<figure><img src="image_1.png" /></figure>

HAL
open science

Cell metabolism

Hadrien Delattre, Wolfram Liebermeister, Elad Noor, Herbert Sauro, Orkun

Soyer, Robert West

► To cite this version:

Hadrien Delattre, Wolfram Liebermeister, Elad Noor, Herbert Sauro, Orkun Soyer, et al. Cell metabolism. Economic Principles in Cell Biology, The Economic Cell Collective, 2025, 10.5281/zen-odo.8156823. hal-05195020

HAL Id: hal-05195020

https://hal.inrae.fr/hal-05195020v1

Submitted on 31 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Cell metabolism

Hadrien Delattre<fnref n="1" /> ID, Wolfram Liebermeister² ID, Elad Noor³ ID, Herbert M. Sauro⁴ ID, Orkun S. Soyer<fnref n="1" /> ID, and Robert West<fnref n="1" /> ID

## Abstract

Metabolism is a dynamical process, in which a myriad of chemicals are converted via biochemical reactions inside the cell over time. This chapter motivates the dynamical nature of metabolism and introduces mathematical approaches for modeling it. Basic models of biochemical reaction rates are introduced and their assumptions and relations to each other are explained. Using these rate models, simple metabolic systems are modeled and their dynamical behavior analyzed. Experimental findings on metabolic dynamics are summarized and possible explanatory models of such dynamics are discussed.

**Keywords:** metabolism – metabolic pathway – metabolite – enzyme kinetics – metabolic network – biochemical reactions – thermodynamics – steady state

**Contributions:** The chapter was originally drafted and written by Orkun Soyer, Hadrien Delattre, and Robert West. Later, sections were added by Wolfram Liebermeister, Elad Noor, and Herbert Sauro. We thank Francesco Moro, Markus Arthur Köbis, and Maarten Droste for comments on an earlier version of this chapter. We acknowledge funding from the UK's Biotechnology and Biological Sciences Research Council (BBSRC) grant ID BB/T010150/1 and from the Gordon and Betty Moore Foundation grant ID GBMF9200.

**To cite this chapter:** Hadrien Delattre, Elad Noor, Herbert M. Sauro, Orkun S. Soyer, Robert West. Cell metabolism (Version July 2025). doi: 10.5281/zenodo.8156823. Chapter from: The Economic Cell Collective (2025). Economic Principles in Cell Biology. No commercial publisher | Online open access book | doi: 10.5281/zenodo.8156386

The authors are listed in alphabetical order.

<figure><img src="image_3.png" /></figure>

This is a chapter from the open textbook "Economic Principles in Cell Biology".
Free download from [principlescellphysiology.org/book-economic-principles/](http://principlescellphysiology.org/book-economic-principles/).
Lecture slides for this chapter are available on the website.

© 2025 The Economic Cell Collective.

Licensed under Creative Commons License CC-BY-SA 4.0.
An online open access book. No publisher has been paid.
doi: 10.5281/zenodo.8156386

<figure><img src="image_4.png" /></figure>

<footnote n="1">¹ School of Life Sciences, University of Warwick, Gibbet Hill Road, Coventry, UK; ² Université Paris-Saclay, INRAE, MalAGE, 78350 Jouy-en-Josas, France; ³ Department of Plant and Environmental Sciences, Weizmann Institute of Science, 76100 Rehovot, Israel; ⁴ University of Washington, Seattle, WA, US; ⁵ Université Paris-Saclay, INRAE, MalAGE, 78350 Jouy-en-Josas, France.</footnote>
## Chapter overview

This chapter introduces cell metabolism as a dynamical system. While the previous chapter gave an overview of the constituents of this system, i.e. enzymes, metabolites, etc., this chapter focuses on conceptual abstraction of the metabolic system as a whole and how to model its dynamics over time. The key areas introduced are:

* Conceptualizing cell metabolism as a dynamical system

* Dynamics and regulation of metabolism

* Toolbox for modeling dynamics of metabolism - biochemical reaction rate laws and their derivations

* Dynamics of metabolism: Examples of experimental evidence and model-based explanations

* Mathematical derivations and example models

### 3.1. Conceptualizing cell metabolism as a dynamical system

Cell metabolism is a dynamical process that converts available metabolites from the environment into biomass and other products. The metabolism of a typical cell involves thousands of biochemical reactions and metabolites. What would be a useful way to think about such a complex, dynamical system? We need a conceptual picture of metabolism to help us formulate more specific ideas about how it functions, how it can be manipulated, or even how it has evolved. Here, we first highlight a few such 'pictures', or ways of thinking about metabolism.

Below we will switch back-and-forth between a high-level view on metabolism, considering all of it, and a more focused, low-level view focusing on modeling individual reactions or small sets of reaction systems (e.g. pathways or motifs). These two viewpoints constitute two ends of a wide spectrum, and our aim in jumping back-and-forth between them is to allow the reader to obtain the skills to model dynamics of reaction systems that make up metabolism, while at the same time to invite them to think about the overall function of the metabolic system.

#### 3.1.1. Metabolism as a collection of pathways

The common and historical view of a metabolic system stems from pioneering biochemical studies from the 1930s onwards, which identified collections of reactions as so-called pathways [1]. Known mostly through the names of their discoverers, these include the EntnerDoudoroff (ED), EmbdenMeyerhofParnas (EMP) and pentose-phosphate (PP) pathways involved in glucose uptake and conversion into pyruvate, and the Krebs pathway (a.k.a. tricarboxylic acid cycle, TCA) involved in the conversion of pyruvate into amino acid, nucleotides, and biomass precursors [2]. This pathway-centric view of cell metabolism lends itself readily to an assembly line analogy and the notion of (linearly) connected pathways (see Economic analogy 3.A).

**Pathways, yes, but not so linear!** The identification of well-established pathways and the subsequent focus upon them gives the false impression that cell metabolism consists of a series of neatly organized and serially connected pathways. This impression is facilitated by pictures of isolated linear pathways, common in textbooks and even research papers. In reality, these pathways are highly interconnected with other pathways.

Part of these interconnections within metabolism arise from co-substrates and specific metabolite pairs that participate in many reactions. For example, co-substrates such as ATP and NADH link many parts of metabolism through reactions in which they are generated or consumed (Fig. 3.2), while the glutamate - α-ketoglutarate pair is involved in the TCA cycle as well as acting as a group donor in all amino acid biosynthesis pathways.

The pathway view provides a useful starting point to think about metabolism, but a complete understanding of metabolism dynamics and metabolic phenotypes requires us to come to terms with the highly connected nature of these pathways (see below, Box 3.C).
<figure><img src="image_5.png" /><figcaption>Figure 3.1: A map of central metabolism in *Escherichia coli* bacteria – The diagram shows the reactions, metabolites, co-factors, and enzymes, as well as a few selected carbon sources and their catabolic pathways.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

Economics analogy 3.A Metabolism as an assembly line

We can make an analogy that presents metabolism as an assembly line in a factory. Metabolites enter the line from outside the cell and are processed i.e. acted upon by enzymes to create new metabolites that are ultimately incorporated into cellular biomass. This picture is reinforced by the common textbook illustration of metabolism as a set of isolated pathways that are placed upstream or downstream of each other, and that produce or consume outputs for each other. A key shortcoming of this analogy is that it conveys a picture in which events are strictly linear and progressive in their nature, ignoring the cyclic and inter-connected nature of metabolism. Despite this shortcoming, this analogy captures the point that the flux of materials through the system can attain a steady-state of equal in- and out-flux across individual reactions (see further discussion of the steady-state concept in the main text). One important difference however between an assembly line and metabolism is that the rate at a given assembly stage in a factory is not a function of how many units are waiting to be processed because factory machines tend run at fixed rates. In metabolism, the rate of a reaction is a function of the substrate concentration until saturated. This leads to distinctive behavior not found in factory assembly lines. Another important difference with a factory assembly line is that unlike an assembly line, metabolism in some cases is able to in both directions along the line. The most well known of these is the bidirectionality of the glycolytic and gluconeogenic pathways.

3.1.2. Coarse grained views of metabolism

The highly connected nature of metabolism makes it difficult to understand its overall dynamics just from individual pathways. It also makes it hard to conceptualize metabolism as a single, linear process, or as serially connected pathways. Here, a coarse-grained viewpoint, focusing on the overall function of cell metabolism, might prove helpful. There have been several such views developed, with two highlighted here.

**Metabolism as biomass generator.** A widely applied coarse-grained view of metabolism considers it as a vehicle to biomass production. In this view, metabolism is considered as two coupled processes, one producing energy and compounds that can act as building blocks (e.g. amino acids), and one that uses these to create larger macro molecules (e.g. proteins and lipids) needed to make a new cell. These two processes are called catabolic and anabolic

<figure><img src="image_7.png" /><figcaption>Figure 3.2: A simplified map of central metabolism, particularly highlighting interconnections among different pro- cesses (i.e. pathways) through the NAD(P)⁺ / NAD(P)H co-substrate pair.</figcaption></figure>
metabolism respectively, and their coupling presents the whole cell metabolism (Fig. 3.3 A). This coarse-grained model is widely used (e.g. [2, 3]). However, it is not always clear how to partition various pathways and reactions as anabolic and catabolic, and the notion of metabolism organized solely to satisfy for biomass production does not capture certain metabolic phenotypes, such as no-growth states or excretion of high-energy metabolites (i.e. metabolic overflow).

**Metabolism as electron flow.** An alternative coarse-grained view of metabolism is obtained from a more chemical standpoint. When one writes down an overall reaction for cellular metabolism, considering compounds taken up from the environment and created at the end of various metabolic processes, one realizes that this is a redox reaction, a type of reaction where electrons are exchanged between participating reactants (see Fig. 3.3 B and Box 3.B). This means that the actual reactions within metabolism that enable this overall reaction must compose also of some redox reactions. In other words, we can argue that metabolism consists of (besides other reactions) a series of redox reactions that enable flow of electrons. Metabolism is thus an inter-connected system of reactions that allows flow of electrons from readily oxidized compounds (electron rich compounds with low or negative reduction potentials) towards readily reduced compounds (electron poor compounds with positive reduction potentials) [4, 5]. (Fig. 3.3 B). As the Nobel laureate Albert Szent-Györgyi (1893–1986), who studied the TCA cycle and discovered vitamin C biosynthesis pathways, once said, "Life is an electron looking for a place to rest."

Emphasizing its redox reactions, the metabolic system can be visualized on a reduction potential chart, which is sometimes called a 'redox ladder' (Fig. 3.4 and box 3.B). This potential chart shows reduction potential of redox half reactions (usually in reduction direction) and allows us to readily visualize the thermodynamic feasibility of redox reaction pairs. The chart is ordered in such a way that any reduction half reaction can be paired with any other placed below it, resulting in a thermodynamically feasible redox reaction, but not with those above it. We notice that cell metabolism, in order to maintain electron flows, needs to maintain thermodynamic feasibility of the overall and all intermediate reactions. The key requirement for this is to have access to electron donors (e.g. carbohydrates) and terminal electron acceptors (e.g. oxygen). One must also note that the redox ladder depicted in Fig. 3.4 is derived for standard concentrations of metabolites, whereas the reduction potentials would depend on actual concentrations in the cell.

### 3.1.3. Keeping flows in a system of interconnected fluxes

It is noticeable that both coarse-grained views presented above involve interconnected fluxes that ultimately enable an overall flux. In the biomass-based view, the flux between catabolism and anabolism is connected to enable flux into biomass. In the electron-flow based view, there is again a set of interconnected flows to enable the overall electron flow from initial donors (e.g. glucose) to final acceptors (e.g. oxygen).

The interconnection of fluxes in metabolism is most clearly visible in reactions involving co-substrates, such as NAD(P)⁺ / NAD(P)H and ADP/ATP pairs (see below, Box 3.C). The NAD(P)⁺ / NAD(P)H pairs form either the oxidation or reduction half-reaction in various redox reactions thereby enabling the aforementioned electron flows within the metabolic system. The ATP⁺/ADP pair forms an energy carrier, providing driving energy to reactions that would be thermodynamically infeasible (see section 3.2.1 below on what we mean by this). This pair is seen as forming the flux connection between catabolism and anabolism, where the former is considered to result in ATP production, and the latter is considered to consume this.

Co-substrates are thus essential in connecting different fluxes, and therefore processes, within metabolism and their dynamics must be important to keep overall metabolic flow. It is tempting to speculate that key co-substrates might be an evolutionary outcome that ensures stable electron flows in the face of changing conditions. While this possibility is difficult to prove or disprove, it is interesting to note that the NAD(P)H/NAD(P)⁺ pairs can attain a broad range of reduction potentials that could enable their redox partnering with many of the different reaction types found in
<figure><img src="image_8.png" /><figcaption>Figure 3.3: Coarse-grained models of cell metabolism – (A) A conceptual drawing of cell metabolism as provider of precursors (catabolism) and generator of biomass from those (anabolism). (B) A conceptual drawing of cell metabolism as enabling an abstract redox reaction between a pair of electron donors and acceptors. The electron donor can at the same time be the carbon source for biomass generation, or there can be a separate ‘carbon-donor’. This overall redox reaction is an abstraction, in the sense that in real metabolism electrons are not directly transferred from the original donor to biomass precursors but rather there are many intermediary redox reactions such as those involving key carrier co-substrate metabolite pairs NAD(P)⁺/NAD(P)H.</figcaption></figure>

cell metabolism [7] - in other words, these two redox pairs seem to be a versatile tool to connect a wide range of redox reactions to each other and ensure electron flows.

### 3.1.4. Metabolic system and recurring motifs

Within the highly inter-connected system that is metabolism, specific reaction arrangements seem to recur frequently, so-called “reaction motifs”. We have already mentioned the cyclic reaction systems, involving co-substrates as one such motif. Other reaction motifs that have been highlighted include autocatalytic cycles [8] and branch points [9]. As we will discuss below, these reaction motifs can give rise to specific nonlinear dynamics and act in auto-regulatory capacity or create constraints on the metabolic system. In general, however, it is difficult to ascertain the evolutionary significance of reaction motifs. While automated approaches, involving graph theoretical analysis of metabolic systems represented as networks, highlighted certain metabolic motifs as significant compared to random networks, it was subsequently shown that this result is dependent both on the original network representation used and the randomized networks used for comparison [10].

## 3.2. Reaction thermodynamics and enzyme kinetics

Independent of our conceptual views on metabolism, the fact remains that the metabolic system involves flux of matter. A myriad of metabolites are combined, converted, broken apart, and re-assembled. These biochemical reactions are catalyzed by enzymes so to improve kinetic rates, and the entire system must obey the laws of thermodynamics (more on these later in section 3.2.1). In summary, metabolism constitutes a ‘system’ of metabolites and their reactions, together with enzymes. Its dynamics over time ensures fluxes of matter.
<figure><img src="image_9.png" /></figure>

Physics box 3.B The redox ladder in metabolism

We can highlight the overall redox reaction implemented by the cellular metabolism further, by writing it as two separate reactions consisting of an oxidation reaction (involving a molecule releasing electrons) and a reduction reaction (involving a molecule accepting electrons) (see Fig. 3.3). The feasibility of the paired, overall redox reaction can be measured by the Gibbs free energy, or the closely related reduction potential, where a positive reduction potential (or a negative Gibbs free energy) indicates a thermodynamically feasible reaction. Thus, a redox reaction with a positive reduction potential implies electrons ‘flowing’ from a molecule with high reduction potential towards that with a low reduction potential a point that can be visualized using a “reduction ladder”, a chart of reduction potentials (Fig. 3.4). Notice that considering redox reactions as composed of individual reduction and oxidation reactions is merely a conceptualization, however, this provides a useful analogy in which we can view a metabolic system as enabling the flux of electrons across many reactions, and between an initial electron donor and a final electron acceptor [1]. While glucose and oxygen are possibly the most well-known electron donor and acceptor pairs, cells, especially microbial cells, can use a wide-range of donors and acceptors, including nitrogen and sulfur containing compounds, thereby contributing significantly to biogeochemical cycles of these compounds [6].

<figure><img src="image_10.png" /></figure>

3.2.1. Biochemical reactions and thermodynamics

Metabolism consists of individual biochemical reactions of the form:

$$n_a A + n_b B \rightleftharpoons n_c C + n_d D \quad (3.1)$$

where $n_i$ are the so-called stoichiometric coefficients, determining the number of molecules of the i'th chemical species taking part in the reaction (Box 3.D). While these reactions are catalyzed by enzymes, they still need to obey thermodynamic laws. We will not provide a full treatise of the thermodynamics of chemical reactions here - we refer the reader to excellent books on physical chemistry for this (e.g. [16]) and also to books for a conceptual introduction to thermodynamics (e.g. [17]). Here, it suffices for us to define the key thermodynamic equation, the Gibbs free energy of reaction, involving the chemical potential of substrates and products. Chemical potentials are related to concentrations, where the relation depends on the ionic strength of the solution. Assuming an ideal solution, we will
<figure><img src="image_11.png" /><figcaption>Figure 3.4: Metabolism on a redox ladder – Cartoon representation highlighting the role of electron flows through redox reactions for a functioning metabolism, and a reduction potential chart listing key redox reactions found in cellular metabolism. Notice that the reduction potential chart shows reduction potentials of half-reactions in the reduction direction and using metabolite concentrations under standard conditions, hence the actual potentials would be different and dynamically changing within the cells. A thermodynamically feasible reaction would need to combine one half reaction (run in reverse, oxidation direction) with another one lying below it (i.e. at a higher reduction potential). Two example feasible redox pairs are shown with the blue and red data points.</figcaption></figure>

<figure><img src="image_12.png" /></figure>

Physics box 3.C Cyclic reaction motifs

The involvement of co-substrate and key metabolites results in the coupling of many different parts of the metabolism and in the emergence of cyclic reaction systems - for example, by connecting different parts of the metabolism, the NAD(P)H/NAD(P)⁺ pairs result in cycling between their different forms. This means that in order to capture the concentration of all the other molecules involved in these reactions, we need to consider dynamics of a series of intertwined cyclic reaction systems, rather than linear pathways akin to an assembly line. Indeed, it has been argued that cyclic reaction motifs should form the basis of developing a dynamic understanding of cell metabolism [11]. It must also be noted that co-substrates, and possibly other key metabolites, can have conserved concentrations in the time scales of metabolic flux dynamics. In other words, these metabolites form conserved moieties within the system, similar to enzymes, such that altering of the total pool size of these co-substrates or the ratio of their different forms (e.g. the NAD⁺/NADH ratio) can possibly affect the flux distribution across different pathways that they are connected to [12, 13, 8, 14, 11, 15].

write here the Gibbs free energy of reaction directly in terms of concentrations:

$$ \Delta_{\mathrm{r}}G' = \Delta_{\mathrm{r}}G'^{\circ} + R \cdot T \cdot \ln \underbrace{\frac{c^{n_{\mathrm{c}}} \cdot d^{n_{\mathrm{d}}}}{a^{n_{\mathrm{a}}} \cdot b^{n_{\mathrm{b}}}}}_{\Gamma}, \qquad (3.2) $$

where the small letters indicate the concentrations of the substrates and products as given in the above reaction. Notice that specifying ‘products’ and ‘substrates’ automatically specifies a ‘forward’ direction to the reaction (Box 3.D). In the above expression, the term in the natural logarithm is the ratio of the concentration of the products to the concentration of the substrates (considering the forward direction of the reaction) and is commonly denoted as the mass action ratio, Γ. The term Δₜₐₗ is the difference between the standard Gibbs free energy of formation of products and substrates.
The Gibbs free energy of a reaction is the key thermodynamic equation we introduce here, as it is this equation that determines whether a reaction would run in the forward direction or not. If the Gibbs free energy of reaction, for a given set of substrates and products concentration, is negative ($\Delta_r G' < 0$), the reaction will be spontaneous in the forward direction as it is written (i.e. in the way the 'substrates' and 'products' are defined). In other words, chemical reactions proceed in the direction of lower energy - they minimize the internal energy of the system. We will see later (in section 3.2.3) that Gibbs free energy will also feature in rate laws for biochemical reactions.

It is important to introduce here the concept of thermodynamic equilibrium, which is attained when $\Delta_r G' = 0$. Re-arranging equation 3.2 under this condition, we can obtain:

$$ \Delta_r G'^{\circ} = -R \cdot T \cdot \ln \frac{c_{eq}^{n_c} \cdot d_{eq}^{n_d}}{a_{eq}^{n_a} \cdot b_{eq}^{n_b}}, \quad (3.3) $$

where the subscript "eq" denotes the concentrations of each species at the thermodynamic equilibrium. The ensuing ratio is known as the equilibrium constant, $K_{eq} = \frac{c_{eq}^{n_c} \cdot d_{eq}^{n_d}}{a_{eq}^{n_a} \cdot b_{eq}^{n_b}}$. Re-arranging equation 3.3, we can derive an expression for $K_{eq}$ as follows:

$$ K_{eq} = e^{\frac{-\Delta_r G'^{\circ}}{R \cdot T}} \quad (3.4) $$

Notice that $K_{eq}$ depends only on $\Delta_r G'^{\circ}$, which is the difference between the standard Gibbs free energy of formation of products and substrates involved in a reaction, and which can be calculated from tabulated values (where available). A good source of $K_{eq}$ values of many biochemical reactions is the eQuilibrator tool (equilibrator.weizmann.ac.il) [18, 19].

This thermodynamic treatment, showing that the equilibrium state of a reaction is captured by a constant relating to the ratios of product and substrate concentrations at that state, is fully supported by seminal experimental works from the second half of 1800s conducted on chemical reactions by Peter Waage (1833 - 1900) and Cato Guldberg (1836 - 1902), and their contemporaries. These works were concerned with the equilibrium, or steady-state, of chemical reactions attained under different conditions and when initiated from various starting concentrations of substrates. The key contribution of these studies was the finding that the equilibrium state in a reaction, that is the ratio of the concentration of substrates and products at steady-state, is characterized by a constant [20].

This finding, referred to as the "mass action law", later gave rise to the notion (rather erroneously) that reaction rate of a chemical reaction at constant temperature is 'proportional to the product of the concentrations of the reacting substances' [21]. This derived statement actually is not a law but presents a possible rate model that would be compatible with the experimentally observed equilibrium state (i.e. with the mass action law of equilibrium) [20, 21] (see Box 3.D and the Appendix 3.6).

### 3.2.2. Enzymes as catalysts of biochemical reactions

We mentioned many biochemical reactions to be catalyzed by enzymes. It is therefore worth briefly explaining enzymes. Enzymes are proteins, chains of amino acids, that fold in the cell in various 3D structures. For our purposes, we do not need to understand all the intricacies of how enzymes are made or how they fold into their structures (the reader is directed to excellent books on these subjects [22, 23]). Suffice to say that in their folded-state, enzymes can bind a set of target metabolites in such a way that puts these metabolites in a specific physio-chemical environment and physical orientation, where their specific biochemical reaction is facilitated. Thus, enzymes are catalysts that facilitate a chemical reaction among metabolites. As we will discuss further below, modeling of biochemical reactions catalyzed by enzymes requires developing a 'mechanistic' picture of how enzymes function. Such models can be developed based on numerous studies on enzyme structure and function. Here, we will only state that a generally accepted model involves enzymes binding their substrates - thereby forming an enzyme-substrate complex - and then transitioning to a state enabling catalysis. We can expand this model by also considering so-called allosteric binding
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

Thermodynamic interpretation

Gibbs free energy of reaction:

$$ \Delta_r G' = \Delta_r G'^{\circ} + R \cdot T \cdot \ln \frac{c_{eq}^{n_c} \cdot d_{eq}^{n_d}}{a^{n_a} \cdot b^{n_b}} $$

At equilibrium:

$$ \Delta_r G'^{\circ} = -R \cdot T \cdot \ln \frac{c_{eq}^{n_c} \cdot d_{eq}^{n_d}}{a^{n_a} \cdot b^{n_b}} \\ e^{\frac{-\Delta_r G'^{\circ}}{R \cdot T}} = \frac{c_{eq}^{n_c} \cdot d_{eq}^{n_d}}{a^{n_a} \cdot b^{n_b}} = K_{eq} $$

Kinetic interpretation

Backward reaction rate:

$$ k_{-} \cdot c^{n_{c}} \cdot d^{n_{d}} $$

Forward reaction rate:

$$ k_{+} \cdot a^{n_{a}} \cdot b^{n_{b}} $$

At equilibrium:

$$ k_{+} \cdot a_{eq}^{n_a} \cdot b_{eq}^{n_b} = k_{-} \cdot c_{eq}^{n_c} \cdot d_{eq}^{n_d} \\ \frac{k_{+}}{k_{-}} = \frac{c_{eq}^{n_c} \cdot d_{eq}^{n_d}}{a_{eq}^{n_a} \cdot b_{eq}^{n_b}} = K_{eq} $$

Cartoon representation of Gibbs free energy of reaction and the thermodynamic equilibrium – As a chemical reaction proceeds, the concentrations of substrates and products change, which in turn affects the 'energy in the chemical system'. We can, thus, capture the reaction advancement in a graph, where the x-axis represents the reaction advancement (i.e. the concentrations of substrates and products at different times in the reaction course) and the y-axis the internal energy of the system. The Gibbs free energy of reaction, in a way, indicates the position of the system in this graphical representation, where the thermodynamic equilibrium would be the energy minima. At equilibrium, reaction Gibbs free energy would be zero, allowing us to derive the relation between substrate and product concentrations at that point and their free energy of formation. This relation is known as the equilibrium constant of the reaction. The same relation can be derived using a rate model to describe the forward and backward reactions that make up the overall reaction. The thermodynamic result (or derivation) shows that a given reaction (under a given temperature) would always have the same substrate and product concentrations at equilibrium, a point that is empirically verified by experiments and that is known as the “mass action law”. The rate-based interpretation of this thermodynamic result (or law) is known as the “mass action rate model” and assumes that rate of a given reaction is proportional to the concentrations of substrates and products to the power of their stoichiometry, and adjusted by a rate constant (shown as $k_+$ and $k_-$ above).

sites, where specific molecules (including sometimes the enzyme's own substrate or product) can bind and alter the kinetics of either enzyme-substrate binding or catalytic activity. These allosteric sites, thus, provide a mechanism for regulation of enzymatic reactions (Fig. 3.5).
<figure><img src="image_15.png" /><figcaption>Figure 3.5: Enzymes and flux regulation – (A) Schematic representation of a biochemical reaction, highlighting the involvement of a catalyzing enzyme. For such enzyme-catalyzed reactions, the flux has an upper limit relating to total enzyme concentration and kinetic parameters of the enzyme (see section 3.2.3 and Appendix 3.7 for enzyme catalyzed reaction rate models). (B) Cartoon representation of enzyme structure and possible mechanisms of allosteric or competitive regulation. Such regulation can emerge either by the substrate of the enzyme or other metabolites binding the enzyme (left), and altering its overall reaction rate (either through competition with the substrate or by altering the enzyme structure and affecting its kinetic parameters, right).</figcaption></figure>

### 3.2.3. Modeling reaction fluxes - reaction rate models

Metabolic reactions can involve diverse biophysical mechanisms (uncatalyzed, enzyme-catalyzed, etc.) and can take place under diverse biophysical conditions inside a cell (membrane-bound, cytosolic, extracellular, coupled across membranes, etc.). As such, mechanistically complete, biophysical representation of all metabolic reactions in dynamic, mathematical models might never be possible [24]. Dynamical models of metabolic systems, as with all mathematical models, must therefore balance abstraction of real mechanistic features of a system with achieving a still useful and insight-providing model. At the core of all dynamical metabolic models are rate laws that aim to capture the kinetics of biochemical reactions.

**Non-enzymatic reactions - the reversible and irreversible mass action rate models** All rate models used in metabolic modeling are based on the so-called 'mass action law' described in Box 3.D above. As discussed in that section, the "mass action law", which is derived from thermodynamic principles, is compatible with a rate model that assumes reaction rate of a chemical reaction at constant temperature to be 'proportional to the product of the concentrations of the reacting substances' [21, 20] (see Box 3.D). This 'mass action rate model' is commonly used, especially in the context of elementary reactions (i.e. reactions involving one single step), and has been shown empirically to apply in the case of some non-elementary reactions [20]. According to the mass action model, the net rate of any reaction of the form given in Eq. (3.1) is given by;

$$v = k_{+} \cdot a^{n_{a}} \cdot b^{n_{b}} - k_{-} \cdot c^{n_{c}} \cdot d^{n_{d}}, \qquad (3.5)$$

where small letters denote concentration of the relevant species of the same letter, $n_i$ denote the stoichiometric coefficient for species i (as introduced above), and $k_+$ and $k_-$ denote kinetic rate constants relating substrate concentrations to reaction rate.

The mass action rate expression is such that if the first term is larger than the second then $v > 0$, and more reactant will convert to product than product converting to reactant (Box 3.D). This situation will continue until some point, where the second term will be larger than the first, and the opposite will occur. Consequently, this expression makes the system converge towards an equilibrium point, or steady-state, where $v = 0$. As long as the reagents are free to move, they will collide and interconvert (in both directions) at the microscopic level, even when the equilibrium is reached. However, at equilibrium, the amount of reactant converting to product equals the amount of product converting to reactant per unit of time, therefore there is no net consumption and production of metabolites (Box 3.D). When we have the concentrations that lead to the thermodynamic equilibrium of the reaction,
i.e. equilibrium concentrations, we will have;

$$v = 0 = k_{+} \cdot a^{n_{a}} \cdot b^{n_{b}} - k_{-} \cdot c^{n_{c}} \cdot d^{n_{d}}$$

$$\frac{k_{+}}{k_{-}} = \frac{c^{n_{c}} \cdot d^{n_{d}}}{a^{n_{a}} \cdot b^{n_{b}}}$$

This ratio is known as the reaction's equilibrium constant $K_{eq}$ and hence the 'mass action rate model' is consistent with the empirical observations of Waage and Guldberg. As we have shown in Eq. (3.4) above, the equilibrium constant is equivalent to the reaction's Gibbs free energy under standard conditions. Note that when considering a biochemical system (rather than a chemical one), it is customary to report Gibbs free energies for standard conditions adjusted for a pH of 7, and denoted with superscript o'. Thus, we can write;

$$\frac{k_{+}}{k_{-}} = K_{eq} = e^{-\frac{\Delta_{r}G^{\prime o}}{R \cdot T}} \qquad (3.6)$$

where $\Delta_{r}G^{\prime o}$ is the Gibbs free energy under biological standard conditions, and R and T denote the molar gas constant<fnref n="1" /> and temperature (in Kelvin) respectively (see Box 3.D). It is important to note here that, given $K_{eq}$ is a constant determined by thermodynamics, the parameters $k_{+}$ and $k_{-}$ cannot be chosen independently, i.e. $k_{-} = K_{eq}/k_{+}$.

Following on from this last point, it is important to consider a reaction with large $K_{eq}$, i.e. a reaction for which $\Delta_{r}G^{\prime o}$ is highly negative. In this case, the value of $k_{-}$ can become small to the extent that the reverse reaction can be negligible. In this case the reaction could be considered as effectively irreversible and the rate model can be approximated by;

$$v = k_{+} \cdot a^{n_{a}} \cdot b^{n_{b}} \qquad (3.7)$$

**Enzymatic reactions** The mass action rate discussed above forms also the basis of modeling enzymatic reactions. This approach is justified by considering each enzymatic reaction as a series of 'elementary steps', each obeying the mass action rate model. To this end, many alternative elementary steps, or 'enzyme mechanisms', can be considered to 'capture' an enzymatic reaction and subsequently many alternative assumptions can be made to simplify the resulting system of steps. It is also possible to include allosteric regulation or other types of inhibition or activation steps within these elementary steps, allowing generation of a rich variety of enzymatic models and rate laws. Here, we will cover some of the most common of such models, noticing that the construction of these models follows the same general principles of (i) drawing up elementary reactions, (ii) writing down mass action based kinetic rates for the system, and (iii) simplifying the system with assumptions on kinetic parameters (see Appendix 3.6). The reader can consult additional books (e.g. [25]) for more specific, elaborate enzymatic reaction schemes, or can attempt them as an exercise.

**Single substrate, irreversible enzymatic rate model (Michaelis-Menten model)** A possible representation of an enzyme mediated reaction consisting in the conversion of a reactant S to a product P could be the following reaction scheme:

$$S + E \underset{k_2}{\overset{k_1}{\rightleftharpoons}} ES \underset{k_{cat}}{\overset{k_{cat}}{\rightleftharpoons}} P + E.$$

This reaction scheme is rather specific, for example, it ignores the possibility that substrate bound enzyme can be converted into product, while remaining bound on the enzyme. Thus, the above reaction scheme is derived from a more complete and more complex reaction scheme through application of several assumptions relating to individual reactions. The resulting rate model from the above scheme is usually known as the Michaelis-Menten model, named after the biochemists Leonor Michaelis and Maud Menten who studied enzyme kinetics in the early 1900's, but several

<footnote n="1">¹The molar gas constant (also known simply as the gas constant) is the molar equivalent to the Boltzmann constant, expressed in units of energy per temperature increment per amount of substance (quantified in moles rather than single particles). Its value is about 8.31 J·K⁻¹·mol⁻¹.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 3.6: Michaelis-Menten rate law – The x- and y-axis show the substrate concentration (normalized by $K_M$) and reaction flux (normalized by $v^{max}$), i.e. $e_{tot} \cdot k_{cat}$.</figcaption></figure>

studies of that time and afterwards arrived at a similar model using different assumptions. Implementation of the specific assumptions, as we detailed in Appendix 3.7, allows one to arrive at the above reaction system, which can be represented by a reduced ODE system, compared to the full system. In this reduced ODE system, the ODE describing the rate of formation of the product, which is equivalent to reaction rate, becomes:

$$v = \frac{s \cdot e_{tot} \cdot k_{cat}}{K_M + s} \qquad (3.8)$$

where $e_{tot}$ represents the total enzyme concentration, $k_{cat}$ is known as the catalytic rate of an enzyme, and $K_M$ is known as the Michaelis-Menten coefficient of the enzyme and is equal to $(k_2 + k_{cat})/k_1$ (we note that depending on the assumptions used, the expression for $K_M$ can vary). Plotting the above rate of formation of product against increasing substrate concentration (see Figure 3.6) shows that the rate is a ‘saturating function’ of substrate, i.e. the rate approaches a threshold point - given by $v^{max} = e_{tot} \cdot k_{cat}$ as substrate concentration increases. Thus, we can see that the enzymatic nature of the reaction introduces a limiting factor on the reaction rate that depends on $v^{max}$, i.e. total enzyme concentration and enzyme’s catalytic rate. This fact underpins the regulation of metabolic flux through regulation of enzyme levels or enzyme’s catalytic rate, and is a key conceptual point for the constraint-based methods discussed later in this book.

**Single substrate, reversible enzymatic rate model (Haldane model)** Considering that all chemical reactions are — at least, in theory — reversible, it is also possible to express the rate of an enzyme-mediated reaction as a function of the concentration of both substrate and product. A method to do so has been introduced by Haldane [26]. It considers the following reaction scheme:

$$S + E \underset{k_2}{\overset{k_1}{\rightleftharpoons}} ES \underset{k_4}{\overset{k_3}{\rightleftharpoons}} EP \underset{k_6}{\overset{k_5}{\rightleftharpoons}} P + E.$$

Deriving the rate law for this reaction scheme is slightly more involved, but it follows the same strategy as explained above, of creating elementary steps, treating them as obeying mass action rate, and making additional simplifying assumptions. As shown in Appendix 3.6, we can follow this strategy to derive the reversible rate law as follows:

$$v = e_{tot} \cdot \frac{k_{cat}^{+}}{K_{S}} \cdot \frac{s - p \cdot \frac{k_{cat}^{-}}{K_{P}}}{1 + \frac{p}{K_{P}} + \frac{s}{K_{S}}} \qquad (3.9)$$

where $K_S$ and $K_P$ are composite constants relating to the substrate and product binding to the enzyme, and $k_{cat}^+$ and $k_{cat}^-$ are Haldane coefficients (again, composite parameters of other kinetic constants) describing catalytic rate
of the enzyme (see Appendix 3.6 for further details of these parameters).

As done in the above section on kinetics of the non-enzymatic reversible reaction, we can consider the equilibrium condition for this enzymatic reversible reaction. This would allow us to derive the corresponding relation between $K_{eq}$ and reaction Gibbs free energy. Recognizing the relation between the Haldane composite parameters and $K_{eq}$ (see Appendix 3.6) and the flux-force relation (see below), we can then re-formulate the reversible rate law as:

$$v = e_{\text{tot}} \cdot k_{\text{cat}}^{+} \cdot \frac{s/K_{\text{S}}}{1 + p/K_{\text{P}} + s/K_{\text{S}}} \cdot \left(1 - e^{\frac{\Delta_{\text{r}}G'}{R \cdot T}}\right) \quad (3.10)$$

where $\Delta_{\text{r}}G'$ is the Gibbs free energy of reaction for a given substrate and product levels under biological conditions and considering the forward direction of the reaction. This rate law shows that forward reaction rate will be independent of thermodynamics, when the reaction free energy is highly negative (i.e. when the reaction is far from thermodynamic equilibrium, $\Delta_{\text{r}}G' \ll 0$). However, as the reaction Gibbs free energy gets close to zero, the reaction rate will decrease, and as such, there will be a dependence of reaction rate on reaction free energy.

Another way of writing equation 3.10 is this one:

$$v = e_{\text{tot}} \cdot k_{\text{cat}}^{+} \cdot \frac{s/K_{\text{S}} \cdot \left(1 - e^{\frac{\Delta_{\text{r}}G'}{R \cdot T}}\right)}{1 + s/K_{\text{S}} \cdot \left(1 + \frac{k_{\text{cat}}^{+}}{k_{\text{cat}}^{-}} \cdot e^{\frac{\Delta_{\text{r}}G'}{R \cdot T}}\right)} \quad (3.11)$$

where we replace $p/K_{\text{P}}$ with an expression that depends on $s$ and $\Delta_{\text{r}}G'$. This alternative expression, developed in the context of modeling microbial metabolism [27, 28], can be useful because it shows us that when the reaction is far from equilibrium ($\Delta_{\text{r}}G' \ll 0$), the term $e^{\Delta_{\text{r}}G'/(R \cdot T)}$ will approach zero and the above formula can be approximated by the irreversible Michaelis-Menten rate law (Equation 3.8). In this case, we further notice that the Haldane coefficient $K_{\text{S}}$ becomes equivalent to $K_{\text{M}}$ introduced above in the irreversible reaction scheme (see section 3.2.3).

It is important to note that many reactions within cell metabolism are experimentally shown to be reversible, indicating that they operate close to thermodynamic equilibrium [29, 30, 18].

**Rate models for representing allosteric effects** Rate models for representing allosteric effects, i.e. binding of additional molecules - or their own substrates - on the enzyme and affecting the enzyme-mediated reaction rate, can be created either by adjusting the rate laws given above empirically, or by considering the additional binding events at ‘allosteric sites’ of the enzyme and deriving a new ‘mechanistic’ rate model. To give an example of the former strategy, we can consider a Michaelis-Menten rate model adjusted for an inhibitory effect of the substrate on the enzymatic reaction rate. This adjusted rate model can be expressed as:

$$v = \frac{\nu^{\text{max}} \cdot s}{K_{\text{M}} + s + s^2/K_{\text{I}}} \quad (3.12)$$

where $K_{\text{I}}$ represents the saturation coefficient for the binding of the substrate at an allosteric site on the enzyme. Notice that we used such a model in the small multi-stable system example introduced above (section 3.3.3) and discussed in Appendix 3.7.

For the same example, the alternative approach (the latter case mentioned above) would be to develop a mechanistic model involving multiple binding reaction on an enzyme. The resulting elementary reactions and their mass action implementation can be then carried out. This process would result in a set of ODEs, which can then be further simplified to draw a rate model for the proposed allosteric regulation. An example of this type model is developed in the context of multi-substrate binding enzymes, and shown to lead to multi-stability under certain parameter conditions [31].
<figure><img src="image_17.png" /><figcaption>Figure 3.7: The ratio of net forward flux (J) to forward reaction rate ($v_+$) as a function of the negative reaction Gibbs free energy</figcaption></figure>

**Flux-force relationship** All chemical reactions, including biochemical reactions, must obey thermodynamic laws. This fact manifests itself in several ways in dynamical modeling. Firstly, reaction direction (or, rather, feasibility) is determined by the sign of the reaction Gibbs free energy. Second, the kinetic constants associated with the elemental reaction steps are constrained by thermodynamics (section 3.2.3). To see the third relation arising from thermodynamics, we consider again the simple non-enzymatic mass action model we used above – reaction schematic given in Eq. (3.1) and the reaction Gibbs free energy given by Eq. (3.2).

We now re-consider the net rate of reaction as given above in Eq. (3.5), and break this into its components of forward reaction rate (or flux) and reverse reaction rate (or flux), which are given by;

$$v_+ = k_+ \cdot a^{n_a} \cdot b^{n_b}$$

$$v_- = k_- \cdot c^{n_c} \cdot d^{n_d}$$

and then, we can express the net forward flux (J) as:

$$J = v_+ - v_- = v_+ \cdot \left(1 - \frac{v_-}{v_+}\right) = v_+ \cdot \left(1 - \frac{k_- \cdot c^{n_c} \cdot d^{n_d}}{k_+ \cdot a^{n_a} \cdot b^{n_b}}\right) = v_+ \cdot \left(1 - \frac{k_-}{k_+} \cdot \Gamma\right)$$

In this re-organized form of the net forward flux, we notice that the expression in parentheses on the right hand side can be re-expressed in terms of reaction free energy (using Eq. (3.6)) as follows:

$$J = v_+ \cdot \left(1 - \frac{k_-}{k_+} \cdot \Gamma\right) = v_+ \cdot \left(1 - \frac{\Gamma}{K_{eq}}\right) = v_+ \cdot \left(1 - e^{\frac{\Delta_r G'}{R \cdot T}}\right)$$

Thus, we find that the net forward flux of the reaction is given by the forward reaction rate multiplied by a thermodynamic factor. When the reaction is energetically favored, i.e. has large negative Gibbs free energy, the thermodynamic factor diminishes and the net forward flux is fully determined by forward reaction rate alone (see Figure 3.7). When the reaction is closer to equilibrium, i.e. small negative or near-zero Gibbs free energy, then the net forward flux will be determined by a combination of forward and reverse flux rates. This relation between net forward flux and thermodynamics is referred to as the flux-force relation [32, 33] and holds also for the enzymatic reversible reaction model described above (see section 3.2.3).

**A note on choosing a reaction rate model** In the above sections, we have introduced several biochemical reaction rate models. These models fall into two main categories, namely those that model enzyme action (i.e. enzymatic models) and those that ignore the enzyme action (i.e. non-enzymatic models). Notice that derivation of both categories of models rely on the mass action law. In the non-enzymatic case, we model reactions as single-step
forward and backward reactions using mass action, while in the enzymatic case, we consider multi-step reaction mechanisms, but still use the mass action for each individual step. For each category, we can consider the reaction thermodynamics and model reactions as reversible, but as we discussed above we can also choose to approximate reactions as ‘irreversible’ when the overall reaction’s Gibbs free energy is very negative (i.e. when $K_{eq}$ is large).

In a given modeling context and metabolic system, it would be a valid question to ask – which model should one use? This question can be answered in parts. In the first instance, we can make a decision about the use of reversible or irreversible rate models. As already mentioned, this decision should be based on the value of $K_{eq}$ – a reaction with a very large $K_{eq}$ can be modeled as irreversible, as long as the product concentrations are known not to reach very high levels (in a cell). However, to represent a metabolic reaction as irreversible is not without consequences even if the reaction always runs in the same direction (notice that the assumption of irreversible reaction means that the reaction rate cannot go negative). Reversible kinetics can capture the negative feedback of reaction products on reaction rate, and irreversible reaction models would lose this feature [34]. A recent study by Shen et al [35] showed how important it can be to include product inhibition to create a predictive metabolic model.

In the case of lower $K_{eq}$ value – in combination with a consideration of possible product concentration – the modeler should opt for the reversible rate models, which are thermodynamically consistent. The decision about use of enzymatic or non-enzymatic reaction models can be made in a practical manner. If the enzyme associated with the modeled reaction has measured kinetic rates, it would be sensible to opt for a enzymatic model (noting that in vivo enzyme kinetics might differ from those measured in vitro and that many enzyme kinetics studies use parameter derivations assuming an irreversible Michaelis-Menten model). Consequently, it may not be possible to find all the required parameters in the literature, so to model a reaction using reversible rate model. In the absence of measured enzyme parameters, the modeler can use ‘guesstimated’ parameters, based – for example – on the distribution of known enzyme kinetic parameters, or alternatively, use the non-enzymatic model.

Given the discussion in the preceding paragraph, it is a useful exercise to consider when the non-enzymatic and enzymatic models might behave in the same way. We have introduced above the concept of flux-force relationship, where we have shown that the net flux in a reversible reaction would be given by the forward flux multiplied by a thermodynamic factor:

$$J = v_{+} \cdot \left( 1 - \frac{\Gamma}{K_{eq}} \right)$$

If we consider this equation for the reversible non-enzymatic and enzymatic models, we would notice that the thermodynamic factor would show the same behavior for both models, depending only on reaction $K_{eq}$ value and substrate and product concentrations. Where the models would differ, would be in the behavior of the $v_+$ term, which takes the form:

For the reversible enzymatic case:

$$v_{+} = e_{tot} \cdot k_{cat}^{+} \cdot (s/K_{S})/(1 + s/K_{S} + p/K_{P})$$

And, for the reversible non-enzymatic case:

$$v_{+} = s \cdot k_{+}$$

Where $k_{cat}$, $K_S$, and $K_P$ are the enzyme kinetic parameters for the enzymatic model and $k_+$ is the forward reaction rate coefficient for the non-enzymatic model. Thus, the two models would behave in a similar way, when there is correspondence between these two terms, which are sometimes referred to as “saturation terms” [33]. By re-arranging
the above terms, we can show that correspondence between the two models can be expressed as:

$$e_{\text{tot}} \cdot k_{\text{cat}}^{+} \cdot (1/K_{\text{S}})/(1 + s/K_{\text{S}} + p/K_{\text{P}}) \approx k_{+}$$

We can see that in the regime, where $s \ll K_S$ and $p \ll K_P$, both models would behave in a linear fashion and their behavior would correspond exactly with the right choice of parameters (i.e. assuming $(e_{\text{tot}} \cdot k_{\text{cat}}^{+}/K_S) = k_+$). Outside of this regime, correspondence would be dependent on both parameters and concentration of S and P. One interesting case to consider is when total amount of S and P would be conserved, for example, with cycling reaction schemes. In this case, we can introduce a new parameter C to describe the total pool of the cycled metabolite (e.g. $C = S + P$) and the correspondence would be expressed as:

$$(e_{\text{tot}} \cdot k_{\text{cat}}^{+}/K_{\text{S}})/(1 + (s \cdot K_{\text{S}} - s \cdot K_{\text{P}})/(K_{\text{S}} \cdot K_{\text{P}}) + C/K_{\text{P}}) \approx k_{+}$$

Thus, in this case of the sum of substrate and product concentrations being conserved, we can have correspondence between the non-enzymatic and enzymatic models when S is small or when $K_S = K_P$.

## 3.3. Dynamics and regulation of metabolism

As explained so far in this chapter, cell metabolism involves biochemical reactions involving metabolites (and often catalyzed by enzymes). Thus, understanding metabolism involves studying the dynamics of this system, trying to predict how metabolite levels will go up or down, or settle to a steady state as cell physiology changes in response to external or internal processes (e.g. cells encountering glucose or undergoing division). Obtaining such understanding requires us to develop models of biochemical reaction systems and predict the ‘dynamics’ of those systems. In the previous section, we learned how to model one biochemical reaction. Now we will see how we can readily expand these models to capture multi-reaction systems. The ‘art’ of developing and analyzing dynamical models falls under the branch of mathematics known as calculus and nonlinear dynamics. Many introductory books to these subjects are available, but we find that two particularly useful ones are those by Silvanus Thompson on calculus [36] and by Steven Strogatz on nonlinear dynamics [37]. Here, we will not re-introduce these topics but focus solely on various reaction rate models for metabolic systems that have been developed based on ODEs. We will highlight relations between these models and reaction thermodynamics and explore their possible limitations and applications in different cases. There are also books that are solely dedicated to models of biochemical reaction kinetics and enzyme kinetics more broadly - the reader is advised to further explore the topic with the help of such books, particularly [22, 25, 38]

### 3.3.1. Stoichiometric matrix and ordinary differential equations

As mentioned above, metabolic systems consists of many reactions. When describing multiple reactions in a biochemical 'system', it is convenient to represent the stoichiometries of individual reactions in a compact form called the stoichiometric matrix, N. The rows and columns of this matrix corresponds to m species (i.e. the metabolites), and to n reactions, found in the system respectively:

$$N \text{ is a } m \times n \text{ matrix}$$

The intersection of a row and column in the matrix indicates whether the species represented by that row takes part in the particular reaction represented by that column, or not. The sign of the element determines whether there is a net loss or gain of substance, and the magnitude describes the relative quantity of substance taking part in the reaction. It is important to appreciate that the elements of the stoichiometry matrix *do not* concern themselves with the rate of reaction, and just indicate the quantities taking part in the reaction.

A full description of a biochemical network, including the time-varying, dynamical behavior of metabolite concentra-
(A) Thermodynamic steady state

<figure><img src="image_18.png" /></figure>

(B) Dynamic steady states – non-equilibrium thermodynamics

<figure><img src="image_19.png" /><figcaption>Figure 3.8: Illustration of thermodynamic equilibrium and dynamical steady state – (A) Thermodynamic steady state. (B) Dynamic steady states – non-equilibrium thermodynamics. While the former happens only at chemical equilibrium, the latter can arise in systems that are far from chemical equilibrium. A cartoon of a flowing water through a tank and a reaction involving co-substrate cycling are shown as examples of systems that can attain dynamical steady states.</figcaption></figure>

tions, will augment the stoichiometry matrix with a rate vector, v, forming a so-called system equation:

$$ \frac{ds}{dt} = Nv(s) \qquad (3.13) $$

This equation represents a system of ordinary differential equations (ODEs) that describe the time evolution of the species, s. In other words, the ODE for species s describes the rate of change in the concentration of s with a given (infinitesimal) change in time. The ODEs can be solved numerically (i.e. simulated) by computer or studied analytically.

Notice that in mathematics, the time varying entities in a dynamical systems - in our context, the concentrations of chemical species - are known as 'variables', while any elements of the system that stay constant over time are known as 'parameters'. For an insightful and accessible mathematical treatment of differential equations and system dynamics, the reader is referred to these two excellent books [36, 37], while for a metabolic view of variables and parameters, the article on the Control of Flux, by Kacser and Burns, offers a valuable perspective [39].

### 3.3.2. Dynamic steady state

As stated above, the ODEs describe the time evolution of all variables s in the system. An informative approach to any dynamical system is to consider its steady state, a state where consuming and generating processes on each variable would have the same rate, i.e. the ODEs are equal to zero, and there would be no change in the variable amounts. For example, a water tank filling at a constant rate but emptying at a rate proportional to the height of water in the tank will eventually reach a steady-state where the output flow equals the inflow of water (Fig. 3.8). Under these conditions the height of water remains constant, or at a steady state.
<figure><img src="image_20.png" /><figcaption>Figure 3.9: Cartoon of a simple pathway that features allosteric enzyme regulation and that can show multiple steady-state solutions (see Appendix 3.6). The metabolite ‘x’ positively regulates the first step, $v_1$. The resulting positive feedback can result in a bistable system under a certain parameter regime.</figcaption></figure>

It is important to note that the thermodynamic equilibrium mentioned above is also a type of steady-state, but this does not mean that steady-state is only attained at thermodynamic equilibrium. In other words, there can be a steady-state where the system is out of thermodynamic equilibrium but the concentrations of metabolites are not changing. An example of this would be a linear metabolic pathway of connected reactions, with influx and outflux of an initial and endpoint metabolite (as seen in Fig. 3.8). In such a system, we can readily consider a scenario where there is influx of the first metabolite, outflux of the last metabolite, and forward flux through each of the reactions in the pathway. Thus, we would have a situation where all reactions are out of thermodynamic equilibrium, but all metabolite concentrations in the pathway attain a dynamic steady-state, where their influx and outflux are equal (Fig. 3.8). The distinction between systems that are both at steady-state and thermodynamic equilibrium, and those that are at steady-state but out of thermodynamic equilibrium, is an important one. It has been shown that complex dynamics, such as bistability and oscillations (as discussed below) are only possible in the latter case [40, 32, 41].

Mathematically speaking, the steady-state is defined when the ODE system, i.e. the system equation, is set to zero:

$$ \frac{ds}{dt} = N v(s) = 0 \qquad (3.14) $$

For simple systems, such as a tank of water filling and emptying, there is only one unique steady-state. This is perhaps better illustrated with a simple biochemical example. Consider a two step pathway where the first step has a constant rate $k_1$ and the second step a variable rate determined by a first-order reaction rate, $k_2$.

$$ \xrightarrow{v_0 = k_1} S \xrightarrow{v_1 = k_2 \cdot s} \qquad (3.15) $$

The differential equation describing the concentration of metabolite S is given by:

$$ \frac{ds}{dt} = k_1 - k_2 \cdot s \qquad (3.16) $$

Setting this equation to zero and solving for s yields the steady-state level of S:

$$ s = \frac{k_1}{k_2} \qquad (3.17) $$

This solution indicates there is only a single steady-state for this system dependent on the parameters $k_1$ and $k_2$.

### 3.3.3. Multiple steady-states and oscillations

In the previous section it was shown that a simple two step pathway admitted a single steady-state. There can be, however, metabolic systems that can show multiple steady states. As a simple example, consider the system shown in Figure 3.9. This shows a linear pathway of two reactions, with the first reaction activated by the species x.

Under certain parameter and model choices, such a system can admit three steady-states. Details of a model that can be simulated can be found in Appendix 3.7). Other examples of metabolic systems with multiple steady-states will be given below. In bi-, or multi-stable systems, there can be multiple sets of steady state concentrations and flux rates that the system can settle at. Which set of steady-states is realized is usually determined by initial concentrations or can be caused by a change in one of the concentrations or parameters. Thus, the system can change its steady-state
value abruptly at a threshold value of a specific parameter of the system. For a metabolic system displaying bistability, we can expect a rapid switch in multiple fluxes with changes in the concentration of one or few metabolites [37]. Furthermore, when bistability is combined with noise in some parameters (e.g. enzyme expression level) there can be a multi-modal distribution of flux states across genetically identical cells (e.g. see [42, 43] and section 3.5).

## 3.4. Concluding remarks

In this chapter, we set out to introduce cellular metabolism as a dynamical system. We have seen that metabolism comprises many biochemical reactions, that are historically cataloged and described into pathways. These pathways are usually not linear, composing of serial conversions of metabolites, but rather display branching points and interconnections through metabolites participating in many reactions. This inter-connected nature of metabolic systems, together with the large numbers of participating metabolites and reactions, makes them a complex system to study and conceptualize.

We have introduced both simplified, coarse-grained viewpoints for describing metabolism, and mechanistic approaches for detailed dynamical modeling of it at the level of single reactions. The former can be used to guide specific ideas on how to study metabolism, or to develop analogies to other disciplines, while the latter can provide a toolbox for constructing dynamical models of small or large metabolic systems. We have provided specific examples of such dynamical models and shown how they can allow us to relate system behavior - steady state or temporal behavior - to specific reaction mechanisms or parameters (e.g. allosteric interactions between metabolites and enzymes, cyclic reaction schemes, branching points).

There are many challenges remaining in the analysis and understanding of metabolism as a dynamical system. Recent studies found for example that many fluxes, where measured, are lower than predicted from a enzymatic irreversible reaction rate model (introduced in Eq. (3.8)) [44], and changes in flux patterns with changing conditions cannot be explained by enzyme levels alone [45]. These findings lead to the question on what determines/limits reaction fluxes and how reaction fluxes are regulated besides regulation of enzyme levels. There are several possible answers, including effects relating to allosteric interactions between metabolites and enzymes, reaction thermodynamics, and substrate-related effects. The experimental study and model incorporation of these possibilities is ongoing in systems biology, with increasing interest to include also more of the physico-chemical aspects of the cellular environment into the study of metabolism - such as diffusion of molecules, involvement of radical chemistry (especially generation of oxygen radicals in respiration) and membrane potential [46, 24]. As such, we are increasingly hoping to move from metabolic reactions studied in isolation, to cell-scale models and physico-chemical concepts that unite cell metabolism and physiology. Some of this emerging movement is captured in subsequent chapters of this book.

## Recommended readings

### Enzyme kinetics and reaction models

* "Enzymes" by J. B. S. Haldane [26]. Historically important book on enzyme kinetics and enzymatic reaction models.

* "Fundamentals of Enzyme Kinetics" by A. Cornish-Bowden [25]. General introductory book on enzymes and enzyme catalysis.

* "Enzyme Kinetics for Systems Biology" (2012) by H. M Sauro [38]. In addition to covering enzyme kinetics, this book also discusses stochastic kinetics and the kinetics of gene regulatory systems with an emphasis on systems biology models.

* "Structure and Mechanism in Protein Science: Guide to Enzyme Catalysis and Protein Folding" by A. Fersht [22]. General introductory book on enzymes and enzyme catalysis.

### Thermodynamics and physical chemistry
* "Understanding thermodynamics" by H. C. van Ness [17]. An excellent book that de-mystifies thermodynamics. It provides a conceptual treatise, leaving the mathematics to the side and focusing on what actually the thermodynamic laws mean.

* "Principles and Problems in Physical Chemistry for Biochemists" by N. C. Price [16]. An introductory book on thermodynamics, physical chemistry, and biochemistry.

## Metabolic system dynamics

* "Energy metabolism of the cell: a theoretical treatise" by J. G. Reich and E. E. Sel'kov [11]. Provides an early view of the importance of reaction dynamics as a 'self-regulatory' element in metabolism. Emphasizes the importance of cyclic reaction schemes and interconnections among metabolic processes.

* "Chemical Biophysics: Quantitative Analysis of Cellular Systems" by D. A. Beard and H. Qian [46]. Provides a rare approach of attempting to combine co-study the more physical aspects of cell physiology, including membrane potential and compartmentalization, with metabolism dynamics.

* "Systems Biology: An Introduction to Metabolic Control Analysis" (2018) by H. M Sauro [47]. Discusses biochemical network dynamics from the perspective of metabolic control analysis.

## Problems

Computer exercises for this chapter can be found on the book website.

### Problem 3.1 An irreversible reaction with simultaneous binding

(a) Write the reaction scheme for an irreversible enzymatic reaction with two substrates. Assume both substrates bind the enzyme simultaneously (forming one complex $ES_1S_2$), and both products are released simultaneously from this complex (i.e. without intermediary $EP_1P_2$ stage).

(b) Find the rate of product production for this system.

### Problem 3.2 A reversible reaction

(a) Write the reaction scheme for a reversible enzymatic reaction with two substrates. Assume both substrates bind the enzyme simultaneously (forming one complex $ES_1S_2$), and both products are released/absorbed simultaneously from/into this complex (i.e. without intermediary $EP_1P_2$ stage).

(b) Find the rate of product production for this system.

### Problem 3.3 An irreversible reaction with sequential binding

(a) Write the reaction scheme for an irreversible enzymatic reaction with two substrates. Assume the substrates bind sequentially (forming complexes $ES_1$ and $ES_1S_2$), and both products are released simultaneously from $ES_1S_2$ (i.e. without intermediary $EP_1P_2$ stage).

(b) Find the rate of product production for this system.

### Problem 3.4 An irreversible reaction with random-order binding

(a) Write the reaction scheme for an irreversible enzymatic reaction with two substrates. Assume the substrates bind the enzyme in any order (forming complexes $ES_1$, $ES_2$ and $ES_1S_2$), and both products are released simultaneously from this $ES_1S_2$ (i.e. without intermediary $EP_1P_2$ stage).

(b) Find the rate of product production for this system. Note that symbolic math tools such as MATHEMATICA, MAPLE or the SYMPY Python library will be helpful for this question (though not essential).

## Appendix sections
<figure><img src="image_21.png" /></figure>

Figure 3.10: Respiration-fermentation switch – Cartoon representation of the respiration and fermentation pathways. There is an observed switch of metabolic fluxes between these pathways, from full respiration to fermentation or a combination of respiration and fermentation, as glucose levels (and consequently growth rate) increases. It is postulated that this relates to a limitation in the respiration and ETC system, but the molecular basis of the switch is not fully clear with several, equally plausible hypotheses postulated.

**3.5. Dynamics of metabolism: experimental evidence and model-based explanations**

The high-level of connectivity among reactions, together with the plurality of molecular level mechanisms that can arise in enzyme-mediated reactions, gives metabolic systems the capacity to display rich dynamic behaviors [11, 41]. Here, we highlight some of the illustrative experimental observations on these metabolic dynamics, and their possible model-based explanations.

**3.5.1. Flux switching / regulation**

We have introduced above the redox-based, electron flow view of metabolism. A common electron donor in metabolism is glucose, while a common electron acceptor is oxygen together with the associated, membrane-bound electron transport chain (ETC). The ensuing metabolic pathway linking glucose oxidation to oxygen reduction is termed as respiration, resulting in formation of CO₂ (from full oxidation of glucose) and water (from reduction of oxygen, cf. Fig. ??). However, it is possible for cell metabolism to stop the sequential oxidation of glucose (or other sugars) at an intermediate level. In this case, the ensuing metabolism is termed as 'fermentation due to production of partially oxidized carbon compounds such as acetate and ethanol (Fig. 3.10) [1].

One of the earliest observation on metabolic dynamics is a shift from pure respiration into fermentation or respiro-fermentation with changing conditions. This shift, known as contre-Pasteur, Warburg, or, Crabtree effect, is described initially in yeast and mammalian cells, especially cancerous cells. The respiration to fermentation shift happens under lack of electron acceptors or with increasing growth rate. While a shift into fermentative pathways due to lack of electron acceptors can be intuitively understood as the only route to sustain electron flow, a similar shift due to increased carbon availability or growth rate are non-intuitive as they occur under the continued presence of strong electron acceptors such as oxygen.

A dominant concept to explain the switch to respiro-fermentation has been the idea of 'overflow metabolism. It postulates that this switch should be seen as an overflow, arising due to limitations in respiration not being sufficient in sustaining metabolic fluxes in the face of increasing substrate availability. The dynamic regulation and origin of
this respiro-fermentation switch is still a focus of significant systems biology research. Hence, this topic is discussed further in other chapters of this book with several alternative models presented for its underlying causes.

It must also be noted that, while respiro-fermentation switch is commonly referred to as 'overflow metabolism' (due to excretion of fermentation products such as acetate, lactate and ethanol), the phenomenon of overflow, i.e. excretion of energy rich compounds is not limited to fermentation. Excretion of amino acids and vitamins seem particularly common, and it is not clear in these cases what type of metabolic flux switching happens or how it happens.

**Flux switching / regulation - flux sensors, branch point dynamics and dynamical flux regulation** How can we understand cells switching their metabolic fluxes with changing external or internal conditions. As discussed above, one possibility is that cells alter the expression levels of their various enzymes, so to achieve a re-distribution of fluxes (given that enzyme levels are directly involved in the determination of fluxes, see Eq. (3.8)). This kind of enzyme-level regulation can be mediated through regulation of transcription factors by specific signaling molecules, including metabolites. The latter case is explored in models of central metabolism, and it was shown that fructose-1,6-bisphosphate could act as a "flux sensor", conveying information about the glycolytic flux onto key transcription factors regulating glycolytic enzymes.

While regulation of enzyme levels can alter flux levels, this type of regulation can be made more sensitive if the coupling between enzyme and flux levels can be made more nonlinear. It has been shown that such nonlinear coupling of flux and enzyme levels can arise at branching points in metabolism [9]. In particular, branch points (as a metabolic motif) have been shown to give rise to ultrasensitivity - a system dynamics feature that describes the situation when a given input to a system results in more than a proportional change in its output. When enzymes at the two branches of a branching point have highly differing affinities for the substrate (i.e. different $K_M$ values, see Eq. (3.8)), then alterations of the maximal rate of one enzyme with higher affinity to the substrate can result in a nonlinear effect on the flux into the other branch of the branching point [9]. Thus, branching points can be one structural motif that can result in switch-like, nonlinear flux changes within metabolism.

It is clear that changing of enzyme levels can regulate fluxes, and can do so in an abrupt, switch-like fashion through structural motifs such as branching points. However, regulation of enzyme levels via transcription factors is found to not capture all observed flux changes in experiments. This suggests that cells might be able to regulate fluxes by other means as well. Recently, one such possible mechanisms is proposed to be the co-substrate pools [15]. For example, in a metabolic branch point, where the two branches involve different co-substrates, regulating the pool sizes of those two co-substrates can induce flux switching at the branch point.

### 3.5.2. Bistability

Bistability is introduced above, and refers to a dynamical system having three steady state, two of which are dynamically stable and can be attained by the system. When bistable systems exist in cell metabolism, their combination with population level variance (i.e. noise) in enzyme levels or activity can lead to bimodal distribution of metabolic fluxes (i.e phenotypes) in isogenic population of cells. In this context, it is notable that significant level of variance is seen in several metabolic parameters, including sugar uptake, ATP levels, and expression levels of the enzymes involved in glycolysis and the TCA cycle.

Bistability in metabolic responses is experimentally implicated in the context of respiration to fermentation switch, and when carbon metabolism is initiated on glucose or when it switches from glucose to other carbon sources. In particular, the latter studies found sub-populations, within isogenic populations (i.e. no mutations), that show different metabolic responses to changing conditions. Experiments with isotope labeled carbon indicated that these sub-populations emerged at the time of the shift in carbon source is induced, i.e. in response to changing conditions, and in a manner dependent on the concentration of the new carbon source. This suggests that the metabolic system implements bistable dynamics, such that changes in external glucose concentrations can lead some cells to shift to a new metabolic steady-state flux distribution, while others remain at their original steady state.
**Bistability - negative feedback via substrate inhibition** There have been many theoretical studies indicating the possibility of bistability within simple enzymatic reaction systems. Bistability is shown to be possible even in a single enzymatic reaction, involving allosteric regulation, or in a system of few coupled enzymatic reactions. A particular 'reaction motif' that has been studied extensively is a two-enzyme cyclic reaction system, where a substrate is converted into a product and then back again, with both forward and backward reactions usually involving different enzymes (see Fig. 3.11). It is common, in these models, that the enzyme catalyzing the forward reaction is assumed to be regulated by substrate inhibition, or by substrate inhibition coupled with product activation. This motif is found in several locations within metabolism, particularly around dehydrogenases, such as lactate dehydrogenase, and kinase/phosphatase pairs, such as those involved around fructose-6-phosphate, that can convert different metabolites back and forth, using the NAD+/NADH or ADP/ATP pairs as reaction partners.

These theoretical findings are supported by several *in vitro* re-constitution experiments that confirmed bistability experimentally using enzyme preparations of pyruvate kinase, lactate dehydrogenase, and isocitrate dehydrogenase enzymes and their corresponding partners resulting in cyclic reaction schemes.

It is notable that many of these models incorporated negative feedback via empirical alteration of a Michaelis-Menten type reaction rate model (i.e. one of the approaches we mentioned in the paragraph above on allosteric rate models, see 3.12). This raises the question about the actual biochemical mechanisms that can lead to bistability in a enzyme-mediated reaction model. In a recent study, it was shown that the presence of multiple enzyme-substrate complexes, as would be the case in an enzyme with multiple substrate binding sites, creates a potential in the reaction system for bistability [31] (see Fig. 3.11). Thus, multi-site enzymes could be points of multistability generation in metabolic systems and any larger models featuring such enzymes or inherently including feedback regulation can demonstrate bistability.

### 3.5.3. Oscillations

Sustained and damped oscillations are common dynamics in nonlinear systems and can arise from a combination of positive and negative feedbacks [37]. In metabolic systems both types of oscillations are seen *in vivo* or *in situ*, with cell extracts, where concentrations of all observed metabolites are found to oscillate over time. In the case of experiments involving cell extracts, these metabolite oscillations had a period ranging from few to tens of minutes. In these experiments, oscillations are verified not to be due to artificial changes in ATP dynamics arising in the cell extract preparations, and oscillations could be entrained by controlled glucose additions. This shows that there is an inherent ability for oscillatory dynamics in the underpinning enzymatic reaction system. This ability is suggested to be linked to the enzyme phosphofructokinase (PFK), which catalyzes the phosphorylation of fructose-6-phosphate into fructose-diphosphate in the glycolysis pathway and is allosterically regulated.

*In vivo*, oscillatory dynamics are observed to occur within the central carbon pathways and displaying a phase of tens of minutes up to several hours. Metabolic oscillations were demonstrated at single cell level and are found to be autonomous of, but coupled with, the cell cycle oscillations. Additional studies across cell populations found that cells can synchronize metabolic oscillations under some conditions, and proposed several possible mediators for such synchronization, including acetaldehyde, hydrogen sulphide, carbon dioxide, and media pH.

**Oscillations - intertwined negative and positive feedbacks** Several mathematical models of the reaction catalyzed by the enzyme phosphofructokinase (PFK) in the glycolysis pathway has shown that oscillations are possible to arise from the dynamics of this reaction alone. These models incorporate some of the observed allosteric regulation of PFK both by its substrates and products, resulting in intertwined negative and positive feedbacks.

It must be noted that some of these models, and others, use the same basic models that show bistable behavior (as discussed above) and extend them with in- and out-fluxes of involved metabolites, to display oscillations. While these theoretical demonstrations of specific enzymatic schemes leading to oscillations have not been explored in detail experimentally, metabolic oscillations are readily observed both *in vivo* and *in vitro*, as discussed above. Models,
(A) Allosteric enzyme model

<figure><img src="image_22.png" /></figure>

At steady state:

$$ \underbrace{\frac{V_1 \cdot [S]}{K_1 + [S] + [S]^2/K_3}}_{\text{P production}} = \underbrace{\frac{V_2 \cdot (C - [S])}{K_2 + (C - [S])} + \alpha \cdot (S_0 - [S])}_{\text{P consumption}} $$

(B) Multi-site enzyme model

<figure><img src="image_23.png" /></figure>

Figure 3.11: Cartoon representations and brief analysis results of two enzymatic models capable of bistability. (A) Allosteric enzyme model. The first model considers an enzyme that can convert a substrate (S) into a product (P) and that is allosterically regulated by its own substrate. This regulation takes the form of inhibition and is implemented mathematically in the rate of the enzyme - black colored equation. This model results in a nonlinear curve for the relation between rate of production of P at steady state and the total concentration of substrate and product in the system, $S_{tot}$ (black curve on the top right panel). The intersections of this curve with the linear curve for the relation between rate of consumption of P at steady state and $S_{tot}$ (red curve on top right panel). We can see that the model is capable of resulting in three intersections, i.e. three steady states of the system. (B) Multi-site enzyme model. The second model considers instead of allostery, an enzyme that binds multiple substrates. This results in several enzyme-substrate complexes depending on the number of binding sites - 3 sites in the model shown. The resulting model can be solved for the steady state values of flux through each enzyme complex against $S_{tot}$ (shown in red and blue colors on the bottom right panel). The sum of these gives the rate of production of P at steady state (black curve on the bottom right panel). This model can also result in a non-linear production curve and three steady states. For further discussion of these models, see relevant citations.

involving some of these proposed synchronization molecules, were also developed and could reproduce experimental findings.

## 3.6. Derivation of enzymatic rate laws

Enzymatic reactions can be modeled using a mechanistic model of enzyme binding and catalysis. The general approach is to develop a 'cartoon' model of the physical steps in a reaction. This cartoon model usually takes the form of a series of reactions, involving either binding / unbinding events or chemical conversions. Once a model is developed one can write down ordinary differential equations (ODEs) based on these reactions, and assuming each reaction to
be governed by mass action kinetics (see Section 3.2.3). The ODEs can be simplified using certain assumptions, or sometimes just kept as is, before applying a quasi steady-state assumption (which states the enzyme-substrate complexes to be in steady-state). This assumption would allow us to solve the ODE for the enzyme-substrate complex($e_s$) at steady-state. We then enter these solutions into the ODE for the product, so to obtain a reduced system and a specific rate law for product formation. This approach forms the basis of obtaining simplified rate laws, that is, a reduced ODE for the rate of product formation, for enzymatic reactions.

### 3.6.1. Derivation of the single substrate, irreversible rate law

This is the most generic model of an enzymatic reaction that has been developed/studied by Leonor Michaelis (1875 1947) and Maud Leonora Menten (1879 1960), and their contemporaries. It involves the following reaction scheme, where a substrate binds to an enzyme to form a enzyme-substrate complex, gets converted into a product, and then released from the enzyme:

$$S + E \underset{k_2}{\overset{k_1}{\rightleftharpoons}} ES \underset{k_4}{\overset{k_3}{\rightleftharpoons}} EP \underset{k_6}{\overset{k_5}{\rightleftharpoons}} P + E. \qquad (3.18)$$

We can simplify this reaction system by assuming that (1) the transition between enzyme complexes ES and EP are instantaneous and are therefore considered as a single entity, e.g. $e_s$, and (2) that the release of product and enzyme is irreversible. The scheme now becomes:

$$S + E \underset{k_2}{\overset{k_1}{\rightleftharpoons}} ES \underset{k_3}{\overset{k_3}{\rightleftharpoons}} P + E. \qquad (3.19)$$

We can now write a set of ODEs to describe the dynamics of this reaction system - using mass action kinetics. The ODEs are as follows:

$$\frac{ds}{dt} = -s \cdot e \cdot k_1 + e_s \cdot k_2$$

$$\frac{de}{dt} = -s \cdot e \cdot k_1 + e_s \cdot (k_2 + k_3)$$

$$\frac{dc}{dt} = s \cdot e \cdot k_1 - e_s \cdot (k_2 + k_3)$$

$$\frac{dp}{dt} = e_s \cdot k_3$$

where we used the small letter notation to represent the concentration of each species, e.g. "e" for the concentration of the enzyme, E, and "e_s" for the concentration of the enzyme-substrate complex, ES. At this stage, we can see that if we can formulate "e_s" as a function of "s", we can provide a simpler rate model that relates production of the product, P, to the level of the substrate, s. To achieve this we make several additional assumptions. First, we will assume that the total level of the enzyme is conserved, i.e. $e + e_s = C$, where C is a constant (referred to as $e_{tot}$ in the main text). This assumption effectively means that total enzyme levels are fixed in the timescale of reaction dynamics. This assumption already allows us to re-define the ODEs and reduce their number to three from four - since, we can now express e, as a function of $e_s$. The new ODEs look like this:

$$\frac{ds}{dt} = -s \cdot (C - e_s) \cdot k_1 + e_s \cdot k_2$$

$$\frac{de_s}{dt} = s \cdot (C - e_s) \cdot k_1 - e_s \cdot (k_2 + k_3)$$

$$\frac{dp}{dt} = e_s \cdot k_3$$

Second, we will assume that the binding/unbinding of substrate to the enzyme happens much faster than release of product from the enzyme-substrate complex. This assumption, together with the additional assumption that enzyme
levels are much lower than substrate levels, allows us to consider the enzyme-substrate complex to remain constant throughout the reaction. In other words, we consider the enzyme-substrate complex to be in a 'quasi steady-state'. This allows us to solve the second ODE from above for steady-state:

$$
\begin{align*}
\frac{\mathrm{d}e_s}{\mathrm{d}t} &= 0 = s \cdot (C - e_s) \cdot k_1 - e_s \cdot (k_2 + k_3) \\
e_s \cdot (k_2 + k_3) &= s \cdot (C - e_s) \cdot k_1 \\
e_s \cdot (k_2 + k_3) &= sC \cdot k_1 - s \cdot e_s \cdot k_1 \\
e_s \cdot (k_2 + k_3 + s \cdot k_1) &= s \cdot C \cdot k_1 \\
e_s &= \frac{s \cdot C \cdot k_1}{(k_2 + k_3 + s \cdot k_1)}
\end{align*}
$$

We have now an expression for “e_s”, which we can simply introduce to the ODE system. We have effectively reduced our ODE system from a three variable system into a two variable one:

$$
\begin{align*}
\frac{ds}{dt} &= -s \cdot \left(C - \frac{s \cdot C \cdot k_1}{(k_2 + k_3 + s \cdot k_1)}\right)k_1 + \frac{s \cdot C \cdot k_1}{(k_2 + k_3 + s \cdot k_1)} \cdot k_2 \\
\frac{dp}{dt} &= \frac{s \cdot C \cdot k_1}{(k_2 + k_3 + s \cdot k_1)} \cdot k_3
\end{align*}
$$

The second ODE describes the rate of change in product, P, as a function of substrate, S. It is a rate model for this enzymatic reaction, and holds under the assumptions we made in its derivation. It is known as the Michaelis-Menten kinetic rate model and is commonly expressed as:

$$
v = \frac{s \cdot e_{\text{tot}} \cdot k_{\text{cat}}}{K_{\text{M}} + s}
$$

where e<sub>tot</sub> is equal to C and represents total enzyme concentration, k<sub>cat</sub> is equal to k<sub>3</sub> and is known as the maximal catalytic rate of an enzyme, and K<sub>M</sub> is equal to (k<sub>2</sub> + k<sub>3</sub>)/k<sub>1</sub> and is known as the Michaelis-Menten coefficient of the enzyme. Plotting this rate against increasing substrate concentration would show that the rate is a 'saturating function' of s, i.e. the rate approaches a threshold point - given by v<sup>max</sup> = e<sub>tot</sub> · k<sub>3</sub> as substrate increases. The enzymatic nature of the reaction introduces a limiting factor on the reaction rate! This saddle point is actually a underpinning point for some of the constraint-based methods discussed in this book.

3.6.2. Derivation of a two substrate, irreversible rate law

See Problem 3.2

3.6.3. Derivation of the single substrate, reversible rate law

We now return to the reaction scheme we considered in the above section:

$$
S + E \underset{k_2}{\stackrel{k_1}{\rightleftharpoons}} ES \underset{k_4}{\stackrel{k_3}{\rightleftharpoons}} EP \underset{k_6}{\stackrel{k_5}{\rightleftharpoons}} P + E.
$$

The corresponding ODE system, written only for the key variables *e*<sub>s</sub>, *e*<sub>p</sub>, and *p*, is as follows:

$$
\begin{align*}
\frac{\mathrm{d}e_s}{\mathrm{d}t} &= e \cdot s \cdot k_1 + e_p \cdot k_4 - e_s \cdot (k_2 + k_3) \\
\frac{\mathrm{d}e_p}{\mathrm{d}t} &= e \cdot p \cdot k_6 + e_s \cdot k_3 - e_p \cdot (k_4 + k_5) \\
\frac{\mathrm{d}p}{\mathrm{d}t} &= e_p \cdot k_5 - e \cdot p \cdot k_6
\end{align*}
$$
As above, we will now introduce the assumptions of (1) total enzyme being conserved, and (2) the quasi steady-state, but this time for both of the enzyme-substrate and enzyme-product complexes. We will denote total enzyme concentration as C, as before, and use these two assumptions to express e_s and ep in terms of each other, and the other variables. Let us first proceed with e_s ;

$$ \frac{\mathrm{d}e_s}{\mathrm{d}t} = 0 = e \cdot s \cdot k_1 + ep \cdot k_4 - e_s \cdot (k_2 + k_3) $$

$$ e_s \cdot (k_2 + k_3) = (C - e_s - ep) \cdot s \cdot k_1 + ep \cdot k_4 $$

$$ e_s \cdot (k_2 + k_3 + s \cdot k_1) = (C - ep) \cdot s \cdot k_1 + ep \cdot k_4 $$

$$ e_s = \frac{C \cdot s \cdot k_1 + ep \cdot (k_4 - s \cdot k_1)}{(k_2 + k_3 + s \cdot k_1)} $$

We carry the same derivation for ep ;

$$ \frac{\mathrm{d}ep}{\mathrm{d}t} = 0 = e \cdot p \cdot k_6 + e_s \cdot k_3 - ep \cdot (k_4 + k_5) $$

$$ ep \cdot (k_4 + k_5) = (C - e_s - ep) \cdot p \cdot k_6 + e_s \cdot k_3 $$

$$ ep \cdot (k_4 + k_5 + p \cdot k_6) = (C - e_s) \cdot p \cdot k_6 + e_s \cdot k_3 $$

$$ ep = \frac{C \cdot p \cdot k_6 + e_s \cdot (k_3 - p \cdot k_6)}{(k_4 + k_5 + p \cdot k_6)} $$

We see that we have a symmetry in the expressions for e_s and ep, in that the two expressions can be derived from each other by a replacement of variables (k_1, k_4, k_2, s) → (k_6, k_3, k_5, p). Keeping this symmetry in mind, we now attempt to eliminate one of the complexes from the equation for the other:

$$ ep \cdot (k_4 + k_5 + p \cdot k_6) = C \cdot p \cdot k_6 + e_s \cdot (k_3 - p \cdot k_6) $$

$$ ep \cdot (k_4 + k_5 + p \cdot k_6) = C \cdot p \cdot k_6 + \frac{C \cdot s \cdot k_1 + ep \cdot (k_4 - s \cdot k_1)}{(k_2 + k_3 + s \cdot k_1)} \cdot (k_3 - p \cdot k_6) $$

$$ ep \cdot (k_4 + k_5 + p \cdot k_6) = C \cdot p \cdot k_6 + \frac{C \cdot s \cdot k_1 k_3 - C \cdot s \cdot k_1 \cdot p \cdot k_6 + ep \cdot (k_4 - s \cdot k_1) \cdot (k_3 - p \cdot k_6)}{(k_2 + k_3 + s \cdot k_1)} $$

$$ ep \cdot (k_4 + k_5 + p \cdot k_6) \cdot (k_2 + k_3 + s \cdot k_1) = C \cdot p \cdot k_6 \cdot (k_2 + k_3 + s \cdot k_1) + C \cdot s \cdot k_1 k_3 - C \cdot s \cdot k_1 \cdot p \cdot k_6 + $$

$$ ep \cdot (k_4 - s \cdot k_1) \cdot (k_3 - p \cdot k_6) $$

$$ ep \cdot (k_4 + k_5 + p \cdot k_6) \cdot (k_2 + k_3 + s \cdot k_1) = C \cdot p \cdot k_6 k_2 + C \cdot p \cdot k_6 k_3 + C \cdot s \cdot k_1 k_3 + ep \cdot (k_4 - s \cdot k_1) \cdot (k_3 - p \cdot k_6) $$

$$ ep \cdot ((k_4 + k_5 + p \cdot k_6) \cdot (k_2 + k_3 + s \cdot k_1) - (k_4 - s \cdot k_1) \cdot (k_3 - p \cdot k_6)) = C \cdot p \cdot k_6 k_2 + C \cdot p \cdot k_6 k_3 + C \cdot s \cdot k_1 k_3 $$

$$ ep = \frac{C \cdot p \cdot k_6 \cdot (k_2 + k_3) + C \cdot s \cdot k_1 k_3}{(k_4 + k_5 + p \cdot k_6) \cdot (k_2 + k_3 + s \cdot k_1) - (k_4 - s \cdot k_1) \cdot (k_3 - p \cdot k_6)} $$

$$ ep = \frac{C \cdot p \cdot (k_6 k_2 + k_6 k_3) + C \cdot s \cdot k_1 k_3}{(k_4 k_2 + k_5 k_2 + k_5 k_3 + s \cdot k_1 k_3 + s \cdot k_1 k_4 + s \cdot k_1 k_5 + p \cdot k_6 k_2 + p \cdot k_6 k_3 + p \cdot k_6 k_4)} $$

Note that, in the above equation set, we have dropped the dot notation from multiplication of parameters for simplicity of expression. Based on the above argument of symmetry, or by following the same steps for "e_s", we can show that we will have a similar expression with different parameters in the numerator:

$$ e_s = \frac{C \cdot s \cdot (k_1 k_5 + k_1 k_4) + C \cdot p \cdot k_6 k_4}{(k_4 k_2 + k_5 k_2 + k_5 k_3 + s \cdot k_1 k_3 + s \cdot k_1 k_4 + s \cdot k_1 k_5 + p \cdot k_6 k_2 + p \cdot k_6 k_3 + p \cdot k_6 k_4)} $$
With these expressions for $e_s$ and $ep$ at hand, we can now derive an expression for $e$:

$$e = C - e_s - ep$$

$$e = C - \frac{C \cdot s \cdot (k_1 k_5 + k_1 k_4) + C \cdot p \cdot k_6 k_4}{(k_4 k_2 + k_5 k_2 + k_5 k_3 + s \cdot k_1 k_3 + s \cdot k_1 k_4 + s \cdot k_1 k_5 + p \cdot k_6 k_2 + p \cdot k_6 k_3 + p \cdot k_6 k_4)} - \frac{C \cdot p \cdot (k_6 k_2 + k_6 k_3) + C \cdot s \cdot k_1 k_3}{(k_4 k_2 + k_5 k_2 + k_5 k_3 + s \cdot k_1 k_3 + s \cdot k_1 k_4 + s \cdot k_1 k_5 + p \cdot k_6 k_2 + p \cdot k_6 k_3 + p \cdot k_6 k_4)}$$

$$e = C - C \cdot \frac{s \cdot (k_1 k_3 + k_1 k_5 + k_1 k_4) + p \cdot (k_6 k_2 + k_6 k_3 + k_6 k_4)}{(k_4 k_2 + k_5 k_2 + k_5 k_3 + s \cdot k_1 k_3 + s \cdot k_1 k_4 + s \cdot k_1 k_5 + p \cdot k_6 k_2 + p \cdot k_6 k_3 + p \cdot k_6 k_4)}$$

$$e = C \cdot \frac{k_3 k_5 + k_2 k_5 + k_2 k_4}{(k_4 k_2 + k_5 k_2 + k_5 k_3 + s \cdot k_1 k_3 + s \cdot k_1 k_4 + s \cdot k_1 k_5 + p \cdot k_6 k_2 + p \cdot k_6 k_3 + p \cdot k_6 k_4)}$$

We are now ready to substitute all these expressions into the ODE for the product, so to obtain our rate law:

$$\frac{dp}{dt} = \frac{C \cdot p \cdot (k_6 k_2 + k_6 k_3) + C \cdot s \cdot k_1 k_3}{(k_4 k_2 + k_5 k_2 + k_5 k_3 + s \cdot k_1 k_3 + s \cdot k_1 k_4 + s \cdot k_1 k_5 + p \cdot k_6 k_2 + p \cdot k_6 k_3 + p \cdot k_6 k_4)} \cdot k_5 - C \cdot \frac{k_3 k_5 + k_2 k_5 + k_2 k_4}{(k_4 k_2 + k_5 k_2 + k_5 k_3 + s \cdot k_1 k_3 + s \cdot k_1 k_4 + s \cdot k_1 k_5 + p \cdot k_6 k_2 + p \cdot k_6 k_3 + p \cdot k_6 k_4)} \cdot p \cdot k_6$$

$$\frac{dp}{dt} = C \cdot \frac{s \cdot k_1 k_3 k_5 - p \cdot k_2 k_4 k_6}{(k_4 k_2 + k_5 k_2 + k_5 k_3 + s \cdot k_1 k_3 + s \cdot k_1 k_4 + s \cdot k_1 k_5 + p \cdot k_6 k_2 + p \cdot k_6 k_3 + p \cdot k_6 k_4)}$$

We can somewhat simplify this expression by defining the following composite rate constants:

$$K_S = \frac{k_3 k_5 + k_2 k_5 + k_2 k_4}{k_1 \cdot (k_3 + k_4 + k_5)}$$

$$K_P = \frac{k_3 k_5 + k_2 k_5 + k_2 k_4}{k_6 \cdot (k_2 + k_3 + k_4)}$$

$$k_{\text{cat}}^{+} = \frac{k_3 k_5}{k_3 + k_4 + k_5}$$

$$k_{\text{cat}}^{-} = \frac{k_2 k_4}{k_2 + k_3 + k_4}$$

and substituting them into the rate expression from above, to get:

$$\frac{dp}{dt} = v = C \cdot \frac{k_{\text{cat}}^{+}}{K_S} \cdot \frac{s - p \cdot \frac{k_{\text{cat}}^{-}}{K_P}}{1 + \frac{k_{\text{cat}}^{+}}{K_S} + \frac{s}{K_S}}$$

This reaction rate is referred to as the Haldane kinetic rate law, named after Jack Burden Sanderson Haldane (5 November 1892 1 December 1964). It can be re-expressed by recognizing the fact that the fraction entering as a multiplier for the product concentration is actually equivalent to the equilibrium constant of the reaction scheme drawn above, at the beginning of this section, when we assume the reaction proceeding in the forward direction, i.e.
towards product formation:

$$ \frac{k_{\text{cat}}^{-}}{K_{\text{P}}} = \frac{k_2 k_4 k_6}{k_1 k_3 k_5} = 1/K_{\text{eq}} $$

This allows us to re-express the Haldane rate law as:

$$ v = C \cdot k_{\text{cat}}^{+} \cdot \frac{s/K_{\text{S}}}{1 + \frac{p}{K_{\text{P}}} + \frac{s}{K_{\text{S}}}} \cdot \left(1 - \frac{p/s}{K_{\text{eq}}}\right) $$

This re-arranged expression is interesting because we can recognize that the last term is related to the thermodynamic Gibbs free energy of the reaction, allowing us to finally derive:

$$ v = C \cdot k_{\text{cat}}^{+} \cdot \frac{s/K_{\text{S}}}{1 + p/K_{\text{P}} + s/K_{\text{S}}} \cdot \left(1 - e^{\Delta_{\text{r}} G' / RT}\right) $$

where $\Delta_{\text{r}}G'$ is the Gibbs free energy of reaction for given substrate and product levels, considering forward direction, and R and T stand for the gas constant and temperature respectively. This rate law shows that forward reaction rate will be independent of thermodynamics, when the reaction free energy is highly negative (i.e. thermodynamically highly favored), but the reaction rate will decrease as Gibbs free energy gets close to zero.

A second, faster derivation of this rate law is found by noting that the ODEs for $\frac{de_s}{dt}$ and $\frac{de_p}{dt}$ are linear in $e$, $e_s$ and $e_p$, and can therefore be solved with linear matrix algebra. One may write:

$$ \begin{pmatrix} s & k_1 & -(k_2 + k_3) & k_4 \\ p & k_6 & k_3 & -(k_4 + k_5) \\ 1 & 1 & 1 & 1 \end{pmatrix} \begin{pmatrix} e \\ e_s \\ e_p \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ C \end{pmatrix}, \qquad (3.20) $$

where the first two rows of the matrix correspond to $\frac{de_s}{dt} = 0$ and $\frac{de_p}{dt} = 0$, and the last row represents conservation of total enzyme concentration. The equilibrium concentrations of $e$, $e_s$ and $e_p$ are then found by left-multiplying both sides of the equation by the inverse of this matrix. The obtained results are the same as given above.

### 3.6.4. Derivation of two-substrate, reversible rate law with simultaneous binding

The two-substrate case is described by the following reaction scheme:

$$ S_1 + S_2 + E \underset{k_2}{\stackrel{k_1}{\rightleftharpoons}} ES_1S_2 \underset{k_4}{\stackrel{k_3}{\rightleftharpoons}} EP_1P_2 \underset{k_6}{\stackrel{k_5}{\rightleftharpoons}} P_1 + P_2 + E, $$

Where we assume that binding and unbinding of the substrates and products occurs simultaneously. Proceeding as above we let $e$, $e_{s1}s_2$, $ep_1p_2$, $s_1$, $s_2$, $p_1$ and $p_2$ denote the concentrations of E, $ES_1S_2$, $EP_1P_2$, $S_1$, $S_2$, $P_1$ and $P_2$ respectively. The differential equations for $e_{s1}s_2$, $ep_1p_2$ and $p_1 + p_2$ are:

$$ \frac{des_1s_2}{dt} = e \cdot s_1 \cdot s_2 \cdot k_1 + ep_1p_2 \cdot k_4 - es_1s_2 \cdot (k_2 + k_3) $$

$$ \frac{dep_1p_2}{dt} = e \cdot p_1 \cdot p_2 \cdot k_6 + es_1s_2 \cdot k_3 - ep_1p_2 \cdot (k_4 + k_5) $$

$$ \frac{d(p_1 + p_2)}{dt} = ep_1p_2 \cdot k_5 - e \cdot p_1 \cdot p_2 \cdot k_6. $$
Proceeding as in the single substrate case, we note that the ODEs for $\frac{des_1s_2}{dt}$ and $\frac{dep_1p_2}{dt}$ are linear in $e$, $es_1s_2$ and $ep_1p_2$, and that the total enzyme concentration $e + es_1s_2 + ep_1p_2$ is constant, denoted $C$.

$$ \begin{pmatrix} s_1s_2k_1 & -(k_2+k_3) & k_4 \\ p_1p_2k_6 & k_3 & -(k_4+k_5) \\ 1 & 1 & 1 \end{pmatrix} \begin{pmatrix} e \\ es_1s_2 \\ ep_1p_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ C \end{pmatrix}. \quad (3.21) $$

We therefore see that the results for the two-substrate case are the same as for the single substrate case, with $s$ replaced by $s_1s_2$ and $p$ replaced by $p_1p_2$. This result is dependent on the assumption that binding/unbinding of substrates/products occurs simultaneously.

## 3.7. A simple model illustrating product activation

This model demonstrates that allosteric regulation of an enzymatic reaction by its product can create a bistable system. In this simple example, we consider enzymatic production of a metabolite (labelled 'x') and its non-enzymatic consumption. It is assumed that the metabolite allosterically regulates the enzyme that produces it. The listing uses the Antimony format [48] which can be easily converted into SBML [49]. An online converter can be found at sys-bio.github.io/makesbml/

1. // The following model admits three steady-states at:

2. // x = 0.325, x = 1.671, and x = 0.873

3. // The first reaction step '-> x' uses a rate law that models

4. // positive feedback via the product x. The constant 0.2

5. // is to ensure that the lower steady-state is non-zero.

6. // The statement 'ext Xo' indicates that the species Xo is fixed.

7. ext Xo

8. Xo -> x; (vo*x^n)/(1 + x^n) + 0.2

10. x ->; k1*x

11.

12. k1 = 0.65

13. n = 4; vo = 1

14. x = 0

<caption>Listing 3.1: Model illustrating bistability</caption>

1. # Equivalent model as a differential equation in python:

2. def ode(x, t):

3. vo = 1

4. n = 4

5. k1 = 0.65

6. return [(vo*x**n)/(1 + x**n) + 0.2) - k1*x]

<caption>Listing 3.2: Equivalent model as a differential equation in python</caption>

# Solutions to problems

## Problem 3.1 (An irreversible reaction with simultaneous binding)
$$ (a) \qquad E + S_1 + S_2 \underset{k_2}{\overset{k_1}{\rightleftharpoons}} ES_1 S_2 \xrightarrow{k_3} E + P_1 + P_2 \qquad (3.22) $$

$$ (b) \qquad \frac{dp}{dt} = k_3 \frac{s_1 s_2 C}{s_1 s_2 + \frac{k_2 + k_3}{k_1}}, \qquad (3.23) $$

where $p = [P_1 + P_2]$ and $C = [E] + [ES_1S_2]$.

### Problem 3.2 (A reversible reaction with simultaneous binding)

$$ (a) \qquad E + S_1 + S_2 \underset{k_2}{\overset{k_1}{\rightleftharpoons}} ES_1 S_2 \underset{k_4}{\overset{k_3}{\rightleftharpoons}} E + P_1 + P_2 \qquad (3.24) $$

$$ (b) \qquad \frac{dp}{dt} = k_3 \frac{C(s_1 s_2 - \frac{k_2 k_4}{k_1 k_3} p)}{s_1 s_2 + \frac{k_4}{k_1} p + \frac{k_2 + k_3}{k_1}} \qquad (3.25) $$

where $p = [P_1 + P_2]$ and $C = [E] + [ES_1S_2]$.

### Problem 3.3 (An irreversible reaction with sequential binding)

$$ (a) \qquad E + S_1 \underset{k_2}{\overset{k_1}{\rightleftharpoons}} ES_1 \\ ES_1 + S_2 \underset{k_4}{\overset{k_3}{\rightleftharpoons}} ES_1 S_2 \xrightarrow{k_5} E + P_1 + P_2 \qquad (3.26) $$

$$ (b) \qquad \frac{dp}{dt} = k_5 \frac{s_1 s_2 C}{s_1 s_2 + s_1 \frac{k_4+k_5}{k_3} + s_2 \frac{k_5}{k_3} + \frac{k_2}{k_1 k_3} (k_1 + k_5)}, \qquad (3.27) $$

where $p = [P_1 + P_2]$ and $C = [E] + [ES_1] + [ES_2]$.

### Problem 3.4 (An irreversible reaction with random-order binding)

$$ (a) \qquad E + S_1 \underset{k_2}{\overset{k_1}{\rightleftharpoons}} ES_1 \\ ES_1 + S_2 \underset{k_4}{\overset{k_3}{\rightleftharpoons}} ES_1 S_2 \\ E + S_2 \underset{k_6}{\overset{k_5}{\rightleftharpoons}} ES_2 \\ ES_2 + S_1 \underset{k_8}{\overset{k_7}{\rightleftharpoons}} ES_1 S_2 \\ ES_1 S_2 \xrightarrow{k_9} E + P_1 + P_2 \qquad (3.28) $$

$$ (b) \qquad \frac{dp}{dt} = k_9 \frac{C s_1 s_2 (k_1 k_3 (k_6 + k_7 s_1) + k_5 k_7 (k_2 + k_3 s_2))}{s_1 A(s_1) + s_2 B(s_2) + s_1 s_2 C(s_1, s_2) + D} \qquad (3.29) $$
where $p = [P_1 + P_2]$, $C = [E] + [ES_1] + [ES_2] + [ES_1S_2]$, and

$$A(s_1) = k_1k_6(k_4 + k_8 + k_9) + k_7(k_0 + k_4)(k_2 + k_1s_1)$$

$$B(s_2) = k_2k_5(k_4 + k_8 + k_9) + k_3(k_0 + k_8)(k_6 + k_5s_2)$$

$$C(s_1, s_2) = k_1k_3(k_6 + k_8 + k_7s_1) + k_5k_7(k_2 + k_4 + k_3s_2) + k_3k_7k_9$$

$$D = k_2k_6(k_4 + k_8 + k_9)$$
Bibliography

[1] G. Gottschalk. *Bacterial Metabolism*. Springer Series in Microbiology. Springer, New York, 1985. ISBN 978-0-387-96153-8. doi: https://doi.org/10.1007/978-1-4612-1072-6.

[2] F. C. Neidhardt, J. L. Ingraham, and M. Schaechter. *Physiology go the bacterial cell: A molecular approach*. Sinauer Associates, 1990. ISBN 0878936084.

[3] H. V. Westerhoff, K. J. Hellingwerf, and K. Van Dam. Thermodynamic efficiency of microbial growth is low but optimal for maximal growth rate. *Proc Natl Acad Sci U S A*, 80(1):305-9, 1983. ISSN 0027-8424 (Print) 1091-6490 (Electronic) 0027-8424 (Linking). doi: 10.1073/pnas.80.1.305.

[4] E. Branscomb and M. J. Russell. Turnstiles and bifurcators: the disequilibrium converting engines that put metabolism on the road. *Biochim Biophys Acta*, 1827(2):62-78, 2013. ISSN 0006-3002 (Print) 0006-3002 (Linking). doi: 10.1016/j.bbabio.2012.10.003.

[5] C. Zerfass, M. Asally, and O. S. Soyer. Interrogating metabolism as an electron flow system. *Curr Opin Syst Biol*, 13:59-67, 2019. ISSN 2452-3100 (Print) 2452-3100 (Linking). doi: 10.1016/j.coisb.2018.10.001.

[6] W. H. Schlesinger and E. S. Bernhardt. *Biogeochemistry: an analysis of global change*. Academic Press, 2013. ISBN 0123858747.

[7] A. Jinich, A. Flamholz, H. Ren, S. J. Kim, B. Sanchez-Lengeling, C. A. R. Cotton, E. Noor, A. Aspuru-Guzik, and A. Bar-Even. Quantum chemistry reveals thermodynamic principles of redox biochemistry. *PLoS Comput Biol*, 14 (10):e1006471, 2018. ISSN 1553-7358 (Electronic) 1553-734X (Linking). doi: 10.1371/journal.pcbi.1006471.

[8] U. Barenholz, D. Davidi, E. Reznik, Y. Bar-On, N. Antonovsky, E. Noor, and R. Milo. Design principles of autocatalytic cycles constrain enzyme kinetics and force low substrate saturation at flux branch points. *Elife*, 6, 2017. ISSN 2050-084X (Electronic) 2050-084X (Linking). doi: 10.7554/eLife.20667.

[9] D. C. LaPorte, K. Walsh, and Jr. Koshland, D. E. The branch point effect. ultrasensitivity and subsensitivity to metabolic control. *J Biol Chem*, 259(22):14068-75, 1984. ISSN 0021-9258 (Print) 0021-9258 (Linking).

[10] M. E. Beber, C. Fretter, S. Jain, N. Sonnenschein, M. Muller-Hannemann, and M. T. Hutt. Artefacts in statistical analyses of network motifs: general framework and application to metabolic networks. *J R Soc Interface*, 9(77): 3426-35, 2012. ISSN 1742-5662 (Electronic) 1742-5662 (Linking). doi: 10.1098/rsif.2012.0490.

[11] J.G. Reich and E.E. Selkov. *Energy metabolism of the cell : a theoretical treatise*. Academic Press, London ; New York, 1981. ISBN 0125859201.

[12] J. H. Hofmeyr, H. Kacser, and K. J. van der Merwe. Metabolic control analysis of moiety-conserved cycles. *Eur J Biochem*, 155(3):631-41, 1986. ISSN 0014-2956 (Print) 0014-2956 (Linking). doi: 10.1111/j.1432-1033.1986.tb09534.x.

[13] H. M. Sauro. Moiety-conserved cycles and metabolic control analysis: problems in sequestration and metabolic channelling. *Biosystems*, 33(1):55-67, 1994. ISSN 0303-2647 (Print) 0303-2647 (Linking). doi: 10.1016/0303-2647(94)90061-2.
[14] T. S. Hatakeyama and C. Furusawa. Metabolic dynamics restricted by conserved carriers: Jamming and feedback. *PLoS Comput Biol*, 13(11):e1005847, 2017. ISSN 1553-7358 (Electronic) 1553-734X (Linking). doi: 10.1371/journal.pcbi.1005847.

[15] Robert West, Hadrien Delattre, Elad Noor, Elisenda Feliu, and Orkun Soyer. Dynamics of co-substrate pools can constrain and regulate metabolic fluxes. *Elife*, 12:e84379, 2023.

[16] N. C. Price, R. A. Dwek, R. G. Ratcliffe, and M. R. Wormald. *Principles and Problems in Physical Chemistry for Biochemists*. Oxford University Press, 1997. ISBN 0198792816.

[17] H. C. V. Ness. *Understanding thermodynamics*. Dover Publications, 1983. ISBN 0486632776.

[18] E. Noor, A. Bar-Even, A. Flamholz, Y. Lubling, D. Davidi, and R. Milo. An integrated open framework for thermodynamics of reactions that combines accuracy and coverage. *Bioinformatics*, 28(15):2037–44, 2012. ISSN 1367-4811 (Electronic) 1367-4803 (Print) 1367-4803 (Linking). doi: 10.1093/bioinformatics/bts317.

[19] M. E. Beber, M. G. Gollub, D. Mozaffari, K. M. Shebek, A. I. Flamholz, R. Milo, and E. Noor. eQuilibrator 3.0: a database solution for thermodynamic constant estimation. *Nucleic Acids Res*, 50(D1):D603–D609, 2022. ISSN 1362-4962 (Electronic) 0305-1048 (Print) 0305-1048 (Linking). doi: 10.1093/nar/gkab1106.

[20] K. J. Mysels. Textbook errors: VII. the laws of reaction rates and of equilibrium. *Journal of Chemical Education*, 33(4), 1956.

[21] E. A. Guggenheim. Textbook errors: IX. more about the laws of reaction rates and of equilibrium. *Journal of Chemical Education*, 33(11), 1956.

[22] A. Fersht. *Structure and Mechanism in Protein Science: Guide to Enzyme Catalysis and Protein Folding*. W. H. Freeman & Co., 3rd edition, 1998. ISBN 0716732688.

[23] C. I. Branden and C. Tooze. *Introduction to Protein Structure*. Garland Science, 2nd edition, 1999.

[24] D. A. Beard. Simulation of cellular biochemical system kinetics. *Wiley Interdiscip Rev Syst Biol Med*, 3(2): 136–46, 2011. ISSN 1939-005X (Electronic) 1939-5094 (Print) 1939-005X (Linking). doi: 10.1002/wsbm.116.

[25] A. Cornish-Bowden. *Fundamentals of Enzyme Kinetics*. Wiley-Blackwell, 4th edition, 2012. ISBN 9783527330744.

[26] J. B. S. Haldane. *Enzymes*. MIT Press, 1965. ISBN 0262580039.

[27] C. Y. Hoh and R. Cord-Ruwisch. A practical kinetic model that considers endproduct inhibition in anaerobic digestion processes by including the equiibrium constant. *Biotechnol Bioeng*, 51, 1996.

[28] T. Grosskopf and O. S. Soyer. Microbial diversity arising from thermodynamic constraints. *ISME J*, 10(11): 2725–2733, 2016. ISSN 1751-7370 (Electronic) 1751-7362 (Linking). doi: 10.1038/ismej.2016.49.

[29] André B. Canelas, Cor Ras, Angela ten Pierick, Walter M. van Gulik, and Joseph J. Heijnen. An in vivo data-driven framework for classification and quantification of enzyme kinetics and determination of apparent thermodynamic data. *Metabolic Engineering*, 13(3):294–306, 2011. doi: 10.1016/j.ymben.2011.02.005.

[30] J. O. Park, S. A. Rubin, Y. F. Xu, D. Amador-Noguez, J. Fan, T. Shlomi, and J. D. Rabinowitz. Metabolite concentrations, fluxes and free energies imply efficient enzyme usage. *Nat Chem Biol*, 12(7):482–9, 2016. ISSN 1552-4469 (Electronic) 1552-4450 (Linking). doi: 10.1038/nchembio.2077.

[31] C. Hayes, E. Feliu, and O. S. Soyer. Multisite enzymes as a mechanism for bistability in reaction networks. *ACS Synth Biol*, 11(2):596–607, 2022. ISSN 2161-5063 (Electronic) 2161-5063 (Linking). doi: 10.1021/acssynbio.1c00272.
[32] D. A. Beard and H. Qian. Relationship between thermodynamic driving force and one-way fluxes in reversible processes. *PLoS One*, 2(1):e144, 2007. ISSN 1932-6203 (Electronic) 1932-6203 (Linking). doi: 10.1371/journal.pone.0000144.

[33] E. Noor, A. Flamholz, W. Liebermeister, A. Bar-Even, and R. Milo. A note on the kinetics of enzyme action: a decomposition that highlights thermodynamic effects. *FEBS Lett*, 587(17):2772-7, 2013. ISSN 1873-3468 (Electronic) 0014-5793 (Linking). doi: 10.1016/j.febslet.2013.07.028.

[34] A. Cornish-Bowden and M. L. Cardenas. Information transfer in metabolic pathways. effects of irreversible steps in computer models. *Eur J Biochem*, 268(24):6616-24, 2001. ISSN 0014-2956 (Print) 0014-2956 (Linking). doi: 10.1046/j.0014-2956.2001.02616.x.

[35] L. Shen, M. Kohlhaas, J. Enoki, R. Meier, B. Schonenberger, R. Wohlgemuth, R. Kourist, F. Niemeyer, D. van Niekerk, C. Brasen, J. Niemeyer, J. Snoep, and B. Siebers. A combined experimental and modelling approach for the Weimberg pathway optimisation. *Nat Commun*, 11(1):1098, 2020. ISSN 2041-1723 (Electronic) 2041-1723 (Linking). doi: 10.1038/s41467-020-14830-y.

[36] P. T. Silvanus. *Calculus made easy*. Macmillan, 1914.

[37] S. H. Strogatz. *Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering*. CRC Press, 2000. ISBN 0738204536.

[38] H. M. Sauro. *Enzyme Kinetics for Systems Biology*. Ambrosius Publishing, 2012. ISBN 0982477333.

[39] H. Kacser and J. A. Burns. The control of flux. *Biochem Soc Trans*, 23(2):341-66, 1995. ISSN 0300-5127 (Print) 0300-5127 (Linking). doi: 10.1042/bst0230341.

[40] I. Prigogine, R. Lefever, A. Goldbeter, and M. Herschkowitz-Kaufman. Symmetry breaking instabilities in biological systems. *Nature*, 223(5209):913-6, 1969. ISSN 0028-0836 (Print) 0028-0836 (Linking). doi: 10.1038/223913a0.

[41] A. Goldbeter. Dissipative structures in biological systems: bistability, oscillations, spatial patterns and waves. *Philos Trans A Math Phys Eng Sci*, 376(2124), 2018. ISSN 1471-2962 (Electronic) 1364-503X (Linking). doi: 10.1098/rsta.2017.0376.

[42] E. Simsek and M. Kim. The emergence of metabolic heterogeneity and diverse growth responses in isogenic bacterial cells. *ISME J*, 12(5):1199-1209, 2018. ISSN 1751-7370 (Electronic) 1751-7362 (Linking). doi: 10.1038/s41396-017-0036-2.

[43] A. Z. Rosenthal, Y. Qi, S. Hormoz, J. Park, S. H. Li, and M. B. Elowitz. Metabolic interactions between dynamic bacterial subpopulations. *Elife*, 7, 2018. ISSN 2050-084X (Electronic) 2050-084X (Linking). doi: 10.7554/eLife.33099.

[44] Dan Davidi, Elad Noor, Wolfram Liebermeister, Arren Bar-Even, Avi Flamholz, Katja Tummler, Uri Barenholz, Miki Goldenfeld, Tomer Shlomi, and Ron Milo. Global characterization of in vivo enzyme catalytic rates and their correspondence to in vitro $k_{cat}$ measurements. *Proceedings of the National Academy of Sciences*, 113 (12):3401-3406, 2016. doi: 10.1073/pnas.1514240113.

[45] L. Gerosa, B. R. Haverkorn van Rijsewijk, D. Christodoulou, K. Kochanowski, T. S. Schmidt, E. Noor, and U. Sauer. Pseudo-transition analysis identifies the key regulators of dynamic metabolic adaptations from steady-state data. *Cell Syst*, 1(4):270-82, 2015. ISSN 2405-4712 (Print) 2405-4712 (Linking). doi: 10.1016/j.cels.2015.09.008.
[46] D. A. Beard and H. Qian. *Chemical Biophysics: Quantitative Analysis of Cellular Systems*. Cambridge University Press, 2008. ISBN 0521870704.

[47] H. M. Sauro. *Systems Biology: An Introduction to Metabolic Control Analysis*. Ambrosius Publishing, 2018. ISBN 0982477368.

[48] Lucian P Smith, Frank T Bergmann, Deepak Chandran, and Herbert M Sauro. *Antimony: a modular model definition language*. *Bioinformatics*, 25(18):2452–2454, 2009.

[49] Michael Hucka, Andrew Finney, Herbert M Sauro, Hamid Bolouri, John C Doyle, Hiroaki Kitano, Adam P Arkin, Benjamin J Bornstein, Dennis Bray, Athel Cornish-Bowden, et al. The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. *Bioinformatics*, 19(4): 524–531, 2003.