# Hopf Bifurcation Analysis for the Model of the Chemostat with One Species of Organism

Authors: Haiyun Bai; Yanhui Zhai
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We research the dynamics of the chemostat model with time delay. The conclusion confirms that a Hopf bifurcation occurs due to the existence of stability switches when the delay varies. By using the normal form theory and center manifold method, we derive the explicit formulas determining the stability and direction of bifurcating periodic solutions. Finally, some numerical simulations are given to illustrate the effectiveness of our results.


---

## Body

## 1. Introduction

Since late 60s, many researchers have been devoted to studying the chemostat, which is considered as an important laboratory set used for breeding microorganism and studying biological systems. In [1] Li et al. propose some new ideals by modifying the basic chemostat model with one species of organism and studying the Hopf bifurcations and stability of the modified one. In [2] Li et al. study the chemostat model with two time delays. They only research the stability of the equilibrium and the existence of the local Hopf bifurcation. However, some subtle mathematical questions on the behavior of solutions of the model are far from completely answered, for example, the bifurcating direction and stability of periodic solutions. Based on this, the main purpose of this study is to provide an insight into these unexplored aspects of the model by using the theory of the center manifold and the normal forms method.Now we consider the basic model of the chemostat with one species of organism (see [3]):
(1)x˙(t)=-(d-p(y(t-τ)))x(t),y˙(t)=(y0-y(t))d-1γp(y(t))x(t),
where x(t) is the concentration of the organism at time t, y(t) is the concentration of the nutrient at time t, p(y)=cy/(k+y) (c, k are positive constants) is the growing rate of y, γ>0 is the ratio of the mass of organism formed and the mass of substrate used, y0>0 is the concentration of the input nutrient, τ>0 is time lag of digestion, and d>0 is flowing rate.


## 2. Stability and Local Hopf Bifurcation

To consider the meaning of the biology, in the section we only focus on investigating the local stability of the interior equilibrium for the system (1). We know that if the equilibrium of system (1) is stable when τ=0 and the characteristic equation of (1) has no purely imaginary roots for any τ>0, it is also stable for any τ>0. On the other hand, if the equilibrium of system (1) is stable when τ=0 and there exist some positive values τ such that the characteristic equation of (1) has a pair of purely imaginary roots, there exists a domain concerning τ such that the equilibrium of system (1) is table in the domain.Wheny(t)<y0 for ∀t>0 and c>d, the system (1) has a unique interior equilibrium. We denote this unique interior equilibrium by (x*,y*).Then it satisfies(2)(d-cy*k+y*)x*=0,(y0-y*)d-cx*y*γ(k+y*)=0.That is(3)(x*,y*)=(γcy0-dy0-dkc-d,dkc-d).Letx1(t)=x(t)-x* and y1(t)=y(t)-y* and still denote x(t), y(t), respectively.Then system (1) becomes
(4)x˙(t)=-d(x(t)+x*)+c(y(t-τ)+y*)(x(t)+x*)k+y(t-τ)+y*,y˙(t)=d(y0-y(t)-y*)-c(y(t)+y*)(x(t)+x*)γ(k+y(t)+y*).The linearization of (4) around (x*,y*) is
(5)x˙(t)=(n-d)x(t)+my(t-τ),y˙(t)=-nγx(t)-(mγ+d)y(t),
where m=kcx*/(k+y*)2 and n=cy*/(k+y*), whose characteristic equation is
(6)λ2+(mγ-n+2d)λ-(n-d)(mγ+d)+mnγe-λτ=0.By setting(m/γ)+d=A, n-d=B and (mn/γ)=K, the characteristic equation (6) can be rewritten as
(7)λ2+(A-B)λ-AB+Ke-λτ=0.Lemma 1.
Whenτ=0 and
(*)A-B>0
are met, the equilibrium (x*,y*) of system (1) is asymptotically stable.Proof.
Whenτ=0, (7) becomes
(8)λ2+(A-B)λ-AB+K=0,
whose characteristic value is
(9)λ1,2=-(A-B)±(A-B)2+4(AB-K)2.
Obviously, when(*) holds, the real parts of λ1,2 are negative.
This completes the proof.In the following, we investigate the distribution of the eigenvalues of the characteristic equation (7).Lemma 2.
Assume that(**)K>AB
is satisfied. Then (7) has a pair of purely imaginary roots ±iω0 when τ=τj, where
(10)ω0=(-(A2+B2)+(A2+B2)2-4(A2B2-K2)1/2)1/2,τj=1ω0[arccosω02+ABK+2jπ],j=0,1,2,….Proof.
Letiω (ω>0) be a root of (7). Then
(11)-ω2+i(A-B)ω-AB+K(cosωτ-isinωτ)=0.
The separation of the real and imaginary parts yields(12)-ω2-AB+Kcosωτ,(A-B)ω-Ksinωτ.
Hence(13)ω2=-(A2+B2)±(A2+B2)2-4(A2B2-K2)2.
Obviously,(**) implies that
(14)ω0=(-(A2+B2)+(A2+B2)2-4(A2B2-K2)2)1/2,
and, hence,
(15)τ0=1ω0arccosω02+ABK.
Defineτj=τ0+(2jπ/ω0), j=0,1,2,…. Then (τj,ω0) solves (12).
This means thatiω0 is a root of (7) when τ=τj, j=0,1,2….
This completes the proof.Lemma2 shows that there exist some positive values τj such that the characteristic equation (7) has a pair of purely imaginary roots.Lemma 3.
Letλ(τ)=a(τ)+ib(τ) be the root of (7) with a(τj)=0 and b(τj)=ω0. When (**) holds, a′(τj)>0.Proof.
By differentiating both sides of (7) with respect to τ, we obtain
(16)2λdλdτ+(A-B)dλdτ-Ke-τλ(λ+τdλdτ)=0.
Then(17)dλdτ=Kλe-τλ2λ+(A-B)-Kτe-τλ.
Substitutingτ=τj and λ=iω0 into (17), we obtain
(18)dλdτ|τ=τj=Kiω0e-τjω0i2iω0+(A-B)-Kτje-τjω0i.
According to (12),
(19)dλdτ|τ=τj=iω0(ω02-i(A-B)ω0-AB)2iω0+(A-B)-τj(ω02-i(A-B)ω0-AB).
Hence(20)a′(τj)=Redλdτ|τ=τj=ω02(2ω02+A2+B2)[A+B-τjω02+τjAB]2+[2ω0+(A+B)τjω0]2>0.
The conclusion is completed.Lemma3 explains that the real parts a′(τ) are monotonously increased in a small neighbourhood concerning τj. In other words, the root of (7) crosses the imaginary axis from the left to the right as τ continuously varies from a number less than τj to one greater than τj.Lemma 4.
When(*) and (**) hold, then there exist τ0<τ1<τ2<⋯ such that all the roots of (7) have negative real parts when τ∈[0,τ0) and (7) has at least one root with positive real parts when τ∈(τj,τj+1), j=1,2,3,…, where τj is defined as in (10).In fact, according to Lemmas1, 2, and 3, it is easy to obtain the results.Applying Lemmas2, 3, and 4 and the Hopf bifurcation theorem (see [4]), we have the following results.Theorem 5.
If(*) and (**) are satisfied, then the equilibrium (x*,y*) is asymptotically stable for τ∈[0,τ0) and unstable for τ>τ0. System (4) undergoes a Hopf bifurcation at (x*,y*) when τ=τj, j=0,1,2…, where τj is defined as in (10).


## 3. Direction and Stability of the Bifurcating Periodic Solutions

Throughout the following section,C([-1,0];R+2) is a phase space, and A stands for an operator, which is different from A in Section 2.In Section3, we will research the stability and direction of the bifurcating periodic solutions of system (1). For convenience, let t=sτ and still denote t. Then the system (5) can be rewritten as (21)x˙(t)=τ(n-d)x(t)+τmy(t-1),y˙(t)=-τnγx(t)-τ(mγ+d)y(t).By using the Taylor series and lettingτ=τ0+μ, we have (22)x˙(t)=(τ0+μ)[kc(k+y*)2kc(k+y*)2(n-d)x(t)+my(t-1)+kc(k+y*)2x(t)y(t-1)+O(x(t)2+y(t-1)2)kc(k+y*)2],y˙(t)=(τ0+μ)[kcγ(k+y*)2-nγx(t)-(mγ+d)y(t)+kcγ(k+y*)2x(t)y(t)+Ox(t)2+y(t)2(mγ+d)(mγ+d)].Clearly,μ=0 is the Hopf bifurcation value of system (21).LetC1=C([-1,0],R+2). For ∀ϕ=(ϕ1,ϕ2)T∈C1, let
(23)Lμ(ϕ)=(τ0+μ)Aϕ(0)+(τ0+μ)Bϕ(-1),(24)f(μ,ϕ)=(τ0+μ)×(kc(k+y*)2ϕ2(-1)ϕ1(0)+O(ϕ22(-1)+ϕ12(0))kcγ(k+y*)2ϕ1(0)ϕ2(0)+O(ϕ12(0)+ϕ22(0))),
where A=(n-d0-(n/γ)-((m/γ)+d)), B=(0m00).By Riesz’s representation theorem, there exists a matrix whose components are bounded variation functionsη(θ,μ) in θ∈[-1,0), such that
(25)Lμ(θ)=∫-10‍d(η(θ,μ)ϕ(θ))
for ϕ∈C.Forϕ∈C([-1,0],R2), we define the operators A and R as (26)A(μ)ϕ(θ)={dϕ(θ)dθ,θ∈[-1,0),∫-10‍d(η(t,μ)ϕ(t)),θ=0,R(μ)ϕ(θ)={0,θ∈[-1,0),f(μ,θ),θ=0.Then the system (21) is equivalent to the following abstract differential equation [5]:
(27)u˙t=A(μ)ut+R(μ)ut,
where u=(u1,u2)T, ut(θ)=u(t+θ) for θ∈[-1,0).As in [6], the bifurcating periodic solutions x(t,μ) of system (21) are indexed by a small parameter ɛ. A solution x(t,μ(ɛ)) has amplitude O(ɛ), period T(ɛ), and  nonzero Floquet exponent β(ɛ) with β(0)=0. Under the present assumptions, μ, T, and β have expansions
(28)μ=μ2ɛ2+μ4ɛ4+⋯,T=2πω(1+T2ɛ2+T4ɛ4+⋯),β=β2ɛ2+β4ɛ4+⋯.The sign ofμ2 determines the direction of bifurcation: if μ2>0(<0), then the Hopf bifurcation is forward (backward). β2 determines the stability of the bifurcating periodic solutions: asymptotically orbitally stable (unstable) if β2<0(>0). And T2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T2>0(<0).Next, we only compute the coefficientsμ2, T2, β2 in these expansions.We define the adjoint operatorA* of A as
(29)A*ψ(s)={-dψ(s)ds,s∈(0,1],∫-10‍d(ηT(t,0)ψ(-t)),s=0.Forϕ∈C[-1,0] and ψ∈C[0,1], define a bilinear form
(30)〈ψ,ϕ〉=ψ-T(0)ϕ(0)-∫-10‍∫0θ‍ψ-T(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0).To determine the normal form of operatorA, we need to calculate the eigenvectors q(θ) and q1*(s) of A and A* corresponding to iτ0ω0 and -iτ0ω0, respectively.Proposition 6.
Assume thatq(θ) and q*(s) are the eigenvectors of A and A* corresponding to iτ0ω0 and -iτ0ω0, respectively, satisfying 〈q*,q〉=1 and 〈q*,q-〉=0.
Then(31)q(θ)=(α,β)Teiτ0ω0θ=(iω0+mγ+d,-nγ)Teiτ0ω0θ,q*(s)=D(α*,β*)Teiτ0ω0s=D(-nγ,iω0-n+d)Teiτ0ω0s,
where D-=1/((m/γ)+2d-n-mτ0e-iτ0ω0).Proof.
Without loss of generality, we just consider the eigenvectorq(θ).
Firstly, whenθ∈[-1,0), by the definition of A and q(θ), we obtain the form q(θ)=(α,β)Teiω0τ0 (here, α, β are unknown parameters).
In what follows, notice thatq(0)=(α,β)T, and Aq(0)=∫-10‍d(η(t,μ)ϕ(t))=iω0τ0q(0). We have
(32)α=iω0+mγ+d,β=-nγ.
Finally, by〈q*,q〉=1, we obtain the parameter D-.
The proof is completed.Now we construct the coordinates of the center manifoldC0 at μ=0.Let(33)z(t)=〈q*,ut〉,W(t,θ)=ut(θ)-2Re{z(t)q(θ)}.On the center manifoldC0, we have
(34)W(t,θ)=W(z(t),z(t)-,θ),
where
(35)W(z,z-,θ)=W20(θ)z22+W11(θ)zz-+W02z-22+W30z36⋯.z and z- are local coordinates for the center manifold C0 in the direction of q and q*, respectively. Since μ=0, we have
(36)z′(t)=iτ0ω0z(t)+〈q*(θ),f(W+2Re{z(t)q(θ)})〉=iτ0ω0z(t)+q*(0)-f(W(z,z-,0)+2Re{z(t)q(0)})≜iτ0ω0z(t)+q*(0)-f0(z,z-),
where
(37)f0(z,z-)=fz2z22+fz-2z-22+fzz-zz-+⋯.We rewrite this as(38)z′(t)=iτ0ω0z+g(z,z-),
with
(39)g(z,z-)=q*-(0)f0(z,z-)=g20z22+g11zz-+g02z-22+g21z2z-2+⋯.By (27) and (36), we obtain
(40)W˙=u˙t-z˙q-z-˙q={AW-2Re{q-*(θ)f0q(θ)},θ∈[-1,0)AW-2Re{q-*(0)f0q(0)}+f0,θ=0≜AW+H(z,z-,θ),
where
(41)H(z,z-,θ)=H20z22+H11zz-+H02z-22+⋯.
By comparing the coefficients of the previous series, we obtain
(42)(A-2iτ0ω0I)W20(θ)=-H20(θ),AW11(θ)=-H11(θ),….Noticingut=(u1t,u2t)=zq(θ)+z-q(θ)-+W(t,θ), it follows that
(43)y(t-1)=βe-iτ0ω0z+β-eiτ0ω0z-+W20(2)(-1)z22+W11(2)(-1)zz-+⋯,u1t(0)=αz+α-z-+W20(1)(0)z22+W11(1)(0)zz-+W02(1)(0)z-22+⋯,u2t(0)=βz+β-z-+W20(2)(0)z22+W11(2)(0)zz-+W02(2)(0)z-22+⋯.Thus, from (39), we have
(44)g20=2D-kcτ0(k+y*)2(αα*βe-iτ0ω0+1γαα*β),g11=D-kcτ0(k+y*)2(1γα-ββ*α-α*βe-iτ0ω0+αα*β-eiτ0ω0+1γαβ*β-+1γα-ββ*),g02=2D-kcτ0(k+y*)2(α-α*β-eiτ0ω0+1γα-β-β*),g21=2D-kcτ0(k+y*)2(12γα*βe-iτ0ω0W11(1)(0)+αα*W11(2)(-1)+α*β-eiτ0ω0W20(1)(0)+12α-α*W20(2)(-1)+1γαβ*W11(2)(0))+12γββ*W20(1)(0)+1γββ*W11(1)(0)+12γα-β*W20(2)(0).Since there areW20, W11 in g21, we still need to compute them.Forθ∈[-1,0),
(45)H(z,z-,θ)=-q*(0)f0q(θ)-q*(0)f-0q-(θ)=-g(z,z-)q(θ)-g(z,z-)q-(θ),
and by comparing coefficients with (41),we obtain
(46)H20(θ)=-g20q(θ)-g02q-(θ),H11(θ)=-g11q(θ)-g-11q-(θ).By substituting these relations into (42), we can derive the following equation:
(47)W20′(θ)=2iτ0ω0W20(θ)+g20q(θ)+g-02q-(θ).By solving forW20(θ), W11(θ), we obtain
(48)W20(θ)=ig20τ0ω0q(0)eiτ0ω0θ+ig-023τ0ω0q-(0)e-iτ0ω0θ+E1e2iτ0ω0θ,W11(θ)=-ig11τ0ω0q(0)eiτ0ω0θ+ig-11τ0ω0q-(0)e-iτ0ω0θ+E2,
where E1=(E1(1),E1(2)), E2=(E2(1),E2(2)).From the definition ofA and (46), we obtain
(49)(2iτ0ω0I-∫-10‍e2iτ0ω0θdη(θ))E1=2kcτ0(k+y*)2(αβe-iω0τ0,1γαβ)T.Thus(50)E1(1)=2kcαβ(k+y*)2(e-iω0τ0(2iω0+mγ+d)+mγe-2iω0τ0)×((2iω0-n+d)(2iω0+mγ+d)+mnγe-2iω0τ0)-1,E1(2)=2kcαβ(k+y*)2(2iω0-n+d-nγe-iω0τ0)×((2iω0-n+d)(2iω0+mγ+d)+mnγe-2iω0τ0)-1.Similarly, we have(51)E2(1)=4kc(k+y*)2(Re{α-βe-iω0τ0}(mγ+d)+Re{α-β}mγe-2iω0τ0)×((d-n)(mγ+d)+mnγe-2iω0τ0)-1,E2(2)=4kc(k+y*)2(1γ(n-d)Re{α-β}+nγRe{α-βe-iω0τ0})×((d-n)(mγ+d)+mnγe-2iω0τ0)-1.Thus, we can compute the parametersW(0) and W(-1).In conclusion, we have computed all the coefficients in (39): g20, g11, g02, and g21.Next, we can compute the following quantities:(52)C1(0)=i2τ0ω0(g20g11-2|g11|-13|g02|2)+g212,μ2=-Re{C1(0)}Re{λ′(τ0)},β2=2Re{C1(0)},T2=-Im{C1(0)}+μ2(Im{λ′(τ0)})ω0.From the discussion in Section2, we know that Re{λ′(τ0)}>0. We therefore have the following result.Theorem 7.
IfRe{C1(0)}<0(>0), the direction of the Hopf bifurcation of the system (1) at the equilibrium (x*,y*) when τ=τ0 is forward (backward) and the bifurcating periodic solutions are orbitally asymptotically stable (unstable).


## 4. Numerical Simulation

In this section, we give a particular example to illustrate the effectiveness of our results. We take the coefficientsy0=0.65, c=0.7, d=0.56, γ=2, k=0.1 in (1). By simple computing, we have the equilibrium (x*,y*)=(0.4,5.4), ω0≐0.2540, τ0≐2.018. Further, we obtain the numerical results directly by means of the software Matlab:(53)g20≐-0.1839+0.0979i,g11≐0.0629+0.0164ig02≐0.0397-0.1067i,g21≐-0.0271+0.0643i.Re{λ′(τ0)}≐1.1591×10-5>0.Re{C1(0)}≐-0.0276<0.From the previous arithmetic, the equilibrium(0.4,5.4) is asymptotically stable when τ=2<τ0=2.018 (see Figure 2). When τ=2.1>τ0=2.018 and Re{C1(0)}≐-0.0276<0, the stable periodic solutions occur from the equilibrium (0.4,5.4) (see Figure 1).Figure 1
Whenτ=2.1, the equilibrium (0.4,5.4) is unstable and periodic solution occurs around them.Figure 2
Whenτ=2, the equilibrium (0.4,5.4) becomes locally stable.Thus, the numerical simulation clarifies the effectiveness of our results.


## 5. Conclusions

In this paper, we have discussed the chemostat model with one species of organism. Firstly, we get the stable domain of equilibrium, and by regarding the delaysτ as the bifurcation parameters and applying the theorem of Hopf bifurcation, we draw the sufficient conditions of the Hopf bifurcation. Further, by using the center manifold and the normal form method, we research the Hopf bifurcating direction and the stability of the model when τ=τ0. Our analysis indicates that the dynamics of the model of the chemostat with one species of organism can be much more complicated than we may have expected. It is interesting to describe the global dynamics of the model by means of the local properties of the interior equilibrium.



---

*Source: 829045-2013-04-17.xml*