page 1

///General notes

// I used some java or c shorthand

// A += C is the same as A = A + C

// A *= C is the same as A = A * C

//The model uses 3 loci (unused loci fitnesses are set to 1.0), and I 

tend to think of it as a cube of fitnesses and

//frequencies

//when you see nested for loops like:

 for (L1 = SS to RR)

 for (L2 = SS to RR)

 for (L3 = SS to RR)

...

//It means the code is looping through the whole cube perfoming some 

operation

/**

*

* @author Mike Caprio

*/

ModelGuts 

 public int ModelRun(LocusmodelView lmv) 

 stop = false

 generation = 0

 Input fitnessL1[SS]

Input dominanceL1

 fitnessL1[RR] = 1.0

 fitnessL1[RS] = (fitnessL1[RR] - fitnessL1[SS]) * dominanceL1 + 

fitnessL1[SS]

 Input fitnessL2[SS]

Input dominanceL2

 fitnessL2[RR] = 1.0

 fitnessL2[RS] = (fitnessL2[RR] - fitnessL2[SS]) * dominanceL2 + 

fitnessL2[SS]

 Input fitnessL3[SS]

Input dominanceL3

 fitnessL3[RR] = 1.0

 fitnessL3[RS] = (fitnessL3[RR] - fitnessL3[SS]) * dominanceL3 + 

fitnessL3[SS]

Input refugeProp

Input refugeMortality

Input dailyRefugeSurvival

 Input initialResistanceAlleleFrequency

Input baseLarvalMovement

Input daysOfLarvalMovement

Input maxGenerations

TITLE: Model Guts Pseudo

AUTHOR: Mike Caprio, Mississippi State University

DATE: November 1, 2010
page 2

Input includeLarvalMovement //Difference between structured refuge 

or seed-mix

Include asymmetricalMovement //Turns on asymmetrical movement 

- boolean

Input includeL2 //Changes between single or dual locus 

model

Input includeL3 //Includes the third locus

 Input proportionNonCompliant

Input noncompliantMixing

Input nonRandomMating

if not(includeL2)

 for j=SS to RR

 fitnessL2[j] = 1.0

if not(includeL3)

 for j=SS to RR

 fitnessL3[j] = 1.0

//Now convert rates to daily rates - we do this by taking the dayth 

root of not moving

//If 40% of larvae move during the selection phase, we take the 

proportion that don't move

//(movement isn't multiplicative, not moving is) and take the 

daysOfSelection root of that number

//to get the proportion that don't move per day.

//allMovement includes base larval movement and asymmetrical 

movement off of Bt plants

 allMovement = 1.0 - (1.0 - asymmetricalMovement) * (1.0 -

baseLarvalMovement)

 proportionalMovement = 1.0 - exp(log((1.0 -

allMovement)/dayOfSelection))

 baseLarvalMovement = 1.0 - exp(log(1.0 - baseLarvalMovement) / 

daysOfSelection)

//Overall survivorship is the oberseved mortality of wild type 

insects

//The expect surivorship assuming multiplicativity is fitnessL1[SS] 

* fitnesssL2[SS] *fitnessL3[SS]

 nonAdditiveInteraction = OverallSurvivorship/(fitnessL1[SS] * 

fitnesssL2[SS] *fitnessL3[SS]);

//Because there are three two locus interactions, we take the third 

root of the total non-additive interaction

 //and use that to determine the survivorship of the individuals 

susceptible to two loci and susceptible to the third

 nonAdditiveInteraction2L = exp(log(nonAdditiveInteraction)/3);
page 3

//This is the wild type individual...it's fitness is the observed 

survivorship on the three gene pyramid

//Survivorships is a 3x3x3 cube with the fitness for each 3 locus 

genotype

 Survivorships[SS][SS][SS] = OverallSurvivorship

 Survivorships[RR][RR][RR] = 1.0

//Next are the three genotypes where we have to include the 

interactions among susceptibility

//to two loci while the third locus is RR

 Survivorships[RR][SS][SS] = fitnessL2[SS] * fitnessL3[SS] * 

nonAdditiveInteraction2L

//Sanity checks. Fitness can't go up if you add resistance alleles

 if (Survivorships[RR][SS][SS] > min(fitnessL2[SS],fitnessL3[SS]))

 Survivorships[RR][SS][SS] = min(fitnessL2[SS],fitnessL3[SS])

 if(Survivorships[RR][SS][SS] < OverallSurvivorship) 

 Survivorships[RR][SS][SS] = OverallSurvivorship

 Survivorships[SS][RR][SS] = fitnessL1[SS]*fitnessL3[SS] * 

nonAdditiveInteraction2L

 if (Survivorships[SS][RR][SS] > min(fitnessL1[SS],fitnessL3[SS]))

 Survivorships[SS][RR][SS] = min(fitnessL1[SS],fitnessL3[SS])

 if(Survivorships[SS][RR][SS] < OverallSurvivorship) 

 Survivorships[SS][RR][SS] = OverallSurvivorship

 Survivorships[SS][SS][RR] = fitnessL1[SS] * fitnessL2[SS] * 

nonAdditiveInteraction2L

 if (Survivorships[SS][SS][RR] > min(fitnessL1[SS],fitnessL2[SS]))

 Survivorships[SS][SS][RR] = min(fitnessL1[SS],fitessL2[SS])

 if(Survivorships[SS][SS][RR] < OverallSurvivorship) 

 Survivorships[SS][SS][RR] = OverallSurvivorship

 

 //Now cases where just one locus is susceptible...this collapses to 

the single toxin data provided

//Note that we already have assigned Survivorships[RR][RR][RR]

 Survivorships[SS][RR][RR] = fitnessL1[SS]

 Survivorships[RS][RR][RR] = (fitnessL1[RR] - fitnessL1[SS]) * 

dominancel1 + fitnessL1[SS]

 Survivorships[RR][SS][RR] = fitnessL2[SS]

 Survivorships[RR][RS][RR] = (fitnessL2[RR] - fitnessL2[SS]) * 

dominancel2 + fitnessL2[SS]

 Survivorships[RR][RR][SS] = fitnessL3[SS]
page 4

 Survivorships[RR][RR][RS] = (fitnessL3[RR] - fitnessL3[SS]) * 

dominancel3 + fitnessL3[SS]

 Survivorships[SS][RS][RR] = (Survivorships[SS][RR][RR] -

Survivorships[SS][SS][RR]) * dominancel2 + Survivorships[SS][SS][RR]

 Survivorships[RS][SS][RR]= (Survivorships[RR][SS][RR] -

Survivorships[SS][SS][RR]) * dominancel1 + Survivorships[SS][SS][RR]

 Survivorships[RR][RS][SS] = (Survivorships[RR][RR][SS] -

Survivorships[RR][SS][SS]) * dominancel2 + Survivorships[RR][SS][SS]

 Survivorships[RS][RR][SS] = (Survivorships[RR][RR][SS] -

Survivorships[SS][RR][SS]) * dominancel1 + Survivorships[SS][RR][SS]

 Survivorships[RR][RR][RS] = (Survivorships[RR][RR][RR] -

Survivorships[RR][RR][SS]) * hl3 + Survivorships[RR][RR][SS]

 Survivorships[SS][RR][RS] = (Survivorships[SS][RR][RR] -

Survivorships[SS][RR][SS]) * hl3 + Survivorships[SS][RR][SS]

//Genotypes with just one resistance allele. I use 

OverallSurvivorship, which is actually Survivorships[SS][SS][SS]

 Survivorships[RS][SS][SS] = (Survivorships[RR][SS][SS] -

OverallSurvivorship) * dominancel1 + OverallSurvivorship

 Survivorships[SS][RS][SS] = (Survivorships[SS][RR][SS] -

OverallSurvivorship) * dominancel2 + OverallSurvivorship

 Survivorships[SS][SS][RS] = (Survivorships[SS][SS][RR] -

OverallSurvivorship) * dominancel3 + OverallSurvivorship

//Double heterozygotes...these could all be calculated two 

different ways

 Survivorships[RR][RS][RS] = (Survivorships[RR][RS][RR] -

Survivorships[RR][RS][SS]) * hl3 + Survivorships[RR][RS][SS]

 Survivorships[RS][SS][RS] = (Survivorships[RS][SS][RR] -

Survivorships[RS][SS][SS]) * hl3 + Survivorships[RS][SS][SS]

 Survivorships[RS][RR][RS] = (Survivorships[RS][RR][RR] -

Survivorships[RS][RR][SS]) * hl3 + Survivorships[RS][RR][SS]

 Survivorships[SS][RS][RS] = (Survivorships[SS][RS][RR] -

Survivorships[SS][RS][SS]) * hl3 + Survivorships[SS][RS][SS]

 Survivorships[RS][RS][SS] = (Survivorships[RR][RS][SS] -

Survivorships[SS][RS][SS])*(hl1+hl2)/2 + Survivorships[SS][RS][SS]

 Survivorships[RS][RS][RR] = (Survivorships[RR][RS][RR] -

Survivorships[SS][RS][RR])*(hl1+hl2)/2 + Survivorships[SS][RS][RR]

//The triple heterozygotes..could be calculated 3 different ways

 Survivorships[RS][RS][RS] = (Survivorships[RR][RS][RS] -

Survivorships[SS][RS][RS])*(hl1+hl2)/2 + Survivorships[SS][RS][RS]

//Convert Daily survivorships from Survivorships using dayth root 

again

 for (L1 = SS to RR)

 for (L2 = SS to RR)

 for (L3 = SS to RR)

 DailySurv[L1][L2][L3] = 

exp(log(Survivorships[L1][L2][L3])/daysOfSelection)
page 5

 //Initialize frequency arrays to Hardy-Weinberg proportions

//p = (1.0 - initialResistanceAlleleFrequency

//q = initialResistanceAlleleFrequency 

 //locus1 are the allele frequencies in compliant fields

//locus1nc are the allele frequencies in noncompliant fields

 locus1[SS] = (1.0 - initialResistanceAlleleFrequency) * (1.0 -

initialResistanceAlleleFrequency) * (1 - proportionNoncompliant)

 locus1nc[SS] = (1.0 - initialResistanceAlleleFrequency) * (1.0 -

initialResistanceAlleleFrequency) * (proportionNoncompliant)

 locus1[RS] = 2 * initialResistanceAlleleFrequency * (1.0 -

initialResistanceAlleleFrequency) * (1 - proportionNoncompliant)

 locus1nc[RS] = 2 * initialResistanceAlleleFrequency * (1.0 -

initialResistanceAlleleFrequency) * (proportionNoncompliant)

 locus1[RR] = initialResistanceAlleleFrequency * 

initialResistanceAlleleFrequency * (1 - proportionNoncompliant)

 locus1nc[RR] = initialResistanceAlleleFrequency * 

initialResistanceAlleleFrequency * (proportionNoncompliant)

//Copy these into allele frequencies for each of the other loci, 

they are identical

 for (j1 = SS to RR)

 locus2[j1] = locus1[j1]

 locus3[j1] = locus1[j1]

 locus2nc[j1] = locus1nc[j1]

 locus3nc[j1] = locus1nc[j1]

//This is the main loop of the program. We stop when we exceed 

maxGenerations or all allele frequencies under selection exceed 50%

 while (!stop){

//Bt[] is the frequencies of genotypes in the Bt portion of 

Bt fields

//Refuge[] is the frequencies of the genotypes in the refuge 

portion of the Bt fields

//Btnc[] is the frequencies of the genotypes in the 

noncompliant fields

//We need to rescale each locus in each population (compliant 

and noncompliant) to sum to a probability of 1.0

//to calculate the joint probability of each genotyp

 tc1=0; tc2 = 0; tc3 = 0; tnc1=0; tnc2 = 0; tnc3=0;

 for (j1 = SS to RR)

 tc1 += locus1[j1]

 tc2 += locus2[j1]

 tc3 += locus3[j1]

 tnc1 += locus1nc[j1]

 tnc2 += locus2nc[j1]
page 6

 tnc3 += locus3nc[j1]

 for (L1 =SS to RR)

 for (L2 = SS to RR)

 for (L3 = SS to RR) 

 Bt[L1][L2][L3] = locus1[L1]/tc1 * locus2[L2]/tc2 * 

locus3[L3]/tc2

 Refuge[L1][L2][L3] = locus1[L1]/tc1 * 

locus2[L2]/tc2 * locus3[L3]/tc3 

 Btnc[L1][L2][L3] = locus1nc[L1]/tnc1 * 

locus2nc[L2]/tnc2 * locus3nc[L3]/tnc3

//Now we rescale back to the original proportions in each 

habitat

 for (L1 = SS to RR)

 for (L2 = SS to RR)

 for (L3 = SS to RR)

 Bt[L1][L2][L3] *= 

(tc1+tc2+tc3)/(tc1+tc2+tc3+tnc1+tnc2+tnc3) * (1.0 - refugeProp)

 Refuge[L1][L2][L3] *= 

(tc1+tc2+tc3)/(tc1+tc2+tc3+tnc1+tnc2+tnc3) * (refugeProp)

 Btnc[L1][L2][L3] *= 

(tnc1+tnc2+tnc3)/(tc1+tc2+tc3+tnc1+tnc2+tnc3)

 

//Now we do selection with larval movement (if there's larval 

movement)

 if(includeLarvalMovement) 

 for (j4 = 0 to daysOfSelection) 

 for (L1 = SS to RR)

 for (L2 = SS to RR)

 for (L3 = SS to RR)

 sum = 0;

 Bt[L1][L2][L3] *= DailySurv[L1][L2][L3]

 Refuge[L1][L2][L3] *= (dailyRefugeSurvival)

 if (AsymmetricalMovement) 

 btMovement = proportionalMovement

 else

 btMovement = baseLarvalMovement

//Sum up all movement

 sum += Bt[L1][L2][L3] * btMovement

sum += Refuge[L1][L2][L3] * baseLarvalMovement

//Remove migrants from their original 

habitat

 Bt[L1][L2][L3] -= (Bt[L1][L2][L3] * btMovement)

 Refuge[L1][L2][L3] -= Refuge[L1][L2][L3] * 

baseLarvalMovement
page 7

//Return them to new habitats based on 

the refuge proportion

 Refuge[L1][L2][L3] += sum * refugeProportion

 Bt[L1][L2][L3] += sum * (1.0 - refugeProportion)

 //This section is used for the structured refuge

 else 

 for (L1 = SS to RR){

 for (L2 = SS to RR){

 for (L3 = SS to RR){

 //Because it all happens at once, we don't use the 

daily rates

 Bt[L1][L2][L3] *= Survivorships[L1][L2][L3]

Refuge[L1][L2][L3] *= (1.0 - refugeMortality)

 Btnc[L1][L2][L3] *= Survivorships[L1][L2][L3]

 //Now we use the final compartment, the 

non-random mating compartment, that contains

 //individuals from Bt fields that do not 

mate at random with the refuge adults

 nrMate[L1][L2][L3] = Bt[L1][L2][L3] * nr_mating

 Bt[L1][L2][L3] -= nrMate[L1][L2][L3]

 sum = 0.0 //Sum of Bt fields (include refuges), with the 

exception of the nonrandom mating compartment

 sumnr = 0.0 //Sum of population in nonrandom mating portion of Bt 

fields

 sumnc = 0.0 //Sum of population in noncompliant fields

 for (L1 = SS to RR) 

 for (L2 = SS to RR) 

 for (L3 = SS to RR) 

 sum += Bt[L1][L2][L3] + Refuge[L1][L2][L3]

 sumnc += Btnc[L1][L2][L3]

 sumnr += nrMate[L1][L2][L3]

//Get the inverse of the sums, unless zero, then assign 0 to 

inverse

 if (sum > 0.0)

 isum = 1.0/sum

 else

 isum = 0.0;

 if (sumnc > 0.0)

 isumnc = 1.0/sumnc

 else

 isumnc = 0.0

 if (sumnr > 0.0)

 isumnr = 1.0/sumnr

 else

 isumnr = 0.0
page 8

 

//pc is the proportion of the total population in the 

randomly mating portions of the compliant fields

 pc = sum / (sum + sumnc + sumMix + sumnr)

 //pnc is the proportion of the total population in the 

noncomliant fields

 pnc = sumnc / (sum + sumnc + sumMix + sumnr)

//pnr is the proportion of the total population in the 

nonrandomly mating portions of the compliant fields

 pnr = sumnr / (sum + sumnc + sumMix + sumnr)

 //Now we convert each group sum to 1 to do the mating

//Hardy-Weinberg requires the sum of the genotypes to be 1.0, 

and we rescale

//each compartment to 1.0

//Note: I could have add the Refuge and Bt compartments, as 

they are now adults and

//form a panmictic population (but not nrMate, the 

nonrandomly mating compartment)

 for (L1 = SS to RR

 for (L2 = SS to RR)

 for (L3 = SS to RR)

 Bt[L1][L2][L3] *= isum

Refuge[L1][L2][L3] *= isum

Btnc[L1][L2][L3] *= isumnc

nrMate[L1][L2][L3] *= isumnr

 //This gets the sum in each of the four compartments

 for (L1 = SS to RR)

 for (L2 = SS to RR)

 for (L3 = SS to RR)

 tsumbt += Bt[L1][L2][L3]

tsumref += Refuge[L1][L2][L3]

tsumnc += Btnc[L1][L2][L3]

tsumnr += nrMate[L1][L2][L3]

 //Zero out the frequency arrays.

 for (j1 = SS to RR)

 locus1[j1] = 0.0

 locus2[j1] = 0.0

 locus3[j1] = 0.0

 locus1nc[j1] = 0.0

 locus2nc[j1] = 0.0

 locus3nc[j1] = 0.0

 locus1nr[j1] = 0.0

 locus2nr[j1] = 0.0

 locus3nr[j1] = 0.0
page 9

 //Now we sum the frequencies of of the genotypes at each locus in 

each compartment

 for (L1 = SS to RR)

 for (L2 = SS to RR)

 for (L3 = SS to RR)

 locus1[L1] += Bt[L1][L2][L3]

locus1[L1] += Refuge[L1][L2][L3]

locus1nc[L1] += Btnc[L1][L2][L3]

locus2[L2] += Bt[L1][L2][L3]

locus2[L2] += Refuge[L1][L2][L3]

locus2nc[L2] += Btnc[L1][L2][L3]

 locus3[L3] += Bt[L1][L2][L3]

 locus3[L3] += Refuge[L1][L2][L3]

locus3nc[L3] += Btnc[L1][L2][L3];

locus1nr[L1] += nrMate[L1][L2][L3]

locus2nr[L2] += nrMate[L1][L2][L3]

 locus3nr[L3] += nrMate[L1][L2][L3]

 //Now we calculate allele frequencies. We have dealt with 

genotypes until now

//l1A are the alleles frequencies at locus 1 in compliant 

fields

//l1Anc are the alleles frequencies at locus 1 in 

noncompliant fields

//l1Anr are the alleles frequencies at locus 1 in compliant, 

but nonrandomly mating fields

//Because these are allele freqeuncies, there are only two 

array indices, S and R

 l1A[S] = locus1[SS] + 0.5 * locus1[RS]

 l1A[R] = 1.0 - l1A[S]

 l2A[S] = locus2[SS] + 0.5 * locus2[RS]

 l2A[R] = 1.0 - l2A[S]

 l3A[S] = locus3[SS] + 0.5 * locus3[RS]

 l3A[R] = 1.0 - l3A[S]

 l1Anc[S] = locus1nc[SS] + 0.5 * locus1nc[RS]

 l1Anc[R] = 1.0 - l1Anc[S]

 l2Anc[S] = locus2nc[SS] + 0.5 * locus2nc[RS]

 l2Anc[R] = 1.0 - l2Anc[S]

 l3Anc[S] = locus3nc[SS] + 0.5 * locus3nc[RS]

 l3Anc[R] = 1.0 - l3Anc[S]

 l1Anr[S] = locus1nr[SS] + 0.5 * locus1nr[RS]

 l1Anr[R] = 1.0 - l1Anr[S]

 l2Anr[S] = locus2nr[SS] + 0.5 * locus2nr[RS]

 l2Anr[R] = 1.0 - l2Anr[S]

 l3Anr[S] = locus3nr[SS] + 0.5 * locus3nr[RS]

 l3Anr[R] = 1.0 - l3Anr[S]
page 10

//Now we complete mating in each compartment using HardyWeinberg and rescale to the original population sizes

//the first line is p^2 * pc

//the fourth line is 2pq * pc

//the seventh line is q^2 * pc

 locus1[SS] = l1A[S] * l1A[S] * pc

 locus1nc[SS] = l1Anc[S] * l1Anc[S] * pnc

 locus1nr[SS] = l1Anr[S] * l1Anr[S] * pnr

 locus1[RS] = 2 * l1A[S] * l1A[R] * pc

 locus1nc[RS] = 2 * l1Anc[S] * l1Anc[R] * pnc

 locus1nr[RS] = 2 * l1Anr[S] * l1Anr[R] * pnr

 locus1[RR] = l1A[R] * l1A[R] * pc;

 locus1nc[RR] = l1Anc[R] * l1Anc[R] * pnc

 locus1nr[RR] = l1Anr[R] * l1Anr[R] * pnr

 locus2[SS] = l2A[S] * l2A[S] * pc

 locus2nc[SS] = l2Anc[S] * l2Anc[S] * pnc

 locus2nr[SS] = l2Anr[S] * l2Anr[S] * pnr

 locus2[RS] = 2 * l2A[S] * l2A[R] * pc

 locus2nc[RS] = 2 * l2Anc[S] * l2Anc[R] * pnc

 locus2nr[RS] = 2 * l2Anr[S] * l2Anr[R] * pnr

 locus2[RR] = l2A[R] * l2A[R] * pc

 locus2nc[RR] = l2Anc[R] * l2Anc[R] * pnc

 locus2nr[RR] = l2Anr[R] * l2Anr[R] * pnr

 locus3[SS] = l3A[S] * l3A[S] * pc

 locus3nc[SS] = l3Anc[S] * l3Anc[S] * pnc

 locus3nr[SS] = l3Anr[S] * l3Anr[S] * pnr

 locus3[RS] = 2 * l3A[S] * l3A[R] * pc

 locus3nc[RS] = 2 * l3Anc[S] * l3Anc[R] * pnc

 locus3nr[RS] = 2 * l3Anr[S] * l3Anr[R] * pnr

 locus3[RR] = l3A[R] * l3A[R] * pc

 locus3nc[RR] = l3Anc[R] * l3Anc[R] * pnc

 locus3nr[RR] = l3Anr[R] * l3Anr[R] * pnr

 //Now we need to combine the nrMate and refuge+Bt alleles for 

random oviposition

 for (j1 = SS to RR)

 locus1[j1] += locus1nr[j1]

 locus2[j1] += locus2nr[j1]

 locus3[j1] += locus3nr[j1]

 //Now move some of the eggs back and forth

 //This simulates mated adults moving between compliant and 

noncompliant fields

 for (j1 = SS to RR)

 temp = locus1[j1] + locus1nc[j1] * complianceMix * (1 -

propNC) -

 locus1[j1] * complianceMix * propNC

 locus1nc[j1] += locus1[j1] * complianceMix * propNC -

 locus1nc[j1] * complianceMix * (1 - propNC)
page 11

 locus1[j1] = temp

 temp = locus2[j1] + locus2nc[j1] * complianceMix * (1 -

propNC) -

 locus2[j1] * complianceMix * propNC

 locus2nc[j1] += locus2[j1] * complianceMix * propNC -

 locus2nc[j1] * complianceMix * (1 - propNC)

 locus2[j1] = temp

 temp = locus3[j1] + locus3nc[j1] * complianceMix * (1 -

propNC) -

 locus3[j1] * complianceMix * propNC

 locus3nc[j1] += locus3[j1] * complianceMix * propNC -

 locus3nc[j1] * complianceMix * (1 - propNC)

 locus3[j1] = temp

//Default is to stop. If one or more allele frequencies at 

selected loci are less than 0.5

//we set stop to false and continue the simulation

 stop = true

 if (l1A[1] < 0.5 OR l1Anc[1] < 0.5)

 stop = false

 else if (includeL2 AND (l2A[1] < 0.5 OR l2Anc[1] < 0.5))

 stop = false

 else if (includeL3 AND (l3A[1] < 0.5 OR l3Anc[1] < 0.5))

 stop = false

 if (generation >= MaxGenerations)

 stop = true

 

 generation += 1

//This concludes the main while loop of the program.

//Now just output the final results

//These are the final resistance allele frequencies at the three 

loci

Output L1A[R]

if (includeLocus2)

Output L2A[R]

if (includeLocus3)

Output L3A[R]

Output generations