Study of the sterile insect release technique for a two-sex mosquito population model

1.   Introduction

As an insect control method, the sterile insect technique (SIT) has made a big success in pest management and disease vector control over the past few decades ^[1,2,3,4,5]. In most instances, sterile males are released into an area to compete with wild males. If a female mates with a sterile one, it is unable to produce offspring, then the insect population will be eventually controlled or even eliminated by the reproductive mate attrition. SIT was verified to work against mosquitoes by field trials in the 1970s and 1980s ^[6], and then developed rapidly for emerging approaches from SIT and new rearing techniques ^[7,8,9,10]. In addition to the classic SIT, the biological control methods already used in the laboratory or field also include the genetic approaches and the Wolbachia driven mosquito control technique.

To investigate the impact and effectiveness of these biological control measures, various mathematical models have been formulated and analyzed. K. R. Fister et al. ^[11] proposed an optimal control framework to explore the effect of sterile mosquito releases on reducing the incidence of mosquito-borne diseases. S. M. White et al. explored the mechanism how an insect fitness cost affects different control policies by constructing a stage-structured mathematical model for the mosquito Aedes aegypti in ^[12]. J. Li et al. ^[13,14,15] investigated multiple policies of sterile mosquito release by formulating discrete and continuous dynamical systems for the interaction between two mosquito populations. While L. Cai et al. ^[16] studied the impact of the SIT on disease transmission by constructing dynamical systems incorporating constant, proportional and Holling-Ⅱ type release rates. Y. Dumont and J. M. Tchuenche proposed mathematical models of SIT to exploit the control of an epidemic of Chikungunya in ^[17], and pulsed periodic releases was especially studied. While in ^[8], M. Strugarek et al. investigated the application of SIT in reducing and eliminating wild mosquitoes by a simplified population model for Aedes, and several release modes were considered and necessary conditions to guarantee elimination in each case were obtained. There are also lots of works focusing on the dynamical analysis of mosquito population models with different characteristics ^[18,19,20].

It was pointed out in ^[21] that the duration of the release process each time is relatively short and it usually takes multiple releases to make the mosquito population under control. Therefore, multiple pulsed releases may be a realistic assumption. The impulsive release has been studied in several cases ^{[8,17,21,22,23,24,25,26,27]}. The authors in ^[17] explained the pulsed release of sterile males by a mosquito-human epidemiological model for Chikungunya, and investigated the impact of periodic pulsed releases on the disease transmission. In ^[21], different two-dimensional models with periodic impulsive releases and state feedback impulsive releases are constructed, and the authors exploited the influence of different release strategies on the population development of wild mosquitoes. And P. A. Bliman et al.^[22] also investigated impulsive release of sterile male mosquitoes, and they studied the periodic impulsive releases under open-loop control, state feedback impulsive releases under closed-loop control and a mixed release strategy that combines open-loop and closed-loop controls. While research in ^[8] established necessary conditions which can guarantee the eventual extinction of the wild mosquito population.

For the Wolbachia driven mosquito control technique, J. Yu ^[23] introduced a model of differential equations with a time delay to study the suppression dynamics of wild mosquitoes intervened by the releasing of Wolbachia-infected males. Unlike many studies, the population suppression in this work tried to avoid releasing infected females, just released living infected males, and aimed for eliminating the whole population of mosquitoes. Then J. Yu et al. introduced the sexual lifespan of sterile mosquitoes and assumed that the interaction happens only when the sterile mosquitoes are still sexually active ^{[24,25,26,27]}. They investigated the impact of the sexual lifespan of sterile mosquitoes on mosquito population suppression based on delay differential equations and gave a lot of important results. J. Yu and B. Zheng ^[28] specially studied Wolbachia persistence by extra releases of Wolbachia-infected mosquitoes based on difference equations, and obtained a maximal maternal leakage rate threshold such that infected mosquitoes can persist. While M. Huang et al. ^[29] introduced a system of delay differential equations, including both the adult and larval stages of wild mosquitoes, interfered by Wolbachia infected males. They explored its global dynamics and determined a threshold level of infected male releasing which can ensure that the wild population is suppressed completely.

Although rearing techniques are steadily updated, mosquito mass rearing is still one of the major obstacles preventing the application of the SIT against mosquitoes in large scale. Artificial rearing of sterile mosquitoes comes at an economic cost which cannot be neglected. Optimal control method can be used to balance the conflict between the control level of the wild mosquitoes and the economic cost.

Optimal control problem for impulsive dynamical systems has its special peculiarities, and there is a technical difficulty compared to continuous dynamical systems due to the dependence of the state of variables on uncertain pulse effects. Many researchers have been making contributions to overcome the difficulty and providing available control methodologies. These methodologies have been used in the optimal impulsive management of population ^{[11,13,30,31,32,33,34]}. For example, in ^[33] the authors provided multiple kinds of optimal policies for an eco-epidemiological model with impulsive interferences. An pest management system incorporating impulsive release natural enemies is studied in ^[34], and three optimal release strategies are given.

In this study, we did not consider the effect of the sexual lifespan and propose a two-sex mosquito population model with stage structure and impulsive releases of sterile males, and an Allee effect is also incorporated to describe the scarcity of the available mating area in the field. We firstly study the large-scale time control based on aims of the extinction of wild mosquitoes, and then investigate limited-time optimal control of wild mosquitoes to gain a suitable control strategy by selecting optimal release parameters.

The structure of paper is as follows: In Section 2, we establish a hybrid dynamical system for the large-scale time control of wild mosquitoes with impulsive releases of sterile males, and then study its dynamical properties and exploit threshold conditions whether or not the wild mosquito population is eliminated. In Section 3, we take into account the control effect of mosquito population level and the economic input, and raise three limited-time optimal control problems for the impulsive release strategies. By using a time rescaling technique, the gradients of cost function with respect to all control parameters are obtained. Then numerical simulations are performed in Section 4 to determine the optimal values of the release timing and release amount. Finally, a brief conclusion is presented in Section 5.

2.   Release strategies to eliminate wild mosquitos for large-scale time control

2.1.   Model formulation

The model proposed in this section aims to investigate release tactics of sterilizing males which can effectively reduce and eventually eliminate wild mosquitoes in a field.

According to the life habits of many species of mosquitoes (for example, Aedes genus), eggs produced by fertile females may keep unhatched for a long time to wait for rainy seasons when the natural breeding sites are available. So there may be large egg stocks in a given field which has to be taken into account when built our mathematic model. Besides, in many cases the available mating area in the field is relatively scarce and fertile females need overcome difficulties to get successful fertilization. Therefore, we add an Allee effect in the wild female population as some researchers have done. Based on the model assumption in ^[8], we propose the following population development model for wild mosquitoes with impulsive releases of sterile males

with $K > W_E(0) > 0, W_M(0) > 0, W_F(0) > 0$ and $G_M(0) > 0$. $W_E(t), W_M(t)$ and $W_F(t)$ represent the densities of eggs, fertile males and fertile females of wild mosquitoes at time $t$, while $G_M(t)$ is the density of sterilizing males released in the field. The logistic term $\beta W_F (1-\frac{W_E}{K})$ describes the "skip oviposition" behavior of fertile females for they are capable to avoid depositing eggs in an area which has supported too many larvae. $\beta$ measures the effective fecundity and $K$ is the environmental carrying capacity. $\mu_i (i = 1, 2, 3, 4)$ denote the death rates, $\rho$ is the hatching rate and $\theta$ represents the sex ratio of wild mosquitoes. $b$ stands for the insemination rate of emerging females, while parameters $\gamma$ and $\alpha$ measure the strength of Allee effect which involves a female's mating likelihood and the mating competitiveness of sterile males, respectively. Besides, $\omega$ is the release period and $\delta$ is the amount released each time.

In ^[8], the Allee effect caused by mating limitation is modeled by a negative exponential function $1-\exp(-k W_M)$, while in this work a rectangular hyperbola function $\frac{b W_M}{\gamma +W_M+\alpha G_M}$ is used. This form of Allee effect is also widely used to describe mating limitation in sexually reproducing organisms ^{[14,15,16,35,36]}. We hope to explore this kind of Allee effect in this work and to see if it will cause any changes.

In the next subsection, we will investigate dynamical behaviors of system (2.1), and exploit the release tactics for the large-scale time control of wild mosquitoes aiming to wipe out wild mosquitoes from the field eventually.

2.2.   Existence and stability of wild mosquito-free periodic solution

We firstly discuss the existence of the wild mosquito-free periodic solution of system (2.1), then determine conditions of its global stability. We will provide theoretical analysis and practical method for the selection of the release amount $\delta$ and the release period $\omega$ so that the wild mosquitoes can be wiped out from the field.

Let $X(t) = (W_E(t), W_M(t), W_F(t), G_M(t))^T$ be any solution of system (2.1). Obviously, $X(t)$ is piecewise continuous and $X(k\omega^+) = \lim \limits_{\varepsilon\rightarrow 0^+}X(k\omega+\varepsilon)$ exists. Since the right hand side of system (2.1) is locally lipschitz continuous on $R_+^4$, system (2.1) has a unique solution ^[37,38].

The positivity and boundedness of the solution of system (2.1) are firstly investigated.

Proposition 1. Solutions of system (2.1) are always non-negative if the initial values are non-negative.

Proof. Let $(W_E(t), W_M(t), W_F(t), G_M(t))$ be a solution of system (2.1) with non-negative initial conditions. From the forth equation of system (2.1), we can get ${d G_M(t)}/{dt} = 0$ if $G_M(t) = 0$, so we have $G_M(t) > 0$, $t\geq 0$ for $G_M(0) > 0$.

Assume that there exists $t > 0$ satisfying $W_E(t) < 0$. Denote $\tau_1 = \inf\{t:W_E(t) < 0\}$. From this, we have $W_E(\tau_1) = 0$ and $W_E'(\tau_1)\leq0$. Substitute them into the first equation of system (2.1), then we obtain $W_E'(\tau_1) = \beta W_F(\tau_1)\leq0$.

If $W_F(t)$ is non-negative for all $t > 0$, then $W_F(\tau_1)\geq0$ and we get $W_E'(\tau_1)\geq 0$. It follows from $W_E'(\tau_1)\leq 0$ that $W_E'(\tau_1) = 0$, then we have $W_E'(\tau_1) = W_E(\tau_1) = W_F(\tau_1) = W_F'(\tau_1) = 0$. According to the definition of $\tau_1$, there must be a sufficiently small constant $\epsilon_1 > 0$ such that $W_E(\tau_1+\epsilon_1) < 0$ and $W_E'(\tau_1+\epsilon_1) < 0$. Since $W_F(\tau_1+\epsilon_1)\geq0$, we have $W_E'(\tau_1+\epsilon_1) = \beta W_F(\tau_1+\epsilon_1) (1-\frac{W_E(\tau_1+\epsilon_1)}{K})-(\rho +\mu_1)W_E(\tau_1+\epsilon_1) > 0$, which leads to a contradiction. Therefore, there must exist $t > 0$ such that $W_F(t) < 0$. Denote $\tau_2 = \inf\{t:W_F(t) < 0\}$, then $W_F(\tau_2) = 0$ and $W_F'(\tau_2)\leq0$. We can easily obtain $\tau_2 < \tau_1$ and $W_E(\tau_2) > 0$.

In fact, if $\tau_1 = \tau_2$, there is $W_E'(\tau_1) = \beta W_F(\tau_1) = \beta W_F(\tau_2) = 0$, and as discussed above we can get a contradiction. If $\tau_1 < \tau_2$, there is $W_F(\tau_1) > 0$, and then we have $W_E'(\tau_1) = \beta W_F(\tau_1) > 0$, which also leads to a contradiction.

Further more, if $W_M(t)$ is non-negative for all $t > 0$, then we obtain $W_M(\tau_2)\geq0$. By the third equation of system (2.1), we have

It follows from $W_F'(\tau_2)\leq 0$ that $W_F'(\tau_2) = 0$, then we have $W_F'(\tau_2) = W_F(\tau_2) = W_M(\tau_2) = 0, W_M'(\tau_2) = \theta \rho W_E(\tau_2) > 0$. According to the definition of $\tau_2$, there must be a sufficiently small constant $\epsilon_2 > 0$ such that $W_F(\tau_2+\epsilon_2) < 0$, $W_F'(\tau_2+\epsilon_2) < 0$ and $W_M(\tau_2+\epsilon_2) > 0$. Since

which leads to a contradiction. Thus there must exist some $t > 0$ such that $W_M(t) < 0$. Denote $\tau_3 = \inf\{t:W_M(t) < 0\}$, then we have $W_M(\tau_3) = 0$ and $W_M'(\tau_3)\leq0$. According to the second equation of (2.1) yields $W_M'(\tau_3) = \theta \rho W_E(\tau_3)\leq0$. If $\tau_1 > \tau_3$, then there is $W_E(\tau_3) > 0$ and $W_M'(\tau_3) = \theta \rho W_E(\tau_3) > 0$, which also leads to a contradiction. Thus we have $\tau_1\leq\tau_3$.

By the above discussion, we have $\tau_2 < \tau_1\leq\tau_3$ and $W_E(\tau_2) > 0, W_M(\tau_2) > 0, W_F(\tau_2) = 0$. However,

which also contradicts the definition of $\tau_2$.

To sum up, we have $W_E(t)$ is non-negative for all $t > 0$. Then there are ${d W_M(t)}/{dt} = 0$ if $W_M(t) = 0$ and ${d W_F(t)}/{dt}\geq0$ if $W_F(t) = 0$, and all solutions of system (2.1) with non-negative initial conditions are always non-negative. The proof is completed.

Denote

Proposition 2. $\Omega$ is a forward invariant and globally attracting set of system (2.1) in $R^4_+$.

Proof. Suppose $X(t) = (W_E(t), W_M(t), W_F(t), G_M(t))^T$ is a solution of system (2.1) with $X(0) = (W_E(0), W_M(0), W_F(0), G_M(0))^T \in R^4_+$. Firstly, we prove that $\Omega$ is a forward invariant set. By the forth equation and the impulsive conditions of system (2.1), we get

It is obvious that the dynamics of sterilizing males is completely unaffected by that of the wild mosquitoes. System (2.3) admits a unique positive periodic solution

which is globally asymptotically stable. Besides, we can easily get

and

If $X(0)\in \Omega$, then we have $G_M(0)\leq L_3 = \tilde{G}_M(0)$ and $G_M(t)\leq\tilde{G}_M(t)\leq\tilde{G}_M(0) = L_3, t\geq0$.

Besides, from the first equation of system (2.1), we get $\frac{dW_E}{dt}|_{X(t)\in \Omega}\leq \beta L_2(1-\frac{W_E}{K})$, thus we can easily obtain that $0\leq W_E(t)\leq K, t\geq0$ for any $X(0)\in \Omega$. According to the second equation of system (2.1), we have $\frac{dW_M}{dt}|_{X(t)\in \Omega}\leq \theta \rho K-\mu_2W_M = \theta \rho K(1-\frac{W_M}{L_2})$, thus there is $0\leq W_M(t)\leq L_2, t\geq0$ for any $X(0)\in \Omega$. Similarly, by the third equation of system (2.1), we have $\frac{dW_F}{dt}|_{X(t)\in \Omega}\leq (1-\theta)\rho K b-\mu_3W_F = (1-\theta)\rho K b(1-\frac{W_M}{L_3})$, and then there is $0\leq W_F(t)\leq L_3, t\geq0$ for any $X(0)\in \Omega$. Therefore, for any $X(0) = (W_E(0), W_M(0), W_F(0), G_M(0))^T \in \Omega$, there is

That is to say, $\Omega$ is a forward invariant set of system (2.1) in $R^4_+$.

In the following, we prove that $\Omega$ is globally attractive. For any initial point $X(0) = (W_E(0), W_M(0), W_F(0), G_M(0))^T \notin \Omega$, we study the trajectory trend of $X(t, 0, X(0))$ with the increase of time $t$. If $W_E(0)\leq K$, similar to the above discussion above the invariant set, we can get $0\leq W_E(t)\leq K, t\geq0$. If $W_E(0) > K$, from the first equation of system (2.1), we get $\frac{dW_E}{dt}|_{W_E\geq K} < -(\rho +\mu_1)K < 0$, and there must exist a time $t_1 > 0$ such that $0\leq W_E(t)\leq K, t\geq t_1$. When $t\geq t_1$, by the second equation of system (2.1) and $W_E(t)\leq K$, we have $\frac{dW_M}{dt}\leq \theta \rho K-\mu_2W_M$. If $W_M(t_1)\leq L_1$, similar to the above discussion above the invariant set, we can get $0\leq W_M(t)\leq L_1, t\geq t_1$. If $W_M(t_1) > L_1$, we have $\frac{dW_M}{dt}|_{t\geq t_1}\leq \theta \rho K-\mu_2W_M,$ then there must exist a time $t_2 > t_1$ such that $0\leq W_M(t)\leq L_1, t\geq t_2$. When $t\geq t_2$, by the third equation of system (2.1) and $W_E(t)\leq K$, we have $\frac{dW_F}{dt}\leq (1-\theta)\rho K b-\mu_3W_F$. If $W_F(t_2)\leq L_2$, similar to the above discussion above the invariant set, we can get $0\leq W_F(t)\leq L_2, t\geq t_2$. If $W_F(t_2) > L_2$, we have $\frac{dW_F}{dt}|_{t\geq t_2}\leq (1-\theta)\rho K b-\mu_3W_F,$ then there must exist a time $t_3 > t_2$ such that $0\leq W_F(t)\leq L_2, t\geq t_3$. Further more, since $\lim \limits_{t\rightarrow +\infty}G_M(t) = \tilde{G}_M(t)$ and $0\leq \tilde{G}_M(t)\leq L_3, t\geq 0$, we deduce that $\Omega$ is a globally attracting set of system (2.1) in $R^4_+$. This completes the proof.

Theorem 1. System (2.1) has a wild mosquito-free periodic solution $(0, 0, 0, \tilde{G}_M(t))$, where

and

Proof. According to the positive periodic solution of model (2.3) obtained in Proposition 2, we can easily get $(0, 0, 0, \tilde{G}_M(t))$ is a periodic solution of system (2.1) which implies the eradication of wild mosquitoes. The proof is completed.

We first prove the local stability of the $\omega-$period solution $(0, 0, 0, \tilde{G}_M(t))$.

Theorem 2. The wild mosquito-free periodic solution $(0, 0, 0, \tilde{G}_M(t))$ of system (2.1) is locally asymptotically stable.

Proof. In order to study the local stability of the wild mosquito-free periodic solution $(0, 0, 0, \tilde{G}_M(t))$, we consider the following subsystem

Obviously, system (2.6) has a trivial equilibrium $(0, 0, 0)$. Computing the Jacobian matrix of system(2.6) at $(0, 0, 0)$, we get

and it has three real and negative eigenvalues $\lambda_i, i = 1, 2, 3$. Then the Floquet multipliers of the corresponding monodromy matrix $e^{J_O \omega}$ are $e^{\lambda_i \omega} < 1, i = 1, 2, 3$. It follows from the Floquet theorem that the trivial equilibrium $(0, 0, 0)$ of system (2.6) is always locally stable, which also implies that the wild mosquito-free periodic solution $(0, 0, 0, \tilde{G}_M(t))$ of system (2.1) is locally asymptotically stable. The proof is completed.

In the following, we determine conditions under which the wild mosquito-free periodic solution $(0, 0, 0, \tilde{G}_M(t))$ is also globally attractive.

Theorem 3. The wild mosquito-free periodic solution $(0, 0, 0, \tilde{G}_M(t))$ of system (2.1) is a global attractor provided one of the following conditions holds

(i) $\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho} \geq b$;

(ii) $\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho} < b$ and $\frac{1}{\alpha} \big(\frac{(b-\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho})^2}{4b\frac{\mu_2}{\theta \rho K}\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho}} -\gamma\big) < \frac{\delta \exp(-\mu_4\omega)}{(1-\exp(-\mu_4\omega))}$.

Proof. Since $\lim \limits_{t\rightarrow +\infty}G_M(t) = \tilde{G}_M(t)$, to discuss the global attractivity of $(0, 0, 0, \tilde{G}_M(t))$, we only need to find conditions that can guarantee the global stability of $(0, 0, 0)$ for system (2.6).

According to Theorem 2, the trivial equilibrium $(0, 0, 0)$ of system (2.6) is locally stable. We now exploit conditions under which it is globally attractive.

By letting the sterilizing males $\tilde{G}_M(t)$ be a constant, that is, $\tilde{G}_M(t)\triangleq G_M^{cst}\geq 0$, we construct a comparison system as follows

We calculate the Jacobian matrix of system (2.7) as follows

Obviously, all the off-diagonal elements are non-negative on the set $\Omega_1 = \big\{(W_E, W_M, W_F)\in R_+^3:W_E\leq K\big\}$, which implies that system (2.7) is monotone on $\Omega_1$ in the sense of the monotone systems theory ^[39]. Besides, it has a trivial equilibrium $O(0, 0, 0)$ which is locally stable. To verify the existence of positive steady state, we need to solve the following algebraic equations

By direct calculation, we get

For simplicity, denote

then the number of positive steady states of system (2.7) equals the number of positive roots of the following equations with respect to $x$:

For the quadratic equation in system (2.9), it is straightforward to show that it has no positive root if $(i) B_2\geq b$ or $(ii) B_2 < b$ and $\frac{1}{\alpha} \big(\frac{(b-B_2)^2}{4bB_1B_2} -\gamma\big) < G_M^{cst}$ holds.

If $B_2 < b$ and $\frac{1}{\alpha} \big(\frac{(b-B_2)^2}{4bB_1B_2} -\gamma\big) > G_M^{cst}$, then the quadratic equation in system (2.9) has positive roots

which must satisfy the second inequality in system (2.9). That is to say, system (2.7) must have two positive equilibria, one of which is locally stable.

Due to the monotonicity of the system, the trivial equilibrium $O(0, 0, 0)$ of system (2.7) is globally asymptotically stable if it is the unique steady state. According to the above analysis, when $B_2 < b$, there is a critical value $G_M^{crit} = \frac{1}{\alpha} \big(\frac{(b-B_2)^2}{4bB_1B_2} -\gamma\big)$ for $G_M^{cst}$ and $G_M^{cst} > G_M^{crit}$ can ensure that $O(0, 0, 0)$ is a global attractor.

Based on the analytical expression of $\tilde{G}_M$ in system (2.4), we can easily get its upper and lower bounds

If $\tilde{G}_M^L = \frac{\delta \exp(-\mu_4\omega)}{1-\exp(-\mu_4\omega)} > G_M^{crit}$, then we have $\tilde{G}_M(t) > G_M^{crit}, t\geq 0$. By the monotonicity and the relation between systems (2.6) and (2.7), we know that the trivial equilibrium $O(0, 0, 0)$ is globally asymptotically stable for system (2.6) if it is globally asymptotically stable for system (2.7).

Thus, the trivial equilibrium $O(0, 0, 0)$ is globally asymptotically stable for system (2.6) if $(i) B_2\geq b$ or $(ii) B_2 < b$ and $\tilde{G}_M^L > G_M^{crit}$ holds, then wild mosquito-free periodic solution $(0, 0, 0, \tilde{G}_M(t))$ of system (2.1) is globally asymptotically stable under the same conditions. The proof is completed.

Remark 1. In Theorem 3, $\frac{1}{B_2} = \frac{\beta (1-\theta)\rho}{(\rho +\mu_1)\mu_3}$ involves with the fecundity of the wild mosquitoes in the field. According to the results in Theorem 3, if the fertility is weak enough, that is, $b\frac{1}{B_2}\leq 1$, then the wild mosquitoes will eventually go extinct even without human intervention. If the fertility is relatively strong, for example, $b\frac{1}{B_2} > 1$, we can also eliminate the wild mosquitoes in the long run by adjusting the intensity of releases of sterile males.

2.3.   Large-scale time control for wild mosquitoes

In the following, we study the large-scale time control strategies for wild mosquitoes by theoretical analysis and give practical methods for selecting the release amount $\delta$ and release period $\omega$ so that the wild mosquitoes can be wiped out from the field.

Based on the first kind of condition listed in Theorem 3, if the fertility in a given field is weak $(\frac{\beta (1-\theta)\rho}{(\rho +\mu_1)\mu_3}\leq \frac{1}{b})$, wild mosquitoes will always go extinct without any human interventions. So we mainly consider a more common case when wild mosquitoes has a relatively strong fertility and we need release sterilizing males reasonably into the field to wipe out the wild ones.

According to the second kind of conditions listed in Theorem 3, we discuss the control strategies if $\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho} < b$ holds.

Obviously, if $\frac{(b-\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho})^2}{4b\frac{\mu_2}{\theta \rho K}\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho}} \leq\gamma$, then the conditions listed in $(ii)$ of Theorem 3 are valid and the wild mosquitoes will eventually go extinct without human intervention. If $\frac{(b-\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho})^2}{4b\frac{\mu_2}{\theta \rho K}\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho}} > \gamma$, we denote $\mathcal{N}(\delta, \omega) : = \frac{\delta \exp(-\mu_4\omega)}{(1-\exp(-\mu_4\omega))}$ and consider the equation

For any given release period $\omega^*$, $\mathcal{N}(\delta, \omega^*)$ is monotonically increasing with respect to $\delta$. Besides, by simple calculation, we can get $\mathcal{N}(0, \omega^*) = 0$ and $\mathcal{N}(+\infty, \omega^*) = +\infty$. Then there is a unique $\tilde{\delta} \ge 0$ satisfying $\mathcal{N}(\delta, \omega^*) = \frac{1}{\alpha} \big(\frac{(b-\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho})^2}{4b\frac{\mu_2}{\theta \rho K}\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho}} -\gamma\big)$. Hence it follows from Theorem 3 that the wild mosquito-free periodic solution $(0, 0, 0, \tilde{G}_M(t))$ of (2.1) is globally stable provided $\delta > \tilde{\delta}$ holds.

Similarly, for a given release amount $\delta^*$, $\mathcal{N}(\delta^*, \omega)$ is monotonically decreasing with respect to $\omega$. Since it is obvious to have $\mathcal{N}(\delta^*, 0) = +\infty$ and $\mathcal{N}(\delta^*, +\infty) = 0$, there exists a unique $\tilde{\omega} \ge 0$ such that $\mathcal{N}(\delta^*, \omega) = \frac{1}{\alpha} \big(\frac{(b-\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho})^2}{4b\frac{\mu_2}{\theta \rho K}\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho}} -\gamma\big)$. Thus the wild mosquito-free periodic solution $(0, 0, 0, \tilde{G}_M(t))$ of system (2.1) is globally stable provided $\omega < \tilde{\omega}$ according to Theorem 3.

In the following, we will investigate the release tactics for large-scale time control of wild mosquitoes by numerical simulations. Most model parameters are chosen from ^[8] and ^[17] (refer to Table 1). While the environmental carrying capacity $K$ and the strength of Allee effect $\gamma$ remain pending for their values often change with environments.

In this paper, we consider parameters as follows

For the wild mosquito extinction induced by a strong Allee effect, we select $\gamma = 1500$ and there is $\frac{(b-\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho})^2}{4b\frac{\mu_2}{\theta \rho K}\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho}} = 1240.1868 < \gamma = 1500$. Then the wild mosquitoes eventually go extinct without sterile males deliveries (see Figure 1).

Keep other parameters the same as in Eq (2.11) but change $\gamma = 1500$ to $\gamma = 200$, then we have $\frac{(b-\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho})^2}{4b\frac{\mu_2}{\theta \rho K}\frac{(\rho +\mu_1)\mu_3}{\beta (1-\theta)\rho}} > \gamma$. By direct calculation, we get

We first fix the release period $\omega^* = 3$ and through simple calculation we obtain that the unique positive root of $\mathcal{N}(\delta, \omega^*) = 1040.1868$ is $\tilde{\delta} = 1161.8887$. To verify our theoretical results, we compare two release amounts $\delta = 1200$ and $\delta = 800$ in Figures 2 and 3, respectively. We can see that if the release amount $\delta = 1200 > \tilde{\delta}$, then the wild mosquito-free periodic solution of system (2.1) is globally stable (see Figure 2), while if the release amount $\delta = 800 < \tilde{\delta}$, then there exists a locally stable positive coexistence period solution in addition to the wild mosquito-free periodic solution (see Figure 3).

We then fix the release amount $\delta^* = 1000$ and get that the unique positive root of $\mathcal{N}(\delta^*, \omega) = 1040.1868$ is $\tilde{\omega} = 2.6946$. Similarly, we compare two release periods $\omega = 2.5$ and $\omega = 4$ in Figures 4 and 5, respectively. We see that if the release period $\omega = 2.5 < \tilde{\omega}$, then the wild mosquito-free periodic solution of system (2.1) is globally stable (see Figure 4), while if the release period $\omega = 4 > \tilde{\omega}$, then a locally stable positive coexistence period solution coexists with the wild mosquito-free periodic solution (see Figure 5).

In this section, we focus on the asymptotic behaviors of system (2.1), and verify the theoretical results by numerical simulations. From Figures 1–5, we can see that although the sequential impulsive releases of sterile mosquitoes can make the wild mosquitoes go extinct, there are disadvantages in terms of cost control. After the wild population entered the basin of attraction of the extinction solution, the extinction of wild mosquitoes is a foregone conclusion, and the releases of sterile mosquitoes only speed up the process. This is unreasonable from the perspective of cost control. Because the cost control is not involved in this part of the study, the conclusion of this part is more limited to the theoretical discussion, and ignores the rationality of practical application.

3.   Optimal release strategies to control wild mosquitos in a limited time

As pointed out in Section 2, regardless of the state of wild mosquito population in the environment, blind sequential releasing sterile mosquitoes in large-scale time can ensure the extinction of wild mosquitoes, but it will cause unnecessary cost waste in practice, which is not desirable. Furthermore, the prevalence of mosquito borne diseases describes obviously seasonal and regional characteristics. To this end, in this section we investigate limited-time optimal control of wild mosquitoes incorporating both of the population control level of wild mosquitoes and the economic input, and study three different release strategies: optimal release amount for periodic releases, optimal release timing for a fixed release amount and mixed control with optimal release timing and amount each time.

3.1.   Optimization by mixed control

Firstly, we consider a relatively complicated scenario when both release timings and release amounts are chosen as control parameters. Assume that we need to make mosquito population under control in a predefined finite interval $[0, T]$ and $N-1$ releases of sterilizing males are planed. Suppose the release amount and release timing of the $i$th release are $\delta_i$ and $t_i\in [0, T]$, respectively, with $i = 1, 2, \cdots, N-1$. Then a limited-time control system is proposed in the following

with initial conditions

$T$, the length of control time, can be converted to be dependent on the transmission of mosquito borne infectious diseases in this situation. $t_i, i = 1, 2, \cdots, N-1$, the release timings, are assumed to satisfy $0 = t_0\leq t_1\leq t_2\leq \cdots \leq t_{N-1}\leq t_N = T$. Denote $T_i = t_i-t_{i-1}$ which represents the time interval between the $(i-1)$th and $i$th release. According to practical meaning, this time interval cannot be too long or too short. Thus we give constraints

where $\tau_i^1$ and $\tau_i^2$ are given constants. Similarly, we also give constraint for the release amount $\delta_i$

where $\delta_i^1$ and $\delta_i^2$ are given constants which represent the minimum and maximum amount allowed in the $i$th release. Denote $\Gamma = (T_1, T_2, \cdots, T_N)^T$ and $\delta = (\delta_1, \delta_2, \cdots, \delta_{N-1})^T$, where $T_i$ and $\delta_i$ meet the stated constraints Eqs (3.3) and (3.4), respectively. Let $\Psi_1$ and $\Psi_2$ be sets of all $\Gamma \in R^N, \delta\in R^{N-1}$ satisfying Eqs (3.3) and (3.4), respectively.

Since the right hand side of system (3.1) are differentiable, system (3.1) with initial condition Eq (3.2) has a unique solution $(W_E(t), W_M(t), W_F(t), G_M(t))^T$ corresponding to each pair $(\Gamma, \delta)\in (\Psi_1, \Psi_2)$ ^[37,38].

Based on our key considerations of the mosquito population control, we define a cost function as follows

Here, $r_0$ stands for the per unit cost of rearing sterilizing mosquitoes.

For the optimal control problem in this scenario, we can state it formally as follows.

${(P1)}$ Subject to the dynamical system (3.1) with initial condition Eq (3.2), find a feasible parameter vector pair $(\Gamma, \delta)$ such that the cost function $\mathfrak{J}(\pi, \delta)$ is minimized over $(\Psi_1, \Psi_2)$.

Since the state of variables $W_E(t), W_M(t), W_F(t)$ and $G_M(t)$, functions of $t, \Gamma$ and $\delta$, depends on uncertain release timings and uncertain release amounts, the optimal control problem ${(P1)}$ cannot be solved directly by general optimization techniques. In this paper, we apply a time rescaling method which has been used in several studies ^{[33,34,40,41]} to transform these uncertain pulse time points into fixed ones. By this method, the optimal problem ${(P1)}$ is turned into an equivalent optimal parameter selection problem which is described as a series of ordinary differential equations with periodic initial conditions. And the new equivalent problem can be solved by utilizing a gradient-based optimization technique.

To this end, let $t = \sum_{j = 1}^{i-1}T_j+T_i s$, $t\in(\sum_{j = 1}^{i-1}T_j, \sum_{j = 1}^{i}T_j]$, and denote

Then system (3.1) with initial condition Eq (3.2) is converted into $N$ subsystems

with

The cost function Eq (3.5) can also be redefined as follows

Based on the above transformations, we can restate the problem ${(P1)}$ as follows

${(P2)}$ Subject to the dynamical system (3.7) with initial condition Eq (3.8), find a feasible parameter vector pair $(\Gamma, \delta)$ such that the cost function defined in Eq (3.9) is minimized over $(\Psi_1, \Psi_2)$.

Using Theorem 6.1 in ^[42], we define Hamiltonian functions $H_i, i = 1, 2, \cdots, N$ in the following

Here $\lambda^i(s) = (\lambda_1^i(s), \lambda_2^i(s), \lambda_3^i(s), \lambda_4^i(s))$ are costate variables which are governed by the following costate equations

with

Define

and from system (3.7) there is

From Theorems 4.1 and 4.2 in ^[41], we have

Proposition 3. The gradients of the cost functional $\mathfrak{J}_1$ with respect to $\delta$ and $\Gamma$ are given by

and

respectively.

By straightforward calculation, we get

Theorem 4. The gradients of the cost function $\mathfrak{J}_1(\Gamma, \delta)$ with respect to the release timing $T_i$ and release amount $\delta_l$ are given by

for $i = 1, 2, \cdots, N$, and

for $l = 1, 2, \cdots, N-1$, respectively.

3.2.   Optimization by release amount

In this subsection, we study a relatively simple but common scenario when sterile male mosquitoes are released into the field periodically with a fixed release amount. To determine an optimal release amount for this mode, we suppose that sterile mosquitoes are periodically released with a constant release amount $\delta_d$ in the limited time $[0, T]$ and $N-1$ times of releases are totally planed, that is to say, the release period is $\omega = \frac {T}{N}$. Then the limited-time control system is proposed as follows

with initial conditions Eq (3.2).

Just like we did in section 3.1, we also give a constraint for the release amount $\delta_d$

where $\delta_{low}$ and $\delta_{up}$ are given constants which represent the minimum and maximum amount that are allowed in each release.

In this scenario, we can define the cost function for the control problem $(P1)$ as follows

Here, $\delta_d$ is the only control parameter. That is, we need to determine a release amount $\delta_d$ such that $\mathcal{\bar{J}}(\delta_d)$ is minimized over $[\delta_{low}, \delta_{up}]$.

Apply the time rescaling technique again and let $t = (i-1)\omega+s\omega$ for $i = 1, 2, \cdots, N$, then the system (3.15) with initial condition Eq (3.2) is converted into the following $N$ subsystems

with the same initial conditions Eq (3.8).

Then the cost function Eq (3.17) can be redefined as

while the optimal control problem can be restated as: determine a $\delta_d$ such that $\mathfrak{\bar{J}}_1(\delta_d)$ is minimized over $[\delta_{low}, \delta_{up}]$.

Then by similar discussion, we get the following result.

Theorem 5. The gradient of $\mathfrak{\bar{J}}_1(\delta_d)$ with respect to the release amount $\delta_d$ is

3.3.   Optimization by release timing and uniform release amount

In this subsection, the release timings are added as new control parameters on the basis of the preceding scenario. That is, sterile mosquitoes are released at irregular moments $0\leq t_1\leq t_2\leq\cdots\leq t_{N-1}\leq T$ with a same release amount $\delta_d$. To find a set of optimal release timings and an optimal release amount for this case, the following limited-time control system is proposed

with the same initial conditions listed in Eq (3.2). Besides, release timings $t_i, i = 1, 2, \cdots, N-1$ and release amount $\delta_d$ meet the stated constraints Eqs (3.3) and (3.16), respectively.

We define the cost function of control problem $(P1)$ as follows

with $\Gamma = (T_1, T_2, \cdots\cdots, T_N)^T, T_i = t_i-t_{i-1}$.

Use the time rescaling technique and let $t = \sum_{j = 1}^{i-1}T_j+T_i s$ for $i = 1, 2, \cdots, N$, then the system (3.21) is turned into

with the same initial conditions listed in Eq (3.8).

Accordingly, the cost function Eq (3.22) is transformed into an equivalent form

Then we obtain the following result.

Theorem 6. The gradients of $\mathfrak{\hat{J}}_1(\Gamma, \delta_d)$ with respect to the release timing $T_k$ and release amount $\delta_d$ are given by

and

respectively. Here $\nabla_{T_k}\mathfrak{J}_1(\Gamma, \delta)$ is the same as Eq (3.13).

4.   Numerical simulations for the optimal control

To determine the optimal values of the release timings and release amounts for the three limited-time optimal control problems in the preceding section, a series of numerical simulations are performed in this section. These three optimal release strategies will also be compared in different ways.

Before proceeding further, we need to introduce the calculation method of the cost function and its gradients with respect to the control parameters in detail. This method has been used in ^[34] and ^[41], and we will explain it in the following by the case presented in Section 3.1.

(ⅰ) We firstly solve straightforward the differential equations (3.7) with initial conditions Eq (3.8) to obtain $W_E^i(s), W_M^i(s), W_F^i(s), G_M^i(s), s\in [0, 1]$ for $i = 1, 2, \cdots, N$.

(ⅱ) Using $W_E^i(s), W_M^i(s), W_F^i(s), G_M^i(s)$ obtained in last step, we first solve backwards the costate equations (3.11) with boundary conditions Eq (3.12) and obtain the costate variables $\lambda_1^i(s), \lambda_2^i(s), \lambda_3^i(s)$ and $\lambda_4^i(s)$ for $i = N$. Then we obtain the costate variables for $i = N-1$ in the same way, and continue until we obtain the costate variables for $i = 1$.

(ⅲ) Based on the expression in Eq (3.9), we compute the cost function $\mathfrak{J}_1(\Gamma, \delta)$ by using $W_M^N(s)$ and $W_F^N(s)$ and release amounts $\delta_i$.

(ⅳ) Applying $W_E^i(s), W_M^i(s), W_F^i(s), G_M^i(s), \lambda_1^i(s), \lambda_2^i(s), \lambda_3^i(s)$ and $\lambda_4^i(s)$ obtained in step $(i)$ and $(ii)$, we calculate $\nabla_{T_i}\mathfrak{J}_1(\Gamma, \delta)$ for $i = 1, 2, \cdots, N$ and $\nabla_{\delta_l}\mathfrak{J}_1(\Gamma, \delta)$ for $l = 1, 2, \cdots,$$N-1$.

Keep the parameter values in Eq (2.11) and consider the Allee effect in a relatively high level with $\gamma = 200$. Besides, suppose that the average cost of rearing per unit of sterilizing mosquitoes $r_0 = 0.01$. When doing numerical simulation, we measure the time in days and take 20 days as the total control time, that is to say, $T = 20$, and 4 releases of sterilizing males are planed, then these 20 days should be divided into $N = 5$ parts according to different rules.

In the following, we will study three different optimal strategies in impulsive control by numerical simulations. Surely there is no guarantee that the optimal solution we find numerically is unique, so we just present some optimal ones with special initial release periods and amounts. Specifically, with the help of nonlinear optimization-Matlab function library, we will find the optimal parameters by Matlab according to the objective function and the correlation gradient calculated in the above steps.

Example 1. Optimal release amount for periodic releases

It is obvious that the release period is $\omega = \frac{T}{N} = 4$. Starting with a initial release amount $\delta_d = 600$, if no optimal control is taken and only simple impulsive releases are employed, we can obtain that after five periods the cost value is $\mathfrak{J}_0 = 371.9008$ and the total fertile wild mosquito population is $W_M(T)+W_F(T) = 347.9008$ at $T = 20$.

Under the constraint $0\leq \delta_d \leq 1000$, we solve the corresponding optimal problem numerically by using the algorithm listed above in Matlab. We get an optimal release amount $\delta^*_d = 559.54$ and the corresponding cost value $\mathfrak{J}^* = 371.8468$, while the total fertile wild mosquito population is $W_M^*(T)+W_F^*(T) = 349.4652$ at time $T = 20$. We plot the time series diagrams of fertile wild mosquitoes for this kind of optimal control, non-control and simple impulsive control in Figure 6(a). By these three curves, we find that although the simple impulsive control has an obvious superiority in reducing fertile mosquito population (in most time of the control process, the total fertile wild mosquito population of the optimal one is a little higher than that of the simple impulsive control one), it cost more sterile mosquitoes to achieve such an effect (see Table 2).

In addition, we investigate the influence of release amount on the cost function and the number of the fertile wild mosquito population at time $T$ (see Figure 6(b)), and find that when the release amount varies in the interval $0\leq \delta_d \leq 1000$, the cost function $\mathfrak{J}(\delta_d)$ admits a minimum point, which also verifies the optimum result we obtained above. Furthermore, we notice that the increase of the release amount leads to the population reduction of fertile wild mosquitoes at the terminal time $T$. However, the wild mosquito cannot be eliminated regardless of the releasing amount of sterile mosquitoes each time.

Example 2. Optimal release timing for a fixed release amount

To be consistent with Example 1, we choose the same initial release amount $\delta_d = 600$ and select $T_1 = 2.5, T_2 = T_3 = T_4 = 4, T_5 = 5.5$ as the initial release intervals. In order to determine optimal time intervals $T_i$ and optimal releasing amount $\delta_d$ which can minimize the cost function $\mathfrak{J}$, we consider constraint conditions

and $0\leq \delta_d\leq 1000$.

Then by using the Matlab program, we solve this optimal problem numerically and obtain the following optimal release intervals

and an optimal release amount

In addition, we get the minimum cost value $\mathfrak{J}^* = 371.0832$ and the total fertile wild mosquito population $W_M(T)+W_F(T) = 348.1168$ at $T = 20$.

We plot the time series diagrams of fertile wild mosquitoes for the optimal release timing control, optimal release amount control and non-control in Figure 7(a). From these curves, we find that the optimal release timing control has a better control effect with a relatively low cost function value. Besides, for every $\delta_d \in [0, 1000]$, we solve the corresponding optimal time intervals under the restriction Eq (4.1) and then calculate the value of cost function and the amount of total fertile wild mosquitoes at $T = 20$. We find that, see Figure 7(b), when the release amount varies in the interval $0\leq \delta_d \leq 1000$, the cost function $\mathfrak{J}(\Gamma, \delta_d)$ also admits a minimum point, which agrees with the optimal values we have obtained. Similarly, the increase of the release amount reduces the fertile wild mosquito population, but the wild mosquitoes cannot be eliminated even if the release amount reaches the upper constrained bound and the cost value is very high.

Example 3. Optimal release timing and release amounts

Keep the initial release intervals $T_1 = 2.5, T_2 = T_3 = T_4 = 4, T_5 = 5.5$ and choose the initial release amounts $\delta_1 = \delta_2 = \delta_3 = \delta_4 = 600$, we deal with the optimal problem with constraints Eq (4.1) and $0\leq \delta_i\leq 1000, i = 1, 2, 3, 4$. Solving this optimal problem numerically in Matlab, we obtain the optimal release amounts

and the optimal release intervals

This release strategy is showed in Figure 8(b). Besides, we obtain the minimum cost value $\mathfrak{J}^* = 371.0803$ and the total fertile wild mosquito population $W_M(T)+W_F(T) = 347.6869$ at $T = 20$. The time series diagrams of the number of total fertile wild mosquitoes under four types of control modes are plotted in Figure 8(a) from which we can see that mixed optimal control produces the best control effect.

Finally, we compare these three optimal release strategies (refer to Table 3 and Figure 9). We find that the optimal release timing control is superior to the optimal release amount control, while the mixed control produces the best integrated control effect since the lowest fertile wild mosquito level is reached at the minimal cost function value. From Figure 9(b), we also see that the mixed optimal control releases the most sterile mosquitoes in the whole control process, but its cost function value is not large because the least wild mosquitoes stayed in the field at the terminal time. Although it can be seen from Figure 9(b) that the largest accumulated amount of sterile mosquitoes are released during the whole control process in the mixed optimal control, its cost function value is the smallest due to the smallest final population size of wild mosquitoes in the field.

5.   Conclusions

The SIT has been a hot topic in the research field of mosquito-borne infectious disease control in recent years, and lots of field trails have been conducted all over the world and large numbers of researchers have been devoting themselves in this area and have already achieved many progresses. Mathematical model, as an important tool, plays a significant role in the research process. However, most of these models in previous studies are constructed by continuous or discrete dynamical systems, which cannot describe the release process accurately. Most works focused on the asymptotic behavior of the system in infinite time. However, the control of mosquitoes in most cases should be a shorter-term action.

Release of sterile mosquitoes has been used to reduce or eliminate the wild mosquito population in order to control vector-borne infectious diseases. In this paper, we proposed and studied a stage-structured two-sex mosquito population model with an Allee effect and impulsive releases of sterile males. By adjusting different types of control parameters, both large-scale time control and limited-time optimal control of wild mosquitoes were investigated.

We firstly studied the large-scale time control aiming to wipe out wild mosquitoes. By using the monotone system theory and the comparison theorem, we showed the existence, uniqueness and globally stability of the wild mosquito-free periodic solution. For fixed release period $\omega^*$ (or release amount $\delta^*$), we established threshold value for release amount $\tilde{\delta}$ (or release period $\tilde{\omega}$) which determines the extinction or persistence of the wild mosquito population.

Then for the limited-time optimal control of wild mosquitoes, we took into account both of the population control level of wild mosquitoes and the economic cost, and investigated three different release tactics: optimal release amount for periodic releases, optimal release timing for a fixed release amount and a combination of optimal release timing and release amounts. A time rescaling technique was applied to overcome the technical difficulty that the state of variables depends on uncertain pulse effects. We obtained the optimal release amounts and release timings numerically for each release strategy. Numerical simulations indicate that the optimal release timing control is a more effective strategy than the optimal release amount control. However, simultaneous optimal selection of release amount and release timing leads to the best control performance.

In the limited-time control, we construct cost functions by referring to the pest control in agriculture and only focus on terminal control but ignore process control. And in our future work, we will consider both terminal control and process control in the limited time control, so as to ensure that the number of wild mosquitoes cannot be too large during the control process.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (12071407, 11901502 and 11871415), Training plan for young backbone teachers in Henan Province (2019GGJS157), Foundation of Henan Educational Committee under Contract (21A110022), Program for Science $&$ Technology Innovation Talents in Universities of Henan Province (21HASTIT026), Scientific and Technological Key Projects of Henan Province (212102110025), Nanhu Scholars Program of XYNU and Nanhu Scholars Program for Young Scholars of XYNU.

Conflict of interest

The authors declare that they have no competing interests.