New revised simple models for interactive wild and sterile mosquito populations and their dynamics

ABSTRACT

Based on previous research, we formulate revised, new, simple models for interactive wild and sterile mosquitoes which are better approximations to real biological situations but mathematically more tractable. We give basic investigations of the dynamical features of these simple models such as the existence of equilibria and their stability. Numerical examples to demonstrate our findings and brief discussions are also provided.

KEYWORDS:

• Allee effects
• density dependence
• stability
• sterile mosquito releases
• vector-borne diseases

AMS SUBJECT CLASSIFICATION:

• 34D20
• 92D25
• 92D30
• 92D40

1. Introduction

Mosquito-borne diseases, such as malaria and dengue fever, are considerable public health concerns worldwide. These diseases are transmitted between humans by blood-feeding mosquitoes. No vaccines are available and an effective way to prevent these mosquito-borne diseases is to control mosquitoes. Among the mosquitoes control measures, the sterile insect technique (SIT) 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,21].

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 incorporating different strategies in releasing sterile mosquitoes have been formulated and studied in [9,16]. The wild and sterile mosquitoes are assumed homogeneous in these models which makes the model reasonably simple and the analysis tractable. Nevertheless, since mosquitoes undergo complete metamorphosis going through four distinct stages of development during a lifetime: egg, pupa, larva, and adult [8], simple stage-structured models are formulated in [17] where the wild mosquitoes are divided into two groups one of which is composed of the three aquatic stages. 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,18]. Density dependence then is included in both the homogeneous and stage-structured models in [9,16,17] to make all model systems nonlinear.

Comparing the homogeneous and stage-structured models, the model formulations in [17] are more reasonable and realistic. However, since our goal is, after having a fundamental understanding of the dynamics for the interactive mosquitoes, to have the mosquito models incorporated into disease transmission models for the mosquito-borne diseases [14,15], the three dimensional models for the mosquitoes make the analysis mathematically more difficult.

We also notice 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 density dependence for the sterile mosquitoes, which are basically adult mosquitoes, can be ignored, and the density dependence for the wild mosquitoes need to be focused on the newborns survivals.

Based on the better understanding we obtained from the model analysis in [17], we revise the models studied in [9] and formulate new simple models in this paper. To make the models as simple as possible, we still assume homogeneous populations such that no stage distinction is concerned. We focus the density dependence on the larvae development process and assume the death rates for wild and sterile mosquitoes are density independent [19,20,22]. We first give general modelling descriptions in Section 2. We then consider three strategies of releases, similar to those in [9,17] to formulate different models, and give complete mathematical analyses for the model dynamics in Sections 3, 4, and 5, respectively. Numerical examples are provided as well to confirm our analytic results, and we finally give brief discussions on our findings in Section 6.

2. Model formulation

Based on the works in [2–6], homogeneous models for the interactive wild and sterile mosquitoes were formulated in [9]. Let $w\left(t\right)$ be the number of wild mosquitoes and $g\left(t\right)$ the number of sterile mosquitoes at time t. Density dependence is assumed for the death rates of both wild and sterile mosquitoes. The model, in a more general setting, has the form (1) $\begin{array}{rl}\frac{\mathrm{d}w}{\mathrm{d}t}& =\left(C\left(N\right)\frac{aw}{w+g}-({\mu }_1+{\xi }_1(w+g\left)\right)\right)w,\\ \frac{\mathrm{d}g}{\mathrm{d}t}& =B(\cdot )-({\mu }_2+{\xi }_2(w+g\left)\right)g,\end{array}$(1) where $C\left(N\right)$ is the number of matings per individual, per unit of time, with N=w + g; a is the number of wild offspring produced per mate; ${\mu }_i$ and ${\xi }_i$, i=1,2, are the density independent and dependent death rates of the wild and sterile mosquitoes, respectively, and $B(\cdot )$ is the release rate of the sterile mosquitoes. With different mating rates, possible Allee effects, and three strategies of releases, fundamental analysis was carried out and results were obtained. System (1(1) $\begin{array}{rl}\frac{\mathrm{d}w}{\mathrm{d}t}& =\left(C\left(N\right)\frac{aw}{w+g}-({\mu }_1+{\xi }_1(w+g\left)\right)\right)w,\\ \frac{\mathrm{d}g}{\mathrm{d}t}& =B(\cdot )-({\mu }_2+{\xi }_2(w+g\left)\right)g,\end{array}$(1) ) is based on the classical Lotka–Volterra model. The populations are homogeneous and density dependence is assumed in both the wild and sterile mosquitoes death rates.

Considering the distinct metamorphic stages of mosquitoes, stage-structured models were formulated and studied in [17]. The three aquatic stages of wild mosquitoes were grouped into one class, called the larvae class and denoted by J, and adult stage was kept as a separate class and denoted by A in the models. Introducing the number of sterile mosquitoes at time t, denoted by $g\left(t\right)$, into the two classes of wild mosquitoes, an interactive model, in more general setting, was described by the following system (2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{A+g}A-\frac{\alpha J}{1+J}-(d_0+d_{1}J)J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\frac{\alpha J}{1+J}-{\mu }_{1}A,\\ \frac{\mathrm{d}g}{\mathrm{d}t}& =B(\cdot )-{\mu }_{2}g,\end{array}$(2) where σ is the number of offspring produced per wild mosquito through all matings, per unit of time, $d_0$ is the density independent death rate and $d_1$ is the density dependent death factor for larvae, ${\mu }_i$, i=1,2, are the adults death rates which are assumed density independent. Based on three strategies of releases, fundamental investigations were performed. Stage structure is included in the model. Note that the sterile mosquitoes released are adult mosquitoes and the competition between adult mosquitoes are relatively weak. Moreover, while the intraspecific competition, or the effect of crowding, represents a major density dependent source in the population dynamics, the crowding basically takes place in water. The death rates for adult wild mosquitoes and sterile mosquitoes are then assumed density independent.

Each of the two models has its strong as well as weak points. In the mean time, from the relatively full investigations, it was shown that the dynamic features of the two models (1(1) $\begin{array}{rl}\frac{\mathrm{d}w}{\mathrm{d}t}& =\left(C\left(N\right)\frac{aw}{w+g}-({\mu }_1+{\xi }_1(w+g\left)\right)\right)w,\\ \frac{\mathrm{d}g}{\mathrm{d}t}& =B(\cdot )-({\mu }_2+{\xi }_2(w+g\left)\right)g,\end{array}$(1) ) and (2(2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{A+g}A-\frac{\alpha J}{1+J}-(d_0+d_{1}J)J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\frac{\alpha J}{1+J}-{\mu }_{1}A,\\ \frac{\mathrm{d}g}{\mathrm{d}t}& =B(\cdot )-{\mu }_{2}g,\end{array}$(2) ) are similar. As mentioned in Section 1, our goal is, after gaining a better understanding of the interactive dynamics of the wild and sterile mosquitoes, to have the mosquitoes models incorporated into the transmission models of diseases. Thus, it is beneficial and even necessary in practice to keep the mosquitoes models as simple as possible, but more accurate and realistic. To this end, learning from the strong points of the two models in systems (1(1) $\begin{array}{rl}\frac{\mathrm{d}w}{\mathrm{d}t}& =\left(C\left(N\right)\frac{aw}{w+g}-({\mu }_1+{\xi }_1(w+g\left)\right)\right)w,\\ \frac{\mathrm{d}g}{\mathrm{d}t}& =B(\cdot )-({\mu }_2+{\xi }_2(w+g\left)\right)g,\end{array}$(1) ) and (2(2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{A+g}A-\frac{\alpha J}{1+J}-(d_0+d_{1}J)J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\frac{\alpha J}{1+J}-{\mu }_{1}A,\\ \frac{\mathrm{d}g}{\mathrm{d}t}& =B(\cdot )-{\mu }_{2}g,\end{array}$(2) ), we propose in this paper new revised models for the interactive wild and sterile mosquitoes. We keep using the two dimensional homogeneous models based on system (1(1) $\begin{array}{rl}\frac{\mathrm{d}w}{\mathrm{d}t}& =\left(C\left(N\right)\frac{aw}{w+g}-({\mu }_1+{\xi }_1(w+g\left)\right)\right)w,\\ \frac{\mathrm{d}g}{\mathrm{d}t}& =B(\cdot )-({\mu }_2+{\xi }_2(w+g\left)\right)g,\end{array}$(1) ) but assume that the density dependence only affects the survivals of the larvae and hence appears only in the birth term of the wild mosquitoes.

Let $S_v$ be the number of wild mosquitoes at time t. In the absence of sterile mosquitoes, the dynamics of $S_v$ are governed by (3) $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=P\left(S_v\right)-{\mu }_v{S}_v,$(3) where $P\left(S_v\right)$ is the birth function, and constant ${\mu }_v$ is the death rate of the wild mosquitoes.

We assume that the birth function has the form $P\left(S_v\right)=C\left(S_v\right)a{S}_v(1-{\xi }_v{S}_v),$ where $C\left(S_v\right)$ is the mating rate, a is the maximum number of survived offspring reproduced per individual, and ${\xi }_v$ is the carrying capacity parameter such that $1-{\xi }_v{S}_V$ describes the density dependent survival probability.

If the mating rate is constant, written as $C\left(S_v\right):=C$, Equation (3(3) $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=P\left(S_v\right)-{\mu }_v{S}_v,$(3) ) becomes the classical logistic equation $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=Ca{S}_v(1-{\xi }_v{S}_v)-{\mu }_v{S}_v=\left(\right(Ca-{\mu }_v)-Ca{\xi }_v{S}_v)S_v.$

Considering the possible difficult for mosquitoes to find their mates when the population size of mosquitoes is small, we assume an Allee effect such that the mating rate takes a form of Holling-II type with $C\left(S_v\right)=c{S}_v/(1+S_v)$. Then Equation (3(3) $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=P\left(S_v\right)-{\mu }_v{S}_v,$(3) ) becomes (4) $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=\frac{a{S}_v}{1+S_v}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,$(4) where we merge c and a but still denote it as a.

Clearly, the trivial equilibrium $S_v=0$ is a locally asymptotically stable equilibrium for Equation (4(4) $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=\frac{a{S}_v}{1+S_v}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,$(4) ). Set the right-hand side of Equation (4(4) $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=\frac{a{S}_v}{1+S_v}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,$(4) ) equal to zero. A positive equilibrium of system (4(4) $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=\frac{a{S}_v}{1+S_v}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,$(4) ) satisfies $-\frac{S_v{Q}_w\left(S_v\right)}{1+S_v}=0,$ where (5) $Q_w\left(S_v\right):=a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v,$(5) with $r_0:=a-{\mu }_v$, the intrinsic growth rate of wild mosquitoes in the absence of sterile mosquitoes.

Function $Q_w\left(S_v\right)$ has real roots (6) $S_{vw}^{\pm }=\frac{r_0\pm \sqrt{r_0^2-4a{\xi }_v{\mu }_v}}{2a{\xi }_v},$(6) provided $r_0^2>4a{\xi }_v{\mu }_v$. Thus, there exists no, one, or two positive equilibria of Equation (4(4) $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=\frac{a{S}_v}{1+S_v}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,$(4) ) if $r_0<2\sqrt{a{\xi }_v{\mu }_v}$, $r_0=2\sqrt{a{\xi }_v{\mu }_v}$, or $r_0>2\sqrt{a{\xi }_v{\mu }_v}$, respectively.

Suppose $r_0>2\sqrt{a{\xi }_v{\mu }_v}$ such that there exist two positive equilibria. It follows from Equations (4(4) $\frac{\mathrm{d}{S}_v}{\mathrm{d}t}=\frac{a{S}_v}{1+S_v}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,$(4) ) – (6(6) $S_{vw}^{\pm }=\frac{r_0\pm \sqrt{r_0^2-4a{\xi }_v{\mu }_v}}{2a{\xi }_v},$(6) ) that, at a positive equilibrium, $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}{S}_v}\left(\frac{a{S}_v}{1+S_v}(1-{\xi }_v{S}_v)-{\mu }_v\right)& =\frac{a(1-{\xi }_v{S}_v)}{(1+S_v{)}^2}-\frac{a{\xi }_v{S}_v}{1+S_v}=\frac{a(1-{\xi }_v{S}_v)}{{\mu }_v(1+S_v{)}^2}\left({\mu }_v-a{\xi }_v{S}_v^2\right)\\ & =\frac{a(1-{\xi }_v{S}_v)}{{\mu }_v(1+S_v{)}^2}\left({\mu }_v-(r_0{S}_v-{\mu }_v)\right)\\ & =\frac{a(1-{\xi }_v{S}_v)}{{\mu }_v(1+S_v{)}^2}\left({\mu }_v-\frac{r_0^2}{2a{\xi }_v}\mp \frac{r_0\sqrt{r_0^2-4a{\xi }_v{\mu }_v}}{2a{\xi }_v}\right)\\ & =\frac{a(1-{\xi }_v{S}_v)}{{\mu }_v(1+S_v{)}^2}\left(4a{\xi }_v{\mu }_v-r_0^2\mp r_0\sqrt{r_0^2-4a{\xi }_v{\mu }_v}\right)\\ & =\frac{a(1-{\xi }_v{S}_v)\sqrt{r_0^2-4a{\xi }_v{\mu }_v}}{{\mu }_v(1+S_v{)}^2}\left(-\sqrt{r_0^2-4a{\xi }_v{\mu }_v}\mp r_0\right).\end{array}$ Thus $\begin{array}{rl}& \frac{\mathrm{d}}{\mathrm{d}{S}_v}\left(\frac{a{S}_v}{1+S_v}(1-{\xi }_v{S}_v)-{\mu }_v\right)\left(S_{vw}^{-}\right)>0,\\ & \frac{\mathrm{d}}{\mathrm{d}{S}_v}\left(\frac{a{S}_v}{1+S_v}(1-{\xi }_v{S}_v)-{\mu }_v\right)\left(S_{vw}^{+}\right)<0,\end{array}$ and then positive equilibrium $S_{vw}^{-}$ is unstable, and $S_{vw}^{+}$ is asymptotically stable.

After sterile mosquitoes are released into the wild mosquito population, we let $S_g$ be the number of sterile mosquitoes at time t, and assume that the interactive birth rate of wild mosquitoes follows the harmonic mean. Then, if there is no Allee effect included, the birth function has the form $a{S}_v/(S_v+S_g)$ and the interactive dynamics are governed by the following system (7) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =B\left(S_v\right)-{\mu }_v{S}_g.\end{array}$(7) If an Allee effect is considered, we then have the following system (8) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =B\left(S_v\right)-{\mu }_v{S}_g.\end{array}$(8)

3. Constant releases

In this section, we suppose that the sterile mosquitoes are constantly released such that $B(\cdot ):=b$, a constant. Then no mating difficulty will be concerned and based on system (7(7) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =B\left(S_v\right)-{\mu }_v{S}_g.\end{array}$(7) ), we have (9) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b-{\mu }_v{S}_g.\end{array}$(9)

Define set ${\mathrm{\Omega }}_1$ as ${\mathrm{\Omega }}_1:=\left\{(S_v,S_g):0\le S_v\le \frac{r_0}{a{\xi }_v},0\le S_g\le \frac{b}{{\mu }_v}\right\},$ where we assume $r_0>0$. Then ${\mathrm{\Omega }}_1$ is a positive invariant set for system (9(9) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b-{\mu }_v{S}_g.\end{array}$(9) ). We assume $(S_v,S_g)\in {\mathrm{\Omega }}_1$ in this section.

System (9(9) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b-{\mu }_v{S}_g.\end{array}$(9) ) has a boundary equilibrium $E_0^c:=(0,g^0)$, where $g^0=b/{\mu }_v$.

A positive equilibrium of system (9(9) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b-{\mu }_v{S}_g.\end{array}$(9) ) satisfies $\frac{a{\mu }_v{S}_v}{{\mu }_v{S}_v+b}(1-{\xi }_v{S}_v)-{\mu }_v:=-\frac{{\mu }_v{Q}_1\left(S_v\right)}{{\mu }_v{S}_v+b}=0,$ where (10) $Q_1\left(S_v\right):=a{\xi }_v{S}_v^2-r_0{S}_v+b,$(10) which has two solutions (11) $S_{vb}^{\pm }=\frac{r_0\pm \sqrt{r_0^2-4a{\xi }_{v}b}}{2a{\xi }_v}.$(11) Define a threshold value of releases for system (9(9) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b-{\mu }_v{S}_g.\end{array}$(9) ) by (12) $b_c^c:=\frac{r_0^2}{4a{\xi }_v}.$(12) Then there exist no, one, or two positive solutions given in (11(11) $S_{vb}^{\pm }=\frac{r_0\pm \sqrt{r_0^2-4a{\xi }_{v}b}}{2a{\xi }_v}.$(11) ) to $Q_1=0$ in Equation (10(10) $Q_1\left(S_v\right):=a{\xi }_v{S}_v^2-r_0{S}_v+b,$(10) ), that is, no, one, or two positive equilibria of system (9(9) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b-{\mu }_v{S}_g.\end{array}$(9) ), if $b>b_c^c$, $b=b_c^c$, or $b<b_c^c$ , respectively.

We next investigate the stability of the equilibria. The Jacobian matrix at an equilibrium has the form $J^{\left(1\right)}:=\left(\begin{array}{cc}{J}_{11}^{\left(1\right)}& -\frac{a{S}_v^2}{(S_v+g^0{)}^2}\\ 0& -{\mu }_v\end{array}\right),$ where $J_{11}^{\left(1\right)}=-{\mu }_v$ at the boundary equilibrium. Thus the boundary equilibrium $E_0^c$ is always asymptotically stable. Furthermore, if there exists no positive interior equilibrium, since set ${\mathrm{\Omega }}_1$ is positively invariant, according to the Poincaré-Bendixson Theorem, $E_0^c$ is globally asymptotically stable.

At a positive equilibrium, $S_v+S_g=\frac{a{S}_v(1-{\xi }_v{S}_v)}{{\mu }_v},$ and thus, in addition to following from Equation (10(10) $Q_1\left(S_v\right):=a{\xi }_v{S}_v^2-r_0{S}_v+b,$(10) ), $\begin{array}{rl}{J}_{11}^{\left(1\right)}& =S_v{\left.\left(\frac{a{S}_v(1-{\xi }_v{S}_v)}{S_v+S_g}-{\mu }_v\right)\right|}_{S_v}=a{S}_v\left(\frac{S_g(1-{\xi }_v{S}_v)}{(S_v+S_g{)}^2}-\frac{{\xi }_v{S}_v}{S_v+S_g}\right)\\ & =\frac{a{S}_v}{(S_v+S_g{)}^2}\left(S_g\right(1-{\xi }_v{S}_v)-{\xi }_v{S}_v(S_v+S_g\left)\right)=\frac{a{S}_v(1-{\xi }_v{S}_v)}{{\mu }_v(S_v+S_g{)}^2}(b-a{\xi }_v{S}_v^2)\\ & =\frac{a{S}_v(1-{\xi }_v{S}_v)}{{\mu }_v(S_v+S_g{)}^2}(2b-r_0{S}_v)=\frac{a{S}_v(1-{\xi }_v{S}_v)}{{\mu }_v(S_v+S_g{)}^2}\left(2b-\frac{r_0^2}{2a{\xi }_v}\mp \frac{r_0\sqrt{r_0^2-4a{\xi }_{v}b}}{2a{\xi }_v}\right)\\ & =\frac{a{S}_v(1-{\xi }_v{S}_v)}{{\mu }_v(S_v+S_g{)}^{2}2a{\xi }_v}(4a{\xi }_{v}b-r_0^2\mp r_0\sqrt{r_0^2-4a{\xi }_{v}b})\\ & =\frac{a{S}_v(1-{\xi }_v{S}_v)\sqrt{r_0^2-4a{\xi }_{v}b}}{{\mu }_v(S_v+S_g{)}^{2}2a{\xi }_v}(-\sqrt{r_0^2-4a{\xi }_{v}b}\mp r_0).\end{array}$

Then, at positive equilibria $E_c^{-}:=(S_{vb}^{-},g^0)$ and $E_c^{+}:=(S_{vb}^{+},g^0)$, $J_{11}^{\left(1\right)}\left(S_{vb}^{-}\right)=\frac{a{S}_v(1-{\xi }_v{S}_v)\sqrt{r_0^2-4a{\xi }_{v}b}}{{\mu }_v(S_v+S_g{)}^{2}2a{\xi }_v}(-\sqrt{r_0^2-4a{\xi }_{v}b}+r_0)>0,$ and $J_{11}^{\left(1\right)}\left(S_{vb}^{+}\right)=\frac{a{S}_v(1-{\xi }_v{S}_v)\sqrt{r_0^2-4a{\xi }_{v}b}}{{\mu }_v(S_v+S_g{)}^{2}2a{\xi }_v}(-\sqrt{r_0^2-4a{\xi }_{v}b}-r_0)<0,$ which implies that $E_c^{-}$ is a saddle point and $E_c^{+}$ is a locally asymptotically stable node. In summary, we have

Theorem 3.1:

Assume $r_0=a-{\mu }_v>0$. System (9(9) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b-{\mu }_v{S}_g.\end{array}$(9) ) has a boundary equilibrium $E_c^0=(0,g^0)$ with $g^0=b/{\mu }_v,$ which is globally asymptotically stable if there exists no positive equilibrium, and locally asymptotically stable if a positive equilibrium exists. System (9(9) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b-{\mu }_v{S}_g.\end{array}$(9) ) has no, one, or two positive equilibria if $b>b_c^c,$ $b=b_c^c,$ or $b<b_c^c,$ respectively, where the release threshold $b_c^c$ is given in (12(12) $b_c^c:=\frac{r_0^2}{4a{\xi }_v}.$(12) ). If there exist two positive equilibria, $E_c^{-}=(S_{vb}^{-},g^0)$ and $E_c^{+}=(S_{vb}^{+},g^0),$ where $S_{vb}^{\pm }$ are given in (11(11) $S_{vb}^{\pm }=\frac{r_0\pm \sqrt{r_0^2-4a{\xi }_{v}b}}{2a{\xi }_v}.$(11) ), $E_c^{-}$ is an unstable saddle point and $E_c^{+}$ is a locally asymptotically stable node.

We give an example to illustrate the results in Theorem 3.1 as follows.

Example 1:

Parameters are given as (13) $a=10,{\mu }_v=0.5,{\xi }_v=0.3,$(13) and thus the release threshold is $b_c^c=7.5208$.

Boundary equilibrium $E_c^0=(0,12)$ is asymptotically stable. With $b=6<b_c^c$, there exist two positive equilibria $E_c^{-}=(0.8713,12)$ and $E_c^{+}=(2.2953,12)$. Equilibrium $E_c^{-}$ is a saddle point and $E_c^{+}$ is a locally asymptotically stable node. Solutions approach either $E_c^0$ or $E_c^{+}$ depending on their initial values, as shown in the left side of Figure 1. For $b=9>b_c^c$, there exists no positive equilibrium and $E_c^0$ is globally asymptotically stable. All solutions approach $E_c^0$ as shown in the right side of Figure 1.

Figure 1. Parameters are given in (13(13) $a=10,{\mu }_v=0.5,{\xi }_v=0.3,$(13) ) and the release threshold is $b_c^c=7.5208$. For $b=6<b_c^c$, there are two positive equilibria $E_c^{-}=(0.8713,12)$ and $E_c^{+}=(2.2953,12)$. Equilibrium $E_1^c$ is a saddle point and $E_c^{+}$ is a locally asymptotically stable node, solutions approach $E_c^{+}$ or the locally asymptotically stable boundary equilibrium $E_c^0$ depending on their initial values, as shown in the left-side figure. For $b=9>b_c^c$, no positive equilibrium exists and $E_c^0$ becomes globally asymptotically stable. All solutions approach $E_c^0$ as shown in the right-side figure.

4. Releases proportional to the wild mosquitoes

As in [9,17], to have a more optimal and economically effective strategy of releasing sterile mosquitoes in an area where the population size of wild mosquitoes is relatively small, we may consider to keep closely sampling or surveillance of the wild mosquito population and let the rate of releases be proportional to the population size of the wild mosquitoes such that $B\left(S_v\right)=b{S}_v$ where b is a constant. In the meantime, in a situation where a small mosquito population size may cause possible mating difficulty, we assume an Allee effect and thus the interactive system for the wild and sterile mosquitoes, based on system (8(8) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =B\left(S_v\right)-{\mu }_v{S}_g.\end{array}$(8) ), has the form (14) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b{S}_v-{\mu }_v{S}_g.\end{array}$(14)

Define set ${\mathrm{\Omega }}_2$ as ${\mathrm{\Omega }}_2:=\left\{(S_v,S_g):0\le S_v\le S_{vw}^{+},0\le S_g\le \frac{b}{{\mu }_v}{S}_{vw}^{+}\right\},$ where $S_{vw}^{+}$ is given in (6(6) $S_{vw}^{\pm }=\frac{r_0\pm \sqrt{r_0^2-4a{\xi }_v{\mu }_v}}{2a{\xi }_v},$(6) ). Then ${\mathrm{\Omega }}_2$ is a positive invariant set for system (14(14) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b{S}_v-{\mu }_v{S}_g.\end{array}$(14) ). We assume $(S_v,S_g)\in {\mathrm{\Omega }}_2$ in this section.

The origin $(0,0)$ is a trivial equilibrium and is always locally asymptotical stable. A positive equilibrium of system (14(14) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b{S}_v-{\mu }_v{S}_g.\end{array}$(14) ) satisfies the following equation $\frac{a{S}_v(1-{\xi }_v{S}_v)}{1+S_v+S_g}-{\mu }_v=-\frac{{\mu }_v{Q}_2\left(S_v\right)}{{\mu }_v+({\mu }_v+b)S_v}=0,$ where (15) $Q_2\left(S_v\right):=a{\xi }_v{S}_v^2-(r_0-b)S_v+{\mu }_v,$(15) which has two roots (16) $S_{vp}^{\pm }=\frac{r_0-b\pm \sqrt{(r_0-b{)}^2-4a{\xi }_v{\mu }_v}}{2a{\xi }_v}.$(16) If $b\ge r_0$, then there are no positive roots to $Q_2\left(S_v\right)=0$, that is, no positive equilibria to system (14(14) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b{S}_v-{\mu }_v{S}_g.\end{array}$(14) ).

Assume $b<r_0$. Then $(r_0-b{)}^2-4a{\xi }_v{\mu }_v\ge 0$, if and only if $0<b\le r_0-2\sqrt{a{\xi }_v{\mu }_v}$. Define the threshold value of releases as (17) $b_p^c:=r_0-2\sqrt{a{\xi }_v{\mu }_v}.$(17) Then system (4) has no, one, or two positive equilibria of system (9(9) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b-{\mu }_v{S}_g.\end{array}$(9) ) if $b>b_p^c$, $b=b_p^c$, or $b<b_p^c$, respectively.

Assume $b<b_p^c$ so that there exist two positive equilibria, denoted by $E_p^{-}:=(S_{vp}^{-},b{S}_{vp}^{-}/{\mu }_v)$ and $E_p^{+}:=(S_{vp}^{+},b{S}_{vp}^{+}/{\mu }_v)$. We investigate the stability of the positive equilibria as follows.

The Jacobian matrix at an equilibrium has the form $J^{\left(2\right)}:=\left(\begin{array}{cc}{J}_{11}^{\left(2\right)}& -\frac{a{\mu }_v^2{S}_v^2(1-{\xi }_v{S}_v)}{({\mu }_v+({\mu }_v+b)S_v{)}^2}\\ b& -{\mu }_v\end{array}\right),$ where $\begin{array}{rl}{J}_{11}^{\left(2\right)}& =S_v{\left.\left(\frac{a{S}_v(1-{\xi }_v{S}_v)}{1+S_v+S_g}-{\mu }_v\right)\right|}_{S_v}=\frac{a{S}_v}{(1+S_v+S_g{)}^2}\left(-{\xi }_v{S}_v^2+(1-2{\xi }_v{S}_v)(1+S_g)\right)\\ & =-\frac{a{\mu }_v{S}_v}{({\mu }_v+({\mu }_v+b)S_v{)}^2}\left(\right({\mu }_v+2b){\xi }_v{S}_v^2+(2{\mu }_v{\xi }_v-b)S_v-{\mu }_v).\end{array}$ Thus, writing ${\mathrm{\Delta }}_1:=\frac{({\mu }_v+({\mu }_v+b)S_v{)}^2}{a{\mu }_v^2{S}_v}det{J}^{\left(2\right)},$ we have (18) $\begin{array}{rl}{\mathrm{\Delta }}_1& =({\mu }_v+2b){\xi }_v{S}_v^2+(2{\mu }_v{\xi }_v-b)S_v-{\mu }_v+b{S}_v(1-{\xi }_v{S}_v)\\ & =({\mu }_v+b){\xi }_v{S}_v^2+2{\mu }_v{\xi }_v{S}_v-{\mu }_v.\end{array}$(18) It follows from Equation (15(15) $Q_2\left(S_v\right):=a{\xi }_v{S}_v^2-(r_0-b)S_v+{\mu }_v,$(15) ) that, at an equilibrium, ${\xi }_v{S}_v=\frac{1}{a}\left(\right(r_0-b)S_v-{\mu }_v).$ Substituting it into Equation (18(18) $\begin{array}{rl}{\mathrm{\Delta }}_1& =({\mu }_v+2b){\xi }_v{S}_v^2+(2{\mu }_v{\xi }_v-b)S_v-{\mu }_v+b{S}_v(1-{\xi }_v{S}_v)\\ & =({\mu }_v+b){\xi }_v{S}_v^2+2{\mu }_v{\xi }_v{S}_v-{\mu }_v.\end{array}$(18) ) yields (19) $\begin{array}{rl}a{\mathrm{\Delta }}_1& =({\mu }_v+b)\left(\right(r_0-b)S_v-{\mu }_v)+2a{\mu }_v{\xi }_v{S}_v-a{\mu }_v\\ & =\left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)S_v-(a+{\mu }_v+b){\mu }_v.\end{array}$(19) Then substituting (16(16) $S_{vp}^{\pm }=\frac{r_0-b\pm \sqrt{(r_0-b{)}^2-4a{\xi }_v{\mu }_v}}{2a{\xi }_v}.$(16) ) into Equation (19(19) $\begin{array}{rl}a{\mathrm{\Delta }}_1& =({\mu }_v+b)\left(\right(r_0-b)S_v-{\mu }_v)+2a{\mu }_v{\xi }_v{S}_v-a{\mu }_v\\ & =\left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)S_v-(a+{\mu }_v+b){\mu }_v.\end{array}$(19) ), we have $\begin{array}{rl}2a^2{\xi }_v{\mathrm{\Delta }}_1& =\left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)(r_0-b)-2a{\xi }_v(a+{\mu }_v+b){\mu }_v\\ & \pm \left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)\sqrt{(r_0-b{)}^2-4a{\xi }_v{\mu }_v}.\end{array}$ Notice that $\begin{array}{rl}& \left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)(r_0-b)-2a{\xi }_v(a+{\mu }_v+b){\mu }_v\\ & =({\mu }_v+b)(r_0-b{)}^2+2a{\mu }_v{\xi }_v(r_0-b-(a+{\mu }_v+b))\\ & =({\mu }_v+b)(r_0-b{)}^2-4a{\mu }_v{\xi }_v({\mu }_v+b)=({\mu }_v+b\left)\right((r_0-b{)}^2-4a{\mu }_v{\xi }_v).\end{array}$ Then $\begin{array}{rl}2a^2{\xi }_v{\mathrm{\Delta }}_1& =({\mu }_v+b)\left(\right(r_0-b{)}^2-4a{\mu }_v{\xi }_v)\\ & \pm \left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)\sqrt{(r_0-b{)}^2-4a{\xi }_v{\mu }_v},\end{array}$ and thus we have the determinant of the Jacobian at $E_p^{+}$, $det{J}^{\left(2\right)}\left(S_{vp}^{+}\right)>0$.

At $E_p^{-}$, writing ${\mathrm{\Delta }}_2:=\frac{2a^2{\xi }_v{\mathrm{\Delta }}_1}{\sqrt{(r_0-b{)}^2-4a{\xi }_v{\mu }_v}},$ we have $\begin{array}{rl}{\mathrm{\Delta }}_2& =({\mu }_v+b)\sqrt{(r_0-b{)}^2-4a{\mu }_v{\xi }_v}-\left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)\\ & =\frac{({\mu }_v+b{)}^2((r_0-b{)}^2-4a{\mu }_v{\xi }_v)-\left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v{)}^2}{({\mu }_v+b)\sqrt{(r_0-b{)}^2-4a{\mu }_v{\xi }_v}+\left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)}\\ & =-\frac{4a{\mu }_v{\xi }_v\left(\right({\mu }_v+b{)}^2+({\mu }_v+b)(r_0-b))+4a^2{\mu }_v^2{\xi }_v^2}{({\mu }_v+b)\sqrt{(r_0-b{)}^2-4a{\mu }_v{\xi }_v}+\left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)}\\ & =-\frac{4a{\mu }_v{\xi }_v\left(\right({\mu }_v+b\left)\right({\mu }_v+r_0\left)\right)+a{\mu }_v{\xi }_v)}{({\mu }_v+b)\sqrt{(r_0-b{)}^2-4a{\mu }_v{\xi }_v}+\left(\right({\mu }_v+b\left)\right(r_0-b)+2a{\mu }_v{\xi }_v)}<0,\end{array}$ which implies that the determinant of $J^{\left(2\right)}$ at $E_p^{-}$ has $det{J}^{\left(2\right)}\left(S_{vp}^{-}\right)<0$. Hence $E_p^{-}$ is an unstable saddle point.

Moreover, at $E_p^{+}$, $\begin{array}{rl}\mathrm{tr}{J}^{\left(2\right)}\left(S_v\right)& =-\frac{a{\mu }_v{S}_v}{({\mu }_v+({\mu }_v+b)S_v{)}^2}\left(\right({\mu }_v+2b){\xi }_v{S}_v^2+(2{\mu }_v{\xi }_v-b)S_v-{\mu }_v)-{\mu }_v\\ & =-\frac{{\mu }_v{\mathrm{\Delta }}_3}{({\mu }_v+({\mu }_v+b)S_v{)}^2},\end{array}$ where (20) $\begin{array}{rl}{\mathrm{\Delta }}_3:=& a{S}_v\left(\right({\mu }_v+2b){\xi }_v{S}_v^2+(2{\mu }_v{\xi }_v-b)S_v-{\mu }_v)+({\mu }_v+({\mu }_v+b)S_v{)}^2\\ & =a{S}_v\left(\right({\mu }_v+2\left){\xi }_v{S}_v^2+2{\mu }_v{\xi }_v{S}_v-{\mu }_v\right)+b{S}_{v}a{\xi }_v{S}_v^2-ab{S}_v^2+\left({\mu }_v+\right({\mu }_v+b)S_v{)}^2\\ & =a{S}_v{\mathrm{\Delta }}_1+b{S}_{v}a{\xi }_v{S}_v^2-ab{S}_v^2+({\mu }_v+({\mu }_v+b)S_v{)}^2.\end{array}$(20) It follows from Equation (15(15) $Q_2\left(S_v\right):=a{\xi }_v{S}_v^2-(r_0-b)S_v+{\mu }_v,$(15) ) that $a{\xi }_v{S}_v^2=(r_0-b)S_v-{\mu }_v.$ Substituting it into Equation (20(20) $\begin{array}{rl}{\mathrm{\Delta }}_3:=& a{S}_v\left(\right({\mu }_v+2b){\xi }_v{S}_v^2+(2{\mu }_v{\xi }_v-b)S_v-{\mu }_v)+({\mu }_v+({\mu }_v+b)S_v{)}^2\\ & =a{S}_v\left(\right({\mu }_v+2\left){\xi }_v{S}_v^2+2{\mu }_v{\xi }_v{S}_v-{\mu }_v\right)+b{S}_{v}a{\xi }_v{S}_v^2-ab{S}_v^2+\left({\mu }_v+\right({\mu }_v+b)S_v{)}^2\\ & =a{S}_v{\mathrm{\Delta }}_1+b{S}_{v}a{\xi }_v{S}_v^2-ab{S}_v^2+({\mu }_v+({\mu }_v+b)S_v{)}^2.\end{array}$(20) ) then yields $\begin{array}{rl}{\mathrm{\Delta }}_3=& a{S}_v{\mathrm{\Delta }}_1+b{S}_v(r_0-b)S_v-b{S}_v{\mu }_v-ab{S}_v^2+{\mu }_v^2+2{\mu }_v({\mu }_v+b)S_v+({\mu }_v+b{)}^2{S}_v^2\\ & =a{S}_v{\mathrm{\Delta }}_1+b{S}_v^2(r_0-b-a)+{\mu }_v^2+{\mu }_v(2{\mu }_v+b)S_v+({\mu }_v+b{)}^2{S}_v^2\\ & =a{S}_v{\mathrm{\Delta }}_1-b{S}_v^2({\mu }_v+b)+{\mu }_v^2+{\mu }_v(2{\mu }_v+b)S_v+b({\mu }_v+b)S_v^2+{\mu }_v({\mu }_v+b)S_v^2\\ & =a{S}_v{\mathrm{\Delta }}_1+{\mu }_v^2+{\mu }_v(2{\mu }_v+b)S_v+{\mu }_v({\mu }_v+b)S_v^2>0.\end{array}$ Thus $\mathrm{tr}{J}^{\left(2\right)}\left(S_{vp}^{+}\right)=-{\mathrm{\Delta }}_3<0$ and hence positive equilibrium $E_p^{+}$ is a locally asymptotically stable node or spiral.

In summary, we have the following results.

Theorem 4.1:

The origin $(0,0)$ is a trivial equilibrium of system (14(14) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b{S}_v-{\mu }_v{S}_g.\end{array}$(14) ), which is always locally asymptotically stable. Assume $r_0>2\sqrt{a{\xi }_v{\mu }_v}$. Then system (14(14) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =b{S}_v-{\mu }_v{S}_g.\end{array}$(14) ) has no, one, or two positive equilibria if $b>b_c^p,$ $b=b_c^p,$ or $b<b_c^p,$ respectively, where the release threshold $b_c^p$ is given in (17(17) $b_p^c:=r_0-2\sqrt{a{\xi }_v{\mu }_v}.$(17) ). If there exist two positive equilibria, $E_p^{-}=(S_{vp}^{-},b{S}_{vp}^{-}/{\mu }_v)$ and $E_p^{+}=(S_{vp}^{+},b{S}_{vp}^{+}/{\mu }_v),$ where $S_{vp}^{\pm }$ are given in (16(16) $S_{vp}^{\pm }=\frac{r_0-b\pm \sqrt{(r_0-b{)}^2-4a{\xi }_v{\mu }_v}}{2a{\xi }_v}.$(16) ), equilibrium $E_p^{-}$ is an unstable saddle point and $E_p^{+}$ is a locally asymptotically stable node or spiral.

A numerical example for the proportional releases of sterile mosquitoes is given below.

Example 2:

With the same parameters as in Example 1, that is, (21) $a=10,{\mu }_v=0.5,{\xi }_v=0.3,$(21) the release threshold is $b_p^c=7.0505$.

When $b=6<b_p^c$, there exist two positive equilibria $E_p^{-}=(0.6667,2)$ and $E_p^{+}=(1,12)$. Equilibrium $E_p^{-}$ is a saddle point and $E_p^{+}$ is a locally asymptotically spiral. Solutions approach either $E_p^{+}$ or the origin depending on their initial values as shown in the left side of Figure 2. For $b=9>b_p^c$, there exists no positive equilibrium and the origin is globally asymptotically stable. All solutions approach the origin as shown in the right side of Figure 2.

Figure 2. Parameters are given in (21(21) $a=10,{\mu }_v=0.5,{\xi }_v=0.3,$(21) ) and the release threshold is $b_p^c=7.0505$. For $b=6<b_p^c$, there are two positive equilibria $E_p^{-}=(0.6667,2)$ and $E_p^{+}=(1,12)$. Equilibrium $E_p^{-}$ is a saddle point and $E_p^{+}$ is a locally asymptotically spiral. Solutions approach either $E_p^{+}$ or the origin depending on their initial values as shown in the left-side figure. For $b=9>b_p^c$, there exists no positive equilibrium. The origin is globally asymptotically stable. All solutions approach the origin as shown in the right-side figure.

5. Proportional releases with saturation

A new strategy to compromise the two strategies studied in Sections 3 and 4 was proposed in [9,17], where the release rate is proportional to $S_v$ as $S_v$ small, but is saturated and approaches a constant as $S_v$ sufficiently large. Let $B\left(S_v\right)=b{S}_v/(1+S_v)$, and the interactive dynamics of the wild and sterile mosquitoes are governed by the following system (22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22)

Define set ${\mathrm{\Omega }}_3$ as ${\mathrm{\Omega }}_3:=\left\{(S_v,S_g):0\le S_v\le S_{vw}^{+},0\le S_g\le \frac{b}{{\mu }_v}\right\},$ where $S_{vw}^{+}$ is given in (6(6) $S_{vw}^{\pm }=\frac{r_0\pm \sqrt{r_0^2-4a{\xi }_v{\mu }_v}}{2a{\xi }_v},$(6) ). Then ${\mathrm{\Omega }}_3$ is a positive invariant set for system (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ). We assume $(S_v,S_g)\in {\mathrm{\Omega }}_3$ in this section.

The origin $(0,0)$ is a trivial equilibrium and is always locally asymptotical stable. A positive equilibrium of system (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ) satisfies the following equation $\frac{a{S}_v(1-{\xi }_v{S}_v)}{1+S_v+\frac{b{S}_v}{{\mu }_v(1+S_v)}}-{\mu }_v=-\frac{{\mu }_v{Q}_3\left(S_v\right)}{{\mu }_v(1+S_v{)}^2+b{S}_v}=0,$ where $Q_3:=a{\xi }_v{S}_v^3+(a{\xi }_v-r_0)S_v^2+(b+{\mu }_v-r_0)S_v+{\mu }_v.$ Clearly, the positive equilibria of (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ) correspond to the positive roots of $Q_3\left(S_v\right)=0$. Since there exist at most two positive roots of $Q_3\left(S_v\right)=0$, there exist at most two positive equilibria of (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ).

Rewrite $Q_3\left(S_v\right)=Q_w\left(S_v\right)(1+S_v)+b{S}_v,$ where $Q_w$ is given in (5(5) $Q_w\left(S_v\right):=a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v,$(5) ). For fixed parameters, if there is no positive root for $Q_w\left(S_v\right)=0$, that is, no intersection between the curve of $Q_w\left(S_v\right)$ and the positive $S_v$-axis, the curve of $Q_3\left(S_v\right)$ has no intersection with the positive $S_v$-axis for $b\ge 0$. Hence there is no positive equilibrium of system (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ) if $r_0\le 2\sqrt{a{\xi }_v{\mu }_v}$.

Suppose $r_0>2\sqrt{a{\xi }_v{\mu }_v}$. The curve of $Q_3\left(S_v\right)$ for b=0, that is the curve of $Q_w\left(S_v\right)$, intersects the positive $S_v$-axis twice. It is clear that the curve of $Q_3\left(S_v\right)$ moves up as b increases and is completely above the positive $S_v$-axis as b sufficiently large. Hence, there exists a threshold value of b, denoted by $b_H^c$, such that the curve is tangent to the positive $S_v$-axis, for all of other parameters fixed. Thus system (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ) has no, one, or two, positive equilibria, if $b>b_H^c$, $b=b_H^c$, or $b<b_H^c$, respectively. Indeed, the threshold value $b_H^c$ can be determined as follows.

At a positive critical point $(S_v^H,b_H^c)$, the following two equations are satisfied simultaneously: (23) $\begin{array}{rl}{Q}_3\left(S_v\right)& =Q_w\left(S_v\right)(S_v+1)+b{S}_v=(a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v)(S_v+1)+b{S}_v=0,\\ Q_3^{\prime }\left(S_v\right)& =Q_w^{\prime }\left(S_v\right)(S_v+1)+Q_w\left(S_v\right)+b\\ & =(2a{\xi }_v{S}_v-r_0)(S_v+1)+a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v+b=0,\end{array}$(23) which results in a unique $S_v^H$ satisfying (24) $2a{\xi }_v(S_v^H{)}^3+a{\xi }_v(S_v^H{)}^2-{\mu }_v=0,$(24) and then (25) $b_H^c:=r_0-{\mu }_v-3a{\xi }_v(S_v^H{)}^2+2(r_0-a{\xi }_v)S_v^H.$(25)

We next work on the stability of the positive equilibria.

The Jacobian matrix at an equilibrium has the form $J^{\left(3\right)}:=\left(\begin{array}{cc}{J}_{11}^{\left(3\right)}& J_{12}^{\left(3\right)}\\ \frac{b}{(1+S_v{)}^2}& -{\mu }_v\end{array}\right),$ where $\begin{array}{rl}{J}_{11}^{\left(3\right)}& =S_v{\left.\left(\frac{a{S}_v(1-{\xi }_v{S}_v)}{1+S_v+S_g}-{\mu }_v\right)\right|}_{S_v}=S_v\left(\frac{a(1-{\xi }_v{S}_v)}{1+S_v+S_g}-\frac{a{S}_v(1+{\xi }_v(1+S_g\left)\right)}{(1+S_v+S_g{)}^2}\right)\\ & ={\mu }_v-\frac{a{S}_v^2(1+{\xi }_v(1+S_g\left)\right)}{(1+S_v+S_g{)}^2},\end{array}$ and $J_{12}^{\left(3\right)}=S_v{\left.\left(\frac{a{S}_v(1-{\xi }_v{S}_v)}{1+S_v+S_g}-{\mu }_v\right)\right|}_{S_g}=-\frac{a{S}_v^2(1-{\xi }_v{S}_v)}{(1+S_v+S_g{)}^2}=-\frac{{\mu }_v{S}_v}{1+S_v+S_g}.$

It follows from Equations (23(23) $\begin{array}{rl}{Q}_3\left(S_v\right)& =Q_w\left(S_v\right)(S_v+1)+b{S}_v=(a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v)(S_v+1)+b{S}_v=0,\\ Q_3^{\prime }\left(S_v\right)& =Q_w^{\prime }\left(S_v\right)(S_v+1)+Q_w\left(S_v\right)+b\\ & =(2a{\xi }_v{S}_v-r_0)(S_v+1)+a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v+b=0,\end{array}$(23) ) that, at a positive equilibrium, (26) $a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v=a{S}_v({\xi }_v{S}_v-1)+{\mu }_v(1+S_v)=-\frac{b{S}_v}{1+S_v}.$(26) Since $S_g=\frac{b{S}_v}{{\mu }_v(1+S_v)},$ then Equation (26(26) $a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v=a{S}_v({\xi }_v{S}_v-1)+{\mu }_v(1+S_v)=-\frac{b{S}_v}{1+S_v}.$(26) ) leads to $S_g=\frac{a{S}_v(1-{\xi }_v{S}_v)}{{\mu }_v}-(1+S_v).$ Thus (27) $1+S_v+S_g=\frac{a{S}_v(1-{\xi }_v{S}_v)}{{\mu }_v},$(27) and $1+{\xi }_v(1+S_g)=1+{\xi }_v\left(\frac{a{S}_v(1-{\xi }_v{S}_v)}{{\mu }_v}-S_v\right)=\frac{1-{\xi }_v{S}_v}{{\mu }_v}(a{\xi }_v{S}_v+{\mu }_v).$ We then have (28) $\begin{array}{rl}{J}_{11}^{\left(3\right)}& ={\mu }_v-\frac{a{S}_v^2(1+{\xi }_v(1+S_g\left)\right)}{(1+S_v+S_g{)}^2}={\mu }_v-\frac{a{S}_v^2(1-{\xi }_v{S}_v)(a{\xi }_v{S}_v+{\mu }_v){\mu }_v^2}{a^2{S}_v^2(1-{\xi }_v{S}_v{)}^2{\mu }_v}\\ & ={\mu }_v-\frac{(a{\xi }_v{S}_v+{\mu }_v){\mu }_v}{a(1-{\xi }_v{S}_v)}=\frac{{\mu }_v}{a(1-{\xi }_v{S}_v)}\left(a\right(1-{\xi }_v{S}_v)-a{\xi }_v{S}_v-{\mu }_v)\\ & =\frac{{\mu }_v}{a(1-{\xi }_v{S}_v)}\left(a\right(1-2{\xi }_v{S}_v)-{\mu }_v).\end{array}$(28)

On the other hand, it follows from Equations (23(23) $\begin{array}{rl}{Q}_3\left(S_v\right)& =Q_w\left(S_v\right)(S_v+1)+b{S}_v=(a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v)(S_v+1)+b{S}_v=0,\\ Q_3^{\prime }\left(S_v\right)& =Q_w^{\prime }\left(S_v\right)(S_v+1)+Q_w\left(S_v\right)+b\\ & =(2a{\xi }_v{S}_v-r_0)(S_v+1)+a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v+b=0,\end{array}$(23) ) and (26(26) $a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v=a{S}_v({\xi }_v{S}_v-1)+{\mu }_v(1+S_v)=-\frac{b{S}_v}{1+S_v}.$(26) ) that $\begin{array}{rl}{Q}_3^{\prime }& =(2a{\xi }_v{S}_v-r_0)(S_v+1)-\frac{b{S}_v}{1+S_v}+b=\left(a\right(2{\xi }_v{S}_v-1)+{\mu }_v)(S_v+1)+\frac{b}{1+S_v}\\ & =-\left(a\right(1-2{\xi }_v{S}_v)-{\mu }_v)(S_v+1)+\frac{b}{1+S_v}.\end{array}$ Then $a(1-2{\xi }_v{S}_v)-{\mu }_v=\frac{1}{1+S_v}\left(\frac{b}{1+S_v}-Q_3^{\prime }\right),$ and thus $J_{11}^{\left(3\right)}=\frac{{\mu }_v}{a(1-{\xi }_v{S}_v)(1+S_v)}\left(\frac{b}{1+S_v}-Q_3^{\prime }\right).$ Hence, we have $\begin{array}{rl}det{J}^{\left(3\right)}& =-{\mu }_v{J}_{11}^{\left(3\right)}-\frac{b}{(1+S_v{)}^2}{J}_{12}^{\left(3\right)}\\ & =\frac{{\mu }_v^2{Q}_3^{\prime }}{a(1-{\xi }_v{S}_v)(1+S_v)}-\frac{b{\mu }_v^2}{a(1-{\xi }_v{S}_v)(1+S_v{)}^2}+\frac{ba{S}_v^2(1-\xi S_v)}{(1+S_v{)}^2(1+S_v+S_g{)}^2},\end{array}$ and it follows from Equation (27(27) $1+S_v+S_g=\frac{a{S}_v(1-{\xi }_v{S}_v)}{{\mu }_v},$(27) ) that (29) $\begin{array}{rl}det{J}^{\left(3\right)}& =\frac{{\mu }_v^2{Q}_3^{\prime }}{a(1-{\xi }_v{S}_v)(1+S_v)}-\frac{b{\mu }_v^2}{a(1-{\xi }_v{S}_v)(1+S_v{)}^2}+\frac{b{\mu }_v^2}{(1+S_v{)}^{2}a(1-{\xi }_v{S}_v)}\\ & =\frac{{\mu }_v^2{Q}_3^{\prime }}{a(1-{\xi }_v{S}_v)(1+S_v)}.\end{array}$(29)

Suppose there exist two positive equilibria $E_1^H:=(S_v^{\left\{1\right\}},S_g^{\left\{1\right\}})$ and $E_2^H:=(S_v^{\left\{2\right\}},S_g^{\left\{2\right\}})$ with $S_v^{\left\{1\right\}}<S_v^c<S_v^{\left\{2\right\}}$. Then $Q_3^{\prime }\left(S_v^{\left\{1\right\}}\right)<0,Q_3^{\prime }\left(S_v^{\left\{2\right\}}\right)>0,$ which implies, from Equation (29(29) $\begin{array}{rl}det{J}^{\left(3\right)}& =\frac{{\mu }_v^2{Q}_3^{\prime }}{a(1-{\xi }_v{S}_v)(1+S_v)}-\frac{b{\mu }_v^2}{a(1-{\xi }_v{S}_v)(1+S_v{)}^2}+\frac{b{\mu }_v^2}{(1+S_v{)}^{2}a(1-{\xi }_v{S}_v)}\\ & =\frac{{\mu }_v^2{Q}_3^{\prime }}{a(1-{\xi }_v{S}_v)(1+S_v)}.\end{array}$(29) ) that $det{J}^{\left(3\right)}\left(S_v^{\left\{1\right\}}\right)<0,det{J}^{\left(3\right)}\left(S_v^{\left\{2\right\}}\right)>0.$ Thus, $E_1^H$ is an unstable saddle.

Moreover, we have, from Equation (28(28) $\begin{array}{rl}{J}_{11}^{\left(3\right)}& ={\mu }_v-\frac{a{S}_v^2(1+{\xi }_v(1+S_g\left)\right)}{(1+S_v+S_g{)}^2}={\mu }_v-\frac{a{S}_v^2(1-{\xi }_v{S}_v)(a{\xi }_v{S}_v+{\mu }_v){\mu }_v^2}{a^2{S}_v^2(1-{\xi }_v{S}_v{)}^2{\mu }_v}\\ & ={\mu }_v-\frac{(a{\xi }_v{S}_v+{\mu }_v){\mu }_v}{a(1-{\xi }_v{S}_v)}=\frac{{\mu }_v}{a(1-{\xi }_v{S}_v)}\left(a\right(1-{\xi }_v{S}_v)-a{\xi }_v{S}_v-{\mu }_v)\\ & =\frac{{\mu }_v}{a(1-{\xi }_v{S}_v)}\left(a\right(1-2{\xi }_v{S}_v)-{\mu }_v).\end{array}$(28) ), $\begin{array}{rl}\mathrm{tr}{J}^{\left(3\right)}& =J_{11}^{\left(3\right)}-{\mu }_v=\frac{{\mu }_v}{a(1-{\xi }_v{S}_v)}\left(a\right(1-2{\xi }_v{S}_v)-{\mu }_v)-{\mu }_v\\ & =\frac{{\mu }_v}{a(1-{\xi }_v{S}_v)}\left(a\right(1-2{\xi }_v{S}_v)-{\mu }_v-a(1-{\xi }_v{S}_v\left)\right)\\ & =\frac{{\mu }_v}{a(1-{\xi }_v{S}_v)}(-a{\xi }_v{S}_v)-{\mu }_v)<0.\end{array}$ Therefore, positive equilibrium $E_2^H$ is a locally asymptotically stable node or saddle.

In summary, we have

Theorem 5.1:

The origin $(0,0)$ is a trivial equilibrium of system (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ), which is always locally asymptotically stable. If $r_0\le 2\sqrt{a{\xi }_v{\mu }_v},$ system (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ) has no positive equilibrium. If $r_0>2\sqrt{a{\xi }_v{\mu }_v},$ there exists a threshold value $b_H^c=b_H^c\left(S_v^H\right)$ given in (25(25) $b_H^c:=r_0-{\mu }_v-3a{\xi }_v(S_v^H{)}^2+2(r_0-a{\xi }_v)S_v^H.$(25) ) with $S_v^H$ being the unique positive solution to Equation (24). System (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ) has no, one, or two positive equilibria if $b>b_H^c,$ $b=b_H^c,$ or $b<b_H^c,$ respectively. If there exist two positive equilibria $E_1^H=(S_v^{\left\{1\right\}},S_g^{\left\{1\right\}})$ and $E_2^H=(S_v^{\left\{2\right\}},S_g^{\left\{2\right\}})$ with $S_v^{\left\{1\right\}}<S_v^H<S_v^{\left\{2\right\}},$ equilibrium $E_1^H$ is an unstable saddle and $E_2^H$ is a locally asymptotically stable node or spiral.

We give a numerical example for system (22(22) $\begin{array}{rl}\frac{\mathrm{d}{S}_v}{\mathrm{d}t}& =\frac{a{S}_v}{1+S_v+S_g}(1-{\xi }_v{S}_v)S_v-{\mu }_v{S}_v,\\ \frac{\mathrm{d}{S}_g}{\mathrm{d}t}& =\frac{b{S}_v}{1+S_v}-{\mu }_v{S}_g.\end{array}$(22) ) as follows.

Example 3:

Given the parameters (30) $a=10,{\mu }_v=0.5,{\xi }_v=0.3,$(30) the solution $S_v^H$ to system (23(23) $\begin{array}{rl}{Q}_3\left(S_v\right)& =Q_w\left(S_v\right)(S_v+1)+b{S}_v=(a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v)(S_v+1)+b{S}_v=0,\\ Q_3^{\prime }\left(S_v\right)& =Q_w^{\prime }\left(S_v\right)(S_v+1)+Q_w\left(S_v\right)+b\\ & =(2a{\xi }_v{S}_v-r_0)(S_v+1)+a{\xi }_v{S}_v^2-r_0{S}_v+{\mu }_v+b=0,\end{array}$(23) ) is $S_v^H=0.31899$ and then $b_H^c=12.231$ in (25(25) $b_H^c:=r_0-{\mu }_v-3a{\xi }_v(S_v^H{)}^2+2(r_0-a{\xi }_v)S_v^H.$(25) ).

With $b=11<b_H^c$, there exist two positive equilibria $E_1^H=(0.60312,8.2768)$ and $E_2^H=(1.7239,13.9232)$, in addition to the trivial equilibrium $(0,0)$ which is a locally asymptotically stable node. Equilibrium $E_1^H$ is a saddle point while $E_2^H$ is a locally asymptotically stable spiral, as shown in the upper left side of Figure 3. As b increases but is less than $b_H^c$, the two positive equilibria move closer. For b=12.07, $E_1^H$ is still a saddle point, $E_2^H$ becomes a locally asymptotically stable node as shown in the upper right side of Figure 3. There is an unstable saddle-node point as $b=b_c^H$, and as $b=13>b_c^H$, there exists no positive equilibrium and all solutions approach the origin as shown in the lower part of Figure 3.

Figure 3. Parameters are given in (30(30) $a=10,{\mu }_v=0.5,{\xi }_v=0.3,$(30) ) and the release threshold is $b_H^c=12.231$. For $b=11<b_H^c$, there are two positive equilibria, one of which is a saddle point and one of which is a locally asymptotically stable spiral, as shown in the upper left-side figure. For b=12.07, the two positive equilibria move closer, where one is still saddle point and the other becomes a locally asymptotically stable node as shown in the upper right-side figure. For $b=13>b_H^c$, there exists no positive equilibrium and all solutions approach the origin as shown in the lower figure.

6. Concluding remarks

Based on the models in [9,17], we revised and formulated new simple models for the interactive wild and sterile mosquitoes in this paper. The main goal is to make the models more appropriately describing the biological situations, and in the meantime mathematically more tractable.

While mosquitoes undergo complete metamorphosis going through distinct stages of development during a lifetime, and hence stage-structured models are needed, the investigations on the stage-structured models in [17] and the simplified homogeneous models without stage structure in [9] showed that they exhibited very similar dynamical features. Since our fundamental goal of these studies is to eventually have the interactive mosquitoes models incorporated into transmission models of mosquito-borne diseases in higher dimensional spaces, a reduction of even one dimension in the model formulation could benefit our further research. Hence, it seems reasonable to keep using two-dimensional homogeneous models for the interaction of wild and sterile mosquitoes but with modifications.

We note that the released sterile mosquitoes are all adults. The competition for resources between the adult mosquitoes are relatively weak and as a result, it seems reasonable to assume density independent death rates for both wild and sterile mosquitoes. In the meantime, the intraspecific competition could represent a major density dependent source for the larvae development, and in particular, the crowding basically takes place in water. Thus, we need to have the density dependence more focused on the birth terms, more specifically, on the survivals of the newborns. Comparing to the models in [9], we thus assume that the birth rate has the form of $a{S}_v(1-{\xi }_v{S}_v)$ with a the maximum reproduction rate. With such revisions, we believe that the models are more biologically appropriate and better fitting our research design.

We fully studied the dynamical features for the models formulated in this paper. We determined the existence of all possible equilibria and completely investigated their stability for the three models with different strategies of releases in Sections 3,4, and 5. We established release threshold values of sterile mosquitoes for each model. As the release rate of sterile mosquitoes exceeds the threshold value, there exists no positive equilibrium, and the boundary equilibrium for the constant releases and the trivial equilibrium for the other two strategies are globally asymptotically stable, which means all wild mosquitoes go extinct eventually. On the other hand, it the release rate is less than the threshold value, all of the three models have two positive interior equilibria one of which is always a saddle point, and one of which is locally asymptotically stable. The stable equilibrium for the constant releases is a node, but it can be either a node or a spiral for the other strategies of releases. This is different from the model dynamics in [9,17].

The model dynamics are relatively simple compared to the models in [9,17], but again since the focus of modelling mosquitoes interactions is on our fundamental goal of having the mosquito models incorporated into disease transmission models, this is exactly what we need. Further studies on the impact of releases of sterile mosquitoes on disease transmissions, based on the interactive mosquito models in this paper, will benefit from such revisions, and will appear in the near future. We would also like to point out that the new models in this paper can be applied to other interacting species as well.

Acknowledgments

The author thanks the handling editor and two anonymous reviewers for their valuable comments and suggestions.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was supported partially by U.S. National Science Foundation grant [DMS-1118150].