Modeling Wolbachia Diffusion in Mosquito Populations by Discrete Competition Model

Dengue fever is caused by dengue virus and transmitted by Aedes mosquitoes. A promising avenue to control this disease is to infect the wild Aedes population with the bacteriumWolbachia 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.


Introduction
Dengue fever is a mosquito-borne infectious disease caused by dengue virus, mainly transmitted by Aedes 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 bacterium Wolbachia 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].
is 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. e strategies of using Wolbachia for mosquito vector control mainly include population suppression and population replacement. e 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 of Wolbachia 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].
ey 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][15][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 of Wolbachia in mosquito populations from the perspective of competition, which is rarely mentioned in the existing literature. e most frequently used equation for studying the discrete model is the Beverton-Holt equation: where b > 0 denotes the birth rate and c 11 > 0 denotes the intraspecific competition coefficient. If b > 1, then all solutions of the equation with x 0 > 0 tend monotonically to the equilibrium x � (b − 1)/c 11 . If x t 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: where b 1 , b 2 > 0 denote the birth rate of x t , y t , respectively. c 11 , c 22 > 0 denote the intraspecific competition coefficient and c 12 , c 21 > 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 x t represents the number of mosquito population infected with Wolbachia at t and y t represents the number of uninfected mosquito population at t, and their birth rates are b 1 and b 2 , respectively. Since the mosquitoes are overlapping generation populations, we introduce d 1 , d 2 as the mortality rates of x t and y t , respectively. en, their survival rates are, respectively, 1 − d 1 and 1 − d 2 . Because Wolbachia often induces fitness cost [10], and it is reasonable to assume that 0 < d 2 < d 1 < 1. Furthermore, considering that both x t and y t 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 c 11 � c 12 � α and c 21 � c 22 � β 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 y t is replaced by b 2 y t /(x t + y t ); hence, the discrete model of competition between two species is as follows: where all the coefficients are positive, and we mainly discuss this model. e structure of this paper is arranged as follows. In Section 2, 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 E 0 , boundary equilibrium E 1 , E 2 , and positive equilibrium E 3 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.

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 Let σ � 1 − d 1 , then we have x t+1 − σx t ≤ (b 1 /α). Furthermore, we can obtain a series of inequalities Add them up, and we get As σ ∈ (0, 1), it is easy to obtain the following inequality for any given initial value x 0 : lim sup Similarly, for y t , we have lim sup 2 Discrete Dynamics in Nature and Society erefore, is 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: us, (3) can be modified as us, E 0 (0, 0) can be regarded as an equilibrium of (3). e other equilibria of (3) satisfy the following algebraic equations: ere are two nonnegative boundary equilibria of model (12), Besides, a positive equilibrium (3), we have the following conclusions.

Theorem 1. For model
Proof. e Jacobian of model (3) is e Jacobians evaluated at E 1 and E 2 read as respectively. eir 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 ⟶ ∞.
the system is monotonically decreasing, and if Obviously, x t+1 � x t , y t+1 � y t if and only if x t � y t � 0 for all t ≥ 0. Otherwise, we claim that 0 ≤ x t+1 < x t and 0 ≤ y t+1 < y t for all t ≥ 0. It follows that lim t⟶+∞ x t and lim t⟶+∞ y t exist and must be zero. erefore, E 0 is globally asymptotically stable. For the remaining cases, i.e., b 1 When y t � 0, from the first equation of model (3), we have . It follows that the solution emanating from the X-axis with 0 < x 0 < (1/α) ((b 1 /d 1 ) − 1) will be away from E 0 along the X-axis, so E 0 is unstable. e eigenvalues of J 1 belong to the interval (0, 1), and hence E 1 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 E 1 . erefore, E 1 is globally asymptotically stable.
In this case, we can prove the results in a similar way to case (b), so we omit it.
□ From now on, we assume that b 1 > d 1 , b 2 > d 2 . By a similar argument, we know that E 0 is unstable. erefore, in the sequel, we will mainly study the boundary equilibria E 1 , E 2 and positive equilibrium E 3 .
We claim that E 3 exists if and only if We analyze it graphically, from the first equation of (12), we can obtain the perpendicular isoclines denoted by L 1 , which is a segment of lines with slope − 1 and intercept (1/α)((b 1 /d 1 ) − 1). From the second equation of (12), we can obtain the horizontal isoclines denoted by L 2 , which is a parabola, and it goes through three fixed points (0, 0), (− (1/β), 0), and (0, (1/β)((b 2 /d 2 ) − 1)). We distinguish two cases depending on the position of L 1 and L 2 as shown in Figure 1. e inequalities satisfied by the coefficients that correspond to these cases can be listed as follows: In this case, there exist three equilibria E 0 , E 1 , and E 2 , but E 3 does not exist.
In this case, there exist four equilibria E 0 , E 1 , E 2 , and E 3 , the positive equilibrium E 3 is the intersection of L 1 and L 2 .
For the above two cases, a feasible biological interpretation of the conditions is that if (1/β)((b 2 /d 2 ) − 1) < (1/α)((b 1 /d 1 ) − 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/β)((b 2 /d 2 ) − 1) > (1/α)((b 1 /d 1 ) − 1), then the environment is less favorable for infected mosquitoes, and we claim that the infection has a fitness cost [10]. e stability analysis of the equilibria in both cases will be shown in the following theorem.
en, we have the following conclusions.
(a) In the fitness benefit case , both E 1 and E 2 are locally asymptotically stable while E 3 is a saddle.

Proof
(a) Case A: from J 2 (see (16)), we know that E 2 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 E 2 . By Liu and Elaydi's theorem, all solutions must approach E 1 . Besides, the eigenvalues of J 1 (see (15)) are positive and less than 1, so E 1 is locally asymptotically stable. erefore, E 1 is globally asymptotically stable.
(b) Case B: all of the eigenvalues of J 1 (see (15)) and J 2 (see (16)) are positive and less than 1, so both E 1 and E 2 are locally asymptotically stable. e stability of E 3 is discussed below. e Jacobian at E 3 takes the form According to the Jury criteria [27], the equilibrium E 3 is locally asymptotically stable if the Jury condition holds, where tr J 3 is the trace and det J 3 is the determinant of J 3 . If at least one of these inequalities is reversed, then the equilibrium is unstable. e inequalities are equivalent to the following three inequalities.
Discrete Dynamics in Nature and Society e first part of eorem 2 predicts that if the environment favors Wolbachia-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 E 3 , whose stable curves separate the attraction regions of E 1 and E 2 . We illustrate our results by numerical simulations in the next section.

Numerical Results
In order to show the dynamic stability of the model (3), we use MATLAB technical computing software for numerical simulations. e simulations for (a) and (b) in eorem 2 is mainly conducted by using different parameters. Some results are shown as follows.
From Figure 2 we see that in Case A, namely, the fitness benefit case, E 1 is globally asymptotically stable. e result of (a) in eorem 2 is verified. Whatever the initial value of Wolbachia-infected mosquito population x t and uninfected mosquito population y t are, y t will extinct and x t 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 b 1 and reducing the competition coefficient α and the death rate d 1 will be beneficial to this result.
From Figure 3, we see that in Case B, namely, the fitness cost case, both E 1 and E 2 are locally asymptotically stable, which is consistent with the result of (b) in eorem 2. Each equilibrium has its own attractive regions and the initial values of Wolbachia-infected mosquitoes x t and uninfected mosquitoes y t will affect the trend of the solution; if the solution approaches E 1 , then y t will become extinct and x t 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 E 2 , then x t will go to extinct and y t will persist, namely, Wolbachia diffusion fails.
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 of E 3 . eorem 2 shows that E 3 is a saddle whose stable manifold separates the attractive regions of E 1 and E 2 . 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. rough figure analysis, we claim that the stable curve can be expressed as a smooth function that starts at E 0 , strictly increasing in the positive direction, and its image biased towards the X-axis. e stable manifold is composed of two separation curves with E 3 as the dividing point. erefore, the attraction region of E 2 is larger than that of E 1 .

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 employ Wolbachia, 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 x t and noninfected mosquitoes population y t . We showed the global asymptotic properties of the four equilibria of the model through elaborate analysis. First of all, the trivial equilibrium E 0 is globally asymptotically stable when b 1 ≤ d 1 and b 2 ≤ d 2 , mosquito populations will go to extinct, otherwise E 0 is unstable. en, we obtain a complete result of E 1 and E 2 . If b 1 > d 1 , b 2 ≤ d 2 , then E 1 is globally asymptotically stable, which means that Wolbachia-infected mosquitoes population x t will persist and noninfected mosquitoes y t will become extinct, namely, Wolbachia successfully spreads. On the contrary, if b 1 ≤ d 1 , b 2 > d 2 , then E 2 is globally asymptotically stable, namely, Wolbachia fails to spread. Furthermore, we study the case when b 1 > d 1 , b 2 > d 2 . In the fitness benefit case, namely, (1/β)((b 2 /d 2 ) − 1) < (1/α) ((b 1 /d 1 ) − 1), and then E 2 is a saddle while E 1 is globally asymptotically stable and E 3 is not exist; that is, Wolbachia successfully spreads. In the fitness cost case, namely, (1/β)((b 2 /d 2 ) − 1) < (1/α)((b 1 /d 1 ) − 1), both E 1 and E 2 are locally asymptotically stable while E 3 is a saddle, E 1 and E 2  Figure 4: e parameters and some conclusions are given in Figure 3. E 3 is the saddle, and its stable manifolds separate the attractive regions of E 1 and E 2 . e solutions with different initial values will approach to different equilibria. e initial value (x 0 , y 0 ) � (1, 1) and (x 0 , y 0 ) � (5, 4) will lead to the solution approaches E 1 as t ⟶ ∞, shown in (a) and (c), while the initial value (x 0 , y 0 ) � (1, 2) and (x 0 , y 0 ) � (6, 5) will lead to the solution approaches E 2 as t ⟶ ∞, shown in (b) and (d).
have their own attractive regions, and the stable manifold of E 3 is separated from their attracting domain. e 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 condition b 1 > d 1 , b 2 > d 2 , 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 E 3 , we showed its approximate position through simulations.
ere 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 assume c 11 � c 12 � α and c 21 � c 22 � β to obtain a simpler model. It is obvious that if c 11 ≠ c 12 , c 21 ≠ c 22 , 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.

Data Availability
No data were used to support the study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.