# Optimization of microalgae biosynthesis via controlled algal-bacterial symbiosis

Rand Asswad, Walid Djema, Olivier Bernard, Jean-Luc Gouzé, Eugenio Cinquemani

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

## To cite this version:

Rand Asswad, Walid Djema, Olivier Bernard, Jean-Luc Gouzé, Eugenio Cinquemani. Optimization of microalgae biosynthesis via controlled algal-bacterial symbiosis. CDC 2024 - 63rd IEEE Conference on Decision and Control, IEEE CSS, Dec 2024, Milan, Italy. pp.1-6, 10.1109/CDC56724.2024.10886300. hal-04727571

HAL Id: hal-04727571

https://inria.hal.science/hal-04727571v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimization of microalgae biosynthesis via controlled algal-bacterial symbiosis*

R. Asswad<fnref n="1" />, W. Djema<fnref n="2" />, O. Bernard<fnref n="2" />, J.-L. Gouzé<fnref n="3" />,†, E. Cinquemani<fnref n="1" />,†

*Abstract*— We investigate optimization of an algal-bacterial consortium, where an exogenous control input modulates bacterial resource allocation between growth and synthesis of a resource that is limiting for algal growth. Maximization of algal biomass synthesis is pursued in a continuous bioreactor, with dilution rate as an additional control variable. We formulate optimal control in the two variants of static and dynamic control problems, and address them by theoretical and numerical tools. We explore convexity of the static problem and uniqueness of its solution, and show that the dynamic problem displays a solution with bang-bang control actions and singular arcs that result in cyclic control actions. We finally discuss the relation among the two solutions and show the extent to which dynamic control can outperform static optimal solutions.

## I. INTRODUCTION

Microbial consortia, that is, communities of several microbial species in interaction, are common in nature as the result of co-evolution [1]. Such natural symbioses are increasingly exploited to reduce the needs in fertilizers, vitamins or micro nutriments, by having these elements produced within the consortium instead of adding them artificially [2]. A chief example is that of mixed cultures of algae and bacteria. The potential of microalgae for the production of food, feed, cosmetics or pharmaceutics molecules has been highlighted in the last decade [3]. It has been shown that co-culturing algae with probiotic bacteria could lead to a cost decrease, higher robustness of the culture and even growth enhancement [4].

The drawback of using several microorganisms in place of single species is the increased complexity of the resulting dynamics, with challenging issues for process control and optimization [5]. State-of-the-art modelling and experimental capabilities sharpened the mathematical characterization of bioprocess dynamics [6]. While a range of methods for online estimation, optimization and control have been developed over the past decades [7], they are often based on models of the growth of a single species in a bioreactor. Control methods dedicated to microbial consortia are still scarce [8], [9] and much remains to be explored.

In this paper, we focus on the optimal control of an algal-bacterial consortium growing in a continuous bioreactor. We consider the scenario where algal growth is dependent on a

limiting resource, such as vitamins, synthesized by a bacterial population (see also [10] and references therein for consortia displaying similar interactions). Bacterial synthesis of the vitamins is regulated by an exogenous control input (e.g. by the use of optogenetics [11]), and we make the realistic assumption that, as in analogous resource re-allocation scenarios [12], increased synthesis comes at the expense of reduced bacterial growth. Motivated by a synthetic consortium under development comprising microalga *Chlorella* and bacterium *Escherichia coli*, this provides a rich scenario to study control of microbial consortia and maximization of productivity of algal biomass, a proxy for optimized biosynthesis of added value compounds.

We consider two alternative Optimal Control Problems (OCPs): a Static OCP (SOCP), where productivity is to be maximized at equilibrium relative to dilution and vitamin synthesis rate parameters, and a Dynamic OCP (DOCP), where total productivity is to be maximized over a finite time period with time-varying rate functions. For the SOCP, we will first explore equilibria and stability of the algal-bacterial consortium, notably in the conditions for species coexistence. We then investigate problem convexity and, in the same spirit as [10], we numerically illustrate the control tradeoffs and the uniqueness of the solution. For the DOCP, we seek the solution via geometric control theory ([13], [14]) and the Pontryagin's Maximum Principle (PMP, see e.g., [13], [15] and related applications in [16], [17], [18]). We establish a solution in the form of bang-bang control actions and singular arcs. Then, on a simulated case study, we use the numerical optimization software Bocop [19] to illustrate the theoretical results, extend them by showing in particular the cyclic nature of the solution, and relate the solution with that of the corresponding SOCP.

The paper is organized as follows. In Section II we introduce and analyze the algal-bacterial consortium model. The OCPs of interest are stated in Section III. Section IV is dedicated to the analysis and numerical solution of the SOCP, while Section V focuses on the theoretical solution of the DOCP. Simulation results comparing the solutions of the SOCP and of the DOCP are developed in Section VI. Conclusions and perspectives are in Section VII.

## II. MICROBIAL CONSORTIUM MODEL AND ANALYSIS

We present a model that describes the co-culturing of *E. coli* bacteria and *Chlorella* microalgae in a continuously stirred-tank bioreactor. Bacterial growth dynamics are described by the well-known Monod model under glucose limitation [20]. Bacteria are engineered so as to synthesize

<footnote>*Work supported in part by project Ctrl-AB [ANR-20-CE45-0014]</footnote>

<footnote>†These authors contributed equally to the work</footnote>

<footnote n="1">¹MICROCOSME, Centre Inria de l'Université Grenoble Alpes rand.asswad@inria.fr ; eugenio.cinquemani@inria.fr</footnote>

<footnote n="2">²BIOCORE, Centre Inria d'Université Côte d'Azur walid.djema@inria.fr; olivier.bernard@inria.fr</footnote>

<footnote n="3">³MACBES, Centre Inria d'Université Côte d'Azur jean-luc.gouze@inria.fr</footnote>
a vitamin that is a growth-limiting resource for the algal strain. Vitamin synthesis is controlled through optogenetics [11], and it is assumed that allocation of resources to vitamin synthesis slows down bacterial growth ([12] and own experimental evidence). Algal growth dynamics are described by the Caperon-Droop model [21], treating vitamins as the limiting growth substrate. Also known as variable yield model, this model accounts for the import of substrate contributing to an internal quota and subsequent utilization of this quota for algal biomass synthesis.

Let s, e, v, and c be the glucose, bacterial biomass, secreted vitamin and algal biomass concentration in the (fixed volume) bioreactor, in the same order. Let q be the internal algal vitamin quota. Defining the system state $\mathbf{x}(t) = (s(t), e(t), v(t), q(t), c(t))^{\top}$, its dynamics are given by

$$ \dot{s} = -\frac{1}{\gamma}\varphi(s)e + d(t)(s_{\text{in}} - s) \quad (1) $$

$$ \dot{e} = (1 - \alpha(t))\varphi(s)e - d(t)e \quad (2) $$

$$ \dot{v} = \alpha(t)\beta\varphi(s)e - \rho(v)c - d(t)v \quad (3) $$

$$ \dot{q} = \rho(v) - \mu(q)q \quad (4) $$

$$ \dot{c} = \mu(q)c - d(t)c \quad (5) $$

with functions $\varphi$ and $\rho$ defined over $[0, \infty)$, and $\mu$ defined over $[q_{\text{min}}, \infty)$, given by

$$ \varphi(s) = \frac{\varphi_{\text{max}}s}{k_s + s}, \quad \rho(v) = \frac{\rho_{\text{max}}v}{k_v + v}, \quad \mu(q) = \mu_{\text{max}} \left(1 - \frac{q_{\text{min}}}{q}\right). $$

State x takes values in the space $\Omega = \mathbb{R}_+ \setminus \{q < q_{\text{min}}\}$.

We consider throughout the article a constant substrate input flow $s_{\text{in}}$. Constant $\gamma$ represents the bacterial growth yield, $\varphi(s)$ and $\mu(q)$ are the specific growth rates of bacteria and microalgae, respectively. $\rho(v)$ is the vitamin uptake rate, $d(t)$ is the dilution rate and $\beta$ is the vitamin production yield. All the model parameters are strictly positive. Functions $\varphi$ and $\rho$ are called Monod functions as they have the form $t \mapsto at/(b+t)$ defined on $\mathbb{R}_+$ with a and b strictly positive. Factors $1 - \alpha(t)$ and $\alpha(t)$ in Eq. (2)-(3) represent the allocation of resources toward biomass or vitamin synthesis in response to optogenetic control $\alpha(t)$. For simplicity we do not model the optogenetic system response dynamics to light induction (the actual exogenous control input). We simply let $\alpha(t)$ be the control variable and assume we can set it to any value in $[0, 1]$ at any time $t$.

For numerical simulations only, we will refer to the specific choice of parameter values in Table I (ratios of similar units, e.g. $g \cdot g^{-1}$, refer to distinct molecules, e.g. grams of vitamins per gram of biomass). Parameter values pertaining the algal dynamics ($\rho_{\text{max}}$, $q_{\text{min}}$ and $\mu_{\text{max}}$) were determined based on experimental data obtained with collaborators at the Laboratoire d'Océanographie de Villefranche (IMEV, Villefranche-sur-Mer, France). Bacterial growth parameters $\gamma$, $\varphi_{\text{max}}$ and $k_s$ were borrowed from [10], with $\varphi_{\text{max}}$ reduced so as to account for the growth of bacteria at a lower-than-optimal temperature (well below 37° for optimal algal growth). Yield $\beta$ was borrowed from [22].

<table><caption>TABLE I: Model parameters</caption><tr><td>k<sub>v</sub></td><td>0.57</td><td>mg·L<sup>-1</sup></td><td>k<sub>s</sub></td><td>0.09</td><td>g·L<sup>-1</sup></td></tr><tr><td>ρ<sub>max</sub></td><td>27.3</td><td>mg·g<sup>-1</sup>·day<sup>-1</sup></td><td>φ<sub>max</sub></td><td>6.48</td><td>day<sup>-1</sup></td></tr><tr><td>q<sub>min</sub></td><td>2.7628</td><td>mg·g<sup>-1</sup></td><td>γ</td><td>0.44</td><td>g·g<sup>-1</sup></td></tr><tr><td>μ<sub>max</sub></td><td>1.0211</td><td>day<sup>-1</sup></td><td>β</td><td>23</td><td>mg·g<sup>-1</sup></td></tr></table>

### A. Steady states and asymptotic behavior

The system can be studied in two parts: the dynamics of bacterial growth $\mathbf{x}^B(t) = (s(t), e(t))^{\top} \in \Omega_B = \mathbb{R}_+^2$, and those of algal growth $\mathbf{x}^A(t) = (v(t), q(t), c(t))^{\top} \in \Omega_A = \mathbb{R}_+^3 \setminus \{q < q_{\text{min}}\}$. While both systems have been studied separately [20], [21], [23], in this paper a novel consortium system is presented and studied. The model is a cascade autonomous dynamical system

$$ \dot{\mathbf{x}}^B(t) = f_B(\mathbf{x}^B(t)) \quad (6) $$

$$ \dot{\mathbf{x}}^A(t) = f_A(\mathbf{x}^B(t), \mathbf{x}^A(t)) \quad (7) $$

where the bacterial dynamics (6) do not depend on $\mathbf{x}^A(t)$.

1) **Bacterial steady states:** The system (6) has a trivial equilibrium point $\mathbf{x}_0^B = (s_{\text{in}}, 0)^{\top}$ referred to as the bacterial washout equilibrium. A nontrivial equilibrium exists $\mathbf{x}_1^B = (s^*, e^*)^{\top}$ if $d < (1 - \alpha)\varphi(s_{\text{in}})$ with

$$ s^* = \varphi^{-1} \left( \frac{d}{1-\alpha} \right) \quad \text{and} \quad e^* = (1-\alpha)\gamma(s_{\text{in}} - s^*). \quad (8) $$

If $\mathbf{x}_1^B$ exists, it is globally asymptotically stable (GAS) over $\Omega_B \setminus \{e=0\}$ and the washout is unstable, otherwise the latter is GAS over $\Omega_B$ [24].

Let $v_{\text{in}}(\mathbf{x}^B) = \alpha\beta\varphi(s)e/d$, it follows that $v_{\text{in}} = 0$ at bacterial washout and if $\mathbf{x}_1^B$ exists then

$$ v_{\text{in}}^* = v_{\text{in}}(\mathbf{x}_1^B) = \alpha\beta\gamma(s_{\text{in}} - s^*) > 0. \quad (9) $$

This allows reformulating the system (7) in the usual form of a Droop model.

2) **Algal steady states:** Algal equilibria can be investigated by exploring the zero dynamics of (7):

$$ 0 = -\rho(v)c + d(v_{\text{in}} - v) \implies c = d(v_{\text{in}} - v)/\rho(v) \quad (10) $$

$$ 0 = \rho(v) - \mu(q)q \implies \rho(v) = \mu(q)q \quad (11) $$

$$ 0 = \mu(q)c - dc \implies c = 0 \text{ or } \mu(q) = d \quad (12) $$

□ The case of $v_{\text{in}} = 0$ leads to a degenerate Droop model that has a single equilibrium at $(0, q_{\text{min}}, 0)^{\top}$.

□ In the nondegenerate case of $v_{\text{in}} > 0$, the algal system has up to two distinct equilibria that we explore at the bacterial equilibrium $\mathbf{x}_1^B$ for which the vitamin feed is $v_{\text{in}}^*$.

◆ If $c = 0$ the system is at an algal washout equilibrium $\mathbf{x}_0^A = (v_{\text{in}}^*, q_0, 0)^{\top}$ with $q_0 = q_{\text{min}} + \rho(v_{\text{in}}^*)/\mu_{\text{max}}$.

◆ For $c > 0$, the distinct equilibrium $\mathbf{x}_1^A = (v^*, c^*, q^*)^{\top}$ might exist under suitable conditions. From (12), $q^* = \mu^{-1}(d)$ if $d < \mu_{\text{max}}$. Additionally, if $\mu(q^*)q^* < \rho_{\text{max}}$ (i.e. $q^* < q_{\text{max}} = q_{\text{min}} + \rho_{\text{max}}/\mu_{\text{max}}$), then (11) has a unique solution $v^* = \rho^{-1}(\mu(q^*)q^*)$. Injecting in (10) gives $c^* = (v_{\text{in}}^* - v^*)/q^*$ that is well-defined if $v^* < v_{\text{in}}^*$.

For convenience, we introduce the function $\psi$ such that $\psi^{-1}(y) = \rho^{-1}(y \cdot \mu^{-1}(y))$. Hence, $v^* = \psi^{-1}(d)$. It can be shown that $\psi(v) = \psi_{\text{max}}v/(k_c + v)$ with $\psi_{\text{max}} =$
$\mu(q_{\max}) = \mu_{\max}\rho_{\max}/(\rho_{\max} + q_{\min}\mu_{\max})$ and $k_c = k_vq_{\min}\mu_{\max}/(\rho_{\max} + q_{\min}\mu_{\max})$, $\psi$ is therefore a Monod function. Consequently, $\psi$ is increasing and bounded above by $\psi_{\max}$. Thus, $\psi^{-1}(d) = v^* < v_{\min}^* \implies d < \psi(v_{\min}^*) < \psi_{\max} < \mu_{\max}$, this in turn guarantees $q^* < q_{\max}$. This establishes $d < \psi(v_{\min}^*)$ as a necessary and sufficient condition for the existence of $x_1^A$ defined at

$$v^* = \psi^{-1}(d), \quad q^* = \mu^{-1}(d), \quad c^* = \frac{v_{\min}^* - v^*}{q^*}. \qquad (13)$$

Although the vitamin feed $v_{\min}^*$ is meaningful for the algal model and convenient for defining $x_1^A$ and its stability, it carries less significance in the consortium model as it depends on the parameters $\alpha$ and $s_{\min}$. In order to relate the existence of $x_1^A$ to these parameters, we inject (13) and (9) in the inequality $v^* < v_{\min}^*$, which gives $\psi^{-1}(d) < \alpha\beta\gamma(s_{\min} - \varphi^{-1}(d/(1-\alpha)))$. Introducing the increasing function $\psi_{\alpha}$ defined such that

$$\psi_{\alpha}^{-1}(y) = \varphi^{-1}\left(\frac{y}{1-\alpha}\right) + \frac{\psi^{-1}(y)}{\alpha\beta\gamma} \qquad (14)$$

allows rewriting the last inequality as $\psi_{\alpha}^{-1}(d) < s_{\min}$. Subsequently, $x_1^A$ exists if and only if $d < \psi_{\alpha}(s_{\min})$. When $x_1^A$ exists, it is GAS over $\Omega_A \setminus \{c=0\}$ and the washout is unstable. Otherwise, the washout point is GAS over $\Omega_A$. The stability of these equilibria is extensively studied in [23].

3) Consortium's steady states: The algal-bacterial system has up to three equilibrium points: A washout equilibrium for both species always exists at $x_0 = (s_{\min}, 0, 0, q_{\min}, 0)^{\top}$, an algal washout point at $x_{1,0} = (s^*, e^*, v_{\min}^*, q_0, 0)^{\top}$ if $d < (1-\alpha)\varphi(s_{\min})$, and an equilibrium with both species at $x_{1,1} = (s^*, e^*, v^*, q^*, c^*)^{\top}$ if $d < \psi(v_{\min}^*)$. Stability of these equilibria (Table II) results from the stability of (6-7) and the boundedness of all orbits of the full model [25]. Stability is not addressed at bifurcation values for its practical insignificance [24]. From now on, $x_{1,1}$ is referred to as the functional equilibrium, denoted $x^*$, being the equilibrium of interest with nonzero bacteria and algae.

<table><caption>TABLE II: Stability of equilibria over $\Omega \setminus \{e=0, c=0\}$ with $d_1 = \psi_{\alpha}(s_{\min})$ and $d_2 = (1-\alpha)\varphi(s_{\min})$</caption><thead><tr><th></th><th>0 &lt; d &lt; d<sub>1</sub></th><th>d<sub>1</sub> &lt; d &lt; d<sub>2</sub></th><th>d<sub>2</sub> &lt; d</th></tr></thead><tbody><tr><td><b>x<sub>0</sub></b></td><td>unstable</td><td>unstable</td><td>GAS</td></tr><tr><td><b>x<sub>1,0</sub></b></td><td>unstable</td><td>GAS</td><td>-</td></tr><tr><td><b>x<sub>1,1</sub></b></td><td>GAS</td><td>-</td><td>-</td></tr></tbody></table>

### III. OPTIMAL CONTROL PROBLEMS

The objective of our study is to maximize the production of microalgal biomass in a photobioreactor system through two control variables: the vitamin production activation $\alpha(t)$ and the dilution rate $d(t)$. Let $\mathbf{u}(t) = (\alpha(t), d(t))^{\top}$. The produced microalgae are harvested from the bioreactor output which is the function $f_0(\mathbf{x}(t), \mathbf{u}(t)) = d(t)c(t)$. The system dynamics (1)-(5) can be expressed as $\dot{\mathbf{x}} = f(\mathbf{x}, \mathbf{u})$.

The optimal control problem (OCP) for the system in section II is approached in two ways: the static OCP in which the control inputs are chosen to maximize the produced

microalgae biomass at equilibrium, and the dynamical OCP where the total harvested microalgae biomass is optimized over a finite time period. In both cases, the problem consists of choosing $\alpha \in \mathcal{A}$ and $d \in \mathcal{D}$, where $\mathcal{U} = \mathcal{A} \times \mathcal{D}$ is the relevant set of admissible controls, in order to maximize the harvested microalgae from the bioreactor. The two problems are formalized next.

The definition of the dynamic and static OCPs that we propose in this section closely follows the generic formulation found in [26], [27], [28].

#### A. Static Optimal Control Problem (SOCP)

The SOCP consists of choosing admissible control inputs $\mathbf{u} = (\alpha, d)^{\top} \in \mathcal{U}$ that maximize $f_0(\mathbf{x}, \mathbf{u}) = d \cdot c$ at the equilibrium over the set of reachable states $\mathcal{R} \subset \Omega$, i.e.,

$$\max_{(\mathbf{x}, \mathbf{u}) \in \mathcal{R} \times \mathcal{U}} f_0(\mathbf{x}, \mathbf{u}). \qquad (\text{SOCP})$$

#### B. Dynamic Optimal Control Problem (DOCP)

We aim to find optimal control function $\mathbf{u} \in \mathcal{U} = \mathcal{A} \times \mathcal{D}$ to maximize the total production of microalgal biomass over a finite time horizon $[0, t_f]$, i.e.,

$$\max_{\mathbf{u} \in \mathcal{U}} \int_{0}^{t_f} f_0(\mathbf{x}(t), \mathbf{u}(t)) dt \qquad (\text{DOCP}) \\ \text{subject to } \dot{\mathbf{x}}(t) = f(\mathbf{x}(t), \mathbf{u}(t)), \quad \mathbf{x}(t) \in \Omega$$

and verifying the boundary condition $\mathbf{x}(0) \in \Omega$, where $\mathcal{A}$ and $\mathcal{D}$ are defined here as

$$\mathcal{A} = \{\alpha \in L^{\infty}([0, t_f]) | 0 \le \alpha(t) \le 1, \forall t \in [0, t_f]\},$$

$$\mathcal{D} = \{d \in L^{\infty}([0, t_f]) | 0 \le d(t) \le d_{\max}, \forall t \in [0, t_f]\}.$$

### IV. ANALYSIS AND NUMERICAL SOLUTION OF SOCP

As explored in Section II, the system has an asymptotically stable equilibrium $\mathbf{x}^*$ with strictly positive microalgal biomass if and only if $d < \psi_{\alpha}(s_{\min})$. Considering $\mathcal{A} = (0, 1)$ and $\mathcal{D}_{\alpha} = (0, \psi_{\alpha}(s_{\min}))$ for $\alpha \in \mathcal{A}$, we define $f_0^*(\mathbf{u}) = f_0(\mathbf{x}^*(\mathbf{u}), \mathbf{u}) = d \cdot c^*(\alpha, d)$ over $\mathcal{U} = \mathcal{A} \times \mathcal{D}_{\alpha}$ which relaxes the condition $f(\mathbf{x}, \mathbf{u}) = 0$. By injecting the expression of $c^*$ from (8), (9), (13), and (14) we obtain the equivalent problem

$$\max_{(\alpha, d) \in \mathcal{U}} f_0^*(\alpha, d) = \frac{d(v_{\min}^*(\alpha, d) - \psi^{-1}(d))}{\mu^{-1}(d)} \qquad (15)$$

$$= (s_{\min} - \psi_{\alpha}^{-1}(d)) \left( \frac{\alpha\beta\gamma d}{\mu^{-1}(d)} \right) \qquad (16)$$

that is a constrained maximization problem.

#### A. Convexity properties of the problem

A maximization problem is convex if its domain is a convex set and its objective function $g$ is concave (i.e. $-g$ is convex) [29]. The complexity of the expression $f_0^*$ and of the definition of the domain $\mathcal{U}$ makes the problem's convexity analysis cumbersome. However, we can prove weaker properties for (15-16) that still guarantee the existence of global maxima. In particular, we will rely on the notion of biconvexity [30].
<figure><img src="image_3.png" /><figcaption>Fig. 1: SOCP objective function $f_0^*(\alpha, d)$ contours, its maxima along $\alpha$ and $d$, and its global maximum.</figcaption></figure>

**Remark 1.** Every Monod function is strictly concave and its inverse is strictly convex. Both a Monod function and its inverse are strictly positive and strictly increasing.

**Lemma 2.** $f_0^*$ is strictly concave with respect to $\alpha$.

*Proof.* From (15), it is clear that $f_0^*$ is concave with respect to $\alpha$ if and only if $v_{in}^*(\alpha, d) = \alpha\beta\gamma \left(s_{in} - \varphi^{-1}\left(\frac{d}{1-\alpha}\right)\right)$ is concave since $f_0^*$ is affine with respect to $v_{in}^*$ with a positive coefficient $d/\mu^{-1}(d)$. Moreover, the function $\alpha \mapsto d/(1-\alpha)$ is strictly positive, increasing, and convex. Since the composition of an increasing function (here $\varphi^{-1}$) with a strictly convex increasing function is also strictly convex increasing function [29], it follows that $h(\alpha) = s_{in} - \varphi^{-1}(d/(1-\alpha))$ is strictly concave and decreasing. Finally, $\alpha \mapsto \alpha h(\alpha)$ is strictly concave for $\alpha > 0$ because its second derivative $\alpha \mapsto 2h'(\alpha) + \alpha h''(\alpha)$ is strictly negative given $h$ has strictly negative first and second derivatives. Hence, $v_{in}^*(\alpha, d)$ is strictly concave with respect to $\alpha$ and then so is $f_0^*$. ■

**Lemma 3.** $f_0^*$ is strictly log-concave w.r.t. $d$.

*Proof.* Function $\varphi_\alpha^{-1}(d)$ defined in (14) is strictly convex with respect to $d$ for all $\alpha \in \mathcal{A}$ since $\varphi^{-1}$ is increasing and convex, and so is $\psi^{-1}$. It follows that $s_{in} - \varphi_\alpha^{-1}(d)$ is strictly concave, and positive for all $d \in \mathcal{D}_\alpha$ thus strictly log-concave. Moreover, $\alpha\beta\gamma d/\mu^{-1}(d)$ is a concave parabola and strictly positive for $d \in \mathcal{D}_\alpha$, then it is strictly log-concave. From (16), $f_0^*$ is the product of two strictly log-concave functions, hence it is strictly log-concave w.r.t. $d$ [29]. ■

**Proposition 4.** $f_0^*$ is strictly log-biconcave over $\mathcal{U}$.

*Proof.* Since $f_0^*$ is positive and strictly concave with respect to $\alpha$ over $\mathcal{A}$ (Lemma 2), it is strictly logarithmically concave w.r.t. $\alpha$. From Lemma 3, $f_0^*$ is strictly logarithmically concave w.r.t $d$ over $\mathcal{D}_\alpha$. Moreover, $\mathcal{A}$ and $\mathcal{D}_\alpha$ are open intervals therefore trivially convex sets. Subsequently, $\mathcal{U} = \mathcal{A} \times \mathcal{D}_\alpha$ is a biconvex set and $\ln f_0^*$ is a biconvex function. ■

## B. SOCP solution

Proposition 4 establishes the logarithmic biconvexity of the objective function $f_0^*$ over $\mathcal{U}$, which proves the existence, although not the uniqueness, of global maximizers [30]. Note that maximizers of $\ln f_0^*(\mathbf{u})$ are also maximizers of $f_0^*(\mathbf{u})$.

In practice, sequential maximization (w.r.t. $\alpha$ then w.r.t. $d$) produces a single global maximum in Figure 1 at $\bar{\alpha} \approx 0.8251$ and $\bar{d} \approx 0.4409$ per day, which gives approximately 330.786 $mg \cdot L^{-1} \cdot day^{-1}$ of harvested microalgae for the parameters in Table I, with $s_{in} = 0.5 \cdot g \cdot L^{-1}$.

## V. DOCP ANALYSIS: APPLICATION OF THE PMP

In this section, we derive necessary conditions on the optimal controls using the PMP.

### A. Hamiltonian and switching functions

For $\mathbf{x} = (s, e, v, q, c)^\top \in \Omega$, $\mathbf{\lambda} = (\lambda_s, \lambda_e, \lambda_v, \lambda_q, \lambda_c)^\top \in \mathbb{R}^5$, $\mathbf{u} \in \mathcal{U}$, and a scalar $\lambda_0$, we define the Hamiltonian $H(\mathbf{x}, \mathbf{\lambda}, \lambda_0, \mathbf{u}) = \langle \mathbf{\lambda}, f(\mathbf{x}, \mathbf{u}) \rangle + \lambda_0 f_0(\mathbf{x}, \mathbf{u})$, that is,

$$ H = \lambda_0 d(t)c + \lambda_s \left[ -\frac{1}{\gamma} \varphi(s)e + d(t)(s_{in} - s) \right] \\ + \lambda_e \left[ (1 - \alpha(t))\varphi(s) - d(t)]e \right. \\ + \lambda_v \left[ \alpha(t)\beta\varphi(s)e - \rho(v)c - d(t)v \right] \\ + \lambda_q \left[ \rho(v) - \mu(q)q \right] + \lambda_c [\mu(q) - d(t)]c. \quad (17) $$

Note that this can be written as

$$ H = \tilde{H}(\mathbf{x}, \mathbf{\lambda}) + \zeta_{\alpha}(\mathbf{x}, \mathbf{\lambda})\alpha(t) + \zeta_{d}(\mathbf{x}, \mathbf{\lambda}, \lambda_0)d(t) \quad (18) $$

where $\tilde{H}(\mathbf{x}, \mathbf{\lambda}) = -\frac{\lambda_s}{\gamma}\varphi(s)e + \lambda_e\varphi(s)e - \lambda_v\rho(v)c + \lambda_q[\rho(v) - \mu(q)q] + \lambda_c\mu(q)c$ and

$$ \zeta_{\alpha}(\mathbf{x}, \mathbf{\lambda}) = (\beta\lambda_v - \lambda_e) \varphi(s)e, \quad (19) $$

$$ \zeta_d(\mathbf{x}, \mathbf{\lambda}, \lambda_0) = \lambda_s(s_{in} - s) - \lambda_e e - \lambda_v v + (\lambda_0 - \lambda_c)c. \quad (20) $$

### B. The co-state dynamics and transversality conditions

From the PMP, along extremal solutions, there exists an adjoint state $\mathbf{\lambda} : [0, t_f] \to \mathbb{R}^5$ absolutely continuous and a $\lambda_0 \ge 0$ such that $(\mathbf{\lambda}, \lambda_0)$ is non-trivial, satisfying

$$ \dot{\mathbf{\lambda}} = -\partial H/\partial \mathbf{x}. \quad (21) $$

Next, since the final state $\mathbf{x}(t_f)$ is free, the transversality conditions are given by

$$ \mathbf{\lambda}(t_f) = 0. \quad (22) $$

These conditions, combined with the co-state dynamics, form a boundary value problem that must be solved simultaneously with the state equations to determine the optimal trajectories of the controls $\alpha(t)$ and $d(t)$ over $[0, t_f]$.

### C. The PMP maximization condition

The PMP states that for the DOCP, the control strategy maximizes the Hamiltonian $H$ over $[0, t_f]$. Thus the PMP maximization condition requires that, for all $t \in [0, t_f]$,

$$ H(\mathbf{x}, \mathbf{\lambda}, \lambda_0, \mathbf{u}) = \max_{\mathbf{v} \in \mathcal{U}} H(\mathbf{x}, \mathbf{\lambda}, \lambda_0, \mathbf{v}). \quad (23) $$

Given the form of (18), we establish the following result.

**Proposition 5.** For a fixed $t_f > 0$, the optimal controls $\alpha$ and $d$ satisfy, for almost all $t \in [0, t_f]$,

$$ \alpha(t) = \begin{cases} 0 & \text{if } \zeta_{\alpha}(t) < 0, \\ 1 & \text{if } \zeta_{\alpha}(t) > 0, \\ \alpha_{\text{singular}}(t) & \text{if } \zeta_{\alpha}(t) = 0, \ t \in [t_1^{\alpha}, t_2^{\alpha}], \end{cases} \quad (24) $$
$$d(t) = \begin{cases} 0 & \text{if } \zeta_d(t) < 0, \\ d_{\max} & \text{if } \zeta_d(t) > 0, \\ d_{\text{singular}}(t) & \text{if } \zeta_d(t) = 0, \ t \in [t_1^d, t_2^d]. \end{cases} \quad (25)$$

**Proof.** Expressions (24)-(25) follow directly from the application of the PMP maximization condition (23) in the light of (18), which is affine w.r.t. both controls $\alpha$ and $d$, such that the signs of the switching functions $\zeta$ determine the bang controls. If a switching function is zero over a time interval $[t_1, t_2]$ a singular arc may occur. The derivation of the analytic expressions of $\alpha_{\text{singular}}$ and $d_{\text{singular}}$ (which are both of order 1) is not reported due to space limitations. ■

### D. The optimal controls at the final time

We have the following result.

**Proposition 6.** At the final time $t_f$,

$$\zeta_{\alpha}(t_f) = 0 \quad \text{and} \quad \dot{\zeta}_{\alpha}(t_f) = 0, \quad (26)$$

$$\zeta_d(t_f) = \lambda_0 c(t_f) > 0. \quad (27)$$

Moreover, there exists a strictly positive $\varepsilon > 0$, such that for all $t \in [t_f - \varepsilon, t_f]$, $d(t) = d_{\max}$.

**Proof.** The expressions in (26) and (27) follow from those of $\zeta_{\alpha}$ and $\zeta_d$, the costate dynamics (21), and the transversality conditions (22). The inequality in (27) follows from $\lambda_0 > 0$ and $c(t)$ being positive for every $c(0) > 0$. The second statement follows from (27) and the continuity of $c$. ■

Proposition 6 asserts that the optimal dilution strategy $d(t)$ ends with a constant bang arc $d_{\max}$. This makes sense from a biotechnological standpoint since it corresponds to harvesting the remaining algae from the bioreactor. The final stretch of $\alpha(t)$ cannot be characterized in the same manner since a singular arc occurs only if $\zeta_{\alpha}$ is zero on a non-empty interval, which was not confirmed nor excluded.

## VI. NUMERICAL SOLUTION OF THE DOCP AND COMPARISON WITH SOCP

In this section, we employ a direct optimization method that converts the DOCP into a nonlinear programming problem (NLP), wherein the state variables as well as the two control functions $\alpha$ and $d$ are discretized. We implement the DOCP over a large time interval ($t_f = 20$ days) and solve it in Bocop [19] with the following settings: Gauss II discretization (implicit Gauss-Legendre 2-stage scheme of order 4); 7000 time steps; tolerance (Ipot option) $10^{-14}$. In the numerical illustration, we utilize the parameter values of Table I, with $s_{\text{in}} = 0.5 \cdot g \cdot L^{-1}$, and we consider the experimental initial condition $\mathbf{x}_{\text{init}} = (s_0, e_0, v_0, q_0, c_0)^{\top} = (0.1629, 0.0487, 0.0003, 17.7, 0.035)^{\top}$ (in the relevant units).

### A. The optimal optogenetic control $\alpha(t)$

The numerically optimized control $\alpha$ (Figure 2, top) is primarily bang-bang, with short singular arcs between the recurring bangs. The corresponding switching function $\zeta_{\alpha}$ is obtained by injecting the state and costate trajectories into (19). The figure shows correspondence between the sign of the switching function $\zeta_{\alpha}$ and the optimal $\alpha$, with occurrence

<figure><img src="image_4.png" /><figcaption>Fig. 2: Optimal controls $\alpha(t)$ and $d(t)$ [red] and switching functions $\zeta_{\alpha}(t)$ and $\zeta_d(t)$ [blue] over $t \in [0, t_f]$, along with the optimal static controls $\bar{\alpha}$ and $\bar{d}$ [dashed green].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3: The optimal state trajectories plotted against time over 20 days, under the optimal dynamic control $u(t)$ [solid lines] and under the static optimal control $\bar{u}$ [dashed lines]. Concentrations are given in $g \cdot L^{-1}$ and the internal cell quota is given in $g \cdot g^{-1}$</figcaption></figure>

of short singular arcs ($\zeta_{\alpha} = 0$ over short time-intervals). This confirms the properties of the control (24) obtained from the PMP. We also observe $\zeta_{\alpha}(t_f) = 0$ as per Proposition 6. Concerning the comment on $\alpha$ after the proposition, here a singular arc near $t_f$ is not apparent.

### B. The optimal dilution rate control $d(t)$

The numerically optimized control $d$ (Figure 2, bottom) has three major phases: a minimal bang arc $d(t) = 0$ for all $t \in [0, t_1^d)$, then a singular arc $d_{\text{singular}}(t)$ for $t \in [t_1^d, t_2^d]$, followed by a maximal bang arc $d(t) = d_{\max}$ for $t \in (t_2^d, t_f]$ as per Proposition 6. This illustrates the properties of the control (25) obtained from the PMP. The singular arc $d_{\text{singular}}(t)$ exhibits a cyclic-like behavior that appears to be paced by the optimal control $\alpha(t)$. This behaviour will be object of further investigation.

### C. SOCP and DOCP comparison

The static optimal control law $\bar{u}$ obtained in Section IV for the equilibrium state $\mathbf{x}^*$ is compared with the dynamic optimal control law solving the DOCP with initial condition
$x(0) = x_{\text{init}}$ (see Figure 2 for a comparison of the control laws). First, we apply both laws to the system starting from initial condition $x(0) = x_{\text{init}}$. The resulting state trajectories are illustrated in Figure 3 for $t \in [0, t_f]$ and $t_f = 20$ days. In a second test, we apply both laws to the system starting from initial condition $x(0) = x^*$. The performance of the two strategies in terms of total harvested microalgae over the period $[0, t_f]$ is compared in Table III. When the system starts from $x(0) = x_{\text{init}}$ (first row of the table), the dynamic strategy is optimal and it improves upon the static one by approximately 16.85%, which is remarkable for biotechnological applications. The dynamic strategy outperforms the static one (by approximately 5.88%) even when the system starts from $x(0) = x^*$ (second row), the scenario for which the SOCP was solved. We draw the conclusion that the cyclic-like regime of the DOCP solution significantly improves the objective criterion, a phenomenon known as *overyielding* (see [31], [32] and references therein).

<table><caption>TABLE III: Total harvested microalgae over $t_f = 20$ days.</caption><thead><tr><th>x(0)</th><th>static optimal control</th><th>dynamic optimal control</th></tr></thead><tbody><tr><td>x<sub>init</sub></td><td>4.665953 (g/L)</td><td>5.45218 (g/L)</td></tr><tr><td>x<sup>*</sup></td><td>6.615718 (g/L)</td><td>7.005 (g/L)</td></tr></tbody></table>

## VII. CONCLUSIONS

In this study, we discussed static and dynamical optimal control of a microbial consortium model, with the objective to enhance productivity of algal biomass in the interest of biotechnological applications. In particular, our investigation reveals that the dynamic control involves *bang-bang* and singular phases, leading to a significant improvement of the control performance achieved with a static optimal solution at equilibrium. In our future work, we aim to explore the design of robust feedback controls for practical implementation in bioreactors.

## ACKNOWLEDGMENT

We thank Juan-Carlos Arceo-Luzanilla for contributing to establish parameter values, Antoine Sciandra for experimental data, and Inria project-team MICROCOSME for modelling discussions.

## REFERENCES

[1] J. L. Sachs, R. G. Skophammer, and J. U. Regus, "Evolutionary transitions in bacterial symbiosis," *PNAS*, vol. 108, no. supplement 2, pp. 10800-10807, 2011.

[2] K. M. Rapp, J. P. Jenkins, and M. J. Betenbaugh, "Partners for life: building microbial consortia for the future," *Curr. Opin. Biotechnol.*, vol. 66, pp. 292-300, 2020.

[3] M. Rizwan, G. Mujtaba, S. A. Memon, K. Lee, and N. Rashid, "Exploring the potential of microalgae for new biotechnology applications and beyond: A review," *Renew. Sustain. Energy Rev.*, vol. 92, pp. 394-404, 2018.

[4] O. A. Palacios, B. R. López, and L. E. de Bashan, "Microalga growth-promoting bacteria (mgpb): A formal term proposed for beneficial bacteria involved in microalgal-bacterial interactions," *Algal Res.*, vol. 61, p. 102585, 2022.

[5] A. R. Zomorrodi and D. Segrè, "Synthetic ecology of microbes: mathematical models and applications," *J. Mol. Biol.*, vol. 428, no. 5, pp. 837-861, 2016.

[6] A. J. Lopatkin and J. J. Collins, "Predictive biology: modelling, understanding and harnessing microbial complexity," *Nat. Rev. Microbiol.*, vol. 18, no. 9, pp. 507-520, 2020.

[7] R. Simutis and A. Lübbert, "Bioreactor control improves bioprocess performance," *Biotechnol. J.*, vol. 10, no. 8, pp. 1115-1130, 2015.

[8] C. Aditya, F. Bertaux, G. Batt, and J. Ruess, "A light tunable differentiation system for the creation and control of consortia in yeast," *Nat. Commun.*, vol. 12, no. 1, p. 5829, 2021.

[9] D. Salzano, D. Fiore, and M. di Bernardo, "Ratiometric control of cell phenotypes in monoen strain microbial consortia," *J. R. Soc. Interface.*, vol. 19, no. 192, p. 20220335, 2022.

[10] M. Mauri, J.-L. Gouzé, H. De Jong, and E. Cinquemani, "Enhanced production of heterologous proteins by a synthetic microbial community: Conditions and trade-offs," *PLoS Comput. Biol.*, vol. 16, no. 4, p. e1007795, 2020.

[11] M. Benisch, S. K. Aoki, and M. Khammash, "Unlocking the potential of optogenetics in microbial applications," *Curr. Opin. Microbiol.*, vol. 77, p. 102404, 2024.

[12] A. Y. Weiße, D. A. Oyarzún, V. Danos, and P. S. Swain, "Mechanistic links between cellular trade-offs, gene expression, and growth," *PNAS*, vol. 112, no. 9, pp. E1038-E1047, 2015.

[13] H. Schättler and U. Ledzewicz, *Geometric optimal control: theory, methods and examples*. Springer, 2012, vol. 38.

[14] R. B. Vinter and R. Vinter, *Optimal control*. Springer, 2010, vol. 2, no. 1.

[15] L. S. Pontryagin, *Mathematical theory of optimal processes*. Routledge, 2018.

[16] W. Djema, T. Bayen, and O. Bernard, "Optimal darwinian selection of microorganisms with internal storage," *Processes*, vol. 10, no. 3, p. 461, 2022.

[17] I. Yegorov, F. Mairet, H. De Jong, and J.-L. Gouzé, "Optimal control of bacterial growth for the maximization of metabolite production," *J. Math. Biol.*, vol. 78, pp. 985-1032, 2019.

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

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

[20] G. Bastin and D. Dochain, *On-line estimation and adaptive control of bioreactors*, ser. Process measurement and control. Amsterdam ; New York: Elsevier, 1990, no. 1.

[21] H. L. Smith and P. Waltman, "Competition for a Single Limiting Resource in Continuous Culture: The Variable-Yield Model," *SIAM J. Appl. Math.*, vol. 54, no. 4, pp. 1113-1131, Aug. 1994.

[22] Z. Lin, Z. Xu, Y. Li, Z. Wang, T. Chen, and X. Zhao, "Metabolic engineering of Escherichia coli for the production of riboflavin," *Microb. Cell. Fact.*, vol. 13, no. 1, p. 104, Jul. 2014.

[23] K. Lange and F. J. Oyarzun, "The attractiveness of the droop equations," *Math. Biosci.*, vol. 111, no. 2, pp. 261-278, Oct. 1992.

[24] J. Harmand, C. Lobry, A. Rapaport, and T. Sari, *The Chemostat: Mathematical Theory of Microorganism Cultures* | Wiley. Wiley-ISTE, Aug. 2007.

[25] P. Seibert and R. Suarez, "Global stabilization of nonlinear cascade systems," *Syst. Control Lett.*, vol. 14, no. 4, pp. 347-352, Apr. 1990.

[26] E. Trélat and E. Zuazua, "The turnpike property in finite-dimensional nonlinear optimal control," *J. Differ. Equ.*, vol. 258, no. 1, pp. 81-114, 2015.

[27] J.-B. Caillau, W. Djema, J.-L. Gouzé, S. Maslovskaya, and J.-B. Pomet, "Turnpike property in optimal microbial metabolite production," *J. Optim. Theory Appl.*, vol. 194, no. 2, pp. 375-407, 2022.

[28] W. Djema, L. Giraldi, S. Maslovskaya, and O. Bernard, "Turnpike features in optimal selection of species represented by quota models," *Automatica*, vol. 132, p. 109804, 2021.

[29] S. Boyd and L. Vandenberghe, *Convex Optimization*, 1st ed. Cambridge University Press, Mar. 2004.

[30] J. Gorski, F. Pfeuffer, and K. Klamroth, "Biconvex sets and optimization with biconvex functions: A survey and extensions," *Math. Methods Oper. Res.*, vol. 66, no. 3, pp. 373-407, Dec. 2007.

[31] T. Bayen, A. Rapaport, and F. Z. Tani, "Optimal periodic control for scalar dynamics under integral constraint on the input," *Math. Control Relat. Fields.*, vol. 10, no. 3, pp. 547-571, 2020.

[32] —, "Improvement of performances of the chemostat used for continuous biological water treatment with periodic controls," *Automatica*, vol. 121, p. 109199, 2020.