Population Genetics

The model has four di-allelic, autosomal, resistance genes, one for each type of insecticidal traits, including Cry1A, Cry1F, Cry2Ab, and Vip3A, which are designated based on their mode of action and/ or site of action. For the purposes of this model, traits Cry1Ab, Cry1Ac, and Cry1A.105 are all considered to have the same mode of action and thus share the same resistance gene. We assumed that there is no linkage between the four resistance genes and that no mutations occur after the start of simulations.

Fig. 1.

Hectares planted with corn in 2013 in nine states.  http://usda.mannlib.cornell.edu/usda/current/Acre/Acre-06-30-2014.pdf

Landscape

We have attempted to represent the complex, heterogeneous crop landscape in the southern United States consisting of corn, cotton, soybean, peanut, and sorghum (Figs. 1 and 2). The relative proportions of crops in Fig. 2 and information about non-Bt corn refuge (described later in the text) allowed us to categorize the states. In our analysis, we focused on GA, Texas (TX), and MS. Mississippi seemed to be a conservative choice relative to IRM, because it was the state with the lowest predicted durability for Bt cotton, according to Gustafson et al. (2006).

Corn earworm has up to six generations in the southern United States, but we only modeled the generations under Bt selection and did not model the generations in early spring or late fall when corn earworm feeds on wild plants, weeds, and volunteer corn (Babu et al. 2014). During the first two generations of corn earworm under Bt selection, the landscape in the model mainly consists of corn. During Generations 3 and 4, the landscape consists of cotton and unstructured refuge (i.e., non-Bt cotton, soybean, peanut, and sorghum).

Insecticidal cotton has been planted in the southern United States for many years, first with a single trait and now as a pyramid of two traits. No refuge is required for Bt cotton because of the availability of refuge in other crops (Onstad and Knolhoff 2014). According to a 2007 survey of insect control management practices of farmers performed by the U.S. Department of Agriculture (USDA), 96–97% of cotton in GA and MS and 48% of cotton in TX is Bt cotton (USDA 2014). However, more recent estimates from 2012 indicate even higher proportions of Bt cotton in the three states. Based on the report by Williams (2013), we use 0.01, 0.02, and 0.16 as the proportions of non-Bt cotton in GA, MS, and TX, respectively. We assume that a pyramid of two to three insecticidal traits, including Vip3A, will be planted in the future: half with Cry1Ac and Cry2Ab2 and half expressing Cry1Ac and Cry1F. In a second scenario, in which corn has three traits, we included Vip3A in both cotton crops.

Corn Refuge

Corn refuge consists of structured and unstructured refuge. Structured refuge is required by companies and the U.S. Environmental Protection Agency to be planted along with Bt corn. Unstructured refuge consists of non-Bt corn that is planted for any other reason, including the choice to plant more than that required. (In our model, blended refuge produced by mixing seed in the commercial bags is not unstructured refuge.) The USDA and several states in the southern United States, including MS, do not collect data on area planted with non-Bt or Bt corn (USDA-NASS 2013). Thus, our estimates for the total corn refuge in the southern United States were obtained from two sources.

Fig. 2.

Hectares planted with cotton, soybean, peanut, and sorghum in 2013 in nine states.  http://usda.mannlib.cornell.edu/usda/current/Acre/Acre-06-302014.pdf

First, we used observations by Babu (2013) working under the direction of Dr. Fred Musser at Mississippi State University. Of 48 cornfields tested in MS in 2012, 85% (41 fields) were positive for at least one Bt protein: 42% (20 fields) were positive for Cry1A and Cry2A, 42% (20 fields) were positive for Cry1F, and 2% (1 field) was positive for Cry1A alone. Thus, 15% of blocks visited were non-Bt corn.

According to a 2010 USDA survey of insect control management practices of farmers, 13% (SE = 4%) of corn acreage in GA was planted with Bt corn, indicating that 87% of corn is refuge (USDA 2014). For TX and North Carolina, the 2010 survey estimated 35 and 55% corn refuge, respectively (USDA 2014). Another standard USDA survey found 31% non-Bt corn in TX in 2012 (USDA-NASS 2013).

In a 2013 survey of corn growers performed at the request of the Agricultural Biotechnology Stewardship Committee, an industry group, at least 46% of corn growers in cotton-growing regions reported that they complied with the refuge size requirement for corn borers in the southern United States (Smith and Smith 2014). The survey also learned that 25% partially met the size requirement and 29% did not plant a refuge. We assumed 50% compliance with block refuge size requirements in all simulations.

Total refuge was calculated by combining both structured and unstructured corn refuge. First, we assumed that 50% of growers comply with the requirement for a 20% block refuge (the other 50% do not plant a required refuge). Thus, 10% of proportional area of Bt corn should be structured refuge (0.1 × 0.85, 0.1 × 0.13 and 0.1 × 0.65 for MS, GA, and TX, respectively). Unstructured refuge is the total proportional area of non-Bt corn minus the calculated structured refuge. Therefore, unstructured corn refuge is 0.15 − 0.085 = 0.065, 0.87 − 0.013 = 0.857, and 0.35 − 0.065 = 0.285, respectively, for MS, GA, and TX.

Calculation of Refuge Effectiveness

Effective refuge is the percent of refuge (non-Bt crops and non-crop hosts) that effectively contributes to susceptible insect production. The effective refuge value adjusts land area numbers within a given insect generation to account for adoption of Bt crop, insecticide use in the same crop without Bt traits, and amount and quality of alternative hosts.

We found information about insecticide use in cotton, soybean, peanut, and sorghum in various sources. In non-Bt cotton, 64% of the fields are treated and 64.6% of the corn earworm larvae are killed by insecticide (Gustafson et al. 2006). In soybean, 38% of the fields are treated with insecticides lethal to corn earworm ( http://quickstats.nass.usda.gov/results/114C2D53-D65E-3B50-A14F-17E 81DA4DA99; data for MS accessed May 9, 2014). In treated soybean fields, mortality rate is 86% of larvae 6–7 d after treatment (Akin et al. 2012; Herbert et al. 2012a, b). For sorghum, 10% of the fields are treated, resulting in 82% insecticide-related mortality 4–6 d after treatment ( http://quickstats.nass.usda.gov/results/ A6E9EBC1-DD9A-3A70-BD9E-D70D701DDAEF; data for TX accessed July 22, 2014; Adams et al. 2013a, b). In GA, 85% of the peanut fields are treated with insecticides lethal to corn earworm and mortality rate is 90% (M. Abney, University of Georgia, personal communication).

The calculation of the proportion of effective refuge for a specific crop is presented in the following equations. First, the proportions of land area planted to Bt, P_B, and non-Bt, P_N, are calculated for each crop:

where AD is Bt crop adoption level (proportional area) and Ref is the product of the proportion of required refuge and compliance rate.

Gustafson et al. (2006) point out that, when corn earworm feeds on cotton, soybean, sorghum, and peanut in late summer, soybean produces only a tenth as many adults as would a similar area of untreated non-Bt cotton. We assumed that peanut has the same production of corn earworm as soybean. Based on knowledge of corn earworm on sorghum in TX, we assume that sorghum produces one-fifth the number of moths per ha as cotton when neither is treated with insecticide. When we calculate the effective refuge contribution of soybean, peanut, and sorghum for the model, we multiply the landscape proportion by the proportional productivity coefficient. Based on the total acreage for each crop in a state (Fig. 2), we calculated the proportional land area, LP, for each crop. We then calculated the total cotton refuge for Generations 3 and 4 by applying the following equation:

where i for non-corn crops. A(i) is the productivity relative to cotton (1 for cotton, 0.2 for sorghum, and 0.1 for soybean and peanut). T(i) is the proportion of non-Bt crop treated with insecticide lethal to corn earworm, and LS(i) is the proportion of larvae surviving this treatment. And P_B(i) is 0 for soybean, sorghum, and peanut, as there is no Bt variety of these crops in southern United States. If we do not consider alternative hosts, then the proportion of effective refuge, Ref_eff, is calculated by dividing the productive non-Bt part by the productive total crop, which is the special case of equation 2.

For GA and MS, equation 2 calculated effective refuge levels for Bt cotton of 0.032 and 0.354, respectively. For MS, soybean is important for effective refuge. For GA and MS, P_N for corn is equal to the effective refuge for the first two generations in benchmark simulations, because we assume that no insecticide lethal to corn earworm is used exclusively on non-Bt corn.

We model the system in TX differently from those in GA and MS by making the first two generations of corn earworm in TX infest both corn and sorghum. Thus, we use half the acreage planted to sorghum (Fig. 2) as refuge for Bt corn and half for refuge for Bt cotton (in final two generations).

We assume that sorghum produces equal numbers of corn earworm moths per ha as non-Bt corn that has not been sprayed. The young sorghum crop is not treated with insecticides in either of the first two generations. Non-Bt corn is not sprayed in the first generation, but 20% of the crop is treated in the second generation. For blended corn refuge scenarios, the unstructured refuge is sprayed and none of the blended refuge is sprayed. For block refuge scenarios in TX, the standard simulation assumed 20% of all non-Bt corn is sprayed. (We also modeled the assumption that farmers do not spray any cornfields, which is the same as we modeled for MS and GA.) We assume 30% survival of corn earworm on treated non-Bt corn. Therefore, every ha of non-Bt corn produces 0.8 + 0.2 × 0.3 = 0.86 moth for every moth produced by sorghum in the second generation in TX.

The modeling of the third and fourth generations in TX is similar to that for the other two states. The productivity (oviposition preference and immature survival) of sorghum is one-fifth that of cotton, as noted previously. After accounting for insecticide use, the weighted average production of moths is 1 from sorghum for every 3 from non-Bt cotton. For TX, we calculated an effective refuge level for Bt cotton of 0.144 (equation 2).

Insect and Crop Synchrony

For simplicity, we assume that the generations of the corn earworm are discrete in time and that each generation infests a well-defined set of crop species that are available for the entire period of activity for immature stages. Within each of the three most common host plants, corn earworm prefers to oviposit on 1) tasselling and silk stage corn, 2) blooming and early pod set soybean, and 3) flowering cotton, which is preferred over preflowering cotton (squares). The plant species are not equally attractive, as our discussion of effective refuge suggested. In 2002–2003, Jackson et al. (2008) observed no overlap between cotton and corn in five southern states. From 2002 to 2006 in the southern United States, Head et al. (2010) observed that circa 2% of moths flying during “corn” generations had fed as larvae on cotton. The results of Head et al. (2010) also indicate that the natal fields of corn earworm are not soybean until after the generations on corn. Hutchison and Storer (2010) modeled the first two generations on corn only and the third and fourth generations on cotton and soybean. Thus, our benchmark assumption is that corn earworm infests corn for two generations, followed by two generations on cotton, soybean, peanut, and sorghum.

In the sensitivity analysis for MS, we allow corn earworm to exist on both corn and cotton during one generation. When corn is infested for three generations, we make the relative productivity of corn in the third generation equal to one-fifth the production of corn earworm on cotton (Gustafson et al. 2006). When cotton is infested by the second to fourth generations, we make the relative productivity of corn earworm in its second generation on cotton equal to one-fifth the production on corn. For this scenario, soybean, peanut, and sorghum are still only available as refuge for the last two generations. Thus, in this scenario, effective cotton refuge in MS is only 0.012 in the second corn earworm generation (first on cotton).

When we simulate the overlap of one generation of corn earworm on both corn and cotton, corn should be added into equation 2 as a fifth crop. For the third generation, mainly feeding on cotton, A(i) is 0.2, 0.2, 1, and 0.1 for corn, sorghum, cotton, and other crops, respectively.

Calculation of Survival in a Given Patch in Landscape

We used the model to evaluate block refuge and blended refuge (mixture of Bt and non-Bt seed). In both cases, we account for survival before and after a single movement event, probability of larval movement, and survival during movement. We calculate survival for both young and old larvae and assume that survival for entire larval stage is the product of the two. Young larvae are those that do not move off of their natal plant (Burkness et al. 2015). Because of a lack of data, the same rates are used for movement off of Bt and non-Bt plants resulting in symmetrical movement. Burkness et al. (2015) did not observe any early larval movement, but they found that the movement rate of old larvae is 15% with range between 0 and 33%. Thus, our benchmark value for the proportion of old larvae moving off of a corn plant was set at 0.15. For simplicity, we assumed that survival during late larval movement was 1.

The following equations calculate the survival of larvae in genotype g and in patch f for a given generation. Survival of young larvae is simply S(g,f)^q, where S(g,f)¹ is observed survival for a larva in genotype g that remains on a single plant in patch f for entire larval stage. The exponent q was chosen to cause most of mortality in young larvae (q > 0.5). The proportion of older larvae, OL, of genotype g on refuge plant surviving to adulthood i:

Where YL is the young larvae in genotype g in patch f, either refuge, r, or Bt plants, b. Parameter M is the probability of moving off of plant in patch f and P is the proportional area of patch. The terms in the square brackets in equation 4 represent those leaving plants not in patch f and moving to a plant in patch f, those larvae staying on natal plant, and those leaving natal plant but reaching similar plant in patch f. For refuge plants, survival S(g,r) = 1, which means survival is 1 for both young and old larvae regardless of q value. For survival on blended refuge corn plants owing to cross-pollination, we used reduced survival owing to cross-pollination for S(g,r)^(1-q) and kept S(g,b)^q = 1. The proportion of older larvae on Bt plant surviving to adulthood is:

Again exponent q distributes mortality across two larval stages. A comparison of two studies of young and old larvae by Herbert et al. (2012a, b) leads us to believe that 0.6 < q < 0.65 when larvae are exposed to insecticides. We used equations 4 and 5 for larvae infesting both corn and cotton. We used the same values of M and q in both Bt corn and Bt cotton, but the S values are different.

Survival on Corn Plants

The relative survival rates on reproductive-stage corn are based on data from field studies performed in three locations in GA and MS (Crespo et al. 2015). The data for each replicate are shown in the Suppl Appendix (online only). We calculated the mean survival by dividing the mean for all replicates for each treatment by the mean for all replicates observed in refuge corn. The means are similar to those used by Edwards et al. (2013), who based their parameters on the opinions of eight entomologists. For the three insecticidal traits acting alone, our means and standard deviations and the estimates of Edwards et al. (2013) are (0.77 ± 0.16, 0.69) for Cry1F, (0.54 ± 0.25, 0.58) for Cry1Ab, and (0.02 ± 0.01, 0.07) for Vip3A, respectively. For corn with MON810 × DAS-01507-1 × MIR162, we observed a mean survival of 0.001 with an SD of 0.003 in 2012 (Suppl Appendix [online only]). Therefore, we use the values observed in the field in 2012 to parameterize survival on reproductive-stage corn. For vegetative-stage corn, productivity is believed to be less than that on reproductive-stage corn (Gustafson et al. 2006, Edwards et al. 2013). Therefore, we use the following survival rates for homozygous susceptible larvae on vegetative corn: 0.49 on Cry1F and 0.36 on Cry1Ab (both based on expert opinions reported by Edwards et al. 2013), and 0.01 on Vip3A, which is less than the 0.02 observed in the field in reproductive-stage corn. We assume redundant killing (multiplicative survival) in all Bt corn products.

For refuge plants, we assume 100% survival unless cross-pollination reduces survival in reproductive stage corn in blended refuge. For refuge plants pollinated by adjacent Bt plants, we base our choice of survival rate for older larvae (those feeding on ears) on several sources of information. We do not model direct mortality caused by Bt pollen consumption by larvae on refuge plants, or at least we do not consider this a separate rate of survival. However, because older corn earworm larvae commonly feed on corn kernels, we do model the influence of cross-pollination from Bt plants on toxicity of refuge-plant kernels. This reduction in survival of older, homozygous susceptible larvae on refuge plants occurs in our model only on corn in the second generation when a seed blend is used for refuge.

Burkness et al. (2015) performed field experiments in Minnesota during 2012–2013 using artificial infestations to determine movement behaviors of corn earworm larvae. In these experiments, only the central plant was initially infested. By comparing the relative number of individuals surviving on the central refuge plant in pure stand refuge plots and blended plots, any negative effects of cross-pollination on survival could be inferred. For blended plots, the central refuge plant was surrounded by corn plants with MON810 × DAS-01507-1 × MIR162. Burkness et al. (2015) observed either no difference in the number of late-instar larvae surviving on central refuge plants or more larvae on the central refuge plant in blended plots than in pure stand refuge plots.

In GA and MS in 2012, D.B. and A.C. (unpublished data) performed field studies with natural infestations in a three-plant cluster planted with refuge corn and corn with MON810 ×DAS-01507-1 × MIR162. Note that in this study, the total number of larvae is the sum of the number collected and the number of exit holes in the corn ears. In MS, the refuge plant in pure stand refuge clusters had higher numbers of larvae on average compared with the blended refuge. The relative survival in blended refuge was 0.61. In GA, the relative survival was 0.95, but there was no significant difference between the number of larvae found on refuge plants in pure stand refuge versus refuge plants in blends.

In TX in 2014, P.P. (unpublished data) performed a field study with natural infestations in a 95:5 seed blend containing corn with MON810 × DAS-01507-1 × MIR162. Porter observed similar numbers of corn earworm adults surviving on seed blend refuge ears and ears in a pure stand of non-Bt corn. The relative survival rate in this study was 1 in the blended refuge.

Thus, in four out of five data sets concerning MON810 × DAS-01507-1 × MIR162 distinguished by year and/or location, there was no evidence of an effect of cross-pollination on corn earworm larvae. Therefore, our benchmark survival rate for older larvae on refuge corn in blends during the second generation is the mean of the five results described previously (1 + 1 + 1 + 0.95 + 0.61)/5 = 0.91. This rate is implemented by taking the square or cubed root of 0.91, which is 0.95 or 0.97, respectively, and assigning this as the survival rate for each toxin in model. In a sensitivity analysis, we varied survival owing to cross-pollination from 0 to 1. In all cases, the dominance of resistance for this case is the same as for all other cases, the dominance equals the survival rate for homozygous-susceptible.

It is possible that asymmetrical movement away from Bt corn could cause movement to counteract any losses owing to cross-pollination on refuge plants, but this is a more complicated scenario to model and no data exist to test this hypothesis. In cotton, movement away from Bt and non-Bt plants may be different or the same depending on variety (Jackson et al. 2010).

Survival on Bt Cotton

Edwards et al. (2013) obtained opinions from seven entomologists to select their survival rates on Bt cotton. For Cry1Ac, Cry1Ac x Cry2Ab, Cry1Ab × Cry1F, and Cry1Ab × Vip3A, the estimated survival rates were 0.46, 0.14, 0.20, and 0.09, respectively. Sivasupramaniam et al. (2008) observed survival on Cry1Ac cotton of 0.40 and Cry1Ac × Cry2Ab2 of 0.08. Dividing the Cry2Ab pyramid estimates in both publications by the respective survival on Cry1Ac (0.14/0.46 and 0.08/0.40) allows us to find the Cry2Ab contribution to survival, which is 0.30 and 0.20 in each study. Brévault et al. (2013) measured survival on plant material in the laboratory. They observed survival of 0.70 on non-Bt cotton, 0.04 on Cry1Ac (lower than expected), and 0.01 on Cry1Ac × Cry2Ab. This means that, relative to that observed on non-Bt corn, survival owing to Cry1Ac was 0.057 and that owing to Cry2Ab was 0.25. Thus, the Cry2Ab survival seems to be 0.25 based on all three studies. However, the survival rate of 0.057 on Cry1Ac observed by Brévault et al. (2013) seems similar to that observed by Anilkumar et al. (2009) but much lower than that observed by others (Kurtz 2005, Sivasupramaniam et al. 2008), including the eight experts cited by Edwards et al. (2013). We chose to use the survival rate observed by Sivasupramaniam et al. (2008), which is 0.40 on Cry1Ac cotton. Dividing survival rate of 0.2 on Cry1Ab × Cry1F cotton proposed by experts cited in Edwards et al. (2013) by the 0.4 value just chosen for Cry1A, we obtain a survival rate of 0.5 for Cry1F cotton.

Bommireddy and Leonard (2008) observed 7–28% survival on Vip3A cotton plant tissues after 4 d of exposure. The authors noted that survival would likely be less over a longer period. Therefore, we believe that a benchmark value for survival of the entire larval stage on Vip3A cotton should be 0.07. For the other three survival rates, assuming multiplicative survival on pyramids, we chose the following survival rates to be benchmark values: 0.40 on Cry1Ac, 0.25 on Cry2Ab, and 0.50 on Cry1F.

Brévault et al. (2013) concluded that inheritance of resistance on Cry1Ac cotton was not completely recessive. The value of h was 0.25. Brévault et al. (2013) found no fitness costs for resistant larvae on non-Bt cotton. They concluded that resistant phenotypes of corn earworm have incomplete resistance on at least one insecticidal trait. However, to be conservative in our analysis, we chose to assume complete resistance (100% survival) for homozygous-resistant larvae on all Bt cotton. We also assume redundant killing for larvae on pyramided Bt cotton.

Calculation of Genotypic Survival

We assume there is complete resistance (100% survival) when resistance is homozygous at a single locus. Heterozygote survival is calculated in the typical manner. Given Wss and Wrr survival for the two homozygotes, Wrs = hWrr + (1-h)Wss, where h is dominance of resistance. When h = 0, resistance is completely recessive. As demonstrated by Tabashnik et al. (2004), we modeled functional dominance based on dose of insecticide. We assume that dominance, h, equals the proportional survival rate for homozygous-susceptible larvae over the entire larval period. Although Brévault et al. (2013) calculated a value of h = 0.25 for expression of the allele conferring resistance to Cry1Ac, we do not use this value because it would actually be lower than all of the values derived from our rule for Cry1Ab and Cry1Ac. We use h to calculate survival of heterozygotes and then use the exponent q to calculate survival for each larval stage (equations 4 and 5).

Our model allows cross-resistance to be modeled as complete or partial and symmetrical or asymmetrical. The parameter CR_1-2 is defined as the cross-resistance for the Bt protein 2 when resistance occurs to Bt protein 1. When CR_1-2 = 0, there is no cross-resistance. When CR_1-2 = × (0 < × < 1), there is partial cross-resistance between Bt proteins 1 and 2. In other words, CR_1-2 is a variable that denotes cross-resistance and shows that when 100% survival occurs on Bt protein 1, there will be x survival on Bt protein 2. Although there is no evidence of cross-resistance in corn earworm with the common insecticidal traits, we chose a conservative value of 0.1 for cross-resistance for all resistance genes.

The formulas used for calculating survival for two-toxin pyramids are shown in Table 1. Formulas similar to these were created in our computer program (R language in Rx64 V3.02) for a pyramid of three or more Bt proteins. For example, given wild-type susceptible alleles A, X, and M and resistance alleles B, Y, and N, we can calculate the survival rate of the genotype ABXYMN on a pyramid of three Bt proteins by separately calculating the survival rates of genotype AB on Bt protein 1, XY on Bt protein 2, and MN on Bt protein and multiplying these three survival rates to get the overall survival rate. For example, to calculate the survival rate S_sr2 of het-erozygote AB on Bt protein 1, the program performs the following steps. First, it calculates the AB survival rate without cross-resistance, S_sr = (1-h₁) × S₁ + h₁ × s1. Then, it calculates the survival rate based on cross-resistance from the Y allele, which is CR_2-1 × h₂ × s2. The survival rate of AB with cross-resistance of Y is S_sr1 = S_sr+(1-S_sr) × (CR_2-1 × h₂ × s₂). Third, it calculates the survival rate based on cross-resistance from the N allele, which is CR_3-1 × h₃ × s_3. The survival rate of AB with cross-resistance caused by Y and N is S_sr2 = S_sr1+(1- S_sr1) × (CR_3-1 × h₃ × s₃), where s₃ and h₃ are survival rate of NN on Bt protein 3 and associated dominance value, and CR_3-1 is partial cross-resistance between Bt proteins 1 and 3, respectively. Note that the value of CR_3-1 is the survival on Bt protein 1 when 100% survival occurs on Bt protein 3. This algorithm can be applied to other genotypes. We apply the cross-resistance factor to the separate survival rates for young and old larvae (before and after movement).

Table 1.

Example of formulas for survival rates for larvae feeding on a pyramid of Bt proteins, with subscripts representing resistance genes

Required Corn Refuge

We evaluated IRM strategies that require 5–20% refuge either as a block or blend. We assume that block refuges are planted with a 50% compliance rate.

Calculating Changes in Resistance Allele Frequency

The final step is calculating the landscape-level change in allele frequency by determining mean survival of each genotype. Survival in each part of landscape is multiplied by the proportion of the landscape (after all adjustments described previously) that consists of that plant type. After the proportion of each genotype in the landscape is calculated, the allele frequencies for the adults are determined. The genotypes of the offspring are derived from the adult gene pool by assuming Hardy–Weinberg ratios. Edwards et al. (2013) assumed that initial resistance allele frequencies for Cry1A and Cry2 were 0.002 and 0.0002, respectively. We take a more conservative approach and assume that the benchmark values of initial resistance allele frequencies for all loci are 0.005.

Computer Program Verification

We verified the logic of the computer program (R language in R ×64 V3.02) in several ways. The program was compared with two sets of previous model results. A simple, abstract, two-gene model was used to test two-gene scenarios by eliminating larval movement factors. The results from Mallet and Porter (1992) were also verified by studying one major resistance gene and only incorporating late larval movement. In addition, the genotype survival rates were manually calculated under three-gene scenarios to match the results derived from model. We used a general two-gene model developed in the R programming language to test the effective-refuge and larval-movement algorithms.