involving the augmented matrix $\bar{\mathbf{A}}$ and the non-negative vector $\mathbf{y}$. Correspondingly, we have to transform the vector

$$ \mathbf{c} = \begin{pmatrix} \mathbf{c}^{\circ} \\ \mathbf{c}^{*} \end{pmatrix}. $$

Indeed,

$$ \mathbf{c}^{\top}\mathbf{x} = \bar{\mathbf{c}}^{\top}\mathbf{y} \quad \text{with} \quad \bar{\mathbf{c}} := \begin{pmatrix} \mathbf{c}^{\circ} \\ \mathbf{c}^{*} \\ -\mathbf{c}^{*} \end{pmatrix}. $$

As a result, we obtain an LP in standard form,

$$ \begin{array}{ll} \max_{\mathbf{y}} & \bar{\mathbf{c}}^{\top}\mathbf{y} \\ \text{s.t.} & \bar{\mathbf{A}}\mathbf{y} \leq \mathbf{b} \\ & \text{and } \mathbf{y} \geq 0, \end{array} \qquad (5.2) $$

which is equivalent to the original LP (5.1) in the following sense:

* On the one hand, if $\mathbf{x}$ is a feasible solution of (5.1), that is, $\mathbf{A} \mathbf{x} \leq \mathbf{b}$, with objective value $\mathbf{c}^{\top}\mathbf{x}$, then every $\mathbf{y}$ with $\mathbf{y}^{\circ} = \mathbf{x}^{\circ}$ and $\mathbf{y}^{+} - \mathbf{y}^{-} = \mathbf{x}^{*}$ is a feasible solution of (5.2), that is, $\bar{\mathbf{A}} \mathbf{y} \leq \mathbf{b}$ and $\mathbf{y} \geq 0$, with the same objective value $\bar{\mathbf{c}}^{\top}\mathbf{y} = \mathbf{c}^{\top}\mathbf{x}$.

In particular, we can choose $\mathbf{y}^{+}$ and $\mathbf{y}^{-}$ as follows: for every free variable $x_i$, if $x_i \geq 0$, then we set $y_i^{+} = x_i$ and $y_i^{-} = 0$, and if $x_i \leq 0$, then we set $y_i^{+} = 0$ and $y_i^{-} = -x_i$. Equivalently, we can choose $\mathbf{y}^{+}$ and $\mathbf{y}^{-}$ such that $\mathbf{y}^{+} - \mathbf{y}^{-} = \mathbf{x}^{*}$ and $y_i^{+} \cdot y_i^{-} = 0$ for all indices $i \in \mathcal{I}^{*}$.

*Note.* The latter choice involves if-then statements or nonlinear equations. Hence, it cannot be written in terms of linear inequalities, that is, it cannot be directly incorporated into an LP. Still, the choice can be used to obtain equivalent (but simpler) formulations of LPs. See the next subsection, in particular, the formulation of enzyme constraints in FBA.

* On the other hand, if $\mathbf{y}$ is a feasible solution of (5.2), then $\mathbf{x}$ with $\mathbf{x}^{\circ} = \mathbf{y}^{\circ}$ and $\mathbf{x}^{*} = \mathbf{y}^{+} - \mathbf{y}^{-}$ is a feasible solution of (5.1) and $\mathbf{y}$ and $\mathbf{x}$ have the same objective value.

### 5.2.2 Metabolic modeling via linear programs aka Flux balance analysis

Flux balance analysis (FBA) studies metabolic models via linear programs. The feasible solutions (fluxes) are specified by

- homogeneous equations arising from mass-balance constraints,

- homogeneous inequalities (nonnegativity conditions) arising from irreversible reactions,

- and extra inhomogeneous (in-)equalities encoding lower/upper bounds for individual fluxes or enzyme constraints involving weighted sums of fluxes.

(See also the previous chapter.)

The objective function can be an individual flux or a weighted sum. Hence, a general FBA problem can be written as

$$ \begin{array}{ll} \max_{\mathbf{v}} & \mathbf{c}^{\top}\mathbf{v} \\ \text{s.t.} & \mathbf{N}\mathbf{v} = 0, \mathbf{v} \to \geq 0, \\ & \text{and } \mathbf{G}\mathbf{v} \geq \mathbf{h}. \end{array} \qquad (5.3) $$

*Two types of extra constraints.* As already mentioned, the inhomogeneous inequalities may involve lower/upper bounds for (all) individual fluxes,

$$ \ell \leq \mathbf{v} \leq \mathbf{u}, \qquad (5.4) $$

so-called box constraints. Componentwise, $\ell_i \leq v_i \leq u_i$ for reaction $i$. (If there is no lower or upper bound, one may set $\ell_i = -\infty$ or $u_i = \infty$.) Of course, the bounds must be consistent with the irreversibility
The linear problem can be described by:

$$
\begin{array}{ll}
\text{minimize} & \displaystyle\sum_{i} v_{i} \\
\text{subject to} & \mathbf{N}\mathbf{v} = 0 \\
& \forall i \quad 0 \le v_{i} \le \beta \\
& v_{\text{glycolysis}} = -1
\end{array}
\qquad (6.2)
$$

where v is the flux variable, and N is comprised of the universal stoichiometric matrix , and in addition one reaction (whose flux is denoted vglycolysis) which has the stoichiometry of Eq. (6.1). The constraint vglycolysis = -1 ensures that the sum of all active reactions except for vglycolysis will together form a full glycolysis pathway, since their net reaction has to balance the stoichiometry of vglycolysis given the mass balance constraint Nv = 0. β given the upper bound on the flux for all reactions. For simplicity, we assume that all fluxes are positive and that reversible reactions are split into their two opposing directionalities . β is a tunable parameter that is an upper bound on all the fluxes in the solution pathways. Setting it too low would exclude solutions with complex stoichiometries. On the other hand, a very high value would increase the complexity of the search and lead to very long run-times. Typically, we choose β = 10 which is a good balance between the two extremes. Finally, we set the objective function (Σi vᵢ) to minimize the sum of fluxes. As we will explain shortly, we can iterate through all possible solutions and therefore the objective will only determine the order at which we find them.

To find all possible glycolysis pathways comprising known enzymatic reactions, Ng et al. [214] iteratively introduced
constraints in order to exclude all previous solutions and find the next optimal one [217]: to exclude a solution,
they add an integer cut, which is an inequality constraint ensuring that the number of active reactions is strictly
larger than the sum over their indicator variables (boolean variables that are equal to 1 if the reaction is active, i.e.
carries a nonzero flux). Therefore, at least one of those reactions must be inactive in all future solutions. This is
quite similar to constrained Minimal Cut Sets (cMCS) which were introduced in Chapter 4.4 as a way of exploring
the flux space.

Formally, if {P₀, P₁... Pₘ} are the set of solutions already discovered by our algorithm (where ∀j Pⱼ ⊆ {0, ..., n},
i.e. each solution is a set of integers which correspond to indices of active reactions) then the added constraints
will be:

$$
\begin{align*}
\forall i & \quad z_i \in \{0, 1\} \\
\forall i & \quad v_i - \beta z_i \le 0 \\
\forall j & \quad \sum_{i \in P_j} z_i < \|P_j\| \tag{6.3}
\end{align*}
$$

where ||P_j|| is the length of pathway j (i.e. the number of reactions). The z_i are boolean reaction indicators, i.e. z_i
must be equal to 1 if a reaction is active (v_i > 0). The final set of constraints eliminate P_j and any pathway which
is a superset of P_j from the solution space. Using this extra set of constraints iteratively, each time generating the
next pathway and adding it to the excluded list, will eventually go through all possible solutions (by increasing order
of their sum of fluxes). It is important to note that using integer cuts requires switching to an MILP (Mixed-Integer
Linear Program) solver, which is computationally much more demanding and typically requires a commercial license.

The objective set by the linear problem (6.2) is minimizing the sum of fluxes, which corresponds to pathways with fewer
reactions and low fluxes in each one. As discussed in 5.2, this objective is only a crude proxy for the efficiency of a
pathway, and its only purpose is to get the pathway solutions in a relatively logical order. Although we have discussed
global enzyme constraints in previous chapters (such as molecular crowding and proteome allocation), when comparing
pathways we will focus only on the efficiency of the pathway itself. This will allow us to compare pathways without
thinking about the rest of the cell or a specific metabolic context. But how can one quantify the efficiency of a pathway?
constraints.

Additionally/alternatively, one may consider enzyme constraints. First, for every reaction $i \in \mathcal{R} \to \cup \mathcal{R}$, one introduces the concentration $e_i$ of the corresponding enzyme. Then, for every irreversible reaction $i \in \mathcal{R} \to$, enzyme kinetics implies the inequality $0 \le v_i \le k_i^{\text{cat}, \to} e_i$. Similarly, for every reversible reaction $i \in \mathcal{R} \rightleftharpoons$, one obtains

$$ \left\{ \begin{array}{ll} v_i \le k_i^{\text{cat}, \to} e_i & \text{if } v_i \ge 0, \\ -v_i \le k_i^{\text{cat}, \leftarrow} e_i & \text{if } v_i \le 0. \end{array} \right. \qquad (5.5a) $$

Finally, using nonnegative weights $\omega_i$, one formulates a capacity constraint

$$ \sum_{i \in \mathcal{R} \to \cup \mathcal{R} \rightleftharpoons} \omega_i e_i = \omega^{\top} e \le \Omega \qquad (5.5b) $$

(or several such constraints). Clearly, the inequalities (5.5a) involve if-then statements. Hence, they cannot be directly incorporated into an LP. To resolve this problem, we perform reaction splitting, which ultimately leads to an LP in standard form. In the following, we describe this transformation in detail.

**Reaction splitting.** As described in the previous chapter, for every reversible reaction $i \in \mathcal{R} \rightleftharpoons$ with net reaction rate $v_i$, we define a forward reaction "rate" $w_i^{\to} \ge 0$ and a backward reaction "rate" $w_i^{\leftarrow} \ge 0$ such that $v_i = w_i^{\to} - w_i^{\leftarrow}$. In vector notation, $\mathbf{v}^{\to} = \mathbf{w}^{\to} - \mathbf{w}^{\leftarrow}$ with $\mathbf{w}^{\to}, \mathbf{w}^{\leftarrow} \ge 0$. Again, we 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. Further, for every irreversible reaction $i \in \mathcal{R} \to$, we write $v_i = w_i^{\to} \ge 0$ to obtain a uniform notation. That is, $\mathbf{v}^{\to} = \mathbf{w}^{\to} \ge 0$. We introduce the augmented stoichiometric matrix $\bar{\mathbf{N}}$ and the corresponding non-negative flux vector $\mathbf{w}$,

$$ \bar{\mathbf{N}} := (\mathbf{N} \to \quad \mathbf{N} \rightleftharpoons \quad -\mathbf{N} \leftarrow) \quad \text{and} \quad \mathbf{w} := \begin{pmatrix} \mathbf{w}^{\to} \\ \mathbf{w}^{\leftarrow} \\ \mathbf{w}^{\leftarrow} \end{pmatrix}, $$

and find

$$ \mathbf{N} \mathbf{v} = \bar{\mathbf{N}} \mathbf{w}. $$

Analogously, we introduce the augmented matrix $\bar{\mathbf{G}}$ and the (objective) vector $\bar{\mathbf{c}}$,

$$ \bar{\mathbf{G}} := (\mathbf{G} \to \quad \mathbf{G} \rightleftharpoons \quad -\mathbf{G} \leftarrow) \quad \text{and} \quad \bar{\mathbf{c}} := \begin{pmatrix} \mathbf{c}^{\to} \\ \mathbf{c}^{\leftarrow} \\ -\mathbf{c}^{\to} \end{pmatrix}, $$

and obtain the LP

$$ \begin{array}{ll} \max_{\mathbf{w}} & \bar{\mathbf{c}}^{\top} \mathbf{w} \\ \text{s.t.} & \bar{\mathbf{N}} \mathbf{w} = 0, \mathbf{w} \ge 0, \\ & \text{and } \bar{\mathbf{G}} \mathbf{w} \ge \mathbf{h}, \end{array} \qquad (5.6) $$

which is equivalent to the original LP (5.3).

**Choice of w.** To incorporate the enzyme constraints (5.5a), we choose $\mathbf{w}^{\to}$ and $\mathbf{w}^{\leftarrow}$ as follows: for every reversible reaction $i \in \mathcal{R} \to$, if $v_i \ge 0$, then we set $w_i^{\to} = v_i$ and $w_i^{\leftarrow} = 0$, and if $v_i \le 0$, then we set $w_i^{\to} = 0$ and $w_i^{\leftarrow} = -v_i$. We find

$$ \left\{ \begin{array}{ll} w_i^{\to} = w_i^{\leftarrow} + v_i \le v_i \le k_i^{\text{cat}, \to} e_i & (\text{and } w_i^{\leftarrow} = 0) \quad \text{if } v_i \ge 0, \\ w_i^{\leftarrow} = w_i^{\to} - v_i \le -v_i \le k_i^{\text{cat}, \leftarrow} e_i & (\text{and } w_i^{\to} = 0) \quad \text{if } v_i \le 0 \end{array} \right. $$

and hence

$$ w_i^{\to} / k_i^{\text{cat}, \to} + w_i^{\leftarrow} / k_i^{\text{cat}, \leftarrow} \le e_i, \qquad (\rightleftharpoons) $$

independently of the sign of $v_i$ (that is, not involving an if-then statement). To summarize, the if-then
The next section will be dedicated to exactly this question.

## 6.2. Pathway efficiency - some basic notions and thoughts

For glycolysis alone, Ng et al. [214] found 11,916 alternatives that produce at least one mole of ATP per mole of glucose. These include, of course, the EMP pathway. Although evolution can explore these options, natural selection typically converges on one or a few efficient variants. This does not mean that every single pathway observed in nature must be optimal, but we generally expect cells hosting highly inefficient pathways to eventually become extinct. Iacometti et al. [218] tested this experimentally by knocking out the EMP pathway from E. coli and forcing the cells to use the alternatives that naturally exist in this bacterium. In all cases, growth rates were slower than in the wild-type.

Before we discuss other examples for metabolic pathways, we need to define what we mean by “efficiency”. There are several criteria one should consider:

* Low consumption rate of the substrate

* High generation rate of the product

* High regeneration rate or low consumption rate of the co-factor

* Small number of steps [219]

* Higher thermodynamic forces [220, 221]

* High enzyme turnover numbers

* High enzyme saturation levels

Some of these criteria refer to the cost (or investment) of the pathway, while others reflect the benefit (or profit) to the cell. By considering two common scenarios – single nutrient limitation or exponential growth in rich media – we can focus on two simple criteria which provide good measures of efficiency.

When the availability of a single nutrient is limiting growth, maximizing the molar yield (i.e. the number of moles of product generated for each mole of the nutrient) becomes the important feature. Yield is rather straightforward to calculate, as it is a direct outcome of the stoichiometry of the pathway. For example, anaerobic fermentation is often compared to respiration and deemed inefficient since it yields two moles of ATP per glucose, instead of ≈30 [222].

On the other hand, when conditions are good, such as during exponential growth in rich media, minimizing the total number of proteins required is often the objective which determines growth rate. Here, we will be using the enzyme demand (e.g. in grams of protein) per unit of flux (typically, in mmol per hour per gram of cell dry weight). In fact, the enzyme demand per flux, as an objective, takes into consideration both the cost (protein) and the benefit (flux). Importantly, these two criteria scale linearly with respect to each other: doubling the amount of all enzymes without changing any of the metabolite concentrations would directly double the flux in the pathway. Therefore, this measure of efficiency is independent of the magnitude of the flux in the pathway. But, as we will see shortly, enzyme demand is a non-linear function, making it trickier to compute compared to other constraint-based problems such as ones we've seen in previous chapters.

Notably, these two measures of efficiency are not only useful for evolutionary processes, but for bioengineering as well. Obviously, the molar yield has economical implications when, for example, producing ethanol from sugar. However, the rate of a bioprocess is important as well due to the costs involved, e.g. for maintaining an operational bioreactor. One can imagine a computational model that accurately predicts the enzyme demand per flux of a pathway. Choosing the pathways with the lowest demand would be a good strategy for increasing the overall rate of bioproduction [223].

We define the enzyme demand per unit flux as the total amount of enzyme (in grams of protein) that is required to catalyze all of the pathway reactions at their required rates. We start by deriving a formula for the demand of a single enzymatic reaction. Consider an enzyme-catalyzed reaction:

$$S \rightleftharpoons P \qquad (6.4)$$
statement (5.5a) for v together with our choice for w (another if-then statement) imply the inequality ($\Leftrightarrow$) which can be incorporated into the LP (5.6). In fact, we further rewrite (=$\Leftrightarrow$) as

$$w_i^- / k_i^{\text{cat}, \to} \le e_i^-, \quad w_i^- / k_i^{\text{cat}, \leftarrow} \le e_i^-, \quad \text{and} \quad e_i^- + e_i^+ \le e_i,$$

where we introduce the forward and backward enzyme "concentrations", $e_i^-$ and $e_i^+$, respectively.

The treatment of the other extra constraints is simple. First, for every irreversible reaction $i \in \mathcal{R}^\to$, we get $w_i^\to = v_i \le k_i^{\text{cat}, \to} e_i$ or, equivalently, $w_i^\to / k_i^{\text{cat}, \to} \le e_i$. Altogether,

$$w_i \le k_i^{\text{cat}} \bar{e}_i, \quad \text{for } i \in \mathcal{R}^\to \cup \mathcal{R}^\leftarrow \cup \mathcal{R}^\leftarrow, \qquad (5.7a)$$

where we use the vectors

$$\mathbf{k}^{\text{cat}} := \begin{pmatrix} \mathbf{k}^{\text{cat}, \to} \\ \mathbf{k}^{\text{cat}, \leftarrow} \\ \mathbf{k}^{\text{cat}, \leftarrow} \end{pmatrix} \quad \text{and} \quad \bar{\mathbf{e}} := \begin{pmatrix} \mathbf{e}^\to \\ \mathbf{e}^\leftarrow \\ \mathbf{e}^\leftarrow \end{pmatrix}$$

to streamline the notation. Second, we adapt the capacity constraint (5.5b). Using the augmented vector

$$\bar{\omega} := \begin{pmatrix} \omega^\to \\ \omega^\leftarrow \\ \omega^\leftarrow \end{pmatrix},$$

we obtain

$$\sum_{i \in \mathcal{R} \to \cup \mathcal{R} \leftarrow \cup \mathcal{R}^-} \bar{\omega}_i \bar{e}_i = \sum_{i \in \mathcal{R} \to} \omega_i e_i + \sum_{i \in \mathcal{R}^=} \omega_i \underbrace{(e_i^- + e_i^+)}_{\le e_i} \\ \le \sum_{i \in \mathcal{R} \to \cup \mathcal{R}^=} \omega_i e_i \\ \le \Omega,$$

in short,

$$\sum_{i \in \mathcal{R} \to \cup \mathcal{R} \leftarrow \cup \mathcal{R}^-} \bar{\omega}_i \bar{e}_i = \bar{\omega}^{\top} \bar{e} \le \Omega. \qquad (5.7b)$$

Of course, box constraints (5.4) for v can also be written in terms of w: for the irreversible reactions, we get $\ell^- \le w^- \le u^- \quad (\text{with } \ell^- \ge 0)$, and for the reversible reactions, we get $\ell^+ \le w^- - w^- \le u^+ \quad (\text{with } \ell^+ < 0 \text{ and } u^+ > 0)$. The latter inequalities can be further simplified using our choice of w.

**Summary.** We have shown that an FBA problem (5.3) involving box constraints (5.4) and enzyme constraints (5.5) can be transformed into an equivalent LP problem (5.6) in standard form. In particular, for every feasible v (and e) in $\mathbb{N}v = 0$, $v^- \ge 0$, and (5.5), there is a feasible w (and $\bar{e}$) in $\overline{\mathbb{N}}w = 0$, $w \ge 0$, and (5.7), namely, our choice of w. Conversely, for every feasible w (not necessarily in the form of our choice) in the standard LP problem, there is a feasible v in the original FBA problem.

### 5.2.3 Optimization problems for metabolic fluxes

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.1) and (4.16). 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), \quad \text{such that} \quad \mathbf{A}v \ge \mathbf{b}, \qquad (5.8)$$
<figure><img src="image_43.png" /><figcaption>Figure 6.2: Enzyme cost in metabolism. (a) Enzyme-specific flux depends on a number of physical factors. Under ideal conditions, an enzyme molecule catalyses its reaction at a maximal rate given by the enzyme's forward catalytic constant (blue). The rate is reduced by microscopic reverse fluxes (magenta) and by incomplete saturation with substrate, causing waiting times between reaction events, or by enzyme inhibition or incomplete activation (red). (b-c) On a logarithmic scale, catalytic rates and enzyme demand can be split into sums of efficiency terms. With lower catalytic rates, larger amounts of enzyme are required for realizing the same metabolic flux.</figcaption></figure>

where $s$ and $p$ will be the concentrations of the substrate (S) and product (P) respectively, and $E$ the concentration of the enzyme which catalyzes this reaction (for simplicity, we drop the tot subscript from $E_{tot}$). Here, we will be using the factorized rate law (Eq. 3.15), but other kinetic rate laws would produce similar results. The rate of a reaction is given by:

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

where $k_{\text{cat}}^{+}$ is the forward turnover rate, $K_{\text{S}}$ and $K_{\text{P}}$ are the Michaelis-Menten constants for the S and product P, and $\Delta G'_{r}$ is the Gibbs free energy. So, the minimal amount of enzyme that is required for reaching a given rate $v$ is:

$$q \equiv v \cdot h \cdot \frac{1}{k_{\text{cat}}^{+}} \cdot \frac{1 + p/K_{\text{P}} + s/K_{\text{S}}}{s/K_{\text{S}}} \cdot \left(1 - e^{\Delta G'_{r}/RT}\right)^{-1}, \quad (6.6)$$

where $h$ is a number converting enzyme concentration $e$ into enzyme amount $q$ (for example, the enzyme molecular mass). For an illustration, see Figure 6.2. Summing up the demand across all the reactions in the pathway (each with its own rate, kinetic parameters, and substrate/product concentrations) will produce the total enzyme demand. Looking at this function, we can already make some interesting observations. First, the kinetic parameters ($k_{\text{cat}}^{+}$, $K_{\text{P}}$, and $K_{\text{S}}$) can be treated as constants since they change only in evolutionary timescales, and we often assume that existing enzymes already have near-optimal kinetics (although that's not always the case). Since we care about the demand *per pathway* flux one can, without loss of generality, set $v$ to 1. However, if the pathway requires a non-trivial ratio between some reactions, the value of $v$ can be different based on the stoichiometry. Finally, the thermodynamic term, i.e. $1 - e^{\Delta G'_{r}/RT}$ (which we will discuss in more detail in the following section, 6.3), is a function of the metabolite concentrations and the $K_{\text{eq}}$, which is another constant. So, generally speaking, enzyme demand is defined by a set of constants that are unique to each pathway, and variables that represent the metabolite concentrations. Since these concentrations are subject to change depending on the growth conditions, we often treat them as optimization variables and try to find the minimal demand possible within certain constraints. In Section 6.4, we will see a general method for finding the minimal value using convex optimization.

Most of the proposed criteria for good pathways have either to do with material investments (such as substrate, cofactor, or energy demand) or with “machine investments”, that is, enzyme demands. Enzyme demands, in turn, depend on
<figure><img src="image_56.png" /></figure>

Economics analogy 5.A Optimization and economic planning in the Soviet Union

Linear programming as an algorithmic approach to solving constrained linear optimization problems was first developed by soviet mathematician and economist Leonid Kantorovich in the 1930s [134, 135]. Kantorovich was tasked with helping to optimize production in the soviet plywood industry, but soon discovered that the underlying problems 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 [136]. 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).

with

$$A = \begin{pmatrix} N_{\text{int}} \\ -N_{\text{int}} \\ \mathbf{I} \\ \mathbf{G} \end{pmatrix}, \quad b = \begin{pmatrix} 0 \\ 0 \\ 0 \\ \mathbf{h} \end{pmatrix}. \qquad (5.9)$$

Recall that $N_{\text{int}}v = 0$ models the steady-state assumption, while the multiplication with the identity matrix $(I_{n \times n}v \ge 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 \ge 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.10)$$

where coefficients $c_i$ weigh the relevance of the different reaction rates in the objective function. Problems of the form (5.8), (5.9), and (5.10) 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 [137, 138]. 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
Mathematical details 6.C: Factorized rate laws and enzyme cost function

According to Eq. (6.13), reversible rate laws can be factorized into five terms that depend on metabolite concentrations in different ways [102]. For a reaction S ⇌ P with reversible Michaelis-Menten kinetics Eq. (6.11), a driving force θ = −Δ<sub>r</sub>G'/RT, and a prefactor for non-competitive inhibition, the rate law can be written as

$$
v = E \cdot k_{\text{cat}}^{+} \cdot \underbrace{[1 - e^{-\theta}]}_{\eta^{\text{for}}} \cdot \frac{s/K_{\text{S}}}{1 + s/K_{\text{S}} + p/K_{\text{P}}} \cdot \underbrace{\frac{1}{1 + x/K_{\text{I}}}}_{\eta^{\text{reg}}}
$$

$$
\text{Rate} = \begin{array}{c} \text{enzyme} \\ \text{level} \end{array} \cdot \begin{array}{c} \text{forward catalytic} \\ \text{constant} \end{array} \cdot \begin{array}{c} \text{thermodynamic} \\ \text{factor} \end{array} \cdot \begin{array}{c} \text{saturation} \\ \text{factor} \end{array} \cdot \begin{array}{c} \text{regulation} \\ \text{factor} \end{array}
$$

with inhibitor concentration $x$. The product of the first two terms, $E$ and $k_{\text{cat}}^{+}$, represents the maximal velocity, i.e. the rate at full substrate-saturation without backward flux and without enzyme inhibition. The following factors decrease this velocity for different reasons: $\eta^{\text{for}}$ describes a decrease due to backward fluxes, $\eta^{\text{sat}}$ – the decrease due to incomplete substrate saturation, and $\eta^{\text{reg}}$ – the decrease due to small-molecule regulation (see Figure b). While $k_{\text{cat}}^{+}$ is an enzyme-specific constant (yet, dependent on conditions such as pH, ionic strength, or molecular crowding in cells), the efficiency factors are concentration-dependent, unitless, and can vary between 0 and 1. The thermodynamic factor $\eta^{\text{for}}$ depends on the driving force (and thus, indirectly, on metabolite concentrations), and the equilibrium constant is required for its calculation. The saturation factor $\eta^{\text{sat}}$ depends directly on metabolite levels and contains the $K_{\text{M}}$ values as parameters. Enzyme regulation by small molecules yields additive or multiplicative terms in the rate law denominator, which in our example and can be captured by a separate factor $\eta^{\text{reg}}$. The enzyme cost for a flux $v$, with an enzyme burden $h_e$, can be written as

$$
q = h_e \cdot E = h_e \cdot v \cdot \frac{1}{k_{\text{cat}}^{+}} \cdot \underbrace{\frac{1}{[1 - e^{-\theta}]}}_{1/\eta^{\text{for}}} \cdot \underbrace{\frac{1 + s/K_{\text{S}} + p/K_{\text{P}}}{s/K_{\text{S}}}}_{1/\eta^{\text{sat}}} \cdot \underbrace{\frac{[1 + x/K_{\text{I}}]}{1/\eta^{\text{reg}}}}_{\eta^{\text{reg}}}
$$

and contains the terms from the rate law in inverse form. The first factors, $h_e v / k_{\text{cat}}^{+}$, define a minimum enzyme cost, which is then increased by the following efficiency factors. By omitting some of these factors, one can construct simplified enzyme cost functions with higher specific rates, or lower enzyme demands (compare Figure 6.2b). For a closer approximation, the factors may be substituted with constant numbers between 0 and 1.

pathway length, enzyme masses, and enzyme efficiency, and therefore on rate laws (where $k_{\text{cat}}$ values, thermodynamic forces, and metabolite concentrations come into play).

In fact, many criteria which we discussed earlier as indicators of efficiency are actually an approximation of the enzyme demand under certain assumptions. For example, the number of steps is proportional to the total demand if all enzymes have exactly the same $k_{\text{cat}}^{+}$, saturation, and thermodynamics. Therefore, it is quite a useful rule-of-thumb in case not much else is known about the enzymes themselves.

A better approximation, denoted *Pathway Specific Activity*, was used by [215] to compare CO₂ fixation cycles. If we assume that all enzymes are fully saturated and irreversible, the demand would be a direct function of the individual enzyme specific activities (specifically, proportional to the sum of all their reciprocal values).

But even if we know nothing about the enzyme kinetic parameters, thermodynamics alone can provide us with useful information with which to grade pathways. Specifically, the $K_{\text{eq}}$ of a reaction is a universal constant that is not affected by enzymes, but rather determined solely by the chemical structures of the substrates and products.

In the following sections, we will focus on enzyme use efficiency as a main objective and consider a thermodynamic approximation, relating enzyme demands to thermodynamic forces. For linear metabolic pathways, optimal enzyme profiles (and the associated metabolite profiles and enzyme costs) can be computed with closed formulae. We will also
<figure><img src="image_57.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 $n_1$ and $n_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 $n_X$ molecules of ATP.</figcaption></figure>

$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 $n_X$ molecules of ATP with flux value $v_X$. We also introduce as a linear objective function the total rate of ATP production, $vATP$. Since in this example, reactions $v_1$ and $v_3$ produce ATP with stoichiometric coefficients $n_1$ and $n_3$, respectively, the total rate of ATP production is given by $v_{ATP} = n_1v_1 + n_3v_3 - n_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 $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}
\displaystyle \max_{v} v_{ATP} = n_{1}v_{1} + n_{3}v_{3} - n_{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} \ge v_{0}, \\
v_{4}^{ub} \ge v_{4}, \\
v_{0}, v_{1}, v_{2}, v_{3}, v_{4} \ge 0.
\end{array} 
\qquad (5.11) 
$$

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 $n_3 = n_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}$.
discuss a way to compute optimal enzyme profiles numerically, for networks of any shape and size, as long as the flux mode is known.

## 6.3. The role of thermodynamics

In general, when considering larger metabolic networks, thermodynamic feasibility can play an important or even crucial role in determining which pathways are used. In this section we will discuss this role more explicitly and see how thermodynamics can still give us useful insights about pathway efficiency even when no other kinetic data is available.

Why are thermodynamic driving forces a meaningful criterion for good pathways? In brief, the driving forces, defined as $\theta \equiv -\Delta G'/RT$, play a double role: first, they determine whether or not a pathway flux is feasible at all, given the metabolite concentrations at the pathway boundary (i.e. the metabolites that form connections to the broader metabolic network); and second, in case the pathway is feasible, driving forces can affect enzyme efficiency and, consequently, the enzyme demand for a given desired pathway flux. In Chapter 3, we learned that $\Delta G'$, and hence the driving force $\theta$, depends on the equilibrium constant $K_{eq}$ of the reaction and on the substrate and product concentrations. We also learned that for a flux in forward direction, the driving force must be positive. Beyond that, the efficiency of an enzyme is proportional to $\eta^{\text{for}}(\theta) = 1 - e^{-\theta}$, a function that ranges between 0 (for $\theta = 0$, reactions in thermodynamic equilibrium) and 1 ($\theta \gg 1$, reactions far from equilibrium). Let us now see how this non-equilibrium relation affects pathway efficiency.

### 6.3.1. Kinetics and driving forces

We should remind ourselves some of the lessons learned in Chapter 3. Specifically, recall the factorized rate law [102] with a reversibility term that is an explicit function of the Gibbs energy (Eq. 3.15):

$$v = E \cdot k_{\text{cat}}^{+} \cdot \frac{\prod_{i} s_{i}^{\nu_{i}} / K_{s}}{1 + \prod_{j} p_{j}^{\nu_{j}} / K_{p} + \prod_{i} s_{i}^{\nu_{i}} / K_{s}} \cdot (1 - e^{\Delta G'_{r} / RT}). \quad (6.7)$$

The enzyme mechanism behind this formula assumes fast binding and unbinding of substrate and product, and a slow reversible conversion step (of bound substrate into bound product). Note that here we generalize the rate law for cases with more than one substrate and one product, where $\nu_i$ and $\nu_j$ are the stoichiometric coefficients of substrates and products, respectively<fnref n="1" />. This generalization is one out of many, and corresponds to the assumption that all reactants bind independently to the enzyme (and at random order). We focus on this rate law because it is one of the simplest, but the theoretical results in this chapter apply to most other generalizations as well (e.g. convenience kinetics [224]).

According to the definition of $k_{\text{cat}}^{+}$, and also by noticing that the middle and rightmost terms in Eq. (6.7) are each smaller than 1, the rate of an enzymatic reaction is bounded by $v \le E \cdot k_{\text{cat}}^{+}$ (see Mathematical Details Box 6.C for a detailed explanation). However, the additional terms are often much lower than 1, which means that the rate does not reach its maximum. If we try to measure the apparent catalytic rate by dividing the rate by the enzyme abundance ($k_{\text{app}} = v/E$) we would typically get a value that is lower than $k_{\text{cat}}^{+}$, while only in rare “ideal” cases, $k_{\text{app}}$ would approach the $k_{\text{cat}}^{+}$. In fact, this reasoning was used by Davidi et al. [31] to estimate the $k_{\text{cat}}^{+}$ values of more than 100 enzymes in E. coli, where they sampled many growth conditions and took the maximum $k_{\text{app}}$ as the estimate.

As discussed in Section 3.3.2, the factorized rate law has a thermodynamic perspective based on the flux-force relationship, where we view the reversibility term as a “penalty” for the fact that by lowering the energy barrier, enzymes must catalyze reactions in both directions. When the driving force ($\theta$) is low, the reverse reaction flux can become significant and lower the net flux. On the other hand, if the driving force is large enough, this term can be ignored and the rate law

<footnote n="1">¹In general, reaction stoichiometries can be arbitrarily scaled. For example, instead of a reaction $2 A \rightarrow B$, we may write $A \rightarrow \frac{1}{2} B$ for convenience, which will only lead to a scaling factor in the reaction rate. However, this holds only if reaction stoichiometries are used to describe mass-balance. In cases like Eq. (6.7), where stoichiometries appear in kinetic rate laws or in thermodynamic balances, we do not have this choice. In these cases, the stoichiometries must reflect the molecularities, that is, the actual number of reactant molecules involved in the enzymatic reaction.</footnote>
## 5.3 Choice of objective functions in Flux Balance Analysis

Solving the constraint-based optimization problem of (5.8) 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 [139, 140, 141].

**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 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[142]. 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 [143, 144, 145].

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 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 10 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
<figure><img src="image_44.png" /><figcaption>Figure 6.3: The thermodynamic efficiency term $\eta^{\text{for}}$ and some approximations. (a) In a given reaction, the thermodynamic efficiency term $\eta^{\text{for}} = 1 - e^{-\theta}$ (solid line) can vary between 0 and 1 depending on the driving force $\theta$. Small driving forces make the enzyme inefficient, since $\eta^{\text{for}} \to 0$, while for large forces, thermodynamics does not play a role as $\eta^{\text{for}} \to 1$. The dashed lines show two linear approximation that hold always as bounds, but can also be used as good approximations for small or large $\theta$ values, respectively: $(1 - e^{-\theta}) < \theta$ and $(1 - e^{-\theta}) < 1$. (b) The reciprocal value $1/\eta^{\text{for}}$ is one of the factors determining enzyme demand. The solid line shows the thermodynamic demand factor $1/\eta^{\text{for}}$, while the dashed lines show the resulting approximations $1/\eta^{\text{for}} > 1/\theta$ and $1/\eta^{\text{for}} > 1$, corresponding respectively to the enzyme demand approximations $E \ge \frac{v}{k_{\text{cat}}\theta}$ and $E \ge \frac{v}{k_{\text{cat}}}$.</figcaption></figure>

resembles irreversible kinetics.

So far we’ve seen that increasing the driving force of a single reaction translates to a better enzyme efficiency and lower demand. If we consider whole pathways, ones whose overall driving force is larger have more of it to distribute among the reactions and therefore should also have higher efficiencies overall. However, using “too much” driving force can also have downsides. Using a larger amount of the Gibbs energy to drive the pathway reactions means that less of that energy would go for building biomass or currency metabolites such as ATP. An example for this trade-off between the efficiency of single enzymes (in terms of backward rates) and the overall pathway efficiency (in terms of ATP yield) was demonstrated by Flamholz et al. [225] who analyzed two versions of the famous glycolytic pathway (see Figure 6.1 below).

### 6.3.2. Small driving forces should be avoided

With the factorized rate law 6.7, we can approximate the reaction rates by $v \le E k_{\text{cat}} (1 - e^{-\theta})$ (where we assume positive fluxes by convention). The thermodynamic efficiency $\eta^{\text{for}} = 1 - e^{-\theta}$ plays a prominent role. As shown in Figure 6.3, this formula yields two important approximations: for small forces $\theta$, that is, close to equilibrium, we obtain $\eta^{\text{for}} \approx \theta$, while for large forces, that is, for strongly forward-driving reactions, we obtain $\eta^{\text{for}} \approx 1$. In fact, both approximations also serve as upper bounds across all $\theta$ values. What does this mean? Far from equilibrium, the thermodynamic term does not play a role and can be ignored. Close to equilibrium, in contrast we obtain a simple approximation for fluxes

$$v < E \cdot k_{\text{cat}}^{+} \cdot (1 - e^{-\theta}) < E \cdot k_{\text{cat}}^{+} \cdot \theta \quad (6.8)$$

and hence for the enzyme demand

$$E > \frac{v}{k_{\text{cat}}^{+} \cdot (1 - e^{-\theta})} > \frac{v}{k_{\text{cat}}^{+} \cdot \theta} \quad (6.9)$$

As $\theta$ goes to zero, the enzyme demand (for a given desired flux) goes to infinity. We already know the reason from Chapter 3: the driving force determines the ratio of forward and reverse one-way fluxes, $\frac{v_{+}}{v_{-}} = e^{\theta}$. If $\theta$ comes close to zero, their relative difference becomes very small, and in order to obtain a given net flux $v = v_{+} - v_{-}$, both $v_{+}$ and $v_{-}$
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 maximized. 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 [146]. 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.4 Enzyme-constrained FBA

In its most simple form, flux balance analysis requires a stoichiometric matrix, an objective function, and at least one flux constraint to ensure that the problem is bounded. Solving an FBA problem allows for the prediction of intracellular fluxes and essential gene knockouts, given measured uptake and secretion rates.

However, whilst classical FBA can capture the effects of essential gene knockouts well, it falls down when it comes to non-lethal knockdowns, and the prediction of growth phenotypes. For example, overflow metabolism in E. coli, and similarly the Crabtree and Warburg effects in S. cerevisiae and cancer cells respectively, cannot be captured in FBA models without ad-hoc flux constraints being imposed. These names refer to the seemingly wasteful strategy of cells at high growth rates using a combination of respiration and fermentation, despite the higher ATP-yield of respiration.

It has been proposed that overflow metabolism results from optimal protein allocation in the cell [147, 148]. In FBA models, we capture this by imposing total proteome constraints to perform enzyme-constrained flux balance analysis (ecFBA). The usual formulation for ecFBA can be written as follows [149]:

$$ 
\begin{array}{ll}
\max_{v,e} & v_r \\
\text{s.t.} & (C1) \quad N v = 0, \\
& (C2) \quad v_i \ge 0 \forall i, \\
& (C3) \quad v_i \le k_i^{\text{cat}} \cdot e_i \quad \forall i \in R \\
& (C4) \quad \sum_{i \in R_k} e_i \le E_k \quad \forall k \in [1, \dots, K]
\end{array}
\qquad (5.12) $$

Here, we wish to maximise flux through the objective reaction $v_r$, subject to four conditions. The matrix $N$ is our stoichiometric matrix, with all reversible reactions split into a forward and a reverse reaction, and the condition $C2$ ensures that all fluxes are positive. In this formulation, we give all metabolic reactions an associated catalysing enzyme, and stipulate that the flux through a reaction is equal to the concentration of this enzyme multiplied by the apparent turnover number ($k^{\text{cat}}$) value (C3). Finally, we constrain total proteome constraints, in the form of enzyme pools. The total concentration of enzymes in the $k$-th enzyme pool must not exceed the constraint $E_k$ (C4).

Predictions using ecFBA do not rely on the input of flux constraints, but rather good estimates for the total protein in different cellular compartments (for example the membrane and the cytosol), as well as the $k_{\text{cat}}$-values (for details, see Chapter 10).

The ecFBA formulation ensures that all metabolic enzymes have an associated cost, relative to the gene
<figure><img src="image_45.png" /></figure>

Figure 6.4: Thermodynamic forces, enzyme efficiency, and enzyme demand in a linear chain of reactions. The plot in the center represents two possible profiles of the thermodynamic driving forces (blue and red). The curves describe the cumulative $\Delta G'$ values: while the total $\Delta G'$ is fixed (and determined by external metabolite concentrations), the shape of the profile can vary. In the optimal profile (in red), small driving forces are avoided. The driving forces determine the ratios of forward and backward one-way fluxes (red arrows), and at a given net flux (black arrows) the enzyme demands. In the suboptimal blue curve, in contrast, the last three reactions show lower forces, and therefore relatively high reverse fluxes (blue arrows); to obtain the same net flux, forward and backward fluxes have to be strongly increased, which increases the enzyme demand.

must grow enormously, which would require a large amount of enzyme. This effect concerns only very small $\theta$ values - for $\theta$ much larger than 1 (or $\Delta G'$ much smaller than -RT), it can be neglected. Therefore, redistributing driving forces between reactions, to avoid very small forces, can save enzyme costs. The relation between driving forces, enzyme efficiency enzyme demand is shown in more detail in Figure 6.4.

If small driving forces should be avoided (to prevent enzyme costs from going infinity), how can this happen in practice? The driving forces themselves depend on metabolite levels, which can vary over several orders of magnitude. While the true metabolite concentrations are usually unknown, we hypothesize that selection favors concentration profiles that prohibit very small driving forces, in order to escape the ensuing large enzyme demands. Of course, completely avoiding small driving forces may be impossible, as there is always a trade-off: if a metabolite concentration decreases, the driving forces of all reactions producing it will increase, but the driving forces of all reactions consuming it will decrease simultaneously. So, all else being equal, the optimal metabolite profile is one that distributes its driving forces as evenly as possible.

### 6.3.3. Max-min Driving Force method

Previously in Chapter 4.3.2, we discussed adding thermodynamic constraints to constraint-based models in order to comply with the second law of thermodynamics. We can extend that approach in order to implement the idea of avoiding small driving forces. When we talk about the thermodynamic profile of a metabolic pathway, we usually try to visualize it by the cumulative Gibbs energy of reaction: we start at 0 and at each step add the $\Delta G'$ of the next reaction, which, assuming the pathway is feasible, is a negative number. The profile therefore has a shape of a downhill slope. The end point represents the total Gibbs energy and depends only on the concentrations of the metabolites that are part of the net reaction. Intermediate metabolites do not affect it, but they do determine the shape of the profile itself (see Figure 6.4). Specifically, each intermediate metabolite typically affects the driving force of two reactions – the one producing it and the one consuming it – with opposite signs. Therefore, changing the concentration of an intermediate can help increase the driving force of one reaction, but always at the expense of another reaction. This strong coupling
<figure><img src="image_58.png" /><figcaption>Figure 5.2: We show a cartoon of the solution space of a metabolic network, the so-called flux cone, with respectively (A) zero, (B) one and (C) 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>

product molar mass and turnover number. A simplified technique to provide a proxy for these costs is parsimonious flux balance analysis (pFBA). Central to pFBA is the assumption that cells minimize their total enzyme usage. Here, an optimal objective value is calculated via standard FBA, and the sum of the gene-associated reaction fluxes is then minimized. pFBA significantly reduces flux variability compared to standard FBA, but still does not typically capture overflow metabolism [150].

### 5.4.1 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 $e_i$ is the i-th EFM and $\lambda_i \ge 0$ its coefficient in v. 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 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.10). 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 con-
<figure><img src="image_46.png" /><figcaption>Figure 6.5: Max-Min Driving force method (MDF): an optimality problem in metabolite space. (a) Example pathway with given equilibrium constants and fixed concentrations of the external metabolites X and Y. What are the most favorable concentrations of the internal metabolites A and B? Assuming that small driving forces should be avoided in all reactions, MDF determines the metabolite profile that optimizes a worst case: it maximizes the worst (that is, smallest) driving force among all three reactions. (b) Driving force in reaction 1, as a function of the logarithmic concentrations of A and B, called $\ln a$ and $\ln b$. Higher concentrations of A (the reaction product) lead to smaller driving forces. Above a critical value (where X and A are in equilibrium), the driving force becomes negative, and a forward flux is impossible (grey region). The concentration of B, which does not participate in the reaction, does not play a role. (c) Driving force for reaction 2. Here, it is the ratio $b/a$ that counts. The lower the ratio (lower right), the higher the driving force. If the ratio is higher than the equilibrium constant, the driving force becomes negative (grey region). (d) Driving force for reaction 3. (e) By overlaying the contours in (b), (c), and (d) and taking the minimum value, we obtain the minimal driving force $\theta^{\min}$ among all three reactions. $\theta^{\min}$ is a piecewise linear function of $\ln a$ and $\ln b$ within the feasible range, yielding positive forces in all three reactions. The maximum point of this function is the optimum metabolite profile predicted by MDF. In the example shown, the feasible concentration space is entirely defined by the driving forces themselves, given the external concentrations. In general, physiological concentration ranges for all metabolites could further decrease the solution space and shift the optimum point (not shown).</figcaption></figure>

between $\Delta_r G'$ is why it is not trivial to find the optimal thermodynamic profile of a pathway.

The Max-Min driving force method (MDF) [226] is a method for predicting metabolite concentrations, based on the principle of evenly distributed driving forces. All fluxes are fixed and given, and assumed to be positive. It assumes that each metabolite concentration must remain in a predefined range, converts each choice of metabolite concentrations into the corresponding pattern of driving forces, and determines the smallest resulting driving force in the network. If this smallest driving force is negative, the flux distribution cannot be realized thermodynamically. Otherwise, the larger this smallest driving force, the better the overall metabolite profile. Hence, among all possible metabolite profiles, MDF predicts the one that maximizes the value of the minimal driving force across the network. Mathematically, this leads to a linear optimization problem: in the space of logarithmic metabolite concentrations, a lower bound on all driving forces (denoted B) is maximized (Eq. 6.10). An illustrative example is shown in Figure 6.5.

$$
\begin{array}{ll}
\text{Maximize}_{x,B} & B \\
\text{Subject to} & -(\Delta_r G' + RT \cdot N^\top x) \geq B \\
& \ln(C_{\min}) \leq x \leq \ln(C_{\max})
\end{array}
\qquad (6.10)
$$

MDF is easy to apply: it is based on a simple Linear Programming problem and requires only the following input data: (i) the stoichiometric network; (ii) the flux directions; (iii)) the known equilibrium constants (or equivalently, the standard reaction Gibbs free energies); (iv) physiological ranges for metabolite concentrations. Based on these data
straint. In the previous chapter, inhomogeneous constraint were written in the general form

$$ \sum_{i} \omega_{i}^{p} v_{i} \leq h_{p}, \quad p = 1, \dots, P \qquad (5.13) $$

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[151, 152]: only a limited number of macromolecules (proteins, ribosomes, etc) fit inside a cell. Since these molecules catalyze reactions, reaction rates are proportional to their concentrations:

$$ v_i = e_i k_i(\mathbf{s}), \qquad (5.14) $$

where $e_i$ is the concentration of the enzyme that catalyzes reaction $i$, and $k_i(\mathbf{s})$ 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} \omega_{i} e_{i} \leq 1, \qquad (5.15) $$

where $\omega_i$ are weights that determine how much of the resources are taken up by one unit of the ith enzyme; these 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{\omega_{i}}{k_{i}} v_{i} \leq 1, \qquad (5.16) $$

such that these resource-allocation constraints again fit the form presented in Equation (5.13). A well-known example of a modeling framework that uses such a constraint is FBA with macromolecular crowding (FBAwMC, [153]) 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 [154]. Moreover, we may have multiple inhomogenous constraints on flux values as Equation (5.13) 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 [154]. 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 [155], 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.5 Multi-objective flux analysis and flux variability

## 5.5.1 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
alone, metabolite concentrations and driving forces (or $\Delta G'$ values) are predicted. An example application can be found in Hädicke et al. [227], where the potential of CO₂ fixation in E. coli via endogenous pathways was analyzed using MDF.

A theoretical insight from MDF is the notion of distributed bottlenecks. A simple bottleneck would consist of a single reaction whose driving force cannot be increased because the substrates are at their upper concentration bounds and the products are at their lower concentration bounds. Given the fixed equilibrium constant, nothing can be done to increase the driving force in this reaction. A distributed bottleneck is more complicated: it consists of a series of reactions that all share the same low driving force, which, because of all the concentration constraints in the system, cannot be further increased (e.g. as in Figure 6.4). Even though each single reaction looks “harmless” because its own driving force could still be increased, this increase would happen at the expense of other driving forces.

### 6.3.4. The roles of thermodynamics for metabolic states

In summary, thermodynamics provides important clues both about the feasibility of pathways fluxes and about their enzyme demand. To use this knowledge, fluxes need to be considered together with metabolite concentrations (to obtain the possible driving forces), but no detailed knowledge of enzyme kinetics is required. Thermodynamics alone yields an upper bound on fluxes (and hence, a lower bound on enzyme demands) that holds for any kinetic rate laws. The only required data (except for the metabolic network itself) are equilibrium constants (or equivalently, standard Gibbs free energies of reactions $\Delta G'^{\circ}$), which can be obtained from *equilibrator.weizmann.ac.il* [228] as well as physiological bounds on metabolite concentrations. Given this information, and given a feasible choice of metabolite concentrations, we can compute the driving forces of all reactions, and from the factorized rate law (and assuming positive fluxes by convention) we can then approximate the reaction rates by $v \le E k_{cat} (1 - e^{-\theta})$.

We also remember from Chapter 3 that driving forces are not independent between reactions, but depend on the metabolite concentrations, which creates trade-offs: in a chain $A \xrightarrow{R_1} B \xrightarrow{R_2} C$, a lower concentration of B will increase the driving force in $R_1$, but decrease the driving force in $R_2$. For high enzyme efficiency (low enzyme demand), all driving forces should in principle be high, but this is most important for low $\theta$ values (while for $\theta \gg 1$ it does not even matter). Therefore we may conclude that, to save enzyme, a cell should rearrange its metabolite levels within physiological bounds such that small $\theta$ are avoided. Implementing this as an optimality problem, we obtain MDF.

In conclusion, we described (i) a general rule of thumb that poor thermodynamics makes reactions costly; (ii) simple approximations of enzyme cost; and (iii) practical methods (MDF) to obtain metabolite profiles with favorable thermodynamic properties.

## 6.4. Enzyme cost minimization

### 6.4.1. Enzyme cost minimization

The problem of minimizing the total enzyme demand (or cost) for a given pathway can be solved numerically, thanks to the fact that they are always convex [229]. Finding the minimum of the convex objective (the total enzyme cost) in a convex set (the set of admissible metabolite profiles, a convex polytope in log-metabolite space) can be done efficiently. In contrast to general optimality problems, such problems have a unique local optimum, which can be found by simple numerical methods. In this section, we demonstrate it with a simple example, the same three-reaction pathway that you already saw in Section 6.3 above.

### 6.4.2. Enzyme cost landscape of a metabolic pathway

Given the fluxes, kinetics, and concentration bounds in a metabolic pathway model, one can predict the enzyme demand by assuming that cells minimize the enzyme cost in that pathway. In the Enzyme Cost Minimization method A reaction
<figure><img src="image_59.png" /><figcaption>Figure 5.3: Phenotypic phase plane of the metabolic model shown in Figure 5.1, calculated as a function of the uptake rates of oxygen and glucose.</figcaption></figure>

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, an 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 $n_X = 10$, $n_1 = 1$, $n_3 = 4$, $v_4 = 0$, $v_0^{ub} = 10$, and $v_{O_2}^{ub} = 15$ mmol g⁻¹ h⁻¹ (cell dry mass).

At zero oxygen uptake, Figure 5.3.A shows that growth is possible reaching a specific growth rate of 1 1/h 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 mmol g⁻¹ h⁻¹ (cell dry mass) is 0.1 g mmol⁻¹ (cell dry mass). 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.17) $$

where z is the objective function optimal value (specific growth rate in this case) and $b_i^v$ corresponds to the violation of a mass balance constraint and is equivalent to the uptake reaction of the i-th metabolite (glucose in this example)[156]. 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 1/h. 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
rate $v = e \cdot f(c)$ depends on enzyme level $e$ and metabolite concentrations $c_i$ through the enzymatic rate law, $f(c)$. If the metabolite concentrations were known, we could directly compute enzyme demands $e = v/f(c)$ from fluxes, and similarly calculate the flux-specific enzyme demand $e/v = 1/f(c)$. However, metabolite concentrations are usually unknown and vary between experimental conditions. Therefore, there can be many solutions for $e$ and $c$ realizing one flux distribution. To select one of them, we employ an optimality principle: we define an enzyme cost function (for instance, total enzyme mass) and choose the enzyme profile with the lowest cost while restricting the metabolite levels to physiological ranges and imposing some thermodynamic constraints. As we shall see below, the solution is in many cases unique.

Let us demonstrate this procedure with a simple example (Figure 6.6 (a)). In the pathway $X \rightleftharpoons A \rightleftharpoons B \rightleftharpoons Y$, the external metabolite levels [X] and [Y] are fixed and given, while the intermediate levels [A] and [B] need to be found. As rate laws for each of the three reactions, we use reversible Michaelis-Menten (MM) kinetics

$$v = E \frac{k_{\text{cat}}^{+} s / K_{\text{S}} - k_{\text{cat}}^{-} p / K_{\text{P}}}{1 + s / K_{\text{S}} + p / K_{\text{P}}} \qquad (6.11)$$

with enzyme level $E$, substrate and product levels $s$ and $p$, turnover rates $k_{\text{cat}}^{+}$ and $k_{\text{cat}}^{-}$, and Michaelis constants $K_{\text{S}}$ and $K_{\text{P}}$. In kinetic modeling, steady-state concentrations would usually be obtained from given enzyme levels and initial conditions through numerical integration. Here, instead, we fix a desired pathway flux $v$ and compute the enzyme demand as a function of metabolite concentrations:

$$E(s, p, v) = v \frac{1 + s / K_{\text{S}} + p / K_{\text{P}}}{k_{\text{cat}}^{+} s / K_{\text{S}} - k_{\text{cat}}^{-} p / K_{\text{P}}}. \qquad (6.12)$$

Figure 6.6 shows how the enzyme demand in each reaction depends on the logarithmic reactant concentrations. To obtain a positive flux, substrate levels $s$ and product levels $p$ must be restricted: for instance, to allow for a positive flux in reaction 2, the rate law numerator $k_{\text{cat}}^{+} [A] / K_{\text{S}} - k_{\text{cat}}^{-} [B] / K_{\text{P}}$ must be positive. This implies that $[B]/[A] < K_{\text{eq}}$ where the reaction's equilibrium constant $K_{\text{eq}}$ is determined by the Haldane relationship, $K_{\text{eq}} = (k_{\text{cat}}^{+} / k_{\text{cat}}^{-}) \cdot (K_{\text{P}} / K_{\text{S}})$. With all model parameters set to 1, we obtain the constraint $[B]/[A] < 1$, i.e. $\ln[B] - \ln[A] < 0$, putting a straight boundary on the feasible region (Figure 6.6 (c)). Close to chemical equilibrium ($[B]/[A] \approx K_{\text{eq}}$), the enzyme demand $e_2$ approaches infinity. Beyond that ratio ($[B]/[A] > K_{\text{eq}}$) no positive flux can be achieved (grey region). Such a threshold exists for each reaction (see Figure 6.6 (b)-(d)). The remaining feasible metabolite profiles form a triangle in log-concentration space, which we call *metabolite polytope* $\mathcal{P}$ (Figure 6.6 (e)), and Eq. (6.12) yields the total enzyme demand $E_{\text{tot}} = E_1 + E_2 + E_3$, as a function on the metabolite polytope. The demand increases steeply towards the edges and becomes minimal in the center. The minimum point marks the optimal metabolite profile, and via Eq. (6.12) we obtain the resulting optimal enzyme profile.

The metabolite polytope and the large enzyme demand at its boundaries follow directly from thermodynamics. To see this, we consider the unitless *thermodynamic driving force* $\Theta = -\Delta_{\text{r}}G'/RT$ [88] derived from the reaction Gibbs free energy $\Delta_{\text{r}}G'$. The thermodynamic force can be written as $\Theta = \ln \frac{K_{\text{eq}}}{[B]/[A]}$, i.e. the driving force is positive whenever $[B]/[A]$ is smaller than $K_{\text{eq}}$, and it vanishes if $[B]/[A] = K_{\text{eq}}$. How is this force related to enzyme cost? A reaction's net flux is given by the difference $v = v^{+} - v^{-}$ of forward and backward fluxes, and the ratio $v^{+}/v^{-}$ depends on the driving force as $v^{+}/v^{-} = e^{\Theta}$. Thus, only a fraction $v/v^{+} = 1 - e^{-\Theta}$ of the forward flux acts as a net flux, while the remaining forward flux is partially canceled by the backward flux. Close to chemical equilibrium, where the mass-action ratio $[B]/[A]$ approaches the equilibrium constant $K_{\text{eq}}$, the driving force goes to zero, the reaction's backward flux increases, and the flux per enzyme level drops. This is what happens at the triangle edges in Figure 6.6: a reaction approaches chemical equilibrium, the driving force $\Theta$ goes to zero, and large enzyme amounts are needed for compensation. Exactly on the edge, the driving force vanishes and no enzyme level, no matter how large, can support a positive flux. The quantitative cost depends on model parameters: for example, by lowering a $k_{\text{cat}}$ value, the increase in enzyme cost at the boundary becomes steeper and the optimum point is shifted away from the boundary (see Figure 6.6 (f)).
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 [156]. 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.2 Non-uniqueness of optimal metabolic states: possible reasons

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 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.3 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 \text{), such that:} \\
\displaystyle\sum_{j \in J} N_{i,j} v_j = 0, \quad \forall i \in I, \\
v_j^{lb} \le v_j \le v_j^{ub}, \\
v_{\text{biomass}} = f \cdot v_{\text{biomass}}^*, \\
v_j \in \mathbb{R}, \quad \forall j \in J.
\end{array} 
\qquad (5.18) 
$$

Hence, 2N optimization problems need to be solved if there are N unconstrained fluxes. The results of the FVA 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 FVA. Instead, FVA inscribes this polytope in the smallest possible “box” [157]. 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 FVA) 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 FVA, could help in the identification of dead-end metabolites, and
<figure><img src="image_47.png" /><figcaption>Figure 6.6: Enzyme demand in a metabolic pathway. (a) Pathway with reversible Michaelis-Menten kinetics (equilibrium constants, catalytic constants, and $K_M$ values are set to values of 1, [A] and [B] denote the variable concentrations of intermediates A and B in mM). The external metabolite concentrations [X] and [Y] are fixed. Plots (b)-(d) show the enzyme demand of reactions 1, 2, and 3 at given flux $v = 1$ according to Eq. (6.12). Grey regions represent infeasible metabolite profiles. At the edges of the feasible region (where A and B are close to chemical equilibrium), the thermodynamic driving force goes to zero. Since small forces must be compensated by high enzyme levels, edges of the feasible region are always dark blue. For example, in reaction 1 (panel (b)), enzyme demand increases with the level of A (x-axis) and goes to infinity as the mass-action ratio [A]/[X] approaches the equilibrium constant (where the driving force vanishes). (e) Total enzyme demand, obtained by summing all enzyme levels. The metabolite polytope – the intersection of feasible regions for all reactions – is a triangle, and enzyme demand is a cup-shaped function on this triangle. The minimum point defines the optimal metabolite concentrations and optimal enzyme levels. (f) As the $k_{cat}$ value of the first reaction is lowered by a factor of 5, states close to the triangle edge of reaction 1 become more expensive and the optimum point is shifted away from the edge. (g) The same model with a physiological upper bound on the concentration [A]. The bound defines a new triangle edge. Since this edge is not caused by thermodynamics, it can contain an optimum point, in which driving forces are far from zero and enzyme costs are kept low. Please note the resemblance to the MDF problem for the same pathway, shown in Figure 6.5.</figcaption></figure>

### 6.4.3. Enzyme cost as a function of metabolite profiles

The prediction of optimal metabolite and enzyme levels can be extended to models with general rate laws and complex network structures. In general, enzyme demand depends not only on driving forces and $k_{cat}$ values, but also on the kinetic rate law, which includes $K_M$ values and small-molecule regulation. We can conveniently model or approximate these factors by using factorized rate laws. Let us write this rate laws here again in a general form to see the different factors at play. As we learned in Section 6.2, the rate of a reaction depends on enzyme level $e$, forward catalytic constant $k_{cat}^+$ (i.e. the maximal possible forward rate per unit of enzyme, in s⁻¹), driving force (i.e. the ratio of forward and backward fluxes), and on kinetic effects such as substrate saturation or small-molecule regulation. If all active fluxes are positive, reversible rate laws like the Michaelis-Menten kinetics in Eq. (6.11) can be factorized as [102]:

$$v = E \cdot k_{cat}^{+} \cdot \eta^{\text{for}} \cdot \eta^{\text{sat}} \cdot \eta^{\text{reg}}. \qquad (6.13)$$

Negative fluxes, which would complicate our formulae, can be avoided by orienting the reactions in the direction of fluxes.

Enzyme demand can be quantified as a concentration (e.g. enzyme molecules per volume) or mass concentration (where enzyme molecules are weighted by their molecular weights). If rate laws, fluxes, and metabolite concentrations are known, the enzyme demand of a single reaction $l$ follows from Eq. (6.13) as

$$E_l(\mathbf{c}, v_l) = v_l \cdot \frac{1}{k_{cat,l}^+} \cdot \frac{1}{\eta_l^{\text{for}}(\Theta(\mathbf{c}))} \cdot \frac{1}{\eta_l^{\text{sat}}(\mathbf{c})} \cdot \frac{1}{\eta_l^{\text{reg}}(\mathbf{c})}. \qquad (6.14)$$
<figure><img src="image_60.png" /></figure>

Philosophical remark 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 actually is. The problem is that, since we have many different parameters to choose from, many different models can explain the same metabolic observations [158].

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

in the long run, in model improvement.

5.6 Concluding remarks

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

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
To determine the enzyme demand of an entire pathway, we sum over all reactions: $E_{tot}^{path} = \sum_{l} e_{l}$. Based on its enzyme demands $E_{l}$, we can associate each metabolic flux with an enzyme cost $q = \sum_{l} h_{E_{l}} E_{l}$, describing the effort of maintaining the enzymes. The burdens $h_{E_{l}}$ of different enzymes represent, e.g. differences in molecular mass, post-translational modifications, enzyme maintenance, overhead costs for ribosomes, as well as effects of misfolding and non-specific catalysis. The enzyme burdens $h_{E_{l}}$ can be chosen heuristically, for instance, depending on enzyme sizes, amino acid composition, and lifetimes. Setting $h_{E_{l}} = m_{l}$ (protein mass in grams per mole), $q$ will be in gram protein per gram cell dry weight. Considering the specific amino acid composition of enzymes, we can also assign specific costs to the different amino acids. Alternatively, an empirical cost per protein amount can be established by the level of growth impairment that an artificial induction of protein would cause [60, 230]. Thus, each reaction flux $v_{l}$ is associated with an enzyme cost $q_{l}(v_{l}, c) \equiv h_{E_{l}} E_{l}(c, v_{l})$ of flux and metabolite concentrations. From now on, we refer to log-scale metabolite concentrations $s_{i} = \ln c_{i}$ to obtain simple optimality problems below. From the factorized rate law Eq. (6.14), we obtain the enzyme cost function

$$q(\mathbf{s}, \mathbf{v}) \equiv \sum_{l} h_{E_{l}} E_{l}(v_{l}, \mathbf{s}) = \sum_{l} h_{E_{l}} \cdot v_{l} \cdot \frac{1}{k_{\text{cat},l}^{+}} \cdot \frac{1}{\eta_{l}^{\text{for}}(\mathbf{s})} \cdot \frac{1}{\eta_{l}^{\text{sat}}(\mathbf{s})} \cdot \frac{1}{\eta^{\text{reg}}(\mathbf{s})} \quad (6.15)$$

for a given pathway flux $v$. If the fluxes are fixed and given, our enzyme cost becomes, at least formally, a function of the metabolite levels. The cost function is defined on the metabolite polytope $P$, a convex polytope in log-concentration space containing the feasible metabolite profiles. Like the triangle in Figure 6.6, the polytope is defined by physiological and thermodynamic constraints.

Beyond minimizing the total enzyme cost, one can also use Enzyme Cost Minimization to analyze the individual enzyme demands. When the metabolite levels are known, the demand can be directly calculated and each efficiency factor ($\eta$) in Eq. (6.15). By omitting some factors or replacing them by constant numbers $0 < \eta \le 1$, simplified enzyme cost functions with fewer parameters can be obtained. For example, $\eta^{\text{for}} = 1$ would imply an infinite driving force $\Theta \to \infty$ and a vanishing backward flux, $\eta^{\text{sat}} = 1$ implies full substrate saturation, and $\eta^{\text{reg}} = 1$ implies full enzyme activation and no enzyme inhibition (or no small-molecule regulation at all). In these limiting cases, enzyme activity will not be reduced, and enzyme demand will be given by the capacity-based estimate $v/k_{\text{cat}}^{+}$, a lower estimate of the actual demand. Instead of omitting an efficiency factor, it can also be set to a constant value between 0 and 1. Such simplifications and the resulting enzyme cost functions with fewer parameters can be practical if kinetic constants are unknown.

### 6.4.4. General lessons from Enzyme Cost Minimization

Enzyme cost minimization not only provides numerical solutions, but also some general insights.

* **Convexity** Enzyme Cost Minimization shows again the importance of the metabolite polytope. The usage of logarithmic metabolite concentrations not only leads to a good search space for feasible metabolite profiles (as in MDF), but also facilitates optimization because enzyme cost is a convex function of the metabolite log-concentrations [231]. Convexity makes this optimization tractable and scalable – unlike a direct optimization in enzyme space. Convexity holds for a wide range of rate laws and for extended versions of the problem, e.g. including bounds on the sum of (non-logarithmic) metabolite concentrations or bounds on weighted sums of enzyme fractions.

* **Factorized rate laws disentangle individual enzyme cost effects** To see how metabolic states are shaped by different physical factors, we considered factorized rate laws. The different terms in these functions represent specific physical factors and require different kinetic and thermodynamic data for their calculation. By neglecting some terms, one obtains different approximations of the true enzyme cost. By comparing the different scores, we can estimate the enzyme cost that cells “pay” for running reactions at small driving forces (to save Gibbs free energy) or for keeping enzymes beneath substrate-saturation (e.g., to dampen fluctuations in metabolite levels).

* **Relationship to other optimality approaches** Beyond their practical advantages, factorized enzyme cost functions also allow us to easily compare our method to earlier modeling and optimization approaches. These approaches typically focused on only one or two of the factors that are taken into account in Enzyme Cost Minimization, and
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.

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.

## Recommended readings and tools

Escher FBA Escher FBA (sbrg.github.io/escher-fba/) is a nicely illustrative tool for 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

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

### Problem 5.1 Flux distribution with constraint

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

(a) Using a spreadsheet and its associated linear programming optimizer.

(b) Using an LP solver in Python such as `linprog` available in `scipy/optimize`.

(c) 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 Choice between phenotypes

The metabolic network illustrated in Figure 5.1, adapted from [156], was designed to include four pheno-
many of them can be reformulated as approximations of this method [226, 232, 220].

* **Enzyme cost is related to thermodynamics** In FBA, thermodynamic constraints and flux costs appear as completely unrelated aspects of metabolism (as is explained in Chapter 5). Thermodynamics is used to restrict flux directions, and to relate them to metabolite bounds, while flux costs are used to suppress unnecessary fluxes. In Enzyme Cost Minimization, thermodynamics and flux cost appear as two sides of the same coin. Like in FBA, flux profiles are thermodynamically feasible if they lead to a non-empty metabolite polytope, allowing for positive forces in all reactions. However, the values of these forces also play a role in shaping the enzyme cost function on that polytope. Together, metabolite polytope and enzyme cost function (as in Figure 6.6) summarize all relevant information about flux cost.

Many pathways are regulated, for instance by feedback inhibition of enzymes via the end product. While this may stabilise the dynamics and adapt it to current demands, such enzyme regulation comes at a cost, which we can estimate by following the logic of Enzyme Cost Minimization. Many enzymes are regulated by small molecules that act as competitive or allosteric inhibitors [233], an effective way to implement feedback control, for example to adapt the flux in biosynthesis pathways to current needs. In order for such a regulation to work, the enzyme needs to be partially inhibited on average (because only then, its activity can be increased on demand, by alleviating the inhibition). Therefore, the enzyme efficiency goes down, and the cell needs to provide more enzyme to catalyze the same flux than without the inhibition. How much will this regulation cost the cell as part of the enzyme budget? From the perspective of Enzyme Cost Minimization, where we start from desired fluxes and compute the enzyme demand, this question is easy to answer: in the inhibited enzyme case, the lower efficiency will be described by a factor $\eta^{\text{reg}} \in [0, 1]$ (Mathematical Details Box 6.C). In the same reaction, the enzyme demand increases by a factor $1/\eta^{\text{reg}}$, so the extra cost is simply $1/\eta^{\text{reg}} - 1$ times the “baseline” cost of this enzyme (without inhibition). Specifically, a non-competitive inhibitor, with efficiency factor $\eta^{\text{reg}} = \frac{1}{1+c/K_I}$ yields a cost factor $1+c/K_I$. If the metabolite concentrations are fixed, this corresponds to an extra enzyme demand $\Delta E_l = \frac{E_l c_i}{K_{I,li}}$. Similarly, an enzyme activation with efficiency factor $\eta^{\text{reg}} = \frac{c/K_A}{1+c/K_A}$ in the rate laws yields a cost factor $\frac{1+c/K_A}{c/K_A} = 1 + K_A/c$ in the formulae for enzyme demands. If the metabolite concentrations are fixed, this corresponds to an extra enzyme demand $\Delta E_l = \frac{E_l K_{A,li}}{c_i}$ (where $l$ and $i$ denote the regulated reaction and the regulating metabolite, respectively). As usually in Enzyme Cost Minimization, an optimal rearrangement of enzyme and metabolite concentrations must be taken into account, which will then slightly reduce the overall cost.

The predictions of optimal states by Enzyme Cost Minimization rely on two main inputs: a metabolic model that relates metabolite concentrations, enzyme levels, and fluxes, and an optimality principle based on the assumption that cells realize their production fluxes at a minimal total enzyme cost. To test whether this optimality principle holds at all, Noor et al. [229] compared the predictions from Enzyme Cost Minimization to predictions from the same metabolic model and the same flux distribution, but with randomly sampled metabolite profiles (and the corresponding enzyme profiles). In comparison, metabolite profiles sampled close to the Enzyme Cost Minimization optimum yielded significantly better enzyme level predictions than metabolite profiles sampled more broadly. This strongly supports the idea that *E. coli* metabolism, in the conditions studied, is at least partially optimized for low enzyme cost, and thus supports cost-optimality as a principle in living cells.

## 6.5. Comparison of alternative pathways

Having clarified our main functional criteria for pathways (substrate productivity and enzyme productivity) and how they depend on pathway details (including outer concentrations), we can now compare alternative pathways by their substrate and enzyme demand per production flux (an example of “cost per benefit”) and see which one scores better.

### 6.5.1. A tale of two glycolyses

One of the canonical examples discussed throughout this book is how cells choose between respiration and fermentation for making their ATP. However, having a precise kinetic model for respiration is difficult, since it involves electron transfer
<table><caption>Table 5.1: Stoichiometry of the metabolic network for problem 5.2. Adapted from [160] after [156].</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td><i>qA</i> &rarr; A</td><td>ATP <i>R</i><sub>ft</sub></td></tr><tr><td>A + ATP <i>R</i><sub>1</sub> &rarr; B</td><td>C + 10 ATP <i>R</i><sub>z</sub> &rarr; Biomass</td></tr><tr><td>B <i>R</i><sub>2</sub> &rarr; 2 ATP + 3 NADH + C</td><td><i>qO</i><sub>2</sub> &rarr; O<sub>2</sub></td></tr><tr><td>0.2 C <i>R</i><sub>3</sub> &rarr; 2 NADH</td><td>C <i>C</i><sub>out</sub></td></tr><tr><td>C <i>R</i><sub>4</sub> &rarr; ATP + 3 D</td><td>D <i>D</i><sub>out</sub></td></tr><tr><td>C + 2 NADH <i>R</i><sub>5</sub> &rarr; 3 E</td><td>E <i>E</i><sub>out</sub></td></tr><tr><td>NADH + O<sub>2</sub> <i>R</i><sub>Res</sub> &rarr; 2 ATP</td><td></td></tr></tbody></table>

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

(a) If the uptake of the carbon source A is bounded between 0 and 10 mmol g⁻¹ h⁻¹ and no restrictions on the oxygen uptake are imposed, prepare a plot showing the biomass, C, D and E fluxes attained at different uptakes of A.

(b) Repeat the preceding analysis, but limit the maximum uptake rate of oxygen to 10 mmol g⁻¹ h⁻¹.

(c) If substrates uptakes are bounded between 0 and 10 mmol g⁻¹ h⁻¹ for A and 0 and 20 mmol g⁻¹ 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.
and membrane-bound reactions. Therefore, it is challenging to calculate the enzyme cost of respiration using models like those discussed in this chapter. Flamholz et al. [225] analyzed a similar but simpler case by comparing between the EMP and ED variants of glycolysis, since all the required enzymes are soluble and expressed in the cytoplasm and/or the periplasm and many of their kinetic parameters are measured. The common description of glycolysis ends in pyruvate (e.g., as depicted in Figure 6.1). This means that the pathway is not neutral in terms of redox, since the oxidation state of pyruvate is higher than glucose. In order to simplify the comparison and focus only on ATP yield (rather than NADH), the EMP and ED pathways were extended to end in lactate by including lactate dehydrogenase (ldh) as an extra step, making them redox neutral. These could be thought of as the more relevant versions of the pathways in anaerobic conditions.

Although EMP-based fermentation is usually described in textbooks as less efficient than respiration, since it produces only 2 moles of ATP per mole glucose instead of ≈ 30, the ED pathway has an even lower yield – 1 mole of ATP. Nevertheless, the ED pathway is quite common among the bacteria. For example, *Zymomonas mobilis* – the bacterium used in fermenting pulque (a.k.a., agave wine [234]) and a promising platform for bio-production [235] – lacks key enzymes from the EMP pathway and uses the ED pathway exclusively to metabolize sugars. These bacteria don't seem to be bothered by the low ATP yield and can achieve high growth rates [236]. This already suggests to us that the ED pathway is probably superior to EMP in other aspects, such as the enzyme demand. Another clue was provided by a study which found that the ED pathway improves *E. coli* growth during glucose up-shifts and that the flux through it increases by 130% [237] (see Economic Analogy Box 6.D).

To see if indeed the models provide predictions that are consistent with the experimental evidence, Flamholz et al. [225] first used the MDF method to compare the two pathways. The ED pathway was found to be substantially more thermodynamically favorable, with a much higher score than the EMP pathway (8.0 versus 4.8 kJ/mol, see Figure 6.7 upper row).

Although the EMP pathway is clearly more favorable, we can still argue that an MDF of 4.8 kJ/mol is good enough, as it means $\theta > 1.9$ for each one of the pathway reactions. In this case, $\eta^{\text{for}} > 0.85$ (see Figure 6.3) and therefore it might be a small price to pay for double the ATP yield. But, as discussed earlier, the efficiency of a pathway is affected by other factors besides the thermodynamics. Flamholz et al. [225] tried to see whether ED is superior to EMP also in terms of the enzyme cost using the Enzyme Cost Minimization method. Indeed, they found that the ED pathway would require ≈5 times less protein compared to EMP for catalyzing the same flux (see Figure 6.7 bottom row). So, although the ATP yield of the ED pathways is half that of EMP, one can still generate ATP at a higher rate using the same amount of protein, according to the model.

The comparison of EMP and ED provided some insight as to a trade-off that can exists between the yield of a pathway and its cost, or enzyme burden. However, one can expand the question and ask if there are any other theoretically possible glycolysis pathways that might be able to break this trade-off and be more efficient than EMP and ED in both aspects. Ng et al. [214] tried to address this question with an algorithm they called *optStoic* that generates all biochemically feasible routes between glucose and pyruvate, with various ATP/glucose yields. They then ran pathway analysis on all 11,916 options and found that indeed both EMP and ED are both (nearly) Pareto-optimal. This suggests that evolution may indeed select for features such as high yield and low enzyme cost, where one might be more important than the other depending on the context.

### 6.5.2. Metabolic engineering

Besides the quest for understanding the evolution of existing biochemical pathways, pathway analysis methods like MDF and Enzyme Cost Minimization have also been used by metabolic engineers in order to rank and prioritize different alternative designs. For example, Volpers et al. [238] used the MDF algorithm and the *Pathway Specific Activity* measure to compare between designs of photo-electro-autotrophic strains. Similarly, Löwe and Kremling [239] used the Enzyme Cost Minimization algorithm to predict the enzyme demand of both natural and synthetic carbon fixation cycles.
# Chapter 6

# The enzyme cost of metabolic fluxes

Wolfram Liebermeister and Elad Noor

## Chapter overview

* 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 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 features that might appear strange at first glance [161]: 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
<figure><img src="image_48.png" /><figcaption>Figure 6.7: Comparing two metabolic pathways using the MDF and the Enzyme Cost Minimization methods. For the MDF analysis (top row), the dashed green line represents the cumulative Gibbs energy along the pathway if all metabolite concentrations were 1 mM. The MDF solution is presented as a grey line, where the bottleneck reactions are marked in red. For the Enzyme Cost Minimization analysis (bottom row), we used the same kinetic parameters for all enzymes in both pathways ($k_{cat}$ = 200 s⁻¹, K_M = 200 µM, same as in [225]). However, here we used an updated version of Enzyme Cost Minimization with the factorized rate law, therefore the results are not identical. The script for generating the figure can be found in the following Jupyter notebook.</figcaption></figure>

### 6.5.3. Predicting concentrations

So far, the examples given in this section focused on analyzing and comparing pathway alternatives in isolation, outside of the context of actual living organisms. However, we should not forget that the motivation for optimization goals such as enzyme demand are derived from physiological and evolutionary principles. Therefore, the optimal solutions coming from MDF and Enzyme Cost Minimization might be good predictions for the actual metabolic state that exists in naturally evolved organisms.

For example, a few years after the *in silico* analysis of the ED pathway [225], Jacobson et al. [240] measured the intracellular concentrations ED intermediates in *Z. mobilis*, and used them to calculate the Gibbs energies of the pathway's reactions. Indeed, they found that they closely fit the predicted values from the MDF solution. Similarly, measured values of enzyme and metabolite concentrations in *E. coli* correlate with predicted values from Enzyme Cost Minimization (when empricial reaction fluxes were obtained from ¹³C-MFA measurements, Figure 6.8) [229]. In a related paper, Wortel et al. [241] expanded the idea of this method to explore the entire flux polytope.

These results suggest that indeed the optimization process that occurs throughout evolution is somewhat similar to the (much simplified) models presented here. Of course, improving the accuracy of the inputs and accounting for other effects that impact fitness could improve the predictions further. On the other hand, it might be naïve to expect natural systems to be optimal, which would mean that using basic principles to *precisely* predict phenotypes is an impossible task.
<figure><img src="image_61.png" /></figure>

Philosophical remark 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.

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 [162] – 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.

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 gen-
eral term for pathways that convert glucose to pyruvate while producing ATP [161]. One variant of glycol-
ysis, named after Gustav Embden, Otto 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 [162]. 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). An-
<figure><img src="image_49.png" /><figcaption>Figure 6.8: Validation of metabolite and enzyme concentrations, predicted by Enzyme Costs Minimization, in *E. coli*'s central carbon metabolism.</figcaption></figure>

Economic analogy 6.D: The push for fast growth

The ED pathway seems to be useful as a quick response to a sudden increase in abundance of resources (glucose), but less efficient than EMP when the environment is steady. This is somewhat analogous to start-up companies, which burn large amounts of venture capital in order to grow rapidly. However, after reaching a certain scale, the dynamic nature of start-ups often becomes a burden, where overhead costs pile up and signal that it is time to join a larger corporation.

6.6. Can we now predict which pathways will be selected by evolution?

Coming back to our initial question, what have we learned from theory about the choice between possible pathways? The “choice between pathways” in a larger network is actually a choice between (network-wide) flux distributions that use different alternative pathways. Here we discussed how to score the usefulness of given flux distributions, which can also be used to score single pathways.

Importantly, flux distributions are scalable (by scaling all enzyme levels proportionally, and keeping all metabolite levels constant). If we scale the fluxes, this will scale both the flux benefit (for instance, the production of a desired product or biomass) and the required resources (substrates consumed, enzyme budget invested, or toxic byproducts produced). Because of this scaling property, our “quality criteria” mostly have the form of ratios between an output flux (as the benefit) and some (limited) resource (the cost). Such ratios are called “productivities”, where in Chapter 4-5 we focused mostly on substrate productivity (or yield on substrate) and in this chapter on enzyme productivity (or enzyme-specific rate) as important criteria. Why these criteria? On the one hand, they are closely related to some big objectives of the entire cell – depending on the type of competition it is facing. On the other hand, they are easy to link to some concrete criteria about metabolic pathways such as product yield, pathway length, $k_{cat}$ values, thermodynamic forces, etc.

Since yield on substrate depends only on the shape of the flux distribution, it can be studied by methods like FBA (see chapters 4 and 5). In this chapter, we focused on the more difficult case, enzyme productivity, where thermodynamics, enzyme kinetics, and the arrangement of metabolite and enzyme concentrations come into play. The factorized law in Eq. (6.7) shows us how the enzyme demand of a flux distribution can be computed if metabolite concentrations are known, and how the demand depends on forward $k_{cat}$, the thermodynamic force, and enzyme saturation. The only difficulty is that the thermodynamic forces and metabolite concentrations are usually not known. Here we considered some best-case scenarios, assuming that the cell will realize the concentration arrangements that optimize pathway performance. When considering thermodynamics alone (and making some further simplifications), this led to the MDF method. For the full problem, the solution is provided by Enzyme Costs Minimization. This method is directly related to
<figure><img src="image_62.png" /><figcaption>Figure 6.1: Two natural variants of the glycolysis pathway, named after their discoverers: Embden-Meyerhof-Parnas (EMP) and Entner-Doudoroff (ED) – The pathways are part of the core metabolism shown in Figure 3.1.</figcaption></figure>

other common variant was discovered in 1952 by Nathan Entner and Michael Doudoroff [163] (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 [93].

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. [164] 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. [165] 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 [166]) 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. [164] used this algorithm with the stoichiometry from 6.1 to find all possible glycolysis pathways comprising known enzymatic reactions (see Box 6.B).

The objective set by the linear problem (6.2) is minimizing the sum of fluxes, which corresponds to pathways with fewer reactions and low fluxes in each one. As discussed in 5.2, this objective is only a crude proxy for the efficiency of a pathway, and its only purpose is to get the pathway solutions in a relatively logical order. Although we have discussed global enzyme constraints in previous chapters (such as molecular crowding and proteome allocation), when comparing pathways we will focus only on the efficiency of the pathway itself. This will allow us to compare pathways without thinking about the rest of the cell or a specific metabolic context. But how can one quantify the efficiency of a pathway? The next section will be dedicated to exactly this question.
### Economic analogy 6.E: Two important assumptions: homogeneity and stationarity

In the models described in this chapter, we generally assume that our system (for example, a metabolic pathway in a cell) is spatially and temporally homogeneous, and that it shows stable stationary states. This is clearly a simplification: in reality, cells are inhomogeneous, with compartments, with enzymes unequally distributed across the cell, and with enzymes forming complexes or dedicated compartments like the glycosome (an organelle in some organisms that contains the glycolytic enzymes), which changes (average) enzyme kinetics. Cells are also dynamic on various time scales (chemical noise, metabolic dynamics, protein expression dynamics), which also may change (average) enzyme kinetics. If we ignore this in our models—assuming a timeless steady state—this will not only cause approximation errors in our metabolic model, but much more importantly, we ignore the fact that the cell can exploit spatial inhomogeneity (e.g. compartments or channeling) and non-steady states (e.g. metabolic oscillations, or adaptation to fluctuations in the environment) to further improve its fitness (as compared to a steady-state, constant enzyme model).

Interestingly, classical economic theory makes similar assumptions – e.g. about markets in equalibrium– which ignore the spatio-temporal, dynamic side of real economic systems, which – as in the case of metabolic models – is likely to lead to wrong results.

the different pathway criteria we discussed initially (including pathway length, thermodynamic forces, and $k_{cat}$ values) and thus shows how these different factors determine enzyme demand. As a numerical method, it is relatively easy to use because it is a convex optimization problem. But if little data is available, simpler methods such as MDF, with their lower demand for parameters, may be useful tools to predict pathway usage.

### Recommended readings

**A search for efficient pathways, based on different criteria:** Arren Bar-Even, Elad Noor, Nathan E. Lewis, and Ron Milo. Design and analysis of synthetic carbon fixation pathways. Proceedings of the National Academy of Sciences, 107(19):8889–8894, 2010. doi: 10.1073/pnas.0907176107.

**The max-min driving force method:** Elad Noor, Arren Bar-Even, Avi Flamholz, Ed Reznik, Wolfram Liebermeister, and Ron Milo. Pathway thermodynamics highlights kinetic obstacles in central metabolism. PLoS Comput. Biol., 10(2):e1003483, February 2014. doi: 10.1371/journal.pcbi.1003483.

**Enzyme cost minimization:** Elad Noor, Avi Flamholz, Arren Bar-Even, Dan Davidi, Ron Milo, and Wolfram Liebermeister. The protein cost of metabolic fluxes: Prediction from enzymatic rate laws and cost minimization. PLoS Comput. Biol., 12(11):e1005167, November 2016. doi: 10.1371/journal.pcbi.1005167.

### Problems

**Problem 6.1 Pathway efficiencies** Estimate pathway efficiencies (i.e. product production rates per total enzyme concentration) from simple back-of the envelope calculations and plausible numbers (refer to the BioNumbers database for realistic values). (a) From pathway length (assuming reasonable apparent $k_{cat}$ values); (b) from given apparent $k_{app}$ values (or given $k_{cat}$ values and $\Delta G$). (c) Convert the results into growth rates (assuming realistic estimates of the total protein density; the proteome fraction of metabolic enzymes; the biomass production rate etc). Assume plausible numbers in all cases.

**Problem 6.2 Efficiency – dependence on substrate** Compute the reduction of pathway efficiency in a linear chain when decreasing the external substrate concentration (no constraints on metabolite levels).

**Problem 6.3 ATP yield in glycolysis** Derive the optimal ATP yield in glycolysis under a linear flux-force relationship
<figure><img src="image_63.png" /></figure>

Math box 6.B Integer cuts for iterating all possible pathway variants

The linear problem can be described by:

$$
\begin{array}{ll}
\text{minimize} & \displaystyle\sum_{i} v_i \\
\text{subject to} & \mathbf{N}_{\text{int}} \mathbf{v} = 0 \qquad (6.2) \\
& \forall i \quad 0 \le v_i \le \beta \\
& v_{\text{glycolysis}} = -1
\end{array}
$$

where **v** is the flux variable, and **N** is comprised of the universal stoichiometric matrix , and in addition one
reaction (whose flux is denoted v<sub>glycolysis</sub>) which has the stoichiometry of Eq. (6.1). The constraint v<sub>glycolysis</sub> = -1
ensures that the sum of all active reactions except for v<sub>glycolysis</sub> will together form a full glycolysis pathway, since
their net reaction has to balance the stoichiometry of v<sub>glycolysis</sub> given the mass balance constraint N<sub>int</sub> v = 0. β
given the upper bound on the flux for all reactions. For simplicity, we assume that all fluxes are positive and
that reversible reactions are split into their two opposing directionalities . β is a tunable parameter that is an
upper bound on all the fluxes in the solution pathways. Setting it too low would exclude solutions with complex
stoichiometries. On the other hand, a very high value would increase the complexity of the search and lead to
very long run-times. Typically, we choose β = 10 which is a good balance between the two extremes. Finally,
we set the objective function (∑i v<sub>i</sub>) to minimize the sum of fluxes. As we will explain shortly, we can iterate
through all possible solutions and therefore the objective will only determine the order at which we find them.

To find all possible glycolysis pathways comprising known enzymatic reactions, Ng et al. [164] iteratively intro-
duced constraints in order to exclude all previous solutions and find the next optimal one [167]: to exclude a
solution, they add an integer cut, which is an inequality constraint ensuring that the number of active reactions
is strictly larger than the sum over their indicator variables (boolean variables that are equal to 1 if the reaction
is active, i.e. carries a nonzero flux). Therefore, at least one of those reactions must be inactive in all future so-
lutions. This is quite similar to constrained Minimal Cut Sets (cMCS) which were introduced in Chapter 4.5 as a
way of exploring the flux space.

Formally, if {P₀, P₁... Pₘ} are the set of solutions already discovered by our algorithm (where ∀j Pⱼ ⊆ {0, ..., n}, i.e. each solution is a set of integers which correspond to indices of active reactions) then the added constraints will be:

$$
\begin{align*}
& \forall i \quad z_i \in \{0, 1\} \\
& \forall i \quad v_i - \beta z_i \le 0 \\
& \forall j \quad \sum_{i \in P_j} z_i < \|P_j\| \tag{6.3}
\end{align*}
$$

where ||P_j|| is the length of pathway j (i.e. the number of reactions). The z_i are boolean reaction indicators,
i.e. z_i must be equal to 1 if a reaction is active (v_i > 0). The final set of constrains eliminate P_j and any path-
way which is a superset of P_j from the solution space. Using this extra set of constraints iteratively, each time
generating the next pathway and adding it to the excluded list, will eventually go through all possible solutions
(by increasing order of their sum of fluxes). It is important to note that using integer cuts requires switching to
an MILP (Mixed-Integer Linear Program) solver, which is computationally much more demanding and typically
requires a commercial license.

6.2 Pathway efficiency - some notions and thoughts

For glycolysis alone, Ng et al. [164] found 11,916 alternatives that produce at least one mole of ATP per mole of glucose. These include, of course, the EMP pathway. Although evolution can explore these options, natural selection typically converges on one or a few efficient variants. This does not mean that every single pathway observed in nature must be optimal, but we generally expect cells hosting highly inefficient pathways to eventually become extinct. Iacometti et al. [168] tested this experimentally by knocking out the EMP pathway from E. coli and forcing the cells to use the alternatives that naturally exist in this bacterium. In all cases, growth rates were slower than in the wild-type.

Before we discuss other examples for metabolic pathways, we need to define what we mean by “efficiency”.
There are several criteria one should consider:
**Problem 6.4 MDF method** Implement the MDF method in a programming language of your choice.

**Problem 6.5 MDF and enzyme cost** The optimality principle of MDF (avoiding small thermodynamic driving forces) can be justified by assuming that low driving forces would entail high enzyme demands. Do you expect that MDF solutions are also Enzyme Costs Minimization solutions (or vice versa)? Otherwise, can you think of an approximation of the Enzyme Costs Minimization problem, such that MDF provides the correct solution? Show how the Enzyme Costs Minimization objective could be approximated step by step, and illustrate this with an example.

**Problem 6.6 Cycle of chemical reactions** Assume a cycle of chemical reactions $A \leftrightarrow B \leftrightarrow C \leftrightarrow A$ without co-factors or external inputs/outputs. (a) Show that there is no stationary, thermodynamically feasible flux distribution except for the (trivial) vanishing flux. (b) Explain why, if there were a flux, this would be a perpetuum mobile.

**Problem 6.7 Optimal enzyme levels in two-reaction chain** Consider a chain of two reactions $S \leftrightarrow X \leftrightarrow P$ with enzymes $E_1$ and $E_2$, $v_1 = E_1(k_{+1}S - k_{-1}X)$, $v_2 = E_2(k_{+2}X - k_{-2}P)$. Compute the steady state flux given $E_1, E_2$. Let $E_1+E_2 = E_{tot}^{path}$ be fixed. Determine $E_1, E_2$ such that the flux is maximal. Use Lagrange multipliers. Hint: Assume forward flux where $P/S < (k_{+1}k_{+2})/(k_{-1}k_{-2}) = q_1q_2$.

**Problem 6.8 Flux maximization in a linear pathway** Prove that the function:

$$f(\vec{E}) = \frac{1}{\sum_{i}(A_{i}E_{i})^{-1}} \qquad (6.16)$$

for a fixed $\vec{A}$ and under the constraint $\sum_i E_i = E_{tot}$, is at its maximum when:

$$E_i = E_{tot} \cdot \frac{A_i^{-1/2}}{\sum_i A_i^{-1/2}}$$

**Problem 6.9 Haldane kinetic rate law** Haldane described an enzyme-catalyzed reaction by three steps, each following a mass-action rate law:

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

The ODE system describing the change in time of each species is:

$$\begin{aligned}
\frac{d[ES]}{dt} &= [E] \cdot [S] \cdot k_1 + [EP] \cdot k_4 - [ES] \cdot (k_2 + k_3) \\
\frac{d[EP]}{dt} &= [E] \cdot [P] \cdot k_6 + [ES] \cdot k_3 - [EP] \cdot (k_4 + k_5) \\
\frac{d[E]}{dt} &= -[E] \cdot [S] \cdot k_1 + [ES] \cdot k_2 + [EP] \cdot k_5 - [E] \cdot [P] \cdot k_6
\end{aligned}
\qquad (6.18)$$

Prove that at quasy-steady-state (where the total enzyme concentration is fixed, and the concentration of each species doesn't change over time), the rate in which $[S]$ is converted to $[P]$ is governed by the following rate law:

$$v = [E_0] \frac{k_{\text{cat}}^{+}[S]/K_S - k_{\text{cat}}^{-}[P]/K_P}{1 + [S]/K_S + [P]/K_P} \qquad (6.19)$$
* Low consumption rate of the substrate

* High generation rate of the product

* High regeneration rate or low consumption rate of the co-factor

* Small number of steps [169]

* Higher thermodynamic forces [170, 171]

* High enzyme turnover numbers

* High enzyme saturation levels

Some of these criteria refer to the cost (or investment) of the pathway, while others reflect the benefit (or profit) to the cell. By considering two common scenarios – single nutrient limitation or exponential growth in rich media – we can focus on two simple criteria which provide good measures of efficiency.

When the availability of a single nutrient is limiting growth, maximizing the molar yield (i.e. the number of moles of product generated for each mole of the nutrient) becomes the important feature. Yield is rather straightforward to calculate, as it is a direct outcome of the stoichiometry of the pathway. For example, anaerobic fermentation is often compared to respiration and deemed inefficient since it yields two moles of ATP per glucose, instead of ≈30 [172].

On the other hand, when conditions are good, such as during exponential growth in rich media, minimizing the total number of proteins required is often the objective which determines growth rate. Here, we will be using the enzyme demand (e.g. in grams of protein) per unit of flux (typically, in mmol per hour per gram of cell dry weight). In fact, the enzyme demand per flux, as an objective, takes into consideration both the cost (protein) and the benefit (flux). Importantly, these two criteria scale linearly with respect to each other: doubling the amount of all enzymes without changing any of the metabolite concentrations would directly double the flux in the pathway. Therefore, this measure of efficiency is independent of the magnitude of the flux in the pathway. But, as we will see shortly, enzyme demand is a non-linear function, making it trickier to compute compared to other constraint-based problems such as ones we've seen in previous chapters.

Notably, these two measures of efficiency are not only useful for evolutionary processes, but for bioengineering as well. Obviously, the molar yield has economical implications when, for example, producing ethanol from sugar. However, the rate of a bioprocess is important as well due to the costs involved, e.g. for maintaining an operational bioreactor. One can imagine a computational model that accurately predicts the enzyme demand per flux of a pathway. Choosing the pathways with the lowest demand would be a good strategy for increasing the overall rate of bioproduction [173].

We define the enzyme demand per unit flux as the total amount of enzyme (in grams of protein) that is required to catalyze all of the pathway reactions at their required rates. We start by deriving a formula for the demand of a single enzymatic reaction. Consider an enzyme-catalyzed reaction:

$$S \rightleftharpoons P \qquad (6.4)$$

where s and p will be the concentrations of the substrate (S) and product (P) respectively, and E the concentration of the enzyme which catalyzes this reaction (for simplicity, we drop the tot subscript from Etot). Here, we will be using the factorized rate law (Eq. 3.10), but other kinetic rate laws would produce similar results. The rate of a reaction is given by:

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

where $k_{\text{cat}}^+$ is the forward turnover rate, $K_S$ and $K_p$ are the Michaelis-Menten constants for the S and product P, and $\Delta_{\text{r}}G'$ is the Gibbs free energy. So, the minimal amount of enzyme that is required for reaching a given rate v is:

$$q \equiv v \cdot h \cdot \frac{1}{k_{\text{cat}}^{+}} \cdot \frac{1 + p/K_{\text{P}} + s/K_{\text{S}}}{s/K_{\text{S}}} \cdot \left(1 - e^{\Delta_{\text{r}}G'_{\text{R}}/RT}\right)^{-1}, \qquad (6.6)$$
where:

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

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

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

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

**Problem 6.10 The factorized rate law** Use the Haldane relationship:

$$\frac{k_{\text{cat}}^{+} K_{P}}{k_{\text{cat}}^{-} K_{S}} = \frac{k_{1} k_{3} k_{5}}{k_{2} k_{4} k_{6}} = K_{\text{eq}} \qquad (6.20)$$

and the definition of Gibbs free energy:

$$\begin{aligned} \Delta G'_{r} &= -R \cdot T \cdot \ln K_{\text{eq}} \\ \Delta G'_{r} &= \Delta G'_{r}^{\circ} + R \cdot T \cdot \ln ([P]/[S]) \end{aligned} \qquad (6.21)$$

to prove that Eq. (6.19) is equivalent to the following factorized rate law:

$$v = [E_0] k_{\text{cat}}^{+} \cdot (1 - e^{\Delta_r G' / RT}) \cdot \frac{[S]/K_S}{1 + [S]/K_S + [P]/K_P} \qquad (6.22)$$
<figure><img src="image_64.png" /></figure>

Figure 6.2: Enzyme cost in metabolism – (A) Enzyme-specific flux depends on a number of physical factors. Under ideal conditions, an enzyme molecule catalyses its reaction at a maximal rate given by the enzyme’s forward catalytic constant (blue). The rate is reduced by microscopic reverse fluxes (magenta) and by incomplete saturation with substrate, causing waiting times between reaction events, or by enzyme inhibition or incomplete activation (red). (B-C) On a logarithmic scale, catalytic rates and enzyme demand can be split into sums of efficiency terms. With lower catalytic rates, larger amounts of enzyme are required for realizing the same metabolic flux.

where *h* is a number converting enzyme concentration *e* into enzyme amount *q* (for example, the enzyme molecular mass). For an illustration, see Figure 6.2. Summing up the demand across all the reactions in the pathway (each with its own rate, kinetic parameters, and substrate/product concentrations) will produce the total enzyme demand. Looking at this function, we can already make some interesting observations. First, the kinetic parameters (*k*<sub>cat</sub><sup>+</sup>, *K*<sub>p</sub>, and *K*<sub>s</sub>) can be treated as constants since they change only in evolutionary timescales, and we often assume that existing enzymes already have near-optimal kinetics (although that's not always the case). Since we care about the demand *per pathway flux* one can, without loss of generality, set *v* to 1. However, if the pathway requires a non-trivial ratio between some reactions, the value of *v* can be different based on the stoichiometry. Finally, the thermodynamic term, i.e. $1 - e^{\Delta_r G' / RT}$ (which we will discuss in more detail in the following section, 6.3), is a function of the metabolite concentrations and the *K*<sub>eq</sub>, which is another constant. So, generally speaking, enzyme demand is defined by a set of constants that are unique to each pathway, and variables that represent the metabolite concentrations. Since these concentrations are subject to change depending on the growth conditions, we often treat them as optimization variables and try to find the minimal demand possible within certain constraints. In Section 6.4, we will see a general method for finding the minimal value using convex optimization.

Most of the proposed criteria for good pathways have either to do with material investments (such as
substrate, cofactor, or energy demand) or with "machine investments", that is, enzyme demands. Enzyme
demands, in turn, depend on pathway length, enzyme masses, and enzyme efficiency, and therefore on
rate laws (where kcat values, thermodynamic forces, and metabolite concentrations come into play). In
fact, many criteria which we discussed earlier as indicators of efficiency are actually an approximation of
the enzyme demand under certain assumptions. For example, the number of steps is proportional to the
total demand if all enzymes have exactly the same kcat, saturation, and thermodynamics. Therefore, it
is quite a useful rule-of-thumb in case not much else is known about the enzymes themselves. A better
approximation, denoted Pathway Specific Activity, was used by [165] to compare CO₂ fixation cycles. If we
assume that all enzymes are fully saturated and irreversible, the demand would be a direct function of
the individual enzyme specific activities (specifically, proportional to the sum of all their reciprocal values).
But even if we know nothing about the enzyme kinetic parameters, thermodynamics alone can provide
# Chapter 7

## Metabolism in states of maximal enzyme efficiency

Andreas Kremling, Wolfram Liebermeister, Elad Noor and Meike T. Wortel

### Chapter highlights

* Enzyme-efficient states are metabolic states that realize a given flux objective at a minimal enzyme cost.

* In models without further constraints, enzyme-efficient states are Elementary Flux Modes (EFMs).

* Elementary Flux Modes can be used to find enzyme-efficient states in networks that would be too large to optimize metabolic states "by brute force".

* Biomass per enzyme efficiency can be converted to into cell growth rate.

* As growth conditions are changing, the flux profile either changes continuously (and metabolite and enzyme concentrations as well) or fluxes change discontinuously together with metabolite and enzyme concentrations.

### 7.1. Introduction

In a simple economic picture of cells, we assume that cells adjust their metabolic state in each environment to obtain a maximal fitness advantage. This may be impossible in reality, but it remains an interesting question what this best metabolic state would look like, according to our knowledge of cells. So what is the best metabolic state overall (comprising metabolic fluxes, metabolite concentrations and enzyme levels)? What pathways should a cell use, which enzymes should be induced or repressed, and how should this change in a new environment? To answer this, we need to remember that all metabolic variables (fluxes, metabolite levels, enzyme levels, and enzyme efficiencies) depend on each other. Physically, fluxes depend on metabolite concentrations through kinetics and enzyme regulation (e.g. competitive inhibition) and metabolites are produced and consumed by the fluxes until a steady state is reached. Hence, if we think in terms of cellular economics (treating enzymes as control variables), then all metabolic variables must be optimized together.

In the previous chapters we saw some ways to predict optimal metabolic fluxes, metabolite concentrations and enzyme levels separately: in Flux Balance Analysis (FBA, Chapter 5), we optimized fluxes (while ignoring kinetics or assuming constant enzyme efficiencies); in Enzyme Costs Minimization (Chapter 6), in contrast, fluxes were fixed and given and concentrations (and enzyme efficiencies) were optimized. Each of these methods is based on a strong assumption: FBA requires measured flux ranges and/or apparent catalytic rates and assumes enzyme saturation effects can be neglected, while enzyme cost minimization requires a given flux distribution. But what if we don't know any of the variables in advance? How can we predict all of them from first principles?
Σ Math box 6.C Factorized rate laws and enzyme cost function

According to Eq. (6.13), reversible rate laws can be factorized into five terms that depend on metabolite concentrations in different ways [82]. For a reaction S ⇌ P with reversible Michaelis-Menten kinetics Eq. (6.11), a driving force θ = -ΔrG'/RT, and a prefactor for non-competitive inhibition, the rate law can be written as

$$v = e \cdot k_{\text{cat}}^{+} \cdot \underbrace{[1 - e^{-\theta}]}_{\eta^{\text{for}}} \cdot \underbrace{\frac{s/K_{\text{S}}}{1 + s/K_{\text{S}} + p/K_{\text{P}}} \cdot \frac{1}{1 + x/K_{\text{I}}}}_{\eta^{\text{sat}}}$$

Rate = enzyme level · forward catalytic constant · thermodynamic factor · saturation factor · regulation factor

with inhibitor concentration x. The product of the first two terms, E and $k_{cat}^+$, represents the maximal velocity, i.e. the rate at full substrate-saturation without backward flux and without enzyme inhibition. The following factors decrease this velocity for different reasons: $η^{for}$ describes a decrease due to backward fluxes, $η^{sat}$ – the decrease due to incomplete substrate saturation, and $η^{reg}$ – the decrease due to small-molecule regulation (see Figure 6.2 b). While $k_{cat}^+$ is an enzyme-specific constant (yet, dependent on conditions such as pH, ionic strength, or molecular crowding in cells), the efficiency factors are concentration-dependent, unitless, and can vary between 0 and 1. The thermodynamic factor $η^{for}$ depends on the driving force (and thus, indirectly, on metabolite concentrations), and the equilibrium constant is required for its calculation. The saturation factor $η^{sat}$ depends directly on metabolite levels and contains the $K_M$ values as parameters. Enzyme regulation by small molecules yields additive or multiplicative terms in the rate law denominator, which in our example and can be captured by a separate factor $η^{reg}$. The enzyme cost for a flux v, with an enzyme burden $h_e$, can be written as

$$q = h_e \cdot e = h_e \cdot v \cdot \frac{1}{k_{\text{cat}}^{+}} \cdot \underbrace{\frac{1}{[1 - e^{-\theta}]} \cdot \frac{1 + s/K_{\text{S}} + p/K_{\text{P}}}{s/K_{\text{S}}} \cdot \underbrace{[1 + x/K_{\text{I}}]}_{1/\eta^{\text{reg}}}}_{1/\eta^{\text{for}}}$$

and contains the terms from the rate law in inverse form. The first factors, $h_e v/k_{cat}^+$, define a minimum enzyme cost, which is then increased by the following efficiency factors. By omitting some of these factors, one can construct simplified enzyme cost functions with higher specific rates, or lower enzyme demands (compare Figure 6.2b). For a closer approximation, the factors may be substituted with constant numbers between 0 and 1.

us with useful information with which to grade pathways. Specifically, the $K_{eq}$ of a reaction is a universal constant that is not affected by enzymes, but rather determined solely by the chemical structures of the substrates and products.

In the following sections, we will focus on enzyme use efficiency as a main objective and consider a thermodynamic approximation, relating enzyme demands to thermodynamic forces. For linear metabolic pathways, optimal enzyme profiles (and the associated metabolite profiles and enzyme costs) can be computed with closed formulae. We will also discuss a way to compute optimal enzyme profiles numerically, for networks of any shape and size, as long as the flux mode is known.

## 6.3 The role of thermodynamics

In general, when considering larger metabolic networks, thermodynamic feasibility can play an important or even crucial role in determining which pathways are used. In this section we will discuss this role more explicitly and see how thermodynamics can still give us useful insights about pathway efficiency even when no other kinetic data is available.

Why are thermodynamic driving forces a meaningful criterion for good pathways? In brief, the driving forces, defined as $θ = -ΔrG'/RT$, play a double role: first, they determine whether or not a pathway flux is feasible at all, given the metabolite concentrations at the pathway boundary (i.e. the metabolites that form connections to the broader metabolic network); and second, in case the pathway is feasible, driving forces can affect enzyme efficiency and, consequently, the enzyme demand for a given desired pathway
Before thinking about this, let us briefly step back: what do we actually mean by an "optimal state"? What quantity should be maximized in metabolism? There could be very different aims (e.g. production in biotechnology, vs offspring and survival in a wild-type cell). However, in both cases an important aim is cell growth – or at least, avoiding strong growth deficits. Below we will see that cell growth depends, to a good approximation, on biomass/enzyme efficiency, that is, biomass production per total enzyme invested. Hence, whenever fast growth is important, cells should maximize this efficiency.

Thus, we will consider the following optimality problem: maximize biomass/enzyme efficiency, defined as the production flux per invested enzyme with respect to all metabolic variables (metabolites, enzymes and fluxes) and under all constraints (steady state, enzyme kinetics, etc). Solutions to this problem are considered optimal states.

## 7.2. The enzyme-efficient metabolic states are Elementary Flux Modes

The optimization problem in this chapter is to reach maximal objective flux(es) for minimal enzyme investment. The biological interpretation is that this would lead to the highest growth rate, because it optimizes the ratio between gains (flux(es)) and costs (enzymes). When we solve this optimization problem with mathematical tools, it is convenient to either find the minimal enzyme investment for a certain flux, or the maximum flux for a fixed enzyme investment. Although one could think of different biological explanations for those two ways to state the optimization problem, mathematically they are equivalent. For the outline of the proof that optimal states are Elementary Flux Modes, it is convenient to fix the objective flux to an arbitrary value (we choose 1) and then minimize the enzyme investment. This leads to the following optimization problem:

$$
\begin{array}{ll}
\underset{\mathbf{v},\mathbf{e},\mathbf{c}}{\text{minimize}} & \displaystyle\sum_{i=1}^{r} h_i e_i \\
\text{subject to:} & \mathbf{N} \cdot \mathbf{v} = \mathbf{0} \\
& \forall i: v_i = e_i f_i(\mathbf{c}) \\
& \mathbf{e}, \mathbf{c} \ge \mathbf{0} \\
& v_r = 1 \\
& \mathbf{c} \le c_{\max}
\end{array}
\qquad (7.1)
$$
steady state
enzyme kinetics
positive concentrations
fixed objective flux
metabolite bounds

where $\forall i$ means for all reactions $i$, $r$ is the number of reactions (with the last the objective), $h_i$ are the weights, and bold face items denote vectors. This optimization problem states that by adjusting the fluxes ($v$), metabolite concentrations ($c$) and enzyme concentrations ($e$), the cost (sum of the costs ($h_i e_i$) for every reaction) is minimized, while keeping the objective flux constant. The weights ($h_i$) can be thought of as the size or production costs of the enzymes (measured in e.g. molecular weight or gene length), which might be different for different enzymes. We need certain constraints to be met for this solution: 1) the metabolic network needs to be in steady state to avoid built-up of intermediates, 2) enzyme kinetics: the flux of each reaction ($v_i$) has to be equal to the enzyme concentration ($e_i$) times a metabolite dependent term ($f_i(c)$), 3) all enzyme metabolite concentrations have to be positive, and 4) the objective flux is equal to 1. Optionally, we can add bounds on the metabolite levels, which is mostly necessary with irreversible kinetics. Reversible kinetics usually lead to bounded metabolite levels, because very high concentrations of products will inhibit the reaction forming the products.

In this section we will explain why the optimal state is reached at an Elementary Flux Mode (EFM). One important starting point is that, as we have seen before in Chapter 4, convex optimization problems with only positivity or equality constraints (no other inequalities) lead to an optimal solution at an extreme points of the feasible solution space, and those extreme points are Elementary Flux Modes (EFMs). However, the optimization problem (7.1) is not convex, mainly due to the hyperbolic dependency of reaction rates on the metabolite concentrations ($f_i(c)$ is usually not linear). There
flux. In Chapter 3, we learned that $\Delta_r G'$, and hence the driving force $\theta$, depends on the equilibrium constant $K_{eq}$ of the reaction and on the substrate and product concentrations. We also learned that for a flux in forward direction, the driving force must be positive. Beyond that, the efficiency of an enzyme is proportional to $\eta^{\text{for}}(\theta) = 1 - e^{-\theta}$, a function that ranges between 0 (for $\theta = 0$), reactions in thermodynamic equilibrium) and 1 ($\theta \gg 1$, reactions far from equilibrium). Let us now see how this non-equilibrium relation affects pathway efficiency.

### 6.3.1 Enzyme kinetics and driving forces

We should remind ourselves some of the lessons learned in Chapter 3. Specifically, recall the factorized rate law [82] with a reversibility term that is an explicit function of the Gibbs energy (Eq. 3.10):

$$v = e \cdot k_{\text{cat}}^{+} \cdot \frac{\prod_{i} s_{i}^{n_{i}} / K_{s}}{1 + \prod_{j} p_{j}^{n_{j}} / K_{p} + \prod_{i} s_{i}^{n_{i}} / K_{s}} \cdot (1 - e^{\Delta_{r} G' / RT}). \quad (6.7)$$

The enzyme mechanism behind this formula assumes fast binding and unbinding of substrate and product, and a slow reversible conversion step (of bound substrate into bound product). Note that here we generalize the rate law for cases with more than one substrate and one product, where $n_i$ and $n_j$ are the stoichiometric coefficients of substrates and products, respectively<fnref n="1" />. This generalization is one out of many, and corresponds to the assumption that all reactants bind independently to the enzyme (and at random order). We focus on this rate law because it is one of the simplest, but the theoretical results in this chapter apply to most other generalizations as well (e.g. convenience kinetics [174]).

According to the definition of $k_{\text{cat}}^{+}$, and also by noticing that the middle and rightmost terms in Eq. (6.7) are each smaller than 1, the rate of an enzymatic reaction is bounded by $v \le e \cdot k_{\text{cat}}^{+}$ (see Mathematical Details Box 6.C for a detailed explanation). However, the additional terms are often much lower than 1, which means that the rate does not reach its maximum. If we try to measure the apparent catalytic rate by dividing the rate by the enzyme abundance ($k_{\text{app}} = v/e$) we would typically get a value that is lower than $k_{\text{cat}}^{+}$, while only in rare “ideal” cases, $k_{\text{app}}$ would approach the $k_{\text{cat}}^{+}$. In fact, this reasoning was used by Davidi et al. [24] to estimate the $k_{\text{cat}}^{+}$ values of more than 100 enzymes in E. coli, where they sampled many growth conditions and took the maximum $k_{\text{app}}$ as the estimate.

As discussed in Section 3.2.3, the factorized rate law has a thermodynamic perspective based on the flux-force relationship, where we view the reversibility term as a “penalty” for the fact that by lowering the energy barrier, enzymes must catalyze reactions in both directions. When the driving force ($\theta$) is low, the reverse reaction flux can become significant and lower the net flux. On the other hand, if the driving force is large enough, this term can be ignored and the rate law resembles irreversible kinetics.

So far we've seen that increasing the driving force of a single reaction translates to a better enzyme efficiency and lower demand. If we consider whole pathways, ones whose overall driving force is larger have more of it to distribute among the reactions and therefore should also have higher efficiencies overall. However, using “too much” driving force can also have downsides. Using a larger amount of the Gibbs energy to drive the pathway reactions means that less of that energy would go for building biomass or currency metabolites such as ATP. An example for this trade-off between the efficiency of single enzymes (in terms of backward rates) and the overall pathway efficiency (in terms of ATP yield) was demonstrated by Flamholz et al. [175] who analyzed two versions of the famous glycolytic pathway (see Figure 6.1 below).

### 6.3.2 Driving forces should not be too small

With the factorized rate law 6.7, we can approximate the reaction rates by $v \le e k_{\text{cat}} (1 - e^{-\theta})$ (where we assume positive fluxes by convention). The thermodynamic efficiency $\eta^{\text{for}} = 1 - e^{-\theta}$ plays a prominent role. As shown in Figure 6.3, this formula yields two important approximations: for small forces $\theta$, that

<footnote n="1">¹In general, reaction stoichiometries can be arbitrarily scaled. For example, instead of a reaction $2 A \rightarrow B$, we may write $A \rightarrow \frac{1}{2} B$ for convenience, which will only lead to a scaling factor in the reaction rate. However, this holds only if reaction stoichiometries are used to describe mass-balance. In cases like Eq. (6.7), where stoichiometries appear in kinetic rate laws or in thermodynamic balances, we do not have this choice. In these cases, the stoichiometries must reflect the molecularities, that is, the actual number of reactant molecules involved in the enzymatic reaction.</footnote>
<figure><img src="image_50.png" /><figcaption>Figure 7.1: Translation from flux to enzyme space retains EFMs as extreme rays. The top left panel shows the feasible flux space with the steady state constraints, all fluxes positive (using splitting of fluxes, as explained in the text, if necessary) and a fixed objective flux. The extreme points here are points where one flux becomes 0 and are Elementary Flux Modes (see Chapter 5). Here we show that when we have assumed metabolite concentrations, such as when we keep them at an optimal solution, we get a linear transformation and the extreme rays are maintained. Different metabolite levels, for example solutions to different environmental conditions, can lead to different transformations and therefore different optima (minimal total enzyme), but those are always located at an EFM.</figcaption></figure>

are several ways to proof that the solution of this optimization problem is an EFM, of which some are outlined in the papers by Wortel et al. [242] and Müller et al. [243]. In this book we will outline the proof that a solution that would not be an EFM, will lead to a contradiction.

Assume we have some optimal state that is not an EFM. Any optimal solution is associated with a set of fluxes, and enzyme and metabolite concentrations. Now we set the metabolite concentrations to the concentrations of the assumed optimal state. Then all metabolite dependent terms ($f_i(c)$) become constants and we return to a convex problem. As explained in Chapter 9 and Figure 7.1, the optimum (now in terms of enzyme concentrations and fluxes) of this problem is a flux distribution that is an EFM. But that is a contradiction: we starting with assuming the optimal state from which we took the set of enzyme concentrations was not an EFM. We conclude that the optimal state has to be an EFM.

## 7.3. Illustration with an example network

To illustrate the proof, we study a simple network that we have seen previously in Chapter 5 (Figure 7.2), in which we will now include enzyme kinetics. We will use the factorized rate law as in Chapter 6, but then generalized for $n_s$ substrates and $n_p$ products (also compare Eq. (3.15) in Chapter 3):

$$v = e \cdot k_{\text{cat}}^{+} \cdot \frac{\prod_{n_s}^{j=1} s_j / K_{S,j}}{1 + \prod_{n_p}^{k=1} p_k / K_{P,k} + \prod_{n_s}^{j=1} s_j / K_{S,j}} \cdot (1 - e^{\Delta G'_{r}/RT}) \quad (7.2)$$
<figure><img src="image_65.png" /><figcaption>Figure 6.3: The thermodynamic efficiency term $\eta^{\text{for}}$ and some approximations – (A) In a given reaction, the thermodynamic efficiency term $\eta^{\text{for}} = 1 - e^{-\theta}$ (solid line) can vary between 0 and 1 depending on the driving force $\theta$. Small driving forces make the enzyme inefficient, since $\eta^{\text{for}} \to 0$, while for large forces, thermodynamics does not play a role as $\eta^{\text{for}} \to 1$. The dashed lines show two linear approximation that hold always as bounds, but can also be used as good approximations for small or large $\theta$ values, respectively: $(1 - e^{-\theta}) < \theta$ and $(1 - e^{-\theta}) < 1$. (B) The reciprocal value $1/\eta^{\text{for}}$ is one of the factors determining enzyme demand. The solid line shows the thermodynamic demand factor $1/\eta^{\text{for}}$, while the dashed lines show the resulting approximations $1/\eta^{\text{for}} > 1/\theta$ and $1/\eta^{\text{for}} > 1$, corresponding respectively to the enzyme demand approximations $e \ge \frac{v}{k_{\text{cat}}\theta}$ and $e \ge \frac{v}{k_{\text{cat}}}$.</figcaption></figure>

is, close to equilibrium, we obtain $\eta^{\text{for}} \approx \theta$, while for large forces, that is, for strongly forward-driving reactions, we obtain $\eta^{\text{for}} \approx 1$. In fact, both approximations also serve as upper bounds across all $\theta$ values. What does this mean? Far from equilibrium, the thermodynamic term does not play a role and can be ignored. Close to equilibrium, in contrast we obtain a simple approximation for fluxes

$$v < e \cdot k_{\text{cat}}^{+} \cdot (1 - e^{-\theta}) < e \cdot k_{\text{cat}}^{+} \cdot \theta \qquad (6.8)$$

and hence for the enzyme demand

$$e > \frac{v}{k_{\text{cat}}^{+} \cdot (1 - e^{-\theta})} > \frac{v}{k_{\text{cat}}^{+} \cdot \theta} \qquad (6.9)$$

As $\theta$ goes to zero, the enzyme demand (for a given desired flux) goes to infinity. We already know the reason from Chapter 3: the driving force determines the ratio of forward and reverse one-way fluxes, $\frac{v^+}{v^-} = e^\theta$. If $\theta$ comes close to zero, their relative difference becomes very small, and in order to obtain a given net flux $v = v^+ - v^-$, both $v^+$ and $v^-$ must grow enormously, which would require an a large amount of enzyme. This effect concerns only very small $\theta$ values - for $\theta$ much larger than 1 (or $\Delta_r G'_{\text{G}}$ much smaller than -RT), it can be neglected. Therefore, redistributing driving forces between reactions, to avoid very small forces, can save enzyme costs. The relation between driving forces, enzyme efficiency enzyme demand is shown in more detail in Figure 6.4.

If small driving forces should be avoided to prevent enzyme costs from going infinity, how can this happen in practice? The driving forces themselves depend on metabolite levels, which can vary over several orders of magnitude. While the true metabolite concentrations are usually unknown, we hypothesize that selection favors concentration profiles that prohibit very small driving forces, in order to escape the ensuing large enzyme demands. Of course, completely avoiding small driving forces may be impossible, as there is always a trade-off: if a metabolite concentration decreases, the driving forces of all reactions producing it will increase, but the driving forces of all reactions consuming it will decrease simultaneously. So, all else being equal, the optimal metabolite profile is one that distributes its driving forces as evenly as possible.
<table><caption>Table 7.1: Outcomes of the optimization of the example network with standard kinetics, parameter values and external concentrations (see Box 7.4) for varying levels of [G<sub>ex</sub>].</caption><thead><tr><th>[G<sub>ex</sub>]</th><th>E<sub>tot</sub></th><th>v<sub>0</sub></th><th>v<sub>1</sub></th><th>v<sub>2</sub></th><th>v<sub>3</sub></th><th>v<sub>4</sub></th><th>v<sub>BM</sub></th><th>e<sub>0</sub></th><th>e<sub>1</sub></th><th>e<sub>2</sub></th><th>e<sub>3</sub></th><th>e<sub>4</sub></th><th>e<sub>BM</sub></th><th>[G]</th><th>[P]</th><th>[ATP]</th><th>[ADP]</th></tr></thead><tbody><tr><td>0.01</td><td>156.2</td><td>5</td><td>5</td><td>0</td><td>9</td><td>0</td><td>1</td><td>54.4</td><td>4.4</td><td>0</td><td>94.4</td><td>0</td><td>2.9</td><td>0.08</td><td>15.14</td><td>0.05</td><td>20.09</td></tr><tr><td>0.1</td><td>91.3</td><td>50</td><td>50</td><td>99</td><td>0</td><td>0</td><td>1</td><td>61.3</td><td>11.3</td><td>14.2</td><td>0</td><td>0</td><td>4.4</td><td>0.13</td><td>4.55</td><td>0.11</td><td>20.09</td></tr><tr><td>1</td><td>36.2</td><td>50</td><td>50</td><td>99</td><td>0</td><td>0</td><td>1</td><td>13.0</td><td>8.0</td><td>12.5</td><td>0</td><td>0</td><td>2.7</td><td>0.60</td><td>7.65</td><td>0.11</td><td>20.09</td></tr></tbody></table>

See Box 7.4 for all detailed rate laws of the example networks. We can simplify this equation by combining the forward catalytic constant, the thermodynamic efficiency factor, the saturation efficiency factor, and the regulation efficiency factor (if that exists) in a function *f(c)*, which only depends on the metabolites, and not on the enzyme concentrations. We will below write *f* for *f(c)*.

$$v_i = e_i \cdot f_i \qquad (7.3)$$

Now we take $v_{BM} = 1$ and optimize all fluxes, enzymes concentrations and metabolite concentrations to minimize the enzyme costs, while satisfying the constraints posed in Equations (7.1), for different levels of external glucose and standard levels of the other external metabolites. We see that for different external glucose concentrations, we get different optimal fluxes, enzyme levels and metabolite levels (Table 7.1).

We notice that the total enzyme needed for the biomass flux of one decreases with the increasing glucose levels, as we expect. Also, the optimal level of internal glucose increases with increasing external glucose. Moreover, the fluxes of the solutions do follow an EFM (see Figure 7.2b). We can now reformulate the problem for only the flux and enzyme levels, while keeping the metabolite levels as they are in the table. With the metabolite levels in the first row of the table, we can linearly relate the enzyme and flux levels (with the factors $f_i$), and thereby the extreme rays of the enzyme and flux space will be equal and EFMs, as pointed out above (see also Chapter 5, Figure 7.1). Optimization in this space will lead to the optimal flux distributions following an EFM. If the fluxes in the table would not follow an EFM, this will lead to an inconsistency; if there is an optimum, this should still be the outcome if we keep some variables in the optimum constant and optimize over the rest of the variables.

We point out two important aspects, using the network (Figure 7.2) as an example. First, it is convenient to split reversible reactions, such that fluxes are always positive. In this case that means that the reversible reaction from P to P<sub>1</sub> is split into the forward reaction v<sub>2</sub> and the reverse reaction v<sub>4</sub>, which both can have only positive flux. This splitting makes sure that EFMs are the extreme rays of the flux space (see Chapter 5). This splitting is purely a mathematical conveniency, we still assume this to be one reaction in the biological sense, and therefore the kinetic equations of both the forward and the backward reaction will be exactly the same. Only, depending on in which direction the flux goes, either one of the reactions will be positive. Any solution with both reactions positive is infeasible, but minimizing enzyme levels will never lead to such a solution, therefore we do not need to set additional constrains. Second, the feasibility of EFMs can depend on external concentrations. In this network the biomass reaction (v<sub>BM</sub>) is the objective flux and there are three EFMs leading to the production of biomass: EFM1 consisting of v<sub>0</sub>, v<sub>1</sub>, v<sub>2</sub> and v<sub>BM</sub>, EFM2 consisting of v<sub>0</sub>, v<sub>1</sub>, v<sub>3</sub> and v<sub>BM</sub> and EFM3 consisting of v<sub>4</sub>, v<sub>3</sub> and v<sub>BM</sub>. However, if P<sub>1</sub> is absent in the environment, the uptake flux v<sub>4</sub> will always be 0 and therefore EFM3 will not be feasible.

## 7.4. Computation of the optimal state

We can now use the result that states of maximal enzyme efficiency are reached at an Elementary Flux Mode to calculate these states. We can follow the following steps to calculate these states:

1. Enumerate the Elementary Flux Modes

2. Calculate the minimal enzyme for each EFM scaled to an objective flux of 1
<figure><img src="image_66.png" /><figcaption>Figure 6.4: Thermodynamic forces, enzyme efficiency, and enzyme demand in a linear chain of reactions – The plot in the center represent two possible profiles of the thermodynamic driving forces (blue and red). The curves describe the cumulative $\Delta_r G'$ values: while the total $\Delta_r G'$ is fixed (and determined by external metabolite concentrations), the shape of the profile can vary. In the optimal profile (in red), small driving forces are avoided. The driving forces determine the ratios of forward and backward one-way fluxes (red arrows), and at a given net flux (black arrows) the enzyme demands. In the suboptimal blue curve, in contrast, the last three reactions show lower forces, and therefore relatively high reverse fluxes (blue arrows); to obtain the same net flux, forward and backward fluxes have to be strongly increased, which increases the enzyme demand.</figcaption></figure>

### 6.3.3 Max-Min driving force method

Previously in Chapter 4.4.3, we discussed adding thermodynamic constraints to constraint-based models in order to comply with the second law of thermodynamics. We can extend that approach in order to implement the idea of avoiding small driving forces. When we talk about the thermodynamic profile of a metabolic pathway, we usually try to visualize it by the cumulative Gibbs energy of reaction: we start at 0 and at each step add the $\Delta_r G'$ of the next reaction, which, assuming the pathway is feasible, is a negative number. The profile therefore has a shape of a downhill slope. The end point represents the total Gibbs energy and depends only on the concentrations of the metabolites that are part of the net reaction. Intermediate metabolites do not affect it, but they do determine the shape of the profile itself (see Figure 6.4). Specifically, each intermediate metabolite typically affects the driving force of two reactions – the one producing it and the one consuming it – with opposite signs. Therefore, changing the concentration of an intermediate can help increase the driving force of one reaction, but always at the expense of another reaction. This strong coupling between $\Delta_r G'$ is why it is not trivial to find the optimal thermodynamic profile of a pathway.

The Max-Min driving force method (MDF) [176] is a method for predicting metabolite concentrations, based on the principle of evenly distributed driving forces. All fluxes are fixed and given, and assumed to be positive. It assumes that each metabolite concentration must remain in a predefined range, converts each choice of metabolite concentrations into the corresponding pattern of driving forces, and determines the smallest resulting driving force in the network. If this smallest driving force is negative, the flux distribution cannot be realized thermodynamically. Otherwise, the larger this smallest driving force, the better the overall metabolite profile. Hence, among all possible metabolite profiles, MDF predicts the one that maximizes the value of the minimal driving force across the network. Mathematically, this leads to a linear optimization problem: in the space of logarithmic metabolite concentrations, a lower bound on all driving forces (denoted $B$) is maximized (Eq. 6.10). An illustrative example is shown in Figure 6.5.
<figure><img src="image_51.png" /><figcaption>Figure 7.2: Maximal efficiency states in example model. (a) Example network from Chapter 5 with added stoichiometry. (b) Three elementary flux modes of this network. (c) Calculated enzyme investment needed for a biomass flux of 1. There is a single concentration of $G_{ex}$ for which the optimal EFM switches from EFM2 to EFM1. (d) Specific fluxes (flux divided by total enzyme) associated with the optimal EFM for different levels of $G_{ex}$. The rates show a discontinuity when the optimal EFM switches from EFM1 to EFM2.</figcaption></figure>

### 3. Compare the EFMs and select the one with minimal enzyme demands

Step 2 is a convex optimization problem as we have seen in Chapter 6. Step 1 is possible for relatively large networks, although usually not for genome scale metabolic networks. The method of these three steps is called Enzyme Flux Cost Minimization, because it is similar to Enzyme Cost Minimization, but while that is focussed on fixes fluxes, Enzyme Flux Cost Minimization simultaneously finds the optimal fluxes. In this section we will show the method on the example network of Figure 7.2.

First, we describe the network with the Stoichiometric matrix (N) and the concentration vector (c):
<figure><img src="image_67.png" /><figcaption>Figure 6.5: Max-Min Driving force method (MDF): an optimality problem in metabolite space – (A) Example pathway with given equilibrium constants and fixed concentrations of the external metabolites X and Y. What are the most favorable concentrations of the internal metabolites A and B? Assuming that small driving forces should be avoided in all reactions, MDF determines the metabolite profile that optimizes a worst case: it maximizes the worst (that is, smallest) driving force among all three reactions. (B) Driving force in reaction 1, as a function of the logarithmic concentrations of A and B, called ln a and ln b. Higher concentrations of A (the reaction product) lead to smaller driving forces. Above a critical value (where X and A are in equilibrium), the driving force becomes negative, and a forward flux is impossible (gray region). The concentration of B, which does not participate in the reaction, does not play a role. (C) Driving force for reaction 2. Here, it is the ratio b/a that counts. The lower the ratio (lower right), the higher the driving force. If the ratio is higher than the equilibrium constant, the driving force becomes negative (grey region). (C) Driving force for reaction 3. (E) By overlaying the contours in (B), (C), and (D) and taking the minimum value, we obtain the minimal driving force θmin among all three reactions. θmin is a piecewise linear function of ln a and ln b within the feasible range, yielding positive forces in all three reactions. The maximum point of this function is the optimum metabolite profile predicted by MDF. In the example shown, the feasible concentration space is entirely defined by the driving forces themselves, given the external concentrations. In general, physiological concentration ranges for all metabolites could further decrease the solution space and shift the optimum point (not shown).</figcaption></figure>

$$
\begin{array}{ll}
\text{Maximize}_{x,B} & B \\
\text{Subject to} & -(\Delta_r \mathbf{G}^{\prime\circ} + RT \cdot \mathbf{N}^{\top} \mathbf{x}) \geq B \\
& \ln(\mathbf{C}_{\min}) \leq \mathbf{x} \leq \ln(\mathbf{C}_{\max})
\end{array}
\qquad (6.10)
$$

MDF is easy to apply: it is based on a simple Linear Programming problem and requires only the following input data: (i) the stoichiometric network; (ii) the flux directions; (iii)) the known equilibrium constants (or equivalently, the standard reaction Gibbs free energies); (iv) physiological ranges for metabolite concentrations. Based on these data alone, metabolite concentrations and driving forces (or ΔᵣG' values) are predicted. An example application can be found in Hädicke et al. [177], where the potential of CO₂ fixation in E. coli via endogenous pathways was analyzed using MDF.

A theoretical insight from MDF is the notion of distributed bottlenecks. A simple bottleneck would consist of a single reaction whose driving force cannot be increased because the substrates are at their upper concentration bounds and the products are at their lower concentration bounds. Given the fixed equilib-
$$
\begin{equation}
\begin{aligned}
\mathbf{N} &= \begin{pmatrix}
1 & -1 & 0 & 0 & 0 & 0 \\
0 & 2 & -1 & -1 & 1 & -1 \\
0 & 2 & 0 & 10 & 0 & -100 \\
0 & -2 & 0 & -10 & 0 & 100
\end{pmatrix} \\
\mathbf{c} &= \begin{pmatrix}
[G] \\
[P] \\
[ATP] \\
[ADP]
\end{pmatrix}
\end{aligned}
\tag{7.4}
\end{equation}
$$

And with the stoichiometric matrix we can describe the steady state constraints:

$$
\dot{\mathbf{c}} = \mathbf{N}\mathbf{v} = \begin{pmatrix} 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 2 & -1 & -1 & 1 & -1 \\ 0 & 2 & 0 & 10 & 0 & -100 \\ 0 & -2 & 0 & -10 & 0 & 100 \end{pmatrix} \begin{pmatrix} v_0 \\ v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_{BM} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \tag{7.5}
$$

Now we find the EFM5 (for example with EFMtool [244]). It is easily cheched that the following EFM5 are in the
nullspace of the Stoichiometric matrix:

$$
\text{EFM1} = \begin{pmatrix} 5 \\ 5 \\ 0 \\ 9 \\ 0 \\ 1 \end{pmatrix}, \quad \text{EFM2} = \begin{pmatrix} 50 \\ 50 \\ 99 \\ 0 \\ 0 \\ 1 \end{pmatrix}, \quad \text{EFM3} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 10 \\ 11 \\ 1 \end{pmatrix} \tag{7.6}
$$

The next step is to do the convex optimization over the metabolite levels for all three EFM5. Therefore we convert the enzyme levels to ratio of the flux and the function f(c), using Equation 7.3. By taking the sum of those we make a function for the total enzyme as a function of fluxes, metabolite concentrations and parameters. We use the standard parameters (Box 7.4) and for each EFM input the fluxes. We are then left with a convex optimization over the metabolite levels, an Enzyme Cost Minimization problem as in Chapter 6. For [Gex] = 0.05 we obtain a total enzyme of 111.1 for EFM1 and of 146.3 for EFM2. That means that for this conditions we will conclude that EFM1 is optimal, and we obtain the metabolite concentrations from the optimization of [G] = 0.08, [P] = 3.93, [ATP] = 0.11 and [ADP] = 20.09. We can next use the rate equations to calculate the enzyme levels from the fluxes and metabolite levels, using the values for the parameters and external concentrations.

We can repeat this procedure for different levels of external concentrations and see that the optimal EFM can change
depending on the external concentration (Figure 7.2c). When the optimum shifts to using a different EFM, there is a
discontinuity in the fluxes at the external metabolite concentration (Figure 7.2d). Many cells show shifts in metabolic
strategies depending on the external conditions and with Enzyme Flux Cost Minimization those shift could be explained.

In conclusion, we have now found the metabolic state of maximum enzyme efficiency. Even though we obtained the enzyme concentrations last, this is how the cell controls the system. When the cells makes enzymes in the concentrations we calculated, a steady state will set in with the fluxes and metabolite levels that will lead to this optimal state.
rium constant, nothing can be done to increase the driving force in this reaction. A distributed bottleneck is more complicated: it consists of a series of reactions that all share the same low driving force, which, because of all the concentration constraints in the system, cannot be further increased (e.g. as in Figure 6.4). Even though each single reaction looks “harmless” because its own driving force could still be increased, this increase would happen at the expense of other driving forces.

### 6.3.4 The role of thermodynamics for metabolic states

In summary, thermodynamics provides important clues both about the feasibility of pathways fluxes and about their enzyme demand. To use this knowledge, fluxes need to be considered together with metabolite concentrations (to obtain the possible driving forces), but no detailed knowledge of enzyme kinetics is required. Thermodynamics alone yields an upper bound on fluxes (and hence, a lower bound on enzyme demands) that holds for any kinetic rate laws. The only required data (except for the metabolic network itself) are equilibrium constants (or equivalently, standard Gibbs free energies of reactions $Δ_rG°$), which can be obtained from the eQuilibrator tool (equilibrator.weizmann.ac.il) [178, 68, 69] as well as physiological bounds on metabolite concentrations. Given this information, and given a feasible choice of metabolite concentrations, we can compute the driving forces of all reactions, and from the factorized rate law (and assuming positive fluxes by convention) we can then approximate the reaction rates by $v ≤ e k_{cat} (1−e⁻⁰)$.

We also recall from Chapter 3 that driving forces are not independent between reactions, but depend on the metabolite concentrations, which creates trade-offs: in a chain $A \xrightarrow{R_1} B \xrightarrow{R_2} C$, a lower concentration of B will increase the driving force in $R_1$, but decrease the driving force in $R_2$. For high enzyme efficiency (low enzyme demand), all driving forces should in principle be high, but this is most important for low $θ$ values (while for $θ ≫ 1$ it does not even matter). Therefore we may conclude that, to save enzyme, a cell should rearrange its metabolite levels within physiological bounds such that small $θ$ are avoided. Implementing this as an optimality problem, we obtain MDF.

In conclusion, we described (i) a general rule of thumb that poor thermodynamics makes reactions costly; (ii) simple approximations of enzyme cost; and (iii) practical methods (MDF) to obtain metabolite profiles with favorable thermodynamic properties.

## 6.4 Enzyme cost minimization

### 6.4.1 Enzyme cost minimization

The problem of minimizing the total enzyme demand (or cost) for a given pathway can be solved numerically, thanks to the fact that they are always convex [179]. Finding the minimum of the convex objective (the total enzyme cost) in a convex set (the set of admissible metabolite profiles, a convex polytope in log-metabolite space) can be done efficiently. In contrast to general optimality problems, such problems have a unique local optimum, which can be found by simple numerical methods. In this section, we demonstrate it with a simple example, the same three-reaction pathway that you already saw in Section 6.3 above.

### 6.4.2 Enzyme cost landscape of a metabolic pathway

Given the fluxes, kinetics, and concentration bounds in a metabolic pathway model, one can predict the enzyme demand by assuming that cells minimize the enzyme cost in that pathway. In the Enzyme Cost Minimization method A reaction rate $v = e \cdot κ(s)$ depends on enzyme level $e$ and metabolite concentrations $s_i$ through the enzymatic rate law, $κ(s)$. If the metabolite concentrations were known, we could directly compute enzyme demands $e = v/κ(s)$ from fluxes, and similarly calculate the flux-specific enzyme demand $e/v = 1/κ(s)$. However, metabolite concentrations are usually unknown and vary between experimental conditions. Therefore, there can be many solutions for $e$ and $s$ realizing one flux distribution. To select one of them, we employ an optimality principle: we define an enzyme cost function (for instance, total enzyme mass) and choose the enzyme profile with the lowest cost while restricting the metabolite levels to physiological ranges and imposing some thermodynamic constraints. As we shall see below, the solution is in many cases unique.
<figure><img src="image_52.png" /><figcaption>Figure 7.3: Translation of enzyme-specific biomass rate to growth rate. A Both from experimental data and a cell-optimization point of view, the ribosomal fraction of the proteome increases with the growth rate, while the metabolic fraction decreases. B This leads to a hyperbolic dependency of the growth rate on the biomass production rate.</figcaption></figure>

Box 7.A: Kinetics of the example network

The detailed kinetic equations for the example model (Figure 7.2) using the factorized rate law (see Equation (7.2) and Chapters 3 and 6) are:

$$
\begin{align}
v_0 &= e_0 \cdot k_{\text{cat},0}^+ \cdot \frac{[\text{G}_{\text{ex}}]/K_{\text{G}_{\text{ex}}}}{1 + [\text{G}]/K_{\text{G}} + [\text{G}_{\text{ex}}]/K_{\text{G}_{\text{ex}}}} \cdot \left(1 - e^{\Delta G'_{r,0}/RT}\right) \nonumber \\
v_1 &= e_1 \cdot k_{\text{cat},1}^+ \cdot \frac{([\text{G}]/K_{\text{G}})/([\text{ADP}]/K_{\text{ADP}})}{1 + ([\text{P}]/K_{\text{P}})/([\text{P}]/K_{\text{P}})/([\text{ATP}]/K_{\text{ATP}}) + ([\text{G}]/K_{\text{G}})/([\text{ADP}]/K_{\text{ADP}})} \cdot \left(1 - e^{\Delta G'_{r,1}/RT}\right) \nonumber \\
v_2 &= e_2 \cdot k_{\text{cat},2}^+ \cdot \frac{([\text{P}]/K_{\text{P}})}{1 + [\text{P}_1]/K_{\text{P}_1} + [\text{P}]/K_{\text{P}}} \cdot \left(1 - e^{\Delta G'_{r,2}/RT}\right) \nonumber \\
v_3 &= e_3 \cdot k_{\text{cat},3}^+ \cdot \frac{([\text{P}]/K_{\text{P}})/([\text{ADP}]/K_{\text{ADP}})/([0_2]/K_{0_2})}{1 + ([\text{C0}_2]/K_{\text{C0}_2})/([\text{ATP}]/K_{\text{ATP}}) + ([\text{P}]/K_{\text{P}})/([\text{ADP}]/K_{\text{ADP}})/([0_2]/K_{0_2})} \cdot \left(1 - e^{\Delta G'_{r,3}/RT}\right) \tag{7.7} \\
v_4 &= e_4 \cdot k_{\text{cat},4}^+ \cdot \frac{([\text{P}_1]/K_{\text{P}_1})}{1 + [\text{P}_1]/K_{\text{P}_1} + [\text{P}]/K_{\text{P}}} \cdot \left(1 - e^{\Delta G'_{r,4}/RT}\right) \nonumber \\
v_{\text{BM}} &= e_{\text{BM}} \cdot k_{\text{cat,BM}}^+ \cdot \frac{([\text{P}]/K_{\text{P}})/([\text{ATP}]/K_{\text{ATP}})}{1 + ([\text{BM}]/K_{\text{BM}})/([\text{ADP}]/K_{\text{ADP}}) + ([\text{P}]/K_{\text{P}})/([\text{ATP}]/K_{\text{ATP}})} \cdot \left(1 - e^{\Delta G'_{r,0}/RT}\right) \nonumber
\end{align}
$$

Note that P is a product twice in $v_1$, as $v_1$ produces 2P. Note that $v_2$ and $v_4$ have a very similar rate equation, but in the different direction, and both should only be used in the positive direction.

The standard set of parameters we used for the toy model is all $k_{cat,i} = 10$ except $k_{cat,3} = 0.1$, all $\Delta G'_{r,i} = -1000$, $R = 8.3$, $T = 293$ and all $K_M = 1$. For the external metabolites $[P_1] = 1$, $[G_{ex}] = 0.05$, $[0_2] = 0.1$, $[BM] = 1$ and $[C0_2] = 10$ unless mentioned otherwise.

7.5. Translating enzyme efficiency into cell growth rate

In the section above, we learned how to optimize metabolic states for a maximal overall enzyme efficiency. Why is this quantity relevant? One reason is that overall enzyme efficiency, according to some simple reasoning, determines the cell growth rates. If microbes compete by growing fast, their fitness is largely determined by their momentary growth rate in their respective environment. In such environments, the biomass/enzyme efficiency will be under selection, which makes it one of the important objective functions in this book. If higher enzyme efficiency means higher growth rate, and if we have a conversion formula for this, we can put "growth rate" instead of "overall enzyme efficiency" on the axes of our plots.

Enzyme-efficient metabolic states allow us to compute specific biomass production rates, i.e. the rate of biomass production per metabolic enzyme invested. If biomass consisted only of enzymes, the ratio "enzyme production rate / enzyme concentration" would give us directly the growth rate. However, biomass does not only consist of metabolic enzymes,