EPA/600/R-10/055
May 2010
Simulation Models Evaluation of Pest
Resistance Development to Refuge in the Bag
Concepts Related to Pioneer Submission
Letter Report
By
Michael Caprio1, John A. Glaser2
^isssissippi State University, Department of Entomology and Plant Pathology, Mississippi State, MS
39762-9775, 2US Environmental Protection Agency, National Risk Management Research Laboratory,
Sustainable Technology Division, 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 Caprio
-------
Introduction
The USEPA, under its administration of the Federal Insecticide Fungicide and Rodenticide Act
(FIFRA), requires the registration of all pesticides and pesticidal materials. The GM crops containing
pesticidal traits are subject to FIFRA registration requirements. Concerns relating to future
environmental effects of these crops can be investigated 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 broad spectrum pesticides that have
significant environmental footprints.
Background
Resistance management (RM) simulation models are designed as deterministic, stochastic, and
spatially explicit analytical methodologies. Simulation models can provide a realistic assessment of the
risk of resistance evolution given the allele frequencies used as initial conditions. New concepts of PIP
crop deployment such as refuge-in-the-bag require a detailed 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 International, Inc. This research specifically analyzed the
relative risks related to the control of western corn rootworm in block refuge deployment and 5% seed-
mixture scenarios.
Research Design
To understand the relative risks related to the evolution of resistance in western corn rootworm
between the currently mandated block refuge deployment and a 5% seed-mixture deployment as
submitted Pioneer Hi-Bred International, Inc., two different models were developed and evaluated. The
two models consisted of a modification of a spatially-explicit stochastic model (Caprio et al. 2006) and
a simpler, frequency-based deterministic model. The latter could be run in a single simulation mode
with a graphical interface to enter parameters or in a risk-assessment mode capable of running
thousands of simulations to estimate the effects of parameter uncertainty. Realistic estimates for
-------
several parameters using in the deterministic model were developed through use of the stochastic
model.
Impacts and Outcomes
This research was designed to provide substantive information about the operation and
capabilities of different resistance management models to assist the regulatory expert in its proper use
and interpretation of results.
Models Employed
EPAPOPGEN (Population Genetics Simulation Models) are a series of generalized, flexible
population genetics modeling systems emphasizing the evolution of resistance 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 non-
compliance, 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 modification of a previous written stochastic, spatially
explicit, dual locus simulation model designed to evaluate the population dynamics of crop pest insects
relating to the development of resistance to pesticidal controls installed in the crop. (Caprio et al 2006).
Materials and Methods
Stochastic Model
The stochastic model was a modification of a previous 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. These changes are discussed
below. The stochastic model did not incorporate larval movement and was used only for evaluations of
-------
different block refuge scenarios.
Mortality rates: Mortality rates on transgenic corn were adopted from Pioneer's submission.
We used two different dominance values (//), 0.05 and 0.2. If wss is the likelihood that a susceptible
individual will emerge from transgenic corn relative to non-Bt corn, and we standardize Wn- as 1, then
the fitness of heterozygotes (wrs) is estimated as (WIT- wss)*/7+wss. When /z=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 wss = 0.01. Those rates decrease as the fitness of susceptible homozygotes increase (a reflection
of the difficulty in using dominance rather than the wrs:wss ratio).
Developmental delays: To replicate the developmental 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 group 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 mean developmental delay of 2.3 days.
Block size: The simulated universe consisted of a 40x40 matrix of 30m x 30m blocks, or 144
ha. Four different types of simulations were run with different field sizes. Because the current refuge
strategy requires that refuges be associated 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/field, with most fields in the 30-50 ha/field range.
In each case, a 20% refuge was planted along one edge of the field. With the 5x5 block fields, this
refuge consisted of a single row of 5 blocks along one edge. In the 10x10 block 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 fields 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 rectangular (length:width ratio
increased), refuges become thinner but the maximum distance between any Bt-plant and the refuge
decreases. Under these conditions a square field (as we used) maximizes the distance between refuge
and Bt plants and can be seen as a worst case. Pioneer's model (or at least my interpretation of their
4
-------
model), maximizes the distance between refuge and Bt plants, but is not consistent 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 distributions based on the frequency of beetles moving
different distances. In short, if the total height of the bars in the above graphs summed to one, then the
individual bar heights would represent the likelihood that a beetle would move that distance. One
change was made. The Pioneer data includes beetles that did not disperse. The zero movement class
consisted of 45% of the beetles, so the stochastic model assumes that 45% beetles do not disperse (on a
daily basis) while the remainder move according to the Nowatski data. The terminology difference
may be important. Pioneer's interpretation of the data is that beetles 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 the 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 Nowatski 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 dispersal kernal matched the Nowatski data, a
simulation was run in which one block was populated with adults and the remaining blocks were
empty. The number of beetles along a transect through the fields was counted every 3 days (Table 1,
Fig. 1), corrected for the number of cells at a distance 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 approximates 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 dispersal. 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 process, a simple spreadsheet was constructed in which
10 individuals are simulated that all move 14.5m/day. At the start of each day a new direction is
randomly drawn and the beetles move 14.5m in that direction. The right-hand column shows estimates
of the dispersal rate using the Nowatski method. Only the estimates from the first day accurately
-------
estimate dispersal rate/day (because all the individuals are moving away from the release point).
Estimates based on beetles displacement on subsequent days tend to underestimate daily movement.
After three days of movement the estimated dispersal rate/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 capture distributions for beetles of different ages, we
cannot correct for this problem. In retrospect, analysis of the data using a diffusion model may have
been more appropriate. If we assume that Table 3 can show the extent of underestimation of the
dispersal, and that equal numbers of beetles were caught on days 1,2 and 3, a "ballpark" figure is
indicates Nowatski underestimated the daily dispersal rate by a factor of 0.68. Using a definition of
estimating the dispersal as the distance a beetle moves when it moves (this is purely semantics but
critical to define), and correct by the above factor, the resultant dispersal rate is ca. 41m/day. It should
be noted that, because a different definition of dispersal was used, this rate is not the same as the
45m/day dispersal rate used in the sensitivity analysis in the Pioneer submission.
Refuge movement: Several scenarios were explored for yearly refuge deployment. One
scenario, used for the transect experiments presented below, used refuges that were moved every
generation. All growers planted their refuge on the west side of their 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 first determines if a field will be
compliant, and if so, randomly chooses which side (east 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 non-
compliance, 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 and 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 refuge size. As noted in the Pioneer submission, a critical
-------
component of their modeling is that refuges are moved every generation. As mentioned before, while
we believe that this may overestimate movement of refuges by growers, this was incorporated into the
model by defining the refuge size not as the physical area devoted to the refuge, but rather the
proportion of the eggs in a field that start out underneath the new refuge. We refer to this as the
effective refuge size to distinguish it from physical refuge size (the proportion 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 swapping them every year). The
worst case scenario was based on planting the refuge in the area of the field where the egg density was
the lowest. Because the egg density distribution across fields has never been measured, it is unlikely
that growers could accomplish this, but it does represent and theoretical lower bound for effective
refuge size. The worst case effective refuge sizes for 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 effective refuge size of 0.055% for 144 ha fields. Other
estimates for dispersal are higher and may depend on the region of the country. To test the effect of
higher dispersal rates we doubled the vector of daily movement (to 29m/day) and repeated the large
field size simulation (fig. le). 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 range 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 genotypes in the field. In
this case, a larvae leaving a transgenic plant had a 95% probability of moving onto another transgenic
plant. As with other rates in the model, this movement rate is converted into a daily rate (depending on
the days of larval movement parameter).
-------
To incorporate our uncertainty in the parameters used in the deterministic model, we developed
distributions based on data for the lowest and maximum possible values we felt the parameter might
take, as well as a most likely value. These values 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 model (Table 5). Each model was run 2000 times, each
time drawing parameter values from the relevant PERT distributions and using the sampled parameters
as inputs to the model. The time it took for the resistance allele to exceed 50% was recorded for each
simulation. We measured the mean and 25th and 75th 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 minimal 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 refuge model (Table 4) and those specific to the seed-blend model (Table
5).
Common parameters
1. Fitness of susceptible homozygotes on transgenic plants (wss). The Pioneer estimate for
this parameter was 0.0125. In the past, however, as transgenic events have been
exposed to a wide variety of environments, wss 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 associated with events that are not a true
high-dose is greater than the Pioneer estimate. Our estimate was centered on the
Pioneer estimate but included the possibility that it could be 5x less or 5x greater.
3. Initial resistance allele frequency (IRF). While many estimates of the initial resistance
allele frequency for loci providing resistance to Bt toxins 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-random 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 refuges
-------
between years and on corn rootworm dispersal rates. It is likely that this parameter is
responsible for many of the differences 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 rate 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 effects 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 related 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 accidentally non-compliant in one or more
fields. Different monitoring groups suggested 5-7% of growers were persistently non-
compliant, so 93% was chosen as a most likely value, 95% as the maximum and 85% as
a conservative minimum value.
4. Non-random mating. This parameter is an estimate of the proportion of the population
emerging in the transgenic fields that mates among themselves with no contribution
from refuge insects. Based on the simulated 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 refuges 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 non-
transgenic plants in the absence on high larval densities (which we find unlikely in a 5%
seed-mixture). The range for this parameter was from 1%, to 5%, with a most likely
value of 3%
2. Days of larval movement. Based on Mallet and Porter (1992), larval movement effects
the evolution of resistance by shifting dominance. This shift occurs gradually over a
number of days. The range of time over which larval movement occurred ranged from 1
to 18 days, with a most likely value of 10 days. 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 bioassys in many species suggest that larvae
will tend to taste and then move from transgenic tissue at a greater rate than non-Bt
tissue. Hibbard et al (2005) found, in a two year study, that this rate varied in those two
seasons between 0 and 30%. Because of the limited data on this parameter, the range of
this parameter was 0 to 50%, with a most likely value of 30%.
Two locus deterministic model
One final set of resistance simulations was conducted for two locus block refuge and seed-
mixture models. In this case the fitness 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 longevity of two genes stacked
-------
into the same plant versus using the same two genes sequentially in single gene plants. This is relevant
because introducing a single gene event into a landscape with dual gene events may lead to rapid
resistance at the first locus and essentially turn all dual gene events that share the 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 are presented in figure 4. Because the bars each
represent a year, the density of the bars is the proportion of simulations that resulted in resistance
evolving in that generation. The mean time to resistance 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 the simulations exceeded 10 years.
The results from the two locus block refuge deterministic 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 fold longer than if the two toxins
had been used sequentially in a block refuge design. Similarly, for the two locus seed-mixture model,
the dual gene product lasted a mean of 52.6 years (Fig. 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 single gene event selects rapidly for resistance to one gene, and once
this occurs the dual gene events are effectively rendered single gene events.
Stochastic Model
The stochastic model was used to evaluate 4 different sizes and two different dispersal rates
(Fig. 6). Because the largest field size took the entire simulated universe, it was not possible to
incorporate non-compliance. To address this issue, 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 were used in Fig. 6 and in the data analysis. Analysis of variance
showed no significant impact of dispersal rate, but field size was highly significant (p « 0.001). When
the simulation was run with a single 144 ha field and no non-compliance, the times to resistance were
much longer (mean = 3292 d). It was not expected that doubling the dispersal rate would have so little
10
-------
impact on the time to resistance, and it suggests the model is relatively insensitive 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 relevant 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 been 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 again 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 fields might be non-compliant. In a larger system
the mean number of non-compliant fields would remain more constant over time (Law of Larger
Numbers) and there might be a better opportunity to investigate how dispersal interacts with refuge size
and placement. It is possible that dispersal has less effect than might have been expected because there
are several different situations that arise in the model because of yearly random refuge placement and
non-compliance. When a refuge is not moved, lower dispersal 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 the location of where a refuge had been planted in a
neighboring field in the previous year, lower rates of dispersal might also maximize effective refuge
size. In contrast, when refuges are moved to areas relatively distant from the locations 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 sensitive than expected with respect to dispersal rates.
Discussions over whether the Nowatski data over or underestimate movement rates for the entire corn
belt are, perhaps, less critical than expected.
Discussion
11
-------
Pioneer presents in Table 1 (MRID PHI-2009-197) the results from their model using the
benchmark parameters listed in Table 2. Using the deterministic model presented here with similar
parameter values, similar but not identical durabilities were found. In the case of no refuge, Pioneer's
model suggested resistance would evolve in 5 years, while this model suggested durability would be 3
years. The stochastic model predicted a durability of 3.5y in the absence of refuges. For comparison's
sake, we measure other rates relative to these baseline numbers. For the block refuge model, Pioneer
suggested resistance might evolve in 7 years, or 1.4x the time without a refuge. Assuming our most
likely value of a 15% effective refuge (and otherwise using Pioneer's benchmark values), we obtain a
predicted durability of 14 years, 4.7 fold longer than the simulation without a refuge. One possible
reason for the large discrepancy is the effective refuge 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 distribution for this parameter. This parameter was set as the minimum
because it would require large fields (on the order of 150 ha), that the dispersal rates used by Pioneer
were not underestimated and are typical of rootworm dispersal rates throughout the corn belt (there is
some data to suggest otherwise), and that growers consistently move refuges to regions of their fields
with the lowest egg density. The overall risk profile from our uncertainty 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 uncertainty 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 distribution in time. The inclusion of the possibility of higher levels of
dominance, in particular, could lead to shorter predicted durability estimates. 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-compliance was weighted more heavily towards lower
values (with a most likely value of 10% non-compliance) than reported by CSPI (Center for Science in
the Public Interest), based on EPA's analysis that report overestimated non-compliance when measured
on a field by field basis. We assumed that most non-compliant growers were "accidentally" non-
compliant, perhaps due to miscalculations or mistakes during planting. Only 5% of non-compliant
growers were assumed to be persistently non-compliant year after year.
For the larval movement model, the current model 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 their results run without a refuge. The values both compare quite well to
the mean of our risk profile, 8.1 years.
12
-------
These data suggest there are minor differences between the larval 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 that non-random mating is important in their results, we did not find
that to be the case in our simulations (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 xlOOm), or 40 ha. It
was unclear if this is a single field or represented multiple smaller fields, but given their results, a
single field is likely. The exact placement of the refuge in the 4x10 cell universe 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 rows 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 fact that 20% of the cells by planting 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 move across the intervening 800m. In actual field systems there are
other fields surrounding the simulated field with refuges that, on average, would be adjacent to the
current field, and the mean distance between the furthest from a refuge and the refuge would be 400m.
Most simulations are designed as toruses (an insect moving off the left side of the system arrives at the
right side, as if the system were a donut) to overcome 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 were 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 dispersal kernel based on the Pioneer dispersal data
ignore the potential of longer distance dispersal (Storer et al 2003). Although the contribution from
13
-------
such moths is likely to be small, in areas of fields were there are essentially no eggs (see Fig 3d), the
contribution of these few individuals could be significant.
A number of explanations for the differences in the predicted longevity of the block refuge
strategy between the Pioneer simulations 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 lOOxlOOm 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 submission and reduce predicted durability of the block
refuge design.
Pioneer provided little detail in the spatial 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, where was it moved to? Was it moved to a random
location, or was it always moved to the opposite end of the 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 insect 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 furthest from a refuge
and the nearest refuge. If Pioneer chose not to use a torus design, and simulated a single field with a
refuge at either the west or east 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 creating a long string of identical fields. This design correctly models the
14
-------
mean distance, but there is no variance. The simulations 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 that in some cases the refuges 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 indicates that rather than non-random 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 adjacent 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 could 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 long distance distances.
15
-------
References
Caprio, M A, T Nowatzki, BD Siegfried, LJ Meinke, 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. Evaluating 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 Evaluation 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 Cry3Bbl 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 Adaptation 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 transgenic corn within three
generations of on-plant greenhouse selection. Proc.Nat. Acad. Sci. 105: 19177-19182.
Pan, Z., DW Onstad, BH Stanley, JL Flexner, Q Chen, and SALefko, Addendum to: A Spatially
Explicit Computer Simulation Study of Corn Rootworm, Diabrotica spp. [Coleoptera: Chrysomelidae],
Pioneer Hi-Bred International, 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 Computer Simulation Study of Corn
Rootworm, Diabrotica spp. [Coleoptera: Chrysomelidae], Pioneer Hi-Bred International, Inc., Study
Number: PHI-2009-197
Vose, 2000. Risk Analysis: A quantitative guide. Wiley, New York.
16
-------
Beetle Position By Day After Release
JD
M—
O
JD
E
5000
500
50
0.5
0 30 60 90 120 150 180 210 240 270 300 330 360
Distance (m)
Figure 1 Beetle Dispersal Time Series
17
-------
Non-compliant
Bt
Bt
Compliant
Refu<
Figure 2. Schematic of the deterministic frequency based model for both seed-blend and block refuge deployments.
18
-------
6000
Figure 3a Small 2.25 hectare field size with a 30m row of refuge along the edge of a 120m wide field oftransgenic 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 would get 19.8% of the overwintering eggs. However, if they planted in the area with the
lowest density, the refuge would 15% of the overwintering eggs
25000
7
8
10
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.
19
-------
40000
35000
30000
25000
20000
15000
10000
5000
10 11 12 13 14 15 16 17 18 19 20
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).
25000
20000
15000
10000
5000
111
.ll
1 23456789 10111213141516171819202122232425262728293031323334353637383940
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).
20
-------
40000
35000
30000
25000
20000
15000
10000
5000
40 112 184 256 328 400 472 544 616 688 760 832 904 976 10481120119212641336140814801552
4 76 148 220 292 364 436 508 580 652 724 796 868 940 101210841156122813001372144415161588
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.
21
-------
S -
Generations
Figure 4a. Block refuge model
~ -,
O.I ' TE !^r
OOOOOOOOOO
20
ne ratio ns
Figure 4b. Seed-mixture model
Figure 4. The distributions of the time to resistance
22
-------
100
200
Dataset$time
300
400
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.
23
-------
Plot of Means
o
o
10
LU
CO
F
O
o
in
o
o
o
225
9
FieldSizec-$ D ispe rsa
36
Field Size
144
Figure 6. Means of four simulations at each field size (given in hectares). Simulations used for the dark line used the
Nowatski data, simulations used in the red line doubled that dispersal rate
24
-------
Days
cells
0
4604
675
224
110
58
44
20
1
30
542
386
169
103
61
31
26
8
60
119
182
99
60
39
31
20
16
90
20
64
72
57
35
26
13
24
120
2
19
37
30
25
28
13
32
150
3
12
15
17
9
18
9
40
Distance
180
0
4
9
13
7
5
4
48
210
0
3
2
7
5
0
3
56
240
0
3
0
0
2
1
1
64
270
0
2
0
0
2
4
1
72
300
0
0
4
1
0
3
1
80
330
0
0
0
0
0
1
0
88
360
0
0
0
1
1
1
0
96
Table 1. Beetle counts over time in a simulated dispersal experiment. Distance is the distance in meters from the initial release point.
Days
0
0
0
0
0
0
0
0
30
130080
92640
40560
24720
14640
7440
6240
60
114240
174720
95040
57600
37440
29760
19200
90
43200
138240
155520
123120
75600
56160
28080
120
7680
72960
142080
115200
96000
107520
49920
Distance
150 180
18000
72000
90000
102000
54000
108000
54000
0
34560
77760
112320
60480
43200
34560
210
0
35280
23520
82320
58800
0
35280
240
0
46080
0
0
30720
15360
15360
270
0
38880
0
0
38880
77760
19440
300
0
0
96000
24000
0
72000
24000
330
0
0
0
0
0
29040
0
360
0
0
0
34560
34560
34560
0
Sum
313200H
705360
720480
675840
50112oM
580800
286080
n/day
27.22
17.64
13.66
11.09
9.28
8.98
6.81
Table 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 time
multiplied by the distance. These are summed, averaged by beetle numbers and divided by the number of days for the sample to obtain the mean displacement in
m/day.
25
-------
Day 0
Day 0
Day 0
Day 0
Day 0
Day 0
Day 0
Day 0
Day 0
Day 0
m/day
X-position
Y-position
direction
displacement x
displacement y
total displacement
direction
displacement x
displacement y
total displacement
direction
displacement x
displacement y
direction
displacement x
displacement y
total displacement
0
0
Day 1
2.43
-10.94
9.51
Distance
14.5
Day 2
1.95
-16.38
22.96
Distance
28.2
Day 3
0.99
-8.48
35.12
Distance
36.13
Day 4
3.05
-22.92
36.38
Distance
43
0
0
Day 1
1.06
7.03
12.68
Distance
14.5
Day 2
5.11
12.7
-0.67
Distance
12.71
Day 3
1.43
14.71
13.69
Distance
20.1
Day 4
1.13
20.84
26.83
Distance
33.98
0
0
Day 1
3.26
-14.39
-1.77
Distance
14.5
Day 2
0.28
-0.45
2.23
Distance
2.28
Day 3
6.06
13.67
-1.03
Distance
13.71
Day 4
5.24
20.92
-13.59
Distance
24.95
0
0
Day 1
3.6
-13.03
-6.36
Distance
14.5
Day 2
1.48
-11.77
8.09
Distance
14.28
Day 3
1.88
-16.25
21.88
Distance
27.25
Day 4
1.71
-18.27
36.24
Distance
40.58
0
0
Day 1
2.08
-7.03
12.68
Distance
14.5
Day 2
2.13
-14.71
24.98
Distance
28.99
Day 3
5.86
-1.47
19.08
Distance
19.14
Day 4
1.73
-3.74
33.4
Distance
33.61
0
0
Day 1
3.12
-14.5
0.25
Distance
14.5
Day 2
4.59
-16.26
-14.14
Distance
21.55
Day 3
5.41
-6.94
-25.25
Distance
26.18
Day 4
5.9
6.5
-30.68
Distance
31.36
0
0
Day 1
5.85
13.14
-6.13
Distance
14.5
Day 2
5.52
23.57
-16.2
Distance
28.6
Day 3
6.25
38.06
-16.71
Distance
41.57
Day 4
0.45
51.1
-10.35
Distance
52.13
0
0
Day 1
0.52
12.56
7.25
Distance
14.5
Day 2
1.24
17.28
20.96
Distance
27.16
Day 3
1.48
18.54
35.4
Distance
39.97
Day 4
0.26
32.55
39.16
Distance
50.92
0
0
Day 1
1.08
6.81
12.8
Distance
14.5
Day 2
4.42
2.57
-1.06
Distance
2.78
Day 3
0.91
11.5
10.36
Distance
15.48
Day 4
5.69
23.52
2.25
Distance
23.62
0
0
Day 1
2.01
-6.13
13.14
Distance
14.5
Day 2
3.89
-16.73
3.25
Distance
17.05
Day 3
4.29
-22.63
-9.99
Distance
24.74
Day 4
1.85
-26.63
3.94
Distance
26.92
14.5
9.18
9.03
Table 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 mean
displacement (m/day) is an average from 10 individuals and will chance depending on the random directions drawn.
26
-------
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 simulations
27
-------
Coefficients: Estimate Std. Error lvalue 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 regression on the results from the block refuge model.
Coefficients: Estimate Std. Error lvalue
(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.
28
-------
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 under was used to deve;op our estimate. 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 distribution is derived from Pioneer's estimate of h=0.05,
but that seems low compared to many estimates based on LCSOs. LCSOs, 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 impression from dealing with CRW folks is that this
value has tended to be higher in empirical estimates than we have seen for other 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 being in a refuge (perhaps after being moved), and the non-random
oviposition patterns can be addressed using this parameter. If more individuals 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 year. 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 fact 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 population. The most likely value is a literature value
(CSPI). This parameter should be based on a field by field basis and incorporate some estimate of
amount of non-compliance. A person planting a 3% refuge 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 determines how much Bt-corn is planted without any
refuge.
Non-compliance mixing - This determines how much movement there is between compliant and non-
compliant fields. If non-compliant growers are accidentally 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 other hand, if non-compliant growers are growers who
29
-------
purposefully don't plant refuges and continue this practice year-by-year, then this parameter should be
low. How low would depend on how much movement you think naturally occurs between non-
compliant and compliant fields. Inevitably this leads 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 parameter addresses the proportion of the population in the Bt-portion of
compliant fields that is far enough away from the refuges 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 decline is found to be a greater problem than density-dependent factors
such as mortality or dispersal. The model does include 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 resistant homozygotes would only move at the base larval
movement rate. Of course, an individual leaving a transgenic plant only has a 5% change of moving
onto a non-Bt plant. Various estimates 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 relatively unimportant compared to other parameters and a
precise estimate would be a misapplication of resources.
Days of larval movement - Movement and larval development are assumed to occur over this period of
days. Movement occurs once per day. Daily survivorship 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 dispersal) , and of the remaining 64.5%, 36.5% will
disperse (see above). Of these dispersants, only 5% will move onto non-Bt plants, increasing the
refuge size from 5% to ca. 6%. If the base larval movement is low (remember it is spread over days as
well), then this process can lead to an increase in the refuge population and actually decrease the rate of
resistance evolution. I put this into the model just to see what effect it might have. The value was
varied between 5, 10 and 15 and it had very little effect on the results.
30
------- |