Introduction

For bacterial cells steadily growing under fixed, substrate-rich culture conditions, several quantitative laws pertaining to the physiology of growing bacteria have been identified. These include the linear association between growth rate and the ribosomal fraction in total proteomes Schaechter et al. (1958); Scott et al. (2010), the linear relationship linking substrate uptake to growth rate Pirt (1965), and the correlations between transcriptional responses to various sublethal stresses Kaneko et al. (2015).

In recent years, the dynamical aspects of bacterial physiology has garnered significant attention. A wealth of intriguing phenomena under such conditions has been gradually revealed. For instance, cells can “remember” the duration of starvation and the nature of environmental changes—such as the rate at which nutrients deplete—and this can dramatically alter the distribution of lag times Kaplan et al. (2021); Himeoka and Kaneko (2017); Levin-Reisman et al. (2010); Himeoka et al. (2022a). Both the proteome and metabolome of starved cells display characteristically slow dynamics Radzikowski et al. (2016). The metabolites used for storage play a pivotal role in fluctuating environments Sekar et al. (2020). Further, there is a trade-off between the growth rate under nutrient-rich conditions and other physiological parameters: namely, the lag time upon nutrient shifts Basan et al. (2020) and the death rate under starved conditions Biselli et al. (2020).

To achieve a comprehensive understanding of physiological responses, the development of quantitative theories that elucidate the dynamics of cellular physiology is indispensable. Metabolism, often a primary trigger for physiological transitions, stands at the forefront of efforts to forge a unified perspective on the dynamics of bacterial physiology. To this end, extensive research has been undertaken to understand the dynamic responses of metabolic systems using mathematical frameworks; The metabolic control analysis (MCA) Heinrich and Rapoport (1974), the biochemical systems theory (BST) especially S-systems Savageau (1988), and the linear stability analysis of massbalance kinetic models Chakrabarti et al. (2013); Lee et al. (2014). The exploration of the study on the stability and responsiveness of the metabolic system is limited only to the linear region except S-systems. However, the S-system approach omits the mass-balance. The mass-balance condition, or in other words, the stoichiometric constraints, is known to have strong impact on the dynamics of chemical reaction systems Feinberg (2019). For the development of theoretical foundations of metabolic dynamics, it is essential to explore the dynamic responses of mass-balancing metabolism beyond the linear regime.

In this manuscript, we examine the mass-balancing kinetic models of Escherichia coli’s central carbon metabolism and ask how the metabolic state responds to perturbations with which the linear approximation is no longer valid. We study three independently proposed kinetic models Chassagnole et al. (2002); Khodayari et al. (2014); Boecker et al. (2021), and explore the features shared among the models. We demonstrate that the three models exhibit strong responses to perturbations on the metabolic states; The effects of these perturbations are magnified, leading to significant deviations in metabolite concentrations from their steady-state values, depending on which and how metabolite concentrations are perturbed. The observed responsiveness is a hallmark of the realistic metabolic network as the toy metabolic models lack such a strong response. Computational analysis revealed that cofactors, such as ATP and ADP, play a crucial role in the strong response to perturbations. We also discovered the importance of the sparse structure of the metabolic network by adding virtual reactions to it. To validate our findings, we developed a simple, minimal model of metabolic reactions, and we confirmed that cofactor dynamics and network sparsity are indeed the key ingredients of the strong responses.

Results

Model

In the present study, we investigate the dynamic response of metabolic states to perturbations in metabolite concentrations using kinetic models of E.coli central carbon metabolism. In kinetic modeling approaches, the temporal evolution of metabolites’ concentrations is modeled using ordinary differential equations to capture behaviors in out-of-steady-state metabolism.

We employ three kinetic models of E.coli’s central carbon metabolism, emphasizing the common features exhibited by the models. A limitation of kinetic models is their specificity; unlike constraint-based modeling, kinetic modeling necessitates extensive biochemical information, such as the function form of the reaction rate equation and the parameter values for each reaction. By examining a single kinetic model, conclusions may apply only to that specific model. Therefore, our focus lies on the shared features of the three models to circumvent the limitations of individual model specificity. We also believe there is potential for the results obtained commonly among the models to apply to real biological systems as well.

The models under consideration are those proposed by Chassagnole et al. Chassagnole et al. (2002), Khodayari et al. Khodayari et al. (2014), and Boecker et al. Boecker et al. (2021). While all three models incorporate the glycolytic pathway, only the models by Khodayari et. al. and Chassagnole et. al. feature the pentose phosphate (PP) pathway. Notably, the Chassagnole model excludes the tricarboxylic acid (TCA) cycle. A graphical summary of the metabolic modules in these models is presented in Fig. 1A.

Figure 1.

The Boecker and Chassagnole models utilize glucose as the sole carbon source, whereas the Khodayari model can assimilate other carbon sources, such as fructose, formate, and acetate depending on differences between the extracellular- and intracellular concentrations. The Boecker model explicitly models biomass formation, considering several metabolites as biomass precursors that are consumed with a fixed stoichiometry to produce biomass. Conversely, the Khodayari and Chassagnole models do not incorporate the biomass formation reaction. Instead, biomass precursor metabolites are consumed independently, implying that the stoichiometry of biomass formation varies depending on the concentrations. The dilution effect due to the volume growth is modeled in both the Boecker and Chassagnole models. In the Boecker model, the dilution rate corresponds proportionally to the biomass formation reaction rate, while it remains constant in the Chassagnole model. The Khodayari model omits the dilution effect. In addition to the metabolic reaction implemented in each model, the specific forms of the reaction rate equations also differ among the models.

Within these models, the dynamics of gene regulation mediated by transcription factors are not modeled; that is, the total concentration of each metabolic enzyme is constant while the substralevel regulations are incorporated into the models (The lists of substrate-level regulations are provided in Table. S1-S3. For a detailed description, see Materials and Methods).

Perturbation-response simulation

In this section, we present the results of a series of simulations, which we call perturbation-response simulations. These simulations highlight the unique responses of the metabolic models to perturbations and shed light on the dynamic behavior of metabolic systems. It may be useful for the reader to first describe the perturbation-response simulation procedure. The procedure is consistently applied across all three models used in our study and can be summarized as follows:

1. Compute the growing steady-state attractor.

2. Generate N_ini initial points by perturbing the metabolite concentrations from the attractor.

3. Simulate the model dynamics starting from each initial point.

1. Computing the attractor

First, we identify the steady-state attractor where the production/consumption of each metabolite is balanced. For the Khodayari and Boecker models, the attractor corresponds to the steady-state studied in their original papers. For the Chassagnole model, we numerically determined the steadystate attractor (details provided in Materials and Methods).

2. Generating the initial points by perturbation

Once the steady-state attractor has been computed, we establish a set of initial points for subsequent computations by perturbing the metabolite concentrations. This perturbation is a proxy for stochastic fluctuations within cells. We postulate that the perturbation source in the variability in protein concentrations due to the inherent randomness of transcription and translation processes. Also, the cell division would be another source of perturbations of metabolites’ concentrations while the cell division is not explicitly considered in the study. Based on a stochastic model of transcription and translation with a biologically relevant setup, the relative size of the concentration fluctuations is estimated to be several tens of percent (for a detailed derivation, see Materials and Methods). Consequently, the perturbed concentration of chemical n as , where r_n is a uniformly distributed random number with a maximum strength of 40%, i.e. ranging from 0.6 to 1.4. The perturbation strength is beyond the linear region that the relaxation dynamics is approximated by the linearlized model around the steady state. The choice of uniform distribution is not essential for the following result. Also, the characteristic features of the response studied later remain consistent even when different perturbation strengths are employed (see Fig. S2).

3. Simulating the dynamics

We compute the model dynamics using the N_ini points generated by the perturbation as initial points to explore the metabolic response to the perturbations. As far as we have tried, all initial points have returned to the original growing steady state, and thus, we hereafter study the relax-ation behavior.

Strong responses in the kinetic models

After executing the perturbation-response simulations for the models, we observed distinct relaxation dynamics in each. In Fig. 1B-D, the top and bottom panels display typical examples of dynamics from each model that either respond strongly or weakly to perturbations. In the weakly responding dynamics (bottom panel), the effects of the perturbations diminish almost monotonically. Conversely, in the dynamics that respond strongly (top panel), the initial displacements of concentrations are amplified over time, leading some metabolites’ concentrations to overshoot or undershoot (significant concentration changes in metabolites are highlighted in Fig.S3).

For more in-depth analysis, we introduce a measure of the responsiveness to the perturbation that we call the response coefficient χ, defined as

where ‖⋅‖ is the Euclidean norm. The response coefficient distributions for each model are plotted in the bottom panels of Fig. 1E-G. All three models have the peak at χ = 1 (note that this is the minimum value of χ). These peaks correspond to the trajectories where the Euclidean distance between the state and the attractor decreases monotonically with time (i.e. max_t ‖ ln x(t) − ln x^st ‖ = ‖ ln x(0) − ln x^st ‖). In addition, the distribution of the Khodayari model exhibits another peak in a larger χ region while the other two models have plateaus adjacent to the primary peak. We have computed the response coefficient by linearizing each model. As shown in Fig. S1, the responsive-ness of the original model cannot be described by the linearized model.

For comparison purposes, we conducted perturbation-response simulations for the random catalytic reaction networks (RCRN) model, a toy model of the metabolic network Furusawa and Kaneko (2003, 2012). In one instance of the RCRN model, the “metabolites” are interconnected through a random network, and every metabolite is considered a catalyst as well as a reactant. We generated 128 instances of the RCRN model, each having the same number of metabolites and the same distribution of the reaction rate constant. For instance, to draw a comparison with the Khodayari model, we generated the model instances with 778 metabolites and 3112 reactions, where the reaction rate constants are randomly chosen from that of the Khodayari model (detailed construction of the RCRN model is provided in Materials and Methods). The top panel of Fig. 1E-G presents the response coefficient distribution for each RCRN model instance, represented by different lines. Additionally, the “average” response coefficient distribution, which is the mean of the distributions plotted in the top panel, is overlaid (green dotted line). As a fact, the emergence of the additional peak and plateaus, indicative of strong responsiveness, are unique traits of the realistic metabolic reaction network model and are seldom seen in toy representations of the metabolic network.

Key metabolites on the responsiveness

What factors contribute to the strong responses of the metabolic state? In this section, we investigate which metabolites play pivotal roles in inducing these strong responses in the model. To understand the role of each metabolite’s dynamics, we conducted perturbation-response simulations using model equations derived by fixing each metabolite’s concentration to its steady-state value, i.e, set and dx_met/dt = 0 for the concentration of the metabolite met whose con-centration is fixed. The impact of holding each metabolite’s concentration constant is quantified by the change in the average response coefficient:

where x denotes the metabolites’ IDs and χ_x represents the response coefficient of the model with the concentration of the metabolite x fixed to constant. Note that for each metabolite in each model, either the or is defined, depending on whether the modified model (with the metabolite x as a held constant) exhibits an average response coefficient larger (for) or smaller (for) than the original model.

Fig. 2 displays values for representative metabolites (comprehensive results for all metabolites can be found in Fig. S4). While most metabolites have a bilateral effect on the average response coefficient depending on the original model (or exert a small unilateral effect), ATP and ADP (we refer to them collectively as AXPs) consistently demonstrate a marked impact in reducing the average response coefficient across models. This observed significance of AXPs in metabolic dynamics aligns with previous findings Himeoka and Mitarai (2022). The common importance of AXPs are shown also when we quantify the impact of their dynamics on the responsiveness by utilizing the Sobol’ total sensitivity index (Fig. S5) Homma and Saltelli (1996).

Figure 2.

Interestingly, halting the dynamics of several metabolites can amplify the model’s response coefficient. Across the models, metabolites that augment the response coefficient when their concentrations are fixed tend to be allosteric regulators associated with negative regulation. Freezing these metabolites’ concentrations effectively removes the negative feedback loop present in the metabolic systems. Typically, negative feedback stabilizes system behavior. In the present context, it assists the system in returning to its original state after perturbations. Hence, fixing these metabolites’ concentrations allows the system to respond more vigorously to perturbations compared to the original model. We checked the effect of feedback itself on the responsiveness rather than the concentration dynamics by weakening the strength of substrate-level regulation of OAA, FUM, AKG in the Boecker model (Fig. S6). It was shown that weakening the substrate-level regulation of CSICD by AKG and MDH by OAA increases the responsiveness. On the other hand, weakening the strength of substrate-level regulation of neither FRD by OAA nor MDH by FUM alters the responsiveness, and in addition there is no substrate-level regulation by SUC. The increase in responsiveness by fixing the concentrations of OAA and AKG would be explained by substrate-level regulation of the corresponding metabolites. However, substrate-level regulation is not sufficient to describe the effect of fixing the concentration of FUM and SUC. Higher order effects may be required.

Pyruvate is an exception to the above description, which is not the regulator in the models. The strong influence of pyruvate on the increase in response coefficients is attributed to the phosphotransferase system (PTS). The external glucose is taken up by converting phosphoenolpyruvate into pyruvate through the PTS. This reaction elevates pyruvate concentration, subsequently decelerating glucose uptake. By holding pyruvate concentration constant, this negative feedback effect is negated, enhancing the autocatalytic nature of the reaction system formed by the PTS and glycolysis (GLC_ex + PEP → … → 2PEP).

The role of the sparsity of the networks

The key role of AXPs in changing the kinetic models’ responsiveness was shown in the previous section. Nevertheless, this does not imply that all kinetic models of chemical reaction systems with cofactors such as AXPs always exhibit a strong response to perturbations. The strong responses demonstrated by the three models (as seen in Fig. 1B-D) are absent in abstract toy models of metabolism, which consist of random reaction networks incorporating such cofactor chemicals Kondo and Kaneko (2011). This discrepancy prompts us to investigate other potential determinants causing metabolic models to display strong responses shown in the three models.

A noteworthy attribute of metabolic networks is their inherent sparsity. While cofactors are linked to numerous reactions, the backbone networks—networks excluding cofactors—are distinctly sparse. As an example, the glycolytic pathway predominantly mirrors a linear sequence of reactions, resulting in many metabolites within this pathway participating in merely two reactions.

Does network structure play a role in the dynamics of metabolic systems? For this question, pre-vious research in the opposite direction—reaction dynamics on dense networks—provides useful insights. There are several studies using the random networks as an abstract model of metabolic networks Furusawa and Kaneko (2003); Kaneko et al. (2015); Awazu and Kaneko (2009); Himeoka et al. (2022b). In such types of models, the effect of the initial perturbations on metabolite concentrations typically decays monotonically over time Awazu and Kaneko (2009); Himeoka et al. (2022b). This is in stark contrast to the behavior observed in Fig. 1B-D.

Considering the inherent tendency of random reaction networks to display dynamics characterized by weak responses, we are compelled to investigate the influence of network sparseness on the dynamic behavior of metabolic models. To this end, we increase the reaction density of the metabolic network in the three models and assess the resultant changes in responsiveness. Our approach entails initially introducing N_add random reactions into the network and subsequently al-locating kinetic parameters to these additions from the distribution of parameter values computed from the original model. We then execute a perturbation-response simulation on this extended model to determine the response coefficient, which in turn helps ascertain the distribution and average of these coefficients. This process is schematically depicted in Fig. 3A.

Figure 3.

Fig. 3B illustrates the alterations in responsiveness resulting from the model extensions. We generated extended metabolic networks and accompanying ODE models by adding N_add random reactions, and executed the perturbation-response simulation. The average response coefficient is then computed for each model and the distribution of the average response coefficient is presented as the violin plots. In the perturbation-response simulation, the extended models with only a single attractor are subjected to the response coefficient analysis (comprehensive procedural details available in the Materials and Methods section). Note that in this procedure, metabolites are not newly introduced to the model, but only the reactions are. Thus, the reaction density increases.

As evident from the figures, a consistent trend emerges across the three models: the introduction of random reactions weakens responsiveness. However, the magnitude of this reduction varies among the models. It is important to note that the reactants of the added reactions are selected randomly, thus not based on any specific chemical rationale. Nevertheless, this diminished responsiveness persists even when the additional reactions are confined solely to those cataloged in the metabolic model database, as visualized in Fig. S7.

Such attenuation aligns with expectations shaped by prior studies on dense reaction networks Furusawa and Kaneko (2003); Kaneko et al. (2015); Awazu and Kaneko (2009); Himeoka et al. (2022b). Evidently, the inherent sparsity of the backbone reaction network plays a critical role, enabling the kinetic model of metabolic systems to respond strongly to perturbations.

A minimal model for disrupted homeostasis

So far, we have studied the kinetic models of the central carbon metabolism of E.coli. Our findings indicate that cofactors and the sparsity of the backbone network (network devoid of cofactors) are the keys to the responsiveness of the metabolic dynamics. In this section, we develop a simple minimal model to investigate if these two factors are sufficient to induce strong responses in reaction systems.

Our minimal model is designed in a manner that allows us to adjust the sparsity of the backbone network and the proportion of cofactor-coupled reactions. First, we explain the construction of the backbone networks. These networks are comprised of N chemical species, C₁, C₂, …, C_N, and R re-actions (R ≥ N −1). For simplicity, we construct the backbone network using only uni-uni reactions. The minimum number of uni-uni reactions required to interlink N chemical species is N − 1. To assemble a backbone network with N chemicals and R reactions, we first connect N reactions using N − 1 reactions and subsequently introduce random reactions for the pairs of chemicals that lack reactions, totaling R−(N − 1). We incorporate exchange reactions with extracellular environments for 5% of the metabolites (this percentage is not critical to the results). The exchange reactions are always introduced to the two endpoints metabolites of the sequential linkage to ensure there are no dead-ends on reaction networks. We posit that one of the endpoint metabolites from the sequential linkage is the nutrient metabolite, and we set its external concentration at 100, while for the remainder, it is set to unity.

As cofactors, we introduce the three forms of the activation for mimicking AXP, A^**, A^* and A, which are not included in the N chemicals in the backbone network. We introduce the coupling fraction f to modulate the number of reactions with cofactor coupling. With the coupling fraction f, the number of reactions with cofactor coupling is given as Round(fR) with Round(⋅) as the function to round the value to the nearest integer. When a reaction is selected for cofactor coupling, two forms of the cofactor (chosen from A^**, A^* and A) are randomly selected. The reaction then gets coupled to the conversion between these two cofactor forms, with the direction being randomly set.

Overall, the differential equation governing the minimal model can be expressed as

In this expression, the concentrations of the chemicals from the backbone network, as well as the cofactors, are combined in a vertical stack, i.e., n = 1, 2, …, N, N + 1, N + 2, N + 3. x_n with n = 1, …, N represents the concentrations of the nth chemical in the backbone network, and x_N+1, x_N+2 and x_N+3 denote the concentrations of A^**, A^* and A, respectively. S is the stoichiometric matrix of the reaction network. D_n is unity if the nth chemical has an influx/efflux term, otherwise, it is zero. X_n is the external concentration of chemical n, which is a constant parameter.

Mass action kinetics is employed for the reaction rate function J_r. The rate function has the form

where Sub(r) and Prd(r) are the set of substrates and products, respectively, of the rth reaction. If the reaction is coupled with the cofactor, |Sub(r) | = |Prd(r) | = 2 and otherwise |Sub(r) | = |Prd(r) | = 1. We neither allowed the self-loops in the backbone network (C_n ⇌ C_n) nor the catalytic activity of the cofactors (e.g., C_n + A^* ⇌ C_m + A^*).

The parameters v_r and are the reaction rate and irreversibility, respectively. The v_r values are randomly set as with the uniformly distributed random number, u_r ∼ U(−3.66, 7.13). The range for the kinetic parameter is derived from the Khodayari model. is calculated using the Arrhenius equation based on the standard chemical potential of each chemical. We randomly set the standard chemical potential μ_n for each chemical, where μ_n follows the uniform distribution U(0, 1), while the standard chemical potentials of A^**, A^* and A are fixed to the constant values, 1.0, 0.5, and 0.0, respectively. To satisfy the detailed balance condition, are then given by with Δμ as the standard chemical potential difference between the products and the substrates, Δμ_r = ∑_n∈Prd(r) μ_n − ∑_n∈Sub(r) μ_n, and β is the inverse temperature.

We generate M networks for several choices of reaction number R and coupling fraction f. We performed the perturbation-response simulation to determine the response coefficient distribution of each network. In Fig. 4, the average response coefficient, averaged over the networks, is depicted as functions of R and f (see the rectangles labeled as “Cofactor”). These figures clearly illustrate that the average responsiveness tends to decrease as either R increases or f decreases. The coupling with the cofactors leads to the two-body reactions, and thus, increases the nonlinearity of the model equation. It is noteworthy that the increase in responsiveness with f cannot be fully attributed to the increase in the model’s nonlinearity. To validate this, we employed the model (3) without the cofactors, but incorporated catalytic reactions. In this model, each reaction proceeds with the catalytic ability of an enzyme, represented as C_n +C_l ⇌ C_m +C_l (here, autocatalytic reactions are permissible, meaning l can be either n or m).

Figure 4.

In terms of the numbers of constant, linear, and quadratic terms in the model equation, the catalytic reaction version of the model corresponds to the f = 1 case of the original minimal model. If the increase in responsiveness with f is dictated solely by nonlinearity, the catalytic reaction model should have comparable responsiveness with the “cofactor” model with high f. Yet, as shown in the rectangles above the main plot in Fig. 4 (labeled as “Catalytic”), the responsiveness metrics are notably lower than those observed in the minimal model.

Discussion

In the present manuscript, we have studied three models of E. coli’s central carbon metabolism. Our findings indicate that these models commonly exhibit strong and weak responses to perturbations of the steady-state concentration. In the weak response, the effect of the perturbation dwindles almost monotonically over time. In contrast, during a strong response, initial deviations from the steady state get amplified. Note that the linearized models cannot capture the response dynamics. We introduced the response coefficient, χ, as a metric to quantify the strength of these responses and undertook comprehensive computational analyses. The strong responsiveness were observed not only in the three models studied in the manuscript, but also in a model of glycolytic pathway with multiple sets of parameter choices (Fig. S8) Haiman et al. (2021).

We aim to identify key metabolites whose dynamics play a crucial role in the emergence of strong response. We systematically fixed the concentration of each metabolite to its steady-state value and assessed the influence of individual metabolite dynamics on the overall response strength. For a majority of metabolites, the impact on the response coefficient is bidirectional: halting the dynamic concentration of a metabolite can either amplify or diminish the response, depending upon the base model. Even if their effects consistently amplify or diminish the response, the magnitude of change in the average response coefficients remains relatively modest.

However, three metabolites—ATP, ADP, and pyruvate—stood out. They consistently and significantly influenced the average response coefficient across models. Pyruvate’s pronounced impact can be attributed to the phosphotransferase system, whereas the effects of ATP and ADP seem to arise synergistically across multiple reactions.

One of the authors has recently reported a kinetic model of E.coli metabolism to exhibit anomalous relaxations and the role of AXPs in the relaxations Himeoka and Mitarai (2022). In the paper, they reduced the metabolic network into a coarse-grained minimal network. The minimal network revealed that, depending on the initial conditions, a futile cycle conversing ATP into AMP be overly active. Then, the futile cycle can deplete ATP and ADP, leading to a significant slowdown or even cessation of several ATP/ADP-coupled reactions, pushing the entire metabolic network into a “choked” state. The extensive connectivity of AXPs is a common trait in metabolic models. Hence, it is plausible that the impact of AXPs’ dynamics on the entire metabolic system is a common feature irrespective of specific kinetic models.

Several experimental studies underscore the pivotal role of AXPs in metabolic homeostasis. For example, Teusink et al. demonstrated that in yeast, a sudden increase in glucose levels in the culture medium can activate a positive feedback loop in ATP production in the glycolytic pathway, leading to growth arrest van Heerden et al. (2014); Teusink et al. (1998). This phenomenon of growth arrest, triggered by a nutrient level surge, has been observed across various bacterial species, including E.coli Strange and Hunter (1966); Postgate and Hunter (1964); Calcott and Postgate (1972), though its molecular mechanism has been most extensively studied in S. cerevisiae. Moreover, a rapid decrease in the total adenyl cofactor concentration as a cellular response to external glucose up-shift has been documented Theobald et al. (1997); Kresnowati et al. (2006); Rizzi et al. (1997); Walther et al. (2010). Walther et al. found that inhibiting this transient reduction in the adenyl cofactor pool increased the number of cells transitioning to a growth-arrested state upon glucose up-shift Walther et al. (2010).

Next, we focused on the structure of the backbone network as suggested by the relatively simple dynamics typically observed in the random network metabolic models Furusawa and Kaneko (2003); Kaneko et al. (2015); Awazu and Kaneko (2009); Himeoka et al. (2022b). To this end, we increased the density of the backbone network in the three models by randomly selecting metabolites and introducing hypothetical reactions between them. The perturbation-response simulation on the extended models indicated that as more reactions were integrated, the system’s response weakened.

Finally, we devised a minimal model to evaluate our hypothesis that the dynamics of cofactors and the sparsity of the backbone network are pivotal determinants of responsiveness. This model successfully showed a strong response to perturbations when the fraction of reactions to be coupled with cofactors becomes high and the backbone network is sparse. It is worth noting that the heightened response due to cofactor coupling is not solely a result of the model equation’s increased nonlinearity. Interestingly, the model’s catalytic reaction variant demonstrated attenuated responses to perturbations, even though its nonlinearity mirrored that of the minimal model where all reactions were cofactor-coupled.

The present study showed that there are similar metabolic responses to perturbations across three distinct kinetic models of E.coli metabolism. We like to underscore that these models were proposed by different research groups and within distinct contexts. Furthermore, the specifics of each model diverge considerably: both the Khodayari and Chassagnole models are tailored for aerobic conditions, whereas the Boecker model is for anaerobic conditions. Certain reactions and metabolites present in the Khodayari model are absent in both the Boecker and Chassagnole models. The reaction rate functions and parameter values differ across these models because each model has been developed to recapitulate different experimental results. Yet, in spite of these pronounced disparities, the models exhibit converging characteristics to a certain degree, notably the pivotal role of AXPs and the diminishing responsiveness attributed to the addition of reactions. Such congruence in outcomes might hint at the presence of features that remain resilient against variations in model specifics, potentially being shaped by the foundational network structures inherent to core metabolic processes.

From a theoretical standpoint, the robustness of metabolic systems, as well as the interplay between network structure and reaction dynamics, have been explored in the chemical reaction network theory and related fields Feinberg (2019); Gutenkunst et al. (2007); Himeoka et al. (2022b); Kobayashi et al. (2022). However, the practical relevance of these theories to metabolic systems remains to be ascertained since the theories are mathematically rigorous and the strong assumptions narrow down the applicable systems. The cross-talk between bacterial physiology and chemical reaction network theory should enhance the developments of each field.

In the following paragraphs, we discuss potential implications of strong responsiveness in the context of biological sciences. The unique feature of strong responsiveness, as observed in our realistic metabolic model, naturally prompts the question: what could be its biological function? An initial hypothesis suggests its role in sensing ambient condition changes and fluctuations. Wellknown design principles of signal transduction systems emphasize the amplification of external stimuli Ferrell and Ha (2014). Given that certain transcription factors are known to sense concentrations of metabolites in central carbon metabolism Kotte et al. (2010); Okano et al. (2019); Kochanowski et al. (2013), strong responsiveness could serve as an amplifier for external cues.

Exploring metabolic dynamics also holds promise for synthetic biology and metabolic engineering Shimizu et al. (2001); Schwander et al. (2016); Bar-Even et al. (2010); Doerr et al. (2019); Claassens et al. (2019); Kurisu et al. (2023); Opgenorth et al. (2017, 2014). The construction of artificial metabolic reaction networks often relies on pathway prediction algorithms, constraint-based models, and machine learning Moriya et al. (2010); Araki et al. (2015); Orth et al. (2010); Varma and Palsson (1993); Segler et al. (2018). While these algorithms consider stoichiometric feasibility, they frequently overlook dynamics. This oversight can result in proposed network structures underperforming Schwander et al. (2016); Bar-Even et al. (2010); Doerr et al. (2019); Claassens et al. (2019), especially in terms of cofactor recycling Opgenorth et al. (2017, 2014). Notably, in silico reconstructed metabolic systems often disrupt cofactor balance Okamoto et al. (1980); Okamoto and Hayashi (1983), posing challenges to bottom-up synthesis in artificial cells. We hope our study offers foundational insights for the development of quantitative metabolic dynamics theories.

Acknowledgments

The authors thank Namiko Mitarai and Taro Furubayashi for fruitful discussions. This work is supported by JSPS KAKENHI Grant Numbers JP 21K20626, JP 22K15069, and JP 22H05403 to YH; JP 22K21344 to CF;, the Japan Science and Technology Agency (JST) ERATO (JPMJER1902 to CF), and GteX Program (JPMJGX23B4 to YH).

Code availability

Codes for generating the computational models studied in the manuscript are deposited on GitHub (https://github.com/yhimeoka/Perturbation-Response-Analysis.git).

Materials and methods

Models and modifications

In Khodayari et al. Khodayari et al. (2014), the kinetic parameter values are estimated by using the ensemble modeling Khodayari et al. (2014); Tran et al. (2008); Tan and Liao (2012); Rizk and Liao (2009) so that the steady flux distribution of the metabolic reactions becomes consistent with the experimental measurements. However, the steady-state dealt with in the paper turned out to be unstable, i.e., the maximum eigenvalue of the Jacobi matrix of the linearized system at the steady state was positive. For carrying out the perturbation-response simulation, we need to make the steady-state stable. To make the steady-state stable, we fixed the periplasmic hydrogen concentration to the steady-state value and removed the competitive inhibition of phosphoenolpyruvate carboxylase (Those two factors were manually found). Also, we set the concentrations of the metabolites in the culture media (extracellular metabolites) with the assumption that the model cell is growing in a sufficiently large flask. In addition, the small molecules, H₂O, O₂, CO₂, NH₄ and phosphate, are set to be constant with further assuming that the exchange rate of the small chemicals between the extracellular culture is sufficiently fast.

For the dynamics of ATP and ADP concentrations, in the Chassagnole model Chassagnole et al. (2002), we used the reactions that coupled with ATP and ADP. Since the kinetic rate equations of those reactions are defined with the dependency of ATP and ADP concentrations in the original model, we utilized the kinetic rate equations as they are. The reactions for the adenine base synthesis are not modeled in the model, and thus, the total concentration of ATP and ADP are the conserved quantity in the model. By assuming that the decreases of the ATP and ADP due to the growth dilution are compensated, the growth dilution term is not introduced to the model equa-tions of ATP and ADP concentrations. Lastly, we introduced the non-growth associated ATP hydrolysis reaction (ATP→ADP) since the model fails to balance ATP and ADP and relaxes to a steady state where the concentration of some metabolites is almost zero. We utilized a simple MichaelisMenten rate equation for the hydrolysis reaction J_AH = vx_ATP/(K_M + x_ATP), where v = 0.1 (sec⁻¹) and K_M = 1.0 (mM)

The Boecker model Boecker et al. (2021) is used without modification.

Computation of dynamics

All the ODE simulations in the present manuscript are done using ode15s() function in MATLAB (MathWorks) Inc. (2022). For computing the attractor of the Chassagnole model, we used the initial concentrations registered in chassagnole4 on JWS online Olivier and Snoep (2004) as the initial concentrations of the modified Chassagnole model to run the dynamics. The steady-state concentrations of the metabolites reach from the initial concentrations and are utilized as the attractor.

Estimate of the fluctuation strength

For estimating the strength of fluctuations, we utilize the following simple metabolic reaction

The concentration of the enzyme P and Q follows the stochastic transcription and translation model Friedman et al. (2006); Paulsson and Ehrenberg (2000). The temporal change of the concentration of the metabolite X here is given by a simple model with the mass-action kinetics,

where x, p, q are the concentration of X, P, and Q, respectively. v and u are the kinetic parameters of the reactions. According to the time-scale separation of metabolic reactions and transcription/translation, we suppose that x is the fast variable and its dynamics slave to p and q. We assume that the timescales in the change of the metabolite concentrations are sufficiently smaller than that of the proteins and apply the quasi-steady-state approximation: we solve the steady-state of the deterministic equation of (5) Kaneko (1981); Risken (1996). Then, we calculate the mean and variance of the solution x(t) = vp(t)/uq(t), which is a stochastic variable. Since the stochastic transcription and translation model leads to the Gamma distribution of the protein concentration, the average and the variance are given by

where a and b are the parameters set by the rate of the transcription, translation, mRNA degradation, and protein degradation.

The coefficient of variation is then given by

a is the ratio between the transcription rate and the protein degradation rate. According to the large-scale proteomic analysis of E.coli Taniguchi et al. (2010), the mean value of a among the essential genes is ≈ 6.82 ± 2.34 (only with the essential metabolic enzymes listed in the iML1515 Monk et al. (2017), the average value is ≈ 6.72 ± 2.34). By substituting a = 6.82 into (8), we obtain the relative noise strength as ≈ 62%. Since in the estimate, we suppose the irreversible reactions lead to larger noise amplitudes than reversible reactions, we regard 62% as the maximum relative noise strength and use a bit smaller value, 40%.

The random catalytic reaction network model

The random catalytic reaction network (RCRN) model is a simple, mass-balancing toy model of cellular metabolism Furusawa and Kaneko (2003); Kaneko et al. (2015). The ordinary differential equation dictating the temporal evolution of the ith metabolite’s concentration x_i is given by

where v_ijk, D_i and μ are the rate constant of the reaction C_j → C_i catalyzed by C_k (C_i denotes the ith metabolite), the substrate update rate of C_i, and the specific growth rate of the model cell, respectively. T (i) denotes the transporter metabolite (enzyme) for the uptake of C_i, and is the external concentration of C_i. In the RCRN model, the total volume of the cell is usually set to be equal to the total amount of the metabolites. This is because each “metabolite” in the model is interpreted as an enzyme as well as a reactant of the metabolic reactions. Hence, the specific growth rate of the model cell is set to the total uptake rate of the metabolites,.

In the construction of an instance of the RCRN model, we first generate a random, connected network among N metabolites, and then, we assign a single catalyst on each reaction edge. If the reaction C_j → C_i catalyzed by C_k exists in the network, v_ijk is non-zero, while otherwise, it is set to zero. In the present manuscript, we suppose all the reactions are reversible, and thus, v_jik is non-zero if and only if v_ijk is non-zero.

For making the comparison of the responsiveness of the realistic model to the RCRN model, we generated 128 instances of the RCRN model where the number of metabolites N is set to be equal to the realistic model that we want to make a comparison. The number of reactions R is set to 3.5N for the Boecker and Chassagnole model while 4N for the Khodayari model. The ratio between N and R is set so that the Erdős–Rényi random graph of given size becomes connected with sufficiently high probability. The reaction rate constants v_ijk’s are sampled from the model to compare. For the sake of simplicity, we suppose only the 1st metabolite is the nutrient metabolite which is taken up from the external environment and the Nth metabolite is the transporter of it.

D_i and are set to unity for i = 1, and for the else, those parameters are set to zero.

For each generated instance of the RCRN model, we performed the perturbation response simulation where the attractor(s) are computed by simulating the model dynamics from randomly generated 32 initial concentrations x_ini ∈ (10⁻³, 10³)^N. The RCRN models exhibiting multistability or computational failure due to numerical issues such as heavy stiffness are not subjected to further analysis (Less than 10% of the model instances showed multistability for corresponding models. The computational failure is ≈ 20% for the RCRN model with the parameter distribution of Chassagnole model, ≈ 50% for that of the Khodayari model, while less than 1% for that of the Boecker model). For quantifying the response coefficient, we computed 128 trajectories for each of those models.

Random addition of the metabolic reactions

First, we have add uni-uni reactions where the reactants are randomly picked from the non-cofactor metabolites. Next, we decide whether the additional reaction to be coupled with the adenyl cofactor metabolites (ATP, ADP, or AMP) with a propability p. The probability p equals to the fraction of adenyl cofactor-coupled reactions in the original model. The addition reaction is either the uni-uni reaction (if the reaction is not coupled with AXPs) or the bi-bi reaction (if the reaction is coupled with the adenyl cofactors). We computed the parameter value distribution for each reaction scheme and the parameter values for the additional reactions are sampled from the corresponding parameter value distribution.

For the Khodayari model, the additional reaction is modeled using the elementary reaction decomposition (see the next section), i.e, a reaction A ⇌ B is decomposed into the three elementary reaction steps, E + A ⇌ EA, EA ⇌ EB, and EB ⇌ E + B. The mass-action kinetics is adopted as the reaction rate function for each elementary reaction step.

The elementary reaction decomposition is not applied to the Boecker model and the Chassagnole model. Also, several types of reaction rate functions are used in the two models. Thus, we decided to use the simplest enzymatic reaction kinetics used in each paper. The reversible Michaelis-Menten kinetics

is used for uni-uni reactions, and for bi-bi reactions, the ordered bi-bi reaction kinetics

is adopted in the Boecker model, a reaction kinetics

is used in the Chassagnole model. [S_i], [P_i], and K_* (* is either S₁, S₂, P₁ or P₂) represent the concentration of the ith substrate and product, and the dissociation constant of the corresponding chemical, respectively.

After constructing an extended model, we compute the steady-state attractor by simulating the dynamics from 128 random initial concentrations. These initial conditions are generated by applying a 40% relative perturbation to the steady-state concentrations of the original model. If all the initial conditions converge to a single attractor, then we generate another set of 128 initial conditions by applying a 40% relative perturbation to the steady-state concentrations of the extended model in order to perform the perturbation-response simulation. For the Boecker model, we require the extended model to have a steady growth rate greater than half of that of the original model. This is because growth dilution is typically the slowest process in the model, and thus, a significant slowdown of growth dilution can lead to an artificial overestimation of the model’s responsiveness.

Elementary reaction decomposition

Here, a brief description of the elementary reaction decomposition (ERD) is provided. Suppose that we have the following bi-to-uni reaction

and that the reaction is catalyzed by an enzyme E. For implementing this reaction into the kinetic model, the Michaelis-Menten type kinetics of chemical reaction is often adopted. With the Michaelis-Menten type kinetics, the rate of the reaction J is given by

where we supposed that the binding of the substrate to the enzyme occurs in the order of A → B.

The Michaelis-Meneten type reaction kinetics ((10)) is obtained by the pseudo-equilibrium ap-proximation or the steady-state approximation of the following reactions

In the ERD scheme, each step ((11)-(14)) of the overall reaction ((9)) is modeled using the massaction kinetics. Also, the forward- and the backward reactions are dealt as different reactions. For instance, the rates of the reaction A + E ⇌ EA is given by J_A,E→EA = v_A,E→EA[E][A] and J_EA→A,E = v_EA→A,E[EA].

The advantage of the ERD is that every reaction rate function is given either by a linear or quadratic function of the metabolites’ concentrations. This feature allows us, in the present manuscript, to randomize the reaction network in a unified manner. Also, the ERD has the advantage of reduc-ing the computational cost for the parameter estimation based on omics data. For more details, see original papers Khodayari et al. (2014); Tran et al. (2008); Tan and Liao (2012); Rizk and Liao (2009).

The generative AIs and other software

ChatGPT (OpenAI, September 25, 2023 version) and Grammarly (Grammarly, Inc.) are used to improve clarity and to brush up on English grammar.

Supplementary tables and figures

Table S1.

Table S2.

Table S3.

Table S4.

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Figure S8.