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)$ (Fig. 8.4D). 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 (Fig. 8.4C). 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 Growth laws derived 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.

### 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 (Fig. 8.3B). This model consists of two reactions and two proteome sectors. According to the proteome allocation constraint, the sum of the
A

$$x_0 \rightarrow x_1 \rightarrow x_2 \rightarrow x_3$$

B

<figure><img src="image_66.png" /><figcaption>Figure 10.1: Example of optimal enzyme expression in an unbranched metabolic pathway. A. A simple scheme of the metabolic pathway. B. Time evolution of the optimal enzyme expression $u_i$ and metabolite concentration $x_i$. For the simulations, the functions $g_i$ are Michaelis-Menten with constants $k = (1, 2, 4, 3)$ s$^{-1}$, $K = 1$ mM, $V = 0.2$ mM s$^{-1}$ and $x_0 = 5$ mM. Enzymatic weights are set to $\alpha_i = 1$ mM$^{-1}$ s and maximum enzyme availability $E_{tot} = 1$ mM. Resulting activation times are $t_0 = 1.59$ s, $t_1 = 2.2$ s and $t_f = 2.55$ s.</figcaption></figure>

In this section we have reviewed some optimal control approaches for the optimization of unbranched metabolic pathways. While differing in their formulations and solution strategies, these approaches provide substantial computational evidence that some temporal patterns observed in metabolic dynamics can be understood as the solution of an optimal control problem. In the next section we focus on approaches that go beyond individual pathways and include additional components and processes of the cellular machinery.

## 10.4. Dynamic optimization of resource allocation in coarse-grained models of cellular growth

In the previous section, we considered models that were essentially limited to metabolic pathways. The optimization problems were formulated in terms of the allocation of enzymes to the different reactions in the pathway. In this section, we generalize the perspective by increasing the scope of the models from metabolism to protein synthesis and growth. The optimization problems concern the allocation of resources to the synthesis of enzymes catalyzing different metabolic reactions, but also to the synthesis of ribosomes in charge of the production of proteins. Growth is explicitly defined in terms of the increase of protein mass, and leads to growth dilution of all cellular components. The models are very similar to those considered in Chapter 8, but the optimization problems are dynamic rather than static. That is, instead of searching an allocation of cellular resources to the synthesis of different classes of proteins that is optimal at steady state, during balanced growth, we are interested in finding a time-varying resource allocation strategy optimizing an objective defined over an interval of time, e.g., during a transition between two states of balanced growth.

We consider the class of models with dynamics given by Eq. 10.5, where the input $u$ is interpreted as the (time-varying) resource allocation strategy. Among the cellular components $x$, we distinguish between metabolites and proteins, with concentrations $c$ and $p$, respectively. Accordingly, the concentration vector can be written as $x = [c, p]'$. We also distinguish between enzymatic reactions and protein synthesis reactions. While the former have metabolites as substrates and products, the latter convert metabolites (amino acids) into proteins. An enzymatic reaction $i$ has the following
<figure><img src="image_94.png" /><figcaption>(A) Conservation of mass and steady state</figcaption></figure>

$\frac{d \text{ cell component}}{dt} = \text{incoming flux} - \text{outgoing flux} = 0$

<figure><img src="image_95.png" /><figcaption>(C) Volume and surface area</figcaption></figure>

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

$P_{\text{GlucoseUptake}}V_{\text{GlucoseUptake}} + \dots + s_0 = \text{cell surface area}$

<figure><img src="image_97.png" /><figcaption>(D) Reaction fluxes</figcaption></figure>

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 \rightarrow y$ is carried out by proteome sector $p_{x \rightarrow y}$. The maximal rate is reached for saturating substrate concentrations and is determined by the size of the proteome sector.

proteome sectors must sum to one (according to Eq. (8.7)):

$$p_{n \rightarrow x} + p_{x \rightarrow 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 \rightarrow x} = p_{n \rightarrow x} \beta_{n \rightarrow x}; \quad J_{x \rightarrow B} = p_{x \rightarrow B} \beta_{x \rightarrow 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 \rightarrow x} = J_{x \rightarrow B}. \qquad (8.14)$$
reaction rate function:

$$v_i(t) = k_i p_j(t) h_i(c), \qquad (10.21)$$

where $k_i$ is a catalytic constant, $p_j$ the concentration of protein j, and $h_i$ a function describing enzyme saturation. Enzyme saturation is determined by the substrates, products, and activators/inhibitors of the reaction. Typical rate functions $v_i$ follow mass-action kinetics or (ir)reversible Michaelis-Menten kinetics. The synthesis of protein i is associated with the reaction-rate function

$$v_i(t) = u_i(t) v_R(t), \qquad (10.22)$$

where $v_R$ is the total protein synthesis rate defined by

$$v_R(t) = k_R p_R(t) h_R(c(t)), \qquad (10.23)$$

with $k_R$ the maximum protein synthesis rate, $p_R$ the concentration of ribosomes, and $h_R$ a function describing the saturation of ribosomes by their substrate, that is, amino acids (or more precisely, tRNAs charged with amino acids). The function $u_i$ in Eq. 10.22 is a time-varying resource allocation function, describing the fraction of the total protein synthesis rate allocated to the synthesis of protein i. The fractions are non-negative and sum to 1, that is, for every time t,

$$\sum_i u_i(t) = 1, \quad \text{and} \quad u_i(t) \ge 0, \quad \text{for all } i. \qquad (10.24)$$

In most models, the biomass of a growing cell population is equated with the mass of proteins, the most abundant cellular component (Chapter 2). Under the further assumption that the biomass density is constant, it follows that the total protein concentration $p_{tot}$ must be constant, where

$$p_{tot} = \sum_i p_i(t), \qquad (10.25)$$

with the index i running over all proteins. Moreover, the growth rate reduces to the relative (or specific) increase of the protein mass, which leads to

$$\mu(t) = \frac{v_R(t)}{p_{tot}} = \frac{k_R p_R(t) h_R(c(t))}{p_{tot}}. \qquad (10.26)$$

The above model couples metabolism, protein synthesis, and growth in a single formalism, in the spirit of the small resource allocation models discussed in Chapter 8.

Figure 10.2 gives an example of a resource allocation model, describing a simple self-replicatory microbial system [305, 306] inspired by the model of Scott et al. [307] (see Chapter ?? related models). The model divides the proteome into three categories: ribosomes, enzymes, and housekeeping proteins, with concentrations $p_Q$, $p_R$, and $p_M$, respectively. In addition to the three categories of protein, we add a metabolite representing the precursors for protein synthesis, with concentration c. The precursors are produced from nutrients in the environment at a rate $v_M$, a macroreaction catalyzed by the enzymes. Protein synthesis occurs at a rate $v_R$, catalyzed by the ribosomes. The resource allocation functions $u_Q$, $u_R$, and $u_M$ determine the fraction of the protein synthesis rate assigned to each of the three protein categories, where $u_Q$ is assumed to be a constant, growth-rate-independent fraction. The rate equations for the metabolic and protein synthesis reactions follow irreversible Michaelis-Menten kinetics, where the substrate concentration in the medium is assumed saturating.

The resource allocation functions in the model are not explicitly specified by regulatory mechanisms, but assumed to
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.* [203] and later confirmed by Bremer *et al.* [204]. 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$. Fig. 8.3C. 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 [196].

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. 8.18 - 8.20 reveals that
<figure><img src="image_67.png" /><figcaption>Figure 10.2: Example of optimal resource allocation strategy in a coarse-grained model of microbial growth. **A.** Representation of simple self-replicator model of microbial growth. **B.** Model and optimization problem for the self-replicator shown in panel A, as discussed in the text. **C.** Optimal dynamic resource allocation strategy, in terms of the fraction of resources attributed to ribosome synthesis ($u_R$). **D.** Time-varying protein mass fractions corresponding to the optimal solution shown in panel C. The parameter values used for the simulation are $k_M = 0.5$, $k_R = 1$, $K_R = 0.5$ and $u_Q = 0.6$.</figcaption></figure>

follow a dynamic pattern optimizing an objective criterion. In many cases, the objective criterion is based on the hypothesis that microorganisms have evolved to maximize growth, that is, the accumulation of biomass. While this hypothesis can be criticized on theoretical and empirical reasons, it is a reasonable choice in well-mixed environments and provides an interesting baseline in other environments. In the model framework considered here, this gives rise to the following objective function:

$$ \max_{u \in U} J(x(t), u(t), 0, T) = \int_{0}^{T} k_{R} p_{R}(t) h_{R}(c(t)) dt, \quad (10.27) $$

where like in the general case of Eq. 10.1, U denotes the set of admissible profiles for the resource allocation functions u. Note that the maximization of growth over the time-interval $[0, T]$ amounts to taking the integral of the instantaneous
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 maximization* – 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.

### Example 3 - Multiple energy generating pathways

In this example, we focus on a classic question in cell physiology known as overflow metabolism [212, 213]. Within the cell, two primary energy-generating pathways exist: the oxygen-requiring respiration pathway and the oxygen-independent fermentation pathway. It is established that, in the presence of oxygen, the respiration pathway fully oxidizes available nutrients, rendering it more nutrient-efficient in contrast to the fermentation pathway [214]. Utilization of the fermentation pathway is marked by the secretion of byproducts, such as acetate in E. coli or ethanol in yeast, making it inherently wasteful. Intriguingly, experimental observations reveal a counterintuitive phenomenon: even under oxygen-rich conditions, cells often opt for the less efficient fermentation pathway. Under growth rates surpassing a critical threshold, the secretion rate of byproducts, indicating an increased reliance on the fermentation pathway, exhibits a linear rise [147, 215, 216]. This counterintuitive preference for fermentation has long presented a profound question in bacterial physiology.

Based on previous studies [147], we present a coarse-grained model to elucidate this observed phenomenon (Fig. 8.3D). The model postulates steady-state growth on a single nutrient source, denoted as n. This nutrient is taken up from the environment, and channeled towards biomass through the proteome sector $p_{n \to B}$. Additionally, it serves as a precursor for energy generation, either through the respiration pathway catalyzed by proteome sector $p_{n \to r}$ or the fermentation pathway catalyzed by proteome sector $p_{n \to f}$. Thus, according to the proteome allocation constraint (Eq. 8.7), we have:

$$ p_{n \to B} + p_{n \to r} + p_{n \to f} = 1. \quad (8.22) $$

Diverging from earlier models presented in this chapter, our model necessitates two precursors for biomass generation: energy and a carbon precursor. Carbon assimilation is coarse-grained into the biomass generation pathway $n \to B$, while energy is generated through the energy-producing pathways of respiration $n \to r$ and fermentation $n \to f$. Consequently, two mass balance equations are requisite – one for carbon flux and another one for energy flux. The carbon mass balance equates the carbon uptake rate coming from nutrient uptake $J_{in}^C$ to the carbon fluxes utilized for cell biosynthesis $J_{n \to B}^C$, fermentation $J_{n \to f}^C$ and respiration $J_{n \to r}^C$:

$$ J_{in}^C = J_{n \to B}^C + J_{n \to f}^C + J_{n \to r}^C. \quad (8.23) $$

Similarly, the energy balance equation asserts that the energy generated by fermentation $J_{n \to f}^E$ and respiration $J_{n \to r}^E$ equals the energy utilized for the biomass synthesis reaction $J_{n \to B}^E$:

$$ J_{n \to B}^E = J_{n \to f}^E + J_{n \to r}^E. \quad (8.24) $$

Consistent with prior examples in this chapter, we maintain a linear correlation between reaction rates and proteome sector fractions (as per Eq. 8.8).
growth rate over that time-interval, defined by Eq. 10.26. This objective does not generally reduce to maximizing the
instantaneous growth rate at every time-point of this interval.

The question can be asked, for the microbial self-replicator in Figure 10.2, how the cells redistribute their resources
over the different protein categories after a change in environment, in particular a shift of the cells from a poor to a
rich carbon source. In the case of *E. coli*, for example, such as shift might involve a change from minimal medium
with acetate to minimal medium with glucose. Given that *E. coli* grows faster on glucose than on acetate, and that a
higher growth rate requires an increased proportion of resources to be allocated to ribosomes according to the growth
law (Chapter 10), one expects $u_R$ to increase after the shift. Since $u_Q$ is assumed constant, and the resource allocation
functions must sum to 1 at every time-point, this overall increase of $u_R$ must be balanced by a decrease of $u_M$. These
expectations concern resource allocation before the shift (balanced growth on acetate) and a long time after the shift
(balanced growth on glucose), but the growth law provides no information on the pattern of adaptation immediately
after the shift.

In order to investigate the optimal adaptation pattern of $u_R$ immediately after the growth transition, we solve the
dynamic optimization problem specified in Figure 10.2. For the simple example considered here, the optimal solution can
be characterized analytically [308, 305, 306]. This is not possible for more complicated examples, however, which require
the optimal solution to be constructed numerically, using one of the tools discussed in Appendix ]. Figure 10.2C-D show
a typical solution for parameter values estimated from experimental data [306]. Starting from a low value of $u_R$ during
balanced growth on acetate, the optimal resource allocation scheme consists of a sequence of switches between $u_R = 1$
(maximal ribosome synthesis) and $u_R = 0$ (no ribosome synthesis), until an intermediate value of $u_R$ for balanced
growth on glucose is attained. The value of $u_R$ during balanced growth on glucose is higher than that for balanced
growth on acetate, as expected from the growth law.

The sequence of on-off switches followed by the intermediate steady-state value is called a bang-bang-singular solution
in optimal control theory [305, 306]. The solution reflects a dynamic trade-off between the two different functions
contributing to growth: metabolism and protein synthesis. When, due to growth dilution, the ribosome concentration
falls to a level that is limiting for maximal protein synthesis, the synthesis of ribosomal proteins is switched on ($u_R = 1$),
leading to an increase of the ribosome concentration. Switching on the synthesis of ribosomal proteins causes the
synthesis of metabolic enzymes to be switched off. When, due to growth dilution, the concentration of metabolic
enzymes next falls to a level that the precursors produced by the latter become limiting, the synthesis of metabolic
enzymes is switched on ($u_R = 0$) to replenish the precursor pool, and so on.

Optimal solutions with a similar bang-bang pattern were already encountered in the previous section. They also occur
in a model with a more detailed description of different precursor (amino acid) synthesis pathways under the objective
under the minimal time of adaptation after a shift from a medium supplemented with amino acid to a medium lacking
amino acids [309]. In a different type of problem, the development of intestinal crypts, the minimal time to mature
crypts was found to depend on the on-off control of the proliferation of stem and non-stem cells [310]. There is no
convincing experimental evidence that the adaptation of ribosomal synthesis after a nutrient upshift from a poor to a
carbon source follows a bang-bang singular pattern. The interpretation of proteomics data after a nutrient upshift in *E.
coli* shows that the simple upregulation of ribosomal resource allocation to the steady-state value for growth on the rich
nutrient captures the ribosomal protein expression data well [311].

This example serves to emphasize that, while the optimality assumption may lead to thought-provoking predictions, these
need to be confronted with experimental data. In case the optimal solutions do not agree with the data, several revisions
of the problem could be considered. While growth optimization was chosen as the objective criterion in the example of
Figure 10.2, there is evidence that during balanced growth, microorganisms find a trade-off between maximizing growth
rate in a given environment and minimizing necessary adjustments to other environments [296]. The problem could
therefore be generalized to a multi-criteria optimization problem. An example is the analysis of a model similar to that
considered here under the objectives of biomass maximization and minimal adaptation time after a nutrient shift [312].
The formulation of the optimization problem in Figure 10.2 does not put any constraints on valid optimal resource
allocation strategies, except that the individual functions $u_i$ components need to sum to 1 (Eq. 10.24). Bearing in
Both fermentation and respiration reactions utilize a carbon substrate and produce energy, with a key distinction lying in their nutrient utilization efficiency. The ratio of carbon utilized in these reactions to energy generated is expressed as:

$$J_{n \to r}^{E} = \epsilon_{n \to r} J_{n \to r}^{C}; \quad J_{n \to f}^{E} = \epsilon_{n \to f} J_{n \to f}^{C} \qquad (8.25)$$

Given that the respiration pathway exhibits higher nutrient efficiency than the fermentation pathway: $\epsilon_{n \to r} > \epsilon_{n \to f}$.

Concluding the model description, we incorporate the cellular requirements for growth precursors (energy and carbon) and the proteome. Under carbon limitation, the proteome fraction dedicated to cell biosynthesis $p_{n \to B}$ exhibits a linear growth rate dependence [147, 217, 205, 218]:

$$p_{n \to B} = p_0 + \sigma_{n \to B} \lambda. \qquad (8.26)$$

The growth rate correlates with the flux of growth precursors, adhering to a fixed stoichiometry of the metabolic network [219, 220]:

$$J_{n \to B}^{E} = \sigma_{E} \lambda; \quad J_{n \to B}^{C} = \sigma_{C} \lambda. \qquad (8.27)$$

Another key assumption of the model posits that, while the respiration pathway is more nutrient-efficient, utilizing less nutrients per energy unit generated, the fermentation pathway is more proteome-efficient, requiring a smaller proteome fraction per energy unit generated. This assumption is embodied in the efficiency parameters of the reaction fluxes: $\beta_{n \to f} > \beta_{n \to r}$.

To validate the efficacy of our model in capturing the experimentally observed linear increase in acetate secretion with high growth rates, we endeavored to predict acetate secretion as a function of growth rate. The acetate secretion rate is governed by the flux through the fermentation pathway, represented by $J_{ac} = S_{ac}J_{n \to f}^{C}$, where $S_{ac}$ is determined by the involved stoichiometry. Solving Eqs 8.22 - 8.27 for acetate secretion yields an expression that increases linearly with the growth rate:

$$J_{ac} = \frac{S_{ac}}{\epsilon_{n \to f}} \beta_{E} (p_{E} - \lambda(\sigma_{x \to B} + \frac{\sigma_{E}}{\beta_{x \to r}})). \qquad (8.28)$$

where $\beta_{E} = \frac{\beta_{n \to r} \beta_{n \to f}}{\beta_{n \to r} - \beta_{n \to f}}$ and $p_{E} = 1 - p_{0}$. The negative value of $\beta_{E}$, arising from the higher proteome efficiency of the fermentation pathway, results in a positive slope and a negative intercept on the $J_{ac}$-axis. The model provides a good quantitative fit to the experimental observation [147]. The critical growth rate $\lambda_{cr}$, signifying the growth rate at which the cell activates the fermentation pathway, occurs when $J_{ac} = 0$, giving $\lambda_{ac} = \frac{p_{E}}{\sigma_{n \to B} + \sigma_{E}/\beta_{n \to r}}$.

It is crucial to highlight the key assumption underlying this solution, which lies in the relative efficiencies of the energy-generating pathways. At high growth rates, the cell encounters inhibition not only in its ability to rapidly extract energy from the nutrient but, more significantly, it is constrained by the available proteome. Consequently, the cell shifts to utilize the more efficient fermentation pathway.

It is also noteworthy to identify the assumptions overlooked by the model. For instance, the model excludes the proteome sector for nutrient uptake, coarsely integrating it into the biomass biosynthesis and energy generation pathways. While this assumption is reasonable for growth on a single nutrient, a model considering multiple nutrients with varying uptake efficiencies necessitates the inclusion of proteome sectors for nutrient uptake. Further analysis of the model can be found in [147, 221].
mind that the regulatory mechanisms underlying a resource allocation come with a cost, and need to respect certain
physical constraints, the predicted resource allocation strategy may not be feasible. When such constraints are taken
into account, the optimal solution may no longer be bang-bang singular, but resemble the observed adaptation pattern
[313, 306].

In summary, the dynamical optimization approach for studying microbial growth presented here provides a way to test the consequences of hypothesized objective functions in combination with simple resource allocation models. The predictions can be confronted with experimental data, but may also inform the redesign of microbial strains for metabolic engineering purposes [314] (Box 10.A).

Box 10.A: Time-dependent optimization for biotechnological applications and process design

As explained in Section 10.2, time-dependent optimization problems are defined by an objective function, expressing the criterion that microorganisms presumably optimize. In the case of microorganisms growing in natural conditions, the choice of a particular objective function is difficult to make and several functions may qualify. For example, microorganisms could be assumed to maximize their biomass over a given interval of time or minimize the time to adapt to their new environment after a change in growth conditions. The choice may be somewhat arbitrary and in many cases it makes sense to consider a multi-criteria optimization problem. Even more fundamentally, the idea that microorganisms have evolved to the point that they optimize one or several objective functions, is controversial.

In bioengineering applications, however, the formulation of an objective function is less problematic. In this context, the objective function is not assumed to have evolved through natural selection, but is rather stipulated by the metabolic engineer in an *a-priori* manner, in agreement with a practical objective. Possible objective functions are the maximal amount of fermentation product that can be obtained from a given amount of substrate (maximal yield) or the minimal time to produce a given amount of fermentation product (maximal productivity). The use of optimal control methods for process design in bioengineering is well-known [315, 316]. Most of these methods, however, treat microbial growth as a black box and do not provide much detail about the underlying cellular processes, contrary to the formalisms discussed in this chapter. Opening up the black box of microbial growth allows the use of control variables that go beyond standard process parameters of the bioreactor and represent directed perturbations of specific cellular processes.

One example is the use of coarse-grained models of microbial growth for the design of optimal operating conditions for the so-called growth switch [314]. The growth switch is a synthetic regulatory circuit allowing growth of *E. coli* to be arrested in order to passively reorient the resources thus becoming available towards the production of a metabolite of interest [317]. The maximal production of this metabolite from a given amount of substrate, within a given interval of time, was formulated as an optimal control problem. Its solution showed that the optimal solution consists of two phases: a first phase of maximal biomass production followed by a second phase of maximal product synthesis [314]. The conclusion that this two-phase procedure is optimal corresponds well with established practice in biotechnology [318]. Very similar conclusions were attained in related work by Jeanne *et al.* [319].

10.5. Dynamic flux balance analysis (dFBA) of metabolic networks

Dynamic Flux Balance Analysis (dFBA) is an extension of Flux Balance Analysis (FBA) as described in Chapter 5, that can simulate the interactions between the metabolism of an organism and its dynamic environment. In contrast to the constant, steady-state flux solutions that are generated by classical FBA, dFBA yields flux solutions that may dynamically depend on concentrations of extracellular metabolites, such as sugars or other carbon sources, dissolved oxygen, or secreted waste metabolites. Applying these fluxes to the concentration balance of extracellular metabolites also permits to capture dynamic changes in these concentrations due to the cells' metabolic activity, and track the resulting overall biomass growth. It is noted that a basic assumption of dFBA is that organisms rapidly reach intracellular steady state
<figure><img src="image_98.png" /><figcaption>Figure 8.5: Schematic of the dynamic growth model – The model focuses on key cellular processes: nutrient uptake, transcription and translation. Enzymes (shown in blue and dark green) import and metabolize extracellular nutrient (shown in orange), which yields energy (yellow). Availability of energy impacts transcription and translation, however, it is assumed that energy consumption is dominated by translation. The different species of mRNA compete for ribosomes (light green), and their translation consumes energy. Assuming that biomass is dominated by protein, the total rate of translation determines the rate of growth (lower right). Four classes of proteins are modelled: ribosomes, nutrient transporters, enzymes and other house-keeping proteins (red).</figcaption></figure>

## 8.4 Mechanistic links between cellular trade-offs, gene expression, and growth

This section presents a coarse-grained cell model that describes the dynamic adaptation of global mechanisms driving the growth of bacterial cells. Compared to the models previously described in this chapter, this model is dynamic, i.e. not based on steady-state assumptions, and it has a higher level of granularity. It is also based on explicit mechanisms, which allows extension with additional mechanisms of interest, for example, the effects of antibiotics or of heterologous gene expression on cellular growth.

Energy metabolism and protein production are the main pillars of biomass production and cell growth, and form the basis of the growth model. A set of ordinary differential equations describes the dynamic interplay of (i) nutrient internalization and catabolism, (ii) transcription, and (iii) and translation (see Fig. 8.5). A key assumption of the model is that biomass is dominated by proteins, and so the cellular growth rate corresponds to the total rate of protein synthesis via translation. All processes are part of a feedback loop in which the final protein products act as catalyzers of the model reactions, creating a self-replicating system.

In its basic form, the growth model includes 14 intracellular variables: internal nutrient, $s_i$; energy molecules, $a$; and four types of proteins along with their corresponding free ($m_x$) and ribosome-bound mRNAs ($c_x$). Of the four types of proteins considered, there are three groups of catalyzing molecules: transporters ($e_t$), metabolic enzymes ($e_m$) and ribosomes ($r$), and one group of housekeeping proteins ($q$). As the model does not assume steady state, the different reactions are defined in terms of reaction rates instead of reaction fluxes. A simplified description of the main reaction rates of the model is shown in Table 8.1. For details on all reactions and parameters, readers are referred to the supplementary information of [222]. In what follows, the focus will be on the conceptual aspects underlying the prediction of cellular growth rate, and some examples of model applications.
in response to extracellular perturbations, and on the long run no metabolite can accumulate or deplete.

In general, a dFBA model comprises three main parts as demonstrated in Figure 10.3: the dynamic equations, in the form of differential equations, for biomass and extracellular metabolites, constraints on the fluxes as in the FBA model, and an optimization objective that determines how to choose the optimal fluxes.

We first consider the dynamic equations used for dFBA.

The biomass dynamics are given by

$$ \dot{X} = \mu X, \qquad (10.28) $$

where $X$ denotes the biomass concentration, typically measured as dry mass in g/L, and $\mu$ denotes the growth rate, typically measured in 1/h. In principle, this equation follows the equations for balanced growth. However, instead of using simple models, like a Monod equation for the growth rate, the growth rate is taken from the value of the biomass reaction in an FBA model (check in Chapter 5!).

Denoting the concentrations of the extracellular metabolites that are modelled dynamically as the vector $c$, the dynamics for these metabolites can be formulated as the differential equation

$$ \dot{c} = S_{exch}vX. \qquad (10.29) $$

Here, $v$ is the flux vector for the complete metabolic network, including uptake and production reactions for exchange metabolites, and $S_{exch}$ is the stoichiometric matrix that links these reaction fluxes to the metabolite concentrations which are balanced dynamically. Multiplication with the biomass $X$ is necessary, since the flux values in the FBA model are determined relative to biomass, whereas the concentrations $c$ of the dynamic metabolites are relative to the system volume. The equations are not yet closed, because the fluxes $v$ (including the growth rate $\mu$ as one element of the flux vector) still need to be determined by optimization.

As constraints, two types of constraints are used in dFBA models. A flux balance constraint as in steady state FBA models is applied to the concentrations of all metabolites that are not dynamically balanced in (10.29), e.g., intracellular metabolites. This steady state constraint is given by

$$ S_{int}v = 0, \qquad (10.30) $$

where $S_{int}$ is the stoichiometric matrix that links reaction fluxes in the vector $v$ to the steady state metabolites. Further, upper and lower bounds need to be put on the individual reaction fluxes. In contrast to classical FBA, where these bounds are constant, in dFBA flux bounds can depend on concentrations of metabolites in the vector $c$. This is mostly applied to uptake reactions for nutrients, and often as Michelis-Menten kinetics. For example, if $c_i$ is the concentration of a sugar substrate, and $v_i$ is the uptake reaction for this substrate (conventionally negative in FBA models), bounds of the form

$$ -\frac{V_{i,max}c_i}{K_M + c_i} \le v_i \le 0 \qquad (10.31) $$

would be used, where $V_{i,max}$ and $K_M$ are the common parameters of the Michaelis-Menten kinetics (Chapter 3 [MET]).

In recent years, dFBA is increasingly applied for the simulation of dynamic biological systems, especially due to the promising use of GSMMs for interpreting cell physiology and evolution, as well as for guiding metabolic engineering and bioprocess design and optimization [320, 321, 322]. The dFBA applications based on GSMMs include the bacteria *Escherichia coli* [323, 324, 325, 326, 327, 328, 329, 330] and *Lactococcus lactis* [331], as well as the yeast species *Saccharomyces cerevisiae* [332, 333, 334, 335, 336, 324, 327, 337, 338, 330] and *Scheffersomyces (Pichia) stipites* [327]. However, the majority of dFBA applications use small-scale metabolic models, most of which include less than 100 reactions. Such applications include models of bacteria, like *Escherichia coli* [339, 340, 341, 327, 342, 343, 344, 345], *Corynebacterium glutamicum* [346, 326, 347, 329] and *Bordetella pertussis* [348], models of yeast for wine fermentation [349] and *Saccharomyces cerevisiae* [350, 338, 351, 352], but also plant and animal models, such as a model for the photosynthetic metabolism of C3 plants [353], a four-tissue (leaf, root, seed, and stem) model of the core metabolism
Building on the assumptions of mass balance and proteome allocation described in Section 8.2 of this chapter, the model centers around three fundamental constraints, namely (i) a finite pool of cellular energy that fuels protein biosynthesis, (ii) a finite pool of ribosomes for which mRNAs compete for translation, and (iii) a finite cell mass. As a result, the model predicts the dynamic allocation of internal resources and its emergent impact on cellular growth rate without the need to assume growth rate maximisation.

<table><caption>Table 8.1: Summary of main model reactions and their accompanying rates. The four proteins represented in the model are denoted in the reactions by $x, x \in r, e_t, e_m, q, \gamma(a)$ is the rate of translational elongation, defined as $\frac{\gamma_{max}a}{K_{\gamma+a}}$, and $n_x$ is the average length of a protein molecule in amino acids. The parameter $n_s$ represents nutrient quality and determines the yield of energy per catabolized nutrient.</caption><thead><tr><th>Description</th><th>Reaction</th><th>Reaction rate</th></tr></thead><tbody><tr><td>Nutrient internalisation</td><td>s &rarr; s<sub>i</sub></td><td>e<sub>t</sub> &middot; <sup>v<sub>t</sub>s</sup>&frasl;<sub>(K<sub>t</sub>+s)</sub></td></tr><tr><td>Nutrient metabolism</td><td>s<sub>i</sub> &rarr; n<sub>s</sub>a</td><td>e<sub>m</sub> &middot; <sup>v<sub>m</sub>s<sub>i</sub></sup>&frasl;<sub>(K<sub>m</sub>+s<sub>i</sub>)</sub></td></tr><tr><td>Transcription</td><td>&Oslash; &rarr; m<sub>x</sub></td><td>&omega;<sub>x</sub> &middot; <sup>a</sup>&frasl;<sub>(&theta;<sub>x</sub>+a)</sub></td></tr><tr><td>Ribosome binding</td><td>m<sub>x</sub> + r &leftrightarrow; c<sub>x</sub></td><td>k<sub>b</sub> &middot; m<sub>x</sub>r, k<sub>u</sub> &middot; c<sub>x</sub></td></tr><tr><td>Translation</td><td>c<sub>x</sub> + n<sub>x</sub>a &rarr; x + m<sub>x</sub> + r</td><td>c<sub>x</sub> &middot; <sup>&gamma;(a)</sup>&frasl;<sub>n<sub>x</sub></sub></td></tr></tbody></table>

### 8.4.1 Model definitions

**Growth rate and biomass synthesis** Based on the assumption that biomass is dominated by protein, and other contributions are negligible, the biomass $B$ of a cell can be calculated by summing over the coarse-grained proteome,

$$B = \sum_{x} n_{x}x + n_{r} \sum_{x} c_{x}, \quad x \in r, e_{t}, e_{m}, q, \qquad (8.29)$$

which sums over all proteins ($x$) and mRNA-bound ribosomes ($c_x$), with $n_r$ and $n_x$ denoting the lengths of proteins in terms of amino acids. Equation (8.29) is equivalent to the mass balance assumption described in section 1.2.1 of this chapter. As a consequence, the proteome allocations, defined by $\phi_x = x/B$ for $x \in \{e_m, e_t, r, q\}$ sum to 1, i.e. $\sum_x \phi_x = 1$.

Similar to the previous examples in this chapter (Section 8.3), the model correlates the growth rate with biomass production, which depends on translating ribosomes and their translation elongation rate $\gamma(a)$. Importantly, the rate of elongation depends on the energy produced in the catabolic processes described in the model, which dynamically couples protein synthesis with metabolism. Defining the number of translating ribosomes $R_t = \sum_x c_x$, the change in cellular biomass over time becomes

$$\frac{dB}{dt} = \gamma(a)R_t - \lambda B. \qquad (8.30)$$

The second term, $\lambda B$, accounts for dilution via redistribution of mass to daughter cells at division. In homeostatic conditions, that is when $B$ is in steady state and so $\frac{dB}{dt} = 0$, it then follows that $\lambda^*$ is proportional to the rate of protein synthesis. To define growth dynamically,

$$\lambda(t) := \frac{\gamma(a)R_t}{B_0}, \quad B_0 > 0. \qquad (8.31)$$

Setting $B_0$ to the typical biomass of a cell in mid-exponential growth ensures that cells will have a steady-state biomass of $B^* = B_0$.

**Rate of translation** In actively growing bacteria, protein synthesis, and in particular translation-associated processes, account for a major part of the energy budget. The model assumes a simplified mechanism to derive the dependence of the translation rates on the energy levels of the cell. It is assumed that each elongation step of translation consumes a fixed amount of energy (Figure 8.6), and further that intermediate reactions are in quasi-steady state. It can then be shown that the net rate of translation elongation
<figure><img src="image_68.png" /><figcaption>Figure 10.3: Schematic representation of dFBA. As in FBA, the intracellular environment in dFBA is represented by a linear programming (LP) optimization problem that describes the metabolism of the microorganism based on its genome-scale metabolic model (GSMM). FBA assumes that all intracellular metabolite concentrations remain constant while the cells optimally distribute their metabolic fluxes to maximize their growth rate and hence, an LP can calculate the growth rate, as well as the intracellular and exchange fluxes of the GSMM. The calculated growth rate and exchange fluxes can be used to update the extracellular environment. The extracellular environment in dFBA is represented by ordinary differential equations (ODEs) that describe the mass balance equations for biomass and metabolites found outside of the cell. Moreover, the intracellular GSMM and the extracellular mass balance equations can be linked through kinetic rules for substrate uptake, like the Michaelis-Menten equations, that can raise concentration-dependent constraints for exchange fluxes and predict growth rate dependencies on substrate concentrations.</figcaption></figure>

of *Arabidopsis thaliana* [354], a model for fatty acid metabolism and lipid accumulation in rat hepatocytes [355], and a model for energy metabolism in myocardial cells [356].

It is noted that most of the dFBA applications for microorganisms simulate microbial fermentations under batch or fed-batch conditions. Since dFBA can be used for the analysis, control and optimization of biochemical processes, many dFBA applications focus on either dynamic metabolic engineering or optimal control of bioreactors, or both simultaneously. Dynamic metabolic engineering studies can predict the effect of strain gene insertion and deletion on the dynamic behavior and productivity of a bioprocess [332, 323, 333, 334], while optimal control of batch or fed-batch operation of bioreactors is important for the production of desired chemicals [332, 352]. Finally, dFBA has also been expanded for the study of microbial communities, where each microorganism is represented by an LP that is solved independently [357, 358, 359]. Co-culture simulations with dFBA can predict possible consortia compositions, as well as metabolic engineering approaches to improve the productivity of the consortia, but they are out of the scope of this chapter.

Coming to the mathematical formulation of dFBA models, dFBA is an optimization problem coupled with a system of ordinary differential equations, that can be solved with the help of various mathematical and numerical techniques. Even though dFBA was first introduced in 1994 [339], it was not formalized until 2002 [340]. The existing formalized solution approaches that are going to be discussed here involve the static optimization approach (SOA), the dynamic optimization approach (DOA), and the direct approach (DA). More recently, reformulation approaches and surrogate models for the optimization problem have also been proposed in order to ease the computational complexity of dFBA simulations. This complexity arises from several characteristics of dFBA. More specifically, the solution of dFBA problems faces challenges in terms of:

1. problem size and scalability: As the size of the metabolic network increases, the computational cost increases. For
<figure><img src="image_99.png" /><figcaption>Figure 8.6: Mechanistic derivation of the translational elongation rate – The model assumes that each elongation step consumes a fixed amount of energy. In a first step, energy reversibly binds the mRNA-ribosome complex, upon which elongation takes place. Once the peptide reaches its final length, the protein is released and ribosome and mRNA are freed up.</figcaption></figure>

takes the form

$$ \gamma(a) = \frac{\gamma_{\max} a}{K_{\gamma} + a}. \qquad (8.32) $$

Here, $\gamma_{\max}$ denotes the maximal rate of translation elongation per ribosome and $K_{\gamma}$ the energy threshold of half-maximal elongation. For any protein $x$, the rate of its translation is then given by

$$ \nu_x(c_x, a) = c_x \frac{\gamma(a)}{n_x}, \qquad (8.33) $$

where $c_x$ denotes ribosomes bound to mRNA of type $x$ and division by $n_x$ accounts for the number of elongation steps to take place for the production of one $p_x$.

**Rate of transcription** The model assumes that transcription is energy-dependent, but that its consumption is negligible compared to that of translation. Analogous to translation, under the assumption of fixed energy consumption per elongation step, the rate of transcription takes the same shape and is defined by

$$ \omega_x(a) = \frac{\omega_x a}{\theta_x + a}, \quad x \in r, e_t, e_m. \qquad (8.34) $$

Here, the energy threshold of half-maximal transcription, $\theta_x$, is specific for each proteome sector $x$, which dynamically links the proteome allocations $\phi_x$ with different growth conditions. In particular, $\theta_r \gg \theta_x$ for $x \neq r$ ensures that the ribosomal sector increases in rich growth conditions (cf. growth laws in Fig. 8.1C).

In addition, the model assumes that the transcription of household genes is negatively auto-regulated to maintain near constant levels across different conditions. Therefore

$$ \omega_q(q, a) = \frac{w_q a}{\theta_q + a} \cdot \mathcal{I}(q), \quad \text{with} \quad \mathcal{I}(q) := \frac{1}{1 + (q/K_q)^{h_q}}, \qquad (8.35) $$

where $\mathcal{I}$ is the auto-inhibition function with threshold $K_q$ and Hill-coefficient $h_q$.

### 8.4.2 Model predictions

The model recovers the bacterial growth laws through the automodulation of finite cellular resources in response to changing environments. It robustly fits empirical data (Fig. 8.7), suggesting the growth laws are an emerging property of the constraints integrated into the modeling approach.

The model predicts a hyperbolic dependence of the growth rate on nutrient availability as described by Monod's law (Fig. 8.7 inset), derived using the conservation of mass assumption and when $\phi_r \ll \phi_q$. Energy is created from the metabolism of internalized nutrients and determines the rates of transcription ($\omega_x(a)$) and translation ($\gamma(a)$). In the absence of antibiotics, the latter is proportional to the growth rate of the cell as described in Eq. (8.31). As the nutrient quality is increased, more energy will be available and therefore more transcription will occur. Due to the relationship between transcription thresholds ($\theta_r \gg \theta_x$), the transcription of ribosomes is increased comparatively more, leading to an increase in the ribosomal mass fraction as seen in Fig. 8.7.

In a fixed nutrient condition, inhibiting translation by the addition of an antibiotic increases intracellular energy levels as fewer ribosomes can translate. Again, with $\theta_r \gg \theta_x$, this energy increase leads to a
this reason, simulations that involve large genome-scale metabolic models or multispecies microbial communities
are limited.

2. stiffness: The stiff behavior of dFBA has been observed in many cases, such as the simulation of the diauxic growth in *E. coli* [340].

3. nonlinearity: The presence of nonlinear constraints or objective functions can significantly increase the computa-tional cost.

4. feasibility: The intracellular optimization problem can become infeasible and lead to failure of the integration of the extracellular ODEs.

5. differentiability: The optimal value of the intracellular optimization problem may not be continuously differentiable, which poses an obstacle when dFBA is used for optimal control or parameter estimation.

6. non-unique solutions: The solution of the intracellular optimization problem is usually not unique which can make fluxes unrealistically “jump” between different optimal solutions.

**Static Optimization Approach (SOA)** divides the total time horizon of the dFBA simulation into several smaller time intervals. The optimization problem is solved to obtain the flux distribution at the beginning of each time interval, and then the ODEs are integrated over the time interval with this fixed flux distribution. The dynamics calculated from this time step are used to constrain the optimization problem solved at the beginning of the next time interval, and the process is repeated until the end of the simulation time is reached. SOA can be implemented easily with the use of an Euler scheme for integrating the system and a suitable existing LP solver for solving the FBA at each time step. SOA is also implemented in the constraint-based reconstruction and analysis (COBRA) toolbox for MATLAB [360] which can perform dFBA simulations. Since its implementation is relatively simple, SOA has been widely used in studies for the diauxic [340, 342, 361], aerobic and fermentative [339, 362, 363, 323] growth of *E. coli*, for *S. cerevisiae* fermentations [350, 334, 335, 336], as well as for the growth of other bacterial [331] and plant organisms [354]. Many of these applications include larger-scale or genome-scale metabolic networks, due to the scalability of SOA. However, the main drawback is that SOA is inefficient and can become computationally expensive because it has to solve the optimization problem at each time step. This can be challenging for most dFBA problems which are stiff and require small time steps to ensure accuracy, convergence, and stability of the solution.

**Dynamic Optimization Approach (DOA)** follows closely the general dynamic optimization framework described in Section 10.2: an objective function that depends on the dynamic states of the system over the complete time horizon of interest is formulated, and the dynamics (10.28)–(10.29) and algebraic constraints (10.30)–(10.31) are added as optimization constraints. In other words, DOA discretizes the total time horizon of the dFBA simulation, and then transforms the dynamic optimization problem into a non-linear programming (NLP) problem, which is solved once by simultaneously optimizing over the entire time of the simulation. In this way, DOA obtains the time profiles of fluxes and metabolite concentrations in the system, and allows the formulation of a dynamic objective function, which could provide useful information about the design of genetically modified metabolic networks or the maximization of bioprocess productivity. Because of this characteristic, DOA is often used in dynamic metabolic engineering, parameter estimation and optimal control applications. For example, DOA has been used for simulating the diauxic growth of *E. coli* [340, 364, 365], as well as the growth of engineered *E. coli* strains on glucose [366, 367], and the growth of various eukaryotic organisms, such as *S. cerevisiae* [351, 368], plant [353] and animal [355, 356] cells under genetic and environmental perturbations. On the downside, even though the optimization problem does not need to be repeatedly solved like in SOA, the single NLP of DOA can become easily intractable, as its dimension increases with the fineness of time discretization. Additionally, DOA has been mainly limited to small-scale metabolic networks, since it cannot be easily applied to genome-scale metabolic networks due to the large number of variables and constraints that are introduced in the NLP as the size of the network increases.

**Direct Approach (DA)** has been formulated more recently than SOA and DOA, and directly includes the LP solver for the FBA in the right-hand side evaluator function of the ODEs. In this way, it can take advantage of existing ODE
<figure><img src="image_100.png" /></figure>

Figure 8.7: Mechanistic cell model – Experimental data (coloured circles) and model simulations (lines) depicting the relationship between growth rate and cellular composition. The data describes the ribosomal fraction of the proteome $\phi_r$ in different growth conditions. Each colour represents a different media composition, with increasing drug-free growth going from red to green. The numbers within the circles indicate the addition of the antibiotic chloramphenicol to the growth media at a certain concentration [in $\mu M$]. Although this antibiotic inhibits translation, an increase in $\phi_r$ can be observed through all media compositions. The model fit to the experimental data demonstrates the capacity of this model to describe two of the growth laws. (Inset) Model simulation. Besides the composition, varying the amount of external nutrient in the growth media increases the steady-state growth rate up to a saturation point. This reproduces Monod's growth law.

proportionally larger increase in transcription of ribosomal mRNAs and so to a larger $\phi_R$. In contrast to the scenario without antibiotics, fewer ribosomes can actively translate and therefore the growth rate will be lower. Consequently, a negative dependence of $\phi_R$ and growth rate arises.

### 8.4.3 Applications of the model

Due the coarse-grained modeling of mechanisms and the use of non-steady state dynamics, the model lends itself to modular extension for a range of applications. For example, to reproduce the negative correlation between growth rate and ribosome content amid translational inhibition (Fig. 8.7), the model was extended to account for inhibitory actions of the antibiotic chloramphenicol on ribosomes. Similarly, mechanisms that account for drugs with other modes of action could also be included. Further, in [222], it was shown that the model can be extended to study a number of applications:

Firstly, the model was extended to account for expression of a heterologous gene circuit and predict constraints between heterologous circuit expression, circuit function, and the growth of the host. This has applications in areas such as chemical production in biotechnology, where host-circuit interactions are not understood and where synthetic circuits have to operate robustly in different growth conditions. In this context, the model can serve to quantify host-circuit interactions for a more host-aware design of synthetic gene circuits.

In another application, the model's ability to dynamically predict growth rate emergently from intracellular mechanisms was used as a proxy for evolutionary 'fitness' to study when gene regulation was evolutionarily stable. This was done by augmenting the cell model with population growth, assuming that all cells of a population are identical, and modeling competitive interactions between a resident and mutant strain.

Finally, in [222] it was shown how to use the model to study specific mechanisms within a wider cellular context. With the example of gene-dosage compensation, where the effects of a gene deletion can be reduced by increasing the expression of a paralogous gene, it was shown how and when global regulatory mechanisms caused compensation. The example showed that the constraints underpinning the growth
integrators with adaptive step size and error control that can reduce the number of integration steps and provide better
solution accuracy compared to SOA. DA has been implemented in the ORCA toolbox [369], which complements the
constraint-based reconstruction and analysis (COBRA) toolbox for MATLAB [360]. Furthermore, DA has been used for
studying the diauxic [370], aerobic and anaerobic[341, 371] growth of wild type and engineered *E. coli* strains, the aerobic
growth of *Corynebacterium glutamicum* on glucose and xylose in biorefinery simulations [347], as well as the aerobic and
anaerobic growth of wild type and engineered *S. cerevisiae* strains [332, 333, 335]. Some of these applications involve
dynamic metabolic engineering for product maximization, and many of them include genome-scale metabolic networks,
since DA is relatively easily scalable like SOA.

However, DA requires the LP to be resolved at least once, every time the right-hand side of the ODEs is evaluated [372]. This can make DA computationally demanding, especially for larger metabolic networks. Another major challenge is that when evaluating the right-hand side of the ODEs close to the boundary of feasibility, the LP can become infeasible and make the dFBA simulation fail. The LP can become infeasible either because it is really infeasible and the simulation should be terminated, or because the ODE integrator becomes unable to evaluate the right-hand side of the ODEs and the simulation is discontinued, or erroneous death phase messages are being displayed. The latter can happen as dFBA simulations involve discrete events that correspond to switches in the active set of the LP solution. More specifically, different bases for the optimal solution of the LP can emerge at each time step. Moreover, at the points of change of the active set, the dFBA model is not differentiable, since the optimal value of the LP as a function of the right-hand side of the constraints is not continuously differentiable. This is a problem because the first and second derivatives of the model must be computed when dFBA is used for optimal control or parameter estimation applications. Finally, another drawback emerges due to the primal multiplicity of the LP. As it is well-known, FBA is formulated as an underdetermined problem and therefore, the LP does not have a unique solution [373]. Non-unique optimal reaction fluxes can lead different ODE integrators to different results.

Mathematical details 10.B: dFBAlab

In order to address some of the computational challenges of dFBA, Höffner, Harwood, and Barton proposed a
simulator for dFBA, which was initially coded in FORTRAN [374], but gained popularity when implemented in
MATLAB with the name Dynamic Flux Balance Analysis laboratory (DFBAlab) [359], and more recently in Python
[375]. It is noted that the DFBAlab is compatible with the COBRA toolbox [360]. Based on this dFBA simulator,
it is not necessary to resolve the LP each time the right-hand side of the ODEs is evaluated and consequently,
the solution process becomes faster. This is possible because the FBA solution at an initial time can be used
to compute future optimal solutions by detecting changes in the active set [329] or by computing the optimal
basis [376] of the FBA. Unfortunately, such formulations need to continuously monitor the active set of the LP,
which increases with the size of the metabolic network, or need to choose a basis for the optimal solution that is
most likely to remain optimal as the simulation proceeds [376]. The latter is challenging since the optimal basis
can be non-unique even for a unique optimal solution. Nevertheless, DFBAlab manages to reduce the number
of times that the LP is resolved, and also avoids obtaining infeasible LPs and numerical failure by using the LP
feasibility problem and the Karush-Kuhn-Tucker (KKT) optimality conditions of the FBA problem (see below)
[324, 326, 329]. In addition, the differentiability problem could be solved with the help of non-smooth analysis
which provides optimality conditions in terms of sub-gradients or generalized gradients, for convex and non-convex
functions respectively [325]. Furthermore, to tackle the issue of primal multiplicity of the FBA problem, DFBAlab
performs lexicographic optimization [324, 325, 326, 329].

**10.6. Concluding remarks**

One approach for understanding the response of microorganisms to changes in their environment is to assume that
this response has been optimized by evolution. That is, the regulatory mechanisms controlling the response optimize
laws can also cause global negative feedbacks on proteins affecting growth.

## 8.5 Concluding remarks

In this chapter, we delved into the intricate world of coarse-grained modeling of microbial growth. We began by describing key experimental evidence that has led to what is known as bacterial growth laws. These laws are derived from growth measurements and are deemed to be conserved for various organisms. We then mathematically described the fundamental assumptions necessary to model bacterial growth. Using basic modeling systems, we showed how to analyze such a system and derive fundamental conclusions for bacterial growth. These models reproduce the bacterial growth laws, providing a link between theoretical models and experimental results. Finally, we introduced a more complex model that includes various cell processes such as translation, transcription, and the cellular growth process. Overall, this chapter highlights the power of coarse-grained modeling in unraveling the complexities of microbial growth and offers a framework for exploring a wide range of biological questions.

While this chapter lays a foundation for research on various topics in biology, many areas remain to be explored. For example, the effects of changing environmental conditions such as dynamic changes in nutrient availability, acidity, or temperature are not discussed. Furthermore, various cellular processes such as protein degradation and membrane assembly are not covered in the chapter. Including these processes in a coarse-grained model could potentially lead to the discovery of other growth laws.

In the next chapter, you will explore models that further refine the biological cell and bridge between coarse-grained models and genome-scale models. These models incorporate several of the assumptions discussed here but utilize more knowledge of the metabolic network.

## Recommended readings

### Growth laws in microbiology:

* Monod, The growth of bacterial cultures, *Annual Review of Microbiology*, 1949 [200]. Classical reference for the quantitative modeling of microbial growth.

* Schaechter, Maaløe, and Kjeldgaard, Dependency on medium and temperature of cell size and chemical composition during balanced growth of *Salmonella typhimurium*, *Microbiology*, 1958 [203]. Classical article that introduced the growth law for ribosomes.

* Scott, Gunderson, Mateescu, Zhang, and Hwa. Interdependence of cell growth and gene expression: origins and consequences, *Science*, 2010 [205]. Article that renewed interest in growth laws for the quantitative study of microbial physiology.

* Jun, Si, Pugatch, and Scott, Fundamental principles in bacterial physiology – history, recent progress, and the future with focus on cell size control: a review, *Reports on Progress in Physics*, 2018 [198]. A very complete review of growth laws in microbiology.

### Coarse-grained modeling of microbial growth:

* Hinshelwood, On the chemical kinetics of autosynthetic systems, *Journal of the Chemical Society*, 1952 [223]. Historical reference for coarse-grained modeling of microbial growth.

* Kafri, Metzl-Raz, Jonas and Barkai, Rethinking cell growth models, *FEMS Yeast Research*, 2016 [224]. Review of coarse-grained models of microbial growth.

* de Jong *et al.*, Mathematical modeling of microbes: metabolism, gene expression and growth, *J. R. Soc. Interface*, 2017 [225]. Review comparing coarse-grained models of microbial growth with other modeling frameworks.

* Bruggeman, Planqué, Molenaar, and Teusink, Searching for principles of microbial physiology, *FEMS Microbiology Reviews*, 2020 [226]. Review summarizing biological insights obtained from coarse-grained models.

### Examples of coarse-grained models:

* Molenaar, van Berlo, de Ridder and Teusink, Shifts in growth strategies reflect tradeoffs in cellular eco-
Mathematical details 10.C: Lexicographic Optimization

Lexicographic or hierarchical optimization involves the solution of a series of LPs with auxiliary objectives ranked
in priority order. The use of auxiliary objectives reduces the feasible space and leads to a unique optimal solution,
while the auxiliary objectives can have specific biological meanings and can be selected based on prior knowledge
about the organism [377]. Some of the most used auxiliary objectives are based on the assumption that evolution
leads to the exclusion of inefficient pathways so that cells can biosynthesize the smaller possible number of enzymes.
Examples of such auxiliary objectives include the minimization of enzyme cost [378], the minimization of the total
reaction flux [379], and the minimization of the number of active reactions [380]. However, it has been shown that
such objectives may not be suitable for some engineered cells [348]. In general, it is not trivial to find a series of
auxiliary objectives that are consistent with experimental data, assure uniqueness and preserve continuity of the
optimal solution. In some cases, even when all auxiliary objectives have been used, hierarchical optimization cannot
ensure the uniqueness of the dFBA solution. Apart from the use of auxiliary objectives, auxiliary rules or auxiliary
parameters have also been proposed to address the primal multiplicity of FBA. For example, geometric methods
have been proposed to identify a unique distribution of reaction fluxes for FBA [381], even though there is no
biological evidence to justify such methods.

an objective, or a trade-off between competing objectives, subject to a variety of physical and biochemical constraints.
This approach gives rise to dynamic optimization problems (10.1)-(10.4) that can be solved by techniques from optimal
control theory. Three examples of such problems were considered in this chapter: dynamic optimization of enzyme
expression in metabolic pathways, dynamic optimization of coarse-grained models of cellular growth, and dynamic flux
balance analysis. This does not exhaust the range of possible problems that can be considered. One example is the
combination of the resource allocation perspective with dynamic flux balance analysis [364].

Some of the predictions obtained by means of dynamic optimization seem to be supported by available experimental data,
such as the time-ordering of enzyme expression in a linear pathway. Other predictions cannot currently be tested or may
not be consistent with the available experimental data, such as the dynamical ribosomal protein synthesis pattern. The
contradiction between a predicted optimal response and the observations is interesting, because it indicates that some of
the assumptions underlying the problem need to be revised. The model may not account for all relevant processes taking
place in the cell, important constraints may have been ignored, or the objective may not capture the actual processes
taking place.

The solution of a dynamic optimization problem may be different from the concatenation of the solutions of repeated
static optimization problems defined over short, consecutive time intervals making up the time horizon. For instance, the
optimal pattern of resource allocation over a time horizon may involve the accumulation of a reserve of unused resources
that, while being wasteful in the short run, is beneficial when the whole time interval of interest is considered. One
example was given in the introduction of this chapter, concerning the expression of maltose enzymes in the presence of
lactose [292]. Another example is the accumulation of glycogen in cyanobacteria during daylight, providing the energetic
resources for maintenance metabolism in the night time [289].

The analysis of the growth of microorganisms using dynamic optimization critically depends on the choice of an ap-
propriate objective function. In the examples discussed above, optimization was performed with respect to a single
objective, e.g., the maximal accumulation of biomass over a time interval or the minimal time to deplete a given amount
of substrate. It is plausible, however, that microorganisms have evolved under the necessity to satisfy several objectives
simultaneously. This can be taken into account by formulating a weighted sum of the different objective functions, such
as the simultaneous minimization of throughput time and investment in enzymes in Eq. 10.18. Another approach is to
generalize the optimization problem to a multi-objective optimization problem, with sets of Pareto optimal solutions,
each providing a trade-off between mutually conflicting objectives. One example of such a multi-objective optimization
problem is given by a generalization of the dynamic optimization of enzyme expression in metabolic pathways in Sec-
tion 10.3, with the double objective of minimizing the time to consume a given amount of substrate and minimizing the
nomics, *Molecular Systems Biology*, 2009 [227]. Influential article illustrating the explanatory capacity of coarse-grained models.

○ Weiße *et al.*, Mechanistic links between cellular trade-offs, gene expression, and growth, *Proceedings of the National Academy of Sciences of the USA*, 2015 [222]. Article describing how growth laws for ribosomes can be recovered from coarse-grained model of microbial growth.

○ Basan *et al.*, Overflow metabolism in *Escherichia coli* results from efficient proteome allocation, *Nature*, 2015 [147]. Article describing how proteome allocation constraints can account for overflow metabolism in bacteria.

○ Zavřel *et al.*, Quantitative insights into the cyanobacterial cell economy, *eLife*, 2019 [228]. Example of the use of coarse-grained models for explaining physiological principles underlying growth of less-studied (photosynthetic) microorganisms.

## Problems

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

### Problem 8.1 Linear chain model

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 Linear chain model with Michaelis-Menten rate laws

Solve example 1 when the nutrients are not available in excess. Use Michaelis-Menten relations for both reactions. First, derive the concentration of metabolite $x$ as function of catabolic sector proteome size. What is the minimal size for the catabolic sector? What happens if the catabolic sector is smaller than that? Next, determine the proteome allocation that maximizes the growth rate.

### Problem 8.3 Linear chain model with Michaelis-Menten rate law for the catabolic reaction

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

### Problem 8.4 Simple model with allosteric regulation of catabolic reaction [229]

A metabolic system is growing in an environment with one nutrient available. The system allosterically regulates its catabolic reaction according to the concentration of metabolite $x$. Assume Michaelis-Menten kinetics for all reactions. What is the growth rate as function of catabolic sector proteome size? This is a complex solution, don't solve it analytically and plot a numerical solutions instead. What is the catabolic sector proteome size that maximizes the growth rate?

### Problem 8.5 Growth on a single nutrient that is degraded to both energy and biomass precursors

Consider the model from section 1.3, example 3. Solve the model for the nutrient uptake rate as function of growth rate for:

(a) Growth rates above the onset of acetate secretion

(b) Growth rates below the onset of acetate secretion

### Problem 8.6 Simulating models numerically

Simple coarse-grained models can generally be solved analytically. However, for models with a higher level of granularity, like the one presented in this section, reaching an analytical solution to the model equations is highly complex. Computational approaches that allow numerically solving high-dimensional systems are of great value.

(a) With the help of the provided code and following the detailed description of the ODE system in the SI of [222], implement and solve the system of ODEs. Using this implementation, reproduce Monod's law, as seen in the inset of Figure 8.1.

(b) The nutrient composition of the growth media is the main driver of increasing growth rates. Sim-
<figure><img src="image_69.png" /></figure>

B

Exchange constraints

$$
\begin{align*}
-10.5 \frac{\text{mmol}}{\text{gDW h}} \frac{c_G}{2.7 \frac{\text{mg}}{\text{L}} + c_G} &\le v_G \le 0 \\
-30 \frac{\text{mmol}}{\text{gDW h}} \frac{c_{O_2}}{20 \frac{\text{mg}}{\text{L}} + c_{O_2}} &\le v_{O_2} \le 0 \\
0 &\le v_E \\
-30 \frac{\text{mmol}}{\text{gDW h}} \frac{c_A}{100 \frac{\text{mg}}{\text{L}} + c_A} &\le v_A
\end{align*}
$$

C

Dynamic equations

$$
\begin{align*}
\dot{X} &= \mu X \\
\dot{c}_{O_2} &= v_{O_2} m_{O_2} X + k_L a \left(10 \frac{\text{mg}}{\text{L}} - c_{O_2}\right) \\
\dot{c}_G &= v_G m_G X \\
\dot{c}_E &= v_E m_E X \\
\dot{c}_A &= v_A m_A X
\end{align*}
$$

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

E

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

F

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

Figure 10.4: Example result of a DFBA simulation for the E. coli core model [382] performed with DFBAlab. **A:** Illustration of dynamic metabolites and relevant exchange fluxes. Extracellular metabolites use mass concentration (g/L), exchange fluxes are in molar amount per dry biomass and time (mmol/gDW h). **B:** Concentration-dependent constraints applied to the exchange fluxes during the DFBA simulation. **C:** Differential equation model for biomass concentration X and metabolite mass concentrations. Growth rate μ and exchange fluxes v* (in red) are optimal values from the underlying FBA model. mO2, mG, mE, mA are molar masses of O2, glucose, ethanol, and acetate, respectively. kLa = 8.5 h-1 is the volumetric mass transfer coefficient for oxygen. **D:** Simulation results for concentrations of biomass and metabolites. We can observe four growth phases: aerobic growth on glucose with production of acetate, anaerobic growth on glucose with production of ethanol and acetate, aerobic growth on acetate, and a stationary phase. **E:** Simulation result for oxygen concentration in the liquid medium. Oxygen is depleted in the second growth phase due to mass transfer limitations, but replenishes at the start of the third phase. **F:** Penalty function time course. Increases of the penalty function indicate periods where the underlying FBA model is infeasible. Here, this occurs when the ATP maintenance constraint cannot be satisfied due to a lack of substrates, and happens in this simulation during a brief period where the switch from glucose to acetate as a substrate takes place, since oxygen needs to be replenished first, and in the stationary phase.

concentration of (possibly toxic) intermediate metabolites [383].

Instead of making *a-priori* assumptions about the objectives presumably optimized by microorganisms, one could try to infer the latter from the experimental data. This inverse optimization approach leverages the large amounts of time-course data on the dynamic response of microorganisms to environmental perturbations that have accumulated in the past decade. Inverse optimization requires the solution of complex inverse optimal control problems that have been little explored until now [295].
ulate the model to steady state for different values of nutrient qualities. What model species are
most impacted by an increase in nutrient quality?

(c) As seen in Figure 8.1, the addition of a drug that inhibits protein synthesis results in an upregulation of the ribosomal fraction $\phi_R$. Reproduce Figure 8.1. How do the observed results relate to your answer in question 2?
Chapter 11

Metabolic diversity in cell populations:
probability densities over the flux polytope

Andrea De Martino¹,² and Marcelo Rivas-Astroza³

Chapter highlights

Even in clonal populations, cells appear to be strongly heterogeneous in terms of, e.g., protein levels, RNA levels, sizes at birth or division, interdivision times and elongation rates. Part of this variability is likely due to the inherent stochasticity of gene expression at the level of single cells. It is however known that heterogeneous populations may possess an evolutionary advantage, for instance in variable environments or under stress. Despite appearing to be at odds with the idea of optimality presented in the previous Chapters, metabolic diversity can be described and modeled within the constraint-based framework introduced in the previous Chapters. Specifically, a statistical representation of heterogeneous populations can be obtained by defining suitable probability distributions on the flux polytope. This Chapter addresses

* the different sources of variation that affect microbial metabolism along with the mechanisms that may favor higher variability,

* the methods devised to represent heterogeneous microbial populations within the framework of constraint-based models, and

* how these approaches connect to the optimality scenario presented in the previous Chapters.

11.1. Introduction

The theory of cellular metabolism developed up to this point through constrained-based models (CBMs) relies crucially on some type of optimality assumption: among all viable flux states encoded in the flux polytope by mass-balance, thermodynamic and regulatory constraints, cells strive for those that maximise a physiologically motivated objective function. For *E. coli* cells growing on carbon-limited substrates, for instance, it is reasonable to take such a function to be the (specific) growth rate. At the very least, these optimal states provide reference points to gauge cellular behaviour. In this respect, having a good grasp of what makes a configuration of fluxes through the network 'optimal' with respect to a certain objective is rather important from a theoretical viewpoint. On the other hand, it is not easy to prove directly in an experiment that a certain function is *actually* being optimised (in any physical system, let alone in a microbe or a microbial population). The optimality assumption can usually be corroborated *a posteriori*, e.g. by comparing optimality predictions to experimentally measured fluxes or growth rates [384, 385], or indirectly, e.g. by showing that, in a given growth medium, certain metabolic enzymes are expressed at just the level ensuring maximal growth [386]. By looking
# Chapter 9

## Universal features of autocatalytic
systems

David Lacoste, Barnabé Ledoux

### Chapter overview

* A stoichiometric theory of autocatalysis is outlined, which is based on the notion of productivity (either economical or chemical). The framework is applied to the von Neumann universal constructor model as an example. New methods of identifying autocatalytic subnetworks in complex chemical networks follow from this approach.

* The expanding economic model of the von Neumann model is a linear model of a circular and productive economy also based on stoichiometric considerations. As shown in [230], a notion of growth factor introduced in this model has valuable applications for the characterization of autocatalytic chemical reaction networks.

* A special case of the von Neumann model is the Leontief model, from which the Leontief production function can be derived. This approach is useful in modeling certain features of cell metabolism, such as cell growth laws and inhibition of bacterial growth by various types of antibiotics.

### 9.1 Stoichiometric versus dynamical autocatalysis

#### 9.1.1 Stoichiometric autocatalysis

According to IUPAC, "an autocatalytic reaction is a chemical reaction in which a product (or a reaction intermediate) also functions as a catalyst. In such a reaction, the observed rate of reaction is often found to increase with time from its initial value". While this definition provides a sound kinetic characterization of autocatalysis, it is not easy to use it to identify autocatalysis in situations in which kinetics is poorly known. Such situations arise frequently when trying to analyze complex chemical mixtures, such as astrophysical samples analyzed by researchers studying the origin of life research, or man-made prebiotic systems such as complex interacting RNA networks, in which most of the species and reactions are unknown.

On the theory side, the concept of autocatalytic sets was introduced by S. Kauffman in 1971, and played an important role in his early investigations of the Origin of order in living systems. In 2004, W. Hordijk and M. Steel expanded this original work by introducing the concept of reflexively autocatalytic food-generated networks (RAFs), namely self-sustaining networks that collectively catalyze all their reactions using only compounds from the food sets [231]. This formalism is based on the assumption that any compound (or a fraction of them) involved in randomly picked reactions has a certain probability to be catalytic [232].
at the behaviour of individual cells in a population, however, one cannot help but notice a salient feature: their diversity. Individual cells are macroscopically heterogeneous in terms of parameters like interdivision times, elongation rates, sizes at birth or division, etc. This suggests that a corresponding diversity is present at the level of intracellular processes like cell cycle, gene expression and, of course, metabolism. Quantitative experiments probing populations at single-cell resolution (see Experimental Methods Box 11.A) can nowadays characterise such a diversity in some detail. Among the remarkable outcomes of these studies is that, when analysed with a lens that accounts for diversity, bacterial growth displays signatures of universality [387, 388, 389], suggesting the existence of general, system- and condition-independent control mechanisms (e.g. of cell division and growth) that do not change with specifics like strain, quality of medium, etc. Identifying these mechanisms yields robust insight (and predictive capacity) into the physiology of microbial systems (see also Chapter on *Control of cell division and coordination with other cell-cycle processes*).

Experimental methods 11.A: Quantitative methods for single-cell analysis of microbial systems

At the very minimum, quantitative experimental characterisation of cell-to-cell diversity in microbial populations requires (i) the possibility of achieving steady-state cell growth in controlled environments, and (ii) the possibility of identifying individual cells within a population. The two setups that are most important for the present Chapter (and most widely used in general for the study of cell-to-cell heterogeneity in microbial systems) are the following.

* High-resolution optical microscopy of bacteria growing on agarose pads. Optical microscopy is the first and still most used technique to address cellular individuality [390]. Besides giving direct information about the macroscopic growth dynamics of individual cells [391, 389], it can be used in conjunction with gene expression reporters like fluorescent proteins to quantify diversity in gene expression levels [392] and dynamics [393]. Optical means usually allow to reliably follow the expression of a relatively small number of genes. In addition, however, they can also provide information about many other aspects of bacterial physiology, like motility, chemotaxis or the spatial self-organisation of colonies.

* Microfluidic 'lab-on-a-chip' devices. In essence, these techniques allow to confine single cells or small lineages thereof in controlled environments for long-term data acquisition [394]. A well-known example is the 'mother machine' [395]. In a mother machine cells grow in narrow (ca. 1 µm) microfluidic dead-end channels such that (a) all cells in the same channel are daughters of a mother cell stuck at the closed side of the channel; (b) a main feeding channel carries away cells that grow out of the length of dead-end channels (which suffice to contain a few cells, usually 5 to 10); and (c) nutrient in-flow and waste out-flow from the feeding channel ensure a constant medium in all dead-end channels via diffusion. This setup effectively keeps the population size fixed. Growing bacteria can then be imaged and analysed by standard means like time-lapse microscopy to obtain the statistics of quantities like the interdivision time or the size at birth at stationarity [388].

The setup of mother machines has the advantage that cells can be followed for many more generations than on agarose pads, since the latter tend to become overcrowded after a limited number of rounds of divisions. On the other hand, agarose pads offer a more natural environment for cell division.

In addition to these, a host of other techniques are being increasingly refined and used to probe single-cell properties and behaviour in bacterial populations, including single-cell metabolomics by mass spectrometry [396], nanoscale secondary ion mass spectrometry (nanoSIMS) [397], and single-cell transcriptomics [398].

It is not hard to guess why a bunch of identical cells sharing the same medium would, say, elongate at different rates. For one, gene expression has a stochastic component, from e.g. the random diffusion of transcription factors to targets to the thermal noise driving the on/off dynamics of transcription events. We also know that the cell cycle can be highly variable [399]. And other 'natural' sources of variance can be found in the dynamics of expression in genetic circuits, ageing, asymmetric partitioning of cellular resources at division, inter-cellular interactions, and epigenetic modifications [400]. In other terms, a degree of variability across a population is to be expected. The question, however, is, how can variability be reconciled with the optimality picture? And related to this: can we explain cell-to-cell differences in terms of some other, perhaps more involved, optimality criterion? Are there cases in which variability is optimised? Can we describe quantitatively a microbial population in ways that account for inter-cellular diversity? Note that cell-to-cell
Although very nice results follow from this assumption, such as the existence of a phase transition controlled by the connectivity of the network, this assumption is a bit problematic, because a given species can act as a catalyst or not depending on the presence of other molecules and depending on the reactions it is part of. In other words, the probability for a reaction to be catalytic is context dependent and strongly constrained by the topology of the network itself. To address both issues, namely the lack of available data on the kinetics and the shortcomings of the RAF formalism, a better starting point is to define autocatalysis from stoichiometry rather than from kinetics. While alternate definitions of autocatalysis are possible [233], we now detail the definition of [234], which encompasses autocatalytic sets and RAFs as particular cases.

Given a complete chemical network, autocatalysis is defined at the level of a subnetwork of stoichiometric matrix $N$, i.e. for a subset of species and reactions of the full network. Species which are not part of the autocatalytic subnetwork can still play an important role, for instance, food species or building blocks. The autocatalysis of this subnetwork requires two essential properties: *autonomy* (i) and *productivity* (ii). The *autonomy* condition is there to exclude direct injection or loss of species from the environment within the set of autocatalysts. More precisely, this condition requires that each species is produced by at least one reaction (each row of the stoichiometric matrix has at least one positive entry) and each reaction must at least consume one species (each column has at least one negative entry). This notion of autonomy is a close analog of the notion of circular economy in the economic context, provided that species are replaced by goods, and reactions are replaced by industry sectors of an economy. *Productivity* (ii) means the absence of mass-like conservation laws and is related to the notion of a productive economy. More precisely, this condition requires the existence of a non-zero reaction vector $v$, such that

$$ \Delta n = N \cdot v > 0, \qquad (9.1) $$

element-wise. The existence of the flux vector $v$ guarantees that there is a set of species and reactions such that all species of the set are produced by the autocatalytic subnetwork. The proof that this property is equivalent to the absence of mass-like conservation law is the content of Gordan's theorem of linear algebra, which is explained in the box 9.A.

A third more technical condition introduced in Ref. [234] is that $N$ should be *non-ambiguous* (iii) which means that a species can not be both a reactant and a product of the same reaction. This condition (iii) is less essential but is convenient as it ensures that catalytic steps can be distinguished at the level of stoichiometric matrix. Indeed, otherwise the stoichiometric matrix would be ambiguous in the sense that it would not be possible to separate the contribution of the consumption of reactants from the production of products. In practice it is always possible to transform an ambiguous reaction into a non-ambiguous one provided intermediates in the reaction are added. Another way to get around the issue of ambiguity without having to transform the network, is to define the chemical network from the start by two stoichiometric

Σ Math box 9.A Gordan's theorem

Gordan's theorem is the following result of linear algebra which takes the form of an alternative:

$$ \exists v \text{ s.t. } N \cdot v > 0, \qquad (9.2) $$

or

$$ \exists \rho > 0 \text{ s.t. } N^{\top} \cdot \rho = 0. \qquad (9.3) $$

The first side of the alternative in 9.2 corresponds to the stoichiometric definition of autocatalysis, the second side of the alternative in 9.3 corresponds to the existence of a so called mass-like mass conservation law, i.e. a conservation law with only strictly positive entries, in which case no autocatalysis is present. Note that when autocatalysis is present, only mass-like conservation are forbidden but not general conservation laws in which the entries of $\rho$ are positive and negative. An example of an autocatalytic network with non-mass like conservation law is given in [235], which also provides a nice geometric interpretation of conservation laws as manifolds. Note also that a stronger condition for autocatalysis is the absence of any conservation law, which mathematically means $ker N^{\top} = 0$.
<figure><img src="image_73.png" /><figcaption>Figure 11.1: One versus multiple optima in the flux polytope. (a) A two-dimensional flux polytope (shaded area) with non-negative fluxes and the defining constraints shown as dashed green lines. (b) The linear objective function represented by the red line has a unique maximum (magenta dot). (c) The linear objective function represented by the blue line has infinitely many maxima, lying on the highlighted magenta segment.</figcaption></figure>

variability is inherently a population-level concept. Addressing it therefore requires a framework that is capable of clearly distinguishing single-cell properties from population-level ones.

It is definitely possible to explain cell-to-cell variability in the optimality framework (see Chapter on *Solutions of constraint-based metabolic models*). For example, one could say that, in appropriate conditions, all microbes in a population are optimal, but the optima are slightly different for different cells. As a matter of fact, optimal states in CBMs need not be isolated points belonging to the flux polytope. There can in fact be infinitely many flux vectors that maximise an objective function (this happens, for instance, when an objective function attains its maxima on one of the edges or faces of the polytope, see Figure 11.1). This implies that identical cells subject to the same constraints and sharing the same objective may end up having different metabolic profiles despite carrying the same value for the objective function. In this scenario, diversity is induced by a very special feature of the objective function and, unless some other ingredient is brought into the game to lift the degeneracy, all optimal states would be equally likely for cells. If having an objective function of this type seems unlikely in a high-dimensional setup such as metabolism, one may imagine a scenario in which all cells optimise the same objective but with slightly different constraints (i.e. in a slightly different polytope, e.g. due to small variations in regulatory constraints, energy demands, or nutrient uptakes). In this case, each cell would solve its own optimisation problem, ending up having, along with a different metabolic profile, a slightly different value of the objective function. Metabolic diversity is therefore induced by variability in the constraints. But it is also possible that, if cells are subject to fluctuating exogenous constraints (e.g. variable nutrient levels), they would prefer to maximise their, say, growth rate *averaged over conditions*, especially if fluctuations occur on faster timescales than those over which metabolic reactions equilibrate. In such a case, the average growth rate would be maximum (given the external variability), but other than that every cell could carry a different growth rate and a different metabolic profile. In this respect, one can say that diversity is now being optimally adapted to external conditions. Or one may even think that different cells have different objective functions. This scenario, possibly unrealistic for growing microbial populations but not for other cell types (think for instance of the mixture of neurons with high energy demands and glia with low energy demands in the brain), would also lead to heterogeneous flux profiles and objectives. And so on.

It is clear from these examples that, in order to represent heterogeneity within CBMs, one must, first and foremost, clarify the origin of heterogeneity as much as possible. Next, it is necessary to shift from the language of individual flux vectors belonging to the flux polytope to that of *ensembles* of flux vectors or, more reasonably for large populations, of *probability densities* defined on the flux polytope. This transition is less trivial and more momentous than it sounds and, together with the causes of variability, is the core subject of the present Chapter. We shall begin by giving a more precise characterisation of the different types and sources of diversity that can be considered when modelling metabolic networks. Next, we shall introduce probability densities on the flux polytope and briefly discuss a few simple examples. We shall then address the general problem of using probability densities to model heterogeneity and uncertainty, most notably that seen in empirical data. Finally, we will show how these ideas can be used to generalise the notion of
<figure><img src="image_101.png" /><figcaption>Figure 9.1: Five minimal autocatalytic motifs (figure taken from [234]). Yellow circles represent species, black lines connecting them represent reactions and the orange squares indicate the locations where further reactions could be added while preserving the motif type.</figcaption></figure>

metric matrices instead of one as we will do in section 9.2.

Remarkably, this mathematical definition is enough to guarantee the existence of a small number of minimal autocatalytic motifs called autocatalytic cores. The minimality of these cores means that they cannot contain smaller cores in them. In [234], it was found that with the above assumptions, only five minimal motifs could exist, which are represented in Fig. 9.1. It also follows from that construction that an orientation of the reactions in a core must exist such that v has only non-negative components [236].

As a simple illustration of this framework, let us consider the universal constructor which von Neumann introduced in 1940 [237]. This is an idealized machine U, that would be able to construct any object including itself when given some set of instructions I. von Neumann identified a potential recursion issue related to the self-replication of the machine together with its instructions and he also understood that for such a machine to evolve without compromising its replication, an additional player in addition to U and I was needed. This additional player would be a universal copy machine X that would copy the instructions without translating them [238]. With this remarkable insight, von Neumann foresaw the essential mechanism of the DNA based translation-transcription machinery that we know today. We can summarize the reactions in which the universal constructor U, the universal copy machine X and the instructions I are involved by the following simple chemical network :

$$ 
\begin{align*}
& X + I \rightarrow 2I + X, \\
& U + I \rightarrow 2U + I, \\
& U + I \rightarrow U + I + X.
\end{align*}
\qquad (9.4) $$

When written in this way, the sub-network of these three reactions with the set {$U, I, X$} contains several catalytic steps, so that neither conditions (i) nor (iii) are satisfied. To solve the issue, let us introduce food species $F_1$, $F_2$ and $F_3$, intermediate species $XIF_1$, $UIF_2$ and $UIF_3$ and their corresponding reactions so that the network is no longer explicitly catalytic :

$$ 
\begin{align*}
& F_1 + X + I \rightarrow XIF_1, \quad XIF_1 \rightarrow 2I + X, \\
& F_2 + U + I \rightarrow UIF_2, \quad UIF_2 \rightarrow 2U + I, \\
& F_3 + U + I \rightarrow UIF_3, \quad UIF_3 \rightarrow U + I + X.
\end{align*}
\qquad (9.5) $$

Now the subnetwork of species {$U, I, X$} satisfies the conditions (i) (autonomy) and (iii) (non-ambiguity), and also condition (ii) (productivity) because the reaction vector that corresponds to summing all the reactions produces all the species of the set, according to the overall reaction $X + 2U + 3I \rightarrow 2X + 3U + 4I$. Therefore this network is autocatalytic from the point of view of stoichiometry. Note also that all reactions have been assumed to be irreversible here, which is allowed within the framework of [234], since thermodynamic compatibility is not considered.

A significant benefit of the stoichiometric definition of autocatalysis outlined above is that linear programming algorithms, which are a classic tool of analysis of optimization problems, can be used to search effectively for autocatalytic subnetworks within large chemical networks [239, 230, 236]. In carrying out this
optimality to heterogeneous populations.

## 11.2. Sources of variability and uncertainty in metabolism

Metabolic heterogeneity is widespread among clonal populations of prokaryotic and eukaryotic cells. Populations of *Escherichia coli* display diverse cell-to-cell conversion yields of glucose into final products, such as fatty acids and tyrosine [401]; not surprisingly, the intracellular concentration of co-factors, including ATP, also vary significantly between cells [402]. *Saccharomyces cerevisiae* metabolic states have been observed to divert over time for each cell. For instance, a single budding yeast does uptake oxygen before the duplicating its genetic material, but it changes to an anaerobic metabolism once DNA duplication starts to prevent mutations related to free-radicals [403, 404]. Animal cells within a single tissue also show heterogeneous metabolisms. Non-small cell lung cancer show a remarkable diversity of preferred carbon sources. Within the tumour, some cells consume glucose and produce lactate, whereas others divert their metabolism to consume lactate as a carbon source [405].

The root cause of this metabolic heterogeneity is manifold, including uneven distribution of nutrients in the environment, asymmetric cell partitioning at division, and noise in gene expression [406, 407]. These effects are stochastic, preventing the determination of each cell's metabolic state in advance. This type of uncertainties are rooted in the nature of metabolism itself, thus we say it is objective uncertainty.

There is another type of uncertainty, the one that comes from our models of metabolism. There can be missing reactions in the metabolic network reconstruction [408], some reactions' directionality can not be anticipated under *in vivo* conditions [409], and experimental errors in the techniques used to measure metabolic quantities –such as exchange fluxes, or the weights of the biomass reaction– does not provide precise values but approximated intervals [410]. Even when using a *bona fide* metabolic network conditioned by precisely measured parameters, the use of optimality principles –e.g. maximisation of biomass growth– are typically implemented as under-determined mathematical formulations, whose output is not a single flux vector but rather a reduction of the viable polytope [411] (see Chapter on *Solutions of constraint-based metabolic models*). This is exemplified in Fig. 11.2.A, where the maximisation of the network fourth flux, $v_4$, can only reduce the fluxome space to the subspace defined by the line shown in Fig. 11.2.B. Uncertainties that are the result of modelling uncertainties can be categorized as subjective, as they arise solely from an observer's lack of knowledge.

**Problem 11.1** Use the model in Fig. 11.2 to study different objective functions, specifically combinations of fluxes. Can you find cases in which the optimum is not unique?

Although objective and subjective uncertainties arise from different sources, both can be modelled using probability theory as we will see in the following section.

## 11.3. Probability densities over the flux polytope

In what follows, we shall denote the convex flux polytope by $\mathcal{P}$ and a generic flux configuration in $\mathcal{P}$ by $\mathbf{v} = \{v_i\}_{i=1}^N$. A probability density $p$ defined on $\mathcal{P}$ is any non-negative function such that

$$ \int_{\mathcal{P}} p(\mathbf{v})dv_1 \cdots dv_N \equiv \int_{\mathcal{P}} p(\mathbf{v})d\mathbf{v} = 1. \qquad (11.1) $$

Notice that the integral over $\mathcal{P}$ implicitly encodes two types of constraints: mass-balance equations (i.e. $S\mathbf{v} = 0$) and ranges of variability of the form $v_{i,min} \le v_i \le v_{i,max}$ (see Chapter on *Solutions of constraint-based metabolic models*). The quantity $\int_{\mathcal{P}} d\mathbf{v}$ represents therefore the *a priori* volume of $\mathcal{P}$ (which, understandably, is far from simple to calculate for high-dimensional polytopes like those corresponding to genome-scale metabolic network reconstructions). As usual, $p(\mathbf{v})$ can be interpreted as the relative likelihood of flux configuration $\mathbf{v}$: if we imagine that a cell is assigned a flux configuration by “randomly sampling it from $\mathcal{P}$” using the rule described by $p$, then $p(\mathbf{v})d\mathbf{v}$ represents the probability
program, the authors of [230] found that the number of autocatalytic subnetworks typically grow exponentially with system size. In practice, a much smaller number of subnetworks is expected to be relevant dynamically in large networks, an issue which we address now.

### 9.1.2 Dynamical autocatalysis

Historically, autocatalysis as studied in classic chemistry has been related to a certain type of kinetic pattern. To distinguish this classic definition from the previous based on stoichiometry, we will call this form of autocatalysis, dynamical autocatalysis. In this view, dynamical autocatalysis is defined as a chemical process in which one of the products catalyzes its own formation according to

$$ \frac{dx_i}{dt} = k(\mathbf{X}) \cdot x_i^n + f(\mathbf{X}), \quad \text{for } k > 0, n > 0, |k| \gg |f|, \qquad (9.6) $$

where $\mathbf{X}$ is the vector of all the concentrations $x_i$, the term $k(\mathbf{X}) \cdot x_i^n$ describes the contribution from autocatalysis while the function $f$ describes the contribution coming from the rest of the chemical system [240]. It follows from this definition that diverse forms of autocatalysis are possible depending on $k(\mathbf{X})$ and $n$, and that these forms could remain concealed if $f(\mathbf{X})$ is too large. It also follows from this definition that when $k(\mathbf{X})$ is constant, autocatalysis can generate exponential growth for $n = 1$, but also over-exponential $n > 1$ or sub-exponential $n < 1$ evolution. In practice, the exponential growth regime is anyway limited to an intermediate time window either because the reaction eventually runs out of substrates (at long times) or due to product inhibition, which can trigger sub-exponential behavior. Further, $k$ and $n$ are not independent factors since they arise from shared physical factors (availability and diffusion of ligands, size and flexibility of the molecules...) resulting in a trade-off between these two parameters that is relevant for designing catalysts or autocatalysts from bottom up [241].

From the definition in Eq. (9.6), it follows that dynamical autocatalysis has the potential to destabilize a dynamical state, which would otherwise have remained stable. Recently, the connection between the topology of the autocatalytic reaction network and its dynamical stability has been explored in two separate works that address different sides of that issue. In the first one, it was proven that for fully connected dilute systems with no degradation, the stoichiometric definition of autocatalysis leads to dynamical autocatalysis, characterized by a strictly positive Lyapunov exponent [242]. In the second one, for a certain class of parameter-rich kinetics, it was shown that the stoichiometric definition of autocatalysis implies a choice of reaction rates for which an unstable fixed point necessarily exists [243]. In a nutshell, the first work provides explicit results regarding the relation between stoichiometric and dynamical aspects of autocatalysis but the results are limited to the diluted regime, while the second work does not have this limitation, but is only a proof of existence: it does not provide an explicit method to obtain the reaction rate stated in the result.

Several recent studies have explored growing systems from the point of view of non-equilibrium thermodynamics [244]. In particular, in [235], Kamimura et al. have built a comprehensive chemical thermodynamic theory of open systems which are also self-replicating. This approach clarifies the thermodynamic conditions under which growth is possible in a system in which the volume is also growing [235]. To include this change of volume, these researchers developed an extension of traditional chemical thermodynamics theory. The growth of the volume is an important feature of autocatalysis, which manifests itself in certain experiments such as that of [245]. In this work, small-molecule autocatalytic reactions occur in compartments made of water in oil droplets. Small molecules, which act as fuel in these reactions, can diffuse between compartments while the large molecules which are produced inside the compartments cannot diffuse across compartments. This work provides a stunning demonstration that autocatalysis can drive compartment growth, competition and reproduction.

The question of how to connect stoichiometric, kinetic and thermodynamic features of autocatalysis is an ongoing active area of research which is pursued by several groups theoretically [236, 246, 235, 247]. For instance, in [246], a framework has been proposed to derive structural and thermodynamic bounds for autocatalytic chemical networks assuming mass-action law kinetics. The term structural means that these bounds depend on the topology of the network but not on the value of the rate constants. Another
<figure><img src="image_74.png" /><figcaption>Figure 11.2: Minimal metabolic network with multiple optima. Toy network where the top metabolite has one producing reaction flux $v_1$ and two reactions fluxes, $v_2$ and $v_3$, that convert it into the lower metabolite, which is excreted via $v_4$. When fixing $v_1 = 10$, the maximisation of $v_4$ –a proxy of biomass growth rate– results in $v_4 = 10$. Even in this case, there are infinitely many points that are coherent with this solutions, among them the colour-coded as red blue and green (A). The subspace of solutions form a line (B).</figcaption></figure>

that the cell's flux configuration will lie in a small volume $dv$ around $v$. It is clear then that probability densities on $P$ provide a mathematically convenient way of describing the metabolic state of large populations (or ensembles) of cells at a given time, provided one can assume that cells have the same metabolic network and are subject to the same constraints, so that $P$ is the same for all of them. For the population of cells described by $p$, the probability density clearly contains all the statistics of metabolic fluxes, from mean values to variances to correlations. For instance, by integrating $p$ over all fluxes except the $i$-th, one obtains the marginal probability density of flux $v_i$, i.e.

$$ \int_{P} p(v) dv \backslash i = p_i(v_i) , \qquad (11.2) $$

where the index $\backslash i$ corresponds to 'except for the flux of index $i'$ (so $dv \backslash i = dv_1 \cdots dv_{i-1} dv_{i+1} \cdots dv_N$). And from $p_i$ we can immediately retrieve the statistical features of flux $v_i$ (e.g. mean value, variance, etc).

Let us make a few simple examples.

* If we assume that all cells in the population maximize the same objective function, and that there is no degeneracy in the optimal state, then

$$ p(v) = \delta(v - v^*) , \qquad (11.3) $$

where $v^*$ denotes the (unique) objective-maximising flux vector and $\delta(x)$ denotes Dirac's $\delta$-distribution.

* If we can make no assumption on the cells' metabolic activity other than it has to be compatible with the constraints encoded by $P$, then any flux vector $v \in P$ is equally likely to occur in a population. This means that $p$ is constant on $P$. Specifically, its value must be equal to the inverse of the volume of $P$:

$$ p(v) = \left( \int_{P} dv' \right)^{-1} \quad (v \in P) . \qquad (11.7) $$

For any given flux polytope, this distribution can be sampled in practice using the methods described in the Chapter *The space of metabolic flux distributions*.

* Imagine having a dataset derived from a $^{13}$C labelling experiment (mass spectrometry) that gives the mean value $\bar{v}_i$ of every flux in the network (the average being over the population of cells used in the experiment), together with an experimental error $\sigma_i$ (which likely conflates different sources of uncertainty of which we may know very
important structural property of biochemical networks is for instance robust perfect adaptation (RPA). This property means that there exist some subnetworks, with specific topological features, whose parameters are irrelevant to the steady-state properties of the rest of the network [248].

In [236], the authors proved there always exists a well-defined CRN corresponding to an autocatalytic core, where by well-defined, we mean that there exists a list of chemical species with finite concentrations and a list of reversible chemical reactions obeying mass action law kinetics that realize this network dynamically. Within these assumptions, they also showed that thermodynamic constraints prevent certain associations of autocatalytic cores, which suggests that only a restricted number of cores is relevant for the dynamics of a given autocatalytic network.

## 9.2 von Neumann's model of an expanding economy

Besides his theory of the universal constructor mentioned above, J. von Neumann made another essential contribution to our topic by proposing in 1945 a model for an expanding economy. The model assumes the economy to be circular, which means that products (or goods) are produced from other goods and from building blocks using a number of processes with a certain intensity $v_j$ [249]. There are $n$ goods and $m$ processes with $n < m$, which are characterized by constant ratios of inputs to outputs. The model is formulated in terms of an output matrix $B$ and an input matrix $A$. Since the total amount of good produced must match the internal and external demand, described by the positive vector $d$, we have the equation

$$B \cdot v = A \cdot v + d. \qquad (9.7)$$

This economic model can be directly mapped onto a chemical reaction network, if goods are interpreted as chemical species, the vector $v$ represent chemical fluxes [250] and $d$ could represent a dilution or degradation. To formalize this analogy, it is convenient to split the stoichiometric matrix $N$ into the part that concerns the production of products denoted $N^+$ and the part that concerns the consumption of goods or species $N^-$, so that we can use $B = N^+$, $A = N^-$ and $N = N^+ - N^-$. We say that an economy is productive when there exists a non-negative vector $v$ such that $B \cdot v > A \cdot v$. This condition maps exactly to the notion of productivity introduced in Eq. (9.1) for stoichiometric autocatalysis, while the condition of circular economy maps to the notion of autonomy introduced at the same time.

Obviously, we must require the positivity of the vectors $v$ and the condition $\sum_i v_i > 0$. Another condition is that the economic system is irreducible, which means that it cannot be decomposed into isolated independent sub-parts (a sub-part being here a subset of goods which does not require goods from outside the subset to produce all the goods of the subset). To show the existence and unicity of the dynamic equilibrium for such economies, J. von Neumann introduced the function

$$\alpha(\mathbf{v}) = \min_i \frac{\sum_j N_{ij}^+ v_j}{\sum_j N_{ij}^- v_j}, \qquad (9.8)$$

where the minimum is taken over all goods $i$, and he proved that the function $\alpha$ is uniquely defined from the vector $v$, which is itself part of the solution. The quantity

$$\alpha = \max_{v} \alpha(v), \qquad (9.9)$$

represent a growth factor of the economic system and the value of the vector $v$ at the maximum describes the set of goods that defines a dynamic economic equilibrium.

In [250], Blanco et al. applied this framework to the study of autocatalytic reaction networks. They underlined the importance of $\alpha$ for chemical reaction networks, which they call the maximum growth factor (MGF). Note that this maximum growth factor is different from the dynamical growth rate of species within the network. The MGF does not have dimensions of a growth rate, instead it is a dimensionless factor, which can be evaluated based only stoichiometry. Blanco et al. proved that this MGF is strictly larger
Mathematical details 11.B: Dirac's δ-distribution

For our purposes, the defining property of the δ-distribution in one dimension is the following: if a variable x is
δ-distributed around the finite value x*, then, for any continuous function f,

$$
\int_{-\infty}^{+\infty} f(x) \delta(x - x^*) dx = f(x^*) . \qquad (11.4)
$$

This means that, intuitively, $\delta(x - x^*) = 0$ everywhere on the real axis except at $x^*$, where its value is $+\infty$. Such a function only makes sense under an integral sign. In this respect, (11.3) should be seen as an abuse of notation, albeit a convenient one. There are however several ways to represent the $\delta$-distribution that comply with the above requirement. For example, one can define

$$
\begin{align}
\int_{-\infty}^{+\infty} f(x) \delta(x - x^*) dx &:= \lim_{\sigma \to 0} \int_{-\infty}^{+\infty} f(x) \frac{1}{\sqrt{2\pi}\sigma^2} e^{-\frac{(x-x^*)^2}{2\sigma^2}} dx \nonumber \\
&= \lim_{\sigma \to 0} \int_{-\infty}^{+\infty} f(x^* + \sigma y) \frac{1}{\sqrt{2\pi}} e^{-\frac{y^2}{2}} dy = f(x^*) . \tag{11.5}
\end{align}
$$

The generalisation to $n > 1$ dimensions is obtained by straightforwardly assuming $\delta(\mathbf{x} - \mathbf{x}^*) = \prod_{i=1}^{n} \delta(x_i - x_i^*)$,
so that

$$
\int_{\mathbb{R}^n} f(\mathbf{x}) \delta(\mathbf{x} - \mathbf{x}^*) d\mathbf{x} = f(\mathbf{x}^*) . \quad (11.6)
$$

Because the δ-distribution effectively has non-zero probability mass only at a single point, it is reasonable to expect
(11.6) to hold also if the integral is carried out over a compact domain D, provided x* belongs to D. This is
indeed the case, although the proof requires some work. For a quick guide to the many other interesting and useful
properties of the δ-distribution that are beyond our current scopes, see [412].

little, if anything at all), such that the experimental population-level estimate of $v_i$ is $\bar{v}_i \pm \sigma_i$. Let us assume that we know enough about the experiment to be able to define a flux polytope for the cell type ($\mathcal{P}$), and that all empirically measured (averages and errors) are in $\mathcal{P}$. Then, if we want to describe the population by a density in $\mathcal{P}$ that is uniform over the domain defined by experimental estimates, we can set

$$
p(\mathbf{v}) = \prod_{i=1}^{N} \frac{\theta(\bar{v}_i + \sigma_i - v_i) \theta(v_i - \bar{v}_i + \sigma_i)}{2\sigma_i} \quad (\mathbf{v} \in \mathcal{P}) , \qquad (11.8)
$$

where $\theta(x)$ denotes the Heaviside (step) function.

**Mathematical details 11.C: The step function**

The Heaviside function is defined as

$$
\theta(x) = \begin{cases} 1 & \text{for } x > 0 \\ 0 & \text{for } x < 0 \end{cases} . \qquad (11.9)
$$

**Problem 11.2** Show that, for a real variable *x*, a continuous function *f* and upon integration over *R*, *d/dx*θ(*x*) = *δ*(*x*).
**Hint.** Use the fact that *d/dx*[θ(*x*)*f*(*x*)] = *θ'*(*x*)*f*(*x*) + *θ*(*x*)*f'*(*x*).

* (Boltzmann distribution) Let $f(\mathbf{v})$ denote a generic function of the flux vector, such as $f(\mathbf{v}) = \sum_{i=1}^{N} c_i v_i$, with $c_i$ prescribed constants. The Boltzmann distribution is defined as

$$
p(\mathbf{v}) = \frac{1}{Z(\beta)} e^{\beta f(\mathbf{v})} \quad (\mathbf{v} \in \mathcal{P}) , \qquad (11.10)
$$
<figure><img src="image_102.png" /><figcaption>Figure 9.2: Autonomous subnetworks with highest value of MGF are shown for the formose network (A) and for the *E. coli* metabolism (B). The intensity of fluxes in each subnetwork is shown with the color scale. The figure is reproduced from [230] with permission from the authors.</figcaption></figure>

than one, if the network is autocatalytic [250]. Conversely, if an autonomous network (where autonomy is now defined at the level of $N^-$ and $N^+$) has an MGF strictly larger than one, then it is autocatalytic. As an illustration, the autocatalytic nature of the network of Eq. (9.1.1) can be assessed without the need of introducing intermediates, because one can show that its MGF is strictly larger than one and that it is autonomous.

Blanco et al. also developed efficient algorithms to identify the strongest, maximal and minimal autocatalytic subnetworks, and used them to study the formose and E. coli reaction networks. The formose network is an important prebiotic network, which is known to be autocatalytic. Their results are markedly different for the two networks as shown in Fig. 9.2. From an analysis of MGF of subnetworks, they found that in the formose network, a single subnetwork dominates all the others but is fragile with respect to perturbations. In contrast, the *E. coli* network is built from interlinked cores, which together forms a robust structure.

### 9.2.1 Leontief's production function

Wassily W. Leontief won the Nobel prize in economics for his work on input-output relations in economic systems, which he started in 1936 [254]. He later developed this mathematical tool to study the American economy, and in particular the interdependency between industries.

The Leontief model is the special case of the von Neumann model, in which each productive activity has a single output (no joint products) whereas there may be many activities producing the same output and each good is produced by at least one industry. As in the von Neumann model, it is also assumed in the Leontief model that goods are produced with fixed ratios of production factors. This assumption leads to the Leontief production function discussed in the box 9.B.

We can think of this production function as the outcome of a supply chain where production factors have to be assembled in fixed proportions in order to form a product. For instance, in order to build one bike, you need two wheels, one saddle, two pedals, etc. Those ratios are fixed, or equivalently the elasticity of substitution is $\sigma = 0$ as discussed in the box 9.B. If you want to form a product *P*, for which you need to assemble $n_1$ units of $R_1$, $n_2$ units of $R_2$, ... up to $n_N$ units of $R_N$, the production rate will be limited by the smaller value of $R_j/n_j$, that is the number of sets of resource *j* required to produce *P*. If the minimum
<figure><img src="image_75.png" /><figcaption>Figure 11.3: Boltzmann distribution on the flux polytope. The Boltzmann distribution, Eqn (11.10), morphs from a uniform probability density to a $\delta$-distribution concentrated on the flux vector that maximises the function $f$ as $\beta$ varies from 0 to $+\infty$.</figcaption></figure>

where $\beta$ is a constant and $Z$ is a factor ensuring normalisation (i.e. (11.1)), namely $Z(\beta) = \int_P e^{\beta f(\mathbf{v})} d\mathbf{v}$. The behaviour of $p$ is simple to grasp in three limits.

1. For $\beta \to 0$, (11.10) reduces to (11.7): in other words, $p$ becomes uniform over $\mathcal{P}$ (and therefore insensitive to $f$).

2. For $\beta \to +\infty$, $p$ effectively concentrates on the flux vector $\mathbf{v}^*$ that maximises $f$ (which for simplicity we assume to be unique). To see this at a heuristic level, it suffices to notice that, for any $\mathbf{v} \neq \mathbf{v}^*$, the ratio

$$ \frac{p(\mathbf{v}^*)}{p(\mathbf{v})} = e^{\beta[f(\mathbf{v}^*) - f(\mathbf{v})]} \qquad (11.11) $$

increases exponentially as $\beta$ increases. Because densities are normalised, when this ratio becomes large, $p(\mathbf{v})$ must become very small. Hence, when integrated over $\mathcal{P}$, the larger is $\beta$, the closer to $\mathbf{v}^*$ must flux vectors be in order to give a significant contribution to the integral. For $\beta \to +\infty$, the only relevant contribution comes from $\mathbf{v}^*$, so that, effectively, $p(\mathbf{v}) \approx \delta(\mathbf{v} - \mathbf{v}^*)$. This conclusion can be reached more precisely using Laplace's method (a.k.a. saddle-point approximation) to evaluate integrals of the form $\int_{\mathbb{R}^n} e^{\beta g(\mathbf{x})} d\mathbf{x}$ in the limit $\beta \to \infty$ for fixed $n$ (see e.g. [413], Ch. 27).

3. By a similar reasoning, for $\beta \to -\infty$ the only relevant contribution to integrals involving $p$ comes from the (unique, by assumption) flux vector $\mathbf{v}_*$ that minimizes $f$, so that, effectively, $p(\mathbf{v}) \approx \delta(\mathbf{v} - \mathbf{v}_*)$.

When $\beta$ varies, things depend strongly on the form of $f$ and can become rather complicated when $f$ is non-linear, especially if it involves terms that involve the product of two or more fluxes ('high-order interactions'). However, in the simple case in which $f$ is linear (as outlined above), then the probability density gradually morphs from a uniform distribution over $\mathcal{P}$ to a $\delta$-distribution around the maximum of $f$ as $\beta$ increases from 0 to $+\infty$ as shown in Fig. 11.3 (and likewise when $\beta$ decreases from 0 to $-\infty$). In this respect, the parameter $\beta$ can be seen simply as a 'degree of optimisation': the closer a population is to optimising $f$, the higher the value of $\beta$. For reasons that will become clear in the next section, the Boltzmann distribution plays an especially important role in this Chapter.

**Problem 11.3 Well- versus ill-defined flux spaces.** Using the sampling methods introduced in the Chapter (The space of metabolic flux distributions) and a linear objective function of your choice, write a program that will sample a toy two-dimensional flux polytope according to (11.10), and check the outcome for a few values of $\beta$. Then try changing the shape of the polytope in different ways by changing the constraints. What features of the polytope can make sampling harder and/or less accurate (i.e. require a larger number of samples)? Can you work out a modification of the sampling algorithms that alleviates these problems?

* In Constrained Allocation FBA [385] (see the Chapter Large cell models), one considers an ensemble of growth-rate
<figure><img src="image_103.png" /></figure>

Economics analogy 9.B Production functions

In economics, production functions relate the quantities of outputs to that of inputs in a system. There are various production functions commonly used in economics, one of the most well-known is the Cobb-Douglas production function [251, 252], which can be written in a general way:

$$y = c \prod_{i} x_{i}^{\alpha_{i}}, \qquad (9.10)$$

where $x_i$ are production factors (which could be labor, capital or other factors), $y$ is an amount of product and $c, \alpha_i$ are positive coefficients. Mathematically, this law bears similarities with the mass action law in chemistry, where the rate of production of a species is related to the product of the concentrations of the reactants to the power of their associated stoichiometric coefficient. A more general production function is the constant elasticity of substitution (CES) production function [253]. This function accounts for the fact that one product may be substituted with another one and can be written in the form:

$$y = c \left( \sum_{i} \alpha_{i} x_{i}^{\rho} \right)^{\frac{1}{\rho}}. \qquad (9.11)$$

Here, $\rho$ is the coefficient of substitution and $\alpha_i$ is the weight of the production factor $i$ in the total production ($\sum_i \alpha_i = 1$). Instead of $\rho$, one often uses the elasticity of substitution $\sigma$ which is such that $\rho = \frac{\sigma-1}{\sigma}$.

Here it is assumed that the elasticity of substitution is the same for all pairs of production factors. Note that $\sigma$ could also be defined as an elasticity coefficient that compares the change in the ratio of inputs to changes in the ratio of marginal products. Indeed, the marginal product with respect to factor $j$ is $y_j = \frac{\partial y}{\partial x_j} = c\rho\alpha_j x_j^{\rho-1} (\sum_i \alpha_i x_i^{\rho})^{\frac{1-\rho}{\rho}}$ and represents the sensitivity of the product to the amount of production factor $j$. Therefore, $y_j/y_k = \alpha_j x_j^{\rho-1} / \alpha_k x_k^{\rho-1}$, and we recover that:

$$\sigma = \frac{\partial \ln(x_j/x_k)}{\partial \ln(y_k/y_j)}. \qquad (9.12)$$

Therefore, $\sigma$ measures precisely how the ratio of two production factors is modified when the ratio of two corresponding marginal products is modified. Now, three famous cases can be considered [253]:

* $\sigma \to \infty$: In this case, any modification in the ratio of marginal products would require an infinite modification in the ratio of production factors. This means that the ratio of marginal products remains constant whatever the modification in the ratio of production factor. This is the so called perfect factor substitution limit. Indeed, the CES production function becomes linear ($\rho = 1$), $y = c \sum_i \alpha_i x_i$, and any production factor can be substituted by another (even if $\exists i, x_i = 0$, it can be replaced by any other $x_j$).

* $\sigma = 1$: In this case, a modification in the ratio of marginal products translates to the same modification in the ratio of production factors. If we want the dependency of the production in production factor $j$, we need to double the amount of production factor $j$. This limit corresponds to $\rho \to 0$, in which case we recover the Cobb-Douglas production function $y = c \prod_i x_i^{\alpha_i}$. This means that the level of substitution of a any production factor is null (if one of the $x_i = 0$, then $y = 0$). Interestingly, for all values of $\sigma \in [0, 1]$, there is no possible substitution between production factors.

* $\sigma = 0$: This case is particularly interesting, it means that the ratios of the production factors $x_j/x_k$ are fixed whatever the ratios $y_j/y_k$. In this limit we recover the so called Leontief function with fixed factor proportions. There is still no substitution between the products (if one $x_i = 0$, then $y = 0$), but this also means that forming one unit of product always requires the production factors in the same proportions. This corresponds to $\rho \to -\infty$, which gives $y = c \min(x_i)$. Interestingly, this is a way to get rid of the coefficients $\alpha_i$ which represent the share of each production factor in the total production. This boils down to saying that production is limited by the scarcest resource.

time to produce one unit of $P$ is $\tau_P$, and the minimum time to use resource $R_j$ in order to produce one $P$ is $\tau_i$, we can write

$$\frac{dP}{dt} = \frac{1}{\tau_P} \min \left( \frac{\tau_P R_1}{\tau_1 n_1}, \frac{\tau_P R_2}{\tau_2 n_2}, \dots, \frac{\tau_P R_N}{\tau_N n_N} \right). \qquad (9.13)$$

Indeed, $\tau_P/\tau_j n_j$ is the number of products which can be simultaneously produced from one resource
maximisation problems constructed by sampling (from a prescribed probability density) a family of random variables representing the proteome fraction to be invested in each metabolic enzyme per unit flux of the corresponding reaction. The idea in CAFBA is that different sets of parameters effectively correspond to different cells, reflecting the cell-to-cell variability in e.g. transcription levels and protein abundances. The population-level behaviour is then obtained by averaging over different choices of these parameters (i.e. over a population of heterogeneous cells). An alternative interpretation is however possible, namely that different parameters reflect the different environmental conditions that a species can encounter over its life process history. By averaging over parameters one obtains a growth strategy that levels out this environmental variability. Such a strategy may be the one that cells prefer to implement e.g. when environmental fluctuations are fast (faster than regulatory timescales). In either case, in CAFBA, randomness in a family of parameters related to the optimisation problem induces randomness in the solutions, and therefore a probability density over the feasible space. This probability is unfortunately hard to write down explicitly in the case of CAFBA due to the complexity of the optimization problem. Its marginal distributions are however easy to calculate numerically. Two of them, specifically for the single-cell growth rate and acetate excretion fluxes, are shown in Fig. 2 in [385].

We could provide more examples but the key message of this section should already be visible: probability densities on the flux polytope are useful (a) when one wants to explicitly represent how uncertainties, experimental knowledge (with errors), or variability in parameters impact our knowledge of what part of the flux space $P$ is occupied by the metabolic states that occur in a true microbial population; and (b) when one is interested in representing an *optimal* (in some sense) population in a way that explicitly accounts for heterogeneity. If one has data (with errors), a probability density can provide a representation of the data, as in (11.8). It can likewise describe the solution to a population-level optimisation problem, and therefore a purely theoretical prediction, as in (11.3). Or the solution to an optimisation problem with uncertainty, i.e., partial knowledge or variability in some of the parameters, in which case it represents an ‘informed’ theoretical prediction (as in the CAFBA example, where the ‘information’ injected into the problem comes from the probability density from which parameters are sampled). Or it can simply be a tool to interpolate between extreme cases when we are unsure about how well a certain function is being optimised (as in (11.10)). Notice how, in our examples, different motivations activate different theoretical routes, all of which lead to working with probability densities that have a *priori* different origins and meanings even though they can be formally the same.

The two broad motivations for working with probability densities on $P$ outlined above [i.e. (a) representing uncertainty and (b) representing optimal heterogeneous populations], pose fundamentally different modelling challenges. In the first case, the key question is one of model selection: given some empirical knowledge, what is the probability density on $P$ that best represents our residual uncertainty? For instance: how good of a choice for $p$ is (11.8) given the data we had? Are there criteria that can guide our choice of a probability density? We will briefly consider these issues in the upcoming Sec. 11.4. When attempting to model optimal heterogeneous populations at the theoretical level, instead, one basically has to generalise the problem tackled by CBMs models like FBA to the case in which an optimal probability density is searched for instead of an optimal flux configuration. We will see how this can be done in Sec. 11.5.

## 11.4. Representing heterogeneity and uncertainty

### 11.4.1. ML, MAP and Bayesian inference

We have seen that probability densities on $P$ can represent, under certain assumptions, populations of microbes whose metabolism can be described by the same flux polytope, and that different probability densities can be surmised to model the distribution of $v \in P$ when some external information (e.g. experimental data) is available. Here, we will address the following question: how can one choose the $p$ that best represents our knowledge about the metabolic state of a population in presence of these external data?

To summarise the huge and highly involved set of problems behind the above (very general) question [413] in a way that is useful for the purposes of this Chapter, we can start by assuming we have a *priori* chosen a form of $p$ that depends
j. This result holds true if resources are fully allocated to the production of P, but if several products P<sub>i</sub> need to be produced in parallel, one resource may be used by different production chains simultaneously, meaning that a fraction α<sub>i,j</sub> of total available resources R<sub>j</sub> must be used for the specific product P<sub>i</sub>, so that:

$$ \frac{dP_i}{dt} = \frac{1}{\tau_P} \min \left( \alpha_{i,1} \frac{\tau_P R_1}{\tau_1 n_1}, \alpha_{i,2} \frac{\tau_P R_2}{\tau_2 n_2}, \dots, \alpha_{i,N} \frac{\tau_P R_N}{\tau_N n_N} \right). \quad (9.14) $$

Then the prefactor before R<sub>i</sub> represents the maximal number of copies of the product that you can produce simultaneously from one unit of resource i. The fact that resources may not be substituted with other resources has important consequences for cell metabolism [255].

In the problem ??, we study a single resource - single product industry or workstation and we show that a Leontief production function emerges from mass-action law kinetics when a certain time scale separation holds. This example is important because it illustrates that the Leontief production function not only involves fluxes associated to reactions or industries but can typically also include stocks associated to goods or metabolites, which are only available in finite amounts.

Further applications of this formalism to describe for instance the ability of a metabolic network to switch from one behavior to another one (as in the crabtree or Warburg effects in biology and in the Giffen behavior in economics) is studied in [256] and in a coming chapter to be written.

### 9.2.2 Liebig's law

Interestingly, the idea that the production rate could be limited by the scarcest resource is present in the field of agronomy under the name of Liebig's law of minimum [257, 258, 259, 260]. It was initially used to describe plant growth, which requires various resources, and where it is observed that varying the amount of fully available resources did not modify the final production. This suggests that only scarce resources will limit production and translates to the principle that when a population is growing using various resources, the scarcest will set the growth rate, and the others will be consumed accordingly. This law can be used to model the growth of an organism in an environment where resources are constrained. The link between mass action laws and Liebig's law of minimum has been studied [259]. The interest of this method is to obtain equations that are easier to solve on domains where one resource is scarcest. One main difference between the Leontief production function and Liebig's law of the minimum is that for the latter, the minimum is not necessary taken between the numbers of each production factor, but between the yields of those production factors (which can be non-linear functions) [261]. In particular, this means that the rate of production is set by the minimum of the yields of each production factor:

$$ \frac{dY}{dt} = \frac{1}{\tau} \min (\{f_i(x_i)\}), \quad (9.15) $$

where f<sub>i</sub> are functions of the production factors x<sub>i</sub>. To model the requirements of plants in nutrients, yields given by Michaelis-Menten kinetics can be used [259]. Instead of directly comparing the numbers of each production factors, it consists in comparing the yields. However, the idea that one resource will be limiting remains the same.

### 9.2.3 Application to metabolism

The law of minimum was used to build simplified models of metabolism as an ensemble of coupled autocatalytic cycles [255]. Metabolism is seen as a supply chain, where production factors must first be produced and then assembled in fixed proportions to form a product. We call P the number of proteins of one type, produced by translating mRNA (of number m) with ribosomes (of number R) and substrates (of number S). Using the Leontief production function, we can write (following the method of [255]):

$$ \left. \frac{dP}{dt} \right|_{\text{prod}} = \frac{1}{\tau_P} \min \left( \alpha_P R, \alpha_P S, \frac{\tau_P}{\tau_e s_R} m, \dots \right) \quad (9.16) $$
on certain free parameters and ask how to tailor parameters so that $p$ 'optimally' matches the empirical evidence. To be concrete, let us denote by $\psi$ the vector of parameters of $p$, and by $\mathbf{W} = \{\mathbf{w}^1, \mathbf{w}^2, ..., \mathbf{w}^R\}$ a set of $R$ experimental samples of $\mathbf{v}$. Each measurement, $\mathbf{w}$, is a vector of metabolic fluxes that ideally should include all the reactions of a metabolic network. In practice, a $\mathbf{w}$ typically spans only a subset of all the reactions of the metabolic network, those that are amenable to $^{13}$C labeling (TCA, glycolysis, and pentose phosphate pathways) or that correspond to exchange fluxes that can be reliably measured (glucose and oxygen consumption, or lactate and ethanol, to name a few). According to Bayes' rule (we assume all variables to be continuous), the quantities

* $p(\psi|\mathbf{W})$: the conditional probability density of the parameters given the observations (a.k.a. the *posterior*);

* $p(\mathbf{W}|\psi)$: the conditional probability density of the observations given the parameters (a.k.a. the *likelihood*);

* $p(\psi)$: the prior probability density of parameters (a.k.a. the *prior*);

* $p(\mathbf{W})$: the (marginal) probability density of observations (a.k.a. the *evidence*).

are related by the formula

$$p(\psi|\mathbf{W}) = \frac{p(\mathbf{W}|\psi)p(\psi)}{p(\mathbf{W})} . \qquad (11.12)$$

Ideally, what one would like to know in order to 'optimally' set the parameters of $p$ is how likely a parameter set is given the data, i.e. the full posterior $p(\psi|\mathbf{W})$, as it allows to quantify our uncertainty on the model itself. One may however also consider different (less ambitious) ways to choose parameters. The three best known methods are the following:

* Maximum Likelihood (ML) inference aims at finding the parameter vector that maximises the likelihood:

$$\psi_{\text{ML}} = \arg \max_{\psi} p(\mathbf{W}|\psi) . \qquad (11.13)$$

In standard cases, this produces a single 'optimal' vector $\psi$ (hence it is called a 'point estimator'), resulting in a $p$ that models -in a context-specific manner- the metabolic heterogeneity within the cellular population.

* Maximum a Posteriori (MAP) inference aims instead at finding the parameter vector that maximises the posterior:

$$\psi_{\text{MAP}} = \arg \max_{\psi} p(\psi|\mathbf{W}) = \arg \max_{\psi} p(\mathbf{W}|\psi)p(\psi) , \qquad (11.14)$$

where the last equality follows from the fact that $p(\mathbf{W})$ does not depend on $\psi$. Again, the MAP estimator is a point estimator.

* Bayesian inference aims finally at computing the full posterior distribution $p(\psi|\mathbf{W})$. It is therefore a 'distribution estimator' rather than a point estimator.

The following problem should clarify the way in which point estimators differ from (and are less informative than) distribution estimators in practice.

**Problem 11.4 MAP inference versus Bayesian inference.** Consider a Bernoulli random variable with parameter $\psi$, i.e. such that the probability of having $k$ successes in $n$ trials given $\psi$ is

$$p(k|\psi) = \binom{n}{k} \psi^k (1-\psi)^{n-k} , \qquad (11.15)$$

and assume that the prior for $\psi$ is a $\beta$-distribution with parameters $a$ and $b$, i.e.

$$p(\psi) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \psi^{a-1} (1-\psi)^{b-1} . \qquad (11.16)$$

Calculate the full posterior $p(\psi|k)$ and the MAP estimator for $\psi$ as a function of $k$, $n$, $a$ and $b$. Then assume $a = b = 2$ and compare the following situations: (i) a Bernoulli process that returned 2 successes in 3 trials; (ii) a Bernoulli process that returned 20 successes in 33 trials. Show that the MAP estimator for $\psi$ is 60% for both (i) and (ii) (so the two processes are indistinguishable to MAP), while the posterior is different. Knowing the posterior, which process would you pick if you were asked to point to the one that is more likely to have $\psi = 0.6$?
<figure><img src="image_104.png" /></figure>

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

Figure 9.3: (A) Scheme of coupled autocatalytic networks interacting with a toxic agent. The orange box linking two arrows represents the Leontief production function. $B_1$ represents active ribosomes; $C_1$ active RNA polymerases; similarly $B_2, ..., B_{N-1}$ and $C_2, ..., C_{K-1}$ are intermediates, $R_N, R_K$ are building blocks. We suppose that "toxic" inhibiting agents in numbers $A$ can bind to one of the autocatalysts (chosen here to be $B_1$ for simplicity). (B) Illustration of the growth laws when varying either the amount of antibiotics or the nutrient quality linked to pre-exposure growth rate $\lambda_0$ displayed on the right scale. This figure is reproduced from [263] with permission.

Ribosomes and substrates have to be used simultaneously to produce different types of proteins, $\alpha_P$ is the fraction of the total population of ribosomes (and substrates) used to produce the particular protein $P$. The minimum time to produce one $P$ is $\tau_P$, the minimum time to elongate the polymer by one amino acid is $\tau_e$, and $s_R$ is the size of the domain on the mRNA that has to be dedicated to the production of one polymer $P$ at a given time. As explained in [255], $\tau_P/\tau_e s_R$ is then the maximum number of ribosomes that can translate simultaneously one mRNA, and thus the maximum number of copies of $P$ you can produce simultaneously from one unit of mRNA. This is indeed what was predicted from Leontief's model for input-output systems: the prefactor before every amount of resources is the maximum number of copies of the protein you can produce simultaneously from this specific resource. This term is that of production of the protein, now the protein could be consumed to form a product, or degraded.

To use ribosomes, you first need to form ribosomes by assembling proteins and ribosomal RNA. Similarly, to use mRNA you must first polymerize RNA. This suggests that a minimal autocatalytic network that begins to capture the structure of the transcription-translation machinery is made of two coupled autocatalytic cycles associated respectively to RNA and ribosomes. These two cycles will then be described by coupled equations with production terms described by Leontief's production function [255].

From such an framework [255], one can derive the various growth laws that characterize the cell metabolism, which have been discussed in detail in various chapters of this book. Let us mention briefly two recent works that follow this line of research : In the first one carried out in [262], the authors noticed that RNA polymerase, and mRNA levels correlate in experiments with growth rates in contrast to the belief that ribosomes should be the sole drivers of the growth rate. To explain these observations, the authors de- veloped a theoretical framework building on [255], which account for the joint role of all these factors in the observed growth rate.

Another recent application of the above framework concerns a model for the inhibition of bacterial growth by antibiotics [263]. In that work, the cell metabolism is modeled as two coupled autocatalytic cycles, in which one cycle describes the production of ribosomes, while the other describes RNA-polymerase pro- duction as shown in Fig. 9.3. It is assumed that the antibiotic inhibits one of these two essential auto- catalytic cycles by targeting some essential metabolites in them. Growth laws can be recovered from the model as shown in Fig. 9.3B. A first law describes the increase of ribosome fraction as a function of growth rate when nutrient quality is increased (solid magenta curve) while a second growth law describes the up-
Mathematical details 11.D: Inference in metabolic network modeling

ML is the most commonly used point estimation method. As said above, the estimated parameters, $\hat{\psi}$, are computed as the argument that maximizes the likelihood of the observed data, i.e.

$$ \hat{\psi} = \arg \max_{\psi} p(\mathbf{W}|\psi) = \arg \max_{\psi} \prod_{i=1}^{R} p(\mathbf{w}^{(i)}|\psi) . \quad (11.17) $$

ML takes a familiar form if one follows, for instance, Theorell et al. [414] in modelling v according to a multivariate normal distribution:

$$ v \sim N(v|\psi) = \frac{1}{\sqrt{(2\pi)^N |\Sigma|}} e^{(-\frac{1}{2}(\bar{v}-v)^T \Sigma^{-1} (\bar{v}-v))} \quad (11.18) $$

The parameters encompass the mean flux vector, $\bar{v}$, and the covariance matrix, $\Sigma$. This is, $\psi = [\bar{v}, \Sigma]$. Accordingly, each of the i terms of the rightmost term of Eq. (11.17) can be expressed as

$$ p(\mathbf{w}^{(i)}|\psi) = N(\mathbf{w}^{(i)}|\psi) = \frac{1}{\sqrt{(2\pi)^N |\Sigma|}} e^{(-\frac{1}{2}(\bar{v}-\mathbf{w}^{(i)})^T \Sigma^{-1} (\bar{v}-\mathbf{w}^{(i)}))} , \quad (11.19) $$

so that their product can be transformed into a summation by applying the log function (without changing the final outcome):

$$ \hat{\psi} = \arg \max_{\psi} \sum_{i=1}^{R} \log[p(\mathbf{w}^{(i)}|\psi)] \quad (11.20) $$

$$ = \arg \max_{\bar{v}, \Sigma} \sum_{i=1}^{R} \left[ -\frac{1}{2}(\bar{v} - \mathbf{w}^{(i)})^T \Sigma^{-1} (\bar{v} - \mathbf{w}^{(i)}) - \log \left( \sqrt{(2\pi)^N |\Sigma|} \right) \right] \quad (11.21) $$

The second term of the summation of Eq. (11.21) is a constant and can be omitted from the maximisation. We are therefore left with

$$ \hat{\psi} = [\hat{\bar{v}}, \hat{\Sigma}] = \arg\min_{\bar{v}, \Sigma} \sum_{i=1}^{R} (\bar{v} - \mathbf{w}^{(i)})^T \Sigma^{-1} (\bar{v} - \mathbf{w}^{(i)}) . \quad (11.22) $$

Eq. (11.22) corresponds to the well-known weighted least squares estimator. The resulting $\hat{\bar{v}}$ corresponds to the mode of p, which can be interpreted as the flux vector with the highest frequency within P. With $\hat{\bar{v}}$ and $\hat{\Sigma}$, the frequency of any flux configuration can be computed from $N(v|\hat{\bar{v}}, \hat{\Sigma})$.

Standard techniques, such as confidence intervals, can be applied to assess the precision of $\hat{\psi}$. The larger the number of samples, R, the smaller the uncertainty in $\hat{\psi}$. Unfortunately, measurements of flux vectors typically include only a handful of points [415], which may lead to $\hat{\psi}$ over-fitted to the sample set. One way to overcome limited sample sizes is to regularise the estimation procedure by incorporating prior information on $\psi$ via the MAP estimation method (11.14). The evidence $p(\psi)$ in MAP can be used to encode the distribution of v values observed in previous experiments or formulated as a plausible non-informative probability distribution. For example, Heinonen et al. [416] formulated $p(\psi)$ as a multivariate normal distribution with mean values equal to zero, and variances for each flux adjusted to prevent fluxes extending beyond their lower and upper bounds defined in P. MAP estimation can be considered as an ML estimation whose objective function has been augmented by the prior distribution of $p(\psi)$. In this sense, MAP estimation is a 'regularised' ML estimation, which helps prevent overfitting.

MAP estimation however does not exploit the capacity of Bayes' theorem to explore the full set of values that the parameters can achieve. By producing a distribution estimation of the parameters, Bayesian inference allows quantifying the parameters' variability. Compared to point estimation methods, though, Bayesian inference is computationally expensive as it requires to asses how different values of $p(\psi)$ affect $p(\mathbf{W}|\psi)$. Fortunately, some families of p are susceptible to methods such as Gibbs sampling or Markov Chain Monte Carlo that offer an efficient way to compute the posterior numerically [417]. This is the case, for instance, for the truncated multivariate normal
regulation of ribosomes as a result of the inhibition of translation by ribosome inhibitors (colored solid lines). In addition, the model successfully describes the experimental dependence of the growth rate on various types of antibiotics and confirms the existence of growth bistability, namely a regime in which two possible values of the growth rates are possible in the same range of physical parameters.

## 9.3 Concluding remarks

In this chapter, we have established a connection between the universal constructor model and the expanding economic model, which were both introduced by von Neumann. Interestingly, von Neumann himself did not discuss the relation between these two works, which we make in this chapter. This expanding economic model of von Neumann (1945) and the input-output model of Leontief (1936, 1941) laid the foundation of a modern framework for economic analysis. Their framework turned out to be essential to quantify the relative interdependency of various parts of an economy and the nature and the structure of economic equilibria. Remarkably, these tools continue to inspire developments in other fields as recent studies of autocatalytic chemical networks show. Thanks to linear programming methods, we are now able to identify autocatalytic subnetworks efficiently and characterize them using notions such as the maximum growth factor or using thermodynamic bounds. More work is needed to understand the interactions between autocatalytic networks, and their role in the emergence of chemical complexity linked to pre-Darwinian evolution.

## Acknowledgements

We thank R. Pugatch, Y. De Decker, A. Despons and P. Gagrani for insightful comments and a careful reading of this chapter. We thank R. Pugatch for providing us with a nice logo for this chapter.

## Recommended readings

This book chapter is mainly based on these two recent articles: "A unifying autocatalytic network-based framework for bacterial growth laws" by Roy et al. [255], and "Universal motifs and the diversity of autocatalytic systems" by Blokhuis et al. [234].

## Problems

### Problem 9.1 Production function of a workstation with a single resource and single product

We model a workstation with a single resource and a single product [261]. Let $x$ be the amount of available resource or stock, $r$ the supply rate of this stock, $q$ the specific rate with which the stock is being lost or degraded. Further, let $u$ (resp. $v$) be the number of idle (resp. busy) machines which can process the stock. The processing time is $1/\beta$, $b$ is the number of output products per machine, $a$ is rate of capture of stock by machines, $\alpha$ is how many machines get involved per unit resource/stock.

(a) Assuming mass action law for the processing of the stock, show that the equations of the problem are

$$ \dot{x} = r - a u x - q x, $$

$$ \dot{u} = \beta v - \alpha u x, $$

$$ \dot{v} = -\beta v + \alpha u x, $$

$$ \dot{j} = bv. $$

(b) Show that the total number of machines is a constant denoted $u_0$. Derive the steady state number of busy machines $\bar{v}$ in terms of the steady state amount of stock $\bar{x}$. Comment on the form of the output production that you find.

(c) Introduce the variables $\gamma = q/(a u_0)$ and $\rho = \alpha r/(a \beta u_0)$. Derive the expression of $\bar{v}$ in the limit $\gamma \ll 1$ and $\gamma \ll |\rho - 1|$ and show that it has a Leontief form. Interpret this form in terms of how the stock is handled in the regime $\rho < 1$ and $\rho > 1$.
distributions that Heinonen et al. [416] used for the likelihood and prior functions appearing in (11.12). The posterior can then be used to derive statistical features of quantities that depend on $\psi$, e.g. metabolic fluxes.

In practice, most parameters underlying the mechanisms that govern cellular metabolism -e.g., enzymes' allosteric regulation or the local conditions within cells' organelles- remain unknown. Various hypotheses can be advanced to close this knowledge gap. Alas, it is not uncommon to have conflicting scenarios. For instance, to explain overflow metabolism in S. cerevisiae and E. coli [418, 419, 420], numerous plausible explanations have been pushed forward, including ATP savings for the production of non-oxidative enzymes (which by being smaller, compared to their oxidative counterparts, require less ATP in their synthesis) [421, 422], limited uptake rates capacity [423], and an upper limit on the dissipation of Gibbs energy [424]. (See [425] for an excellent review of optimisation-based explanations.) Because each mechanism can be encoded through a different prior, it is clear that the choice of the prior is a delicate matter in Bayesian inference. Generally speaking, the choice of the prior becomes less and less problematic the more data we have, i.e. the better sampling we have of the state space of the system. However, if data is scant, the prior will leave an important imprint on the resulting posterior. In these cases, a careful selection of the prior is paramount. Among the methods most commonly employed are (a) the construction of empirical priors (namely priors that encode previous knowledge about parameters), (b) the use of so-called “non-informative priors” (i.e. priors that reflect ‘vague knowledge’ about parameters, like the fact that a certain parameter is non-negative) [426], and (c) the selection of priors based on the Maximum Entropy principle (see below) [427, 428].

### 11.4.2. MaxEnt inference

According to the principle of Maximum Entropy (MaxEnt) [429, 430], among all probability densities that are consistent with given prior knowledge or data, the one having the largest value of the entropy

$$H[p] = - \int_{\mathcal{P}} p(\mathbf{v}) \ln p(\mathbf{v}) d\mathbf{v} \qquad (11.23)$$

is the one that best represents our knowledge about the system. A classical intuitive justification of the MaxEnt principle is most easily given for discrete variables [431].

#### Mathematical details 11.E: Entropy

Consider $N$ cells, each of which can be found in any of $K$ states (what precisely defines a state is immaterial for this reasoning). Let an assignment $\mathbf{n} = \{n(i)\}$ be given, such that $n(i)$ denotes the number of cells in state $i$ (with $\sum_{i=1}^K n(i) = N$). Because we can always exchange the states of two cells without changing $\mathbf{n}$, there are multiple ‘microscopic’ ways to realise an assignment $\mathbf{n}$. Combinatorics tells us that the number of different microscopic realisations of an assignment $\mathbf{n}$ is given by

$$\mathcal{N}(\mathbf{n}) = \frac{N!}{\prod_{i=1}^{K} n(i)!} \qquad (11.24)$$

If all $n(i)$'s are large enough, we can use Stirling's approximation ($n! \approx (n/e)^n$) to see that

$$\mathcal{N}(\mathbf{n}) \approx e^{NH(\mathbf{n})} \quad , \quad H(\mathbf{n}) = -\sum_{i=1}^{K} \frac{n(i)}{N} \ln \frac{n(i)}{N} \equiv -\sum_{i=1}^{K} p(i) \ln p(i) \equiv H(\mathbf{p}) \quad , \quad (11.25)$$

where $p(i)$ denotes the fraction of cells in state $i$ (or, equivalently for us, the probability to find a cell in state $i$). $H$ is the entropy of the assignment $\mathbf{n}$, and is in essence a measure of the microscopic degeneracy that underlies a macroscopic arrangement. The distribution $\mathbf{p} = \{p(i)\}$ carrying the largest entropy subject to certain constraints is therefore the one having the largest underlying microscopic degeneracy given those constraints. So, if one were to randomly pick a microscopic state given those constraints, the most likely macroscopic state would be the maximum
(d) Show that for small but non-zero $\gamma$, the input-output function falls below the line boundary defined by the Leontief function.

### Problem 9.2 UPF model

The UPF model is a toy model of metabolism made of two coupled autocatalytic cycles [255]. In the model, a fraction $\alpha$ of machines of type $U$ catalyze themselves. The remaining machines synthetize another type of machines $P$. The $P$ machines convert an external substrate $f$ to an internal substrate $F$, which is used by $U$ to make more copies of itself and new $Ps$. To simplify, we assume that to make one more unit of $U$ by the first reaction, one unit of $F$ is needed, and one unit to make one unit of $P$ by the second reaction.

The model is defined by the following set of equations

$$ 
\begin{aligned}
U + F &\xrightarrow{\alpha} 2U, \\
F + U &\xrightarrow{1-\alpha} P + U, \\
f + P &\rightarrow F + P, \\
U &\rightarrow \emptyset.
\end{aligned
$$

(a) Let $n_i$ be the number of molecules of type $i$, where $i \in \{U, P, F, f\}$. The life type of a machine of type $U$ is $\tau_L$, the incorporation time of $f$ into $F$ is $\tau_F$ and the incorporation time of one unit $F$ is $\tau_a$. Show that within the framework of Leontief production functions, the equations of the model are

$$ 
\begin{aligned}
\frac{dn_U}{dt} &= \frac{\alpha \min(n_U, n_F)}{\tau_a} - \frac{n_U}{\tau_L} \\
\frac{dn_P}{dt} &= \frac{(1-\alpha) \min(n_U, n_F)}{\tau_a} \\
\frac{dn_F}{dt} &= \frac{\min(n_P, n_f)}{\tau_F} - \frac{\min(n_U, n_F)}{\tau_a}.
\end{aligned}
\qquad (9.17) $$

(b) Discuss the four limiting regimes of the model: (i) $n_F \gg n_U$ and $n_f \gg n_P$, (ii) $n_F \gg n_U$ and $n_f \ll n_P$, (iii) $n_F \ll n_U$ and $n_f \gg n_P$ and (iv) $n_F \ll n_U$ and $n_f \ll n_P$. Simplify the equations for each regime by introducing a common growth rate $\mu$.

(c) Summarize your results by deriving the growth laws of the various regimes in a single plot representing $\alpha$ vs. $\mu$.

### Problem 9.3 MGF of two simple networks

(a) Calculate the MGF introduced in Eq. (9.9) for the following two simple networks:

$$ 
\begin{aligned}
A &\rightarrow B \\
B &\rightarrow 2A
\end{aligned}
\qquad (9.18) $$

and

$$ 
\begin{aligned}
2A &\rightarrow B \\
B &\rightarrow A.
\end{aligned}
\qquad (9.19) $$

(b) Comment on the values obtained in the previous question.
entropy distribution. In other terms, the MaxEnt distribution is the least biased distribution compatible with the constraints, as any other distribution satisfying the same constraints would correspond to a smaller underlying degeneracy, thereby neglecting some feasible (i.e. constraint-satisfying) microscopic configurations. In this respect, a MaxEnt distribution requires the least information besides prior knowledge (i.e. constraints). (A more detailed justification for using the MaxEnt principle as an inference tool is given e.g. in [431].)

If for instance cells are assigned to states in a completely random way, the MaxEnt distribution is the solution of

$$ \max_{\mathbf{p}} - \sum_{i=1}^{K} p(i) \ln p(i) \quad \text{subject to} \quad \sum_{i=1}^{K} p(i) = 1, \qquad (11.26) $$

which can be found via the method of Lagrange multipliers.

**Problem 11.5** Show that the solution of the above maximisation problem is the uniform distribution $p(i) = 1/K$ for all i.

If other constraints are imposed, though, the MaxEnt distribution will clearly change.

**Problem 11.6** MaxEnt distribution in different cases. Assume that a certain real variable $x$ takes values $x(i)$ in the $K$ states (one can for instance think of $x(i)$ as the growth rate of cells in state $i$). Show that the MaxEnt distributions for constraints imposed on (i) normalisation of the distribution, (ii) normalisation and mean value of $x$, (iii) normalisation, mean value of $x$ and second moment of $x$, and (iv) normalisation and mean of the logarithm of $x$, are, respectively, uniform, exponential, Gaussian, and power-law.

For our purposes, the continuous case with entropy given by (11.23) can be seen as a straightforward generalisation of the discrete one.

To get some grasp of the scenario that the MaxEnt rule provides within CBMs, let us work out one especially noteworthy case, namely the MaxEnt probability density of flux configurations with a given mean value of a generic function $f$ of the fluxes. This probability density is the solution of

$$ \max_{p(\mathbf{v})} - \int_{\mathcal{P}} p(\mathbf{v}) \ln p(\mathbf{v}) d\mathbf{v} \quad \text{subject to} \quad \int_{\mathcal{P}} p(\mathbf{v}) d\mathbf{v} = 1 \quad \text{and} \quad \int_{\mathcal{P}} f(\mathbf{v}) p(\mathbf{v}) d\mathbf{v} = \bar{f}. \qquad (11.27) $$

To find it, we construct the functional

$$ \mathcal{L}[p] = H[p] + \alpha \left[ \int_{\mathcal{P}} p(\mathbf{v}) d\mathbf{v} - 1 \right] + \beta \left[ \int_{\mathcal{P}} f(\mathbf{v}) p(\mathbf{v}) d\mathbf{v} - \bar{f} \right], \qquad (11.28) $$

where $\alpha$ and $\beta$ are Lagrange multipliers for the normalisation and the mean-value-of-$f$ constraints, respectively. Variation of $\mathcal{L}$ with respect to $p$ yields the maximum condition

$$ 1 - \ln p(\mathbf{v}) + \alpha + \beta f(\mathbf{v}) = 0. \qquad (11.29) $$

Solving for $p$ results in

$$ p(\mathbf{v}) = \frac{e^{\beta f(\mathbf{v})}}{e^{-(\alpha+1)}}. \qquad (11.30) $$

The normalisation condition however implies

$$ \int_{\mathcal{P}} e^{\beta f(\mathbf{v})} d\mathbf{v} = e^{-(\alpha+1)} \equiv Z(\beta), \qquad (11.31) $$

so that, finally,

$$ p(\mathbf{v}) = \frac{1}{Z(\beta)} e^{\beta f(\mathbf{v})} \quad (\mathbf{v} \in \mathcal{P}). \qquad (11.32) $$
<figure><img src="image_76.png" /><figcaption>Figure 11.4: MaxEnt modelling of single-cell growth rate distributions. Empirical distributions are reproduced by a MaxEnt assumption where the mean growth rate is constrained, leading to a Boltzmann distribution over the flux polytope (Eq. (11.34)).</figcaption></figure>

The value of $\beta$ must be determined from second constraint, namely

$$ \frac{1}{Z(\beta)} \int_{\mathcal{P}} f(\mathbf{v}) e^{\beta f(\mathbf{v})} d\mathbf{v} = \bar{f} . \qquad (11.33) $$

Notice that the result is nothing but Boltzmann's distribution (11.10). We have therefore found that (11.10) is the MaxEnt distribution for a given mean value of the function $f$. This means that if we have a dataset returning the empirical mean value of an observable $f$ over a population of cells, our knowledge is best represented by assuming that $p(\mathbf{v})$ is of the form (11.10), with $\beta$ ensuring the matching of empirical and theoretical means.

This suggests a possible way to represent single-cell growth-rate distributions [432], such as the E. coli populations growing in rich media studied e.g. in [433, 388, 389] (see Figure 11.4). Let us assume that all cells in the population can be described by the same flux polytope $\mathcal{P}$ and let $\lambda(\mathbf{v})$ denote the growth rate associated to flux configuration $\mathbf{v}$. We can ask the following question: what is the $p(\mathbf{v})$ on $\mathcal{P}$ that best represents our knowledge that the mean growth rate of cells is $\bar{\lambda}$ (empirical)? The answer is

$$ p(\mathbf{v}) = \frac{1}{Z(\beta)} e^{\beta \lambda(\mathbf{v})} \quad (\mathbf{v} \in \mathcal{P}) , \qquad (11.34) $$

where $Z(\beta) = \int_{\mathcal{P}} e^{\beta \lambda(\mathbf{v})} d\mathbf{v}$, and where $\beta$ is set so that the empirical mean growth rate ($\bar{\lambda}$) matches the theoretical mean, i.e.

$$ \frac{1}{Z(\beta)} \int_{\mathcal{P}} \lambda(\mathbf{v}) e^{\beta \lambda(\mathbf{v})} d\mathbf{v} = \bar{\lambda} . \qquad (11.35) $$

We can therefore solve the above equation (numerically) and analyse the resulting distribution. One sees from (11.34) that $\beta$ has a 'natural' unit given by $\lambda_{max}$, the maximum growth rate achievable in $\mathcal{P}$ (which is easily computed by LP). In the populations analysed in [432], the value of $\beta$ that ensures the matching condition ranges from $190/\lambda_{max}$ to $300/\lambda_{max}$, suggesting that indeed the degree of optimisation of $\lambda$ is significant. The most remarkable result, however is that the marginal distribution of the growth rate computed from (11.34), namely

$$ p(\lambda) = \int_{\mathcal{P}} \delta(\lambda - \lambda(\mathbf{v})) p(\mathbf{v}) d\mathbf{v} , \qquad (11.36) $$

matches the overall empirical growth-rate distributions. In other words, if one adjusts the parameter of (11.34) so that the theoretical mean growth rate and the empirical one coincide, then (11.36) reproduces the entire empirical growth-rate distribution. This observation confirms the empirical fact that the variance of single-cell growth-rate distributions is a function of the mean, such that, if growth rates are re-scaled by the mean, distributions roughly collapse on 'universal curves' that turn out to only depend on the degree of optimisation (i.e. on $\beta$). In addition, the analysis of [434] has shown that predictions for individual fluxes obtained from (11.34) (i.e. mean values plus standard deviations) provide a
Chapter 10

Resource allocation in complex cell
models

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

Chapter overview

○ Whole-cell resource allocation models on a genomic scale combine a detailed, FBA-like description of metabolism with a model of macromolecule synthesis, formulated as linear constraint-based problems.

○ Resource allocation models of cells can be built based on three basic constraints: stationary fluxes (balancing production and consumption fluxes, uptake and excretion fluxes, as well as compound dilution by cell growth); catalytic constraints relating fluxes to the amounts of catalyzing enzymes (or other machines); and density constraints, limiting molecule amounts in cell compartments, or molecule concentrations.

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

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

**10.1 Detailed resource allocation models of cells**

In the previous chapters, we have saw two principal approaches to modeling production processes in cells. To keep the number of variables low, but with intention of well-parametrizing the model, one can construct small, coarse-grained models of growing cells (Chapter 8). On the contrary, Flux Balance Analysis (FBA) models can accommodate a very large number of variables (Chapter 5), making them an excellent choice to model metabolic networks at genome-scale.

Small, coarse-grained models are very suitable for investigating base principles of life. Likely the best example to illustrate this is the work of Douwe Molenaar and co. [264], where a self-replicator model was used to proposed that the low-yield, or substrate-inefficient ("wasteful") metabolic strategies are adopted as a consequence of these pathways being more efficient in terms of protein use, compared to the high-yield pathways. In other terms, the growth output of the "wasteful" strategy per unit protein is higher than the "efficient" one. Thus we now believe that fermentation of glucose, often called under an umbrella term "overflow metabolism", will take place in many organisms if the substrate in their environments is abundant enough.

However, the chemistry of life is extremely diverse, and even such a familiar concept as fermentation can
become complicated. Take three representatives of the tree of life: a bacterium *Escherichia coli*, budding
better fit to experimentally measured fluxes than growth-rate maximising fluxes obtained from FBA. This is especially
important as it suggests that, despite the relatively high degree of optimisation, the cell-to-cell variability underlied by
(11.34) is biologically relevant.

In the following section we will use this observation as a springboard to study optimal heterogeneous populations.

11.5. Representing optimal populations

Let us start from a rather abstract question. Suppose that an organism is actually maximising a certain function F,
unknown to us, which depends on metabolic fluxes v as well as on a set of other variables w that are not part of
metabolism: F ≡ F(v, w). We shall denote by (v*, w*) the (supposedly unique) configuration of variables where F
attains its maximum. Let's furthermore say that we have a guess for what the organism's objective function might be,
and that this guess is only a function of metabolic fluxes, which we denote by f ≡ f(v). If we trust our guess, and if f
is maximised by the (supposedly unique) flux vector ̂v, our prediction for the fluxes would be ̂v. Question: what is the
probability that ̂v is the true optimum, i.e. that ̂v = v*? Note that f(v*) ≡ f* < ̂f ≡ f(̂v) (i.e. at the 'true' optimum
the value of f is bound to be smaller than the maximum value of f).

The answer goes like this: according to the MaxEnt principle, the probability density p(v) for any flux configuration v to be the true state of the system (i.e. the true optimum) should be undetermined other than by our knowledge that the real optimum has some value of f below ̂f. What is the correct constraint to enforce (besides normalisation) if we are to look for such a p(v)? We could impose that allowed configurations strictly have some fixed value of f < ̂f. This choice would lead to a uniform density over all states with a given value of f. In this way, however, we are assuming that states with a different value of f are inaccessible, a rather strong hypothesis. However, if we impose that only the mean value of f is constrained, MaxEnt will return a probability density with the exact same mean value as the uniform density just described (by construction) but a much larger entropy, just because –intuitively– it will assign a non-zero probability to many more states. Hence, as long as we have no other information, the best prediction we can make for p(v) is given by the probability density that maximises the entropy H[p] subject to the constraint ⟨f⟩ ≡ ∫P p(v)f(v)dν = f*, i.e. by the solution of

$$
\max_{p} - \int_{\mathcal{P}} p(\mathbf{v}) \ln p(\mathbf{v}) d\mathbf{v} \quad \text{subject to} \quad \int_{\mathcal{P}} p(\mathbf{v}) d\mathbf{v} = 1 \quad \text{and} \quad \int_{\mathcal{P}} p(\mathbf{v}) f(\mathbf{v}) d\mathbf{v} = f^* . \qquad (11.37)
$$

We now know the result to be given by (11.10), i.e.

$$
p(\mathbf{v}) = \frac{1}{Z(\beta)} e^{\beta f(\mathbf{v})} \quad (\mathbf{v} \in \mathcal{P}) , \qquad (11.38)
$$

where $\beta$ is the Lagrange multiplier enforcing the constraint $\langle f \rangle = f^*$. What this means in practice is this: if one is
modelling a microbe's metabolism and is unsure about the objective function but has a guess ($f$), information theory
suggests that the best one can do is to assume that metabolic flux configurations are selected according to (11.38).
Ideally, the value of $\beta$ for which one obtains the best agreement between predictions based on sampling (11.38) and
experiments is the 'degree' to which the system optimises $f$. If $f$ is the true objective function, then the agreement
between theory and experiments will get better and better as $\beta$ increases. It is important to keep in mind that (i) while
we have assumed that the organism is actually maximising something, we didn't really use the fact that $F$ is maximised
at $(\mathbf{v}^*, \mathbf{w}^*)$ (only that the true state of the system has some value of $f$ below $\hat{f}$); (ii) this is a totally ideal situation
(for instance, experimental data have errors, so whether comparisons between theory and experiments are informative
doesn't only depend on the theory but also on the quality of the data).

The fact that (11.38) is ‘optimal’ in a rather fundamental sense (a priori different from the sense in which f-maximising populations are optimal) encourages to view distributions described by (11.34) through a different lens. When we maximise the entropy at fixed mean growth rate, in practice, we are looking for the ‘broadest’ probability density (i.e. the most variable population) on $\mathcal{P}$ that is compatible with the given mean. In other terms, we are saying that, given a
yeast *Saccharomyces cerevisiae*, and mammalian, say, human cells. All three exhibit overflow metabolism - even when enough oxygen is available in the environment - yet the underlying biochemistry tells us that *E. coli* ferments glucose into acetic acid, *S. cerevisiae* - into ethanol, and human cells - into lactate. Bringing more contrasts on the table, there might be extreme differences in a single taxon already: some yeasts, for instance, will never produce ethanol when oxygen is present; some of them have lost the ability to do respiration at all over the course of evolution. This might sound like playing a trivia game, but in many cases, meaningful modeling of complex biological systems requires both taking and making biochemical insight. Therefore, when we aim not only to uncover the underlying principles, but also to learn biochemistry, more detailed models such as FBA model have an upper hand.

Yet we already know from previous chapters that the predictions of canonical FBA models are limited to substrate-efficient metabolic states. Continuing with the example of the overflow metabolism, FBA models would predict *E. coli* or *S. cerevisiae* to respire on minimal medium with glucose as the main carbon source – regardless of the maximal flux of glucose into the cell. Thus the prediction of substrate-inefficient metabolism using FBA over the years used to rely on introducing additional, mainly empirical (e.g. maximal oxygen uptake), constraints onto the system [158]. Moreover, we can impose only linear constraints in FBA models, and this greatly reduces our options.

Overall, we often seek to take the advantageous points of both "schools of modeling", however, this is where we need to start doing compromises. In an ideal world, the self-replicator models from Chapter 8 would be much more detailed, and would be extended with explicit kinetics and thermodynamic constraints to obtain a detailed cell model. However, the number of variables would increase tremendously, and non-linear optimization is very inefficient already past even small systems. On the contrary, we could try to advance on existing FBA-type models by introducing the concepts of protein economy (Chapter 7) at genome-scale, as well as self-replication. Following our best understanding, these, again, would constitute non-linear relationships (e.g. enzyme kinetics), yet large-scale non-linear programming is not a viable option either. Thus simplifications are necessary to keep linearity (and convexity) to solve optimization problems for large-scale models.

So can we make large-scale models tractable? If we 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 of metabolites in the growing cell volume, or simply - the growth rate $\mu$ of the cell.

### 10.1.1 Replacing enzyme kinetics by linear catalytic constraints

To obtain large, detailed cell model that we can actually solve, all relationships between models variables have to be linearized. This concerns, most importantly, all catalyzed processes: we assume a linear dependence between a catalyzed flux and its catalyst (enzyme or machine) concentration, but ignore the dependence on the concentrations of substrates, products, 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(s)$: the rate $v$ is proportional to enzyme level $e$ and enzyme efficiency $k$, which is given by a kinetic rate law $k(s)$, a nonlinear function of the metabolite concentrations. Depending on the context, $k$ is also called apparent $k_{\text{cat}}$. The kinetic rate laws $k(s)$ have typical shapes, as described in Chapter 3.

To linearize the expression for $v$, while keeping the dependence on $e$, we need to replace the ratio $k = v/e$ by a fixed number, and so $k$ becomes 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$ [24], which is feasible for a limited number of enzymes, however. Obviously, in reality, neither $v$ 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).
<figure><img src="image_77.png" /><figcaption>Figure 11.5: Fitness-information bound (general form). The black line encodes the maximum mean growth rate achievable for any given value of the information content $I$ (Eq. (11.39)) of a metabolic flux distribution (or the minimum value of $I$ required to achieve any given mean growth rate).</figcaption></figure>

mean growth rate, the optimal population is the one that has the largest possible variability. To quantify variability in a more readily understandable way, it is convenient to transform it into a measure of the amount of information encoded in $p$. One can reason as follows: if no prior information is available about the population, uncertainty is maximal and all flux vectors in $P$ must be considered to be equally likely. This means that, for such a population, the probability density over $P$ is uniform (see (11.7)). We shall denote the entropy of the uniform distribution over $P$ by $H(0)$. When we inject information into the problem (e.g. the fact that the population has a certain mean growth rate), then the probability density is no longer uniform but given, say, by (11.34). The uncertainty is therefore reduced by $H(0) - H(\beta)$, where $H(\beta)$ is the entropy of (11.34). (Clearly, $H(0)$ is just the entropy of (11.34) for $\beta = 0$.) The quantity

$$I = \frac{H(0) - H(\beta)}{\ln 2} \qquad (11.39)$$

denotes the amount of information (in bits, hence the factor $\ln 2$) injected by a non-zero value of $\beta$. Re-formulating our population-level optimisation, we can say that, for any fixed mean growth rate $\langle \lambda \rangle$, the optimal population is the one carrying the smallest value of $I$. A short calculation shows that $\langle \lambda \rangle$ and $I$ are related by

$$\beta \langle \lambda \rangle = I \ln 2 + \int_{0}^{\beta} \langle \lambda \rangle d\beta' , \qquad (11.40)$$

where it should be noted that $\langle \lambda \rangle$ is an increasing function of $\beta$ (as $\beta$ increases, the density gets more and more concentrated around the growth-rate maximising flux vector, thereby leading to an increase of $\langle \lambda \rangle$).

**Problem 11.7** Retrieve formula (11.40).

The curve $\langle \lambda \rangle$ versus $I$ described by (11.40) can therefore be computed numerically for any metabolic network reconstruction (as the only ingredients required are encoded in the flux polytope $P$) [432]. The resulting line (see Figure 11.5) separates the $(\langle \lambda \rangle, I)$ plane in a viable (achievable) region and a forbidden region where the mean growth rates are too large for the amount of information encoded in the population. This ‘phase diagram’ yields, first and foremost, a general prediction linking the mean growth rate (fitness) of a microbial population to its metabolic heterogeneity: all populations must have fitness-heterogeneity values in the viable region. Recent work relying on an advanced statistical inference framework has shown that actual microbial populations indeed lie in the viable part of the plane [435]. In addition, it provides a quantitative definition of an optimal population that accounts for variability: optimal populations have fitness-heterogeneity pairs that lie on the boundary between the viable and the forbidden region. In this respect, results from [432, 434, 435] can be summarised by saying that heterogeneous, faster-growing E. coli populations (mean growth rate larger than roughly 1/h, richer growth media) are very close to optimality, while slower-growing ones tend to have
Under this assumption, we can replace all kinetic constraints by two linear constraints on the enzyme. If we consider coefficients $k$ and $k'$ to approximate enzyme kinetics in the forward and backward direction, respectively, the flux the enzyme $e$ catalyzes should satisfy $-e k' \le v \le e k$. We set $k' = 0$ for irreversible reactions, and, for simplicity reasons, we usually assume $k = k'$ for reversible reactions, unless kinetic measurements are available that suggest otherwise. This relationship can be formulated as enzyme capacity constraints in order to replace the kinetic rate laws in the FBA model. By writing down such constraints for each enzyme in the model, we can couple the metabolic fluxes with the demand for enzymes, needed to operate these fluxes.

By linearizing all formulae as described above, it is possible to build very large models, describing resource allocation on genome scale. What we commonly refer to as "resource allocation models" therefore formalize the mathematical relationships defining the interactions and allocation of resources between the cellular processes to describe optimal resource allocation using constraint-based models. All these relationships take the form of linear, growth-rate dependent equalities and inequalities, and, when linearized, form a convex feasibility problem [265, 266, 267].

### 10.1.2 Overview of existing FBA extensions

By itself, the idea of constraining metabolic models to represent limited metabolic capacity of cells is not new. There are two ways to approach this budgeting problem. The first possibility is "protein budgeting", where a fixed amount of protein needs to be partitioned in the optimal manner (e.g. to maximize growth). The second, alternative method is "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.

Some enzyme-constrained variants of FBA account for empirical constraints on the total concentration of metabolic enzymes (FBA with molecular crowding, or FBAwMC [268]), or on proteome sectors (Constrained-Allocation FBA, or CAFBA [269]). While these types of models 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. Therefore, resource allocation models extend and embed the ideas of proteome partitioning beyond frameworks like CAFBA and GECKO [270], or representing metabolic capacity limitations beyond FBAwMC.

Currently, there are three main implementations of large-scale resource allocation models: Resource Balance Analysis (RBA) [271], Models of Metabolism and Macromolecular Expression (ME-models) [272] and proteome-constrained models (pc-models) [39]. All these implementations are formalized as LP feasibility problems at a fixed growth rate, where the growth rate can then be maximized in an additional optimization loop. 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.

### 10.1.3 Why maximize the growth rate?

Under the assumption of the balanced growth, the copy number of each cell component is doubled between two consecutive cell divisions. If metabolites are described by their concentration, we can think of cell growth as dilution by which the concentrations of all compounds would go down if their amounts remain the same. For a given compounds, dilution by growth can be effectively modeled of every metabolite by a "consuming reaction", with a flux given by $v_{dil} = \mu s$, 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 s$ 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 $s$ exists. Furthermore, the
mean growth rates and information contents that get more and more sub-optimal the less rich is the growth medium. (Of course, this notion of optimality refers to the growth rate and information content as the key parameters to evaluate a population's performance. It may well be, and this is an issue definitely worth exploring, that slower-growing population are optimal with respect to some other parameter(s).) At any rate, the above definition of optimality coincides with the standard one (growth-rate maximisation) for $\beta \to +\infty$, in which case variability goes strictly speaking to zero as all cells collapse on the same flux configuration. And we now understand how it generalises it: by stressing the way in which heterogeneous populations can be optimal despite growing at sub-maximal rates.

For later convenience, note that, because the entropy is a convex functional, the solution to the MaxEnt problem is the same as the solution to

$$ \max_p \int_{\mathcal{P}} p(\mathbf{v})\lambda(\mathbf{v})d\mathbf{v} \quad \text{subject to} \quad \int_{\mathcal{P}} p(\mathbf{v})d\mathbf{v} = 1 \quad \text{and} \quad -\int_{\mathcal{P}} p(\mathbf{v}) \ln p(\mathbf{v})d\mathbf{v} = H^*. $$

The above problem has perhaps a more direct interpretation: the optimal population is the one that has the largest mean growth rate at fixed variability (entropy) or, equivalently, at fixed information content.

Before moving on, we notice that, in the above setting, optimality of heterogeneous populations has a rather simple mechanistic interpretation in terms of how populations 'occupy' the flux polytope. If one considers the uniform distribution on $\mathcal{P}$, Eq. (11.7), and calculates the marginal distribution for the growth rate (i.e. (11.36)), one finds that the growth-rate landscape in which populations grow is extremely skewed towards slow growth rates: the overwhelming majority of metabolic flux configurations corresponds to slow-growing cells, i.e. with growth rates roughly two orders of magnitude below $\lambda_{max}$. This implies that, whatever flux vector we are in, a small random change to it is overwhelmingly more likely to decrease our growth rate than increase it. In this respect, slow states have an 'entropic' advantage over fast states. On the other hand, by definition, fast-growing flux configurations replicate faster than slow-growing ones, and therefore have a replicative advantage. It is therefore tempting to interpret the probability density (11.34) as resulting from the balance between these two tendencies. One can for instance imagine that a microbial population grows and evolves in time in $\mathcal{P}$ due to (i) replication events, and (ii) small random changes of the flux vector (due e.g. to gene expression noise). Ref. [432] has indeed shown that such a population naturally evolves toward a distribution very close to (11.34), where the role of $\beta$ is played by the inverse rate of diffusion of the population in $\mathcal{P}$, that is, by the inverse of the rate at which small random changes occur: fast rate implies small $\beta$, and vice versa. (As the mathematical analysis of this scenario requires the toolbox of non-linear Fokker-Planck equations, it is beyond the scopes of this Chapter.)

The above theory can be extended in various directions. We shall limit ourselves to one example here, namely that of optimal populations in fluctuating environments [436]. The basic assumption we make is that the growth rate $\lambda$ is a function of both the flux vector $\mathbf{v}$ and of a single (for simplicity) exogenous variable $s \ge 0$ representing the stress level to which the population is subject: $\lambda \equiv \lambda(\mathbf{v}, s)$. We furthermore assume that $s$ is a random variable with probability density $P(s)$. For any value of $s$, $\lambda$ will be maximised by a certain flux vector $\mathbf{v}^* \equiv \mathbf{v}^*(s)$. If fluctuations of $s$ are sufficiently slow, then cells may be able to perfectly adapt their metabolic response to every value of $s$ they encounter. But this is unlikely to be possible in rapidly fluctuating environments. In the latter case, it is instead reasonable to assume that cells will try to maximise their average growth rate, where the average is taken over the distribution of $s$. And the relevant distribution to consider is now the probability to observe a certain flux configuration $\mathbf{v}$ given that the state of the environment is $s: p(\mathbf{v}|s)$. So the objective function (to be maximised over $p(\mathbf{v}|s)$) is just

$$ \langle \lambda \rangle = \int ds P(s) \int_{\mathcal{P}} p(\mathbf{v}|s)\lambda(\mathbf{v}, s)d\mathbf{v}. \qquad (11.41) $$

Now we must specify the constraints. One is simple and concerns normalisation: $\int_{\mathcal{P}} p(\mathbf{v}|s)d\mathbf{v}$ should be equal to one for all $s$. To introduce the second one, we note that, because one expects $\mathbf{v}$ to encode information about the environment, it is convenient to constrain the mutual information between $\mathbf{v}$ and $s$, i.e.

$$ I(\mathbf{v}; s) = \int ds P(s) \int_{\mathcal{P}} p(\mathbf{v}|s) \log_2 \frac{p(\mathbf{v}, s)}{p(\mathbf{v})P(s)} d\mathbf{v}, \qquad (11.42) $$