Metabolic dynamics restricted by conserved carriers: Jamming and feedback

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.


Original full model
Ordinary differential equations (ODEs) of the simplest model of the carrier cycling cascade (CCC) are given as equations of mass action kinetics, as follows: where m i represents the i-th metabolite, c and c * are the active and inactive carriers, respectively, and [x] denotes the concentration of x. m 0 and m 1 can form complexes with the active and inactive carriers as cm 0 and c * m 1 with the association rates k b0 and k b1 , respectively, and these complexes can decompose with the dissociation rates k d0 and k d1 , respectively. Active carriers are consumed with the rate k c when m 1 is transformed from m 0 and are produced with rate k p when m 2 is transformed from m 1 . m 0 is supplied and diluted with rates k in and k leak , respectively, and m 2 is diluted with rate k out . If we assume that k c , k p ≪ k bi , k di , these equations can be reduced to five ODEs, and we calculated the five-variable model. Here, we use the parameters presented in S1

Reversible model
ODEs of the reversible model are given as follows: where k r1 and k r2 are the speed of reactions from m 1 to m 0 and from m 2 to m 1 , respectively. K r1 and K r2 are the dissociation constants between m 1 and c * and between m 2 and c, respectively. Here, we use the parameters presented in S2 Table. TableS. 2: Parameters used in the reversible model.

Parameter Value
The reduced CCC model is given as: Here, we obtained k th in analytically. When k leak is zero, a nullcline for [m 0 ] is analytically calculated as: This nullcline converges into [m 1 ] = c sum − k in /k c for [m 0 ] → ∞. Although a nullcline for [m 1 ] becomes too complicated, analysis of the limit of [m 0 ] → ∞ is helpful. In these limits, the first term of the right side in Eq.(3b) becomes because the maximum speed can be obtained as a multiplication of the turnover rate of the enzyme and the active carrier concentration. Thus, the nullcline in the limit of [m 0 ] → ∞ is obtained as: where α = k c K 1 /(k c + k p ). Therefore, from Eqs.(4) and (5), k th in is obtained as: For the limit of K 1 → 0, i.e., in the case where m 1 can perfectly bind to c * and never unbind, Eq.(6) becomes:

Analytical calculation of the frequency response of the CCC model
To obtain the frequency response analytically, we calculated the nullclines using the large [m 0 ] limit. Using the limit of [m 0 ] → ∞, the speed of the active coenzyme-consuming reaction becomes respectively. A nullcline for [m 1 ] in the large [m 0 ] and small K 1 limit is calculated from Eq.(5) as Therefore, from Eqs. (9) and (10), the fixed point value of [m 0 ] is given as: When k leak is small, two nullclines are close to the fixed point, and then the speed of change in [m 0 ] on the nullcline for [m 1 ] is much slower than that approaching the nullcline and is proportional to the distance between two nullclines. Following this, the two-dimensional dynamics (Eqs.(3a) and (3b)) can be reduced into one-dimensional dynamics of [m 0 ] around the fixed point as When k in (t) is given as a sinusoidal function k in (t) = A in cos(2πf t) + k 0 in , [m 0 ](t) at the steady state is calculated as Here, the above approximations are feasible when [m 0 ] is higher than [m 0 ] th = k th in /k c due to saturation, which is the maximal m 0 concentration processed in the first reaction per unit of time. When [m 0 ] becomes lower than [m 0 ] th due to changes in the influx rate, the fixed point is drastically changed and the concentrations of [m 1 ], [c] t , and [c * ] t are altered. Therefore, the cut-off frequency is given as the frequency at which the minimal value of Eq.(13) is the same as [m 0 ] th , as follows: The estimated cut-off frequency fits well with the result of our simulation (Fig. 4B in the main text), suggesting that the complex dynamics can be reduced into one-dimensional dynamics due to saturation, and the need for robustness against external fluctuation is determined by the conditions underlying this saturation.

Internal fluctuation of the metabolite concentration
We considered the stochastic dynamics of the number of m 1 molecules, n. First, we calculated the steady-state distribution of n in a non-saturated condition, i.e., the m 0 concentration is lower than the maximal coenzyme concentration c max . Here, we consider the following conditions: c sum < c pool and c max is the same as c sum . Using the limits of K 0 → 0 and K 1 → 0, i.e., the metabolites bind to coenzymes perfectly, the production rate of m 1 becomes k c [m 0 ], and the consumption rate of m 1 becomes k p n when the number of m 1 molecule is n. Here, only the consumption rate is proportional to n, while the production rate is not, so that the master equation can be given as follows: where p n is the probability that the number of m 1 molecule is n. The steadystate distribution is the Poisson distribution: Therefore, both < n > and σ 2 are given as k c [m 0 ]/k p and σ 2 / < n > becomes 1, which is similar to a previously reported condition [1]. Under the saturated condition, i.e., the m 0 concentration is higher than c max , because of coenzyme conservation, the production rate becomes k c (c max − n) while the consumption rate remains the same as the non-saturated condition. Here, both the production and consumption rates are proportional to n. The master equation is thus given as follows: The steady-state distribution is the binomial distribution: where K is k c /k p . The moment-generating function is given as: and then, Therefore, the fluctuation can be reduced depending on the k c and k p values (see Fig. 5B in the main text). The carrier cycling can improve the signalto-noise ratio by feedback regulation through the conserved concentration of a carrier, and this effect does not depend on the concentration of a coenzyme as long as the metabolite is saturated before the coenzyme-consuming step. This suggests that the feedback in the CCC model can reduce the fluctuations of the active and inactive carrier concentrations because of the conserved quantities between [c] t and [m 1 ].
To investigate the differences in the CCC from ordinal Michaelis-Menten type reactions, we considered a double Michaelis-Menten model; i.e., reactions from m 0 to m 1 and m 1 to m 2 are catalyzed by different catalysts. Under the non-saturated condition, the behavior of the double Michaelis-Menten model is similar to that of the CCC model, and the Fano factor is approximately 1. However, under the saturated condition, the concentration of m 1 shows a random walk (Fig. S6) when a similar parameter set as used for the CCC model was used, as demonstrated in Fig. 5A in the main text. As shown for the saturated condition, the Fano factor of the m 1 concentration never decreases in the double Michaelis-Menten model. Furthermore, in the double Michaelis-Menten model, under the nonsaturated condition, the master equation is the same as that of the CCC model, and a steady-state distribution is given as the Poisson distribution.
However, under the saturated condition, the master equation is given as where c 1 and c 2 represent the catalysis for the first and second step reactions, respectively. The steady-state distribution is given as Here, if k c c 1 /k p c 2 is higher than 1, the time evolution of the number of m 1 is governed by the asymmetric random walk and the number of m 1 will diverge. If k c c 1 /k p c 2 is 1, the time evolution of the number of m 1 is governed by the symmetric random walk from 0 to ∞, as shown in S6 Fig. If k c c 1 /k p c 2 is lower than 1, the steady-state distribution can be given as a geometric distribution: The average and variance are: Therefore, in the double Michaelis-Menten model, σ 2 / < n > should be higher than 1 in all cases.

Coupled carrier cycling cascade (CCCC) model
We considered the condition in which two carrier cycling cascades are coupled through a common coenzyme pool. Association and dissociation reactions between the coenzyme and substrates are eliminated adiabatically, as in the single cascade model, and the ODEs can be given as: where i represents an index of cascade numbers. Although we set N to 2, i.e., two CCCs are coupled (S7 Fig), the obtained results do not change for higher N when the following conditions are satisfied. The CCCC shows a qualitatively similar response to the environmental changes as the CCC when one cascade is a major pathway. Here, we considered two examples: 1. c and m 2 0 with a larger dissociation constant than c and m 1 0 . 2. The influx rate (leak rate) of m 2 0 is smaller (larger) than that of m 1 0 . In both cases, cascade 1 is the major cascade, and the CCCC displays responses to environmental changes in a similar way as observed for the CCC (S8 Fig). Therefore, our findings should be applicable to more complicated metabolic networks as well.
For these, we used the parameter set presented in S3 Table, unless otherwise specified.