# Dynamical analysis and optimization of a generalized resource allocation model of microbial growth

Agustín Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé, Hidde de Jong, Francis Mairet

► To cite this version:

Agustín Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé, Hidde de Jong, Francis Mairet. Dynamical analysis and optimization of a generalized resource allocation model of microbial growth. SIAM Journal on Applied Dynamical Systems, 2022, 21 (1), pp.137-165. 10.1137/21M141097X. hal-03251044v2

HAL Id: hal-03251044

https://inria.hal.science/hal-03251044v2

Submitted on 9 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamical analysis and optimization of a generalized resource allocation model of microbial growth*

Agustín Gabriel Yabo†, Jean-Baptiste Caillau‡, Jean-Luc Gouzé†, Hidde de Jong§, and Francis Mairet¶

**Abstract.** Gaining a better comprehension of the growth of microorganisms is a major scientific challenge, which has often been approached from a resource allocation perspective. Simple mathematical self-replicator models based on resource allocation principles have been surprisingly effective in accounting for experimental observations of the growth of microorganisms. Previous work, using a three-variable resource allocation model, predicted an optimal resource allocation scheme for the adaptation of microbial cells to a sudden nutrient change in the environment. We here propose an extended version of this model considering also proteins responsible for basic housekeeping functions, and we study their impact on predicted optimal strategies for resource allocation following changes in the environment. A full dynamical analysis of the system shows there is a single globally attractive equilibrium, which can be related to steady-state growth conditions of bacteria observed in experiments. We then explore the optimal allocation strategies using optimization and optimal control theory. We show that the solutions to this dynamical problem have a complicated structure that includes a second-order singular arc given in feedback form, and characterized by i) Fuller's phenomenon, and ii) the turnpike effect, producing a very particular asymptotic behaviour towards the solution of the static optimization problem. Our work thus provides a generalized perspective on the analysis of microbial growth by means of simple self-replicator models.

**Key words.** systems biology, bacterial growth laws, resource allocation, nutritional shifts, optimal control, turnpike

**AMS subject classifications.** 37N25, 49K15, 92C42

**1. Introduction.** The growth of microorganisms is a paradigm example of self-replication in Nature. Microbial cells are capable of transforming nutrients from the environment into new microbial cells astonishingly fast and in a highly reproducible manner [1]. The biochemical reaction network underlying microbial growth has evolved under the pressure of natural selection, a process that has retained changes in the network structure and dynamics increasing fitness, i.e., favoring the ability of the cells to proliferate in their environment. Gaining a better comprehension of the growth of microorganisms in the context of evolution is a major scientific challenge [2], and the ability to externally control growth is critical for a wide range of applications, such as in combating antibiotics resistance, food preservation, and biofuel

<footnote>*Submitted to the editors 8th April 2021.</footnote>

<footnote>**Funding:** This work was partially supported by ANR project Maximic (ANR-17-CE40-0024-01), Inria IPL Cosy and Labex SIGNALIFE (ANR-11-LABX-0028-01). We acknowledge the support of the FMJH Program PGMO and the support to this program from EDF-THALES-ORANGE.</footnote>

<footnote>†Université Côte d'Azur, Inria, INRAE, CNRS, Sorbonne Université, Biocore Team, Sophia Antipolis, France (agustin.yabo@inria.fr, www.agustinyabo.com.ar)</footnote>

<footnote>‡Université Côte d'Azur, CNRS, Inria, LJAD, France (jean-baptiste.caillau@univ-cotedazur.fr)</footnote>

<footnote>§Université Grenoble Alpes, Inria, 38000 Grenoble, France (hidde.de-jong@inria.fr)</footnote>

<footnote>¶Ifremer, Physiology and Biotechnology of Algae laboratory, rue de l'Ile d'Yeu, 44311 Nantes, France (francis.mairet@ifremer.fr)</footnote>
production [3, 4, 5].

A fruitful perspective on microbial growth is to view it as a resource allocation problem [6]. Microorganisms must assign their available resources to different cellular functions, including the uptake and conversion of nutrients into molecular building blocks of proteins and other macromolecules (metabolism), the synthesis of proteins and other macromolecules from these building blocks (gene expression), and the detection of changes in the environment and the preparation of adequate responses (signalling and regulation). It is often assumed that microorganisms have evolved resource allocation strategies so as to maximize their growth rate, as this would allow them to outgrow competing species.

Simple mathematical models based on resource allocation principles have been surprisingly effective in accounting for experimental observations of the growth and physiology of microorganisms [6, 7, 8, 9, 10, 11, 12]. Instead of providing a detailed description of the entire biochemical reaction network, these models include a limited number of macroreactions responsible for the main growth-related functions of the cell. The models usually take the form of nonlinear ODE systems, typically 3-10 equations with parameters obtained from the experimental literature or estimated from published data. The models have been instrumental in explaining a number of steady-state relations between the growth rate and the cellular composition, in particular the concentration of ribosomes, protein complexes that are responsible for the synthesis of new proteins [8, 6, 13, 10, 14]. Moreover, they have brought out a trade-off between the rate and yield of alternative metabolic pathways that produce energy-carrying molecules, necessary for driving forward many cellular reactions, such as those involved in the synthesis of proteins and other macromolecules [8, 15, 16].

In previous work, using a three-variable resource allocation model, it was possible to predict an optimal resource allocation scheme for the response of microbial cells to a sudden nutrient change in the environment [10]. The prediction was based on the Infinite Horizon Maximum Principle, a generalization of the well-known PMP (Pontrjagin Maximum Principle) [17, 18]. A feedback control strategy inspired by a known regulatory mechanism for growth control in the bacterial cell was shown to give a quasi-optimal approximation of the optimal solution. Strategies for optimal control were also explored for an extension of the model, inspired by recent experimental work [19], which comprises a pathway for the production of a metabolite of biotechnical interest as well as an external signal allowing growth to be switched off [20, 21, 22, 23]. We showed by a combination of analytical and computational means that the optimal solution for the targeted metabolite production problem consists of a phase of growth maximization followed by a phase of product maximization, in agreement with strategies proposed in metabolic engineering. Optimal control approaches have also been used for studying other dynamic optimization problems in biology (see [24] for a review). A classical example is the determination of optimal activation patterns of metabolic pathways, such as to minimize the transition time of metabolites or minimize enzyme costs [25, 26].

The resource allocation model that lies at the basis of the above-mentioned work [10] has a number of limitations. First, the biomass of the cell was assumed to consist of two classes of proteins, enzymes catalyzing metabolic reactions and ribosomes responsible for protein synthesis, whose relative proportions vary with the growth rate. However, experimental data show that a large fraction of the total protein contents of the cell is growth rate-independent [27]. This suggests the introduction of a third protein category, dedicated mainly to ba-
sic housekeeping functions of the cell. The proportion of these proteins is independent of the growth rate and thus constrains the variations in the other two, growth rate-dependent categories [6, 13]. Second, the concentration of ribosomes and enzymes, the two protein categories included in the original model, have both a growth rate-dependent and a growth rate-independent component [6, 27]. This implies that the protein synthesis rate, and thus the growth rate, does not depend on the total ribosome concentration, as in the original model, but only on its growth rate-dependent fraction [13].

In the present manuscript, we revise the above modeling assumptions and study their impact on predicted optimal strategies for resource allocation following changes in the environment of different nature (i.e., changes in the nutrient concentration or stress responses). This leads to a number of interesting problems in mathematical analysis and control, which are addressed using tools from dynamical systems analysis and optimal control theory. A full dynamical analysis of the system shows there is a single globally attractive equilibrium, which can be related to steady-state growth conditions of bacteria observed in experiments. In spite of the simplicity of the presented model, the solutions of the associated biomass maximization problems exhibit quite interesting features. Notably, the second-order singular arc is characterized by a) the Fuller's phenomenon at its junctions, yielding an infinite set of switching points in a finite-time window, and b) the turnpike effect, which produces very particular asymptotic behaviors towards the solution of the static optimization problem. We provide a full description of the singular arc in terms of the state, as well as an explicit proof of the presence of the turnpike effect. While the predicted (optimal) control dynamics does not change much qualitatively in comparison with the previous model, the more realistic modeling assumptions offer a more general perspective of the biological problem. For example, in contrast with the previous model where the absence of growth-rate independent protein yields a constant singular arc equal to the solution of the static optimization problem, the singular arc of the new model is not constant, but governed by a turnpike phenomenon.

In Section 2, we describe the model used in this study, followed by a global dynamical analysis of the model in Section 3. In Section 4, we calibrate the model from literature data using the equilibrium of interest for an optimal steady-state allocation parameter, and in Section 5 we formulate an optimal control problem and prove properties of the optimal solutions. In Section 6, we show that the general analysis can be applied to two different cases of environmental changes related to nutrient shifts and stress responses.

**2. Model definition.** We define a self-replicator system composed of the mass of precursor metabolites P, the gene expression machinery R (ribosomes, RNA polymerase, ...) and the metabolic machinery M (enzymes, transporters, ...), as shown in Figure 1. Essentially, the ribosomal proteins R are responsible for the fabrication of new proteins, and the metabolic proteins M are in charge of the uptake of nutrients for building precursor metabolites P. Following Scott et al. [6], we also introduce a class Q of proteins whose functions fall outside the range of tasks performed by M and R. This sector comprises mainly growth rate-independent proteins such as housekeeping proteins responsible for the maintenance of certain basic cellular functions. Needless to say, the synthesis of Q proteins draws resources away from the pathways to M and R, and consequently imposes an upper bound on the fraction of resources dedicated to self-replication and nutrient uptake. This constraint appears in the model through a constant
<figure><img src="image_2.png" /><figcaption>Figure 1: Coarse-grained self-replicator model. The external substrate S is consumed by bacteria and transformed into precursor metabolites P by the metabolic machinery M. The precursors are used to produce macromolecules of classes R, M and Q, with proportions $\gamma\alpha$, $\gamma(1 - \alpha)$, and $1 - \gamma$, respectively. Solid lines indicate the macroreactions with their respective synthesis rates, and dashed lines denote a catalytic effect.</figcaption></figure>

$\gamma \in [0, 1]$, and it indicates the maximum fraction of the protein synthesis rate available for making ribosomes and metabolic enzymes. The overall allocation process can be represented by the biochemical macroreactions

$$S \xrightarrow{V_M} P,$$

$$P \xrightarrow{V_R} \gamma\alpha R + \gamma(1-\alpha) M + (1-\gamma) Q.$$

The first reaction describes the transformation of external substrate S into precursor metabolites P at a rate $V_M$. The second reaction represents the conversion of precursors into macromolecules R, M, and Q at a rate $V_R$. The roles of the enzymes M in the uptake and metabolization of nutrients and the ribosomal proteins R in the production of proteins are represented through catalytic effects, indicated with dotted arrows in Figure 1. In this context, protein A catalyzes reaction B means that the rate of reaction B is proportional to the cellular concentration of A, but the reaction itself does not consume A. The natural resource allocation strategy is modeled through the time-varying function $\alpha(t) \in [0, 1]$. Thus, the proportion of the total synthesis rate of proteins dedicated to the gene expression machinery R is $\gamma\alpha$, while that of the metabolic machinery M is $\gamma(1 - \alpha)$. In particular, the allocation parameter does not influence the synthesis rate of Q, with constant proportion $1 - \gamma$, as the synthesis of proteins in this class is auto-regulated through mechanisms not relevant in this study. From a biological perspective, the function $\alpha(t)$ represents the naturally-evolved allocation strategy of the cell which is, a priori, unknown. In the context of control theory, and throughout this paper, $\alpha$ is treated as the control input of the system.
**2.1. Self-replicator system.** Generalizing upon Giordano et al. [10], a mass balance analysis yields the dynamical system

$$ 
\left\{
\begin{aligned}
\dot{P} &= V_M - V_R, \\
\dot{R} &= \gamma\alpha V_R, \\
\dot{M} &= \gamma(1 - \alpha)V_R, \\
\dot{Q} &= (1 - \gamma)V_R.
\end{aligned}
\right.
$$

where mass quantities $P$, $M$, $R$ and $Q$ are described in grams (g), the synthesis rates $V_M$ and $V_R$ in grams per hour, and $\alpha$ is the dimensionless allocation parameter. In what follows, we will assume that the proteins of classes $R$, $M$ and $Q$ are responsible for most of the bacterial mass [1], and so we define the bacterial volume $V$ measured in liter units (L) as

$$ (2.1) \qquad V = \beta(R + M + Q), $$

where $\beta$ corresponds to a density constant relating mass and bacterial volume [28], such that the total biomass in grams is given by $V/\beta$. The above assumption implies that the mass of precursor metabolites represents a negligible fraction of the total biomass, (in other words, $P \ll V/\beta$). We define the intracellular concentrations in grams per liter as

$$ (2.2) \qquad p_v \doteq \frac{P}{V}, \quad r_v \doteq \frac{R}{V}, \quad m_v \doteq \frac{M}{V}, \quad q_v \doteq \frac{Q}{V}. $$

Using (2.1) and (2.2), we to obtain the relation

$$ (2.3) \qquad r_v + m_v + q_v = \frac{1}{\beta}. $$

We also define the rates of mass flow per unit volume, which we assume to be functions of the available concentrations, as

$$ v_M(s, m_v) \doteq \frac{V_M}{V}, \quad v_R(p_v, r_v) \doteq \frac{V_R}{V}, $$

where $s$ corresponds to the extracellular concentration of substrate measured in grams per liter. The growth rate of the bacterial population is defined as the relative change of the bacterial volume

$$ \mu \doteq \frac{\dot{V}}{V} = \frac{\beta V_R}{V} = \beta v_R(p_v, r_v). $$
We write the system in terms of the concentrations as

$$
\begin{equation}
\begin{aligned}
\dot{p}_{\nu} &= v_M(s, m_{\nu}) - (1 + \beta p_{\nu}) v_R(p_{\nu}, r_{\nu}), \\
\dot{r}_{\nu} &= (\gamma \alpha - \beta r_{\nu}) v_R(p_{\nu}, r_{\nu}), \\
\dot{m}_{\nu} &= (\gamma(1 - \alpha) - \beta m_{\nu}) v_R(p_{\nu}, r_{\nu}), \\
\dot{q}_{\nu} &= ((1 - \gamma) - \beta q_{\nu}) v_R(p_{\nu}, r_{\nu}), \\
\dot{\nu} &= \beta v_R(p_{\nu}, r_{\nu}) \nu.
\end{aligned}
\tag{2.2}
\end{equation}
$$

2.2. Kinetic definition. We define the kinetics of the reaction system by taking into account that a minimal concentration of ribosomal proteins r<sub>V,min</sub> ∈ (0, γ/β) is required for protein synthesis to take place. In other words, a part of the bacterial volume is occupied by ribosomal proteins which do not directly contribute to growth [13]. Such behavior can be modeled as

$$
v_R(p_\nu, r_\nu) \doteq w_R(p_\nu) (r_\nu - r_{\nu,\min})^+, \quad \text{with } (r_\nu - r_{\min})^+ = \begin{cases} r_\nu - r_{\nu,\min} & \text{if } r_\nu \ge r_{\nu,\min}, \\ 0 & \text{if } r_\nu < r_{\nu,\min}. \end{cases}
$$

Later on, we will see that there is no need to define $v_R(p_ν, r_ν)$ for $r_ν < r_{ν,min}$ if the initial conditions lie in a particular region of the state space. The rate of nutrient uptake is defined as

$$
v_M(s, m_v) \doteq w_M(s) m_v.
$$

We will make the following assumption for functions $w_R(p_ν)$ and $w_M(s)$.

Hypothesis 2.1. Function $w_i(x) : \mathbb{R}_+ \to \mathbb{R}_+$ is

• Continuously differentiable w.r.t. $x$,

• Null at the origin: $w_i(0) = 0$,

189 • Strictly increasing: $w'_{i}(x) > 0, \forall x \ge 0,$

190 • Strictly concave: $w''_{i}(x) < 0, \forall x \ge 0,$

191 • Upper bounded: $\lim_{x \to \infty} w_i(x) = k_i > 0$.

The classical Michaelis-Menten kinetics satisfies Hypothesis 2.1. While most of the math-
ematical results are based on this general definition, for the calibration of the model and
numerical simulations, we will resort to the particular case where the functions are defined as

$$
(2.4) \qquad w_R(p_\nu) \doteq k_R \frac{p_\nu}{K_R + p_\nu}, \quad w_M(s) \doteq k_M \frac{s}{K_S + s},
$$

where $k_R$ and $k_M$ are the maximal reaction rates in h<sup>-1</sup>, and $K_M$ and $K_R$ are the half-saturation constants of the synthesis rates in g L<sup>-1</sup>. For the general case introduced in Hypothesis 2.1 we will define

$$
k_R \doteq \lim_{p_\nu \to \infty} w_R(p_\nu).
$$
**2.3. Constant environmental conditions.** We assume that the availability of the substrate in the medium is constant over the time-window analyzed. The latter can be modeled by setting *s* constant, and thus removing the dynamics of *s* from the system.

**Hypothesis 2.2.** The flow of substrate can be expressed as $w_M(s) = e_M$ with $e_M > 0$ constant.

Using this assumption, the dynamical equation of $p_V$ becomes

$$ \dot{p}_V = e_M m_V - (1 + \beta p_V) w_R(p_V) (r_V - r_{V,\min})^+. $$

The constant $e_M$ models the substrate availability of the medium, but it is also related to the quality of the nutrient and the efficiency of the macroreaction that produces precursor metabolites.

**2.4. Mass fraction formulation and non-dimensionalization.** We define mass fractions of the total bacterial mass as

$$ p \doteq \beta p_V, \quad r \doteq \beta r_V, \quad r_{\min} \doteq \beta r_{V,\min}, \quad m \doteq \beta m_V, \quad q \doteq \beta q_V, $$

which, replacing in (2.3), yields the relation

$$ (2.5) \qquad r + m + q = 1. $$

We also define the non-dimensional time variable $\hat{t} \doteq k_R t$, and the non-dimensional growth rate

$$ (2.6) \qquad \hat{\mu}(p, r) \doteq \frac{\mu(p_V, r_V)}{k_R} = \hat{w}_R(p)(r - r_{\min}), $$

with $\hat{w}_R(p) : \mathbb{R}_+ \to [0, 1)$ defined as $\hat{w}_R(p) \doteq w_R(p_V)/k_R$, and $E_M \doteq e_M/k_R$. For the sake of simplicity, let us drop all hats from the current notation. Then, the model becomes

$$ \left\{ \begin{array}{l} \dot{p} = E_M m - (p+1)w_R(p)(r-r_{\min})^+, \\ \dot{r} = (\gamma\alpha - r)w_R(p)(r-r_{\min})^+, \\ \dot{m} = (\gamma(1-\alpha) - m)w_R(p)(r-r_{\min})^+, \\ \dot{V} = w_R(p)(r-r_{\min})^V, \\ m+r \le 1, \end{array} \right. $$

(S)

where *q* has been removed since it can be expressed in terms of the other concentrations through (2.5); and the constraint $m + r \le 1$ is required to comply with $q \ge 0$. The model differs from that of Giordano et al. by the addition of the category of housekeeping proteins (*q*) and a minimum concentration of ribosomes for protein synthesis (*r*<sub>min</sub>). In what follows, we will systematically investigate how these differences affect the asymptotic behavior and optimal resource allocation strategies.
**3. Asymptotic behavior.** In the present section, we will study the asymptotic behavior of the reduced system representing the intracellular dynamics

$$ (3.1) \qquad \left\{ \begin{array}{l} \dot{p} = E_M m - (p+1)w_R(p)(r - r_{\min})^+, \\ \dot{r} = (\gamma\alpha - r)w_R(p)(r - r_{\min})^+, \\ \dot{m} = (\gamma(1-\alpha) - m)w_R(p)(r - r_{\min})^+, \\ m + r \le 1, \end{array} \right. $$

where $\mathcal{V}$ has been removed since none of the remaining states explicitly depends on it, and it only reaches a steady state when there is no bacterial growth (otherwise, $\dot{\mathcal{V}} > 0$). We will start by stating the invariant set of interest.

**Lemma 3.1.** *The set*

$$ \Gamma = \{(p, r, m) \in \mathbb{R}^3 : p \ge 0, \gamma \ge r \ge r_{\min}, \gamma \ge m \ge 0, m + r \le 1\} $$

*is positively invariant by (3.1).*

*Proof.* This can be easily verified by evaluating the differential equations of system (3.1) over the boundaries of $\Gamma$. As for the condition $m+r \le 1$, we can define a variable $z \doteq m+r$ that obeys the dynamics

$$ \dot{z} = (\gamma - z)w_R(p)(r - r_{\min})^+ $$

which, when evaluated at $z = 1$ yields $\dot{z} \le 0$, as $r_{\max} < 1$, which proves its invariance. ■

This Lemma states that $\gamma \ge r \ge r_{\min}$ for any trajectory with initial conditions in $\Gamma$. As a consequence, there is no need to define the flow $v_R(p, r)$ for values of $r$ under $r_{\min}$. The same thing can be said for the constraint $m+r \le 1$, which is valid for every trajectory starting in $\Gamma$. Additionally, since $\gamma$ represents the maximal ribosomal mass fraction, we will define the following parameter.

**Definition 3.2.** *The maximal ribosomal mass fraction is $r_{\max} \doteq \gamma$.*

Then, we will reduce the study of the system to this set and so, using Definition 3.2, we redefine (3.1) as

$$ (S') \qquad \left\{ \begin{array}{l} \dot{p} = E_M m - (p+1)w_R(p)(r - r_{\min}) \\ \dot{r} = (r_{\max}\alpha - r)w_R(p)(r - r_{\min}) \\ \dot{m} = (r_{\max}(1-\alpha) - m)w_R(p)(r - r_{\min}) \end{array} \right. $$

where $(r-r_{\min})^+$ has been replaced by $r-r_{\min}$, and the constraint $m+r \le 1$ has been removed. Furthermore, we will define the minimum constant allocation parameter $\alpha_{\min}^*$ necessary to allow steady-state self-replication, given by

$$ \alpha_{\min}^* \doteq \frac{r_{\min}}{r_{\max}}. $$

Its importance will be analyzed throughout the current section.
## 3.1. Local stability.

Theorem 3.3. System (S') has the equilibria

* $E_1 \doteq (p^*, r^*, m^*)$, locally stable if $\alpha^* > \alpha_{\min}^*$.

* $E_2 \doteq (p, r_{\min}, 0)$, locally unstable if $\alpha^* > \alpha_{\min}^*$.

* $E_3 \doteq (0, r, 0)$, locally unstable if $r \neq r_{\min}$.

with

$$ (3.2) \qquad \begin{aligned} p^* &\doteq \{p \in \mathbb{R}_+ : (p+1)w_R(p) = \frac{E_M m^*}{r^* - r_{\min}}\}, \\ r^* &\doteq r_{\max} \alpha^*, \\ m^* &\doteq r_{\max}(1 - \alpha^*). \end{aligned} $$

*Proof.* The general Jacobian matrix of the system (S') is

$$ (3.3) \qquad \begin{bmatrix} -(w_R(p) + (p+1)w'_R(p))(r-r_{\min}) & -(p+1)w_R(p) & E_M \\ (r_{\max}\alpha - r)w'_R(p)(r-r_{\min}) & (r_{\max}\alpha - 2r + r_{\min})w_R(p) & 0 \\ (r_{\max}(1-\alpha) - m)w'_R(p)(r-r_{\min}) & (r_{\max}(1-\alpha) - m)w_R(p) & -w_R(p)(r-r_{\min}) \end{bmatrix}. $$

We first see that, if $\alpha^* > \alpha_{\min}^*$, the value $p^*$ is unique since $(p+1)w_R(p)$ is a monotone increasing function satisfying $w_R(0) = 0$ and $\lim_{p \to \infty} (p+1)w_R(p) = \infty$ (as stated in Hypothesis 2.1), and $E_M m^* / (r^* - r_{\min}) > 0$, so the set (3.2) yields a unique solution. For $\alpha^* < \alpha_{\min}^*$, the equation for $p^*$ in (3.2) has no valid solution as $E_M m^* / (r^* - r_{\min})$ becomes negative, and therefore the equilibrium does not exist. The Jacobian (3.3) for $E_1$ becomes

$$ J_1 = \begin{bmatrix} -(w_R(p^*) + (p^* + 1)w'_R(p^*)) (r^* - r_{\min}) & -(p^* + 1)w_R(p^*) & E_M \\ 0 & -(r^* - r_{\min})w_R(p^*) & 0 \\ 0 & 0 & -w_R(p^*)(r^* - r_{\min}) \end{bmatrix} $$

and so the local stability of the equilibrium is given by the signs of the roots of the characteristic polynomial, which are $\lambda = -(w_R(p^*) + (p^*+1)w'_R(p^*)) (r^* - r_{\min})$, $\lambda = -(p^*+1)w_R(p^*)$, and $\lambda = -w_R(p^*)(r^* - r_{\min})$. As the three roots are negative, we conclude that, if the equilibrium exists, it is locally stable. For the second equilibrium $E_2$, the Jacobian is

$$ J_2 = \begin{bmatrix} 0 & -w_R(p) & E_M \\ 0 & (r^* - r_{\min})w_R(p) & 0 \\ 0 & r_{\max}(1 - \alpha^*)w_R(p) & 0 \end{bmatrix} $$

with characteristic polynomial

$$ P_2(\lambda) = \lambda^2(\lambda - (r^* - r_{\min})w_R(p)). $$
296 If $\alpha^* > \alpha_{\min}^*$, then $J_2$ has one positive eigenvalue and $E_2$ becomes locally unstable. As for $E_3$, 297 the Jacobian is

$$J_3 = \begin{bmatrix} -w'_R(0)(r - r_{\min}) & 0 & E_M \\ (r_{\max}\alpha - r)w'_R(0)(r - r_{\min}) & 0 & 0 \\ r_{\max}(1 - \alpha)w'_R(0)(r - r_{\min}) & 0 & 0 \end{bmatrix}$$

300 with characteristic polynomial

$$P_3(\lambda) = \lambda^2(\lambda + w'_R(0)(r - r_{\min})) - E_M r_{\max}(1 - \alpha)w'_R(0)(r - r_{\min})\lambda.$$

303 One root is $\lambda = 0$, and the two remaining roots can be found by solving the equation

$$\lambda^2 + \lambda w'_R(0)(r - r_{\min}) - E_M r_{\max}(1 - \alpha)w'_R(0)(r - r_{\min}) = 0.$$

306 By the Routh-Hurwitz criterion, the two remaining roots are in the open left half plane if and only if $w'_R(0)(r - r_{\min}) > 0$ and $E_M r_{\max}(1 - \alpha)w'_R(0)(r - r_{\min}) < 0$, which is never true. Consequently, for $r \neq r_{\min}$, there is at least one positive root, and so the equilibrium is 308 unstable. 309 ■

310 **3.2. Global behavior.** We will study the global behavior of system (S') for the initial 311 conditions

$$(IC) \qquad p(0) > 0, \quad r(0) \in (r_{\min}, r_{\max}), \quad m(0) \in (0, r_{\max}), \quad r(0) + m(0) \le 1.$$

314 and for a given constant allocation parameter

$$\alpha(t) = \alpha^* \in (\alpha_{\min}^*, 1).$$

317 Under this constraint, we see that the dynamics of $r$ and $m$ become

$$\dot{r} = (r^* - r)w_R(p)(r - r_{\min}), \quad \dot{m} = (m^* - m)w_R(p)(r - r_{\min}),$$

320 which means that, if $p > 0$ and $r > r_{\min}$, the signs of $\dot{r}$ and $\dot{m}$ are given by the signs of $r^* - r$ 321 and $m^* - m$, respectively (and both $\dot{r}$ and $\dot{m}$ are zero if $p = 0$ or $r = r_{\min}$). Then, let us 322 divide $\Gamma$ into the subsets

$$\begin{align*} \mathcal{R}^{-} &\doteq \{(p, r, m) \in \Gamma : r \in (r_{\min}, r^{*})\}, & \mathcal{M}^{-} &\doteq \{(p, r, m) \in \Gamma : m \in (0, m^{*})\}, \\ \mathcal{R}^{+} &\doteq \{(p, r, m) \in \Gamma : r \in (r^{*}, r_{\max})\}, & \mathcal{M}^{+} &\doteq \{(p, r, m) \in \Gamma : m \in (m^{*}, r_{\max})\}, \end{align*}$$

325 such that $\Gamma = \bar{\mathcal{R}}^{-} \cup \bar{\mathcal{R}}^{+} = \bar{\mathcal{M}}^{-} \cup \bar{\mathcal{M}}^{+}$. In these sets, the following holds.

326 **Lemma 3.4.** For $\alpha(t) = \alpha^* \in (\alpha_{\min}^*, 1)$, the closed sets $\bar{\mathcal{R}}^{-}$, $\bar{\mathcal{R}}^{+}$, $\bar{\mathcal{M}}^{-}$ and $\bar{\mathcal{M}}^{+}$ are invariant 327 by (S'), and

$$\left\{ \begin{array}{ll} \dot{r} \ge 0 & \text{if } (p, r, m) \in \mathcal{R}^{-}, \\ \dot{r} \le 0 & \text{if } (p, r, m) \in \mathcal{R}^{+}, \end{array} \right. \qquad \left\{ \begin{array}{ll} \dot{m} \ge 0 & \text{if } (p, r, m) \in \mathcal{M}^{-}, \\ \dot{m} \le 0 & \text{if } (p, r, m) \in \mathcal{M}^{+}. \end{array} \right.$$

330 Again, the invariance of the sets can be checked by evaluating the vector field over the 331 boundaries of the sets. Now we state a first result.
**Proposition 3.5.** For $\alpha(t) = \alpha^* \in (\alpha_{\min}^*, 1)$ and initial conditions (IC), system (S') has a lower bound

$$ (p, r, m) \geq (p_{\text{low}}, r_{\text{low}}, m_{\text{low}}) \text{ for all } t \geq 0, $$

with

$$ \begin{gathered} r_{\text{low}} \doteq \min(r(0), r^*), \quad m_{\text{low}} \doteq \min(m(0), m^*), \\ p_{\text{low}} \doteq \left\{ p \in \mathbb{R}_+ : (p+1)w_R(p) = \frac{E_M m_{\text{low}}}{r_{\max} - r_{\min}} \right\}. \end{gathered} \tag{3.4} $$

*Proof.* For a trajectory emanating from $\mathcal{R}^-$ (respectively, $\mathcal{R}^+$), it follows that $\dot{r} \geq 0$ (respectively, $\dot{r} \leq 0$) for all $t$ (according to Lemma 3.4), and so $r \geq r(0)$ (respectively, $r \geq r^*$) for all $t$. This proves that $r \geq \min(r(0), r^*) > r_{\min}$ for all $t$ (depending on whether the trajectory starts in $\mathcal{R}^-$ or $\mathcal{R}^+$). Similarly, a trajectory starting in $\mathcal{M}^-$ (respectively, $\mathcal{M}^+$) meets $\dot{m} \geq 0$ (respectively, $\dot{m} \leq 0$) for all $t$, and so $m \geq m(0)$ (respectively, $m \geq m^*$) for all $t$. Then, it follows that $m \geq \min(m(0), m^*)$ for all $t \geq 0$. The equation for $p$ can thus be lower-bounded to

$$ \dot{p} \geq E_M m_{\text{low}} - (p+1)w_R(p)(r_{\max} - r_{\min}), $$

which means $p \geq p_{\text{low}}$ for all $t \geq 0$, with $p_{\text{low}}$ the solution of (3.4), which is unique by the same arguments as those used in Theorem 3.3. ■

A lower bound on system (S') is a stronger condition than the classical persistence for biological populations, as the bound is imposed not only for $t \to \infty$ but for the whole trajectory. As a consequence, the growth rate never vanishes, as it meets $\mu(p, r) \geq w_R(p_{\text{low}})(r_{\text{low}} - r_{\min}) > 0$ for all $t \geq 0$. Then, the global stability of the system is straightforward.

**Theorem 3.6.** For $\alpha(t) = \alpha^* \in (\alpha_{\min}^*, 1)$ and initial conditions (IC), every solution of (S') converges to the equilibrium $E_1$.

*Proof.* Since $p \geq p_{\text{low}} > 0$ and $r \geq r_{\text{low}} > r_{\min}$ for all $t \geq 0$, we have that $\text{sign}(\dot{r}) = \text{sign}(r^* - r)$ and $\text{sign}(\dot{m}) = \text{sign}(m^* - m)$, showing that $r$ and $m$ converge asymptotically to $r^*$ and $m^*$, respectively. Consequently, the dynamical equation of $p$ becomes $\dot{p} = E_M m^* - (p + 1)w_R(p)(r^* - r_{\min})$ and so $\text{sign}(\dot{p}) = \text{sign}(p^* - p)$, which means that $p$ converges asymptotically to the steady-state value $p^*$. ■

**Remark 3.7.** For the case over the invariant plane given by $r(0) = r_{\min}$ and $m(0) > 0$, concentrations $m$ and $r$ are constant along the whole trajectory, and $p$ increases linearly with time (as $\dot{p} = E_M m(0)$). This is a degenerate case that contradicts the assumption $p \ll 1$, and lacks biological relevance.

**3.3. Maximum steady-state growth rate.** A classical hypothesis in the literature is to suppose bacterial populations in steady-state regimes maximize their growth rate ([10] and references therein). We are interested in finding the static allocation strategy $\alpha^*$ that produces this situation. Since the only equilibrium that admits bacterial growth is $E_1$, we will express
369 the static optimization problem as

$$370 \quad \max_{\alpha^* \in [\alpha_{\min}^*, 1]} \mu(p^*, r^*), \\ 371$$

372 which can be rewritten as $\mu(p^*, r^*) = w_R(p^*)(r^* - r_{\min})$. It is possible to express $\alpha^*$ in terms of $p^*$ through the relation

$$374 \quad (3.5) \qquad \alpha^*(p^*) = \frac{E_M + (p^* + 1)w_R(p^*)\alpha_{\min}^*}{E_M + (p^* + 1)w_R(p^*)}.$$

376 Moreover, since the above function $\alpha^*(p^*) : \mathbb{R}_+ \to (\alpha_{\min}^*, 1]$ is monotone decreasing, it is possible to write the optimization problem in terms of $p^*$ instead of $\alpha^*$. The growth rate in terms of $p^*$ writes

$$379 \quad (3.6) \qquad w_R(p^*)(r^* - r_{\min}) = (r_{\max} - r_{\min}) \left( \frac{E_M w_R(p^*)}{E_M + (p^* + 1)w_R(p^*)} \right).$$

381 We differentiate w.r.t. $p^*$ and we get the relation $w_R(p^*)^2 = E_M w_R'(p^*)$, which has a unique solution since, according to Hypothesis 2.1, $w_R(p)^2$ is a monotone increasing function satisfying $w_R^2(0) = 0$ and $\lim_{p \to \infty} w_R^2(p) = 1$, and $w_R'(p)$ is a monotone decreasing function satisfying $w_R'(0) > 0$ and $\lim_{p \to \infty} w_R'(p) = 0$ (as $w_R(p)$ is a strictly increasing upper-bounded function). Then, the condition for optimality can be expressed as

$$386 \quad (3.7) \qquad \frac{w_R(p_{\text{opt}}^*)^2}{E_M w_R'(p_{\text{opt}}^*)} = 1.$$

388 Thus, the optimal allocation parameter $\alpha^*$ is obtained by replacing $p_{\text{opt}}^*$ in (3.5), and the maximal static growth rate can be calculated using (3.6). From (3.7), it can be seen that $p_{\text{opt}}^*$ depends neither on $r_{\min}$ nor on $r_{\max}$, suggesting that the steady-state precursor concentration is independent of the housekeeping protein fraction $q$ and of the growth rate-independent ribosomal fraction. Conversely, the precursor concentration is rather determined by the environmental conditions and by the nature of the function $w_R(p)$. It can be proven that the latter result is not a consequence of assumption (2.1): when considering a definition of the bacterial volume as $\beta(P + R + M + Q)$, which takes into account the mass $P$, the optimal precursor concentration amounts to $p_{\text{opt}}^*/(1 + p_{\text{opt}}^*)$.

397 In addition, from $\dot{p} = 0$ in (S'), we get:

$$398 \quad \frac{r^* - r_{\min}}{m^*} = \frac{E_M}{(p^* + 1)w_R(p^*)}.$$

400 This shows that, for the optimal steady state, the concentration ratio of the active gene expression machinery over the metabolic machinery does not depend on $r_{\max}$ either. Thus, a cellular strategy regulating the precursor concentration and the balance between gene expression and metabolism could lead to the optimal equilibrium, regardless of the demand for $Q$.
**4. Model calibration.** Whereas the parameter values do not affect the results above and the optimal control analysis in the next section, they are nevertheless important for simulations illustrating the dynamics and optimal allocation strategies of system (S'). Below, we derive such parameters for the model bacterium *Escherichia coli*, using published sources. The $\beta$ constant used in the definition of the bacterial volume (2.1) corresponds to the inverse of the protein density, which is set to 0.003 [L g$^{-1}$] based on [10]. According to [6], the ribosomal fraction of the proteome<fnref n="1" /> can vary between 6% and 55%. In more recent studies [27], this sector is divided into growth-rate dependent and independent fractions. The maximal growth-rate dependent ribosomal fraction of the proteome is estimated to be 41%, and the growth rate-independent fraction is 9%. Based on these experimental estimations, we set $r_{max} = 0.41 + 0.09 = 0.5$. We performed further calibrations using data sets from [29, 30, 6, 27] containing measurements of various strains of *E. coli* growing in different media. The data sets are composed essentially of data points (*growth rate*, *RNA/protein mass ratio*) measured at steady state. Most RNA is ribosomal RNA found overwhelmingly in ribosomes, the main constituent of the gene expression machinery. In order to adjust the measurements to model (S'), the observed RNA/protein ratios can be converted to mass fractions $r$ through multiplication with a conversion factor $\rho = 0.76 \mu g$ of protein/$\mu g$ of RNA [6]. As a result, we have $n$ measurements of form $(\tilde{\mu}_k, \tilde{r}_k)$ which are assumed to follow a linear relation [6], as seen in Figure 2a. From the vertical intercept of the linear regression performed using the data points, we obtain $r_{min} = 0.07$, in agreement with previous studies [6, 13, 27]. Each data point, composed of an observed growth rate and its associated ribosomal mass fraction, can be related to an optimal steady state of system (S') for a certain environmental condition $e_M$. Thus, each $k$th pair $(\tilde{\mu}_k, \tilde{r}_k)$ of the $n$ measurements should yield a constant environmental condition $e_{M,k}$, and all pairs should simultaneously adjust the rate constant $k_R$. Such fitting can be done by resorting to the Michaelis-Menten kinetic form introduced in (2.4). Based on [10], we fix the half-saturation constant of protein synthesis $K_R = 1 \text{ g L}^{-1}$. We then define the parameter vector $\theta = (k_R, e_{M,1}, \dots, e_{M,n})$ which is computed by solving a least-squares regression problem. Using the relation (2.6), the cost function to minimize is

$$ \min_{\theta \in \mathbb{R}_{+}^{n+1}} \sum_{k=1}^{n} (\tilde{\mu}_k - \mu_{\text{opt}}^*(k_R, e_{M,k}))^2 + (\tilde{r}_k - r_{\text{opt}}^*(k_R, e_{M,k}))^2, $$

where the non-dimensional growth rate $\mu_{\text{opt}}^*$ is calculated using (3.6), and the optimal steady state $(r_{\text{opt}}^*, p_{\text{opt}}^*, m_{\text{opt}}^*)$ is expressed in terms of $\alpha_{\text{opt}}^*$ (using Theorem 3.3) which is, at the same time, a function of $k_R$ and $e_{M,k}$. The numerical solution yields $k_R = 6.23 \text{ h}^{-1}$, and different values of $e_M$ matching different nutrients from the dataset (see Figure 2b). We can validate these results by computing the maximal growth rate $k_R(r_{max} - r_{min}) = 2.68 \text{ h}^{-1}$ based on the adjusted parameters, which is a value that corresponds well with literature values of the maximal growth rate of *E. coli* in rich media [30].

**5. Optimal resource allocation.**

**5.1. Problem definition.** In this section we formulate the dynamic optimization problem under the hypothesis that microbial populations have evolved resource allocation strategies

<footnote n="1">¹The proteome is the total amount of protein in the cell</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Experimental data from [6, 29, 30] plotted in (a) shows a linearity of $r^2 = 0.9739$ (dashed line, fitted to data) with a vertical intersect $r_{min} = 0.07$ and slope $k_R = 6.23 \text{ h}^{-1}$. In (b), steady-state growth rate curves $\mu^*$ are shown in terms of the mass fraction $r^* \in (r_{min}, r_{max})$ for different fitted values of $e_M$. Each optimal pair $(\mu_{opt}^*, r_{opt}^*)$ marked with color circles corresponds to a sample from the data set of Scott et al. denoted in (a) with circles of matching colors.</figcaption></figure>

enabling them to maximize their biomass [31, 32]. This is represented by an OCP (Optimal Control Problem) where the objective is to maximize the final volume at time $T$ given by $V(T)$. For the sake of convenience, we propose to maximize the quantity $\log V(T)$ (since $\log$ is an increasing function) given by

$$ \log V(T) = \int_{0}^{T} \mu(p, r) dt + \log V(0). $$

As the initial condition $V(0)$ is fixed, we define the cost function

$$ J(u) \doteq \int_{0}^{T} \mu(p, r) dt = \int_{0}^{T} w_{R}(p)(r - r_{\min}) dt. $$

Since $V$ appears neither in the dynamics nor in the cost function, the optimal problem will be written considering the reduced system introduced in (S') with initial conditions given by
(IC). We write the optimal control problem

$$
\begin{array}{ll}
\text{(OCP)} & \left\{
\begin{array}{ll}
\text{maximize} & J(u) = \int_{0}^{T} w_{R}(p)(r - r_{\min}) dt \\
\text{subject to} & \text{dynamics (S'),} \\
& \text{initial conditions (IC),} \\
& \alpha(\cdot) \in \mathcal{U},
\end{array}
\right. \\
&
\end{array}
$$

with $\mathcal{U}$ the set of admissible controllers, which are Lebesgue measurable real-valued functions defined on the interval $[0, T]$ and satisfying the constraint $\alpha(t) \in [0, 1]$. Problem (OCP) has neither final state constraints nor path constraints. In the context of dynamical optimization, the use of path constraints can be useful to restrict the solutions to those meeting certain physical and biological limitations, especially when dealing with more complex models. While enforcing additional constraints on the OCP increases the dimension of the problem, standard optimal control solvers are able to handle such formulations. In this work, imposing initial conditions (IC) guarantees that every trajectory of the system stays within the set $\Gamma$ defined in Lemma 3.1, which ensures that the solutions are consistent with the biological assumptions. In principle, this formulation of the problem resembles the optimal control problem proposed in [10]: the objective is to maximize the accumulation of a certain quantity within the system during a fixed time interval $[0, T]$. The main difference lies in the dynamics of the system, as the introduction of the protein Q increases the system dimension by one, which yields a more relevant (and more complex) associated OCP. We will see in following sections that the problem raised in this work can be solved by a generalization of Giordano et al.'s approach.

**5.2. Pontrjagin Maximum Principle.** Existence of a solution for this class of OCPs is rather trivial. Given that there are no terminal constraints, there is no controllability issue. Moreover, the dynamics is affine in the control with the latter included in a compact and convex set (a closed interval), and one can easily check that every finite-time trajectory remains bounded. So existence is guaranteed by Filippov's theorem [33]. Then, for an optimal control problem (OCP) with state $\varphi \in \mathbb{R}^n$, Pontrjagin maximum principle (PMP) ensures that there exist $\lambda^0 \le 0$ and a piecewise absolutely continuous mapping $\lambda(\cdot) : [0, T] \to \mathbb{R}^n$, with $(\lambda(\cdot), \lambda^0) \ne (0, 0)$, such that the extremal $(\varphi, \lambda, \lambda^0, \alpha)$ satisfies the generalized Hamiltonian system

$$
\begin{array}{ll}
\text{(PMP)} & \left\{
\begin{array}{l}
\dot{\varphi} = \frac{\partial}{\partial \lambda} H(\varphi, \lambda, \lambda^0, \alpha), \\
\dot{\lambda} = -\frac{\partial}{\partial \varphi} H(\varphi, \lambda, \lambda^0, \alpha), \\
H(\varphi, \lambda, \lambda^0, \alpha) = \max_{\alpha \in [0,1]} H(\varphi, \lambda, \lambda^0, \alpha),
\end{array}
\right.
\end{array}
$$

for almost every $t \in [0, T]$. For our particular case, we have the state vector $\varphi \doteq (p, r, m)$ and
adjoint vector $\lambda \doteq (\lambda_p, \lambda_r, \lambda_m)$ and the Hamiltonian given by

$$
(5.1) \qquad H(\varphi, \lambda, \lambda^0, \alpha) = \lambda^0 w_R(p)(r - r_{\min}) + \langle \lambda, F(\varphi, u) \rangle,
$$
where $F$ represents the right-hand side of system (S'). Given that in (OCP) there is no terminal condition on the state $\varphi(T)$, the transversality condition for the adjoint state is $\lambda(T) = 0$, and we can discard abnormal extremals from the analysis. In other words, any extremal $(\varphi, \lambda, \lambda^0, \alpha)$ satisfying PMP is normal, so $\lambda^0 \neq 0$. Developing (5.1) yields the Hamiltonian

$$H = (E_M m - (p+1)w_R(p)(r-r_{\min}))\lambda_p + (r_{\max}\alpha - r)w_R(p)(r-r_{\min})\lambda_r \\ + (r_{\max}(1-\alpha) - m)w_R(p)(r-r_{\min})\lambda_m - \lambda^0 w_R(p)(r-r_{\min}),$$

and the adjoint system is

$$ (5.2) \qquad \left\{ \begin{aligned} \dot{\lambda}_p &= w_R(p)(r-r_{\min})\lambda_p + (p+1)w'_R(p)(r-r_{\min})\lambda_p - (r_{\max}\alpha - r)w'_R(p)(r-r_{\min})\lambda_r \\ &\quad -(r_{\max}(1-\alpha) - m)w'_R(p)(r-r_{\min})\lambda_m + \lambda^0 w'_R(p)(r-r_{\min}), \\ \dot{\lambda}_r &= (p+1)w_R(p)\lambda_p + w_R(p)(r-r_{\min})\lambda_r - (r_{\max}\alpha - r)w_R(p)\lambda_r \\ &\quad -(r_{\max}(1-\alpha) - m)w_R(p)\lambda_m + \lambda^0 w_R(p), \\ \dot{\lambda}_m &= -E_M\lambda_p + w_R(p)(r-r_{\min})\lambda_m. \end{aligned} \right. $$

Since the Hamiltonian is linear in the control $\alpha$, we rewrite it in the input-affine form $H = H_0 + \alpha H_1$ with

$$ H_0 = (E_M m - (p+1)w_R(p)(r-r_{\min}))\lambda_p - rw_R(p)(r-r_{\min})\lambda_r \\ + (r_{\max}-m)w_R(p)(r-r_{\min})\lambda_m - \lambda^0 w_R(p)(r-r_{\min}), $$

$$ (5.3) \qquad H_1 = r_{\max}w_R(p)(r-r_{\min})(\lambda_r - \lambda_m). $$

The constrained optimal control $\alpha$ should maximize the Hamiltonian, so the solution is

$$ \alpha(t) = \begin{cases} 0 & \text{if } H_1 < 0, \\ 1 & \text{if } H_1 > 0, \\ \alpha_{\text{sing}}(t) & \text{if } H_1 = 0, \end{cases} $$

where $\alpha_{\text{sing}}(t)$ is called a singular control, showing that any optimal control is a concatenation of bangs $(\alpha = \pm 1)$ and singular arcs, depending on the sign of the switching function $H_1$. As obtained in [21, 23], a bang arc $\alpha = 0$ (respectively $\alpha = 1$) corresponds to a pure allocation strategy where the production of $R$ (respectively $M$) is completely switched off. While a full description of the optimal control is often difficult to obtain through PMP, there are certain analyses that can be performed to help understand its structure. We will first see that the final bang of the optimal control is an upper bang $\alpha = 1$.

**Lemma 5.1.** There exists $\epsilon$ such that the optimal control solution of (OCP) is $\alpha(t) = 1$ for the interval of time $[T - \epsilon, T]$.
Proof. We define $\lambda_z = \lambda_r - \lambda_m$, where its dynamics can be obtained from (5.2). It can be seen that, when evaluating its dynamics at final time, we get

$$ \dot{\lambda}_z(T) = \lambda^0 w_R(p(T)) < 0, $$

due to the whole adjoint state being null at final time except for $\lambda^0$. As $\lambda_z(T)$ also vanishes due to the transversality conditions, we have $\lambda_z(T - \epsilon) > 0$ for a certain $\epsilon$. Then, $H_1 > 0$ for the interval $[T - \epsilon, T]$, which corresponds to a bang arc $\alpha = 1$. ■

A control $\alpha = 1$ implies a strategy in which all resources are allocated to ribosome synthesis, thus favoring the synthesis of proteins. An intuitive interpretation of Lemma 5.1 is that, when approaching the final time $T$, the most efficient strategy is to exploit as much as possible the available precursors. This is achieved by maximizing the proteins catalyzing $v_R$, at the expense of arresting the uptake of nutrients $v_M$ from the environment. In order to further describe the optimal control, we can analyze the singular extremals. A singular arc occurs when the switching function $H_1$ vanishes over a subinterval of time. A detailed description of the singular arcs can be done by differentiating successively the switching function $H_1$ until the singular control $\alpha_{\text{sing}}$ can be obtained as a function of the state $\varphi$ and the adjoint state $\lambda$.

### 5.3. Study of the singular arcs.

**5.3.1. Introduction.** We assume $H_1$ vanishes on a whole sub-interval $[t_1, t_2] \subset [0, T]$, so the extremal belongs to the singular surface

$$ \Sigma \doteq \{(\varphi, \lambda) \in \mathbb{R}^6 : H_1(\varphi, \lambda) = 0\}. $$

Since $H_1$ vanishes identically, so does its derivative with respect to time. Differentiating along an extremal $(\varphi, \lambda)$ amounts to taking a Poisson bracket<fnref n="2" /> with the Hamiltonian $H$ [33]. Indeed, along the singular arc,

$$ 0 = \dot{H}_1 = \frac{\partial H_1}{\partial \varphi} \dot{\varphi} + \frac{\partial H_1}{\partial \lambda} \dot{\lambda} = \sum_{i=1}^{n} \left( \frac{\partial H}{\partial \lambda_i} \frac{\partial H_1}{\partial \varphi_i} - \frac{\partial H}{\partial \varphi_i} \frac{\partial H_1}{\partial \lambda_i} \right) = \{H, H_1\} = \{H_0, H_1\}. $$

The first derivative $\dot{H}_1 = H_{01} \doteq \{H_0, H_1\}$ is equal to $\langle \lambda, F_{01} \rangle$, where $F_{01}$ corresponds to the Lie bracket of the vector fields $F_0$ and $F_1$. Differentiating again we obtain

$$ 0 = \dot{H}_{01} = H_{001} + \alpha H_{101}. $$

Again, $H_{001} \doteq \langle \lambda, F_{001} \rangle$ where, with the same notation as before, $F_{001}$ is the Lie bracket of $F_0$ with $F_{01}$. If, on the set

$$ \Sigma' \doteq \{(\varphi, \lambda) \in \mathbb{R}^6 : H_1(\varphi, \lambda) = H_{01}(\varphi, \lambda) = 0\}, $$

$$ \{f, g\} = \sum_{i=1}^{n} \left( \frac{\partial f}{\partial \lambda_i} \frac{\partial g}{\partial \varphi_i} - \frac{\partial f}{\partial \varphi_i} \frac{\partial g}{\partial \lambda_i} \right). $$

<footnote n="2">²The Poisson bracket $\{f, g\}$ of two functions $f$ and $g$ along an extremal $(\varphi, \lambda)$ is defined as</footnote>
the bracket $H_{101}$ is also zero, the control disappears from the previous equality, and one has to differentiate at least two more times to retrieve the control: $H_{0001}$ is also zero, and

$$ (5.4) \qquad 0 = H_{00001} + \alpha H_{10001}. $$

If the length-five bracket $H_{10001}$ is not zero, the singular arc is of order two. When $H_{101}$ vanishes not only on $\Sigma'$ but on all $\mathbb{R}^6$, the order is said to be intrinsic and connections between bang and singular arcs can only occur through an infinite number of switchings [34], the so-called Fuller phenomenon. Otherwise, the order is termed local, and Fuller phenomenon may or may not occur. Using (5.4), the singular control $u_s$ is obtained as a function of both the state $\varphi$ and the adjoint state $\lambda$ as

$$ \alpha_s(\varphi, \lambda) \doteq -\frac{H_{00001}}{H_{10001}}. $$

In our low-dimensional situation, there exists the possibility that the singular control is in feedback form, that is, as a function of the state only. The latter can be verified by rewriting the system in dimension four (Mayer optimal control formulation where the final volume is maximized), in terms of $\tilde{\varphi} \doteq (p, r, m, \mathcal{V})$ and its adjoint $\tilde{\lambda} \doteq (\lambda_p, \lambda_r, \lambda_m, \lambda_\mathcal{V})$. The dynamics is affine in the control,

$$ \dot{\tilde{\varphi}} = \tilde{F}_0(\tilde{\varphi}) + \alpha \tilde{F}_1(\tilde{\varphi}), $$

and so is the Hamiltonian:

$$ \tilde{H}(\tilde{\varphi}, \tilde{\lambda}, \alpha) = \tilde{H}_0 + \alpha \tilde{H}_1, $$

with $\tilde{H}_i = \langle \tilde{\lambda}_i, \tilde{F}_i \rangle$, $i = 0, 1$. The same computation as before leads to the following relations along a singular arc of order two:

$$ 0 = \dot{\tilde{H}}_1 = \tilde{H}_{01} = \tilde{H}_{001} = \tilde{H}_{0001}, $$

and

$$ 0 = \tilde{H}_{00001} + \alpha \tilde{H}_{10001}. $$

**Proposition 5.2.** Assume that, for all $\varphi$, $\tilde{F}_1$, $\tilde{F}_{01}$, and $\tilde{F}_{001}$ are independent. Then, an order two singular control depends only on the state $\tilde{\varphi}$, and can be expressed as

$$ \alpha_s(\tilde{\varphi}) = -\frac{\det(\tilde{F}_1, \tilde{F}_{01}, \tilde{F}_{001}, \tilde{F}_{0001})}{\det(\tilde{F}_1, \tilde{F}_{01}, \tilde{F}_{001}, \tilde{F}_{10001})}. $$

*Proof.* The previous relations imply that $\tilde{\lambda}$ is orthogonal to $\tilde{F}_1$, $\tilde{F}_{01}$, $\tilde{F}_{001}$, and also to $\tilde{F}_{00001} + \alpha \tilde{F}_{10001}$. If these four vector fields were independent at some point along the singular
arc, $\tilde{\lambda} \in \mathbb{R}^4$ would vanish: for a problem in Mayer form, this would contradict the maximum principle. So their determinant must vanish everywhere along the arc and

$$ \det (\tilde{F}_1, \tilde{F}_{01}, \tilde{F}_{001}, \tilde{F}_{00001}) + \alpha \det (\tilde{F}_1, \tilde{F}_{01}, \tilde{F}_{001}, \tilde{F}_{10001}) = 0. $$

If the second determinant was zero, given the rank assumption on the first three vector fields, $F_{10001}$ would belong to their span; but this is impossible since it would imply $H_{10001} = 0$, contradicting the fact that the singular is of order two. ■

Going back to the three-dimensional formulation, one can explicit the computations by successively differentiating the expression (5.3).

**5.3.2. Singular arc in feedback form.** The condition $H_1 = 0$ could be a consequence of the growth rate $w_R(p)(r - r_{\min})$ vanishing over the whole interval $[t_1, t_2]$. We will see this is not possible given the dynamics of the system.

**Proposition 5.3.** *The growth rate $\mu(p, r) = w_R(p)(r - r_{\min})$ cannot vanish along the optimal solution of (OCP).*

*Proof.* For any trajectory of (S') with initial conditions (IC), control $\alpha(\cdot) \in \mathcal{U}$ and $t \in [0, T]$, we have $\dot{p} \le E_M r_{\max}$, which means $p \le p_{\max}^T = E_M r_{\max} T + p(0)$. Then, $\dot{r} \ge -r_{\max} w_R(p_{\max}^T)(r - r_{\min})$. Additionally, since $w_R(p)$ is continuously differentiable, there exists $c$ such that $cp \ge w_R(p)$, which means that $\dot{p} \ge -cp(p_{\max}^T + 1)(r_{\max} - r_{\min})$. Then, at worst, the state $p$ (respectively, $r$) decays exponentially towards the value 0 (respectively, $r_{\min}$), which cannot be attained in finite time. ■

As a consequence of Proposition 5.3, the condition $H_1 = 0$ becomes

$$ (\text{Condition 1}) \qquad \lambda_r - \lambda_m = 0. $$

We define the quantity $\phi(\varphi, \lambda) \doteq (r_{\max} - m - r)\lambda_r - (p+1)\lambda_p - \lambda^0$, so that the time derivative of (Condition 1) is

$$ (\text{Condition 2}) \qquad \phi(\varphi, \lambda)w_R(p) - E_M\lambda_p = 0. $$

Along a singular arc, the Hamiltonian can be rewritten as

$$ (5.5) \qquad H = E_M m \lambda_p + \phi(\varphi, \lambda) w_R(p) (r - r_{\min}), $$

and, using (Condition 1) and (Condition 2), the adjoint system becomes

$$ \left\{ \begin{aligned} \frac{d\lambda_p}{dt} &= w_R(p)(r - r_{\min})\lambda_p - \phi(\varphi, \lambda)w'_R(p)(r - r_{\min}), \\ \frac{d\lambda_r}{dt} &= w_R(p)(r - r_{\min})\lambda_r - \phi(\varphi, \lambda)w_R(p). \end{aligned} \right. $$

**Proposition 5.4.** *Neither $\phi(\varphi, \lambda)$ nor $\lambda_p$ can vanish along a singular arc.*
Proof. According to (Condition 2), if either $\phi(\varphi, \lambda)$ or $\lambda_p$ are null, then both of them are null. Then, if $\phi(\varphi, \lambda) = \lambda_p = 0$, equation (5.5) would imply that the Hamiltonian vanishes in $\Sigma$, and therefore it would vanish for the whole interval $[0, T]$ (as it is constant along the solution). However, one can see in (5.1) that the Hamiltonian evaluated at final time is $-\lambda^0 w_R(p(T))(r(T) - r_{\min})$ which cannot be 0 due to Proposition 5.3 and $\lambda^0 \neq 0$. ■

We differentiate (Condition 2) w.r.t. time and we get $\dot{\phi}(\varphi, \lambda)w_R(p)+\phi(\varphi, \lambda)w'_R(p)\dot{p}-E_M\dot{\lambda}_p = 0$. Replacing the latter and using Proposition 5.4 allows us to reduce the expression to

$$ (Condition 3) \qquad -(r_{\max} - r_{\min})w_R(p)^2 + E_M(m + r - r_{\min})w'_R(p) = 0, $$

which allows us to express $m + r$ in terms of $p$.

Lemma 5.5. Along a singular arc over the interval $[t_1, t_2]$,

$$ m + r = x(p) $$

with $x(p) : \mathbb{R}_+ \to [r_{\min}, \infty)$ defined as

$$ x(p) \doteq (r_{\max} - r_{\min}) \frac{w_R(p)^2}{E_M w'_R(p)} + r_{\min}, $$

which, using (3.7), yields $x(p_{\text{opt}}^*) = r_{\max}$.

The fact that the control does not show up in (Condition 3)—which is obtained by differentiating (Condition 1) twice—means that the singular arc is at least of order two. We differentiate (Condition 3) and we get

$$ (Condition 4) \qquad (r_{\max} - x(p) + (p+1)x'(p))w_R(p)(r - r_{\min}) - E_M m x'(p) = 0. $$

We define the function

$$ (5.6) \qquad y(p) \doteq w_R(p)(r_{\max} - x(p) + (p+1)x'(p)). $$

Using (Condition 3) and (5.6) in (Condition 4) yields

$$ (x(p) - r_{\min})y(p) - (E_M x'(p) + y(p))m = 0, $$

which means we can express $m$ and $r$ in terms of $p$ along the singular arc.

Lemma 5.6. Along a singular arc over the interval $[t_1, t_2]$,

$$ (5.7) \qquad m = (x(p) - r_{\min}) \frac{y(p)}{E_M x'(p) + y(p)}, $$

$$ (5.8) \qquad r = x(p) - (x(p) - r_{\min}) \frac{y(p)}{E_M x'(p) + y(p)}. $$
We differentiate (Condition 4) and we get

$$ (5.9) \qquad \begin{aligned} & -(r_{\max}(1-\alpha) - m)w_R(p)(r - r_{\min}) + x'(p) \frac{y(p)}{E_M x'(p) + y(p)} \dot{p} \\ & +(x(p) - r_{\min}) \left( \frac{y'(p)}{E_M x'(p) + y(p)} - \frac{y(p)}{(E_M x'(p) + y(p))^2} (E_M x''(p) + y'(p)) \right) \dot{p} = 0, \end{aligned} $$

meaning that we can express

$$ \alpha_{\text{sing}}(p) = 1 - \frac{m}{r_{\max}} \left( \left( \frac{x'(p)}{x(p) - r_{\min}} + \frac{y'(p)}{y'(p)} - \frac{E_M x''(p) + y'(p)}{E_M x'(p) + y(p)} \right) \frac{\dot{p}}{w_R(p)(r - r_{\min})} + 1 \right). $$

While (Condition 3) showed that the order of the singular arc is at least two, the latter relation proves that it is exactly two. Indeed, the coefficient before $\alpha$ in (5.9) is $-r_{\max}w_R(p)(r - r_{\min})$, which cannot vanish as proven in Proposition 5.3. The singular arc is said to be locally of order two, as the coefficient of $\alpha$ in (Condition 3) is zero along the singular arc, but not everywhere on the cotangent bundle [34]. In this case, the presence of the Fuller phenomenon (i.e., the junctions between bang and singular arcs constituting an infinite number of switchings) is not guaranteed. However, this turns out to be the case as it will be shown in the numerical computations. Besides, in accordance with Proposition 5.2, the order two singular control can be expressed in feedback form, i.e., as a function of the state only. We performed a numerical rank test using Singular Value Decomposition, which confirmed that the rank condition is fulfilled. More precisely, the actual computation proves that the singular control can be expressed as a function of $p$ only (Lemma 5.6 entails that $r, m$ and therefore $\dot{p}$ can be expressed in terms of $p$), which allows to retrieve the turnpike behaviour as described in the following section.

**5.3.3. The turnpike phenomenon.** Using (5.7) and (5.8), we see that, along a singular arc, the dynamical equation of $p$ becomes

$$ \dot{p} = E_M w_R(p) \frac{x(p) - r_{\min}}{E_M x'(p) + y(p)} (r_{\max} - x(p)), $$

which is only equal to 0 when $r_{\max} = x(p)$. This is only true at $p = p_{\text{opt}}^*$, and so

$$ \text{sign}(\dot{p}) = \text{sign}(p_{\text{opt}}^* - p), $$

meaning that, in a singular arc over the interval $[t_1, t_2]$, the concentration $p$ converges asymptotically to the optimal value $p_{\text{opt}}^*$. This means that $m$ and $r$ would also converge to the optimal values $m_{\text{opt}}^*$ and $r_{\text{opt}}^*$, respectively, and the singular control $\alpha_{\text{sing}}$ to $\alpha_{\text{opt}}^*$. We formalize this in the following theorem.

**Theorem 5.7.** On a singular arc, the system states and singular control tend asymptotically to

$$ (p, r, m) = (p_{\text{opt}}^{*}, r_{\text{opt}}^{*}, m_{\text{opt}}^{*}), $$

$$ \alpha_{\text{sing}}(t) = \alpha_{\text{opt}}^{*}. $$
The above theorem is an explicit proof of the presence of the turnpike property: an optimal control characterized by a singular arc that stays exponentially close to the steady-state solution of the static optimal control problem [35]. This phenomenon has been considerably studied in econometry [36], and more recently in biology [37, 10, 20]. It has been shown that, for large final times, the trajectory of the system spends most of the time near the optimal steady state, and that in infinite horizon problems, it converges to this state.

**5.4. Numerical results.** The computations of the optimal trajectories were performed with Bocop [38], which solves the optimal control problem through a direct method. An online version of the numerical computations can be visualized and executed on the gallery of the ct (Control Toolbox) project<fnref n="3" />. The time discretization algorithm used is Lobato IIIC (implicit, 4-stage, order 6) with 2000 time steps. Figure 3 shows an optimal trajectory with $r(0) + m(0) < r_{max}$, where most of the bacterial mass corresponds to class Q proteins. The obtained optimal control confirms the conclusions of the latter section: a large part of the time, the optimal control remains near the optimal steady-state allocation $\alpha_{opt}^*$, according to the turnpike theory (Theorem 5.7). The solution presents chattering after and before the singular arc, as expected in the presence of Fuller's phenomenon (even if only a finite number of bangs is computed by the numerical method), and the final bang corresponds to $\alpha = 1$ (Lemma 5.1). In order to verify the optimality of the singular arc, we performed a numerical computation of the derivatives of $H_1$, which is shown in Figure 4. The fact that the factor of $\alpha$ in the fourth derivative is different from 0 confirms that the singular arc is of order 2. Moreover, its negativity complies with the generalized Legendre-Clebsch condition given by

$$ (5.10) \qquad (-1)^k \frac{\partial}{\partial \alpha} \left( \frac{d^{2k}}{dt^{2k}} H_1 \right) < 0, $$

along the singular arc, which is a necessary condition for optimality. As we state in [23], even if there exist no available sufficient condition to verify local optimality of extremals with Fuller arcs, a check of the Legendre-Clebsch condition along the singular arc can ensure that the extremal obtained is not a too crude local minimizer. For the second-order singular arc case, the condition corresponds to the case $k = 2$. The initial conditions used in Figure 3 were only chosen to confirm the theoretical results found throughout this section, by emphasizing the main features of the solution. However, from a biological perspective, a situation where $r + m$ is significantly different from its steady-state value $r_{max}$ is not to be expected: a common assumption in these classes of coarse-grained models is that the transcription of Q proteins is autoregulated around stable levels [39], which translates into a constant $q = 1 - r_{max}$ (and therefore $m + r = r_{max}$) for the whole interval $[0, T]$. We will see in next section that this hypothesis produces a very particular structure of the optimal control solution.

**6. Biologically relevant scenarios.** Despite their simplicity, self-replicator models have been capable of accounting for a number of observable phenomena during steady-state microbial growth, under the assumption that bacteria allocate their resources in such a way as to maximize growth. Here, we apply the general optimal allocation strategy derived in the previous section to predict the bacterial response to certain environmental changes. We

<footnote n="3">³https://ct.gitlabpages.inria.fr/gallery/bacteria/bacteria.html</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3: Numerical simulation of (OCP) obtained with Bocop, for the parameter values derived in Section 4. Initial state is $p(0) = 0.03$, $r(0) = 0.1$, $m(0) = 0.2$ with $E_M = 0.6$. As predicted, the optimal control $\alpha$ involves chattering after and before the singular arc. The mass fraction $q$ converges to $1 - r_{max}$ and $m + r$ to $r_{max}$. Moreover, along the singular arc, the states $(p^*, r^*, m^*)$ converge asymptotically to $(p_{opt}^*, r_{opt}^*, m_{opt}^*)$.</figcaption></figure>

consider two situations that commonly affect bacteria: changes in the nutrient concentration in the medium, and changes in the environment submitting the cell to a particular stress.

**6.1. Nutrient shift.** Bacteria are known to traverse different habitats throughout their lifetime, experiencing fluctuating nutrient concentrations in the medium. In [10], we explored how bacteria dynamically adjust their allocation strategy when facing a nutrient upshift. In this work, we show that considering a class of growth rate-independent proteins in the model refines these previous results. We consider the optimal control problem with the initial state being the optimal steady state for a low value of $E_M$, and we set a higher $E_M$ for the time interval $[0, T]$, representing a richer medium. Setting initial conditions at steady state has an impact on the singular arc of the optimal control: it holds that $m + r = r_{max}$ and $q = 1 - r_{max}$ for the whole trajectory, which yields a constant singular arc.

**Theorem 6.1.** If $r(0) + m(0) = r_{max}$ (i.e., $q$ starts from a steady-state value), then any singular arc over the interval $[t_1, t_2]$ of the optimal control corresponds to the optimal steady state.
<figure><img src="image_5.png" /><figcaption>Figure 4: Factors of $\alpha$ in the derivatives of $H_1$ evaluated over the trajectory plotted in Figure 3. The intervals where the functions vanish are marked in red. As expected, all functions vanish along the singular arc except for the factor in the fourth derivative (highlighted in green) which is negative according to the Legendre-Clebsch condition (5.10).</figcaption></figure>

Proof. The dynamical equation for $q$ is $\dot{q} = ((1 - r_{\max}) - q)w_R(p)(r - r_{\min})$, where it can be seen that the set $q = 1 - r_{\max}$ is invariant. This means that, for any trajectory emanating from a steady state, $q$ remains constant even under changes of the nutrient quality $E_M$. Then, by using the relation (2.5), we obtain

$$ (6.1) \qquad m + r = r_{\max}. $$

Along the singular arc, it holds that $m + r = x(p)$, which, using (6.1), implies that $p = p_{\text{opt}}^*$, meaning that the precursor concentration along the singular arc is constant and optimal. Then, $\alpha_{\text{sing}} = \alpha_{\text{opt}}^*$, $m = m_{\text{opt}}^*$ and $r = r_{\text{opt}}^*$ for the whole singular arc. ■

A numerical simulation of this scenario is shown in Figure 5. As expected, the increase in $E_M$ produces a higher ribosomal mass fraction $r$, which translates into an increase of the growth rate, stabilizing at the maximal steady-state growth rate $\mu_{\text{opt}}^*$ through an oscillatory phase. It is noteworthy that, in comparison to Giordano et al.'s model, the relative changes in mass fractions $r$ and $m$ are much lower, which corresponds well with the relative changes observed in [6]. Additionally, while the presence of $r_{\min}$ does not noticeably affect the solution
<figure><img src="image_6.png" /></figure>

Figure 5: Numerical simulation of the optimal control problem starting from a steady state. The initial state corresponds to the optimal steady state for $E_M = 0.3$ (poor medium), and the new environmental constant is fixed to $E_M = 0.7$ (rich medium). As predicted, $m + r (= 1 - q)$ remains constant, even if they vary individually, in opposition to the previous case. Naturally, an increase in the nutrient quality produces a higher steady-state ribosomal mass fraction $r^*$, which yields an increased steady-state growth rate $\mu_{opt}^*$ with respect to the growth rate before the upshift.

of the optimal control problem, it contributes to a model that more accurately reproduces the experimental data (Figure 2a), representing a significant improvement from the modeling point of view.

**6.2. Bacterial response to stress.** The other scenario of interest is an environmental change imposing a certain stress on the microbial population, which is counteracted through the synthesis of a stress response protein W. This protein is also growth rate-independent like Q, and its production can be triggered by many different situations. For instance, when subject to extreme temperatures, the production of so-called molecular chaperones helps bacteria counter the effect of protein unfolding [40, 14]. Likewise, the production of other proteins is known to protect bacteria like *E. coli* against acid stress [41]. Another possible scenario is the response to metabolic load imposed by the induced overexpression of a heterologous protein [42]. All of these situations are known to reduce the resources available for growth-associated proteins (Figure 6), consequently decreasing the maximal growth rate attainable. Here, we
model a general stress response through the production of the W protein that takes up a fraction w of the proteome, thus reducing $r_{max}$ to a certain $r_{max}^w < r_{max}$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Left: original case. Right: new proposed case, where q remains unchanged, but the maximal allocation $m + r$ is restricted to a $r_{max}^w < r_{max}$.</figcaption></figure>

As before, we assume q takes up a constant fraction $1 - r_{max}$ of the proteome, but the proportions of resources allocated to M and R are now $r_{max}^w\alpha$ and $r_{max}^w(1 - \alpha)$ respectively. By construction, we have $w = r_{max} - m - r$, which means we can express

$$ \dot{w} = (r_{max} - r_{max}^w - w)w_R(p)(r - r_{min}), $$

showing that the mass fraction w converges asymptotically to the difference $r_{max} - r_{max}^w$. The remaining mass fractions p, r and m obey the dynamics of system (S'), so the application of the optimal solution found in last section is straightforward. An example is shown in Figure 7. As predicted, $m + r$ converges to the reduced $r_{max}^w$, q remains constant at $1 - r_{max}$ and w converges to $r_{max}^w - r_{max}$. The reduction of resources available for growth-associated proteins (M and R) causes the growth rate to drop, as was shown experimentally [6].

**7. Conclusion.** In this work, we proposed a dynamical self-replicator model of bacterial growth based on the work of [10], which introduces a growth rate-independent class of protein. As a consequence, the proteome of the bacterial cell can be divided into the metabolic machinery M, the gene expression machinery R, and the housekeeping machinery Q. While Q is growth rate-independent, this is also the case for a fraction of R required for cell replication to occur. As a consequence of this hypothesis, a maximum ribosomal concentration $r_{max}$ appears in the model kinetics, limiting the allocation of resources to M and R. We studied the asymptotic behavior of the system, showing that, under certain conditions, all solutions converge towards the only globally attractive equilibrium. We then explored the optimal dynamic allocation strategies that consider maximizing the bacterial population volume in terms of the resource allocation parameter $\alpha$. This involved a study of the static and dynamic aspects of optimal strategies. For the first one, we showed there is a unique optimal steady state, which corresponds to experimental observations of growing cultures of E. coli [29, 30, 6, 27]. The dynamic problem is approached through optimal control theory, by application of the Pontrjagin's Maximum Principle. The obtained optimal control has a Fuller-singular-Fuller structure with a non-constant singular arc, in contrast to the constant singular arc obtained in Giordano et al.'s approach. We performed a detailed analysis of the OCP in both analytic and numerical ways. In particular, the singular arc of the optimal solution is characterized
<figure><img src="image_8.png" /><figcaption>Figure 7: Numerical simulation of an optimal trajectory where the initial conditions are the optimal steady state for $E_M = 0.7$ and $r_{max} = 0.5$. A certain stress is induced at $t = 0$, which triggers the synthesis of the growth rate-independent protein *w*, reducing the fraction $r_{max}$ to $r_{max}^w = 0.3$. As a result, the steady-state growth rate is significantly reduced.</figcaption></figure>

by i) its feedback form (i.e., being expressed as a function of the state only), ii) being exactly of order 2, and iii) the turnpike phenomenon (where the state trajectory and optimal control converge asymptotically towards the optimal steady state and control). Moreover, we showed that, when the mass fraction of class *Q* proteins is at steady state, the singular arc of the optimal solution corresponds to the optimal steady state. Additionally, we showed that the dynamical approach can be used to predict the behavior of the system when subject to stress. The latter is modeled through a reduction of the fraction of growth rate-dependent protein synthesis as the production of a *w* protein that reduces $r_{max}$.

While the main features of Giordano *et al.'s* work are present in this approach, our generalization shows a better agreement with the experimental data given by the introduction of the parameters $r_{max}$ and $r_{min}$ in the model. Additionally, the proposed partitioning of the proteome in a dynamic setting can account for certain natural phenomena known to reduce the fraction of growth rate-dependent proteins in the cell. These modifications yield interesting optimal control problems, which could potentially help understand the internal decision-making mechanisms evolved by bacteria.

Our approach was built on the joint exploitation of theoretical and numerical results.
806 When tackling more complex problems as proposed, e.g., in Tsiantis & Banga [43], a PMP
807 perspective tends to yield very complicated mathematical formulations. Using direct methods
808 has the advantage of avoiding these issues, but it often requires some knowledge to initialize
809 the optimization algorithm or to check the validity of the solutions. In order to investigate
810 complex biological systems, we advocate the development and theoretical analysis of simple
811 models, in line with the question to be investigated, coupled with numerical exploration of
812 optimal solutions (using larger models if necessary).

813 **Acknowledgments.** We would like to acknowledge the help of S. Psalmon and B. Schall
814 from Polytech Nice Sophia for the numerical simulations.

REFERENCES

816 [1] M. SCHAECHTER, J. L. INGRAHAM, AND F. C. NEIDHARDT, *Microbe*, ASM Press, Washington, DC, 2006.

817 [2] F. NEIDHARDT, *Bacterial growth: constant obsession with dN/dt*, Journal of Bacteriology, 181 (1999), pp. 7405-7408.

820 [3] S. BRUL, S. VAN GERWEN, AND M. ZWIETERING, *Modelling Microorganisms in Food*, Elsevier, 2007.

821 [4] A. BRAUNER, O. FRIDMAN, O. GEFEN, AND N. Q. BALABAN, *Distinguishing between resistance, tolerance and persistence to antibiotic treatment*, Nature Reviews Microbiology, 14 (2016), pp. 320-330.

823 [5] J. C. LIAO, L. MI, S. PONTRELLI, AND S. LUO, *Fuelling the future: microbial engineering for the production of sustainable biofuels*, Nature Reviews Microbiology, 14 (2016), pp. 288-304.

825 [6] M. SCOTT, C. W. GUNDERSON, E. M. MATEESCU, Z. ZHANG, AND T. HWA, *Interdependence of cell growth and gene expression: origins and consequences*, Science, 330 (2010), pp. 1099-1102.

827 [7] H. VAN DEN BERG, Y. KISELEV, AND M. ORLOV, *Optimal allocation of building blocks between nutrient uptake systems in a microbe*, Journal of Mathematical Biology, 44 (2002), pp. 276-296.

829 [8] D. MOLENAAR, R. VAN BERLO, D. DE RIDDER, AND B. TEUSINK, *Shifts in growth strategies reflect tradeoffs in cellular economics*, Molecular Systems Biology, 5 (2009), p. 323.

831 [9] A. Y. WEISSE, D. A. OYARZÚN, V. DANOS, AND P. S. SWAIN, *Mechanistic links between cellular trade-offs, gene expression, and growth*, Proceedings of the National Academy of Sciences, 112 (2015), pp. E1038-E1047.

834 [10] N. GIORDANO, F. MAIRET, J.-L. GOUZÉ, J. GEISELMANN, AND H. DE JONG, *Dynamical allocation of cellular resources as an optimal control problem: novel insights into microbial growth strategies*, PLoS Computational Biology, 12 (2016), p. e1004802.

837 [11] D. W. ERICKSON, S. J. SCHINK, V. PATSALO, J. R. WILLIAMSON, U. GERLAND, AND T. HWA, *A global resource allocation strategy governs growth transition kinetics of Escherichia coli*, Nature, 551 (2017), pp. 119-123.

840 [12] H. DOURADO AND M. J. LERCHER, *An analytical theory of balanced cellular growth*, Nature Communications, 11 (2020), pp. 1-14.

842 [13] M. SCOTT, S. KLUMPP, E. MATEESCU, AND T. HWA, *Emergence of robust growth laws from optimal regulation of ribosome synthesis*, Molecular Systems Biology, 10 (2014), p. 747.

844 [14] F. MAIRET, J.-L. GOUZÉ, AND H. DE JONG, *Optimal proteome allocation and the temperature dependence of microbial growth laws*, npj Systems Biology and Applications, 7 (2021), pp. 1-11.

846 [15] M. BASAN, S. HUI, H. OKANO, Z. ZHANG, Y. SHEN, J. WILLIAMSON, AND T. HWA, *Overflow metabolism in Escherichia coli results from efficient proteome allocation*, Nature, 528 (2015), pp. 99-104.

849 [16] A. MAITRA AND K. DILL, *Bacterial growth laws reflect the evolutionary importance of energy efficiency*, Proceedings of the National Academy of Sciences of the USA, 112 (2015), pp. 406-411.

851 [17] D. A. CARLSON, A. B. HAURIE, AND A. LEIZAROWITZ, *Infinite Horizon Optimal Control*, Springer, Berlin, Heidelberg, 1991.

853 [18] I. YEGOROV, F. MAIRET, AND J.-L. GOUZÉ, *Optimal feedback strategies for bacterial growth with degra-
dation, recycling, and effect of temperature, Optimal Control Applications and Methods, 39 (2018),
pp. 1084–1109.

[19] J. IZARD, C. G. BALDERAS, D. ROPERS, S. LACOUR, X. SONG, Y. YANG, A. LINDNER, J. GEISELMANN, AND H. DE JONG, A synthetic growth switch based on controlled expression of RNA polymerase, Molecular Systems Biology, 11 (2015), p. 840.

[20] I. YEGOROV, F. MAIRET, H. DE JONG, AND J.-L. GOUZÉ, Optimal control of bacterial growth for the maximization of metabolite production, Journal of Mathematical Biology, 78 (2019), pp. 985–1032.

[21] A. G. YABO, J.-B. CAILLAU, AND J.-L. GOUZÉ, Singular regimes for the maximization of metabolite production, in 2019 IEEE 58th Conference on Decision and Control (CDC), IEEE, 2019, pp. 31–36.

[22] A. G. YABO AND J.-L. GOUZÉ, Optimizing bacterial resource allocation: metabolite production in continuous bioreactors, IFAC-PapersOnLine, 53 (2020), pp. 16753–16758.

[23] A. G. YABO, J.-B. CAILLAU, AND J.-L. GOUZÉ, Optimal bacterial resource allocation: metabolite production in continuous bioreactors, Mathematical Biosciences and Engineering, 17 (2020), pp. 7074–7100.

[24] J. EWALD, M. BARTL, AND C. KALETA, Deciphering the regulation of metabolism with dynamic optimization: an overview of recent advances, Biochemical Society Transactions, 45 (2017), pp. 1035–1043.

[25] E. KLIPP, R. HEINRICH, AND H.-G. HOLZHÜTTER, Prediction of temporal gene expression: Metabolic optimization by re-distribution of enzyme activities, European Journal of Biochemistry, 269 (2002), pp. 5406–5413.

[26] D. A. OYARZÚN, B. P. INGALLS, R. H. MIDDLETON, AND D. KALAMATIANOS, Sequential activation of metabolic pathways: a dynamic optimization approach, Bulletin of Mathematical Biology, 71 (2009), p. 1851.

[27] S. HUI, J. M. SILVERMAN, S. S. CHEN, D. W. ERICKSON, M. BASAN, J. WANG, T. HWA, AND J. R. WILLIAMSON, Quantitative proteomic analysis reveals a simple strategy of global resource allocation in bacteria, Molecular Systems Biology, 11 (2015), p. 784.

[28] M. BASAN, M. ZHU, X. DAI, M. WARREN, D. SÉVIN, Y.-P. WANG, AND T. HWA, Inflating bacterial cells by increased protein synthesis, Molecular Systems Biology, 11 (2015), p. 836.

[29] J. FORCHHAMMER AND L. LINDAHL, Growth rate of polypeptide chains as a function of the cell growth rate in a mutant of Escherichia coli 15, Journal of Molecular Biology, 55 (1971), pp. 563–568.

[30] H. BREMER, P. P. DENNIS, ET AL., Modulation of chemical composition and other parameters of the cell by growth rate, Escherichia coli and Salmonella: Cellular and Molecular Biology, 2 (1996), pp. 1553–69.

[31] J. S. EDWARDS, R. U. IBARRA, AND B. O. PALSSON, In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data, Nature Biotechnology, 19 (2001), pp. 125–130.

[32] N. E. LEWIS, K. K. HIXSON, T. M. CONRAD, J. A. LERMAN, P. CHARUSANTI, A. D. POLPITIYA, J. N. ADKINS, G. SCHRAMM, S. O. PURVINE, D. LOPEZ-FERRER, ET AL., Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models, Molecular Systems Biology, 6 (2010), p. 390.

[33] A. A. AGRACHEV AND Y. SACHKOV, Control Theory from the Geometric Viewpoint, vol. 87, Springer Science & Business Media, 2013.

[34] M. I. ZELIKIN AND V. F. BORISOV, Theory of Chattering Control: with Applications to Astronautics, Robotics, Economics, and Engineering, Springer Science & Business Media, 2012.

[35] E. TRÉLAT AND E. ZUAZUA, The turnpike property in finite-dimensional nonlinear optimal control, Journal of Differential Equations, 258 (2015), pp. 81–114.

[36] D. CASS, Optimum growth in an aggregative model of capital accumulation, The Review of Economic Studies, 32 (1965), pp. 233–240.

[37] J.-M. CORON, P. GABRIEL, AND P. SHANG, Optimization of an amplification protocol for misfolded proteins by using relaxed control, Journal of Mathematical Biology, 70 (2015), pp. 289–327.

[38] I. S. TEAM COMMANDS, BOCOP: an open source toolbox for optimal control. http://bocop.org, 2017.

[39] S. KLUMPP, Z. ZHANG, AND T. HWA, Growth rate-dependent global effects on gene expression in bacteria, Cell, 139 (2009), pp. 1366–1375.

[40] D. RONCARATI AND V. SCARLATO, Regulation of heat-shock genes in bacteria: from signal sensing to gene expression output, FEMS Microbiology reviews, 41 (2017), pp. 549–574.
908

909

910

911

912

913

914

[41] N. GUAN AND L. LIU, *Microbial response to acid stress: mechanisms and applications*, Applied Microbiology and Biotechnology, 104 (2020), pp. 51-65.

[42] H. DONG, L. NILSSON, AND C. G. KURLAND, *Gratuitous overexpression of genes in Escherichia coli leads to growth inhibition and ribosome destruction.*, Journal of Bacteriology, 177 (1995), pp. 1497-1504.

[43] N. TSIANTIS AND J. R. BANGA, *Using optimal control to understand complex metabolic pathways*, BMC Bioinformatics, 21 (2020), pp. 1-33.