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

  1. Robert West
  2. Hadrien Delattre
  3. Elad Noor
  4. Elisenda Feliu  Is a corresponding author
  5. Orkun S Soyer  Is a corresponding author
  1. School of Life Sciences, University of Warwick, United Kingdom
  2. Department of Plant and Environmental Sciences, Weizmann Institute of Science, Israel
  3. Department of Mathematics, University of Copenhagen, Denmark

Abstract

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.

Editor's evaluation

This manuscript presents an important mathematical analysis of metabolic "co-substrates" and how their cycling can affect metabolic fluxes. Through mathematical analysis of simple network motifs, it shows the impact of co-substrate 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 co-substrate cycling. The work will be of interest to researchers who study microbial metabolism and metabolic engineering.

https://doi.org/10.7554/eLife.84379.sa0

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.

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 (Heinrich et al., 1991). 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 (kappmax) 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 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).

Motif, time-series and threshold in a single co-substrate involving reaction.

(A) Cartoon representation of a single irreversible reaction with co-substrate cycling (see Appendices for other reaction schemes). The co-substrate is considered to have two forms A0 and A1. (B) Concentrations of the metabolites M0 (red) and M1 (green), and the A0/A1 ratio (blue) are shown as a function of time. At t=500, the parameters are switched from the white dot in panel (C) (where a steady state exists) to the black dot (where we see continual build-up of M0 and decline of A0 without steady state). (C & D) Heatmap of the steady state concentration of M0 as a function of the total co-substrate pool size (Atot) and inflow flux (kin). White area shows the region where there is no steady state. On both panels, the dashed line indicates the limitation from the primary enzyme, kin<Vmax,E0, and the solid line indicates the limitation from co-substrate cycling, kin<AtotVmax,Ea/(KM,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 M0 into M1) are picked from within a physiological range (see Supplementary file 1) and are set to: Etot=0.01 mM, kcat=100/s, KM,E0=KM,Ea=50μM, while kout is set to 0.1/s. The Etot and kcat 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 kcat of the co-substrate cycling enzyme, which are set to 0.7 times those for the primary enzyme.

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 M0 to M1 and involves co-substrate cycling (Figure 1A). For co-substrate cycling, we consider additional ‘background’ enzymatic reactions that are independent of M0 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, M0 and M1 (denoted as m0 and m1) are given by:

(1) m0=α kinKM,E0(Vmax,E0kin)(Vmax,EaAtotkin(KM,Ea+Atot))m1=kinkout

where kin and kout denote the rate of in-flux of M0, and out-flux of M1, 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 A0 and A1 in Figure 1A). The parameters Vmax,E0 and Vmax,Ea are the maximal rates (i.e. Vmax=kcatEtot) for the enzymes catalysing the conversion of A0 and M0 into A1 and M1 (enzyme E0), and the turnover of A1 into A0 (enzyme Ea), respectively, while the parameters KM,E0 and KM,Ea are the individual or combined Michaelis-Menten coefficients for these enzymes’ substrates (i.e. for A0 and M0 and A1, respectively). The term α is (in this case where all reactions are irreversible) equal to Vmax,Ea-kin, and in general is a positive expression comprising kin, 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 A0 and A1, 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:

(2) kin<Vmax,E0andkin<AtotVmax,EaKM,Ea+Atot.

Additionally, the ‘shape’ of (1) indicates a ‘threshold effect’ on the steady state value of m0, where it would rise towards infinity as kin 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 M0 is given by kin, while its consumption rate is a function of the Vmax,E0 and the concentration of A0. In turn, the concentration of A0 is determined by its re-generation rate (which is a function of KM,Ea and Vmax,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 (E0 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 kin<Vmax,E0, while considering it as a non-enzymatic reaction with co-substrate cycling only, leads to only the constraint kin<AtotVmax,Ea/(KM,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.

Motif, time-series and thresholds for the linear pathway model with n=4.

(A) Cartoon representation of a chain of reversible reactions with co-substrate cycling occurring solely inter-pathway. The co-substrate is considered to have two forms A0 and A1. (B) Heatmap of the steady state concentration of M0 as a function of the total metabolite pool size (Atot) and inflow rate constant (kin). White area shows the region where there is no steady state. The dashed and solid lines indicate the limitations arising from primary enzyme (E1 in this case) and co-substrate cycling, respectively, as in Figure 1. (C) Concentrations of M0-4, and A0/A1 ratio as a function of time (with colours as indicated in the inset). At t=1000 s, the parameters are switched from the white dot in panel (B) (where a steady state exists) to the black dot (where we see build-up of all substrates that are produced before the first co-substrate cycling reaction, and continued decline of A0). The parameters used are picked from within a physiological range (see Supplementary file 1) and are set to: Etot=0.01 mM, kcat=100/s, KM=50μm, for all reactions, and kout=0.1/s.

We find the same kind of threshold dynamics as in the single reaction case. When kin is above a threshold value, the metabolite M0 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 kin (Figure 2C and Appendix 3—figure 1). In other words, pathway output cannot be increased further by increasing kin beyond the threshold. Finally, as kin 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 kin, 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 kin 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.

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 kcat 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 kin (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-kcat 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.

Measured and FBA-predicted flux values are typically lower than the calculated primary enzyme threshold.

(A) Measured and FBA-predicted flux values (from 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.

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 flux – v/(E0kcat) (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 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.

Motif, heatmap and time-series for the branching pathway model.

(A) Cartoon representation of two branching pathways from the same upstream metabolite. The two branches are linked to separate co-substrate pools, A and B. Note that pathway independent turnover of the co-substrates is included in the model (see Figure 4—source code 1). (B) The pathways’ flux ratio (i.e. flux into M2,2 divided by flux into M2,1) shown in colour mapping, against the ratio of co-substrate pool sizes, Atot and Btot, and the influx rate, kin, into the upstream metabolite. In the block colour areas, the system has no steady state and the indicated metabolite(s) M0 and one of the metabolites M1,2 or M1,1 accumulate towards infinity. (C) Concentrations of upstream and branch-endpoint metabolites over time, coloured as shown in the inset of the panel. The solid lines show results using parameters indicated by the white dot in panel (B), where Btot>Atot, while the dashed lines show results using parameters indicated by the black dot in panel (B), where Atot>Btot. For both simulations, all kinetic parameters are arbitrarily set to 1, apart from the pathway-independent co-substrate recycling (Vmax,Ea and Vmax,Eb) that is set to 10 (see Figure 4—source code 1).

Figure 4—source code 1

Python implementation of branched pathway model, presented in Figure 4.

https://cdn.elifesciences.org/articles/84379/elife-84379-fig4-code1-v2.zip

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 kin values, but the threshold effect is still present with high enough kin 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 pool size (see upper corner regions on 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 kin. 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 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.

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 A0 and A1. It is converted from A0 to A1 on the lower pathway, and from A1 to A0 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, M1,2 and M1,1, as a function of the total amount of co-substrate (Atot) and the extent of fluctuations in the two pathway influxes, kin,1 and kin,2. The influx fluctuation is characterised by a waiting time that is exponentially distributed with mean τ, after which the log ratio of the kin values is drawn from a standard normal distribution. The mean of the kin 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 M1,2 (green) and M1,1 (magenta), pathway influx ratio (pink), and A0/A1 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 kin fluctuations is on a similar timescale to the other reactions. The system is largely unresponsive to the noise, and the M1,2 and M1,1 are fully correlated (i.e. the green and magenta curves overlap). (D) Concentrations of metabolites M1,2 (green) and M1,1 (magenta), pathway influx ratio (pink), and A0/A1 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 kin values occur at a much lower rate than the other reactions. The system is responsive to the noise, yet the M1,2 and M1,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 pathway-independent co-substrate recycling (Vmax,Ea) that is set to 0.01 (see Figure 5—source code 1).

Figure 5—source code 1

Python implementation of connected pathway model, presented in Figure 5.

https://cdn.elifesciences.org/articles/84379/elife-84379-fig5-code1-v2.zip

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 Vmax,Ea=0.01 are representative of those for Vmax,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, kin, 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 nitrogen-source concentrations) causing higher kin values across various pathways, their metabolic systems can be close to the threshold presented here. While both kin 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

Request a detailed protocol

The metabolic system shown in Figure 1A comprises the following biochemical reactions:

(3) kinM0M1koutM0+A0M1+A1A0A1

where metabolites are denoted by Mi and the different forms of the co-substrate are denoted by Ai. We assume additional conversion between A1 and A0, mediated through other enzymatic reactions. The parameters kin, and kout denote the in- and out- flux of M0 and M1 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:

(4) dm0dt=kinVmax,E0 a0m0KM,E0+a0m0dm1dt=Vmax,E0 a0m0KM,E0+a0m0koutm1da0dt=Vmax,Ea a1KM,Ea+a1Vmax,E0 a0m0KM,E0+a0m0da1dt=Vmax,E0 a0m0KM,E0+a0m0Vmax,Ea a1KM,Ea+a1

where mi and ai denote the concentrations of Mi and Ai respectively, KM 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 a0+a1=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.

Symbolic and numerical computations

Request a detailed protocol

For 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 protocol

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).

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 A0, A1. 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, A0 and A1, 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 Mi-1Mi, while in the latter case, they are shown as Mi-1Mi. 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:

A0k6k5A1.

We talk about irreversible co-substrate conversion, if k5=0 or k6=0, that is, only conversion in one direction is considered. We talk about no free co-substrate conversion, if k5=k6=0, that is, the co-substrate cycling is only related through the reactions in the metabolic pathway.

Appendix 1—figure 1
Cartoon representation of select metabolic pathways involving co-substrate cycling.

(A) Cartoon representation of upper glycolysis pathway. Note that stoichiometric balance across the pathway changes after F1,6P and metabolites and co-substrates highlighted with gray background have a stoichiometry of 2. (B & C) Cartoon representation of nitrogen assimilation via glutamate and involvement of glutamate cycling in amino acid biosynthesis. (D) Cartoon representation of metabolite cycles involved in one-carbon metabolism around tetrahydorfolate. Enzyme and metabolite name abbreviations are: Glc – Glucose, G6P – Glucose-6-phosphate, F6P – Fructose-6-phosphate, F1,6P – Fructose-1,6-biphosphate, G3P – glyceraldehyde-3-phosphate, G1,3P – 1,3-biphospho-D-glycerate, 3 PG – 3-phospho-D-glycerate, PEP – Phosphoenolpyruvate, Pyr – Pyruvate, Lac – Lactate, DHF – Dihydrofolate, THF - Tetrahydrofolate, AICAR - 5-amino-4-imidazolecarboxamide ribotide, FAICAR - 5’-phosphoribosyl-formamido-carboxamide, hk – hexokinase, pfk – phosphofructokinase, gapdh – glyceraldehyde-3-phosphate dehydrogenase, pgk – phosphoglycerate kinase, pyk – phophoenolpyruvate kinase, ldh – lactate dehydrogenase, glud – glutamate dehydrogenase, GOGAT – glutamate synthase, DHFR – dihydrofolate reducatase.

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 a0,a1 denote the concentrations of the co-substrate pair A0 and A1, respectively and mi to denote the concentration of Mi, 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:

Mi1+A0+EkoffkonEMi1A0koff,1kon,1EMiA1koff,2kon,2Mi+A1+E.

For this reversible reaction scheme, the rate of production of Mi takes the form

v=EiLi1mi1a0FiKimia1KiLi1+Kimia1+Li1mi1a0.

Likewise, for the reversible enzymatic conversion Mi-1Mi, we have the following reaction scheme:

Mi1+EkoffkonEMi1koff,1kon,1EMikoff,2kon,2Mi+E.

The rate of production of Mi is given by

v=EiLi1mi1FiKimiKiLi1+Kimi+Li1mi1.

In both of these reversible rate equations, the parameters K and L are equivalent to the Haldane coefficients KS and KP, respectively and are given by

(5) Ki=kon,1kon,2+koffkon,2+koffkoff,1kon(kon,2+koff,1+kon,1)andLi1=k on,1kon,2+koffkon,2+koffkoff,1koff,2(koff,1+k on,1+koff).

When there are two substrates that take part in the reaction, the kon and koff,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 KS and KP.

The parameters E and F correspond to the Haldane coefficients kcat+ and kcat-, multiplied by the total enzyme concentration (denoted Etot, below), and are given by

(6) Ei=Etotkon,1kon,2kon,2+koff,1+kon,1andFi=Etotkoffkoff,1koff+koff,1+kon,1.

For the irreversible enzymatic reaction, the reaction schemes simplify to:

Mi1+A0+EkoffkonEMi1A0kcatMi+A1+E.
Mi1+EkoffkonEMi1kcatMi+E.

And the rate of production for the two cases are given by

v=Eimi1a0Ki+mi1a0,andv=Eimi1Ki+mi1,

where Ei is the catalytic rate coefficient of the i-th enzyme multiplied by its total concentration, and Ki is its Michaelis-Menten coefficient. Again, when there are two substrates, the kon 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 Ki 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 ++-, +--, -+-, + 0-, 0+-.

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 A0 and A1, with biosynthesis and degradation of both forms.

0k2k1A00k4k3A1A0k6k5A1.

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 a0=[A0] and a1=[A1], the differential equations for these concentrations are:

da0dt=k1(k2+k5)a0+k6a1da1dt=k3(k4+k6)a1+k5a0.

Since all the terms are linear or constant, the steady state values are the solutions of the linear equation:

((k2+k5)k6k5(k4+k6))(a0a1)=(k1k3).

The steady states are then found to be:

a0=k1k4+k1k6+k3k6k4k5+k2k4+k2k6,a1=k2k3+k3k5+k1k5k4k5+k2k4+k2k6.

If we consider the case where the synthesis and degradation rates of the different forms of the co-substrate (i.e. cycled metabolite) are the same, i.e. k1 = k3=ks and k2 = k4=kd, these equations simplify to:

a0=ks(kd+2k6)kd(kd+k5+k6),a1=ks(kd+2k5)kd(kd+k5+k6),

and the eigenvalues of the Jacobian of the system evaluated at this steady state are always real and negative. When kd 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:

a0=2k6ks/kdk5+k6,a1=2k5ks/kdk5+k6,

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. k1=k2=k3=k4=0). In that case, the steady states are a0=Tk6/(k5+k6) and a1=Tk5/(k5+k6), 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 co-substrates 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 k1=k3=ks and k2=k4=kd, the ODEs for the system are:

da0dt=kskda0a0EaLaa1FaKaa0La+a1Ka+KaLa,da1dt=kskda1+a0EaLaa1FaKaa0La+a1Ka+KaLa.

The only real and positive steady state is now found to be:

a0=(KaLakd+Ka(Fa+3ks)+La(Eaks)4Kaks(KaLa)(2(Fa+ks)+kdLa)+(FaKa+3Kaks+EaLa+KakdLaksLa)2)/(2kd(KaLa)),a1=(KaLakdKa(Fa3ks)La(Ea+3ks)+4Kaks(KaLa)(2(Fa+ks)+kdLa)+(FaKa+3Kaks+EaLa+KakdLaksLa)2)/(2kd(KaLa)).

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 a0+a1 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:

0kinM0M0+A0M1+A1M1kout0A0A1.

This corresponds to the motif depicted in Figure 1A of the main text, and the resulting ODEs are:

dm0dt=kinE1L0m0a0F1K1m1a1K1m1a1+L0a0m0+K1L0dm1dt=E1L0m0a0F1K1m1a1K1m1a1+L0a0m0+K1L0koutm1da0dt=E1L0m0a0F1K1m1a1K1m1a1+L0a0m0+K1L0EaLaa0FaKaa1Kaa1+Laa0+KaLada1dt=E1L0m0a0F1K1m1a1K1m1a1+L0a0m0+K1L0+EaLaa0FaKaa1Kaa1+Laa0+KaLa.

This ODE system has one conservation law, namely the sum of a0 and a1 is constant:

a0+a1=Atot.

The steady states of the system are:

m0=K1kin(E1kin)(FaAtotkin(La+Atot))α,m1=kinkout,a0=Ka(FaAtotkin(La+Atot))Ka(Fakin)+La(Ea+kin),a1=La(Kakin+T(Ea+kin))Ka(Fakin)+La(Ea+kin).

where, we introduced the composite parameter α, as follows:

(7) α=koutL0(Ka(Fakin)+La(Ea+k))+(F1+kin)La(Kakin+kinAtot+EaAtot)KaL0kout.

For the steady state equations given above to be positive, the following conditions must be satisfied:

(8) kin<E1andkin<FaAtotLa+Atot.

Note that as FaAtotLa+Atot<Fa, the second condition readily implies kin<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:

dm0dt=kinE1m0a0K1+m0a0,dm1dt=E1m0a0K1+m0a0koutm1,da0dt=E1m0a0K1+m0a0EaLaa0FaKaa1Kaa1+Laa0+KaLa,da1dt=E1m0a0K1+m0a0+EaLaa0FaKaa1Kaa1+Laa0+KaLa.

The steady state of the system is:

m0=K1kin(E1kin)(FaAtotkin(La+Atot))α,m1=kinkout,a0=Ka(FaAtotkin(La+Atot))Ka(Fakin)+La(Ea+kin),a1=La(Kakin+Atot(Ea+kin))Ka(Fakin)+La(Ea+kin).

where the composite parameter α is defined (differently to the reversible case) as:

(9) α=Ka(Fakin)+La(Ea+kin)Ka.

For this to be positive the same conditions as in the reversible case must be satisfied:

(10) kin<E1andkin<FaAtotLa+Atot.

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 A1 to A0, the ODEs describing the system dynamics are:

dm0dt=kinE1m0a0K1+m0a0,dm1dt=E1m0a0K1+m0a0koutm1,da0dt=E1m0a0K1+m0a0+Faa1La+a1,da1dt=E1m0a0K1+m0a0Faa1La+a1.

The steady state of the system is:

m0=K1kin(E1kin)(FaAtotkin(La+Atot))α,m1=kinkout,a0=AtotkinLaFakin,a1=kinLaFakin.

where the composite parameter α is defined as:

(11) α=Fakin

For these steady states to be positive the same conditions as in the other two cases must be satisfied:

(12) kin<E1andkin<FaAtotLa+Atot.

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:

0kinM0M0+A0M1+A1M1kout0A0k6k5A10k8k7A00k10k9A1.

The resulting system of ODEs is:

dm0dt=kinE1L0m0a0F1K1m1a1K1m1a1+L0a0m0+K1L0dm1dt=E1L0m0a0F1K1m1a1K1m1a1+L0a0m0+K1L0koutm1da0dt=E1L0m0a0F1K1m1a1K1m1a1+L0a0m0+K1L0(k5+k8)a0+k6a1+k7da1dt=E1L0m0a0F1K1m1a1K1m1a1+L0a0m0+K1L0+k5a0(k6+k10)a1+k9.

The steady state concentrations are then given by:

a0=k10(k7kin)+k6(k7+k9)k10k5+k10k8+k6k8a1=k8(k9+kin)+k5(k7+k9)k10k5+k10k8+k6k8m0=K1kin((F1+kin)(k8(k9+kin)+k5(k7+k9))+koutL0(k10k5+k10k8+k6k8))L0kout(E1kin)(k10(k7kin)+k6(k7+k9))m1=kinkout.

These expressions are positive if and only if

kin<k7+k6k10(k7+k9)andkin<E1.

If the main path is irreversible, similarly to the previous section, the ODEs describing the system’s dynamics, are:

dm0dt=kinE1m0a0K1+m0a0dm1dt=E1m0a0K1+m0a0koutm1da0dt=E1m0a0K1+m0a0(k5+k8)a0+k6a1+k7da1dt=E1m0a0K1+m0a0+k5a0(k6+k10)a1+k9.

The steady state in this case is:

a0=k10(k7kin)+k6(k7+k9)k10k5+k10k8+k8k6a1=k8(k9+kin)+k5(k7+k9)k10k5+k10k8+k8k6m0=K1kin(k10k5+k10k8+k8k6)(E1kin)(k10(k7kin)+k6(k7+k9))m1=kinkout.

These expressions are positive if and only if

kin<k7+k6k10(k7+k9)andkin<E1.

Comparing these results with those of Subsection C.1, we see some similar themes. Firstly, the total amount of co-substrate a0+a1 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 kin<E1, and any other conditions take the form kin<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:

0kinM0M0k2k1M1M1kout0.

The resulting system of ODEs is:

dm0dt=kinE1L0m0F1K1m1K1m1+L0m0+K1L0dm1dt=E1L0m0F1K1m1K1m1+L0m0+K1L0koutm1.

In this case, there is no conservation law, and the steady states of the system are:

m0=K1kin(F1+kin+koutL0)koutL0(E1kin),m1=kinkout.

These expressions are positive if and only if kin<E1 and the steady state is asymptotically stable when this holds.

If the main path is irreversible, the ODEs describing system dynamics are:

dm0dt=kinE1m0K1+m0dm1dt=E1m0K1+m0koutm1,

and the steady state is:

m0=K1kinE1kin,m1=kinkout.

Again, expressions are positive if and only if kin<E1 and the steady state is asymptotically stable when this holds.

From this we see that the condition kin<E1 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:

0kinM0M0+A0k2k1M1A1M1kout0A0k4k3A1.

The resulting system of ODEs, describing system dynamics, is:

dm0dt=kink1m0a0+k2m1a1dm1dt=k1m0a0k2m1a1koutm1da0dt=k1m0a0+k2a1m1k3a0+k4a1da1dt=k1m0a0k2m1a1+k3a0k4a1

There is the conservation law

a0+a1=Atot.

The steady state of the system is:

m0=kin(k2(kin+k3Atot)+kout(k3+k4))k1kout(k4Atotkin)m1=kinkout,a0=k4Atotkink3+k4,a1=k3Atot+kink3+k4.

These expressions are positive if and only if kin<k4Atot, and when this is satisfied the steady state is asymptotically stable.

If the system reactions are irreversible, the ODEs describing the system dynamics are:

dm0dt=kink1m0a0dm1dt=k1m0a0koutm1da0dt=k1m0a0+k4a1da1dt=k1m0a0k4a1

and the steady state is:

m0=k4kink1(k4Atotkin),m1=kinkout,a0=k4Atotkink4,a1=kink4.

As in the reversible case, these expressions are positive when kin<k4Atot 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.

Appendix 3—figure 1
Concentrations of non-accumulating metabolites vs kin in the models presented in Figures 1 and 2 of the main text.

(A) Concentration of non-accumulating metabolite M1 from Figure 1. Red line shows the critical kin value, after which the concentration of M1 remains constant (B) Concentrations of non-accumulating metabolites from the models presented in Figure 2. Blue line shows the threshold value for kin. Once the threshold value for kin is reached the concentrations of the non-accumulating metabolites do not change. Atot=10-1 in panel (A) and 10-6 in panel (B). All other parameters are the same as their counterparts in the main text.

Appendix 4

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:

(13) 0kinM0M0+A0M1+A1Mn1+A1Mn+A0M1M2Mi1MiMi+1Mn1Mnkout0.

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:

(14) 0kinM0M0+A0M1+A1M1M2M3kout0M2+A1M3+A0.

Considering reversible reaction kinetics, the evolution of the concentrations of the species in time is modelled by the following ODE system:

dm0dt=kinE1L0m0a0F1K1m1a1K1m1a1+L0m0a0+K1L0dm1dt=E1L0m0a0F1K1m1a1K1m1a1+L0m0a0+K1L0E2L1m1F2K2m2K2m2+L1m1+K2L1dm2dt=E2L1m1F2K2m2K2m2+L1m1+K2L1E3L2m2a1F3K3m3a0K3m3a0+L2m2a1+K3L2dm3dt=E3L2m2a1F3K3m3a0K3m3a0+L2m2a1+K3L2koutm3da0dt=E1L0m0a0F1K1m1a1K1m1a1+L0m0a0+K1L0+E3L2m2a1F3K3m3a0K3m3a0+L2m2a1+K3L2da1dt=E1L0m0a0F1K1m1a1K1m1a1+L0m0a0+K1L0E3L2m2a1F3K3m3a0K3m3a0+L2m2a1+K3L2.

The system has two conservation laws:

(15) a0+a1=Atot,m1+m2+a0=W.

Solving recursively the steady state equations given by dm0dt++dm3dt=0, dm3dt=0, dm2dt+dm3dt=0, dm1dt+dm2dt+dm3dt=0 for m3,m2,m1,m0 and the first conservation law, we obtain

(16) m0=K1((F1+kin)m1a1+L0kin)a0L0(E1kin),m1=K2((F2+kin)m2+L1kin)L1(E2kin),m2=kinK3((F3+kin)a0+L2kout)kouta1L2(E3kin),m3=kinkout,a1=Atota0.

By substituting recursively m2 in m1 and m1 in m0, we obtain expressions of all variables at steady state in terms of a0. For these to be positive, it needs to hold that 0<a0<Atot, and further

(17) kin<min(E1,E2,E3).

We substitute (16) into the remaining conservation law, and obtain a polynomial in a0 whose roots in the interval (0,T) are in one-to-one correspondence with the positive steady states, provided (17) holds. The polynomial is

p(a0)=L1L2kout(E3kin)(E2kin)a02+(L1L2kout(Atot+W)(E3kin)(E2kin)K2L1L2kinkout(E3kin)+K3kin(F3+kin)(E2L1,+F2K2)+K3kin2(K2L1)(F3+kin,))a0+L2K2α,

where

(18) α=kin(kin+F2)K2K3+kinAtot(E3kin)K2L1+kin(E2kin)K3L1AtotW(E3kin)(E2kin)L1.

As a0=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 kin=0, α=AtotW(E3kin)(E2kin)L1<0. Hence, for kin small enough, α<0 and also (17) holds, and the system has a positive steady state. The steady state value of a0 is found by solving the quadratic equation p(a0)=0 and considering the smallest root.

We note that α is a polynomial in kin 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 kin=E2 or kin=E3, then at this value of kin, α 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

kin<min(E1,B).

If the minimum is attained at B, then, as kin approaches B, the smallest positive root of p(a0) approaches 0 as α approaches zero. Using this in (16), we obtain that

(19) m0kinB ,m1kinB Bkout(F2K3+L1Atot(E3B)+BK3)L1Atot(E3B)(E2B),m2kinBBK3Atot(E3B)m3kinB Bkout.

If E1<B, then kin needs to be smaller than E1. When kin approaches E1, a0 approaches the root of p when kin=E1, which is not zero. Then m0 still approaches as it has E1-kin in the denominator, while m1,m2 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.

(20) 0kinM0M0M1M1MM+A0M+1+A1M+1Mρ1MρMρ+A1Mρ+1+A0Mρ+1Mn1MnMnkout0.

For the system to make sense biochemically, we require 0<ρ1n1.

To write the equations in a simple format, we write

xi=mi,i,ρ,x=ma0,xρ=mρa1

and

yi=mi,i+1,ρ+1,y+1=m+1a1,yρ+1=mρ+1a0.

We start with the reversible reaction kinetics. Then the associated ODE system becomes:

dm0dt=kinE1L0x0F1K1y1K1L0+L0x0+K1y1dmidt=EiLi1xi1FiKiyiKiLi1+Kiyi+Li1xi1Ei+1LixiFi+1Ki+1yi+1Ki+1Li+Ki+1yi+1+Lixii=1,,n1dmndt=EnLn1xn1FnKnynKnLn1+Knyn+Ln1xn1koutmnda0dt=E+1LxF+1K+1y+1K+1L+Lx+K+1y+1+Eρ+1LρxρFρ+1Kρ+1yρ+1Kρ+1Lρ+Lρxρ+Kρ+1yρ+1da1dt=E+1LxF+1K+1y+1K+1L+Lx+K+1y+1Eρ+1LρxρFρ+1Kρ+1yρ+1Kρ+1Lρ+Lρxρ+Kρ+1yρ+1.

This ODE system has precisely two conservation laws:

a0+a1=Atot,m+1++mρ+a0=W.

Note that we have the following equalities:

(21) dm0dt++dmidt=kinEi+1LixiFi+1Ki+1yi+1Ki+1Li+Ki+1yi+1+Lixi,i=0,,n1
(22) dm0dt++dmndt=kinkoutmn

For i=1,,n-2, by solving (21) for xi, we obtain the following recursive formulas:

(23) xi=Ki+1(kinLi+(kin+Fi+1)yi+1)(Ei+1kin)Li,i=0,,n1.

Finally, from (22) and the conservation law a0+a1=Atot, we obtain

(24) xn=kinkout,a1=Atota0.

These expressions are all positive if and only if 0<a0<Atot and kin<Ei for i=1,,n. Note that the value of mi can now be found recursively from mn using (23), as we show next.

Recall that 0<ρ1n1. Then, it holds

yi=xiϵi,ϵi=1,i,+1,ρ,ρ+1,ϵ=ϵρ+11=1a0,ϵ+1=ϵρ1=a1.

We write for shortness

xi=zi+bixi+1ci,i=0,,n1,

where for i=0,,n-1,

(25) zi=kinKi+1Libi=ϵi+1(kin+Fi+1)Ki+1ci=(Ei+1kin)Li.

Then we claim that

(26) xi=(j=in1zj(bibj1)(cj+1cn1))+(bibn1)xncicn1,i=1,,n1.

We prove this by descending induction on i. Note that products over an empty index equal 1. For instance, bjbj1. 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 xi. For n1>i10, we have

xi1=zi1+bi1xici1=zi1+bi1(j=in1zj(bibj1)(cj+1cn1))+(bibn1)xncicn1ci1=zi1cicn1+bi1((j=in1zj(bibj1)(cj+1cn1))+(bibn1)xn)ci1cicn1=zi1cicn1+(j=in1zj(bi1bj1)(cj+1cn1))+(bi1bn1)xnci1cn1=(j=i1n1zj(bi1bj1)(cj+1cn1))+(bi1bn1)xnci1cn1.

This is (26) for i-1. Hence, the formula holds. Finally, we can use that xn=mn=kinkout to obtain:

(27) xi=kout(j=in1zj(bibj1)(cj+1cn1))+(bibn1)kinkout cicn1,i=0,,n1.

Let

b¯u,j:=i=uj(kin+Fi+1)Ki+1,c¯u,j:=cucj=i=uj(Ei+1kin)Li.

This gives:

(28) mi=kout(j=in1zj(bibj1)c¯j+1,n1)+(bibn1)kinkout c¯i,n1,i,ρ,
(29) m=kout(j=in1zj(bibj1)c¯j+1,n1)+(bibn1)kina0kout c¯,n1,
(30) mρ=kout(j=in1zj(bibj1)c¯j+1,n1)+(bibn1)kin(Ta0)kout c¯ρ,n1

Remember that b,b+1,bρ,bρ+1 depend on a0,a1, while the rest of b’s do not. For a product of the form bubj with uj, we have the following:

If u-1:

 bubj={i=uj(kin+Fi+1)Ki+1j<1, or ρj1a0i=uj(kin+Fi+1)Ki+1j=1, or j=ρ1a1a0i=uj(kin+Fi+1)Ki+11<j<ρ1.

If u=:

bubj={a1i=uj(kin+Fi+1)Ki+1j<ρ1i=uj(kin+Fi+1)Ki+1j=ρ1a0i=uj(kin+Fi+1)Ki+1ρj.

If <uρ1:

bubj={i=uj(kin+Fi+1)Ki+1j<ρ11a1i=uj(kin+Fi+1)Ki+1j=ρ1a0a1i=uj(kin+Fi+1)Ki+1ρj.

If u=ρ:

 bubj=a0i=uj(kin+Fi+1)Ki+1.

If u>ρ:

 bubj=i=uj(kin+Fi+1)Ki+1.

Summarising bubj, equals

{i=uj(kin+Fi+1)Ki+1u1,j<1, or u1,ρj, or u=,j=ρ1, or <uρ1,j<ρ1, or u>ρ1a0i=uj(kin+Fi+1)Ki+1j=1, or u1,j=ρ1a1a0i=uj(kin+Fi+1)Ki+1u1,1<j<ρ1a1i=uj(kin+Fi+1)Ki+1u=,j<ρ1a0i=uj(kin+Fi+1)Ki+1u=,ρj, or u=ρ1a1i=uj(kin+Fi+1)Ki+1<uρ1,j=ρ1a0a1i=uj(kin+Fi+1)Ki+1<uρ1,ρj.

Observe that for all i, the summand of j=in-1zj(bibj-1)(cj+1cn-1) corresponding to i=j is

 zi(ci+1cn1),

which does not depend either on a0 or a1 .

In particular, we deduce that mi for i-1, the term (bibn-1) does not depend on a0,a1, and in the sum j=in-1zj(bibj-1)(cj+1cn-1) there are summands involving 1a0, for eample when j=. Hence mi, for i-1, goes to infinity as a0 goes to zero. Note that a1 goes to T when a0 goes to zero. When i=, the denominator of mi is multiplied by a0. As the numerator has at least one term that is not multiple of a0, m goes to infinity as well. We conclude that

 mi a00 ,0i.

When i, then no summand in the numerator of mi in (27–29) that involves 1a0, and neither the denominator has a0 as factor. As the numerator has at least one term that is not multiple of a0, mi goes to a finite value as a0 goes to zero.

 mi a00 number,i.

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+1++mρ+a0=W, to determine the value of a0. We have

i=+1ρmi+a0W=i=+1ρ1kout(j=in1zj(bibj1)c¯j+1,n1)+(bibn1)kinkout c¯i,n1+kout(j=ρn1zj(bρbj1)c¯j+1,n1)+(bρbn1)kina1kout c¯ρ,n1+a0W=().

We have (bibn-1)=a0a1b¯i,n-1 when <i<ρ and (bρbn-1)=a0b¯ρ,n-1. Also

j=ρn1zj(bρbj1)c¯j+1,n1=zρc¯ρ+1,n1+j=ρ+1n1a0zjb¯ρ,j1c¯j+1,n1.

Finally, for +1iρ-1,

j=in1zj(bibj1)c¯j+1,n1=j=iρ1zj(bibj1)c¯j+1,n1+zρ(bibρ1)c¯ρ+1,n1+j=ρ+1n1zj(bibj1)c¯j+1,n1=j=iρ1zjb¯i,j1c¯j+1,n1+zρ1a1b¯i,ρ1c¯ρ+1,n1+a0a1j=ρ+1n1zjb¯i,j1c¯j+1,n1.

This gives

()=i=+1ρ1kout(j=iρ1zjb¯i,j1c¯j+1,n1+zρ1a1b¯i,ρ1c¯ρ+1,n1+a0a1j=ρ+1n1zjb¯i,j1c¯j+1,n1)+a0a1b¯i,n1kinkout c¯i,n1+kout(zρc¯ρ+1,n1+j=ρ+1n1a0zjb¯ρ,j1c¯j+1,n1)+a0b¯ρ,n1kina1kout c¯ρ,n1+a0W=i=+1ρ1kout(a1j=iρ1zjb¯i,j1c¯j+1,n1+zρb¯i,ρ1c¯ρ+1,n1+a0j=ρ+1n1zjb¯i,j1c¯j+1,n1)+a0b¯i,n1kina1kout c¯i,n1+kout(zρc¯ρ+1,n1+j=ρ+1n1a0zjb¯ρ,j1c¯j+1,n1)+a0b¯ρ,n1kina1kout c¯ρ,n1+a0W.

Hence (*) vanishes if and only if

(31) 0=i=+1ρ1kout(0,0,1a1j=iρ1zjb¯i,j1c¯+1,i1c¯j+1,n1+zρb¯i,ρ1c¯+1,i1c¯ρ+1,n1+0,0,1a0j=ρ+1n1zjb¯i,j1c¯+1,i1c¯j+1,n1)+0,0,1a0b¯i,n1c¯+1,i1kin+kout(zρc¯+1,ρ1c¯ρ+1,n1+j=ρ+1n10,0,1a0zjb¯ρ,j1c¯+1,ρ1c¯j+1,n1)+0,0,1a0b¯ρ,n1c¯+1,ρ1kin+(0,0,1a0W)koutc¯+1,n10,0,1a1.

As a1=Atot-a0, this is a degree 2 polynomial in a0. The leading term comes from the term (0,0,1a0W)koutc¯+1,n10,0,1a1, and is

koutc¯+1,n1a02,

which is negative. The independent term is

i=+1ρ1kout(Atotj=iρ1zjb¯i,j1c¯+1,i1c¯j+1,n1+zρb¯i,ρ1c¯+1,i1c¯ρ+1,n1)+koutzρc¯+1,ρ1c¯ρ+1,n1AtotWkoutc¯+1,n1.

We divide by kout and define

α=i=+1ρ1c¯+1,i1(Atotj=iρ1zjb¯i,j1c¯j+1,n1+zρb¯i,ρ1c¯ρ+1,n1)+zρc¯+1,ρ1c¯ρ+1,n1AtotWc¯+1,n1,

where recall from (25) that

 zi=kinKi+1Libi=ϵi+1(kin+Fi+1)Ki+1ci=(Ei+1kin)Li.

When a0=Atot, all terms multiplying a1 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 kin=0 this inequality holds, as all terms with zj vanish. When kin approaches one of Ei with +1<in, 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

α=AtotW(E2kin)L1(E3kin)L2(E4kin)+kinK4(E2kin)L1(E3kin)L2+AtotkinkoutL1(E4kin)(E3kin)L2+AtotkinK3L2(E4kin)(F2+kin)kout+kinK4(F3+kin)K3(F2+kin)kout+kinK3L2(E2kin)L1(E4kin)Atot+kinK4(E2kin)L1(F3+kin)K3.

Let B be the smallest positive root of α=0 as a polynomial in kin, 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 Ej is attained for some Ei with +1<in, then α is positive at kin=Ei. In that case, as α is negative at kin=0, there is at least one value of kin<Ei for all i such that α=0 and hence B is finite. This shows that min(E1,,En,B)=min(E1,,E+1,B).

Putting it all together, we have shown that for all

(32)  kin<min(E1,,E+1,B)

the system has positive steady states, and if

kinmin(E1,,E+1,B),orkin<min(E1,,E+1),B<kin<B,

Taking condition (32), if the minimum is attained at B, then when kin 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 a0, we see that a0 approaches zero, and the mi converge to the values given above. Specifically,

mj kinB ,for 0j,mj kinB number,for j.

By the comment above, (32) cannot be attained at Ei with +1<in. If (32) is attained at Ei with i+1, then as kin approaches this minimum, a0 converges to a number. In this case, the concentrations that tend to infinity are those with Ei-kin in the denominator:

mj kinEi ,for 0j<i,mj kinEi number,for ji.

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:

(33) 0kinM0M0+ATPM1+ADPM1+NADM2+NADHM4kout0M2+ADPM3+ATPM3+NADHM4+NAD

We write A0= ATP, A1= ADP, A2= NAD, A3= NADH, for simplicity. The ODE system governing the dynamics of the network is:

dm0dt=kinE1L0a0m0F1K1a1m1K1a1m1+L0a0m0+K1L0dm1dt=E1L0a0m0F1K1a1m1K1a1m1+L0a0m0+K1L0E2L1a2m1F2K2a3m2K2a3m2+L1a2m1+K2L1dm2dt=E2L1a2m1F2K2a3m2K2a3m2+L1a2m1+K2L1E3L2a1m2F3K3a0m3K3a0m3+L2a1m2+K3L2dm3dt=E3L2a1m2F3K3a0m3K3a0m3+L2a1m2+K3L2E4L3a3m3F4K4a2m4K4a2m4+L3a3m3+K4L3dm4dt=E4L3a3m3F4K4a2m4K4a2m4+L3a3m3+K4L3koutm4da0dt=E1L0a0m0F1K1a1m1K1a1m1+L0a0m0+K1L0E3L2a1m2F3K3a0m3K3a0m3+L2a1m2+K3L2da1dt=E1L0a0m0F1K1a1m1K1a1m1+L0a0m0+K1L0E3L2a1m2F3K3a0m3K3a0m3+L2a1m2+K3L2da2dt=E2L1a2m1F2K2a3m2K2a3m2+L1a2m1+K2L1+E4L3a3m3F4K4a2m4K4a2m4+L3a3m3+K4L3da3dt=E2L1a2m1F2K2a3m2K2a3m2+L1a2m1+K2L1E4L3a3m3F4K4a2m4K4a2m4+L3a3m3+K4L3.

There are four conservation laws:

 a0+a1=Atot,a2+a3=Btot,m1+m2+a0=W,m2+m3+a2=M.

We consider the equations dm0dt++dm4dt=0, dm4dt=0, dm3dt+dm4dt=0, dm2dt+dm3dt+dm4dt=0, dm1dt+dm2dt+dm3dt+dm4dt=0, and solve them iteratively for m4,m3,m2,m1,m0 and obtain:

m4=kinkout,m3=kinK4((F4+kin)a2+L3kout)kouta3L3(E4kin),m2=K3((F3+kin)a0m3+L2kin)a1L2(E3kin),m1=K2((F2+kin)a3m2+L1kin)a2L1(E2kin),m0=K1((F1+kin)a1m1+L0kin)a0L0(E1kin),

and a1=Atot-a0, a3=Btot-a2. As usual, a necessary condition for positive steady states is

(34) kin<min(E1,E2,E3,E4).

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 m1+m2+a0=W, and m2+m3+a2=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

M=1,W=1,F1=1F2=1,F3=1,F4=2,K1=1K2=0.1,K3=0.1,K4=2,L0=1L1=1L2=2,L3=10,Atot=3,Btot=2,E1=3,E2=3,E3=12,E4=12,kin=2,kout=10,

the polynomial system becomes

0=20000a22a060000a22+180560a260315.2a2a095200+32120a00=20000a22a0240000a2a0279996.64a22a0+164086.88a2a07928a0+59720a22131680a2+24480

and the solutions are:

a00.69,a20.679,

which are in the desired interval.

If instead we replace L3 by 1, we obtain a system whose solutions are

a02.978,a20.323,

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 5—figure 1
Motif, time series and threshold for the linear pathway model with multiple co-substrates.

(A) Cartoon representation of a chain of reversible reactions with metabolite cycling of two different metabolites, as shown in Equation 33. Each cycled metabolite has two forms, and there is no pathway independent cycling. (B) Heatmap of of the steady state concentration of M0 as a function of the metabolite pool sizes (Atot=Btot) and the inflow flux (kin). White area shows the region where there is no steady state. (C) Concentrations of M0-4, A0/A1 and B0/B1 ratios as a function of time. At t=1000 s, parameters are switched from the white dot in panel (B) (where a steady state exists) to the black dot (where we see build up of M0), and continued decline of A0. The B0/B1 ratio remains constant however as these are still cycled by reactions after the build up. In (B) and (C) the other parameters are kcat=100s-1, Etot=0.01mM, Km=50μM and kout=0.1s-1.

Appendix 6

Different stoichiometries for co-substrate cycling along a single pathway – mimicking the case seen in upper glycolysis

In this section, we consider a scenario of intra-pathway 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:

(35) 0kinM0M0+ATPM1+ADPM1+ATP2M2+ADPM2+ADPM3+ATPM3kout0.

We let A0=ATP and A1=ADP as usual.

With Michaelis-Menten kinetics, the ODE system is

dm0dt=kinE1L0m0a0F1K1m1a1K1m1a1+L0m0a0+K1L0dm1dt=E1L0m0a0F1K1m1a1K1m1a1+L0m0a0+K1L0E2L1m1a0F2K2m22a1K2m22a1+L1m1a0+K2L1dm2dt=2(E2L1m1a0F2K2m22a1)K2m22a1+L1m1a0+K2L1E3L2m2a1F3K3m3a0K3m3a0+L2m2a1+K3L2dm3dt=E3L2m2a1F3K3m3a0K3m3a0+L2m2a1+K3L2koutm3da0dt=E1L0m0a0F1K1m1a1K1m1a1+L0m0a0+K1L0E2L1m1a0F2K2m22a1K2m22a1+L1m1a0+K2L1+E3L2m2a1F3K3m3a0K3m3a0+L2m2a1+K3L2da1dt=E1L0m0a0F1K1m1a1K1m1a1+L0m0a0+K1L0+E2L1m1a0F2K2m22a1K2m22a1+L1m1a0+K2L1E3L2m2a1F3K3m3a0K3m3a0+L2m2a1+K3L2.

The system has two conservation laws:

a0+a1=Atot,m1+m2+a0=W.

By considering the equation 0=2dm0dt+2dm1dt+dm2dt+dm3dt=koutm3-2kin, we obtain

 m3=2kinkout.

Upon substitution of this value of m3 into dm3dt=0, dm2dt+dm3dt=0 and 2dm1dt+dm2dt+dm3dt=0, and solving recursively for M2,m1,m0, we obtain the following steady state relations:

m2=2kinK3((F3+2kin)a0+L2kout)kouta1L2(E32kin),m1=K2((F2+kin)a1m22+L1kin)L1a0(E2kin),m0=K1((F1+kin)a1m1+L0kin)a0L0(E1kin),

and recall a1=Atot-a0. We see that if 0<a0<Atot, these expressions are all positive if and only if

(36) kin<min(E1,E2,E3/2).

As usual, we consider the remaining conservation law, m1+m2+a0=W, together with these expressions, to find a polynomial p(a0) 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 a0=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

α=L2kout(L1L2(E32kin)2kout(AtotW(E2kin)+K2kin)2K3L1L2(E2kin)(E32kin)kinkout8K2K32kin2(F3+2kin)(F2+kin)),

which depends on Atot and W. To summarize for positive steady states to exist, we need

 kin<E1,kin<E2,kin<E32,α>0.

If kin is small enough, then the conditions hold as α>0 at kin=0.

However, these conditions are not sufficient. To see this, by inspecting the term α, one can see that if K2,K3 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,F1=1,F2=2,F3=3,L0=1L1=1,L2=2,K1=1kin=1,kout=1,E1=2,E2=2,E3=3,

and note that (35) holds. With K2=K3=0.1, the polynomial of interest becomes p(a0)=4a03-14.300a02+7.360a0-0.448, and so α>0. The polynomial has two positive roots under Atot=1, namely 0.07 and 0.537. For K2=0.2,K3=0.3, the polynomial is p(a0)=4a03-23.400a02+2.080a0-1.664, and although α>0, it has no root in the interval (0,1).

An analogous reasoning holds for the irreversible system.

Appendix 6—figure 1
Motif, time series and threshold for a model with differential co-substrate stoichiometry, as seen in glycolysis.

(A) Reaction system modelled in Appendix 6. Branching arrowhead from M1 to M2 indicates that two M2 molecules are produced/used in the forward/backward reactions. (B) Heatmap of the steady state M0 concentration. In the white area there is no steady state due to continual build up of M0. Dashed line shows the smallest limit imposed by Equation 35. (C) Time series of M04 and A0/A1. At t=1000 s, parameters are switched from the white dot in panel (B) (where a steady state exists) to the black dot (where we see build up of M0). Note that in the build up regime, M2 reduces as well as A0/A1, as M2 production is dependent on the presence of A0. In (B) and (C) the other parameters are kcat=100s-1, Etot=0.01mM, Km=50μM and kout=0.1s-1.

Appendix 7

Co-substrate cycling along with metabolite cycling – mimicking the case seen in nitrogen assimilation

In this section, we consider a scenario of intertwined, co-substrate cycling within a pathway. In particular, we consider a simplified version of the case seen in nitrogen assimilation where NAD(P)H cycling co-occurs together with a cycling from alpha-ketoglutarate to glutamate (mediated by glutamate dehydrogenase) and back from glutamate to alpha-ketoglutarate and glutamine (mediated by glutamate synthase). The simplified reaction scheme we consider comprises:

(37) 0kinNH3M0+NADH+NH3M1+NADM2kout0M1+NADM0+M2+NADH.

Denoting A0= NADH, A1= NAD, A2= NH3, the ODE system is:

dm0dt=E1L0a0a2m0F1K1a1m1L0a0a2m0+K1a1m1+K1L0+E2L1a1m1F2K2a0m0m2K2a0m0m2+L1a1m1+K2L1dm1dt=E1L0a0a2m0F1K1a1m1L0a0a2m0+K1a1m1+K1L0F2K2a0m0m2+E2L1a1m1K2a0m0m2+L1a1m1+K2L1dm2dt=F2K2a0m0m2+E2L1a1m1K2a0m0m2+L1a1m1+K2L1koutm2da0dt=E1L0a0a2m0F1K1a1m1L0a0a2m0+K1a1m1+K1L0+F2K2a0m0m2+E2L1a1m1K2a0m0m2+L1a1m1+K2L1da1dt=E1L0a0a2m0F1K1a1m1L0a0a2m0+K1a1m1+K1L0F2K2a0m0m2+E2L1a1m1K2a0m0m2+L1a1m1+K2L1da2dt=kinE1L0a0a2m0F1K1a1m1L0a0a2m0+K1a1m1+K1L0.

There are three conservation laws:

 a0+m1=M,a0+a1=Atot,m0+a1=W.

From da2dt+dm1dt+dm2dt=0 we get

 m2=kinkout

as expected. From dm2dt=0 and da2dt=0 we get

 m1=kinK2(F2a0m0+kina0m0+L1kout)kouta1L1(E2kin),a2=K1(F1a1m1+a1kinm1+L0kin)m0a0L0(E1kin).

From the second and third conservation laws, we get

 a1=Atota0,m0=Wa1=WAtot+a0.

This gives that for a positive steady state, we need kin<E1,kin<E2 and max(AtotW,0)<a0<Atot. Note that Atot-W is the constant value of a0-m0 along trajectories.

Plugging these expressions into the first conservation law, we have that steady states are in one-to-one correspondence with the solutions to

(38)  M=a0+kinK2(F2a0(a0Atot+W)+a0kin(a0Atot+W)+L1kout)kout(a0+Atot)L1(E2kin).

We first note that the derivative of the right hand side of (38) with respect to a0 is:

1+kinK2(F2(a0A tot+W)+F2a0+kin(a0Atot+W)+a0kin)kout(a0+Atot)L1(E2kin)+kinK2(F2a0(a0Atot+W)+a0kin(a0Atot+W)+L1kout)kout(a0+Atot)2L1(E2kin)

As Atota0>0 and a0AtotW>0, this function is positive in the interval of interest. Therefore, the right hand side of (38) is an increasing function of a0 when max(AtotW,0)<a0<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(AtotW,0)<a0<Atot of the following polynomial

p(a0)=(L1(E2kin)kout+K2kin(F2+kin))a02+(L1(E2kin)kout(M+T)+(WAtot)K2kin(F2+kin))a0L1kout(MAtot(E2kin)K2kin).

When a0=Atot, the polynomial is positive.

Case 1: Atot-W0. In this case we want 0<a0<Atot. If the independent term of p is negative, then p has exactly one solution in (0,Atot). So, if MT(E2kin)K2kin<0, we have one positive steady state. This is equivalent to

 kin<MAtotE2MAtot+K2.

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: AtotW>0. In this case, we want AtotW<a0<Atot. We find that p(Atot-W) is

 p(AtotW)=L1kout(W(MAtot+W)(E2kin)K2kin).

If this is negative, then there is one positive steady state, as p is positive at Atot. The condition is

 kin<WE2(MAtot+W)W(MAtot+W)+K2.

Note that M+W-Atot=a0+m1+m0+a1-a0-a1=m0+m1, 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 kin<E1 and

 kin<MAtotE2MAtot+K2andAtotW0,

or

 kin<WE2(MAtot+W)W(MAtot+W)+K2,AtotW>0,andM+WAtot>0.
Appendix 7—figure 1
Motif, time series and threshold for model with both co-substrate and metabolite cycling, mimicking nitrogen assimilation.

(A) Reaction system modelled in Appendix 7. (B) Heatmap of the steady state A2 concentration. In the white area there is no steady state due to continual build up of A2. Dashed line shows the limit kin<MTE2MT+K2, which is the smallest limit for these parameters. (C) Time series of M02, A2 and A0/A1. At t=1000 s, parameters are switched from the white dot in panel (B) (where a steady state exists) to the black dot (where there is continual build up of A2). In (B) and (C) the other parameters are kcat=100s-1, Etot=0.01mM, Km=50μM and kout=0.1s-1.

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 8—figure 1
Measured flux values (from Davidi et al., 2016; Gerosa et al., 2015) plotted against the calculated primary enzyme kinetic threshold (first part of Eq. (1) of the main text) adjusted by substrate affinity of the enzyme.

Note that the flux data shown here is a subset of the flux data presented in (Figure 3 of the main text), focusing only on those where the main substrate concentration was experimentally measured and the relevant KM is known. For both panels, the solid line indicates the equivalence of the two values and the dashed lines indicate 10-fold change interval on this, as a guide to the eye. Point colour indicates the nature of co-substrate involved and fill state indicates the data source (as shown on the inset).

Appendix 8—figure 2
Cumulative distribution of the ratio of observed flux (measured or FBA-predicted) to enzyme kinetic based limit (Vmax) (left panel) or to enzyme kinetic based limit accounting for substrate affinity (right panel).
Appendix 8—figure 3
Measured flux values (from Davidi et al., 2016; Gerosa et al., 2015) plotted against ΔGm which is the standard Gibbs energy of the reaction assuming all metabolites are at 1 mM concentrations (note this can be different than the usual definition of ΔG where concentrations are set to 1 M – a concentration that is not reflective of standard physiological conditions).

(A) shows a scatter plot of the ratio between the observed flux and the enzyme kinetic based limit (Vmax) against the ΔGm. The points highlighted in orange are for the malate dehydrogenase reaction. (B) A scatter plot of the enzyme kinetic based limit against the measured flux (same data as in Figure 3 of the main text), where the ΔGm is shown using a colourmap.

Appendix 8—figure 4
Measured flux values under different experimental conditions (from Davidi et al., 2016; Gerosa et al., 2015) for select reactions plotted against the corresponding co-substrate pool size.

See Supplementary file 1 for short notations for reactions. Each point on each panel is a separate flux measurement under a different environmental condition, where the co-substrate pool size is also measured. Point colours represent co-substrate type, with red for AMP +ADP + ATP, blue for NAD +NADH, and green for NADP +NADPH. Normalised RMSE of the best linear fit and the p-value are shown in the panel title.

Appendix 8—table 1
Correlation coefficients between observed flux (measured or FBA-predicted) and co-substrate pool size.

For reaction ID and descriptions, see Supplementary file 1.

ReactionFlux sourceNumber of conditionsPearson-rp-Value
25MDHGerosa70.8743100.010042
33PGKGerosa40.9735470.026453
03OAR140Davidi70.7943770.032850
7DHFSDavidi70.7877190.035438
46UAMASDavidi70.7837670.037025
45UAMAGSDavidi70.7837670.037025
44UAAGDSDavidi70.7837670.037025
48UGMDDSDavidi70.7837670.037025
2AIRC2Davidi70.7835230.037125
39PRASCSiDavidi70.7835230.037125
27NADKDavidi70.7835230.037125
19HSKDavidi70.7830420.037321
4CTPS2Davidi70.7823030.037624
41PTPATiDavidi70.7817950.037833
35PNTKDavidi70.7817950.037833
8DTMPKDavidi70.7801770.038501
43TMPKDavidi70.7774540.039643
34PMPKDavidi70.7774540.039643
1ADSKDavidi70.7774180.039658
15GLU5KDavidi70.7744610.040918
37PRAGSrDavidi70.7739000.041159
40PRFGSDavidi70.7739000.041159
38PRAISDavidi70.7739000.041159
14GK1Davidi70.7729180.041584
6DBTSDavidi70.7719000.042026
5CYTK1Davidi70.7621250.046410
20ICDHyrDavidi70.7430030.055681
47UAPGRDavidi70.7027990.078196
9G5SDDavidi70.6928300.084415
28P5CRDavidi70.6928300.084415
12GAPDDavidi7–0.6812030.091988
11G6PDH2rGerosa70.6299820.129423
10G6PDH2rDavidi70.6179090.139204
16GNDDavidi70.6179090.139204
22IMPDDavidi7–0.5399180.210941
23IPMDDavidi7–0.5359730.214953
36PPNDDavidi7–0.5205160.231016
26MTHFR2Davidi7–0.5190490.232569
32PGCDDavidi7–0.5174750.234241
3AKGDHGerosa70.4986660.254643
17GNDGerosa70.4986360.254676
42PYKGerosa4–0.7188210.281179
18HISTDDavidi7–0.4693230.288020
30PFKDavidi50.5372600.350444
24MDHDavidi70.3114710.496502
13GAPDGerosa4–0.4518170.548183
21ICDHyrGerosa7–0.1695100.716347
29PDHGerosa70.0980150.834403
31PFKGerosa2nannan

Appendix 9

Pathway branching into two pathways with independent co-substrates

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:

(39) 0kinM0M0M1,2M0M1,1M1,2+A0M2,2+A1M1,1+B0M2,1+B1A1A0B1B0M2,2kout,20M2,1kout,10.

The ODE system is:

(40) dm0dt=kinE1,1L0m0F1,1K1,1m1,1K1,1L0+K1,1m1,1+L0m0E1,2L0m0F1,2K1,2m1,2K1,2L0+K1,2m1,2+L0m0dm1,2dt=E1,2L0m0F1,2K1,2m1,2K1,2L0+K1,2m1,2+L0m0E2,2L1,2m1,2a0F2,2K2,2m2,2a1K2,2L1,2+K2,2m2,2a1+L1,2m1,2a0dm1,1dt=E1,1L0m0F1,1K1,1m1,1K1,1L0+K1,1m1,1+L0m0E2,1L1,1m1,1b0F2,1K2,1m2,1b1K2,1L1,1+K2,1m2,1b1+L1,1m1,1b0dm2,2dt=E2,2L1,2m1,2a0F2,2K2,2m2,2a1K2,2L1,2+K2,2m2,2a1+L1,2m1,2a0m2,2kout,2dm2,1dt=E2,1L1,1m1,1b0F2,1K2,1m2,1b1K2,1L1,1+K2,1m2,1b1+L1,1m1,1b0m2,1kout,1da0dt=EaLaa0FaKaa1KaLa+Kaa1+Laa0E2,2L1,2m1,2a0F2,2K2,2m2,2a1K2,2L1,2+K2,2m2,2a1+L1,2m1,2a0da1dt=EaLaa0FaKaa1KaLa+Kaa1+Laa0+E2,2L1,2m1,2a0F2,2K2,2m2,2a1K2,2L1,2+K2,2m2,2a1+L1,2m1,2a0db0dt=EbLbb0FbKbb1KbLb+Kbb1+Lbb0E2,1L1,1m1,1b0F2,1K2,1m2,1b1K2,1L1,1+K2,1m2,1b1+L1,1m1,1b0db1dt=EbLbb0FbKbb1KbLb+Kbb1+Lbb0+E2,1L1,1m1,1b0F2,1K2,1m2,1b1K2,1L1,1+K2,1m2,1b1+L1,1m1,1b0.

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, kin along with the ratio of the pool sizes, Atot/Btot in the main text (Figure 4 and Appendix 9—figure 1), while varying kin with the ratio of moiety back-cycling rates ka/kb is presented in Appendix 9—figure 2.

Appendix 9—figure 1
Effect of varying Atot/Btot and kin in the branching system shown in Figure 4A, on the ratio of flux into M2,2/M2,1 for reversible and irreversible dynamics and different Vmax,A and Vmax,B values.

Downstream metabolite has higher flux when it’s respective co-substrate has a higher concentration. (A) Reversible reactions with Vmax,A=Vmax,B=0.1 (i.e. one order of magnitude smaller than the other reaction rates). (B) Irreversible reactions with Vmax,A=Vmax,B=0.1 (i.e. one order of magnitude smaller than the other reaction rates). (C) Reversible reactions with Vmax,A=Vmax,B=10 (i.e. one order of magnitude larger than the other reaction rates).

Appendix 9—figure 2
Effect of varying Vmax,B/Vmax,A (i.e maximum back-cycling) ratio and kin in the branching system shown in Figure 4A.

(A) Motif showing the dynamics. (B) Heatmap showing the ratio of flux into M2,2/M2,1. Downstream metabolite has higher flux when the back-cycling rate of its co-substrate is higher, with a greater effect for higher kin values. (C) Time series of downstream metabolites and shared precursor. Solid/dashed lines show parameters for black/white dot in (B).

Appendix 10

Independent pathways coupled by co-substrate 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 A0A1 in the ‘forward’ direction of the first pathway, while A1A0 in the ‘forward’ direction of the other. The metabolite reactions are reversible, and there is pathway-independent cycling of the shared conserved moiety. The motif is shown in Figure 5A in the main text, and the reactions are as follows:

(41) 0kin,1M0,1M0,1+A0M1,1+A1M1,1kout,100kin,2M0,2M0,2+A1M1,2+A0M1,2kout,20A0A1.

The resulting system of ODEs is:

(42) dm0,1dt=kin,1E1,1L0,1m0,1a0F1,1K1,1m1,1a1K1,1L0,1+K1,1m1,1a1+L0,1m0,1a0dm1,1dt=E1,1L0,1m0,1a0F1,1K1,1m1,1a1K1,1L0,1+K1,1m1,1a1+L0,1m0,1a0kout,1m1,1dm0,2dt=kin,2E1,2L0,2m0,2a1F1,2K1,2m1,2a0K1,2L0,2+K1,2m1,2a0+L0,2m0,2a1dm1,2dt=E1,2L0,2m0,2a1F1,2K1,2m1,2a0K1,2L0,2+K1,2m1,2a0+L0,2m0,2a1kout,2m1,2da0dt=E1,2L0,2m0,2a1F1,2K1,2m1,2a0K1,2L0,2+K1,2m1,2a0+L0,2m0,2a1E1,1L0,1m0,1a0F1,1K1,1m1,1a1K1,1L0,1+K1,1m1,1a1+L0,1m0,1a0EaLaa0FaKaa1KaLa+Kaa1+Laa0da1dt=E1,2L0,2m0,2a1F1,2K1,2m1,2a0K1,2L0,2+K1,2m1,2a0+L0,2m0,2a1+E1,1L0,1m0,1a0F1,1K1,1m1,1a1K1,1L0,1+K1,1m1,1a1+L0,1m0,1a0+EaLaa0FaKaa1KaLa+Kaa1+Laa0.

This system has one conservation law

 a0+a1=Atot,

giving that at steady state, a1=Atot-a0. By solving the steady state equations dm0,1dt=dm1,1dt=dm0,2dt=dm1,2dt=0 for m0,1,m0,2,m1,1,m1,2, we obtain

m0,1=K1,1kin,1((F1,1+kin,1)a1+L0,1kout,1)kout,1L0,1a0(E1,1kin,1),m1,1=kin,1kout,1,m0,2=K1,2kin,2((F1,2+kin,2)a0+L0,2kout,2)kout,2L0,2a1(E1,2kin,2),m1,2=kin,2kout,2.

The expressions are positive provided E1,1>kin,1, E1,2>kin,2, and a0,a1>0.

Finally, we use da0dt=0 to solve for a0 and obtain

(43)  a0=Ka((kin,2kin,1)(La+Atot)+AtotFa)(kin,2kin,1)(KaLa)+EaLa+FaKa.

For all quantities to be positive at steady state, we require 0<a0<Atot. By subtracting from Atot the value of a0 at steady state, we obtain

(44)  Atota0=La((Atot+Ka)(kin,1kin,2)+AtotEa).

To summarize, there is a positive steady state if and only if E1,1>kin,1, E1,2>kin,2, and (44) and (43) are positive. For(44) to be positive,

(45) AtotEaAtot+Ka>kin,2kin,1,that isAtotEaAtot+Ka<kin,1kin,2,

and either

(46) AtotFaAtot+La>kin,1kin,2,EaLa+FaKa>(kin,1kin,2)(KaLa)

or

(47) AtotFaAtot+La<kin,1kin,2,EaLa+FaKa<(kin,1kin,2)(KaLa)

By analysing these cases, we obtain all possible scenarios for a positive steady state to exist, and these are dictated by the difference kin,1-kin,2.

For example, if kin,1=kin,2, then (45) and (46) hold directly. If Ka=La, then we require AtotFaAtot+La>kin,1kin,2>AtotEaAtot+Ka to hold.

If Ka>La, the conditions lead to either

min(AtotFaAtot+La,EaLa+FaKa(KaLa))>kin,1kin,2>AtotEaAtot+Ka

or

max(AtotFaAtot+La,EaLa+FaKa(KaLa),AtotEaAtot+Ka)<kin,1kin,2.

If Ka<La, the conditions lead to either

AtotFaAtot+La>kin,1kin,2>max(AtotEaAtot+Ka,EaLa+FaKa(KaLa))

or

max(AtotFaAtot+La,AtotEaAtot+Ka)<kin,1kin,2<EaLa+FaKa(KaLa).

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 kin 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.

Appendix 10—figure 1
Comparison of the onset of instability in the coupled pathway case vs the single pathway case.

(A) Pathway enzyme and cycling enzyme limits in the single and coupled pathway cases (see legend, pathway enzyme limit is the same for each case). For the coupled pathway, kin,2 is set to.0.9kin,1 Since the limiting factor depends on the difference between kin,1 and kin,2, higher kin values are possible without build up of upstream metabolites if all other parameters are kept the same. (B) Time series of M0, M1 and A0/A1 in the single pathway case where the parameters for Atot and kin are set by the black dot in (A). Here we see build up M0 because kin>AtotFa/(Atot+La). (C) Time series of M0,1, M0,2, M1,1, M1,2 and A0/A1 in the coupled pathway case where the parameters for Atot, kin,1 and kin,2 are set by the black dot in (A). Here we see the system admits a steady state because.kin,1kin,2<AtotFa/(Atot+La) The parameters in all panels are: AtotFaAtotLa>kin,1kin,2>max(AtotEaAtotKa,EaLa+FaKa(KaLa)),K1=L0=K1,1=K1,2=L0,1=L0,2=1M,kout=kout,1=kout,2=0.1s-1,Ea=Fa=15Ms-1.Ka=La=1 M.

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 log-ratio (i.e.log10kin,1/kin,2) from a standard normal distribution with mean μ=0 and variance σ2=1. A new log-ratio is drawn after waiting a time that is drawn from an exponential distribution with mean τ. Example time-series are shown in Figure 5(C and D) of the main text. The log-ratio 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 discrete-time analogue of the Ornstein-Uhlenbeck process: it has the same steady state distribution, mean, variance and correlation function as its continuous-time 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 τ and total pool size Atot 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 anti-correlated. Furthermore, this behaviour is observed whether the metabolite reactions are irreversible or reversible. Since the transition from correlated to anti-correlated is so sharp, relatively small changes in the total amount of cycled co-substrate, combined with noisy influx rates, could lead to the downstream metabolites changing from being correlated to anti-correlated, or vice-versa.

Appendix 10—figure 2 also further demonstrates the effect of changing Ea, τ and Atot on the time series of the products M1,1 and M1,2 and the co-substrate ratio A0/A1. For parameters represented by the black dot in the first row (panels (E) - (H)), the products are correlated when Ea1, and become uncorrelated when Ea=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 Atot. For parameters represented by the grey dot in the first row (panels (I) - (L)), the products switch from correlated to anti-correlated as Ea increases, and remain responsive to the noise. The variation in A0/A1 reduces from left to right because Ea increases. For parameters represented by the white dot in the first row (panels (M) - (P)), the products remain anti-correlated and responsive to the noise. The variation in A0/A1 reduces from left to right but is not as pronounced for small Ea because Atot is very large.

Appendix 10—figure 2
Correlation between products of two coupled pathways, using a model with irreversible reaction rates.

(A) - (D) show results for Ea=Fa=0,0.01,1 and 10, respectively. As in the Figure 5 of the main text, all other parameters are set to 1, apart from the kout values that are each 0.5. As Ea increases, the products become anti-correlated at smaller Atot values. This is also true for the model using reversible reaction rates (see Figure 5 of the main text). (E) - (P) show time series for the products M1,1 and M1,2 (purple and green respectively) on the left axis and the ratios A0/A1 and kin,2/kin,1 (blue and pink respectively) on the right axis. Panels in the same column share the same Ea value as the heatmap in the top (e.g. Ea=0 in (E), (I) and (M)). Panels (E) - (H), (I) - (L) and (M) - (P) have Atot and τ given by the black, grey and white dots respectively. See text.

Appendix 11

Table of parameters

Appendix 11—table 1
List of parameters together with the biological meanings and units.

M and t stand for Molar and time, respectively. Where different notation is used between this and the main text, the equivalent parameters have been made clear.

ParameterBiological meaning and notesUnits
Ei, FiCatalytic rate coefficients for forward and backward reactions respectively. In the main text we only use the forward rate and denote it Vmax,Ei.Mt-1
Ki, Li-1Michaelis-Menten coefficients for forward and backward reactions, respectively. In the main text, we only use the forward rate and denote it KM,Ei.Mn, where n is the number of different substrates involved in the reaction.
kinInflux rate.Mt-1
koutOutflow rate.t-1

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/OSS-Lab/CoSubstrateDynamics/tree/v1.0.0.

The following data sets were generated
    1. West R
    2. Delattre H
    3. Noor E
    4. Feliu E
    5. Soyer O
    (2023) Zenodo
    Co-substrate dynamics analysis software and data.
    https://doi.org/10.5281/zenodo.7565439

References

  1. Book
    1. Reich JG
    2. Sel‘kov EE
    (1981)
    Energy Metabolism of the Cell: A Theoretical Treatise
    Academic press.

Article and author information

Author details

  1. Robert West

    School of Life Sciences, University of Warwick, Warwick, United Kingdom
    Contribution
    Conceptualization, Resources, Software, Formal analysis, Writing – review and editing
    Contributed equally with
    Hadrien Delattre
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9348-0258
  2. Hadrien Delattre

    School of Life Sciences, University of Warwick, Warwick, United Kingdom
    Contribution
    Conceptualization, Investigation, Writing – review and editing
    Contributed equally with
    Robert West
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5488-8370
  3. Elad Noor

    Department of Plant and Environmental Sciences, Weizmann Institute of Science, Rehovot, Israel
    Contribution
    Resources, Formal analysis, Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8776-4799
  4. Elisenda Feliu

    Department of Mathematics, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Formal analysis, Methodology, Writing – review and editing
    For correspondence
    efeliu@math.ku.dk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7205-6511
  5. Orkun S Soyer

    School of Life Sciences, University of Warwick, Warwick, United Kingdom
    Contribution
    Conceptualization, Supervision, Investigation, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    O.Soyer@warwick.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9504-3796

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

  1. Preprint posted: September 6, 2022 (view preprint)
  2. Received: October 22, 2022
  3. Accepted: February 16, 2023
  4. Accepted Manuscript published: February 17, 2023 (version 1)
  5. 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

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. Robert West
  2. Hadrien Delattre
  3. Elad Noor
  4. Elisenda Feliu
  5. Orkun S Soyer
(2023)
Dynamics of co-substrate pools can constrain and regulate metabolic fluxes
eLife 12:e84379.
https://doi.org/10.7554/eLife.84379

Share this article

https://doi.org/10.7554/eLife.84379

Further reading

    1. Computational and Systems Biology
    Qianmu Yuan, Chong Tian, Yuedong Yang
    Tools and Resources

    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 genome-scale applications. Besides, these methods haven’t fully explored the geometry of the protein structures. Here, we propose GPSite, a multi-task 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 end-to-end manner. Experiments demonstrate that GPSite substantially surpasses state-of-the-art sequence-based and structure-based approaches on various benchmark datasets, even when the structures are not well-predicted. The low computational cost of GPSite enables rapid genome-scale 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://bio-web1.nscc-gz.cn/app/GPSite.

    1. Cell Biology
    2. Computational and Systems Biology
    Thomas Grandits, Christoph M Augustin ... Alexander Jung
    Research Article

    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 speed-ups 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 root-mean-square 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.