Modeling diffusive movement of sterile insects released along aerial flight lines

ABSTRACT The redistribution of insects that are released from an airplane is described using a diffusion equation to derive optimal spacing of flight lines and time interval between flights to achieve a reasonably uniform spatial distribution of released insects and to minimize costs. This optimization is done based on relative costs of sterile males and of flying time. An example is presented using data on the Mediterranean fruit fly, Ceratitis capitata (Weidmann), from the Moscamed program in Guatemala. For the parameter values used, the cost (sterile males plus flights) is minimized when insects are released at intervals of approximately every 2 days (for daily mortality μ = 0.24) to 10 days (for μ = 0.04) and flights are spaced at 300 m (for μ = 0.24) to 600 m (for μ = 0.04) apart, depending on daily medfly mortality estimates, which vary widely in the literature. A simpler approximate method of optimization is then presented based on a relationship observed in the optimization results between flight-line separations and the standard deviation of the distribution.


Introduction
Certain methods of insect pest control, such as the inundative or augmentative release of predators or parasitoids as well as the sterile insect technique (SIT), rely on released insects moving and finding the wild insects that are the target of the control program. Predators and parasitoids must find their prey or hosts and often use kairomones, which may be host pheromones or other host odors. The SIT also relies on the ability of the sterile males to find and mate with wild females in order to be successful in controlling the target population. Part of the efficiency of the SIT involves the movement of the released insects to cover the control area. If individuals move, they may intermix with the wild population, allowing sterile individuals to encounter wild females, which is critical to SIT. In this paper, we focus on movement of sterile males, although the results also apply to other categories of insects.
Various models have been proposed to describe the movements of animals. The model we use is a diffusion equation. Originally derived to describe the diffusion of heat along a metal rod (Fourier 1822), diffusion models have also been used to describe the movement of discrete entities, such as atoms, genes and, more recently, insects (Skellam 1951;Turchin 1998). Starting with seminal papers by Fisher (1937), Dobzhansky and Wright (1947) and Skellam (1951), diffusion equations have become widely used in ecology (Okubo 1980;Kareiva 1983;Shigesada & Kawasaki 1997;Turchin 1998) and, more recently, in pest control (Lewis & van den Driessche 1993;Bouyer et al. 2007; Barclay et al. 2011). They have the advantage of simplicity, as the equations deal with only one process: random diffusion in all directions (but see Holmes 1993). They also appear to be of potential use in describing the movement of insects after being released from airplanes, provided wind is not a major factor and movement of the airplane does not result in much displacement while the insects are falling. However, in their simplest form, they have the disadvantage of having only a single parameter, D, the diffusion coefficient, thus disallowing an evaluation of the effects of other relevant factors.
In field operations involving SIT, insects are released from airplanes over flat terrain or from helicopters over mountainous terrain (Vargas et al. 1995), flying on specified parallel flight lines over areas to be protected. Ignoring convective effects of wind, we want to find the density of insects at the midline between two flight lines and optimize insect delivery at minimal cost after a certain number of days following the release. The major variables that are expected to affect the results are (1) number of insects released per flight, (2) interval between releases (days), (3) distance between flight lines (km), (4) rate of diffusion of insects, and (5) daily mortality of the insects. We expect interactions among all these variables. In addition, altitude of the airplane may be important, as releases at higher elevations above ground level tend to spread out more before landing on the ground (Vargas et al. 1995). Since the rate of diffusion and mortality are species characteristics and are not readily under the control of the pest manager, the optimizations will be done for several values of each.
Specific objectives in this paper are (1) to quantify the extent to which insects released from the air distribute themselves throughout a randomly distributed wild population and (2) to optimize the flight-line separations and flight intervals, accordingly. Our use of a diffusion equation describes passive diffusion in the absence of wind and other disturbing factors. If wind is present, but gentle and not turbulent, then the predictions of models with no wind will simply undergo a translation in space in the direction of the wind and consistent with its speed (Turchin 1998).

Modeling description and results
2.1. Diffusion models of dispersal 2.1.1. The simple model with random deaths of insects Let U(x,y,t) be a density function describing the relative density of insects at a point (x,y) at time t, and having units "number of insects per unit area". Then the two-dimensional diffusion equation with fixed instantaneous death rate m and diffusion coefficient D is @Uðx; y; tÞ=@t D D ð@ 2 Uðx; y; tÞ=@x 2 C @ 2 Uðx; y; tÞ=@y 2 ÞÀmUðx; y; tÞ: (1) Assuming one insect is released at a point (0,0) in space at time zero, Uðx; y; 0Þ D 0; for ðx; yÞ 6 ¼ 0.
There are no boundary conditions as diffusion is assumed to take place over the entire (x,y) plane. Equation (1) then has the following solution: Uðx; y; tÞ D ð1=ð4pDtÞÞ Â ½expf ¡ mt À ððx 2 C y 2 Þ=4DtÞg for t > 0: (2) Equation (2) is a modification of a bivariate normal distribution and has variance 2Dt and is centered on the point of release. The two parameters are (1) D, the diffusion coefficient that describes the rate of diffusion and has units distance 2 /time, and (2) m, the instantaneous rate of mortality per unit time. The relationship between m and the discrete rate of mortality per unit time, q x , for t D 1, as calculated from daily measurements, is q x D 1:0 Àe ¡ m ; or inverted, m D ¡ ln ð1:0 À q x Þ where q x is proportional mortality per unit time as displayed in a life table (Poole 1974), and ln is the natural logarithm. For small values of m (say, m < 0.1) and t D 1, q x % m. In addition, if an accumulated mortality (1.0 ¡ Q) is documented over several time periods, t (e.g. days), where Q is the cumulative survivorship, then m can be estimated from 2.1.2. Release of insects from a flight line of length 2L If the initial distribution of insects is along a line of length 2L km, the diffusion equation governing this case is the same as Equation (1), but has the initial condition: one insect is released along each unit length of the flight line of length 2L km, which is centered at y D 0, and Uðx; y; 0Þ D 0, for x 6 ¼ 0. The solution then is and U(x,y,t) has units of insects/km 2 . The parameter U K is a constant defined as one insect per linear kilometer of flight line (i.e. U K has a value of 1.0). At time t D 0, all the insects are concentrated on the line at the frequency of one insect per unit length (e.g. km). For a derivation of the solution to this model, see Appendix 1. This is also given in more general form by Haberman (2004). As Equation (4) is a modification of the normal distribution, it is not integrable with a finite number of terms, but can be solved numerically in a variety of ways, such as using Simpson's rule (Dorn & McCracken 1972) or Gauss' error function of the normal distribution (see Appendix 2). We use flight lines of length 2L, centered on y D 0, and solve the integral by Simpson's method using a step size of 0.0002 in the dummy variable ξ. Six data files were produced that each contained values of U(x,y,t) for 10 values of daily mortality (m D 0.02À0.2 in steps of 0.02); 321 values of x (0.0À16.0 in steps of 0.05) and 120 values of t (1À120), totaling 385,200 values for each file. Each file had one of six values of D (0.002, 0.005, 0.02, 0.05, 0.2 and 0.5); these values were chosen to span the range of values of D for medflies (0.005) and tsetse (0.46). These files were used as catalogs of function U(x,y,t) values to construct the graphs in Figures 1À7 as well as in the Supplementary Material. The use of 321 values of x facilitates inclusion of eight flight lines on each side of a midline of interest where the flight lines are between 0.1 and 2.0 km apart. Values of m from 0.02 to 0.2 likely include the majority of daily sterile male mortalities encountered in the field (being less than 0.04 for tsetse (Hargrove 1981) and at least 0.26 (Wong et al. 1982); daily mortalities of greater than 0.2 would likely mean that most sterile flies would die before mating), and 120 values of t facilitate including eight previous flights separated by as much as 14 days each, for inclusion in computing the sum of contributions to the total density at the midline from several flight lines and several flights.

Effects of flight-line length and diffusion coefficient on patterns of diffusion on a flight line
The diffusion coefficient, D, has been estimated by Plant and Cunningham (1991) to be D D 0.006 km 2 /day for Mediterranean fruit flies, Ceratitis capitata (Weidmann), and by Bouyer et al. (2007) to be D D 0.29 or 0.46, depending on how it is measured, for the riverine tsetse, Glossina palpalis gambiensis Vanderplank. For the treatment below, we use D D 0.005 and 0.5 km 2 /day as being extremes, likely to include most values encountered for insect species to be controlled by SIT. For flight lines whose half length, L, is much greater than three standard deviations of U(x,y,t) (i.e. 3 (2Dt) 1/ 2 ), the integral in Equation (4) is very close to 1.0. The resulting distribution of density along a line parallel to the flight line is shown in Figure 1; for longer lines, the density is flat for most of the length of the line, unless diffusion is quite high.
Since diffusion occurs in all directions, there will be some losses of released insects at the ends of the flight lines. These are expected to be greater with high diffusion coefficients than low ones, and such is seen to be Figure 1. Height of the density function U(x,y,t) along the line x D 0.5 parallel to the flight line (20 km) at times 1, 5 and 20 days after the insect release (t); m D 0.05 in both panels. Top panel: D D 0.005 and there is a slight loss of insects from the ends of the line to the area outside the control area. On day 1, t D 1, very few insects have yet reached the mid-line at x D 0.5, and the density is very close to zero. At both t D 5 and t D 20, the density is still increasing as a result of the low diffusion rate of D D 0.005. This diffusion rate approximates a measured value of 0.006 for Mediterranean fruit flies (Plant and Cunningham, 1991). Bottom panel: D D 0.5. There is a considerable loss of insects from the ends of the line to the area outside the control area. At t D 1, many insects have reached the midline at x D 0.5, and their density is the highest of any shown on the graph. By t D 5 and t D 20, the density is already decreasing as a result of the high diffusion rate. This diffusion rate approximates a measured value of 0.46 for tsetse flies (Bouyer et al. 2007). This illustrates the feature that high rates of diffusion require long flight lines to be efficient and minimize losses to outside the control area. the case (Figure 1) where there are few losses for D D 0.005, but a substantial "spillover" effect at the ends of the flight lines for D D 0.5. This could be countered by releasing sterile males beyond the ends of the flight lines to maintain the insect density at the required value for control at the ends of the flight lines, so that a small region outside the control area also receives sterile males (Cunningham et al. 1980).

Diffusion of insects from flight lines
2.3.1. Sum of density functions, U(x,y,t), between two flight lines Insects will be diffusing in all directions from two adjacent flight lines, both into the space between the two flight lines and also away from the area between the flight lines. If the values of the density functions between two adjacent flight lines are added, we obtain the sum of the contributions from both lines at all points between the lines. The variable x measures the distance to the right of the left-hand line, and the sums of the two density functions are noted every 50 m between the flight lines that are v km apart. Thus this sum is Uðx; 0; tÞ C Uðv ¡ x; 0; tÞ for values of t t, the time interval between flights. At the midline between two flight lines, x D v=2, so this sum simplifies to 2Uðv=2; 0; tÞ.
2.3.2. Contributions from several flight lines to the midline between two adjacent flight lines In addition to contributions of insects to a given midline from present releases along adjacent lines, there may be contributions from past flights unless diffusion is low, daily mortality is high and flight-line separations are wide. We consider the densities, Uðx; 0; tÞ, just prior to a flight, as this should represent the lowest insect density achieved at that point. The time between flights is t days, and the distance of the midline from either flanking flight line is v/2, so we consider the density of insects Uðv=2; 0; tÞ as being the contribution from each of the two flight lines adjacent to the midline just prior to the next flight, where v is the flight-line separation. Assuming symmetry in x, the contributions from n flight lines on each side of the midline from a single flight t days ago are 2Uðv=2; 0; tÞ C 2Uðv C v=2; 0; tÞ Similarly, the contributions from the two flight lines adjacent to the midline for the last m flights are 2Uðv=2; 0; tÞ C 2Uðv=2; 0; 2tÞ C 2Uðv=2; 0; 3tÞ C Á Á Á C 2Uðv=2; 0; mtÞ: Thus the total density contributions, T U , from n flight lines on each side of the midline, each of which contributes insects from the previous m flights, are Here the standard deviation is 1.0, and the inner two density functions cross at their points of inflection, at one standard deviation. The horizontal scaling is distance of flight lines from the midline in the centre. This represents the coverage, on a previously untreated field of width 30 km, of one single flight, t days ago, where the flightline separation of the recommended v D 2ð2DtÞ 1=2 has been used. Note that the coverage in the troughs, which lie on the midlines, is more than 97% of the coverage at the peaks, which lie along the flight lines, as the sum cycles with a maximum of 0.3400 and a minimum of 0.3304. T U D 2Uðv=2; 0; tÞ C 2Uðv C v=2; 0; tÞ C 2Uð2 v C v=2; 0; tÞ C Á Á Á C 2Uððn ¡ 1Þ v C v=2; 0; tÞ C 2Uðv=2; 0; 2tÞ C 2Uðv C v=2; 0; 2tÞ C 2Uð2 v C v=2; 0; 2tÞ C Á Á Á C 2Uððn ¡ 1Þv C v=2; 0; 2tÞ C 2Uðv=2; 0; 3tÞ C 2Uðv C v=2; 0; 3tÞ C 2Uð2v C v=2; 0; 3tÞ C Á Á Á C 2Uððn ¡ 1Þv C v=2; 0; 3tÞ . . . C 2Uðv=2; 0; mtÞ C 2Uðv C v=2; 0; mtÞ C 2Uð2v C v=2; 0; mtÞ C Á Á Á C 2Uððn ¡ 1Þv C v=2; 0; mtÞ: More compactly, this sum of density functions with units of insects per km 2 is where i assumes values from 0 to n ¡ 1 and k assumes values from 1 to m.
It can be shown that if v is twice the standard deviation of T U (i.e. v D 2ð2DtÞ 1=2 Þ; then the midline is at the point of inflection of the U functions between adjacent flight lines, and the sum of the eight values of U from flight lines on each side of a midline and the eight previous flights is nearly flat (Figure 2). This represents a configuration of maximal conservation of sterile males, as they are almost equally abundant everywhere between the two flight lines. Ideally, the heights of the sum of the two densities at all points between the two lines should be close to flat for maximum efficiency of the SIT, and this can be shown to be approximately true if D 0.05 and flight lines are 1.0 km or less apart, with multiple flight lines as described below (Figure 3). They will be very close to being flat, if v D 2ð2DtÞ 1=2 .
Then T U can be used as the density to compute the actual number, N K , of insects to be released per kilometer of flight line and the number, N SK , of insects to be released per km 2 . If N M is the minimum number of sterile males per km 2 required to ensure control, i.e. to result in population growth being stopped, and T U is the total obtained from summing the density functions, U(x,y,t), from 8 flight lines on each side of a midline (i.e. 16 flight lines) from each of the most recent 8 flights (henceforth called "triple 8 fl"), then we can compute the value N K as and Equations (6) and (7) are used as intermediates in the calculation of costs (in Equations (8) and (9)). For most parameter values, the density U(x,y,t) is very close to zero for flight lines greater than four distant from a midline and flown longer ago than the most recent four flights. For large values of D and narrow separations between flight lines, a few of the functions, U(x,y,t), are still positive but small, for eight flight lines distant and eight flights ago. As flight lines will probably have wide spacing for species with large D, for the parameter combinations used here, eight flight lines appear to be sufficient to include almost all contributions to the sum of the U(x,y,t) functions at any midline. This is made more specific in the medfly example given below.
The parameter N M is mentioned above as being the density of sterile males required to achieve control and it must be determined for each control program prior to using the methods given in this paper. The equation is given by Barclay (2001Barclay ( , 2005 to determine the sterile release rate for a population of size M for the adult male component, a rate of increase per generation of λ, a sterile male release rate per generation of S and a competitive coefficient, c, for the ability of sterile males to find and fertilize wild females compared with that of wild males. Then the sterile release rate to hold the population at equilibrium, S Ã , is given for unequal sterile and wild male mating competitive ability. This is exactly the source of N M , because N M D S Ã in this context, perhaps inflated by some factor to ensure success. The value of S Ã (and thus N M ) for a given pest control program depends on the values of M, λ and c, which will in general not be known without some prior investigation, and will vary from one control program to another. In practice, other features such as residual fertility of irradiated males, sterile male survivorship and immigration will also need to be considered in computing N M . In some control programs, this algebra is side-stepped and a simple multiple of the existing wild male population is used instead. For example, in controlling medfly using SIT in the Moscamed program in Guatemala, a multiple of 200 is used as a multiple of the existing wild male fly density (M) to estimate the intended density of sterile males, so that N M is assigned the value of 200 M.

Density function totals at points between the two central flight lines from triple 8 fl
The sum of total contributions, T X , to a point x between two flight lines and from triple 8 fl is ktÞ for i going from 0 to 7 (the first summation) and k from 1 to 8 (the second summation). Values of the sum of densities are shown in Figure  The only case that noticeably violated the assumption above was the case of D D 0.005, m D 0.04, v D 1.0 and t D 7, in which the sum of densities was lowest for t D 1, then for t D 3, t D 5 and highest for t D 7. For D D 0.05 and 0.5, the assumption held.

Computation of costs
We assume here that the major costs of an SIT program are the costs associated with purchase of sterile males and flying time. Other costs are assumed to be factored into the costs of sterile males and flying time, so that we need only consider these two. Table 1 defines several components of the two major costs, and these can be used to obtain the two quantities: C FD is the cost of flying per km 2 per day C SD is the cost of sterile males per km 2 per day is The total cost, C T , per km 2 per day is the sum of C FD and C SD : The information required to evaluate C T requires knowledge of the following quantities: A, B, L, S, C H , C S , N M, m, D, T U , T FL , N FL, v and t (see Table 1 for definitions).

Optimization procedure
The optimization procedure consists of calculating the sums (T U ) of densities (U(x,y,t)) for all combinations of flight-line separation (v) and days between flights (t) and the calculation of the number of sterile males that need to be released per km 2 (N SK ) for each combination of diffusion coefficient (D) and mortality (m), using Equations (4) and (5) for T U . Once this is done, the cost of each combination of v, t and N SK is calculated using Equation (10) for each combination of D and m. Then for each combination of D and m, the values of v, t and N SK that minimize the cost are identified by inspection. Figure 4 shows the optimal flight-line separations, v, and flight intervals, t, for 3 values of D (0.005, 0.05 and 0.5) and 10 values of m (0.02À0.2 in steps of 0.02).
In the absence of specific information to the contrary, the same relative costs of sterile males (US$250 per million flies) and flying time (US$1000 per hour) were used for all three values of D as are assumed for the medflies below. The optimal flight-line separations increase with D and decrease with m, whereas the flight intervals decrease with m but are almost unaffected by D.
2.4.3. Effects of relative costs of sterile males and air time The optimization described above assumes known relative costs for sterile males and for flying time, and the results depend on the relative sizes of these costs. Figure 5 shows the optimal flight-line separation and costs for a range of mortality values (m D 0.04À0.24 per day) for a population in which D D 0.005. In this graph, the relative costs are varied by adjusting the required density of sterile males per km 2 to be 200,000, 1,000,000 (default) and 5,000,000. This shows the differences in optimal outcomes if costs are varied and LÀhalf the length of a flight line BÀwidth of the control area (km) A À total area being controlled (km 2 ) D 2L¢B v À separation between two adjacent flight lines (km) t À interval (days) between successive flights N FL À number of flight lines; N FL D B=v S À speed of the airplane (km/hour) T FL À time to fly one flight line (hours) D 2L/S T FA À time to fly the total area (hours) D T FL N FL D 2LB=vS T FD À time to fly per km 2 per day D T FA =At D ð2LB=vSÞ=2LBt D 1=Svt x, y À horizontal and vertical spatial coordinates t À time since release of insects (days) D À diffusion coefficient (km 2 /day) m À instantaneous death rate of insects (1/day) q x À daily death rate as measured on each day N M À the minimum density of sterile males (number per km 2 ) required for eradication of the pest population, and it is determined independently from models of SIT and other considerations. For the purpose of optimization in this paper, N M is assigned the value of 10 6 /km 2 . U(x,y,t) À function giving relative frequency of insects at spatial coordinates x and y and at time t following release U K À a constant defined as one insect released per kilometer of flight line T U À sum of densities, U, for eight flight lines on each side of a midline for the previous eight flights N K À number of sterile males to be released per kilometer of flight line D U K N M =T U N SK À number of sterile males to be released per km 2 per day D N K =tv D U K N M =ðT U vtÞ λ À required number of flight lines on each side of a midline to obtain adequate accuracy for T U ' À required number of prior flights to obtain adequate accuracy for T U C H À cost of air time per hour C FD À cost of flying per km 2 per day D T FD C H D C H =ðSvtÞ C S À cost per million sterile males C SD À cost of sterile males per km 2 per day D illustrates the importance of knowing these relative costs. The values of 200,000, 1,000,000 and 5,000,000 are not derived from any particular control program and are only used to determine the sensitivity of the system to changes in the relative costs of sterile males and flying time. Any other three (or more) values could have been used for this. The assumption of one million sterile males being required per km 2 (above) thus will affect the results of the optimization. The graphs show the predictable pattern that when the cost of sterile males is relatively low and the cost of flying is unchanged, then the optimization procedure yields wider flight-line separations, whereas when the cost of sterile males is increased, then narrow flight-line separations are optimal. This is a conservative estimate of the change of optimal flight-line separation and flight intervals, as throughout most of the comparisons the cost of flying was only in the order of 10% of the cost of sterile males, and with this imbalance of costs this comparison would be expected to be less sensitive than if the costs of insects and flying time were approximately equal.

Staggered flight lines
In some SIT programs, the flight lines alternate from one flight to the next, so that on one flight, lines 1, 3, 5, … are flown and on the following flight (t/2 days later), lines 2, 4, 6, … are flown. This complicates the computations only slightly, as the contributions to midline densities alternate between these two sets of lines flight by flight. Thus on one flight the contributions from the eight flight lines 2v apart that were flown t days ago are Uðv=2; 0; tÞ C Uðv C v=2; 0; tÞ C Uð2v C v=2; 0; tÞ C Á Á Á C Uð7v C v=2; 0; tÞ: Similarly, from the eight flight lines 2v apart flown t/2 days ago, where v is taken as the distance between adjacent flight lines flown t and t/2 days ago. The sum from the previous 8 flights of the contributions to a midline of interest from 8 flights flown t days ago and eight flights flown t / 2 days ago is for i D 0À7 and k D 0À7. In each summation in this formula we include flight lines from both sides of the midline, whereas in Equation (5), we sum on one side only and then double the value because of symmetry. A comparison of the resulting optimization curves is shown in Figure 6 using ratios of T U for staggered flight lines and T U for regular flight lines. It is seen that staggered flight lines are more efficient than the regular pattern, especially for low values of v and high values of D, m and t.

Example: Mediterranean fruit fly
In order to optimize the flight-line separations (v) and flight intervals (t), the daily sterile male mortality and the diffusion coefficient for the population must be determined. These may be approximated by published values for the species, but may also depend strongly on local habitat conditions. Also needed are the costs of purchase of sterile males and of flying time. The relationships in the Supplementary Material can then be used to estimate the optimal values of v and t, and then Equation (7) will give the number of sterile males that need to be released per km 2 .
2.5.1. Sterile male daily mortality Data on daily medfly mortality can easily be obtained in the laboratory, as was done by Carey (1982) and by Vargas et al. (1984), but these data are generally not applicable to field conditions, as they are collected under ideal conditions with no external sources of mortality from inclement weather, predators, parasitoids, etc. Studies involving mark-recapture are an important source of data, although difficult with species that have a high daily mortality, as a result of low recaptures of marked individuals (Barry et al. 2004). Thus, most field estimates are both of low accuracy and often specific to local habitat conditions. We have tapped several sources for our estimates, and they vary widely, as expected. In order to estimate daily mortality from capture data, we had to assume that proportions of the total captures that were captured each day were accurate reflections of corresponding proportions of the actual population and that insect response to traps does not vary with time or age. We use the relationship given in Equation (3) to obtain m from the graphs given in the publications cited below. Data from the Moscamed program in Guatemala indicate that there is a 50% mortality of sterile males by day 3, which translates into a value of m % 0.23. Cunningham et al. (1980) stated that the weekly survivorship of medflies was 0.65, which put the daily mortality, m, at 0.06. Wong et al. (1982) listed weekly catches at 7.2 for week 1 and 1.2 for week 2, which yields a daily survivorship of 0.26, similar to that of the Moscamed flies. Barry et al. (2002) published a survivorship graph for catches of two strains of medflies at all trap distances; using the catches at days zero and 14, the values of m for the two strains were 0.17 and 0.20. A graph published by Andress et al. (2013) indicates that daily mortality was about 0.06. These values come from different strains and different release habitats and make generalization difficult. The high values of apparent mortality of the sterile males may be a result of starvation, desiccation, being blown away by wind, drowned by rain, eaten by predators, not having   Plant and Cunningham (1991) constant responses to traps over time or other sources. It is important to estimate this parameter accurately, in view of the variability noted above.

Diffusion coefficient
Estimates of the diffusion coefficient, D, for medflies are even less available in the ecological literature. Such estimates as exist also may be strongly biased toward particular habitats and seasons. The diffusion coefficient has been estimated to be D D 0.0058 km 2 /day for Mediterranean fruit flies, which are relatively poor fliers and typically disperse only short distances (Plant and Cunningham 1991). Crude estimates of diffusion can be obtained from published graphs of trapping data. Assuming that insect response to traps does not change with time or distance and that the distribution is approximately normal, then the point at which the trapping data approach zero (on the vertical axis) can be used to estimate three standard deviations from the mean, and then the formula for variance, s 2 D 2Dt; can be used to estimate D if the graph is for only one time value after release. Such graphs are given by Cunningham and Couey (1986, Figure 1), Plant and Cunningham (1991, Figure 3) and Andress et al. (2013, Figure 4, for days 1 and 3). We have visually estimated the point representing three standard deviations for each of these and computed D; these are given in Table 2. The value of D obtained from estimating the third standard deviation from the graph by Plant and Cunningham (0.0054) was very close to the estimate they gave (0.0058) using a more comprehensive estimation method. Estimates of diffusion may turn out to be just as variable and difficult to assess accurately as estimates of mortality.

Example from the Moscamed program in Guatemala
We use an example from a sterile medfly release block located in a coffee plantation in Guatemala to illustrate the computations of cost for a given set of parameter values. The flight lines were 20 km in length and 0.5 km apart, and releases occurred once a week, with staggered flights occurring over lines 1, 3, 5, … one week and lines 2, 4, 6, … on alternate weeks. Thus on a given flight, the distance between flight lines flown that day is 1.0 km and those flight lines are flown once every two weeks; however, since midlines receive contributions from flights on both weeks, the flight-line separation is 0.5 km. The densities at the midline were found as in Equation (4)  Thus the total cost/km 2 /day is C T D C FD C C SD D $0.71 C $109.82 D $110.53/km 2 /day, or 2000(110.53) D $221,060 for the control area per day.

Effects of travel time to control site
So far, we have been assuming that the travel time to the control site is minimal and that the control site is very close to the airport. This is not likely to be true for most control sites. Travel to and from the control site adds to the total air time and thus can be included in our calculations by simply increasing the cost per hour of air time in calculating the time to fly the flight lines and cost of doing so. For example, if the control site is half an hour from the airport, then one hour is added to the flying time. If the aircraft has sufficient fuel to fly safely for five hours and travel time is one hour, then the cost of flying is increased by 25%. If this is inserted into the C H of Equation (8), then the optimization results may be different from those without considering travel time (as in Figure 5).

Optimizing for parallel linear aerial flight lines
To optimize flight procedure for given values of D and m, we use the method immediately above. The cost is calculated for each combination of v and t, with N M determined from field data and considerations of sterile male effectiveness, and N SK calculated from N M , v, and T U using Equations (6) and (7). We then choose the combination with the lowest total cost by visual inspection of computed values. For this, we calculated values of U(x,y,t) for x D 0.02À4.0 in steps of 0.02, for t D 1À180 in steps of 1.0 and for m D 0.04À0.24 in steps of 0.02. We used only D D 0.005 for this as the parameters were chosen to illustrate results applicable to medflies. The results are given in Table 3 and show a considerable range of cost, flight interval (t) and sterile males required (N SK ), depending on daily mortality. At the low end of m (0.04), longer flight intervals are suggested, and fewer sterile males need to be released, and the high end of m (0.24), shorter flight intervals are suggested and higher numbers of sterile males need to be released. It is notable that the optimal number of sterile males to be released per km 2 (and therefore cost of sterile males) is approximately linearly related to the daily mortality, m (Table 3). In addition, for the medfly case, flight-line separation, v, is almost equal to twice the standard deviation of the distribution, 2(2Dt) 1/2 for all values of m, thus ensuring a flat curve for T U , and optimizing efficiency of the use of sterile males.

Approximate method for optimizing
The close relationship between v and the standard deviation of the distribution, U(x,y,t), that is evident in Table 3, can be used to derive a simpler approximate method for optimizing the values of v and t. This method contains two assumptions: (1) the relationship between v and the standard deviation of the distribution observed in Table 3 is assumed to hold for other values of the parameters, and this is supported by Table A1; (2) the midline being used for computations is assumed to be well within the control area so that edge effects are not evident and diffusion is assumed to be in equilibrium, with equal numbers of insects entering and leaving the area of computation. This approximate method is easier to use than the full estimation method outlined above.
Given the values of D and m, take v D 2 ð2DtÞ 1=2 for each possible value of t, t D 1, 2, 3, …, k. The original motivation for this v value was to ensure that the points of inflection of the normal distribution of insects from two adjacent flight lines coincide (Figure 2). A calculation then shows that from a single flight t days later, the ratio between the highest density and the lowest density in the field is 50.7À49.3. The density from all previous flights is essentially uniform. For n flights (Table 3 was computed using n D 8), so N SK ðexpð ¡ mtÞ C expð ¡ 2 mtÞ C Á Á Á C expð ¡ nmtÞÞ D N M D 10 6 insects in the medfly example are required.
Using the formula for the sum of a geometric series with n terms, the equation gives the number of steriles to be released per km 2 as Now, to find the optimal value of t, we must minimize the formula given in (10): Total cost, C T , per km 2 per day is C T D C H = ðSv tÞ C U K ðC S = 10 6 Þ ½N M = ðt v T U Þ: Replacing the terms U K N M /(vT U ) U K N M =ðvTUÞ as in Equation (7) by N SK as in Equation (12), gives Table 3 has been computed from a cost of $250 for 10 6 insects, and a total flying cost that is an average of $5 per kilometer actually flown over the field, once the transit to and from the field is included. Using these Table 3. Optimal values of flight-line separation (v, km), flight interval (t, days), sums of densities at the midline, T U , and numbers of sterile flies to be released as well as associated costs (US$/km 2 /day), for values of m from 0.04 to 0.24 with D D 0.005 throughout. Here travel time to and from the control site is one hour and total flying time is limited to five hours per trip, thus increasing computed flying time by 25%. (a) These results were calculated using the methods outlined prior to Equation(10). (b) These results were calculated using Equation (14) and the plots of C T (t) against t. "Inclusion" indicates the number of flight lines on each side of a midline and (in parentheses) number of prior flights to include 99% and 95% of total T U . Given values of D and t, 2ð2DtÞ 1=2 is two standard deviations of the distribution of insects from a single flight line. Note that the values for steriles/km 2 /day are very accurately predicted by the approximate method, which were calculated in the section on medflies. C T ðtÞ D 5=ð2ð2DÞ 0:5 t 1:5 Þ C 250ðexpðmtÞ ¡ 1Þ=ðt ð1 ¡ expð ¡ nmtÞÞÞ : For medflies, D D 0.005, so with n D 8 flights, When C T (t) from Equation (14) is plotted against t, this is a convex (i.e. concave up) curve, with a minimum value that represents the optimal values of t (Figure 7). A hand-held graphing calculator, or graphplotting software on a computer, can easily be used to obtain the minimum value, and the t (which must be a positive integer) at which it occurs. The results using this method are given in Table 3 (panel b) and give excellent agreement with Table 3 (panel a) for the number of sterile males to be released, the total daily cost, and for t, as well as reasonable agreement for v.
As an example of the approximate method of calculation for how the m D 0.24 column of Table 3 (panel  b) is produced, using a basic calculator, from the formula in Section 2.5.6, we use Equation (15) and obtain the following: The first instance when the cost increases, signals to stop calculating, as the cost will only increase beyond this point. In our example, the best value of t is 2 days. This then gives v D 2 sqrt(2 Ã 0.005 Ã 2) D 0.283. The number of steriles to be released per day is then 1; 000; 000 £ ðexpð0:48Þ ¡ 1Þ= ð2 Ã ð1 ¡ expð ¡ 3:84ÞÞ D 315; 000: 2.5.7. Determination of required number of flight lines and prior flights to estimate T U The development above uses eight flight lines on each side of a focal midline and eight prior flights to ensure adequate estimation of T U . This was determined from visual inspection to be satisfactory for most parameter combinations. Certain parameter combinations, while perhaps seldom encountered in reality, require many more flight lines or prior flights, such as large D coupled with small m, v and t. Determinations of T U were made for numbers of flight lines on each side of a midline from 1 to a maximum of 200 and also for numbers of prior flights from 1 to a maximum of 200 to obtain asymptotic estimates of T U . The accuracy was ensured by requiring that at least five consecutive estimates of T U be within 0.0001 of each other before accepting the result. The required number of flight lines (λ) on each side of a midline and prior flights (') to obtain 95% and 99% of the asymptotic value of T U were then noted. The optimization in Table 3 indicates that for medfly parameters, optimal values of tau vary inversely with m, while optimal values of v vary inversely with the square root of m and this has the effect of rendering the values of both λ and ' to be relatively constant with changes in m, v and t. The last two lines in Table 3 show the values of λ (not in parentheses) and of ' (in parentheses) to obtain 95% and 99% of asymptotic T U using staggered flight lines. The same patterns were Figure 7. Curves for optimizing t using the approximation given in Equations (13) and (14). In both curves, the value of D is 0.005, while in the top curve m is 0.04 and m D 0.24 in the bottom curve. In both cases, eight flights have been used. The location of the optimum value of t is well defined (t D 2) for m D 0.24, but not so well defined (t D 10) for m D 0.04. For m D 0.24, the cost would be notably greater if t were assigned to 1 or 3, whereas the difference in cost would be minimal if t was misassigned to be 8 or 9, or 11 or 12 when the optimum was 10 (as with m is 0.04).
found to hold for regular flight lines. Note that for the parameter values used, the values of λ are very low, because v is relatively large for such a small diffusion coefficient, whereas the values of ' are larger (between 6 and 11) because t is large only when m is small, and this results in survival of insects over several flights.

Assessment of the models used
The modeling approach used here is simple in concept, and the results are easy to understand and incorporate into a pest control program. Diffusion in two dimensions is symmetric in the model used here, and this is both an asset and a liability. Our model can be seen as a basic model that adequately describes insect movement in the absence of disturbing factors. However, wind is ubiquitous and will change the patterns described here in ways that are not easily predicted (Baker et al. 1986). If the wind is mild and relatively non-turbulent, then the effect is mostly just to displace the distribution of insect numbers in the direction of the wind and by an amount consistent with its speed (Okubo 1980;Turchin 1998); this can be described by an advectionÀdiffusion equation. However, if the wind is turbulent, then the situation is less predictable. Insects released from the air are likely to be swept away by wind and may leave the control area entirely. In addition, there are only two parameters embedded within the model, the diffusion coefficient and mortality, and this limits the flexibility of the approach. If it is possible to plan sterile releases with some flexibility, then the adverse effects of wind, rain and temperature can likely be minimized by planning releases to avoid the effects of inclement weather.
There have been many modifications to the original diffusion equation to make the model more useful for ecological purposes (Okubo 1980;Turchin 1998). Some of these could be used to expand the usefulness of the approach given here. For example, wind could be addressed by using advection-diffusion equations (for laminar air movement) and by using the models for turbulence (Okubo & Levin 2001), while aggregation of wild and released insects could be addressed with the model for attraction and population pressure and of correlated movement (Okubo & Levin 2001). Plant and Cunningham (1991) used a model of settled and non-settled portions of their medfly population and got a better fit to the data than by assuming homogeneity of the population. If this turns out to be true for most species, then it would be a useful approach to improve the general treatment presented here.

Data requirements
The data requirements for the optimization were stated near Equation (10). The information required to evaluate the total cost, C T , is A, B, L, S, C H , C S , N M , m, D and T U . From these we can obtain T FL , N FL and N SK, and then from the optimization we can obtain t and v. This gives the five quantities required to set up a SIT program: C T , N SK , T U , t and v (i.e. total cost per km 2 per day, number of sterile males that need to be released per km 2 , the relative density, flight interval (days) and flight-line separation (km)). These appear to be rather stringent data requirements, but A, B, L, S, C H , C S and N M , are all known or can be determined ahead of the establishment of the SIT program, T U is available from this paper, including the Supplementary Material (or can be derived using Equation (4) as was done here), T FL , N FL and N SK are easily derivable from the foregoing, and only m and D must be determined empirically. Alternatively N SK can be computed by the simpler approximate method given above. Although m and D are biological characteristics of the target species, they are also affected by the local habitat and perhaps also by inbreeding or other characteristics of the rearing facility. Mortality, m may be modified to some extent (McInnis et al. 2002) and perhaps D as well, but their limits are determined genetically, and these parameters should be measured prior to the onset of a SIT program. If the breeding facility puts limits on the quality of the sterile insects, then these too may be improved (Calkins & Parker 2005).

Requirement that insect density must be greater than the minimum everywhere
The optimization method outlined above depends, among other things, on N M , the minimum number of sterile males per km 2 needed to be present to ensure control of the pest. This number should be equaled or exceeded everywhere and at all times to ensure that all parts of the wild population have an over-flooding ratio sufficient to cause a decline throughout the control area. This assumes a uniform spatial distribution of the wild population, and it may result in considerable wastage of sterile males that are in a position, as a result of diffusion and mortality, to be in excess of this minimum, although for higher values of D the sums of the T U function are nearly flat. It is possible that some lower value of T U would suffice, so that at some locations the local population would be decreasing as a result of excess sterile males and at some locations the population would be increasing, resulting in a net decrease. However, this might result in a clumped spatial distribution. It is, therefore, advisable to try to exceed the minimum value throughout the control area.
The second factor to consider is that the optimization predictions depend on the relative costs of the sterile males and the flights (shown in Figure 5). If an extremely large over-flooding ratio is desired, then that can be taken into account in allocating N M in Equation (6) to allow an appropriate optimization to be achieved.

Assumption of minimal density being at the end of the flight interval, t
We chose to calculate our required sterile male release rates based on the sum of the heights of the density functions at the midline between two flight lines at a time just prior to the next flight; this was on the assumption that the sum would be minimized at that time. This turns out usually to be the case, but not always. However, even when that sum is not the lowest at the midline, the difference is likely not great and the approach still appears to be appropriate. The main violation of this assumption was for D D 0.005 and m D 0.04, a combination that may turn out to be uncommon.

Patterns of dispersal
The question of whether sterile males disperse in a similar manner to the wild individuals needs addressing for each SIT program. Sterilization has been shown repeatedly to affect the survivorship, mobility, mating competitiveness and behavior of many species, although improvements in sterilization techniques have been made (Calkins & Parker 2005). Vreysen et al. (2011) have demonstrated that in their study tsetse showed similar patterns of dispersal for steriles and wild males. By releasing marked sterile Mexican fruit fly (Anastrepha ludens, Loew) from a single point, data in Enkerlin (1987) demonstrated that under the conditions of a non-suitable habitat, sterile flies disperse in an eccentric fashion following the direction of dominant winds (northwest in that case), whereas, in a suitable homogeneous habitat, the dispersal pattern was random. In the same study, using a simple method proposed by Hamada (1980), the mean dispersal distance of sterile flies released in the non-suitable habitat was found to be 1.3-fold greater than flies released in a suitable habitat. These findings are useful in assessing fruit fly dispersal ability under different environmental conditions and thus for planning sterile fly releases. If sterile insects move in patterns similar to wild insects, then they too will move toward (or away from) clumps and control becomes easier; however, if they do not disperse in patterns similar to the wild insects, then control will be more difficult, as they will not encounter wild females as often as the wild males will. When chilled sterile flies are released by airplane, during the first moments after being released the population may maintain a relatively uniform distribution. Once these flies reach the ground and warm up and start moving according to ecological and environmental factors present in the release site, they may tend to clump. However, on homogeneous areas (e.g. forest or large areas covered with a single host orchard), the sterile population may maintain a more uniform or random distribution. These patterns of movement of sterile and wild insects should be determined in order to assess their consequences. This has been done with some species (e.g. tsetse by Vreysen et al. 2011), but movement may vary from one location to another, so it would be useful to do these estimations for each SIT program (Enkerlin 1987). If it is determined that sterile males do not disperse in a manner that facilitates intermixing with the wild population, then narrower flightline separations would be advisable.

Aerial release of insects: altitude
While many estimates of diffusion are based on ground releases from a point, when insects are released from an airborne vessel, they must travel downward some distance before reaching the ground. On the way down, it seems probable that they spread out to some extent, so that when they land they cover some area on both sides of the flight line. In this way they would get a head start on insects that were released at a point or along a line on the ground. Thus it is reasonable to displace the time by a certain time interval to allow the diffusion model's predicted insect distribution to match that where the insects land. This could account for the large value of the estimate of the diffusion coefficient obtained from the data of Andress et al. (2013), as they released their sterile medflies from the air, whereas the others in Table 2 released the insects from the ground. It may be significant that the value of three standard deviations for day 1 was almost identical to that for day 3 (Andress et al. 2013), as most of the spread in each case may have resulted from spreading out during the fall from the aircraft. One way of allowing for this would be to fit the data to the model, varying t until a good fit is obtained at some starting value of t, t s , then using t s as the initial value of the model to make predictions. The value of the starting time would presumably be a function of altitude, condition of the insects when released, wind, rain, etc., and these would need to be involved in any estimation of t s . It appears reasonable that the value of t s would increase with altitude of release, but so would the displacing effects of wind, so that a tradeoff between having the insects spread out and losing them to wind would have to be considered on each flight.

The Medfly case
Using values of m and D obtained from the literature, the recommendations for control of medfly by SIT are to use flight-line separations closer than are presently used, ranging from 300 to 600 m (Table 3), and to use flight intervals of 2 (for high mortality) to 10 days (for low mortality). These are not too different from those being used in a number of action programs that integrate SIT. The FAO and IAEA guideline (Enkerlin 2007) states that, for chilled fruit fly adult release, a range of 268 m (880 feet) to 536 m (1760 feet) apart is normally used and flights occur every seven days, although every three days is recommended due to high mortality rate. Unfortunately, many programs release at 500 m apart and every seven days because of high operational costs. The balance between flight-line separation and flight interval has been shown here to depend on the relative costs of sterile males and flight time, so that where the cost of release of insects is much less than the cost of production of the insects, it makes economic sense to have more flights and fewer insects being released, both per flight and also in total. Although more flights means that release becomes more expensive, fewer insects need be released, because the local density of insects is maintained at a higher level as a result of densities at the midlines between flight lines being higher (Figure 3). Insect mortality also plays a major role in the release frequency, with more frequent flights and narrower separation between flight lines being recommended (Table 3) with high mortality. In fact it appears to be wasteful of reared insects to have wide separation between flight lines and base the numbers released on the minimum sum of densities, when most of the area between the flight lines will be much higher than the minimum required.
Data from the Moscamed program (Ponciano et al. 2015, in preparation) confirm that there is more uniform distribution of insects between the flight lines when they are 500 m apart than when they are 1000 m apart. On the other hand, the measures of "flies per trap per day", the sterile to wild ratio, and percent recovery in monitoring traps were not statistically different between those two flight-line separations. However, the better coverage was enough to convince Moscamed to maintain flight lines at 500 m apart. Also, the added expense of twice as many flight lines (56 minutes flight time at 500 m separation compared with 28 minutes at 1000 m separation, at US$935 per hour) should be offset somewhat by the lower number of flies that need to be released as a result of a greater sum of fly densities during fly diffusion.

The simpler approximate method
This method is based on (1) the observation above that when v D 2ð2DtÞ 1=2 the T U curve between two flight lines is nearly flat and (2) the derivation surrounding Equation (12) for an alternative method in calculating N SK . The predictions for the optimal values of t and v are in remarkable agreement for the two methods of calculating N SK , i.e. Equations (7) and (12), for 3 values of D and 10 values of m. Clearly it appears feasible to use the simpler method for calculating optimal values of sterile male releases for the range of parameter values considered here.

Conclusions
1. Predicted optimal flight-line separations depend on flight intervals and the diffusion coefficient whereas flight intervals and overall costs depend sensitively on sterile male mortality (m) and relative costs of sterile males and air time (Table 3 and Figure 4). 2. For medflies, with low values of D (0.005) and high values of m (0.24), narrow spacing of flight lines and frequent flights are recommended (Table 3). 3. For the same number of flight lines per unit time, flying staggered flight lines on alternate flights is slightly more efficient than flying all the flight lines at the same time (Equation (11) and Figure 6). 4. The most important biological information to obtain prior to the aerial release of sterile males is sterile male mortality and diffusive ability. Both appear to be highly variable for a given species, likely depending on habitat and weather (Table 2). 5. The simpler approximate method of optimization gives an excellent estimate of the number of sterile males to be released (Table 3 (panel a) and (panel b)) as well as the optimal values of t and v.
From Equation (4) it follows that Uðx; y; tÞ D U K ð4pDtÞ ¡ 1=2 exp ¡ mt ¡ ðx 2 =4DtÞ I;where I is the integral in Equation (4), and using Gauss' error function, I becomes I D 0:5½erf ðz high =ð2Þ 1 2 Þ ¡ erf ðz low =ð2Þ 1=2 with z high D ðL ¡ yÞ=ð2DtÞ 1=2 and z low D ¡ ðL C yÞ= ð2DtÞ 1=2 Appendix 3. Table for accuracy of the approximation using the standard deviation method Table A1. Optimal values of t and v for two values of D and seven values of m with v being two standard deviations ðv D 2ð2DtÞ 1=2 Þ using the optimized t to compute v. Also, the number (nfl) of flight lines to obtain 99% of the total T U was 2 for all values of m except for m D 0.02, for D D 0.05, and 5 flight lines for D D 0.5. The number of flights (n) required for 99% of T U (in parentheses) was much greater than the required number of flight lines in all cases.