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

HAL
open science

# Singular regimes for the maximization of metabolite production

Agustín Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé

► To cite this version:

Agustín Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé. Singular regimes for the maximization of metabolite production. CDC 2019 - 58th IEEE Conference on Decision and Control, Dec 2019, Nice, France. pp.31-36. hal-02418444

HAL Id: hal-02418444

https://inria.hal.science/hal-02418444

Submitted on 18 Dec 2019

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

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

Agustín Gabriel Yabo<fnref n="1" />, Jean-Baptiste Caillau<fnref n="2" /> and Jean-Luc Gouzé<fnref n="1" />

*Abstract*— In nature, microorganisms are continuously facing nutrient availability changes in the environment, and thus they have evolved to dynamically adapt their physiology to cope with this phenomenon, by dynamically allocating resources to different cellular functions. In order to study their behaviours, the fitness of such microorganisms can be represented as a dynamical growth maximization strategy, which is formulated as an OCP (Optimal Control Problem) in coarse-grained self-replicator systems. This study inspired the use of biotechnological engineering to maximize the production of a metabolite of interest in *E. coli* by means of both analytical and computational techniques. Motivated by this approach, we incorporate the metabolite production scheme in a CSTR (Continuous Stirred-Tank Reactor) Bioreactor, which can be interpreted as a general case of the preceding models. We then derive two particular cases and study the associated OCP, so as to stress the importance of singular regimes and chattering arcs in optimal solutions. From a biological point of view, our results show that the natural allocation of resources of bacteria has to be modified in order to achieve optimal metabolite production. Finally, we go over the computations of the second order singular arc, and provide a numerical check of the Legendre-Clebsch condition.

## I. INTRODUCTION

In nature, microorganisms are continuously facing nutrient availability changes in the environment, and thus they have evolved to dynamically adapt their physiology to cope with this phenomenon. This is achieved through reorganization of the gene expression machinery, by dynamically allocating resources to different cellular functions. Among all possible allocation strategies, only few will guarantee survival when competing for nutrients, leading to complex and highly optimized organisms. In the specific case of *Escherichia coli*, studies have shown that, under certain conditions, bacterial populations achieve nearly maximal growth-rate, suggesting that this feature is indeed a design objective resulting from evolutionary processes [1]. These experimental results have triggered a large number of studies where the growth-rate maximization strategy is a central assumption when approaching resource allocation problems [2]. However, most of these works consider the resource allocation problem in steady-state conditions, which do not represent the natural environment of bacterial populations, thus motivating a dynamical approach to the problem. Such dynamical behaviours can be modeled through the so-called self-replicator models,

widely used in bacterial growth representations for its simplicity and its capacity to reproduce observed experimental behaviours [3].

The starting point for this line of research is [4], where the authors addressed the problem of dynamical allocation of cellular resources, showing that maximizing the steady-state growth is a sub-optimal strategy under changing environments. The dynamical growth-rate maximization can be interpreted as a biomass maximization problem during a fixed time period. Thus, it is possible to reformulate the question as an OCP (Optimal Control Problem), to be solved by means of the PMP (Pontryagin's Maximum Principle). This theoretical approach can provide gold standard allocation strategies, that can be then compared to feasible growth control implementations in bacterial cells. These results have provided a baseline understanding upon which it is possible to re-engineer the naturally-evolved behaviors of the cell in order to improve certain productivity measures. In particular, we consider the problem of producing a certain metabolite of interest, as considered in [5], for its relevance in biotechnological processes. In this regard, Optimal Control theory can not only shed light on the natural bacterial intracellular behaviours, but also help enhance industrial processes, as well as provide guidance in biotechnological research.

In this work, we present a general coarse-grained model for a self-replicating system extended with the metabolite production pathway based on [4], [5]. We start, in Section II, by considering the case of the CSTR (Continuous Stirred-Tank Reactor) Bioreactor scheme as our baseline, and we show that it is possible to derive the models previously analyzed in the field (fed-batch, constant substrate, no production of metabolites) as an initial step towards a full analysis of the new system. Then, in Section III, we focus on two dynamical problems: 1) biomass maximization when there is no production of metabolites, a feature assumed to be achieved by living organisms through evolution; and 2) product maximization under constant environmental conditions (corresponding to the fed-batch bioreactor), an artificial objective stated purely for biotechnological purposes. Ultimately, the comparison between this two problems should help to understand how to dynamically disrupt the natural allocation process in order to prioritize the metabolite production pathway instead of the population's growth rate, which is the natural behavior of bacteria. From the biological point of view, our results show that in order to optimally produce the artificial compound, the dynamical allocation of resources should be progressively altered to allocate more resources to the metabolic machinery of the cell population. Both problems are tackled through Optimal Control theory,

<footnote>*This work was partially supported by ANR project Maximic (ANR-17-CE40-0024-01), Inria IPL Cosy and Labex SIGNALIFE (ANR-11-LABX-0028-01).</footnote>

<footnote n="1">¹Université Côte d'Azur, Inria, INRA, CNRS, Sorbonne Université, Bio-core Team, Sophia Antipolis, France. agustin.yabo@inria.fr, jean-luc.gouze@inria.fr</footnote>

<footnote n="2">²Université Côte d'Azur, CNRS, Inria, LJAD, France. jean-baptiste.caillau@univ-cotedazur.fr</footnote>
and then solved using PMP. The solutions of both OCPs turn out to be singular controls, characterized by the existence of a singular arc along the solution. Consequently, we proceed in Section IV to characterize the singular trajectories in both cases, detailing the computations required to obtain the second order singular arc, and providing a numerical check of the suitable Legendre-Clebsch condition.

## II. SELF-REPLICATOR MODELS

### A. CSTR Bioreactor model with metabolite production

As previously stated, the problem of resource allocation in bacteria can be studied through the so-called self-replicator models. We consider a self-replicating system in a CSTR Bioreactor of volume $V_{ext}$. The cell is composed of the gene expression machinery (R) and the metabolic machinery (M), as seen in Figure 1a. Based on the extension introduced in [5], a metabolic pathway for the production of a certain metabolite of interest X is included (Figure 1b).

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

The system is described by three chemical macroreactions,

$$ 
\begin{aligned}
S & \xrightarrow{V_M} P, \\
P & \xrightarrow{V_R} \alpha R + (1-\alpha) M, \\
P & \xrightarrow{V_X} X.
\end{aligned}
$$

The first reaction transforms an external substrate (S) into precursor metabolites (P) and is catalyzed by M. The second one converts precursors into macromolecules, and is catalyzed by R. Finally, a third reaction transforms precursors into the product X, and is also catalyzed by M. The parameter $\alpha \in [0, 1]$ represents the resource allocation choice, and determines for each time instant the proportion of precursor allocated to the gene expression machinery, while $1 - \alpha$ indicates the proportion allocated to the metabolic machinery. The rates at which these reactions occur are $V_M$, $V_R$ and $V_X$ [g h⁻¹]. The system is subject to a constant volumetric flow rate F [L h⁻¹] which generates both an inflow of fresh medium rich in substrate, and an outflow of biomass and metabolites [6]. Moreover, the scheme is extended with the growth switch described in [7] that allows to shut off the production of ribosomes and other components of the gene expression machinery. Then, the new resource allocation variable becomes

$$ u(t) = I(t) \alpha(t), \quad u \in [0, 1], \qquad (1) $$

where I is the external control and $\alpha$ the natural allocation mechanism used in [4]. While these two control functions are supposed to act independently, we are interested in obtaining the optimal combination of them. Thus, in this work, we restrict the analysis to calculate optimal control u, without decoupling the individual controls. Then, the time evolution of the mass of each component can be written as

$$ 
\left\{
\begin{aligned}
\dot{S} &= V_{S_{in}} - V_M - V_{S_{out}}, \\
\dot{P} &= V_M - V_R - V_X - V_{P_{out}}, \\
\dot{M} &= (1-u) V_R - V_{M_{out}}, \\
\dot{R} &= u V_R - V_{R_{out}}, \\
\dot{X} &= V_X - V_{X_{out}},
\end{aligned}
\right.
\qquad (2) $$

where the inflow/outflow rates are defined as $V_{S_{out}} = DS$, $V_{P_{out}} = DP$, $V_{M_{out}} = DM$, $V_{R_{out}} = DR$, $V_{X_{out}} = DX$, $V_{S_{in}} = FS_{in}$, $s_{in}$ [g L⁻¹] being the concentration of the nutrient input, and D [h⁻¹] the dilution rate given by the relation $F/V_{ext}$. Under the assumption that the cytoplasmic density of the cells is constant throughout the population, we define the volume of the cell population $V$ [L] as

$$ V \doteq \beta(M + R), \qquad (3) $$

where $\beta$ [L g⁻¹] corresponds to the inverse of the cytoplasmic density. The definition (3) is based on the experimental fact that macromolecules explain most of the biomass in microbial cells [8]. Then, for the sake of convenience, quantities of the system are expressed as concentrations,

$$ p \doteq \frac{P}{V}, \quad r \doteq \frac{R}{V}, \quad m \doteq \frac{M}{V}, \quad s \doteq \frac{S}{V_{ext}}, \quad x \doteq \frac{X}{V} $$

where p, r and m [g L⁻¹] are intracellular concentrations of precursor metabolites, ribosomes (and other components of the gene expression machinery) and metabolic enzymes respectively (with respect to the cell population volume), s [g L⁻¹] is the extracellular concentration of substrate with respect to a constant external volume $V_{ext}$ [L]; and x [g L⁻¹] a concentration without a precise biological interpretation (since the metabolite is excreted from the cell population). As a result, it is possible to exclude the dynamics of m from the analysis since, by construction, $r + m = 1/\beta$. Replacing with concentrations leads to the following system,

$$ 
\left\{
\begin{aligned}
\dot{s} &= D(s_{in} - s) - v_M \frac{V}{V_{ext}}, \\
\dot{p} &= v_M - v_R - v_X - \mu p, \\
\dot{r} &= u v_R - \mu r, \\
\dot{m} &= (1 - u) v_R - \mu m, \\
\dot{x} &= v_X - \mu x \\
\dot{V} &= (\mu - D) V,
\end{aligned}
\right.
\qquad (4) $$

where $v_M$, $v_R$ and $v_X$ [g L⁻¹ h⁻¹] are the mass fluxes per unit volume obtained from dividing the rates $V_M$, $V_R$ and $V_X$ by $V$; and $\mu$ [h⁻¹] is the growth rate of the self-replicator
system that, using (2), is defined as

$$ \frac{\dot{V}}{V} = \frac{\dot{M} + \dot{R}}{M + R} = \beta v_R - D \quad \rightarrow \quad \mu \doteq \frac{\dot{V}}{V} \bigg|_{F=0} = \beta v_R. \qquad (5) $$

The latter basically means that the growth rate is defined as the relative variation of cell volume ($\dot{V}/V$) when there is no volumetric flow rate. The synthesis rates are modeled as Michaelis-Menten kinetics,

$$ \begin{aligned} v_M \doteq k_M m \frac{s}{K_M + s}, \quad v_R \doteq k_R r \frac{p}{K_R + p}, \\ v_X \doteq k_X m \frac{p}{K_X + p}. \end{aligned} \qquad (6) $$

with rate constants $k_M$, $k_R$, $k_X$ [h⁻¹] and half-saturation constants $K_M$, $K_R$, $K_X$ [g L⁻¹]. Moreover, quantities and time-scale are nondimensionalized in order to simplify the analysis, by defining appropriate new variables and constants,

$$ \begin{aligned} \dot{t} \doteq k_R t, \quad \dot{p} \doteq \beta p, \quad \dot{r} \doteq \beta r, \quad \dot{x} \doteq \beta x \quad \dot{X} \doteq \beta X, \\ K \doteq \beta K_R, \quad K_1 \doteq \beta K_X, \quad K_2 \doteq \beta K_M \\ k_1 \doteq \frac{k_X}{K_R}, \quad k_2 \doteq \frac{k_M}{K_R}, \end{aligned} $$

By replacing all variables and dropping all hats, system (4) becomes

$$ \left\{ \begin{aligned} \dot{s} &= D(s_{in} - s) - k_2 \frac{(1-r)s}{K_2+s} \frac{V}{V_{ext}}, \\ \dot{p} &= k_2 \frac{s(1-r)}{K_2+s} - k_1 \frac{p(1-r)}{K_1+p} - (p+1) \frac{pr}{K+p}, \\ \dot{r} &= (u-r) \frac{pr}{K+p}, \\ \dot{x} &= k_1 \frac{p(1-r)}{K_1+p} - \frac{pr}{K+p} \hat{x}, \\ \dot{V} &= \left( \frac{pr}{K+p} - D \right) V. \end{aligned} \right. \qquad (7) $$

where the dynamical expression of $m$ has been removed since $m = 1 - r$. The analysis of the model (7) can be greatly simplified by resorting to the particular case where $D = 0$, meaning that there is no substrate inflow, and no biomass output.

### B. Substrate depletion

In this case, the substrate is not replenished from the outside, a situation that can describe batch cultivation. The dynamical equation for $s$ becomes

$$ \dot{s} = -k_2 \frac{s(1-r)}{K_2 + s} \frac{V}{V_{ext}}. \qquad (8) $$

The optimal product maximization problem was already analyzed for this particular case in [5] where, due to the complexity of the computations in the PMP, most of the analysis was performed through a numerical approach. To allow an analytical study, we further simplify the system by assuming there is substrate in excess and the depletion occurs slowly enough (if, for example, $V \ll V_{ext}$). Thus, it is possible to exclude the dynamics of $s$ from the analysis.

### C. Constant environmental conditions

The environmental conditions can be modeled as constant over time as a result of $s$ being constant due some external

regulation of the variable, but it can also describe an environment with abundant substrate, where $s \gg K_M$ in (6). A constant value $E_M \doteq k_2 \frac{s}{K_2+s}$ is defined, such that model (8) becomes

$$ \left\{ \begin{aligned} \dot{p} &= E_M(1-r) - k_1 \frac{p(1-r)}{K_1+p} - (p+1) \frac{pr}{K+p}, \\ \dot{r} &= (u-r) \frac{pr}{K+p}, \\ \dot{x} &= k_1 \frac{p(1-r)}{K_1+p} - \frac{pr}{K+p}x, \\ \dot{V} &= \frac{pr}{K+p}V. \end{aligned} \right. \qquad (9) $$

This assumption can also represent fed-batch cultivation where the nutrient concentration is maintained high enough in order to achieve exponential growth rate. For this particular case study, we state the problem of product maximization as an OCP, and characterize the singular arcs of the solution.

### D. Allocation problem with no metabolite production

By overriding the production of the compound $X$, the model becomes the coarse-grained self-replicator depicted in Figure 1a. This is a particular case of the model (9), when $k_1 = 0$, and the external control $I(t)$ introduced in (1) is overridden (so that $u(t) = \alpha(t)$),

$$ \left\{ \begin{aligned} \dot{p} &= E_M(1-r) - (p+1) \frac{pr}{K+p}, \\ \dot{r} &= (\alpha - r) \frac{pr}{K+p}. \end{aligned} \right. \qquad (10) $$

The problem of biomass maximization by natural mechanisms has been extensively analyzed for this model in [4], so in this paper we merely recall its dynamics and the OCP associated, in order to compare the solution with that of the metabolite production problem.

## III. OPTIMAL CONTROL PROBLEM

### A. Dynamical biomass maximization in Model (10)

In model (10), the control input is the natural allocation $\alpha(t)$. Let $\mathcal{U}$ be the set of admissible controllers, which are Lebesgue measurable real-valued functions defined on the interval $[0, T]$ and satisfying the constraint $\alpha \in [0, 1]$. The maximization of the biomass can be formulated as an OCP,

$$ \begin{array}{ll} \text{maximize} & J_{\mu}(\alpha) = \int_{0}^{T} \mu \, dt \\ \text{subject to} & \alpha(\cdot) \in \mathcal{U} \end{array} $$

with $\mu = \frac{pr}{K+p}$ being the growth rate defined in (5). Given the state $\varphi = (p, r)$, and according to Pontryagin Maximum Principle, the Hamiltonian is defined as

$$ H_A(\varphi, \lambda, \alpha) \doteq \lambda^0 \mu + \langle \lambda, F(\varphi, \alpha) \rangle, \qquad (11) $$

where $F$ denotes the right-hand side of (10) and where $\lambda = (\lambda_p, \lambda_r)$ is the adjoint state. We assume that the process $(\alpha, \varphi, \lambda)$ satisfying PMP conditions is a normal extremal and set $\lambda^0 = -1$, which yields

$$ H_A = \lambda_p E_M (1 - r) - \frac{pr^2}{K+p} [\lambda_p (1 + p) + \lambda_r (r - \alpha) - 1]. $$
Since the Hamiltonian depends linearly on the control $\alpha$, the
optimal solution is,

$$
\alpha_{opt}(t) = \begin{cases} 0, & \text{if } \phi(\cdot) < 0, \\ 1, & \text{if } \phi(\cdot) > 0. \end{cases} \tag{12}
$$

being $\phi(\cdot) = \lambda_r \frac{pr}{K+p}$ the switching function.

B. Dynamical product maximization in Model (9)

The problem can be formulated as,

$$
\begin{array}{ll}
\text{maximize} & J_X(u) = X(T) - X_0 \\
\text{subject to} & u(\cdot) \in \mathcal{U},
\end{array}
$$

where $X_0 \doteq X(0)$ is the initial mass of metabolite. Using
the dynamical equation for $X$ in (9), the criterion can be
written as

$$
J_X(u) = \int_0^T k_1 \frac{p(1-r)}{K_1 + p} \nu \, dt.
$$

It can be proved that if $X_0 > 0$, the cost function can be
equivalently reformulated as

$$
J_X(u) = \int_0^T \frac{k_1 p(1-r)}{x K_1 + p} dt.
$$

For the sake of convenience, we do a change of variables
$y \doteq \ln x$ so that the criterion becomes

$$
J_X(u) = \int_0^T k_1 e^{-y} \frac{p(1-r)}{K_1 + p} dt,
$$

with new initial condition $y(0) = y_0 = \ln x(0)$. Given the
state $\varphi = (p, r, y)$ and adjoint state $\lambda = (\lambda_p, \lambda_r, \lambda_y)$, we
again assume that the extremal triple $(\varphi, u, \lambda)$ satisfies PMP
in normal form, and set $\lambda^0 = -1$. The Hamiltonian becomes

$$
H_B = E_M \lambda_p (1-r) - \frac{((\lambda_y+1)e^{-y} + \lambda_p)k_1 p (1-r)}{K_1 + p} - \frac{(\lambda_p(p+1) + \lambda_r(r-u) + \lambda_y)pr}{K+p}, \quad (13)
$$

and the adjoint system $\dot{\lambda} = -\frac{\partial H_B}{\partial \varphi}$ writes

$$
\left\{
\begin{aligned}
\dot{\lambda}_p &= \frac{((\lambda_y+1)e^{-y}-\lambda_p)k_1(1-r)\left(\frac{p}{K_1+p}-1\right)}{K_1+p} + \frac{\lambda_p pr}{K+p} \\
&\quad - \frac{(\lambda_p(p+1)+\lambda_r(r-u)+\lambda_y)r\left(\frac{p}{K+p}-1\right)}{K+p}, \\[1em]
\dot{\lambda}_r &= E_M \lambda_p + \frac{((\lambda_y+1)e^{-y}-\lambda_p)k_1 p}{K_1+p} \\
&\quad + \frac{(\lambda_p(p+1)+\lambda_r(2r-u)+\lambda_y)p}{K+p}, \\[1em]
\dot{\lambda}_y &= -\frac{k_1(\lambda_y+1)p(r-1)e^{-y}}{K_1+p}.
\end{aligned}
\right.
$$

As in the previous model, the Hamiltonian for this OCP still
depends linearly on $u$ and so the solution is again as (12),
with the same switching function $\phi(\cdot) = \lambda_r \frac{pr}{K+p}$.

C. Characterization of singular regimes

Singular trajectories play a major role in optimal control
theory [9], [10], and their characterization is a necessary
step towards a complete description of the optimal solution.

In both explored OCPs, and as expected in linear optimiza-
tion problems, maximization of the Hamiltonian gives no
information when the switching function vanishes during a
whole interval of time. In this case, the OCP has a singular
arc, that can be obtained by computing the successive time
derivatives of the switching function until it is possible to
obtain an explicit expression of the control. In the case of the
simple resource allocation problem solved in III-A, vanishing
of the switching function implies than $\lambda_r = 0$, so $\dot{\lambda}_r = 0$
along the singular arc. Since the Hamiltonian is preserved
along an extremal trajectory, $\lambda_p$ is constant along a singular
arc, which means that such arc corresponds to a steady-state
of the system. Moreover, it has been shown in [4] that the
singular arc is in fact the optimal steady-state<fnref n="1" />, and so the
optimal control is $\alpha(t) = \alpha_{opt}$. Additionally, the arc is of
order two, which implies that it has to be entered and exited
through a chattering arc [11], i.e. an arc with an infinite
number of switchings. We now focus on the model (9) and
provide a detailed computation of the corresponding singular
arcs. The Hamiltonian (13) can be expressed as

$$
H_B = H_0 + u H_1
$$

where,

$$
\begin{align*}
H_0 &= E_m \lambda_p (1-r) + \frac{((\lambda_y + 1)e^{-y} - \lambda_p)k_1 p (1-r)}{K_1 + p} \\
    &\quad - \frac{(\lambda_p (p+1) + \lambda_r r + \lambda_y)pr}{K + p}, \\
H_1 &= \phi(t) = \frac{\lambda_r pr}{K + p}.
\end{align*}
$$

Assume $H_1$ vanishes on a whole sub-interval $\mathcal{I} = [t_1, t_2] \subset [0, T]$. The switching surface is the set (here $n = 3$)

$$
\Sigma \doteq \{(\varphi, \lambda) \in \mathbb{R}^{2n} \mid H_1 = 0\}.
$$

The time derivative of $H_1$ along the extremal is equal to the
Poisson bracket $\{H_0, H_1\}$, denoted as $H_{01}$, since

$$
\begin{align*}
\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_B}{\partial \lambda_i} \frac{\partial H_1}{\partial \varphi_i} - \frac{\partial H_B}{\partial \varphi_i} \frac{\partial H_1}{\partial \lambda_i} \right) \\
&= \{H_B, H_1\} = \{H_0 + uH_1, H_1\} = \{H_0, H_1\} = H_{01}.
\end{align*}
$$

Using the latter, we proceed to compute the subsequent
derivatives of $H_1$ with respect to time, which should also
vanish along the singular arc, until it is possible to compute
the singular control $u$:

$$
\begin{equation}
\begin{split}
0 = H_{01} = -\frac{pr^2}{(K+p)^2} \lambda_r p + \frac{\lambda_r p^2 r^2}{(K+p)^2} & \\
-\left(E_m(r-1) - \frac{k_1 p(r-1)}{K_1+p} + \frac{(p+1)pr}{K+p}\right) \left(1 - \frac{pr}{K+p}\right) & \frac{\lambda_r r}{K+p} \\
+\left(E_m \lambda_p + \frac{((\lambda_y+1)e^{-y}-\lambda_p)k_1 p}{K_1+p} + \frac{(\lambda_p(p+1)+\lambda_r r+\lambda_y)p}{K+p}\right) & \frac{pr}{K+p},
\end{split}
\tag{14}
\end{equation}
$$

<footnote n="1">¹By optimal steady-state we mean the steady-state with constant input $\alpha_{opt}^*$ that maximizes the integrand of the cost function.</footnote>
and evaluating (14) for $z \in \Sigma$ (which, in this case, is just setting $\lambda_r = 0$) yields

$$0 = H_{01}(z) = \frac{\left(E_m \lambda_p + \frac{((\lambda_y+1)e^{-y}-\lambda_p)k_1 p}{K_1+p} + \frac{(\lambda_p(p+1)+\lambda_y)p}{K+p}\right)pr}{K+p}. \quad (15)$$

From (15) we have the following new condition:

$$E_m \lambda_p + \frac{((\lambda_y + 1)e^{-y} - \lambda_p)k_1 p}{K_1 + p} + \frac{(\lambda_p(p + 1) + \lambda_y)p}{K + p} = 0,$$

which defines a subset of the switching surface $\Sigma$,

$$\Sigma' \doteq \{(\varphi, \lambda) \in \mathbb{R}^{2n} | H_{01} = 0\} \cap \Sigma.$$

Similarly, the second derivative can be computed as

$$\begin{aligned} \ddot{H}_1 = \dot{H}_{01} &= \{H_B, H_{01}\} \\ &= \{H_0, H_{01}\} + u\{H_1, H_{01}\} = H_{001} + uH_{101}, \end{aligned}$$

and so one has

$$\begin{aligned} H_{101} = &-pr\phi_7 \frac{\lambda_r r}{(K+p)^2} - \frac{\phi_3 r \left(\frac{p}{K+p}-1\right)}{(K+p)^2} \lambda_r p \\ &+ \frac{\left(E_m \lambda_p + \frac{k_1 p \phi_2}{K_1+p} + \frac{p \phi_1}{K+p}\right)}{(K+p)^2} p^2 r + \lambda_r p \frac{p^2 r^2}{(K+p)^3} \\ &+ pr\phi_7 \frac{\lambda_r pr}{(K+p)^3} + \frac{\lambda_r \phi_3 \left(\frac{p}{K+p}-1\right)}{(K+p)^2} pr \\ &+ \frac{\lambda_r \phi_7 \left(\frac{p}{K+p}-1\right)}{(K+p)^2} pr^2 \end{aligned}$$

where

$$\begin{aligned} \phi_1 &= \lambda_p(p+1) + 2\lambda_r r + \lambda_y \\ \phi_2 &= (\lambda_y+1)e^{-y} - \lambda_p \\ \phi_3 &= E_m(r-1) - \frac{k_1 p(r-1)}{K_1+p} + \frac{(p+1)pr}{K+p} \\ \phi_7 &= \frac{\phi_3 + E_m - \frac{k_1 p}{K_1+p}}{r}. \end{aligned}$$

Evaluating in $z \in \Sigma'$ yields

$$\begin{aligned} H_{101}(z) &= \frac{\left(E_m \lambda_p + \frac{((\lambda_y+1)e^{-y}-\lambda_p)k_1 p}{K_1+p} + \frac{(\lambda_p(p+1)+\lambda_y)p}{K+p}\right)p^2 r}{(K+p)^2} \\ &= H_{01}(z) \frac{p}{K+p} = 0 \end{aligned}$$

on $\Sigma'$, entailing that singular arcs must be at least of order two. Since $\ddot{H}_1$ should vanish along the singular arc, we calculate $H_{001}$ and evaluate it in $z \in \Sigma'$:

$$0 = H_{001}(z) = \frac{k_1 p \phi_2 \phi_7 \left(\frac{p}{K_1+p}-1\right)}{(K+p)(K_1+p)} + \frac{k_1 (\lambda_y+1) p^3 e^{-y}}{(K+p)^2 (K_1+p)} - \frac{\left(\frac{k_1 \phi_2 \left(\frac{p}{K_1+p}-1\right)}{K_1+p} - \frac{\lambda_p p}{K+p} + \frac{\phi_1 \left(\frac{p}{K_1+p}-1\right)}{K+p}\right) \left(E_m - \frac{k_1 p}{K_1+p}\right) p}{K+p},$$

which defines the set

$$\Sigma'' \doteq \{(\varphi, \lambda) \in \mathbb{R}^{2n} | H_{001} = 0\} \cap \Sigma'.'$$

As previously stated, the procedure should be repeated until $u$ appears explicitly. For this problem, going up to the fourth derivative of $H_1$, one is able to retrieve the control using

Poisson brackets of length five,

$$0 = H_{0001} + uH_{10001},$$

provided $H_{10001}$ is not zero. This turns out to be the case, and one has to check the Kelley (or generalized Legendre-Clebsch) condition (here for $k=2$)

$$(-1)^k \frac{\partial}{\partial u} \left[ \frac{d^{2k}}{dt^{2k}} \left( \frac{\partial H_B}{\partial u} \right) \right] < 0, \quad t \in \mathcal{I},$$

along the singular arc. This condition is necessary for optimality and we devise a numerical check in next section. As for problem III-A, having an order two singular control implies that singular arcs can only be entered and exited through chattering (Fuller phenomenon), as illustrated by the numerical computations presented next.

## IV. NUMERICAL SOLUTION

The numerical verifications were performed with Bocop [12], using the following parameters: $\beta = 0.003$ L g$^{-1}$, $E_M = 1$, $k_R = 1.6$ h$^{-1}$, $K_R = 1$ h$^{-1}$, $k_X = 1$ h$^{-1}$, $K_X = 1$ h$^{-1}$ and $V_{ext} = 8$ L. The number of time steps was fixed to 5000, and the tolerance to $1e-14$. The discretization method used was the sixth-order Lobatto III C$^2$. In both cases, initial conditions were set to $p(0) = 0.024$, $r(0) = 0.18$ and $V(0) = 1$; and in the product maximization problem $y(0) = 0$. As expected, both optimal control solutions are

<figure><img src="image_3.png" /><figcaption>(a) Biomass maximization (III-A) (b) Product maximization (III-B)</figcaption></figure>

<caption>Fig. 2: Solutions of the OCPs computed using Bocop.</caption>

characterized by a singular arc of order two that is entered and exited through chattering (Figure 2). Moreover, the solution for the biomass maximization problem $\alpha_{opt}$ matches the optimal steady-state input $\alpha_{opt}^*$ along the circular arc (Figure 2a), as predicted in the computations. However, it is not the case for the product maximization problem, where the singular solution moves away from the optimal steady-state over time (Figure 2b) showing that a constant input is a sub-optimal control strategy. This result suggests that, in order to maximize the production of the metabolite $X$, it is necessary to induce the microbial population to increasingly allocate more resources to the metabolic machinery as time passes by. Thus, the external control $I$ should act on the natural allocation $\alpha(t)$ along the singular arc to match the

<footnote>$^2$Bocop definition files are available from the authors on request.</footnote>
lower value $u(t)$ depicted in Figure 2b. In both cases, the optimal steady-state inputs $\alpha_{opt}^*$ and $u_{opt}^*$ were computed numerically by selecting the steady-state input that maximizes the integrand of each criterion. The fact that, for the product maximization problem, the singular arc is of order two, can be verified by evaluating the derivatives of $H_1$ over the optimal trajectories (Figure 3). In this figure, all subsequent derivatives vanish along the sub-interval $\mathcal{I}$, except for the fourth one. We eventually provide a numerical check of the generalized Legendre-Clebsch condition.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Simulation of the OCP III-B (Product maximization) with *Bocop*. The intervals where the functions vanish are highlighted in light red. All functions vanish along the singular arc but $\frac{\partial}{\partial u}\left[\frac{d^4}{dt^4}\left(\frac{\partial}{\partial u}H_B\right)\right]$, highlighted in green, which is negative as expected.</figcaption></figure>

## V. CONCLUSIONS

Self-replicator models of bacterial growth are capable of reproducing the growth laws of certain organisms at steady-state. However, in [4], it has been shown that these models can also account for dynamical environments. In this context, natural biomass maximization was achieved through a bang-singular strategy. In [5], authors showed that the same kind of strategy is necessary in order to maximize the synthetic production of a metabolite. In line with this work, we have proposed a metabolite production scheme in a CSTR Bioreactor, with the particularity that it encompasses previously studied models in the field. We have focused on the most relevant particular cases to emphasize the importance of singular regimes and chattering arcs on optimal control solutions.

We have provided an analytical derivation of the results, as well as a numerical characterization of the singular regimes. Contrary to the natural biomass maximization process, maximizing the production of a metabolite is accomplished by a singular solution that does not contain any steady-state, showing that a time-varying action is indeed required to achieve maximization. For this particular problem, an optimal external control should increasingly allocate more resources to the metabolic machinery, while reducing the allocation to the gene expression machinery. Results raise interesting questions from the optimal control point of view as it can provide ideal control scenarios, but also about how to implement these open-loop strategies in real environments. Indeed, in order to do so, further analysis is required for the newly presented CSTR Bioreactor model, which is our current objective.

## REFERENCES

[1] R. U. Ibarra, J. S. Edwards, and B. O. Palsson, "Escherichia coli k-12 undergoes adaptive evolution to achieve in silico predicted optimal growth," *Nature*, vol. 420, no. 6912, p. 186, 2002.

[2] E. Bosdriesz, D. Molenaar, B. Teusink, and F. J. Bruggeman, "How fast-growing bacteria robustly tune their ribosome concentration to approximate growth-rate maximization," *The FEBS journal*, vol. 282, no. 10, pp. 2029-2044, 2015.

[3] A. L. Koch, "Why can't a cell grow infinitely fast?," *Canadian journal of microbiology*, vol. 34, no. 4, pp. 421-426, 1988.

[4] N. Giordano, F. Mairet, J.-L. Gouzé, J. Geiselmann, and H. De Jong, "Dynamical allocation of cellular resources as an optimal control problem: novel insights into microbial growth strategies," *PLoS computational biology*, vol. 12, no. 3, p. e1004802, 2016.

[5] I. Yegorov, F. Mairet, H. De Jong, and J.-L. Gouzé, "Optimal control of bacterial growth for the maximization of metabolite production," *Journal of mathematical biology*, pp. 1-48, 2018.

[6] H. L. Smith and P. Waltman, *The theory of the chemostat: dynamics of microbial competition*, vol. 13. Cambridge university press, 1995.

[7] J. Izard, C. D. G. 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*, vol. 11, no. 11, p. 840, 2015.

[8] H. Bremer, P. P. Dennis, et al., "Modulation of chemical composition and other parameters of the cell by growth rate," *Escherichia coli and Salmonella: cellular and molecular biology*, vol. 2, no. 2, pp. 1553-169, 1996.

[9] B. Bonnard and M. Chyba, *Singular trajectories and their role in control theory*, vol. 40. Springer Science & Business Media, 2003.

[10] Y. Chitour, F. Jean, and E. Trélat, "Singular trajectories of control-affine systems," *SIAM Journal on Control and Optimization*, vol. 47, no. 2, pp. 1078-1095, 2008.

[11] V. Borisov, "Fuller's phenomenon," *Journal of Mathematical Sciences*, vol. 100, no. 4, pp. 2311-2354, 2000.

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