To uncover the processes and mechanisms of cellular physiology, it first necessary to gain an understanding of the underlying metabolic dynamics. Recent studies using a constraint-based approach succeeded in predicting the steady states of cellular metabolic systems by utilizing conserved quantities in the metabolic networks such as carriers such as ATP/ADP as an energy carrier or NADH/NAD+ as a hydrogen carrier. Although such conservation quantities restrict not only the steady state but also the dynamics themselves, the latter aspect has not yet been completely understood. Here, to study the dynamics of metabolic systems, we propose adopting a carrier cycling cascade (CCC), which includes the dynamics of both substrates and carriers, a commonly observed motif in metabolic systems such as the glycolytic and fermentation pathways. We demonstrate that the conservation laws lead to the jamming of the flux and feedback. The CCC can show slow relaxation, with a longer timescale than that of elementary reactions, and is accompanied by both robustness against small environmental fluctuations and responsiveness against large environmental changes. Moreover, the CCC demonstrates robustness against internal fluctuations due to the feedback based on the moiety conservation. We identified the key parameters underlying the robustness of this model against external and internal fluctuations and estimated it in several metabolic systems.
In recent decades, mathematical modeling of metabolic systems has been intensively explored in the field of systems biology [1–6]. Several studies have demonstrated that metabolic systems can be quantitatively predicted using mathematical models [7–10]. A particularly important feature for the mathematical modeling of a metabolic system is that some conserved quantities characterize the system. The total number of atoms remains unchanged before and after the reactions, allowing for mass conservation. Furthermore, several coenzymes act as the carriers of molecules and energy through various reactions, e. g. , ATP/ADP as an energy carrier and NADH/NAD+ as a hydrogen carrier [11,12]. The total concentrations of the carriers represent the conserved quantities in the steady state, which is referred to as moiety conservation [11], by which various coenzyme-related reactions have to be balanced. These constraints, i. e. , the law of mass conservation and moiety conservation, restrict the solution space and facilitate the analysis of complicated metabolic networks. One of the most successful approaches used for the modeling of metabolic systems to date is the constraint-based analysis of metabolic fluxes based on the stoichiometry that describes conservation laws in metabolic fluxes [2,5]. In this approach, the metabolic fluxes are assumed to be in a steady state, meaning that, for any metabolite pool, the fluxes governing its synthesis and degradation are balanced. Owing to its predictive power and advantage of not requiring detailed information on the kinetic parameters, a steady state-based metabolic modeling approach has been applied in various studies analyzing the characteristics of metabolic systems. However, in contrast to the steady-state solution, the characteristics of the dynamic behavior of metabolic systems have not been elucidated. In a fluctuating environment, a metabolic system is not in a steady state, and thus its dynamic behavior should be analyzed over time. In a biological system, the producing and decomposition reactions of such carriers are also faster than cellular growth. Since the conservation laws are considered to be maintained for a longer time than several cell cycles, the law of mass conservation and moiety conservation restrict both the steady state and the transient dynamics. Thus, the constraints of metabolic networks should be taken into account when studying metabolic dynamics. In particular, the moiety conserved coenzymes should have particularly strong effects on the various pathways in analyzed metabolic systems due to their recycling [11]. For example, in the glycolytic pathway, ATP is transformed into ADP by phosphofructokinase, while ADP is recycled to ATP in downstream reactions. During the anaerobic growth of microorganisms, NAD+ is transformed into NADH by glyceraldehyde-3-phosphate dehydrogenase in upper glycolysis, while NADH is recycled to NAD+ through the fermentation pathways [1]. These reactions involving carriers must be balanced not only in the steady state but also in the dynamic states, and yet the effect of such balance remains unknown. Several previous studies have analyzed the characteristics of metabolic pathways with moiety conservation. Reich and Sel’kov [11] derived simple dynamical system models of metabolic system, and pointed out that recycling of the moiety conservations represents the skeleton of energy metabolism. Such a moiety-conserved cycle was formulated using mass-action kinetics and was analyzed with dynamical systems theory. They found that the cycle of moiety conservation works as a positive feedback to produce the high-energy carrier autocatalytically if its concentration is low, whereas a higher concentration of the high-energy carrier limits its own production. Following these pioneering works, the behavior of the moiety-conserved metabolic cycle has been extensively studied [13,14]. Although the majority of these studies focused only on the analysis of steady states, some researchers addressed the dynamic behavior of the moiety-conserved cycle after environmental changes [15,16]. These studies demonstrated transient switching behavior based on a simple model with mass-action kinetics without considering complex formation among enzymes, substrates, and cofactors. Recently, some theoretical studies demonstrated that complex formation in a chemical reaction system can cause transient abnormal behavior, which could not be obtained using simple mass-action kinetics [17]. This finding suggested that the complex dynamic behavior of metabolic systems, especially those including moiety conservation, can also be captured by considering complex formation. Following this idea, in the present study, we focused on metabolic reaction pathways with carrier recycling (also known as a “turbo design” [18]). To investigate the dynamic behavior of nonlinear metabolic systems with the law of mass conservation and moiety conservation, we here explore a simple motif designated as a carrier cycling cascade (CCC), which simultaneously considers the dynamics of substrates and coenzymes. We formulate the CCC with complex formation to analyze the dynamics of both substrates and carriers simultaneously. The CCC shows slow relaxation into the steady state against external perturbations around a critical point where the relaxation speed follows the power law. We analyze the origin of the slow relaxation using dynamical systems theory. Additionally, the CCC maintains robustness of a metabolic state with intrinsic noise. We demonstrate that the negative feedback through the conserved carrier provides such robustness, i. e. , the moiety-conserved cycle can reduce fluctuations in the concentrations of metabolites. We further discuss the relationship between microscopic feedback and macroscopic dynamics, and estimate the effect of this feedback using experimentally determined parameters.
The steady-state characteristics in a moiety-conserved cycle were intensively studied by Reich and Sel’kov [11]. However, most of these proposed models are too complicated to effectively extract the essence of the dynamic features. The metabolite flow is usually branched, and the same coenzyme can be utilized at multiple steps. Hence, to investigate the effects of a cycling carrier on the dynamics of metabolic systems, we here focus on the simplest cascade, the CCC, which models ATP-ADP cycling in the glycolysis pathway and NADH-NAD+ cycling in the fermentation pathway (Fig 1). Here, we refer to the high-energy and low-energy carriers as the active and inactive carriers, respectively. The CCC consists only of active carrier-consuming and -producing steps. Our model is similar to the simplest skeleton proposed by Reich and Sel’kov [11]. In contrast to the previous study, which considered the mass-action dynamics without complex formation, we here include the complex-formation process. We assume that the enzyme, substrate, and coenzyme form the complex, while the enzyme is saturated under the ordinal metabolic condition [19,20] and does not appear in the equation explicitly. Original ordinary differential equations (ODEs) are given as seven mass-action kinetic equations (see section 1 in S1 Text), which can be reduced to five ODEs by eliminating the association and dissociation reactions between coenzymes and substrates adiabatically. d [ m 0 ] d t = k in - k c [ m 0 ] [ c ] K 0 + [ c ] - k leak [ m 0 ], (1a) d [ m 1 ] d t = k c [ m 0 ] [ c ] K 0 + [ c ] - k p [ m 1 ] [ c * ] K 1 + [ c * ], (1b) d [ m 2 ] d t = k p [ m 1 ] [ c * ] K 1 + [ c * ] - k out [ m 2 ], (1c) d [ c ] t d t = - k c [ m 0 ] [ c ] K 0 + [ c ] + k p [ m 1 ] [ c * ] K 1 + [ c * ], (1d) d [ c * ] t d t = k c [ m 0 ] [ c ] K 0 + [ c ] - k p [ m 1 ] [ c * ] K 1 + [ c * ], (1e) [ c ] t = [ c ] + [ m 0 ] [ c ] K 0 + [ c ], (1f) [ c * ] t = [ c * ] + [ m 1 ] [ c * ] K 1 + [ c * ], (1g) where mi is the i-th metabolite, c and c* are active and inactive carriers, respectively, and [x] denotes the concentration of x. m0 is supplied and diluted with rates kin and kleak, respectively, and m2 is diluted with rate kout. c and m0 make complex cm0, and c* and m1 make complex c* m1. The active carrier is consumed with rate kc when m1 is transformed from m0, and is produced with rate kp when m2 is transformed from m1. [c]t and [c*]t represent the total concentration of active and inactive carriers as [c]t = [c] + [cm0] and [c*]t = [c*] + [c*m1], respectively. K0 and K1 represent the dissociation constants between c and m0 and between c* and m1, respectively. Here, we considered the condition in which the total concentrations of the carriers are conserved. There are multiple reactions between the active carrier-producing and -consuming reactions in the actual metabolic networks, while these reactions are reversible and not rate-limiting under ordinal metabolic conditions [16,18,20]. Then, multi-step reactions can be reduced to the present form. The conserved quantities in the CCC model are as follows: Note that the number of independent conserved quantities is two because of c pool = c sum + c diff *. In metabolic networks, some types of inputs will be in a CCC, e. g. , changes in the influx and efflux rates of the first metabolite, changes in the ratio between active and inactive carriers due to the environmental changes, and the new synthesis and degradation of a carrier. Here, we consider the condition in which the total carrier concentration is conserved so that the cpool does not change. Therefore, we consider changes in kin and the ratio between active and inactive carriers as inputs. Initially, we changed the influx rate. We allowed the system to approach a steady state, and then changed kin at time = 0 (Fig 2A). When kin is altered, the concentrations of all molecules do not change in the timescales of the active carrier consuming and producing reactions (kc and kp are set as 1. 0). Accordingly, the changing of [m0] is slower than the timescale of kc and kp, while the concentrations of other components do not change as in the quasi-steady state. When [m0] decreases to become sufficiently small, the concentrations of the others drastically change with a fast timescale, i. e. , the timescale of relaxation in [m0] is approximately one thousand times that of kc and kp. Consequently, the system reaches the true steady state, i. e. , the other components relax into the steady-state values in the timescale of kc and kp. When we change the ratio between [c] and [c*], which is considered the change in csum, a similar type of behavior is observed (S1 Fig). Although the concentrations of molecules without m0 change when the ratio changes, the system can reach the quasi-steady state quickly. Then, [m0] changes slowly and the other components subsequently change rapidly. The timescale of the relaxation depends on the input strength. We defined the relaxation time τ as the time when the sum of differences of the concentrations of all molecules from the previous time point falls below the threshold (10−7). When we change kin, τ diverges at the point where k in = k in th (Fig 2B, inset). The relaxation time is proportional to the inverse of the difference between kin and k in th (Fig 2B) so that the relaxation time critically slows down around the critical point. Note that although the relaxation time is also critically prolonged in the vicinity of saddle-node bifurcation, the exponent is 1/2 [21] and differs from the obtained result. When the relaxation time slows down, the relaxation manner of [m0] becomes linear with time rather than exponential (S2 Fig). This can be observed in both cases where kin and csum change in the system (Fig 2A and S1 Fig). If both the active carrier-producing and -consuming reactions are reversible, similar slow dynamics appears, while the switch from the quasi-steady state to the true steady state is smooth for small kout values (see S3 Fig and section 2 in S1 Text). This may be due to inhibition of the active carrier-producing reaction by binding of m2 to c*. Therefore, the slow relaxation is considered to be a general behavior of the CCC in response to different environmental changes. To elucidate the process of the slow relaxation, we reduced the number of arguments in our model by using conserved quantities. Although our model had five ODEs even when the association and dissociation reactions are eliminated adiabatically, [m2] is involved only in the equation for [m2] and not in other equations. Moreover, two equations can be eliminated because of two independent conserved quantities. Hence, the model can be reduced to the following two ODEs: d [ m 0 ] d t = k in - k c [ m 0 ] [ c ] K 0 + [ c ] - k leak [ m 0 ], (2a) d [ m 1 ] d t = k c [ m 0 ] [ c ] K 0 + [ c ] - k p [ m 1 ] [ c * ] K 1 + [ c * ], (2b) [ c ] = - β + { β 2 + 4 K 0 (c sum - [ m 1 ]) } 1 / 2 2, (2c) [ c * ] = - γ + { γ 2 + 4 K 1 (c pool - c sum + [ m 1 ]) } 1 / 2 2, (2d) where β = [m0] + K0 − csum + [m1] and γ = K1 + csum − cpool. The nullclines for [m0] and [m1] are presented in Fig 3, in the case of kleak = 0. By altering kin, the vertical position of the nullcline for [m0] is changed as well. The CCC behavior changes drastically in the vicinity of the critical point. When kin is smaller than k in th, a stable fixed point appears (Fig 3A), and when kin is close to k in th, relaxation to the fixed point is critically slowed down due to the approach of two nullclines. In this case, a slow manifold is located on the nullcline for [m1], which is nearly parallel to the nullcline for [m0]. In the region where two nullclines are parallel, the relaxation speed across the slow manifold is proportional to the distance from the nullcline for [m0]; i. e. , the speed of change in [m0] is given as a constant that is proportional to k in th - k in. Hence, the degradation of m0 is the rate-limiting process when the distance between the nullclines is sufficiently small, and the relaxation time is proportional to the inverse of k in th - k in (Fig 2B), while the relaxation manner of m0 depends linearly on time (S2 Fig). However, when kin is larger than k in th, nullclines do not cross anywhere and [m0] diverges (Fig 3C). In particular, if kin is equal to k in th, the points on the overlapped nullclines become neutrally stable fixed points (Fig 3B). When kin is fixed and csum is changed, both of the nullclines move. Moreover, in this case, the nullclines approach each other for the critical csum value (see S4 Fig), and the slow relaxation occurs. Here, k in th can be obtained analytically for kleak = 0 (see section 3 in S1 Text). k in th = k c k p 2 (k c + k p) { c pool + c sum + α + (c pool - c sum + α) 2 + 4 c sum α }. (3) where α = kc K1/ (kc + kp). For the limit of K1 → 0, i. e. , m1 can perfectly bind to c* and never dissociate, Eq (3) becomes: k in th = { k c k p c sum (k c + k p) (c pool > c sum), k c k p c pool (k c + k p) (c pool < c sum). (4) Eq (4) represents the maximal capacity of the flux of CCC. Therefore, when the influx exceeds the capacity, the CCC can become jammed, which leads to the appearance of slow dynamics. If kleak is not zero, the nullcline for [m0] is tilted, and a fixed point is obtained for large kin (S5 Fig). For a system with a finite kleak value, k in th cannot be defined. However, in the parameter range where kin is close to the k in th for the system with no kleak, the fixed point value changes drastically and slow relaxation appears, which cannot be represented by bifurcation. For organisms, robustness against small fluctuations in nutrient uptake is important for maintenance of the intracellular environment. However, in the case of a sudden decrease in the concentration of nutrients, which may lead to starvation, stress-resistant systems should be activated [22,23]. Here, we analyzed the frequency responses of the CCC using a cyclic nutrient uptake rate kin (t) with different amplitudes (Fig 4). The CCC demonstrates the low-pass filter characteristics with a sharp cut-off frequency for both weak and strong inputs (Fig 4A). For the case when kin (t) is higher than k in th, m0 accumulates and the concentrations of the others does not change over time. For the case when kin (t) is lower than k in th, m0 decreases slowly, while the concentrations of the other components remains the same. Here, m0 concentration represents a buffer, and the sharp cut-off frequency can be achieved. The timescale of the decrease depends on the distance between kin and k in th, and thus the filter characteristics depend on the amplitude of inputs (Fig 4B). For inputs with a smaller amplitude, the cut-off frequency becomes lower, but when a larger amplitude of inputs is used, the CCC demonstrates the response against the higher frequency of inputs. For larger concentrations of m0, the rates of enzymatic reactions, i. e. , the active coenzyme-consuming and -producing reactions, are considered to be constant due to saturation. Therefore, the nullcline for [m1] can be considered to be nearly constant, while that for [m0] can be regarded as a linear equation of [m0] due to a leak term (S5 Fig). When two nullclines are close, the flow rate along the nullclines can be approximated as the distance between the two nullclines, which is slower than the flow rate approaching the nullclines. Accordingly, the two-dimensional dynamics can be reduced into one-dimensional dynamics along the nullcline for [m1]. We analytically obtained the frequency response of the dynamics. When kleak is not as large, the two-dimensional dynamics (Eqs (2a) and (2b) ) can be reduced into one-dimensional dynamics of [m0] (see section 4 in S1 Text). d [ m 0 ] d t = - k leak k c [ m 0 ] - k c c sum - k in (t) k c + k c c sum (k c + k p). (5) For kin (t) as a sinusoidal k in (t) = A in cos (2 π f t) + k in 0, when [m0] is saturated, the cut-off frequency is given as: 2 π f = k leak k c A in k c k in 0 - (k c + k leak) k in th - k leak k c. (6) The estimated cut-off frequency fits well with the simulation result (red dashed line in Fig 4B), suggesting that the complex dynamics can be reduced to one-dimensional dynamics due to the conserved quantities, and that the necessity for robustness against external fluctuation is determined by the conditions underlying this saturation. Note that Eq (5) is independent of a form of function of the influx rate. Hence, if the influx rate is regulated by downstream products, the dynamics of the CCC are reduced into one-dimensional dynamics. Although the filter characteristics (Eq (6) ) depend on the regulation, robustness by the slow dynamics will always be achieved. To investigate the robustness of the CCC against the intrinsic noise caused by the stochasticity of biochemical reactions, we calculated the stochastic dynamics of the original full model with consideration of complex formation using the Gillespie algorithm [24], which can generate a statistically possible trajectory of the solution of stochastic equations and is often used for modeling biochemical reactions. We calculated the long trajectories of the solution and several statistics from the given trajectories. When the cpool is larger than the critical value, the average number of m1 is almost constant for different cpool values. In this condition, the Fano factor, which is the ratio of the variance to the average, is approximately 1, as in similar non-catalytic reactions. However, below the critical point where the flux is limited by the concentration of the carrier, i. e. , the first metabolite is saturated, the Fano factor of the number of m1 decreases below 1 (Fig 5A). This suggests that the carrier cycling can reduce the variance of the concentration of intermediate metabolites. To investigate the mechanisms underlying the decrease in the intrinsic noise, we analyzed the probabilistic dynamics of the number of m1 molecules, n (see section 5 in S1 Text). Under the not-saturated condition, i. e. , when the m0 concentration is lower than the maximal coenzyme concentration cmax, which is the same as csum when csum < cpool, only the consumption rate of m0 is proportional to n but the production rate is not. Therefore, the steady-state distribution of the number of m1 is given as a Poisson distribution [25] in the limit of K0 → 0 and K1 → 0; i. e. , the metabolites can bind to coenzymes perfectly. Thus, the Fano factor becomes 1, which is similar to the previously reported condition [26]. In the saturated condition, i. e. , the m0 concentration is higher than cmax, the production rate becomes kc (cmax − n) due to conservation of the coenzyme, while the consumption rate remains the same as above. Here, both the production and consumption rates are proportional to n, which is considered the feedback-regulated production through the conserved quantity. The steady-state distribution represents the binomial distribution [25] and the Fano factor is given as: σ 2 < n > = 1 - k c k c + k p. (7) Therefore, the fluctuation is reduced depending on the kc and kp values (Fig 5B), and the recycling of carriers improves the signal-to-noise ratio by feedback regulation via the moiety conservation. This effect does not depend on the concentration of the coenzyme as long as the metabolite is saturated before the coenzyme consumption step. The feedback can also reduce the fluctuation in the active and inactive carrier concentrations. Note that when both the production and consumption rates of m1 obey the Michaelis-Menten reaction with different coenzymes, the number of m0 shows a random walk in the range of 0 to ∞, or becomes zero, or diverges, depending on the parameter settings. In any case, the Fano factor never falls below 1 (see section 5 in S1 Text).
Here, we proposed a minimal model of metabolic systems that includes recycling of the moiety-conserved carriers with complex formation between carriers and other metabolites. In transient dynamics, the conserved quantities constrain a dimension of orbits moving on to the steady state, similar to the steady-state solution space, and dynamics in the restricted dimension demonstrate various phenomena. These effects can be summarized by two properties: 1) jamming and 2) feedback, which are followed by the slow relaxation and robustness against internal and external fluctuations. We demonstrated that the relaxation dynamics in the CCC are decelerated by the jamming when the nutrient uptake rate is close to the capacity of the cascade. Such slow relaxation to the steady state has also been discussed in other enzymatic networks [17]. From the viewpoint of dynamical systems theory, this jamming is due to a restriction of the phase space and the closure of nullclines by the conservation of carrier concentration. In this situation, the concentrations of internal metabolites, the final product, and active and inactive carriers are almost constant and can be drastically different from the concentrations in the steady state. Hence, if a metabolic system seems to be in a steady state on a short timescale after an environmental change, this system will not always be in the true steady state but rather in a quasi-stable state. Consequently, the previous theoretical studies considering steady metabolic states were not able to analyze the cellular metabolism in these quasi-stable states. Such jamming in the metabolic flow is likely to be observed in several metabolic systems, including the glycolytic cascades. In fact, some mutants of the budding yeast exhibit growth defects, which were explained due to the abnormal accumulation of intermediate metabolites based on computer simulations [18]. Such time evolution is consistent with our results caused by the jamming, while we further uncovered the origin of this phenomenon using dynamical systems theory. This study suggests the possibility that the native metabolic cascade resides in the vicinity of the critical point and can be easily jammed by genetic perturbations, which will be validated in future studies. Furthermore, we suggest that the jamming mechanism may represent the mediator between molecular- and organism-level timescales [27]; i. e. , the jamming slows down the fast enzymatic turnover and may determine the timescale of physiological behaviors. The slow relaxation can help organisms maintain their cellular condition against changes in the nutrient condition and will work as a memory of past environmental change. For example, the slow timescale might be related to a slow response of the relaxed strain against carbon and amino acid starvation [22,23]. We have demonstrated the robustness of metabolic systems against external fluctuation, which opens the door for further theoretical studies of the quasi-steady state of metabolic dynamics. The concentrations of some metabolites and carriers are maintained throughout small fluctuations in the influx rates. At the same time, the CCC is responsible for the large alterations in the influx rate. Therefore, the CCC shows both robustness against small fluctuations and responsiveness to large fluctuations; i. e. , the CCC does not respond to small and short-term changes in the nutrient uptake rate, but it can rapidly respond to larger alterations in nutrient uptake rates. This robustness and responsiveness are most likely due to the wide range of the quasi-steady states in the phase space and the input-dependent timescales of relaxation. These two properties depend on the non-linear dynamics owing to the complex formation that is not in the vicinity of a fixed point, which has not been investigated using a constraint-based approach or with steady-state analysis of the previous model. A relationship between robustness and responsiveness has also been studied in other systems [28] and will be further investigated in various biological systems. To measure quasi-steady state characteristics, the dynamics of metabolites should be measured non-invasively for a long time; however, the real-time measurement techniques of metabolites are insufficient at present. Thus, we expect that the continuous development of new measurement techniques will help to validate our predictions experimentally. We demonstrated that an internal metabolite regulates its production via feedback mechanisms by the moiety conservation. The noise in the concentrations of metabolites and carrier can be reduced drastically. This feedback mechanism is important for stabilization of the metabolite and carrier concentrations at the finite constant values during the slow relaxation. When this feedback is missing, this noise cannot be reduced microscopically. Therefore, the concentration of the internal molecule should decrease to zero or diverge in the saturated condition (S6 Fig and section 5 in S1 Text). This suggests that the feedback mechanism involving the carrier cycling may be responsible for the existence of the quasi-steady states. There are several important parameters contributing to the robustness against both external and internal fluctuations. In particular, kc, the turnover rate of an enzyme in the active carrier-consuming reaction, and kp, that in the active carrier-producing reaction, are important. We estimated these turnover rates using data from previously published studies: the glycolytic pathway of Escherichia coli during continuous aerobic cultivation [3,29] and the lactic acid fermentation pathway of Lactococcus lactis during anaerobic cultivation with an 80% lactic acid yield [30–32]. In both pathways, the enzyme levels necessary for each reaction are considered to be sufficient and not a rate-limiting factor, which we assumed in our model, and allowed us to consider the kinetics of substrate and carrier alterations subsequently. Based on our estimation, kc and kp in the glycolytic pathway are given as 50 s−1 and 12 s−1, respectively, while the Fano factor of the internal metabolite is estimated at nearly 0. 2. For the lactic acid fermentation pathway, kc and kp are given as 1. 2 s−1 and 4. 9 s−1, respectively, and the Fano factor is estimated at 0. 85 (see Table 1). Our estimations are too simplified to allow for a quantitative discussion about the cellular metabolic process, because actual metabolic processes are more complicated than represented by our model. However, the properties observed with our model are preserved even if the details of the model change. Indeed, if two CCCs are coupled through a common carrier pool, the slow relaxation to the steady state is observed (see S7 and S8 Figs, and section 6 in S1 Text). This suggests that the moiety conservation pool can underlie the jamming and the feedback processes that determine the dynamics of more complicated metabolic networks than the CCC. If there are some branches in the CCC, our results would still be qualitatively reproduced, although some quantities would no longer be conserved but rather given as quasi-steady-state values. It has been reported that such branches will change the steady-state characteristics [11], e. g. , the bistability and hysteresis, and the slow dynamics might appear in transient dynamics approaching each stable fixed point. We expect that further non-trivial dynamic phenomena will be observed in more complicated metabolic networks with moiety conservation via the jamming and feedback. To investigate the effects of conserved carriers in more complicated metabolic networks quantitatively, both theoretical and experimental investigations are required in the future.
Our analyses suggest that enzyme turnover rates in metabolic pathways are the essential contributors to the robustness of metabolic dynamics. We estimated the parameters described in actual metabolic pathways, i. e. , the glycolytic pathway of Escherichia coli during continuous aerobic cultivation [29] and the lactic acid fermentation pathway of Lactococcus lactis during anaerobic cultivation with an 80% lactic acid yield [30], and calculated the effects of carrier cycling. We approximated the kinetics of each reaction using the Michaelis-Menten equation. Hence, the concentrations of a substrate and carrier, the dissociation constant, and the speed of flux are necessary for the estimation of kc and kp. For kc of phosphofructokinase (PKF) in the glycolytic pathway, the concentrations of fructose 6-phosphate (F6P) as a substrate and ATP as a carrier used are 3. 21 × 10−2 mM and 8. 28 × 10−1 mM, respectively [29]. The dissociation constant is given as 0. 16 mM [3]. The flux was estimated from the glucose uptake rate, which is 1. 8 mM s−1 under the same conditions described previously [29]. However, the entire amount of glucose is not catalyzed by PKF, and 20–40% of glucose is considered to be involved in the pentose phosphorylation pathway [29]. Here, we considered this leak to be 20%, and the flux is estimated as 1. 8 × 0. 8 = 1. 44 mM s−1. Subsequently, kc is estimated to be approximately 50 s−1, using the Michaelis-Menten form. For kp of pyruvate kinase (PK), the concentrations of phosphoenolpylvate (PEP) and ADP are 1. 49 × 10−1 mM and 9. 65 × 10−1 mM, respectively [29]. The dissociation constant is 0. 26 mM [3]. We assumed that the flux of PK is similar to that of PKF, and kp is estimated to be approximately 12 s−1. From the estimated kc and kp, σ2/<n> of the internal metabolite is determined to be 0. 2. In the same manner, kc and kp in the lactic acid fermentation pathway can be estimated. For kc of glyceraldehyde phosphate dehydrogenase (GAPDH), the concentrations of glyceraldehyde 3-phosphate (GAP) and NAD+ are 6. 0 mM and 8. 4 mM, respectively [30]. The dissociation constant is 0. 2 mM [31] and the flux is 5 mM s−1 [30]. Therefore, kc is estimated to be 0. 85 s−1. For kp of lactate dehydrogenase (LDH), the concentration of pyruvate and NADH are estimated as 1 mM and 0. 7 mM, respectively [30]. The dissociation constant is 0. 08 mM [32], while the flux is given as 4. 4 mM s−1 [30]. Therefore, kp is estimated to be 4. 9 s−1, and the σ2/<n> of the internal metabolite is determined to be 0. 85.
Although a metabolic shift is essential for the adaptation of cells or organisms to environmental changes, the transient behaviors of metabolic systems are poorly understood. When describing the time development of metabolic systems, there are several conserved quantities to consider due to balances of metabolic reactions, e. g. , the cycling of coenzymes. Such conserved quantities limit the possible changes in the metabolic state and can generate non-trivial dynamical behaviors. We here propose a minimal motif of metabolic reactions that includes coenzyme recycling to investigate the effect of conserved quantities on metabolic dynamics. We demonstrate that the dynamics with this motif intrinsically show slow relaxation to the steady state after environmental changes. Moreover, this motif can maintain robustness against external and internal fluctuations owing to the conservation of coenzymes. Overall, these results suggest that the complex metabolic dynamics generated by coenzyme recycling are beneficial to organisms.
relaxation time
coenzymes
enzymes
metabolic networks
enzymology
systems science
mathematics
metabolites
network analysis
enzyme metabolism
enzyme chemistry
relaxation (physics)
computer and information sciences
proteins
metabolic pathways
dynamical systems
physics
biochemistry
biology and life sciences
cofactors (biochemistry)
physical sciences
metabolism
Metabolic dynamics restricted by conserved carriers: Jamming and feedback