One step proximal point schemes for monotone vector field inclusion problems

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 ^[1]. In this regard, optimality of bacterial organisms has been a central hypothesis in several research allocation studies ^[2].

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 ^[3]. This motivated a series of works with a dynamical-oriented approach to the problem ^[4,5], 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 internal 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 ^[6]. 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) ^[7].

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 ^[8,9]. 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 ^[10,11,12]. In this approach, a control loop is developed through a bacterial growth switch that allows to externally modify the natural resource allocation decision ^[13]. 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 ^[14], 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 ^[10] 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). 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 ^[15]. Examples of self-replicator models in continuous bioreactors can be found in some recent works: In ^[16], 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 ^{[17,18,19,20,21]}. 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 ^[22]. 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) ^[23]. 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 ^[4], and taking into account the constraints for the existence of the equilibrium of interest.

2.   Model definition

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 $\mathcal{V}_\mathrm{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 ^[10] (Figure 1).

The model represents three chemical macroreactions,

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 ^[24]. As already stated, we include in our scheme the growth switch described in ^[13] to externally affect the allocation decision $\alpha$ by varying the inducer concentration in the medium. The mechanism is modeled as

where $I$ is the external signal, which acts independently from the allocation parameter $\alpha$. The agammaegated form of Eq (2.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 = \alpha$). In this work, we limit the analysis to calculating the control input $u$, without decoupling both individual controls $I$ and $\alpha$ (refer to Section 5 in ^[10] 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,

The inflow/outflow rates are defined as

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

We define the volume of the cell population $\mathcal{V}$ [L] as

where $\beta$ [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 ^[25]. Thus, the mass of precursors $P$ is excluded from the volume $\mathcal{V}$ as a simplifying assumption. Then, we express the quantities of the system as concentrations with respect to the volumes,

where $p$, $r$ and $m$ [g L$^{-1}$] are intracellular concentrations of precursor metabolites, components of the gene expression machinery and metabolic enzymes respectively; and $s$ and $x$ [g L$^{-1}$] 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}_\mathrm{ext}$. Then, using definition (2.2), we obtain that $m + r = 1/\beta$. We define the growth rate of the bacterial population $\mu$ [h$^{-1}$] 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

where $v_M(s, m)$, $v_R(p, r)$ and $v_X(p, m)$ [g L$^{-1}$ h$^{-1}$] 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

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

which yields the relation

and we define the non-dimensional synthesis rates

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

where the dynamical expression of $m$ has been omitted since it can be computed from $r$, as shown in Eq (2.5). 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 Eq (2.5). 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:

● $v_i(\cdot)$ bounded w.r.t. $y$: $\lim_{y\rightarrow \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 ^[26]. 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)

where $c(r) : (0, 1) \rightarrow \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 ^[27]. Notably, for fixed values of $r$, both synthesis rates are simply proportional. This assumption reduces $S_1$ to

We note that both Assumptions 1 and 2 are formulated for the non-dimensional synthesis rates given in definition (2.6), 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}] \rightarrow \mathbb{R}^+$).

2.2.   Asymptotic behavior

The asymptotic behavior of system $S_1$ 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_1$ 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_1$ converges to the equilibria of its limiting system. Let us start the analysis of system $S_1$ by defining its invariant region in the following lemma.

Lemma 1. The set

is positively invariant for the initial value problem.

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

We will study the initial value problem of system $S_1$ with initial conditions

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.

2.2.1.   Mass conservation

System $S_1$ can be rewritten as

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

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

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

Moreover, when fixing $u(t)$ to $\bar u \in (0, 1)$, the quantity $w_2$ also obeys the same Eq (2.10), 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$ given by Eq (2.8) lies in the hyperplane

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

Further on, we will use Lemma 2 to analyze system $S_1$ through its limiting system.

2.2.2.   Limiting systems

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

meaning that, as $t \rightarrow \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_1$ becomes

where

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

2.2.3.   Local stability

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

Definition 2. We define the function

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

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) \geq 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

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

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

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

meaning that $p = p_w$ is either repulsive or invariant, so a new invariant set $\Gamma' \subset \Gamma$ can be defined,

Equilibrium $E_i$ The existence of the interior equilibrium $E_i$ is given by the boundedness of the flows,

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 (2.17) 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 (2.18) is met, $s_i$ should be unique. As is standard in continuous bioreactors, inequality (2.17) implies that the maximal growth rate of the bacterial population should be bigger than the dilution rate $D$. At the same time, the inequality (2.18) requires the maximal uptake flow to be bigger than the flows responsible for bacterial growth and metabolite production. In $\Gamma'$, these conditions become

where the second condition (2.18) is included in the inequality (2.19), as $max \; \bar v_M(s) = \bar v_M(s_{in}) = \bar f(p_w)$, and the condition $\bar f(p_w) \geq \bar f(p_i)$ is true if and only if $\mu(p_w) \geq D$. The Jacobian matrix is given by

and the characteristic polynomial is

Using Definition 2, it can be seen that $-D$ is an eigenvalue by replacing in expression (2.20)

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

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

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

with characteristic polynomial

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

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

Proof. The quantity

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

Using Lemma 2, we can further reduce system $S_1'$ by expressing

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

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

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

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

● $E_w^* \doteq (p_w, 0)$ 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:

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

Through the theory of asymptotically autonomous systems ^[22], 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 (2.9) 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) \geq D$.

Proof. We resort to a more particular case of the general theory of asymptotically autonomous systems (introduced in Appendix F of ^[24]), 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 Eqs (2.10) and (2.21) 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 Eq (2.23).

(A5) $S_1''$ has no periodic orbits, as shown in Eq (2.23).

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.

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 ^[4] 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)$.

3.1.   Kinetic's definition

Michaelis-Menten kinetics are defined as

where the biological constants and their values ^[4,10] are described in Table 1.

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

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

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

which, using definitions (2.3), can be expressed as $D \mathcal{V}_\mathrm{ext} x$. Then, the total metabolite production over an interval $T$ amounts to

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

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 $\mathcal{V}_\mathrm{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 3. 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 (2.8) 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 (2.16) is invariant), and since

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 ^[28].

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_ \mathcal{V})$ associated to the state vector $\varphi \doteq (s, p, r, x, \mathcal{V})$, and we write the Hamiltonian

where $F$ represents the right-hand side of system $S_1$ given by Eq (2.8). Developing the expression we obtain

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

where

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

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 differentiating 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), showing that the singular regime should be entered and leaved through Fuller phenomenon or “chattering” ^[29] (an infinite number of switchings between the control bounds $0$ and $1$ over a finite time interval). The theorem hereunder formalizes this statement.

Theorem 4. 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

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

where

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

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

as $\left\{ u H_1, H_1 \right\} = u \left\{ H_1, H_1 \right\} = 0$. In this context, we will use the notation $H_{01}$ to refer to $\left\{ H_0, H_1 \right\}$, and so forth. The second derivative of $H_1$ is

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

where

Replacing the expression (3.5) in Eq (3.4) yields

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.

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 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 ^[31] that relates the singular arc of the dynamic optimization problem with the solution of the static one. See, e.g., ^[32] for a preliminary analysis on a similar case.

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

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 simulations of such strategy are shown in Figure 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 control problems ^[31]. 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 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 ^[34]) condition

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 5. Although 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, numerically 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 a very good approximation of the solution.

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 $\mathcal{V}$ and of metabolite $X$. The static biomass maximization problem (SP$_{ \mathcal{V}}$) can be written as

Analogously, the product maximization problem can be defined as

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 $\mathcal{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$

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

Proposition 5. Every solution $(D, \bar u)$ of the static optimization problems (SP$_ \mathcal{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 ^[10]. Proposition 5 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 Eq (2.11). Figure 6 illustrates a numerical analysis of the static problems. It is interesting to notice in both subfigures 6a and 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 ^[35]. 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 ^[4,10,11]. 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 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 $\mathcal{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 $D \mathcal{V}_i$ is attained through a fairly high growth rate ($\approx 85\%$ of the maximal growth rate).

Product maximization objective Results for this case are shown in Figure 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$ catalyzes the synthesis of metabolite $X$. It is noteworthy that the optimal point $D X_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 between 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 compounds of the metabolic machinery (by artificially lowering the allocation parameter $\bar u$) without going to washout. Figure 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 7b—for each problem: for the (SP$_X$) case, there is a reduction in the biomass $\mathcal{V}$ of only $20\%$ w.r.t. the (SP$_ \mathcal{V}$) case. However, we can see an increase on the amount of metabolite $X$ of about $5$ times that of the (SP$_ \mathcal{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 metabolite $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.   Discussions

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 production regime. Based on previous dynamical systems approaches ^[4,10,12], 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 $\mathcal{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 control 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 allocation problem ^[36]. 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 provide guidance for developing online feedback solutions based on real-time measuring of the process.

Acknowledgments

This work was partially supported by ANR project Maximic (ANR-17-CE40-0024-01), Inria IPL Cosy and Labex SIGNALIFE (ANR-11-LABX-0028-01). We acknowledge the support of the FMJH Program PGMO and the support to this program from EDF-THALES-ORANGE. We also thank Hidde de Jong (INRIA Grenoble - Rhône-Alpes) for discussions on the manuscript.

Conflict of interest

All authors declare no conflicts of interest in this paper.