GLOBAL DYNAMICS IN A CHEMOSTAT WITH MULTIPLE RESOURCES AND TOXIN

A mathematical model for the dynamics of the bacterial population in response to toxin challenge is presented. Then it is developed the theory by analyzing two extreme cases (no control and constant control) and the intermediate case. This leads to a local sufficient condition for treatment success of waste water. Also we give a condition for global treatment, supported by numerical simulations of the model.


Introduction
Water pollution by industrial chemicalwaste represents a relevant topic with a growing awareness nowadays, due to its continuous release in environment, which is a potential cause of social and environmental problems [1].In biological wastewater treatment, bacteria grow on and consume various constituents (e.g.organic carbon, nitrogen, phosphorus) in the water and are subsequently removed as sludge.The bacteria can grow suspended in the water and/or attached to solid media.In the attached growth processes, bacteria can interact spatially by setting up zones with different environmental conditions.For example, nitrifying bacteria may grow on the aerobic surface and denitrifying bacteria in the anoxic interior of a biofilm.In continuous flow systems (e.g.conventional activated sludge), the bacteria grow in "mixed liquor" reactor and are settled out in subsequent clarifier, and a portion of the settled biomass is recycled back to the biological reactor.In batch systems (e.g.sequencing batch reactor), the various processes are carried out in a single reactor by altering the conditions (i.e.mixing, aeration) [2].
It is important and useful to model some control strategies against the water pollution by using numerical simulations.Gujer [3] developed a relatively simple individual-based models (IBM) for continuous and batch suspended growth systems, based on the Activated Sludge Model (ASM), that accounts for one bacteria functional group and microbial storage products (MSP).Schuler [4][5][6] developed a more complex and realistic model for the continuous treatment scheme also based on ASM, which accounts for multiple functional groups and MSPs, as well as variable biomass.
Recently, many papers have been devoted to the analysis of mathematical models describing the growth of the organism [7][8][9][10][11].In [9] transient oscillations in cell number density are experimentally observed during cell growth, but considering a simple autonomous model with single nutrient.Li and Smith [10] study a chemostat model that describes competition between two species for two essential resources based on storage, but ignoring the effect of the inhibitor (toxicant).Uncontrolled contribution of pollutant to the environment has led many species to extinction.In order to use and regulate toxic substance wisely, we must assess the risk of the population exposed to toxicant [7].Therefore, it is important to study the effects of toxicant on populations and to find proper control strategies against the water pollution.
As motivated by the above-mentioned brief literature survey, in this paper, we will construct a chemostat with multiple nutrients and toxin.For mathematical simplicity, we make several assumptions in the models.Details are shown in next section.
The paper is organized as follows.We begin by describing the model for the dynamics of the bacterial population in response to toxin challenge.Then we give a local sufficient condition for treatment success of waste water by analyzing two extreme cases (no control and constant control) and the intermediate case.Finally, the condition for global treatment is obtained, supported by numerical simulations of the model.

Model formulation
The equations describe a population of bacteria in waste water.The variable X 1 is the bacterial biomass which is sensitive to toxin and X 2 is one which is insensitive to toxin.
where the function F j (t) ≥ 0 is the flux of nutrient, the constant c j > 0 is the fraction of bacterial biomass consisting of nutrient S j (∑ c j = 1), S = (S 1 , S 2 , • • • , S m ) T .The parameters d 1 > 0 and d 2 > 0 are removal rates of sensitive bacteria and insensitive bacteria, respectively.We expect that the insensitive bacteria have an advantage in lower removal rates: In the special case that resources are supplied at a constant rate F j (t) = dS 0 j , where S 0 j > 0 is the concentration of S j in the feed to a well-mixed continuous culture with dilution rate d and assuming that sensitive bacteria and insensitive bacteria are removed at rate d (ignoring death).
Our system becomes mathematically tractable because various conservation principles become available.Then the simplified system becomes (2.2) The following assumptions for model (2.2) are made: (H1) The functions k 1 (t) and k 2 (t) are non-negative, time-varying functions which describe the effect of toxin on the population.
is the killing rate of the sensitive population.Note in particular that the killing rate is proportional to the growth rate of the bacteria.It is positive when both toxin and nutrient are present, but zero when either one is missing.
(ii) k 2 (t) is the (per capita) rate at which insensitive bacteria revert to the sensitive state when toxin is absent (it is zero when toxin is present).
(iii) For simplicity, assume that the functions k 1 (t) and k 2 (t) are τ−periodic (τ > 0) and of the following the bang-bang type [12]: The rate k f (s)B s is the loss rate of sensitive bacteria that become insensitive bacteria, regardless of whether or not the toxin is administered.
(v) Assume that the constant loss of sensitive bacteria to the insensitive compartment does not prevent a net positive growth of the sensitive class, that is 1 − k > 0. But the additional effect of the toxin is lethal to the sensitive bacteria, that is 1 − k 1 − k < 0.
(H2) Nutrient uptake is governed by the nutrient uptake rate functions f (s), which is assumed to satisfy the following hypotheses: be the portion of the ray through S 0 in the direction −c belonging to R m + .(i) Assume that is continuous and C 1 on a neighborhood of L and satisfies ∇ f (S) • c > 0 on L.
(ii) Assume that f (S) − d < 0 for S ≤ S 0 − c(min 1≤ j≤m S 0 j /c j ).Less general conditions, involving only f , which imply those above for every choice of S 0 and c > 0 are: The main purpose of this paper is to investigate how the behavior of system (2.2) changes both qualitatively and quantitatively in terms of p.

Main results
Lemma 3.1.System (2.2) has R m+2 + as a forward invariant set, and its solutions are ultimately bounded by a uniform bound.
Proof.The first assertion is obvious.The second follows from consideration of the dynamics of given by and hence Note that solutions of (2.2) are attracted to the set where S ≤ S 0 − c(X 1 + X 2 ).The domain of (3.1) is the triangular region Proof.(ii) Case p = 0 (no control ).Assume The eradication steady state E 0 is unstable, and there exists a unique positive steady state which is locally stable.
Proof.(i) If p = 1, then k 1 (t) ≡ k d and k 2 (t) ≡ 0 for all t.Then system (3.1)becomes (3. 3) The Jacobian matrix of (3.3) at a point (X 1 , X 2 ) is given by Obviously, Jacobian J at E 0 is where f is evaluated at E 0 .By (2.4) and together with −d < 0, the eradication steady state E 0 is locally stable.
(ii) If p = 0, then k 1 ≡ 0 and k 2 (t) ≡ k 2 for all t.Then system (3.1)becomes (3.4) The Jacobian matrix of (3.4) at a point (X 1 , X 2 ) is given by where f , f X 2 and f X 1 are evaluated at (X 1 , X 2 ).Obviously, the Jacobian J at E 0 is whose determinant is negative because of (3.2).Hence the Jacobian matrix at E 0 has one positive and one negative eigenvalues, and thus it is unstable.
Then let us establish the existence of the unique positive steady state E * .Define A calculation shows that the determinant of the coefficient matrix of system (3.1) is −dG(S), where .
And since h(min j S 0 j /c j ) < 0 by hypothesis there is a unique To show that it is locally stable we determine the determinant and the trace of the Jacobian matrix evaluated at E * .The calculation is straightforward and yields the following formulas: and We therefore conclude that the 'coexistence' steady state E * , when it exists, is asymptotically stable.
Proof.Using Dulac's function 1/(X 1 X 2 ), if p = 1, we have If p = 0, the following inequality is valid.Proof.(i) If p = 1, Lemma 3.3 implies that E 0 is locally stable.It is globally attractive by the Poincaré-Bendixson theory since there is no other equilibrium in the positively invariant set Γ.
(ii) If p = 0 and (3.2), G(S 0 ) > 0 and by Lemma 3.3 E 0 is a saddle whose stable manifold lies outside Γ.Since Γ is a bounded and positively invariant, the Poincaré-Bendixson theory implies that there exists another equilibrium.In addition, by Lemma 3.4, the nontrivial equilibrium E * is the global attractor.
Next we deal with the τ-periodic model (3.1), assuming that p ∈ (0, 1).We will also assume that the dilution rate is not too large, as in condition (3.2) of Lemma 3.3.
Notice that system (3.1) has a steady state E 0 , regardless of the value of p.To determine its stability properties, we define x = (X 1 , X 2 ) and calculate the τ−periodic variational equation at Using (2.2), these Floquet multipliers are the eigenvalues of the following matrix: where and are quasimonotone matrices.Notice also that all eigenvalues of A 1 are in the open half-plane by (v) of hypothesis (H1), and that A 2 has one negative and one positive eigenvalue, as shown in the proof of Lemma 3.3.
Since A 1 and A 2 are quasi-monotone, it follows that their matrix exponentials are (entry-wise) non-negative matrices and then their product Φ is a (entry-wise) positive matrix whose spectral radius ρ(φ ) is an eigenvalue by the Perron-Frobenius theorem [13].Consequently, to determine stability of E 0 , we need to establish whether or not ρ(Φ) is inside the unit circle: If ρ(Φ) < 1, then E 0 is locally asymptotically stable.If ρ(Φ) > 1, then E 0 is unstable.
Our main concern is knowing how ρ(φ ) varies as a continuous function of p (this variation is continuous since eigenvalues of a matrix are continuous functions of its entries, and clearly the entries of φ are continuous in p).For p = 0 (never using toxin), and hence also for p near 0 by continuity of ρ(φ ), we have that ρ(φ ) = ρ(e τA 2 ) > 1 because A 2 has a positive eigenvalue.This is in accordance with Theorem 1 where it was shown that the microbial population exists and E * is globally asymptotically stable.For p = 1 (using toxin continuously) we have by (v) of hypothesis (H1) that ρ(φ ) = ρ(e τA 1 ) = e −dτ < 1.The inequality also holds for p near 1.
This is in accordance with Theorem 3.1 as well because it was shown there that all solutions converge to E 0 in this case.

Simulations and discussion
In the following, we describe results of numerical simulations and present some conclusions.
These numerical experiments illustrate the extreme treatment.The numerical values are shown in Table 1, and it can be easily verified that these choices satisfy the condition (v).The per capita growth rate is of Michaelis-Menten type: Let us assume for definiteness, although one can argue more generally, that there are only two resources S 1 and S 2 .Consider the function f (S 1 , S 2 ) = f 1 (S 1 ) f 2 (S 2 ) [14,15] as the resourcedependent growth function.
Table 1 m 2 The first two scenarioes is computed to check dynamic consistency between the theoretical results obtained in Theorem 3.1 and the numerical simulation of the model.If p = 1, the eradication steady state E 0 of system (3.1) is globally asymptotically stable (Fig. 1).If p = 0 and (3.2), the positive steady state E * of system (3.1) is globally asymptotically stable (Fig. 2).
The third and forth scenario is simulated varying the control duration µ while keeping the other parameters fixed.This numerical simulation allows to observe the effect of this paramter on the treatment of waste water.When p = 0, the level of the bacterial boomass increases if µ = 0.6 (Fig. 3), while the level of the bacterial boomass decreases if p = 0.2 (Fig. 4).
Futhermore, in Theorem 3.2, we show that when p ∈ (0, 1), the eradication steady state E 0 of system (3.1) with (2.2) is locally stable if ρ(φ ) < 1, but unstable if ρ(φ ) > 1.Though we do not describe results of numerical simulations of the intermediate cases, it is easy to see the behavior of system (2.2) changes both qualitatively and quantitatively in terms of p.Therefore, it is important to study the effects of toxicant on populations and to find proper control strategies against the water pollution.
Clearly, our model is a simplification of reality because the concentration of toxin is not expected to be of the bang-bang type.In a more realistic model, the functions k 1 (t) and k 2 (t) would be replaced by functions depending on (at least) a new state variable for the concentration of the antibiotic, and the periodicity would arise through a periodic forcing term in the equation for this new variable.We leave the study of such a model to the future.

FIGURE 1 .FIGURE 2 .
FIGURE 1. Phase plane trajectories corresponding to different initial levels.The figure clearly indicate that the eradication equilibrium E 0 is asymptotically stable when p = 1 and µ = 0.4.