feasible combinations ($\mu$, $v$, $s$) form a convex set, with possible solutions ($v$, $s$) for low values of $\mu$ and no solutions above a critical value $\mu_{max}$, the maximal possible growth rate for our model. Finding this critical value as well as the corresponding optimal fluxes $v$ and compound concentrations $s$ is relatively easy, and can be done by bisection: solving a series of Linear Programming problems (checking for potential solutions ($v$, $s$) for different values of $\mu$).

## 10.2 The basic constraints in resource allocation models

As mentioned above, fine-grained resource allocation models build on genome-scale metabolic models (GEMs) to encompass all the reactions that can be employed in a metabolic network. The technical advance, when constructing such models, is to impose sets of additional constraints onto GEMs to couple the metabolic fluxes with investment into metabolic pathways (production of enzymes). To the date, different implementations of this concept were proposed to predict optimal resource allocation in different microorganisms [273].

Although the precise formulations vary, resource allocation models build on three principal types of constraints (Figure 10.1):

(1) Mass-conservation constraints

(2) Flux coupling constraints

(3) Compartment capacity, or protein density, constraints

The general description of these constraints in fact is the same as for small, coarse-grained self-replicator models, only the number of individual constraints increases. Moreover, every of the constraints described can be split into a number of constraints, considering only a subset of fluxes in the model (e.g., fluxes taking place in a certain cell compartment).

Alongside these three major types of constraints, there is another set of constraints, which we may call "environment" constraints - these correspond to, e.g. the composition of growth medium, biomass composition at at given growth rate $\mu$, etc. They are implemented by setting target values for amounts and/or fluxes defining a viable cell in a given (or several) environmental conditions, but they are not structural constraints. These constraints usually are added *ad hoc* and do not need to bear any functional meaning *per se*. We will now expand on the three types of constraints used in resource allocation models; note that the description is not exhaustive and peculiarities may vary among different formulations.

### 10.2.1 Steady-state and mass-conservation constraints

Mass-conservation constraints define the metabolic network (stoichiometry and relation between fluxes). The initial building blocks of these extended models are GEMs, and thus the metabolic network stoichiometry is already there; what remains to be defined are the protein turnover processes. We consider four types of protein turnover reactions in fine-grained resource allocation models: protein synthesis, folding, degradation and dilution-by-growth. So, for every protein present in such a model, we add these four reactions: two of them, translation and degradation, include the stoichiometry of amino acids needed for its translation and released during degradation based on the protein sequence. The reactions which represent either protein folding modeled as the conversion of the "unfolded" protein species into the "folded" ones, and the dilution-by-growth is modeled as a sink for the "folded" protein species ("folded" $\rightarrow \emptyset$).

### 10.2.2 Catalytic constraints

Next, the flux coupling constraints couple the metabolic fluxes with protein usage: usually, the usage scales with the catalytic turnover value $k_{cat}$ of the enzyme. In this step we have to collect the kinetic information (in most cases, $k_{cat}$ values), which are used as model parameters. We establish the coupling between fluxes and protein synthesis by setting $v = k_{cat} e \eta$, where $e$ is the enzyme concentration and $0 < \eta \le 1$ is an efficiency term summarizing the effects of reaction thermodynamics, enzyme saturation, and possibly small-molecule regulation. The value for $\eta$ can be either assumed or fitted from experimental data, and when $\eta = 1$, the enzyme is considered to operate at its maximal rate. Coupling constraints are introduced to couple both (i) the metabolic reactions with enzyme usage (as described above) and (ii) protein
where $p(\mathbf{v}, s) = P(s)p(\mathbf{v}|s)$ is the joint distribution of $\mathbf{v}$ and $s$, while $p(\mathbf{v}) = \int ds P(s)p(\mathbf{v}|s)$. Clearly, $I = 0$ if $p(\mathbf{v}, s)$ factorises over $\mathbf{v}$ and $s$ and it gets larger and larger as $\mathbf{v}$ and $s$ become more and more correlated. Putting these pieces together, we can write the cell's optimisation problem as

$$ \max_{p(\mathbf{v}|s)} \int ds P(s) \int_{\mathcal{P}} p(\mathbf{v}|s) \lambda(\mathbf{v}, s) d\mathbf{v} \quad \text{subject to} \quad \int_{\mathcal{P}} p(\mathbf{v}|s) d\mathbf{v} = 1 \quad (\forall s) \\ \text{and} \quad \int ds P(s) \int_{\mathcal{P}} p(\mathbf{v}|s) \log_2 \frac{p(\mathbf{v}|s)}{p(\mathbf{v})} d\mathbf{v} = I^* . \quad (11.43) $$

A comparison with (11.5) should clarify how the above generalises the previously discussed optimisation framework. Again using the method of Lagrange multipliers one finds that the optimal probability density is now given by

$$ p(\mathbf{v}|s) = \frac{p(\mathbf{v})}{Z(s, \beta)} e^{\beta \lambda(\mathbf{v}, s)}, \quad (11.44) $$

where

$$ Z(s, \beta) = \int_{\mathcal{P}} d\mathbf{v} p(\mathbf{v}) e^{\beta \lambda(\mathbf{v}, s)}, \quad (11.45) $$

while $\beta$ is a Lagrange multiplier.

**Problem 11.8** Retrieve formula (11.44) (hard).

The meaning of (11.44) is straightforward: when $\beta \to 0$, the metabolic flux configuration $\mathbf{v}$ becomes independent of $s$, implying $I = 0$. As $\beta$ increases, $\mathbf{v}$ and $s$ get more and more correlated, while $p^*(\mathbf{v}|s)$ tends to get more and more sharply peaked around $\mathbf{v}^*(s)$. In the limit $\beta \to +\infty$ cells respond to each value of $s$ by selecting the exact flux configuration that maximises $\lambda$. To achieve this, maximal $I$ is required. A detailed study of the optimal probability density emerging in this case within a highly coarse-grained model of metabolism has been carried out in [436], showing how complex metabolic strategies (including the coexistence of slow-growing, persistent states with fast-growing ones) arise as optimal responses to a fluctuating environment.

## 11.6. Discussion

Metabolic variability in cell populations has, as we have discussed, multiple origins, both rooted in unavoidable stochastic effects and (possibly) in the fact that, in certain cases, being heterogeneous can be optimal for a microbial population. Models can account for variability by representing (sufficiently large) populations via probability densities defined on the flux polytope. Two main (different) goals can be achieved. First, one can look for the probability density that yields the best (in a precise sense) description of a set of empirical data. Methods like Maximum Likelihood and Maximum Entropy provide different, albeit related, approaches to this task. Second, one can formulate optimization problems for populations, whose general solution is a probability density rather than a single flux configuration. Solutions to these problems can highlight how fitness and variability are related in optimal populations, providing useful theoretical benchmarks for real microbial systems. While possibly more demanding from a mathematical viewpoint (and certainly more demanding from a computational viewpoint), these approaches significantly expand the scope of CBMs, including in terms of predictive power. In addition, they can refine the notion of optimality and provide insights into the fundamental principles that govern the organization of metabolism across populations. The question of whether variability confers an advantage to microbial populations is however very general, and goes beyond the metabolic level of CBMs on which we focused here. A broader discussion of these aspects is presented in the Chapter *Cell behavior in the face of uncertainty*.

### Economic analogy 11.F: Maximum Entropy economic equilibrium

Most of economic theory relies on the assumption that markets are capable of allocating resources optimally, i.e. so that the utilities of each of the participating agents is maximised (an assumption can be seen as the analog
<figure><img src="image_106.png" /><figcaption>Figure 10.1: Overview of biological components and mathematical constraints in large-scale resource allocation models – A Resource Balance Analysis (RBA) model is shown as an example. (A) Typically, an RBA model describes metabolisms and macromolecule production in a growing cell (yellow blocks). Precursors from metabolism are needed to produce macromolecules, and some macromolecules serve as enzymes to catalyze metabolic reactions. In addition, macromolecules are diluted and are localized in cell compartments. (B) Sets of mathematical constraints. The variables and processes described by an RBA model must satisfy a number of constraints, include mass-balance constraints (between production, degradation, and dilution of compounds); capacity constraints (relating process velocities to the concentrations of catalysts); density constraints (on the total amount of compounds in a cell compartment); and possibly empirical physiological constraints on any types of "target variables", to ensure realistic models.</figcaption></figure>

turnover reactions with the respective macromolecular machinery (e.g. sum demand of ribosomes for protein translation, $v_{\text{translation}} = [\text{Ribosome}] \times k_{\text{cat,ribosome}}$). The sheer number of the kinetic parameters needed for formulating the coupling constraints in the fine-grained models requires the modeler to consider different assumptions and simplifications when building and parameterizing these models, as briefly discussed below.

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

The next issue is kinetic parametrization of these fine-grained models. We currently can use only very simplified kinetics in the models (flux coupling $v = k_{\text{cat}} e^{-\eta t}$), and simplify such factors as enzyme saturation and thermodynamic driving force into a single value of factor $\eta$. Two approaches are used to deal with
of each cell in a population maximising its growth rate). In order to achieve optimal states (called 'equilibria' in economics), agents endowed with a *priori* different preferences, resources and goals identify the actions that maximise their utilities (e.g. transactions, trade, or production) and carry them out. This process however can become more and more demanding as the number of agents that take part in the market gets larger and larger, because (in short) the set of viable transactions for each agent can become exceedingly large. How can one describe the equilibria that arise from these situations?

A possible approach, used at least since [437] (see also [438, 439, 440]), is based on the Maximum Entropy principle. The idea, in short, is the following. Once every agent has somehow chosen their preferred actions (i.e. once a system-wide 'configuration of individual actions' has been selected), the market as a whole presents a set of transactions to be carried out that aggregate the choices of individual agents. When looked at the aggregate level, though, each set of transactions can correspond to more than one configuration of individual actions. (This can happen, for instance, because agents have a degree overlap in their characteristics which makes them indistinguishable from an economic perspective.) If one assumes that agents choose their actions at random from their set of viable transactions, then some sets of transactions are bound to be more likely than others, simply because they can be realised in more 'microscopic' ways (for instance, by interchanging agents of the same type). It is then reasonable to think that the likelihood of any particular set of transactions will be larger, the larger the number of microscopic ways in which it can be realised. Taking entropy as a measure of multiplicity, the most likely set of transactions, then, is the one that maximises the entropy.

A model of market where the above program is worked out in detail is found in [437]. The 'statistical equilibrium' theory that follows from the use of the Maximum Entropy principle generalises the standard competitive equilibrium discussed in microeconomics by providing a description of optimality in large markets with heterogeneous participants. This line of work has also inspired further developments that explicitly included agents' heterogeneity into the theory of competitive equilibria [441, 442, 443]. To the best of our knowledge, a similar approach has not yet been used to model heterogeneous microbial systems.
Box 10.A Protein abundance versus concentration in resource allocation models

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

The predicted amount of protein in cells can be compared to experimental measurements in two ways. First option is to convert abundance to concentration using the relationship between the cell volume and dry weight (e.g. $V_{gDW} = 1.7 \text{ mL gDW}^{-1}$ in *Saccharomyces cerevisiae*, [17]). Alternatively, proteome mass fractions are a popular unit in label-free mass spectrometry-based protein quantification, a popular method in quantitative microbiology. Respectively, predicted proteome mass fractions can be inferred by converting protein abundance in *mmol* to *g*, and scaling to the protein content in dry cell biomass. Here, it is important to consider the conversion factors (protein content in dry biomass). *E. coli* maintains rather constant protein content in dry weight across growth rates (ca. 0.55 (*g* protein) *gDW*⁻¹) [275, 276]. On the contrary, the protein content is known to vary in *S. cerevisiae* as a function of growth rate [17].

this, as a large fraction of parameters are not even available. First, condition-dependent kinetic param-
eters ("apparent catalytic constants", k App ) are fitted from experimental (mostly quantitative proteomics)
data (setting k eff = k cat α, where 0 < α ≤ 1) with a value α chosen to match predicted enzyme abundance
and experimental measurements. Otherwise, for the enzymes with measured k cat values, we can assume
that enzymes work at their maximal rate, i.e. the saturation function η = 1. Then the model computes
the minimal protein requirement to sustain the flux through the metabolic reactions. The comparison of
minimal predicted vs. observed protein abundance can represent the "apparent saturation", or "overcapac-
ity" of enzymes. For instance, it is common in yeast S. cerevisiae that the flux and not protein expression
varies across conditions, and the relationship between predicted and measured expression can suggest
the nature of the observed protein expression [274].

10.2.3 Protein density constraints

The final layer of information in our resource allocation models is a set of protein density constraints in
each cell. These constraints put an upper limit on the amounts machines driving the cellular processes,
e.g. a maximal protein capacity of a compartment. These constraints are formulated as weighted sums
of protein abundance, usually with weights proportional to the proteins' molecular weight. Usually, the
density constraints are expressed in terms of (usually maximal) mass, area, and volume of the compartment
(e.g. "what is the maximal mass the mitochondrial proteins can take up in gDW of cells?"). Based on the
biological interpretation of the constraints, we formulate the weighing multipliers to represent either of
the metrics (mass/area/volume) that every protein occupies.

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

10.2.4 Interpreting the consequences of the additional constraints

We have briefly discussed what types of additional constraints need to be implemented to extend FBA
models to account for cellular resource allocation, and now let us recap on what these sets of rules mean
in biological terms. The constraints described above shall couple the metabolic fluxes with the production
of enzymes that operate these functions, so the model has to produce amino acids and generate ATP in
order to use them for protein translation. Moreover, the enzyme demand will be coupled with the produc-
tion of the macromolecular machines required to produce, fold, and degrade these enzymes (ribosomes,
Chapter 12

Control of cell division and coordination
with other cell-cycle processes

Mattia Corigliano, Marco Cosentino Lagomarsino, Jacopo Grilli, and Gabriele Micali

Chapter highlights

* Cells require coordination of growth and division, as well as coordination of cell-cycle progression with several essential sub-tasks, such as chromosome replication and segregation.

* Single-cell dynamics data offer correlation patterns that can be used to understand these decisional processes.

* The cell-cycle progression and cell-division decisional process can be described by continuous-time and discrete-time stochastic processes.

* There are quantitative relationships that connect growth, cell-cycle progression, and resource allocation.

* There are differences and common points in the decisional processes by which single cells of different organisms commit to divide (sizers, adders, accumulators, titration-dilutors, etc.)

12.1. Introduction: the decision to divide illustrated through single-cell *E. coli* data.

As nicely put by the Nobel prize winner François Jacob, "the dream of every cell is to become two cells". Achieving this dream often requires multiple steps, such as growing by a certain size, replicating DNA, and dividing. The previous chapters have addressed cell growth as a consequence of optimization of catabolic and biosynthetic fluxes through optimally regulated resource allocation; this chapter deals with the decision to divide (and to progress the cell cycle), based on growth and other important cellular processes and cues. Clearly this decision to divide or progress the cell cycle must be based on a set of inputs (growth, production processes such as DNA replication and cell-wall biosynthesis, partitioning processes, etc.) and entails several outputs, prominently cell division, but also intermediate key cell-cycle substeps, such as initiation of DNA replication or construction of a "divisome" organelle. The questions that we will consider concern the characterization of the known aspects of this decisional process and its coupling to cell size, to cell growth, and to the chromosome cycle. We will use throughout the chapter *E. coli* as an example. This section provides a description of the main problem through an introduction to the data, based on *E. coli* bacteria. Sections 2-5 start from a mathematical toolbox of models that are useful in this context and compare them with data. Finally, section 6 describes applications to other organisms than *E. coli*.

Capturing the key processes regulating cell division is a fundamental question in biology, which remains open despite a history of more than 60 years. During the years, scientists have learned a great deal about the size and shape of bacteria in different nutrient conditions, what most of the molecular players involved in the division process are, how the DNA
chaperones, and proteases, respectively), requiring the same building blocks (see Chapter 2). These con-
straints therefore formalize a self-replicating molecular system in balanced growth subject to different
structural constraints:

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

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

3. The intracellular density of compartments and the occupancy of membranes must not exceed the de-fined limits.

As highlighted before, the biological interpretation of the additional constraints discussed above is rather
universal for different implementations of resource allocation models, with minor deviations in terminol-
ogy and/or formulation. To illustrate how resource allocation models are built from conventional GEMs,
and how the respective models are formalized in mathematical terms, in the following we will consider
one of the popular formulations of resource allocation models in more depth.

**10.3 Resource Balance Analysis: model construction and simulation**

Resource Balance Analysis (RBA) is a flexible and generic modeling framework that describes the function-
ing of an organism using linear equality and inequality constraints, as described in general terms in Section
10.2. As a consequence, an RBA model includes all known metabolic reactions coupled to relevant cell
processes with major protein investments (production of biomass precursors; including, but not limited to
protein translation, protein folding, protein transmembrane transport, and protein degradation). Where
applicable, circumstantial information can be included into the model to establish the dependency of en-
zyme activity on metal ions, vitamins, and/or cofactors. Which metabolic reactions and cell processes are
regarded as relevant may vary between organisms and is a modeler's choice.

**10.3.1 Building a draft RBA model**

The software package RBApy [277] contains all the routines needed to build and simulate RBA models. In order to build a new RBA model, it takes as an input a genome-scale metabolic network in SBML format [278], together with additional information to formulate all the constraints described above. Different types of biological data are needed to build an RBA model for a given type of cell:

○ Amino acid sequences for metabolic enzymes and macromolecular machines (e.g. ribosomes and chap-
erones),

○ If applicable, stoichiometry of known cofactors (e.g. metal ions),

○ Efficiencies of metabolic enzymes,

O Molecular weights and localization of proteins (for density constraints),

○ Any empirical constraints on concentrations or fluxes ("targets", see previous section).

The software first extends the input GEM with a description of protein production and dilution in the cell.
To do so, it extracts information from the input files on (i) protein sequences and cofactors, (ii) the subunit
stoichiometry of protein complexes, and (iii) protein localization (using information from public databases
such as UniProt). Using this information, reactions corresponding to protein synthesis, folding, degrada-
tion, and dilution by growth are added automatically. Finally, the software maps enzymes to the reactions
they catalyze and to the proteins they consist of. The output of the routine is a draft (uncalibrated) RBA
model.

**10.3.2 Mathematical description of a RBA problem**

**Notation.** Below $A^T$ refers to the transpose of the matrix $A$. $\mathbb{R}_{>0}^n \triangleq \{x \in \mathbb{R}^n | x_i > 0 \text{ for all } i \in \{1, \dots, n\}\}$, $\mathbb{R}_{>0} \triangleq \mathbb{R}_{>0}^1$, $\mathbb{R}_{\ge 0}^n \triangleq \{x \in \mathbb{R}^n | x_i \ge 0 \text{ for all } i \in \{1, \dots, n\}\}$ and $\mathbb{R}_{\ge 0} \triangleq \mathbb{R}_{\ge 0}^1$.
<figure><img src="image_78.png" /><figcaption>Figure 12.1: Salient quantitative features of cell-division control, explained through *E. coli* data. A: *E. coli* cells are rod-like. Within a condition they grow by increasing their length, and they divide symmetrically. Following single-cell lineages, growth in length or volume is close to exponential. B: Size-growth plots quantify the strength of division control. For a timer, multiplicative growth quantified by $G = \log(s_d/s_0)$ is uncoupled to birth size, for a sizer, it is maximally coupled. The single-cell data show an intermediate trend. C: Since $G = \log(s_d/s_0) = \alpha\tau$, the size-growth plot can be split into contributions correlating birth size to growth rate (top) and/or interdivision time. The data show that *E. coli* bacteria only compensate by modulating interdivision times. D: Equivalent quantifications of the strength of the division control by relating birth and division size (left) or added size and birth size (right). The intermediate control strategy adopted by *E. coli* adds a size that is independent from the initial size ("adder"). This strategy is sufficient to achieve size homeostasis.</figcaption></figure>

replication machinery is formed and how it proceeds along the chromosome, how the septum and the new cell wall are synthesized. However, the vast majority of these data are based on population averages, out of which it turns out to be impossible to extract any direct and/or causal link between the different processes involved in cell growth that set cell division [444]. Today, a new generation of data has the potential to answer several open questions [444, 445, 446]. These data differ from the previous generation in the ability to measure single bacterial cells over multiple division events in controlled conditions. At the same time, the expression of a specific gene or the concentration of specific proteins of interest can be monitored using fluorescent reporters. For example, fluorescent tags on the proteins involved in replication are used to score the initiation of replication in each cell cycle. Single-cell data allow for validating mathematical models and thus bring insights into the causal link between the several processes a cell need to complete before dividing.

By following lineages of cells over multiple generations under controlled environmental conditions, scientists collected different important pieces of evidence (Figure 12.1): First, within a cell cycle, the cell size $s(t)$ is well-described by a single exponential in time<fnref n="1" /> [449, 450]: $s(t) = s_0 \exp(\alpha t)$, where $s_0$ is the size at birth, $\alpha$ is the growth rate, and $t$ is the time since cell birth.

If division occurs at time $\tau_d$, a simple relationship connects the size at division $s_d$ with the other cell properties: $s_d = s_0 \exp(\alpha\tau_d)$. All the four parameters of this equation are subject to stochasticity in time and vary across single cells, even when they grow in controlled conditions. Second, in steady growth, the size distribution of newborn cells does not change over time, an observation that is referred to as cell-size homeostasis [447]. Equivalently, cells show specific correlation patterns between size at growth and size at division, which are related to their cell-division strategy [451, 447].

Let us try to understand more in detail how single-cell correlation patterns can be used to understand cell-division behaviors. The observation of near-exponential growth immediately suggests a change of variables that is useful to formulate mathematical models and to understand how single cells control cell division. Indeed, if we can assume that growth is exponential, we can use logarithmic sizes instead of linear sizes. One robust observation, is that the elongation $G = \log(s_d/s_0) = \alpha\tau$ depends on the size at birth $s_0$ (Figure 12.1B). This allows us to generate so-called "size-growth"

<footnote n="1">¹Note that most of the studies today use cell length as a proxy for size. However, different choices are possible such as volume or mass, and the differences are not fully characterized [447, 448].</footnote>
In a standard RBA model, we consider balanced growth (see Chapter ??), that is, the average state of a cell in a cell bacterial population growing exponentially at the specific (constant) growth rate $\mu \ge 0$, i.e. the amount of produced biomass per biomass per cell per unit of time. Our simulated average cell is composed of different molecule species:

1. $n_y$ types of molecular machines, which can be subdivided further into $n_e$ enzymes and transporters involved in the metabolic network $\mathbb{E} \triangleq (E_1, \dots, E_{n_e})$ at the concentrations $\mathbf{e} \triangleq (e_1, \dots, e_{n_e})^T$ and metabolic fluxes $\mathbf{v} \triangleq (v_1, \dots, v_{n_e})^T$; and $n_m$ macromolecular machines $\mathbb{M} \triangleq (M_1, \dots, M_{n_m})$ involved in non-metabolic cellular processes, such as the translation apparatus, at the concentrations $\mathbf{m} \triangleq (m_1, \dots, m_{n_m})^T$;

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

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

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

For a given cell growth rate $\mu \ge 0$, the RBA optimization problem (named $\mathcal{P}_{\text{rba}}(\mu)$) can be formalized mathematically as follows. For a fixed vector of concentrations $\mathbf{p} \in \mathbb{R}_{>0}^{N_{n_p}}$ and the given growth rate $\mu \ge 0$,

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

where all the inequalities are defined component-wise and:

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

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

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

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

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

○ $\bar{\mathbf{d}}$ is a vector of size $n_c$, where $n_c$ is the number of compartments (compartment membrane and/or compartment interior for which density constraints are considered. $\bar{\mathbf{d}}^i$ is the density of molecular entities within the volume or surface area. Densities are typically expressed as a number of amino-acid
plots (Figure 12.1B), in which the log-multiplicative growth during a cell cycle of a single cell is plotted as a function of the logarithmic size at birth [451]. Different mechanisms of size control predict different slopes for this plot. A cell division set by a “timer”, for instance, would predict no relation between $G$ and size. Since $G = \log s_d - \log s_0$, if instead $\log s_d$ were independent of the initial size, a “sizer”, one would predict a slope $= -1$. The E. coli data typically fall half way in between these two predictions, a negative slope of about 0.5 (Figure 12.1B).

By noticing that the overall logarithmic growth $G$ during a cell cycle is the product of the single-cell growth rate and inter-division time ($G = \alpha\tau$), we can ask the question of which one of these variables is responsible for the correlation. This analysis disentangles the contributions to cell division control due to growth rate and inter-division timing (Figure 12.1C). In other words, the dependency of $G$ on initial size can be further decomposed on the dependency of growth rate $\alpha$ and division time $\tau$. In E. coli, when growth rate and interdivision times are plotted separately as a function of the logarithmic size at birth, the negative slope is only observed in the interdivision-time plot, suggesting that cell control size by adjusting the single-cell interdivision time rather than their single-cell growth rates. Hence, E. coli data indicate that $\tau$ does depend strongly on initial size, while the growth rate has only a weak dependency [450].

One can visualize and quantify the mutual dependencies between cell sizes and growth properties in other equivalent ways (Figure 12.1D). For example, in E. coli data, the scatter plot relating size at division in the y-axis to size at birth in the x-axis for single cells has a slope of around 1 (and once again this observation holds true for different strains and under different environmental conditions). In this plot, a slope of 0 would suggest that cells on average need to reach a threshold in size upon division, a sizer. More technically, the division size $s_d$ is independent on the initial size $s_0$ in the case of a sizer. Instead, a slope of 2 in this plot would suggest that cells on average need to wait a fixed time upon division, a timer. The observed intermediate slope of 1 can also be understood using the equivalent plot in which the added size between birth and division is used on the y-axis, studying the dependency of the added size $s_d - s_0$ on $s_0$. This latter way to plot the data is particularly popular, given that, for many datasets it shows no dependency, suggesting that adding a constant added size is the mechanism of size control effectively in place. Indeed, for E. coli the experimentally observed slope is always close to 0 [452, 453, 447], an observation that goes under the name of “adder” behavior since cells appear to add on average a constant size during the cell cycle (Figure 12.1BD).

It is fairly simple to rationalize why, for exponentially growing cells, a cell division strategy based on a timer does not achieve a homeostatic size. In order to do this, we can call $q(i) = \log(s_0(i))$ the logarithmic cell size at birth of cell-cycle $i$, and look at its dynamics through subsequent cell cycles. Since $s(\tau) = s_0 \exp(\alpha\tau)$, and $\langle \alpha\tau \rangle = \log 2$, and assuming that cells divide perfectly in two halves, one immediately gets that

$$q(i+1) - q(i) = \nu(i)$$

where $\nu(i)$ is a zero-average random variable independent for each cell-cycle, arising from the size-independent fluctuations of inter-division times (hence, in technical jargon, we can model $\nu$ as a discrete-time Markovian random process). Since the jumps in logarithmic size between subsequent cell cycles are random and independent, cell size at birth makes a discrete-time multiplicative random walk, hence, within a population, the distribution of cell sizes at birth tends to get wider and wider across divisions. The following two sections will explain how size homeostasis can be achieved by size-coupled cell divisions.

**Problem 12.1** Show that for cells that grow linearly in time an adder and a timer are the same.

## 12.2. Hazard rate approach to cell division

As we have seen in the previous section, E. coli cells grow roughly exponentially. Hence, we can describe their growth by a trajectory for size $s$ (measured as cell mass or volume) of the kind $s(t) = s_0 \exp(\alpha t)$, where $t$ is time from cell birth. While experimentally the growth rate $\alpha$ fluctuates with time, we will neglect its variability and assume for the moment that it is constant. As a consequence, the cell grows as a simple exponential function of time. We will address different
residues by volume or surface area.

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

For given growth rate and medium composition, all equalities and inequalities in our RBA problem $P_{rba}(\mu)$ is linear in the decision variables $(y, v)$ and is proven to be convex [265, 267]. At given $\mu$, $P_{rba}(\mu)$ is a feasibility optimization problem, where constraints $(C_1-C_3)$ define the feasibility domain. The feasibility domain can be empty or non-empty. If there exists a solution $(y, v)$ to $P_{rba}(\mu)$ -i.e. the feasibility domain is non-empty-, then there exists a feasible resource distribution compatible with the given growth rate. In other words, the cell can grow at this growth rate value. By construction, the feasibility domain of $P_{rba}(\mu)$ corresponds to the set of all possible phenotypes of the cell at a growth rate $\mu \ge 0$.

We conclude this with some remarks:

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

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

3. In [266, 279], an RBA model was built for *Bacillus subtilis*. It integrates two macromolecular processes in constraint $C_{2a}$, the translation and chaperoning of proteins, and two density constraints, the limitation of the cytosolic density and of the membrane occupancy. An RBA model can be refined by integrating for instance other cellular processes and molecular machines, such as the transcription machinery, the protein secretion apparatus (see [279, 277]), or molecule turnover [280], as well as other types of constraints.

### 10.3.3 Simulation and analysis of RBA models

**How to incorporate the medium composition.** We represent the medium composition in two aspects, namely (i) qualitatively, by allowing exchange of the medium metabolites in the model ($UB_{Exchange,n} > 0$). Note that some metabolites, although not explicitly represented by the growth media, should also adhere to this rule (e.g. oxygen, water, and protons). The (ii) quantitative composition of the growth medium is determined by extracellular concentrations, which, in turn, dictate the efficiencies of metabolic transporters via Michaelis-Menten-like rate laws (as nonlinear $k(c)$ functions; see section 10.1.3). For an extracellular nutrient $S_{ext,i}$ with concentration $s_{ext,i} \ge 0$, the efficiency of the corresponding metabolic transporter(s) is given by $k_E(s_{ext,i}) = \frac{k_{cat}s_{ext,i}}{K_M+s_{ext,i}}$, with parameters $k_{cat}$ and $K_M$ for the turnover number and the affinity of the transporter, respectively.

**Obtaining the RBA solution for a given parameter set.** For an RBA problem with given parameters, there exists a maximal growth rate $\mu^* \ge 0$, such that for any $\mu$, $P_{rba}(\mu)$ is feasible if and only if $\mu \le \mu^*$ [265, 267]. For a given medium composition, the maximal growth rate $\mu^*$ can be computed by using a bisection algorithm, in which a series of LP problems are solved to narrow down the exact growth rate at which the problem becomes infeasible. A real-life example would be simulating growth in glucose-limited chemostat cultures under different dilution rates $D$. With increasing $D$, the glucose availability increases, and a set of $n$ different glucose uptake rates $q_{Glc}(q_{Glc,1}, q_{Glc,2}, \dots, q_{Glc,n})$ can be subjected to an RBA model to obtain a set of optimal metabolic states $(\mu_1^*, \mu_2^*, \dots, \mu_n^*)$.

Together with the maximal feasible growth rate one obtains the optimal cell configuration maximizing growth $(\mu^*, y^*, v^*)$. The principle of optimal performance, in this case, that a cell phenotype should maximize growth rate, in fact, coincides with the principle of parsimonious resource allocation between cellular processes.

**Exploration of the feasibility domain.** Although RBA models inherently reduce the solution space due to principle of parsimonious resource allocation, the solutions obtained might still contain considerable flux variability. In the same vein as Flux Variability Analysis ([281], see Chapter 5), the feasibility domain can be
<figure><img src="image_79.png" /><figcaption>Figure 12.2: Illustration of the inverse hazard rate approach on data. Data from many lineages of dividing cells can be used to estimate the cumulative distribution of non-divided cells, which can also be conditioned on different variables. The drawn example refers to the case where the tested variable is the added size $s - s_0$. In this case, the formalism allows to extract mathematically the hazard rate $h_d(s - s_0)$ from this distribution. Experimental E. coli data are consistent with this adder scenario, with an hazard rate that peaks at a characteristic added size, after which the division control weakens.</figcaption></figure>

hypotheses regarding this point in the later sections.

A simple way to describe the decision processes leading to division (or other cell cycle progression events) is the so-called “hazard rate” model [454, 450, 453]. In this framework, as the cell cycle progresses, each cell has a certain probability to divide, and we call $h_d$ the rate of cell division. In principle, this rate can be a function of many different internal cellular parameters, all the processes that contribute setting cell division. However, since we have in mind experiments measuring cell size versus time and recording cell divisions, the most general “empirically accessible” $h_d$ can depend on $s, t, s_0, \alpha$ with the constraint that $s/s_0 = \exp(\alpha t)$. This means that there can be at most three free parameters. We can also consider simplified models, such as $h_d = h_d(s)$ or $h_d = h_d(s, t)$. Empirically, the lack of correlation between $\alpha$ and birth size suggests a smaller role for this parameter. It is important to realize that this formalism is very powerful, as it can be applied more widely to any sub-cell cycle decision (for example, entry into a specific phase, such as initiation of DNA replication, mitosis, etc.), and to measurements of different relevant cell-cycle processes (for example chromosome configurations or the expression of cell-cycle proteins or other factors), which the hazard rate may depend on.

Given a model for the hazard rate, we are interested in the cumulative probability $F(t|s_0, \alpha)$ that a cell born at $t=0$ has not divided at time $t$, given that its initial size is $s_0$ and its exponential growth rate $\alpha$. The Mathematical Detail Box 12.A discusses the mathematical formalism to obtain this probability.

The considerations we made so far are sufficient to produce “forward models” where a hazard rate is assumed, and one explores its consequences on the division dynamics. The simulation of such a model is straightforward. For each discretized time increment $dt$, the cell will grow by the prescribed dynamics $s(t)$ and will divide with hazard rate $h_d$. If a division occurs, the mother’s cell size will halve, and go from $s$ to $s/2$ (we assume for simplicity perfect binary divisions, but this assumption can easily be relaxed). What is a “sizer” in this framework? We can define it as a model where $h_d = h_d(s)$ [455]. Equally, a timer is a model where $h_d = h_d(t)$, and an adder has $h_d = h_d(s - s_0)$. At this stage, it is only intuitive, but not formally grounded, that the scatter plots of the previous section correspond precisely to these models. This problem will be discussed in section 12.3. Note that not all the choices of hazard rates will guarantee a steady-state cell size distribution. As a particular case, one can consider a constant division rate $h_d(t) = r$, which is a simple Poisson process (see the problem above). This is a pure timer and we expect that it will not maintain a steady-state cell size distribution (the reader can verify it, e.g. by simulations).

Beyond the forward approach, we would like to recognize the trends in the data that favor one model rather than another. In particular, we can ask which model best describes the E. coli data, presented in the first section of this chapter. This question is a “reverse problem”, and is equivalent to the inference of the hazard rate $h_d$ from data (Figure 12.2). It is a very common reverse problem for the literature, used for example in the so-called “survival analysis” in clinical studies [456]. In that case, the hazard rate typically corresponds to a one-time negative outcome (death of the patient) and the process is not repeated along lineages as in the case of cell divisions. However, the mathematical ingredients are very similar. Consequently, there are many regression methods available in the literature, which can be transferred to our case. One of the most famous is Cox regression [457]. However, most of these regression methods need an ansatz
explored at optimal ($\mu^*$) or sub-optimal ($\mu \le \mu^*$) growth rates. For one decision variable $y_i$ (resp. $v_i$), two LP problems are solved, where (i) constraints $C_1, C_2$ and $C_3$ remain unchanged; (ii) the decision variable $y_i$ (resp. $v_i$) is maximized (LP 1) and minimized (LP 2). This operation is repeated for each decision variable to obtain *in fine* the feasibility domain of all decision variables.

It was proven that the feasibility domain becomes smaller with increasing growth rate [265, 267], so it might be worthwhile to probe the solution space at slow-growth regimes. In practice, at the optimum, the cell configuration ($\mu^*, \mathbf{y}^*, \mathbf{v}^*$) is often unique. Indeed, non-unique solutions will exist if two alternative metabolic pathways have exactly the same cost in resources. Since all enzymes have different amino acid sequences, use different cofactors, are differently localization, etc, this is highly unlikely. A caricatural example of a model with non-unique solutions would be one in which an enzyme pool is arbitrarily split into two, and the two new "enzyme species" are given different names, although they are physically exactly the same.

### 10.3.4 Calibration of model parameters

An RBA model may contain a high number of model parameters. First, the global parameters to be estimated are related to cell composition: (i) the concentrations of bulk biomass components $\bar{b}$, which is usually deduced from the biomass reaction of the genome-scale metabolic network of the organism. Using quantitative proteomics data [282], one can infer (ii) the protein densities in different compartments ($\bar{d}$), and (iii) the abundance of housekeeping (unspecified) proteins ($\mathbf{p}$).

The next set of parameters we need to collect concerns the efficiencies of molecular machines ($K_E$, $K_{E1}$, $K_{T1}$). As we learned in Chapters 2 and 3, the rate of an enzymatic reaction $v$ depends on the enzyme efficiency or "apparent catalytic rate", given by $v = ek_{app}$, with $k_{app} = \kappa(\mathbf{c}) = k_{cat}^+ \cdot \eta^{\text{rev}}(\mathbf{c}) \cdot \eta^{\text{sat}}(\mathbf{c}) < k_{cat}^+$. The $k_{app}$ values are always below the $k_{cat}$ value, but may vary from state to state depending on metabolite concentrations. Since internal metabolite concentrations $s$ are unknown and difficult to measure at genome-scale, we cannot estimate $k_{app}$ from the explicit kinetic law $\kappa(\mathbf{c})$. We need to obtain these $k_{app}$ parameters empirically, for example by measuring the flux $v$ and the protein abundance $e$ in one condition and taking their ratio.

Hence, for a given environmental condition, efficiency parameters can be estimated using quantitative proteomics in combination with fluxomics [279] or FBA to estimate the flux distribution [277]. To account for variable enzyme efficiencies, one may make the simplifying assumption that enzyme efficiencies depend mostly on growth rate. By estimating the enzyme efficiencies at different growth rates and interpolating between them, one obtains empirical relationships between efficiency and the growth rate [279] to be used in $P_{rba}(\mu)$. For instance, several estimates of enzymatic efficiencies obtained in contrasting growth conditions will provide a relationship $K_E(\mu)$ instead of a constant $K_E$ value.

### 10.3.5 Enzyme efficiencies: use of -omics data-informed $k_{app}$ vs. naïve $k_{cat}$ values

The three most popular formalisms of fine-grained resource allocation models, RBA [271], ME-models [272], and pc-models [39], are variations on the same theme, implementing the four major constraints discussed in Section 10.2. Thus most of the ideas, concepts, and constraints are equivalent (or at least highly similar) in their biological interpretation. Most of the differences arise from the approach taken towards parametrization of these models, and consequently, interpretation of model output. Here we will discuss an example where implementations differ significantly.

In resource allocation models, two types of constraints define the proteome capacity at given growth rate $\mu$, the protein density vector $\bar{b}$, and the fraction of housekeeping proteins $\mathbf{p}$ in the proteome. The remaining proteome space is to be distributed among the proteins that are explicitly defined in the model. The RBA formalism requires to formulate the function $k_{app}(\mu)$ (or $K_E(\mu)$ in the RBA problem, Section 10.3.2) for every protein in the model using -omics data (see Section 10.3.3), and the fraction of the "housekeeping" proteins in the proteome is determined from data for each simulation.

Conversely, the formulation of pc-models [39] allows more flexibility to the "unspecified" protein *UP*, represented by a single artificial protein of average size and amino acid composition. Instead of setting a
Mathematical details 12.A: Probability distribution of (un)divided cells

This Mathematical Detail Box derives the probability distribution of (un)divided cells from the hazard rate. The probability that a cell divides between time $t$ and $t + dt$ is the probability of not having divided so far times the probability of dividing between $t$ and $t + dt$, in turn given by the product of the hazard rate and the time interval $dt$, $F(t|s_0, \alpha)h_d(s(t), t, s_0, \alpha)dt$. During the same time interval, the cumulative probability of not having divided will decrease by the same amount. Hence, we can write

$$F(t + dt|s_0, \alpha) = F(t|s_0, \alpha)[1 - h_d(s(t), t, s_0, \alpha)dt]. \quad (12.1)$$

In the limit of $dt \to 0$ we obtain a differential equation, which governs the evolution of our system

$$\frac{d}{dt}F(t|s_0, \alpha) = -h_d(s(t), t, s_0, \alpha)F(t|s_0, \alpha), \quad (12.2)$$

and whose formal solution is (for $t \ge 0$)

$$F(t|s_0, \alpha) = e^{-\int_0^t dz h_d(s(z), z)}. \quad (12.3)$$

Since we said that the probability of a cell division event in the time interval $[t, t + dt]$ is $P(t|s_0, \alpha)dt = F(t|s_0, \alpha)h_d dt$, the corresponding probability density is

$$P(t|s_0, \alpha) = h_d(s, t)e^{-\int_0^t dz h_d(s(z), z)} = -\frac{d}{dt}F(t|s_0, \alpha). \quad (12.4)$$

Alternatively, the size $s$ can be used as a coordinate, considering for $s > s_0$,

$$F(s|s_0, \alpha) = e^{-\int_{s_0}^{s} dz h_d^*(z, t(z))}, \quad (12.5)$$

while $F(s|s_0, \alpha) = 0$ for $s < s_0$. Here, $h_d^*(s, t(s))dx$ is the probability of cell division in the size range between $s$ and $s + ds$. The two rates are simply related by $h_d^*(s, t(s))ds = h_d(s(t), s)dt$, where $ds/dt = h_q(s) = \alpha s$ is the rate of growth.

**Problem 12.2** Analyze the consequences of a constant per-size hazard rate $h_d^* = 1/\tilde{s}$ and compare them to the consequence of a constant per-time $h_d = r$ (a Poisson process).

**Problem 12.3** Analyze the forward hazard rate model for cell division where $h_d(s) = (s/\tilde{s}^2)$ by simulation and/or analytical calculations.

for the parameterization of the model, which might be a nuisance, as it would require some previous knowledge. Here we consider a simpler, direct inference method, which does not need any parameterization (but is effective only with a sufficient amount of data, i.e., for many cell divisions).

Suppose for simplicity we deal with a sizer. In this case, it is possible generate an estimator for the functional form of $h_d(s)$ using Eq. (12.5). By inversion, we obtain

$$h_d(s) = -\alpha s \frac{d}{ds} \log[F(s|s_0)], \quad (12.6)$$

where $F$ can easily be estimated from data, from the cumulative fraction of undivided cells at size $s$ with initial size $s_0$. In our case, we can use the mean value of the growth rate $\langle \alpha \rangle$, since we are neglecting fluctuations in growth rate.

Since we do not know whether our assumption of a sizer applies to data, we can first combine the data and the inference to falsify the assumption [450]. In order to do this, we can further condition our histograms in order to fix $s_0$. If $h_d$
fixed amount allocated to **p** which changes across conditions, one can determine the minimal fraction of this protein in proteome $UP_{min}$, and formulate the demand to produce $UP$ as an inequality constraint $UP \ge UP_{min}$. Interestingly, in *Saccharomyces cerevisiae*, the proteome mass fraction occupied by non-metabolic proteins is relatively constant under different glucose-limited conditions, as determined by quantitative proteomics data (see [39], Fig. S1 for a plot).

This inequality constraint can be interpreted as the upper limit of available protein space, i.e., under fixed protein density $y + p = \text{const.}$, the proteome not occupied by $y \triangleq e + m$ is allocated to $p$. Since now the model can distribute the proteome among explicitly-defined vs. unspecified protein freely, the procedure of fitting $k_{\text{app}}$ values is no longer a prerequisite. Using $k_{\text{cat}}$ values, collected from literature/databases/own experimental measurements, rather than apparent $k_{\text{app}}$ values, has consequences both for predictions and the data use: first, the model prediction on the protein use is the "demand" of the enzyme and is strictly coupled to the flux through the enzyme (equivalent to the ECM1 layer of enzyme costs in the *Enzyme Cost Minimization* method, Chapter 6). Second, the condition-dependent quantitative proteomics data can be used as validation dataset for model predictions instead [283], as the predicted protein abundance is not dependent on these datasets.

Using less data for parameter fitting and redirecting these data-rich datasets towards validation of model prediction strengthens the argument for using resource allocation models for learning new biology, and already has real-life examples. For instance, the discrepancies in predicted vs. observed levels of glycolytic enzymes at glucose-scarce conditions in [39] inspired the same team to revisit the question whether the high levels of glycolytic enzymes represent the optimal expression given very low thermodynamic driving force and undersaturation of glycolytic enzymes. Comparing predictions of *Enzyme Cost Minimization* models with the results of the pc-model and experimental data, [274] proposed that *S. cerevisiae* expresses genuine excess of glycolytic enzymes in glucose-limited conditions, meant to amply consume any glucose as soon as it appears in the environment.

## 10.4 Biomass composition as a constraint or as a prediction

Cell models describe, among other things, what a cell is composed of (see Chapter 2). In FBA, specifically, "biomass" refers to the proportions of different molecule classes (e.g. lipids, protein, DNA, RNA, cofactors) in 1 gram dry weight of cells, and biomass composition needs to be defined prior to optimization. Since, at least for FBA models of microbes, biomass production usually is the optimization objective, the literature frequently refers to the mathematical description of cell composition as "biomass objective function" (BOF). In most cases, it is assumed that the proportions of biomass constituents are fixed, only the total production (flux through BOF) changes.

For the predictions of FBA models to be reliable, a realistic BOF is a must (see Chapter 5). Therefore, there is a sustained effort to experimentally determine biomass composition, even for *E. coli* [284]; for more details on the usual experimental measurement methods, see the box in Chapter 2. In case supporting data are available, the cell composition in the BOF may be described in a more fine-grained manner for individual molecule types (e.g., individual lipids, proteins, mRNA species, etc), or even in terms of atomic composition (which in turn gives clues about the amounts of molecule classes). So, overall, biomass composition acts as a global, and one of the most stringent, constraint on the predicted solution space in FBA-based models.

However, cell composition may greatly vary not only between (micro-)organisms, or different cell types within the same organism, but also for the same organism/cell type across different conditions. The budding yeast *S. cerevisiae*, for instance, exhibits rather linear relationships between the proportions of bulk biomass constituents as a function of growth rate in glucose-limited cultures [3]. This variable composition often poses a challenge for models: just like the uptake rates, the varying biomass composition reflects complex global rearrangements of resources (for instance, different ribosome content at different growth rates [35] leads to changes in RNA-to-protein ratio in the cells), and choices between metabolic strategies (e.g. depletion of storage carbohydrates in glucose-fermenting *S. cerevisiae* [17]).

A main advance of resource allocation models, compared to conventional FBA models, is that only a part of the biomass composition is given as input information just like in FBA ($\bar{b}$ in RBA). The proteome com-
depends solely on $s$, then the inferred function $\tilde{h}_d$ should not change with varying $s_0$. This is indeed the case if the procedure is applied to simulated data. However, when we apply the same procedure to the experimental data shown in the previous section, the inferred $h_d(s)$ changes if it is inferred for different bins of birth size $s_0$. Hence, we conclude that our E. coli data do not behave as a sizer, in the sense of the hazard rate. Instead, if we consider the adder ansatz for the hazard rate $h_d(s - s_0)$, and we repeat the procedure, we find that further conditioning by birth size or time from birth does not change our inferred hazard rate [453]. Hence, we can conclude that a hazard-rate analysis of the data supports an adder (or at least that the data cannot falsify this simple model).

How does the inferred $h_d$ depend on size? Curiously, for any fixed $s_0$, $h_d$ increases superlinearly for small cell sizes, then reaches a maximum after which it decreases. In other words, some cells may “miss” a cell division event and keep growing until they find a better occasion to divide. This process is called “filamentation” (because the cells that miss one or more division elongate and end up looking like filaments), and is typically the consequence of stress, but also present in stress-free growth conditions. experimental observations show that E. coli forms filaments in response to DNA damage, antibiotics, host immune systems, temperature, starvation, and many other stresses. As a consequence, size plasticity may be in many cases an adaptive strategy. The quantitative division rules of filamentous E. coli cells have been studied experimentally [458], but we lack a comprehensive mathematical model.

One very robust observation of cell division statistics, in E. coli and beyond [459, 453, 460, 461], is that the distributions of size at birth, size at division, and division times measured across conditions, collapse onto the same curve when rescaled by their mean. For instance, the distributions around these values are clearly non-overlapping: the single-cell birth-size distribution in glucose $p_{glu}(s_0)$ strongly differ from the one in TSB medium $p_{TSB}(s_0)$. In particular, the typical size at birth for E. coli growing in glucose $\langle x_0 \rangle_{glu}$ is about 2/3 the size of E. coli growing in TSB $\langle s_0 \rangle_{TSB}$ and the average division time $\langle \tau_d \rangle_{TSB}$ is TSB is half the one of E. coli in glucose $\langle \tau_d \rangle_{glu}$. This appears to be valid across different environmental conditions (e.g., nutrient quality, temperature, pH, etc.). The remarkable empirical observation is that, when comparing two conditions, the rescaled distribution is universal. If we introduce the rescaled size $\tilde{s}_0 = s_0/\langle s_0 \rangle_c$, the distribution of $\tilde{s}_0$ is universal, independent of the condition. This observation applies also to size at division, added size between divisions, interdivision time, and, to a certain extent, growth rate [460].

An obvious question that follows from this observation is how the size-scaling properties of cell-size at birth constrain the mechanisms of homeostasis and the properties of stochasticity at the single-cell level. A necessary consequence of the distribution collapse is that the processes leading to single-cell heterogeneity and homeostasis must have common underlying properties across conditions. Conditions differ because they are characterized by different dimensional scales, but, phenomenologically, division control is governed by the same underlying principles (although the key molecular players may vary). The collapse of all the distributions, when the variables are rescaled by the mean has another, stronger, consequence: whatever the division control mechanism is, it depends on only two scales, a size-scale (setting the typical cell size) and a temporal scale (setting growth rate and division time).

This constraint has strict consequences on the variability of the hazard rate across conditions. In particular, it implies that the hazard rate must take the mathematical form [462]

$$h_d(s(t), s_0, t\alpha) = \alpha \tilde{h} \left( \frac{s(t)}{\langle s \rangle_c}, \frac{s_0}{\langle s \rangle_c} \right), \quad (12.7)$$

where the function $\tilde{h}(\cdot, \cdot)$ is the same across conditions. The dependency on $\alpha$ and $t$ disappears, as the scaling of division time, implies the existence of a unique time scale. Since $\tilde{h}(\cdot, \cdot)$ is by definition adimensional, it can only depend on the product $\alpha t$, which can always be re-expressed as a function of $s$ and $s_0$, as $\alpha t = \log(s(t)/s_0)$. While this is a powerful observation, as it allows to naturally connect division mechanisms across conditions, it does not provide any evidence to a particular decisional mechanism enforcing cell division, which is encoded in the function $\tilde{h}(\cdot, \cdot)$. Addressing this question needs further experimental details.
position, on the contrary, becomes a genuine prediction of the optimization procedure. Unlike small self-replicator models (see the models in Chapter 8), this prediction is very detailed, as the predicted proteome composition is represented by the sum of individual protein abundances. Moreover, if proteins require trace elements or cofactors (e.g. iron in iron-containing proteins) for function, the demand and contribution to the overall biomass of these metabolites will also be predicted by the model (as it will vary with the expression level of those proteins).

In theory, the abundance of biomass constituents other than proteome could be formulated in the way they become predictions of the resource allocation models, rather than hardcoded inputs. Following the idea implemented in the small, coarse-grained models of [264], one could set relationships between, e.g., protein density in the cells and production of lipids (in [264], the biological interpretation was to maintain the surface area-to-volume ratio constant). Currently this is not widely accepted as a standard practice, and, as we can see from the example above, requires comprehensive experimental evidence, which, by itself, could be interpreted still as "input to the model".

## 10.5 Concluding remarks

In this chapter we considered whole-cell resource allocation models that couple metabolic networks with a description of the macromolecular machinery that is required to operate them. Compared to FBA models, these models contain a large number of additional reactions, metabolites, constraints, and model parameters, and, overall, offer a fine-grained representation of cellular economy. Many of the kinetic parameters cannot be accurately measured for individual enzymes, and/or are condition-dependent. The quantitative nature of the predictions of resource allocation models (and the most cellular decisions/phenotype shifts), however, is largely due to global constraints: for instance, when the protein density $g gDW^{-1}$ in a compartment reaches its upper limit (=that compartment is fully packed with protein), the cells switch from fully-respiratory to respiro-fermentative growth (see [272] for E. coli, or [39] for S. cerevisiae). Unlike the kinetic parameters in single reactions, which are rather uncertain, these "global", cell-wide constraints are based on more trustworthy evidence.

Thus these models still retain a reasonable compromise concerning numerical tractability and model complexity, and can accurately predict complex adaptations, which cannot be captured by GEMs in an autonomous way, i.e. without the addition of empirical constraints on fluxes. A successful use case of using resource allocation models is dissecting iron economy, using RBA models: some proteins require iron for their function, and the cell growth can become iron-limited in some conditions. The RBA model was used to predict cell behavior under iron starvation, and the predictions suggested couple of scenarios, (i) the cell may *increase* the import of iron, but also (ii) *avoid* using proteins that contain iron (and the pathways in which they operate) [279, 285].

As with the biomass composition, another aspect of resource allocation models (and FBA-based models in general) with some duality in its interpretation is the objective function. Although its validity has been always debated since conception, maximization of instantaneous growth rate as the optimization objective has shown incredible success in predicting microbial physiology. The current approach we apply for resource allocation models still remains the FBA-based assumption that the desired cell phenotypes are the ones maximizing instantaneous growth rate $\mu$. This time, however, the $\mu$ is also a model variable, so we have to apply bisection to obtain the optimal solution for each parameter set we use in resource allocation models.

It is becoming more and more evident that many cell phenotypes (and microbial species!), which we try to predict, do not actually maximize instantaneous growth rate. For instance, most experimental research on microbial physiology has been focused on carbon-limited (C-limited) cultures, especially the yeast work in Delft, the Netherlands (see [3, 17] for examples). It seems that the principle of growth rate maximization works very well in C-limited case, and the success of resource allocation models to quantitatively capture these phenotypes [272, 279, 39] affirms this assumption. But is C-limitation descriptive of natural environments? Let us continue the argument with yeasts as an example.

Yeasts in the wild, for instance, are often subjected to feast-famine cycles in terms of carbon availability,
12.3. Cell-division control as discrete-time linear response process

In the previous section, we have seen how the cell-division control mechanism can be mathematically defined using the hazard-rate framework. This approach uses as a fundamental ingredient the probability per unit time of cell division $h_d$, which is, a-priori, a function of many internal cellular parameters. This approach is, in some sense, very general, as it allows to characterize any complex cellular decision process. However, this generality limits the tractability and interpretability of the model. In this section, we introduce an alternative discrete-time mathematical framework which greatly simplifies the parameterization and the interpretation of a cell-division control model [463, 464], and easily maps to the empirical parameters discussed in Figure 12.1.

Specifically, instead of tracking the division rate at different stages of the cell-cycle, it is often convenient to model directly the cell size at birth across different generations. In this case we can, in full generality, write

$$s_0^{i+1} = f(s_0^i, \alpha, \dots) + \eta^i(s_0^i, \alpha, \dots). \qquad (12.8)$$

where $s_0^i$ is the birth size of the cell at generation $i$. The function $\eta(\cdot)$ represents a random variable whose mean is equal to 0 and having, a priori, arbitrary probability distribution. The function $f(\cdot)$ described the control over cell division. Specifically, the function $f(\cdot)$ can be simply (almost tautologically) defined as the conditional average of the size at birth at generation $i+1$ given all the variables that contribute to cell division control (the previous size at birth, the growth rate, and others),

$$f(s_0^i, \alpha, \dots) := \langle s_0^{i+1} \rangle_{s_0^i, \alpha, \dots} \qquad (12.9)$$

The random variable $\eta(s_0^i, \alpha, \dots)$ characterizes the fluctuations around this conditionally averaged birth size.

This formulation of the process is as general as the hazard-rate formalism as it allows to express any division probability $F(s|s_0, \alpha, \dots)$. Eq. (12.8) simply isolates the contribution of the (conditional) average size at division from the deviations from this average. This separation is useful because it allows a clear interpretation of the mechanism of division control, and because the conditional average size at division is typically accessible from single-cell experiments. For instance, a timer corresponds to $f(s_0^i, \alpha) \propto s_0^i$, where the proportionality constant equals $\exp(\alpha\tau_d)/2$. A sizer corresponds to $f(\cdot)$ being a constant, independent of the initial size $s_0^i$. Along the same lines, an adder is defined as $f(s_0^i, \alpha) = (s_0^i + \Delta(\alpha))/2$, where $\Delta(\alpha)$ corresponds to the (average) added size. The formalism also shows how there is a continuum of possible intermediate behaviors besides these three limit cases.

Given the facts that growth is exponential, and the distribution of sizes at birth is approximately Lognormal [453, 460], it is once again convenient to introduce the logarithmic size $q_0^i = \log s_0^i$. One can derive the dynamics of the variable $q_0^i$ as a function of the dynamics defined in Eq. (12.8) [462]. This procedure is described in Mathematical Detail Box 12.B. Since the fluctuations of this variable are small, this dynamics is fully specified by a set of linear-response parameters $\lambda_{ab}$ relating the main observables (i.e. in our case each of the variables $a, b$ can be $q_0, \alpha, \tau, G$).

Mathematical details 12.B

This Mathematical Detail Box derives the dynamic equations of the birth size $q_0^i$ across generations (indexed by $i$) in the discrete-time formalism. We define $\langle q_0 \rangle_\alpha$ as the average value of $q_0$, and the log size deviation $\delta q_0^i := q_0^i - \langle q_0 \rangle_\alpha$. The dynamics for the log-size deviation takes the form

$$\delta q_0^{i+1} = g(\delta q_0^i, \alpha) + \zeta^i(\delta q_0^i, \alpha), \qquad (12.10)$$

where $\zeta^i(\delta q_0^i, \alpha)$ is a random variable with zero mean. This equation has the same degree of generality of Eq. (12.8) and can express any arbitrary division control model (or equivalently any shapes of the hazard rate function). In order to make further mathematical (and biological) progress, we need to simplify the equation and make the comparison with data possible. There are several possible choices. In the following, for simplicity, we first neglect the fluctuation of the growth rate $\alpha$. Assume that the size at birth is the only variable influencing cell division
and one could argue that in the famine phase of the cycle, these yeasts should act as if they were glucose-limited. Yet the current opinion in the yeast ecology seems to see feast-famine cycles as a continuous, although reduced, supply of carbon, and steer towards embracing a higher role of nitrogen (N) limitation in natural environments instead. Currently, our understanding of N-limited growth is not very comprehensive, and N-limitation is also a case where the instantaneous growth rate maximization breaks down: the pc-models of *S. cerevisiae* cannot quantitatively capture the cell behavior under N-limited conditions.

So the selection of a suitable optimization objective can be a choice followed by huge success, but also, the optimal solution might end up contradicting the existing knowledge. How can we try to mitigate that? One huge advance of resource allocation models is that at any condition, the available solution space is greatly reduced, compared to conventional FBA. We can argue that we have introduced a whole new set, a whole new type of constraints into the model by accepting assumptions stemming from the metabolism-molecular machinery coupling. In theory, we should be able to reason further regarding any additional (even empirical/ad hoc) constraints and/or additional optimization objectives which would bring our model predictions closer to observed biology. Just remember: fitting models is not a sin; but nontransparent/reckless fitting is! After all, modeling is an art, and there is no one cookbook that represents the ground truth: we should be free to explore the secrets of biology, as unrealistic as our assumptions are at times.

A final remark on modeling being an art. In this book, we have explored several types of cell models of different size, detail, and assumptions behind. This whole hierarchy and diversity of different implementations and formalisms might seem overcomplicated and unnecessary, but it is a mere reflection that "one size does not fit all". Hence, and we invite (future) modelers to be creative, mix, match, and tailor different models (and modeling types) to cover more and more biological knowledge. The compromise between fine-grained but linear modeling vs. complex kinetics that materialized into resource allocation models is an inspiring example of how one can push bounds of different methods to create mathematical representation of our mental pictures of cells.

## Recommended readings

**RBA website** Website rba.inrae.fr for further details on RBA. Under Tools, you can find example models and Jupyter notebooks for running them.

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

## Problems

### Problem 10.1 The role of metabolite concentrations

The available cell space for proteins depends on the assumed space occupied by small metabolites.

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

(b) In what way would a cell, in reality, profit from a lower small metabolite content? Can we assume that the ratio between small metabolites and proteins is optimized? Describe possible aspects of this compromise! For inspiration, see [286].
$(g(\cdot)$ is a function of $\delta q_0^i$ only) will allow us to introduce a linear-response framework. We will then describe how to consider the heterogeneity of multiple growth parameters.

The main empirical observation that comes to our help is the fact that the coefficient of variation of $q_0^i$ is small (typically around 0.15) [445, 453, 463, 460, 462]. The small value of the coefficient of variation strongly suggests the possibility of Taylor-expanding the function $g(\delta q_0^i)$ around $\delta q_0^i = 0$ [463]. In this limit, the function $g(\delta q_0^i)$ is approximately linear and the random variable $\zeta^i(\delta q_0^i, \alpha)$ can be well approximated by a Gaussian random variable with zero mean and constant variance [462]. The resulting equation reads

$$ \delta q_0^{i+1} = (1 - \lambda)\delta q_0^i + \sigma\xi^i, \qquad (12.11) $$

where $\xi^i$ is a Gaussian random variable with zero mean and unit variance. The two parameters $\lambda$ and $\sigma$ encode, respectively, the relevant information about the mechanism of size control and the level of stochasticity. The parameter $\sigma$ simply corresponds to $\zeta^i(0, \alpha)$. The parameter $\lambda$, which quantifies the strength of size control, has a direct relationship with the mechanism at its origin. It is defined as $\lambda = 1 - g'(0, \alpha)$. For instance, the sizer corresponds to $\lambda = 1$ and an adder to $\lambda = 1/2$. The case $\lambda = 0$ does not lead to a stationary process and corresponds to a timer. Consequently, this parameter can easily be inferred from the plots in Figure 12.1.

Eq. (12.11) can be solved analytically [462]. In particular one can show that the conditional probability of observing a log-size deviation $\delta q_0^i$ from the average at generation $i$ given a deviation at generation 0, is a Gaussian with mean

$$ \langle \delta q_0^i \rangle \delta q_0^0 = (1 - \lambda)^i \delta q_0^0. \qquad (12.12) $$

This result clearly shows how different mechanisms correspond to different strengths of cell-size homeostasis, leading to fluctuations persisting across a different number of generations. For a sizer, $\lambda = 1$, the expected deviation of the daughter cell is independent of the mother cell fluctuations. A timer, with $\lambda = 0$, does not lead to homeostasis, as the expected deviation of size at birth of a daughter cell is the same as the deviation of the mother. The adder, $\lambda = 1/2$, leads on average to a halving of the size at birth deviation at each generation, as approximately observed in experiments [453].

One can generalize the linear-response framework to consider fluctuations of different growth parameters [464]. In general, one can assume that the size at birth of the daughter cell depends on both size at birth of the mother and her individual growth rate fluctuations.

$$ \delta q_0^{i+1} = (1 - \lambda_{qq})\delta q_0^i - \lambda_{q\alpha}\delta\alpha^i + \xi_q^i. \qquad (12.13) $$

Along the same lines, one can assume that the growth rate fluctuations obey a similar equation

$$ \delta\alpha^{i+1} = -\lambda_{q\alpha}\delta q_0^i - \lambda_{\alpha\alpha}\delta\alpha^i + \xi_{\alpha}^i. \qquad (12.14) $$

This kind of equation can be written in multiple forms, i.e. including multiple variables. For example, one can write an equation explicitly for the elongation rate between divisions $\delta G := \delta q_0^{i+1} - \delta q_0^i$ or for the division time. Since the linear-response equations assume that the fluctuations around the means of these variables are small, all these choices turn out to be mathematically equivalent. This is also the reason why the different plots in Figure 12.1 are equivalent. While a linear dependency of growth rate $\alpha$ and division time $\tau_d$ on (log-)size at birth $q_0$ would induce a non linear dependency of the elongation $G = \alpha\tau_d$ on the initial size, such non-linearities can be neglected in the limit of small fluctuations, leading always to linear dependencies [463, 464].

The values of the parameters $\lambda_{ab}$ can be easily inferred using the standard tools of linear regression. Notably, the best (maximum likelihood) estimates of these parameters can be directly obtained from the variable covariances [462, 464]. For instance, $\langle \delta q^{i+1} \delta q^i \rangle = \lambda_{qq}\sigma_q^2 + \lambda_{q\alpha}\langle \delta\alpha^i \delta q^i \rangle$. By writing the expressions for other correlations (e.g., $\langle \delta q^{i+1} \delta\alpha^i \rangle$ or $\langle \delta\alpha^{i+1} \delta q^{i+1} \rangle$) one can map the coefficient $\lambda_{ab}$ with the measured covariances.
**Problem 12.4** Find the hazard rate corresponding to the process defined by Eq. (12.11).

**Problem 12.5** Write an explicit expression of the four parameters $\lambda_{ab}$ appearing in Eq. (12.14) and Eq. (12.13) as a function of the covariances between the fluctuations of growth rates and log-size at the same or different generations.

**Problem 12.6** Prove that the adder strategy rapidly achieves cell size homeostasis (that is, a controlled cell size at birth) after a few cell generations, independently of the starting initial size. Prove that convergence to homeostasis and loss of memory of the initial cell size is exponential in the number of cell cycles. Write down a simple numerical code to simulate this process and verify your analytical predictions. What is the role of noise in setting the inter-division added size?

The linear-response framework offers a flexible and analytically tractable tool to formulate and explore different models of division control. The models can be constrained by correlation patterns measured in data, quantified for example by covariances, which relate to the coupling parameters $\lambda_{ab}$. However, the question remains of whether such models are consistent with data. For *E. coli* data, the linear-response framework predicts the correct consistency relations between experimental measurements, thereby confirming its usefulness to characterize empirical data [464]. A second, more biologically relevant, question is identifying the biological mechanism reproducing the observed dependency patterns. As already discussed, the observation that $\lambda_{qq} \sim 0.5$ is a strong indication of adder-like size-control mechanisms [463, 453, 462, 464]. Interestingly, one can show that the non-zero correlation between growth rate and log-initial size $\langle\delta\alpha^{i+1}\delta q_0^i\rangle$ can be explained because of the correlation between mother and daughter single cell growth rates (the presence of a non-zero value $\lambda_{\alpha\alpha}$ and a dependency of the division size on the growth rate (a non-zero term $\lambda_{q\alpha}$). Such a relation between parameters point to some dependency on the size at division on the single cell-growth rate. For *E. coli*, it is possible [464] to reproduce the empirical values of these coupling parameters by assuming an adder model where the added size depends exponentially on the single-cell growth rate, following the same dependency it has on the population growth rate (this behavior will be discussed in more detail below, and is sometimes termed Schaechter's Law [465]).

12.4. Coordination of cell division with different cell-cycle processes

In the previous sections, we learned that *E. coli* single-cell dynamic data reveal the adder size-control behavior, which allows bacterial cells to maintain size homeostasis. We also discussed a mathematical framework that describes how size control is achieved, and, in particular, how the key measured variables (logarithmic size at birth, interdivision time, growth rate, and total growth during a cell cycle) are connected. Here, we introduce a joint description of the DNA replication cycle, which at the modeling level makes it necessary to partition the cell cycle into sub-periods. We then present the key elements and observations around the debate on whether and how DNA replication and genome segregation is limiting cell division in *E. coli*. In presenting this debate we aim to (i) highlight the positive and innovative aspects of some of the cornerstone studies of recent years, (ii) provide the reader with robust tools necessary to compare mathematical models against data. Finally, we conclude the section by underlying a few open questions.

It is a classic question in biology [467, 468] how cells achieve the precise coordination of the cell cycle with chromosome replication and segregation is necessary for cell survival. DNA replication defines a way to subdivide the cell cycle into sub-periods. In *E. coli*, the period between cell division and initiation of DNA replication is normally referred as the B-period. The C-period is the period needed to complete replication. Bacterial DNA is organized in circular chromosomes which replicate starting from a well-defined “origin” region (called *ori* locus). The replication machinery moves bi-directionally, and the two “replication forks” proceed approximately at the same speed and terminate in a “terminus” region of the chromosome called *ter* locus [469, 470, 471]. For *E. coli* cells dividing at mean interdivision times from about 20 minutes to about one hour, the replication speed is approximately constant, resulting in an approximately constant C period of around 40 minutes [472]. The D-period is the period that lasts from the end of replication to the next division which thus includes segregation and septum formation. Note that the inter-division time, i.e. the time between two consecutive division events, can be as short as 20 minutes in *E. coli*. How can a cell with a division
Chapter 11

Optimal cell behavior in time

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

Chapter overview

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

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

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

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

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

11.1 Introduction

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

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

In many situations, the precise functioning of the molecular mechanisms regulating the adaptation of mi-
crobial physiology to changes in the environment is not or only qualitatively understood. This precludes
their inclusion in quantitative models that accurately predict the dynamic response of the cell in a variety
<figure><img src="image_80.png" /><figcaption>Figure 12.3: Comparison of different cell-cycle models including chromosome sub-periods proposed in the literature for *E. coli*. (A) The DNA replication-segregation cycle divides of the cell cycle into sub-periods. The B-period is the period between cell birth and initiation of DNA replication; the C-period is the period needed for completing DNA replication; and the D-period is the period between the termination of DNA replication and cell division. Finally, the I-period is the period between two consecutive initiations of DNA replication, which usually spans two generations. (B) Scheme of the 'replication-centric' class of models in which DNA replication-segregation sets division (first column). These models usually assume that the CD and the I periods are adders (blue lines in the third and fourth column, respectively), in agreement with data (red lines in the same panels). The G-period correlation pattern is a prediction of the model in general agreement with data (yellow vs red lines in the second column). (C) Schematic for the 'replication-agnostic' class of models in which a process starting at cell birth drive division (first column). These models assume the G and I periods to be adders (blue lines in the second and fourth panels, respectively). The C+D period correlation pattern is a prediction of this model which does not agree with the available data (yellow vs red lines in the third panel). (D) Schematic for the 'concurrent cycles' class of models in which two processes compete to set division through an AND gate (first column). These models assume the I periods to be an adder (blue lines in the fourth column) and using additional parameters predict both adders in the G and C+D periods (yellow lines in the second and third column). (E) Plotting the slope of the G versus the C+D-period allows to compare the different models with data. Schematic similar Figure 4 in [466].</figcaption></figure>

time shorter than the C-period duration have at least two copies of the DNA? Classical studies have shown that *E. coli* and other bacteria can set up multiple overlapping rounds of replication, as summarized by Cooper and Helmstetter in 1968 [468]. For example, a cell at birth is already replicating DNA and has two forks. During the cell cycle, two new initiation events take place, which will only terminate in the daughter cells [473]. We will refer to the "G-period" and the "I-period" as the periods between two consequent division and initiation events, respectively.

As briefly mentioned in the introduction of this chapter, the recent single-cell experiments allow to score initiation and termination of DNA replication by fluorescently tagging proteins involved in the formation of the replication forks or directly the *ori* locus [472, 474, 476, 477, 479]. The scoring of initiation and termination makes it possible to produce the size-growth, and the equivalent adder, plots for any of the sub-periods BCD<fnref n="2" /> as well as for the G- and I-periods (jointly). In the remainder of this section, we will refer to the slope of the size-growth plot of a sub-period X (X= B, C, D, G, or I) as $\lambda_X$, and to the slope of the corresponding adder plot as $\zeta_X$. The two slopes are linked by the equation $(1 - \lambda_X) = \frac{\zeta_X+1}{Q_X}$, where $Q_X = \exp(\langle \text{growth during X} \rangle)$ (see Mathematical Detail Box 12.C).

<footnote n="2">²Note that under fast-growing conditions the termination is experimentally harder to score reliably and hence in many studies the C and D periods of single cells are considered together as a "CD period".</footnote>
of conditions. The lack of mechanistic, quantitative information can be bypassed by making appropriate assumptions about the regulatory systems, in particular that the latter have evolved under the selection pressure of the environment to optimize the response to external perturbations. More precisely, it is assumed that microorganisms have developed mechanisms that allocate limiting resources (proteins, fluxes, ...) to cellular processes so as to maximize or minimize some objective function, or combination of objective functions, over the time-interval of environment changes.

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

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

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

For each of the biological problems and corresponding models considered, we give the precise definition of the modeling formalism and the optimization problem, a small example as an illustration, a discussion of the solution of the problem, and a brief description of more realistic applications and the insights they have given into the functioning of cellular networks. The chapter does not give a detailed explanation of the mathematical methods that are used for solving different classes of optimal control problems, because this would require knowledge of specialized mathematical concepts with which the average reader of the book may not be familiar. Moreover, these methods have been the subject of dedicated textbooks [294, 295]. Rather, we focus on the definition of the dynamic optimization problems and the interpretation of the
Having formally defined sub-periods for the cell cycle and the corresponding linear-response formalism, we now proceed by discussing a schematic overview of the experimental observations in E. coli that any mathematical model should reproduce:

* The G-period shows an adder behavior, ($\lambda_G = -0.5$, $\zeta_G = 0$) [452, 453].

* The C-period duration is approximately constant across cells and experimental conditions with, a tendency to increase for slow growth rates and the C-period generally shows a timer behavior<fnref n="3" /> ($\lambda_C = 0$, $\zeta_C = Q_C - 1$) [478, 479, 480, 474, 470].

* The I-period shows an adder behavior, ($\lambda_G = -0.5$, $\zeta_G = 0$) [481, 477, 476].

* The CD-period shows an adder behavior ($\lambda_{CD} = \frac{Q_{CD}-1}{Q_{CD}}$, $\zeta_{CD} = 0$) [477, 482].

* The single-cell growth rate and the duration of the CD period are inversely proportional [475].

Other interesting observations that are considered in the mathematical models we will present shortly are

* E. coli cells divide symmetrically with a narrow distribution of division length with CV = 0.05 [450]. Note that this CV is lower than the CV of both the growth-rate distribution (CV ≈ 0.1) and interdivision time distribution (CV ≈ 0.2).

* The growth rate of the mother cell is correlated positively with the growth rate of the daughter cells, with a Pearson correlation of around 0.5 [449].

Mathematical details 12.C

This Mathematical Detail Box shows how to translate the linear response ("$\lambda$-formalism") to an equivalent formalism based on the slopes of adder plots ("$\zeta$-formalism"). The interested reader can find more information in [463, 464, 481, 483, 466]. As discussed previously, Eq. (12.11) makes it possible to estimate the linear-response parameter $\lambda$ in experimental data from the covariance of log-size fluctuations between subsequent generations, by noticing that $(1 - \lambda_G) = \frac{\langle \delta q_0^{i+1} \delta q_0^i \rangle}{\sigma_{q_0}^2}$, where we refer to $\lambda$ in Eq. (12.11) as $\lambda_G$, to highlight the fact that this equation refers to the G-period. Exponential growth dictates that $2s_0^{i+1} = s_0^i e^{\alpha^i \tau^i}$, where $s_0^i$, $\alpha^i$, and $\tau^i$ are the size at birth, the growth rate and the interdivision time, respectively. For the cell cycle $i$ one can expand the logarithmic growth $G_G^i := \alpha^i \tau^i$ around its average value ($\langle G_G \rangle \approx \log 2$) in terms of variations around the logarithmic size at birth $q_0^i := \log s_0^i$. Following this procedure, the cell size at birth of generation $i + 1$ within a lineage can be expressed as a function of the parameters of generation $i$, as follows,

$$2s_0^{i+1} = Q_G (s_0^i)^{1-\lambda_G} \langle s_0 \rangle^{\lambda_G} + \nu_0^i, \quad (12.15)$$

where $Q_G = e^{\langle G_G \rangle} = \exp \langle \log s_d/s_0 \rangle$, $s_d$ is the cell size at division and $\nu_0^i$ is a discrete-time Gaussian noise with mean zero and standard deviation $\sigma_{s_0}$. Expanding around the average size, for small fluctuations we obtain a mapping between added size and slope of the size-growth plot,

$$2s_0^{i+1} = Q_G \langle s_0 \rangle + (1 - \lambda_G) Q_G \delta s_0^i + \nu_0^i \\ \delta \Delta_G^i = +[(1 - \lambda_G) Q_G - 1] \delta s_0^i + \nu_0^i.$$

Here $\Delta_G^i = s_f^i - s_0^i$ is the added size during a cell cycle, and $\delta \Delta_G^i = \Delta_G^i - \langle \Delta_G^i \rangle$ is its fluctuation. Hence, by definition, the term in square brackets must be the slope of the adder plot

$$\zeta_G := (1 - \lambda_G) Q_G - 1. \quad (12.16)$$

Solving the equation for $\lambda_G$, we get

$$(1 - \lambda_G) = \frac{(\zeta_G + 1)}{Q_G}, \quad (12.17)$$

<footnote n="3">³Given the difficulty in observing the C-period in single cells, this last question requires further experimental investigation.</footnote>
solutions returned by available numerical solvers of optimal control problems.

## 11.2 Mathematical form of dynamic optimization problems

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

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

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

such that

$$ \begin{aligned} \frac{dx}{dt} &= f(x(t), u(t)), \quad x(0) = x_0, \\ 0 &\ge c_1(x(t), u(t)), \\ 0 &\ge c_2(x(0), x(T)). \end{aligned} \qquad (11.2) $$

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

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

$$ \frac{dx}{dt} = N v(x(t), u(t)) - \mu(t) x(t), \quad x(0) = x_0, \qquad (11.3) $$

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

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

Many methods for solving optimal control problems (11.1)-(11.2) exist. While some optimal control problems can be solved analytically, most of the problems considered in the examples below require numerical approximations to be solved. All examples developed in the sections below have been solved by means of freely available solvers.
which can be used (assuming as usual small fluctuations) to convert the slope $\zeta_G$ of the adder plot into the slope of the size-growth plot $\lambda_G$, and vice-versa.

**Problem 12.7** Write the equivalent of Eq. (12.11) for the I-period and for sub-periods B and CD, and prove the following relationships:

$$ (1 - \lambda_I) = \frac{\langle \delta q_I^{i+1} \delta q_I^i \rangle}{\sigma_{q_I}^2}, \quad (1 - \lambda_B) = \frac{\langle \delta q_I^i \delta q_0^i \rangle}{\sigma_{q_0}^2}, \quad (1 - \lambda_{CD}) = \frac{\langle \delta q_0^{i+1} \delta q_I^i \rangle}{\sigma_{q_I}^2}, $$

where the log-size fluctuation at initiation for the cell cycle i is $\delta q_I^i := q_I^i - \langle q_I \rangle \approx \log(s_I^i / \langle s_I \rangle)$, with $s_I^i$ the cell size at initiation.

**Problem 12.8** Write the equivalent of Eq. (12.17) for the I-period and for sub-periods B and CD.

The mathematical models proposed in the literature can all be described with the general framework we provided so far. However, they are different in terms of ingredients and relevant variables (Fig. 12.3). Specifically, they can be grouped into two broad classes with fundamentally different views on the role of DNA replication, its impact on cell division control, and ultimately on how the cell division and replication cycles are coupled [483, 484, 470, 476, 477]. A class of 'replication-centric' models see the completion of DNA replication as the crucial checkpoint for cell-cycle progression, which fundamentally limits division and initiation events [475, 477]. Instead, 'replication-agnostic' models assume that cell division is limited by a cell cycle-related process such as septum or cell wall formation and not by DNA replication [485, 476].

The linear-response theory over sub-periods coupled with the new-generation experimental observations on single cells gives us a powerful tool to compare the different models (see Mathematical Detail Box 12.D). Crucially, while the slopes of the size-growth plots are ultimately correlation patterns, the interpretation of the causal link between them changes across different models. For instance, the replication-centric models generally assume that two parameters among $\lambda_I$, $\lambda_B$, $\lambda_{CD}$ are input variables, fixed by an underlying molecular mechanism, while $\lambda_G$ is an output of the model, i.e. an emergent correlation pattern predicted by the model. In contrast, the replication-agnostic models assume a mechanism for the G-period ($\lambda_G$ is fixed), and the other correlation patterns are outputs of the model. Hence, the observed relationships between linear-response constants across conditions can be used to select a specific model. In the following, we present replication-agnostic theories first, then replication-centric models, then we introduce a class of models that find a solution of this dichotomy.

The replication-centric models are in line with the classic views on the E. coli cell cycle, but they are challenged by recent findings [468, 488, 483, 485, 478]. The 1968 Cooper and Helmstetter model was based only on the available population-average data at that time. The model posits that cell division happens within a defined period (CD) of time after initiation. Shortly after, Donachie [488] combined the Cooper and Helmstetter observation of a constant (population average) CD period with the even older observation that population-average cell size increases with the growth rate with a trend that is compatible with an exponential (the "Schaechter law" [465], which we mentioned above) and postulated that the population-average mass-per-origins is constant with the growth rate. Crucially, the classic paradigm by which replication limits division rested on indirect conclusions based on population averages, but these assumptions needed to be verified by single-cell data, which showed that things are much more complex [444].

In recent times, Ho and Amir [486] were the first to connect the Cooper-Helmstetter-Donachie ideas with the new observation of adder correlation patterns over the G-period. The authors assumed an adder mechanism during the I-period and a timer mechanism during the CD period. This model produces (in the limit of small noise in the timing of the CD period) an adder behavior in the G-period. Note that in this model $\lambda_I = -0.5$ and $\lambda_{CD} = 0$ are inputs while $\lambda_G \approx -0.5$ is an output of the model. This model, by definition, fails in reproducing the adder behavior in the CD period (which was not known at the time). Although it turned out to be an oversimplification, this work has the merit of connecting the old theories with new single-cell data into a simple and elegant replication-centric model.
## 11.3 Dynamic optimization of enzyme expression in metabolic pathways

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

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

$$ 
\begin{aligned}
\frac{dx_0}{dt} &= -k_1 e_1 x_0, \\
\frac{dx_i}{dt} &= k_i e_i x_{i-1} - k_{i+1} e_{i+1} x_i, \\
\frac{dx_n}{dt} &= k_n e_n x_{n-1},
\end{aligned}
\qquad (11.4) $$

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

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

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

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

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

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

The first experimental demonstration of the just-in-time principle was presented by Zaslaver and colleagues [299]. This work employed luminescent and fluorescent reporters to measure the temporal adaptation of *Escherichia coli* upon withdrawal of amino acids from the growth media. Clear just-in-time patterns of enzyme expression were found in the serine, methionine and arginine biosynthetic pathways. To better understand such patterns, the authors studied a model for an unbranched pathway with three
Mathematical details 12.D: Comparing cell-cycle sub-periods models with data

This box describes the quantitative tools necessary to systematically compare cell-cycle sub-periods models with data using the linear-response formalism and size-growth plots. Since the formalism may become very heavy, to avoid complications we will present the case of slow-growth conditions, in which there are no overlapping replication rounds. In addition, we will assume that the growth rate is a constant parameter and we will assume perfectly symmetric division.

Replication-centric models assume $\lambda_{CD}$ and either $\lambda_B$ or $\lambda_I$ to be input parameters in the model. Here, we focus on the case in which $\lambda_{CD}^*$ and $\lambda_I^*$ are fixed, which is the case for the Cooper and Helmstetter, Ho and Amir, and Witz et al models [468, 486, 487]. In these models, one has that $\delta q_I^{i+1} = (1 - \lambda_I^*)\delta q_I^i + \alpha \nu_I^i$ and $\delta q_0^{i+1} = (1 - \lambda_{CD}^*)\delta q_I^i + \alpha \nu_{CD}^i$, where $q_0^i$ and $q_I^i$ are the logarithmic sizes at birth and initiation of the cell cycle $i$, respectively; $\alpha$ is the growth rate, and $\nu_I^i$ and $\nu_{CD}^i$ are the white noise contribution related to the I and CD periods, respectively. In this class of models, $\lambda_G$ and $\lambda_B$ are mathematically linked to $\lambda_{CD}^*$ and $\lambda_I^*$, which provides predictions that can be validated or falsified with data:

$$ (1 - \lambda_G) := \frac{\langle \delta q_0^{i+1} \delta q_0^i \rangle}{\sigma_{q_0}^2} = \frac{(1 - \lambda_{CD}^*)^2 (1 - \lambda_I^*) \sigma_{qI}^2}{\sigma_{q_0}^2}, \qquad (12.18) $$

$$ (1 - \lambda_B) := \frac{\langle \delta q_I^i \delta q_0^i \rangle}{\sigma_{q_0}^2} = \frac{(1 - \lambda_{CD}^*) (1 - \lambda_I^*) \sigma_{qI}^2}{\sigma_{q_0}^2}. \qquad (12.19) $$

Note that by combining (12.18) with (12.19), we also get the relationship

$$ (1 - \lambda_G) = (1 - \lambda_{CD}) (1 - \lambda_B) . \qquad (12.20) $$

**Problem 12.9** Write the predicted $\lambda_G$ and $\lambda_I$ for a model in which $\lambda_{CD}^*$ and $\lambda_B^*$ are input parameters of the model. Does Eq. (12.20) still hold?

**Problem 12.10** Extend the previous models for:

* Overlapping rounds of DNA replication. This case is more difficult to address analytically, but can be easily simulated.

* The $\zeta$-formalism (without overlapping rounds). Use the model to answer the question: can an adder in the I- and CD-period provide the adder behavior in the G-period$^a$?

The first studies measuring the initiation of DNA replication in single cells [474, 475] brought two new experimental pieces of evidence into the field: they observed the duration of the CD period was inversely proportional to the single-cell growth rate and that the C period does not display any size compensation. Based on their data, Wallden and coworkers proposed a replication-centric model with a sizer in the B-period ($\zeta_B = -1$), which was later falsified [481, 477, 476]. A subsequent study by a different group [477] measured consecutive initiation events in single cells and observed three adders in the G, I, and CD periods. They then designed an improved version of the Ho-Amir model (already proposed for mycobacteria [489]) in which the initiation of DNA replication triggers both the next initiation and a division event with an adder mechanism. In this model, the adder in the G-period is an output of the model, which emerges from the adder in I and CD when the growth rate is a random variable and a sufficiently skewed asymmetry in cell division is added into the model. This replicaiton-centric model is unable to capture the growth rate – CD period inverse relationship discovered by Wallden and coworkers. However, it has the merit of improving the Ho and Amir model accounting for both adders in I and CD and introducing a debate over the importance of asymmetric division.

The replication-agnostic models entered the debate more recently. Based on dynamic cell-wall and cell-geometry mea-

<footnote>$^a$Note that the adder behavior can be recovered introducing asymmetric divisions [477]</footnote>
enzymatic steps and Michaelis-Menten kinetics:

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

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

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

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

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

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

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

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

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

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

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

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

In problem (11.11), the constant $F_{\text{goal}}$ is a prescribed production flux that the pathway should achieve at steady state. Minimization of the objective in (11.11) accounts for the activation of the pathway from an "off" state until it reaches a prescribed flux $F_{\text{goal}}$. This formulation differs from the previous example [301] in two important ways. First, it accounts for cellular resources in the objective function itself. The first term of the objective quantifies the total amount of repressor produced through the optimization horizon, and thus relates to the amount of cellular resources required to activate the pathway. Second, the decision variables are the regulatory parameters, not the temporal profiles of the molecular species. Therefore, strictly speaking, this is not an optimal control problem but rather a static optimization problem subject to dynamic constraints encapsulated by the pathway ODE model. Through numerical solutions for different values of the protein cost weight $a$ and optimization horizon $T$, the authors determined conditions under
surements, Harris and Theriot proposed a model in which the completion of the division septum, and not the chromosome, was the limiting factor for cell division [484, 485]. This model proposes a simple molecular mechanism for the adder based on three main ingredients: (i) a crucial factor involved in setting division is produced at a rate proportional to the cell size; (ii) this factor needs to reach a threshold in the number in order the cell to divide; (iii) the factor in the next generation has to be reset, with no history dependencies on the previous cell cycle (in the case of the septum, this is natural, as a new septum needs to be produced from zero at every cell cycle). This model structure is still the basis for different mechanistic models explaining the adder during the G period, but the mechanistic factor was also proposed to be a protein [490, 491, 476]. Further evidence in favor of a replication-agnostic view came from experiments performed by the Jun lab [476] aiming to perturb independently the adder correlation pattern in the G-period, while maintaining intact the adder pattern over the I-period, and viceversa. The perturbations were achieved by inducing oscillating levels of the FtsZ protein, which forms a contractile ring structure at the future cell-division site and of the DnaA protein, responsible for the initiation of replication, respectively. The authors interpreted the results of these experiments as a proof that the replication and division cycles are independently regulated, and in particular that completion of DNA replication and segregation is not a limiting factor for cell division. Additionally, the authors re-interpreted the ‘molecular adder’ model proposed by Harris and Theriot, suggesting that the FtsZ may be the “adder protein” setting division. This work has the merit of providing precious experimental information. However, the model fails to explain the adder behavior over the CD period, as well as the correlation patterns related to how the replication and the division cycles are coordinated [481, 483, 466].

The replication-centric and replication-agnostic views have been firmly opposing each other in recent years (see e.g. [492, 487, 493]). However, a standpoint that is gaining consensus is that neither of these views is able to capture the full complexity of the correlation patterns in the data [481, 483, 466, 478, 470, 494]. The recently proposed “concurrent-cycles” scenario [481, 483, 466] bridges the two opposing views and is in better agreement with the data compared to all the above models. The key innovative element in this theoretical framework lies in the assumption that there is no unique process limiting cell division. Rather a set of competing processes have to be completed before division, and some “downstream control” module (modelled as a logic gate) has to process the input from these processes. In its original formulation [481, 483], based on the available data the competing processes are the DNA replication processes defined by an adder in the I-period, a timer in the replication-segregation period cycle, and a cell division process that adds constant size between two consecutive divisions (division-related cycle). The division is decided by an AND gate, which triggers when both of two actions are completed, the interdivision period is complete and the replication-segregation period is complete. Therefore, the AND gate selects the slowest of the two random processes (which vary across single cells) to set the timing. Note that in this framework the CD period can be set by the intrinsic replication-segregation period of this is the slowest process, or by the interdivision period in case this other process is the slowest one. The concurrent-cycles framework makes precise predictions on how the sub-periods correlations of size change when either the replication-related or the division-related cycles are perturbed. Recently, experiments in which cell wall insertion is delayed confirmed the prediction of the model [466]. Other recent studies proposed similar frameworks, adding mechanistic details, where the onset of constriction at the divisome [494] and/or a “progression control complex” including the chromosome and the divisome play the role of the gate deciding cell division [478, 470]. Technically, concurrent cycle models need an additional set of parameters compared to the replication-centric and agnostic models (see Mathematical Detail Box 12.E). These parameters are ultimately summarized by one extra relevant parameter, which can be expressed as the probability that the division-related process to sets division (in a given cell cycle). Thus, the replication-centric and replication-agnostic models can be seen as limit cases of the concurrent-cycles framework, where this probability is zero or one respectively.

Despite the large improvement that the concurrent-cycles framework provides in the agreement with data, many questions remain open. For example, we do not know the probability of either of the concurrent processes limiting division varies under different conditions. Recent surveys of the available data [466, 494] suggest that the probability of a chromosome-agnostic cycle increases with increasing growth rate. At very slow growth (interdivision times of 300 minutes or more), it has been suggested that replication-segregation is the limiting process. Additionally, we currently do not know what tunes such probability and what the role of the growth rate may be. We also do not know how many concurrent processes there are and which precisely are the relevant players at the molecular level. Finally, the regulation of initiation
which the optimal solutions showed two features of the just-in-time property, namely:

$$ \tau_1 < \tau_2 < \tau_3, \quad \max_t e_1 > \max_t e_2 > \max_t e_3, \qquad (11.13) $$

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

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

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

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

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

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

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

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

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

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

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

where $F_{\text{goal}}$ is a (constant) target pathway flux, similar as in (11.11). The terminal constraint ensures that the pathway reaches a steady state at the final time $T$. Using Pontryagin's Minimum Principle [295], the authors showed that the optimal enzyme concentrations follow a bang-bang temporal profile that matches the order in which they act on the pathway. This result was shown to be independent of the weight $\alpha$, the number of enzymatic steps, and valid for a wide range of enzyme kinetics satisfying the assumptions
of DNA replication could also be set by a "gate" integrating a set of processes, a hypothesis that remains underexplored in the literature.

## Mathematical details 12.E: The concurrent-cycles framework

This Mathematical Detail Box provides the mathematical relationships that correspond to the ones appearing in Mathematical Detail Box 12.D for the more general concurrent-cycles framework. Given the complexity of this model, we restrict to the case of no overlapping rounds. In particular, we will show how Eq. (12.20) is no longer valid in the concurrent-cycles framework (without the need to include additional ingredients such as asymmetric division or mother-daughter growth rate correlations).

In the concurrent-cycles model, cell division is determined by the slowest of two processes. The first process is an interdivision, (chromosome-agnostic) cycle that is concluded, for generation $i$, at a log-size $q_H^i$, which is expressed as $q_H^i = q_H^* + (1 - \lambda_H^*) (q_0^i - (q_H^* - \log 2)) + \alpha \nu_H^i$, with $\lambda_H$ size control parameter of this process. The second process is a chromosome replication-segregation cycle (replication-centric), that is concluded, for generation $i$, at a log-size $q_R^i$, which is expressed as $q_R^i = q_R^* + \delta q_I^i + \alpha \nu_I^i$. Note that this equation assumes a timer for this process, $\lambda_{CD'}^* = 0$, where $CD'$ identify the time needed for completing DNA replication, which is identical to the measurable CD-period only when this second cycle sets division. The cell size at division is determined by the slowest process, i.e. $q_d^i = \max(q_H^i, q_R^i)$. The initiation of DNA replication decides the next initiation independently on the size at birth or division, generating the fluctuation around the logarithmic size at initiation that we already found in Mathematical Detail Box 12.D, $\delta q_I^{i+1} = (1 - \lambda_I^*) \delta q_I^i + \alpha \nu_I^i$.

To calculate the fluctuations of the logarithmic size at division, we assume that the replication-centric process sets the division of generation $i$ with probability $p_H$ independently on $q_0^i$ and $q_I^i$. With this assumption, and considering $\lambda_H^*$, $\lambda_I^*$ and $\lambda_{CD'}^* = 0$, the model predicts the following values for the strength of the size-growth plots in the B-, CD- and G-period,

$$ (1 - \lambda_B) = (1 - \lambda_{CD})(1 - \lambda_I) \frac{\sigma_{q_I}^2}{\sigma_{q_0}^2} \qquad (12.21) $$

$$ (1 - \lambda_{CD}) = (1 - p_H) + p_H (1 - \lambda_H^*) (1 - \lambda_B) \frac{\sigma_{q_I}^2}{\sigma_{q_0}^2}, \qquad (12.22) $$

$$ (1 - \lambda_G) = (1 - p_H)(1 - \lambda_B) + p_H(1 - \lambda_H^*) \qquad (12.23) $$

Overall, the concurrent-cycles model allows to match the experimental trends in the size-growth plots with an additional parameter ($p_H$). In particular, it allows to break the relationship in Eq. (12.20) without including asymmetric divisions or mother-daughter correlations in growth rates [481, 483, 466].

## 12.5. Protein sectors and cell division

This chapter focuses on quantitative descriptions of the cell cycle and cell division control, and it is natural to wonder whether and how these consideration relate to the topic of previous Chapters 8 and 9 which deal with resource allocation models where cell growth is set by catabolism and biosynthesis. There is a strong link between regulation of growth and cell-cycle progression, which remains a largely open area of investigation both in biology and in quantitative biology / physics of living systems. This section discusses some recent models aimed to describe some specific aspects of the coordination between cell growth and cell-cycle progression. We will start by presenting the main questions that we want to address with the aid of mathematical models. Then we will discuss the main ideas and ingredients behind the models that address these questions, and present some relevant predictions that can be tested and validated against experimental data.

The maintenance of an interplay between cell growth and cell cycle is crucial for the correct functioning of the cell.
<figure><img src="image_107.png" /></figure>

<figure><img src="image_108.png" /><figcaption>Figure 11.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 (11.15), thus extending the original finding in [301] to a larger class of pathways. Figure 11.1 shows a numerical example of the optimal activation pattern obtained for an unbranched metabolic pathway of length three (see also Exercise 11.1).

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.

## 11.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.
Specifically, a cell has to adapt both growth and division rates concertedly when either one is perturbed. For example, the response and adaptation to environmental stresses, such as sudden shifts in nutrient conditions or exposure to drugs or toxins, requires the ability to reprogram in a coordinated way cell growth and cell division. Consequently, cells across all kingdoms of life have developed specific mechanisms to precisely coordinate cell cycle progression with cell growth and biosynthesis [196, 495, 447, 496, 451, 497]. There are many mechanisms involved in this coordination, and we lack a complete and coherent quantitative understanding of how this coordination works in different contexts. Sometimes we even lack simple ways to frame questions concerning the effects on cell cycle progression of cell growth perturbations/inhibitions, or the effects of cell growth of cell-cycle perturbations (such as cell cycle arrest).

To formulate and address these questions quantitatively, we would need a theoretical framework where both growth physiology (as in “how does a cell grow?”) and cell-cycle decisions/progression (as in “how does a cell decide when to divide?”) aspects are allowed to play a role and influence each other. However, while both cell growth and cell cycle progression alone have been subject of intense study in the past (especially in bacteria [445], but see ref. [498] for a recent review of these themes in eukaryotes), comparatively little effort has been directed so far toward the development of such unified framework. Nonetheless, recent work has advanced our quantitative understanding of the cross-talk between cell growth and cell cycle progression in bacteria. The remainder of this section will focus on discussing these aspects.

Relatively to the bacterium *E. coli*, recent and current efforts aimed at integrating already existing coarse-grained models of cell physiology and cell cycle control. More precisely, several studies have extended the classic proteome allocation theory, (presented in chapters 8 and 9), which has proven successful in describing several physiological laws, to include also a cell-division proteome sector “X”, whose dynamics should implement cell-division control (or cell-cycle progression control) strategies at a phenomenological or molecular level (Fig. 12.4). The current models for *E. coli* usually include a threshold accumulation process for cell division, i.e., proteins of the division sector accumulate during cell cycle progression up to a threshold level that triggers cell division. The previous section has mentioned some candidate molecular players for this accumulation (the FtsZ protein and the cell wall insertion).

Let us take a closer look at the ingredients of this modeling framework. The two main ingredients are (i) the standard proteome allocation theory extended to include a division sector *X*, alongside to the standard main sectors (see Chapter 8), *Q* (house-keeping), *R* (ribosomes), *P* (catabolism and transport), together with (ii) a threshold-accumulation division strategy to set the decision to divide (Fig.12.4A). Note that the fact that the division factor *X* is a protein is an implicit assumption in these framework and experimentally things could be more complex. Crucially, the fact that cell division is a proteome sector couples the rates of cellular growth and division, by controlling the synthesis of division proteins. Specifically, the models encode a trade-off between ribosomes and division protein synthesis, which as we will see determines many salient predictions.

The following Mathematical Detail Box shows how these ideas and ingredients can be translated into a mathematical model. The framework that we are now going to discuss is consistent with different models recently developed in the literature [499, 500, 501, 502, 496].

## Mathematical details 12.F

The model consists of two different layers of dynamical equations, and one relationship connecting them. The first set of equations describes cell growth and division as cellular processes

$$ \frac{ds}{dt} = \lambda s, \quad \frac{dX}{dt} = k_X s - \frac{d_X}{m_X} X, \quad X(\tau_d) = X_{th} \implies \begin{cases} s(\tau_d) \to s(\tau_d)/2 \\ X(\tau_d) \to 0 \end{cases}, \quad (12.24) $$

where cell size *s* (mass or volume) grows exponentially at a rate $\lambda$ ($[\lambda] = [T]^{-1}$), while division proteins *X*, of mass $m_X$ being synthesized and degraded at rates ($k_X$ ($[k_X] = [s]^{-1}[T]^{-1}$), $d_X$ ($[d_X] = [M][T]^{-1}$)), accumulate until a threshold amount of them is reached and cell division occurs, after that cell size is divided exactly in half and division proteins number is reset to zero.
We consider the class of models with dynamics given by Eq. (11.3), 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 reaction rate function:

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

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 (11.20)$$

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

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

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. (11.20) 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 (11.22)$$

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 (11.23)$$

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 (11.24)$$

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 11.2 gives an example of a resource allocation model, describing a simple self-replicatory microbial system [306, 307] (see Chapter ?? for 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
**Problem 12.11** Run numerical simulations of Eqs. (12.24). Prove that in order to obtain an adder, the ingredients of a size-specific (rather than constant) production rate of the division protein $k_X$ and a reset to zero (rather than partitioning in half in the two daughter cells) of the division factor $X$ turn out to be essential.

The second set of equations describes the dynamical allocation of the proteome and the biosynthesis layer underlying cell growth, as follows

$$ \begin{aligned} \frac{dA}{dt} &= \frac{1}{m_a} \left( k_n P - ak_t R f_a + \sum_{P_i \in \{Q, P, R, X\}} d_{P_i} P_i \right), \\ \frac{dP_i}{dt} &= \frac{1}{m_{P_i}} (ak_t f_{P_i} R f_a - d_{P_i} P_i) . \quad P_i \in \{Q, P, R, X\}. \end{aligned} \quad (12.25) $$

According to Eq. (12.25), free amino-acids (A) are produced from import/catalysis of nutrients at a rate $k_n$ ($[k_n] = [M][T]^{-1}$) per number of catabolic/transport proteins P, and from protein degradation, occurring at a rate $d_{P_i}P_i$ (where $d_{P_i}$ ($[d_{P_i}] = [M][T]^{-1}$) is the degradation rate) for each specific sector. Free amino-acids are taken up to synthesise each proteome sector $P_i$ at a rate equal to the number of active ribosomes ($Rf_a$), times the fraction of ribosomes synthesising the specific sector $f_{P_i}$, times an overall protein translation rate, which in this particular model is equal to a constant translation rate per ribosomes $k_t$ ($[k_t] = [s][T]^{-1}$) times the concentration of free amino-acids $a \equiv (m_a A) / ([a] = [M][s]^{-1})$.

**Problem 12.12** Rewrite the system of equations (12.25) in terms of protein fractions, either defined as protein mass fractions $\phi_i \equiv M_i/M_{prot}$ or protein number fraction $\psi_i \equiv P_i/\sum_i P_i$, where $M_{prot} = m_Q Q + M_P P + m_R R + m_X X = M - m_A A$. In both cases one has the obvious constraint $\sum_i \psi_i = 1 = \sum_i \phi_i$. Find the connection between $\psi_i$ and $\phi_i$. What can be generally said about the stationary composition of the proteome? How does the scenario change if degradation can be neglected?

Finally, there must be a connection between the two levels of description, in the sense that cellular rates should be regarded as the result of the underlying biosynthesis dynamics. To make this connection explicit, we write the equation

$$ s = \gamma M = \gamma (m_A A + m_P P + m_R R + m_Q Q + m_X X), \quad (12.26) $$

representing mass conservation (if "size" stands for "mass" $s = M$), or the assumption of constant density (if "size" stands for "volume" $s = V$), verified in E. coli for population averages but not for single cells, or for certain perturbations [448, 503].

Together, Eqs.12.24, 12.25 and 12.26 fully specify the model.

**Problem 12.13** For the mathematically curious readers, show that the model described so far can be written in more general mathematical terms as

$$ \begin{aligned} \frac{dX_i}{dt} &= f_i(\mathbf{X}); \qquad \frac{dZ}{dt} = h(\mathbf{X}, Z) \\ V(\mathbf{X}, Z) &= \sum_{i=1}^{N} v_i X_i + v_Z Z \end{aligned} \quad (12.27) $$

where $V$ is the volume of the cell and $X_i$, $Z$ its chemical constituents. Identify the functions $f_i$s and $h$. Show that the $f_i$s satisfy the property of homogeneity, $f_i(\beta \mathbf{X}) = \beta f_i(\mathbf{X})$. The predictions of this model have been studied in the wider framework of dynamical systems theory [504, 505].

In order to exemplify how this framework can generate relevant predictions, the following Mathematical Detail Box shows some concrete examples from the literature. The mathematical derivations are not exhaustive, but aimed to give
<figure><img src="image_109.png" /><figcaption>Figure 11.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>

is assumed saturating.

The resource allocation functions in the model are not explicitly specified by regulatory mechanisms, but assumed to 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 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 (11.25) $$

where like in the general case of Eq. (11.1), U denotes the set of admissible profiles for the resource allocation functions u. Note that maximization of biomass over the time-interval $[0, T]$ amounts to taking the integral of the instantaneous growth rate over that time-interval, defined by Eq. (11.24). This objective

### Dynamical system

$$ \begin{aligned} \dot{c} &= v_M - \mu c \\ \dot{p}_Q &= u_Q v_R - \mu p_Q \\ \dot{p}_R &= u_R v_R - \mu p_R \\ \dot{p}_M &= u_M v_R - \mu p_M \end{aligned} $$

### Synthesis rates

$$ \begin{aligned} v_M &= k_M p_M \\ v_R &= k_R p_R \frac{c}{c + K_R} \\ \mu &= \frac{v_R}{P_Q + P_R + P_M} \\ u_Q + u_R + u_M &= 1 \end{aligned} $$

### Objective function

$$ J(u) = \int_{0}^{T} \mu(t) dt $$

### Optimal control

$$ u_{\text{opt}} = \arg \max_{u \in \mathcal{U}} J(u) $$
<figure><img src="image_81.png" /><figcaption>Figure 12.4: Ingredients and predictions of modeling frameworks integrating sector models with cell-division control. (A) The framework unifies growth and cell division by extending the standard proteome allocation model to include a division sector X, implementing a threshold-accumulation process setting the decision to divide. (B) The model is naturally suited to uncover general relationships and growth laws involving proteome composition and growth rate, as well as trade-offs between different proteome sectors. The inclusion of a division protein sector X regulating cell division allows the model to make predictions on cell size control and study the transient dynamics in nutrient shifts.</figcaption></figure>

the reader a feeling of the "recipe" followed to obtain a given prediction starting from the model's ingredients. The interested reader should have sufficient information to work out the mathematical calculations autonomously or follow the complete derivation in the cited references (for example by Serbanescu et al. [499, 500]).

## Mathematical details 12.G

**Growth laws and trade-offs between protein sectors.** Prototypical predictions are the so-called "growth laws", general quantitative relationships linking proteome composition and rates of cellular processes. The reason why relationships of the kind $\lambda = \lambda(\phi_R, \phi_X, \dots)$ and $k_X(\phi_R, \phi_X, \dots)$ naturally emerge in the framework is due to cell growth and division rates being coupled to proteome allocation dynamics.

**Growth law for the ribosome sector.** For example, the first growth law, stating that the ribosome mass fraction increases linearly with the nutrient-imposed growth rate, that is $\lambda = \lambda(\phi_R) = K(\phi_R - \phi_R^{min})$, is obtained straightforwardly by noting that upon differentiation of Eq. (12.26) with respect to time and substitution of Eq. (12.24) and Eq. (12.25) one finds the dynamical relation $\lambda(t) = \frac{k_n P(t)}{M}$, which at equilibrium reads (neglecting degradation)

$$ \lambda^* = \frac{k_n P^*}{M} = \frac{ak_t R^* f_a}{M} = \frac{ak_t M_{prot}}{m_R} (\phi_R - \phi_R^{min}), \quad (12.28) $$

since at equilibrium the amino-acid import flux $k_n P^*$ matches the biosynthesis flux $ak_t R^* f_a$ ($dA/dt = 0$ in absence of degradation). Note that we have used the definitions $\phi_i \equiv (m_i P_i)/M_{prot} = (m_i P_i)/(M - M_a)$ and $Rf_a = R^{active} = R - R^{inactive}$ and we have identified $\phi_R^{inactive} = \phi_R^{min}$.

**Trade-Offs between Ribosomes and Division Protein Synthesis.** Following Refs. [499, 500], we re-write Eq. (12.28) as $k_n = \frac{ak_t}{m_R} m_P \frac{\phi_R - \phi_R^{min}}{\phi_P}$ and use the constraint $\phi_R^{max} = 1 - \phi_Q = \phi_R + \phi_P + \phi_X$ to obtain

$$ \phi_X = - \frac{K_n + K_t}{K_n} \phi_R + \frac{K_t \phi_R^{min} + K_n \phi_R^{max}}{K_n}, \quad (12.29) $$
does not generally reduce to growth rate maximization, that is, maximization of the instantaneous growth rate at every time-point of this interval. In particular, there may be situations where a lower-than-optimal growth rate over some sub-interval of $[0, T]$ may prepare for a higher growth rate after a sudden change in conditions, and thus turn out to be beneficial for the total accumulation of biomass over $[0, T]$. The dynamic perspective on microbial growth in this chapter thus entails a generalization of the objective criterion in comparison with previous chapters, where balanced growth of microorganisms at steady state was considered.

The question can be asked, for the microbial self-replicator in Figure 11.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 11), 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 11.2. For the simple example considered here, the optimal solution can be characterized analytically [308, 306, 307]. This is not possible for more complicated examples, however, which require the optimal solution to be constructed numerically. Figure 11.2C-D show a typical solution for parameter values estimated from experimental data [307]. 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 [306, 307]. 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 (see Exercise 11.2).

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 of ensuring 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 11.2, there is evidence that during balanced
where $K_n \equiv k_n/m_P$ ($[K_n] = [T]^{-1}$) and $K_t \equiv ak_t/m_R$ ($[K_t] = [T]^{-1}$). Eq. (12.29) shows a negative correlation between the ribosome and division sectors under nutrient or translational perturbations, in agreement with recent published data [506]. Also, since the rates of growth and division protein synthesis are respectively proportional to the ribosome and the division sector, this negative correlation reflects a trade-offs between allocating ribosomal resources towards growth or division (see Fig.1F in Ref. [499]).

**Growth law for the division sector.** So, the larger the fraction of ribosomes making division proteins the smaller the fraction of ribosomes making ribosomes. In other words, there is a negative correlation between the growth rate and the division protein sector. Indeed, the ribosome sector is related to the growth rate via the first growth law $\phi_R = \frac{M\lambda}{M_{prot}K_t} + \phi_R^{min}$, but it is also related to $\phi_X$ via Eq. (12.29) $\phi_R = \frac{K_t\phi_R^{min}+K_n\phi_R^{max}}{K_n+K_t} - \frac{K_n}{K_n+K_t}\phi_X$. Equating the two terms yields

$$ \lambda = \frac{K_n K_t}{K_n + K_t} \frac{M_{prot}}{M} (\phi_R^{max} - \phi_R^{min} - \phi_X), \qquad (12.30) $$

which is Eq. (9) in Ref. [499].

We now discuss how two known steady-growth size-related behaviors emerge in the unified framework from the interplay between cell growth and cell division.

**Adder mechanism.** As we discussed, E. coli cells regulate their size by adding a constant volume between consecutive cell divisions (adder mechanism). In a previous problem, we investigated with numerical simulations the range of validity of this property. In the following one, we instead show analytically that the adder property is naturally embedded in the unified framework.

**Problem 12.14** By directly integrating Eq. (12.24), derive the following expression for the threshold number of division proteins $X_{th} \equiv X(\tau_d)$

$$ X(t) = \frac{k_X s_0}{\lambda + \frac{d_X}{m_X}} \left( 2^{\frac{t}{\tau_d}} - 2^{-\frac{d_X}{m_X \lambda} \frac{t}{\tau_d}} \right) \implies X_{th} = \frac{k_X}{\lambda + \frac{d_X}{m_X}} \left( s_d - s_0 2^{-\frac{d_X}{m_X \lambda}} \right). \qquad (12.31) $$

It can be seen then that whenever $\lambda \gg d_X/m_X$ (e.g. fast growth conditions), $\Delta s_{1cycle} \approx \frac{\lambda}{k_X} X_{th} = \text{const}$ which is the adder property. Notably, in increasingly slower growth conditions, where degradation becomes with the growth rate, deviations from the adder are predicted, up to the point $\lambda \ll d_X/m_X$ where $s_d \approx X_{th}d_X/(k_X m_X) = \text{const}$.

**"Schaechter-Maaloe-Kjeldgaard" (SMK) growth law.** According to this law, the population-averaged cellular size scales with growth rate in an approximately exponential fashion [507]. Interestingly, deviations from the exponential trend have recently been reported, particularly at slow growth, leading to a different proposition [478]. Notably, deviations from this law are accounted in our framework. Indeed, in an exponentially expanding population the average cell size can be expressed as $\langle s \rangle = 2 \log 2 \langle s_0 \rangle$ [445], which, combined with Eq. (12.31) and $\langle s_d \rangle = 2 \langle s_0 \rangle$ leads to

$$ \langle s \rangle = \frac{\lambda + \frac{d_X}{m_X}}{\tilde{k}_X \left( 2 - 2^{-\frac{d_X}{m_X \lambda}} \right)} \qquad (12.32) $$

where, following Ref. [499], we have defined $\tilde{k}_X \equiv k_X/(2 \log 2 X_{th})$. Note that since $\lambda \propto \phi_R$ and $k_X \propto -\phi_R$ the average cell size increases with ribosome abundance, a trend observed in experiments. Notably, upon determining the model parameters and making full explicit the growth rate dependence, the authors in Ref. [499] with no further fitting showed that Eq. (12.32) recapitulates the experimental data [478, 499], a remarkable achievement of the unified framework.

**Non-steady relationships.** Finally, we contextualize within the unified framework some predictions of a model recently proposed to unify cell division and growth in non-steady growth conditions [491]. As we saw, although there is consensus on an inter-division adder at the phenomenological level, the mechanisms regulating cell division
growth, microorganisms find a trade-off between maximizing growth rate in a given environment and minimizing necessary adjustments to other environments [297]. 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 11.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. (11.22)). Bearing in 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, 307].

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 (Box 11.A).

Box 11.A Temporal optimization for biotechnological applications and process design

As explained in Section 11.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 [314]. 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 [315]. 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 [316]. 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 [315]. The conclusion that this two-phase procedure is optimal corresponds well with established practice in biotechnology [317]. Very similar conclusions were attained in related work by Jeanne et al. [318].

## 11.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 metabolic activity of the cells, and track the resulting overall biomass growth. It is noted that a basic assumption of dFBA is that organisms rapidly reach intracellular steady state in response to extracellular perturbations, and on the long run no metabolite can accumulate or deplete.
dynamics in the bacterium *E. coli* are still widely debated. In particular several mechanistic models based on different mechanisms for division control were proposed for the adder [485, 422, 476, 475]. In order to help selecting different scenarios, experiments beyond steady-state growth help comparing the specific causal relationships underlying different models with data. Following this philosophy, and aiming to shed more light on cell division dynamics, Panlilio et al. [491] ran multiple long-term *E. coli* microfluidics experiments jointly monitoring size-division dynamics and reporters of ribosomal and constitutive genes through nutritional up-shifts. The fluorescent reporters can be seen as proxies for the dynamics of the *R* and *P* sectors during the shift. Remarkably, in their experiments they observed highly-complex multiple-timescale dynamics in different cell-division variables (particularly inter-division time, division rate, added volume and added-to-initial volume ratio) during the nutritional up-shift. Notably, in spite of this complex dynamics, they found the division control strategy to be unaffected by the shift. The transient observed division dynamics in their shift data falsifies several scenarios, such as the Harris-Theoriot septum-limited division and the classic scenario of replication-limited division. Instead, the authors found that a threshold accumulation model such as the one described by Eq. (12.24) could not be falsified,

$$ \frac{ds(t)}{dt} = \alpha(t)s(t) \qquad \frac{dN(t)}{dt} = r_X(t)s(t). \qquad (12.33) $$

This the usual scenario where a constitutive *X*-sector protein accumulates to a threshold value $N^*$ and at that point triggers cell division. The regulation of cell division from a constitutive sector is coherent with the observation that ppGpp is a cell size and cell division regulator [508]. These results are also in line with independent conclusions based on steady-state data [500, 476, 490] and isolate FtsZ as a likely candidate cell-division trigger, although the previous section has clarified how the complexity of the decision to divide is likely higher than described by the chromosome-agnostic cell-division models that are used in integrated frameworks. Future efforts will have to integrate this complexity in a description that also accounts for the interplay of different processes relevant for cell cycle progression with cell growth.

## 12.6. Control of cell division across species and kingdoms

The concepts described in the previous sections are widely applicable, but there are many relevant species-specific aspects, so that different crucial assumptions that we have taken so far might break down for different species and kingdoms. Additionally, it should be noted that the approach described here is purely phenomenological, while a biological investigation might be concerned with the detailed molecular players responsible for the cell division and cell-cycle progression decisions. Even in this case, the approach is useful and is being applied in recent work. For example, if the goal is to understand how the size control phenomenon is regulated, the phenomenological analyses can quantify how the phenomenology of size correction behaves under different mutants and perturbations, helping to identify molecular players and their effects on cell-cycle decisions.

Let us consider briefly some important variations of the approach used so far, relevant for the understanding of different species-specific behaviors. First, it is not granted that single cells grow exponentially, or even that exponential growth is a good approximate description. Even in the cases where exponential growth appears to be a good average description, these averages may emerge from more complex behaviors at the single-cell level or in cell cycle sub-periods. For bacteria, most studies conclude that exponential growth is a sufficiently good description, although recent accounts show deviations [509, 510]. In budding yeast (*S. cerevisiae*), the average growth rate was reported to change at regulatory checkpoints with the cell-cycle phase [511, 512, 513]. In the fission yeast *S. pombe*, a systematic study of single-cell growth concludes that the majority of growth trajectories are best described by a bi-linear growth [514]. In cell lines of animal cells, most studies suggest that, on average, cells grow exponentially until a certain saturation size after which they slow down, but this mean behavior hides many details [447]. For example, it seems that cells in the G1 phase of the cell cycle grow at a slightly slower rate than in later stages of the cell cycle [515].

A second important aspect to consider is whether division is symmetric or not. In *E. coli*, cells divide symmetrically,
In general, a dFBA model comprises three main parts as demonstrated in Figure 11.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 (11.26) $$

where $X$ denotes the biomass concentration, typically measured as dry mass in g L⁻¹, and $\mu$ denotes the growth rate, typically measured in 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 (11.27) $$

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 (11.27), e.g., intracellular metabolites. This steady state constraint is given by

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

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 (11.29) $$

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

In recent years, dFBA is increasingly applied for the simulation of dynamic biological systems, especially due to the promising use of genome-scale metabolic models (GSMMs) for interpreting cell physiology and evolution, as well as for guiding metabolic engineering and bioprocess design and optimization [319, 320, 321]. The dFBA applications based on GSMMs include the bacteria *Escherichia coli* [322] and *Lactococcus lactis* [323], as well as the yeast *Saccharomyces cerevisiae* [324, 325, 326, 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* [328, 329, 330, 331] and *Corynebacterium glutamicum* [332], models of yeast [333], but also plant and animal models, such as a model of the core metabolism of *Arabidopsis thaliana* [334].

It is noted that most of the dFBA applications for microorganisms simulate microbial fermentations un-
giving rise to two daughter cells that are nearly equal in size, with a precision of a few percent [450]. However, different species use very different strategies for cell division, which increase variability or explicitly aim for asymmetry. For example, S. cerevisiae reproduces through budding (hence the term “budding yeast”). The parent cell creates a small outgrowth that eventually becomes a daughter cell. Both division strategies are common among unicellular organisms (many filamentous fungi grow via budding). Budding creates a parent/offspring distinction in which age-related aspects are not transmitted equally. Since aging may correspond to a decrease in fitness/growth rate, it can also create diversity along lineages. A third important aspect to consider is that the growth rate may be coupled to size and enforce size homeostasis. In other words, homeostasis can be achieved by modulating cell-cycle duration based on size at birth, but also if large-born cells grow slower than small-born ones.

As an example of how different issues can be analyzed with extensions of the phenomenological approaches discussed so far, it is instructive to discuss in more detail how one can use the linear-response framework to detect indications of growth-based size homeostasis. As we mentioned previously, the overall multiplicative growth of a cell in one cycle is quantified by $G = q_f - q_0 = \log \frac{s_f}{s_0} =: \alpha\tau$. The slope $\lambda$ of the size-growth plot is equivalent to considering the conditional average of G over logarithmic size $q$,

$$ \langle G \rangle_q = \langle G \rangle - \lambda \delta_q \qquad (12.34) $$

As we have seen in Fig 12.1, we can consider the separate contributions of timing and growth to the coupling by taking separate scatter plots with growth rate and cell division time. We can give a more formal quantification of their contributions as follows. We call $\theta$ the coupling strength derived from the slope the first plot quantifying control by modulation of interdivision time,

$$ \langle \tau \rangle_q = \langle \tau \rangle - \langle \tau \rangle \theta \delta_q, \qquad (12.35) $$

and $\gamma$ the slope quantifying modulation of growth rate based on birth size,

$$ \alpha - \langle \alpha \rangle = -\langle \alpha \rangle (\gamma \delta_q) + \nu_{\alpha}. \qquad (12.36) $$

For positive values of $\gamma$, cells that are born larger than average can correct their sizes by growing with a slower growth rate, and cells that are born with a smaller size than average can correct by growing at a faster rate. Conversely, for negative values of $\gamma$, birth-size related specific growth rate variations increases systematically size variability.

Intuitively, we can understand that $\theta \gamma$ and $\lambda$ must be related. First, the overall homeostasis must be the result of the one enforced by growth-rate modulation and the one enforced by interdivision-time modulation. More formally, the slopes of the correlation plots illustrated in Fig. 12.1 for G, $\alpha$ and $\tau$ versus logarithmic birth size must be related, because $G = \alpha\tau$.

Using the linear response approach defined in section 12.3, one can derive the following equation

$$ \lambda = \theta \langle \alpha \rangle \langle \tau \rangle + \gamma \langle \alpha \rangle \langle \tau \rangle. \qquad (12.37) $$

Eq. (12.37) states that the overall correction to size over a cell cycle has to be the sum of a correction due to modulation of timing and a correction due to the modulation of specific growth rate based on size at birth. For example, if the overall strength is an adder, and the size coupling of the duration of the cell cycle is already an adder, the growth rate must be uncoupled from initial size.

Going back to the data, one can use Eq. (12.37) to evaluate the different strategies, by evaluating the couplings $\theta$, $\gamma$ and $\lambda$ from the different scatter plots. Importantly, the constraint imposed by Eq. (12.37) is realized in data from several bacterial species and growth conditions, indicating that the framework is sufficient to describe the data. Work on different bacteria shows widespread adder correlations [445], hence $\lambda \simeq 0.5$. What is more surprising is that adder behavior has been reported for in budding yeast and cultured human cells. Hence, for many species, the inter-division correlation patterns are nearly always close to an adder. One interesting exception is the fission yeast S. pombe, discussed below. The widespread adder patterns may suggest common general principles underlying the division control
<figure><img src="image_110.png" /><figcaption>Figure 11.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>

der 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 [324, 322, 325], while optimal control of batch or fed-batch operation of bioreactors is important for the production of desired chemicals [324]. Finally, dFBA has also been expanded for the study of microbial communities, where each microorganism is represented by a linear program (LP) that is solved independently [335, 336, 337]. 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 [328], it was not formalized until 2002 [329]. 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 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 [329].

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

4. feasibility: The intracellular optimization problem can become infeasible and lead to failure of the integration of the extracellular ODEs.
of microorganisms and cultured single mammalian cells. Considering the couplings $\theta, \gamma$ shows a different scenario, with a clear distinction between microorganisms and cultured mammalian cells. In the studied unicellular microbes, the inter-division adder is always due to the modulation of cell-cycle duration. Instead, cultured mammalian cells also rely on growth rate modulation to correct their size. In particular, this rejects the hypothesis that adder behavior may be favored by common underlying mechanisms. Additionally, for budding yeast and mammalian cells, the overall adder behavior emerges from homeostatic regulations acting close to the initiation of replication (G1/S transition) during the cell cycle, and from a weaker regulation of the subsequent parts of the cell cycle [447]. Cell growth outside of G1 is critical in setting the average cell size but appears to be less significant for the size homeostasis effect setting cell-to-cell variability in birth size. This is not the case in bacteria, where we have seen that key questions regarding the specific events in the cell cycle where homeostasis is exerted are still under debate.

The fission yeast *S. pombe* is an interesting case to discuss. This rapidly dividing microorganism is a yeast but uses symmetric division (hence it is sometimes called “fission yeast”), and was the central model system in pioneering studies of the cell cycle. Its size-correction mechanism is the strongest observed in nature, because it can correct size fluctuations in a single cell cycle. Its inter-division size pattern is close to a sizer, but recently the study of mutants with different cell widths has shown that the mechanism that triggers division is based on a surface-area sensor, triggered at a critical cell surface. The molecular effector of this sensing, a protein called Cdr2, has been identified [516]. Curiously, genetic knockout of this protein does not lead to an ablation of size homeostasis. Rather, fission yeast cells fall back to a volume-based mechanism, suggesting that multiple biochemical circuits play a role in the decision to divide.

Finally, since cells of different species and in different conditions use a range of ways to control cell division, for example sizers or adders. An important question is why a particular species would implement one particular strategy. One possibility is that this trait is under selection, and the fitness of individual cells decreases away from the optimal size. In this case sizers would be favored, because they can compensate for deviations in one cell cycle and minimize fluctuations. A second, more likely, possibility is that intrinsic physiological constraints linking cell cycle and growth are important in determining cell division control. For example, it has been argued that in bacteria size control is a result of a cell's attempt to exert a tight control over the initiation of DNA replication – rather than cell division [517].

## 12.7. Concluding remarks

This chapter focused on modeling the cell cycle. The reader should have acquired an overview of some of the key recent experimental results in this area, as well as the basic mathematical toolbox to address biological questions motivated by single-cell dynamic data, concerning (i) decisional processes during the cell cycle and primarily the decision to divide, (ii) coordination between different cell-cycle processes, and primarily the chromosome cycle with cell division and (iii) the coordination of cell cycle progression with growth.

This chapter is connected with Chapters 8 and 9 describing resource allocation models used here to describe growth, to Chapter 11, describing models of growth rate variability, because it provides a framework to include a description of the division rate variability.

## Recommended readings

* Osella M, Tans SJ, Cosentino Lagomarsino M. Step by Step, Cell by Cell: Quantification of the Bacterial Cell Cycle. Trends Microbiol. 2017 Apr;25(4):250-256. doi: 10.1016/j.tim.2016.12.005. Epub 2017 Jan 13. PMID: 28094092.

* Willis L, Huang KC. Sizing up the bacterial cell cycle. Nat Rev Microbiol. 2017 Oct;15(10):606-620. doi: 10.1038/nr-micro.2017.79. Epub 2017 Aug 14. PMID: 28804128.

* Cadart, C., Venkova, L., Recho, P. et al. The physics of cell-size regulation across timescales. Nat. Phys. 15, 993–1004 (2019). https://doi.org/10.1038/s41567-019-0629-y