# Stability and Sensitive Analysis of a Model with Delay Quorum Sensing

Authors: Zhonghua Zhang; Yaohong Suo; Juan Zhang
Journal: Abstract and Applied Analysis
 (2015)

---

## Abstract

This paper formulates a delay model characterizing the competition between bacteria and immune system. The center manifold reduction method and the normal form theory due to Faria and Magalhaes are used to compute the normal form of the model, and the stability of two nonhyperbolic equilibria is discussed. Sensitivity analysis suggests that the growth rate of bacteria is the most sensitive parameter of the threshold parameterR0 and should be targeted in the controlling strategies.


---

## Body

## 1. Introduction

Quorum sensing is a process that enables bacteria to communicate using secreted signaling molecules calledautoinducers [1]. It makes bacteria regulate their gene expression collectively and control their behaviors on community scale. Quorum sensing was initially observed in the marine bacteriumVibrio fischeri about 30 years ago [2, 3]. Now, many other species are observed to exhibit quorum sensing behavior, including major human pathogens such asStaphylococcus aureus andPseudomonas aeruginosa. Quorum sensing has received more and more attention (see [4–15] and the references therein) and some models are formulated to investigate its effect on the transmission of disease. Braselton and Waltman [4] formulated the dynamically allocated inhibitor production. Dockery and Keener [5] were devoted to developing and studying an ODE and a PDE mathematical models for quorum sensing inPseudomonas aeruginosa and found that quorum sensing works because of a biochemical switch between two stable steady solutions, one with low levels of autoinducer and one with high levels of autoinducer. Koerber et al. [6] presented a mathematical model for the early stages of the infection process byPseudomonas aeruginosa in burn wounds which accounts for the quorum sensing and the diffusion of signalling molecules in the burn-wound environment, and the effects of important parameters on the dynamic properties of the model are discussed in detail. They gave some sufficient conditions for the global asymptotic stability of two boundary equilibria which, respectively, correspond to the survival of theallelochemical producer species or the susceptible one. Fergola et al. [7] formulated anallelopathic competition model in which a distributed delay term simulates quorum sensing which regulates the delay production process ofallelochemicals, and they proved the unique existence of the positive solution and the stability of biologically meaningful steady-state solutions. Anguige et al. [9] constructed a multiphase mathematical model of quorum sensing in a maturingPseudomonas aeruginosa biofilm to investigate the effect of antiquorum sensing and antibiotic treatments on theexopolysaccharide concentration, signal level, bacterial numbers, and biofilm growth rate. The above articles leave out the immune response to the bacterial invasion. However, the immune status of the hosts has a significant impact on the transmission of an infection in a population. Literatures [11–15] employed a quadratic function to describe the quorum sensing of bacteria and formulate some models to characterize the competition between bacteria and the immune system, in which the existence of periodical solution, chaotic motion, and subharmonic bifurcation, the properties of Hopf bifurcation, and the stability of equilibrium et al. were investigated. As a novelty of this paper, a cubic function is used to express quorum sensing. This makes the model own more general nonlinearity, which results in much wider set of outcomes including the coexistence of multiple positive equilibria and the existence of critical equilibrium with simple zero singularity.


## 2. Model Formulation

We denote byXU(t) the concentration of the uninfected target cells, XI(t) the concentration of the infected target cells, B(t) the concentration of the bacteria, IR(t) the concentration of the innate cells, and IA(t) the concentration of the adaptive cells. The dynamic relations among them are as follows: the uninfected target cells have a natural turnover SU and a half-life μXU and they are infected by bacteria with mass-action term α1XUB; the infected target cells are cleared by half-life μXI or adaptive immune cells with mass action term α2XIIA; both the innate and the adaptive immune cells have a source term and a half-life time; for the innate immunity, the source term SIR includes a wide range of cells involved in the first wave of defense of the host such as natural killer cells,polymorphonuclear cells,macrophages, anddendritic cells, and for the adaptive immunity, the source term SIA represents the memory cells, derived from a previous infection or vaccination, a zero source means the first infection with this pathogen and there are no memory cells; both of the two kinds of cells are increased by the signals captured by the bacteria load; the bacteria population has a net growth term represented by a logistic function α20B(1-B/σ) and it is cleared by the innate immunity with mass action term α3BIR. Here, we use a function(1)α20B0Bt2B(t-τ)to formulize the bacteria that compete with the immune cells at time t, which receive signal molecules τ time units ago. B0 is a positive constant, α20 is the growth rate of bacteria, and σ is the effective carrying capacity of the environment. Consequently, the vital dynamics are governed by(2)dBtdt=α20Bt1+B2t-τB0-Btσ-α3BtIRt-α4BtIAt,dXU(t)dt=SU-α1XU(t)B(t)-μXUXU(t),dXI(t)dt=α1XU(t)B(t)-α2IA(t)XI(t)-μXIXI(t),dIR(t)dt=SIR+β1B(t)-μIRIR(t),dIA(t)dt=SIA+β2B(t)-μIAIA(t).Remark 1.
The first equation of system (2) suggests that bacteria are controlled and increased by quorum sensing except for their net growth and they are cleared by the innate immune cells. The second equation of system (2) characterizes the dynamics of the uninfected target cells, and the third one reflects the dynamics of the infected target cells. The uninfected target cells are infected by bacteria in mass action law and they have their own constant input flow, and the infected target cells are killed by the adaptive immune cells. The last two equations of system (2) show that each kind of the immune cells has a special source term and their responses are enhanced by the bacteria load. The target cells and the immune cells have their own half-life terms.


## 3. The Existence and Stability of Equilibria

We introduce(3)R0=α3SIRα20μIR+α4SIAα20μIA,R1=4α202(1-R0)B0α20/σ+α3β1/μIR+α4β2/μIA2.Theorem 2.
System (2) always admits a bacteria-free equilibrium E0=(0,SU/μXU,0,SIR/μIR,SIA/μIA). If R0>1, then system (2) admits a unique positive equilibrium E1. If R0=1, then system (2) admits a unique positive equilibrium E2. If R0<1 and R1>1, then system (2) admits no positive equilibrium. If R0<1 and R1=1, then system (2) admits a unique positive equilibrium E4. If R0<1 and R1<1, then system (2) has two positive equilibria E1 and E3.
Specifically,(4)Ej=Bj*,SUα1Bj*+μXU,kkkkklα1SUBj*μIAα2SIA+β2B4*+μXIμIAα1B4*+μXU,kkkkkkBj*,SUα1Bj*+μXU,SIR+β1Bj*μIR,SIA+β2Bj*μIA,B1,3*=B02α20α20σ+α3β1μIR+α4β2μIA1±1-R1,kkkkkkkkkklkkllB2*=B0α20α20σ+α3β1μIR+α4β2μIA,kkkkkkkkkkklkB4*=B02α20α20σ+α3β1μIR+α4β2μIA,kkkkkkkkkkkkkkkkkkkkkkkklllllllkkj=1,2,3,4.Theorem 3.
IfR0>1, the bacteria-free equilibrium E0 is asymptotically stable, while if R0<1, E0 is unstable.The proofs for Theorems2 and 3 are trivial, so omit them.In the sequel, we study the stability of the positive equilibriumEj. First of all, we transfer it to the origin and get(5)dx(t)dt=Mx(t)+Nx(t-τ)+f(x(t)),where(6)fx=f1x-α1x1x2α1x1x2-α2x3x500,N=2α20Bj*2B0000000000000000000000000,M=-α20Bj*σ00-α3Bj*-α4Bj*-α1SUα1Bj*+μXU-μXU-α1Bj*000α1SUα1Bj*+μXUα1Bj*-μXI-α2(SIA+β2Bj*)μIA0Θβ100-μIR0β2000-μIA,f1x=α20Bj*B02x1+x1t-τx1(t-τ)-x1(α3x4+α4x5)-α20x1x1σ-x12t-τB0,Θ=-α1α2SUμIABj*μXU+α1Bj*μXIμIA+α2SIA+β2Bj*.The characteristic equation of (5) at the origin is(7)λ+μXI+α2SIA+β2Bj*μIAλ+μXU+α1Bj*·λ3+p1λ2+p2λ+p3+q1λ2+q2λ+q3e-λτ=0,where(8)p1=μIR+μIA+Bj*α20σ,p2=μIRμIA+μIR+μIABj*α20σ+β1α3Bj*+β2α4Bj*,p3=μIAμIRBj*α20σ+β1α3Bj*μIA+β2α4μIRBj*,q1=-2α20Bj*2B0,q2=-2α20Bj*2(μIA+μIR)B0,q3=-2α20Bj*2μIAμIRB0.Equation (7) has two negative roots -μXI-α2(SIA+β2Bj*)/μIA, -μXU-α1Bj*, and the other roots can be obtained by solving the following equation:(9)λ3+p1λ2+p2λ+p3+(q1λ2+q2λ+q3)e-λτ=0.It can be seen that(10)p3+q3=Bj*B0σB0(α20μIRμIA+β1α3σμIA+β2α4σμIR)-2α20Bj*μIRμIAσkkkkkkkkkkkll-2α20Bj*μIRμIAσ,p3+q3<0,j=1,2;p3+q3>0,j=3;p3+q3=0,j=4.Theorem 4.
Forτ=0, both E1 and E2 are unstable when they exist.Obviously,E3 does not own zero eigenvalue singularity; namely, zero is not a eigenvalue of the Jacobi matrix at E3. By using the Routh-Hurwitz stability criterion, E3 is locally asymptotically stable if p1+q1>0, p2+q2>0, and (p1+q1)(p2+q2)>(p3+q3) are satisfied together.Theorem 5.
Forτ=0, E3 is asymptotically stable if R0<1, R1<1, and 1-R0 are small sufficiently.Proof.
Letk=1-1-R1 and ε=1-R0. From Theorem 2, E3 exists if and only if R0<1 and R1<1, which leads to 0<1-R0<1 and 0<1-R1<1; that is, 0<ε<1 and 0<k<1. By using R1=4α202ε/B0(α20/σ+α3β1/μIR+α4β2/μIA)2, it follows that k approaches 0 if ε does so. Because (11)B3*2B0=B3*σ+α3α20IR*+α4α20IA*-1,IR*=SIR+β1B3*μIR,IA*=SIA+β2B3*μIA,we have(12)p1+q1=2α20ε-α20σ+2β1α3μIR+2β2α4μIAB3*+μIR+μIA,p2+q2=2α20μIR+μIAε+μIRμIAp2+q2=+α20σμIR+μIA+β1α3+β2α4+2β1α3μIAμIR+2β2α4μIRμIAp2+q2=kkk+2β1α3μIAμIR+2β2α4μIRμIAB3*,p3+q3=2α20μIRμIAεp3+q3=-α20σμIRμIA+β1α3μIA+β2α4μRB3*.Note that B3*=(B0/2α20)(α20/σ+α3β1/μIR+α4β2/μIA)k. We have(13)p1+q1p2+q2-p3+q3=-B02α202α20σ(μIR+μIA)+β1α3μIR(3μIR+4μIA)·μIR+μIA+β22α42μIA25μIA+8μIRkkkkkkkkkllk+2α20σ(μIR+μIA)+β1α3μIR(3μIR+4μIA)β2α4μIA·(3μIA+4μIR)ε+α20B0σ2kkkkkkklkll·α20σ(μIR+μIA)+β1α3μIR(4μIR+5μIA)kkkkkkkklkllll+β2α4μIA(4μIA+5μIR)kkkkkkkklll+B0σβ12α32μIR2(5μIR+8μIA)+13β1β2α3α4μIRμIAkkkkkkkkkkkkllkl·μIR+μIA+β22α42μIA25μIA+8μIR·α20σ+α3β1μIR+α4β2μIAk+μIR+μIAε2+α20B02σ2μIR+μIA+β1α3B02σμIR3μIR+4μIA+3β1β2α3α4B0α20μIRμIA(μIR+μIA)kkkkklkl+β2α4B02σμIA(3μIA+4μIR)+μIR+μIA2kkkklkkl+β12α32B0α20μIR2μIR+2μIA+β22α42B0α20μIA2μIA+2μIRkkkkkkll+3β1β2α3α4B0α20μIRμIA(μIR+μIA)ε+(μIR+μIA)μIRμIA.Seen from above formulae, p1+q1, p2+q2, p3+q3, and (p1+q1)(p2+q2)-(p3+q3) are positive when ε is small, which results in the locally asymptotical stability of E3.Clearly, the left side of (9) is continuous in τ and has roots with positive real parts if and only if it has purely imaginary roots. We will determine whether (9) has purely imaginary roots or not, from which we then will be able to get conditions for all eigenvalues to have negative real parts.Denote the eigenvalue of the characteristic equation (5) by λ=ρ(τ)+iω(τ), where ρ(τ), ω(τ) continually depend on the delay τ. Under the same conditions as Theorem 5, we have ρ(0)<0. Since λ is continuous in τ, one still has ρ(τ)<0 and E3 remains stable if τ is sufficiently small. If there exists a positive value τ0 satisfying ρ(τ0)=0, that is, λ=iω(τ0) is a purely imaginary root of (9), then E3 loses its stability and eventually becomes unstable when ρ(τ) becomes positive. On the other hand, if such a ρ(τ0) does not exist E3 is always stable.Obviously, (9) has a purely imaginary root iω, ω>0, if and only if(14)-iω3-p1ω2+ip2ω+p3+(-q1ω2+iq2ω+q3)e-iωτ=0.Separating the real and imaginary parts of (14) and adding up the squares of them lead to(15)z3+A1z2+A2z+A3=0,where z=ω2, A1=p12-2p2-q12,  A2=p22-2p1p3+2q1q3-q22, and A3=p32-q32.It can be verified thatA3<0 and (15) has positive roots. Without loss of generality, one assumes that (15) has three positive roots defined by z1, z2, and z3, respectively. Then, (14) has three positive roots ωj=zj, j=1,2,3. And then, we have(16)cos⁡(ωkτ)=(p1ωk2-p3)(q3-q1ωk2)+q2ωk2(ωk2-p2)q22ωk2+q3-q1ωk22.Thus, if we denote(17)τkj=1ωkarccos⁡p1ωk2-p3q3-q1ωk2+q2ωk2ωk2-p2q22ωk2+q3-q1ωk22kkkkkkkkklkp1ωk2-p3q3-q1ωk2+q2ωk2ωk2-p2q22ωk2+q3-q1ωk22+2jπωk,where k=1,2,3, j=0,1,2,…, then ±iωk is a pair of purely imaginary roots of (9). Define τ0=min⁡⁡τk0,k=1,2,3. We have the following.Theorem 6.
Under the same conditions as Theorem5, E3 is asymptotically stable for τ∈[0,τ0).Theorem 7.
If3ωk4+2A1ωk2+A2≠0, system (2) undergoes Hopf bifurcation at the positive equilibrium E3 when τ=τkj.


## 4. Normal Forms on the Center Manifold

From the discussions in the above section, it can be seen that the Jacobi matrix atE4 has a uniquely simple zero eigenvalue if R1=1 and τ≠(p2+q2)/q3. To determine the dynamic properties of E4, we have to compute the normal forms on the center manifold. The method used is based on the center manifold reduction and normal form theory due to Faria and Magalhaes; see [16, 17].By means ofR1=1, it obtains(18)a2α32+a1α3+a0=0,where(19)a2=α20μIRμIA2B0β12σ2,a1=2α20μIR2μIA2σ(2μIAσα20SIR+B0α20β1μIA+σB0β1β2α4),a0=α20μIR3μIAB0α202μIA2+B0α42β22σ2+2B0α20μIAα4β2σkkkkkkkkkkkk+4α20σ2μIAα4SIA-4α202σ2μIA2.It is seen from (19) that a0<0 if SIA and B0 are small enough, which means there exists a unique positive α3 solving (18). Denote it by d and define α3 as(20)α3=d+μ,where μ is a small parameter. Obviously, R1=1 if μ=0. Next, we transfer E4 to the origin by(21)x1t=Bt-B4*,x2t=XUt-SUα1B4*+μXU,x3(t)=XI(t)-α1SUB4*μIAα2SIA+β2B4*+μXIμIAα1B4*+μXU,x4(t)=IR(t)-SIR+β1B4*μIR,x5(t)=IA(t)-SIA+β2B4*μIA.Normalizing the delay by t→t/τ, denoting xi(τt) by ui(t), and neglecting the higher order terms O(μ2), we have(22)du(t)dt=A(μ)u(t)+B(μ)u(t-1)+f(u(t),μ),where(23)fu,μ=τf1u,μ-α1u1u2α1u1u2-α2u3u500,B(μ)=τ2α20b0B0(b(μ)+b1μ)000000000000000000000000,Aμ=τα20σbμ00-d+μbμ-α4bμ-α1dμ-μXU-α1bμ000α1dμα1bμ-μXI-α2cμ0α2kμβ100-μIR0β2000-μIA,f1u,μ=u12α20b(μ)B0u1(t-1)-(d+μ)u4kkkkk2α20b(μ)B0u1(t-1)-(d+μ)u4-α4u5+α20B0u12(t-1)-α20σu1+α20b(μ)B0u12(t-1),b0=B02α20α20σ+dβ1μIR+α4β2μIA,b1=B0β12α20μIR,b(μ)=b0+μb1,c0=SIA+β2b0μIA,c1=β2b1μIA,cμ=c0+μc1,d0=SUμXU+α1b0,d1=-SUα1b1μXU+α1b02,dμ=d0+μd1,k(μ)=k0+μk1,k0=-α1d0b0μXI+α2c0,k1=d12(α2c1b0(μXU+b0α1)-μXU(μXI+α2c0))b1μXI+α2c02.LetC=C([-1,0],R5) be the Banach space of continuous functions from [-1,0] into R5 with supremum norm. Define zt∈C as zt(θ)=z(t+θ), θ∈[-1,0]. Equation (22) can be written as the functional differential equation(24)z˙(t)=L(μ)(zt)+F(zt,μ),where V is a neighborhood of zero in space of real numbers, L:C×V→R5 is a parameterized family of bounded linear operators, and F:C×V→R5 is a function with F(0,μ)=0, ∂F(0,μ)/∂z=0 for all μ∈R. They have the following respective forms: (25)L(μ)(φ)=L(0)φ+L1(μ)φ,F(φ,μ)=F2(φ,μ)+F3(φ,μ),where(26)L0φ=τ-α20b0σφ10-2α20b0B0φ1-1-b0dφ40-α4b0φ50-α1d0φ10-μXU+α1b0φ20α1d0φ10+α1b0φ20-α2c0+μXIφ30-α1k0φ50β1φ10-μIRφ40β2φ10-μIAφ50,L1μφ=τμ-α20b1σφ10+4α20b0b1B0φ1-1-b0+db1φ40-α4b1φ50-α1d1φ10-α1b1φ20α1d1φ10+α1b1φ20+α2c1φ30-α2kφ5000,F2φ,μ=τφ102α20b0B0μφ1-1-dφ40-α4φ50-α20σφ10+α20b1B0φ12-1-α1φ10φ20α1φ10φ20-α2φ30φ5000,F3φ,μ=τφ102b1α20B0μφ1-1-μφ40+α20B0φ12-1+α20b1B0μφ12-10000.From the Riesz representation theorem the linear mapL can be expressed in integral form as follows:(27)L(μ)(φ)=∫-10dημ(θ)φ(θ),where ημ is a bounded variation matrix-valued function on [-1,0]. In fact, we can define ημ as(28)ημ(θ)=A(μ)δ(θ)+B(μ)δ(θ+1),where δ(·) is the Dirac delta function.LetR5* be the 5-dimensional vector space of row vectors, and denote C*=C([-1,0],R5*). Then, the adjoint bilinear form on C*×C can be defined as (29)ψs,φθ=ψ0φ0-∫-10∫0θψ(ξ-θ)dημ(θ)φ(ξ)dξ,where φ∈C and ψ∈C*.LetA(μ) be the infinitesimal generator of the flow for the linear system(30)z˙(t)=L(μ)ztwith spectrum σ[A(μ)]. The adjoint operator A*(μ) is defined as the infinitesimal generator for the solution operator of the adjoint equation in C*: (31)w˙(t)=-∫-10w(t-θ)dημ(θ).It is well known that the eigenvalues ofA(μ) with zero real parts play an important role in the bifurcation theory of RFDES. Denote A(0) by A0, and let Λ0={λ∈σ(A0)∣Reλ=0}. We have Λ0={0}.Using the formal adjoint theory for FDEs in [18], the phase space C can be decomposed by Λ0 as C=P⊕Q, where P is the generalized eigenspace associated with the eigenvalues in Λ0, Q={φ∈C∣〈ψ,φ〉=0forallψ∈P*}, and the dual space P* is the generalized eigenspace for A*(0) associated with the eigenvalues in Λ0. Assume that Φ and Ψ are the respective dual bases of P and P* and satisfy 〈Ψ(s),Φ(θ)〉=1. We might as well choose Φ and Ψ as follows:(32)Φθ=1,-α1d0μXU+α1b0,α1μXI+α2c0d0-b0d0μXU+α1b0-k0β2μIA,β1μIR,β2μIAkkKKKlll·d0-b0d0μXU+α1b0-k0β2μIA,β1μIR,β2μIAT,kkkkkkkkkkkkkkklkkkkkkkkKKkkkklkl-1≤θ≤0,Ψ(s)=1,0,0,-b0dμIR,-α4b0μIAl,0≤s≤1,where ·T is the transpose of · and l=1/((1-b0dβ1/μIR2-α4b0β2/μIA2)+τ(2α20b0/B0)).LetB=0. Then, the following equations hold simultaneously: (33)Φ˙=ΦB,Ψ˙=-BΨ.As shown in [16, 17], an appropriate phase space for considering normal forms of (24) is the Banach spaceBC of functions from [-1,0] into R5, which are uniformly continuous on [-1,0) and with a jump discontinuity at 0. Then, the elements ofBC have the form φ+X0ρ, where φ∈C, ρ∈R5, and (34)X0(θ)=I,θ=0,0,-1≤θ<0,so thatBC is identified with C×R5 with the norm |φ+X0ρ|=|φ|C+|ρ|R5.Letπ:BC→P denote the projection (35)π(φ+X0ρ)=ΦΨ,φ+Ψ0ρ,φ∈C,ρ∈R5,and then the decomposition C=P⊕Q yields a decomposition ofBC by Λ0 as the topological direct sumBC=P⊕Ker⁡π with the property Q⊂Ker⁡π, where Q is an infinite dimensional complementary subspace of P and C as shown above. Now, we decompose zt∈C1 in (24) as zt=Φx(t)+y, where x(t)∈R and y∈Q1=Q∩C1, C1 is the subset of C consisting of continuously differentiable functions.Next, we rewrite (24) as follows: (36)z˙(t)=L0zt+L1(μ)zt+F(zt,μ),And, then, under the composition zt=Φx(t)+y, (24) can be decomposed as a system of ODEs in R×Ker⁡π as follows:(37)x˙=Bx+f21x,y,μ+f31x,y,μ,y˙=AQ1y+f22(x,y,μ)+f32(x,y,μ),where AQ1=y˙+X0(L0(y)-y˙(0)) is the restriction of A as an operator from Q1 into Ker⁡π, and(38)f21x,y,μ=Ψ0L1μΦx+y+F2Φx+y,μ,f22x,y,μ=I-πX0L1μΦx+y+F2Φx+y,μ,f31x,y,μ=Ψ0F3Φx+y,μ,f32(x,y,μ)=(I-π)X0F3(Φx+y,μ).As for autonomous ODEs inR5, the normal forms are obtained by a recursive process of changes of variables. At a step j, the terms of order j=2 are computed from the terms of the same order and from the terms of lower orders already computed in previous steps. Assuming that steps of orders 2,3,…,j-1, have already been performed leads to(39)x˙=Bx+∑l≥2j-1‍gj1(x,y,μ)+f~j1(x,y,μ)+h.o.t,y˙=AQ1y+∑l≥2j-1‍gj2(x,y,μ)+f~j2(x,y,μ)+h.o.t,where f~j∈(f~j1,f~j2) is the terms of order j in (x,y,μ) after the previous transformations of variables and h.o.t stands for the higher order terms. Following the algorithm of   [16, 17] at step j, using a change of variables of the form (40)(x,y)=(x^,y^)+Uj(x^,μ)≡(x^,y^)+Uj1x^,μ,Uj2x^,μ,where x,x^∈R,y,y^∈Q1, and Uj1:R2→R, Uj2:R2→Q1 are homogeneous polynomials of degree j in x^ and ϵ, after dropping the hats for simplification of notations, (37) can be put into the normal form(41)x˙=Bx+∑j≥2‍gj1x,y,μ,y˙=AQ1y+∑j≥2‍gj2(x,y,μ),where(42)gj1x,y,μ=f~j1x,y,μ-DxUj1x,μBx-BUj1x,gj2(x,y,μ)=f~j2(x,y,μ)-[DxUj2(x,μ)Bx-AQ1Uj2(x)].It can be verified that (24) satisfies nonresonance conditions; see [16, 17]. Then, the locally invariant manifold for (24) tangent to P at zero must be y=0 and the flow on this manifold is given by 1-dimensional ODE (43)x˙=Bx+g21x,0,μ+g31x,0,μ+h.o.t.The nonlinear terms in (41) are in normal form in the classical sense with respect to matrix B. In application, gj1(x,0,μ) usually can be determined by the following procedure.Theorem 8.
Forj≥2, let Mj denote the operator defined in Vj(R2×Ker⁡π), with values in the same place, by(44)Mjh1,h2=Mj1h1,Mj2h2,Mj1h1x,μ=Dxh1x,μBx-Bh1x,μ,(Mj2)(x,μ)=Dxh2(x,μ)Bx-AQ1[h2(x,μ)]with domain D(Mj)=Vj2(R2)×Vj2(Q1). Here, we use the natation Vj2(Y) to denote the space of homogeneous polynomials of degree j in 2 variables (x,ϵ)∈R2, with coefficients in a Banach space Y.According to [16, 17], we derive(45)Uj(x)=Mj-1PI,jf~j(x,0,μ)∈Ker⁡Mjc,and then(46)gj1(x,0,μ)=(I-PI,j)f~j1(x,0,μ)∈ImMj1c,where PI,j=(PI,j1,PI,j2) is the projection of  Vj2(R2)×Vj2(Ker⁡π) on Im(Mj1)×Im(Mj2).SinceB=0, it can be checked that [Im(M21)]c=span⁡{x2,xμ,μ2}. Then, we obtain(47)g21x,0,μ=lμτ4α20b0b1B0-dβ1μIR-b0+db1β1μIR-α4b1β2μIAx+lτ(3b0α20B0-dβ1μIR-α4β2μIA-α20σ+α20b0B0)x2,and further then the normal form of (24) on the invariant local center manifold y=0 is given by(48)x˙=lμτ4α20b0b1B0-dβ1μIR-b0+db1β1μIR-α4b1β2μIAx+lτ3b0α20B0-dβ1μIR-α4β2μIA-α20σx2+f3(x,0,μ)+h.o.t.If R0<1, by using the definitions, we can verify that(49)b0α20B0>dβ1μIR+α4β2μIA+α20σ.Theorem 9.
IfR1=1, the positive equilibrium E4 is unstable for any τ>0.ForR0=1, E0 is also an equilibrium with simple zero singularity. To discuss its stability, we employ the following perturbation form:(50)α20=d¯+μ¯,where d¯=(α3SIR/μIR)+(α4SIA/μIA) and μ¯ is a small parameter. Then, the normal form of (24) near the bacteria-free equilibrium is as follows: (51)x˙=τμ¯x-τd¯σ+α3β1μIR+α4β2μIAx2+h.o.t.Theorem 10.
IfR0=1, the bacteria-free equilibrium is unstable for any τ>0.


## 5. Sensitive Analysis

Sensitivity indices allow us to measure the relative change in a variable when a parameter changes. The normalized forward sensitivity index of a variable to a parameter is the ratio of the relative change in the variable to the relative change in the parameter. When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives. Here, one adopts the following definition as described by Chitnis et al. [19].Definition 11.
The normalized forward sensitivity index of a variable,u, that depends differentiably on a parameter, p, is defined as γpu=(∂u/∂p)×(p/u).To clear bacteria in the body, we must take measures to makeR0>1 hold. For this end, it is important to determine how crucial each parameter is to R0. Table 1 exhibits the analytical sensitivity indices of R0 to the parameters α20, α3, α4, μIR, μIA, SIR, and SIA.Table 1
Sensitivity indices ofR0 to its parameters.
Parameterj
α
20
α
3
α
4
μ
I
R
μ
I
A
S
I
R
S
I
A
γ
j
R
0
-
1
l
1
l
2
-
l
1
-
l
2
l
1
l
2
l
1
=
S
I
R
μ
I
A
α
3
/
(
α
3
S
I
R
μ
I
A
+
α
4
S
I
A
μ
I
R
) and l2=SIAμIRα4/(α3SIRμIA+α4SIAμIR).The negative sign of the sensitivity index forR0 implies that increase in the relevant parameter leads to the decrease in R0. Note that 0<l1,l2<1. The most sensitive parameter is the growth rate α20, which has a high impact on R0 and should be targeted by intervention strategies. To reduce α20 to ensure R0>1, we may develop an inhibitor to control the reproduction of bacteria or to kill bacteria individual, which agrees with an intuitive expectation.To understand howα20 affects R0, we might as well assume α4=kα3, SIA=kSIR, and μIA=kμIR, and then gets R0=(α3SIR/α20μIR)(1+k). Clearly, for k=0, only the innate immune cell competes with bacteria. Figure 1 exhibits the relationship between R0 and α20 with α3=μIR=0.02, SIR=0.04, and k=0,2. It can be seen from the plot that as α20 decreases R0 increases faster when k=2 than the case of k=0, which reveals that vaccination or other strategies adopted to stimulate immunity of the body are beneficial to the clearance of bacteria.Figure 1
Plot ofR0 against the the effective reproductive rate α20∈(0,20] for k=0,2.Figure2 shows the solutions of system (2) starting from different initial values chosen arbitrarily together with α20=0.01,0.5,2, α3=0.02,  α4=0.04, SIR=1, SIA=2, μIR=0.02, μIA=0.04, B0=1, σ=1/3, β1=0.01, β2=0.02, α1=0.01, α2=0.01, μXU=0.02, and μXI=0.04. It is easy to see that the bacteria concentration B(t) tends to zero, and the smaller growth rate leads to the faster convergence.Figure 2
Simulations ofB(t) for α20=2,0.5,0.01 in each figure.
(a)
(b)
(c)


## 6. Conclusions

This paper formulates the competition between bacteria and immune system by DDEs. Then, the qualitative properties of the model are analyzed. Specially, by virtue of the center manifold reduction and normal form theory due to Faria and Magalhaes [16, 17], the normal form of system (2) associated with zero eigenvalue is computed, from which one deduces that the bacteria-free equilibrium E0 and the positive equilibrium E4 are unstable under the conditions of R0=1 and R1=1, respectively. Next, sensitivity analysis and numerical simulations indicate that the effective reproductive rate α20 is the most sensitive parameter to R0. Theorem 3 suggests the strategies target the decrease of the growth rate which can be successful in disease elimination.On the biological viewpoint, the termsα3SIR/μIR and α4SIA/μIA measure the respective strengths of the innate and adaptive immune system defense against the bacterial challenge, while the factor α20 measures the bacteria’s offensive strength. So with (α3SIR)/(α20μIR)+(α4SIA)/(α20μIA), we can compare the strength of the immune system against the bacterial offensive. Thus, Theorems 2 and 3 have the biological explications: in the domain of attraction of E0, bacteria will be cleared if R0>1; that is, the strength of the immune system defense against the bacteria challenge is not weaker than the bacteria’s offensive strength; in the domain of attraction of E3, bacteria coexist with immune cells when R0<1 and the bacterial challenge is weaker than bacteria’s offensive strength.



---

*Source: 964735-2015-03-05.xml*