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 initialResistanceAlleleFreq uency

Input baseLarvalMovement

Input daysOfLarvalMovement

Input maxGenerations 

### TITLE : Model Guts Pseudo

AUTHOR : Mike Caprio, Mississippi State University

DATE : November 1, 2010page 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

fit nessL2[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

all Movement = 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] *fi tnessL3[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

nonA dditiveInteraction2L = 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

Survivors hips[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]) *

d ominancel2 + 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]

//Geno types 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][S S] = (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 + Survivor ships[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 us ing 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 no ncompliant 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) * (pr oportionNoncompliant)

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 ge notypes 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 *

locu s3[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 la rval

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 thei r original

habitat

Bt[L1][L2][L3]

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

Refuge[L1][L2][L3]

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

baseLarvalMovementpage 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[L 1][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 m ating 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.0page 8 //pc is the proportion of the to tal 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[L 1][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.0page 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]

l 1A[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

- l1 Anc[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 Hardy

Weinberg 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;

locu s1nc[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

l ocus3nr[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 m ated 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] += locus 1[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 pr ogram.

//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