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

HAL
open science

# Optimal resource allocation in bacterial growth: theoretical study and applications to metabolite production

Agustín Gabriel Yabo

► To cite this version:

Agustín Gabriel Yabo. Optimal resource allocation in bacterial growth: theoretical study and applications to metabolite production. Optimization and Control [math.OC]. Université Côte d'Azur, 2021. English. NNT: 2021COAZ4095. tel-03636842v2

HAL Id: tel-03636842

https://theses.hal.science/tel-03636842v2

Submitted on 11 Apr 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.
<figure><img src="image_2.png" /></figure>

# THÈSE DE DOCTORAT

Allocation optimale de ressources dans la croissance bactérienne : étude théorique et applications à la production de métabolites

## Agustín Gabriel Yabo

Institut National de Recherche en Sciences et Technologies du Numérique (Inria)
Équipe Projet Biocore / McTao

**Présentée en vue de l'obtention du grade de docteur en Automatique, Traitement du Signal et des Images d'Université Côte d'Azur**

**Dirigée par:** Jean-Luc Gouzé

**Co-dirigée par:** Jean-Baptiste Caillau

**Soutenue le:** 9 decembre 2021

**Devant le jury, composé de:**

Julio R. Banga, Professeur, Instituto de Investigaciones

Marinas, CSIC

Jean-Baptiste Caillau, Professeur, Université Côte d'Azur

Madalena Chaves, Directrice de recherche, Inria Sophia

Antipolis-Méditerranée

Jean-Luc Gouzé, Directeur de recherche, Inria Sophia

Antipolis-Méditerranée

Hidde de Jong, Directeur de recherche, Inria Grenoble -

Rhône-Alpes

Alain Rapaport, Directeur de recherche, INRAE Montpellier

Mohab Safey El Din, Professeur, Sorbonne Université
# Allocation optimale de ressources dans la croissance bactérienne : étude théorique et applications à la production de métabolites

## Optimal resource allocation in bacterial growth: theoretical study and applications to metabolite production

### Jury:

#### Président du jury

Madalena Chaves, Directrice de recherche, Inria Sophia Antipolis-Méditerranée (France)

#### Rapporteurs

Julio R. Banga, Professeur, Instituto de Investigaciones Marinas, CSIC (Espagne)

Alain Rapaport, Directeur de recherche, INRAE Montpellier (France)

#### Examinateurs

Jean-Baptiste Caillau, Professeur, Université Côte d'Azur (France)

Jean-Luc Gouzé, Directeur de recherche, Inria Sophia Antipolis-Méditerranée (France)

Hidde de Jong, Directeur de recherche, Inria Grenoble - Rhône-Alpes (France)

Mohab Safey El Din, Professeur, Sorbonne Université (France)
# Allocation optimale de ressources en croissance bactérienne : étude théorique et application à la production de métabolite

Agustín Gabriel Yabo

## Résumé

Les micro-organismes évoluent sous la pression de la sélection naturelle, améliorant leur capacité à proliférer dans leur environnement en développant des réseaux métaboliques optimisés. Des études ont montré que les populations bactériennes peuvent atteindre un taux de croissance presque maximal dans certaines conditions, ce qui leur permet de supplanter les espèces concurrentes. Considérer l'auto-réplication microbienne comme un problème d'allocation de ressources est une nouvelle approche qui a répondu avec succès à certaines des questions sous-jacentes dans le domaine. Dans cette approche, les ressources cellulaires disponibles sont affectées dynamiquement à différentes fonctions comme le métabolisme ou la synthèse des protéines. Ce cadre a également motivé de nombreuses applications à la production artificielle de métabolites d'intérêt ; l'objectif principal est alors de détourner les ressources cellulaires des voies natives vers une voie hétérologue dans le but de synthétiser efficacement un composé spécifique (par exemple des agents antitumoraux, des antibiotiques, de l'insuline, des agents immunosuppresseurs, etc.) À cette fin, des techniques biotechnologiques récentes permettent de contrôler de manière externe la croissance bactérienne en interrompant l'expression de l'ARN polymérase. Cette thèse se concentre sur les aspects mathématiques d'une certaine classe de modèles d'auto-réplicateurs basés sur les principes d'allocation des ressources susmentionnés. Ces modèles, basés sur des hypothèses minimales, sont étonnamment efficaces pour rendre compte des lois de croissance empiriques bien étudiées des cultures microbiennes. Tout au long du manuscrit, nous revisitons certains des cadres industriels les plus pertinents pour la croissance bactérienne parmi lesquels la culture par lots et les bioréacteurs continus, ainsi que d'autres modèles simplifiés où la concentration en nutriments reste constante. L'idée est de comparer les stratégies d'allocation des ressources évoluant naturellement, où l'objectif est de maximiser la biomasse de la population bactérienne, avec les stratégies artificielles visant à maximiser la production d'un métabolite voulu. L'étude implique une analyse dynamique et une optimisation des modèles proposés, et nous avons recours à la théorie du contrôle optimal pour déterminer les stratégies d'allocation conformes à ces objectifs, tant d'un point de vue analytique que numérique. Ces stratégies optimales constituent des références de choix pour le développement de stratégies de rétroaction basées sur la mesure en temps réel des processus industriels.

Mots clés : biologie des systèmes, croissance bactérienne, contrôle optimal, biotechnologie, systèmes dynamiques non linéaires.
# Optimal resource allocation in bacterial growth: theoretical study and applications to metabolite production

Agustín Gabriel Yabo

## Abstract

Microorganisms have evolved under the pressure of natural selection, improving their capacity to proliferate in the environment by developing highly optimized metabolic networks. Studies showed that bacterial populations can achieve nearly maximal growth rate under certain conditions, allowing them to outgrow competing species. Considering microbial self replication as a resource allocation problem is a novel approach that successfully answered some of the underlying questions in the field. This perspective considers the problem of dynamically assigning the available cellular resources to different cellular functions, such as metabolism and protein synthesis. The framework has also motivated numerous application to the artificial production of metabolites of interest, where the main objective is to divert the cellular resources from the native pathways into a heterologous pathway in order to efficiently synthesize a specific compound (e.g. antitumour agents, antibiotics, insulin, immunosuppressive agents, etc.). To this end, recent biotechnological techniques allow to externally control bacterial growth by shutting off the expression of RNA polymerase. This thesis focuses on the mathematical aspects of a certain class of self-replicator models based on the aforementioned resource allocation principles. These models, based on minimal assumptions, are surprisingly effective in accounting for well-studied empirical growth laws of microbial cultures. Throughout the manuscript, we revisit some of the most relevant industrial frameworks for bacterial growth, such as batch cultivation and continuous bioreactors, and other simplifying approximations where the nutrient concentration remains constant. The objective is to compare the naturally-evolved resource allocation strategies, where the objective is to maximize the biomass of the bacterial population, with the artificial strategies aiming to maximize the production of the metabolite of interest. The study involves dynamical analysis and optimization of the proposed models, and we resort to Optimal Control theory to find the allocation strategies complying with these objectives, both from analytical and numerical points of view. Ultimately, such optimal strategies may provide guidance for developing feedback strategies based on real-time measuring of industrial processes.

Keywords: systems biology, bacterial growth, optimal control, biotechnology, nonlinear dynamical systems.
Acknowledgements

I would like to start by thanking both my supervisors Jean-Luc Gouzé and Jean-Baptiste Caillau, without whom this manuscript would not have been written. I am specially grateful for their trust, the freedom I was given to make my own path and to explore my own ideas, and for pushing me to learn "la langue de Molière".

I am also very grateful to the members of the jury in the PhD defense, in particular to Alain Rapaport and Julio R. Banga for the effort of thoroughly going through this PhD manuscript, and the rich feedbacks they gave me. I also appreciate the participation of Hidde de Jong and Madalena Chaves, who have been part of the project in many different levels (as members of both the Personal Monitoring Committee and the PhD jury), and with whom we shared very fruitful scientific discussions. Additionally, I am grateful to Hidde and the whole Microcosme team from Inria Grenoble - Rhône-Alpes (specially to Johannes Geiselmann, Eugenio Cinquemani and Antrea Pavlou) for welcoming me and helping me during the (brief but) exciting experimental phase of my PhD. Finally, I thank Mohab Safey El Din for taking the time to participate in the jury and for the ongoing joint work.

I am more than grateful to the teams Biocore and McTao for the excellent research environment we shared these years, and for the friends I have done during that period, with whom we shared lots of excellent moments around a coffee/mate/beer. I particularly enjoyed the scientific works we have done with my friend Nicolas Augier. I also acknowledge the inspiring time I spent at the Ctrl-A and Datamove teams in Grenoble, during my master's internship.

I would like to thank the Université Côte d'Azur for funding the project, as well as its Inria research center and the ANR project Maximic (ANR-17-CE40-0024-01).

Agradezco fuertemente el apoyo que me brindó Céline, qui m'a accompagné pendant la dernière année de ma thèse, et qui m'a offert la meilleure journée de soutenance de thèse qu'un humble doctorant puisse souhaiter.

No puedo dejar de agradecer a quien me viera dar mis primeros pasos en el mundo de la investigación, Damián E. Oliva, que supo orientarme cuando más lo necesité; y a los miembros del Departamento de Ciencia y Tecnología de la Universidad Nacional de Quilmes.

Finalemente, les debo un infinito agradecimiento a mis padres Andrea y Gabriel, y a mi hermana Anita, quienes siempre fueron mi modelo a seguir. Gracias por brindarme tantas herramientas y motivarme siempre a hacer lo que me apasiona.
# Contents

## 1 Introduction 6

1.1 Mathematical modelling 7

1.2 Optimality in nature 9

1.3 Self-replicator models 11

1.4 A dynamical perspective 12

1.5 Metabolite production 13

1.6 Studied problems 14

## 2 Constant substrate inflow 18

2.1 Introduction 18

2.2 Self-replicator models 20

2.3 Dynamical analysis of model (S) 25

2.4 Optimal control problem 33

2.5 Numerical solution 42

2.6 Conclusions 43

## 3 Substrate depletion 46

3.1 Introduction 46

3.2 Model definition 47

3.3 Model analysis 51

3.4 The biomass maximization case 55

3.5 The product maximization case 61

3.6 Discussion 66

## 4 Continuous bioreactor 68

4.1 Introduction 68

4.2 Model definition 71
<table><tr><td>4.3 Metabolite production</td><td>86</td></tr><tr><td>4.4 Discussion</td><td>98</td></tr></table>

5 Model-predictive control schemes 101

<table><tr><td>5.1 Introduction</td><td>101</td></tr><tr><td>5.2 Model definition</td><td>103</td></tr><tr><td>5.3 Artificial metabolite production</td><td>108</td></tr><tr><td>5.4 Numerical results</td><td>113</td></tr><tr><td>5.5 Conclusion</td><td>113</td></tr></table>

6 A generalized resource allocation model of microbial growth 116

<table><tr><td>6.1 Introduction</td><td>116</td></tr><tr><td>6.2 Model definition</td><td>119</td></tr><tr><td>6.3 Asymptotic behavior</td><td>125</td></tr><tr><td>6.4 Model calibration</td><td>132</td></tr><tr><td>6.5 Optimal resource allocation</td><td>134</td></tr><tr><td>6.6 Biologically relevant scenarios</td><td>145</td></tr><tr><td>6.7 Conclusion</td><td>149</td></tr></table>

7 Conclusion 152

A Time-optimal control of piecewise affine bistable gene-regulatory networks 156

<table><tr><td>A.1 Introduction</td><td>156</td></tr><tr><td>A.2 Bistable-switch model</td><td>160</td></tr><tr><td>A.3 Time-optimal transfer</td><td>165</td></tr><tr><td>A.4 Main results</td><td>171</td></tr><tr><td>A.5 Proof of the main results</td><td>172</td></tr><tr><td>A.6 Lower bound on the minimal time</td><td>180</td></tr><tr><td>A.7 Numerical results</td><td>182</td></tr><tr><td>A.8 Conclusion</td><td>187</td></tr></table>
Thesis-related publications

[AYT1] Yabo, A., J.B. Caillau, J.L. Gouzé, Hidde de Jong, and Francis Mairet. Dynamical analysis and optimization of a generalized resource allocation model of microbial growth (accepted for publication). *SIAM Journal on Applied Dynamical Systems*, 2021.

[AYT2] Agustín Gabriel Yabo, Jean-Baptiste Caillau, and Jean-Luc Gouzé. Optimal bacterial resource allocation: metabolite production in continuous bioreactors. *Mathematical Biosciences and Engineering*, 17(6):7074-7100, 2020.

[AYT3] Yabo, A., J.B. Caillau, and J.L. Gouzé. Hierarchical MPC applied to bacterial resource allocation and metabolite synthesis (accepted for publication). In *The 60th IEEE conference on Decision and Control*, 2021.

[AYT4] Agustín Gabriel Yabo and Jean-Luc Gouzé. Optimizing bacterial resource allocation: metabolite production in continuous bioreactors. *IFAC-PapersOnLine*, 53(2):16753-16758, 2020.

[AYT5] Agustín Gabriel Yabo, Jean-Baptiste Caillau, and Jean-Luc Gouzé. Singular regimes for the maximization of metabolite production. In *2019 IEEE 58th Conference on Decision and Control (CDC)*, pages 31-36. IEEE, 2019.

[AYT6] Agustín Yabo, Jean-Baptiste Caillau, and Jean-Luc Gouzé. Bacterial growth strategies as optimal control problems: maximizing metabolite production. In *FGS'2019-19th French-German-Swiss conference on Optimization*, 2019.
Other publications

[AYO1] N. Augier and Yabo, A. Time-optimal control of piecewise affine bistable gene-regulatory networks (accepted for publication). *International Journal of Robust and Nonlinear Control*, 2021.

[AYO2] N. Augier and Yabo, A. Time-optimal control of piecewise affine bistable gene-regulatory networks: preliminary results. In *7th IFAC Conference on Analysis and Design of Hybrid Systems*, 2021.

[AYO3] C. Djuikem, Yabo, A., F. Grognard, and S Touzeau. Mathematical modelling and optimal control of the seasonal coffee leaf rust propagation. In *7th IFAC Conference on Analysis and Design of Hybrid Systems*, 2021.
Online contributions

All the works presented in this manuscript contain a numerical section includ-
ing simulations of the dynamical systems and trajectories of the optimal control
problems. In order to guarantee the usability and reproducibility of the results, a
series of online notebooks have been developed, which allow readers to execute
the simulations and optimal control algorithms. The examples are part of the
gallery (https://ct.gitlabpages.inria.fr/gallery/) of the ct (control toolbox)
project, which contains several examples of optimal control problems solved numer-
ically using a python wrapper of bocop3 (https://www.bocop.org/) that can be
executed on a Binder live web environment (https://mybinder.org/). The note-
books are listed below:

* Optimal resource allocation in a generalized model of microbial growth
https://ct.gitlabpages.inria.fr/gallery/bacteria/bacteria.html

* Bacterial biomass maximization with substrate depletion
https://ct.gitlabpages.inria.fr/gallery/substrate/depletion.html

* Bistable gene-regulatory networks
https://ct.gitlabpages.inria.fr/gallery/bistable/bistable.html

* Seasonal coffee leaf rust propagation
https://ct.gitlabpages.inria.fr/gallery/clr/clr.html
# Chapter 1

## Introduction

Microbiological research, while recent in comparison with the study of other living creatures, has been the main task of many scientists across the globe, even when the available technological means were scarce. There is a vast number of compelling reasons that explain this fact. For instance, understanding bacteria is a first step towards understanding—and potentially controlling—the diseases caused by them. In addition, most of the antibiotics manufactured by the pharmaceutical industry are produced by bacteria as a *secondary metabolite*, so named for not being involved in the growth and reproduction of these organisms. Numerous compounds essential in medical practices are also synthesized by bacteria, such as antitumor agents, immunosuppressive agents and insulin. The latter is also the case in the food industry, for products such as yogurt, cheese and pickles [7].

Whereas investigating such microorganisms usually starts with an experimental phase in a laboratory, a fundamental part of the discipline involves the exploitation of the obtained data, and the challenge of drawing knowledge from it. Due to its complexity, this task brings together scientists from very different backgrounds, such as mathematicians, biologists, statisticians and computer scientists. While the theoretical and practical approaches used in each discipline are quite diverse, scientists agree on a very common, interdisciplinary tool: mathematical models, widely used in the domain—and in science in general—for their capacity to capture and predict key behaviours of these living systems.
## 1.1 Mathematical modelling

A powerful approach used to study bacterial growth is mathematical modelling, for its proved success in representing complex networks of biochemical reactions behind internal cellular processes such as metabolism and gene expression [8]. Due to the complexity of even the smallest living organisms, simple features such as the growth rate of a bacterial population is usually governed by a large number of internal reactions forming these complex networks, and thus it is not expected to be explained by simple deterministic laws. However, through continuous experimental analysis of growing microbial culture in controlled environments, it is possible to observe certain recurring phenomena, and to describe them through empirical mathematical expressions. For instance, during the exponential growth phase of a bacterial culture, the rate of growth is constant, and thus it is reasonable—and customary—to consider all relative concentrations in the cell as constants [9]. At this stage, the number of cells in the population can be modelled using a simple dynamical equation

$$ \dot{N} = \mu N, \qquad (1.1) $$

where $N$ is the number of cells, and $\mu$ the constant population growth rate. Thus, the amount of cells at time $t$ is given by

$$ N(t) = N_0 e^{\mu t}. $$

In the same line, a major milestone in microbiology is the work of Monod [10], for successfully being able to describe the relation between microbial growth and its medium. The main contribution was the concept of limiting nutrient: an empirical relationship between the exponential steady-state growth rate $\mu$ of a growing culture of *Escherichia coli* and the concentration $s$ of available substrate in the medium. The latter can be represented by an hyperbolic function

$$ \mu(s) = \mu_{\max} \frac{s}{K_s + s}, $$
where $\mu_{max}$ is the maximal possible bacterial growth rate, and $K_s$ a constant. Then, equation (1.1) can be written as

$$ \dot{N} = \mu_{max} \frac{s}{K_s + s} N $$

The latter can be related to classical Michaelis-Menten dynamics [11], an equation proposed more than a hundred years ago predicting a proportionality relation between the rate of an enzyme-catalyzed reaction and the concentration of the substrate consumed in the reaction.

Almost ten years after Monod's work, experimental studies found that the concentration of RNA (ribonucleic acid) in microbial cells was directly linked to the growth rate [12]. More precisely, the results showed an empirical linear relation between the RNA/total protein ratio and the growth rate at steady state in bacteria [1], as illustrated in Figure 1.1. This study was one of the first relating growth rate and the cell macromolecular composition, and motivated further investigation on these empirical laws, currently known in the literature as *growth laws*. These

<figure><img src="image_3.png" /><figcaption>Figure 1.1: The relation between steady-state growth rate and RNA content of *A. aerogenes*, strain 5-P14 [1].</figcaption></figure>

results have been instrumental in the way we understand microorganisms, and are
still valid nowadays [2], as seen in Figure 1.2.

<figure><img src="image_4.png" /><figcaption>Figure 1.2: Correlation of the RNA/protein ratio *r* with growth rate *λ* for various strains of E. coli. The growth rate is modulated by changing the quality of nutrients [2].</figcaption></figure>

An essential machinery of bacteria is the metabolism, defined as the ensemble of internal chemical reactions responsible for the production of energy from nutrients, growth and replication. There has been numerous works focusing on the prediction of metabolic networks during exponential growth (i.e. steady-state conditions) aiming to better understand the internal resources distribution [13, 14]. These studies aim to understand and predict how bacteria adapt to different mediums, such as acetate, glucose, and glycerol. The obtained results support an assumption that has been adopted in many previous works: under certain conditions, bacteria adapt the way they distribute their internal resources so as to maximize their growth rate.

## 1.2 Optimality in nature

A common hypothesis when studying bacterial growth is that bacteria seek to maximize their growth rate [15, 2, 16, 17]. From an evolutionary perspective, the latter can be viewed as a result of the pressure of natural selection: in environments with scarce nutrient, a higher growth rate represents an evolutionary advantage allowing
species to outgrow competitors. This concept is known as fitness, and is usually represented in mathematical models as optimality with respect to a specific criteria (i.e. growth rate, growth yield, ...) [18]. Different modelling techniques have been employed to predict the optimality principles behind bacterial growth. A classical example is FBA (Flux Balance Analysis), an approach that aims to predict metabolic flux distributions in a quantitative way based on certain optimization criteria, while complying with a number of physical and chemical constraints such as mass balance and energy balance [19]. As a result, this technique is able to predict (sometimes multiple) steady-state solutions of internal fluxes inside a cell based solely on the stoichiometry of the metabolic networks. However, in certain cases, this technique can require imposing very strong simplifying assumptions in order to predict maximal growth rate, as the approach tends to exclude some important biomass constituents from the model [8]. Other similar approaches such as RBA (Resource Balance Analysis) consider the optimization problem together with the production of proteins and enzyme-related cost functions [20].

More recent studies focus on the problem of adapting to changing environments from a resource allocation perspective. The approach is based on the idea that the mechanisms behind the allocation of internal resources to different cellular functions have been optimized through evolution to maximize the population growth rate [21]. The advantage of this framework is the capacity to yield models encompassing mechanistic trade-offs related to physical and biochemical limitations. An example can be seen in [16], where the model is constructed taking into account constraints such as the limited size of the pools of energy and ribosomes in a cell. Intuitively, producing proteins to tackle a specific task implies diverting resources from other tasks, which creates numerous inherent trade-offs between cellular functions. There are well-studied examples in the literature, such as the trade-off between nutrient uptake and protein synthesis [2], or growth and stress response [22]. Through this methodology, it has been possible to link the optimal regulation of ribosome synthesis in the cell to a well-known natural regulation mechanism based on the presence of guanosine tetraphosphate (also known as the ppGpp alarmone) [23]. Recent experimental results were able to confirm such theoretical predictions by showing how the artificial increase or reduction of the concentration of the ppGpp molecule in *E. coli* leads to non-optimal resource allocation strategies [24].

A natural way of representing bacterial growth is through self-replicator models,
which are becoming increasingly popular for their simplicity and their potential to predict the empirical growth laws above mentioned. In this context, the impact of simple natural trade-offs in bacterial growth can be studied through the question: what is the resource allocation strategy that maximizes the rate of replication of a self-replicating system?.

## 1.3 Self-replicator models

The study of self-replicator systems is a domain that has drawn the attention of scientists for more than 70 years. The questions of whether a system is able to self replicate, and how, have been of interest even before the discovery of the DNA and the current available knowledge in biological research. Its interest does not only come from living beings, but also from its applications to space exploration, nanotechnology, and computer science in general [25]. However, cellular replication has been a major driver in the study of self replication, for its potential to understand the biological mechanisms behind it, to transfer this knowledge to artificial systems, and to optimize the existing biological processes.

Bacterial growth is not strictly a replication process, as each new individual in the culture is subject to a strong genetic variability, producing mutations from one generation to another. However, representing microbial growth as a process of self replication has proven surprisingly helpful in elucidating multiple natural mechanisms of unicellular organisms [26]. A very simple example of a real self-replicating biological system consists of an RNA enzyme that synthesises itself from a precursor [27]. As briefly explained in previous sections, if the substrate in the medium is constant, the population is expected to undergo exponential growth. Self-replicator models have also been useful in studying fundamental biological mechanisms as homeostasis and circadian rhythms [28]. Other examples were able to reproduce a metabolic strategy adopted by unicellular organisms known as overflow metabolism, which consists of an energetic inefficiency in bacteria growing at a high growth rate [26].

Most of these models consider the process of self replication in a steady-state setting, which is rarely a condition attained by bacteria in nature. While these works can be useful to study phenomena that occur in reduced timescales, or in controlled environments such as a laboratory, some questions regarding the adaptation of living
systems to changing environments require an explicit formulation of such dynamism,
which entail a different mathematical approach to the problem.

1.4 A dynamical perspective

The changing nature of the environments in which bacteria inhabit motivates a dy-
namical perspective of the resource allocation problem [29]. A recent paper has
focused in the dynamical allocation of cellular resources facing nutrient upshifts
[30], where the bacterial cell is represented through a coarse-grained self-replicator
dynamical model described by a system of ODEs. The approach considers a sim-
plified proteome (i.e. the entire set of proteins in a cell) divided into metabolic
enzymes (responsible for nutrient uptake) and ribosomal proteins (responsible for
the synthesis of macromolecules) and tackles the question of how to allocate re-
sources between these two cellular functions so as to optimally transition between
two different environments. The main hypothesis behind the study is that the re-
source allocation strategy has been optimized by natural selection to maximize the
bacterial growth rate. In a system of ODEs, the mathematical problem of finding
the time-varying function that maximizes a certain criterion is called an OCP (Op-
timal Control Problem), and it can be solved through PMP (Pontrjagin's Maximum
Principle) [31]. Results show that the optimal allocation strategy consists of a tran-
sient period followed by a steady-state phase matching previous static results [26].
The transient is composed of oscillations between the two limit allocations: maximal
enzyme synthesis and maximal ribosomal production. Surprisingly, the natural reg-
ulatory mechanism behind the ppGpp molecule produces a similar control structure
to that of the optimal allocation found through the optimal control approach.

The approach of Giordano et al. is a prime example of how simple natural
trade-offs can be studied through abstract mathematical models. A similar work
investigates the trade-off between storage and enzymatic compounds in cells [32].
By hypothesis, the natural objective is to cumulate biomass, which is described by
both compounds, but the metabolic capacity of the bacterial cell is supposed to be
affected only by enzymes. The paper also tackles the dynamical optimization prob-
lem using a PMP approach. A quantitative model of dynamical resource allocation
describing transitions between different nutrients has been proposed in [33], which is
independent of the kinetic parameters behind each reaction, thus providing robust-
ness of the results to parametric uncertainty. A resource allocation model was used to describe bacterial growth during nutritional upshifts [34], suggesting that, in unstable environments, ribosomal sub-saturation can represent a selective advantage for allowing faster growth rate increases after nutrient upshifts. A more recent work proposed a minimalistic whole-cell coarse-grained dynamical model representing the proteome composition, growth rate and size of the cell, that takes into account a class of proteins dedicated to cell division [35]. While the model considers that the growth rate is affected by environmental conditions such as changes in the nutrient quality or over-expression of growth rate-independent proteins, it is not based on the maximal growth rate principle, in contrast with most of the studies in the literature.

Dynamical approaches based on minimal models have been key in understanding the regulation mechanisms behind natural trade-offs, but it have also proved fruitful in analyzing trade-offs between natural and artificially engineered resource pathways as, for instance, those targeting the synthesis of value-added compounds.

## 1.5 Metabolite production

A classic strategy for producing value-added metabolites and other chemical compounds is to artificially engineer the production of an heterologous enzyme into a cell host, which will thus synthesize the targeted metabolite. For instance, *E. coli* can be modified to produce glucuronic and glucaric acid (used in the production of detergents) from glucose by reenginering a synthetic pathway into the microbial cell [36]. However, the limitation of cellular resources yields a competition between the native pathways of the cell (dedicated to growth) and the synthetic circuits (dedicated to metabolite synthesis). Numerous works have studied this trade-off, aiming to improve the current industrial bioprocesses. Bacterial growth can be inhibited through growth inhibitors, thus increasing the yield of metabolites, as shown in recent experimental results [37]. While this is an effective technique, it can also be expensive, and can produce toxicity to the host cell [38]. A less costly, innovative tool is the optogenetic regulation of gene expression, which has been successfully implemented in a real-time feedback loop in continuous liquid cultures [39].

A promising alternative is the engineering of synthetic gene circuits. An example can be seen in [40], where a synthetic control circuit is able to compensate for the demand of resources from both the native and engineered pathways by regulating
the expression of enzymes. A method for engineering metabolic control circuits is proposed in [41], which tackles the burden of deviating resources from growth into the heterologous pathway through a multiobjective optimization problem, towards a framework for automated design of genetic circuits. Other approaches, based on coarse-grained models of the proteome, focuses on the loss of fitness of individuals in the population produced by the burden of the heterologous pathway, and on how to overcome such effect [42].

The capacity to externally affect gene expression through external real-time controllers, combined with the increasing understanding of natural resource allocation strategies, has motivated approaches based on self-replicator models [43]. The latter focuses on how to regulate the expression of RNA polymerase through a synthetic growth switch [44] in order to maximize the production of a metabolite of interest. The model represents the trade-off biomass/compound synthesis through a resource allocation problem similar to the one studied in Giordano et al., and the results emphasizes the differences between the control strategies depending on the particular objective.

Based on these works, we propose a series of coarse-grained self-replicator dynamical models representing the proteome of a bacterial cell, and we study the trade-offs arising from inherent natural allocation compromises, as well as those related to the artificial production of compounds of interest. The mathematical tools employed for this task range from dynamic systems analysis of the given models, to optimization and optimal control theory for studying the associated static and dynamical maximization objectives.

## 1.6 Studied problems

In this section, a brief description of the nature of the models used in this manuscript is given, all devised from a resource allocation perspective, that consider some of the most relevant industrial production frameworks. Each chapter is written in a self-contained manner, as they are based on (already published or in production) scientific papers, and thus there are substantial differences in the biological hypotheses and the overall notation. However, some common definitions can be provided to describe the organization of the manuscript.

The time-varying quantities represented in the models are divided into the extra-
<figure><img src="image_5.png" /><figcaption>Figure 1.3: Schema of the system.</figcaption></figure>

cellular concentrations (i.e. with respect to a bioreactor of volume $V_{ext}$) and intracellular concentrations (i.e. with respect to the bacterial volume $V$). The extracellular variables are $s$ and $x$, corresponding to the substrate and metabolite concentrations with respect to the external volume, respectively. The intracellular variables are $p$, $r$ and $m$, which correspond to the concentrations of precursor metabolites, proteins of the gene expression machinery, and proteins of the metabolic machinery, respectively. The main reaction rates are $v_M$, $v_R$ and $v_X$, and are functions of the concentrations. The precursor synthesis rate $v_M(s, m)$ is increasing with respect to the concentration of substrate in the medium $s$ and catalyzed by $m$. The macromolecule synthesis rate $v_R(p, r)$ is increasing with respect to the concentration of precursors $p$ and catalyzed by $r$. Finally, $v_X(p, m)$ is the metabolite synthesis rate, also increasing with respect to the concentration of precursors $p$, and catalyzed by $m$. A schema of the system is depicted in Figure 1.3. While the derivations of each particular model are detailed in each chapter throughout the manuscript, and thus they are omitted in this introduction, we describe here the obtained differential equations governing the time evolution of the defined quantities. The dynamics of the pool of precursors is described by the differential equation

$$ \dot{p} = v_M(s, m) - v_X(p, m) - v_R(p, r)(p+1). $$

The dynamics of the concentrations of macromolecules in the cell, described by
ribosomal and enzymatic proteins, are given by

$$ 
\begin{aligned}
\dot{r} &= (u - r)v_R(p, r), \\
\dot{m} &= (1 - u - m)v_R(p, r),
\end{aligned}
$$

where $u \in [0, 1]$ is the time-varying allocation strategy, which can represent either the natural allocation or an artificial human-driven allocation resulting from the effect of an external control, depending on the objective to be studied. The dynamics of the bacterial volume and the metabolite concentration are

$$ 
\begin{aligned}
\dot{\mathcal{V}} &= v_R(p, r)\mathcal{V}, \\
\dot{x} &= v_X(p, m)\frac{\mathcal{V}}{\mathcal{V}_{ext}}.
\end{aligned}
$$

The simplified case, where the dynamics of $s$ are omitted—representing the case where the substrate concentration is constant for the system—was first studied in [3], and is further investigated in Chapter 2. In this work, $s$ is assumed to be constant, and so the dynamics of $p$ becomes

$$ \dot{p} = v_M(m) - v_X(p, m) - v_R(p, r)(p + 1), $$

where $v_M$ only depends on $m$. Intuitively, this case is expected to undergo exponential bacterial growth after a certain transient. An extension of this work is presented in Chapter 5, where a hierarchical MPC (Model-Predictive Control) loop is proposed, based on the studied open-loop optimal control.

The case where an initial amount of substrate $s_0$ in the bioreactor is gradually consumed represents a batch process, and is studied in Chapter 3. The latter is modelled through the dynamics of $s$ given by the equation

$$ \dot{s} = -v_M(s, m)\frac{\mathcal{V}}{\mathcal{V}_{ext}}, $$

and is expected to eventually reach a steady state when all the available substrate in the bioreactor is depleted. A scheme of metabolite production in continuous bioreactors is studied in Chapter 4. In this case, there is a continuous constant inflow of substrate of concentration $s_{in}$ to the bioreactor, as well as an outflow (of the same volumetric flow rate) of substrate, biomass and metabolites. This flow is
characterized by the constant $D$ denoted as the dilution rate. Thus, the dynamics are modified with additional terms

$$ \dot{s} = Ds_{in} - v_{M}(s, m)\frac{V}{V_{ext}} - Ds, $$

$$ \dot{V} = v_{R}(p, r)V - DV, $$

$$ \dot{x} = v_{X}(p, m)\frac{V}{V_{ext}} - Dx. $$

The latter system is also expected to reach a steady state, depending on the dilution rate $D$ and on the allocation variable $u$. An experimental stage of this project, consisting of E. coli cultures in continuous bioreactors, was conducted in collaboration with Ibis team (Inria Grenoble - Rhône-Alpes). However, it was not possible to perform the required amount of experiments due to the pandemic, and so these results were not included in the manuscript.

Chapter 6 describes a self-replicator model that considers growth rate-independent proteins in the proteome. The latter aims to model a wild-type bacteria with no metabolite production, and so the synthesis rate of metabolites $v_X(p, m) = 0$. Growth rate-independent proteins are divided into ribosomal proteins and housekeeping proteins. The fact that there are proteins of the gene expression machinery that do not contribute to growth produces a minimal ribosomal concentration $r_{min}$ required to have bacterial growth. Thus, the synthesis rate of macromolecules becomes $v_R(p, r - r_{min})$. Additionally, since a part of the proteome is dedicated to housekeeping proteins $q$, there is a maximal ribosomal concentration $r_{max}$. The equations for the intracellular concentrations in the proteome are then given by

$$ \dot{r} = (r_{max}u - r)v_{R}(p, r - r_{min}), $$

$$ \dot{m} = (r_{max}(1 - u) - m)v_{R}(p, r - r_{min}), $$

$$ \dot{q} = ((1 - r_{max}) - q)v_{R}(p, r - r_{min}). $$

Finally, an additional section in Appendix A describes supplementary work carried out on gene-regulatory networks, produced as part of a collaborative work during my PhD. The results fall within the realm of Systems Biology, but they have been excluded from the main body of the manuscript for being out of the specific scope of the thesis.
# Chapter 2

## Constant substrate inflow

This chapter reproduces [AYT5], published in the 58th Conference on Decision and Control (IEEE CDC 2019), and some unpublished results on the stability of the system of interest and the uniqueness of the solution of the static optimization problem.

### 2.1 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 [14]. 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 [21]. 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 [45].

The starting point for this line of research is [30], 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 [43], 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 [30, 43]. We
start, in Section 2.2, 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 2.4, 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 com-
parison between this two problems should help to understand how to dynamically
disrupt the natural allocation process in order to prioritize the metabolite produc-
tion 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 opti-
mally 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, 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 2.5 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.

## 2.2 Self-replicator models

### 2.2.1 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 2.1a. Based on the extension introduced in [43], a metabolic pathway for the production of a certain metabolite of interest $X$ is included (Figure 2.1b).

<figure><img src="image_6.png" /><figcaption>Figure 2.1: Self-replicator models of bacterial growth.</figcaption></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$^{-1}$]. The system is subject to a constant volumetric flow rate F [L h$^{-1}$] which generates both an inflow of fresh medium rich in substrate, and an outflow of biomass and metabolites [46]. Moreover, the scheme is extended with the growth switch described in [44] 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 (2.1)$$

where I is the external control and $\alpha$ the natural allocation mechanism used in [30]. 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.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$^{-1}$] being the concentration of the nutrient input, and D [h$^{-1}$] 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 (2.3)$$

where $\beta$ [L g⁻¹] corresponds to the inverse of the cytoplasmic density. The definition (2.3) is based on the experimental fact that macromolecules explain most of the biomass in microbial cells [6]. 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 = \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(s, m)\frac{V}{V_{ext}}, \\
\dot{p} &= v_M(s, m) - v_R(p, r) - v_X(p, m) - \mu(t)p, \\
\dot{r} &= uv_R(p, r) - \mu(t)r, \\
\dot{m} &= (1 - u)v_R(p, r) - \mu(t)m, \\
\dot{x} &= v_X(p, m) - \mu(t)x \\
\dot{V} &= (\mu(t) - D)V,
\end{aligned}
\right.
\qquad (2.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(t)$ [h⁻¹] is the growth rate of the self-replicator system that, using (2.2), is defined as

$$\frac{\dot{V}}{V} = \frac{\dot{M} + \dot{R}}{M + R} = \beta v_R(p, r) - D \quad \rightarrow \quad \mu(t) \doteq \frac{\dot{V}}{V} \bigg|_{F=0} = \beta v_R(p, r). \qquad (2.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(s, m) &\doteq k_M m \frac{s}{K_M + s}, \\
v_R(p, r) &\doteq k_R r \frac{p}{K_R + p}, \\
v_X(p, m) &\doteq k_X m \frac{p}{K_X + p},
\end{aligned} 
\qquad (2.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, & \hat{p} &\doteq \beta p, & \hat{r} &\doteq \beta r, & \hat{x} &\doteq \beta x & \hat{X} &\doteq \beta X, \\
K &\doteq \beta K_R, & K_1 &\doteq \beta K_X, & K_2 &\doteq \beta K_M \\
k_1 &\doteq \frac{k_X}{k_R}, & k_2 &\doteq \frac{k_M}{k_R}.
\end{aligned}
$$

By replacing all variables and dropping all hats, system (2.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} x, \\
\dot{V} &= \left( \frac{pr}{K + p} - D \right) V.
\end{aligned}
\right. 
\qquad (2.7) $$

where the dynamical expression of $m$ has been removed since $m = 1 - r$. The analysis of model (2.7) is part of a more recent work [AYT2], which is described in this thesis in Chapter 4. In the next section, we present the particular case where $D = 0$, meaning that there is no substrate inflow, and no biomass output.
### 2.2.2 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 (2.8) $$

The optimal product maximization problem was partly analyzed for this particular case in [43] where, due to the complexity of the computations in the PMP, most of the analysis was performed through a numerical approach. In Chapter 3, a more detailed study of this model and its associated optimal control problems is performed. To allow an analytical study, we 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, which yields the model of interest in this paper.

### 2.2.3 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 (2.6). A constant value

$$ E_M \doteq k_2 \frac{s}{K_2+s} $$

is defined, such that model (2.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 (S) $$

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.

### 2.2.4 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 2.1a. This is a particular case of the model (S), when $k_1 = 0$, and the external control $I(t)$ introduced in (2.1) is overridden (so that $u(t) = \alpha(t)$),

$$ \begin{cases} \dot{p} = E_M(1-r) - (p+1)\frac{pr}{K+p}, \\ \dot{r} = (\alpha-r)\frac{pr}{K+p}. \end{cases} \qquad (2.9) $$

The problem of biomass maximization by natural mechanisms has been extensively analyzed for this model in [30], 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.

## 2.3 Dynamical analysis of model (S)

In this section, we perform a local and global analysis of model (S). Given the nullcline $p$ that satisfies the equation

$$ 0 = E_M(1-r) - w_X(p)(1-r) - w_R(p)(p+1)r, $$

with

$$ w_R(p) = \frac{p}{K+p}, \quad w_X(p) = \frac{k_1p}{K_1+p}, $$

we can define the function

$$ r_n(p) \doteq \frac{E_M - w_X(p)}{E_M - w_X(p) + w_R(p)(p+1)} \in (0, 1], \qquad (2.10) $$
which is monotonically decreasing w.r.t. $p$, such that the equation for $p$ can be rewritten as

$$ \frac{dp}{dt} = (r_n(p) - r)[E_M - w_X(p) + w_R(p)(p+1)], $$

leading to the system

$$ \left\{ \begin{aligned} \frac{dp}{dt} &= (r_n(p) - r)[E_M - w_X(p) + w_R(p)(p+1)], \\ \frac{dr}{dt} &= (u - r)w_R(p)r, \\ \frac{dx}{dt} &= w_X(p)(1 - r) - w_R(p)rx, \\ \frac{dV}{dt} &= w_R(p)rV. \end{aligned} \right. \quad (S_1) $$

### 2.3.1 Global behavior

Given that the volume $V$ admits no steady state except the trivial $V = 0$, we will study the dynamical behavior of the system without considering this variable. Additionally, since the dynamics of variables $(p, r)$ do not depend on $x$, we will resort to well-known decomposition techniques presented in [47] to study the stability properties of $S_1$ through its subsystems. More precisely, we will separate the system into two subsystems arranged in a hierarchical form

$$ \begin{aligned} z_1 &= \begin{bmatrix} p \\ r \end{bmatrix}, \quad \frac{dz_1}{dt} = f_1(z_1), \\ z_2 &= \begin{bmatrix} x \end{bmatrix}, \quad \frac{dz_2}{dt} = f_2(z_1, z_2), \end{aligned} $$

where $f_1$ depends only on $z_1$. Therefore, in the present section, we will study the stability of the subsystem given by variables $(p, r)$

$$ \left\{ \begin{aligned} \frac{dp}{dt} &= (r_n(p) - r)[E_M - w_X(p) + w_R(p)(p+1)] \\ \frac{dr}{dt} &= (u - r)w_R(p)r \end{aligned} \right. \quad (S'_1) $$
to conclude, towards the end of the section, on the stability of the global system $(S_1)$.

**Lemma 2.3.1.** *The set,*

$$ \Gamma \doteq \{(p, r) \in \mathbb{R}^2 : p \ge 0, 1 \ge r \ge 0\}, $$

is positively invariant in $(S'_1)$ for the initial value problem.

*Proof.* Let us analyze the boundaries of $\Gamma$,

$$ \dot{p}|_{p=0} = E_M(1-r) \ge 0 \quad \Rightarrow \quad p=0 \text{ is repulsive or invariant.} $$

$$ \dot{r}|_{r=0} = 0 \quad \Rightarrow \quad r=0 \text{ is invariant.} $$

$$ \dot{r}|_{r=1} = (u-1)w_R(p) \le 0 \quad \Rightarrow \quad r=1 \text{ is repulsive or invariant.} $$

□

In order to study the steady-state behavior of the system, we fix $u(t) = \bar{u} \in (0, 1)$.

**Lemma 2.3.2.** *The system admits two equilibria:*

* *The interior equilibrium:* $E_i = (p_i, \bar{u})$

* *The no-growth equilibrium:* $E_w = (0, 1)$

with $p_i$ solution of the equation

$$ E_M(1 - \bar{u}) - (p+1)w_R(p)\bar{u} - w_X(p)(1 - \bar{u}) = 0 \quad (2.11) $$

*Proof.* To prove uniqueness of $E_i$, it is sufficient to prove that the value of $p_i$ is unique in (2.11). We define the constant value,

$$ c_1 \doteq E_M(1 - \bar{u}) \ge 0, $$

replace them in (2.11) and define the function $f$ as

$$ f(p) \doteq w_X(p)(1 - \bar{u}) + (p+1)w_R(p)\bar{u}, $$
with the following properties

$$f(0) = 0, \quad \lim_{p \to \infty} f(p) = \infty, \quad \frac{\partial f(p)}{\partial p} > 0.$$

Therefore, $f(p) = c_1$ has a unique solution $p_i$, which also leads to a unique value $x_i$ and thus the equilibrium $E_i$ is unique. □

### 2.3.2 Local stability of equilibria

We first write the Jacobian matrix

$$J = \begin{bmatrix} -(p+1)w'_{R}(p)r - w_{R}(p)r - w'_{X}(p)(1-r) & -E_{M} - w_{R}(p)(p+1) + w_{X}(p) \\ w'_{R}(p)r(\bar{u} - r) & -w_{R}(p)r + (\bar{u} - r)w_{R}(p) \end{bmatrix}, \quad (2.12)$$

where

$$w'_{X}(p) \doteq \frac{\partial w_{X}(p)}{\partial p}, \quad w'_{R}(p) \doteq \frac{\partial w_{R}(p)}{\partial p}.$$

#### Equilibrium $E_i$

Replacing the values of $E_i$ in (2.12)

$$J_i = \begin{bmatrix} -(p+1)w'_{R}(p)\bar{u} - w_{R}(p)\bar{u} - w'_{X}(p)(1-\bar{u}) & -E_{M} - w_{R}(p)(p+1) - w_{X}(p) \\ 0 & -w_{R}(p)\bar{u} \end{bmatrix},$$

with eigenvalues

$$\lambda = -(p+1)w'_{R}(p)\bar{u} - w_{R}(p)\bar{u} - w'_{X}(p)(1-\bar{u}), -w_{R}(p)\bar{u},$$

and so the equilibrium is locally stable.
## Equilibrium $E_w$

In $E_w$, one has $w_R(p) = w_X(p) = 0$ and $w'_R(p)r = \frac{1}{K}$, and so the Jacobian matrix becomes

$$J_w = \begin{bmatrix} -\frac{1}{K} & -E_M \\ -\frac{1}{K}(1-\bar{u}) & 0 \end{bmatrix},$$

with characteristic polynomial

$$P(\lambda) = \left(\lambda + \frac{1}{K}\right) \lambda - E_M \frac{1}{K} (1 - \bar{u}) = \lambda^2 + \lambda \frac{1}{K} - E_M \frac{1}{K} (1 - \bar{u}).$$

Then, the eigenvalues are given by

$$\lambda_1 = \frac{-\frac{1}{K} - c_1}{2} < 0, \quad \lambda_2 = \frac{-\frac{1}{K} + c_1}{2} > 0, \quad \text{with } c_1 = \sqrt{\frac{1}{K}^2 + 4E_M \frac{1}{K}(1-\bar{u})}$$

showing that the equilibrium is a saddle point. Now, let us find the eigenvector $v_1 = [p_1, r_1]^T$ associated to the stable eigenvalue $\lambda_1 < 0$,

$$\begin{bmatrix} \lambda_1 + \frac{1}{K} & E_M \\ \frac{1}{K}(1-\bar{u}) & \lambda_1 \end{bmatrix} \begin{bmatrix} p_1 \\ r_1 \end{bmatrix} = 0, \quad \Rightarrow \quad \begin{cases} (\lambda_1 + w'_R(p)r)p_1 + E_M r_1 = 0 \\ w'_R(p)r(1-\bar{u})p_1 + \lambda_1 r_1 = 0 \end{cases}$$

Solving the system yields the relations

$$p_1 = -r_1 \frac{E_M}{\lambda_1 + w'_R(p)r}, \quad p_1 = -r_1 \frac{\lambda_1}{w'_R(p)r(1-\bar{u})},$$

where it can be seen that both equations are equivalent

$$\frac{E_M \frac{1}{K}(1-\bar{u})}{(\lambda_1 + \frac{1}{K})\lambda_1} = \frac{4E_M \frac{1}{K}(1-\bar{u})}{-(\frac{1}{K}-c_1)(\frac{1}{K}+c_1)} = \frac{4E_M \frac{1}{K}(1-\bar{u})}{(\frac{1}{K})^2 + 4E_M \frac{1}{K}(1-\bar{u}) - (\frac{1}{K})^2} = 1,$$

and so, since $\lambda_1 < 0$, $p_1$ and $r_1$ have the same sign, which means that the eigenvector associated to the stable eigenvalue points to the first and third quadrants, which is outside the invariant set $\Gamma$.
### 2.3.3 Global analysis of $S'_1$

The analysis of the global behavior consists on dividing $\Gamma$ in four regions delimited by the system nullclines (as seen in Figure A.2) and studying how the system behaves in each region separately. We start by stating what has been showed in Section 2.3.2:

**Lemma 2.3.3.** A trajectory originating in $\Gamma \setminus \{E_w\}$ can't converge to the unstable equilibria $E_w$.

*Proof.* As shown in Section 2.3.2, $E_w$ is a saddle point, and the eigenvector associated to the stable eigenvalue points to the first and third quadrant w.r.t. $E_w = (0, 1)$. As this is the only possible direction through which the saddle point can be attained, the equilibrium is not an attractor in $\Gamma$, except for trajectories starting in $E_w$. $\square$

<figure><img src="image_7.png" /><figcaption>Figure 2.2: Phase plane of system $S'_1$.</figcaption></figure>

**Lemma 2.3.4.** Every solution of the initial value problem starting in the sets,

$$ \Gamma_A \doteq \{(p, r) \in \mathbb{R}^2 : p_i \ge p > 0, r_n(p) \ge r \ge \bar{u}\}, $$

$$ \Gamma_C \doteq \{(p, r) \in \mathbb{R}^2 : p \ge p_i, \bar{u} \ge r \ge r_n(p)\}, $$

converges asymptotically to the equilibrium $E_i$.
*Proof.* First, let us show that $\Gamma_A$ is positively invariant to the system

$$ \dot{p}|_{p=0} = E_M(1-r) > 0, \forall r \in (\bar{u}, 1) \quad \Rightarrow \quad p=0 \text{ is repulsive.} $$

$$ \dot{r}|_{r=r_n(p)} = (\bar{u} - r_n(p))w_R(p)r_n(p) < 0, \forall p \in (0, p_i) \quad \Rightarrow \quad r = r_n(p) \text{ is repulsive.} $$

$$ \dot{r}|_{r=\bar{u}} = 0 \quad \Rightarrow \quad r = \bar{u} \text{ is invariant.} $$

as well as $\Gamma_C$

$$ \dot{r}|_{r=r_n(p)} = (\bar{u} - r_n(p))w_R(p)r_n(p) > 0, \forall p > p_i \quad \Rightarrow \quad r = r_n(p) \text{ is repulsive.} $$

$$ \dot{r}|_{r=\bar{u}} = 0 \quad \Rightarrow \quad r = \bar{u} \text{ is invariant.} $$

Then, asymptotic stability of $E_i$ can be proved using the following simple quadratic Lyapunov function

$$ V_{A,C}(t) \doteq \frac{1}{2}(p_i - p)^2 + \frac{1}{2}(\bar{u} - r)^2 > 0, \quad \forall (p, r) \in (\Gamma_A \cup \Gamma_C) \setminus \{E_i\}, $$

such that

$$ \frac{d}{dt} V_{A,C}(t) = -(p_i - p)\dot{p} - (\bar{u} - r)^2 w_R(p) < 0, \quad \forall (p, r) \in (\Gamma_A \cup \Gamma_C) \setminus \{E_i\}, $$

which verifies the inequality given that in $\Gamma_A$, $p_i > p$ and $\dot{p} > 0$, and in $\Gamma_C$, $p_i < p$ and $\dot{p} < 0$. The equilibrium $E_i$ is the only point where the derivative $\frac{d}{dt}V_{A,C}(t) = 0$, as well as the largest invariant set, which confirms through LaSalle theorem [48] that every solution starting in $\Gamma_A$ and $\Gamma_C$ approaches $E_i$ as $t \to \infty$. $\square$

**Lemma 2.3.5.** Every solution of the initial value problem starting in the sets,

$$ \Gamma_B \doteq \{(p, r) \in \mathbb{R}^2 : p > 0, 1 \ge r > \max\{\bar{u}, r_n(p)\}\}, $$

$$ \Gamma_D \doteq \{(p, r) \in \mathbb{R}^2 : p > 0, \min\{\bar{u}, r_n(p)\} > r > 0\}, $$

converges asymptotically to the equilibrium $E_i$.

*Proof.* Every point in $\Gamma_B$ meets $r_n(p) < r$ and $r > \bar{u}$, meaning that the vector field has signs $\dot{p} < 0$, $\dot{r} < 0$. Analogously, every point in $\Gamma_D$ meets $r_n(p) > r$ and $r < \bar{u}$, so the vector field has signs $\dot{p} > 0$, $\dot{r} > 0$. As stated in Lemma 2.3.3, no trajectory can converge towards $E_w$, so it can go whether towards the equilibrium $E_i$, or towards any of both nullclines, where eventually it will converge to $E_i$ as
proven in Lemma 2.3.4.

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

Lemmas 2.3.4 and 2.3.5 lead us to the following theorem.

**Theorem 2.3.1.** Every solution of the initial value problem with initial conditions $p(0) > 0$ and $r(0) \in (0, 1)$ converges asymptotically to the equilibrium $E_i$ as $t \to \infty$.

**Connection with the original system**

As the global stability of the equilibrium $E_i$ in subsystem $S'_1$ has been established, we can conclude on the global behavior of the original system $S_1$. For that, we will state a simplified version of the theorem presented in [47] that relates the stability of a hierarchical system

$$S_z : \begin{cases} \frac{dz_1}{dt} = f_1(z_1) \\ \frac{dz_2}{dt} = f_2(z_1, z_2) \end{cases}$$

with equilibrium $\bar{z} = [\bar{z}_1, \bar{z}_2]$, with the stability of its so-called isolated subsystems

$$S_{z_1} : \frac{dz_1}{dt} = f_1(z_1), \quad S_{z_2} : \frac{dz_2}{dt} = f_2(\bar{z}_1, z_2),$$

where, in the second system $S_{z_2}$, $z_1$ is set to $\bar{z}_1$.

**Theorem 2.3.2.** [47] Given $f_1$ and $f_2$ continuous, autonomous and continuously differentiable with respect to all variables, satisfying

$$\sup_{t>0} \sup_{z_1} ||\nabla_{z_1} f_1(z_1)|| < \infty, \\ \sup_{t>0} \sup_{(z_1,z_2)} ||\nabla_{(z_1,z_2)} f_2(z_1,z_2)|| < \infty,$$

and $f_1(\bar{z}) = f_2(\bar{z}) = 0$, then $z = \bar{z}$ is a globally asymptotically stable equilibrium point of $S_z$ if and only if $z_1 = \bar{z}_1$ and $z_2 = \bar{z}_2$ are globally asymptotically stable equilibrium points of isolated subsystems $S_{z_1}$ and $S_{z_2}$ respectively.

In the present case, as the asymptotic stability of the equilibrium $E_i$ in $S'_1$ has been already proven, and given that all terms in the Jacobian matrix (2.12) are
bounded for every value of $\Gamma$, we proceed to evaluate the dynamical expression of $x$ in $(p, r) = E_i$

$$ \frac{dx}{dt} = w_X(p_i)(1 - \bar{u}) - w_R(p_i)\bar{u}x, $$

showing that the $x = x_i$ is globally asymptotically stable, with

$$ x_i \doteq \frac{w_X(p_i)(1 - \bar{u})}{w_R(p_i)\bar{u}}. $$

The latter is formalized in the following theorem.

**Theorem 2.3.3.** In every solution of the initial value problem in system $S_1$ starting in $p(0) > 0$, $r(0) \in (0, 1)$, $x(0) \ge 0$ and $\mathcal{V}(0) > 0$, when $t \to \infty$,

* $(p, r)$ converge asymptotically to the equilibrium $E_i$.

* The biomass volume $\mathcal{V}$ grows exponentially with rate

$$ \mu(p_i, \bar{u}) = w_R(p_i)\bar{u} = \frac{p_i\bar{u}}{K + p_i} $$

* The quantity $x$ converges asymptotically to

$$ x_i = \frac{w_X(p_i)(1 - \bar{u})}{w_R(p_i)\bar{u}} $$

as the total mass of metabolite X also grows exponentially.

## 2.4 Optimal control problem

### 2.4.1 Dynamical biomass maximization in Model (2.9)

In model (2.9), 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(p, r) \, dt \\
\text{subject to} & \alpha(\cdot) \in \mathcal{U}
\end{array}
$$

with

$$
\mu(p, r) = \frac{pr}{K + p}
$$

being the growth rate defined in (2.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, \quad (2.13)
$$

where *F* denotes the right-hand side of (2.9) 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{2.14}
$$

being

$$
\phi(\cdot) = \lambda_r \frac{pr}{K+p}
$$

the switching function.

2.4.2 Product maximization in Model (S)

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 (S), the criterion can be written as

$$J_X(u) = \int_0^T k_1 \frac{p(1-r)}{K_1 + p} \mathcal{V} 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. \qquad (2.15)$$

## Static problem

We are interested in the steady state that maximizes the metabolite production in terms of the constant allocation $\bar{u}$. The latter is represented by the integrand of the cost function (2.15). The problem can be stated as

$$\begin{array}{ll} \text{maximize} & J(\bar{u}, p_i(\bar{u})) = \frac{w_X(p)(1 - \bar{u})}{x_i}, \\ \text{subject to} & \bar{u} \in (0, 1). \end{array} \qquad (2.16)$$

While an explicit calculation of the solution of the latter problem is rather hard to obtain, we will focus on showing the existence and uniqueness of such solution. First, we notice that the cost can be equivalently rewritten as

$$J(\bar{u}, p_i(\bar{u})) = w_R(p)\bar{u} = \frac{p_i(\bar{u})\bar{u}}{K + p_i(\bar{u})},$$

which shows that maximizing the metabolite production at steady state is equivalent to maximizing the biomass production. Recalling the function $r_n(p)$ defined in (2.10), we can express $\bar{u}$ in terms of $p$ and maximize $J(p, \bar{u}(p))$ with respect to $p$ (as the steady-state value $p_i(\bar{u})$ is strictly decreasing w.r.t. $\bar{u}$). Moreover, for the sake of simplicity, rather than maximizing the growth rate $w_R(p)$, we will minimize its inverse, which yields the optimization problem

$$\begin{array}{ll} \text{minimize} & \hat{J}(p, \bar{u}(p)) \doteq \frac{1}{J(p, \bar{u}(p))} = \frac{K+p}{p \bar{u}(p)}, \\ \text{subject to} & p \in (0, \infty). \end{array} \qquad (2.17)$$
The objective can be further developed by replacing $\bar{u}(p)$ by its expression, which yields

$$ \hat{J}(p) = \frac{K}{p} + 1 + \frac{(p+1)(K_1 + p)}{E_M K_1 + p(E_M - k_1)}. $$

The following assumption is required to ensure the well-posedness of the problem:

**Assumption 2.4.1.** $E_M > k_1$.

The latter assumption implies that the maximal substrate intake is strictly greater than the maximal metabolite production. Note that this assumption also includes model (2.9), where $k_1 = 0$, as a particular case of the current analysis. For notation purposes, we define the constants

$$ c_1 \doteq \frac{E_M K_1}{E_M - k_1} > 0, \quad c_2 \doteq \frac{(1 - c_1)(K_1 - c_1)}{E_M - k_1}, \qquad (2.18) $$

such that

$$ K_1 - c_1 = K_1 \left( 1 - \frac{E_M}{E_M - k_1} \right) < 0. $$

Then, the objective function becomes

$$ \hat{J}(p) = \frac{K}{p} + 1 + \frac{p+1}{E_M - k_1} + \frac{K_1 - c_1}{E_M - k_1} + \frac{c_2}{p + c_1}. \qquad (2.19) $$

It can be seen from (2.19) that, at the boundaries, the function $\hat{J}(p)$ verifies

$$ \lim_{p \to 0} \hat{J}(p) = \infty, \quad \lim_{p \to \infty} \hat{J}(p) = \infty, $$

which is coherent with the fact that the growth rate $\mu(p, r) = 0$ at the boundaries. The latter, and the continuity of $\hat{J}(p)$, implies the existence of, at least, one solution. Additionally, it means that the optimal value $p$ is interior, so in order to prove uniqueness of the solution, it suffices to show that the derivative is cancelled in only one value of $p$ by solving

$$ \frac{\partial}{\partial p} \hat{J}(p) = -\frac{K}{p^2} + \frac{1}{E_M - k_1} - \frac{c_2}{(p + c_1)^2} = 0. \qquad (2.20) $$

For that, there are two possible cases to explore in terms of the sign of $c_2$.
**Case $c_2 \ge 0$**

It can be seen from (2.18) that if $1 - c_1 \le 0$, then $c_2 \ge 0$ and so we can rewrite (2.20) as

$$ \frac{1}{E_M - k_1} = \underbrace{\frac{K}{p^2} + \frac{c_2}{(p + c_1)^2}}_{f_1(p)} $$

where

$$ \lim_{p \to 0} f_1(p) = \infty, \quad \lim_{p \to \infty} f_1(p) = 0, \quad \frac{\partial}{\partial p} f_1(p) = -\frac{2K}{p^3} - \frac{2c_2}{(p + c_1)^3} < 0, $$

which means that the equation $f_1(p) = \frac{1}{E_M - k_1}$ has a unique solution. Then, there is a single value of $p$ that cancels (2.20), and so the minimization problem (2.17) has a unique solution.

**Case $c_2 < 0$**

In the case where $c_2 < 0$ due to $1 - c_1 > 0$, we can rewrite the equality (2.20) as

$$ K = \underbrace{\frac{p^2}{E_M - k_1} - c_2 \left( \frac{p}{p + c_1} \right)^2}_{f_2(p)} $$

where

$$ f_2(0) = 0, \quad \lim_{p \to \infty} f_2(p) = \infty, \quad \frac{\partial}{\partial p} f_2(p) = \frac{2p}{E_M - k_1} - \frac{2c_2 p}{(p + c_1)^2} \left( 1 - \frac{p}{p + c_1} \right) > 0, $$

showing that the equation $f_2(p) = K$ has a unique solution, which again means that the minimization problem (2.17) has a unique solution. The latter can be formalized in the following theorem:

**Theorem 2.4.1.** Under assumption 2.4.1, there exists a unique constant allocation $\bar{u} \in (0, 1)$ that maximizes the steady-state metabolite production of Problem (2.16) (equivalent to the steady-state biomass production).
## Dynamic problem

For the sake of convenience, we perform a change of variables $y = \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 (2.21)$$

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}, \\
\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}, \\
\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 (2.14), with the same switching function $\phi(\cdot) = \lambda_r \frac{pr}{K+p}$.

### 2.4.3 Characterization of singular regimes

Singular trajectories play a major role in optimal control theory [49, 50], and their characterization is a necessary step towards a complete description of the optimal solution. In both explored OCPs, and as expected in linear optimization 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 2.4.1, vanishing of the switching function implies that $\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 [30] that the singular arc is in fact the optimal steady state solution of the static problem, 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 [51], i.e. an arc with an infinite number of switchings. We now focus on the model (S) and provide a detailed computation of the corresponding singular arcs. The Hamiltonian (2.21) can be expressed as

$$H_B = H_0 + u H_1$$

where

$$H_0 = 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 + \lambda_y)pr}{K + p}, \\ H_1 = \phi(t) = \frac{\lambda_r pr}{K + p}.$$

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} | 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{aligned} \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{aligned}$$
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{2.22}
\end{equation}
$$

and evaluating (2.22) 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 (2.23)
$$

From (2.23) 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 Σ given by

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

Similarly, the second derivative can be computed as

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

and so one has

$$
\begin{equation}
\begin{split}
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{split}
\tag{2.24}
\end{equation}
$$
where

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

Evaluating in $z \in \Sigma'$ yields

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

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)} \\ - \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 + p} - 1 \right)}{K + p} \right) \left( E_m - \frac{k_1 p}{K_1 + p} \right) \frac{p}{K + p},$$

which defines the set

$$\Sigma'' \doteq \{(\varphi, \lambda) \in \mathbb{R}^{2n} \mid 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_{00001} + 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. Having a singular control of intrinsic order two (that is such that $H_{101}$ vanishes identically on the whole cotangent space) implies that singular arcs can only be entered and exited through chattering (Fuller phenomenon); see [69]. Although the singular is only of local order two in our case ($H_{101}$ vanishes when $H_1$ and $H_{01}$ do, not identically), the numerical simulations below show that there is indeed chattering in and out to enter and exit the singular arc.

## 2.5 Numerical solution

The numerical verifications were performed with Bocop [52], 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$^1$. 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 characterized by a singular arc of order two that is entered and exited through chattering (Figure 2.3). Moreover, the solution for the biomass maximization problem $\alpha_{opt}$ matches the optimal steady-state input $\alpha_{opt}^*$ along the circular arc (Figure 2.3a), 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 2.3b) 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 lower value $u(t)$ depicted in Figure

<footnote>$^1$Bocop definition files are available from the authors on request.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 2.3: Solutions of the OCPs computed using *Bocop*.</figcaption></figure>

2.3b. 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, which corresponds to the solution of the static problem, already shown to exist and to be unique in previous section. 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 2.4). 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.

## 2.6 Conclusions

Self-replicator models of bacterial growth are capable of reproducing the growth laws of certain organisms at steady state. However, in [30], 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 [43], 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. Con-
<figure><img src="image_10.png" /><figcaption>Figure 2.4: Simulation of the OCP 2.4.2 (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>

trary 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.
# Chapter 3

## Substrate depletion

This chapter describes unpublished work in progress.

### 3.1 Introduction

The study of living microorganisms through resource allocation models has increas-
ingly become relevant for its capacity to elucidate natural behaviors of microbia
through very simple dynamical models. The core idea is to represent the distribution
of internal resources through optimal control strategies, based on the assumption
that evolutionary processes have tuned these endogenous allocation strategies to
attain nearly-optimal levels. Numerous problems arise in this context, one of them
being the optimal production of metabolites regulated by an external control capa-
ble of arresting the bacterial growth. To this end, a resource allocation approach
can help understand how to modify the naturally-evolved allocation strategies so as
to efficiently produce such metabolite.

This is the subject of this paper, which tackles the problem of batch processing
from a resource allocation perspective. A simple coarse-grained self-replicator model
is introduced, based on [3] and [AYT5], with minimal biological assumptions. Using
mass conservation laws, it is possible to analyse the asymptotic behavior and stabil-
ity of the dynamical system, showing that for every possible allocation strategy, all
component of the system are transformed either into proteins or into metabolites, a
condition later defined as *Full depletion*. Then, two main studies are performed: the
biomass maximization case, representing the natural objective of wild-type micro-
bial cultures; and the metabolite maximization case, using a bacterial model that
includes a pathway for metabolite synthesis for industrial purposes. Both problems
are analyzed in infinite time and in finite time, the latter stated as an Optimal
Control Problem (OCP) which is solved through the application of the Pontryagin's
Maximum Principle (PMP). The problem presented here is similar to the one stud-
ied in [3], but approached from an analytical perspective, and including the case
with no metabolite synthesis as a starting point.

**3.2 Model definition**

**3.2.1 Self-replicator model**

We define a self-replicator model describing the dynamics of a microbial population
growing inside a closed bioreactor of fixed volume $V_{ext} > 0$ liters. At the beginning
of the experience, there is an initial mass of substrate S inside the bioreactor, that
is gradually consumed by the bacterial population, transforming it into precursor
metabolites P. Precursors are transformed into components of the gene expression
machinery R, into enzymes that makes up the metabolic machinery M, and into
a metabolite of interest X that is excreted from the cell. While the production of
proteins M and R is catalyzed by the enzymes R, the absorption of S and synthesis
of X are both catalyzed by M (dashed arrows in Figure 3.1). The proportion of

<figure><img src="image_11.png" /><figcaption>Figure 3.1: Self-replicator model</figcaption></figure>
precursors dedicated to the production of M and R is decided by the allocation control $u(t) \in [0, 1]$: $u = 0$ means no production of R, while $u = 1$ means no production of M. The control u can represent the natural allocation used by bacteria, as modelled in [30], as well as the artificially modified allocation modelled in [3], depending on the objective to be analyzed. The dynamics of the system are described by

$$ 
\left\{ 
\begin{aligned}
\dot{S} &= -V_M, \\
\dot{P} &= V_M - V_X - V_R, \\
\dot{M} &= (1-u)V_R, \\
\dot{R} &= uV_R, \\
\dot{X} &= V_X,
\end{aligned}
\right. 
$$

where the variables $S(t)$, $P(t)$, $M(t)$, $R(t)$ and $X(t)$ represent the masses (in grams) of substrate, precursors metabolites, the metabolic machinery, the gene expression machinery and the metabolites of interest at time $t$, respectively; $V_M$ [$g/h$], $V_R$ [$g/h$] and $V_X$ [$g/h$] are the reaction rates of the system, and $u(t)$ is the allocation control previously defined. We define the structural volume of the bacterial population inside the bioreactor $V_L(t)$ measured in liters as

$$ V_L \doteq \beta(M + R), \qquad (3.1) $$

where $\beta$ is a constant relating density and volume. This allows us to define concentrations with respect to this volume

$$ p_v \doteq \frac{P}{V_L} \left[ \frac{g}{L} \right], \quad m_v \doteq \frac{M}{V_L} \left[ \frac{g}{L} \right], \quad r_v \doteq \frac{R}{V_L} \left[ \frac{g}{L} \right], \qquad (3.2) $$

and define relative rates involved in the processes as

$$ v_M(s_v, m_v) \doteq \frac{V_M}{V_L} \left[ \frac{g}{Lh} \right], \quad v_R(p_v, r_v) \doteq \frac{V_R}{V_L} \left[ \frac{g}{Lh} \right], \quad v_X(p_v, m_v) \doteq \frac{V_X}{V_L} \left[ \frac{g}{Lh} \right]. $$

Likewise, we define the concentrations related to the external volume

$$ s_v = \frac{S}{V_{ext}} \left[ \frac{g}{L} \right], \quad x_v = \frac{X}{V_{ext}} \left[ \frac{g}{L} \right]. $$
From (3.1) and (3.2) we have that

$$m_{\nu} + r_{\nu} = \frac{M + R}{\mathcal{V}_L} = \frac{1}{\beta}. \qquad (3.3)$$

We define the growth rate $\mu = \mu(t)$ as

$$\mu \doteq \frac{\dot{\mathcal{V}}_L}{\mathcal{V}_L} = \beta v_R(p_\nu, r_\nu).$$

and that the dynamical system can be expressed in terms of the concentrations as

$$\left\{
\begin{aligned}
\dot{s}_{\nu} &= -v_M(s_{\nu}, m_{\nu}) \frac{\mathcal{V}_L}{\mathcal{V}_{ext}}, \\
\dot{p}_{\nu} &= v_M(s_{\nu}, m_{\nu}) - v_X(p_{\nu}, m_{\nu}) - v_R(p_{\nu}, r_{\nu})(\beta p_{\nu} + 1), \\
\dot{r}_{\nu} &= (u - \beta r_{\nu}) v_R(p_{\nu}, r_{\nu}), \\
\dot{m}_{\nu} &= (1 - u - \beta m_{\nu}) v_R(p_{\nu}, r_{\nu}), \\
\dot{\mathcal{V}}_L &= \beta v_R(p_{\nu}, r_{\nu}) \mathcal{V}_L, \\
\dot{x}_{\nu} &= v_X(p_{\nu}, m_{\nu}) \frac{\mathcal{V}_L}{\mathcal{V}_{ext}}.
\end{aligned}
\right.$$

### 3.2.2 Kinetics definition

We model the kinetics of the system by supposing they are linear in $m_\nu$ and $r_\nu$ [2]. Thus, they can be expressed as

$$v_M(s_\nu, m_\nu) = w_M(s_\nu)m_\nu,$$

$$v_R(p_\nu, r_\nu) = w_R(p_\nu)r_\nu,$$

$$v_X(p_\nu, m_\nu) = \gamma w_R(p_\nu)m_\nu,$$

where $\gamma > 0$ is a proportionality constant, under the assumption that $v_X(p_\nu, m_\nu)/v_R(p_\nu, r_\nu)$ does not depend on $p_\nu$ but on the ratio $m_\nu/r_\nu$. The functions $w_i$ are assumed to have the following behaviour:

**Assumption 3.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 monotonically increasing: $w'_i(x) > 0, \forall x \ge 0$,

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

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

For numerical simulations, we will resort to the particular case where the functions follow Michalis-Menten kinetics. For that case, we define

$$w_R(p_V) \doteq k_R \frac{p_V}{K_R + p_V}, \quad w_M(s_V) \doteq k_M \frac{s_V}{K_M + s_V},$$

where the values of the constants $k_R$, $K_R$, $k_M$ and $K_M$ are based on the literature [30]. For the general case introduced in Hypothesis 3.2.1 we will define

$$k_R \doteq \lim_{p_V \to \infty} w_R(p_V), \quad k_M \doteq \lim_{s_V \to \infty} w_M(s_V).$$

### 3.2.3 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 m \doteq \beta m_V, \quad s \doteq \beta s_V, \quad x \doteq \beta x_V,$$

and the biomass fraction of the bioreactor

$$V \doteq \frac{V_L}{V_{ext}},$$

which, replacing in (3.3), yields

$$m + r = 1. \qquad (3.4)$$

We define the non-dimensional time $\hat{t} \doteq k_R t$ and the non-dimensional functions

$$\hat{w}_R(p) = \frac{w_R(p_V)}{k_R}, \quad \hat{w}_M(s) = \frac{w_M(s_V)}{k_R}$$
so that $\lim_{p \to \infty} \hat{w}_R(p) = 1$. Using (3.4), $m$ can be expressed in terms of $r$, and so the dynamical equation of $m$ can be removed from the system. For the sake of simplicity, let us drop all hats from the current notation. Thus, system becomes

$$ \left\{ \begin{aligned} \dot{s} &= -w_M(s)(1-r)\mathcal{V}, \\ \dot{p} &= w_M(s)(1-r) - \gamma w_R(p)(1-r) - w_R(p)r(p+1), \\ \dot{r} &= (u-r)w_R(p)r, \\ \dot{\mathcal{V}} &= w_R(p)r\mathcal{V}, \\ \dot{x} &= \gamma w_R(p)(1-r)\mathcal{V}. \end{aligned} \right. \qquad (S) $$

## 3.3 Model analysis

**Lemma 3.3.1.** *The set,*

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

is positively invariant for the initial value problem.

Thus, we fix initial conditions

$$ s(0) = s_0 > 0, \quad p(0) = p_0 > 0, \quad r(0) = r_0 \in (0, 1), \quad \mathcal{V}(0) = \mathcal{V}_0 > 0, \quad x(0) = 0. \qquad (\text{IC}) $$

Some relations are immediate from the dynamics: as $\dot{s} \le 0$ and $\dot{\mathcal{V}} \ge 0$ for all $t$, we have

$$ s(t) \le s_0, \quad \mathcal{V}(t) \ge \mathcal{V}_0, \qquad (3.5) $$

representing the fact that the substrate can only be consumed (and not replenished), and the biomass can only grow.

### 3.3.1 Total available mass

As typically occurs in batch processes, there is neither inflow nor outflow of mass in the bioreactor, which is reflected in the dynamics of the system though a mass
conservation law. We define the constant

$$ \Sigma \doteq s_0 + (p_0 + 1)\mathcal{V}_0, $$

representing the initial mass concentration in the system. It can be seen that the total mass concentration

$$ z \doteq s + (p + 1)\mathcal{V} + x $$

is constant for all $t$ (as $\dot{z} = 0$). This means that

$$ s(t) + (p(t) + 1)\mathcal{V}(t) + x(t) = \Sigma, \qquad (\text{First integral 1}) $$

for all $t$. Variables $\mathcal{V}$ and $x$ are maximal when the remaining variables are equal to 0, and so they are upper bounded. In particular, both $\mathcal{V}(t)$ and $x(t)$ are decreasing w.r.t. $s(t)$ and $p(t)$. As neither $s$ nor $p$ can be negative, we have that

$$ \mathcal{V}(t) + x(t) = \Sigma \qquad (3.6) $$

when $s(t) = p(t) = 0$. This condition means that all the available substrate and precursor metabolites have been depleted and transformed into biomass and metabolites, which is intuitively what one would expect from system (S) for $t$ sufficiently large. Additionally, using (3.5) and (First integral 1), we can obtain the following result.

**Proposition 3.3.2.** $\mathcal{V}(t) \in [\mathcal{V}_0, \Sigma]$, $x(t) \in [0, \Sigma - \mathcal{V}_0]$ and $p(t) \in [0, p_{max}]$ for all $t$, with $p_{max} = \Sigma/\mathcal{V}_0 - 1$.

### 3.3.2 Infinite-time full depletion

Dynamics (S) shows that, under initial conditions (IC), $s(t)$ and $p(t)$ can only vanish asymptotically, that is, when $t \to \infty$. The latter can be proved by seeing that the dynamical equations of $s$ and $p$ can be bounded to

$$ \dot{s} \geq -w_M(s)\Sigma, \quad \dot{p} \geq -w_R(p)p_{max}(p_{max} + 1 + \gamma), $$

which means that, at worst, $s$ and $p$ decay exponentially (as functions $w_i(x)$ can be upper bounded by linear functions $w_i(x) \leq c_i x$), and thus $s(t) = p(t) = 0$ cannot be
attained in finite time. Thus, we define the depletion of $s$ and $p$ in an infinite time horizon.

**Definition 3.3.3.** System (S) achieves Full depletion when all the substrate and the precursors are asymptotically depleted, i.e.

$$ \lim_{t \to \infty} s(t) = \lim_{t \to \infty} p(t) = 0, \qquad (\text{Full depletion}) $$

### 3.3.3 Asymptotic behaviour

Now, we will study the system dynamics for an infinite time $t \to \infty$.

#### Constant allocation $u^*$

**Theorem 3.3.1.** For any trajectory of system (S) with initial conditions (IC) and constant allocation $u(t) = u^*$, it follows that

$$ r(t) = u^* - (u^* - r_0) \frac{\mathcal{V}_0}{\mathcal{V}(t)}. \qquad (\text{First integral 2}) $$

*Proof.* Under a constant allocation $u(t) = u^*$, the dynamics of $r$ becomes

$$ \dot{r} = (u^* - r)\mu(p)r. $$

Using dynamics (S), it is possible to see that the quantity $R_u = (u^* - r)\mathcal{V}$ is constant (as $\dot{R}_u = 0$), meaning that $(u^* - r(t))\mathcal{V}(t) = (u^* - r_0)\mathcal{V}_0$ for all $t$, which yields (First integral 2). $\square$

#### General allocation $u(t)$

Due to the boundedness of $\mathcal{V}$ stated in Lemma 3.3.2, and the relation between $\mathcal{V}$ and $r$ shown in (First integral 1), we can see that any constant control $u^*$ yields a bounded ribosomal concentration $r$. We will extend this notion to any function $u(t)$.

**Lemma 3.3.4.** For any trajectory of system (S) with initial conditions (IC) and any control $u(t)$, the ribosomal concentration has bounds $r(t) \in [r_{\min}, r_{\max}]$ for all $t$, with

$$ r_{\min} \doteq r_0 \frac{\mathcal{V}_0}{\Sigma} > 0, \qquad r_{\max} \doteq 1 - (1 - r_0) \frac{\mathcal{V}_0}{\Sigma} < 1. $$
**Proof.** Let us extend system (S) by defining variables $r_{low}(t)$ and $r_{up}(t)$ with dynamics

$$ 
\begin{aligned}
\dot{r}_{\text{low}} &= -r_{\text{low}}w_R(p)r \le 0, & \dot{r}_{\text{up}} &= (1-r_{\text{up}})w_R(p)r \ge 0, \\
r_{\text{low}}(0) &= r_0, & r_{\text{up}}(0) &= r_0,
\end{aligned}
$$

which correspond to the dynamics of $r$ with $u = 0$ and $u = 1$ respectively, and which satisfy

$$ r_{\text{low}}(t) \le r(t) \le r_{\text{up}}(t) $$

for all $t$. The latter can be easily proved by showing that the time-varying differences

$$ \Delta_{\text{low}}(t) = r(t) - r_{\text{low}}(t), \quad \Delta_{\text{up}}(t) = r_{\text{up}}(t) - r(t) $$

with dynamics

$$ \dot{\Delta}_{\text{low}} = (u - \Delta_{\text{low}})w_R(p)r, \quad \dot{\Delta}_{\text{up}} = (1 - u - \Delta_{\text{up}})w_R(p)r $$

are always non-negative: they satisfy $\Delta_{\text{low}}(0) = \Delta_{\text{up}}(0) = 0$ and are repulsive or (at worst) invariant at 0. Then, based on the same principle used to obtain (First integral 2), we define the quantities $R_{\text{low}} = r_{\text{low}}V$ and $R_{\text{up}} = (1 - r_{\text{up}})V$ which are constant (as $\dot{R}_{\text{low}} = \dot{R}_{\text{up}} = 0$), and so

$$ r_{\text{low}}(t) = r_0 \frac{V_0}{V(t)}, \quad r_{\text{up}}(t) = 1 - (1 - r_0) \frac{V_0}{V(t)}, $$

for all $t$. As $V_0 \le V(t) \le \Sigma$ for all $t$, we have

$$ r_{\text{low}}(t) \in \left[ r_0 \frac{V_0}{\Sigma}, r_0 \right], \quad r_{\text{up}}(t) \in \left[ r_0, 1 - (1 - r_0) \frac{V_0}{\Sigma} \right] $$

which shows that $r_{\min} \le r(t) \le r_{\max}$ for all $t$. $\square$

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

Lemma 3.3.4 states that, for any control $u(t)$, the ribosomal concentration never reaches the bounds $r = 0$ and $r = 1$, and thus neither the substrate intake nor the protein synthesis is arrested. We will see that this feature produces (Full depletion) in the following theorem.

**Theorem 3.3.2.** Any trajectory of system (S) with initial conditions (IC) and any control $u(t)$ achieves (Full depletion).
*Proof.* Using Lemma 3.3.4, it is easy to see that

$$ \dot{s} \leq -w_M(s)(1 - r_{\max})\mathcal{V}_0, $$

which means that $s(t)$ converges to 0 as $t \to \infty$. Then, this means that

$$ \dot{p} \leq -\gamma w_R(p)(1 - r_{\max}) - w_R(p)r_{\min}, $$

and so $p(t)$ also converges to 0 as $t \to \infty$. □

## 3.4 The biomass maximization case

For this section, we will write the problem of maximizing the biomass for both infinite time and finite time. The latter is a mathematical representation of the naturally-evolved resource allocation strategy used by bacteria in nature. For that, we will assume that no metabolite is produced, as the pathway responsible for its production is artificially engineered, and thus not present in wild-type bacteria. This is simply modeled through $\gamma = 0$. The resulting system is

$$ \left\{ \begin{array}{l} \dot{s} = -w_M(s)(1-r)\mathcal{V}, \\ \dot{p} = w_M(s)(1-r) - w_R(p)r(p+1), \\ \dot{r} = (u-r)w_R(p)r, \\ \dot{\mathcal{V}} = w_R(p)r\mathcal{V}, \end{array} \right. \quad (S_V) $$

### 3.4.1 Infinite-time problem

#### Problem formulation

We first write the biomass maximization problem for an infinite time horizon. The problem can be expressed as

$$ \max \lim_{u(t) \atop t \to \infty} \mathcal{V}(t). $$
Since $\mathcal{V} \in [\mathcal{V}_0, \Sigma]$, applying (Full depletion) in (First integral 1) yields the condition

$$ \lim_{t \to \infty} \mathcal{V}(t) = \Sigma. $$

meaning that, in infinite time, the biomass is maximized for every control $u(t)$. As a consequence, using Lemma 3.3.1, we have the following result for constant allocations.

**Corollary 3.4.1.** For any trajectory of system $(S_\mathcal{V})$ with initial conditions (IC) and constant control $u(t) = u^*$,

$$ \lim_{t \to \infty} r(t) = u^* - (u^* - r_0) \frac{\mathcal{V}_0}{\Sigma} $$

### Numerical simulations

Examples of trajectories confirming the analytical results are shown in Figure 3.2 and Figure 3.3, where we see that the system approaches (Full depletion) asymptotically in every case, thus approaching the maximal biomass value $\mathcal{V}(t) = \Sigma$. Figure 3.2 shows the resulting trajectories associated to the same initial conditions, when varying the allocation parameter $u$. On the other hand, Figure 3.3 illustrates the trajectories for different values of $r_0$.

## 3.4.2 Finite-time problem

### Problem formulation

Let us fix a final time $t_f > 0$, and write the OCP maximizing the final bacterial volume $\mathcal{V}(t_f)$ with initial conditions (IC). The latter writes

$$ \left\{ \begin{array}{ll} \text{maximize} & \mathcal{V}(t_f), \\ \text{subject to} & \text{dynamics of } (S_\mathcal{V}), \\ & \text{initial conditions (IC)}, \\ & u(\cdot) \in \mathcal{U}. \end{array} \right. \qquad (\text{OCP}_\mathcal{V}) $$

For this class of optimal control problem, where there are no terminal constraints, there is no controllability issues. Additionally, the dynamics is affine in the control,
<figure><img src="image_13.png" /><figcaption>Figure 3.2: Simulation of ($S_v$) with $s_0 = 0.3$, $p_0 = 0.001$, $r_0 = 0.8$, $V_0 = 0.003$ and different allocation functions $u$.</figcaption></figure>

with the latter included in a compact and convex set (a closed interval), and it can be checked that every finite-time trajectory remains bounded. Thus, existence of a solution is guaranteed by Filippov's theorem [53]. Then, for a problem (OCP$_V$) with state $\varphi \in \mathbb{R}^n$, PMP ensures that there exist $\lambda^0 \le 0$ and a piecewise absolutely continuous mapping $\lambda(\cdot) : [0, t_f] \to \mathbb{R}^n$, with $(\lambda(\cdot), \lambda^0) \ne (0, 0)$, such that the extremal $(\varphi, \lambda, \lambda^0, u)$ satisfies the generalized Hamiltonian system

$$ 
\left\{ 
\begin{array}{l}
\dot{\varphi} = \frac{\partial}{\partial \lambda} H(\varphi, \lambda, \lambda^0, u), \\
\\
\dot{\lambda} = -\frac{\partial}{\partial \varphi} H(\varphi, \lambda, \lambda^0, u), \\
\\
H(\varphi, \lambda, \lambda^0, u) = \max_{u \in [0,1]} H(\varphi, \lambda, \lambda^0, u),
\end{array}
\right.
\qquad (\text{PMP})
$$
<figure><img src="image_14.png" /><figcaption>Figure 3.3: Simulation of ($S_V$) with $s_0 = 0.3$, $p_0 = 0.001$, $V_0 = 0.003$, $u = 0.5$ and different values of $r_0$.</figcaption></figure>

for almost every $t \in [0, t_f]$. We define the adjoint states for this particular case as
$\lambda = (\lambda_s, \lambda_p, \lambda_r, \lambda_V)$, and we write the Hamiltonian

$$H = -w_M(s)(1-r)\mathcal{V}\lambda_s + (w_M(s)(1-r) - w_R(p)r(p+1))\lambda_p + w_R(p)r\mathcal{V}\lambda_V + (u-r)w_R(p)r\lambda_r,$$

and the adjoint system as

$$\left\{
\begin{aligned}
\frac{d\lambda_s}{dt} &= w'_M(s)(1-r)(\mathcal{V}\lambda_s - \lambda_p), \\
\frac{d\lambda_p}{dt} &= (w'_R(p)r(p+1) + w_R(p)r)\lambda_p - w'_R(p)r\mathcal{V}\lambda_V - (u-r)w'_R(p)r\lambda_r, \\
\frac{d\lambda_r}{dt} &= -w_M(s)(\mathcal{V}\lambda_s - \lambda_p) + w_R(p)(p+1)\lambda_p - w_R(p)\mathcal{V}\lambda_V - (u-2r)w_R(p)\lambda_r, \\
\frac{d\lambda_V}{dt} &= w_M(s)(1-r)\lambda_s - w_R(p)r\lambda_V.
\end{aligned}
\right.$$
Given that there are no terminal conditions on the state, the transversality conditions for the adjoint state are $\lambda(t_f) = (0, 0, 0, \lambda_0)$, where the presence of $\lambda_0$ is due to the cost function being $\mathcal{V}(t_f)$. Since the Hamiltonian is linear in the control $u$, we rewrite it in the input-affine form $H = H_0 + uH_1$ with

$$ 
\begin{aligned}
H_0 &= -w_M(s)(1-r)\mathcal{V}\lambda_s + (w_M(s)(1-r) - w_R(p)r(p+1))\lambda_p \\
    &\quad + w_R(p)r\mathcal{V}\lambda_v - w_R(p)r^2\lambda_r, \\
H_1 &= w_R(p)r\lambda_r.
\end{aligned}
$$

The constrained optimal control $u$ should maximize the Hamiltonian, so the solution of $(\text{OCP}_v)$ is

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

where $u_{\text{sing}}(t)$ is called a singular control, showing that any optimal control is a concatenation of bangs ($u = 0$ and $u = 1$) and singular arcs, depending on the sign of the switching function $H_1$.

## Singular arcs

The singular arc is produced when $H_1$ vanished (as well as its successive derivatives w.r.t. time) on a whole sub-interval $[t_1, t_2] \subset [0, t_f]$. As $r$ is bounded and $p$ cannot vanish in finite time, along the singular arc one has

$$ \phi_0 = 0, \quad \text{with } \phi_0 = \lambda_r, \qquad (C1) $$

and thus, along a singular arc, the Hamiltonian becomes

$$ H = -w_M(s)(1-r)\mathcal{V}\lambda_s + (w_M(s)(1-r) - w_R(p)r(p+1))\lambda_p + w_R(p)r\mathcal{V}\lambda_v. $$

We differentiate (C1) and we get

$$ \phi_1 = 0, \quad \text{with } \phi_1 = -w_M(s)(\mathcal{V}\lambda_s - \lambda_p) + w_R(p)(p+1)\lambda_p - w_R(p)\mathcal{V}\lambda_v. \qquad (C2) $$
Then, the Hamiltonian can be expressed as

$$H = -w_M(s)(\mathcal{V}\lambda_s - \lambda_p) - r\phi_1 + (u-r)H_1 = c$$

and so $\phi_1$ becomes

$$\phi_1 = c + w_R(p) ((p+1)\lambda_p - \mathcal{V}\lambda_\mathcal{V}).$$

We differentiate (C2) and we get

$$\phi_2 = 0, \quad \text{with } \phi_2 = \frac{w_R^2(p)}{w_M(s)w_R'(p)} - 1. \qquad (C3)$$

Given that $u$ does not appear in (C3), we have the following result.

**Lemma 3.4.2.** Any singular arc is of order two and, along it, $s$ can be expressed in terms of $p$ (provided that $w_M(s)$ is invertible) through the equation $w_M(s) = x(p)$, where $x(p) : \mathbb{R}^+ \to \mathbb{R}^+$ is defined as $x(p) \doteq w_R^2(p)/w_R'(p)$.

We differentiate (C3) and we get

$$\phi_3 = 0, \quad \text{with } \phi_3 = w'_M(s)w_M(s)(1-r)\mathcal{V} + x'(p)\left(w_M(s)(1-r) - w_R(p)r(p+1)\right), \qquad (C4)$$

which shows another relation between states. The latter implies:

**Lemma 3.4.3.** Along the singular arc, $r$ can be expressed in terms of $p$ and $\mathcal{V}$ as $r = \phi(p, \mathcal{V})$ with $\phi(p, \mathcal{V}) : \mathbb{R}^+ \times \mathbb{R}^+ \to (0, 1)$ defined as

$$\phi(p, \mathcal{V}) \doteq \frac{x(p)\left(w'_{M}(s)\mathcal{V} + x'(p)\right)}{x(p)\left(w'_{M}(s)\mathcal{V} + x'(p)\right) + x'(p)w_{R}(p)(p+1)}$$

Now we have:

**Theorem 3.4.1.** The singular optimal control $u(t) = u(p, \mathcal{V})$ (i.e. is in feedback form).

*Proof.* Derivating the expression $r = \phi(p, \mathcal{V})$ and solving for $u$ yields

$$u = \frac{1}{w_R(p)r} \frac{d\phi(p, \mathcal{V})}{dt} + r.$$
As $\phi$ only depends on $p$ and $V$, its derivative will only depend on the state. Then, since $s$ and $r$ can be expressed in terms of $p$ and $V$, the control $u$ becomes a function of these variables. $\square$

## Numerical simulations

The optimal trajectories were computed with Bocop [4], which solves the OCP through a direct method. The time discretization algorithm used is Lobato IIIC (implicit, 4-stage, order 6) with 2000 time steps. Figures 3.4 and 3.5 show optimal trajectories for the same set of initial conditions and different values of $t_f$. Using the mass conservation law (First integral 1), the quantities are represented in the plots as fractions of the total mass in the bioreactor $\Sigma$. The optimal control is characterized by the presence of the chattering phenomenon after and before the singular arc, as expected in singular arcs of order two. From a biological point of view, both allocation strategies prioritize the synthesis of proteins of the metabolic machinery M (red in both Figures): the singular arc takes rather small values, and, depending on the choice of $t_f$, a large proportion of the optimal control at the end of the trajectory is composed of a bang arc $u = 0$. The latter strategy promotes nutrient uptake, which results in a faster depletion of the substrate. By assumption, the concentration of precursor metabolites in the bioreactor remains negligible in comparison with the other quantities. Figure 3.6 illustrates the trajectory represented in Figure 3.4 in the sp-plane, showing that the trajectory approaches the curve $\phi_2(p, s) = 0$ obtained from the singular surface, slides along it, and then follows a trajectory obtained from the $u = 0$ arc that goes towards (Full depletion).

## 3.5 The product maximization case

As done in the previous section, we approach the product maximization objective in infinite time and finite time using the full model (S) where $\gamma \in \mathbb{R}^+$.
<figure><img src="image_15.png" /><figcaption>Figure 3.4: Simulation with $s_0 = 0.1$, $p_0 = 0.003$, $r_0 = 0.1$, $V_0 = 0.003$ and $t_f = 30$. The final volume $V(t_f)$ is at 95% of $\Sigma$ (knowing that $mV + rV = V$).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 3.5: Simulation with $s_0 = 0.1$, $p_0 = 0.003$, $r_0 = 0.1$, $V_0 = 0.003$ and $t_f = 40$. The final volume $V(t_f)$ is at 99.8% of $\Sigma$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 3.6: Same trajectory as Figure 3.4, in the sp-plane. The state approaches the curve $\phi_2 = 0$ and slides along it.</figcaption></figure>

### 3.5.1 Infinite-time problem

The problem of maximizing the product concentration at infinite time is given by the expression

$$ \max \lim_{u^* \to \infty} x(t), $$

which, using (3.6), can be rewritten as

$$ \min \lim_{u^* \to \infty} \mathcal{V}(t) $$

indicating that maximizing the metabolite concentration at infinite time equates to minimizing the biomass. While, in the previous section, the conditions (Full depletion) and (First integral 1) were sufficient to determine the asymptotic behavior of the system, the presence of $x$ in this particular problem does not allow a similar resolution. Thus, we write a relaxed version of the problem in terms of a constant ribosomal concentration $r^*$, which can provide an insight into the optimization problem.
Constant ribosomal concentration

The system with constant ribosomal concentration $r^*$ writes

$$ egin{cases} 
\dot{s} = -w_M(s)(1-r^*)\mathcal{V}, \\ 
\dot{p} = w_M(s)(1-r^*) - \gamma w_R(p)(1-r^*) - w_R(p)r^*(p+1), \quad (S_r) \\ 
\dot{V} = w_R(p)r^*\mathcal{V}. 
\end{cases} $$

It can be seen that the study of the asymptotic behavior of system (S) applies to $(S_r)$ as the latter is a particular case of the original one (S) where $r(0) = u^* = r^*$. We then maximize the final product $x^*$ in terms of the constant ribosomal concentration $r^* \in [r_{min}, r_{max}]$. We can see that the quantity

$$ z = s + (p+1)\mathcal{V} + \gamma \frac{1-r^*}{r^*} \mathcal{V} $$

is constant. Thus,

$$ \mathcal{V}^* + \gamma \frac{1 - r^*}{r^*} (\mathcal{V}^* - \mathcal{V}_0) = \Sigma $$

which, using (3.6), yields

$$ x^* = \gamma \frac{1 - r^*}{r^*} (\mathcal{V}^* - \mathcal{V}_0). $$

Using the fact that $\mathcal{V}^* + x^* = \Sigma$ from (3.6), we see that $x^*$ is monotone decreasing w.r.t. $r^*$, and so the value maximizing the product is $r^* = r_{min}$. This is what one would expect intuitively in an infinite time horizon, as $r^* = r_{min}$ favours the production of M, which catalyzes the production of X without arresting the production of biomass (given by the case $r^* = 0$, which cannot be attained in trajectories starting in $\Gamma$). However, in a finite horizon, a first phase dedicated to bacterial growth could also foster the production of X, which depends directly on the concentration of bacteria in the bioreactor.
### 3.5.2 Finite-time problem

#### Problem formulation

In current section, we study the metabolite production objective, in which the final concentration of metabolite in the bioreactor $x(t_f)$ is maximized. Maximizing the final biomass $V(t_f)$ was already done for model ($S_v$) representing a wild-type bacteria. However, it is likely that the presence of the heterologous pathway responsible for the production of $x$ might affect the results already obtained. Thus, the two objectives will be compared in model (S) from a numerical perspective. Given a fixed final time $t_f > 0$, the OCP maximizing $cV(t_f) + (1-c)x(t_f)$ (with $c = 0$ or $c = 1$ depending on the objective) with initial conditions (IC) writes

$$ \left\{ \begin{array}{l} \text{maximize} \quad cV(t_f) + (1-c)x(t_f), \\ \text{subject to} \quad \text{dynamics of (S),} \\ \text{and} \quad u(\cdot) \in \mathcal{U}, \end{array} \right. \qquad (\text{OCP}_x) $$

One can easily see that, given the dynamics of the system, applying PMP would yield a Hamiltonian linear in the control for both values of $c$, which means that the solution of $(\text{OCP}_x)$ is similar to that of $(\text{OCP}_v)$, given by expression (3.7). A numerical analysis of this results is provided in next section.

#### Numerical simulations

The optimal trajectories were obtained following the same procedure as in the biomass maximization case. Figures 3.7 and 3.8 illustrate optimal trajectories for the same set of initial conditions and same final time $t_f$. Figure 3.7 is the solution of $(\text{OCP}_x)$ where the objective is the product maximization $x(t_f)$. As expected, and similar to the results obtained for $(\text{OCP}_v)$, the optimal control takes the value $u = 0$ for most of the interval, which promotes the synthesis of proteins of the metabolic machinery M, catalyzing the production of $x$. Additionally, according to our simulations, these results do not depend on the final time $t_f$: the final bang $u = 0$ of the optimal control is always predominant in the control strategy. The finite-time numerical results match the results obtained for the infinite-time case, in which the ribosomal sector of the cell should be minimized to maximize the production of $x$. Figure 3.8 shows an optimal trajectory solution of $(\text{OCP}_x)$ with cost function $V(t_f)$.
In this case, the allocation strategy is described by an initial bang $u = 1$ followed by a singular arc that takes up most of the optimal solution, with values near to an intermediate strategy $u = 0.5$; and a small bang $u = 0$ at the end. Such strategy leads to a bacterial composition much more balanced between ribosomal and enzymatic proteins, in opposition to the metabolite production case, where most of the bacterial proteins were dedicated to the metabolic machinery. The latter behavior illustrates a natural trade-off between two opposed strategies: maximizing the number of ribosomes to prioritize the synthesis of macromolecules over the production of $x$ and, at the same time, maximizing the enzymatic activity in order to consume the substrate in the medium as fast as possible.

<figure><img src="image_18.png" /><figcaption>Figure 3.7: Solution of (OCP$_x$) for the metabolite maximization case $x(t_f)$, with $s_0 = 0.1$, $p_0 = 0.001$, $r_0 = 0.1$, $V_0 = 0.003$ and $t_f = 60$. The final product concentration $x(t_f)$ is at 62% of $\Sigma$, while the final volume $V(t_f)$ is only at 23% of $\Sigma$.</figcaption></figure>

## 3.6 Discussion

This chapter described a preliminary study of resource allocation in bacteria for batch processing. A model considering the production of a value-added chemical compound is proposed, and a dynamical study of the system based on mass con-
<figure><img src="image_19.png" /><figcaption>Figure 3.8: Solution of (OCP$_x$) for the biomass maximization case $V(t_f)$, with $s_0 = 0.1$, $p_0 = 0.001$, $r_0 = 0.1$, $V_0 = 0.003$ and $t_f = 60$. The final product concentration $x(t_f)$ is at 27% of $\Sigma$, while the final volume $V(t_f)$ is at 56% of $\Sigma$.</figcaption></figure>

servation laws shows that, under all possible resource allocation strategies, all the
substrate in the medium is consumed. Then, the particular case of a wild-type
bacteria with no metabolite production is analyzed, showing that the optimal al-
location propitious for biomass production should yield a very high $m/r$ ratio in
the cell. Paradoxically, for the metabolite production case, these kind of strategies
would rather maximize the production of $x$, while for maximizing the biomass under
the presence of the heterologous pathway, a more balanced allocation is required.
This work will be readily completed with a more detailed analysis of the impact of
the final time in the allocation decision.

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 4

## Continuous bioreactor

This chapter reproduces [AYT4] and [AYT2], published in the 21th IFAC World Congress and the Mathematical Biosciences and Engineering journal respectively.

### 4.1 Introduction

Microorganisms continuously have to contend with environmental changes in nature, and so they have evolved to accordingly adapt their physiology to cope with this unsteadiness. This is done by reorganizing the gene expression machinery, which is accomplished by dynamically allocating resources to different cellular functions. Microorganisms like bacteria exhibit great genetic variability, which is the main driver of natural selection, a phenomenon that depends on the continuous mutations in living populations. For this reason, they have achieved highly optimized resource management strategies. Such is the case of *E. Coli*: For specific environmental conditions, they seek to maximize their growth rate, which is a naturally evolved characteristic of this bacterium, and well known in the scientific community [14]. In this regard, optimality of bacterial organisms has been a central hypothesis in several research allocation studies [21].

Nevertheless, most of previous works in the literature consider the resource allocation problem in steady-state, without taking into account the changing conditions of natural environments [54]. This motivated a series of works with a dynamical-oriented approach to the problem [30, 8], that aimed to shed light on how bacteria tune their allocation mechanisms when facing changes on the nutrient concentration of the medium. Using Optimal Control theory, they investigate how their inter-
nal pathways can be dynamically readjusted so as to maximize their growth rate, obtaining different strategies that can be related to well known natural regulation mechanisms (such as the ppGpp-mediated sensing of the pool of aminoacids). The approach is based on a widely used modeling technique in systems biology: The so-called coarse-grained self-replicator models, used in bacterial growth representations for their simplicity and their capacity to reproduce observed experimental behaviors [45]. Although this kind of single-cell models are somewhat limited when predicting complex phenomena, they can accurately account for bacterial culture growth laws under the right assumptions (e.g., homogeneity of the culture) [55].

From an industry point of view, a natural question triggered by these studies is: How can we divert the natural allocation strategies of bacteria to improve current biotechnological production schemes? Such is the case of the artificial synthesis of metabolites or proteins of interest. The synthesis of such compounds is highly relevant for its wide range of applications: Production of antitumor agents, insulin, antibiotics, immunosuppressive agents and insecticides, among others [56, 7]. Motivated by the increasing understanding of the biosynthetic properties of certain microorganisms, research on this area can potentially lead to more efficient and sustainable production schemes. This is the matter addressed in recent work using a strain of *Escherichia coli* that includes an artificially engineered heterologous pathway for the production of a certain metabolite of interest [3, 74, AYT5]. In this approach, a control loop is developed through a bacterial growth switch that allows to externally modify the natural resource allocation decision [44]. The mechanism is implemented by re-engineering the transcriptional control of the expression of RNA polymerase, a key component of the gene expression machinery. This way, it is possible to optimize the productivity of the bioprocess by channeling resources into this new heterologous pathway.

At the same time, the synthesis of these metabolites draws resources from the native pathways used for producing biomass in bacteria, thus leading to an inherent compromise between these two objectives. One possible approach to this trade-off is to model it through different cost functions, thus obtaining multi-objective optimization problems. This is the case of [41], where the authors aim to maximize the production of a metabolite while minimizing the genetic burden caused by pathway expression. In contrast to this method, the work of [3] models the main trade-offs behind the process through a single decision parameter, which considerably reduces
the complexity of the optimization problem.

In a similar vein, we address a classic production scheme: the Continuous Stirred-Tank Reactor (CSTR), also known as chemostat [57]. While resource allocation in bacteria has been vastly studied in constant environments (e.g., under the assumption that there is always enough substrate in the medium), how this goes in continuous bioreactors is not trivial, since a feedback occurs from the physiology of the cell to the environmental conditions given by the interaction bacteria-medium [26]. Examples of self-replicator models in continuous bioreactors can be found in some recent works: In [58], authors use a coarse-grained self-replicator kinetic model of *Saccharomyces cerevisiae's* in a continuous bioreactor to investigate the trade-off between respiratory and fermentative metabolism, showing that optimal strategies are ‘pure’ metabolic strategies (e.g., either respiration of fermentation, but not respiro-fermentation). Likewise, it is rather classical in the continuous bioreactor scheme to maximize a certain performance measure (i.e., biomass production) in terms of the operational parameters related to the setup, such as dilution rate and/or concentration of the substrate inflow [59, 60, 61, 62, 63]. An example can be seen in [64], where the infinite-time optimal control problem of maximizing the average biogas production in terms of the dilution rate is studied. However, incorporating the aforementioned external control—that can disrupt the natural allocation process of the whole culture—provides an extra degree of freedom, which can, in turn, contribute to further improve the classical production scheme.

Based on the presented works, we show novel results addressing the problem of bacterial resource allocation in the CSTR framework. Our approach is based on a coarse-grained self-replicator dynamical model that accounts for the microbial culture growth inside the continuous bioreactor, and incorporates the external allocation control previously described. The novelty of the approach lies in the combination of the resource allocation control scheme, and the capacity to regulate the bacterial growth rate through the dilution rate of the continuous bioreactor. Further on, we study its asymptotic behavior using dynamical systems theory, and we provide conditions for the persistence of the bacterial population. The analysis is carried out by studying the local and global behavior of its limiting system, and relating its convergence to the original model through the theory of asymptotically autonomous systems [65]. Then, we pose the problem of maximizing the synthesis of the metabolite of interest during a fixed interval of time in terms of the resource
allocation decision. The latter is expressed as an OCP (Optimal Control Problem) which is then explored using PMP (Pontryagin's Maximum Principle) [31]. We analyze the solution of the problem and propose a sub-optimal control strategy given by a constant allocation decision, which eventually takes the system to the optimal steady-state production regime. On this basis, we study and compare the two most significant steady-state production objectives in CSTRs: Biomass production and metabolite production. For this last purpose, and in addition to the allocation parameter, we control the constant volumetric flow of the bioreactor (or dilution rate), and we analyze the results through a numerical approach. The resulting two-dimensional optimization problem is defined in terms of Michaelis-Menten kinetics with the parameter values of [30], and taking into account the constraints for the existence of the equilibrium of interest.

## 4.2 Model definition

### 4.2.1 Self-replicator model

In this section, we define the coarse-grained self-replicator model including the continuous bioreactor scheme and the allocation parameter. We consider a growing bacterial population in a CSTR bioreactor of constant volume $V_{ext}$ [L]. The self-replicating system that models the culture is composed of the metabolic machinery $M$ (transporters, enzymes...) and the gene expression machinery $R$ (RNA polymerase, ribosomes...), both responsible of the cell growth. The validity of this single-cell model as a representation of a growing bacterial culture depends on a number of simplifying assumptions, one of them being that individual cells share the same macromolecular composition. For the production of the metabolite of interest $X$, we consider the artificially engineered pathway introduced in [3] (Figure 4.1).

The model represents three chemical macroreactions,

$$S \xrightarrow{V_M} P$$

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

$$P \xrightarrow{V_X} X$$
<figure><img src="image_20.png" /><figcaption>Figure 4.1: Extended coarse-grained self replicator model introduced in [3]. The substrate in the bioreactor ($S_1$) is consumed by the bacterial culture and transformed into precursors $P$ through the action of the metabolic machinery $M$. Then, precursors are used to make macromolecules of the gene expression machinery $R$ and the metabolic machinery $M$, with proportions $\alpha$ and $1 - \alpha$ respectively; and to synthesize the metabolite $X$, which is excreted from the cell to the bioreactor. The external control $I$ affects how the precursors $P$ are distributed between both cellular functions $M$ and $R$.</figcaption></figure>

The first reaction is catalyzed by $M$ and describes the transformation of external substrate $S$ into precursor metabolites $P$ at a rate $V_M$. The second one represents the conversion of precursors into macromolecules $R$ and $M$ and is catalyzed by $R$, at a rate $V_R$. Finally, the third reaction describes the transformation of precursors $P$ into product $X$ at a rate $V_X$, and catalyzed by $M$. The natural resource allocation parameter is modeled through the dimensionless parameter $\alpha \in [0, 1]$, that represents the proportion of precursors allocated to the gene expression machinery $R$, while $1 - \alpha$ indicates that of the metabolic machinery $M$. The model describes all mass quantities $S$, $P$, $M$, $R$ and $X$ in grams, and the rates in grams per hour. A constant volumetric flow rate $F$ [L h$^{-1}$] produces an inflow of fresh medium rich in substrate, and an outflow of bacterial culture and metabolites [46]. As already stated, we include in our scheme the growth switch described in [44] to externally affect the allocation decision $\alpha$ by varying the inducer concentration in the medium. The mechanism is modeled as

$$u(t) = I(t) \alpha(t), \quad u \in [0, 1] \qquad (4.1)$$
where I is the external signal, which acts independently from the allocation pa- rameter α . The aggregated form of equation (4.1) accounts for the fact that the synthetic switch is capable of adjusting bacterial growth between zero (by setting I = 0, which yields u = 0) and the maximal growth rate supported by the medium (given by I = 1, which yields u = α ). In this work, we limit the analysis to calculat- ing the control input u , without decoupling both individual controls I and α (refer to Section 5 in [3] for more details on the implementation of this external controller). Then, we write the system of differential equations describing the evolution of the mass of each component,

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

The inflow/outflow rates are defined as

$$
V_{S_{out}} = DS, \quad V_{P_{out}} = DP, \quad V_{M_{out}} = DM, \quad V_{R_{out}} = DR, \quad V_{X_{out}} = DX, \quad V_{S_{in}} = F s_{in},
$$

where s$_{in}$ [g L − $^{1}$] is the nutrient concentration of the inflow of fresh medium, and D [h − $^{1}$] the dilution rate defined as

$$
D \doteq \frac{F}{V_{\text{ext}}}.
$$

We define the volume of the cell population V [L] as

$$
\nu \doteq \beta(M + R), \tag{4.2}
$$

where β [L g − $^{1}$] is the inverse of the cytoplasmic density. The above definition is based on the assumption that the cytoplasmic density of the cells is constant for the whole culture, and it also takes into account the experimental results showing that macromolecules are responsible for most of the biomass in microbial cells [6]. Thus, the mass of precursors P is excluded from the volume V as a simplifying assumption.
Then, we express the quantities of the system as concentrations with respect to the volumes,

$$p \doteq \frac{P}{\mathcal{V}}, \quad r \doteq \frac{R}{\mathcal{V}}, \quad m \doteq \frac{M}{\mathcal{V}}, \quad s \doteq \frac{S}{\mathcal{V}_{\text{ext}}}, \quad x \doteq \frac{X}{\mathcal{V}_{\text{ext}}}, \qquad (4.3)$$

where $p, r$ and $m$ [g L⁻¹] are intracellular concentrations of precursor metabolites, components of the gene expression machinery and metabolic enzymes respectively; and $s$ and $x$ [g L⁻¹] the extracellular concentrations of substrate and metabolite. It is worth stressing that intracellular concentrations are defined with respect to the bacterial volume $\mathcal{V}$, while extracellular concentrations with respect to the volume of the bioreactor $\mathcal{V}_{\text{ext}}$. Then, using definition (4.2), we obtain that $m + r = 1/\beta$. We define the growth rate of the bacterial population $\mu$ [h⁻¹] as the relative variation of cell volume $\dot{\mathcal{V}}/\mathcal{V}$ without considering the effect of the outflowing biomass. Replacing with concentrations leads to the system

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

where $v_M(s, m)$, $v_R(p, r)$ and $v_X(p, m)$ [g L⁻¹ h⁻¹] are the mass fluxes per unit volume obtained from dividing the rates $V_M$, $V_R$ and $V_X$ by $\mathcal{V}$, function of the concentrations of system (S₁). In this new system, the growth rate becomes

$$\mu(p, r) \doteq \left. \frac{\dot{\mathcal{V}}}{\mathcal{V}} \right|_{F=0} = \left. \frac{\dot{M} + \dot{R}}{M + R} \right|_{F=0} = \beta v_R(p, r)$$

showing that the bacterial growth rate is proportional to the macromolecule synthesis rate. We propose a change of variables that simplifies the expressions of the
system

$$ \hat{s} = \beta s, \quad \hat{p} = \beta p, \quad \hat{r} = \beta r, \quad \hat{m} = \beta m, \quad \hat{x} = \beta x, \quad \hat{V} = \frac{V}{V_{\text{ext}}}, $$

which yields the relation

$$ \hat{r} + \hat{m} = 1, \qquad (4.4) $$

and we define the non-dimensional synthesis rates

$$ \hat{v}_M(\hat{s}, \hat{m}) = \beta v_M(s, m), \quad \hat{v}_R(\hat{p}, \hat{r}) = \beta v_R(p, r), \quad \hat{v}_X(\hat{p}, \hat{m}) = \beta v_X(p, m), \qquad (4.5) $$

and the non-dimensional substrate concentration $\hat{s}_{in} = \beta s_{in}$. Then, dropping all hats yields the following system

$$ \left\{ \begin{aligned} \dot{s} &= D(s_{in} - s) - v_M(s, 1-r)V, \\ \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{x} &= v_X(p, 1-r)V - Dx, \\ \dot{V} &= (\mu(p, r) - D)V, \end{aligned} \right. \qquad (S_1) $$

where the dynamical expression of $m$ has been omitted since it can be computed from $r$, as shown in equation (4.4). It can also be seen that both concentrations $m$ and $r$ are limited to the interval $[0, 1]$ due to physical constraints from the relation in equation (4.4). For this latter, and due to the nature of the reactions involved in the studied problem, we will make some assumptions on the synthesis rates $v_M(s, m)$, $v_R(p, r)$ and $v_X(p, m)$:

**Assumption 1.** Functions $v_M(s, m)$, $v_R(p, r)$ and $v_X(p, m)$ meet

* $v_i(y, z) : \mathbb{R}^+ \times [0, 1] \rightarrow \mathbb{R}^+$

* $v_i(y, z)$ continuously differentiable w.r.t. both variables
* $v_i(0, z) = v_i(y, 0) = 0$

* $v_i(\cdot)$ strictly monotonically increasing:

$$ \frac{\partial v_i}{\partial y}(y, z) > 0, \forall (y, z) \in \mathbb{R}_{>0} \times (0, 1], \quad \frac{\partial v_i}{\partial z}(y, z) > 0, \forall (y, z) \in \mathbb{R}_{>0} \times (0, 1] $$

* $v_i(\cdot)$ bounded w.r.t. $y$: $\lim_{y \to \infty} v_i(y, z) = v_{i,max}(z)$.

Assumption 1 encompasses all general monotone increasing kinetics models used in biological models, such as Michaelis-Menten or Hill equation-based kinetics [66]. In this first work, we will focus on a particular kind of systems where the synthesis rate related to the metabolite production depends on the growth rate:

**Assumption 2.** For $r \in (0, 1)$, the metabolite synthesis rate $v_X(p, 1-r)$ can be expressed in terms of the macromolecule synthesis rate $v_R(p, r)$ (i.e., the growth rate)

$$ v_X(p, 1-r) = c(r) v_R(p, r), $$

where $c(r) : (0, 1) \to \mathbb{R}^+$ is a positive continuously differentiable function.

As previously described, the reaction $v_X(p, m)$ is catalyzed by $m$, and the reaction $v_R(p, r)$ is catalyzed by $r$, meaning that the ratio between $M$ and $R$ in the microbial culture determines whether the resources are being allocated to the production of biomass or metabolite. This represents the trade-off described in the introduction of this paper, which is here modeled through the function $c(r)$. The fact that $c$ does not depend on the concentration of precursors implies that the host cell has the same affinity to synthesize both biomass and metabolite from the precursors, even when the reactions are not expected to consume the precursors in the same proportion. For the particular case of Michaelis-Menten kinetics, this phenomenon is captured by the half-saturation constant [67]. Notably, for fixed values of $r$, both
synthesis rates are simply proportional. This assumption reduces (S₁) to

$$ 
\left\{ 
\begin{aligned}
\dot{s} &= D(s_{in} - s) - v_M(s, 1 - r)\mathcal{V}, \\
\dot{p} &= v_M(s, 1 - r) - \mu(p, r)(p + c(r) + 1), \\
\dot{r} &= (u - r)\mu(p, r), \\
\dot{x} &= c(r)\mu(p, r)\mathcal{V} - Dx, \\
\dot{\mathcal{V}} &= (\mu(p, r) - D)\mathcal{V}.
\end{aligned}
\right. 
\quad (S_1) 
$$

We note that both Assumptions 1 and 2 are formulated for the non-dimensional synthesis rates given in definition (4.5), but they also hold for the original functions $v_M$, $v_R$ and $v_X$ (taking into account that, for these functions, the domain is defined as $\mathbb{R}^+ \times [0, \frac{1}{\beta}] \to \mathbb{R}^+$).

### 4.2.2 Asymptotic behavior

The asymptotic behavior of system (S₁) describes the “open-loop” operation mode of the continuous bioreactor, where the resource allocation control $u(t)$ is fixed to $\bar{u} \in (0, 1)$. In the present section we propose a series of mass conservation laws that allow to reduce (S₁) to a 3-dimensional limiting system. Then, we study the local stability of their equilibria, and we show, using the theory of asymptotically autonomous systems, that the full system (S₁) converges to the equilibria of its limiting system. Let us start the analysis of system (S₁) by defining its invariant region in the following lemma.

**Lemma 1.** *The set*

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

is positively invariant for the initial value problem.
*Proof.* Let us analyze the boundaries of Γ,

$$
\begin{align*}
\dot{s}|_{s=0} &= Ds_{in} > 0 & \Rightarrow & s = 0 \text{ is repulsive.} \\
\dot{s}|_{s=s_{in}} &= -v_M(s, 1-r)\mathcal{V} \le 0 & \Rightarrow & s = s_{in} \text{ is repulsive/invariant.} \\
\dot{p}|_{p=0} &= v_M(s, 1-r) \ge 0 & \Rightarrow & p = 0 \text{ is repulsive/invariant.} \\
\dot{r}|_{r=0} &= 0 & \Rightarrow & r = 0 \text{ is invariant.} \\
\dot{r}|_{r=1} &= (u-1)\mu(p, 1) < 0 & \Rightarrow & r = 1 \text{ is repulsive.} \\
\dot{x}|_{x=0} &= v_X(p, m)\mathcal{V} \ge 0 & \Rightarrow & x = 0 \text{ is repulsive/invariant.} \\
\dot{V}|_{\mathcal{V}=0} &= 0 & \Rightarrow & \mathcal{V} = 0 \text{ is invariant.}
\end{align*}
$$

□

We will study the initial value problem of system (S₁) with initial conditions

$$s_{in} \geq s(0) \geq 0, \quad p(0) \geq 0, \quad 1 \geq r(0) > 0, \quad x(0) \geq 0, \quad \mathcal{V}(0) > 0, \quad (4.6)$$

where two cases have been excluded for being trivial to the analysis: $\mathcal{V}(0) = 0$, since it is necessary to have an initial amount of biomass to have bacterial growth; and $r(0) = 0$, since an empty gene expression machinery pool implies null growth rate $\mu(p, 0) = 0$, and therefore it is not possible to self-replicate from that point.

## Mass conservation

System (S₁) can be rewritten as

$$\begin{cases}
\dot{\varphi} = D(s_{in}\vec{v}_{in} - \vec{v}_{out}) + N\vec{v}_{i} - \vec{v}_{\mu}\mu(p, r), \\
\dot{V} = (\mu(p, r) - D)\mathcal{V},
\end{cases}$$

where

* $\varphi \doteq [s, p, r, m, x]^T$ is the state vector of concentrations in the bioreactor.

* N is the stoichiometry matrix of the macroreactions.

* $v_i$ is the vector of internal synthesis rates.

* $v_{in}$ and $v_{out}$ are the vectors of inflows and outflows respectively, associated to the continuous bioreactor setup.
• $v_\mu$ is the vector modeling the dilution effect due to variation of the bacterial volume.

which are defined as

$$ N \doteq \begin{bmatrix} -\mathcal{V} & 0 & 0 \\ 1 & -1 & -1 \\ 0 & u & 0 \\ 0 & 1-u & 0 \\ 0 & 0 & \mathcal{V} \end{bmatrix}, \quad \vec{v}_i \doteq \begin{bmatrix} v_M(s, 1-r) \\ v_R(p, r) \\ v_X(p, 1-r) \end{bmatrix}, \quad \begin{array}{l} \vec{v}_{in} \doteq [1, 0, 0, 0]^T, \\ \vec{v}_{out} \doteq \operatorname{diag}(\varphi) [1, 0, 0, 0, 1]^T, \\ \vec{v}_{\mu} \doteq \operatorname{diag}(\varphi) [0, 1, 1, 1, 0]^T. \end{array} $$

By analyzing the left null space of $N$, it can be seen that there are two mass conservation laws related to the total mass inside the bioreactor.

**Definition 1.** We define the quantities

$$ \begin{aligned} w_1 &\doteq s + (p + m + r) \mathcal{V} + x = s + (p + 1) \mathcal{V} + x, \\ w_2 &\doteq s + \left(p + \frac{r}{\bar{u}}\right) \mathcal{V} + x. \end{aligned} $$

The first quantity tends asymptotically to $w_1 = s_{in}$ as $t \to \infty$ for every input $u(t)$, as it obeys the dynamical equation

$$ \dot{w}_1 = D(s_{in} - w_1). \qquad (4.7) $$

Moreover, when fixing $u(t)$ to $\bar{u} \in (0, 1)$, the quantity $w_2$ also obeys the same equation (4.7), meaning that this second quantity also converges to $s_{in}$, which greatly simplifies the analysis of the asymptotic behavior of the system.

**Lemma 2.** The $\omega$-limit set of any solution of system (S$_1$) lies in the hyperplane

$$ \Omega_1 \doteq \{(s, p, r, x, \mathcal{V}) \in \mathbb{R}^5 : s + (p+1)\mathcal{V} + x = s_{in}\}. $$

Moreover, under constant input $u(t) = \bar{u}$, this is also true for the hyperplane

$$ \Omega_2 \doteq \{(s, p, r, x, \mathcal{V}) \in \mathbb{R}^5 : s + \left(p + \frac{r}{\bar{u}}\right) \mathcal{V} + x = s_{in}\}. $$
Further on, we will use Lemma 2 to analyze system (S₁) through its limiting system.

## Limiting systems

Lemma 2 presents two mass conservation laws that can be used to reduce subsystem (S₁) by two dimensions. We will first analyze the asymptotic behavior of concentration r: when $t \to \infty$, quantities $w_1 = w_2 = s_{in}$, so

$$s + (p + 1)\mathcal{V} + x = s + \left(p + \frac{r}{\bar{u}}\right)\mathcal{V} + x \quad \Rightarrow \quad r = \bar{u}$$

meaning that, as $t \to \infty$, r will converge to the value $\bar{u}$. We can also express $x = s_{in} - s - (p + 1)\mathcal{V}$, so that the limiting system of (S₁) becomes

$$\left\{
\begin{aligned}
\dot{s} &= D(s_{in} - s) - \bar{v}_M(s)\mathcal{V}, \\
\dot{p} &= \bar{v}_M(s) - \bar{\mu}(p)(p + \bar{c} + 1), \\
\dot{\mathcal{V}} &= (\bar{\mu}(p) - D)\mathcal{V},
\end{aligned}
\right.
\quad (S'_1)$$

where

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

Details on the convergence of the limiting system (S'₁) to the original one (S₁) will be addressed later in the article. In next section, we will fully describe the asymptotic behavior of (S'₁).

## Local stability

In the interest of simplifying the notation, we define the following function.

**Definition 2.** We define the function

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

Function f meets $\bar{f}(p) > 0, \bar{f}'(p) > 0, \forall p \in \Gamma$ (positive and monotonically in-
creasing).

The main result of local stability study is summarized in Theorem 1.

**Theorem 1.** The local stability of equilibria is given by the following criterion.

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

In the above criterion, the equilibria are defined as follows:

* The interior equilibrium $E_i \doteq (s_i, p_i, \mathcal{V}_i)$, with

$$p_i : \{p \in \mathbb{R}^+ : \bar{\mu}(p) = D\}, \qquad (4.8)$$

$$s_i : \{s \in \mathbb{R}^+ : \bar{v}_M(s) = \bar{f}(p_i)\}, \qquad (4.9)$$

$$\mathcal{V}_i \doteq \frac{D(s_{in} - s_i)}{\bar{v}_M(s_i)}. \qquad (4.10)$$

* The washout equilibrium $E_w \doteq (s_{in}, p_w, 0)$, with

$$p_w : \{p \in \mathbb{R}^+ : \bar{f}(p) = \bar{v}_M(s_{in})\}. \qquad (4.11)$$

*Proof.* First, we will prove the boundedness of $p$. Using definition (4.11), it is possible to define a constant upper-bound on $p$ by analyzing the time-varying upper bound $p_{up}(t)$ with dynamical equation

$$\dot{p}_{up} \doteq \bar{v}_M(s_{in}) - \bar{f}(p_{up}). \qquad (4.12)$$

In equation (4.12), $p_{up}$ converges to the value $p_w$ satisfying equation (4.11). Additionally, the vector field at $p = p_w$ is always negative (or null when $s = s_{in}$)

$$\dot{p}|_{p=p_w} = \bar{v}_M(s) - \bar{f}(p_w) \le 0$$
meaning that $p = p_w$ is either repulsive or invariant, so a new invariant set $\Gamma' \subset \Gamma$ can be defined,

$$ \Gamma' \doteq \{(s, p, \mathcal{V}) \in \mathbb{R}^3 : s_{in} \ge s > 0, p_w \ge p > 0, \mathcal{V} \ge 0\}. \quad (4.13) $$

**Equilibrium E<sub>i</sub>** The existence of the interior equilibrium E<sub>i</sub> is given by the boundedness of the flows,

$$ \max \bar{\mu}(p) \ge D, \quad (4.14) $$

$$ \max \bar{v}_M(s) \ge \bar{f}(p_i), \quad (4.15) $$

and its uniqueness can be proved through monotonicity arguments: In $\bar{\mu}(p) = D$, $\bar{\mu}(p)$ is strictly monotonically increasing w.r.t. $p$ so, if inequality (4.14) is met, $p_i$ should be unique. Similarly, in $\bar{v}_M(s) = \bar{f}(p_i)$, $\bar{v}_M(s)$ is again strictly monotonically increasing so, if inequality (4.15) is met, $s_i$ should be unique. As is standard in continuous bioreactors, inequality (4.14) implies that the maximal growth rate of the bacterial population should be bigger than the dilution rate $D$. At the same time, the inequality (4.15) requires the maximal uptake flow to be bigger than the flows responsible for bacterial growth and metabolite production. In $\Gamma'$, these conditions become

$$ \bar{\mu}(p_w) \ge D, \quad (4.16) $$

where the second condition (4.15) is included in the inequality (4.16), as $\max \bar{v}_M(s) = \bar{v}_M(s_{in}) = \bar{f}(p_w)$, and the condition $\bar{f}(p_w) \ge \bar{f}(p_i)$ is true if and only if $\mu(p_w) \ge D$. The Jacobian matrix is given by

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

and the characteristic polynomial is

$$ \begin{aligned} P_i(\lambda) &= (\lambda + D + \bar{v}'_M(s)\mathcal{V}_i)(\lambda + \bar{f}'(p))\lambda + \bar{v}_M(s)\bar{v}'_M(s)\bar{\mu}'(p)\mathcal{V}_i \quad (4.17) \\ &= \lambda^3 + \lambda^2(D + \bar{v}'_M(s)\mathcal{V}_i + \bar{f}'(p)) + \lambda(D + \bar{v}'_M(s)\mathcal{V}_i)\bar{f}'(p) + \bar{v}_M(s)\bar{v}'_M(s)\bar{\mu}'(p)\mathcal{V}_i. \end{aligned} $$
Using Definition 2, it can be seen that $-D$ is an eigenvalue by replacing in expression 4.17

$$P_i(-D) = -D\bar{v}_M'(s)\mathcal{V}_i(p + \bar{c} + 1)\bar{\mu}'(p) + \bar{v}_M(s)\bar{v}_M'(s)\bar{\mu}'(p)\mathcal{V}_i = \bar{v}_M'(s)\mathcal{V}_i\bar{\mu}'(p)\dot{p} = 0.$$

Then, dividing $P_i(\lambda)$ by $\lambda + D$ yields

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

which, by Routh-Hurwitz criteria, implies that all eigenvalues have negative real part. Then, if it exists, the equilibrium is always stable.

**Equilibrium $E_w$** The washout equilibrium $E_w$ exists for all values of $\Gamma'$, since the only condition for existence is given by

$$\max \bar{f}(p) \geq \bar{v}_M(s_{in}),$$

which is always true in $\Gamma'$ since $\bar{v}_M(s_{in}) = \bar{f}(p_w)$ and so the inequality becomes $\bar{f}(p_w) \geq \bar{f}(p_w)$. Again, as $\bar{f}(p)$ is strictly monotonically increasing, there is a unique solution $p_w$. Its stability is given by

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

Eigenvalues $\lambda = -D$ and $\lambda = -\bar{f}'(p)$ are always real and negative, and so the stability criterion becomes

$E_w$ stable if and only if $\bar{\mu}(p_w) < D$.

□
## Global analysis

In the current section, we show how the limiting system ($S'_1$) can be further reduced using another mass conservation law given by the fact that $\bar{c}$ is constant. This is formalized in the following lemma.

**Lemma 3.** The $\omega$-limit set of any solution of the limiting system ($S'_1$) lies in the hyperplane

$$ \Omega_3 \doteq \{(s, p, \mathcal{V}) \in \mathbb{R}^3 : s + (p + \bar{c} + 1) \mathcal{V} = s_{in}\} $$

*Proof.* The quantity

$$ w_3 \doteq s + (p + \bar{c} + 1) \mathcal{V} \quad (4.18) $$

obeys the dynamical equations $\dot{w}_3 = D(s_{in} - w_3)$, which means that the system converges asymptotically to the limit set $\Omega_3$. $\square$

Using Lemma 2, we can further reduce system ($S'_1$) by expressing

$$ s = s_{in} - \mathcal{V}(p + \bar{c} + 1), \quad (4.19) $$

for the values $(p, \mathcal{V})$ that meet $s_{in} - \mathcal{V}(p + \bar{c} + 1) > 0$, and so the new limiting system becomes

$$ \left\{ \begin{array}{l} \dot{p} = \bar{v}_M(s(\cdot)) - \bar{\mu}(p)(p + \bar{c} + 1), \\ \dot{\mathcal{V}} = (\bar{\mu}(p) - D)\mathcal{V}. \end{array} \right. \quad (S_1'') $$

Since ($S_1''$) is a 2-dimensional continuous system, its global behavior (illustrated in Figure 4.2) can be studied through the Poincaré-Bendixson trichotomy.

**Lemma 4.** Every solution of the limiting system ($S_1''$) with initial conditions (4.6) converges to

$$ \bullet E_i^* \doteq (p_i, \mathcal{V}_i) \text{ if } \bar{\mu}(p_w) \ge D $$

$$ \bullet E_w^* \doteq (p_w, 0) \text{ if } \bar{\mu}(p_w) < D $$
**Proof.** The Poincaré-Bendixson trichotomy ensures that every non-empty compact $\omega$-limit set of $(S_1'')$ is either a fixed point, a periodic orbit or a cycle of equilibria. Then, by applying Poincaré-Bendixson theorem through the Dulac criterion, we can discard periodic orbits and cycles of equilibria:

$$ 
\begin{aligned}
\frac{\partial}{\partial p}\dot{p} + \frac{\partial}{\partial \mathcal{V}}\dot{\mathcal{V}} &= \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,
\end{aligned}
\quad (4.20) $$

as $\partial s(\cdot)/\partial p = -\mathcal{V}$ from equation (4.19). This ensures that the new limiting system should converge to one of the stable equilibria as $t \to \infty$ which, according to Theorem 1, is known to be unique. $\square$

<figure><img src="image_21.png" /><figcaption>Figure 4.2: Phase plane of the limiting system $S_1''$ showing: The case where $\bar{\mu}(p_w) \ge D$ (left), so that the equilibrium $E_i^*$ exists and attracts all solutions; and the case where $\bar{\mu}(p_w) < D$ (right) and $E_w^*$ attracts all solutions.</figcaption></figure>

Through the theory of asymptotically autonomous systems [65], we can relate the asymptotic behavior of the 2-dimensional limiting system $S_1''$ to that of the full 5-dimensional system $(S_1)$. The main result is established in Theorem 2.

**Theorem 2.** Every solution of system $(S_1)$ with initial conditions (4.6) converges to,

* The extended interior equilibrium $\hat{E}_i \doteq (s_i, p_i, \bar{u}, x_i, \mathcal{V}_i)$ if it exists, with $x_i \doteq \bar{c}\mathcal{V}_i$.
• The extended washout equilibrium $\hat{E}_w \doteq (s_{in}, p_w, \bar{u}, 0, 0)$ if $\hat{E}_i$ does not exist. where the condition for the existence of the extended interior equilibrium is $\bar{\mu}(p_w) \ge D$.

*Proof.* We resort to a more particular case of the general theory of asymptotically autonomous systems (introduced in Appendix F of [46]), that requires a certain number of hypotheses to be met related to the limiting system $S_1''$, its equilibria and the mass conservation equations defined in Lemma 2:

(A1) The dynamical equations (4.7) and (4.18) of quantities $w_1, w_2$ and $w_3$ are stable.

(A2) $S_1''$ has 2 rest points $E_i^*$ and $E_w^*$, which are hyperbolic.

(A3) If $E_i^*$ exists, $dim(M^+(E_i^*)) = 2$ and $dim(M^+(E_w^*)) = 1$. If $E_i^*$ does not exist, $dim(M^+(E_w^*)) = 2$.

(A4) $S_1''$ has no cycles of rest points, as shown in equation (4.20).

(A5) $S_1''$ has no periodic orbits, as shown in equation (4.20).

Then, almost all trajectories of the original system ($S_1$) converge to one of the asymptotically stable equilibria of the limiting system, which is always unique. $\square$

## 4.3 Metabolite production

In this section, we aim to maximize the production of the metabolite of interest $X$. We start by modeling the synthesis rates as explicit functions of the concentrations of the system. While this latter is not required for most of the results, it enables the numerical simulations both for the dynamical and for the static optimization problems. In particular, we resort to the Michaelis-Menten kinetics defined in previous works [30] and the same biological constants. Then, we pose the problem of maximizing the production of $X$ over a fixed period of time, and we approximate the optimal solution with a constant allocation strategy. This finally motivates the static optimization approach, that is solved in terms of the tuple $(D, \bar{u})$.
### 4.3.1 Kinetic's definition

Michaelis-Menten kinetics are defined as

$$ 
\begin{aligned}
v_M(s, m) &\doteq k_M m \frac{s}{K_M + s} = k_M (1 - r) \frac{s}{K_M + s}, \\
v_R(p, r) &\doteq k_R r \frac{p}{K_R + p}, \\
v_X(p, m) &\doteq k_X m \frac{p}{K_X + p} = k_X (1 - r) \frac{p}{K_X + p},
\end{aligned}
$$

where the biological constants and their values [30, 3] are described in Table 4.1.

<table><caption>Table 4.1: Relevant biological system variables. RT and HSC stand for "Rate constant" and "Half-saturation constant" respectively.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>β</td><td>Inverse of the cytoplasmic density</td><td>L g<sup>-1</sup></td><td>0.003</td></tr><tr><td>k<sub>M</sub></td><td>RT of the metabolic macroreaction</td><td>h<sup>-1</sup></td><td>4.32</td></tr><tr><td>K<sub>M</sub></td><td>HSC of the metabolic macroreaction</td><td>g L<sup>-1</sup></td><td>33.33</td></tr><tr><td>k<sub>R</sub></td><td>RT of the protein synthesis reaction</td><td>h<sup>-1</sup></td><td>3.6</td></tr><tr><td>K<sub>R</sub></td><td>HSC of the protein synthesis reaction</td><td>g L<sup>-1</sup></td><td>1</td></tr><tr><td>k<sub>X</sub></td><td>RT of the metabolite synthesis reaction</td><td>h<sup>-1</sup></td><td>0.5</td></tr><tr><td>K<sub>X</sub></td><td>HSC of the metabolite synthesis reaction</td><td>g L<sup>-1</sup></td><td>1</td></tr></tbody></table>

Assumption 2 is implemented by setting the half-saturation constant $K_X = K_R$, such that the function $c(r)$ becomes

$$ c(r) = \frac{v_X(p, 1 - r)}{v_R(p, r)} = \frac{k_X (1 - r)}{k_R r} $$

The substrate concentration of the inflow $s_{in}$ is set to 0.4 g/L.

### 4.3.2 Dynamic optimization problem

The problem of maximizing the production of the metabolite X during a fixed interval of time T is explored in this section. As it is classical in the continuous bioreactor framework, the instantaneous production of metabolite is described by the quantity

$$ DX \ [g h^{-1}] \qquad (4.21) $$
which, using definitions (4.3), can be expressed as $DV_{\text{ext}}x$. Then, the total metabo-
lite production over an interval $T$ amounts to

$$
\int_{0}^{T} D V_{\text{ext}} x \mathrm{d}t \quad [g]. \qquad (4.22)
$$

We will be first interested in dynamically adjusting the aggregated control $u(t)$ so as to maximize the quantity (4.22) for a fixed dilution rate $D$. The problem is tackled through an optimal control approach, formulated as

$$
\left\{
\begin{array}{ll}
\text{maximize} & J_x(u) = D V_{\text{ext}} \int_0^T x(t) \, dt \\
\text{subject to} & \text{dynamics of } (S_1), \\
& u(\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 $u(t) \in [0, 1]$. We first note that, since $D$ and $V_{\text{ext}}$ are constants, the problem reduces to maximizing the integral of the concentration $x$. Thus, the nature of the OCP allows us to obtain a first result on the existence of the solutions.

**Proposition 4.3.1.** The dynamic maximization problem has at least one solution.

*Proof.* Since there are no terminal conditions, the set of admissible controls is not empty (any constant control within the prescribed bounds is admissible). For bounded controls ($u(t)$ belongs to $[0, 1]$) and for a fixed final time $T > 0$, the dynamics $(S_1)$ cannot blow up in finite time so all trajectories remain in a fixed compact. Indeed, all state variables but $\mathcal{V}$ are bounded (as set (4.13) is invariant), and since

$$
\dot{\mathcal{V}} = (\mu(p, r) - D)\mathcal{V}
$$

with $\mu(p, r)$ bounded, the volume has at worst an exponential rate and is also bounded. As the dynamics is affine in the control, the set of velocities is convex and existence holds by Filippov's theorem [68]. $\square$

In order to further explore the solution of (OCP), we define the adjoint state vector $\lambda \doteq (\lambda_s, \lambda_p, \lambda_r, \lambda_x, \lambda_\nu)$ associated to the state vector $\varphi \doteq (s, p, r, x, \mathcal{V})$, and
we write the Hamiltonian

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

where $F$ represents the right-hand side of system (S$_1$) given by equations (S$_1$). Developing the expression we obtain

$$H = \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)) \\ + \lambda_r (u-r)\mu(p, r) + \lambda_x (v_X(p, 1-r)\mathcal{V} - Dx) + \lambda_\nu (\mu(p, r) - D)\mathcal{V} - \lambda_0 x,$$

which shows that the Hamiltonian is linear in the control $u$, so it can be rewritten in the affine form

$$H = H_0 + uH_1,$$

where

$$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)) \\ - r\lambda_r\mu(p, r) + \lambda_x (v_X(p, 1-r)\mathcal{V} - Dx) + \lambda_\nu (\mu(p, r) - D)\mathcal{V} - \lambda_0 x,$$

$$H_1 = \lambda_r \mu(p, r).$$

In the absence of terminal constraints, there are no abnormal extremals and $\lambda_0$ can be set to $-1$. Since the constrained optimal control $u$ should maximize the Hamiltonian, one has

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

These alternatives account for the possibility of bang and singular arcs, an extremal being in general an arbitrary concatenation of these. Bang controls correspond to pure allocation strategies $u = 0$ and $u = 1$ (i.e., purely geared towards either the metabolic machinery $M$ or the gene expression machinery $R$ respectively), while singular control $u_s(t)$ occurs when the function $H_1$ identically vanishes over some subinterval. Such singular arcs can be further described by successively differentiat-
ing the switching function $H_1$ until the singular control $u$ can be explicitly computed. For the present case, this occurs when differentiating four times ("order two" singular arc). More precisely, the singular arc is of local (not intrinsic) order two (see [69]). In this case, although chattering (a.k.a. Fuller phenomenon) in and out is not necessary to enter and exit the singular arc, it is indeed what will be observed on the numerical simulations.

**Theorem 3.** Any singular arc $u_s(t)$ of a normal extremal process solution of (OCP) is at least of order two.

*Proof.* We assume that the process $(\varphi, u, \lambda)$ that satisfies PMP is a normal extremal, and we set $\lambda_0 = -1$. In order to further describe the singular arc, we assume $H_1$ vanishes on a whole sub-interval $\tau = [t_1, t_2] \subset [0, T]$. Then, the switching surface is the set

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

where, for this case, $n = 5$. The first derivative is

$$ \dot{H}_1 = \dot{\lambda}_r \mu(p, r) + \lambda_r \left( \mu_r(p, r) \dot{r} + \mu_p(p, r) \dot{p} \right), \quad (4.23) $$

where

$$ \mu_r(p, r) = \frac{\partial \mu(p, r)}{\partial r}, \quad \mu_p(p, r) = \frac{\partial \mu(p, r)}{\partial p}. $$

Along the singular arc $H_1$ is identically zero, so expression (4.23) also vanishes. In order to compute the successive derivatives of $H_1$, we will resort to the Poisson bracket operator [70],

$$ \{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). $$

Applying the latter definition to the derivative of $H_1$ we obtain

$$ \begin{aligned} \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 + uH_1, H_1\} = \{H_0, H_1\} \end{aligned} $$
as {$uH_1, H_1$} = $u\{H_1, H_1\} = 0$. In this context, we will use the notation $H_{01}$ to refer to {$H_0, H_1$}, and so forth. The second derivative of $H_1$ is

$$ \ddot{H}_1 = \dot{H}_{01} = \{H, H_{01}\} = \{H_0, H_{01}\} + u\{H_1, H_{01}\} = H_{001} + uH_{101}. $$

Showing that $H_{101} = 0$ along the switching surface $\Sigma$ ensures that the singular arc is at least of order 2 (otherwise, a singular arc of order 1 could be computed as $u = -H_{001}/H_{101}$). Since $H_1$ depends only on $\lambda_r, p$ and $r$, the computation of the Poisson bracket reduces to

$$ \begin{aligned} H_{101} = \{H_1, H_{01}\} &= \frac{\partial H_1}{\partial \lambda_r} \frac{\partial H_{01}}{\partial r} - \frac{\partial H_1}{\partial r} \frac{\partial H_{01}}{\partial \lambda_r} + \frac{\partial H_1}{\partial \lambda_p} \frac{\partial H_{01}}{\partial p} - \frac{\partial H_1}{\partial p} \frac{\partial H_{01}}{\partial \lambda_p} \\ &= \mu(p, r) \frac{\partial H_{01}}{\partial r} - \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} \end{aligned} \quad (4.24) $$

where

$$ \frac{\partial H_{01}}{\partial r} = \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). \quad (4.25) $$

Replacing the expression (4.25) in equation (4.24) yields

$$ \begin{aligned} 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} \end{aligned} $$

which is also equal to 0 along the switching surface since every term is multiplied either by $\lambda_r$ or $\dot{\lambda}_r$, both identically zero along a singular arc (see $H_1$ expression), showing that the singular arc is at least of order 2. $\square$

It is noteworthy that the above proof holds for all general flows $v_M(s, 1 - r)$, $v_X(p, 1 - r)$ and $v_R(p, r)$ considered in Assumption 1, with no need to use the defined Michaelis-Menten kinetics. Numerical results shown in Figure 4.3 confirm that, when the conditions for the existence of the interior equilibrium $E_i$ are met, the optimal control strategy consists on a series of bang-bang arcs, and a singular arc that is entered and left through the chattering phenomenon. While, as already mentioned, it is compulsory to enter and leave order two singular arcs through chattering, a more precise description of the structure of the extremal would require a
deeper analysis (for instance to prove that there is only one singular arc). Moreover, when the simulation time is long enough, we observe that the singular control $u_s(t)$ converges to a constant value, eventually taking the system to steady-state. This is related to the so-called turnpike phenomenon [71] that relates the singular arc of the dynamic optimization problem with the solution of the static one. See, e.g., [72] for a preliminary analysis on a similar case.

<figure><img src="image_22.png" /><figcaption>Figure 4.3: Results of the numerical simulations using BOCOP [4], showing the optimal control input $u$ and the state variables. The simulation has been done in a time horizon $T = 80$, with 5000 time steps and Midpoint discretization method. The Ipopt (interior point nonlinear optimizer) arguments are: max_iter = 1000, tol = $1.0e-14$, The initial conditions are fixed to: $s(0) = 0.1$, $p(0) = 0.024$, $r(0) = 0.1$, $x(0) = 0$, $V(0) = 0.2$.</figcaption></figure>

Additionally, we observe that the steady-state approached by the system during the singular regime maximizes the integrand of the cost function, which is in fact the instantaneous production of metabolite described in expression (4.21) (otherwise, it would not be the optimal strategy). Therefore, the static optimal control to which the dynamic one converges along the singular arc is the solution of the optimization
problem

$$
\left\{
\begin{array}{ll}
\underset{\text{subject to}}{\text{maximize}} & J_x(\bar{u}) \doteq x \\
& \text{dynamics of } (S_1), \\
& \dot{\varphi} = 0, \\
& 0 \le \bar{u} \le 1.
\end{array}
\right.
\qquad (\text{SP})
$$

The dynamical optimal solution acts as the gold-standard in terms of what can
be achieved through control techniques. However, implementing such strategy is,
in practice, unfeasible. Consequently, it is possible to design a static suboptimal
strategy consisting of a constant allocation $u_{sp}$ that takes the system to the same
steady-state which maximizes the integrand of the cost function. Numerical sim-
ulations of such strategy are shown in Figure 4.4, where it can be seen that the
area below curves $x_{sp}$ and $x$ in OCP only differ slightly ($< 1\%$ for this particular
simulation). The constant control $u_{sp}$ basically disregards the initial and final bang
arcs, as well as the chattering phenomena, which constitutes a small fraction of the
complete time horizon $T$. Moreover, this fraction gets smaller as $T$ becomes large,
which is a typical feature of the turnpike effect for a specific class of optimal con-
trol problems [71]. This way, the difference between strategies becomes marginal in
long-term production schemes.

Additionally, we provide a numerical computation of the successive derivatives
of the switching function $H_1$ in Figure 4.5. It can be seen that the fourth derivative
$H_{10001} \neq 0$ over the interval where $H_1 = 0$, which shows that the singular arc is
of order 2. Furthermore, we verify that the Kelley (or generalized Legendre-Clebsch
[73]) condition

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

is met along the singular arc, which in this case is equivalent to $H_{10001} < 0, \forall t \in \mathcal{I}$.
A check of this condition, necessary for optimality, is also shown in Figure 4.5. Al-
though there is no available sufficient condition to test local optimality of extremals
with Fuller arcs, verifying the Legendre-Clebsch condition along the singular arc
only ensures that we do not compute a too crude local minimizer. Besides, numeri-
cally only a small finite number of bang arcs are retrieved by the optimizer for the
chattering parts before and after the singular arc. This is usually sufficient to give
<figure><img src="image_23.png" /><figcaption>Figure 4.4: Results of the numerical simulations using BOCOP, comparing the control $u$ solution of the OCP to the solution of the static optimization problem $u_{sp}$ (both with same initial conditions). The parameters for the numerical simulation match the ones used in Figure 4.3. The last plot emphasizes the area below the curves $x_{sp}$ and $x$ in OCP, as they are proportional to the total mass of metabolite produced, which is the quantity to be maximized.</figcaption></figure>

a very good approximation of the solution.

### 4.3.3 Static optimization problem

We have shown that a constant allocation decision $\bar{u}$ represents a simplified alternative to the optimal control solution, a strategy composed of bang arcs, a time-varying singular arc and the chattering artifact. In this section, we will further explore the static optimization problem by adding a second degree of freedom to the problem: the dilution rate $D$. In addition to that, we investigate two objectives: the production of biomass $V$ and of metabolite $X$. The static biomass maximization problem
<figure><img src="image_24.png" /><figcaption>Figure 4.5: Successive derivatives of $H_1$ obtained with BOCOP. The intervals where the functions vanish are highlighted in light red. All functions vanish along the singular arc except for $H_{10001}$, highlighted in green, which is negative as required by Kelley condition (4.26).</figcaption></figure>

$(SP_V)$ can be written as

$$ \left\{ \begin{array}{ll} \text{maximize} & J_V(\bar{u}, D) \doteq DV \\ \text{subject to} & \text{dynamics of } (S_1), \\ & \dot{\varphi} = 0, \\ & 0 \le \bar{u} \le 1. \end{array} \right. \qquad (SP_V) $$

Analogously, the product maximization problem can be defined as

$$ \left\{ \begin{array}{ll} \text{maximize} & J_X(\bar{u}, D) \doteq DX = DV_{\text{ext}}x \\ \text{subject to} & \text{dynamics of } (S_1), \\ & \dot{\varphi} = 0, \\ & 0 \le \bar{u} \le 1. \end{array} \right. \qquad (SP_X) $$

Since we look for the steady-states that maximize each objective, the washout equilibrium $E_w$ can be excluded from the analysis since, as shown in Theorem 2, the
equilibrium corresponds to the steady-state values $V = 0$ and $x = 0$. Therefore, the static problems are reduced to finding the equilibria $E_i$ in terms of the pair $(D, \bar{u})$ that maximize each objective function. Moreover, it can be shown that the optimal solution cannot belong to the boundary of the equilibrium $E_i$

$$ \Theta \doteq \{(\bar{u}, D) \in \mathbb{R}^2 : \bar{\mu}(p_w) = D\}. $$

Indeed, using the definitions (4.8), (4.9) and (4.10), it can be seen that on $\Theta$ there is no bacterial population, as $s_i = s_{in}$, which means that $V_i = x_i = 0$. We formalize the latter reasoning in the following proposition.

**Proposition 4.** Every solution $(D, \bar{u})$ of the static optimization problems $(SP_V)$ and $(SP_X)$ is in the region of existence of $E_i$ given by the condition $\bar{\mu}(p_w) > D$.

When considering the same self-replicator scheme under constant environmental conditions (which could describe fed-batch cultivation) the solution for both static problems corresponds to the steady-state with maximal growth rate [3]. Proposition 4 shows that this is not the case in continuous bioreactors, where maximal growth rate is attained at the boundary set of existence of the equilibrium $E_i$. Two interesting particular cases within $\Theta$ are the pure static allocation strategies $\bar{u} = 0$ and $\bar{u} = 1$: A pure metabolic strategy $\bar{u} = 0$ will lead to a bacterial population with no RNA polymerase (i.e., $r = 0$), which will eventually stop the production of biomass (as $v_R(p, 0) = 0$), leading to washout in the bioreactor. Analogously, allocating all resources to the gene expression machinery will finally empty the metabolic machinery $m$, halting the absorption of substrate from the environment (as $v_M(p, 0) = 0$) and depleting the bacteria from precursors, which will also lead to washout. From this analysis, we can conclude that any optimal steady-state allocation $\bar{u}$ should belong to $(0, 1)$.

We recall that, in continuous bioreactors, the growth rate $\bar{\mu}(p_i)$ is fixed by the dilution rate $D$, as shown in equation (4.8). Figure 4.6 illustrates a numerical analysis of the static problems. It is interesting to notice in both subfigures 4.6a and 4.6b that the model accounts for the classical quasi-linear relation between the maximal growth rate (that lies on the boundary $\Theta$), and the control $\bar{u}$ which regulates the RNA/protein mass ratio of the bacterial population [2]. The latter is a phenomenon which has been first observed experimentally, and later used to develop dynamical self-replicator models for natural and biotechnological purposes
[30, 3, 74]. However, in the present case, the growth rate is fixed through D, so it is not a result of the nutrient quality in the environment, which enables the multivariate approach.

**Biomass maximization objective** Results for this problem are shown in Figure 4.6a. For this case, the curve $\bar{u}_{opt}$ remains over 0.6 for all values of the growth rate, tending to $\bar{u}_{opt} = 1$ as the growth rate goes to 0. This shows that, in order to maximize $V$, the allocation strategy should prioritize the synthesis of macromolecules of the gene expression machinery $R$ over the metabolic machinery $M$, which catalyzes the production of biomass and not the synthesis of metabolites. Moreover, the maximum $DV_i$ is attained through a fairly high growth rate (≈ 85% of the maximal growth rate).

<figure><img src="image_25.png" /><figcaption>(a) Results associated to (SP$_V$) with $J_V = (b)$ Results associated to (SP$_X$) with $J_X = DV$. $DX$.</figcaption></figure>

Figure 4.6: Numerical results for both static problems. The values for the objective functions $J_V$ and $J_X$ are represented through a qualitative colormap. The set $\Theta$ of maximum growth rate delimits the region of existence of the interior equilibrium $E_i$. Additionally, curves $\bar{u}_{opt}(D)$ show the optimal allocation $\bar{u}$ in terms of the dilution rate $D$.

**Product maximization objective** Results for this case are shown in Figure 4.6b. We can see that, in opposition to the first case, the optimal solution requires allocating as much precursors as possible into the metabolic machinery $M$, no matter the value of the dilution rate. In other words, the allocation control $\bar{u}$ should be as low as possible within the region of existence of the equilibrium $E_i$ (but not
in $\Theta$). This result is consistent with the fact that the metabolic machinery $M$ cat-
alyzes the synthesis of metabolite $X$. It is noteworthy that the optimal point $DX_i$ is
accomplished through a rather lower dilution rate $D$ in comparison to the biomass
production case, resulting in a continuous production at a rate of about 35% of
the maximal growth rate. The latter might appear counter-intuitive, as it is well
established in the literature that high dilution rates in continuous bioreactors imply
high production rates. This characteristic can be attributed to the compromise be-
tween allocating resources to the metabolic machinery $M$ and increasing the dilution
rate, linked to the maximum value of the function $c(r)$. In other words, the lower
the dilution rate of operation, the wider the interval of existence of the equilibrium
$(u_{\min}, u_{\max})$, which enables the possibility of further promoting the synthesis of com-
pounds of the metabolic machinery (by artificially lowering the allocation parameter
$\bar{u}$) without going to washout. Figure 4.7a illustrates the resulting synthesis rates for
each solution. In $(SP_X)$, the synthesis rate of the metabolite $v_X$ is increased at the
expense of reducing both precursor and biomass synthesis rates $v_M$ and $v_R$, which is
in large part due to the reduction of the dilution rate $D$. This difference in the flux
distribution impacts directly on the mass quantities inside the bioreactor at steady-
state—depicted in Figure 4.7b—for each problem: for the $(SP_X)$ case, there is a
reduction in the biomass $V$ of only 20% w.r.t. the $(SP_V)$ case. However, we can see
an increase on the amount of metabolite $X$ of about 5 times that of the $(SP_V)$ case.
As already said, the allocation parameter $\bar{u}$ becomes quite lower in the metabolite
synthesis objective. In turn, this inhibits the synthesis of macromolecules, which
translates into a decrease of the bacterial population’s growth rate. To compensate
this effect, the steady-state pool of precursors $P$ required for producing the metabo-
lite $X$ becomes considerably bigger than that of the biomass production objective
(about 10 times), therefore increasing the rate of fabrication of biomass $v_R(p, r)$.

## 4.4 Discussion

The objective in this paper was to synthesize a certain metabolite of interest by
re-adjusting the natural allocation of resources in bacteria. This is done by drawing
resources from the native pathways of the cell originally used for producing biomass,
and allocating them into the production of the compound X. We addressed the
problem in the continuous bioreactor framework, which allows a steady-state pro-
<figure><img src="image_26.png" /><figcaption>Figure 4.7: Numerical results for both static problems.</figcaption></figure>

duction regime. Based on previous dynamical systems approaches [30, 3, AYT5], we proposed a coarse-grained self-replicator model capable of accounting for well-studied bacterial growth laws in a simplified way, and we studied its asymptotic behavior. We tackled two different production objectives: Biomass maximization $V$ and metabolite maximization $X$, and we compared results in order to better understand the potential control strategies required to that effect. We rely on novel bio-engineering techniques capable of delivering groundbreaking control schemes: A synthetic growth switch that allows to control the transcription of RNA polymerase through an inducible promoter. In addition to this regulation mechanism, we include the dilution rate as a control input, which yields a multi-variable optimization problem. We concluded by showing very contrasting results, but in accordance with our previous understanding of microbial resource allocation. The biomass-oriented strategy involves an almost maximal dilution rate, and prioritizes investing resources into the gene expression machinery. Conversely, producing the compound of interest requires a rather low value of the dilution rate, which allows an allocation strategy more geared towards the synthesis of components of the metabolic machinery (i.e., a lower value of $\bar{u}$). The latter shows there exist a compromise between augmenting the dilution rate and artificially boosting the production of enzymes and transporters.

The capacity to interfere with the natural bacterial behaviors at molecular level represents a promising tool to improve biotechnological processes. We are interested in further exploring the details of the implementation of such control schemes, by considering different models of the external growth switch. A next step in this regard
would involve taking into account the nature of the external signal *I* and its physical
constraints, in order to adapt our current results towards a more implementable con-
trol loop. Additionally, our approach is based on very simplified representations of
bacterial growth. Such biological processes usually involve numerous reactions and
variables. In our case, we clustered all macromolecules into only 2 different classes,
and we purposely omitted a number of known phenomena in bacteria, such as cell
division, protein degradation and the influence of temperature. However, some of
these effects have been proven to affect only marginally the results regarding the al-
location problem [75]. Thus, in our case, a simplified representation becomes useful
to emphasize the effects of optimally dealing with the internal resource distribution
of bacteria in industrial frameworks. Eventually, such optimal strategies could pro-
vide guidance for developing online feedback solutions based on real-time measuring
of the process.
# Chapter 5

## Model-predictive control schemes

This chapter reproduces [AYT3], accepted for publication in the 60th Conference on Decision and Control (IEEE CDC 2021).

### 5.1 Introduction

Microorganisms have evolved over millions of years under natural selection, submit-
ted to a continuous optimisation process that has improved their capacity to prolifer-
ate in nature. Thus, they have developed highly optimized distribution mechanisms
of their internal resources to cellular functions enabling them to face changing envi-
ronments. Unraveling these internal mechanisms has always been of great interest
for the scientific community, not only from a pure biological point of view, but also
for biotechnological purposes. In this context, being able to understand and control
the growth process is key for several industrial applications, such as in combating
antibiotics resistance, food preservation, and biofuel production [7].

Considering the microbial self-replication process as a resource allocation prob-
lem is a novel approach that has successfully answered some of the underlying ques-
tion in the field [30]. The latter has also motivated numerous applications to the
artificial production of metabolites of interest [3, AYT5, AYT4, AYT2]. These stud-
ies aim to find how to divert the cell internal resources into a heterologous pathway
in order to efficiently synthesize a specific protein. This is done through an external
control that is able to disrupt the cellular allocation process of a growing culture by
reengineering the transcriptional control of the expression of RNA polymerase [44].
In a dynamical systems framework, the problem can be posed as an Optimal Con-
trol Problem (OCP), which can be approached through the well-known Pontrjagin's
Maximum Principle (PMP).

Model-based optimal control studies are essential in understanding the overall
allocation process, as they are able to provide the gold-standard strategies, i.e. the
best that can be achieved from a theoretical point of view. However, in most cases, it
is impossible to obtain a closed-loop control strategy: the obtained optimal control
often depends on the so-called adjoint state, which hinders its implementation (as
it is the case in [3]). Additionally, such approaches depend on the accuracy of
the model and the precision of its parameters, which often tend to be limited for
most biochemical and biological processes. At the same time, the existing industrial
applications that allow a closed-loop implementation are mainly based on general
schemes such as non-linear Model Predictive Controllers (NMPC), which tend to
disregard the structure of each particular problem [76].

Motivated by the lack of synergy between pure theoretical approaches and very
general implementations, in this work we revisit the metabolite production prob-
lem. We summarize the open-loop optimal allocation strategies found in the litera-
ture, which are characterized by sharing the same simple structure and a common
parametrisation. Then, we propose a hierarchical NMPC scheme designed on the
basis of these open-loop optimal controllers. In particular, we resort to the shrinking
horizon NMPC (sh-NMPC) [77], an approach targeted to control processes of known
time duration, such as batch processes [78]. In contrast to typical receding horizon
approaches, in the sh-NMPC, the final time of the process is fixed, and so the time
window considered in the optimisation problem (i.e. the prediction horizon) shrinks
at each step.

Our approach uses an optimal control-based input parametrisation that takes
into account the structure of the open-loop natural allocation strategy of the cell.
Thus, we first implement an MPC loop which creates a closed-loop natural alloca-
tion, followed by a second MPC that computes the external control maximising the
production of a metabolite of interest. Similar approaches have been proposed to
control batch and semi-batch processes through the manipulation of the feedrate
[77, 79], which is a standard scheme in the bioreactor framework. The novelty of
this work resides in a hierarchical control scheme that aims to affect the internal
pathways of the cells in a bacterial growing culture by affecting the expression of
RNA polymerase.
<figure><img src="image_27.png" /><figcaption>Figure 5.1: Coarse-grained self-replicator model. The external substrate S is consumed by bacteria and transformed into precursor metabolites P through the action of the metabolic machinery M. These precursors are used to produce macro-molecules of the gene expression machinery R, the metabolic machinery M, the housekeeping machinery Q, and metabolites X. The external control I is able to externally affect the natural allocation parameter $\alpha$ in order to channel resources into the production of metabolites of interest.</figcaption></figure>

We start the paper by defining the model, and the naturally-evolved resource allocation strategy used by the bacteria. Then, we propose a suboptimal parametrization intended to emulate the open-loop strategy. In Section 5.3, we introduce the open-loop metabolite maximisation problem, and the closed-loop hierarchical scheme. Finally, we provide a numerical simulation of the approach and a comparison with the optimal case, followed by a conclusion.

## 5.2 Model definition

### 5.2.1 Self-replicator system

Based on [30], we define the self-replicator system composed of the mass (in grams) of: precursor metabolites P, the gene expression machinery R, the metabolic machinery M, the housekeeping machinery Q, and a metabolite of interest X. As illustrated
in Figure 5.1, substrate S is taken from the environment and transformed into P at rate $V_M$ through a reaction catalyzed by M. Then, the precursors P are transformed into M, Q, R and X at rates $r_{max}(1-\alpha)V_R$, $(1-r_{max})V_R$, $r_{max}\alpha V_R$, and $V_X$, respectively. While the reactions that produce M, Q and R are catalyzed by R, the reaction synthesizing X is catalyzed by M. In short, the ribosomal proteins R are responsible of producing new proteins, and the metabolic proteins M are responsible for the uptake of nutrients into the cell, and the production of metabolites X. The latter represents a classical trade off in synthetic biology, and is modeled through the parameter $\alpha$, defined as a time function with bounds $\alpha(t) \in [0, 1]$. The dynamical system is

$$ 
\begin{cases}
\dot{S} = V_{in} - V_M \\
\dot{P} = V_M - V_X - V_R, \\
\dot{R} = r_{max}\alpha V_R, \\
\dot{M} = r_{max}(1 - \alpha)V_R, \\
\dot{Q} = (1 - r_{max})V_R, \\
\dot{X} = V_X.
\end{cases}
$$

where the time variable $t$ is measured in hours. The bacterial volume is defined as $V \doteq \beta(R + M + Q)$, and the growth rate given by $\mu \doteq \dot{V}/V$. We define the intracellular concentrations

$$ p = \frac{P}{V}, \quad r = \frac{R}{V}, \quad m = \frac{M}{V}, \quad q = \frac{Q}{V} $$

and the extracellular concentration of substrate $s$. Using the definition of bacterial volume, we obtain the relation $\beta(r + m + q) = 1$. Then, following [16], we assume the transcription of proteins $Q$ to be internally autoregulated to a constant value, such that

$$ \beta(r + m) = r_{max}, \quad \beta q = q_{max} \doteq 1 - r_{max}. \qquad (5.1) $$

We define the rates of mass flow per unit volume, which we assume to be functions of the concentrations $s$, $m$ and $r$, as $v_M(s, m) \doteq V_M/V$ and $v_R(p, r) \doteq V_R/V$. In this new system, the growth rate becomes $\mu = \beta V_R/V = \beta v_R(p, r)$. Taking into account that a minimal concentration of ribosomes $r_{min}$ is required in order for bacteria to self replicate, we define the kinetics of the problem as $v_M(s, m) \doteq w_M(s)m$,
$v_R(p, r) \doteq w_R(p) (r - r_{\min})$ and $v_X(p, m) \doteq w_X(p) m$, with $w_R(p) \doteq k_R p/(K_R + p)$, $w_M(s) \doteq k_M s/(K_S + s)$ and $w_X(p) \doteq k_X p/(K_X + p)$. We will model a production process in which the substrate remains constant. This could be the result of an external control regulating through an inflow of fresh medium to the bioreactor, or due to high availability in the medium. Thus, we replace $w_M(s) = e_M$, with $e_M > 0$ constant. We define the non-dimensional timescale $\hat{t} = k_R t$, as well as the mass fractions of the total volume $\mathcal{V}$: $\hat{p} \doteq \beta p$, $\hat{r} \doteq \beta r$, $\hat{r}_{\min} \doteq \beta r_{\min}$, $\hat{m} \doteq \beta m = r_{\max} - r$, $\hat{q} \doteq \beta q = 1 - r_{\max}$. Additionally, we define non-dimensional synthesis rates $\hat{w}_R(p) = w_R(p)/k_R$, $\hat{w}_X(p) = w_X(p)/k_R$, and parameter $E_M \doteq e_M/k_R$. Then, dropping all hats, the model becomes

$$ 
\left\{ 
\begin{aligned}
\dot{p} &= E_M (r_{\max} - r) - w_X(p)(r_{\max} - r) \\
&\quad - (p + 1)w_R(p)(r - r_{\min}), \\
\dot{r} &= (r_{\max}\alpha - r)w_R(p)(r - r_{\min}), \\
\dot{X} &= w_X(p)(r_{\max} - r)\mathcal{V}, \\
\dot{\mathcal{V}} &= w_R(p)(r - r_{\min})\mathcal{V},
\end{aligned}
\right. 
\qquad (5.2) 
$$

where $q$ and $m$ have been removed using equations (5.1). The parameter values of the kinetics and of bounds $r_{\min}$ and $r_{\max}$ are fixed based on previous studies [30, 3, 23].

### 5.2.2 Naturally-evolved resource allocation strategy

A common assumption in biology is that microorganisms have evolved resource allocation strategies that maximise their growth rate, which allow them to outgrow competing organisms. Such assumption can be represented by an OCP, in which the objective is to maximise the synthesis of biomass in an interval of time $T$ given by $\Delta\mathcal{V}(T) = \mathcal{V}(T) - \mathcal{V}(0)$. This defines the cost function

$$ J_N(\alpha) = \int_0^T w_R(p)(r - r_{\min})\mathcal{V} dt. $$
Thus, as neither the states nor the cost function depend on variable $X$, we will define the OCP for the reduced state $(p, r, \mathcal{V})$ with dynamics

$$ egin{cases} 
\dot{p} = E_M (r_{\max} - r) - w_X(p)(r_{\max} - r) - (p+1)w_R(p)(r - r_{\min}), \\ 
\dot{r} = (r_{\max}\alpha - r)w_R(p)(r - r_{\min}), \\ 
\dot{\mathcal{V}} = w_R(p)(r - r_{\min})\mathcal{V}, 
\end{cases} 
\quad (S_N) $$

and initial conditions

$$ p(0) = p_0, \quad r(0) = r_0 \quad \mathcal{V}(0) = \mathcal{V}_0. \quad (\text{IC}) $$

with $p_0 > 0$, $r_0 \in (r_{\min}, r_{\max})$ and $\mathcal{V}_0 > 0$. The OCP is then defined as

$$ egin{cases} 
\underset{\alpha}{\text{maximise}} & \text{biomass production } J_N(\alpha), \\ 
\text{subject to} & \text{dynamics } (S_N), \\ 
& \text{initial conditions } (\text{IC}), \\ 
& \alpha(\cdot) \in \mathcal{U}, \\ 
& t \in [0, T]. 
\end{cases} 
\quad (\text{OCP}_N) $$

where $\mathcal{U}$ is the set of admissible controllers, which are Lebesgue measurable real-valued functions defined on the time interval $[0, T]$ and satisfying $\alpha(t) \in [0, 1]$. In [30], the particular case where $w_X(p) = 0$, $r_{\max} = 1$ and $r_{\min} = 0$ has been studied, and similar analyses have been carried out in [3, AYT5]. By application of the PMP, it is possible to show that the optimal control that solves $\text{OCP}_N$ has bang and singular arcs, where the singular arc corresponds to the solution of the static optimal control problem obtained through the addition of the constraint $(\dot{p}, \dot{r}, \dot{\mathcal{V}}) = 0$. The solution is characterized by the presence of the Fuller phenomenon before and after the constant singular arc, which produces an infinite number of bangs (also known as *chattering*), a feature that, due to obvious physical limitations, is not possible to implement (nor expected to be found in nature). An example of this kind of structure is shown in Figure 5.2.

**Remark 5.2.1.** Problem $\text{OCP}_N$ can be further simplified by considering the cost function $\ln \mathcal{V}(T)$ instead of $\mathcal{V}(T)$. Then, neither the dynamics nor the cost function depend on $\mathcal{V}$, so the problem can be rewritten in terms of the state $(p, r)$ [30].
<figure><img src="image_28.png" /><figcaption>Figure 5.2: Optimal control $\alpha$ obtained with Bocop [4]. Simulation in a rich medium with $e_M = k_R$, meaning the medium enables the maximum growth rate. Initial conditions are $p_0 = 0.024$, $r_0 = 0.2$, and $V_0 = 0.003$, and the simulation time is set to $T = 15$.</figcaption></figure>

### 5.2.3 MPC parametrisation of the natural allocation

In order to incorporate the natural allocation strategy of the cell into the MPC loop, we propose a sub-optimal parametric form of $\alpha$ given by

$$ \alpha_{so}(\theta, t) = \begin{cases} b_1 & \text{if } t < t_1, \\ \alpha^* & \text{if } t_1 \le t \le t_2, \\ b_2 & \text{if } t > t_2, \end{cases} $$

with the set of parameters $\theta = (b_1, b_2, t_1, t_2, \alpha^*)$ subject to

$$ \begin{array}{l} b_1 \in \{0, 1\}, \quad b_2 \in \{0, 1\}, \\ t_f \ge t_2 \ge t_1 \ge 0, \quad 1 \ge \alpha^* \ge 0, \end{array} \qquad (5.3) $$

where $b_1$ and $b_2$ are Boolean parameters. The suboptimal parametric allocation $\alpha_{so}$ deliberately neglects the chattering artifact from the optimal control $\alpha$, replacing it by pure bang controls during the intervals $[0, t_1)$ and $(t_2, T]$. In order to compare the performance of the proposed controller, we write an optimisation problem with the same biomass production objective $J_N$. At each time instant of the control loop, the algorithm finds the vector of parameters $\theta$ that maximises the final volume of biomass $V(T)$. The latter amount to solving four optimisation problems in terms of $(t_1, t_2, \alpha^*)$ given by all possible combinations of Boolean parameters $(b_1, b_2)$. At each
iteration k, the control loop starts by measuring the system and getting an estimation ($\tilde{p}_k, \tilde{r}_k, \tilde{V}_k$) of the system state. Thus, the optimisation problem at iteration k is formulated with initial conditions

$$ (p(k\tau), r(k\tau), \mathcal{V}(k\tau)) = (\tilde{p}_k, \tilde{r}_k, \tilde{V}_k). \qquad (5.4) $$

This defines the optimisation problem

$$
\left\{
\begin{array}{ll}
\underset{\theta}{\text{maximise}} & \text{biomass production } J_N(\theta) \\
\text{subject to} & \text{dynamics (S_N),} \\
& \text{initial conditions (5.4)} \\
& \alpha(\cdot) = \alpha_{so}(\theta, t), \\
& \text{input constraints (5.3)} \\
& t \in [k\tau, T]
\end{array}
\right.
\qquad (\text{OP}^k)
$$

which is solved at each instant $k\tau$. The scheme proposes a closed-loop form of the open-loop optimal control found in Subsection 5.2.2, with the purpose of implementing it in the hierarchical MPC loop.

### 5.2.4 Numerical example

Figure 5.3 shows a comparison of the optimal control $\alpha$ and the proposed suboptimal control $\alpha_{so}$. The initial and final Fuller arcs are approximated by pure bang arcs (which are $\alpha = 1$ for this particular case), and the parameter $\alpha^*$ of the suboptimal control takes exactly the same value of the static optimal control $\alpha_{opt}^*$. The difference between both control functions is minor, which translates into an imperceptible difference in the trajectories.

## 5.3 Artificial metabolite production

The artificial metabolite production problem is to maximise the synthesis of X over a fixed interval of time $[0, T]$, which is equal to $\Delta X(T) = X(T) - X(0)$, and can be expressed as

$$ J_X(u) = \int_0^T w_X(p)(r_{\max} - r)\mathcal{V} dt. $$
<figure><img src="image_29.png" /><figcaption>Figure 5.3: Comparison of the optimal control $\alpha(t)$ solution of OCP$_N$ and the MPC scheme parametrized with the suboptimal control $\alpha_{so}(\theta, t)$. Initial conditions are set to $p_0 = 0.024$, $r_0 = 0.2$, and $V_0 = 0.003$. The scheme is executed with time step $\tau = 0.3$. The quantity $\Delta X$ amounts to $X(T) - X(0)$.</figcaption></figure>

As neither the states nor the cost function depend on variable $X$, the reduced system ($S_N$) can be used.

### 5.3.1 Optimal Control Problem

In the original approach [3], the naturally-evolved resource allocation parameter $\alpha$ is overridden by the external control $u$, and so the dynamical equation of $r$ becomes

$$ \dot{r} = (r_{\max}u - r)w_R(p)(r - r_{\min}). $$