Mathematical Modelling of the Sterile Insect Technique Using Different Release Strategies

Abstract

We study simple mathematical models for the dynamics of interactive wild and sterile insect populations. As well as being mathematically tractable, these models can be used as first approximations to real situations occurring with the Sterile Insect Technique (SIT) in which sterile males are released to reduce or eradicate a pest population. This is a method of biological control which can effectively help contain the spread of many pest insects such as the Red Palm Weevil (RPW). Models formulated in this paper are continuous-time, include a strong Allee effect that captures extinction events, and incorporate different strategies of releasing sterile insects. We perform basic studies of dynamical features of these models, with an emphasis on the condition of excitation, and the impact of the different release methods is investigated. Our findings are also demonstrated with some numerical examples.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In recent years, there has been a rapid rise in the use of biological methods for the control of insect pests. One tool, which has proved effective in the area-wide control of various insects, is the SIT. This method, introduced initially by Knipling [1], consists in releasing high numbers of sterilized males into the environment. Such a technique constitutes a biological control process that disturbs the natural reproduction of the insect pests. This is carried out by using chemical or physical or other radical procedures to treat male insect pests to make them infertile, so they become unable to reproduce regardless of their sex drive. The infertile males are then introduced to the environment and compete for mates with fertile males, such that interaction between sterilized males and any female wild insect pests will not lead to any insect reproduction, thereby disrupting the natural reproductive process of the population. Despite the fact that frequent release of treated insect pests in large amounts will eventually eradicate the wild pest population completely, it is more practical to control the wild insect pest population instead of eliminating it completely. Then, if the number of released sterile males is high enough and is repeated over a sufficient period of time, the average fertility of the target population could be reduced leading to the control or even the eradication of the pest population from large areas. The first successful SIT operation was against the screwworm population in Florida in the late 1950s. Later, this technique has been applied to combat a number of pests and disease vectors, such as the Mediterranean fruit fly (medfly), the RPW in coconut and date palm gardens, and the tsetse fly in Africa (see [2] for an overall review of the SIT and its applications). Moreover, diseases like dengue fever and malaria that are transmitted to humans by blood-feeding mosquitoes present a significant health concern for people. Roughly around one to three million people every year succumb to malaria as indicated by the World Health Organization (WHO). Malaria vastly hits Africa and South America, majorly taking the lives of children and pregnant women. As there are no vaccines available to prevent mosquito-borne diseases, the only way to prevent these diseases is to control the mosquitoes.

On the other hand, for vector control in particular, new approaches with similar working principles as the SIT have been developed. Those include, on the one hand, genetically modified control methods, such as the RIDL (Release of Insects with Dominant Lethality) technique, and, on the other hand, the Wolbachia technique. The former involve the release of genetically engineered insects (that have a lethal gene in their genome in the RIDL strategy), while the latter utilizes the Cytoplasmic Incompatibility (CI) property of the Wolbachia bacterium [3–5]. Indeed, these bacteria have the property to alter the sperm of infected males making it unable to fertilize uninfected eggs. This is the principle of the Incompatible Insect Technique (IIT) [6–11]. Moreover, the CI property raises considerably the progeny of infected females. And since Wolbachia is maternally inherited, releasing high numbers of W-females into a target population may lead to a Population Replacement (PR) by Wolbachia-infected insects and eventually to the elimination of the wild population (see, for example, [12] for a recent review of the Wolbachia-based PR strategy). Note that PRs and invasions have been observed in natural populations, such as with the Californian Culex pipiens [13] and the Australian Aedes aegypti [14].

Motivated by the issue of controlling a pest population by means of an SIT-like method, numerous theoretical studies, especially on the mathematical modelling of the classical SIT, the IIT, and the Wolbachia-PR, have been conducted (see, for example, [15–26] for SIT/IIT and [27–30] for PR and references therein). As a matter of fact, mathematical models have proven valuable in understanding various important issues in population dynamics, such as suppression mechanisms and the success or failure of different strategies. Thus, various classes of models have been formulated, including deterministic, stochastic, continuous-time, discrete-time, hybrid approaches, and temporal and spatiotemporal models.

In this paper, we study the dynamics of the interactive wild and sterile insects with a particular focus on the impact of the strategy adopted in releasing sterile individuals. Three release methods are then incorporated based on works [25, 31, 32]. The sterile-fertile interaction is assumed to be a one-sided competition that affects only the wild-type population. To reflect the need of a critical threshold density for the persistence of the wild population, a strong Allee effect is included. Moreover, to keep the model reasonably simple, we consider homogeneous insect populations such that no male-female or stage distinction is made and death rates for sterile and fertile insects are assumed to be density-independent and equal.

The paper is organized as follows. In Section 2, we present our general modelling assumptions. In Sections 3, 4, and 5, respectively, we consider three submodels, each with a different strategy of release: the first involves a constant release rate, the second assumes a release rate proportional to the size of the wild population, and the third uses a release rate of Holling-II type. We carry out detailed mathematical analysis of these models and discuss their dynamical features, especially the existence of equilibria and their stability. We also illustrate our analytical findings with numerical calculations. In the final section, we give a brief conclusion.

2. The General Model

We consider a two-dimensional one-stage model that involves density dependence solely in the growth term of pest insects. We assume that the birth rate of the sterile insects is their release rate and that the sterile-fertile competition affects only the wild population.

Let $Nt$ be the number of wild insects at time $t$. In the absence of sterile insects, $Nt$ is assumed to evolve according to the following equation:$$\begin{array}{}\text{(1)}& \frac{\text{d}N}{\text{d}t}=GN-\mu N,\hspace{1em}t\ge 0,\end{array}$$where $GN$ is the growth function of the wild population and $\mu$ is its death rate taken to be density independent. Furthermore, considering that the subsistence of the wild population requires some critical threshold density, we assume an Allee effect such that the birth function takes the form [33, 34]$$\begin{array}{}\text{(2)}& GN=\alpha rNK-NN-\epsilon ,\end{array}$$where $K$ is the carrying capacity of the niche and $ϵ$ is the Allee effect parameter $ϵ<K$, and we have set $\alpha =4/{K-ϵ}^2$, so that the parameter $r$ is the maximum value of the growth rate as in the usual logistic model. Note that the Allee effect has been attracting a lot of interest recently due to its strong potential impact on population dynamics (see, for instance, [35]). This effect may arise from different causes, but the most obvious of them is the difficulty of finding mates at low population sizes [36].$$\begin{array}{}\text{(3)}& \frac{\text{d}N}{\text{d}t}=\alpha rNK-NN-\epsilon -\mu N.\end{array}$$

This equation has a trivial equilibrium point at $N_0^{\ast }=0$, which is locally asymptotically stable. A positive equilibrium of equation (3) should be a zero of the function $QN$ defined by$$\begin{array}{}\text{(4)}& QN≔K-NN-\epsilon -\tilde{\mu },\end{array}$$with $\tilde{\mu }=\mu /\alpha r$. This function has two roots:$$\begin{array}{}\text{(5)}& N_{\pm }^{\ast }=\frac{K+\epsilon \pm K-\epsilon \sqrt{1-\mu /r}}{2},\end{array}$$which are real under the condition $r>\mu$, i.e., provided that the mortality rate remains below the maximum growth rate. Hence, the model given by equation (3) may have no, one, or two positive equilibria depending on whether $r>\mu$, $\mu =r$, or $r<\mu$, respectively.

In the case when the condition $r>\mu$ is fulfilled, at a positive equilibrium $N^{\ast }$ we have$$\begin{array}{}\text{(6)}& \frac{\text{d}}{\text{d}N}\alpha rNK-NN-\epsilon -\mu N{N}^{\ast }=\alpha r{N}^{\ast }K+\epsilon -2N^{\ast }.\end{array}$$

Thus,$$\begin{array}{}\text{(7)}& \frac{\text{d}}{\text{d}N}\alpha rNK-NN-\epsilon -\mu N{N}_{\pm }^{\ast }=\mp N_{\pm }^{\ast }K-\epsilon \sqrt{\frac{1-\mu }{r}},\end{array}$$which means that the equilibrium $N_{-}^{\ast }$ is unstable, whereas the equilibrium $N_{+}^{\ast }$ is asymptotically stable. From now, we assume that the condition $r>\mu$ is satisfied.

Now, after sterile insects are released throughout the wild population, its reproductive success will be reduced. We shall assume that the birth rate of wild insects is affected so that it follows the harmonic mean. On the other hand, as the sterile-fertile interaction is admitted to be a one-sided competition (for mates), sterile insects are not affected by the presence of fertile individuals. Thus, if we denote by $Mt$ the number of sterile insects at time $t$, the interactive dynamics are then governed by the following system:$$\begin{array}{c}\text{(8)}& \frac{\text{d}N}{\text{d}t}=\alpha r\frac{N^2}{N+M}K-NN-\epsilon -\mu N,\frac{\text{d}M}{\text{d}t}=RN-\mu M,\end{array}$$

where we have assumed that sterile insects have the same survivability as wild insects. $RN$ is a functional characterizing the release strategy, assumed to be independent of $t$. Subsequently, we will be considering three different forms of $RN$, along with the assumption $r>\mu$.

3. Constant Release Rate

We consider here the situation where sterile insects are constantly released so that $RN=\gamma$, with $\gamma$, a positive constant. Then, system (8) takes the form$$\begin{array}{c}\text{(9)}& \frac{\text{d}N}{\text{d}t}=\alpha r\frac{N^2}{N+M}K-NN-\epsilon -\mu N,\frac{\text{d}M}{\text{d}t}=\gamma -\mu M.\end{array}$$

For this model, it is easy to check that the rectangle of the phase plane$$\begin{array}{}\text{(10)}& {\Omega }_1≔N,M\in {ℝ}^2:0\le N\le N_{+}^{\ast },\hspace{1em}0\le M\le \eta ,\end{array}$$with $\eta ≔\gamma /\mu$, is a positive invariant, and we assume in this section that $N,M\in {\mathrm{\Omega }}_1$.

System (9) has a first equilibrium $E_0^c≔0,\eta$ lying on the boundary of ${\mathrm{\Omega }}_1$. It corresponds to the situation of no wild insects, but only a constant population of sterile individuals is present. Besides that, if it exists, a positive equilibrium should be a solution of the cubic equation:$$\begin{array}{}\text{(11)}& NQN-\tilde{\gamma }=0,\end{array}$$with $\tilde{\gamma }=\gamma /\alpha r$. By Descartes’ rule of signs, this equation has always one real negative solution (that we will disregard) and either no, one, or two real positive solutions. And one or two positive roots occur if the curve of the function $WN≔NQN$ intercepts the $\tilde{\gamma }$-axis once or twice, respectively. Furthermore, an elementary study of the variations of $WN$ shows that $WN<0$ for $N<N_{-}^{\ast }$. Then, as $N$ increases from $N_{-}^{\ast }$, this function also increases from zero to attain a certain maximum at the point$$\begin{array}{}\text{(12)}& N_h≔\frac{K+\epsilon +\sqrt{K\epsilon +{K-\epsilon }^{2}1-3\mu /4r}}{3},\end{array}$$before lowering indefinitely. Note here that $N_h$ is always real and positive as we have presumed that $r>\mu$. We then deduce that the condition for the existence of two positive solutions for equation (11) is$$\begin{array}{}\text{(13)}& N_{h}Q{N}_h>\tilde{\gamma }.\end{array}$$

This defines a threshold value of the release rate $\gamma$ as$$\begin{array}{}\text{(14)}& {\gamma }_c^c≔\alpha r{N}_{h}Q{N}_h,\end{array}$$at which equation (11) admits a unique positive root. Note that this definition of ${\gamma }_c^c$ tacitly requires that $Q{N}_h$ be positive, and we can easily show that this is still the case as $r>\mu$.

Thus, we infer that the model given by equation (9) will have no, one, or two positive equilibria if $\gamma >{\gamma }_c^c$, $\gamma ={\gamma }_c^c$, or $\gamma <{\gamma }_c^c$, respectively. Note also that all the equilibria of the model lie on the border of the rectangle ${\mathrm{\Omega }}_1$ at which $M=\eta$.

Next, we address the stability of the equilibria. The Jacobian matrix at an equilibrium point $N^{\ast },\eta$ has the form$$\begin{array}{}\text{(15)}& J^c≔\begin{array}{cc}{J}_{11}^c,& -\frac{\mu N^{\ast }}{N^{\ast }+\eta },\\ 0,& -\mu .\end{array}\end{array}$$

At the boundary equilibrium $E_0^c$, we have $J_{11}^c=-\mu$; therefore, $E_0^c$ is always asymptotically stable. Moreover, in the case where no positive equilibrium exists as ${\mathrm{\Omega }}_1$ is a positive invariant set, we conclude by the Poincaré–Bendixson theorem that $E_0^c$ is globally asymptotically stable. Thus, for a release rate of sterile insects $\gamma >{\gamma }_c^c$, the wild population evolves to extinction whatever its initial number is. This is a very important outcome of the model since estimating the optimal number of sterile insects to be released in order to eradicate the pest population remains one of the ultimate goals of the SIT modelling.

At a positive equilibrium $N^{\ast },\eta$, the matrix element $J_{11}^c$ can be written as$$\begin{array}{}\text{(16)}& J_{11}^c=\frac{\text{d}}{\text{d}N}\frac{\alpha rN}{N+\eta }WN-\tilde{\gamma }{N}^{\ast }.\end{array}$$

Then, keeping in mind that $W{N}^{\ast }=\tilde{\gamma }$, we can re-express $J_{11}$ as$$\begin{array}{}\text{(17)}& J_{11}^c=\frac{\alpha r{N}^{\ast }}{N^{\ast }+\eta }{W}^{\prime }{N}^{\ast },\end{array}$$with $W^{\prime }=\text{d}W/\text{d}N$. On the other hand, if we denote by $N_1^{\ast ,c}$ and $N_2^{\ast ,c}$ the two positive equilibria of the model (with $N_2^{\ast ,c}>N_1^{\ast ,c}$), we can deduce from the aforesaid behavior of the function $WN$ that $N_{-}^0<N_1^{\ast ,c}<N_h$ and $N_2^{\ast ,c}>N_h$. Moreover, since $WN$ is strictly increasing over $N_{-}^0,N_h$, then, in particular, $W^{\prime }{N}_1^{\ast ,c}$ is strictly positive and thus $J_{11}{N}_1^{\ast ,c}>0$. Similarly, as $WN$ is strictly decreasing for $N>N_h$, it follows that, in particular, $W^{\prime }{N}_2^{\ast ,c}$ is strictly negative and so $J_{11}{N}_2^{\ast ,c}<0$. This eventually implies that the equilibrium $E_1^c≔N_1^{\ast ,c},\eta$ is a saddle point and that $E_2^c≔N_2^{\ast ,c},\eta$ is a locally asymptotically stable node.

The above results are summed up in the following theorem.

Theorem 1.

Assume that $\mu <r$; then, system (9) possesses a boundary equilibrium $E_0^c≔0,\eta$, which is globally asymptotically stable if no positive equilibrium occurs and locally asymptotically stable if it is not the case. System (9) has either no, one, or two positive equilibria depending on whether $\gamma >{\gamma }_c^c$, $\gamma ={\gamma }_c^c$, or $\gamma <{\gamma }_c^c$, respectively, where ${\gamma }_c^c$ is defined by equation (14). In that case where two positive equilibria exist, $E_1^c≔N_1^{\ast ,c},\eta$ and $E_2^c≔N_2^{\ast ,c},\eta$ with $N_{-}^{\ast }<N_1^{\ast ,c}<N_h$ and $N_2^{\ast ,c}>N_h$, $E_1^c$ is a saddle point, and $E_2^c$ is locally asymptotically stable node.

These results are illustrated with numerical examples presented in Figure 1.

[figure omitted; refer to PDF]

4. Release Rate Proportional to the Wild Population

In general, a constant release rate of sterile insects is not optimal and better strategies can be adopted by adjusting the release rate to the size of the wild population [37–39]. One choice is to let the release rate be proportional to the number of pest insects if the latter is relatively small [25, 31]. Of course, a close monitoring of the pest population will be required, and particularly, its smallness remains critical for the economy of this choice. Within this strategy, the release function is given by $RN=\gamma N$, where $\gamma$ is a positive constant (hereafter referred to as the relative release rate). Then, system (8) takes the form$$\begin{array}{c}\text{(18)}& \frac{\text{d}N}{\text{d}t}=\alpha r\frac{N^2}{N+M}K-NN-\epsilon -\mu N,\frac{\text{d}M}{\text{d}t}=\gamma N-\mu M.\end{array}$$

Define the region of the phase plan ${\mathrm{\Omega }}_2$ as$$\begin{array}{}\text{(19)}& {\Omega }_2≔N,M\in {ℝ}^2:0\le N\le N_{+}^{\ast },\hspace{1em}0\le M\le \eta N_{+}^{\ast },\end{array}$$and then we can easily verify that ${\mathrm{\Omega }}_2$ is a positive invariant set for system (18), and we only consider $N,M\in {\mathrm{\Omega }}_2$ in this section.

The model (18) has a trivial equilibrium at the origin $\mathrm{0,0}$, and it is always locally asymptotically stable. A possible positive equilibrium verifies the quadratic equation$$\begin{array}{}\text{(20)}& K-NN-\epsilon -\tilde{\mu }+\tilde{\gamma }=0.\end{array}$$

It is worth noting that equation (20) is completely similar to the equation $QN=0$, giving the positive equilibria of the model in the absence of sterile insects, but with $\tilde{\mu }$ substituted now by ${\tilde{\mu }}_e:=\tilde{\mu }+\tilde{\gamma }$. This term highlights the effect of the sterile-fertile interactivity within the present model and may be seen as the effective death rate after sterile insects are released. Then, equation (4) has two roots:$$\begin{array}{}\text{(21)}& N_{\pm }^{\ast ,p}≔\frac{K+\epsilon \pm K-\epsilon \sqrt{1-{\tilde{\mu }}_e/r}}{2}.\end{array}$$

Thus, if ${\tilde{\mu }}_e<r$, both solutions $N_{p,\pm }^{\ast }$ are real and system (18) has two positive equilibria. But if ${\tilde{\mu }}_e>r$, that is, the effective mortality rate exceeds the maximum growth rate, no positive equilibrium occurs. This clearly defines a threshold value of the relative release rate: ${\gamma }_c^p≔r-\mu$, such that the model (18) has either no, one, or two positive equilibria if either $\gamma >{\gamma }_c^p$, $\gamma ={\gamma }_c^p$, or $\gamma <{\gamma }_c^p$, respectively.

Moreover, for $\gamma >{\gamma }_c^p$ the unique equilibrium $\mathrm{0,0}$ is globally asymptotically stable since ${\mathrm{\Omega }}_2$ is a positive invariant set for system (18). Thus, if the relative release rate of sterile insects is greater than $\gamma >{\gamma }_c^p$, the wild population ends up being eradicated whatever its initial size is.

In the case where $\gamma <{\gamma }_c^p$, the equilibria are given by $E_{-}^p≔N_{-}^{\ast ,p},\eta N_{-}^{\ast ,p}$ and $E_{+}^p≔N_{+}^{\ast ,p},\eta N_{+}^{\ast ,p}$. The stability of these equilibria is investigated below.

The general form of the Jacobian matrix at an equilibrium is$$\begin{array}{}\text{(22)}& J^p≔\begin{array}{cc}{J}_{11}^p,& J_{12}^p,\\ \gamma ,& -\mu .\end{array}.\end{array}$$

At a positive equilibrium $N^{\ast },\gamma N^{\ast }/\mu$, the matrix element $J_{11}^c$ can be written as$$\begin{array}{}\text{(23)}& J_{11}^p=\frac{\partial }{\partial N}\frac{\alpha rN}{N+M}WN-\tilde{\mu }M\frac{N^{\ast },\gamma N^{\ast }}{\mu }.\end{array}$$

Then, since $W{N}^{\ast }=\tilde{\gamma }{N}^{\ast }$, we readily find$$\begin{array}{}\text{(24)}& J_{11}^p=\frac{\alpha r}{1+\eta }\tilde{\gamma }+N^{\ast }{Q}^{\prime }{N}^{\ast }=\frac{\alpha r}{1+\eta }\tilde{\gamma }+N^{\ast }K+\epsilon -2N^{\ast }.\end{array}$$

The calculation of the element $J_{12}^p$ is also straightforward and gives$$\begin{array}{}\text{(25)}& J_{12}^p=-\frac{\mu }{1+\eta }.\end{array}$$

This leads to$$\begin{array}{}\text{(26)}& \det J^p=-\frac{\alpha r\mu }{1+\eta }{N}^{\ast }K+\epsilon -2N^{\ast },\end{array}$$wherein $N^{\ast }$ is either $N_{-}^{\ast ,p}$ or $N_{+}^{\ast ,p}$. It follows from equation (26) that$$\begin{array}{}\text{(27)}& \det J^p{N}_{\pm }^{\ast ,p}=\mp \frac{\alpha r\mu }{1+\eta }{N}_{\pm }^{\ast ,p}K-\epsilon \sqrt{\frac{1-{\tilde{\mu }}_e}{r}}.\end{array}$$

This means firstly that the equilibrium $E_{-}^p$ is a saddle point. For the equilibrium $E_{+}^p$, we need once more to see the sign of $tr{J}^p$. Then, we easily verify that$$\begin{array}{}\text{(28)}& tr{J}^p{N}_{+}^{\ast ,p}=-\frac{\alpha r\mu }{1+\eta }{N}_{+}^{\ast ,p}\tilde{\mu }+K-\epsilon \sqrt{\frac{1-{\tilde{\mu }}_e}{r}}<0.\end{array}$$

Equation (27) and equation (28) allow concluding that $E_{+}^p$ is a locally asymptotically stable node or spiral.

We sum up the results of this section in the following theorem.

Theorem 2.

The origin $\mathrm{0,0}$ is a locally asymptotically stable node for system (18). Furthermore, assume that $\mu <r$; if $\gamma >{\gamma }_c^p$ where ${\gamma }_c^p=r-\mu$, the origin is the unique equilibrium and it is globally asymptotically stable. And if $\gamma <{\gamma }_c^p$, system (9) possesses two positive equilibria, $E_{-}^p=N_{-}^{\ast ,p},\eta N_{-}^{\ast ,p}$ and $E_{+}^p=N_{+}^{\ast ,p},\eta N_{+}^{\ast ,p}$, with $N_{+}^{\ast ,p}$ given by equation (21). Moreover, $E_{-}^p$ is a saddle point and $E_{+}^p$ is a locally asymptotically stable node or spiral.

We illustrate the above results with numerical examples as given in Figure 2.

[figure omitted; refer to PDF]

5. Saturating Proportional Release Rate

As noted in the previous section, the proportional release rate may turn out to be a very costly process if the wild population becomes big sized, since the number of releases should be great as well. Then, a new strategy, compromising the two previous strategies, has been proposed in [31]. It consists in adjusting the release rate so that it is proportional to $N$ for a small $N$ but tends towards saturation for sufficiently large $N$. Following [31], in this section we consider a release function of Holling-II type: $RN=\gamma N/1+N$ with $\gamma$ a positive constant. Then, system (8) takes the form$$\begin{array}{c}\text{(29)}& \frac{\text{d}N}{\text{d}t}=\alpha r\frac{N^2}{N+M}K-NN-\epsilon -\mu N,\frac{\text{d}M}{\text{d}t}=\frac{\gamma N}{1+N}-\mu M.\end{array}$$

Then, it can be easily checked that ${\mathrm{\Omega }}_1$ defined in equation (19) is a positive invariant set for system (29), and we assume $N,M\in {\mathrm{\Omega }}_1$ in this section.

The model (29) has a trivial equilibrium at the origin $\mathrm{0,0}$, which is always locally asymptotically stable. For a positive equilibrium, the following cubic equation should be satisfied:$$\begin{array}{}\text{(30)}& N+1QN-\tilde{\gamma }=0,\end{array}$$or equivalently$$\begin{array}{}\text{(31)}& N^3-K+\epsilon -1N^2-K+\epsilon -K\epsilon -\tilde{\mu }N-K\epsilon +\tilde{\mu }+\tilde{\gamma }=0.\end{array}$$

The situation here is very similar to that of constant release rate. According to Descartes’ rule of signs, equation (31) may have either no, one, or two real positive solutions. Note that the non-existence of positive roots requires both conditions $K+ϵ<1$ and $K+ϵ-Kϵ<\tilde{\mu }$, which cannot be fulfilled simultaneously because $r>\mu$. Furthermore, a positive root occurs when the curve of the function $PN≔N+1QN$ intercepts the $\tilde{\gamma }$-axis. Furthermore, the study of the variations of $PN$ shows that $PN<0$ for $N<N_{-}^{\ast }$. Then, as $N$ increases from $N_{-}^{\ast }$, the function $PN$ rises to reach a maximum at the point$$\begin{array}{}\text{(32)}& N_h^{\prime }≔\frac{K+\epsilon -1+\sqrt{K+1\epsilon +1+{K-\epsilon }^{2}1-3\mu /4r}}{3},\end{array}$$with $N_h^{\prime }$ always real and positive since $r>\mu$. And for $r>\mu$, $PN$ is strictly decreasing. Hence, the condition for the occurrence of two positive solutions for equation (30) is$$\begin{array}{}\text{(33)}& N_h^{\prime }+1Q{N}_h^{\prime }>\tilde{\gamma }.\end{array}$$

This clearly defines a threshold value of $\gamma$ as$$\begin{array}{}\text{(34)}& {\gamma }_c^s≔\alpha r{N}_h^{\prime }+1Q{N}_h^{\prime },\end{array}$$such that equation (30) admits only one positive root if $\gamma ={\gamma }_c^s$. Therefore, we conclude that system (29) has no, one, or two positive equilibria if $\gamma >{\gamma }_c^s$, $\gamma ={\gamma }_c^s$, or $\gamma <{\gamma }_c^s$, respectively. Moreover, for $\gamma >{\gamma }_c^s$ the origin $\mathrm{0,0}$ is globally asymptotically stable as ${\mathrm{\Omega }}_1$ is a positive invariant set for system (29). For $\gamma <{\gamma }_c^s$, the additional positive equilibria of the model are given by $E_1^s≔N_1^{\ast ,s},\eta N_1^{\ast ,s}/1+N_1^{\ast ,s}$ and $E_2^s≔N_2^{\ast ,s},\eta N_2^{\ast ,s}/1+N_2^{\ast ,s}$, where $N_1^{\ast ,s}$ and $N_2^{\ast ,s}$ are positive solutions of equation (30).

Next, we investigate the stability of the equilibria $E_1^s$ and $E_2^s$. The Jacobian matrix at an equilibrium point has the form$$\begin{array}{}\text{(35)}& J^s≔\begin{array}{cc}{J}_{11}^s,& -\frac{\mu 1+N^{\ast }}{N^{\ast }+1+\eta },\\ \frac{\gamma }{{1+N^{\ast }}^2},& -\mu ,\end{array},\end{array}$$wherein $N^{\ast }$ is either $N_1^{\ast ,s}$ or $N_2^{\ast ,s}$. In addition, using the relation $W{N}^{\ast }=\tilde{\gamma }{N}^{\ast }/N^{\ast }+1$, $J_{11}^s$ can be written as$$\begin{array}{}\text{(36)}& J_{11}^s=\frac{\alpha r{N}^{\ast }+1}{N^{\ast }+1+\eta }{W}^{\prime }{N}^{\ast }.\end{array}$$

This entails$$\begin{array}{}\text{(37)}& \det J^s=-\frac{\mu \alpha r}{N^{\ast }+1+\eta }{N^{\ast }+1}^2{W}^{\prime }{N}^{\ast }-\tilde{\gamma }.\end{array}$$

On the other hand, it follows from the relation $WN=NPN/N+1$ and the fact that $P{N}^{\ast }=\tilde{\gamma }$ that$$\begin{array}{}\text{(38)}& W^{\prime }{N}^{\ast }=\frac{N^{\ast }}{N^{\ast }+1}{P}^{\prime }{N}^{\ast }+\frac{\tilde{\gamma }}{{N^{\ast }+1}^2}.\end{array}$$

This leads to$$\begin{array}{}\text{(39)}& \det J^s=-\mu \alpha r\frac{N^{\ast }{N}^{\ast }+1}{N^{\ast }+1+\eta }{P}^{\prime }{N}^{\ast }.\end{array}$$

Furthermore, from the variations of the function $PN$, we can infer that $N_{-}^0<N_1^{\ast ,s}<N_h$ and $N_2^{\ast ,s}>N_h$. Subsequently, as $PN$ is strictly increasing over $N_{-}^0,N_h$, then $P^{\prime }{N}_1^{\ast ,s}$ is certainly strictly positive and thus $\det J^s{N}_1^{\ast ,s}>0$. Similarly, as $PN$ is strictly decreasing for $N>N_h$, $P^{\prime }{N}_2^{\ast ,s}$ is consequently strictly negative and thus $\det J^s{N}_2^{\ast ,s}<0$. This implies firstly that the equilibrium $E_1^s$ is a saddle point. Next, to identify the nature of the equilibrium $E_2^s$ let us write $\text{tr}{J}^s$. We have$$\begin{array}{}\text{(40)}& tr{J}^s=\alpha r\frac{N^{\ast }+1}{N^{\ast }+1+\eta }{W}^{\prime }{N}^{\ast }-\tilde{\mu }-\frac{\tilde{\gamma }}{N^{\ast }+1}.\end{array}$$

Then, using equation (38) we find$$\begin{array}{}\text{(41)}& \text{tr}{J}^s=\alpha r\frac{N^{\ast }}{N^{\ast }+1+\eta }{P}^{\prime }{N}^{\ast }-\frac{\tilde{\mu }{N}^{\ast }+1}{N^{\ast }}-\frac{\tilde{\gamma }}{N^{\ast }+1},\end{array}$$

And since $P^{\prime }{N}_2^{\ast ,s}<0$, we conclude that $E_2^s$ is a locally asymptotically stable node.

The results of this section are summed up in the following theorem.

Theorem 3.

System (29) has a locally asymptotically stable node at the origin $\mathrm{0,0}$. Moreover, assume that $\mu <r$, if $\gamma >{\gamma }_c^s$ where ${\gamma }_c^s$ is given in equation (34), the origin is the unique equilibrium and it is globally asymptotically stable. And if $\gamma <{\gamma }_c^s$ system two positive equilibria occur for (29): $E_1^s≔N_1^{\ast ,s},\eta N_1^{\ast ,s}/1+N_1^{\ast ,s}$ and $E_2^s≔N_2^{\ast ,s},\eta N_2^{\ast ,s}/1+N_2^{\ast ,s}$, with $N_{-}^{\ast }<N_1^{\ast ,s}<N_h$ and $N_2^{\ast ,s}>N_h$. In addition, $E_1^s$ is a saddle point and $E_2^s$ is a locally asymptotically stable node.

The results for this model are illustrated with numerical examples shown in Figure 3.

[figure omitted; refer to PDF]

6. Conclusion

In summary, this work studies relatively simple mathematical models describing the dynamics of interactive wild and sterile insect populations, occurring within the SIT. The latter is a method of biological control, in which sterile males are released to reduce or eradicate a pest population, which can effectively help contain the spread of many pest insects such as the Red Palm Weevil (RPW). Modelling assumptions adopted in this paper allow for substantial simplifications of the SIT dynamics but in the meantime yield a model that can bias the results in comparison with real biological situations. As a matter of fact, we assumed homogeneous insect populations such that no male-female or stage distinction has been made. Here, we note that a recent comparison made between stage structured models and homogeneous models revealed that they share very similar dynamical features [25]. Moreover, death rates for sterile and fertile insects are assumed density-independent and equal. This seems reasonable since released insects are all adults so that the competition between them for natural resources is relatively weak. On the other hand, the competition for mates is the mechanism emphasised by SIT. Such competition obviously does not affect the sterile population. The sterile-fertile interaction is then assumed to be a one-sided competition that affects only the wild-type population. Moreover, to account for the need of a critical threshold density in order that the wild population could persist, a strong Allee effect has been included in the growth term of the wild population. Subsequently, we considered three submodels, each characterized by a different strategy of release: the first involves a constant release rate, the second assumes a release rate proportional to the size of the wild population, and the third uses a release rate of Holling-II type. We have carried out complete mathematical analysis of these submodels and discussed their dynamical features, especially the existence of equilibria and their stability. In particular, we demonstrated the existence of release threshold for all the three strategies. Thus, if the release rate is below the threshold value, each submodel admits two positive interior equilibria: a saddle point and a stable node (or spiral for the second submodel). However, as soon as the release rate exceeds the threshold value, positive equilibria no longer exist. Each of the three submodels possesses a unique equilibrium with a vanishing number of wild insects. In this situation, the wild population evolves to extinction whatever its initial number is. Finally, our analytical findings for all submodels have been illustrated with numerical examples.

Authors’ Contributions

The study was carried out in collaboration among all authors. All authors read and approved the final manuscript.

Acknowledgments

The third author would like to thank his Professors/Scientists Prof. Mohamed Haiour, Prof. Ahmed-Salah Chibi, and Prof. Azzedine Benchettah at Annaba University in Algeria for the important content of Master’s and PhD courses in pure and applied mathematics which he received during his studies. Moreover, he thanks them for the additional help they provided to him during office hours in their office about the few concepts/difficulties he had encountered, and he appreciates their talent and dedication for their postgraduate students currently and previously. In addition, the authors gratefully acknowledge Qassim University, represented by the Deanship of Scientific Research, on the material support for this research under the number 6766-alrasscac-bs-2019-2-2-I during the academic year 1441 AH/2019 AD.

References

[1] E. F. Knipling, "Possibilities of insect control or eradication through the use of sexually sterile Males1," Journal of Economic Entomology, vol. 48 no. 4, pp. 459-462, DOI: 10.1093/jee/48.4.459, 1955.

[2] V. A. Dyck, J. Hendrichs, A. S. Robinson, The Sterile Insect Technique, Principles and Practice in Area-wide Integrated Pest Management, 2006.

[3] K. Bourtzis, Wolbachia-based Technologies for Insect Pest Population Control: Advances in Experimental Medicine and Biology, 2008.

[4] S. P. Sinkins, "Wolbachia and cytoplasmic incompatibility in mosquitoes," Insect Biochemistry and Molecular Biology, vol. 34 no. 7, pp. 723-729, DOI: 10.1016/j.ibmb.2004.03.025, 2004.

[5] S. Zabalou, M. Riegler, M. Theodorakopoulou, C. Savakis, K. Bourtzis, "Wolbachia-induced cytoplasmic incompatibility as a means for insect pest population control," Proceedings of the National Academy of Sciences, vol. 101 no. 42, pp. 15042-15045, DOI: 10.1073/pnas.0403853101, 2004.

[6] M. Marranzano, R. Ragusa, M. Platania, G. Faro, M. A. Coniglio, "Knowledge, attitudes and practices towards patients with HIV/AIDS in staff Nurses in one university hospital in Sicily," Ital. Journal of Public Health, vol. 10 no. 1, 2013.

[7] K. Nikolouli, F. Sassù, L. Mouton, C. Stauffer, K. Bourtzis, "Combining sterile and incompatible insect techniques for the population suppression of Drosophila suzukii," Journal of Pest Science, vol. 93 no. 2, pp. 647-661, DOI: 10.1007/s10340-020-01199-6, 2020.

[8] R. Ragusa, S. Russo, L. Villari, G. Schilirò, "Hodgkin’s disease as a second malignant neoplasm in childhood: report of a case and review of the literature," Pediatric Hematology and Oncology, vol. 18 no. 6, pp. 407-414, DOI: 10.1080/088800101316922038, 2001.

[9] G. Reale, G. I. Russo, M. Di Mauro, "Association between dietary flavonoids intake and prostate cancer risk: a case-control study in Sicily," Complementary Therapies in Medicine, vol. 39, pp. 14-18, DOI: 10.1016/j.ctim.2018.05.002, 2018.

[10] D. Stauffer, R. S. Lees, Z. Xi, "Combining the sterile insect technique with Wolbachia-based approaches: II a safer approach to Aedes albopictus population suppression programmes, designed to minimize the consequences of inadvertent female release," PLoS One, vol. 10 no. 8,DOI: 10.1371/journal.pone.0135194, 2015.

[11] D. Zhang, X. Zheng, Z. Xi, K. Bourtzis, J. R. Gilles, "Combining the sterile insect technique with the incompatible insect technique: I-impact of wolbachia infection on the ftness of triple- and double-infected strains of Aedes albopictus," PLoS One, vol. 10 no. 4,DOI: 10.1371/journal.pone.0121126, 2015a.

[12] P.-S. Yen, A.-B. Failloux, "A review: wolbachia-based population replacement for mosquito control shares common points with genetically modified control approaches," Pathogens, vol. 9 no. 5, pp. 404-421, DOI: 10.3390/pathogens9050404, 2020.

[13] J. L. Rasgon, T. W. Scott, "Wolbachia and cytoplasmic incompatibility in the California culex pipiens mosquito species complex: parameter estimates and infection dynamics in natural populations," Genetics, vol. 165 no. 4, pp. 2029-2038, 2003.

[14] A. A. Hoffmann, B. L. Montgomery, J. Popovici, "Successful establishment of wolbachia in aedes populations to suppress dengue transmission," Nature, vol. 476 no. 7361, pp. 454-457, DOI: 10.1038/nature10356, 2011.

[15] S. Boulaaras, M. Haiour, "The finite element approximation of evolutionary Hamilton-Jacobi-Bellman equations with nonlinear source terms," Indagationes Mathematicae, vol. 24 no. 1, pp. 161-173, DOI: 10.1016/j.indag.2012.07.005, 2013.

[16] S. Boulaaras, M. Haiour, "L∞-asymptotic behavior for a finite element approximation in parabolic quasi-variational inequalities related to impulse control problem," Applied Mathematics and Computation, vol. 217 no. 14, pp. 6443-6450, DOI: 10.1016/j.amc.2011.01.025, 2011.

[17] S. Boulaaras, "Some existence results for a new class of elliptic kirchhoff equation with logarithmic source terms," Journal of Intelligent & Fuzzy Systems, vol. 37 no. 6, pp. 8335-8344, DOI: 10.3233/jifs-190885, 2019.

[18] S. Boulaaras, "Some new properties of asynchronous algorithms of theta scheme combined with finite elements methods for an evolutionary implicit 2-sided obstacle problem," Mathematical Methods in the Applied Sciences, vol. 40 no. 18, pp. 7231-7239, DOI: 10.1002/mma.4525, 2017.

[19] S. Boulaaras, "An optimal error estimate of finite element method for parabolic quasi-variational inequalities with non linear source terms," Asymptotic Analysis, vol. 100 no. 3-4, pp. 193-208, DOI: 10.3233/asy-161392, 2016.

[20] S. Boulaaras, "Existence of positive solutions for a new class of Kirchhoff parabolic systems," Rocky Mountain Journal of Mathematics, vol. 50 no. 2, pp. 445-454, DOI: 10.1216/rmj.2020.50.445, 2020.

[21] C. Dufourd, Y. Dumont, "Modeling and simulations of mosquito dispersal. the case of Aedes albopictus," Biomath, vol. 1 no. 2, 2012.

[22] C. Dufourd, Y. Dumont, "Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control," Computers & Mathematics with Applications, vol. 66 no. 9, pp. 1695-1715, DOI: 10.1016/j.camwa.2013.03.024, 2013.

[23] Y. Dumont, J. M. Tchuenche, "Mathematical studies on the sterile insect technique for the chikungunya disease and Aedes albopictus," Journal of Mathematical Biology, vol. 65 no. 5, pp. 809-855, DOI: 10.1007/s00285-011-0477-6, 2012.

[24] M. Huang, X. Song, J. Li, "Modelling and analysis of impulsive releases of sterile mosquitoes," Journal of Biological Dynamics, vol. 11 no. 1, pp. 147-171, DOI: 10.1080/17513758.2016.1254286, 2017.

[25] J. Li, L. Cai, Y. Li, "Stage-structured wild and sterile mosquito population models and their dynamics," Journal of Biological Dynamics, vol. 11 no. 1, pp. 79-101, DOI: 10.1080/17513758.2016.1159740, 2017.

[26] J. Li, Z. Yuan, "Modelling releases of sterile mosquitoes with different strategies," Journal of Biological Dynamics, vol. 9 no. 1,DOI: 10.1080/17513758.2014.977971, 2015.

[27] J. Z. Farkas, P. Hinow, "Structured and unstructured continuous models for wolbachia infections," Bulletin of Mathematical Biology, vol. 72 no. 8, pp. 2067-2088, DOI: 10.1007/s11538-010-9528-1, 2010.

[28] A. Fenton, K. N. Johnson, J. C. Brownlie, G. D. D. Hurst, "Solving the wolbachia paradox: modeling the tripartite interaction between host, wolbachia, and a natural enemy," The American Naturalist, vol. 178 no. 3, pp. 333-342, DOI: 10.1086/661247, 2011.

[29] J. Z. Farkas, S. A. Gourley, R. Liu, A.-A. Yakubu, "Modelling Wolbachia infection in a sex-structured mosquito population carrying west nile virus," Journal of Mathematical Biology, vol. 75 no. 3, pp. 621-647, DOI: 10.1007/s00285-017-1096-7, 2017.

[30] G. Nadin, M. Strugarek, N. Vauchelet, "Hindrances to bistable front propagation: application to wolbachia invasion," Journal of Mathematical Biology, vol. 76 no. 6, pp. 1489-1533, DOI: 10.1007/s00285-017-1181-y, 2018.

[31] A. Bendhahbi, Y. Chargui, S. Boulaaras, S. Ben Khalifa, "A one-sided competition mathematical model for the sterile insect technique," Complexity, vol. 2020,DOI: 10.1155/2020/6246808, 2020.

[32] J. Li, "New revised simple models for interactive wild and sterile mosquito populations and their dynamics," Journal of Biological Dynamics, vol. 11 no. 2, pp. 316-333, DOI: 10.1080/17513758.2016.1216613, 2017.

[33] F. Courchamp, T. Clutton-Brock, B. Grenfell, "Inverse density dependence and the Allee effect," Trends in Ecology & Evolution, vol. 14 no. 10, pp. 405-410, DOI: 10.1016/s0169-5347(99)01683-3, 1999.

[34] Y. Gruntfest, R. Arditi, Y. Dombrovsky, "A fragmented population in a varying environment," Journal of Theoretical Biology, vol. 185 no. 4, pp. 539-547, DOI: 10.1006/jtbi.1996.0358, 1997.

[35] J. A. Walter, K. L. Grayson, D. M. Johnson, "Variation in allee effects: evidence, unknowns, and directions forward," Population Ecology, vol. 59 no. 2, pp. 99-107, DOI: 10.1007/s10144-017-0576-3, 2017.

[36] F. Courchamp, L. Berec, J. Gascoigne, Allee Effects in Ecology and Conservation, 2008.

[37] C. P. Ferreira, H. M. Yang, L. Esteva, "Assessing the suitability of sterile insect technique applied toaedes aegypti," Journal of Biological Systems, vol. 16 no. 04, pp. 565-577, DOI: 10.1142/s0218339008002691, 2008.

[38] J. Kim, H. Lee, C. Lee, S. Lee, "Assessment of optimal strategies in a two-patch dengue transmission model with seasonality," PLoS One, vol. 12 no. 3,DOI: 10.1371/journal.pone.0173673, 2017.

[39] L. Multerer, T. Smith, N. Chitnis, "Modeling the impact of sterile males on an Aedes aegypti population with optimal control," Mathematical Biosciences, vol. 311, pp. 91-102, DOI: 10.1016/j.mbs.2019.03.003, 2019.