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}}}} \quad (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}}}. \quad (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$ [81] 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)).

### 6.4.3 Enzyme cost as a function of metabolite concentrations

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_{\text{cat}}$ values, but also on the kinetic rate law, which includes $K_{\text{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_{\text{cat}}^{+}$ (i.e. the maximal possible forward rate per unit of enzyme, in $s^{-1}$), 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
<figure><img src="image_53.png" /><figcaption>Figure 7.4: Optimal growth rates of the two EFMs for different levels of the external metabolite G<sub>ex</sub> (using Equation 7.8)</figcaption></figure>

but includes ribosomal enzymes, RNA, DNA, lipids, and other compounds. Therefore we need a formula for converting
biomass/enzyme efficiency into cellular growth rate.

Mathematically, a cell's growth rate is given by $\mu = v_{BM}/c_{BM}$, where $v_{BM}$ is the biomass production rate (biomass produced per cell volume and time) and $c_{BM}$ is the biomass amount per cell volume. If a cell consisted only of metabolic enzymes (more precisely, of the enzymes described in our model) the biomass/enzyme efficiency $K_{BM} = v_{BM}/h_{enz}$ would directly describe the cellular growth rate. Since that is not the case, we need to convert between $h_{enz}$ and $c_{BM}$. The metabolic protein fraction decreases with the growth rate, leading to a hyperbolic dependency of the growth rate on the biomass production rate (Figure 7.3). We may use the empirical approximation $h_{enz}/c_{BM} = f_{prot}(a - b \mu)$, where $f_{prot} = 0.5$ is the fraction of protein mass within the cell dry mass and the parameters $a = 0.27$ and $b = 0.2$ h were fitted to describe the metabolic enzyme fraction in proteomics data, assuming a linear dependence on growth rate [18]. This yields the conversion formula (see also [241]):

$$ \mu = \frac{a f_{prot} v_{BM}}{h_{enz} + b f_{prot} v_{BM}}. \qquad (7.8) $$

This formula has been used to convert the minimal enzyme cost per biomass flux for different external concentrations in the toy model (Figure 7.2c) to the maximal growth for each EFM (Figure 7.4).

7.6. Conclusions

In this chapter we considered the metabolic network of a cell - and enzyme levels, metabolite concentrations, and fluxes as the state variables - and studied its maximally efficient states. Finding such states can be difficult because fluxes, metabolite concentrations, and enzyme levels are tightly coupled: metabolite concentrations determine enzyme efficiencies, enzyme efficiencies determine optimal enzyme levels, and enzyme levels determine fluxes and metabolite concentrations, which in turn determine enzyme efficiencies. To find an optimal state, all variables need to be optimized at the same time, which is a non-linear optimality problem with (possibly) many local optima. In small toy models, solutions can be found numerically, but for large detailed models, the computational effort becomes enormous. Instead of simplifying the problem (as in the previous chapters) we here used the insight that (in models without extra flux bounds) the optimal solutions must be EFMs. Biologically, it is enzyme levels rather than fluxes that are regulated directly (by transcriptional regulation), while metabolite concentrations and fluxes respond dynamically. But from a functional point
<figure><img src="image_68.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>

rate laws like the Michaelis-Menten kinetics in Eq. (6.11) can be factorized as [82]:

$$v = e \cdot k_{\text{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{s}, v_l) = v_l \cdot \frac{1}{k_{\text{cat},l}^+} \cdot \frac{1}{\eta_l^{\text{for}}(\theta(\mathbf{s}))} \cdot \frac{1}{\eta_l^{\text{sat}}(\mathbf{s})} \cdot \frac{1}{\eta_l^{\text{reg}}(\mathbf{s})}. \qquad (6.14)$$

To determine the enzyme demand of an entire pathway, we sum over all reactions: $e_{\text{tot}}^{\text{pw}} = \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
of view, we may see this differently: we may first think of a task (a flux to be realized) and then ask how this flux can be performed optimally.

Recommended readings

M.T. Wortel, H. Peters, J. Hulshof, B. Teusink, and F.J. Bruggeman. Metabolic states with maximal specific rate carry flux through an elementary flux mode. FEBS Journal, 281(6):1547-1555, 2014.

S. Müller, G. Regensburger, and R. Steuer. Enzyme allocation problems in kinetic metabolic networks: Optimal solutions are elementary flux modes. Journal of Theoretical Biology, 347:182-190, 2014.
ribosomes, as well as effects of misfolding and non-specific catalysis. The enzyme burdens $h_{el}$ can be chosen heuristically, for instance, depending on enzyme sizes, amino acid composition, and lifetimes. Setting $h_{el} = 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 [40, 180]. Thus, each reaction flux $v_l$ is associated with an enzyme cost $q_l$, which can be written as a function $q_l(v_l, s) \equiv h_{el} e_l(s, v_l)$ of flux and metabolite concentrations. From now on, we refer to log-scale metabolite concentrations $s_i = \ln s_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_{el} e_l(v_l, \mathbf{s}) = \sum_{l} h_{el} \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 $\mathbf{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 $\mathcal{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.

1. **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 [181]. 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.

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

3. **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 many of them can be reformulated as approximations of this method [176, 182, 170].
4. 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 [183], 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+\frac{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 l_i}$. 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. [179] 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,
# Chapter 8

## Models of growing cells

Ohad Golan, Hollie J. Hindley, Hidde de Jong, Markus Köbis, Elena Pascual Garcia, and Andrea Weiße

### Chapter highlights

* A comprehensive description of fundamental growth laws in microbial growth, elucidating the core principles that govern biological growth patterns.

* A detailed exploration of the contrasts between coarse-grained and fine-grained modeling is presented, offering insights into the varying levels of detail that each approach encompasses.

* A thorough breakdown of the key assumptions in the modeling of metabolic systems is provided, underlining the foundational premises that are crucial for accurately representing these complex systems.

* The process of deriving fundamental growth laws by modeling key assumptions is meticulously demonstrated, enabling a clear understanding of how theoretical constructs translate into biological realities.

### 8.1. Introduction

A key feature of living systems is that they are able to grow and reproduce. The reproductive success in a given environment defines the fitness of a living system. The study of the growth of bacteria and other microorganisms is crucial for better understanding their capacity to cause diseases in humans or for better exploiting their use in biotechnological or environmental processes. Beyond their interest for a variety of applications, bacteria and other microorganisms have shown themselves ideal model systems for investigating fundamental questions on the relation between growth, fitness and characteristics of the environment.

One of the first to systematically and quantitatively study the growth of bacterial cultures was Jacques Monod in the 1940s. He performed so-called diauxic growth experiments, in which bacteria were cultured in a medium containing two different limiting carbon sources. He showed that the bacteria first deplete one carbon source before starting to assimilate the second carbon source. The order in which the primary and secondary carbon source were consumed was determined by the growth rate they support: the preferred carbon source allows the culture to grow at a higher rate. Further work on the molecular basis of diauxic growth led to the discovery that cells inhibit the expression and activity of functions for the use of secondary carbon sources when a preferred carbon source is present, a global regulatory mechanism known as carbon catabolite repression [245, 246].

Monod characterized bacterial growth by means of batch culture experiments in a well-defined growth medium allowing bacteria to reach a state of balanced growth, where the accumulation of biomass can be described by a single constant, the exponential growth rate. Together with the chemostat, a device allowing continuous culture of microorganisms at a predefined growth rate [247], these methods have become standard in microbial physiology. They notably underlie
<figure><img src="image_69.png" /><figcaption>Figure 6.7: Comparing two metabolic pathways using the Max-min Driving Force (MDF) method. The light blue line represents the cumulative Gibbs free energy along the pathway if all metabolite concentrations were 1 mM. The MDF solution is presented as a gray line, where the bottleneck reactions are highlighted in red.</figcaption></figure>

since it involves electron transfer 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. [175] 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 [184]) and a promising platform for bio-production [185] – 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 [186]. 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% [187] (see Economic Analogy Box 6.D)

To see if indeed the models provide predictions that are consistent with the experimental evidence, Flamholz et al. [175] 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).

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. [175] 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
<figure><img src="image_54.png" /><figcaption>Figure 8.1: Growth laws of bacterial growth. A. Monod growth law: growth rate dependency on nutrient availability [Data from [250]]. B. Correlation between growth rate and nutrient uptake rate [Data from [256]]. C. Correlation between growth rate and cellular composition [Data from [254]]</figcaption></figure>

the discovery of a number of so-called growth laws, relating the growth rate to a variety of properties of the physiology of growing bacteria. The growth laws are conserved across different organisms and a broad range of experimental conditions. Here, we list three well-known growth laws [248, 249]:

1. Dependency of the growth rate on nutrient availability [250]: In his characterization of bacterial growth, Monod discovered the first growth law. He observed that the growth rate of bacteria depends upon the nutrient concentration in the medium in a hyperbolic fashion.

2. Correlation between growth rate and nutrient uptake rate [251]: In continuous cultures, the growth rate was shown to vary linearly with the nutrient uptake rate. The slope of this linear relation is called the biomass yield and the offset the 'maintenance energy', as it is assumed to be derived from the energy spent on processes required to maintain the basic processes of the cell, in the absence of growth [252].

3. Correlation between growth rate and cellular composition [253, 254]: In 1959, Schaechter, Maaløe and Kjeldgaard showed that RNA, DNA and the number of nuclei in *Salmonella typhimurium* linearly correlate with the growth rate. Later, it was further shown that other physiological parameters, such as the mass fraction of ribosomes in growing populations, alsolinearly correlate with the growth rate [254]. Initially, it was believed that the correlation between ribosomal mass fraction and growth was strictly positive, however, Scott et al. [255] showed that when growth is inhibited through translation-inhibiting drugs, growth rate and ribosomal mass fraction exhibit a negative (near-)linear relation.

The conserved nature of the growth laws has led scientists to ask whether there are fundamental principles governing bacterial growth. To answer this question, different types of mathematical models have been developed. One approach aims at integrating all known molecular constituents of the cell and the reactions involving these constituents into a big model, an *in-silico* copy, or 'digital twin', of the cell. Such models, known as fine-grained models, can be useful to predict emergent phenotypes, but they are difficult to construct and maintain, and their complexity makes it hard to grasp certain principles that underpin growth. In this chapter, we will focus on coarse-grained models of bacterial growth. Rather than assembling individual reactions in a bottom-up manner, these models are based on the top-down definition of a limited number of basic cellular functions or processes involved in growth, described by appropriate macro-reactions. Coarse-grained models are smaller and therefore easier to construct and analyze. The lack of molecular detail can make their predictions less accurate, but their simplicity allows a focus on how basic cellular functions and their interactions shape bacterial growth. How much detail is included in a model depends on the specific scientific question asked, and similarly, models may vary in their underlying assumptions. Oftentimes, assumptions are based on biochemical principles governing intracellular reactions, on physical limitations faced by cells, on optimality principles, or on a combination of these.

In this chapter, we show how to understand and, ultimately, how to develop coarse-grained models of cellular growth. We present a number of coarse-grained models with increasing levels of granularity. The models have been chosen
<figure><img src="image_70.png" /></figure>

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

<figure><img src="image_72.png" /><figcaption>Figure 6.8: Comparing two metabolic pathways using the Enzyme Cost Minimization (ECM) method. We used the same kinetic parameters for all enzymes in both pathways ($k_{cat}^{+} = 200 s^{-1}$, $K_{M} = 200 \mu M$, same as in [175]). However, here we used an updated version of Enzyme Cost Minimization with the factorized rate law, therefore the results are not identical. A Jupyter notebook for generating the figure can be found on the book website.</figcaption></figure>

for catalyzing the same flux (see Figure 6.8). 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. [164] 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. [188] used the MDF algorithm and the Pathway Specific Activity measure to compare between designs of photo-electro-autotrophic strains. Similarly, Löwe and Kremling [189] used the Enzyme Cost Minimization algorithm to predict the enzyme demand of both natural and synthetic carbon fixation cycles.

## 6.5.3 Predicting the metabolite 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 [175], Jacobson et al. [190] measured the intracellular concentrations ED intermediates in *Z. mobilis*, and used them to calculate the Gibbs
<figure><img src="image_55.png" /><figcaption>Figure 8.2: Coarse grained modelling of cellular growth. Compared to genome-scale FBA and whole-cell models, coarse grained models zoom out of the molecular detail and focus on key processes.</figcaption></figure>

to also represent a variety of commonly used assumptions, for example, based on growth rate maximization or on phenomenological or mechanistic constraints. Despite these differences, however, models we discuss generally recover the basic growth laws, and we show how the latter can be derived from solving two of the simplest coarse-grained models. The goals of this chapter are:

1. To enable the reader to understand and analyze any model of microbial growth from the literature.

2. To enable the reader to develop their own coarse-grained model of a metabolic system that is directed at their specific scientific question.

3. To provide the reader with a new perspective on modeling of complex systems and specifically the biological cell.

## 8.2. Fundamental modeling assumptions of microbial growth

The models of microbial growth we consider here are based on fundamental assumptions that follow from biochemical and biophysical constraints. In this section, we discuss and mathematically define assumptions that are found, explicitly or implicitly, in most coarse-grained models of microbial growth. The assumptions are formulated in an abstract manner to hold for any self-replicating biological system, irrespective of the specifics of the underlying molecular mechanisms. In the next section, we use these assumptions to construct increasingly complex models of microbial growth and show how the latter can be used to derive the experimentally observed growth laws presented in the introduction of this chapter.

The growth of microorganisms consists of the uptake of nutrients from the environment and the conversion of these nutrients into new microbial cells through a number of coupled metabolic processes. This description brings out the self-replicating or autocatalytic nature of microbial growth: cells transform nutrients from the environment into new cells. In what follows, we consider growth on the population level, that is, an increase in the total amount of cells or, equivalently in many situations, an increase of the biomass of the population. This leads to the well-known model of microbial growth, where the change in biomass over time is proportional to the amount of biomass ():

$$ \frac{dB}{dt} = \lambda B, \qquad (8.1) $$

where $t$ [h] denotes time, $B$ in gram dry weight [gDW] the biomass and $\lambda$ [1/h] the population growth rate.

If the growth rate is constant, the solution to Eq. (8.1) describes exponential growth of the biomass:

$$ B = B_0 e^{\lambda t}, \qquad (8.2) $$

where $B_0$ [gDW] is the initial biomass at $t = 0$.
<figure><img src="image_73.png" /><figcaption>Figure 6.9: The *E. coli* model of central metabolism used in Noor et al. [179].</figcaption></figure>

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, Figures 6.9 and 6.10) [179]. In a related paper, Wortel et al. [191] 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_74.png" /></figure>

### Economics 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.
<figure><img src="image_56.png" /><figcaption>Figure 8.3: Coarse grained models of metabolic systems with increasing complexity. A. A self replicating system. B. The simplest description of a metabolic system: coupled catabolic and anabolic reactions. C. A metabolic system that can catabolize two different nutrient sources. D. A catabolic system requiring two substrates to grow: x and E.</figcaption></figure>

The growth rate is a key parameter that is often used as a proxy for the fitness of microorganisms. It is dependent on the metabolic processes, that is, how a cell utilizes the nutrients to synthesize new biomass (self-replication). The simplest description of metabolism is that it takes up a nutrient, breaks it down into metabolites (catabolism), and then utilizes these metabolites to produce new biomass (anabolism). Catabolic and anabolic processes comprise a variety of biochemical reactions that are carried out by different sets of proteins and enzymes. The reaction rates of these processes are limited biochemically and biophysically. We formulate these limitations as modeling assumptions and define them as mathematical constraints, four of which we briefly review below.

### 8.2.1. Conservation of mass and quasi-steady-state assumption

Dry biomass is often a more readily measurable quantity than cell volume. The latter relates absolute abundances of cell components to their intracellular concentrations. Yet, because bacterial cells have been observed to maintain approximately constant cell density across various growth conditions [257, 258] (though transient exceptions have been observed at the single-cell level [259]), biomass can be regarded a proxy for volume and is therefore assumed to be proportional to cell volume in many growth models. All models considered in this chapter are based on the assumption of constant cell density and approximate the concentration x of a cellular component x (we use normal font for cell components and italic font for their concentrations) by its absolute abundance divided by the cell mass.

According to the law of mass conservation, the change of mass is equal to the inflow minus the outflow of mass. As a consequence, the change in concentration of a cell component, for example a metabolite pool, is determined by the sum of the rates of the reactions consuming and producing this cell component. The mass balance for any cell component x is given by the following equation:

$$ \frac{dx}{dt} = \sum_{y} r_{y \to x} - \sum_{k} r_{x \to k}, \qquad (8.3) $$

where $r_{y \to x}$ denotes the rate of the reaction converting cell component y into cell component x (production of x), and $r_{x \to k}$ the rate of the reaction converting cell component x into cell component k (consumption of x). Typically, cell component concentrations have units mg/gDW or mmol/gDW, so that rates of metabolic reactions are expressed in units mg/(gDW h) or mmol/(gDW h).

In the simple system of, there are two reactions: one converting the nutrient source n into metabolite x and one
(A)

**Metabolite concentrations**

predicted [M]

measured [M]

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

(B)

**Enzyme concentrations**

predicted [M]

measured [M]

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

(C)

Total mass density = 59.3 [g/L]

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

Figure 6.10: Validation of metabolite and enzyme concentrations, predicted by Enzyme Cost Minimization, in the central carbon metabolism of *E. coli* (Fig 6.9) – (A) Comparing predicted and measured metabolic concentrations. The thin diagonal line marks *x* = *y*, i.e. where the predictions match the measurements. Full blue points are for all metabolites whose allowed concentration range was set to 1µM - 10 mM. Hollow blue points represent co-factors whose concentration was fixed in the analysis and therefore are not actually predicted – these were omitted from the statistics. The Root Mean Squared Error (RMSE, in log₁₀ scale) was 0.62, the r² (Pearson correlation) was 0.50, and *p*-value was 2.2 × 10⁻⁴. (B) Comparing predicted and measured enzyme concentrations. Here, the RMSE was 0.47, r² = 0.54, and *p* = 5.3 × 10⁻⁶. (C) A pie chart showing the distribution of the predicted absolute mass-concentrations for both enzymes (orange) and metabolites (blue) together. Note that *aconitase* (catalyzing the reactions acn1 and acn2) has a lower specific activity than *glyceraldehyde-3P dehydrogenase* (catalyzing gap), and therefore occupies a higher fraction of the mass-concentration even though the required concentration of the latter enzyme is higher. Labels of enzymes and metabolites that occupy the smallest fractions of the biomass are omitted due to lack of space. Data and model predictions taken from [179].

**6.6 Concluding remarks**

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”
utilizing the metabolite for the synthesis of biomass. According to (8.3), the flux balance of metabolite pool x is given by $dx/dt = r_{n→x} - r_{x→B}$.

A key assumption is that intracellular concentrations are in quasi-steady state. This means that cell component pools remain constant:

$$ \frac{dx}{dt} = 0, \text{ for all cell components } x. \qquad (8.4) $$

The quasi-steady-state assumption simplifies the mathematical analysis of the system significantly and holds for balanced growth of the microbial population. In this chapter, we focus mostly on situations in which the quasi-steady-state assumption applies, but also give an example of a model with metabolic dynamics. In metabolic modeling, the rates of reactions at steady state are called fluxes, denoted by the symbol $J$. With the quasi-steady-state assumption, Eq. (8.3) becomes

$$ \sum_{y} J_{y→x} = \sum_{k} J_{x→k} \qquad (8.5) $$

that is, for every cell component, the sum of production fluxes equals the sum of consumption fluxes. In the example system, we have $J_{n→x} = J_{x→B}$.

### 8.2.2. Proteome allocation assumption

The biochemical reactions breaking down nutrients into intracellular metabolites, and the reactions utilizing these metabolites for the synthesis of new biomass, do not occur spontaneously. The reactions are catalyzed mostly by proteins complexes, in particular metabolic enzymes and ribosomes. In coarse-grained models, well-defined sets of biochemical reactions are grouped together into macro-reactions. The cell components that are necessary to catalyze the individual steps of a macro-reaction are grouped together into a corresponding so-called proteome sector. A proteome sector includes mostly proteins that catalyze metabolic reactions but also ribosomes catalyzing the reaction of protein biosynthesis. Proteins constitute most of the biomass of the cell [260]. As a first approximation, therefore, the sum of the proteome sectors equals the total biomass of the growing population measured in units of $g$:

$$ \sum_{r \in \{x \to y\}} P_r = B, \qquad (8.6) $$

where $P_{x→y}$ is the proteome sector catalyzing the macro-reaction that transforms cell component x into cell component y. The proteome sectors as defined above are extensive quantities, summed over the entire growing population, like the total biomass $B$. For the models, we are rather interested in intensive quantities, the amount of a proteome sector relative to the total amount of biomass (protein), corresponding to protein concentrations or protein fractions. Dividing the left-hand and right-hand sides of Eq. (8.6) by $B$, we thus obtain:

$$ \sum_{r \in \{x \to y\}} p_r = 1 \qquad (8.7) $$

where $p_{x→y}$ is the fraction of the proteome converting x into y, defined by $p_{x→y} = P_{x→y}/B$. Proteome fractions are dimensionless and sum to one.

In the simple example of, we distinguish two macro-reactions: a catabolic reaction and an anabolic reaction (biomass synthesis). We therefore define two proteome sectors, corresponding to enzymes and ribosomes, respectively, with fractions $p_{n→x}$ and $p_{x→B}$, respectively. In later examples in the chapter, the catabolic and anabolic processes are further broken down into smaller macro-reactions and so are the proteome sectors.

### 8.2.3. Mathematical description of reaction fluxes

The rate at which a reaction is converting one cell component, e.g., a metabolite, into another is determined by the proteome fraction, the concentrations of the substrates of the reaction and possible regulation by other cell components
<figure><img src="image_78.png" /></figure>

Economics analogy 6.E Two central 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 equilibrium– 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.

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 Cost Minimization. This method is directly related to 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.

In order to predict optimal metabolic states, we started in the previous chapter with models that optimize the fluxes in an entire network. However, to keep the models linear, kinetics and concentrations were largely ignored. In FBA with molecular crowding, a connection between fluxes and enzyme levels was made via empirical parameters, the apparent catalytic rates or "enzyme efficiencies". We now saw that these parameters are not at all constant parameters, but emerge from kinetics and given concentration profiles, and we also saw how optimal concentration profiles can be computed for a given flux distribution. This means: we now know how to predict optimal fluxes from known enzyme efficiencies, and we know how to predict optimal concentrations (and therefore enzyme efficiencies) from known fluxes. In the next chapter we will put these two things together, in order to predict all variables in the system – fluxes, metabolite concentrations, enzyme efficiencies, and enzyme levels – from a single principle of maximal overall enzyme efficiency.
in the system. While mass-action kinetics provide a principled framework to develop rate equations for biochemical reactions, in practice, various approximations based on mechanistic assumptions are often used to obtain simplified equations [261]. Below are a few examples of rate laws defining the fluxes in coarse-grained models:

1. Excess substrate and no allosteric interactions. The simplest relation of the flux to the relevant proteome sector is linear, such that

$$J_{x \to y} = p_{x \to y} \beta_{x \to y}, \qquad (8.8)$$

where $\beta_{x \to y}$ is a parameter describing the efficiency of proteome sector $p_{x \to y}$ in generating a flux from $x$ to $y$. This expression assumes substrate $x$ is in excess and disregards any regulation of the flux by allosteric interactions of the enzymes and other cell components.

2. Limited substrate and allosteric interactions. A more complex relation is obtained when the substrate is in excess and allosteric interactions involving a cell component n play a role in the modulation of the flux. The expression of the flux is multiplied by two regulatory functions $f(x)$ and $g(n)$ describing the modulation of the flux by the substrate and the allosteric cell component, respectively:

$$J_{x \to y} = p_{x \to y} \beta_{x \to y} f(x) g(n). \qquad (8.9)$$

It is important to note that that both $f(x)$ and $g(n)$ return values between 0 and 1, and that the flux remains linear in the proteome fraction. Typically, a Michaelis-Menten relation is taken for the effect of the concentration of substrate $x$ on the flux, such that $f(x) = x/(k_{x \to y} + x)$. When the concentration $x$ is in excess, such that $x \gg k_{x \to y}$, the function $f(x)$ becomes approximately 1. Other types of regulatory functions can be used depending on the macroreactions concerned and the growth conditions.

### 8.2.4. Volume and surface area assumptions

The intracellular volume as well as the surface area of the cell are limited. Obviously, the total volume occupied by the components of the cell, in particular proteins, cannot be larger than the cell volume. As such, the total volume of the cell is larger than the sum of the volume of the proteome sectors that are functioning inside the cell plus some constant volume taken up by other cell components such as DNA. This gives the following constraint:

$$\text{Cell volume} \geq \sum_{r \in \{x \to y\}} p_r v_r + v_0 \qquad (8.10)$$

where $v_{x \to y}$ is the volume of proteome sector $p_{x \to y}$ and $v_0$ is some constant volume filled by other cell components. Similarly, the total surface occupied by proteins and lipids making up the cell membrane has to equal the surface area of the cell. This constraint gives:

$$\text{Cell surface area} \geq \sum_{r \in \{x \to y\}} p_r s_r + l_0 \qquad (8.11)$$

where $s_{x \to y}$ is the surface area of proteome sector $p_{x \to y}$ and $l_0$ is the surface area of the lipids in the cell membrane.

### 8.3. Derivation of growth laws from basic modeling assumptions

In the following section, we will build upon the fundamental assumptions discussed earlier to construct models of microbial metabolism with increasing complexity. We will introduce additional assumptions as necessary to solve each model, and use them to derive one of the growth laws presented in the introduction of this chapter that have been experimentally observed in microorganisms.
## 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, 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, 2016. doi: 10.1371/journal.pcbi.1005167.

## Problems

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

### 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_r 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 a glycolysis model with a linear flux-force relationship

### 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 Cost Minimization solutions (or vice versa)? Otherwise, can you think of an approximation of the Enzyme Cost Minimization problem, such that MDF provides the correct solution? Show how the Enzyme Cost 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_{\text{tot}}^{\text{pw}}$ 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
<figure><img src="image_57.png" /></figure>

Conservation of mass and steady state

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

$$ \frac{d \text{ cell component}}{dt} = \text{incoaming flux} - \text{outoing flux} = 0 $$

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

Proteome allocation

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

$$ p_{\text{catabolism}} + p_{\text{biosynthesis}} + \cdots + p_{\text{uptake}} = 1 $$

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

Volume and surface area

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

$$ p_{\text{catabolism}}v_{\text{catabolism}} + \cdots + v_0 = \text{cell volume} $$

$$ p_{\text{Glucose uptake}}S_{\text{Glucose uptake}} + \cdots + s_0 = \text{cell surface area} $$

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

Reaction fluxes

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

$$ J_{x \to y} = p_{x \to y} \frac{s}{s + k_m} $$

**Figure 8.4:** Fundamental assumptions in the modeling of microbial growth. A. Conservation of mass and steady-state assumption: The change in concentration of a cell component is equal to the incoming flux minus the outgoing flux. At steady state, the concentration of the cell component is constant. B. Proteome allocation assumption: the proteome is divided into different proteome sectors. The number of proteome sectors in a model depends on the model granularity. The sum of all the proteome sectors always equals 1. C. Volume and surface area assumption: The volume of the cell is limited and is filled with intracellular cell components such as proteins. The sum of the volumes of the intracellular cell components is equal to the cell volume. Similarly, the surface area of the cell is limited and contains membrane cell components such as lipids. The sum of the surface areas of membrane cell components is equal to the cell surface area. D. Example of flux assumption according to Michaelis-Menten kinetics: the reaction $x \to y$ is carried out by proteome sector $p_{x \to y}$. The maximal rate is reached for saturating substrate concentrations and is determined by the size of the proteome sector.

**Example 1 - Basic metabolic system with saturating substrate concentrations** In this example, we will use the basic metabolic model to derive the relationship between the concentration of ribosomes and the growth rate in microorganisms. The most basic metabolic model involves the uptake of a single nutrient from the environment, the catabolism of that nutrient into a metabolite x, and the use of this metabolite in anabolic processes to synthesize biomass. This model consists of two reactions and two proteome sectors. According to the proteome allocation constraint, the sum of the proteome sectors must sum to one (according to Eq. (8.7)):

$$ p_{n \to x} + p_{x \to B} = 1. \qquad (8.12) $$

For simplicity, we assume that the rate of each reaction is proportional to the allocation of the proteome to that reaction (according to Eq. (8.8)), so that:

$$ J_{n \to x} = p_{n \to x} \beta_{n \to x}; \quad J_{x \to B} = p_{x \to B} \beta_{x \to B}. \qquad (8.13) $$

The mass conservation constraint, with the assumption of a steady state for metabolite x, gives (according to Eq. (8.5)):

$$ J_{n \to x} = J_{x \to B}. \qquad (8.14) $$
Prove that the function:

$$f(\mathbf{e}) = \frac{1}{\sum_{i} (A_i e_i)^{-1}} \qquad (6.16)$$

for a fixed A and under the constraint $\sum_i e_i = e_{\text{tot}}$, is at its maximum when:

$$e_i = e_{\text{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_{\text{S}} - k_{\text{cat}}^{-} [P] / K_{\text{P}}}{1 + [S] / K_{\text{S}} + [P] / K_{\text{P}}} \qquad (6.19)$$

where:

$$K_{\text{S}} = \frac{k_2 k_4 + k_2 k_5 + k_3 k_5}{k_1 (k_3 + k_4 + k_5)}; \quad K_{\text{P}} = \frac{k_2 k_4 + k_2 k_5 + k_3 k_5}{k_6 (k_2 + k_3 + k_4)}; \quad k_{\text{cat}}^{+} = \frac{k_3 k_5}{k_3 + k_4 + k_5}; \quad 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_{\text{cat}}^{-}} \frac{K_{\text{P}}}{K_{\text{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_{\text{r}} G^{\prime\circ} &= -R \cdot T \cdot \ln K_{\text{eq}} \\ \Delta_{\text{r}} G' &= \Delta_{\text{r}} G^{\prime\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_{\text{r}} G' / RT}) \cdot \frac{[S] / K_{\text{S}}}{1 + [S] / K_{\text{S}} + [P] / K_{\text{P}}} \qquad (6.22)$$
Finally, due to conservation of mass, the biomass synthesis flux equals the growth rate:

$$ \lambda = J_{x \to B}. \qquad (8.15) $$

Solving equations (8.12)-(8.15) gives a prediction for the growth rate:

$$ \lambda = \frac{\beta_{x \to B} \beta_{n \to x}}{\beta_{x \to B} + \beta_{n \to x}}. \qquad (8.16) $$

Solution (8.15) for the growth rate is based solely on *mechanistic assumption* - that is, assumptions that are based on the mechanistic properties of the biochemical reactions in the cell. In this case, that is that the fluxes are linear to the relevant proteome sector. Because we have taken a steady state approximation and the rates of the two reactions must be equal, the growth rate is determined by the relative values of the catalytic constants.

Using this model, we can now derive the relationship between the concentration of ribosomes and the growth rate. Combining Eq. ((8.12)) and ((8.14)) gives:

$$ \lambda = p_{x \to B} \beta_{x \to B} \qquad (8.17) $$

This shows that the growth rate is linearly proportional to the anabolic sector. Given that the anabolic sector is composed mostly of ribosomes, this fits well with the experimentally observed linear relationship between the concentration of ribosomes and the growth rate, which was first described by Schaechter et al.[253] and later confirmed by Bremer et al.[254]. It is important to notice that this relation is due to the assumption that the biomass synthesis flux is linear in the ribosomal proteome sector.

In summary, we have derived the linear relationship between the concentration of ribosomes and the growth rate using only basic assumptions about the properties of the biochemical reactions in the cell and the conservation of mass. This relationship is one of the experimentally observed growth laws in microbial systems.

**Example 2 - Growth on two nutrient sources** In this example, we consider a metabolic system that grows on two different nutrient sources, $n_1$ and $n_2$. We use the fundamental assumptions outlined in Section 1.2 and an additional assumption of growth-rate maximization to demonstrate how cells may exhibit catabolite repression - a phenomenon in which cells utilize only one nutrient even when multiple nutrients are available in the environment [246].

The metabolic system in this example catabolizes both nutrient sources to the same metabolite x, but at different efficiencies. The anabolic reaction is the same as in Example 1. There are now three proteome sectors in this model: two for catabolism of the nutrients and one for anabolism. Thus, according to the proteome allocation constraint (Eq. 8.7), we have:

$$ p_{n_1 \to x} + p_{n_2 \to x} + p_{x \to B} = 1. \qquad (8.18) $$

As before, we assume a linear correlation between reaction rates and proteome sector fractions (according to Eq. (8.8)). The different efficiencies of the catabolic sectors is represented as $\beta_{n_2} > \beta_{n_1}$. Applying the mass conservation assumption for metabolite x, combined with the steady state assumption, gives

$$ J_{n_1 \to x} + J_{n_2 \to x} = J_{x \to B}. \qquad (8.19) $$

The growth rate is again equal to biomass synthesis flux, as in Example 1:

$$ \lambda = J_{x \to B} \qquad (8.20) $$

Given that there are more variables than constraints in this example, solving eqs reveals that there is no unique solution
for the growth rate, but rather a solution space with one free variable $p_{n_1 \to x}$:

$$ \lambda = \frac{\beta_{x \to B} \beta_{n_2 \to x}}{\beta_{x \to B} + \beta_{n_2 \to x}} + p_{n_1 \to x} \left( \frac{\beta_{n_1 \to x} - \beta_{n_2 \to x}}{\beta_{x \to B} + \beta_{n_2 \to x}} \right) \beta_{x \to B} \quad (8.21) $$

The solution shows that the metabolic system has a decision to make regarding how much of the proteome to invest in sector $p_{n_1 \to x}$. To solve this system, we introduce an additional assumption of *growth rate optimization* – that is, to maximize its fitness, the metabolic system maximizes the growth rate in a given condition. In this example, to maximize the growth rate, the cell uses only the more efficient catabolic system, setting $p_{n_1 \to x} = 0$ and the solution for the growth rate is as in example 1. The model predicts that the cells will only utilize the nutrient source with the higher efficiency, even if both nutrient sources are available in the environment. This solution fits the catabolic repression experimental result presented in the introduction in which in which the metabolic system represses the use of a less efficient nutrient source in favor of a more efficient one.

## Problems

**Problem 8.1** A system is composed of a set of 2 linear reactions: nutrient -> metabolite $x_1$ -> metabolite $x_2$ -> biomass. Using the same approximations as in example 1, solve for the growth rate. What would be the solution for a system composed of $N$ reactions? Show that the least efficient reaction determines the growth rate.

**Problem 8.2** Solve Example 2 when the nutrients are not available in excess. Use Michaelis-Menten relations for the catabolic reactions. At what point does the metabolic system switch to use the other nutrient source?

**Problem 8.3** A simple model with allosteric regulation of catabolic reaction: A metabolic system is growing in an environment with one nutrient available. The system regulates its catabolic reaction according to the concentration of metabolite $x$. What is the growth rate as function of catabolic sector proteome size? What is the catabolic sector proteome size that maximizes the growth rate? (
# Chapter 7

## Optimization of metabolic states

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

### Chapter overview

* Optimal metabolic states in this chapter refer to enzyme-efficient states, which are metabolic states that realize a given objective flux at a minimal enzyme cost.

* In models without further flux constraints, flux distributions of 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 cell growth rate by simple approximate formulae.

* The Elementary Flux Mode that is realized in an enzyme-efficient state depends on the external conditions.

* As growth conditions change, either the flux profile changes continuously (together with metabolite and enzyme concentrations), or fluxes change discontinuously, implying jumps also in 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 these questions, 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 by maximizing an objective function (typically biomass) while in Enzyme Cost Minimization [175, 179] (Chapter 6) metabolite concentrations were optimized by minimizing cost (or, equivalently, maximizing the enzyme efficiencies). 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?

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, versus number of 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.

In conclusion, 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 Enzyme-efficient metabolic states use elementary flux modes

The optimization problem in this chapter is to reach maximal objective flux with minimal enzyme investment. The biological interpretation is that this would lead to the highest growth rate, because it optimizes the ratio between gains (fluxes) 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 over the fluxes (v), enzymes levels (e) and internal metabolite concentrations (s):

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

where $h_i$ are the weights, $\mathbf{N}$ is the stoichiometry matrix, $i$ is the index of the reactions (ranging from 1 to $r$), with the last reaction (with index $r$) representing the objective. This optimization problem states that by adjusting the fluxes (v), metabolite concentrations (s) and enzyme concentrations (e), the total cost (sum of costs – $h_i e_i$ – for every reaction) is minimized, while keeping the objective flux constant (any arbitrary constant can be chosen, here we chose 1). The weights ($h_i$) can be thought of as the size or production costs of the enzymes (measured, for example, in molecular weight or gene length) We require certain constraints: (i) the metabolic network needs to be in steady state to avoid built-up of intermediates, (ii) enzyme kinetics – the flux of each reaction ($v_i$) has to be equal to the enzyme concentration ($e_i$) times a metabolite dependent (e.g., saturation) term ($\kappa_i(s)$), (iii) all enzyme and metabolite concentrations have to be positive, (iv) the objective flux is equal to 1, and (v) the metabolite concentrations are within their given bounds. The latter constraint is optional and is mostly necessary when dealing with irreversible kinetics. Reversible kinetics usually lead to bounded metabolite levels because very high concentrations of products inhibit the reaction that forms 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 point of the feasible solution space, and those extreme points are Elementary Flux Modes. However, the
# Chapter 9

## Large resource allocation models of cells

Hugo Dourado, Anne Goelzer, Pranas Grigaitis, Wolfram Liebermeister, and Elad Noor

### Chapter highlights

* Resource balance analysis models are cell models based on three basic constraints formulated at genome-scale: stationary fluxes (balancing production and consumption fluxes, uptake and excretion fluxes, as well as compound dilution by cell growth); flux coupling constraints relating fluxes to the amounts of catalyzing enzymes (or other machines); and density constraints, limiting molecule amounts in cell compartments, or molecule concentrations.

* These constraints narrow down the solution space predicted by FBA towards more physiological solutions

* Large resource allocation models build on the same principles, and have been implemented as different variations (RBA models in a narrow sense, ME-models, and pc-models).

### 9.1. Introduction: Quantifying cell investments at genome-scale

#### 9.1.1. Can we build detailed models of entire growing cells?

In the previous chapter on growing cells, we derived possible growth rates and corresponding cell compositions and metabolic strategies based on resource allocation ideas. We asked what amounts of cellular machines will allow a cell to self-replicate, and possibly, to self-replicate at a maximal possible speed. However, the models described were rather simple, with reactions coarse-grained into proteome sectors. Following the same idea, would it be possible to model an entire cell by a one-to-one mapping of all of its reactions, including realistic kinetics and find the configuration for this entire system that would allow the cell to optimally self-replicate? Such models would in principle be capable of describing the detailed resource allocation into each reaction in the cell: they would extend genome-scale models of metabolism (such as in FBA) to models covering all cell subsystems, at the same time accounting for all important constraints on growth, in particular realistic kinetics.

So why don't we directly describe metabolism and macromolecule production in one single, detailed model? If we do that and attempt to find the optimal state, we quickly see that the optimality problems for genome-size models are hard to solve, and maybe unsolvable with current means. The reason is simple: due to non-linearities in the rate laws (where rates are typically dependent on substrates and machines, and are therefore at least quadratic), the optimality problem may have multiple local optima. Since the space for optimization is very high-dimensional (with the number of dimensions on the order of the number of compounds and reactions considered) and non-convex (because certain regions of the initially convex solution space are excluded by thermodynamic constraints), the global optimum may be hard (or in practice,
optimization problem (7.1) is not convex, mainly due to the hyperbolic dependence of reaction rates on the concentrations of metabolites ($\kappa_i(s)$ is usually not linear in the internal metabolite concentrations).

There are several ways to prove that the solution of this optimization problem is an EFM, of which some are outlined in the papers by Wortel et al. [192] and Müller et al. [193]. Here we will outline a proof by contradiction: assuming a solution to the optimization problem that is not an EFM and showing that this leads to a contradiction.

**Theorem 2.** *The flux distribution that maximizes an objective flux over the total enzyme cost in a metabolic network without additional constraints is an Elementary Flux Mode.*

*Proof.* Assume we have some optimal state where the flux distribution is not an EFM. Any optimal solution is associated with a set of fluxes, enzyme concentrations and metabolite concentrations. Now we set the metabolite concentrations to the concentrations of the assumed optimal state. Then all metabolite-dependent terms ($\kappa_i(s)$) become constants, and we return to a convex problem. As explained in Chapter 10 and Figure 7.1, the optimum of this problem (now in terms of enzyme concentrations and fluxes) is a flux distribution that is an EFM. But this contradicts our initial assumption that the optimal state from which we took the set of metabolite concentrations was not an EFM. The proof by contradiction shows that the optimal state must be an EFM. □

## 7.3 Enzyme-efficient states in an example network

To illustrate the proof, we study a simple network representing growth on glucose and pyruvate that we have seen previously in Chapter 5 (Figure 7.2). We use G and P for glucose and pyruvate in the equations, we use the subscript ex when a metabolite is extracellular and square brackets to denote a concentration. For the use in this chapter, we add enzyme kinetics to this network. 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.10) 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_r G' / RT}) \quad (7.2)$$

See Box 7.A for all detailed rate laws of the example network. 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 $\kappa(s)$, which only depends on the metabolites, and not on the enzyme concentrations. We will below write $\kappa$ for $\kappa(s)$.

$$v_i = e_i \cdot \kappa_i \quad (7.4)$$

Now we take $v_{BM} = 1$ and optimize all fluxes, enzyme concentrations and metabolite concentrations to minimize the enzyme costs ($e_{tot} = \sum_i e_i$), 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 concentrations of external glucose, lead to different optimal fluxes, enzyme levels and metabolite levels (Table 7.1).

The table shows that the total enzyme needed for a biomass flux of one decreases with increasing glucose levels, as we expect. In addition, the optimal level of internal glucose increases with increasing external glucose. This is because a higher external glucose allows for a higher internal glucose while still maintaining a steady glucose influx, and a higher internal glucose allows fewer enzymes to drive further metabolism. Moreover, the fluxes of the solutions 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
impossible) to find.

We can also see this more formally. While kinetic self-replicator models could be formulated, the resulting problems would be non-convex [242, 262, 263] and potentially hard to solve. How do we know this? In the optimal solution, the flux distribution must be an elementary flux mode [262]. Since no better method is known, to find the globally optimal metabolic fluxes, protein and metabolite abundances all elementary flux modes must be enumerated. For large-scale metabolic networks this is hard, and in fact this problem is computationally \#P-complete.[264], which prevents its resolution for large-scale models. However, recent attempts of integration of thermodynamics constraints in EFM computation drastically reduced their number [265], which may provide significant help in future developments to solve large self-replicator models.

So in order to model a cell in high resolution (including individual metabolic reactions and their enzymes), what can we do? In previous chapters, we approached this problem from two directions. Flux Balance Analysis (FBA) models, as described in Chapter 5, capture metabolic fluxes in detail, but ignore the rest of the cell (such as macromolecule production; and also the fact that all cellular compounds are diluted by cell growth). Self-replicator models, as described in Chapter 8, describe a growing cell directly and consider whole-cell aspects such as dilution by growth. However, they contain only very few reactions. Larger cell processes are lumped into a few global variables, ignoring the details. This makes it impossible to model the effects of detailed changes (e.g. the fate of individual metabolites, the inhibition of single proteins, or deletions of single genes).

So why don't we just combine the two approaches? In Chapter 6 we already did this, to an extent: we considered a metabolic (FBA-like) model to compute biomass/enzyme productivity, and then used a formula derived from bacterial growth laws (see Chapter 8) to convert biomass/enzyme productivity into cell growth rate. But there were still some disadvantages, most importantly, that the biomass composition has to be known. Biological facts such as biomass composition should not just come as empirical knowledge, but should be explained *ab initio*. Subdividing a system into metabolism and other systems (with "fixed interfaces", like a predefined biomass composition) makes it impossible to understand details of complex behavior. For example, if some proteins are needed to retrieve iron via metabolism and other proteins (or possibly, the same proteins) also contain iron, this will create a complicated protein allocation response to changing iron supply. Iron depletion is expected lead to an upregulation of the first type of proteins, and a downregulation of the second type. However, FBA can only capture the first effect, because its biomass function is fixed.

### 9.1.2. Linearization: replacing kinetics by catalytic constraints

So how can large-scale models be made tractable? If we radically linearize all formulae, then instead of a biconvex or convex/concave problem, we obtain a linear problem (a bit like FBA); more precisely, a system of linear equalities and inequalities that define a set of feasible states. This set is a polytope, and linear optimality problems on this set can be solved easily. More specifically, to model metabolism in a growing cell, we need to consider dilution in the growing cell volume. If compounds are described by concentration, dilution can be effectively modeled by "consuming reactions", with a flux given by $v^{\text{dil}} = \mu c$, the compound concentration multiplied by the growth rate. By adding these hypothetical dilution reactions to the metabolic network, we obtain a new stationarity condition $N v = \mu c$ that connects the vectors of fluxes and compound concentrations, and in which the growth rate $\mu$ appears as a parameter. For each choice of the parameter $\mu$, we can ask whether a feasible steady growth state – i.e. a feasible combination of $v$ and $c$ exists, Furthermore, the feasible combinations $(\mu, v, c)$ form a convex set, with possible solutions $(v, c)$ for low values of $\mu$ and no solutions above a critical value $\mu^{\text{max}}$. Finding this critical value – the maximal possible growth rate for our model – as well as the corresponding optimal fluxes $v$ and compound concentrations $c$ is relatively easy, and can be done by solving a series of Linear Programming problems (checking for potential solutions $(v, c)$ for different values of $\mu$).

The main downside of this approach is that all relationships between models variables need to be linearized. This concerns, most importantly, all catalyzed processes, in which we assume a linear dependence between catalyzed flux and catalyst (enzyme or machine) concentration, but ignore the dependence on the concentrations of substrates, products,
<figure><img src="image_79.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>

the enzyme and flux levels (with the factors $κ_i$ that have become constants now we have set the internal metabolite concentrations), and thus the extreme rays of the enzyme and flux space will be equal and EFMs, as pointed out above (see also Chapter 5 and Figure 7.1). Optimization in this space will lead to the optimal flux distributions following an EFM (see Box 7.B for the detailed calculations). As fixing part of the optimal solution should lead to the same optimal solution, this required the flux distribution of the optimization over all variables to follow an EFM, as was indeed the case.

We point out two important aspects, using the network (Figure 7.2) as an example. First, it is convenient to

<table><caption>Table 7.1: Outcomes of the optimization of the example network with standard kinetics, parameter values and external concentrations (see Box 7.A) 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>
cofactors, or additional regulators. What does this mean in practice? As we know from Chapter 3, typical enzymatic rate laws have the form $v = e k(c)$: the rate $v$ is proportional to enzyme level $e$ and enzyme efficiency $k$, which is given by a kinetic rate law $k(c)$, a nonlinear function of the metabolite concentrations. Depending on the context, $k$ is also called catalytic rate or apparent $k_{cat}$. The kinetic rate laws $k(c)$ have typical shapes, as described in Chapter 3. To linearize the expression for $v$, while keeping the dependence on $e$, we need to replace $k$ by a fixed number, a model parameter. If the metabolite concentrations were known (experimentally, or from kinetic models under optimality assumptions, see Chapter 6), the value of $k$ could be computed. Otherwise, it can also be determined experimentally, by measuring $v$ and $e$ and setting $k = v/e$. Obviously, in reality, neither $c$ nor $k$ will be fixed and given, but for our linearized model, we need to assume this. This holds both for metabolic reactions (with enzymes as catalysts) and for macromolecular reactions (with molecular machines as catalysts). Under this assumption, we can replace all kinetic constraints by two linear constraints on the enzyme to fulfill its function as $-e k' \le v \le e k$, where coefficients $k$ and $k'$ are the approximation coefficients of the enzyme kinetics in the forward and backward direction respectively. These constraints, which replace the kinetic rate laws, are called capacity constraints and resemble the linear flux/enzyme relation in FBA with molecular crowding. Finally, there is one important exception: like in satFBA. some kinetic rate laws can remain, namely those in which the metabolite concentrations themselves are not variables (or disregarded completely), but are given as model parameters. This concerns, in particular, all transport reactions whose substrate is an extracellular compound with a known concentration in the growth medium. Nonlinearities can be kept as long as they do not depend on decision variables (even if the nonlinear problem can still sometimes be solved for a few nonlinearities if convexity is kept).

### 9.1.3. Towards fine-grained resource allocation models

Some FBA variants account for extra empirical constraints on the total concentration of metabolic enzymes (FBA with molecular crowding, or FBAwMC), or on proteome sectors (Constrained-Allocation FBA, or CAFBA). While these can predict metabolic states more reliably, the empirical constraints come as model assumptions and thus cannot be understood by the models themselves. In these models, the primary assumption is that the cell phenotype is obtained by genetic regulations, and the main goal and utility of genetic regulation can be interpreted as ways of saving resources. Thus in many cases when we predict cell phenotype maximizing growth, we find predictions in good agreement with the experimental observations.

To obtain an ideal detailed cell model, the self-replicator models from the previous chapter would have to be extended with explicit kinetics and thermodynamic constraints. Since we know that such problems are hard to solve, simplifications are currently necessary to keep linearity (and convexity) to solve optimization problems for large-scale models.

## 9.2. Additional constraints in resource allocation models

The canonical FBA/genome-scale metabolic models (GEMs), as described above (and in previous chapters), hold a number of assumptions which allow metabolic states irregardless of the corresponding resource allocation pattern. Therefore, we need to impose sets of additional constraints on metabolic networks. The fine-grained models build on the genome-scale metabolic models to encompass all the reactions that potentially could happen in a metabolic network, and different extensions of GEMs were proposed to predict optimal resource allocation in different microorganisms [266].

Resource allocation models describe optimal resource allocation in constraint-based, genome-scale whole-cell models and formalizes the mathematical relationships defining the interactions and allocation of resources between the cellular processes. All these relationships take the form of linear, growth-rate dependent equalities and inequalities, and form a convex feasibility problem [267, 268, 269]. We will now expand on the three types of constraints used in resource allocation models; note that the description is not exhaustive and peculiarities may vary among different formulations.
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{Gex}]/K_{\text{Gex}}}{1 + [\text{G}]/K_{\text{G}} + [\text{Gex}]/K_{\text{Gex}}} \cdot (1 - e^{\Delta_r G'_{0}/RT}) \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 (1 - e^{\Delta_r G'_{1}/RT}) \nonumber \\
v_2 &= e_2 \cdot k_{\text{cat},2}^+ \cdot \frac{[\text{P}]/K_{\text{P}}}{1 + [\text{Pex}]/K_{\text{Pex}} + [\text{P}]/K_{\text{P}}} \cdot (1 - e^{\Delta_r G'_{2}/RT}) \nonumber \\
v_3 &= e_3 \cdot k_{\text{cat},3}^+ \cdot \frac{([\text{P}]/K_{\text{P}})([\text{ADP}]/K_{\text{ADP}})([\text{O}_2]/K_{\text{O}_2})}{1 + ([\text{CO}_2]/K_{\text{CO}_2})([\text{ATP}]/K_{\text{ATP}}) + ([\text{P}]/K_{\text{P}})([\text{ADP}]/K_{\text{ADP}})([\text{O}_2]/K_{\text{O}_2})} \cdot (1 - e^{\Delta_r G'_{3}/RT}) \tag{7.3} \\
v_4 &= e_4 \cdot k_{\text{cat},4}^+ \cdot \frac{[\text{Pex}]/K_{\text{Pex}}}{1 + [\text{Pex}]/K_{\text{Pex}} + [\text{P}]/K_{\text{P}}} \cdot (1 - e^{\Delta_r G'_{4}/RT}) \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 (1 - e^{\Delta_r G'_{5}/RT}) \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$ are the same reaction, but defined in the opposite direction. The standard set of parameters we used for the toy model is all $k_{\text{cat},i}^+ = 10 \text{ s}^{-1}$ except $k_{\text{cat},3}^+ = 0.1 \text{ s}^{-1}$, all $\Delta_r G'_{i}/RT = -440$ and all $K_{\text{M}} = 1 \text{ mM}$. For the external metabolites $[\text{Pex}] = 1 \text{ mM}$, $[\text{Gex}] = 0.05 \text{ mM}$, $[\text{O}_2] = 0.1 \text{ mM}$, $[\text{BM}] = 1 \text{ mM}$ and $[\text{CO}_2] = 10 \text{ mM}$ unless mentioned otherwise.

split reversible reactions such that fluxes are always positive. In this case, that means that the reversible reaction from P to P<sub>ex</sub> is split into the forward reaction v<sub>2</sub> and the reverse reaction v<sub>4</sub>, both of which 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 convenience; we still assume this to be one reaction in the biological sense, and therefore the kinetic equations of both the forward and the backward reactions will be exactly the same. Depending on in which direction the flux goes, either one of the reactions will be positive and the other zero. 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 to avoid this. 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>ex</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 Calculation of optimal states

We can now use the result that states of maximal enzyme efficiency are reached at an elementary flux mode to calculate optimal states in a metabolic network using the following steps:

1. Enumerate the elementary flux modes that include the objective flux

2. Calculate the minimal enzyme for each EFM scaled to an objective flux of 1

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

Step 1 is possible for relatively large networks, although usually not for genome scale metabolic networks.

Step 2 is a convex optimization problem as we have seen in Chapter 6 and Step 3 is straightforward.

These three steps together are called Enzyme Flux Cost Minimization, because it is similar to Enzyme Cost Minimization, but while that is focused on fixed fluxes, Enzyme Flux Cost Minimization simultaneously finds the optimal fluxes, enzyme and metabolite levels. In this section we will show the method on the example network of Figure 7.2.
Constraining metabolic models to represent limited metabolic capacity of cells is not a new idea. There are two ways to approach this budgeting problem: "protein budgeting", where a fixed amount of protein needs to be partitioned in the optimal manner (maximizing growth), and "resource budgeting", where models include both the protein budgeting and the descriptions of demands for protein synthesis. However, "protein budgeting" problems assume that investments in protein production follow the budget, and not vice versa. Therefore, resource allocation models extend and embed the ideas of proteome partitioning beyond frameworks like CAFBA [270] and GECKO [271], or representing metabolic capacity limitations beyond FBAwMC [272].

Currently, there are three main implementations of large-scale resource allocation models: Resource Balance Analysis (RBA) [273], Models of Metabolism and Macromolecular Expression (ME-models) [274] and proteome-constrained models (pc-models) [59]. All these implementations are formalized as LP feasibility problems at fixed growth rate.

Originally, ME-models were considered as an extension of M-models, by including predictions for mRNA, protein, and ribosome levels. Importantly, they do not consider density constraints that, for instance, RBA does. Therefore, limitations on the capacity of exchange fluxes (as in FBA) are necessary to obtain a solution.

### 9.2.1. Types of constraints

The general description of these constraints in fact is the same as for small, coarse-grained self-replicator models, only the number of individual constraints increases. Although the precise formulations vary, resource allocation models build on three principal types of constraints:

* Mass-conservation constraints

* Flux coupling constraints

* Compartment capacity constraints

Alongside these three major classes there is another set of constraints, which we could call "environment" constraints - these correspond to, e.g. the composition of growth medium, biomass composition at a given growth rate $\mu$, etc. They are implemented by setting target values for amounts and/or fluxes defining a viable cell in a given (or several) environmental conditions, but they are not structural constraints. These constraints usually are added *ad hoc* and do not need to bear any functional meaning *per se*.

The *mass-conservation* constraints define the metabolic network (stoichiometry and relation between fluxes). The initial building blocks of these extended models are GEMs, and thus the metabolic network stoichiometry is already there; what remains to be defined are the protein turnover processes. We consider 4 types of protein turnover reactions in fine-grained resource allocation models: protein synthesis, folding, degradation and dilution-by-growth. So, for every protein present in such a model, we add these four reactions: two of them, translation and degradation, include the stoichiometry of amino acids needed for its translation and released during degradation based on the protein sequence. Protein folding and dilution are modeled by the same types of reactions for different proteins (a series of two steps, "unfolded" → "folded" → ∅).

Next, the *flux coupling* constraints couple the metabolic fluxes with protein usage: usually, the usage scales with the catalytic turnover value $k_{\text{cat}}$ of the enzyme. In this step we have to collect the kinetic information (in most cases, $k_{\text{cat}}$ values), which are used as model input. We establish the coupling between fluxes and protein synthesis by setting $v = k_{\text{cat}} e^{-\eta}$, where $e$ is the enzyme concentration and $0 < \eta \le 1$ is an efficiency term summarizing the effects of reaction thermodynamics, enzyme saturation, and possibly small-molecule regulation. Coupling constraints are introduced to couple both (i) the metabolic reactions with enzyme usage (as described above) and (ii) protein turnover reactions with the respective macromolecular machinery (e.g. sum demand of ribosomes for protein translation, $v_{\text{translation}} = [\text{Ribosome}] \times k_{\text{cat,ribosome}}$). The sheer number of the kinetic parameters needed for formulating the coupling constraints in the fine-grained models requires the modeller to consider different assumptions and simplifications when building and parametrizing these models, as briefly discussed below.
Box 7.B Optimal metabolic states in the example network

We minimize the enzyme investment for $v_{BM} = 1$ with $[P_{ex}] = 0$ (and therefore $v_4 = 0$ and EFM3 is not feasible) for the network in Figure 7.2 (the optimization problem in Equation (7.1)). Assuming all $h_i = 1$, the objective function $\sum_{i=1}^{r} h_i e_i = e_0 + e_1 + e_2 + e_3 + e_{BM}$. The constraints $v_{BM} = 1$ and $e, s \ge 0$ in Eq. (7.1) are straightforward. The steady state of all internal metabolites (G, P, ADP and ATP) leads to the following equalities (the steady states of ADP and ATP lead to the same equality):

$$ \text{Steady state ATP} \implies 100 v_{BM} = 2 v_1 + 10 v_3 $$

$$ \text{Steady state P} \implies 2 v_1 + v_4 = v_2 + v_3 + v_{BM} $$

$$ \text{Steady state G} \implies v_0 = v_1 $$

Substituting $v_{BM} = 1$ and $v_4 = 0$ and solving this set of linear equations, we can write all fluxes as functions of $v_2$: $v_0 = v_1 = 5 + \frac{5}{11}v_2$ and $v_3 = 9 - \frac{1}{11}v_2$ (there is only one independent flux in this system). This means we can draw the feasible flux space on the $v_2$ line and we can express the objective function in terms of $v_2$:

$$ \begin{aligned} \sum_{i=1}^{r} h_i e_i &= e_0 + e_1 + e_2 + e_3 + e_{BM} \\ &= v_0/\kappa_0 + v_1/\kappa_1 + v_2/\kappa_2 + v_3/\kappa_3 + v_{BM}/\kappa_{BM} \\ &= (5+5/11v_2)/\kappa_0 + (5+5/11v_2)/\kappa_1 + v_2/\kappa_2 + (9-1/11v_2)/\kappa_3 + 1/\kappa_{BM} \quad (7.5) \\ &= \underbrace{(5/\kappa_0 + 5/\kappa_1 + 9/\kappa_3 + 1/\kappa_{BM})}_{\alpha} + \underbrace{(5/(11\kappa_0) + 5/(11\kappa_1) + 1/\kappa_2 - 1/(11\kappa_3))}_{\beta} v_2 \\ &= \alpha + \beta v_2 \end{aligned} $$

The kinetic functions ($\kappa_i$) depend on several parameters (external metabolite levels [$G_{ex}$], $O_2$, [CO$_2$] and [P$_{ex}$], catalytic constants, Michaelis constants and Gibbs free energies) and the variables [G], [P], [ATP] and [ADP]. That means that once we have a set of internal metabolite concentrations $s$, the enzyme levels in the objective function can be written as a constant times the flux: $e_i = v_i/\kappa_i$, with $\kappa_i$ a constant. For a set of parameters, $\alpha$ and $\beta$ are positive or negative depending on the choice of $s$. It is clear that when we minimize this objective function by adjusting $v_2$, we will always have an optimum at $v_2 = 0$ (when $\beta$ is positive) or $v_2 = 99$ (when $\beta$ is negative). $v_2 = 99$ is the maximum of $v_2$ because then $v_3 = 9 - \frac{1}{11}v_2 = 0$, and higher values of $v_2$ would lead to negative values for $v_3$.

In conclusion, the optimum cannot be at a value of $0 < v_2 < 99$. If there would be an optimum with $0 < v_2 < 99$, we can determine $s$ and calculate whether $\beta > 0$ to find a lower objective value at $v_2 = 0$ or $v_2 = 99$, contradicting that we started with an optimum. Only if $\beta = 0$ there is a range of optima, but this requires very precise parameter values. $v_2 = 0$ and $v_2 = 99$ correspond to EFMs of this network (Figure 7.2).

First, we describe the network with the stoichiometric matrix ($N$) and the concentration vector ($s$):

$$ 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}, \quad s = \begin{pmatrix} [G] \\ [P] \\ [ATP] \\ [ADP] \end{pmatrix} \quad (7.6) $$

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

$$ \frac{d}{dt} s = N 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} \quad (7.7) $$

Now we find the EFMs (for example with EFMtool [194]). It can easily be checked that the following EFMs
The number of processes described in a fine-grained manner directly translates to the number of reactions and metabolites in the model. For instance, transcription is modelled explicitly in the ME-models [274]. The modellers' decision is key here: under assumption that transcription and translation form a linear pathway with fixed scaling factors (i.e. there is a fixed ratio of peptides translated per mRNA transcribed), the flux through mRNA translation reaction can be computed post-optimization based on the flux through the protein translation reaction. Explicit modelling of transcription would require describing processes of mRNA transcription, processing, export from nucleus, and then cytosolic degradation after the mRNA is translated – for each of the transcripts, with precise stoichiometry and a new set of coupling constraints.

The next issue is kinetic parametrization of these fine-grained models. We currently can use only very simplified kinetics in the models (flux coupling $v = k_{cat} e^{\eta}$), and simplify such factors as enzyme saturation and thermodynamic driving force into a single value of factor $\eta$. Two approaches are used to deal with this, as a large fraction of parameters are not even available. First, condition-dependent kinetic parameters ("apparent catalytic constants", $k_{app}$) are fitted from experimental (mostly quantitative proteomics) data (setting $k_{eff} = k_{cat}$ $\alpha$, where $0 < \alpha \le 1$) with a value $\alpha$ chosen to match predicted enzyme abundance and experimental measurements. Otherwise, for the enzymes with measured $k_{cat}$ values, we can assume that enzymes work at their maximal rate, i.e. the saturation function $\eta = 1$. Then the model computes the minimal protein requirement to sustain the flux through the metabolic reactions. The comparison of minimal predicted vs. observed protein abundance can represent the "apparent saturation", or "overcapacity" of enzymes. For instance, it is common in yeast *S. cerevisiae* that the flux and not protein expression varies across conditions, and the relationship between predicted and measured expression can suggest the nature of the observed protein expression [275].

The final layer of information in the fine-grained resource allocation models is a set of capacity constraints. These constraints describe the [upper] limit of cellular process(es), e.g. maximal protein capacity of a compartment. The capacity constraints are formulated as weighted sums of protein abundance, and usual weighing multipliers are proportional to the molecular weight of the protein. Usually, the capacity constraints are expressed in terms of mass, area, and volume of the compartment (e.g. "what is the maximal mass the mitochondrial proteins can take up in gDW of cells?"). Based on the biological interpretation of the constraints, we formulate the weighing multipliers to represent either of the metrics (mass/area/volume) that every protein occupies.

The capacity constraints can be both equality and inequality constraints: more frequent are the latter (usually defining the "upper limit" of, e.g. amount of protein targeted to mitochondria). However, some cell properties should be described through equality constraints: one of these is the protein density of biomass, defining the "target" protein translation per gram dry cell biomass.

### Box 9.B: Protein abundance versus concentration in fine-grained resource allocation models

Here we would like to include a relevant note for interpretation of the output of the fine-grained resource allocation models. Both the classical FBA and these extensions do not consider "metabolite concentration" as a concept: optimization variables are all *fluxes*. Frameworks discussed in this chapter model protein synthesis from amino acids and energy equivalents explicitly, with a typical flux dimension of *mmol gDW⁻¹ h⁻¹* (as for any other fluxes). To compute the amount of protein that has to be produced in the steady-state growth, we should consider the flux balance for the protein *e*: $v_{synthesis,e} = v_{degradation,e} + v_{dilution,e}$, or, rewritten with the respective parameters, $v_{synthesis,e} = (k_{deg,e} + \mu) e$. Here, $k_{deg,e}$ is the degradation rate for the protein *e*, and $\mu$ is the specific growth (= dilution-by-growth) rate. The [*e*] in the rewritten equation holds dimension of *mmol gDW⁻¹*, which is protein abundance, rather than concentration.

The predicted amount of protein in cells can be compared to experimental measurements in two ways. First option is to convert abundance to concentration using the relationship between the cell volume and dry weight (e.g. $V_{gDW} = 1.7 \text{ mL gDW}^{-1}$ in *Saccharomyces cerevisiae*, [8]). Alternatively, proteome mass fractions are a popular unit in label-free mass spectrometry-based protein quantification, a popular method in quantitative microbiology. Respectively, predicted proteome mass fractions can be inferred by converting protein abundance in *mmol* to *g*,
<figure><img src="image_80.png" /><figcaption>Figure 7.2: States of maximal efficiency in an 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. At a very low concentration of extracellular glucose ([Gex]), EFM3 has the lowest cost. But as we move along the x-axis, at around [Gex] = 0.02 there is a switch to EFM1 and later, at around [Gex] = 0.07, EFM2 becomes the one with the lowest cost. (D) Specific fluxes (flux divided by total enzyme) associated with the optimal EFM for different levels of Gex. Note that v1 is not shown as it is always equal to v0. The rates show a discontinuity when there is a switch from one optimal EFM to another.</figcaption></figure>

(denoted by vectors $f^{(i)}$) are in the nullspace of the stoichiometric matrix:

$$f^{(1)} = \begin{pmatrix} 5 \\ 5 \\ 0 \\ 9 \\ 0 \\ 1 \end{pmatrix}, \quad f^{(2)} = \begin{pmatrix} 50 \\ 50 \\ 99 \\ 0 \\ 0 \\ 1 \end{pmatrix}, \quad f^{(3)} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 10 \\ 11 \\ 1 \end{pmatrix} \qquad (7.8)$$
and scaling to the protein content in dry cell biomass. Here, it is important to consider the conversion factors (protein content in dry biomass). Fast-growing bacteria, e.g. *E. coli*, maintains rather constant protein content in dry weight across growth rates (ca. 0.55 (*g* protein) *g*DW⁻¹) [276, 277]. On the contrary, the protein content is known to vary in *S. cerevisiae* as a function of growth rate [8].

### 9.2.2. Interpretation and outcomes of the additional constraints

For cells growing in exponential phase at a growth rate $\mu$, the constraints listed above formalize different structural constraints:

(I) the metabolic network has to produce all metabolic precursors necessary for biomass production and mass conservation must hold for all intracellular molecule species - i.e. intracellular metabolites and molecular machines.

(II) the capacity of each type of molecular machine must be sufficient to ensure its function, i.e. to catalyze chemical conversions at a sufficient rate;

(III) the intracellular density of compartments and the occupancy of membranes must not exceed the defined limits.

Metabolism and protein production depend on each other in both ways!

* Possible states are defined by constraints:

(i) stationarity, including dilution

(ii) kinetics/capacity constraints

(iii) density constraints

(iv) other physiological bounds ("targets")

* Select an optimal state by an assumed objective (e.g. maximal growth; protein production; in chemostat, growth at a given growth rate and at minimal external substrate level)

* A change in one cell component may lead to small perturbations or global rearrangements

* Relaxing one constraint can lead to changes in the entire system!

## 9.3. Resource balance analysis

Resource Balance Analysis (RBA) has been developed as (and is considered) a flexible and generic modeling framework which describes the functioning of an organism using the most relevant set of linear equality and inequality constraints. Importantly, RBA does not refer to a *specific* unique set of constraints, but it provides the framework in which such constraints can be taken into account. Here we distinguish the general RBA framework from a specific type of RBA model, i.e. the model of a specific organism, which contains by definition a set of specific constraints. The general structure of such models, as well as the typical constraints, are shown in Figure 9.1.

### 9.3.1. Mathematical description

**Notation.** Below $A^T$ refers to the transpose of the matrix $A$. $\mathbb{R}_{>0}^n \triangleq \{x \in \mathbb{R}^n | x_i > 0 \text{ for all } i \in \{1, \dots, n\}\}$, $\mathbb{R}_{>0} \triangleq \mathbb{R}_{>0}^1$, $\mathbb{R}_{\ge 0}^n \triangleq \{x \in \mathbb{R}^n | x_i \ge 0 \text{ for all } i \in \{1, \dots, n\}\}$ and $\mathbb{R}_{\ge 0} \triangleq \mathbb{R}_{\ge 0}^1$.

In a standard RBA model, we consider balanced growth, that is, the average state of a cell in a cell bacterial population growing exponentially at the specific (constant) growth rate $\mu \ge 0$, i.e. the amount of produced biomass per biomass per cell per unit of time. Our simulated average cell is composed of different molecule species:

(i) $n_y$ types of molecular machines, which can be subdivided further into $n_e$ enzymes and transporters involved in the metabolic network $E \triangleq (E_1, \dots, E_{n_e})$ at the concentrations $e \triangleq (e_1, \dots, e_{n_e})^T$ and metabolic fluxes
<figure><img src="image_81.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 per amount of enzymes.</figcaption></figure>

The next step is to perform the convex optimization over the metabolite levels for each one of the three EFMs. Therefore, we express the enzyme levels as a ratio of the flux and the function $f(s)$, using Equation 7.4. Summing over all enzymes, we get a function for the total enzyme cost (level) as a function of fluxes, metabolite concentrations and parameters:

$$e_{\text{tot}} = \sum_{i} e_{i} = \sum_{i} \frac{v_{i}}{\kappa_{i}(\text{s})}. \qquad (7.9)$$

We use the standard parameters (Box 7.A) and replace $v_i$ by the values given by each EFM. We are then left with a convex optimization over the metabolite levels, an Enzyme Cost Minimization problem as in Chapter 6. For $[G_{ex}] = 0.05$ we obtain a total enzyme of 111.1 for EFM1, of 146.3 for EFM2 and 136.5 for EFM3. That means that for these conditions we will conclude that EFM1 is optimal. From the optimization we obtain the metabolite concentrations: $[G] = 0.08$, $[P] = 3.93$, $[ATP] = 0.11$ and $[ADP] = 20.09$ (that the internal glucose concentration is higher than the external is because we described the transport with regular enzyme kinetics instead of transporter enzyme kinetics, which would have been more realistic). 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 Enzyme Flux Cost Minimization is one way of explaining those shifts.

Above, Enzyme Cost Flux Minimization was used to find the metabolic state with the maximum enzyme efficiency. Although in our calculation we obtain the enzyme concentrations last, it is by enzyme concentrations that cells actually control metabolism. If cells produce enzymes at the concentrations we calculated and reach a steady state, this state will realize the fluxes and metabolite levels that lead to our optimal state.

## 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's growth rate. 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 plot the "growth rate" of the different EFMs instead of "overall enzyme efficiency".
<figure><img src="image_65.png" /><figcaption>Figure 9.1: Resource balance analysis model: overview of biological components and mathematical constraints. (a) Typically, an RBA model describes metabolisms and macromolecule production in a growing cell (yellow blocks). Precursors from metabolism are needed to produce macromolecules, and some macromolecules serve as enzymes to catalyze metabolic reactions. In addition, macromolecules are diluted and are localized in cell compartments. (b) Mathematical constraints. The variables and processes described by an RBA model must satisfy a number of constraints, include mass-balance constraints (between production, degradation, and dilution of compounds); capacity constraints (relating process velocities to the concentrations of catalysts); density constraints (on the total amount of compounds in a cell compartment); and possibly empirical physiological constraints on any types of "target variables", to ensure realistic models.</figcaption></figure>

$\nu \triangleq (\nu_1, \dots, \nu_{n_e})^T$; and $n_m$ macromolecular machines $M \triangleq (M_1, \dots, M_{n_m})$ involved in non-metabolic cellular processes, such as the translation apparatus, at the concentrations $\mathbf{m} \triangleq (m_1, \dots, m_{n_m})^T$;

(ii) $n_p$ proteins $\mathbb{P} \triangleq \{P_1, \dots, P_{n_p}\}$ belonging to unspecified cellular processes. $\mathbf{p} \triangleq (p_1, \dots, p_{n_p})^T$ denotes the set of concentrations of $\mathbb{P}$;

(iii) $n_s$ intracellular and mass-balanced metabolites $\mathbb{S} \triangleq (S_1, \dots, S_{n_s})$. Within the set $\mathbb{S}$, we distinguish a subset $\mathbb{B} \triangleq (B_1, \dots, B_{n_b})$ of abundant metabolites which have fixed growth-independent concentrations $\bar{\mathbf{b}} \triangleq (\bar{b}_1, \dots, \bar{b}_{n_b})^T$ (and usually coincide with biomass macro-components such as DNA, cell wall or plasmic membrane). We also consider a set of extracellular metabolites $\mathbb{S}_{\text{ext}} \triangleq (S_{\text{ext},1}, \dots, S_{\text{ext},n_{\text{ext}}})$ of concentrations $\mathbf{s}_{\text{ext}} \triangleq (s_{\text{ext},1}, \dots, s_{\text{ext},n_{\text{ext}}})^T$ that are not mass-balanced.

Finally, let us introduce the vector $\mathbf{y}^T \triangleq (\mathbf{e}^T, \mathbf{m}^T)$ of concentrations of molecular machines of size $n_y$. Typical units of concentrations e, m and p are in millimoles per gram of cell dry weight, and fluxes $\nu$ in millimoles per gram of cell dry weight per unit of time.

For a given cell growth rate $\mu \ge 0$, the RBA optimization problem (named $\mathcal{P}_{\text{rba}}(\mu)$) can be formalized mathematically as follows.
<figure><img src="image_82.png" /><figcaption>Figure 7.4: Optimal growth rates of the two EFMs for different levels of the external metabolite Gex, computed using Equation (7.10) from the enzyme demands (at a unit biomass production rate) shown in Figure 7.2 (C).</figcaption></figure>

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 per total enzyme demand" would give us directly the growth rate. However, biomass does not only consist of metabolic enzymes, but includes ribosomal enzymes, RNA, DNA, lipids, and other compounds. Therefore we need a formula for converting biomass/enzyme efficiency into cellular growth rate.

Mathematically, a cell's growth rate is given by $\mu = v_{BM}/s_{BM}$, where $v_{BM}$ is the biomass production rate (biomass produced per cell volume and time) and $s_{BM}$ is the biomass amount per cell volume. If a cell contained nothing but metabolic enzymes (more precisely, the enzymes described in our model), the biomass/enzyme efficiency $\kappa_{BM} = v_{BM}/h_{enz}$ would directly describe the cellular growth rate. Since that is not the case, we need to convert $h_{enz}$ to $s_{BM}$. The metabolic protein fraction decreases with the growth rate, leading to a hyperbolic dependency of the growth rate on the biomass production rate (Figure 7.3). We may use the empirical approximation $h_{enz}/s_{BM} = \kappa_{prot}(a - b\mu)$, where $\kappa_{prot} = 0.5$ is the fraction of protein mass within the cell dry mass and the parameters $a = 0.27$ and $b = 0.2$ h were fitted to describe the metabolic enzyme fraction in proteomics data, assuming a linear dependence on growth rate [19]. This yields the conversion formula (see also [191]):

$$ \mu = \frac{a \kappa_{prot} v_{BM}}{h_{enz} + b \kappa_{prot} v_{BM}}. \qquad (7.10) $$

This formula has been used to convert the minimal enzyme cost per biomass flux for different external concentrations in the toy model (Figure 7.2c) to the maximal growth for each EFM (Figure 7.4).

## 7.6 Application to central metabolism in E. coli bacteria

In the previous sections, we saw that finding enzyme-efficient metabolic states can be done by iterating through all possible EFMs and performing the enzyme cost minimization on each one. We demonstrated it on a toy model comprising only 3 EFMs. In Wortel et al. [191], this method was scaled up and applied to a more realistic model covering the central metabolic network, as shown in Figure 7.5A. For this larger network, there are 1566 biomass-generating EFMs. Each reaction is assigned to a single enzyme along with its molecular weight, $k_{cat}^{+}$, $K_{M}$, and $\Delta_{r}G^{\circ}$, and follows the generalized factorized rate law as in Equation (7.2). These parameters are listed in Appendix section B.1, and the full procedure for obtaining them is described in Wortel et al. [191], along with other model parameters.
For a fixed vector of concentrations $\mathbf{p} \in \mathbb{R}_{>0}^{N_{n_p}}$ and the given growth rate $\mu \ge 0$,

$$
\begin{array}{ll}
\text{find possible cell states} & \mathbf{y} \in \mathbb{R}_{\ge 0}^{n_y}, \nu \in \mathbb{R}^{n_e}, \\
\text{subject to} & (C_1) \qquad -\Omega\nu + \mu(\mathbf{C}_Y^S \mathbf{y} + \mathbf{C}_B^S \bar{\mathbf{b}} + \mathbf{C}_P^S \mathbf{p}) = 0 \\
& (C_{2a}) \qquad \mu(\mathbf{C}_Y^M \mathbf{y} + \mathbf{C}_P^M \mathbf{p}) - \mathbf{K}_T \mathbf{y} \le 0 \\
& (C_{2b}) \qquad -\mathbf{K}'_E \mathbf{y} \le \nu \le \mathbf{K}_E \mathbf{y} \\
& (C_3) \qquad \mathbf{C}_Y^D \mathbf{y} + \mathbf{C}_P^D \mathbf{p} - \bar{\mathbf{d}} \le 0
\end{array}
$$

where all the inequalities are defined componentwise and:

* $\Omega$ is the stoichiometry matrix of the metabolic network of size $n_s \times n_e$, where $\Omega_{ij}$ corresponds to the stoichiometry of metabolite $S_i$ in the j-th enzymatic reaction;

* $\mathbf{C}_Y^S$ (resp. $\mathbf{C}_P^S$) is an $n_s$ $n_y$ (resp. $n_s$ $n_p$) matrix where each coefficient $\mathbf{C}_{Y_{ij}}^S$ corresponds to the number of metabolite $S_i$ consumed (or produced) for the synthesis of one machine $Y_j$ (resp. $P_j$); $\mathbf{C}_{Y_{ij}}^S$ is then positive, negative or null if $S_i$ is produced, consumed or not involved in the synthesis of one machine $Y_j$ (resp. $P_j$);

* $\mathbf{C}_B^S$ is an $n_s \times n_b$ matrix in which each coefficient $\mathbf{C}_{B_{ij}}^S$ corresponds to a metabolite $S_i$ consumed (or produced) for the synthesis of one $B_j$;

* $\mathbf{K}_T$ ($\mathbf{K}_E$ and $\mathbf{K}'_E$, respectively) are matrices of size $n_m \times n_y$ ($n_e \times n_y$, respectively) in which each coefficient $k_{T_i}$ ($k_{E_i}$ and $k'_{E_i}$, respectively) is positive and corresponds to the efficiency of molecular machine $M_i$, i.e. the rate of the process per amount of the catalyzing molecular machine, (the efficiency of the enzyme $E_i$ in forward and backward sense, respectively);

* $\mathbf{C}_Y^M$ (resp. $\mathbf{C}_P^M$) is an $n_m \times n_y$ (resp. $n_m \times n_p$) matrix in which each coefficient $\mathbf{C}_{Y_{ij}}^M$ typically corresponds to the length in amino acids of the machine $Y_j$ (resp. $P_j$). In some cases (for instance for the constraints on protein chaperoning), the length in amino acids can be multiplied by a coefficient, such as the fraction of the whole proteome that necessitates chaperoning;

* $\bar{\mathbf{d}}$ is a vector of size $n_c$, where $n_c$ is the number of compartments (compartment membrane and/or compartment interior for which density constraints are considered. $\bar{\mathbf{d}}^i$ is the density of molecular entities within the volume or surface area. Densities are typically expressed as a number of amino-acid residues by volume or surface area.

* $\mathbf{C}_Y^D$ (resp. $\mathbf{C}_P^D$) is an $n_c \times N_y$ (resp. $n_c \times n_p$) matrix in which each coefficient $\mathbf{C}_{Y_{ij}}^D$ corresponds to the density of one machine $Y_j$ (resp. $P_j$) in the compartment $i$. By construction, we have one unique localization per machine.

At given $\mu$, $P_{rba}(\mu)$ is a feasibility optimization problem, where constraints $(C_1-C_3)$ define the feasibility domain. The feasibility domain can be empty or non-empty. If there exists a solution $(\mathbf{y}, \nu)$ to $P_{rba}(\mu)$ -i.e. the feasibility domain is non-empty-, then there exists a feasible resource distribution compatible with the given growth rate. In other words, the cell can grow at this growth rate value.

Remarks:

1. In practice, the vector $\bar{\mathbf{b}}$ contains non-zero values only for the concentrations of macro-components such as DNA, cell wall, and lipid membranes, and for a few set of metabolites. These values are usually extracted from the biomass formation reaction used in FBA models (see Chapter 5).

2. To model reversible enzymes, we introduced two diagonal matrices containing the enzyme efficiencies, i.e. $\mathbf{K}_E$ and $\mathbf{K}'_E$, describing the capacity constraints of enzymes in both directions. If an enzyme $E_i$ is considered irreversible, $k'_{E_i}$ is set to 0.

3. In [268, 278], an RBA model was built for *Bacillus subtilis*. It integrates two macromolecular processes in constraint $C_{2a}$, the translation and chaperoning of proteins, and two density constraints, the limitation of the cytosolic density and of the membrane occupancy. An RBA model can be refined by integrating for instance other cellular processes and molecular machines, such as the transcription machinery, the protein secretion apparatus (see [278, 279]), or molecule turnover [280], as well as other types of constraints.
<figure><img src="image_83.png" /><figcaption>Figure 7.5: Model of central metabolism in *E. coli* bacteria. (A) The metabolic network of the *E. coli* model used by Wortel et al. [191]. Note that only for the purpose of visualization, the network shown here has been condensed by lumping consecutive reactions that are fully coupled (e.g., the reactions between DHAP and PEP are now represented by a single arrow). Furthermore, some groups of metabolites have been merged to a single node: H6P – representing the hexose phosphates G6P, F6P, and FBP; T3P – representing the triose phosphates G3P and DHAP; P5P – representing the pentose phosphates R5P, X5P, and Ru5P. The metabolites that are direct substrates of the biomass reaction are marked in bold. (B) A Venn diagram showing statistics of biomass-producing EFMs in the model and their reliance on oxygen.</figcaption></figure>

First, Wortel et al. [191] wanted to study the effect of environmental conditions on the growth rate of *E. coli*, and see whether the model would be able to recapitulate empirical phenomena. The external glucose concentration was set to 100 mM and oxygen levels were varied between 1 µ and 10 mM. They selected 4 flux modes as representatives (the EFMs *max-gr*, *ana-lac*, and *aero-ace*, and *exp*, which is based on experimentally measured fluxes; the flux distributions are shown in Figure 7.6), and calculated their predicted growth rates in each condition, using Equation (7.10). The results are shown in Figure 7.7. When focusing on a single flux mode, one can see that as the oxygen level increases so does the growth rate. The increase saturates at some point, which depends on the flux modes and on the kinetic parameters in the model. Indeed, it has long been known that growth rate dependence on a limiting substrate concentration has this specific shape – a relationship generally called the *Monod curve*.

In this specific example, it is interesting to see the Monod curves of the different EFMs, and try to understand the differences. First, the EFM called *ana-lac* (red curve), is a flat line. This makes sense because cells that use this EFM do not utilize the oxidative phosphorylation system and therefore do not require oxygen at all for growth. *max-gr*, on the other hand, is very sensitive to the level of oxygen mainly because
**How to incorporate the medium composition.** We represent the medium composition by extracellular concentrations, which determine the efficiencies of metabolic transporters via Michaelis-Menten like rate laws (as nonlinear $k(c)$ functions; see section 9.1.2). For an extracellular nutrient $S_{ext,i}$ with concentration $s_{ext,i} \ge 0$, the efficiency of the corresponding metabolic transporter(s) is given by $k_E(s_{ext,i}) = \frac{k_{cat}s_{ext,i}}{K_m+s_{ext,i}}$, with parameters $k_{cat}$ and $K_m$ for the turnover number and the affinity of the transporter, respectively.

**Convexity.** For given growth rate and medium composition, all equalities and inequalities in our RBA problem $P_{rba}(\mu)$ is linear in the decision variables $(y, v)$ and is proven to be convex [267, 269]. By construction, the feasibility domain of $P_{rba}(\mu)$ corresponds to the set of all possible phenotypes of the cell at a growth rate $\mu \ge 0$.

**Parsimonious resource allocation.** The principle of parsimonious resource allocation between cellular processes coincides with principles of optimal performance, for example, that a cell phenotype should maximize growth rate. Importantly, for an RBA problem with given parameters, there exists a maximal growth rate $\mu^* \ge 0$, such that for any $\mu$, $P_{rba}(\mu)$ is feasible if and only if $\mu \le \mu^*$ [267, 269]. For a given medium composition, the maximal growth rate $\mu^*$ can be computed by using a bisection algorithm, in which a series of LP problems are solved to narrow down the exact growth rate at which the problem becomes infeasible. Together with the maximal feasible growth rate one obtains the optimal cell configuration maximizing growth $(\mu^*, y^*, v^*)$.

**Exploration of the feasibility domain.** It was proven that the feasibility domain contracts with increasing growth rate [267, 269]. In the same vein as Flux Variability Analysis ([281], see Chapter 5), the feasibility domain can be explored at optimal $(\mu^*)$ or sub-optimal $(\mu \le \mu^*)$ growth rates. For one decision variable $y_i$ (resp. $v_i$), two LP problems are solved, where (i) constraints $C_1, C_2$ and $C_3$ remain unchanged; (ii) the decision variable $y_i$ (resp. $v_i$) is maximized (LP 1) and minimized (LP 2). This operation is repeated for each decision variable to obtain *in fine* the feasibility domain of all decision variables. In practice, at the optimum, the cell configuration $(\mu^*, y^*, v^*)$ is often unique. Indeed, non-unique solutions will exist if two alternative metabolic pathways have exactly the same cost in resources. Since all enzymes have different amino acid sequences, use different cofactors, are differently localized, etc, this is highly unlikely. A caricatural example of a model with non-unique solutions would be one in which an enzyme pool is arbitrarily split into two, and the two new "enzyme species" are given different names, although they are physically exactly the same.

### 9.3.2. How to build and calibrate RBA models?

For a given organism, an RBA model can include all known metabolic reactions and all (potentially relevant) cell processes with major protein investments (including in particular metabolism, transport, ribosomes, chaperones). Which metabolic reactions and cell processes are regarded as relevant may vary between organisms and is a modeling choice. RBA integrates the description of the stoichiometry of the metabolic network, the composition of some macromolecular components as DNA, mRNA, lipids, cell wall components, but also a detailed description of molecular machines including the precise protein and nucleic sequences, their dependency on metallic ions, vitamins, cofactors for activity, and their cellular localization.

**Requirement for RBA model creation.** Several pieces of information, that is, several types of biological data, are needed to build an RBA model for an organism:

* Metabolic network and reaction stoichiometry

* Mapping from reactions to enzymes and further to proteins

* Protein amino acid sequences and cofactors (e.g. metal ions)

* Efficiencies of metabolic enzymes,

* Types of cellular machines (e.g. ribosomes and chaperones), their efficiencies and cofactor demands

* Molecule sizes and localisations (for density constraints)

* Other constraints on concentrations or fluxes ("targets") defining a viable cell

In practice, RBA models can be built and simulated using a software called RBApy [279]. It takes as input (i) a genome-scale metabolic network in SBML format [282] and having a high-quality mapping between genes and reactions; (ii) the
<figure><img src="image_84.png" /></figure>

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

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

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

Figure 7.6: The four flux modes chosen for drawing the Monod curves in Figure 7.7: (A) *max-gr* – the EFM with the highest growth rate under the standard conditions chosen in the study, (B) *aero-ace* – an EFM which mixes between respiration and acetate fermentation, (C) *ana-lac* – an EFM that does not require oxygen (i.e. anaerobic) and uses lactate fermentation, and (D) *exp* – which is not an elementary flux mode, but rather one based on the measured flux distribution for *E. coli* growing on minimal media and glucose. The active reactions are highlighted in color (with the flux direction indicated by the arrowhead). The magnitude of each flux is not shown here but can be found in the Supplementary section of Wortel et al. [191]. The biomass reaction is not shown here due to space limitations, but is always active.

of the high flux going through oxidative phosphorylation. It is also the EFM with the highest growth rate
in standard oxygen levels (0.21 mM), even when taking all the other ~1500 EFMs into account (not shown
here).

Instead of screening only external oxygen levels, we can also screen several model parameters and compute "winning EFMs", their enzyme demands, and the resulting growth rates for our parameter combination. By screening glucose and oxygen concentrations, we obtain the Monod landscape shown in Figure
description of molecular machines involved in macromolecular processes (e.g ribosomes); (iii) the NCBI taxon of the
organism of interest. Then it automatically extracts information on protein sequences and cofactors, subunit structures,
localization from public repositories such as Uniprot, maps enzymes to the reactions they catalyze and to the proteins
they consist of, and builds a draft uncalibrated RBA model.

**Calibration of RBA model parameters.** An RBA model may contain a high number of model parameters. How can we obtain them? Remember that the main parameters to be estimated are the abundance of housekeeping proteins (p), densities (d̄) and efficiencies of molecular machines ($K_E, K'_E, K_T)^1$. As we learned in Chapter 3, the rate of an enzymatic reaction, v = e kapp, depends on the enzyme's efficiency or "apparent catalytic rate", given by kapp = f(c) = kcat + ηrev(c) · ηsat(c) < kcat. The kapp values are always below the kcat value, but may vary from state to state depending on metabolite concentrations. Since internal metabolite concentrations c are unknown and difficult to measure at genome-scale, we cannot estimate kapp from the explicit kinetic law f(c). We need to obtain these kapp parameters empirically, for example by measuring the flux v and the protein abundance e in one condition and taking their ratio. Hence, for a given environmental condition, RBA parameters can be estimated using quantitative proteomics [283] for p and d̄, and a combination of quantitative proteomics and fluxomics (if available) as in [278] or FBA to estimate the flux distribution as in [279] for efficiency parameters. To account for variable enzyme efficiencies, one may make the simplifying assumption that enzyme efficiencies depend mostly on growth rate. By estimating the enzyme efficiencies at different growth rates and interpolating between them, one obtains empirical relationships between them and the growth rate [278] that can be further used in Prba(μ). For instance, several estimates of enzymatic efficiencies obtained in several contrasted growth conditions will provide a relationship such as KE(μ) instead of a constant KE value.

9.3.3. Examples of RBA applications

Despite its assumptions and simplifications, RBA is able to predict realistic complex adaptations (e.g. some proteins contain iron; under iron starvation, the cell may *increase* the import of iron, but also *avoid* using proteins that contain iron, and the pathways in which they operate), even in complex growth conditions [278, 284]. These effects cannot be captured by standard CBM methods in an autonomous way, i.e. without the addition of specific constraints on fluxes. Thus altogether, RBA provides a reasonable compromise between numerical tractability, cell phenotype prediction accuracy and model complexity.

9.4. Prediction of variable biomass composition

Cell models describe, among other things, what a cell is composed of (see Chapter 2). In such models, cell composition
may either serve as input information (e.g. putting a certain amount of lipids as a constraint) or model results
(to be validated by data). Cell composition may be described broadly by molecule classes (e.g. the overall amounts of
lipids, protein, DNA, RNA, and so on), or more precisely for individual molecule types (individual lipids, proteins, mRNA
species, and so on), or even in terms of atomic composition (which in turn gives clues about the amounts of molecule
classes). Another measure of biomass composition is the amounts of precursors that are made in metabolism in order to
produce macromolecules. In FBA, specifically, "biomass" refers to the proportions of these precursor amounts (including
cofactors), assuming that the total production may change, but the proportions stay fixed.

Cell composition may greatly vary between cell types, and also for a given cell type between conditions (for instance, because of the different ribosome content of a cell at different growth rates). This variable composition often poses a challenge for models. Just like the uptake rates, the varying biomass composition reflects complex global rearrangements of resources, and choices between metabolic strategies. In classical FBA, uptake rates and biomass composition need to be predefined; if they are not known from experiments, predictions may not be reliable (see Chapter 5).

<footnote n="1">¹$\bar{b}$ is deduced from the biomass reaction of the genome-scale metabolic network of the organism</footnote>
<figure><img src="image_88.png" /><figcaption>Figure 7.7: Monod curves (cell growth rate as a function of oxygen level) computed using the model shown in Figure 7.5 – Each curve was computed using one of the EFMs and the associated (oxygen-dependent) enzyme demands. The ana-lac strategy (anaerobic growth with lactate secretion) does not use oxygen, therefore it curve is flat.</figcaption></figure>

7.8 (A). Just like in Figure 7.7, there are distinct parameter regions in which optimal growth is reached by specific EFMs. While the *max-gr* EFM remains best when glucose and oxygen levels are high, at low oxygen levels we see a large number of different EFMs, one of them *ana-lac* (see the EFM phase diagram in Figure 7.8 (B)). More results for this model (fluxes plotted in the EFM phase diagram and in flux space, as well as ideal and real enzyme costs for all EFMs), are shown in Appendix Section B.2.

## 7.7 Concluding remarks

In this chapter we considered the metabolic network of a cell - and enzyme levels, metabolite concentrations, and fluxes as the state variables - and studied its maximally efficient states. Finding such states can be difficult because fluxes, metabolite concentrations, and enzyme levels are tightly coupled: metabolite concentrations determine enzyme efficiencies, enzyme efficiencies determine optimal enzyme levels, and enzyme levels determine fluxes and metabolite concentrations, which in turn determine enzyme efficiencies. To find an optimal state, all variables need to be optimized at the same time, which is a non-linear optimality problem with (possibly) many local optima. In small toy models, solutions can be found numerically, but for large detailed models, the computational effort becomes enormous. Instead of simplifying the problem (as in the previous chapters) we here used the insight that (in models without extra flux bounds) the optimal solutions must be EFMs.

### Recommended readings

* M.T. Wortel, H. Peters, J. Hulshof, B. Teusink, and F.J. Bruggeman. Metabolic states with maximal specific rate carry flux through an elementary flux mode. FEBS Journal, 281(6):1547-1555, 2014.

* S. Müller, G. Regensburger, and R. Steuer. Enzyme allocation problems in kinetic metabolic networks: Optimal solutions are elementary flux modes. Journal of Theoretical Biology, 347:182-190, 2014.

* M.T. Wortel, E. Noor, M. Ferris, F.J. Bruggeman, and W. Liebermeister. Metabolic enzyme cost explains variable trade-offs between microbial growth rate and yield. PLoS Computational Biology, 14 (2):e1006010, 2018.

### Problems

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

#### Problem 7.1 Effect of oxygen concentration
A main achievement of RBA models is that the biomass composition is not fixed and given as input information as in FBA, but is predicted by the model itself. The model predicts the optimal macromolecule composition of the cell and – unlike in small self-replicator models (see the models in Chapter 8) – this prediction is very detailed: it follows from the predicted proteome and may comprise details like trace elements (e.g. iron, which appears in some of the proteins) and whose contribution to the overall biomass will vary with the expression level of those proteins. How does this work? In RBA models, in contrast to FBA, uptake rates and biomass composition are usually not predefined (although this is possible, if required, and may be used to fit models to data), but obtained as model predictions from resource allocation ideas. For example, in RBA, the proportions of biomass precursors would reflect the types and amounts of macromolecules (to be made from these precursors) in the optimal state, which in turn (among other things) reflect the amounts of enzymes needed to make the precursors.

## Recommended readings

* Website rba.inrae.fr for further details on RBA

* A recent review on large-scale resource allocation models: K. de Becker et al. "Using resource constraints derived from genomic and proteomic data in metabolic network models" Curr Opin Syst Biol 2022, 29:100400

## Problems

### Problem 9.1 The role of metabolite concentrations

* What if the metabolite content of the cell has been underestimated? Assume that the amount of small metabolites in cells is currently underestimated. What problems in model predictions would arise from the fact? In what way would predictions (by FBA or other methods) be distorted?

* In what way would a cell, in reality, profit from a lower small metabolite content? Can we assume that the ratio between small metabolites and proteins is optimized? Describe possible aspects of this compromise! For inspiration, see [285].
<figure><img src="image_89.png" /></figure>

Figure 7.8: Monod landscape – (A) Similar to the 1-dimensional Monod curve (Figure 7.7), the graphics shows the cell growth rate as a function of external glucose and oxygen concentrations, predicted from the *E. coli* model in Figure 7.5. The growth rate of the “winning” EFM – i.e. the one with the highest growth rate under the glucose and oxygen levels matching the *x* and *y* values – determines the height of each point. Each color represents the region in which a certain EFM is the “winning” one. (B) EFM phase diagram. The same plot as in (A), seen from above. The “winning EFMs” form a sort of phase diagram. At the boundary between every two regions, the two EFMs lead to the same growth rate (similar to the intersections between curves in Figure 7.7). The EFMs from Figure 7.6 are marked by their names. Note that the colors in this figure do not match the previous colors marking these select EFMs.

Consider the model in Figure 7.2. What would be the qualitative effect of a change in oxygen concentration on the enzyme cost of the three EFMs and on the choice of the optimal strategy?

**Problem 7.2 Effect of external metabolites**

Consider the model in Figure 7.2 under standard conditions (Box 7.A and [$G_{ex}$] = 1, such that EFM2 is optimal, and EFM1 second best (remember that the higher the enzyme cost, the less optimal the EFM). What might happen when we gradually increase the concentration [P<sub>ex</sub>]? What is the qualitative effect on the enzyme cost of the three EFMs?

**Problem 7.3 States of maximal growth rate**

Consider the following small toy network:

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

We want to optimize the specific pathway flux for the production of P (which is v<sub>3</sub>/e<sub>tot</sub>, where we assume all enzymes to have equal costs: i.e. e<sub>tot</sub> = e<sub>1</sub> + e<sub>2</sub> + e<sub>3</sub>) at steady state. We assume mass-action kinetics, meaning the rate is the enzyme concentration times the forward rate constant times the substrate minus the backward rate constant times the product: v = e(k<sup>+</sup>s - k<sup>-</sup>p). Unless mentioned otherwise, we use the values s<sub>2</sub> = 10, k<sub>1</sub><sup>+</sup> = 2, k<sub>1</sub><sup>-</sup> = 1, k<sub>2</sub><sup>+</sup> = 3, k<sub>2</sub><sup>-</sup> = 1, k<sub>3</sub><sup>+</sup> = 1, k<sub>3</sub><sup>-</sup> = 0.1, p = 0 (concentrations are denoted by lower case letters).

(a) Write out the rate equations for all three rates in terms of the parameters and the concentrations.

(b) Give an expression of the total enzyme concentration in terms of fluxes and the metabolite concentrations s<sub>1</sub> and x.
(c) Find the concentration of X for which the specific flux $v_3/e_{tot}$ is maximal for $e_1 = 0$ and $s_2 = 10$, and also give the corresponding value of $v_3$. HINT: Is it easiest to set $e_{tot} = 1$ and maximize $v_3$, replacing $e_3$ using the equation for the total enzyme cost and the steady state assumption.

(d) Find the concentration of X for which the specific flux $v_3/e_T$ is maximal for $e_2 = 0$ and $s_1 = 10$, and also give the corresponding value of $v_3/e_{tot}$.

(e) Find the concentration of X for which the enzyme cost is minimal for $e_1 = e_2$ and $s_1 = s_2 = 10$, and also give the corresponding value of $v_3/e_{tot}$.

(f) What was the best distribution of enzymes from the three options above for $s_1 = 10$?

(g) Find the concentration of X for which the enzyme cost is minimal for $e_1 = 0$ and $s_1 = 50$, and also give the corresponding value of $v_3/e_{tot}$.

(h) Find the concentration of X for which the enzyme cost is minimal for $e_2 = 0$ and $s_1 = 50$, and also give the corresponding value of $v_3/e_{tot}$.

(i) Find the concentration of X for which the enzyme cost is minimal for $e_1 = e_2$ and $s_1 = 50$, and also give the corresponding value of $v_3/e_{tot}$.

(j) What was the best distribution of enzymes from the three options above for $s_1 = 50$?

(k) Interpret the results from this problem in light of the proof shown in this chapter about the optimal specific flux being attained at an EFM.
# Chapter 10

## Optimal cell behavior in time

Hidde de Jong, Dafni Giannari, Diego A. Oyarzún, Steffen Waldherr, and Agustín G. Yabo

### Chapter highlights

* Microorganisms live in continually changing environments, which require them to develop adaptation strategies.

* These strategies have been profitably studied under the assumption that microorganisms have evolved to optimize one or several aspects of their adaptive response.

* The mathematical formalization of this assumption leads to dynamic optimization problems that can be solved by means of techniques from optimal control theory.

* The chapter discusses three example problems: dynamic optimization of enzyme expression in metabolic pathways, dynamic optimization of coarse-grained models of cellular growth, and dynamic flux balance analysis.

* The results obtained for these problems illustrate the interest of studying adaptation strategies from the perspective of dynamic optimization, and the strengths and weaknesses of this approach.

### 10.1. Introduction

The study of microorganisms in the laboratory has often focused on the creation of stable conditions enabling balanced, reproducible growth of the population. Such conditions are almost never found in nature. Microorganisms live in continually changing environments in which nutrients are only intermittently available and in which the cells are submitted to a variety of other temporally varying stresses (acidity, heat and cold, drought, ...). In order to survive in these conditions, microorganisms have developed a range of molecular mechanisms to detect changes in the environment, or signals announcing such changes, and to adapt their functioning accordingly.

A well-studied example of the dynamic response of bacteria to changes in their environment is the phenomenon of diauxic growth, discovered by Jacques Monod ([286] (see also Chapter 8). When *Escherichia coli* is grown in a medium containing a mixture of two carbon sources, e.g., glucose and lactose, the cells generally first deplete the carbon source supporting the highest growth rate (glucose) before starting to assimilate the other carbon source (lactose). A variety of mechanisms are involved in this switch from a preferred to a secondary carbon source, including the release of the repression of enzymes necessary for lactose utilization, the release of the inhibition of lactose transporters, and the global regulation of a large number of other genes [287, 288].

In many situations, the precise functioning of the molecular mechanisms regulating the adaptation of microbial physiology to changes in the environment is not or only qualitatively understood. This precludes their inclusion in quantitative models that accurately predict the dynamic response of the cell in a variety of conditions. The lack of mechanistic, quantitative
# Chapter 8

## Principles of cell growth

Ohad Golan, Hollie J. Hindley, Hidde de Jong, Markus Köbis, Elena Pascual Garcia, and Andrea Weiße

### Chapter overview

* A comprehensive description of fundamental growth laws in microbial growth, elucidating the core principles that govern biological growth patterns.

* A detailed exploration of the contrasts between coarse-grained and fine-grained modeling is presented, offering insights into the varying levels of detail that each approach encompasses.

* A thorough breakdown of the key assumptions in the modeling of metabolic systems is provided, underlining the foundational premises that are crucial for accurately representing these complex systems.

* The process of deriving fundamental growth laws by modeling key assumptions is meticulously demonstrated, enabling a clear understanding of how theoretical constructs translate into biological realities.

### 8.1 Introduction

A key feature of living systems is that they are able to grow and reproduce. The reproductive success in a given environment defines the fitness of a living system. The study of the growth of bacteria and other microorganisms is crucial for better understanding their capacity to cause diseases in humans or for better exploiting their use in biotechnological or environmental processes. Beyond their interest for a variety of applications, bacteria and other microorganisms have shown themselves ideal model systems for investigating fundamental questions on the relation between growth, fitness and characteristics of the environment.

One of the first to systematically and quantitatively study the growth of bacterial cultures was Jacques Monod in the 1940s. He performed so-called diauxic growth experiments, in which bacteria were cultured in a medium containing two different limiting carbon sources. He showed that the bacteria first deplete one carbon source before starting to assimilate the second carbon source. The order in which the primary and secondary carbon source were consumed was determined by the growth rate they support: the preferred carbon source allows the culture to grow at a higher rate. Further work on the molecular basis of diauxic growth led to the discovery that cells inhibit the expression and activity of functions for the use of secondary carbon sources when a preferred carbon source is present, a global regulatory mechanism known as carbon catabolite repression [195, 196].

Monod characterized bacterial growth by means of batch culture experiments in a well-defined growth medium allowing bacteria to reach a state of balanced growth, where the accumulation of biomass can be described by a single constant, the exponential growth rate. Together with the chemostat, a device
information can be bypassed by making appropriate assumptions about the regulatory systems, in particular that the
latter have evolved under the selection pressure of the environment to optimize the response to external perturbations.
More precisely, it is assumed that microorganisms have developed mechanisms that allocate limiting resources (proteins,
fluxes, ...) to cellular processes so as to maximize or minimize some objective function, or combination of objective
functions, over the time-interval of environment changes.

The use of an optimality assumption to make up for missing or incomplete information was already exploited with success in Chapter 5 of this book. The difference with those approaches is that here we are interested in cases where the optimality criterion is defined over an interval of time rather than at steady state, and thus we need to consider dynamic instead of static optimization. Moreover, some methods take into account that cells may vary the allocation of limiting resources to cellular processes over the time interval in which the environmental changes occur, instead of only considering a constant response in a stable environment. This generalization of the problem enormously increases its complexity. It may also lead to nontrivial dynamical effects that are not found in the case of static optimization, such as the accumulation of resource buffers to anticipate future changes in the environment [289].

The classical argument motivating the optimality assumption in the case of microorganisms is that mutants of genes coding for enzymes in central metabolism often have a lower growth rate than the wild-type strain, where growth rate is interpreted as indicating fitness [290]. This argument, however, derives from observations of balanced growth in a stable environment. Is there any evidence that, in the case of changing environments, microorganisms have evolved to perform dynamic optimization? Some circumstantial evidence is provided by the observed capacity of microorganisms to anticipate changes in their environment. For example, when moving along the digestive tract, *E. coli* cells are exposed first to lactose and then to maltose, thus requiring the ability to switch from growth on lactose to growth on maltose (reminiscent of diauxic growth in the laboratory) [291]. Interestingly, reporter gene studies found that the enzymes required for maltose assimilation are expressed at a much higher level in the presence than in the absence of lactose, in otherwise identical conditions [292]. This suggests a specific effect of the presence of lactose on the expression of maltose enzymes, preparing the cells for the expected future availability of maltose. This and other examples of anticipatory behavior are not conclusive in themselves, but they suggest that dynamic optimization is a plausible working hypothesis that may be useful in practice.

The aim of this chapter is to show how microbial physiology can be studied by means of dynamic optimization, by combining a specific objective function, or combination of objective functions, with models of different scope and granularity, while taking into account a number of biophysical and biochemical constraints. We first provide a general definition of dynamic optimization problems in the mathematical framework of optimal control. We then instantiate this general definition for three types of biological problems, each giving rise to a specific class of models. In particular, we discuss (i) dynamic optimization of enzyme expression in metabolic pathways, (ii) dynamic optimization of resource allocation in coarse-grained models of cellular growth, and (iii) dynamic flux balance analysis (dFBA) of metabolic networks. Across the different examples, the scope of the models varies from metabolic pathways (i) to metabolic networks (iii) to the entire cell (ii). The increase in scope is sometimes traded against a lower granularity of the description of cellular process (ii). Some of the models provide a kinetic description of the rates of the individual reactions (i and ii), whereas other models only provide constraints on the reaction rates (iii). In every case, different objective functions are tried, for example the minimization of the time to produce a given compound or the maximization of the amount of biomass produced.

For each of the biological problems and corresponding models considered, we give the precise definition of the modeling formalism and the optimization problem, a small example as an illustration, a discussion of the solution of the problem, and a brief description of more realistic applications and the insights they have given into the functioning of cellular networks. The chapter does not give a detailed explanation of the mathematical methods that are used for solving different classes of optimal control problems, because this would require knowledge of specialized mathematical concepts with which the average reader of the book may not be familiar. Moreover, these methods have been the subject of dedicated textbooks [293, 294]. Rather, we focus on the definition of the dynamic optimization problems and the interpretation of the solutions returned by available numerical solvers of optimal control problems.
<figure><img src="image_91.png" /><figcaption>Figure 8.1: Bacterial growth laws – (A) Monod growth law: growth rate dependency on nutrient availability (data from [200]). (B) Correlation between growth rate and nutrient uptake rate (data from [206]). (C) Correlation between growth rate and cellular composition (data from [204]).</figcaption></figure>

allowing continuous culture of microorganisms at a predefined growth rate [197], these methods have become standard in microbial physiology. They notably underlie the discovery of a number of so-called growth laws, relating the growth rate to a variety of properties of the physiology of growing bacteria. The growth laws are conserved across different organisms and a broad range of experimental conditions. Here, we list three well-known growth laws [198, 199]:

1. *Dependency of the growth rate on nutrient availability* [200]: In his characterization of bacterial growth, Monod discovered the first growth law. He observed that the growth rate of bacteria depends upon the nutrient concentration in the medium in a hyperbolic fashion (Fig. 8.1A).

2. *Correlation between growth rate and nutrient uptake rate* [201]: In continuous cultures, the growth rate was shown to vary linearly with the nutrient uptake rate (Fig. 8.1B). The slope of this linear relation is called the biomass yield and the offset the ‘maintenance energy’, as it is assumed to be derived from the energy spent on processes required to maintain the basic processes of the cell, in the absence of growth [202].

3. *Correlation between growth rate and cellular composition* [203, 204]: In 1959, Schaechter, Maaløe and Kjeldgaard showed that RNA, DNA and the number of nuclei in *Salmonella typhimurium* linearly correlate with the growth rate. Later, it was further shown that other physiological parameters, such as the mass fraction of ribosomes in growing populations, also linearly correlate with the growth rate [204] (Fig. 8.1C). Initially, it was believed that the correlation between ribosomal mass fraction and growth was strictly positive, however, Scott et al. [205] showed that when growth is inhibited through translation-inhibiting drugs, growth rate and ribosomal mass fraction exhibit a negative (near-)linear relation.

The conserved nature of the growth laws has led scientists to ask whether there are fundamental principles governing bacterial growth. To answer this question, different types of mathematical models have been developed. One approach aims at integrating all known molecular constituents of the cell and the reactions involving these constituents into a big model, an *in-silico* copy, or ‘digital twin’, of the cell. Such models, known as fine-grained models, can be useful to predict emergent phenotypes, but they are difficult to construct and maintain, and their complexity makes it hard to grasp certain principles that underpin growth. In this chapter, we will focus on coarse-grained models of bacterial growth. Rather than assembling individual reactions in a bottom-up manner, these models are based on the top-down definition of a limited number of basic cellular functions or processes involved in growth, described by appropriate macro-reactions (Fig. 8.2). Coarse-grained models are smaller and therefore easier to construct and analyze. The lack of molecular detail can make their predictions less accurate, but their simplicity allows a focus on how basic cellular functions and their interactions shape bacterial growth.

How much detail is included in a model depends on the specific scientific question asked, and similarly, models may vary in their underlying assumptions. Oftentimes, assumptions are based on biochemical principles governing intracellular reactions, on physical limitations faced by cells, on optimality principles,
## 10.2. Mathematical formalization of dynamic optimization problems

The models of cellular processes considered in this chapter have the form of systems of ordinary differential equations (ODEs) (Chapter 3). Dynamic optimization problems for such systems take the form of so-called optimal control problems, which have their roots in physics and engineering [293, 294].

Let $\mathbf{x}(t)$ be the (time-varying) state of the dynamical system, typically concentrations of (intracellular or extracellular) metabolites or proteins, and let $\mathbf{f}(\cdot)$ describe the (linear or nonlinear) dynamics of the state. $\mathbf{u}(t)$ denotes the (time-varying) control variables, e.g., fluxes allocated to specific reactions or protein fractions allocated to specific enzymes. The time-points $0$ and $T > 0$ indicate the bounds of the interval over which the behavior of the system is optimized, with respect to an objective function $J$. The behavior of the system, given the control exerted by $\mathbf{u}(t)$, is subject to constraints $\mathbf{c}_1(\cdot)$ and $\mathbf{c}_2(\cdot)$ on the admissible control inputs at specific time-points $t$ or over the entire time-interval $[t_0, t_e]$, respectively. The constraints express physical limitations, such as the intracellular density of molecular constituents (Chapter 2), or biochemical limitations, such as the maximum protein synthesis rate. Combining the above elements, we obtain the following definition of dynamic optimization problems:

$$ \max_{u \in U} J(\mathbf{x}(t), \mathbf{u}(t), 0, T), \qquad (10.1) $$

such that

$$ \frac{d\mathbf{x}}{dt} = \mathbf{f}(\mathbf{x}(t), \mathbf{u}(t)), \quad \mathbf{x}(0) = \mathbf{x}_0, \qquad (10.2) $$

$$ 0 \ge \mathbf{c}_1(\mathbf{x}(t), \mathbf{u}(t)), \qquad (10.3) $$

$$ 0 \ge \mathbf{c}_2(\mathbf{x}(0), \mathbf{x}(T)). \qquad (10.4) $$

In summary, the problem consists in finding controls that, given the dynamics of the system, maximize the objective function and satisfy the constraints [295].

The above definition makes no specific assumptions about the dynamics of the system under consideration. Given that we deal with biochemical reaction systems, the dynamics can be refined to

$$ \frac{d\mathbf{x}}{dt} = \mathbf{N}\mathbf{v}(\mathbf{x}(t), \mathbf{u}(t)) - \mu(t)\mathbf{x}(t), \quad \mathbf{x}(0) = \mathbf{x}_0, \qquad (10.5) $$

where $\mathbf{N}$ represents the stoichiometry matrix and $\mu$ is the (time-varying) growth rate. The principles of describing the structure of biochemical reactions systems by means of a stoichiometry matrix were described in Chapter 3 above.

The problem definition assumes that there is only a single objective function to be optimized. This may not be appropriate, since microorganisms seem to optimize several criteria in parallel, for example growth rate and survival under stress [296]. In many situations, it is therefore more appropriate to generalize the above problem to the case where $\mathbf{J}(\cdot \cdot \cdot)$ represents a vector of $n$ objective functions $\mathbf{J} = [J_1, \dots, J_n]$. Thus generalized, the problem does not usually have a single solution, but rather an infinite set of solutions located on a so-called Pareto surface [297]. Solutions on the Pareto surface have the property that every alternative solution improving the performance with respect to some objective necessarily degrades the performance with respect to at least one of the other objectives. In the problems developed in the sections below, we principally consider optimality in the case of a single, possibly composite objective.

Many methods for solving optimal control problems (10.1)-(10.4) exist. While some optimal control problems can be solved analytically, most of the problems considered in the examples below require numerical approximations to be solved. All examples developed in the sections below have been solved by means of freely available solvers.
<figure><img src="image_92.png" /><figcaption>Figure 8.2: Coarse-grained modeling of cellular growth – Compared to genome-scale FBA and whole-cell models (Chapters 4, 5, 10, coarse-grained models zoom out of the molecular detail and focus on key processes.</figcaption></figure>

or on a combination of these.

In this chapter, we show how to understand and, ultimately, how to develop coarse-grained models of cellular growth. We present a number of coarse-grained models with increasing levels of granularity. The models have been chosen to also represent a variety of commonly used assumptions, for example, based on growth rate maximization or on phenomenological or mechanistic constraints. Despite these differences, however, models we discuss generally recover the basic growth laws, and we show how the latter can be derived from solving two of the simplest coarse-grained models. The goals of this chapter are:

1. To enable the reader to understand and analyze any model of microbial growth from the literature.

2. To enable the reader to develop their own coarse-grained model of a metabolic system that is directed at their specific scientific question.

3. To provide the reader with a new perspective on modeling of complex systems and specifically the biological cell.

## 8.2 Fundamental modeling assumptions of microbial growth

The models of microbial growth we consider here are based on fundamental assumptions that follow from biochemical and biophysical constraints. In this section, we discuss and mathematically define assumptions that are found, explicitly or implicitly, in most coarse-grained models of microbial growth. The assumptions are formulated in an abstract manner to hold for any self-replicating biological system, irrespective of the specifics of the underlying molecular mechanisms. In the next section, we use these assumptions to construct increasingly complex models of microbial growth and show how the latter can be used to derive the experimentally observed growth laws presented in the introduction of this chapter.

The growth of microorganisms consists of the uptake of nutrients from the environment and the conversion of these nutrients into new microbial cells through a number of coupled metabolic processes (Fig. 8.2). This description brings out the self-replicating or autocatalytic nature of microbial growth: cells transform nutrients from the environment into new cells. In what follows, we consider growth on the population level, that is, an increase in the total amount of cells or, equivalently in many situations, an increase of the biomass of the population. This leads to the well-known model of microbial growth, where the change in biomass over time is proportional to the amount of biomass (Fig. 8.3A):

$$ \frac{dB}{dt} = \lambda B, \qquad (8.1) $$

where $t$ [h] denotes time, $B$ in gram dry weight [gDW] the biomass and $\lambda$ [1/h] the population growth rate. If the growth rate is constant, the solution to Eq. (8.1) describes exponential growth of the biomass:

$$ B = B_0 e^{\lambda t}, \qquad (8.2) $$

where $B_0$ [gDW] is the initial biomass at $t = 0$.
10.3. Dynamic optimization of enzyme expression in metabolic pathways

A number of experimental works suggest that metabolic regulation encodes temporal patterns in enzyme expression that may be beneficial for cell physiology [298, 299]. Since the timing of gene expression can directly control resource expenditure, several authors have attempted to rationalize such patterns as solutions of optimal control problems defined as in (10.1)-(10.4). The general idea is to optimize the temporal evolution of enzyme concentrations using objective functions that are representative of cellular goals. This provides a rationale to reverse-engineer optimality principles that underlie the expression patterns observed in experiments. In this section, we briefly describe results obtained for unbranched metabolic pathways, the basic building blocks of the metabolic networks of the cell.

Dynamic optimization of enzymatic concentrations was first considered by Klipp and co-workers [300]. The problem under study was the minimal-time activation of an unbranched network from an “off” state, where only the precursor is present, to a state where all substrate has been converted into product. To this end, the authors considered an unbranched pathway with n enzymes and (n + 1) metabolites:

$$
\begin{align}
\frac{dx_0}{dt} &= -k_1 e_1 x_0, \nonumber \\
\frac{dx_i}{dt} &= k_i e_i x_{i-1} - k_{i+1} e_{i+1} x_i, \tag{10.6} \\
\frac{dx_n}{dt} &= k_n e_n x_{n-1}, \nonumber
\end{align}
$$

with a given initial condition $x_0(0) \neq 0$ and $x_i(0) = 0$ for $i = 1, 2, \dots, n$, and where all enzymatic reactions are assumed to follow mass-action kinetics with rate constant $k_i$. To model the “off” state prior to pathway activation, the initial conditions can be set to $x_0(0) = s$, where s is the concentration of precursor at $t = 0$, and $x_i(0) = e_i(0) = 0$ for all $i = 1, \dots, n$. The goal was to determine a vector of optimal enzyme concentrations $e(t)$ that solve the following problem:

$$
\mathbf{e}^*(t) = \arg \min_{e \in \mathcal{U}} \frac{1}{s} \int_{0}^{\infty} (s - x_n(t)) \, dt, \quad (10.7)
$$

subject to the dynamic model in (10.6) and constraint set U as in (10.1) defined by a limited overall enzyme abundance
over the optimization horizon:

$$
\sum_{i=1}^{n} e_i(t) = e_{\text{tot}}, \tag{10.8}
$$

where $e_{\text{tot}}$ is a constant amount of total enzyme concentration. The objective function in (10.7) is called the *transition time* of the pathway and quantifies the time needed to convert all precursor into product. Note that the optimization problem (10.6)-(10.8) falls within the general class of problems defined by (10.1)-(10.4).

Numerical solutions of the optimization problem reveal that the enzyme profiles have a temporal sequence that matches the order in which the enzymes appear in the pathway. Crucially, such pattern resembles the “just-in-time” strategies widely studied in operations research [301], whereby costly resources are deployed only when needed in a production line. In the context of cellular metabolism, such a strategy implies that minimal time activation tends to express biosynthetic enzymes only when their substrates have been built up to sufficiently high concentrations, and thus avoid wasteful protein expression.

The first experimental demonstration of the just-in-time principle was presented by Zaslaver and colleagues [298]. This work employed luminescent and fluorescent reporters to measure the temporal adaptation of *Escherichia coli* upon withdrawal of amino acids from the growth media. Clear just-in-time patterns of enzyme expression were found in the serine, methionine and arginine biosynthetic pathways. To better understand such patterns, the authors studied a model
<figure><img src="image_93.png" /><figcaption>Figure 8.3: Coarse-grained models of metabolic systems with increasing complexity – (A) A self-replicating system. (B) The simplest description of a metabolic system: coupled catabolic and anabolic reactions. (C) A metabolic system that can catabolize two different nutrient sources. (D) A catabolic system requiring two substrates to grow: x and E.</figcaption></figure>

The growth rate is a key parameter that is often used as a proxy for the fitness of microorganisms. It is dependent on the metabolic processes, that is, how a cell utilizes the nutrients to synthesize new biomass (self-replication). The simplest description of metabolism is that it takes up a nutrient, breaks it down into metabolites (catabolism), and then utilizes these metabolites to produce new biomass (anabolism) (Fig. 8.3B). Catabolic and anabolic processes comprise a variety of biochemical reactions that are carried out by different sets of proteins and enzymes. The reaction rates of these processes are limited biochemically and biophysically. We formulate these limitations as modeling assumptions and define them as mathematical constraints, four of which we briefly review below.

### 8.2.1 Conservation of mass and quasi-steady-state assumption

Dry biomass is often a more readily measurable quantity than cell volume. The latter relates absolute abundances of cell components to their intracellular concentrations. Yet, because bacterial cells have been observed to maintain approximately constant cell density across various growth conditions [207, 208] (though transient exceptions have been observed at the single-cell level [209]), biomass can be regarded a proxy for volume and is therefore assumed to be proportional to cell volume in many growth models. All models considered in this chapter are based on the assumption of constant cell density and approximate the concentration x of a cellular component x (we use normal font for cell components and italic font for their concentrations) by its absolute abundance divided by the cell mass.

According to the law of mass conservation, the change of mass is equal to the inflow minus the outflow of mass. As a consequence, the change in concentration of a cell component, for example a metabolite pool, is determined by the sum of the rates of the reactions consuming and producing this cell component (Fig. 8.4A). The mass balance for any cell component x is given by the following equation:

$$ \frac{dx}{dt} = \sum_{y} r_{y \to x} - \sum_{k} r_{x \to k}, \qquad (8.3) $$

where $r_{y \to x}$ denotes the rate of the reaction converting cell component y into cell component x (production of x), and $r_{x \to k}$ the rate of the reaction converting cell component x into cell component k (consumption of x). Typically, cell component concentrations have units mg/gDW or mmol/gDW, so that rates of metabolic reactions are expressed in units mg/(gDW h) or mmol/(gDW h).

In the simple system shown in Fig. 8.3B, there are two reactions: one converting the nutrient source N into a metabolite X and one utilizing the metabolite for the synthesis of biomass. According to (8.3), the flux balance of metabolite pool x is given by $dx/dt = r_{n \to x} - r_{x \to B}$.

A key assumption is that intracellular concentrations are in quasi-steady state. This means that cell component pools remain constant:

$$ \frac{dx}{dt} = 0, \text{ for all cell components } x. \qquad (8.4) $$
for an unbranched pathway with three enzymatic steps and Michaelis-Menten kinetics:

$$ \frac{dx_i}{dt} = k_{\text{cat},i} e_i \frac{x_{i-1}}{x_{i-1} + K_{M,i}} - k_{\text{cat},i+1} e_{i+1} \frac{x_i}{x_i + K_{M,i}} - \mu x_i, \quad i = 1, \dots, 3, \qquad (10.9) $$

with given initial conditions $x_1(0) \neq 0$, $x_2(0) = x_3(0) = 0$, and where ($k_{\text{cat},i}$, $K_{M,i}$) are the enzyme turnover rate and Michaelis-Menten constants of each enzyme, respectively. The precursor concentration $x_0$ is assumed to be constant. The model also includes a dilution term that accounts for dilution by cell growth at rate $\mu$. In contrast to previous works, this model also includes an explicit description of enzyme expression:

$$ \frac{de_i}{dt} = \frac{\beta_i}{1 + r/\kappa_i} - \mu e_i, \quad i = 1, \dots, 3, \qquad (10.10) $$

where the first term is a lumped model of enzyme expression controlled by a time-varying (active) repressor concentration $r(t)$, with maximal expression rate $\beta_i$, and $\kappa_i$ being the concentration of (active) repressor required for half-maximal expression. Moreover, since bacterial amino acid pathways are often subject to end-product feedback, the model assumed that the repressor gets activated by the pathway product:

$$ r(t) = r_{\text{T}}(t) \frac{x_3(t)}{K_r + x_3(t)}, \qquad (10.11) $$

where $r_{\text{T}}(t)$ denotes the total (active and inactive) repressor concentration. The model also included negative autoregulation of the repressor itself:

$$ \frac{dr_{\text{T}}}{dt} = \frac{\beta_0}{1 + r/\kappa_0} - \mu r_{\text{T}}, \qquad (10.12) $$

where $\beta_0$ and $\kappa_0$ define the strength of autoregulation similarly as in the lumped model for enzyme expression in (10.10).

The authors constructed an optimization problem so as to study the relation between optimality, and the strength of the regulatory parameters $\mathbf{k} = (k_1, k_2, k_3)$ and $\beta = (\beta_1, \beta_2, \beta_3)$. To this end, they defined the optimization problem

$$ \min_{\mathbf{k},\beta} a \cdot \sum_{i=1}^{3} \underbrace{\int_{0}^{T} \frac{\beta_0}{1+r(t)/\kappa_0} dt}_{\text{total amount of repressor}} + \underbrace{\int_{0}^{T} |F-F_{\text{goal}}| dt}_{\text{deviation from steady state}}, \qquad (10.13) $$

where $a$ is a scalar weight accounting for the protein costs, $T$ is the optimization horizon, and $F$ is the rate of product synthesis:

$$ F = k_{\text{cat},3} e_3 \frac{x_2}{x_2 + K_{M,3}}. \qquad (10.14) $$

In problem (10.13), the constant $F_{\text{goal}}$ is a prescribed production flux that the pathway should achieve at steady state. Minimization of the objective in (10.13) accounts for the activation of the pathway from an “off” state until it reaches a prescribed flux $F_{\text{goal}}$. This formulation differs from the previous example [300] in two important ways. First, it accounts for cellular resources in the objective function itself. The first term of the objective quantifies the total amount of repressor produced through the optimization horizon, and thus relates to the amount of cellular resources required to activate the pathway. Second, the decision variables are the regulatory parameters, not the temporal profiles of the molecular species. Therefore, strictly speaking, this is not an optimal control problem but rather a static optimization problem subject to dynamic constraints encapsulated by the pathway ODE model. Through numerical solutions for different values of the protein cost weight $a$ and optimization horizon $T$, the authors determined conditions under which the optimal solutions showed two features of the just-in-time property, namely:

$$ \tau_1 < \tau_2 < \tau_3, \quad \max_t e_1 > \max_t e_2 > \max_t e_3, \qquad (10.15) $$
The quasi-steady-state assumption simplifies the mathematical analysis of the system significantly and holds for balanced growth of the microbial population. In this chapter, we focus mostly on situations in which the quasi-steady-state assumption applies, but also give an example of a model with metabolic dynamics. In metabolic modeling, the rates of reactions at steady state are called fluxes, denoted by the symbol $J$. With the quasi-steady-state assumption, Eq. (8.3) becomes

$$ \sum_{y} J_{y \to x} = \sum_{k} J_{x \to k} \qquad (8.5) $$

that is, for every cell component, the sum of production fluxes equals the sum of consumption fluxes. In the example system, we have $J_{n \to x} = J_{x \to B}$.

### 8.2.2 Proteome allocation assumption

The biochemical reactions breaking down nutrients into intracellular metabolites, and the reactions utilizing these metabolites for the synthesis of new biomass, do not occur spontaneously. The reactions are catalyzed mostly by proteins complexes, in particular metabolic enzymes and ribosomes. In coarse-grained models, well-defined sets of biochemical reactions are grouped together into macro-reactions. The cell components that are necessary to catalyze the individual steps of a macro-reaction are grouped together into a corresponding so-called proteome sector. A proteome sector includes mostly proteins that catalyze metabolic reactions but also ribosomes catalyzing the reaction of protein biosynthesis. Proteins constitute most of the biomass of the cell [210]. Therefore, as a first approximation, the sum of the proteome sectors equals the total biomass of the growing population measured in units of $g$ (Fig. 8.4B):

$$ \sum_{r \in \{x \to y\}} P_r = B, \qquad (8.6) $$

where $P_{x \to y}$ is the proteome sector catalyzing the macro-reaction that transforms cell component x into cell component y. The proteome sectors as defined above are extensive quantities, summed over the entire growing population, like the total biomass $B$. For the models, we are rather interested in intensive quantities, the amount of a proteome sector relative to the total amount of biomass (protein), corresponding to protein concentrations or protein fractions. Dividing the left-hand and right-hand sides of Eq. (8.6) by $B$, we thus obtain:

$$ \sum_{r \in \{x \to y\}} p_r = 1 \qquad (8.7) $$

where $p_{x \to y}$ is the fraction of the proteome converting x into y, defined by $p_{x \to y} = P_{x \to y}/B$. Proteome fractions are dimensionless and sum to one.

In the simple example system in Fig. 8.3B, we distinguish two macro-reactions: a catabolic reaction and an anabolic reaction (biomass synthesis). We therefore define two proteome sectors, corresponding to enzymes and ribosomes, respectively, with fractions $p_{n \to x}$ and $p_{x \to B}$, respectively. In later examples in the chapter, the catabolic and anabolic processes are further broken down into smaller macro-reactions and so are the proteome sectors.

### 8.2.3 Mathematical description of reaction fluxes

The rate at which a reaction is converting one cell component, e.g., a metabolite, into another is determined by the proteome fraction, the concentrations of the substrates of the reaction and possible regulation by other cell components in the system. While mass-action kinetics provide a principled framework to develop rate equations for biochemical reactions, in practice, various approximations based on mechanistic assumptions are often used to obtain simplified equations [211]. Below there are a few examples of rate laws defining the fluxes in coarse-grained models:

1. Excess substrate and no allosteric interactions. The simplest relation of the flux $J$ to the relevant proteome
where $\tau_i$ is the response time, i.e. the time to reach 50% of maximal concentration, and $\max_t e_i$ is the peak concentration of each enzyme. This theoretical model was designed to mimic the architecture of gene regulation in such pathways, whereby the end product commonly represses the expression of upstream enzymes, and thus gave both experimental and computational evidence that just-in-time patterns may be the result of optimality principles underlying the regulation of metabolic pathways.

Further experimental evidence of temporal patterns in enzyme expression have been found in other pathways [302] and organisms [299], and number of subsequent works have explored their optimality in more detail; we refer the reader to the review in [303] for a detailed discussion on such approaches. Oyarzún and colleagues [304], in particular, gave the first mathematical proof that just-in-time dynamics are a general property in models of unbranched metabolic pathways. Using a cost-benefit objective function that balances the speed of response against the cost of expressing pathway enzymes, they showed that the just-in-time patterns emerge in pathways of arbitrary length and with minimal assumptions on the enzyme kinetics. Specifically, they considered a general model for an unbranched pathway with $n + 1$ reactions:

$$ \frac{dx_i}{dt} = g_{i-1}(x_{i-1}) e_{i-1} - g_i(x_i) e_i, \quad i = 1, \dots, n, \qquad (10.16) $$

with initial conditions $x_i(0) = 0$ for $i = 1, 2, \dots, n$, and the precursor $x_0$ assumed to be at a constant concentration. The functions $g_i$ represent a general kinetic turnover rate satisfying the following conditions:

$$ \begin{aligned} g_i(0) &= 0, \\ \frac{\partial g_i(x_i)}{\partial x_i} &> 0. \end{aligned} \qquad (10.17) $$

The above assumptions are generally satisfied by most enzyme kinetic functions, as catalytic rates are typically a monotonic function of the substrate concentration. In particular, the assumptions in (10.17) are met by common kinetics such as mass-action, Michaelis-Menten and Hill equations. The optimization problem considered in [304] corresponds to a free final-time optimal control problem:

$$ e^*(t) = \arg \min_{e \in U} \int_0^T (1 + \alpha' e(t)) dt, \qquad (10.18) $$

where $e(t)$ is the vector of enzyme concentration, $\alpha$ is an $(n+1)$-dimensional vector of tuneable weights, $T$ is a free optimization horizon, and $U$ is a constraint set as in (10.1). The first term in the objective function (10.18) accounts for the total time taken to activate the pathway from the “off” state up to a steady state flux, while the second term weighs the cost of pathway activation. To account for limited availability of cellular resources, the authors also included a temporal constraint on the enzyme concentrations:

$$ \sum_{i=0}^{n} e_i(t) \le e_{\text{tot}}, \qquad (10.19) $$

which is a relaxation of the constraint originally employed by Klipp et al in (10.8), as well as a terminal constraint of the form:

$$ e_i(t) = \frac{F_{\text{goal}}}{g_i(x_i(T))}, \quad \text{for } t \ge T, \qquad (10.20) $$

where $F_{\text{goal}}$ is a (constant) target pathway flux, similar as in (10.13). The terminal constraint ensures that the pathway reaches a steady state at the final time $T$. Using Pontryagin's Minimum Principle [294], the authors showed that the optimal enzyme concentrations follow a bang-bang temporal profile that matches the order in which they act on the pathway. This result was shown to be independent of the weight $\alpha$, the number of enzymatic steps, and valid for a wide range of enzyme kinetics satisfying the assumptions in (10.17), thus extending the original finding in [300] to a larger class of pathways. Figure 10.1 shows a numerical example of the optimal activation pattern obtained for an unbranched metabolic pathway of length three.