# Simulation Models Evaluation of Pest Resistance Development to Refuge in the Bag Concepts Related to Pioneer Submission

Letter Report

By

Michael Caprio 1 , John A. Glaser 2

1 Misssissippi State University, Depa rtment of Entomology and Plant Pathology, Mississippi State, MS 39762-9775, 2 US Environmental Protection Agency, Nationa l Risk Management Research Laboratory, Sustainable Technology Di vision, Cincinnati, OH 45268

January 22, 2010 Version 1.0

Sustainable Technology Division National Risk Management Research Laboratory US Environmental Protection Agency Cincinnati, OH 45268

Keywords: PIP corn variety, PIP crop registration, resistance management, simulation modeling, refuge in bag

A portion of the work outlined in this report was accomplished under contract with Dr Michael Capri o2 Introduction

The USEPA, under its administration of the Fede ral Insecticide Fungicide and Rodenticide Act (FIFRA), requires the registration of all pesticides and pesticidal ma terials. The GM crops containing pesticidal traits are subject to FIFRA registra tion requirements. Concer ns relating to future environmental effects of these crops can be i nvestigated only using simulation models. Model application in this instance offers an analysis of the useful lifetime for the pesticidal traits. USEPA has clearly supported the use of these traits as replacements of broa d spectrum pesticides that have significant environmental footprints.

Background

Resistance management (RM) simulation models ar e designed as deterministic, stochastic, and spatially explicit analytic al methodologies. Simulation models can provide a realistic assessment of the risk of resistance evolution given the allele frequenc ies used as initial conditions. New concepts of PIP crop deployment such as refuge-in-the-bag require a de tailed evaluation to ensure the claims associated with each concept offer the sustainable protection of the crop. Simulation models offer a means to evaluate the threat of crop loss through the development of resistance in the near future.

Purpose

This collaborative research effort was designed to evaluate the relative merit of refuge-in-bag concepts as submitted by Pioneer Hi-Bred Internationa l, Inc. This research specifically analyzed the

relative risks related to the control of western corn rootworm in block refuge deployment and 5% seedmixture scenarios.

Research Design

To understand the relative risks rela ted to the evolution of resist ance in western corn rootworm between the currently mandated block refuge depl oyment and a 5% seed-mixture deployment as submitted Pioneer Hi-Bred International, Inc., two di fferent models were developed and evaluated. The two models consisted of a modification of a spatially -explicit stochastic mode l (Caprio et al. 2006) and a simpler, frequency-based deterministic model. Th e latter could be run in a single simulation mode with a graphical interface to en ter parameters or in a risk-a ssessment mode capable of running

thousands of simulations to estimate the effects of parameter uncertainty. Realistic estimates for3 several parameters using in the deterministic mode l were developed through use of the stochastic

model.

Impacts and Outcomes

This research was designed to provide s ubstantive information about the operation and capabilities of different resistance ma nagement models to assist the re gulatory expert in its proper use and interpretation of results.

Models Employed

EPA POPGEN (Population Genetics Simulation Mode ls) are a series of generalized, flexible population genetics modeling systems emphasizing the evolution of resi stance to conventional pesticides as well as transgenic crops. The programs are capable of simulating complex landscapes that vary spatially as well as temporally over the course of a growing season.

POPGEN-D The deterministic model POPGEN-D is a simple frequency-based model that incorporates noncompliance, non-random mating in compliant fields and larval movement. The model includes a switch to vary between a block-refuge design and a seed-mixture design. Biological parameters were primarily derived from that model with the exception of mortality rates, block size, dispersal rates and spatial aspects of refuge placement and movement.

POPGEN-S2 The stochastic model POPGEN-S2 is a modificati on of a previous written stochastic, spatially explicit, dual locus simulation model designed to eval uate the population dynamics of crop pest insects relating to the development of resist ance to pesticidal contro ls installed in the cr op. (Caprio et al 2006).

Materials and Methods

Stochastic Model

The stochastic model was a modification of a pr evious western corn rootworm model (Caprio et al 2006). Biological parameters were primarily derived from that model with the exception of mortality rates, developmental delays, block size and dispersal rates. Th ese changes are discussed

below. The stochastic model did not incorporate larval movement a nd was used only for evaluations of4 different block refuge scenarios.

Mortality rates: Mortality ra tes on transgenic corn were adopted from Pioneer's submission. We used two different dominance values ( h ), 0.05 and 0.2. If w ss is the likelihood th at a susceptible individual will emerge from transgenic corn relative to non-Bt corn, and we standardize w rr as 1, then the fitness of heterozygotes (w rs ) is estimated as (w rr – w ss )* h +w ss . When h =0.05, heterozygotes were ca. 5x more likely to survive than susceptibles , and were ca. 20x more likely to survive when h =0.2

when w ss = 0.01. Those rates decrease as the fitness of susceptible homozygotes increase (a reflection of the difficulty in using dominance rather than the w rs :w ss ratio).

Developmental delays: To replicate the developm ental delay of larvae on transgenic plants in the Pioneer model, we added the ability to delay insects by not moving a portion of any age group by a fraction. Normally, all the individuals in an age gr oup would age one day in the model per time step. In the current model, 5% were held back from each age group (though they did experience mortality) and did not age. This resulted in a m ean developmental delay of 2.3 days.

Block size: The simulated universe consiste d of a 40x40 matrix of 30m x 30m blocks, or 144 ha. Four different types of simula tions were run with different field sizes. Because the current refuge strategy requires that refuges be a ssociated with each field, the four different scenarios were field sizes of 5x5 blocks (each block is 30m x 30m, for a total of 2.25 ha), 10x10 blocks (9 ha field size), 20x20 blocks (36 ha field size) and finally a single 40x40 block field of 144 ha. Mean field sizes in the U.S. range from as low as 2.5 ha/field to as high as 150 ha/f ield, with most fields in the 30-50 ha/field range. In each case, a 20% refuge was planted along one e dge of the field. With the 5x5 block fields, this refuge consisted of a single row of 5 blocks along one edge. In the 10x10 bloc k fields, the refuge was 2 rows wide and 10 blocks long. The refuge widths for the remaining scenarios were, of course, 4 and 8 rows wide, maintaining in each case a 20% refuge (by area). Pioneer maintained that a square field is

a best-case scenario (Flexner 2009), but we feel the opposite is true . It is difficult to determine from the Pioneer submission, but it appears that they simulated rectangular fi elds with refuges planted at the short ends of these fields. Most growers, however , plant lengthwise along such fields to minimize the number of turns they must make. It is most common to plant the refuge along one long edge of the field. If dispersal of adults is in question, then as fields become more rect angular (length:width ratio increased), refuges become thinner but the maximu m distance between any Bt-plant and the refuge

decreases. Under these conditions a square field (as we used) maxi mizes the distance between refuge and Bt plants and can be seen as a worst case. Pi oneer's model (or at least my interpretation of their5 model), maximizes the distance between refuge and Bt plants, but is not consiste nt with typical grower practices.

Adult dispersal: Dispersal was modeled based on Pioneer's dispersal data. We combined the data from Pioneer Fig. 3 and built probability dist ributions based on the frequency of beetles moving different distances. In short, if th e total height of the bars in the a bove graphs summed to one, then the individual bar heights would repr esent the likelihood that a beetle would move that distance. One change was made. The Pioneer data includes beetles that did not disp erse. The zero movement class consisted of 45% of the beetles, so the stochastic model assumes that 45% beet les do not disperse (on a daily basis) while the remainder move according to the Nowatski data. The terminology difference may be important. Pioneer's interp retation of the data is that beetle s move an average of 14.5m/day. Using the same data, but only looking at beetles that move, the average movement rate is ca. 29m/day. This is purely a difference based on definition. In th e our stochastic model, each day a beetle that was dispersing would be assigned a random position in the 30x30m block, chose a random direction, determine the distance it moved (based on the Nowats ki data) and move it to a new block. As beetles disperse, they may move only a short distance and not change the 30x30m block where they are currently located.

To determine how well the stochastic model di spersal kernal matched the Nowatski data, a

simulation was run in which one bl ock was populated with adults and the remaining blocks were empty. The number of beetles al ong a transect through the fields wa s counted every 3 days (Table 1, Fig. 1), corrected for the number of cells at a dist ance and multiplied by the distance they had dispersed to estimate the mean dispersal rate (per Nowatski) at different times (Table 2). The dispersal rate estimated after 3 days is 17.6m/day, which approximat es the Nowatski estimate. Given the model uses actual field data to parametrize the dispersal, it was felt that this was a suitable match. Interestingly,

the data from Table 2 demonstrate that the estimate of the mean dispersal rate drops over time. This

occurs because beetles don' t move in a straight line day after day. Brownian diffusion (based on a random walk) is usually assumed for insect disper sal. By confounding measurements over 3 days and using that data in a daily time-step model, Pioneer may have underestimated the daily dispersal rate of western corn rootworm. To demonstrate this proce ss, a simple spreadsheet was constructed in which

10 individuals are simulated that all move 14.5m/day. At the star t of each day a new direction is randomly drawn and the beetles move 14.5m in that di rection. The right-hand column shows estimates of the dispersal rate using the Nowatski method. Only the estimates from the first day accurately6 estimate dispersal rate/day (because all the indivi duals are moving away from the release point). Estimates based on beetles displacement on subseque nt days tend to underestimate daily movement. After three days of movement the estimated dispersal ra te/day is less then half the actual value (because the spreadsheet uses random numbers, values may differ). Unfortunately, because we do not have access to the original data and do not know the captur e distributions for beetles of different ages, we cannot correct for this problem. In retrospect, anal ysis of the data using a diffusion model may have been more appropriate. If we assume that Tabl e 3 can show the extent of underestimation of the dispersal, and that equal numbers of beetles we re caught on days 1,2 and 3, a “ballpark” figure is indicates Nowatski underestimated the daily dispersa l rate by a factor of 0.68. Using a definition of estimating the dispersal as the distance a beetle m oves when it moves (this is purely semantics but critical to define), and correct by th e above factor, the resultant dispersa l rate is ca. 41m/day. It should be noted that, because a different definition of disp ersal was used, this rate is not the same as the 45m/day dispersal rate used in the sensitiv ity analysis in the Pioneer submission.

Refuge movement: Several scenarios were e xplored for yearly refuge deployment. One scenario, used for the transect experiments presen ted below, used refuges that were moved every generation. All growers planted thei r refuge on the west side of th eir fields in one generation and

moved them to the east side in the next. Because the east side of one field is adjacent to the west side of a neighboring field, this scenario resulted in refuges always being planted adjacent to a refuge from the previous year. Ultimately, random placement of refuges each year was chosen as the most likely grower behavior. In this case, at the start of year, the simulation fi rst determines if a field will be compliant, and if so, randomly chooses which side (eas t or west) of the field to plant the refuge. This means that occasionally refuges are not moved, while in other cases they are moved to relatively

isolated areas.

Deterministic model

The deterministic model is a relatively simple frequency-based model that incorporates noncompliance, non-random mating in compliant fields and larval movement. It includes a switch to vary between a block-refuge design and a seed-mixture design. Appendix A includes descriptions of the different parameters. A schematic of the model is shown in fig. 2. The primary difference between the Bt population in the non-compliant population an d the Bt-population in the non-random mating population is that no refuge was planted for the non-compliant fields. Non-compliant populations reduce both non-random mating and re fuge size. As noted in the Pioneer submission, a critical7 component of their modeling is that refuges are moved ever y generation. As mentioned before, while we believe that this may overestimate movement of refuges by growers, this wa s incorporated into the model by defining the refuge size not as the physi cal area devoted to the refuge, but rather the proportion of the eggs in a field th at start out underneath the new refuge. We refer to this as the effective refuge size to distinguish it from physical refuge size (the pr oportion of an area planted to a refuge). These small refuge sizes are independent of non-compliance and result from growers moving their refuge to areas of the field were few eggs were deposited the previous year. To estimate this important parameter, we used the stochastic model to estimate egg densities at the start of a year across our four different field sizes (Fig 3a-e). In each case, we estimated a best and worst scenario. For the best case scenario, the grower moved their refuge to the opposite side of the field, next to where the refuge would have been planted in the neighboring field the previous year (these simulations were run with all growers planting refuges to the same sides of their fields and swappi ng them every year). The worst case scenario was based on plan ting the refuge in the area of th e field where the egg density was the lowest. Because the egg densit y distribution across fields has ne ver been measured, it is unlikely that growers could accomplish this, but it does re present and theoretical lower bound for effective refuge size. The worst case effective refuge sizes fo r the four different field sizes were (from smallest to largest field size): 15%, 11%, 3% and 0.055%. As noted by Pioneer, few insects are moving into the center of large fields as reflected by the effec tive refuge size of 0.055% for 144 ha fields. Other estimates for dispersal are higher and may depend on th e region of the country. To test the effect of higher dispersal rates we doubled the vector of da ily movement (to 29m/day) and repeated the large field size simulation (fig. 1e). In this case the density of eggs in the center of the field was fifty times greater under the worst case scenario and the effective refuge size was 1% . Most fields in the corn belt, however, are smaller, in the rang e of 30-50 ha, and we decided to use a lower limit of 1% for this

parameter.

Asymmetrical movement was included in the model, though the amount movement of susceptible larvae from transgenic corn is uncertain (Hibbard et al 2005). Asymmetrical movement was defined as the probability that larvae that have survived exposure on a Bt-plant will move. When larvae move, they move onto plant genotypes according to the ratio of those genot ypes in the field. In this case, a larvae leaving a transgenic plant had a 95% probability of movi ng onto another transgenic plant. As with other rates in the model, this movement rate is conve rted into a daily rate (depending on the days of larval movement parameter).8 To incorporate our uncertainty in the parameters used in the deterministic model, we developed

distributions based on data for th e lowest and maximum possible valu es we felt the parameter might take, as well as a most likely value. These valu es were used in PERT distributions based on a beta distribution with the most likely value weighted by 4 (Vose 2001). This was done for both the block refuge model (Table 4) and larval movement mode l (Table 5). Each model was run 2000 times, each time drawing parameter values from the relevant PE RT distributions and using the sampled parameters as inputs to the model. The time it took for the resi stance allele to exceed 50% was recorded for each

simulation. We measured the mean and 25 th and 75 th percentiles of the time to resistance. We determined the proportion of simulations in which the time to resistance exceeded 10 years (i.e., 11 years or longer). At a cutoff point of 11 years, the gene has already been in use for 4 years and previous EPA-SAPs have suggested 15 years as a reasonable minima l goal for the lifetime of single gene products.

Uncertainty distributions for the deterministic model.

The uncertainty distributions are divided into those that are common to both models (Tables 4- 5), those specific to the block ref uge model (Table 4) and those specific to the seed-blend model (Table 5).

Common parameters 1. Fitness of susceptible homozygot es on transgenic plants (w ss ). The Pioneer estimate for this parameter was 0.0125. In the past, how ever, as transgenic events have been exposed to a wide variety of environments, w ss has often been found to exceed initial estimates. The minimum value used was 0.0125, the maximum value was 0.08 and the most likely value was 0.04. 2. Dominance (h). Occasionally, dominance associ ated with events that are not a true high-dose is greater than the Pioneer estim ate. Our estimate was centered on the Pioneer estimate but included the possibility that it coul d be 5x less or 5x greater. 3. Initial resistance allele freque ncy (IRF). While many estimates of the initial resistance allele frequency for loci pr oviding resistance to Bt toxi ns are in the range of 0.001, recent studies with western corn rootworm have suggested they may be more common, perhaps in the 0.02 range. These two values were used as our maximum and minimum values, with 0.005 used as the most likely value.

Block refuge parameters 1. Effective refuge size. Without non-compliance, the physical refuge size was always assumed to be 20% of the area. Non-ra ndom oviposition and refuge placement could, however, lead to that 20% of the area having more or less than 20% of eggs initially (see Figs 3a-e). This is a complex parameter that depends on the size and shape of fields, how growers physically place the refuges in the fields, how they move the refuges9 between years and on corn rootworm dispersal rate s. It is likely that this parameter is responsible for many of the di fferences between the results of this model and the model used by Pioneer. The effective refuge size varied from 1% to 20%, with 15% as the most likely value. 2. Non-compliance. Non-compliance must be measured on a field by field basis, and EPA's interpretation that a non-compliance ra te of 30% is an overestimate seems reasonable. There may be localized areas, however, where non-compliance may be greater than 30%, and to incorporate the ef fects of these areas the maximum value of non-compliance was set to 60%. The minimum value was 5% with a most likely value of 10%. 3. Non-compliance mixing. This parameter is re lated to the rate of mixing per generation between compliant and non-compliant fields and may be seen as a measure of persistent non-compliance versus growers that an acci dentally non-compliant in one or more fields. Different monitoring groups suggested 5-7% of growers we re persistently noncompliant, so 93% was chosen as a most lik ely value, 95% as th e maximum and 85% as a conservative minimum value. 4. Non-random mating. This parameter is an estimate of the propor tion of the population emerging in the transgenic fields that ma tes among themselves with no contribution from refuge insects. Based on the simulate d dispersal experiment , it appears as though as many as 50% of these individuals would be unlikely to mate with refuge individuals in the largest field sizes. Because re fuges are moved each year, these non-random populations are merged with the random mating populations for reproduction. The maximum value used was 50%, with a minimum of 5% and a most likely value of 30%.

Seed-blend model 1. Base larval movement. Data suggest there is limited larval movement among nontransgenic plants in the absence on high larval densities (which we find unlikely in a 5% seed-mixture). The range for this paramete r was from 1%, to 5%, with a most likely value of 3% 2. Days of larval movement. Based on Mallet and Porter (1992), larv al movement effects the evolution of resistance by shifting dominance. This sh ift occurs gradually over a number of days. The range of time over wh ich larval movement occurred ranged from 1 to 18 days, with a most likely value of 10 da ys. The mortality rate was assumed to be constant over this period, with the overall mortality rate equal to 1 – w, the relative fitness of the genotype. 3. Asymmetrical movement. Food choice bioassy s in many species suggest that larvae will tend to taste and then move from transg enic tissue at a great er rate than non-Bt tissue. Hibbard et al (2005) found, in a two year study, that th is rate varied in those two seasons between 0 and 30%. Because of the lim ited data on this para meter, the range of this parameter was 0 to 50%, w ith a most likely value of 30%.

Two locus deterministic model

One final set of resistance simulations was conducted for two locus block refuge and seedmixture models. In this case the f itness of the susceptibles and the dominance of resistance at each of two loci were independently determined according to the PERT distributions from the single locus

model. This distribution is informative because it tells us the relative longev ity of two genes stacked10 into the same plant versus using the same two genes se quentially in single gene pl ants. This is relevant because introducing a single gene event into a lands cape with dual gene events may lead to rapid resistance at the first locus and essentially turn all dua l gene events that share th e toxin into single gene

events.

Results and Discussion

Deterministic Model

The distributions of the time to resistance for the block refuge strategy and the larval movement strategy with and without asymmetrical movement ar e presented in figure 4. Because the bars each represent a year, the density of the bars is the proportion of simulati ons that resulted in resistance evolving in that generation. The mean time to resi stance for the block refuge strategy was 11.95 years, and 56.55% of the simulations exceeded 10 years. For the larval movement model, the mean time to resistance was 8.1 years and 13.5% of th e simulations exceeded 10 years.

The results from the two locus block refuge de terministic model (Fig. 5a) show that, though the second toxin has identical uncertainty distributions as the first toxin, the mean expected lifetime of the dual gene product is 132 generations in a block refuge deployment, 6 fo ld longer than if the two toxins had been used sequentially in a bl ock refuge design. Similarly, for the two locus seed-mixture model, the dual gene product lasted a mean of 52.6 years (Fi g. 5b), 3.25 fold longer than if the two toxins had been used sequentially. Single gene products that share toxins with dual gene events can lead to a pseudo sequential use because the sing le gene event selects rapidly fo r resistance to one gene, and once this occurs the dual gene events are ef fectively rendered si ngle gene events.

Stochastic Model

The stochastic model was used to evaluate 4 di fferent sizes and two different dispersal rates (Fig. 6). Because the largest field size took th e entire simulated universe, it was not possible to incorporate non-compliance. To address this issu e, the model was run in an 80x80 cell mode (6400 30x30m cells) for the largest field size, allowing simulation of four 144 ha fields. These larger

simulations for the 144 ha field size we re used in Fig. 6 and in the data analysis. Analysis of variance showed no significant impact of disp ersal rate, but field size was hi ghly significant (p << 0.001). When the simulation was run with a single 144 ha field a nd no non-compliance, the times to resistance were much longer (mean = 3292 d). It was not expected that doubling the disper sal rate would have so little11 impact on the time to resistance, and it suggests the m odel is relatively insensitiv e to this parameter. Using Pioneer's 40 ha field size as a benchmark, the stochastic model predicted an increase in durability of 3.94 fold compared to simulations with no refuge (see discussion). This is much higher than the 1.4 fold reported for this field size in the Pioneer simulations. The simulations presented predict that a block refuge system with 36 ha fields would last 13.9 y, almost double the time predicted by Pioneer. The 144 ha simulations showed much lower durability (8.4 y, or 2.4 fold longer than without a refuge). It may be rele vant that the dimensions of these fields had a 240m wide refuge and a 960 meter wide strip of Bt-corn. Pioneer's simulations may have b een run with a 200m wide refuge and an 800 wide strip of Bt-corn (this occurs in a 40 ha field size in the Pioneer model but a 144 ha field in the current model).

To determine further the sensitivity of the model to the rate of dispersal, the 36 ha field size was run with 2x, 4x and 8x dispersal. Dispersal was ag ain not significant in an ANOVA (p=0.058). The 36 ha field size has only 4 refuges in the default 144 ha system, and this increases variance when in any given year each field has a 30% of being non-compliant. There is a 24% chance that no fields will be non-compliant in any year and a 0.8% that all four fi elds might be non-compliant. In a larger system the mean number of non-compliant fields would rema in more constant over time (Law of Larger Numbers) and there might be a better opportunity to i nvestigate how disp ersal interacts with refuge size and placement. It is possible that dispersal has less ef fect than might have been expected because there

are several different situations that arise in the mo del because of yearly random refuge placement and non-compliance. When a refuge is not moved, lower di spersal rates result in a larger effective refuge. Simulations suggest this frequently delays resistance evolution even though it reduces the probability

that an individual beetle will move and mate at random with beetles emerging in Bt-fields (Caprio 2001). Similarly, if the refuge is moved next to th e location of where a refuge had been planted in a neighboring field in the previous ye ar, lower rates of dispersal might also maximize effective refuge size. In contrast, when refuges ar e moved to areas relativel y distant from the loca tions of the previous years refuges, higher rates of dispersal would maximize the effective refuge size. It is clear, however, that these results suggest that the model is less sensitiv e than expected with respect to dispersal rates.

Discussions over whether the Nowatski data over or underestimate move ment rates for the entire corn belt are, perhaps, less critical than expected.

Discussion12 Pioneer presents in Table 1 (MRID PHI-2009-19 7) the results from their model using the benchmark parameters listed in Table 2. Using th e deterministic model presented here with similar parameter values, similar but not identical durabilit ies were found. In the case of no refuge, Pioneer's model suggested resistance would evolve in 5 years, while this model suggest ed durability would be 3 years. The stochastic model predicted a durability of 3.5y in the absence of re fuges. For comparison's sake, we measure other rates relative to these base line numbers. For the block refuge model, Pioneer suggested resistance might evolve in 7 years, or 1.4x the time wit hout a refuge. Assuming our most likely value of a 15% effective refuge (and otherwis e using Pioneer's benchmark values), we obtain a predicted durability of 14 years, 4.7 fold longer th an the simulation without a refuge. One possible reason for the large discrepancy is the effective re fuge size. To duplicate the 1.4x durability increase noted by Pioneer, the effective refuge size would need to be decreased to 1%. This was the minimum value in our uncertainty distribut ion for this parameter. This parameter was set as the minimum because it would require large fields (on the order of 150 ha), that th e dispersal rates used by Pioneer

were not underestimated and are t ypical of rootworm dispersal rates throughout the corn belt (there is some data to suggest otherwise), and that growers c onsistently move refuges to regions of their fields with the lowest egg density. The ove rall risk profile from our uncerta inty simulations for block refuges is impacted by other parameters that deviated from the Pioneer benchmark values as well. Although relative fitness, initial frequency and dominance used the same uncer tainty distributions for both the block and seed-blend simulations are were likely to affect both models in a similar manner, these parameters could shift both distribu tion in time. The inclusion of the possibility of higher levels of dominance, in particular, could lead to shorter predicted durability es timates. Other parameters specific to the block refuge model could lead to the longer durability predicted in these simulations compared to Pioneer's simulations. The distribution for non-comp liance was weighted more heavily towards lower values (with a most likely value of 10% non-complian ce) than reported by CSPI (Center for Science in the Public Interest), based on EPA's analysis th at report overestimated non-compliance when measured on a field by field basis. We assumed that mo st non-compliant growers were “accidentally” noncompliant, perhaps due to miscalculations or mist akes during planting. Only 5% of non-compliant growers were assumed to be persisten tly non-compliant year after year.

For the larval movement model, the current mode l predicted a durability of 9 years, or a 3-fold increase in durability compared to no refuge. This compares reasonably to the Pioneer estimate of 10 years, or a 2-fold benefit over thei r results run without a refuge. Th e values both compare quite well to the mean of our risk profile, 8.1 years.13 These data suggest there are mi nor differences between the larv al movement results generated by the Pioneer model and the model presented here, but the major differences lie in the block refuge model. Although Pioneer suggests th at non-random mating is important in their results, we did not find

that to be the case in our simulati ons (Table 6). The parameter most likely to be responsible for the

differences between the two models is the effective refuge size. We agree with the Pioneer assertion that it is possible that areas of large fields could have very low egg densities, and should growers manage to locate and plant refuges over those areas , the block refuge strategy could be severely compromised. Understanding of the differences in the results is limited by details missing in the Pioneer data submission. They simulated a size of 4x10 cells (each cell as 100m x100m), or 40 ha. It was unclear if this is a single field or represente d multiple smaller fields, but given their results, a single field is likely. The exact placement of the refuge in the 4x10 cell unive rse is important and not stated. Based on an ESA presentation (Flexner 2009) , it appears as though the refuge was located to one “thin” edge of this system (perhaps 2 4 cell ro ws to one edge), which would maximize the furthest distance between transgenic plants and the refuge, but would seem to be at odds with common grower behavior (this is supported by the fa ct that 20% of the cells by planti ng along the long edge in this size is not possible in this geometry). Again, based on the same ESA presentation (Flexner 2009), it appears as though the only way for insects to reach the plants to the right of the system (when the refuge was planted on the left side) was to m ove across the intervening 800m. In actual field systems there are other fields surrounding the simulate d field with refuges that, on aver age, would be adjacent to the current field, and the mean distance between the furthe st from a refuge and the refuge would be 400m. Most simulations are designed as toru ses (an insect moving off the left side of the system arrives at the right side, as if the system were a donut) to overco me this problem, but the animation (Flexner 2009) did not appear to be designed in this manner. This model minimizes the effects of dispersal and would tend to lower Pioneer's estimate of the durability of a block refuge system. In the current stochastic model, we simulated multiple fields arranged in a large torus. Refuges were placed randomly along one long edge of fields, and non-compliant fields we re also chosen randomly each year. In these simulations, occasionally refuges in adjacent fields would be side by side, while in other cases they would lie at opposite ends. In the latter case there would clearly be an area with low egg densities between the two refuges, but that is apparently offset by refuges that must logically be more closely spaced in other areas. The Pioneer model and our di spersal kernel based on the Pioneer dispersal data ignore the potential of longer dist ance dispersal (Storer et al 2003). Although the contribution from14 such moths is likely to be small, in areas of fiel ds were there are essentia lly no eggs (see Fig 3d), the

contribution of these few indi viduals could be significant.

A number of explanations for the differences in the predicted longevity of the block refuge strategy between the Pioneer simulati ons and those presented here are possible. Important details of the Pioneer model are missing and it is difficult to determine how certain key elements were modeled. The Pioneer model used a single 40 ha block of 100x100m cells arranged in a rectangular fashion (4x10). In discussions, corn entomologists expected most growers would plant lengthwise in such a field, and the refuge would be planted along either one of the long edges (here termed the east and west edges). Simulations of refuges planted to the short ends of these fields (the north and south edges) would not seem to be the normal grower practice and would tend to increase the mean distance between refuges and Bt plants. This choice, if made, would increase the likelihood of the holes in random mating mentioned by Pioneer in their submi ssion and reduce predicted durability of the block refuge design.

Pioneer provided little detail in the spatia l aspects of their simulation. Was a single field simulated with all the refuge gathered in one place, or were smaller fields simulated with multiple refuges? When the refuge was moved each, wher e was it moved to? Was it moved to a random location, or was it always moved to the opposite end of th e field from it's previous location? The latter choice might lead to a significant decrease in the effective refuge size and more rapid resistance development.

How were boundary conditions treated? When an in sect moved from the field to the west side, for example, what happened to it? Did it die, was it reflected back into the system, or was a torus simulated? The choice here may be important as it can again impact the mean distance between refuge and Bt plants and ultimately effective refuge size. As an example, consider a string of 4x10 fields, each 1000m long, each with it's own 200m wide refuge. The maximum distance between any two refuges would be 1600m, and the maximum distance between a refuge and a Bt plant would be 800m. The average distance would be 800m between refuges and 400m between Bt plant fu rthest from a refuge and the nearest refuge. If Pioneer chose not to use a torus design, and simulate d a single field with a refuge at either the west or eas t end, the maximum distance between a Bt plant and the nearest refuge would be doubled (800m). In a simulation of a torus, insects leaving the east edge of the field end up in the west edge, in essence creati ng a long string of identical fields. This design correctly models the15 mean distance, but there is no variance. The simula tions described here included multiple fields with random placement of refuges on a yearly basis to either the east or west edge of a field, which seems to best mimic grower behavior. This design meant th at in some cases the re fuges would be at opposite ends of adjacent fields, while in other cases they might be side by side. This design choice allowed

simulation of both the correct mean distance between refuges, but the variance in this distribution.

Our interpretation of the data in dicates that rather than non-rand om mating being the threat to a block refuge system, it is non-random oviposition and refuge placement that are most important. The conclusion is supported the results of the stochastic model run with just one 144 ha field. The time to resistance in this case was over five years longer than any other scenario, despite the fact that substantial areas of non-random mating must have existed in the center of these blocks (see Fig. 3d). What did occur is the field was always moved adj acent to a previous years refuge, so there was a relatively high density of eggs under the new refuge and the effective refuge size was high. This suggests a potentially plausible method that Pioneer might use to bolster their claims that failure of the block refuge system is imminent. If the relative density of eggs coul d be estimated in transects across fields of different sizes (similar to Fig. 3a-e), and if more accurate estimates of how growers plant and move their refuges over years, then major differences between the two models might be resolved. Such a study might provide an estimate of the contribution of beetles dispersing l ong distance distances.16 References

Caprio, M A, T Nowatzki, BD Siegfried, LJ Mein ke, RJ Wright, and LD Chandler, 2006, Assessing Risk of Resistance to Aerial Applications of Methyl-Parathion in Western Corn Rootworm (Coleoptera: Chrysomelidae), J. Econ. Entomol. 99(2): 483-493.

Caprio, M. A., D. Sumerford & S. Sims. 2006b. Evalua ting transgenic plants for suitability in pest and resistance management programs, 2nd edition. L. Lacey & H. Kaya [eds]. Field Manual of Techniques for the Application and Eval uation of Entomopathogens.

Hibbard, BE, MI Higdon, DP Duran, YM Schweikert, and MR Ellersieck. 2004. Role of egg density on establishment and plant-to-plant movement by western corn rootworm larvae (Coleoptera: Chrysomelidae). J. Econ. Entomol. 97(3): 871-882

Hibbard, BE, TT Vaughn, IO Oyediran, TL Clark, MR Ellersieck. 2005. Effect of Cry3Bb1 expressing transgenic corn on plant-to-plant movement by western corn rootworm larvae (Coleoptera: Chrysomeldiae). J. Econ. Entomol. 98(4): 1126-1138

Mallet, J, P. Porter. 1992 Preventing Insect Adaptati on to Insect-Resistant Crops : Are Seed Mixtures or Refugia the Best Strategy?, Proc R Soc London B 250, 165-169.

Meihls, LN, ML Higdon, BD Siegfried, TA Spencer, NJ Miller, TW Sappington, MR Ellersieck, and BE Hibbard. 2008. Increased survival of western corn rootworm on transgen ic corn within three generations of on-plant greenhouse selec tion. Proc.Nat. Acad. Sci. 105: 19177-19182.

Pan, Z., DW Onstad, BH Stanley, JL Flexner, Q Chen, and SA Lefko, Addendum to: A Spatially Explicit Computer Simulation Study of Corn Rootworm, Diabrotica spp. [Coleoptera: Chrysomelidae], Pioneer Hi-Bred Internationa l, Inc., Study ID: PHI-2009-125

Rood, T, Z Pan, D Onstad, B Stanley, JL Flexner, Q Chen, T Nowatzki, M Ziegler, and S Lefko, Additional Addendum to MRID 473567- 08: A Spatially Explicit Com puter Simulation Study of Corn Rootworm, Diabrotica spp. [Coleoptera: Chrysomelidae], Pioneer Hi-Bred International, Inc., Study Number: PHI-2009-197

Vose, 2000. Risk Analysis: A quantit ative guide. Wiley, New York.17 Figure 1 Beetle Dispersal Time Series 0 30 60 90 120 150 180 210 240 270 300 330 360 0.5

5 50 500 5000

1 4 7 10 13 16 19 

### Beetle Position By Day After Release

1 4 7 10 13 16 19

Distance (m)Number of beetles

Da y P

o s t - Re

l e a s e18 Figure 2. Schematic of the deterministic frequency based model for both seed-blend and block refuge deployments.



## Non-compliant Compliant Bt Refuge Bt

Non-random mating Random mating Selection

Mating

Dispersal19 12345 0 1000 2000 3000 4000 5000 6000

Figure 3a Small 2.25 hectare field size with a 30m row of refuge along the edge of a 120m wide field of transgenic corn. The refuge is always planted along an edge and is switched every year. Previous year’s refuge was in strip 5, while the current one is in strip 1. Growers always plant to the same side of their fields, so the reason there is a larger egg population under strip 1 than in strip two is because there was a refuge to the left in another field. In this case, if growers planted in strip one, the new refuge woul d get 19.8% of the overwintering eggs. How ever, if they planted in the area with the lowest density, the refuge would 15% of the overwintering eggs

1 2 3 4 5 6 7 8 9 10 0 5000 10000 15000 20000 25000Egg Density/900m2

Figure 3b 9 hectare field. The highest density is under the previous years refuge, but continued “drift” from neighboring fields is observed. In this case 18.6% of the eggs would be under the refuge in the default planting pattern, but the refuge would only get 11% if growers planted the refuge in the area with the lowest density of eggs. The refuge in this case was 2 strips (60m) wide.20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 0 5000 10000 15000 20000 25000 30000 35000 40000

Figure 3c 36 hectare field (about an average field for many parts of the U.S.). With the default planting pattern, the refuge would receive 20% of the eggs, but only 3% if growers planted over the area with the lowest density. The refuge in this case was 4 strips wide (120m).

1 2 3 4 5 6 7 8 9 1011 1213 1415 1617 18 1920 2122 2324 25 2627 2829 3031 32 3334 3536 3738 39 40 0 5000 10000 15000 20000 25000

Figure 3d 144 hectare field size. At this size, the default planting would result in 23% of the eggs in the new refuge, but the worst case drops to 0.055% (if growers plant in the middle). The refuge is 8 strips wide (240m).21 4 292 580 868

40 76 112 148

184

220

256 328 364

400

436

472

508

544 616 652

688

724

760

796

832 904 940

976

1012

1048

1084

1120

1156

1192

1228

1264

1300

1336

1372

1408

1444

1480

1516

1552

1588 0 5000 10000 15000 20000 25000 30000 35000 40000

Figure 3e Doubled dispersal distance in a 144 hectare field. This figure contains a column for each 30x30m block rather than an average of 3 transects as in the other four figures

.

Figure 3. Simulated egg densities in transects across Bt-corn and associated refuges

.22 Figure 4a. Block refuge model

Figure 4b. Seed-mixture model

Figure 4. The distributions of the time to resistance23 Figure 5. The distribution of time to resistance from 2000 simulations of a block refuge deployment with two resistance loci. The mean time is 132 generations and 100% of the simulations persisted beyond 26 years.24 Figure 6. Means of four simulations at each field size (given in hectares). Simulations us ed for the dark line used the Nowatski data, simulations used in the red line doubled that dispersal rate25DistanceDays 0 30 60 90 120 150 180 210 240 270 300 330 3601 4604 542 119 20 2 3 0 0 0 0 0 0 04 675 386 182 64 19 12 4 3 3 2 0 0 07 224 169 99 72 37 15 9 2 0 0 4 0 010 110 103 60 57 30 17 13 7 0 0 1 0 113 58 61 39 35 25 9 7 5 2 2 0 0 116 44 31 31 26 28 18 5 0 1 4 3 1 119 20 26 20 13 13 9 4 3 1 1 1 0 0cells 1 8 16 24 32 40 48 56 64 72 80 88 96Table 1.Beetle counts over time in a simulated dispersal experiment. Distance is the distance in meters from the initial release point.DistanceDays 0 30 60 90 120 150 180 210 240 270 300 330 360 Sum m/day1 0 130080 114240 43200 7680 18000 0 0 0 0 0 0 0 313200 27.224 0 92640 174720 138240 72960 72000 34560 35280 46080 38880 0 0 0 705360 17.647 0 40560 95040 155520 142080 90000 77760 23520 0 0 96000 0 0 720480 13.6610 0 24720 57600 123120 115200 102000 112320 82320 0 0 24000 0 34560 675840 11.0913 0 14640 37440 75600 96000 54000 60480 58800 30720 38880 0 0 34560 501120 9.2816 0 7440 29760 56160 107520 108000 43200 0 15360 77760 72000 29040 34560 580800 8.9819 0 6240 19200 28080 49920 54000 34560 35280 15360 19440 24000 0 0 286080 6.81Table 2.Estimate of dispersal rates over time in a simulated dispersal experiment. Each cell is the product of the corrected number of beetles at a distance in timemultiplied by the distance. These are summed, averaged by beetlenumbers and divided by the number of days for the sample to obtain the mean displacement inm/day.26Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 Day 0 m/dayX-position 0 0 0 0 0 0 0 0 0 0Y-position 0 0 0 0 0 0 0 0 0 0Day 1 Day 1 Day 1 Day 1 Day 1 Day 1 Day 1 Day 1 Day 1 Day 1direction 2.43 1.06 3.26 3.6 2.08 3.12 5.85 0.52 1.08 2.01displacement x -10.94 7.03 -14.39 -13.03 -7.03 -14.5 13.14 12.56 6.81 -6.13displacement y 9.51 12.68 -1.77 -6.36 12.68 0.25 -6.13 7.25 12.8 13.14Distance Distance Distance Distance Distance Distance Distance Distance Distance Distancetotal displacement 14.5 14.5 14.5 14.5 14.5 14.5 14.5 14.5 14.5 14.5 14.5Day 2 Day 2 Day 2 Day 2 Day 2 Day 2 Day 2 Day 2 Day 2 Day 2direction 1.95 5.11 0.28 1.48 2.13 4.59 5.52 1.24 4.42 3.89displacement x -16.38 12.7 -0.45 -11.77 -14.71 -16.26 23.57 17.28 2.57 -16.73displacement y 22.96 -0.67 2.23 8.09 24.98 -14.14 -16.2 20.96 -1.06 3.25Distance Distance Distance Distance Distance Distance Distance Distance Distance Distancetotal displacement 28.2 12.71 2.28 14.28 28.99 21.55 28.6 27.16 2.78 17.05 9.18Day 3 Day 3 Day 3 Day 3 Day 3 Day 3 Day 3 Day 3 Day 3 Day 3direction 0.99 1.43 6.06 1.88 5.86 5.41 6.25 1.48 0.91 4.29displacement x -8.48 14.71 13.67 -16.25 -1.47 -6.94 38.06 18.54 11.5 -22.63displacement y 35.12 13.69 -1.03 21.88 19.08 -25.25 -16.71 35.4 10.36 -9.99Distance Distance Distance Distance Distance Distance Distance Distance Distance Distance36.13 20.1 13.71 27.25 19.14 26.18 41.57 39.97 15.48 24.74 8.81Day 4 Day 4 Day 4 Day 4 Day 4 Day 4 Day 4 Day 4 Day 4 Day 4direction 3.05 1.13 5.24 1.71 1.73 5.9 0.45 0.26 5.69 1.85displacement x -22.92 20.84 20.92 -18.27 -3.74 6.5 51.1 32.55 23.52 -26.63displacement y 36.38 26.83 -13.59 36.24 33.4 -30.68 -10.35 39.16 2.25 3.94Distance Distance Distance Distance Distance Distance Distance Distance Distance Distancetotal displacement 43 33.98 24.95 40.58 33.61 31.36 52.13 50.92 23.62 26.92 9.03Table 3.Example of estimate daily dispersal rate over time from insects moving in a random walk at 14.5m/day. Direction is given in radians. The meandisplacement (m/day) is an average from 10 individuals andwill chance depending on the random directions drawn.27 Min Most Likely Max Wss 0.01250 0.04000 0.08000 Dominance 0.01000 0.05000 0.20000 IRF 0.00100 0.00500 0.02000 Effective refuge 0.01000 0.15000 0.20000 Non-compliance 0.05000 0.10000 0.60000 Non-compliance mixing 0.85000 0.93000 0.95000 Non-random mating 0.05000 0.30000 0.50000

Table 4. Uncertainty parameters used for the block refuge simulations

Min Most Likely Max Wss 0.01250 0.04000 0.08000 Dominance 0.01000 0.05000 0.20000 IRF 0.00100 0.00500 0.02000 Refuge 5.00000 5.00000 5.00000 Base larval movement 0.01000 0.03000 0.05000 Days of larval movement 1.00000 10.00000 18.00000 Asymmetrical movement 0.00000 0.30000 0.50000

Table 5. Uncertainty distributions used for the seed-mixture simulations28 Coefficients: Es timate Std. Error t va lue Pr(>|t|) (Intercept) 11.6958 2.2321 5.240 1.78e-07 *** dominance -84.9810 1.3097 -64.887 < 2e-16 *** IRF -629.7121 14.7731 -42.626 < 2e-16 *** Non.compliance -8.8690 0.4426 -20.037 < 2e-16 *** non.compliance.mixing 1.5415 2.4028 0.642 0.521 non.random.mating -0.7173 0.5336 -1.344 0.179 Refuge 54.6024 1.2589 43.374 < 2e-16 *** wss 61.5433 3.4862 17.653 < 2e-16 ***

--- Signif. codes: 0 '***' 0.001 '* *' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.797 on 1992 degrees of freedom Multiple R-squared: 0.8095, Adjusted R-squared: 0.8088 F-statistic: 1209 on 7 and 1992 DF, p-value: < 2.2e-16

Table 6. Regression coefficients for the regressi on on the results from the block refuge model.

Coefficients: Estimate Std. Error t va lue Pr(>|t|) (Intercept) 13.25428 0.20471 64.748 <2e-16 *** BaseLarvalMovement 7.43119 4.18984 1.774 0.0763 . dominance -63.84439 0.91981 -69.410 <2e-16 *** IRF -466.03593 10.13300 -45.992 <2e-16 *** SelectionDays 0.09145 0.01017 8.992 <2e-16 *** wss 57.75906 2.49046 23.192 <2e-16 ***

--- Signif. codes: 0 '***' 0.001 '* *' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Table 7. Regression coefficients for the results from the seed-mixture model with asymmetrical movement.29 Appendix A. Deterministic parameter description

Wss – The fitness of homozygous susceptible genotypes on the transgenic plants. Pioneer used a value of 0.0125. Without better knowledge of the Pioneer field tests it's difficult to develop a range, but a specific factor over and unde r was used to deve;op our estimat e. The sorts of questions to be asked might be what is the lowest/highest performance under field conditions. Instead of 98.75% mortality, could it be 95% (as an example)?

Dominance – The relative fitness of heterozygotes. I use the traditional population genetics equation: Wrs = (Wrr – Wss) * h + Wss

where h is the dominance and Wrr = 1.0. Our distri bution is derived from Pi oneer's estimate of h=0.05, but that seems low compared to many estimates ba sed on LC50s. LC50s, however, are based on a log scale, however this research has used a linear scale. In general, for products with lower effectiveness, dominance tends to be higher (see Caprio et al. 2006b).

IRF – The initial resistance allele frequency. My im pression from dealing with CRW folks is that this value has tended to be higher in empirical estimates th an we have seen for othe r insects. I believe I used Pioneer's default value as my most likely value.

Block refuge model parameterization:

Values cited above are shared with the larval movement model and changes in would effect both the block and larval movement models to a similar degree.

Refuge proportion – Initially set this to 20%. However, it can be viewed as the proportion of eggs in a compliant field that end up bei ng in a refuge (perhaps after being moved), and the non-random oviposition patterns can be addressed using this parameter. If more i ndividuals survive in refuges, and they tend to lay eggs near that area, then the number of eggs in the soil to start the new growing season will be greatest near the refuge from the previous y ear. If the refuge is moved, it stands to reason the new refuge, even if it covers 20% of the area, will receive less than 20% of the eggs oviposited in the field. Based on this modeling research, this f act and the likelihood that refuges are moved are identified as key assumptions of the Pioneer model. Unfortunately, there is little empirical data on the non-random oviposition patterns of CRW.

Non-compliance – This is the proportion of fields that are planted without a refuge. Individuals emerging from this area mate as a separate populati on. The most likely value is a literature value (CSPI). This parameter should be based on a fiel d by field basis and incorporate some estimate of amount of non-compliance. A person planting a 3% re fuge in a field is not the same as someone not planting any refuge. Alternatively, it shouldn't measure if a field is compliant (as yes./no answer), but how close to compliant a non-compliant field is. General comments have been encountered which indicated that 30% non-compliance was too low. It is not obvious how non-compliance was estimated. This parameter may be important because it dete rmines how much Bt-corn is planted without any refuge.

Non-compliance mixing – This determines how much movement there is between compliant and noncompliant fields. If non-compliant growers are accid entally non-compliant (i.e. they never got around to planting the refuge because of weather), then they will generally be compliant in the next year, and this parameter should be close to 1.0. On the ot her hand, if non-compliant growers are growers who30 purposefully don't plant refuges and c ontinue this practice year-by-year, then this parameter should be low. How low would depend on how much move ment you think naturally occurs between noncompliant and compliant fields. Inevitably this l eads to discussion of the clustering of non-compliant fields. Matt Carrol has done some work on this with the Monsanto model.

Non-random mating – This paramete r addresses the proportion of th e population in the Bt-portion of compliant fields that is far enough away from the ref uges that they do not mate at random. These differ from non-compliant fields in that refuge is planted, but it was too far away. The non-random oviposition that results from this non-random mating (and limited dispersal capability) is covered by the reductions in refuge size discussed above.

Larval movement model parameters

The refuge proportion is fixed at 0.05 in the larval movement model. This was based on Pioneer's last data that indicated they had reduced the variance considerably.

Base larval movement – This is the amount of movement that might occur in a non-transgenic field with populations low enough so that significant density -dependent dispersal is not a factor. In the stochastic models, population declin e is found to be a greater problem than density-dependent factors such as mortality or dispersal. The model does in clude asymmetrical movement, with the movement of individuals on Bt plants increased by an amount related to their daily mortality rates. Susceptibles therefore move more than heterozygotes, and resist ant homozygotes would only move at the base larval movement rate. Of course, an individual leaving a transgenic plant only ha s a 5% change of moving onto a non-Bt plant. Various estimat es of larval dispersal rates have been encountered but mostly on the low side. Generally, this research assumes that populations will be low in seed-mix fields and there will be little to no density-dependent larval dispersal. This could be an area we need help with, or we might find by selecting alternate numbers is relativel y unimportant compared to other parameters and a precise estimate would be a misapplication of resources.

Days of larval movement – Movement and larval de velopment are assumed to occur over this period of days. Movement occurs once per day. Daily survivorsh ip is constant and is determined as the rate that raised to the days power results in the overall fitness listed for the genotype in the original table. In the case of Pioneer, that was 0.0125 for susceptible larvae. If the days of larval movement is 10, then the daily survivorship rate is 0.645. This means that on day 1, 36.5% of the SS larvae on Bt die (I assume conservatively that mortality occurs before di spersal) , and of the remaining 64.5%, 36.5% will disperse (see above). Of these dispersants, on ly 5% will move onto non-Bt plants, increasing the refuge size from 5% to ca. 6%. If the base la rval movement is low (remember it is spread over days as well), then this process can lead to an increase in the refuge population and act ually decrease the rate of resistance evolution. I put this in to the model just to see what ef fect it might have. The value was varied between 5, 10 and 15 and it had very little effect on the results.