Dynamics of co-substrate pools can constrain and regulate metabolic fluxes

Cycling of co-substrates, whereby a metabolite is converted among alternate forms via different reactions, is ubiquitous in metabolism. Several cycled co-substrates 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 co-substrates in different parts of central metabolism. Here, we develop a mathematical framework to analyse the effect of co-substrate cycling on metabolic flux. In the cases of a single reaction and linear pathways, we find that co-substrate 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 co-substrate. 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 co-substrate cycling. We discuss how the limitations imposed by co-substrate cycling provide experimentally testable hypotheses on specific metabolic phenotypes. We conclude that measuring and controlling co-substrate dynamics is crucial for understanding and engineering metabolic fluxes in cells.


Introduction
Dynamics of cell metabolism directly influences individual and population-level 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;Diaz-Ruiz et al., 2009;Carmona-Fontaine et al., 2013) and cell-to-cell cross-feeding (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 . The resulting theory predicted a trade-off 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 max app ) 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 ( kcat ) (Davidi et al., 2016). However, individual estimates from each condition (i.e. individual kapp 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 so-called co-substrate 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 Cornish-Bowden, 2000). For example, the total pool of ATP eLife 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 co-substrates, 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 co-substrate 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 co-substrate 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 co-substrates.
They found that co-substrate 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 co-substrates in the cell: changing the number of co-substrates 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 co-substrates in different cell types and conditions. 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 switch-like and bistable dynamics (Okamoto and Hayashi, 1983;Hervagault and Cimino, 1989), and that metabolite cycling motifs introduce total co-substrate 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 'turbo-design', 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 co-substrate 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 co-substrate cycling in cell metabolism dynamics, we undertook here an analytical and simulation-based mathematical study together with analyses of measured fluxes. We created models of enzymatic reaction systems featuring co-substrate cycling, abstracted from real metabolic systems such as glycolysis, nitrogen-assimilation, and central carbon metabolism. We found that co-substrate 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 co-substrate. Analysing measured fluxes, we find that several of the co-substrate 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 co-substrate cycling. In addition to its possible constraining role, we show that co-substrate cycling can also act as a regulatory element, where control of co-substrate pool size can allow control of flux dynamics across connected or branching pathways. Together, these findings show that co-substrate 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 co-substrate pool sizes and turnover dynamics, and can be expanded to explain dynamic measurements of metabolite concentrations in different perturbation experiments.

Results
Co-substrate cycling represents a ubiquitous motif in metabolism with co-substrate 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 co-substrate 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 co-substrates. Examples for the latter include NADPH, FADH2, GTP, and Acetyl-CoA and their corresponding alternate forms, while examples for the former include the tetrahydrofolate (THF) / 5,10-Methylene-THF and glutamate / α -ketoglutarate (akg) pairs involved in one-carbon 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 genome-scale 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 co-substrate involving cycling reactions can be abstracted as a simplified motif as shown in ( Figure 1A). This abstract representation highlights the fact that the total pool-size 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 co-substrates, 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 co-substrate 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 inter-conversions among pools rather than through biosynthesis (Walther et al., 2010).

Co-substrate cycling introduces a limitation on reaction flux
To explore the effect of co-substrate 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 co-substrate Heatmap of the steady state concentration of M 0 as a function of the total co-substrate pool size ( Atot ) and inflow flux (k in ). White area shows the region where there is no steady state. On both panels, the dashed line indicates the limitation from the primary enzyme, k in < V max,E 0 , and the solid line indicates the limitation from co-substrate cycling, k in < AtotV max,Ea /(K M,Ea + Atot) . In panel (C), there is a range of Atot values for which the first limitation is more severe than the second. In contrast, in panel (D), the second limitation is always more severe than the first. In (B & C) the parameters used for the primary enzyme (for the reaction converting M 0 into M 1 ) are picked from within a physiological range (see Supplementary file 1) and are set to: Etot = 0.01 mM, kcat = 100/s , K M,E 0 = K M,Ea = 50 µM , while k out is set to 0.1/s. The Etot and k cat for the co-substrate cycling enzyme are 1.2 times those for the primary enzyme. In panel (D) the parameters are the same except for the Etot and k cat of the co-substrate cycling enzyme, which are set to 0.7 times those for the primary enzyme. cycling ( Figure 1A). For co-substrate cycling, we consider additional 'background' enzymatic reactions that are independent of M 0 and can also convert the co-substrate (denoted Ea 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 in-flux of M 0 , and out-flux of M 1 , either in-and-out of the cell or from other pathways, and Atot 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 max,E0 and V max,Ea are the maximal rates (i.e. Vmax = kcatEtot ) 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 Ea ), respectively, while the parameters K M,E0 and K M,Ea are the individual or combined Michaelis-Menten coefficients for these enzymes' substrates (i.e. for A 0 and M 0 and A 1 , respectively). The term α is (in this case where all reactions are irreversible) equal to V max,Ea − k in , and in general is a positive expression comprising k in , and the Michaelis-Menten coefficients and the Vmax 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 co-substrate cycling, the production rate of M 0 is given by k in , while its consumption rate is a function of the V max,E0 and the concentration of A 0 . In turn, the concentration of A 0 is determined by its re-generation rate (which is a function of K M,Ea and V max,Ea ) and the pool size ( Atot ). This explains the inequalities given in (2) and shows that a cycled co-substrate 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 co-substrate cycling leads to only the constraint k in < V max,E0 , while considering it as a non-enzymatic reaction with co-substrate cycling only, leads to only the constraint k in < AtotV max,Ea /(K M,Ea + Atot) becoming the sole limitation on the system (see Appendix 3). In other words, the two limitations act independently.
To conclude this section, we re-iterate its main result. The flux of a reaction involving co-substrate cycling is limited either by the kinetics of the primary enzyme mediating that reaction, or by the turnover rate of the co-substrate. The latter is determined by the co-substrate pool size and the kinetics of the enzyme(s) mediating its turnover.
Co-substrate 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 co-substrate cycling, for example as seen in upper glycolysis (Appendix 1- figure 1A). In this model, we only consider intra-pathway metabolite cycling, i.e. the co-substrate is consumed and re-generated 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 co-substrate 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 Atot 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) intra-pathway 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 re-generation 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 co-substrate pool size and the enzyme kinetics. substrates that are produced before the first co-substrate cycling reaction, and continued decline of A 0 ). The parameters used are picked from within a physiological range (see 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. Vmax ) (Davidi et al., 2016). This finding is based on calculating Vmax 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 co-substrate pool size and turnover dynamics.
To explore this possibility, we re-analysed 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 enzyme-flux-k cat values with full experimental data and 259 combinations with only FBA-predicted 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 co-factor 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 Michaelis-Menten 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 ∆rG ′ is calculated assuming that all reactants are at 1 mM and pH = 7) against the normalized fluxv/(E 0 · kcat) (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 co-substrate pool sizes. Co-substrate pool sizes do change between different conditions, and we note that such changes cannot be due to flux changes in co-substrate utilising reactions. But, on the other hand, changes in pool size can affect flux in those reactions, where co-substrate dynamics is limiting (as predicted by the theory). For both measured and FBA-predicted fluxes, we find that several reactions show significant correlation between flux and co-substrate pool size (see Figure 3B-D, see also Appendix 8-table 1 and Appendix 8- figure 4). In the case of FBA-predicted fluxes, however, we note that these results can be confounded due to additional, flux-to-flux 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 glucose-6-phosphate dehydrogenase (G6PDH) linked with the NADPH pool.
In summary, these results show that for reactions involving co-substrate cycling (1) measured fluxes are lower than those predicted by kinetics of the primary enzyme (i.e. enzyme involved in substrate  Davidi et al., 2016;Gerosa et al., 2015) plotted against the calculated primary enzyme kinetic threshold (first part of eq. (1)). Notice that there are 7 points for each reaction, corresponding to the different experimental conditions under which measurements or FBA modelling was done (see Supplementary file 1 for data, along with reaction names and metabolites involved). (B-D) Measured flux values under different experimental conditions (from Gerosa et al., 2015) for select reactions plotted against the corresponding co-substrate pool size. Panels B to D show reactions for phosphoglycerate kinase (PGK), malate dehydrogenase (MDH), and glucose-6-phosphate dehydrogenase (G6PDH). Each point on these panels is a separate flux measurement under a different environmental condition, where the co-substrate pool size is also measured. Error bars represent standard deviations of flux and metabolite measurements as they appear in the dataset from Gerosa et al. Point colours represent co-substrate type and are as shown in the legend to panel A. Lines show the best linear fit with the corresponding normalised RMSE shown in the panel title. conversion) alone, and (2) there is -for some reactions -a correlation between flux and co-substrate pool size. Both observations could indicate co-substrate pool sizes and/or co-substrate 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 co-substrate cycling within the study of metabolic system dynamics.

Co-substrate cycling allows regulation of branch point fluxes
In addition to its possible constraining effects on fluxes, we wondered if co-substrate dynamics can offer a regulatory element in cellular metabolism. In particular, co-substrate 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 co-substrate, 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 co-substrate pool sizes, Atot and Btot , could allow regulation of the fluxes on the two branches. To test this hypothesis, we ran numerical simulations with different co-substrate 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 Atot to Btot ( 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 build-up can affect either branch depending on which co-substrate has the lower The two branches are linked to separate co-substrate pools, A and B . Note that pathway independent turnover of the cosubstrates is included in the model (see Figure 4-source code 1). (B) The pathways' flux ratio (i.e. flux into M 2,2 divided by flux into M 2,1 ) shown in colour mapping, against the ratio of co-substrate pool sizes, Atot and Btot , and the influx rate, k in , into the upstream metabolite. In the block colour areas, the system has no steady state and the indicated metabolite(s) M 0 and one of the metabolites M 1,2 or M 1,1 accumulate towards infinity.  Figure 4B). There is also a regime of only the upstream, branch point metabolite building-up, but this happens only when all reactions are considered as reversible and the extent of it depends on turnover rates of the two co-substrates (Appendix 9-figure 1).
In the no-build-up, steady state regime, changing the pool size ratio of the two co-substrates 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).
Inter-pathway co-substrate 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 co-substrate 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 Figure 5. Motif, heatmap and time-series for the coupled pathway model with noisy influxes. (A) Cartoon representation of two pathways coupled via the same co-substrate cycling. The two forms of the co-substrate are indicated as A 0 and A 1 . It is converted from A 0 to A 1 on the lower pathway, and from A 1 to A 0 in the upper pathway. The presented results are for a model with reversible enzyme kinetics, while the results from a model with irreversible enzyme kinetics are shown in Appendix 10-figure 2. (B) Correlation coefficient of the two pathway product metabolites, M 1,2 and M 1,1 , as a function of the total amount of co-substrate ( Atot ) and the extent of fluctuations in the two pathway influxes, k in,1 and k in,2 . The influx fluctuation is characterised by a waiting time that is exponentially distributed with mean τ , after which the log ratio of the k in values is drawn from a standard normal distribution.
The mean of the k in values is set to be 0.1 and the pathway-independent cycling occurs at a much lower rate compared to the other reactions (see Figure 5-source code 1). (C) Concentrations of metabolites M 1,2 (green) and M 1,1 (magenta), pathway influx ratio (pink), and A 0 /A 1 ratio (blue) as a function of time. The simulation is run with parameters corresponding to the grey dot in (B) where the products are correlated, and the rate of k in fluctuations is on a similar timescale to the other reactions. The system is largely unresponsive to the noise, and the M 1,2 and M 1,1 are fully correlated (i.e. the green and magenta curves overlap). (D) Concentrations of metabolites M 1,2 (green) and M 1,1 (magenta), pathway influx ratio (pink), and A 0 /A 1 ratio (blue) as a function of time. The simulation is run with parameters corresponding to the black dot in (B) where the products are correlated, but the fluctuations in k in values occur at a much lower rate than the other reactions. The system is responsive to the noise, yet the M 1,2 and M 1,1 are fully correlated (i.e. the green and magenta curves overlap). For both simulations, all kinetic parameters are arbitrarily set to 1, apart from the pathwayindependent co-substrate recycling ( V max,Ea ) that is set to 0.01 (see Figure 5-source code 1).
The online version of this article includes the following source code for figure 5: Source code 1. Python implementation of connected pathway model, presented in Figure 5. phosphate pathway and some amino acid biosynthesis pathways (notably Methionine), through NADPH generation and consumption respectively. We hypothesised that such inter-pathway co-substrate cycling might cause; (1) the co-substrate 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 pathway-independent turn-over of the co-substrate (see Appendix 10). In other words, for two pathways coupled via co-substrate cycling, the cycling-dependent 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 steady-state outputs of the pathways with random fluctuations in their influx ( Figure 5B). We found that there is either a high level of anti-correlation 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 anti-correlation 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 sub-regimes. 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 non-responsive 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 co-substrate cycling dynamics. Because the turnover of co-substrate 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 max,Ea = 0.01 are representative of those for V max,Ea = 0 (see Appendix 10-figure 2). As we increase the rate of the assumed background reaction, that is pathway-independent turnover of the co-substrate, we find that these results remain qualitatively the same, but the transition point from anti-correlation in outputs to correlation, shifts to lower Atot values (Appendix 10- figure 2).
These results show that coupling by co-substrate 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 co-substrate 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 co-substrate (Stucki, 1980).

Discussion
We presented a mathematical analysis of metabolic systems featuring co-substrate 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 co-substrate cycling, increase towards infinity and the system cannot reach steady state. Specifically, for reactions involving co-substrates, 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 co-substrate and the pool size of that co-substrate.
This second, additional constraint arising from co-substrate cycling can be directly relevant for cell physiology. We particularly note that this threshold can be highly dynamic, and condition-and cell-dependent. When cells have a permanently or occasionally lowered total co-substrate pool size (i.e. lower Atot ), 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 Atot can be adjusted in the long term, for example by reducing substrate transporter expression or increasing co-substrate 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 co-substrate linked reactions could indeed by limited by co-substrate 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 co-substrate 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 co-substrate linked reactions. In this quest, we expect that the expansion of measurements to eukaryotic cells to be particularly challenging due to organelle-specific pools, but some progress is being made to achieve at least mitochondrial and cytosolic measurements (Chen et al., 2016). Despite the current limitations, our data-based analyses highlighted three key reactions, that are possibly limited by co-substrate dynamics, and that are mediated by phosphoglycerate kinase (PGK), malate dehydrogenase (MDH), and glucose-6-phosphate dehydrogenase (G6PDH) and linked to ATP, NADH, and NADPH pools. Possible flux limitation of these reactions by co-substrate 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 'substrate-induced death' (van Heerden et al., 2014) or 'overflow metabolism'. The latter usually refers to a respiration-to-fermentation switch under respiratory conditions (e.g. the Warburg and Crabtree effects [Warburg, 1956;Diaz-Ruiz 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 respiration-to-fermentation switch (Diaz-Ruiz 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 substrate-induced death and metabolite excretions could relate to increasing substrate influx rate reaching close to the limits imposed by co-substrate dynamics. There is experimental support for this hypothesis in the case of both observations. Substrate-induced 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 respiration-to-fermentation 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, sulfur-compound excretions are linked to alterations in the NADPH pool through changes in the amino acid metabolism (Olin-Sandoval et al., 2019;Green et al., 2020).
Dynamical thresholds relating to co-substrate pools would be relevant for all co-substrates, 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 co-substrate 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 co-substrate consumption and production, or by altering co-substrate 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 co-substrate consumption or production (e.g. ATPases, oxidases, or other) and controlling their expression. It would also be possible to monitor co-substrate 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 co-substrate 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 inter-connected pathways, thereby ensuring correlation between outputs of different metabolic processes, or regulating flux across branch points. Regulation of fluxes through co-substrate 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, S-adenosylmethionine, 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 self-regulation and the pool sizes of key co-substrates could be seen as 'tuning points' controlling a more complex metabolic system. We thus propose further experimental analyses focusing on co-substrate pool sizes and turnover dynamics to understand and manipulate cell physiology.

Materials and methods
Model of a single reaction with co-substrate cycling The metabolic system shown in Figure 1A comprises the following biochemical reactions: where metabolites are denoted by M i and the different forms of the co-substrate 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 Michaelis-Menten 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 Michaelis-Menten coefficients of the enzyme for its substrates, and Vmax is the total enzyme concentration times its catalytic rate (i.e. Vmax = kcatEtot ) (see Appendix 11-table 1 for a list of parameters and their units). We further have the conservation relation a 0 + a 1 = Atot , where Atot 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 re-arrangements 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.

Reaction fluxes and enzyme kinetic parameters
To support the model findings on co-substrate pools acting as a possible limitation on reaction fluxes, we analysed measured and FBA-derived flux data collated previously (Davidi et al., 2016;Gerosa et al., 2015). We focused our analyses on reactions involving co-substrates only. We compared measured (or FBA-derived) 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 proteomics-based 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ía-Contreras et al., 2012). When comparing flux levels against co-substrate pool sizes, we used the matching, measured pool-sizes 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).

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.

Additional files
Supplementary files • Supplementary file 1. Enzyme kinetics, flux, metabolite concentration, and enzyme abundance data associated with flux analyses.
The following dataset was generated: Author (

Appendix 1
The mathematical approach and modelling setting As explained in the main text, we are interested in understanding the effect of co-substrate cycling on the flux through metabolic pathways, such as those shown in Appendix 1-figure 1.
In these Appendices, we use a generic co-substrate pair, denoted as A 0 , A 1 . We consider the synthesis or degradation of the co-substrate pair, or consider it as a conserved moiety, i.e. having a fixed total concentration. Our generic co-substrate pair, A 0 and A 1 , can be taken as representing a specific co-substrate, such as NAD(H), but note that the mathematical analyses presented would be applicable to any co-substrate 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 i−1 ⇌M i , while in the latter case, they are shown as M i−1 −→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 co-substrate 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 co-substrate considered here, we have: We talk about irreversible co-substrate conversion, if k 5 = 0 or k 6 = 0 , that is, only conversion in one direction is considered. We talk about no free co-substrate 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 Michaelis-Menten (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 co-substrate 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 co-substrate 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 i−1 ⇌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 K i = k on,1 k on,2 + k off k on,2 + k off k off,1 kon(k on,2 + k off,1 + k on,1 ) and L i−1 = k on,1 k on,2 + k off k on,2 + k off k off,1 k off,2 (k off,1 + k on,1 + k off ) .
When there are two substrates that take part in the reaction, the k on and k 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 kcat+ and k cat− , multiplied by the total enzyme concentration (denoted Etot , below), and are given by E i = Etot k on,1 k on,2 k on,2 + k off,1 + k on,1 and F i = Etot k off k off,1 k off + k off,1 + k on,1 .
For the irreversible enzymatic reaction, the reaction schemes simplify to: Appendix 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 Michaelis-Menten 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 non-enzymatic 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 Michaelis-Menten 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 well-known Descartes' rule of signs.
Lemma 1. Let p(x) 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 (0, T) 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 τ of sign changes in the sequence of coefficients ignoring zero coefficients, and the difference between τ 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 + + − , If the independent term is positive or zero, then the sign sequence is one of + + − , + − − , + 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 +∞ , the root must be in the interval (T, +∞) .
If the independent term is negative, then the sign sequence is − + − . The polynomial is negative both at 0 and at +∞ . As p(T) > 0 , there must be a positive root in (0, T) and one in (T, +∞) , and there cannot be more by the Descartes' rule of signs. From this, the statement of the lemma follows.

Appendix 2 Considering co-substrate pool size
We first consider cycling of a generic co-substrate 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 = [ A 0 ] and a 1 = [ A 1 ] , 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: , a 1 = ks(k d + 2k 5 ) k d (k d + k 5 + k 6 ) , 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 co-substrate conversion rates, it can be safely neglected in the brackets, resulting in the expression of steady state formulas as: a 0 = 2k 6 ks/k d k 5 + k 6 , a 1 = 2k 5 ks/k d k 5 + k 6 , 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 = Tk 6 /(k 5 + k 6 ) and a 1 = Tk 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 co-substrates is an enzymatic reaction and make the same simplifying assumptions as above that k 1 = k 3 = ks 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:  (Ka − La)).
This is stable as long as all parameters are positive. Note that in the case of Ka = La , 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 non-enzymatic case presented above. Thus, whether the metabolite cycling is considered as a non-enzymatic or enzymatic (i.e. following Michaelis-Menten kinetics) reaction, the co-substrates 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 co-substrate cycling or not, as presented in the main text.

Enzymatic reaction with co-substrate cycling
We first reconsider the case where all reactions, including the off-pathway cycling, are enzymatic (hence are modelled with Michaelis-Menten 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 α , as follows: For the steady state equations given above to be positive, the following conditions must be satisfied: Note that as FaAtot La+Atot < Fa , the second condition readily implies k in < Fa and α is positive. When there is a positive steady state, then it is asymptotically stable.
If the main pathway is irreversible, but the co-substrate reaction is still reversible, the ODEs describing the system dynamics are: The steady state of the system is: where the composite parameter α 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, α is positive, and the positive steady state is asymptotically stable. If the main pathway is irreversible, and the co-substrate 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 α 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 co-substrate cycling, biosynthesis and degradation
To extend the previous analysis we consider the same simple scenario where a pathway involves a single co-substrate consuming reaction and a back reaction re-generating the co-substrate, with the addition of synthesis and degradation of the co-substrate forms. This system is comprised of the following reactions: The resulting system of ODEs is: The steady state concentrations are then given by: a 0 = k 10 (k 7 − k in ) + k 6 (k 7 + k 9 ) k 10 k 5 + k 10 k 8 + k 6 k 8 a 1 = k 8 (k 9 + k in ) + k 5 (k 7 + k 9 ) k 10 k 5 + k 10 k 8 + k 6 k 8 m 0 = K 1 k in ( (F 1 + k in )(k 8 (k 9 + k in ) + k 5 (k 7 + k 9 )) + koutL 0 (k 10 k 5 + k 10 k 8 + k 6 k 8 ) ) These expressions are positive if and only if k in < k 7 + k 6 k 10 (k 7 + k 9 ) and k in < E 1 .
Comparing these results with those of Subsection C.1, we see some similar themes. Firstly, the total amount of co-substrate 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 in < E 1 , and any other conditions take the form k in < f , where f is a function of the parameters controlling synthesis, degradation and the turnover of the co-substrate. Thus, in the analysis that follows, we only consider the case of conserved co-substrate, as adding synthesis and degradation only affects the quantitative results, not the qualitative behaviour.

Enzymatic reaction without co-substrate cycling
Considering an enzymatic reaction without co-substrates, 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 in < 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 in < E 1 and the steady state is asymptotically stable when this holds. From this we see that the condition k in < E 1 for stability in the enzymatic system with co-substrate cycling is a result of the reaction being enzymatic.

Non-enzymatic reaction with co-substrate cycling
Consider a theoretical case of non-enzymatic 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 in < k 4 Atot , 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 in < k 4 Atot and the steady state is asymptotically stable when this holds. Hence, we see that introducing co-substrate cycling always introduces a new constraint into the system, even in this simple case where the cycling is not enzymatic. Linear, arbitrary length, pathway model with co-substrate 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 dm0 dt + · · · + dm3 dt = 0 , dm3 dt = 0 , dm2 dt + dm3 dt = 0 , dm1 dt + dm2 dt + dm3 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 < Atot , 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 one-to-one correspondence with the positive steady states, provided (17) holds. The polynomial is As a 0 = Atot , we have p(Atot) > 0 , and hence, by Lemma 1, the system has positive steady states if and only if α < 0 . Note that at k in = 0 , α = − Atot W(E 3 − k in )(E 2 − k in )L 1 < 0 . Hence, for k in small enough, α < 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 α is a polynomial in k 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 in = E 2 or k in = E 3 , then at this value of k in , α is positive independently of the rest of the parameters. Hence, in the region of interest α is negative only in an interval of the form (0, B) , where B is the smallest positive root of α . The bound for positive steady states is now given as If the minimum is attained at B , then, as k in approaches B , the smallest positive root of p(a 0 ) approaches 0 as α approaches zero. Using this in (16), we obtain that If E 1 < B , then k in needs to be smaller than E 1 . When k in approaches E 1 , a 0 approaches the root of p when k in = E 1 , which is not zero. Then m 0 still approaches ∞ as it has E 1 − k 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 co-substrate recycling in arbitrary places.
To write the equations in a simple format, we write We start with the reversible reaction kinetics. Then the associated ODE system becomes: Kρ+1Lρ+Lρxρ+Kρ+1yρ+1 .
This ODE system has precisely two conservation laws: Note that we have the following equalities: For i = 1, . . . , n − 2 , by solving (21) for x i , we obtain the following recursive formulas: Finally, from (22) and the conservation law a 0 + a 1 = Atot , we obtain These expressions are all positive if and only if 0 < a 0 < Atot and k in < E i for i = 1, . . . , 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 ≤ ℓ < ρ − 1 ≤ n − 1 . Then, it holds We write for shortness where for i = 0, . . . , n − 1 , 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 · · · b j−1 . For i = n − 1 , (26) agrees with (23). Assume the formula is true for some , and we prove it for i − 1 . To do so, we use (23) and the induction hypothesis for x i . For n − 1 > i − 1 ≥ 0 , we have This is (26) for i − 1 . Hence, the formula holds. Finally, we can use that xn = mn = kin kout to obtain: This gives: Remember that b ℓ , b ℓ+1 , bρ, b ρ+1 depend on a 0 , a 1 , while the rest of b 's do not. For a product of the form bu · · · b j with u ≤ j , we have the following: Observe that for all i , the summand of ∑ n−1 j=i z j (b i · · · b j−1 )(c j+1 · · · c n−1 ) 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 ≤ ℓ − 1 , the term (b i · · · b n−1 ) does not depend on a 0 , a 1 , and in the sum ∑ n−1 j=i z j (b i · · · b j−1 )(c j+1 · · · c n−1 ) there are summands involving 1 a0 , for eample when j = ℓ . Hence m i , for i ≤ ℓ − 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 = ℓ , 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 ℓ goes to infinity as well. We conclude that When i ≥ ℓ , then no summand in the numerator of m i in (27-29) that involves 1 a0 , 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  and is a function of the parameters of the system, not involving W .
As a 1 = Atot − a 0 , this is a degree 2 polynomial in a 0 . The leading term comes from the term (0, 0, 1a 0 − W)koutc ℓ+1,n−1 0, 0, 1a 1 , and is −koutc ℓ+1,n−1 a 2 0 , which is negative. The independent term is We divide by kout and define where recall from (25) that When a 0 = Atot , 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 α < 0.
Note that when k in = 0 this inequality holds, as all terms with z j vanish. When k in approaches one of E i with ℓ + 1 < i ≤ n , the negative term of α approaches zero, and the inequality does not hold. For example, for n = 3 , ℓ = 0 and ρ = 2 , α was found in (18). For n = 4 , ℓ = 0 and ρ = 2 , we have Let B be the smallest positive root of α = 0 as a polynomial in k in , if it exists, or take B = ∞ if not. Similarly, let B ′ be the second such root, if it exists, or B ′ = ∞ if not. If the smallest of E j is attained for some E i with ℓ + 1 < i ≤ n , then α is positive at k in = E i . In that case, as α is negative at k in = 0 , there is at least one value of k in < E i for all i such that α = 0 and hence B is finite. This shows that min(E 1 , . . . , En, B) = min(E 1 , . . . , E ℓ+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 in approaches B , the first positive root of the polynomial in (3031, ) approaches zero (as α 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 ℓ + 1 < i ≤ n . If (32) is attained at E i with i ≤ ℓ + 1 , then as k 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 in in the denominator:

Appendix 5
Multiple co-substrate cycling along a single pathway -mimicking the case seen in glycolysis, combined with fermentation In this section, we consider a scenario of intra-pathway cycling with two different co-substrates. 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: a 0 + a 1 = Atot, a 2 + a 3 = Btot, m 1 + m 2 + a 0 = W, m 2 + m 3 + a 2 = M.
Upon substitution of this value of m 3 into dm3 dt = 0 , dm2 dt + dm3 dt = 0 and 2 dm1 dt + dm2 dt + dm3 dt = 0 , and solving recursively for M 2 , m 1 , m 0 , we obtain the following steady state relations: and recall a 1 = Atot − a 0 . We see that if 0 < a 0 < Atot , these expressions are all positive if and only if k in < min(E 1 , E 2 , E 3 /2).
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, Atot) are in one-to-one 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 = Atot . With this information, we cannot immediately conclude that there is a positive root in (0,Atot) . 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 Atot and W . To summarize for positive steady states to exist, we need If k in is small enough, then the conditions hold as α > 0 at k in = 0 . However, these conditions are not sufficient. To see this, by inspecting the term α , 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 Atot = 1, W = 2, F 1 = 1, F 2 = 2, F 3 = 3, L 0 = 1 L 1 = 1, L 2 = 2, K 1 = 1 k in = 1, kout = 1, E 1 = 2, E 2 = 2, E 3 = 3, and note that (35) holds. With K 2 = K 3 = 0.1 , the polynomial of interest becomes p(a 0 ) = 4a 3 0 − 14.300a 2 0 + 7.360a 0 − 0.448 , and so α > 0 . The polynomial has two positive roots under Atot = 1 , namely 0.07 and 0.537. For K 2 = 0.2, K 3 = 0.3 , the polynomial is p(a 0 ) = 4a 3 0 − 23.400a 2 0 + 2.080a 0 − 1.664 , and although α > 0 , it has no root in the interval (0, 1) . An analogous reasoning holds for the irreversible system. From da2 dt + dm1 dt + dm2 dt = 0 we get m 2 = k in kout as expected. From dm2 dt = 0 and da2 dt = 0 we get m 1 = k in K 2 ( F 2 a 0 m 0 + k in a 0 m 0 + L 1 kout ) From the second and third conservation laws, we get This gives that for a positive steady state, we need k in < E 1 , k in < E 2 and max(Atot − W, 0) < a 0 < Atot . Note that Atot − 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 oneto-one 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 Atot − a 0 > 0 and a 0 − Atot − 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 max(Atot − W, 0) < a 0 < Atot . It follows that if (38) has a solution, it has exactly one. Rewriting (38) as a polynomial equation, steady states are in one-to-one correspondence with the roots in the interval max(Atot − W, 0) < a 0 < Atot of the following polynomial p(a 0 ) = ( −L 1 (E 2 − k in )kout + K 2 k in (F 2 + k in ) ) a 2 0 + (L 1 (E 2 − k in )kout(M + T) When a 0 = Atot , the polynomial is positive. Case 1: Atot − W ≤ 0 . In this case we want 0 < a 0 < Atot . If the independent term of p is negative, then p has exactly one solution in (0,Atot) . So, if MT(E 2 − k in ) − K 2 k 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 Atot . Therefore, no positive steady states in this case. Case 2: Atot − W > 0 . In this case, we want Atot − W < a 0 < Atot . We find that p(Atot − W) is If this is negative, then there is one positive steady state, as p is positive at Atot . The condition is Note that M + W − Atot = 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(Atot − W) ≥ 0 , then we are in the situation where p is nonnegative both at Atot − W and Atot , so, if there are roots in the interval (Atot − W, Atot) , 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 in < E 1 and
For all quantities to be positive at steady state, we require 0 < a 0 < Atot . By subtracting from Atot the value of a 0 at steady state, we obtain Atot − a 0 = La ( (A tot + Ka)(k in,1 − k in,2 ) + AtotEa ) .
To summarize, there is a positive steady state if and only if E 1,1 > k in,1 , E 1,2 > k in,2 , and (44) and (43) are positive. For(44) to be positive, and either or AtotFa Atot + La < k in,1 − k in,2 , EaLa + FaKa < (k in,1 − k in,2 )(Ka − La) By analysing these cases, we obtain all possible scenarios for a positive steady state to exist, and these are dictated by the difference k in,1 − k in,2 . For example, if k in,1 = k in,2 , then (45) and (46) < k in,1 − k in,2 .
If Ka < La , the conditions lead to either AtotFa Atot + La > k in,1 − k in,2 > max 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 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 Atot where this applies. An example is shown in Appendix 10-figure 1.