# Dynamics and Stability Analysis of a Brucellosis Model with Two Discrete Delays

Authors: Paride O. Lolika; Steady Mushayabasa
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

We present a mathematical model for brucellosis transmission that incorporates two discrete delays and culling of infected animals displaying signs of brucellosis infection. The first delay represents the incubation period while the second account for the time needed to detect and cull infectious animals. Feasibility and stability of the model steady states have been determined analytically and numerically. Further, the occurrence of Hopf bifurcation has been established. Overall the findings from the study, both analytical and numerical, suggest that the two delays can destabilize the system and periodic solutions can arise through Hopf bifurcation.


---

## Body

## 1. Introduction

Brucellosis is one of the neglected zoonotic diseases that remains a major public health problem world over, especially in Middle Eastern countries, southern Europe and North Africa, countries in South and Central Asia, sub-Saharan Africa, Mexico, the Caribbean, and countries in South and Central America [1], with an annual occurrence of more than 500 000 cases [2].In animals, brucellosis is usually transmitted through direct contact between a susceptible and an infectious animal or indirectly, i.e., when a susceptible animal ingest contaminated forages or the excrement containing large quantities of bacteria, generally discharged by infected animals [3]. In humans, however, the majority of the infections result from direct or indirect exposure to infected animals or consumption of raw animal products such as unpasteurized milk or cheese [4]. Since human-to-human transmission of the disease is extremely rare [5], the ultimate management of human brucellosis can be achieved through effective control of brucellosis in livestock. Some researchers postulate that eradication of brucellosis in animals can be attained by combining vaccination with test-and-slaughter programs [1].Mathematical models have proved to be essential guiding tools for epidemiologists, biologists, and policy makers. Models can provide solutions to phenomena which are difficult to measure practically. Recently, a number of mathematical models have been proposed to study the spread and control of brucellosis (see, for example, [3, 6–15] and references therein). A limitation of these previous studies, however, is the noninclusion of the time taken before an infectious animal is detected and culled, despite the fact that in many countries where the disease is endemic lack of financial and human resources often results in delay of detection and culling of infectious animals. The size of this delay may play an important role in minimizing the spread of the disease in the community.It is therefore essential to gain a better and more comprehensive understanding of the effects of time delay on brucellosis transmission and control. Prior studies have shown that epidemic models with time delay often exhibit periodic solutions and as a consequence understanding the nature of these periodic outbreaks plays a crucial role in designing policies that can successfully control the disease. In fact, a recent analysis of brucellosis dataset in countries where the disease is endemic has shown that the disease incidences exhibit a strong periodic behavior with mortality and morbidity of the disease concentrated in a few months each year [16, 17]. Understanding the impact of such seasonal variations is crucial in managing the spread of the disease in the community.Our main goal in this study is to explore the dynamics and stability analysis of a brucellosis model with two discrete delays. Hence we formulated a mathematical model that incorporates two discrete delays. The first delay represents the incubation period while the second delay accounts for the time taken to detect and cull infectious animals. In addition, we subdivide the total animal population into classes of susceptible, vaccinated, infectious undetected, and infectious detected animals. In certain situations immediate slaughter of detected animals may not be feasible and more often these animals are isolated from the rest. However, due to lack of financial and human resources, in addition to lack of knowledge and attitude of farmers, isolation of detected animals has not been a successful practice in most developing nations where animal infections are rampant. Thus in our modelling process we assume that a proportion of detected animals that are not immediately culled are also responsible for disease transmission. Utilizing both analytical and numerical results we have demonstrated that the two delays can destabilize the system and lead to Hopf bifurcation.The article is organised as follows. The model description is given in Section2. Analytical and numerical results are given in Sections 3 and 4, respectively. We end with Section 5 of conclusions.


## 2. Mathematical Model

We subdivide the total animal populationN(t) into compartments of susceptible S(t), vaccinated V(t), undetected infectious animals I1(t), and infectious detected and unculled I2(t). Although brucellosis can be transmitted indirectly (environmental transmission), prior studies [6, 8] suggest that indirect transmission plays a relatively small role on the spread of brucellosis, and as such we have ignored this aspect in our study. Brucellosis dynamics in this study are governed by the following autonomous system:(1)dStdt=A-βI1t+1-pI2tSt-μ+σSt+κVt,dVtdt=σSt-γβI1t+1-pI2tVt-μ+κVt,dI1tdt=βI1t-τ1+1-pI2t-τ1St-τ1+γVt-τ1-α+μ+dI1t,dI2tdt=αI1t-τ2-μ+c+dI2t,whereA is the recruitment rate through birth, μ is the natural death rate, β is the disease direct transmission rate, p is fraction of detected animals that have been culled, σ is the vaccination rate, κ is the vaccination waning rate, τ1 represents the latency delay, τ2 is the delay in detection, γ is the modification factor, α is the rate at which animals are detected and quarantined, c is the culling rate of detected animals, and d is disease-induced death.


## 3. Analytical Results

### 3.1. Initial Conditions

The appropriate space for system (1) is X=C([-τ,0],R+4) the Banach space of continuous functions mapping the interval [-τ,0] into R+4 equipped with subnorm where τ=maxτ1,τ2. From the standard results of functional differential equations [18] it follows that, given any initial conditions x0∈X, there exists a unique solution ϕ(t,x0)=(S(t,x0),V(t,x0),I1(t,x0),I2(t,x0)) of system (1), which satisfies ϕ0=x0, the initial conditions are given by(2)Sθ=x01θ,Vθ=x02θ,I1θ=x03θ,I2θ=x04θ,θ∈-τ,0,wherex0=(x01,x02,x03,x04)∈X, with x0i(θ)≥0,(θ∈[-τ,0],i=1,2…,4) and x03(0),x04(0)>0.


### 3.2. The Basic Reproduction Number

By direct calculation, we find that system (1) has a disease-free equilibrium E0, given by E0=(S0,V0,0,0), with (3)S0=Aμ+κμμ+σ+κ,V0=Aσμμ+σ+κ,S0+γV0=Aμ+κ+γσμμ+σ+κ.By utilizing the next-generation matrix method [19], one can deduce that the basic reproduction number of model (1) is(4)R0=βS0+γV0α1-p+k2k1k2,with k1=(μ+α+d), and k2=(μ+c+d).


### 3.3. Stability of the Disease-Free Equilibrium

In this section, we study the local and global stability of the disease-free equilibrium.Theorem 1.
The disease-free equilibriumE0 of model (1) is locally asymptotically stable when R0<1 and unstable when R0>1.Proof.
To study the local stability of the disease-free equilibriumE0, we linearized system (1) about this point and obtained the characteristic equation, given by the following determinant:(5)-μ+σ-λκ-βS0-β1-pS0σ-μ+κ-λ-γβV0-γβ1-pV000S0+γV0βe-λτ1-k1-λβS0+γV01-pe-λτ100αe-λτ2-k2-λ=0,where λ is the eigenvalue.
From (5) the characteristic equation is(6)λ+μλ+σ+κ+μλ+k2λ+k1-βS0+γV0e-λτ1-1-pαβS0+γV0e-λτ1+τ2=0.
Clearly,-μ and -(σ+κ+μ) are eigenvalues and the other two can be determined from the following equation:(7)λ+k2λ+k1-βS0+γV0e-λτ1-1-pαβS0+γV0e-λτ1+τ2=0.Let(8)gλ,τ1,τ2=λ+k2λ+k1-βS0+γV0e-λτ1-1-pαβS0+γV0e-λτ1+τ2.
Through direct calculation one can easily verify thatg(λ,τ1,τ2) is an analytic function and it follows that(9)g0,τ1,τ2=k1k21-R0,gλ,0,0=λ+k2λ+k1-βS0+γV0-1-pαβS0+γV0.
Now we proceed to investigate the distribution of the solutions of (7) in the following cases:(a)
IfR0<1, then g(0,τ1,τ2)>0. Since the derivative gλ′(λ,τ1,τ2)>0 for λ≥0, τ1>0 and τ2>0, (7) has no zero root and positive real roots for all positive τ1 and τ2. Now we assume that the solution of (7) does not have any purely imaginary roots λ=iω, (ω>0) for some τ1>0, τ2=0. Then by computation, ω must be positive real root of(10)ω4+k12+k22-βS0+γV02ω2+k1k22-βS0+γV0k2+α1-p2=0.IfR0<1, (10) has no positive roots. Hence (7) does not have any purely imaginary roots. We can easily see that the roots of g(λ,0,0)=0 all have negative real parts when R0<1. By the implicit function theorem and the continuity of g(λ,τ1,τ2), we know that all roots of (7) have negative real parts for positive τ1 and τ2=0 which implies that E0 is stable.(b)
IfR0=1 then g(0,τ1,τ2)>0. Since the derivative g′(λ,τ1,τ2)>0 for λ≥0, τ1>0 and τ2>0, (7) has a simple zero and no positive root for all positive τ1 and τ2. By the same argument in case (a), we can obtain that all roots of (7) have negative real parts for positive τ1 and τ2=0 except a zero root. Thus E0 is a degenerate equilibrium of codimension and is stable except in one direction.(c)
IfR0>1, then g(0,τ1,τ2)<0. Since we have limλ→∞g(λ,τ1,τ2)=∞ and g′(λ,τ1,τ2)>0 for λ≥0,τ1>0 and τ2>0, (7) has a unique positive real root for all positive τ1 and τ2=0 and E0 is unstable.Theorem 2.
The disease-free equilibrium of model (1) is globally asymptotically stable when R0≤1 and unstable when R0>1.Proof.
We denote byxt the translation of the solution of the system (1); that is, xt=(S(t+θ),V(t+θ),I1(t+θ),I2(t+θ)), where θ∈[-τ,0] and τ=max{τ1,τ2}. We consider the following Lyapunov functional:(11)Uxt=βα1-p+k2k1k2I1t+β1-pk2I2t+βα1-pk2∫t-τ2tI1θdθ+β2α1-p+k2k1k2∫t-τ1tI1θ+1-pI2θSθ+γVθdθ.Taking the derivative of U along the solutions of (1) gives(12)dUxtdt=βα1-p+k2k1k2βI1t-τ1+1-pI2t-τ1St-τ1+γVt-τ1-βα1-p+k2k2I1t+β1-pk2αI1t-τ2-β1-pI2t+βα1-p+k2k1k2βI1t+1-pI2tSt+γVt-βα1-p+k2k1k2βI1t-τ1+1-pI2t-τ1St-τ1+γVt-τ1+β1-pk2αI1t-β1-pk2αI1t-τ2=βα1-p+k2k1k2βSt+γVtI1t+1-pI2t-βI1t+1-pI2t≤ββα1-p+k2S0+γV0k1k2-1I1t+1-pI2t=βR0-1I1t+1-pI2t.Therefore, U˙<0 holds if R0<1. Furthermore, U˙=0 if R0=1. Therefore, the largest invariant set of U˙ is a singleton {E0} such that S(t)=S0, V(t)=V0, I1(t)=I2(t)=0. It follows from LaSalle’s invariance principle [20] that the disease-free equilibrium of system (1) denoted by E0 is globally asymptotically stable whenever R0≤1. This completes the proof of Theorem 2.


### 3.4. Disease Persistence

System (1) is said to be uniformly persistent if there exists a constant η0>0 such that any solution (S(t),V(t),I1(t),I2(t)) of (1) satisfies(13)liminft→∞St≥η0,liminft→∞Vt≥η0,liminft→∞I1t≥η0,liminft→∞I2t≥η0.Now we give a result on the uniform persistence of system (1). To proceed we introduce the following notation and terminology. Denote by P(t),t≥0 the family of solution operators corresponding to (1). ω-limit set ω(x) of x consists of y∈X such that there exists a sequence tn→∞ as n→∞ with P(tn)x→y as n→∞.Theorem 3.
System (1) is uniformly persistent, if it satisfies R0>1.Proof.
LetX0={x0∈X:x03(0)>0,x04(0)>0}, ∂X=X∖X0={x0∈X:x03(0)=0orx04(0)=0} which is relatively closed in X.
In what follows we demonstrate thatX0 is positively invariant for P(t). From the third and fourth equations of (1) we have(14)dI1tdt≥-α+μ+dI1t,dI2tdt≥-μ+c+dI2t.
SinceI1(0,x0)=x03(0)>0, we have I2(0,x0)=x04(0)>0, and it follows from (14) that(15)I1t,x0≥x030.e-α+μ+dt,I2t,x0≥x040.e-μ+c+dt,∀t≥0.
ThusX0 is positively invariant for P(t).
We set(16)M∂=x0∈X:ϕtx0satisfies1andϕtx0∈∂X,∀t≥0.
We claim that(17)M∂=S,V,0,0.
Assumingϕ(t)∈M∂, ∀t≥0; it suffices to show that I1(t)=I2(t)=0, ∀t≥0. If it is not true, then there exists t0>0 such that either (a) I1(t0)>0, I2(t0)=0; or (b) I1(t0)=0, I2(t0)>0. For case (a), from the fourth equation of (1), we have(18)dI2dtt=t0=αI1t0-τ2>0.
Hence there isϵ0>0 such that I2(t)>0, ∀t∈(t0,t0+ϵ0). On the other hand, from I1(t)>0 there exists ϵ1(0<ϵ1<ϵ0) such that I1(t)>0, ∀t∈(t0,t0+ϵ1). Thus we have I1(t)>0, I2(t)>0, ∀t∈(t0,t0+ϵ1) which contradicts the assumption that (S(t),V(t),I1(t),I2(t))∈M∂,∀t≥0. Similarly, we can obtain a contradiction for case (b). This proves the claim (17).
LetF=⋂x∈Fbω(x), where Fb is the global attractor of P(t) restricted to ∂X. We show that F={E0}. In fact, from F⊆M∂ and the first and second equation of (1), we have limt→∞S(t)=S0 and limt→∞V(t)=V0. Thus, {E0} is the isolated invariant set in X.
Next we show thatWS(E0)∩X0=∅. If this is not true, then there exists a solution (St,Vt,I1t,I2t)∈X0 such that(19)limt→∞St=Aμ+κμμ+σ+κ,limt→∞Vt=Aσμμ+σ+κ,limt→∞I1t=0,limt→∞I2t=0,.
For any sufficiently small constantϵ>0, there exists a positive constant T0=T0(ϵ) such that S(t)>S0-ϵ>0, V(t)>V0-ϵ>0, ∀t>T0. For the constant ϵ given above, it follows from the third and fourth equations of (1) that(20)dI1tdt≥βI1t-τ1+1-pI2t-τ1S0-ϵ+γV0-ϵ-α+μ+dI1t,dI2tdt=αI1t-τ2-μ+c+dI2t,t≥T0+τ
IfI1(t),I2(t)→∞, then, by a standard comparison argument and the nonnegativity, the solution (I~1(t),I~2(t)) of the following monotone system(21)dI~1tdt=βI~1t-τ1+1-pI~2t-τ1S0-ϵ+γV0-ϵ-α+μ+dI~1t,dI~2tdt=αI~1t-τ2-μ+c+dI~2t,t≥T0+τwith initial condition I~1(t)=I1(t),I~2(t)=I2(t),∀t∈[T0,T0+τ] converges to (0,0) as well. Thus limt→∞W^(t)=0, where W^(t)>0 is defined by(22)W^t=βα1-p+k2k1k2I~1t+β1-pk2I~2t+βα1-pk2∫t-τ2tI~1ζdζ+β2α1-p+k2S0-ϵ+γV0-ϵk1k2∫t-τ1tI~1ζ+1-pI~2ζdζ.Differentiating W^(t) with respect to time t yields(23)dW^tdt12=β2α1-p+k2S0-ϵ+γV0-ϵk1k2-βI~1t+1-pI~2t.
SinceR0>1, we have β2α(1-p)+k2(S0-ϵ)+γ(V0-ϵ)/k1k2-β>0 for a sufficiently small ϵ. Therefore W^(t) goes to either infinity or a positive number as t→∞, which leads to a contradiction with limt→∞W^(t)=0. Thus we have WS(E0)∩X0=∅. Define m:X→R+ by m(x0)=min{x03(0),x04(0)},∀x0∈X. It is clear that X0=m-1(0,∞) and ∂X=m-1(0). Thus by [21] theorem 3 we have liminft→∞(I1(t),I2(t)≥(η1,η1) for some constant η1>0. Let η0=min{η1,ϵ} where ϵ is the constant such that liminft→∞S(t)≥ϵ>0, liminft→∞V(t)≥ϵ>0. We showed that liminft→∞S(t)≥η0, liminft→∞V(t)≥η0, liminft→∞I1(t)≥η0, liminft→∞I2(t)≥η0. This completes the proof of Theorem 3.


### 3.5. Existence of the Endemic Equilibrium

Theorem 4.
IfR0>1, model (1) admits a unique endemic equilibrium.Proof.
The endemic equilibriumE∗=(S∗,V∗,I1∗,I2∗) of model (1) is determined by the following:(24)A-βI1t+1-pI2tSt-μ+σSt+κVt=0,σSt-γβI1t+1-pI2tVt-μ+κVt=0,βI1t+1-pI2tSt+γVt-k1I1t=0,αI1t-k2I2t=0.From the last equation in (24) we have(25)I2=αI1k2.The first two equations in (24) give(26)S=Aγβ1+α1-p/k2I1+μ+κβ1+α1-p/k2I1+μ+σγβ1+α1-p/k2I1+μ+κ-κσ,V=Aσβ1+α1-p/k2I1+μ+σγβ1+α1-p/k2I1+μ+κ-κσ.
ForI1≠0, substituting (25) into the third equation in (24) gives(27)S+γV=k1k2βα1-p+k2.
Substituting (26) into (27) yields(28)FI1=Aγβ1+α1-p/k2I1+μ+κ+γσβ1+α1-p/k2I1+μ+σγβ1+α1-p/k2I1+μ+κ-κσ-k1k2βα1-p+k2.Direct calculations show(29)F′I1=-Aβ21+α1-p/k22γ2βI121+α1-p/k2+2γγσ+μ+κI1+Mβ1+α1-p/k2I1+μ+σγβ1+α1-p/k2I1+μ+κ-κσ2<0,
where(30)M=Aβ1+α1-pk2γσ2κ+μ+γ2σμ+σ+μ+κ2.
then functionF(I1) is monotonic decreasing for I1>0, and then we can define the following function:(31)F0=k1k2βα1-p+k2R0-1.
Therefore, by monotonicity of a functionF(I1) there exists a unique endemic equilibrium E∗=(S∗,V∗,I1∗,I2∗)


### 3.6. Stability of the Endemic Equilibrium

In this section, we will investigate the local and global stability of the endemic equilibrium point.Theorem 5.
The endemic equilibriumE∗ of the system (1) is locally asymptotically stable if R0>1 and conditions (36) are satisfied.Proof.
The characteristics equation of system (1) on the infected equilibrium E∗ is given by the following determinant:(32)-μ+σ-λκ-βS∗-β1-pS∗σ-μ+κ-λ-γβV∗-γβ1-pV∗r31r32r33r3400αe-λτ2-μ+c+d-λ=0,
with(33)r31=βI1∗+1-pI2∗e-λτ1,r32=γβI1∗+1-pI2∗e-λτ1,r33=S∗+γV∗βe-λτ1-α+μ+d-λ,r34=βS∗+γV∗1-pe-λτ1.
After some algebraic manipulations one can show that the characteristic equation has the following form:(34)λ4+a1λ3+a2λ2+a3λ+a4=0with(35)a1=4μ+κ+α+σ+2d+c+βγ+1I1∗+1-pI2∗-S∗+γV∗βe-λτ1,a2=β2S∗+γ2V∗I1∗+1-pI2∗e-λτ1+μμ+κ+σ+γβ2I1∗+1-pI2∗2+βγμ+γσ+κ+μI1∗+1-pI2∗-αβ1-pS∗+γV∗e-λτ1+τ2+μ+α+d-S∗+γV∗βe-λτ12μ+κ+σ+βγ+1I1∗+1-pI2∗+μ+c+d3μ+κ+σ+α+d+βγ+1I1∗+1-pI2∗-βS∗+γV∗e-λτ1,a3=βS∗+γV∗κ+γσ+β2γS∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗βe-λτ1+μS∗+γ2V∗I1∗+1-pI2∗β2e-λτ1+μμ+κ+σμ+α+d-S∗+γV∗βe-λτ1+γβ2I1∗+1-pI2∗2μ+α+d-S∗+γV∗βe-λτ1+μμ+c+dμ+κ+σ+βγμ+γσ+μ+κI1∗+1-pI2∗μ+α+d-S∗+γV∗βe-λτ1+γβ2μ+c+dI1∗+1-pI2∗2+β2μ+c+dS∗+γ2V∗I1∗+1-pI2∗e-λτ1+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗+αβ21-pS∗+γ2V∗I1∗+1-pI2∗e-λτ1+τ2+μ+c+d2μ+κ+σ+βγ+1I1∗+1-pI2∗μ+α+d-βS∗+γV∗e-λτ1-αβ1-pS∗+γV∗2μ+κ+σ+βγ+1I1∗+1-pI2∗e-λτ1+τ2,a4=μ+c+dβS∗+γV∗κ+γσ+γβ2S∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗βe-λτ1+μ+c+dμS∗+γ2V∗I1∗+1-pI2∗β2e-λτ1+μμ+c+dκ+μ+σμ+α+d-S∗+γV∗βe-λτ1+γβ2μ+c+dI1∗+1-pI2∗2μ+α+d-S∗+γV∗βe-λτ1+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗μ+α+d-S∗+γV∗βe-λτ1+αγβ21-pS∗σ+κV∗I1∗+1-pI2∗e-λτ1+τ2+αγβ21-pV∗μ+σ+βI1∗+1-pI2∗I1∗+1-pI2∗e-λτ1+τ2+αβ21-pS∗μ+κ+γβI1∗+1-pI2∗I1∗+1-pI2∗e-λτ1+τ2-αβ1-pS∗+γV∗μκ+μ+σ+βγμ+γσ+κ+μI1∗+1-pI2∗e-λτ1+τ2-αβ1-pS∗+γV∗γβ2I1∗+1-pI2∗2e-λτ1+τ2.
By the Routh-Hurwitz criterion, all roots of the characteristics equation (34) have negative real parts and the endemic equilibrium E∗ of system (1) is locally asymptotically stable if τ1=τ2=0, if and only if ai>0(i=1,2,3,4), M1=a1>0,(36)M2=a1a31a2>0,M3=a1a3a51a2a40a1a3>0,M4=a1a3001a2a400a1a3001a2a4>0.
Now, we wish to explore if there is a possibility of having complex roots with positive real part for(a)τ1>0, τ2=0 and (b)τ1=0, τ2>0. We now proceed to explore the above cases as follows.
(a) Ifτ1>0, τ2=0, then the characteristics equation (34) becomes(37)λ4+a11λ3+a21λ2+a31λ+a41=e-λτ1m11λ3+m21λ2+m31λ+m41,
with(38)a11=4μ+κ+α+σ+2d+c+βγ+1I1∗+1-pI2∗,a21=μμ+κ+σ+γβ2I1∗+1-pI2∗2+βγμ+γσ+κ+μI1∗+1-pI2∗+μ+α+d2μ+κ+σ+βγ+1I1∗+1-pI2∗+μ+c+d3μ+κ+σ+α+d+βγ+1I1∗+1-pI2∗,a31=μμ+κ+σμ+α+d+γβ2I1∗+1-pI2∗2μ+α+d+μμ+c+dμ+κ+σ+βγμ+γσ+μ+κI1∗+1-pI2∗μ+α+d+γβ2μ+c+dI1∗+1-pI2∗2+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗+μ+c+d2μ+κ+σ+βγ+1I1∗+1-pI2∗μ+α+d,a41=μμ+c+dκ+μ+σμ+α+d+γβ2μ+c+dI1∗+1-pI2∗2μ+α+d+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗μ+α+d,m11=S∗+γV∗β,m21=-β2S∗+γ2V∗I1∗+1-pI2∗+αβ1-pS∗+γV∗+S∗+γV∗β2μ+κ+σ+βγ+1I1∗+1-pI2∗+μ+c+dβS∗+γV∗,m31=-βS∗+γV∗κ+γσ+β2γS∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗β-μS∗+γ2V∗I1∗+1-pI2∗β2+μμ+κ+σS∗+γV∗β+γβ2I1∗+1-pI2∗2S∗+γV∗+βγμ+γσ+μ+κI1∗+1-pI2∗S∗+γV∗β-β2μ+c+dS∗+γ2V∗I1∗+1-pI2∗-αβ21-pS∗+γ2V∗I1∗+1-pI2∗+μ+c+d2μ+κ+σ+βγ+1I1∗+1-pI2∗βS∗+γV∗+αβ1-pS∗+γV∗2μ+κ+σ+βγ+1I1∗+1-pI2∗,m41=-μ+c+cβS∗+γV∗κ+γσ+γβ2S∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗β-μ+c+dμS∗+γ2V∗I1∗+1-pI2∗β2+μμ+c+dκ+μ+σS∗+γV∗β+γβ2μ+c+dI1∗+1-pI2∗2S∗+γV∗β+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗S∗+γV∗β-αγβ21-pS∗σ+κV∗I1∗+1-pI2∗-αγβ21-pV∗μ+σ+βI1∗+1-pI2∗I1∗+1-pI2∗-αβ21-pS∗μ+κ+γβI1∗+1-pI2∗I1∗+1-pI2∗Now we need to show that all roots of (37) have negative real parts for all τ1∈(0,τ∗). To do so, we show that (37) does not have any purely imaginary roots for all τ1∈(0,τ∗). We assume that λ=iω with ω>0 is a root of (37). Then ω must satisfy the following system: (39)ω4-a21ω2+a41=m41-m21ω2cos⁡ωτ1+m31ω-m11ω3sin⁡ωτ1,a31ω-a11ω3=m31ω-m11ω3cos⁡ωτ1-m41-m21ω2sin⁡ωτ1.Now, we square both sides of each equation above and add the resulting equations; thus ω must be a positive root of(40)ω8+b1ω6+b2ω4+b3ω2+b4=0,where(41)b1=a112-2a21-m112,b2=a212+2a41-a11a31+m11m31-m212,b3=a312+2m21m41-a21a41-m312,b4=a412-m412.
Letz=ω2; then (40) becomes(42)Fz=z4+b1z3+b2z2+b3z+b4=0.
One can observe that ifbi≥0, (i=1,2,3,4), then (42) has no positive roots. Therefore (37) does not have any purely imaginary roots for all τ1∈(0,τ∗) so that all roots of the characteristic equation (37) have negative real parts and the endemic equilibrium E∗ of (1) is stable for all τ1∈(0,τ∗).
(b) Ifτ2>0, τ1=0 then the characteristics equation (34) becomes(43)λ4+α11λ3+α21λ2+α31λ+α41=e-λτ2n11λ3+n21λ2+n31λ+n41(44)α11=4μ+κ+α+σ+2d+c+βγ+1I1∗+1-pI2∗-S∗+γV∗β,α21=β2S∗+γ2V∗I1∗+1-pI2∗+μμ+κ+σ+γβ2I1∗+1-pI2∗2+βγμ+γσ+κ+μI1∗+1-pI2∗+μ+α+d-S∗+γV∗β2μ+κ+σ+βγ+1I1∗+1-pI2∗+μ+c+d3μ+κ+σ+α+d+βγ+1I1∗+1-pI2∗-βS∗+γV∗,α31=βS∗+γV∗κ+γσ+β2γS∗+γV∗I1∗+1-pI2I1∗+1-pI2∗β+μS∗+γ2V∗I1∗+1-pI2∗β2+μμ+κ+σμ+α+d-S∗+γV∗β+γβ2I1∗+1-pI2∗2μ+α+d-S∗+γV∗β+μμ+c+dμ+κ+σ+βγμ+γσ+μ+κI1∗+1-pI2∗μ+α+d-S∗+γV∗β+γβ2μ+c+dI1∗+1-pI2∗2+β2μ+c+dS∗+γ2V∗I1∗+1-pI2∗+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗+μ+c+d2μ+κ+σ+βγ+1I1∗+1-pI2∗μ+α+d-βS∗+γV∗,α41=μ+c+cβS∗+γV∗κ+γσ+γβ2S∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗β+μ+c+dμS∗+γ2V∗I1∗+1-pI2∗β2+μμ+c+dκ+μ+σμ+α+d-S∗+γV∗β+γβ2μ+c+dI1∗+1-pI2∗2μ+α+d-S∗+γV∗β+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗μ+α+d-S∗+γV∗β,n11=0,n21=αβ1-pS∗+γV∗n31=-αβ21-pS∗+γ2V∗I1∗+1-pI2∗αβ1-pS∗+γV∗2μ+κ+σ+βγ+1I1∗+1-pI2∗,n41=-αγβ21-pS∗σ+κV∗I1∗+1-pI2∗-αγβ21-pV∗μ+σ+βI1∗+1-pI2∗I1∗+1-pI2∗-αβ21-pS∗μ+κ+γβI1∗+1-pI2∗I1∗+1-pI2∗+αβ1-pS∗+γV∗μκ+μ+σ+βγμ+γσ+κ+μI1∗+1-pI2∗+αβ1-pS∗+γV∗γβ2I1∗+1-pI2∗2.
Using the same discussion as in the above case then (43) can be written as(45)hz=z4+c1z3+c2z2+c3z+c4=0
with(46)c1=α112-2α21,c2=α212+2α41-α11α31-n212,c3=α312+2n21n41-α21α41-n312,c4=α412-n412.
It follows that all roots of (43) have negative real parts for τ2(0,τ2∗) when cj≥0, j=1,2,3,4 and this implies that endemic equilibrium is locally asymptotically stable for τ2∈(0,τ2∗). This completes the proof.We now explore the global stability of the endemic equilibrium.Theorem 6.
IfR0>1, then E∗ is globally asymptotically stable.Proof.
Let us consider the Lyapunov function:(47)Wt=W1t+W2t+W3.
Here,(48)W1t=St-S∗-S∗ln⁡StS∗+Vt-V∗-V∗ln⁡VtV∗+I1t-I1∗-I1∗ln⁡I1tI1∗+β1-pS∗+γV∗I2∗αI1∗I2t-I2∗-I2∗ln⁡I2tI2∗,W2t=βS∗I1∗∫0τ1I1t-ωSt-ωS∗I1∗-1-lnI1t-ωSt-ωS∗I1∗dω+βγV∗I1∗∫0τ1I1t-ωVt-ωV∗I1∗-1-lnI1t-ωVt-ωV∗I1∗dω+β1-pS∗I2∗∫0τ1I2t-ωSt-ωS∗I2∗-1-lnI2t-ωSt-ωS∗I2∗dω+βγ1-pV∗I2∗∫0τ1I2t-ωVt-ωV∗I2∗-1-lnI2t-ωVt-ωV∗I2∗dω,W3t=β1-pS∗+γV∗I2∗∫0τ2I1t-ωI1∗-1-lnI1t-ωI1∗dω.The derivatives of W1(t) are given by (49)dW1tdt=1-S∗StdSdt+1-V∗VtdVdt+1-I1∗I1tdI1dt+β1-pS∗+γV∗I2∗αI1∗1-I2∗I2tdI2dt.
Substituting the appropriate differentials from (1) we have(50)dW1tdt=1-S∗SA-βI1t+1-pI2tSt-μ+σSt+κVt+1-V∗VσSt-γβI1t+1-pI2tVt-μ+κVt+1-I1∗I1βI1t-τ1+1-pI2t-τ1St-τ1+γVt-τ1-k1I1t+β1-pS∗+γV∗I2∗αI1∗1-I2∗I2αI1t-τ2-k2I2t.
At endemic equilibrium, we have(51)A=βI1∗+1-pI2∗S∗+μ+σS∗-κV∗,μ+κV∗=σS∗-γβI1∗+1-pI2∗V∗,k1I1∗=βI1∗+1-pI2∗S∗+γV∗,k2I2∗=αI1∗.
Using the above constants, we have(52)dW1tdt=μS∗2-StS∗-S∗St+κV∗2-StS∗.V∗Vt-S∗St.VtV∗+μV∗3-S∗St-VtV∗-StS∗.V∗Vt+βI1∗S∗2-StS∗.I1tI1∗-S∗St+β1-pS∗I2∗2-StS∗.I2tI2∗-S∗St-I1tI1∗+βγV∗I1∗3-S∗St-I1tI1∗.VtV∗-StS∗.V∗Vt+βγ1-pV∗I2∗3-S∗St-StS∗.V∗Vt-I2tI2∗.VtV∗-I1tI1∗+βI1t-τ1St-τ11-I1∗I1t+β1-pI2t-τ1St-τ11-I1∗I1t+βγI1t-τ1Vt-τ11-I1∗I1t+βγ1-pI2t-τ1Vt-τ11-I1∗I1t+β1-pI1t-τ2S∗I2∗I1∗-I2∗I2t.I2∗I1∗+βγ1-pI1t-τ2V∗I2∗I1∗-I2∗I2t.I2∗I1∗.
The derivatives ofW2 are given by(53)
d
W
2
t
d
t
=
β
S
∗
I
1
∗
d
d
t
∫
0
τ
1
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
d
ω
+
β
γ
V
∗
I
1
∗
d
d
t
∫
0
τ
1
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
d
ω
+
β
1
-
p
S
∗
I
2
∗
d
d
t
∫
0
τ
1
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
d
ω
+
β
γ
1
-
p
V
∗
I
2
∗
d
d
t
∫
0
τ
1
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
d
ω
,
=
β
S
∗
I
1
∗
∫
0
τ
1
d
d
t
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
d
ω
+
β
γ
V
∗
I
1
∗
∫
0
τ
1
d
d
t
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
d
ω
+
β
1
-
p
S
∗
I
2
∗
∫
0
τ
1
d
d
t
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
d
ω
+
β
γ
1
-
p
V
∗
I
2
∗
∫
0
τ
1
d
d
t
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
d
ω
,
=
-
β
S
∗
I
1
∗
∫
0
τ
1
d
d
ω
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
d
ω
-
β
γ
V
∗
I
1
∗
∫
0
τ
1
d
d
ω
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
d
ω
-
β
1
-
p
S
∗
I
2
∗
∫
0
τ
1
d
d
ω
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
d
ω
-
β
γ
1
-
p
V
∗
I
2
∗
∫
0
τ
1
d
d
ω
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
d
ω
=
β
S
∗
I
1
∗
I
1
t
S
t
S
∗
I
1
∗
-
I
1
t
-
τ
1
S
t
-
τ
1
S
∗
I
1
∗
+
l
n
I
1
t
-
τ
1
S
t
-
τ
1
I
1
t
S
t
+
β
γ
V
∗
I
1
∗
I
1
t
V
t
V
∗
I
1
∗
-
I
1
t
-
τ
1
V
t
-
τ
1
V
∗
I
1
∗
+
l
n
I
1
t
-
τ
1
V
t
-
τ
1
I
1
t
V
t
+
β
1
-
p
S
∗
I
2
∗
I
2
t
S
t
S
∗
I
2
∗
-
I
2
t
-
τ
1
S
t
-
τ
1
S
∗
I
2
∗
+
l
n
I
2
t
-
τ
1
S
t
-
τ
1
I
2
t
S
t
+
β
γ
1
-
p
V
∗
I
2
∗
I
2
t
V
t
S
∗
I
2
∗
-
I
2
t
-
τ
1
V
t
-
τ
1
V
∗
I
2
∗
+
l
n
I
2
t
-
τ
1
S
t
-
τ
1
I
2
t
V
t
=
β
I
1
∗
S
∗
.
S
t
S
∗
.
I
1
t
I
1
∗
+
β
γ
I
1
∗
V
∗
.
V
t
V
∗
.
I
1
t
I
1
∗
+
β
1
-
p
I
2
∗
S
∗
.
S
t
S
∗
.
I
2
t
I
2
∗
+
β
γ
1
-
p
I
2
∗
V
∗
.
V
t
V
∗
.
I
2
t
I
2
∗
-
β
I
1
t
-
τ
1
S
t
-
τ
1
-
β
γ
I
1
t
-
τ
1
V
t
-
τ
1
-
β
1
-
p
I
2
t
-
τ
1
S
t
-
τ
1
-
β
γ
1
-
p
I
2
t
-
τ
1
V
t
-
τ
1
+
β
S
∗
I
1
∗
ln
⁡
I
1
t
-
τ
1
S
t
-
τ
1
I
1
t
S
t
+
β
γ
V
∗
I
1
∗
ln
⁡
I
1
t
-
τ
1
V
t
-
τ
1
I
1
t
V
t
+
β
1
-
p
S
∗
I
2
∗
ln
⁡
I
2
t
-
τ
1
S
t
-
τ
1
I
2
t
S
t
+
β
γ
1
-
p
V
∗
I
2
∗
ln
⁡
I
2
t
-
τ
1
V
t
-
τ
1
I
2
t
V
t
.
The derivatives ofW3(t) are given by(54)dW3tdt=β1-pS∗+γV∗I2∗ddt∫0τ2I1t-ωI1∗-1-lnI1t-ωI1∗dω=β1-pS∗+γV∗I2∗∫0τ2ddtI1t-ωI1∗-1-lnI1t-ωI1∗dω=-β1-pS∗+γV∗I2∗∫0τ2ddωI1t-ωI1∗-1-lnI1t-ωI1∗dω=β1-pS∗+γV∗I2∗I1tI1∗-I1t-τ2I1∗+ln⁡I1t-τ2I1t=β1-pS∗I2∗.I1tI1∗+βγ1-pV∗I2∗.I1tI1∗-β1-pS∗I2∗I1t-τ2.1I1∗-βγ1-pV∗I2∗I1t-τ2.1I1∗+β1-pS∗I2∗ln⁡I1t-τ2I1t+βγ1-pV∗I2∗ln⁡I1t-τ2I1t.
Combining the derivatives ofWj(t), for j=1,2,3, we have(55)dWtdt=μS∗2-StS∗-S∗St+κV∗2-StV∗S∗Vt-S∗VtStV∗+μV∗3-S∗St-VtV∗-StV∗S∗Vt+βI1∗S∗2-S∗St-St-τ1I1t-τ1S∗I1t+ln⁡I1t-τ1St-τ1I1tSt+β1-pS∗I2∗2-S∗St-St-τ1I2t-τ1I1∗S∗I2∗I1t-I1t-τ2I2∗I1∗I2t+ln⁡I2t-τ1St-τ1I1t-τ2I2tStI1t+βγV∗I1∗3-S∗St-StV∗S∗V-Vt-τ1I1t-τ1V∗I1t+ln⁡I1t-τ1Vt-τ1I1tVt+βγ1-pV∗I2∗3-S∗St-StV∗S∗Vt-Vt-τ1I1∗I2t-τ1V∗I1tI2∗-I1t-τ2I2∗I1∗I2t+lnI2t-τ1Vt-τ1I1t-τ2I2tVtI1t.
Note that(56)2≤StS∗+S∗St,2≤StV∗S∗Vt+S∗VtStV∗,3≤S∗St+VtV∗+StV∗S∗Vtfor all S(t)>0 and V(t)>0, because the arithmetic mean is greater than or equal to the geometric mean. Further, note that G(t)=1-g(t)+ln⁡g(t) is always nonpositive for any function g(t)>0, and g(t)=0 if and only if g(t)=1. Hence, it follows that W(t)≤0 and, consequently, W˙(t)≤0. Moreover, the largest invariant set of W˙(t)=0 is a singleton, where S(t)≡S∗, V(t)≡V∗, I1(t)≡I1∗, and I2(t)≡I2∗. Using LaSalle’s invariance principle [20], we conclude that the endemic equilibrium point E∗ is globally asymptotically stable if R0>1.


### 3.7. Hopf Bifurcation Analysis

In this section we determine criteria for Hopf bifurcation to occur using the time delayτ1 and τ2 as the bifurcation parameters to find the interval in which the infected equilibria is stable and unstable out of the same margins. Now to consider Hopf bifurcation we consider the cases (a)τ1=τ10>0, τ2=0 and (b)τ2=τ20>0 and τ1=0. Our analysis is as follows:(a) Whenτ1=τ10>0 and τ2=0 we need to show that dReλ(τ10)/dτ1>0 differentiating (37) with respect to τ1 we get(57)4λ3+3a11λ2+2a21λ+a31dλdτ1=-τ1e-λτ1m11λ3+m21λ2+m31λ+m41+e-λτ13m11λ2+2m21λ+m31dλdτ1-λeλτ1m11λ3+m21λ2+m31λ+m41.This gives(58)dλdτ1-1=4λ3+3a11λ2+2a21λ+a31-λe-λτ1m11λ3+m21λ2+m31λ+m41+3m11λ2+2m21λ+m31λm11λ3+m21λ2+m31λ+m41-τ1λ=3λ4+2a11λ3+a21λ2-a41-λ2λ4+a11λ3+a21λ2+a31λ+a41+2m11λ3+m21λ2-m41λ2m11λ3+m21λ2+m31λ+m41-τ1λ.Thus,(59)signdReλdτ1λ=iω0=sign⁡Redλdτ1-1λ=iω0=sign⁡Re3λ4+2a11λ3+a21λ2-a41-λ2λ4+a11λ3+a21λ2+a31λ+a41λ=iω0+sign⁡Re2m11λ3+m21λ2-m41λ2m11λ3+m21λ2+m31λ+m41λ=iω0=sign⁡Re3ω04-2a11ω03i-a21ω02-a41ω02ω04-a11ω03i-a21ω02+a31ω0i+a41+sign⁡Re-2m11ω03i-m21ω02-m41-ω02-m11ω03i-m21ω02i+m31ω0i+m41=sign⁡3ω08+2a112-2a21ω06+a212+2a41-a11a31ω04-a412ω02ω04-a21ω02+a412+a31ω0-a11ω032+sign⁡m412-2m112ω06-m212-2m11m11ω04ω02m41-m21ω022+m31ω0-m11ω032=sign⁡3ω08+2a112-m112-2a21ω06+a212-m212+2a41+m11m31-a11a31ω04-a412+m412ω02ω04-a21ω02+a412+a31ω0-a11ω032=sign⁡4ω06+3a112-m112-2a21ω04+2a212-m212+2a41+m11m31-a11a31ω02+k0ω04-a21ω02+a412+a31ω0-a11ω032.with(60)k0=a312-m312+2m21m41-a21a41.Lemma 7.
Suppose thatxi, i=1,2,3,4, are the roots of equation g(x)=x4+ϑ1x3+ϑ2x2+ϑ3x+ϑ4=0(ϑ3<0) and x4 is the largest positive root; then(61)dgxdxx=x4>0.In our case considering F(z)=z4+b1z3+b2z2+b3z+b4=0 defined in (42), and assuming b3<0 and ω02 as the largest positive root, we have(62)dReλdτ1=dFz/dzω04-a21ω02+a412+a31ω0-a11ω032>0.The above analysis can be summarized into the following theorem.Theorem 8.
Suppose that(a)R0>1. If either (b)b4<0 or (c)b4≥0 and b3<0 is satisfied, and ω0 is the largest positive simple root of (42) then the infected equilibrium E∗ of model (1) is locally asymptotically stable when τ1<τ10 and unstable when τ1>τ10 where(63)τ10=1ω0arccos⁡m41-m21ω02ω04-a21ω02+a41+m31ω0-m11ω03a31ω0-a11ω03m41-m21ω022+m31ω0-m11ω032,when τ1=τ10, a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcates from E∗ as τ1 passes through the critical value τ10.(b) Whenτ2=τ20>0 and τ1=0 we also need to show that dReλ(τ20)/dτ2>0 differentiating (43) with respect to τ2 we get(64)4λ3+3α11λ2+2α21λ+α31dλdτ2=-τ2e-λτ2n21λ2+n31λ+n41+e-λτ22n21λ+n31dλdτ2-λeλτ2n31λ+n41.This gives(65)dλdτ2-1=4λ3+3α11λ2+2α21λ+α31-λe-λτ2n21λ2+n31λ+n41+2n21λ+n31λn21λ2+n31λ+n41-τ2λ=3λ4+2α11λ3+α21λ2-α41-λ2λ4+α11λ3+α21λ2+α31λ+α41+n21λ2-n41λ2n21λ2+n31λ+n41-τ2λ.Thus,(66)signdReλdτ2λ=iω0=signRedλdτ2-1λ=iω0=sign⁡Re3λ4+2α11λ3+α21λ2-α41-λ2λ4+α11λ3+α21λ2+α31λ+α41λ=iω0+sign⁡Ren21λ2-n41λ2n21λ2+n31λ+n41λ=iω0=sign⁡Re3ω04-2α11ω03i-α21ω02-α41ω02ω04-α11ω03i-α21ω02+α31ω0i+α41+sign⁡Re-n21ω02-n41-ω02-n21ω02i+n31ω0i+n41=sign⁡3ω08+2α112-2α21ω06+α212+2α41-α11α31ω04-α412ω02ω04-α21ω02+α412+α31ω0-α11ω032+sign⁡n412-n212ω04ω02n41-n21ω022+n312ω02=sign⁡3ω08+2α112-2α21ω06+α212-n212+2α41-α11α31ω04-α412+n412ω02ω04-α21ω02+α412+α31ω0-α11ω032=sign⁡4ω06+3α112-2α21ω04+2α212-n212+2α41-α11α31ω02+α312-n312+2n21n41-α21α41ω04-α21ω02+α412+α31ω0-α11ω032.Lemma 9.
Suppose thatxi, i=1,2,3,4, are the roots of equation g(x)=x4+φ1x3+φ2x2+φ3x+φ4=0(φ3<0) and x4 is the largest positive root; then(67)dgxdxx=x4>0.In our case considering h(z)=z4+c1z3+c2z2+c3z+c4=0 defined in (45), and assuming c3<0 and ω02 as the largest positive root, we have(68)dReλdτ2=dhz/dzω04-α21ω02+α412+α31ω0-α11ω032>0.The above analysis can be summarized into the following theorem.Theorem 10.
Suppose that(a)R0>1. If either (b)c4<0 or (c)c4≥0 and c3<0 is satisfied, and ω0 is the largest positive simple root of (45) then the infected equilibrium E∗ of model (1) is locally asymptotically stable when τ2<τ20 and unstable when τ2>τ20 where(69)τ20=1ω0arccos⁡n41-n21ω02ω04-α21ω02+α41+n31ω0α31ω0-α11ω03n41-n21ω022+n312ω02,when τ2=τ20, a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcates from E∗ as τ2 passes through the critical value τ20.From the analysis above, we can deduce that Hopf bifurcations may arise if conditions in Theorems8 and 10 are satisfied. Thus, the introduction of time delay in system (1) can destabilize the system.


## 3.1. Initial Conditions

The appropriate space for system (1) is X=C([-τ,0],R+4) the Banach space of continuous functions mapping the interval [-τ,0] into R+4 equipped with subnorm where τ=maxτ1,τ2. From the standard results of functional differential equations [18] it follows that, given any initial conditions x0∈X, there exists a unique solution ϕ(t,x0)=(S(t,x0),V(t,x0),I1(t,x0),I2(t,x0)) of system (1), which satisfies ϕ0=x0, the initial conditions are given by(2)Sθ=x01θ,Vθ=x02θ,I1θ=x03θ,I2θ=x04θ,θ∈-τ,0,wherex0=(x01,x02,x03,x04)∈X, with x0i(θ)≥0,(θ∈[-τ,0],i=1,2…,4) and x03(0),x04(0)>0.


## 3.2. The Basic Reproduction Number

By direct calculation, we find that system (1) has a disease-free equilibrium E0, given by E0=(S0,V0,0,0), with (3)S0=Aμ+κμμ+σ+κ,V0=Aσμμ+σ+κ,S0+γV0=Aμ+κ+γσμμ+σ+κ.By utilizing the next-generation matrix method [19], one can deduce that the basic reproduction number of model (1) is(4)R0=βS0+γV0α1-p+k2k1k2,with k1=(μ+α+d), and k2=(μ+c+d).


## 3.3. Stability of the Disease-Free Equilibrium

In this section, we study the local and global stability of the disease-free equilibrium.Theorem 1.
The disease-free equilibriumE0 of model (1) is locally asymptotically stable when R0<1 and unstable when R0>1.Proof.
To study the local stability of the disease-free equilibriumE0, we linearized system (1) about this point and obtained the characteristic equation, given by the following determinant:(5)-μ+σ-λκ-βS0-β1-pS0σ-μ+κ-λ-γβV0-γβ1-pV000S0+γV0βe-λτ1-k1-λβS0+γV01-pe-λτ100αe-λτ2-k2-λ=0,where λ is the eigenvalue.
From (5) the characteristic equation is(6)λ+μλ+σ+κ+μλ+k2λ+k1-βS0+γV0e-λτ1-1-pαβS0+γV0e-λτ1+τ2=0.
Clearly,-μ and -(σ+κ+μ) are eigenvalues and the other two can be determined from the following equation:(7)λ+k2λ+k1-βS0+γV0e-λτ1-1-pαβS0+γV0e-λτ1+τ2=0.Let(8)gλ,τ1,τ2=λ+k2λ+k1-βS0+γV0e-λτ1-1-pαβS0+γV0e-λτ1+τ2.
Through direct calculation one can easily verify thatg(λ,τ1,τ2) is an analytic function and it follows that(9)g0,τ1,τ2=k1k21-R0,gλ,0,0=λ+k2λ+k1-βS0+γV0-1-pαβS0+γV0.
Now we proceed to investigate the distribution of the solutions of (7) in the following cases:(a)
IfR0<1, then g(0,τ1,τ2)>0. Since the derivative gλ′(λ,τ1,τ2)>0 for λ≥0, τ1>0 and τ2>0, (7) has no zero root and positive real roots for all positive τ1 and τ2. Now we assume that the solution of (7) does not have any purely imaginary roots λ=iω, (ω>0) for some τ1>0, τ2=0. Then by computation, ω must be positive real root of(10)ω4+k12+k22-βS0+γV02ω2+k1k22-βS0+γV0k2+α1-p2=0.IfR0<1, (10) has no positive roots. Hence (7) does not have any purely imaginary roots. We can easily see that the roots of g(λ,0,0)=0 all have negative real parts when R0<1. By the implicit function theorem and the continuity of g(λ,τ1,τ2), we know that all roots of (7) have negative real parts for positive τ1 and τ2=0 which implies that E0 is stable.(b)
IfR0=1 then g(0,τ1,τ2)>0. Since the derivative g′(λ,τ1,τ2)>0 for λ≥0, τ1>0 and τ2>0, (7) has a simple zero and no positive root for all positive τ1 and τ2. By the same argument in case (a), we can obtain that all roots of (7) have negative real parts for positive τ1 and τ2=0 except a zero root. Thus E0 is a degenerate equilibrium of codimension and is stable except in one direction.(c)
IfR0>1, then g(0,τ1,τ2)<0. Since we have limλ→∞g(λ,τ1,τ2)=∞ and g′(λ,τ1,τ2)>0 for λ≥0,τ1>0 and τ2>0, (7) has a unique positive real root for all positive τ1 and τ2=0 and E0 is unstable.Theorem 2.
The disease-free equilibrium of model (1) is globally asymptotically stable when R0≤1 and unstable when R0>1.Proof.
We denote byxt the translation of the solution of the system (1); that is, xt=(S(t+θ),V(t+θ),I1(t+θ),I2(t+θ)), where θ∈[-τ,0] and τ=max{τ1,τ2}. We consider the following Lyapunov functional:(11)Uxt=βα1-p+k2k1k2I1t+β1-pk2I2t+βα1-pk2∫t-τ2tI1θdθ+β2α1-p+k2k1k2∫t-τ1tI1θ+1-pI2θSθ+γVθdθ.Taking the derivative of U along the solutions of (1) gives(12)dUxtdt=βα1-p+k2k1k2βI1t-τ1+1-pI2t-τ1St-τ1+γVt-τ1-βα1-p+k2k2I1t+β1-pk2αI1t-τ2-β1-pI2t+βα1-p+k2k1k2βI1t+1-pI2tSt+γVt-βα1-p+k2k1k2βI1t-τ1+1-pI2t-τ1St-τ1+γVt-τ1+β1-pk2αI1t-β1-pk2αI1t-τ2=βα1-p+k2k1k2βSt+γVtI1t+1-pI2t-βI1t+1-pI2t≤ββα1-p+k2S0+γV0k1k2-1I1t+1-pI2t=βR0-1I1t+1-pI2t.Therefore, U˙<0 holds if R0<1. Furthermore, U˙=0 if R0=1. Therefore, the largest invariant set of U˙ is a singleton {E0} such that S(t)=S0, V(t)=V0, I1(t)=I2(t)=0. It follows from LaSalle’s invariance principle [20] that the disease-free equilibrium of system (1) denoted by E0 is globally asymptotically stable whenever R0≤1. This completes the proof of Theorem 2.


## 3.4. Disease Persistence

System (1) is said to be uniformly persistent if there exists a constant η0>0 such that any solution (S(t),V(t),I1(t),I2(t)) of (1) satisfies(13)liminft→∞St≥η0,liminft→∞Vt≥η0,liminft→∞I1t≥η0,liminft→∞I2t≥η0.Now we give a result on the uniform persistence of system (1). To proceed we introduce the following notation and terminology. Denote by P(t),t≥0 the family of solution operators corresponding to (1). ω-limit set ω(x) of x consists of y∈X such that there exists a sequence tn→∞ as n→∞ with P(tn)x→y as n→∞.Theorem 3.
System (1) is uniformly persistent, if it satisfies R0>1.Proof.
LetX0={x0∈X:x03(0)>0,x04(0)>0}, ∂X=X∖X0={x0∈X:x03(0)=0orx04(0)=0} which is relatively closed in X.
In what follows we demonstrate thatX0 is positively invariant for P(t). From the third and fourth equations of (1) we have(14)dI1tdt≥-α+μ+dI1t,dI2tdt≥-μ+c+dI2t.
SinceI1(0,x0)=x03(0)>0, we have I2(0,x0)=x04(0)>0, and it follows from (14) that(15)I1t,x0≥x030.e-α+μ+dt,I2t,x0≥x040.e-μ+c+dt,∀t≥0.
ThusX0 is positively invariant for P(t).
We set(16)M∂=x0∈X:ϕtx0satisfies1andϕtx0∈∂X,∀t≥0.
We claim that(17)M∂=S,V,0,0.
Assumingϕ(t)∈M∂, ∀t≥0; it suffices to show that I1(t)=I2(t)=0, ∀t≥0. If it is not true, then there exists t0>0 such that either (a) I1(t0)>0, I2(t0)=0; or (b) I1(t0)=0, I2(t0)>0. For case (a), from the fourth equation of (1), we have(18)dI2dtt=t0=αI1t0-τ2>0.
Hence there isϵ0>0 such that I2(t)>0, ∀t∈(t0,t0+ϵ0). On the other hand, from I1(t)>0 there exists ϵ1(0<ϵ1<ϵ0) such that I1(t)>0, ∀t∈(t0,t0+ϵ1). Thus we have I1(t)>0, I2(t)>0, ∀t∈(t0,t0+ϵ1) which contradicts the assumption that (S(t),V(t),I1(t),I2(t))∈M∂,∀t≥0. Similarly, we can obtain a contradiction for case (b). This proves the claim (17).
LetF=⋂x∈Fbω(x), where Fb is the global attractor of P(t) restricted to ∂X. We show that F={E0}. In fact, from F⊆M∂ and the first and second equation of (1), we have limt→∞S(t)=S0 and limt→∞V(t)=V0. Thus, {E0} is the isolated invariant set in X.
Next we show thatWS(E0)∩X0=∅. If this is not true, then there exists a solution (St,Vt,I1t,I2t)∈X0 such that(19)limt→∞St=Aμ+κμμ+σ+κ,limt→∞Vt=Aσμμ+σ+κ,limt→∞I1t=0,limt→∞I2t=0,.
For any sufficiently small constantϵ>0, there exists a positive constant T0=T0(ϵ) such that S(t)>S0-ϵ>0, V(t)>V0-ϵ>0, ∀t>T0. For the constant ϵ given above, it follows from the third and fourth equations of (1) that(20)dI1tdt≥βI1t-τ1+1-pI2t-τ1S0-ϵ+γV0-ϵ-α+μ+dI1t,dI2tdt=αI1t-τ2-μ+c+dI2t,t≥T0+τ
IfI1(t),I2(t)→∞, then, by a standard comparison argument and the nonnegativity, the solution (I~1(t),I~2(t)) of the following monotone system(21)dI~1tdt=βI~1t-τ1+1-pI~2t-τ1S0-ϵ+γV0-ϵ-α+μ+dI~1t,dI~2tdt=αI~1t-τ2-μ+c+dI~2t,t≥T0+τwith initial condition I~1(t)=I1(t),I~2(t)=I2(t),∀t∈[T0,T0+τ] converges to (0,0) as well. Thus limt→∞W^(t)=0, where W^(t)>0 is defined by(22)W^t=βα1-p+k2k1k2I~1t+β1-pk2I~2t+βα1-pk2∫t-τ2tI~1ζdζ+β2α1-p+k2S0-ϵ+γV0-ϵk1k2∫t-τ1tI~1ζ+1-pI~2ζdζ.Differentiating W^(t) with respect to time t yields(23)dW^tdt12=β2α1-p+k2S0-ϵ+γV0-ϵk1k2-βI~1t+1-pI~2t.
SinceR0>1, we have β2α(1-p)+k2(S0-ϵ)+γ(V0-ϵ)/k1k2-β>0 for a sufficiently small ϵ. Therefore W^(t) goes to either infinity or a positive number as t→∞, which leads to a contradiction with limt→∞W^(t)=0. Thus we have WS(E0)∩X0=∅. Define m:X→R+ by m(x0)=min{x03(0),x04(0)},∀x0∈X. It is clear that X0=m-1(0,∞) and ∂X=m-1(0). Thus by [21] theorem 3 we have liminft→∞(I1(t),I2(t)≥(η1,η1) for some constant η1>0. Let η0=min{η1,ϵ} where ϵ is the constant such that liminft→∞S(t)≥ϵ>0, liminft→∞V(t)≥ϵ>0. We showed that liminft→∞S(t)≥η0, liminft→∞V(t)≥η0, liminft→∞I1(t)≥η0, liminft→∞I2(t)≥η0. This completes the proof of Theorem 3.


## 3.5. Existence of the Endemic Equilibrium

Theorem 4.
IfR0>1, model (1) admits a unique endemic equilibrium.Proof.
The endemic equilibriumE∗=(S∗,V∗,I1∗,I2∗) of model (1) is determined by the following:(24)A-βI1t+1-pI2tSt-μ+σSt+κVt=0,σSt-γβI1t+1-pI2tVt-μ+κVt=0,βI1t+1-pI2tSt+γVt-k1I1t=0,αI1t-k2I2t=0.From the last equation in (24) we have(25)I2=αI1k2.The first two equations in (24) give(26)S=Aγβ1+α1-p/k2I1+μ+κβ1+α1-p/k2I1+μ+σγβ1+α1-p/k2I1+μ+κ-κσ,V=Aσβ1+α1-p/k2I1+μ+σγβ1+α1-p/k2I1+μ+κ-κσ.
ForI1≠0, substituting (25) into the third equation in (24) gives(27)S+γV=k1k2βα1-p+k2.
Substituting (26) into (27) yields(28)FI1=Aγβ1+α1-p/k2I1+μ+κ+γσβ1+α1-p/k2I1+μ+σγβ1+α1-p/k2I1+μ+κ-κσ-k1k2βα1-p+k2.Direct calculations show(29)F′I1=-Aβ21+α1-p/k22γ2βI121+α1-p/k2+2γγσ+μ+κI1+Mβ1+α1-p/k2I1+μ+σγβ1+α1-p/k2I1+μ+κ-κσ2<0,
where(30)M=Aβ1+α1-pk2γσ2κ+μ+γ2σμ+σ+μ+κ2.
then functionF(I1) is monotonic decreasing for I1>0, and then we can define the following function:(31)F0=k1k2βα1-p+k2R0-1.
Therefore, by monotonicity of a functionF(I1) there exists a unique endemic equilibrium E∗=(S∗,V∗,I1∗,I2∗)


## 3.6. Stability of the Endemic Equilibrium

In this section, we will investigate the local and global stability of the endemic equilibrium point.Theorem 5.
The endemic equilibriumE∗ of the system (1) is locally asymptotically stable if R0>1 and conditions (36) are satisfied.Proof.
The characteristics equation of system (1) on the infected equilibrium E∗ is given by the following determinant:(32)-μ+σ-λκ-βS∗-β1-pS∗σ-μ+κ-λ-γβV∗-γβ1-pV∗r31r32r33r3400αe-λτ2-μ+c+d-λ=0,
with(33)r31=βI1∗+1-pI2∗e-λτ1,r32=γβI1∗+1-pI2∗e-λτ1,r33=S∗+γV∗βe-λτ1-α+μ+d-λ,r34=βS∗+γV∗1-pe-λτ1.
After some algebraic manipulations one can show that the characteristic equation has the following form:(34)λ4+a1λ3+a2λ2+a3λ+a4=0with(35)a1=4μ+κ+α+σ+2d+c+βγ+1I1∗+1-pI2∗-S∗+γV∗βe-λτ1,a2=β2S∗+γ2V∗I1∗+1-pI2∗e-λτ1+μμ+κ+σ+γβ2I1∗+1-pI2∗2+βγμ+γσ+κ+μI1∗+1-pI2∗-αβ1-pS∗+γV∗e-λτ1+τ2+μ+α+d-S∗+γV∗βe-λτ12μ+κ+σ+βγ+1I1∗+1-pI2∗+μ+c+d3μ+κ+σ+α+d+βγ+1I1∗+1-pI2∗-βS∗+γV∗e-λτ1,a3=βS∗+γV∗κ+γσ+β2γS∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗βe-λτ1+μS∗+γ2V∗I1∗+1-pI2∗β2e-λτ1+μμ+κ+σμ+α+d-S∗+γV∗βe-λτ1+γβ2I1∗+1-pI2∗2μ+α+d-S∗+γV∗βe-λτ1+μμ+c+dμ+κ+σ+βγμ+γσ+μ+κI1∗+1-pI2∗μ+α+d-S∗+γV∗βe-λτ1+γβ2μ+c+dI1∗+1-pI2∗2+β2μ+c+dS∗+γ2V∗I1∗+1-pI2∗e-λτ1+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗+αβ21-pS∗+γ2V∗I1∗+1-pI2∗e-λτ1+τ2+μ+c+d2μ+κ+σ+βγ+1I1∗+1-pI2∗μ+α+d-βS∗+γV∗e-λτ1-αβ1-pS∗+γV∗2μ+κ+σ+βγ+1I1∗+1-pI2∗e-λτ1+τ2,a4=μ+c+dβS∗+γV∗κ+γσ+γβ2S∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗βe-λτ1+μ+c+dμS∗+γ2V∗I1∗+1-pI2∗β2e-λτ1+μμ+c+dκ+μ+σμ+α+d-S∗+γV∗βe-λτ1+γβ2μ+c+dI1∗+1-pI2∗2μ+α+d-S∗+γV∗βe-λτ1+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗μ+α+d-S∗+γV∗βe-λτ1+αγβ21-pS∗σ+κV∗I1∗+1-pI2∗e-λτ1+τ2+αγβ21-pV∗μ+σ+βI1∗+1-pI2∗I1∗+1-pI2∗e-λτ1+τ2+αβ21-pS∗μ+κ+γβI1∗+1-pI2∗I1∗+1-pI2∗e-λτ1+τ2-αβ1-pS∗+γV∗μκ+μ+σ+βγμ+γσ+κ+μI1∗+1-pI2∗e-λτ1+τ2-αβ1-pS∗+γV∗γβ2I1∗+1-pI2∗2e-λτ1+τ2.
By the Routh-Hurwitz criterion, all roots of the characteristics equation (34) have negative real parts and the endemic equilibrium E∗ of system (1) is locally asymptotically stable if τ1=τ2=0, if and only if ai>0(i=1,2,3,4), M1=a1>0,(36)M2=a1a31a2>0,M3=a1a3a51a2a40a1a3>0,M4=a1a3001a2a400a1a3001a2a4>0.
Now, we wish to explore if there is a possibility of having complex roots with positive real part for(a)τ1>0, τ2=0 and (b)τ1=0, τ2>0. We now proceed to explore the above cases as follows.
(a) Ifτ1>0, τ2=0, then the characteristics equation (34) becomes(37)λ4+a11λ3+a21λ2+a31λ+a41=e-λτ1m11λ3+m21λ2+m31λ+m41,
with(38)a11=4μ+κ+α+σ+2d+c+βγ+1I1∗+1-pI2∗,a21=μμ+κ+σ+γβ2I1∗+1-pI2∗2+βγμ+γσ+κ+μI1∗+1-pI2∗+μ+α+d2μ+κ+σ+βγ+1I1∗+1-pI2∗+μ+c+d3μ+κ+σ+α+d+βγ+1I1∗+1-pI2∗,a31=μμ+κ+σμ+α+d+γβ2I1∗+1-pI2∗2μ+α+d+μμ+c+dμ+κ+σ+βγμ+γσ+μ+κI1∗+1-pI2∗μ+α+d+γβ2μ+c+dI1∗+1-pI2∗2+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗+μ+c+d2μ+κ+σ+βγ+1I1∗+1-pI2∗μ+α+d,a41=μμ+c+dκ+μ+σμ+α+d+γβ2μ+c+dI1∗+1-pI2∗2μ+α+d+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗μ+α+d,m11=S∗+γV∗β,m21=-β2S∗+γ2V∗I1∗+1-pI2∗+αβ1-pS∗+γV∗+S∗+γV∗β2μ+κ+σ+βγ+1I1∗+1-pI2∗+μ+c+dβS∗+γV∗,m31=-βS∗+γV∗κ+γσ+β2γS∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗β-μS∗+γ2V∗I1∗+1-pI2∗β2+μμ+κ+σS∗+γV∗β+γβ2I1∗+1-pI2∗2S∗+γV∗+βγμ+γσ+μ+κI1∗+1-pI2∗S∗+γV∗β-β2μ+c+dS∗+γ2V∗I1∗+1-pI2∗-αβ21-pS∗+γ2V∗I1∗+1-pI2∗+μ+c+d2μ+κ+σ+βγ+1I1∗+1-pI2∗βS∗+γV∗+αβ1-pS∗+γV∗2μ+κ+σ+βγ+1I1∗+1-pI2∗,m41=-μ+c+cβS∗+γV∗κ+γσ+γβ2S∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗β-μ+c+dμS∗+γ2V∗I1∗+1-pI2∗β2+μμ+c+dκ+μ+σS∗+γV∗β+γβ2μ+c+dI1∗+1-pI2∗2S∗+γV∗β+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗S∗+γV∗β-αγβ21-pS∗σ+κV∗I1∗+1-pI2∗-αγβ21-pV∗μ+σ+βI1∗+1-pI2∗I1∗+1-pI2∗-αβ21-pS∗μ+κ+γβI1∗+1-pI2∗I1∗+1-pI2∗Now we need to show that all roots of (37) have negative real parts for all τ1∈(0,τ∗). To do so, we show that (37) does not have any purely imaginary roots for all τ1∈(0,τ∗). We assume that λ=iω with ω>0 is a root of (37). Then ω must satisfy the following system: (39)ω4-a21ω2+a41=m41-m21ω2cos⁡ωτ1+m31ω-m11ω3sin⁡ωτ1,a31ω-a11ω3=m31ω-m11ω3cos⁡ωτ1-m41-m21ω2sin⁡ωτ1.Now, we square both sides of each equation above and add the resulting equations; thus ω must be a positive root of(40)ω8+b1ω6+b2ω4+b3ω2+b4=0,where(41)b1=a112-2a21-m112,b2=a212+2a41-a11a31+m11m31-m212,b3=a312+2m21m41-a21a41-m312,b4=a412-m412.
Letz=ω2; then (40) becomes(42)Fz=z4+b1z3+b2z2+b3z+b4=0.
One can observe that ifbi≥0, (i=1,2,3,4), then (42) has no positive roots. Therefore (37) does not have any purely imaginary roots for all τ1∈(0,τ∗) so that all roots of the characteristic equation (37) have negative real parts and the endemic equilibrium E∗ of (1) is stable for all τ1∈(0,τ∗).
(b) Ifτ2>0, τ1=0 then the characteristics equation (34) becomes(43)λ4+α11λ3+α21λ2+α31λ+α41=e-λτ2n11λ3+n21λ2+n31λ+n41(44)α11=4μ+κ+α+σ+2d+c+βγ+1I1∗+1-pI2∗-S∗+γV∗β,α21=β2S∗+γ2V∗I1∗+1-pI2∗+μμ+κ+σ+γβ2I1∗+1-pI2∗2+βγμ+γσ+κ+μI1∗+1-pI2∗+μ+α+d-S∗+γV∗β2μ+κ+σ+βγ+1I1∗+1-pI2∗+μ+c+d3μ+κ+σ+α+d+βγ+1I1∗+1-pI2∗-βS∗+γV∗,α31=βS∗+γV∗κ+γσ+β2γS∗+γV∗I1∗+1-pI2I1∗+1-pI2∗β+μS∗+γ2V∗I1∗+1-pI2∗β2+μμ+κ+σμ+α+d-S∗+γV∗β+γβ2I1∗+1-pI2∗2μ+α+d-S∗+γV∗β+μμ+c+dμ+κ+σ+βγμ+γσ+μ+κI1∗+1-pI2∗μ+α+d-S∗+γV∗β+γβ2μ+c+dI1∗+1-pI2∗2+β2μ+c+dS∗+γ2V∗I1∗+1-pI2∗+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗+μ+c+d2μ+κ+σ+βγ+1I1∗+1-pI2∗μ+α+d-βS∗+γV∗,α41=μ+c+cβS∗+γV∗κ+γσ+γβ2S∗+γV∗I1∗+1-pI2∗I1∗+1-pI2∗β+μ+c+dμS∗+γ2V∗I1∗+1-pI2∗β2+μμ+c+dκ+μ+σμ+α+d-S∗+γV∗β+γβ2μ+c+dI1∗+1-pI2∗2μ+α+d-S∗+γV∗β+βμ+c+dγμ+γσ+κ+μI1∗+1-pI2∗μ+α+d-S∗+γV∗β,n11=0,n21=αβ1-pS∗+γV∗n31=-αβ21-pS∗+γ2V∗I1∗+1-pI2∗αβ1-pS∗+γV∗2μ+κ+σ+βγ+1I1∗+1-pI2∗,n41=-αγβ21-pS∗σ+κV∗I1∗+1-pI2∗-αγβ21-pV∗μ+σ+βI1∗+1-pI2∗I1∗+1-pI2∗-αβ21-pS∗μ+κ+γβI1∗+1-pI2∗I1∗+1-pI2∗+αβ1-pS∗+γV∗μκ+μ+σ+βγμ+γσ+κ+μI1∗+1-pI2∗+αβ1-pS∗+γV∗γβ2I1∗+1-pI2∗2.
Using the same discussion as in the above case then (43) can be written as(45)hz=z4+c1z3+c2z2+c3z+c4=0
with(46)c1=α112-2α21,c2=α212+2α41-α11α31-n212,c3=α312+2n21n41-α21α41-n312,c4=α412-n412.
It follows that all roots of (43) have negative real parts for τ2(0,τ2∗) when cj≥0, j=1,2,3,4 and this implies that endemic equilibrium is locally asymptotically stable for τ2∈(0,τ2∗). This completes the proof.We now explore the global stability of the endemic equilibrium.Theorem 6.
IfR0>1, then E∗ is globally asymptotically stable.Proof.
Let us consider the Lyapunov function:(47)Wt=W1t+W2t+W3.
Here,(48)W1t=St-S∗-S∗ln⁡StS∗+Vt-V∗-V∗ln⁡VtV∗+I1t-I1∗-I1∗ln⁡I1tI1∗+β1-pS∗+γV∗I2∗αI1∗I2t-I2∗-I2∗ln⁡I2tI2∗,W2t=βS∗I1∗∫0τ1I1t-ωSt-ωS∗I1∗-1-lnI1t-ωSt-ωS∗I1∗dω+βγV∗I1∗∫0τ1I1t-ωVt-ωV∗I1∗-1-lnI1t-ωVt-ωV∗I1∗dω+β1-pS∗I2∗∫0τ1I2t-ωSt-ωS∗I2∗-1-lnI2t-ωSt-ωS∗I2∗dω+βγ1-pV∗I2∗∫0τ1I2t-ωVt-ωV∗I2∗-1-lnI2t-ωVt-ωV∗I2∗dω,W3t=β1-pS∗+γV∗I2∗∫0τ2I1t-ωI1∗-1-lnI1t-ωI1∗dω.The derivatives of W1(t) are given by (49)dW1tdt=1-S∗StdSdt+1-V∗VtdVdt+1-I1∗I1tdI1dt+β1-pS∗+γV∗I2∗αI1∗1-I2∗I2tdI2dt.
Substituting the appropriate differentials from (1) we have(50)dW1tdt=1-S∗SA-βI1t+1-pI2tSt-μ+σSt+κVt+1-V∗VσSt-γβI1t+1-pI2tVt-μ+κVt+1-I1∗I1βI1t-τ1+1-pI2t-τ1St-τ1+γVt-τ1-k1I1t+β1-pS∗+γV∗I2∗αI1∗1-I2∗I2αI1t-τ2-k2I2t.
At endemic equilibrium, we have(51)A=βI1∗+1-pI2∗S∗+μ+σS∗-κV∗,μ+κV∗=σS∗-γβI1∗+1-pI2∗V∗,k1I1∗=βI1∗+1-pI2∗S∗+γV∗,k2I2∗=αI1∗.
Using the above constants, we have(52)dW1tdt=μS∗2-StS∗-S∗St+κV∗2-StS∗.V∗Vt-S∗St.VtV∗+μV∗3-S∗St-VtV∗-StS∗.V∗Vt+βI1∗S∗2-StS∗.I1tI1∗-S∗St+β1-pS∗I2∗2-StS∗.I2tI2∗-S∗St-I1tI1∗+βγV∗I1∗3-S∗St-I1tI1∗.VtV∗-StS∗.V∗Vt+βγ1-pV∗I2∗3-S∗St-StS∗.V∗Vt-I2tI2∗.VtV∗-I1tI1∗+βI1t-τ1St-τ11-I1∗I1t+β1-pI2t-τ1St-τ11-I1∗I1t+βγI1t-τ1Vt-τ11-I1∗I1t+βγ1-pI2t-τ1Vt-τ11-I1∗I1t+β1-pI1t-τ2S∗I2∗I1∗-I2∗I2t.I2∗I1∗+βγ1-pI1t-τ2V∗I2∗I1∗-I2∗I2t.I2∗I1∗.
The derivatives ofW2 are given by(53)
d
W
2
t
d
t
=
β
S
∗
I
1
∗
d
d
t
∫
0
τ
1
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
d
ω
+
β
γ
V
∗
I
1
∗
d
d
t
∫
0
τ
1
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
d
ω
+
β
1
-
p
S
∗
I
2
∗
d
d
t
∫
0
τ
1
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
d
ω
+
β
γ
1
-
p
V
∗
I
2
∗
d
d
t
∫
0
τ
1
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
d
ω
,
=
β
S
∗
I
1
∗
∫
0
τ
1
d
d
t
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
d
ω
+
β
γ
V
∗
I
1
∗
∫
0
τ
1
d
d
t
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
d
ω
+
β
1
-
p
S
∗
I
2
∗
∫
0
τ
1
d
d
t
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
d
ω
+
β
γ
1
-
p
V
∗
I
2
∗
∫
0
τ
1
d
d
t
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
d
ω
,
=
-
β
S
∗
I
1
∗
∫
0
τ
1
d
d
ω
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
S
t
-
ω
S
∗
I
1
∗
d
ω
-
β
γ
V
∗
I
1
∗
∫
0
τ
1
d
d
ω
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
-
1
-
l
n
I
1
t
-
ω
V
t
-
ω
V
∗
I
1
∗
d
ω
-
β
1
-
p
S
∗
I
2
∗
∫
0
τ
1
d
d
ω
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
S
t
-
ω
S
∗
I
2
∗
d
ω
-
β
γ
1
-
p
V
∗
I
2
∗
∫
0
τ
1
d
d
ω
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
-
1
-
l
n
I
2
t
-
ω
V
t
-
ω
V
∗
I
2
∗
d
ω
=
β
S
∗
I
1
∗
I
1
t
S
t
S
∗
I
1
∗
-
I
1
t
-
τ
1
S
t
-
τ
1
S
∗
I
1
∗
+
l
n
I
1
t
-
τ
1
S
t
-
τ
1
I
1
t
S
t
+
β
γ
V
∗
I
1
∗
I
1
t
V
t
V
∗
I
1
∗
-
I
1
t
-
τ
1
V
t
-
τ
1
V
∗
I
1
∗
+
l
n
I
1
t
-
τ
1
V
t
-
τ
1
I
1
t
V
t
+
β
1
-
p
S
∗
I
2
∗
I
2
t
S
t
S
∗
I
2
∗
-
I
2
t
-
τ
1
S
t
-
τ
1
S
∗
I
2
∗
+
l
n
I
2
t
-
τ
1
S
t
-
τ
1
I
2
t
S
t
+
β
γ
1
-
p
V
∗
I
2
∗
I
2
t
V
t
S
∗
I
2
∗
-
I
2
t
-
τ
1
V
t
-
τ
1
V
∗
I
2
∗
+
l
n
I
2
t
-
τ
1
S
t
-
τ
1
I
2
t
V
t
=
β
I
1
∗
S
∗
.
S
t
S
∗
.
I
1
t
I
1
∗
+
β
γ
I
1
∗
V
∗
.
V
t
V
∗
.
I
1
t
I
1
∗
+
β
1
-
p
I
2
∗
S
∗
.
S
t
S
∗
.
I
2
t
I
2
∗
+
β
γ
1
-
p
I
2
∗
V
∗
.
V
t
V
∗
.
I
2
t
I
2
∗
-
β
I
1
t
-
τ
1
S
t
-
τ
1
-
β
γ
I
1
t
-
τ
1
V
t
-
τ
1
-
β
1
-
p
I
2
t
-
τ
1
S
t
-
τ
1
-
β
γ
1
-
p
I
2
t
-
τ
1
V
t
-
τ
1
+
β
S
∗
I
1
∗
ln
⁡
I
1
t
-
τ
1
S
t
-
τ
1
I
1
t
S
t
+
β
γ
V
∗
I
1
∗
ln
⁡
I
1
t
-
τ
1
V
t
-
τ
1
I
1
t
V
t
+
β
1
-
p
S
∗
I
2
∗
ln
⁡
I
2
t
-
τ
1
S
t
-
τ
1
I
2
t
S
t
+
β
γ
1
-
p
V
∗
I
2
∗
ln
⁡
I
2
t
-
τ
1
V
t
-
τ
1
I
2
t
V
t
.
The derivatives ofW3(t) are given by(54)dW3tdt=β1-pS∗+γV∗I2∗ddt∫0τ2I1t-ωI1∗-1-lnI1t-ωI1∗dω=β1-pS∗+γV∗I2∗∫0τ2ddtI1t-ωI1∗-1-lnI1t-ωI1∗dω=-β1-pS∗+γV∗I2∗∫0τ2ddωI1t-ωI1∗-1-lnI1t-ωI1∗dω=β1-pS∗+γV∗I2∗I1tI1∗-I1t-τ2I1∗+ln⁡I1t-τ2I1t=β1-pS∗I2∗.I1tI1∗+βγ1-pV∗I2∗.I1tI1∗-β1-pS∗I2∗I1t-τ2.1I1∗-βγ1-pV∗I2∗I1t-τ2.1I1∗+β1-pS∗I2∗ln⁡I1t-τ2I1t+βγ1-pV∗I2∗ln⁡I1t-τ2I1t.
Combining the derivatives ofWj(t), for j=1,2,3, we have(55)dWtdt=μS∗2-StS∗-S∗St+κV∗2-StV∗S∗Vt-S∗VtStV∗+μV∗3-S∗St-VtV∗-StV∗S∗Vt+βI1∗S∗2-S∗St-St-τ1I1t-τ1S∗I1t+ln⁡I1t-τ1St-τ1I1tSt+β1-pS∗I2∗2-S∗St-St-τ1I2t-τ1I1∗S∗I2∗I1t-I1t-τ2I2∗I1∗I2t+ln⁡I2t-τ1St-τ1I1t-τ2I2tStI1t+βγV∗I1∗3-S∗St-StV∗S∗V-Vt-τ1I1t-τ1V∗I1t+ln⁡I1t-τ1Vt-τ1I1tVt+βγ1-pV∗I2∗3-S∗St-StV∗S∗Vt-Vt-τ1I1∗I2t-τ1V∗I1tI2∗-I1t-τ2I2∗I1∗I2t+lnI2t-τ1Vt-τ1I1t-τ2I2tVtI1t.
Note that(56)2≤StS∗+S∗St,2≤StV∗S∗Vt+S∗VtStV∗,3≤S∗St+VtV∗+StV∗S∗Vtfor all S(t)>0 and V(t)>0, because the arithmetic mean is greater than or equal to the geometric mean. Further, note that G(t)=1-g(t)+ln⁡g(t) is always nonpositive for any function g(t)>0, and g(t)=0 if and only if g(t)=1. Hence, it follows that W(t)≤0 and, consequently, W˙(t)≤0. Moreover, the largest invariant set of W˙(t)=0 is a singleton, where S(t)≡S∗, V(t)≡V∗, I1(t)≡I1∗, and I2(t)≡I2∗. Using LaSalle’s invariance principle [20], we conclude that the endemic equilibrium point E∗ is globally asymptotically stable if R0>1.


## 3.7. Hopf Bifurcation Analysis

In this section we determine criteria for Hopf bifurcation to occur using the time delayτ1 and τ2 as the bifurcation parameters to find the interval in which the infected equilibria is stable and unstable out of the same margins. Now to consider Hopf bifurcation we consider the cases (a)τ1=τ10>0, τ2=0 and (b)τ2=τ20>0 and τ1=0. Our analysis is as follows:(a) Whenτ1=τ10>0 and τ2=0 we need to show that dReλ(τ10)/dτ1>0 differentiating (37) with respect to τ1 we get(57)4λ3+3a11λ2+2a21λ+a31dλdτ1=-τ1e-λτ1m11λ3+m21λ2+m31λ+m41+e-λτ13m11λ2+2m21λ+m31dλdτ1-λeλτ1m11λ3+m21λ2+m31λ+m41.This gives(58)dλdτ1-1=4λ3+3a11λ2+2a21λ+a31-λe-λτ1m11λ3+m21λ2+m31λ+m41+3m11λ2+2m21λ+m31λm11λ3+m21λ2+m31λ+m41-τ1λ=3λ4+2a11λ3+a21λ2-a41-λ2λ4+a11λ3+a21λ2+a31λ+a41+2m11λ3+m21λ2-m41λ2m11λ3+m21λ2+m31λ+m41-τ1λ.Thus,(59)signdReλdτ1λ=iω0=sign⁡Redλdτ1-1λ=iω0=sign⁡Re3λ4+2a11λ3+a21λ2-a41-λ2λ4+a11λ3+a21λ2+a31λ+a41λ=iω0+sign⁡Re2m11λ3+m21λ2-m41λ2m11λ3+m21λ2+m31λ+m41λ=iω0=sign⁡Re3ω04-2a11ω03i-a21ω02-a41ω02ω04-a11ω03i-a21ω02+a31ω0i+a41+sign⁡Re-2m11ω03i-m21ω02-m41-ω02-m11ω03i-m21ω02i+m31ω0i+m41=sign⁡3ω08+2a112-2a21ω06+a212+2a41-a11a31ω04-a412ω02ω04-a21ω02+a412+a31ω0-a11ω032+sign⁡m412-2m112ω06-m212-2m11m11ω04ω02m41-m21ω022+m31ω0-m11ω032=sign⁡3ω08+2a112-m112-2a21ω06+a212-m212+2a41+m11m31-a11a31ω04-a412+m412ω02ω04-a21ω02+a412+a31ω0-a11ω032=sign⁡4ω06+3a112-m112-2a21ω04+2a212-m212+2a41+m11m31-a11a31ω02+k0ω04-a21ω02+a412+a31ω0-a11ω032.with(60)k0=a312-m312+2m21m41-a21a41.Lemma 7.
Suppose thatxi, i=1,2,3,4, are the roots of equation g(x)=x4+ϑ1x3+ϑ2x2+ϑ3x+ϑ4=0(ϑ3<0) and x4 is the largest positive root; then(61)dgxdxx=x4>0.In our case considering F(z)=z4+b1z3+b2z2+b3z+b4=0 defined in (42), and assuming b3<0 and ω02 as the largest positive root, we have(62)dReλdτ1=dFz/dzω04-a21ω02+a412+a31ω0-a11ω032>0.The above analysis can be summarized into the following theorem.Theorem 8.
Suppose that(a)R0>1. If either (b)b4<0 or (c)b4≥0 and b3<0 is satisfied, and ω0 is the largest positive simple root of (42) then the infected equilibrium E∗ of model (1) is locally asymptotically stable when τ1<τ10 and unstable when τ1>τ10 where(63)τ10=1ω0arccos⁡m41-m21ω02ω04-a21ω02+a41+m31ω0-m11ω03a31ω0-a11ω03m41-m21ω022+m31ω0-m11ω032,when τ1=τ10, a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcates from E∗ as τ1 passes through the critical value τ10.(b) Whenτ2=τ20>0 and τ1=0 we also need to show that dReλ(τ20)/dτ2>0 differentiating (43) with respect to τ2 we get(64)4λ3+3α11λ2+2α21λ+α31dλdτ2=-τ2e-λτ2n21λ2+n31λ+n41+e-λτ22n21λ+n31dλdτ2-λeλτ2n31λ+n41.This gives(65)dλdτ2-1=4λ3+3α11λ2+2α21λ+α31-λe-λτ2n21λ2+n31λ+n41+2n21λ+n31λn21λ2+n31λ+n41-τ2λ=3λ4+2α11λ3+α21λ2-α41-λ2λ4+α11λ3+α21λ2+α31λ+α41+n21λ2-n41λ2n21λ2+n31λ+n41-τ2λ.Thus,(66)signdReλdτ2λ=iω0=signRedλdτ2-1λ=iω0=sign⁡Re3λ4+2α11λ3+α21λ2-α41-λ2λ4+α11λ3+α21λ2+α31λ+α41λ=iω0+sign⁡Ren21λ2-n41λ2n21λ2+n31λ+n41λ=iω0=sign⁡Re3ω04-2α11ω03i-α21ω02-α41ω02ω04-α11ω03i-α21ω02+α31ω0i+α41+sign⁡Re-n21ω02-n41-ω02-n21ω02i+n31ω0i+n41=sign⁡3ω08+2α112-2α21ω06+α212+2α41-α11α31ω04-α412ω02ω04-α21ω02+α412+α31ω0-α11ω032+sign⁡n412-n212ω04ω02n41-n21ω022+n312ω02=sign⁡3ω08+2α112-2α21ω06+α212-n212+2α41-α11α31ω04-α412+n412ω02ω04-α21ω02+α412+α31ω0-α11ω032=sign⁡4ω06+3α112-2α21ω04+2α212-n212+2α41-α11α31ω02+α312-n312+2n21n41-α21α41ω04-α21ω02+α412+α31ω0-α11ω032.Lemma 9.
Suppose thatxi, i=1,2,3,4, are the roots of equation g(x)=x4+φ1x3+φ2x2+φ3x+φ4=0(φ3<0) and x4 is the largest positive root; then(67)dgxdxx=x4>0.In our case considering h(z)=z4+c1z3+c2z2+c3z+c4=0 defined in (45), and assuming c3<0 and ω02 as the largest positive root, we have(68)dReλdτ2=dhz/dzω04-α21ω02+α412+α31ω0-α11ω032>0.The above analysis can be summarized into the following theorem.Theorem 10.
Suppose that(a)R0>1. If either (b)c4<0 or (c)c4≥0 and c3<0 is satisfied, and ω0 is the largest positive simple root of (45) then the infected equilibrium E∗ of model (1) is locally asymptotically stable when τ2<τ20 and unstable when τ2>τ20 where(69)τ20=1ω0arccos⁡n41-n21ω02ω04-α21ω02+α41+n31ω0α31ω0-α11ω03n41-n21ω022+n312ω02,when τ2=τ20, a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcates from E∗ as τ2 passes through the critical value τ20.From the analysis above, we can deduce that Hopf bifurcations may arise if conditions in Theorems8 and 10 are satisfied. Thus, the introduction of time delay in system (1) can destabilize the system.


## 4. Numerical Results

In order to explore the behavior of system (1) and illustrate the stability of the equilibria solutions, we numerically solve system (1) using MATLAB and parameter values adopted from Table 1.Table 1
Model parameters and variables and their baseline values.
Symbol Definition Value Unit Source d Elimination rate due to brucellosis 0.15 year-1 [8] p Fraction of infectious animals culled upon detection 0.5 - [22] β Direct transmission rate 3.844 × 1 0 - 6 animal-1year-1 [6] κ Vaccination waning rate 0.4 year-1 [8] μ Natural elimination rate 0.25 year-1 [6] γ Modification factor 0.18 - [8] A Recruitment rate 76434 animalsyear-1 [6] σ Vaccination rate 0.316 year-1 [8] α Detection rate Varied year-1 - c Culling rate 0.15 year-1 [8]In Figure1 we illustrate the effects of varying the delay (τ1=τ2) on the dynamics of system (1). Figures 1(a) and 1(b) demonstrate that the system approaches the stable disease-free or endemic equilibrium for R0<1 and R0>1, respectively. One should note that, according to Theorems 2 and 6, the stability of the model steady states does not depend on the value of the time delays, but rather on the basic reproduction number R0, only. In addition, we observe that the range of values for the two time delays does not lead to periodic solutions but an increase in both delays translate to an increase in the infectious population, both detected and undetected.Figure 1
Stability of the infected and free-infected equilibrium of model system (1) showing plots of I1(t) and I2(t) with varying delay (τ1=τ2). The direction of the arrow depicts an increase in delay with a step size of 2.0 starting from 2.0 to 10. The blue patterns in both (a) and (b) highlight brucellosis dynamics when R0<1 while the red pattern is for R0>1. Initial population levels were assumed as follows: S(0)=1000 animals, V(0)=500 animals, I1(0)=500 animals, and I2(0)=0 animals.
(a) (b)Figure2 depicts the numbers of infectious undetected and infectious detected animals over time with varying delays. The results clearly show that the incubation related delay (τ1) has more influence on shaping the dynamics of brucellosis compared to the culling related delay (τ2). More precisely, the incubation period delay significantly increases the infectious population (both detected and undetected) for 0<t<20 and thereafter its impact will be the same as that of detection (τ2).Figure 2
Numerical solutions of model system (1) illustrating the effects of different time delay on brucellosis infection level in the community. Initial population levels were assumed as follows: S(0)=1000 animals, V(0)=500 animals, I1(0)=500 animals, and I2(0)=0 animals.
(a) (b)In Figure3 we illustrate the stability of the disease-free equilibrium E0 with τ1=30 and τ2=5 (note that R0=0.686281). As we can observe, for certain parameter values and initial population levels, system (1) exhibits some periodic oscillation. Precisely, we note that the infected population (I1(t) and I2(t)) oscillates with a reduced amplitude from the start till when t is approximately 400; thereafter the oscillations dies off the solutions converges to the disease-free equilibrium. These simulation results demonstrate the occurrence of periodic solutions through Hopf bifurcation for delay values τ1=30 and τ2=5. In contrast, we can observe that there are no periodic oscillations for the uninfected populations S(t) and V(t).Figure 3
Numerical solutions demonstrating the stability ofE0 equilibrium of model system (1) with R0=0.686281. We set τ1=30, τ2=5, β=6.844×10-6animal-1year-1, γ=0.2, α=0.15year-1, A=16434 animals year-1, and the remainder retained the baseline values in Table 1. Further, we set the initial conditions as follows: S(0)=100 animals, V(0)=0 animals, I1(0)=10 animals, and I2(0)=0 animals.
(a) (b) (c) (d)In Figure4, we demonstrate the dynamic for model system (1) with respect to the stability of infection-free equilibrium for different pair of delay values (τ1,τ2) and from the simulation results we can conclude that both delays do not have a huge influence on the stability of disease-free equilibrium.Figure 4
Dynamics of model system (1) for different values of (τ1,τ2), which illustrate the stability of infection-free equilibrium E0 at R0=0.686281. We set β=6.844×10-6animal-1year-1, γ=0.2, α=0.15year-1, A=16434 animals year-1, and the remainder retained the baseline values in Table 1. Further, we set the initial conditions as follows: S(0)=100 animals, V(0)=0 animals, I1(0)=10 animals, and I2(0)=0 animals.
(a) (b) (c) (d)In Figure5 we observe that, for certain parameter values and initial population levels, system (1) may admit periodic oscillations when R0>1. As we can observe, when R0=3.77333 both the solutions of the infected and uninfected populations exhibit periodic oscillation for a certain period, before stability at endemic point is attained.Figure 5
Stability ofE∗ equilibrium of model system (1) with R0=3.77333. The time delay τ1 was fixed to be 60 and τ2 was fixed to be 1. We set the model parameters and variables as follows: β=6.844×10-6animal-1year-1, γ=0.2, α=0.015year-1, S(0)=100 animals, V(0)=0 animals, I1(0)=10 animals, and I2(0)=0 animals while the other parameter values are as in Table 1.
(a) (b) (c) (d)In Figure6, we illustrate the dynamics for model system (1) with respect to the stability of endemic equilibrium for several pair of delay values (τ1,τ2). The results confirm that the incubation related delay τ1 has more influence on shaping the dynamic of brucellosis compared to the culling related delay τ2.Figure 6
Numerical results of model system (1) for different values of (τ1,τ2), which demonstrate the stability of infected equilibrium E∗ at R0=3.77333. We set the model parameters and variables as follows: β=6.844×10-6animal-1year-1, γ=0.2, α=0.015year-1, S(0)=100 animals, V(0)=0 animals, I1(0)=10 animals, and I2(0)=0 animals while the other parameter values are as in Table 1.
(a) (b) (c) (d) (e) (f)To explore influence of model parameters on the reproduction numberR0, we perform a local sensitivity analysis of the basic reproduction number following the approach in [23]. The local sensitivity analysis will be useful on identifying parameters with greatest influence to change R0. To this end, denoting by Φ the generic parameter of system (1), we evaluate thenormalised sensitivity index:(70)SΦ=ΦR0∂R0∂Φ,which indicates how sensitiveR0 is to a change of parameter Φ. Model parameters with positive index increase the value of R0 whenever they are increased while those with a negative index decrease the value of R0 whenever they are increased. We consider the parameter values in Table 1, and we set α=0.015 in order to evaluate the normalized sensitivity index and the results are depicted in Figure 7. Here, we observe that parameters A, β, κ, γ have a positive correlation with R0, such that increasing these parameters will increase R0. However, it is the increase of A and β that has the greatest influence to change R0. Precisely, increasing either A or β by 50% will increase R0 by 50%. We also note that increasing parameters, c, μ, σ, d, p, and α, will lower the reproduction number.Figure 7
Sensitivity index forR0 with respect to model parameters that define it.


## 5. Conclusion

Zoonotic brucellosis remains a major public health problem in many developing nations. This is mainly attributed to several challenges associated with effective disease control in these nations. The challenges for effective control of brucellosis in developing nations range from inadequate veterinary personnel and vaccines as well as the failure by farmers to adhere to some of the aforementioned brucellosis control and eradication program activities. Furthermore, these challenges often lead to delay in detection and culling of infectious animals. In this study, we developed and analysed a mathematical model for brucellosis infection that incorporates two discrete delays. The first delay accounts for the latent period and the second delay represents the time taken to detect infectious animals. We computed the basic reproduction number and demonstrated that it is an important threshold quantity for stability of equilibria. By constructing suitable Lyapunov functionals, it has been shown that the model has a globally asymptotically stable infection-free equilibrium whenever the reproduction is less than unity. Further, it has been demonstrated that whenever the model reproduction number is greater than unity then the model has a unique endemic equilibrium point which is globally asymptotically stable. Numerical simulations are carried out to illustrate the main results.Although culling of symptomatic animals is a relatively easy strategy to implement, it is worth noting that some studies suggest that culling of both infected and susceptible animal may be more effective [24, 25]. The rationale is that by decreasing the host density, then the number of contacts per unit time between animals is low, thereby reducing disease transmission. In [25] it was demonstrated that culling of both susceptible and symptomatic animals only can be effective whenever the number of infected host is above a certain critical level [25]. We expect to improve this study in our future work by developing brucellosis model(s) with time delay that will enable us to compare aforementioned aspects.



---

*Source: 6456107-2018-08-07.xml*