I e-mailed Dr. Matthew Forister about getting a copy of the R code he used for his paper that we are discussing today entitled: “Ant association facilitates the evolution of diet breadth in a lycaenid butterfly“. About 1 hour later I got the following response:

*Hi Thomas,*

*Sure, I’m happy to provide the code! See attached. It’s been a while since I looked at this, but I just verified that it does work, and I added a bit of annotation at the top that should help. Please let me know if you or your group has any questions.*

*thanks for your interest,*

*Matt*

I specifically asked if I could post the code on our blog so our members could run the code and test the simulations. I have posted the code below. You can find out more about Dr. Forister’s research on his faculty webpage, or his lab webpage.

`####### Code for running simulations in Forister et al. 2011, "Ant association facilitates..."`

## Below there are two functions: runAntsPlants and plotAntsPlants

## After these functions have been loaded, you can run the simulation with defaults thus:

## result <- runAntsPlants()

## Of course parameters can be changed by adding specifications in those parentheses.

## And then time-plots of the results can be generated with:

## plotAntsPlants(result)

##

## Note that there is one external package that needs to be loaded (Hmisc) which provides

## an error bar function for the graphing.

library(Hmisc)

runAntsPlants<-function(nRep=1000,nGen=100,elSurvMin=0.15,elSurvMax=0.65,lpSurvMin=0.1,lpSurvMax=0.55,

paSurvMin=0.35,paSurvMax=0.8,aeMin=10,aeMax=180,AFmin=1,AFmax=1,MFmin=1,MFmax=1,

Tend=1,freqMed=0.5,eggK=50000,startPropK=0.75,astragFlowerFactor=1.95,medQualFactor=0.053,

antFactor=2.32,flowerFactor=7.53,medFecFactorWF=0.67,medFecFactorNF=0.17){

cat(“ants and plants model running \n”)

## set up data objects

eggsMed<-array(dim=c(nRep,nGen))

eggsAst<-array(dim=c(nRep,nGen))

larvaMed<-array(dim=c(nRep,nGen))

larvaAst<-array(dim=c(nRep,nGen))

pupaMed<-array(dim=c(nRep,nGen))

pupaAst<-array(dim=c(nRep,nGen))

adultMed<-array(dim=c(nRep,nGen))

adultAst<-array(dim=c(nRep,nGen))

## row1=median, row2=ub, row3=lb

summaryEM<-array(dim=c(3,nGen))

summaryEA<-array(dim=c(3,nGen))

summaryLM<-array(dim=c(3,nGen))

summaryLA<-array(dim=c(3,nGen))

summaryPM<-array(dim=c(3,nGen))

summaryPA<-array(dim=c(3,nGen))

summaryAM<-array(dim=c(3,nGen))

summaryAA<-array(dim=c(3,nGen))

## vectors for lambda

curLambda<-numeric(nGen-1)

lambda<-numeric(nRep)

## begin simulation and cycle through reps and generations

for (i in 1:nRep){

cat(“rep”,i,”\n”)

## draw AF and MF for initialization

MF<-MFmin + rbeta(1,2,3) * (MFmax – MFmin)

AF<-AFmin + rbeta(1,2,3) * (AFmax – AFmin)

## figure out proportion of eggs on Medicago

regOut<-lm(c(1/astragFlowerFactor,1) ~ c(0,1))

prefMed<-regOut$coefficients[2] * AF + regOut$coefficients[1]

propMed<-freqMed * prefMed

## startPropK gives the ratio of the initial egg population to carrying capacity for eggs

eggsMed[i,1]<-eggK * startPropK * propMed

eggsAst[i,1]<-eggK * startPropK *(1 – propMed)

## begin looping through generations

for (j in 1:nGen){

## draw AF and MF for current generation

MF<-MFmin + rbeta(1,2,3) * (MFmax – MFmin)

AF<-AFmin + rbeta(1,2,3) * (AFmax – AFmin)

## figure out proportion of eggs that will be laid this gen on each host plant

prefMed<-regOut$coefficients[2] * AF + regOut$coefficients[1]

propMed<-freqMed * prefMed

## figure out relatvie larval survival on medicago (based on ants and flower, once per gen)

medSurvival<-medQualFactor * (antFactor * Tend + 1 – Tend) * (flowerFactor * MF + 1 – MF)

## if there is at least one egg go ahead

if (eggsMed[i,j]>=1 | eggsAst[i,j]>=1){

## survival of eggs to newly hatched larvae, host, flowers and ants don’t matter here

if (eggsMed[i,j]>=1) larvaMed[i,j]<-eggsMed[i,j] * (elSurvMin + rbeta(1,2,3) * (elSurvMax – elSurvMin))

else larvaMed[i,j:nGen]<-rep(0,nGen-j+1)

if (eggsAst[i,j]>=1) larvaAst[i,j]<-eggsAst[i,j] * (elSurvMin + rbeta(1,2,3) * (elSurvMax – elSurvMin))

else larvaAst[i,j:nGen]<-rep(0,nGen-j+1)

## survival of newly hatched larvae to pupation, depends on host, and if med, ants and flowers

if (larvaMed[i,j]>=1) pupaMed[i,j]<-medSurvival *

(larvaMed[i,j] * (lpSurvMin + rbeta(1,2,3) * (lpSurvMax – lpSurvMin)))

else pupaMed[i,j:nGen]<-rep(0,nGen-j+1)

if (larvaAst[i,j]>=1) pupaAst[i,j]<-larvaAst[i,j] * (lpSurvMin + rbeta(1,2,3) * (lpSurvMax – lpSurvMin))

else pupaAst[i,j:nGen]<-rep(0,nGen-j+1)

## survival of pupae to adults, host, flowers and andts don’t matter

if (pupaMed[i,j]>=1) adultMed[i,j]<-pupaMed[i,j] * (paSurvMin + rbeta(1,2,3) * (paSurvMax – paSurvMin))

else adultMed[i,j:nGen]<-rep(0,nGen-j+1)

if (pupaAst[i,j]>=1) adultAst[i,j]<-pupaAst[i,j] * (paSurvMin + rbeta(1,2,3) * (paSurvMax – paSurvMin))

else adultAst[i,j:nGen]<-rep(0,nGen-j+1)

## reproduction, requires at least two adults, and doesn’t occur in the last generation

if ((adultMed[i,j] + adultAst[i,j]) >= 2){

if (j < nGen){

astragEggs<-aeMin + rbeta(((adultAst[i,j])/2),2,3) * (aeMax – aeMin)

medWFeggs<-aeMin + rbeta(((adultMed[i,j])/2 * MF),2,3) * ((aeMax * medFecFactorWF) – aeMin)

medNFeggs<-aeMin + rbeta(((adultMed[i,j])/2 * (1 – MF)),2,3) * ((aeMax * medFecFactorNF) – aeMin)

totalEggs<-sum(c(astragEggs,medWFeggs,medNFeggs))

if (totalEggs > eggK) totalEggs<-eggK

eggsMed[i,j+1]<-totalEggs * propMed

eggsAst[i,j+1]<-totalEggs * (1 – propMed)

}

}

else {

eggsMed[i,j:nGen]<-rep(0,nGen-j+1)

eggsAst[i,j:nGen]<-rep(0,nGen-j+1)

larvaMed[i,j:nGen]<-rep(0,nGen-j+1)

larvaAst[i,j:nGen]<-rep(0,nGen-j+1)

pupaMed[i,j:nGen]<-rep(0,nGen-j+1)

pupaAst[i,j:nGen]<-rep(0,nGen-j+1)

adultMed[i,j:nGen]<-rep(0,nGen-j+1)

adultAst[i,j:nGen]<-rep(0,nGen-j+1)

}

}

else {

eggsMed[i,j:nGen]<-rep(0,nGen-j+1)

eggsAst[i,j:nGen]<-rep(0,nGen-j+1)

larvaMed[i,j:nGen]<-rep(0,nGen-j+1)

larvaAst[i,j:nGen]<-rep(0,nGen-j+1)

pupaMed[i,j:nGen]<-rep(0,nGen-j+1)

pupaAst[i,j:nGen]<-rep(0,nGen-j+1)

adultMed[i,j:nGen]<-rep(0,nGen-j+1)

adultAst[i,j:nGen]<-rep(0,nGen-j+1)

}

if (j != 1){

if ((adultMed[i,j-1] + adultAst[i,j-1]) > 0) curLambda[j-1]<-(adultMed[i,j] + adultAst[i,j])/(adultMed[i,j-1] + adultAst[i,j-1])

else curLambda[i]<-NA

}

## end of generation loop

}

lambda[i]<-exp(mean(log(curLambda),na.rm=TRUE))

## end of rep loop

}

## calculate median and quantiles for each generation, I can make this mean and sd if prefered

lb<-0.25

ub<-0.75

for (j in 1:nGen){

summaryEM[,j]<-quantile(eggsMed[,j],prob=c(lb,0.5,ub),na.rm=TRUE)

summaryEA[,j]<-quantile(eggsAst[,j],prob=c(lb,0.5,ub),na.rm=TRUE)

summaryLM[,j]<-quantile(larvaMed[,j],prob=c(lb,0.5,ub),na.rm=TRUE)

summaryLA[,j]<-quantile(larvaAst[,j],prob=c(lb,0.5,ub),na.rm=TRUE)

summaryPM[,j]<-quantile(pupaMed[,j],prob=c(lb,0.5,ub),na.rm=TRUE)

summaryPA[,j]<-quantile(pupaAst[,j],prob=c(lb,0.5,ub),na.rm=TRUE)

summaryAM[,j]<-quantile(adultMed[,j],prob=c(lb,0.5,ub),na.rm=TRUE)

summaryAA[,j]<-quantile(adultAst[,j],prob=c(lb,0.5,ub),na.rm=TRUE)

}

## determine proportion of replicates ending in extinction

propExtinct<-sum((adultMed[,nGen] + adultAst[,nGen]) == 0)/nRep

cat(“Proportion of replicates extinct: “,propExtinct,”\n”)

cat(“lamda: median=”, quantile(lambda,probs=0.5),” lb=”, quantile(lambda,probs=0.025),” ub=”,quantile(lambda,probs=0.975),”\n”)

## make list to be returned

resultsSummary<-list(summaryEM,summaryEA,summaryLM,summaryLA,summaryPM,summaryPA,summaryAM,summaryAA)

names(resultsSummary)<-c(“eggsMed”,”eggsAst”,”larvaMed”,”larvaAst”,”pupaMed”,”pupaAst”,”adultMed”,

“adultAst”)

out<-list(resultsSummary,eggsMed,eggsAst,larvaMed,larvaAst,pupaMed,pupaAst,adultMed,adultAst,

propExtinct,lambda)

names(out)<-c(“summary”,”eggsMed”,”eggsAst”,”larvaMed”,”larvaAst”,”pupaMed”,”pupaAst”,”adultMed”,

“adultAst”,”propExtinct”,”lambda”)

return(out)

## end of function

}

## Hmisc needed for errorbar plots, library location should be adjusted

library(Hmisc,lib.loc=”~/Library/R/2.8/library”)

## this plots the summary data, which is probably the most useful for now

plotAntsPlants<-function(results=NULL,pdf=FALSE,myFile=”./antsPlants.pdf”,myfile2=”./antsPlantsLam.pdf”){

## make errbar plots

if (pdf == TRUE) pdf(file=myFile,width=16,height=11)

else quartz()

par(mfrow=c(2,4))

par(pty=”s”)

errbar(1:dim(results$eggsMed)[2],results$summary$eggsMed[2,],results$summary$eggsMed[3,],results$summary$eggsMed[1,],

xlab=”Generation number”,ylab=”No. eggs on Medicago”)

axis(1)

axis(2)

box()

errbar(1:dim(results$eggsMed)[2],results$summary$eggsAst[2,],results$summary$eggsAst[3,],results$summary$eggsAst[1,],

xlab=”Generation number”,ylab=”No. eggs on Astragalus”)

axis(1)

axis(2)

box()

errbar(1:dim(results$eggsMed)[2],results$summary$larvaMed[2,],results$summary$larvaMed[3,],results$summary$larvaMed[1,],

xlab=”Generation number”,ylab=”No. larva on Medicago”)

axis(1)

axis(2)

box()

errbar(1:dim(results$eggsMed)[2],results$summary$larvaAst[2,],results$summary$larvaAst[3,],results$summary$larvaAst[1,],

xlab=”Generation number”,ylab=”No. larva on Astragalus”)

axis(1)

axis(2)

box()

errbar(1:dim(results$eggsMed)[2],results$summary$pupaMed[2,],results$summary$pupaMed[3,],results$summary$pupaMed[1,],

xlab=”Generation number”,ylab=”No. pupa on Medicago”)

axis(1)

axis(2)

box()

errbar(1:dim(results$eggsMed)[2],results$summary$pupaAst[2,],results$summary$pupaAst[3,],results$summary$pupaAst[1,],

xlab=”Generation number”,ylab=”No. pupa on Astragalus”)

axis(1)

axis(2)

box()

errbar(1:dim(results$eggsMed)[2],results$summary$adultMed[2,],results$summary$adultMed[3,],results$summary$adultMed[1,],

xlab=”Generation number”,ylab=”No. adults from Medicago”)

axis(1)

axis(2)

box()

errbar(1:dim(results$eggsMed)[2],results$summary$adultAst[2,],results$summary$adultAst[3,],results$summary$adultAst[1,],

xlab=”Generation number”,ylab=”No. adults from Astragalus”)

axis(1)

axis(2)

box()

#if (pdf == TRUE) dev.off()

## lambda plot

#if (pdf == TRUE) pdf(file=myFile2,width=16,height=11)

#else quartz()

#boxplot(results$lambda,ylab=”lambda”)

#if(pdf == TRUE) dev.off()