# ModelingWolbachia Diffusion in Mosquito Populations by Discrete Competition Model

Authors: Yijie Li; Zhiming Guo; Yanyuan Xing
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

Dengue fever is caused by dengue virus and transmitted byAedes mosquitoes. A promising avenue to control this disease is to infect the wild Aedes population with the bacterium Wolbachia driven by cytoplasmic incompatibility (CI). To study the invasion of Wolbachia into wild mosquito population, we formulate a discrete competition model and analyze the competition between released mosquitoes and wild mosquitoes. We show the global asymptotic properties of the trivial equilibrium, boundary equilibrium, and positive equilibrium and give the conditions for the successful invasion of Wolbachia. Finally, we verify our findings by numerical simulations.


---

## Body

## 1. Introduction

Dengue fever is a mosquito-borne infectious disease caused by dengue virus, mainly transmitted byAedes aegypti [1]. More than 2.5 billion people in more than 100 countries are currently at risk from this disease, making mosquito control a major public health priority. Conventional measures, such as mass spraying of insecticides, are limited and unsafe, and worse, no effective vaccine that against dengue has yet been developed [2].An innovative and effective method of mosquito control is using specific endosymbiotic bacteriumWolbachia to block the transmission of dengue fever. With infection, Wolbachia can reduce the mosquito’s dengue transmission potential. Besides, Wolbachia often induces cytoplasmic incompatibility (CI) which leads to early embryonic death when Wolbachia-infected males mate with uninfected females, so the release of infected males helps to reduce mosquito populations [3]. This bacterium infects many arthropods, including some mosquito species in nature, which was first identified in 1924 [4]. However, Wolbachia is not present in some important mosquito vectors. In 2005, Xi et al. established stable mosquito strains carrying Wolbachia for the first time by injecting Wolbachia into Aedes aegypti [5, 6], which is the basis of using Wolbachia to control mosquito-borne diseases. The strategies of using Wolbachia for mosquito vector control mainly include population suppression and population replacement. The former only releases infected males, and produces CI effect after mating with female mosquitoes in the field to suppress the number of mosquito vectors in the field [7, 8]; the latter releases both male and female mosquitoes infected with Wolbachia using vertical maternal transmission of Wolbachia and certain advantages to let Wolbachia spread in the mosquito population to stop spreading disease [9].In recent years, the spreading dynamics ofWolbachia has become a hot research topic in academia. In 2014, Zheng et al. have established a delay differential equation model to study Wolbachia infection dynamics [10]. They gave a precise threshold for the infection frequency, and their numerical simulation is well fitting with the experimental data; then, Huang et al. proposed the corresponding reaction diffusion equations with homogeneous Neumann boundary condition and obtained similar results [11, 12]. In 2019, Zheng et al. have carried out a mathematical model to study combining incompatible and sterile insect techniques (IIT-SIT) to eliminate the population of Aedes mosquitoes in the wild, and they proved the feasibility of regional control of mosquito vector population by combining IIT-SIT through field experiments [13]. Besides, some interesting mathematical models have also been developed, such as the time-delay continuous model [8, 14–16], the discrete model [17, 18], and the stochastic systems of differential equations [19].Motivated by those valuable research studies and since the experimental data of the detection of mosquito population in the field are discrete, we managed to establish a discrete model and make it as reasonable as possible. In this paper, a discrete model was established to study the spread ofWolbachia in mosquito populations from the perspective of competition, which is rarely mentioned in the existing literature.The most frequently used equation for studying the discrete model is the Beverton-Holt equation:(1)xt+1=b11+c11xtxt,where b>0 denotes the birth rate and c11>0 denotes the intraspecific competition coefficient. If b>1, then all solutions of the equation with x0>0 tend monotonically to the equilibrium x=b−1/c11. If xt is considered as the number of a species at time t, it is a discrete model similar to the classical logistic continuous model [20]. If we consider the discrete model of competition between two species, the Beverton-Holt equation needs to be modified, which leads to the Leslie-Gower discrete competition model:(2)xt+1=b111+c11xt+c12ytxt,yt+1=b211+c21xt+c22ytyt,where b1,b2>0 denote the birth rate of xt,yt, respectively. c11,c22>0 denote the intraspecific competition coefficient and c12,c21>0 denote the interspecific competition coefficient, respectively. Model (2) has rich dynamic properties and satisfies the competitive exclusion principle [21, 22]. Based on (2), we establish a discrete competition model describing Wolbachia spread in mosquito population. We assume that xt represents the number of mosquito population infected with Wolbachia at t and yt represents the number of uninfected mosquito population at t, and their birth rates are b1 and b2, respectively. Since the mosquitoes are overlapping generation populations, we introduce d1,d2 as the mortality rates of xt and yt, respectively. Then, their survival rates are, respectively, 1−d1 and 1−d2. Because Wolbachia often induces fitness cost [10], and it is reasonable to assume that 0<d2<d1<1. Furthermore, considering that both xt and yt denote the numbers of Aedes aegypti, for simplicity, we ignore the difference between intraspecific competition and interspecific competition and focus on the effect of competitions on their birth rates, so we assume c11=c12=α and c21=c22=β denote the competition coefficient within (or between) species, respectively. Finally, considering the infection Wolbachia CI mechanism of mosquito population and vertical maternal transmission of Wolbachia [23], the birth rate of yt is replaced by b2yt/xt+yt; hence, the discrete model of competition between two species is as follows:(3)xt+1=b1xt1+αxt+yt+1−d1xt,yt+1=b2yt1+βxt+ytytxt+yt+1−d2yt,where all the coefficients are positive, and we mainly discuss this model.The structure of this paper is arranged as follows. In Section2, firstly, we prove the nonnegativity and boundedness of the solution of model (3) and obtain the positive invariant set. Secondly, we give four equilibrium points, including the trivial equilibrium E0, boundary equilibrium E1,E2, and positive equilibrium E3 as well as their existence conditions, and we also introduce the global attraction conditions for each equilibrium by stability analysis. In Section 3, some numerical simulations are presented to illustrate the behavior of the model and we interpret out main results biologically. Finally, the conclusion and discussion will be given in Section 4.


## 2. Model Analysis

We first consider the positivity and boundedness of the solution of model (3), positivity is clear, and the boundedness can be obtained as follows. From the first equation of model (3), we find(4)xt+1≤b1α+1−d1xt.Letσ=1−d1, then we have xt+1−σxt≤b1/α. Furthermore, we can obtain a series of inequalities(5)xt−σxt−1≤b1α,σxt−1−σ2xt−2≤σb1α,σ2xt−2−σ3xt−3≤σ2b1α,…σt−1x1−σtx0≤σt−1b1α.Add them up, and we get(6)xt−σtx0≤b1α1+σ+σ2+⋯+σt−1.Asσ∈0,1, it is easy to obtain the following inequality for any given initial value x0:(7)limsupt⟶∞xt≤b1α11−σ=b1αd1.Similarly, foryt, we have(8)limsupt⟶∞yt≤b2βd2.Therefore,(9)Γ=x,y∈R20≤x≤b1αd1,0≤y≤b2βd2is a positive invariant.Next, we analyze the existence and stability of the equilibria to (3). Note that (3) is not defined at the origin. In order to study the dynamical behavior of solutions to (3) near the origin, we define an auxiliary function:(10)fx,y=b2y1+βx+yyx+y+1−d2y,x,y≠0,0,0,x,y=0,0.Thus, (3) can be modified as(11)xt+1=b1xt1+αxt+yt+1−d1xt,yt+1=fxt,yt.Thus,E00,0 can be regarded as an equilibrium of (3). The other equilibria of (3) satisfy the following algebraic equations:(12)x=b1x1+αx+y+1−d1x,y=b2y1+βx+yyx+y+1−d2y.There are two nonnegative boundary equilibria of model (12), E11/αb1/d1−1,0 and E20,1/βb2/d2−1, provided b1>d1 and b2>d2. Besides, a positive equilibrium(13)E31αb1d1−11−d2b21+βαb1d1−1,d2αb2b1d1−11+βαb1d1−1exists if b1>d1,b2>d2, and 1/βb2/d2−1>1/αb1/d1−1.Theorem 1.
For model (3), we have the following conclusions.(a)
Ifb1≤d1,b2≤d2, then E0 is globally asymptotically stable.(b)
Ifb1>d1,b2≤d2, then E0 is unstable and E1 is globally asymptotically stable.(c)
Ifb1≤d1,b2>d2, then E0 is unstable and E2 is globally asymptotically stable.Proof.
The Jacobian of model (3) is(14)J=b11+αy1+αx+y2+1−d1−αb1x1+αx+y2−b2y21+2βx+y1+βx+y2x+y2b2y2x+y+2b2βxyx+y1+βx+y2x+y2+1−d2.
The Jacobians evaluated atE1andE2 read as(15)J1=d12b1+1−d1d12b1−d101−d2,(16)J2=b11+α/βb2/d2−1+1−d10d22b2−2d2d22b2+1−d2,respectively. Their eigenvalues appear along the diagonals. Recall that an equilibrium is globally asymptotically stable if it is locally asymptotically stable and all orbits tend to the equilibrium as t⟶∞.(a)
Ifb1≤d1,b2≤d2, there is no nonnegative equilibrium except for E0. When b1<d1,b2<d2, from (3), we have(17)0≤xt+1≤b1−d1+1xt,0≤yt+1≤b2−d2+1yt,the system is monotonically decreasing, and ift⟶∞, then xt⟶0,yt⟶0, so that E0 is globally asymptotically stable. Furthermore, when b1=d1,b2=d2, from (3), we have(18)xt+1=b11+αxt+yt+1−d1xt,yt+1=b21+βxt+ytytxt+yt+1−d2yt.Obviously,xt+1=xt, yt+1=yt if and only if xt=yt=0 for all t≥0. Otherwise, we claim that 0≤xt+1<xt and 0≤yt+1<yt for all t≥0. It follows that limt⟶+∞xt and limt⟶+∞yt exist and must be zero. Therefore, E0 is globally asymptotically stable. For the remaining cases, i.e., b1<d1,b2=d2 and b1=d1,b2<d2, the conclusion can be verified similarly.(b)
Ifb1>d1,b2≤d2, there are two equilibria, E0 and E1. When yt=0, from the first equation of model (3), we have(19)xt+1=b11+αxt+1−d1xt.Clearly,b1/1+αxt+1−d1>1 if and only if xt<1/αb1/d1−1. It follows that the solution emanating from the X-axis with 0<x0<1/αb1/d1−1 will be away from E0 along the X-axis, so E0 is unstable. The eigenvalues of J1 belong to the interval 0,1, and hence E1 is locally asymptotically stable. In addition, according to a theorem due to Liu and Elaydi [24], all bounded orbits will eventually converge to an equilibrium; thus, all orbits will converge to E1. Therefore, E1 is globally asymptotically stable.(c)
Ifb1≤d1,b2>d2, there are two equilibria, E0 and E2. In this case, we can prove the results in a similar way to case (b), so we omit it.
From now on, we assume thatb1>d1,b2>d2. By a similar argument, we know that E0 is unstable. Therefore, in the sequel, we will mainly study the boundary equilibria E1, E2 and positive equilibrium E3.
We claim thatE3 exists if and only if(20)b1>d1,b2>d2,1βb2d2−1>1αb1d1−1.
We analyze it graphically, from the first equation of (12), we can obtain the perpendicular isoclines(21)y=−x+1αb1d1−1,denoted by L1, which is a segment of lines with slope −1 and intercept 1/αb1/d1−1. From the second equation of (12), we can obtain the horizontal isoclines(22)βx2+2βxy+βy2+x+1−b2/d2y=0,denoted by L2, which is a parabola, and it goes through three fixed points 0,0,−1/β,0, and 0,1/βb2/d2−1. We distinguish two cases depending on the position of L1 and L2 as shown in Figure 1. The inequalities satisfied by the coefficients that correspond to these cases can be listed as follows:Case A(23)b1>d1,b2>d2,1βb2d2−1<1αb1d1−1.In this case, there exist three equilibriaE0,E1, and E2, but E3 does not exist.Case B(24)b1>d1,b2>d2,1βb2d2−1>1αb1d1−1.
In this case, there exist four equilibriaE0,E1,E2, and E3, the positive equilibrium E3 is the intersection of L1 and L2.
For the above two cases, a feasible biological interpretation of the conditions is that if1/βb2/d2−1<1/αb1/d1−1, then the environment is more favorable for infected mosquitoes, and we say that Wolbachia infection has a fitness benefit. On the contrary, if 1/βb2/d2−1>1/αb1/d1−1, then the environment is less favorable for infected mosquitoes, and we claim that the infection has a fitness cost [10]. The stability analysis of the equilibria in both cases will be shown in the following theorem.Figure 1
The two possible cases ofL1 and L2. (a) Case A. (b) Case B.
(a)
(b)Theorem 2.
Supposeb1>d1,b2>d2. Then, we have the following conclusions.(a)
In the fitness benefit case1/βb2/d2−1<1/αb1/d1−1, E2 is a saddle and E1 is globally asymptotically stable.(b)
In the fitness cost case1/βb2/d2−1>1/αb1/d1−1, both E1 and E2 are locally asymptotically stable while E3 is a saddle.Proof.
(a)
Case A: from J2 (see (16)), we know that E2 is a saddle and is unstable and the stable manifold is the Y-axis. From the stable manifold theory and Hartman-Grobman theorems [25, 26], it follows that no solution can approach E2. By Liu and Elaydi’s theorem, all solutions must approach E1. Besides, the eigenvalues of J1 (see (15)) are positive and less than 1, so E1 is locally asymptotically stable. Therefore, E1 is globally asymptotically stable.
(b)
Case B: all of the eigenvalues of J1 (see (15)) and J2 (see (16)) are positive and less than 1, so both E1 and E2 are locally asymptotically stable.
The stability ofE3 is discussed below. The Jacobian at E3 takes the form(25)J3=d121+d2/b2b1/d1−11+β/αb1/d1−1b1+1−d1d12b1−d11−d2b21+βαb1d1−1−d221+2β/αb1/d1−1b2−d221+2β/αb1/d1−1b2+1+d2.
According to the Jury criteria [27], the equilibrium E3 is locally asymptotically stable if the Jury condition(26)trJ3<1+detJ3<2,holds, where trJ3 is the trace and detJ3 is the determinant of J3. If at least one of these inequalities is reversed, then the equilibrium is unstable. The inequalities are equivalent to the following three inequalities.(1)
1+detJ3<2,(2)
−1−detJ3<trJ3,(3)
trJ3<1+detJ3.
By a direct calculation, it shows that(27)detJ3=d121+d21+d2/b2b1/d1−11+β/αb1/d1−1b1−d221+2β/αb1/d1−1b2+1−d11+d2,trJ3=d121+d2/b2b1/d1−11+β/αb1/d1−1b1−d221+2β/αb1/d1−1b2+2−d1+d2.
Using these formulas, we find that inequality3 is equivalent to(28)d12d21+d2/b2b1/d1−11+β/αb1/d1−1b1−d1d2>0,and the reverse inequality holds in Case B, so that E3 is unstable.
In what follows, we prove thatE3 is a saddle.
Inequality2 is equivalent to(29)d12+d2d1b11+d2b2b1d1−11+βαb1d1−1−1+2d21−d2b21+2βαb1d1−1+4>0.
We will prove that it holds inCase B.
Let’s analyze the two complicated parts of the above equation. Obviously, the inequality(30)−1<d1b11+d2b2b1d1−11+βαb1d1−1−1<0,holds in Case B. Also, we have(i)
If1/βb2/d2−1<2/βb1/d1−1, then 1−d2/b21+2β/αb1/d1−1<0.(ii)
If1/βb2/d2−1=2/αb1/d1−1, then 1−d2/b21+2β/αb1/d1−1=0.(iii)
If1/βb2/d2−1>2/αb1/d1−1, then 1−d2/b21+2β/αb1/d1−1>0.
If1/βb2/d2−1>2/αb1/d1−1 or 1/βb2/d2−1=2/αb1/d1−1, then(31)d12+d2d1b11+d2b2b1d1−11+βαb1d1−1−1+2d21−d2b21+2βαb1d1−1+4≥d12+d2d1b11+d2b2b1d1−11+βαb1d1−1−1+4>−d12+d2+4>0.
Consequently, we only need to prove−1−detJ3<trJ3 holds when 1/αb1/d1−1<1/βb2/d2−1<2/αb1/d1−1. By using this condition, we have(32)d12+d2d1b11+d2b2b1d1−11+βαb1d1−1−1+2d21−d2b21+2βαb1d1−1+4≥d12+d2d1b11+d2b2b1d1−11+β2βb2d2−1−1+2d21−d2b21+2ββb2d2−1+4≥d12+d2d1b11+d2b2b1d1−112b2d2+1−1+2d21−d2b21+2b2d2−1+4≥d12+d212d1b1b1d1+1−1+2d21−d2b22b2d2−1+4>−12d12+d2−2d2+4>0.
The last inequality employs the condition0<d2<d1<1.
Now, we consider the characteristic quadratic polynomial ofJ3, namely,(33)pλ=λ2−trJ3λ+detJ3.
Clearly,p1=1−trJ3+detJ3<0 because the inequality 3 is reversed in Case B. Moreover, p−1=1+trJ3+detJ3>0 because the inequality 2 holds, and p+∞=+∞. Thus, pλ has a positive root λ>1 and a root between −1 and 1. Therefore, E3 is a saddle.
The first part of Theorem 2 predicts that if the environment favorsWolbachia-infected populations, then successful invasion is ensured for any initial release size [28]. Biologically, the other case is more feasible when Wolbachia-infected populations have a fitness cost. In this case, model (3) possesses a saddle point E3, whose stable curves separate the attraction regions of E1 and E2. We illustrate our results by numerical simulations in the next section.


## 3. Numerical Results

In order to show the dynamic stability of the model (3), we use MATLAB technical computing software for numerical simulations. The simulations for a and b in Theorem 2 is mainly conducted by using different parameters. Some results are shown as follows.From Figure2 we see that in Case A, namely, the fitness benefit case, E1 is globally asymptotically stable. The result of a in Theorem 2 is verified. Whatever the initial value of Wolbachia-infected mosquito population xt and uninfected mosquito population yt are, yt will extinct and xt will persist, this means that the result of the competition is that Wolbachia successfully spreads throughout the mosquito population. In addition, it can be shown from the conditions that increasing the birth rate b1 and reducing the competition coefficient α and the death rate d1 will be beneficial to this result.Figure 2
The red solid and the black dotted curves represent the numbers of infected and uninfected mosquitoes inWolbachia, respectively. The parameter values are α=β=1,b1=3,b2=2,d1=0.4,andd2=0.3. The initial conditions are given by x0,y0=1,2 in (a), x0,y0=3,1 in (b), and x0,y0=10,10 in (c), and all of them show E1 is globally asymptotically stable in Case A.
(a)
(b)
(c)From Figure3, we see that in Case B, namely, the fitness cost case, both E1 and E2 are locally asymptotically stable, which is consistent with the result of b in Theorem 2. Each equilibrium has its own attractive regions and the initial values of Wolbachia-infected mosquitoes xt and uninfected mosquitoes yt will affect the trend of the solution; if the solution approaches E1, then yt will become extinct and xt will persist, that is to say the result of the competition is that Wolbachia successfully spreads throughout the mosquito population. On the contrary, if the solution approaches E2, then xt will go to extinct and yt will persist, namely, Wolbachia diffusion fails.Figure 3
The red solid and the black dotted curves represent the numbers of infected and uninfected mosquitoes inWolbachia, respectively. The parameter values are α=β=1,b1=2,b2=3,d1=0.4,andd2=0.3. It satisfies the condition of Case B, both E1 and E2 are locally asymptotically stable, and they have their own attractive regions. The solutions with different initial values will approach to different equilibria. The initial value x0,y0=2,1 and x0,y0=10,5 will lead to the solution approach E1 as t⟶∞, shown in (a) and (b), while the initial value x0,y0=1,3 and x0,y0=10,10 will lead to the solution approach E2 as t⟶∞, shown in (c) and (d).
(a)
(b)
(c)
(d)From above, we find the factors that affect the competition results are the competition coefficient within (or between) species, the birth rate, the death rate, and the initial value of the population. Reducing the competition coefficient and the death rate of the population and increasing the birth rate and the initial value of the population will be beneficial to the lasting survival of the population. It also makes sense in biological terms.Finally, we give an explanation on the stable manifold ofE3. Theorem 2 shows that E3 is a saddle whose stable manifold separates the attractive regions of E1 and E2. However, the specific position of the stable manifold cannot be determined. Here, we use numerical simulation to roughly estimate the position of the stable curve by using different initial values, as shown in Figure 4, and observing which equilibrium that the solution will tend to. Through figure analysis, we claim that the stable curve can be expressed as a smooth function that starts at E0, strictly increasing in the positive direction, and its image biased towards the X-axis. The stable manifold is composed of two separation curves with E3 as the dividing point. Therefore, the attraction region of E2 is larger than that of E1.Figure 4
The parameters and some conclusions are given in Figure3. E3 is the saddle, and its stable manifolds separate the attractive regions of E1 and E2. The solutions with different initial values will approach to different equilibria. The initial value x0,y0=1,1 and x0,y0=5,4 will lead to the solution approaches E1 as t⟶∞, shown in (a) and (c), while the initial value x0,y0=1,2 and x0,y0=6,5 will lead to the solution approaches E2 as t⟶∞, shown in (b) and (d).
(a)
(b)
(c)
(d)


## 4. Conclusion and Discussion

Dengue fever is a mosquito-borne disease, which has great harm to human society. Due to the lack of vaccines and efficient clinical cures, we need to control mosquito population to block the spread of the disease. An innovative and effective method to control mosquitoes is to employWolbachia, which has led to a growing number of researchers building models to study the dynamics of Wolbachia transmission. Considering that the collection data of mosquitoes in the wild are discrete, we established a discrete competition model to study the conditions for Wolbachia to successful spread in mosquitoes.Model (3) is a discrete competition model of overlapping generations of Wolbachia-infected mosquitoes population xt and noninfected mosquitoes population yt. We showed the global asymptotic properties of the four equilibria of the model through elaborate analysis. First of all, the trivial equilibrium E0 is globally asymptotically stable when b1≤d1 and b2≤d2, mosquito populations will go to extinct, otherwise E0 is unstable. Then, we obtain a complete result of E1 and E2. If b1>d1,b2≤d2, then E1 is globally asymptotically stable, which means that Wolbachia-infected mosquitoes population xt will persist and noninfected mosquitoes yt will become extinct, namely, Wolbachia successfully spreads. On the contrary, if b1≤d1,b2>d2, then E2 is globally asymptotically stable, namely, Wolbachia fails to spread. Furthermore, we study the case when b1>d1,b2>d2. In the fitness benefit case, namely, 1/βb2/d2−1<1/αb1/d1−1, and then E2 is a saddle while E1 is globally asymptotically stable and E3 is not exist; that is, Wolbachia successfully spreads. In the fitness cost case, namely, 1/βb2/d2−1<1/αb1/d1−1, both E1 and E2 are locally asymptotically stable while E3 is a saddle, E1 and E2 have their own attractive regions, and the stable manifold of E3 is separated from their attracting domain. The size of the initial population value will affect its own persistence.Numerical simulations are also provided to demonstrate these theoretical results. We mainly showed the simulation under conditionb1>d1,b2>d2, and we found that the simulation results are consistent with the theoretical results. In particular, since we cannot determine the exact position of the stable manifold of E3, we showed its approximate position through simulations.There are some limitations of the model presented in this paper. When we built the model, we only took into account the competitive factors affecting the birth rate of the population, not the death rate. In fact, both intraspecific and interspecific competition can affect mortality in populations. Also, we assumec11=c12=α and c21=c22=β to obtain a simpler model. It is obvious that if c11≠c12,c21≠c22, the model will be more reasonable and can better describe the actual situation of competition between two mosquito populations. However, it presents great challenges in equilibrium analysis.



---

*Source: 8987490-2020-02-27.xml*