Dynamics of cosubstrate pools can constrain and regulate metabolic fluxes
Abstract
Cycling of cosubstrates, whereby a metabolite is converted among alternate forms via different reactions, is ubiquitous in metabolism. Several cycled cosubstrates are well known as energy and electron carriers (e.g. ATP and NAD(P)H), but there are also other metabolites that act as cycled cosubstrates in different parts of central metabolism. Here, we develop a mathematical framework to analyse the effect of cosubstrate cycling on metabolic flux. In the cases of a single reaction and linear pathways, we find that cosubstrate cycling imposes an additional flux limit on a reaction, distinct to the limit imposed by the kinetics of the primary enzyme catalysing that reaction. Using analytical methods, we show that this additional limit is a function of the total pool size and turnover rate of the cycled cosubstrate. Expanding from this insight and using simulations, we show that regulation of these two parameters can allow regulation of flux dynamics in branched and coupled pathways. To support these theoretical insights, we analysed existing flux measurements and enzyme levels from the central carbon metabolism and identified several reactions that could be limited by the dynamics of cosubstrate cycling. We discuss how the limitations imposed by cosubstrate cycling provide experimentally testable hypotheses on specific metabolic phenotypes. We conclude that measuring and controlling cosubstrate dynamics is crucial for understanding and engineering metabolic fluxes in cells.
Editor's evaluation
This manuscript presents an important mathematical analysis of metabolic "cosubstrates" and how their cycling can affect metabolic fluxes. Through mathematical analysis of simple network motifs, it shows the impact of cosubstrate cycling on constraining metabolic fluxes. The combination of mathematical modeling and comparisons with existing data from previous studies offers convincing support for the potential biological relevance of cosubstrate cycling. The work will be of interest to researchers who study microbial metabolism and metabolic engineering.
https://doi.org/10.7554/eLife.84379.sa0eLife digest
Metabolism powers individual cells and ultimately the body. It comprises a sequence of chemical reactions that cells use to break down substances and generate energy. These reactions are catalyzed by enzymes, which are proteins that speed up the rate of the reaction. Many reactions also involve cosubstrates, which are themselves transformed by individual reactions but are eventually converted back into their original form in a series of steps. This process is known as cosubstrate cycling.
Scientists have long been interested in understanding what controls the rate at which metabolic reactions and metabolic pathways convert a substance into a final product. This is a difficult subject to study because of the complexity of the metabolic pathways, with their branched, linear or coupled structures. In the past, researchers have looked at the influence of enzymes on the rate of a metabolic pathway, but less has been known about the effect of cosubstrate cycling.
To find out more, West, Delattre et al. developed a series of mathematical models to describe different types of metabolic pathways in terms of the number of metabolites that enter and leave it, including the influence of cosubstrates.
They found that cosubstrate cycling, when involved in a metabolic reaction, limits the speed with which the reaction happens. This is distinct from the limit that enzymes impose on the speed of the reaction. It depends on the total amount of cosubstrates in the cell: changing the number of cosubstrates in the cell influences the speed at which the metabolic reaction takes place.
This study has increased our understanding of how metabolic pathways work, and what controls the speed at which reactions take place. It opens up a new potential method for explaining how cells control metabolic reaction rates and how metabolic substrates can be directed across different pathways. This research is likely to inspire future research into the influence of cosubstrates in different cell types and conditions.
Introduction
Dynamics of cell metabolism directly influences individual and populationlevel cellular responses. Examples include metabolic oscillations underpinning the cell cycle (Papagiannakis et al., 2017; Murray et al., 2007) and metabolic shifts from respiration to fermentation, as observed in cancer phenotypes (Warburg, 1956; DiazRuiz et al., 2009; CarmonaFontaine et al., 2013) and celltocell crossfeeding (Ponomarova et al., 2017; Campbell et al., 2015; Großkopf et al., 2016). Predicting or conceptualising these physiological responses using dynamical models is difficult due to the large size and high connectivity of cellular metabolism. Despite this complexity, however, it is possible that cellular metabolism features certain ‘design principles’ that determine the overall dynamics. There is ongoing interest in finding such simplifying principles.
A key concept for understanding the dynamics of any metabolic system is that of ‘reaction flux’, which is a measure of the rate of biochemical conversion in a given reaction. To identify possible limitations on reaction fluxes, early studies focused on linear pathways involving ATP production and studied their dynamics under the optimality assumption of maximisation of overall pathway flux under limited enzyme levels available to the pathway (Heinrich et al., 1991). The resulting theory predicted a tradeoff between pathway flux vs. yield (i.e. rate of ATP generation vs. amount of ATP generated per metabolite consumed by the pathway) in linear pathways (Heinrich and Hoffmann, 1991). This theory is subsequently used to explain the emergence of different metabolic phenotypes (Pfeiffer et al., 2001). In related studies, models pertaining to flux optimisation and enzyme levels being a key limitation are used to explain the structure of different metabolic pathways (Flamholz et al., 2013), and the metabolic shifting from respiration to fermentative pathways under increasing glycolysis rates (Großkopf et al., 2016; Basan et al., 2015; Majewski and Domach, 1990). There are, however, increasing number of studies suggesting that enzyme levels alone might not be sufficient to explain observed flux levels. For example, it was shown that the maximal value of the apparent activities (${k}_{app}^{max}$) of an enzyme, derived using measured enzyme levels and fluxes under different conditions, was a good estimate for the specific activity of that enzyme in vitro (${k}_{cat}$) (Davidi et al., 2016). However, individual estimates from each condition (i.e. individual ${k}_{app}$ values) were commonly lower than the specific activity – suggesting that the flux is limited by something other than enzyme levels under those conditions. Other studies have shown that metabolic flux changes, caused by perturbations in media conditions, are not explained solely by changes in expression levels of enzymes (Chubukov et al., 2013; Gerosa et al., 2015).
Another conceptual framework emphasized the importance of cyclic reaction motifs, particularly those involving socalled cosubstrate pairs, such as ATP / ADP or NAD(P)H / NAD(P)${}^{+}$, as a key to understanding metabolic system dynamics (Reich and Sel‘kov, 1981). This framework is linked to the idea of considering the supply and demand structures around specific metabolites as regulatory blocks within metabolism (Hofmeyr and CornishBowden, 2000). For example, the total pool of ATP and its derivates (the ‘energy charge’) is suggested as a key determinant of physiological cell states (Atkinson, 1968). Inspired by these ideas, theoretical studies have shown that metabolic systems featuring metabolite cycling together with allosteric regulation can introduce switchlike and bistable dynamics (Okamoto and Hayashi, 1983; Hervagault and Cimino, 1989), and that metabolite cycling motifs introduce total cosubstrate level as an additional control element in metabolic control analysis (Hofmeyr et al., 1986; Sauro, 1994). Specific analyses of ATP cycling in the glycolysis pathway, sometimes referred to as a ‘turbodesign’, and metabolite cycling with autocatalysis, as seen for example in the glyoxylate cycle, have shown that these features constrain pathway fluxes (Koebmann et al., 2002; Teusink et al., 1998; van Heerden et al., 2014; Hatakeyama and Furusawa, 2017; Barenholz et al., 2017; Kurata, 2019). Taken together, these studies indicate that metabolite cycling, in general, and cosubstrate cycling specifically, could provide a key ‘design feature’ in cell metabolism, imposing certain constraints or dynamical properties to it.
Towards better understanding the role of cosubstrate cycling in cell metabolism dynamics, we undertook here an analytical and simulationbased mathematical study together with analyses of measured fluxes. We created models of enzymatic reaction systems featuring cosubstrate cycling, abstracted from real metabolic systems such as glycolysis, nitrogenassimilation, and central carbon metabolism. We found that cosubstrate cycling introduces a fundamental constraint on reaction flux. In the case of single reaction and short linear pathways, we were able to derive a mathematical expression of the constraint, showing that it relates to the pool size and turnover rate of the cosubstrate. Analysing measured fluxes, we find that several of the cosubstrate featuring reactions in central carbon metabolism carry lower fluxes than expected from the kinetics of their primary enzymes, suggesting that these reactions might be limited by cosubstrate cycling. In addition to its possible constraining role, we show that cosubstrate cycling can also act as a regulatory element, where control of cosubstrate pool size can allow control of flux dynamics across connected or branching pathways. Together, these findings show that cosubstrate cycling can act both as a constraint and a regulatory element in cellular metabolism. The resulting theory provides testable hypotheses on how to manipulate metabolic fluxes and cell physiology through the control of cosubstrate pool sizes and turnover dynamics, and can be expanded to explain dynamic measurements of metabolite concentrations in different perturbation experiments.
Results
Cosubstrate cycling represents a ubiquitous motif in metabolism with cosubstrate pools acting as ‘conserved moieties’
Certain pairs of metabolites can be interconverted via different reactions in the cell, thereby resulting in their ‘cycling’. This cycling creates interconnections within metabolism, spanning either multiple reactions in a single, linear pathway, or multiple pathways that are independent or are branching from common metabolites. For example, in glycolysis, ATP is consumed in reactions mediated by the enzymes glucose hexokinase and phosphofructokinase, and is produced by the downstream reactions mediated by phosphoglycerate and pyruvate kinase (Appendix 1—figure 1A). In the nitrogen assimilation pathway, the NAD^{+} / NADH pair is cycled by the enzymes glutamine oxoglutarate aminotransferase and glutamate dehydrogenase (Appendix Dynamics of cosubstrate pools can constrain and regulate metabolic fluxes  Appendix 1—figure 1B). Many other cycling motifs can be identified, involving either metabolites from the central carbon metabolism or metabolites that are usually referred to as cosubstrates. Examples for the latter include NADPH, FADH2, GTP, and AcetylCoA and their corresponding alternate forms, while examples for the former include the tetrahydrofolate (THF) / 5,10MethyleneTHF and glutamate / $\mathrm{\alpha}$ketoglutarate (akg) pairs involved in onecarbon transfer and in amino acid biosynthesis pathways, respectively (Appendix 1—figure 1C and D). For some of these metabolites, their cycling can connect many reactions in the metabolic network. Taking ATP (NADH) as an example, there are 265 (118) and 833 (601) reactions linked to the cycling of this metabolite in the genomescale metabolic models of Escherichia coli and human respectively models iJO1366 (Orth et al., 2011) and Recon3d (Brunk et al., 2018).
We notice here that many of the cosubstrate involving cycling reactions can be abstracted as a simplified motif as shown in (Figure 1A). This abstract representation highlights the fact that the total poolsize involving all the different forms of a cycled metabolite can become a conserved quantity. This would be the case even when we consider biosynthesis or environmental uptake of cosubstrates, as the total concentration of a cycled metabolite across its different forms at steady state would then be given by a constant defined by the ratio of the influx and outflux rates (see Appendices 2 and 3). In other words, the cycled metabolite would become a ‘conserved moiety’ for the rest of the metabolic system and can have a constant ‘pool size’. Supporting this, temporal measurement of specific cosubstrate pool sizes shows that ATP and GTP pools are constant under stable metabolic conditions, but can rapidly change in response to external perturbations, possibly through interconversions among pools rather than through biosynthesis (Walther et al., 2010).
Cosubstrate cycling introduces a limitation on reaction flux
To explore the effect of cosubstrate cycling on pathway fluxes, we first consider a didactic case of a single reaction. This reaction converts an arbitrary metabolite M_{0} to M_{1} and involves cosubstrate cycling (Figure 1A). For cosubstrate cycling, we consider additional ‘background’ enzymatic reactions that are independent of M_{0} and can also convert the cosubstrate (denoted ${E}_{a}$ on Figure 1A). We use either irreversible or reversible enzyme dynamics to build an ordinary differential equation (ODE) kinetic model for this reaction system and solve for its steady states analytically (see Methods and Appendix 3). In the case of using irreversible enzyme kinetics, we obtain that the steady state concentration of the two metabolites, M_{0} and M_{1} (denoted as m_{0} and m_{1}) are given by:
where k_{in} and k_{out} denote the rate of influx of M_{0}, and outflux of M_{1}, either inandout of the cell or from other pathways, and ${A}_{\mathrm{tot}}$ denotes the total pool size of the cycled metabolite (with the different forms of the cycled metabolite indicated as A_{0} and A_{1} in Figure 1A). The parameters ${V}_{\mathrm{max},{E}_{0}}$ and ${V}_{\mathrm{max},{E}_{a}}$ are the maximal rates (i.e. ${V}_{\mathrm{max}}={k}_{\mathrm{cat}}{E}_{\mathrm{tot}}$) for the enzymes catalysing the conversion of A_{0} and M_{0} into A_{1} and M_{1} (enzyme E_{0}), and the turnover of A_{1} into A_{0} (enzyme ${E}_{a}$), respectively, while the parameters ${K}_{M,{E}_{0}}$ and ${K}_{M,{E}_{a}}$ are the individual or combined MichaelisMenten coefficients for these enzymes’ substrates (i.e. for A_{0} and M_{0} and A_{1}, respectively). The term $\alpha $ is (in this case where all reactions are irreversible) equal to ${V}_{\mathrm{max},{E}_{a}}{k}_{\mathrm{in}}$, and in general is a positive expression comprising k_{in}, and the MichaelisMenten coefficients and the ${V}_{\mathrm{max}}$ parameters of the background enzymes in the model (see Appendix 3, Equations 7; 9; 11). The steady states for the model with all enzymatic conversions being reversible, and for a model with degradation and synthesis of A_{0} and A_{1}, are given in Appendix 3. The steady state solutions of these alternative models are structurally akin to (1), and do not alter the qualitative conclusions we make in what follows.
A key property of (1) is that it contains terms in the denominator that involve a subtraction. The presence of these terms introduces a limit on the parameter values for the system to attain a positive steady state. Specifically, we obtain the following conditions for positive steady states to exist:
Additionally, the ‘shape’ of (1) indicates a ‘threshold effect’ on the steady state value of m_{0}, where it would rise towards infinity as k_{in} increases towards the lower one among the limits given in (2) (see Figure 1B).
Why does (1) show this specific form, leading to these limits? We find that this is a direct consequence of the steady state condition, where metabolite production and consumption rates need to be the same at steady state. In the case of cosubstrate cycling, the production rate of M_{0} is given by k_{in}, while its consumption rate is a function of the ${V}_{\mathrm{max},{E}_{0}}$ and the concentration of A_{0}. In turn, the concentration of A_{0} is determined by its regeneration rate (which is a function of ${K}_{M,{E}_{a}}$ and ${V}_{\mathrm{max},{E}_{a}}$) and the pool size (${A}_{\mathrm{tot}}$). This explains the inequalities given in (2) and shows that a cycled cosubstrate creates the same type of limitation (mathematically speaking) on the flux of a reaction it is involved in, as that imposed by the enzyme catalysing that reaction (E_{0} in this example) (see Figure 1C & D). We also show that considering the system shown in Figure 1A as an enzymatic reaction without cosubstrate cycling leads to only the constraint $k}_{\mathrm{i}\mathrm{n}}<{V}_{\mathrm{m}\mathrm{a}\mathrm{x},{E}_{0}$, while considering it as a nonenzymatic reaction with cosubstrate cycling only, leads to only the constraint ${k}_{\mathrm{i}\mathrm{n}}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}{V}_{\mathrm{m}\mathrm{a}\mathrm{x},{E}_{a}}/({K}_{M,{E}_{a}}+{A}_{\mathrm{t}\mathrm{o}\mathrm{t}})$ becoming the sole limitation on the system (see Appendix 3). In other words, the two limitations act independently.
To conclude this section, we reiterate its main result. The flux of a reaction involving cosubstrate cycling is limited either by the kinetics of the primary enzyme mediating that reaction, or by the turnover rate of the cosubstrate. The latter is determined by the cosubstrate pool size and the kinetics of the enzyme(s) mediating its turnover.
Cosubstrate cycling causes a flux limit on linear metabolic pathways
We next considered a generalised, linear pathway model with $n+1$ metabolites and arbitrary locations of reactions for cosubstrate cycling, for example as seen in upper glycolysis (Appendix 1—figure 1A). In this model, we only consider intrapathway metabolite cycling, i.e. the cosubstrate is consumed and regenerated solely by the reactions of the pathway. Here, we show results for this model with 5 metabolites as an illustration (Figure 2A), while the general case is presented in Appendix 4.
We find the same kind of threshold dynamics as in the single reaction case. When k_{in} is above a threshold value, the metabolite M_{0} accumulates towards infinity and the system does not have a steady state (Figure 2B). A numerical analysis, as well as our analytical solution, reveals that the accumulation of metabolites applies to all metabolites upstream of the first reaction with cosubstrate cycling (Figure 2C and Appendix 4). Additionally, metabolites downstream of the cycling reaction accumulate to a steady state level that does not depend on k_{in} (Figure 2C and Appendix 3—figure 1). In other words, pathway output cannot be increased further by increasing k_{in} beyond the threshold. Finally, as k_{in} increases, the cycled metabolite pool shifts towards one form and the ratio of the two forms approaches zero (Figure 2C).
An analytical expression for the threshold for k_{in}, like shown in (2), could not be derived for linear pathways with $n>3$, but our analytical study indicates that (i) the threshold is always linked to ${A}_{\mathrm{tot}}$ and enzyme kinetic parameters, and (ii) the concentration of all metabolites upstream (downstream) to the reaction coupled to metabolite cycling will accumulate towards infinity (a fixed value) as k_{in} approaches the threshold (see Appendix 4). In Figure 2, we illustrate these dynamics with simulations for a system with $n=4$.
We also considered several variants of this generalised linear pathway model, corresponding to biologically relevant cases as shown in Appendix 1—figure 1. These included (i) intrapathway cycling of two different metabolites, as seen with ATP and NADH in combined upper glycolysis and fermentation pathways (Appendix 5), (ii) different stoichiometries for consumption and regeneration reactions of the cycled metabolite, as seen in upper glycolysis (Appendix 6), and (iii) cycling of one metabolite interlinked with that of another, as seen in nitrogen assimilation (Appendix 7). The results in the Appendices confirm that all these cases display similar threshold dynamics, where the threshold point is a function of the cosubstrate pool size and the enzyme kinetics.
Cycled metabolite related limit could be relevant for specific reactions from central metabolism
Based on flux values that are either experimentally measured or predicted by flux balance analysis (FBA), many reactions from the central carbon metabolism of the model organism Escherichia coli are shown to have lower flux than expected from the kinetics of their immediate enzymes (i.e. ${V}_{\mathrm{max}}$) (Davidi et al., 2016). This finding is based on calculating ${V}_{\mathrm{max}}$ from in vitro measured k_{cat} values of specific enzymes and their in vivo levels based on proteomics studies in E. coli (see Materials and methods). The flux and enzyme concentration data were from other studies which measured them during the exponential phase in E. coli growing on minimal media supplemented with various carbon sources (Schmidt et al., 2016; Gerosa et al., 2015). If we consider measured fluxes for each reaction as a proxy for k_{in} (notice that these two would be equal at steady state), we can conclude from the fact that there were no observed substrate accumulation in these reactions, as an indication for the analysed reactions carrying fluxes below the first limit identified above in (2). There could be several explanations for this observation of measured fluxes being lower than the limit set by measured enzyme kinetics and level. One simple explanation could be that there is a discrepancy between in vitro measured enzyme kinetics and in vivo realised ones. Alternatively, this discrepancy can be low, but the lower flux could be arising because there are additional limiting factors other than the enzymes mediating the main reaction. Among such additional limiting factors, substrate limitation and thermodynamic effects are shown to partially explain observed lower fluxes in some reactions (Davidi et al., 2016; see also below results). Here, we highlight that the presented theory shows that an additional possible limitation could be the cosubstrate pool size and turnover dynamics.
To explore this possibility, we reanalysed the flux values compiled previously (Davidi et al., 2016; Gerosa et al., 2015) and focused solely on reactions that are linked to ATP, NADH, or NADPH pools (see Materials and methods and Supplementary file 1). The resulting dataset contained fluxes, substrate concentrations, and enzyme levels for 45 different reactions determined under 7 different conditions along with turnover numbers and kinetic constants of the corresponding enzymes. In total, we gathered 49 combinations of enzymefluxk_{cat} values with full experimental data and 259 combinations with only FBApredicted flux values. We compared the flux values that would be expected from the primary enzyme limit identified above, under all conditions analysed (Figure 3A), and in addition checked whether the saturation effect of the primary substrate could explain the difference (Appendix 8—figure 1). We found that in both cases, about 80% of these reactions carry flux lower than what is expected from enzyme kinetics (Appendix 8—figure 2), suggesting that the limits imposed by cofactor dynamics might be constraining the flux further. The low number of the cases where the flux exceeds the limit might be due to uncertainties in measurement of flux, enzyme or substrate level.
It is also possible that observed lower fluxes are due to thermodynamic limitations. This is very difficult to analyse without more data, as calculating reaction thermodynamics requires knowledge of concentrations for all substrates and products, as well as enzyme MichaelisMenten constants in both forward and backward directions. This information is currently not available except for few of the reactions among the ones we analysed. Nevertheless, to give as much insight as possible on the thermodynamic effect, we analysed the physiological Gibbs free energy (the ${\mathrm{\Delta}}_{r}{G}^{\prime}$ is calculated assuming that all reactants are at 1 mM and pH = 7) against the normalized flux – $v/({E}_{0}\cdot {k}_{cat})$ (Appendix 8—figure 3). This shows that although in few cases, such as malate dehydrogenase (MDH), the normalised flux seems to be greatly reduced by the thermodynamic barrier, the general picture is that there is little correlation between reaction flux and thermodynamics.
We have also checked the relation between fluxes and cosubstrate pool sizes. Cosubstrate pool sizes do change between different conditions, and we note that such changes cannot be due to flux changes in cosubstrate utilising reactions. But, on the other hand, changes in pool size can affect flux in those reactions, where cosubstrate dynamics is limiting (as predicted by the theory). For both measured and FBApredicted fluxes, we find that several reactions show significant correlation between flux and cosubstrate pool size (see Figure 3B–D, see also Appendix 8—table 1 and Appendix 8—figure 4). In the case of FBApredicted fluxes, however, we note that these results can be confounded due to additional, fluxtoflux correlations and correlations between pool sizes and growth rate. Among reactions with measured fluxes, the three reactions with high correlation to pool size are those mediated by malate dehydrogenase (MDH), linked with NADH pool, phosphoglycerate kinase (PGK), linked with the ATP pool, and glucose6phosphate dehydrogenase (G6PDH) linked with the NADPH pool.
In summary, these results show that for reactions involving cosubstrate cycling (1) measured fluxes are lower than those predicted by kinetics of the primary enzyme (i.e. enzyme involved in substrate conversion) alone, and (2) there is – for some reactions – a correlation between flux and cosubstrate pool size. Both observations could indicate cosubstrate pool sizes and/or cosubstrate cycling dynamics being a limiting factor for flux. We can not state this as a certainty, however, as there are possibly other factors acting as the extra limitation, including thermodynamic effects. These points call for further experimental analysis of cosubstrate cycling within the study of metabolic system dynamics.
Cosubstrate cycling allows regulation of branch point fluxes
In addition to its possible constraining effects on fluxes, we wondered if cosubstrate dynamics can offer a regulatory element in cellular metabolism. In particular, cosubstrate cycling can commonly interconnect two independent pathways, or pathways branching from the same upstream metabolite, where it could influence flux distributions among those pathways. To explore this idea, we considered a model of a branching pathway, with each branch involving a different cosubstrate, $A$ and $B$ (Figure 4A and Appendix 1). This scenario is seen in the synthesis of certain amino acids that start from a common precursor but utilise NADH or NADPH, for example Serine and Threonine.
We hypothesised that regulating the two cosubstrate pool sizes, ${A}_{\mathrm{tot}}$ and ${B}_{\mathrm{tot}}$, could allow regulation of the fluxes on the two branches. To test this hypothesis, we ran numerical simulations with different cosubstrate pool sizes and influx rates into the branch point. We found that the ratio of fluxes across the two branches can be regulated by changing the ratio of ${A}_{\mathrm{tot}}$ to ${B}_{\mathrm{tot}}$ (Figure 4B). The regulation effect is seen with a large range of k_{in} values, but the threshold effect is still present with high enough k_{in} values leading to loss of steady state and metabolite build up. In that case, the resulting metabolite buildup can affect either branch depending on which cosubstrate has the lower pool size (see upper corner regions on Figure 4B). There is also a regime of only the upstream, branch point metabolite buildingup, but this happens only when all reactions are considered as reversible and the extent of it depends on turnover rates of the two cosubstrates (Appendix 9—figure 1).
In the nobuildup, steady state regime, changing the pool size ratio of the two cosubstrates causes a change in fluxes and metabolite levels, The change in flux ratio is of the same order as the change in pool size ratio (Figure 4C & D), while the change in the ratio of metabolite levels is in general less. This relation between pool size ratio and flux ratio on each branch is unaffected by the value of k_{in}. We also evaluated the level of regulation that can be achieved by varying the turnover rates of $A$ and $B$. The flux regulation effect in this case is weaker, unless the difference in the turnover rates is large and the influx rate is close to the threshold (Appendix 9—figure 2).
Interpathway cosubstrate cycling limits maximum influx difference and allows for correlating pathway outfluxes despite influx noise
We next considered a simplified model of two independent pathways interconnected by a single cosubstrate pool (Figure 5A and Appendix 10). This model can represent several different processes in metabolism, for example the coupling between the TCA cycle and the respiratory electron transfer chain, through NADH generation and consumption respectively, or the coupling between the pentose phosphate pathway and some amino acid biosynthesis pathways (notably Methionine), through NADPH generation and consumption respectively. We hypothesised that such interpathway cosubstrate cycling might cause; (1) the cosubstrate related limit to relate to difference in pathway influxes, rather than input into one pathway, and (2) a coupling of the pathway output fluxes against influx fluctuations, such that the output fluxes remain correlated to each other, despite differences in influx levels.
To address the first hypothesis, we used analytical methods and explored the relation between the systems’ ability to reach steady state and system parameters. We found that, indeed, for this coupled system, the ability to reach steady state depends on the influx difference between two pathways, as well as the total pool size and the kinetic parameters relating to pathwayindependent turnover of the cosubstrate (see Appendix 10). In other words, for two pathways coupled via cosubstrate cycling, the cyclingdependent flux limit for each is not determined by their own influx but rather on how different this is to the coupled pathways’ influx (Appendix 10—figure 1).
To test the second hypothesis about the output coupling, we considered the correlation of the steadystate outputs of the pathways with random fluctuations in their influx (Figure 5B). We found that there is either a high level of anticorrelation or correlation between pathway outputs for all pool sizes tested (blue and yellow regions in Figure 5B). As the pool size decreases, the system reaches a point where there is a transition from anticorrelation to high correlation in product outputs (blue to yellow region in Figure 5B). At low pool sizes, pathway outputs are fully correlated despite significant fluctuation in pathway influx (Figure 5C & D). Within this correlated regime, we identified two different subregimes. The first is a regime where the metabolite concentrations stay relatively constant despite the influx noise (Figure 5C). This regime arises because the influx fluctuations are occurring at a much faster rate than the pathway kinetics and the system is rather nonresponsive to influx noise. In a second regime, the influx noise is at time scales comparable to pathway kinetics. Here, the metabolite concentrations can readily change with the influx changes, and the system is ‘responsive’, yet the output levels are still correlated (Figure 5D). This regime is directly a result of cosubstrate cycling dynamics. Because the turnover of cosubstrate is essentially coupling the two pathways, their outputs become directly correlated. This effect does not depend on whether pathway reactions are modelled as reversible or irreversible, and the results here for ${V}_{\mathrm{max},{E}_{a}}=0.01$ are representative of those for ${V}_{\mathrm{max},{E}_{a}}=0$ (see Appendix 10—figure 2). As we increase the rate of the assumed background reaction, that is pathwayindependent turnover of the cosubstrate, we find that these results remain qualitatively the same, but the transition point from anti correlation in outputs to correlation, shifts to lower ${A}_{\mathrm{tot}}$ values (Appendix 10—figure 2).
These results show that coupling by cosubstrate cycling can introduce a limit on influxes of independent pathways or metabolic processes. Furthermore, such coupling can allow high correlation in the pathway outputs, despite significant noise in the inputs of those pathways. These effects are most readily seen where the turnover of the coupling cosubstrate by other processes is low. We note that an example case for such a scenario is the coupling of respiration and oxidative phosphorylation, where transmembrane proton cycling takes the role of the cycled cosubstrate (Stucki, 1980).
Discussion
We presented a mathematical analysis of metabolic systems featuring cosubstrate cycling and showed that such cycling introduces a threshold effect on system dynamics. As the pathway’s influx rate, k_{in}, approaches a threshold value, the steady state concentrations of metabolites that are upstream of a reaction linked to cosubstrate cycling, increase towards infinity and the system cannot reach steady state. Specifically, for reactions involving cosubstrates, there are two thresholds on influx rate, one relating to the kinetics of the enzyme directly mediating that reaction, and another relating to the kinetics of the enzymes mediating the turnover of the cosubstrate and the pool size of that cosubstrate.
This second, additional constraint arising from cosubstrate cycling can be directly relevant for cell physiology. We particularly note that this threshold can be highly dynamic, and condition and celldependent. When cells have a permanently or occasionally lowered total cosubstrate pool size (i.e. lower ${A}_{\mathrm{tot}}$), or when they are placed under challenging conditions (e.g. high carbon or nitrogensource concentrations) causing higher k_{in} values across various pathways, their metabolic systems can be close to the threshold presented here. While both k_{in} and ${A}_{\mathrm{tot}}$ can be adjusted in the long term, for example by reducing substrate transporter expression or increasing cosubstrate biosynthesis, there can be short term impact on cells experiencing significant flux limitations and metabolite accumulations.
Comparing measured flux data against estimated flux values based on measured enzyme levels from proteomics and enzyme kinetics from in vitro studies, we have provided support that fluxes in cosubstrate linked reactions could indeed by limited by cosubstrate pool dynamics under physiological conditions. This analysis was based on the model organism E. coli and is limited even for this organism due to limited flux and proteomics data. For example, the data compiled here contained 14 cosubstrate reactions with experimentally measured fluxes, but only half of these could be used due to lack of measurement on enzyme concentrations. We hope that the presented theory will provide motivation to further expand the available data sets, especially for reactions relating to cosubstrate linked reactions. In this quest, we expect that the expansion of measurements to eukaryotic cells to be particularly challenging due to organellespecific pools, but some progress is being made to achieve at least mitochondrial and cytosolic measurements (Chen et al., 2016). Despite the current limitations, our databased analyses highlighted three key reactions, that are possibly limited by cosubstrate dynamics, and that are mediated by phosphoglycerate kinase (PGK), malate dehydrogenase (MDH), and glucose6phosphate dehydrogenase (G6PDH) and linked to ATP, NADH, and NADPH pools. Possible flux limitation of these reactions by cosubstrate dynamics can also be subjected to further experimental study – as we discuss further below.
Overall, the presented theoretical results could contribute to our understanding of two commonly observed metabolic dynamics that arise under increasing or high substrate concentrations, and that are shown to cause either ‘substrateinduced death’ (van Heerden et al., 2014) or ‘overflow metabolism’. The latter usually refers to a respirationtofermentation switch under respiratory conditions (e.g. the Warburg and Crabtree effects [Warburg, 1956; DiazRuiz et al., 2009; Basan et al., 2015; Meyer et al., 1984]), but other types of overflow metabolism, involving excretion of amino acids and vitamins, has also been observed (Ponomarova et al., 2017; Jiang et al., 2018). Several arguments have been put forward to explain these observations, including osmotic effects arising from high substrate concentrations causing cell death and limitations in respiratory pathways or cell’s protein resources causing a respirationtofermentation switch (DiazRuiz et al., 2009; Majewski and Domach, 1990; Basan et al., 2015). Notwithstanding the possible roles of these processes, the presented theory leads to the hypothesis that both substrateinduced death and metabolite excretions could relate to increasing substrate influx rate reaching close to the limits imposed by cosubstrate dynamics. There is experimental support for this hypothesis in the case of both observations. Substrateinduced death and associated mutant phenotypes are linked to the dynamics associated with ATP regeneration in glycolysis (Teusink et al., 1998; Koebmann et al., 2002; van Heerden et al., 2014). Based on that finding, it has been argued that cells aim to avoid the threshold dynamics through allosteric regulation of those steps of the glycolysis that involve ATP consumption (Teusink et al., 1998). In the case of respirationtofermentation switch, it has been shown that the glucose influx threshold, at which fermentative overflow starts, changes upon introducing additional NADH conversion reactions in both yeast and E. coli populations (Vemuri et al., 2006; Vemuri et al., 2007). In another supportive case, sulfurcompound excretions are linked to alterations in the NADPH pool through changes in the amino acid metabolism (OlinSandoval et al., 2019; Green et al., 2020).
Dynamical thresholds relating to cosubstrate pools would be relevant for all cosubstrates, and not just for ATP or NADH, which have been the focus of most experimental studies to date. We would expect that altering kinetics of enzymes involved in cosubstrate cycling can have direct impact on cell physiology, and in particular on metabolic excretions. This prediction can be tested by exploring the effect of mutations on enzymes linked to cosubstrate consumption and production, or by altering cosubstrate pool sizes and assessing effects of such perturbations on the dynamics of metabolic excretions. These tests can be experimentally implemented by introducing additional enzymes specialising in cosubstrate consumption or production (e.g. ATPases, oxidases, or other) and controlling their expression. It would also be possible to monitor cosubstrate pool sizes in cells in real time by using fluorescent sensors on key metabolites such as ATP or glutamate, or by measuring autofluorescence of certain pool metabolites, such as NAD(P)H, under alterations to influx rate of glucose or ammonium.
Besides acting as a flux constraint, we find that cosubstrate pools can also allow for regulation of pathway fluxes through regulation of pool size or turnover dynamics. We find that such regulation can take the form of balancing interconnected pathways, thereby ensuring correlation between outputs of different metabolic processes, or regulating flux across branch points. Regulation of fluxes through cosubstrate pools can act to adjust metabolic fluxes at time scales shorter than possible via gene regulation, and possibly at similar time scales as with allosteric regulation – especially when considering pool size alterations through exchange among connected pools. Possibility of such a regulatory role has been indicated experimentally, where total ATP pool size is found to change when yeast cells are confronted with a sudden increase in glucose influx rate (Walther et al., 2010). In that study, the change in the ATP pool is found to link to the purine metabolism pathways, which are linked to several conserved moieties; GTP, ATP, NAD, NADP, Sadenosylmethionine, and Coenzyme A. These findings suggest that cells could dynamically alter pool sizes associated with different parts of metabolism, limiting flux through some pathways, while allowing higher flux in others, and thereby shifting the metabolites from the latter to the former. This could provide a dynamic selfregulation and the pool sizes of key cosubstrates could be seen as ‘tuning points’ controlling a more complex metabolic system. We thus propose further experimental analyses focusing on cosubstrate pool sizes and turnover dynamics to understand and manipulate cell physiology.
Materials and methods
Model of a single reaction with cosubstrate cycling
Request a detailed protocolThe metabolic system shown in Figure 1A comprises the following biochemical reactions:
where metabolites are denoted by ${M}_{i}$ and the different forms of the cosubstrate are denoted by ${A}_{i}$. We assume additional conversion between A_{1} and A_{0}, mediated through other enzymatic reactions. The parameters k_{in}, and k_{out} denote the in and out flux of M_{0} and M_{1} respectively, from and to other pathways or across cell boundary. The ordinary differential equations (ODEs) for the system shown in (3) (and Figure 1A), using irreversible MichaelisMenten enzyme kinetics would be:
where m_{i} and a_{i} denote the concentrations of ${M}_{i}$ and ${A}_{i}$ respectively, ${K}_{M}$ denotes a composite parameter of the MichaelisMenten coefficients of the enzyme for its substrates, and ${V}_{\mathrm{max}}$ is the total enzyme concentration times its catalytic rate (i.e. ${V}_{\mathrm{max}}={k}_{\mathrm{cat}}{E}_{\mathrm{tot}}$) (see Appendix 11—table 1 for a list of parameters and their units). We further have the conservation relation ${a}_{0}+{a}_{1}={A}_{\mathrm{tot}}$, where ${A}_{\mathrm{tot}}$ is a constant. This assumption would be justified when influx of any form of the cycled metabolite into the system is independent of the rest of the metabolic system (see further discussion and analysis in Appendix 2). The steady states of (4) can be found by setting the left side equal to zero and performing algebraic rearrangements to isolate each of the variables. The resulting analytical expressions for steady state metabolite concentration are shown in (2), and in Appendix 3 for this model with reversible enzyme kinetics, as well as for other models.
Symbolic and numerical computations
Request a detailed protocolFor all symbolic computations, utilised in finding steady state solutions and deriving mathematical conditions on rate parameters, we used the software Maple 2021, as well as theoretical results presented in Torres and Feliu, 2021. To run numerical simulations of select systems, we used Python packages with the standard solver functions. All numerical simulations were performed in the Python environment. The main model simulation files relating to Figures 4 and 5 are provided as , while all remaining simulation and analysis scripts are made available through a dedicated repository (West et al., 2023).
Reaction fluxes and enzyme kinetic parameters
Request a detailed protocolTo support the model findings on cosubstrate pools acting as a possible limitation on reaction fluxes, we analysed measured and FBAderived flux data collated previously (Davidi et al., 2016; Gerosa et al., 2015). We focused our analyses on reactions involving cosubstrates only. We compared measured (or FBAderived) fluxes to flux thresholds based on enzyme kinetics (i.e. first condition in Eq. 2). To calculate the latter, we used data on enzyme kinetics and levels as collated in Davidi et al., 2016, which is based on the BRENDA database (Chang et al., 2021) and proteomicsbased measurements (Schmidt et al., 2016). We note that most available kinetic constants for enzymes have been obtained under in vitro conditions, which can be very different from those of the cytosol (GarcíaContreras et al., 2012). When comparing flux levels against cosubstrate pool sizes, we used the matching, measured poolsizes from Gerosa et al., 2015. All the data used in this analysis is provided in the Supplementary file 1, and through a dedicated repository (West et al., 2023).
Appendix 1
The mathematical approach and modelling setting
As explained in the main text, we are interested in understanding the effect of cosubstrate cycling on the flux through metabolic pathways, such as those shown in Appendix 1—figure 1.
In these Appendices, we use a generic cosubstrate pair, denoted as A_{0}, A_{1}. We consider the synthesis or degradation of the cosubstrate pair, or consider it as a conserved moiety, i.e. having a fixed total concentration. Our generic cosubstrate pair, A_{0} and A_{1}, can be taken as representing a specific cosubstrate, such as NAD(H), but note that the mathematical analyses presented would be applicable to any cosubstrate pair in natural metabolic pathways (as discussed in the main text).
For our analyses, we consider a generalised model of a linear metabolic pathway, as well as additional metabolic pathway structures. Throughout the presented analyses, we consider reactions to be either enzyme mediated or not, and when they are enzyme mediated, we consider them either to be reversible or irreversible. In the former case, the enzymatic conversions are shown as ${M}_{i1}\rightleftharpoons {M}_{i}$, while in the latter case, they are shown as ${M}_{i1}\u27f7{M}_{i}$. These notations do not show enzyme complexes explicitly, but we use enzymatic rates derived from reaction schemes accounting for enzyme complexes (see below).
In certain models, we consider some, or all, cycling reactions of the cosubstrate to occur independently of the enzymatic reactions involved in the metabolic pathway, for example due to hydrolysis reactions. We refer to this type of recycling as free conversion, for example in the case of a generic cosubstrate considered here, we have:
We talk about irreversible cosubstrate conversion, if ${k}_{5}=0$ or ${k}_{6}=0$, that is, only conversion in one direction is considered. We talk about no free cosubstrate conversion, if ${k}_{5}={k}_{6}=0$, that is, the cosubstrate cycling is only related through the reactions in the metabolic pathway.
Enzyme kinetics
Each metabolic pathway is modelled using either MichaelisMenten (irreversible case) or Haldane (reversible case) enzyme kinetics, for the individual reactions it comprises. The general kinetics can be expressed as follows, where we let ${a}_{0},{a}_{1}$ denote the concentrations of the cosubstrate pair A_{0} and A_{1}, respectively and m_{i} to denote the concentration of ${M}_{i}$, the $i$th metabolite in the pathway.
In the case of a reversible, enzymatic reaction involving a cosubstrate and assuming simultaneous binding of both substrates to the enzyme, we have the following reaction scheme:
For this reversible reaction scheme, the rate of production of ${M}_{i}$ takes the form
Likewise, for the reversible enzymatic conversion ${M}_{i1}\rightleftharpoons {M}_{i}$, we have the following reaction scheme:
The rate of production of ${M}_{i}$ is given by
In both of these reversible rate equations, the parameters $K$ and $L$ are equivalent to the Haldane coefficients ${K}_{S}$ and ${K}_{P}$, respectively and are given by
When there are two substrates that take part in the reaction, the k_{on} and ${k}_{\mathrm{off},2}$ parameters are composite parameters composed of two binding coefficients, one for each substrate. This does not affect the derivations, so for convenience we use ${K}_{S}$ and ${K}_{P}$.
The parameters $E$ and $F$ correspond to the Haldane coefficients ${k}_{\mathrm{cat}+}$ and ${k}_{\mathrm{cat}}$, multiplied by the total enzyme concentration (denoted ${E}_{\mathrm{tot}}$, below), and are given by
For the irreversible enzymatic reaction, the reaction schemes simplify to:
And the rate of production for the two cases are given by
where ${E}_{i}$ is the catalytic rate coefficient of the $i$th enzyme multiplied by its total concentration, and ${K}_{i}$ is its MichaelisMenten coefficient. Again, when there are two substrates, the k_{on} parameter is a composite parameter composed of two binding coefficients, one for each substrate. As in the reversible case, this does not affect the derivations, so we use ${K}_{i}$ for convenience. Influx and outflux follow nonenzymatic reaction kinetics, with reaction rate constants as indicated by the labels of the reactions.
Analytical approach
Our mathematical analysis is primarily concerned with finding conditions on the kinetic parameters, if any, that imply that the system has a positive steady state. This is different than system reduction, for example as done in the analyses leading to MichaelisMenten kinetics. Our analysis distinctively solves the entire system for steady states and determines conditions on kinetic parameters to satisfy the steady state equations.
Thus, for each of the metabolic pathway motifs we consider, we build the ODEs defining the rates of change of variables, find the conservation laws among variables, and consider a system of equations whose solutions are the steady states of the ODEs constrained by the conservation laws. We then follow one of two strategies. We first attempt to solve all equations for all concentrations. For some systems, we readily get an expression in terms of the parameters of the system. For other systems, this approach is not possible. In this case, using all equations in the system but one, we solve for the steady states of all concentrations but one. This gives all concentrations in terms of the remaining concentration, say $x$. Plugging these expressions in the remaining equation of the system, we obtain a final equation whose solutions characterize the steady states of the system. We need then to study when the solutions obtained this way are positive.
We are also interested in proving if a given system has a positive steady state for all parameter combinations, and that this steady state is stable. When there is one positive steady state, we find the Hurwtiz determinants associated with the characteristic polynomial of the Jacobian of the system of ODEs, evaluated at the steady state. If these are all positive, then the steady state is asymptotically stable (Torres and Feliu, 2021).
To decide on the existence of a steady state, throughout the analysis, we will use repeatedly the following lemma, which is a consequence of the wellknown Descartes’ rule of signs.
Lemma 1. Let $p\mathit{}\mathrm{(}x\mathrm{)}$ be a univariate polynomial of degree two, with negative leading term. If at some value $T$, we have $p(T)>0$, then $p$ has a root in the interval $\mathrm{(}\mathrm{0}\mathrm{,}T\mathrm{)}$ if and only if the independent term of $p$ is negative.
Proof. The Descartes’ rule of signs establishes that the number of positive roots of a polynomial cannot exceed the number $\tau $ of sign changes in the sequence of coefficients ignoring zero coefficients, and the difference between $\tau $ and the number of positive roots is an even number. As the polynomial $p$ in the statement attains positive values, it must have some positive coefficient. Furthermore, as the degree two polynomial has negative leading term, the sequence of the sign of terms (when terms are ordered from lowest exponent to highest) is one of the following $++$, $+$, $+$, $+\mathrm{\hspace{0.17em}0}$, $0+$.
If the independent term is positive or zero, then the sign sequence is one of $++$, $+$, $+\mathrm{\hspace{0.17em}0}$, $0+$. In this case, there is one sign change in the sequence, and it follows that the polynomial has exactly one positive root. As $p(0)>0$, $p(T)>0$ and $p$ becomes negative as $x$ goes to $+\mathrm{\infty}$, the root must be in the interval $(T,+\mathrm{\infty})$.
If the independent term is negative, then the sign sequence is $+$. The polynomial is negative both at 0 and at $+\mathrm{\infty}$. As $p(T)>0$, there must be a positive root in $(0,T)$ and one in $(T,+\mathrm{\infty})$, and there cannot be more by the Descartes’ rule of signs. From this, the statement of the lemma follows.
Appendix 2
Considering cosubstrate pool size
We first consider cycling of a generic cosubstrate A_{0} and A_{1}, with biosynthesis and degradation of both forms.
We suppose that the biosynthesis occurs at a constant rate, while degradation and cycling are proportional to the concentration of the relevant chemical species. Writing ${a}_{0}=\left[{\mathrm{A}}_{0}\right]$ and ${a}_{1}=\left[{\mathrm{A}}_{1}\right]$, the differential equations for these concentrations are:
Since all the terms are linear or constant, the steady state values are the solutions of the linear equation:
The steady states are then found to be:
If we consider the case where the synthesis and degradation rates of the different forms of the cosubstrate (i.e. cycled metabolite) are the same, i.e. k_{1} = k_{3}=k_{s} and k_{2} = k_{4}=k_{d}, these equations simplify to:
and the eigenvalues of the Jacobian of the system evaluated at this steady state are always real and negative. When k_{d} is sufficiently small compared to cosubstrate conversion rates, it can be safely neglected in the brackets, resulting in the expression of steady state formulas as:
We can compare the above expressions with those obtained from the case, where we assume a constant pool size of the cycled metabolite (i.e. ${k}_{1}={k}_{2}={k}_{3}={k}_{4}=0$). In that case, the steady states are ${a}_{0}=T{k}_{6}/({k}_{5}+{k}_{6})$ and ${a}_{1}=T{k}_{5}/({k}_{5}+{k}_{6})$, where $T$ is the total pool size. Thus, under the limit of degradation rates being much smaller than conversion rates, the two cases will be identical and cosubstrates will act as a conserved moiety for the rest of the metabolic system.
If we now assume that the cycling of cosubstrates is an enzymatic reaction and make the same simplifying assumptions as above that ${k}_{1}={k}_{3}={k}_{s}$ and ${k}_{2}={k}_{4}={k}_{d}$, the ODEs for the system are:
The only real and positive steady state is now found to be:
This is stable as long as all parameters are positive. Note that in the case of ${K}_{a}={L}_{a}$, the steady state solutions converge to a real number less than infinity by l’Hopital’s Rule. Also, note that the sum ${a}_{0}+{a}_{1}$ is constant as in the nonenzymatic case presented above. Thus, whether the metabolite cycling is considered as a nonenzymatic or enzymatic (i.e. following MichaelisMenten kinetics) reaction, the cosubstrates will act as a conserved moiety for the rest of the metabolic system in both cases.
Appendix 3
Single reaction models
In this section, we derive results for a single reaction of two metabolites, involving cosubstrate cycling or not, as presented in the main text.
Enzymatic reaction with cosubstrate cycling
We first reconsider the case where all reactions, including the offpathway cycling, are enzymatic (hence are modelled with MichaelisMenten kinetics). The reactions are:
This corresponds to the motif depicted in Figure 1A of the main text, and the resulting ODEs are:
This ODE system has one conservation law, namely the sum of a_{0} and a_{1} is constant:
The steady states of the system are:
where, we introduced the composite parameter $\alpha $, as follows:
For the steady state equations given above to be positive, the following conditions must be satisfied:
Note that as $\frac{{F}_{a}{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}}{{L}_{a}+{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}}<{F}_{a}$, the second condition readily implies $k}_{\text{in}}<{F}_{a$ and $\alpha $ is positive. When there is a positive steady state, then it is asymptotically stable.
If the main pathway is irreversible, but the cosubstrate reaction is still reversible, the ODEs describing the system dynamics are:
The steady state of the system is:
where the composite parameter $\alpha $ is defined (differently to the reversible case) as:
For this to be positive the same conditions as in the reversible case must be satisfied:
When these are satisfied, $\alpha $ is positive, and the positive steady state is asymptotically stable.
If the main pathway is irreversible, and the cosubstrate reaction only flows from A_{1} to A_{0}, the ODEs describing the system dynamics are:
The steady state of the system is:
where the composite parameter $\alpha $ is defined as:
For these steady states to be positive the same conditions as in the other two cases must be satisfied:
When these are satisfied, the positive steady state is asymptotically stable.
Enzymatic reaction with cosubstrate cycling, biosynthesis and degradation
To extend the previous analysis we consider the same simple scenario where a pathway involves a single cosubstrate consuming reaction and a back reaction regenerating the cosubstrate, with the addition of synthesis and degradation of the cosubstrate forms. This system is comprised of the following reactions:
The resulting system of ODEs is:
The steady state concentrations are then given by:
These expressions are positive if and only if
If the main path is irreversible, similarly to the previous section, the ODEs describing the system’s dynamics, are:
The steady state in this case is:
These expressions are positive if and only if
Comparing these results with those of Subsection C.1, we see some similar themes. Firstly, the total amount of cosubstrate ${a}_{0}+{a}_{1}$ at steady state is a constant, even when it is not explicitly considered to be a conserved moiety. Secondly, one of the conditions for a positive steady state is $k}_{\text{in}}<{E}_{1$, and any other conditions take the form ${k}_{\text{in}}<f$, where $f$ is a function of the parameters controlling synthesis, degradation and the turnover of the cosubstrate. Thus, in the analysis that follows, we only consider the case of conserved cosubstrate, as adding synthesis and degradation only affects the quantitative results, not the qualitative behaviour.
Enzymatic reaction without cosubstrate cycling
Considering an enzymatic reaction without cosubstrates, the reactions are:
The resulting system of ODEs is:
In this case, there is no conservation law, and the steady states of the system are:
These expressions are positive if and only if $k}_{\mathrm{i}\mathrm{n}}<{E}_{1$ and the steady state is asymptotically stable when this holds.
If the main path is irreversible, the ODEs describing system dynamics are:
and the steady state is:
Again, expressions are positive if and only if $k}_{\mathrm{i}\mathrm{n}}<{E}_{1$ and the steady state is asymptotically stable when this holds.
From this we see that the condition $k}_{\mathrm{i}\mathrm{n}}<{E}_{1$ for stability in the enzymatic system with cosubstrate cycling is a result of the reaction being enzymatic.
Nonenzymatic reaction with cosubstrate cycling
Consider a theoretical case of nonenzymatic reactions for the main reaction:
The resulting system of ODEs, describing system dynamics, is:
There is the conservation law
The steady state of the system is:
These expressions are positive if and only if $k}_{\mathrm{i}\mathrm{n}}<{k}_{4}{A}_{\mathrm{t}\mathrm{o}\mathrm{t}$, and when this is satisfied the steady state is asymptotically stable.
If the system reactions are irreversible, the ODEs describing the system dynamics are:
and the steady state is:
As in the reversible case, these expressions are positive when $k}_{\mathrm{i}\mathrm{n}}<{k}_{4}{A}_{\mathrm{t}\mathrm{o}\mathrm{t}$ and the steady state is asymptotically stable when this holds.
Hence, we see that introducing cosubstrate cycling always introduces a new constraint into the system, even in this simple case where the cycling is not enzymatic.
Appendix 4
Linear, arbitrary length, pathway model with cosubstrate cycling
We consider next a linear, generic pathway of $n+1$ metabolites, and comprising the following reaction mechanism:
Steady states of the linear pathway model: the case with n=3
We first consider the dynamics of model (13) with $n=3$, as this is the minimal pathway length where the system makes biochemical sense. This system has the form:
Considering reversible reaction kinetics, the evolution of the concentrations of the species in time is modelled by the following ODE system:
The system has two conservation laws:
Solving recursively the steady state equations given by $\frac{d{m}_{0}}{dt}+\mathrm{\cdots}+\frac{d{m}_{3}}{dt}=0$, $\frac{d{m}_{3}}{dt}=0$, $\frac{d{m}_{2}}{dt}+\frac{d{m}_{3}}{dt}=0$, $\frac{d{m}_{1}}{dt}+\frac{d{m}_{2}}{dt}+\frac{d{m}_{3}}{dt}=0$ for ${m}_{3},{m}_{2},{m}_{1},{m}_{0}$ and the first conservation law, we obtain
By substituting recursively m_{2} in m_{1} and m_{1} in m_{0}, we obtain expressions of all variables at steady state in terms of a_{0}. For these to be positive, it needs to hold that $0<{a}_{0}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}$, and further
We substitute (16) into the remaining conservation law, and obtain a polynomial in a_{0} whose roots in the interval $(0,T)$ are in onetoone correspondence with the positive steady states, provided (17) holds. The polynomial is
where
As ${a}_{0}={A}_{\mathrm{tot}}$, we have $p({A}_{\mathrm{t}\mathrm{o}\mathrm{t}})>0$, and hence, by Lemma 1, the system has positive steady states if and only if $\alpha <0$.
Note that at ${k}_{\text{in}}=0$, $\alpha =\phantom{\rule{thinmathspace}{0ex}}{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}\phantom{\rule{thinmathspace}{0ex}}W({E}_{3}{k}_{\text{in}})({E}_{2}{k}_{\text{in}}){L}_{1}<0$. Hence, for ${k}_{\text{in}}$ small enough, $\alpha <0$ and also (17) holds, and the system has a positive steady state. The steady state value of a_{0} is found by solving the quadratic equation $p({a}_{0})=0$ and considering the smallest root.
We note that $\alpha $ is a polynomial in ${k}_{\text{in}}$ of degree 2 with negative independent term. The sign of the leading term depends on the choice of parameters. If the minimum in (17) is attained at ${k}_{\text{in}}={E}_{2}$ or ${k}_{\text{in}}={E}_{3}$, then at this value of ${k}_{\text{in}}$, $\alpha $ is positive independently of the rest of the parameters. Hence, in the region of interest $\alpha $ is negative only in an interval of the form $(0,B)$, where $B$ is the smallest positive root of $\alpha $. The bound for positive steady states is now given as
If the minimum is attained at $B$, then, as ${k}_{\text{in}}$ approaches $B$, the smallest positive root of $p({a}_{0})$ approaches 0 as $\alpha $ approaches zero. Using this in (16), we obtain that
If ${E}_{1}<B$, then ${k}_{\text{in}}$ needs to be smaller than E_{1}. When ${k}_{\text{in}}$ approaches E_{1}, a_{0} approaches the root of $p$ when ${k}_{\text{in}}={E}_{1}$, which is not zero. Then m_{0} still approaches $\mathrm{\infty}$ as it has ${E}_{1}{k}_{\text{in}}$ in the denominator, while ${m}_{1},{m}_{2}$ converge to some positive value.
Steady states of the linear pathway model: the general case
We consider the dynamics of model (13) with $n$ taking any positive value strictly larger than 1, and positioning the cosubstrate recycling in arbitrary places.
For the system to make sense biochemically, we require $0\le \ell <\rho 1\le n1$.
To write the equations in a simple format, we write
and
We start with the reversible reaction kinetics. Then the associated ODE system becomes:
This ODE system has precisely two conservation laws:
Note that we have the following equalities:
For $i=1,\mathrm{\dots},n2$, by solving (21) for x_{i}, we obtain the following recursive formulas:
Finally, from (22) and the conservation law ${a}_{0}+{a}_{1}={A}_{\mathrm{tot}}$, we obtain
These expressions are all positive if and only if $0<{a}_{0}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ and $k}_{\text{in}}<{E}_{i$ for $i=1,\mathrm{\dots},n$. Note that the value of m_{i} can now be found recursively from m_{n} using (23), as we show next.
Recall that $0\le \ell <\rho 1\le n1$. Then, it holds
We write for shortness
where for $i=0,\mathrm{\dots},n1$,
Then we claim that
We prove this by descending induction on $i$. Note that products over an empty index equal 1. For instance, $b}_{j}\cdots {b}_{j1$. For $i=n1$, (26) agrees with (23). Assume the formula is true for some , and we prove it for $i1$. To do so, we use (23) and the induction hypothesis for x_{i}. For $n1>i1\ge 0$, we have
This is (26) for $i1$. Hence, the formula holds. Finally, we can use that ${x}_{n}={m}_{n}=\frac{{k}_{\text{in}}}{{k}_{\text{out}}}$ to obtain:
Let
This gives:
Remember that ${b}_{\mathrm{\ell}},{b}_{\mathrm{\ell}+1},{b}_{\rho},{b}_{\rho +1}$ depend on ${a}_{0},{a}_{1}$, while the rest of $b$’s do not. For a product of the form ${b}_{u}\mathrm{\cdots}{b}_{j}$ with $u\le j$, we have the following:
If $u\le \mathrm{\ell}1$:
If $u=\mathrm{\ell}$:
If $\ell <u\le \rho 1$:
If $u=\rho $:
If $u>\rho$:
Summarising $b}_{u}\cdots {b}_{j$, equals
Observe that for all $i$, the summand of ${\sum}_{j=i}^{n1}{z}_{j}({b}_{i}\mathrm{\cdots}{b}_{j1})({c}_{j+1}\mathrm{\cdots}{c}_{n1})$ corresponding to $i=j$ is
which does not depend either on a_{0} or a_{1} .
In particular, we deduce that m_{i} for $i\le \mathrm{\ell}1$, the term $({b}_{i}\mathrm{\cdots}{b}_{n1})$ does not depend on ${a}_{0},{a}_{1}$, and in the sum ${\sum}_{j=i}^{n1}{z}_{j}({b}_{i}\mathrm{\cdots}{b}_{j1})({c}_{j+1}\mathrm{\cdots}{c}_{n1})$ there are summands involving $\frac{1}{{a}_{0}}$, for eample when $j=\mathrm{\ell}$. Hence m_{i}, for $i\le \mathrm{\ell}1$, goes to infinity as a_{0} goes to zero. Note that a_{1} goes to $T$ when a_{0} goes to zero. When $i=\mathrm{\ell}$, the denominator of m_{i} is multiplied by a_{0}. As the numerator has at least one term that is not multiple of a_{0}, ${m}_{\mathrm{\ell}}$ goes to infinity as well. We conclude that
When $i\ge \mathrm{\ell}$, then no summand in the numerator of m_{i} in (27–29) that involves $\frac{1}{{a}_{0}}$, and neither the denominator has a_{0} as factor. As the numerator has at least one term that is not multiple of a_{0}, m_{i} goes to a finite value as a_{0} goes to zero.
The number can be found using (Equations 28–30) and is a function of the parameters of the system, not involving $W$.
We consider now the remaining equation, namely the conservation law ${m}_{\mathrm{\ell}+1}+\mathrm{\cdots}+{m}_{\rho}+{a}_{0}=W$, to determine the value of a_{0}. We have
We have $({b}_{i}\mathrm{\cdots}{b}_{n1})=\frac{{a}_{0}}{{a}_{1}}{\overline{b}}_{i,n1}$ when $\ell <i<\rho$ and $({b}_{\rho}\mathrm{\cdots}{b}_{n1})={a}_{0}{\overline{b}}_{\rho ,n1}$. Also
Finally, for $\mathrm{\ell}+1\le i\le \rho 1$,
This gives
Hence $(*)$ vanishes if and only if
As ${a}_{1}={A}_{\mathrm{tot}}{a}_{0}$, this is a degree 2 polynomial in a_{0}. The leading term comes from the term $(0,0,1{a}_{0}W){k}_{\text{out}}\phantom{\rule{thinmathspace}{0ex}}{\overline{c}}_{\ell +1,n1}\phantom{\rule{thinmathspace}{0ex}}0,0,1{a}_{1}$, and is
which is negative. The independent term is
We divide by ${k}_{\text{out}}$ and define
where recall from (25) that
When ${a}_{0}={A}_{\mathrm{tot}}$, all terms multiplying a_{1} vanish, and then the polynomial is a sum of positive terms, hence positive. By Lemma 1, the system has a positive steady state if and only if
Note that when ${k}_{\text{in}}=0$ this inequality holds, as all terms with z_{j} vanish. When ${k}_{\text{in}}$ approaches one of ${E}_{i}$ with $\ell +1<i\le n$, the negative term of $\alpha $ approaches zero, and the inequality does not hold.
For example, for $n=3$, $\mathrm{\ell}=0$ and $\rho =2$, $\alpha $ was found in (18). For $n=4$, $\mathrm{\ell}=0$ and $\rho =2$, we have
Let $B$ be the smallest positive root of $\alpha =0$ as a polynomial in ${k}_{\text{in}}$, if it exists, or take $B=\mathrm{\infty}$ if not. Similarly, let ${B}^{\prime}$ be the second such root, if it exists, or ${B}^{\prime}=\mathrm{\infty}$ if not. If the smallest of ${E}_{j}$ is attained for some ${E}_{i}$ with $\ell +1<i\le n$, then $\alpha $ is positive at ${k}_{\text{in}}={E}_{i}$. In that case, as $\alpha $ is negative at ${k}_{\text{in}}=0$, there is at least one value of $k}_{\text{in}}<{E}_{i$ for all $i$ such that $\alpha =0$ and hence $B$ is finite. This shows that $\mathrm{min}({E}_{1},\mathrm{\dots},{E}_{n},B)=\mathrm{min}({E}_{1},\mathrm{\dots},{E}_{\mathrm{\ell}+1},B)$.
Putting it all together, we have shown that for all
the system has positive steady states, and if
Taking condition (32), if the minimum is attained at $B$, then when ${k}_{\text{in}}$ approaches $B$, the first positive root of the polynomial in (3031, ) approaches zero (as $\alpha $ goes to zero). As this determines the steady state value of a_{0}, we see that a_{0} approaches zero, and the m_{i} converge to the values given above. Specifically,
By the comment above, (32) cannot be attained at ${E}_{i}$ with $\ell +1<i\le n$. If (32) is attained at ${E}_{i}$ with $i\le \mathrm{\ell}+1$, then as ${k}_{\text{in}}$ approaches this minimum, a_{0} converges to a number. In this case, the concentrations that tend to infinity are those with ${E}_{i}{k}_{\text{in}}$ in the denominator:
Appendix 5
Multiple cosubstrate cycling along a single pathway – mimicking the case seen in glycolysis, combined with fermentation
In this section, we consider a scenario of intrapathway cycling with two different cosubstrates. This is a simplified version of the case seen in the combined pathways of glycolysis and fermentation. The reaction scheme we consider comprises:
We write ${A}_{0}=$ ATP, ${A}_{1}=$ ADP, ${A}_{2}=$ NAD, ${A}_{3}=$ NADH, for simplicity. The ODE system governing the dynamics of the network is:
There are four conservation laws:
We consider the equations $\frac{d{m}_{0}}{dt}+\mathrm{\cdots}+\frac{d{m}_{4}}{dt}=0$, $\frac{d{m}_{4}}{dt}=0$, $\frac{d{m}_{3}}{dt}+\frac{d{m}_{4}}{dt}=0$, $\frac{d{m}_{2}}{dt}+\frac{d{m}_{3}}{dt}+\frac{d{m}_{4}}{dt}=0$, $\frac{d{m}_{1}}{dt}+\frac{d{m}_{2}}{dt}+\frac{d{m}_{3}}{dt}+\frac{d{m}_{4}}{dt}=0$, and solve them iteratively for ${m}_{4},{m}_{3},{m}_{2},{m}_{1},{m}_{0}$ and obtain:
and ${a}_{1}={A}_{\mathrm{tot}}{a}_{0}$, ${a}_{3}={B}_{\mathrm{tot}}{a}_{2}$. As usual, a necessary condition for positive steady states is
For this system, we are left with two conservation laws that we have not used, and two steady state values that are still free. By plugging the expressions above into ${m}_{1}+{m}_{2}+{a}_{0}=W$, and ${m}_{2}+{m}_{3}+{a}_{2}=M$, we obtain a system of two polynomial equations in two variables. For some combinations of parameter values, there will be positive steady states, and for others, not, as the following examples shows.
By choosing
the polynomial system becomes
and the solutions are:
which are in the desired interval.
If instead we replace L_{3} by 1, we obtain a system whose solutions are
and hence there are no positive steady states.
This shows that there is a condition for positive steady states to exist, and although (34) is necessary, it is not sufficient.
Appendix 6
Different stoichiometries for cosubstrate cycling along a single pathway – mimicking the case seen in upper glycolysis
In this section, we consider a scenario of intrapathway cycling with varying stoichiometry in the pathway. In particular, we consider a simplified version of the case seen in upper glycolysis. The reaction scheme we consider comprises:
We let ${A}_{0}=\mathrm{ATP}$ and ${A}_{1}=\mathrm{ADP}$ as usual.
With MichaelisMenten kinetics, the ODE system is
The system has two conservation laws:
By considering the equation $0=2\frac{d{m}_{0}}{dt}+2\frac{d{m}_{1}}{dt}+\frac{d{m}_{2}}{dt}+\frac{d{m}_{3}}{dt}={k}_{\text{out}}{m}_{3}2{k}_{\text{in}},$ we obtain
Upon substitution of this value of m_{3} into $\frac{d{m}_{3}}{dt}=0$, $\frac{d{m}_{2}}{dt}+\frac{d{m}_{3}}{dt}=0$ and $2\frac{d{m}_{1}}{dt}+\frac{d{m}_{2}}{dt}+\frac{d{m}_{3}}{dt}=0$, and solving recursively for ${M}_{2},{m}_{1},{m}_{0}$, we obtain the following steady state relations:
and recall ${a}_{1}={A}_{\mathrm{tot}}{a}_{0}$. We see that if $0<{a}_{0}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}$, these expressions are all positive if and only if
As usual, we consider the remaining conservation law, ${m}_{1}+{m}_{2}+{a}_{0}=W$, together with these expressions, to find a polynomial $p({a}_{0})$ whose roots in $(0,{A}_{\mathrm{tot}})$ are in onetoone correspondence with positive steady states, provided (36) holds. The polynomial has degree 3, positive leading term, negative independent term and term of degree 2, and is positive at ${a}_{0}={A}_{\mathrm{tot}}$. With this information, we cannot immediately conclude that there is a positive root in $(0,{A}_{\mathrm{tot}})$. But for positive steady states to exist, we need the term of degree 1 to be positive (this follows from Descartes’ rule of signs, as usual), and this sets an extra condition on the parameters. When this happens, there will be two positive steady states.
Specifically, the term of degree 1 is
which depends on ${A}_{\mathrm{tot}}$ and $W$. To summarize for positive steady states to exist, we need
If ${k}_{\text{in}}$ is small enough, then the conditions hold as $\alpha >0$ at ${k}_{\text{in}}=0$.
However, these conditions are not sufficient. To see this, by inspecting the term $\alpha $, one can see that if ${K}_{2},{K}_{3}$ are small enough, then there will be two positive steady states, while if they are larger, there will be none. We have verified that both scenarios occur. For example, fix the parameter values to be
and note that (35) holds. With ${K}_{2}={K}_{3}=0.1$, the polynomial of interest becomes $p({a}_{0})=4{a}_{0}^{3}14.300{a}_{0}^{2}+7.360{a}_{0}0.448$, and so $\alpha >0$. The polynomial has two positive roots under ${A}_{\mathrm{tot}}=1$, namely 0.07 and 0.537. For ${K}_{2}=0.2,{K}_{3}=0.3$, the polynomial is $p({a}_{0})=4{a}_{0}^{3}23.400{a}_{0}^{2}+2.080{a}_{0}1.664$, and although $\alpha >0$, it has no root in the interval $(0,1)$.
An analogous reasoning holds for the irreversible system.
Appendix 7
Cosubstrate cycling along with metabolite cycling – mimicking the case seen in nitrogen assimilation
In this section, we consider a scenario of intertwined, cosubstrate cycling within a pathway. In particular, we consider a simplified version of the case seen in nitrogen assimilation where NAD(P)H cycling cooccurs together with a cycling from alphaketoglutarate to glutamate (mediated by glutamate dehydrogenase) and back from glutamate to alphaketoglutarate and glutamine (mediated by glutamate synthase). The simplified reaction scheme we consider comprises:
Denoting ${A}_{0}=$ NADH, ${A}_{1}=$ NAD, ${A}_{2}=$ NH3, the ODE system is:
There are three conservation laws:
From $\frac{d{a}_{2}}{dt}+\frac{d{m}_{1}}{dt}+\frac{d{m}_{2}}{dt}=0$ we get
as expected. From $\frac{d{m}_{2}}{dt}=0$ and $\frac{d{a}_{2}}{dt}=0$ we get
From the second and third conservation laws, we get
This gives that for a positive steady state, we need $k}_{\text{in}}<{E}_{1},{k}_{\text{in}}<{E}_{2$ and $\mathrm{m}\mathrm{a}\mathrm{x}({A}_{\mathrm{t}\mathrm{o}\mathrm{t}}W,0)<{a}_{0}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}$. Note that ${A}_{\mathrm{tot}}W$ is the constant value of ${a}_{0}{m}_{0}$ along trajectories.
Plugging these expressions into the first conservation law, we have that steady states are in onetoone correspondence with the solutions to
We first note that the derivative of the right hand side of (38) with respect to a_{0} is:
As ${A}_{\mathrm{t}\mathrm{o}\mathrm{t}}{a}_{0}>0$ and ${a}_{0}{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}W>0$, this function is positive in the interval of interest. Therefore, the right hand side of (38) is an increasing function of a_{0} when $\mathrm{m}\mathrm{a}\mathrm{x}({A}_{\mathrm{t}\mathrm{o}\mathrm{t}}W,0)<{a}_{0}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}$. It follows that if (38) has a solution, it has exactly one.
Rewriting (38) as a polynomial equation, steady states are in onetoone correspondence with the roots in the interval $\mathrm{m}\mathrm{a}\mathrm{x}({A}_{\mathrm{t}\mathrm{o}\mathrm{t}}W,0)<{a}_{0}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ of the following polynomial
When ${a}_{0}={A}_{\mathrm{tot}}$, the polynomial is positive.
Case 1: ${A}_{\mathrm{tot}}\mathrm{}W\mathrm{\le}\mathrm{0}$. In this case we want $0<{a}_{0}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}$. If the independent term of $p$ is negative, then $p$ has exactly one solution in $(0,{A}_{\mathrm{tot}})$. So, if $MT({E}_{2}{k}_{\text{in}}){K}_{2}{k}_{\text{in}}<0$, we have one positive steady state. This is equivalent to
If the independent term is positive or zero, then we note that the degree 1 term is also positive. So either $p$ has all coefficients positive and no positive roots, or the leading term is negative, in which case there is one root larger than ${A}_{\mathrm{tot}}$. Therefore, no positive steady states in this case.
Case 2: ${A}_{\mathrm{t}\mathrm{o}\mathrm{t}}W>0$. In this case, we want $A}_{\mathrm{t}\mathrm{o}\mathrm{t}}W<{a}_{0}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}$. We find that $p({A}_{\mathrm{tot}}W)$ is
If this is negative, then there is one positive steady state, as $p$ is positive at ${A}_{\mathrm{tot}}$. The condition is
Note that $M+W{A}_{\mathrm{tot}}={a}_{0}+{m}_{1}+{m}_{0}+{a}_{1}{a}_{0}{a}_{1}={m}_{0}+{m}_{1}$, and hence needs to be positive. This is assumed here.
If $p({A}_{\mathrm{tot}}W)\ge 0$, then we are in the situation where $p$ is nonnegative both at ${A}_{\mathrm{tot}}W$ and ${A}_{\mathrm{tot}}$, so, if there are roots in the interval $({A}_{\mathrm{tot}}W,{A}_{\mathrm{tot}})$, there must be two. This contradicts that we already showed that there was at most one. So, in this case, no steady states.
To summarize, there is one positive steady state if and only if $k}_{\text{in}}<{E}_{1$ and
or
Appendix 8
Analysis of existing flux data against predicted limits
To support the presented theory, we have analysed existing flux data – compiled from experiments and using flux balance analysis modelling – against predicted limits. The details and main results of this analysis are presented in the main text, under the results and methods sections. Here, we provide further analysis results
Appendix 9
Pathway branching into two pathways with independent cosubstrates
We consider a scenario where two pathways share a common upstream metabolite, as shown in the motif in Figure 4A in the main text. Each branch has its own conserved moiety that is cycled and all reactions are reversible. The reaction system is as follows:
The ODE system is:
While it is not possible to solve this system directly, we examine through simulations the effect of varying the influx of the shared upstream metabolite, k_{in} along with the ratio of the pool sizes, ${A}_{\mathrm{tot}}/{B}_{\mathrm{tot}}$ in the main text (Figure 4 and Appendix 9—figure 1), while varying k_{in} with the ratio of moiety backcycling rates ${k}_{a}/{k}_{b}$ is presented in Appendix 9—figure 2.
Appendix 10
Independent pathways coupled by cosubstrate cycling
We consider a scenario with two pathways that when isolated take the form examined in Appendix 3. However, we now suppose that they are coupled by the conserved moiety, that is ${A}_{0}\to {A}_{1}$ in the ‘forward’ direction of the first pathway, while ${A}_{1}\to {A}_{0}$ in the ‘forward’ direction of the other. The metabolite reactions are reversible, and there is pathwayindependent cycling of the shared conserved moiety. The motif is shown in Figure 5A in the main text, and the reactions are as follows:
The resulting system of ODEs is:
This system has one conservation law
giving that at steady state, ${a}_{1}={A}_{\mathrm{tot}}{a}_{0}$. By solving the steady state equations $\frac{d{m}_{0,1}}{dt}=\frac{d{m}_{1,1}}{dt}=\frac{d{m}_{0,2}}{dt}=\frac{d{m}_{1,2}}{dt}=0$ for ${m}_{0,1},{m}_{0,2},{m}_{1,1},{m}_{1,2}$, we obtain
The expressions are positive provided $E}_{1,1}>{k}_{\mathrm{i}\mathrm{n},1$, $E}_{1,2}>{k}_{\mathrm{i}\mathrm{n},2$, and ${a}_{0},{a}_{1}>0$.
Finally, we use $\frac{d{a}_{0}}{dt}=0$ to solve for a_{0} and obtain
For all quantities to be positive at steady state, we require $0<{a}_{0}<{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}$. By subtracting from ${A}_{\mathrm{tot}}$ the value of a_{0} at steady state, we obtain
To summarize, there is a positive steady state if and only if $E}_{1,1}>{k}_{\mathrm{i}\mathrm{n},1$, $E}_{1,2}>{k}_{\mathrm{i}\mathrm{n},2$, and (44) and (43) are positive. For(44) to be positive,
and either
or
By analysing these cases, we obtain all possible scenarios for a positive steady state to exist, and these are dictated by the difference ${k}_{\mathrm{in},1}{k}_{\mathrm{in},2}$.
For example, if ${k}_{\mathrm{in},1}={k}_{\mathrm{in},2}$, then (45) and (46) hold directly. If ${K}_{a}={L}_{a}$, then we require $\frac{{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}{F}_{a}}{{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}+{L}_{a}}>{k}_{\mathrm{i}\mathrm{n},1}{k}_{\mathrm{i}\mathrm{n},2}>\frac{{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}{E}_{a}}{{A}_{\mathrm{t}\mathrm{o}\mathrm{t}}+{K}_{a}}$ to hold.
If $K}_{a}>{L}_{a$, the conditions lead to either
or
If $K}_{a}<{L}_{a$, the conditions lead to either
or
A key consequence of these conditions is that coupled pathways can admit higher influx values without upstream metabolite build up, due to the cycling enzyme condition now depending on the difference between the ${k}_{\text{in}}$ values rather than the values themselves. This occurs when the limit due to the pathway enzyme is larger than that of the cycling enzyme, so there is always a range of ${A}_{\mathrm{tot}}$ where this applies. An example is shown in Appendix 10—figure 1.
To complement this steady state analysis with dynamics of this system, we use numerical simulations to study the system. In particular, we consider the effect of randomly fluctuating influxes on the downstream metabolites. The analysis is achieved by fixing the average of the influxes, while drawing the logratio (i.e.${\mathrm{log}}_{10}{k}_{\mathrm{in},1}/{k}_{\mathrm{in},2}$) from a standard normal distribution with mean $\mu =0$ and variance ${\sigma}^{2}=1$. A new logratio is drawn after waiting a time that is drawn from an exponential distribution with mean $\tau $. Example timeseries are shown in Figure 5(C and D) of the main text. The logratio is chosen as the variable instead of simply the ratio as it allows us to examine large variations, while keeping the effect on each pathway symmetric. This random process can be thought of as the discretetime analogue of the OrnsteinUhlenbeck process: it has the same steady state distribution, mean, variance and correlation function as its continuoustime counterpart, while being much easier and faster to implement as part of a larger system where the other equations are deterministic.
The effect of different values of $\tau $ and total pool size ${A}_{\mathrm{tot}}$ is shown in Figure 5 of the main text, and is further explored in Appendix 10—figure 2. These analyses show that the effect of increasing the pathway independent moiety cycling rate is to reduce the critical pool size above which the downstream metabolites are anticorrelated. Furthermore, this behaviour is observed whether the metabolite reactions are irreversible or reversible. Since the transition from correlated to anticorrelated is so sharp, relatively small changes in the total amount of cycled cosubstrate, combined with noisy influx rates, could lead to the downstream metabolites changing from being correlated to anticorrelated, or viceversa.
Appendix 10—figure 2 also further demonstrates the effect of changing ${E}_{a}$, $\tau $ and ${A}_{\mathrm{tot}}$ on the time series of the products ${M}_{1,1}$ and ${M}_{1,2}$ and the cosubstrate ratio ${A}_{0}/{A}_{1}$. For parameters represented by the black dot in the first row (panels (E)  (H)), the products are correlated when ${E}_{a}\le 1$, and become uncorrelated when ${E}_{a}=10$. The system is still responsive to the noise, but with less variation compared to the other rows because it is limited by the low ${A}_{\mathrm{tot}}$. For parameters represented by the grey dot in the first row (panels (I)  (L)), the products switch from correlated to anticorrelated as ${E}_{a}$ increases, and remain responsive to the noise. The variation in ${A}_{0}/{A}_{1}$ reduces from left to right because ${E}_{a}$ increases. For parameters represented by the white dot in the first row (panels (M)  (P)), the products remain anticorrelated and responsive to the noise. The variation in ${A}_{0}/{A}_{1}$ reduces from left to right but is not as pronounced for small ${E}_{a}$ because ${A}_{\mathrm{tot}}$ is very large.
Appendix 11
Table of parameters
Data availability
All used data and models for cases shown in Figure 4 and 5 are made available via a dedicated repository (https://doi.org/10.5281/zenodo.7565439) and the following Github page: https://github.com/OSSLab/CoSubstrateDynamics/tree/v1.0.0.
References

Recon3D enables a threedimensional view of gene variation in human metabolismNature Biotechnology 36:272–281.https://doi.org/10.1038/nbt.4072

Brenda, the elixir core data resource in 2021: new developments and updatesNucleic Acids Research 49:D498–D508.https://doi.org/10.1093/nar/gkaa1025

Tumor cell energy metabolism and its common features with yeast metabolismBiochimica et Biophysica Acta 1796:252–265.https://doi.org/10.1016/j.bbcan.2009.07.003

Metabolic dynamics restricted by conserved carriers: jamming and feedbackPLOS Computational Biology 13:e1005847.https://doi.org/10.1371/journal.pcbi.1005847

Kinetic parameters of enzymatic reactions in states of maximal activity; an evolutionary approachJournal of Theoretical Biology 151:249–283.https://doi.org/10.1016/s00225193(05)803637

Mathematical analysis of enzymic reaction systems using optimization principlesEuropean Journal of Biochemistry 201:1–21.https://doi.org/10.1111/j.14321033.1991.tb16251.x

Dynamic behaviors of an open substrate cycle: a graphical approachJournal of Theoretical Biology 140:399–416.https://doi.org/10.1016/S00225193(89)800955

Metabolic control analysis of moietyconserved cyclesEuropean Journal of Biochemistry 155:631–640.https://doi.org/10.1111/j.14321033.1986.tb09534.x

The glycolytic flux in Escherichia coli is controlled by the demand for atpJournal of Bacteriology 184:3909–3916.https://doi.org/10.1128/JB.184.14.39093916.2002

Selfreplenishment cycles generate a threshold responseScientific Reports 9:17139.https://doi.org/10.1038/s41598019535891

Simple constrainedoptimization view of acetate overflow in E. coliBiotechnology and Bioengineering 35:732–738.https://doi.org/10.1002/bit.260350711

Dynamic behavior of cyclic enzyme systemsJournal of Theoretical Biology 104:591–598.https://doi.org/10.1016/00225193(83)902473

A comprehensive genomescale reconstruction of Escherichia coli metabolism2011Molecular Systems Biology 7:535.https://doi.org/10.1038/msb.2011.65

The quantitative and conditiondependent Escherichia coli proteomeNature Biotechnology 34:104–110.https://doi.org/10.1038/nbt.3418

The optimal efficiency and the economic degrees of coupling of oxidative phosphorylationEuropean Journal of Biochemistry 109:269–283.https://doi.org/10.1111/j.14321033.1980.tb04792.x

The danger of metabolic pathways with turbo designTrends in Biochemical Sciences 23:162–169.https://doi.org/10.1016/s09680004(98)012055

Symbolic proof of bistability in reaction networksSIAM Journal on Applied Dynamical Systems 20:1–37.https://doi.org/10.1137/20M1326672

Overflow metabolism in Escherichia coli during steadystate growth: transcriptional regulation and effect of the redox ratioApplied and Environmental Microbiology 72:3653–3661.https://doi.org/10.1128/AEM.72.5.36533661.2006

Control of ATP homeostasis during the respiro‐fermentative transition in yeastMolecular Systems Biology 6:344.https://doi.org/10.1038/msb.2009.100
Article and author information
Author details
Funding
Biotechnology and Biological Sciences Research Council (BB/T010150/1)
 Robert West
Novo Nordisk (F18OC0052483)
 Elisenda Feliu
Gordon and Betty Moore Foundation (GBMF9200)
 Orkun S Soyer
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This project is funded by the Biotechnology and Biological Sciences Research Council (BBSRC) (grant BB/T010150/1). EF acknowledges funding from the Novo Nordisk Foundation (grant NNF18OC0052483), while OSS acknowledges support from the Gordon and Betty Moore Foundation (grant https://doi.org/10.37807/GBMF9200).
We would like to thank Wenying Shou for constructive comments on an earlier version of this manuscript, and Dan Davidi for his help with datasets of reaction fluxes and enzyme abundances.
Version history
 Preprint posted: September 6, 2022 (view preprint)
 Received: October 22, 2022
 Accepted: February 16, 2023
 Accepted Manuscript published: February 17, 2023 (version 1)
 Version of Record published: March 20, 2023 (version 2)
Copyright
© 2023, West, Delattre et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 792
 views

 175
 downloads

 1
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
Revealing protein binding sites with other molecules, such as nucleic acids, peptides, or small ligands, sheds light on disease mechanism elucidation and novel drug design. With the explosive growth of proteins in sequence databases, how to accurately and efficiently identify these binding sites from sequences becomes essential. However, current methods mostly rely on expensive multiple sequence alignments or experimental protein structures, limiting their genomescale applications. Besides, these methods haven’t fully explored the geometry of the protein structures. Here, we propose GPSite, a multitask network for simultaneously predicting binding residues of DNA, RNA, peptide, protein, ATP, HEM, and metal ions on proteins. GPSite was trained on informative sequence embeddings and predicted structures from protein language models, while comprehensively extracting residual and relational geometric contexts in an endtoend manner. Experiments demonstrate that GPSite substantially surpasses stateoftheart sequencebased and structurebased approaches on various benchmark datasets, even when the structures are not wellpredicted. The low computational cost of GPSite enables rapid genomescale binding residue annotations for over 568,000 sequences, providing opportunities to unveil unexplored associations of binding sites with molecular functions, biological processes, and genetic variants. The GPSite webserver and annotation database can be freely accessed at https://bioweb1.nsccgz.cn/app/GPSite.

 Cell Biology
 Computational and Systems Biology
Computer models of the human ventricular cardiomyocyte action potential (AP) have reached a level of detail and maturity that has led to an increasing number of applications in the pharmaceutical sector. However, interfacing the models with experimental data can become a significant computational burden. To mitigate the computational burden, the present study introduces a neural network (NN) that emulates the AP for given maximum conductances of selected ion channels, pumps, and exchangers. Its applicability in pharmacological studies was tested on synthetic and experimental data. The NN emulator potentially enables massive speedups compared to regular simulations and the forward problem (find drugged AP for pharmacological parameters defined as scaling factors of control maximum conductances) on synthetic data could be solved with average rootmeansquare errors (RMSE) of 0.47 mV in normal APs and of 14.5 mV in abnormal APs exhibiting early afterdepolarizations (72.5% of the emulated APs were alining with the abnormality, and the substantial majority of the remaining APs demonstrated pronounced proximity). This demonstrates not only very fast and mostly very accurate AP emulations but also the capability of accounting for discontinuities, a major advantage over existing emulation strategies. Furthermore, the inverse problem (find pharmacological parameters for control and drugged APs through optimization) on synthetic data could be solved with high accuracy shown by a maximum RMSE of 0.22 in the estimated pharmacological parameters. However, notable mismatches were observed between pharmacological parameters estimated from experimental data and distributions obtained from the Comprehensive in vitro Proarrhythmia Assay initiative. This reveals larger inaccuracies which can be attributed particularly to the fact that small tissue preparations were studied while the emulator was trained on single cardiomyocyte data. Overall, our study highlights the potential of NN emulators as powerful tool for an increased efficiency in future quantitative systems pharmacology studies.