1. $v \in \mathcal{F}C$,

2. for all $v' \in \mathcal{F}C$, if $\text{supp}(v') \subseteq \text{supp}(v)$ then either $v' = 0$ or $v' = \alpha v$ for some $\alpha > 0$.

This means that $v$ is an EFM only if there is no non-zero flux vector in the flux cone that uses only a subset of the reactions that are active in $v$. This also means that if any of the flux-carrying reactions in an EFM is deleted, the flux through the resulting reactions must violate the mass-balance constraints and can therefore not occur in steady-state metabolism; the EFMs are thus minimal in the sense that they cannot be reduced further.

One may now wonder how it is possible that we have seen two definitions of EFMs. First, we introduced them as the extreme rays of the flux cone; then, we introduced them as support-minimal metabolic subnetworks. Indeed, the beautiful thing about EFMs is that these two characterizations are equivalent (see the Mathematical Interlude-box for a proof). These two definitions of EFMs are complementary. Understanding EFMs as minimal metabolic subnetworks enables us to interpret an EFM in terms of its biological function; an EFM can be seen as a metabolic strategy that a cell can use to obtain steady-state metabolism, and which it can combine with other strategies to reach its purpose. The extreme ray-interpretation on the other hand, allows us to write an arbitrary flux vector $v \in \mathcal{F}C$ as a combination of EFMs, as is done in (4.7). This also means that we can learn something about all flux vectors $v$ by learning something about all EFMs. For example, if we know that there is no EFM that produces compound $Y$ without using reaction $r$, this immediately implies that there is no flux vector at all that can do this, and that reaction $r$ is thus essential for the production of $Y$. In addition, as we will discuss below, the identification of EFMs as extreme rays enables the use of efficient computational algorithms to enumerate them.

We note that there is currently no limit on the amount of flux that a particular EFM may carry, since $\lambda e^k$ is also an EFM for any $\lambda > 0$ when $e^k$ is an EFM, and consequently the absolute value of any flux vector $v$ in $\mathcal{F}C$ remains unbounded. However, we will see in the next section that this is not necessarily true when additional constraints are introduced.

### 4.2.3. Irreversibility constraints

In principle, all reactions in a metabolic reaction network are able to run in both directions, but in many practical examples certain thermodynamic arguments can be used to justify treating a subset of reactions as *irreversible*, meaning that for all modeling purposes they can only run in one direction. This stems from their activation energy being so large for realistic levels of products and substrates that the rate of the reverse reaction is effectively negligible compared to the rate of the forward. The choice of which reactions to assume irreversible is a decision to be taken by the modeler and may affect the results of the downstream constraint-based analysis. Returning to the simple model for central carbon metabolism provided in Figure 4.2, the reactions importing glucose into the cell, converting it to pyruvate, and the production of product $P_2$ are assumed to be irreversible. A more detailed representation of central carbon metabolism may allow for the first two reactions to be reversible such that ‘gluconeogenesis’ becomes possible. This would result in the introduction of a fourth EFM associated with the conversion, via pyruvate, of the fermentation product to glucose that is then exported out of the cell.

More generally, in our mathematical description of the metabolic network (4.4), we decomposed each reaction into a forward and a reverse reaction, treating each as individual degrees of freedom, and imposing additional irreversibility constraints is thus as simple as removing a reverse reaction from the network. When this is done, the EFMs of the resulting network do not need to be re-calculated, but can be obtained by simply removing all the EFMs that use the removed directions. This nicely shows the complementarity of Definitions 1 and 2 of EFMs: although it is not easy to see that the set of extreme rays of the new flux cone will be a subset of the previous set of extreme rays, it is easy to see that each vector that is a support-minimal flux vector after the removal of some reactions, must have been a support-minimal flux vector before that. The new flux cone is therefore made up of all flux vectors that are expressed as conical combinations of the remaining EFMs, and remains non-bounded.
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 [94, 74]. 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 [76]. Historically important book on enzyme kinetics and enzymatic reaction models.

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

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

## Thermodynamics and physical chemistry

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

## Metabolic system dynamics

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

# Problems

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

## Problem 3.1 An irreversible reaction with simultaneous binding

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

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

## Problem 3.2 A reversible reaction

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

In this box we will prove (following Gagneur and Klamt [159]) that the two characterizations of elementary flux modes that we gave in this chapter are equivalent. To prove this, we should first make precise what we mean when we define EFMs as 'extreme rays' of the flux cone. Let us start generally: a pointed polyhedral cone $\mathcal{P}$ is defined by a full rank matrix $\mathbf{A}$ such that

$$ \mathcal{P} = \{\mathbf{x} \in \mathbb{R}^n : \mathbf{A}\mathbf{x} \geq \mathbf{0}\}, $$

(as illustrated in Figure 4.1). Such a cone is thus defined as all points that satisfy a certain set of linear inequalities, which can be seen as a region bounded off by halfspaces. Note that the definition of the flux cone as given in (4.4) can be written in this form by taking $\mathbf{A} = [\mathbf{N}, -\mathbf{N}, \mathbf{I}_{n \times n}]^T$, where $\mathbf{I}_{n \times n}$ is the identity matrix, so $\mathcal{F}\mathcal{C}$ is indeed a pointed polyhedral cone.

A vector, $\mathbf{r}$, is a ray of $\mathcal{P}$ if $\alpha\mathbf{r} \neq \mathbf{0}$ and for all $\alpha > 0$ we have $\mathbf{r} \in \mathcal{P}$. Each ray $\mathbf{r} \in \mathcal{P}$ has a zero set defined as $Z(\mathbf{r}) = \{i : \mathbf{A}_{i\bullet}\mathbf{r} = \mathbf{0}\}$. Thus, the zero set is the index set of inequalities that are met with equality by the ray $\mathbf{r}$. We call $\mathbf{r}$ an extreme ray when for all $\mathbf{r}' \in \mathcal{P}$ if $Z(\mathbf{r}) \subseteq Z(\mathbf{r}')$ then $\mathbf{r}' = \alpha\mathbf{r}$ for some $\alpha > 0$. In other words, a ray is called extreme if the set of inequalities that it satisfies with equality cannot be increased. With this we are ready for specifying our second definition of EFM, after which we can prove the equivalency of the two definitions.

**Definition 2.** A vector $\mathbf{v}$ is called an elementary flux mode if it is an extreme ray of the flux cone $\mathcal{F}\mathcal{C}$.

**Lemma 1.** In a metabolic network captured by stoichiometric matrix $\mathbf{N}$ in which all reactions are irreversible, the definitions of elementary flux modes as the extreme rays of the flux cone (Def. 2) and as support-minimal steady-state flux vectors (Def. 1) are equivalent.

*Proof.* Let $\mathbf{v}$ be an elementary flux mode according to Definition 1. The first requirement in this definition immediately implies that $\mathbf{v}$ is a ray of the flux cone, where we can define the flux cone as all $\mathbf{x} \in \mathbb{R}^n$ such that $\mathbf{A}\mathbf{x} \geq \mathbf{0}$ with: $\mathbf{A} = [\mathbf{N}, -\mathbf{N}, \mathbf{I}_{n \times n}]^T$. To show that it is also an extreme ray, let us assume that there is another ray $\mathbf{v}'$ such that $Z(\mathbf{v}) \subseteq Z(\mathbf{v}')$. Since all rays of $\mathcal{F}\mathcal{C}$ must satisfy the first $2m$ inequalities, this specifically means that whenever $v_i = 0$ for some $1 < i < n$, also $v'_i = 0$, i.e. $\text{supp}(\mathbf{v}') \subseteq \text{supp}(\mathbf{v})$, but according to Property 2 of Definition 1 we must then have $\mathbf{v}' = \alpha\mathbf{v}$. This implies that $\mathbf{v}$ is indeed an extreme ray, so it is also an EFM according to Definition 2.

For the converse, let $\mathbf{v}$ be an EFM according to Definition 2. Again, this immediately shows that $\mathbf{v} \in \mathcal{F}\mathcal{C}$, so we should now show that it is support-minimal. For that, let $\mathbf{v}' \in \mathcal{F}\mathcal{C}$ such that $\text{supp}(\mathbf{v}') \subseteq \text{supp}(\mathbf{v})$. This means that whenever $v'_i = 0$, also $v_i = 0$. Since $\mathbf{v}$ and $\mathbf{v}'$ both satisfy the first $2m$ inequalities of $\mathbf{A}\mathbf{v} \geq \mathbf{0}$ with equality, and this shows that whenever $\mathbf{v}'$ saturates one of the last $n$ inequalities, then also $\mathbf{v}$ does this, we conclude that $Z(\mathbf{v}) \subseteq Z(\mathbf{v}')$. Using Definition 2 this implies that $\mathbf{v}' = \alpha\mathbf{v}$ for some $\alpha > 0$. This shows that $\mathbf{v}$ is indeed support-minimal, and is thus an EFM according to Definition 1. □

4.2.4. Practical uses of elementary flux modes

EFMs represent the full set of possible metabolic capacities of an organism, which can therefore make EFM analysis a useful tool for biology. To this end, application of EFM analysis to bioengineering has been proposed to guide the genetic manipulation of microorganisms to perform desirable properties such as synthesis of a bio-compound or efficient production of a recombinant protein (e.g. [160, 161]). From a more theoretical point of view, EFMs have also been used in attempts to quantify cellular robustness [162], in particular regarding robustness under genetic perturbations [163]. The relevance of elementary flux mode analysis to cellular robustness stems from the fact that there is rarely a unique conical combination of elementary flux modes for any given flux vector, which implies there are multiple combinations of minimal metabolic pathways to achieve the same desired effect. This redundancy can be interpreted as a measure for the metabolic robustness of an organism, in terms of preserving essential metabolic functionalities under loss of a gene, for example.
(b) Find the rate of product production for this system.

**Problem 3.3 An irreversible reaction with sequential binding**

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

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

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

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

(b) Find the rate of product production for this system. Note that symbolic math tools such as Mathematica, Maple or the SymPy Python library will be helpful for this question (though not essential).
There have also been several ways that EFM analysis has been incorporated into analysis of multi-omics data. For example, on the basis of transcriptomic profiling of microorganisms, metabolic pathways associated with elementary flux modes have been scored according to their probability of carrying flux [164]. The principle here is that, although levels of RNA often serve as a poor proxy for flux through the reaction associated with that particular enzyme's gene, by creating a gene set associated with an entire EFM there might be a better chance of concretely assessing whether the metabolic pathway as a whole is likely to carry flux. The study [164] suggested that the integration of EFM analysis with gene expression data enabled the identification of certain metabolic pathways activated during stress conditions, and that the organization of elementary flux mode utilization in *Saccharomyces cerevisiae* involves a disparate combination of highly specialized and multi-tasking roles. Beyond transcriptomic profiling, isotope tracing experiments in principle provide a much more direct insight into quantifying metabolic flux. To interpret isotope tracing data, an extension of the concept of an EFM was introduced in [165].

### 4.2.5. Computational challenges for elementary flux mode analysis

Enumerating EFMs for large networks can be computationally challenging if not impossible. In principle, EFMs can be found by removing one reaction at a time and solving the resulting mass-balance constraint problem until it is no longer possible to remove a reaction and still obtain a flux vector that satisfies the steady state conditions. However, the one-to-one identification of EFMs with extreme rays of the flux cone that we described before enables the use of algorithms that are specialized in the efficient enumeration of extreme rays of polyhedral cones, such as the double description method [166]. Various tools have been developed for elementary flux mode enumeration based on this algorithm (e.g. EFMTOOL [167] or MetaTool [168]). However, when the size of the metabolic reaction network grows, the number of EFMs scales disproportionately, leading to a combinatorial explosion that effectively makes enumeration impossible for genome-scale networks containing several thousands of reactions [169]. Currently, EFM analysis is therefore restricted to medium-scale reconstructions containing on the order of several hundreds of reactions, and results in the identification of several hundred million EFMs (e.g. enumeration based on the *Escherichia coli* core model results in approximately 272 million EFMs).

Approaches to reduce the complexity of dealing with so many EFMs even for metabolic reaction networks of modest size have also been proposed. These include invoking transcriptional regulatory constraints to eliminate most of the EFMs to be considered in downstream analysis. Imposing additional constraints based on thermodynamic conditions similarly reduces the set of EFMs considerably. A problem with these approaches is evidently that they do still depend on an initial calculation of all EFMs, and so do not solve the problem of enumeration complexity. A rigorous study of the complexity of EFM mode enumeration was performed by Acuña and colleagues [170]. They showed that the decision problem if there exists an EFM containing two specific reactions is NP-complete whilst the complexity of enumerating all EFMs remains open.

Later in this chapter we will explore some alternatives to EFM enumeration that reduce the difficulty of enumeration, cf. Section 4.4.

## 4.3. Additional constraints and flux polyhedra

### 4.3.1. Inhomogeneous linear flux constraints

We have so far been working exclusively with mass-conservation and irreversiblity constraints, which are captured entirely by the stoichiometric matrix where each row is associated with a metabolite concentration at steady state. We also saw that these considerations alone result in a flux cone that is by definition unbounded, meaning that a flux vector in this space is allowed to take on any absolute value (i.e. multiplying a flux vector in the flux cone by an arbitrarily large positive number again returns a flux vector in the flux cone). However, there are physical constraints limiting the magnitude of flux vectors, especially on the values of flux through exchange reactions that may depend on concentrations
Chapter 4

Metabolic flux distributions

Daan de Groot, Wolfram Liebermeister, Maxime Mahout, Stefan Müller, David Ruckerbauer, Felipe Scott,
and David Tourigny

Chapter overview

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

o 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 and 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 (for a review of the various methods, see [96]). 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, conver-
sion, and production of metabolites, as described by the metabolic fluxes. The rate at which a chemical
of extracellular substrates, numbers of transporter molecules in the membrane, or for which we might have direct experimental measurements. Typically, such bounds on flux values are imposed using inequality constraints of the form $v_i^{lb} \le v_i \le v_i^{ub}$ where $v_i^{lb}$ and $v_i^{ub}$ are lower and upper bounds, respectively, for the flux through the $i$th reaction. When reactions have been decomposed into forward and reverse directions, both upper and lower bounds are non-negative where the latter is usually zero.

In the example from Figure 4.2 one may impose an upper bound on the flux value $v_0$, suggesting that there is a maximal rate at which the cell or organism can import glucose from the extracellular environment. In this case the total set of constraints on the flux vector $v$ take the form

$$Nv = 0, \quad v \ge 0, \quad v_0 \le v_0^{ub} \qquad (4.10)$$

where $v_0^{ub}$ is the maximal glucose uptake rate. It is important to note that the new constraint is of a different kind than the mass-balance and irreversibility constraints: the right-hand side of the constraint is non-zero. Constraints that involve a non-zero are called *inhomogeneous constraints*. We can write these constraints in matrix form as

$$Av \ge b \qquad (4.11)$$

with

$$A = \begin{pmatrix} N \\ -N \\ I \\ G \end{pmatrix}, \quad b = \begin{pmatrix} 0 \\ 0 \\ 0 \\ h \end{pmatrix}, \qquad (4.12)$$

where in this particular case

$$G = (-1 \quad 0 \quad 0 \quad 0 \quad 0), \quad h = (-v_0^{ub}). \qquad (4.13)$$

In general, the matrix $G$ will have $P$ rows corresponding to $P$ linear, inhomogeneous constraints of the form

$$\sum_i w_i^p v_i \le h_p, \quad p = 1, \dots, P \qquad (4.14)$$

where each $h_p$ corresponds to a component of the $P$-dimensional vector $h$ and $n$ weights $w_i^p$ ($i = 1, \dots, n$) are supplied for each constraint. Many constraints can be written in this general form, for example, one might imagine modeling a bound on the total flux that a cell can catalyze, which would be captured by setting all weights equal to 1. In this form the set of constraints on $v$ define a *general polyhedron* that is necessarily contained within the flux cone which was based only on the homogeneous constraints: $Nv = 0$ and $v \ge 0$. The additional inhomogeneous constraints serve to sequentially close up the cone such that various (if not all) dimensions become bounded, thus bounding the total magnitude of the flux vector $v$.

In the example described above, bounding the extracellular glucose uptake rate puts an upper bound on the weights of EFM's $e^1, e^3$ illustrated in Figure 4.3, whose support includes the glucose transport reaction. However, the weight of the EFM associated with uptake and oxidation of the fermentation product ($e^2$) can remain unbounded. In short, this restricts the values of weights $\lambda_1, \lambda_2, \lambda_3$ in the representation (4.7) to satisfy

$$\lambda_1, \lambda_2, \lambda_3 \ge 0, \quad \lambda_1 + \lambda_3 \le v_0^{ub}.$$

Recalling that each EFM is associated with an extreme ray of the flux cone coming from mass-balance constraints, for this simple example it is quite straightforward to interpret the geometric consequences of the maximal glucose uptake rate. Any flux vector $v$ in the resulting *flux polytope* is now constrained to lie within the bounded convex plane represented in Figure 4.4, but remains free to take on any value along the remaining extreme ray parameterized by $\lambda_2$. The maximal glucose uptake has therefore closed up the flux cone along two directions while leaving the third untouched, and the geometry of the resulting flux polytope is represented in Figure 4.5. Imposing an upper bound on the uptake rate of
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 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.

For the moment, we therefore make some simplifying assumptions in order to obtain tractable models:

1. **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 macro-molecule synthesis) that happen "outside" our metabolic network are ignored.

2. **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".

3. **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!

4. **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.

5. **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".

6. **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
<figure><img src="image_34.png" /><figcaption>Figure 4.4: Representation of the bounded convex plane within which a flux vector satisfying the mass-balance and maximal glucose uptake constraints must lie. The possible combinations of EFM weights $\lambda_1$ and $\lambda_3$ are contained within or on the line in blue given by the equation $\lambda_1 + \lambda_3 = v_0^{ub}$.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure 4.5: Geometry of the flux polytope (blue) containing any flux vector that satisfies the mass-balance and maximal glucose uptake rate constraints. While bounded in the directions parameterized by $\lambda_1, \lambda_3$, it remains unbounded in the direction parameterized by $\lambda_2$.</figcaption></figure>

the fermentation product, of the form $v_4 \le v_4^{ub}$, will serve to bound this remaining direction of the polytope such that weights of the EFMs are then restricted to the space defined by

$$ \lambda_1, \lambda_2, \lambda_3 \ge 0, \quad \lambda_1 + \lambda_3 \le v_0^{ub}, \quad \lambda_2 \le v_4^{ub}. \qquad (4.15) $$

More general constraints, for larger metabolic reaction networks will be more difficult to interpret and visualize in such simple geometric terms. Quite quickly the combinatorial complexity associated with combinations of multiple constraints and EFMs will become unmanageable. The intuitive treatment of inhomogeneous linear constraints is partially assisted using the concept of elementary flux vectors on which we will add a section in a later version of this book, but both geometrically and biologically these objects are nowhere near as easy to interpret as their EFM counterparts. We shall see that alternative computational methods for exploring flux space therefore become imperative.

As a final remark, we clarify once more that the general form of constraints (4.14) is by no means restricted to sums on the left hand side that involve just a single reaction and can of course include constraints on weighted sums of flux values for different reactions. These weighted sums are often associated with particular biological interpretations: for example, in the metabolic reaction network from Figure 4.2 we might want to restrict our search of flux space to those flux vectors $v$ that produce adenosine triphosphate (ATP) at a rate of at least $v^{ATP}$. Although a more elaborate model would of course include ATP as one of the metabolites, in this example we can use our biological understanding of central carbon metabolism to see that ATP is produced in reactions $v_1$ and $v_3$. A lower bound on ATP production would thus be a lower bound on a combination of $v_1$ and $v_3$ with weights determined by stoichiometry which depends on the organism under investigation. We could write such a constraint as

$$ w_1v_1 + w_3v_3 \ge v^{ATP} \qquad (4.16) $$
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.

7. **Constrain solutions by modeling limited resources** 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 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)$-
with appropriate weights $w_1, w_3$. Such a constraint forms an additional row of the matrix $G$ and we leave it as an exercise for the reader to explore how this affects the geometry of the flux polytope for various values of the weights, minimal ATP production rate and maximal glucose and fermentation product uptake rates. Particular combinations of constraints will be impossible to satisfy simultaneously (i.e. when the minimal rate of ATP production is impossible to achieve under the given bounds on glucose and fermentation product uptake rates), resulting in a flux polytope that is empty. In such cases the set of constraints on $v$ are called *infeasible*.

### 4.3.2. Thermodynamic constraints

In Chapter 3 the basic concepts of chemical thermodynamics were introduced, in particular, the Gibbs free energy of a metabolic reaction was defined in terms of the concentrations of its products and substrates. For a metabolic reaction network with stoichiometric matrix $N$, the vector of Gibbs free energies (one for each reaction in the network) $\Delta_r G'$ can be written in matrix form as

$$ \Delta_r G' = \Delta_r G'^o + RT \cdot N^T \cdot \ln(c) \quad (4.17) $$

where $R$ is the Gas constant, $T$ the temperature and $c$ the vector of metabolite concentrations at steady state. The components of the vector $\Delta_r G'^o$ are the changes in standard Gibbs free energy for each corresponding reaction. Typically, these values are not known precisely for reactions in the network, but can be estimated or approximated from experimental data using methods beyond the scope of this chapter. Similarly, although it is often difficult to accurately measure all metabolite concentrations, in principle the vector $c$ can be obtained experimentally. However, in practice experimental data on $c$ and $\Delta_r G'^o$ are almost never available. Various methods have therefore been developed to combine estimation of $\Delta_r G'^o$ (sometimes with partial measurements of $c$) with advanced computational techniques that allow simultaneous optimization (see next chapter) or sampling (see below) of $v$ and $c$ (or equivalently: $\Delta_r G'$).

The second law of thermodynamics applied to chemical reaction networks can be summarized by saying that every component of the metabolic flux vector $v$ must satisfy the condition

$$ \operatorname{sign}(v_i) = -\operatorname{sign}(\Delta_r G'_i) \quad (4.18) $$

where $v_i$ and $\Delta_r G'_i$ are the $i$th components of $v$ and $\Delta_r G'$, respectively, and $\operatorname{sign}(x)$ denotes the sign of a variable $x$, and $\operatorname{sign}(0) = 0$. It is important to point out that this notation is different to that used previously, where we had assumed all $v_i$ to be non-negative by decomposing each reaction into irreversible forward and backward reactions. Returning to this reversible notation simplifies the inclusion of thermodynamic constraints into constraint-based models and also their interpretation. According to the second law, a reaction can only proceed in a direction where the change in Gibbs free energy is negative. Thus, to be consistent with mass-balance and the second law of thermodynamics, a flux vector $v$ must simultaneously satisfy both (4.1) and (4.18), with $\Delta_r G'$ defined in (4.17). The consequence of these additional constraints on the geometry of the space of metabolic flux distributions is to exclude quadrants incompatible with the signs of $\Delta_r G'$. Equivalently, imposing the second law of thermodynamics on metabolic flux distributions removes regions of the space that are associated with combinations of thermodynamically-infeasible reaction directionalities.

The resulting space of feasible flux vectors is almost always non-convex, which means more advanced computational methods are required to explore it efficiently. The intuitive reason for this is that imposing thermodynamic constraints on top of the mass-balance constraint is usually done in terms of Boolean variables, which breaks the linearity of the problem that we had and exploited so far. Relating this to the EFMs that were discussed previously, it for example becomes clear that any EFM representing an internal cycle—not including any exchange reactions—will never be thermodynamically feasible. Thus, thermodynamic constraints also reduce the set of EFMs that are possible in a metabolic network. Interestingly, it turns out that any thermodynamically-feasible metabolic flux vector can be expressed solely in terms of thermodynamically-feasible EFMs [171], but the converse statement is not true: a linear combination of thermodynamically-feasible elementary flux mode does not necessarily satisfy the thermodynamic constraints. This shows how the workable properties of convex spaces break down as the mathematical models become more complex, in
dimensional stoichiometric matrix N where each of the m rows corresponds to a metabolite and each of the n columns corresponds to a reaction. The entry N<sub>ij</sub> is the coefficient of the i-th metabolite in the j-th chemical reaction. Then, we can gather all n net reaction rates in an n-dimensional flux vector: v = (v<sub>1</sub>,...,v<sub>n</sub>)<sup>T</sup>. This is convenient because the multiplication N v 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: ds/dt = N v. 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 for v,

$$
\mathbf{N} \mathbf{v} = \mathbf{0}. \tag{4.1}
$$

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 *nullspace* of N.

4.2.2 Irreversibility constraints

In principle, all reactions in a metabolic reaction network are able to run in both directions, but in many
practical examples certain thermodynamic arguments can be used to justify treating a subset of reactions
as irreversible, meaning that in a given model they can run in only one direction. The choice of which
reactions to assume irreversible depends on the experimental conditions and affects the results of the
downstream constraint-based analysis.

Due to microscopic reversibility, the net reaction rate $v_i$ (of reaction $i$) is the difference of the forward and backward reaction rates, that is, $v_i = v_i^- - v_i^-$ (with both $v_i^- > 0$ and $v_i^+ > 0$ if all reactants are present), and $v_i$ can be either positive, zero, or negative. As stated above, thermodynamics may determine the direction of certain reactions, that is, the sign of the net reaction rate. In this sense, if a reaction proceeds in the forward reaction, one adds the nonnegativity constraint $v_i \ge 0$. (Conversely, if a reaction proceeds in the backward reaction, one redefines the reaction by exchanging forward and backward and again adds $v_i \ge 0$.) For a compact notation, let $\mathcal{R}^- \subseteq \{1, \dots, n\}$ be the index set of the irreversible reactions (and $\mathcal{R}^+ \subseteq \{1, \dots, n\}$ be the reversible reactions). We require $v^- := v_{\mathcal{R}^-} \ge 0$, that is, $v_i \ge 0$ if $i \in \mathcal{R}^-$.

4.2.3 The flux cone

Mass balance and irreversibility constraints together define the flux cone

$$
C = \{ \mathbf{v} \mid \mathbf{N}\mathbf{v} = \mathbf{0}, \mathbf{v} \xrightarrow{\mathbf{r}} \mathbf{0} \}. \qquad (4.2)
$$

Elements of the flux cone are called *flux modes*. The flux cone $\mathcal{C}$ is called an *s-cone* (subspace cone) in Müller and Regensburger (2016) [97], since it arises from a linear subspace and nonnegativity constraints.

To provide a concrete example, we consider the simple representation of central carbon metabolism pre-
sented in Figure 4.1. In this example there are four external metabolites, G<sub>ex</sub>, O, P<sub>1</sub>, P<sub>2</sub> and two internal
metabolites: G and P. In our model we only require mass-balance for internal metabolites, such that the
steady-state constraint can be written as

$$
\mathbf{N}\mathbf{v} = \begin{pmatrix} 1 & -1 & 0 & 0 \\ 0 & 2 & -1 & -1 \end{pmatrix} \begin{pmatrix} v_1 \\ v_2 \\ v_3 \\ v_4 \end{pmatrix} = \mathbf{0}, \quad (4.3)
$$

where each column thus corresponds to one of the four (reversible or irreversible) reactions, and where
the rows correspond to *G* and *P* respectively. The entry 1 in the first row of the first column thus corre-
sponds to the import of one glucose molecule *G*. The mass-balance equations

$$
v_1 - v_2 = 0, \quad 2v_2 - v_3 - v_4 = 0, \tag{4.4a}
$$
this case by accounting for thermodynamics.

## 4.4. Alternative methods for flux space exploration

As we described above, exploration of all possible flux distributions using EFMs can become very complex for larger models. A genome-scale model, which comprises all metabolic reactions that an organism can catalyze, typically contains thousands of reactions, which prohibits the enumeration of EFMs. At the moment, it is unclear whether, even if we would have an enormously fast computer that could compute all EFMs, the number of EFMs would not be so large that we cannot store the EFMs anywhere, nor analyze it in any meaningful way. Here we discuss several alternatives for exploring the metabolic capabilities of a cell that try to avoid the combinatorial complexity that hinders EFM analysis.

### 4.4.1. Elementary conversion modes

If we are interested in the metabolic capabilities of an organism, is it always necessary to know all possible flux vectors? For example, what if we want to lab-culture an organism of which we have a reconstructed metabolic network, but no idea what nutrients it needs to grow. Then we only need to know from what combinations of nutrients it can make all its cell components. Or, what if we want to model the possible cross-feeding interactions between several microbial species? Then we are mostly interested in what each of them can consume and produce, and not really in how they do that. Elementary conversion modes (ECMs), introduced in 2005 by Urbanczik and Wagner [172], capture all possible overall conversions from nutrients to products that an organism can catalyze, while ignoring which individual reactions are used for this.

ECMs focus on the net results of metabolism, i.e. on the uptake and production of compounds external to the metabolic network, such as sugars, nitrogen sources, fermentation products but also 'biomass'. To get information about these compounds we need to extend our metabolic network by including the external compounds as rows in the stoichiometry matrix; this is in general easy to do since we already had exchange reactions (reactions where an external compound was imported or exported) so we only have to find the stoichiometric coefficient in which the external compound was involved in these reactions. Let us denote the original stoichiometry matrix by $N_{\text{int}}$ and the submatrix that we add by $N_{\text{ext}}$; together they form $N_{\text{tot}}$. We can then define the conversion cone:

$$C = \{ \dot{c} = N_{\text{ext}}v \mid N_{\text{int}}v = 0, v \ge 0 \}. \qquad (4.19)$$

If we look carefully at this definition we can see that the flux vectors $v$ need to satisfy exactly the same constraints as in the flux cone (Eq. (4.4)). The only difference between flux and conversion cones is that we are either interested in the fluxes themselves, or rather in the conversions that they induce: $\dot{c} = N_{\text{ext}}v$.

**Definition 3.** The set of ECMs is the minimal set of conversions {$ecm^1, \dots, ecm^K$} (where $ecm_k^i$ is the amount of metabolite $k$ produced in the $i$th elementary conversion mode), such that

1. all conversions $\dot{c} \in C$ can be written as a positive sum of these elementary conversion modes: $\dot{c} = \sum_i \lambda_i ecm^i$, with $\lambda_i \ge 0$,

2. without the production of any metabolite being canceled in that sum, i.e. for all metabolites $k$: we either have for all $\lambda_i > 0$ that $ecm_k^i \ge 0$ or for all $\lambda_i > 0$ that $ecm_k^i \le 0$.

We will explain both parts of this definition below, but let us first remark that the definition is in fact perfectly analogous to the definition of EFMs: EFMs the *elementary vectors* (or precisely: conformally non-decomposable vectors) of the flux cone, and ECMs of the conversion cone. The reason that the definition of ECMs has an additional requirement (2.) is just that the analogous requirement was automatically satisfied for EFMs because we assumed all reactions to be irreversible.

In Figure 4.6A we show a small metabolic network with external metabolites A, B and BM, and internal metabolites
Σ Math box 4.A Generators of a polyhedral cone

For every polyhedral cone, and hence for every subspace cone such as the flux cone $C$, there exists a finite, minimal set of *generators* (minimal in the sense that no proper subset forms a generating set). In particular, for the flux cone $C$, there exists a finite set {$f^{(1)}, \dots, f^{(\ell)}$} of $n$-dimensional vectors such that

$$C = \left\{ \mathbf{v} \mid \mathbf{v} = \sum_{k=1}^{\ell} \lambda_k f^{(k)} \text{ with } \lambda_k \ge 0 \right\},$$

that is, any flux vector $\mathbf{v}$ in the flux cone $C$ can be expressed as a conical (non-negative) linear combination of generators {$f^{(1)}, \dots, f^{(\ell)}$}.

*Remark.* The generators $f^{(k)}$ can be multiplied with scalars, that is, any $\lambda f^{(k)}$ with $\lambda > 0$ could replace $f^{(k)}$ in the set of generators.

*Remark.* For a general polyhedral cone, there is no *unique* minimal generating set. However, there is a unique minimal set of *conformal* generators; see [97, Section 3.4]. For subspace cones (such as the flux cone), these are the support-minimal vectors (EFMs); see Theorem 1 below.

and the non-negativity conditions

$$v_1, v_2, v_3 \ge 0, \qquad (4.4b)$$

induced by the irreversible reactions 1, 2, 3, define the flux cone $C$ as the space of all flux vectors $\mathbf{v}$ that satisfy all of these constraints simultaneously.

## 4.3 Elementary flux modes

Equation (4.2) gives a mathematical definition of the flux cone (via equations and inequalities). Here, we will provide an equivalent characterization of this space (via generators, see Math box 4.A). Note that definition (4.2) makes it easy to check if a given steady-state flux distribution $\mathbf{v}$ lies in $C$. However, it is not clear how to generate the flux cone. As a set of generators, we will introduce “minimal” flux distributions, called *elementary flux modes* (EFMs), that can be combined to generate all possible flux distributions in $C$. These EFMs generate the flux cone, similar to how basis vectors generate a linear subspace (but with non-negative coefficients).

In order to define EFMs formally, we introduce the *support* of a vector $\mathbf{v}$ as the index set $\text{supp}(\mathbf{v}) = \{i \mid v_i \neq 0\}$, that is, the support of a flux vector is the set of reactions that have a nonzero rate.

**Definition 1.** A nonzero vector $\mathbf{v} \in C$ is an EFM if it is *support-minimal*, that is, if $\text{supp}(\mathbf{v}') \subseteq \text{supp}(\mathbf{v})$ for any nonzero vector $\mathbf{v}' \in C$ implies $\text{supp}(\mathbf{v}') = \text{supp}(\mathbf{v})$.

*Remark.* If $\mathbf{v}$ is an EFM and $\text{supp}(\mathbf{v}') = \text{supp}(\mathbf{v})$, then further $\mathbf{v}' = \lambda \mathbf{v}$ for some scalar $\lambda$.

Definition 1 states that $\mathbf{v}$ is an EFM if there is no nonzero flux vector in the flux cone that uses only a strict subset of the reactions that are active in $\mathbf{v}$. This also means that if any of the flux-carrying reactions in an EFM is deleted, the flux through the remaining reactions must violate the mass-balance constraints and can therefore not occur in steady-state metabolism; the EFMs are thus minimal in the sense that they cannot be reduced further.

To illustrate the concept of EFMs, we return to the simple representation of central carbon metabolism presented in Figure 4.1 with the stoichiometric matrix

$$\mathbf{N} = \begin{pmatrix} 1 & -1 & 0 & 0 \\ 0 & 2 & -1 & -1 \end{pmatrix} \qquad (4.5)$$
<figure><img src="image_36.png" /><figcaption>Figure 4.6: A) Small toy network with three ECMs shown in blue, yellow and red. Note that the red mode can be decomposed as a positive combination of the blue and yellow elementary conversion modes, but that would cancel the production of *B* so this is not allowed. B) The conversion cone is shown in gray, and the blue and yellow arrow correspond to the blue and yellow ECMs are the extreme rays. The red ECM needs to be added because it is on the intersection with the $\dot{B}=0$-plane.</figcaption></figure>

C, D and E. We can find 9 EFMs in this network: one that goes from A to B, four that produce BM starting from A and four that produce BM from B. We get four EFMs to go from A to BM because there are two ways of going from C to D and again two for converting D into E. This makes clear that having a number of modules of alternative reactions can quickly give rise to large numbers of EFMs, even though the overall conversion from nutrients to products remains the same. In contrast, we will explain that we only get three ECMs.

In Figure 4.6B we see the conversion cone in gray. Note that this cone does not live in flux space, but rather in the space of external metabolite changes, or conversions. We recognize that the cone can be spanned by two extreme rays, which correspond to converting *A* into *B* (blue) and to using 2*B* to produce *BM* (yellow), so these rays correspond to elementary conversion modes following the first part of Definition 3. Now why do we have a third ECM, when the blue and yellow one already span the whole conversion cone? Indeed, the third vector in Figure 4.6B can be obtained by summing the yellow vector and two times the blue vector: $2(-1, 1, 0) + (0, -2, 1) = (-2, 0, 1)$. However, note that the production of metabolite *B* would cancel in this sum, which is not allowed according to the second part of Definition 3. The reason that this second part of the definition is important, is that the elementary conversion modes are intended to capture all metabolic capabilities of an organism, so taking only the first two modes would not be enough: we also want to account for the possibility of making *BM* from *A* even if we decide that the elementary conversion mode from *B* to *BM* is not possible in the current environment, for example because *B* is not present as a nutrient in the medium.

Because many EFMs result in the same overall conversion, the exploration of metabolic capabilities can now be done in larger networks, at the cost of ignoring information about which reactions are used [173]. This way of thinking can be pushed even further: what if one is not interested in the conversions between all nutrients and products, but only between a subset of these? In that case, one would want to compute the ECMs only between the external metabolites of the most interest. This can be done with a small trick. Say that we are not interested in the production of external metabolite *X*. Before we start the enumeration algorithm we add a virtual reaction to the network that consumes and produces *X* from nothing, i.e. we add $X \rightleftharpoons \emptyset$, and then we change *X* from an external metabolite to an internal metabolite. Consequently, it now has to satisfy the mass-balance constraint (which can always be done trivially using the added virtual reaction), and will thus never show up in the computed elementary conversions. In this way it was possible to compute all ECMs between glucose, oxygen and biomass for a real genome-scale network of *E. coli*.
<figure><img src="image_48.png" /><figcaption>Figure 4.1: A simple representation of central carbon metabolism as a metabolic network. (A) Extracellular glucose, $G_{ex}$, is imported into the cell via reaction 1, and intracellular glucose, $G$, is converted to pyruvate, $P$, via reaction 2, having stoichiometric coefficients of two pyruvate molecules for one glucose molecule. Pyruvate is then either converted to a fermentation product, $P_1$, via reaction 4 or, in the presence of oxygen, $O$, converted to an oxidative phosphorylation (OXPHOS) terminal product $P_2$ via reaction 3. The fermentation product $P_1$ can also be converted back to pyruvate via the backward reaction of 4. (B) EFMs $f^{(1)}, f^{(2)}, f^{(3)}$. From our understanding of central carbon metabolism, $f^{(1)}$ represents glycolytic fermentation, $f^{(2)}$ the oxidative metabolism of glucose, and $f^{(3)}$ the oxidative metabolism of the fermentation product.</figcaption></figure>

and the flux vector $v = (v_1, v_2, v_3, v_4)^{\top}$, where $v_1, v_2, v_3 \ge 0$. As it turns out, the set of EFMs is given by

$$f^{(1)} = \begin{pmatrix} 1 \\ 1 \\ 0 \\ 2 \end{pmatrix}, \quad f^{(2)} = \begin{pmatrix} 1 \\ 1 \\ 2 \\ 0 \end{pmatrix}, \quad f^{(3)} = \begin{pmatrix} 0 \\ 0 \\ 1 \\ -1 \end{pmatrix}, \qquad (4.6)$$

and these are depicted in Figure 4.1 (B). From our understanding of central carbon metabolism, we see that these three EFMs represent the “minimal” metabolic pathways of ($f^{(1)}$) glycolytic fermentation, ($f^{(2)}$) oxidative metabolism of glucose, and ($f^{(3)}$) oxidative metabolism of the fermentation product.

In Math box 4.A, we have characterized a polyhedral cone (the flux cone $C$) in terms of its generators (the EFMs). In our toy carbon metabolism network, this means that any flux vector $v$ can be viewed as a conical 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*.

In Math box 4.A, we have also mentioned that EFMs need not form the *unique* minimal set of generators, but they form the *unique* minimal set of *conformal* generators. We first motivate conformality by thermodynamic arguments and then provide a formal definition. For every reaction, Gibbs free energy determines its direction, and hence for every flux, it determines its sign $(-, 0, +)$. Now, since any flux vector is the conical combination of EFMs, the signs of the flux vector determine the signs of the EFMs. In particular, if a certain flux component is zero, then this flux component is zero in all EFMs. (Zero flux cannot arise from a cancellation of positive and negative fluxes.) If a certain flux component is nonzero, then this flux component has the same sign or is zero in all EFMs. (Zero flux can arise from a zero enzyme concentration and hence is thermodynamically sound.)

The above argument can be formalized as follows: For a vector $v \in \mathbb{R}^n$, we obtain the sign vector $\text{sign}(v) \in \{-, 0, +\}^n$ by applying the sign function componentwise, that is, $\text{sign}(v)_i = \text{sign}(v_i)$. In order to capture “conformal signs”, we define the partial order $0 < -$ and $0 < +$ on $\{-, 0, +\}$, which implies the inequalities $0 \le 0$ (zero flux conforms to zero flux), $+, 0 \le +$ (positive or zero flux conforms to positive flux), and $-, 0 \le -$ (negative or zero flux conforms to negative flux). The partial order on $\{-, 0, +\}$ induces a partial order on $\{-, 0, +\}^n$: for two sign vectors $\sigma, \tau \in \{-, 0, +\}^n$, we write $\sigma \le \tau$ if the inequality holds componentwise, and we say that $\sigma$ conforms to $\tau$. If $\sigma \le \tau$ (and $\tau_i$ is given), then $\sigma_i = \tau_i$ or $\sigma_i = 0$. To summarize, if $\sigma$ conforms to $\tau$, then it has the same entries or some more zeros.
### 4.4.2. Flux sampling

In addition to the computational complexity of EFM enumeration for large metabolic networks, these objects are not necessarily related to experimentally-derived flux measurements. This is because when a vector of experimentally-measured flux values $v$ would be decomposed into EFMs, this generally does not give a unique solutions because it can be done in many ways. Flux sampling methods can be employed to solve both the computational and interpretability problems simultaneously, exploring the set of flux vectors (i.e. directly measurable in principle) by computationally sampling from the flux space. The goal of flux sampling in general terms is to produce a sequence of flux vectors that satisfy the steady state constraints until enough samples have been generated to provide an approximate representation of the entire flux space. The flux polyhedra defined by mass-balance and additional inhomogenous linear constraints are convex, and therefore uniform sampling of these flux spaces can be achieved using variants of an algorithm developed for convex analysis called the coordinate hit-and-run (CHR) algorithm [174]. Briefly, the most basic implementation of the CHR algorithm generates a Markov chain of flux vectors by starting in a random position within the flux polytope, picking a direction at random (uniform), and moving a random distance (uniform) in that direction from the current point. The resulting point is returned as a flux vector instance and the process repeats from there. It has been proven that the CHR algorithm converges to a stationary distribution of the Markov chain that is a uniform distribution in the flux space. Alternatives to uniform sampling (i.e. alternative distributions across the flux polytope) can also be achieved using variants of the CHR algorithm.

As highlighted previously in Section 4.3.2, mass-balance and inhomogeneous linear constraints alone often do not contain enough information to sufficiently reduce the space of biologically-feasible flux vectors. For example, thermodynamic constraints on flux vectors are important for ruling out a large proportion of the sampled flux vectors as infeasible, but this may disproportionately dominate the resulting sampling distributions. Unfortunately, for mathematical reasons too deep to go into here, simply removing these infeasible flux distributions post-sampling will not result in a uniform distribution over the thermodynamically-feasible portion of flux space. In fact, this relevant subset of flux space cannot be defined explicitly, and is usually neither convex nor connected meaning that no Markov chain methods exist for sampling. As an alternative, a recent method [175] has been developed to combine thermodynamic constraints, physiological observations and estimated thermodynamic parameters, with mass-balance and inhomogeneous linear constraints to provide a probabilistic thermodynamic analysis of metabolic reaction networks. Advances such as these will almost certainly aid a more complete characterization of flux space as data and methods become available.

### 4.4.3. Minimal cut sets

A minimal cut set (MCS) is a set of reactions that, when disabled, disables a set of modes, which in turn can represent a biological function, such as the secretion of a side product. This enables the prediction of gene deletion targets, given that the genes coding for the involved reactions are known. A cut set is minimal if the removal of one or more reactions from the set leads to at least one of the targeted modes not being disabled.

In order to avoid also disabling desired functionalities, such as product secretion and growth, the concept of constrained minimal cut sets (cMCSs) has been developed. cMCSs enable targeting a set of modes while at the same time making sure that some elements of another set of modes will remain active.

**Motivation for (constrained) Minimal Cut Sets** The concept of MCSs was introduced by Klamt and Gilles in 2004 [176] and subsequently generalized and improved [177, 178, 179]. As briefly outlined above, the idea is to define a set of EFMs which should be disabled, for example because they generate an unwanted side product or because they don't generate the product of interest with a sufficiently high yield. Since EFMs are minimal, removing a single reaction will disable it. A cut set is a set of reactions of which at least one is active in each of the EFMs in the targeted group. Thus, disabling the reactions contained in the cut set will disable all of the targeted EFMs, and each cut set therefore represents the prediction of a set of gene deletions. Since it would be pointless to remove reactions which only target
Now, we can refine the characterization of a flux cone in terms of generators, as given in Math Box 4.A. Indeed, we have the following conformal sum theorem, see [97, Theorem 3].

**Theorem 1.** Let $\mathcal{C}$ be the flux cone and $\{f^{(1)}, \dots, f^{(\ell)}\}$ be the set of EFMs. Then,

$$ \mathcal{C} = \left\{ \mathbf{v} \mid \mathbf{v} = \sum_{k=1}^{\ell} \lambda_k f^{(k)} \text{ with } \lambda_k \ge 0 \text{ and } \operatorname{sign}(f^{(k)}) \le \operatorname{sign}(\mathbf{v}) \right\}. $$

That is, any flux vector $\mathbf{v}$ in the flux cone $\mathcal{C}$ can be expressed as a conformal sum of EFMs $\{f^{(1)}, \dots, f^{(\ell)}\}$.

Again, we illustrate the theoretical concepts in the simple representation of central carbon metabolism. The flux distribution $\mathbf{v} = (1, 1, 1, 1)^T$ lies in the flux cone, cf. (4.4), and hence can be written as a conical linear combination of EFMs (in a non-unique way):

$$ \begin{align} \mathbf{v} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} &= f^{(1)} + f^{(3)} = \begin{pmatrix} 1 \\ 1 \\ 0 \\ 2 \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ 1 \\ -1 \end{pmatrix} \tag{4.7} \\ &= \frac{1}{2} f^{(1)} + \frac{1}{2} f^{(2)} = \begin{pmatrix} \frac{1}{2} \\ \frac{1}{2} \\ 0 \\ 1 \end{pmatrix} + \begin{pmatrix} \frac{1}{2} \\ \frac{1}{2} \\ 1 \\ 0 \end{pmatrix}. \notag \end{align} $$

Note that the first sum is not conformal: The fourth component of $\mathbf{v}$ is positive, whereas the corresponding component of $f^{(3)}$ is negative. That is, the contributing EFMs have different signs in the net reaction rates of the fourth reaction, which leads to cancellation and is not meaningful thermodynamically. (Gibbs free energy determines reaction directions, see Section 4.4.3.) Still, the second sum is conformal: no cancellation occurs, and the decomposition is thermodynamically meaningful. Theorem 1 states that a decomposition as a conformal sum is always possible.

On the one hand, we introduced EFMs as the support-minimal vectors of the flux cone, corresponding to minimal metabolic subnetworks. On the other hand, EFMs form the (unique minimal) set of (conformal) generators of the flux cone. Indeed, the beautiful thing about EFMs is that they have several equivalent (but complementary) definitions, see Math box 4.B for examples and proofs.

Viewing EFMs as minimal metabolic subnetworks enables us to interpret an EFM in terms of its biological function; an EFM can be seen as a metabolic strategy that a cell can use to obtain steady-state metabolism, and which it can combine with other strategies to reach its purpose. The interpretation as conformal generators allows us to write an arbitrary flux vector $\mathbf{v} \in \mathcal{C}$ as a combination of EFMs in a thermodynamically meaningful way, see Theorem 1. This also means that we can learn something about all flux vectors $\mathbf{v}$ by learning something about all EFMs. For example, if we know that there is no EFM that produces compound $Y$ without using reaction $r$, this immediately implies that there is no flux vector at all that can do this, and that reaction $r$ is thus essential for the production of $Y$.

Finally, after reaction splitting, as described in Section 4.3.2, the flux cone is contained in the non-negative orthant and hence is pointed. Then, EFMs agree with the extreme vectors and can be computed via algorithms based on the double-description method, as discussed in Section 4.4.4.

So far, we did not consider a limit on the amount of flux that a particular EFM may carry, since $\lambda f^{(k)}$ is an EFM for any $\lambda > 0$ and any EFM $f^{(k)}$, and consequently the absolute value of any flux vector $\mathbf{v}$ in $\mathcal{C}$ is unbounded. In Section 4.4, we will see that this is not necessarily true when additional constraints are introduced.

### 4.3.1 Practical relevance of EFMs

EFMs represent the full set of possible metabolic capacities of an organism, which can therefore make EFM analysis a useful tool for biology. To this end, application of EFM analysis to bioengineering has been
EFMs that were already targeted by other reactions, cut sets are required to be minimal. This means that removing a single reaction from the cut set would lead to one or more of the targeted EFMs to survive the intervention and also that adding a single reaction to the cut set would have no additional effect on the set of target EFMs.

The pitfall when using MCSs is that while they guarantee the elimination of the targeted EFMs, all other EFMs may be affected as well. This means that modes with desired phenotypes, such as high growth and/or high product yield, may become impossible. Therefore, cMCSs were developed [180]. In this extension of the concept of MCSs it is now possible to additionally define a set of EFMs which are desired, i.e. which can not be disabled by the cMCSs. This is usually implemented by the requirement that at least a specified minimum number of EFMs of the desired set need to remain active. Summarizing, cMCSs are sets of reactions which guarantee that (i) the full set of target EFMs is disabled and (ii) a certain minimum of desired EFMs has to remain unaffected. The drawback, with both MCSs and cMCSs, is that the target (and desired) EFMs need to be defined. This is generally achieved by defining cut-offs in terms of product yield and growth, which is, however, ultimately arbitrary.

**Calculation of (constrained) Minimal Cut Sets** Since minimal cut sets in a metabolic network are EFMs in a dual network [181], methods used for calculating EFMs can be used to calculate MCSs. Among other approaches [182] one based on binary integer programming has been developed [183, 184]. While it requires that the EFMs are calculated before it can be applied, the advantage is that the algorithm is very intuitive. After having calculated the modes, each is represented as a binary vector which is zero for reactions with zero flux and one otherwise. The EFMs are then divided into either targeted or desired. A binary vector, corresponding to the cMCSs being calculated is introduced. It will have a one if the corresponding reaction remains active and zero if the reaction is disabled. The first requirement is that cMCS needs to disable all target modes and thus the vector must have zero elements such that each target EFM must have at least one corresponding non-zero element. The second requirement is that at least a defined minimum of desired modes must remain active. This is achieved by introducing a second binary vector. This vector has an element for each EFM and is calculated so that it has a zero when the mode is disabled by the cMCS and one otherwise. By adding the constraint that the number of ones in this vector must at least equal the previously defined minimum, the second requirement is met. Maximizing the vector corresponding to the cMCS yields the first solution. The next solution can be found by adding constraints to make sure that the current one is excluded.

## 4.5. Conclusions of our flux cone explorations

In this chapter we studied how the individual reactions that an organism can catalyze together give rise to the overall conversion of nutrients into cell components and secretion products. For that, we studied the cell's metabolism under a number of simplifying assumptions, most notably, we model metabolism in steady-state. Given this steady-state constraint, we explained how all feasible flux distributions form a space of a specific type: a pointed polyhedral cone. By exploring this 'flux cone' we can chart the metabolic capabilities of an organism.

We have seen that an exhaustive charting of these metabolic capabilities is the computation of all *elementary flux modes*: minimal subnetworks that can individually give rise to steady-state flux distributions, and that may be interpreted as minimal metabolic strategies. An especially important use of EFM analysis can be found in the prediction of the effect of gene knockouts: when all EFMs that produce compound *Y* use reaction *r*, then the organism cannot make this compound when the gene is knocked out that codes for the enzyme that catalyzes *r*. And conversely, sometimes gene knockouts can be found such that the cell cannot grow anymore without producing a certain compound of interest. Clearly, these analyses can be very useful for the design of organisms in bio-industry.

On the other hand, we also saw that for large models the computation of all EFMs becomes impossible. There are simply too many of these minimal subnetworks. We presented several alternatives. One could use *elementary conversion modes* if one still desires an exhaustive list of the metabolic capabilities of the cell. The ECMs are easier to enumerate because one can choose to focus only on all possible conversions between (a subset of) the nutrients and products, instead of requiring all information about which reactions are used to get these conversions. For the design of gene
<figure><img src="image_49.png" /></figure>

Math box 4.B Equivalent Definitions of Elementary Flux Modes (EFMs)

In the main text, we have introduced EFMs as the support minimal vectors of the flux cone, see Definition 1. In fact, EFMs can be defined as the support-minimal, support-wise non-decomposable, sign-minimal, sign-wise non-decomposable, and conformally non-decomposable vectors of the flux cone; cf. [97]. Here, we consider the latter definition for three reasons: (i) it matches Theorem 1 on the decomposition of flux distributions into conformal sums of EFMs, (ii) it also applies to general polyhedral cones (not just s-cones such as the flux cone) and even to polyhedra and polytopes, and (iii) it establishes a link to the case when the flux cone is contained in the negative orthant. (In the latter case, the cone is pointed and generated by the extreme vectors.)

**Definition 2.** A nonzero vector **v** ∈ **C** is conformally non-decomposable if **v** = **v**¹ + **v**² for any nonzero vectors **v**¹, **v**² ∈ **C** with sign(**v**¹), sign(**v**²) ≤ sign(**v**) implies **v**¹ ∼ **v**² (that is, **v**¹ = λ**v**²).

As stated above, EFMs can be defined as the conformally non-decomposable vectors of the flux cone. Indeed,
we have the following equivalence.

**Proposition 1.** *A nonzero vector* **v** ∈ **C** *is conformally non-decomposable if and only if it is support-minimal.*

*Proof.* Assume that **v** ∈ **C** is conformally decomposable, that is, **v** = **v**¹ + **v**² for nonzero **v**¹, **v**² ∈ **C** with sign(**v**¹), sign(**v**²) ≤ sign(**v**) and **v**¹ ≁ **v**². Then also **v**¹ ≁ **v**, and there exists a largest λ > 0 such that the nonzero vector **v**' = **v** − λ**v**¹ fulfills sign(**v**') ≤ sign(**v**). For this λ, **v**' ∈ **C** (that is, N**v**' = 0 and **v**' → ≥ 0) and sign(**v**') < sign(**v**) (in particular, v'_i = 0 and v_i ≠ 0 for some i). Hence, supp(**v**') ⊂ supp(**v**), that is, **v** is not support-minimal.

Conversely, assume that **v** ∈ **C** is not support-minimal, that is, supp(**v') ⊂ supp(**v**) for a nonzero **v'** ∈ **C**.
Then, there exists a largest λ > 0 such that the nonzero vectors **v**¹ = ½**v** + λ**v'** and **v**² = ½**v** - λ**v'** fulfill
sign(**v**¹), sign(**v**²) ≤ sign(**v**). For this λ, either sign(**v**¹) < sign(**v**) or sign(**v**²) < sign(**v**); in any case,
**v**¹, **v**² ∈ **C** and **v**¹ ≁ **v**². Clearly, **v** = **v**¹ + **v**², that is, **v** is conformally decomposable. □

Conformally non-decomposable vectors are closely related to extreme (or non-decomposable) vectors.

**Definition 3.** A nonzero vector **v** ∈ **C** is extreme if **v** = **v**¹ + **v**² for any nonzero vectors **v**¹, **v**² ∈ **C** implies **v**¹ ∼ **v**².

If the flux cone is contained in the non-negative orthant (in particular, after reaction splitting, as described in
Section 4.3.2) and hence is pointed, EFMs can be defined as the extreme vectors.

**Proposition 2.** Let $\mathcal{C} \subseteq \mathbb{R}^n$. A nonzero vector $\mathbf{v} \in \mathcal{C}$ is extreme if and only if it is conformally non-decomposable.

*Proof.* If **v**, **v**¹, **v**² ∈ **C** ⊆ ℝ≥ⁿ, then **v** = **v**¹ + **v**² implies sign(**v**¹), sign(**v**²) ≤ sign(**v**), and Definitions 2 and 3 agree. □

proposed to guide the genetic manipulation of microorganisms to perform desirable properties such as
synthesis of a bio-compound or efficient production of a recombinant protein (e.g. [98, 99]). From a more
theoretical point of view, EFMs have also been used in attempts to quantify cellular robustness [100], in
particular regarding robustness under genetic perturbations [101]. The relevance of elementary flux mode
analysis to cellular robustness stems from the fact that there is rarely a unique conical combination of
elementary flux modes for any given flux vector, which implies there are multiple combinations of minimal
metabolic pathways to achieve the same desired effect. This redundancy can be interpreted as a measure
for the metabolic robustness of an organism, in terms of preserving essential metabolic functionalities
under loss of a gene, for example.

There have also been several ways that EFM analysis has been incorporated into analysis of multi-omics
data. For example, on the basis of transcriptomic profiling of microorganisms, metabolic pathways asso-
ciated with elementary flux modes have been scored according to their probability of carrying flux [102].
The principle here is that, although levels of RNA often serve as a poor proxy for flux through the reaction
associated with that particular enzyme's gene, by creating a gene set associated with an entire EFM there
might be a better chance of concretely assessing whether the metabolic pathway as a whole is likely to carry
flux. The study [102] suggested that the integration of EFM analysis with gene expression data enabled
the identification of certain metabolic pathways activated during stress conditions, and that the organiza-
tion of elementary flux mode utilization in *Saccharomyces cerevisiae* involves a disparate combination of
<figure><img src="image_37.png" /><figcaption>Figure 4.7: Adapted from [187]</figcaption></figure>

knockouts specifically, *minimal cut sets* may be used. Finally, we discussed that the flux cone can be sampled randomly to characterize the flux cone, if this characterization does not need to be exhaustive.

In many cases we have additional information that determines that part of the flux cone is infeasible. For example, some metabolic fluxes may have been measured so that these reaction rates can be fixed to their observed value. In other cases, one may want to use thermodynamic properties to prohibit reactions from occurring that would violate the second law of thermodynamics. These additional constraints can be imposed on top of the mass-balance constraint to further bound the space of feasible flux distributions; each correctly-imposed constraint narrows down the space of feasible fluxes, and thus increases our knowledge of the metabolic state of the cell.

All explorations of the space of feasible flux distributions show one unavoidable conclusion: the metabolic network is incredibly flexible. Even when several constraints are imposed, a genome-scale metabolic model will allow for an almost incomprehensible number of modes in which the metabolic network can function. Consequently, to predict the metabolic state of a cell in more detail we need to make an additional assumption. In the following chapter, we will study what predictions we can make when we assume that the metabolic state is optimized to perform a certain function.

## Recommended readings

A nice paper, called "Elementary flux modes in a nutshell: Properties, calculation and applications", which gives an understandable introduction of elementary flux mode analysis and its applications, can be found here: doi.org/10.1002/biot.201200269

*Elementary Flux Vectors* were introduced as an analog of *Elementary Flux Modes* in the case that the flux mode is further bound by at least one inhomogeneous constraint. A nice review of these EFVs is can be found in [185].

## Problems

**Problem 4.1** *Spirallus insilicus*, a completely fictional organism [186], is characterized by the metabolic network depicted in Figure 4.7:

$X$, $S$ and $P$ represent the biomass, one substrate and one product, while metabolites A to E denote intracellular metabolites. One directional arrows indicate irreversible reactions (all but $v_4$)
highly specialized and multi-tasking roles. Beyond transcriptomic profiling, isotope tracing experiments in principle provide a much more direct insight into quantifying metabolic flux. To interpret isotope tracing data, an extension of the concept of an EFM was introduced in [103].

### 4.3.2 Reaction splitting for EFM computation

The computation of EFMs via the double description (DD) method as well as the solution of linear programs (LPs) via the simplex algorithm assume that the flux cone is given in certain standard forms. (Note, however, that the computation of EFMs via lexicographic reverse search (Irs) does not involve such an assumption.)

Recall that the flux cone is given by the mass-balance constraints $N\mathbf{v} = 0$ and the irreversibility constraints $\mathbf{v} \ge 0$, whereas standard forms are given by $A\mathbf{v} = 0$ (for DD) or $A\mathbf{v} \ge 0$ and $\mathbf{v} \ge 0$ (for LP) with a matrix $A$ of appropriate dimensions. To bring the flux cone into standard form, we will split reversible reactions into irreversible forward and backward reactions. First, we order reactions such that

$$N\mathbf{v} = (N \to \ N \rightleftharpoons \ \begin{pmatrix} \mathbf{v} \\ \mathbf{v} \end{pmatrix}), \qquad (4.8)$$

where the superscripts $\to$ and $\rightleftharpoons$ refer to the irreversible and reversible reactions, $R \to$ and $R \rightleftharpoons$, respectively. Next, for every reversible reaction $i \in R \rightleftharpoons$ with net reaction rate $v_i$, we define a forward reaction with "rate" $w_i \to \ge 0$ and a backward reaction with "rate" $w_i \leftarrow \ge 0$ such that $v_i = w_i \to - w_i \leftarrow$. (In vector form, $\mathbf{v} \rightleftharpoons = \mathbf{w} \to - \mathbf{w} \leftarrow$.) Note that the "rates" $w_i \to$, $w_i \leftarrow$ do not denote the (microscopic) forward and backward reaction rates $v_i \to$, $v_i \leftarrow$ that determine the net reaction rate $v_i = v_i \to - v_i \leftarrow$. They are auxiliary quantities, and only their difference $w_i \to - w_i \leftarrow = v_i$ has a biochemical meaning (and is the subject of constraint-based metabolic modeling). Further, for every irreversible reaction $i \in R \to$, we write $v_i = w_i \to$ to obtain a uniform notation. (In vector form, $\mathbf{v} \to = \mathbf{w} \to$.) Now,

$$N\mathbf{v} = (N \to \ N \rightleftharpoons \ -N \rightleftharpoons) \begin{pmatrix} \mathbf{w} \\ \mathbf{w} \\ \mathbf{w} \end{pmatrix}. \qquad (4.9)$$

By introducing the augmented stoichiometric matrix $\overline{\mathbf{N}}$ and the corresponding non-negative flux vector $\mathbf{w}$, we can write

$$N\mathbf{v} = \overline{\mathbf{N}}\mathbf{w} \quad \text{with} \quad \overline{\mathbf{N}} = (N \to \ N \rightleftharpoons \ -N \rightleftharpoons), \quad \mathbf{w} = \begin{pmatrix} \mathbf{w} \\ \mathbf{w} \\ \mathbf{w} \end{pmatrix}. \qquad (4.10)$$

As a consequence, the augmented flux cone is given by

$$\bar{C} = \{\mathbf{w} \mid \overline{\mathbf{N}}\mathbf{w} = 0, \mathbf{w} \ge 0\} \qquad (4.11)$$

or, in LP standard form, by

$$\bar{C} = \{\mathbf{w} \mid A\mathbf{w} \ge 0, \mathbf{w} \ge 0\} \quad \text{with} \quad A = \begin{pmatrix} \overline{\mathbf{N}} \\ -\overline{\mathbf{N}} \end{pmatrix}, \qquad (4.12)$$

after writing equations as non-strict inequalities.

Obviously, $\bar{C}$ is contained in the non-negative orthant and hence pointed. As an important consequence, EFMs can be defined as the extreme vectors of the flux cone, see the 'Math box', and be computed by algorithms based on the DD method.

For examples of pointed polyhedral cones (in the non-negative orthant), see Figure 4.2. Note that the cone in Figure 4.2 (A) is not an s-cone and hence not a flux cone. In particular, its generators/extreme vectors lie in the interior of the non-negative orthant. On the other hand, the cone in Figure 4.2 (B) is an s-cone. Its generators/support-minimal vectors/EFMs arise from the intersection of a subspace (the nullspace of
1. How many intracellular metabolites, intracellular reactions and transport reactions are involved in the model?

2. Obtain the stoichiometric matrix ($N$) and the vector of fluxes. How many elements are in the product $Nv$ and what do they represent?

3. Is the matrix $N$ of full rank? How many fluxes should be specified to have a unique solution?

4. Transform the set of constraints so that they define a pointed cone. Determine the number of variables (fluxes) and constraints.

**Problem 4.2** Consider the following small metabolic network:

$$ 
\begin{aligned}
S_e &\xrightarrow{v_0} S_c \\
S_c &\xrightarrow{v_1} P_c \\
P_c &\xrightarrow{v_2} C_c \\
P_c &\xrightarrow{v_3} D_c \\
P_c + 2C_c &\xrightarrow{v_4} X
\end{aligned}
 $$

Metabolites with a $c$ subscript are located in the cytosol (intracellular) while $e$ stands for extracellular and $X$ represent biomass. All fluxes are positive.

1. Represent the model as a reaction network (a sketch with metabolites and reactions)

2. Obtain the stoichiometric matrix ($N$) and list the variables of the metabolic model ($v$)

3. Show that there is no solution to the mass balance equation $Nv = 0$ producing metabolite $D$. Identify why this is so and modify the model so the production of $D$ is allowed ($v_3 > 0$)

**Problem 4.3** Assume reaction $v_4$ is irreversible from $A$ to $D$ in *Spirallus insilicus* (Problem 4.1). Calculate all the Elementary Flux Modes.

1. By hand.

2. Using a software of your choice (e.g. https://pypi.org/project/efmtool/)

**Problem 4.4 Elementary Flux Modes** Consider the following metabolic network

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

with stoichiometric matrix

$$ N = \begin{pmatrix} 1 & -2 & 0 & 0 & 0 & 0 & -2 \\ 0 & 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 2 & 0 & -1 & -2 & 0 & 0 \\ 0 & 0 & 1 & 0 & 1 & -1 & 1 \end{pmatrix} $$

1. Red points denote carbon atoms. Check that carbon atoms are conserved in all reactions.

2. All metabolites are treated as internal, that is, they need to be mass-balanced. Find all EFMs (by pure reasoning or by using a software). Determine all EFMs in which all fluxes are in forward direction, i.e. along the “conventional
<figure><img src="image_50.png" /><figcaption>Figure 4.2: Pointed polyhedral cones. (A) A pointed polyhedral cone that is not a flux cone; all its generators lie in the interior of the non-negative orthant. (B) A pointed polyhedral cone that is a flux cone; its generators arise from the intersection of a subspace with the boundaries of the non-negative orthant.</figcaption></figure>

the stoichiometric matrix) with the boundaries of the non-negative orthant.

Again, we return to the simple representation of central carbon metabolism presented in Figure 4.1. After reaction splitting, the mass-balance constraint can be written as

$$ \bar{\mathbf{N}} \mathbf{w} = \begin{pmatrix} 1 & -1 & 0 & 0 \\ 0 & 2 & -1 & -1 \end{pmatrix} \begin{pmatrix} w_1 \\ w_2 \\ w_3 \\ w_4 \end{pmatrix} = \mathbf{0}. \qquad (4.13) $$

In particular, the reversible fourth reaction has been split into irreversible forward and backward reactions with reaction vectors $(\begin{pmatrix} 0 \\ -1 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix})$ and "rates" $w_4^-, w_4^-$, see Equation (4.5). Now, algorithms based on the DD method can be applied to the mass-balance and irreversibility constraints in standard form, $\bar{\mathbf{N}} \mathbf{w} = \mathbf{0}$ and $\mathbf{w} \ge \mathbf{0}$. As it turns out, the set of EFMs (support-minimal vectors) is given by

$$ g^{(1)} = \begin{pmatrix} 1 \\ 1 \\ 0 \\ 2 \\ 0 \end{pmatrix}, \quad g^{(2)} = \begin{pmatrix} 1 \\ 1 \\ 2 \\ 0 \\ 0 \end{pmatrix}, \quad g^{(3)} = \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \\ 1 \end{pmatrix}, \quad \text{and} \quad g^{(4)} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \\ 1 \end{pmatrix}. \qquad (4.14) $$

EFMs $g^{(1)}, g^{(2)}, g^{(3)}$ correspond to EFMs $f^{(1)}, f^{(2)}, f^{(3)}$ before reaction splitting, see Equation (4.6). Just recall $v_4 = w_4^- - w_4^-$ . However, EFM $g^{(4)}$ corresponds to zero flux. More specifically, it represents the fourth reaction having equal forward and backward "rates" and hence zero net reaction rate. Such EFMs are artifacts of reaction splitting and need to be discarded when translating the EFMs of the augmented flux cone back to the EFMs of the original flux cone.

## 4.4 Extra constraints and flux polyhedra

### 4.4.1 Inhomogeneous linear flux constraints

We have so far been working exclusively with mass-conservation and irreversibility constraints, which are captured entirely by the stoichiometric matrix where each row is associated with a metabolite concentration at steady state. We also saw that these considerations alone result in a flux cone that is by definition unbounded, meaning that a flux vector in this space is allowed to take on any absolute value (i.e. multiplying a flux vector in the flux cone by an arbitrarily large positive number again returns a flux vector in
directions" indicated by arrows.

3. Which of the EFMs are thermodynamically realizable? Explain why.
the flux cone). However, there are physical constraints limiting the magnitude of flux vectors, especially on the values of flux through exchange reactions that may depend on concentrations of extracellular substrates, numbers of transporter molecules in the membrane, or for which we might have direct experimental measurements. Typically, such bounds on flux values are imposed using inequality constraints of the form $v_i^{lb} \le v_i \le v_i^{ub}$ where $v_i^{lb}$ and $v_i^{ub}$ are lower and upper bounds, respectively, for the flux through the $i$th reaction. When reactions have been decomposed into forward and reverse directions, both upper and lower bounds are non-negative where the latter is usually zero.

In the example from Figure 4.1, one may impose an upper bound on the flux value $v_1$, suggesting that there is a maximal rate at which the cell or organism can import glucose from the extracellular environment. In this case the total set of constraints on the flux vector $v$ take the form

$$N v = 0, \quad v \ge 0, \quad v_1 \le v_1^{ub}, \qquad (4.15)$$

where $v_1^{ub}$ is the maximal glucose uptake rate. It is important to note that the new constraint is of a different kind than the mass-balance and irreversibility constraints: the right-hand side of the constraint is nonzero. Constraints that involve a nonzero are called *inhomogeneous constraints*. We can write these constraints in matrix form as

$$G v \ge h, \qquad (4.16)$$

where in this particular case

$$G = (-1 \quad 0 \quad 0 \quad 0), \quad h = (-v_1^{ub}). \qquad (4.17)$$

In general, the matrix $G$ will have $\ell$ rows corresponding to $\ell$ inhomogeneous linear constraints of the form

$$\sum_i G_{ji} v_i \le h_j, \quad j = 1, \dots, \ell. \qquad (4.18)$$

That is, for constraint $j$, there are $n$ entries $G_{ji}$ ($i = 1, \dots, n$) of the matrix $G$ and the component $h_j$ of the $\ell$-dimensional vector $h$. Many constraints can be written in this general form. For example, after reaction splitting, one may impose a bound on the total flux that a cell can catalyze, by setting all entries (in the corresponding row of $G$) to 1.

Altogether, the constraints on $v$ define a *flux polyhedron* that is necessarily contained within the flux cone given by the homogeneous constraints $N v = 0$ and $v \ge 0$. The additional inhomogeneous constraints serve to further restrict the cone such that various (if not all) dimensions become bounded, thus bounding the total magnitude of the flux vector $v$.

### 4.4.2 From the flux polyhedron to the EFM weight polyhedron

Via the conical sum $v = \sum_{k=1}^{\ell} \lambda_k f^{(k)}$, constraints on the fluxes $v$ define constraints on the EFM weights $\lambda$ and hence a corresponding *EFM weight polyhedron*. Whereas elements $v$ of the flux polyhedron have entries $v_i$ for every reaction $i \in \mathcal{R}^-$, elements $\lambda$ of the EFM weight polyhedron have entries $\lambda_k$ for every EFM $f^{(k)}$, $k = 1, \dots, \ell$ (and hence can be very high-dimensional).

In the example from Figure 4.1, let $\lambda_1, \lambda_2, \lambda_3 \ge 0$ be the weights of the corresponding EFMs in the (conformal) sum $v = \sum_{i=1}^3 \lambda_i f^{(i)}$. Bounding the extracellular glucose uptake rate puts an upper bound on the weights of EFMs $f^{(1)}, f^{(2)}$ (involving the glucose uptake reaction),

$$\lambda_1 + \lambda_2 \le v_1^{ub}, \qquad (4.19)$$

see also Figure 4.1 (B). However, the weight of EFM $f^{(3)}$ (associated with uptake and oxidation of the fermentation product) can remain unbounded.

For this simple example, it is quite straightforward to interpret the geometric consequences of the maximal glucose uptake rate. Any flux vector $v$ in the resulting flux polyhedron now corresponds to a point $(\lambda_1, \lambda_2)$ in the (projected) EFM weight polyhedron depicted in Figure 4.3 (A). However, the weight $\lambda_3$ remains free, and the (full) EFM weight polyhedron is depicted in Figure 4.3 (B). In terms of the flux polyhedron, the
<figure><img src="image_51.png" /><figcaption>Figure 4.3: Feasible regions in the space of EFM weights - (A) Possible combinations of EFM weights $\lambda_1$ and $\lambda_2$, given by the inequality $\lambda_1 + \lambda_2 \le v_1^{ub}$ (and $\lambda_1, \lambda_2 \ge 0$). (B) Geometry of the EFM weight polyhedron (blue) representing any flux vector that satisfies the mass-balance, irreversibility, and maximal glucose uptake rate constraints. While bounded in $\lambda_1, \lambda_2$, it is unbounded in $\lambda_3$.</figcaption></figure>

maximal glucose uptake has restricted the flux cone along $v_1, v_2$ while leaving $v_3, v_4$ unbounded.

In order to obtain a bounded flux polyhedron (a flux polytope), we impose an upper bound on the uptake rate of the fermentation product, that is, $-v_4 \le v_4^{ub}$. In terms of the EFM weights, we obtain the bound $-2\lambda_1 + \lambda_3 \le v_4^{ub}$. Since conformal sums are sufficient to generate the flux cone, this simplifies to $\lambda_3 \le v_4^{ub}$. Altogether, the EFM weight polyhedron is given by

$$ \lambda_1, \lambda_2, \lambda_3 \ge 0, \quad \lambda_1 + \lambda_2 \le v_1^{ub}, \quad \lambda_3 \le v_4^{ub}. \qquad (4.20) $$

Indeed, all EFM weights and hence all fluxes are bounded.

More general constraints, for larger metabolic reaction networks will be more difficult to interpret and visualize in such simple geometric terms. Quite quickly the combinatorial complexity associated with combinations of multiple constraints and EFMs will become unmanageable. The intuitive treatment of inhomogeneous linear constraints is partially assisted using the concept of elementary flux vectors on which we will add a section in a later version of this book, but both geometrically and biologically these objects are nowhere near as easy to interpret as their EFM counterparts. We shall see that alternative computational methods for exploring flux space therefore become imperative.

As a final remark, we clarify once more that the general form of constraints (4.16) is by no means restricted to sums on the left hand side that involve just a single reaction and can of course include constraints on weighted sums of flux values for different reactions. These weighted sums are often associated with particular biological interpretations: in the example from Figure 4.1, we might want to restrict our search of flux space to those flux vectors $v$ that produce adenosine triphosphate (ATP) at a rate of at least $v^{ATP}$. Although a more elaborate model would of course include ATP as one of the metabolites, in this example we can use our biological understanding of central carbon metabolism to see that ATP is produced in reactions $v_2$ and $v_3$. A lower bound on ATP production would thus be a lower bound on a combination of $v_2$ and $v_3$ with coefficients determined by stoichiometry (depending on the organism under investigation). We could write such a constraint as

$$ \alpha_1 v_1 + \alpha_3 v_3 \ge v^{ATP} \qquad (4.21) $$

with appropriate coefficients $\alpha_1, \alpha_3$. Such a constraint forms an additional row of the matrix $G$ and we leave it as an exercise for the reader to explore how this affects the geometry of the flux polytope for various values of the coefficients, minimal ATP production rate and maximal glucose and fermentation product uptake rates. Particular combinations of constraints will be impossible to satisfy simultaneously (i.e. when the minimal rate of ATP production is impossible to achieve under the given bounds on glucose and fermentation product uptake rates), resulting in a flux polytope that is empty. In such cases the set of
# Chapter 5

## Optimal solutions of constraint-based
metabolic models

Daan de Groot, Felipe Scott, and David Tourigny

### Chapter highlights

* An optimization objective can be added to constraint-based models to make more specific predictions.

* Different purposes can be served by choosing different optimization objectives and constraints

* The optimal solutions can be understood in terms of elementary flux modes

### 5.1. Can we use an optimality assumption to predict metabolic behavior?

In the previous chapter we began to characterize an organism's metabolism based on the individual chemical reactions it can potentially catalyze by expressing enzymes encoded within its genome. To understand how the genome constrains patterns of metabolic flux we needed to make several simplifying assumptions. The first important assumption was that intracellular metabolism is at steady-state, i.e., that the production and consumption of all metabolites is balanced such that their concentrations are constant in time. These resulted in the mass-balance constraints on the flux vector v. The 'flux cone' of all flux vectors satisfying the mass-balance constraints could be further reduced by additional constraints on v, based on extra physical and biological assumptions about the magnitude and directionality of certain reactions within the network. We introduced several ways in which the entire flux space could be explored.

When applied to very large metabolic networks, the flux space will often contain an infinite number of flux vectors v that simultaneously satisfy all constraints. From a mathematical perspective, this implies that the constraints do not include enough information to uniquely specify a flux vector v. This makes sense biologically, since if we imagine constraints are related to experimental observations it is very unlikely that we will ever be able to make enough to fully account for every reaction encoded within the entire genome of an organism (no matter how simple it might be). Often, however, researchers do want to further narrow down the set of flux vectors that they think biologically relevant to the organism and conditions they are studying, perhaps even to a unique v imagined to describe the metabolic state of an organism at a given moment in time. One popular approach for doing so is to provide an additional assumption (or set thereof) in the form of an objective function: it is assumed that the metabolic state of an organism is such that some function of v (e.g. growth rate) is maximized to satisfy some criteria (e.g. evolutionary selective advantage). The computational problem then becomes one of constrained-optimization: find a flux vector v that is optimal in terms of the objective function(s) that simultaneously satisfies all constraints. The resulting space of optimal flux vectors (sometimes containing just one unique vector) is often considerably smaller than the space of those that satisfy only the constraints.
constraints on v are called infeasible.

### 4.4.3 Thermodynamic constraints

In Chapter 3 the basic concepts of chemical thermodynamics were introduced, in particular, the Gibbs free energy of a metabolic reaction was defined in terms of the concentrations of its products and substrates. For a metabolic reaction network with stoichiometric matrix N, the vector of Gibbs free energies (one for each reaction in the network) Δ<sub>r</sub>G' can be written in matrix form as

$$ \Delta_{\mathrm{r}}G' = \Delta_{\mathrm{r}}G'^{\circ} + RT \cdot N^{\top} \ln(s) \qquad (4.22) $$

where R is the gas constant, T the temperature and s the vector of metabolite concentrations at steady state. The components of the vector Δ<sub>r</sub>G'° are the changes in standard Gibbs free energy for each corresponding reaction. Typically, these values are not known precisely for reactions in the network, but can be estimated or approximated from experimental data using methods beyond the scope of this chapter. Similarly, although it is often difficult to accurately measure all metabolite concentrations, in principle the vector s can be obtained experimentally. However, in practice experimental data on s and Δ<sub>r</sub>G'° are almost never available. Various methods have therefore been developed to combine estimation of Δ<sub>r</sub>G'° (sometimes with partial measurements of s) with advanced computational techniques that allow simultaneous optimization (see next chapter) or sampling (see below) of v and s (or equivalently: Δ<sub>r</sub>G').

The second law of thermodynamics applied to chemical reaction networks can be summarized by saying that every component of the metabolic flux vector v must satisfy the condition

$$ \mathrm{sign}(v_i) = -\mathrm{sign}(\Delta_{\mathrm{r}}G'_i) \qquad (4.23) $$

where v<sub>i</sub> and Δ<sub>r</sub>G'<sub>i</sub> are the ith components of v and Δ<sub>r</sub>G', respectively, and sign(x) denotes the sign of a variable x, and sign(0) = 0. It is important to point out that this notation is different to that used previously, where we had assumed all v<sub>i</sub> to be non-negative by decomposing each reaction into irreversible forward and backward reactions. Returning to this reversible notation simplifies the inclusion of thermodynamic constraints into constraint-based models and also their interpretation. According to the second law, a reaction can only proceed in a direction where the change in Gibbs free energy is negative. Thus, to be consistent with mass-balance and the second law of thermodynamics, a flux vector v must simultaneously satisfy both (4.1) and (4.23), with Δ<sub>r</sub>G' defined in (4.22). The consequence of these additional constraints on the geometry of the space of metabolic flux distributions is to exclude quadrants incompatible with the signs of Δ<sub>r</sub>G'. Equivalently, imposing the second law of thermodynamics on metabolic flux distributions removes regions of the space that are associated with combinations of thermodynamically-infeasible reaction directionalities.

The resulting space of feasible flux vectors is almost always non-convex, which means more advanced computational methods are required to explore it efficiently. The intuitive reason for this is that imposing thermodynamic constraints on top of the mass-balance constraint is usually done in terms of Boolean variables, which breaks the linearity of the problem that we had and exploited so far. Relating this to the EFMs that were discussed previously, it for example becomes clear that any EFM representing an internal cycle—not including any exchange reactions—will never be thermodynamically feasible. Thus, thermodynamic constraints also reduce the set of EFMs that are possible in a metabolic network. Interestingly, it turns out that any thermodynamically-feasible metabolic flux vector can be expressed solely in terms of thermodynamically-feasible EFMs [104], but the converse statement is not true: a linear combination of thermodynamically-feasible elementary flux mode does not necessarily satisfy the thermodynamic constraints. This shows how the workable properties of convex spaces break down as the mathematical models become more complex, in this case by accounting for thermodynamics.

### 4.4.4 Computational challenges for EFM analysis

Enumerating EFMs for large networks can be computationally challenging if not impossible. In principle, EFMs can be found by removing one reaction at a time and solving the resulting mass-balance constraint
In this chapter, we will study metabolic models based on constrained-optimization. We will introduce a selection of commonly used objective functions and the computational methods used to solve the associated constrained-optimization problem. We will also characterize optimal solutions that we get in terms of the minimal metabolic strategies that we identified in the previous chapter: elementary flux modes. Finally, we will explain how we can handle the cases where the solutions are, even after optimization, not unique.

## 5.2. Metabolic models based on linear optimization problems

In the previous chapter, we described how linear homogeneous and inhomogeneous constraints arising from biological and physical knowledge can be combined into matrix and vector notation and written in the general form presented in Equations (4.11) and (4.12). The resulting space of all flux vectors $v$ satisfying these constraints is called the flux polyhedron. The flux polyhedron can remain high-dimensional and, as explained above, an objective function $f$ can be used to narrow down the set of flux vectors to only those that are optimal (i.e., maximize the objective function). The general form in which we can write the resulting constraint-based optimization problem is therefore:

$$ \max_{v} f(v), \text{ such that } Av \geq b, \qquad (5.1) $$

with

$$ A = \begin{pmatrix} N \\ -N \\ I \\ G \end{pmatrix}, \quad b = \begin{pmatrix} 0 \\ 0 \\ 0 \\ h \end{pmatrix}. \qquad (5.2) $$

Recall that $Nv = 0$ models the steady-state assumption, while the multiplication with the identity matrix $(I_{n \times n} v \geq 0)$ captures the fact that we forced all reactions to be irreversible by splitting reversible reactions into a forward and a backward reaction. Finally, $Gv \geq h$ can be used to impose additional 'inhomogeneous' constraints that can be used to input additional biological knowledge such as an experimentally measured upper bound on the uptake rate of a certain nutrient.

In many cases, the objective function is chosen to be a linear function of the fluxes, i.e.,

$$ f(v) = \sum_{i} c_{i}v_{i}, \qquad (5.3) $$

where coefficients $c_i$ weigh the relevance of the different reaction rates in the objective function. Problems of the form (5.1), (5.2), and (5.3) in general are called *linear programming problems* and as the name suggests can be solved using *linear programming*. Applied to metabolic models, linear programming is called *Flux Balance Analysis* (FBA). Linear programming problems are well studied, such that FBA is perhaps the most popular approach to genome-scale metabolic models [188, 189]. FBA problems are relatively easy to solve using specialized optimization software, which have been highly developed due to the general applicability of linear programming in economics, logistics, and many other fields also. In the following subsections we will briefly describe various choices that can be made for the linear objective function $f(v)$ in FBA.

As an example FBA problem, in Figure 5.1 we have extended the minimal example from the previous chapter to include ATP and biomass ($X$) production, assuming the latter is produced from pyruvate using a single reaction that consumes $s_X$ molecules of ATP with flux value $v_X$. We also introduce as a linear objective function the total rate of ATP production, $v_{ATP}$. Since in this example, reactions $v_1$ and $v_3$ produce ATP with stoichiometric coefficients $s_1$ and $s_3$, respectively, the total rate of ATP production is given by $v_{ATP} = s_1v_1 + s_3v_3 - s_Xv_X$. The FBA problem is then given by simply maximizing $v_{ATP}$ subject to $v_0, v_{O_2}, v_1, v_2, v_3, v_4, v_X$ satisfying the mass-balance constraints but, as we will see in the next subsection, this would result in a problem that is unbounded: the flux vectors and resulting optimal value of $v_{ATP}$ could be indefinitely large. Biologically, this is because there are no bounds on the uptake rates of glucose
problem until it is no longer possible to remove a reaction and still obtain a flux vector that satisfies the steady state conditions. However, the equivalence of EFMs and extreme vectors of the flux cone (after re-action splitting) described in Section 4.3.2 enables the use of algorithms that are specialized in the efficient enumeration of extreme rays of polyhedral cones, such as the double description method [105]. Various tools have been developed for elementary flux mode enumeration based on this algorithm (e.g. EFMTOOL [106] or MetaTool [107]). However, when the size of the metabolic reaction network grows, the number of EFMs scales disproportionately, leading to a combinatorial explosion that effectively makes enumeration impossible for genome-scale networks containing several thousands of reactions [108]. Currently, EFM analysis is therefore restricted to medium-scale reconstructions containing on the order of several hundreds of reactions, and results in the identification of several hundred million EFMs (e.g. enumeration based on the *Escherichia coli* core model results in approximately 272 million EFMs).

Approaches to reduce the complexity of dealing with so many EFMs even for metabolic reaction networks of modest size have also been proposed. These include invoking transcriptional regulatory constraints to eliminate most of the EFMs to be considered in downstream analysis. Imposing additional constraints based on thermodynamic conditions similarly reduces the set of EFMs considerably. A problem with these approaches is evidently that they do still depend on an initial calculation of all EFMs, and so do not solve the problem of enumeration complexity. A rigorous study of the complexity of EFM mode enumeration was performed by Acuña and colleagues [109]. They showed that the decision problem if there exists an EFM containing two specific reactions is NP-complete whilst the complexity of enumerating all EFMs remains open.

Later in this chapter we will explore some alternatives to EFM enumeration that reduce the difficulty of enumeration, cf. Section 4.5.

### 4.4.5 Reducing combinatorics of EFMs computation

In order to reduce the combinatorics of EFM computation to a feasible order, the search space may be limited to the biologically relevant EFMs only. This can be done by considering additional biological constraints before, during, or after the computation of EFMs. One way to restrict the search space is to remove all 'irrelevant' reactions in a metabolic network, that is,

* reactions that are not essential for the cell (not part of the core metabolism),

* reactions that are not performed for *chemo-physical*, *kinetic*, or *thermodynamic* reasons,

* reactions that are too expensive in terms of enzymatic resource allocation,

* reactions that transport metabolites which are not present in the growth medium ('environmental regulation'),

* reactions that are catalyzed by enzymes whose expression is inhibited by transcriptional regulation.

The purpose of incorporating biological constraints, from the perspective of a modeler, is to reduce the number of pathways the biologist needs to analyze. Additionally, the computation of EFMs becomes much more efficient because fewer solutions need to be computed.

Below we are going to illustrate the last two types of constraints: environmental and transcriptional regulation. Both types can be expressed using Boolean constraints. A Boolean constraint is a Boolean function $f: \mathbb{B}^k \to \mathbb{B}$, where $\mathbb{B} = \{0, 1\}$, which takes in $k$ Boolean inputs $z \in \mathbb{B}^k$ and produces a Boolean output $b \in \mathbb{B}$ such that $b = f(z)$. In our case, Boolean functions determine whether reactions are allowed or not in EFMs based on biological conditions. To this end, reactions are associated with a Boolean indicator. The value of this indicator (either 0 or 1) determines whether that reaction can participate in an EFM.

The following relationship, for a set of reactions $R$ with corresponding fluxes $v$ and indicators $z$, determines how Boolean regulation affects the presence of reactions in EFMs:

$$ \forall r \in R : (z_r = 0) \implies (v_r = 0) $$

As an example, we consider the following small metabolic model from [111, 110], which involves transcriptional and environmental regulation. The network contains 18 reactions and 18 metabolites (10 internal
<figure><img src="image_39.png" /><figcaption>Figure 5.1: A simple representation of the metabolic reaction network for central carbon metabolism. Extracellular glucose is imported into the cell via a reaction with flux $v_0$ and converted via intracellular glucose, G, to pyruvate, P, via the reaction with flux $v_1$ that has a stoichiometric 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_2$ imported via $v_{O_2}$, converted to an oxidative phosphorylation (OXPHOS) terminal product $P_2$ via the reaction with flux $v_3$. It can also be converted to biomass X with rate $v_X$. The reactions with flux values $v_1$ and $v_3$ produce ATP from ADP (in red) with stoichiometry $s_1$ and $s_3$, respectively, which can vary between species. The production of 1.0 grams of new cells, in a dry weight basis, requires one molecule of pyruvate and $s_X$ molecules of ATP.</figcaption></figure>

$v_0^{ub}$ and the fermentation product $v_4^{ub}$. Thus, if we re-impose these bounds as in the last chapter, the result is an FBA problem that is bounded and therefore has a finite objective value:

$$ 
\begin{array}{l}
\max_{v} v_{ATP} = s_{1}v_{1} + s_{3}v_{3} - s_{X}v_{X}, \quad \text{such that:} \\
0 = v_{0} - v_{1}, \\
0 = v_{O_{2}} - v_{3}, \\
0 = 2v_{1} - v_{2} - v_{3} + v_{4} - v_{X}, \\
v_{0}^{ub} \geq v_{0}, \\
v_{4}^{ub} \geq v_{4}, \\
v_{0}, v_{1}, v_{2}, v_{3}, v_{4} \geq 0.
\end{array}
\qquad (5.4) $$

To illustrate a particular instance of this FBA problem, we consider the very simple case where $v_4^{ub} = 0$, $v_0^{ub} > 0$ and $s_3 = s_1 = 1$. It can be checked by hand that an optimal solution is given by $v_0 = v_1 = v_2/2 = v_0^{ub}$, with $v_2 = v_4 = v_X = 0$. The optimal objective value is given by $v_{ATP} = 3v_0^{ub}$.

### 5.2.1. Types of linear objective functions used in FBA

Solving the constraint-based optimization problem of (5.1) will reduce the set of flux vectors to those that are optimal (maximize the objective function), but the biological validity of this prediction is critically dependent on the particular choice of $f$. Consequently, there has been a lot of consideration and debate among researchers working on FBA about the appropriate objective functions to use in different contexts and how best interpret the results. Below, we will provide some popular examples, but for a more systematic comparison of different objective functions we refer the reader to [193, 194, 195].

**Evolutionary justifications for objective functions: the rate of biomass production**

Objective functions are often based on evolutionary arguments: the objective is chosen to capture some proxy for the evolutionary fitness of an organism. The motivation behind this is that cells with a metabolic state that scores well on this fitness-proxy would
<figure><img src="image_52.png" /><figcaption>Figure 4.4: Metabolic model from [110]. Transcriptional regulation shown in Figure 4.5.</figcaption></figure>

and 8 external) and has 80 EFMs. For an illustration, see Figure 4.4. The formulae describing reaction stoichiometries and regulation rules are shown in Figure 4.5. This model makes for a good basis for studying the effect of Boolean constraints on a small scale: out of 80 EFMs, only 26 are consistent with the regulation in the most permissive growth medium – and even fewer are found when the growth medium gets restricted [110].

As mentioned above, we distinguish two types of Boolean functions. First, *environmental regulation* applies to uptake transporters and is automatically constructed from the defined growth medium. For example, the oxygen transport reaction can only be active if external oxygen is present in the growth medium ($\neg m_{\text{oxygen}} \Rightarrow \neg z_{\text{tox}}$). Second, *transcriptional regulation* is reconstructed from a literature review and curated by the modeler. For instance, r7 is regulated by the level of metabolite B in the cell, its enzyme cannot be expressed at the same time as B is being produced by r2 ($z_{\text{r2b}} \Rightarrow \neg z_{\text{r7}}$). Some individual constraints mimic the behavior of *E. coli*: the activation of respiration reaction r5a, r5b, rres depends on the presence of oxygen (motivated by the transcriptional factors ArcA and FNR); tc2 is deactivated when faced with carbon1, mimicking the behavior of glucose catabolite repression by CRP. Ultimately, these transcriptional and environmental constraints serve to filter out EFMs. For instance, the elementary mode {r2b, r3, r4, r5b, r8b, rres, th, tox, growth} is not consistent with regulation. Indeed, we have: $z_{\text{r5b}} \Rightarrow \neg m_{\text{oxygen}}$ and $z_{\text{tox}} \Rightarrow m_{\text{oxygen}}$, a contradiction.

Regulation Boolean constraints could be incorporated into the EFM computation by the method *regEFM-Tool*, as well as in the tools *SMTTool* and *aspefm* [112, 113, 114]. These constraints lend themselves naturally to logical encoding, making logic programming such as Answer Set Programming (ASP) well suited to this type of problem. Unlike traditional double description methods, which struggle with the combinatorial explosion of EFMs by the number of reactions and do not inherently handle regulatory constraints, ASP allows for an intuitive representation of Boolean constraints and efficient pruning of infeasible solutions early in the computation. In the simplest cases reactions that cannot respect the regulation constraints are directly deactivated in pre-processing.

Adding *environmental regulation* and restricting the analysis to a limited growth medium is crucial for reducing the computational load of the analysis. The software *regEFMTool* from Jungreuthmayer et al. was tested on Orth, Fleming and Palsson's *E. coli* core model [112, 115], a central carbon metabolic model of
Economic analogy 5.A: Linear programming and economic planning

Linear programming as an algorithmic approach to solving constrained linear optimization problems was first devel-
oped by soviet mathematician and economist Leonid Kantorovich in the 1930s [190, 191]. Kantorovich was tasked
with helping to optimize production in the soviet plywood industry, but soon discovered that the underlying prob-
lems could not be solved using analytical methods. He instead developed a method for solving linear optimization
problems using an iterative process through which a solution is continuously improved until an optimum is reached.
Kantorovich argued that this could be used to make soviet economic planning more efficient.

Soviet planning was primarily based on material balancing, which aimed to create a consistent plan with regards to the inputs and outputs of various industries. For example, the input requirement of steel consuming industries ought not to exceed steel production targets. In a balanced plan the input requirements for steel would match the production of steel. But a balanced plan is not necessarily an optimal one. There can well be several consistent plans of which some lead to higher overall production output than others. Kantorovich observed that productive resources were often not used where they could yield the greatest benefit. By using linear programming, planners could in principle calculate a plan that made the best use of economic resources and maximized production output. One of the problems that needed to be overcome by Kantorovich was that optimization always aims to optimize a singular objective function. However, there was no obvious way of measuring the output of qualitatively distinct products on a single scale. Without prior valuation of the products (for example through market prices) it is not clear whether 3 tanks and 10 trucks should be counted as more than 4 tanks and 8 trucks. Kantorovich circumvented this problem by assuming that outputs ought to be produced at given proportions. For example, it might be specified that 2 trucks ought to be produced for every tank. Linear programming can then be used to calculate the plan that maximizes output at these proportions. Unlike most contemporary economic applications of linear programming, this does not depend on a monetary objective function. So, what's being maximized is not monetary value. Instead, the objective function measures purely physical quantities (such as number of trucks or tons of steel).

In the context of economic planning, constraints are used to represent limits to available economic resources (such as fertile land). A plan that uses more resources than are available will not be feasible and must thus be excluded. Constraints can also be used to fix the proportions at which distinct outputs ought to be produced [192]. While it was first developed for economic planning, the fundamental principles of linear programming can also be applied to other problems (for example in biology).

come to dominate the cell-population because they outgrow their competitors. Proxies for fitness are in principle very
hard to choose since evolutionary fitness is mostly related to the average net reproduction rate of a cell over a very long
time [196]. Therefore, to know the metabolic objective that aligns with the maximization of fitness would require us to
know what the cell has been selected for in its evolutionary history. This is a non-trivial question, for example, is an
*E. coli* cell growing in the human gut selected for the same metabolic objective as a muscle cell in your body?

An objective that is used very often is the maximization of a biomass production rate, because this is used as a proxy for maximizing growth rate. It is indeed arguable that unicellular organisms with high growth rates are selected, since in stationary conditions these cells will come to dominate the population. Indeed, FBA models in which the biomass production rate is optimized seem to predict metabolic states reasonably well [197, 198, 199].

But what exactly do we mean by "biomass"? This is extensively discussed in the Chapter 2, but for our purposes it
is sufficient to say that it is the entirety of all components that constitute a new cell. In metabolic models, however,
"biomass" refers to all precursors that are outputs of the model and that are needed to produce a new cell. This has
two consequences. First, biomass in our model does not only consist of the components of which the cell is built, but
also of components needed to do the building itself, such as a certain amount of ATP. Second, what is contained in
biomass will depend on where we draw a line around the metabolic network - all necessary cell components that are
not inside are regarded as biomass. In practice, biomass appears in metabolic networks in the form of a virtual *biomass
reaction* that consumes all necessary precursor molecules in the right proportions and produces one unit of "standard
(A) Variables

Internal metabolites: $M = \{A, B, C, D, E, F, G, H, O_2, ATP, NADH\}$

External metabolites: $E = \{\text{carbon1}, \text{carbon2}, \text{D}_{\text{ext}}, \text{E}_{\text{ext}}, \text{F}_{\text{ext}}, \text{H}_{\text{ext}}, \text{oxygen}, \text{biomass}\}$

Boolean indicators for growth medium: $G = \{m_{\text{carbon1}}, m_{\text{carbon2}}, m_{\text{H}}, m_{\text{F}}, m_{\text{oxygen}}\}$

Reactions: $R = \{\text{r1}, \text{r2a}, \text{r2b}, \text{r3}, \text{r4}, \text{r5a}, \text{r5b}, \text{r6}, \text{r7}, \text{r8a}, \text{r8b}, \text{rres}, \text{tc1}, \text{tc2}, \text{td}, \text{te}, \text{tf}, \text{th}, \text{tox}, \text{growth}\}$

Flux vector: $v = \{v_{\text{r1}}, v_{\text{r2a}}, v_{\text{r2b}}, \dots\}$

Boolean indicators for reactions: $z = \{z_{\text{r1}}, z_{\text{r2a}}, z_{\text{r2b}}, \dots\}$

(B) Stoichiometry

Internal reactions

$$
\begin{align*}
\text{r1:} \quad & A + \text{ATP} \rightarrow B \\
\text{r2a, r2b:} \quad & B \leftrightarrow 2 \text{ ATP} + 2 \text{ NADH} + C \\
\text{r3:} \quad & B \rightarrow F \\
\text{r4:} \quad & C \rightarrow G \\
\text{r5a:} \quad & G \rightarrow 0.8 C + 2 \text{ NADH} \\
\text{r5b:} \quad & G \rightarrow 0.8 C + 2 \text{ NADH} \\
\text{r6:} \quad & C \rightarrow 2 \text{ ATP} + 3 \text{ D} \\
\text{r7:} \quad & C + 4 \text{ NADH} \rightarrow 3 \text{ E} \\
\text{r8a, r8b:} \quad & G + \text{ATP} + 2 \text{ NADH} \leftrightarrow H \\
\text{rres:} \quad & \text{NADH} + \text{O}_2 \rightarrow \text{ATP}
\end{align*}
$$

Transport reactions

$$
\begin{align*}
\text{tc1:} \quad & \text{carbon1} \rightarrow A \\
\text{tc2:} \quad & \text{carbon2} \rightarrow A \\
\text{td:} \quad & D \rightarrow D_{\text{ext}} \\
\text{te:} \quad & E \rightarrow E_{\text{ext}} \\
\text{tf:} \quad & F_{\text{ext}} \rightarrow F \\
\text{th:} \quad & H_{\text{ext}} \rightarrow H \\
\text{tox:} \quad & \text{oxygen} \rightarrow \text{O}_2 \\
\text{growth:} \quad & C + F + H + 10 \text{ ATP} \rightarrow \text{biomass}
\end{align*}
$$

(C) Regulation

Transcriptional

$$
\begin{align*}
z_{\text{r2b}} &\implies \neg z_{\text{r2a}} \\
\neg m_{\text{oxygen}} &\implies \neg z_{\text{r5a}} \\
m_{\text{oxygen}} &\implies \neg z_{\text{r5b}} \\
z_{\text{r2b}} &\implies \neg z_{\text{r7}} \\
m_{\text{H}} &\implies \neg z_{\text{r8a}} \\
\neg m_{\text{oxygen}} &\implies \neg z_{\text{rres}} \\
m_{\text{carbon1}} &\implies \neg z_{\text{tc2}}
\end{align*}
$$

Environmental

$$
\begin{align*}
\neg m_{\text{carbon2}} &\implies \neg z_{\text{tc2}} \\
\neg m_{\text{carbon1}} &\implies \neg z_{\text{tc1}} \\
\neg m_{\text{H}} &\implies \neg z_{\text{th}} \\
\neg m_{\text{F}} &\implies \neg z_{\text{tf}} \\
\neg m_{\text{oxygen}} &\implies \neg z_{\text{tox}}
\end{align*}
$$

Figure 4.5: Formulae for the metabolic model from Figure 4.4. Stoichiometry is given for reactions and metabolites; simple arrow or double arrow represent reversibility, for instance reaction r1 consumes one A and one ATP to produce one B. The names r2a, r2b and r8a, r8b denote the forward and backward directions of the respective reactions, while r5a, r5b represent isozymes. Boolean inputs for the Boolean functions can be either growth medium metabolites or reactions.

95 reactions containing a complete transcriptional regulation network. The analysis was performed with all uptake reactions allowed. The total number of EFMs was reduced from 226.3 million to 2 million EFMs after post-processing.

Using aspefm, Mahout *et al* applied environmental regulation, transcriptional regulation, as well as thermo-dynamic constraints in order to further reduce that set to a subset of only 10³ EFMs for post-processing analysis of optimal uptake rates [114, 116]. In general, we therefore recommend to routinely incorporate basic regulation constraints checking in order to drastically reduce the complexity of search of EFMs on metabolic models. This is particularly true for genome-scale models, which number of reactions reach thousands and number of EFMs reach billions. Ideally the procedure should be done in pre-processing,
biomass". Maximizing the biomass production rate thus takes the very simple form of just maximizing the rate through the biomass reaction.

The use of such a fixed biomass reaction represents an important assumption, because in reality the biomass composition will be condition dependent. For example, if a cell grows faster and contains more ribosomes, this increases the cellular fraction of proteins and polynucleotides, and hence the need for the respective precursors (amino acids and nucleotides). Moreover, biomass composition can even depend on the choice of metabolic strategy. If a pathway includes enzymes that contain a lot of iron, then depending on the flux solution (which uses this pathway or not), more or less iron will be contained in the biomass. So, the flux solution must be known to know the biomass composition, but the biomass composition must also be known to get to a flux solution. To resolve this, we would need a model of the entire cell, including the synthesis reactions of all enzymes. Such models will be discussed later, in the Chapter 9 on large cell models.

**Evolutionary justifications for objective functions: alternative fitness-proxies** In some cases, modeling the maximization of the instantaneous growth rate through the biomass reaction is an unrealistic proxy of the evolutionary fitness. For example, in multicellular organisms each cell performs a task that contributes to the fitness of the whole organism, but this is not related to the reproduction rate of the individual cells. In those cases, we may still try to capture an evolutionary objective when we know the main task of the cell-type. For example, beta-cells in the pancreas have as their main task to produce insulin, and we may thus model their metabolism by maximizing the production of insulin.

In other cases, our metabolic model is focused only on a very small part of the true metabolic network, and therefore does not model the production of all biomass precursors. In such cases, energy production rate in the form of ATP production rate is often maximised. Yet other objective functions that are sometimes used and have a (somewhat vague) evolutionary motivation are the minimization of overall ATP usage and the minimization of overall fluxes.

**Synthetic design-oriented objective functions** Metabolic modeling can also be used to identify metabolic states that lead to a certain desired behavior of a microorganism. For example, we may seek to genetically perturb a microbe such that it produces a certain compound of industrial or medicinal interest, while it also retains a certain minimal growth rate[200]. Indeed, it is often desired to retain a certain minimal ability to grow such that the genetically engineered organisms can be lab-grown after which the produced compound of interest can be harvested. In that case, we can combine maximizing the production rate of the compound while imposing an inhomogeneous constraint that sets a lower bound on the biomass production rate. This can even be combined with a calculation in which we solely maximize the biomass production rate: maximizing the biomass production rate is a model for the wild-type cell, whereas maximizing the generation of the compound models the desired phenotype. By comparing the flux distributions between these 'strains', we can search for target genes that should be up- or downregulated.

## 5.3. Optimal metabolism in terms of elementary flux modes

In the previous chapter we introduced elementary flux modes (EFMs) and identified them as the fundamental metabolic pathways that carry flux through the metabolic reaction network. Here, we will show how elementary flux modes also can be very useful for describing optimal metabolic states. We briefly recapitulate the notion of elementary flux modes. All metabolic flux vectors v that satisfy both the mass-balance and irreversibility constraints form a pointed polyhedral cone, called the flux cone. The EFMs are the extreme rays of this cone, so that they can be used to decompose all steady-state flux vectors:

$$v = \sum_{i} \lambda_{i} e_{i},$$

where $\lambda_i \ge 0$ and $e_i$ is the i-th EFM. Moreover, the EFMs turn out to be the minimal metabolic subnetworks that a cell can use in steady-state without needing any other reaction, so that we can view EFMs as minimal metabolic strategies.

In Figure 5.2 we depict the EFMs as black lines, and the region in-between these lines is the steady-state solution space
coupled with network compression.

Instead of inactivating reactions, one might be interested by computing all EFMs containing a specific reaction, such as the biomass, or several reactions, e.g. biomass synthesis and ATP maintenance. This is not a good idea to try to incorporate these constraints directly into the computation as such a constraint adds an hyperplane on the solution space, changing the resulting solutions [117, 118]. As a result, these kind of candidate constraints are best left for post-processing.

## 4.5 Alternative methods for flux space exploration

As we described above, exploration of all possible flux distributions using EFMs can become very complex for larger models. A genome-scale model, which comprises all metabolic reactions that an organism can catalyze, typically contains thousands of reactions, which prohibits the enumeration of EFMs. At the moment, it is unclear whether, even if we would have an enormously fast computer that could compute all EFMs, the number of EFMs would not be so large that we cannot store the EFMs anywhere, nor analyze it in any meaningful way. Here we discuss several alternatives for exploring the metabolic capabilities of a cell that try to avoid the combinatorial complexity that hinders EFM analysis.

### 4.5.1 Elementary conversion modes

If we are interested in the metabolic capabilities of an organism, is it always necessary to know all possible flux vectors? For example, what if we want to lab-culture an organism of which we have a reconstructed metabolic network, but no idea what nutrients it needs to grow. Then we only need to know from what combinations of nutrients it can make all its cell components. Or, what if we want to model the possible cross-feeding interactions between several microbial species? Then we are mostly interested in what each of them can consume and produce, and not really in how they do that. Elementary conversion modes (ECMs), introduced in 2005 by Urbanczik and Wagner [119], capture all possible overall conversions from nutrients to products that an organism can catalyze, while ignoring which individual reactions are used for this.

ECMs focus on the net results of metabolism, i.e. on the uptake and production of compounds external to the metabolic network, such as sugars, nitrogen sources, fermentation products but also 'biomass'. To get information about these compounds we need to extend our metabolic network by including the external compounds as rows in the stoichiometry matrix; this is in general easy to do since we already had exchange reactions (reactions where an external compound was imported or exported) so we only have to find the stoichiometric coefficient in which the external compound was involved in these reactions. Let us denote the original stoichiometry matrix by $N_{\text{int}}$ and the submatrix that we add by $N_{\text{ext}}$; together they form $N_{\text{tot}}$. We can then define the conversion cone:

$$C = \left\{ \frac{ds}{dt} = N_{\text{ext}}v \mid N_{\text{int}}v = 0, v \ge 0 \right\}. \qquad (4.24)$$

If we look carefully at this definition we can see that the flux vectors $v$ need to satisfy exactly the same constraints as in the flux cone (Eq. (4.2)). The only difference between flux and conversion cones is that we are either interested in the fluxes themselves, or rather in the conversions that they induce: $ds/dt = N_{\text{ext}}v$.

**Definition 4.** The set of ECMs is the minimal set of conversions $\{ecm^1, \dots, ecm^l\}$ (where $ecm_k^i$ is the amount of metabolite $k$ produced in the $i$th elementary conversion mode), such that

1. all conversions $ds/dt \in C$ can be written as a positive sum of these elementary conversion modes: $ds/dt = \sum_i \lambda_i ecm^i$, with $\lambda_i \ge 0$,

2. without the production of any metabolite being canceled in that sum, i.e. for all metabolites $k$ we either have for all $\lambda_i > 0$ that $ecm_k^i \ge 0$ or for all $\lambda_i > 0$ that $ecm_k^i \le 0$.

We will explain both parts of this definition below, but let us first remark that the definition is in fact perfectly analogous to the definition of EFMs: EFMs are the elementary vectors (or precisely: conformally
<figure><img src="image_40.png" /><figcaption>Figure 5.2: We show a cartoon of the solution space of a metabolic network, the so-called flux cone, with respectively zero, one and two constraints. With one constraint, the optimal solution for any linear objective can be attained in a vertex of the space, which means that it can be attained in a single EFM. With two constraints, we need to combine at most two EFMs to describe the optimal solution.</figcaption></figure>

that is spanned by the EFMs. Note that this illustration is great simplification, usually the flux cone is a high-dimensional object that can only be visualized in trivial toy examples. In fact, the flux cone is a subspace of $\mathbb{R}^n$ where $n$ (the number of reactions) can be in the thousands for a typical genome-scale metabolic network. Moreover, the number of extreme rays of the cone would be overwhelming, due to the complexity issues associated with EFM enumeration as described in the previous chapter.

Figure 5.2a also shows that there is a direction in which the objective increases fastest. This direction is determined by the choice of objective function, to be specific: the direction of maximal increase of the objective is given by the vector of coefficients, $[c_1 \cdots c_n]^T$, appearing in the linear objective function (5.3). However, as long as we do not impose an inhomogeneous constraint, the flux cone is unbounded, so that we can usually reach infinite values. This makes sense when we think of the metabolic states in terms of elementary flux modes: when we have an EFM that reaches some nonzero objective value, we can always multiply it by any positive scalar. This multiplication will increase the objective value, while the steady-state and irreversibility constraints will not be affected.

Metabolism, however, is never unconstrained, so we will always have at least one inhomogeneous constraint. In the previous chapter, inhomogeneous constraint were written in the general form

$$ \sum_i w_i^p v_i \le h_p, \quad p = 1, \dots, P \qquad (5.5) $$

where each $h_p$ corresponds to a component of the $P$-dimensional vector $\mathbf{h}$ and $n$ weights $w_i^p$ ($i = 1, \dots, n$) are supplied for each of the $P$ constraints. The second panel of Figure 5.2 shows how a single inhomogeneous constraint (i.e. the case $P = 1$) can constrain the flux cone and theory dictates an optimal flux vector is found at a vertex of the resulting flux polyhedron, which geometrically corresponds to the intersection of the flux cone and the hyperplane of the inhomogeneous constraint. One particular biological argument for such a constraint is related to resource allocation[201, 202]: only a limited number of macromolecules (proteins, ribosomes, etcetera) fit inside a cell. Since these molecules catalyze reactions, reaction rates are proportional to their concentrations:

$$ v_i = e_i g_i(\mathbf{x}_{\text{metab}}), \qquad (5.6) $$

where $e_i$ is the concentration of the enzyme that catalyzes reaction $i$, and $g_i(\mathbf{x}_{\text{metab}})$ is a function that describes enzyme kinetics in a non-linear way that is for most reactions unknown. The resource-allocation constraint then takes the form

$$ \sum_i w_i e_i \le 1, \qquad (5.7) $$

where $w_i$ are weights that determine how much of the resources are taken up by one unit of the $i$th enzyme; these
<figure><img src="image_53.png" /><figcaption>Figure 4.6: Elementary conversion modes – (A) Small toy network with three ECMs shown in blue, yellow and red. Note that the red mode can be decomposed as a positive combination of the blue and yellow elementary conversion modes, but that would cancel the production of B so this is not allowed. (B) The conversion cone is shown in gray, and the blue and yellow arrow correspond to the blue and yellow ECMs are the extreme rays. The red ECM needs to be added because it is on the intersection with the dB/dt = 0-plane.</figcaption></figure>

non-decomposable vectors) of the flux cone, and ECMs of the conversion cone. The reason that the def-
ition of ECMs has an additional requirement (2.) is just that the analogous requirement was automatically
satisfied for EFMs because we assumed all reactions to be irreversible.

In Figure 4.6A we show a small metabolic network with external metabolites *A*, *B* and *BM*, and internal
metabolites *C*, *D* and *E*. We can find 9 EFMs in this network: one that goes from *A* to *B*, four that produce
*BM* starting from *A* and four that produce *BM* from *B*. We get four EFMs to go from *A* to *BM* because
there are two ways of going from *C* to *D* and again two for converting *D* into *E*. This makes clear that
having a number of modules of alternative reactions can quickly give rise to large numbers of EFMs, even
though the overall conversion from nutrients to products remains the same. In contrast, we will explain
that we only get three ECMs.

In Figure 4.6B we see the conversion cone in gray. Note that this cone does not live in flux space, but rather in the space of external metabolite changes, or conversions. We recognize that the cone can be spanned by two extreme rays, which correspond to converting *A* into *B* (blue) and to using *2B* to produce *BM* (yellow), so these rays correspond to elementary conversion modes following the first part of Definition 4. Now why do we have a third ECM, when the blue and yellow one already span the whole conversion cone? Indeed, the third vector in Figure 4.6B can be obtained by summing the yellow vector and two times the blue vector: $2(-1, 1, 0) + (0, -2, 1) = (-2, 0, 1)$. However, note that the production of metabolite *B* would cancel in this sum, which is not allowed according to the second part of Definition 4. The reason that this second part of the definition is important, is that the elementary conversion modes are intended to capture all metabolic capabilities of an organism, so taking only the first two modes would not be enough: we also want to account for the possibility of making *BM* from *A* even if we decide that the elementary conversion mode from *B* to *BM* is not possible in the current environment, for example because *B* is not present as a nutrient in the medium.

Because many EFMs result in the same overall conversion, the exploration of metabolic capabilities can now be done in larger networks, at the cost of ignoring information about which reactions are used [120]. This way of thinking can be pushed even further: what if one is not interested in the conversions between all nutrients and products, but only between a subset of these? In that case, one would want to compute the ECMs only between the external metabolites of the most interest. This can be done with a small trick. Say that we are not interested in the production of external metabolite *X*. Before we start the enumeration algorithm we add a virtual reaction to the network that consumes and produces *X* from nothing, i.e. we add *X* ⇄ ∅, and then we change *X* from an external metabolite to an internal metabolite. Consequently,
weights can for example be proportional to the volume, the mass, or the number of amino acids of the enzyme. Making a change of variables to express the constraint in terms of fluxes gives:

$$ \sum_{i} \frac{w_i}{g_i} v_i \leq 1, \qquad (5.8) $$

such that these resource-allocation constraints again fit the form presented in Equation (5.5). A well-known example of a modeling framework that uses such a constraint is FBA with macromolecular crowding (FBAwMC, [203]) where such a constraint arises due to a physical limitation on the number of enzymes contained within the cell.

It is not necessarily always the case that an inhomogenous constraint applies to all EFMs. For example, in a metabolic model of an organism able to grow on multiple carbon sources, many EFMs may remain unbounded. For treatment of these cases, the reader is referred to [204]. Moreover, we may have multiple inhomogenous constraints on flux values as Equation (5.5) suggests. The third panel of Figure 5.2 illustrates how a second inhomogenous constraint can further constrain the solution space where theory implies an optimal flux vector is found on a vertex lying on the edge between two EFMs (as shown in the example in the figure). Imposing additional inhomogenous constraints can therefore lead to the superposition of additional EFMs in the solution. In general, if we consider a constraint-based model with K inhomogeneous constraints it can be proved that an optimal flux vector will be built out of at most K EFMs[204]. We therefore see another important property of EFMs: not only do they form the minimal building blocks that span all metabolic capabilities of the cell, they are also optimal building blocks. When metabolism is optimized, only few of these EFMs are used. As a result, solutions to linear constraint-based optimizations can usually be rationalized in terms of the properties of the available EFMs[205], for example, a flux balance analysis with only one constraint on a nutrient uptake will just return the EFM with the highest 'yield', i.e. the highest efficiency of making biomass per nutrient.

## 5.4. Phenotypic phase plane analysis

The analysis of the metabolic response to environmental changes is often sought assuming that there is only one substrate limiting growth (or other metabolic reaction). For example, we could be interested in the growth and ethanol production by S. cerevisiae under oxygen limitation in a chemostat. In this experimental setup, every other substrate should be provided in excess, including the carbon and energy source. If no oxygen is supplied, ATP must be produced only using oxidative phosphorylation reactions and a fermentation product, such as ethanol, will be produced. On the other extreme, if enough oxygen is available, a fraction of the carbon source will be completely oxidized, producing ATP via respiration. In both cases, the fraction of the carbon and energy source not used for energy generation will be used for the production of biomass at a specific growth rate equal to the dilution rate of the chemostat.

This behavior can be analyzed using the phenotypic phase plane analysis. To calculate a phenotypic phase plane (PPhP), the uptake fluxes values under analysis, typically the uptakes of oxygen and the carbon source are discretized between their upper and lower values and used to construct a meshgrid containing the 2-D grid coordinates based on the coordinates contained in the discretized vectors of oxygen and carbon uptake fluxes. At each tuple in the 2-D grid a FBA problem is solved after fixing the lower and upper bounds of the corresponding fluxes to the values in the tuple. Figure 5.3.A shows the PPhP of the metabolic network presented in Figure 5.1 with $s_X = 10$, $s_1 = 1$, $s_3 = 4$, $v_4 = 0$, $v_0^{UB} = 10$, and $v_{O_2}^{UB} = 15$ mmolg$^{-1}$h$^{-1}$.

At zero oxygen uptake, Figure 5.3.A shows that growth is possible reaching a specific growth rate of $1 h^{-1}$ at the maximum glucose uptake. Notice that the slope of the line connecting the origin of coordinates and the point of the highest growth rate at a glucose uptake of $10$ mmolg$^{-1}$h$^{-1}$ is $0.1$ gCDW(mmol)$^{-1}$. Biologically, this slope corresponds to the biomass yield on glucose under anaerobic conditions, and in terms of linear programming to the negative of the shadow price defined as:

$$ \gamma_i = \frac{-dz}{db_i^v}, \qquad (5.9) $$

where z is the objective function optimal value (specific growth rate in this case) and $b_i^v$ corresponds to the violation of
it now has to satisfy the mass-balance constraint (which can always be done trivially using the added virtual reaction), and will thus never show up in the computed elementary conversions. In this way it was possible to compute all ECMs between glucose, oxygen and biomass for a real genome-scale network of *E. coli*.

### 4.5.2 Flux sampling

In addition to the computational complexity of EFM enumeration for large metabolic networks, these objects are not necessarily related to experimentally-derived flux measurements. This is because when a vector of experimentally-measured flux values $v$ would be decomposed into EFMs, this generally does not give a unique solutions because it can be done in many ways. Flux sampling methods can be employed to solve both the computational and interpretability problems simultaneously, exploring the set of flux vectors (i.e. directly measurable in principle) by computationally sampling from the flux space. The goal of flux sampling in general terms is to produce a sequence of flux vectors that satisfy the steady state constraints until enough samples have been generated to provide an approximate representation of the entire flux space. The flux polyhedra defined by mass-balance and additional inhomogenous linear constraints are convex, and therefore uniform sampling of these flux spaces can be achieved using variants of an algorithm developed for convex analysis called the coordinate hit-and-run (CHR) algorithm [121]. Briefly, the most basic implementation of the CHR algorithm generates a Markov chain of flux vectors by starting in a random position within the flux polytope, picking a direction at random (uniform), and moving a random distance (uniform) in that direction from the current point. The resulting point is returned as a flux vector instance and the process repeats from there. It has been proven that the CHR algorithm converges to a stationary distribution of the Markov chain that is a uniform distribution in the flux space. Alternatives to uniform sampling (i.e. alternative distributions across the flux polytope) can also be achieved using variants of the CHR algorithm.

As highlighted previously in Section 4.4.3, mass-balance and inhomogeneous linear constraints alone often do not contain enough information to sufficiently reduce the space of biologically-feasible flux vectors. For example, thermodynamic constraints on flux vectors are important for ruling out a large proportion of the sampled flux vectors as infeasible, but this may disproportionately dominate the resulting sampling distributions. Unfortunately, for mathematical reasons too deep to go into here, simply removing these infeasible flux distributions post-sampling will not result in a uniform distribution over the thermodynamically-feasible portion of flux space. In fact, this relevant subset of flux space cannot be defined explicitly, and is usually neither convex nor connected meaning that no Markov chain methods exist for sampling. As an alternative, a recent method [122] has been developed to combine thermodynamic constraints, physiological observations and estimated thermodynamic parameters, with mass-balance and inhomogeneous linear constraints to provide a probabilistic thermodynamic analysis of metabolic reaction networks. Advances such as these will almost certainly aid a more complete characterization of flux space as data and methods become available.

### 4.5.3 Minimal cut sets

A minimal cut set (MCS) is a set of reactions that, when disabled, disables a set of modes, which in turn can represent a biological function, such as the secretion of a side product. This enables the prediction of gene deletion targets, given that the genes coding for the involved reactions are known. A cut set is minimal if the removal of one or more reactions from the set leads to at least one of the targeted modes not being disabled.

In order to avoid also disabling desired functionalities, such as product secretion and growth, the concept of constrained minimal cut sets (cMCSs) has been developed. cMCSs enable targeting a set of modes while at the same time making sure that some elements of another set of modes will remain active.

**Motivation for (constrained) Minimal Cut Sets** The concept of MCSs was introduced by Klamt and Gilles in 2004 [123] and subsequently generalized and improved [124, 125, 126]. As briefly outlined above, the idea is to define a set of EFMs which should be disabled, for example because they generate an unwanted side product or because they don't generate the product of interest with a sufficiently high yield. Since EFMs are minimal, removing a single reaction will disable it. A cut set is a set of reactions of which at
<figure><img src="image_41.png" /><figcaption>Figure 5.3: Phenotypic phase plane of the metabolic model shown in Figure 5.1, calculated as a function of the uptakes of oxygen and glucose</figcaption></figure>

a mass balance constraint and is equivalent to the uptake reaction of the *i*-th metabolite (glucose in this example)[206]. Figure 5.3.C shows that the glucose shadow price is equal to -0.1 at every point in the feasible region of the problem. Figure 5.3.D shows the shadow price values for oxygen uptake. For every unit increase in the oxygen uptake flux, the biomass specific growth rate increases by 0.4 $h^{-1}$. Thus, the plane of increasing growth rate values in Figure 5.3.A can be described by the equation $0.1v_G + 0.4v_{O_2}$. Concomitantly, as the oxygen uptake increases, the flux of product P1 decreases as more ATP is generated in reaction $v_3$. For every constant glucose uptake flux, the specific growth rate increases and the production of P1 decreases until the optimally line (red line) is reached in Figure 5.3.A. This line represents the optimal relation between the two metabolic fluxes in the PhPP [206]. In this example, the optimally line represents the combinations of glucose and oxygen uptake fluxes leading to a complete oxidation of the substrate, and thus supporting the maximal biomass yield. Finally, increasing the oxygen consumption beyond the optimally line, at a constant glucose uptake, leads to an infeasible problem since there is no further glucose to be oxidized.

## 5.5. Non-uniqueness of the optimal metabolic state

Although the optimization of some objective function strongly reduces the number of solutions, it is still possible that many different metabolic states satisfy the constraints and reach the same maximal value for the objective. In that case, we are again undecided on which of the solutions gives the most useful information about the biological problem. This non-uniqueness of the optimum can be explained in terms of the elementary flux modes. In the second panel of
least one is active in each of the EFMs in the targeted group. Thus, disabling the reactions contained in the cut set will disable all of the targeted EFMs, and each cut set therefore represents the prediction of a set of gene deletions. Since it would be pointless to remove reactions which only target EFMs that were already targeted by other reactions, cut sets are required to be minimal. This means that removing a single reaction from the cut set would lead to one or more of the targeted EFMs to survive the intervention and also that adding a single reaction to the cut set would have no additional effect on the set of target EFMs.

The pitfall when using MCSs is that while they guarantee the elimination of the targeted EFMs, all other EFMs may be affected as well. This means that modes with desired phenotypes, such as high growth and/or high product yield, may become impossible. Therefore, cMCSs were developed [127]. In this extension of the concept of MCSs it is now possible to additionally define a set of EFMs which are desired, i.e. which can not be disabled by the cMCSs. This is usually implemented by the requirement that at least a specified minimum number of EFMs of the desired set need to remain active. Summarizing, cMCSs are sets of reactions which guarantee that (i) the full set of target EFMs is disabled and (ii) a certain minimum of desired EFMs has to remain unaffected. The drawback, with both MCSs and cMCSs, is that the target (and desired) EFMs need to be defined. This is generally achieved by defining cut-offs in terms of product yield and growth, which is, however, ultimately arbitrary.

**Calculation of (constrained) Minimal Cut Sets** Since minimal cut sets in a metabolic network are EFMs in a dual network [128], methods used for calculating EFMs can be used to calculate MCSs. Among other approaches [129] one based on binary integer programming has been developed [130, 131]. While it requires that the EFMs are calculated before it can be applied, the advantage is that the algorithm is very intuitive. After having calculated the modes, each is represented as a binary vector which is zero for reactions with zero flux and one otherwise. The EFMs are then divided into either targeted or desired. A binary vector, corresponding to the cMCSs being calculated is introduced. It will have a one if the corresponding reaction remains active and zero if the reaction is disabled. The first requirement is that cMCS needs to disable all target modes and thus the vector must have zero elements such that each target EFM must have at least one corresponding non-zero element. The second requirement is that at least a defined minimum of desired modes must remain active. This is achieved by introducing a second binary vector. This vector has an element for each EFM and is calculated so that it has a zero when the mode is disabled by the cMCS and one otherwise. By adding the constraint that the number of ones in this vector must at least equal the previously defined minimum, the second requirement is met. Maximizing the vector corresponding to the cMCS yields the first solution. The next solution can be found by adding constraints to make sure that the current one is excluded.

## 4.6 Concluding remarks

In this chapter we studied how the individual reactions that an organism can catalyze together give rise to the overall conversion of nutrients into cell components and secretion products. For that, we studied the cell's metabolism under a number of simplifying assumptions, most notably, we model metabolism in steady-state. Given this steady-state constraint, we explained how all feasible flux distributions form a space of a specific type: a pointed polyhedral cone. By exploring this 'flux cone' we can chart the metabolic capabilities of an organism.

We have seen that an exhaustive charting of these metabolic capabilities is the computation of all *elementary flux modes*: minimal subnetworks that can individually give rise to steady-state flux distributions, and that may be interpreted as minimal metabolic strategies. An especially important use of EFM analysis can be found in the prediction of the effect of gene knockouts: when all EFMs that produce compound *Y* use reaction *r*, then the organism cannot make this compound when the gene is knocked out that codes for the enzyme that catalyzes *r*. And conversely, sometimes gene knockouts can be found such that the cell cannot grow anymore without producing a certain compound of interest. Clearly, these analyses can be very useful for the design of organisms in bio-industry.

On the other hand, we also saw that for large models the computation of all EFMs becomes impossible. There are simply too many of these minimal subnetworks. We presented several alternatives. One could use *elementary conversion modes* if one still desires an exhaustive list of the metabolic capabilities of the
Figure 5.2 we saw that the optimal solution was located on the vertex that was as far as possible in the optimization direction. One can imagine, however, that the flux cone can be located in the space such that there are two vertices that both reach out equally far into that direction. In that case, the two corresponding elementary flux modes perform equally well, and consequently, all convex combinations of these elementary flux modes also reach the same objective value. In metabolic modeling we often work in a high-dimensional space with constraints that concern only few of those dimensions (for example a bound only on a nutrient uptake rate). In such cases it is very likely that many elementary flux modes perform equally well, so that there is a whole subspace of equivalent solutions.

### 5.5.1. Flux Variability Analysis

The equality and inequality constraints of the FBA problem form a polytope where the problem is feasible, a cone if the problem is written in canonical form. The optimal solutions of the LP problem can lay on a vertex of the polytope, and be unique, or be non-unique solutions if the objective function hyperplane is parallel to a facet of the constraint polytope at the solution. This means that one or several variables can change their values without affecting the value of the objective function. These variables can be identified using flux variability analysis (FVA), where each flux of the reactions in the metabolic network (the set of *J* reactions with *N* elements and *I* metabolites) maximized and minimized, one at a time, while fixing the value of the objective function to a fraction of the optimal value obtained in the original FBA problem.

$$ 
\begin{array}{l}
\max_{v} v_j (\text{or } \min_{v} v_j), \quad \text{such that:} \\
\displaystyle\sum_{j \in J} S_{i,j} v_j = 0, \quad \forall i \in I, \\
LB_j \le v_j \le UB_j, \\
v_{\text{biomass}} = f \cdot v_{\text{biomass}}^*, \\
v_j \in \mathbb{R}, \quad \forall j \in J.
\end{array} 
\qquad (5.10) 
$$

Hence, 2N optimization problems need to be solved if there are N unconstrained fluxes. The results of the VFA analysis should be carefully interpreted. Since the maximum and minimum fluxes are calculated one at a time, and although changes in this flux might not affect the objective function, this typically requires changes in the remaining fluxes. Therefore, the polytope that describes all alternate optimal solutions is not captured by VFA. Instead, FVA inscribes this polytope in the smallest possible “box” [207]. Besides being useful for the identification of alternative solutions, FVA can be utilized to identify blocked reactions under a given growth condition. These reactions are characterized by minimum and maximum flux values (as calculated by VFA) equal to zero and arise due to regulatory constraints imposed to the FBA or due to network gaps, for example, metabolites lacking a consumption or production pathways for whom a steady-state mass balance is impossible. Thus VFA, could help in the identification of dead-end metabolites, and in the long run, in model improvement.

## 5.6. Limitations of constraint-based metabolic models

In this and the previous chapter, we have introduced constraint-based analysis of genome-scale metabolic models. We started by pointing out many of the simplifying assumptions that are associated with the study of large metabolic reaction networks. For example, we only considered systems in chemical steady state with their environment, we ignored the effects of metabolite dilution, and we made semi-informed choices for which intracellular molecules are contained in our model or summarized in a biomass reaction. All these assumptions can be relaxed, at the cost of making models more complex. Although it is tempting to think that the more complex a model the more realistic it will be, there is not much use to adding additional complexity if we don't have the data to support it. Constraint-based analysis therefore
cell. The ECMs are easier to enumerate because one can choose to focus only on all possible conversions between (a subset of) the nutrients and products, instead of requiring all information about which reactions are used to get these conversions. For the design of gene knockouts specifically, *minimal cut sets* may be used. Finally, we discussed that the flux cone can be sampled randomly to characterize the flux cone, if this characterization does not need to be exhaustive.

In many cases we have additional information that determines that part of the flux cone is infeasible. For example, some metabolic fluxes may have been measured so that these reaction rates can be fixed to their observed value. In other cases, one may want to use thermodynamic properties to prohibit reactions from occurring that would violate the second law of thermodynamics. These additional constraints can be imposed on top of the mass-balance constraint to further bound the space of feasible flux distributions; each correctly-imposed constraint narrows down the space of feasible fluxes, and thus increases our knowledge of the metabolic state of the cell.

All explorations of the space of feasible flux distributions show one unavoidable conclusion: the metabolic network is incredibly flexible. Even when several constraints are imposed, a genome-scale metabolic model will allow for an almost incomprehensible number of modes in which the metabolic network can function. Consequently, to predict the metabolic state of a cell in more detail we need to make an additional assumption. In the following chapter, we will study what predictions we can make when we assume that the metabolic state is optimized to perform a certain function.

## Recommended readings

**Elementary flux modes** and their applications are introduced in an intuitive way in: J. Zanghellini, D. E. Ruckerbauer, M. Hanscho, C. Jungreuthmayer (2013). *Elementary flux modes in a nutshell: Properties, calculation and applications*. Biotechnology Journal 8 (9), 1009. doi: 10.1002/biot.201200269

**Elementary Flux Vectors** were introduced as an analog of Elementary Flux Modes in the case that the flux mode is further bound by at least one inhomogeneous constraint. A nice review of these EFVs is can be found in: S. Klamt, G. Regensburger, M. P. Gerstl, C. Jungreuthmayer, S. Schuster, R. Mahadevan, J. Zanghellini, and S. Müller (2017). *From elementary flux modes to elementary flux vectors: Metabolic pathway analysis with arbitrary linear flux constraints*. PLoS Computational Biology, 13(4):e1005409, doi: 10.1371/journal.pcbi.1005409.

## Problems

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

### Problem 4.1 A small metabolic network (1)

*Spirallus insilicus*, a completely fictional organism [133], is characterized by the metabolic network depicted in Figure 4.7 *X*, *S* and *P* represent the biomass, one substrate and one product, while metabolites *A* to *E* denote intracellular metabolites. One directional arrows indicate irreversible reactions (all but $v_4$)

(a) How many intracellular metabolites, intracellular reactions and transport reactions are involved in the model?

(b) Obtain the stoichiometric matrix ($N$) and the vector of fluxes. How many elements are in the product $Nv$ and what do they represent?

(c) Is the matrix $N$ of full rank? How many fluxes should be specified to have a unique solution?

(d) Transform the set of constraints so that they define a pointed cone. Determine the number of variables (fluxes) and constraints.

### Problem 4.2 A small metabolic network (2)
provides one way to study metabolism at genome-scale when data are limited. In the following chapters we will study
the consequences of lifting one or more of these simplifying assumptions.

In constraint-based analysis, one considers reaction rates (fluxes) as the variables in the model, giving the illusion that these are directly set by the cell to regulate its metabolic state. In reality, however, the reaction fluxes are the combined consequences of enzyme expression, regulation and metabolite concentrations. If we wish to model metabolism in more detail, we should build models that incorporate gene expression and metabolite concentrations systematically. Some of the next chapters attempt this, but we have described that FBA is useful when experimental data are limited. Certain extensions of FBA discussed in later chapters also move beyond the steady state assumption, allowing the environment to change with time. One example is the method dynamic FBA, which will also be discussed in a later chapter.

Philosophical remarks 5.B: Qualities of a model

When have we made a good model? Is the quality of a model determined by whether it fits all experimental
observations? What is the ideal size of a model? Is the purpose of a model that it predicts, or rather that it
provides insight into the biological processes?

The answers to these questions are as common as it is unsatisfying: 'it depends'. Sometimes a model can be
very useful if it just predicts, and does not explain, as witnessed by the undebatable success that machine learning
models have across the sciences. However, only true understanding of the studied process can lead to hypotheses
and predictions on phenomena that are far away from the currently available data. The more a model is fitted to
a specific dataset, the less we are able to extrapolate it beyond this dataset.

These questions are very relevant in the context of metabolic modeling. Metabolic models have many unknown
parameters, stemming from our ignorance of the biological process: What is the true objective? What constraints
are relevant for determining metabolism? It is a deceptive trap to view the success of the model in reproducing
the observed data as a validation that the right parameters, objective and constraints were chosen. A successful
model only indicates that the modeled mechanism can be similar to the true biological mechanism, but it does not
show it is. The problem is that, since we have many different parameters to choose from, many different models
can explain the same metabolic observations[208].

An especially important question is whether metabolism is truly optimized for some evolutionary function. It is now
an attractive option to view the success of optimization-based models as proof that the cells are indeed optimized,
but this would be wrong because we can also explain the data with models that do not require optimization.
To really quantify whether metabolism is optimized we should therefore devise quantitative tests that distinguish
between randomly chosen and optimized metabolic states. An interesting approach for describing the metabolic
outcome of cells, relying on statistical mechanics rather than on a selected objective function, has already been
introduced [209].

5.7. Conclusions of our flux cone optimizations

In this chapter we built upon the exploration of flux spaces derived from constraints by imposing optimality criteria in terms of an objective function. The choice of the objective function(s) and the constraints depend on the modeling purpose. We will summarize some of the possible choices by listing three purposes that this type of models can have.

First, constraint-based optimization can be used to collect, integrate and extrapolate data on the metabolism of a
specific organism. In this case, as much experimental information as possible can be used to refine the model. For
example, measured fluxes can be fixed with constraints, measured metabolite concentrations can be used to determine
the thermodynamically feasible direction of reactions, and transcriptome information can be used to exclude some
reactions because the corresponding genes are not expressed. One of the applications is then that unknown variables
can be inferred such that they are in accordance with the metabolic network and all the measured variables.
<figure><img src="image_54.png" /><figcaption>Figure 4.7: *Spirallus insilicus* network, adapted from [132]</figcaption></figure>

Consider the following small metabolic network:

$$
\begin{align*}
S_e &\xrightarrow{v_0} S_c \\
S_c &\xrightarrow{v_1} P_c \\
P_c &\xrightarrow{v_2} C_c \\
P_c &\xrightarrow{v_3} D_c \\
P_c + 2C_c &\xrightarrow{v_4} X
\end{align*}
$$

Metabolites with a *c* subscript are located in the cytosol (intracellular) while *e* stands for extracellular and *X* represent biomass. All fluxes are positive.

(a) Represent the model as a reaction network (a sketch with metabolites and reactions)

(b) Obtain the stoichiometric matrix (N) and list the variables of the metabolic model (v)

(c) Show that there is no solution to the mass balance equation N v = 0 producing metabolite D.
Identify why this is so and modify the model so the production of D is allowed ($v_3 > 0$)

**Problem 4.3 Elementary Flux Modes (1)**

Assume reaction *v*₄ is irreversible from *A* to *D* in *Spirallus insilicus* (Problem 4.1). Calculate all the Elementary Flux Modes.

(a) By hand.

(b) Using a software of your choice (e.g. pypi.org/project/efITOOL/)

**Problem 4.4 Elementary Flux Modes (2)**

Consider the following metabolic network

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

$$
N = \begin{pmatrix}
1 & -2 & 0 & 0 & 0 & 0 & -2 \\
0 & 1 & -2 & 1 & 0 & 0 & 0 \\
0 & 2 & 0 & -1 & -2 & 0 & 0 \\
0 & 0 & 1 & 0 & 1 & -1 & 1
\end{pmatrix}
$$
Second, hypotheses can be tested on why the studied organism attains its metabolic state. By choosing an objective function we can propose what *drives* the metabolic behavior and by choosing the constraints we propose what *limits* the metabolic behavior. If the model is then in accordance with the experimental observations, we know that at least the hypotheses were not proven wrong. On the other hand, we must be careful not to conclude from this that the hypotheses must be right, as we discussed in the box with philosophical remarks.

Third, we may use these models to search for a metabolic state that results in a certain desired behavior, for example in the secretion of a product that is useful for industrial or medicinal reasons. In this case, the objective function is picked such that exactly the desired behavior is maximized, often while requiring that some biomass production is still possible because the cells need to be able to grow before the harvesting of the product can start.

Despite these useful purposes, we have also identified several limitations of the FBA-type models that we described here, such as ignoring metabolite concentrations, enzyme kinetics, and the assumption of a stationary metabolic state. The reason that these models are still very popular is their computational simplicity: as long as the objective function and constraints are linear in the reaction rates, the optimal solution is relatively easy to find using linear programming. This makes it feasible to make and run these models on genome-scale metabolic networks, which are networks that comprise all the metabolic enzymes for which the genome encodes, and can include thousands of reactions.

Understanding the solutions of such large models can also be very difficult due to their dimensionality. This is made easier when one uses elementary flux modes: we have seen that a solution is always a combination of a relatively small number of EFMs. More precisely, the number of EFMs that are active in the optimal solution cannot exceed the number of imposed constraints. This means that to understand the solution, we only need to understand which EFMs are selected and why. As such, we can interpret optimal solutions in terms of the EFMs, i.e. the minimal metabolic strategies, that are used.

Finally, we explained that the optimal solution is often still not unique, and in a future version of this book, we will describe several methods that allow us to explore these equivalent optimal solutions.

## Recommended readings and tools

* Escher FBA (https://sbrg.github.io/escher-fba/) is nicely illustrative. It does FBA on an *E. coli* core model. Bounds on all reactions can be changed and different objectives can be explored. The resulting flux distribution is shown graphically.

## Problems

**Problem 5.1** Augment the metabolic network of *Spirallus insilicus* (Problem 4.1) by adding the in-homogeneous constraint $v_{upt} \le 10 \frac{\text{mmol}}{\text{gDWh}}$ and calculate the flux distribution if biomass is the objective function (maximize $v_5$).

1. Using a spreadsheet and its associated linear programming optimizer.

2. Using an LP solver in Python such as `linprog` available in `scipy.optimization`.

3. Is the flux distribution unique? Calculate the maximum and minimum values of each flux (except for the uptake of substrate and biomass production) if $v_5$ should be equal to its optimal value ($v_5^*$) and if this constraint is relaxed to $v_5 \ge 0.9 v_5^*$.

**Problem 5.2** The metabolic network illustrated in Figure 5.1, adapted from [206], was designed to include four phenotypes that can be reached depending on the ratios of the oxygen and carbon source (A) uptake, defining zones of single nutrient and dual nutrient limitation.

1. If the uptake of the carbon source A is bounded between 0 and $10 \frac{\text{mmol}}{\text{gDWh}}$ and no restrictions on the oxygen
Please note that some stoichiometric coefficients in *N* are different from 1 (not shown in the graphics).

(a) In the network drawing, gray dots denote carbon atoms. Check that carbon atoms are conserved in all reactions. What's the carbon content of the byproduct (not shown) of the reaction from A to D?

(b) All metabolites are treated as internal, that is, they need to be mass-balanced. Find all EFMs (by pure reasoning or by using a software). Determine all EFMs in which all fluxes are in forward direction, i.e. along the "conventional directions" indicated by arrows.

(c) Which of the EFMs are thermodynamically realizable? Explain why.
uptake are imposed, prepare a plot showing the biomass, C, D and E fluxes attained at different uptakes of A.

2. Repeat the preceding analysis, but limit the maximum uptake rate of oxygen to 10 $\frac{\text{mmol}}{\text{gDW h}}$.

3. If substrates uptakes are bounded between 0 and 10 $\frac{\text{mmol}}{\text{gDW h}}$ for A and 0 and 20 $\frac{\text{mmol}}{\text{gDW h}}$ for oxygen, calculate the phenotype phase plane. In each region of the phase plane (defined by a different slope), pick a combination of A and oxygen uptakes and analyze the fluxes of C, D and E.

<table><caption>Table 5.1: Stoichiometry of the metabolic network for problem 5.2. Adapted from [210] after [206].</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>$\frac{q_A}{\rightarrow} A$</td><td>$\frac{ATP}{\rightarrow} R_{ft}$</td><td></td></tr><tr><td>$A + ATP \xrightarrow{R_1} B$</td><td>$C + 10ATP \xrightarrow{R_z} \text{Biomass}$</td><td></td></tr><tr><td>$B \xrightarrow{R_2} 2ATP + 3NADH + C$</td><td>$\frac{q_{O_2}}{\rightarrow} O_2$</td><td></td></tr><tr><td>$0.2C \xrightarrow{R_3} 2NADH$</td><td>$C \xrightarrow{C_{out}}$</td><td></td></tr><tr><td>$C \xrightarrow{R_4} ATP + 3D$</td><td>$D \xrightarrow{D_{out}}$</td><td></td></tr><tr><td>$C + 2NADH \xrightarrow{R_5} 3E$</td><td>$E \xrightarrow{E_{out}}$</td><td></td></tr><tr><td>$NADH + O_2 \xrightarrow{R_{Res}} 2ATP$</td><td></td><td></td></tr></tbody></table>
# Chapter 6

## Rating metabolic pathways by enzyme efficiency

Wolfram Liebermeister and Elad Noor

### Chapter highlights

* In this chapter we discuss why certain pathway designs have been selected by evolution, by hypothesizing that some are more beneficial than others - based on several possible criteria and optimization goals: minimizing the number of reactions, maximizing product yield, increasing reaction turnover rates, and avoiding small thermodynamic driving forces.

* It turns out that all these criteria are related to a single objective: minimizing enzyme demand per product production rate or, equivalently, maximizing "enzyme productivity".

* We first focus on simple unbranched pathways with predefined flux distributions. We discuss several feasibility and optimality problems where metabolite concentrations are independent variables and solve for the minimal enzyme demand. In this setting, we see how enzyme productivity can be assessed or predicted and how it depends on different system parameters such as kinetics, thermodynamics, and concentrations of enzymes and metabolites.

* We discuss the difference between growth rate and yield. We then illustrate it by comparing between pathway options for glycolysis.

### 6.1. What guides evolution to select one pathway over another?

In the previous chapters, we asked what flux distributions are possible in a network, and which are most profitable for a certain task. Now we shall ask, more specifically, what led to the choice of existing pathways, or what makes a pathway variant favorable over another one that exists, or may have existed, in evolution. Of course, the same question plays also an important role in metabolic engineering, when new pathways are added to an organism, typically with the goal of achieving a maximal production, while imposing the smallest possible burden on the cell.

The chemical space is vast and many options exist for the same process, even if we consider only reactions with known enzyme mechanisms and impose thermodynamic constraints. Hence, while evolution had a choice between many pathway variants, only a tiny fraction of these possible variants is actually realized in nature, and a core part of central metabolism almost always follows the exact same design. The few exceptions that exist actually prove the rule, such as the two natural variants of glycolysis discussed later in this chapter. How can we understand why a certain variant is used in a certain organism or situation? And why are many variants not used at all? Moreover, some very successful pathways show
Chapter 5

Optimization of metabolic fluxes

Hester Chapman, Jan Pilipp Dapprich, Daan de Groot, Stefan Müller, Felipe Scott, and David Tourigny

Chapter overview

○ An optimization objective can be added to constraint-based models to make more specific predictions.

○ Different purposes can be served by choosing different optimization objectives and constraints

○ The optimal solutions can be understood in terms of elementary flux modes

5.1 Can optimality principles help us predict metabolic behavior?

In the previous chapter, we characterized an organism's metabolism by listing all the biochemical reactions that can be catalyzed by the enzymes encoded within the organism's genome. To understand how the genome constrains patterns of metabolic flux we needed to make several simplifying assumptions. The first important assumption was that intracellular metabolism is at steady-state, i.e., that the production and consumption of all metabolites is balanced such that their concentrations are constant in time. These resulted in the mass-balance constraints on the flux vector **v**. The flux cone of all flux vectors satisfying the mass-balance constraints could be further reduced by additional constraints on **v**, based on extra physical and biological assumptions about the magnitude and directionality of certain reactions within the network. We introduced several ways in which the entire flux space could be explored.

When applied to very large metabolic networks, the flux space will often contain an infinite number of flux
vectors **v** that simultaneously satisfy all constraints. From a mathematical perspective, this implies that
the constraints do not include enough information to uniquely specify a flux vector **v**. This makes sense
biologically, since if we imagine constraints are related to experimental observations it is very unlikely that
we will ever be able to make enough to fully account for every reaction encoded within the entire genome
of an organism (no matter how simple it might be). Often, however, researchers do want to further narrow
down the set of flux vectors that they think biologically relevant to the organism and conditions they are
studying, perhaps even to a unique **v** imagined to describe the metabolic state of an organism at a given
moment in time. One popular approach for doing so is to provide an additional assumption (or set thereof)
in the form of an objective function: it is assumed that the metabolic state of an organism is such that some
function of **v** (e.g. growth rate) is maximized to satisfy some criteria (e.g. evolutionary selective advantage).
The computational problem then becomes one of constrained-optimization: find a flux vector **v** that is
optimal in terms of the objective function(s) that simultaneously satisfies all constraints. The resulting
space of optimal flux vectors (sometimes containing just one unique vector) is often considerably smaller
than the space of those that satisfy only the constraints.

In this chapter, we will study metabolic models based on constrained-optimization. We will introduce a
features that might appear strange at first glance [211]: in glycolysis, an initial investment of ATP is required, and only
later it is recovered in higher amounts leading to a net gain. Is this just an evolutionary accident, i.e. a case where the
pathway that evolved first is the one that stuck around although it is not necessarily better than all the alternatives? Or,
rather, evolution did manage to find the optimal solution and therefore we should try to explain what the advantages of
these "engineered" features are?

In this chapter, we assume that it was a selection for functional features, not chance, that determined these pathway "choices", and ask: what guides evolution to select one pathway over another? What are the criteria that make pathways "efficient" or "profitable" for a cell or, alternatively, for a metabolic engineer? To compare pathways, we assume that each pathway comes with a predefined flux distribution, and therefore a predefined product yield, and alternative pathways (yielding the same product) are compared at equal product production rates.

When people talk about natural ecosystems, diversity is usually the first topic discussed. Indeed, evolution through
natural selection is almost guaranteed to create diversity where species evolve to occupy biological niches while exploring
the vast space of possible phenotypes. Similarly, the world of biochemistry is a vast space of possible reactions. Metabolic
enzymes participate in a network of pathways that supply cells with energy, and building blocks for biomass. Scientists
have been studying these biochemical reactions for nearly 300 years [212] – so far tens of thousands such reactions have
been classified; certainly many more exist in nature. Here are a few online databases where biochemical reaction data
are collected or predicted: MetaNetX, KEGG, MetaCyc, BiGG, ModelSEED, ATLAS of biochemistry.

To study the choice between pathways variants, we consider alternative pathways leading from A to B (or having a certain net sum formula) and their respective advantages and disadvantages. For simplicity, let us focus on biosynthesis pathways whose main task is more or less clear: producing a precursor molecule. Thus, the theoretical question would be: if a cell needs to make B from A, which pathway should it use? More specifically, how should the metabolic reactions be chosen and in what order? What should their kinetics and how should they be regulated?

If the pathway variant found in nature is due to selection for “good functioning”, then what are the features that make existing pathway designs successful? In short, what are criteria for “good” pathways? One possible criterion seems to be simplicity, that is, choosing a short route from pathway substrate and pathway product.

Philosophical remarks 6.A: What do we mean by a pathway?

The notion of “pathways” is common in cell biology to describe a set of reactions, proteins, or processes that form
a functional unit. However, there is no general definition: in practice, a pathway is often just a subregion of interest
within a larger network. In metabolism, “pathways” often lead from some important substrate to some important
product, with a simple and predefined flux distribution that consumes substrate(s), generates product(s), and may
or may not make use of co-factors. Considering fluxes in specific pathways (instead of flux distributions in the
entire network) is often a practical choice and, importantly, a choice that assumes that we can model, understand,
manipulate, or engineer such a pathway without strongly affecting the rest of the cell. This has a number of
benefits: (i) Instead of studying a huge network, we can look at pathways separately; (ii) there are reasons to
believe that the flux distributions in enzyme-efficient metabolic states must be elementary flux modes (see Chapter
4). Since EFMs often entail discrete choices between different pathways, it can make sense to study these pathways
separately (iii) once we understand the costs and benefits of single pathways (with a single, scalable flux mode),
we can apply the same thinking to analysing flux distribution on the entire metabolic network. Thus, in the rest
of this chapter, all results about “pathways” will also hold generally for entire networks, as long as a (scalable) flux
mode is given. Instead of comparing alternative pathways, we can compare alternative flux modes. In the following
chapter, we use this for optimizing over the set of all possible flux modes that a given network can support.

In contrast to the huge diversity that is allowed by the catalytic capabilities of enzymes, a few metabolic pathways are
extremely ubiquitous and exist virtually in every living cell. For example, glycolysis is a general term for pathways that
convert glucose to pyruvate while producing ATP [211]. One variant of glycolysis, named after Gustav Embden, Otto
selection of commonly used objective functions and the computational methods used to solve the associated constrained-optimization problem. We will also characterize optimal solutions that we get in terms of the minimal metabolic strategies that we identified in the previous chapter: elementary flux modes. Finally, we will explain how we can handle the cases where the solutions are, even after optimization, not unique.

## 5.2 Metabolic modeling based on linear optimization

We begin by introducing linear programs and transforming them into standard form. (For the history of linear programming, see Box 5.A.) Then, we discuss two types of inhomogeneous constraints in Flux Balance Analysis (FBA), the application of linear programming to metabolic modeling. Finally, we illustrate FBA in an example.

### 5.2.1 Linear programming problems

Linear programs (LPs) are specified for a vector $x$ (of real variables) and involve a vector $c$ (defining the objective function) as well as a matrix $A$ and a corresponding vector $b$ (defining the linear inequalities),

$$ \begin{array}{ll} \max_{x} & c^{\top}x \\ \text{s.t.} & A x \leq b. \end{array} \qquad (5.1) $$

In general, an individual inequality $a^\top x \leq b$ (where the vector $a$ denotes a row of the matrix $A$ and $b$ denotes the corresponding component of $b$) constrains a weighted sum of the variables. However, if $a$ has only one nonzero component, e.g. $a_i = 1$ and $a_j = 0$ for all $j \neq i$, then it restricts a single variable, $x_j \leq b$.

*Towards LP standard form.* Some variables may have nonnegativity constraints. Hence, we order the variables such that

$$ x = \begin{pmatrix} x^{\circ} \\ x^{*} \end{pmatrix}, $$

where the superscripts $\circ$ and $*$ refer to the index sets of nonnegative and “free” variables, $I^\circ$ and $I^*$, respectively. Now, we can formulate the nonnegativity constraints as $x^\circ \geq 0$. To obtain nonnegativity constraints for all variables, we will consider free variables as differences of nonnegative variables. Given the order on the variables, we first write

$$ A x = (A^{\circ} \quad A^{*}) \begin{pmatrix} x^{\circ} \\ x^{*} \end{pmatrix}. $$

Next, for every $i \in I^*$ (for every free variable $x_i$), we introduce $y_i^+, y_i^- \geq 0$ such that $y_i^+ - y_i^- = x_i$. In vector notation, we have

$$ y^+ - y^- = x^*. $$

Note that this transformation is not one-to-one. There are infinitely many pairs $y^+, y^-$ that have the same difference $x^*$. Analogously, for every $i \in I^\circ$ (for every nonnegative variable $x_i$), we introduce $y_i^\circ = x_i \geq 0$ to obtain a uniform notation. That is,

$$ y^{\circ} = x^{\circ}. $$

Now, we can write

$$ A x = (A^{\circ} \quad A^{*} \quad -A^{*}) \begin{pmatrix} y^{\circ} \\ y^{+} \\ y^{-} \end{pmatrix}. $$

In short,

$$ A x = \bar{A} y \quad \text{with} \quad \bar{A} := (A^{\circ} \quad A^{*} \quad -A^{*}) \quad \text{and} \quad y := \begin{pmatrix} y^{\circ} \\ y^{+} \\ y^{-} \end{pmatrix}, $$
<figure><img src="image_42.png" /><figcaption>Figure 6.1: Two natural variants of the glycolysis pathway, named after their discoverers: Embden-Meyerhof-Parnas (EMP) and Entner-Doudoroff (ED)</figcaption></figure>

Fritz Meyerhof, and Karol Parnas (or the EMP pathway for short, see Figure 6.1), was the first metabolic pathway to be discovered by scientists [212]. Often, the pyruvate is reduced to lactate or ethanol, which makes the pathway redox balanced. Therefore, it one of the most common way for producing ATP anaerobically (i.e. without oxygen to serve as an electron acceptor). Another common variant was discovered in 1952 by Nathan Entner and Michael Doudoroff [213] (ED for short). For example, E. coli is capable of metabolizing glucose through both the EMP or the ED variants, and often does so simultaneously [121].

More generally, the overall reaction describing glycolysis is:

$$ \text{Glucose} + 2 \text{NAD(P)}^{+} + n \text{ADP} + n \text{Phosphate} \rightarrow 2 \text{Pyruvate} + 2 \text{NAD(P)}\text{H} + n \text{ATP} + n \text{H}_{2}\text{O} \quad (6.1) $$

where the value of n for the EMP pathway is 2. Ng et al. [214] explored the space of all possible glycolyses (with different values of n), by exhaustively enumerating all glycolytic pathway variants. In order to generate the variants, they adapted a computational method first introduced by Bar-Even et al. [215] for finding alternative carbon fixation cycles – metabolic cycles whose net reaction converts CO₂ into organic compounds. You start by collecting a database of known biochemical reactions (e.g. from a database such as KEGG [216]) and then use a linear-programming algorithm to identify the set of reactions with the minimal sum of fluxes that conform to the predefined net reaction (e.g. 6.1). The objective is somewhat arbitrary, but since solving the LP requires setting an objective, we chose the min-flux as a reasonable proxy for the simplicity of the pathway. In any case, we will soon see how one can iterate through all possible solutions. Ng et al. [214] used this algorithm with the stoichiometry from 6.1 to find all possible glycolysis pathways comprising known enzymatic reactions (see Mathematical Details Box 6.B).

Mathematical details 6.B: Integer cuts can be used for iterating all possible solutions