## R Scripts reproducing the simulations described in the Supplemtary Note.
## These scripts reproduce the simulations and results 
## i) on the evolution of genetic diversity during outbreaks (Figures S8-S9)
## ii) on the fixation of mutations during successive outbreaks (Figure S10). 
## They require the devel version of adegenet (1.3-0 or later), 
## which can be installed from http://adegenet.r-forge.r-project.org/,
## as well as R version 2.9.2 or later. Created objects are saved in the working directory. 


#### GENETIC DIVERSITY DURING AN OUTBREAK (Figures S8-S9) ####

## INSTALL AND LOAD PACKAGE ##
install.packages("adegenet", repos="http://R-Forge.R-project.org")
library(adegenet) 


## TRACKING PAIRWISE DISTANCES ## 
## BUILD A REFERENCE POPULATION 
ref.pd <- haploPopDiv(1e4, max.pop.size=100, max.nb.pop=1, birth.func=1.1, ini.pop=1e2, regen=TRUE, p.new.pop=0, 
mu=1e-6, death=0, track="div") 
save(ref.pd, file="ref.pd.100.RData")


## OUTBREAK SIMULATIONS 
NREP <- 50 

## replicates R=1.1 
x.1.1 <- replicate(NREP, haploPopDiv(1000, max.pop.size=1e5, max.nb.pop=1, birth.func=1.1, ini.pop=100, 
regen=FALSE, p.new.pop=0, mu=1e-6, death=0, track="div"), simplify=FALSE) 
save(x.1.1, file="x.1.1.RData") 

## replicates R=1.5 
x.1.5 <- replicate(NREP, haploPopDiv(1000, max.pop.size=1e5, max.nb.pop=1, birth.func=1.5, ini.pop=100, 
regen=FALSE, p.new.pop=0, mu=1e-6, death=0, track="div"), simplify=FALSE) 
save(x.1.5, file="x.1.5.RData") 

## replicates R=2 
x.2 <- replicate(NREP, haploPopDiv(1000, max.pop.size=1e5, max.nb.pop=1, birth.func=2, ini.pop=100, regen=FALSE, 
p.new.pop=0, mu=1e-6, death=0, track="div"), simplify=FALSE) 
save(x.2, file="x.2.RData") 


## TRACKING ALLELE FREQUENCIES ## 
## REFERENCE POPULATION 
ref.f <- haploPopDiv(1e4, max.pop.size=1e2, max.nb.pop=1, birth.func=1.1, ini.pop=1e2, regen=TRUE, p.new.pop=0, 
mu=1e-6, death=0, track="freq") 
save(ref.f, file="ref.f.100.RData") 


## OUTBREAKS 
f.1.1 <- replicate(NREP, haploPopDiv(1000, max.pop.size=1e5, max.nb.pop=1, birth.func=1.1, ini.pop=100, 
regen=FALSE, p.new.pop=0, mu=1e-6, death=0, track="freq"), simplify=FALSE) 
save(f.1.1, file="f.1.1.RData") 

## replicates R=1.5 
f.1.5 <- replicate(NREP, haploPopDiv(1000, max.pop.size=1e5, max.nb.pop=1, birth.func=1.5, ini.pop=100, 
regen=FALSE, p.new.pop=0, mu=1e-6, death=0, track="freq"), simplify=FALSE) 
save(f.1.5, file="f.1.5.RData") 

## replicates R=2 
f.2 <- replicate(NREP, haploPopDiv(1000, max.pop.size=1e5, max.nb.pop=1, birth.func=2, ini.pop=100, regen=FALSE, 
p.new.pop=0, mu=1e-6, death=0, track="freq"), simplify=FALSE) 
save(f.2, file="f.2.RData") 


## GRAPHICS ## 
## AUXILIARY FUNCTIONS 
getOneDiff <- function(x, lag=10){ 
   peak <- which.max(x$popSize) 
   gen2 <- x$div[[peak]] 
   gen1 <- x$div[[peak-lag+1]] 
   res <- mean(gen2) - mean(gen1) 
   return(res) 
} 

getRefDiff <- function(x, lag=10){ 
   L <- length(x$div) 
   if(L< 1000+lag) stop("simulation too short") 
   posi <- sample(1000:(L-lag), 1) 
   gen1 <- x$div[[posi]] 
   gen2 <- x$div[[posi+lag]] 
   res <- mean(gen2) - mean(gen1) 
   return(res) 
} 

transp <- function(col, alpha=.5){ 
   res <- apply(col2rgb(col),2, function(c) rgb(c[1]/255, c[2]/255, c[3]/255, alpha)) 
   return(res) 
} 

res <- lapply(list(x.1.1, x.1.5, x.2), function(e) sapply(e, getOneDiff )) 
res <- c(list(replicate(1000, getRefDiff(ref.pd))) , res) 
R <- factor(rep(c("equ","1.1","1.5","2"), sapply(res, length)), levels= c("equ","1.1","1.5","2")) 


## DRAFT OF FIGURE S9 
boxplot(unlist(res)~R, xaxt="n",yaxt="n", col=c("grey","gold1","red","deepskyblue2"), horizontal=TRUE, outline=FALSE) 


## DRAFT OF FIGURE S10 
par(mfrow=c(2,2)) 

## ref 
par(mar=c(4.1,4.1,.6,.6)) 
temp.ref <- unlist(lapply(ref.f$div[1001:2000], function(e) e/100)) 
hist(temp.ref, col="grey", xlab="Allele frequencies", ylab="Relative frequency", proba=TRUE,main="", ylim=c(0,15)) 

## f.1.1 
par(mar=c(4.1,4.1,.6,.6)) 
temp.1.1 <- unlist(lapply(f.1.1, function(sim) unlist(mapply(function(f,s) f/s, sim$div, sim$popSize)) )) 
hist(temp.1.1, col="gold1", xlab="Allele frequencies", ylab="Relative frequency", proba=TRUE, main="") 
hist(temp.ref, col=transp("black",.2), xlab="Allele frequencies", ylab="Relative frequency", proba=TRUE, add=TRUE, 
border=NA) 

## f.1.5 
par(mar=c(4.1,4.1,.6,.6)) 
temp.1.5 <- unlist(lapply(f.1.5, function(sim) unlist(mapply(function(f,s) f/s, sim$div, sim$popSize)) )) 
hist(temp.1.5, col="red", xlab="Allele frequencies", ylab="Relative frequency", proba=TRUE, main="") 
hist(temp.ref, col=transp("black",.2), xlab="Allele frequencies", ylab="Relative frequency", proba=TRUE, add=TRUE, 
border=NA) 

## f.2 
par(mar=c(4.1,4.1,.6,.6)) 
temp.2 <- unlist(lapply(f.2, function(sim) unlist(mapply(function(f,s) f/s, sim$div, sim$popSize)) )) 
hist(temp.2, col="deepskyblue2", xlab="Allele frequencies", ylab="Relative frequency", proba=TRUE, main="") 
hist(temp.ref, col=transp("black",.2), xlab="Allele frequencies", ylab="Relative frequency", proba=TRUE, add=TRUE, 
border=NA) 


## STUDENT TESTS FOR THE MEAN GENETIC DISTANCE 
lapply(res, t.test, mu=0, alter='greater') 



#### FIXATION OF MUTATIONS IN SUCCESSIVE OUTBREAKS (Figure S10) ####

## DEFINE VARIABLES ##
R = 1.1
K <- 5e5
MU <- 5*(2.3e-8/365)
L <- 4e6
N.REPLI <- 10
N.PROPAG <- 30


## FUNCTION PERFORMING SIMULATIONS ##
makeSuccOut <- function(n.outbreak, n.propag=N.PROPAG, ...){
 n.gen <- which.max(haploPopDiv(1e6, regen=FALSE,p.new.pop=0, ini.pop=n.propag, mu=0, haplo.l=10, track="nb", birth=R, death=0, max.pop=K, quiet=TRUE)$popSize)-1

 res <- numeric()

 temp <- haploPop(n.gen, regen=FALSE, p.new.pop=0, ini.pop=n.propag, mu=MU, haplo.l=L, birth=R, death=0, max.pop=K, quiet=TRUE)
 propag <- sample.haploPop(temp, n.propag)
 res <- c(res, sum(table(unlist(propag$pop[[1]])) == n.propag))

    for(i in 2:n.outbreak){
        cat("\n== i : ",i,"==")
        temp <- haploPop(n.gen, regen=FALSE, p.new.pop=0, ini.obj=propag, mu=MU, haplo.l=L, birth=R, death=0, max.pop=K, quiet=TRUE)
        propag <- sample.haploPop(temp, n.propag)
        res <- c(res, sum(table(unlist(propag$pop[[1]])) == n.propag))
    }

    return(res)
}


## PERFORM SIMULATIONS ##
nbMutFix.100rep <- lapply(1:100, function(i) makeSuccOut(n.out=100)) # takes about 2 days
save(nbMutFix.100rep, file="nbMutFix.100rep.RData")


## BUILD DATA FRAME OF RESULTS ##
X <- data.frame(nbMutFix.100rep)-1
colnames(X) <- paste("rep",1:100, sep=".")


## COMPUTE QUANTILES / MEAN
lmin <- apply(X, 1, quantile, 0.01)
lmax <- apply(X, 1, quantile, 0.99)
lmean <- apply(X, 1, mean)


## COMPUTE ESPERANCE OF NB MUT FIX
## equals MU*L for each generation (77 generations per outbreak)
theo <- cumsum(rep(MU*L*77, 100))


## MAKE FIGURE (DRAFT OF FIGURE S10) ##
matplot(X, lwd=2, xlab="Years", ylab="Number of fixed mutations", col="transparent", xaxt="n", cex.lab=1.3)
polygon(c(1:100, 100:1), c(lmin, rev(lmax)), col="black", border="black")
matplot(X,type="l", lwd=1, add=TRUE, col=rainbow(102)[1:100], lty=1:3)
lines(1:100, lmin, lwd=4, col="black")
lines(1:100, lmax, lwd=4, col="black")
lines(1:100, lmean, lwd=5, col="red")
lines(1:100, theo, lwd=5, col="blue", lty=2)

xax.at <- seq(0, 110, by= 10*(365/385)) # one tick mark every 10 years
axis(side=1, at=xax.at, label=seq(0, le=length(xax.at), by=10))