# Periodic Oscillations in a Chemostat Model with Two Discrete Delays

Authors: Tiansi Zhang; Xuehui Ji; Bo Li
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

Periodic oscillations of solutions of a chemostat-type model in which a species feeds on a limiting nutrient are considered. The model incorporates two discrete delays representing the lag in nutrient recycling and nutrient conversion. Through the study of characteristic equation associated with the linearized system, a unique positive equilibrium is found and proved to be locally asymptotically stable under some conditions. Meanwhile, a Hopf bifurcation causing periodic solutions is also obtained. Numerical simulations illustrate the theoretical results.


---

## Body

## 1. Introduction

The chemostat is very important to imitate ecosystem and is used extensively in microbiology and population biology for its measurable parameters and maneuverable experiments; see [1, 2], for example. Many models have been constructed based on chemostat-type equations such as the growth of plankton of unicellular algae and plasmid-bearing and plasmid-free competition; see [3–5] and the references therein. Amongst them, regeneration of nutrient decomposed by bacterial played an indispensable role for the stability and development of the nature and thus was taken into account; see [6–8], but no one treated it with a time delay. In fact, a delay in nutrient recycling on ecosystem always exists and is influenced by the change of temperature [9]. In view of this, Bischi [10] embedded this time delay and discussed its stabilizing effect for a chemostat model with oscillations. After that, He, Ruan and Wang et al. studied the global behaviors of a chemostat model with delayed nutrient recycling and periodically pulsed input in [11] and [12], respectively. Freedman and Xu [13] turned to consider a competition model with delayed nutrient recycling which generalized the problem in [14] and gave some persistence and extinction conditions for the competing populations. For some other excellent work about the delayed nutrient recycling, one may refer to [15–20].On the other hand, environment is a vital factor that affects the growth of microorganisms. When it changes, a lag phase of the growth would appear. In 1969, when studying the growth ofIsochrysis galbana, Caperon [21] used a time lag in his model with alterable environment and his research well met the corresponding experiment results. Besides it, Ellermeyer [22] and Xia et al. [23] also investigated the delayed response in growth and got global asymptotic behaviors of a competition and transient oscillations, respectively.From the above we see that time delays have biological significance since they can alter dynamics of solutions of a system and display more transient oscillations especially in numerical simulation. So in this paper, we confine us to a single species model with a limited nutrient at a constant rate together with two delays in nutrient recycling and nutrient conversion, which takes the following form:(1)dSdt=D(S0-S)-aP(S)x+γD1x(t-τ2),dxdt=ae-Dτ1PSt-τ1xt-τ1-Dx-D1x,
with initial value conditions
(2)S(θ)=φ(θ)≥0,x(θ)=ψ(θ)≥0,
where S(t) denotes the limiting nutrient concentration and x(t) denotes the concentration of the microorganism at time t. φ,ψ∈BC[-max⁡{τ1,τ2},0], the Banach space of all continuous bounded functions, and all parameters are positive constants. S0 is the input concentration of the limiting nutrient, D is the washout rate of the chemostat, D1 is the death rate of the microorganisms, γ∈(0,1) is the fraction of nutrient conversion, and a is the maximal nutrient uptake rate for the microorganisms. There are two delays of τ1 and τ2 due to nutrient conversion and nutrient recycling, respectively. P(S) is the nutrient uptake function and has the following properties:
(3)P(0)=0,dSdt>0,lim⁡S→∞P(S)=1.
This is the first time to study such a problem. The method used here is established by Ruan and Wei [24]. Our analysis will develop from discussing roots of the associated characteristic equation of the linearized system, whose coefficients depend heavily on delays.The paper is organized as follows. In the next section, we investigate the existence of a positive equilibrium. Then in Section3, we study the asymptotic stability of the equilibrium by first linearizing the system and then analyzing the roots of the associated characteristic equation, which is a second degree exponential polynomial. This process is divided into three cases of which two delays are zero, only one delay is zero, and no delays are zero. Finally a Hopf bifurcation causing the appearance of periodic solutions is found to exist. The last two sections are to numerically illustrate our findings and to summarize the discussion in the paper, respectively.


## 2. The Boundedness and Existence of Nonnegative Equilibria of Model1

Firstly, by studying model1 with 3, we have the following result about the boundedness of solutions.Theorem 1.
All solutions of system1 with 3 are bounded.Proof.
Define a Lyapunov function(4)Vt,S,x=S(t)+eDτ1x+∫t-τ1t‍aP(S(u))x(u)du+γD1∫t-τ2t‍x(u)du.
Obviously, V≥0 and V→+∞ as S,x→+∞. Differentiate V along trajectories of 1; that is,
(5)V˙=D(S0-S)-aP(S)x+γD1x(t-τ2)+eDτ1(ae-Dτ1P(S(t-τ1))x(t-τ1)-Dx-D1x)+aPSx-PSt-τ1xt-τ1+γD1xt-xt-τ2=D(S0-S)-((D+D1)eDτ1-γD1)x.
It is easy to see that the coefficient of x is negative since γ∈(0,1). So outside the region bounded by the positive orthant and the line
(6)x=D(S0-S)(D+D1)eDτ1-γD1,V˙ is negative. This completes the proof.Secondly, note thatE0=(S0,0) is always an equilibrium for system 1, and from the following equations
(7)D(S0-S*)=(aP(S*)-γD1)x*,aP(S*)x*=(D+D1)eDτ1x*,
system 1 has a positive equilibrium E*(S*,x*) given by
(8)S*=P-1D+D1aeDτ1,x*=D(S0-S*)aP(S*)-γD1,
provided that
(9)D+D1aeDτ1<P(S0),
where P-1 means the inverse function of P. Moreover, since the function P(S) is increasing, the equilibrium E* is unique. So we can conclude the following.Theorem 2.
System1 has a unique trivial equilibrium E0, and has a unique positive equilibrium E* defined by 8 if condition 9 holds.


## 3. Asymptotic Stability ofE*

In this section, we are going to discuss the asymptotic stability of the positive equilibriumE* of system 1. Throughout this section, we assume that condition 9 holds. To begin we linearize 1 at E* and then construct the associated characteristic equation.Settingz1=S-S* and z2=x-x* and putting them into system 1, we have
(10)dz1dt=-D+ax*P′S*z1-aP(S*)z2+γD1z2(t-τ2),dz2dt=-(D+D1)z2+ax*P′S*e-Dτ1z1(t-τ1)+aP(S*)e-Dτ1z2(t-τ1).
Denote
(11)A0=-D+ax*P′S*-aP(S*)0-(D+D1),A1=00ax*P′S*e-Dτ1aP(S*)e-Dτ1,A2=0γD100.
Then the characteristic equation associated with 10 is given by
(12)det⁡(λI-A0-e-λτ1A1-e-λτ2A2)=0,
or equivalently by the second degree exponential polynomial in λ,
(13)λ2+D+ax*P′S*+D+D1λ+D+D1D+ax*P′S*-aP(S*)e-Dτ1(λ+D)e-λτ1-aγD1e-Dτ1x*P′S*e-λ(τ1+τ2)=0.From now we restrict our attention to discuss the local asymptotic stability ofE* by analyzing the sign of the real parts of roots of the characteristic equation 13. Note that the coefficients of 13 depend heavily on the time delays τ1 and τ2; there is not a direct way to go, so by use of the method proposed by Ruan and Wei in [24], we consider three cases: τ1=0,τ2=0; τ1=0,τ2>0; and τ1>0,τ2>0.


### 3.1. The Caseτ1=τ2=0

Assume in13 that τ1=τ2=0; then we have
(14)λ2+D+ax*P′S*λ+ax*P′S*aPS*-γD1=0.
According to the Routh-Hurwith criterion, the real parts of all eigenvalues of 14 are negative if and only if 9 is valid, which implies the asymptotic stability of E*. Hence, we have the following result.Theorem 3.
Assume that condition9 and τ1=τ2=0 hold. Then all roots of 13 have negative real parts, and the positive equilibrium E* of system 1 is locally asymptotically stable.


### 3.2. The Caseτ1=0, τ2>0

Settingτ1=0 in 13, since aP(S*)=D+D1, the characteristic equation is now
(15)λ2+Aλ+B-Ce-λτ2=0,
where
(16)A=D+ax*P′S*,B=a2x*PS*P′S*,C=γD1ax*P′S*.
It is not hard to see that A>0 and B>C>0.We know from Theorem3 that the E* is locally asymptotically stable if set τ2=0 in 15, while when τ2 increases from zero, the equilibrium could lose its stability if pure imaginary roots of 13 appear. Therefore it is convenient to look for pure imaginary roots of 15.Putλ=iω, ω≠0∈R into 15 and separate the real and the imaginary parts of the equation; then we have
(17)cos⁡(ωτ2)=B-ω2C,sin(ωτ2)=-AωC.
Taking the sum of squares in 17 yields
(18)ω4+(A2-2B)ω2+B2-C2=0.
It has real solutions
(19)ω±2=122B-A2±2B-A22-4(B2-C2),
provided
(20)A2-2B2>4(B2-C2)>0>A2-2B.
Notice that if ω is a solution of 17, so is -ω. Accordingly 15 has two pure imaginary solutions λ±=iω± with ω+>ω->0.Now differentiate both sides of15 with respect to τ2,
(21)(2λ+A+Cτ2e-λτ2)dλdτ2=-Cλe-λτ2,
and then
(22)dλdτ2-1=-2λ+Aeλτ2Cλ-τ2λ,eλτ2=Cλ2+Aλ+B.
Therefore,
(23)sign⁡dReλdτ2λ=iω=sign⁡Redλdτ2-1λ=iω=sign⁡Re-2λ+Aλλ2+Aλ+Bλ=iω=sign⁡A2-2B-ω2A2ω2+B-ω22=sign⁡A2-2B+2ω2.
Clearly the sign of d(Reλ)/dτ2 is positive at ω+ or is negative at ω-.Furthermore, we see that there are two sets ofτ2 satisfying 17:
(24)τ2,1n=θ1ω++2nπω+,
where 0≤θ1<2π, n=1,2,3,…, and
(25)cos⁡θ1=B-ω+2C,sinθ1=-Aω+C,τ2,2n=θ2ω-+2nπω-,
where 0≤θ2<2π, and
(26)cos⁡θ2=B-ω-2C,sinθ2=-Aω-C.
To sum up, we have the following.Lemma 4.
Letτ2,i* be an element of the sequences {τ2,in}, i=1,2. For τ2=τ2,i*, the characteristic equation 15 has a pair of simple conjugate pure imaginary roots ±iω±, satisfying
(27)sign⁡d(Reλ)dτ2τ=τ2,i*=sign⁡{A2-2B+2ω±2}.From the above analysis, the stability of the positive equilibriumE* of system 1 for τ1=0 is obvious.Theorem 5.
Assume that conditions9 and 20 hold and τ1=0. Let
(28)τ2*=min⁡i=1,2;n∈N{τ2,in}.
Then the equilibrium E* of system 1 is locally asymptotically stable when τ<τ2* and a Hopf bifurcation occurs at E* when τ=τ2* if and only if sign⁡{A2-2B+2ω±2}>0.


### 3.3. The Caseτ1,τ2>0

The caseτ1,τ2>0 can be regarded as a general one of that in the last section, so our discussion here is based on the result in Section 3.2. The asymptotic stability of E* can be derived from Theorem 5. To state it, we first introduce a lemma about the sign of the real parts of roots of the characteristic equation 13.Lemma 6.
If all roots of15 have negative real parts for τ2>0, then there exists a τ1*(τ2)>0, such that all roots of 13 have negative real parts for τ1∈[0,τ1*(τ2)).Proof.
Suppose that15 has no roots with nonnegative real parts for τ2>0; that is, all the roots of 13 have negative real parts for τ1=0 and τ2>0. According to Theorem  2.1 in [21], on the left side of 13, the sum of the multiplicity of zeros may change if a zero locates on or passes through the imaginary axis in the open right half-plane when τ1 varies.
Recall that13 with τ1=0 has no root with nonnegative real part; then there must exist a τ1*(τ2)>0 such that all roots of 13 have negative real parts for τ1<τ1*(τ2).Theorem 7.
Assume that conditions9 and 20 hold. Then for any τ2∈[0,τ2*), where τ2* is defined in Theorem 5, there exists a τ1*(τ2)>0 such that the positive equilibrium E* of system 1 is locally asymptotically stable when τ1∈[0,τ1*(τ2)).Proof.
From Theorem5, it follows that all roots of 15 have negative real parts when τ2∈[0,τ2*). Then Lemma 6 indicates the validity of the conclusion.


## 3.1. The Caseτ1=τ2=0

Assume in13 that τ1=τ2=0; then we have
(14)λ2+D+ax*P′S*λ+ax*P′S*aPS*-γD1=0.
According to the Routh-Hurwith criterion, the real parts of all eigenvalues of 14 are negative if and only if 9 is valid, which implies the asymptotic stability of E*. Hence, we have the following result.Theorem 3.
Assume that condition9 and τ1=τ2=0 hold. Then all roots of 13 have negative real parts, and the positive equilibrium E* of system 1 is locally asymptotically stable.


## 3.2. The Caseτ1=0, τ2>0

Settingτ1=0 in 13, since aP(S*)=D+D1, the characteristic equation is now
(15)λ2+Aλ+B-Ce-λτ2=0,
where
(16)A=D+ax*P′S*,B=a2x*PS*P′S*,C=γD1ax*P′S*.
It is not hard to see that A>0 and B>C>0.We know from Theorem3 that the E* is locally asymptotically stable if set τ2=0 in 15, while when τ2 increases from zero, the equilibrium could lose its stability if pure imaginary roots of 13 appear. Therefore it is convenient to look for pure imaginary roots of 15.Putλ=iω, ω≠0∈R into 15 and separate the real and the imaginary parts of the equation; then we have
(17)cos⁡(ωτ2)=B-ω2C,sin(ωτ2)=-AωC.
Taking the sum of squares in 17 yields
(18)ω4+(A2-2B)ω2+B2-C2=0.
It has real solutions
(19)ω±2=122B-A2±2B-A22-4(B2-C2),
provided
(20)A2-2B2>4(B2-C2)>0>A2-2B.
Notice that if ω is a solution of 17, so is -ω. Accordingly 15 has two pure imaginary solutions λ±=iω± with ω+>ω->0.Now differentiate both sides of15 with respect to τ2,
(21)(2λ+A+Cτ2e-λτ2)dλdτ2=-Cλe-λτ2,
and then
(22)dλdτ2-1=-2λ+Aeλτ2Cλ-τ2λ,eλτ2=Cλ2+Aλ+B.
Therefore,
(23)sign⁡dReλdτ2λ=iω=sign⁡Redλdτ2-1λ=iω=sign⁡Re-2λ+Aλλ2+Aλ+Bλ=iω=sign⁡A2-2B-ω2A2ω2+B-ω22=sign⁡A2-2B+2ω2.
Clearly the sign of d(Reλ)/dτ2 is positive at ω+ or is negative at ω-.Furthermore, we see that there are two sets ofτ2 satisfying 17:
(24)τ2,1n=θ1ω++2nπω+,
where 0≤θ1<2π, n=1,2,3,…, and
(25)cos⁡θ1=B-ω+2C,sinθ1=-Aω+C,τ2,2n=θ2ω-+2nπω-,
where 0≤θ2<2π, and
(26)cos⁡θ2=B-ω-2C,sinθ2=-Aω-C.
To sum up, we have the following.Lemma 4.
Letτ2,i* be an element of the sequences {τ2,in}, i=1,2. For τ2=τ2,i*, the characteristic equation 15 has a pair of simple conjugate pure imaginary roots ±iω±, satisfying
(27)sign⁡d(Reλ)dτ2τ=τ2,i*=sign⁡{A2-2B+2ω±2}.From the above analysis, the stability of the positive equilibriumE* of system 1 for τ1=0 is obvious.Theorem 5.
Assume that conditions9 and 20 hold and τ1=0. Let
(28)τ2*=min⁡i=1,2;n∈N{τ2,in}.
Then the equilibrium E* of system 1 is locally asymptotically stable when τ<τ2* and a Hopf bifurcation occurs at E* when τ=τ2* if and only if sign⁡{A2-2B+2ω±2}>0.


## 3.3. The Caseτ1,τ2>0

The caseτ1,τ2>0 can be regarded as a general one of that in the last section, so our discussion here is based on the result in Section 3.2. The asymptotic stability of E* can be derived from Theorem 5. To state it, we first introduce a lemma about the sign of the real parts of roots of the characteristic equation 13.Lemma 6.
If all roots of15 have negative real parts for τ2>0, then there exists a τ1*(τ2)>0, such that all roots of 13 have negative real parts for τ1∈[0,τ1*(τ2)).Proof.
Suppose that15 has no roots with nonnegative real parts for τ2>0; that is, all the roots of 13 have negative real parts for τ1=0 and τ2>0. According to Theorem  2.1 in [21], on the left side of 13, the sum of the multiplicity of zeros may change if a zero locates on or passes through the imaginary axis in the open right half-plane when τ1 varies.
Recall that13 with τ1=0 has no root with nonnegative real part; then there must exist a τ1*(τ2)>0 such that all roots of 13 have negative real parts for τ1<τ1*(τ2).Theorem 7.
Assume that conditions9 and 20 hold. Then for any τ2∈[0,τ2*), where τ2* is defined in Theorem 5, there exists a τ1*(τ2)>0 such that the positive equilibrium E* of system 1 is locally asymptotically stable when τ1∈[0,τ1*(τ2)).Proof.
From Theorem5, it follows that all roots of 15 have negative real parts when τ2∈[0,τ2*). Then Lemma 6 indicates the validity of the conclusion.


## 4. Numerical Simulations

In this section, we give some numerical tests to support our findings in the previous section. Besides this, we observe the existence of periodic solutions arising from the Hopf bifurcation.In system1, we assume that the nutrient uptake function P(S) is of Michaelis-Menten type
(29)P(S)=Sk+S,
where k>0 is the half-saturation constant [25, 26]. We take the parameters D=0.01, D1=0.5, γ=0.8, k=0.1, a=1, and S0=1. Figure 1 shows the existence and the stability of E* for the case τ1=τ2=0.Figure 1
Asymptotic stability. The steady stateE* is locally asymptotically stable. The parameters are D=0.01, D1=0.5, γ=0.8, k=0.1, a=1, and S0=1. The time delays are τ1=0 and τ2=0.Forτ1=0 and τ2<τ2*, Figure 2 implies the asymptotical stability of E* as stated in Theorem 5. Moreover there also exists the phenomenon of damped oscillations about the equilibrium.Figure 2
Asymptotic stability. The steady stateE* of 1 is locally asymptotically stable, even though damped oscillations can be observed. The parameters are D=0.01, D1=0.5, γ=0.8, k=0.1, a=1, and S0=1. The time delays are τ1=0 and τ2=10.Notice that, in Theorem5, a Hopf bifurcation exists for τ1=0 and τ2=τ2*. It is natural to find that, for τ1=0 and τ2>τ2*, the oscillating solutions emerge; please see Figure 3. Besides, we can see that the periodic oscillation is stable.Figure 3
Oscillation solutions withτ1=0. Parameters are D=0.01, D1=0.5, γ=0.8, k=0.1, a=1, and S0=1. The time delay is τ2=12.217.Meanwhile, we can still expect to obtain the existence of the steady stateE* by taking τ1≤57.8 under condition 9.Figure4 presents the local asymptotic stability of E* for τ1=0.9 and τ2=12.217 which is in accordance with the result in Theorem 7.Figure 4
Asymptotic stability. The steady stateE* of 1 is locally asymptotically stable, even though damped oscillations can be observed. Parameters are D=0.01, D1=0.5, γ=0.8, k=0.1, a=1, and S0=1. The time delays are τ1=0.9 and τ2=12.217.


## 5. Discussion

The most simple idealization for an ecosystem is probably the chemostat, of which the imitation is very close to that in the natural environment such as lakes. But the real environment is more complicated than a chemostat and the continuous culture of microorganisms can be affected by multiple factors, for instance, temperature, climate, food, and so forth. In addition, the nutrient recycling and absorption are also crucial to microorganisms in nature particularly with variable environment. Based on these, many models of chemostat systems have been developed; see [16, 17] for the models with delayed nutrient recycling and [19, 27, 28] with delayed response in growth, respectively.In this paper, we have investigated a chemostat model of a single species with inhibitory exponential substrate uptake and two discrete delays. One indicates the lag caused by gestation; the other models the partially recycled nutrient decomposed by bacteria after the death of biomass.The characteristic equation of the linearized system at the nonzero equilibrium is defined by an exponential polynomial of second degree with delay-dependent coefficients. Because of the existence of two delays, general analytical methods to study the roots of this equation could not be used directly. In the paper, we carry out a rigorous analysis of the stability of the model. Firstly we concentrate on the case when one of the delays,τ1, equals zero and obtains a critical value for the delay τ2: when τ2<τ2* all roots of the characteristic equation have negative parts and when τ2=τ2* purely imaginary roots appear. Then fix the range τ2<τ2*; we show that there exists a number τ1*(τ2) such that all roots of the characteristic equation have negative real parts for τ1∈[0,τ1*(τ2)]. Consequently, the stability results for the chemostat model with two independent delays are obtained and the existence of a Hopf bifurcation that produces periodic solutions is also proven. To test our findings, some numerical simulations are demonstrated which present the periodic oscillations.



---

*Source: 306302-2015-03-03.xml*