OPTIMAL STERILE INSECT RELEASE FOR AREA-WIDE INTEGRATED PEST MANAGEMENT IN A DENSITY REGULATED PEST POPULATION

. To determine optimal sterile insect release policies in area-wide integrated pest management is a challenge that users of this pest control method inevitably confront. In this note we provide approximations to best policies of release through the use of simulated annealing. The discrete time model for the population dynamics includes the eﬀects of sterile insect release and density dependence in the pest population. Spatial movement is introduced through integrodiﬀerence equations, which allow the use of the stochastic search in cases where movement is described through arbitrary dispersal kernels. As a byproduct of the computations, an assessment of appropriate control zone sizes is possible.

1. Introduction.The sterile insect technique (SIT) is used worldwide for control and eradication of insect pest invasions.An innovative application of radioactivity to improve agriculture, was successfully implemented by Edward J. Knipling in the late 30's, [14,17].As originally conceived, the goal is to induce the production of males' damaged gametes through a massive, but controlled, exposure of insects to gamma radiation.After irradiation, males are released in the fields, where they appear indistinguishable to females.As a consequence, sterile males' impaired reproductive cells cause inhibition of mitosis and death of fertilized eggs or embryos [1].For many years it has been observed how the appropriate and repeated application of this technique eventually drives pest populations to extinction, or supresses them to small sizes susceptible of extermination by alternative methods [14].The technique artificially induces female fertilization disruption, creating an Allee effect [9,16] that might drive pest extinction.The understanding of some practical aspects of SIT has continuously benefited from an increasing development of mathematical modeling elaborated during the last decades, which systematically included biological details with potential impact in the success or failure of control programs.Deterministic [5,6], stochastic [8] and computational [27] models have been considered with the aim to predict or explain pest population behavior under constrains, like incomplete sterility [4], decreased competitivity in sterile males [2] or the aggregation of insects [3] among other relevant features.

LUIS F. GORDILLO
Extensive evidence about SIT efficiency eliminating or suppressing pests and its clear benefits over the use of insecticides (it does not contaminate the environment, does not target other insect populations and is safe for farmers health) situates the technique among one of the most preferred, although its implementation can be excessively expensive.The technique was used to eradicate the New World screwworm (Cochliomyia hominivorax ) in Florida in 1958, followed by successful programmes covering the Southwest of the United States, Mexico and Central America [14].SIT has been also applied worldwide to control Tephritid fruit flies [11] and Tsetse flies in Africa [14].However, it might confront several issues when applied against lepidopteran or coleopteran species, including a more difficult (and costly) large-scale rearing of insects than that of flies [14].
In general, one of the main concerns in current eradication programs that use SIT is to determine optimal sterile release policies that guarantee economic and technical feasibility, a task that frequently faces difficulties.The lack of appropriate modeling for area-wide integrated pest management might leave most of the decisions to be taken only by guessing.As a consequence, some projects might be abandoned due to economic constraints and, in other cases, the use of SIT simply fails due to a poor design.Recent efforts to find adequate models to support control policies integrating area-wide integrated pest management have been made in [7], where reaction-difussion approximations are suggested to compute appropriate size of control areas.
In this paper we show how to compute optimal SIT application in spatial contexts by combining a discrete time model, which includes density regulation on the larval survival to the adult stage [22], and dispersal using integrodifference equations [16].This allows us to approximate optimal policies through the use of stochastic search optimization.The approach is general and admits the treatment of cases where pest movement is associated with different types of dispersal.
The paper is organized as follows.In Section 2 we briefly review basic models for SIT release in discrete time.In Section 3 dispersal is included and the effect of sterile release on traveling waves is observed.The case of variable spatial sterile release is considered in Section 4, where approximate optimal sterile release policies are computed in an example.
2. Knipling's paradigm in discrete time.In this section basic ideas of SIT mathematical modeling are briefly reviewed without consideration of spatial effects.The reader is referred to [5] for an extensive account.Consider an insect population that has non overlapping generations.Let F t and M t denote the female and fertile male densities in the population at generation t (we assume all females are fertile), and λ the female's rate of increase per generation.Assuming that the sex ratio is 1:1 and denoting with S the density of sterile males then the fraction M t /(M t + S) approximately determines the female population density in the next generation through the relation the last equality justified by the sex ratio assumption.This equation is the basic Knipling's model [15] for SIT, which has 0 and a positive value F * as stable and unstable equilibria, respectively.The positive steady state determines the critical release rate S * : if S > S * = (λ − 1)F * then the pest population will decline, otherwise it will rise.The model can be extended to integrate positive density dependent regulation on the larval stage as a variation of the Beverton-Holt model [22], a curve that has been found to be consistent with the dynamics of several insect populations [16].If N denotes the carrying capacity for the insect population in the environment then The curve defined by y = f (x) has an increasing sigmoidal shape for x > 0 that is tangential to the horizontal axis at the origin.Thus, x = 0 is a solution of f (x) = x that is always stable and, for relatively large values of S, the only solution.However, two additional positive equilibria appear if the density of released sterile males is smaller than the critical value S * = (λ − 1)N/4.Therefore, to assure extinction, S has to be larger than S * .Normalizing the density with respect to N , with u t = F t /N and s = S/N , equation ( 2) is rewritten as The impact of SIT on spatial invasions has been usually addressed in continuous time models through the inclusion of diffusion terms.The spatial models for SIT studied in [18], for instance, are based on the non-spatial framework proposed by [6] and provide threshold conditions for the existence of traveling waves that invade or retreat, either when the release of sterile insects is kept at a constant level or when it is considered variable.
3. Dispersal and pest traveling fronts.We integrate pest dispersal through redistribution kernels, which are probability distributions for the distances that insects move.For insect species, movement is generally well described with dispersal kernels that have fat tails [16,20,24] or with Gaussian kernels, being the latter associated also with Fickian diffusion models in continuous time.
Let K be the redistribution kernel corresponding to the pest population, that is, we assume that K(x) ≥ 0, K(x) = K(−x) and ∞ −∞ K(x)dx = 1.Then, the spatial dispersal density under sterile males release effort is given by with g defined as in (3).Next we address the question about a traveling wave on invasion: under what conditions does an established pest population invade a contiguous area initially free of these insects?Traveling waves integrating the release of sterile insects were studied in [18] for continuous time.One of the main achievements there is to find threshold values for the sterile population that determine waves of invasion or extinction.Our goal is to find similar criteria for the discrete time model (4).We notice first, from equation (3), that the release of sterile males induces a (strong) Allee effect in the non spatial case.It is then plausible to determine the behavior of traveling fronts for the integrodifference equation (4) using results from [26].The conditions required on the recruitment function g are: (i) it has at least two equilibrium points, (ii) g belongs to C ∞ [0, 1], g > 0, g ≥ 0, g(0) = 0 and g(1) = 1, where it is admissible to have another fixed point between 0 and 1.These conditions are satisfied in our problem if 1. s is selected smaller than (λ − 1)/4, If the speed is negative, that is the front moves in the opposite direction, the invasion process fails to get established.
2. g (as defined in equation ( 3)) is normalized with respect to its larger equilibrium u * = 1/2 + 1/4 − s/(λ − 1).After rescaling, g is rewritten as We can assume that the spatial distribution, u(x), of the pest population satifies lim x→∞ u t (x) = 0, lim x→−∞ u t (x) = 1 and u t (x) < 0, connecting the largest and the trivial equilibria of the recuitment function g, see Figure 1.According to Theorem 2.1 in [26] (see the Appendix), a successful invasion front appears with positive velocity c > 0 if and only if Recall now that, according to the non-spatial model (3), successful invasions should occur in the region determined by the points (λ, s) that satisfy s < (λ − 1)/4, see Figure 2. In the spatial case, it is the integral (5) instead that determines the region in the parameter space (λ, s) where successful invasions (c > 0) may occur.This region corresponds to the points located under the curve determined by Φ, see Figure 2. The discrepancy between these two regions evidences the existence of parameter values for which traveling waves retreat (c < 0), and invasion ceases, even if the non-spatial model predicts the opposite, in analogy with the results in [18] for continuous time.This happens because female population reproduction at a particular point in space is not fast enough to compensate the dispersal from that point.

4.
Variable release of sterile insects.Up to this point it has been assumed that the release of sterile insects is done uniformly on space, consistently with the fact that sterile insects are often released by aircrafts and thus their spatial distribution is approximately homogeneous.However, if aerial release is neither cost-effective nor efficient, a frequent scenario in small areas, ground release is the preferred method [12].This leads to the problem of finding an appropriate non-uniform distribution for the release of sterile insects.We start considering an established core area, that is, an area that has not yet been contaminated, or perhaps a space where the insect pest has been expunged by other methods.The core area does not need to have crops in it, maybe it is convenient to keep it free of pest if there is danger of it being transformed into pest reservoir.This area is adjacent to a control (buffer) zone that separates the core and the area with high pest density, see Figure 3.Control measures are applied within the buffer zone to restrict pest insect movement into the core area.
Implicit in the optimal spatial release of sterile insects is the problem of finding the control zone width necessary to halt the invasion.This problen has been studied for continuous time in [7] through the use of reaction-diffusion equations.In general, however, the relation between the control zone size and the sterile release pressure is not straightforward.As we will see, the numerical approach introduced offers, as byproduct, results in this direction that might be helpful.

4.1.
In search of the best policy: A numerical approach.For the sake of simplicity we assume only one spatial dimension and choose a segment I that will represent the control zone.We assume that at each point a certain number of sterile males has to be released, but as previously noticed, the releases should not necessarly be made with the same intensity at distinct points.After the core area is determined, it becomes desirable to find the optimal release pressure within a minimum control zone extension where the invasion is halted.
Depending on the context confronted, the implementation of SIT generates costs that include those of production, packing and transport to the site where sterile insects are actually released.We define q as the estimated total production cost per sterile insect and generate a reasonable fine partition (p 1 , ..., p n ) of I.If S i denotes the density of sterile insects released at location p i in one generation then qS i represents the cost, in one generation, of sterile release at the corresponding location.We set as our goal to protect the core area for a fixed time (expressed in number of insect generations) at a minimal total cost, which can be expressed as where we restrict the possible values of S i to be chosen from a convenient discretization that depends on the particular application.To illustrate the numerical procedure, the following assumptions are used in the example below: (i) initially pest insects are uniformly distributed at the carrying capacity in a limited patch, (ii) sterile release is maintained constant at each generation, (iii) insects disperse with a Gaussian kernel.We assume the natural reproductive female growth per generation to be λ = 15, which is among the typical values found for instance in the light brown apple moth (Epiphyas postvittana), [10].Although it has been observed that under normal conditions the apple moth might present dipersal distances about 0.1 Km each generation [23], for purposes of illustration we assume a worse scenario, with dispersal speed of 1.0 Km/generation.We remark again that, if necessary, the scheme is general enough to allow the use of leptokurtic kernels, i.e. kernels with thick tails, commonly used for describing dispersal capabilities associated to several insect species [16,20,24].The goal is to find the optimal among the policies that do not allow insects to reach the control zone boundaries, a condition expressed as u = 0 at the endpoints of the segment.(7) We notice that the requirement might be applied to other locations within the control zone, which evidently would induce different results.We minimize the production cost functional (6) using simulated annealing (see the Appendix for a brief description).The feasible states for the stochastic search are the vectors (S 1 , ..., S n ) that induce condition (7) through equation ( 4) at all times.The resulting approximated optimal sterile release policy, guarantees that pest insects do not cross to the core zone, halting pest dispersal or extinguishing its presence in the control zone.Results are presented in Figure 4 where sterile release is suggested to overflow pest density within the first two kilometers after which it jumps to a lower value and smoothly decreases along the control zone.The computations suggest that the control zone length might be reduced to 8 Km.

5.
Discussion.Integrating dispersal with sterile release and density dependence in the pest population suggests the possibility of overestimating the minimum value of sterile pressure necessary for pest extinction: here we determine the existence of a region in the (λ, s)−parameter space where the traveling wave of invasion retreats, even if the sterile pressure is less than the threshold predicted for successful control.This observation directly leads to the problem of searching for optimal policies that produce distributions of sterile release over a delimited space and avoid the possibility of insufficient sterile insects in the field.In this note we solve the problem through the use of stochastic search optimization, as illustrated in Section 4, where an overflow of sterile insects is suggested only over a limited extension around the pest original location.Incidentally to the computations of the sterile release policy distribution, we are able to determine an appropriate control zone size from an initial guess that could be conveniently resized (reduced or expanded).
In the example of Section 4, we initially consider a control zone that extends 9.5 Km beyond the pest boundary location but the computations suggest it could be reduced.It is worth noting that, in contrast with classical diffusion approaches, integrodifference equation models permit the treatment of cases where long-range dispersal is observed, a characteristic of the movement of several insect species [16,24,25].The computation of approximately optimal release policies and appropriate control zone sizes would be of use for the control phase of benefit/cost analysis carried on scenarios subject the sterile insect technique [19], reducing the risk of making poor decisions.The model presented here in one spatial dimension is amenable to extensions that may include the combination with other control measures.These characteristics provide considerable flexibility to test management options in a complete economic analysis.
We emphasize the fact that we are considering conceptual models, which contain inherent limitations.Nevertheless, they are particularly useful to generate a theoretical framework for supporting management decisions.The practical implementation of these tools requires careful parameter estimation for the model and the dispersal kernels, in addition to the consideration and cautious analysis of many other factors that might affect the insect population survival and spread, like wind, rainfall and seasonal temperatures [7,21].
Simulated annealing.For the sake of completeness we include a brief description of how simulated annealing works.For more detailed accounts the reader is encouraged to look at the reference [13].Essentially, the algorithm runs a search in a space of states for an element that globally minimize (or maximize) the value of a function.In our case, the space of states is given by a multidimensional lattice where each point is a vector with the feasible controls (space units and sterile insect pressure) at each time as entries.The grid is chosen in such a way that it has practical significance and computational feasibility.
Once the cost functional and the finite set of states have been defined, the following algorithm can be used to find the minimum.First choose an initial state and define an inhomogeneous Markov chain X(k) that evolves as follows: if X(k) is in state i, i.e.X(k) = i, choose randomly a neighbor state j that can be reached in one transition step.If J(j) ≤ J(i) then set X(k + 1) = j.Otherwise, if J(j) > J(i) then set X(k + 1) = j with probability q = exp − J(j) − J(i) T (k) , (the parameter T (k) > 0 controls the algorithm speed) and X(k + 1) = i with probability 1 − q.It can be shown [13] that the probability of choosing an element from state space that minimizes J approaches 1 when the temperature parameter tends to 0. The choice of an appropriate sequence T (k), called cooling schedule, is important: if the cooling is too fast, the chain may get stuck in a local minimum.So there is a balance that has to be reached, where you want a reasonable time of computation to observe convergence, but slow enough to avoid convergence to an element that is not a global minimizer.

Figure 1 .
Figure1.A feasible spatial distribution u(x) for the invasive pest amenable to travel to the right or left with speed c.The horizontal axis is the one dimensional space.If the front has positive speed eventually every location x will saturate with invading organisms.If the speed is negative, that is the front moves in the opposite direction, the invasion process fails to get established.

Figure 2 .
Figure2.In the region s < (λ − 1)/4 the existence of two nontrivial equilibria is guaranteed in the non-spatial model.However, the curve Φ(λ, s) = 0 crosses it defining two subregions where the invasive front speed, c, has opposite signs.Thus, in the subregion where the invasive traveling wave retreats, that is c < 0, eradication is anyway sucessful.

Figure 3 .
Figure3.After the core area is determined, it becomes desirable to find the optimal release pressure within a minimum control zone extension where the invasion is halted.

Figure 4 .
Figure 4.The control zone is initially guessed to be the segment [−10 4 , 10 4 ], measured in meters.(a) Evolution of pest dispersal assuming no sterile release with initial uniform distribution of 10 3 m around 0 (not shown).After an initial decay pest insects saturate the zone at the sixth generation.(b) Computed optimal sterile release.(c) Evolution of pest dispersal with the optimal control from (b).Notice the changes of axis orientations (time is reversed) for better visualization.No insect crosses the boundary control zone.The policy is applied constantly in every generation.For this example we consider female reproduction growth equal to 15 per individual per generation, together with Gaussian dispersal with speed of 10 3 m/generation.