Cell metabolism

Cell metabolism Hadrien Delattre1, Wassili Dimitriew2, Wolfram Liebermeister3, Elad Noor4, Herbert M. Sauro5, Orkun S. Soyer1, and Robert West1 1School of Life Sciences, University of Warwick, Gibbet Hill Road, Coventry, UK; 2Bioinformatics, Friedrich-Schiller University Jena, Germany; 3Université Paris-Saclay, INRAE, MaIAGE, 78350 Jouy-en-Josas, France; 4Department of Plant and Environmental Sciences, Weizmann Institute of Science, 76100 Rehovot, Israel; 5University of Washington, Seattle, WA, US; 6Université Paris-Saclay, INRAE, MaIAGE, 78350 Jouy-en-Josas, France. Abstract Metabolism is a dynamical process, in which a myriad of chemicals are converted via biochemical reactions inside the cell over time. This chapter gives a brief overview of cell metabolism, motivates the dynamical nature of metabolism and introduces mathematical approaches for modeling it. Basic models of biochemical reaction rates are introduced and their assumptions and relations to each other are explained. Using these rate models, simple metabolic systems are modeled and their dynamical behavior analyzed. Experimental findings on metabolic dynamics are summarized and possible explanatory models of such dynamics are discussed. Keywords: cell metabolism – metabolic pathway – enzyme kinetics – metabolic network – biochemical reactions – steady state – catabolism – anabolism – redox reaction – electron flow – Michaelis-Menten rate law – thermodynamic equilibrium Contributions: The chapter was originally drafted and written by Orkun Soyer, Hadrien Delattre, and Robert West. Later, sections were added by Wassili Dimitriew, Wolfram Liebermeister, Elad Noor, and Herbert Sauro. We thank Francesco Moro, Markus Arthur Köbis, and Maarten Droste for comments on an earlier version of this chapter. We acknowledge funding from the UK’s Biotechnology and Biological Sciences Research Council (BBSRC) grant ID BB/T010150/1 and from the Gordon and Betty Moore Foundation grant ID GBMF9200. To cite this chapter: Hadrien Delattre, Wassili Dimitriew, Elad Noor, Herbert M. Sauro, Orkun S. Soyer, Robert West. Cell metabolism (Version January 2026). doi: 10.5281/zenodo.8156823. Chapter from: The Economic Cell Collective (2026). Economic Principles in Cell Biology. No commercial publisher | Online open access book | doi: 10.5281/zenodo.8156386 The authors are listed in alphabetical order. This is a chapter from the open textbook “Economic Principles in Cell Biology”. Free download from principlescellphysiology.org/book-economic-principles/. Lecture slides for this chapter are available on the website. c 2026 The Economic Cell Collective. Licensed under Creative Commons License CC-BY-SA 4.0. An online open access book. No publisher has been paid. doi: 10.5281/zenodo.8156386 1 Chapter overview ◦Metabolism is a dynamical process, in which a myriad of chemicals are converted via biochemical reactions inside the cell over time. ◦This chapter gives a brief overview of cell metabolism, motivates the dynamical nature of metabolism and introduces mathematical approaches for modeling it. ◦Basic models of biochemical reaction rates are introduced and their assumptions and relations to each other are explained. ◦Using these rate models, simple metabolic systems are modeled and their dynamical behavior analyzed. ◦Experimental findings on metabolic dynamics are summarized and possible explanatory models of such dynamics are discussed. 3.1. Cell metabolism The previous chapter gave an overview of all the cell components and their amounts in a cell. The goal of this chapter is to describe interactions between these components and the resulting dynamics. Living systems need to constantly maintain themselves, replacing parts that degrade over time, or make new ones in order to grow. We will take a closer look at how all these components are being produced. Cell metabolism is the umbrella term for all the processes involved in the constant inflow of compounds taken up and converted into other compounds through a large network of biochemical reactions. We will zoom in on this important part of the cell and show how it can be described by models. 3.1.1. What is cell metabolism? Metabolism encompasses the totality of biochemical reactions that occur within a living organism, enabling it to maintain life, and underlies the growth, reproduction, repair, and adaptation capabilities of all living organisms. For example, an E. coli bacterium has an overall elemental composition of approximately C4H6O2N. However, these atoms form thousands of different biomolecules or "metabolites", which have to be imported or synthesized, converted into each other, and possibly exported or degraded. The system achieving all this, cellular metabolism, consists of an integrated network where thousands of reactions are coordinated to achieve metabolic functions such as building new cellular mass or storing energy in the form of ATP, as well as responding to environmental cues. These reactions are often subdivided into groups of sequential steps called metabolic pathways that manage the conversion of nutrients into energy and building blocks for biosynthesis, and the generation of molecules necessary for cellular signaling and other functions. At the molecular level, metabolism supports the synthesis of the molecules that form the basis of life. These building blocks have been described in detail in Section 2.2.2 in [1] Biological molecules, and we only provide a brief recap here. The molecular components of a cell include lipids, amino acids, and nucleotides. Lipids serve as key components of cell membranes, bilayer structures that form the boundaries of cells and organelles, and participate in signal transduction and energy storage. Twenty types of amino acids can be assembled into proteins, which perform a vast array of cellular tasks, including catalysis, structural support, signaling, and transport. Nucleotides are the building blocks of DNA and RNA, which carry genetic information and play central roles in transcription and translation, the two main steps in producing proteins. To synthesize these different components and to maintain continuous growth and essential cellular processes - like the production, modification and degradation of enzymes, formation of temporary organelles, reaction to external stimuli etc. - cells require an uptake of nutrients. These nutrients, including carbohydrates, amino acids, fatty acids, and some inorganic molecules, can be made and interconverted via complex, regulated pathways. Catabolic pathways break down nutrients to release energy and generate intermediate metabolites, while anabolic pathways use these 2 catabolism anabolism carbohydrate / fat / protein byproducts NADPH NADP+ ATP ADP building blocks respiration NADH NAD+ growth H2O O2 Substrate-level phosphorylation NADH NAD+ NADPH NADP+ ATP ADP Figure 3.1: Scheme of central carbon metabolism and cofactor balance – In catabolism, a carbon source in the form of sugar, another carbohydrate, lipid (fat), or protein is degraded by catabolic reactions in order to regenerate cofactor molecules such as ATP, NADH and NADPH, from their paired forms ADP, NAD+, and NADP+, and producing important building blocks as well as byproducts such as CO2, ethanol, lactate, or acetate. The building blocks and cofactor molecules are then used by anabolic reactions to produce cellular components needed for growth. Cellular respiration, which oxidizes NADH, is the most effective path for replenishing ATP, but typically requires molecular oxygen to serve as the electron acceptor. However, substrate-level phosphorylation, i.e. when ATP is regenerated independently of oxygen by breaking up sugars, is an important alternative that allows some cells to grow anaerobically. intermediates to build complex macromolecules (see Figure 3.1). Regulation of these pathways allows cells to allocate their resources dynamically, responding to internal requirements and external changes, thus making metabolism a central determinant of cellular physiology. If we zoom into any of these metabolic pathways, we’ll see that they are comprised of simpler biochemical reactions, each involving one or more substrates, which are converted into products through the catalytic action of enzymes (see Box 3.A). The product of one reaction serves as a substrate for subsequent reactions, creating a cascade of molecular transformations in which each reaction can influence multiple other processes. In addition, there are sometimes alternative pathways that enable the cell to reroute fluxes under different physiological conditions. Put together, all these reactions synthesize the building blocks and energy carriers for the production of macromolecules. Carbon atoms are derived from a variety of sources, usually sugars, and are funneled into biosynthesis. Nitrogen-containing intermediates are used in the production of amino acids and nucleotides. Once synthesized, the building blocks are polymerized and assembled into functional macromolecular structures: lipids form membranes, amino acids form proteins, and nucleotides assemble into nucleic acids (DNA and RNA). This molecule synthesis from a set of small building blocks connects small-molecule metabolism to the production of macromolecules, and further to complex cellular structures. 3.1.2. Main metabolic pathways of cells – central carbon metabolism Metabolic networks consist of thousands of interconnected biochemical reactions. Since their overall architecture may be overwhelming, it is helpful to focus on smaller network regions with specific functions. Quite generally, one may distinguish between catabolism, comprising processes that typically break down larger molecules to smaller ones, and anabolism, comprising biosynthetic processes that typically build larger molecules from their smaller building 3 Biology box 3.A The machinery behind metabolism Enzymes Proteins and protein complexes that act as biochemical catalysts are called enzymes. Without them, most metabolic reactions would proceed too slowly to sustain life. By lowering the activation energy of a reaction, enzymes enable it to proceed more rapidly, even under mild physiological conditions. Moreover, enzymes can also serve as regulators within metabolic pathways: by modulating their activity, metabolic conversion rates can be increased or decreased, allowing the cell to adapt to metabolic demands and environmental conditions. The energy carrier ATP Adenosine triphosphate (ATP) and Adenosine diphosphate (ADP) are two nucleotides with the same basic structure but different numbers of phosphate groups. When an ATP molecule is hydrolyzed into ADP, one phosphate group is released together with about 28 kJ/mol of free energy (BNID: 100777), which can be used to drive other reactions. Aside from anabolic reactions, which are discussed in this chapter, processes such as substrate phosphorylation and ion transport - which, for example, is needed for muscle contraction and nerve impulse propagation in animals - rely on ATP as an energy source. However, other nucleotide triphosphates are of much importance as well. For example, GTP (Guanosine Tri-Phosphate) is an important energy source for protein biosynthesis. Nevertheless, ATP is the most abundant energy carrier in most organisms that exist today and is often called the energy currency of a cell. And like most currencies, cells must earn back whatever they spend and maintain a balanced budget, i.e. constant ATP and ADP concentrations. Therefore, cells are required to continuously perform the reverse process of replenishing their ATP from ADP by adding a phosphate group. Since the conversions between them in the entire network must be balanced, overall and on average, the fluxes of different pathways - including ATP-consuming processes outside of metabolism - are effectively coupled. Reduction equivalents NADH and NADPH In many chemical reactions, electrons can be released (oxidation) or added (reduction) to one of the substrates. From an electron flow perspective, aerobic metabolism involves the transfer of electrons from reduced substrates (e.g., glucose) to terminal electron acceptors (e.g., oxygen), via various redox reactions. Electron transfer between these molecules is essential for ATP production and redox homeostasis. Nicotinamide Adenine Dinucleotide (NADH) is an important molecule that acts as a mediator for transferring these electrons in many reactions. It exists in two forms - oxidized (NAD+) and reduced (NADH) - which form a cofactor pair, similar to ATP and ADP. Other redox carriers exist and fulfill similar roles, for example NADPH (a phosphorylated version of NADH). Like in the case of ATP, the interconversion NADH ↔ NAD+ (and NADPH ↔NADP+) in the two directions must be balanced across the network. Cofactors Cofactor pairs, such as ATP/ADP, NADH/NAD+, or NADPH/NADP+ are involved in many reactions and in different parts of the network. Being interconverted in multiple pathways, they couple metabolic fluxes in these pathways across the entire network. blocks (see Figures 3.1 and 3.2). But each of these processes can be divided into even smaller parts, called metabolic pathways, that perform specific chemical transformations. A metabolic pathway is a series of sequential, enzymecatalyzed reactions that convert an initial substrate (or sometimes, several substrates) into one or several final products, while (potentially) interconverting different cofactors on the way. Pathways are typically linear, but some are branched or cyclic. They are subject to regulation to ensure efficient resource usage and to meet the cell’s needs in varying conditions. Central metabolism contains prominent pathways such as glycolysis, the TCA cycle, the pentose phosphate pathway, the various amino acid synthesis pathways, and fatty acid synthesis. ◦Glycolysis [2] is a central metabolic pathway that breaks down a glucose molecule into two pyruvate molecules. Occurring in the cytoplasm, glycolysis consists of ten enzyme-catalyzed steps and generates a net yield of 2 ATP and 2 NADH molecules per glucose molecule. It also provides essential intermediates for other pathways and is the primary route of carbohydrate catabolism. ◦The tricarboxylic acid (TCA) cycle, also known as the Krebs cycle or citric acid cycle oxidizes acetyl-CoA derived from carbohydrates, fats, and proteins into CO2. The cycle produces high-energy molecules: NADH, NADPH, 4 carbohydrates sugars glycolysis glucose glyceraldehyde-3-P pyruvate acetyl-CoA TCA cycle oxidative phosphorylation amino acids fatty acids Figure 3.2: Pathways in central metabolism – Catabolic and anabolic parts of central metabolism and GTP, which are essential for ATP generation via respiration, in an oxygen-consuming process called oxidative phosphorylation. In eukaryotes, the cycle operates in the mitochondrial matrix, while in prokaryotes it operates in the cytoplasm. In prokaryotes (bacteria and archaea), all compounds are contained in a single volume called the cytoplasm, except for those that are part of the membrane and maybe a cell wall around it. Eukaryotic cells, in contrast, contain spatial compartments called organelles enclosed by membranes, which allow metabolic activities to be segregated. For example, most of the catabolism operates in the cytosol while respiration takes place in the mitochondria. Transport proteins facilitate the movement of metabolites and proteins across membranes, ensuring the coordination of distributed processes in the cell. 3.1.3. Metabolic strategies in central carbon metabolism Glycolysis and TCA cycle enable cells to produce ATP in different ways, called fermentation and respiration (see Figure 3.5). When a cell changes its metabolic strategy from fermentation to respiration, metabolic fluxes are rerouted. But there are also changes of strategy in which a metabolic flux flips its direction. An example is glycolysis, which can also run in the other direction. This reverse process, called gluconeogenesis, is necessary when cells feed on carbon sources that enter the pathway from the other end - for example, it provides a way to synthesize glucose instead of breaking it down, or to convert pyruvate into nucleotides needed as RNA or DNA precursors. 5 D-fructose-1,6P D-xylulose-5P D-sedoheptulose-7P D-ribulose-5P D-glucono-lactone-6P D-gluconate-6P glycerone-P D-glucose L-malate succinate citrate D-isocitrate cis-aconitate fumarate D-glycerate D-xylose  D-glycerate-2P D-gluconate  glyoxylate 2-keto-3-deoxy-gluconate-6P D-lactate D-xylulose acetaldehyde ethanol D-glycerate-1,3BP D-fructose-6P D-glucose-6P D-glyceraldehyde-3P D-ribose-5P D-erythrose-4P D-glycerate-3P phosphoenol pyruvate pyruvate oxaloacetate acetyl-CoA succinyl-CoA 2-ketoglutarate acetyl-P acetate glycerol glycerol-P formate 17% 3% 1% 1% 10% 10% 3% 13% 18% 20% 4% * * * zwf pgi pfk fbp glk pgl gnd rpi rpe tkt tal tkt fba gap pgk tpi gpm eno pykF pdh gntK glpK glt ppc acn acn mdh fum sdh sucCD icd sucAB ppck eda mae pps aceA glcB,aceB eda edd frd glpD ldhA xylA xylB pfl adh adh pta ackA acs garK  ATP   ADP   ATP   ADP  Pi  ATP   ADP  CO2  ATP   ADP   CoA   CoA   CoA, ATP  ADP, Pi Pi  ATP   AMP, Pi   ADP   ATP   ADP   ATP   ATP   ,ADP   Pi   ATP   ADP   CoA   CoA   CoA   Pi   ATP   ADP   ATP, CoA  AMP, PPi  ATP   ADP  CO2, 2e2e2e2e2e2e2e2e2e2e-  CoA  CO2, 2eCO2, 2eCO2, 2eCO2 CO2, 2eCO2 Figure 3.3: A map of core 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. The choice between the two flux directions – glycolysis or gluconeogenesis – depends on the substrate and product concentrations, and on the upregulation and downregulation of certain enzymes with different usage of cofactors. The change in substrate and product concentrations, as well as the different usage of cofactors, will change the direction of the thermodynamic driving forces along the pathway, which is necessary to flip the flux direction. 6 Figure 3.4: Central metabolism in E. coli bacteria. Aside from core metabolism, as shown in Figure 3.3 and here in the central upper part of the map, this map also contains the biosynthesis pathways for all main macromolecule building blocks: amino acids (for proteins), nucleotides (for RNA and DNA), and fatty acids (for lipids). Colors show a metabolic flux distribution for aerobic growth on glucose. The network represents a computational model called iCH360 [3], from which the flux distribution was computed using parsimonious FBA (see Chapter 5 in [1]). 3.1.4. Enzyme kinetics and regulation Enzyme kinetics describe how reaction rates depend on the concentrations of substrates, products, and other small molecules that act as regulatory effectors. Typically, the rates of enzymatic reactions are proportional to the enzyme level. In general, higher substrate concentrations increase a reaction rate while higher product concentrations may even decrease it. Effector molecules, called activators or inhibitors, increase or decrease the enzyme activity, respectively. Enzyme kinetics are often modeled using the Michaelis-Menten equation which already appeared in the previous chapter and is described in more detail below. A property of this rate law - which it shares with many others - 7 Respiration (high yield, low rate) vs. fermentation (low yield, high rate) Substrate Respiration TCA (Krebs) cycle Glucose Glycolysis Pyruvate … Fermentation NADH NAD+ ETC (respiration) ATPase ADP ATP akg NH4 Glutamate O2 H2O Ferment. Figure 3.5: Fermentation and respiration in cells – In fermentation, pyruvate is converted into organic acids or alcohols (e.g., lactate, ethanol), regenerating NAD+ from NADH. Fermentation does not use oxygen and allows cells to produce ATP via glycolysis only, for example when oxygen is unavailable. Cellular respiration happens by running the TCA cycle (which produces reduction equivalents) and oxidative phosphorylation, in which electrons from NADH and FADH2 flow through the electron transport chain to oxygen, the terminal electron acceptor, driving ATP synthesis via ATP synthase. In eukaryotes, oxidative phosphorylation takes place in the inner mitochondrial membrane. is that it describes enzyme saturation: as the substrate level increases, the enzyme becomes increasingly saturated with substrate and the reaction rate approaches a maximum (called Vmax). Microscopically, a positive reaction rate can be seen as the difference of a rate in forward direction (from substrates to products) and a smaller rate in backward direction (from products to substrates). The interplay between these rates depends on thermodynamics. A consequence of this is the Haldane relationship, which links kinetic parameters to the reaction equilibrium constant, thus providing thermodynamic constraints on reaction reversibility and efficiency. Metabolic dynamics itself tends to be stabilizing. Due to reaction kinetics and thermodynamics, a momentary increase of a metabolic concentration will lead to a decreased production and increased consumption of the metabolite, which has a general stabilizing effect on metabolic states. However, the rates of enzymatic reactions can also be regulated in several other ways. The regulation mechanisms include: ◦Direct enzyme regulation involving small molecules that interact with enzymes. Competitive inhibitors, which chemically resemble a substrate, can bind to the enzyme’s active site without forming a covalent bonds. By preventing substrate molecules from binding, they slow down the reaction. Allosteric regulator molecules bind to other parts of the enzyme, causing conformational changes that modulate the enzymes activity. ◦Post-translational protein modifications such as phosphorylation - the chemical addition of phosphate groups - can activate or inactivate enzymes, providing another means to modulate the amount of active enzyme at any moment in time. ◦Transcriptional regulation operates on a much slower time scale. By adjusting the production rate of the proteins, enzyme levels can be adjusted over time. This often occurs across multiple enzymatic steps at any given time. These different regulation mechanisms act on different time scales: ◦Short-term regulation by allosteric regulators and covalent modification is used to ensure that production and consumption of essential molecules such as ATP or monomer units for biosynthesis is balanced. For example, Glycolysis has an intricate set of allosteric regulators to ensure that the supply of ATP matches consumption of ATP elsewhere in the cell. ◦On a longer time scale, enzyme levels are controlled through transcriptional regulation of enzyme production. Changing enzyme levels – or even a complete repression of certain enzymes – can be used to rewire metabolism, 8 to help it adapt to changing external nutrient supplies or when a cell approaches major transition points such as cell division or stress responses. In the interplay of enzyme production and degradation, a quick change in enzyme production will lead to slower, smooth changes of enzyme levels, on a time scale determined by the time scale of enzyme degradation. In growing cells, proteins are diluted, which has a similar effect as degradation, and in microbes, the typical time scale of protein turnover is therefore often the cell cycle time. 3.1.5. Metabolic networks and flux distributions As we have already seen, a metabolic network is a complex web of biochemical reactions. Metabolites produced in one pathway may serve as inputs for others, and the flow of matter across the network is controlled by enzyme activities and regulatory signals. Some pathways are linear, progressing through a fixed sequence of steps (e.g. glycolysis), while others are interconnected or cyclic (e.g. TCA cycle), with multiple inputs and outputs. Alternative pathways may provide redundancy and adaptability to changing conditions. In models, it is often assumed that metabolism operates in a steady state, that is, a state in which concentrations of internal metabolites remain constant over time despite ongoing flux. If the production and consumption rate of a metabolite are very similar, molecules are constantly turned over, maybe at a high rate, but the metabolite concentration remains almost the same. Assuming an exact steady state simplifies analysis and allows for the prediction of metabolic behavior based on network structure and constraints. A steady state resembles a flowing river: while water continuously flows in and out, the overall water level remains stable. Similarly, in metabolism, the input and output fluxes of each metabolite are balanced , maintaining a constant concentration of the metabolite. The same also goes for cofactor pairs: for example, the sums of fluxes in ATP-producing and ATP-consuming reactions must be in balance, which effectively couples fluxes in the entire network. In all this, enzymes act as adjustable parameters or "knobs" in the metabolic network. By modulating enzyme abundance or activity, cells can redistribute fluxes, optimize resource use, and tune metabolite concentrations to achieve homeostasis or adapt to new environments. A "metabolic strategy" is a flux distribution consuming and producing compounds of biological relevance, possibly together with the metabolite and enzyme concentrations that support it. The same type of cells can employ different metabolic strategies depending on environmental conditions (e.g. supply of different nutrients) and physiological goals (e.g. growing versus showing a stress response). Different strategies may involve different choices of nutrients or different alternative pathways that vary in energy efficiency, speed, or biosynthetic output. The environment, especially nutrient availability and oxygen presence, plays a key role in determining which pathways can be active at all, and which pathways are the most economical in terms of enzyme demand. For example, cells may switch between oxidative phosphorylation and fermentation based on oxygen levels. Aerobic respiration requires oxygen as the terminal electron acceptor. In its absence, cells must resort to anaerobic pathways such as fermentation, which are less efficient but still support ATP production. However, the choice between pathways may also depend on the achievable flux per enzyme, resource availability, or regulatory constraints. For instance, while respiration yields more ATP per molecule of glucose (see Figure 3.5), fermentation can allow for a higher ATP production rate at a given amount of enzymes. Therefore, fermentation can produce ATP rapidly, which may be advantageous under high demand, and is used by many cells even in oxygen-rich environments. Some organisms, like yeast, also employ respiro-fermentation: a hybrid strategy where fermentation and respiration occur together in the presence of oxygen, allowing rapid growth despite a lower energy yield compared to pure respiration. The choice between metabolic strategies can be modeled by optimality principles (see Chapters 6 in [1] and 7 in [1]). 3.2. 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 15 Notice that Keq depends only on ∆rG0◦, 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 Keq values of many biochemical reactions is the eQuilibrator tool (equilibrator.weizmann.ac.il) [21,22]. 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 [23]. 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’ [24]. 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) [23,24] (see Box 3.E and the Appendix 3.7). 3.3.2. Enzymes as catalysts of biochemical reactions We mentioned many biochemical reactions to be catalyzed by enzymes. As we saw above, 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 [25,26]). 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 a 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 enzymesubstrate binding or catalytic activity. These allosteric sites, thus, provide a mechanism for regulation of enzymatic reactions (Fig. 3.9). 3.3.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 [27]. Dynamical models of metabolic systems, as with all mathematical models, must therefore balance abstraction of real mechanistic features of a system with achieving a still useful and insight-providing model. At the core of all dynamical metabolic models are rate laws that aim to capture the kinetics of biochemical reactions. Non-enzymatic reactions - the reversible and irreversible mass action rate models All rate models used in metabolic modeling are based on the so-called ‘mass action law’ described in Box 3.E 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 16 Physics box 3.E Mass action law for chemical reactions naA+nbB | {z } substrates k+ −−* )−− k− ncC+ndD | {z } products ξ∗ ∆rG0<0 ∆rG0= 0 ∆rG0>0 ξ(Reaction adv.) Internal energy Thermodynamic interpretation Gibbs free energy of reaction: ∆rG0= ∆rG0◦ +R·T·ln cnc·dnd ana·bnb At equilibrium: ∆rG0◦ =−R·T·ln cnc eq ·dnd eq ana eq ·bnb eq e −∆rG0◦ R·T=cnc eq ·dnd eq ana eq ·bnb eq =Keq (3.5) Kinetic interpretation Backward reaction rate: k−·cnc·dnd Forward reaction rate: k+·ana·bnb At equilibrium: k+·ana eq ·bnb eq =k−·cnc eq ·dnd eq k+ k− =cnc eq ·dnd eq ana eq ·bnb eq =Keq (3.6) 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). the concentrations of the reacting substances’ [24,23] (see Box 3.E). 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 [23]. According to the mass action model, the net 17 (A) (B) substrate([S]) flux: v≈etot·kcat −−−−−−−−−−→ enzyme: etot product([P]) enzyme substrate binding site allosteric site enzyme S enzyme S S, P, I kcat altered kcat Figure 3.9: 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.3.3 and Appendix 3.8 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). rate of any reaction of the form given in Eq. (3.1) is given by; v=k+·ana·bnb−k−·cnc·dnd,(3.7) where small letters denote concentration of the relevant species of the same letter, nidenote the stoichiometric coefficient for species i(as introduced in Equation 3.1), 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.E). 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 a 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.E). When we have the concentrations that lead to the thermodynamic equilibrium of the reaction, i.e. equilibrium concentrations, we will have v= 0 = k+·ana·bnb−k−·cnc·dnd k+ k− =cnc·dnd ana·bnb(3.8) This ratio is known as the reaction’s equilibrium constant Keq 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 ◦0. Thus, we can write; k+ k− =Keq = e−∆rG0◦ R·T(3.9) 18 where ∆rG0◦ is the Gibbs free energy under biological standard conditions, and Rand Tdenote the molar gas constant1and temperature (in Kelvin) respectively (see Box 3.E). It is important to note here that, given Keq is a constant determined by thermodynamics, the parameters k+and k−cannot be chosen independently, i..e k−=Keq/k+. Following on from this last point, it is important to consider a reaction with large Keq, i.e. a reaction for which ∆rG0◦ 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+·ana·bnb(3.10) Enzymatic reactions The mass action rate discussed above forms also the basis of modeling enzymatic reactions. This approach is justified by considering each enzymatic reaction as a series of ‘elementary steps’, each obeying the mass action rate model. To this end, many alternative elementary steps, or ‘enzyme mechanisms’, can be considered to ‘capture’ an enzymatic reaction and subsequently many alternative assumptions can be made to simplify the resulting system of steps. It is also possible to include allosteric regulation or other types of inhibition or activation steps within these elementary steps, allowing generation of a rich variety of enzymatic models and rate laws. Here, we will cover some of the most common of such models, noticing that the construction of these models follows the same general principles of (i) drawing up elementary reactions, (ii) writing down mass action based kinetic rates for the system, and (iii) simplifying the system with assumptions on kinetic parameters (see Appendix 3.7). The reader can consult additional books (e.g. [28]) 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 k1 −−* )−− k2 ES kcat −−→ 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 3.8, 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=s·etot ·kcat KM+s(3.11) where etot represents the total enzyme concentration, kcat is known as the catalytic rate of an enzyme, and KMis known as the Michaelis-Menten coefficient of the enzyme and is equal to (k2+kcat)/k1(we note that depending on the assumptions used, the expression for KMcan vary). Plotting the above rate of formation of product against increasing substrate concentration (see Figure 3.10) shows that the rate is a ‘saturating function’ of substrate, i.e. the rate approaches a threshold point - given by νmax =etot ·kcat 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 νmax, i.e. total enzyme concentration and enzyme’s catalytic rate. This fact underpins the regulation of metabolic flux 1The 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−1 ·mol−1. 19 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 s/KM v/νmax Figure 3.10: Michaelis-Menten rate law – The xand y-axis show the substrate concentration (normalized by KM) and reaction flux (normalized by νmax) respectively. The dashed horizontal line corresponds to νmax, i.e. etot ·kcat. 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 [29]. It considers the following reaction scheme: S+E k1 −−* )−− k2 ES k3 −−* )−− k4 EP k5 −−* )−− k6 P+E. Deriving the rate law for this reaction scheme is slightly more involved, but it follows the same strategy as explained above, of creating elementary steps, treating them as obeying mass action rate, and making additional simplifying assumptions. As shown in Appendix 3.7, we can follow this strategy to derive the reversible rate law as follows: v=etot ·k+ cat KS · s−p·k− cat/KP k+ cat/KS 1 + p KP +s KS (3.12) where KSand KPare composite constants relating to the substrate and product binding to the enzyme, and k+ cat and k− cat are Haldane coefficients (again, composite parameters of other kinetic constants) describing catalytic rate of the enzyme (see Appendix 3.7 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 Keq and reaction Gibbs free energy. Recognizing the relation between the Haldane composite parameters and Keq (see Appendix 3.7) and the flux-force relationship (see below), we can then re-formulate the reversible rate law as: v=etot ·k+ cat ·s/KS 1 + p/KP+s/KS ·1−e∆rG0 R·T(3.13) where ∆rG0is 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, ∆rG00). However, as the reaction Gibbs free energy gets close to zero, the reaction rate will decrease, and as such, there will be a dependence of reaction rate on reaction free energy. 20 Another way of writing equation 3.13 is this one: v=etot ·k+ cat · s/KS·1−e∆rG0 R·T 1 + s/KS·1 + k+ cat k− cat ·e∆rG0 R·T(3.14) where we replace p/KPwith an expression that depends on sand ∆rG0. This alternative expression, developed in the context of modeling microbial metabolism [30,31], can be useful because it shows us that when the reaction is far from equilibrium (∆rG00), the term e∆rG0/(R·T)will approach zero and the above formula can be approximated by the irreversible Michaelis-Menten rate law (Equation 3.11). In this case, we further notice that the Haldane coefficient KSbecomes equivalent to KMintroduced above in the irreversible reaction scheme (Equation 3.11). 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 [32,33,21]. 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=νmax ·s KM+s+s2/KI (3.15) where KIrepresents the saturation coefficient for the binding of the substrate at an allosteric site on the enzyme. Notice that we use such a model in the small multi-stable system example introduced below (section 3.4.3) and discussed in Appendix 3.8. 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 [34]. 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.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.7), and break this into its components of forward reaction rate (or flux) and reverse reaction rate (or flux), which are given by; v+=k+·ana·bnb v−=k−·cnc·dnd(3.16) 21 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 −∆rG0/RT J/v+ Figure 3.11: The ratio of net forward flux (J) to forward reaction rate (v+) as a function of the negative reaction Gibbs free energy and then, we can express the net forward flux (J) as: J=v+−v−=v+·1−v− v+=v+·1−k−·cnc·dnd k+·ana·bnb=v+·1−k− k+ ·Γ(3.17) 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.9)) as a flux-force relationship: J=v+·1−k− k+ ·Γ=v+·1−Γ Keq =v+·1−e∆rG0 R·T(3.18) In this equation, we used the fact that the disequilibrium ratio ρ= Γ/Keq can also be written as a direct function e∆rG0 R·Tof the Gibbs free energy of reaction ∆rG0. For convenience, we also define the thermodynamic driving force θ=−∆rG0/RT and sometimes write ρ= Γ/Keq = e−θ. Eq. 3.18 tells us 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.11). 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 relationship [35,36] and holds also for the enzymatic reversible reaction model described above (see section 3.3.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 Keq 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 Keq – a reaction with a very large Keq 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 22 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 [37]. A recent study by Shen et al [38] showed how important it can be to include product inhibition to create a predictive metabolic model. In the case of lower Keq value – in combination with a consideration of possible product concentration – the modeler should opt for the reversible rate models, which are thermodynamically consistent. The decision about use of enzymatic or non-enzymatic reaction models can be made in a practical manner. If the enzyme associated with the modeled reaction has measured kinetic rates, it would be sensible to opt for a enzymatic model (noting that in vivo enzyme kinetics might differ from those measured in vitro and that many enzyme kinetics studies use parameter derivations assuming an irreversible Michaelis-Menten model). Consequently, it may not be possible to find all the required parameters in the literature, so to model a reaction using reversible rate model. In the absence of measured enzyme parameters, the modeler can use ‘guesstimated’ parameters, based, for example, on the distribution of known enzyme kinetic parameters (see Figure 2.7 in [1] from Chapter 2 in [1]), 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+·1−Γ Keq (3.19) 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 Keq 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+=etot ·k+ cat ·(s/KS)/(1 + s/KS+p/KP)(3.20) And, for the reversible non-enzymatic case: v+=s·k+(3.21) Where kcat,KS, and KPare 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” [36]. By re-arranging the above terms, we can show that correspondence between the two models can be expressed as: etot ·k+ cat ·(1/KS)/(1 + s/KS+p/KP)≈k+(3.22) We can see that in the regime, where sKSand pKP, both models would behave in a linear fashion and their behavior would correspond exactly with the right choice of parameters (i.e. assuming (etot ·k+ cat/KS) = k+). Outside of this regime, correspondence would be dependent on both parameters and concentration of sand p. One interesting case to consider is when total amount of sand pwould be conserved, for example, with cycling reaction schemes. In this case, we can introduce a new parameter cto describe the total pool of the cycled metabolite (e.g. c=s+p) 23 and the correspondence would be expressed as: (etot ·k+ cat/KS)/(1 + s·(KP−KS)/(KS·KP) + c/KP)≈k+(3.23) Thus, in this case when the sum of substrate and product concentrations is conserved, we can have correspondence between the non-enzymatic and enzymatic models when sis small or when KS=KP. 3.4. Dynamics and regulation of metabolism As explained so far in this chapter, cell metabolism involves biochemical reactions involving metabolites (and often catalyzed by enzymes). Thus, understanding metabolism involves studying the dynamics of this system, trying to predict how metabolite levels will go up or down, or settle to a steady state as cell physiology changes in response to external or internal processes (e.g. cells encountering glucose or undergoing division). Obtaining such understanding requires us to develop models of biochemical reaction systems and predict the ‘dynamics’ of those systems. In the previous section, we learned how to model one biochemical reaction. Now we will see how we can readily expand these models to capture multi-reaction systems. The ‘art’ of developing and analyzing dynamical models falls under the branch of mathematics known as calculus and nonlinear dynamics. Many introductory books to these subjects are available, but we find that two particularly useful ones are those by Silvanus Thompson on calculus [39] and by Steven Strogatz on nonlinear dynamics [40]. 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 [25,28,41] 3.4.1. Stoichiometric matrix and 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 mspecies (i.e. the metabolites), and to nreactions, found in the system respectively: Nis a m×nmatrix 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: ds dt=N v(s)(3.24) 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 sdescribes the rate of change in the concentration of swith a given (infinitesimal) change in time. The ODEs can be solved numerically (i.e. simulated) by computer or studied analytically. 24 (A) Thermodynamic steady state naA+nbB−−* )−− ncC+ndD ξ∗ ∆G= 0 ξ(Reaction adv.) Internal energy (B) Dynamic steady states – non-equilibrium thermodynamics M0M1 kin kout ATP ADP Figure 3.12: 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. 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 [39,40], while for a metabolic view of variables and parameters, the article on the Control of Flux, by Kacser and Burns, offers a valuable perspective [42]. 3.4.2. Dynamic steady state As stated above, the ODEs describe the time evolution of all variables sin 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.12). 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.12). 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 31 (A) Allosteric enzyme model 0 20 40 60 80 100 1.0 1.5 2.0 2.5 3.0 S P E1 E2 Consumption Production Stot Flux α αα α At steady state: V1·[S] K1+ [S]+[S]2/K3 | {z } P production =V2·(C−[S]) K2+ (C−[S]) +α·(S0−[S]) | {z } P consumption (B) Multi-site enzyme model S P E12 ES1 ES2 ES1,2 3 complexes ES1,2 n ES1,2,…n … ES1,3 2n 2n1 complexes E2 E1 E12 E2 E1 Pproduction flux (total) E12 E2 E1 Flux Stot α α 0 2 4 6 8 10 0.2 0.4 0.6 0.8 1.0 Figure 3.15: Cartoon representations and brief analysis results of two enzymatic models capable of bistability. (A) Allosteric enzyme model. The first model considers an enzyme that can convert a substrate (S) into a product (P) and that is allosterically regulated by its own substrate. This regulation takes the form of inhibition and is implemented mathematically in the rate of the enzyme - black colored equation. This model results in a nonlinear curve for the relation between rate of production of P at steady state and the total concentration of substrate and product in the system, Stot (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 Stot (red curve on top right panel). We can see that the model is capable of resulting in three intersections, i.e. three steady states of the system. (B) Multi-site enzyme model. The second model considers instead of allostery, an enzyme that binds multiple substrates. This results in several enzyme-substrate complexes depending on the number of binding sites - 3 sites in the model shown. The resulting model can be solved for the steady state values of flux through each enzyme complex against Stot (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. Oscillations - intertwined negative and positive feedbacks Several mathematical models of the reaction catalyzed by the enzyme phosphofructokinase (PFK) in the glycolysis pathway has shown that oscillations are possible to arise from the dynamics of this reaction alone. These models incorporate some of the observed allosteric regulation of PFK both by its substrates and products, resulting in intertwined negative and positive feedbacks. It must be noted that some of these models, and others, use the same basic models that show bistable behavior (as discussed above) and extend them with inand out-fluxes of involved metabolites, to display oscillations. While these theoretical demonstrations of specific enzymatic schemes leading to oscillations have not been explored in detail experimentally, metabolic oscillations are readily observed both in vivo and in vitro, as discussed above. Models, 32 involving some of these proposed synchronization molecules, were also developed and could reproduce experimental findings. 3.7. Derivation of enzymatic rate laws Enzymatic reactions can be modeled using a mechanistic model of enzyme binding and catalysis. The general approach is to develop a ‘cartoon’ model of the physical steps in a reaction. This cartoon model usually takes the form of a series of reactions, involving either binding / unbinding events or chemical conversions. Once a model is developed one can write down ordinary differential equations (ODEs) based on these reactions, and assuming each reaction to be governed by mass action kinetics (see Section 3.3.3). The ODEs can be simplified using certain assumptions, or sometimes just kept as is, before applying a quasi steady-state assumption (which states the enzyme-substrate complexes to be in steady-state). This assumption would allow us to solve the ODE for the enzyme-substrate complex(es) at steady-state. We then enter these solutions into the ODE for the product, so to obtain a reduced system and a specific rate law for product formation. This approach forms the basis of obtaining simplified rate laws, that is, a reduced ODE for the rate of product formation, for enzymatic reactions. 3.7.1. Derivation of the single substrate, irreversible rate law This is the most generic model of an enzymatic reaction that has been developed/studied by Leonor Michaelis (1875 1947) and Maud Leonora Menten (1879 1960), and their contemporaries. It involves the following reaction scheme, where a substrate binds to an enzyme to form a enzyme-substrate complex, gets converted into a product, and then released from the enzyme: S+E k1 −−* )−− k2 ES k3 −−* )−− k4 EP k5 −−* )−− k6 P+E.(3.29) We can simplify this reaction system by assuming that (1) the transition between enzyme complexes ES and EP are instantaneous and are therefore considered as a single entity, e.g. es, and (2) that the release of product and enzyme is irreversible. The scheme now becomes: S+E k1 −−* )−− k2 ES k3 −−→ P+E.(3.30) We can now write a set of ODEs to describe the dynamics of this reaction system - using mass action kinetics. The ODEs are as follows: ds dt=−s·e·k1+es·k2 de dt=−s·e·k1+es·(k2+k3) dc dt=s·e·k1−es·(k2+k3) dp dt=es·k3 where we used the small letter notation to represent the concentration of each species, e.g. “e” for the concentration of the enzyme, E, and “es” for the concentration of the enzyme-substrate complex, ES. At this stage, we can see that if we can formulate “es” as a function of “s”, we can provide a simpler rate model that relates production of the product, P, to the level of the substrate, s. To achieve this we make several additional assumptions. First, we will assume that the total level of the enzyme is conserved, i.e. e+es=C, where Cis a constant (referred to as etot in the main text). This assumption effectively means that total enzyme levels are fixed in the timescale of reaction dynamics. This assumption already allows us to re-define the ODEs and reduce their number to three from four - 33 since, we can now express e, as a function of es. The new ODEs look like this: ds dt=−s·(C−es)·k1+es·k2 des dt=s·(C−es)·k1−es·(k2+k3) dp dt=es·k3 Second, we will assume that the binding/unbinding of substrate to the enzyme happens much faster than release of product from the enzyme-substrate complex. This assumption, together with the additional assumption that enzyme levels are much lower than substrate levels, allows us to consider the enzyme-substrate complex to remain constant throughout the reaction. In other words, we consider the enzyme-substrate complex to be in a ‘quasi steady-state’. This allows us to solve the second ODE from above for steady-state: des dt= 0 = s·(C−es)·k1−es·(k2+k3) es·(k2+k3) = s·(C−es)·k1 es·(k2+k3) = sC ·k1−s·es·k1 es·(k2+k3+s·k1) = s·C·k1 es=s·C·k1 (k2+k3+s·k1) We have now an expression for “es”, which we can simply introduce to the ODE system. We have effectively reduced our ODE system from a three variable system into a two variable one: ds dt=−s·(C−s·C·k1 (k2+k3+s·k1))k1+s·C·k1 (k2+k3+s·k1)·k2 dp dt=s·C·k1 (k2+k3+s·k1)·k3 The second ODE describes the rate of change in product, P, as a function of substrate, S. It is a rate model for this enzymatic reaction, and holds under the assumptions we made in its derivation. It is known as the Michaelis-Menten kinetic rate model and is commonly expressed as: v=s·etot ·kcat KM+s where etot is equal to Cand represents total enzyme concentration, kcat is equal to k3and is known as the maximal catalytic rate of an enzyme, and KMis equal to (k2+k3)/k1and is known as the Michaelis-Menten coefficient of the enzyme. Plotting this rate against increasing substrate concentration would show that the rate is a ‘saturating function’ of s, i.e. the rate approaches a threshold point - given by νmax =etot ·k3as substrate increases. The enzymatic nature of the reaction introduces a limiting factor on the reaction rate! This saddle point is actually a underpinning point for some of the constraint-based methods discussed in this book. 3.7.2. Derivation of a two substrate, irreversible rate law See Problem 3.4 34 3.7.3. Derivation of the single substrate, reversible rate law We now return to the reaction scheme we considered in the above section: S+E k1 −−* )−− k2 ES k3 −−* )−− k4 EP k5 −−* )−− k6 P+E. The corresponding ODE system, written only for the key variables es,ep, and p, is as follows: des dt=e·s·k1+ep·k4−es·(k2+k3) dep dt=e·p·k6+es·k3−ep·(k4+k5) dp dt=ep·k5−e·p·k6 As above, we will now introduce the assumptions of (1) total enzyme being conserved, and (2) the quasi steadystate, but this time for both of the enzyme-substrate and enzyme-product complexes. We will denote total enzyme concentration as C, as before, and use these two assumptions to express esand ep in terms of each other, and the other variables. Let us first proceed with es; des dt= 0 = e·s·k1+ep ·k4−es·(k2+k3) es·(k2+k3) = (C−es−ep)·s·k1+ep ·k4 es·(k2+k3+s·k1) = (C−ep)·s·k1+ep ·k4 es=C·s·k1+ep ·(k4−s·k1) (k2+k3+s·k1) We carry the same derivation for ep; dep dt= 0 = e·p·k6+es·k3−ep ·(k4+k5) ep ·(k4+k5) = (C−es−ep)·p·k6+es·k3 ep ·(k4+k5+p·k6) = (C−es)·p·k6+es·k3 ep =C·p·k6+es·(k3−p·k6) (k4+k5+p·k6) We see that we have a symmetry in the expressions for esand ep, in that the two expressions can be derived from each other by a replacement of variables (k1, k4, k2, s)→(k6, k3, k5, p). Keeping this symmetry in mind, we now 35 attempt to eliminate one of the complexes from the equation for the other: ep ·(k4+k5+p·k6) = C·p·k6+es ·(k3−p·k6) ep ·(k4+k5+p·k6) = C·p·k6+C·s·k1+ep ·(k4−s·k1) (k2+k3+s·k1) ·(k3−p·k6) ep ·(k4+k5+p·k6) = C·p·k6+C·s·k1k3−C·s·k1·p·k6+ep ·(k4−s·k1)·(k3−p·k6) (k2+k3+s·k1) ep ·(k4+k5+p·k6)·(k2+k3+s·k1) = C·p·k6·(k2+k3+s·k1) + C·s·k1k3−C·s·k1·p·k6+ ep ·(k4−s·k1)·(k3−p·k6) ep ·(k4+k5+p·k6)·(k2+k3+s·k1) = C·p·k6k2+C·p·k6k3+C·s·k1k3+ep ·(k4−s·k1)·(k3−p·k6) ep ·((k4+k5+p·k6)·(k2+k3+s·k1)−(k4−s·k1)·(k3−p·k6)) = C·p·k6k2+C·p·k6k3+C·s·k1k3 ep =C·p·k6·(k2+k3) + C·s·k1k3 (k4+k5+p·k6)·(k2+k3+s·k1)−(k4−s·k1)·(k3−p·k6) ep =C·p·(k6k2+k6k3) + C·s·k1k3 (k4k2+k5k2+k5k3+s·k1k3+s·k1k4+s·k1k5+p·k6k2+p·k6k3+p·k6k4) Note that, in the above equation set, we have dropped the dot notation from multiplication of parameters for simplicity of expression. Based on the above argument of symmetry, or by following the same steps for “es”, we can show that we will have a similar expression with different parameters in the numerator: es=C·s·(k1k5+k1k4) + C·p·k6k4 (k4k2+k5k2+k5k3+s·k1k3+s·k1k4+s·k1k5+p·k6k2+p·k6k3+p·k6k4) With these expressions for esand ep at hand, we can now derive an expression for e: e=C−es−ep e=C−C·s·(k1k5+k1k4) + C·p·k6k4 (k4k2+k5k2+k5k3+s·k1k3+s·k1k4+s·k1k5+p·k6k2+p·k6k3+p·k6k4) −C·p·(k6k2+k6k3) + C·s·k1k3 (k4k2+k5k2+k5k3+s·k1k3+s·k1k4+s·k1k5+p·k6k2+p·k6k3+p·k6k4) e=C−C·s·(k1k3+k1k5+k1k4) + p·(k6k2+k6k3+k6k4) (k4k2+k5k2+k5k3+s·k1k3+s·k1k4+s·k1k5+p·k6k2+p·k6k3+p·k6k4) e=C·k3k5+k2k5+k2k4 (k4k2+k5k2+k5k3+s·k1k3+s·k1k4+s·k1k5+p·k6k2+p·k6k3+p·k6k4) We are now ready to substitute all these expressions into the ODE for the product, so to obtain our rate law: dp dt=C·p·(k6k2+k6k3) + C·s·k1k3 (k4k2+k5k2+k5k3+s·k1k3+s·k1k4+s·k1k5+p·k6k2+p·k6k3+p·k6k4)·k5 −C·k3k5+k2k5+k2k4 (k4k2+k5k2+k5k3+s·k1k3+s·k1k4+s·k1k5+p·k6k2+p·k6k3+p·k6k4)·p·k6 dp dt=C·s·k1k3k5−p·k2k4k6 (k4k2+k5k2+k5k3+s·k1k3+s·k1k4+s·k1k5+p·k6k2+p·k6k3+p·k6k4) 36 We can somewhat simplify this expression by defining the following composite rate constants: KS=k3k5+k2k5+k2k4 k1·(k3+k4+k5) KP=k3k5+k2k5+k2k4 k6·(k2+k3+k4) k+ cat =k3k5 k3+k4+k5 k− cat =k2k4 k2+k3+k4 and substituting them into the rate expression from above, to get: dp dt=v=C·k+ cat KS · s−p· k− cat KP k+ cat KS 1 + p KP+s KS This reaction rate is referred to as the Haldane kinetic rate law, named after Jack Burden Sanderson Haldane (5 November 1892 1 December 1964). It can be re-expressed by recognizing the fact that the fraction entering as a multiplier for the product concentration is actually equivalent to the equilibrium constant of the reaction scheme drawn above, at the beginning of this section, when we assume the reaction proceeding in the forward direction, i.e. towards product formation: k− cat KP k+ cat KS =k2k4k6 k1k3k5 = 1/Keq This allows us to re-express the Haldane rate law as: v=C·k+ cat ·s/KS 1 + p KP+s KS ·(1 −p/s Keq ) This re-arranged expression is interesting because we can recognize that the last term is related to the thermodynamic Gibbs free energy of the reaction, allowing us to finally derive: v=C·k+ cat ·s/KS 1 + p/KP+s/KS ·(1 −e∆rG0/RT ) where ∆rG0is the Gibbs free energy of reaction for given substrate and product levels, considering forward direction, and Rand Tstand for the gas constant and temperature respectively. This rate law shows that forward reaction rate will be independent of thermodynamics, when the reaction free energy is highly negative (i.e. thermodynamically highly favored), but the reaction rate will decrease as Gibbs free energy gets close to zero. A second, faster derivation of this rate law is found by noting that the ODEs for des dtand dep dtare linear in e,esand ep, and can therefore be solved with linear matrix algebra. One may write:    s k1−(k2+k3)k4 p k6k3−(k4+k5) 1 1 1      e es ep   =   0 0 C   ,(3.31) 37 where the first two rows of the matrix correspond to des dt= 0 and dep dt= 0, and the last row represents conservation of total enzyme concentration. The equilibrium concentrations of e,esand epare then found by left-multiplying both sides of the equation by the inverse of this matrix. The obtained results are the same as given above. 3.7.4. Derivation of two-substrate, reversible rate law with simultaneous binding The two-substrate case is described by the following reaction scheme: S1+ S2+ E k1 −−* )−− k2 ES1S2 k3 −−* )−− k4 EP1P2 k5 −−* )−− k6 P1+ P2+ E, Where we assume that binding and unbinding of the substrates and products occurs simultaneously. Proceeding as above we let e,es1s2,ep1p2,s1,s2,p1and p2denote the concentrations of E,ES1S2,EP1P2,S1,S2,P1and P2 respectively. The differential equations for es1s2,ep1p2and p1+p2are: des1s2 dt=e·s1·s2·k1+ep1p2·k4−es1s2·(k2+k3) dep1p2 dt=e·p1·p2·k6+es1s2·k3−ep1p2·(k4+k5) d(p1+p2) dt=ep1p2·k5−e·p1·p2·k6. Proceeding as in the single substrate case, we note that the the ODEs for des1s2 dtand dep1p2 dtare linear in e,es1s2 and ep1p2, and that the total enzyme concentration e+es1s2+ep1p2is constant, denoted C.    s1s2k1−(k2+k3)k4 p1p2k6k3−(k4+k5) 1 1 1      e es1s2 ep1p2   =   0 0 C   .(3.32) We therefore see that the results for the two-substrate case are the same as for the single substrate case, with s replaced by s1s2and preplaced by p1p2. This result is dependent on the assumption that binding/unbinding of substrates/products occurs simultaneously. 3.8. A simple model illustrating product activation This model demonstrates that allosteric regulation of an enzymatic reaction by its product can create a bistable system. In this simple example, we consider enzymatic production of a metabolite (labelled ’x’) and its non-enzymatic consumption. It is assumed that the metabolite allosterically regulates the enzyme that produces it. The listing uses the Antimony format [51] which can be easily converted into SBML [52]. An online converter can be found at sysbio.github.io/makesbml. 1// The following model admits three steady-states at: 2// x = 0.325, x = 1.671, and x = 0.873 3// The first reaction step ‘-> x’ uses a rate law that models 4// positive feedback via the product x. The constant 0.2 5// is to ensure that the lower steady-state is non-zero. 6// The statement ‘ext Xo’ indicates that the species Xo is fixed. 7 8ext Xo 9Xo -> x; (vo*x^n)/(1 + x^n) + 0.2 10 x ->; k1*x 11 38 12 k1 = 0.65 13 n = 4; vo = 1 14 x = 0 Listing 3.1: Model illustrating bistability 1#Equivalent model as a differential equation in python: 2def ode (x, t): 3vo = 1 4n=4 5k1 = 0.65 6return [((vo*x**n)/(1 + x**n) + 0.2) - k1*x] Listing 3.2: Equivalent model as a differential equation in python Solutions to problems Problem 3.1 (Elemental composition of bacteria) BioNumbers (BNID 106666) gives mass fractions for E. coli: 47% carbon (47%/12=3.91%); 23% oxygen (23/16=1.43); 14% nitrogen (14/14=1); 6% hydrogen (6/1=6); 10% phosphorus plus sulphur (10/32=0.3). Dividing by atom weights yields the numbers for carbon (47/12=3.91); oxygen (23/16=1.43); nitrogen (14/14=1); hydrogen (6/1=6); phosphorus plus sulphur (10/32=0.3), which can then be approximated by integer numbers. A closer approximation (with a more precise value of the oxygen fraction) is C8H12O3N2. Problem 3.2 (Alien life form) We leave it to you to speculate. What is sure is that most biochemical processes in such an organisms would be much slower, if would be much harder to privilege certain "intended reaction" over concurrent side reactions, and it would not be possible to regulate the branching ratios of fluxes by regulating the catalyst. All this would probably put much stronger constraints on what set of chemical compounds can be used, because the compounds would determine more directly all the reaction rates. In addition, Darwinian evolution could not act on catalysts. However, alternative mechanisms, resembling catalysis, might be used to speed up or regulate reactions, for example, a spatial co-localization of compounds. Problem 3.3 (An irreversible reaction with simultaneous binding) (a) E+S1+S2 k1 −−* )−− k2 ES1S2 k3 −−→ E+P1+P2(3.33) (b) dp dt =k3 s1s2C s1s2+k2+k3 k1 ,(3.34) where p= [P1+P2]and C= [E]+[ES1S2]. Problem 3.4 (A reversible reaction with simultaneous binding) (a) E+S1+S2 k1 −−* )−− k2 ES1S2 k3 −−* )−− k4 E+P1+P2(3.35) (b) dp dt=k3 C(s1s2−k2k4 k1k3p) s1s2+k4 k1p+k2+k3 k1 (3.36) 39 where p= [P1+P2]and C= [E]+[ES1S2]. Problem 3.5 (An irreversible reaction with sequential binding) (a) E+S1 k1 −−* )−− k2 ES1 ES1+S2 k3 −−* )−− k4 ES1S2 k5 −−→ E+P1+P2(3.37) (b) dp dt=k5 s1s2C s1s2+s1k4+k5 k3+s2k5 k3+k2 k1k3(k1+k5),(3.38) where p= [P1+P2]and C= [E]+[ES1]+[ES1S2] Problem 3.6 (An irreversible reaction with random-order binding) (a) E+S1 k1 −−* )−− k2 ES1 ES1+S2 k3 −−* )−− k4 ES1S2 E+S2 k5 −−* )−− k6 ES2 ES2+S1 k7 −−* )−− k8 ES1S2 ES1S2 k9 −−→ E+P1+P2(3.39) (b) dp dt=k9 Cs1s2k1k3(k6+k7s1) + k5k7(k2+k3s2) s1A(s1) + s2B(s2) + s1s2C(s1, s2) + D(3.40) where p= [P1+P2],C= [E]+[ES1]+[ES2]+[ES1S2], and A(s1) = k1k6(k4+k8+k9) + k7(k0+k4)(k2+k1s1) B(s2) = k2k5(k4+k8+k9) + k3(k0+k8)(k6+k5s2) C(s1, s2) = k1k3(k6+k8+k7s1) + k5k7(k2+k4+k3s2) + k3k7k9 D=k2k6(k4+k8+k9) Bibliography [1] The Economic Cell Collective, editor. Economic Principles in Cell Biology. Free online book, 2023. doi: 10.5281/zenodo.8156386. [2] N.M. Grüning, F. Agostini, C. Caldana, J. Hartl, M. Heinemann, M.A. Keller, J.L. Krüsemann, C. Lamperti, C.L. Linster, S.N. Lindner, J. Muenzner, J. Nielsen, Z. Nikoloski, B. Siebers, J.L. Snoep, H. Tenenboim, B. Teusink, S.J. Williams, M.M.C. Wamelink, and M. Ralser. The return of metabolism: biochemistry and physiology of glycolysis. Biol Rev Camb Philos Soc, 2025. doi: 10.1111/brv.70104. [3] Marco Corrao, Hai He, Wolfram Liebermeister, Elad Noor, and Arren Bar-Even. A compact model of escherichia coli core and biosynthetic metabolism. PLoS Comp Biol, 21(10):e1013564, 2025. doi: 10.1371/journal.pcbi. 1013564. [4] G. Gottschalk. Bacterial Metabolism. Springer Series in Microbiology. Springer, New York, 1985. ISBN 9780-387-96153-8. doi: https://doi.org/10.1007/978-1-4612-1072-6. [5] F. C. Neidhardt, J. L. Ingraham, and M. Schaechter. Physiology go the bacterial cell: A molecular approach. Sinauer Associates, 1990. ISBN 0878936084. [6] H. V. Westerhoff, K. J. Hellingwerf, and K. Van Dam. Thermodynamic efficiency of microbial growth is low but optimal for maximal growth rate. Proc Natl Acad Sci U S A, 80(1):305–9, 1983. ISSN 0027-8424 (Print) 1091-6490 (Electronic) 0027-8424 (Linking). doi: 10.1073/pnas.80.1.305. [7] E. Branscomb and M. J. Russell. Turnstiles and bifurcators: the disequilibrium converting engines that put metabolism on the road. Biochim Biophys Acta, 1827(2):62–78, 2013. ISSN 0006-3002 (Print) 0006-3002 (Linking). doi: 10.1016/j.bbabio.2012.10.003. [8] C. Zerfass, M. Asally, and O. S. Soyer. Interrogating metabolism as an electron flow system. Curr Opin Syst Biol, 13:59–67, 2019. ISSN 2452-3100 (Print) 2452-3100 (Linking). doi: 10.1016/j.coisb.2018.10.001. [9] W. H. Schlesinger and E. S. Bernhardt. Biogeochemistry: an analysis of global change. Academic Press, 2013. ISBN 0123858747. [10] A. Jinich, A. Flamholz, H. Ren, S. J. Kim, B. Sanchez-Lengeling, C. A. R. Cotton, E. Noor, A. Aspuru-Guzik, and A. Bar-Even. Quantum chemistry reveals thermodynamic principles of redox biochemistry. PLoS Comput Biol, 14 (10):e1006471, 2018. ISSN 1553-7358 (Electronic) 1553-734X (Linking). doi: 10.1371/journal.pcbi.1006471. [11] U. Barenholz, D. Davidi, E. Reznik, Y. Bar-On, N. Antonovsky, E. Noor, and R. Milo. Design principles of autocatalytic cycles constrain enzyme kinetics and force low substrate saturation at flux branch points. Elife, 6, 2017. ISSN 2050-084X (Electronic) 2050-084X (Linking). doi: 10.7554/eLife.20667. [12] D. C. LaPorte, K. Walsh, and Jr. Koshland, D. E. The branch point effect. ultrasensitivity and subsensitivity to metabolic control. J Biol Chem, 259(22):14068–75, 1984. ISSN 0021-9258 (Print) 0021-9258 (Linking). [13] M. E. Beber, C. Fretter, S. Jain, N. Sonnenschein, M. Muller-Hannemann, and M. T. Hutt. Artefacts in statistical analyses of network motifs: general framework and application to metabolic networks. J R Soc Interface, 9(77): 3426–35, 2012. ISSN 1742-5662 (Electronic) 1742-5662 (Linking). doi: 10.1098/rsif.2012.0490.