Optimal Darwinian Selection of Microorganisms with Internal Storage

: In this paper, we investigate the problem of species separation in minimal time. Droop model is considered to describe the evolution of two distinct populations of microorganisms that are in competition for the same resource in a photobioreactor. We focus on an optimal control problem (OCP) subject to a ﬁve-dimensional controlled system in which the control represents the dilution rate of the chemostat. The objective is to select the desired species in minimal-time and to synthesize an optimal feedback control. This is a very challenging issue, since we are are dealing with a ten-dimensional optimality system. We provide properties of optimal controls allowing the strain of interest to dominate the population. Our analysis is based on the Pontryagin Maximum Principle (PMP), along with a thorough study of singular arcs that is crucial in the synthesis of optimal controls. These theoretical results are also extensively illustrated and validated using a direct method in optimal control (via the Bocop software for numerically solving optimal control problems). The approach is illustrated with numerical examples with microalgae, reﬂecting the complexity of the optimal control structure and the richness of the dynamical behavior.


Introduction
The interaction between species coexisting in an ecosystem is complex and affected by external factors. Depending on their environment, some species will dominate, while others, less adapted, will progressively decline. This Darwinian pressure, when it can be manipulated [1], provides the opportunity of guiding the evolution of species of interest. This concept can be applied to artificial ecosystems to select individuals with a desired trait. Here, we focus on microalgae, unicellular photosynthetic microorganisms with promising potential for industrial applications [2,3]. The great biodiversity of microalgae opens the door for a large range of applications [4]. They are grown for their pigments, antioxidants or essential fatty acids [5], and, over the longer term, their efficient way of producing proteins, bricks for green chemistry, biofuel and CO 2 mitigation [2,6,7]. To date, microalgae do not have the place they deserve in biotechnology (see, e.g., [8][9][10]) and many optimization steps must be carried out to improve the economic and environmental performances of these processes at a large scale [6,11]. Currently, only wild organisms sampled in nature are used on an industrial scale. One of the key challenges is to improve the productivity of these strains.
Species in agriculture have been improved after centuries of selection and hybridization. The objective of this work is to develop an alternative approach adapted to microorganisms to select, on a shorter time scale, more productive microalgae strains, by Darwinian where q i is the quota storage of the i-th species and s in is the input substrate concentration. The dilution rate D(·) is a bounded non-negative control function such that D(t) ∈ [0, D max ], where D max > 0 is the maximal admissible value of the dilution rate, above the maximum actual growth rates [22] (this will be made more precise in Section 2.3) of the two species as shown in Figure 1. In addition, for i = 1, 2, ρ i is a non-negative function representing the rate of substrate absorption, i.e., the uptake rate of the free nutrient s, and µ i is also a non-negative function representing the growth rate of the i-th species (see [25]).
Following for instance [25], we suppose that the uptake rates ρ i (s) are expressed as, which corresponds to Michaelis-Menten's kinetics. Here, the parameters K i and ρ m i are positive, i = 1, 2. In addition, we assume that, for i = 1, 2, there is k i > 0 such that cell division does not occur if k < k i . Concerning the Droop model, the kinetics µ i are defined by In addition, for i = 1, 2, let M i stand for, and letq 1 ,q 2 be such that µ i (q i )q i = M i , i = 1, 2 (observe thatq 1 ,q 2 are uniquely defined). Thus, one has, System (1) satisfies the following invariance property.

Proposition 1.
For every q m 1 ≥q 1 , and for every q m 2 ≥q 2 , the set is forward invariant by (1).
Proof. First, observe that, for i = 1, 2, x i never vanishes whenever x 0 i = x i (0) > 0. Now, (0, s in ) is clearly invariant by the dynamics of s(·) sinceṡ ≥ 0 (resp.ṡ ≤ 0) whenever s = 0 (resp. s = s in ). Similarly, for i = 1, 2: where the last inequality follows from the choice of M i and the fact that q m i ≥q i , i = 1, 2. This ends the proof.
The parameter q m i represents the maximum internal storage quota. Since Ω is invariant by (1) (Proposition 1), we notice thatq i stands for the effective maximum internal storage quota for s ∈ (0, s in ). Thus, in the sequel, we consider without loss of generality that q m i =q i for i = 1, 2. In the sequel, we also assume that ρ 1 , ρ 2 fulfill the following hypothesis: Assumption 1. The affinity of species 1 for the substrate is higher than the one of species 2, i.e., K s 2 > K s 1 , or equivalently: From (2), we deduce that, for s ≥ 0, It means that species 1 with the lower K i absorbs nutrients slightly faster.
We are now in a position to formulate the OCP of interest.

Statement of the Optimal Control Problem (OCP)
In this work, we suppose that the first species (with biomass concentration x 1 ) is the one of interest. Our aim is to compute the best feeding strategy, that is, the optimal (dilution rate) control function D(·), such that x 1 becomes predominant in the photobioreactor in minimal-time. This can be formulated and quantified in terms of the ratio between the two competing species. Intuitively, we wish to find an adequate control strategy (if possible optimal) D(·) for which, at the end of the process, we have x 1 Firstly, the set of admissible controls is defined as, where L ∞ loc (R + ) is the space of locally integrable functions on every compact on R + and D max > 0 is the maximum pump feeding capacity. In practice, D max is designed above the maximum growth rates of the coexisting species (see Section 2.3).
To handle the selection process between the two species, let us define a subset T of Ω as, We choose the parameter ε > 0 such that ε 1 in such a way to quantify the contamination rate of the interesting strain x 1 . Whenever a trajectory reaches the target set T , this means that the biomass of the first species is significantly greater than the other one when reaching the target T at the terminal time (if possible). Objective 1. The optimal control problem (OCP) can then be stated: determine a dilution-based control strategy D(·) in such a way that trajectories of (1) starting from an initial condition within the set Ω reach the target set T in minimal-time, i.e., where X(·) is the unique solution of (1) associated with D(·) ∈ D such that X(0) = X 0 ∈ Ω and t D f ∈ (0, +∞] is the first entry time of X(·) into the target set. In the sequel, we will use the simpler notation t f instead of t D f .
In other words, for every positive initial conditions X 0 = (s 0 , q 0 1 , x 0 1 , q 0 2 , x 0 2 ) such that q 0 i > k i , we are seeking an admissible control strategy D = D X 0 ∈ D, steering the trajectory X(t) of the system (1) from X 0 to the target set T in minimal-time, for a fixed D max ( Figure 1) and a given contamination rate ε 1. Note that, if one is able to synthesize such an optimal control for every X 0 ∈ Ω, then one is able to construct an optimal feedback control over Ω as X 0 → D X 0 (0). Such an optimal control problem falls into the class of minimal-time control problems governed by a mono-input affine controlled system, for which the synthesis of an optimal feedback control, thanks to geometric control theory, is a crucial (but also delicate) issue. In particular, handling the high dimension of the Droop model in competition and its resulting optimality system is challenging. Note also that the linearity of the problem w.r.t. D (in contrast, for instance, with strictly convex cost functionals) leads to technicalities because singular arcs usually occur in this setting, see Section 4.

Basic Properties
We now introduce the so-called actual growth rates. These key functions will have an important role in the optimal separation strategy. For that, let us firstly start with the following observations: is well-defined and is also one-to-one; and is one-to-one. From these observations, one can immediately check that the mappings δ i , i = 1, 2 and δ −1 i are increasing. Indeed, for i = 1, one can write δ −1 1 (s) = c 1 + c 2 s+c 3 with c 1 , c 3 > 0 and which implies that δ −1 i is increasing. The actual growth rate of species i is then defined as the mapping The resulting generic functions are illustrated in Figure 1: Let us now define, Throughout the paper, we suppose that ∆ satisfies the following assumption.
Taking into account that ∆(ŝ) = 0, the inequalities satisfied by ∆ according to Assumption 2 can also be written: If we assume that s is regulated to s * (t) = s c using an appropriate control D, we notice that the q-variables are regulated to some unique q ic ∈ [k i , q mi ], for i = 1, 2. The unique point (s c , q 1c , q 2c ), where s c ∈ [0, s in ], plays a crucial role in the optimal control strategy of (OCP) as discussed in Section 5. Finally, D max (the maximum dilution rate) is assumed to be large enough in order to drive competition between the two species. More precisely, we assume that D max satisfies the hypothesis: Assumption 3. The maximal value of the dilution rate D max satisfies ∀s ∈ [0, s in ], D max > max µ 1 (δ −1 1 (s)), µ 2 (δ −1 2 (s)) .
These assumptions will ensure reachability (as detailed in the next section) of the target T , and establishes a generic framework where both species may win the competition for sufficiently large time (considering for instance various constant control parameters D that favor species 1 or 2). Thus, under these considerations, we ensure the well-posedness of the optimal control problem of interest. In this general framework, the objective is then to determine the optimal D that steers the trajectories in minimal-time to the desired target.  Figure 1. The illustrated absorption functions ρ i and growth rates µ i lead to a generic form of the actual growth rates s → µ i (δ i (s)) where both species may win the competition. The maximum dilution rate D max is a fixed constant value above the maximum of the functions s → µ i (δ i (s)) for i = 1, 2, as stated in Assumption 3.

Reachability of the Target
Our next aim is to show that the target is reachable from every initial condition. First, let us recall that the set is an invariant and attractive manifold for (1) for a given persistently exciting control (i.e., an admissible control function D(·) such that +∞ 0 D(t) dt = +∞).

Proposition 2.
For every initial condition X 0 ∈ Ω, there exists an admissible control D(·) and a time t e ≥ 0 such that X(t e ) ∈ T , where X(·) is the unique solution of (1), starting from X 0 , associated with D(·).
Proof. Let s † ∈ (0,ŝ) and X 0 ∈ Ω. Without any loss of generality, we may assume that s(0) = s † . Indeed, observe that s = s † is not a steady-state ofṡ whenever D = 0 over R + or D = D max over R + . Thus, if we apply D = 0 (in that caseṡ < 0) or D = D max (in that casė s > 0), then s(t) = s † is reached in a finite horizon. Consider the feedback control function, in such a way that the unique solution of (1) associated with this control satisfies s(t) = s † for every time t ≥ 0 (Cauchy-Lipschitz's Theorem). We claim that there exists t 1 ≥ 0 large enough such that D † (x 1 (t), x 2 (t)) ∈ [0, D max ] for every time t ≥ t 1 . Indeed, when t → +∞, taking into account (8), it follows: Now, when t → +∞, it is easy to see that q 1 (t) converges to ρ 1 (s † ) . At steady-state, we thus have In conclusion, when t → +∞, or, equivalently, using that, at steady state, Thanks to Assumption 3, this last expression is upper bounded by D max for t large enough, which proves our claim. Finally, posit y i := ln(x i ) and observe thaṫ This ends the proof.

Motivation of Studying the OCP
Thanks to Proposition 2, the target set is reachable from any initial condition, thus the existence of an optimal control of (6) is standard (namely because the dynamics is affine w.r.t. the control): it is an application of the Fillipov Theorem, see, e.g., [26][27][28]. Considering such a control as in the proof of Proposition 2 then indeed allows the system to let the species of interest dominate the reactor, but this process can be long (see, for instance, Example 1). Another possible strategy is to use a constant control D. Following [12], depending on the value of D, species 1 may win the competition, i.e., In that case, this (simple) strategy indeed allows for reaching the target. However, this convergence is asymptotic and depends on the value of D as in the competitive exclusion principle. Roughly speaking, if D > µ 1 (q 1 ) (whereq 1 and s † are such thatq 1 , then species 2 wins the competition, whereas, if D < µ 1 (q 1 ), species 1 wins the competition. We refer to [12] for more details about the asymptotic behavior of (1) for a constant control D. Thus, the target set may not always be reachable with a constant control D. Our objective in this paper is precisely to propose a methodology to compute a control strategy to reach the target set T faster, playing on the control D(·) as illustrated in Figure 2.
Example 1. Let us consider the following parameters: ρ 1m = 0.8, ρ 2m = 0.95, K s 1 = 1, K s 2 = 1.4, The contamination rate is fixed to ε = 0.05. The initial conditions are given by: s 0 = 2, q 0 i = 2.5 and x 0 i = 1. As illustrated in Figure 2, the target T is reached after t f = 46.774 days using the control D(t) = D opt , while x 1 dominates the culture after t f = 62.68 days using the constant control D = 0.48. Let us also point out that, using an arbitrary constant control D ∈ [0, D max ], we are not even sure that x 1 wins the competition (trajectories do not reach the target in that case).  Figure 2. It is possible to find a constant feeding strategy D ∈ [0, D max ] that could favor one species or the other. However, this is a risky time-consuming process, while the target is reached much faster using an optimized control D opt , thus saving several days of costly cultures.

Necessary Conditions on Optimal Controls
We start this section by generalizing previously obtained results [22,29] characterizing optimal solutions of (6). For that, we apply the PMP which allows for obtaining necessary conditions satisfied by optimal controls of (6). We denote by X = (s, q 1 , x 1 , q 2 , x 2 ) and λ = (λ s , λ q 1 , λ x 1 , λ q 2 , λ x 2 ), respectively, the state and adjoint variables (also called co-state or covector). The Hamiltonian associated with the optimal control problem

is given by
Let X 0 ∈ Ω\T and let (X(·), u(·)) be an optimal pair such that X(·) reaches the set T in a time t f ≥ 0. Thanks to the PMP, there exist an absolutely-continuous map The covector satisfieṡ • The Hamiltonian maximization condition writes • At the terminal time, the transversality condition writes: . (11) Here, N T (X) = {p ∈ R 5 ; ∀Y ∈ T , p · (Y − X) ≤ 0} stands for the normal cone to T at some point X ∈ T , see [28]. The adjoint Equation (9) is equivalent to: We define an extremal as a quadruplet (X(·), λ(·), λ 0 , D(·)) such that (λ(·), λ 0 ) is non zero and such that (1) and (9)-(11) is verified. Whenever λ 0 = 0, we say that the extremal is abnormal ; if λ 0 = 0, then we say that the extremal is normal. Because (6) is autonomous (i.e., the system does not depend explicitly on time), the Hamiltonian computed along an extremal is constant. In addition, since the terminal time is not fixed, we classically obtain, following optimal control theory, that H = 0.
The transversality condition is crucial for obtaining properties on optimal controls by reasoning backward in time from the terminal time t = t f . We shall next extend earlier results [22] by taking into account explicitly the fact that T is a half-space of R 5 and exploiting that X(t f ) belongs to the set E := {X ∈ R 5 ; x 2 − εx 1 = 0} (the boundary of the target set). Then, condition (11) can be transformed more explicitly as follows. At X(t f ), the normal cone to T writes Therefore, inclusion (11) is then equivalent to and the inequalities λ Using the constancy of H, we would also obtain λ 0 = 0 contradicting the PMP. We can then conclude that Throughout the paper, we suppose that only normal extremals occur, i.e., λ 0 < 0, and, without any loss of generality, we may assume that λ 0 = −1 (up to a renormalization of the necessary conditions that are linear w.r.t. λ).

Remark 1.
Abnormal extremals are not generic. They correspond to the optimal path reaching the target set in some particular subset of the target set and are such that λ 0 = 0 (and thus λ(·) = 0). From the conservation of H, this implies that µ 1 (q 1 (t f )) = µ 2 (q 2 (t f )) in such a way that λ(t f ) is not uniquely defined in contrast with the normal case (see below). At such a singular point, the value function (the minimal time as a function of X 0 ) is also non-differentiable.
Going back to the normal case, i.e., λ 0 = −1, the covector λ at t = t f can be completely determined (thanks to the conservation of H) as follows: Notice that the quantity µ 1 (q 1 (t f )) − µ 2 (q 2 (t f ))) is non-zero along a normal extremal. As a consequence, the transversality condition (11) coupled with the conservation of H are equivalent to (16). The computation of λ at t = t f is useful to integrate the state adjoint system backward in time from the target set.
We now wish to exploit the Hamiltonian condition (10). It is of particular interest to introduce the switching function, which allows us to determine the optimal control D according to the sign of the switching function. For that, let us denote byΦ the switching function: associated with the control function D(·).
In the case whereΦ > 0, resp.Φ < 0 over a time interval [t 1 , t 2 ], we say that the optimal control u is of bang type (denoted by B + , resp. B − ). When the control D(·) is non-constant in every neighborhood of a time t c ∈ (0, t c ), we say that t c is a switching time, and one must haveΦ(t c ) = 0. Next, when the switching functionΦ vanishes over a time-interval [t 1 , t 2 ], we state that a singular arc occurs. In this case, the corresponding trajectory is singular over [t 1 , t 2 ], and such an arc will be denoted by S. Singular arcs are essential to optimize the time to steer an initial condition to the target set. Now, we are ready to state some main features of the switching functionΦ, and then investigate properties of the singular paths.

Remark 2.
Following the formalism of geometric control theory,˙Φ never involves D explicitly, but D is present in the expression ofΦ (2k) , k ≥ 1, see [27]. If k ≥ 1 is the first integer for which the control is present in the expression ofΦ (2k) , we usually say that the singular arc is of order k.
At this step, an optimal control is a concatenation of bang and singular arcs: with possibly infinitely many crossing times (in particular, if there is a singular arc of order 2 [27]). The occurrence and properties of singular arcs as well as the various (possible) structure for an optimal control of (6) will be precisely the matter of the next section. The goal is to reduce (if possible) the number of possible structures for an optimal control.

Legendre-Clsebsch's Necessary Condition and Computation of the Singular Control
The analysis of singular arcs requires to computeΦ. Indeed, the Hamiltonian condition does not give any information about an optimal control during a singular phase. Thanks to the next computations, we shall also be able to deduce an expression of the singular control during a singular arc.
Doing so, let us define θ 1 , θ 2 : [0, T] → R as Hereafter, to simplify the layout, we do not write the dependency of θ i , s, x i , λ s , λ q i w.r.t. the time. To shorten the notation, we also do not write explicitly the dependency of certain functions w.r.t. some variables. Using (12)- (18), one can writėΦ Lemma 2.
Note that these computations have been verified thanks to a symbolic computation software. The next step is to establish whether Legendre-Clebsch's condition is verified or not along a singular arc. Recall that this condition is necessary for optimality and that it can be stated as follows (see, e.g., [27,30,31]).
Using the expression of the derivativeΦ given in (22), we provide in the next lemma the expression of the second derivative,Φ | D . Lemma 3. Let I = [t 1 , t 2 ] be such that the trajectory is singular over [t 1 , t 2 ]. Then, one has: Proof. In (22), the only term involving the control D is related toṡ. We obtain (24) using that θ 1 ≡ 0 over [t 1 , t 2 ]. Proposition 3. Along a singular arc that occurs over a time interval [t c , t f ], it holds that:

over [t c , t f ] and Legendre-Clebsch's condition (with a strict inequality) is fulfilled over [t c , t f ).
Proof. Because the trajectory is singular over [t c , t f ], one has θ 1 ≡ 0 over this interval; thus, λ s (·) satisfies:λ s = Dλ s , λ s (t f ) = 0.
It follows that λ s ≡ 0 over [t c , t f ]. In a left neighborhood of t = t f , one has from (12) λ q 1 < 0; thus, since λ q 1 vanishes at t = t f , we necessarily have λ q 1 > 0 in a left neighbor- Then, one must haveλ q 1 (t ) ≥ 0 since λ q 1 > 0 over (t , t f ). However, at t = t , the adjoint equation implies thaṫ . This is a contradiction and thus λ q 1 does not vanish over t c , t f . Assumption 1 implies that We can then conclude that Legendre-Clebsch's condition (with a strict inequality) is fulfilled over the whole interval [t c , t f ).
A consequence of the previous proposition is that, when a singular arc occurs over some time interval [t c , t f ], then it is of order 1. Based on this proposition, we shall only consider singular arcs of first order in the remaining of the paper. If Legendre-Clebsch's condition holds true, the singular arc is said to be of turnpike type [26]. The expression defining the singular control can then be derived using (22). Next, let ς(X, λ) be defined by: We now give an expression of the singular control as a feedback of the state and covector.

Proposition 4.
Suppose that an extremal is singular over [t 1 , t 2 ] and that (23) is verified over [t 1 , t 2 ] with a strict inequality. Then, the singular control D s is given by where we recall that and ς is given by (25).

Remark 3.
(i) Note that Legendre-Clebsch's condition (with a strict inequality) is equivalent to θ 2 > 0 (and that this condition is always verified over [t c , t f ) with t c close to t f . (ii) In view of the general expression giving the singular control, see (26), there is no guarantee a priori that the singular control D s is always with values in [0, D max ], i.e., that the singular arc is always admissible (even if Legendre-Clebsch's condition is verified). This can bring additional difficulties; however, we may discard this point by choosing D max large enough. (iii) Notice that (27) is at least active at t = t f in the case where a singular arc D s steers the model trajectories towards the target T , since the transversality conditions ensure that λ s (t f ) = 0.

About the Occurrence of a Terminal Singular Arc at the Terminal Time
The aim of this section is to discuss the possibility of having a singular arc over some time interval [t f − τ, t f ] (with τ > 0) and the structure of optimal controls. Our main questioning is as follows: Does any optimal trajectory contain a singular arc over some time interval [t f − τ, t f ]?
To analyze this point, let us summarize properties of the switching function at t = t f (that are consequences of transversality conditions associated with the codimension 1 target): • The switching function and its derivative vanish at • The second derivative of the switching function satisfies: • In addition, Legendre-Clebsch's condition (23) is always satisfied along a singular arc defined in a left neighborhood of the terminal time t = t f . The necessary conditions (28) and (29) are a very good indication for the occurrence of a singular arc and are thus strong arguments to answer positively to the above question. Thus, we could now wonder whether or not conditions (28) and (29) are sufficient to ensure the occurrence of a singular arc in some time interval [t f − τ, t f ]. It appears that this question is complex and falls into the setting of geometric optimal control theory. As far as we know, such conditions are not equivalent to the occurrence of a singular arc over some time interval [t f − τ, t f ] (this may depend, in particular, on the initial condition). It is, however, worth mentioning that these conditions (in particular (28)) are commonly used numerically to implement a singular arc in shooting methods [32]. In our context of Droop model, it is very interesting to notice that singular arcs are the cornerstone of the optimal control, in particular for a large set of initial conditions that are biologically meaningful (typically for heterogeneous cultures where x 0 1 /x 0 2 ≈ 1). However, the answer to the above question is not always true and depends on the initial condition (as it has been confirmed using direct optimization methods, see Section 5). Indeed, as illustrated in Example 2-Section 5, when the initial conditions x 0 i are taken very close to the target (x 2 (t f )/x 1 (t f ) = ε), and µ 1 (q 0 1 ) − µ 2 (q 0 2 ) > 0, the singular arc does not appear or appears marginally at t = t f to satisfy the transversality conditions.
Recall thatΦ(t f ) =˙Φ(t f ) = 0. Hence, the sign ofΦ depends onΦ(t f ) that is computed in the next lemma.

Lemma 4. At the terminal time, one has
In addition, the second derivative of the switching function exists at t = t f and: Proof. Inequality (30) follows from the expression of λ(t f ) and the transversality condition (16). The expression ofΦ(t f ) in (31) follows from (22).
We can now define the function: From the previous lemma, we deduce the behavior of an optimal path near the terminal time: • First, the target set can only be reached at some point X(t f ) such that (30) is fulfilled. • In addition, if ξ(t f ) = 0, then, in a left neighborhood of t = t f , the optimal control D(·) is of bang type and satisfies D(t) = sign(ξ(t)).
• If a singular arc occurs in a left neighborhood of t = t f , then one must have ξ(t f ) = 0, i.e., a singular arc reaches the target in the subset of T defined as:

Toward an Optimal Synthesis Characterizing the Optimal Solutions
Reducing the number of switching times is in general non-tractable for nonlinear optimal control problems governed by a system in dimension greater than three. Nevertheless, thanks to the properties of singular arcs, we obtained previously and of the switching function at t = t f , we can expect a limited number of possible structures for an optimal control as we formulate in the next conjecture.

Conjecture 1. Every initial condition
in Ω is steered optimally to the target set via a control D that has a finite number of switching times. In addition, for almost every initial condition, an optimal control presents the following structure: For a large set of initial conditions in some subset Ω † ⊂ Ω (far from the target), there is a single bang arc and a terminal singular arc, whereas, for some initial conditions close to the target set, no singular arc occurs (i.e., S is of zero duration).
This conjecture has been verified numerically for a large number of initial conditions (see Section 5). Our argumentation to confirm this conjecture is as follows. • From the PMP, we have seen that, for every initial condition in Ω, an optimal control is a concatenation of bang arcs B ± and singular arcs S. Moreover, since the switching functionΦ satisfies the strong requirements˙Φ(t f ) = Φ(t f ) = 0,Φ | D (t f ) = 0 (from the transversality conditions) as well as Legendre-Clebsch's condition, we conjecture that, for almost all initial conditions, an optimal control is singular in a left neighborhood of the terminal time. This implies in particular that the number of switchings is finite since we proved that any terminal singular arc is of first order. In addition, the number of switchings is minimal in general (apart when chattering occurs, see [27], which is not the case here). Thus, an optimal control should be of type B ± B ∓ S with one or two bang arcs before the terminal singular arc. As we have seen, we also must have ξ(t f ) = 0, which only involves state variables at the terminal time. This surprising condition mixing the first derivative of the basic functions in the Droop model, with the s variable on one side and the q i variables on the other side, is, however, hard to interpret biologically. • For some marginal--but admissible-initial conditions outside of Ω † , the structure of an optimal control of (6) may be of bang type for almost all t ∈ [0, t f ), or [0, t f − τ], with very small τ > 0 (see, e.g., Example 2 in the next section). This is the case when typically x 0 1 x 0 2 , i.e., the initial condition is very close to the target set T , with in addition µ 1 (q 0 1 ) − µ 2 (q 0 2 ) > 0. Thus, the requirement µ 1 (q 1 (t f )) − µ 2 (q 2 (t f )) > 0 is easily satisfied. Thus, in this particular situation, it comes as no surprise that the fastest path to reach the target T is the one exploiting the fact thatẋ 1 (0) It is worth noticing that, when no singular arc occurs, the strategy mainly consists of "pushing" x 1 and x 2 as quickly as possible towards the target T using D = 0, when the initial conditions x 0 1 and x 0 2 are very close to T . Nevertheless, this strategy may not be the optimal one whenever q 0 1 and q 0 2 are "far" from satisfying µ 1 (q 1 (t f )) − µ 2 (q 2 (t f )) > 0. This is typically the case illustrated in Example 3 in the next section.
For an optimal control of type B ± S, the occurrence of a singular arc is related to the socalled turnpike phenomenon that we now explain in this framework. For a large subset of initial conditions S † ⊂ Ω that are biologically the most relevant, the structure of the optimal control is bang-singular B ± S. The singular arc is the control D s given in (27) that reaches the target T . Moreover, this singular phase coincides with optimal trajectories (s(t), q 1 (t), q 2 (t)) that stay most of the time close to the critical point (s c , q 1c , q 2c ) defined in Section 2.3 (related to the actual growth rates and the function ∆(s) = µ 1 (δ −1 1 (s)) − µ 2 (δ −1 2 (s))). Observe, for instance, the trajectories s, q 1 and q 2 in Figure 3a.We also believe that the concatenation of bang arcs before the major singular phase exclusively aims at moving (s 0 , q 0 1 , q 0 2 ) towards (s c , q 1c , q 2c ). Then, the singular arc D s takes over at a switching-time instant t = t s and ensures that the associated singular trajectory, denoted (s s (t), q 1s (t), q 2s (t)), satisfies the so-called turnpike inequality (see, e.g., [33]), This is, for instance, the case for optimal controls illustrated in Figure 4 (of type B − S) and in Figure 5a (of type B − B + S). The inequality (32) usually holds when the time interval [0, t f ] is not excessively short [33][34][35], which is the generic case in the Droop model (1) associated with (6). Indeed, in practice, the most significant biological experiments aim to separate species and select x 1 starting from an homogeneous culture (a well-balanced initial culture with x 0 1 /x 0 2 ≈ 1) or even from x 0 2 x 0 1 with the challenging issue of selecting the minority species (x 1 ), which is not naturally promoted. In these cases, Droop's kinetics ensure that the minimum selection time t f cannot be excessively short and therefore singular arcs as well as the turnpike-type behavior appear systematically in the optimal strategy of (6).  Figure 5. (a) The optimal state trajectories s, q 1 and q 2 (associated with the optimal control D in Figure 6) get closer over time to the critical point (s c , q 1c , q 2c ). The target T , with ε = 0.08, is reached quickly, without resorting to the singular arc. The transversality conditions are satisfied, and in particular λ  6. Evolution of the quantity ρ 2 (s(t))µ 2 (q 2 (t)) − ρ 1 (s(t))µ 1 (q 1 (t)) along the optimal trajectories given in Figure 3a.

Direct Optimization and Numerical Results
In this section, a direct-optimization approach is performed in order to solve (6) and illustrate the different cases discussed in Section 4.3. The numerical direct methods that we use throughout this paper are implemented in Bocop [36] (an optimal control toolbox), and they have the characteristic to transform the studied (6) into a nonlinear programming problem (NLP) in finite-dimension [37], through the discretization step of the control and the state variables [38]. Numerical results are organized as follows: • An optimal control of type B − S is developed throughout Example 1. • An optimal control of type B − is developed throughout Example 2. • An optimal control of type B − B + S is developed throughout Example 3. In all the numerical examples, we consider the model parameters given in Table 1, with the settings in Table 2. Assumption 3 is verified when D max = 1.5, namely because D max is precisely chosen above the maximum actual growth rates of the species. The contamination rate in all the examples is fixed to a significantly small value, ε = 0.08.
In Bocop, the state variables (and even the time, in minimal-time OCPs) of the Droop model (1) are discretized with a Lobatto scheme based on Runge-Kutta methods of type Lobatto-IIIC of order 6, which uses an implicit trapezoidal rule. The main settings used in Bocop are given in Table 3.  Table 1 and Figure 7, associated with the settings in Tables 2 and 3, and the initial conditions given in Table 4.   Table 1, and used throughout the Examples 1-3. We note that the value of s that maximizes the function ∆(s) is given s c = 0.092. This point maximizes the difference between the actual growth rates (as discussed in Section 2.3) and defines the unique static solution (s c , q 1c , q 2c ). All the assumptions that ensure the well-posedness of the generic (6) are satisfied in this case. In particular, we highlight that ∆(s) can be positive and negative, and both species may win the competition using an appropriate control D (see Section 2.4).
In this example, we have a well-balanced initial culture since x 0 1 /x 0 2 = 1 (see Section 4.3). The direct optimization method allows us to determine the optimal control D, given in Figure 4, that steers the model trajectories towards T (with ε = 0.08) in minimal-time t f = 18.3526 days.
We check and analyze the evolution of the switching functionΦ, its derivatives, and the co-state of the substrate s ( Figure 8) in order to characterize the switching instant t s ∈ (0, t f ). This time-instant coincides with λ s = 0 (since the singular arc is the one reaching the target T ),Φ =˙Φ = 0 (thus activating D s , according to the PMP). We also notice that the conditioñ Φ(t f ) = 0 is also satisfied. The optimal state and co-state trajectories are depicted in Figure 3, where we notice that s(t), q 1 (t) and q 2 (t) evolve around the static critical point (s c , q 1c , q 2c ) for almost all t ∈ [t s , t f ], see the turnpike-like property discussed in Section 4.3.
The behavior of the optimal control and optimal trajectories described in Example 2 is definitively the most compelling one (with bang(0)-singular or bang(D max )-singular arcs) from a biological standpoint, since it is the one that systematically appears when the final time is not extremely short. Indeed, in practice, initial conditions start more commonly sufficiently "far" from the target T , leading to a final time that allows the singular arc and the turnpike-like behaviors to hold. Example 3. Now, let us consider the initial conditions in Table 5.  It is worth noticing that the initial conditions in Example 3 are intuitively favourable for reaching the target T in a very short time, since µ 1 (q 0 1 ) − µ 2 (q 0 2 ) > 0 and x 0 2 /x 0 1 = 0.083, while ε = 0.08 (very close to the target). The optimal control in this case is mainly a bang(0) over time, as illustrated in Figure 10 (see Section 4.3 for more details). The optimal state trajectories and co-state trajectories (that satisfy the transversality conditions) are illustrated in Figure 5.

Example 4.
In the last example, let us consider the initial conditions in Table 6.  In this situation, we notice that initial conditions are still favourable for reaching the target T in a very short time because x 0 2 /x 0 1 = 0.0818 (with ε = 0.08, so x 0 i are very close to the target). However, we also note that µ 1 (q 1 (0)) − µ 2 (q 2 (0)) < 0, while µ 1 (q 1 (t f )) − µ 2 (q 2 (t f )) should be positive at t = t f (see Section 3). Thus, the issue of minimal-time separation is slightly more complex than Example 3, and we obtain a structure for the optimal control of bang(0)-bang(1)-singular type, as illustrated in Figure 11. The corresponding model trajectories and optimal co-states are provided in Figure 12.  Figure 11. The optimal control given in (a) is bang(0)-bang(D max )-singular over [0, t f ), where t f = 4.5541 days. A first switching instant from bang(0) to bang(D max ) occurs around 0.3 days, then t s occurs when λ s = 0, as illustrated in (b), starting the singular arc that steers the model trajectories towards T , with ε = 0.08.  Figure 12. The optimal model trajectories (a), and the optimal co-state trajectories (b) that satisfy the transversality conditions.

Conclusions
The minimal-time OCP for the competition between two microbial species accumulating nutrients is a key issue. Progressing along this problem will definitely help for experiments that nowadays last more than 6 months [14,15]. However, the competition described by the Droop model turns out to be significantly more complicated than for the Monod model in dimension 2 [17]. We have improved the preliminary results about this problem to be found in [22,29] in order to provide an optimal synthesis depending on the initial condition. In particular, we applied the PMP, we discussed the structure of the optimal control, and we identified the singular arc steering optimal paths to the desired target set in a minimal amount of time. This study also highlights the turnpike behavior [33] although an exact verification in our case is not possible in our setting since the problem is affine w.r.t. entries in contrast with [35] that, in general, requires coercive hypotheses on the Hamiltonian w.r.t. entries. As usual in optimal control problems that are affine in the control, the study of singular arcs via geometric methods is a crucial issue.
This study also raised a mathematical (open) question outside the scope of this paper on the existence or not of a terminal singular arc whenever the switching function, its derivative, and second derivative w.r.t. the input vanish on a terminal manifold of codimension 1 (see, e.g., [39]).
Future work will focus on the determination of closed loop (sub-optimal) controllers to be applied for bioreactor control subject to uncertainties that are inherently present in biological systems. Finally, the proposed strategy must now be tested experimentally to assess the gain in experimental time it can offer.