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.

KEYWORDS:

• Allee effects
• stability
• stage structure
• sterile mosquito releases
• vector-bornediseases

AMS SUBJECT CLASSIFICATIONS:

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

1. 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–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.

2. 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 $\beta (\cdot )$; the rate of emergence from larvae to adults be a function of the larvae with the form of $\alpha (1-\kappa (J\left)\right)$, where $\alpha >0$ is the maximum emergence rate, $0\le \kappa \left(J\right)\le 1$, with $\kappa \left(0\right)=0$, ${\kappa }^{\prime }\left(J\right)>0$, and $\underset{J\to \mathrm{\infty }}{lim}\kappa \left(J\right)=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 ${\mu }_1$. Then we arrive at, in the absence of sterile mosquitoes, the following system of equations: (1) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\beta (\cdot )A-\alpha (1-\kappa (J\left)\right)J-(d_0+d_{1}J)J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\alpha (1-\kappa (J\left)\right)J-{\mu }_{1}A.\end{array}$(1) We further assume a functional response for $\kappa \left(J\right)$, as in [16], in the form $\kappa \left(J\right)=\frac{J}{1+J}.$ System (1(1) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\beta (\cdot )A-\alpha (1-\kappa (J\left)\right)J-(d_0+d_{1}J)J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\alpha (1-\kappa (J\left)\right)J-{\mu }_{1}A.\end{array}$(1) ) then becomes $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\beta (\cdot )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.\end{array}$

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.

2.1. 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 $\beta (\cdot ):=\beta$. The interactive dynamics for the wild mosquitoes are governed by the following system: (2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\beta 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.\end{array}$(2)

It is easy to see that system (2(2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\beta 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.\end{array}$(2) ) has no closed orbits by the Bendixson–Dulac theorem. Define the intrinsic growth rate for the mosquito population in system (2(2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\beta 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.\end{array}$(2) ) as (3) ${\overline{r}}_0:=\frac{\beta \alpha }{(\alpha +d_0){\mu }_1}.$(3) The dynamics of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\beta 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.\end{array}$(2) ) can be summarized as follows.

Theorem 2.1

Theorem 3.1 in [16]

If${\overline{r}}_0\le 1,$ where ${\overline{r}}_0$ is defined in Equation (3(3) ${\overline{r}}_0:=\frac{\beta \alpha }{(\alpha +d_0){\mu }_1}.$(3) ), the trivial equilibrium $(0,0)$ of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\beta 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.\end{array}$(2) ) is a globally asymptotically stable node, and there exists no positive equilibrium. If${\overline{r}}_0>1,$ the trivial equilibrium $(0,0)$ is unstable, and there exists a unique positive equilibrium $J^{\left(0\right)}=\frac{\left(\right(d_0+d_1{)}^2-4d_1(\alpha +d_0)(1-r_0){)}^{1/2}}{2d_1}-\frac{d_0+d_1}{2d_1},A^{\left(0\right)}=\frac{\alpha J^{\left(0\right)}}{{\mu }_1(1+J^{\left(0\right)})},$ which is a globally asymptotically stable node.

2.2. 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 $\beta (\cdot )=\frac{\beta A}{\gamma +A},$ where γ is the Allee effect constant. The stage-structured wild mosquito population model then has the following form: (4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4)

System (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) has a trivial equilibrium $(0,0)$ which is always asymptotically stable. A positive equilibrium of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) satisfies (5) $\begin{array}{rl}\frac{\beta A}{\gamma +A}A& =\frac{\alpha J}{1+J}+(d_0+d_{1}J)J,\\ \frac{\alpha J}{1+J}& ={\mu }_{1}A,\end{array}$(5) which leads to $\frac{\beta \alpha J}{\gamma {\mu }_1+(\gamma {\mu }_1+\alpha )J}\frac{\alpha }{{\mu }_1(1+J)}=\frac{\alpha }{1+J}+d_0+d_{1}J,$ that is ${\mathrm{\Phi }}_0\left(J\right):=\frac{\beta {\alpha }^{2}J}{\gamma {\mu }_1+(\gamma {\mu }_1+\alpha )J}-{\mu }_1(\alpha +(d_0+d_{1}J\left)\right(1+J\left)\right)=\frac{-{\mathrm{\Psi }}_0\left(J\right)}{\gamma {\mu }_1+(\gamma {\mu }_1+\alpha )J}=0,$ where (6) ${\mathrm{\Psi }}_0\left(J\right):=K_1{J}^3+K_2{J}^2+K_{3}J+\gamma {\mu }_1^2(d_0+\alpha ),$(6) with $\begin{array}{rl}{K}_1& :={\mu }_1{d}_1(\gamma {\mu }_1+\alpha ),\\ K_2& :={\mu }_1(\gamma {\mu }_1{d}_1+(\gamma {\mu }_1+\alpha \left)\right(d_0+d_1\left)\right),\\ K_3& :={\mu }_1\left(\gamma {\mu }_1\right(d_0+d_1)+(\gamma {\mu }_1+\alpha \left)\right(d_0+\alpha \left)\right)(1-{\delta }_0).\end{array}$ Here ${\delta }_0$ is defined as (7) ${\delta }_0:=\frac{\beta {\alpha }^2}{{\mu }_1\left(\gamma {\mu }_1\right(d_0+d_1)+(\gamma {\mu }_1+\alpha \left)\right(d_0+\alpha \left)\right)}.$(7)

The J component in a positive equilibrium of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) is then equivalently a positive root of ${\mathrm{\Psi }}_0\left(J\right)=0$. Notice that $K_1>0$ and $K_2>0$. Thus if ${\delta }_0\le 1$, there exists no positive equilibrium of system (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ).

Suppose ${\delta }_0>1$ and thus $K_3<0$. It follows from (8) ${\mathrm{\Psi }}_0^{\prime }\left(J\right)=3K_1{J}^2+2K_{2}J+K_3$(8) that there exists a unique positive root of ${\mathrm{\Psi }}_0^{\prime }\left(J\right)=0$ given by (9) $J_0:=\frac{\sqrt{K_2^2-3K_1{K}_3}-K_2}{3K_1}.$(9) Then $\mathrm{\Psi }\left(J\right)$ (and hence $\mathrm{\Phi }\left(J\right)$) has no, one, or two positive roots; that is, system (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) has no, one, or two positive equilibria if ${\mathrm{\Psi }}_0\left(J_0\right)>0$, ${\mathrm{\Psi }}_0\left(J_0\right)=0$, or ${\mathrm{\Psi }}_0\left(J_0\right)<0$, respectively.

It follows from Equations (6(6) ${\mathrm{\Psi }}_0\left(J\right):=K_1{J}^3+K_2{J}^2+K_{3}J+\gamma {\mu }_1^2(d_0+\alpha ),$(6) ) and (8(8) ${\mathrm{\Psi }}_0^{\prime }\left(J\right)=3K_1{J}^2+2K_{2}J+K_3$(8) ) that $\begin{array}{rl}{\mathrm{\Psi }}_0\left(J_0\right)& =\frac{J_0}{3}(3K_1{J}_0^2+2K_2{J}_0+K_3)+\frac{1}{3}{K}_2{J}_0^2+\frac{2}{3}{K}_3{J}_0+\gamma {\mu }_1^2(d_0+\alpha )\\ & =\frac{1}{3}{K}_2\left(-\frac{2K_2}{3K_1}{J}_0-\frac{K_3}{3K_1}\right)+\frac{2}{3}{K}_3{J}_0+\gamma {\mu }_1^2(d_0+\alpha )\\ & =\left(\frac{2}{3}{K}_3-\frac{2K_2^2}{9K_1}\right)J_0-\frac{K_2{K}_3}{9K_1}+\gamma {\mu }_1^2(d_0+\alpha )\\ & =\frac{2}{9K_1}(3K_1{K}_3-K_2^2)J_0-\frac{K_2{K}_3}{9K_1}+\gamma {\mu }_1^2(d_0+\alpha )\\ & =-\frac{2}{27K_1^2}(K_2^2-3K_1{K}_3)\left(\sqrt{K_2^2-3K_1{K}_3}-K_2\right)-\frac{K_2{K}_3}{9K_1}+\gamma {\mu }_1^2(d_0+\alpha ).\end{array}$ Write (10) ${\mathrm{\Delta }}_0:=\frac{27K_1^2\gamma {\mu }_1^2(d_0+\alpha )}{2(K_2^2-3K_1{K}_3)(\sqrt{K_2^2-3K_1{K}_3}-K_2)+3K_1{K}_2{K}_3}.$(10) Then system (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) has no, one, or two positive equilibria if ${\mathrm{\Delta }}_0>1$, ${\mathrm{\Delta }}_0=1$, or ${\mathrm{\Delta }}_0<1$, respectively.

We next investigate the stability for the positive equilibria of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) as follows.

The Jacobian matrix at a positive equilibrium has the form ${\mathrm{\Lambda }}_1:=\left(\begin{array}{cc}-d_0-2d_{1}J-\frac{\alpha }{(1+J{)}^2}& \frac{\beta A}{\gamma +A}\left(1+\frac{\gamma }{\gamma +A}\right)\\ \frac{\alpha }{(1+J{)}^2}& -{\mu }_1\end{array}\right).$

The trace of ${\mathrm{\Lambda }}_1$ is $\mathrm{tr}{\mathrm{\Lambda }}_1=-d_0-2d_{1}J-\frac{\alpha }{(1+J{)}^2}-{\mu }_1<0,$ and the determinant of ${\mathrm{\Lambda }}_1$ is $\begin{array}{rl}det{\mathrm{\Lambda }}_1& ={\mu }_1\left(d_0+2d_{1}J+\frac{\alpha }{(1+J{)}^2}\right)-\frac{\alpha }{(1+J{)}^2}\frac{\beta A}{\gamma +A}\left(1+\frac{\gamma }{\gamma +A}\right)\\ & ={\mu }_1\left(d_0+d_{1}J+\frac{\alpha }{(1+J{)}^2}+d_{1}J\right)-\frac{{\mu }_{1}A}{(1+J)J}\frac{\beta A}{\gamma +A}\left(1+\frac{\gamma }{\gamma +A}\right)\\ & =\frac{{\mu }_1}{(1+J)J}\left((d_0+d_{1}J)J(1+J)+\frac{\alpha J}{1+J}+d_1(1+J)J^2\right)\\ & -\frac{{\mu }_{1}A}{(1+J)J}\frac{\beta A}{\gamma +A}\left(1+\frac{\gamma }{\gamma +A}\right).\end{array}$ It follows from Equation (5(5) $\begin{array}{rl}\frac{\beta A}{\gamma +A}A& =\frac{\alpha J}{1+J}+(d_0+d_{1}J)J,\\ \frac{\alpha J}{1+J}& ={\mu }_{1}A,\end{array}$(5) ) that (11) $\begin{array}{rl}det{\mathrm{\Lambda }}_1& =\frac{{\mu }_1}{(1+J)J}\left(\frac{\beta A}{\gamma +A}A+(d_0+d_{1}J+d_1(1+J\left)\right)J^2\right)\\ & -\frac{{\mu }_{1}A}{(1+J)J}\frac{\beta A}{\gamma +A}\left(1+\frac{\gamma }{\gamma +A}\right)\\ & =\frac{{\mu }_1}{(1+J)J}\left((d_0+d_{1}J+d_1(1+J\left)\right)J^2-\frac{\beta \gamma A^2}{(\gamma +A{)}^2}\right),\end{array}$(11) and from Equation (5(5) $\begin{array}{rl}\frac{\beta A}{\gamma +A}A& =\frac{\alpha J}{1+J}+(d_0+d_{1}J)J,\\ \frac{\alpha J}{1+J}& ={\mu }_{1}A,\end{array}$(5) ) again that (12) $\frac{\beta \gamma A^2}{(\gamma +A{)}^2}=\frac{\beta \gamma {\alpha }^2{J}^2}{(\gamma {\mu }_1+(\gamma {\mu }_1+\alpha )J{)}^2}.$(12) Substituting Equation (12(12) $\frac{\beta \gamma A^2}{(\gamma +A{)}^2}=\frac{\beta \gamma {\alpha }^2{J}^2}{(\gamma {\mu }_1+(\gamma {\mu }_1+\alpha )J{)}^2}.$(12) ) into Equation (11(11) $\begin{array}{rl}det{\mathrm{\Lambda }}_1& =\frac{{\mu }_1}{(1+J)J}\left(\frac{\beta A}{\gamma +A}A+(d_0+d_{1}J+d_1(1+J\left)\right)J^2\right)\\ & -\frac{{\mu }_{1}A}{(1+J)J}\frac{\beta A}{\gamma +A}\left(1+\frac{\gamma }{\gamma +A}\right)\\ & =\frac{{\mu }_1}{(1+J)J}\left((d_0+d_{1}J+d_1(1+J\left)\right)J^2-\frac{\beta \gamma A^2}{(\gamma +A{)}^2}\right),\end{array}$(11) ) yields $\frac{(1+J)J}{{\mu }_1{J}^2}det{\mathrm{\Lambda }}_1=d_0+d_{1}J+d_1(1+J)-\frac{\beta \gamma {\alpha }^2}{(\gamma {\mu }_1+(\gamma {\mu }_1+\alpha )J{)}^2}.$ It follows from ${\mathrm{\Phi }}^{\prime }\left(J\right)=\frac{{\mu }_1\beta \gamma {\alpha }^2}{(\gamma {\mu }_1+(\gamma {\mu }_1+\alpha )J{)}^2}-{\mu }_1(d_0+d_1+2d_{1}J)$ that $\frac{(1+J)J}{{\mu }_1{J}^2}det{\mathrm{\Lambda }}_1=-\frac{{\mathrm{\Phi }}^{\prime }\left(J\right)}{{\mu }_1}.$ If there exist two positive equilibria of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ), $E_1=(J_1,A_1)$ and $E_2=(J_2,A_2)$, then ${\mathrm{\Psi }}^{\prime }\left(J_1\right)<0$ and ${\mathrm{\Psi }}^{\prime }\left(J_2\right)>0$. Since ${\mathrm{\Phi }}^{\prime }\left(J\right)=-\frac{{\mathrm{\Psi }}^{\prime }\left(J\right)}{\gamma {\mu }_1+(\gamma {\mu }_1+\alpha )J}$ at a positive equilibrium, we have ${\mathrm{\Phi }}^{\prime }\left(J_1\right)>0$ and ${\mathrm{\Phi }}^{\prime }\left(J_2\right)<0$, and hence $det{\mathrm{\Lambda }}_1\left(E_1\right)<0,det{\mathrm{\Lambda }}_1\left(E_2\right)>0.$ 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(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) has no closed orbits. We summarize the results for system (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) as follows.

Theorem 2.2.

System (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) has a trivial equilibrium $(0,0)$ which is always locally asymptotically stable. The existence of positive equilibria for system (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) is summarized in the following table.

Table

${\delta }_0\le 1$	${\delta }_0>1$
No positive equilibrium	${\mathrm{\Delta }}_0>1$	${\mathrm{\Delta }}_0=1$	${\mathrm{\Delta }}_0<1$
	No positive equilibrium	Unique positive equilibrium	2 positive equilibria

Here ${\delta }_0$ is given in Equation (7(7) ${\delta }_0:=\frac{\beta {\alpha }^2}{{\mu }_1\left(\gamma {\mu }_1\right(d_0+d_1)+(\gamma {\mu }_1+\alpha \left)\right(d_0+\alpha \left)\right)}.$(7) ), and ${\mathrm{\Delta }}_0$ is given in Equation (10(10) ${\mathrm{\Delta }}_0:=\frac{27K_1^2\gamma {\mu }_1^2(d_0+\alpha )}{2(K_2^2-3K_1{K}_3)(\sqrt{K_2^2-3K_1{K}_3}-K_2)+3K_1{K}_2{K}_3}.$(10) ). If there exist two positive equilibria $E_1=(J_1,A_1)$ and $E_2=(J_2,A_2)$ with $J_1<J_0<J_2,$ where $J_0$ is given in Equation (9(9) $J_0:=\frac{\sqrt{K_2^2-3K_1{K}_3}-K_2}{3K_1}.$(9) ), then $E_1$ is an unstable saddle posit and $E_2$ is a locally asymptotically node or spiral. System (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ) has no closed orbits.

3. Mosquito population model with sterile mosquitoes

We now consider that sterile mosquitoes are released to a wild mosquito population, and let $g\left(t\right)$ 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(\cdot )$. 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 ${\mu }_2$. The dynamics of sterile mosquitoes are then determined by the following simple equation: (13) $\frac{\mathrm{d}g}{\mathrm{d}t}=B(\cdot )-{\mu }_{2}g.$(13)

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 (14) $\beta =\frac{\sigma A}{A+g},$(14) where σ is the number of offspring produced per wild mosquito through all matings, per unit of time. Based on systems (1(1) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\beta (\cdot )A-\alpha (1-\kappa (J\left)\right)J-(d_0+d_{1}J)J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}& =\alpha (1-\kappa (J\left)\right)J-{\mu }_{1}A.\end{array}$(1) ), (13(13) $\frac{\mathrm{d}g}{\mathrm{d}t}=B(\cdot )-{\mu }_{2}g.$(13) ), and (14(14) $\beta =\frac{\sigma A}{A+g},$(14) ), the dynamics of the interactive wild and sterile mosquitoes are described by the following system: (15) $\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}$(15)

It is easy to see that the first quadrant is positive invariant under the flow of system (15(15) $\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}$(15) ) and the intrinsic growth rate of the wild mosquito population in the absence of sterile mosquitoes is (16) $r_0:=\frac{\sigma \alpha }{(\alpha +d_0){\mu }_1}.$(16)

In the case where Allee effects are included, based on systems (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ), (13(13) $\frac{\mathrm{d}g}{\mathrm{d}t}=B(\cdot )-{\mu }_{2}g.$(13) ), and (14(14) $\beta =\frac{\sigma A}{A+g},$(14) ), the dynamics of the interactive wild and sterile mosquitoes are described by the following system: (17) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}$(17)

We consider three different cases with different release strategies as in [9, 18] in the following sections.

4. 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(\cdot ):=b$. Then there is no difficulty for mosquitoes to find mates and based on system (15(15) $\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}$(15) ), the interactive dynamics are governed by the following system: (18) $\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-{\mu }_{2}g.\end{array}$(18)

There exists a trivial boundary equilibrium $E_0:=(J,A,g)=(0,0,g^0)$, where $g^0=b/{\mu }_2$, for system (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ). Simple linear stability analysis shows that it is locally asymptotically stable.

At a positive equilibrium, we have $g=b/{\mu }_2$. Substituting it into the first equation for γ in Equation (5(5) $\begin{array}{rl}\frac{\beta A}{\gamma +A}A& =\frac{\alpha J}{1+J}+(d_0+d_{1}J)J,\\ \frac{\alpha J}{1+J}& ={\mu }_{1}A,\end{array}$(5) ) in Section 2.2, we have the following equation to determine the positive equilibria: ${\mathrm{\Phi }}_c\left(J\right)=-\frac{{\mathrm{\Psi }}_c\left(J\right)}{{\mu }_{1}b+({\mu }_{1}b+{\mu }_2\alpha )J},$ where (19) ${\mathrm{\Psi }}_c\left(J\right)=B_1{J}^3+B_2{J}^2+B_{3}J+{\mu }_1^2(d_0+\alpha )b,$(19) and (20) $\begin{array}{rl}{B}_1& :={\mu }_1{d}_1({\mu }_{1}b+{\mu }_2\alpha ),\\ B_2& :={\mu }_1\left({\mu }_1\right(d_0+2d_1)b+\alpha {\mu }_2(d_0+d_1\left)\right),\\ B_3& :={\mu }_1\left({\mu }_1\right(2d_0+d_1+\alpha )b+{\mu }_2\alpha (d_0+\alpha \left)\right)-{\mu }_2\sigma {\alpha }^2={\mu }_1^2(2d_0+d_1+\alpha )(b-b_c).\end{array}$(20) Here we define a threshold value of releases as (21) $b_c:=\frac{\alpha {\mu }_2(\sigma \alpha -(\alpha +d_0\left){\mu }_1\right)}{{\mu }_1^2(\alpha +2d_0+d_1)}=\frac{\alpha {\mu }_2(\alpha +d_0)(r_0-1)}{{\mu }_1(\alpha +2d_0+d_1)}.$(21) Clearly, if $r_0\le 1$, there exists no positive equilibrium for system (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ) for any $b\ge 0$. If $r_0>1$, there exists no positive equilibrium for system (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ) if $b\ge b_c>0$ since then $B_3>0$.

Assume $r_0>1$ and $b<b_c$. Similarly as in Section 2.2, if we define (22) ${\mathrm{\Delta }}_c:=\frac{27{\mu }_1^2(d_0+\alpha )B_1^{2}b}{2(B_2^2-3B_1{B}_3)(\sqrt{B_2^2-3B_1{B}_3}-B_2)+3B_1{B}_2{B}_3},$(22) then there exist no, one, or two positive equilibria for system (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ) if ${\mathrm{\Delta }}_c>1$, ${\mathrm{\Delta }}_c=1$, or ${\mathrm{\Delta }}_c<1$, respectively.

The expression of ${\mathrm{\Delta }}_c$ is depending on b. We can further determine a threshold value $b_s^c$ for the existence of positive equilibria as follows.

Let function $G_c\left(b\right)$ be defined as $G_c\left(b\right):={\mathrm{\Psi }}_c\left(J_c\right(b),b)=B_1\left(b\right)J_c^3\left(b\right)+B_2\left(b\right)J_c^2\left(b\right)+B_3\left(b\right)J_c\left(b\right)+{\mu }_1^2(d_0+\alpha )b,$ where ${\mathrm{\Psi }}_c$ is given in Equation (19(19) ${\mathrm{\Psi }}_c\left(J\right)=B_1{J}^3+B_2{J}^2+B_{3}J+{\mu }_1^2(d_0+\alpha )b,$(19) ), and $J_c$ satisfies ${\mathrm{\Psi }}_c^{\prime }\left(J\right)=0$ and is given as follows: (23) $J_c:=\frac{\sqrt{B_2^2-3B_1{B}_3}-B_2}{3B_1}.$(23) Since ${\mathrm{\Psi }}_c^{\prime }\left(J_c\right)=0$ and $B_i^{\prime }\left(b\right)>0$, i=1,2,3, it follows from $\begin{array}{rl}{G}_c^{\prime }\left(b\right)& ={\mathrm{\partial }}_{J_c}{\mathrm{\Psi }}_c{J}_c^{\prime }\left(b\right)+{\mathrm{\partial }}_b{\mathrm{\Psi }}_c\\ & ={\mathrm{\partial }}_b{\mathrm{\Psi }}_c=B_1^{\prime }\left(b\right)J_c^3\left(b\right)+B_2^{\prime }\left(b\right)J_c^2\left(b\right)+B_3^{\prime }\left(b\right)J_c\left(b\right)+{\mu }_1^2(d_0+\alpha )>0\end{array}$ that ${\mathrm{\Psi }}_c\left(J_c\right(b\left)\right)$ is an increasing function of b, which implies that ${\mathrm{\Delta }}_c\left(b\right)$ is monotone increasing as b increases. Therefore, there exists a threshold value $b_s^c>0$ such that as $b>b_s^c$, $b=b_s^c$, or $b<b_s^c$, ${\mathrm{\Delta }}_c>1$, ${\mathrm{\Delta }}_c=1$, or ${\mathrm{\Delta }}_c<1$, that is, system (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ) has no, one, or two positive equilibria, respectively.

For the asymptotic dynamics of system (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ), it is easy to see that the JS-plane is a global attractor for system (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ), and the stability analysis for Equation (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ) is similar to that for Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\beta A}{\gamma +A}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.\end{array}$(4) ). The results are summarized as follows.

Theorem 4.1.

System (18(18) $\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-{\mu }_{2}g.\end{array}$(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(18) $\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-{\mu }_{2}g.\end{array}$(18) ) are summarized in the following table.

Table

$r_0\le 1$	$r_0>1$
	$b\ge b_c$	$b<b_c$
		$b>b_s^c$	$b=b_s^c$	$b<b_s^c$
No positive equilibrium	No positive equilibrium	No positive equilibrium	Unique positive equilibrium	2 positive equilibria

Here $r_0$ is the intrinsic growth rate of the wild mosquito population in the absence of sterile mosquitoes given in Equation (16(16) $r_0:=\frac{\sigma \alpha }{(\alpha +d_0){\mu }_1}.$(16) ), the threshold value of the existence of positive equilibria $b_c$ is given in Equation (21(21) $b_c:=\frac{\alpha {\mu }_2(\sigma \alpha -(\alpha +d_0\left){\mu }_1\right)}{{\mu }_1^2(\alpha +2d_0+d_1)}=\frac{\alpha {\mu }_2(\alpha +d_0)(r_0-1)}{{\mu }_1(\alpha +2d_0+d_1)}.$(21) ), and the threshold value $b_s^c$ is determined by ${\mathrm{\Delta }}_c\left(b_s^c\right)=1$ with ${\mathrm{\Delta }}_c$ given in Equation (22(22) ${\mathrm{\Delta }}_c:=\frac{27{\mu }_1^2(d_0+\alpha )B_1^{2}b}{2(B_2^2-3B_1{B}_3)(\sqrt{B_2^2-3B_1{B}_3}-B_2)+3B_1{B}_2{B}_3},$(22) ). If there exist two positive equilibria $E_1=(J_1,A_1,g_1)$ and $E_2=(J_2,A_2,g_2)$ with $J_1<J_c<J_2,$ where $J_c$ is given in Equation (23(23) $J_c:=\frac{\sqrt{B_2^2-3B_1{B}_3}-B_2}{3B_1}.$(23) ) and $B_i,$ $i=1,2,3,$ are given in Equation (20(20) $\begin{array}{rl}{B}_1& :={\mu }_1{d}_1({\mu }_{1}b+{\mu }_2\alpha ),\\ B_2& :={\mu }_1\left({\mu }_1\right(d_0+2d_1)b+\alpha {\mu }_2(d_0+d_1\left)\right),\\ B_3& :={\mu }_1\left({\mu }_1\right(2d_0+d_1+\alpha )b+{\mu }_2\alpha (d_0+\alpha \left)\right)-{\mu }_2\sigma {\alpha }^2={\mu }_1^2(2d_0+d_1+\alpha )(b-b_c).\end{array}$(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(18) $\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-{\mu }_{2}g.\end{array}$(18) ) in Example 1.

Example 1.

Given the parameters (24) $\sigma =18,\alpha =7,d_0=0.35,d_1=0.3,{\mu }_1=0.25,{\mu }_2=0.196,$(24) we have (25) $r_0=68.5714,b_c=340.7019,b_s^c=\mathrm{183.8706.}$(25) For $b=190>b_s^c=183.8706$, there exists no positive equilibrium. All solutions go to the origin as shown in the left figure in Figure 1. For $b=170<b_s^c$, there exist two positive equilibria $E_1=(1.1213,14.8005,867.3469)$ and $E_2=(2.9408,20.8948,867.3469)$. Equilibrium $E_1$ is an unstable saddle and $E_2$ is a stable node. Solutions approach either the origin or $E_2$ depending on their initial values. We only show a solution approaching $E_2$ in the right figure in Figure 1.

Figure 1. Parameters are given in Equation (24(24) $\sigma =18,\alpha =7,d_0=0.35,d_1=0.3,{\mu }_1=0.25,{\mu }_2=0.196,$(24) ). For $b=190>b_c^s$, there exists no positive equilibrium. Solutions approach the origin as shown in the left figure. For $b=170<b_s^c$, there exist two positive equilibria $E_1$ and $E_2$. Solutions approach either the origin or $E_2$ depending on their initial values. A solution approaching $E_2$ is shown in the right figure. Here only the curves for wild larvae and adults are presented.

5. 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(\cdot )=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(17) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}$(17) ), are governed by the following system: (26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26)

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(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ), we have $g=\frac{b}{{\mu }_2}A.$ Then it follows from $\frac{\sigma A}{1+\left(1+\frac{b}{{\mu }_2}\right)A}A=\frac{{\mu }_2\sigma A}{{\mu }_2+({\mu }_2+b)A}A=\frac{\alpha J}{1+J}+(d_0+d_{1}J)J$ that the corresponding function for the existence of positive equilibria has the form $\begin{array}{rl}{\mathrm{\Phi }}_p& :=\frac{{\mu }_2\sigma {\alpha }^{2}J}{{\mu }_1{\mu }_2+\left({\mu }_2\right({\mu }_1+\alpha )+\alpha b)J}-{\mu }_1(\alpha +(d_0+d_{1}J\left)\right(1+J\left)\right)\\ & =-\frac{{\mathrm{\Psi }}_p\left(J\right)}{\left({\mu }_2\right({\mu }_1+\alpha )+\alpha b)J+{\mu }_1{\mu }_2},\end{array}$ where (27) ${\mathrm{\Psi }}_p\left(J\right):=C_1{J}^3+C_2{J}^2+C_{3}J+{\mu }_1^2{\mu }_2(d_0+\alpha ),$(27) with (28) $\begin{array}{rl}{C}_1& :={\mu }_1{d}_1(\alpha b+{\mu }_2({\mu }_1+\alpha \left)\right),\\ C_2& :={\mu }_1\left(\alpha \right(d_0+d_1)b+{\mu }_2({\mu }_1(d_0+2d_1)+\alpha (d_0+d_1)\left)\right),\\ C_3& :={\mu }_1\alpha (d_0+\alpha )(b-b_p).\end{array}$(28) Here the threshold value for the releases, $b_p$, is defined as (29) $\begin{array}{rl}{b}_p& :=\frac{{\mu }_2(\sigma {\alpha }^2-{\mu }_1({\mu }_1(2d_0+d_1+\alpha )+\alpha (d_0+\alpha )\left)\right)}{{\mu }_1\alpha (d_0+\alpha )}\\ & =\frac{{\mu }_1{\mu }_2\left({\mu }_1\right(2d_0+d_1+\alpha )+\alpha (d_0+\alpha \left)\right)(r_p-1)}{{\mu }_1\alpha (d_0+\alpha )}\end{array}$(29) and (30) $r_p:=\frac{\sigma {\alpha }^2}{{\mu }_1\left({\mu }_1\right(2d_0+d_1+\alpha )+\alpha (d_0+\alpha \left)\right)}.$(30) System (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) clearly has no positive equilibrium if $r_p\le 1$ for any $b\ge 0$, or $r_p>1$ and $b\ge b_p$.

Suppose $r_p>1$ and $b<b_p$. We define (31) ${\mathrm{\Delta }}_p:=\frac{27{\mu }_1^2{\mu }_2(d_0+\alpha )C_1^2}{2(C_2^2-3C_1{C}_3)(\sqrt{C_2^2-3C_1{C}_3}-C_2)+3C_1{C}_2{C}_3}.$(31) Then system (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) has no, one, or two positive equilibria provided ${\mathrm{\Delta }}_p>1$, ${\mathrm{\Delta }}_p=1$, or ${\mathrm{\Delta }}_p<1$, respectively.

Similarly as in Section 4, since ${\mathrm{\Delta }}_p$ depends on b, we can determine a threshold value of b for the existence of positive equilibria of system (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) as follows.

Let function $G_p\left(b\right)$ be defined as $G_p\left(b\right):={\mathrm{\Psi }}_p\left(J_p\right(b),b)=C_1\left(b\right)J_p^3\left(b\right)+C_2\left(b\right)J_p^2\left(b\right)+C_3\left(b\right)J_p\left(b\right)+{\mu }_1^2{\mu }_2(d_0+\alpha ),$ where (32) $J_p:=\frac{\sqrt{C_2^2-3C_1{C}_3}-C_2}{3C_1}.$(32) Then since ${\mathrm{\Psi }}_p^{\prime }\left(J_p\right)=0$ and $C_i^{\prime }\left(b\right)>0$, i=1,2,3, it follows from $G_p^{\prime }\left(b\right)={\mathrm{\partial }}_{J_c}{\mathrm{\Psi }}_p{J}_p^{\prime }\left(b\right)+{\mathrm{\partial }}_b{\mathrm{\Psi }}_p={\mathrm{\partial }}_b{\mathrm{\Psi }}_p=C_1^{\prime }\left(b\right)J_p^3\left(b\right)+C_2^{\prime }\left(b\right)J_p^2\left(b\right)+C_3^{\prime }\left(b\right)J_p\left(b\right)>0$ that ${\mathrm{\Psi }}_p\left(J_p\right(b\left)\right)$ is an increasing function of b, which implies that ${\mathrm{\Delta }}_p\left(b\right)$ is monotone increasing as b increases. Therefore, there exists a threshold value $b_s^p>0$ such that as $b>b_s^p$, $b=b_s^p$, or $b<b_s^p$, ${\mathrm{\Delta }}_p>1$, ${\mathrm{\Delta }}_p=1$, ${\mathrm{\Delta }}_p<1$, and hence system (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) has no, one, or two positive equilibria, respectively.

We next study the stability of the positive equilibria. The Jacobian matrix of system (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) has the form $\begin{array}{rl}{\mathrm{\Lambda }}_2& :=\left(\begin{array}{ccc}-d_0-2d_{1}J-\frac{\alpha }{(1+J{)}^2}& \frac{\sigma A}{1+A+g}\left(1+\frac{1+g}{1+A+g}\right)& -\frac{\sigma A^2}{(1+A+g{)}^2}\\ \frac{\alpha }{(1+J{)}^2}& -{\mu }_1& 0\\ 0& b& -{\mu }_2\end{array}\right)\\ & :=\left(\begin{array}{ccc}-a_{11}& a_{12}& -a_{13}\\ a_{121}& -{\mu }_1& 0\\ 0& b& -{\mu }_2\end{array}\right).\end{array}$ The characteristic equation of the linearized system of (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) at an equilibrium then is given by $H\left(\lambda \right):={\lambda }^3+P_1{\lambda }^2+P_2\lambda +P_3=0,$ where (33) $\begin{array}{rl}{P}_1& :=-\mathrm{tr}{\mathrm{\Lambda }}_2={\mu }_1+{\mu }_2+a_{11},\\ P_2& :=\sum 2\times 2\mathrm{principleminors}={\mu }_1{\mu }_2+a_{11}({\mu }_1+{\mu }_2)-a_{21}{a}_{12},\\ P_3& :=-det{\mathrm{\Lambda }}_2={\mu }_1{\mu }_2{a}_{11}+b{a}_{13}{a}_{21}-{\mu }_2{a}_{12}{a}_{21}.\end{array}$(33)

Notice that the J component of a positive equilibrium satisfies Equation (27(27) ${\mathrm{\Psi }}_p\left(J\right):=C_1{J}^3+C_2{J}^2+C_{3}J+{\mu }_1^2{\mu }_2(d_0+\alpha ),$(27) ). Similarly as in Section 2.2, if there are two positive equilibria, denoted by $E_1=(J_1,A_1,g_1)$ and $E_2=(J_2,A_2,g_2)$, where $J_1<J_p<J_2,$ and $J_p$ is given in Equation (32(32) $J_p:=\frac{\sqrt{C_2^2-3C_1{C}_3}-C_2}{3C_1}.$(32) ), then we have ${\mathrm{\Psi }}_p^{\prime }\left(J_1\right)<0$ and ${\mathrm{\Psi }}_p^{\prime }\left(J_2\right)>0$, and hence ${\mathrm{\Phi }}_p^{\prime }\left(J_1\right)>0$ and ${\mathrm{\Phi }}_p^{\prime }\left(J_2\right)<0$. By direct calculation, we have $P_3=-\frac{{\mu }_1{\mu }_{2}J}{1+J}{\mathrm{\Phi }}_p^{\prime }\left(J\right).$ (Details are given in Appendix 1). Thus $P_3\left(E_1\right)<0$ and $P_3\left(E_2\right)>0$, and since $P_3=-det{\mathrm{\Lambda }}_2$, we have $det{\mathrm{\Lambda }}_2\left(E_1\right)>0,det{\mathrm{\Lambda }}_2\left(E_2\right)<0,$ which implies that $E_1$ is unstable. Furthermore, at $E_2$, direct but tedious calculation shows that if ${\mu }_2\ge {\mu }_1$, then $P_1{P}_2-P_3>0$. (Details are given in Appendix 2.) Since $P_3\left(E_2\right)>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(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) is always locally asymptotically stable. The existence results for positive equilibria of system (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) are summarized in the following table.

Table

$r_p\le 1$	$r_p>1$
	$b\ge b_p$	$b<b_p$
		$b>b_s^p$	$b=b_s^p$	$b<b_s^p$
No positive equilibrium	No positive equilibrium	No positive equilibrium	Unique positive equilibrium	2 positive equilibria

Here the threshold values $b_p$ and $r_p$ are defined in Equations (29(29) $\begin{array}{rl}{b}_p& :=\frac{{\mu }_2(\sigma {\alpha }^2-{\mu }_1({\mu }_1(2d_0+d_1+\alpha )+\alpha (d_0+\alpha )\left)\right)}{{\mu }_1\alpha (d_0+\alpha )}\\ & =\frac{{\mu }_1{\mu }_2\left({\mu }_1\right(2d_0+d_1+\alpha )+\alpha (d_0+\alpha \left)\right)(r_p-1)}{{\mu }_1\alpha (d_0+\alpha )}\end{array}$(29) ) and (30(30) $r_p:=\frac{\sigma {\alpha }^2}{{\mu }_1\left({\mu }_1\right(2d_0+d_1+\alpha )+\alpha (d_0+\alpha \left)\right)}.$(30) ), respectively. The threshold value $b_s^p$ is determined by ${\mathrm{\Delta }}_p\left(b\right)=1$ with ${\mathrm{\Delta }}_p$ given (31(31) ${\mathrm{\Delta }}_p:=\frac{27{\mu }_1^2{\mu }_2(d_0+\alpha )C_1^2}{2(C_2^2-3C_1{C}_3)(\sqrt{C_2^2-3C_1{C}_3}-C_2)+3C_1{C}_2{C}_3}.$(31) ). If there exist two positive equilibria $E_1=(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(23) $J_c:=\frac{\sqrt{B_2^2-3B_1{B}_3}-B_2}{3B_1}.$(23) ), equilibrium $E_1$ is an unstable saddle posit and $E_2$ is a locally asymptotically node or spiral, provided ${\mu }_2\ge {\mu }_1$.

The dynamics of system (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) have been determined in Theorem 5.1 for ${\mu }_2\ge {\mu }_1$. However, the model dynamics can be complex if ${\mu }_2<{\mu }_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(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) in Example 2.

Example 2.

Let parameters be given by (34) $\sigma =18,\alpha =7,d_0=0.35,d_1=0.3,{\mu }_1=0.25,{\mu }_2=\mathrm{0.196.}$(34) Note that ${\mu }_2=0.196<{\mu }_1=0.25$. It follows from Equations (29(29) $\begin{array}{rl}{b}_p& :=\frac{{\mu }_2(\sigma {\alpha }^2-{\mu }_1({\mu }_1(2d_0+d_1+\alpha )+\alpha (d_0+\alpha )\left)\right)}{{\mu }_1\alpha (d_0+\alpha )}\\ & =\frac{{\mu }_1{\mu }_2\left({\mu }_1\right(2d_0+d_1+\alpha )+\alpha (d_0+\alpha \left)\right)(r_p-1)}{{\mu }_1\alpha (d_0+\alpha )}\end{array}$(29) ) and (30(30) $r_p:=\frac{\sigma {\alpha }^2}{{\mu }_1\left({\mu }_1\right(2d_0+d_1+\alpha )+\alpha (d_0+\alpha \left)\right)}.$(30) ) that $b_p=13.2364,r_p=66.0056,$ and then $b_s^p=13.0521$.

For $b=12.57<b_s^p$, ${\mathrm{\Delta }}_p=0.0797<1$. There exist two positive equilibria $E_1=(0.010699,0.29641,19.00956)$ and $E_2=(0.47332,8.99533,576.89438)$. Equilibrium $E_1$ is an unstable saddle and $E_2$ is a stable spiral. There exists an unstable closed orbit as shown in the left figure in 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.

Figure 2. Parameters are given in Equation (34(34) $\sigma =18,\alpha =7,d_0=0.35,d_1=0.3,{\mu }_1=0.25,{\mu }_2=\mathrm{0.196.}$(34) ). With b=12.57, there exist positive equilibrium $E_1$ which is an unstable saddle and positive equilibrium $E_2$ which is a stable spiral. There exists a unstable closed orbit as shown in the left figure. Solutions initially started near $E_2$ spiral towards $E_2$ and solutions initially started away from the closed orbit approach the origin as both shown in the right figure. Here only the curves for the wild adults are presented. Notice that the speed of the convergence to $E_2$ is slow.

With the same parameters but b=12.5724 which is still less than $b_s^p$, 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.

Figure 3. With b=12.57, there are two positive equilibria $E_1$ which is an unstable saddle and $E_2$ which is an stable spiral. There exists a unstable closed orbit as shown in the left figure. Sustained oscillations occur. Solutions initially started near $E_2$ spiral towards the bistable closed orbit and solutions initially started away from the closed orbit still approach the origin as both shown in the right figure. Here also only the curves for the wild adults are presented.

6. 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(\cdot )=bA/(1+A)$. Then the interactive system becomes (35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35)

The origin $(0,0,0)$ is a trivial equilibrium and is clearly locally asymptotically stable. At a positive equilibrium, we have $J=\frac{A}{\rho -A},0<A<\rho ,g=\frac{bA}{{\mu }_2(1+A)},$ where $\rho =\alpha /{\mu }_1$. Substituting them into the first equation for positive equilibria in Equation (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ), we have $\frac{{\mu }_2\sigma A^2(1+A)}{{\mu }_2(1+A{)}^2+bA}=\frac{\alpha J}{1+J}+(d_0+d_{1}J)J.$ Then the equation to determine positive equilibria of Equation (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ) is equivalent to (36) $\begin{array}{rl}{\mathrm{\Phi }}_H\left(A\right)& :=\frac{{\mu }_2\sigma A(1+A)}{{\mu }_2(1+A{)}^2+bA}-{\mu }_1-\frac{\rho d_0+(d_1-d_0)A}{(\rho -A{)}^2}\\ & =\sigma -{\mu }_1-\sigma \frac{{\mu }_2+({\mu }_2+b)A}{{\mu }_2(1+A{)}^2+bA}-\frac{\rho d_0+(d_1-d_0)A}{(\rho -A{)}^2}=0.\end{array}$(36)

Simple calculation yields $\begin{array}{rl}{\mathrm{\Phi }}_H^{\prime }\left(A\right)& =\sigma \frac{{\mu }_2^2+2{\mu }_2^{2}A+{\mu }_2({\mu }_2+b)A^2}{\left({\mu }_2\right(1+A{)}^2+bA{)}^2}-\frac{(d_1-d_0)(\rho -A)+2(d_{1}A+d_0(\rho -A\left)\right)}{(\rho -A{)}^3}\\ & =\sigma \frac{{\mu }_2^2+2{\mu }_2^{2}A+{\mu }_2({\mu }_2+b)A^2}{\left({\mu }_2\right(1+A{)}^2+bA{)}^2}-\frac{(d_0+d_1)\rho +(d_1-d_0)A}{(\rho -A{)}^3}\end{array}$ and ${\mathrm{\Phi }}_H^{″}\left(A\right)=-2{\mu }_2^2\sigma \frac{b(1+A^3)+{\mu }_2(1+A{)}^3}{\left({\mu }_2\right(1+A{)}^2+bA{)}^3}-2\frac{(d_0+2d_1)\rho +(d_1-d_0)A}{(\rho -A{)}^4}.$ It follows from $0<A<\rho$ that $(d_0+2d_1)\rho +(d_1-d_0)A>3d_{1}A>0.$ Hence ${\mathrm{\Phi }}_H^{″}\left(A\right)<0$ for $0\le A<\rho$, and there exists a unique $A_H$, with $0<A_H<\rho$, such that ${\mathrm{\Phi }}_H^{\prime }\left(A_H\right)=0$.

Notice that ${\mathrm{\Phi }}_H\left(0\right)=-{\mu }_1-d_0/\rho$ and $\underset{A\to \rho -}{lim}{\mathrm{\Phi }}_H\left(A\right)=-\mathrm{\infty }$, ${\mathrm{\Phi }}_H\left(A_H\right)=\underset{0<A<\rho }{max}{\mathrm{\Phi }}_H\left(A\right)$. Thus ${\mathrm{\Phi }}_H\left(A\right)$ has no, one, or two positive roots, for $A<\rho$, if ${\mathrm{\Phi }}_H\left(A_H\right)<0$, ${\mathrm{\Phi }}_H\left(A_H\right)=0$, or ${\mathrm{\Phi }}_H\left(A_H\right)>0$, respectively.

Similarly as in Sections 4 and 5, we define $G_H\left(b\right):={\mathrm{\Phi }}_H\left(A_H\right(b),b)$. Then it follows from $G_H^{\prime }\left(b\right)={\mathrm{\partial }}_{A_H}{G}_H\left(A_H\right)A_H^{\prime }\left(b\right)+{\mathrm{\partial }}_b{G}_H\left(b\right)={\mathrm{\partial }}_b{G}_H\left(b\right)=-\frac{{\mu }_2\sigma A^2(1+A)}{\left({\mu }_2\right(1+A{)}^2+bA{)}^2}<0$ that $\mathrm{\Phi }\left(A_H\right)$ is a decreasing function of b and there exists a unique threshold value $b_s^H>0$ such that then if $b>b_s^H$, $b=b_s^H$, or $b<b_s^H$, then ${\mathrm{\Phi }}_H\left(A_H\right)<0$, ${\mathrm{\Phi }}_H\left(A_H\right)=0$, or ${\mathrm{\Phi }}_H\left(A_H\right)>0$, and hence system (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ) has no, one, or two positive equilibria, respectively.

We next investigate the stability of the positive equilibria.

At a positive equilibrium, the Jacobian matrix of system (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ) has the following form: $\begin{array}{rl}{\mathrm{\Lambda }}_3& :=\left(\begin{array}{ccc}-d_0-2d_{1}J-\frac{\alpha }{(1+J{)}^2}& \frac{\sigma A}{1+A+g}\left(1+\frac{1+g}{1+A+g}\right)& -\frac{\sigma A^2}{(1+A+g{)}^2}\\ \frac{\alpha }{(1+J{)}^2}& -{\mu }_1& 0\\ 0& \frac{b}{(1+A{)}^2}& -{\mu }_2\end{array}\right)\\ & :=\left(\begin{array}{ccc}-c_{11}& c_{12}& -c_{13}\\ c_{21}& -{\mu }_1& 0\\ 0& c_{32}& -{\mu }_2\end{array}\right),\end{array}$ and hence the characteristic equation is ${\lambda }^3+Q_1{\lambda }^2+Q_2\lambda +Q_3=0,$ where (37) $\begin{array}{rl}{Q}_1& ={\mu }_1+{\mu }_2+c_{11},\\ Q_2& ={\mu }_1{\mu }_2+c_{11}({\mu }_1+{\mu }_2)-c_{21}{c}_{12},\\ Q_3& ={\mu }_1{\mu }_2{c}_{11}+c_{32}{c}_{13}{c}_{21}-{\mu }_2{c}_{12}{c}_{21}.\end{array}$(37)

Similarly as in Section 5, direct calculation yields $Q_3=-\frac{{\mu }_1{\mu }_2{A}^2}{(1+J)J}{\mathrm{\Phi }}_H^{\prime }\left(A\right).$ (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 ${\mathrm{\Phi }}_H^{\prime }\left(E_1\right)>0$ and ${\mathrm{\Phi }}_H^{\prime }\left(E_2\right)<0$ that $Q_3\left(E_1\right)<0$ and $Q_3\left(E_2\right)>0$. Hence $det{\mathrm{\Lambda }}_3\left(E_1\right)>0,det{\mathrm{\Lambda }}_3\left(E_2\right)<0,$ which implies the instability of $E_1$. At $E_2$, direct but tedious calculation shows that if ${\mu }_2\ge {\mu }_1$, then $Q_1{Q}_2-Q_3>0$. (Details are given in Appendix 4.) Since $Q_3\left(E_2\right)>0$, positive equilibrium $E_2$ is locally asymptotically stable. The results are summarized as follows.

Theorem 6.1.

System (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ) has a trivial equilibrium $(0,0,0)$ which is always locally asymptotically stable. There exists a threshold value $b_s^H$ determined by ${\mathrm{\Phi }}_H\left(A\right)=0,$ where ${\mathrm{\Phi }}_H\left(A\right)$ is given in Equation (36(36) $\begin{array}{rl}{\mathrm{\Phi }}_H\left(A\right)& :=\frac{{\mu }_2\sigma A(1+A)}{{\mu }_2(1+A{)}^2+bA}-{\mu }_1-\frac{\rho d_0+(d_1-d_0)A}{(\rho -A{)}^2}\\ & =\sigma -{\mu }_1-\sigma \frac{{\mu }_2+({\mu }_2+b)A}{{\mu }_2(1+A{)}^2+bA}-\frac{\rho d_0+(d_1-d_0)A}{(\rho -A{)}^2}=0.\end{array}$(36) ), such that system (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ) has no, one, or two positive equilibria if $b>b_s^H,$ $b=b_s^H,$ or $b<b_s^H,$ for $A<\alpha /{\mu }_1,$ respectively. If there exist two positive equilibria $E_1=(J_1,A_1,g_1)$ and $E_2=(J_2,A_2,g_2),$ with $A_1<A_2,$ equilibrium $E_1$ is an unstable saddle point and $E_2$ is locally asymptotically provided ${\mu }_2\ge {\mu }_1$.

We note that the stability condition ${\mu }_2\ge {\mu }_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.

7. 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–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(18) $\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-{\mu }_{2}g.\end{array}$(18) ), (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ), and (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(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(18) $\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-{\mu }_{2}g.\end{array}$(18) ) and (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ). The threshold value for system (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ) was also implicitly obtained. We completely determined the local stability conditions of the equilibria for system (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ), and sufficient stability conditions for the positive equilibria for systems (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) and (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ) were obtained. The dynamics of system (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) and (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ) are complicated. The existence of stable or unstable closed orbits for system (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) is shown in Example 2. While we have obtained sufficient conditions for the stability of the positive equilibria for system (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(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.

In principle. the dynamics of the three model systems with different release strategies in this paper are similar to the dynamics of the model systems in [9] where no mosquito metamorphic stages are included. That is, based on threshold release values and under certain conditions, there exist two positive equilibria, one of which is unstable and one of which is locally asymptotically stable for all of the three model systems. System (18(18) $\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-{\mu }_{2}g.\end{array}$(18) ) has a boundary equilibrium $(0,0,g^0)$ and systems (26(26) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =bA-{\mu }_{2}g.\end{array}$(26) ) and (35(35) $\begin{array}{rl}\frac{\mathrm{d}J}{\mathrm{d}t}& =\frac{\sigma A}{1+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}& =\frac{bA}{1+A}-{\mu }_{2}g.\end{array}$(35) ) have the trivial equilibrium $(0,0,0)$ other than the positive equilibria. Solutions approach either the boundary (or the trivial) equilibrium, or the stable positive equilibrium, depending on their initial values. The comparisons between the three release strategies are also similar to those in [9]. Nevertheless, the analysis is more difficult for the three-dimensional systems in this paper than the two-dimensional systems in [9]. As is illustrated above, the inclusion of the mosquitoes' metamorphic stages and the density dependence in larvae's death rates is necessary from the modelling perspective. On the other hand, we seem to have also learned from this study, as has been well described in many other studies as well, that simplified models may not necessarily lose key features that the more complicated models exhibit. Therefore, it would probably be always a good idea to consider to start with relatively simple models when we work on new real world problems.

Acknowledgments

The authors thank two anonymous reviewers for their careful reading and 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 [DMS-1118150], the National Natural Science Foundation of China [11371305], China Scholarship Council [201308410212] and Nanhu Scholars Program for Young Scholars XYNU. A part of the work was performed while Liming Cai was visiting the University of Alabama in Huntsville.

Appendix 1. Calculation of P₃

(38) $\begin{array}{rl}{P}_3& ={\mu }_1{\mu }_2{a}_{11}+b{a}_{13}{a}_{21}-{\mu }_2{a}_{12}{a}_{21}\\ & ={\mu }_1{\mu }_2\left(d_0+2d_{1}J+\frac{\alpha }{(1+J{)}^2}\right)+\frac{b\sigma A^2}{(1+A+g{)}^2}\frac{\alpha }{(1+J{)}^2}\\ & -\frac{{\mu }_2\sigma A}{1+A+g}\left(1+\frac{1+g}{1+A+g}\right)\frac{\alpha }{(1+J{)}^2}\\ & =\frac{{\mu }_1{\mu }_2}{(1+J)J}\left((d_0+d_{1}J)J(1+J)+\frac{\alpha J}{(1+J)}+d_1(1+J)J^2\right)+{\mu }_1{\mu }_{2}g\frac{\sigma A^2}{(1+A+g{)}^2}\frac{1}{(1+J)J}\\ & -{\mu }_2\frac{\sigma A}{1+A+g}\left(1+\frac{1+g}{1+A+g}\right)\frac{{\mu }_{1}A}{(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2}{(1+J)J}\left(\frac{\sigma A^2}{1+A+g}+(d_0+d_{1}J+d_1(1+J\left)\right)J^2\right)+{\mu }_1{\mu }_{2}g\frac{\sigma A^2}{(1+A+g{)}^2}\frac{1}{(1+J)J}\\ & -{\mu }_2\frac{\sigma A}{1+A+g}\left(1+\frac{1+g}{1+A+g}\right)\frac{{\mu }_{1}A}{(1+J)J}\\ & =\frac{{\mu }_1{\mu }_{2}J}{1+J}(d_0+d_1+2d_{1}J)+{\mu }_1{\mu }_{2}g\frac{\sigma A^2}{(1+A+g{)}^2}\frac{1}{(1+J)J}-{\mu }_1{\mu }_2\frac{\sigma A^2(1+g)}{(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_{2}J}{1+J}(d_0+d_1+2d_{1}J)-\frac{{\mu }_1{\mu }_2\sigma A^2}{(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_{2}J}{1+J}\left(d_0+d_1+2d_{1}J-\frac{\sigma {\alpha }^2{\mu }_2^2}{({\mu }_1{\mu }_2+({\mu }_1{\mu }_2+b\alpha +\alpha {\mu }_2)J{)}^2}\right)\\ & =-\frac{{\mu }_1{\mu }_{2}J}{1+J}{\mathrm{\Phi }}_p^{\prime }\left(J\right).\end{array}$(38)

Appendix 2. Proof of P₁P₂ - P₃ > 0

It follows from Equation (33(33) $\begin{array}{rl}{P}_1& :=-\mathrm{tr}{\mathrm{\Lambda }}_2={\mu }_1+{\mu }_2+a_{11},\\ P_2& :=\sum 2\times 2\mathrm{principleminors}={\mu }_1{\mu }_2+a_{11}({\mu }_1+{\mu }_2)-a_{21}{a}_{12},\\ P_3& :=-det{\mathrm{\Lambda }}_2={\mu }_1{\mu }_2{a}_{11}+b{a}_{13}{a}_{21}-{\mu }_2{a}_{12}{a}_{21}.\end{array}$(33) ) that $\begin{array}{rl}{P}_1{P}_2-P_3& =({\mu }_1+{\mu }_2+a_{11})({\mu }_1{\mu }_2+a_{11}({\mu }_1+{\mu }_2)-a_{21}{a}_{12})-{\mu }_1{\mu }_2{a}_{11}-b{a}_{13}{a}_{21}+{\mu }_2{a}_{12}{a}_{21}\\ & =({\mu }_1+{\mu }_2)({\mu }_1{\mu }_2+a_{11}({\mu }_1+{\mu }_2\left)\right)+a_{11}^2({\mu }_1+{\mu }_2)-({\mu }_1+a_{11})a_{21}{a}_{12}-b{a}_{13}{a}_{21}\\ & ={\mu }_1{\mu }_2({\mu }_1+{\mu }_2)+({\mu }_1+{\mu }_2{)}^2{a}_{11}+a_{11}(({\mu }_1+{\mu }_2)a_{11}-a_{12}{a}_{21})-{\mu }_1{a}_{12}{a}_{21}\\ & -b{a}_{13}{a}_{21}.\end{array}$ Write ${\hat{a}}_{11}:=\frac{\sigma A^2}{(1+A+g)(1+J)J},{\overline{a}}_{11}:=\frac{J}{1+J}(d_0+d_1+2d_{1}J).$ Then $a_{11}=\frac{1}{(1+J)J}\left(\frac{\sigma A^2}{1+A+g}+(d_0+d_1+2d_{1}J)J^2\right)={\hat{a}}_{11}+{\overline{a}}_{11},$ and we have $\begin{array}{rl}{P}_1{P}_2-P_3& ={\mu }_1{\mu }_2({\mu }_1+{\mu }_2)+({\mu }_1+{\mu }_2{)}^2{\overline{a}}_{11}+{\mu }_1(({\mu }_1+{\mu }_2){\hat{a}}_{11}-a_{12}{a}_{21})\\ & +{\mu }_2({\mu }_1+{\mu }_2){\hat{a}}_{11}-b{a}_{13}{a}_{21}+a_{11}\left(\right({\mu }_1+{\mu }_2)a_{11}-a_{12}{a}_{21}).\end{array}$ It follows from, at a positive equilibrium, $a_{21}=\frac{{\mu }_{1}A}{(1+J)J},b=\frac{{\mu }_{2}g}{A}$ that if we assume ${\mu }_2\ge {\mu }_1$, then $\begin{array}{rl}& {\mu }_1\left(\right({\mu }_1+{\mu }_2){\hat{a}}_{11}-a_{12}{a}_{21})+{\mu }_2({\mu }_1+{\mu }_2){\hat{a}}_{11}-b{a}_{13}{a}_{21}\\ & =\frac{{\mu }_1({\mu }_1+{\mu }_2)\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1^2\sigma A^2}{(1+A+g)(1+J)J}\left(1+\frac{1+g}{1+A+g}\right)\\ & +\frac{{\mu }_2({\mu }_1+{\mu }_2)\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1{\mu }_2\sigma A^{2}g}{(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1^2\sigma A^2(1+g)}{(1+A+g{)}^2(1+J)J}\\ & +\frac{{\mu }_2({\mu }_1+{\mu }_2)\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1{\mu }_2\sigma A^{2}g}{(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2\sigma A^2}{(1+A+g)(1+J)J}\left(1-\frac{g}{1+A+g}\right)\\ & +\frac{\sigma A^2}{(1+A+g)(1+J)J}\left({\mu }_2^2+{\mu }_1\left({\mu }_2-{\mu }_1\frac{1+g}{1+A+g}\right)\right)>0\end{array}$ and (39) $\begin{array}{rl}({\mu }_1+{\mu }_2)a_{11}-a_{12}{a}_{21}& =({\mu }_1+{\mu }_2){\hat{a}}_{11}-a_{12}{a}_{21}+({\mu }_1+{\mu }_2){\overline{a}}_{11}\\ & =\frac{({\mu }_1+{\mu }_2)\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1\sigma A^2}{(1+A+g)(1+J)J}\\ & \cdot \left(1+\frac{1+g}{1+A+g}\right)+({\mu }_1+{\mu }_2){\overline{a}}_{11}\\ & =\frac{\sigma A^2}{(1+A+g)(1+J)J}\left({\mu }_2-{\mu }_1\frac{1+g}{1+A+g}\right)+({\mu }_1+{\mu }_2){\overline{a}}_{11}>0.\end{array}$(39) Thus $P_1{P}_2-P_3>0$, if ${\mu }_2\ge {\mu }_1$.

Appendix 3. Calculation of Q₃

It follows from, at a positive equilibrium, $c_{21}={\mu }_{1}A/(1+J)J$ and $b={\mu }_2(1+A)g/A$ that $c_{32}={\mu }_{2}g/A(1+A)$. Then (40) $\begin{array}{rl}{Q}_3& ={\mu }_1{\mu }_2{c}_{11}+c_{32}{c}_{13}{c}_{21}-{\mu }_2{c}_{12}{c}_{21}\\ & ={\mu }_1{\mu }_2\left(d_0+2d_{1}J+\frac{\alpha }{(1+J{)}^2}\right)+\frac{{\mu }_2\sigma Ag}{(1+A)(1+A+g{)}^2}\frac{\alpha }{(1+J{)}^2}\\ & -\frac{{\mu }_2\sigma A}{1+A+g}\left(1+\frac{1+g}{1+A+g}\right)\frac{\alpha }{(1+J{)}^2}\\ & =\frac{{\mu }_1{\mu }_2}{(1+J)J}\left((d_0+d_{1}J)J(1+J)+\frac{\alpha J}{(1+J)}+d_1(1+J)J^2\right)\\ & +\frac{{\mu }_1{\mu }_2\sigma A^{2}g}{(1+A)(1+A+g{)}^2(1+J)J}-{\mu }_2\frac{\sigma A}{1+A+g}\left(1+\frac{1+g}{1+A+g}\right)\frac{{\mu }_{1}A}{(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2}{(1+J)J}\left(\frac{\sigma A^2}{1+A+g}+(d_0+d_{1}J+d_1(1+J\left)\right)J^2\right)+\frac{{\mu }_1{\mu }_2\sigma A^{2}g}{(1+A)(1+A+g{)}^2(1+J)J}\\ & -{\mu }_2\frac{\sigma A}{1+A+g}\left(1+\frac{1+g}{1+A+g}\right)\frac{{\mu }_{1}A}{(1+J)J}\\ & =\frac{{\mu }_1{\mu }_{2}J}{1+J}(d_0+d_1+2d_{1}J)+\frac{{\mu }_1{\mu }_2\sigma A^{2}g}{(1+A)(1+A+g{)}^2(1+J)J}-\frac{{\mu }_1{\mu }_2\sigma A^2(1+g)}{(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2}{(1+J)J}\left(\frac{(d_0\rho +(d_1-d_0\left)A\right)A}{(\rho -A{)}^2}+\frac{(\rho d_0+(d_1-d_0\left)A\right)A}{(\rho -A{)}^2}+\frac{(d_1-d_0)A(\rho -A)}{(\rho -A{)}^2}\right)J\\ & +\left(\frac{g}{1+A}-1-g\right)\frac{{\mu }_1{\mu }_2\sigma A^2}{(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2}{(1+J)J}\left(\frac{d_1\rho A}{(\rho -A{)}^2}+\frac{(d_0\rho +(d_1-d_0\left)A\right)A}{(\rho -A{)}^2}\right)J\\ & +\left(\frac{g}{1+A}-1-g\right)\frac{{\mu }_1{\mu }_2\sigma A^2}{(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2}{(1+J)J}\left(\frac{\left(\right(d_1+d_0)\rho +(d_1-d_0\left)A\right)A^2}{(\rho -A{)}^3}\right)+\left(\frac{g}{1+A}-1-g\right)\frac{{\mu }_1{\mu }_2\sigma A^2}{(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2{A}^2}{(1+J)J}\left(\frac{(d_1+d_0)\rho +(d_1-d_0)A}{(\rho -A{)}^3}+\sigma \left(\frac{bA}{{\mu }_2(1+A)}-\frac{{\mu }_2(1+A)+bA}{{\mu }_2}\right)\right.\\ & \left.\cdot \frac{{\mu }_2^2(1+A)}{\left({\mu }_2\right(1+A{)}^2+bA{)}^2}\right)\\ & =\frac{{\mu }_1{\mu }_2{A}^2}{(1+J)J}\left(\frac{(d_1+d_0)\rho +(d_1-d_0)A}{(\rho -A{)}^3}-\frac{{\mu }_2^2(1+A{)}^2+b{\mu }_2{A}^2}{\left({\mu }_2\right(1+A{)}^2+bA{)}^2}\right)\\ & =-\frac{{\mu }_1{\mu }_2{A}^2}{(1+J)J}{\mathrm{\Phi }}_H^{\prime }\left(A\right).\end{array}$(40)

Appendix 4. Proof of Q₁Q₂ - Q₃ > 0

It follows from Equation (37(37) $\begin{array}{rl}{Q}_1& ={\mu }_1+{\mu }_2+c_{11},\\ Q_2& ={\mu }_1{\mu }_2+c_{11}({\mu }_1+{\mu }_2)-c_{21}{c}_{12},\\ Q_3& ={\mu }_1{\mu }_2{c}_{11}+c_{32}{c}_{13}{c}_{21}-{\mu }_2{c}_{12}{c}_{21}.\end{array}$(37) ) that $\begin{array}{rl}{Q}_1{Q}_2-Q_3& =({\mu }_1+{\mu }_2+c_{11})({\mu }_1{\mu }_2+c_{11}({\mu }_1+{\mu }_2)-c_{21}{c}_{12})-{\mu }_1{\mu }_2{c}_{11}-c_{32}{c}_{13}{c}_{21}+{\mu }_2{c}_{12}{c}_{21}\\ & =({\mu }_1+{\mu }_2)({\mu }_1{\mu }_2+c_{11}({\mu }_1+{\mu }_2\left)\right)+c_{11}^2({\mu }_1+{\mu }_2)-({\mu }_1+c_{11})c_{21}{c}_{12}-c_{32}{c}_{13}{c}_{21}\\ & ={\mu }_1{\mu }_2({\mu }_1+{\mu }_2)+({\mu }_1+{\mu }_2{)}^2{c}_{11}+c_{11}(({\mu }_1+{\mu }_2)c_{11}-c_{12}{c}_{21})-{\mu }_1{c}_{12}{c}_{21}\\ & -c_{32}{c}_{13}{c}_{21}.\end{array}$ For convenience, we set ${\hat{c}}_{11}={\hat{a}}_{11}$ and ${\overline{c}}_{11}={\overline{a}}_{11}$. Thus $c_{11}={\hat{c}}_{11}+{\overline{c}}_{11}$ and we have $\begin{array}{rl}{Q}_1{Q}_2-Q_3& ={\mu }_1{\mu }_2({\mu }_1+{\mu }_2)+({\mu }_1+{\mu }_2{)}^2{\overline{c}}_{11}+{\mu }_1(({\mu }_1+{\mu }_2){\hat{c}}_{11}-c_{12}{c}_{21})\\ & +{\mu }_2({\mu }_1+{\mu }_2){\hat{c}}_{11}-c_{32}{c}_{13}{c}_{21}+c_{11}\left(\right({\mu }_1+{\mu }_2)c_{11}-c_{12}{c}_{21}).\end{array}$

If we assume ${\mu }_2\ge {\mu }_1$, then $\begin{array}{rl}& {\mu }_1\left(\right({\mu }_1+{\mu }_2){\hat{c}}_{11}-c_{12}{c}_{21})+{\mu }_2({\mu }_1+{\mu }_2){\hat{c}}_{11}-c_{32}{c}_{13}{c}_{21}\\ & =\frac{{\mu }_1({\mu }_1+{\mu }_2)\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1^2\sigma A^2}{(1+A+g)(1+J)J}\left(1+\frac{1+g}{1+A+g}\right)\\ & +\frac{{\mu }_2({\mu }_1+{\mu }_2)\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1{\mu }_2\sigma A^{2}g}{(1+A)(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1^2\sigma A^2(1+g)}{(1+A+g{)}^2(1+J)J}\\ & +\frac{{\mu }_2({\mu }_1+{\mu }_2)\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1{\mu }_2\sigma A^{2}g}{(1+A)(1+A+g{)}^2(1+J)J}\\ & =\frac{{\mu }_1{\mu }_2\sigma A^2}{(1+A+g)(1+J)J}\left(1-\frac{g}{(1+A)(1+A+g)}\right)\\ & +\frac{\sigma A^2}{(1+A+g)(1+J)J}\left({\mu }_2^2+{\mu }_1\left({\mu }_2-{\mu }_1\frac{1+g}{1+A+g}\right)\right)>0.\end{array}$ Same as in Equation (A2(39) $\begin{array}{rl}({\mu }_1+{\mu }_2)a_{11}-a_{12}{a}_{21}& =({\mu }_1+{\mu }_2){\hat{a}}_{11}-a_{12}{a}_{21}+({\mu }_1+{\mu }_2){\overline{a}}_{11}\\ & =\frac{({\mu }_1+{\mu }_2)\sigma A^2}{(1+A+g)(1+J)J}-\frac{{\mu }_1\sigma A^2}{(1+A+g)(1+J)J}\\ & \cdot \left(1+\frac{1+g}{1+A+g}\right)+({\mu }_1+{\mu }_2){\overline{a}}_{11}\\ & =\frac{\sigma A^2}{(1+A+g)(1+J)J}\left({\mu }_2-{\mu }_1\frac{1+g}{1+A+g}\right)+({\mu }_1+{\mu }_2){\overline{a}}_{11}>0.\end{array}$(39) ), we have, if ${\mu }_2\ge {\mu }_1$, $({\mu }_1+{\mu }_2)c_{11}-c_{12}{c}_{21}=({\mu }_1+{\mu }_2){\hat{c}}_{11}-c_{12}{c}_{21}+({\mu }_1+{\mu }_2){\overline{c}}_{11}>0.$ Therefore, $Q_1{Q}_2-Q_3>0$, if ${\mu }_2\ge {\mu }_1$.