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.

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 [35]. 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) [611]. 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, [1526] for SIT/IIT and [2730] 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 be the number of wild insects at time . In the absence of sterile insects, is assumed to evolve according to the following equation:where is the growth function of the wild population and 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]where is the carrying capacity of the niche and is the Allee effect parameter , and we have set , so that the parameter 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].

This equation has a trivial equilibrium point at , which is locally asymptotically stable. A positive equilibrium of equation (3) should be a zero of the function defined bywith . This function has two roots:which are real under the condition , 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 , , or , respectively.

In the case when the condition is fulfilled, at a positive equilibrium we have

Thus,which means that the equilibrium is unstable, whereas the equilibrium is asymptotically stable. From now, we assume that the condition 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 the number of sterile insects at time , the interactive dynamics are then governed by the following system:

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

3. Constant Release Rate

We consider here the situation where sterile insects are constantly released so that , with , a positive constant. Then, system (8) takes the form

For this model, it is easy to check that the rectangle of the phase planewith , is a positive invariant, and we assume in this section that .

System (9) has a first equilibrium lying on the boundary of . 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:with . 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 intercepts the -axis once or twice, respectively. Furthermore, an elementary study of the variations of shows that for . Then, as increases from , this function also increases from zero to attain a certain maximum at the pointbefore lowering indefinitely. Note here that is always real and positive as we have presumed that . We then deduce that the condition for the existence of two positive solutions for equation (11) is

This defines a threshold value of the release rate asat which equation (11) admits a unique positive root. Note that this definition of tacitly requires that be positive, and we can easily show that this is still the case as .

Thus, we infer that the model given by equation (9) will have no, one, or two positive equilibria if , , or , respectively. Note also that all the equilibria of the model lie on the border of the rectangle at which .

Next, we address the stability of the equilibria. The Jacobian matrix at an equilibrium point has the form

At the boundary equilibrium , we have ; therefore, is always asymptotically stable. Moreover, in the case where no positive equilibrium exists as is a positive invariant set, we conclude by the Poincaré–Bendixson theorem that is globally asymptotically stable. Thus, for a release rate of sterile insects , 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 , the matrix element can be written as

Then, keeping in mind that , we can re-express aswith . On the other hand, if we denote by and the two positive equilibria of the model (with ), we can deduce from the aforesaid behavior of the function that and . Moreover, since is strictly increasing over , then, in particular, is strictly positive and thus . Similarly, as is strictly decreasing for , it follows that, in particular, is strictly negative and so . This eventually implies that the equilibrium is a saddle point and that is a locally asymptotically stable node.

The above results are summed up in the following theorem.

Theorem 1. Assume that ; then, system (9) possesses a boundary equilibrium , 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 , , or , respectively, where is defined by equation (14). In that case where two positive equilibria exist, and with and , is a saddle point, and is locally asymptotically stable node.

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

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 [3739]. 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 , where is a positive constant (hereafter referred to as the relative release rate). Then, system (8) takes the form

Define the region of the phase plan asand then we can easily verify that is a positive invariant set for system (18), and we only consider in this section.

The model (18) has a trivial equilibrium at the origin , and it is always locally asymptotically stable. A possible positive equilibrium verifies the quadratic equation

It is worth noting that equation (20) is completely similar to the equation , giving the positive equilibria of the model in the absence of sterile insects, but with substituted now by . 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:

Thus, if , both solutions are real and system (18) has two positive equilibria. But if , 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: , such that the model (18) has either no, one, or two positive equilibria if either , , or , respectively.

Moreover, for the unique equilibrium is globally asymptotically stable since is a positive invariant set for system (18). Thus, if the relative release rate of sterile insects is greater than , the wild population ends up being eradicated whatever its initial size is.

In the case where , the equilibria are given by and . The stability of these equilibria is investigated below.

The general form of the Jacobian matrix at an equilibrium is

At a positive equilibrium , the matrix element can be written as

Then, since , we readily find

The calculation of the element is also straightforward and gives

This leads towherein is either or . It follows from equation (26) that

This means firstly that the equilibrium is a saddle point. For the equilibrium , we need once more to see the sign of . Then, we easily verify that

Equation (27) and equation (28) allow concluding that is a locally asymptotically stable node or spiral.

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

Theorem 2. The origin is a locally asymptotically stable node for system (18). Furthermore, assume that ; if where , the origin is the unique equilibrium and it is globally asymptotically stable. And if , system (9) possesses two positive equilibria, and , with given by equation (21). Moreover, is a saddle point and is a locally asymptotically stable node or spiral.

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

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 for a small but tends towards saturation for sufficiently large . Following [31], in this section we consider a release function of Holling-II type: with a positive constant. Then, system (8) takes the form

Then, it can be easily checked that defined in equation (19) is a positive invariant set for system (29), and we assume in this section.

The model (29) has a trivial equilibrium at the origin , which is always locally asymptotically stable. For a positive equilibrium, the following cubic equation should be satisfied:or equivalently

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 and , which cannot be fulfilled simultaneously because . Furthermore, a positive root occurs when the curve of the function intercepts the -axis. Furthermore, the study of the variations of shows that for . Then, as increases from , the function rises to reach a maximum at the pointwith always real and positive since . And for , is strictly decreasing. Hence, the condition for the occurrence of two positive solutions for equation (30) is

This clearly defines a threshold value of assuch that equation (30) admits only one positive root if . Therefore, we conclude that system (29) has no, one, or two positive equilibria if , , or , respectively. Moreover, for the origin is globally asymptotically stable as is a positive invariant set for system (29). For , the additional positive equilibria of the model are given by and , where and are positive solutions of equation (30).

Next, we investigate the stability of the equilibria and . The Jacobian matrix at an equilibrium point has the formwherein is either or . In addition, using the relation , can be written as

This entails

On the other hand, it follows from the relation and the fact that that

This leads to

Furthermore, from the variations of the function , we can infer that and . Subsequently, as is strictly increasing over , then is certainly strictly positive and thus . Similarly, as is strictly decreasing for , is consequently strictly negative and thus . This implies firstly that the equilibrium is a saddle point. Next, to identify the nature of the equilibrium let us write . We have

Then, using equation (38) we find

And since , we conclude that 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 . Moreover, assume that , if where is given in equation (34), the origin is the unique equilibrium and it is globally asymptotically stable. And if system two positive equilibria occur for (29): and , with and . In addition, is a saddle point and is a locally asymptotically stable node.

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

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.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

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.