Discrete-time model for malaria transmission with constant releases of sterile mosquitoes

ABSTRACT

In this study, we first formulate a baseline discrete-time mathematical model for malaria transmission where the survival function of mosquitoes is of Beverton–Holt type. We then introduce sterile mosquitoes to the baseline model to explore the transmission dynamics with sterile mosquitoes. We derive formulas for the reproductive number $R_0$ of infection and determine the existence and uniqueness of endemic fixed points as well, for the models with or without sterile mosquitoes. We then study the impact of the releases of sterile mosquitoes on the disease transmissions by investigating the effects of varying the release rates of the sterile mosquitoes. We use a numerical example to illustrate our results for all cases and finally give brief discussions of our findings.

KEYWORDS:

• Mathematical modelling
• vector-borne diseases
• malaria transmission
• discrete-time models
• Beverton–Holt survivability
• sterile mosquitoes

AMS SUBJECT CLASSIFICATIONS:

• 39A10
• 39A60
• 92D25
• 92D30
• 92D40

1. Introduction

Mosquito-borne diseases, such as malaria, are a big concern for the public health worldwide. Malaria is a leading cause of death in many developing countries. There are 3.2 billion people, half of the world's population, living in areas at risk of malaria transmission in 106 countries and territories. In 2016, estimated 216 million cases of malaria occurred worldwide and 445,000 people died, mostly children in the African Region. About 1700 cases of malaria are diagnosed in the United States each year [13, 23, 30].

Malaria is transmitted indirectly not from a human to a human but through vectors as blood-feeding mosquitoes. The effective way to prevent the transmission of such diseases is to control mosquitoes and among the biological methods in control of mosquitoes is the sterile insect technique (SIT) [7], in which the natural reproductive process of the target population is disrupted. By chemical or physical methods, male mosquitoes are genetically modified to be sterile despite being sexually active. These sterile male mosquitoes are then released into the environment to mate with the wild female mosquitoes. A wild female mosquito that mates with a sterile male mosquito will either not reproduce or produce eggs that do not hatch. Repeated releases of genetically modified mosquitoes or the releases of a significantly large number of sterile mosquitoes may eventually wipe out a wild mosquito population, although it is, in practice, often more useful to consider controlling the population rather than eradicating it [5, 8, 29].

Mathematical models have played an important role in providing insights into the transmission dynamics of diseases and in influencing the decision making for preventing the diseases. There are many models in the literature for the study of vector-borne diseases and models incorporating sterile mosquitoes, are also formulated for the disease transmission dynamics [1, 9, 17, 18, 20, 28]. However, most of the models are of continuous time, based on differential equations, to take advantage of the rich theoretical foundation of the well-developed dynamical systems theory and bifurcation theory. As time steps are sufficiently small or population sizes are sufficiently large, these models are valid. The results from those continuous-time models have shown their significance in helping understand the transmission dynamics and make useful strategies for control and prevention of the disease.

On the other hand, the timescales for the dynamics of the human and mosquito populations are significantly different. Then, discrete-time models, based on difference equations, become equally important and even more appropriate [3, 10, 11, 14]. Moreover, discrete-time models also have many strengths in the theory of population dynamics and mathematical epidemiology [2, 10, 11, 27, 31].

In this paper, we first formulate a discrete-time susceptible-exposed-infective-recovered (SEIR) compartmental model for humans and a susceptible-exposed-infective (SEI) compartmental model for mosquitoes, based on the model in [21], as our baseline model in Section 2. We derive a formula for the reproductive number $R_0$ of infection and study the existence of endemic fixed points. We then incorporate the sterile mosquitoes into the baseline model such that the susceptible mosquitoes consist of wild and sterile mosquitoes in Section 3. We consider the case where the release rate of sterile mosquitoes is constant and derive a formula for the reproductive number. We investigate the impact of the releases of sterile mosquitoes on the malaria transmission based on the reproductive number and the endemic fixed points, and give a numerical example to verify our theoretical results in Section 4. Brief discussions of our findings are provided in Section 5.

2. Baseline model for malaria transmission

To study the impact of interactive wild and sterile mosquitoes on the malaria transmission, based on the discrete-time malaria model in [21], we formulate a new discrete-time model as our baseline model in this section, where the survivability of mosquitoes assumes the Beverton–Holt type of nonlinearity instead of the Ricker type in [21].

Since the mosquitoes lifespan is shorter than their infective period, we assume that mosquitoes are unable to recover after infection. Then we divide the mosquito population into groups of susceptible, exposed or incubating, and infective individuals, and denote their numbers, at generation n, by $S^v\left(n\right)$, $E^v\left(n\right)$ and $I^v\left(n\right)$, respectively [6, 21].

For the human population, we divide it into groups of susceptible, exposed or incubating, infective and recovered individuals. We let $S^h\left(n\right)$ be the number of susceptible humans, $E^h\left(n\right)$ the number of exposed or incubating humans who are infected but not infectious yet, $I^h\left(n\right)$ the number of infective humans who are infected and also infectious, $R^h\left(n\right)$ the number of humans who are recovered from infection but partly lose their immunity [19, 24–26], and $N^h\left(n\right)=S^h\left(n\right)+E^h\left(n\right)+I^h\left(n\right)+R^h\left(n\right)$, the total number of humans at time n.

Then, the model equations for humans and mosquitoes described in Figure 1 are given by (2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) where the parameters are given in Table 1. We assume that the infection has no effects on the birth and that the intraspecific competition only affects the survivability of the susceptible newborns. We further assume that the survivability of the susceptible newborns, ${\alpha }_b^v\left(N^v\right(n\left)\right)$, has the Beverton–Holt-II type of nonlinearity as (2.2) ${\alpha }_b^v\left(N^v\right(n\left)\right):=\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}.$(2.2)

Figure 1. Schematic diagram for the malaria transmission. The mosquito population is divided into three groups, $S^v$, $E^v$ and $I^v$ and the human population is divided into four groups, $S^h$, $E^h$, $I^h$ and $R^h$. The mosquito infection rate ${\lambda }^v\left(n\right)$ is given in (2.3(2.3) ${\lambda }^v\left(n\right)={\beta }^{h}r\frac{I^h\left(n\right)}{N^h\left(n\right)},$(2.3) ), and the human infection rate ${\lambda }^h\left(n\right)$ is given in (2.5(2.5) ${\lambda }^h\left(n\right)=G\left(L\right(n\left)\right),$(2.5) ).

Table 1. The parameters for the malaria transmission model.

Λ:	The constant input flow of the susceptible humans
${\lambda }^h\left(n\right)$:	The human infection rate with $0\le {\lambda }^h\left(n\right)\le 1$
${\alpha }_k^h$:	$k=1,\dots ,8$. Survival probabilities of susceptible, exposed or incubating, infective and recovered humans, respectively, which are all constants
${\gamma }^h$:	The developing rate of incubating humans becoming infective such that $1/{\gamma }^h$ is the incubating period
${\theta }^h$:	The rate of partial immunity loss, which is a constant
${\lambda }^v\left(n\right)$:	The mosquito infection rate with $0\le {\lambda }^v\left(n\right)\le 1$
${\gamma }^v$:	The rate of incubating individuals becoming infective, which is a constant with $0\le {\gamma }^v\le 1$
$b^v$:	The per capita birth rate of mosquitoes
${\alpha }_k^v$:	$k=1,\dots ,5$. Survival probabilities of susceptible, incubating and infective mosquitoes, respectively, which are all constants
$k_b^v$:	The maximum survival probability for newborns with $0\le k_b^v\le 1$
${\eta }_b^v$:	The density-dependant factor
${\gamma }^v$:	The developing rate of incubating mosquitoes becoming infective. $1/{\gamma }^v$ is the incubating period
${\beta }^h$:	The transmission probability per bite to a susceptible mosquito from an infective human
${\beta }^v$:	The transmission probability per bite to a susceptible human from an infective mosquito
r:	The number of average bites by a single mosquito on all human hosts

To determine the infection rates ${\lambda }^v\left(n\right)$ and ${\lambda }^h\left(n\right)$, we let r be the average number of bites that a single mosquito makes on all human hosts. Then the total number of bites made by susceptible mosquitoes is $r{S}^v$, among which, a proportion of $I^h/N^h$ goes to the infected hosts, causing $r\cdot S^v\cdot \left(I^h/N^h\right)\cdot {\beta }^h={\beta }^{h}r\left(I^h/N^h\right)\cdot S^v$ new infections for mosquitoes, where ${\beta }^h$ is the transmission probability per bite to a susceptible mosquito. Thus the infection rate for mosquitoes is given by (2.3) ${\lambda }^v\left(n\right)={\beta }^{h}r\frac{I^h\left(n\right)}{N^h\left(n\right)},$(2.3) where we assume ${\beta }^{h}r<1$.

Similarly, let $r_1$ be the average number of bites a human host receives from all mosquitoes. Then the total number of bites received by susceptible human hosts is $r_1{S}^h$, among which, a proportion of $I^v/N^v$ is from infected mosquitoes, causing $r_1{S}^h\cdot \left(I^v/N^v\right)\cdot {\beta }^v={\beta }^v{r}_1\left(I^v/N^v\right)S^h$ new infections for hosts, where ${\beta }^v$ is the transmission probability per bite to a susceptible human. Here, however, we need to note that the number of infected mosquitoes can be sufficiently larger than the total number of humans, that is, $I^v\left(n\right)/N^h\left(n\right)>1$. Then, it may lead to the human infection rate greater than 1 even when ${\beta }^{v}r<1$, and thus $S^h(n+1)<0$ from (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ), if we directly define ${\lambda }^h\left(n\right)$ in the same way as for ${\lambda }^v\left(n\right)$. Thus we consider the following factor for the infection rate for humans: (2.4) $L\left(n\right):={\beta }^v{r}_1\frac{I^v\left(n\right)}{N^v\left(n\right)}.$(2.4) Moreover, the total number of bites that all mosquitoes make is the same as the total number of bites all human hosts receive. We then need the following balance constraint: $r{N}^v\left(n\right)=r_1{N}^h\left(n\right),$ and this balance constraint and (2.4(2.4) $L\left(n\right):={\beta }^v{r}_1\frac{I^v\left(n\right)}{N^v\left(n\right)}.$(2.4) ) immediately lead to $L\left(n\right)={\beta }^v\frac{r{N}^v\left(n\right)}{N^h\left(n\right)}\frac{I^v\left(n\right)}{N^v\left(n\right)}={\beta }^{v}r\frac{I^v\left(n\right)}{N^h\left(n\right)}.$ Since the number of infected mosquitoes can be sufficiently larger than the total number of humans, which may lead to the factor $L\left(n\right)$ much greater than 1, we assume that the infection rate for humans is given by (2.5) ${\lambda }^h\left(n\right)=G\left(L\right(n\left)\right),$(2.5) where G is a positive function of L, satisfying the following conditions [10]: (2.6) $G:[0,\mathrm{\infty })\to [0,1],G\left(0\right)=0,G^{\prime }\left(L\right)>0,G^{″}\left(L\right)<0,\mathrm{\forall }L\ge 0,\underset{L\to \mathrm{\infty }}{lim}G\left(L\right)=1.$(2.6)

2.1. The reproductive number

We derive a formula for the reproductive number by investigating the local stability of the infection-free fixed point. The Jacobian matrix of system (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ) at the infection-free fixed point $(E^v,I^v,E^h,I^h,S^v,S^h,R^h)=(0,0,0,0,S_0^v,S_0^h,0)$ with (2.7) $S_0^v=\frac{k_b^v{b}^v-(1-{\alpha }_1^v)}{{\eta }_b^v(1-{\alpha }_1^v)},S_0^h=\frac{\mathrm{\Lambda }}{1-{\alpha }_1^h},$(2.7) has the form (2.8) $J_1:=\left(\begin{array}{cc}F+T& 0\\ \cdot & C\end{array}\right),$(2.8) where F is the fertility matrix given by $F=\left(\begin{array}{cccc}0& 0& 0& r{\beta }^h{\alpha }_2^v\frac{s_0^v}{s_0^h}\\ 0& 0& 0& 0\\ 0& G^{\prime }\left(0\right)r{\beta }^v{\alpha }_2^h& 0& 0\\ 0& 0& 0& 0\end{array}\right),$ T is the transition matrix given by $T=\left(\begin{array}{cccc}(1-{\gamma }^v){\alpha }_3^v& 0& 0& 0\\ {\gamma }^v{\alpha }_4^v& {\alpha }_5^v& 0& 0\\ 0& 0& (1-{\gamma }^h){\alpha }_3^h& 0\\ 0& 0& {\gamma }^h{\alpha }_4^h& (1-{\eta }^h){\alpha }_5^h\end{array}\right),$ and matrix C has the form $C=\left(\begin{array}{ccc}{c}_{11}& 0& 0\\ 0& {\alpha }_1^h& {\theta }^h{\alpha }_8^h\\ 0& 0& (1-{\theta }^h){\alpha }_7^h\end{array}\right),$ where $c_{11}=\frac{k_b^v{b}^v}{1+{\eta }_b^v{S}_0^v}-\frac{k_b^v{b}^v{\eta }_b^v{S}_0^v}{(1+{\eta }_b^v{S}_0^v{)}^2}+{\alpha }_1^v<1,$ since $c_{11}-1=(1-{\alpha }_1^v)\left(\frac{1-{\alpha }_1^v}{k_b^v{b}^v}-1\right)<0.$ Then, the next generation matrix of system (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ) [4, 12, 15, 16] is given by $Q=F(I-T{)}^{-1}=\left(\begin{array}{cccc}0& 0& q_{13}& q_{14}\\ 0& 0& 0& 0\\ q_{31}& q_{32}& 0& 0\\ 0& 0& 0& 0\end{array}\right),$ where $\begin{array}{rl}{q}_{13}& :=\frac{r{\beta }^h{\alpha }_4^h{\alpha }_2^v{S}_0^v}{S_0^h(1-(1-{\gamma }^h\left){\alpha }_3^h\right)(1-(1-{\eta }^h\left){\alpha }_5^h\right)},q_{14}:=\frac{r{\beta }^h{\alpha }_2^v{S}_0^v}{S_0^h(1-(1-{\eta }^h\left){\alpha }_5^h\right)},\\ q_{31}& :=\frac{G^{\prime }\left(0\right)r{\alpha }_2^h{\beta }^v{\gamma }^v{\alpha }_4^v}{(1-{\alpha }_5^v)(1-(1-{\gamma }^v\left){\alpha }_3^v\right)},q_{32}:=\frac{G^{\prime }\left(0\right)r{\alpha }_2^h{\beta }^v}{1-{\alpha }_5^v}.\end{array}$

The characteristic polynomial of Q is $\det (\rho I-Q)={\rho }^2({\rho }^2-q_{13}{q}_{31})={\rho }^2({\rho }^2-q_{13}{q}_{32}),$ which leads to the eigenvalues as ${\rho }_1=0$ and ${\rho }_2=\sqrt{q_{13}{q}_{32}}$. According to [16], the infection-free fixed point is locally asymptotically stable if and only if the eigenvalues ${\rho }_{1,2}<1$. Thus we define the net reproductive number of infection as $R_0:={\rho }_2=\sqrt{q_{13}{q}_{32}}$.

Hence, (2.9) $R_0:=\sqrt{\frac{G^{\prime }\left(0\right)r{\beta }^h{\gamma }^h{\alpha }_2^h{\alpha }_4^h}{(1-(1-{\gamma }^h\left){\alpha }_3^h\right)(1-(1-{\eta }^h\left){\alpha }_5^h\right)}\cdot \frac{r{\beta }^v{\gamma }^v{\alpha }_2^v{\alpha }_4^v}{(1-{\alpha }_5^v)(1-(1-{\gamma }^v\left){\alpha }_3^v\right)}\cdot \frac{S_0^v}{S_0^h}},$(2.9) with $S_0^v$ and $S_0^h$ given in (2.7(2.7) $S_0^v=\frac{k_b^v{b}^v-(1-{\alpha }_1^v)}{{\eta }_b^v(1-{\alpha }_1^v)},S_0^h=\frac{\mathrm{\Lambda }}{1-{\alpha }_1^h},$(2.7) ).

We have the following results.

Theorem 2.1

Define the reproductive number of infection, $R_0,$ in (2.9(2.9) $R_0:=\sqrt{\frac{G^{\prime }\left(0\right)r{\beta }^h{\gamma }^h{\alpha }_2^h{\alpha }_4^h}{(1-(1-{\gamma }^h\left){\alpha }_3^h\right)(1-(1-{\eta }^h\left){\alpha }_5^h\right)}\cdot \frac{r{\beta }^v{\gamma }^v{\alpha }_2^v{\alpha }_4^v}{(1-{\alpha }_5^v)(1-(1-{\gamma }^v\left){\alpha }_3^v\right)}\cdot \frac{S_0^v}{S_0^h}},$(2.9) ) for system (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ). Then, the infection-free fixed point is locally asymptotically stable if $R_0<1,$ and is unstable if $R_0>1$.

2.2. The endemic fixed point

We next explore the existence of endemic fixed points of system (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ).

The components of an endemic fixed point need to satisfy the following equations: (2.10a) $S^h=\mathrm{\Lambda }+(1-{\lambda }^h)S^h{\alpha }_1^h+{\theta }^h{R}^h{\alpha }_8^h,$(2.10a) (2.10b) $E^h={\lambda }^h{S}^h{\alpha }_2^h+(1-{\gamma }^h)E^h{\alpha }_3^h,$(2.10b) (2.10c) $I^h={\gamma }^h{E}^h{\alpha }_4^h+(1-{\eta }^h)I^h{\alpha }_5^h,$(2.10c) (2.10d) $R^h={\eta }^h{I}^h{\alpha }_6^h+(1-{\theta }^h)R^h{\alpha }_7^h,$(2.10d) (2.10e) $S^v=\frac{k_b^v{b}^v{N}^v}{1+{\eta }_b^v{N}^v}+(1-{\lambda }^v)S^v{\alpha }_1^v,$(2.10e) (2.10f) $E^v={\lambda }^v{S}^v{\alpha }_2^v+(1-{\gamma }^v)E^v{\alpha }_3^v,$(2.10f) (2.10g) $I^v={\gamma }^v{E}^v{\alpha }_4^v+I^v{\alpha }_5^v.$(2.10g)

Solving (2.10a(2.10a) $S^h=\mathrm{\Lambda }+(1-{\lambda }^h)S^h{\alpha }_1^h+{\theta }^h{R}^h{\alpha }_8^h,$(2.10a) )–(2.10d(2.10d) $R^h={\eta }^h{I}^h{\alpha }_6^h+(1-{\theta }^h)R^h{\alpha }_7^h,$(2.10d) ), we have (2.11) $\begin{array}{rl}{E}^h& ={\omega }_1^h{\lambda }^h{S}^h,I^h={\omega }_1^h{\omega }_2^h{\lambda }^h{S}^h,\\ R^h& ={\omega }_1^h{\omega }_2^h{\omega }_3^h{\lambda }^h{S}^h,S^h=\frac{\mathrm{\Lambda }}{1-{\alpha }_1^h+A_1^h{\lambda }^h},\end{array}$(2.11) where $\begin{array}{rl}{\omega }_1^h& :=\frac{{\alpha }_2^h}{1-(1-{\gamma }^h){\alpha }_3^h},{\omega }_2^h:=\frac{{\gamma }^h{\alpha }_4^h}{1-(1-{\eta }^h){\alpha }_5^h},\\ {\omega }_3^h& :=\frac{{\eta }^h{\alpha }_6^h}{1-(1-{\theta }^h){\alpha }_7^h},A_1^h:={\alpha }_1^h-{\theta }^h{\alpha }_8^h{\omega }_1^h{\omega }_2^h{\omega }_3^h.\end{array}$

Then we have (2.12) $N^h=\frac{\mathrm{\Lambda }(1+A_2^h{\lambda }^h)}{1-{\alpha }_1^h+A_1^h{\lambda }^h},$(2.12) where $A_2^h:={\omega }_1^h+{\omega }_1^h{\omega }_2^h+{\omega }_1^h{\omega }_2^h{\omega }_3^h$.

Solving (2.10e(2.10e) $S^v=\frac{k_b^v{b}^v{N}^v}{1+{\eta }_b^v{N}^v}+(1-{\lambda }^v)S^v{\alpha }_1^v,$(2.10e) )–(2.10g(2.10g) $I^v={\gamma }^v{E}^v{\alpha }_4^v+I^v{\alpha }_5^v.$(2.10g) ), we have (2.13) $S^v=\frac{k_b^v{b}^v{N}^v}{(1-(1-{\lambda }^v\left){\alpha }_1^v\right)(1+{\eta }_b^v{N}^v)},E^v={\omega }_1^v{\lambda }^v{S}^v,I^v={\omega }_1^v{\omega }_2^v{\lambda }^v{S}^v,$(2.13) where (2.14) ${\omega }_1^v:=\frac{{\alpha }_2^v}{1-(1-{\gamma }^v){\alpha }_3^v},{\omega }_2^v:=\frac{{\gamma }^v{\alpha }_2^v{\alpha }_4^v}{(1-{\alpha }_5^v)(1-(1-{\gamma }^v\left){\alpha }_3^v\right)}.$(2.14) Substituting (2.11(2.11) $\begin{array}{rl}{E}^h& ={\omega }_1^h{\lambda }^h{S}^h,I^h={\omega }_1^h{\omega }_2^h{\lambda }^h{S}^h,\\ R^h& ={\omega }_1^h{\omega }_2^h{\omega }_3^h{\lambda }^h{S}^h,S^h=\frac{\mathrm{\Lambda }}{1-{\alpha }_1^h+A_1^h{\lambda }^h},\end{array}$(2.11) )–(2.13(2.13) $S^v=\frac{k_b^v{b}^v{N}^v}{(1-(1-{\lambda }^v\left){\alpha }_1^v\right)(1+{\eta }_b^v{N}^v)},E^v={\omega }_1^v{\lambda }^v{S}^v,I^v={\omega }_1^v{\omega }_2^v{\lambda }^v{S}^v,$(2.13) ) into (2.3(2.3) ${\lambda }^v\left(n\right)={\beta }^{h}r\frac{I^h\left(n\right)}{N^h\left(n\right)},$(2.3) ) and (2.5(2.5) ${\lambda }^h\left(n\right)=G\left(L\right(n\left)\right),$(2.5) ), respectively, we obtain (2.15) ${\lambda }^v={\beta }^{h}r\frac{I^h}{N^h}={\beta }^{h}r\frac{{\omega }_1^h{\omega }_2^h{\lambda }^h}{1+A_2^h{\lambda }^h}$(2.15) and $\begin{array}{rl}L\left({\lambda }^h\right)& ={\beta }^{v}r\frac{I^v}{N^h}=\frac{{\beta }^{v}r{\omega }_1^v{\omega }_2^v{\alpha }_2^v{S}^v(1-{\alpha }_1^h+A_1^h{\lambda }^h)}{\mathrm{\Lambda }(1+A_2^h{\lambda }^h)}{\lambda }^v\\ & =\frac{{\beta }^{v}r{\omega }_1^v{\omega }_2^v{\alpha }_2^v{S}^v(1-{\alpha }_1^h+A_1^h{\lambda }^h)}{\mathrm{\Lambda }(1+A_2^h{\lambda }^h)}\cdot \frac{{\beta }^{h}r{\omega }_1^h{\omega }_2^h{\lambda }^h}{1+A_2^h{\lambda }^h}\\ & =R_0^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h{)}^2}\frac{S^v\left({\lambda }^h\right)}{S_0^v}{\lambda }^h.\end{array}$ Hence, there exists an endemic fixed point if and only if there is a positive solution to equation ${\lambda }^h=G\left(L\right({\lambda }^h\left)\right)$, or equivalently, $H\left({\lambda }^h\right):=G\left(L\right({\lambda }^h\left)\right)-{\lambda }^h$ for $0\le {\lambda }^h\le 1$.

Notice that $L\left(0\right)=0$ and hence $G\left(L\right(0\left)\right)=0$ by (2.6(2.6) $G:[0,\mathrm{\infty })\to [0,1],G\left(0\right)=0,G^{\prime }\left(L\right)>0,G^{″}\left(L\right)<0,\mathrm{\forall }L\ge 0,\underset{L\to \mathrm{\infty }}{lim}G\left(L\right)=1.$(2.6) ). Then $H\left(0\right)=0$ and $H\left(1\right)<0$ since $G^{\prime }\left(L\right)>0$ and $\underset{L\to \mathrm{\infty }}{lim}G\left(L\right)=1$. Moreover, it follows from $S^v\left(0\right)=S_0^v$ and $L^{\prime }\left(0\right)=R_0^2/G^{\prime }\left(0\right)$ that $H^{\prime }\left(0\right)=G^{\prime }\left(0\right)L^{\prime }\left(0\right)-1=R_0^2-1.$ Thus if $R_0>1$, there exists a positive solution to $H\left({\lambda }^h\right)=0$, that is, to ${\lambda }^h=G\left(L\right({\lambda }^h\left)\right)$. In other words, there exists an endemic fixed point if $R_0>1$.

Theorem 2.2

For system (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ), if the reproductive number $R_0>1,$ there exits an endemic fixed point.

3. The interactive transmission model with sterile mosquitoes

Now suppose that sterile mosquitoes are released into a wild mosquito population and we let $B\left(n\right)$ be the number of the sterile mosquitoes released at time n. Since sterile mosquitoes do not reproduce, there is no maturation process from larvae to adults for sterile mosquitoes. Hence, the number of sterile mosquitoes at n is just the number of released sterile mosquitoes, and the size of total mosquitoes is $N^v\left(n\right)+B\left(n\right)$. We assume sterile mosquitoes are constantly released such that $B\left(n\right):=b>0$ is a positive constant. After the sterile mosquitoes are released, the mating interaction between the wild and sterile mosquitoes takes place. The number of susceptible offspring produced per mating is $\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}.$ Then the model for the interactive mosquitoes and the human population can be described as (3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) where the human part remains the same as described in system (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ).

Assume that all of the survival probabilities of mosquitoes ${\alpha }_i^v$, $i=1,2,3,4,5$, in system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ) are the same, denoted by ${\alpha }^v$. Then the interactive dynamics between the total wild mosquitoes and the constant sterile mosquitoes are governed by the following equation: (3.2) $N^v(n+1)=\frac{k_b^v{b}^v\left(N^v\right(n){)}^2}{(1+{\eta }_b^v{N}^v(n\left)\right)\left(N^v\right(n)+b)}+{\alpha }^v{N}^v\left(n\right).$(3.2) It follows from [22] that the results for (3.2(3.2) $N^v(n+1)=\frac{k_b^v{b}^v\left(N^v\right(n){)}^2}{(1+{\eta }_b^v{N}^v(n\left)\right)\left(N^v\right(n)+b)}+{\alpha }^v{N}^v\left(n\right).$(3.2) ) can be stated as follows.

Lemma 3.1

Theorem 4.1 [22]

For given $r_0=k_b^v{b}^v,$ system (3.2(3.2) $N^v(n+1)=\frac{k_b^v{b}^v\left(N^v\right(n){)}^2}{(1+{\eta }_b^v{N}^v(n\left)\right)\left(N^v\right(n)+b)}+{\alpha }^v{N}^v\left(n\right).$(3.2) ) has no, one or two positive fixed points if $b>b_c,$ $b=b_c$ or $b<b_c,$ respectively, where the release threshold $b_c$ is given by (3.3) $b_c:=\frac{(\sqrt{k_b^v{b}^v}-\sqrt{1-{\alpha }^v}{)}^2}{{\eta }_b^v(1-{\alpha }^v)}.$(3.3) If $b>b_c,$ there exists no positive fixed point and $N_0^{vb}=0$ is the only fixed point, which is globally asymptotically stable. Thus the wild mosquito population goes extinct regardless of the initial population size. If $b=b_c,$ there exists a unique positive fixed point which is unstable and thus the trivial fixed point $N_0^{vb}=0$ is also globally asymptotically stable. If $b<b_c,$ there exist two positive fixed points (3.4) $\begin{array}{rl}& N^{vb(\pm )}\\ & =\frac{k_b^v{b}^v-(1+b{\eta }_b^v)(1-{\alpha }^v)\pm \sqrt{(k_b^v{b}^v-(1+b{\eta }_b^v\left)\right(1-{\alpha }^v){)}^2-4b{\eta }_b^v(1-{\alpha }^v{)}^2}}{2{\eta }_b^v(1-{\alpha }^v)},\end{array}$(3.4) where $N^{vb(-)}$ is unstable and $N^{vb(+)}$ is locally asymptotically stable with a basin of attraction $(N^{vb(-)},\mathrm{\infty }),$ while the trivial fixed point $N_0^{vb}=0$ is locally asymptotically stable with a basin of attraction $(0,N^{vb(-)})$.

Notice that if $b>b_c$, for system (3.2(3.2) $N^v(n+1)=\frac{k_b^v{b}^v\left(N^v\right(n){)}^2}{(1+{\eta }_b^v{N}^v(n\left)\right)\left(N^v\right(n)+b)}+{\alpha }^v{N}^v\left(n\right).$(3.2) ), there exists no positive fixed point and the trivial fixed point $N_0^{vb}=0$ is the only fixed point, which is globally asymptotically stable. Correspondingly, for system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ), all wild mosquitoes eventually go extinct so that there will be no malaria transmission. Thus we only consider the case of $b<b_c$, and we assume the initial size of the wild mosquitoes is greater than $N^{vb(-)}$, as shown in Figure 2.

Figure 2. Solutions approach $N^v=0$ or $N^{vb(+)}$ depending on their initial values.

3.1. The reproductive number and disease spread

We derive a formula for the reproductive number of infection after the sterile mosquitoes are released into the wild mosquito population with $b<b_c$ and the initial size in $(N^{vb(-)},\mathrm{\infty })$ as follows.

Consider the following infection-free fixed point of system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ): $(E^v,I^v,E^h,I^h,S^v,S^h,R^h)=(0,0,0,0,S_0^{vb},S_0^{hb},0).$ It follows from (2.7(2.7) $S_0^v=\frac{k_b^v{b}^v-(1-{\alpha }_1^v)}{{\eta }_b^v(1-{\alpha }_1^v)},S_0^h=\frac{\mathrm{\Lambda }}{1-{\alpha }_1^h},$(2.7) ) that $S_0^{hb}=S_0^h$ and from (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ) that (3.5) $S_0^{vb}=\frac{k_b^v{S}_0^{vb}}{1+{\eta }_b^v{S}_0^{vb}}\cdot \frac{b^v{S}_0^{vb}}{S_0^{vb}+b}+S_0^{vb}{\alpha }_1^v.$(3.5) We then define function $P\left(b\right):=\mathrm{\Phi }\left(S_0^{vb}\right(b),b)$ by (3.6) $\mathrm{\Phi }\left(S_0^{vb}\right(b),b):=(1-{\alpha }_1^v){\eta }_b^v{S_0^{vb}}^2+\left[(1-{\alpha }_1^v)(1+b{\eta }_b^v)-k_b^v{b}^v\right]S_0^{vb}+b(1-{\alpha }_1^v).$(3.6) Clearly, there exists no, one or two positive roots of (3.6(3.6) $\mathrm{\Phi }\left(S_0^{vb}\right(b),b):=(1-{\alpha }_1^v){\eta }_b^v{S_0^{vb}}^2+\left[(1-{\alpha }_1^v)(1+b{\eta }_b^v)-k_b^v{b}^v\right]S_0^{vb}+b(1-{\alpha }_1^v).$(3.6) ) if $b>b_c$, $b=b_c$ or $b<b_c$, respectively, where $b_c$ is given in (3.3(3.3) $b_c:=\frac{(\sqrt{k_b^v{b}^v}-\sqrt{1-{\alpha }^v}{)}^2}{{\eta }_b^v(1-{\alpha }^v)}.$(3.3) ). Therefore, in the case of $b<b_c$, there exist two infection-free fixed points.

The Jacobian matrix $J_b$ evaluated at an infection-free fixed point has the following form: $J_b:=\left(\begin{array}{cc}{F}_b+T_b& 0\\ \cdot & C_b\end{array}\right),$ where $F_b$ is the fertility matrix $F_b=F$, $T_b$ is the transition matrix $T_b=T$, with F and T given in (2.8(2.8) $J_1:=\left(\begin{array}{cc}F+T& 0\\ \cdot & C\end{array}\right),$(2.8) ), and $C_b=\left(\begin{array}{ccc}{c}_{11}^b& 0& 0\\ 0& {\alpha }_1^h& {\theta }^h{\alpha }_8^h\\ 0& 0& (1-{\theta }^h){\alpha }_7^h\end{array}\right)$ with $\begin{array}{rl}{c}_{11}^b& =\frac{k_b^v{b}^v{S}_0^{vb}\left[2\right(1+{\eta }_b^v{S}_0^{vb}\left)\right(b+S_0^{vb})-S_0^{vb}(1+2{\eta }_b^v{S}_0^{vb}+b{\eta }_b^v\left)\right]}{(1+{\eta }_b^v{S}_0^{vb}{)}^2(b+S_0^{vb}{)}^2}\\ & =2-\frac{S_0^{vb}(1+2{\eta }_b^v{S}_0^{vb}+b{\eta }_b^v)}{(1+{\eta }_b^v{S}_0^{vb})(b+S_0^{vb})}\\ & =2-\frac{(1-{\alpha }_1^v)(1+2{\eta }_b^v{S}_0^{vb}+b{\eta }_b^v)}{k_b^v{b}^v}.\end{array}$ It follows from (3.6(3.6) $\mathrm{\Phi }\left(S_0^{vb}\right(b),b):=(1-{\alpha }_1^v){\eta }_b^v{S_0^{vb}}^2+\left[(1-{\alpha }_1^v)(1+b{\eta }_b^v)-k_b^v{b}^v\right]S_0^{vb}+b(1-{\alpha }_1^v).$(3.6) ) that $\begin{array}{rl}{\mathrm{\partial }}_{S_0^{vb}}\mathrm{\Phi }\left(S_0^{vb}\right(b),b)& =2(1-{\alpha }_1^v){\eta }_b^v{S}_0^{vb}+(1-{\alpha }_1^v)(1+b{\eta }_b^v)-k_b^v{b}^v\\ & =(1-{\alpha }_1^v)(1+2{\eta }_b^v{S}_0^{vb}+b{\eta }_b^v)-k_b^v{b}^v,\end{array}$ which implies that, for the two positive roots ${S_0^{vb}}^{(-)}<{S_0^{vb}}^{(+)}$, ${\mathrm{\partial }}_{{S_0^{vb}}^{(-)}}\mathrm{\Phi }<0,{\mathrm{\partial }}_{{S_0^{vb}}^{(+)}}\mathrm{\Phi }>0;$ that is, $\frac{(1-{\alpha }_1^v)(1+2{\eta }_b^v{S_0^{vb}}^{(-)}+b{\eta }_b^v)}{k_b^v{b}^v}<1,\frac{(1-{\alpha }_1^v)(1+2{\eta }_b^v{S_0^{vb}}^{(+)}+b{\eta }_b^v)}{k_b^v{b}^v}>1.$ Thus we have $c_{11}^b>1$ for ${S_0^{vb}}^{(-)}$, and hence ${S_0^{vb}}^{(-)}$ is unstable.

On the other hand, it follows from $c_{11}^b<1$ for ${S_0^{vb}}^{(+)}$ that all eigenvalues of matrix C are inside the unit circle, and thus the local stability of this infection-free fixed point is determined by matrix $F_b+T_b$. It follows then from Section 2.1 that we can define the reproductive number of infection for system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ) as (3.7) $R_0^b\left(b\right):=\sqrt{\frac{G^{\prime }\left(0\right)r{\beta }^h{\gamma }^h{\alpha }_2^h{\alpha }_4^h}{(1-(1-{\gamma }^h\left){\alpha }_3^h\right)(1-(1-{\eta }^h\left){\alpha }_5^h\right)}\cdot \frac{r{\beta }^v{\gamma }^v{\alpha }_2^v{\alpha }_4^v}{(1-{\alpha }_5^v)(1-(1-{\gamma }^v\left){\alpha }_3^v\right)}\cdot \frac{S_0^{vb(+)}\left(b\right)}{S_0^h}},$(3.7) where $S_0^{vb(+)}\left(b\right)$ equals $N^{vb(+)}$ given in (3.4(3.4) $\begin{array}{rl}& N^{vb(\pm )}\\ & =\frac{k_b^v{b}^v-(1+b{\eta }_b^v)(1-{\alpha }^v)\pm \sqrt{(k_b^v{b}^v-(1+b{\eta }_b^v\left)\right(1-{\alpha }^v){)}^2-4b{\eta }_b^v(1-{\alpha }^v{)}^2}}{2{\eta }_b^v(1-{\alpha }^v)},\end{array}$(3.4) ).

By using b as a variable, $S_0^{vb(+)}\left(b\right)$ is a function of b and so is $R_0^b$. When b=0, it is clear that $S_0^{vb(+)}\left(b\right)=S_0^v$ and $R_0^b=R_0$. Then $R_0^b\left(b\right)=\sqrt{\frac{S_0^{vb(+)}\left(b\right)}{S_0^v}}{R}_0,$ and, for $0<b<b_c$, $0<S_0^{vb(+)}\left(b\right)<S_0^v$. Thus $R_0^b<R_0$.

It follows from (3.6(3.6) $\mathrm{\Phi }\left(S_0^{vb}\right(b),b):=(1-{\alpha }_1^v){\eta }_b^v{S_0^{vb}}^2+\left[(1-{\alpha }_1^v)(1+b{\eta }_b^v)-k_b^v{b}^v\right]S_0^{vb}+b(1-{\alpha }_1^v).$(3.6) ) that $P^{\prime }\left(b\right)={\mathrm{\partial }}_{S_0^{vb}}\mathrm{\Phi }\left(S_0^{vb}\right(b),b)\cdot {S_0^{vb}}^{\prime }\left(b\right)+(1-{\alpha }_1^v)(1+{\eta }_b^v{S}_0^{vb})=0,$ and we have (3.8) ${S_0^{vb}}^{\prime }\left(b\right)=-\frac{(1-{\alpha }_1^v)(1+{\eta }_b^v{S}_0^{vb})}{{\mathrm{\partial }}_{S_0^{vb}}\mathrm{\Phi }\left(S_0^{vb}\right(b),b)}<0.$(3.8) Since ${\mathrm{\partial }}_{{S_0^{vb}}^{(+)}}\mathrm{\Phi }>0$, ${S_0^{vb}}^{(+)}$ is monotone decreasing with respect to b. Notice that function $R_0^b$ is monotone increasing with respect to ${S_0^v}^{(+)}$. That is to say, the composed function $R_0^b$ is monotone decreasing with respect to b. Hence, there is a unique threshold $\overline{b}$ such that $R_0^b\left(\overline{b}\right)=1$ and $R_0^b\left(b\right)\left\{\begin{array}{ll}>1,& \mathrm{if}b<\overline{b},\\ <1,& \mathrm{if}b>\overline{b}.\end{array}\right.$ In fact, the threshold $\overline{b}$ can be explicitly solved as follows.

Since $1=\sqrt{S_0^{vb(+)}\left(\overline{b}\right)/S_0^v}{R}_0$, we have $S_0^{vb(+)}\left(\overline{b}\right)=\frac{S_0^v}{R_0^2}.$ It follows from (3.5(3.5) $S_0^{vb}=\frac{k_b^v{S}_0^{vb}}{1+{\eta }_b^v{S}_0^{vb}}\cdot \frac{b^v{S}_0^{vb}}{S_0^{vb}+b}+S_0^{vb}{\alpha }_1^v.$(3.5) ) that (3.9) $\overline{b}=\frac{k_b^v{b}^v{S}_0^{vb}}{(1-{\alpha }_1^v)(1+{\eta }_b^v{S}_0^{vb})}-S_0^{vb}=\frac{k_b^v{b}^v{S}_0^v}{(1-{\alpha }_1^v)(R_0^2+{\eta }_b^v{S}_0^v)}-\frac{S_0^v}{R_0^2}.$(3.9) Thus, the infection-free fixed point of system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ) is locally asymptotically stable if $b>\overline{b}$ and unstable if $b<\overline{b}$. We also notice that if $b=b_c$, the infection-free fixed point with the unique positive component $S_0^{vb(+)}$ is unstable and if $b>b_c$, there exists no such positive component $S_0^{vb(+)}$. Thus all wild mosquitoes will be wiped out if $b\ge b_c$ and hence $\overline{b}<b_c$. We summarize the results as follows.

Theorem 3.1

Assume sterile mosquitoes are released into the wild mosquito population constantly with the number of releases b>0. Define the two threshold values of releases $b_c$ and $\overline{b}$ in (3.3(3.3) $b_c:=\frac{(\sqrt{k_b^v{b}^v}-\sqrt{1-{\alpha }^v}{)}^2}{{\eta }_b^v(1-{\alpha }^v)}.$(3.3) ) and (3.9(3.9) $\overline{b}=\frac{k_b^v{b}^v{S}_0^{vb}}{(1-{\alpha }_1^v)(1+{\eta }_b^v{S}_0^{vb})}-S_0^{vb}=\frac{k_b^v{b}^v{S}_0^v}{(1-{\alpha }_1^v)(R_0^2+{\eta }_b^v{S}_0^v)}-\frac{S_0^v}{R_0^2}.$(3.9) ), respectively. Then we have the following results.

• If $b>b_c,$ there exists no positive fixed point for the interactive mosquitoes system (3.2(3.2) $N^v(n+1)=\frac{k_b^v{b}^v\left(N^v\right(n){)}^2}{(1+{\eta }_b^v{N}^v(n\left)\right)\left(N^v\right(n)+b)}+{\alpha }^v{N}^v\left(n\right).$(3.2) ) and the only trivial fixed point $N_0^{vb}=0$ of system (3.2(3.2) $N^v(n+1)=\frac{k_b^v{b}^v\left(N^v\right(n){)}^2}{(1+{\eta }_b^v{N}^v(n\left)\right)\left(N^v\right(n)+b)}+{\alpha }^v{N}^v\left(n\right).$(3.2) ) is globally asymptotically stable. Correspondingly, for system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ), all wild mosquitoes are wiped out and there will be no infection.

• If $b=b_c,$ the unique positive fixed point of system (3.2(3.2) $N^v(n+1)=\frac{k_b^v{b}^v\left(N^v\right(n){)}^2}{(1+{\eta }_b^v{N}^v(n\left)\right)\left(N^v\right(n)+b)}+{\alpha }^v{N}^v\left(n\right).$(3.2) ) is unstable and the trivial fixed point $N_0^{vb}=0$ of system (3.2(3.2) $N^v(n+1)=\frac{k_b^v{b}^v\left(N^v\right(n){)}^2}{(1+{\eta }_b^v{N}^v(n\left)\right)\left(N^v\right(n)+b)}+{\alpha }^v{N}^v\left(n\right).$(3.2) ) is globally asymptotically stable. Correspondingly, for system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ), all wild mosquitoes are wiped out as well and hence there is no infection.

• If $\overline{b}<b<b_c,$ the sterile and wild mosquitoes coexist, but the reproductive number $R_0^b<1$ and the infection-free fixed point of system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ) associated with the locally asymptotically stable positive fixed point $N^{vb(+)}$ given in (3.4(3.4) $\begin{array}{rl}& N^{vb(\pm )}\\ & =\frac{k_b^v{b}^v-(1+b{\eta }_b^v)(1-{\alpha }^v)\pm \sqrt{(k_b^v{b}^v-(1+b{\eta }_b^v\left)\right(1-{\alpha }^v){)}^2-4b{\eta }_b^v(1-{\alpha }^v{)}^2}}{2{\eta }_b^v(1-{\alpha }^v)},\end{array}$(3.4) ) of system (3.2(3.2) $N^v(n+1)=\frac{k_b^v{b}^v\left(N^v\right(n){)}^2}{(1+{\eta }_b^v{N}^v(n\left)\right)\left(N^v\right(n)+b)}+{\alpha }^v{N}^v\left(n\right).$(3.2) ) is locally asymptotically stable. Thus the infection will eventually go extinct.

• If $b<\overline{b}$ such that $R_0^b>1,$ then the infection-free fixed point of system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ) is unstable. The disease spreads when the initial size of wild mosquito is greater than $N^{vb(-)}$ given in (3.4(3.4) $\begin{array}{rl}& N^{vb(\pm )}\\ & =\frac{k_b^v{b}^v-(1+b{\eta }_b^v)(1-{\alpha }^v)\pm \sqrt{(k_b^v{b}^v-(1+b{\eta }_b^v\left)\right(1-{\alpha }^v){)}^2-4b{\eta }_b^v(1-{\alpha }^v{)}^2}}{2{\eta }_b^v(1-{\alpha }^v)},\end{array}$(3.4) ).

3.2. Endemic fixed point

Similarly as in Section 2.2, we determine the existence of endemic fixed points for system (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ) as follows.

The components of an endemic fixed point for the human part are the same as those in Section 2.2 for the baseline model. The components for the wild mosquitoes at an endemic fixed point satisfy the following system: (3.10a) $S^{vb}=\frac{k_b^v}{1+{\eta }_b^v{N}^{vb}}{N}^{vb}\frac{b^v{N}^{vb}}{N^{vb}+b}+(1-{\lambda }^v)S^{vb}{\alpha }_1^v,$(3.10a) (3.10b) $E^{vb}={\lambda }^v{S}^{vb}{\alpha }_2^v+(1-{\gamma }^v)E^{vb}{\alpha }_3^v,$(3.10b) (3.10c) $I^{vb}={\gamma }^v{E}^{vb}{\alpha }_4^v+I^{vb}{\alpha }_5^v,$(3.10c) which leads to $N^{vb}=\frac{k_b^v{b}^v{N^{vb}}^2}{(1+{\eta }_b^v{N}^{vb})(b+N^{vb})}+\alpha N^{vb}$ with ${\alpha }_i^v=\alpha$ for $i=1,\dots ,5$. Then $N^{vb}=S_0^{vb}=N^{vb(+)}$ given in (3.4(3.4) $\begin{array}{rl}& N^{vb(\pm )}\\ & =\frac{k_b^v{b}^v-(1+b{\eta }_b^v)(1-{\alpha }^v)\pm \sqrt{(k_b^v{b}^v-(1+b{\eta }_b^v\left)\right(1-{\alpha }^v){)}^2-4b{\eta }_b^v(1-{\alpha }^v{)}^2}}{2{\eta }_b^v(1-{\alpha }^v)},\end{array}$(3.4) ).

Solving (3.10a(3.10a) $S^{vb}=\frac{k_b^v}{1+{\eta }_b^v{N}^{vb}}{N}^{vb}\frac{b^v{N}^{vb}}{N^{vb}+b}+(1-{\lambda }^v)S^{vb}{\alpha }_1^v,$(3.10a) )–(3.10c(3.10c) $I^{vb}={\gamma }^v{E}^{vb}{\alpha }_4^v+I^{vb}{\alpha }_5^v,$(3.10c) ), we have (3.11) $\begin{array}{rl}{S}^{vb}& =\frac{k_b^v{b}^v{N^{vb}}^2}{(1+{\eta }_b^v{N}^{vb})(b+N^{vb})(1-(1-{\lambda }^v\left){\alpha }_1^v\right)},\\ E^{vb}& =\frac{{\lambda }^v{\alpha }_2^v}{1-(1-{\gamma }^v){\alpha }_3^v}{S}^{vb}={\omega }_1^v{\lambda }^v{S}^{vb},\\ I^{vb}& =\frac{{\gamma }^v{\alpha }_4^v}{1-{\alpha }_5^v}\cdot \frac{{\lambda }^v{\alpha }_2^v}{1-(1-{\gamma }^v){\alpha }_3^v}{S}^{vb}={\omega }_2^v{\lambda }^v{S}^{vb},\end{array}$(3.11) where ${\omega }_1^v$ and ${\omega }_2^v$ are defined in (2.14(2.14) ${\omega }_1^v:=\frac{{\alpha }_2^v}{1-(1-{\gamma }^v){\alpha }_3^v},{\omega }_2^v:=\frac{{\gamma }^v{\alpha }_2^v{\alpha }_4^v}{(1-{\alpha }_5^v)(1-(1-{\gamma }^v\left){\alpha }_3^v\right)}.$(2.14) ). Then we obtain (3.12) $S^{vb}=\frac{N^{vb}}{1+{\lambda }^v({\omega }_1^v+{\omega }_2^v)}.$(3.12) Substituting (3.11(3.11) $\begin{array}{rl}{S}^{vb}& =\frac{k_b^v{b}^v{N^{vb}}^2}{(1+{\eta }_b^v{N}^{vb})(b+N^{vb})(1-(1-{\lambda }^v\left){\alpha }_1^v\right)},\\ E^{vb}& =\frac{{\lambda }^v{\alpha }_2^v}{1-(1-{\gamma }^v){\alpha }_3^v}{S}^{vb}={\omega }_1^v{\lambda }^v{S}^{vb},\\ I^{vb}& =\frac{{\gamma }^v{\alpha }_4^v}{1-{\alpha }_5^v}\cdot \frac{{\lambda }^v{\alpha }_2^v}{1-(1-{\gamma }^v){\alpha }_3^v}{S}^{vb}={\omega }_2^v{\lambda }^v{S}^{vb},\end{array}$(3.11) ) and (2.15(2.15) ${\lambda }^v={\beta }^{h}r\frac{I^h}{N^h}={\beta }^{h}r\frac{{\omega }_1^h{\omega }_2^h{\lambda }^h}{1+A_2^h{\lambda }^h}$(2.15) ) into $L\left(t\right)$ in (2.5(2.5) ${\lambda }^h\left(n\right)=G\left(L\right(n\left)\right),$(2.5) ) yields $\begin{array}{rl}L\left({\lambda }^h\right)& ={\beta }^{v}r\frac{I^{vb}}{N^h}=\frac{{\beta }^{v}r{\omega }_1^v{\alpha }_2^v{S}^{vb}(1-{\alpha }_1^h+A_1^h{\lambda }^h)}{\mathrm{\Lambda }(1+A_2^h{\lambda }^h)}{\lambda }^v\\ & =\frac{{\beta }^{v}r{\omega }_1^v{\alpha }_2^v{S}^{vb}(1-{\alpha }_1^h+A_1^h{\lambda }^h)}{\mathrm{\Lambda }(1+A_2^h{\lambda }^h)}\cdot \frac{{\beta }^{h}r{\omega }_1^h{\omega }_2^h{\lambda }^h}{1+A_2^h{\lambda }^h}\\ & ={R_0^b}^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h{)}^2}\frac{S^{vb}\left({\lambda }^h\right)}{S_0^{vb}}{\lambda }^h.\end{array}$ It follows from (3.12(3.12) $S^{vb}=\frac{N^{vb}}{1+{\lambda }^v({\omega }_1^v+{\omega }_2^v)}.$(3.12) ) that $\frac{S^{vb}\left({\lambda }^h\right)}{S_0^{vb}}=\frac{1}{1+{\lambda }^v({\omega }_1^v+{\omega }_2^v)}.$ Then (3.13) $\begin{array}{rl}L\left({\lambda }^h\right)& ={R_0^b}^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h{)}^2}\frac{1}{1+{\lambda }^v({\omega }_1^v+{\omega }_2^v)}{\lambda }^h\\ & ={R_0^b}^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h)(1+D{\lambda }^h)}{\lambda }^h,\end{array}$(3.13) where $D=A_2^h+{\beta }^{h}r{\omega }_1^h{\omega }_2^h({\omega }_1^v+{\omega }_2^v)>0$ and $A_2^h>0$. Hence, there exists an endemic fixed point if and only if there is a positive solution to equation ${\lambda }^h=G\left(L\right({\lambda }^h\left)\right)$, or equivalently, (3.14) $H^b\left({\lambda }^h\right):=G\left(L\right({\lambda }^h\left)\right)-{\lambda }^h,$(3.14) for $0\le {\lambda }^h\le 1$.

Notice that $L\left(0\right)=0$ and hence $G\left(L\right(0\left)\right)=0$ by (2.6(2.6) $G:[0,\mathrm{\infty })\to [0,1],G\left(0\right)=0,G^{\prime }\left(L\right)>0,G^{″}\left(L\right)<0,\mathrm{\forall }L\ge 0,\underset{L\to \mathrm{\infty }}{lim}G\left(L\right)=1.$(2.6) ). Then $H^b\left(0\right)=0$ and $H^b\left(1\right)<0$ since $G^{\prime }\left(L\right)>0$ and $\underset{L\to \mathrm{\infty }}{lim}G\left(L\right)=1$. Moreover, it follows from $L^{\prime }\left(0\right)=R_0^{b2}/G^{\prime }\left(0\right)$ that ${H^b}^{\prime }\left(0\right)=G^{\prime }\left(0\right)L^{\prime }\left(0\right)-1=R_0^{b2}-1.$ Thus if $R_0^b>1$, there exists a positive solution ${\lambda }^{\ast }$ to $H^b\left({\lambda }^h\right)=0$, that is, to ${\lambda }^{\ast }=G\left(L\right({\lambda }^{\ast }\left)\right)$. Then there exists an endemic fixed point if $R_0^b>1$.

Next, we prove the uniqueness of this endemic fixed point.

We rewrite the equation $H^b\left({\lambda }^h\right)=G\left(L\right({\lambda }^h\left)\right)-{\lambda }^h$ as $H^b\left(x\right)=G\left(L\right(x\left)\right)-x$ for convenience and take the first derivative with respect to variable x. Then we have ${H^b}^{\prime }\left(x\right)=G^{\prime }\left(L\right(x\left)\right)L^{\prime }\left(x\right)-1,$ where $G^{\prime }\left(L\right(x\left)\right)>0$. It follows from (3.13(3.13) $\begin{array}{rl}L\left({\lambda }^h\right)& ={R_0^b}^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h{)}^2}\frac{1}{1+{\lambda }^v({\omega }_1^v+{\omega }_2^v)}{\lambda }^h\\ & ={R_0^b}^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h)(1+D{\lambda }^h)}{\lambda }^h,\end{array}$(3.13) ) that $L^{\prime }\left(x\right)=\frac{{R_0^b}^2}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)}\cdot \frac{l_1\left(x\right)}{(1+A_2^{h}x{)}^2(1+Dx{)}^2},$ where $l_1\left(x\right):=k{x}^2+2A_1^{h}x+(1-{\alpha }_1^h)$ and we write (3.15) $k:=A_1^h(A_2^h+D)-A_2^{h}D(1-{\alpha }_1^h).$(3.15) If k<0, there exists a unique positive solution $x^{\ast }$ satisfying $l_1\left(x^{\ast }\right)=0$. Then $L^{\prime }\left(x\right)>0$ if $x\in (0,x^{\ast })$ and $L^{\prime }\left(x\right)<0$ which leads to ${H^b}^{\prime }\left(x\right)<0$ if $x>x^{\ast }$. Thus $H^b\left(x\right)$ is monotone decreasing if $x>x^{\ast }$.

We then prove that $L^{\prime }\left(x\right)$ is monotone decreasing for $x<x^{\ast }$, which is equivalent to $L^{″}\left(x\right)<0$ for $x<x^{\ast }$.

It follows from (3.13(3.13) $\begin{array}{rl}L\left({\lambda }^h\right)& ={R_0^b}^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h{)}^2}\frac{1}{1+{\lambda }^v({\omega }_1^v+{\omega }_2^v)}{\lambda }^h\\ & ={R_0^b}^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h)(1+D{\lambda }^h)}{\lambda }^h,\end{array}$(3.13) ) that $L^{″}\left(x\right)=\frac{2{R_0^b}^2}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)}\cdot \frac{l_2\left(x\right)}{(1+A_2^{h}x{)}^3(1+Dx{)}^3},$ where (3.16) $l_2\left(x\right)=[-k{A}_2^{h}D]x^3-3A_1^h{A}_2^{h}D{x}^2-3A_2^{h}D(1-{\alpha }_1^h)x+A_1^h-(A_2^h+D)(1-{\alpha }_1^h).$(3.16) Since $A_1^h<A_2^{h}D(1-{\alpha }_1^h)/\left(A_2^h+D\right)$ from (3.15(3.15) $k:=A_1^h(A_2^h+D)-A_2^{h}D(1-{\alpha }_1^h).$(3.15) ) with k<0, the constant part of function $l_2\left(x\right)$ becomes $\begin{array}{rl}{A}_1^h-(A_2^h+D)(1-{\alpha }_1^h)& <\frac{A_2^{h}D(1-{\alpha }_1^h)}{A_2^h+D}-(A_2^h+D)(1-{\alpha }_1^h)\\ & <\frac{1-{\alpha }_1^h}{A_2^h+D}\cdot \left(-(A_2^h-\frac{D}{2}{)}^2-\frac{3D^2}{4}\right)<0.\end{array}$ Hence there exists at most one positive solution $x_0$ satisfying $l_2\left(x_0\right)=0$ from Descartes' rule of sign. Thus $L^{″}\left(x\right)<0$ if $x<x_0$ and $L^{″}\left(x\right)>0$ if $x>x_0$.

In order to compare $x^{\ast }$ with $x_0$, we assume $x^{\ast }>x_0$. Then, for any $x\in (x_0,x^{\ast })$, $L^{\prime }\left(x\right)$ is increasing and $L^{\prime }\left(x\right)>0$ but $L^{\prime }\left(x^{\ast }\right)=0$, which is a contradiction. Thus, $x^{\ast }<x_0$, and then $L^{″}\left(x\right)<0$ for $x<x^{\ast }$. It follows from (3.17) ${H^b}^{″}\left(x\right)=G^{″}\left(L\right(x\left)\right)\cdot \left(L^{\prime }\right(x){)}^2+G^{\prime }(L\left(x\right))\cdot L^{″}(x),$(3.17) where $G^{″}\left(L\right(x\left)\right)<0$, that ${H^b}^{″}\left(x\right)<0$ for $x<x^{\ast }$. Therefore, function $H^b\left(x\right)$ is concave down for $x<x^{\ast }$ and then decreasing for $x>x^{\ast }$ with ${H^b}^{\prime }\left(0\right)>0$, which implies that the endemic fixed point is unique.

If k>0, then $A_1^h>0$ and thus $L^{\prime }\left(x\right)>0$ for all x. If the constant part of $l_2\left(x\right)$ from (3.16(3.16) $l_2\left(x\right)=[-k{A}_2^{h}D]x^3-3A_1^h{A}_2^{h}D{x}^2-3A_2^{h}D(1-{\alpha }_1^h)x+A_1^h-(A_2^h+D)(1-{\alpha }_1^h).$(3.16) ) is negative, $l_2\left(x\right)<0$ for all $x\in (0,1]$ and thus $L^{″}\left(x\right)<0$, which implies that ${H^b}^{″}\left(x\right)<0$ from (3.17(3.17) ${H^b}^{″}\left(x\right)=G^{″}\left(L\right(x\left)\right)\cdot \left(L^{\prime }\right(x){)}^2+G^{\prime }(L\left(x\right))\cdot L^{″}(x),$(3.17) ). Since function $H^b\left(x\right)$ is concave down for all $x\in (0,1]$, the endemic fixed point is unique.

If the constant part of $l_2\left(x\right)$ from (3.16(3.16) $l_2\left(x\right)=[-k{A}_2^{h}D]x^3-3A_1^h{A}_2^{h}D{x}^2-3A_2^{h}D(1-{\alpha }_1^h)x+A_1^h-(A_2^h+D)(1-{\alpha }_1^h).$(3.16) ) is positive, there exists a unique solution $x_1$ such that $L^{″}\left(x_1\right)=0$. Then $L^{″}\left(x\right)>0$ if $x<x_1$ and $L^{″}\left(x\right)<0$ if $x>x_1$. Hence ${H^b}^{\prime }\left(x\right)$ is monotone increasing if $x<x_1$ and monotone decreasing if $x>x_1$ with ${H^b}^{\prime }\left(0\right)={R^b}^2-1>0$, where there exists a unique solution $\overline{x}$ such that ${H^b}^{\prime }\left(\overline{x}\right)=0$. Therefore, function $H^b\left(x\right)$ is monotone increasing for $x<\overline{x}$ and monotone decreasing for $x>\overline{x}$. Thus the endemic fixed point is unique.

In summary, we have shown that when $R_0^b>1$, the infection-free fixed point is unstable and there exists a unique endemic fixed point.

In the following, we prove that if $R_0^b<1$, there exists no endemic fixed point.

If k<0, it follows from the analysis above that $H^b\left(x\right)$ is concave down if $x<x^{\ast }$ and monotone decreasing if $x>x^{\ast }$. Thus there exists no intersection with x-axis since ${H^b}^{\prime }\left(0\right)={R^b}^2-1<0$.

If k>0, it follows from (3.16(3.16) $l_2\left(x\right)=[-k{A}_2^{h}D]x^3-3A_1^h{A}_2^{h}D{x}^2-3A_2^{h}D(1-{\alpha }_1^h)x+A_1^h-(A_2^h+D)(1-{\alpha }_1^h).$(3.16) ) that $H^b\left(x\right)$ is concave down with ${H^b}^{\prime }\left(0\right)<0$ if the constant part of $l_2\left(x\right)$ is negative, and hence there exists no endemic fixed point. On the other hand, if the constant part of $l_2\left(x\right)$ is positive, ${H^b}^{\prime }\left(x\right)$ is monotone increasing if $x<x_1$ and monotone decreasing if $x>x_1$ with ${H^b}^{\prime }\left(0\right)={R^b}^2-1<0$. Thus ${H^b}^{″}\left(x_1\right)=0$ and $L^{\prime }\left(x_1\right)=0$ since $L^{″}\left(x_1\right)=0$ from (3.17(3.17) ${H^b}^{″}\left(x\right)=G^{″}\left(L\right(x\left)\right)\cdot \left(L^{\prime }\right(x){)}^2+G^{\prime }(L\left(x\right))\cdot L^{″}(x),$(3.17) ). Then ${H^b}^{\prime }\left(x_1\right)=G^{\prime }\left(L\right(x_1\left)\right)\cdot L^{\prime }\left(x_1\right)-1=-1<0$, which implies that ${H^b}^{\prime }\left(x\right)<0$ for all $x\in (0,1]$. Hence $H^b\left(x\right)$ is monotone decreasing for all $x\in (0,1]$ and there exists no endemic fixed point.

4. Impact of releases of sterile mosquitoes

To study the impact of the releases of sterile mosquitoes on the malaria transmission dynamics, first, it is clear that if $b\ge \overline{b}$, there is no infection. We then consider the interval $(0,\overline{b})$, where $\overline{b}$ is given in (3.9(3.9) $\overline{b}=\frac{k_b^v{b}^v{S}_0^{vb}}{(1-{\alpha }_1^v)(1+{\eta }_b^v{S}_0^{vb})}-S_0^{vb}=\frac{k_b^v{b}^v{S}_0^v}{(1-{\alpha }_1^v)(R_0^2+{\eta }_b^v{S}_0^v)}-\frac{S_0^v}{R_0^2}.$(3.9) ). For each $b\in (0,\overline{b})$, the corresponding reproductive number $R_0^b>1$ and there exists a unique endemic fixed point associated with ${\lambda }^{\ast }$ to equation $H^b\left(\lambda \right)=0$.

With this ${\lambda }^{\ast }\left(b\right)$, we have ${\left.\frac{\mathrm{\partial }{H}^b}{\mathrm{\partial }\lambda }\right|}_{\lambda ={\lambda }^{\ast }\left(b\right)}=G^{\prime }\left({\lambda }^{\ast }\right)L^{\prime }\left({\lambda }^{\ast }\right)-1<0.$ Then it follows from (3.13(3.13) $\begin{array}{rl}L\left({\lambda }^h\right)& ={R_0^b}^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h{)}^2}\frac{1}{1+{\lambda }^v({\omega }_1^v+{\omega }_2^v)}{\lambda }^h\\ & ={R_0^b}^2\frac{1-{\alpha }_1^h+A_1^h{\lambda }^h}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^h)(1+D{\lambda }^h)}{\lambda }^h,\end{array}$(3.13) ) that $L\left(b\right)=L\left(R_0^b\right(b),{\lambda }^{\ast }(b\left)\right)$ since both $R_0^b$ and ${\lambda }^{\ast }$ can be regarded as functions of variable b. Thus (4.1) $\begin{array}{rl}{L}^{\prime }\left(b\right)& =\frac{\mathrm{\partial }L\left(b\right)}{\mathrm{\partial }{R}_0^b\left(b\right)}\cdot {R_0^b}^{\prime }\left(b\right)+\frac{\mathrm{\partial }L\left(b\right)}{\mathrm{\partial }{\lambda }^{\ast }\left(b\right)}\cdot {{\lambda }^{\ast }}^{\prime }\left(b\right)\\ & =\frac{2R_0^b{\lambda }^{\ast }\left(b\right)(1-{\alpha }_1^h+A_1^h{\lambda }^{\ast }(b\left)\right)}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^{\ast }(b\left)\right)(1+D{\lambda }^{\ast }(b\left)\right)}\cdot {R_0^b}^{\prime }\left(b\right)+L^{\prime }\left({\lambda }^{\ast }\right(b\left)\right)\cdot {{\lambda }^{\ast }}^{\prime }\left(b\right),\end{array}$(4.1) and by taking the derivative of equation $H^b\left(\lambda \right(b\left)\right)=0$ with respect to b on both sides, $\begin{array}{rl}{H^b}^{\prime }\left(b\right)& =G^{\prime }\left(L\right(b\left)\right)L^{\prime }\left(b\right)-{{\lambda }^{\ast }}^{\prime }\left(b\right)\\ & =\frac{2R_0^b{G}^{\prime }\left({\lambda }^{\ast }\right){\lambda }^{\ast }\left(b\right)(1-{\alpha }_1^h+A_1^h{\lambda }^{\ast }(b\left)\right)}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^{\ast }(b\left)\right)(1+D{\lambda }^{\ast }(b\left)\right)}\\ & \cdot {R_0^b}^{\prime }\left(b\right)+\left(G^{\prime }\right({\lambda }^{\ast }\left)L^{\prime }\right({\lambda }^{\ast })-1)\cdot {{\lambda }^{\ast }}^{\prime }\left(b\right)\\ & =0.\end{array}$ Equivalently, (4.2) ${{\lambda }^{\ast }}^{\prime }\left(b\right)\left(G^{\prime }\right({\lambda }^{\ast }\left)L^{\prime }\right({\lambda }^{\ast })-1)=-\frac{2R_0^b{G}^{\prime }\left({\lambda }^{\ast }\right){\lambda }^{\ast }\left(b\right)(1-{\alpha }_1^h+A_1^h{\lambda }^{\ast }(b\left)\right)}{G^{\prime }\left(0\right)(1-{\alpha }_1^h)(1+A_2^h{\lambda }^{\ast }(b\left)\right)(1+D{\lambda }^{\ast }(b\left)\right)}\cdot {R_0^b}^{\prime }\left(b\right),$(4.2) where $G^{\prime }\left({\lambda }^{\ast }\right(b\left)\right)L^{\prime }\left({\lambda }^{\ast }\right(b\left)\right)-1<0$ and $R_0^{b^{\prime }}<0$. Hence ${{\lambda }^{\ast }}^{\prime }\left(b\right)<0$, which implies that ${\lambda }^{\ast }$ is monotone decreasing with respect to b.

Moreover, it follows from (2.11(2.11) $\begin{array}{rl}{E}^h& ={\omega }_1^h{\lambda }^h{S}^h,I^h={\omega }_1^h{\omega }_2^h{\lambda }^h{S}^h,\\ R^h& ={\omega }_1^h{\omega }_2^h{\omega }_3^h{\lambda }^h{S}^h,S^h=\frac{\mathrm{\Lambda }}{1-{\alpha }_1^h+A_1^h{\lambda }^h},\end{array}$(2.11) ) that $I^h\left(b\right)=\mathrm{\Lambda }{\omega }_1^h{\omega }_2^h\frac{{\lambda }^{\ast }\left(b\right)}{1-{\alpha }_1^h+A_1^h{\lambda }^{\ast }\left(b\right)},$ and then $I^{h^{\prime }}\left(b\right)=\frac{{\omega }_1^h{\omega }_2^h\mathrm{\Lambda }(1-{\alpha }_1^h)}{(1-{\alpha }_1^h+A_1^h{\lambda }^{\ast }{)}^2}{\lambda }^{{\ast }^{\prime }}<0.$ Thus component $I^h\left(b\right)$ of the unique endemic fixed point is monotone decreasing with respect to b. In other words, we can reduce the population of the infected humans by increasing the amount of sterile mosquitoes.

Furthermore, it follows from (3.11(3.11) $\begin{array}{rl}{S}^{vb}& =\frac{k_b^v{b}^v{N^{vb}}^2}{(1+{\eta }_b^v{N}^{vb})(b+N^{vb})(1-(1-{\lambda }^v\left){\alpha }_1^v\right)},\\ E^{vb}& =\frac{{\lambda }^v{\alpha }_2^v}{1-(1-{\gamma }^v){\alpha }_3^v}{S}^{vb}={\omega }_1^v{\lambda }^v{S}^{vb},\\ I^{vb}& =\frac{{\gamma }^v{\alpha }_4^v}{1-{\alpha }_5^v}\cdot \frac{{\lambda }^v{\alpha }_2^v}{1-(1-{\gamma }^v){\alpha }_3^v}{S}^{vb}={\omega }_2^v{\lambda }^v{S}^{vb},\end{array}$(3.11) ) and (3.12(3.12) $S^{vb}=\frac{N^{vb}}{1+{\lambda }^v({\omega }_1^v+{\omega }_2^v)}.$(3.12) ) that $I^v\left(b\right)={\omega }_2^v\cdot \frac{{\lambda }^v\left(b\right)N^{vb+}\left(b\right)}{1+{\lambda }^v\left(b\right)({\omega }_1^v+{\omega }_2^v)}.$ Taking the derivative with respect to b, we obtain $I^{v^{\prime }}\left(b\right)=\frac{{\omega }_2^v}{(1+{\lambda }^v(b\left)\right({\omega }_1^v+{\omega }_2^v){)}^2}\cdot \left(N^{vb+}{\lambda }^{v^{\prime }}+(1+({\omega }_1^v+{\omega }_2^v\left){\lambda }^v\right)N^{vb{+}^{\prime }}\left(b\right)\right).$ In addition, it follows from (2.15(2.15) ${\lambda }^v={\beta }^{h}r\frac{I^h}{N^h}={\beta }^{h}r\frac{{\omega }_1^h{\omega }_2^h{\lambda }^h}{1+A_2^h{\lambda }^h}$(2.15) ) that ${\lambda }^{v^{\prime }}\left(b\right)=\frac{{\beta }^{h}r{\omega }_1^h{\omega }_2^h{\lambda }^{{\ast }^{\prime }}}{(1+A_2^h{\lambda }^{\ast }{)}^2}<0,$ and $N^{vb{+}^{\prime }}\left(b\right)<0$ from (3.8(3.8) ${S_0^{vb}}^{\prime }\left(b\right)=-\frac{(1-{\alpha }_1^v)(1+{\eta }_b^v{S}_0^{vb})}{{\mathrm{\partial }}_{S_0^{vb}}\mathrm{\Phi }\left(S_0^{vb}\right(b),b)}<0.$(3.8) ). Consequently, $I^{v^{\prime }}\left(b\right)<0$, which implies that component $I^v\left(b\right)$ of the unique endemic fixed point is also monotone decreasing with respect to b. As a result of increasing the releases of the sterile mosquitoes, we can possibly wipe out all of the wild mosquitoes eventually.

Therefore, for $b\in (0,\overline{b})$, even though $R_0^b\left(b\right)>1$ such that the disease spreads and goes to a positive steady state for $n\to \mathrm{\infty }$, we can still reduce the components of the infected humans and mosquitoes by increasing the releases of the sterile mosquitoes to make the transmission under control.

We provide an example below to demonstrate our findings.

Example 4.1

We use the following parameters for the transmission model: (4.3) $\begin{array}{cccccccc}& \mathrm{\Lambda }=6;& k_b^v=0.5;& {\eta }_b^v=1.5;& b^v=80;& {\gamma }^v=0.8;& {\beta }^v=0.1;& {\beta }^h=0.15;\\ & {\alpha }_1^h=0.3;& {\alpha }_2^h=0.7;& {\alpha }_3^h=0.5;& {\alpha }_4^h=0.8;& {\alpha }_5^h=0.8;& {\alpha }_6^h=0.8;& {\alpha }_7^h=0.8;\\ & {\alpha }_8^h=0.8;& {\gamma }^h=0.7;& {\eta }^h=0.25;& {\theta }^h=0.5;& r=6;\end{array}$(4.3) and ${\alpha }_1^v={\alpha }_2^v={\alpha }_3^v={\alpha }_4^v={\alpha }_5^v=0.5$.

Before the sterile mosquitoes are released, the reproductive number of infection for system (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ) is $R_0=1.3039>1$ and hence the disease spreads. After the releases of sterile mosquitoes, we have the existence threshold value $b_c=42.0743$ such that if b=45>42.0743, there exists no positive fixed point and thus the infection eventually dies out, as shown in Figure 3. If b<42.0743, there exist two positive fixed points with components $N^{vb(\pm )}\left(b\right)=26.3333-0.5000b\pm 0.6667\sqrt{(39.5-0.75b{)}^2-1.500b},$ where $N^{vb(-)}$ is unstable and $N^{vb(+)}$ is locally asymptotically stable with initial value larger than $N^{vb(-)}$. $N^{vb(+)}$ is a decreasing function with respect to release value b.

Figure 3. With the parameters given in (4.3(4.3) $\begin{array}{cccccccc}& \mathrm{\Lambda }=6;& k_b^v=0.5;& {\eta }_b^v=1.5;& b^v=80;& {\gamma }^v=0.8;& {\beta }^v=0.1;& {\beta }^h=0.15;\\ & {\alpha }_1^h=0.3;& {\alpha }_2^h=0.7;& {\alpha }_3^h=0.5;& {\alpha }_4^h=0.8;& {\alpha }_5^h=0.8;& {\alpha }_6^h=0.8;& {\alpha }_7^h=0.8;\\ & {\alpha }_8^h=0.8;& {\gamma }^h=0.7;& {\eta }^h=0.25;& {\theta }^h=0.5;& r=6;\end{array}$(4.3) ), the threshold value is $b_c=42.0743$. If $b=45>b_c$, the trivial fixed point is globally asymptotically stable, which leads to the infection eventually extinct as shown above.

However, the threshold value for the disease spread is $\overline{b}=21.2330$ such that the reproductive number $R_0^b\left(b\right)<1$ if b=30>21.2330 and thus the disease dies out eventually, which is shown in Figure 4 (right). If b<21.2330, the reproductive number $R_0^b\left(b\right)>1$, which makes the infection-free fixed point unstable and thus the disease spreads. Moreover, $R_0^b\left(b\right)$ is monotone decreasing with respect to b, as shown in Figure 4 (left).

Figure 4. With the same parameters in (4.3(4.3) $\begin{array}{cccccccc}& \mathrm{\Lambda }=6;& k_b^v=0.5;& {\eta }_b^v=1.5;& b^v=80;& {\gamma }^v=0.8;& {\beta }^v=0.1;& {\beta }^h=0.15;\\ & {\alpha }_1^h=0.3;& {\alpha }_2^h=0.7;& {\alpha }_3^h=0.5;& {\alpha }_4^h=0.8;& {\alpha }_5^h=0.8;& {\alpha }_6^h=0.8;& {\alpha }_7^h=0.8;\\ & {\alpha }_8^h=0.8;& {\gamma }^h=0.7;& {\eta }^h=0.25;& {\theta }^h=0.5;& r=6;\end{array}$(4.3) ), the threshold values are $\overline{b}=21.2330$ and $b_c=42.0743$, respectively. By using b as a variable, the horizontal axis is for b and the vertical axis is for $R_0^b\left(b\right)$ (left). The reproductive number $R_0^b\left(0\right)=R_0=1.3039$ at b=0. For b=30, the infection-free fixed point is stable since $R_0^b\left(b\right)<1$. Then the infected humans and infected mosquitoes eventually go extinct (right).

For b<21.2330, $R_0^b\left(b\right)>1$, the corresponding ${\lambda }^h\left(b\right)$ determined in (3.14(3.14) $H^b\left({\lambda }^h\right):=G\left(L\right({\lambda }^h\left)\right)-{\lambda }^h,$(3.14) ) is a positive and decreasing function as shown in Figure 5 (left). Corresponding to ${\lambda }^h\left(b\right)$, there exists a unique endemic fixed point whose components $I^v\left(b\right)$ and $I^h\left(b\right)$ are monotone decreasing with respect to b as shown in Figure 5 (right), which indicates that increasing the releases of sterile mosquitoes reduces the disease spread.

Figure 5. With still the same parameters in (4.3(4.3) $\begin{array}{cccccccc}& \mathrm{\Lambda }=6;& k_b^v=0.5;& {\eta }_b^v=1.5;& b^v=80;& {\gamma }^v=0.8;& {\beta }^v=0.1;& {\beta }^h=0.15;\\ & {\alpha }_1^h=0.3;& {\alpha }_2^h=0.7;& {\alpha }_3^h=0.5;& {\alpha }_4^h=0.8;& {\alpha }_5^h=0.8;& {\alpha }_6^h=0.8;& {\alpha }_7^h=0.8;\\ & {\alpha }_8^h=0.8;& {\gamma }^h=0.7;& {\eta }^h=0.25;& {\theta }^h=0.5;& r=6;\end{array}$(4.3) ), the threshold values are $\overline{b}=21.2330$ and $b_c=42.0743$, respectively. The curve on the left figure is for ${\lambda }^h\left(b\right)$ at the endemic fixed point for each b. The upper curve and lower curve on the right figure are for $I^h\left(b\right)$ and $I^v\left(b\right)$ at the endemic fixed point for each b, respectively. Clearly, ${\lambda }^h\left(b\right)$, $I^h\left(b\right)$ and $I^v\left(b\right)$ are all negative if $b>\overline{b}$, which implies that no endemic fixed point exists although positive $N^{vb(+)}$ exists for $\overline{b}<b<b_c$.

To show the dynamical impact of the releases of the sterile mosquitoes, we also present the solutions of the disease transmission system in Figure 6. When no sterile mosquitoes are released, the reproductive number $R_0=1.3039>1$ and the disease spreads as shown in Figure 6 (left). For $b=15<\overline{b}=21.2330$, the reproductive number $R_0^b=1.0988>1$ such that the infection still exists but the number of infected humans and infected mosquitoes are reduced significantly as shown in Figure 6 (right).

Figure 6. (Left) With again the parameters given in (4.3(4.3) $\begin{array}{cccccccc}& \mathrm{\Lambda }=6;& k_b^v=0.5;& {\eta }_b^v=1.5;& b^v=80;& {\gamma }^v=0.8;& {\beta }^v=0.1;& {\beta }^h=0.15;\\ & {\alpha }_1^h=0.3;& {\alpha }_2^h=0.7;& {\alpha }_3^h=0.5;& {\alpha }_4^h=0.8;& {\alpha }_5^h=0.8;& {\alpha }_6^h=0.8;& {\alpha }_7^h=0.8;\\ & {\alpha }_8^h=0.8;& {\gamma }^h=0.7;& {\eta }^h=0.25;& {\theta }^h=0.5;& r=6;\end{array}$(4.3) ), the reproductive number is $R_0=1.3039>1$ and hence the disease spreads in the absence of sterile mosquitoes. (Right) After the sterile mosquitoes are released, even the released amount is less than the threshold with $b=15<\overline{b}$, the infected components $I_h$ and $I_v$ are reduced so that we can have the infection under control.

5. Concluding remarks

To have a better understanding of the impact of releasing sterile mosquitoes on malaria transmission, we first formulated a simple compartmental SEIR model for the malaria transmission, the survivability of the susceptible newborns has the Beverton–Holt-II type of nonlinearity, in (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ) as our baseline model in Section 2. We derived a formula for the reproductive number of infection $R_0$ in (2.9(2.9) $R_0:=\sqrt{\frac{G^{\prime }\left(0\right)r{\beta }^h{\gamma }^h{\alpha }_2^h{\alpha }_4^h}{(1-(1-{\gamma }^h\left){\alpha }_3^h\right)(1-(1-{\eta }^h\left){\alpha }_5^h\right)}\cdot \frac{r{\beta }^v{\gamma }^v{\alpha }_2^v{\alpha }_4^v}{(1-{\alpha }_5^v)(1-(1-{\gamma }^v\left){\alpha }_3^v\right)}\cdot \frac{S_0^v}{S_0^h}},$(2.9) ) and showed that the infection-free fixed point of the baseline model is asymptotically stable if $R_0<1$ and unstable if $R_0>1$. We also showed that if $R_0>1$, there exists an endemic fixed point for the baseline model.

We then introduced sterile mosquitoes into our baseline model and formulated the interactive model in Section 3. The human components are the same as those in (2.1(2.1) $\begin{array}{rl}{S}^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\\ S^v(n+1)& =b^v{N}^v\left(n\right){\alpha }_b^v\left(N^v\right(n\left)\right)+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\end{array}$(2.1) ) but the mosquito components are given in (3.1(3.1) $\begin{array}{rl}{S}^v(n+1)& =\frac{k_b^v}{1+{\eta }_b^v{N}^v\left(n\right)}{N}^v\left(n\right)\frac{b^v{N}^v\left(n\right)}{N^v\left(n\right)+b}+(1-{\lambda }^v(n\left)\right)S^v\left(n\right){\alpha }_1^v,\\ E^v(n+1)& ={\lambda }^v\left(n\right)S^v\left(n\right){\alpha }_2^v+(1-{\gamma }^v)E^v\left(n\right){\alpha }_3^v,\\ I^v(n+1)& ={\gamma }^v{E}^v\left(n\right){\alpha }_4^v+I^v\left(n\right){\alpha }_5^v,\\ S^h(n+1)& =\mathrm{\Lambda }+(1-{\lambda }^h(n\left)\right)S^h\left(n\right){\alpha }_1^h+{\theta }^h{R}^h\left(n\right){\alpha }_8^h,\\ E^h(n+1)& ={\lambda }^h\left(n\right)S^h\left(n\right){\alpha }_2^h+(1-{\gamma }^h)E^h\left(n\right){\alpha }_3^h,\\ I^h(n+1)& ={\gamma }^h{E}^h\left(n\right){\alpha }_4^h+(1-{\eta }^h)I^h\left(n\right){\alpha }_5^h,\\ R^h(n+1)& ={\eta }^h{I}^h\left(n\right){\alpha }_6^h+(1-{\theta }^h)R^h\left(n\right){\alpha }_7^h,\end{array}$(3.1) ). We only consider the case of constant releases of sterile mosquitoes. We derived a formula for the reproductive number $R_0^b$, presented in (3.7(3.7) $R_0^b\left(b\right):=\sqrt{\frac{G^{\prime }\left(0\right)r{\beta }^h{\gamma }^h{\alpha }_2^h{\alpha }_4^h}{(1-(1-{\gamma }^h\left){\alpha }_3^h\right)(1-(1-{\eta }^h\left){\alpha }_5^h\right)}\cdot \frac{r{\beta }^v{\gamma }^v{\alpha }_2^v{\alpha }_4^v}{(1-{\alpha }_5^v)(1-(1-{\gamma }^v\left){\alpha }_3^v\right)}\cdot \frac{S_0^{vb(+)}\left(b\right)}{S_0^h}},$(3.7) ), and showed $R_0^b<R_0$ for any release amount b. We proved that the infection-free fixed point is asymptotically stable if $R_0^b<1$ and is unstable if $R_0^b>1$. Using the rate of releases b as a variable, we determined threshold value $b_c$ for the existence of a unique endemic fixed point for the interactive wild and sterile mosquitoes and threshold value $\overline{b}$ that determines whether $R_0^b<1$ or $R_0^b>1$, that is, whether the disease dies out or spreads.

We studied the impact of the releases of sterile mosquitoes on the transmission dynamics in Section 4 by investigating the variation of the reproductive number, $R_0^b\left(b\right)$, and the infected components, $I^h\left(b\right)$ and $I^v\left(b\right)$, induced by ${\lambda }^h\left(b\right)$, as b varies, respectively. Based on the threshold values $b_c$ and $\overline{b}$, we provided Example 4.1 to confirm and demonstrate our findings. If the rate of releases b is greater than threshold $b_c$, the trivial fixed point for the interactive wild and sterile mosquitoes is locally asymptotically stable. All wild mosquitoes are wiped out and thus there is no infection. On the other hand, if $b<b_c$, while the wild mosquitoes cannot be wiped out and the two types of mosquitoes coexist, the disease can still go extinct when there are sufficient sterile mosquitoes with $\overline{b}<b<b_c$, which leads to $R_0^b<1$. In the case even if we are unable to release enough sterile mosquitoes with $b<\overline{b}$, which leads to $R_0^b>1$ and thus the disease spreads when the initial wild mosquitoes is larger than $N^{vb(-)}$ given in (3.4(3.4) $\begin{array}{rl}& N^{vb(\pm )}\\ & =\frac{k_b^v{b}^v-(1+b{\eta }_b^v)(1-{\alpha }^v)\pm \sqrt{(k_b^v{b}^v-(1+b{\eta }_b^v\left)\right(1-{\alpha }^v){)}^2-4b{\eta }_b^v(1-{\alpha }^v{)}^2}}{2{\eta }_b^v(1-{\alpha }^v)},\end{array}$(3.4) ), the infected components $I^h\left(b\right)$ and $I^v\left(b\right)$ are monotone decreasing with respect to b; that is, the infection can be reduced. All of the results are demonstrated by numerical simulations in Example 4.1.

In conclusion, based on the discrete-time interactive malaria model we formulated and the analysis we accomplished in this study, we showed that the releases of sterile mosquitoes have significant impacts on controlling the transmission of malaria. It will be useful and helpful in considering effective strategies for control and prevention of mosquito-borne diseases as well.

Disclosure statement

No potential conflict of interest was reported by the authors.