URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu GLOBAL STABILITY OF A DELAYED MOSQUITO-TRANSMITTED DISEASE MODEL WITH STAGE STRUCTURE B

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu

GLOBAL STABILITY OF A DELAYED MOSQUITO-TRANSMITTED DISEASE MODEL

WITH STAGE STRUCTURE

B. G. SAMPATH ARUNA PRADEEP, WANBIAO MA

Abstract. This article presents a new eco-epidemiological deterministic de- lay differential equation model considering a biological controlling approach on mosquitoes, for endemic dengue disease with variable host (human) and vari- able vector (Aedes aegypti) populations, and stage structure for mosquitoes.

In this model, predator-prey interaction is considered by using larvae as prey and mosquito-fish as predator. We give a complete classification of equilibria of the model, and sufficient conditions for global stability/global attractivity of some equilibria are given by constructing suitable Lyapunov functionals and using Lyapunov-LaSalle invariance principle. Also, numerical simulations are presented to show the validity of our results.

1. Introduction

Recently, many scholars have proposed and investigated various kinds of epi- demic models in order to understand and describe the dynamics of infectious dis- ease. Most of the mathematical models pertinent to epidemiology are dependents of the baseline SIR (Susceptible, Infectious and Recovered) model which was pre- sented by Kermark and Mackendrick in 1927 based on ODE [23] with the concept of “compartment modelling”. By referring to the classical books [3, 5, 28, 32, 47], the readers can find not only the history of mathematical epidemiology but also the theories related delayed incorporated to biological systems. After Kermark and Mackendrick’s primus model numerous number of models emerged with time delay (see for example [1, 9, 10, 19, 18, 27, 31, 38, 42]), without time delay [8, 17, 46], epidemic model with stage structure [21, 45] and SVIR model with stochastic per- turbation [7].

Many serious epidemics such as AIDS and dengue (here, we mentioned few of them for more details in [4, 39]) can be transmitted horizontally as well as verti- cally. In addition, diseases also spread due to their parental genes [4]. The disease transmission via the vectors for examples West Nile fever, malaria, dengue and Rift vally fever has been studied by many researchers [11, 14, 15, 40]. With regard to dengue, one of the most spread out flavivirus disease propagated by adult female

2000Mathematics Subject Classification. 92D25, 34D23, 92B05, 93C23.

Key words and phrases. Epidemic model; time delay; Lyapunov functional;

global stability; nonlinear response.

c

2015 Texas State University - San Marcos.

Submitted April 22, 2014. Published January 7, 2015.

1

Aedes aegyptimosquito species (predominantly by this type, although there are oth- ers). Prevalent throughout the year in tropical countries but transmission reaches its peak when the highest rainfall takes place. It is known that this virus spreads among host (human) after an infectious matured femaleAedes aegypti(vector) hav- ing a blood meal from a susceptible host. On the other hand, the susceptible vector is infected after taking a blood meal from an infectious host. Although the disease spreads vertically among vectors, there is very low possibility of getting infected a new comer ([12] and references there in) which discourage us to include vertical transmission to our model. The most appropriate way to eradicate the transmis- sion of this viral disease is to control winged stage female mosquitoes and aquatic stage of mosquitoes because in the near future, it cannot be anticipated a vaccine to prevent from dengue fever ([36] and references there in). However, information regarding the possibility of vaccine and review of the development of a vaccine is found in [30]. Contemporary, eradication and control methods of mosquitos are similar to those arranged over half a century back. In the academic article [13], author has exhibited one of the controlling strategies named Sterile Insect Tech- nique (SIT) for the control of Aedes aegypti mosquitoes. Further, RIDL (Release of Insects Carrying a Dominant Lethal) based on new genetic sexing system for male mosquitoes is introduced by which allow only to born of male mosquitoes by blocking of female production of Aedes aegypti [20]. As a cost effective methods most countries use high toxic chemicals such as Malathion and insecticides to con- trol mosquito population which are very dangerous for public health. Places with immovable and clear water are available; the femaleAedes aegypti mosquitoes use those places for oviposition. By continuous awareness programs, the public can be made aware to avoid building up (source reduction) such places. Yet, it is hard to control aquatic stage of mosquitoes in the places such as lakes and ponds. As a biological control method, we can introduce mosquito-fish (predator) into wa- ter body in which immature mosquitoes (prey) are usually inhabited. In [35] the prey-dependent consumption predator prey model has been considered and some valuable results have been obtained.

Subsequently in the in 1920s, Lotka and Volterra introduced baseline model for predator-prey interaction, since then wide variety of modified and developed predator-prey models were seen in the literature. Further, in article [16], the authors have considered predator-prey model with infectious disease, models like SI, SIS and SIR with mass action incidence have been applied. In addition to that maturation time taken for the conversion process or gestation time delay for the predator or the prey in predator-prey model has been used to upgrade the model coherence with the natural world [26, 29, 37, 41, 44, 48]. It is ecologically important to look on the effects on stage-structure of immature to become mature of a certain species as all beings in the real world encounter to the stage structure. In the monograph [2], Aiello and Freedman suggested and analyzed stage structure model with constant maturation time delay for single species, for more examples one can refer [6, 24, 33, 34]. The authors in [43] considered a model with maturation delay and completely studied the stability properties and bifurcation analysis.

Motivated by the above works and as predator involvement is positively effect on controlling spread of vector transmission disease. We concern to construct a new echo-epidemic model for vector spread disease with both predator-prey interaction

and stage structure. Nonlinear functional response for predator-prey system is applied here. Biological appearance of our model is shown in Figure 1.

Figure 1. Mosquito life-cycle with predator and disease transmis- sion among human

We form a model which is described by a system of differential equations, with the aid of schematic diagram shown in Figure 2.

Let Sh(t), Ih(t) and Rh(t) represent the classes for the human population, in- dicating the number of susceptible, infective and recovered individuals at time t, respectively. Λ_h, γ_h, µ_h, q_v, τ, µ_l,µ_v, β_v, β_p, k_v, γ_v, µ_p, λ_landλ_pare all positive con- stants. Here,µ_h, µ_v, µ_l, andµ_p denote the per capita mortality rate of the human, vector, larvae and predator populations respectively. Λ_h indicates the recruitment of human to susceptible class while γ_h represents the recovered rate. Further, the vector population sub-divided into two classes, namely,S_v(t) is the number of sus- ceptible andIv(t) is the number of infective. It is assumed that vectors and human populations are mixing homogeneously. It is also accepted that the infected vectors will never be recovered and they transmit the virus in their entire life-span. The birth to the immature population (larvae) is proportional to the currently available number of matured population (vector) andqv is the conversion rate. In addition, λlis the rate of encounter andλp is the conversion efficiency. Λv represents the re- cruitment of immigrated mosquitoes to the susceptible class (see for example [22]).

We assumed that larvae are the only available food for predators. The density dependent mortality rate is denoted byβv (see, for example [2]) and the crowding rate is denoted β_p. The state variables L_v(t) and P(t) show up the number of larvae and number of predator at time t, respectively. The larvae who was born

Figure 2. Transfer diagram of the disease in human and vectors with predator

at t−τ and still survive at timet, transforming from larvae (immature stage) to susceptible vector (matured stage) is given by the term bqve^{−τ µl}Sv(t−τ) and to infected vector (matured stage) by aqve^{−τ µl}Iv(t−τ), where a and b are positive constants such thata+b= 1. The number of recovered human at timetis denoted byRh(t) which has not appeared in the other equations of system (1.1). Therefore, it is excluded from further consideration. The formulated model is given below.

S˙h(t) = Λh−βvhIv(t)Sh(t)−µhSh(t), I˙h(t) =khβvhIv(t)Sh(t)−(µh+γh)Ih(t),

S˙v(t) = Λv+bqve^{−τ µl}Sv(t−τ)−βhvIh(t)Sv(t)−µvSv(t)−βvS_v²(t), I˙v(t) =aqve^{−τ µl}Iv(t−τ) +kvβhvIh(t)Sv(t)−(µv+γv)Iv(t),

L˙_v(t) =b[S_v(t)−e^{−τ µl}S_v(t−τ)] +a[I_v(t)−e^{−τ µl}I_v(t−τ)]

−µlLv(t)−λlf(Lv(t))P(t),

P(t) =˙ P(t)(−µp−βpP(t)) +λpf(Lv(t))P(t).

(1.1)

In system (1.1) it is assumed that at timetthe number ofβvhIv(t)Sh(t) is removed from susceptible human class and simultaneously the number of khβvhIv(t)Sh(t) enters to infected human class. It is further assumed that at timet the number of βhvIh(t)Sv(t) is removed from susceptible mosquito class and simultaneously the number ofkvβhvIh(t)Sv(t) enters to infected mosquito class, whereβvh,βhv,khand k_v are positive constants. f(L_v) is a function such that monotonically increasing, positive and differentiable for all L_v > 0 andf(0) = 0. More general model for

system (1.1) is formulated as follows:

S˙h(t) = Λh−βvhIv(t)Sh(t)−µhSh(t),

I˙h(t) =khβvhe^−µhσIv(t−σ)Sh(t−σ)−(µh+γh)Ih(t), S˙v(t) = Λv+bqve^{−τ µl}Sv(t−τ)−βhvIh(t)Sv(t)−µvSv(t)−βvS_v²(t), I˙_v(t) =aq_ve^{−τ µl}I_v(t−τ) +k_vβ_hve^−µvωI_h(t−ω)S_v(t−ω)−(µ_v+γ_v)I_v(t),

L˙_v(t) =b[S_v(t)−e^{−τ µl}S_v(t−τ)] +a[I_v(t)−e^{−τ µl}I_v(t−τ)]

−µlLv(t)−λlf(Lv(t))P(t),

P˙(t) =P(t)(−µp−β_pP(t)) +λ_pf(L_v(t−ρ))P(t−ρ),

(1.2) where σ >0 is the latent delay for human,ω >0 is the latent delay for mosquito and ρ >0 is the predation delay and others have the same biological meaning as in system (1.1).

This work is organized as follows. In next section, we state some lemmas that are important for our discussion and show the existence, boundedness of the solutions of system (1.1) with the initial condition (2.1). Moreover, conditions are given for existence of all kinds of equilibria of system (1.1); the reproduction number is simultaneously calculated. In Section 3, stability properties of some equilibria are established by means of Lyapunov functionals, and instability of equilibria is given by using characteristic equations. Further, the global attractiveness of other equilibria is also considered. Numerical simulations of the results are presented in Section 4. The paper ends with a brief discussion.

2. Boundedness of solutions and analysis of equilibria

In this section, we consider the boundedness of solutions and existence of equi- libria for system (1.1). First, the initial conditions of system (1.1) are

S_h(θ) =ϕ₁(θ), I_h(θ) =ϕ₂(θ), S_v(θ) =ϕ₃(θ),

Iv(θ) =ϕ4(θ), Lv(θ) =ϕ5(θ), P(θ) =ϕ6(θ), (2.1) where (−τ ≤ θ ≤ 0) and ϕi(θ) (i = 1,2, . . . ,6) belong to Banach space C = C([−τ,0], R+) of continuous functions mapping from the interval [−τ,0] intoR+:=

[0,+∞), equipped with the supremum norm.

Using the basic theory of delay differential equations (see, for example, [24]), it is not difficult to show that, for the initial conditions given above, the solution (S_h(t),I_h(t),S_v(t),I_v(t),L_v(t),P(t)) of system (1.1) exists and unique for all time t≥0. The following lemma is used to obtain our results.

Lemma 2.1. Consider the delay differential equationx(t) =˙ α+βx(t−τ)−γx(t)−

δx²(t)whereα, β, γ >0,δ≥0,τ≥0 are constants, the initial functionϕ∈Cand ϕ(θ)>0.

(i) If δ >0 andγ≥β then unique positive equilibrium x^∗=β−γ+p

(β−γ)²+ 4αδ

2δ ,

is globally asymptotically stable.

(ii) If δ = 0 and γ > β, then unique positive equilibrium x^∗ = α/(γ−β) is globally asymptotically stable.

Proof. Let us consider case (i). It is not difficult to verify that the solutionx(t) with any initial functionϕof this equation is positive. Fort≥0, let us define

V =x−x^∗−x^∗ln x x^∗ +β

Z t

t−τ

x(θ)−x^∗−x^∗lnx(θ) x^∗

dθ.

By considering the time derivative along the solution we have that V˙ = (γ−β)x^∗

2− x x^∗ −x^∗

x

+βx^∗

1−x(t−τ)

x + lnx(t−τ) x

− δ

x(x−x^∗)²(x+x^∗).

If γ ≥ β, it is clear that ˙V ≤ 0. Therefore, it follows easily from Corollary 5.2 of Kuang [24] that x^∗ is globally asymptotically stable. Proof of part (ii) can be

obtained by using similar method as in part (i).

For biological reasons, throughout this paper we discuss the dynamical behavior of system withµv≥max[aqve^{−τ µl}, bqve^{−τ µl}].

For boundedness of the solutions to (1.1), we have the following result.

Theorem 2.2. Every solution (Sh(t), Ih(t), Sv(t), Iv(t), Lv(t), P(t)) of (1.1) has the following properties:

lim sup

t→+∞

S_h(t) + 1 k_hI_h(t)

≤M^Sh, lim sup

t→+∞

I_h(t)≤M^Ih, lim sup

t→+∞

S_v(t)≤M^Sv lim sup

t→+∞

I_v(t)≤M^Iv, lim sup

t→+∞

Lv(t)≤M^Lv, lim sup

t→+∞

P(t)≤M^P, where

M^Sv =S_v^e0, M^Ih =k_hβ_vhS^e_h⁰ µh+γh

M^Iv, M^Iv = khkvβhvS_h^e0

µv+γv−aqve^{−τ µl}M^Sv, M^Lv =bM^Sv +aM^Iv µl

, M^P =

(_λ

pf(M^Lv)−µp

β_p λpf(M^Lv)> µp, 0, λpf(M^Lv)≤µp,

M^Sh =S_h^e0 =Λh

µh

, S_v^e0= 1

2βv

q_vbe^{−τ µl}−µ_v+p

(q_vbe^{−τ µl}−µ_v)²+ 4β_vΛ_v . Proof. From the first two equations of (1.1), fort≥0, we have

Sh(t) + 1 k_hIh(t)0

≤Λh−µh

Sh(t) + 1 k_hIh(t)

, from which, we have

lim sup

t→+∞

(Sh(t) +Ih(t)/kh)≤S_h^e0, lim sup

t→+∞

Ih(t)≤khS_h^e0. Hence, from (1.1) it follows that fort≥0,

S˙_v(t)≤Λ_v+bq_ve^{−τ µl}S_v(t−τ)−µ_vS_v(t)−β_vS_v²(t).

The well-known comparison theorems for delay differential equations (see, for ex- ample [25]) and Lemma 2.1 imply that lim sup_t→+∞ Sv(t)≤M^Sv.

For any sufficiently smallε >0, there exists some sufficiently largeT, such that, fort≥T, it has from system (1.1) that

I˙v(t)≤aqve^{−τ µl}Iv(t−τ) +kvβhv(khS_h^e0+ε)(M^Sv +ε)−(µv+γv)Iv(t).

Hence, from Lemma 2.1, it is easy to obtain that lim sup

t→+∞

Iv(t)≤k_vβ_hv(k_hS_h^e0+ε)(M^Sv+ε) µv+γv−aqve^{−τ µl} . Letε→0⁺, it has that lim sup_t→+∞Iv(t)≤M^Iv.

With similar arguments as above, for any sufficiently small ε >0, there exists some sufficiently largeT1, such that fort≥T1, from system (1.1), we have

I˙h(t)≤khβvh(M^Iv+ε)(S_h^e0+ε)−(µh+γh)Ih(t), L˙v(t)≤b(M^Sv+ε) +a(M^Iv+ε)−µlLv(t), from which it is easy obtain that

lim sup

t→+∞

I_h(t)≤k_hβ_vh(M^Iv+ε)(S_h^e0+ε)/(µ_h+γ_h), lim sup

t→+∞

L_v(t)≤(b(M^Sv+ε) +a(M^Iv +ε))/µ_l. Further, by lettingε→0⁺, one has that

lim sup

t→+∞

Ih(t)≤M^Ih, lim sup

t→+∞

Lv(t)≤M^Lv.

For a sufficiently smallε >0, there exists some sufficiently large T2, such that fort≥T₂, from system (1.1), one has

P(t)˙ ≤P(t)(−µp−βpP(t)) +λpf(M^Lv +ε)P(t)

= [λ_pf(M^Lv +ε)−µ_p−β_pP(t)]P(t).

Ifλpf(M^Lv)> µp, for sufficient smallε >0, it has thatλpf(M^Lv+ε)−µp>0.

Hence, we can obtain that lim sup_t→+∞P(t)≤ λ_pf(M^Lv+ε)−µ_p

/β_p. By letting ε→0⁺, one has that lim sup_t→+∞P(t)≤ λpf(M^Lv)−µp

/βp.

Further, if λpf(M^Lv) ≤ µp, it has that for t ≥ T2, ˙P(t) ≤ (ε−βpP(t))P(t), which implies that lim sup_t→+∞P(t)≤ε/βp. By letting ε → 0⁺, we have that lim sup_t→+∞P(t) =0. Therefore, it proves that lim sup_t→+∞P(t)≤M^P. Next, we study the existence of all possible nonnegative equilibria of system (1.1). We have following four cases to be considered.

(i) There exists boundary equilibrium (disease-free and predator-free ) E0 = (S_h^e0,0, S_v^e0,0, L^e_v⁰,0), where

L^e_v⁰= b µl

(1−e^{−τ µl})S_v^e0,

(ii) Iff(L^e_v¹)> µp/λpholds, there exists boundary equilibrium (disease-free with predator)E1 = (S_h^e1,0, S_v^e1,0, L^e_v¹, P^e1), where S_h^e1 =S_h^e0, S_v^e1 =S^e_v⁰. From last two equations of system (1.1), we have thatP^e1= λ_pf(L^e_v¹)−µ_p

/β_p and G(L_v) =b(1−e^{−τ µl})S_v^e1−µ_lL_v−λ_lf(L_v)P^e1 = 0,

whereL_v is any value which satisfy the fifth equation.

It is easy to show thatG(0) =b(1−e^{−τ µl})S_v^e1 >0 andG(L1) =−λlf(L1)P^e1<

0. Therefore, there exists a uniqueL^e_v¹∈(0, L1) whereL1=b(1−e^{−τ µl})S_v^e1/µl. (iii) IfS_v^e0< S_v^e2R0holds, there exists boundary equilibrium (predator free with disease)E₂= (S_h^e2, I_h^e2, S_v^e2, I_v^e2, L^e_v²,0), where

S_h^e2 =S_v^e0S_h^e0 S^ev²R0

, I_h^e2= k_hΛ_h µh+γh

1− S_v^e0 R0Sv^e2

, I_v^e2 = kvβhvS_v^e2I_h^e2

(µ_v+γ_v−aq_ve^{−τ µl}), L^e_v² =(bS_v^e2+aI_v^e2)

µ_l (1−e^{−τ µl}).

S_v^e2= −K+p

K²+ 4βvM

/2βv is given byβv(S_v^e2)²+KS_v^e2−M = 0, where K= β_hvk_hΛ_h

µh+γh

+µ_v−bq_ve^{−τ µl}, M = Λ_v+(µ_v+γ_v−aq_ve^{−τ µl})µ_h kvβvh

, and

R0= khβvh

µv+γv−aqve^{−τ µl}S_h^e0 kvβhv

µh+γh

S_v^e0 is the basic reproduction number of system (1.1).

(iv) IfS_v^e0 < S^∗_vR0andf(L^∗_v)> µp/λp hold, there exists a unique positive equi- librium (disease with predator)E^∗= (S_h^∗, I_h^∗, S_v^∗, I_v^∗, L^∗_v,P^∗), whereS_h^∗=S_h^e2, I_h^∗= I_h^e2, S_v^∗=S^e_v² andI_v^∗=I_v^e2.

From last two equations of system (1.1), we have that P^∗= (λpf(L^∗_v)−µp)/βp

and

H(Lv) = bS_v^∗+aI_v^∗

(1−e^{−τ µl})−µlLv−λlf(Lv)P^∗= 0,

whereLvis any value which satisfy the fifth equation. It is not difficult to show that H(0) = (bS_v^∗+aI_v^∗)(1−e^{−τ µl})S_v^∗ >0 andH(L2) = −λlf(L2)P^∗ <0. Therefore, there exist a uniqueL^∗_v∈(0, L2) whereL2= (bS_v^∗+aI_v^∗)(1−e^{−τ µl})/µl.

3. Stability of equilibria

In this section, we analyze the stability properties of each equilibrium of sys- tem (1.1). The characteristic equation of system (1.1) at any equilibrium E = (Sh, Ih, Sv, Iv, Lv, P) has the form

F(λ, τ) (3.1)

=

λ+a1 0 0 βvhSh 0 0

−khβvhIv λ+b1 0 −khβvhSh 0 0

0 βhvSv λ+q(λ, τ) 0 0 0

0 −kvβhvSv −kvβhvIh λ+c(λ, τ) 0 0

0 0 g(λ, τ) j(λ, τ) λ+d λlf(Lv)

0 0 0 0 −λpf⁰(Lv)P λ+e

= 0,

where

a1=βvhIv+µh, b1=µh+γh, c(λ, τ) =µv+γv−aqve^{−(λ+µl)τ}, d=µl+λlf⁰(Lv)P, e=µp+ 2βpP−λpf(Lv),

g(λ, τ) =b(e^{−(λ+µl)τ}−1), q(λ, τ) =h−qvbe^{−(λ+µl)τ}, h=βhvIh+µv+ 2βvSv, j(λ, τ) =a(e^{−(λ+µl)τ}−1).

In next theorem, we establish stability properties of the equilibriumE₀.

Theorem 3.1. The following conclusions hold for any time delayτ ≥0.

(a) IfR0≤1 andf(L^e_v⁰)≤µp/λp thenE1, E2 andE^∗ are not existent, andE0

is globally asymptotically stable,

(b) EitherR₀>1 orf(L^e_v¹)> µ_p/λ_p holds, thenE₀ unstable.

Proof. The characteristic equation (3.1) atE0 is reduced to

(λ+a1)(λ+q(λ, τ))(λ+d)(λ+e)∆1(λ, τ) = 0, (3.2) where

∆₁(λ, τ) = (λ+b₁)(λ+c(λ, τ))−k_vβ_hvS_v^e0k_hβ_vhS_h^e0

=λ²+ (µh+γh+µv+γv−aqve^{−(λ+µl)τ})λ + (µh+γh)(µv+γv−aqve^{−(λ+µl)τ})(1−R0).

By following a similar procedure as the one used in [24] it easy to show for the characteristics equation (3.2) that if f(L^e_v⁰) < µp/λp and R0 ≤ 1 hold, E0 of system (1.1) is locally asymptotically stable.

Define a Lyapunov functional as W1= k1

β_vhS_h^e0V1+ k1

k_hβ_vhS_h^e0Ih+ 1 β_hvSv^e0

V2+ 1 k_vβ_hvSv^e0

V3+U, wherek1 is determined later, and

V1=Sh−S^e_h⁰−S_h^e0ln Sh

S^e_h⁰, V2=Sv−S_v^e0−S^e_v⁰ln Sv

Sv^e0

, V3=Iv+aqve^{−τ µl}

Z t

t−τ

Ivdt, U = qvbe^{−τ µl} βhvSv^e0

Z t

t−τ

[Sv−S_v^e0−S_v^e0ln Sv

Sv^e0

]dt.

By considering the derivative along the solution, we have W˙1= k1

β_vhS_h^e0(1−S_h^e0

S_h)(Λh−βvhIvSh−µhSh) + k₁

khβvhS^e_h⁰(k_hβ_vhI_vS_h−(µ_h+γ_h)I_h)

+ 1

βhvSv^e0

(1−S_v^e0 Sv

)(Λ_v+bq_ve^{−τ µl}S_v(t−τ)−β_hvI_hS_v−µ_vS_v−β_vS_v²)

+ 1

kvβhvSv^e0

(aq_ve^{−τ µl}I_v(t−τ) +k_vβ_hvI_hS_v−(µ_v+γ_v)I_v) + aq_ve^{−τ µl}

kvβhvSv^e0

(I_v(t)−I_v(t−τ)) +q_vbe^{−τ µl}

βhvSv^e0

(Sv−Sv(t−τ) +S_v^e0lnS_v(t−τ) Sv

).

By choosing k1 = ^µv−aq_k ^{ve−τ µl+γv}

vβ_hvS_v^e0 and noting that Λv = (µv −bqve^{−τ µl})S_v^e0 + βv(S^e_v⁰)², Λh=µhS_h^e0, we have that

W˙1= k1µh

β_vh (2− Sh

S_h^e0 −S^e_h⁰

S_h) + (1− 1

R₀)Ih+(µv−qvbe^{−τ µl})

β_hv (2− Sv

Sv^e0

−S_v^e0 S_v )

− βv

βhvSvSv^e0

(S_v−S^e_v⁰)²(S_v+S_v^e0)

+qvbe^{−τ µl} βhv

1−Sv(t−τ) Sv

+ lnSv(t−τ) Sv

.

It can be shown that ifR₀ ≤1, then ˙W₁≤0. Define the subsetE={ϕ= (ϕ₁, ϕ2, ϕ3,ϕ4,ϕ5, ϕ6)|W˙1(ϕ) = 0}. Further, letM be the largest invariant set inE with respect to system (1.1). Let us further show thatM ={E0}. Denote

Sht(θ) =Sh(t+θ), Svt(θ) =Sv(t+θ), I_ht(θ) =I_h(t+θ), I_vt(θ) =I_v(t+θ),

Lvt(θ) =Lv(t+θ), Pt(θ) =P(t+θ) (−τ≤θ≤0).

For any solution (Sht, Iht, Svt, Ivt, Lvt, Pt) with the initial functionϕ = (ϕ1, ϕ2,ϕ3,ϕ4,ϕ5,ϕ6)∈M, one has, from the invariance, that for allt∈R, (Sht,Iht, Svt,Ivt,Lvt,Pt)∈M.

IfR0<1, then ˙W1= 0 if and only if Iht(0) = 0 and S_ht(0)

S_h^e0 =S_vt(0) Sv^e0

=S_ht(−τ) Sv

= 1.

Hence, from invariance of M, we can obtain that Iht(0) = Ih(t) = 0, Sht(0) = Sh(t) = S^e_h⁰, Svt(0) =Sv(t) = S_v^e0. Further from first equation of (1.1), we can obtain thatIvt(0) =Iv(t) = 0.

IfR0= 1, then ˙W1= 0 if and only if S_ht(0)

S_h^e0 =S_vt(0) Sv^e0

= S_vt(−τ) Sv

= 1.

Hence, for all t ∈ R, one has that Sht(0) = Sh(t) = S_h^e0, Svt(0) = Sv(t) = S_v.^e0. From first and third equations of system (1.1), we can show thatIht(0) =Ih(t) = 0 andIvt(0) =Iv(t) = 0 respectively, for allt∈R. Therefore, for allt∈R, in subset M last two equations of system (1.1) are reduced to

L˙v(t) =b(1−e^{−τ µl})S_v^e0−µlLv(t)−λlf(Lv(t))P(t), P(t) =˙ P(t)(−µp−βpP(t)) +λpf(Lv(t))P(t).

Define another Lyapunov functional as W2=Lv−L^e_v⁰−

Z L_v

L^ev⁰

f(L^e_v⁰)

f(θ) dθ+k2P,

wherek2due to be determined later. By taking time derivative along the solution, we have that

W˙2=

1−f(L^e_v⁰) f(Lv)

b(1−e^{−τ µl})S_v^e0−µlLv−λlf(Lv)P +k₂ (−µp−β_pP)P+λ_pf(L_v)P

.

Lettingk2=λl/λp and noting thatb(1−e^{−τ µl})S^e_v⁰ =µlL^e_v⁰, we have that W˙2=µl(L^e_v⁰−Lv)

1−f(L^e_v⁰) f(L_v)

+λlP

f(L^e_v⁰)−µp

λ_p −λl

λ_pβpP².

Iff(L^e_v⁰)≤µp/λpholds, and from properties of the function we have that ˙W2≤0.

Further, W˙₂ = 0 if and only if L_v(t) = L^e_v⁰ and P(t) = 0. This proves that

M ={E0}. By Lyapunov-LaSalle invariance principle from [24], one proves that E0 is globally attractive. Hence,E0 is globally asymptotically stable.

IfR0>1, we have from (3.2) that limλ→+∞∆1(λ, τ) = +∞ and ∆1(0, τ)<0.

Hence, ∆₁(λ, τ) = 0 has at least one positive real root.

Next, we consider the factor λ+e=λ+µ_p−λ_pf(L^e_v⁰) = 0. Clearly, it has a positive real root whenf(L^e_v⁰)> µ_p/λ_p. Remark: At equilibriumE1predators consume some larvae (i.e. f(L^e_v¹)≤f(L^e_v⁰).

Hence, if (f(L^e_v¹)≤)f(L^e_v⁰)≤µp/λp holds, E0 stable. From which, it implies that E1 is not existent. Further, it is clear from Theorem 2.2 and Theorem 3.1 that S_v^e2< S^e_v⁰ and ifR0≤1 holdsE0 stable, respectively. From which, it implies that E2andE^∗are not existent. Therefore, ifE0stable,E1,E2andE^∗are not existent.

In Theorem 3.2, we establish the global stability properties of the equilibrium E₁.

Theorem 3.2. If f(L^e_v¹)> µ_p/λ_p(i.e. E₀is unstable), then following conclusions hold for any time delayτ ≥0.

(a) IfR₀≤1(i.e. E₂andE^∗are not existent), thenE₁is globally asymptotically stable,

(b) IfR0>1 holds, thenE1 is unstable.

Proof. At equilibriumE1the characteristic equation (3.1) becomes

(λ+a1)(λ+q(λ, τ))[(λ+d)(λ+e) +λlλpP^e1f⁰(L^e_v¹)f(L^e_v¹)]∆2(λ, τ) = 0, (3.3) where

∆₂(λ, τ) = (λ+b₁)(λ+c(λ, τ))−k_vβ_hvS_v^e1k_hβ_vhS_h^e1

=λ²+ (µ_h+γ_h+µ_v+γ_v−aq_ve^{−(λ+µl)τ})λ + (µ_h+γ_h)(µ_v+γ_v−aq_ve^{−(λ+µl)τ})(1−R₀).

By following a similar procedure as in [24] it is easy to show, for the characteristics equation (3.3), that if R₀ ≤1 holds, E₁ of system (1.1) is locally asymptotically stable.

We define the same Lyapunov functional (W₁) and applying the same procedure as in proof of Theorem 3.1, we can show thatS_h(t) =S_h^e1,I_h(t) = 0,S_v(t) =S_v^e1, Iv(t) = 0 onM. For all t∈R, in subset M last two equations of system (1.1) are reduced to

L˙v(t) =b(1−e^{−τ µl})S_v^e1−µlLv(t)−λlf(Lv(t))P(t), P(t) =˙ P(t)(−µp−βpP(t)) +λpf(Lv(t))P(t).

Define a Lyapunov functional as W₃=L_v−L^e_v¹−

Z Lv

L^e_v¹

f(L^e_v¹) f(θ) dθ+ λ_l

λp

P−P^e1−P^e1ln P P^e1

. By taking the time derivative along the solution, we have that

W˙₃=

1−f(L^e_v¹) f(Lv)

(b(1−e^{−τ µl})S_v^e1−µ_lL_v−λ_lf(L_v)P) + λl

λ_p 1−P^e1 P

(−µp−βpP)P+λpf(Lv)P .

Noting that b(1−e^{−τ µl})S_v^e1 =µlL^e_v¹ +λlf(L^e_v¹)P^e1 and µp =λpf(L^e_v¹)−βpP^e1, we have that

W˙3=µl(L^e_v¹−Lv)

1−f(L^e_v¹) f(Lv)

+λlf(L^e_v¹)

2−f(L^e_v¹)

f(Lv)−f(Lv) f(L^ev¹)

−λlβp

λp

(P−P^e1)². From the properties of the functionf we have that ˙W3≤0. Further, ˙W3= 0 if and only if Lv(t) = L^e_v¹ and P(t) =P^e1. This proves that M ={E1}. By Lyapunov- LaSalle invariance principle from [24],E1 is globally attractive. It proves that E1

is globally asymptotically stable.

If R0 >1 holds, we can easily show from (3.3) that limλ→+∞∆2(λ, τ) = +∞

and ∆2(λ, τ)<0. Hence, ∆2(λ, τ) = 0 has at least one positive real root.

Theorem 3.3. If S_v^e0 < S^e_v²R0 (i.e. E0,E1 are unstable), then following conclu- sions hold for any time delayτ ≥0.

(a) Iff(L^e_v²)≤µ_p/λ_p (i.e. E^∗ is not existent), thenE₂ is globally attractive, (b) Iff(L^∗_v)> µ_p/λ_p holds, thenE₂ is unstable.

Proof. Define a Lyapunov functional W4= 1

βvhIv^e2S_h^e2V1(t) + 1

khβvhIv^e2S_h^e2V2(t) + 1 βhvI_h^e2Sv^e2

(V3(t) +U1)

+ 1

k_vβ_hvI_h^e2Sv^e2

(V₄(t) +U₂), where

V₁(t) =S_h−S_h^e2−S_h^e2ln S_h

S_h^e2, V₄(t) =I_v−I_v^e2−I_v^e2ln I_v Iv^e2

, V₂(t) =I_h−I_h^e2−I_h^e2ln Ih

I_h^e2, U₁=q_vbe^{−τ µl} Z t

t−τ

[S_v−S_v^e2−S_v^e2ln Sv

Sv^e2

]dt, V3(t) =Sv−S_v^e2−S_v^e2ln Sv

Sv^e2

, U2=aqve^{−τ µl} Z t

t−τ

[Iv−I_v^e2−I_v^e2ln Iv

Iv^e2

]dt.

By considering time derivative along the solution, we have W˙₄= 1

βvhIv^e2S_h^e2(1−S_h^e2 Sh

)(Λ_h−β_vhI_vS_h−µ_hS_h)

+ 1

khβvhIv^e2S^e_h²(1−I_h^e2 Ih

)(khβvhIvSh−(µh+γh)Ih)

+ 1

βhvI_h^e2Sv^e2

(1−S_v^e2 Sv

)(Λv+bqve^{−τ µl}Sv(t−τ)−βhvIhSv−µvSv−βvS_v²) + q_vbe^{−τ µl}

βhvI_h^e2Sv^e2

(S_v−S_v(t−τ) +S^e_v²lnS_v(t−τ) Sv

)

+ 1

kvβhvI_h^e2Sv^e2

(1−I_v^e2 Iv

)(aqve^{−τ µl}Iv(t−τ) +kvβhvIhSv−(µv+γv)Iv) + aq_ve^{−τ µl}

kvβhvI_h^e2Sv^e2

(I_v−I_v(t−τ) +I_v^e2lnI_v(t−τ) Iv

).

Note that

Λh=βvhI_v^e2S_h^e2+µhS_h^e2, Λv= (µv−bqve^{−τ µl})S^e_v²+βhvI_h^e2S^e_v²+βv(S^e_v²)², k_hβ_vhI_v^e2S_h^e2 = (µ_h+γ_h)I_h^e2, k_vβ_hvI_h^e2S^e_v² = (µ_v−aq_ve^{−τ µl}+γ_v)I_v^e2.