Dynamics analysis of a Zika–dengue co-infection model with dengue vaccine and antibody-dependent enhancement

Physica A 522 (2019) 248–273

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Dynamics analysis of a Zika–dengue co-infection model with dengue vaccine and antibody-dependent enhancement Liping Wang, Hongyong Zhao

∗

Department of Mathematics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, People’s Republic of China

highlights • • • • •

Study the effect of antibody-dependent enhancement and dengue vaccine on two diseases. We derive the basic and invasion reproduction numbers for Zika and dengue. The criteria of vaccination against dengue strategies are established for disease eradication. Simulations show that dengue vaccine has a negative effect on the control of Zika. Our model fits well the reported human Zika and dengue cases in Brazil.

article

info

Article history: Received 26 September 2018 Received in revised form 31 December 2018 Available online 31 January 2019 Keywords: Zika Dengue Co-infection Vaccine Antibody-dependent enhancement Stability

a b s t r a c t Mounting clinical and experimental evidence indicates that antibodies against dengue can enhance Zika infection through the antibody-dependent enhancement (ADE) mechanism. ADE is an alarming evolutionary development that makes it difficult to develop a vaccination campaign against dengue. To better understand the effect of dengue vaccine and ADE on a Zika and dengue outbreak, we formulate a Zika–dengue co-infection model that describes the joint dynamics of Zika and dengue. We derive the basic and invasion reproduction numbers, which are threshold values to identify the existence and stability of a disease-free state and two boundary equilibrium points, where only one strain (Zika or dengue) is present, and to identify the persistence of disease. Our analysis and numerical simulations suggest that although vaccination against dengue has a positive effect on the control of dengue, it has a negative effect on the control of Zika, and the increasing level of ADE induces a large number of accumulated Zika cases. Later, in this paper, criteria of vaccination against dengue strategies are established for disease eradication. Finally, the model is applied to simulate the accumulated cases of Zika and dengue in Brazil from February 19, 2016 to April 14, 2018. Based on the parameters estimated, we obtain that in Brazil the dengue epidemic is endemic whereas the Zika epidemic will eventually disappear. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Zika and dengue are mosquito-borne infectious diseases that are part of a global public health crisis. For both viruses, the main vector is Aedes aegypi [1,2]. Dengue is prevalent throughout the tropical and subtropical regions across at least 100 countries in Africa, the Americas and the Eastern Mediterranean [3,4]. It has increased over the past several decades, infecting ∗ Corresponding author. E-mail address: [email protected] (H. Zhao). https://doi.org/10.1016/j.physa.2019.01.099 0378-4371/© 2019 Elsevier B.V. All rights reserved.

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

249

an estimated 390 million people per year [5]. In recent years, Zika has become an emerging mosquito-borne pathogen. It was isolated from rhesus monkeys in the Zika forest, Uganda in 1947 [6]. Few human cases were reported until the first known outbreak, which occurred on Yap Island, Micronesia, in April through July 2007 [7]. Later, an outbreak occurred in French Polynesia between October 2013 and April 2014 [8]. Since then, the Zika virus is no longer a mild infection limited to Africa and Asia [9,10]. It has spread rapidly across continents, swept across Central and Southern America, and arrived in North America in 2016. The spread of the Zika virus led the World Health Organization (WHO) to declare the Zika outbreak a Public Health Emergency of International Concern on February 1, 2016 [11]. In the 2015–16 outbreak, up to 534,553 cumulative Zika cases were reported by countries and territories in the Americas. The dengue virus is known to be endemic in many of the regions involved in the current Zika outbreak. Co-infection cases of Zika and dengue have been reported in New Caledonia in 2014 [12], in Brazil in 2015 [13], in Colombia in 2016 [14,15]. With more than one-third of the world’s population living in areas at risk for infection of dengue virus, the co-circulation of Zika and dengue presents a potentially serious public concern. Furthermore, it is known that Zika shares up to 60% of nucleotide identity with dengue. Mounting clinical and experimental evidence indicates that Zika and dengue are highly cross-reactive [16–21], and that dengue virus sero-cross-reactivity drives antibody-dependent enhancement (ADE) of infection with Zika virus, that is, antibodies against dengue can enhance Zika infection. With respect to vaccine, ADE is an alarming evolutionary development in both viruses. With the increase in dengue infections in recent years, faster development of a vaccine becomes a serious concern [22], but development of a vaccination campaign is challenged by the possibility of ADE. Mathematical models have played an important role in gaining insights into the dynamics of interaction between host and vector [23–28]. Many models have been proposed to study the dynamic behaviors of distinct serotypes of dengue viruses. Feng et al. [29] proposed a two-strain model to simulate the spread of dengue virus between mosquitoes and humans and obtained permit coexistence in competing serotypes. Sriprom et al. [30] introduced a cross-protection dengue model that allowed for both symptomatic and asymptomatic infections. More studies of cross immunity have been applied to multiserotype dengue virus [31–34]. However, the evidence suggests that ADE is a salient feature in multi-serotype dengue virus [35,36]. The effect of ADE on the transmission dynamics of dengue was first studied by Ferguson et al. [37], who showed that ADE could promote the co-existence of strains and lead to large oscillations. Cummings et al. [38] further demonstrated that ADE-driven dynamics could provide a good qualitative approximation to the chaotic desynchronization of strains commonly observed within longitudinal datasets of numbers of clinical cases. Billings et al. [39] proposed an epidemic multistrain model only considering human populations with ADE to examine the effects of single-strain vaccine campaign. Hu et al. [40] used three deterministic models for dengue to gain insight into the roles of ADE and cross-immunity on the dynamics of dengue epidemics. There are very few models to investigate the transmission dynamics of co-infection of Zika and dengue. Tang et al. [41] first proposed an age-structured model of co-infection of Zika and dengue to study the implication of vaccination against dengue for Zika outbreak, but the dengue vaccine is not intuitively embodied in the basic reproduction number, and only numerical simulations were conducted to demonstrate the effects of vaccination against dengue on Zika outbreak. To our knowledge, research addressing the effect of ADE on the control and transmission of Zika virus in the presence of dengue vaccine is still at the preliminary stage. The main aim of this paper is to investigate the effect of both ADE and dengue vaccine on the transmission dynamics and interaction between Zika and dengue, and whether or not vaccination against dengue is beneficial to disease eradication. In our work, we will formulate an epidemiological model to consider the effect of ADE on the control and transmission of Zika virus in the presence of dengue vaccine, and obtain the basic reproduction number in conjunction with the ADE factor and dengue vaccine. In addition, based on the ideas in [42], we introduce the concept of invasion reproduction number in our mosquito-borne model. The invasion reproduction number is a threshold value that determines whether or not one strain can invade in a population when the other is already present. The paper is organized as follows. In Section 2, we formulate an epidemiological model with ADE and dengue vaccine to describe the joint dynamics of Zika and dengue. In addition, the basic properties on the positivity and boundness of solutions are discussed. In Section 3, we derive the existence of equilibria and the basic and invasion reproduction numbers that determine the dynamic outcomes of the interaction between Zika and dengue. A systematic qualitative analysis of the model is provided in Section 4, including stability of equilibria for two subsystems involving Zika or dengue only and for the co-infection system, disease persistence, and effect of dengue vaccine and ADE on the disease transmission. In Section 5, to support our theoretical predictions, some numerical simulations are given. In Section 6, the system is applied to simulate the Zika and dengue data in Brazil. A brief discussion is given in Section 7. 2. Model formulation To study the effect of dengue vaccine and ADE on the spread of Zika and dengue, in this section, we present the formulation of a model that examines the dynamics of the Zika–dengue co-infection in two populations: human host population Nh (t), mosquito vector population Nm (t). The human host population is divided into eleven categories: susceptible human to both Zika and dengue Sh (t), vaccinated human against dengue and permanent immunity for dengue but susceptible to Zika Vh (t), vaccinated human against dengue and permanent immunity for dengue but infected with Zika Iv z (t), infected human with Zika only but still susceptible to dengue Iz (t), infected human with dengue only but still susceptible to Zika Id (t), infected

250

L. Wang and H. Zhao / Physica A 522 (2019) 248–273 Table 1 Parameter descriptions. Parameter

Description

Λh Λm µh µm

Humans recruitment rate (Week−1 ) Mosquitoes recruitment rate (Week−1 ) Natural death rate of humans (Week−1 ) Natural death rate of mosquitoes (Week−1 ) The effective coverage rate of dengue vaccine among humans (Dimensionless) Biting rate of mosquitoes (Week−1 ) Transmission probability from mosquitoes with Zika to humans (Dimensionless) Transmission probability from mosquitoes with dengue to humans (Dimensionless) Transmission probability from humans with Zika to mosquitoes (Dimensionless) Transmission probability from humans with dengue to mosquitoes (Dimensionless) Recovery rate of infected humans with Zika only (Week−1 ) Recovery rate of infected humans with dengue only (Week−1 ) Recovery rate of infected humans with both Zika and dengue Idz (Week−1 ) Modification parameter or amplification factor (Dimensionless) Modification parameter or decrease recovery rate factor (Dimensionless) susceptibility index or ADE factor (Dimensionless)

p a

βmhz βmhd βhmz βhmd γz γd γdz η1 , η2 ϵ1 , ϵ2 φv , φ1 , φ2

human with Zika and dengue Idz (t), recovered human from Zika and permanent immunity for Zika but susceptible to dengue only Rz (t), recovered human from dengue and permanent immunity for dengue but susceptible to Zika only Rd (t), infected human with Zika and immune to dengue Ydz (t), infected human with dengue and immune to Zika Yzd (t) and recovered human from Zika and dengue and permanent immunity for Zika and dengue Rdz (t). Here, the total human population is given by Nh (t) = Sh (t) + Iz (t) + Id (t) + Idz (t) + Rz (t) + Rd (t) + Ydz (t) + Yzd (t) + Vh (t) + Iv z (t) + Rdz (t). The mosquito vector population is divided into three categories: susceptible mosquito to both Zika and dengue Sm (t), infected mosquito with Zika only Imz (t) and infected mosquito with dengue only Imd (t). So, Nm (t) = Sm (t) + Imz (t) + Imd (t). Let the effective coverage rate of dengue vaccine among humans be p. Motivated by the idea of modeling in [41], in this study, it is assumed that vaccinated humans Vh (t) against dengue are recruited into the population at a constant rate pΛh . Then susceptible humans Sh (t) are recruited into the population at a constant rate (1 − p)Λh . Susceptible humans Sh (t) acquire infection with Zika or dengue via the bite of infectious Zika or dengue mosquitoes during probing and feeding at a per I capita rate aβmhz INmz or aβmhd Nmd , and progress to the infectious class Iz (t) or Id (t), respectively. Furthermore, humans infected h

h

with Zika Iz (t) acquire infection with dengue at a per capita rate η1 aβmhd

Imd , Nh

and become individuals with Zika and dengue

co-infection Idz (t). Also, humans in the Idz (t) class are acquired via Id (t) infected with Zika at a per capita rate η2 aβmhz INmz . The h constant ηi > 1, i = 1, 2, is an amplification factor that captures the fact that individuals co-infected with both diseases are highly infectious. Individuals in Idz class either recover with permanent immunity for dengue only and progress to the Ydz class at a rate ϵ1 γd , or recover with permanent immunity for Zika only and enter the Yzd class at a rate ϵ2 γz or recover with permanent immunity for both Zika and dengue and become the Rdz class at a rate γdz . 0 < ϵi < 1, i = 1, 2, captures the fact that co-infection reduces both Zika and dengue recovery rate. Infected humans Iz (t) recover at a rate γz to join the class Rz (t) and have permanent immunity for Zika, but individuals in Rz (t) class are susceptible to dengue. Humans Rz (t) can I acquire infection with dengue at a per capita rate φ1 aβmhd Nmd and enter the class Yzd (t). Individuals in Yzd (t) class recover h at a recovery rate γd and enter the class Rdz (t). Similarly, Infected humans Id (t) enter the class Rd (t) at a recovery rate γd and have permanent immunity for dengue, but individuals in Rd (t) class are susceptible to Zika. Humans Rd (t) are infected with Zika at a per capita rate φ2 aβmhz INmz to join the class Ydz (t). Individuals in Ydz (t) class recover at a recovery rate γz to h

join the class Rdz (t). Vaccinated humans Vh (t) acquire infection with Zika at a per capita rate φv aβmhz INmz , and progress to the h infected humans Iv z (t). Then individuals in Iv z (t) class recover at a recovery rate γz to enter the class Rdz (t). There is a per capita natural mortality rate µh in all humans sub-population. Here, the susceptibility indices φv , φ1 and φ2 are positive real numbers that may mimic either cross-immunity (0 < φi < 1, i = v, 1, 2) or increased susceptibility (φi > 1, i = v, 1, 2) by ADE where is a phenomenon in which viral replication is increased rather than decreased by immune sera. In terms of φi > 1, i = v, 1, 2, the case φv > 1 represents the scenario where dengue vaccine induces an increase in Zika cases due to ADE [43]. The case φ1 > 1 (or φ2 > 1) represents the scenario where the pre-existing antibodies to previous Zika (or dengue) infection cannot neutralize but rather enhance the dengue (or Zika) infection [16]. Throughout this paper, we investigate the case of φi > 1, i = v, 1, 2. Susceptible mosquitoes are recruited at a constant rate Λ( m and move to)the infection ( class with ) Zika or dengue after biting infectious Zika or dengue humans at an average rate aβhmz

Iz +Idz +Ydz +Iv z Nh

or aβhmd

Id +Idz +Yzd Nh

, and later die naturally at a

rate µm . We assume that once a mosquito is infected, it never recovers and it cannot be reinfected with a different strain of virus [29]. A flow diagram of co-infection of Zika and dengue transmission among humans and mosquitoes is depicted in Fig. 1. The meanings of all parameters are described in Table 1.

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

251

Fig. 1. Dynamical transmission flow diagram of Zika–dengue co-infection among humans and mosquitoes.

Based on the flow diagram in Fig. 1, we formulate a transmission model for co-infection of Zika and dengue transmission among humans and mosquitoes

⎧ Sh Sh dSh ⎪ ⎪ = (1 − p)Λh − aβmhz Imz − aβmhd Imd − µh Sh , ⎪ ⎪ dt Nh Nh ⎪ ⎪ ⎪ Iz Sh dIz ⎪ ⎪ ⎪ ⎪ dt = aβmhz Imz Nh − η1 aβmhd Imd Nh − γz Iz − µh Iz , ⎪ ⎪ ⎪ ⎪ dId Sh Id ⎪ ⎪ = aβmhd Imd − η2 aβmhz Imz − γd Id − µh Id , ⎪ ⎪ dt N N ⎪ h h ⎪ ⎪ I I dI ⎪ ⎪ ⎪ dz = η1 aβmhd Imd z + η2 aβmhz Imz d − (ϵ1 γd + ϵ2 γz + γdz )Idz − µh Idz , ⎪ ⎪ dt N N ⎪ h h ⎪ ⎪ dRz Rz ⎪ ⎪ = γ I − φ a β I − µ R , ⎪ z z 1 mhd md h z ⎪ ⎪ dt Nh ⎪ ⎪ ⎪ dR R d d ⎪ ⎪ = γd Id − φ2 aβmhz Imz − µh Rd , ⎪ ⎪ dt Nh ⎪ ⎪ ⎪ Rd ⎪ dYdz ⎪ ⎪ ⎨ dt = ϵ1 γd Idz + φ2 aβmhz Imz N − γz Ydz − µh Ydz , h

dYzd Rz ⎪ = ϵ2 γz Idz + φ1 aβmhd Imd − γd Yzd − µh Yzd , ⎪ ⎪ ⎪ dt N h ⎪ ⎪ dV ⎪ Vh h ⎪ ⎪ = pΛh − φv aβmhz Imz − µh Vh , ⎪ ⎪ dt N ⎪ h ⎪ ⎪ dI V ⎪ ⎪ ⎪ vz = φv aβmhz Imz h − γz Ivz − µh Ivz , ⎪ ⎪ dt N ⎪ h ⎪ ⎪ dRdz ⎪ ⎪ = γ I + γ Y ⎪ dz dz z dz + γd Yzd + γz Iv z − µh Rdz , ⎪ ⎪ dt ( ) ( ) ⎪ ⎪ ⎪ dSm Iz + Idz + Ydz + Iv z Id + Idz + Yzd ⎪ ⎪ = Λm − aβhmz Sm − aβhmd Sm − µm Sm , ⎪ ⎪ Nh ) Nh ⎪ dt ⎪ ( ⎪ ⎪ dImz Iz + Idz + Ydz + Iv z ⎪ ⎪ = aβhmz Sm − µm Imz , ⎪ ⎪ ⎪ dt Nh ) ⎪ ( ⎪ ⎪ dI I + Idz + Yzd ⎪ ⎪ ⎩ md = aβhmd d Sm − µm Imd , dt

Nh

(2.1)

252

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

in which the parameters of the model are described in Table 1. The total population sizes of humans and mosquitoes can be determined by dNh (t) dt

= Λh − µh Nh ,

dNm (t) dt

= Λm − µm Nm .

It is easily to see that both for the humans population and for the mosquitoes population the corresponding total population Λ m . This implies that without loss of generality, sizes are asymptotically constant: limt →∞ Nh (t) = µ h and limt →∞ Nm (t) = Λ µ we can assume Nh (t) =

Λh µh

and Nm (t) =

Λm µm

m

h

for all t ≥ 0. Next, we gives the basic property of system (2.1).

2.1. Basic property System (2.1) describes both humans and mosquitoes population during the course of disease epidemic. So, it is epidemiologically meaningful to explore the positivity and boundness of solutions of system (2.1). Lemma 2.1. Let x(t) = (Sh , Iz , Id , Idz , Rz , Rd , Ydz , Yzd , Vh , Iv z , Rdz , Sm , Imz , Imd ). If initial values x(0) are nonnegative, then all solutions x(t) of system (2.1) with initial values x(0) are nonnegative for all t ≥ 0. In particular, the solution x(t) of system (2.1) is positive for t > 0 in the existence interval of the solution if x(0) > 0. Moreover, the region Γ = {(Sh , Iz , Id , Idz , Rz , Rd , Ydz , Yzd , Vh , Ivz , Rdz , Sm , Imz , Imd ) | Sh ≥ 0, Iz ≥ 0, Id ≥ 0, Idz ≥ 0, Rz ≥ 0, Rd ≥ 0, Ydz ≥ 0, Yzd ≥ 0, Vh ≥ 0, Iv z ≥ 0, Rdz ≥ 0, Sm ≥ 0, Imz ≥ 0, Imd ≥ 0, 0 ≤ Sh + Iz + Id + Idz + Rz + Rd + Ydz + Yzd + Vh + Iv z + Rdz = Λh m , 0 ≤ Sm + Imz + Imd = Λ } is positively invariant for the system (2.1). µ µ m

h

The proof of Lemma 2.1 can be found in the Appendix. 3. Reproduction numbers One of the most important concepts in epidemiology is the reproduction number which determines the ability of an infectious disease invading a population. For mosquito-borne disease, it is defined to be the average number of secondary cases produced in a susceptible population for a specific virus, by one infectious individual (mosquito or human) during its lifetime as infectious. Since the generation of secondary cases requires two transmission processes in mosquito-borne disease, the definition of square root form was given in [44]. In this section, we derive the basic and invasion reproduction numbers for Zika and dengue. 3.1. Basic reproduction number and boundary equilibria The basic reproduction number determines whether both Zika and dengue can be invade a completely susceptible population. Firstly, through calculation, the system (2.1) always has a disease-free equilibrium E0 = (Sh0 , 0, 0, 0, 0, 0, 0, 0, Vh0 , 0, 0 m . Applying the next generation matrix method in [45,46], we can 0, Sm , 0, 0), where Sh0 = (1−µp)Λh , Vh0 = pµΛh , Sm0 = Λ µm h h calculate the basic reproduction number of system (2.1) as follows R0 = max Rz , Rd ,

{

}

(3.1)

where

√ Rz =

a2 βmhz βhmz Λm µh (1 + p(φv − 1))

(3.2)

Λh µ2m (γz + µh )

and

√ Rd =

a2 βmhd βhmd Λm µh (1 − p)

Λh µ2m (γd + µh )

.

(3.3)

The square roots in the expressions in (3.2) and (3.3) account for the two generations needed to complete the transmission cycle of each of the two diseases. Rz represents the number of secondary infections produced by an infected Zika mosquito or human during its lifetime as infectious in a completely susceptible population. The number of new humans infected with Zika is generated by an infected Zika mosquito contacting with susceptible and vaccinated humans at a rate its lifetime

1

µm

aβmhz (Sh0 +φv Vh0 ) Nh

during

, and the number of new mosquitoes infected with Zika is generated by an infected Zika human contacting

with mosquitoes at a rate z

0 aβhmz Sm Nh

during its lifetime

1

γz +µh

d

. Thus, the geometric mean of

aβmhz (Sh0 +φv Vh0 ) Nh

· µ1m and

0 aβhmz Sm Nh

1 · γz +µ

h

gives the form R shown as (3.2). The interpretation of R is similar. The system (2.1) has other two boundary equilibria: Zika-only (dengue-free) boundary equilibrium, Ez , and dengue-only (Zika-free) boundary equilibrium, Ed .

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

253

At the dengue-free boundary equilibrium Ez , the components Id , Idz , Rd , Ydz , Yzd and Imd are zero. Then this equilibrium is given by ∗ ∗ , 0), , Imz Ez = (Sh∗ , Iz∗ , 0, 0, R∗z , 0, 0, 0, Vh∗ , Iv∗z , R∗dz , Sm

(3.4)

where Sh∗ =

(1 − p)Λh

(1 − p)Λh λ∗h

, Iz∗ =

, R∗z =

(1 − p)γz Λh λ∗h

λ∗h + µh (γz + µh )(λ∗h + µh ) µh (γz + µh )(λ∗h + µh ) ∗ φ p Λ λ p Λ v h h h Vh∗ = , I∗ = , φv λ∗h + µh vz (γz + µh )(φv λ∗h + µh ) γz φv pΛh λ∗h Λm Λm λ∗m ∗ = R∗dz = , Sm∗ = ∗ , Imz , ∗ µh (γz + µh )(φv λh + µh ) λm + µm µm (λ∗m + µm )

, (3.5)

in which

λ∗h =

∗ aβmhz Imz

Nh

, λ∗m =

aβhmz (Iz∗ + Iv∗z )

(3.6)

Nh

be per-capita infection rate from mosquitoes to humans and from humans to mosquitoes at steady state Ez , respectively. Substituting (3.5) in (3.6), one can get

λ∗m =

aβhmz µh ((φv λ∗h + µh ) + pµh (φv − 1))λ∗h (γz + µh )(λ∗h + µh )(φv λ∗h + µh )

,

(3.7)

and λ∗h satisfies the following quadratic equation A(λ∗h )2 + Bλ∗h + C = 0,

(3.8)

where A = Λh µm φv (aβhmz µh + µm (γz + µh )) > 0, B = aβhmz Λh µ2h µm (1 + p(φv − 1)) + Λh µ2m µh (γz + µh ) + a2 βhmz βmhz Λm µ2h φv p(φv − 1)

+ Λh µ2m µh φv (γz + µh )(1 − (Rz )2 ),

(3.9)

C = Λh µ2h µ2m (γz + µh )(1 − (Rz )2 ). Assume Rz ≤ 1. Clearly, C ≥ 0 and B > 0. Eq. (3.8) has no positive root. If Rz > 1, that is, C2 < 0. Then Eq. (3.8) has a unique positive root

√ λ∗h =

−B +

B2 − 4AC

2A

.

(3.10)

So, the state Ez exists only for Rz > 1. Similarly, dengue-only (Zika-free) boundary equilibrium Ed can be given by ∗∗ ∗∗ ∗∗ Ed = (Sh∗∗ , 0, Id∗∗ , 0, 0, R∗∗ d , 0, 0, Vh , 0, 0, Sm , 0, Imd ),

(3.11)

where (1 − p)Λh − (γd + µh )Id∗∗ Λ2h µ2m ((Rd )2 − 1) γd ∗∗ , Sh∗∗ = , R∗∗ I , d = aβhmd µh (aβmhd Λm + Λh µm ) µh µh d aβhmd Λm µh Id∗∗ pΛh Λh Λm ∗∗ = , Sm∗∗ = . ∗∗ , Imd = µh Λh µm + aβhmd µh Id µm (Λh µm + aβhmd µh Id∗∗ )

Id∗∗ = Vh∗∗

(3.12)

So the state Ed exists only for Rd > 1. 3.2. Invasion reproduction number Because the basic reproduction number only provides information in a population without any strains, it is important to consider another threshold value, the invasion reproduction number. Note that the invasion reproduction number is the average number of secondary infections caused by one infected individual in a population that is completely susceptible to one strain, but already has the other strain. In this part, we give the invasion reproduction number for each disease. That is, the reproduction number for Zika when dengue is endemic, and conversely, the reproduction number for dengue when Zika is endemic. The concept of invasion reproduction number was given in [42].

254

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

3.2.1. Invasion reproduction number for dengue The invasion reproduction number for dengue is derived in a population where Zika is already established at the endemic equilibrium. Consider the dengue-free boundary equilibrium Ez of system (2.1). It follows from [46] that the matrices Fd and Vd (corresponding to the new infection and remaining transfer terms, respectively) are given by

⎡

aβmhd Imd

Sh

⎤ ⎡

⎢ ⎥ Nh ⎢ ⎥ ⎢ ⎥ Iz ⎢ η1 aβmhd Imd ⎥ ⎢ ⎥ N h ⎢ ⎥, Fd = ⎢ ⎥ Rz ⎢ φ1 aβmhd Imd ⎥ ⎢ ⎥ Nh ⎢ ⎥ ( ) ⎣ ⎦ Id + Idz + Yzd aβhmd Sm

η2 aβmhz Imz

Id

+ γd Id + µh Id

⎤

⎥ ⎢ Nh ⎥ ⎢ ⎥ ⎢ Id ⎢ − η a β I + ( ϵ γ + ϵ γ + γ )I + µ I Vd = ⎢ 2 mhz mz 1 d 2 z dz dz h dz ⎥ . ⎥ N h ⎥ ⎢ ⎦ ⎣ − ϵ2 γz Idz + γd Yzd + µh Yzd µm Imd

Nh

Calculating the derivative of (Id , Idz , Yzd , Imd )T , then substituting dengue-free boundary equilibrium Ez into the variables, we can obtain

⎡ ⎢ ⎢ Fd = ⎢ ⎢ ⎣ ⎡

aβmhd Sh∗

0

0

0

0

0

0

0

0

0

∗ aβhmd Sm Nh

∗ aβhmd Sm Nh

∗ aβhmd Sm Nh

η2 aβmhz Nh

−

⎢ Vd = ⎢ ⎣

I∗

η1 aβmhd Iz∗ Nh

φ1 aβmhd R∗z Nh

∗ η2 aβmhz Imz

Nh

0 0

⎥ ⎥ ⎥, ⎥ ⎦

0

+ γd + µh

mz

⎤

Nh

(3.13)

0

0

0

⎤

ϵ1 γd + ϵ2 γz + γdz + µh −ϵ2 γz

0

γd + µh

0 0

⎥ ⎥. ⎦

0

0

µm

Then, the next generation matrix for dengue can be expressed by

⎡ Fd Vd−1

(

where r1 =

⎢ ⎢ ⎢ =⎢ ⎢ ⎢ ⎣

0

0

0

0

0

0

0

0

0

∗ aβhmd Sm r1 Nh

∗ (γ +µ +ϵ γ ) aβhmd Sm d h 2 z Nh (ϵ1 γd +ϵ2 γz +γdz +µh )(γdz +µh )

∗ aβhmd Sm Nh (γd +µh )

I∗

)

ϵ1 γd +ϵ2 γz +γdz +µh +η2 aβmhz Nmz (γd +µh )+ϵ2 γz η2 aβmhz h ( ) I∗ η2 aβmhz Nmz +γd +µh (ϵ1 γd +ϵ2 γz +γdz +µh )(γd +µh )

∗ Imz Nh

aβmhd Sh∗

µm Nh η1 aβmhd Iz∗ µm Nh φ1 aβmhd R∗z µm Nh

0

⎤ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎦

.

h

Rdz , the spectral radius of the next generation matrix Fd Vd−1 , is given by Rdz = ρ (Fd Vd−1 ) =

√

Rd1 + Rd2 + Rd3 ,

(3.14)

where

Rd1 =

) ) I∗ I∗ ϵ1 γd + ϵ2 γz + γdz + µh + η2 aβmhz Nmz (γd + µh ) + ϵ2 γz η2 aβmhz Nmz h h ( ) , ∗ Imz 2 Λh µm η2 aβmhz N + γd + µh (ϵ1 γd + ϵ2 γz + γdz + µh )(γd + µh )

∗ ∗ a2 βhmd βmhd µ2h Sm Sh

((

h

Rd2 = Rd3 =

∗ a2 βhmd βmhd µ2h Sm η1 Iz∗ (γd + µh + ϵ2 γz )

Λ2h 2

µm (ϵ1 γd + ϵ2 γz + γdz + µh )(γd + µh )

∗ a βhmd βmhd µ2h Sm φ1 R∗z

Λ2h µm (γd + µh )

,

(3.15)

.

∗ ∗ in which Sh∗ , Iz∗ , R∗z , Sm and Imz are defined in (3.5). Rdz is the invasion reproduction number for dengue in a population where Zika is already established at the endemic equilibrium. We know Sh , Iz and Rz are susceptible for dengue. So, the threshold quantity Rdz is associated with the dengue

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

255

disease transmission between humans in the Sh class and mosquitoes, humans in the Iz class and mosquitoes and humans in the Rz class and mosquitoes. The quantity Rd1 is associated with dengue transmission between humans in the Sh class and mosquitoes. Next, we give the biological interpretation of Rd1 . Susceptible mosquitos acquire infection with dengue by an effective contact with either the infected human with dengue only but still susceptible to Zika Id , or infected human with both Zika and dengue Idz or infected human with dengue and immune to Zika Yzd . The total number of infected mosquitoes with dengue is given by infectious human Id , Idz and Yzd . So, the number of new mosquitoes infected with dengue is generated by an infected human with dengue only but still susceptible to ∗ aβ Sm 1 . The number of mosquitoes Zika (Id ) contacting with mosquitoes at a rate hmd during its infectious period I∗ N η2 aβmhz

h

mz Nh

+γd +µh

infected with generated by an infected human with Zika and dengue (Idz ) during infectiousness period is given by the product ∗ aβ Sm , the probability that a human in the Id class will become infected human with both Zika and of the infection rate hmd N h

dengue (Idz )

∗ /N η2 aβmhz Imz h

I∗ η2 aβmhz Nmz h

and the average duration of the infectious lifetime of a human in the Idz class

+γd +µh

1

ϵ1 γd +ϵ2 γz +γdz +µh

.

The number of mosquitoes infected with dengue generated by an infected human with dengue and immune to Zika (Yzd ) ∗ aβ Sm during infectious period is given by the product of the infection rate hmd , the probability that a human in the Idz class will Nh ϵ γz and the average duration of the infectious become infected human with dengue and immune to Zika (Yzd ) ϵ γ +ϵ γ2 +γ +µ 1 d

2 z

dz

h

1 . Thus, in a susceptible mosquito population, the average number of new mosquitoes lifetime of a human in the Yzd class γ +µ d h infected with dengue generated by an infected human with dengue ( Id , Idz and Yzd ) is given by

∗ aβhmd Sm

Nh

·

1

η 2 aβ

aβhmd Sm

∗ Imz mhz N h

Nh

∗ aβhmd Sm

η2 aβmhz Imz /Nh

Nh

∗

∗

+

+ γd + µh

+

·

∗

η2 aβmhz INmz + γd + µh h

·

·

∗ η2 aβmhz Imz /Nh

η 2 aβ

∗ Imz mhz N h

+ γd + µh

·

1

ϵ1 γd + ϵ2 γz + γdz + µh

ϵ2 γz 1 · . ϵ1 γd + ϵ2 γz + γdz + µh γd + µh

(3.16)

The number of new humans infected with dengue is generated by an infected dengue mosquito contacting with susceptible aβmhd S ∗

humans at a rate N h during its lifetime µ1 . Thus, in a susceptible human population, the average number of new humans m h infected with dengue generated by an infected mosquito with dengue is given by aβmhd Sh∗ Nh

·

1

µm

.

(3.17)

Thus, the geometric mean of (3.16) and (3.17) gives the Rd1 in (3.15). The quantity Rd2 is associated with dengue transmission between humans in the Iz class and mosquitoes. When a typical infected with dengue person during the entire period of infection is introduced into a population made up infected humans with Zika only but still susceptible to dengue, secondary dengue infections include the humans in Idz class and Yzd class. So, similar the above analysis, in a susceptible mosquito population, the average number of new mosquitoes infected with dengue generated by an infected human (Idz and Yzd ) is given by ∗ aβhmd Sm

Nh

·

1

ϵ1 γd + ϵ2 γz + γdz + µh

+

∗ aβhmd Sm

Nh

·

ϵ2 γz 1 · . ϵ1 γd + ϵ2 γz + γdz + µh γd + µh

(3.18)

In a population made up infected humans with Zika only but still susceptible to dengue, the average number of humans infected with dengue generated by an infected mosquito with dengue is given by

η1 aβmhd Iz∗ Nh

·

1

µm

.

(3.19)

Thus, the geometric mean of (3.18) and (3.19) gives the Rd2 in (3.15). The quantity Rd3 in (3.15) has similar biological explanation.

3.2.2. Invasion reproduction number for Zika Similarly, consider the Zika-free boundary equilibrium Ed of system (2.1), then Rzd , invasion reproduction number for Zika, is given by Rzd =

√

Rz1 + Rz2 + Rz3 + Rz4 ,

(3.20)

256

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

here

Rz1 =

) ) I ∗∗ I ∗∗ ϵ1 γd + ϵ2 γz + γdz + µh + η1 aβmhd Nmd (γz + µh ) + ϵ1 γd η1 aβmhd Nmd h h ( ) , ∗∗ Imd 2 Λh µm η1 aβmhd N + γz + µh (ϵ1 γd + ϵ2 γz + γdz + µh )(γz + µh )

∗∗ ∗∗ a2 βhmz βmhz µ2h Sm Sh

((

h

Rz2 = Rz3 Rz4

∗∗ a2 βhmz βmhz µ2h Sm η2 Id∗∗ (γz + µh + ϵ1 γd )

Λ2h µm (ϵ1 γd + ϵ2 γz + γdz + µh )(γz + µh ) ∗∗ a2 βhmz βmhz µ2h Sm φ2 R∗∗ d = , 2 Λh µm (γz + µh ) ∗∗ a2 βhmz βmhz µ2h Sm φv Vh∗∗ = , 2 Λh µm (γz + µh )

,

(3.21)

∗∗ ∗∗ ∗∗ in which Sh∗∗ , Id∗∗ , R∗∗ d , Vh , Sm and Imd are defined in (3.12). Rzd is the invasion reproduction number for Zika in a population where dengue is already established at the endemic equilibrium. We know Sh , Vh , Id and Rd are susceptible for Zika. So, the threshold quantity Rzd is associated with the Zika disease transmission between humans (Sh , Vh , Id and Rd ) and mosquitoes. The quantity Rz1 measures Zika transmission between humans in the Sh class and mosquitoes. The quantity Rz2 measures Zika transmission between humans in the Id class and mosquitoes. The quantity Rz3 is associated with Zika transmission between humans in the Rd class and mosquitoes. The quantity Rz4 is associated with Zika transmission between humans in the Vh class and mosquitoes. The biological significance of Rzi in (3.21) is similar to Rdj , i = 1, 2, 3, 4, j = 1, 2, 3.

4. Stability of equilibria and disease persistence In this section, the main focus is to discuss the local and global behavior of system (2.1), including stability a disease-free state, conditions for the invasion of Zika (dengue) into a population in which dengue (Zika) is already established, and disease persistence. In the process, the stability of equilibria for two subsystems involving Zika or dengue only is also obtained. Furthermore, we present an effect of dengue vaccine and ADE on the disease transmission.

4.1. Local stability of disease-free equilibrium By analyzing the eigenvalues of the Jacobian matrix of (2.1) at disease-free equilibrium E0 , we have the following result regarding the stability of E0 . Theorem 4.1. The disease-free equilibrium E0 of system (2.1) is locally asymptotically stable if R0 < 1 and unstable if R0 > 1. Proof. The stability of the disease-free equilibrium E0 of system (2.1) is determine by the eigenvalues of the Jacobian matrix J(E0 ). the eigenvalues of J(E0 ) are −µh (multiplicity 5), −µm , −(ϵ1 γd + ϵ2 γz + γdz + µh ), −(γd + µh ), −(γz + µh ) (multiplicity 2) and the roots of the polynomials

λ2 + b1 λ + b2 = 0

(4.1)

λ2 + c1 λ + c2 = 0,

(4.2)

and

where b1 = µm + γd + µ(h > 0, ) b2 = µm (γd + µh ) 1 − (Rd )2 , c1 = µm + γz + µ(h > 0, ) c2 = µm (γz + µh ) 1 − (Rz )2 , If R0 < 1, i.e., Rd < 1 and Rz < 1, then clearly b2 > 0, and c2 > 0. The coefficients of Eqs. (4.1) and (4.2) are always positive. Therefore, the coefficients of Eqs. (4.1) and (4.2) satisfy Routh–Hurwitz conditions, that is, the two roots of Eq. (4.1) and the two roots of Eq. (4.2) have negative real parts. Thus, all eigenvalues of the Jacobian matrix J(E0 ) have negative real parts when R0 < 1, and Eq. (4.1) (or Eq. (4.2)) has at least one positive real part root when Rd > 1 (or Rz > 1). In conclusion, E0 is locally asymptotically stable if R0 < 1 and unstable if R0 > 1. □

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

257

4.2. Local stability of boundary equilibria

To give the stability of boundary equilibria Ez and Ed , we consider the following two subsystems firstly,

Zika subsystem

⎧ Sh dSh ⎪ ⎪ = (1 − p)Λh − aβmhz Imz − µh Sh , ⎪ ⎪ dt Nh ⎪ ⎪ ⎪ ⎪ Sh ⎪ dIz ⎪ = aβmhz Imz − γz Iz − µh Iz , ⎪ ⎪ ⎪ dt Nh ⎪ ⎪ ⎪ dRz ⎪ ⎪ ⎪ = γz Iz − µh Rz , ⎪ ⎪ dt ⎪ ⎪ ⎪ dVh Vh ⎪ ⎪ = pΛh − φv aβmhz Imz − µh Vh , ⎪ ⎨

dt Nh dI V v z h ⎪ ⎪ = φv aβmhz Imz − γz Ivz − µh Ivz , ⎪ ⎪ dt Nh ⎪ ⎪ ⎪ ⎪ dRdz ⎪ ⎪ ⎪ ⎪ dt = γz Ivz − µh Rdz , ⎪ ⎪ ( ) ⎪ ⎪ dSm Iz + Iv z ⎪ ⎪ Sm − µm Sm , = Λm − aβhmz ⎪ ⎪ ⎪ Nh ⎪ dt

(4.3)

⎪ ( ) ⎪ ⎪ dImz Iz + Iv z ⎪ ⎪ Sm − µm Imz , = aβhmz ⎩ dt

Nh

and

dengue subsystem

⎧ Sh dSh ⎪ ⎪ = (1 − p)Λh − aβmhd Imd − µ h Sh , ⎪ ⎪ dt N ⎪ h ⎪ ⎪ ⎪ dId Sh ⎪ ⎪ = aβmhd Imd − γd Id − µh Id , ⎪ ⎪ ⎪ dt N h ⎪ ⎪ ⎪ dR ⎪ ⎪ ⎨ d = γd Id − µh Rd , dt

(4.4)

dVh ⎪ ⎪ ⎪ = pΛh − µh Vh , ⎪ ⎪ dt ⎪ ⎪ ⎪ Id dSm ⎪ ⎪ = Λm − aβhmd Sm − µm Sm , ⎪ ⎪ ⎪ dt Nh ⎪ ⎪ ⎪ dI I ⎪ ⎪ ⎩ md = aβhmd Sm d − µm Imd . dt

Nh

Firstly, we state some results of Zika subsystem (4.3). Consider the feasible region of Zika subsystem (4.3) Γz =

{(Sh , Iz , Rz , Vh , Ivz , Rdz , Sm , Imz ) | Sh ≥ 0, Iz ≥ 0, Rz ≥ 0, Vh ≥ 0, Ivz ≥ 0, Rdz ≥ 0, Sm ≥ 0, Imz ≥ 0, 0 ≤ Λ m Sh + Iz + Rz + Vh + Iv z + Rdz = µ h , 0 ≤ Sm + Imz = Λ }. It follows from Lemma 2.1 that Γz is positively invariant. µ m

h

0 The subsystem (4.3) always has a disease-free equilibrium E0z = (Sh0 , 0, 0, Vh0 , 0, 0, Sm , 0), where Sh0 = pΛh (1−p)Λh Λm 0 0 z , V = , S = . When R > 1, the subsystem (4.3) has a unique endemic equilibrium, denote Ezz = m h µ µh µm ∗ h ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ (Sh , Iz , Rz , Vh , Iv z , Rdz , Sm , Imz ), where Sh , Iz , Rz , Vh , Iv z , Rdz , Sm and Imz are defined in (3.5). In the following, we

will present global stability of equilibria for subsystem (4.3). Theorem 4.2. If Rz < 1, then E0z is globally asymptotically stable for subsystem (4.3). Proof. We employ Lyapunov function taking advantage of the properties of the function g(x) = x − 1 − ln x. g(x) is positive in R+ except at x = 1, and g(x) = 0 if and only if x = 1. Therefore, consider the Lyapunov function V =

Sh0

γz + µh

( g

Sh Sh0

) +

Iz

γz + µh

+

Vh0

γz + µh

( g

Vh Vh0

) +

Iv z

γz + µh

+

Nh aβhmz

( g

Sm 0 Sm

Obviously, V is a positive definite function and attains zero at E0z . We need to show derivative of V along the solution of the subsystem (4.3) can be written as

) + dV dt

Nh 0 aβhmz Sm

Imz .

is negative definite. The orbital

258

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

⏐ dV ⏐ ⏐ dt ⏐

(Sh − Sh0 ) dSh

(γz + µh )Sh dt

+

1

dIz

γz + µh dt

+

(Vh − Vh0 ) dVh

(γz + µh )Vh dt

+

1

dIv z

+

0 Nh (Sm − Sm ) dSm

Nh

dImz

+ 0S 0 dt γz + µh dt aβhmz Sm aβhmz Sm m dt (4.3) ( ) ( ) (Sh − Sh0 ) Sh 1 Sh = (1 − p)Λh − aβmhz Imz − µ h Sh + aβmhz Imz − (γz + µh )Iz (γz + µh )Sh Nh γz + µh Nh ( ) ( ) (Vh − Vh0 ) Vh 1 Vh pΛh − φv aβmhz Imz − µh Vh + φv aβmhz Imz − (γz + µh )Ivz + (γz + µh )Vh Nh γz + µh Nh ( ( ) ) 0 Nh (Sm − Sm ) Iz + Idz + Ydz + Iv z + Λm − aβhmz Sm − µm Sm 0S aβhmz Sm Nh (m ( ) ) Nh Iz + Iv z + aβhmz Sm − µm Imz 0 aβhmz Sm Nh µh (Sh − Sh0 )2 aβmhz Sh0 Imz µh (Vh − Vh0 )2 φv aβmhz Vh0 Imz µm Nh (Sm − Sm0 )2 µm Nh Imz = − + − + − − 0S 0 (γz + µh )Sh (γz + µh )Nh (γz + µh )Vh (γz + µh )Nh aβhmz Sm aβhmz Sm m 0 2 0 2 0 2 ) µh (Vh − Vh ) µh (Sh − Sh ) µm Nh (Sm − Sm ) µm Nh ( − − − 1 − (Rz )2 Imz . = − 0S 0 (γz + µh )Sh (γz + µh )Vh aβhmz Sm a β S m hmz m ⏐ dV ⏐ z z z If R < 1, then dt (4.3) ≤ 0. So, E0 is locally stable. It can be seen that E0 is the largest positively invariant set contained ⏐ ⏐ in dV = 0. Thus, from the LaSalle Theorem, E0z is globally attractive. Therefore, E0z is globally asymptotically stable for dt (4.3) subsystem (4.3) if Rz < 1. □ =

Theorem 4.3. If Rz > 1, then subsystem (4.3) has a unique globally asymptotically stable endemic equilibrium Ezz in Γz \E0z . Proof. Set

(

( ) ( ) Iz Vh ∗ ∗ , V = I g , V = V g , 2 3 z h Sh∗ Iz∗ Vh∗ ( ) ( ) ( ) Iv z Sm Imz ∗ ∗ ∗ V 4 = Iv z g , Vk = S m g , k = 5, 6, V7 = Imz g ∗ . ∗ ∗

V1 = Sh∗ g

Sh

)

Iv z

Sm

Imz

Differentiating Vi along the solution of subsystem (4.3), we can obtain

⏐ ( )( ) Sh∗ Sh ⏐ = 1− (1 − p)Λh − aβmhz Imz − µh Sh dt ⏐(4.3) Sh Nh ) ( ) ( ∗ ∗ aβmhz Imz Sh Sh∗ aβmhz Imz Sh ∗ = 1− + µh Sh − − µ h Sh Sh Nh Nh ( ) ∗ ∗ ∗ )2 aβmhz Imz Sh S µh ( Imz Sh Imz =− Sh − Sh∗ + 1− h − ∗ ∗ + ∗ Sh Nh Sh Imz Sh Imz ( ) ∗ ∗ aβmhz Imz Sh Imz Sh Imz Sh ≤ ln ∗ − ∗ ∗ + ∗ Nh Sh Imz Sh Imz ( ) ∗ ∗ aβmhz Imz Sh Imz Sh Imz Imz Sh Imz ln ∗ ∗ − ln ∗ − ∗ ∗ + ∗ = Nh I S Imz Imz Sh Imz ( mz h ) ∗ ∗ aβmhz Imz Sh Imz Imz Imz Sh Imz Sh = − ln ∗ − ∗ ∗ + ln ∗ ∗ ∗

dV1 ⏐

Nh

⏐ ⏐ dt ⏐(4.3)

dV2 ⏐

Imz

Imz

Imz Sh

. = a17 G17 , ( )( ) Sh Iz∗ = 1− aβmhz Imz − (γz + µh )Iz Iz

∗

( = 1− =

Iz

≤

Nh

)(

aβmhz Imz Sh

Iz

Nh

∗ ∗ Sh aβmhz Imz

(

aβmhz Imz Sh Nh

. = a21 G21 .

Imz Sh ∗ S∗ Imz h

Nh ∗

Imz Sh

∗

(

Imz Sh ∗ S∗ Imz h

−

−

Iz Iz∗

− ln

∗ ∗ aβmhz Imz Sh Iz

+

Iz∗

Nh Iz∗ Imz Sh ∗ S∗ Iz Imz h

Imz Sh ∗ S∗ Imz h

−

)

Iz Iz∗

) +1 + ln

Iz Iz∗

)

(4.5)

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

259

Fig. 2. The weighted digraph (G , A) constructed for the subsystem (4.3).

Similarly, we can obtain

⏐ ( ) ∗ φv aβmhz Imz Vh∗ Imz Imz Imz Vh Imz Vh . ⏐ − ln − + ln = a37 G37 , ≤ ∗ ∗ ∗ V∗ ∗ V∗ dt ⏐(4.3) Nh Imz Imz Imz I mz h h ⏐ ( ) ∗ φ aβmhz Imz Vh∗ Imz Vh dV4 ⏐ Imz Vh Iv z Iv z . v ⏐ ≤ − ln − + ln = a43 G43 , ∗ V∗ ∗ V∗ dt ⏐(4.3) Nh Imz Imz Iv∗z Iv∗z h h ⏐ ( ) ∗ aβhmz Iz∗ Sm Iz Iz Iz S m Iz S m dVk ⏐ ⏐ ≤ − ln − + ln ∗ ∗ dt ⏐(4.3) Nh Iz∗ Iz∗ Iz∗ Sm Iz∗ Sm ( ) ∗ ∗ aβhmz Iv z Sm Iv z Iv z Iv z Sm Iv z S m . + − ln − + ln = a52 G52 + a64 G64 , k = 1, 2, ∗ ∗ Nh Iv∗z Iv∗z Iv∗z Sm Iv∗z Sm ⏐ ( ) ∗ aβhmz Iz∗ Sm Iz Sm dV7 ⏐ Iz Sm Imz Imz ⏐ ≤ − ln − + ln ∗ ∗ ∗ ∗ dt ⏐(4.3) Nh Iz∗ Sm Iz∗ Sm Imz Imz ( ) ∗ Iv z Sm Iv z S m Imz Imz aβhmz Iv∗z Sm . − ln − + ln = a75 G75 + a76 G76 , + ∗ ∗ ∗ ∗ ∗ ∗

dV3 ⏐

Nh

Iv z Sm

Iv z Sm

aβmhz I ∗ S ∗

Imz

Imz

φv aβmhz I ∗ V ∗

aβ

aβ

I∗S∗

I∗ S∗

mz h mz h z m where a17 = a21 = , a37 = a43 = , a52 = a75 = hmz , a64 = a76 = hmzN vz m . Function Gij can be Nh Nh Nh h defined accordingly. Consider a weight matrix A = (aij )7×7 and denote the corresponding weighted digraph as (G , A) with seven vertices. The associated digraph is shown in Fig. 2. One can readily verify that along each cycle C of such graph one gets ∑ (i,j)∈s(C ) Gij = 0, where s(C ) denotes the arc set of C . So, from Theorem 3.5 in [47], there exist ci ≥ 0, i = 1, 2, . . . , 7, such

∑7

⏐

dV ⏐ that V = i=1 ci Vi is a Lyapunov function for subsystem (4.3) and dt (4.3) ≤ 0. The relation ci can be derived from Theorems + + 3.3 and 3.4 in [47]: d (1) = 1 implies c2 a21 = c1 a17 , d (2) = 1 implies c5 a52 = c2 a21 , d+ (3) = 1 implies c4 a43 = c3 a37 , a d+ (4) = 1 implies c6 a64 = c4 a43 , d+ (5) = 1 implies c7 a75 = c5 a52 and d+ (6) = 1 implies c7 a76 = c6 a64 . Hence, c1 = a75 c7 ,

c2 =

a75 c , a21 7

a76 c , a37 7

c3 =

c4 =

a76 c , a43 7

c5 =

a75 c a52 7

17 ∗ ∗ βhmz Iz∗ Sm βhmz Iv∗z Sm a76 c . So, c1 c2 c4 ∗ S ∗ c7 , c3 ∗ V ∗ c7 and a64 7 βmhz Imz φ β I v mhz mz h h 7 7, that V i=1 ci Vi is a positive definite function. It can be verified that 0. Therefore, Ezz is globally asymptotically stable in the region Γz E0z . □

and c6 =

c5 = c6 = c7 . It follows from ci > 0, i = 1, 2,⏐ . . . , ⏐ = Ezz is the only invariant set in Γz \E0z where dV dt (4.3)

=

=

=

=

=

∑

\

For dengue subsystem (4.4), consider the feasible region Γd = {(Sh , Id , Rd , Vh , Sm , Imd ) | Sh ≥ 0, Id ≥ 0, Rd ≥ 0, Vh ≥ Λ m 0, Sm ≥ 0, Imd ≥ 0, 0 ≤ Sh + Id + Rd + Vh = µ h , 0 ≤ Sm + Imd = Λ }. It follows from Lemma 2.1 that Γd is positively invariant. µ m

h

0 The subsystem (4.4) always has a disease-free equilibrium E0d = (Sh0 , 0, 0, Vh0 , Sm , 0), where Sh0 =

Λm . When µm ∗∗ ∗∗

(1−p)Λh

µh ∗∗

, Vh0 =

pΛh

µh

,

∗∗ R > 1, the subsystem (4.4) has a unique endemic equilibrium, denote = (Sh , Id , Rd , Vh , Sm , Imd ), ∗∗ ∗∗ ∗∗ ∗∗ where Sh , Id , Rd , Vh , Sm and Imd are defined in (3.12). In the following, we will state global stability of equilibria for subsystem (4.4). 0 Sm

=

Edd

d

∗∗

∗∗

∗∗

∗∗

Theorem 4.4. Considering E0d and Edd as equilibria of subsystem (4.4) in Γd , we have the following results. (d1) If Rd < 1, then E0d is globally asymptotically stable for subsystem (4.4). (d2) If Rd > 1, then subsystem (4.4) has a unique globally asymptotically stable endemic equilibrium in the region Γd \E0d . Proof. (d1) Consider the Lyapunov function V = Sh0 g

(

Sh Sh0

)

+ Id + Vh0 g

(

Vh Vh0

) +

Nh (γd + µh ) aβhmd

( g

Sm 0 Sm

) +

Nh (γd + µh ) 0 aβhmd Sm

Imd .

260

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

(d2) Consider the Lyapunov function V =

Nh

( ∗∗ g

aβmhd Imd

Sh

+

(

Nh Imd aβhmd

+

Sh∗∗

∗∗

g I ∗∗ S ∗∗ d

Nh Id∗∗

)

m

∗∗ aβmhd Imd Sh

Imd ∗∗

( ∗∗ g

Imd

)

+ Vh∗∗ g

(

Id

) +

Id∗∗ Vh

)

Vh∗∗

Nh

( ∗∗ g

aβhmd Id

Sm

)

∗∗ Sm

.

The similar proof of Theorem 4.2, we can get that E0d is globally asymptotically stable for subsystem (4.4) if Rd < 1, and Edd is globally asymptotically stable for subsystem (4.4) in the region Γd \E0d if Rd > 1. □ Next, we investigate the dynamic behavior of full system (2.1) in the region Γ . The invasion reproduction numbers given in Section 3.2 and the basic reproduction numbers can be used to describe the stability of boundary equilibria Ez and Ed of system (2.1). Theorem 4.5. Let Rz > 1. Then Ez is locally asymptotically stable if Rdz < 1 and unstable if Rdz > 1. Proof. Let W s (Ez ) be the stable manifold of Ez , and

Γ1 = {(Sh , Iz , Id , Idz , Rz , Rd , Ydz , Yzd , Vh , Ivz , Rdz , Sm , Imz , Imd ) ∈ Γ |Id = Idz = Rd = Yzd = Ydz = Imd = 0}, Γ2 = {(Sh , Iz , Id , Idz , Rz , Rd , Ydz , Yzd , Vh , Ivz , Rdz , Sm , Imz , Imd ) ∈ Γ |Sh = Sh∗ , Iz = Iz∗ , Rz = R∗z , Vh = Vh∗ , ∗ ∗ Iv z = Iv∗z , Rdz = R∗dz , Sm = Sm , Imz = Imz }. It is easy to verify that Γ1 and Γ2 are invariant sets of system (2.1) at the boundary equilibrium Ez . It follows from the center manifold theorem that there exist invariant manifolds, denoted by W1 and W2 , of system (2.1) tangent to the Γ1 and Γ2 , respectively, and dim W1 = dim Γ1 , dim W2 = dim Γ2 . Similarly, Γ1 \E0 is also invariant set of system (2.1). There exists invariant manifolds, denoted by W3 , corresponding, and dim W3 = dim Γ1 \E0 = dim Γ1 . It follows from Theorem 4.3 that Ez is globally stable in Γ1 \E0 . Thus, W3 ⊂ W s (Ez ). Next, consider the linearization of system (2.1) on the subspace Γ2 . Computing the corresponding Jacobin matrix J2 is given by J2 = Fd − Vd , where Fd and Vd are defined by (3.13). Then from Theorem 2 in [46], we have s(J2 ) < 0 (> 0) ⇔ Rdz < 1 (> 1), where s(J2 ) denotes the maximum real part of all eigenvalues of the matrix J2 . If Rdz > 1, i.e., s(J2 ) > 0, then the equilibrium Ez is unstable. If Rdz < 1, then all the real parts of the eigenvalues of the matrix J2 are negative, and the subspace Γ2 is stable. Then the manifold W2 is also stable. Noting that dim W s (Ez ) ≥ dim W3 + dim W2 = 14, it follows from that dim W s (Ez ) = 14. That is, the equilibrium Ez is stable. So, when Rz > 1, Ez is locally asymptotically stable if Rdz < 1 and unstable if Rdz > 1. □ Theorem 4.6. Let Rd > 1. Then Ed is locally asymptotically stable if Rzd < 1 and unstable if Rzd > 1. We omit the proof of Theorem 4.6 as it is similar to proof of Theorem 4.5. 4.3. Disease persistence The analysis of existence of the interior equilibrium is too complicated. Here, we do not analyze it. In fact, the condition under which both Ez and Ed are unstable may provide insights into the coexistence of two pathogens. Next, we will prove the following analytic result regarding the disease persistence. Theorem 4.7. If Rz > 1, Rd > 1, Rdz > 1 and Rzd > 1, the system (2.1) is uniformly persistent, i.e., there exists a positive constant ϵ such that lim inf Sh ≥ ϵ, lim inf Iz ≥ ϵ, lim inf Id ≥ ϵ, lim inf Idz ≥ ϵ, lim inf Rz ≥ ϵ, t →+∞

t →+∞

t →+∞

t →+∞

t →+∞

lim inf Rd ≥ ϵ, lim inf Ydz ≥ ϵ, lim inf Yzd ≥ ϵ, lim inf Vh ≥ ϵ, lim inf Iv z ≥ ϵ, t →+∞

t →+∞

t →+∞

t →+∞

t →+∞

lim inf Rdz ≥ ϵ, lim inf Sm ≥ ϵ, lim inf Imz ≥ ϵ, lim inf Imd ≥ ϵ t →+∞

t →+∞

t →+∞

t →+∞

for system (2.1) with initial value Sh (0) > 0, Iz (0) > 0, Id (0) > 0, Idz (0) > 0, Rz (0) > 0, Rd (0) > 0, Ydz (0) > 0, Yzd (0) > 0, Vh (0) > 0, Iv z (0) > 0, Rdz (0) > 0, Sm (0) > 0, Imz (0) > 0, Imd (0) > 0. Proof. Define X0 = {(Sh , Iz , Id , Idz , Rz , Rd , Ydz , Yzd , Vh , Iv z , Rdz , Sm , Imz , Imd ) ∈ Γ | Iz > 0, Id > 0, Idz > 0, Rz > 0, Rd > 0, Ydz > 0, Yzd > 0, Iv z > 0, Rdz > 0, Imz > 0, Imd > 0},

∂ X 0 = Γ \X 0 .

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

261

It then suffices to show that system (2.1) is uniformly persistent with respect to (X0 , ∂ X0 ). Firstly, it follows from Lemma 2.1 that Γ and X0 are positively invariant for system (2.1), and X0 is relatively closed in Γ . Furthermore, from Lemma 2.1, we know that system (2.1) is point dissipative. For convenience later, denote x(t) = (Sh (t), Iz (t), Id (t), Idz (t), Rz (t), Rd (t), Ydz (t), Yzd (t), Vh (t), Iv z (t), Rdz (t), Sm (t), Imz (t), Imd (t)). Set M∂ = {x(0) | x(t) satisfies system (2.1) and x(t) ∈ ∂ X0 , ∀ t ≥ 0}, and denote

¯ ∂ 0 = {(Sh , 0, 0, 0, 0, 0, 0, 0, Vh , 0, 0, Sm , 0, 0) | Sh ≥ 0, Vh ≥ 0, Sm ≥ 0}, M ¯ ∂ z = {(Sh , Iz , 0, 0, Rz , 0, 0, 0, Vh , Ivz , Rdz , Sm , Imz , 0) | Sh ≥ 0, Vh ≥ 0, Sm ≥ 0, Iz > 0, Rz > 0, Ivz > 0, M Rdz > 0, Imz > 0}, ¯ M∂ d = {(Sh , 0, Id , 0, 0, Rd , 0, 0, Vh , 0, 0, Sm , 0, Imd ) | Sh ≥ 0, Vh ≥ 0, Sm ≥ 0, Id > 0, Rd > 0, Imd > 0}.

¯∂ = M ¯ ∂0 ∪ M ¯ ∂z ∪ M ¯ ∂ d . We will show that Let M ¯ ∂. M∂ = M ¯ ∂ ⊆ M∂ , then we only need to prove M∂ ⊆ M ¯ ∂. It is obvious that M ¯ ∂. Suppose not. Let x(t) be a solution of system (2.1) with the initial condition x(0). then x(t) ∈ M∂ for ∀ t ≥ 0, and x(t) ∈ /M Then x(t) = (Sh (t), Iz (t), Id (t), Idz (t), Rz (t), Rd (t), Ydz (t), Yzd (t), Vh (t), Iv z (t), Rdz (t), Sm (t), Imz (t), Imd (t)) satisfies one of the following conditions: (I) There exists t0 ≥ 0 such that at least one of Idz (t0 ), Ydz (t0 ) and Yzd (t0 ) is not zero; (II) If Idz (t) = Ydz (t) = Yzd (t) = 0, ∀t ≥ 0, then the following two conditions must satisfy: (i) There exists t0 ≥ 0 such that at least one of Iz (t0 ), Iv z (t0 ), Rz (t0 ), Rdz (t0 ) and Imz (t0 ) is not zero; (ii) There exists t0 ≥ 0 such that at least one of Id (t0 ), Rd (t0 ) and Imd (t0 ) is not zero. Assume x(t) satisfies condition (I). Without loss of generality, assume that Idz (t0 ) > 0, Ydz (t0 ) = Yzd (t0 ) = 0. From Lemma 2.1, Idz (t) > 0 for ∀t > t0 . From system (2.1), one can get, for ∀ t > t0 , Imz (t) = Imz (t0 )e−µm (t −t0 ) +

∫

t

(

aβhmz

Iz (s) + Idz (s) + Ydz (s) + Iv z (s) Nh

t

Imd (t) = Imd (t0 )e

−µm (t −t0 )

[ Iz (t) = Iz (t0 )exp −

+

aβmhz Nh

∫ 0t +

aβmhd Nh

∫

aβhmd

Id (s) + Idz (s) + Yzd (s) Nh

t0 t

Imd (s)ds − (γz + µh )(t − t0 )

)

Sm (s)e−µm (t −s) ds > 0,

Sm (s)e−µm (t −s) ds > 0,

]

t0

t

∫

(

)

[ Imz (s)Sh (s)exp −

aβmhd Nh

t

t

∫

]

Imd (ξ )dξ − (γz + µh )(t − s) ds > 0, s

[0 ] ∫ t aβmhz Id (t) = Id (t0 )exp − Imz (s)ds − (γd + µh )(t − t0 ) Nh t0 [ ] ∫ t ∫ t aβmhd aβmhd + Imd (s)Sh (s)exp − Imd (ξ )dξ − (γd + µh )(t − s) ds > 0, Nh Nh t0 s [ ] ∫ t φ1 aβmhd Rz (t) = Rz (t0 )exp − Imd (s)ds − µh (t − t0 ) Nh t0 [ ] ∫ t ∫ φ1 aβmhd t + γz Imd (ξ )ξ − µh (t − s) ds > 0, Iz (s)exp − Nh s t0 [ ] ∫ φ2 aβmhz t Rd (t) = Rd (t0 )exp − Imz (s)ds − µh (t − t0 ) Nh t0 [ ] ∫ ∫ t φ2 aβmhz t Imz (ξ )ξ − µh (t − s) ds > 0, + γd Id (s)exp − t0

Nh

s

262

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

Ydz (t) = Ydz (t0 )e−(γz +µh )(t −t0 ) + Yzd (t) = Yzd (t0 )e−(γd +µh )(t −t0 ) +

∫

∫

t

∫

t0 t

Nh

φ1 aβmhd Nh

φv aβmhz

Iv z (t) = Iv z (t0 )e−(γz +µh )(t −t0 ) + Rdz (t) = Rdz (t0 )e−µh (t −t0 ) +

φ2 aβmhz

∫

Nh

t0 t

Imz (s)Rd (s)e−(γz +µh )(t −s) ds > 0, Imd (s)Rz (s)e−(γd +µh )(t −s) ds > 0.

Imz (s)Vh (s)e−(γz +µh )(t −s) ds > 0,

t0

t

(γdz Idz (s) + γz Ydz (s) + γd Yzd (s) + γz Iv z (s))e−µh (t −s) ds > 0. t0

This means that x(t) ∈ / ∂ X0 for t > t0 . Similarly, It follows from condition (II) that x(t) ∈ / ∂ X0 for t > t0 , which contradicts ¯ ∂ . So M∂ = M ¯ ∂. the assumption x(t) ∈ ∂ X0 . Hence, we get M∂ ⊆ M It is clear that there are three equilibria E0 , Ez and Ed in M∂ , and E0 , Ez and Ed are isolated invariant sets in Γ . Next, we will show that W s (E0 ) ∩ X0 = ∅, W s (Ez ) ∩ X0 = ∅ and W s (Ed ) ∩ X0 = ∅, where W s (Ei ) denotes the stable manifold of Ei i = 0, d, z. That is, there exists a η > 0 such that for any solution Φt (x(0)) of system (2.1) with the initial value x(0) ∈ X0 , we have lim sup d(Φt (x(0)), E0 ) ≥ η, lim sup d(Φt (x(0)), Ez ) ≥ η, lim sup d(Φt (x(0)), Ed ) ≥ η, t →+∞

t →+∞

t →+∞

where d is a distant function in X0 . Here, we only verify lim supt →+∞ d(Φt (x(0)), Ez ) ≥ η. Suppose not, then lim supt →+∞ d(Φt (x(0)), Ez ) < ϵ for ∀ϵ > 0, namely, there exists T > 0, such that Sh∗ − ϵ ≤ Sh ≤ Sh∗ + ϵ , Iz∗ − ϵ ≤ Iz ≤ Iz∗ + ϵ , ∗ R∗z − ϵ ≤ Rz ≤ R∗z + ϵ , Vh∗ − ϵ ≤ Vh ≤ Vh∗ + ϵ , Iv∗z − ϵ ≤ Iv z ≤ Iv∗z + ϵ , R∗dz − ϵ ≤ Rdz ≤ R∗dz + ϵ , Sm − ϵ ≤ Sm ≤ Sm∗ + ϵ , ∗ ∗ Imz − ϵ ≤ Imz ≤ Imz + ϵ , 0 ≤ Id ≤ ϵ , 0 ≤ Idz ≤ ϵ , 0 ≤ Rd ≤ ϵ , 0 ≤ Ydz ≤ ϵ , 0 ≤ Yzd ≤ ϵ and 0 ≤ Imd ≤ ϵ for ∀t > T . From system (2.1), for ∀t > T , we have

⎧ S∗ − ϵ I∗ + ϵ dId ⎪ ≥ aβmhd h Imd − aβmhz mz Id − γd Id − µh Id , ⎪ ⎪ ⎪ dt Nh Nh ⎪ ⎪ ∗ ∗ ⎪ ⎪ dIdz I −ϵ I −ϵ ⎪ ⎪ ≥ aβmhd z Imd + aβmhz mz Id − γdz Idz − µh Idz , ⎨ dt Nh Nh ∗ dY R − ϵ ⎪ zd ⎪ ≥ φ1 aβmhd z Imd − γd Yzd − µh Yzd , ⎪ ⎪ ⎪ dt Nh ⎪ ( ) ⎪ ⎪ ⎪ dI S∗ − ϵ ⎪ ⎩ md ≥ aβhmd m Id + Idz + Yzd − µm Imd , dt Nh

(4.6)

Consider an auxiliary system as follows du dt

= Md (ϵ )u,

(4.7)

where vector u = (u1 , u2 , u3 , u4 )T and

⎡

∗ − aβNmhz (Imz + ϵ ) − γd − µh

0

0

aβmhd (Sh∗ Nh

− ϵ)

−γdz − µh

0

aβmhd ∗ (Iz Nh

− ϵ)

h

⎢ ⎢ ⎢ Md (ϵ ) = ⎢ ⎢ ⎣

aβmhz Nh

∗ (Imz − ϵ)

0 aβhmd ∗ (Sm Nh

−γd − µh

0 aβhmd ∗ (Sm Nh

− ϵ)

− ϵ)

aβhmd ∗ (Sm Nh

− ϵ)

⎤

⎥ ⎥ ⎥ ⎥. φ1 aβmhd ∗ ⎥ (R − ϵ ) z Nh ⎦ −µm

We can get Md (0) = Fd − Vd , where Fd and Vd are defined by (3.13). From Theorem 2 in [46], we have If Rdz > 1, then s(Md (0)) > 0. Since Md (ϵ ) is a continuous for small ϵ , there exists ϵ small enough, such that s(Md (ϵ )) > 0. Thus there is a positive eigenvalue of Md (ϵ ) with a positive eigenvector. Let u(t) = (u1 (t), u2 (t), u3 (t), u4 (t))T be a positive solution of the auxiliary system (4.7), which is strictly increasing with ui (t) → +∞ as t → +∞, i = 1, . . . , 4. Since the system (4.7) is a quasi-monotonic system, applying the comparison principle gives lim Id (t) = +∞, lim Idz (t) = +∞, lim Yzd (t) = +∞, lim Imd (t) = +∞.

t →+∞

t →+∞

t →+∞

t →+∞

From system (2.1), we can get lim Rd (t) = +∞, lim Ydz (t) = +∞.

t →+∞

t →+∞

This contradicts with our assumption. Then Ez is an isolated invariant set in Γ and W s (Ez ) ∩ X0 = ∅. Similarly, we obtain that if Rzd > 1, then Ed is an isolated invariant set in Γ and W s (Ed ) ∩ X0 = ∅, and if Rz > 1 and Rd > 1, then E0 is an isolated invariant set in Γ and W s (E0 ) ∩ X0 = ∅. Therefore, the system is uniformly persistent if Rz > 1, Rd > 1, Rdz > 1 and Rzd > 1. This completes the proof. □

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

263

4.4. Effect of dengue vaccine and ADE on the disease transmission It is clear that the control reproduction number R0 is a very useful quantity in disease control. An effective vaccination strategy should aim to achieve R0 less than one so that the disease will eventually be eradicated in the whole population. Due to the effect of ADE, the level of Zika infection actually get higher when the effective coverage rate of dengue vaccine p is increased. So, it is instructive to determine a criterion of vaccination against dengue in a community. We first give in the absence of vaccination, (i.e., p = 0) denoted by R00 = max{Rz0 , Rd0 }, √ the basic reproduction number √ a2 βmhz βhmz Λm µh

where Rz0 =

Rz = Rz0

Λh µ2m (γz +µh )

and Rd0 =

√

1 + p(φv − 1), Rd = Rd0

a2 βmhd βhmd Λm µh

Λh µ2m (γd +µh )

√

. So, we can rewrite Rz and Rd as follows

1 − p.

(4.8)

It can be shown that R is a increasing function of p, R is a decreasing function of p, R > R and R < R for 0 < p ≤ 1. This means that vaccination against dengue will have a negative effect on the control of Zika although have a positive effect on the control of dengue. To study effects of vaccination against dengue strategies, we assume that disease invades the population, that is R00 > 1. Next, we will give three cases as follows. Case (i) Assume Rz0 < 1 and Rd0 > 1. In order for the disease-free equilibrium E0 to be stable, a dengue vaccine rate need to be found such that Rz < 1, Rd < 1. For this purpose, we need give a critical value of vaccinate rate. Solving for p from Rz = 1 gives p = p∗1 , and from Rd = 1 gives p = p∗2 , where z

∗

p1 =

(

1

φv − 1

d

1 (Rz0 )2

)

− 1 , p∗2 = 1 −

1 (Rd0 )2

.

z

z0

d

d0

(4.9)

It is clear that Rz < 1(≧ 1) for p < p∗1 (≧ p∗1 ) and Rd < 1(≧ 1) for p > p∗2 (≦ p∗2 ). When p∗1 > p∗2 , the disease-free equilibrium E0 is stable if p∗2 < p < p∗1 . It means that both strains can die out. When p∗1 < p∗2 , E0 is unstable for any dengue vaccine rate p and both strains can never die out by a vaccination campaign against dengue. On the other hand, Zika (or dengue) can be eradicated by finding an ADE factor φv that would make the boundary equilibrium Ed (or Ez ) with the dengue (or Zika) strain surviving stable from Theorem 4.6 (or Theorem 4.5). Consider Rzd = Rzd (φv ) and Rdz = Rdz (φv ) as functions of φv , where Rdz and Rzd are defied in (3.14) and (3.20), respectively. Solving for φv from Rzd (φv ) = 1 gives φv = φv 1 , where

φv1 = 1 + in which rz4 =

(1 − (Rz1 + Rz2 + Rz3 + rz4 )) rz4

∗∗ V ∗∗ a2 βhmz βmhz µ2h Sm h

Λ2h µm (γz +µh )

,

(4.10)

, Rz1 , Rz2 and Rz3 are defined in (3.21). It is clear that Rzd < 1(≧ 1) for φv < φv 1 (≧ φv 1 ). It

is too complicated to analysis the curve Rdz (φv ) = 1 analytically. With the help of simulation, we can present that φv = φv 2 satisfies Rdz (φv 2 ) = 1 numerically in Fig. 3, and Rdz < 1(≧ 1) for φv < φv 2 (≧ φv 2 ). If p < p∗1 , then boundary equilibrium Ed is stable if and only if 1 < φv < φv 1 . It means that Zika can die out and dengue is endemic. If p > p∗2 , then boundary equilibrium Ez is stable if and only if 1 < φv < φv 2 . This implies that dengue can die out and Zika is endemic. Furthermore, if p∗1 < p < p∗2 , both Zika and dengue are present if φv > max{φv 1 , φv 2 }. Case (ii) Rz0 > 1, Rd0 < 1. Since Rz > Rz0 , Rz is always more than 1 for any rate of vaccination against dengue. Moreover, dengue can be eradicated and Zika is endemic if 1 < φv < φv 2 . Case (iii) Rz0 > 1, Rd0 > 1. E0 is unstable for any rate of vaccination against dengue p and both strains can never die out by a vaccination campaign against dengue. However, Rd can less than 1 through a vaccination campaign against dengue p > p∗2 , and Ez is stable if ADE factor φv satisfies 1 < φv < φv 2 . Since vaccination against dengue will increase the outbreak of Zika and reduce the outbreak of dengue, Rz0 < 1 and Rd0 > 1 is very interesting in applications. We give stability regions of system (2.1) in (p, φv ) plane based on the analysis of Case (i). We know that there is a unique ADE factor φv 3 such that Rz (φv 3 ) = 1 given by

φv 3 = 1 +

1 p

(

1 (Rz0 )2

) −1 .

(4.11)

It is clear that Rz < 1(≧ 1) for φv < φv 3 (≧ φv 3 ). For convenience analysis later, use φv = φv 3 instead of p = p∗1 to express Rz = 1. To confirm our conjecture, we fix all parameters in Table 2. We plot the parameters regions in (p, φv ) plane which is shown in Fig. 3. Cures L1 : φv (p) = φv 1 (p), L2 : φv (p) = φv 2 (p), L3 : φv (p) = φv 3 (p), and vertical straight line L4 : p = p∗2 divide the (p, φv ) plane into six sub-regions Ω1 , Ω2 , Ω3 , Ω4 , Ω5 and Ω6 . From Theorems 4.1 and 4.5–4.7, we conclude that there is one stable boundary equilibrium Ed in region

Ω1 = {(p, φv ) | φv < φv1 (p), 0 ≤ p < p∗2 }, there is one stable disease-free equilibrium E0 in region

Ω2 = {(p, φv ) | φv < φv3 (p), p∗2 < p ≤ 1},

264

L. Wang and H. Zhao / Physica A 522 (2019) 248–273 Table 2 Parameter descriptions. Parameter Nh Nm

µh µm Λh Λm a

βmhz βhmz βmhd βhmd γz γd γdz φ1 , φ2 , η1 , η2 ϵ1 , ϵ2

Range

Value

Reference

(1, 10) · Nh

2.05 × 108 6.5 × 108

[48] [49] [50] [51] Calculation Calculation [51] [52] [52] [53] [53] [54]

(

7 50

,

7 4

)

(2.1, 7) (0, 0.97) (0, 0.97) (0, 0.97) (0 ( , 0.97) ) 7

( 30 7 30

, ,

7 5 7 6

)

(1, 5) (0, 1)

1 75×52 7 6

1 2.05 × 108 × 75× 52 6.5 × 108 × 76 2.1 0.25 0.3 0.32 0.305 7 6 7 7 7 10

1.5 0.8

[55] Assume [38] Assume

Fig. 3. Stability regions of system (2.1) in (p, φv ) plane. In region Ω1 , Rd > 1, Rzd < 1. In region Ω2 , Rz < 1, Rd < 1. In region Ω3 , Rz > 1, Rdz < 1. In region Ω4 , Rd < 1, Rdz > 1. In region Ω5 , Rz > 1, Rd > 1, Rzd > 1, Rdz > 1. In region Ω6 , Rz < 1, Rzd > 1. Here Rz0 = 0.8777 < 1 and Rd0 = 1.0814 > 1 with all other parameters values shown in Table 2.

there is one stable boundary equilibrium Ez in region

Ω3 = {(p, φv ) | φv3 (p) < φv < φv2 (p), 0 ≤ p ≤ 1}, there is one unstable boundary equilibrium Ez in region

Ω4 = {(p, φv ) | φv > φv2 (p), p∗2 < p ≤ 1}, the system (2.1) is persistence where both Zika and dengue are present in region

Ω5 = {(p, φv ) | φv > max{φv2 (p), φv2 (p), φv3 (p)}, 0 ≤ p < p∗2 } and there is one unstable boundary equilibrium Ed in region

Ω6 = {(p, φv ) | φv1 (p) < φv < φv3 (p), 0 ≤ p < p∗2 }. Fig. 3 can be very helpful for examining the joint effect of vaccination against dengue and ADE factor. If (p, φv ) belongs to region Ω2 , then Zika and dengue cannot spread. Only dengue will break out if (p, φv ) belongs to region Ω1 , and only Zika will be present if (p, φv ) belongs to region Ω3 . Both Zika and dengue will be present when (p, φv ) belongs to region Ω5 . More detailed simulation results for the dynamic in various regions shown in Fig. 3 are discussed in the next section.

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

265

5. Numerical simulations In this section, we implement numerical simulations to confirm and extend the analytic results through the change of two important factors: dengue vaccine rate (p) and ADE factor (φv ), and illustrate the effect of basic and invasion reproduction numbers on the disease dynamics. We will examine the effect of dengue vaccine on the prevalence of disease under a given level of ADE factor. Based on the chosen parameters values in Table 2, we consider the case 1 < φv < max{φv 1 (p), φv 2 (p)} displayed in Fig. 3. Choosing φv = 2.6, we can calculate the critical values of dengue vaccine rate p∗1 = 0.18635 and p∗2 = 0.14491. Firstly, if 0 ≤ p < p∗2 , it follows from Fig. 3 that (p, φv ) = (p, 2.6) lines in region Ω1 in which Ed is stable. Fig. 4 illustrates this pattern where the dengue vaccine rate is chosen as p = 0.13. It means that only dengue is present and Zika cannot spread. As p continuous to increase, if p∗2 < p < p∗1 , (p, φv ) will enter region Ω2 , then the disease-free equilibrium is a unique stable state. For this pattern, we choose p = 0.15 to obtain Fig. 5 in which both dengue and Zika will be die out. When p continuous to increase in excess of p∗1 , (p, φv ) will enter region Ω3 in which Ez is stable. For this pattern, we choose p = 0.2. As depicted in Fig. 6, Zika will break out and dengue will be die out. Figs. 4–6 illustrate vaccination against dengue can induce break out of Zika and decrease spread of dengue. Next, we will present some results regarding the effect of ADE factor under a given level of dengue vaccine rate p = 0.13. We change the value of φv and keep the other parameters unchanged. Consider two different values of ADE factor: φv = 2.6, 3.5. For these simulations, parameter values were chosen such that Rzd = 0.9676 < 1 (i.e., φv = 2.6) and Rzd = 1.0131 > 1 (i.e., φv = 3.5). Figs. 4 and 7 display the two types of dynamical behaviors of system (2.1) correspondingly, and demonstrate the role of the invasion reproduction number for Zika Rzd . Our analytic results suggest that, for Rzd < 1, the boundary equilibrium Ed is stable. This is indeed confirmed by our simulations shown in Fig. 4 which show that dengue converges to a positive level while Zika tends to extinction. The simulations shown in Fig. 7 extend our analytical result presented in Theorem 4.7. While we do not have results on the existence and stability of interior equilibrium of the system (2.1), our numerical studies indicate that stable interior equilibrium may exist exhibited in Fig. 7. It follows from Figs. 4 and 7 that, when dengue is persist, Zika can invade in the population if invasion reproduction number for Zika Rzd is greater than one, while Zika will die out if invasion reproduction number for Zika Rzd is less than one. In Fig. 4, we observe that, when both the basic and invasion reproduction numbers for Zika are less than one, Zika cannot survive in the population. However, when the basic reproduction number is less than one, but the invasion reproduction

Fig. 4. The stable boundary equilibrium Ed of system (2.1) with p = 0.13, φv = 2.6, (p, φv ) ∈ Ω1 , Rd > 1 and Rzd < 1.

Fig. 5. The stable disease-free equilibrium E0 of system (2.1) with p = 0.15, φv = 2.6, (p, φv ) ∈ Ω2 , Rz < 1 and Rd < 1.

266

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

Fig. 6. The stable boundary equilibrium Ez of system (2.1) with p = 0.2, φv = 2.6, (p, φv ) ∈ Ω3 , Rz > 1 and Rdz < 1.

Fig. 7. The existence of interior equilibrium and persistence of system (2.1) with p = 0.13, φv = 3.5, (p, φv ) ∈ Ω5 , Rz > 1, Rd > 1 Rzd > 1 and Rdz > 1.

number is greater than one, it follows from Fig. 8(A) that the survival of Zika depends on whether dengue is already present in the population or not. In this case, when dengue is present, Zika can invade in the population, but when dengue is absent, then Zika dies out. The same results hold for the dengue, as seen in Fig. 8(B).

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

267

Fig. 8. The comparison of full system (2.1) and (A) sub-system (4.3) and (B) sub-system (4.4). (A) Choose p = 0.118, φv = 3.5, (p, φv ) ∈ Ω6 . (B) Choose p = 0.15, φv = 3.5, (p, φv ) ∈ Ω4 .

6. A numerical application to Brazil 6.1. Data fitting There is no systematic data of human infection available to allow comprehensive system studies of the co-infection of Zika and dengue in Brazil. Luckily, time series of weekly Zika and dengue accumulated cases are collected from the Brazil Ministry of Health [56] and the World Health Organization [57]. We choose the seventh week in 2016 as the first observation date and use the weekly reported accumulated cases of Zika and dengue in Brazil from February 19, 2016 to April 14, 2018 to carry out our study. We apply the system (2.1) to fit the accumulated cases in Brazil from February 19, 2016 to April 14, 2018 and take a week as the unit time. So the considered period is about 113 weeks. In this period, the Brazilian population is not vaccinated. So, p = 0. Some other parameters can be get from previous literatures and websites about Zika and dengue for Brazil given in Table 3, while other parameters will be estimated by applying the system (2.1) to simulate the real data by using Markovchain Monte Carlo (MCMC) simulations. Denote unknown parameters as Θ = (Nm , γdz , ϵ1 , ϵ2 ). The relationship between the sizes of reported accumulated Zika and dengue cases and that of system (2.1) accumulated Zika and dengue cases can be written as y1 (t) = y¯1 (t) + ϵ and y2 (t) = y¯2 (t) + ϵ , respectively, where y1 (t) and y2 (t) are theoretic accumulated Zika and

Table 3 Parameter descriptions. Parameter Nh Nm

µh µm Λh Λm a

βmhz βhmz βmhd βhmd γz γd

Range

Value

Reference

(1, 10) · Nh

2.05 × 108 5.2275 × 108

[48] Estimated [50] [51] Calculation Calculation [51] [52] [52] [53] [53] [54]

(

7 50

,

7 4

(2.1, 7) (0, 0.97) (0, 0.97) (0, 0.97) ((0, 0.97) ) 7

( 30 7 30

, ,

7 5 7 6

γdz φ1 , φ2 , η1 , η2 ϵ1 , ϵ2

)

(1, 5) (0, 1)

)

1 75×52

0.3 1 2.05 × 108 × 75× 52 5.2275 × 108 × 0.3 2.1 0.1 0.3 0.205 0.21 7 5 7 8 7 15

1.5 0.8

[55] Estimated [38] Estimated

268

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

Fig. 9. Fitting curve of real data from February 19, 2016 to April 14, 2018 by deterministic system (2.1), where the star denotes the real data of accumulated (A) Zika and (B) dengue cases. Baseline parameters from Table 3 are used with initial conditions Sh (0) = 9.8 × 107 , Iz (0) = 1.6 × 104 , Id (0) = 3.1 × 104 , Idz (0) = 100, Rz (0) = 6.0 × 104 , Rd (0) = 7.5 × 106 , Ydz (0) = 2.0 × 103 , Yzd (0) = 200, Vh (0) = 0, Iv z (0) = 0, Imz (0) = 5.0 × 104 , Imd (0) = 6.0 × 105 .

dengue cases and their change with time are described as dy1 (t) dt dy2 (t) dt

= aβmhz Imz = aβmhd Imd

Sh Nh Sh Nh

+ η2 aβmhz Imz + η1 aβmhd Imd

Id Nh Iz Nh

+ φ2 aβmhz Imz + φ1 aβmhd Imd

Rd Nh Rz Nh

+ φv aβmhz Imz

Vh Nh

,

,

respectively, and y¯1 (t) and y¯2 (t) are surveillance Zika and dengue data in week t. ϵ follows a Gaussian distribution. Using the same method and step in Ref. [49], the vector of parameters Θ = (Nm , γdz , ϵ1 , ϵ2 ) are estimated in Table 3. Through some rational assumptions and parameter estimations in Table 3, numerical simulations of the accumulated human Zika and dengue cases using system (2.1) are shown in Fig. 9. We observe that the model fits well the data over the same time period with the parameters in Table 3. With these parameters, we have that the Zika basic reproduction number is Rz = 0.8962, and the dengue basic reproduction number is Rd = 1.3578. Further, substituting the parameter values shown in Table 3 into the expression (3.20) of Rzd , we calculate invasion reproduction number for Zika Rzd = 0.9933 < 1. From Theorem 4.6, the boundary equilibrium Ed of system (2.1) is stable. This illustrates that dengue epidemic is endemic while Zika epidemic will eventually disappear in Brazil. 6.2. The effect of parameters p and φv on the accumulated cases of Zika and dengue in Brazil To examine the effect of dengue vaccine on the prevalence of Zika and dengue in Brazil, we consider different values of dengue vaccine rates: p = 0, 0.2, 0.4, 0.6, 0.8, 1.0, under a given level of ADE factor φv = 1.35. Fig. 10(A) and (B) show the effect of dengue vaccine on the number of accumulated dengue and Zika, respectively. If 20% of the newborns is vaccinated, the number of accumulated dengue in Brazil will be reduced to 813100 cases which is 1030347 less than that of real reported data at fifteenth week of 2018, while the number of accumulated Zika in Brazil will be increased to 268200 cases which is 33195 more than that of real reported data. Fig. 10(A) illustrates that vaccination against dengue can reduce the number of the accumulated dengue infections, while Fig. 10(B) shows that the number of accumulated Zika infections will increase significantly under a higher rate of vaccination against dengue. Since ADE factor φv > 1 increases transmission probability from Zika infection mosquitoes to vaccinated humans against dengue, we only give the effect of ADE factor on the number of accumulated Zika under a given level of rate of vaccination against dengue p = 0.2. From Fig. 10(C), if φv = 1.4, that is, transmission probability from Zika infection mosquitoes to vaccinated humans against dengue is 1.4 times higher than before, the number of accumulated Zika in Brazil will reach 273900 cases at fifteenth week of 2018. If transmission probability from Zika infection mosquitoes to vaccinated humans against dengue is three times higher than before, the number of accumulated Zika in Brazil will reach 770300 cases. Fig. 10(C) illustrates that the scale of Zika outbreaks in Brazil will get larger when the φv increases.

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

269

Fig. 10. The accumulated cases in Brazil with respect to p and φv . (A) and (B) φv = 1.35, (C) p = 0.2.

Fig. 11. A graph of the regions of effective vaccination rates against dengue p versus the ADE factor φv .

From Fig. 10(C), we obtain that protecting oneself for vaccinated humans against dengue from infection mosquitoes are an effective methods to decrease the size of accumulated Zika case in Brazil. However, Fig. 10(A) and (B) illustrate that higher dengue vaccine rate will increase Zika cases while dengue cases greatly reduce, and lower dengue vaccine rate leads to increasing dengue cases while Zika cases reduce. To this end, we need give a criterion of vaccination against dengue strategy. Fig. 11 suggests that there is a ‘‘low bound’’ for vaccination against dengue, that is, p = p∗2 = 0.4576, above which the dengue can be eradicated when φv < φv∗ = 1.521. If (φv , p) belongs to region D∗ shown in Fig. 11, then Zika and dengue will die out. Fig. 11 provides an ‘‘upper bound’’ for vaccination against dengue , that is, p = p∗1 (φv ), below which the Zika can be eradicated. Here, ADE factor plays an important role in the effects of vaccination against dengue strategies. Thus, when designing vaccination against dengue strategies, one should take into consideration ADE factor. 7. Conclusion and discussion In this paper, we consider a Zika–dengue co-infection model that focuses on investigating the effect of ADE and vaccination programs against dengue on the control and prevention of infectious diseases. One of the main goals of this work is to investigate how the incorporation of both ADE and dengue vaccine influences disease dynamics. Another goal is to determine the role of vaccination against dengue on Zika prevalence. We first calculate and interpret the basic reproduction number, and the invasion reproduction number for the two strains when only one is already present in the population and the other is trying to invade. Then a systematic qualitative analysis including stability of equilibria for the system is obtained. We obtain that although vaccination against dengue has a positive effect on the control of dengue, it has a negative effect on the control of Zika. Hence, the threshold conditions for disease elimination with vaccination against dengue strategies are

270

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

Fig. 12. The graphs of the (A) basic and (B) invasion reproduction numbers when p varies. Here φi = 1, i = 1, 2, v . All other parameters values are shown in Table 2.

derived, that is, p∗2 < p < p∗1 . Otherwise, both Zika and dengue can never die in a vaccination campaign against dengue. On the other hand, one of the two strains can be eradicated by finding an ADE factor φv 1 or φv 2 such that φv < φv 1 or φv < φv 2 . φv can be reduced by protecting oneself for vaccinated humans against dengue from the bite of infection mosquitoes. It suggests that to identify the best disease control vaccination strategies, policymakers must take into account the factor of ADE. This may have significant implications for public health. Our analysis and numerical simulations suggest that elimination and outbreak of disease are generated by changing the factor of ADE under a given level of dengue vaccine rate. In our model, if we do not consider the factor of ADE, then φi = 1, i = 1, 2, v . That is, secondary infections with strain Zika (or dengue) take place as if they were primary infections. In this case, Rz = Rz0 . Rz is not affected by the dengue vaccine. In the numerical simulations (see Fig. 12(B)), Rzd is almost unaffected by the dengue vaccine. From Fig. 12, one can obtain that vaccination against dengue has no effect on the control of Zika but has a positive effect on the control of dengue. Therefore, antibody-dependent enhancement plays an important role in promoting an outbreak of Zika when dengue vaccine is introduced in humans. As an application, we use our system to simulate the reported human Zika and dengue cases in Brazil from February 19, 2016 to April 14, 2018 and obtain reasonable matches. We estimate the Zika basic reproduction number Rz = 0.8962, the dengue basic reproduction number Rd = 1.3578, and invasion reproduction number for Zika Rzd = 0.9933 < 1. From Theorem 4.6, the dengue epidemic is endemic while the Zika epidemic will eventually disappear in Brazil. Vaccination is a cost-effective measure to control dengue, but vaccination against dengue will aggravate the Zika outbreak. To this end, we give a vaccination against dengue strategy. Fig. 11 suggests that there is a ‘‘low bound" for vaccination against dengue p = 0.4576, and an ‘‘upper bound’’ for vaccination against dengue p = p∗1 (φv ) under the level of φv < 1.521. The dengue vaccinate rate is within this range. Zika and dengue will die out. A vaccination and treatment program is seen as an important measure used in infectious disease control. In future work, we intend to study an optimal strategy, considering both the treatment costs of infected individuals and the vaccination costs. By using optimal control theory, we will investigate the cost effectiveness of the control strategies using incremental cost-effectiveness ratio.

Acknowledgments The authors highly appreciate the anonymous reviewers for providing valuable suggestions which helped us to improve the paper. The work is partially supported by National Natural Science Foundation of China under Grant 11571170. Funding of Jiangsu Innovation Program for Graduate Education KYZZ16_0162.

Appendix. Proof of Lemma 2.1 Let x(t) = (Sh , Iz , Id , Idz , Rz , Rd , Ydz , Yzd , Vh , Iv z , Rdz , Sm , Imz , Imd ). It follows from system (2.1), x(0) ≥ 0, Nh (t) = and Nm (t) =

Λm µm

for all t ≥ 0 that

Λh µh

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

271

[ ∫ t( ) ] aβmhz Imz (s) aβmhd Imd (s) Sh (t) = Sh (0)exp − + ds − µh t Nh Nh 0 [ ∫ t( ) ] ∫ t aβmhz Imz (ξ ) aβmhd Imd (ξ ) + (1 − p)Λh exp − + dξ − µh (t − s) ds ≥ 0, Nh Nh 0 s [ ] ∫ t aβmhd Iz (t) = Iz (0)exp − Imd (s)ds − (γz + µh )t Nh 0 [ ] ∫ t ∫ t aβmhd aβmhz Imz (s)Sh (s)exp − Imd (ξ )dξ − (γz + µh )(t − s) ds ≥ 0, + Nh Nh s [0 ] ∫ t aβmhz Id (t) = Id (0)exp − Imz (s)ds − (γd + µh )t Nh 0 [ ] ∫ t ∫ t aβmhd aβmhd + Imd (s)Sh (s)exp − Imd (ξ )dξ − (γd + µh )(t − s) ds ≥ 0, Nh Nh 0 s ) ∫ t( a β I (s)I (s) aβmhz Imz (s)Id (s) mhd md z −(γdz +µh )t + Idz (t) = Idz (0)e + e−(γdz +µh )(t −s) ds ≥ 0, Nh Nh 0 ] [ ∫ φ1 aβmhd t Imd (s)ds − µh t Rz (t) = Rz (0)exp − Nh 0 [ ] ∫ t ∫ φ1 aβmhd t Iz (s)exp − Imd (ξ )ξ − µh (t − s) ds ≥ 0, + γz Nh 0 s [ ] ∫ φ2 aβmhz t Rd (t) = Rd (0)exp − Imz (s)ds − µh t Nh 0 [ ] ∫ t ∫ φ2 aβmhz t Id (s)exp − + γd Imz (ξ )ξ − µh (t − s) ds ≥ 0, Nh 0 s ∫ t φ a β 2 mhz Imz (s)Rd (s)e−(γz +µh )(t −s) ds ≥ 0, Ydz (t) = Ydz (0)e−(γz +µh )t + Nh 0 ∫ φ1 aβmhd t Yzd (t) = Yzd (0)e−(γd +µh )t + Imd (s)Rz (s)e−(γd +µh )(t −s) ds ≥ 0, Nh 0 ] [ ∫ φv aβmhz t Imz (s)ds − µh t Vh (t) = Vh (0)exp − Nh 0 ] [ ∫ t ∫ φv aβmhz t exp − + pΛh Imz (ξ )dξ − µh (t − s) ds ≥ 0, Nh 0 s ∫ t φ a β v mhz Iv z (t) = Iv z (0)e−(γz +µh )t + Imz (s)Vh (s)e−(γz +µh )(t −s) ds ≥ 0, Nh 0 ∫ t Rdz (t) = Rdz (0)e−µh t + (γdz Idz (s) + γz Ydz (s) + γd Yzd (s) + γz Iv z (s))e−µh (t −s) ds ≥ 0, 0 t

∫

−aβhmz

Sm (t) = Sm (0)exp 0

(Iz (s) + Idz (s) + Ydz (s) + Ivz (s)) Nh

ds ×

] (Id (s) + Idz (s) + Yzd (s)) ds − µh t −aβhmd Nh 0 ( ) ∫ t ∫ t Iz (ξ ) + Idz (ξ ) + Ydz (ξ ) + Iv z (ξ ) + Λm exp −aβhmz dξ × Nh 0 s [∫ t ) ] ( Id (ξ ) + Idz (ξ ) + Yzd (ξ ) exp −aβhmd dξ − µh (t − s) ds ≥ 0, Nh s ( ) ∫ t I (s) + Idz (s) + Ydz (s) + Iv z (s) z −µm t Imz (t) = Imz (0)e + aβhmz Sm (s)e−µm (t −s) ds ≥ 0, Nh 0 ( ) ∫ t Id (s) + Idz (s) + Yzd (s) Imd (t) = Imd (0)e−µm t + aβhmd Sm (s)e−µm (t −s) ds ≥ 0. [∫

t

exp

Nh 0 If Sh (0) > 0, Iz (0) > 0, Id (0) > 0, Idz (0) > 0, Rz (0) > 0, Rd (0) > 0, Ydz (0) > 0, Yzd (0) > 0, Vh (0) > 0, Iv z (0) > 0, Rdz (0) > 0, Sm (0) > 0, Imz (0) > 0, Imd (0) > 0, then Sh (t) > 0, Iz (t) > 0, Id (t) > 0, Idz (t) > 0, Rz (t) > 0, Rd (t) > 0, Ydz (t) > 0, Yzd (t) > 0, Vh (t) > 0, Iv z (t) > 0, Rdz (t) > 0, Sm (t) > 0, Imz (t) > 0, Imd (t) > 0 for ∀t > 0.

272

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

Furthermore, from system (2.1), Nh (t) = satisfy the following equations respectively, dNh (t) dt

= Λh − µh Nh = 0,

dNm (t) dt

Λh µh

and Nm (t) =

Λm µm

for all t ≥ 0, the humans and mosquitoes population size

= Λm − µm Nm = 0.

So, Sh + Iz + Id + Idz + Rz + Rd + Ydz + Yzd + Vh + Iv z + Rdz = positively invariant set.

Λh , Sm µh

+ Imz + Imd =

Λm µm

for ∀ t ≥ 0, which implies that Γ is a

References [1] L. Esteva, C. Vargas, A model for dengue disease with variable human population, J. Math. Biol. 38 (3) (1999) 220–240. [2] C.B. Marcondes, M.F.F. Ximenes, Zika virus in Brazil and the danger of infestation by Aedes(Stegomyia) mosquitoes, Rev. Soc. Bras. Med. Trop. 49 (2016) 4–10. [3] G. Chowell, et al., Estimation of the reproduction number of dengue fever from spatial epidemic data, Math. Biosci. 208 (2) (2007) 571–589. [4] L. Esteva, C. Vargas, Coexistence of different serotypes of dengue virus, J. Math. Biol. 46 (1) (2003) 31–47. [5] S. Bhatt, et al., The global distribution and burden of dengue, Nature 496 (7446) (2013) 504. [6] G.W.A. Dick, S.F. Kitchen, A.J. Haddow, Zika virus (I). Isolations and serological specificity, Trans. R. Soc. Trop. Med. Hyg. 46 (5) (1952) 509–520. [7] M.R. Duffy, et al., Zika virus outbreak on Yap Island, Federated States of Micronesia, New England J. Med. 360 (2009) 2536–2543. [8] V.M. Cao-Lormeau, et al., Zika virus, French polynesia, South pacific, 2013, Emerg. Infect. Dis. 20 (6) (2014) 1085. [9] G.S. Campos, A.C. Bandeira, S.I. Sardi, Zika virus outbreak, Bahia, Brazil, Emerg. Infect. Dis. 21 (2015) 1885–1886. [10] M. Roa, Zika virus outbreak: reproductive health and rights in Latin America, Lancet 387 (2016) 843-843. [11] World Health Organization (WHO), WHO statement on the first meeting of the International Health Regulations (2005) Emergency Committee on Zika virus and observed increase in neurological disorders and neonatal malformations, February 1, 2016. http://www.who.int/mediacentre/news/ statements/2016/1st-emergency-committee-zika/en/. (Accessed 26 February 2016). [12] M. Dupont-Rouzeyrol, et al., Co-infection with Zika and dengue viruses in 2 patients, New Caledonia, 2014, Emerg. Infect. Dis. 21 (2) (2015) 381. [13] R. Pessôa, et al., Investigation into an outbreak of dengue-like illness in Pernambuco, Brazil, revealed a cocirculation of Zika, Chikungunya, and dengue virus type 1, Medicine 95 (12) (2016). [14] W.E. Villamil-Gómez, et al., Dengue, chikungunya and Zika co-infection in a patient from Colombia, J. Infect. Publ. Health 9 (5) (2016) 684–686. [15] W.E. Villamil-Gómez, et al., Zika, dengue, and chikungunya co-infection in a pregnant woman from Colombia, Int. J. Infect. Dis. 51 (2016) 135–138. [16] L. Priyamvada, et al., Human antibody responses after dengue virus infection are highly cross-reactive to Zika virus, Proc. Nat. Acad. Sci. 113 (28) (2016) 7852–7857. [17] W. Dejnirattisai, et al., Dengue virus sero-cross-reactivity drives antibody-dependent enhancement of infection with zika virus, Nat. Immunol. 17 (9) (2016) 1102. [18] L.M. Paul, et al., Dengue virus antibodies enhance Zika virus infection, Clin. Transl. Immunol. 5 (12) (2016). [19] P. Castanha, et al., Dengue virus-specific antibodies enhance Brazilian Zika virus infection, J. Infect. Dis. 215 (5) (2017) 781–785. [20] B. Londono-Renteria, et al., A relevant in vitro human model for the study of Zika virus antibody-dependent enhancement, J. Gen. Virol. 98 (7) (2017) 1702–1712. [21] C.J. Reynolds, et al., T cell immunity to Zika virus targets immunocompetent epitomes that show cross-reactivity with other flaviviruses, Sci. Rep. 8 (1) (2018) 672. [22] K.S. Vannice, A. Durbin, J. Hombach, Status of vaccine research and development of vaccines for dengue, Vaccine. 34 (2016) 934–2938. [23] H. Saba, J.G.V. Miranda, M.A. Moret, Self-organized critical phenomenon as a q q math Container Loading Mathjax-exponential decay-Avalanche epidemiology of dengue, Physica A 413 (11) (2014) 205–211. [24] D. He, D. Gao, Y. Lou, S. Zhao, S. Ruan, A comparison study of Zika virus outbreaks in French Polynesia, Colombia and the State of Bahia in Brazil, Sci. Rep. 7 (1) (2017) 273. [25] D.J. Bicout, K. Chalvet-Monfray, P. Sabatier, Infection persistence time of Aedes breeding habitats, Physica A 305 (3) (2002) 597–603. [26] D. Gao, et al., Prevention and control of Zika fever as a mosquito-borne and sexually transmitted disease, Sci. Rep. 6 (2016). [27] F.M.M. Pereira, P.H.T. Schimit, Dengue fever spreading based on probabilistic cellular automata with two lattices, Physica A 499 (2018). [28] D.H. Barmak, C.O. Dorso, M. Otero, Modelling dengue epidemic spreading with human mobility, Physica A 447 (2016) 129–140. [29] Z. Feng, J.X. Velasco-Hernndez, Competitive exclusion in a vector-host model for the dengue fever, J. Math. Biol. 35 (5) (1997) 523–544. [30] M. Sriprom, P. Barbazan, I.M. Tang, Destabilizing effect of the host immune status on the sequential transmission dynamic of the dengue virus infection, Math. Comput Model. 45 (9–10) (2007) 1053–1066. [31] L. Esteva, C. Vargas, Coexistence of different serotypes of dengue virus, J. Math. Biol. 46 (1) (2003) 31–47. [32] N. Nuraini, H. Tasman, E. Soewono, K.A. Sidarto, A with-in host dengue infection model with immune response, Math. Comput. Model. 49 (5–6) (2009) 1148–1155. [33] A. Mishra, S. Gakkhar, The effects of awareness and vector control on two strains dengue dynamics, Appl. Math. Comput. 246 (2014) 159–167. [34] T.T. Zheng, L.F. Nie, Modelling the transmission dynamics of two-strain dengue in the presence awareness and vector control, J. Theoret. Biol. 443 (2018) 82–91. [35] D.M. Morens, S.B. Halstead, Measurement of antibody-dependent infection enhancement of four dengue virus serotypes by monoclonal and polyclonal antibodies, J. Gen. Virol. 71 (12) (1990) 2909–2914. [36] L.C. Katzelnick, et al., Antibody-dependent enhancement of severe dengue disease in humans, Science 358 (6365) (2017) 929–932. [37] N. Ferguson, R. Anderson, S. Gupta, The effect of antibody-dependent enhancement on the transmission dynamics and persistence of multiple-strain pathogens, Proc. Nat. Acad. Sci. 96 (2) (1999) 790–794. [38] D.A. Cummings, I.B. Schwartz, L. Billings, L.B. Shaw, D.S. Burke, Dynamic effects of antibody-dependent enhancement on the fitness of viruses, Proc. Natl. Acad. Sci. USA 102 (42) (2005) 15259–15264. [39] L. Billings, A. Fiorillo, I.B. Schwartz, Vaccinations in disease models with antibody-dependent enhancement, Math. Biosci. 211 (2) (2008) 265–281. [40] K. Hu, C. Thoens, et al., The effect of antibody-dependent enhancement, cross immunity, and vector population on the dynamics of dengue fever, J. Theoret. Biol. 319 (2013) 62–74. [41] B. Tang, Y. Xiao, J. Wu, Implication of vaccination against dengue for Zika outbreak, Sci. Rep. 6 (2016) 35623. [42] Z. Feng, Z. Qiu, Z. Sang, C. Lorenzo, J. Glasser, Modeling the synergy between HSV-2 and HIV and potential impact of HSV-2 therapy, Math. Biosci. 245 (2) (2013) 171–187. [43] J. Cohen, Q & A with Scott Halstead: Zika will subside in 5 years, max, Science Mar. 7 (2016) 7, http://www.sciencemag.org/news/2016/03/qa-scotthalstead-zika-willsubside-5-years-max. [44] N. Chitnis, J.M. Cushing, J.M. Hyman, Bifurcation analysis of a mathematical model for malaria transmission, SIAM J. Appl. Math. 67 (1) (2006) 24–45.

L. Wang and H. Zhao / Physica A 522 (2019) 248–273

273

[45] O. Diekmann, J.A.P. Heesterbeek, J.A. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, J. Math. Biol. 28 (4) (1990) 365–382. [46] P. Van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (1) (2002) 29–48. [47] Z. Shuai, P. Van den Driessche, Global stability of infectious disease models using Lyapunov functions, SIAM J. Appl. Math. 73 (4) (2013) 1513–1532. [48] City Populations Worldwide, Brazil population. http://population.city/brazil/. (Accessed 10 April 2018). [49] L. Wang, H. Zhao, S.M. Oliva, H. Zhu, Modeling the transmission and control of Zika in Brazil, Sci. Rep. 7 (1) (2017) 7721. [50] A. Wiratsudakul, P. Suparit, C. Modchang, Dynamics of Zika virus outbreaks: an overview of mathematical modeling approaches, PeerJ 63 (2018) e4526. [51] M. Andraud, N. Hens, C. Marais, P. Beutels, Dynamic epidemiological models for dengue transmission: a systematic review of structural approaches, PLoS One 7 (2011) 3161–3164. [52] Q. Zhang, et al., Spread of Zika virus in the Americas, Proc. Nat. Acad. Sci. 114 (22) (2017) E4334–E4343. [53] L. Lambrechts, et al., Impact of daily temperature fluctuations on dengue virus transmission by Aedesaegypti, Proc. Nat. Acad. Sci. 108 (18) (2011) 7460–7465. [54] A.C. Gourinat, O. O’Connor, E. Calvez, C. Goarant, M. Dupont-Rouzeyrol, Detection of Zika virus in urine, Emerg. Infect. Dis. 21 (1) (2015) 84. [55] S.B. Halstead, Dengue, Lancet 370 (9599) (2007) 1644–1652. [56] Brazil Ministry of Health, Zika cases from the Brazil Ministry of Health. http://portalms.saude.gov.br/boletins-epidemiologicos. (Accessed 15 April 2018). [57] World Health Organization, Case definitions of Zika virus. http://apps.who.int/iris/bitstream/handle/10665/204381/WHO_ZIKV_SUR_161_por.pdf; jsessionid=D4FD38CA0C8CD27D00310905DC84D644?sequence=5. (Accessed 12 Februry 2016).