# Mathematical Modelling of the Sterile Insect Technique Using Different Release Strategies

Authors: Anis Ben Dhahbi; Yassine Chargui; Salah Mahmoud Boulaaras; Sana Ben Khalifa; Waleed Koko; Faisal Alresheedi
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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, theWolbachia 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 theWolbachia-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 Section2, 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.LetNt 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:(1)dNdt=GN−μN,t≥0,where GN 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](2)GN=αrNK−NN−ε,where K is the carrying capacity of the niche and ϵ is the Allee effect parameter ϵ<K, and we have set α=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].(3)dNdt=αrNK−NN−ε−μN.This equation has a trivial equilibrium point atN0∗=0, which is locally asymptotically stable. A positive equilibrium of equation (3) should be a zero of the function QN defined by(4)QN≔K−NN−ε−μ˜,with μ˜=μ/αr. This function has two roots:(5)N±∗=K+ε±K−ε1−μ/r2,which are real under the condition r>μ, 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>μ, μ=r, or r<μ, respectively.In the case when the conditionr>μ is fulfilled, at a positive equilibrium N∗ we have(6)ddNαrNK−NN−ε−μNN∗=αrN∗K+ε−2N∗.Thus,(7)ddNαrNK−NN−ε−μNN±∗=∓N±∗K−ε1−μr,which means that the equilibrium N−∗ is unstable, whereas the equilibrium N+∗ is asymptotically stable. From now, we assume that the condition r>μ 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 byMt the number of sterile insects at time t, the interactive dynamics are then governed by the following system:(8)dNdt=αrN2N+MK−NN−ε−μN,dMdt=RN−μM,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>μ.


## 3. Constant Release Rate

We consider here the situation where sterile insects are constantly released so thatRN=γ, with γ, a positive constant. Then, system (8) takes the form(9)dNdt=αrN2N+MK−NN−ε−μN,dMdt=γ−μM.For this model, it is easy to check that the rectangle of the phase plane(10)Ω1≔N,M∈ℝ2:0≤N≤N+∗,0≤M≤η,with η≔γ/μ, is a positive invariant, and we assume in this section that N,M∈Ω1.System (9) has a first equilibrium E0c≔0,η lying on the boundary of Ω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:(11)NQN−γ˜=0,with γ˜=γ/α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 γ˜-axis once or twice, respectively. Furthermore, an elementary study of the variations of WN shows that WN<0 for N<N−∗. Then, as N increases from N−∗, this function also increases from zero to attain a certain maximum at the point(12)Nh≔K+ε+Kε+K−ε21−3μ/4r3,before lowering indefinitely. Note here that Nh is always real and positive as we have presumed that r>μ. We then deduce that the condition for the existence of two positive solutions for equation (11) is(13)NhQNh>γ˜.This defines a threshold value of the release rateγ as(14)γcc≔αrNhQNh,at which equation (11) admits a unique positive root. Note that this definition of γcc tacitly requires that QNh be positive, and we can easily show that this is still the case as r>μ.Thus, we infer that the model given by equation (9) will have no, one, or two positive equilibria if γ>γcc, γ=γcc, or γ<γcc, respectively. Note also that all the equilibria of the model lie on the border of the rectangle Ω1 at which M=η.Next, we address the stability of the equilibria. The Jacobian matrix at an equilibrium pointN∗,η has the form(15)Jc≔J11c,−μN∗N∗+η,0,−μ.At the boundary equilibriumE0c, we have J11c=−μ; therefore, E0c is always asymptotically stable. Moreover, in the case where no positive equilibrium exists as Ω1 is a positive invariant set, we conclude by the Poincaré–Bendixson theorem that E0c is globally asymptotically stable. Thus, for a release rate of sterile insects γ>γcc, 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 equilibriumN∗,η, the matrix element J11c can be written as(16)J11c=ddNαrNN+ηWN−γ˜N∗.Then, keeping in mind thatWN∗=γ˜, we can re-express J11 as(17)J11c=αrN∗N∗+ηW′N∗,with W′=dW/dN. On the other hand, if we denote by N1∗,c and N2∗,c the two positive equilibria of the model (with N2∗,c>N1∗,c), we can deduce from the aforesaid behavior of the function WN that N−0<N1∗,c<Nh and N2∗,c>Nh. Moreover, since WN is strictly increasing over N−0,Nh, then, in particular, W′N1∗,c is strictly positive and thus J11N1∗,c>0. Similarly, as WN is strictly decreasing for N>Nh, it follows that, in particular, W′N2∗,c is strictly negative and so J11N2∗,c<0. This eventually implies that the equilibrium E1c≔N1∗,c,η is a saddle point and that E2c≔N2∗,c,η is a locally asymptotically stable node.The above results are summed up in the following theorem.Theorem 1.
Assume thatμ<r; then, system (9) possesses a boundary equilibrium E0c≔0,η, 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 γ>γcc, γ=γcc, or γ<γcc, respectively, where γcc is defined by equation (14). In that case where two positive equilibria exist, E1c≔N1∗,c,η and E2c≔N2∗,c,η with N−∗<N1∗,c<Nh and N2∗,c>Nh, E1c is a saddle point, and E2c is locally asymptotically stable node.These results are illustrated with numerical examples presented in Figure1.Figure 1
Trajectories and equilibria of system (9) for parameters values chosen as K=14.7, ϵ=4.3, r=1.7, and μ=1.1. The release threshold is γcc=5.843. For γ=6.5>γcc, the unique equilibrium is E0c=0,5.91 and it is globally stable (all solutions approach this point). For γ=3.5<γcc, two positive equilibria occur: E1c=7.915,3.182 and E2c=11.686,5.909, with E1c a saddle point and E2c a stable node.


## 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=γN, where γ is a positive constant (hereafter referred to as the relative release rate). Then, system (8) takes the form(18)dNdt=αrN2N+MK−NN−ε−μN,dMdt=γN−μM.Define the region of the phase planΩ2 as(19)Ω2≔N,M∈ℝ2:0≤N≤N+∗,0≤M≤ηN+∗,and then we can easily verify that Ω2 is a positive invariant set for system (18), and we only consider N,M∈Ω2 in this section.The model (18) has a trivial equilibrium at the origin 0,0, and it is always locally asymptotically stable. A possible positive equilibrium verifies the quadratic equation(20)K−NN−ε−μ˜+γ˜=0.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 μ˜ substituted now by μ˜e:=μ˜+γ˜. 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:(21)N±∗,p≔K+ε±K−ε1−μ˜e/r2.Thus, ifμ˜e<r, both solutions Np,±∗ are real and system (18) has two positive equilibria. But if μ˜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: γcp≔r−μ, such that the model (18) has either no, one, or two positive equilibria if either γ>γcp, γ=γcp, or γ<γcp, respectively.Moreover, forγ>γcp the unique equilibrium 0,0 is globally asymptotically stable since Ω2 is a positive invariant set for system (18). Thus, if the relative release rate of sterile insects is greater than γ>γcp, the wild population ends up being eradicated whatever its initial size is.In the case whereγ<γcp, the equilibria are given by E−p≔N−∗,p,ηN−∗,p and E+p≔N+∗,p,ηN+∗,p. The stability of these equilibria is investigated below.The general form of the Jacobian matrix at an equilibrium is(22)Jp≔J11p,J12p,γ,−μ..At a positive equilibriumN∗,γN∗/μ, the matrix element J11c can be written as(23)J11p=∂∂NαrNN+MWN−μ˜MN∗,γN∗μ.Then, sinceWN∗=γ˜N∗, we readily find(24)J11p=αr1+ηγ˜+N∗Q′N∗=αr1+ηγ˜+N∗K+ε−2N∗.The calculation of the elementJ12p is also straightforward and gives(25)J12p=−μ1+η.This leads to(26)detJp=−αrμ1+ηN∗K+ε−2N∗,wherein N∗ is either N−∗,p or N+∗,p. It follows from equation (26) that(27)detJpN±∗,p=∓αrμ1+ηN±∗,pK−ε1−μ˜er.This means firstly that the equilibriumE−p is a saddle point. For the equilibrium E+p, we need once more to see the sign of trJp. Then, we easily verify that(28)trJpN+∗,p=−αrμ1+ηN+∗,pμ˜+K−ε1−μ˜er<0.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 origin0,0 is a locally asymptotically stable node for system (18). Furthermore, assume that μ<r; if γ>γcp where γcp=r−μ, the origin is the unique equilibrium and it is globally asymptotically stable. And if γ<γcp, system (9) possesses two positive equilibria, E−p=N−∗,p,ηN−∗,p and E+p=N+∗,p,ηN+∗,p, with N+∗,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 Figure2.Figure 2
Trajectories and equilibria of system (18) for parameters values chosen as K=17.3, ϵ=4.9, r=2.7, and μ=0.9. The release threshold is γcc=1.9. For γ=3.2>γcc, the unique equilibrium is the origin 0,0 and it is globally stable (all solutions approach the origin). For γ=1.1<γcc, there are two positive equilibria: E−p=7.943,9.708 and E+p=14.257,17.425, where E−p is a saddle point and E+p is a stable node.


## 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=γN/1+N with γ a positive constant. Then, system (8) takes the form(29)dNdt=αrN2N+MK−NN−ε−μN,dMdt=γN1+N−μM.Then, it can be easily checked thatΩ1 defined in equation (19) is a positive invariant set for system (29), and we assume N,M∈Ω1 in this section.The model (29) has a trivial equilibrium at the origin 0,0, which is always locally asymptotically stable. For a positive equilibrium, the following cubic equation should be satisfied:(30)N+1QN−γ˜=0,or equivalently(31)N3−K+ε−1N2−K+ε−Kε−μ˜N−Kε+μ˜+γ˜=0.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ϵ<μ˜, which cannot be fulfilled simultaneously because r>μ. Furthermore, a positive root occurs when the curve of the function PN≔N+1QN intercepts the γ˜-axis. Furthermore, the study of the variations of PN shows that PN<0 for N<N−∗. Then, as N increases from N−∗, the function PN rises to reach a maximum at the point(32)Nh′≔K+ε−1+K+1ε+1+K−ε21−3μ/4r3,with Nh′ always real and positive since r>μ. And for r>μ, PN is strictly decreasing. Hence, the condition for the occurrence of two positive solutions for equation (30) is(33)Nh′+1QNh′>γ˜.This clearly defines a threshold value ofγ as(34)γcs≔αrNh′+1QNh′,such that equation (30) admits only one positive root if γ=γcs. Therefore, we conclude that system (29) has no, one, or two positive equilibria if γ>γcs, γ=γcs, or γ<γcs, respectively. Moreover, for γ>γcs the origin 0,0 is globally asymptotically stable as Ω1 is a positive invariant set for system (29). For γ<γcs, the additional positive equilibria of the model are given by E1s≔N1∗,s,ηN1∗,s/1+N1∗,s and E2s≔N2∗,s,ηN2∗,s/1+N2∗,s, where N1∗,s and N2∗,s are positive solutions of equation (30).Next, we investigate the stability of the equilibriaE1s and E2s. The Jacobian matrix at an equilibrium point has the form(35)Js≔J11s,−μ1+N∗N∗+1+η,γ1+N∗2,−μ,,wherein N∗ is either N1∗,s or N2∗,s. In addition, using the relation WN∗=γ˜N∗/N∗+1, J11s can be written as(36)J11s=αrN∗+1N∗+1+ηW′N∗.This entails(37)detJs=−μαrN∗+1+ηN∗+12W′N∗−γ˜.On the other hand, it follows from the relationWN=NPN/N+1 and the fact that PN∗=γ˜ that(38)W′N∗=N∗N∗+1P′N∗+γ˜N∗+12.This leads to(39)detJs=−μαrN∗N∗+1N∗+1+ηP′N∗.Furthermore, from the variations of the functionPN, we can infer that N−0<N1∗,s<Nh and N2∗,s>Nh. Subsequently, as PN is strictly increasing over N−0,Nh, then P′N1∗,s is certainly strictly positive and thus detJsN1∗,s>0. Similarly, as PN is strictly decreasing for N>Nh, P′N2∗,s is consequently strictly negative and thus detJsN2∗,s<0. This implies firstly that the equilibrium E1s is a saddle point. Next, to identify the nature of the equilibrium E2s let us write trJs. We have(40)trJs=αrN∗+1N∗+1+ηW′N∗−μ˜−γ˜N∗+1.Then, using equation (38) we find(41)trJs=αrN∗N∗+1+ηP′N∗−μ˜N∗+1N∗−γ˜N∗+1,And sinceP′N2∗,s<0, we conclude that E2s 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 0,0. Moreover, assume that μ<r, if γ>γcs where γcs is given in equation (34), the origin is the unique equilibrium and it is globally asymptotically stable. And if γ<γcs system two positive equilibria occur for (29): E1s≔N1∗,s,ηN1∗,s/1+N1∗,s and E2s≔N2∗,s,ηN2∗,s/1+N2∗,s, with N−∗<N1∗,s<Nh and N2∗,s>Nh. In addition, E1s is a saddle point and E2s is a locally asymptotically stable node.The results for this model are illustrated with numerical examples shown in Figure3.Figure 3
Trajectories and equilibria of system (29) for parameters values chosen as K=16.7, ϵ=3.7, r=1.3, and μ=0.7. The release threshold is γcs=6.862. For γ=8>γcs, the origin E0c=0,0 is the unique equilibrium and it is globally stable (all solutions approach this point). For γ=5<γcs, two positive equilibria occur: E1s=8.441,6.386 and E2c=12.877,6.628, with E1s a saddle point and E2s a stable node.


## 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.



---

*Source: 8896566-2020-11-07.xml*