Physica A xxx (xxxx) xxx
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
A periodic malaria model with two delays✩ ∗
Yan Zhang a,b , , Sanyang Liu a , Zhenguo Bai a a b
School of Mathematics and Statistics, Xidian University, Xi’an, 710071, China Xi’an University of Architecture and Technology Huaqing College, Xi’an, 710043, China
article
info
Article history: Received 15 November 2018 Received in revised form 13 April 2019 Available online xxxx MSC: 92D30 37N25 34K13 Keywords: Malaria model Basic reproduction number Time delay Threshold dynamics
a b s t r a c t Malaria is the world’s most prevalent mosquito-borne disease caused by Plasmodium parasites, and responsible for over half a million deaths per year. To understand the effects of the intrinsic and extrinsic incubation periods of the parasite within the humans and mosquitoes, respectively, and the seasonality on disease transmission, we propose a periodic malaria model with delays. The basic reproduction number R0 is derived, and it is shown that R0 is a threshold parameter between the extinction and persistence of the disease. In the case where all the coefficients are constants and the intrinsic incubation period is ignored, we also prove the global attractivity of the endemic equilibrium when R0 > 1. Numerical simulations indicate that prolonging the incubation period in mosquitoes is more effective than prolonging the incubation period in humans for disease control. It is also found that increasing the strength of seasonal forcing will lead to a higher epidemic peak. © 2019 Published by Elsevier B.V.
1. Introduction Malaria is a mosquito-borne disease caused by the Plasmodium parasite and transmitted between humans through bites of female Anopheles mosquitoes [1]. It is a persistent and recurring infectious disease that affects over 100 countries in Africa, Southeast Asia, and the Americas [2]. The disease directly threatens about half of the world’s population, leading to some 600,000 deaths per year [3]. Although mortality rates associated with malaria infection have been reduced, continued and increased efforts to find ways to effectively control the disease are required to achieve elimination and eradication objectives [4]. Mathematical models have been designed and used to understand the transmission dynamics of an infectious disease and thereby, provide guides and suggestions for disease control [5]. The earliest model of malaria transmission is the Ross-MacDonald model [6], which provides insights into the mechanism of malaria transmission and spread through a system of ordinary differential equations. Since then great progress has been made to extend this classical model by considering various aspects related to epidemiological features of malaria, such as [5,7–10]. However, most of these models are investigated under the assumption that the total population is constant. It should be noted that this assumption is relatively valid only when studying diseases of short duration with limited effects on mortality but may no longer hold when dealing with endemic diseases such as malaria, which has a high mortality [11]. It is therefore necessary to develop the model with varying population size to study the malaria transmission dynamics. Apart from considering ✩ This research was supported by the NSF of China (Nos. 11971369, 11671315), the NSF of Shaanxi Province of China (Nos. 2017JM1001, 2019JM-241) and the Fundamental Research Funds for the Central Universities (No. JB190705). ∗ Corresponding author. E-mail address:
[email protected] (Y. Zhang). https://doi.org/10.1016/j.physa.2019.123327 0378-4371/© 2019 Published by Elsevier B.V.
Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
2
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
variable human and mosquito populations, the model studied in this paper also incorporates the following two important biological aspects. The first aspect is the extrinsic and intrinsic incubation periods of parasite in the feedback interactions between mosquitoes and humans. The extrinsic incubation period (EIP) is referred to as the number of days between a mosquito’s infection and when it can yield infectious bites [12]. This period is in a range of 8–14 days depending on malaria parasite species and temperature [13]. The intrinsic incubation period (IIP) is the length of time between the mosquito bite introducing the parasite into the blood stream and the appearance of symptoms, varying from 7 to 30 days [14]. An appropriate candidate to describe the EIP or IIP in mathematics is time delay. There have been some works describing the latencies of the developments of parasites within humans and mosquitoes, see, e.g., [9,10,15] and references therein. The second aspect is the seasonality in malaria transmission. Climate factors, such as temperature, humidity, rainfall, and wind, are known to significantly affect the distribution and seasonal dynamics of mosquito populations, with which substantial implications for disease seasonality and persistence [16,17]. Hence, understanding the relationship between mosquitoes abundance and environmental drivers, such as temperature, is fundamental to the study of the epidemiology of malaria. To achieve this objective, the current study will consider climate factors. The paper is organized as follows. The model is formulated in Section 2. The basic reproduction number R0 and its equivalent characterization is given in Section 3. In Section 4, we establish a threshold result on the disease dynamics in terms of R0 , and for the case of constant coefficients and ignoring the IIP, we also prove the global attractivity of the endemic equilibrium when R0 > 1. Numerical simulations are performed in Section 5 to show the effects of seasonality and delays on disease transmission. Discussions are provided in the final section. 2. The model The model to be designed is motivated by the malaria transmission models in [15,18]. The human population is divided into three epidemiological classes such that at time t, there are Sh (t) susceptible, Ih (t) infectious, and Rh (t) immune humans. The mosquito population is described by two groups such that at time t, there are Sm (t) susceptible, and Im (t) infectious adult female mosquitoes. Denote Nh (t) = Sh (t) + Ih (t) + Rh (t) and Nm (t) = Sm (t) + Im (t) are the total human and mosquito populations at time t, respectively. Note that the subscripts h and m stand for humans and mosquitoes, respectively. Since the climate factor has little impact on human activities, we treat all of the parameters related to humans as constants. Suppose that Λh is the recruitment rate for humans, dh is the natural death rate of humans and αh is the recovery rate of humans. To describe the effects of climate factors on mosquito development, we let Λm (t) be the recruitment rate at which female adult mosquitoes emerge from larvae, and let dm (t) be the mortality rate for female adult mosquitoes. The biting rate β (t) of mosquitoes is the average number of bites per mosquito per unit time at time t. We assume that the total number of bites made by mosquitoes is equal to the total number of bites received by humans, and a mosquito will not bite the same person more than once. Then the number of newly infected humans and newly infected mosquitoes per unit time at time t is, respectively, given by c β (t)
Sh (t) Nh (t)
Im (t)
and bβ (t)
Ih (t) Nh (t)
Sm (t),
where c(b) is the transmission probability from an infectious mosquito (human) to a susceptible human (mosquito) given that contact between the two occurs. The incubation period in a human has duration τ1 . As mentioned in [19], it is possible that some individuals recovered from parasitemia during the IIP, and thus, the probability that a human survives the IIP at time t is e−(dh +αh )τ1 . Then the number of newly infectious humans per unit time at time t is c β (t − τ1 )
Sh (t − τ1 ) Nh (t − τ1 )
Im (t − τ1 )e−(dh +αh )τ1 .
Similarly, let τ2 be the length of the EIP, we can show that bβ (t − τ2 )
Ih (t − τ2 ) Nh (t − τ2 )
∫ − tt−τ dm (s)ds 2
Sm (t − τ2 )e
is the incidence of newly infectious mosquitoes at time t. Consequently, our model under consideration with delays and seasonality can be formulated as follows:
⎧ dSh (t) Sh (t) ⎪ ⎪ = Λh − c β (t) Im (t) − dh Sh (t), ⎪ ⎪ dt N h (t) ⎪ ⎪ ⎪ ⎪ dIh (t) Sh (t − τ1 ) ⎪ ⎪ = c β (t − τ1 ) Im (t − τ1 )e−(dh +αh )τ1 − (dh + αh )Ih (t), ⎪ ⎪ dt Nh (t − τ1 ) ⎪ ⎪ ⎨ dRh (t) = αh Ih (t) − dh Rh (t), ⎪ dt ⎪ ⎪ ⎪ Ih (t) dSm (t) ⎪ ⎪ = Λm (t) − bβ (t) Sm (t) − dm (t)Sm (t), ⎪ ⎪ ⎪ dt N h (t) ⎪ ⎪ ∫ ⎪ ⎪ dI (t) I (t − τ2 ) − t dm (s)ds ⎪ ⎩ m = bβ (t − τ2 ) h Sm (t − τ2 )e t −τ2 − dm (t)Im (t), dt Nh (t − τ2 ) Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
(2.1)
with
two
delays,
Physica
A
(2019)
123327,
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
3
where all constant parameters are positive, and Λm (t), β (t), dm (t) are positive, continuous and periodic function with period ω > 0. Let τ = max{τ1 , τ2 }, C := C ([−τ , 0], R5 ) and C + := C ([−τ , 0], R5+ ). For any φ = (φ1 , φ2 , φ3 , φ4 , φ5 ) ∈ C , denote ∑ ∥φ∥ = 5i=1 ∥φi ∥∞ with ∥φi ∥∞ = max−τ ≤θ≤0 |φi (θ )|. Then, (C , C + ) is an ordered Banach space and C + is a normal cone of C with nonempty interior in C . For any given continuous function u = (u1 , u2 , u3 , u4 , u5 ) : [−τ , σ ) → R5 with σ > 0, we define ut ∈ C by ut (θ ) = (u1 (t + θ ), u2 (t + θ ), u3 (t + θ ), u4 (t + θ ), u5 (t + θ )),
∀θ ∈ [−τ , 0],
for any t ∈ [0, σ ). Lemma 2.1. For any φ ∈ C + , system (2.1) has a unique nonnegative bounded solution u(t , φ ) on [0, ∞) with u0 = φ . Moreover, system (2.1) generates an ω-periodic semiflow Q (t) = ut : C + → C + , i.e., Q (t)φ = ut (φ ), ∀t ≥ 0, and Q := Q (ω) has a global compact attractor. Proof. For any φ ∈ C + , we define
⎛
Λh − c β (t) φ
φ1 (0) 1 (0)+φ2 (0)+φ3 (0)
⎞
φ5 (0) − dh φ1 (0)
⎜ φ1 (−τ1 ) ⎜ c β (t − τ1 ) φ (−τ1 )e−(dh +αh )τ1 − (dh + αh )φ2 (0) ⎜ φ1 (−τ1 )+φ2 (−τ1 )+φ3 (−τ1 ) 5 ⎜ αh φ2 (0) − dh φ3 (0) f (t , φ ) := ⎜ ⎜ φ2 (0) ⎜ Λ (t) − b β (t) φ (0) − dm (t)φ4 (0) m ⎜ φ1 (0)+φ2 (0)+φ3 (0) 4 ⎝ ∫ − t dm (s)ds φ (−τ ) bβ (t − τ2 ) φ (−τ )+φ2 (−τ2 )+φ (−τ ) φ4 (−τ2 )e t −τ2 − dm (t)φ5 (0) 1 2 2 2 3 2
⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎠
φ (s)(or φ (s))
1 2 where φ (s) is defined to be 0 if (φ1 (s), φ2 (s), φ3 (s)) = (0, 0, 0) for some s ∈ [−τ , 0]. Since f (t , φ ) is continuous 1 +φ2 (s)+φ3 (s) in (t , φ ) ∈ R × C + and f (t , φ ) is Lipschitz in φ on each compact subset of R × C + , it then follows from [20, Theorems 2.2.1 and 2.2.3] that system (2.1) has a unique solution u(t , φ ) on its maximal interval [0, σφ ) of existence with u0 = φ . Let φ = (φ1 , φ2 , φ3 , φ4 , φ5 ) be given. If φi (0) = 0 for some i ∈ {1, 2, 3, 4, 5}, then fi (t , φ ) ≥ 0. According to [21, Theorem 5.2.1 and Remark 5.2.1], the unique solution u(t , φ ) of system (2.1) with u0 = φ satisfies ut (φ ) ∈ C + for all t ∈ [0, σφ ). Using the first, fourth and fifth equations of (2.1), we have
dSh (t) dt dSm (t) dt dIm (t) dt
≤ Λh − dh Sh (t),
t ∈ [0, σφ ),
≤ Λm (t) − dm (t)Sm (t), ≤ bβ (t − τ2 )Sm (t − τ2 )e
t ∈ [0, σφ ), ∫ − tt−τ dm (s)ds 2
− dm (t)Im (t),
t ∈ [0, σφ ).
Thus, Sh (t), Sm (t) and Im (t) are bounded on [0, σφ ). As a consequence, we see from the second and third equations of (2.1) that Ih (t) and Rh (t) is bounded on [0, σφ ). Hence, [20, Theorem 2.3.1] implies that σφ = +∞. Define a semiflow Q (t) : C + → C + of (2.1) by Q (t)φ = ut (φ ), ∀φ ∈ C + . By the proof of [22, Lemma 3.5], we can show that {Q (t)}t ≥0 is an ω-periodic semiflow on C + . Let H(t) = Sh (t) + e(dh +αh )τ1 [Ih (t + τ1 ) + Rh (t + τ1 )]. By the first three equations of (2.1) we get dH(t) dt
= Λh − dh H(t). Λh
. Since Sh (t), Ih (t), Rh (t) ≥ 0 for t ≥ 0, it then follows from that for any ε > 0, dh ˆ there exists t > τ1 such that Λh Λh + ε, Ih (t + τ1 ) + Rh (t + τ1 ) ≤ + ε, ∀t > tˆ − τ1 , Sh (t) ≤
It is easy to see that limt →∞ H(t) =
dh
and hence, Ih (t), Rh (t) ≤ one attains
dh
Λh dh
+ ε for all t > tˆ. By applying the similar idea above to the fourth and fifth equations in (2.1),
lim sup(Sm (t) + edm τ2 Im (t + τ2 )) ≤ t →∞
¯m Λ dm
,
¯ m = maxt ∈[0,ω] Λm (t) and dm = mint ∈[0,ω] dm (t). This implies that all solutions of (2.1) are ultimately bounded. where Λ Consequently, the solution semiflow Q (t) : C + → C + is point dissipative. Note that for each t > τ , Q (t) is compact [20, Theorem 3.6.1]. It then follows from [23, Theorem 2.9] that Q := Q (ω) has a strong global attractor in C + . □ Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
4
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
3. Basic reproduction number In this section, we first introduce the basic reproduction number R0 for (2.1) by using the theory developed in [24] and then give a characterization of R0 . We remark that the theory of R0 for time-delayed periodic compartmental models is an extension of earlier ODE models [25,26]. It also applies to abstract delay differential equations (see, e.g., [27,28]). Letting Ih = Im = 0 in (2.1), we get
⎧ dS (t) ⎪ ⎨ h = Λh − dh Sh (t), dt
(3.1)
⎪ ⎩ dSm (t) = Λm (t) − dm (t)Sm (t). dt
∗ One can easily check that system (3.1) has a unique positive ω-periodic solution (Sh∗ , Sm (t)), ∀t ∈ R, where Sh∗ = ∗ (t) = Sm
[∫
t
∫s
Λm (s)e
∫ω
0 dm (r)dr
0
ds +
∫s
Λm (s)e e
0
0 dm (r)dr
∫ω
0 dm (s)ds
ds
] e−
∫t
0 dm (s)ds
−1
Λh dh
and
.
∗ Therefore, E0 (t) := (Sh∗ , 0, 0, Sm (t), 0) is a disease-free periodic solution of (2.1). Linearizing system (2.1) at E0 (t), the Ih and Im equations are given by
dIh (t) dt dIm (t) dt
= c β (t − τ1 )e−(dh +αh )τ1 Im (t − τ1 ) − (dh + αh )Ih (t), dh b
=
Λh
∫ − tt−τ dm (s)ds 2 I
β (t − τ2 )Sm∗ (t − τ2 )e
(3.2) h (t
− τ2 ) − dm (t)Im (t).
Let Z := C ([−τ , 0], R2 ) and Z + := C ([−τ , 0], R2+ ). Then (Z , Z + ) is an ordered Banach space equipped with the maximum norm and the positive cone Z + . For a function z(·) = (z1 (·), z2 (·)) ∈ C ([−τ , ∞), R2 ), we can define zt ∈ Z by zt (θ ) = (z1 (t + θ ), z2 (t + θ )), ∀θ ∈ [−τ , 0]. Let F : R → L(Z , R2 ) be a map and V (t) be a continuous 2 × 2 matrix function on R defined as follows: c β (t − τ1 )e−(dh +αh )τ1 φ2 (−τ1 )
( F (t)φ =
dh b
Λh
β (t − τ2 )Sm∗ (t − τ2 )e
)
∫t
− t −τ dm (s)ds 2 φ
1(
−τ2 )
, V (t) =
(
dh + αh 0
0 dm (t)
)
.
Let Φ (t , s), t ≥ s, be the evolution matrix associated with the following system: du(t) dt
= −V (t)u(t).
That is, for each s ∈ R, the 2 × 2 matrix Φ (t , s) satisfies
∂ Φ (t , s) = −V (t)Φ (t , s), ∀t ≥ s, and Φ (s, s) = I , ∂t where I is the 2 × 2 identity matrix. It then follows that ( −(d +α )(t −s) ) e h h 0 ∫ Φ (t , s) = . − t dm (r)dr 0
e
s
Recall that the exponential growth bound of Φ (t , s) is defined as
ω¯ (Φ ) := inf{ω˜ : ∃ M ≥ 1 such that ∥Φ (t + s, s)∥ ≤ Meω˜ t , ∀s ∈ R, t ≥ 0}. Observing that ∥Φ (t + s, s)∥ ≤ max{e−(dh +αh )t , e−dm t }, we have
ω¯ (Φ ) ≤ − min{dh + αh , dm }. Therefore, F (t) and V (t) satisfy the following assumptions: (H1) Each operator F (t) : Z → R2 is positive in the sense that F (t)Z + ⊆ R2+ ; (H2) The periodic matrix −V (t) is cooperative, and ω ¯ (Φ ) < 0. Let Cω be the Banach space of all ω-periodic functions from R to R2 , equipped with the maximum norm and the positive cone Cω+ := {u ∈ Cω : u(t) ≥ 0, ∀t ∈ R}. Following the procedure in [24, Sect. 2], we define a linear operator L on Cω by
[Lv](t) :=
∞
∫
Φ (t , t − s)F (t − s)v (t − s + ·)ds,
∀t ∈ R, v ∈ Cω .
0
Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
5
Thus, the basic reproduction number for (2.1) is R0 = r(L), the spectral radius of L. Note that R0 is biologically interpreted as the asymptotic per generation growth rate in periodic environments [29]. For structured populations in heterogeneous environments, Inaba [30] introduced a new definition of R0 based on the generation evolution operator, which unifies two definitions in [29,31] and has intuitively clear biological meaning. For any given t ≥ 0, let P(t) be the solution maps of (3.2), that is, P(t)ϕ = zt (ϕ ), t ≥ 0, where z(t , ϕ ) is the unique solution of (3.2) with z0 = ϕ ∈ Z . Then P := P(ω) is the Poincaré map associated with (3.2). Let r(P) be the spectral radius of P. By [24, Theorem 2.1], we have the following result. Lemma 3.1.
R0 − 1 has the same sign as r(P) − 1.
To find the positive solution, with the form of exponential function multiplied by periodic function, of the linear system (3.2) associated with the definition of R0 , we consider the following phase space: Y := C ([−τ2 , 0], R) × C ([−τ1 , 0], R), Y + := C ([−τ2 , 0], R+ ) × C ([−τ1 , 0], R+ ). Given a continuous function v : [−τ2 , ∞) × [−τ1 , ∞) → R2 , denote vt ∈ Y by
vt (θ ) = (v1 (t + θ1 ), v2 (t + θ2 )), ∀θ = (θ1 , θ2 ) ∈ [−τ2 , 0] × [−τ1 , 0], t ≥ 0. Following the same approach as Lemma 2.1, one can derive that for any φ ∈ Y + , system (3.2) admits a unique nonnegative ˆ as the solution map of (3.2) on the space Y , that is, bounded solution v (t , φ ) on [0, ∞) with v0 = φ . Then we define P(t)
ˆ φ = vt (φ ), P(t)
∀φ ∈ Y , t ≥ 0.
ˆ ω) be the Poincaré map of (3.2). Let r(P) ˆ be the spectral radius of P. ˆ By arguments similar to those in [18, Thus Pˆ := P( ˆ Lemma 3.8], we have r(P) = r(P). Now we present a useful result for our next arguments. For this, we first shows that ˆ is eventually strongly positive. the periodic semiflow P(t) Lemma 3.2. For any φ ∈ Y + \ {0}, the solution v (t , φ ) of (3.2) with v0 = φ satisfies vi (t , φ ) > 0 for all t > 2τ , i = 1, 2, ˆ φ ≫ 0 for all t > 3τ . and hence, P(t) Proof. For any given φ = (φ1 , φ2 ) ∈ Y + \ {0}, let v (t , φ ) = (v1 (t , φ ), v2 (t , φ )) be the solution of (3.2) satisfying v0 = φ . For notation simplicity, we denote (v1 (t), v2 (t)) := (v1 (t , φ ), v2 (t , φ )). Without loss of generality, we assume that φ1 > 0. We first show that there exists t0 ∈ [0, τ2 ] such that v2 (t0 ) > 0. Assume, by contradiction, that v2 (t) = 0 for all t ∈ [0, τ2 ]. Then we deduce from the second equation of (3.2) that
φ1 (t − τ2 ) = v1 (t − τ2 ) = 0,
∀t ∈ [0, τ2 ],
which contradicts the assumption that φ1 > 0. In view of dv2 (t) dt
=
dh b
β (t − τ2 )Sm∗ (t − τ2 )e
Λh ≥ −dm (t)v2 (t),
we have that v2 (t) ≥ v2 (t0 )e
−
∫t
v1 (t) = e−(dh +αh )t v1 (0) +
t0 dm (s)ds
∫ − tt−τ dm (s)ds 2 v
1 (t
− τ2 ) − dm (t)v2 (t)
> 0 for all t ≥ t0 . By the integral form of the first equation of (3.2), we get
t
∫
e−(dh +αh )(t +τ1 −s) c β (s − τ1 )v2 (s − τ1 )ds > 0,
∀t ≥ t0 + τ1 .
0
Since t0 ∈ [0, τ2 ], one sees (v1t , v2t ) ≫ (0, 0),
∀t > 3 τ ,
ˆ is strongly positive on Y for each t > 3τ . that is, vt (φ ) ≫ 0 for t > 3τ . This shows that P(t)
□
ˆ 0 ω) is strongly Choose an integer n0 > 0 such that n0 ω > 3τ . From the proof of Lemma 3.2, we see that Pˆ n0 = P(n n0 ˆ positive. Further, [20, Theorem 3.6.1] implies that P is compact. According to the Krein–Rutman theorem, as applied to ˆ n0 , we have λ = r(P) ˆ > 0, where λ is a simple eigenvalue the linear operator Pˆ n0 , together with the fact that r(Pˆ n0 ) = r(P) of Pˆ having a strongly positive eigenvector. A direct application of [32, Lemma 1] yields the following result. Lemma 3.3.
Let µ =
ln r(P)
ω
=
ˆ ln r(P) ω
. Then (3.2) has a solution eµt v ∗ (t) such that v ∗ (t + ω) = v ∗ (t) ≫ 0 for all t ∈ R.
Remark 3.1. Let v¯ (θ ) = eµθ v ∗ (θ ) for θ ∈ [−τ , 0]. By the uniqueness of solutions, one can see that v¯ (t) = eµt v ∗ (t) is a special solution of (3.2) on the state space Z + . Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
6
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
4. Global dynamics In this section, we first establish a threshold-type dynamics of system (2.1) in terms of R0 , and then prove the global attractivity of the positive equilibrium in the case where all the coefficients are constants and the IIP is ignored. 4.1. Threshold dynamics in terms of R0 Let X0 : = {φ = (φ1 , φ2 , φ3 , φ4 , φ5 ) ∈ C + : φ2 (0) > 0 and φ5 (0) > 0},
∂ X0 : = C + \ X0 = {φ = (φ1 , φ2 , φ3 , φ4 , φ5 ) ∈ C + : φ2 (0) = 0 or φ5 (0) = 0}. Then we have the following result. Theorem 4.1. If R0 > 1, then (2.1) admits a positive ω-periodic solution, and there is a real number η > 0 such that the solution (Sh (t), Ih (t), Rh (t), Sm (t), Im (t)) = (Sh (t , φ ), Ih (t , φ ), Rh (t , φ ), Sm (t , φ ), Im (t , φ )) of (2.1) with φ ∈ X0 satisfies lim inf(Ih (t), Im (t)) ⩾ (η, η). t →∞
Proof. In the case where R0 > 1, we have r(P) > 1 due to Lemma 3.1. Recall that Q (t) is the solution maps of (2.1) on C + , i.e. Q (t)φ = ut (φ ), t ≥ 0, and Q := Q (ω) is the Poincaré map associated with (2.1). By the form of (2.1), one easily sees that Q (t)X0 ⊂ X0 for all t ≥ 0. Now we prove that Q is uniformly persistent with respect to (X0 , ∂ X0 ). Let Pδ be the Poincaré map of the following perturbed linear periodic system:
( ⎧ dIh (t) ⎪ ⎪ = c β (t − τ1 ) 1 − ⎨ dt ⎪ dI ⎪ m (t)
⎩
dt
=
dh b
Λh + 3dh δ
4dh δ
)
Λh + 3dh δ
e−(dh +αh )τ1 Im (t − τ1 ) − (dh + αh )Ih (t), ∫ − tt−τ dm (s)ds 2 I
β (t − τ2 )(Sm∗ (t − τ2 ) − δ )e
(4.1) h (t
− τ2 ) − dm (t)Im (t).
Since limδ→0+ r(Pδ ) = r(P) > 1, there exists a sufficiently small δ > 0 such that
δ < min
{
min Sm (t), ∗
t ∈[0,ω]
Λh
}
dh
and r(Pδ ) > 1.
By using Lemma 3.3 and Remark 3.1, there is a positive ω-periodic solution vδ∗ (t) such that v¯ δ (t) = eµδ t vδ∗ (t) is a solution ln r(P ) of (4.1) on Z + , where µδ = ω δ > 0. ∗ ˆ θ ) = 0 for all θ ∈ [−τ , 0]. Then Q (t)M1 = ˆ where Sˆ ∗ (θ ) = S ∗ , S ∗ (θ ) = Sm∗ (θ ) and 0( , 0), Let M1 = (Sˆh∗ , 0ˆ , 0ˆ , Sm0 m0 h h ∗ ˆ , ∀t ≥ 0, and Q (M1 ) = M1 . Since limφ→M (Q (t)φ − Q (t)M1 ) = 0 uniformly for t ∈ [0, ω], there exists (Sh , 0ˆ , 0ˆ , Smt , 0) 1 δ1 = δ1 (δ ) such that
ˆ∗
∥Q (t)φ − Q (t)M1 ∥ ≤ δ,
∀t ∈ [0, ω], ∥φ − M1 ∥ < δ1 .
We further prove the following claim. Claim. lim supn→∞ ∥Q n (φ ) − M1 ∥ ≥ δ1 for all φ ∈ X0 . Suppose the claim is false, then lim supn→∞ ∥Q n (ψ ) − M1 ∥ < δ1 for some ψ ∈ X0 . Hence, there exists an integer n1 ≥ 1 such that ∥Q n (ψ ) − M1 ∥ < δ1 for all n ≥ n1 . For any t ≥ n1 ω, letting t = nω + t ′ with n ≥ n1 and t ′ ∈ [0, ω), we have
∥Q (t)ψ − Q (t)M1 ∥ = ∥Q (t ′ )(Q n (ψ )) − Q (t ′ )M1 ∥ < δ. One immediately sees that for t ≥ n1 ω,
⏐ ⏐ ⏐ ⏐ ⏐Sh (t , ψ ) − Λh ⏐ < δ, ⏐ dh ⏐
0 < Ih (t , ψ ), Rh (t , ψ ) < δ,
∗ Sm (t , ψ ) > Sm (t) − δ.
Then, for any t ≥ n1 ω + τ , Ih (t , ψ ) and Im (t , ψ ) satisfy
( ⎧ dIh (t) ⎪ ⎪ ≥ c β (t − τ1 ) 1 − ⎨ dt ⎪ ⎪ dIm (t)
⎩
dt
≥
dh b
Λh + 3dh δ
4dh δ
)
Λh + 3dh δ
e−(dh +αh )τ1 Im (t − τ1 ) − (dh + αh )Ih (t),
β (t − τ2 )(Sm∗ (t − τ2 ) − δ )e
∫ − tt−τ dm (s)ds 2 I
h (t
− τ2 ) − dm (t)Im (t),
where Ih (t) := Ih (t , ψ ) and Im (t) := Im (t , ψ ). Since Q (t)X0 ⊂ X0 , ∀t ≥ 0, we can choose a sufficiently small k > 0 such that (Ih (t , ψ ), Im (t , ψ )) ≥ kv¯ δ (t) = keµδ t vδ∗ (t), Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
∀t ≥ [n1 ω, n1 ω + τ ]. Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
7
Again, by the comparison theorem [21, Theorem 5.1.1], we can conclude that (Ih (t , ψ ), Im (t , ψ )) ≥ keµδ t vδ∗ (t),
∀ t ≥ n1 ω + τ .
Hence, limt →∞ Ih (t , ψ ) = limt →∞ Im (t , ψ ) = ∞, a contradiction. The above claim implies that M1 is an isolated invariant set for Q in C + and W s (M1 ) ∩ X0 = ∅, where W s (M1 ) is the stable set of M1 for Q . Define M∂ := {φ ∈ ∂ X0 : Q n (φ ) ∈ ∂ X0 , ∀n ≥ 0}. For any given ψ ∈ M∂ , Q n (ψ ) ∈ ∂ X0 , ∀n ≥ 0. Thus, for each n ≥ 0, we have Ih (nω, ψ ) = 0 or Im (nω, ψ ) = 0. In view of dIm (t) ≥ −(dh + αh )Ih (t) and ≥ −dm (t)Im (t), (4.2) dt dt we infer that if there exists some t0 ≥ 0 such that Ii (t0 ) > 0 for some i = h or i = m, then Ih (t) > 0 or Im (t) > 0 for all t ≥ t0 . It then follows that Ih (t , ψ ) ≡ 0 or Im (t , ψ ) ≡ 0 for each t ≥ 0. If Ih (t , ψ ) = 0 for all t ≥ 0, then from the third, ∗ fourth and fifth equations of (2.1), we get that Rh (t , ψ ) → 0, Sm (t , ψ ) − Sm (t) → 0, Im (t , ψ ) → 0 as t → ∞, and hence, Sh (t , ψ ) − Sh∗ → 0 as t → ∞. If Ih (t0 , ψ ) > 0 for some t0 ≥ 0, the second equation implies dIh (t)
Ih (t , ψ ) ≥ Ih (t0 , ψ )e−(dh +αh )(t −t0 ) > 0, ∀t ≥ t0 . Thus, there holds that Im (t , ψ ) = 0, ∀t ≥ t0 . Then from the first and second equation of (2.1), we see that Sh (t , ψ ) − Sh∗ → ∗ 0, Ih (t , ψ ) → 0 as t → ∞, and thus, Rh (t , ψ ) → 0, Sm (t , ψ ) − Sm (t) → 0 as t → ∞. Consequently, we obtain that ω(ψ ) = M1 for any ψ ∈ M∂ , and M1 cannot form a cycle for Q in ∂ X0 . By the acyclicity theorem on uniform persistence for maps (see, e.g., [33, Theorem 1.3.1 and Remark 1.3.1]), it follows that Q : C + → C + is uniformly persistent with respect to X0 . Note that there exists an equivalent norm for C such that for each t > 0, the map Q (t) : C + → C + is an α -contraction on C , where α is the Kuratowski measure of noncompactness (see, e.g., [33, Theorem 3.5.1] and [34, Theorem 4.1.1]). Further, Q n is compact for any integer n with nω > τ . Hence, one can deduce from [33, Theorem 1.3.10] that Q : X0 → X0 admits ⋆ ⋆ (t), Im (t)) is an ω-periodic a global attractor A0 , and Q has a fixed point φ ⋆ ∈ A0 . Then, u(t , φ ⋆ ) = (Sh⋆ (t), Ih⋆ (t), R⋆h (t), Sm solution of (2.1). From the definition of X0 , one easily sees that ⋆ ⋆ (t) ≥ 0, Im (t) > 0. Sh⋆ (t) ≥ 0, Ih⋆ (t) > 0, R⋆h (t) ≥ 0, Sm
We claim that there exists some t¯ ∈ [0, ω] such that Sh⋆ (t¯ ) > 0. If it is not the case, then Sh⋆ (t) ≡ 0 for all t ≥ 0, due to the periodicity of Sh⋆ (t). Then from the first equation of (2.1) we get 0 = Λh > 0, which is a contradiction. Since dSh⋆ (t) dt
(
≥ − c β (t)
⋆ Im (t)
)
Sh⋆ (t) + Ih⋆ (t) + R⋆h (t)
+ dh Sh⋆ (t),
one can derive Sh⋆ (t) > 0 for t ≥ t¯ . The periodicity of Sh⋆ (t) implies that Sh⋆ (t) > 0 for all t ≥ 0. By the third equation ⋆ (t) > 0 for all t ≥ 0. of (2.1) and Ih⋆ (t) > 0 for t ≥ 0, we have that R⋆h (t) > 0 for t ≥ 0. Similarly, we can show that Sm ⋆ ⋆ ⋆ ⋆ ⋆ Therefore, (Sh (t), Ih (t), Rh (t), Sm (t), Im (t)) is a positive ω-periodic solution of (2.1). Moreover, by similar arguments to those in the proof of [32, Theorem 5], we can obtain the practical uniform persistence, that is, there exists η > 0 such that lim inf(Ih (t , φ ), Im (t , φ )) ≥ (η, η). □ t →∞
The following theorem shows that when the incubation period in the host population is equal to zero, the infection will be cleared from the population if R0 < 1. Theorem 4.2. If R0 < 1 and τ1 = 0, then the disease-free periodic solution E0 (t) of system (2.1) is globally attracting in C + . Proof. For the case of R0 < 1 and τ1 = 0, we let Pϵ be the Poincaré map of the following auxiliary system:
⎧ dI (t) h ⎪ = c β (t)Im (t) − (dh + αh )Ih (t), ⎨ dt
⎪ ⎩ dIm (t) = dt
dh b
Λh − dh ϵ
∫ − tt−τ dm (s)ds 2 I
β (t − τ2 )(Sm∗ (t − τ2 ) + ϵ )e
Since limϵ→0+ r(Pϵ ) = r(P) < 1, we fix a sufficiently small ϵ ∈
vϵ∗ (t)
is a positive ω-periodic function ln r(P ) where µϵ = ω ϵ < 0. In view of dNh (t) dt
= Λh − dh Nh (t) and
=
dSm (t) dt
Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
(v1∗ϵ (t), v2∗ϵ (t))
(
(4.3) h (t
0,
− τ2 ) − dm (t)Im (t).
Λh dh
)
such that r(Pϵ ) < 1. In a similar manner, there
such that v¯ ϵ (t) = eµϵ t vϵ∗ (t) is a positive solution of (4.3) on Z + ,
≤ Λm (t) − dm (t)Sm (t), and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
8
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
there exists a sufficiently large integer n2 > 0 with n2 ω ≥ τ such that Nh (t) ≥
Λh dh
− ϵ, Sm (t) ≤ Sm∗ (t) + ϵ, ∀t ≥ n2 ω − τ .
Therefore, for any t ≥ n2 ω, we have
⎧ dI (t) h ⎪ ≤ c β (t)Im (t) − (dh + αh )Ih (t), ⎨ dt
⎪ ⎩ dIm (t) ≤ dt
dh b
Λh − dh ϵ
β (t − τ2 )(Sm∗ (t − τ2 ) + ϵ )e
∫ − tt−τ dm (s)ds 2 I
h (t
− τ2 ) − dm (t)Im (t).
Then the comparison theorem [21, Theorem 5.1.1] implies that (Ih (t), Im (t)) ≤ K v¯ ϵ (t) = Keµϵ t vϵ∗ (t),
∀t ≥ n2 ω,
with constant K > 0 satisfying (Ih (t), Im (t)) ≤ K v¯ ϵ (t) for t ∈ [n2 ω−τ , n2 ω]. As a consequence, limt →∞ (Ih (t), Im (t)) = (0, 0) Λ and thus limt →∞ Rh (t) = 0. Moreover, by using the first and fourth equations of (2.1), we have that limt →∞ Sh (t) = d h h ∗ and limt →∞ (Sm (t) − Sm (t)) = 0. □ Remark 4.1. By applying the perturbation theory of a globally stable fixed point (see [35, Theorem 2.2]) to the Poincaré map of (2.1), we can further prove that the above result still holds for R0 < 1 and sufficiently small delay τ1 > 0. 4.2. Global attractivity in the case of constant coefficients and τ1 = 0 In the case where β (t), Λm (t) and dm (t) are positive constants and τ1 = 0, system (2.1) reduces to the following autonomous system:
⎧ Sh (t) dSh (t) ⎪ ⎪ = Λh − c β Im (t) − dh Sh (t), ⎪ ⎪ dt N h (t) ⎪ ⎪ ⎪ ⎪ Sh (t) dIh (t) ⎪ ⎪ = cβ Im (t) − (dh + αh )Ih (t), ⎪ ⎪ dt N ⎪ h (t) ⎪ ⎨ dRh (t) = αh Ih (t) − dh Rh (t), ⎪ dt ⎪ ⎪ ⎪ dSm (t) Ih (t) ⎪ ⎪ = Λ m − bβ Sm (t) − dm Sm (t), ⎪ ⎪ ⎪ dt N h (t) ⎪ ⎪ ⎪ ⎪ dI (t) I (t − τ2 ) ⎪ ⎩ m = bβ h Sm (t − τ2 )e−dm τ2 − dm Im (t). dt Nh (t − τ2 ) Since Nh (t) = Sh (t) + Ih (t) + Rh (t), system (4.4) is equivalent to ⎧ dN (t) ⎪ ⎪ h = Λh − dh Nh (t), ⎪ ⎪ dt ⎪ ⎪ ⎪ Sh (t) dS ⎪ h (t) ⎪ ⎪ = Λh − c β Im (t) − dh Sh (t), ⎪ ⎪ dt N h (t) ⎪ ⎪ ⎨ dI (t) Sh (t) h = cβ Im (t) − (dh + αh )Ih (t), dt N ⎪ h (t) ⎪ ⎪ ⎪ dSm (t) Ih (t) ⎪ ⎪ = Λ m − bβ Sm (t) − dm Sm (t), ⎪ ⎪ ⎪ dt N h (t) ⎪ ⎪ ⎪ ⎪ dI (t) I (t − τ2 ) ⎪ ⎩ m = bβ h Sm (t − τ2 )e−dm τ2 − dm Im (t). dt Nh (t − τ2 )
(4.4)
(4.5)
0 It is easy to see that system (4.5) has a disease-free equilibrium E 0 = (Nh0 , Sh0 , 0, Sm , 0), where Nh0 = Sh0 = 0
Λh dh
0 and Sm =
Λm dm
.
We consider the infective equations of the linearization of (4.5) at E : dIh (t) dt dIm (t) dt
= c β Im (t) − (dh + αh )Ih (t), =
(4.6)
Λm β dh b −dm τ2 e Ih (t − τ2 ) − dm Im (t). Λh dm
According to [24, Corollary 2.1], the basic reproduction number R0 of system (4.5) can be expressed as
√ R0 = r(Fˆ V −1 ) =
Λm cbdh β 2 e−dm τ2 , Λh d2m (dh + αh )
Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
9
where Fˆ =
(
cβ 0
0 Λm β dh b −dm τ2 e Λ d h m
)
( and
V =
dh + αh 0
0 dm
)
.
Note that system (4.6) is cooperative and irreducible. By using [21, Theorem 5.1 and Corollary 5.2], we know that the linear stability of the trivial equilibrium of (4.6) is equivalent to that of the following ODE system
(
d dt
(
)
Ih (t) Im (t)
=
−(dh + αh )
cβ
Λm β dh b −dm τ2 e Λh dm
−dm
)(
Ih (t) Im (t)
)
( := M
Ih (t) Im (t)
)
.
(4.7)
The characteristic equation of (4.6) is H(λ) := (λ + dh + αh )(λ + dm ) −
Λm cbdh β 2 −dm τ2 −λτ2 e = 0. Λh dm
Therefore, we have max{Reλ : H(λ) = 0} < 0(> 0) if and only if s(M) < 0(> 0), where s(M) := max{Reλ : λ is an eigenvalue of M }. A straightforward computation yields the following result. s(M) < 0(> 0) if and only if R0 < 1(> 1).
Lemma 4.1.
In the case where R0 < 1, we see from Theorem 4.2 that the disease free equilibrium E 0 is globally attractive. For the ⋆ ⋆ case of R0 > 1, system (4.5) has a unique endemic equilibrium E ⋆ = (Nh⋆ , Sh⋆ , Ih⋆ , Sm , Im ) with Nh⋆ = ⋆ Im =
where p = of (4.5) is
Λh dh
qΛm e−dm τ2 Λm Ih⋆
c β dh
Λh
Λh dm (qIh⋆ + dm ) , (pqΛm e−dm τ2 + dh dm q)Ih⋆ + dh d2m dh d2m (R20 − 1) Ih⋆ = , pqΛm e−dm τ2 + dh dm q Sh⋆ =
,
dm (qIh⋆ + dm ) and q =
bβ dh
Λh
,
⋆ Sm =
Λm , qIh + dm ⋆
. Since Nh⋆ is a globally stable equilibrium for the first equation of (4.5), the limiting system
⎧ dSh (t) ⎪ ⎪ = Λh − pSh (t)Im (t) − dh Sh (t), ⎪ ⎪ dt ⎪ ⎪ ⎪ dI (t) ⎪ ⎪ ⎨ h = pSh (t)Im (t) − (dh + αh )Ih (t), dt
(4.8)
dSm (t) ⎪ ⎪ ⎪ = Λm − qIh (t)Sm (t) − dm Sm (t), ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎩ dIm (t) = qIh (t − τ2 )Sm (t − τ2 )e−dm τ2 − dm Im (t). dt
Notice that (4.8) has the same form as system (4.1) with τ1 = 0 in [15]. Hence, [15, Theorem 4.1] implies that the unique ⋆ ⋆ positive equilibrium (Sh⋆ , Ih⋆ , Sm , Im ) of (4.8) is globally asymptotically stable. Next, we will use the theory of internally chain transitive sets (see, e.g., [33, Chapter 1]) to show that E ⋆ is globally attractive for system (4.5). Let W := {φ = (φ1 , φ2 , φ3 , φ4 , φ5 ) ∈ C ([−τ2 , 0], R5+ ) : φ1 ≥ φ2 + φ3 }, and W0 := {φ = (φ1 , φ2 , φ3 , φ4 , φ5 ) ∈ W : φ3 = 0, φ5 = 0}. Then we have the following result. Theorem 4.3. If R0 > 1, then the endemic equilibrium E ⋆ is globally attractive for system (4.5) in W \ W0 . Proof. Let Φ (t) : W → W be the solution semiflow of (4.5), that is,
Φ (t)ψ = (Nht (ψ ), Sht (ψ ), Iht (ψ ), Smt (ψ ), Imt (ψ )),
∀t ≥ 0, ψ ∈ W .
Set Φ = Φ (1), then {Φ }n≥0 define a discrete-time dynamical system on W . For any given ϕ ∈ W , let ω(ϕ ) be the omega limit set of the discrete-time orbit {Φ n (ϕ )}n≥0 . The fact that limt →∞ Nh (t) = Nh⋆ implies that limn→∞ (Φ n (ϕ ))1 = Nˆh⋆ , where ˆ denotes the inclusion R → C ([−τ2 , 0], R) by x → xˆ , xˆ (θ ) = x, θ ∈ [−τ2 , 0]. Thus, there exists a subset ω ˜ of C ([−τ2 , 0], R4+ ) such that ω(ϕ ) = {Nˆh⋆ } × ω ˜. n
Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
10
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
For any φ = (φ1 , φ2 , φ3 , φ4 , φ5 ) ∈ ω(ϕ ), there exists a sequences nk → ∞ such that Φ nk (ϕ ) → φ as k → ∞. Then one has
Φ n |ω(ϕ ) (Nˆh⋆ , φ2 , φ3 , φ4 , φ5 ) = {Nˆh⋆ } × S n |ω˜ (φ2 , φ3 , φ4 , φ5 ),
∀(φ2 , φ3 , φ4 , φ5 ) ∈ ω, ˜ n ≥ 0,
where S n = S(n) and S(t) : C ([−τ2 , 0], R4+ ) → C ([−τ2 , 0], R4+ ) is solution maps of system (4.8). By [33, Lemma 1.2.1], ω(ϕ ) is an internally chain transitive set for Φ , and hence, ω˜ is an internally chain transitive set for S. ˆ } or ω˜ = {(Sˆ ⋆ , Iˆ⋆ , Sˆm⋆ , Iˆm⋆ )}. We now claim that ω˜ ̸= In the case where R0 > 1, it follows that either ω ˜ = {(Sˆh0 , 0ˆ , Sˆm0 , 0) h h 0 ˆ ˆ0 ˆ ˆ }. Then ω(ϕ ) = {(Nˆ 0 , Sˆ 0 , 0ˆ , Sˆm0 , 0) ˆ }. Thus, ˆ {(Sh , 0, Sm , 0)} if ϕ ∈ W \ W0 . Suppose for contradiction that ω˜ = {(Sˆh0 , 0ˆ , Sˆm0 , 0) h h limt →∞ (Ih (t), Im (t)) = (0, 0), and for any ϵ > 0, there exists a T = T (ϵ ) > 0 such that |Nh (t) − Nh0 | < ϵ, |Sh (t) − Sh0 | < ϵ 0 and |Sm (t) − Sm | < ϵ for all t ≥ T . Then for any t ≥ T + τ2 , one has
( ⎧ dIh (t) ⎪ ⎪ ≥ cβ 1 − ⎨ dt ⎪ dI ⎪ m (t)
⎩
dt
2dh ϵ
)
Im (t) − (dh + αh )Ih (t), Λh + dh ϵ dh Λm − dm ϵ −dm τ2 Ih (t − τ2 ) − dm Im (t). ≥ bβ e dm Λh + dh ϵ
(4.9)
Let s(Mϵ ) be the stability modulus of the following matrix:
( Mϵ =
−(dh + αh ) d
bβ d h
m
(
cβ 1 −
Λm −dm ϵ −dm τ2 e Λh +dh ϵ
2dh ϵ Λh +dh ϵ
) ) .
−dm
(
Λ
In light of Lemma 4.1, we have s(M) > 0 due to R0 > 1. Since limϵ→0+ s(Mϵ ) = s(M), we can fix an ϵ ∈ 0, min{ d h , h
such that s(Mϵ ) > 0. For t ≥ T + τ2 , we construct the following comparison linear system:
( ⎧ ¯ dIh (t) ⎪ ⎪ = cβ 1 − ⎨
2dh ϵ
) } dm
Λm
)
I¯m (t) − (dh + αh )I¯h (t), Λh + dh ϵ ¯ ⎪ ⎪ ⎩ dIm (t) = bβ dh Λm − dm ϵ e−dm τ2 I¯h (t − τ2 ) − dm I¯m (t). dt dm Λh + dh ϵ dt
(4.10)
By the same arguments as that for the stability and instability of (4.6) and (4.7) above, one gets that s(Mϵ ) > 0 implies that the trivial solution of (4.10) is unstable. According to (4.9) and the comparison principle, we have that Ii (t) → ∞ as t → ∞, where i = h or i = m, a contradiction to the boundedness of solutions. It then follows that ω ˜ = {(Sˆh⋆ , Iˆh⋆ , Sˆm⋆ , Iˆm⋆ )}, ⋆ ˆ ⋆ ˆ⋆ ˆ ⋆ ˆ⋆ ˆ and hence, ω(ϕ ) = {(Nh , Sh , Ih , Sm , Im )}. This implies that ⋆ ⋆ lim (Nh (t), Sh (t), Ih (t), Sm (t), Im (t)) = (Nh⋆ , Sh⋆ , Ih⋆ , Sm , Im ) . □
t →∞
5. Numerical simulations In this section, we carry out numerical simulations to validate analytical results and explore the effects of key parameters on disease dynamics. We let, using month (about 30.4 days) as the time unit, Λh = 0.8, dh = 0.02, αh = 0.004, b = 0.2, c = 0.011, dm = 1.0032, which are chosen from [1,18]. Motivated by [36], we suppose that the recruitment rate of mosquitoes from larvae is proportional to the biting rate, i.e., Λm (t) = k × β (t), where k = 200 and β (t) = 3.04(1 + 0.6 cos(π t /6)). As mentioned in [13,14], the IIP and EIP for malaria infection are about 7–30 days, and 8–14 days, respectively. Here they are taken to be τ1 = 10/30.4 and τ2 = 12/30.4. To numerically calculate R0 , we use the bisection method with the help of [37, Remark 3.2]. Using the above parameter values, we obtain R0 = 3.3197. In this case, the disease is persistent and a positive periodic solution with period 12 is observed (see Fig. 5.1), which is consistent with Theorem 4.1. If we consider system (2.1) by setting β (t) ≡ 3.04 with other parameters values unchanged, then R0 = 2.9233. This means that the time-averaged system may underestimate the malaria disease risk. Furthermore, taking the same parameters as above except that β (t) = 0.608(1 + 0.6 cos(π t /6)), we get R0 = 0.2969. In this case, the number of infectious humans and infectious mosquitoes eventually tends to 0, and the disease dies out (see Fig. 5.2). Next, we investigate the effects of seasonality on R0 and positive periodic solutions. Let β (t) = 3.04(1 + δ cos(π t /6)), 0 ≤ δ < 1, where δ is the strength of seasonal forcing. Then R0 is strictly increasing with respect to δ from Fig. 5.3. Moreover, the impact of seasonality on malaria transmission is shown in Fig. 5.4. It demonstrates that the peak of disease spread is positively related with the strength of seasonal forcing. Finally, in Fig. 5.5 we show that R0 is a decreasing function of τ1 (or τ2 ) on [0, 1] for fixed τ2 (or τ1 ). This implies that extending the incubation period in the human population or in the vector population with chemical measures will contribute to disease control. However, the impact of the length of the EIP and IIP on malaria transmission are different. As can be seen from Fig. 5.5, R0 is more sensitive to τ2 than τ1 . In addition, it is found that R0 has only a small variance when τ1 varies in [0, 1] for fixed τ2 (see the left plot in Fig. 5.5). This observation can be used for explaining why many malaria models only consider the effects of the EIP on the disease, ignoring the IIP [18,27]. Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
11
Fig. 5.1. Long term behavior of the infectious compartments when R0 > 1.
Fig. 5.2. Long term behavior of the infectious compartments when R0 < 1.
Fig. 5.3. R0 as a function of δ .
Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
12
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
Fig. 5.4. Relationship between the peak of disease compartments and δ .
Fig. 5.5. R0 as functions of τ1 and τ2 . Left: τ2 = 12/30.4; Right: τ1 = 10/30.4.
The impact of the length of the IIP and the EIP on the positive periodic solutions are demonstrated in Figs. 5.6 and 5.7, respectively. Both figures are obtained in the case of R0 > 1. The simulations in Figs. 5.6 and 5.7 indicate that although the periodic solutions for infectious variables vary as τ1 and τ2 vary, respectively, the peak values of periodic solutions are more sensitive to the incubation period τ2 than the incubation period τ1 . For instance, the variance of the periodic solution for infectious mosquitoes arised from the change of the IIP is tiny, see the right plot in Fig. 5.6. 6. Discussion In this paper, taking into account of the incubation periods of the parasites in the humans and mosquitoes, we have proposed a malaria transmission model with two delays and seasonality. Then we defined the basic reproduction number R0 for the model. By applying the persistence theory of dynamical systems, we showed the threshold dynamics of the model in terms of R0 . That is, the disease dies out if R0 < 1 and the delay τ1 is small, whereas when R0 > 1, there exists an endemic periodic solution and the disease will spread in the human and mosquito population. For the autonomous model with constant parameters and τ1 = 0, we obtained the explicit expression of R0 , and further used the theory of internally chain transitive sets to prove the global attractivity of the endemic equilibrium for the case of R0 > 1. Numerically, we explored the effects of the IIP, the EIP and seasonality on the disease transmission. Finally, we remark that there are several questions for future research. Since the EIP is exquisitely temperaturesensitive [38], the EIP delay may be modified to a time-dependent delay. It then would be mathematically much harder to characterize the dynamics of the disease. In addition, as mentioned in [9], the longer incubation period may make exposed humans and mosquitoes spread the parasites to different locations. Therefore, it is not trivial but natural to incorporate the population dispersal into the model. This leads to malaria models in the form of reaction diffusion equations. For such models, many interesting mathematical and biological problems deserved to be investigated, such as traveling wave solutions and spread speed of the disease. Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
13
Fig. 5.6. The impact of the IIP delay on the periodic solution.
Fig. 5.7. The impact of the EIP delay on the periodic solution.
Acknowledgments We are very grateful to two anonymous reviewers for their valuable comments and suggestions. References [1] N. Chitnis, J.M. Hyman, J.M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical mcodel, Bull. Math. Biol. 70 (2008) 1272–1296. [2] J.B. Gutierrez, M.R. Galinski, S. Cantrell, E.O. Voit, From within host dynamics to the epidemiology of infectious disease: Scientific overview and challenges, Math. Biosci. 270 (2015) 143–155. [3] M. Tanner, B. Greenwood, C.J.M. Whitty, et al., Malaria eradication and elimination: views on how to translate a vision into reality, BMC Med. 13 (2015) 167. [4] S.A. Means, R.J. Smith, The impact of human and vector distributions on the spatial prevalence of malaria in sub-Saharan Africa, J. Theoret. Biol. 409 (2016) 70–85. [5] Y. Xiao, X. Zou, Transmission dynamics for vector-borne diseases in a patchy environment, J. Math. Biol. 69 (2014) 113–146. [6] R. Ross, The Prevention of Malaria, John Murray, London, 1911. [7] C. Cosner, J.C. Beier, R.S. Cantrell, D. Impoinvil, L. Kapitanski, M.D. Potts, A. Troyo, S. Ruan, The effects of human movement on the persistence of vector-borne diseases, J. Theoret. Biol. 258 (2009) 550–560. [8] F. Forouzannia, A.B. Gumel, Mathematical analysis of an age-structured model for malaria transmission dynamics, Math. Biosci. 247 (2014) 80–94. [9] S. Ruan, D. Xiao, J.C. Beier, On the delayed Ross-Macdonald model for malaria transmission, Bull. Math. Biol. 70 (2008) 1098–1114. [10] Y. Xiao, X. Zou, On latencies in malaria infections and their impact on the disease dynamics, Math. Biosci. Eng. 10 (2013) 463–481. [11] G.A. Ngwa, W.S. Shu, A mathematical model for endemic malaria with variable human and mosquito populations, Math. Comput. Model. 32 (2000) 747–763.
Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,
14
Y. Zhang, S. Liu and Z. Bai / Physica A xxx (xxxx) xxx
[12] S.E. Bellan, The importance of age dependent mortality and the extrinsic incubation period in models of mosquito-borne disease transmission and control, PLoS One 5 (2010) e10165. [13] L.M. Beck-Johnson, W.A. Nelson, K.P. Paaijmans, A.F. Read, M.B. Thomas, O.N. Bjørnstad, The importance of temperature fluctuations in understanding mosquito population dynamics and malaria risk, R. Soc. Open Sci. 4 (2017) 160969. [14] A.L. Buczak, B. Baugher, E. Guven, L.C. Ramac-Thomas, Y. Elbert, S.M. Babin, S.H. Lewis, Fuzzy association rule mining and classification for the prediction of malaria in South Korea, BMC Med. Inform. Decis. Mak. (2015) 15:47. [15] L. Cai, X.-Z. Li, B. Fang, S. Ruan, Global properties of vector-host disease models with time delays, J. Math. Biol. 74 (2017) 1397–1423. [16] A. Abdelrazec, A.B. Gumel, Mathematical assessment of the role of temperature and rainfall on mosquito population dynamics, J. Math. Biol. 74 (2017) 1351–1395. [17] D.A. Ewing, C.A. Cobbold, B.V. Purse, M.A. Nunn, S.M. White, Modelling the effect of temperature on the seasonal population dynamics oftemperate mosquitoes, J. Theoret. Biol. 400 (2016) 65–79. [18] X. Wang, X.-Q. Zhao, A periodic vector-bias malaria model with incubation period, SIAM J. Appl. Math. 77 (2017) 181–201. [19] D.L. Smith, F.E. McKenzie, Statics and dynamics of malaria infection in Anopheles mosquitoes, Malaria J. (2004) 3–13. [20] J.K. Hale, S.M. Verduyn Lunel, Introduction To Functional Differential Equations, Springer, New York, 1993. [21] H.L. Smith, Monotone Dynamical Systems: An Introduction To the Theory of Competitive and Cooperative Systems, in: Mathematical Surveys and Monographs, vol. 41, American Mathematical Society, 1995. [22] Y. Lou, X.-Q. Zhao, A theoretical approach to understanding population dynamics with seasonal developmental durations, J. Nonlinear Sci. 27 (2017) 573–603. [23] P. Magal, X.-Q. Zhao, Global attractors and steady states for uniformly persistent dynamical systems, SIAM. J. Math. Anal. 37 (2005) 251–275. [24] X.-Q. Zhao, Basic reproduction ratios for periodic compartmental models with time delay, J. Dyn. Differ. Equ. 29 (2017) 67–82. [25] N. Bacaër, S. Guernaoui, The epidemic threshold of vector-borne diseases with seasonality, J. Math. Biol. 53 (2006) 421–436. [26] W. Wang, X.-Q. Zhao, Threshold dynamics for compartmental epidemic models in periodic environments, J. Dyn. Differ. Equ. 20 (2008) 699–717. [27] Z. Bai, R. Peng, X.-Q. Zhao, A reaction–diffusion malaria model with seasonality and incubation period, J. Math. Biol. 77 (2018) 201–228. [28] R. Wu, X.-Q. Zhao, A reaction–diffusion model of vector-borne disease with periodic delays, J. Nonlinear. Sci. 29 (2019) 29–64. [29] N. Bacaër, E.H.A. Dads, On the biological interpretation of a definition for the parameter R0 in periodic population models, J. Math. Biol. 65 (2012) 601–621. [30] H. Inaba, On a new perspective of the basic reproduction number in heterogeneous environments, J. Math. Biol. 65 (2012) 309–348. [31] O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R0 in the models for infectious disease in heterogeneous populations, J. Math. Biol. 28 (1990) 365–382. [32] X. Wang, X.-Q. Zhao, Dynamics of a time-delayed Lyme disease model with seasonality, SIAM J. Appl. Dyn. Syst. 16 (2017) 853–881. [33] X.-Q. Zhao, Dynamical Systems in Population Biology, second ed., Springer, New York, 2017. [34] J.K. Hale, Asymptotic Behavior of Dissipative Systems, in: Mathematical Surveys and Monographs, vol. 25, American Mathematical Society, Providence, RI, 1988. [35] H.L. Smith, P. Waltman, Perturbation of a globally stable steady state, Proc. Amer. Math. Soc. 127 (1999) 447–453. [36] Y. Lou, X.-Q. Zhao, A climate-based malaria transmission model with structured vector population, SIAM J. Appl. Math. 70 (2010) 2023–2044. [37] X. Liang, L. Zhang, X.-Q. Zhao, Basic reproduction ratios for periodic abstract functional differential equations (with application to a spatial model for Lyme disease), J. Dyn. Differ. Equ. 31 (2019) 1247–1278. [38] K.P. Paaijmans, A.F. Read, M.B. Thomas, Understanding the link between malaria risk and climate, Proc. Natl. Acad. Sci. USA 106 (2009) 13844–13849.
Please cite this article as: Y. Zhang, S. Liu https://doi.org/10.1016/j.physa.2019.123327.
and
Z. Bai,
A
periodic
malaria
model
with
two
delays,
Physica
A
(2019)
123327,