THE OPERATING DIAGRAM FOR A MODEL OF COMPETITION IN A CHEMOSTAT WITH AN EXTERNAL LETHAL INHIBITOR

. The inhibition is an important phenomenon, which promotes the stable coexistence of species, in the chemostat. Here, we study a model of two microbial species in a chemostat competing for a single resource in the presence of an external lethal inhibitor. The model is a four-dimensional system of ordinary diﬀerential equations. We give a complete analysis for the existence and local stability of all steady states. We describe the bifurcation diagram which gives the behavior of the system with respect to the operating parameters represented by the dilution rate and the input concentrations of the substrate and the inhibitor. This diagram, is very useful to understand the model from both the mathematical and biological points of view.

In this paper we consider the model of Hsu, Li and Waltman [9]. In this model, two species x and y compete for a single limiting resource S in presence of an external lethal inhibitor p, to which the species x is sensitive and the species y is resistant species. This species is able to detoxify the inhibitor, that is to remove it from the chemostat. Let S(t) denote the concentration of the substrate at time t; let x(t), y(t) denote the concentrations of the competitors respectively; finally, let p(t) denote the concentration of the inhibitor. The model takes the form with initial condition S(0) ≥ 0, x(0) > 0, y(0) > 0 and p(0) ≥ 0. The biological parameters of the model are m 1 , m 2 , δ, K 1 , K 2 , K, β 1 , β 2 , and γ. The meaning and units of these parameters are given in Table 6. In particular the constant γ represents the lethal effect of the inhibitor p on x. These parameters are called biological parameters since they depend on the organisms, substrate and inhibitor considered. These parameters are measurable in the laboratory. In contrast, the operating parameters are the input concentration of the nutrient S 0 , the input concentration of the inhibitor p 0 and the dilution rate D of the chemostat. These parameters are called operating parameters since they are under the control of the experimenter. The complete mathematical analysis of (1), including some global results, was given in [9]. The presence of the inhibitor allows the coexistence of both species. Moreover, the coexistence equilibrium can be unstable, giving the possibility of an asymptotic oscillatory behavior. Due to the importance of this phenomenon, which promotes the coexistence of two species competing for a single resource, the results of [9] were discussed also in the review paper [11].
A useful contribution of the mathematical analysis of a chemostat model, is to give to the engineers the operating diagram which is the bifurcation diagram for which the values of the biological parameters are fixed, and the behavior of the model is discussed with respect to the operating parameters. The operating diagram has the operating parameters as its coordinates and the various regions defined in it correspond to qualitatively different dynamics. As it was noticed by Smith and Waltman [21], p. 252, the operating diagram is probably the most useful answer for the discussion of the behavior of the model with respect to the parameters. This diagram shows how robust or how extensive is the parameter region where coexistence occurs, where the coexistence equilibrium is stable and where it is unstable. This bifurcation diagram is very useful to understand the model from both the mathematical and biological points of view, and is often constructed both in the biological literature [17,23,25,13] and the mathematical literature [1,5,19,20,26].
The parameter space of (1) is twelve dimensional: nine biological parameters and three operating parameters. Exploring all of it is not possible. The approach in [9] consists of rescaling both biological and operating parameters of the model, creating a 'standard' environment in which the operating parameters are normalized S 0 = p 0 = D = 1. This approach is often used in the mathematical literature on the chemostat to reduce the huge number of parameters [21]. Of course, the theory developed in this standard environment potentially permits to present the operating diagrams of the model. Actually an operating diagram is presented in [9]. However, since the computations are made in the standard environment, it is not so easy with this approach to understand how the operating diagram is affected by parameter variations.
In (1), the inhibitor is called lethal since it acts on the death of the species x rather than on its growth, as it was the case in the model of Lenski and Hattingh [13], where the model takes the form with a degree of inhibition f (p) on the growth of species x. In this model the inhibitor is decreasing the growth of the species rather than increasing its death rate. The mathematical analysis of the model (2) was given in [10], and due to its importance, it was discussed also in [11,21]. In [15] the authors assume that the external inhibitor is lethal and can decrease the growth rate of species x, so that the equation for x in (2) is replaced by equation x = m1S K1+S f (p) − D − γp x, the equations for S, y and p being unchanged. The study of (2) was recently revisited in [4] and its operating diagram was described. For (2), the system can be reduced to a three dimensional system through a conservation principle. In the case of (1) this reduction is not possible, so it is necessary to work in the full four dimensional state space. Following [4], our approach to handle the question of robustness with respect to parameter variations, is to present the computations for generic growth and detoxification functions, rather than the specific model (1), and to use the fact that the condition of stability of coexistence equilibrium point of the system is polynomial, of degree 3, with respect to the operating parameter S 0 . The problem is then reduced to the determination of the signs of four real valued functions, which are built with the coefficients of this polynomial. These functions depend only on the operating parameters p 0 and D. Hence, combining mathematical analysis and numerical plots, we can construct the operating diagram of the model. The paper is organized as follows. In Sect. 2, the mathematical model (3), extending the model (1), of competition of two species in presence of a lethal external inhibitor is described. The results of [9] carry over to our extended model with only minor notation changes. We describe these results in Prop. 1. In Sect. 3 we present our main result which gives the conditions of existence and stability of the equilibrium points, explicitly with respect to the operating parameters, see Prop. 2. In Sect. 4 we discuss the properties of the region of instability of the coexistence equilibrium, see Prop. 3. The main result in this section is the definition of the regions Υ i , i = 0, .., 3 of the operating parameter space (D, p 0 ) for which the instability of the coexistence equilibrium can be predicted. In Section 5 we fix the dilution rate D and we provide the operating diagram in the (S 0 , p 0 )-plane, see Fig.  2. In Sect. 6, we make precise our earlier description of the operating diagram, for the model (1) studied in [9,11], see the propositions 5, 6 and 7. The paper concludes with a discussion. Technical proofs and numerical computations are given in the Appendix.
2. Mathematical model . The model of the chemostat with a lethal external inhibitor we consider here is of the form The functional responses f 1 and f 2 represent the specific growth rates of the competitors and the function g represents the absorption rate of the external inhibitor by species y. The model (1) considered in [9], is obtained by letting in (3) In this paper, we consider the general model (3) without restricting ourselves to the special case (4). We suppose only that f 1 , f 2 , and g in (3) are C 1 -functions satisfying the following conditions: (H1) For i = 1, 2, f i (0) = 0 and f i (S) > 0 for all S > 0. (H2) g(0) = 0 and g (p) > 0 for all p > 0. Following [9], we could use a change of variables that reduces the operating parameters D, S 0 and p 0 and the yields coefficients β 1 and β 2 to 1. However, since we are interested in the role of the operating parameters, and since the yields coefficients have their own importance for the biologist, we do not make this reduction. It is easy to prove that for non-negative initial conditions, all solutions of system (3) are bounded and remain non-negative for all t > 0.
• The positive equilibrium E c = (λ 2 , x c , y c , p c ) of coexistence of the species, where p c , y c and x c are given by (11), (12) and (13), respectively. This equilibrium point exists iff λ − < λ 2 < S 0 and λ 2 < λ + and is LES iff, in addition, , and A 4 , are defined by (17). Proof. The proof is similar to the proof given in [9] and is left to the reader.
We observe that E 0 is LES iff E 1 and E 2 do not exist, and E c exists iff E 2 exists and is unstable, and E 1 , if it exists, is unstable. One concludes that there is only one equilibrium which can be stable. The results on existence and local stability of equilibrium points of (3), given by Prop. 1, are summarized in Table 1.

Existence
Local exponential stability Table 1. Existence and stability of equilibrium points E 0 , E 1 , E 2 and E c of (3), given in Prop. 1. Here, λ 2 and λ + are given by (5), λ − is given by (7) and A 1 , A 2 , A 3 and A 4 are given by (17). 3. Existence and stability of equilibrium points with respect to operating parameters. We assume that the biological parameters, that is β 1 , β 2 , γ and the functions f 1 and f 2 , are fixed. The Table 1 gives the conditions of existence and stability of the equilibrium points of (3) implicitly with respect to the operating parameters. We give now our main result which states these conditions more explicitly in the operating parameters. For this purpose, we need the following definitions. We let Note that F 1 is defined on I 2 × R + and F 2 is defined for (D, p 0 ) ∈ J c , where J c is given by (15). We define also the function where A 1 , A 2 , A 3 , and A 4 are defined by (17). Note that F 3 is defined for (D, p 0 , S 0 ) ∈ D c , where D c is given by (16). We have the following description of existence and stability of the equilibrium points of (3).
Proposition 2. Assume that the hypotheses and notations of Prop. 1 hold. The conditions of existence and stability of equilibrium points of (3) can be expressed with respect to the operating parameters D, p 0 and S 0 as in the Table 2, where F 1 (D, p 0 ), F 2 (D, p 0 ) and F 3 (D, p 0 , S 0 ) are defined by (19), (20) and (21), respectively.
Proof. The proof is given in Section A.1 Existence Local exponential stability Table 2. Existence and stability of equilibrium points of (3), with respect to the operating parameters D, S 0 and p 0 . The functions F 1 , F 2 , F 3 are defined by (19), (20), (21), respectively.
4. The region of instability of E c . We give now the necessary and sufficient condition on the operating parameters D, p 0 and S 0 such that the positive equilibrium E c is unstable, that is we discuss the sign of F 3 (D, p 0 , S 0 ), which is a polynomial of degree 3 in x c , thus we chose the following representation With the coefficients a i = a i (D, p 0 ), for i = 0..3. Hence, F 3 given by (22), appears as a third order polynomial in x c whose coefficients are depending only on D and p 0 and not on S 0 . We need to study the positive roots of F 3 . We have the following result. Lemma 1. Let ∆ = ∆(D, p 0 ) be the discriminant of the polynomial F 3 given by (22). Then, one of the four following exclusive cases may occur • F 3 has no positive root if {a 3 > 0 and ∆ < 0} or {a 3 > 0 and ∆ ≥ 0 and a 2 > 0 and a 1 a 2 > a 0 a 3 } (23) • F 3 has one positive (and only one) simple root if • F 3 has two (and only two) positive simple roots if • F 3 has three positive distinct roots if a 3 < 0 and ∆ > 0 and a 2 > 0 and a 1 a 2 < a 0 a 3 Proof. The proof is given in Section A.2 This lemma shows that the stability of E c , that is to say the determination of the sign of F 3 (D, p 0 , S 0 ) defined by (22), is reduced to the determination of the sign of the four real valued functions a 2 , a 3 , ∆ and a 1 a 2 − a 0 a 3 depending only on the operating parameters (D, p 0 ). The method used in the next section consists in fixing the biological parameters, then plotting the curves of the (D, p 0 )-plane where these functions vanish. These curves divide the subset J c of the (D, p 0 )-plane, defined by (15), into four regions, denoted Υ 0 , Υ 1 , Υ 2 and Υ 3 in each of which one and only one of the four cases described in Lemma 1 occurs: where J c is defined by (15). These regions, together with their boundaries, form a partition of J c . The result on the stability of E c can be reformulated as follows: , E c exists and is LES. • Proof. The proof is given in Section A.3 5. Operating diagram. Our aim now is to describe the operating diagram. The boundaries of the regions in the operating diagram are surfaces where bifurcations occur. In order to construct the operating diagram of the system one must compute these boundaries. Using Prop. 2, these boundaries are the following surfaces of the (D, S 0 , p 0 )-space. The subset is the border to which E 1 exists. The subset is the border to which E 2 exists. The subset is the border to which E 1 is stable. The subset is the border to which E 2 is stable. The subsets Γ 3 and Γ 4 are the border to which E c exists. The subset is the border to which E c is stable. It is not easy to represent the regions of existence and stability of the equilibrium points in the three dimensional space (D, p 0 , S 0 ). For this reason we will fix D and show the regions of existence and stability in the operating plane (p 0 , S 0 ).
The intersections of the surfaces Γ 1 , Γ 2 , Γ 3 and Γ 4 with the D-plane where D > 0 is fixed are the curves in the (p 0 , S 0 )-plane described in Table 3. The equations of these curves given in this table are equivalent to those given in (31), (32), (33), (34) and (35) and are obtained by using Lemma 2. In the (p 0 , S 0 )-plane, the curves Γ i , 5 and J S 6 , as illustrated by Fig. 2. In the case where the curve Γ 5 is not empty, two additional regions can appear, J U 5 and J U 6 , such that E c is LES in J S 5 and J S 6 and unstable in J U 5 and J U 6 . As a corollary of Prop. 2 we obtain the following result.
Proposition 4. The behavior of the system in each of the nine regions J 0 , J 1 , J 2 , J 3 , J 4 , J S 5 , J U 5 , J S 6 , and J U 6 is given in Table 4. Some of the regions may be empty as shown on Fig. 2: • Fig. 2(a) corresponds to the case without detoxification, that is g(p) = 0. In this case the line Γ 4 is vertical, so that the regions J S 5 , J S 6 , J U 5 and J U 6 are empty. Therefore, the positive equilibrium E c does not exist. The species cannot coexist. • Fig. 2(b) corresponds to the case where Γ 5 is empty and the tangent of Γ 1 at point (p c , λ 2 ) is under Γ 4 , so that the region J 4 exists. • Fig. 2(c) corresponds to the case where Γ 5 is not empty and the tangent of Γ 1 at point (p c , λ 2 ) is above Γ 4 , so that the region J 4 disappeared, while the regions J U 5 and J U 6 where E c is unstable appeared. Figure 2. Illustrative operating diagrams for D fixed: The curves Γ i , i = 0 · · · 5 defined in the Table 3 divide the operating plane . Some of the regions may be empty. The existence and stability of the equilibrium points in the regions of these diagrams are shown in Table 4.
The operating diagrams shown in Fig. 2 are given only as illustrative examples, showing that our analysis gives a complete description of the behavior of the system. Notice that for plotting operating diagrams we must choose functions f 1 , f 2 , and g in (3) and fix the values of the biological parameters. We illustrate this in the following section for the model considered in [11,9]. Our aim is to show that the plot of the regions Υ i , defined in Section 4, is a useful tool for the understanding of the asymptotic behavior of the model with respect to the operating parameters.
6. Model of Hsu, Li and Waltman [9]. In this section, we consider the model (1), studied by [9]. We illustrate our results for the five sets of biological parameters given in Table 5. The sets of biological parameters considered in lines 1, 2 and 3 of this table were proposed in [9]. We consider also those in lines 4 and 5, since they are exhibiting interesting behaviors. For these sets of biological parameters, the interval I c defined by (14), for which p c (D) exists, is of the form I c = 0, D with D corresponding to the intersection f 1 (S) = f 2 (S). This numerical value which depends only on m 1 , m 2 , K 1 and K 2 is given in the last column of the table.
6.1. Biological parameters in Table 5, Case 1. For these values of the biological parameters, it was shown in [9] that E c exists and is LES for S 0 = p 0 = D = 1, see Fig. 5 in [9] and the text preceding this figure. Using our results on the stability of E c , we show that, actually, for all values of the operating parameters, not only for S 0 = p 0 = D = 1, whenever it exists E c is LES. The proof is based on a numerical computation of the regions Υ i , leading to the following claim: Claim 1. Let the biological parameters be given by Table 5, Case 1. Then J c = Υ 0 .
Proof. Numerical computations show that a 2 > 0, a 3 > 0 and a 1 a 2 > a 0 a 3 for all (D, p 0 ) ∈ J c , where J c is defined by (15). Therefore, only case (23) can occur, so that J c = Υ 0 .  Table 5, Case 1. J c = Υ 0 : therefore E c is LES whenever it exists. The operating diagram for D = 1 shows that p 0 = 1, S 0 = 1 ∈ J S 6 . The existence and stability of the equilibrium points in the regions of this diagram are shown in Table 4 Proposition 5. Let the biological parameters be given by Table 5, Case 1. For all values of the operating parameters, E c is LES, whenever it exists, that is, for all , E c exists and is LES.
Proof. From Claim 1 we have. J c = Υ 0 . From Prop. 3 we deduce that E c is LES whenever it exists.
As an illustration we show in Fig. 3 the operating diagram corresponding to D = 1. This diagram shows that (p 0 = 1, S 0 = 1) ∈ J S 6 , i.e. E c exists and is LES for S 0 = p 0 = D = 1, which agrees with the result of [9].   Table 4. Table 5, Case 2. For these values of the biological parameters, it was shown in [9] that E c exists and is LES for S 0 = p 0 = D = 1, see Figs. 1 and 2 in [9] and the text preceding thess figures. Using Prop. 3, we can show that actually, for D = 1 and all values of the two other operating parameters, not only for S 0 = p 0 = 1, E c is LES whenever it exists. Moreover, there exists a set of values for the operating parameters, such that E c exists and is unstable. The proof is based on a numerical computation of the regions Υ i , leading to the following claim: Claim 2. Let the biological parameters be given by Table 5, Case 2. Then J c = Υ 0 ∪ C ∪ Υ 2 , as shown in Fig. 4, where the blue curve, C, is the Jordan closed component of the curve ∆ = 0.

Biological parameters in
Proof. The proof is given in Section B.1.
Using this claim, which was proved by numerical calculations, we can obtain the following result. Proposition 6. Let the biological parameters be given by Table 5, Case 2. There exist D 1 and D 2 , and p 1 (D) < p 2 (D) defined for D 1 < D < D 2 , such that If D ≤ D 1 or D ≥ D 2 then E c is LES whenever it exists, that is, for all p 0 > p c (D) and all S 0 > λ 2 (D) + yc(D,p 0 ) β2 , E c exists and is LES. If D 1 < D < D 2 then the region of instability of E c in the operating (S 0 , p 0 )-plane is a bounded domain of the form where S 0 i (D, p 0 ), i = 1, 2 are defined in Prop. 3. Proof. The Figure 4 exhibits two bifurcation values D 1 1.83 and D 2 2.65 such that for D 1 ≤ D ≤ D 2 , the boundary between Υ 0 and Υ 2 is defined by p 0 = p 1 (D) As an illustration we show in Fig. 5(a) the operating diagram corresponding to D = 1. This diagram shows that (p 0 = 1, S 0 = 1) ∈ J S 6 , i.e. E c exists and is LES for S 0 = p 0 = D = 1, which agrees with the result obtained in [9]. However, the operating diagram in Fig. 5(a) shows that, for D = 1, and for all operating parameters S 0 and p 0 , not only for S 0 = p 0 = 1, E c is LES, whenever it exists. The operating diagram corresponding to D = 2.2 is shown in Fig. 5(b), illustrating for this set of biological parameters the possibility of instability of E c on J U 6 of the form (36). The particular values p 1 = p 1 (2.2) and p 2 = p 2 (2.2), shown on the p 0 -axis of the figure correspond to the values depicted in Fig. 4. Table 5, Case 3. For these values of the biological parameters, it was shown in [9] that E c exists and is unstable for S 0 = p 0 = D = 1, see Figs. 3 and 4 in [9] and the text preceding these figures. Notice that the parameters are as in Case 2 of the table except that K 1 = 0.03.

Biological parameters in
Using Prop. 3, we can show that actually, E c exists and is unstable for an unbouneded set of values of the operating parameters, not only for S 0 = p 0 = D = 1. The proof is based on a numerical computation of the regions Υ i , leading to the following claim: Claim 3. Let the biological parameters be given by Table 5, Case 3. Then J c = Υ 0 ∪ C 1 ∪ Υ 1 ∪ C 2 ∪ Υ 2 , as shown in Fig. 6, where the blue curve, C 1 , is the Jordan closed component of the curve ∆ = 0 and the red curve, C 2 , is the curve a 3 = 0.
Proof. The proof is given in Section B.2.
Using this claim, which was proved by numerical calculations, we can obtain the following result.  Table 4. Proposition 7. Let the biological parameters be given by Table 5 , then (D, p 0 ) ∈ Υ 1 , and hence, according to Prop. 3, E c is unstable iff S 0 > S 0 1 (D, p 0 ). As an illustration we show in Fig. 7 the operating diagram corresponding to D = 1. Note that D 1 < 1 < D 3 , as shown in Fig. 6(a). Therefore, the operating diagram, for D = 1 exhibits a compact region of instability of E c which is of the form J U 6 as shown in Fig. 7(a). A zoom in Fig. 7(b) shows that (p 0 = 1, S 0 = 1) ∈ J U 6 , i.e. E c is unstable for S 0 = p 0 = D = 1, which agrees with the result obtained in [9]. However, the operating diagram in Fig. 7(a) shows that, for D = 1, and for a large set of operating parameters S 0 and p 0 , not only for S 0 = p 0 = 1, E c is unstable, whenever it exists.
We have claimed at the beginning of this subsection that the instability of E c occurs for an unbounded set of values of the biological parameters. This can be seen from the follwing observation. If D 3 < D < D 4 , then J U 6 becomes unbounded, as predicted by Prop. 7. See the case D = 2.2 in Fig. 8. The particular values p 1 = p 1 (2.2), p 2 = p 2 (2.2), p 3 = p 3 (2.2) and p 4 = p 4 (2.2) shown on the p 0 -axis of the figure correspond to the values depicted in Fig. 6. As predicted by Prop. 7 the vertical lines p 0 = p 3 and p 0 = p 4 correspond to asymptotes of the unbounded branches of the boundary Γ 5 (in green on the figure) of the region J U 6 . The numerical plots of the operating diagrams presented in Fig. 5(a), Fig. 7 and Fig. 8 show that the region of instability of E c , if it exists, is always of the form J U 6 , that is to say, the region J U 5 , depicted in Fig. 2(c), never existed. In the next section we give a set of biological parameter values for which the region J U 5 is non empty. Table 5, Case 4. We see in this section that for a suitable choice of biological parameters, it is possible that J U 5 is not empty. This situation occurs when curves Γ 5 and Γ 1 intersect: both regions of instability J U 5 and J U 6 can exist. Let us consider the parameter values given in Table 5, Case 4. Figure. Figure 10. The biological parameters values are given in Table  5 Table 4.

Biological parameters in
result similar to Prop. 7. We do not give the general result. We emphasize only the following two interesting cases. If D = 0.01, there exist p 1 and p 2 , such that, if p 1 < p 0 < p 2 , then (D, p 0 ) ∈ Υ 2 , see Fig. 9(a). As a consequence of Prop. 3, the region of instability is bounded.
On the other hand for D 1 = 0.013, there exist p 1 , p 2 , p 3 , p 4 , such that, if p 1 < p 0 < p 3 or p 4 < p 0 < p 2 , then (D, p 0 ) ∈ Υ 2 , and if p 3 < p 0 < p 4 then (D, p 0 ) ∈ Υ 1 , see Fig. 9(b). As a consequence of Prop. 3, the region of instability is unbounded.  Table 4. Fig. 10, shows that in both cases, Γ 5 and Γ 1 intersect, so that J U 5 and J U 6 are non empty. For D = 0.001, both J U 5 and J U 6 are bounded, see Fig. 10(a). For D = 0.001, J U 5 is unbounded and J U 6 is bounded, see Fig. 10(b). 6.5. Biological parameters in Table 5, Case 5. When the tangent of Γ 1 at (p c , λ 2 ) is above Γ 4 , the region J 4 exists, see Fig. 2(b). It is shown in [4] that, in this case, for the model (2) where the inhibitor is not lethal, E c is LES whenever it exists. The aim of this section is to show that this property does not hold for the model (1) with a lethal inhibitor. Let us consider the parameter values given in Table 5, Case 5. Fig. 11 shows the corresponding regions Υ i . We emphasize on (a) Figure 13. The biological parameters values are given in Table 5,  Table 4.
the following interesting case. For D = 1.26, there exist p 1 and p 2 , such that, if p 1 < p 0 < p 2 , then (D, p 0 ) ∈ Υ 2 , see Fig. 11(b). As a consequence of Prop. 3, the region of instability J U is bounded, see the region J U 6 in Fig. 12(a). This is a very interesting case where the tangent of Γ 1 at (p c , λ 2 ) is above Γ 4 as shown in the zoom proving that the region J 4 exists, see Fig. 12 This case is also interesting since it depicts a very narrow interval D 1 < D < D 2 , with D 1 1.242 and D 2 1.274, such that the region J U 6 of instability of E c exists and is bounded. This shows that our analysis in the (D, p 0 )-plane and the constructions of the regions Υ i in this plane is a very useful tool for the understanding of the model. A simple numerical exploration of the domain of operating parameters would have not revealed easily this behavior of the model. 6.6. Operating diagrams with respect to S 0 and D . Instead of fixing D and showing the operating diagram in the p 0 , S 0 -plane we can fix p 0 and show the operating diagram in the S 0 , D -plane. This is illustrated in Fig. 13(a) where the biological parameters are given in Case 3 of Table 5 and p 0 = 1. This figure shows that (S 0 = 1, D = 1) ∈ J U 6 which agrees with the result of [9]. Notice that as predicted by the plot of the regions Υ i given in Fig. 6, the region of instability is unbounded, since the horizontal line p 0 = 1 intersects the region Υ 1 in Fig. 6. 6.7. Operating diagrams with respect to p 0 and D. We can fix also S 0 and show the operating diagram in the p 0 , D -plane. This is illustrated in Fig. 14 where the biological parameters are given in Case 3 of Table 5 and S 0 = 1. This figure shows that (p 0 = 1, D = 1) ∈ J U 6 which means that E c exits and is unstable for p 0 = S 0 = D = 1.
As it was said in the introduction, Hsu, Li and Waltman [9] constructed an operating diagram. They fixed S 0 and presented the operating diagram in the p 0 , D -plane, see Fig. 6 in [9]. To show the role of the yields, these authors have chosen the values β 1 = β 2 = 100. In Fig. 15, we use also β 1 = β 2 = 100, and we present the the operating diagram, in the p 0 , D -plane, for S 0 = 1 and the biological parameters given in table 5, Case 3. Our diagram is qualitatively similar Figure 14. The biological parameters values are given in Table  5, Case 3. The operating diagram in the (p 0 , D)-plane for S 0 = 1, showing the instability of E c for S 0 = p 0 = D = 1. The existence and stability of the equilibrium points in the regions of this diagram are shown in Table 4.. Figure 15. The biological parameters values are given in Table 5 Table 4.
to Figure 6 in [9], but has significant quantitative differences, probably due to some differences between our values of biological parameters and those used in Fig. 6 of [9]. Actually, these authors did not give precisely the values of the biological parameters nor that of the operating parameter S 0 , they used in their plot.
In [9] five regions where identified, labeled I = J 0 , corresponding to the stability of E 0 , II = J 1 ∪ J 3 , corresponding to the stability of E 1 , III = J 2 , corresponding to the stability of E 2 , IVa = J U 6 , corresponding to the instability of E c , that is the limiting behavior is oscillatory, and IVb = J S 5 ∪ J S 6 , corresponding to the stability of E c . Maximal growth rates of the competitors 1/time Half saturation constants of the competitors mass/volume δ Maximal growth rate of detoxification 1/time K Half saturation constant of detoxification mass/volume γ Lethal effect of p on x volume/mass β 1 , β 2 Growth yield coefficients dimensionless Table 6. Meanings and units of the variables and parameters of (1).
7. Discussion. In this work we have extended the model (1) of competition in the chemostat with an external lethal inhibitor studied in [9] by considering the model (3) with general growth rate functions of competitors and absorption rate of external inhibitor. Our mathematical analysis of the model has revealed several possible behaviors : Prop. 2 provides a complete theoretical description of the outcome of competition. In order that the results can be useful in practice, one must have a description of the outcome of competition with respect to the operating parameters. The study of bifurcations according to the operating parameters p 0 , D and S 0 is the most meaningful one for the laboratory model, since the experimenter can easily vary these parameters. For instance, the outcome of competition between two species of microorganisms in the chemostat with an external lethal inhibitor has been shown to be quite different from that in a chemostat without an inhibitor, where only competitive exclusion occur. However, the coexistence is possible only if the inhibitor is detoxified from the chemostat. Without detoxification, that is if g = 0, coexistence is not possible as shown in Fig. 2(a).
The previous work [9] on (1) did not addressed the question of the robustness of the results with respect to parameter variations. The robustness is the principal motivation for the present work. Some of the highlights of this paper are the results on the numerical calculations of the regions Υ 0 , Υ 1 and Υ 2 of the operating plane (D, p 0 ), defined by (27), (28) and (29) respectively, and contained in Figs. 4, 6, 9 and 11. We were not able to find a set of biological parameters for (4) such that the region Υ 3 , defined by (30), is not empty. Whether or not the region Υ 3 is always empty for (4) is an interesting question and deserves further attention.
The regions Υ i , together with their boundary curves, are a partition of the first quadrant of (D, p 0 )-plane in parameter regions in which the region of instability is empty, bounded or unbounded. While explicit formulas for the boundary curves are known (see Sect. 4), they are too complicated to be analyzed, yet easy to plot, once the biological parameters are fixed. Figs. 4, 6, 9 and 11 correspond to different sets of values for the biological parameters and illustrate the dramatic effect of a change of the biological parameters on the operating diagram.
Once the biological parameter are fixed, the figure showing the regions Υ i , predicts what should be the outcome of competition as the operating parameters S 0 and p 0 are varied and D is fixed, or S 0 and D are varied and p 0 is fixed. For instance, if the biological parameters values are as given in Table 5, Case 3, then the Fig. 6 shows that if D = 2.2, then the region of instability J U 6 is unbounded, since the vertical line D = 2.2 in Fig. 6 intersects the region Υ 1 . This region is plotted in Fig. 8. Similarly, if p 0 = 1, then the region of instability J U 6 is unbounded, since the horizontal line p 0 = 1 in Fig. 6 intersects the region Υ 1 . This region is plotted in Fig. 13. Therefore, the plots of the operating diagrams depicted in Figs Our propositions 5, 6 and 7, although based on Claims 1, 2 and 3, respectively, whose proofs are obtained by numerical computations of the regions Υ i , predict the outcome of competition in a rigorous way. These propositions extend the results of [9] since they show that: • For the biological parameter values given in Table 5, Case 1, whenever it exists, E c is LES, for all values of the operating parameters, extending the result of [9] obtained for S 0 = p 0 = D = 1.
• For the biological parameter values given in Table 5, Case 2, whenever it exists, E c is LES, for all D ≤ D 1 ≈ 1.83 or D ≥ D 2 ≈ 2.65 and all values of the operating parameters S 0 and p 0 , extending the result of [9] obtained for S 0 = p 0 = D = 1. Moreover, for D 1 < D < D 2 , there exists a bounded set of operating parameters p 0 and S 0 such that E c is unstable whenever it exists.
• For the biological parameter values given in Table 5, Case 3, whenever it exists, E c is unstable, for an unbounded set of values of the operating parameters, extending the result of [9] obtained for S 0 = p 0 = D = 1.
The curves Γ i in the (p 0 , S 0 )-plane (with D fixed) depicted in Figs. 3, 5, 7, 8, 10, 12, divide the first quadrant into parameter regions in which washout occurs, or only one species can survive by itself, or coexistence occurs at a stable or an unstable equilibrium. These curves have explicit equations (see Table 3), except for Γ 5 that can be easily plotted once one of the operating parameter D is fixed and the biological parameters are chosen. For instance, Fig. 7 describes the parameter regions for D = 1 and the biological parameters fixed as in [9] and Proof. Using (5), (H1) and (19), we have On the other hand, using (5), (11) and (H1), we have This completes the proof of the first equivalence in the lemma. Using (5), (7), (11) and (H1), we have Using the fact that ∂W ∂p < 0 and (6), we have This completes the proof of the second equivalence in the lemma. Now we give the proof of Prop. 2. Using (5) and hypothesis (H1), the conditions λ + > S 0 and λ 2 > S 0 of stability of E 0 in Table 1 is equivalent to D > f 1 (S 0 ) − γp 0 and D > f 2 (S 0 ). Similarly, the condition λ + < S 0 of existence of E 1 in Table 1 is equivalent to D < f 1 (S 0 ) − γp 0 . On the other hand, the condition λ 2 < S 0 of existence of E 2 in Table 1 is equivalent to D < f 2 (S 0 ).
Let us consider now the stability of E 1 and E 2 . Using Lemma 2, the condition λ + < λ 2 of stability of E 1 in Table 1 is equivalent to D < F 1 (D, p 0 ) and the condition λ 2 < λ − of stability of E 2 in Table 1 is equivalent to S 0 < F 2 (D, p 0 ).
Let us consider now the existence and stability of E c . Using Lemma 2, we see now that the condition λ − < λ 2 < λ + of existence of E c in Table 1 is equivalent to D > F 1 (D, p 0 ) and S 0 > F 2 (D, p 0 ). Finally, using the definition (21) of the function F 3 , the condition of stability of E c in Table 1 is equivalent to F 3 (S 0 , D, p 0 ) > 0.
A.2. Proof of Lemma 1. We have a 0 = y c bD(ay c + D)(by c + D)A 0 , where A 0 = a(a + b)y 2 c + D(2b + 3a)y c + 2D 2 Thus a 0 > 0. Therefore F 3 (0) = a 0 > 0. If a 3 = 0, then the third order polynomial has three distinct real roots, when ∆ > 0 or one real root and two complex conjugate roots, when ∆ < 0, where ∆ is the discriminant of the polynomial. Recall the Routh-Hurwitz conditions • The roots of polynomial F 3 (x) are all of negative real part iff a 2 > 0 and a 1 a 2 > a 0 a 3 . • The roots of polynomial F 3 (−x) are all of negative real part iff a 2 > 0 and a 1 a 2 < a 0 a 3 . Assume that (23) holds. Since a 3 > 0, F 3 (−∞) = −∞, so that, using F 3 (0) > 0, the polynomial F 3 has at least one negative root. If ∆ < 0, F 3 (x) has no other real root, so that it has no positive root. If ∆ ≥ 0, and the Routh-Hurwitz conditions for F 3 (x) hold, that is to say a 2 > 0 and a 1 a 2 > a 0 a 3 , all the roots have negative real parts. Hence F 3 (x) has no positive root. Therefore, if (23) holds, the polynomial F 3 has no positive root.
Assume that (24) holds. Since a 3 < 0, F 3 (+∞) = −∞, so that, using F 3 (0) > 0, the polynomial F 3 has at least one positive root. If ∆ < 0 then the polynomial F 3 has no other real root. If ∆ ≥ 0, and the Routh-Hurwitz conditions do not hold for F 3 (−x), that is to say a 2 ≤ 0 or a 1 a 2 ≥ a 0 a 3 , there is no other real root. Therefore, if (24) holds, the polynomial F 3 has one and only one positive simple root.
Assume that (25) holds. Since a 3 > 0, F 3 (−∞) = −∞, so that, using F 3 (0) > 0, the polynomial F 3 has at least one negative root. If ∆ > 0, and the Routh-Hurwitz conditions for F 3 (−x) do not hold, that is to say a 2 ≤ 0 or a 1 a 2 ≤ a 0 a 3 , the roots are real, distinct, and two of them are necessarily positive. Therefore, if (25) holds, the polynomial F 3 has two and only two positive roots.
Assume that (26) holds. Since a 3 < 0 then F 3 (+∞) = −∞, so that, using F 3 (0) > 0, the polynomial F 3 has at least one positive root. If ∆ > 0, and the Routh-Hurwitz conditions hold for F 3 (−x), that is to say a 2 > 0 and a 1 a 2 < a 0 a 3 , the roots are distinct, real and positive. Therefore, if (26) holds, the polynomial F 3 has three positive distinct roots.
Appendix B. Construction of the figures 4 and 6. In this section we show how to compute numerically the regions Υ i . Let C be a Jordan curve included in J c . We denote by Int (C), the interior of the Jordan curve C, and by Ext (C 0 ) the intersection of its exterior with the domain J c . Figure 16. The plots of curves p 0 = p c (D), in black, ∆ = 0, in blue, and a 1 a 2 = a 0 a 3 in magenta. The biological parameters values are given in Table 5, Case 2.

BACHIR BAR AND TEWFIK SARI
B.1. Proof of Claim 2. Let the biological parameters be given by Table 5, Case 2. The domain J c , defined by (15) is the region of the (D, p 0 )-plane located above the curve p 0 = p c (D), where 0 ≤ D ≤ D, colored black on Fig. 16. Numerical computations show that a 2 > 0 and a 3 > 0 for all (D, p 0 ) ∈ J c . Therefore we have to consider only the signs of ∆ and a 1 a 2 − a 0 a 3 . The numerical calculation of the curve a 1 a 2 = a 0 a 3 shows that it is a Jordan curve, labeled C 0 and plotted in magenta in Fig. 16. One has a 1 a 2 < a 0 a 3 in Int (C 0 ), and a 1 a 2 > a 0 a 3 on Ext (C 0 ). The numerical calculation of the curve ∆ shows that it has three components, plotted in blue in Fig. 16: a Jordan curve labeled C and located in Int (C 0 ), and two arcs, intersecting at the origin, labeled A 1 and A 2 , and located in Ext (C 0 ). One has ∆ > 0 in Int (C) and ∆ < 0 in Ext (C), excepted in the regions labeled R 1 and R 2 in Fig. 16, where ∆ > 0, and on the curves A 1 and A 2 , where ∆ = 0. On Int (C) one has a 3 > 0 and ∆ > 0 and a 1 a 2 < a 0 a 3 (37) From (37) it is deduced that (25) is satisfied on Int (C). Hence Int (C) ⊂ Υ 2 , where Υ 2 is defined by (29).
Notice that Ext (C) = A ∪ (Ext (C) \ A) where A = R 1 ∪ A 1 ∪ R 2 ∪ A 2 . On A, which is included in Ext (C 0 ), one has a 3 > 0 and ∆ ≥ 0 and a 2 > 0 and a 1 a 2 > a 0 a 3 On Ext (C) \ A, one has a 3 > 0 and ∆ < 0 (39) From (38) and (39) it is deduced that (23) is satisfied on Ext (C). Hence Ext (C) ⊂ Υ 0 , where Υ 0 is defined by (27). Therefore Υ 0 = Ext (C) and Υ 2 = Int (C), as it is claimed in Fig. 4. B.2. Proof of Claim 3. Let the biological parameters be given by Table 5 The plots of curves a 1 a 2 = a 0 a 3 (C 3 ∪ C 4 , in magenta) and a 2 = 0 (C 5 , in cyan). The biological parameters are given in Table  5, Case 3. The plots of curves ∆ = 0 (C 1 ∪ A 1 ∪ A 2 , in blue), a 3 = 0 (C 2 , in red), a 1 a 2 = a 0 a 3 (C 3 ∪ C 4 , in magenta) and a 2 = 0 (C 5 , in cyan). (b): A zoom of the strip 0 < p 0 < 1. The biological parameters are given in Table 5, Case 3. the numerical computations show that all functions a 2 , a 3 , a 1 a 2 − a 0 a 3 and ∆ can change their signs. The signs of these functions are shown On Fig. 17 and summarized on Table 7.