Metabolic pathways and complexity

<figure><img src="image_18.png" /><figcaption>Figure 3.1: (A & B) Upper glycolysis pathway as a linear pathway, with or without co-substrates. Note that the co-substrates 'connect' this pathway to a large number of other reactions that also use these same co-substrates. (C) Upper glycolysis pathway, together with the pentose-phosphate pathway. Notice metabolites participating in both.</figcaption></figure>

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' [63]. Known mostly through the names of their discoverers, these include the Entner-Doudoroff (ED), Embden-Meyerhof-Parnas (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 [64]. 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 (Fig. 3.1).

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.
<figure><img src="image_23.png" /><figcaption>Figure 2.5: Number of molecules per cell in a population (for example protein or mRNA) – The red line is the population mean, which is often the value we use (for modeling). However, the values at a single-cell level can differ several fold.</figcaption></figure>

Finally, we need to zoom out from a single-cell (or average) view of a cell and consider the heterogeneity at the population level. This heterogeneity is often neglected, and we use a single number to describe a concentration of a molecule in a cell/compartment – an average value of the population. However, biological processes are stochastic (noisy), and the actual molecule numbers follow a certain distribution (Figure 2.5), which can be characterized by mean and variance. The effect of the heterogeneity becomes especially important at low copy numbers.

The heterogeneity in molecule copy numbers leads to a heterogeneity in cell phenotypes such as generation time, cell size, stress tolerance, and others. The heterogeneity of the population can affect fitness in a positive or negative way, depending on the conditions. For example, when a cell population encounters an unexpected environment, a certain subpopulation might be better suited to survive. In a different environment, another subpopulation might thrive. We can view this as a microbial “bet-hedging” which increases the chances that at least some part of a population will survive the new conditions. However, when cells try to maximize the growth rate, the variability in the population can decrease fitness because it decreases the average population growth rate [23]. This topic is discussed in detail in Chapter 13.

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

## Experimental methods 2.C Quantification of biomass composition

We can measure biomass composition at different levels of detail – from a coarse elemental or macromolecular composition of an average cell to the quantities of individual molecules in each cellular compartment.

To quantify the main chemical elements (CHNOPS), we can use devices called elemental analyzers. The main macromolecular components – the total protein, lipid, carbohydrate, DNA, and RNA content – can be quantified with simple assays such as detection with fluorescent dyes, chemical reactions that lead to color change, or extraction and weighing of a component. Going into more detail typically requires more sophisticated methods such as liquid or gas chromatography (LC, GC), mass spectrometry (MS) or nuclear magnetic resonance (NMR). For example, for proteins, we can measure an average amino acid composition, and for lipids, the main lipid classes (glycerophospholipids, sphingolipids, sterols, etc.).

If we go down to the level of individual molecules, we enter fields of study collectively termed as *omics*, which aim to characterize and quantify certain pools of biomolecules. Omics methods typically involve high-throughput measurements of hundreds or thousands of different molecules and require a lot of resources (specialized equipment, computational resources) and expertise. The classic omics fields include genomics, transcriptomics, and proteomics which study DNA, RNA and proteins, respectively. Other examples include metabolomics which focuses on small metabolites or fluxomics which measures metabolic fluxes (for example ¹³C metabolic flux analysis).

Combinations of different omics can help us obtain other parameters that are difficult to measure. For example, turnover numbers of enzymes ($k_{cat}$) are notoriously difficult to quantify because the measurements are error-prone and low-throughput. With proteomics and fluxomics data we can calculate apparent turnover numbers ($k_{app}$) at various conditions (see Figure 2.9) and use the maximum value ($k_{app}^{max}$) as an estimate of *in vivo* $k_{cat}$ [24].
Economic analogy 3.A

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 (Fig. 3.1). 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.

<figure><img src="image_19.png" /><figcaption>Figure 3.2: A simplified map of central metabolism, particularly highlighting interconnections among different processes (i.e. pathways) through the NAD⁺ / NADH co-substrate pair.</figcaption></figure>

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, Philosophical Remarks Box 3.C).

3.1.2. Metabolism - coarse grained views

The highly connected nature of metabolism makes it difficult to understand its overall dynamics just from individual
pathways. It also makes it hard to conceptualise 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.
<figure><img src="image_25.png" /></figure>

Experimental methods 2.D Examples of biomass quantification methods

<table><thead><tr><th>Component/parameter</th><th>Examples of quantification methods</th></tr></thead><tbody><tr><td>Cell size</td><td>microscopically</td></tr><tr><td>Dry cell mass</td><td>weighing of a defined amount of dry cells</td></tr><tr><td>Buoyant density</td><td>Percoll gradient</td></tr><tr><td>Protein</td><td>colorimetric (Bradford assay; Lowry assay)</td></tr><tr><td>Lipid</td><td>weighing of extracted and dried lipids</td></tr><tr><td>Carbohydrates</td><td>colorimetric (anthrone assay; phenol-sulphuric acid assay)</td></tr><tr><td>RNA</td><td>fluorimetric (RiboGreen), spectrophotometric</td></tr><tr><td>DNA</td><td>fluorimetric (PicoGreen, Hoechst), spectrophotometric</td></tr><tr><td>Amino acids/lipid classes</td><td>LC/MS, GC/MS</td></tr><tr><td>Genomics</td><td>next-generation sequencing (NGS) - Illumina, PacBio, Nanopore</td></tr><tr><td>Transcriptomics</td><td>NGS (RNA-seq), DNA microarrays</td></tr><tr><td>Proteomics/metabolomics</td><td>LC/MS, GC/MS, NMR</td></tr></tbody></table>

To visualize composition data, consider using Voronoi diagrams instead of the traditional pie charts or bar plots. An online tool is available at [bionic-vis.biologie.uni-greifswald.de](http://bionic-vis.biologie.uni-greifswald.de) for proteomics data, but there is also a tool that works with any type of input data (GitLab repository on the [book website](https://www.wiley.com/go/bionic-vis)).

## 2.5 Macromolecule synthesis and the resources needed

In the previous sections, we have explored the diversity of nature and the abundance of biological molecules. The combination of smaller building blocks into functional units, let it be proteins, membranes, or DNA that conserves information about the organism, is the major stepping stone from an unorganized pack of molecules into what we could call a living system. Therefore, in this section, we will consider the coordination of how cell components are produced, focusing on the biosynthesis of macromolecules.

The overall cell growth can be called *self-replication*: a cell makes a copy of itself by synthesizing macromolecules by using molecules it either produces or takes up from the environment, all in the right amounts and proportions. How, how fast and how big two cells rise from a single parent cell is the question we explore in our field, and often try to complement the experimental observations (*what* and *in what amounts* produced) with computational models of growth (*how* and *why like that*).

In general, three essential types of resources are needed for the synthesizing of macromolecules: (1) small molecule precursors, (2) catalysts, and (3) physical space/volume for the process to happen (Figure 2.6). We will thus discuss how these resources are primed and used for macromolecule synthesis, together with different considerations surrounding each type of these resources. We will first start with discussing the “demand” side of the balance, requirement for the small molecule precursors of the cells, and will continue to zoom out towards the whole-cell economy of volume.

### 2.5.1 Precursors of macromolecules

Biosynthesis of macromolecule precursors (e.g. amino acids, nucleotides, energy equivalents) is an important part of every metabolic network. Many microorganisms can grow on a very limited number of nutrients (in the laboratory context, the so-called minimal media), which usually consist of a single source of carbon, nitrogen, phosphorus, and sulfur. For instance, a minimal growth medium with glucose as the sole carbon source can fully support growth: glucose enters glycolysis as the main energy-harvesting route; however, some of the glycolytic intermediates serve as substrates for, e.g. amino acid, lipid, or nucleotide biosynthesis.

A particularly interesting fact is that metabolic networks can be described as *bow-tie* structures [25]: a large variety of nutrients can be converted into a very small number (usually counted up to 12) essential metabolic intermediates, which give rise to, again, a diverse set of molecules (for a detailed discussion, see Chapter 3). This provides two important insights into metabolic networks. First, this plasticity of the metabolic networks allows organisms to grow in various environments, where different nutrients are avail-
**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 metabolism respectively, and their coupling presents the whole cell metabolism (Fig. 3.3 A). This coarse-grained model is widely used (e.g. [64, 65]). However, it is not always clear how to partition various pathways and reactions as anabolic and catabolic, and the notion of metabolism organised 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) [66, 67]. (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 visualise 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.

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 [63]. 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 [68].
<figure><img src="image_26.png" /></figure>

Economics analogy 2.E Comparing the cell to a bakery

The diversity of metabolic intermediates/end products, stemming from small number of nutrients (e.g. minimal mineral media for yeast growth, containing glucose, ammonium, phosphate and sulphate salts), can be imagined as a bakery. Every pastry starts with a small array of ingredients (flour, water, salt, sugar, ...) and using some machinery (e.g. ovens), one ends up baking bread, pretzels, cookies, muffins etc., which are way diverse in their features, compared to the starting mixture. Likewise, by taking only a handful of compounds, cells, especially microorganisms, can synthesize most of the molecules they need to eventually replicate.

able. Second, because of this organization, the biosynthesis of macromolecule precursors competes for the same starting molecules independently of the initial nutrients.

2.5.2 Catalysts needed for macromolecule synthesis

Many steps of the biosynthesis of macromolecules, as discussed previously, need catalysis to proceed. Therefore, another kind of investment into macromolecule synthesis is expression of necessary proteins and RNAs (in the latter case - ribosomal RNA). Expression of proteins, starting from transcription of messenger RNAs, their translation into proteins, folding, and degradation, involve many steps with energy investment (ATP hydrolysis) and consume large amounts of precursors (nucleotides, amino acids). Talking in energetic terms alone, protein expression accounts for about 40% of energy investments in yeast *S. cerevisiae* [26], and the energy investments for every stage of protein expression are illustrated in Table 2.2 for typical bacterial and eukaryotic cells. This concerted action of several systems, as described above, with substantial investments at every intermediate step, means that these investments thus happen on two levels: investments in the metabolic machinery and in the machinery, producing proteins themselves. We will consider these two levels in the following.

**Metabolic enzymes.** First, metabolic enzymes need to be expressed to convert nutrients into biosynthe-

<figure><img src="image_27.png" /><figcaption>Figure 2.6: The basic building blocks of cells are small molecule precursors – The precursors are needed to make catalysts such as enzymes and machines. In turn, these catalysts synthesize both the precursors and themselves, forming a self-replicating system. These processes need to be coordinated while being constrained by space.</figcaption></figure>

<table><caption>Table 2.2: The estimated energetic costs (units of ATP hydrolysis) of biosynthesis of a gene, as computed by [27]. The estimates are represented as functions of the following parameters: $L_g$, gene length; $N_r$, the steady-state number of mRNAs; $L_{r,pre}$ and $L_{r,mat}$, the length of precursor and mature mRNA, respectively; $\delta_r$, the degradation rate of mRNA; $t$, division time of a cell; $N_p$, the steady-state number of protein molecules; $L_p$, length of the protein chain; $\bar{c}_{AA}$, average cost of an amino acid; $\delta_p$, the degradation rate of proteins.</caption><thead><tr><th>Expression stage</th><th>Bacteria</th><th>Eukaryotes</th></tr></thead><tbody><tr><td>DNA synthesis</td><td>101 L<sub>g</sub></td><td>263 L<sub>g</sub> (×2 for diploids)</td></tr><tr><td>RNA transcription</td><td>2 N<sub>r</sub> L<sub>g</sub>(23 + δ<sub>r</sub>t)</td><td>N<sub>r</sub>(46 × L<sub>r,mat</sub> + 2.17 × δ<sub>r</sub>tL<sub>r,pre</sub>)</td></tr><tr><td>Protein synthesis</td><td></td><td>N<sub>p</sub>L<sub>p</sub>[(c&#773;<sub>AA</sub> - 1) + 5δ<sub>p</sub>t]</td></tr></tbody></table>
Coarse-grained models of cell metabolism

<figure><img src="image_20.png" /><figcaption>Figure 3.3: (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>

Metabolism on a redox ladder!

<figure><img src="image_21.png" /><figcaption>Figure 3.4: 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_28.png" /><figcaption>Figure 2.7: Distributions of the $k_{cat}$ and $K_M$ values (in $s^{-1}$ and $mM$, respectively), collected for *E. coli*, yeast and human enzymes – The vertical solid line depicts the median of each distribution. Values were collected from the BRENDA database, release 2022.1 [29].</figcaption></figure>

sis precursors. Some enzymes are active only in a form of complexes, which also creates a demand to express proteins at defined ratios. Enzymes and their complexes come in different sizes and flavors, and their activity can be described (in very coarse-grained way, for more details see Chapter 3) by two kinetic aspects: the efficacy (represented by the turnover number $k_{cat}$) and substrate specificity (Michaelis constant $K_M$) of an enzyme. Importantly, these two parameters are intertwined: high substrate specificity usually comes at the cost of efficacy and vice versa. Therefore, although some enzymes tend towards extremes in terms of their specificity or efficacy, most of the enzymes land close to the average/median values of these parameters, when considering the distribution of enzyme parameters among different organisms [28] (Figure 2.7).

The metabolic networks need to work in a concerted manner, even though different enzymes need to perform different amounts of “work” (described as metabolite flux through these enzymes, $v$). Thus, even given the similarities in “average” (or “moderate”) enzyme properties, the expression of proteins and the abundance of their substrates span several orders of magnitude. Based on the kinetic interpretation of enzyme kinetic parameters, we can link them to either expression level of the enzyme $e$ ($e \propto \frac{v}{k_{cat}}$) or substrate concentration $s$ (usually, $0.1K_M \le s \le 10K_M$). Note that for substrate concentrations, the suggested range (order-of-magnitude difference from the $K_M$ to each side) is arbitrary, yet supported by empirical observations. On the higher end, the benefit from high substrate concentration becomes negligible (saturation kinetics) as the concentration moves from the order of magnitude of $K_M$ (see Exercises for an example). The lower bound of concentrations is defined through high demand of enzymes: in order to sustain flux, a lot of enzyme would have to be produced. As cells have a finite volume to accommodate proteins, such a strategy works only for a very small number of enzymes. Taken together the limitations on the both sides of the spectrum, enzyme kinetics set the bounds for the concentrations of metabolites in the cells.

To illustrate the diversity of enzyme turnover values $k_{cat}$ and the condition-dependent expression of enzymes (dictated by the flux $v$ these enzymes have to sustain), we can consider the proteome composition of *E. coli* under two conditions: growth medium with the complete supplement of amino acids (all 20 proteogenic amino acids present in medium), in contrast to the supplement with a single amino acid not present in the mix (a “dropout” medium) (Figure 2.8). The growth of *E. coli* in a nutrient-rich medium (glu-
### 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 are ultimately enabling 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, philosophical remarks box 3.C). The NAD(P)⁺/NAD(P)H pair 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 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 cell metabolism [69] - 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.

#### Philosophical remarks 3.C

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 [70]. 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 [71, 72, 73, 74, 70, 75].

### 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 [73] and branch points [76]. 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 [77].
(A) Methionine dropout

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

(B) Complete medium

<figure><img src="image_30.png" /><figcaption>Figure 2.8: Proteome composition of *E. coli*, grown on the growth medium with full amino acid supplement (B) or its version without amino acid methionine (A). Proteome composition data from [30]. Each polygon represents a protein, and its size indicates the amount. Proteins are grouped into sections based on their functions. You can create Proteomaps like this using https://www.proteomaps.net/.</figcaption></figure>

cose + amino acid supplement) is indeed a very fast one (with doubling time of $\tau_{d,rich} = 21.5 \pm 0.4$ vs. $\tau_{d,minimal} = 56.3 \pm 0.5$ minutes). The omission of methionine from the amino acid supplement does increase the doubling time ($\tau_{d,-Met} = 26.5 \pm 1.1$ minutes), yet the growth rate remains high, and so is the methionine biosynthesis demand in these conditions.

Methionine is an amino acid that is energetically the most expensive to make [31], and the final enzymatic reaction in the methionine synthesis pathway is so-called *rate-limiting*, or the reaction which dictates the flux through the whole pathway. Moreover, the enzyme methionine synthase (*MetE*) is a very slow enzyme (Figure 2.8, table on the bottom), thus required at large quantities to provide enough methionine for protein synthesis at high growth. Consequently, it was observed that *MetE* alone could occupy up to ca. 7.5% of the total proteome (by mass) in medium lacking methionine, and growth on a medium, containing methionine, would reduce the proteome fraction by ca. 800-fold, to 0.009% [30]. To contrast this highly condition-dependent expression of *MetE*, we considered a protein in the lower glycolysis, enolase (*Eno*) (Table 2.3). The expression of glycolytic proteins, including *Eno*, was determined to be similar, as both the complete- and the methionine-free media contained glucose as the main carbon source. A noticeable contrast of *Eno* vs. *MetE* is also a ca. 3 orders-of-magnitude higher $k_{cat}$ value compared to the one of *MetE*: having to invest less (per mass) into this enzyme contributes to the ability to sustain a very high flux through enolase when cells grow fast on glucose [30] (see Chapter 6 for a more detailed discussion).

The variable concentrations of metabolic substrates, and their relation to the enzyme parameters ($K_M$ in this case), also bring additional kinetic considerations. The above-introduced turnover value $k_{cat}$ rep-

<table><caption>Table 2.3: Abundance and $k_{cat}$ values of two selected proteins from Figure 2.8: enolase (independent on amino acid supplementation) and methionine synthase (dependent on amino acid supplementation).</caption><thead><tr><th rowspan="2">Pathway</th><th rowspan="2">Enzyme</th><th colspan="2">Proteome mass fraction (%)</th><th rowspan="2">$k_{cat}$ (s<sup>-1</sup>)</th></tr><tr><th>Met dropout</th><th>Complete</th></tr></thead><tbody><tr><td>Glycolysis</td><td>Enolase (Eno)</td><td>0.53</td><td>0.53</td><td>192.95</td></tr><tr><td>Amino acid biosynthesis</td><td>Methionine synthase (MetE)</td><td>7.45</td><td>0.009</td><td>0.12</td></tr></tbody></table>
## 3.2. Dynamics and regulation of metabolism

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 summary, metabolism constitutes a 'system' of metabolites and their reactions, together with enzymes. Its dynamics over time ensures fluxes of matter.

### 3.2.1. Biochemical reactions and thermodynamics

Metabolism consists of individual biochemical reactions of the form:

$$ \nu_a A + \nu_b B \rightleftharpoons \nu_c C + \nu_d D \quad (3.1) $$

where $\nu_i$ are the so-called stoichiometric coefficients, determining the number of molecules of the i'th chemical species taking part in the reaction (Fig. 3.5). 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. [78]) and also to books for a conceptual introduction to thermodynamics (e.g. [79]). 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 write here the Gibbs free energy of reaction directly in terms of concentrations:

$$ \Delta G_r = \Delta G_r^\circ + R \cdot T \cdot \ln \frac{c^{\nu_c} \cdot d^{\nu_d}}{a^{\nu_a} \cdot b^{\nu_b}}, \quad (3.2) $$

where the small letters indicate the concentrations of the reactants and products as given in the above reaction. Notice that specifying 'products' and 'substrates' automatically specifies a 'forward' direction to the reaction (Fig. 3.5). 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, $\Gamma$. The term $\Delta G_r^\circ$ 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 G_r^\circ < 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 that Gibbs free energy will also feature in rate equations for biochemical reactions.

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

$$ \Delta G_r^\circ = -RT \cdot \ln \frac{c_{eq}^{\nu_c} \cdot d_{eq}^{\nu_d}}{a_{eq}^{\nu_a} \cdot b_{eq}^{\nu_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}^{\nu_c} \cdot d_{eq}^{\nu_d}}{a_{eq}^{\nu_a} \cdot b_{eq}^{\nu_b}}$. Re-arranging equation 3.3, we can derive an expression for $K_{eq}$ as follows:

$$ K_{eq} = e^{\frac{-\Delta G_r^\circ}{RT}} \quad (3.4) $$

Notice that $K_{eq}$ depends only on $\Delta G_r^\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
<figure><img src="image_31.png" /><figcaption>Figure 2.9: The relation between the apparent and measured turnover value ($k_{app}$ and $k_{cat}$, respectively) - Factors, leading to low net rate of reaction per unit of protein (e.g. low substrate concentration) lead to $k_{app}$ being significantly lower than the measured $k_{cat}$ value, latter of which corresponds to the maximal rate of the reaction.</figcaption></figure>

resents the highest possible efficacy of the enzyme, where all substrates are accessible in concentration needed to sustain this efficacy (also called *saturating* concentrations). Turnover values are usually measured *in vitro*, with all the substrates highly in excess, thus deliberately minimizing many kinetic effects (enzyme saturation, reversibility of reactions, etc.) that are prevalent in more physiological conditions (see Chapter 3 for details). Therefore, what we usually observe in living cells is not the enzyme efficacy in terms of $k_{cat}$, but rather their *apparent turnover value* $k_{app}$ (Figure 2.9). The ratio of these values ($\frac{k_{app}}{k_{cat}}$) is then called the *enzyme efficiency* and can be used to infer how far the enzyme is from its optimal working conditions. The $k_{app}$ value of an enzyme *in vivo* can be computed as follows: knowing the $k_{cat}$ value, the flux through the reaction, one can calculate the minimal demand (in moles) of the enzyme to maintain that flux. Then, the $k_{app}$ value can be computed by taking the ratio between the predicted minimal enzyme demand and the abundance of enzymes in the cells.

**Macromolecule polymerization.** Moving from metabolic enzymes to macromolecular synthesis machinery, the polymerization of macromolecules (DNA replication, RNA transcription, and protein translation) is catalyzed by large enzyme complexes (and RNA, in the case of ribosomes). DNA and RNA polymerases (DNAP, RNAP) and ribosomes. The resources needed for their expression also contribute significantly to the total costs of macromolecule biosynthesis. For example, the molecular weight of an intact ribosome in *E. coli* is ca. 2.3 MDa (BNID 111560), and the *E. coli* ribosome consists of 62% RNA and 38% protein (in mass %, BNID 109047). Meanwhile, eukaryal ribosomes are even larger, ca. 3.3 MDa for *S. cerevisiae* and ca. 4.3 MDa for human *H. sapiens* (BNID 111560), and have higher protein content [32]. For a comparison, the average length of a protein in *E. coli* is ca. 300 amino acids (BNID 100017) and average amino acid weight is ca. 109 Da (BNID 104877). By multiplying these numbers, the molecular mass of an average protein is ca. 32.7 kDa, roughly 70× lower than the ribosome that synthesizes this protein.

The nature of these large complexes requires an exceptional coordination of resources. The first consideration is the number of individual proteins that form these complexes: the RNA polymerases of *S. cerevisiae* contain up to 17 subunits (BNID 111568), and 79 ribosomal proteins form a fully functional ribosome [33]. Therefore, the assembly of these complexes must be fast and robust: Thus, cells contain a number of assembly factors to facilitate these processes. Next, the coordination also has to be temporal, especially for prokaryotes, where both messenger RNA transcription and protein translation can happen simultaneously. In *E. coli*, this is well illustrated by the 3-fold difference between elongation rates of mRNAs and proteins, ca. 62 nt s⁻¹ and 21 aa s⁻¹, respectively (BNID 103021, 107868). This coordination is essential for coordinated transcription and translation happening at the same time [34], as translation happens in steps, 3 nt each (so-called *triplets*). Even in eukaryote *S. cerevisiae* we observe a similar pattern: mRNA elongation rate of ca. 30 nt s⁻¹ (BNID 103016) and protein chain elongation rate of ca. 10.5 aa s⁻¹ [35], nearly a 3× difference. Also, the polymerization of macromolecules is very tightly connected to the metabolism: different kinds of growth limitations (limiting amounts of nutrients) were shown to create bottlenecks at different stages of protein expression [36], and the optimal regulation of these processes were selected for by the evolution [37, 38].
Law of mass action for a (bio)chemical reaction

$$v_A A + v_B B \rightleftharpoons v_C C + v_D D$$
'reactants' 'products'

<figure><img src="image_22.png" /><figcaption>Figure 3.5: 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. concentrations of substrates and products) 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 derivation of this known as the "mass action law" and the rate model compatible with it is known as the "mass action rate model".</figcaption></figure>

good source of $K_{eq}$ values of many biochemical reactions is the eQuilibrator tool [80, 81].

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 reactants. 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 [82].

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' [83]. 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) [82, 83] (see Fig. 3.5 and the Appendix).
### 2.5.3 Physical proteome space

A final type of asset required for macromolecule synthesis is the physical volume in the cell. As the cells are, again, "bags of things", they possess a finite volume, thus different processes compete for available proteome volume (also called "proteome space" interchangeably). A general trend across microorganisms is that ribosomes occupy larger proteome mass fraction (in the range of 10-40% total proteome) with increasing growth rate [19, 39], with an estimated maximum in *E. coli* of ca. 55% of total proteome mass [19]. Alongside ribosomes, biosynthetic pathways also occupy a substantial share of total proteome (e.g. enzymes, required for amino acid biosynthesis occupy up to 15% of the proteome space in *S. cerevisiae* [39]). Experimentally, the optimal allocation of proteome space can be challenged by, e.g. varying expression of an unneeded (gratuitous) protein. Both for *E. coli* and *S. cerevisiae* it was shown that increasing gratu- itous protein expression directly affects the maximal growth rate on both minimal and rich media [40, 36], suggesting that the decrease in growth rate is not dependent on the nutrient status of the cell.

Numbers provided above were measured for cells, grown in minimal medium, and some of the costs we discussed – not only proteome space, but also precursors and enzymes – could be alleviated by growth in rich medium. Uptake of biosynthetic precursors usually is less costly than biosynthesis, as expression of a single type of transporter can substitute the need of expressing a biosynthetic pathway with tens of enzymes associated. Indeed, transfer of *S. cerevisiae* cells to a amino acid-rich growth medium resulted in an increase of growth rate, caused by increased proteome allocation to ribosomes, in place of the proteins of *de novo* amino acid biosynthesis [41]. In conclusion, the physical space that proteins can occupy is also an asset that the proteins are competing for, and thus the optimal allocation of the available space is key for the cells to grow in the most favorable way under specific conditions.

## 2.6 Physicochemical considerations about cells

### 2.6.1 Cell density

Most cellular parameters we discussed so far – cell size, mass, and composition – vary greatly with the cell type, growth rate, or conditions. However, one quantity does not show such variability – buoyant cell density. Buoyant density is the ratio of cell mass to volume, usually expressed as g mL⁻¹. For most organisms, prokaryotic or eukaryotic, the buoyant cell density is around 1.05-1.15 g mL⁻¹ [42, 14]. This range results from the fact that cells are 70% water which has a density of 1 g mL⁻¹ and that most dry mass is formed by proteins, which have a density of 1.2-1.4 relative to water (BNID 111208, 104272, 101502). Other components range from 1 for lipids (BNID 108142) to 1.4-2 for nucleic acids (BNID 111208). To try the calculation of bacterial density, see Problem 2.4.

For many organisms (*E. coli*, the yeast *Schizosaccharomyces pombe*, Chinese hamster ovary cells, mouse cells), cell density is constant throughout the cell cycle and at different growth rates when growing exponentially. However, it was observed to increase in stationary phase for *E. coli* and *S. pombe* [42, 43]. On the other hand, the density of *S. cerevisiae* fluctuates during the cell cycle, which might be related to a different division mode. The organisms mentioned earlier divide by binary fission – cells divide in the middle and produce two (roughly) identical daughter cells. In contrast, *S. cerevisiae* divides asymmetrically – it grows a bud that breaks away and becomes a smaller daughter cell.

Nevertheless, despite the variability, the range of the observed values is relatively small and similar for most organisms, from bacteria to mammalian cells. There are special cases where cell density deviates from the characteristic values – for example, cells with very high fat content or gas bubbles have lower densities. However, assuming the density of 1.1 g mL⁻¹ is probably a good guess unless you work with a particularly fatty or gassy cell type.

The invariability of cell density suggests that this property is highly regulated and brings us to the next question – is there an optimal density? And what are the constraints that (possibly) determine this optimum? These questions (among others) are discussed in the next section.
### 3.2.2. 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 behaviour of metabolite concentrations, will augment the stoichiometry matrix with a rate vector, v, forming a so-called system equation:

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

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 [84, 85], while for a metabolic view of variables and parameters, the article on the Control of Flux, by Kacser and Burns, offers a valuable perspective [86].

### 3.2.3. 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.6). Under these conditions the height of water remains constant, or at a steady state.

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.6). 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.6). 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 [87, 88, 89].
### 2.6.2 Physical constraints of cell growth

The living cells are constantly subject to a handful of so-called *physical constraints*, which are directly linked to the physics and the chemistry of life. Cells cannot override (evolve to bypass) these limits – only try to cope with them. Thus, sometimes these constraints are also called “hard” constraints. Notice that we consider the “hardness” of these constraints only in the space where conditions can still sustain life: some of these limitations could be relaxed by changing abiotic conditions, but would result in breakdown of biological systems.

One of the abiotic factors would be temperature; however, increased temperatures cause proteins to *denature* (lose their 3D-folded structure, thus functionality) and destabilize biological membranes. Although there are organisms, which live in extremely high temperatures (so-called *thermophiles*), as a rule of thumb, we usually consider the temperature above 393 K (120 °C) to be close to the limit of life. There is an organism known as *Strain 121* (*Geogemma barossii*) which can grow at 121 °C (hence the name), currently the highest temperature known [44]. Next, the suboptimal concentration of inorganic salts (osmolarity) or pH could also drive similar changes, disfavoring life. Here we will consider two prominent physical limits in life: the diffusion and density limits. These two limits describe two aspects of how molecules move in aqueous environments, in our case – living cells.

The diffusion limit describes the state where enzymatic catalysis is so specific and so fast that the reaction speed is determined only by the collisions of substrate molecules to the enzymes, which all result in conversions (i.e. no futile collisions) [45]. Usually, the number of futile collisions vary between 1 and $10^4$ per successful conversion, and thus having as little futile collisions as possible greatly enhances the overall rate of the reaction. Enzymes approaching (operating at) the diffusion limit are also called *perfect* enzymes. Currently there are no enzymes reported which are considerably “above” diffusion limit (see [45] for an in-depth discussion), suggesting the universality of the underlying constraint. Nonetheless, cells do have a strategy to counter the diffusion limit. Consecutive enzymes from a pathway can be placed on a scaffold, which allows the product of one reaction to be channeled directly into the next reaction without diffusing away.

Another aspect to consider is the density, or sum concentration of molecules, of the fluid. As described in previous sections, cell cytosol contains a spectrum of different molecules at different sizes and concentrations. We normally assume that some sort of optimal cell density that maximizes fitness exists, however, the density is known to fluctuate substantially in time and across conditions [46]. One of the most prevalent properties, linked to cytosolic density, is *macromolecular crowding*. As the name suggests, it describes the concentration of biological macromolecules, mainly proteins, in cytosol (thus in bacteria, the genomic DNA also contributes to molecular crowding). For example, the macromolecular crowding is suggested to impose a limit on the protein translation [47], therefore, increased crowding would result in a growth rate decrease. The state of macromolecular crowding is relevant for the cellular function, and is proposed to be in homeostasis (reviewed in [48]): optimal macromolecular crowding corresponds to a state where crowding reduces the path proteins have to diffuse, yet does not substantially decrease the speed of diffusion. In such a way, maintaining high macromolecular crowding is suggested to maximize reaction rates in the cytosol [49].

## 2.7 Concluding remarks

In this chapter, we discussed the properties and the quantities of the main cellular components, how the composition changes in different environmental conditions, and what resources are needed for a cell to replicate itself. It may seem that we already have a vast amount of data, but a lot is still missing. Most available data comes from model organisms such as *E. coli*, *S. cerevisiae*, or humans, but the data for other organisms is still limited. Single-cell data (ideally with subcellular resolution) is also not widely available. Even though we can sequence a genome within a few hours or days, we still do not know the functions of many genes. Many experiments still need to be done, and new high-throughput experimental methods developed to fill the gaps in our knowledge.

Nevertheless, with the basic knowledge from this chapter, we can dive deeper into studying cellular eco-
Steady state

A

Thermodynamic steady state

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

B

Dynamic steady states – non-equilibrium thermodynamics

<figure><img src="image_24.png" /><figcaption>Figure 3.6: Illustration of thermodynamic equilibrium and dynamical steady state. 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>

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

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

$$ v_o = k_1 \qquad S \xrightarrow{v_1=k_2 \cdot s} v_2 \qquad (3.7) $$

The differential equation describing this system is given by:

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

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

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

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

### 3.2.4. 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.7. This shows a linear pathway of two reactions, with the first reaction activated by the species $x$.
nomics and resource allocation with mathematical modeling. How is biomass represented in mathematical models? Often, models only focus on proteome as it is a cell's most abundant and expensive component. However, some models also include other major components (RNA, DNA, lipids, carbohydrates, cofactors, etc.). The components can be modeled at different levels of detail. For example, the cell proteome can be represented simply as a total proteome, divided into protein subgroups (e.g. metabolic, ribosomal, other), or modeled as individual proteins. Finally, there are two contrasting ways to include biomass in mathematical models. On the one hand, some models consider a fixed biomass composition based on measurements or literature (see Chapters 4 and 5). On the other hand, some models predict the biomass composition (i.e. they calculate optimal resource allocation or enumerate all possible compositions, see Chapter 10).

Apart from biomass composition, we can include other cellular properties as constraints or parameters in the models, depending on the type of a model and how detailed it is. For example, we can constrain the transcription/translation rates, enzyme turnover rates, cell surface area or volume.

In conclusion, this chapter introduced the basic building blocks of a cell, processes that make them, how they are coordinated and how they depend on environmental conditions. In the next chapters you will learn how to translate this information into mathematical models and how to use them to gain deeper knowledge of cell biology.

## Recommended readings

### Cell biology

* Alberts, B., et al. (2022). Molecular Biology of the Cell. WW Norton & Co.

### Numbers in cell biology

* Moran, U., Phillips, R., & Milo, R. (2010). SnapShot: key numbers in biology. Cell, 141(7), 1262-1262.

* Milo, R., & Phillips, R. (2015). Cell biology by the numbers. Garland Science.

## Problems

### Problem 2.1 Intuition for biological numbers

Try to answer the following questions, and only then look up the results:

* What is the volume of a cell?

* What is the size of a protein?

* What is bigger, a protein or the mRNA that encoded it?

* How many protein molecules are there in a cell?

* What is the number of genes in a genome?

* How long does it take to transcribe a gene?

* How long does it take to produce a protein molecule?

* What is the minimal doubling time of a cell?

* What other questions come to your mind?

Precise values do not matter here – think about orders of magnitude.

### Problem 2.2 Proteins per cell - estimate one

How many proteins are there in a bacterial/yeast/mammalian cell [14]? Use data from the following table:

<table><tr><td>Protein mass per volume</td><td>0.2 g mL<sup>-1</sup></td></tr><tr><td>Molecular mass of a protein</td><td>40000 g mol<sup>-1</sup></td></tr><tr><td>Avogadro's number</td><td>6 &middot; 10<sup>23</sup> 1/mol</td></tr><tr><td><i>E. coli</i> volume</td><td>1 &mu;m<sup>3</sup></td></tr><tr><td><i>S. cerevisiae</i> volume</td><td>60 &mu;m<sup>3</sup></td></tr><tr><td>Mammalian cell volume</td><td>3000 &mu;m<sup>3</sup></td></tr></table>
<figure><img src="image_25.png" /><figcaption>Figure 3.7: Cartoon of a simple pathway that features allosteric enzyme regulation and that can show multiple steady-state solutions (see Appendix). 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>

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 the chapter appendix (Listing A.2). 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 [85]. 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 [90, 91] and section 3.4).

### 3.2.5. Regulation of fluxes

How does the cell ‘regulate’ the flux of matter in metabolism? How does it decide, for example, to make more of an amino acid or rather more of a lipid? Or do these decisions happen automatically, through system dynamics of the metabolic system? The question of regulation of metabolism is a major research area in its own right. Several hypotheses have been formulated and some have been supported by experimental measurements. It is highly likely that many of these hypotheses are true under some conditions, and actual regulation of metabolism involves multiple mechanisms. Two of the key mechanisms we can highlight here and that we will touch upon in this and other chapters are: flux regulation through control of enzyme levels or enzyme activity. The former is achieved via control of an enzymes’ expression level, while the latter can be achieved via substrate-level allosteric regulation (Fig. 3.8). Notice that the latter case involves regulation of enzyme activity by metabolites, thereby providing a ‘dynamical regulation’ that does not require additional elements (such as gene regulatory factors). Additional examples of such dynamical regulation, which is sometimes referred to as ‘self-regulation’, can also emerge from specific pathway structures and are being proposed and explored continually, e.g. [73, 74, 70, 75]. We will discuss this topic further in the section 3.4.

## 3.3. Toolbox for modelling dynamics 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 this section, we will learn how to model one biochemical reaction, and 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 [84] and by Steven Strogatz on nonlinear dynamics [85]. 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 [92, 93, 94]
### Problem 2.3 Proteins/ribosomes per cell - estimate two

A typical protein has a volume of $25 \text{ nm}^3$ (BNID 101828) and a ribosome $3400 \text{ nm}^3$ (BNID 104919). Given that 70% of a cell volume is water, what is the maximum number of protein/ribosome molecules that fit into a typical *E. coli* cell (see Table 2.2)? Compare your answers with the previous problem/values in BioNumbers database.

### Problem 2.4 Buoyant cell density

Calculate the buoyant density of a typical bacteria using the following data:

<table><thead><tr><td>Component</td><td>Density (g mL<sup>-1</sup>)</td><td>Mass fraction per cell</td></tr></thead><tbody><tr><td>Water</td><td>1</td><td>0.7</td></tr><tr><td>Proteins</td><td>1.3</td><td>0.18</td></tr><tr><td>Nucleic acids</td><td>1.7</td><td>0.08</td></tr><tr><td>Lipids</td><td>1</td><td>0.03</td></tr><tr><td>Carbohydrates</td><td>1.5</td><td>0.01</td></tr></tbody></table>

### Problem 2.5 Concentrations enzymes and substrates

Dourado *et al.* [50] suggested that there is a relationship between the concentrations of enzymes and their substrates in *E. coli*, which is a result of a constraint on the biomass density. They showed that the reaction flux is maximal when the dry mass of each substrate is equal to the dry mass of the unsaturated (free) enzyme. What is the concentration of one enzyme per cell for *E. coli* (in $\text{mol L}^{-1}$)? What would be the optimal concentration of its substrate? Use protein mass and cell volume from Problem 2.2 and the mass of glucose as substrate.

### Problem 2.6 Cell size in different dimensions

Imagine a spherical cell that increases its diameter from 1 to $2\mu\text{m}$. How much do the surface area, volume, and SA/V change? Think about how this could influence the import of nutrients and the diffusion across the cell.

### Problem 2.7 Alien lifeforms

Imagine alien lifeforms. Would they be composed of cells? Why? What features of cells could be completely different? What features are so much dictated by physics that they could not be different in any type of alien cell?

### Problem 2.8 Substrate demand to saturate an enzyme

Take the following rate law: $v = \nu^{\text{max}} \frac{S}{K_M+S}$ (also known as irreversible Michaelis-Menten rate law, see Chapter 3), where $\nu^{\text{max}}$ is the maximal reaction velocity. Plug in the values for $v$ and compare the substrate concentration needed for the reaction rate to increase from (i) 10% to (ii) 90% of the maximal rate $\nu^{\text{max}}$. Hint: express the $S$ in terms of $K_M$ and take the ratio.
Enzymes and flux regulation

<figure><img src="image_26.png" /><figcaption>Figure 3.8: (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.3 and Appendix 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 and altering its overall reaction rate (either through competition with the substrate or by altering the enzyme structure and affecting its kinetic parameters).</figcaption></figure>

3.3.1. Enzymes - a brief note

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 [92, 95]). 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 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.8).
### 3.3.2. 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 [96]. 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 equations 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 section 1.2 above (Fig. 3.5). As discussed in that section, the "law of mass action", 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' [83, 82] (see Fig. 3.5). 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 [82]. 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^{\nu_{a}} \cdot b^{\nu_{b}} - k_{-} \cdot c^{\nu_{c}} \cdot d^{\nu_{d}}, \qquad (3.10)$$

where small letters denote concentration of the relevant species of the same letter, $\nu_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 (Fig. 3.5). 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 (Fig. 3.5). 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^{\nu_{a}} \cdot b^{\nu_{b}} - k_{-} \cdot c^{\nu_{c}} \cdot d^{\nu_{d}}$$

$$\frac{k_{+}}{k_{-}} = \frac{c^{\nu_{c}} \cdot d^{\nu_{d}}}{a^{\nu_{a}} \cdot b^{\nu_{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 G^{\circ}}{RT}} \qquad (3.11)$$

where $\Delta G^{\circ}$ is the Gibbs free energy under biological standard conditions, and R and T denote the gas constant and temperature respectively (see Fig. 3.5). 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_+$.
# Chapter 3

## Cell metabolism

Hadrien Delattre, Wolfram Liebermeister, Elad Noor, Herbert M. Sauro, Orkun S. Soyer, and Robert West

### 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' [51]. Known mostly through the names of their discoverers, these include the Entner-Doudoroff (ED), Embden-Meyerhof-Parnas (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
Following on from this last point, it is important to consider a reaction with large $K_{eq}$, i.e. a reaction for which $\Delta G^{\prime\circ}$ 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^{\nu_{a}} \cdot b^{\nu_{b}} \qquad (3.12)$$

**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 equations. 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). The reader can consult additional books (e.g. [93]) 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 studies of that time and afterwards arrived at a similar model using different assumptions. Implementation of the specific assumptions, as we detailed in the Appendix, 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.13)$$

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.9) 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 [97]. It considers
<figure><img src="image_32.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>

amino acid, nucleotides, and biomass precursors [52]. 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,
<figure><img src="image_27.png" /><figcaption>Figure 3.9: A plot depicting the Michaelis-Menten rate law. The x- and y-axis show the substrate concentration (normalised by $K_M$) and reaction flux (normalised by $v_{max}$) respectively. The dashed horizontal line corresponds to $v_{max}$, i.e. $E_{tot} \cdot k_{cat}$.</figcaption></figure>

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 equation 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 the Appendix, we can follow this strategy to derive the reversible rate equation 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}}} \quad (3.14)$$

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 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) and the flux-force relation (see below), we can then re-formulate the reversible rate equation as:

$$v = E_{tot} \cdot k_{cat}^{+} \cdot \frac{s/K_{S}}{1 + p/K_{P} + s/K_{S}} \cdot (1 - e^{\Delta G'_{r}/RT}) \quad (3.15)$$

where $\Delta G'_{r}$ 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 equation 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 G'_{r} \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 dependency of reaction rate to reaction free energy.
<figure><img src="image_33.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.

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_34.png" /><figcaption>Figure 3.2: A simplified map of central metabolism, particularly highlighting interconnections among different processes (i.e. pathways) through the NAD(P)⁺ / NAD(P)H co-substrate pair.</figcaption></figure>
Another way of writing equation 3.15 is this one:

$$v = E_{tot} \cdot k_{cat}^{+} \cdot \frac{s/K_{S} \cdot (1 - e^{\Delta G'_{r}/RT})}{1 + s/K_{S} \cdot (1 + \frac{k_{cat}^{+}}{k_{cat}} \cdot e^{\Delta G'_{r}/RT})} \quad (3.16)$$

where we replace $p/K_P$ with an expression that depends on $s$ and $\Delta G'_{r}$. This alternative expression, developed in the context of modelling microbial metabolism [98, 99], can be useful because it shows us that when the reaction is far from equilibrium ($\Delta G'_{r} \ll 0$), the term $e^{\Delta G'_{r}/RT}$ will approach zero and the above formula can be approximated by the irreversible Michaelis-Menten rate law (Equation 3.13). In this case, we further notice that the Haldane coefficient $K_s$ becomes equivalent to $K_M$ introduced above in the irreversible reaction scheme.

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 [8, 100, 80].

**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{v_{max} \cdot s}{K_{M} + s + s^{2}/K_{I}} \quad (3.17)$$

where $K_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 and discussed in the Appendix.

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 [101].

**Flux-force relationship** All (bio)chemical 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.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.10, 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^{\nu_{a}} \cdot b^{\nu_{b}}$$

$$v_{-} = k_{-} \cdot c^{\nu_{c}} \cdot d^{\nu_{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^{\nu_{c}} \cdot d^{\nu_{d}}}{k_{+} \cdot a^{\nu_{a}} \cdot b^{\nu_{b}}}\right) = v_{+} \cdot \left(1 - \frac{k_{-}}{k_{+}} \cdot \Gamma\right)$$
### 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 metabolism respectively, and their coupling presents the whole cell metabolism (Fig. 3.3 A). This coarse-grained model is widely used (e.g. [52, 53]). 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) [54, 55]. (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
<figure><img src="image_28.png" /><figcaption>Figure 3.10: The ratio of net forward flux (J) to forward reaction rate ($v_+$) as a function of the negative reaction Gibbs free energy</figcaption></figure>

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.11) 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 G'_c}{RT}} \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.10). 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 [88, 102] and holds also for the enzymatic reversible reaction model described above.

**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 law of mass action. 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 [103]. A recent study by Shen et al [104] 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
<figure><img src="image_35.png" /></figure>

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.

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 cell metabolism [57] - 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 [58] and branch points [59]. 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 [60].

## 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
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” [102]. By re-arranging the above terms, we can show that correspondence between the two models can be expressed as:

$$E_{tot} \cdot k_{cat}^{+} \cdot (1/K_{S})/(1 + s/K_{S} + p/K_{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_{tot} \cdot k_{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_{tot} \cdot k_{cat}^{+}/K_S)/(1 + (s \cdot K_S - s \cdot K_P)/(K_S \cdot K_P) + C/K_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.4. Dynamics of metabolism: Examples of 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 [70, 89]. Here, we highlight some of the illustrative experimental observations on these metabolic dynamics, and their possible model-based
<figure><img src="image_36.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 [51]. 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 [56].

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

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.

### 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. [66]) and also to books for a conceptual introduction to thermodynamics (e.g. [67]). Here, it suffices for us to define
Respiration – fermentation switch

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

Figure 3.11: 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.

explanations.

3.4.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 CO2 (from full oxidation of glucose) and water (from reduction of oxygen, cf. Fig. 3.1). 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.11) [63].

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 [105]. The respiration to fermentation shift happens
under lack of electron acceptors or with increasing growth rate [105, 106, 107, 108, 109, 110, 111]. 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 [112]. 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
<figure><img src="image_38.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_39.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 [61]. 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 [62, 63, 58, 64, 61, 65].

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 write here the Gibbs free energy of reaction directly in terms of concentrations:

$$ \Delta_r G' = \Delta_r G^{\circ} + R \cdot T \cdot \ln \underbrace{\frac{c^{n_c} \cdot d^{n_d}}{a^{n_a} \cdot b^{n_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.
energy rich compounds is not limited to fermentation. Excretion of amino acids and vitamins seem particularly common [113, 114, 115, 116], 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.13). This kind of enzyme-level regulation can be mediated through regulation of transcription factors by specific signalling molecules, including metabolites [117]. The latter case is explored in models of central metabolism [118], and it was shown that fructose-1,6-bisphosphate could act as a "flux sensor", conveying information about glycolysis flux onto key transcription factors regulating glycolytic enzymes [119].

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 [76]. 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.13), 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 [76]. 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 (e.g. [120, 121]). 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 [75]. 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.4.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 [122, 123], ATP levels [124], and expression levels of the enzymes involved in glycolysis and the TCA cycle [91].

Bistability in metabolic responses is experimentally implicated in the context of respiration to fermentation switch [125], and when carbon metabolism is initiated on glucose [126] or when it switches from glucose to other carbon sources [127, 90, 128]. 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 [127]. 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 [70, 89]. A particular
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}}, \qquad (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}} \qquad (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) [68, 69].

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 [70].

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' [71]. 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) [70, 71] (see Box 3.D and the Appendix A.2).

### 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 [72, 73]). 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 sites, where specific molecules (including sometimes the enzyme's own substrate
<figure><img src="image_30.png" /><figcaption>Figure 3.12: Cartoon representations and brief analysis results of two enzymatic models capable of bistability. Top Row: 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. Bottom Row: 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.</figcaption></figure>

'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.12). 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 [129, 130, 131, 132, 133]. This motif is found in several locations within metabolism, particularly around dehydrogenases, such as lactate dehydrogenase [131], and kinase/phosphatase pairs, such as those involved around fructose-6-phosphate [134], 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 [130, 131, 135].

It is notable that many of these models incorporated negative feedback via empirical alteration of a Michealis-Menten type reaction rate model (i.e. one of the approaches we mentioned in the paragraph above on allosteric rate models, see Fig. 3.17). This raises the question about the actual biochemical mechanisms that can lead to bistability in an 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 [101] (see Fig. 3.12). 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
<figure><img src="image_40.png" /></figure>

Physics box 3.D Mass action law for chemical reactions

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

Thermodynamic interpretation

Gibbs free energy of reaction:

$$ \Delta_{\mathrm{r}} G^{\prime} = \Delta_{\mathrm{r}} G^{\prime 0} + R \cdot T \cdot \ln \frac{c^{n_{\mathrm{c}}} \cdot d^{n_{\mathrm{d}}}}{a^{n_{\mathrm{a}}} \cdot b^{n_{\mathrm{b}}}} $$

At equilibrium:

$$ \Delta_{\mathrm{r}} G^{\prime 0} = -R \cdot T \cdot \ln \frac{c_{\mathrm{eq}}^{n_{\mathrm{c}}} \cdot d_{\mathrm{eq}}^{n_{\mathrm{d}}}}{a_{\mathrm{eq}}^{n_{\mathrm{a}}} \cdot b_{\mathrm{eq}}^{n_{\mathrm{b}}}} \\ e^{-\frac{\Delta_{\mathrm{r}} G^{\prime 0}}{R \cdot T}} = \frac{c_{\mathrm{eq}}^{n_{\mathrm{c}}} \cdot d_{\mathrm{eq}}^{n_{\mathrm{d}}}}{a_{\mathrm{eq}}^{n_{\mathrm{a}}} \cdot b_{\mathrm{eq}}^{n_{\mathrm{b}}}} = K_{\mathrm{eq}} $$

Kinetic interpretation

Backward reaction rate:

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

Forward reaction rate:

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

At equilibrium:

$$ k_{+} \cdot a_{\mathrm{eq}}^{n_{\mathrm{a}}} \cdot b_{\mathrm{eq}}^{n_{\mathrm{b}}} = k_{-} \cdot c_{\mathrm{eq}}^{n_{\mathrm{c}}} \cdot d_{\mathrm{eq}}^{n_{\mathrm{d}}} \\ \frac{k_{+}}{k_{-}} = \frac{c_{\mathrm{eq}}^{n_{\mathrm{c}}} \cdot d_{\mathrm{eq}}^{n_{\mathrm{d}}}}{a_{\mathrm{eq}}^{n_{\mathrm{a}}} \cdot b_{\mathrm{eq}}^{n_{\mathrm{b}}}} = K_{\mathrm{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).

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

### 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 [74]. Dynamical models of metabolic systems, as with all mathematical models, must therefore balance abstraction of real mech-
bistability (e.g. [136]).

### 3.4.3. Oscillations

Sustained and damped oscillations are common dynamics in nonlinear systems and can arise from a combination of positive and negative feedbacks [85]. 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 [137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147]. In the case of experiments involving cell extracts, these metabolite oscillations had a period ranging from few to tens of minutes [145, 146, 147]. In these experiments, oscillations are verified not to be due to artificial changes in ATP dynamics arising in the cell extract preparations [146], and oscillations could be entrained by controlled glucose additions [145]. 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 [134].

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

**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 [148, 149, 150].

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 [148, 149, 152, 154, 152]. 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, involving some of these proposed synchronization molecules, were also developed [153, 154, 155] and could reproduce experimental findings.

### 3.5. 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).
<figure><img src="image_42.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 A.3 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>

Anistic 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' [71, 70] (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 [70]. 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}}, \quad (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
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.13) [31], and changes in flux patterns with changing conditions cannot be explained by enzyme levels alone [121]. 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 [156, 96]. 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 [97]. Historically important book on enzyme kinetics and enzymatic reaction models.

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

* "Enzyme Kinetics for Systems Biology" (2012) by H. M Sauro [94]. 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 [92]. General introductory book on enzymes and enzyme catalysis.

### Metabolic system dynamics

* "Energy metabolism of the cell: a theoretical treatise" by J. G. Reich and E. E. Sel'kov [70]. 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 [156]. 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 [157]. Discusses biochemical network dynamics from the perspective of metabolic control analysis.

### Calculus, ODEs, and nonlinear dynamics

* "Calculus made easy" by P. T. Silvanus [84]. A historic, yet excellent, book on introduction to calculus and differential equations. If you are unsure or shy of these mathematical topics, make sure you read this book and you will get an intuitive introduction!

* "Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering" by S. Strogatz [85]. As it says on the tin - an excellent book introducing nonlinear dynamics within the applied sciences context. A brilliant book. If you don't read anything else, read the introduction chapter and be inspired!

### Thermodynamics and physical chemistry
can write;

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

where $\Delta_r G'^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'^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 A.2). The reader can consult additional books (e.g. [75]) 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 studies of that time and afterwards arrived at a similar model using different assumptions. Implementation of the specific assumptions, as we detailed in Appendix A.3, 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

<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>
* "Understanding thermodynamics" by H. C. van Ness [79]. 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 [78]. An introductory book on thermodynamics, physical chemistry, and biochemistry.

# Problems

## Problem 3.1

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

(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

(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).
<figure><img src="image_43.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}$) respectively. The dashed horizontal line corresponds to $v^{max}$, i.e. $e_{tot} \cdot k_{cat}$.</figcaption></figure>

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 [76]. 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 A.2, 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}}} \quad (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 A.2 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 A.2) and the flux-force relation (see below), we can then re-formulate the reversible rate law as:

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

where $\Delta_{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_{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
# Chapter 4

# The space of metabolic flux distributions

Daan de Groot, Felipe Scott, and David Tourigny

## Chapter highlights

* The metabolic capabilities of an organism can be related to the individual chemical reactions it can catalyze

* Elementary flux modes are minimal metabolic strategies that together span all metabolic capabilities.

* When the analysis of elementary flux modes is prohibited by computational limits, alternatives could be used, such as elementary conversion modes, flux sampling and minimal cut sets.

## 4.1. Modeling metabolic fluxes in cells

In the previous chapters we have seen that cells can convert substances from their environment into building blocks for cell components: their metabolism allows cells to grow, reproduce, repair themselves, and produce compounds needed to resist environmental stresses. But how does a cell manage this in detail, and does it have alternative metabolic strategies in case one does not function properly?

The overall metabolic conversion, for example from nutrients + oxygen to all necessary cell components and carbon dioxide, that a cell can use to grow and reproduce is in fact the consequence of many smaller chemical reactions working in concert. All chemical reactions that a cell can catalyze by expressing its enzymes form a very versatile 'metabolic network', which enables a cell to survive and grow, even when the availability of nutrients in its environment changes. There are various (semi-)automatized methods available that can be used to reconstruct this metabolic network from an organism's genome sequence, see Mendoza et al. for a review of the various methods [158]. In this chapter we will zoom in on this metabolic network and study the fluxes (reaction rates) through all individual reactions.

We call the combination of all reaction rates in a cell a 'metabolic flux distribution', and this flux distribution determines if and how a cell succeeds in taking up and converting the right nutrients to sustain itself. For a growing cell, we may ask: what will its flux distribution be, and how does this distribution change when its environment changes? Modeling metabolic fluxes allows us to answer specific questions, for instance about the change of a cell's metabolism after a gene is deleted: will it survive, and if so, will it take up different nutrients or produce different products? In contrast to the previous chapters, in the current and following chapters we are not satisfied with verbal descriptions, but seek predictive models that allow us to compute the state of a cell.

So how can we model metabolism in detail? Our main task is to describe and predict the uptake, conversion, and production of metabolites, as described by the metabolic fluxes. The rate at which a chemical reaction runs depends (through kinetics and thermodynamics) on metabolite concentrations and enzyme activities. Since enzymes are synthesized by the cell itself, the reaction rates are not only controlled by external nutrient supply, but also by gene expression. These
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_P$ with an expression that depends on $s$ and $\Delta_r G'$. This alternative expression, developed in the context of modeling microbial metabolism [77, 78], can be useful because it shows us that when the reaction is far from equilibrium ($\Delta_r G' \ll 0$), the term $e^{\Delta_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_S$ becomes equivalent to $K_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 [17, 79, 68].

**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_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 A.3.

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 [80].

**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)$$
dependencies make this a complicated field of study: the metabolic fluxes depend on the enzyme levels and metabolite concentrations, while the metabolite concentrations are again determined by the balance of fluxes through reactions that produce and consume the metabolites. In turn, enzyme levels are determined by gene expression, which is dependent on both external conditions and internal needs (e.g. the enzyme expression may change when different macromolecules need to be made in different phases of the cell cycle). To make matters even less transparent, most of the parameters (e.g. enzyme kinetic constants and details of enzyme regulation) are unknown.

So for the moment we make some simplifying assumptions in order to obtain tractable models:

* **Focus on small molecules** We focus on a subsystem of the cell, the metabolism of small molecules, which generates macromolecular precursors and energy carriers. All other processes (such as macromolecule synthesis) that happen “outside” our metabolic network are ignored.

* **Ignore spatial structure** We largely ignore the spatial structure of cells: metabolite concentrations and reaction rates are assumed to be homogeneous across the cell. The exception to this rule occurs when there are cell compartments, in which case we describe the metabolites in both compartments as if they were separate compounds (e.g. cytosolic ATP vs mitochondrial ATP), which can be converted in each other through transport “reactions”.

* **Focus on fluxes as the only variables** Instead of considering metabolite concentrations, enzyme levels and metabolic fluxes together, we will only focus on metabolic fluxes. This has important consequences for the mathematical models that we will construct: many variables, and the corresponding equations, will be ignored. Additionally, fluxes cannot be computed through enzyme kinetics, so that we need to find other, non-mechanistic ways to compute the fluxes!

* **Focus on steady-state metabolism** In a simplified picture of balanced growth (see the chapter on Balanced Growth), all metabolic processes are balanced: the rate at which material flows into the cell matches the rate at which it is converted, which again matches the production rate of macromolecule precursors. In addition, we assume that these fluxes are constant, such that the whole metabolic network is in a 'steady-state'. Taken together, we thus assume that the metabolic network can take up and produce external metabolites (e.g. extracellular metabolites and macromolecular precursors), but that all internal metabolites ('inside' the metabolic network) are mass-balanced, that is, for each of these metabolites, production and consumption cancel out.

* **Describe precursor demand by a "biomass reaction"** We assume that cell growth (or: biomass production) requires a fixed set of macromolecule precursors in fixed proportions, corresponding to the average mixture of cell components that are necessary to make a cell. For metabolism, this means that the production of more macromolecule precursors only leads to more biomass production when the production of all precursors is scaled up proportionally. We formally express this by a hypothetical "biomass reaction" that consumes a mix of precursors and energy carriers in the predefined proportions. Hence, in the metabolic models we will describe the term "biomass" has a special meaning: while it usually means "the totality of compounds in a cell", here we use it for "the totality of compounds *outside* our metabolic model, which metabolism needs to produce".

* **Ignore dilution of small molecules** When a cell doubles its size but does not produce a certain metabolite, the concentration of this metabolite will halve. This basic principle is called 'dilution by growth', and in principle affects all compounds in the cell. During balanced growth, the production of macromolecules that are produced but not degraded should balance dilution, i.e. the number of each macromolecule should double when the cell doubles its size. This requires the rate of precursor supply to match the dilution rate, and hence the cell's growth rate. Similarly, small molecules are diluted, but since these are also degraded by consuming reactions, the rate of dilution is usually negligible compared to the production and consumption by metabolic reactions. Therefore, the models below will usually ignore the dilution of such metabolites.

* **Resource allocation** Since each enzyme has a maximal catalytic rate (the $k_{cat}$ value), a reaction flux will require a certain (minimal) amount of enzyme, which takes up cellular space; since cellular space is limited, fluxes cannot increase infinitely since there is always an upper bound on a weighted sum of reaction fluxes. This constraint implies compromises between different reaction fluxes: one flux can only be increased at the expense of others.

With these assumptions, we are converging on a mathematical model: we know which variables to describe (the metabolic
<figure><img src="image_44.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>

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 [81, 82] 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 [83]. A recent study by Shen et al [84] 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
fluxes in steady-state metabolism), which constraints to apply (the balance of production and consumption of all internal metabolites) and what main input information we need (the metabolic network, described by a list of chemical reaction equations). Importantly, the model will be able to describe compromise: for example, with a given carbon influx and assuming mass balance, the carbon atoms can either be used to generate energy or biomass; if one function increases, the other one goes down. To obtain realistic predictions, we may introduce additional constraints, for example known flux directions or experimentally measured uptake rates. All this information will not suffice to predict metabolic fluxes precisely, but it allows us to narrow down the possible flux distributions. Importantly, all formulae in these models are linear, which makes them tractable even for very large model sizes (with thousands or even hundreds of thousands of variables).

Notably, all these assumptions depend only on the list of chemical reaction equations (the stoichiometry of the metabolic network), and nothing needs to be known about enzyme kinetics. So if the networks are already known, what do we gain from this kind of modeling? Even if a metabolic network structure is known reaction by reaction, this does not mean that we understand the network-wide behavior, i.e. which overall flux distributions are possible, and what overall flux distributions are useful for the cell. Our aim here is to make the step from structural information (about the network) to physiological insights about how the network can be used. We can learn, for example, how much biomass can be made from a certain amount of glucose, and whether an enzyme deletion is lethal because a certain precursor cannot be produced anymore.

Metabolic network structures (in the form of stoichiometric matrices) are approximately known for many microbial species, and to some extent for higher organisms. Together with the constraints outlined above, this network determines a range (or “space”) of possible flux distributions. In this chapter we will characterize this space of possible flux distributions according to our assumptions, and since we characterize fluxes entirely by constraints the models will be called “constraint-based models”. We will get to know mathematical tools to characterize this space in a simple way: for instance, to describe all possibilities that a metabolic network provides we can use Elementary Flux Modes (EFMs).

In the next chapter, we will combine such constraint-based models with optimality principles: out of the space of possible flux distributions, specific “optimal” flux distributions will be selected because these are supposedly “most profitable”, either for the cell or for metabolic engineering purposes. Some of the flux prediction methods that we will describe refer also to concentrations; for instance, metabolite concentrations play a role in thermodynamic constraints that exclude certain flux directions, and enzyme concentrations come into play in models that associate fluxes with an enzyme demand. However, in all cases, the connection between fluxes and concentrations is very simple, and real enzyme kinetics are ignored. In later chapters, we will then see how the models change when more and more of the complex details are added about metabolite concentrations, enzyme kinetics, and thermodynamics.

### 4.2. The flux cone

#### 4.2.1. Mass-balance constraints

As described in the introduction, our models will be built on the metabolic network of all chemical reactions that an organism can catalyze. We can conveniently summarize all these chemical reactions as an $(m \times n)$-dimensional stoichiometric matrix $N$ where each of the $m$ rows corresponds to a metabolite at steady state and each of the $n$ columns to a given reaction. The entry $N_{ij}$ is the coefficient of the $i$-th metabolite in the $j$-th chemical reaction equation. Then, we can gather all $n$ reaction rates in an $n$-dimensional flux vector: $v = (v_1, \dots, v_n)^T$. This is convenient because the multiplication $Nv$ now captures the net production and consumption of all $m$ metabolites at this flux distribution, and is therefore equal to the time derivative of the metabolite concentrations: $\dot{c} = Nv$. Therefore, the steady-state assumption combined with the assumption that dilution of metabolites due to growth is negligible, can be mathematically captured in a set of linear equations that we call the *mass-balance constraints* on $v$:

$$ \dot{c} = Nv = 0. \qquad (4.1) $$
for an 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_{\text{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_{\text{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_{\text{tot}} \cdot k_{\text{cat}}^{+} \cdot (s/K_{\text{S}})/(1 + s/K_{\text{S}} + p/K_{\text{P}})$$

And, for the reversible non-enzymatic case:

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

Where $k_{\text{cat}}$, $K_{\text{S}}$, and $K_{\text{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" [82]. 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_{\text{S}}$ and $p \ll K_{\text{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_{\text{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_{\text{S}} = K_{\text{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
Since in a typical metabolic reaction network the number of metabolites is smaller than the number of reactions ($m < n$), the equations for $v$ are under-determined. This means that there are infinitely many solutions, $v$, that satisfy the mass-balance constraints. The space of all such $v$ is called the *null space* of $N$.

In the absence of any additional constraints on $v$, each $v_i$ can take on both positive and negative values, where a negative value would mean that the reaction runs in the reverse direction. However, it will often be more intuitive to think of reaction rates as positive quantities, for example when we want to deduce necessary enzyme levels from the reaction rates by assuming that enzyme levels scale are directly proportional to the catalysed reaction rate: $v_i \propto e_i$. Therefore it is often convenient to replace each reversible reaction by a forward irreversible reaction $v_i^+$ and a backward irreversible reaction $v_i^-$. Mathematically, we thus introduce non-negative variables $v_i^+, v_i^- \ge 0$ such that $v_i = v_i^+ - v_i^-$. The mass-balance constraints in these new variables become

$$0 = Nv = Nv^+ - Nv^- = (N - N) \begin{pmatrix} v^+ \\ v^- \end{pmatrix} \quad (4.2)$$

where $v^+ = (v_1^+, \dots, v_n^+)^T$ and $v^- = (v_1^-, \dots, v_n^-)^T$, respectively. The mass-balance constraints (4.2) combined with the property that $v_i^+, v_i^- \ge 0$ can be expressed in the form

$$A \begin{pmatrix} v^+ \\ v^- \end{pmatrix} \ge 0 \quad (4.3)$$

where

$$A = \begin{pmatrix} N & -N \\ -N & N \\ I & 0 \\ 0 & I \end{pmatrix}.$$

In this form the set of constraints on $(v^+, v^-)^T$ define a *polyhedral cone* and from the condition $v_i^+, v_i^- \ge 0$ we see that the cone is also *pointed*, meaning it contains no complete line and the zero vector is the only vertex (extreme point) of the cone (see Figure 4.1 for an illustration). The space of solutions that satisfies (4.3) is called the *flux cone*. It is important to remember that we only get a pointed polyhedral cone because we have chosen a representation where all reactions are irreversible; it is certainly possible to do useful analyses in the original space with reversible reactions $v$. In what follows however, we will simplify notation by identifying $v$ with $v = (v^+, v^-)^T$ and use $N$ in place of $(N, -N)$ with the implicit understanding that all components of $v$ are non-negative and $N$ accounts for all forward and reverse reactions that exist in the network. In this notation the flux cone is defined as the space

$$\mathcal{FC} = \{v \mid Nv = 0, v \ge 0\} \quad (4.4)$$

where notation $v \ge 0$ demands that each component of $v$ is non-negative.

To provide a concrete example, we consider the simple representation of central carbon metabolism presented in Figure 4.2. In this example there are four external metabolites, $G_{ex}$, $O$, $P_1$, $P_2$ and two internal metabolites: $G_{in}$ and $P$. In our model we only require mass-balance for internal metabolites, such that the steady-state constraint can be written as

$$Nv = \begin{pmatrix} 1 & -1 & 0 & 0 \\ 0 & 2 & -1 & -1 \\ 0 & 0 & -1 & 1 \end{pmatrix} \begin{pmatrix} v_0 \\ \vdots \\ v_4 \end{pmatrix} = 0, \quad (4.5)$$

where each column thus corresponds to one of the five reactions and where the rows correspond to $G_{in}$ and $P$ respectively. The 1 in the first row of the first column thus corresponds to the import of one glucose molecule $G_{in}$. In a small example like this, it is still tractable to write out all separate steady-state equations:

$$0 = v_0 - v_1, \quad 0 = 2v_1 - v_2 - v_3 + v_4. \quad (4.6)$$
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 [85] and by Steven Strogatz on nonlinear dynamics [86]. 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 [72, 75, 87].

### 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 is a *m* × *n* 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 concentrations, 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 [85, 86], while for a metabolic view of variables and parameters, the article on the Control of Flux, by Kacser and Burns, offers a valuable perspective [88].

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

It is important to note that the thermodynamic equilibrium mentioned above is also a type of steady-state,
<figure><img src="image_31.png" /><figcaption>Figure 4.1: A pointed polyhedral cone,</figcaption></figure>

These two mass-balance constraints combine with the non-negativity conditions $v_0, v_1, v_2, v_3, v_4 \ge 0$ to define the flux cone as the space of all flux vectors, $v$, that satisfy all of these constraints simultaneously.

### 4.2.2. Elementary flux modes

Although (4.4) already gives a mathematical description of the flux cone, we will here derive a more useful characterization of this space. One of the main problems of the description in (4.4) is that it does not give us a method to generate (or express in mathematical terms) a steady-state flux distribution, even though it makes it easy to check that any $v$ is in $\mathcal{F}C$. Below, we will instead introduce an exhaustive set of generators: minimal flux distributions that can be combined to make all possible flux distributions in $\mathcal{F}C$, called *elementary flux modes* (EFMs). One can think of these EFM as minimal building blocks that generate the flux cone, similar to how basis vectors generate a linear space, but with a particular rule for combining them as we describe below.

An important property of pointed polyhedral cones such as $\mathcal{F}C$ is that there exists a unique minimal set of $n$-dimensional *generators* $\{e^1, \dots, e^K\}$ such that $\mathcal{F}C$ can also be represented as

$$ \mathcal{F}C = \left\{ \mathbf{v} \mid \mathbf{v} = \sum_{k=1}^{K} \lambda_k e^k, \lambda_k \ge 0 \forall k \right\}. \qquad (4.7) $$

We must remark here that the *generators* $e^k$ are only defined up to scalar multiplication, i.e. any $\alpha e^k$ with $\alpha > 0$ could replace $e^k$ in the set of *generators*. Each *generator* $e^k$ represents one of the $K$ *extreme rays* or “edges” of the pointed cone, and for metabolic reaction networks they turn out to have a particularly useful biological interpretation, which we shall see shortly. In essence, property (4.7) says that any flux vector $v$ in the space $\mathcal{F}C$ can be expressed as a conical combination of *generators* $\{e^1, \dots, e^K\}$.

To illustrate the concept of EFM, we return to the simple representation of central carbon metabolism presented in Figure 4.2 that is associated with the stoichiometric matrix

$$ \mathbf{N} = \begin{pmatrix} 1 & -1 & 0 & 0 & 0 \\ 0 & 2 & -1 & -1 & 1 \end{pmatrix} \qquad (4.8) $$

and corresponding flux vector $v$ with components $v_0, v_1, v_2, v_3, v_4 \ge 0$. A set of EFM that serve as *generators* for the
(A) Thermodynamic steady state

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

(B) Dynamic steady states – non-equilibrium thermodynamics

<figure><img src="image_46.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>

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 [89, 81, 90].

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

$$ v_0 = k_1 \qquad S 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) $$
<figure><img src="image_32.png" /><figcaption>Figure 4.2: A simple representation of the metabolic reaction network for central carbon metabolism. Extracellular glucose, $G_{ex}$, is imported into the cell via reaction with flux $v_0$ and converted via intracellular glucose, $G_{in}$, to pyruvate, $P$, via the reaction with flux $v_1$ that has a stoichiometry coefficient of two pyruvate molecules to each glucose molecule. Pyruvate can then either be converted to a fermentation product, $P_1$, via the reaction with flux $v_2$ or, in the presence of oxygen, $O$, converted to an oxidative phosphorylation (OXPHOS) terminal product $P_2$ via the reaction with flux $v_3$. The fermentation product $P_1$ can also be converted back to pyruvate via the reaction with flux $v_4$.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Figure 4.3: EFMs $e^1, e^2, e^3$ overlaid on the simple metabolic reaction network in Figure 4.2 with concentrations of intracellular glucose and pyruvate assumed to be at steady state. It can be seen from our understanding of central carbon metabolism that $e^1$ represents the glycolytic fermentation pathway, $e^2$ the oxidative metabolism of the fermentation product, and $e^3$ the oxidative metabolism of glucose.</figcaption></figure>

resulting flux cone is given by

$$e^1 = \begin{pmatrix} 1 \\ 1 \\ 2 \\ 0 \\ 0 \end{pmatrix}, \quad e^2 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \\ 1 \end{pmatrix}, \quad e^3 = \begin{pmatrix} 1 \\ 1 \\ 0 \\ 2 \\ 0 \end{pmatrix}, \qquad (4.9)$$

and these are depicted in Figure 4.3. In fact, because we split the reversible reactions into a forward and backward reaction, the combination of reactions $v_2$ and $v_4$ would also be an EFM, but we discard such cycles created by splitting reversible reactions. We see from our understanding of central carbon metabolism that these three EFMs represent the fundamental metabolic pathways of glycolytic fermentation ($e^1$), oxidative metabolism of the fermentation product ($e^2$), and oxidative metabolism of glucose ($e^3$). The definition of the flux cone in terms of EFMs as in (4.7) is in this example equivalent to the statement that any flux vector $v$ in our toy carbon metabolism network can be viewed as a (non-negative) weighted combination of these three minimal metabolic pathways. This interpretation remains true for any metabolic reaction network: *elementary flux modes represent the minimal metabolic pathways through the metabolic reaction network at steady state*.

We can make this more precise with a mathematical definition. First, let us introduce the *support* of a vector as $\text{supp}(v) = \{i | v_i \neq 0\}$, i.e. the support of a flux vector is the set of reactions that have a non-zero rate.

**Definition 1.** A vector $v$ is an EFM if and only if it satisfies the following two properties:
<figure><img src="image_47.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 A.2). 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>

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 A.3). 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 [86]. 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 [91, 92] and section A.1).

## 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 inter-connections 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 an enzymatic irreversible reaction rate model (introduced in Eq. (3.8)) [24], and changes in flux patterns with changing conditions cannot be explained by enzyme levels alone [93]. These findings lead to the question on what determines/limits reaction fluxes and how reaction fluxes are regulated besides regulation