FEEDBACK-MEDIATED COEXISTENCE AND OSCILLATIONS IN THE CHEMOSTAT

We consider a mathematical model that describes the competition of three species for a single nutrient in a chemostat in which the dilution rate is assumed to be controllable by means of state dependent feedback. We consider feedback schedules that are affine functions of the species concentrations. In case of two species, we show that the system may undergo a Hopf bifurcation and oscillatory behavior may be induced by appropriately choosing the coefficients of the feedback function. When the growth of the species obeys Michaelis-Menten kinetics, we show that the Hopf bifurcation is supercritical in the relevant parameter region, and the bifurcating periodic solutions for two species are always stable. Finally, we show that by adding a third species to the system, the two-species stable periodic solutions may bifurcate into the coexistence region via a transcritical bifurcation. We give conditions under which the bifurcating orbit is locally asymptotically stable.


Introduction.
We study the mathematical model of a chemostat: where S is the nutrient concentration and x i are the concentrations of the competing species. The f i are the growth functions of the species and they are assumed to be monotonically increasing. The yield constants, denoted by γ i , reflect that only a fraction of the nutrient of what the different species consume, leads to new biomass.
The two natural control parameters are the input nutrient concentration S 0 and the dilution rate D. If these are constant, then at most one species survives as dictated by the principle of competitive exclusion [17,5,29,25]. But in practice, it often happens that S 0 and/or D change over time, and a well-established literature exists in case these variations are periodic. For results in the case S 0 is periodic, see [23,13,22], and for periodic variations in D see [4,25]. In both cases coexistence is possible under certain conditions, contrasting the competitive exclusion principle. The case when S 0 is a general time-dependent function has been studied in [11].
More recently a program of feedback control in the chemostat has been initiated in [9], where the dilution rate is treated as a feedback variable and made dependent on the concentrations of the competitors in the following simple way: where ǫ and the k i are non-negative design parameters. The feedback approach is perhaps most natural in the lab setting. For instance, optical sensors can be used to measure turbidity, giving a rough estimate of the concentrations of the species. An alternative way to measure concentrations is to use GFP's (Green Fluorescence Proteins), especially since nowadays GFP's can emit light of different color, thus allowing to distinguish between different species. These concentration estimates can be processed by a computer to (online) calculate the dilution rate. The result then determines the speed of the pump -the device that is being actuated-which supplies the reactor with fresh medium.
The introduction of feedback controls can lead to significant changes in the asymptotic behavior of solutions when compared to the case in which the control parameters are fixed. These changes have been well documented for the case n = 2. For example, it was shown in [9] that the system may be rendered coexistent for suitable choices of the feedback parameters. This result was generalized in [12] to include chemostat models with non-monotone growth functions. In both cases, coexistence takes its simplest form: a globally asymptotically stable positive steady state. For other choices of the feedback parameters -but still assuming n = 2 -it was shown in [10] that bistability can occur, meaning that the single-species boundary steady states are stable and attract almost all solutions starting in the interior, with the exception of solutions on the stable manifold of a saddle point in the interior.
Unfortunately, it was pointed out in [9] that feedback-mediated-coexistence using the feedback law above typically does not occur if n > 2. In fact, it was shown in [10] that if there are three competitors, no interior solution is persistent and at least one species disappears.
Here we wish to investigate another type of feedback law for the dilution rate, and see whether or not coexistence is possible when there are more than two competitors. We will show that this is indeed the case when n = 3, using the following feedback law: for certain values of the non-negative parameters ǫ and k i . Notice the subtle difference between the latter and former feedback law, expressed merely by the sign change, and that the concentration of the third species is not used in the feedback law. Now, coexistence will take the form of an interior periodic solution, bifurcating from a periodic solution on the boundary through what is traditionally referred to as a transcritical bifurcation of the associated Poincaré map. Transcritical bifurcations of limit cycles leading to coexistence have been reported in [2], [3] (see also Chapter 3 in [25]), and in [20,1] for different chemostat models. In general, it is not FEEDBACK-MEDIATED OSCILLATORY COEXISTENCE IN THE CHEMOSTAT 3 difficult to establish the existence of the transcritical bifurcation. It considerably more difficult to determine the stability of the bifurcating periodic orbit because the precise location of the boundary periodic solution is typically unknown. In the present case, the stability determination is based upon a perturbation technique used by Smith [24] in the study of a system with two predators and one prey. For a chemostat with three species, is shown that local asymptotically stable periodic orbits result when the growth functions intersect in a specified manner.
As we already alluded to in the last paragraph, we will show along the way that in a chemostat with n = 2 species and controlled by a feedback law of the new type, asymptotically stable periodic solutions are possible (these will be the boundary periodic solutions mentioned in the previous paragraph). Perhaps this result is of interest by itself because this type of dynamical behavior does not occur in chemostats with 2 species which are controlled by the former feedback law, see [10] which will be briefly reviewed in the next section for convenient reference. Finally, this report also investigates the relationship between the divergence criteria [21] and the vague attractor property [19] in determining the stability of a Hopf bifurcation. Specifically, it is shown that the sign of the quadratic form for the divergence coincides with the sign of V ′′′ (0). This result is used in the study of coexistence for three organisms.

2.
Coexistence and bistability for two species. Consider the following chemostat which is an appropriately scaled version of the model presented in the Introduction:Ṡ where (S, x, y) ∈ R 3 + . It is assumed that f and g are nonnegative C 1 functions on R + with f (0) = g(0) = 0 and f ′ , g ′ > 0. Moreover we assume the existence of λ ∈ (0, 1) and D * > 0 such that f (λ) = g(λ) = D * , f ′ (λ) < g ′ (λ) and f (S) > g(S) for S ∈ (0, λ) and f (S) < g(S) for S > λ. We assume that D(x, y) is affine: where k 1 , k 2 and ǫ are nonnegative parameters. Notice that (1)-(2) is well-posed since R 3 + is a forward invariant set. The following coexistence result was obtained in [9].
Suppose that the k i satisfy k 2 <k < k 1 , and k 2 > 0 in case ǫ = 0.
Then system (1) with feedback (2), has a unique steady state , which is globally asymptotically stable with respect to initial conditions in int(R 3 + ). We wish to complete the picture of the dynamical behavior for all ranges of the parameters k i and start with the most interesting case of bistability, see [10]. Theorem 2. Fix ǫ ∈ [0, D * ) and letk := D * −ǫ 1−λ . Suppose that the k i satisfy k 1 <k < k 2 , and k 1 > 0 in case ǫ > 0. Then system (1) with feedback (2), has a unique steady state , which is a saddle. In addition to the washout steady state at E 0 = (1, 0, 0), there are two asymptotically stable single-species boundary steady states Almost all solutions starting in int(R 3 + ) converge to either E 1 or E 2 , except for those solutions starting on the two-dimensional (and hence zero-measure) stable manifold of E.
To complete the global picture we indicate the parameter region for the k i which gives rise to competitive exclusion, see [10] as well.
The results of Theorems 1, 2 and 3 are summarized in Figure 1.
The function D(x, y) is a positive C 3 function to be determined later. The model is well-posed since for V := S+x+y, we have thatV = D(x, y)(1−V ), and so T is forward invariant. Moreover, an obvious reduction argument suggests the study of the following two-dimensional system: where (x, y) ∈ ∆ = {(x, y) ∈ R 2 + | x + y ≤ 1}. Denoting the right-hand side of system (4) by (F (x, y), G(x, y)) and suppressing the argument of D and its partial derivatives, we find that the Jacobian matrix is: For future reference we also calculate the second order derivatives of F and G: and some of their third order derivatives: From now on we will assume that D(x, y) is affine: where k 1 , k 2 and ǫ are non-negative parameters to be determined later. An interior steady state of system (4)-(5) is a solution (x * , y * ) ∈ int(∆) of: Assuming that an interior steady state (x * , y * ) exists, we evaluate the trace and determinant of the variational matrix J(x * , y * ): Notice that if we set k 1 =k 1 and k 2 =k 2 where f ′ (λ) =:k 1 <k 2 := g ′ (λ) (7) in the feedback (5), then tr(J(x * , y * )) = 0 and det(J(x * , y * )) = x * y * (k 2 −k 1 ) 2 > 0, and hence the occurrence of a Hopf bifurcation becomes plausible. This suggests treating (k 1 , k 2 ) as a bifurcation parameter, while using ǫ > 0 to guarantee that 1. (x * , y * ) ∈ int(∆). 2. D(x, y) > 0 in ∆. Explicitly solving (6) for (x * , y * ) shows that and so the first condition is satisfied if a nonempty interval sincek 1 <k 2 . The affine function D(x, y) = −k 1 x −k 2 y + ǫ reaches its minimum in ∆ at the point (0, 1), and so the second condition holds if Notice that both constraints (9) and (10) for ǫ are compatible only if or equivalently,k 2 < D * /λ. The latter inequality is satisfied: by the mean value theorem there is some c ∈ (0, λ) such that g ′ (c) = D * /λ and since g ′ is decreasing as g ′′ < 0, this implies thatk 2 ≡ g ′ (λ) < g ′ (c). Therefore, if we fix an ǫ in the nonempty interval then there is some open neighborhood N of (k 1 ,k 2 ), such that for all (k 1 , k 2 ) ∈ N , D(x, y) > 0 in ∆, and system (4)-(5) has an interior steady state (x * , y * ) (this is because both D(x, y) and the point (x * , y * ) depend continuously on (k 1 , k 2 )). Moreover, the interior steady state (x * , y * ) undergoes a Hopf bifurcation at the bifurcation value (k 1 ,k 2 ). To be more precise, a Hopf bifurcation occurs along any smooth parametric path (k 1 , k 2 ) = (k 1 (σ), k 2 (σ)) such that (k 1 ,k 2 ) = (k 1 (0), k 2 (0)) and x * k ′ 1 (0) + y * k ′ 2 (0) = 0. The first condition ensures that the eigenvalues of J(x * , y * ) are purely imaginary at σ = 0. The second condition implies that this pair of complex conjugate eigenvalues of J(x * , y * ) crosses the imaginary axis transversally. Hence, a Hopf bifurcation occurs at σ = 0. In what follows, we will make use of one particular parametric path (k 1 , k 2 ) = (k 1 + σ,k 2 ) where the value of k 2 is fixed, and k 1 alone is treated as the bifurcation parameter. Such a path clearly satisfies the above conditions for the Hopf bifurcation. To determine the nature of the bifurcation (super-or subcritical), we follow the procedure outlined in [21]. There it is shown that if the sign of the quantity is negative (positive), then a supercritical (subcritical) Hopf bifurcation occurs. Here Q yy = (F * xyy + G * yyy ) + a * x F * yy + a * y (2F * xy + 3G * yy ) (there is no need to calculate Q xy since it is multiplied by F * x in (12), and the latter is 0 in our case). The quantities a * x and a * y are given by: We start by evaluating the first, second and third order derivatives mentioned above at the interior steady state (x * , y * ), using the particular form of the feedback (5) and the choices fork 1 andk 2 , given in (7) (we suppress the argument λ of the derivatives of the functions f and g): Substituting these values into the formulas for a * x and a * y , we obtain Substituting the quantities to see that and Adding the last two equalities, we obtain where we used the fact that x * + y * = 1 − λ. This can be simplified to Notice that the first and last term in the above sum cancel each other, so we conclude that The first term is negative since f ′′′ , g ′′′ < 0 by assumption. We see that if then the second term is negative as well, and hence a supercritical Hopf bifurcation occurs.
If both growth functions are of Michaelis-Menten form, we have a stronger result.
Lemma 1. Let f (s) = m1s a1+s and g(s) = m2s a2+s . Suppose that there exists λ ∈ (0, 1) and Using the relations 1 ai+λ = D * miλ and γ = m2 m1 , we find that We also find that , which can be further simplified substituting the expressions a i = λ( mi D * − 1), so that Replacing m 2 with γm 1 , we obtain the final expression The quantity x * Q xx + y * Q yy can be rewritten as Equivalently, we have Observe that γ > 1 implies M f > M g . The expression for M g can be further simplified to and since 0 < 1 − λ < 1 and D * m1 < 1, it is immediate that M g > 0. Hence, both quantities M f and M g are positive, and the quantity x * f ′′ M f + y * g ′′ M g is negative as long as x * , y * > 0.
, we immediately have: Under the conditions of Theorem 4, system (3)-(5) has an interior steady state in T which undergoes a supercritical Hopf bifurcation when k 1 passes throughk 1 . There exists a δ > 0 such that for all k 1 ∈ (k 1 ,k 1 + δ), (3)-(5) has an asymptotically stable periodic solution having two Floquet multipliers inside the unit circle.
Proof. The statements regarding the occurrence of the Hopf bifurcation are straightforward. Let p(t) = (p 1 (t), p 2 (t)) denote the periodic solution of (4)-(5) of period T , and J r (t) the T -periodic linearization at p(t). Notice that (3)-(5) is equivalent to the following system:V This system has a periodic solution (1, p 1 (t), p 2 (t)), and the linearized system along this periodic solution has the following matrix The fundamental matrix solution associated with J(t) takes the form and Φ r (T ) is the fundamental matrix associated with J r (t). By Theorem 4 the eigenvalues of Φ r (T ) are 1 and r 1 ∈ (0, 1). This implies that two of the Floquet multipliers of (1, p 1 (t), p 2 (t)) are inside the unit circle.
3.2. The divergence criterion and the vague attractor property. In this section, we establish a formal connection between the divergence criterion [20] and the classical method for analyzing the Hopf bifurcation [16,19,6,14]. In the classical approach, the direction of the Hopf bifurcation is determined by checking the vague attractor property at the point of bifurcation. Specifically, if the quantity V ′′′ (0) is negative then the Hopf bifurcation is supercritical, and if V ′′′ (0) is positive then the Hopf bifurcation is subcritical. The divergence criterion allows to transform the original system into one with sign definite divergence near the Hopf point without changing the geometry of the nearby orbits [21]. Knowing that the divergence remains of one sign allows to apply the Dulac criterion and conclude that the Hopf point is a weak attractor/repellor [15,20]. The resulting conclusion is typically weaker than the classical result: one can only infer the existence of periodic solutions with small amplitude on either side of the bifurcation value. In contrast, the fact that V ′′′ (0) = 0 implies the existence of the uniquely determined family of bifurcating limit cycles. The bifurcating limit cycles are stable if V ′′′ (0) < 0 and unstable if The transformation that preserves the orbits generically results in a system whose divergence is asymptotic to a sign definite or semi-definite quadratic form [21]. In the following Theorem, we establish that for such systems the sign of the quadratic form for the divergence coincides with the sign of V ′′′ (0). This result closes the existing gap between the divergence criterion and the classical analysis of the Hopf bifurcation.
Theorem 5. Consider a sufficiently smooth systeṁ that has an equilibrium at (0, 0), such that the variational matrix has pure imaginary eigenvalues. Denote the map of the first return to the positive x-axis as V . Suppose that the divergence of (14) has the form Then the following implications hold (i) If the quadratic form ax 2 + 2bxy + cy 2 = 0 is negative semi-definite, then V ′′′ (0) > 0; (ii) If the quadratic form ax 2 + 2bxy + cy 2 = 0 is positive semi-definite, then V ′′′ (0) < 0.
Proof. We will only prove the implication (i). The proof of (ii) is similar. By interchanging the roles of x and y if necessary, we may assume that r 0 x > 0. Denote the vector field of the system (14) as F = (q, r). Consider the orbit of (14) starting at (r, 0). Let (V (r), 0) be the point of the first return of this orbit to the positive x-axis (see Fig. 2). Let A(r) denote the region enclosed by this orbit as described in Fig. 2. The Green's Theorem implies that where the boundary ∂A(r) is treated as a positively oriented contour with the unit outward normal n and arclength measure ds. Since F · n ≡ 0 along any orbit of the full system, we can express the integral on the l.h.s of (15) as The function V (r) is the map of the first return to the positive x-axis. The elliptical region E(r) represents the interior of the orbit of the linearized system that passes through (r, 0). The region A(r) represents the interior of the region bounded by the part of the orbit of the full system starting at (r, 0) and ending at (V (r), 0), and the part of the x-axis between x = r and x = V (r).
Substituting the Taylor expansions r(x, 0) = r 0 x x+O(x 2 ), and V (r) = r+ 1 6 V ′′′ (0)r 3 + O(r 4 ) [19], and applying the Intermediate Value Theorem, we find that Now we evaluate the leading term of the integral on the r.h.s of (15). Expanding the system near (0, 0), we find that Let z(t) = (x(t), y(t)) denote the solution of the full system with the initial condition z(0) = (r, 0). Let z 1 (t) = (x 1 (t), y 1 (t)) denote the corresponding solution of the linearized system Since the eigenvalues of the variational matrix are pure imaginary, the orbit of z 1 (t) is an ellipse. We denote the interior elliptical region by E(r). For any fixed T > 0, we have the following estimate: Since the length of the elliptical orbit ∂E(r) is proportional to r, and the distiance between the boundaries ∂E(r) and ∂A(r) is at most O(r 2 ), we conclude that the difference between the areas of A(r) and (r) is of the order O(r 3 ). Taking into account the fact that div F = ax 2 + 2bxy + cy 2 + O(r 3 ), we find that The integral on the r.h.s. can be evaluated explicitly. Indeed, let E(r) = {(x, y) : αx 2 + 2βxy + γy 2 ≤ r 2 }, in fact the coefficients α, β, γ can be obtained explicitly from the variational matrix. Now we introduce auxiliary matrices where U > 0 and W ≥ 0. Let Q be an orthogonal matrix such that Q T U Q = diag(e 1 , e 2 ), and let is strictly positive. Finally, we use the polar coordinates to evaluate the integral 1 √ e 1 e 2 Therefore, we conclude that which immediately implies that This observation concludes the proof.
3.3. Three species. Having established the necessary conditions for the existence of stable periodic orbits for two species, we now examine a system with three organisms. Consider the following chemostat model: where (S, x, y, z) ∈ T := {(S, x, y, z) ∈ R 4 + | S + x + y + z ≤ 1}, D(x, y) is given by (5) and f , g, h are Michaelis Menten growth funcitons given by: We assume that there exists substrate concentrationsS andŜ such that 0 <S,Ŝ < 1; f (S) = g(S); and f (Ŝ) = h(Ŝ). For this system the set T is clearly forward invariant and hence (16) is well-posed. As in section 3.1, we introduce the variable V = S + x + y + z and use a reduction argument to study the systeṁ The main result to be derived is given by the following theorem The basic idea of the proof is to first obtain a periodic Hopf orbit in the xy plane and then to bifurcate from this orbit to a limit cycle in the first octant. While this statement of intent is concise, the actual methodology is quite tedious and lengthy. To guide the reader, we first outline the approach and then fill in the necessary details to complete the proof. The development begins with a limiting system for two organisms (3) with growth functions f (S) and g(S) and a diluter rate D(x, y) as described above. By Theorem 4, for sufficiently small positive κ there exists an asymptotically stable Hopf orbit in the x, y plane. This orbit is approximated by a perturbation technique and the system is then expanded by adding a third organism with a growth function h(S, η) where η is the second bifurcation parameter and h(S, 0) = f (S). Next we determine a neutral stability condition for the Hopf orbit with respect to the third organism and establish a neutral stability curve by obtaining a relation between η and the Poincaré-Lindstedt perturbation parameter ε. We define a Poincaré map in a plane parallel to the yz coordinate plane and for the fixed point corresponding to the Hopf orbit test the conditions for a bifurcation from a simple eigenvalue. This bifurcation result provides the existence of parametric family of periodic orbits in the first octant yet does not provide information on the stability of these orbits. We express the new orbits as an expansion in the parameter s. We also expand η as a function of s, η = η * + η 1 s + O(s 2 ), and demonstrate that the stability of the orbits are determined by the sign of η 1 . To find this sign we substitute the s-based parameter expansions into the neutral stability equation and calculate it's derivative with respect to s at s = 0. By combining like-terms of ε in the resulting equations an expression for η 1 is found. Finally, the dependence of the sign of η 1 on the relative positions of the growth functions is determined. Proof.

3.3.1.
A planar Hopf orbit and its approximation. We begin by establishing the existence of a periodic Hopf orbit in the x, y plane and approximating the orbit using a perturbation technique. For the system (3) and the conditions given in the statement of the above theorem, Theorem (4) provides the existence of the Hopf orbit for sufficiently small positive values of the bifurcation parameter κ. An approximation to the orbit is obtained by first varying the parameter κ which shifts the original equilibrium (x,ȳ) (for κ = 0) to the point (x,ŷ). Next, this new system and its equilibrium are perturbed as follows Substituting the above perturbed variables into the system and writing t = ωτ we obtain a sequence of equations by collecting and equating powers of ε. The ε 0 terms correspond to equilibrium and vanish. The ε 1 terms yield In order for this system to have a solution of period 2π we choose ω 0 so that ω 0 = ( √xȳ (k 2 −k 1 )) −1 . Utilizing this new value, a fundamental matrix solution for the above is given by We are free to choose an arbitrary initial value for (x 1 , y 1 ) and choose one such that The ε 2 terms are of the form where f i are quadratic functions of x 1 and y 1 with nontrivial coefficients. In order for the x 2 and y 2 to be 2π periodic, the following secularity condition must be After integrating and simplifying it is found that the secularity conditions can be fulfilled if and only if κ 1 = ω 1 = 0. Thus we may modify the expansions of ω and κ as follows ω = ω 0 + ω 2 ε 2 + ω 3 ε 3 + · · ·, κ = κ 2 ε 2 + κ 3 ε 3 + · · · and write the parameterization of the Hopf orbit as In addition to the ε 1 terms, the constant terms of the functions x 2 (τ ) and y 2 (τ ) will also be required in the next step of method. The general form for these functions is given by The constant terms are obtained as follows Performing this calculation we find

Neutral stability.
We now consider the system represented by (18) and note that the Hopf orbit obtained in the previous section exists in the present system when z = 0. The present focus will be to determine the conditions for the neutral stability of this orbit with respect to the third species. Following the development in Smith [24], by examining the variational equation of (18), the neutral stability condition can be expressed by Expanding the integrand with respect to ε, performing the integration, and rewriting the result in terms of the growth functions we obtain Based on the ε 0 term in G(ε) we define η =S −Ŝ. From the fact that f (Ŝ) = h(Ŝ), it is possible to define m 3 as a function of η and write m 3 (η) = m 1 (a 3 +Ŝ) It is assumed that m 3 (0) = m 1 , so that by continuity, the sign of (m 3 (η) − m 1 ) is constant in some neighborhood of η = 0. The ε 0 term is therefore rewritten asS (m 3 (η) − m 1 )η (a 3 +S)(a 1 +S) and it is noted that the sign of this expression changes as η passes through zero. We define a new function and recognize that H(ε, η) = 0 gives the neutral stability criteria. Using the previous results and definitions, the following are derived .

3.3.3.
A bifurcation result. The next step in the technique is to obtain a parametric family of periodic orbits which bifurcate form the planar Hopf orbit into first octant. We begin my modifying the third equation in (18) to reflect its dependence on η and assume that k 1 and all other parameters are fixed. Letting x(t) = p(t) correspond to the Hopf periodic orbit of the above system with z = 0, it is noted that the orbit it is asymptotically stable in the xy plane. By the appropriate translation of variables the Hopf orbit has the form p(0) = 0 and p ′ (0) = (p ′ 1 (0), 0). Next, φ(t, (x, y, z); η) is defined as a solution to the above system in the translated coordinates with φ(0, (x, y, z); η) = (x, y, z). A Poincaré map P is established in a neighborhood U of 0 in the plane H = {x = 0} where P maps U into V , with V another neighborhood of 0 in H. Corresponding to P is a smooth function τ (y, z; η) which gives the first return time for any point in U to V . It follows that τ (0, 0; η) = τ 0 , the period of the Hopf orbit p(t). Thus P (y, z; η) may be regarded as a projection of φ(τ (y, z; η), 0, y, z; η) onto the yz plane with P (0, 0; η) = (0, 0). We define the displacement mapF : U × R → H byF (y, z; η) = (y, z) − P (y, z; η) and note that periodic solutions to (18) will be zeros ofF . Expressing the component functions of P as P 1 and P 2 and using the invariance of the xy plane for the above system we may write P 2 (y, 0; η) = 0 so that the Jacobian of the Poincaré map is given by Note that since p(t) is asymptotically stable and independent of η, the inequality 0 < ∂P1 ∂y (0, 0; η) < 1 holds for all η. Expressing the component functions of φ as φ 1 , φ 2 , and φ 3 , gives P 2 (0, z; η) = φ 3 (τ (0, z; η), 0, 0, z; η), and taking derivatives yields and since φ 3 (t, 0, 0, 0; η) ≡ 0, the above simplifies to Expressing the third equation in (18) asż = G(x, y, z; η) and integrating the variational equation of the system along p(t), we find If η is an element of the neutral stability curve, then the integral evaluates to zero and ∂P2 ∂z (0, 0; η) = 1. Next we seek a particular branch of solutions ofF (y, z; η) = 0 (i.e., a collection of periodic solutions) and begin by noting .
A similar approach also yields .
Performing a Taylor series expansion of the above about s = 0 provides the equation A similar expansion for the remaining two variables yields From the above series expansions we find the following expressions in which all variables except t have been suppressed for ease of notation Substituting (25) and (26) into (24) yields after simplification dt.

Additional expansions.
In order to use (27) to determine the sign of η 1 , the following perturbation expansions are required (note x 0 (t) and y 0 (t) are chosen so that y ′ 0 (0) = 0): where the constant N is given by Note that the expansion for ∂φ3 ∂z (t) results from taking the partial derivative with respect to z of a term in the variational equation for system (18): where φ 3 (0) = 0 and ∂φ3 ∂z (0) = 1. To evaluate the integrals in (27), we must perform some calculations on certain quantities in the expansions given by (28). First the functions ψ 0 (t), δ 0 (t), γ 0 (t), and ρ 0 (t) must be determined. Considering (φ 1 (t), φ 2 (t), φ 3 (t)) as a solution to the system (18) and taking the derivative of the first two equations with respect to y gives d dt Letting φ 1 (t) = x 0 (t), φ 2 (t) = y 0 (t), φ 3 (t) = 0, substituting from (28), and collecting the ε 0 terms gives the system with the solution given by .
The initial conditions result from the fact that φ 1 (0, x 0 (0), y 0 (0), 0; η) = x 0 (0) and φ 2 (0, x 0 (0), y 0 (0), 0; η) = y 0 (0), which in turn implies ∂φ1 ∂y (0) = 0 and ∂φ2 ∂y (0) = 1 so that we conclude ψ 0 (0) = 0 and δ 0 (0) = 1. Similarly, taking the derivative with respect to z of the first two equations of system (18), substituting the ε expansions, using ∂φ3 ∂z = 1 + O(ε), and collecting the ε 0 terms gives a system with the following solution Additionally certain quantities involving ψ 1 (t), δ 1 (t), γ 1 (t) and ρ 1 (t) are needed. Proceeding as above we collect the ε 1 terms of the expansions and find a system of the form Letting Φ(t) be the fundamental matrix for the homogeneous problem, then the following equations are developed which when solved yields ψ 1 (τ 0 ) = δ 1 (τ 0 ) = 0 and A similar system is obtained for γ 1 (t) and ρ 1 (t) and using the previous techniques gives ρ 1 (τ 0 ) = −ω 0 πh ′ (S). As a final preparatory step, we note that in their text, Marsden and McCracken define a displacement function by V (y) = P 1 (y, 0) − y. They demonstrate that V (0) = V ′ (0) = V ′′ (0) so that by using a Taylor expansion it is possible to write V (y) = 1 3! V ′′′ (0)y 3 + · · · (Note that the dependence of V on the bifurcation parameter has been suppressed). The above implies where y has been replaced by ε since they agree to the first order. (Recall the fixed point of the Poincare map corresponds to (0, 0).) This gives The sign of η 1 . Having performed all of the above calculations, we are in a position to solve (27) for η 1 . We begin by expanding this equation with respect to ε, performing the integrations, simplifying, collecting the like powers of ε, and equating their coefficients to zero. The result of these operations will be an expression for η 1 in terms of the system variables, the growth functions, and their derivatives. To aid in the simplification of (27), we perform some preliminary computations. Based upon (28) and the discussion at the end on section 3.3.5 the following equality holds Using the definition of m 3 (η) and η * , one can easily derive the expansions For the first term of (27), using (28) and (31) the first factor may be expanded as From the development of h(S), if m 3 (η) =m 3 , then h(S) = f (S) so that the ε 0 term vanishes in the above expression. For the second factor, we deduce from section 3.3.5 that ψ 0 (τ 0 ) = ρ 0 (τ 0 ) = γ 0 (τ 0 ) = ψ 1 (τ 0 ) = 0 and upon expansion obtain Combining all of this information, the first term is found to be of order O(ε 4 ). For the second term of (27), the numerator of the integrand can be expanded as a 3 (ε 0m 3 + O(ε 2 ))(ε 0 γ 0 (t) + ε 1 γ 1 (t) + ε 0 ρ 0 (t) + ε 1 ρ 1 (t) + 1), and the denominator can be written as Combining these expressions with (30) yields N ω o x y 2πω0 0 a 3m3 (γ 0 (t) + ρ 0 (t) + 1) From (29) it is clear that γ 0 (t) + ρ 0 (t) + 1 will integrate to zero so the above expression is of order O(ε 4 ). The numerator of the integrand in the third term may be written as a 3 (ε 0m 3 +O(ε 2 ))(ε 0 ψ 0 (t)+ε 1 ψ 1 (t)+ε 0 δ 0 (t)+ε 1 δ 1 (t)+O(ε 2 ))(ε 1 ρ 1 (τ 0 )+ε 2 ρ 2 (τ 0 )+O(ε 3 )) while the denominator may be expanded as Combining this information, the integrand may be expressed as It is clear that the first two terms will integrate to zero while the remaining ε 2 term may be written as After integration and using (28) the third term is given by x y In the fourth term, we note that for the integral, the coefficient of ε 0 is given by 2πω0 0k 1 γ 0 (t) +k 2 ρ 0 (t)dt.
From (29), it is evident that this integral is equal to zero, so that using (30), the fourth term is of order O(ε 4 ). For the fifth term of (27), the integrand may be expanded as which reduces to Each expression which involves ψ 0 (t) or δ 0 (t) integrates to zero so using (28), this term can be expressed as Analysing the final term, we expand the integral as Combining this with (30) and (31) the last term may be written as Using these expansions for the individual terms, we may rewrite (27) as x y Rearranging this equation, evaluating integrals, and substituting we obtain 3.3.7. The stability of bifurcating orbits. As discussed previously, the stability of the bifurcating orbits depends on η 1 , ∂ 2 η ∂ε 2 (0), andS (m3−m1)η (a1+S)(a3+S) . By defining We observe from the statement of the problem as well as the definition of the growth functions that P has the same sign as (a 1 − a 3 ) and since ρ 1 (τ 0 ) is strictly negative, the two terms in (32) are of opposite sign. Having established the above identities we wish to examine Q and to determine the sign of Q as a function of a 3 . Using the fact that f (S) = g(S) = h(S) when η = 0, it follows that a 3 = a i gives m 3 (0) = m i for i ∈ {1, 2}. Thus a 3 = a 1 impliesh(S) = f (S) and a 3 = a 2 impliesh(S) = g(S) (recallh(S) is the growth function with m 3 = m 3 (0)), so that in either case, Q = 0. Using the definitions of the uptake functions it is possible to express Q as a third degree polynomial in a 3 as follows Since, as previously discussed, a 2 > a 1 , the term 2m1m2(a2−a1) (a1+S) 3 (a2+S) 3 (a3+S) 3 is positive, which means that the constant term in the polynomial as well as the coefficient of a 3 3 are both positive. This implies Q has one negative root and the two positive roots given by a 1 and a 2 . From this information it is possible to deduce that Q > 0 for a 3 ∈ (0, a 1 ) ∪ (a 2 , +∞) and Q < 0 for a 3 ∈ (a 1 , a 2 ).
We are now in a position to determine how the interactions of the three growth functions affect the bifurcation. First, Table 1 is constructed based upon the results derived above. From the information in this table as well as the locally quadratic Table 1. The Signs of Essential Quantities as a Function of a 3 .
it is possible to construct the following graphs of the neutral stability curves where U represents the region of instability for the planar Hopf orbit, S represents the region of stability for the planar Hopf orbit, and the arrow points in the direction of the change of η as the value of s becomes positive in the expansion η = η * + sη 1 + O(s 2 ). These graphs demonstrate that for a 3 < a 1 or a 2 < a 3 , the Hopf orbit loses its stability and the bifurcating orbit in the first octant becomes asymptotically stable. The opposite is true for the case a 1 < a 3 < a 2 .
The Figure 4 indicates the relative positions of the growth functions as the value of a 3 is varied. (Recall that g ′ (S) > f ′ (S) implies that f (S) > g(S) for S ∈ (0,S).) Correlating the information in Figures 3 and 4, we deduce that a stable periodic orbit in the first octant results from a bifurcation when h ′ (0) is either greater than or less than both f ′ (0) and g ′ (0) while the bifurcating orbit is unstable otherwise. Thus our initial claim is verified. The proof of Theorem 6 is now complete.

4.
A numerical example. To illustrate the results obtained in Theorem 6, we performed a numerical study of the model (18) with a specific choice of growth  The derivatives of f and g are given byk 1 = f ′ (0.5) = 0.5714, andk 2 = g ′ (0.5) = 1.2 respectively. We set D(x, y) = 1.45 − k 1 x − k 2 y with k 1 =k 1 + 0.05 and k 2 =k 2 so that Lemma 1 warrants the existence of a limit cycle in the (x, y) plane which is asymptotically stable in this plane. Varying the value of the parameter m 3 between 1.1998 and 1.1978, we located the stable periodic solutions by forward numerical integration in Mathematica 5.2 [27]. The orbits of these periodic solutions are shown in Fig. 5. Numerical estimates of the periods T and the Floquet multipliers µ i corresponding to these periodic solutions are presented in Table 2. Near the value m * 3 ≈ 1.1976, the limit cycle in the (x, y) plane undergoes a transcritical bifurcation, and stable positive periodic solutions exist for m 3 > m * 3 . As m 3 increases towards the critical value 1.2, the periods and amplitudes of the periodic solutions decrease monotonically. Another bifurcation occurs at m 3 = 1.2 where the periodic solutions interact with the line of equilibria. The analysis of this bifurcation is a subject for a future study.
5. Discussion. In this paper, we have analyzed the dynamics of microbial competition in the chemostat which is controlled by means of state dependent feedback. We have considered the class of feedbacks where the dilution rate is a positive affine function of the microbial concentrations. In case of two species competing in the chemostat, we showed that the parameters of the feedback function may be varied in such a way that the system undergoes a Hopf bifurcation. We have proved that only supercritical Hopf bifurcations occur when the growth rates of both species  are governed by Michaelis-Menten kinetics (Theorem 4). Thus we established the existence of stable periodic solutions in the two species case. In addition, we demonstrated that under certain conditions on the growth functions, adding a third species results in a transcritical bifurcation of limit cycles, where a periodic solution bifurcates into the region of coexistence for all three species (Theorem 6). We analyzed the direction of the transcritical bifurcation and the stability of the bifurcating periodic solutions. We presented a numerical example illustrating the conclusion of Theorem 6. The importance of the present study is to point the range of dynamic possibilities for the coexistence of species by means of state-dependent control where coexistence is impossible in the absence of such control.
The theoretical question of the existence and stability of stable periodic solutions for growth functions with less stringent constraints remains open. Another open and interesting question is the global dynamics in the two species case. In the present report, we have left out several interesting details such as the existence of multiple saddles on the boundary, possible existence of heteroclinic connections, and resulting multistability. These questions will be addressed in the future reports.