Then, the optimal control problem is defined as

$$ 
\left\{
\begin{array}{ll}
\underset{u}{\text{maximise}} & \text{metabolite production } J_X(u) \\
\text{subject to} & \text{dynamics (S_N)} \\
& \text{initial conditions (IC)} \\
& \quad u(\cdot) \in \mathcal{U}, \\
& t \in [0, T].
\end{array}
\right.
\qquad (\text{OCP}_X)
$$

### 5.3.2 On the solution of the OCP

Applying PMP, we see that the Hamiltonian is affine in the control, so it has the form $H = H_0 + uH_1$, meaning that the solution is bang-singular-bang, given by

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

Examples of optimal trajectories are shown in Figure 5.4 and Figure 5.5, where both problems OCP_N and OCP_X are compared for different environmental conditions representing rich and poor qualities of the nutrient in the medium. While the structures of the optimal control for both problems are similar, the optimal strategy maximizing the production of X is characterized by a non-constant singular arc, which is close to the solution $u_{opt}$ of the static OCP, but deviates from it towards the end. Additionally, the times at which the junctions between bang and singular arc are produced differ, as well as the values of the bangs. In particular, in both Figures, the final bang of the natural control is $\alpha = 1$, while that of the artificial control is $u = 0$. More detailed calculations of the PMP approach can be found in [3].

### 5.3.3 Product maximisation including naturally-evolved al-location

The comparison between OCP_N and OCP_X proves useful to observe that the natural behavior of microbes does not necessarily match the artificial objective of producing a certain metabolite. However, the assumption made in OCP_X is a pure theoretical one, as $\alpha$ cannot be completely substituted by the external control. We then propose
<figure><img src="image_30.png" /><figcaption>Figure 5.4: Optimal control obtained with Bocop. Simulation in a rich medium with $e_M = k_R$, meaning that the substrate enables the maximum growth rate. Initial conditions are $p_0 = 0.024$, $r_0 = 0.1$ and the simulation time is set to $T = 30$.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 5.5: Optimal control obtained with Bocop. Simulation in a poor medium with $e_M = 0.5k_R$. Initial conditions are $p_0 = 0.024$, $r_0 = 0.3$ and the simulation time is set to $T = 30$.</figcaption></figure>

an approach that takes into account two distinct processes: (i) the ability of bacteria to maximise their biomass through the optimal allocation described by the internal control problem $OCP_N$, (ii) the external action of an operator intending to maximise the production of the metabolite of interest X. In practice, the new pathway associated to the dynamics of X in (5.2) is obtained by optogenetic engineering of a strain of bacteria: a light-induced control I is able to externally modulate the natural allocation process. This is modeled by replacing the control $u$ of $OCP_X$ by $u = \alpha(p, r)I$, so that the external control affects in a multiplicative way the internal allocation strategy. The crucial difference with the previous formulation [3] is that the internal
control $\alpha$ of the bacteria now appears in feedback form (thus depending on the two states of OCP$_N$; see Remark 5.2.1), and is modulated by the external light-induced control I. It is noteworthy that no competition occur between the two objectives (biomass vs. metabolite production maximisation). The approach considers the internal control in feedback form, which is mitigated by an external control in relation with a global process that includes the new pathway to produce the metabolite. In this context, the proposed hierarchical approach proves to be more relevant from a biological point of view than a multicriterion one. Additionally, we assume that the feedback $\alpha(p, r)$ is known and smooth. While the latter seems to be a strong assumption from the control point of view (as the solution of OCP$_N$ comprises bang and singular arcs), it is a reasonable assumption in our biological setting, where the kinetics of the involved biochemical reactions prescribe continuous behaviours (see [30] for biologically relevant approximations of the feedback). Thus, the dynamical equation of $r$ becomes

$$ \dot{r} = (r_{\max}\alpha(p, r)I - r)w_R(p)(r - r_{\min}), \quad (5.5) $$

where the new control $I(t)$ is subject to bound constraints $0 \le I(t) \le I_{\max}$, and the cost $J_X$ remains unchanged. This defines problem $\overline{OCP}_X$.

### 5.3.4 Hierarchical MPC for metabolite production

In this approach, we approximate the allocation feedback $\alpha(p, r)$ through the sh-MPC loop described in 5.2.3. At each iteration $k$, solving on $[k\tau, T]$ (where $T$ is, as before, the fixed horizon) yields an approximation of $\alpha(p_k, r_k)$ (based on the suboptimal parametric form $\alpha_{so}$), and of $\alpha$ evaluated at further steps. Then, this suboptimal feedback is injected in the dynamics (5.5) of $\overline{OCP}_X$ so as to find the optimal external control $I$ maximizing $X(T)$. Thus, a second MPC is used "above" the first one. There is a quite large literature on such approaches combining several MPC loops (see, e.g., [80, 81] and references therein). Other relevant matters such as using different time grids for each MPC loop or, more generally, synchronisation issues, are not discussed here (see also the recent paper [82] on convergence of MPC methods in finite horizon). Instead, we focus on the biological application. We note, in particular, that when the feedback $\alpha(p, r)$ is zero (which would be expected for the genuine—though biologically unrealistic—feedback of OCP$_N$ as zero bang arcs
can occur), the external control *I* is not active. In practice, when the allocation is close to zero, the MPC loop would compensate through *I* for this discrepancy. As the external control *I* is bounded, the latter can induce certain performance loss between the results of the ideal model OCP<sub>X</sub> and the more realistic problem $\overline{OCP}_X$. Such comparisons are provided in the next paragraph.

## 5.4 Numerical results

Figure 5.6 shows a comparison between the optimal trajectory solution of OCP<sub>X</sub> and the hierarchical MPC proposed in this paper. Differences between both trajectories are marginal, mainly given by the approximation of the singular arc by a constant control *u*. In Figure 5.7 we see how, in order to match the optimal control *u* solution of OCP<sub>X</sub>, the external signal *I* completely arrests the natural allocation $\alpha$ (around *t* = 12), which implies allocating all the cellular resources to the metabolic machinery M, thus catalyzing the synthesis of X. The latter produces the suboptimal control $\alpha_{so}$ to increasingly compensate until it reaches the value 1. This result shows the difference between the simulated closed-loop behavior of the natural allocation strategy $\alpha_{so}(p, r)$, and the open-loop one (which remains constant almost over the whole interval $[0, T]$).

## 5.5 Conclusion

In this paper, we presented a hierarchical sh-NMPC approach to the problem of optimally producing a metabolite of interest in bacteria. The scheme is based on a parametric version of the naturally-evolved research allocation strategy proposed in [30], which represents a closed-loop alternative to these existing open-loop studies. A second sh-NMPC loop is applied in a hierarchical manner, in order to achieve the metabolite maximisation objective while taking into account the closed-loop natural control. Despite the approach being at an early stage, with no experimental results, it represents a step towards plausible biosynthetic real-time implementations. In future works, we are interested in comparing the natural MPC approximation with alternatives proposed in the literature, such as ppGpp regulation [30]. Other extensions include the possibility of estimating the real value of the natural allocation through online identification techniques.
<figure><img src="image_32.png" /><figcaption>Figure 5.6: Comparison of the optimal control $u(t)$ solution of OCP$_X$ and the hierarchical MPC scheme parametrized that considers the natural allocation as an inner MPC loop. Initial conditions are set to $p_0 = 0.024$, $r_0 = 0.3$, and $V_0 = 0.003$. Final time is set to $T = 30$, the scheme is executed with time step $\tau = 1$ and the environmental constant $e_M = 0.5k_R$.</figcaption></figure>

## Acknowledgments

The authors thank Hidde de Jong (Inria) as well as colleagues from the ANR Maximic project for many inputs on the biological problem.
<figure><img src="image_33.png" /><figcaption>Figure 5.7: Final control *u* and external signal *I* obtained from the MPC loop simulated in Figure 5.6.</figcaption></figure>
# Chapter 6

## A generalized resource allocation model of microbial growth

This chapter reproduces [AYT1], accepted for publication in the SIAM Journal on Applied Dynamical Systems.

### 6.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 [83]. 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 [84], 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 production [85, 86, 87].

A fruitful perspective on microbial growth is to view it as a resource allocation problem [2]. 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 [2, 88, 26, 16, 30, 33, 89]. 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 [26, 2, 23, 30, 90]. 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 [26, 91, 92].

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 [30]. The prediction was based on the Infinite Horizon Maximum Principle, a generalization of the well-known PMP (Pontrjagin Maximum Principle) [93, 43]. 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 [94], 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 [3, AYT5, AYT4, AYT2]. 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 op-
timization problems in biology (see [95] for a review). A classical example is the
determination of optimal activation patterns of metabolic pathways, such as to min-
imize the transition time of metabolites or minimize enzyme costs [96, 97].

The resource allocation model that lies at the basis of the above-mentioned work
[30] 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 [98]. This suggests the introduction
of a third protein category, dedicated mainly to basic 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 [2, 23].
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 [2, 98]. 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 [23].

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 mathemat-
ical 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 charac-
terized 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.

## 6.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 6.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. [2], 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 $\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 precur-
<figure><img src="image_34.png" /><figcaption>Figure 6.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 $γα$, $γ(1 - α)$, and $1 - γ$, respectively. Solid lines indicate the macroreactions with their respective synthesis rates, and dashed lines denote a catalytic effect.</figcaption></figure>

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 6.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 $α(t) \in [0, 1]$. Thus, the proportion of the total synthesis rate of proteins dedicated to the gene expression machinery R is $γα$, while that of the metabolic machinery M is $γ(1 - α)$. In particular, the allocation parameter does not influence the synthesis rate of Q, with constant proportion $1 - γ$, as the synthesis of proteins in this class is auto-regulated through mechanisms not relevant in this study. From a biological perspective, the function $α(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, $α$ is treated as the control input of the system.
### 6.2.1 Self-replicator system

Generalizing upon Giordano et al. [30], a mass balance analysis yields the dynamical system

$$ egin{cases} \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{cases} $$

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 [83], and so we define the bacterial volume $V$ measured in liter units (L) as

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

where $\beta$ corresponds to a density constant relating mass and bacterial volume [99], 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

$$ 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}. \qquad (6.2) $$

Using (6.1) and (6.2), we to obtain the relation

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

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{\mathcal{V}}}{\mathcal{V}} = \frac{\beta V_R}{\mathcal{V}} = \beta v_R(p_{\mathcal{V}}, r_{\mathcal{V}}). $$

We write the system in terms of the concentrations as

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

### 6.2.2 Kinetic definition

We define the kinetics of the reaction system by taking into account that a minimal concentration of ribosomal proteins $r_{\mathcal{V,min}} \in (0, \gamma/\beta)$ 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 [23]. Such behavior can be modeled as

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

Later on, we will see that there is no need to define $v_R(p_{\mathcal{V}}, r_{\mathcal{V}})$ for $r_{\mathcal{V}} < r_{\mathcal{V,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_{\mathcal{V}}) \doteq w_M(s) m_{\mathcal{V}}. $$

We will make the following assumption for functions $w_R(p_{\mathcal{V}})$ and $w_M(s)$.

**Hypothesis 6.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$,

* Strictly increasing: $w'_i(x) > 0, \forall x \ge 0$,

* Strictly concave: $w''_i(x) < 0, \forall x \ge 0$,

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

The classical Michaelis-Menten kinetics satisfies Hypothesis 6.2.1. While most of the mathematical 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

$$w_R(p_v) \doteq k_R \frac{p_v}{K_R + p_v}, \quad w_M(s) \doteq k_M \frac{s}{K_S + s}, \qquad (6.4)$$

where $k_R$ and $k_M$ are the maximal reaction rates in h⁻¹, and $K_M$ and $K_R$ are the half-saturation constants of the synthesis rates in g L⁻¹. For the general case introduced in Hypothesis 6.2.1 we will define

$$k_R \doteq \lim_{p_v \to \infty} w_R(p_v).$$

### 6.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 6.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.
### 6.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 (6.3), yields the relation

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

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

$$\hat{\mu}(p, r) \doteq \frac{\mu(p_v, r_v)}{k_R} = \hat{w}_R(p)(r - r_{\min}), \qquad (6.6)$$

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{aligned}
\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{aligned}
\right.
\qquad (S)$$

where $q$ has been removed since it can be expressed in terms of the other concentrations through (6.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_{\min}$). In what follows, we will systematically investigate how these differences affect the asymptotic behavior and optimal resource allocation strategies.
## 6.3 Asymptotic behavior

In the present section, we will study the asymptotic behavior of the reduced system representing the intracellular dynamics

$$ egin{cases} 
\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{cases} 
\quad (6.7) $$

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 6.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 (6.7).

*Proof.* This can be easily verified by evaluating the differential equations of system (6.7) 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.

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

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 6.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 6.3.2, we redefine (6.7) as

$$ egin{cases} \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{cases} \quad (S') $$

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.

## 6.3.1 Local stability

**Theorem 6.3.1.** 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

$$ egin{aligned} p^* &\doteq \left\{ p \in \mathbb{R}_+ : (p+1)w_R(p) = \frac{E_M m^*}{r^* - r_{\min}} \right\}, \quad (6.8) \\ r^* &\doteq r_{\max}\alpha^*, \\ m^* &\doteq r_{\max}(1 - \alpha^*). \end{aligned} $$
*Proof.* The general Jacobian matrix of the system (S') is

$$ \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}. \quad (6.9) $$

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 6.2.1), and $E_M m^*/(r^* - r_{\min}) > 0$, so the set (6.8) yields a unique solution. For $\alpha^* < \alpha_{\min}^*$, the equation for $p^*$ in (6.8) has no valid solution as $E_M m^*/(r^* - r_{\min})$ becomes negative, and therefore the equilibrium does not exist. The Jacobian (6.9) 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)). $$

If $\alpha^* > \alpha_{\min}^*$, then $J_2$ has one positive eigenvalue and $E_2$ becomes locally unstable.
As for $E_3$, 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}$$

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

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

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

### 6.3.2 Global behavior

We will study the global behavior of system (S') for the initial conditions

$$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. \quad (\text{IC})$$

and for a given constant allocation parameter

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

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}),$$

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$ and $m^* - m$, respectively (and both $\dot{r}$ and $\dot{m}$ are zero if $p = 0$ or $r = r_{\min}$).
Then, let us divide $\Gamma$ into the subsets

$$ \begin{aligned} \mathcal{R}^{-} &\doteq \{(p, r, m) \in \Gamma : r \in (r_{\min}, r^{*})\}, \quad \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})\}, \quad \mathcal{M}^{+} \doteq \{(p, r, m) \in \Gamma : m \in (m^{*}, r_{\max})\}, \end{aligned} $$

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

**Lemma 6.3.3.** For $\alpha(t) = \alpha^{*} \in (\alpha_{\min}^{*}, 1)$, the closed sets $\overline{\mathcal{R}}^{-}$, $\overline{\mathcal{R}}^{+}$, $\overline{\mathcal{M}}^{-}$ and $\overline{\mathcal{M}}^{+}$ are invariant 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. $$

Again, the invariance of the sets can be checked by evaluating the vector field over the boundaries of the sets. Now we state a first result.

**Proposition 6.3.4.** For $\alpha(t) = \alpha^{*} \in (\alpha_{\min}^{*}, 1)$ and initial conditions (IC), system (S') has a lower bound

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

with

$$ \begin{aligned} 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{aligned} \quad (6.10) $$

*Proof.* For a trajectory emanating from $\mathcal{R}^{-}$ (respectively, $\mathcal{R}^{+}$), it follows that $\dot{r} \ge 0$ (respectively, $\dot{r} \le 0$) for all $t$ (according to Lemma 6.3.3), and so $r \ge r(0)$ (respectively, $r \ge r^{*}$) for all $t$. This proves that $r \ge \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} \ge 0$ (respectively, $\dot{m} \le 0$) for all $t$, and so $m \ge m(0)$ (respectively, $m \ge m^{*}$) for all $t$. Then, it follows that $m \ge \min(m(0), m^{*})$ for all $t \ge 0$. The equation for $p$ can thus be lower-bounded to

$$ \dot{p} \ge E_{M}m_{\text{low}} - (p+1)w_{R}(p)(r_{\max} - r_{\min}), $$

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

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) \ge w_R(p_{\text{low}})(r_{\text{low}} - r_{\text{min}}) > 0$ for all $t \ge 0$. Then, the global stability of the system is straightforward.

**Theorem 6.3.2.** 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 \ge p_{\text{low}} > 0$ and $r \ge r_{\text{low}} > r_{\min}$ for all $t \ge 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 6.3.5.** 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.

### 6.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 ([30] 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 the static optimization problem as

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

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

$$ \alpha^*(p^*) = \frac{E_M + (p^* + 1)w_R(p^*)\alpha_{\min}^*}{E_M + (p^* + 1)w_R(p^*)}. \qquad (6.11) $$
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

$$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). \quad (6.12)$$

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

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

Thus, the optimal allocation parameter $\alpha^*$ is obtained by replacing $p_{\text{opt}}^*$ in (6.11), and the maximal static growth rate can be calculated using (6.12). From (6.13), 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 (6.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}}^*)$.

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

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

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$.
## 6.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 (6.1) corresponds to the inverse of the protein density, which is set to 0.003 [L g⁻<fnref n="1" />] based on [30]. According to [2], the ribosomal fraction of the proteome<fnref n="1" /> can vary between 6% and 55%. In more recent studies [98], 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 [5, 6, 2, 98] 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 [2]. As a result, we have $n$ measurements of form ($\tilde{\mu}_k, \tilde{r}_k$) which are assumed to follow a linear relation [2], as seen in Figure 6.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 [2, 23, 98]. 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 (6.4). Based on [30], we fix the half-saturation constant of protein synthesis $K_R = 1$ g L⁻<fnref n="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 (6.6), the

<footnote n="1">¹The proteome is the total amount of protein in the cell</footnote>
<figure><img src="image_36.png" /><figcaption>Figure 6.2: Experimental data from [2, 5, 6] 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>

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^*_{opt}$ is calculated using (6.12), and the optimal steady state $(r^*_{opt}, p^*_{opt}, m^*_{opt})$ is expressed in terms of $\alpha^*_{opt}$ (using Theorem 6.3.1) 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 6.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 [6].
# 6.5 Optimal resource allocation

## 6.5.1 Problem definition

In this section we formulate the dynamic optimization problem under the hypothesis that microbial populations have evolved resource allocation strategies enabling them to maximize their biomass [13, 100]. 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

$$ \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. \qquad (\text{OCP}) $$

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 6.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 [30]: 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.

### 6.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 [53]. 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

$$ \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. \qquad (\text{PMP}) $$

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

$$ H(\varphi, \lambda, \lambda^0, \alpha) = \lambda^0 w_R(p)(r - r_{\min}) + \langle \lambda, F(\varphi, u) \rangle, \qquad (6.14) $$

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 (6.14) 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

$$\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.
\qquad (6.15)$$

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

$$\begin{aligned}
H_0 &= (E_M m - (p+1)w_R(p)(r-r_{\min}))\lambda_p - rw_R(p)(r-r_{\min})\lambda_r \\
&\quad + (r_{\max}-m)w_R(p)(r-r_{\min})\lambda_m - \lambda^0 w_R(p)(r-r_{\min}), \\
H_1 &= r_{\max}w_R(p)(r-r_{\min})(\lambda_r - \lambda_m).
\end{aligned}
\qquad (6.16)$$

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 [AYT5, AYT2], a bang arc $\alpha = 0$ (respec-
tively $\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 6.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 (6.15). 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$. $\square$

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

## 6.5.3 Study of the singular arcs

### 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. Differ-entiating along an extremal $(\varphi, \lambda)$ amounts to taking a Poisson bracket<fnref n="2" /> with the Hamiltonian $H$ [53]. 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 \},$$

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

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

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 [69], the so-called Fuller phenomenon. Otherwise, the order is termed local, and Fuller phenomenon may or may not occur. Using (6.17), 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

$$\{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>
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}, \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 6.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. $\square$
Going back to the three-dimensional formulation, one can explicit the computations by successively differentiating the expression (6.16).

### 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 6.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 \doteq 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. $\square$

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

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

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

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

Along a singular arc, the Hamiltonian can be rewritten as

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

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 6.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 (6.18) 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 (6.14) 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 6.5.3 and $\lambda^0 \neq 0$. $\square$

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 6.5.4 allows us to reduce the expression to

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

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

**Lemma 6.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 (6.13), 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

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

We define the function

$$y(p) \doteq w_R(p) (r_{\max} - x(p) + (p + 1)x'(p)). \quad (6.19)$$
Using (Condition 3) and (6.19) 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 6.5.6.** Along a singular arc over the interval $[t_1, t_2]$,

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

$$ r = x(p) - (x(p) - r_{\min}) \frac{y(p)}{E_M x'(p) + y(p)}. \qquad (6.21) $$

We differentiate (Condition 4) and we get

$$ -(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, \qquad (6.22) $$

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 (6.22) is $-r_{\max}w_R(p)(r - r_{\min})$, which cannot vanish as proven in Proposition 6.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 [69]. 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 6.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 6.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.

## The turnpike phenomenon

Using (6.20) and (6.21), 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

$$ \operatorname{sign}(\dot{p}) = \operatorname{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 6.5.1.** 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 [71]. This phenomenon has been considerably studied in econometry [101], and more recently in biology [102, 30, 3]. 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.

## 6.5.4 Numerical results

The computations of the optimal trajectories were performed with Bocop [4], 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 6.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 6.5.1). 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 6.5.1). In order to verify the optimality of the singular arc, we performed a numerical

<figure><img src="image_37.png" /><figcaption>Figure 6.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>

<footnote n="3">³https://ct.gitlabpages.inria.fr/gallery/bacteria/bacteria.html</footnote>
computation of the derivatives of $H_1$, which is shown in Figure 6.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

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

along the singular arc, which is a necessary condition for optimality. As we state in [AYT2], 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 6.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 [103], 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.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 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.
<figure><img src="image_38.png" /><figcaption>Figure 6.4: Factors of $\alpha$ in the derivatives of $H_1$ evaluated over the trajectory plotted in Figure 6.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 (6.23).</figcaption></figure>

### 6.6.1 Nutrient shift

Bacteria are known to traverse different habitats throughout their lifetime, experiencing fluctuating nutrient concentrations in the medium. In [30], 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.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.

*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 (6.5), we obtain

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

Along the singular arc, it holds that $m + r = x(p)$, which, using (6.24), 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. $\square$

A numerical simulation of this scenario is shown in Figure 6.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 [2]. Additionally, while the presence of $r_{\min}$ does not noticeably affect the solution of the optimal control problem, it contributes to a model that more accurately reproduces the experimental data (Figure 6.2a), representing a significant improvement from the modeling point of view.

### 6.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 [104, 90]. Likewise, the production of other proteins is known to protect bacteria like E. coli against acid stress [105]. Another possible scenario is the response to metabolic load imposed by the induced
<figure><img src="image_39.png" /><figcaption>Figure 6.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.</figcaption></figure>

overexpression of a heterologous protein [106]. All of these situations are known to reduce the resources available for growth-associated proteins (Figure 6.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}$.

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}), $$
<figure><img src="image_40.png" /><figcaption>Figure 6.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>

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 6.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 [2].

## 6.7 Conclusion

In this work, we proposed a dynamical self-replicator model of bacterial growth based on the work of [30], 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* [5, 6, 2, 98]. The dynamic problem is approached through optimal control theory, by
<figure><img src="image_41.png" /><figcaption>Figure 6.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>

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 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. When tackling more complex problems as proposed, e.g., in Tsiantis & Banga [107], a PMP perspective tends to yield very complicated mathematical formulations. Using direct methods has the advantage of avoiding these issues, but it often requires some knowledge to initialize the optimization algorithm or to check the validity of the solutions. In order to investigate complex biological systems, we advocate the development and theoretical analysis of simple models, in line with the question to be investigated, coupled with numerical exploration of optimal solutions (using larger models if necessary).

## Acknowledgments

We would like to acknowledge the help of Sacha Psalmon and Baptiste Schall from Polytech Nice Sophia for the numerical simulations and the production of the online example in the control toolbox gallery.
# Chapter 7

## Conclusion

In this manuscript, we explored a series of problems related to the allocation of resources in mathematical models of bacterial growth, and its applications to the production of compounds of interest.

Chapter 2, 3 and 4 have been dedicated to the production of metabolites in the most relevant cost-effective operating modes used in industry for commercial production. The main difference in these approaches lies in how the rich medium is supplied to the bioreactor and how the end product is extracted, which greatly affects the way the bacterial culture interacts with the environment. In every case, we studied the dynamical behavior of a self-replicator system representing the resource allocation problem in bacteria. An ideal case with constant substrate is analyzed in Chapter 2, representing an environment where bacteria sense no changes in the medium. In Chapter 3, the nutrient is assumed to be consumed by the culture until there is no available substrate in the bioreactor, which represents a batch process. Finally, the continuous bioreactor case is analyzed in Chapter 4, where there are both an inflow and an outflow producing—under certain conditions—a steady-state of the system, which is a well-studied behavior in this framework. The main trade-off intervening in these production schemes arise from the competition between the natural objective of the bacteria that aims to maximize the growth rate, and the biotechnological challenge of maximizing the synthesis of the compound of interest.

The optimal control problems arising from these works are tackled analytically and numerically, and can be considerably instructive in how to cope with such trade-offs in order to accomplish the proposed biosynthetic objectives. However, they frequently yield open-loop allocation strategies, which are often difficult to
implement. This issue is dealt with in Chapter 5, where a model-predictive control loop is proposed. The closed-loop implementation is based on the optimal open-loop strategies found in previous works, and offers a potential alternative to be explored.

Finally, a resource allocation model that considers growth rate-independent proteins is presented in Chapter 6. The latter represents a generalization of the model introduced in [30], and aims to study how bacteria can adapt to changing environments (such as nutrient upshifts and stress) from a dynamical perspective.

The models introduced here are simplified representations of very complex biological systems, which take into account two main cellular processes: metabolism and gene expression. The simplicity of the models allows a qualitative analysis of specific biological and biotechnological problems, and of the trade-offs arising from each objective. By design, several phenomena are not considered in the models, such as protein degradation and cell division. Additionally, the self-replicator systems are intended to represent an average composition of the cells forming the culture, but, in practice, resource distribution in individual cells is not necessarily homogeneous across the population. A major advantage of this methodology is that it allows to obtain analytical results to the proposed mathematical problems, in contrast to other more comprehensive approaches that aim to integrate numerous cellular mechanisms in a single system, yielding rather complex models.

There are numerous aspects to be explored in the future. All of the models are calibrated using data of previous steady-state experiments (or parameters in the literature obtained under the same principle). However, there is a lack of experimental results linking the trajectories of the dynamical models to the real-life bacterial behaviors. As explained in the introduction of this manuscript, a series of experiments with the *E. coli* bacteria were conducted throughout my PhD, with the objective of validating the theoretical results. The experiments were carried out in the context of the ANR Maximic project<fnref n="1" /> with the Ibis team (Inria Grenoble - Rhône-Alpes) with H. de Jong, J. Geiselmann and E. Cinquemani. However, experimentation in the microbiological framework requires multiple iterations in order to consider the obtained data as representative of the studied phenomenon, which was not possible to accomplish due to lack of time. Yet, it remains an important follow-up of this project.

Another potential extension of these works is the inclusion of additional bio-

<footnote n="1">¹https://project.inria.fr/maximic/</footnote>
logical mechanisms and constraints into the bacterial models, and the study of the impact of such modifications in the resulting optimal allocation strategies. This was the objective of the work described in Chapter 6, where the presence of growth rate-independent proteins (representing, for example, stress response proteins) imposes an upper bound on the ribosomal mass fraction of the cell. An important mechanism to contemplate in such models could be the energy metabolism, known to play an important limiting role in protein synthesis, which could refine the results previously obtained. Similarly, changing the optimization criteria can yield very different results in the distribution of resources to different cellular functions. The hypothesis of growth-gate maximization has been validated in many experiments of *E. coli*, and proven to be a selective advantage in certain environments. However, such assumptions are context-dependent, and could lose validity under different environmental conditions.

A modelling assumption that can be reviewed is the action of the external control *I*, modelled in (2.1) (and in the rest of the manuscript) through a multiplicative term affecting directly the natural allocation process. The idea behind this modelling choice is that the external control can arrest the production of RNA polymerase. However, an equally valid hypothesis is that the control can stall general protein synthesis, including enzymes, which would drastically change the control problem. In fact, the latter is part of an ongoing project extending the work on MPC loops presented in this manuscript. The approach is also complemented by a stage of state and parameter estimation, towards an implementation of a real-time bacterial resource allocation controller.

Finally, we recall that the work presented in Chapter 4 relies on a simplifying assumption linking metabolite production and protein synthesis through a certain function of the ribosomal concentration. This hypothesis provides an additional mass conservation law allowing to study the original dynamical system through a lower dimensional limiting system. However, the general case remains unstudied, as even the analysis of the local equilibria through the Routh-Hurwitz stability criterion represents a challenge for such system. This motivated an alternative approach to the problem, in which the stability of the local equilibria is examined by means of real algebraic tools, in collaboration with M. Safey El Din (PolSys team, LIP6, Sorbonne Université).

Manipulation of bacterial cultures for the production of added-value compounds
remains a complex endeavour. Hopefully, theoretical research on the field, boosted
by our growing understanding of microorganisms, will progressively improve the
existing biotechnological techniques towards more efficient and sustainable industrial
schemes.
# Appendix A

## Time-optimal control of piecewise affine bistable gene-regulatory networks

This chapter reproduces [AYO1] and [AYO2], written in collaboration with Nicolas Augier, and accepted for publication in the International Journal of Robust and Nonlinear Control and the 7th IFAC Conference on Analysis and Design of Hybrid Systems respectively. While these works are not in line with the topic of resource allocation in bacterial growth, they have been included in the Appendix as they were produced during my PhD as a result of a joint work. In accordance with the usages in the mathematics field, the authors of this publications are in alphabetical order, regardless of the fact that we equally contributed to the works.

### A.1 Introduction

Understanding complex biological phenomena has become of great interest in the last decades for the scientific community. In the context of synthetic biology, numerous fields of study have been employed to better comprehend and re-engineer the interactions within biological systems [108]. Such is the case of control theory, widely used to explain regulatory mechanisms in nature [30, 109], but also to artificially act upon them for biotechnological purposes [AYT5, AYT2]. A classical example is the metabolism of cells, described by multiple regulatory mechanisms forming complex networks. In this framework, the interaction between genes is a
crucial subject of study [110], whose typical behaviors can be described by positive
and negative feedback loops [111, 112, 113, 114, 115]. The dynamics of these loops
have been extensively analyzed, both from experimental and theoretical perspec-
tives, and are known to present either multistability or oscillatory behaviors. From
a mathematical modeling perspective, such a systems can be modeled through dy-
namical systems of several variables, where the positivity or negativity is given by
the parity of negative interactions forming the loop.

Among all existing patterns, the simplest positive feedback loop is the two-
dimensional bistable system, which is commonly used to represent the so-called
*genetic toggle switch*. The latter is a synthetic flip-flop device first implemented
experimentally in *E. coli* through the genes *lacI* and *tetR* mutually repressing each
other [116]. The state of the device is determined by the concentration of the genes
in the boolean form (low, high) and (high, low). This allows genetic toggle switches
to act as biological memory units capable of storing 1 bit of information, by sus-
taining one of the two possible states through time [117], which offers a biosynthetic
alternative to the classical electronic flip-flop. Since its creation, understanding how
to regulate bistable systems in a reliable manner (e.g. by suppressing undesirable
oscillations [118] or achieving transitions between states [AYO2] [119, 120]) has be-
come highly relevant for their vast implications in biotechnology and biocomputing.

In practice, the state of a genetic toggle switch can be controlled by externally
catalyzing or inhibiting the synthesis rate of the genes. This is done by introduc-
tion of a plasmid, which are essentially small circular DNA molecules that can be
constructed to include an inducible promoter of the studied gene, thus affecting the
synthesis rate of messenger RNA. Thus, the transcription rate can be directly mod-
ified by aggregation of an inducer. In *E. coli*, this is done by externally adding the
diffusible molecules IPTG<fnref n="1" /> and aTc<fnref n="2" />, which are known to repress the *lacI* and *tetR*
genes, respectively [116].

Motivated by this experimental scheme, some authors proposed exact control
strategies based on a piecewise affine model of the bistable system [120]. As discussed
in the work, the importance of studying the phenomenon through qualitative models
arises from the constraints related to the experimental setup, both in measuring the
state and in acting on the system. The proposed model is characterized by the

<footnote n="1">¹isopropyl-β-D-thiogalactopyranoside</footnote>

<footnote n="2">²anhydrotetracycline</footnote>
existence of an "undifferentiated state", where no gene is predominant, and from
which the system can evolve towards one of the two attractors. Mathematically,
the instability of this state appears as a Filippov non-smooth "saddle" singularity
[121]. From a biological point of view, such a state plays a key role in cell decision
making and cell fate differentiation [122]. Its role in fate commitment has also
motivated experimental studies aiming at stabilizing genetic toggle switches around
this undifferentiated point [123, 124].

Whereas most of the theoretical work in the subject has been dedicated to ex-
ternally producing state transfers [125], the time efficiency of state switches has
received little or no attention from the community. Indeed, one of the key issues in
these genetic devices is the time needed to induce a transfer between its two stable
states, due to its importance when studying more complex networks of systems in-
volving different time scales. In particular, the latter becomes a major constraint in
the framework of biological signal processing [126]. Genetic toggle switches operate
at the level of gene transcription and translation, whose duration and timescales
are the main factor delaying the availability of the proteins when facing a switch
between steady states. In this context, the minimization of a state switch, which is
directly linked to the production of the non-expressed protein, becomes highly rel-
evant. Recent works [127] showed the importance of accelerating transitions times
(and minimizing inducer usage) in artificially engineered bistable systems in order
to obtain less costly (and therefore, more sustainable) chemical production schemes.
Thus, in this paper, we investigate the time-optimal control strategies for the afore-
mentioned bistable system. Our aim is to induce transitions between the two stable
steady states in minimal time. Many complex systems are known to involve bistable
processes [128, 129]. Hence, the reduction of the time needed for such transfers could
allow experimentalists to speed up certain chemical reactions or to artificially in-
crease bacterial growth rate, thus improving yield in biotechnological processes. In
a different setting [130, 131], the time efficiency of bistable systems switches in two-
level quantum systems was studied, so as to induce efficient transitions between two
quantum states.

In addition to the biological relevancy of the subject, the resulting OCP (Optimal
Control Problem) yields very interesting results in the framework of Hybrid Optimal
Control. The steady states of the piecewise linear system cannot be reached in finite
time due to the lack of controllability in certain regions, and so one has to consider
a relaxed OCP with "partial targets", that is, driving a given protein to a certain fixed value larger than its corresponding threshold. In this regard, we show that time-optimal strategies for such a problem have a very specific geometric description. When the initial state is far enough from the target, that is below a curve called *separatrix*, we show by an adaptation of the HMP (Hybrid Maximum Principle) to our setting that the optimal control consists in a concatenation of two bang arcs, and the optimal trajectories follow:

* a first phase in which the system reaches the separatrix;

* a second phase where the system slides along this curve, until reaching the "undifferentiated" point of the biological system in finite time;

* a third phase, where the system leaves this curve, slides along a second fixed curve and reaches its target.

These two curves correspond to the stable and unstable manifolds of the undifferentiated saddle-type singularity, and the point where the dynamics achieves its transfer is nothing but the corresponding Filippov equilibrium. The latter behavior can be compared to the turnpike phenomenon [71], where the optimal trajectory for a given OCP for large final times is shown to remain close to a steady-state trajectory solution of the associated static OCP. Besides its specific interest, we expect our method to open new prospects in the study of optimal control of higher dimensional genetic regulatory networks. In particular, it often occurs that trajectories belonging to a given domain may bifurcate into different domains, similarly to what happens in the toggle switch case, and some similar turnpike-like properties may hold in this case.

The paper is organized as follows: in Section A.2, we present both the non-controlled system and the studied controlled system, and we provide some technical results. In Section A.3, we introduce the time-optimal control problem and we adapt the HMP to our setting. In Section A.4, we present the main results, that prove the qualitative features of the optimal trajectories mentioned above. In Section A.6, we give a lower bound for the minimal time, then characterizing the minimal transfer time between the two states of the toggle switch model. In Section A.7, we provide numerical results implemented with Bocop [4], an open-source toolbox for solving OCPs. Additionally, we perform a numerical comparison between the trajectories of the relaxed OCP and the original OCP, that suggests that the results also hold for the original one.
# A.2 Bistable-switch model

## A.2.1 Free dynamics

Consider two variables $x_1$ and $x_2$ which represent two genes mutually inhibiting each other. The individual dynamics, defined in Filippov sense, is the following

$$ \begin{cases} \dot{x}_1 = -\gamma_1 x_1 + k_1 s^-(x_2, \theta_2), \\ \dot{x}_2 = -\gamma_2 x_2 + k_2 s^-(x_1, \theta_1), \end{cases} \quad (A.1) $$

where for $j \in \{1, 2\}$, $x_j \in \mathbb{R}_+$, and for $\theta \in \mathbb{R}$, $s^-(\cdot, \theta) : \mathbb{R} \to \mathbb{R}$ is such that

$$ s^-(x, \theta) = \begin{cases} 1 & \text{if } x < \theta, \\ 0 & \text{if } x > \theta. \end{cases} $$

It is assumed that $s^-(x) \in [0, 1]$ for $x = \theta$. The positive constants $(\gamma_j)_{j \in \{1, 2\}}$, $(k_j)_{j \in \{1, 2\}}$ correspond, respectively, to the degradation and the production rates of each variable. It is classical [120] that the domain $K = [0, k_1/\gamma_1] \times [0, k_2/\gamma_2]$ is forward invariant by the dynamics of Equation (A.1). From now on, we consider only solutions evolving in $K$. Define the regular domains

$$ \begin{aligned} B_{00} &= \{(x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_1 < \theta_1, 0 < x_2 < \theta_2\}, \\ B_{01} &= \{(x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_1 < \theta_1, \theta_2 < x_2 < \frac{k_2}{\gamma_2}\}, \\ B_{10} &= \{(x_1, x_2) \in \mathbb{R}^2 \mid \theta_1 < x_1 < \frac{k_1}{\gamma_1}, 0 < x_2 < \theta_2\}, \\ B_{11} &= \{(x_1, x_2) \in \mathbb{R}^2 \mid \theta_1 < x_1 < \frac{k_1}{\gamma_1}, \theta_2 < x_2 < \frac{k_2}{\gamma_2}\}, \end{aligned} $$

which are defined as open sets in accordance with the HMP approach to be applied in Section A.3.3. Equation (A.1) restricted to a regular domain $B_{ij}$ is an affine dynamical system on $\mathbb{R}^2$ having an asymptotically stable equilibrium, called focal point for system (A.1). Each region $B_{ij}$ for $i, j \in \{0, 1\}$ has a focal point

$$ \phi_{ij} = (\bar{x}_i, \bar{x}_j) $$

corresponding to

$$ \bar{x}_i = \frac{k_i}{\gamma_i} s^-(\bar{x}_j, \theta_j). $$
Thus, system (A.1) has two locally asymptotically stable steady states

$$ \phi_{10} = \left( \frac{k_1}{\gamma_1}, 0 \right) \in \bar{B}_{10}, $$

$$ \phi_{01} = \left( 0, \frac{k_2}{\gamma_2} \right) \in \bar{B}_{01}, $$

and an unstable Filippov equilibrium point at $(\theta_1, \theta_2)$. Figure A.1 illustrates the dynamics of the system for a given set of parameters.

<figure><img src="image_42.png" /><figcaption>Figure A.1: Stream plot with free dynamics given by Equation (A.1). System parameters are $\gamma_1 = 1.1$, $\gamma_2 = 1.7$, $\theta_1 = 0.6$, $\theta_2 = 0.4$, $k_1 = k_2 = 1$.</figcaption></figure>

### A.2.2 Controlled dynamics and some related properties

We write the controlled dynamics assuming that the synthesis rates of each gene can be externally catalyzed or inhibited (e.g. through the introduction of inducible promoters of a given gene), as indicated in the previous section. Mathematically, this is represented by the control input $u$ acting directly on the synthesis rate of each gene, in a multiplicative form. Then, the controlled system, defined in Filippov
sense, is

$$ \left\{ \begin{array}{l} \dot{x}_1 = -\gamma_1 x_1 + u(t)k_1 s^-(x_2, \theta_2), \\ \dot{x}_2 = -\gamma_2 x_2 + u(t)k_2 s^-(x_1, \theta_1), \end{array} \right. \qquad (S) $$

where the control $u(\cdot) \in L^\infty([0, t_f], [u_{\min}, u_{\max}])$, with $0 < u_{\min} < 1 \le u_{\max}$. This system is motivated by the one introduced in [120], which assumed the same control input $u$ for the two variables. The latter aims to model a simple qualitative control, easier to implement in a molecular biology setting than the case with two distinct control variables. We make the following assumptions on the parameters of the system (for more details, see [120]).

**Assumption A.2.1.** The parameters $(\gamma_j)_j$ and $(k_j)_j$ satisfy

$$ \theta_j < \frac{k_j}{\gamma_j}, \quad j \in \{1, 2\}; \quad \frac{\theta_2}{\theta_1} > \frac{k_2 \gamma_1}{k_1 \gamma_2}; \quad \frac{\theta_2}{\theta_1} < \frac{k_2}{k_1}. $$

This assumption is based on intrinsic conditions of the parameters of the non-controlled system, and allows to find a control strategy driving the solution of Equation (S) from $B_{10}$ to $B_{01}$, as well as from $B_{01}$ to $B_{10}$. Note that it implies $\gamma_1 < \gamma_2$, and that the case where $\frac{\theta_2}{\theta_1} > \frac{k_2 \gamma_1}{k_1 \gamma_2}$ and $\frac{\theta_2}{\theta_1} < \frac{k_2}{k_1}$ can be treated analogously by permutation of $x_1$ and $x_2$.

## Separatrix

Now define the *separatrix*, which is a curve playing a fundamental role in the global dynamics of both the open-loop system (A.1) and the controlled system (S). For a fixed value of $u(t) \equiv u \in [u_{\min}, u_{\max}]$, the separatrix ($S_u$) is defined as the stable manifold of the Filippov equilibrium $(\theta_1, \theta_2)$ for Equation (S) restricted to $\overline{B_{00}}$. In the coordinates $(x_1, x_2) \in B_{00}$, for $u \ge 1$, the separatrix ($S_u$) can be written as the curve of equation

$$ x_2 = \alpha(x_1, u) = \frac{k_2 u}{\gamma_2} - \left( \frac{k_2 u}{\gamma_2} - \theta_2 \right) \left( \frac{\frac{k_1 u}{\gamma_1} - x_1}{\frac{k_1 u}{\gamma_1} - \theta_1} \right)^{\frac{\gamma_2}{\gamma_1}}. $$
Using the latter, we define the regions

$$ (S_u)^+ = \left\{ (x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_1 < \theta_1, \alpha(x_1, u) < x_2 < \frac{k_2}{\gamma_2} \right\}, $$

$$ (S_u)^- = \left\{ (x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_2 < \theta_2, \alpha(x_1, u) > x_2, x_1 < \frac{k_1}{\gamma_1} \right\}, $$

such that the domain $\overline{K}$ is divided into

$$ \overline{K} = \overline{(S_u)^+} \cup \overline{(S_u)^-} \cup \overline{B_{11}}, \qquad (A.2) $$

as shown in Figure A.2. The solutions of Equation (S) having initial conditions in

<figure><img src="image_43.png" /><figcaption>Figure A.2: Division of the domain $\overline{K}$ as defined in (A.2), with a vector field defined by a constant control $u < 1$.</figcaption></figure>

$(S_u)^-$ (respectively, $(S_u)^+$) reach $B_{10}$ (respectively, $B_{01}$) in finite time. Moreover, $B_{10}$ (respectively, $B_{10}$) is included in the bassin of attraction of $\phi_{10}$ (respectively, $\phi_{01}$). Notice that, for a fixed value of $u(t) \equiv u \in [u_{\min}, u_{\max}]$, the solutions of Equation (S) having initial conditions in $(S_u)$ reach the Filippov point $(\theta_1, \theta_2)$ in finite time. Once having reached this point, the solution of Equation (S) is then defined by differential inclusion in the Filippov sense. Roughly speaking, there exist
several solutions that will reach either $B_{01}$ or $B_{10}$ (see [120, Appendix] for more precise informations about Filippov solutions of such a system).

## Lower separatrix

Now we define the lower separatrix, which will be useful in Section A.6.

**Definition A.2.2.** For $(\gamma_j)_{j \in \{1,2\}}$ and $(k_j)_{j \in \{1,2\}}$ satisfying Assumption A.2.1, define the lower separatrix $(\tilde{S}_u)$ as the straight line of equation

$$x_2 = \beta(x_1, u) = \frac{k_2 u}{\gamma_1} - \left(\frac{k_2 u}{\gamma_1} - \theta_2\right) \left(\frac{\frac{k_1 u}{\gamma_1} - x_1}{\frac{k_1 u}{\gamma_1} - \theta_1}\right).$$

**Lemma A.2.3.** Let $(\gamma_j)_{j \in \{1,2\}}$ and $(k_j)_{j \in \{1,2\}}$ satisfy Assumption A.2.1 and $u \ge 1$, and let $\bar{x}_1$ be the unique $x_1 \in [0, \theta_1]$ such that $\beta(x_1, u) = 0$. Then for every $x_1 \in [\bar{x}_1, \theta_1]$ and $u \ge 1$, we have $\beta(x_1, u) \le \alpha(x_1, u)$, that is, $(\tilde{S}_u)$ is below $(S_u)$ for every $x_1 \in [\bar{x}_1, \theta_1]$.

*Proof.* For every $x_1 \in [0, \theta_1]$ and $u \ge 1$, define $X(x_1) = \frac{\frac{k_1 u}{\gamma_1} - x_1}{\frac{k_1 u}{\gamma_1} - \theta_1}$. We have easily that $\beta(x_1, u) = 0$ if and only if $x_1 = \bar{x}_1$, where $\bar{x}_1 \in [0, \theta_1]$ is such that $X(\bar{x}_1) = \frac{k_2 u/\gamma_1}{k_2 u/\gamma_1 - \theta_2}$. Evaluating in the expression of $\alpha$ at $x_1 = \bar{x}_1$, we have

$$\alpha(\bar{x}_1, u) = \frac{k_2 u}{\gamma_2} - \left(\frac{k_2 u}{\gamma_2} - \theta_2\right) \left(\frac{k_2 u/\gamma_1}{k_2 u/\gamma_1 - \theta_2}\right)^{\frac{\gamma_2}{\gamma_1}}$$

Using the fact that $\gamma_2 > \gamma_1$ under Assumption A.2.2, one can prove by a direct differentiation that the function

$$x \mapsto \frac{\left(\frac{k_2 u}{\gamma_1} - x\right)^{\frac{\gamma_2}{\gamma_1}}}{\frac{k_2 u}{\gamma_2} - x}$$

is non-decreasing on $[0, \theta_2]$, and that $x_1 \mapsto \alpha(x_1, u)$ is concave for $x_1 \in [0, \theta_1]$. Hence, we have

$$\left( \frac{\frac{k_2 u}{\gamma_1}}{\frac{k_2 u}{\gamma_1} - \theta_2} \right)^{\frac{\gamma_2}{\gamma_1}} \le \frac{\frac{k_2 u}{\gamma_2}}{\frac{k_2 u}{\gamma_2} - \theta_2},$$
and we can deduce $\alpha(\bar{x}_1, u) \ge 0$. Provided that $\alpha(\theta_1, u) = \beta(\theta_1, u) = \theta_2$, we deduce that $\alpha(x_1, u) \ge \beta(x_1, u)$, for every $x_1 \in [\bar{x}_1, \theta_1]$. $\square$

## A.3 Time-optimal transfer

### A.3.1 Problem formulation

The state of a genetic toggle switch is determined by gene expression in the boolean form (low, high) and (high, low), and so the objective in this work is to achieve a transition from one boolean state to the other in minimal time. In the mathematical context, the latter translates into finding trajectories that drive the solution $(x_1(t), x_2(t))$ of Equation (S) towards the steady states $\phi_{01}$ and $\phi_{10}$ of Equation (A.1) in minimum time (where these states correspond to the differentiated states aforementioned). However, due to the lack of controllability in direction $x_1$ (respectively, $x_2$) of Equation (S) restricted to $B_{01}$ (respectively, $B_{10}$), one has to relax the problem. More precisely, the steady state $\phi_{01}$ (respectively, $\phi_{10}$) cannot be reached in finite time, because $u$ does not act on $x_1$ in the domain $B_{01}$ (respectively, $x_2$ in the domain $B_{10}$). Thus, we will be first interested in driving $x_2(t)$ towards an arbitrary value $x_2(t_f) = x_2^f > \theta_2$ (for instance, the value $x_2^f = k_2/\gamma_2$ corresponding to the $x_2$-component of the steady state $\phi_{01}$), with the constraint that at the final time, $x_1(t_f)$ belongs to the interval $[0, \theta_1)$. This target choice ensures that, at the final time, the gene $x_2$ is strongly expressed while the gene $x_1$ is weakly expressed. The symmetric problem, which is equivalent, consists in driving $x_1(t)$ towards an arbitrary value $x_1(t_f) = x_1^f > \theta_1$, with the constraint that at the final time, $x_2(t_f)$ belongs to the interval $[0, \theta_2)$. Fix $x_1^0 \ge \theta_1$, $x_2^0 \le \theta_2$, $x_2^f \ge \theta_2$, and consider the minimization problem

$$ \left\{ \begin{array}{l} \text{minimize } t_f \ge 0, \\ x(t) = (x_1(t), x_2(t)) \text{ is subject to (S)}, \\ x(0) = (x_1^0, x_2^0), \\ x_2(t_f) = x_2^f, \\ x_1(t_f) \in [0, \theta_1), \\ u(\cdot) \in [u_{\min}, u_{\max}]. \end{array} \right. \qquad (OCP) $$
### A.3.2 Reachability of the terminal state

A fundamental aspect of OCPs with fixed terminal state is the existence of a solution. Such a matter is directly linked to the reachability and controllability analysis of the dynamical system, which are often hard to conduct analytically. In this work, we provide sufficient conditions for the feasibility of the proposed trajectory, and we show that a simple piecewise constant control strategy [120] achieves the objective, serving as a candidate to (OCP). This strategy drives asymptotically the system from an initial state in $B_{01}$ to $\phi_{10}$ (or $B_{10}$ to $\phi_{01}$ for the symmetric problem). One can show that, under Assumption A.2.1, there exist $u_{min} < \theta_1\gamma_1/k_1$ and $u_{max} \ge 1$ such that

$$ \Phi^*(u_{min}) \in (S_{u_{max}})^+, $$

with

$$ \Phi^*(u_{min}) \doteq \left( \frac{u_{min}k_1}{\gamma_1}, \frac{u_{min}k_2}{\gamma_2} \right). \qquad (A.3) $$

In previous works [120, Section 3], authors proved the existence of $\bar{u}_{min}$, $\bar{u}_{max}$ such that for every $u_{min}$, $u_{max}$ such that $0 \le u_{min} \le \bar{u}_{min} \le \bar{u}_{max} \le u_{max}$, we have $\Phi^*(u_{min}) \in (S_{u_{max}})^+$. The latter condition is satisfied for $\bar{u}_{min}$, $\bar{u}_{max}$ when there exists $\delta < \frac{\theta_1 k_2}{\theta_2 k_1} - 1$ and $\epsilon > 0$ small enough such that

$$ \bar{u}_{max} > \max \left\{ 1, \frac{(1 + \delta)\gamma_2 \frac{\theta_2}{k_2} - \gamma_1 \frac{\theta_1}{k_1}}{\delta} \right\}, $$

and

$$ 0 < \bar{u}_{min} < \min \left\{ \frac{\gamma_1}{k_1}\theta_1, \frac{\gamma_2}{k_2}\theta_2, (1 - \epsilon)x_1^* \frac{\gamma_1}{k_1} \right\}, $$

where $x_1^* \in (0, \theta_1)$ is the unique solution of $\alpha(x_1^*, u_{max}) = 0$. Concerning the symmetric problem, one can show the existence of another choice of $u_{min} < \theta_1\gamma_1/k_1$ and $u_{max} \ge 1$ such that $\Phi^*(u_{min}) \in (S_{u_{max}})^-$. Due to the symmetry of both problems, we will focus on the first case, and state the following assumption.

**Assumption A.3.1.** Bounds $u_{min}$ and $u_{max}$ are chosen such that $0 \le u_{min} \le \bar{u}_{min} \le \bar{u}_{max}$.
$\bar{u}_{max} \leq u_{max}$, so that $\Phi^*(u_{min}) \in (S_{u_{max}})^+$.

Based on the solution developed in previous works [120], we first propose an input control constrained to two possible values {$u_{min}, u_{max}$} corresponding to the low and high synthesis control. The control law is expressed in terms of the state and time as

$$ u(x,t) = \begin{cases} u_{min} & \text{if } x \in B_{10}, \\ u_{min} & \text{if } t \in [0, t_s), x \in B_{00}, \\ u_{max} & \text{if } t \in [t_s, \infty), x \in B_{00}, \\ u_{max} & \text{if } x \in B_{01}. \end{cases} \quad (A.4) $$

for $t_s > 0$ sufficiently large. During the first phase with $u \equiv u_{min}$, every focal point of the system belongs to $B_{00}$, hence the solution $x(t)$ of Equation (S) converges towards the point $\Phi^*(u_{min}) \in B_{00}$ when $t \to \infty$. During the second phase with $u \equiv u_{max}$, state $x(t)$ reaches $B_{01}$ in finite time, and $x_2(t)$ converges towards $x_2^f$ in finite time. From that point, an open-loop control $u \equiv 1$ drives $x(t)$ to $\phi_{01}$ when $t \to \infty$. An example illustrating this trajectory is shown in Figure A.3, where $x_2^f = k_2/\gamma_2$, matching the coordinate $x_2$ of the point $\phi_{01}$. Indeed, under Assumptions A.2.1 and A.3.1, and by choosing $t_s$ sufficiently large, the control strategy (A.4) ensures that any trajectory starting from $(x_1^0, x_2^0)$ reaches a final point meeting $x_1 \in [0, \theta_1)$ and $x_2 = x_2^f$ in finite time, which shows that the set of admissible controllers for problem (OCP) is non empty.

Notice that, while the latter strategy serves as a candidate, the set of possible controllers is not limited to bang-bang solutions. One could consider, for instance, non bang-bang strategies based on A.4 with intermediate control values (e.g. replacing $u_{max}$ by $\tilde{u}_{max} < u_{max}$) which also achieve the transfer. This implies that there exist several trajectories reaching the target of Problem OCP, which motivates a study from a Pontryagin's Maximum Principle perspective.

### A.3.3 Hybrid optimal control problem with a fixed domain sequence

Consider two compact subsets $M_0$ and $M_1$ of $\mathbb{R}^2$, and assume $M_1$ is reachable from $M_0$ for system (S), that is, such that there exists a time $t_f > 0$, a control $u(\cdot) \in L^\infty([0, t_f], \Omega)$ and $x_0 \in M_0$ such that the solution $x(t)$ of Equation (S), defined in
<figure><img src="image_44.png" /><figcaption>Figure A.3: Optimal trajectory with $x_1^0 = 0.8$, $x_2^0 = 0.3$ and $x_2^f = k_2/\gamma_2$. System parameters are $\gamma_1 = 1.1$, $\gamma_2 = 1.7$, $\theta_1 = 0.6$, $\theta_2 = 0.4$, and $k_1 = k_2 = 1$. Control bounds are set to $u_{min} = 0.4$ and $u_{max} = 1.1$. The control switches from $u \equiv u_{min}$ to $u \equiv u_{max}$ at time $t_s$ (= 3 in this case), after the state $x(t)$ has crossed the separatrix ($S_{u_{max}}$).</figcaption></figure>

the Filippov sense with initial condition $x(0) = x_0$ satisfies $x(t_f) \in M_1$. Consider the problem of steering the system (S) from $M_0$ to $M_1$ in minimal time $t_f$. In order to properly define the problem, one has to choose a sequence $B$ in the set $\{B_{00}, B_{01}, B_{10}, B_{11}\}$ of regular domains, and consider $B$-admissible trajectories of Equation (S), defined as follows.

**Definition A.3.2.** Let $B = (B_j)_{j \in \{1, \dots, k\}}$ be a sequence of regular domains. We say that a solution $x(t)$ of Equation (S) is $B$-admissible if there exists a time $T > 0$, a control $u(\cdot) \in L^\infty([0, t_f], \Omega)$, and times $t_0 = 0 < t_1 < \dots < t_k$ such that $x(t) \in B_j$ for every $t \in \Delta_j$, where $\Delta_j = (t_j, t_{j+1})$.

In particular, the previous definition excludes sliding modes along the frontier between two successive regular domains. Additionally, we require two more assumptions related to the reachability of $B$-admissible solutions for the general case.

**Assumption A.3.3.** $M_0$ (respectively, $M_1$) is included in the adherence $\bar{B}_{jk}$ (respectively, $\bar{B}_{qi}$) of a regular domain, for $j, k, q, i \in \{0, 1\}$.

**Assumption A.3.4.** Assume that there exists a time $T > 0$ and a $B$-admissible solution $(x(t), u(t))$ such that $x(0) \in M_0$ and $x(T) \in M_1$.
Notice that for given sets $M_0, M_1$, the choice of the sequence $B$ is not unique in general. Assume that $M_0, M_1, B$ satisfy the assumptions A.3.3 and A.3.4. For a fixed sequence $B = (B_j)_{j \in \{1, \dots, k\}}$, we can consider Problem (A.4.1) restricted to $B$-admissible trajectories. Necessary conditions of optimality for this problem can be directly derived from the HMP, which we will state in Theorem A.3.1 of the following subsection.

### A.3.4 Hybrid Maximum Principle for time optimal control

In this section, we provide an adaptation of the Hybrid Maximum Principle given by Dmitruk and Kaganovich [132] to the time optimal setting. Let $t_0 < t_1 < \dots < t_\nu$ be real numbers. Denote by $\Delta_k$ the time interval $[t_{k-1}, t_k]$. For continuous functions $x^k : [t_0, t_\nu] \to \mathbb{R}^n$, $k \in \{1, \dots, \nu\}$, define the vector

$$p = (t_0, (t_1, x^1(t_0), x^1(t_1)), \dots, (t_\nu, x^\nu(t_{\nu-1}), x^\nu(t_\nu))) \in \mathbb{R}^d,$$

where $d = 1 + (2n + 1)\nu$. Let $(f_k)_{k \in \{1, \dots, \nu\}}$ be smooth vector fields on $\mathbb{R}^n$, and $(\phi_i)_{i \in \{1, \dots, m\}}$, $(\eta_j)_{j \in \{1, \dots, q\}}$ be two families of smooth functions defined on $\mathbb{R}^{(\nu+1)(n+1)}$. For $t \in [t_0, t_\nu]$ and a collection $(U_k)_{k \in \{1, \dots, \nu\}}$ of subsets of $\mathbb{R}^q$, $q \ge 1$, consider the autonomous hybrid OCP

$$
\left\{
\begin{array}{l}
\displaystyle \text{minimize } t_{\nu} - t_{0}, \\
\displaystyle \dot{x}^{k}(t) = f_{k}(x^{k}(t), u^{k}(t)), \\
\displaystyle u^{k}(t) \in U_{k}, \\
\displaystyle t \in \Delta_{k}, \\
\displaystyle k \in \{1, \dots, \nu\}, \\
\displaystyle \eta_{j}(p) = 0, \quad j = 1, \dots, q, \\
\displaystyle \phi_{i}(p) \le 0, \quad i = 1, \dots, m.
\end{array}
\right.
\qquad (\text{HOCP})
$$

**Definition A.3.5.** For a tuple $w = (t_0; t_k, x^k(t), u^k(t), k = 1, \dots, \nu)$ which is extremal for Problem (HOCP), define:

* the trajectory $(x(t))_{t \in [t_0, t_\nu]}$ which is equal to $x^k(t)$ for every $t \in \Delta_k \setminus \{t_k\}$ and $k \in \{1, \dots, \nu\}$;
* the adjoint trajectory $(\lambda(t))_{t \in [t_0, t_\nu]}$ which is equal to $\lambda^k(t)$ for every $t \in \Delta_k \setminus \{t_k\}$ and $k \in \{1, \dots, \nu\}$;

* the control $(u(t))_{t \in [t_0, t_\nu]}$ which is equal to $u^k(t)$ for every $t \in \Delta_k \setminus \{t_{k-1}\}$ and $k \in \{1, \dots, \nu\}$.

Define, for every $k \in \{1, \dots, \nu\}$, $t \in \Delta_k$,

$$H^k(x^k, \lambda^k, \lambda_0, u^k) = \langle \lambda^k, f_k(x^k, u^k) \rangle - \lambda_0.$$

**Theorem A.3.1.** Assume that $(\tilde{x}(\cdot), \tilde{u}(\cdot), \tilde{p})$ is an optimal solution of Problem (HOCP). Then there exists

$$(\alpha, \beta, \lambda(\cdot), \lambda_0),$$

where $\alpha = (\alpha_1, \dots, \alpha_m) \in \mathbb{R}^m$, $\beta = (\beta_1, \dots, \beta_q) \in \mathbb{R}^q$, $\lambda = (\lambda^1, \dots, \lambda^\nu)$, all $\lambda^k$: $\Delta_k \to \mathbb{R}^n$ for $k \in \{1, \dots, \nu\}$ being Lipschitz functions, and a constant $\lambda_0 \ge 0$ such that:

* $(\lambda_0, \alpha, \beta) \neq 0$;

* For every $i \in \{1, \dots, m\}$, $\alpha_i \ge 0$;

* For every $i \in \{1, \dots, m\}$, $\alpha_i \phi_i(\tilde{p}) = 0$;

* For almost every $t \in \Delta_k$,

$$ \begin{aligned} \dot{x}^k &= \frac{\partial H^k}{\partial \lambda}(x^k, \lambda^k, \lambda_0, \tilde{u}), \\ \dot{\lambda}^k &= -\frac{\partial H^k}{\partial x^k}(x^k, \lambda^k, \lambda_0, \tilde{u}), \\ H^k(x^k, \lambda^k, \lambda_0, \tilde{u}) &= \max_{u \in \Omega} H^k(x^k, \lambda^k, \lambda_0, u) = 0. \end{aligned} \tag{E} $$

Moreover, if we define $L(p) = \lambda_0(t_\nu - t_0) + \sum_{i=1}^m \alpha_i \phi_i(p) + \sum_{j=1}^q \beta_j \eta_j(p)$, then we have the following transversality and discontinuity conditions at times $t = t_0, \dots, t_\nu$:

* At the initial and final times $t_0$ and $t_\nu$, we have

$$\lambda^1(t_0) = \frac{\partial L}{\partial x^1(t_0)}(\tilde{p}),$$

$$\lambda^\nu(t_\nu) = \frac{\partial L}{\partial x^\nu(t_\nu)}(\tilde{p}).$$
• At the crossing times $(t_k)_{k \in \{1, \dots, \nu-1\}}$, we have, for every $k \in \{1, \dots, \nu-1\}$,

$$ \lambda^k(t_{k-1}) = \frac{\partial L}{\partial x^k(t_{k-1})}(\tilde{p}), $$

$$ \lambda^k(t_k) = -\frac{\partial L}{\partial x^k(t_k)}(\tilde{p}). $$

## A.4 Main results

We are interested in solving (OCP) among continuous B-admissible trajectories, as defined in Section A.3.3. We first observe that the regular domain $B_{11}$ is repulsive, and so any B-admissible trajectory with $x(0) \in B_{10}$ and $x(t_f) \in B_{01}$ should pass through $B_{00}$, as the point $(\theta_1, \theta_2)$ cannot be reached from $B_{10}$. Thus, we fix the sequence of regular domains $B = (B_{10}, B_{00}, B_{01})$, with $M_0$ restricted to a point in $K$, and $M_1 = \{(x_1, x_2) \in K \mid x_1 \in [0, \theta_1), x_2 = x_2^f\}$ which has already been proven to be reachable in finite time, verifying assumptions A.3.3 and A.3.4. As previously said, the problem can be further analyzed by applying HMP. The Maximum Principle in the Hybrid framework requires to define functions $(\phi_i)_i$ and $(\eta_j)_j$ that guarantee the continuity of the trajectories and the changes of dynamics at the frontiers $x_1 = \theta_1$ and $x_2 = \theta_2$ [132]. Through its application, we obtain that (OCP) admits an optimal control which can be defined as a very simple feedback.

**Theorem A.4.1.** The optimal strategy $u(x)$ solution of (OCP) for B-admissible trajectories is the feedback control

$$ u(x) = \begin{cases} u_{min} & \text{if } x \in (S_{u_{max}})^-, \\ u_{max} & \text{if } x \in (S_{u_{max}})^+ \cup (S_{u_{max}}). \end{cases} $$

Note that $u(x)$ is not defined in $B_{11}$ due to the lack of control in the region. Figure A.4 illustrates the resulting vector field of (S) with the latter time-optimal control law. As a consequence of the latter theorem, the solutions of (OCP) for B-admissible trajectories are such that:

* the optimal control consists of two bang arcs $u \equiv u_{min}$ and $u \equiv u_{max}$, similar to the suboptimal control (A.4), with the switching between them occurring at the time when the trajectory reaches the separatrix $(S_{u_{max}})$;

* the optimal trajectories passes by the unstable Filippov equilibrium $(\theta_1, \theta_2)$,
<figure><img src="image_45.png" /><figcaption>Figure A.4: Stream plot of the controlled dynamics (S) with the feedback control of Theorem A.4.1. System parameters are $\gamma_1 = 1.1$, $\gamma_2 = 1.7$, $\theta_1 = 0.6$, $\theta_2 = 0.4$, $k_1 = k_2 = 1$. Control bounds are set to $u_{min} = 0.5$ and $u_{max} = 1.5$.</figcaption></figure>

which is reached by its stable manifold corresponding to dynamics of Equation (S) with $u \equiv u_{max}$. Then, the Filippov equilibrium is left by its unstable manifold corresponding to the dynamics of Equation (S) with $u \equiv u_{max}$.

The proof of this result involves showing there are no singular arcs in the optimal control, and thus $u(t)$ can only be a concatenation of bang arcs. Additionally, because of the two-dimensional affine structure in each regular domain, the sign of the switching function in the Hamiltonian can switch at most once throughout the whole interval $[0, t_f]$. Consequently, the optimal control consists of at most two bang arcs ($u_{min}$ or $u_{max}$), and the problem is reduced to finding the optimal switching time between the two arcs. An example of this trajectory and optimal control is shown in Figure A.5.

## A.5 Proof of the main results

In this section, we provide the proof for Theorem A.4.1, which is organized as follows:
<figure><img src="image_46.png" /><figcaption>Figure A.5: Optimal trajectory with $x_1^0 = 0.8$, $x_2^0 = 0.3$ and $x_2^f = 0.7$. System parameters are $\gamma_1 = 1.2$, $\gamma_2 = 1.8$, $\theta_1 = 0.6$, $\theta_2 = 0.4$, and $k_1 = k_2 = 1$. Control bounds are set to $u_{min} = 0.5$ and $u_{max} = 1.5$. Times $t_1$ and $t_2$ are the transition times at which the state meets $x_1(t_1) = \theta_1$ and $x(t_2) = (\theta_1, \theta_2)$.</figcaption></figure>

* In Section A.5.1, we reduce the problem to $B = (B_{00}, B_{01})$-admissible trajectories;

* In Section A.5.2, we prove that any optimal control admits no singular arcs in $B_{00}$;

* In Section A.5.3, we show that the optimal control in $B_{00}$ consists of two bang arcs with a switching time $t_s$ such that $x(t_s) \in (S_{u_{max}})$, and we conclude the proof of Theorem A.4.1 showing that there are no singular arcs in $B_{10}$.

### A.5.1 Reduction of the problem

Let $x(t)$ be the solution of Equation (S) such that $x(0) = x_0$ associated with an arbitrary control $u(t)$, and define the time at which the system crosses the frontier between $B_{10}$ and $B_{00}$ (respectively, between $B_{00}$ and $B_{01}$) as $t_1$ (respectively, $t_2$). We notice that the time needed to achieve a transfer between the point $(x_1(t_2), x_2(t_2)) = (x'_1, \theta_2)$ for $x'_1 < \theta_1$, and the set $\{(x_1, x_2) \in K \mid x_2 = x_2^f, 0 \le x_1 \le \theta_1\}$ does not depend on $x'_1 \le \theta_1$, and by a direct property of Equation (S) restricted to $B_{01}$, we can easily prove that the optimal control strategy for Problem (OCP) is obtained when $u(t) = u_{max}$ for $t \ge t_2$. Moreover, if $u(t) \in [u_{min}, u_{max}]$ is another optimal
control, then we obtain

$$ \int_{t_2}^{t_f} e^{-\gamma_2 s} (u(s) - u_{\max}) ds = 0, $$

hence $u(t) = u_{\max}$ for almost every $t \in [t_2, t_f]$. As a consequence, we can reduce the problem to solving (OCP) with $x_2^f = \theta_2$ among $B = (B_{10}, B_{00})$-admissible trajectories.

### A.5.2 Absence of singular arcs in $B_{00}$

In order to apply Theorem A.3.1, given the choice $B = \{B_{10}, B_{00}\}$, we set $\nu = 2$, and we define the vector fields, for $x = (x_1, x_2) \in \mathbb{R}^2$, $u \in [u_{\min}, u_{\max}]$, by

$$ f_1(x_1, x_2, u) = \begin{pmatrix} -\gamma_1 x_1 + u k_1 \\ -\gamma_2 x_2 \end{pmatrix}, \quad f_2(x_1, x_2, u) = \begin{pmatrix} -\gamma_1 x_1 + u k_1 \\ -\gamma_2 x_2 + u k_2 \end{pmatrix}. $$

The times where changes of regular domains occur for the dynamics are denoted by $t_0 = 0 < t_1$, and the final time is $t_2 = t_f$. Notice that $t_0 = 0$ is assumed to be fixed while $t_1, t_f$ are not fixed quantities a priori. We introduce the following functions $(\eta_j)_{j \in \{1, \dots, 7\}}$, which will guarantee the B-admissibility of the trajectories $x(t)$, which are solutions of Equation (S). In accordance with Definition A.3.5 of Section A.3.4, for a trajectory $x(t)$ which is solution of Equation (S), we define $p = (t_0, (t_1, x^1(t_0), x^1(t_1)), (t_f, x^2(t_1), x^2(t_f)))$. In order to guarantee the B-admissibility and the continuity of the trajectory $x(t)$ at $t = t_1$, we define the functions

$$ \left\{ \begin{aligned} \eta_1(p) &= t_0, \\ \eta_2(p) &= x_1^1(t_0) - x_1^0, \\ \eta_3(p) &= x_2^1(t_0) - x_2^0, \\ \eta_4(p) &= x_1^1(t_1) - \theta_1, \\ \eta_5(p) &= x_1^2(t_1) - \theta_1, \\ \eta_6(p) &= x_2^1(t_1) - x_2^2(t_1), \\ \eta_7(p) &= x_2^2(t_f) - \theta_2. \end{aligned} \right. $$
As in Theorem A.3.1, for $p = (t_0, (t_1, x^1(t_0), x^1(t_1))), (t_f, x^2(t_1), x^2(t_f)))$, $\alpha \in \mathbb{R}$, and $\beta = (\beta_1, \dots, \beta_7) \in \mathbb{R}^7$, define the Lagrangian

$$L(p) = \alpha t_f + \sum_{j=1}^{7} \beta_j \eta_j(p).$$

For $k \in \{1, 2\}$, the Hamiltonian $H^k$ defined in Theorem A.3.1 can be written as $H^k = H_0 + u^k H_1^k$, with $u^k \in [u_{\min}, u_{\max}]$ where, for every $x^k = (x_1^k, x_2^k) \in \mathbb{R}^2$ and $\lambda^k = (\lambda_1^k, \lambda_2^k)$,

$$ \begin{aligned} H_0(x^k, \lambda^k, \lambda_0) &= -\gamma_1 x_1^k \lambda_1^k - \gamma_2 x_2^k \lambda_2^k - \lambda_0, \\ H_1^k(x^k, \lambda^k) &= \xi_1^k k_1 \lambda_1^k + \xi_2^k k_2 \lambda_2^k, \end{aligned} $$

with $\xi_1^1 = 0$, $\xi_1^2 = 1$, $\xi_2^1 = 1$, and $\xi_2^2 = 1$. In this setting, the Adjoint State Equation (E) of Theorem A.3.1 writes

$$ \left\{ \begin{aligned} \dot{\lambda}_1^k &= \gamma_1 \lambda_1^k, \\ \dot{\lambda}_2^k &= \gamma_2 \lambda_2^k, \end{aligned} \right. \qquad (AD) $$

which is independent of $k \in \{1, 2\}$. Then, we can derive conditions from Theorem A.3.1 concerning singular arcs of Equation (OCP) along B-admissible trajectories, as defined in Section A.3.3. For $k \in \{1, 2\}$, extremal singular arcs occur when the variables $(x^k(t), \lambda^k(t), \lambda_0, u^k(t))$ are extremal and satisfy

$$H_1^k(x^k(t), \lambda^k(t)) = 0, \qquad (\text{Sing})$$

for every $t \in [T_1, T_2]$, where $t_1 \le T_1 < T_2 \le t_f$. Along such trajectories, the vanishing condition of the k-th Hamiltonian $H^k$ becomes

$$-\gamma_1 x_1^k(t) \lambda_1^k(t) - \gamma_2 x_2^k(t) \lambda_2^k(t) - \lambda_0 = 0, \qquad (V)$$

for every $t \in [T_1, T_2]$. Define the $B_{00}$ switching function as $\phi(t) = \text{sign}(k_1 \lambda_1^2(t) + k_2 \lambda_2^2(t))$ for $t \in [0, t_f]$. As a direct consequence of Theorem A.3.1, we have the following result.

**Lemma A.5.1.** At times $t_0 = 0$, $t_1$ and $t_f$, we have the following transversality and
discontinuity conditions:

$$ \left\{ \begin{array}{l} \beta_1 = -\alpha = \lambda_0, \\ \lambda_1^1(0) = \beta_2, \\ \lambda_2^1(0) = \beta_3, \\ \lambda_1^1(t_1) = -\beta_4, \\ \lambda_1^2(t_1) = \beta_5, \\ \lambda_2^1(t_1) = \lambda_2^2(t_1) = \beta_6, \\ \lambda_1^2(t_f) = 0, \\ \lambda_2^2(t_f) = -\beta_7. \end{array} \right. \qquad (TD) $$

We can deduce the following property of extremal trajectories of Problem (OCP).

**Lemma A.5.2.** *Extremal trajectories of Problem (OCP) along B-admissible trajectories admit no singular arcs in $B_{00}$, that is, for $t \in [t_1, t_2]$.*

*Proof.* In this case, Condition (Sing) becomes

$$ k_1\lambda_1^2(t) + k_2\lambda_2^2(t) = 0, $$

for $t \in [T_1, T_2]$. Differentiating this equality, we obtain

$$ k_1\gamma_1\lambda_1^2(t) + k_2\gamma_2\lambda_2^2(t) = 0, $$

for $t \in [T_1, T_2]$. Then we get, for $t \in [T_1, T_2]$

$$ \lambda_1^2(t)(\gamma_1 - \gamma_2) = 0. $$

Knowing that $\gamma_1 \neq \gamma_2$ by Assumption A.2.1, we obtain $\lambda_1^2(t) = 0$ for $t \in [T_1, T_2]$, and Condition (V) implies

$$ \lambda_1^2(t) = \lambda_2^2(t) = \lambda_0 = 0. $$

Hence, by Equation (AD), we have $\lambda_1^2(t) = \lambda_2^2(t) = \lambda_0 = 0$ for every $t \in [t_1, t_f]$. Applying Theorem A.3.1 with the functions $(\eta_j)_{j \in \{1, \dots, 7\}}$, and the Lagrangian $L$ as defined as above, we see easily that the transversality and discontinuity conditions (TD) at times $t_0 = 0$, $t_1$ and $t_f$ provide that $\alpha = \beta_j = 0$, for every $j \in \{1, \dots, 7\}$. Indeed, the condition $\lambda_1^2(t) = \lambda_2^2(t) = \lambda_0 = 0$ for every $t \in [t_1, t_f]$
implies that $\alpha = \beta_1 = \beta_5 = \beta_6 = \beta_7 = 0$. By Equation (AD), we get that
$\lambda_2^1(t) = \lambda_2^2(t) = 0$ for every $t \in [0, t_f]$, so that we can deduce $\beta_3 = \beta_6 = 0$. The
null Hamiltonian condition (V) then implies $\lambda_1^1(t) = 0$ for $t \in [0, t_1]$. It follows that
$\beta_2 = \beta_4 = 0$, so that the nontriviality condition $(\alpha, \beta) \neq 0$ of Theorem A.3.1 is
violated.
□

### A.5.3 Optimality of the two bang arcs trajectory for Problem (OCP)

Because of the two-dimensional affine structure in each regular domain, the switching function $\phi$ can switch at most once throughout the whole interval $[0, t_f]$. By reachability considerations, we can deduce the following result.

**Proposition A.5.3.** *Extremal trajectories of Problem (OCP) along B-admissible trajectories are made of two bang arcs in the domain $B_{00}$, that is, there exists $t_s \ge t_1$ such that $u(t) = u_{\min}$ for $t_1 \le t \le t_s$, and $u(t) = u_{\max}$ for $t > t_s$.*

*Proof.* By Equation (AD) and Lemma A.5.2, the switching function $\phi$ switches at most once for $t \in [t_1, t_f]$. Moreover, in order to achieve a transfer between the lines $x_1^2(t_1) = \theta_1$ and $x_2^2(t_f) = \theta_2$, at least one switch is needed. Indeed, a constant control strategy $u \equiv u_{\min}$ is such that the associated solution $x(t)$ of Equation (S) converges towards $\Phi^*(u_{\min})$ when $t \to \infty$, where $\Phi^*(u_{\min})$ is defined as in Equation (A.3), so that $x(t) < \theta_2$ for every $t \ge t_1$. Moreover, a constant control strategy $u \equiv u_{\max}$ is such that $x(t) \in B_{10}$ for $t \ge t_1$, so that $x(t)$ is not B-admissible.
□

Now we prove that that the switching time $t_s$ defined in Proposition A.5.3 for optimal trajectories of Problem (OCP) along B-admissible trajectories is such that $(x_1^2(t_s), x_2^2(t_s)) \in (S_{u_{\max}})$. First notice that, by a direct study of Equation (S) restricted to the domain $B_{00}$, we can define $t^* > 0$ as the unique time at which the solution $(y_1(t), y_2(t))$ of Equation (S) with $u \equiv u_{\min}$ and $y_1(0) = \theta_1$ and $y_2(0) = x_2^2(t_1)$ satisfies $(y_1(t^*), y_2(t^*)) \in (S_{u_{\max}})$. In order to guarantee the conditions $x_1^2(t_f) \in [0, \theta_1]$ and $x_2^2(t_f) = \theta_2$, the time $t_s$ defined in Proposition A.5.3 has to satisfy $t_s \ge t_1 + t^*$. We have the following result.

**Lemma A.5.4.** *Optimal trajectories of Problem (OCP) along B-admissible trajectories are such that $t_s = t_1 + t^*$ and $x_1^2(t_f) = \theta_1$.*
Proof. We have for $t_1 \le t \le t_s$,

$$ 
\begin{aligned}
x_1^2(t) &= (\theta_1 - \frac{k_1 u_{\min}}{\gamma_1})e^{-\gamma_1(t-t_1)} + \frac{k_1 u_{\min}}{\gamma_1}, \\
x_2^2(t) &= (x_2(t_1) - \frac{k_2 u_{\min}}{\gamma_2})e^{-\gamma_2(t-t_1)} + \frac{k_2 u_{\min}}{\gamma_2},
\end{aligned}
$$

and for $t \ge t_s$ we have

$$ 
\begin{aligned}
x_1^2(t) &= (x_1^2(t_s) - \frac{k_1 u_{\max}}{\gamma_1})e^{-\gamma_1(t-t_s)} + \frac{k_1 u_{\max}}{\gamma_1}, \\
x_2^2(t) &= (x_2^2(t_s) - \frac{k_2 u_{\max}}{\gamma_2})e^{-\gamma_2(t-t_s)} + \frac{k_2 u_{\max}}{\gamma_2}.
\end{aligned}
$$

Notice that the condition $x_2^2(t_1) < \theta_2$ implies $x_2^2(t_s) < \theta_2$ and $x_1^2(t_s) < \theta_1$. By direct computations, a time $T \ge t_s$ satisfies $x_2^2(T) = \theta_2$ if and only if

$$ T = T(t_s) \equiv t_s + \frac{1}{\gamma_2} \ln \left( \frac{-x_2^2(t_s) + k_2 u_{\max} / \gamma_2}{-\theta_2 + k_2 u_{\max} / \gamma_2} \right). $$

Notice that the condition $\theta_2 < \frac{u_{\max}}{\gamma_2}$ implies we can define a positive function $T$ : $t_s \mapsto T(t_s)$. Moreover, one can prove that, for every $t_s > 0$,

$$ T'(t_s) = \frac{k_2(u_{\max} - u_{\min})}{-\gamma_2 x_2^2(t_s) + k_2 u_{\max}}. $$

Using the fact that $(x_1^2(t_s), x_2^2(t_s))$ belongs to $B_{00}$, we obtain that $T$ is increasing on $\mathbb{R}^+$, and reaches its minimum in the interval $[t_1 + t^*, +\infty)$ at $t_s = t_1 + t^*$. The result follows from the definitions of $t^*$ and $(S_{u_{\max}})$ (see Section A.2.2). $\square$

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

There remains to understand the structure of an optimal trajectory in the regular domain $B_{10}$, that is, when $t \le t_1$. In the next proposition, we eliminate the possibility of having singular arcs in $B_{10}$ by a direct study of the dynamics of Equation (S) associated with the application of Lemma A.5.4.

**Proposition A.5.5.** *Optimal trajectories of Problem (OCP) along B-admissible trajectories admit no singular arc in $B_{10}$.*

*Proof.* Consider the solution $\bar{x}(t) = (\bar{x}_1(t), \bar{x}_2(t))$ of Equation (S) such that $u \equiv u_{\min}$ while $\bar{x}(t) \in B_{10}$, $u \equiv u_{\min}$ while $\bar{x}(t) \in B_{00} \cap (S_{u_{\max}})^-$, $u \equiv u_{\max}$ while $\bar{x}(t) \in (S_{u_{\max}})$, and the solution $\tilde{x}(t) = (\tilde{x}_1(t), \tilde{x}_2(t))$ of Equation (S) such that $u \equiv \tilde{u}(t)$ while $\tilde{x}(t) \in B_{10}$ for an arbitrary control $t \mapsto \tilde{u}(t) \in [u_{\min}, u_{\max}]$, $u \equiv u_{\min}$ while $\tilde{x}(t) \in B_{00} \cap (S_{u_{\max}})^-$, $u \equiv u_{\max}$ while $\tilde{x}(t) \in (S_{u_{\max}})$, with same initial conditions. Hence
we can define the time $\bar{T} > 0$ (respectively, $\tilde{T}$) at which we have $\bar{x}(\bar{T}) = (\theta_1, \theta_2)$ (respectively, $\tilde{x}(\tilde{T}) = (\theta_1, \theta_2)$). In order to prove that $\bar{T} \ge \bar{T}$, let us first consider the time $\tilde{t}_1 > 0$ (respectively, $\bar{t}_1 > 0$) at which $\tilde{x}(t)$ (respectively, $\bar{x}(t)$) reaches the frontier between $B_{10}$ and $B_{00}$. By a direct property of Equation (S) restricted to the domain $B_{10}$, we have $\tilde{x}_1(t) \ge \bar{x}_1(t)$ and $\tilde{x}_2(t) = \bar{x}_2(t)$ for every $t \le \min(\tilde{t}_1, \bar{t}_1)$. It follows that $\bar{t}_1 \le \tilde{t}_1$ and $\tilde{x}_2(\tilde{t}_1) \le \bar{x}_2(\bar{t}_1)$. Now consider a solution $\tilde{y}(t) = (\tilde{y}_1(t), \tilde{y}_2(t))$ (respectively, $\bar{y}(t) = (\bar{y}_1(t), \bar{y}_2(t))$) of Equation (S) with $u \equiv u_{\min}$ and such that $\tilde{y}_1(0) = \bar{y}_1(0) = \theta_1$, $\tilde{y}_2(0) = \bar{x}_2(\bar{t}_1)$ and $\bar{y}_2(0) = \bar{x}_2(\bar{t}_1)$. By a direct property of Equation (S), the times $\tilde{s}_1$ (respectively, $\bar{s}_1$) at which $\tilde{y}(t)$ (respectively, $\bar{y}(t)$) reaches $(S_{u_{\max}})$ are such that $\tilde{s}_1 \ge \bar{s}_1$. Hence we can deduce $\tilde{T} = \tilde{t}_1 + \tilde{s}_1 \ge \bar{T} = \bar{t}_1 + \bar{s}_1$, and the structure of optimal trajectories in $B_{00}$ given by Lemma A.5.4 allows to prove that an optimal trajectory for Problem (OCP) is such that $u(t) = u_{\min}$ for almost every $t \in [0, t_1]$. In particular, optimal trajectories for Problem (OCP) have no singular arcs in $B_{10}$. $\square$

The latter proposition concludes the proof of Theorem A.4.1. Additionally, as a direct consequence of the previous results, we obtain that for $t \in [0, t_f]$, the optimal control is made of a first bang arc with $u \equiv u_{\min}$ towards $(S_{u_{\max}})$ for $t \in [0, t_s]$, then a second bang arc with $u \equiv u_{\max}$ for $t \in [t_s, t_f]$ so that the system follows $(S_{u_{\max}})$ until reaching $(\theta_1, \theta_2)$.

**Remark A.5.6.** By a direct analysis of the dynamics of Equation (S) in the regular domain $B_{00}$, one can show that the time $t^* = t_s - t_1 > 0$ is the unique non-negative solution of the equation

$$ \left( \frac{k_2 u_{\min}}{\gamma_2} - x_2^2(t_1) \right) e^{-\gamma_2 t} + k_2 \frac{u_{\max} - u_{\min}}{\gamma_2} \qquad (\text{EQ}) \\ = \frac{\left( \frac{k_2 u_{\max}}{\gamma_2} - \theta_2 \right)}{\left( \frac{k_1 u_{\max}}{\gamma_1} - \theta_1 \right)^{\gamma_2/\gamma_1}} \left( \left( \frac{k_1 u_{\min}}{\gamma_1} - \theta_1 \right) e^{-\gamma_1 t} + k_1 \frac{u_{\max} - u_{\min}}{\gamma_1} \right)^{\gamma_2/\gamma_1}. $$

The latter can be obtained by solving $y_2(t^*) = \alpha(y_1(t^*), u_{\max})$, with $y_1(t^*) = x_1^2(t_1+t^*)$ and $y_2(t^*) = x_2^2(t_1+t^*)$. Equation (EQ) is hard to solve explicitly in the general case where Assumption A.2.1 is satisfied, especially because the latter assumption implies $\gamma_1 \ne \gamma_2$.
## A.6 Lower bound on the minimal time

The time required to perform a transition can be minimized to a certain extent, which is imposed by the dynamics of the system and the choice of control bounds, as shown in previous sections. In this section, we show there exist a lower bound to the minimal time. However, an explicit computation requires to solve Equation (EQ) analytically, which is a challenging task. In this section, we give a lower bound on the minimal time of Problem (OCP) in Proposition A.6.3, which is uniform w.r.t. [$u_{min}, u_{max}$] $\subset [0, +\infty)$ and is a function of the parameters $(\gamma_j)_{j \in \{1,2\}}$, $(k_j)_{j \in \{1,2\}}$, $(\theta_j)_{j \in \{1,2\}}$ satisfying Assumption A.2.1. In this purpose, we introduce an additional system which provides a lower bound for Problem (OCP). Let [$u_{min}, u_{max}$] $\subset [0, +\infty)$ be such that Assumption A.3.1 is satisfied. Then, for every $u_{min}, u_{max} \ge 0$ such that [$\bar{u}_{min}, \bar{u}_{max}$] $\subset [u_{min}, u_{max}]$, we have $\Phi^*(u_{min}) \in (S_{u_{max}})^-$. Hence, the optimal control strategy for Problem (OCP) associated with such values of $u_{min}, u_{max}$ is given by Theorem A.4.1.

**Definition A.6.1.** Define the lower trajectories as the solutions of

$$ \begin{aligned} \dot{z}_1 &= -\gamma_1 z_1 + k_1 u(t) s^-(z_2, \theta_2) \\ \dot{z}_2 &= -\gamma_1 z_2 + k_2 u(t) s^-(z_1, \theta_1), \end{aligned} \qquad (\tilde{S}) $$

with $u(t) \equiv u_{min}$ for $z(t) \in (\tilde{S}_{u_{max}})^-$, and $u(t) \equiv u_{max}$ for $z(t) \in (\tilde{S}_{u_{max}})^+ \cup \tilde{S}_{u_{max}}$, where $(\tilde{S}_{u_{max}})$ is defined as in Definition A.2.3.

A direct application of Lemma A.2.3 proves that if a lower trajectory $z(t)$ is such that $z(0) \in (\tilde{S}_{u_{max}})^-$, then $z(t)$ reaches $(\tilde{S}_{u_{max}})$ in finite time $T_{low}(u_{min}, u_{max})$. Moreover, as a consequence of the condition $\gamma_2 > \gamma_1$, we obtain the following lemma.

**Lemma A.6.2.** Consider the solution $x(t)$ of Equation (S) such that $x(0) = x_0 \in (S_{u_{max}})^-$, where $u$ is defined as in Theorem A.4.1, and the solution $z(t)$ of Equation $(\tilde{S})$ such that $z(0) = x_0$. Then we have $x_1(t) \le z_1(t)$ and $x_2(t) \le z_2(t)$, for every $t \in [0, T_{low}(u_{min}, u_{max})]$.

As a direct consequence, we get that the time $t_s$ needed by $x(t)$ in order to reach $(S_{u_{max}})$ (defined as in Proposition A.5.3) is such that $t_s \ge T_{low}(u_{min}, u_{max})$.

Hence, if we denote the minimal time for Problem (OCP) by $T_f(u_{min}, u_{max})$, then we have $T_f(u_{min}, u_{max}) \ge T_{low}(u_{min}, u_{max})$, for every $u_{min}, u_{max}$ be such that $0 \le$
<figure><img src="image_48.png" /><figcaption>Figure A.6: Different trajectories starting from $x_1 = \theta_1$ with fixed control. System parameters are $\gamma_1 = 1.4$, $\gamma_2 = 2$, $\theta_1 = 0.6$, $\theta_2 = 0.4$, and $k_1 = k_2 = 1$. Control is set to $u = u_{min}$ with $u_{min} = 0.5$. Trajectories of $(\tilde{S})$ reach its associated separatrix at the lower-bound time $T_{low}$. Vertical lines at the interception indicate $x_1(T_{low})$.</figcaption></figure>

$u_{min} \le u_{max}$. Furthermore, by definition of Problem (OCP), we have $T_f(u_{min}, u_{max}) \ge T_f(\tilde{u}_{min}, \tilde{u}_{max})$, for every $u_{min}, u_{max}, \tilde{u}_{min}, \tilde{u}_{max}$ such that $0 \le \tilde{u}_{min} \le u_{min} \le u_{max} \le \tilde{u}_{max}$. It follows that for such a choice of $u_{min}, u_{max}, \tilde{u}_{min}, \tilde{u}_{max}$, we have

$$T_f(u_{min}, u_{max}) \ge T_{low}(\tilde{u}_{min}, \tilde{u}_{max}). \qquad (A.5)$$

**Proposition A.6.3.** Set $x_0 = (x_1^0, x_2^0) \in \bar{B}_{10}$ such that $x_2^0 < \theta_2$, and consider $u_{min}, u_{max}$ such that $0 \le u_{min} \le u_{max}$. Let $x(t)$ be the solution of Equation (S) such that $x(0) = x_0$, where $u$ is defined as in Theorem A.4.1. Then we have $T_f(u_{min}, u_{max}) \ge -\frac{1}{\gamma_1} \ln(\frac{\theta_1 k_2 - \theta_2 k_1}{\theta_1 k_2 - x_2^0 k_1}) > 0$.

*Proof.* First assume that $x_1^0 = \theta_1$. Then by an adaptation of the formula given in Remark A.5.6, replacing $\gamma_2$ by $\gamma_1$, the time $T_{low}(u_{min}, u_{max})$ needed by the lower trajectory $z(t)$ to reach $(\tilde{S}_{u_{max}})$ is

$$T_{low}(u_{min}, u_{max}) = -\frac{1}{\gamma_1} \ln \left( \frac{(A(u_{max})k_1 - k_2)(u_{max} - u_{min})}{\gamma_1 \left( \left( \frac{k_2 u_{min}}{\gamma_1} - x_2^0 \right) - A(u_{max}) \left( \frac{k_1 u_{min}}{\gamma_1} - \theta_1 \right) \right)} \right),$$
where $A(u_{\max}) = \frac{k_2 u_{\max} - \theta_2}{\frac{\gamma_1}{k_1 u_{\max}} - \theta_1}$. For every $u_{\min}, u_{\max}, \tilde{u}_{\max}$ such that $0 \le u_{\min} \le u_{\max} \le \tilde{u}_{\max}$, Inequality (A.5) provides

$$T_f(u_{\min}, u_{\max}) \ge T_{\text{low}}(0, \tilde{u}_{\max}) = -\frac{1}{\gamma_1} \ln \left( \frac{(A(\tilde{u}_{\max})k_1 - k_2)\tilde{u}_{\max}}{\gamma_1(-x_2^0 + A(\tilde{u}_{\max})\theta_1)} \right).$$

Noticing that $T_{\text{low}}(0, \tilde{u}_{\max}) \to -\frac{1}{\gamma_1} \ln(\frac{\theta_1 k_2 - \theta_2 k_1}{\theta_1 k_2 - x_2^0 k_1})$ when $\tilde{u}_{\max} \to +\infty$, we deduce that $T_f(u_{\min}, u_{\max}) \ge -\frac{1}{\gamma_1} \ln(\frac{\theta_1 k_2 - \theta_2 k_1}{\theta_1 k_2 - x_2^0 k_1})$, for every $u_{\min}, u_{\max}$ such that $0 \le u_{\min} \le u_{\max}$. Moreover, Assumption A.2.1 and the condition $x_2^0 < \theta_2$ guarantee that $-\frac{1}{\gamma_1} \ln(\frac{\theta_1 k_2 - \theta_2 k_1}{\theta_1 k_2 - x_2^0 k_1}) > 0$. We deduce the general case $x_1^0 \ge \theta_1$ noticing that we have in this case $x_2(t_1) \le x_2^0$, where $t_1 \ge 0$ is the time where $x(t)$ changes regular domain from $B_{10}$ to $B_{00}$, then applying the case $x_1^0 = \theta_1$.

□

## A.7 Numerical results

We illustrate our results with numerical simulations performed with Bocop [4], an open-source toolbox for solving OCPs. In order to guarantee the reproducibility of the numerical results, the computations can be executed through an online version of Bocop<fnref n="3" />. The original problem (OCP) is solved through a direct method, by approximating it by a finite dimensional optimization problem, using a Lobato time discretization method. As the algorithm requires $s^-$ to be regularized to a smooth function, we define, for $x \in \mathbb{R}$ and $k \in \mathbb{N}$, the Hill function

$$\delta(x_i, \theta_i, k) = \frac{\theta_i^k}{x_i^k + \theta_i^k}, \qquad (A.6)$$

which can approximate $s^-$ for large values of $k$ and, when $k \to \infty$, it verifies

$$\lim_{k \to \infty} \delta(x_i, \theta_i, k) = \begin{cases} 1 & x_i < \theta_i, \\ 0 & x_i > \theta_i, \\ 1/2 & x_i = \theta_i. \end{cases}$$

<footnote n="3">³https://ct.gitlabpages.inria.fr/gallery/bistable/bistable.html</footnote>
Replacing $s^{-}$ by Hill functions (A.6) in system (S) yields the non-hybrid system

$$ egin{cases} \dot{x}_1 = -\gamma_1 x_1 + u k_1 \delta(x_2, \theta_2, k), \\ \dot{x}_2 = -\gamma_2 x_2 + u k_2 \delta(x_1, \theta_1, k). \end{cases} $$

System parameters are fixed to $\gamma_1 = 1.2$, $\gamma_2 = 2$, $\theta_1 = 0.6$, $\theta_2 = 0.4$ and $k_1 = k_2 = 1$, which verify Assumption A.2.1; and control bounds are set to $u_{min} = 0.5$ and $u_{max} = 1.5$ satisfying Assumption A.3.1. The parameter $k$ of the Hill function is set to $k = 500$, which proved an acceptable approximation of the $s^{-}$ function. Figure A.5 shows an optimal trajectory representing the transition (high, low) to (low, high). In accordance with the analytical results, the optimal control is a bang-bang control: it consists of a first phase $[0, t_s]$ of low synthesis control $u_{min}$ until $x$ reaches the separatrix ($S_{u_{max}}$), followed by a phase $[t_2, t_f]$ of high synthesis control $u_{max}$ until $x_2$ reaches $x_2^f$. As it is customary when solving OCPs with direct methods, the algorithm does not count on any a priori information of the structure of the optimal control. Yet, the obtained trajectory is in agreement with Theorem A.4.1, which confirms our theoretical results. Moreover, the solver is not restricted to consider only B-admissible trajectories, which suggests that the solution found in this work is optimal not only for Problem (OCP) along B-admissible trajectories but also for the general (OCP), without imposing the domain sequences. Figure A.7 shows different trajectories starting from $(S_{u_{max}})^{+}$ and $(S_{u_{max}})^{-}$. The streamplot represents the closed-loop dynamics for the optimal control defined in Theorem A.4.1. All trajectories starting in $(S_{u_{max}})^{-}$ approach asymptotically the point $\Phi^{*}(u_{min})$ (denoted by a cross) until they reach the separatrix, point at which the state slides over it towards the Filippov equilibrium $(\theta_1, \theta_2)$. The optimal control for trajectories starting in $(S_{u_{max}})^{+}$ consists in $u \equiv u_{max}$ for the whole interval $[0, t_f]$, and do not pass by the Filippov equilibrium.

**Remark A.7.1.** As already mentioned in the introduction, the dynamics is not uniquely defined at the undifferentiated point $(\theta_1, \theta_2)$, and the proposed solution is obtained by making a choice of dynamics at this point. Hence, concerning a biological implementation of our time-optimal strategy, it seems more reasonable to apply $u(t) \equiv u_{min}$ during a slightly longer time $\tilde{t}_s = t_s + \epsilon$ with a small $\epsilon > 0$.

In accordance with Remark A.7.1, one can be interested in comparing the sub-optimal control strategy given in Equation (A.4) with the optimal control given by
<figure><img src="image_49.png" /><figcaption>Figure A.7: Optimal trajectories starting from different initial points, with $x_2^f = 0.7$ and $k = 500$. The streamplot represents the vector field resulting from applying the optimal bang-bang strategy from Theorem A.4.1.</figcaption></figure>

Theorem A.4.1. To this purpose, one can evaluate the time loss when delaying the switch by a time $\delta t > 0$, as $t_s = t_1 + t^* + \delta t$, where $t_1$ and $t^*$ are defined as in Section A.5.3 and depend on the parameters of the system. One can show by simple computations that the difference between the times needed to reach the target $x_2^f$ for the modified trajectory w.r.t. the optimal trajectory is equal to

$$ \frac{1}{\gamma_2} \ln \left( 1 + e^{\gamma_2(\tilde{t} - t^*)} (e^{\gamma_2 \delta t} - 1) \right), $$

where

$$ \tilde{t} = \frac{1}{\gamma_2} \ln \left( \frac{k_2(u_{\min} - u_{\max})}{\gamma_2 \theta_2 - k_2 u_{\max}} \right). $$

### A.7.1 Comparison with the smooth case

In order to explore the differences between the hybrid model studied in this paper and the smooth case, we obtained optimal trajectories for the continuous dynamical model given by (A.6) with lower Hill coefficients. Figure A.8 illustrates the impact of the Hill coefficient $k$ in the functions $\delta(x, \theta, k)$, and how high values of $k$ represent
a suitable approximation of the discrete case. Figure A.9 shows optimal transitions

<figure><img src="image_50.png" /><figcaption>Figure A.8: Different Hill functions with different values of the Hill coefficient.</figcaption></figure>

for $k = 3$ compared to the streamplot obtained from the optimal bang-bang strategy. We can observe that, even for lower values of $k$, the optimal control strategy remains bang-bang, with the switches being produced after the trajectories reach a certain region not necessarily delimited by ($S_{u_{max}}$). Finally, an important difference is that the bang-bang control does not yield trajectories passing by the unstable point ($\theta_1, \theta_2$) in the continuous case. In Figure A.10, three trajectories starting from the same initial conditions are compared for different values of $k$. We observe that, as $k$ is reduced, the trajectory gets closer to the separatrix, and therefore, to the unstable point ($\theta_1, \theta_2$). Additionally, both the final time and the switching time are reduced as $k$ increases towards the idealized hybrid case.

### A.7.2 Supplementary condition $x_1(t_f) < x_1^{\max}$

In bistable systems, a binary switch implies taking the state towards the equilibria $\phi_{10}$ and $\phi_{01}$. However, as stated in Section A.3, (OCP) represents a relaxed version of this problem where $x_1(t_f) > 0$, as it is not possible to control concentration $x_1$ in $B_{01}$. In order to compare the relaxed version with the original one, we investigate
<figure><img src="image_51.png" /><figcaption>Figure A.9: Optimal trajectories starting from different initial points, with $x_2^f = 0.7$ and $k = 3$. The streamplot represents the vector field resulting from applying the optimal bang-bang strategy from Theorem A.4.1.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>Figure A.10: Optimal trajectories starting from the initial point (0.8, 0.3), with $x_2^f = 0.7$ and for different values of $k$. The streamplot represents the vector field resulting from applying the optimal bang-bang strategy from Theorem A.4.1.</figcaption></figure>

numerically the following problem:

$$
\left\{
\begin{array}{l}
\displaystyle \text{minimize } t_f \ge 0 \\
\displaystyle x(t) = (x_1(t), x_2(t)) \text{ is subject to (S)}, \\
\displaystyle x(0) = (x_1^0, x_2^0), \\
\displaystyle x_2(t_f) = x_2^f, \\
\displaystyle x_1(t_f) \in [0, x_1^{\max}), \\
\displaystyle u(\cdot) \in [u_{\min}, u_{\max}].
\end{array}
\right.
\qquad (OCP_2)
$$
In the particular case $x_2^f = k_2/\gamma_2$, solving Problem (OCP$_2$) allows to ensure that the state $x(t)$ is close enough to the steady state $\phi_{01}$ at time $t = t_f$. The main difference with (OCP) is that $x_1(t_f)$ is now constrained to the interval $[0, x_1^{\max}]$ with $x_1^{\max} < \theta_1$. For initial conditions in $B_{01}$ given by $x_1^0 \in (x_1^{\max}, \theta_1]$ and $x_2^0 = \theta_2$, we notice that the time it takes for $x_1(t)$ to reach $x_1^{\max}$ does not depend on the control $u$ (as there is no term depending on the control $u$ in the dynamics of $x_1(t)$). Therefore, the final time $t_f$ does not depend on the control, and so any control driving $x_2(t)$ from $x_2(0) = \theta_2$ to $x_2^f$ in a time $t' \le t_f$ is optimal for Problem (OCP$_2$). Thus, the problem has infinite solutions. Figure A.11 shows different trajectories for different values of $x_1^{\max}$. Among all infinite solutions, the ones found by Bocop depend on the initialization of the optimization algorithm, and have no particular meaning in the regular domain $B_{01}$. However, we verify that, as in (OCP), the switch in the control $u$ occurs at the separatrix $(S_{u_{\max}})$, and then they follow the separatrix until the point $(\theta_1, \theta_2)$. Thus, the simplest bang-bang strategy solution of (OCP$_2$) is

$$ u_1(x) = \begin{cases} u_{min} & \text{if } x \in (S_{u_{max}})^-, \\ u_{max} & \text{if } x \in (S_{u_{max}})^+ \cup (S_{u_{max}}) \text{ and } x_2 < x_2^f, \\ \frac{\gamma_2}{k_2} x_2^f & \text{if } x_2 = x_2^f. \end{cases} $$

where the control $u = x_2^f \gamma_2 / k_2$ is chosen so that $\dot{x}_2 = 0$ in the last phase. In the particular case where the final state is such that $x_2^f = k_2/\gamma_2$ (corresponding to the $x_2$-coordinate of the steady state $\phi_{01}$), the optimal control in the last phase corresponds to the open loop system $u = 1$.

## A.8 Conclusion

This paper addressed the time-optimal control problem of a bistable gene-regulatory network. Through the application of HMP, we showed that any optimal control achieving state transition is a bang-bang control, where its value is a function of the state of the system (i.e. a feedback control). While in previous works [120], the bang-bang nature of the control is imposed as a constraint, we showed that such a characteristic is necessary to produce minimum-time transitions. Results also indicate that optimal trajectories should pass by the Filippov equilibrium $(\theta_1, \theta_2)$, which represents the undifferentiated state, highly relevant from the biological point
<figure><img src="image_53.png" /><figcaption>Figure A.11: Optimal trajectories obtained with Bocop starting from the same initial point (0.8, 0.3), with $x_2^f = 0.7$ and for different values of $x_1^{\max}$. The streamplot represents the vector field resulting from applying the optimal bang-bang strategy from Theorem A.4.1. The first case (with $x_1^{\max} = \theta_1$) is the solution of (OCP).</figcaption></figure>

of view. We showed the existence of a lower bound to the minimal time, by introducing the concept of *lower trajectories*. The numerical simulations obtained through direct methods confirm our analytical results, even when no prior knowledge of the structure of the optimal trajectories is specified. The latter are obtained by approximating the piecewise behavior of the systems with Hill functions, thus simulating a non-hybrid system. Additionally, the numerical results indicate that the trajectories found are optimal not only among B-admissible trajectories, but for all solutions of the hybrid system (S). Finally, we performed a numerical comparison of the trajectories obtained for the relaxed problem (i.e. with a constraint $x_1^f \le \theta_1$) and those of the original one (i.e. with a constraint $x_1^f \le x_1^{\max} < \theta_1$), which suggests that our results are also applicable to the original problem. Our work can be related to other results in the literature. For instance, in [127], an irreversible bistable switch in *E. coli* between the genes FadR and TetR is artificially engineered by replacing an endogenous negative autoregulation loop into a positive one. The control strategy is the feed-in of fatty acid, which is chosen to be bang-bang for simplicity. Our results support such choice by proving it is not only simple, but also time optimal from a mathematical point of view. We expect that our result could be generalized to higher dimensional genetic regulatory networks, where it often occurs that trajec-
tories belonging to a given domain may bifurcate in different domains, similarly to
what happens in the toggle switch case.
Bibliography

[1] Frederick C Neidhardt and Boris Magasanik. Studies on the role of ribonucleic acid in the growth of bacteria. *Biochimica et biophysica acta*, 42:99-116, 1960.

[2] Matthew Scott, Carl W Gunderson, Eduard M Mateescu, Zhongge Zhang, and Terence Hwa. Interdependence of cell growth and gene expression: origins and consequences. *Science*, 330(6007):1099-1102, 2010.

[3] Ivan Yegorov, Francis Mairet, Hidde De Jong, and Jean-Luc Gouzé. Optimal control of bacterial growth for the maximization of metabolite production. *Journal of mathematical biology*, 78(4):985-1032, 2019.

[4] Inria Saclay Team Commands. Bocop: an open source toolbox for optimal control. http://bocop.org, 2017.

[5] Jes Forchhammer and Lasse 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(3):563-568, 1971.

[6] HDPP Bremer, Patrick 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(2):1553-69, 1996.

[7] Liujie Huo, Joachim J Hug, Chengzhang Fu, Xiaoying Bian, Youming Zhang, and Rolf Müller. Heterologous expression of bacterial natural product biosynthetic pathways. *Natural Product Reports*, 36(10):1412-1436, 2019.

[8] Hidde De Jong, Stefano Casagranda, Nils Giordano, Eugenio Cinquemani, Delphine Ropers, Johannes Geiselmann, and Jean-Luc Gouzé. Mathematical modelling of microbes: metabolism, gene expression and growth. *Journal of the Royal Society Interface*, 14(136):20170502, 2017.
[9] Itzhak Fishov, Arieh Zaritsky, and NB Grover. On microbial states of growth. *Molecular microbiology*, 15(5):789-794, 1995.

[10] Jacques Monod. The growth of bacterial cultures. *Annual review of microbiology*, 3(1):371-394, 1949.

[11] Leonor Michaelis, Maud L Menten, et al. Die kinetik der invertinwirkung. *Biochem. z.*, 49(333-369):352, 1913.

[12] Moselio Schaechter, Ole Maaløe, and Niels O Kjeldgaard. Dependency on medium and temperature of cell size and chemical composition during balanced growth of salmonella typhimurium. *Microbiology*, 19(3):592-606, 1958.

[13] Jeremy S Edwards, Rafael U Ibarra, and Bernhard O Palsson. In silico predictions of escherichia coli metabolic capabilities are consistent with experimental data. *Nature biotechnology*, 19(2):125-130, 2001.

[14] Rafael U Ibarra, Jeremy S Edwards, and Bernhard O Palsson. Escherichia coli K-12 undergoes adaptive evolution to achieve in silico predicted optimal growth. *Nature*, 420(6912):186, 2002.

[15] Erez Dekel and Uri Alon. Optimality and evolutionary tuning of the expression level of a protein. *Nature*, 436(7050):588-592, 2005.

[16] Andrea Y Weiße, Diego A Oyarzún, Vincent Danos, and Peter S Swain. Mechanical links between cellular trade-offs, gene expression, and growth. *Proceedings of the National Academy of Sciences*, 112(9):E1038-E1047, 2015.

[17] Benjamin D Towbin, Yael Korem, Anat Bren, Shany Doron, Rotem Sorek, and Uri Alon. Optimality and sub-optimality in a bacterial growth law. *Nature communications*, 8(1):1-8, 2017.

[18] Julio R Banga. Optimization in computational systems biology. *BMC systems biology*, 2(1):1-7, 2008.

[19] Kenneth J Kauffman, Purusharth Prakash, and Jeremy S Edwards. Advances in flux balance analysis. *Current opinion in biotechnology*, 14(5):491-496, 2003.

[20] Anne Goelzer, Vincent Fromion, and Gérard Scorletti. Cell design in bacteria as a convex optimization problem. *Automatica*, 47(6):1210-1218, 2011.
[21] Evert Bosdriesz, Douwe Molenaar, Bas Teusink, and Frank J Bruggeman. How fast-growing bacteria robustly tune their ribosome concentration to approximate growth-rate maximization. *The FEBS journal*, 282(10):2029–2044, 2015.

[22] Oren Shoval, Hila Sheftel, Guy Shinar, Yuval Hart, Omer Ramote, Avi Mayo, Erez Dekel, Kathryn Kavanagh, and Uri Alon. Evolutionary trade-offs, pareto optimality, and the geometry of phenotype space. *Science*, 336(6085):1157–1160, 2012.

[23] M. Scott, S. Klumpp, E.M. Mateescu, and T. Hwa. Emergence of robust growth laws from optimal regulation of ribosome synthesis. *Molecular Systems Biology*, 10:747, 2014.

[24] Manlu Zhu and Xiongfeng Dai. Growth suppression by altered (p) ppgpp levels results from non-optimal resource allocation in escherichia coli. *Nucleic acids research*, 47(9):4684–4693, 2019.

[25] Moshe Sipper. Fifty years of research on self-replication: An overview. *Artificial life*, 4(3):237–257, 1998.

[26] Douwe Molenaar, Rogier Van Berlo, Dick De Ridder, and Bas Teusink. Shifts in growth strategies reflect tradeoffs in cellular economics. *Molecular systems biology*, 5(1):323, 2009.

[27] Tracey A Lincoln and Gerald F Joyce. Self-sustained replication of an rna enzyme. *Science*, 323(5918):1229–1232, 2009.

[28] Christoph Flamm, Lukas Endler, Stefan Müller, Stefanie Widder, and Peter Schuster. A minimal and self-consistent in silico cell model based on macromolecular interactions. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 362(1486):1831–1839, 2007.

[29] Nils Giordano, Francis Mairet, Jean-Luc Gouzé, Johannes Geiselmann, and Hidde De Jong. Dynamic optimisation of resource allocation in microorganisms. In *21st International Symposium on Mathematical Theory of Networks and Systems (MTNS 2014)*, 2014.

[30] Nils Giordano, Francis Mairet, Jean-Luc Gouzé, Johannes Geiselmann, and Hidde De Jong. Dynamical allocation of cellular resources as an optimal control
problem: novel insights into microbial growth strategies. *PLoS computational biology*, 12(3):e1004802, 2016.

[31] Lev Semenovich Pontryagin. *Mathematical theory of optimal processes*. Routledge, 2018.

[32] Steffen Waldherr and Henning Lindhorst. Optimality in cellular storage via the pontryagin maximum principle. *IFAC-PapersOnLine*, 50(1):9889–9895, 2017.

[33] David W Erickson, Severin J Schink, Vadim Patsalo, James R Williamson, Ulrich Gerland, and Terence Hwa. A global resource allocation strategy governs growth transition kinetics of Escherichia coli. *Nature*, 551(7678):119–123, 2017.

[34] Yael Korem Kohanim, Dikla Levi, Ghil Jona, Benjamin D Towbin, Anat Bren, and Uri Alon. A bacterial growth law out of steady state. *Cell reports*, 23(10):2891–2900, 2018.

[35] François Bertaux, Julius Von Kügelgen, Samuel Marguerat, and Vahid Shahrezaei. A bacterial size law revealed by a coarse-grained model of cell physiology. *PLoS computational biology*, 16(9):e1008245, 2020.

[36] Tae Seok Moon, Sang-Hwal Yoon, Amanda M Lanza, Joseph D Roy-Mayhew, and Kristala L Jones Prather. Production of glucaric acid from a synthetic pathway in recombinant escherichia coli. *Applied and environmental microbiology*, 75(3):589–595, 2009.

[37] Songyuan Li, Christian Bille Jendresen, and Alex Toftgaard Nielsen. Increasing production yield of tyrosine and mevalonate through inhibition of biomass formation. *Process Biochemistry*, 51(12):1992–2000, 2016.

[38] Tat-Ming Lo, Si Hui Chng, Wei Suong Teo, Han-Saem Cho, and Matthew Wook Chang. A two-layer gene circuit for decoupling cell growth from metabolite production. *Cell systems*, 3(2):133–143, 2016.

[39] Andreas Milias-Argeitis, Marc Rullan, Stephanie K Aoki, Peter Buchmann, and Mustafa Khammash. Automated optogenetic feedback control for precise and robust regulation of gene expression and cell growth. *Nature communications*, 7(1):1–11, 2016.
[40] Diego A Oyarzún and Guy-Bart V Stan. Synthetic gene circuits for metabolic control: design trade-offs and constraints. *Journal of the Royal Society Interface*, 10(78):20120671, 2013.

[41] Irene Otero-Muras, Ahmad A Mannan, Julio R Banga, and Diego A Oyarzún. Multiobjective optimization of gene circuits for metabolic engineering. *IFAC-PapersOnLine*, 52(26):13–16, 2019.

[42] Axel Nyström, Antonis Papachristodoulou, and Andrew Angel. A dynamic model of resource allocation in response to the presence of a synthetic construct. *ACS synthetic biology*, 7(5):1201–1210, 2018.

[43] Ivan Yegorov, Francis Mairet, Hidde De Jong, and Jean-Luc Gouzé. Optimal control of bacterial growth for the maximization of metabolite production. *Journal of mathematical biology*, pages 1–48, 2018.

[44] Jérôme Izard, Cindy DC Gomez Balderas, Delphine Ropers, Stephan Lacour, Xiaohu Song, Yifan Yang, Ariel B Lindner, Johannes Geiselmann, and Hidde De Jong. A synthetic growth switch based on controlled expression of RNA polymerase. *Molecular systems biology*, 11(11):840, 2015.

[45] Arthur L Koch. Why can't a cell grow infinitely fast? *Canadian journal of microbiology*, 34(4):421–426, 1988.

[46] Hal L Smith and Paul Waltman. *The Theory of the Chemostat: Dynamics of Microbial Competition*, volume 13. Cambridge university press, 1995.

[47] Mathukumalli Vidyasagar. Decomposition techniques for large-scale systems with nonadditive interactions: Stability and stabilizability. *IEEE transactions on automatic control*, 25(4):773–779, 1980.

[48] Joseph LaSalle. Some extensions of Liapunov's second method. *IRE Transactions on circuit theory*, 7(4):520–527, 1960.

[49] Bernard Bonnard and Monique Chyba. *Singular trajectories and their role in control theory*, volume 40. Springer Science & Business Media, 2003.

[50] Yacine Chitour, Frédéric Jean, and Emmanuel Trélat. Singular trajectories of control-affine systems. *SIAM Journal on Control and Optimization*, 47(2):1078–1095, 2008.
[51] VF Borisov. Fuller's phenomenon. *Journal of Mathematical Sciences*, 100(4):2311-2354, 2000.

[52] JF Bonnans, V Grelard, and P Martinon. Bocop, the optimal control solver, open source toolbox for optimal control problems. URL http://bocop.org, 2011.

[53] Andrei A Agrachev and Yuri Sachkov. *Control Theory from the Geometric Viewpoint*, volume 87. Springer Science & Business Media, 2013.

[54] Jan Dirk Van Elsas, Alexander V Semenov, Rodrigo Costa, and Jack T Trevors. Survival of escherichia coli in the environment: fundamental and public health aspects. *The ISME journal*, 5(2):173-183, 2011.

[55] J Shu and ML Shuler. A mathematical model for the growth of a single cell of e. coli on a glucose/glutamine/ammonium medium. *Biotechnology and bioengineering*, 33(9):1117-1126, 1989.

[56] David V Goeddel, Dennis G Kleid, Francisco Bolivar, Herbert L Heyneker, Daniel G Yansura, Roberto Crea, Tadaaki Hirose, Adam Kraszewski, Keiichi Itakura, and Arthur D Riggs. Expression in escherichia coli of chemically synthesized genes for human insulin. *Proceedings of the National Academy of Sciences*, 76(1):106-110, 1979.

[57] Jérôme Harmand, Claude Lobry, Alain Rapaport, and Tewfik Sari. *The chemo-stat: Mathematical theory of microorganism cultures*. John Wiley & Sons, 2017.

[58] Meike T Wortel, Evert Bosdriesz, Bas Teusink, and Frank J Bruggeman. Evolutionary pressures on microbial metabolic strategies in the chemostat. *Scientific reports*, 6:29503, 2016.

[59] David William Spitzer. Maximization of steady-state bacterial production in a chemostat with ph and substrate control. *Biotechnology and bioengineering*, 18(2):167-178, 1976.

[60] Rixa R Lichtl, Michael J Bazin, and David O Hall. The biotechnology of hydrogen production by nostoc flagelliforme grown under chemostat conditions. *Applied microbiology and biotechnology*, 47(6):701-707, 1997.
[61] Marc C D'anjou and Andrew J Daugulis. A rational approach to improving productivity in recombinant pichia pastoris fermentation. *Biotechnology and bioengineering*, 72(1):1-11, 2001.

[62] Térence Bayen and Francis Mairet. Optimization of the separation of two species in a chemostat. *Automatica*, 50(4):1243-1248, 2014.

[63] Carlos Martínez, Olivier Bernard, and Francis Mairet. Maximizing microalgae productivity in a light-limited chemostat. *IFAC-PapersOnLine*, 51(2):735-740, 2018.

[64] Antoine Haddon, Hector Ramirez, and Alain Rapaport. First results of optimal control of average biogas production for the chemostat over an infinite horizon. *IFAC-PapersOnLine*, 51(2):725-729, 2018.

[65] Horst R Thieme. Convergence results and a poincaré-bendixson trichotomy for asymptotically autonomous differential equations. *Journal of mathematical biology*, 30(7):755-763, 1992.

[66] Sylvain Goutelle, Michel Maurin, Florent Rougier, Xavier Barbaut, Laurent Bourguignon, Michel Ducher, and Pascal Maire. The hill equation: a review of its capabilities in pharmacological modelling. *Fundamental & clinical pharmacology*, 22(6):633-648, 2008.

[67] Kenneth A Johnson and Roger S Goody. The original Michaelis constant: translation of the 1913 Michaelis–Menten paper. *Biochemistry*, 50(39):8264-8269, 2011.

[68] Andrei A Agrachev and Yuri Sachkov. *Control theory from the geometric viewpoint*, volume 87. Springer Science & Business Media, 2013.

[69] Michail I Zelikin and Vladimir F Borisov. *Theory of Chattering Control: with Applications to Astronautics, Robotics, Economics, and Engineering*. Springer Science & Business Media, 2012.

[70] AJ Van der Schaft. Symmetries in optimal control. *SIAM journal on control and optimization*, 25(2):245-259, 1987.
[71] Emmanuel Trélat and Enrique Zuazua. The turnpike property in finite-dimensional nonlinear optimal control. *Journal of Differential Equations*, 258(1):81–114, 2015.

[72] Walid Djema, Laetitia Giraldi, Sofya Maslovskaya, and Olivier Bernard. Turnpike features in optimal selection of species represented by quota models. submitted.

[73] HM Robbins. A generalized legendre-clebsch condition for the singular cases of optimal control. *IBM Journal of Research and Development*, 11(4):361–372, 1967.

[74] E Cinquemani, F Mairet, I Yegorov, H de Jong, and J-L Gouzé. Optimal control of bacterial growth for metabolite production: The role of timing and costs of control. In *2019 18th European Control Conference (ECC)*, pages 2657–2662. IEEE, 2019.

[75] Ivan Yegorov, Francis Mairet, and Jean-Luc Gouzé. Optimal feedback strategies for bacterial growth with degradation, recycling, and effect of temperature. *Optimal Control Applications and Methods*, 39(2):1084–1109, 2018.

[76] S Joe Qin and Thomas A Badgwell. A survey of industrial model predictive control technology. *Control engineering practice*, 11(7):733–764, 2003.

[77] Erdal Aydin, Dominique Bonvin, and Kai Sundmacher. Nmpc using pontryagin's minimum principle-application to a two-phase semi-batch hydroformylation reactor under uncertainty. *Computers & Chemical Engineering*, 108:47–56, 2018.

[78] Abhishek S Soni and Robert S Parker. Closed-loop control of fed-batch bioreactors: A shrinking-horizon approach. *Industrial & engineering chemistry research*, 43(13):3381–3393, 2004.

[79] Erdal Aydin, Dominique Bonvin, and Kai Sundmacher. Computationally efficient nmpc for batch and semi-batch processes using parsimonious input parameterization. *Journal of Process Control*, 66:12–22, 2018.

[80] V. Raghuraman, V. Renganathan, T. H. Summers, and J. P. Koeln. Hierarchical mpc with coordinating terminal costs*. In *2020 American Control Conference (ACC)*, pages 4126–4133, 2020.
[81] R. Scattolini and P. Colaneri. Hierarchical model predictive control. In *2007 46th IEEE Conference on Decision and Control*, pages 4803–4808, 2007.

[82] A. L. Dontchev, I. V. Kolmanovsky, M. I. Krastanov, V. M. Veliov, and P. T. Vuong. Approximating optimal finite horizon feedback by model predictive control. *Systems Control Lett.*, 139:104666, 2020.

[83] Moselio Schaechter, John L Ingraham, and Frederick Carl Neidhardt. *Microbe*. ASM Press, Washington, DC, 2006.

[84] F.C. Neidhardt. Bacterial growth: constant obsession with $dN/dt$. *J. Bacteriol.*, 181(24):7405–7408, 1999.

[85] Stanley Brul, Suzanne Van Gerwen, and Marcel Zwietering. *Modelling Microorganisms in Food*. Elsevier, 2007.

[86] Asher Brauner, Ofer Fridman, Orit Gefen, and Nathalie Q Balaban. Distinguishing between resistance, tolerance and persistence to antibiotic treatment. *Nature Reviews Microbiology*, 14(5):320–330, 2016.

[87] James C Liao, Luo Mi, Sammy Pontrelli, and Shanshan Luo. Fuelling the future: microbial engineering for the production of sustainable biofuels. *Nature Reviews Microbiology*, 14(5):288–304, 2016.

[88] H.A. van den Berg, Y.N. Kiselev, and M.V. Orlov. Optimal allocation of building blocks between nutrient uptake systems in a microbe. *Journal of Mathematical Biology*, 44(3):276–296, 2002.

[89] Hugo Dourado and Martin J Lercher. An analytical theory of balanced cellular growth. *Nature communications*, 11(1):1–14, 2020.

[90] Francis Mairet, Jean-Luc Gouzé, and Hidde de Jong. Optimal proteome allocation and the temperature dependence of microbial growth laws. *npj Systems Biology and Applications*, 7(1):1–11, 2021.

[91] M. Basan, S. Hui, H. Okano, Z. Zhang, Y. Shen, J.R. Williamson, and T. Hwa. Overflow metabolism in *Escherichia coli* results from efficient proteome allocation. *Nature*, 528(7580):99–104, 2015.
[92] A. Maitra and K.A. Dill. Bacterial growth laws reflect the evolutionary importance of energy efficiency. *Proceedings of the National Academy of Sciences of the USA*, 112(2):406-411, 2015.

[93] Dean A. Carlson, Alain B. Haurie, and Arie Leizarowitz. *Infinite Horizon Optimal Control*. Springer, Berlin, Heidelberg, 1991.

[94] J. Izard, C. Gomez Balderas, D. Ropers, S. Lacour, X. Song, Y. Yang, A.B. Lindner, J. Geiselmann, and H. de Jong. A synthetic growth switch based on controlled expression of RNA polymerase. *Molecular Systems Biology*, 11(11):840, 2015.

[95] Jan Ewald, Martin Bartl, and Christoph Kaleta. Deciphering the regulation of metabolism with dynamic optimization: an overview of recent advances. *Biochemical Society Transactions*, 45(4):1035-1043, 2017.

[96] Edda Klipp, Reinhart Heinrich, and Hermann-Georg Holzhütter. Prediction of temporal gene expression: Metabolic optimization by re-distribution of enzyme activities. *European journal of biochemistry*, 269(22):5406-5413, 2002.

[97] Diego A Oyarzún, Brian P Ingalls, Richard H Middleton, and Dimitrios Kalamatianos. Sequential activation of metabolic pathways: a dynamic optimization approach. *Bulletin of mathematical biology*, 71(8):1851, 2009.

[98] Sheng Hui, Josh M Silverman, Stephen S Chen, David W Erickson, Markus Basan, Jilong Wang, Terence Hwa, and James R Williamson. Quantitative proteomic analysis reveals a simple strategy of global resource allocation in bacteria. *Molecular systems biology*, 11(2):784, 2015.

[99] Markus Basan, Manlu Zhu, Xiongfeng Dai, Mya Warren, Daniel Sévin, Yi-Ping Wang, and Terence Hwa. Inflating bacterial cells by increased protein synthesis. *Molecular systems biology*, 11(10):836, 2015.

[100] Nathan E Lewis, Kim K Hixson, Tom M Conrad, Joshua A Lerman, Pep Charusanti, Ashoka D Polpitiya, Joshua N Adkins, Gunnar Schramm, Samuel O Purvine, Daniel Lopez-Ferrer, et al. Omic data from evolved e. coli are consistent with computed optimal growth from genome-scale models. *Molecular systems biology*, 6(1):390, 2010.
[101] David Cass. Optimum growth in an aggregative model of capital accumulation. *The Review of economic studies*, 32(3):233–240, 1965.

[102] Jean-Michel Coron, Pierre Gabriel, and Peipei Shang. Optimization of an amplification protocol for misfolded proteins by using relaxed control. *Journal of mathematical biology*, 70(1-2):289–327, 2015.

[103] Stefan Klumpp, Zhongge Zhang, and Terence Hwa. Growth rate-dependent global effects on gene expression in bacteria. *Cell*, 139(7):1366–1375, 2009.

[104] Davide Roncarati and Vincenzo Scarlato. Regulation of heat-shock genes in bacteria: from signal sensing to gene expression output. *FEMS microbiology reviews*, 41(4):549–574, 2017.

[105] Ningzi Guan and Long Liu. Microbial response to acid stress: mechanisms and applications. *Applied microbiology and biotechnology*, 104(1):51–65, 2020.

[106] Hengjiang Dong, Lars Nilsson, and Charles G Kurland. Gratuitous overexpression of genes in Escherichia coli leads to growth inhibition and ribosome destruction. *Journal of bacteriology*, 177(6):1497–1504, 1995.

[107] Nikolaos Tsiantis and Julio R Banga. Using optimal control to understand complex metabolic pathways. *BMC bioinformatics*, 21(1):1–33, 2020.

[108] Han-Yu Chuang, Matan Hofree, and Trey Ideker. A decade of systems biology. *Annual review of cell and developmental biology*, 26:721–744, 2010.

[109] Carlo Cosentino and Declan Bates. *Feedback control in systems biology*. CRC Press, 2019.

[110] Thomas Schlitt and Alvis Brazma. Current approaches to gene regulatory network modelling. *BMC Bioinformatics*, 8(6):S9, 2007.

[111] Murat Arcak and Eduardo D Sontag. A passivity-based stability criterion for a class of biochemical reaction networks. *Mathematical biosciences & engineering*, 5(1):1, 2008.

[112] Etienne Farcot and Jean-Luc Gouzé. A mathematical framework for the control of piecewise-affine models of gene networks. *Automatica J. IFAC*, 44(9):2326–2332, 2008.
[113] Leon Glass and Joel S. Pasternack. Prediction of limit cycles in mathematical models of biological oscillations. *Bull. Math. Biology*, 40(1):27-44, 1978. Papers presented at the Society for Mathematical Biology Meeting (Univ. Pennsylvania, Philadelphia, Pa., 1976).

[114] John Mallet-Paret and Hal L Smith. The poincaré-bendixson theorem for monotone cyclic feedback systems. *Journal of Dynamics and Differential Equations*, 2(4):367-421, 1990.

[115] Eduardo D. Sontag. Passivity gains and the “secant condition” for stability. *Systems & Control Letters*, 55(3):177-183, 2006.

[116] Timothy Gardner, Charles Cantor, and James Collins. Construction of a genetic toggle switch in escherichia coli. *Nature*, 403:339-42, 02 2000.

[117] Mara C Inniss and Pamela A Silver. Building synthetic memory. *Current Biology*, 23(17):R812-R816, 2013.

[118] Lucie Chambon, Ismail Belgacem, and Jean-Luc Gouzé. Qualitative control of undesired oscillations in a genetic negative feedback loop with uncertain measurements. *Automatica*, 112:108642, 2020.

[119] Nicolas Augier, Madalena Chaves, and Jean-Luc Gouzé. Qualitative control strategies for synchronization of bistable gene regulatory networks. https://hal.inria.fr/hal-02953502, 2020.

[120] Madalena Chaves and Jean-Luc Gouzé. Exact control of genetic networks in a qualitative framework: the bistable switch example. *Automatica J. IFAC*, 47(6):1105-1112, 2011.

[121] Aleksei Fedorovich Filippov. Differential equations with discontinuous right-hand side. *Matematicheskii sbornik*, 93(1):99-128, 1960.

[122] Gábor Balázsi, Alexander van Oudenaarden, and James J Collins. Cellular decision making and biological noise: from microbes to mammals. *Cell*, 144(6):910-925, 2011.

[123] Jean-Baptiste Lugagne, Sebastián Sosa Carrillo, Melanie Kirch, Agnes Köhler, Gregory Batt, and Pascal Hersen. Balancing a genetic toggle switch by real-time feedback control and periodic forcing. *Nature communications*, 8(1):1-8, 2017.
[124] Agostino Guarino, Davide Fiore, and Mario Di Bernardo. In-silico feedback control of a mimo synthetic toggle switch via pulse-width modulation. pages 680–685. 2019 18th European Control Conference (ECC), 2019.

[125] Etienne Farcot and Jean-Luc Gouzé. A mathematical framework for the control of piecewise-affine models of gene networks. *Automatica*, 44(9):2326–2332, 2008.

[126] Patrick Hillenbrand, Georg Fritz, and Ulrich Gerland. Biological signal processing with a genetic toggle switch. *PloS one*, 8(7):e68345, 2013.

[127] Ahmad A. Mannan and Declan G. Bates. Designing an irreversible metabolic switch for scalable induction of microbial chemical production. *Nature Communications*, 12(1):3419, 2021.

[128] Jr Loomis, W F and B Magasanik. Glucose-lactose diauxie in escherichia coli. *Journal of bacteriology*, 93(4):1397–1401, 04 1967.

[129] Pierre Salvy and Vassily Hatzimanikatis. Emergence of diauxie as an optimal growth strategy under resource allocation constraints in cellular metabolism. *Proceedings of the National Academy of Sciences*, 118(8), 2021.

[130] Ugo Boscain and Paolo Mason. Time minimal trajectories for a spin 1/2 particle in a magnetic field. *Journal of Mathematical Physics*, 47(6):062101, 2006.

[131] Ugo Boscain, Fredrik Grönberg, Ruixing Long, and Herschel Rabitz. Minimal time trajectories for two-level quantum systems with two bounded controls. *Journal of Mathematical Physics*, 55(6):062106, 2014.

[132] Andrei V Dmitruk and Alexander M Kaganovich. The hybrid maximum principle is a consequence of Pontryagin maximum principle. *Systems & Control Letters*, 57(11):964–970, 2008.