Mathematical modelling and release thresholds of transgenic sterile mosquitoes

ABSTRACT

We formulate simple differential equation models to study the impact of releases of transgenic sterile mosquitoes carrying a dominant lethal on mosquito control based on the modified sterile insects technique. The early acting bisex, late acting bisex, early acting female-killing, and late acting female-killing lethality strategies are all considered. We determine release thresholds of the transgenic sterile mosquitoes, respectively, for these models by investigating the existence of positive equilibria and their stability. We compare the model dynamics, in particular, the thresholds of the models numerically. The late acting lethality strategies are generally more effective than their corresponding early acting lethality strategies, but the comparison between the late acting bisex and early acting female-killing lethality strategies depends on different parameter settings.

KEYWORDS:

• Dominant lethal gene
• positive equilibrium
• sterile insect technique
• release threshold
• transgenic mosquitoes

AMS SUBJECT CLASSIFICATIONS:

• 34D20
• 92D15
• 92D25
• 92D40
• 92D45

1. Introduction

Mosquito-borne diseases are a big concern for the public health and could threaten half the globe by 2050 [1]. The development of vaccines for malaria and dengue has made progress, but is still basically only in the trial stages. An effective way to prevent these diseases is to control mosquitoes. Massive spraying of insecticides or eliminating breeding sites has greatly reduced malaria transmissions in some areas, but the effects are limited. Killing mosquitoes with insecticides can result in an ecological void that could quickly be repopulated with surviving or migrant mosquitoes once control measures are suspended.

Biological controls of mosquitoes have been proven to be biologically and economically more effective. Among those control measures, the sterile insect technique (SIT) has been investigated in detail to reduce or eradicate wild mosquitoes. With SIT, the natural reproductive process of mosquitoes is disrupted. Utilizing radical or chemical methods, male mosquitoes are genetically modified to be sterile and then are released to compete for mates with wild male mosquitoes. Wild female mosquitoes mating with released sterile mosquitoes produce no progeny and then the wild mosquito population tends to decline [2–4].

Conventional SIT programmes have been well operated against major agricultural pests in some areas which had positive control effects and were promising [5]. However, there are currently no large-scale SIT programmes in operation against any mosquito species [6,7]. As further research is conducted, some concerns come up.

Irradiations and chemical methods are used in sterilizing mosquitoes. However, the irradiation of mosquito pupae can damage those pupae, irradiation of adults is operationally more difficult, and using chemicals can lead to trace contamination [7–9]. Moreover, mosquito populations are strongly density-dependent with intraspecific competition, particularly in aquatic stages, and thus are maintained at a relative stable level because of resource limitation. Conventional SIT targeting only the reproductive potential of adult mosquitoes therefore may not have significant impact on mosquito control [7,10].

To overcome these difficulties and improve the utilization of SIT, a modified genetic SIT is used and strategies of releasing mosquitoes carrying a dominant lethal (RIDL) are proposed and developed [5,11–13]. Field trials have also been conducted and have shown that hybrid offspring between the released transgenic and wild mosquitoes are sufficiently robust to be able to reproduce in nature [11].

In the modified genetic SIT, a dominant, repressible, lethal gene is used that induces repressible sex specific lethality. Mosquitoes are not irradiated but genetically sterilized, carrying one or more dominant lethal transgenes. There are potential advantages for genetic control over radiation-based Conventional SIT. Besides the easy sex separation, the ability to engineer the timing or nature of the lethality and thus vary the preferred time of death, or death phase, results in early or late acting lethality. Depending on different acting lethality, the modified genetic SIT can then be differently implemented as follows [10].

If the early acting bisex (EBS) lethality is implemented, the transgenic males carry both dominant female-specific and male-specific lethal genes. Wild females mating with the released transgenic males produce no viable offspring, which is similar to the conventional SIT. If the late acting bisex (LBS) lethality is applied, the transgenic males carry both dominant female-specific and male-specific lethal genes, but the genes act late in development only preventing mosquito larvae from becoming adults. Wild females mating with these transgenic males produce both male and female offspring, but the offspring only survive through the aquatic stage. In such a case, even though the transgenic larvae eventually die prior to adulthood, they increase larval competition reducing wild larvae's chances of survival. In the early acting female-killing (EFK) lethality implementation, the transgenic males carry a dominant female-specific lethal gene. Wild females mating with them produce male but no female offspring. The transgene can be passed on through male progeny. For the late acting female-killing (LFK) lethality implementation, the transgenic males carry a dominant female-specific lethal gene but the gene acts late in development preventing the female larvae from becoming adults. Wild females mating with these transgenic males produce both male and females offspring, but only male offspring, which carry the transgene, survive to adulthood. Similarly to the late acting bisex lethality, although the transgenic female larvae eventually die, they impose larval competition to wild larvae and thus reduce the survivability of wild larvae [5,7,10,13,14].

In order to have better understanding of the impact of utilizing the modified genetic SIT on mosquitoes control, mathematical models have been formulated, based on delayed differential equations, with the density-dependent survivability in the exponential form in [5,7,10,14–19]. Due to the nature of mathematical intractability in those studies, analytic results are limited and the models are basically only numerically solved in those existing studies. In this paper, we formulate homogeneous mosquito population models based on ordinary differential equations. We do not explicitly include metamorphic stages of mosquitoes, but the density-dependent birth functions in the models are based on the population size of the offspring. We respectively incorporate the four transgenic SIT implementations in the models. For simplification, we call the transgenic sterile mosquitoes sterile mosquitoes. As in [20–24], for the released sterile mosquitoes, we only consider and include those who are sexually active so that no separate dynamic equation is needed for the released sterile mosquitoes. We assume that sterile mosquitoes are constantly released with the release rate equal to constant c>0. We investigate the model dynamics and define release thresholds of the sterile mosquitoes for all the implementations to determine whether wild mosquitoes can be eradicated or coexist with the sterile mosquitoes. We discuss our findings and compare the different release strategies with numerical examples.

2. Modelling and release thresholds

We first consider a mosquito population model in the absence of sterile mosquitoes in this section. Let S be the population size of wild mosquitoes or the number of wild adults, and a be the per capitol birth rate such that aS is the total number of offspring or larvae that all wild mosquitoes produce. Because the intraspecific competition essentially occurs in the aquatic stages, we assume that the birth and thus the stage progression process is density-dependent based on the total number of offspring aS, with the logistic type of survivability [25]. The birth function thus takes the form $aS(1-\xi aS),$where $1/\xi$ is the carrying capacity of mosquitoes in aquatic stages, or the larvae. We let the death of wild adults be density-independent and denote the constant death rate by μ. The dynamics of the mosquito population is then governed by the following equation (1) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=a(1-\xi aS)S-\mu S.\end{array}$(1) We assume that the birth rate exceeds the death rate such that $a>\mu$. Equation (1(1) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=a(1-\xi aS)S-\mu S.\end{array}$(1) ) is a typical logistic-type equation. Its zero equilibrium is unstable and the mosquito population has a unique globally asymptotically stable positive equilibrium $S^0=\frac{a-\mu }{a\xi }.$

Remark 2.1

There are a good number of existing studies on interactive dynamics of wild and sterile mosquitoes in the literature where density dependence is assumed, such as [22–24,26–30], to list a few. However, either the density dependence is included in the adults death, or the density-dependent birth is based on the total population size, not the number of offspring.

2.1. Early acting bisex (EBS)

We now consider the early acting bisex lethality strategy. Wild female mosquitoes mating with released male sterile mosquitoes produce no offspring. Assume that the released sterile mosquitoes are maintained at a constant number c. The total number of mosquitoes is S + c and then the fraction of matings of a wild mosquito with wild mosquitoes is $S/(S+c)$. As in (1(1) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=a(1-\xi aS)S-\mu S.\end{array}$(1) ), the intraspecific competition occurs in the aquatic stages, and the total number of wild larvae is $a(S/(S+c))S$ in the density dependent birth term. The model equation then becomes (2) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS}{S+c}S)S-\mu S.\end{array}$(2) The zero equilibrium of (2(2) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS}{S+c}S)S-\mu S.\end{array}$(2) ) is always locally asymptotically stable. To find a positive equilibrium of (2(2) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS}{S+c}S)S-\mu S.\end{array}$(2) ), we need to solve the following equation for a positive root $\frac{aS}{S+c}(1-\xi \frac{aS}{S+c}S)-\mu =0,$that is, equivalently to solve the cubic equation (3) $\begin{array}{r}{F}_1(S):=a^2\xi S^3-(a-\mu )S^2-(a-2\mu )cS+\mu c^2=0\end{array}$(3) for S>0. It is clear that Equation (3(3) $\begin{array}{r}{F}_1(S):=a^2\xi S^3-(a-\mu )S^2-(a-2\mu )cS+\mu c^2=0\end{array}$(3) ) can have at most two positive solutions.

It follows from $F_1^{\prime }(S)=3a^2\xi S^2-2(a-\mu )S-(a-2\mu )c=0$that if $a\ge 2\mu$, $F_1(S)$ has a unique critical point (4) $\begin{array}{r}{\overline{S}}_1=\frac{a-\mu +\sqrt{(a-\mu {)}^2+3a^2\xi (a-2\mu )c}}{3a^2\xi }.\end{array}$(4) At this unique critical point ${\overline{S}}_1$, $\begin{array}{rl}{F}_1({\overline{S}}_1)& =\frac{1}{3}{\overline{S}}_1{F}_1^{\prime }({\overline{S}}_1)-\frac{1}{3}(a-\mu ){\overline{S}}_1^2-\frac{2}{3}(a-2\mu )c{\overline{S}}_1+\mu c^2\\ & =-\frac{a-\mu }{9a^2\xi }{F}_1^{\prime }({\overline{S}}_1)-(\frac{2(a-\mu {)}^2}{9a^2\xi }+\frac{2}{3}(a-2\mu )c){\overline{S}}_1+\mu c^2-\frac{(a-\mu )(a-2\mu )c}{9a^2\xi }\\ & =-2\frac{(a-\mu {)}^2+3a^2\xi (a-2\mu )c}{9a^2\xi }{\overline{S}}_1+\frac{9a^2\xi \mu c-(a-\mu )(a-2\mu )}{9a^2\xi }c.\end{array}$Then Equation (3(3) $\begin{array}{r}{F}_1(S):=a^2\xi S^3-(a-\mu )S^2-(a-2\mu )cS+\mu c^2=0\end{array}$(3) ) has no, one, or two positive equilibria if $F_1({\overline{S}}_1)$ is greater than, equal to, or less than zero.

We now determine the release threshold for the sterile mosquitoes. Setting $F_1({\overline{S}}_1)=0$ and substituting (4(4) $\begin{array}{r}{\overline{S}}_1=\frac{a-\mu +\sqrt{(a-\mu {)}^2+3a^2\xi (a-2\mu )c}}{3a^2\xi }.\end{array}$(4) ) into it, we have ${\overline{S}}_1=\frac{c(9a^2\xi \mu c-(a-\mu )(a-2\mu ))}{2((a-\mu {)}^2+3a^2\xi (a-2\mu )c)}=\frac{a-\mu +\sqrt{(a-\mu {)}^2+3a^2\xi (a-2\mu )c}}{3a^2\xi },$that is, (5) $\begin{array}{rl}& 3a^2\xi c(9a^2\xi \mu c-(a-\mu )(a-2\mu ))-2(a-\mu )((a-\mu {)}^2+3a^2\xi (a-2\mu )c)\\ & =2((a-\mu {)}^2+3a^2\xi (a-2\mu )c)\sqrt{(a-\mu {)}^2+3a^2\xi (a-2\mu )c}.\end{array}$(5) Subtracting the squaring of the right hand side from the square of the left hand side of (5(5) $\begin{array}{rl}& 3a^2\xi c(9a^2\xi \mu c-(a-\mu )(a-2\mu ))-2(a-\mu )((a-\mu {)}^2+3a^2\xi (a-2\mu )c)\\ & =2((a-\mu {)}^2+3a^2\xi (a-2\mu )c)\sqrt{(a-\mu {)}^2+3a^2\xi (a-2\mu )c}.\end{array}$(5) ), we then define a function of c $\begin{array}{rl}{H}_1(c)& :=729a^8{\xi }^4{\mu }^2{c}^4-(486a^6{\xi }^3\mu (a-\mu )(a-2\mu )+108a^6{\xi }^3(a-2\mu {)}^3)c^3\\ & -9a^4{\xi }^2(a-\mu {)}^2(12\mu (a-\mu )+3(a-2\mu {)}^2)c^2\\ & =27a^6{\xi }^2(27a^2{\xi }^2{\mu }^2{c}^2-2\xi (a-2\mu )(2a-\mu )(a+\mu )c-(a-\mu {)}^2)c^2.\end{array}$It is clear that equation $H_1(c)=0$ has a unique positive root (6) $\begin{array}{r}\frac{(a-2\mu )(2a-\mu )(a+\mu )}{27a^2\xi {\mu }^2}+\frac{\sqrt{(a-2\mu {)}^2(2a-\mu {)}^2(a+\mu {)}^2+27a^2{\mu }^2(a-\mu {)}^2}}{27a^2\xi {\mu }^2}.\end{array}$(6) We denote it $c_{EBS}^0$, and define it as the threshold of release of sterile mosquitoes for model (2(2) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS}{S+c}S)S-\mu S.\end{array}$(2) ). Note that $\underset{c\to 0+}{lim}{F}_1(c,{\overline{S}}_1(c))=-\frac{4(a-\mu {)}^3}{27a^3{\xi }^2}<0.$Then $F_1({\overline{S}}_1)$ is less than, equal to, or greater than zero; that is, there exist two, one, or no positive equilibrium, if $0<c<c_{EBS}^0$, $c=c_{EBS}^0$, or $c>c_{EBS}^0$, respectively.

It is also easy to see that if there exists no positive equilibrium, the origin is globally asymptotically stable, and if there exists a unique positive equilibrium $S_{EBS}^{\ast }$, it follows from $\frac{\mathrm{d}S}{\mathrm{d}t}=-\frac{S{F}_1(S)}{(S+c{)}^2}=-\frac{a^2\xi S(S+k_1){(S-S_{EBS}^{\ast })}^2}{(S+c{)}^2},$where $k_1>0$ is a constant, that solutions with initial values less than $S_{EBS}^{\ast }$ approach zero, and solutions with initial values greater than $S_{EBS}^{\ast }$ approach $S_{EBS}^{\ast }$. Thus, $S_{EBS}^{\ast }$ is unstable or semi-stable.

Suppose that there exist two positive equilibria $S_{EBS}^{-}<S_{EBS}^{+}$. Then it follows from $\frac{\mathrm{d}S}{\mathrm{d}t}=-\frac{a^2\xi S(S+k_2)(S-S_{EBS}^{-})(S-S_{EBS}^{+})}{(S+c{)}^2},$where $k_2>0$ is a constant, that solutions approach zero if their initial values are less than $S_{EBS}^{-}$, approach $S_{EBS}^{+}$ if their initial values are between $S_{EBS}^{-}$ and $S_{EBS}^{+}$, or greater than $S_{EBS}^{+}$. Thus, equilibrium $S_{EBS}^{-}$ is unstable and $S_{EBS}^{+}$ is locally asymptotically stable. In summary, we have the following results.

Theorem 2.1

For model (2(2) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS}{S+c}S)S-\mu S.\end{array}$(2) ), the zero solution is always locally asymptotically stable. We assume $a\ge 2\mu$ and define the release threshold of sterile mosquitoes $c_{EBS}^0$ in (6(6) $\begin{array}{r}\frac{(a-2\mu )(2a-\mu )(a+\mu )}{27a^2\xi {\mu }^2}+\frac{\sqrt{(a-2\mu {)}^2(2a-\mu {)}^2(a+\mu {)}^2+27a^2{\mu }^2(a-\mu {)}^2}}{27a^2\xi {\mu }^2}.\end{array}$(6) ). If $0<c<c_{EBS}^0$, there exist two positive equilibria $S_{EBS}^{-}<S_{EBS}^{+}$. Equilibrium $S_{EBS}^{-}$ is unstable and $S_{EBS}^{+}$ is locally asymptotically stable. If $c=c_{EBS}^0$, there exists a unique positive equilibrium $S_{EBS}^{\ast }$, which is unstable or semi-stable. If $c>c_{EBS}^0$, there exists no positive equilibrium. The zero solution is globally asymptotically stable in this case.

Remark 2.2

Notice that while the matings between the wild and sterile mosquitoes lead to producing reduced number of wild offspring, this reduction also results in less crowding and thus weakening the intraspecific competition among the offspring. It therefore increases the survivability of wild larvae. Comparing with many existing works in the literature such as [25] where the lessening crowdedness is not considered, we can see that while the release of sterile mosquitoes plays the same role in reducing the birth of wild mosquitoes, the inclusion of the less crowdedness and thus larger survivability of wild offspring can lead to producing more wild offspring. This can be also seen from $\frac{a{S}^2}{S+c}(1-\xi \frac{aS}{S+c}S)-\frac{a{S}^2}{S+c}(1-\xi aS)=\frac{\xi a^2{S}^3}{S+c}(1-\frac{S}{S+c})>0,$where the first term is the birth function of the models including the reduction of wild offspring due the matings of wild with sterile mosquitoes and the second term is the birth function of the models missing such a reduction.

2.2. Late acting bisex (LBS)

In this section, we assume that the late acting bisex lethality strategy is implemented. Then wild mosquitoes mating with released sterile mosquitoes produce both male and female offspring, but all the offspring only survive through the aquatic stage. Thus, they have no effects on the reproduction but impose larval intraspecific competition and then affect the larvae survivability.

We assume equal sex-ratio for male and female mosquitoes such that the birth rates of male and female offspring are the same. Let b be the number of male or female offspring produced by a wild mosquito mating with a sterile mosquito per unit of time. Since the fraction of sterile mosquitoes in the total mosquito population is $c/(S+c)$, the total number of offspring produced by all wild mosquitoes mating with sterile mosquitoes is $\frac{2bc}{S+c}S.$In the mean time, the total number of wild offspring produced by all wild mosquitoes mating with wild mosquitoes is $a(S/(S+c))S$. Then the total number of all offspring, involving in the larval intraspecific competition and affecting the density dependent survivability, is $\frac{aS+2bc}{S+c}S.$The model dynamics for the case of LBS thus are governed by the following equation (7) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S.\end{array}$(7)

2.2.1. Equal reproduction

We first assume that the matings between wild and sterile mosquitoes produce the same number of offspring such that a = 2b. Then the total number of all offspring is aS and the release of sterile mosquitoes has no effects on the density dependent survivability; that is, it does not reduce the intraspecific competition and therefore does not help wild larvae survival as in the case of EBS. The model equation becomes (8) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi aS)S-\mu S.\end{array}$(8) The origin is locally asymptotically stable. Define the release threshold of sterile mosquitoes by (9) $\begin{array}{r}{c}_{LBS}^0:=\frac{(a-\mu {)}^2}{4a^2\xi \mu }.\end{array}$(9) Then Equation (8(8) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi aS)S-\mu S.\end{array}$(8) ) has no positive equilibrium, a unique positive equilibrium (10) $\begin{array}{r}{S}_{LBS}^{\ast }=\frac{a-\mu }{2a^2\xi },\end{array}$(10) or two positive equilibria (11) $\begin{array}{r}{S}_{LBS}^{\mp }=\frac{a-\mu \mp \sqrt{(a-\mu {)}^2-4a^2\xi \mu c}}{2a^2\xi },\end{array}$(11) if $c>c_{LBS}^0$, $c=c_{LBS}^1$, or $c<c_{LBS}^0$, respectively. The stability analysis is similar and is thus omitted. In summary, we have the following results.

Theorem 2.2

Suppose that wild mosquitoes produce the same number of offspring from matings with wild or sterile mosquitoes such that a = 2b. Define the release threshold $c_{LBS}^0$ in (9(9) $\begin{array}{r}{c}_{LBS}^0:=\frac{(a-\mu {)}^2}{4a^2\xi \mu }.\end{array}$(9) ) for model (8(8) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi aS)S-\mu S.\end{array}$(8) ). If $c>c_{LBS}^0$, there is no positive equilibrium and the zero equilibrium is globally asymptotically stable. If $c=c_{LBS}^0$, the unique positive equilibrium $S_{LBS}^{\ast }$, given in (10(10) $\begin{array}{r}{S}_{LBS}^{\ast }=\frac{a-\mu }{2a^2\xi },\end{array}$(10) ), is semi-stable, and if $c<c_{LBS}^0$, positive equilibria $S_{LBS}^{-}$ and $S_{LBS}^{+}$, given in (11(11) $\begin{array}{r}{S}_{LBS}^{\mp }=\frac{a-\mu \mp \sqrt{(a-\mu {)}^2-4a^2\xi \mu c}}{2a^2\xi },\end{array}$(11) ), are unstable and locally asymptotically stable, respectively.

2.2.2. Nonequal reproduction

We next assume that wild mosquitoes produce less offspring from matings with sterile mosquitoes than with wild mosquitoes due to possible fitness costs such that 2b<a. Then, similarly, to find a positive equilibrium, we define function (12) $\begin{array}{r}{F}_2(S):=a^2\xi S^3-(a-\mu -2ab\xi c)S^2-(a-2\mu )cS+\mu c^2.\end{array}$(12) It follows from $F_2^{\prime }(S)=3a^2\xi S^2-2(a-\mu -2ab\xi c)S-(a-2\mu )c$that $F_2(S)$ has a unique critical point (13) $\begin{array}{r}{\overline{S}}_2:=\frac{a-\mu -2ab\xi c+\sqrt{(a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c}}{3a^2\xi }.\end{array}$(13) At this critical point ${\overline{S}}_2$, $\begin{array}{rl}{F}_2({\overline{S}}_2)& =-\frac{a-\mu -2ab\xi c}{9a^2\xi }(3a^2\xi {\overline{S}}_2^2-2(a-\mu -2ab\xi c){\overline{S}}_2-(a-2\mu )c)\\ & -(\frac{2}{3}(a-2\mu )c+\frac{2(a-\mu -2ab\xi c{)}^2}{9a^2\xi }){\overline{S}}_2\\ & +\mu c^2-\frac{(a-\mu -2ab\xi c)(a-2\mu )c}{9a^2\xi }\\ & =-\frac{2(3a^2\xi (a-2\mu )c+(a-\mu -2ab\xi c{)}^2)}{9a^2\xi }{\overline{S}}_2\\ & +\frac{9a^2\xi c-(a-\mu -2ab\xi )(a-2\mu )c}{9a^2\xi }c.\end{array}$Solving $F_2({\overline{S}}_2)=0$ for ${\overline{S}}_2$ and then following from (13(13) $\begin{array}{r}{\overline{S}}_2:=\frac{a-\mu -2ab\xi c+\sqrt{(a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c}}{3a^2\xi }.\end{array}$(13) ), we have $\begin{array}{rl}{\overline{S}}_2& =\frac{c(9a^2\xi \mu c-(a-\mu -2ab\xi c)(a-2\mu ))}{2((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)}\\ & =\frac{a-\mu -2ab\xi c+\sqrt{(a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c}}{3a^2\xi },\end{array}$that is, (14) $\begin{array}{rl}& 3a^2\xi c(9a^2\xi \mu c-(a-\mu -2ab\xi c)(a-2\mu ))\\ & -2(a-\mu -2ab\xi c)((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\\ & =2((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\sqrt{(a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c}.\end{array}$(14) By subtracting the square of the right-hand side of (14(14) $\begin{array}{rl}& 3a^2\xi c(9a^2\xi \mu c-(a-\mu -2ab\xi c)(a-2\mu ))\\ & -2(a-\mu -2ab\xi c)((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\\ & =2((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\sqrt{(a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c}.\end{array}$(14) ) from the square of the left-hand side of (14(14) $\begin{array}{rl}& 3a^2\xi c(9a^2\xi \mu c-(a-\mu -2ab\xi c)(a-2\mu ))\\ & -2(a-\mu -2ab\xi c)((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\\ & =2((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\sqrt{(a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c}.\end{array}$(14) ), and then dividing it by $c^2$, equation (14(14) $\begin{array}{rl}& 3a^2\xi c(9a^2\xi \mu c-(a-\mu -2ab\xi c)(a-2\mu ))\\ & -2(a-\mu -2ab\xi c)((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\\ & =2((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\sqrt{(a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c}.\end{array}$(14) ) is equivalent to (15) $\begin{array}{rl}& 729a^8{\xi }^4{\mu }^2{c}^2-(486a^6{\xi }^3\mu (a-\mu -2ab\xi c)(a-2\mu )+108a^6{\xi }^3(a-2\mu {)}^3)c\\ & -27a^4{\xi }^2(a-\mu -2ab\xi c{)}^2(4\mu (a-\mu -2ab\xi c)+(a-2\mu {)}^2)=0.\end{array}$(15) Further dividing (15(15) $\begin{array}{rl}& 729a^8{\xi }^4{\mu }^2{c}^2-(486a^6{\xi }^3\mu (a-\mu -2ab\xi c)(a-2\mu )+108a^6{\xi }^3(a-2\mu {)}^3)c\\ & -27a^4{\xi }^2(a-\mu -2ab\xi c{)}^2(4\mu (a-\mu -2ab\xi c)+(a-2\mu {)}^2)=0.\end{array}$(15) ) by $27a^5{\xi }^2$, we obtain the following function of c (16) $\begin{array}{rl}{H}_2(c)& :=27a^3{\xi }^2{\mu }^2{c}^2-2a\xi (a-2\mu )(2a-\mu )(a+\mu )c+36a^{2}b{\xi }^2\mu (a-2\mu )c^2\\ & +32a^2{b}^3{\xi }^3\mu c^3-4a{b}^2{\xi }^2(a^2+8\mu (a-\mu ))c^2\\ & +4b\xi (a-\mu )(a^2+2\mu (a-\mu ))c-a(a-\mu {)}^2\\ & =p_3{c}^3+p_2{c}^2+p_{1}c+p_0,\end{array}$(16) where $\begin{array}{rl}{p}_3& :=32a^2{b}^3{\xi }^3\mu ,\\ p_2& :=a{\xi }^2(27a^2{\mu }^2+36ab\mu (a-2\mu )-4b^2(a^2+8\mu (a-\mu ))),\\ p_1& :=2\xi (2b(a-\mu )(a^2+2\mu (a-\mu ))-a(a-2\mu )(2a-\mu )(a+\mu ))\\ & =2\xi b((a^3+a^2\mu -4a{\mu }^2+2{\mu }^3)-a(2a^3-3a^2\mu -3a{\mu }^2+2{\mu }^3)),\\ p_0& :=-a(a-\mu {)}^2.\end{array}$A positive number c is a solution of (15(15) $\begin{array}{rl}& 729a^8{\xi }^4{\mu }^2{c}^2-(486a^6{\xi }^3\mu (a-\mu -2ab\xi c)(a-2\mu )+108a^6{\xi }^3(a-2\mu {)}^3)c\\ & -27a^4{\xi }^2(a-\mu -2ab\xi c{)}^2(4\mu (a-\mu -2ab\xi c)+(a-2\mu {)}^2)=0.\end{array}$(15) ), or (14(14) $\begin{array}{rl}& 3a^2\xi c(9a^2\xi \mu c-(a-\mu -2ab\xi c)(a-2\mu ))\\ & -2(a-\mu -2ab\xi c)((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\\ & =2((a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c)\sqrt{(a-\mu -2ab\xi c{)}^2+3a^2\xi (a-2\mu )c}.\end{array}$(14) ), if and only if it satisfies $H_2(c)=0$.

Since $p_0<0$, it is clear that if $p_1\le 0$, $H_2(c)=0$ has a unique positive solution. Now we assume $p_1>0$, that is, (17) $\begin{array}{r}a(2a^3-3a^2\mu -3a{\mu }^2+2{\mu }^3)<2b(a^3+a^2\mu -4a{\mu }^2+2{\mu }^3).\end{array}$(17) Since 2b<a, inequality (17(17) $\begin{array}{r}a(2a^3-3a^2\mu -3a{\mu }^2+2{\mu }^3)<2b(a^3+a^2\mu -4a{\mu }^2+2{\mu }^3).\end{array}$(17) ) implies $2a^3-3a^2\mu -3a{\mu }^2+2{\mu }^3<a^3+a^2\mu -4a{\mu }^2+2{\mu }^3,$that is (18) $\begin{array}{r}{a}^2-4a\mu +{\mu }^2<0.\end{array}$(18) Thus $\begin{array}{rl}\frac{1}{a{\xi }^2}{p}_2& =27a^2{\mu }^2+36ab\mu (a-2\mu )-4b^2(a^2+8\mu (a-\mu ))\\ & \ge 4b^2(27{\mu }^2+18\mu (a-2\mu )-(a^2+8\mu (a-\mu ))),\end{array}$and from (18(18) $\begin{array}{r}{a}^2-4a\mu +{\mu }^2<0.\end{array}$(18) ), $\begin{array}{rl}27& {\mu }^2+18\mu (a-2\mu )-(a^2+8\mu (a-\mu ))\\ & =27{\mu }^2+18a\mu -36{\mu }^2-(a^2-4a\mu +{\mu }^2)-12a\mu +9{\mu }^2\\ & =6a\mu -(a^2-4a\mu +{\mu }^2)>0.\end{array}$Hence, if $p_1>0$, coefficient $p_2$ must be positive as well. From Descartes' Rule of Signs, $H_2(c)=0$ has a unique positive solution. We define this unique positive solution as the release threshold of sterile mosquitoes for model (7(7) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S.\end{array}$(7) ) and denote it $c_{LBS}^0$. Similarly to the assertion in Section 2.1, model (7(7) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S.\end{array}$(7) ) has two, one, or no positive equilibrium if $0<c<c_{LBS}^0$, $c=c_{LBS}^0$, or $c>c_{LBS}^0$, respectively.

The stability of the positive equilibrium or equilibria of (7(7) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S.\end{array}$(7) ) can be shown in a way similar to that for (2(2) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS}{S+c}S)S-\mu S.\end{array}$(2) ). We omit the details and summarize the results as follows.

Theorem 2.3

Assume that wild mosquitoes produce less offspring from matings with sterile mosquitoes than with wild mosquitoes due to possible fitness cost such that 2b<a. The zero solution of model (7(7) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S.\end{array}$(7) ) is always locally asymptotically stable. Define the release threshold of sterile mosquitoes $c_{LBS}^0$ by the unique positive solution of $H_2(c)=0$ in (16(16) $\begin{array}{rl}{H}_2(c)& :=27a^3{\xi }^2{\mu }^2{c}^2-2a\xi (a-2\mu )(2a-\mu )(a+\mu )c+36a^{2}b{\xi }^2\mu (a-2\mu )c^2\\ & +32a^2{b}^3{\xi }^3\mu c^3-4a{b}^2{\xi }^2(a^2+8\mu (a-\mu ))c^2\\ & +4b\xi (a-\mu )(a^2+2\mu (a-\mu ))c-a(a-\mu {)}^2\\ & =p_3{c}^3+p_2{c}^2+p_{1}c+p_0,\end{array}$(16) ). Then if $0<c<c_{LBS}^0$, there exist two positive equilibria $S_{LBS}^{-}<S_{LBS}^{+}$. Equilibrium $S_{LBS}^{-}$ is unstable and $S_{LBS}^{+}$ is locally asymptotically stable. If $c=c_{LBS}^0$, there exists a unique positive equilibrium $S_{LBS}^{\ast }$, which is unstable or semi-stable. If $c>c_{LBS}^0$, there exists no positive equilibrium and the zero solution is globally asymptotically stable.

2.3. Early acting female-killing (EFK)

We assume the early acting female-killing lethality strategy in this section. Wild mosquitoes mating with sterile mosquitoes produce male but no female offspring and the transgene can be further passed on to the male progeny. Let a constant proportion γ of sterile male offspring survive and emerge to adults. Then the number of sterile male offspring produced by all wild mosquitoes mating with sterile mosquitoes is $\frac{bc}{S+c}S$, and the number of sterile male adults emerging from the sterile larvae is $\frac{\gamma bc}{S+c}S$. Again, the fraction of sterile mosquitoes in the total mosquito population before the transgenic larvae emerging to adulthood is still $\frac{c}{S+c}$, and the total number of wild mosquitoes produced by all wild mosquitoes mating with wild mosquitoes is $\frac{aS}{S+c}S$. Then the total number of offspring affecting the density dependence is $\frac{aS+bc}{S+c}S,$and the total number of sterile male adults, including those emerging from the sterile larvae, becomes $c+\frac{\gamma bc}{S+c}S.$Here, we simplify the model by ignoring the progression time delay from offspring to adult. Therefore, the model for the early acting female-killing lethality strategy has the following form (19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) The zero equilibrium of (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ) is locally asymptotically stable. We then investigate the existence of positive equilibria which satisfy (20) $\begin{array}{rl}{F}_3(S)& :=\mu (S^2+(2+\gamma b)cS+c^2)+aS(a\xi S^2-(1-b\xi c)S-c))\\ & =a^2\xi S^3-(a-\mu -ab\xi c)S^2-(a-(2+b\gamma )\mu )cS+\mu c^2=0.\end{array}$(20) Equation (20(20) $\begin{array}{rl}{F}_3(S)& :=\mu (S^2+(2+\gamma b)cS+c^2)+aS(a\xi S^2-(1-b\xi c)S-c))\\ & =a^2\xi S^3-(a-\mu -ab\xi c)S^2-(a-(2+b\gamma )\mu )cS+\mu c^2=0.\end{array}$(20) ) has no positive solutions if $a\le (2+\gamma b)\mu \text{and}c\ge \frac{a-\mu }{ab\xi }.$We assume (21) $\begin{array}{r}a>(2+\gamma b)\mu .\end{array}$(21) Equation $F_3^{\prime }(S)=3a^2\xi S^2-2(a-\mu -ab\xi c)S-(a-2\mu -b\gamma \mu )c=0$has a unique positive root (22) $\begin{array}{r}{\overline{S}}_3=\frac{a-\mu -ab\xi c+\sqrt{{(a-\mu -ab\xi c)}^2+3a^2\xi c(a-2\mu -b\gamma \mu )}}{3a^2\xi }.\end{array}$(22) Since $F_3(0)=\mu c^2>0$, then if $F_3({\overline{S}}_3)>0$, $=0$, or $<0$, there exist no, one, or two positive equilibria of Equation (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ).

To determine the release threshold for (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ), it is similar to the manner in Section 2.2.2 that, at ${\overline{S}}_3$, $\begin{array}{rl}{F}_3({\overline{S}}_3)& =-\frac{2(3a^2\xi c(a-2\mu -b\gamma \mu )+(a-\mu -ab\xi c{)}^2)}{9a^2\xi }{\overline{S}}_3\\ & +\frac{9a^2\xi \mu c-(a-\mu -ab\xi c)(a-2\mu -b\gamma \mu )}{9a^2\xi }c.\end{array}$Solving $F_3({\overline{S}}_3)=0$ together from (22(22) $\begin{array}{r}{\overline{S}}_3=\frac{a-\mu -ab\xi c+\sqrt{{(a-\mu -ab\xi c)}^2+3a^2\xi c(a-2\mu -b\gamma \mu )}}{3a^2\xi }.\end{array}$(22) ), we have $\begin{array}{rl}{\overline{S}}_3& =\frac{c(9a^2\xi \mu c-(a-\mu -ab\xi c)(a-2\mu -b\gamma \mu ))}{2(3a^2\xi c(a-2\mu -b\gamma \mu )+(a-\mu -ab\xi c{)}^2)}\\ & =\frac{a-\mu -ab\xi c+\sqrt{{(a-\mu -ab\xi c)}^2+3a^2\xi c(a-2\mu -b\gamma \mu )}}{3a^2\xi },\end{array}$that is, (23) $\begin{array}{rl}& 3a^2\xi c(9a^2\xi \mu c-(a-\mu -ab\xi c)(a-2\mu -b\gamma \mu ))\\ & -2(3a^2\xi c(a-2\mu -b\gamma \mu )(a-\mu -ab\xi c)\\ & +{(a-(2+\gamma b)\mu )}^3)\\ & =2(3a^2\xi c(a-2\mu -b\gamma \mu )+{(a-\mu -ab\xi c)}^2)\\ & \times \sqrt{{(a-\mu -ab\xi c)}^2+3a^2\xi c(a-2\mu -b\gamma \mu )}.\end{array}$(23) Subtracting the squaring of the right hand side from the square of the left hand side of (23(23) $\begin{array}{rl}& 3a^2\xi c(9a^2\xi \mu c-(a-\mu -ab\xi c)(a-2\mu -b\gamma \mu ))\\ & -2(3a^2\xi c(a-2\mu -b\gamma \mu )(a-\mu -ab\xi c)\\ & +{(a-(2+\gamma b)\mu )}^3)\\ & =2(3a^2\xi c(a-2\mu -b\gamma \mu )+{(a-\mu -ab\xi c)}^2)\\ & \times \sqrt{{(a-\mu -ab\xi c)}^2+3a^2\xi c(a-2\mu -b\gamma \mu )}.\end{array}$(23) ) and then dividing it by $108a^4{\xi }^2{c}^2$, we define (24) $\begin{array}{rl}{H}_3(c)& :=27a^4{\xi }^2{\mu }^2{c}^2-18a^2\xi \mu c(a-2\mu -b\gamma \mu )(a-\mu -ab\xi c)\\ & -4a^2\xi c(a-2\mu -b\gamma \mu {)}^3-4\mu {(a-\mu -ab\xi c)}^3\\ & -(a-2\mu -b\gamma \mu {)}^2{(a-\mu -ab\xi c)}^2.\end{array}$(24) Rewrite (25) $\begin{array}{r}{H}_3(c):=k_3{c}^3+k_2{c}^2+k_{1}c+k_0,\end{array}$(25) where $\begin{array}{rl}{k}_3& =4a^3{b}^3{\xi }^3\mu ,\\ k_2& =a^2{\xi }^2(27a^2{\mu }^2+18ab\mu (a-2\mu -b\gamma \mu )-12b^2\mu (a-\mu )-b^2(a-2\mu -b\gamma \mu {)}^2)\\ & =a^2{\xi }^2(27a^2{\mu }^2-b^2(a-2\mu -b\gamma \mu {)}^2+6b\mu (3a(a-2\mu -b\gamma \mu )-2b(a-\mu ))),\\ k_1& =2a\xi (-9a\mu (a-\mu )(a-2\mu -b\gamma \mu )+6b\mu (a-\mu {)}^2-2a(a-2\mu -b\gamma \mu {)}^3\\ & +b(a-\mu )(a-2\mu -b\gamma \mu {)}^2)\\ & =2a\xi (-3\mu (a-\mu )(3a(a-2\mu -b\gamma \mu )-2b(a-\mu ))\\ & -(a-2\mu -b\gamma \mu {)}^2(2a(a-2\mu -b\gamma \mu )-b(a-\mu ))),\\ k_0& =-(a-\mu {)}^2(4\mu (a-\mu )+(a-2\mu -b\gamma \mu {)}^2).\end{array}$If (26) $\begin{array}{r}3a(a-2\mu -b\gamma \mu )\ge 2b(a-\mu ),\end{array}$(26) then $2a(a-2\mu -b\gamma \mu )-b(a-\mu )>0$ and thus $k_1<0$. Hence $H_3(c)=0$ has a unique positive root that we denote $c_{EFK}^0$.

Since $2b\le a$, then $\begin{array}{rl}2(3a(a-2\mu -b\gamma \mu )-2b(a-\mu ))& \ge 3a(2a-4\mu -a\gamma \mu )-2a(a-\mu )\\ & =a((4-3\gamma \mu )a-10\mu ).\end{array}$We assume (27) $\begin{array}{r}a\ge \frac{10\mu }{4-3\gamma \mu }.\end{array}$(27) Then condition (26(26) $\begin{array}{r}3a(a-2\mu -b\gamma \mu )\ge 2b(a-\mu ),\end{array}$(26) ) is satisfied. Moreover, if (28) $\begin{array}{r}a\ge \frac{4\mu }{2-\gamma \mu },\end{array}$(28) it follows from $(2+\gamma b)\mu \le 2\mu +\gamma \mu \frac{a}{2},$that condition (21(21) $\begin{array}{r}a>(2+\gamma b)\mu .\end{array}$(21) ) is satisfied. Since $\frac{10\mu }{4-3\gamma \mu }\ge \frac{4\mu }{2-\gamma \mu },$condition (27(27) $\begin{array}{r}a\ge \frac{10\mu }{4-3\gamma \mu }.\end{array}$(27) ) implies (26(26) $\begin{array}{r}3a(a-2\mu -b\gamma \mu )\ge 2b(a-\mu ),\end{array}$(26) ) and then (21(21) $\begin{array}{r}a>(2+\gamma b)\mu .\end{array}$(21) ). Therefore, if (27(27) $\begin{array}{r}a\ge \frac{10\mu }{4-3\gamma \mu }.\end{array}$(27) ) holds, $H_3(c)=0$ has a unique positive root.

Similarly, since $\underset{c\to 0+}{lim}{F}_3(c,\overline{S}(c))=-\frac{2(a-\mu {)}^3}{27a^4{\xi }^2}<0,$the positive root of $H_3(c)$ is the release threshold of sterile mosquitoes for model (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ).

The stability of the positive equilibrium and equilibria can be determined similarly as in Section 2.2. The details are then omitted. We summarize our results as follows.

Theorem 2.4

Let the birth rate of wild mosquitoes a, the death rate of wild mosquitoes μ, and the emergence rate of sterile larvae γ satisfy condition (27(27) $\begin{array}{r}a\ge \frac{10\mu }{4-3\gamma \mu }.\end{array}$(27) ). Then the release threshold of sterile mosquitoes $c_{EFK}^0$ is given by the positive solution of equation $H_3(c)=0$, where $H_3(c)$ is given in (24(24) $\begin{array}{rl}{H}_3(c)& :=27a^4{\xi }^2{\mu }^2{c}^2-18a^2\xi \mu c(a-2\mu -b\gamma \mu )(a-\mu -ab\xi c)\\ & -4a^2\xi c(a-2\mu -b\gamma \mu {)}^3-4\mu {(a-\mu -ab\xi c)}^3\\ & -(a-2\mu -b\gamma \mu {)}^2{(a-\mu -ab\xi c)}^2.\end{array}$(24) ) or (25(25) $\begin{array}{r}{H}_3(c):=k_3{c}^3+k_2{c}^2+k_{1}c+k_0,\end{array}$(25) ). If $c<c_{EFK}^0$, model (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ) has two positive equilibria, $S_{EFK}^{-}<S_{EFK}^{+}$. Equilibrium $S_{EFK}^{-}$ is unstable and $S_{EFK}^{+}$ is locally asymptotically stable. If $c=c_{EFK}^0$, there exists a unique positive equilibrium $S_{EFK}^{\ast }$, which is unstable or semi-stable. If $c>c_{EFK}^0$, there exists no positive equilibrium and the zero solution is globally asymptotically stable.

2.4. Late acting female-killing (LFK)

In this section, we assume the late acting female-killing lethality strategy. Wild mosquitoes mating with released sterile mosquitoes produce both male and female offspring, but the female offspring die before emerging to adults. The female offspring have no effects on the birth process, but involve in the intraspecific competition and hence the density dependence during the aquatic stages. In this case, the total number of transgenic male and female offspring produced per wild mosquito per unit of time from the matings with sterile mosquitoes is $2b$, and the total number of all offspring produced by all wild mosquitoes mating with sterile mosquitoes per unit of time is $\frac{2bc}{S+c}S$. Since female offspring die before emerging to adults, the total number of sterile mosquitoes is the same as in model (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ). Thus the model equation has the following form (29) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S\\ & =(\frac{aS((S-\xi (aS+bc))S+c)}{(S+c{)}^2+\gamma bcS}-\mu )S.\end{array}$(29) The equation solving for positive equilibrium is (30) $\begin{array}{r}{F}_4(S):=a^2\xi S^3-(a-\mu -2ab\xi c)S^2-(a-2\mu -b\gamma \mu )cS+\mu c^2=0.\end{array}$(30) Analysis similar to those in previous sections can be carried out. Or, if we let $b_1=2b$ and ${\gamma }_1=\gamma /2$, model (29(29) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S\\ & =(\frac{aS((S-\xi (aS+bc))S+c)}{(S+c{)}^2+\gamma bcS}-\mu )S.\end{array}$(29) ) is transformed to (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ) with b and γ replaced by $b_1$ and ${\gamma }_1$, respectively. Thus, we can easily obtain the results shown in the theorem below with detailed proofs omitted.

Theorem 2.5

Let the birth rates of wild mosquitoes mating with wild and sterile mosquitoes a and b, the death rate of wild mosquitoes μ, and the emergence rate of sterile larvae γ satisfy the condition $a\ge \frac{8\mu }{2-3\gamma \mu }.$Then there exists a threshold of releases of sterile mosquitoes $c_{EFK}^0$ given by the positive solution of equation $H_4(c)=0$, where $H_4(c)$ is given by (31) $\begin{array}{r}{H}_4(c):=q_3{c}^3+q_2{c}^2+q_{1}c+q_0,\end{array}$(31) with $\begin{array}{rl}{q}_3& =324a^3{b}^3{\xi }^3\mu ,\\ q_2& =a^2{\xi }^2(27a^2{\mu }^2+36ab\mu (a-2\mu -b\gamma \mu )-48b^2\mu (a-\mu )-4b^2(a-2\mu -b\gamma \mu {)}^2),\\ q_1& =2a\xi (-9a\mu (a-\mu )(a-2\mu -b\gamma \mu )+12b\mu (a-\mu {)}^2-2a(a-2\mu -b\gamma \mu {)}^3\\ & +2b(a-\mu )(a-2\mu -b\gamma \mu {)}^2),\\ q_0& =-(a-\mu {)}^2(4\mu (a-\mu )+(a-2\mu -b\gamma \mu {)}^2).\end{array}$If $c<c_{LFK}^0$, model (29(29) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S\\ & =(\frac{aS((S-\xi (aS+bc))S+c)}{(S+c{)}^2+\gamma bcS}-\mu )S.\end{array}$(29) ) has two positive equilibria, $S_{LFK}^{-}<S_{LFK}^{+}$. Equilibrium $S_{LFK}^{-}$ is unstable and $S_{LFK}^{+}$ is locally asymptotically stable. If $c=c_{LFK}^0$, there exists a unique positive equilibrium $S_{LFK}^{\ast }$, which is unstable or semi-stable. If $c>c_{LFK}^0$, there exists no positive equilibrium, and the zero solution is globally asymptotically stable.

3. Discussions

We have formulated simple models for the modified genetic SIT, specifically for the early acting bisex lethality, late acting bisex lethality, early acting female-killing lethality, and late acting female-killing lethality strategies, respectively. The zero equilibrium is always locally asymptotically stable for all models. The models we formulated in this study are one-dimensional differential equation models. However, with the help of the idea that only sexually active sterile mosquitoes are included, as in [20–24], these models have well described and significantly characterized the features of the four genetic SIT strategies. To determine whether the wild mosquitoes can be suppressed or eradicated, we investigated the existence of positive equilibria and their stability. We then found release thresholds of the transgenic sterile mosquitoes explicitly for model (2(2) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS}{S+c}S)S-\mu S.\end{array}$(2) ), and implicitly for models (7(7) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S.\end{array}$(7) ), (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ), and (29(29) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S\\ & =(\frac{aS((S-\xi (aS+bc))S+c)}{(S+c{)}^2+\gamma bcS}-\mu )S.\end{array}$(29) ).

To demonstrate and compare the results for the four models, we give two numerical examples as follows. For convenience, we simply call the four implementations EBS, LBS, EFK, and LFK.

Example 3.1

In this example, we choose parameters (32) $\begin{array}{r}a=200,b=85,\xi =0.01,\mu =0.43\end{array}$(32) for all models. From Theorems 2.1, 2.3, 2.4, and 2.5, we determine the release threshold and then let the release rate c be less than the threshold and find two positive equilibria for each model. The results are summarized in Table 1. Column 3 shows the thresholds, column 4 lists the release rates we choose, and the corresponding two positive equilibria are shown in column 5 and 6, respectively.

Table 1. Summary table of the results for models (2(2) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS}{S+c}S)S-\mu S.\end{array}$(2) ), (7(7) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S.\end{array}$(7) ), (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ), (29(29) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S\\ & =(\frac{aS((S-\xi (aS+bc))S+c)}{(S+c{)}^2+\gamma bcS}-\mu )S.\end{array}$(29) ).

	γ	Thresholds	c	Positive equilibria
EBS	0	$c_{EBS}=15973.0123$	15900	$S_{EBS}^{-}=49.5045$	$S_{EBS}^{+}=53.5369$
LBS	0	$c_{LBS}=68.0541$	65	$S_{LBS}^{-}=0.2311$	$S_{LBS}^{+}=0.3552$
EFK	0.4	$c_{EFK}=116.3435$	110	$S_{EFK}^{-}=0.4138$	$S_{EFK}^{+}=0.6641$
LFK	0.4	$c_{LFK}=6.8497$	6.5	$S_{LFK}^{-}=0.0156$	$S_{LFK}^{+}=0.5246$

From Table 1, it is clear that the release threshold of LBS is less than the release threshold of EPS, and, the release threshold of LFK is less than the release threshold of EFK. It is also true for any parameter settings that the release thresholds of the transgenic sterile mosquitoes for the late acting lethality strategies are less than those for the corresponding early acting lethality strategies.

Indeed, wild mosquitoes mating with sterile mosquitoes carrying an EBS lethal gene produce no offspring and wild mosquitoes mating with sterile mosquitoes carrying an EFK lethal gene produce only male offspring. In contrast, the matings with sterile mosquitoes carrying a LBS or LFK lethal gene produce both male and female offspring. Even those offspring only survive through the aquatic stage in LBS and LFK, they participate in the intraspecific competition which reduces the survivability of wild larvae. Thus, it is clear that less releases of sterile mosquitoes carrying a late acting lethal gene can possibly obtain the same suppression results as the releases of sterile mosquitoes carrying an early acting lethal gene.

This can be also confirmed mathematically. From the equations for the existence of positive equilibrium $F_1(S)=a^2\xi S^3-(a-\mu )S^2-(a-2\mu )cS+\mu c^2=0$and $F_2(S)=a^2\xi S^3-(a-\mu -2ab\xi c)S^2-(a-2\mu )cS+\mu c^2=0,$$F_2(S)=F_1(S)+2ab\xi c{S}^2$. If $c=c_{EBS}^0$, then since the curve of $F_1(S)$ is tangent to and above the S-axis for $S={\overline{S}}_1$ and $S\ne {\overline{S}}_1$, respectively, $F_2(S)>0$ for all $S\ge 0$. Hence, the release threshold $c_{LBS}$ for LBS is less than $c_{EBS}^0$ for EBS. Similarly, it follows from the equations of the existence of positive equilibrium for EFK and LFK, (20(20) $\begin{array}{rl}{F}_3(S)& :=\mu (S^2+(2+\gamma b)cS+c^2)+aS(a\xi S^2-(1-b\xi c)S-c))\\ & =a^2\xi S^3-(a-\mu -ab\xi c)S^2-(a-(2+b\gamma )\mu )cS+\mu c^2=0.\end{array}$(20) ) and (30(30) $\begin{array}{r}{F}_4(S):=a^2\xi S^3-(a-\mu -2ab\xi c)S^2-(a-2\mu -b\gamma \mu )cS+\mu c^2=0.\end{array}$(30) ), $F_4(S)=F_3(S)+ab\xi c{S}^2$. Then the release threshold $c_{LFK}$ for LFK is less than $c_{EFK}^0$ for EFK for any parameter settings.

Notice that, however, the comparison between the late acting bisex and late female-killing lethality strategies is not straightforward. Both male and female offspring are produced and survive through the aquatic stage in LBS, whereas only male offspring are produced and survive through the aquatic stage in EFK. Then the intraspecific competition among the larvae with the LBS lethality implementation is stronger than that with the EFK lethality strategy. On the other hand, the male offspring only produced with the EFK lethality strategy progress to male adults and then contribute to reducing the birth of wild mosquitoes as well. Thus, the early acting female-killing lethality strategy may or may not be more effective than the late acting bisex lethality strategy, depending on other parameters. We show the comparison in Example 3.2.

Example 3.2

With the parameters given in (32(32) $\begin{array}{r}a=200,b=85,\xi =0.01,\mu =0.43\end{array}$(32) ) and $c_{LBS}=68.0541<c=110<c_{EFK}=116.3435$, model (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ) has two positive equilibria given in Table 1, but there is no positive equilibrium for model (7(7) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S.\end{array}$(7) ) and all wild mosquitoes are wiped out with the LBS lethality strategy. We then choose different parameters (33) $\begin{array}{r}a=185,b=90,\xi =0.01,\mu =0.67,\gamma =\mathrm{0.9.}\end{array}$(33) Now thresholds are $c_{EFK}=37.8252<c_{LBS}=38.065$. If we let c = 37.9, model (7(7) $\begin{array}{r}\frac{\mathrm{d}S}{\mathrm{d}t}=\frac{aS}{S+c}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S.\end{array}$(7) ) has two equilibria $S_{LBS}^{-}=0.2585<S_{LBS}^{+}=0.2949$, but model (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ) has no positive equilibrium and hence all wild mosquitoes are wiped out with the EFK lethality strategy.

Therefore, with the parameters given in (32(32) $\begin{array}{r}a=200,b=85,\xi =0.01,\mu =0.43\end{array}$(32) ), the late acting bisex lethality strategy is more effective in the control of mosquitoes, but the early acting female-killing lethality strategy is more effective with parameters given in (33(33) $\begin{array}{r}a=185,b=90,\xi =0.01,\mu =0.67,\gamma =\mathrm{0.9.}\end{array}$(33) ).

We would like to point out that while our model formulation is significant and we have gained insight into the impact of the releases of different transgenic mosquitoes with the modified genetic SIT on mosquito control, more sophisticated models also need to be formulated. More specifically, we have ignored the time delay from the male larvae carrying a female-killing lethality gene emerging to adults in models (19(19) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+bc}{S+c}S)S-\mu S\\ & =\frac{aS}{(S+c{)}^2+\gamma bcS}(S+c-\xi (aS+bc)S)S-\mu S.\end{array}$(19) ) and (29(29) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\frac{aS}{S+c+\frac{\gamma bc}{S+c}S}(1-\xi \frac{aS+2bc}{S+c}S)S-\mu S\\ & =(\frac{aS((S-\xi (aS+bc))S+c)}{(S+c{)}^2+\gamma bcS}-\mu )S.\end{array}$(29) ). It would be more appropriate to have this time delay included in more realistic models. Moreover, it is clear that we need more modelling effects on the intraspecific competition among the wild and sterile larvae, and more specifically on the density dependence in the aquatic stages. We have implicitly included the larvae stage in the birth process in this study, but models with more explicit stage structure need to be formulated and studied to have better understanding for the mosquito control with the modified genetic SIT even more analytic difficulties and challenging are expected. We plan these projects in the near future.

We would also like to point out that because of again the technique difficulty in mathematical analysis, we have assumed $a\ge 2\mu$ for the case of LBS and (27(27) $\begin{array}{r}a\ge \frac{10\mu }{4-3\gamma \mu }.\end{array}$(27) ) for the case of EFK. While these conditions can be satisfied biologically, we should be able to remove them with different mathematical methods. This is also in our future research plan.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is supported by National Natural Science Foundation of China [grant numbers 12331017, 12071095] and the Program for Changjiang Scholars and Innovative Research Team in University [grant number IRT16R16].