# Optimal bacterial resource allocation

Agustin Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé

► To cite this version:

Agustin Gabriel Yabo, Jean-Baptiste Caillau, Jean-Luc Gouzé. Optimal bacterial resource allocation. DCO 2021 - International Conference on Dynamic Control and Optimization, Feb 2021, Aveiro, Portugal. hal-03528000

HAL Id: hal-03528000

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

Submitted on 17 Jan 2022

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.
# Optimal bacterial resource allocation

Jean-Baptiste Caillau

Université Côte d'Azur, CNRS, Inria, LJAD

Joint work with

Jean-Luc Gouzé

Inria

Agustín Yabo

Université Côte d'Azur, Inria

## Dynamic Control and Optimization

Aveiro, February 2021
Honouring Andrey V. Sarychev's 65

Encounter first with the Portuguese team

<figure><img src="image_1.png" /><figcaption>Geometric control summer school, Bedlewo (2003)</figcaption></figure>

"Even if he falls asleep during a talk, he will ask THE good question!"

Then Trieste, Dijon, Moscow, Paris... (often with Bernard)
Optimal control of microbial cells (ANR Maximic, PI H. de Jong)

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

► Self-replicator model: biomass maximisation (Giordano et al'16)

► Genetic engineering of a strain of bacteria: light-induced control of growth

► New pathway to produce a metabolite of interest (Egorov et al'18)

► Competition between biomass maximisation and metabolite production

Yabo, A.; Caillau, J.-B.; Gouzé, J.-L. Optimal bacterial resource allocation: metabolite production in continuous bioreactors. *Math. Biosci. Eng.* **17** (2020), no. 6, 7074–7100
# Model definition

## Self-replicator model

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

The model represents three chemical macroreactions between substrate (S), precursors (P), genetic machinery/ribosomes (R), metabolic machinery (M) and produced metabolite (X):

$$S \xrightarrow{V_M} P$$

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

$$P \xrightarrow{V_X} X$$

New control: $u(t) := \alpha(t)I(t)$
## Continuous stirred-tank reactor dynamics

In terms of concentrations s, p, r, x (and using $m+r=1$) in the CSTR,

$$ S_1: \left\{ \begin{array}{l} \dot{s} = D(s_{in} - s) - v_M(s, 1-r)\gamma, \\ \dot{p} = v_M(s, 1-r) - v_X(p, 1-r) - \mu(p, r)(p+1), \\ \dot{r} = (u-r)\mu(p, r), \\ \dot{v} = v_X(p, 1-r)\gamma - Dx, \\ \dot{\gamma} = (\mu(p, r) - D)\gamma, \end{array} \right. $$

where $\gamma$ is the volume of the cell population, $D$ is the dilution rate, $s_{in}$ is the substrate inflow, and

$$ \mu(p, r) := \beta v_R(p, r) $$

($\beta$ = inverse of cytoplasmic density) is the growth rate of the bacterial population.

**Assumption.** The synthesis rates $v_M$, $v_R$ and $v_X$ are monotone increasing kinetics (e.g. Michaelis-Menten, $v_R(p, r) := k_R rp/(K_R + p)$) and

$$ v_X(p, 1-r) = c(r) v_R(p, r). $$
# Asymptotic behaviour

## Mass conservation

**Proposition.** The set

$$ \Gamma = \{(s, p, r, x, \mathcal{V}) \in \mathbf{R}^5 : s_{in} \ge s \ge 0, p \ge 0, x \ge 0, 1 \ge r \ge 0, \mathcal{V} \ge 0\} $$

is positively invariant for the initial value problem.

We assume $\mathcal{V}(0) > 0$ and $r(0) > 0$ (no growth otherwise).

**Proposition.** The $\omega$-limit set of any solution of system $S_1$ with constant input $u(t) = \bar{u}$ lies in the set

$$ \Omega := \{(s, p, r, x, \mathcal{V}) \in \mathbf{R}^5 : s + \left(p + \frac{r}{\bar{u}}\right) \mathcal{V} + x = s_{in}\}. $$

This result is used to define a limit system and reduce the dimension by two.
# Asymptotic behaviour

## Limiting system

As $r(t) \to \bar{u}$ when $t \to \infty$ and since $x$ can be eliminated ($\Omega$), we define

$$ S_1' : \begin{cases} \dot{s} = D(s_{in} - s) - \bar{v}_M(s)\nu, \\ \dot{p} = \bar{v}_M(s) - \bar{\mu}(p)(p + \bar{c} + 1), \\ \dot{\nu} = (\bar{\mu}(p) - D)\nu, \end{cases} $$

where

$$ \begin{aligned} \bar{v}_M(s) &:= v_M(s, 1 - \bar{u}), \quad \bar{v}_R(p) := v_R(p, \bar{u}), \quad \bar{v}_X(p) := v_X(p, 1 - \bar{u}), \\ \bar{\mu}(p) &:= \mu(p, \bar{u}), \quad \bar{c} := c(\bar{u}). \end{aligned} $$

We study the local and global stability properties of this reduced system.
# Asymptotic behaviour

## Local stability

**Theorem.** Define

$$ \bar{f}(p) := \bar{v}_R(p) + \bar{v}_X(p) + \bar{\mu}(p)p = \bar{\mu}(p)(p + \bar{c} + 1), $$

and let $p_w$ be the solution of $\bar{f}(p) = \bar{v}_M(s_{in})$. Then,

* If $\bar{\mu}(p_w) \ge D$:

* The interior equilibrium $E_i$ exists, is unique and locally stable.

* The washout equilibrium $E_w$ exists, is unique and locally unstable.

* If $\bar{\mu}(p_w) < D$:

* The interior equilibrium $E_i$ does not exist.

* The washout equilibrium $E_w$ exists, is unique and locally stable.

* The interior equilibrium $E_i := (s_i, p_i, \gamma_i)$ with

$$ \bar{\mu}(p_i) = D, \quad \bar{v}_M(s_i) = \bar{f}(p_i), \quad \gamma_i := \frac{D(s_{in} - s_i)}{\bar{v}_M(s_i)}. $$

* The washout equilibrium is $E_w := (s_{in}, p_w, 0)$.
Asymptotic behaviour

Sketch of proof. If $\bar{\mu}(p_w) \ge D$,

$$
J_i = \begin{bmatrix} -D - \bar{v}'_M(s) \psi_i & 0 & -\bar{v}_M(s) \\ \bar{v}'_M(s) & \bar{f}'(p) & 0 \\ 0 & \bar{\mu}'(p) \psi_i & 0 \end{bmatrix},
$$

and the characteristic polynomial is

$$
P_i(\lambda) = (\lambda + D)(\lambda^2 + \lambda \underbrace{(D + \bar{v}'_M(s)\psi_i + (p+\bar{c}+1)\bar{\mu}'(p))}_{>0} + \bar{v}'_M(s)\psi_i \underbrace{(p+\bar{c}+1)\bar{\mu}'(p)}_{>0})
$$

which, by Routh-Hurwitz criterion, implies that all eigenvalues have negative
real part. If $\bar{\mu}(p_w) < D$,

$$
J_w = \begin{bmatrix} -D & 0 & -\bar{c} \\ \bar{v}'_M(s) & -\bar{f}'(p) & 0 \\ 0 & 0 & \bar{\mu}(p_w) - D \end{bmatrix}
$$

with characteristic polynomial

$$
P_w(\lambda) = (\lambda + D)(\lambda + \bar{f}'(p))(\lambda - \bar{\mu}(p_w) + D). \quad \square
$$
# Asymptotic behaviour

## Global analysis

**Proposition.** The $\omega$-limit set of any solution of the limiting system $S_1'$ lies in the set

$$ \Omega' := \{(s, p, \mathcal{Y}) \in \mathbb{R}^3 : s + (p + \bar{c} + 1) \mathcal{Y} = s_{in}\} $$

As $s = s_{in} - \mathcal{Y}(p + \bar{c} + 1)$, one can further reduce $S_1'$ to

$$ S_1'' := \begin{cases} \dot{p} = \bar{v}_M(s(\cdot)) - \bar{\mu}(p)(p + \bar{c} + 1), \\ \dot{\mathcal{Y}} = (\bar{\mu}(p) - D)\mathcal{Y}. \end{cases} $$

**Proposition.** Every solution of the limiting system $S_1''$ converges to

* $E_i^* := (p_i, \mathcal{Y}_i)$ if $\bar{\mu}(p_w) \ge D$,

* $E_w^* := (p_w, 0)$ if $\bar{\mu}(p_w) < D$.

**Sketch of proof.** Bendixson-Dulac criterion allows to discard periodic orbits and cycles of equilibria as

$$ \begin{aligned} & \frac{\partial}{\partial p} \dot{p} + \frac{\partial}{\partial \mathcal{Y}} \dot{\mathcal{Y}} = \frac{\partial}{\partial p} \bar{v}_M(s(\cdot)) - \bar{\mu}(p) - \bar{\mu}'(p)(p + \bar{c} + 1) + \bar{\mu}(p) - D \\ & = \bar{v}'_M(s(\cdot)) \frac{\partial s(\cdot)}{\partial p} - \bar{\mu}'(p)(p + \bar{c} + 1) - D < 0. \quad \square \end{aligned} $$
Asymptotic behaviour

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

Results on asymptotically autonomous systems (Thieme'92) eventually allow to relate the asymptotic behaviour of the 2-dimensional limiting system $S_1''$ to the behaviour of the full 5-dimensional system $S_1$.

**Theorem.** Every solution of system $S_1$ converges to

► the extended interior equilibrium $\hat{E}_i := (s_i, p_i, \bar{u}, x_i, \mathcal{Y}_i)$, $x_i := \bar{c}\mathcal{Y}_i$, if $\bar{u}(p_w) \ge D$,

► the extended washout equilibrium $\hat{E}_w := (s_{in}, p_w, \bar{u}, 0, 0)$ otherwise.
# Dynamic optimisation

## Metabolite production

Maximisation of the metabolite X is considered over a finite time horizon:

$$ \text{maximize} \quad D \mathcal{V}_{\text{ext}} \int_{0}^{T} x(t) dt, $$

resorting to $(q := (s, p, r, x, \mathcal{V}) \in \mathbb{R}^5)$

$$ \text{(OCP):} \quad \left\{ \begin{array}{ll} \text{maximize} & \displaystyle\int_{0}^{T} x(t) dt \\ \text{subject to} & \dot{q}(t) = F(q(t), u(t)), \\ & q(0) = q_0, \quad q(T) \text{ free,} \\ & u(t) \in [0, 1]. \end{array} \right. $$

**Proposition.** Existence holds.
# Dynamic optimisation

## Extremal flow

We write the (normal) Hamiltonian in terms of state $q := (s, p, r, x, \mathcal{V})$ and costate $\lambda := (\lambda_s, \lambda_p, \lambda_r, \lambda_x, \lambda_\mathcal{V})$:

$$H(q, \lambda, u) = -\lambda^0 x + \langle \lambda, F(q, u) \rangle,$$

where $F$ is the right-hand side of $S_1$. As the dynamics is control-affine, $H = H_0 + uH_1$ with

$$ 
\begin{aligned}
H_0 &= \lambda_s (D(s_{in} - s) - v_M(s, 1-r)\mathcal{V}) + \lambda_p (v_M(s, 1-r) - v_X(p, 1-r) - \mu(p, r)(p+1)) \\
&\quad - r\lambda_r \mu(p, r) + \lambda_x (v_X(p, 1-r)\mathcal{V} - Dx) + \lambda_\mathcal{V} (\mu(p, r) - D)\mathcal{V} - \lambda^0 x, \\
H_1 &= \lambda_r \mu(p, r).
\end{aligned}
$$

By virtue of Pontrjagin's maximum principle, a.e. one has

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

and the control is singular if $H_1$ vanishes along a whole subarc.
# Dynamic optimisation

## Singular arcs

Along a singular arc

$$0 = \dot{H}_1 = \frac{\partial H_1}{\partial q} \dot{q} + \frac{\partial H_1}{\partial \lambda} \dot{\lambda} = \sum_{i=1}^{n} \left( \frac{\partial H}{\partial \lambda_i} \frac{\partial H_1}{\partial q_i} - \frac{\partial H}{\partial q_i} \frac{\partial H_1}{\partial \lambda_i} \right) = \{H, H_1\} = \{H_0, H_1\} := H_{01}$$
and

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

**Proposition.** Singular arcs are at least of order two.

$$H_{101} = \mu(p, r) \left[ \dot{\lambda}_r \mu_r(p, r) + \lambda_r \frac{\partial}{\partial r} \left( \mu_r(p, r) \dot{r} + \mu_p(p, r) \dot{p} \right) \right] - \lambda_r \mu_r(p, r) \frac{\partial H_{01}}{\partial \lambda_r} - \lambda_r \mu_p(p, r) \frac{\partial H_{01}}{\partial \lambda_p}$$

which is also equal to 0 when $H_1 = 0$. $\square$

In the case of an order two singular control ($H_{10001} < 0$ along the arc, generalized Legendre-Clebsch condition)

$$u = - \frac{H_{00001}}{H_{10001}}.$$

Although the order is local ($H_{101}$ is not everywhere zero), one expects connection between bang and singular arcs through Fuller phenomenon.
Dynamic optimisation

Fuller phenomenon

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

Generalized Legendre-Clebsch condition

<figure><img src="image_6.png" /></figure>
## Turnpike property

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

A Lagrange pair $(\bar{q}, \bar{u}, \bar{\lambda})$ of $(SP)_X$ defines a hyperbolic equilibrium of the Hamiltonian singular flow on $\Sigma'' := \{(q, \lambda) | H_1 = H_{01} = H_{001} = H_{00001} = 0\}$,

$$H_s(q, \lambda) := H(q, \lambda, u_s(q, \lambda)), \quad u_s(q, \lambda) := -\frac{H_{00001}}{H_{10001}}(q, \lambda).$$

See Trélat-Zuazua'14, Pighin'16 in the regular case, ongoing work with S. Maslovskaya, W. Djema, J.-B. Pomet in the singular case
# Static optimisation

## Metabolite production vs. biomass maximisation

Maximisation is now performed wrt. the static control $\bar{u}$ plus the dilution rate $D$ of the CSTR:

$$ (SP_Y) : \begin{cases} \text{maximize} \\ \text{subject to} \end{cases} \quad \begin{cases} D\bar{Y} \\ F(\bar{q}, \bar{u}) = 0, \\ 0 \le \bar{u} \le 1, \end{cases} $$

$$ (SP_X) : \begin{cases} \text{maximize} \\ \text{subject to} \end{cases} \quad \begin{cases} D\bar{x} \\ F(\bar{q}, \bar{u}) = 0, \\ 0 \le \bar{u} \le 1. \end{cases} $$

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

**Proposition.** Any solution of the two static problems lies in the region of existence of the interior equilibrium, $\bar{\mu}(p_w) > D$.
# Static optimisation

## Discussion

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

* 20% decrease of biomass (Y) but 500% increase of produced metabolite (X)

* Difference between static and "gold standard" dynamic optimisation negligible for large fixed final time (< 1%)

* Richer models: ubiquity of Fuller and (singular) turnpike phenomena

* Second order conditions in the Fuller case: preliminary approach in Agrachev-Beschastnyi'2020

* Reproducible research: ct (= control toolbox) initiative ct.gitlabpages.inria.fr/gallery
Standing on the shoulders of giants

<figure><img src="image_10.png" /><figcaption>Revaz Valerianovich Gamkrelidze
(born February 4, 1927)</figcaption></figure>

Bon
Anniversaire
Andrey !

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

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

<figure><img src="image_13.png" /><figcaption>Andrey Sarychev</figcaption></figure>