Birth-pulse models of Wolbachia-induced cytoplasmic incompatibility in mosquitoes for dengue virus control

Nonlinear Analysis: Real World Applications 22 (2015) 236–258

Contents lists available at ScienceDirect

Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa

Birth-pulse models of Wolbachia-induced cytoplasmic incompatibility in mosquitoes for dengue virus control Xianghong Zhang a , Sanyi Tang a,∗ , Robert A. Cheke b a

College of Mathematics and Information Science, Shaanxi Normal University, Xi’an, 710062, PR China

b

Natural Resources Institute, University of Greenwich at Medway, Central Avenue, Chatham Maritime, Chatham, Kent, ME4 4TB, UK

highlights • • • •

Birth-pulse models concerning the effects of CI on dengue control are proposed. Two different density dependent death functions are proposed and considered. Strategies of eradicating mosquito and population replacement have been addressed. The two strategies can be realized for different death rate functions.

article

info

Article history: Received 20 June 2014 Received in revised form 5 September 2014 Accepted 12 September 2014 Available online 14 October 2014 Keywords: Dengue fever Wolbachia Cytoplasmic incompatibility Birth-pulse Strong and weak density dependent death rates

abstract Dengue fever, which affects more than 50 million people a year, is the most important arboviral disease in tropical countries. Mosquitoes are the principal vectors of the dengue virus but some endosymbiotic Wolbachia bacteria can stop the mosquitoes from reproducing and so interrupt virus transmission. A birth-pulse model of the spread of Wolbachia through a population of mosquitoes, incorporating the effects of cytoplasmic incompatibility (CI) and different density dependent death rate functions, is proposed. Strategies for either eradicating mosquitoes or using population replacement by substituting uninfected mosquitoes with infected ones for dengue virus prevention were modeled. A model with a strong density dependent death function shows that population replacement can be realized if the initial ratio of number of infected to the total number of mosquitoes exceeds a critical value, especially when transmission from mother to offspring is perfect. However, with a weak density dependent death function, population eradication becomes difficult as the system’s solutions are sensitive to initial values. Using numerical methods, it was shown that population eradication may be achieved regardless of the infection ratio only when parameters lie in particular regions and the initial density of uninfected is low enough. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction The vector-borne diseases, dengue fever and dengue hemorrhagic fever, are emerging globally as the most important arboviral diseases currently threatening human populations. More than approximately 50 million people are estimated to be affected by dengue disease each year [1,2]. The dengue virus is transmitted to humans by mosquitoes, especially Aedes aegypti and Aedes albopictus. As at least four different serotypes of dengue viruses exist, people may be infected with dengue

∗

Corresponding author. Tel.: +86 0 2985310232. E-mail addresses: [email protected], [email protected] (S. Tang).

http://dx.doi.org/10.1016/j.nonrwa.2014.09.004 1468-1218/© 2014 Elsevier Ltd. All rights reserved.

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

237

more than once [3]. Currently, no specific antiviral therapy or vaccines are available against dengue [4,5], so vector population control remains the principal means of dengue prevention. However, traditional control measures, such as the use of insecticides for reducing mosquito populations are often prohibitively expensive, unsustainable, environmentally undesirable and may lead to insecticide resistance [4]. They are also failing to slow the current dengue pandemic, which has stimulated a search for novel technologies based on genetic manipulation of mosquito vectors to break dengue transmission cycles [6,7]. Wolbachia is a maternally transmitted endosymbiotic bacterium that is estimated to infect up to 65% of insect species and approximately 28% of mosquito species that have been surveyed [8,9]. It lives in the testes and ovaries of its hosts and can interfere with their reproductive mechanisms, inducting a variety of mosquito phenotypes, such as those with cytoplasmic incompatibility (CI), parthenogenesis, feminization of genetic males and so on. Moreover, appearances of these phenotypes depend on the host species and on Wolbachia types. The effect of CI causes embryos from females uninfected with Wolbachia to die when they are mated with infected males, whereas infected females are not affected in this manner [10,11]. In mosquito vectors, Wolbachia induced CI and matrilineal inheritance may have opposite effects, such as population extinction, coexistence or all of the uninfected population may be replaced by infected insects depending on the types of Wolbachia involved, the release methods and ratios of infected to uninfected populations. Some Wolbachia cannot only successfully spread within mosquito populations through CI but also stop their mosquito hosts from replicating and transmitting dengue virus [12,13]. In addition to mosquitoes, there are many other economically important insects that harbor Wolbachia that could be exploited in control measures including Simulium squamosum, a vector of onchocerciasis [14], and whiteflies Bemisia tabaci, agricultural pests and vectors of plant pathogens [15]. There are two ways to develop Wolbachia as a biological control agent against the dengue virus. One is a population suppression strategy whereby a population of mosquito vectors die out in small, pilot, experiments after releases of males infected with Wolbachia, such as during the successful suppression of Culex pipiens populations in field tests [16,17]. However, male-only releases would not be practical on the scale required for controlling a large area, or for eradication programmes [18]. This is partly because of the danger that females could be released accidentally as a part of a CI strategy that could permit the Wolbachia infection type to become established in the host field population. Another strategy shifts from population suppression to replacement by use of CI mechanisms and matrilineal inheritance, which can lead to populations of uninfected mosquitoes being replaced by infected ones, despite releases of only small numbers of the latter [12,13,19,20]. At present, two kinds of mathematical models, involving both discrete-time and continuous-time models, have been investigated for the spread of Wolbachia infection [21–27]. The Allee effect and founder control (i.e. one strain survives but another is always driven to extinction) have been demonstrated using a continuous-time model for the behavior of one or two strains of Wolbachia within a single well-mixed population. In addition, a discrete spatial model shows patchy persistence of the two strains but it was habitat characteristics rather than space itself that led to persistence [21,22]. Competitions between mutually incompatible strains were studied by [22,24–26]. A more complex outcome of CI in haplodiploid species was reported by [25]. Reproduction of mosquitoes has been assumed to be continuous throughout the year in most models [22,23,26,27]. However, many insects including mosquito populations exhibit what Cauchley [28] termed as a birth pulse growth pattern. Whilst dengue vectors may be at large in some areas throughout the year, where dengue transmission may thus be continuous, in most places both are highly seasonal with the transmission being interrupted. This is the case in China, where no cases occur between January and April before the outbreak season [29]. Dengue vectors have the ability to survive lean periods as eggs, even being able to hatch after periods of desiccation and so new populations emerge afresh at the start of rainy seasons and then give rise to waves of discrete cohorts. Even where the insects may be perennial, their abundances and dengue transmission are highly seasonal with peaks at the height of the rains or shortly after the rainy season. For instance, this is so in Bangladesh [30], India [31–33] and Thailand [34,35]. A female may survive for a couple of months under ideal temperature and moisture conditions. It needs about 12–18 days to go through four separate stages of its life cycle (egg, larva, pupa, and adult) depending on the environment. Here we focus on adult mosquitoes emerging periodically and ignore the details of their morphological transformations during their life cycles. Moreover fertility of mosquitoes is instantaneous in relation to their life span during tropical rainy seasons. For example every female mosquito oviposits two to three times, on average, and produces up to 500 eggs each year before she finally dies [36,37]. So births are seasonal and take place in a relatively short period at each oviposition. Therefore, it is appropriate to replace continuous breeding models with birth pulse models. Such a phenomenon can be well described by impulsive differential equations (hybrid dynamical systems) [38,39], types of equations which are found in almost every domain of the applied sciences [40–42]. The qualitative properties of these systems are embodied in those of equivalent discrete systems which determine the state after a pulse in terms of the state after the previous pulse. In this study, we investigate an impulsive system at fixed moments with birth-pulse and Wolbachia-induced CI in different density dependent death rate functions. The models represent mosquitoes as the principal vectors of dengue virus, with endosymbiotic Wolbachia bacteria capable of preventing the mosquitoes from reproducing and transmitting the virus. This can be achieved either by mosquito eradication or population replacement. Firstly, we obtain the stroboscopic maps of impulsive systems with different density dependent death functions. Secondly, for impulsive systems with different density dependent death rates, the qualitative properties of these systems are embodied in those of equivalent stroboscopic systems. However, as it is difficult to study the stability of equilibria for the stroboscopic systems directly, we investigate the existence and stability of equilibria for the systems by their equivalent systems. Moreover, we make some comparisons between original systems and equivalent ones under different density dependent death rate functions. Particularly, there

238

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

may be different dynamic behaviors depending on whether transmission from mother to offspring is perfect or not, so we also study this special case. Finally, we discuss some important biological considerations with different death rates. For the model with a strong density dependent death function, if the initial ratio of infected to total mosquitoes exceeds a critical value, strategies for population replacement may be partly realized. In particular, when the transmission is perfect, the periodic solution of impulsive systems derived from the local stability of an equilibrium is globally stable, hence the strategy of population replacement can be successful. However, with a weak density dependent death function, achieving population eradication is difficult due to the solutions of the system being sensitive to initial values. By using numerical methods we show that population eradication can be achieved, but only when parameters lie in particular regions and the initial density of the uninfected mosquitoes is low enough. 2. Mosquito population model with birth-pulse and CI In the absence of CI, we assume that the mosquito population can reproduce every period T which is density dependent, denoted as B(N (t )). Meanwhile we assume that mosquito deaths are continuous, denoted as F (N (t )). Therefore, the size of mosquito vector populations can be described by the following differential equation with birth-pulse



dN (t )

= −F (N (t )), t ̸= nT , dt N (nT + ) = B(N (nT ))N (t ), t = nT ,

(1)

where n (n = 1, 2, . . .) is a positive integer. Considering the density dependent death rate, function F (N (t )) often has the following two forms in biological literature [43,44]: (F1 ): F (N (t )) = δ1 N (t )2 with δ1 > 0; δ1 N ( t ) (F2 ): F (N (t )) = K + with δ1 > 0, K > 0. N (t ) According to the formulae of F1 and F2 we can see that the effect of density of the population on the death functions is much stronger in F1 than F2 . Thus, for convenience, we call F1 a strong density dependent function [43], while F2 is denoted as a weak density dependent function [44]. We aim to discuss how different density dependent death rate functions affect dynamic behaviors of systems and the strategies for preventing the spread of dengue virus. In this paper, the total population of mosquitoes N (t ) is divided into infected and uninfected classes, with the size of each class denoted by I (t ) and U (t ), respectively. We assume that the ratio of Wolbachia-infected males to infected females is the same as that of uninfected males to uninfected females. In general it is assumed that both infected and uninfected individuals have the same natural birth rate b and natural death rate δ (i.e. δ1 in the formulae of F1 and F2 ), but infected individuals may suffer an additional loss of natural fitness, denoted by D, which means that their reproductive capacity is not reduced in infected individuals compared with uninfected ones, while the death rate increases. Note that Wolbachia is primarily transmitted vertically, i.e. it is mostly passed from infected females to their offspring. Moreover, the transmission is usually not perfect in nature, occurring with a probability τ ∈ (0, 1]. The terms τ bI (nT ) and (1 −τ )bI (nT ) represent infected mothers that do and do not transmit Wolbachia bacteria to their offspring, respectively. The effect of the CI mechanism results in zygotic death of potential offspring with a probability q ∈ [0, 1] when an infected male is crossed with an uninfected female, and all other crosses are unaffected. The term b(1 − qI (nT )/N (nT ))U (nT ) refers to the quantity of mosquitoes left after the effect of CI. The infected and uninfected populations are decreased by an amount of (δ + D)I (t )N (t ) (or (δ + D)I (t )/(K + N (t ))) and δ U (t )N (t ) (or δ U (t )/(K + N (t ))), respectively, whenever t ̸= nT . Birth pulses with the CI mechanism occur at t = nT . Therefore, the dynamics of the mosquito population with a strong density dependent death rate function F1 are governed by the following pulse differential equations with fixed moments.

  dI (t )     = −(δ + D ) I ( t ) N ( t ),   dt  t ̸= nT ,     dU (t )  = −δ U (t )N (t ), dt    I (nT + ) = (1 + τ b)I (nT ),       I (nT )  +   U (nT ) = U (nT ) + (1 − τ )bI (nT ) + b 1 − q U (nT ), N (nT )

(2) t = nT .

Similarly, the dynamics of the mosquito population with a weak density dependent death rate function F2 are described as follows:

  dI (t ) I (t )     = −(δ + D ) ,   dt K + N (t )  t ̸= nT ,     dU (t ) = −δ U (t ) ,  dt K + N (t )    I (nT + ) = (1 + τ b)I (nT ),       I (nT )  +   U (nT ) = U (nT ) + (1 − τ )bI (nT ) + b 1 − q U (nT ), N (nT )

(3) t = nT .

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

239

Since we mainly investigate the effects of birth pulses and CI on the influence of control strategies, we let D = 0 for simplicity, which does not change the qualitative behaviors of the above two systems. 3. Stroboscopic map of systems Note that the stroboscopic maps which determine the numbers of infected and uninfected populations immediately after each birth pulse at the discrete times nT play a key role in investigating the dynamics of the proposed models (2) and (3). Thus in the following, we deduce them first. 3.1. Stroboscopic map of system (2) Adding and dividing the first two equations of (2) yields

 dN (t )   = −δ N (t )2 ,  dt

(4)

dU (t ) U (t )   = .  dI (t ) I (t ) Solving the above equations in any birth pulse interval (nT , (n + 1)T ] gives N (nT + )

   N ( t ) =

   ,

N (nT + )δ(t − nT ) + 1 U (nT + )     I (t ), U (t ) =  I (nT + )

t ∈ (nT , (n + 1)T ].

(5)

Owing to N (t ) = I (t ) + U (t ), we have I (nT + )

   I ( t ) =

, δ(t − nT )(U (nT + ) + I (nT + )) + 1 + U (nT )   . U (t ) = δ(t − nT )(U (nT + ) + I (nT + )) + 1

(6)

Denote In = I (nT + ), Un = U (nT + ), and by simplification we have the following discrete equations which is the stroboscopic map of system (2).

(1 + τ b)In , f1 (In , Un ), δ T (In + Un ) + 1   In Un (1 − τ )bIn + 1 + b − bq , g1 (In , Un ). = δ T ( In + U n ) + 1 In + Un δ T (In + Un ) + 1

    In + 1 =   Un+1

(7)

3.2. Stroboscopic map of system (3) In order to obtain the stroboscopic map of system (3), we need to define the Lambert W function and its properties, so we introduce it next. Definition. The Lambert W -function is defined to be a multivalued inverse of the function z → zez satisfying Lambert W(z ) exp(Lambert W(z )) = z. Its derivative with respect to z is LambertW′ (z ) =

LambertW(z ) z (1 + LambertW(z ))

.

(8)

Note that the function z exp(z ) has the positive derivative (z + 1) exp(z ) if z > −1. Define the inverse function of z exp(z ) restricted on the interval [−1, ∞) to be LambertW(0, z ). Similarly, we define the inverse function of z exp(z ) restricted on the interval (−∞, −1] to be LambertW(−1, z ). In the following, for convenience and simplification, we denote LambertW(0, z ) and LambertW(−1, z ) as W (z ) and W (−1, z ), respectively. According to the first two equations of (3), we can get

 dN (t )   = −δ dt  dU (t )



dI (t )

=

N (t )

K + N (t ) U (t ) I (t )

.

, (9)

240

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

The first equation of (9) is independent of N (t ), so we can integrate it directly. In any impulsive interval (nT , (n + 1)T ], the so-called integrated form with the initial value N0 = N (nT + ) is



1

δ

N (nT ) − N (t ) + K ln +



N (nT + )



N (t )

= t − nT .

Rearranging the above equation yields: N (t )



K

ln

δ

 +

N (nT + )

1

δ

N (t ) − N (nT + ) = −(t − nT ),



i.e., N (t )

 ln

 +

N (nT + )

N (t ) K

=

N (nT + ) − δ(t − nT ) K

.

Taking exponents on both sides of the above equation gives N (t ) K

 exp

N (t )

 =

K

N (nT + ) K

 exp

N (nT + ) − δ(t − nT )



K

.

According to the definition of the Lambert W function, we have the following analytical solution for model (9), i.e.

     N (nT + ) N (nT + ) − δ(t − nT )    exp , N (t ) = KW  K

K

t ∈ (nT , (n + 1)T ].   U (nT + )   U ( t ) =  I ( t ). I (nT + ) Based on N (t ) = I (t ) + U (t ), we have     I (nT + ) + U (nT + ) I (nT + ) + U (nT + ) − δ(t − nT ) KI (nT + )   W exp , I ( t ) =  I (nT + ) + U (nT + ) K K     I (nT + ) + U (nT + ) I (nT + ) + U (nT + ) − δ(t − nT ) KU (nT + )  U ( t ) = W exp . I (nT + ) + U (nT + ) K K Denoting In = I (nT + ), Un = U (nT + ), we have the following discrete system  (1 + bτ )KIn   W (z ) , f2 (In , Un ), In+1 = In + Un   In KUn (1 − τ )bKIn   Un+1 = W (z ) + 1 + b − bq W (z ) , g2 (In , Un ), In + U n In + Un In + Un

(10)

with z,

In + Un K

 exp

In + Un − δ T



K

.

4. Equilibria and equivalent systems In order to investigate the existence of periodic solutions of systems (2) and (3), we first get the equilibria of equivalent systems (7) and (10) which are determined by both equations fi (I , U ) − I = 0 and gi (I , U ) − U = 0 (i = 1, 2). Note that both lim

(fi (I , U ) − I ) = 0 and

(I ,U )→(0,0)

lim

(gi (I , U ) − U ) = 0

(I ,U )→(0,0)

hold true, thus we define fi (0, 0) = gi (0, 0) = 0 (i = 1, 2). Clearly both fi and gi are continuous on the closure of R2 = {(I , U )|I ≥ 0, U ≥ 0}. Therefore, trivial equilibria Ei0∗ = (0, 0) (i = 1, 2) always exist in systems (7) and (10). 4.1. Equilibria of systems (7) and (10) For system (7), solving equations f1 (I , U )− I = 0 and g1 (I , U )− U = 0, except for equilibrium (0, 0), there are three other ∗ ∗ equilibria, such as one boundary equilibrium E11 = (0, U11 ) = (0, b/(δ T )), two interior equilibria, denoted by E1j∗ = (I1j∗ , U1j∗ )

(j = 2, 3). Let q2 − 4q(1 − τ ) , ∆ ≥ 0 (i.e., q > 4(1 − τ )), then  √  bτ q − ∆ ∗ ∗ I12 = U13 = , 2qδ T   √ bτ q + ∆ ∗ ∗ I13 = U12 = . 2qδ T

(11)

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

241

Similarly, for system (10), solving equations f2 (I , U ) − I = 0 and g2 (I , U ) − U = 0, except for equilibrium (0, 0), there are also a boundary equilibrium ∗ ∗ E21 = (0, U21 )=



0,

(1 + b)(δ T − K ln(1 + b)) b



,

∗ and interior equilibria E2j = (I2j∗ , U2j∗ ) (j = 2, 3). If ∆ ≥ 0, then



∗ ∗ I22 = U23 =

M (τ ) q − 2q



∗ ∗ I23 = U22 =

√  ∆

M (τ ) q +

√  ∆

2q

, (12)

,

where

(1 + bτ )(δ T − K ln(1 + bτ )) ∗ M (τ ) = I2j + U2j∗ = bτ ∗ is deduced from the equation f2 (I , U )−I = 0 with I ̸= 0 and U ̸= 0. Therefore, equilibrium E21 is positive provided M (1) > 0, ∗ ∗ ∗ i.e., b < exp(δ T /K ) − 1, while E2j = (I2j , U2j ) (j = 2, 3) are positive if and only if M (τ ) > 0, i.e., b < (exp(δ T /K ) − 1)/τ . Taking Poincare transformation xn = Un /In and yn = 1/In , system (10) is changed to  (1 + b)x2n + (2b − bτ − bq + 1)xn + b − bτ  xn+1 = , f3 (xn ), (1 + bτ )(1 + xn ) 1 + xn  yn+1 = , g3 (xn , yn ), (1 + bτ )KW (z1 ) with z1 =

1 + xn Kyn

 exp

1 + xn Kyn

−

δT K



.

Solving equation f3 (x) = x gives x∗1,2 =

q − 2(1 − τ ) ± 2(1 − τ )

√ ∆

.

Based on lim

(x,y)→(x∗1,2 ,0)

(g3 (x, y) − y) = 0,

there are two infinite singular points A1 (x∗1 , 0) and A2 (x∗2 , 0). Therefore, infinity equilibrium E∞ = (+∞, +∞) exists in system (10). 4.2. Equivalent systems of systems (7) and (10) ∗ (i = 1, 2) and E∞ cannot be linearized in systems (7) and (10), hence their local stability On the one hand, equilibria Ei0 cannot be studied directly. On the other hand, as Jacobian matrices of Eij∗ (i = 1, 2, j = 2, 3) have complex forms, it is ∗ difficult to determine whether or not the modulus of the eigenvalues is less than one. So to understand the stability of Ei0 and Eij∗ (i = 1, 2, j = 1, 2, 3), we expand systems (7) and (10) on a whole axis by studying the transformed systems (ξn , Un ) with ξn = In /(In + Un ) ∈ [0, 1] (the ratio of infected individuals to the total population). Then every nonzero equilibrium Eij∗ (i = ∗ 1, 2, j = 1, 2, 3) will possess one corresponding equilibrium in the transformed systems. However Ei0 (i = 1, 2) (or E∞ ) may have more than one corresponding equilibrium (all of them are zero (or infinity) in the second variable Un ) in the transformed systems, which leads to the richness of the dynamic behaviors and many difficulties in the analysis of systems (7) and (10).

4.2.1. The equivalent system of system (7) Based on the above conversion of variables, the equivalent system (ξn , Un ) of system (7) has the following form

(1 + τ b)ξn , H1 (ξn ), 1 + b − bqξn (1 − ξn ) P (ξn )Un   Un+1 = , G1 (ξn , Un ), δ TUn − ξn + 1    ξn+1 =

where P (ξn ) , bqξn2 − (bq + bτ + 1)ξn + b + 1. In order to ensure that it is meaningful for the above system, we first need to prove the following results.

(13)

242

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

Claim. H1 (ξn ) ∈ [0, 1] holds true with ξn ∈ [0, 1] in the parameter space (here 0 ≤ q, τ ≤ 1 and b > 0). Proof. For convenience, we denote ξn as x, then H1 (ξn ) becomes as H1 (x). Due to H1 (x) ≥ 0 with x ∈ [0, 1], the proof of the Claim is transformed to prove inequality H1 (x) ≤ 1, i.e. P (x) ≥ 0, for all of x ∈ [0, 1]. Note that symmetric axis of P (x) is x∗ =

bq + bτ + 1 2bq

.

When x∗ > 1, then P (1) = b(1 − τ ) ≤ P (x) ≤ b + 1 = P (0), so P (x) ≥ 0 with x ∈ [0, 1]. When x∗ ≤ 1, i.e. bτ ≤ bq − 1, then discriminant ∆1 = (bq + bτ + 1)2 − 4bq(b + 1) ≤ 0, so P (x) ≥ 0 with x ∈ [0, 1]. Therefore, H1 (x) ≤ 1 holds true for all of x ∈ [0, 1] in parameter space. We complete the proof of the claim. For system (13), there are six equilibria at most as shown by

   b ∗ ∗ ∗ ∗  ¯ ¯ , E11 = (ξ11 , U11 ) = 0,  E10 = (0, 0),  δT      √ √    ∗ q− ∆ q+ ∆ ∗ ∗ ∗ ∗ ∗ ∗ ∗ E¯ 12 = (ξ12 , U12 ) = , U12 , E¯ 13 = (ξ13 , U13 ) = , U13 , 2q 2q       √ √    q− ∆ q+ ∆  ∗ ∗ ∗ ∗  ¯ ¯ E = (ξ , 0 ) = , 0 , E = (ξ , 0 ) = ,0 .  14 12 15 13 2q

2q

∗ ∗ All of them correspond to the equilibria in system (7). For example, equilibrium E¯ 1j corresponds with E1j (j = 1, 2, 3), ∗ ¯∗ ∗ ∗ whereas the other three equilibria E¯ 10 , E14 and E¯ 15 correspond to E10 .

4.2.2. The equivalent system of system (10) Similarly, the equivalent system of system (10) is

(1 + τ b)ξn , H1 (ξn ), 1 + b − bqξn (1 − ξn )  Un+1 = P (ξn )KW (z2 ) , G2 (ξn , Un ),  ξ

n +1

=

(14)

where z2 ,

Un K (1 − ξn )

 exp

Un − δ T (1 − ξn )



K (1 − ξn )

≥ 0.

There are six finite equilibria at most as follows

 ∗ ∗ ∗ ∗ E¯ 20 = (0, 0), E¯ 21  21 , U21 ) =   (0, M (1)),    = (ξ√ √   q+ ∆ q− ∆ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ , U22 , , U23 , E¯ 23 = (ξ23 , U23 ) = E¯ 22 = (ξ22 , U22 ) =  2q 2q   ¯∗ ∗ ∗ ∗ E24 = (ξ22 , 0), E¯ 25 = (ξ23 , 0). ∗ ∗ Moreover, there are three infinite equilibria E¯ 0∞ = (0, +∞), E¯ 1∞ = (ξ22 , +∞) and E¯ 2∞ = (ξ23 , +∞). All of them corre∗ ¯∗ ∗ ∗ ¯ , E24 spond to the equilibria in system (10). Specifically, equilibrium E2j corresponds to E2j (j = 1, 2, 3). Three equilibria E¯ 20 ∗ ∗ correspond to E20 , while three infinite equilibria E¯ 0∞ , E¯ 1∞ and E¯ 2∞ correspond to E∞ . and E¯ 25

Remark 1. When ∆ > 0, then 0 ≤ ξi2∗ < 1/2 < ξi3∗ ≤ 1 (i = 1, 2); ξi2∗ and ξi3∗ collide together as 1/2 when ∆ = 0; ξi2∗ and ξi3∗ disappear when ∆ < 0. As every equilibrium in system (7) (or (10)) has its corresponding equilibrium in equivalent system (13) (or (14)), the stability of the equilibria in systems (7) and (10) is determined by the stability of equilibria in their equivalent systems. In the following, the stability of equilibria for equivalent systems (13) and (14) is investigated first. 5. Stability of equilibria In the neighborhood of equilibria, the dynamics of discrete systems (13) and (14) are determined by the linearization Xn+1 = JXn

(15)

with Jacobian matrix J as the linear counterpart of (13) and (14), respectively, and X = (ξ , U ). When the modulus of eigenvalues for matrix J is less than one, then the equilibria must be locally stable. For convenience, the two eigenvalues of ¯ ′ij and λ¯ ′′ij (i = 1, 2, j = 0, 1, 2, 3, 4, 5), respectively. equilibrium E¯ ij∗ for systems (13) and (14) are denoted as λ Whether the transmission rate τ is perfect or not may result in different stability and dynamic behaviors, so we will focus on the two cases, i.e., τ ∈ [0, 1) and τ = 1 in the following.

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

243

5.1. The transmission being not perfect, i.e. τ ∈ [0, 1) For equivalent system (13), note that H1 (ξ ) is independent of variable U, taking the derivative of H1 (ξ ) with respect to

ξ , yields

dH1 (ξ )

(1 + τ b)(1 + b − bqξ 2 ) . dξ (1 + b − bqξ (1 − ξ ))2 Solving dH1 (ξ )/dξ = 0 with respect to ξ yields one positive root, denoted by  ξmax =

=

1

bq

+

1 q

(16)

> 1.

Moreover taking the partial derivative of G1 (ξ , U ) with respect to U, we have

∂ G1 (ξ , U ) P (ξ )(1 − ξ ) = . ∂U (δ TU + 1 − ξ )2 Therefore, the Jacobian matrix J at (ξn , Un ) is given by   (1 + bτ )(1 + b − bqξn2 ) 0  (1 + b − bqξ (1 − ξ ))2  , n n J (ξn , Un ) =   P (ξn )(1 − ξ )  ⋆ (δ TU + 1 − ξ )2 where (⋆) is not necessary for the stability analysis, so we omit it. ∗ The Jacobian matrix of system (13) at E¯ 11 takes the form of   1 + τb 0   b+1  1 . ⋆ b+1

(17)

(18)

Since inequalities

     1  1 + τb ′ ′′    <1 ¯ ¯ |λ11 | =  < |λ11 | =  b + 1 b+1  ∗ ∗ hold true, E¯ 11 is locally stable naturally in system (13). This indicates that equilibrium E11 is also locally stable in system (7), which means that given some initial value, Wolbachia bacteria do not spread in the mosquito population, and only an uninfected population exists at the end. ∗ ∗ ∗ To study the stability of equilibria E¯ 1j (j = 2, 3), we substitute values ξ12 and ξ13 into Eq. (16), to get

       dH1 (ξ )  1 + = |λ¯ ′12 | =    dξ ∗  ξ =ξ12        dH1 (ξ )  1 + = |λ¯ ′13 | =    dξ  ξ =ξ ∗ 13



b 4 − q − 4τ +





q2 − 4q(1 − τ ) 

,      b 4 − q − 4τ − q2 − 4q(1 − τ )  .  2(1 + bτ )  2(1 + bτ )

(19)

It is easy to see that inequalities

−4(1 + bτ ) b

< 4 − q − 4τ −



q2 − 4q(1 − τ ) < 0,

and 4 − q − 4τ +



q2 − 4q(1 − τ ) > 0

∗ ¯ ′12 | > 1 and |λ¯ ′13 | < 1. Therefore, E¯ 12 hold true in parameter space, provided that ∆ > 0 (i.e. q > 4(1 − τ )), then we have |λ ∗ for system (13) is unstable. In order to determine the stability of equilibrium E¯ 13 , substituting it into Eq. (17) gives

    ∗ ∗   ∂ G1 (ξ13 , U13 )   1  |λ¯ ′′13 | =  = (20)   1 + bτ  < 1. ∂U ∗ ∗ ∗ Thus equilibrium E¯ 12 is unstable, while E¯ 13 is locally stable for system (13) once it exists. Thereby, for system (7), E12 is

∗ unstable, while E13 is locally stable, which indicates that given some initial values, infected and uninfected mosquitoes may coexist in the end.

244

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

Fig. 1. Two-parameter bifurcation analysis of existence of equilibria for systems (7) and (13).

∗ The Jacobian matrix J of system (13) at E¯ 10 takes the form of

1 + τb 1+b

⋆

0

 .

1+b

Since inequalities

  1 + τb  < 1 and |λ¯ ′′ | = |1 + b| > 1 |λ¯ ′10 | =  10 1+b  ∗ hold true, then E¯ 10 is unstable. ∗ ∗ ∗ ¯ ′14 = λ¯ ′12 and λ¯ ′15 = λ¯ ′13 , then |λ¯ ′14 | > 1 and |λ¯ ′15 | < 1. Substituting E¯ 15 For the stability of E¯ 14 and E¯ 15 , due to λ into (17) gives

λ¯ ′′15 =

∗ P (ξ13 ) ∗ 1 − ξ13

= bτ + 1 > 1 .

∗ ∗ So both E¯ 14 and E¯ 15 are unstable. ∗ ∗ ∗ ¯∗ is unstable, which means that with any , yields E10 Based on the above stability analysis with respect to E¯ 10 , E14 and E¯ 15 initial values, both infected and uninfected mosquitoes cannot die out in the end. Hence, a strategy of aiming to eradicate mosquitoes will fail with a strong density dependent death rate. In order to understand the existence and stability of equilibria of systems (7) and (13) in more detail, we choose τ and q as bifurcation parameters to address the regions of existence and nonexistence of equilibria for these systems as shown in Fig. 1. The three lines (L1 , L2 and L3 ) divide the τ and q parameter space into five regions. Only Ω11 and Ω12 are meaningful ∗ ∗ ∗ ∗ ) coexist in region Ω11 , while all of the equilibria coexist in and E11 (or E¯ 10 and E¯ 11 regions for the existence of equilibria. E10 ∗ ∗ ∗ ∗ ∗ ∗ and E¯ 13 or E¯ 14 and E¯ 15 ) may coincide in line L3 . region Ω12 . E12 and E13 (or E¯ 12 According to the above stability analysis, we have the following main results.

Theorem 1. For system (13) with the transmission being not perfect, i.e. τ ∈ [0, 1), and the rate of zygotic death for CI q ∈ [0, 1], ∗ ¯∗ ¯∗ ¯∗ ¯∗ ∗ there are at most six equilibria, E¯ 10 , E11 , E12 , E13 , E14 and E¯ 15 . ∗ ∗ ∗ (1) If ∆ > 0 (i.e. Ω12 ), then E¯ 11 and E¯ 13 are locally stable, while the other four are unstable. The solution stabilizes at E¯ 11 provided ∗ ∗ ¯ that the initial value (ξ0 , U0 ) ∈ (0, ξ12 ) × (0, +∞); the solution stabilizes at E13 provided that the initial value (ξ0 , U0 ) ∈ ∗ (ξ12 , 1) × (0, +∞). ∗ ∗ ∗ ∗ (2) If ∆ = 0 (i.e. L3 ), then E¯ 12 and E¯ 13 collide together as (1/2, bτ /(2δ T )), while E¯ 14 and E¯ 15 collide together as (1/2, 0), and ∗ only E¯ 11 is locally stable. ∗ ∗ (3) If ∆ < 0 (i.e. Ω11 ), then E¯ 1j (j = 2, 3, 4, 5) disappear, and E¯ 11 is locally stable.

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

A

E

B

F

C

G

D

H

245

Fig. 2. The effects of different initial values on the solutions of systems (7) and (13). The parameters are fixed as follows: τ = 0.9, b = 2, δ = 0.3, T = ∗ ∗ ∗ 2, q = 0.5, then ξ12 = 0.2764 and ξ13 = 0.7236. (A)–(D): The solutions of systems (13) and (7) from initial values (ξ0 , U0 ) = (0.2, 0.4) (ξ0 < ξ12 ) and ∗ (I0 , U0 ) = (0.1, 0.4). (E)–(H): The solutions of systems (13) and (7) from initial values (ξ0 , U0 ) = (0.4, 0.6) (ξ0 > ξ12 ) and (I0 , U0 ) = (0.4, 0.6).

Theorem 2. For system (7) with the transmission being not perfect, i.e. τ ∈ [0, 1), and the rate of zygotic death for CI q ∈ [0, 1], ∗ ∗ ∗ ∗ there are at most four equilibria, E10 , E11 , E12 and E13 . ∗ ∗ ∗ (1) If ∆ > 0 (i.e. Ω12 ), then E11 and E13 are locally stable, while the other two are unstable. The solution stabilizes at E11 provided ∗ ∗ ∗ initial ratio ξ0 ∈ (0, ξ12 ) and any U0 ; the solution stabilizes at E13 provided initial ratio ξ0 ∈ (ξ12 , 1) and any U0 . ∗ ∗ ∗ and E13 collide together as (bτ /(2δ T ), bτ /(2δ T )), and only E11 is locally stable. (2) If ∆ = 0 (i.e. L3 ), then E12 ∗ (3) If ∆ < 0 (i.e. Ω11 ), then two interior equilibria disappear, and E11 is locally stable.

Remark 2. Based on Theorem 1, system (13) clearly demonstrates that there is an Allee effect such that a small level of infection cannot invade an uninfected population as shown in Fig. 3. That is to say low density populations are driven to ∗ extinction, but high density populations can survive. Because when initial value ξ0 is sufficiently small (ξ0 < ξ12 ), function H1 (ξn ) is monotonically decreasing with respect to ξn , and the level of infection will decrease to zero. It is only when the level ∗ of infection reaches some critical threshold ξ12 that function H1 (ξn ) is monotonically increasing with respect to ξn , and the ∗ . The effects of different initial values on the solutions infection will break out, and the level of infection will stabilize at ξ13 of systems (7) and (13) are shown in Fig. 2. From the relations between system (13) and system (2), the existence of equilibria of system (13) indicates the existence of periodic solutions of system (2). When the transmission is not perfect, i.e. τ = 0.9, the periodic solutions of system (2) from different initial values are provided in Fig. 4(A–D). In the following, we investigate the stability of equilibria for systems (10) and (14). For system (14), according to the definition of the Lambert W function, the Jacobian matrix is given by

(1 + bτ )(1 + b − bqξn2 )  (1 + b − bqξ (1 − ξ ))2 n n J (ξn , Un ) =   ⋆





0

 , ∂ W (z2 )  KP (ξn ) ∂ Un

(21)

where

∂ W (z2 ) ∂ W (z2 ) ∂ z2 = ∂ Un ∂ z2 ∂ Un  exp

=

Un −δ T (1−ξn ) K (1−ξn )



n 1 + K (1U−ξ n)



exp(W (z2 ))(1 + W (z2 ))K (1 − ξn )

,

(z2 ≥ 0, i.e. Un ≥ 0).

(22)

246

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

Fig. 3. Cobweb model to show the stability of the first equation of discrete systems (13) and (14). The parameter values are fixed as follows: τ = 0.8, b = 50, q = 0.95.

A

B

C

D

E

F

Fig. 4. The effects of different initial values on the periodic solutions of system (2). The baseline parameter values and initial values are the same as Fig. 2. (A)–(D): When parameter value τ = 0.9, the period solutions of system (2). (E)–(F): When parameter value τ = 1, the periodic solution of system (2) with initial value (I0 , U0 ) = (0.1, 0.4).

Especially, when z2 > 0, i.e. Un > 0, from the properties of the Lambert W function, Eq. (22) can be transformed to exp



Un −δ T (1−ξn ) K (1−ξn )





n 1 + K (1U−ξ W (z2 ) n)

z2 (1 + W (z2 ))K (1 − ξn )

.

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

247

Due to G2 (ξn , Un ) = Un , we have W (z2 ) = Un /(KP (ξn )), then

  Un P (ξn ) K + 1−ξ ∂ W (z2 ) n KP (ξn ) = . ∂ Un KP (ξn ) + Un ∗ ∗ The stability of equilibrium E¯ 2j for system (14) is determined by the eigenvalues of the matrix J (E¯ 2j ) (j = 0, 1, 2, 3, 4, 5). ∗ Substituting equilibrium E¯ 21 into (21) gives the following matrix   1 + bτ 0   ∗ J (E¯ 21 )= 1+b (23) (b + 1)(K + M (1))  . ⋆ K (b + 1) + M (1) The modulus of eigenvalues of the above matrix is

    1 + τb  (b + 1)(K + M (1))  ′ ′′    > 1,  ¯ ¯ |λ21 | =  < 1 and |λ21 | =  1+b  K (b + 1) + M (1)  ∗ ∗ then E¯ 21 is unstable. So E21 for system (10) is also unstable. ∗ ∗ Similarly, for the stability of equilibria E¯ 22 and E¯ 23 , owing to λ′22 = λ′12 and λ′23 = λ′13 , then |λ′22 | > 1 and |λ′23 | < 1. ∗ Furthermore, substituting E¯ 23 into (21) gets another eigenvalue

λ′′23 =

∗ P (ξ23 )(K + M (τ )) ∗ ∗ P (ξ23 )K + U23

> 0.

We have |λ′′23 | > 1, provided that ∗ ∗ P (ξ23 )M (τ ) − U23 =

 √  (1 + bτ )(δ T − K ln(1 + bτ )) q − ∆ 2q

>0

∗ ∗ ∗ ∗ holds true. Then both E¯ 22 and E¯ 23 are unstable when they exist. Thus both interior equilibria E22 and E23 for system (10) are unstable. ∗ For stability of equilibrium E¯ 20 , substituting it into (21) gives

 1 + bτ  1+b

J (E¯ 20 ) = 

⋆



0

(b + 1) exp



 −δ T  .

(24)

K

The modulus of eigenvalues of the above matrix is

       1 + τb  < 1 and |λ¯ ′′ | = (b + 1) exp −δ T  . |λ¯ ′20 | =  20    1+b K ∗ Thus E¯ 20 is locally stable provided that b < exp(δ T /K ) − 1 holds true. ∗ ∗ ∗ ¯ ′24 | > 1 and |λ¯ ′25 | < 1. Substituting E¯ 25 For the stability of E¯ 24 and E¯ 25 , because of λ′24 = λ′12 and λ′25 = λ′13 , we have |λ into (21) gives

λ¯ ′′25 =

  ∗ ) exp −δK T P (ξ25 ∗ 1 − ξ25

= (bτ + 1) exp



−δ T K



> 0.

∗ ∗ is locally stable, provided that So E¯ 24 is unstable, while E¯ 25

b<

  exp δKT − 1 τ

holds true. ∗ ∗ ¯∗ ∗ may be locally stable if b < exp(δ T /K ) − 1 holds true. Therefore, with , E24 and E¯ 25 , E20 From the stability analysis of E¯ 20 a weak density dependent death rate, the strategy of aiming to eradicate mosquitoes may come true. Similarly, in order to understand the existence and stability of equilibria of systems (10) and (14) in more detail, we choose τ , q and b as bifurcation parameters, all others fixed as those in Fig. 5, to determine the regions of existence and stability of equilibria for these systems. Then three curves with respect to parameters τ , q and b are defined as follows:

     δT   Π21 = (τ , q, b)| b = exp −1 ,   K      1 δT  Π = (τ , q , b )| b = exp − 1 , 22   τ K   Π23 = {(τ , q, b)| q − 4(1 − τ ) = 0}.

(25)

248

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

Fig. 5. Three-parameter bifurcation analysis of existences of equilibria for systems (10) and (14). The parameter values are fixed as follows: δ = 0.3, T = 4, K = 1. Table 1 Analysis of existence and stability of equilibria for system (14). Equilibrium

∗ E¯ 20

∗ E¯ 24

∗ E¯ 25

∗ E¯ 21

∗ E¯ 22

∗ E¯ 23

E¯ 0∞

E¯ 1∞

E¯ 2∞

Ω21 Ω22 Ω23

S U U

U U U

S S U

U

U U

U U

×

×

S S S

U U U

S S S

Ω24 Ω25 Ω26

U U S

× × ×

× × ×

× ×

× × ×

× × ×

S S S

× × ×

× × ×

× ×

U

S∼Stable, U∼Unstable, ×∼Nonexistence.

The three curves (Π21 , Π22 and Π23 ) divide τ , q and b parameter space into six regions. The existence, nonexistence or ∗ ¯∗ coexistence of equilibria of systems (10) and (14) in each region are shown in Fig. 5. Moreover, on the curve Π21 , E21 (E21 ) ∗ ∗ ∗ ∗ ∗ ∗ ∗ ¯∗ ¯ ¯ ¯ can coincide with E20 (E20 ); on the curve Π22 , E22 , E23 and E20 can coincide (E22 coincides with E24 , while E23 coincides with ∗ ∗ ∗ ∗ ∗ ∗ and coincide as (1/2, M (τ )/2), while E¯ 24 and E23 can coincide as (M (τ )/2, M (τ )/2) (E¯ 22 and E¯ 23 ); on the curve Π23 , E22 E¯ 25 ∗ ¯E25 coincide as (1/2, 0)). According to the above stability analyses of systems (10) and (14), we have the following main results. ∗ Theorem 3. For system (14) with the transmission being not perfect, i.e. τ ∈ [0, 1), there are at most six finite equilibria E¯ 2j (j =

0, 1, 2, 3, 4, 5), and three infinite equilibria E¯ 0∞ E¯ 1∞ and E¯ 2∞ . ∗ Theorem 4. For system (10) with the transmission being not perfect, i.e. τ ∈ [0, 1), there are at most four finite equilibria E2j (j = 0, 1, 2, 3), and one infinite equilibrium E∞ .

The existence and stability of equilibria of systems (10) and (14) are shown in Tables 1 and 2. ∗ ∗ Remark 3. For system (14), multiple attractors can coexist, such as the pair E¯ 20 and E¯ 25 and another pair E¯ 0∞ and E¯ 2∞ , which ∗ are separated by the line ξ = ξ22 as shown in Figs. 6–8. ∗ Remark 4. According to Tables 1 and 2, it is clear that E20 is locally stable in regions Ω21 , Ω22 and Ω26 , which implies that given some initial values (I0 , U0 ), (In , Un ) → (0, 0) as n → ∞. There are some differences among these regions. Overlook∗ ∗ ing the local stability of infinite equilibria, both E¯ 20 and E¯ 25 are locally stable in Ω21 , so some initial values (I0 , U0 ) make ∗ ∗ (ξn , Un ) → (0, 0), and other initial values make (ξn , Un ) → (ξ23 , 0) as n → ∞. Only E¯ 25 is locally stable in region Ω22 ,

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

249

Table 2 Analysis of existence and stability of equilibria for system (10). Equilibrium

∗ E20

∗ E21

∗ E22

∗ E23

E∞

Ω21 Ω22 Ω23 ,Ω24 ,Ω25 Ω26

S S U S

U

U U

U U

× ×

× ×

S S S S

× × U

S∼Stable, U∼Unstable, ×∼Nonexistence.

A

B

C

D

Fig. 6. The effects of different initial values on the solutions of systems (10) and (14) in region Ω21 . The parameter values are fixed as follows: ∗ ∗ τ = 0.8, b = 2, δ = 0.3, T = 4, q = 0.9, K = 1, interior equilibria of system (14) are ξ22 = 0.3333 and ξ23 = 0.6667.

∗ ∗ so some initial values (I0 , U0 ) make (ξn , Un ) → (ξ23 , 0) as n → ∞; while only E¯ 20 is locally stable in region Ω26 , so some initial values make (ξn , Un ) → (0, 0) as n → ∞.

Remark 5. The stability of infinite equilibrium E∞ cannot be determined, as the second eigenvalues of A1 (x∗1 , 0) and A2 (x∗2 , 0) cannot be calculated. Three infinite equilibria correspond to E∞ of system (10). Unfortunately, none of them can be investigated by the Jacobian matrix. But numerical simulations show local stabilities of E¯ 0∞ and E¯ 2∞ lead to E∞ being locally stable as shown in Figs. 6–8. There are two routes for the local stability of E∞ . Given some initial values (I0 , U0 ), one ∗ ∗ , +∞) as n → ∞. Especially, in regions Ω23 , Ω24 and Ω25 , E20 is is (ξn , Un ) → (0, +∞), and another is (ξn , Un ) → (ξ23 unstable, so E∞ may be globally stable, which make both of the two strategies fail. ∗ ∗ Remark 6. When the two equilibria E¯ 20 and E¯ 21 lie in the curve Π23 , the stabilities of them are not affected. For the stability of (1/2, M (τ )/2) and (1/2, 0), according to Eq. (19) and the forms of λ′′23 and λ′′25 , there are λ′22 = λ′23 = λ′25 = 1,

λ23 = ′′

P

(K + M (τ )) >0 1 M (τ ) P 2 K+ 2

P

1

1 2

and

λ25 = ′′

2

 T exp −δ K

1−

1 2

= (bτ + 1) exp



Note that

  P

1

2

M (τ ) −

M (τ ) 2

=

1 2

(bτ + 1) > 0,

−δ T K



> 0.

250

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

A

B

C

D

Fig. 7. The effects of different initial values on the solutions of systems (10) and (14) in region Ω22 . Except b = 2.5, the other parameter values are the same as Fig. 6. Table 3 Differences between systems (7) and (10). Equilibrium

∗ Ei0

∗ Ei1

∗ Ei2

∗ Ei3

E∞

System (7) (i = 1) System (10) (i = 2)

U S or U

S U or ×

U U or ×

S U or ×

S

×

S∼Stable, U∼Unstable, ×∼Nonexistence.

so λ′′23 > 1. Then (1/2, M (τ )/2) is unstable, while (1/2, 0) is stable, but not asymptotic stability, provided b<

  exp δKT − 1 τ

∗ ∗ holds true. Therefore, the only difference is that instability of E¯ 24 and stability of E¯ 25 are compromised as stability, but not asymptotically stable compared with other regions.

Numerical simulations of solutions of systems (10) and (14) in regions Ω21 , Ω22 and Ω23 (appropriate parameter values chosen from Fig. 5) starting from different initial values are provided in Figs. 6–8. Fig. 6(A–D) show the effects of different initial values on the solutions of systems (10) and (14) in region Ω21 . Fig. 6(A) demonstrates that solutions of system (14) ∗ ∗ ∗ ∗ (ξ0 < ξ22 ) or E¯ 25 (ξ0 > ξ22 ) with different initial values, i.e., (0.2, 0.1), (0.25, 0.1), (0.3, 0.1), (0.35, 0.1), stabilize at E¯ 20 ∗ corresponding to (0.4, 0.1) and (0.45, 0.1), respectively. Fig. 6(B) demonstrates that solutions of system (10) stabilize at E20 ∗ ∗ ¯ ¯ Fig. 6(A). Fig. 6(C) demonstrates that solutions of system (14) stabilize at E0∞ (ξ0 < ξ22 ) or E2∞ (ξ0 > ξ22 ), respectively (initial variate ξ0 being the same as Fig. 6(A) and U0 = 0.3). In Fig. 6(D), solutions of system (10) stabilize at E∞ corresponding to Fig. 6(C). Fig. 7(A–D) show the effects of different initial values on the solutions of systems (10) and (14) in region Ω22 . ∗ ∗ (ξ0 > ξ22 ) from different initial values ((0.35, 0.08), Fig. 7(A) demonstrates that solutions of system (14) stabilize at E¯ 25 ∗ (0.4, 0.08), (0.45, 0.08)). Fig. 7(B) demonstrates that solutions of system (10) stabilize at E20 corresponding to Fig. 7(A). In ∗ ∗ Fig. 7(C), solutions of system (14) stabilize at E¯ 0∞ (ξ0 < ξ22 ) or E¯ 2∞ (ξ0 > ξ22 ) from different initial values, i.e., (0.25, 0.1), (0.2, 0.1), (0.25, 0.5), (0.3, 0.5), (0.35, 0.3), (0.4, 0.2), (0.45, 0.2), respectively. In Fig. 7(D), solutions of system (10) stabilize at E∞ corresponding to Fig. 7(C). Similarly, Fig. 8(A–B) show the effects of different initial values on the solutions of systems (10) and (14) in region Ω23 . ∗ ∗ Now, we compare system (10) (or (14)) with (7) (or (13)) (see Table 3). Firstly, both interior equilibria E22 and E23 of ∗ ∗ ∗ system (10) lose their stability compared with E12 and E13 of system (7). Secondly, E20 may be locally stable by different ways in regions Ω21 , Ω22 and Ω26 (see Remark 4), which means both types of individuals may tend to zero, instead of the

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

A

251

B

Fig. 8. The effects of different initial values on the solutions of systems (10) and (14) in region Ω23 . Except b = 3.5, the parameter values are the same as ∗ ∗ Fig. 6. (A): solutions of system (14) stabilize at E¯ 0∞ (ξ0 < ξ22 ) or E¯ 2∞ (ξ0 > ξ22 ) from different initial values ((0.2, 0.1), (0.25, 0.2), (0.3, 0.3), (0.35, 0.1), (0.4, 0.2), (0.45, 0.3)), respectively. (B): Solutions of system (10) stabilize at E∞ corresponding to (A).

A

B

C

D

Fig. 9. The effects of different initial values on some solutions of system (3). The baseline parameter values are the same as Fig. 6. (A)–(B): The solution of system (3) tends to zero from initial value (I0 , U0 ) = (0.0667, 0.1). (C)–(D): The solution of system (3) tends to infinity from initial value (I0 , U0 ) = (0.15, 0.3).

∗ instability of E10 in system (7). Thirdly, infinity equilibrium E∞ exists and it is locally stable for system (10) by different ways (see Remark 5), thus both types of individuals tend to infinity taking the place of finite coexistence in system (7). In addition, ∗ for given initial value (I0 , U0 ), solutions of system (10) are sensitive to initial values, either tending to equilibria E20 or E∞ as shown in Figs. 6–8. Similarly, if the equilibria of system (14) are locally stable, then the corresponding solutions of system (3) are also locally stable. In Fig. 9, we provide some typical solutions of system (3) initiating from different initial values, where the parameter values are fixed in region Ω21 .

252

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

A

C

B

D

Fig. 10. Simulations of systems (7) and (13) in the case τ = 1. The other parameter values are fixed as follows: b = 2, δ = 0.3, T = 2, q = 0.5. (A)–(B): Solution of system (13) stabilizes at (1, 0) from initial value (ξ0 , U0 ) = (0.2, 0.4). (C)–(D): Solution of system (7) stabilizes at (b/(δ T ), 0) from initial value (I0 , U0 ) = (0.1, 0.4) corresponding to (A) and (B).

As we have addressed, successful strategies for control of dengue virus aim to either eradicate its mosquito vectors or totally replace mosquitoes uninfected with Wolbachia with infected ones. However, with a strong density dependent ∗ death rate, when the transmission rate τ ∈ [0, 1), and the initial ratio of infected exceeds the critical value ξ12 , uninfected mosquitoes are partly replaced by infected ones. In the following, we will consider a special case (τ = 1) with different density dependent death rates such that we can control the spread of dengue diseases. 5.2. The transmission is perfect, i.e. τ = 1 ∗ ∗ ∗ ∗ In this special case, for system (13), E¯ 10 and E¯ 14 merge together as (0, 0); E¯ 11 and E¯ 12 merge together as (0, b/(δ T )); and ∗ ∗ ∗ ∗ ∗ E¯ 13 and E¯ 15 merge together as (1, 0). For system (7), E11 and E12 merge together as (0, b/(δ T )), and E13 becomes (b/(δ T ), 0), ∗ denoted by E14 . Similarly, according to the stability analysis of system (13), (0, 0) is unstable, (0, b/(δ T )) is stable, but not ∗ asymptotically stable and (1, 0) is locally stable. So for system (7), E14 is locally stable as shown in Fig. 10. According to the above analyses, there exists a unique local stability of the periodic solution for system (2) (see ∗ Fig. 4(E)–(F)), denoted by (I14 (t ), 0), with ∗ I14 b ∗ , I14 (t ) = ∗ = I14 δ(t − nT ) + 1 δ(b(t − nT ) + T )

t ∈ (nT , (n + 1)T ].

∗ In order to prove the global stability of periodic solution (I14 (t ), 0), we give some basic properties of the following subsystem

 dI (t )   = −δ I (t )2 , 

t ̸= nT , dt + I (nT ) = (1 + b)I (nT ), t = nT ,    + I (0 ) = I (0).

(26)

∗ ∗ (t )| → 0 Lemma 5. Subsystem (26) has a positive periodic solution I14 (t ), and for any solution I (t ) of (26) we have |I (t ) − I14 as t → ∞, where ∗ I14 (t ) =

b

δ(b(t − nT ) + T )

,

t ∈ (nT , (n + 1)T ].

The proof of Lemma 5 is given in Appendix A. According to this lemma, we have the following main results. Theorem 6. If the transmission rate τ = 1, for system (2) with a strong density-dependent death rate function F1 , then the ∗ periodic solution (I14 (t ), 0) of system (2) is globally asymptotically stable in the first quadrant. ∗ The proof of Theorem 6 is shown in Appendix B. The global stability of periodic solution (I14 (t ), 0) ensures that all uninfected mosquitoes are infected by Wolbachia, and totally replaced by infected ones. So in this case, we fulfill the strategy of population replacement such that dengue diseases will be prevented successfully. In the following, we investigate the stability of equilibria for systems (10) and (14) and consequently the periodic solutions of system (3).

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

253

Table 4 Stability analysis of equilibria with transmission rate τ = 1. System (14) Merge Stability System (10) Merge Stability

∗ ¯∗ E¯ 20 , E24 (0, 0)

∗ ¯∗ E¯ 23 , E25 (1, 0) S

E¯ 0∞ , E¯ 1∞ (0, +∞)

E¯ 2∞ (1, +∞)

U

∗ ¯∗ E¯ 21 , E22 (0, M (1)) U/×

U

S

∗ E20 (0, 0) S

∗ ∗ E21 , E22 (0, M (1)) U/×

∗ E23 (M (1), 0) U/×

E∞

(+∞, +∞)

S

Similarly, according to the stability analysis of (10) and (14), when the transmission rate τ = 1, changes of equilibria and the stability of systems (10) and (14) are shown in Table 4. Although one of the eigenvalues of (1, 0) of system (14) cannot be calculated directly, numerical simulations indicate its local stability. In addition, for system (10), both (0, 0) and E∞ reserve local stability of (1, 0) and (1, +∞) of system (14), respectively, while (M (1), 0) cannot preserve local stability of (1, 0) of system (14). From the stability of (0, 0), there exists a solution for system (3) which makes both infected and uninfected population tend to zero. While from the instability of (M (1), 0), there exists an unstable periodic solution for ∗ system (3), denoted by (I24 (t ), 0), with ∗ I24 (t ) = KW



M (1) K

 exp

M (1) − δ(t − nT ) K



,

t ∈ (nT , (n + 1)T ].

So only the strategy of population eradication will be successful in this special case. 6. Discussion and biological implication Dengue fever and the highly lethal dengue hemorrhagic fever affect more than 50 million people annually, primarily in the tropics and subtropical areas. At present, no treatment or vaccine is available for these diseases, so vector control is currently the primary intervention tool. Mosquitoes are the principal vectors for the spread of the virus, and some endosymbiotic Wolbachia can stop the mosquitoes from replicating and thus interrupt the transmission of dengue virus. Therefore, successful strategies for preventing the spread of dengue aim to either eradicate mosquitoes or to totally replace mosquitoes uninfected with Wolbachia by infected ones. An interesting question posed here was how could Wolbachia-induced CI and birth pulses realize these control aims as far as possible with either a strong or a weak density dependent death rate function or both? 6.1. Control strategies with the strong density dependent death function F1 Based on Theorems 1 and 2, there exists an Allee effect for the first variate ξn in system (13) such that low density populations are driven to extinction, but high density populations can survive. Moreover, when the equilibria are locally stable and horizontal ordinates are nonzero, then the strategy of population replacement can be successful. From stability analyses ∗ ∗ of equilibria, interior equilibrium E13 (or E¯ 13 ) is locally stable, but only when transmission from infected mother to offspring ∗ ∗ (τ ) and the rate of zygotic death from CI (q) are confined in region Ω12 such that E13 (or E¯ 13 ) may exist (see Fig. 1). Therefore, to control dengue virus successfully, appropriate Wolbachia bacteria should be selected for trials to ensure the existence of ∗ ∗ ) (here in region Ω12 ). (or E¯ 13 E13 In order to realize the aim of population replacement as far as possible, some distinct possibilities are proposed to make ∗ ∗ ∗ ∗ ∗ ∗ are sufficiently small. On the one hand, ξ13 (ξ12 ) is an increasing and I13 are sufficiently large, while ξ12 and U13 sure that ξ13 (decreasing) function with respect to q or τ , thus q and τ should be increased greatly in region Ω12 (see Fig. 1). Note that when τ is close to one, for any given initial value, the ratio of infected individuals to the total number of ones (ξn ) and the ∗ ∗ are zero simultaneously, density of infected individuals (In ) will stabilize at one and bτ /(δ T ), respectively, and ξ12 and U13 regardless of the value of q. In this case, all uninfected individuals are infected and replaced by infected ones. So for the strong density dependence death function F1 , we can fulfill the strategy of population replacement such that Wolbachia could totally invade a susceptible population. However, in practice, it is difficult to achieve perfect transmission of Wolbachia. Note that the greater the density of infected individuals is, the more it favors dengue control. So in the following, we fix all parameters except for the key parameters which describe the natural birth rate (b), the probability of transmission from infected mother to offspring (τ ), the natural death rate (δ) and the birth pulses period (T ), and choose different rates of zygotic death for CI (here q ∈ [0, 1]). ∗ Then different parameter sets on the quantity of infected mosquitoes I13 are shown in Fig. 11. Because some domains of ∗ definition of parameter q lead to parameter space q and τ lying in region Ω11 , I13 does not exist at all and there are no ∗ curves in these domains (see Fig. 11). Moreover, the density of infected individuals I13 is increasing with the increase of the rate of zygotic death from CI (q). Fig. 11(B) indicates the effects of the probability of transmission from an infected mother to offspring (τ ), and the rate ∗ of zygotic death for CI (q) on the I13 , and the results indicate that the larger the probability of transmission from infected mother to offspring, the larger the quantity of infected mosquitoes which follows. However, as shown in Fig. 11(C), we found

254

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

A

B

C

D

Fig. 11. The effects of the probability of zygotic death for CI (q), natural birth rate (b), the probability of transmission from infected mother to ∗ offspring (τ ), natural death rate (δ) and birth pulse period (T ) on the density of infected individuals I13 . The baseline parameter values are as follows: ∗ τ = 0.95, b = 2, δ = 0.2, T = 4. (A): The effect of natural birth rate b on I13 . (B): The effect of the probability of transmission from infected mother to ∗ ∗ ∗ offspring τ on I13 . (C): The effect of natural death rate δ on I13 . (D): The effect of natural birth pulse period T on I13 .

∗ an opposite impact of the natural death rate (δ) and the rate of zygotic death from CI (q) on the I13 , i.e. the smaller the natural death rate, the larger the quantity of infected mosquitoes which follows. Similar effects of natural birth rate (b), natural death rate (δ, ) birth pulse period (T ) and the rate of zygotic death for CI ∗ (q) on the density of infected individuals I13 are shown in Figs. 11(A)–(D) and 12. The simulations show that these parameters are all crucial to dengue disease control. They also confirm that the larger the natural birth rate, the shorter the birth pulse period, and the larger rate of zygotic death from CI make dengue control more effective. In practice, when the total population size is low in some areas, a low density of infected individuals can also make the ∗ ratio of infected to total individuals reach above threshold ξ12 . In addition, the immigration of infection from extraneous sources may counteract the Allee effect, such that the ratio of infection exceeds the threshold.

6.2. Control strategies with the weak density dependent death function F2 According to Theorems 3 and 4, with a weak density dependent death function F2 , an Allee effect also exists in system ∗ (10) and two attractors E20 and E∞ may coexist. In this case, it is obvious that strategies for population extinction such that two types of population are reduced as far as possible should be adopted. Based on Tables 1 and 2, equilibrium E∞ is always ∗ locally stable, while E20 is locally stable only in regions Ω21 , Ω22 and Ω26 . Unfortunately, solutions of system (10) may be very sensitive to initial values. Numerical simulations indicate that under the same initial ratio of infected to total individuals, smaller initial U0 makes both infected and uninfected populations become extinct, while a larger U0 makes both infected and uninfected populations tend to infinity for system (10) (see Figs. 6 and 7). Therefore, when the parameter space lies in regions Ω21 , Ω22 and Ω26 , and the initial density of the uninfecteds is low enough, the strategies of population eradication may be realized in spite of the ratio of infection being high, which would prevent the spread of dengue virus. Based on the above discussion, it is worth pointing out that different levels of density dependent death functions F1 and F2 may result in strategies of population replacement and population eradication for dengue virus control respectively even if some conditions, such as the selection of Wolbachia bacteria, effect of CI, release methods and rates of infected mosquitoes and so on are the same. So in practice, it is critical to learn of the mosquito death rate at first so that the appropriate strategies for dengue control can be devised. Individuals of different ages possess different fertility and mortality rates, thus females and males infected with Wolbachia or not may have different fates when they mate due to the CI mechanism and matrilineal inheritance and other sources of

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

A

B

C

D

255

Fig. 12. The effects of natural birth rate (b), the probability of zygotic death for CI (q), the probability of transmission from infected mother to offspring ∗ (τ ), the natural death rate (δ) and the birth pulse period (T ) on the density of infected individuals I13 . The baseline parameter values are as follows: ∗ τ = 0.95, q = 0.8, δ = 0.2, T = 4. (A): The effect of the probability of zygotic death for CI q on I13 . (B): The effect of the probability of transmission from ∗ ∗ ∗ infected mother to offspring τ on I13 . (C): The effect of the natural death rate δ on I13 . (D): The effect of the natural birth pulse period T on I13 .

variation. The latter include different release quantities of infected and uninfected populations and ratios of infected insects to the total population. In future work, we intend to expand our model to include age and gender structures, so that we can analyze how different ages and sexes affect the dynamics of the system and provide improved strategies for dengue control. Acknowledgments This work is supported by the National Natural Science Foundation of China (NSFCs, 11171199, 11471201) and by the Fundamental Research Funds for the Central Universities (GK201305010, GK201401004, and S2014YB01). Appendix A. The proof of Lemma 5 The solution of system (26) at any impulsive interval (nT , (n + 1)T ] is given by I (t ) =

I (nT + )

δ(t − nT )I (nT + ) + 1

,

then 1 I (t )

= δ(t − nT ) +

1 I (nT + )

.

So by induction we have 1 I (t )

= δ(t − nT ) +

δT b

 1−



1

(1 + b)

n

+

1

(1 + b)n I (0)

,

nT < t ≤ (n + 1)T .

(A.1)

256

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

According to the analytic expression of periodic solution I4∗ (t ), we get 1 I4 (t ) ∗

= δ(t − nT ) +

δT b

,

nT < t ≤ (n + 1)T .

(A.2)

It follows from the two Eqs. (A.1) and (A.2) that

       ∗   1    1  1 1    =  I (t ) − I4 (t )  =  δ T − −  I (t ) I ∗ (t )   I (t )I ∗ (t )   b (1 + b)n I (0) (1 + b)n  4 4 → 0 as t → ∞.            δT 1 1 1    =  δ(t − nT ) + δ T δ( t − nT ) + 1 − +  I (t )I ∗ (t )   n n b b (1 + b) (1 + b) I (0)  4 → ∞ as t → ∞.

(A.3)

(A.4)

Based on Eqs. (A.3) and (A.4),

   1  1   |I (t ) − I4 (t )| =  − ∗  I (t )I4∗ (t ) → 0 I (t ) I4 ( t ) ∗

as t → ∞. This completes the proof of Lemma 5. Appendix B. The proof of Theorem 6 ∗ According to Section 5.2, equilibrium E14 of system (7) is locally stable, so the corresponding periodic solution (I4∗ (t ), 0) of system (2) is locally stable. Thus to prove the global stability, we only need to show its global attractivity. It follows from the second and fourth equations of system (2) that we get

 dU (t )   ≤ −δ I (t )U (t ),

t ̸= nT ,

dt

I (nT )U (nT )  U (nT + ) ≤ (1 + b)U (nT ) − bq , I (nT ) + U (nT )

(B.1)

t = nT .

Now we consider the following impulsive differential equation

 dy(t )   = −δ I (t )y(t ),  

t ̸= nT ,

dt

I (nT )y(nT ) y(nT + ) = (1 + b)y(nT ) − bq ,   I ( nT ) + y(nT )   + y(0 ) = U (0).

(B.2)

t = nT ,

In systems (B.1) and (B.2), we assume the items I (t ) > 0, I (nT ) > 0. According to the comparison theorem on impulsive differential equations we get U ((n + 1)T ) ≤ y(nT + ) exp

(n+1)T



−δ I (t )dt



nT



= y(nT ) 1 + b −

bqI (nT )



I (nT ) + y(nT )

(n+1)T

 exp

 −δ I (t )dt ,

nT

= y(0 )η , +

n 1

(B.3)

where



η1 , 1 + b −

bqI (nT ) I (nT ) + y(nT )





(n+1)T

exp



−δ I (t )dt .

(B.4)

nT

According to the first equation of (6), we have I (t ) =

I (nT + )

δ(t − nT )(U (nT + ) + I (nT + )) + 1

.

For system (7), it is impractical to consider the initial value (I0 , U0 ) = (0, 0), any (I0 , U0 ) ∈ R2+ , the solution will be I (nT + ) → b/(δ T ) and U (nT + ) → 0 as n → ∞, then bqI (nT )

η1 =

1 + b − I (nT )+y(nT ) < 1. δ T (U (nT + ) + I (nT + )) + 1

Therefore, U (nT ) ≤ U (0+ )η1n , U (nT ) → 0 as n → ∞, then U (t ) → 0 as t → ∞.

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

257

Next, we prove that I (t ) → I4∗ (t ) as t → ∞. For any ϵ1 > 0, there must exist a t1 such that 0 < U (t ) < ϵ1 /2 for t > t1 . Without loss of generality, we may assume that 0 < U (t ) < ϵ1 /2 hold true for all t > 0, then we have

 ϵ1 2 dI (t ) −δ I (t ) + ≤ ≤ −δ I (t )2 . 2

dt

The right hand inequality follows from the first and third equations of system (2), so we consider the following impulsive differential equation:

 dx(t )   = −δ x(t )2 ,

t ̸= nT , dt +  x(nT+ ) = (1 + b)x(nT ), t = nT , x(0 ) = I (0).

(B.5)

From Lemma 5 and the comparison theorem on impulsive differential equations, we have I (t ) < x(t ) and x(t ) → I4∗ (t ) as t → ∞. Hence I (t ) < x(t ) < I4∗ (t ) + ϵ

(B.6)

for all large t. For simplification we may assume (B.6) holds for all t > 0. For the left hand inequality, we consider the following impulsive differential equation

  dz (t ) ϵ1 2   = −δ z (t ) + , dt + 2  z (nT+ ) = (1 + b)z (nT ), z (0 ) = I (0).

t ̸= nT , (B.7)

t = nT ,

By using the same methods as those for system (B.5), we can easily prove that system (B.7) has a globally stable periodic solution, denoted by z ∗ (t ) and z ∗ (t ) =

z (nT + ) +

δ z (nT + ) + 

ϵ1  2

ϵ1 2

(t − nT ) + 1

−

ϵ1 2

,

t ∈ (nT , (n + 1)T ].

Therefore, for any ϵ2 > 0, there exists a t2 > 0 such that z ∗ (t ) − ϵ2 < I (t ) < I4∗ (t ) + ϵ2 for t > t2 . Let ϵ1 → 0, then we have I4∗ (t ) − ϵ2 < I (t ) < I4∗ (t ) + ϵ2 for t > t2 , which indicates that I (t ) → I4∗ (t ) as t → ∞. This completes the proof of Theorem 6. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

[14] [15] [16] [17] [18] [19] [20] [21]

D.J. Gubler, G. Kuno, Dengue and Dengue Hemorrhagic Fever, Cab International, New York, 1997, p. 478. J.L. Kyle, E. Harris, Global spread and persistence of dengue, Annu. Rev. Microbiol. 62 (2008) 71–92. J.G. Rigau-Pérez, Severe dengue: the need for new case definitions, Lancet Infect. Dis. 6 (2006) 297–302. J. Hemingway, H. Ranson, Insecticide resistance in insect vectors of human disease, Annu. Rev. Entomol. 45 (2000) 371–391. R.F. Qi, L. Zhang, C.W. Chi, Biological characteristics of dengue virus and potential targets for drug design, Acta. Biochim. Biophys. Sin. 40 (2008) 91–101. B.J. Beaty, Genetic manipulation of vectors: a potential novel approach for control of vector-borne diseases, Proc. Natl. Acad. Sci. 97 (2000) 10295–10297. M.M. Pettigrew, S.L. O’Neill, Control of vectorborne disease by genetic manipulation of insect vectors: technological requirements and research priorities, Aust. J. Entomol. 36 (1997) 309–317. J.H. Werren, L. Baldo, M.E. Clark, Wolbachia: master manipulators of invertebrate biology, Nat. Rev. Microbiol. 6 (2008) 741–751. P. Kittayapong, K.J. Baisley, V. Baimai, S.L. ONeill, Distribution and diversity of Wolbachia infections in Southeast Asian mosquitoes (Diptera: Culicidae), J. Med. Entomol. 37 (2000) 340–345. H. Laven, Crossing experiments with Culex strains, Evolution 5 (1951) 370–375. J.H. Yen, A.R. Barr, A new hypothesis of the cause of cytoplasmic incompatibility in Culex pipiens, Nature 232 (1971) 657–658. T. Walker, P.H. Johnson, L.A. Moreira, I. Iturbe-Ormaetxe, F.D. Frentiu, C.J. McMeniman, Y.S. Leong, Y. Dong, J. Axford, P. Kriesner, A.L. Lloyd, S.A. Ritchie, S.L. O’Neill, A.A. Hoffmann, The wMel Wolbachia strain blocks dengue and invades caged Aedes aegypti populations, Nature 476 (2011) 450–453. A.A. Hoffmann, B.L. Montgomery, J. Popovici, I. Iturbe-Ormaetxe, P.H. Johnson, F. Muzzi, M. Greenfield, M. Durkan, Y.S. Leong, Y. Dong, H. Cook, J. Axford, A.G. Callahan, N. Kenny, C. Omodei, E.A. McGraw, P.A. Ryan, S.A. Ritchie, M. Turell, S.L. O’Neill, Successful establishment of Wolbachia in Aedes populations to suppress dengue transmission, Nature 476 (2011) 454–457. J.L. Crainey, M.D. Wilson, R.J. Post, Phylogenetically distinct Wolbachia gene and pseudogene sequences obtained from the African onchocerciasis vector Simulium squamosum, Int. J. Parasitol. 40 (2010) 5690. J.A. Nirgianaki, G.K. Banks, D.R. Frohlich, Z. Veneti, H.R. Braig, T.A. Miller, I.D. Bedford, P.G. Markham, C. Savakis, K. Bourtzis, Wolbachia infections of the Whitefly Bemisia tabaci, Curr. Microbiol. 47 (2003) 93–101. Anon, Oh, New Delhi; Oh, Geneva, Nature 256 (1975) 355–357. H. Laven, Eradication of Culex pipiens fatigans through cytoplasmic incompatibility, Nature 216 (1967) 383–384. R. Pal, WHO, ICMR programme of genetic control of mosquitoes in India, in: R. Pal, M.J. Whitten (Eds.), The Use of Genetics in Insect Control, Elsevier, Amsterdam, 1974, pp. 73–95. M. Turelli, Cytoplasmic incompatibility in population with overlapping generations, Evolution 64 (2010) 232–241. H.L. Yeap, P. Mee, T. Walker, A.R. Weeks, S.L. O’Neill, P. Johnson, S.A. Ritchie, K.M. Richardson, C. Doig, N.M. Endersby, A.A. Hoffmann, Dynamics of the ‘popcorn’ Wolbachia infection in Aedes aegypti in an outbred background, Genetics 187 (2011) 583–595. P.G. Schofield, Spatially explicit models of Turelli–Hoffmann Wolbachia invasive wave fronts, J. Theoret. Biol. 215 (2002) 121–131.

258

X. Zhang et al. / Nonlinear Analysis: Real World Applications 22 (2015) 236–258

[22] M.J. Keeling, F.M. Jiggins, J.M. Read, The invasion and coexistence of competing Wolbachia strains, Heredity 91 (2003) 382–388. [23] J.L. Rasgon, T.W. Scott, Impact of population age structure on Wolbachia transgene driver efficacy: ecologically complex factors and release of genetically modified mosquitoes, Insect Biochem. Mol. Biol. 34 (2004) 707–713. [24] J. Engelstädter, A. Telschow, P. Hammerstein, Infection dynamics of different Wolbachia-types within one host population, J. Theoret. Biol. 231 (2004) 345–355. [25] E. Vautrin, S. Charles, S. Genieys, F. Vavre, Evolution and invasion dynamics of multiple infections with Wolbachia investigated using matrix based models, J. Theoret. Biol. 245 (2007) 197–209. [26] J.Z. Farkas, P. Hinow, Structured and unstructured continuous models for Wolbachia infections, Bull. Math. Biol. 72 (2010) 2067–2088. [27] J.G. Schraiber, A.N. Kaczmarczyk, R. Kwok, et al., Constraints on the use of lifespan-shortening Wolbachia to control dengue fever, J. Theoret. Biol. 297 (2012) 26–32. [28] H. Caswell, Matrix Population Models: Construction, Analysis and Interpretation, Sinauer Associates, Sunderland, 1989. [29] J.Y. Wu, Z.R. Lun, A.A. James, X.G. Chen, Review: dengue fever in mainland China, AM. J. Trop. Med. Hyg. 83 (2010) 664–671. [30] T.U. Ahmed, G.M.S. Rahman, K. Bashar, M. Shamsuzzaman, S. Samajpati, S. Sultana, M.I. Hossain, N.N. Banu, M.S. Rahman, Seasonal prevalence of dengue vector mosquitoes in Dhaka city, Bangladesh, Bangladesh J. Zool. 35 (2007) 205–212. [31] A. Chakravarti, R. Kumaria, Eco-epidemiological analysis of dengue infection during an outbreak of dengue fever, India, Virol. J. 2 (2005) 32. [32] R.S. Sharma, S.M. Kaul, J. Sokhay, Seasonal fluctuations of dengue fever vector, Aedes aegypti (Diptera: Culicidae) in Delhi, India, Southeast Asian J. Trop. Med. Public 36 (2005) 186–190. [33] S. Kaup, J. Sankarankutty, Seroprevalence and seasonal trend of dengue virus infection at a teaching hospital in Tumkur, India, Sch. J. App. Med. Sci. 2 (2014) 922–926. [34] L.M. Bartley, C.A. Donnelly, G.P. Garnett, The seasonal pattern of dengue in endemic areas: mathematical models of mechanisms, T. Roy. Soc. Trop. Med. H. 96 (2002) 387–397. [35] S. Wongkoon, M. Jaroensutasinee, K. Jaroensutasinee, Distribution, seasonal variation & dengue transmission prediction in Sisaket, Thailand, Indian J. Med. Res. 138 (2013) 347–353. [36] N. Becker, D. Petric, M. Zgomba, C. Boase, M. Madon, C. Dahl, A. Kaiser, Mosquitoes and Their Control, Springer Science & Business Media, Berlin, 2010, pp. 9–24. [37] Biological Notes on Mosquitoes, http://www.megacatch.com/mosquitolifecycle.html. [38] D.D. Bainov, P.S. Simeonov, Systems with Impulsive Effect: Stability, Theory and Applications, Kohn Wiley, New York, 1989. [39] A. Lakmeche, O. Arino, Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherapeutic treatment, Contin. Discrete Impulsive Syst. 7 (2000) 165–287. [40] M.G. Roberts, R.R. Kao, The dynamics of an infectious disease in a population with birth pulses, Math. Biosci. 149 (1998) 23–36. [41] S.Y. Tang, R.A. Cheke, State-dependent impulsive models of integrated pest management (IPM) strategies and their dynamic consequences, J. Math. Biol. 50 (2005) 257–292. [42] S.Y. Tang, Y.N. Xiao, L.S Chen, R.A. Cheke, Integrated pest management models and their dynamical behaviour, B. Math. Biol. 67 (2005) 115–135. [43] H.C. Pedersen, H. Steen, L. Kastdalen, H. Brøseth, R.A. Ims, W. Svendsen, N.G. Yoccoz, Weak compensation of harvest despite strong density-dependent growth in willow ptarmigan, Proc. R. Soc. 271 (2004) 381–385. [44] S. Einum, Salmonid population dynamics: stability under weak density dependence? Oikos 110 (2005) 630–633.