Stage-structured wild and sterile mosquito population models and their dynamics

ABSTRACT To study the impact of the sterile insect technique and effects of the mosquitoes' metamorphic stage structure on the transmission dynamics of mosquito-borne diseases, we formulate stage-structured continuous-time mathematical models, based on systems of differential equations, for the interactive dynamics of the wild and sterile mosquitoes. We incorporate different strategies for the releases of sterile mosquitoes in the models and investigate the model dynamics, including the existence of positive equilibria and their stability. Numerical examples are provided to demonstrate the dynamical features of the models.


Introduction
The sterile insect technique (SIT), as one of the mosquitoes control measures, has been applied to reducing or eradicating wild mosquitoes. The SIT is a method of biological control in which the natural reproductive process of mosquitoes is disrupted. Utilizing radical or other chemical or physical methods, male mosquitoes are genetically modified to be sterile which are incapable of producing offspring despite being sexually active. These sterile mosquitoes are then released into the environment to mate with wild mosquitoes that are present in the environment. A wild female that mates with a sterile male will either not reproduce, or produce eggs but the eggs will not hatch. Repeated releases of sterile mosquitoes or releasing a significantly large number of sterile mosquitoes may eventually wipe out or control a wild mosquito population [1,7,23].
Mathematical models have been formulated in the literature to study the interactive dynamics and control of the wild and sterile mosquito populations [2-6, 11, 12]. In particular, models incorporate different strategies in releasing sterile mosquitoes have been formulated and studied in [9,18]. However, the mosquito populations have been assumed homogeneous without distinguishing the metamorphic stages of mosquitoes.
Mosquitoes undergo complete metamorphosis going through four distinct stages of development during a lifetime: egg, pupa, larva, and adult [8]. While interspecific competition and predation may cause larval mortality, intraspecific competition could represent a major density dependent source for the population dynamics, and hence the effect of crowding could be an important factor in the population dynamics of mosquitoes [10,13,20]. Note that the crowding basically takes place in water where the egg, pupa, and larva stages in a mosquito's life cycle present. That is, the mosquito population density dependence mostly exists in the first three aquatic stages. Thus, the assumption of homogenous populations apparently is inappropriate and to have a better understanding of the impact of the releases of sterile mosquitoes, the metamorphic stage structure needs to be included in the model formulations [15-17, 19, 21, 22, 24].
We formulate stage-structured mosquito population models with sterile mosquitoes in this paper. To keep our mathematical models as simple as possible, we follow a line similar to the stage-structured models for mosquitoes in [15][16][17] and group the three aquatic stages of mosquitoes into one class, called larvae, and divide the mosquito population into only two classes, the larvae and adults. We ignore the interspecific competition and predation and assume that the density dependence on the deaths only appears in the larvae. Since the sterile mosquitoes released are adults, the death rate for the sterile mosquitoes is also assumed to be density independent. We first give general modelling descriptions in Section 2. We then study the dynamics of the model, similar to that in [9], where the number of releases of sterile mosquitoes is constant in Section 4. Complete mathematical analysis for the model dynamics is given. We next consider the case where the number of the sterile mosquito releases is proportional to the wild mosquito population size in Section 5. Mathematical analysis and numerical examples are provided to demonstrate the complexity of the model dynamics. To provide a different release strategy, we assume, in Section 6, that the releases are of Holling-II type such that the number of the sterile mosquito releases is proportional to the wild mosquito population size when the wild mosquito population size is small but is saturated and approaches a constant as the wild mosquito population size is sufficiently large. The mathematical analysis for the model dynamics is more complex and partial results are obtained. We finally provide brief discussions on our findings in Section 7.

Model formulation
Consider a wild mosquito population without the presence of sterile mosquitoes. For the simplified stage-structured mosquito population, we group the three aquatic stages into the larvae class, denoted by J, and divide the mosquito population into the larvae class and the adults, denoted by A. We assume that the density dependence exists only in the larvae stage.
We let the birth rate, that is, the oviposition rate of adults be β(·); the rate of emergence from larvae to adults be a function of the larvae with the form of α(1 − κ(J)), where α > 0 is the maximum emergence rate, 0 ≤ κ(J) ≤ 1, with κ(0) = 0, κ (J) > 0, and lim J→∞ κ(J) = 1, is the functional response due to the intraspecific competition [16]. We let the death rate of larvae be a linear function, denoted by d 0 + d 1 J, and the death rate of adults be constant, denoted by μ 1 . Then we arrive at, in the absence of sterile mosquitoes, the following system of equations: We further assume a functional response for κ(J), as in [16], in the form System (1) then becomes We next consider two different cases where there is either or no difficulty for adult mosquitoes to find their mates such that Allee effects are either included or not in the population dynamics.

Wild mosquito population without Allee effects
Suppose mosquito adults have no difficulty to find their mates such that no Allee effects are concerned, and hence the adult birth is constant, simply denoted as β(·) := β. The interactive dynamics for the wild mosquitoes are governed by the following system: It is easy to see that system (2) has no closed orbits by the Bendixson-Dulac theorem. Define the intrinsic growth rate for the mosquito population in system (2) as The dynamics of system (2) can be summarized as follows.
which is a globally asymptotically stable node.

Wild mosquito population with Allee effects
In the case where adult mosquitoes have difficulty in finding their mates, Allee effects are included and the adult birth rate is given by where γ is the Allee effect constant. The stage-structured wild mosquito population model then has the following form: System (4) has a trivial equilibrium (0, 0) which is always asymptotically stable. A positive equilibrium of Equation (4) satisfies which leads to with Here δ 0 is defined as The J component in a positive equilibrium of Equation (4) is then equivalently a positive root of 0 (J) = 0. Notice that K 1 > 0 and K 2 > 0. Thus if δ 0 ≤ 1, there exists no positive equilibrium of system (4).
The Jacobian matrix at a positive equilibrium has the form The trace of 1 is and the determinant of 1 is It follows from Equation (5) that and from Equation (5) again that Substituting Equation (12) into Equation (11) yields It follows from If there exist two positive equilibria of Equation (4), at a positive equilibrium, we have (J 1 ) > 0 and (J 2 ) < 0, and hence Thus E 1 is an unstable saddle point and E 2 is a stable node or spiral. Notice that it follows from the Bendixson-Dulac theorem again that system (4) has no closed orbits. We summarize the results for system (4) as follows. Theorem 2.2: System (4) has a trivial equilibrium (0, 0) which is always locally asymptotically stable. The existence of positive equilibria for system (4) is summarized in the following table.
Here δ 0 is given in Equation (7), and 0 is given in Equation (10). If there exist two positive where J 0 is given in Equation (9), then E 1 is an unstable saddle posit and E 2 is a locally asymptotically node or spiral. System (4) has no closed orbits.

Mosquito population model with sterile mosquitoes
We now consider that sterile mosquitoes are released to a wild mosquito population, and let g(t) be the number of sterile mosquitoes at time t. Since there is no birth for these mosquitoes, the birth rate for sterile mosquitoes is the rate of the releases, denoted by B(·). Because the major density dependence is from the intraspecific competition between larvae, we assume the death rate of sterile mosquitoes is density-independent and is denoted by μ 2 . The dynamics of sterile mosquitoes are then determined by the following simple equation: After the releases of sterile mosquitoes, the interactive mating takes place. By using the harmonic mean, as in [14], in the case where no Allee effects are concerned, the birth function for the wild mosquitoes has the form where σ is the number of offspring produced per wild mosquito through all matings, per unit of time. Based on systems (1), (13), and (14), the dynamics of the interactive wild and sterile mosquitoes are described by the following system: It is easy to see that the first quadrant is positive invariant under the flow of system (15) and the intrinsic growth rate of the wild mosquito population in the absence of sterile mosquitoes is In the case where Allee effects are included, based on systems (4), (13), and (14), the dynamics of the interactive wild and sterile mosquitoes are described by the following system: We consider three different cases with different release strategies as in [9,18] in the following sections.

Constant releases
Assume the sterile mosquitoes are constantly released to the wild mosquito population such that the release rate of sterile mosquitoes is constant, denoted by B(·) := b. Then there is no difficulty for mosquitoes to find mates and based on system (15), the interactive dynamics are governed by the following system: There exists a trivial boundary equilibrium E 0 := (J, A, g) = (0, 0, g 0 ), where g 0 = b/μ 2 , for system (18). Simple linear stability analysis shows that it is locally asymptotically stable.
At a positive equilibrium, we have g = b/μ 2 . Substituting it into the first equation for γ in Equation (5) in Section 2.2, we have the following equation to determine the positive equilibria: and Here we define a threshold value of releases as Clearly, if r 0 ≤ 1, there exists no positive equilibrium for system (18) for any b ≥ 0. If r 0 > 1, there exists no positive equilibrium for system (18) then there exist no, one, or two positive equilibria for system (18) The expression of c is depending on b. We can further determine a threshold value b c s for the existence of positive equilibria as follows. Let function G c (b) be defined as where c is given in Equation (19), and J c satisfies c (J) = 0 and is given as follows: Since c (J c ) = 0 and B i (b) > 0, i = 1,2,3, it follows from (18) has no, one, or two positive equilibria, respectively.
For the asymptotic dynamics of system (18), it is easy to see that the JS-plane is a global attractor for system (18), and the stability analysis for Equation (18) is similar to that for Equation (4). The results are summarized as follows. Theorem 4.1: System (18) has a trivial equilibrium E 0 = (0, 0, g 0 ) which is always locally asymptotically stable. The results for the existence of positive equilibria for system (18) are summarized in the following table. Here r 0 is the intrinsic growth rate of the wild mosquito population in the absence of sterile mosquitoes given in Equation (16), the threshold value of the existence of positive equilibria b c is given in Equation (21), and the threshold value b c s is determined by c (b c s ) = 1 with c given in Equation (22). If there exist two positive equilibria E 1 = (J 1 , A 1 , g 1 ) and where J c is given in Equation (23) and B i , i = 1, 2, 3, are given in Equation (20), equilibrium E 1 is an unstable saddle posit and E 2 is a locally asymptotically node or spiral.
We give a numerical example to illustrate the dynamics of system (18) in Example 1.

Proportional releases
To have a more economically effective strategy for releases of sterile mosquitoes in an area where the population size of wild mosquitoes is relatively small, instead of releasing sterile mosquitoes constantly, we may, by keeping close surveillance of the wild mosquitoes, let the number of releases be proportional to the population size of the wild mosquitos such that B(·) = bA [9,18]. In such a case, there possibly exists difficulty for mosquitoes to find mates when the wild mosquito population size is relatively small. Similarly as in Section 2.2, we then assume Allee effects and that the interactive dynamics between the wild and sterile mosquitoes, based on Equation (17), are governed by the following system: The origin (0, 0, 0) is an equilibrium and is always locally asymptotically stable. We then explore the existence of positive equilibria as follows.
At a positive equilibrium of system (26), we have Then it follows from that the corresponding function for the existence of positive equilibria has the form with C 1 := μ 1 d 1 (αb + μ 2 (μ 1 + α)), d 1 ))), Here the threshold value for the releases, b p , is defined as and System (26) clearly has no positive equilibrium if r p ≤ 1 for any b ≥ 0, or r p > 1 and b ≥ b p . Suppose r p > 1 and b < b p . We define Then system (26) has no, one, or two positive equilibria provided p > 1, p = 1, or p < 1, respectively. Similarly as in Section 4, since p depends on b, we can determine a threshold value of b for the existence of positive equilibria of system (26) as follows.
Let function G p (b) be defined as Then since p (J p ) = 0 and C i (b) > 0, i = 1,2,3, it follows from is an increasing function of b, which implies that p (b) is monotone increasing as b increases. Therefore, there exists a threshold value b p s > 0 such that as b > b p s , b = b p s , or b < b p s , p > 1, p = 1, p < 1, and hence system (26) has no, one, or two positive equilibria, respectively. We next study the stability of the positive equilibria. The Jacobian matrix of system (26) has the form The characteristic equation of the linearized system of (26) at an equilibrium then is given by where P 1 := −tr 2 = μ 1 + μ 2 + a 11 ,  (J 1 , A 1 , g 1 ) and and J p is given in Equation (32), then we have p (J 1 ) < 0 and p (J 2 ) > 0, and hence p (J 1 ) > 0 and p (J 2 ) < 0. By direct calculation, we have (Details are given in Appendix 1). Thus P 3 (E 1 ) < 0 and P 3 (E 2 ) > 0, and since P 3 = − det 2 , we have which implies that E 1 is unstable. Furthermore, at E 2 , direct but tedious calculation shows that if μ 2 ≥ μ 1 , then P 1 P 2 − P 3 > 0. (Details are given in Appendix 2.) Since P 3 (E 2 ) > 0, positive equilibrium E 2 is locally asymptotically stable. In summary, we have the following theorem.
Theorem 5.1: The trivial equilibrium (0, 0, 0) for system (26) is always locally asymptotically stable. The existence results for positive equilibria of system (26) are summarized in the following table.  (J 1 , A 1 , g 1 ) and E 2 = (J 2 , A 2 , g 2 ) with J 1 < J c < J 2 , where J p is given in Equation (23), equilibrium E 1 is an unstable saddle posit and E 2 is a locally asymptotically node or spiral, provided μ 2 ≥ μ 1 . The dynamics of system (26) have been determined in Theorem 5.1 for μ 2 ≥ μ 1 . However, the model dynamics can be complex if μ 2 < μ 1 . While E 1 is always an unstable saddle, E 2 can be locally asymptotically stable or unstable. It can be a node or a spiral. We show the dynamical complexity of system (26) in Example 2.  Figure 2. Solutions initially started near E 2 spiral towards E 2 and solutions initially started away from the closed orbit approach the origin as shown in the right figure in Figure 2.
With the same parameters but b = 12.5724 which is still less than b p s , we again have two positive equilibria E 1 = (0.01074, 0.29752, 19.0844) and E 2 = (0.47174, 8.97486, 575.69148). Equilibrium E 1 is an unstable saddle but E 2 becomes an unstable spiral. There exists a bistable closed orbit as shown in the left figure in Figure 3. Solutions initially started near E 2 spiral towards the closed orbit and solutions initially started away from the closed orbit approach the origin as shown in the right figure in Figure 3. Here only the curves for the wild adults are presented. Notice that the speed of the convergence to E 2 is slow. Here also only the curves for the wild adults are presented.

Proportional releases with saturation
The strategy of proportional releases presented in Section 5, compared to the constant releases, may have an advantage when the size of the wild mosquito population is small since the size of releases is then also small. However, when the wild mosquito population size is large, the release size should presumably also be large, which may exceed our affordability. As in [9,18], we consider a new strategy in which the number of sterile mosquito releases is proportional to the wild mosquito population size when it is small, but is saturated and approaches a constant when the wild mosquito population size increases. Thus, we let the number of the releases be of Holling-II type with the form B(·) = bA/(1 + A). Then the interactive system becomes The origin (0, 0, 0) is a trivial equilibrium and is clearly locally asymptotically stable. At a positive equilibrium, we have where ρ = α/μ 1 . Substituting them into the first equation for positive equilibria in Equation (35), we have and hence the characteristic equation is Similarly as in Section 5, direct calculation yields (Details are given in Appendix 3.) If there exist two positive equilibria, denoted by E 1 = (J 1 , A 1 , g 1 ) and E 2 = (J 2 , A 2 , g 2 ) with A 1 < A 2 , then it follows from H (E 1 ) > 0 and which implies the instability of E 1 . At E 2 , direct but tedious calculation shows that if positive equilibrium E 2 is locally asymptotically stable. The results are summarized as follows.
We note that the stability condition μ 2 ≥ μ 1 for E 2 is only sufficient and very strong. With most sets of parameters, we have unfortunately not been able to find examples where E 2 is unstable or other dynamical features occur.

Concluding remarks
Mosquitoes undergo complete metamorphosis through four distinct stages of development during a lifetime, and the effects of crowding, particularly in the aquatic stages, could be an important factor and thus it is necessary to have the metamorphic stages be included in modelling population dynamics of mosquitoes. Based on the model systems in [9], we introduced the metamorphic stage structure of mosquitoes into dynamical models for interactive wild and sterile mosquitoes to study the impact of the releases of sterile mosquitoes in this paper. We simplified the models by combing the three aquatic metamorphic stages into one group, called larvae as in [15][16][17], and assumed that the density dependence, due to the intraspecific competition, was only on the larvae. We considered three different strategies similarly as in [9,18] for the releases in model systems (18), (26), and (35), respectively. We explored the existence of positive equilibria and studied their stability for all of these model systems. We established threshold release values to determine the existence of positive equilibria for model systems (18) and (26). The threshold value for system (35) was also implicitly obtained. We completely determined the local stability conditions of the equilibria for system (18), and sufficient stability conditions for the positive equilibria for systems (26) and (35) were obtained. The dynamics of system (26) and (35) are complicated. The existence of stable or unstable closed orbits for system (26) is shown in Example 2. While we have obtained sufficient conditions for the stability of the positive equilibria for system (35), the conditions appear too strong and we have been unable to find examples where positive equilibrium E 2 is unstable or where other dynamical phenomena occur when the stability conditions fail.