#GENERAL FUNCTIONS---------------------------------------
#--------------------------------------------------------
superpose.eb <- 
function (x, y, ebl, ebu = ebl, length = 0.08, ...) 
    arrows(x, y + ebu, x, y - ebl, angle = 90, code = 3, 
    length = length, ...)
#---------------------------------------GENERAL FUNCTIONS


#FIGURE 1------------------------------------------------
#--------------------------------------------------------

#Read in raw data----
RPS.AF.df<-read.table("RPS_Anc_Fit.txt", header=TRUE)
attach(RPS.AF.df)

#Create fitness vectors----
wRS<-RPS.AF.df[RPS.AF.df$Competition=="w_R_to_S",]$Fitness
wPR<-RPS.AF.df[RPS.AF.df$Competition=="w_P_to_R",]$Fitness
wSP<-c(0,0,0,0,0,0)

#Calculate means----
WRS<-mean(RPS.AF.df[RPS.AF.df$Competition=="w_R_to_S",]$Fitness)
WPR<-mean(RPS.AF.df[RPS.AF.df$Competition=="w_P_to_R",]$Fitness)
WSP<-0

#Construct mean and sem vectors----
fits<-c(WRS,WPR,WSP)
sem<-c(sd(wRS)/sqrt(length(wRS)),sd(wPR)/sqrt(length(wPR)),sd(wSP)/sqrt(length(wSP)))

#Plot bar graph----
Fit = matrix(fits,1,3)  
eblb = matrix(sem,1,3) 

par(mar=c(7.1,5.7,4.1,2.1))

textsize<-1.6
numsize<-1.5

Fit.Plot <- barplot(fits, col="gray", 
    names.arg=c("","",""), 
    ylim=c(0,1), 
    ylab="relative fitness", 
    cex.lab=textsize, cex.names=textsize, cex.axis=numsize) 
superpose.eb(Fit.Plot, Fit, eblb, col="black", lwd=2, length=0.03)

mtext("Resistant\nrelative to\nSensitive", side=1, at=Fit.Plot[1], line=4.5, cex=textsize)
mtext("Producer\nrelative to\nResistant", side=1, at=Fit.Plot[2], line=4.5, cex=textsize)
mtext("Sensitive\nrelative to\nProducer", side=1, at=Fit.Plot[3], line=4.5, cex=textsize)

lines(c(0,4),c(1,1),lty="dotdash",col="black",lwd=2)

lift<-0.05
text(Fit.Plot[1],fits[1]+sem[1]+lift,"*",cex=3)
text(Fit.Plot[2],fits[2]+sem[2]+lift,"*",cex=3)

#One sample t test for S vs. R competitions----
t.test(wRS, alternative="two.sided",mu=1)

#One sample t test for R vs. P competitions----
t.test(wPR,alternative="two.sided",mu=1)

#------------------------------------------------FIGURE 1

setwd("C:/Users/BEN/My Documents/R Files")

#FIGURE 2------------------------------------------------
#--------------------------------------------------------

#Read in data----
rps.eco.df<-read.table("RPS_ECO.txt", header=TRUE)

#Data processing for restricted migration rps----
P.st <- rps.eco.df[
rps.eco.df$Marker=="P" &
(rps.eco.df$Transfer==0 |
rps.eco.df$Transfer==6 |
rps.eco.df$Transfer==12 |
rps.eco.df$Transfer==18 |
rps.eco.df$Transfer==24 |
rps.eco.df$Transfer==30 |
rps.eco.df$Transfer==36) &
rps.eco.df$Structure &
!rps.eco.df$R_Only &
rps.eco.df$Media == "LB_Tet",]

R.st <- rps.eco.df[
rps.eco.df$Marker=="R" &
(rps.eco.df$Transfer==0 |
rps.eco.df$Transfer==6 |
rps.eco.df$Transfer==12 |
rps.eco.df$Transfer==18 |
rps.eco.df$Transfer==24 |
rps.eco.df$Transfer==30 |
rps.eco.df$Transfer==36) &
rps.eco.df$Structure &
!rps.eco.df$R_Only &
rps.eco.df$Media == "LB_Tet",]

S.st <- rps.eco.df[
rps.eco.df$Marker=="S" &
(rps.eco.df$Transfer==0 |
rps.eco.df$Transfer==6 |
rps.eco.df$Transfer==12 |
rps.eco.df$Transfer==18 |
rps.eco.df$Transfer==24 |
rps.eco.df$Transfer==30 |
rps.eco.df$Transfer==36) &
rps.eco.df$Structure &
!rps.eco.df$R_Only &
rps.eco.df$Media == "LB_Tet",]

aver.P.st <- c(mean(P.st[P.st$Transfer==0,]$Num_cells),
mean(P.st[P.st$Transfer==6,]$Num_cells),
mean(P.st[P.st$Transfer==12,]$Num_cells),
mean(P.st[P.st$Transfer==18,]$Num_cells),
mean(P.st[P.st$Transfer==24,]$Num_cells),
mean(P.st[P.st$Transfer==30,]$Num_cells),
mean(P.st[P.st$Transfer==36,]$Num_cells))

sem.P.st <- c(sd(P.st[P.st$Transfer==0,]$Num_cells)/sqrt(5),
sd(P.st[P.st$Transfer==6,]$Num_cells)/sqrt(5),
sd(P.st[P.st$Transfer==12,]$Num_cells)/sqrt(5),
sd(P.st[P.st$Transfer==18,]$Num_cells)/sqrt(5),
sd(P.st[P.st$Transfer==24,]$Num_cells)/sqrt(5),
sd(P.st[P.st$Transfer==30,]$Num_cells)/sqrt(5),
sd(P.st[P.st$Transfer==36,]$Num_cells)/sqrt(5))

aver.R.st <- c(mean(R.st[R.st$Transfer==0,]$Num_cells),
mean(R.st[R.st$Transfer==6,]$Num_cells),
mean(R.st[R.st$Transfer==12,]$Num_cells),
mean(R.st[R.st$Transfer==18,]$Num_cells),
mean(R.st[R.st$Transfer==24,]$Num_cells),
mean(R.st[R.st$Transfer==30,]$Num_cells),
mean(R.st[R.st$Transfer==36,]$Num_cells))

sem.R.st <- c(sd(R.st[R.st$Transfer==0,]$Num_cells)/sqrt(5),
sd(R.st[R.st$Transfer==6,]$Num_cells)/sqrt(5),
sd(R.st[R.st$Transfer==12,]$Num_cells)/sqrt(5),
sd(R.st[R.st$Transfer==18,]$Num_cells)/sqrt(5),
sd(R.st[R.st$Transfer==24,]$Num_cells)/sqrt(5),
sd(R.st[R.st$Transfer==30,]$Num_cells)/sqrt(5),
sd(R.st[R.st$Transfer==36,]$Num_cells)/sqrt(5))

aver.S.st <- c(mean(S.st[S.st$Transfer==0,]$Num_cells),
mean(S.st[R.st$Transfer==6,]$Num_cells),
mean(S.st[R.st$Transfer==12,]$Num_cells),
mean(S.st[R.st$Transfer==18,]$Num_cells),
mean(S.st[R.st$Transfer==24,]$Num_cells),
mean(S.st[R.st$Transfer==30,]$Num_cells),
mean(S.st[R.st$Transfer==36,]$Num_cells))

sem.S.st <- c(sd(S.st[S.st$Transfer==0,]$Num_cells)/sqrt(5),
sd(S.st[S.st$Transfer==6,]$Num_cells)/sqrt(5),
sd(S.st[S.st$Transfer==12,]$Num_cells)/sqrt(5),
sd(S.st[S.st$Transfer==18,]$Num_cells)/sqrt(5),
sd(S.st[S.st$Transfer==24,]$Num_cells)/sqrt(5),
sd(S.st[S.st$Transfer==30,]$Num_cells)/sqrt(5),
sd(S.st[S.st$Transfer==36,]$Num_cells)/sqrt(5))

#Data processing for unrestricted migration rps----
P.un <- rps.eco.df[
rps.eco.df$Marker=="P" &
(rps.eco.df$Transfer==0 |
rps.eco.df$Transfer==6 |
rps.eco.df$Transfer==12 |
rps.eco.df$Transfer==18 |
rps.eco.df$Transfer==24 |
rps.eco.df$Transfer==30 |
rps.eco.df$Transfer==36) &
!rps.eco.df$Structure &
!rps.eco.df$R_Only &
rps.eco.df$Media == "LB_Tet",]

R.un <- rps.eco.df[
rps.eco.df$Marker=="R" &
(rps.eco.df$Transfer==0 |
rps.eco.df$Transfer==6 |
rps.eco.df$Transfer==12 |
rps.eco.df$Transfer==18 |
rps.eco.df$Transfer==24 |
rps.eco.df$Transfer==30 |
rps.eco.df$Transfer==36) &
!rps.eco.df$Structure &
!rps.eco.df$R_Only &
rps.eco.df$Media == "LB_Tet",]

S.un <- rps.eco.df[
rps.eco.df$Marker=="S" &
(rps.eco.df$Transfer==0 |
rps.eco.df$Transfer==6 |
rps.eco.df$Transfer==12 |
rps.eco.df$Transfer==18 |
rps.eco.df$Transfer==24 |
rps.eco.df$Transfer==30 |
rps.eco.df$Transfer==36) &
!rps.eco.df$Structure &
!rps.eco.df$R_Only &
rps.eco.df$Media == "LB_Tet",]

aver.P.un <- c(mean(P.un[P.un$Transfer==0,]$Num_cells),
mean(P.un[P.un$Transfer==6,]$Num_cells),
mean(P.un[P.un$Transfer==12,]$Num_cells),
mean(P.un[P.un$Transfer==18,]$Num_cells),
mean(P.un[P.un$Transfer==24,]$Num_cells),
mean(P.un[P.un$Transfer==30,]$Num_cells),
mean(P.un[P.un$Transfer==36,]$Num_cells))

sem.P.un <- c(sd(P.un[P.un$Transfer==0,]$Num_cells)/sqrt(5),
sd(P.un[P.un$Transfer==6,]$Num_cells)/sqrt(5),
sd(P.un[P.un$Transfer==12,]$Num_cells)/sqrt(5),
sd(P.un[P.un$Transfer==18,]$Num_cells)/sqrt(5),
sd(P.un[P.un$Transfer==24,]$Num_cells)/sqrt(5),
sd(P.un[P.un$Transfer==30,]$Num_cells)/sqrt(5),
sd(P.un[P.un$Transfer==36,]$Num_cells)/sqrt(5))

aver.R.un <- c(mean(R.un[R.un$Transfer==0,]$Num_cells),
mean(R.un[R.un$Transfer==6,]$Num_cells),
mean(R.un[R.un$Transfer==12,]$Num_cells),
mean(R.un[R.un$Transfer==18,]$Num_cells),
mean(R.un[R.un$Transfer==24,]$Num_cells),
mean(R.un[R.un$Transfer==30,]$Num_cells),
mean(R.un[R.un$Transfer==36,]$Num_cells))

sem.R.un <- c(sd(R.un[R.un$Transfer==0,]$Num_cells)/sqrt(5),
sd(R.un[R.un$Transfer==6,]$Num_cells)/sqrt(5),
sd(R.un[R.un$Transfer==12,]$Num_cells)/sqrt(5),
sd(R.un[R.un$Transfer==18,]$Num_cells)/sqrt(5),
sd(R.un[R.un$Transfer==24,]$Num_cells)/sqrt(5),
sd(R.un[R.un$Transfer==30,]$Num_cells)/sqrt(5),
sd(R.un[R.un$Transfer==36,]$Num_cells)/sqrt(5))

aver.S.un <- c(mean(S.un[S.un$Transfer==0,]$Num_cells),
mean(S.un[R.un$Transfer==6,]$Num_cells),
mean(S.un[R.un$Transfer==12,]$Num_cells),
mean(S.un[R.un$Transfer==18,]$Num_cells),
mean(S.un[R.un$Transfer==24,]$Num_cells),
mean(S.un[R.un$Transfer==30,]$Num_cells),
mean(S.un[R.un$Transfer==36,]$Num_cells))

sem.S.un <- c(sd(S.un[S.un$Transfer==0,]$Num_cells)/sqrt(5),
sd(S.un[S.un$Transfer==6,]$Num_cells)/sqrt(5),
sd(S.un[S.un$Transfer==12,]$Num_cells)/sqrt(5),
sd(S.un[S.un$Transfer==18,]$Num_cells)/sqrt(5),
sd(S.un[S.un$Transfer==24,]$Num_cells)/sqrt(5),
sd(S.un[S.un$Transfer==30,]$Num_cells)/sqrt(5),
sd(S.un[S.un$Transfer==36,]$Num_cells)/sqrt(5))

#Data processing for restricted migration r alone----
R.alone <- rps.eco.df[
rps.eco.df$Marker=="R" &
(rps.eco.df$Transfer==0 |
rps.eco.df$Transfer==6 |
rps.eco.df$Transfer==12 |
rps.eco.df$Transfer==18 |
rps.eco.df$Transfer==24 |
rps.eco.df$Transfer==30 |
rps.eco.df$Transfer==36) &
rps.eco.df$Structure &
rps.eco.df$R_Only &
rps.eco.df$Media == "LB_Tet",]

aver.R.alone <- c(mean(R.alone[R.alone$Transfer==0,]$Num_cells),
mean(R.alone[R.alone$Transfer==6,]$Num_cells),
mean(R.alone[R.alone$Transfer==12,]$Num_cells),
mean(R.alone[R.alone$Transfer==18,]$Num_cells),
mean(R.alone[R.alone$Transfer==24,]$Num_cells),
mean(R.alone[R.alone$Transfer==30,]$Num_cells),
mean(R.alone[R.alone$Transfer==36,]$Num_cells))

sem.R.alone <- c(sd(R.alone[R.alone$Transfer==0,]$Num_cells)/sqrt(5),
sd(R.alone[R.alone$Transfer==6,]$Num_cells)/sqrt(5),
sd(R.alone[R.alone$Transfer==12,]$Num_cells)/sqrt(5),
sd(R.alone[R.alone$Transfer==18,]$Num_cells)/sqrt(5),
sd(R.alone[R.alone$Transfer==24,]$Num_cells)/sqrt(5),
sd(R.alone[R.alone$Transfer==30,]$Num_cells)/sqrt(5),
sd(R.alone[R.alone$Transfer==36,]$Num_cells)/sqrt(5))


#Extra vectors needed----
time<-c(0,6,12,18,24,30,36)
x.poly<-c(0,6,12,18,24,30,36,36,30,24,18,12,6,0)


#Color figure for ecological dynamics of unrestricted migration rps----
plot(time, aver.R.un, ylim=c(10000000,100000000000), xlab="transfer",
ylab="number of cells",log="y",col="white",axes=FALSE,cex.lab=1.5)

axis(1, at = c(0,6,12,18,24,30,36),labels=c(0,6,12,18,24,30,36),cex.axis=1.25)
axis(2, las=1, at = c(10000000,100000000,1000000000,10000000000,100000000000), 
labels=c(expression(10^7),expression(10^8),expression(10^9),expression(10^10),
expression(10^11)),cex.axis=1.25)

polygon(x.poly,c(aver.P.un+sem.P.un,rev(aver.P.un-sem.P.un)),
col="mistyrose",border=FALSE)
lines(time, aver.P.un, col="red",lwd=2,log="y")
polygon(x.poly,c(aver.R.un+sem.R.un,rev(aver.R.un-sem.R.un)),
col="lightgoldenrodyellow",border=FALSE)
lines(time, aver.R.un, col="lightgoldenrod4",lwd=2)
polygon(x.poly,c(aver.S.un+sem.S.un,rev(aver.S.un-sem.S.un)),
col="lightcyan",border=FALSE)
lines(time, aver.P.un, col="red",lty=2)
lines(time, aver.S.un, col="blue",lwd=2)
points(time, aver.P.un, col="red", pch=19,cex=1.3)
points(time, aver.R.un, col="lightgoldenrod4",pch=19,cex=1.3)
points(time, aver.S.un, col="blue",pch=19,cex=1.3)

text(c(37.05),aver.P.un[7],"P",col="red",cex=1.5)
text(c(37.05),aver.S.un[7],"S",col="blue",cex=1.5)
text(c(37.05),aver.R.un[7],"R",col="lightgoldenrod4",cex=1.5)


#Color figure for ecological dynamics of restricted migration rps----
plot(time, aver.R.st, ylim=c(10000000,100000000000), xlab="transfer",
ylab="number of cells",log="y",col="white",axes=FALSE,cex.lab=1.5)

axis(1, at = c(0,6,12,18,24,30,36),labels=c(0,6,12,18,24,30,36),cex.axis=1.25)
axis(2, las=1, at = c(10000000,100000000,1000000000,10000000000,100000000000), 
labels=c(expression(10^7),expression(10^8),expression(10^9),expression(10^10),
expression(10^11)),cex.axis=1.25)

polygon(x.poly,c(aver.P.st+sem.P.st,rev(aver.P.st-sem.P.st)),
col="mistyrose",border=FALSE)
lines(time, aver.P.st, col="red",lwd=2)
polygon(x.poly,c(aver.R.st+sem.R.st,rev(aver.R.st-sem.R.st)),
col="lightgoldenrodyellow",border=FALSE)
lines(time, aver.R.st, col="lightgoldenrod4",lwd=2)
polygon(x.poly,c(aver.S.st+sem.S.st,rev(aver.S.st-sem.S.st)),
col="lightcyan",border=FALSE)
lines(time, aver.P.st, col="red",lty=2)
lines(time, aver.S.st, col="blue",lwd=2)
points(time, aver.P.st, col="red", pch=19,cex=1.3)
points(time, aver.R.st, col="lightgoldenrod4",pch=19,cex=1.3)
points(time, aver.S.st, col="blue",pch=19,cex=1.3)

text(c(37.05),aver.P.st[7],"P",col="red",cex=1.5)
text(c(37.05),aver.S.st[7],"S",col="blue",cex=1.5)
text(c(37.05),aver.R.st[7],"R",col="lightgoldenrod4",cex=1.5)


#Color figure for ecological dynamics of restricted migration r----
plot(time, aver.R.alone, ylim=c(10000000,100000000000), xlab="transfer",
ylab="number of cells",log="y",col="white",axes=FALSE,cex.lab=1.5)

axis(1, at = c(0,6,12,18,24,30,36),labels=c(0,6,12,18,24,30,36),cex.axis=1.25)
axis(2, las=1, at = c(10000000,100000000,1000000000,10000000000,100000000000), 
labels=c(expression(10^7),expression(10^8),expression(10^9),expression(10^10),
expression(10^11)),cex.axis=1.25)



polygon(x.poly,c(aver.R.alone+sem.R.alone,rev(aver.R.alone-sem.R.alone)),
col="lightgoldenrodyellow",border=FALSE)
lines(time, aver.R.alone, col="lightgoldenrod4",lwd=2)
points(time, aver.R.alone, col="lightgoldenrod4",pch=19,cex=1.3)

#polygon(x.poly,c(aver.R.st+sem.R.st,rev(aver.R.st-sem.R.st)),
#col="lightgoldenrodyellow",border=FALSE)
#polygon(x.poly,c(aver.R.un+sem.R.un,rev(aver.R.un-sem.R.un)),
#col="lightgoldenrodyellow",border=FALSE)
#lines(time, aver.R.st, col="gray",lwd=2)
#points(time, aver.R.st, col="gray",pch=19,cex=1.3)
#lines(time, aver.R.un, col="gray",lwd=2)
#points(time, aver.R.un, col="gray",pch=19,cex=1.3)

text(c(37.05),aver.R.un[7],"R",col="lightgoldenrod4",cex=1.5)

#------------------------------------------------FIGURE 2



#FIGURE 3------------------------------------------------
#--------------------------------------------------------

#Read in data----
RPS.df<-read.table("RPS_EVO.txt", header=TRUE)

#Extract fitnesses for the three treatments----
st.all<-c(mean(RPS.df$Fitness[2:9]),mean(RPS.df$Fitness[66:73]),
mean(RPS.df$Fitness[98:105]),mean(RPS.df$Fitness[130:137]),
mean(RPS.df$Fitness[162:169]))

un.all<-c(mean(RPS.df$Fitness[10:17]),
mean(RPS.df$Fitness[74:81]),mean(RPS.df$Fitness[106:113]),
mean(RPS.df$Fitness[138:145]),
mean(RPS.df$Fitness[170:177]))

st.R<-c(mean(RPS.df$Fitness[18:25]),
mean(RPS.df$Fitness[82:89]),
mean(RPS.df$Fitness[114:121]),mean(RPS.df$Fitness[146:153]),
mean(RPS.df$Fitness[178:185]))

#Calculate means and sems----
val<-c(mean(st.all),mean(un.all),mean(st.R))-1
sem<-c(sd(st.all)/sqrt(length(st.all)),sd(un.all)/sqrt(length(un.all)),
sd(st.R)/sqrt(length(st.R)))

Fit = matrix(val,1,3) 
eblb = matrix(sem,1,3) 

#Plot the bargraph----

par(mar=c(7.1,5.7,4.1,2.1))
textsize<-1.6
numsize<-1.5

Fit.Plot <- barplot(val, col="gray", 
    names.arg=c("","",""),
    ylim=c(0,0.8), 
    ylab="relative fitness", 
    cex.lab=textsize, cex.names=textsize, cex.axis=numsize,
	xpd=FALSE, axes=FALSE) 
superpose.eb(Fit.Plot, Fit, eblb, col="black", lwd=2, length=0.03)

mtext("Restricted\nCommunity", side=1, at=Fit.Plot[1], line=3, cex=textsize)
mtext("Unrestricted\nCommunity", side=1, at=Fit.Plot[2], line=3, cex=textsize)
mtext("Restricted\nAlone", side=1, at=Fit.Plot[3], line=3, cex=textsize)


act<-seq(0,0.7,0.1)
axis(2,  at = act, labels=act+1,cex.axis=numsize)

lines(c(0.05,Fit.Plot[3]+Fit.Plot[1]),c(0,0),col="white")



x<-c(rep("a", length(st.all)),rep("b", length(un.all)),rep("c", length(st.R)))
y<-c(st.all,un.all,st.R)
fit.df<-data.frame(ind=x,fit=y)

mod <- aov(fit ~ ind, data=fit.df)
summary(mod)
TukeyHSD(mod)

lift<-(0.04)
text(Fit.Plot[1],val[1]+sem[1]+lift,"a",cex=textsize)
text(Fit.Plot[2],val[2]+sem[2]+lift,"b",cex=textsize)
text(Fit.Plot[3],val[3]+sem[3]+lift,"b",cex=textsize)

#------------------------------------------------FIGURE 3




#FIGURE 4------------------------------------------------
#---------------------------------------------------------

#Start part a---------
#Read in 100x100 simulation at time point 100----
RPS.SIM.df<-read.table("RPS_SIM_t_100.txt", header=TRUE)

summary(RPS.SIM.df)

#Extract fitnesses from RF, Perm, and Shad----
wRF<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & !RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world 
& RPS.SIM.df$number_of_transfers==100 & RPS.SIM.df$Sen>0 & RPS.SIM.df$Pro>0 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30
+RPS.SIM.df$M31+RPS.SIM.df$M32+RPS.SIM.df$M33)>0,]$mean_ending_fitness

wPerm<-RPS.SIM.df[RPS.SIM.df$is_restricted & RPS.SIM.df$permute_world & !RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world
& RPS.SIM.df$number_of_transfers==100 & RPS.SIM.df$Sen>0 & RPS.SIM.df$Pro>0 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30
+RPS.SIM.df$M31+RPS.SIM.df$M32+RPS.SIM.df$M33)>0
,]$mean_ending_fitness

wShad<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & RPS.SIM.df$r_alone & RPS.SIM.df$shadow_world
& RPS.SIM.df$number_of_transfers==100 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30+RPS.SIM.df$M31+RPS.SIM.df$M32
+RPS.SIM.df$M33)>0
,]$mean_ending_fitness

#Compute means of RF, Perm, and Shad----
WRF<-mean(wRF)
WPerm<-mean(wPerm)
WShad<-mean(wShad)

#We want to plot from 1 to the mean fitnesses on our barplot.
#Here we shift our fitnesses down by one so that when the ps file is 
#imported into Illustrator, we don't have massive "boxes"----

fits<-c(WRF,WPerm,WShad)-1
Fit = matrix(fits,1,3)

#Calculate error bars----
sem<-c(sd(wRF)/sqrt(length(wRF)),sd(wPerm)/sqrt(length(wPerm)),sd(wShad)/sqrt(length(wShad)))
eblb = matrix(sem,1,3) 

#Plot barplot----
textsize<-1.6
numsize<-1.5
par(mar=c(7.1,5.7,4.1,2.1))
Fit.Plot <- barplot(fits, col="gray", 
    names.arg=c("","",""),
    ylim=c(0,0.09), 
    ylab="relative fitness", 
    cex.lab=textsize, xpd=FALSE, axes=FALSE) 
superpose.eb(Fit.Plot, Fit, eblb, col="black", lwd=2, length=0.03)

mtext("Restricted\nCommunity", side=1, at=Fit.Plot[1], line=2.85, cex=textsize)
mtext("Restricted\nCommunity with\nPermutation", side=1, at=Fit.Plot[2], line=4.5, cex=textsize)
mtext("Restricted\nAlone with\nShadowing", side=1, at=Fit.Plot[3], line=4.5, cex=textsize)

act<-seq(0,0.09,0.02)
axis(2,  at = act, labels=act+1,cex.axis=numsize)

lines(c(0.05,Fit.Plot[3]+Fit.Plot[1]),c(0,0),col="white")


#Statistical analysis on fitnesses----
x<-c(rep("a", length(wRF)),rep("b", length(wPerm)),rep("c", length(wShad)))
y<-c(wRF,wPerm,wShad)
fit.df<-data.frame(ind=x,fit=y)

mod <- aov(fit ~ ind, data=fit.df)
summary(mod)
TukeyHSD(mod)

#Add appropriate letters on barplot----
lift<-0.003
text(Fit.Plot[1],c(fits[1]+sem[1]+lift),"a",cex=textsize)
text(Fit.Plot[2],c(fits[2]+sem[2]+lift),"b",cex=textsize)
text(Fit.Plot[3],c(fits[3]+sem[3]+lift),"c",cex=textsize)
#End part a---------


#Start part b---------

#Read in 100x100 simulation at time point 400----
RPS.SIM.df<-read.table("RPS_SIM_t_400.txt", header=TRUE)

#Extract fitnesses from RF, Perm, and Shad----
wRF<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & !RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world 
& RPS.SIM.df$number_of_transfers==400 & RPS.SIM.df$Sen>0 & RPS.SIM.df$Pro>0 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30
+RPS.SIM.df$M31+RPS.SIM.df$M32+RPS.SIM.df$M33)>0,]$mean_ending_fitness

wPerm<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$r_alone & RPS.SIM.df$permute_world & !RPS.SIM.df$shadow_world  
& RPS.SIM.df$number_of_transfers==400 & RPS.SIM.df$Sen>0 & RPS.SIM.df$Pro>0 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30
+RPS.SIM.df$M31+RPS.SIM.df$M32+RPS.SIM.df$M33)>0,]$mean_ending_fitness

wShad<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & RPS.SIM.df$r_alone & RPS.SIM.df$shadow_world
& RPS.SIM.df$number_of_transfers==400 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30+RPS.SIM.df$M31+RPS.SIM.df$M32
+RPS.SIM.df$M33)>0,]$mean_ending_fitness

#Compute means of RF, Perm, and Shad----
WRF<-mean(wRF)
WPerm<-mean(wPerm)
WShad<-mean(wShad)

#We want to plot from 1 to the mean fitnesses on our barplot.
#Here we shift our fitnesses down by one so that when the ps file is 
#imported into Illustrator, we don't have massive "boxes"----

fits<-c(WRF,WPerm,WShad)-1
Fit = matrix(fits,1,3)

#Calculate error bars----
sem<-c(sd(wRF)/sqrt(length(wRF)),sd(wPerm)/sqrt(length(wPerm)),sd(wShad)/sqrt(length(wShad)))
eblb = matrix(sem,1,3) 

#Plot barplot----
textsize<-1.6
numsize<-1.5
par(mar=c(7.1,5.7,4.1,2.1))
Fit.Plot <- barplot(fits, col="gray", 
    names.arg=c("","",""),
    ylim=c(0,0.09), 
    ylab="relative fitness", 
    cex.lab=textsize, xpd=FALSE, axes=FALSE) 
superpose.eb(Fit.Plot, Fit, eblb, col="black", lwd=2, length=0.03)

mtext("Restricted\nCommunity", side=1, at=Fit.Plot[1], line=2.85, cex=textsize)
mtext("Restricted\nCommunity with\nPermutation", side=1, at=Fit.Plot[2], line=4.5, cex=textsize)
mtext("Restricted\nAlone with\nShadowing", side=1, at=Fit.Plot[3], line=4.5, cex=textsize)

act<-seq(0,0.09,0.02)
axis(2,  at = act, labels=act+1,cex.axis=numsize)

lines(c(0.05,Fit.Plot[3]+Fit.Plot[1]),c(0,0),col="white")


#Statistical analysis on fitnesses----
x<-c(rep("a", length(wRF)),rep("b", length(wPerm)),rep("c", length(wShad)))
y<-c(wRF,wPerm,wShad)
fit.df<-data.frame(ind=x,fit=y)

mod <- aov(fit ~ ind, data=fit.df)
summary(mod)
TukeyHSD(mod)

#Add appropriate letters on barplot----
lift<-0.003
text(Fit.Plot[1],c(fits[1]+sem[1]+lift),"a",cex=textsize)
text(Fit.Plot[2],c(fits[2]+sem[2]+lift),"b",cex=textsize)
text(Fit.Plot[3],c(fits[3]+sem[3]+lift),"b",cex=textsize)
#----end part b

#------------------------------------------------FIGURE 4




#SI FIGURE 2------------------------------------------------
#-----------------------------------------------------------

#Read in data----

S1.df<-read.table("RPS_Spec_Run_1.txt",header=TRUE)
attach(S1.df)
S2.df<-read.table("RPS_Spec_Run_2.txt",header=TRUE)
attach(S2.df)


#Process sensitive, resistant, and producer growth rates----

for(row in c("B","C","D")) {
    T.high<-2525
    T.low<-225
    Low<-log(0.08)
    High<-log(0.25)
    St.col<-43
    r1.vec<-1:10
    r2.vec<-1:10
    plot(c(T.low,T.high),c(Low,High),col="white")
    for(col in 2:11) {
        well<-paste(row,col,sep="");
        t<-S1.df[(Location==well & Time<T.high & Time>T.low),1]
        y<-log(S1.df[(Location==well & Time<T.high & Time>T.low),3])
        m<-summary(lm(y~t))[[4]][[2,1]]
        b<-summary(lm(y~t))[[4]][[1,1]]
        r1.vec[[col-1]]<-m
        points(t,y,col=colors()[[St.col+col]])
        lines(t,m*t+b,col=colors()[[St.col+col]])
    }
    for(col in 2:11) {
        well<-paste(row,col,sep="");
        t<-S2.df[(Location2==well & Time2<T.high & Time2>T.low),1]
        y<-log(S2.df[(Location2==well & Time2<T.high & Time2>T.low),3])
        m<-summary(lm(y~t))[[4]][[2,1]]
        b<-summary(lm(y~t))[[4]][[1,1]]
        r2.vec[[col-1]]<-m
        points(t,y,col=colors()[[St.col+col]])
        lines(t,m*t+b,col=colors()[[St.col+col]])
    }    
    med1<-c(0,0.0625,0.08333,0.125,0.16667,0.25,0.3333,0.5,0.66667,1)
    med2<-c(0,0.05,0.07,0.1,0.11,0.14,0.2,0.43,0.6,0.75)
    r1.vec<-rev(r1.vec)
    r2.vec<-rev(r2.vec)
    med<-c(med1,med2)
    r.vec<-c(r1.vec,r2.vec)
    model<-nls(r.vec~(0.0001696133*(med))/(k+(med)),start=list(k=.5))
    #v.pred<-summary(model)[[7]][[1,1]]
    k.pred<-summary(model)[[7]][[1,1]]
    pred<-(0.0001696133*(med))/(k.pred+(med))
    
    if(row=="B") {
        s.r<-r.vec;
        #s.v.pred<-v.pred;    
        s.k.pred<-k.pred;
        s.pred<-pred;
    } else {
        if(row=="C") {
            r.r<-r.vec;
            #r.v.pred<-v.pred;    
            r.k.pred<-k.pred;
            r.pred<-pred;
        } else {
            p.r<-r.vec;
            #p.v.pred<-v.pred;    
            p.k.pred<-k.pred;
            p.pred<-pred;
        }
    }
}

par(mar=c(7.1,5.7,4.1,2.1))
textsize<-1.6
numsize<-1.5
plot(med,3600*s.r,col="blue",xlab="proportion LB",
ylab="growth rate", cex.lab=textsize, cex.axis=numsize) 

points(med,3600*r.r,col="goldenrod")
points(med,3600*p.r,col="red")
lines(sort(med),3600*sort(s.pred),col="blue")
lines(sort(med),3600*sort(r.pred),col="goldenrod")
lines(sort(med),3600*sort(p.pred),col="red")
text(.9,.42,"R",font=2,cex=1.75,col="goldenrod")
text(.9,.277,"P",font=2,cex=1.75,col="red")
text(.9,.495,"S",font=2,cex=1.75,col="blue")

#------------------------------------------------SI FIGURE 2




#SI FIGURE 3------------------------------------------------
#-----------------------------------------------------------
#Read in 100x100 simulation at time point 100----
RPS.SIM.df<-read.table("RPS_SIM_100_100_FULL_100.txt", header=TRUE)

#Extract fitnesses from RF and RA----
wRF<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & !RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world 
& RPS.SIM.df$Sen>0 & RPS.SIM.df$Pro>0 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30
+RPS.SIM.df$M31+RPS.SIM.df$M32+RPS.SIM.df$M33)>0,]$mean_ending_fitness

wRA<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world
& (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30+RPS.SIM.df$M31+RPS.SIM.df$M32
+RPS.SIM.df$M33)>0,]$mean_ending_fitness


#Compute means of RF and RA----
WRF<-mean(wRF)
WRA<-mean(wRA)


#We want to plot from 1 to the mean fitnesses on our barplot.
#Here we shift our fitnesses down by one so that when the ps file is 
#imported into Illustrator, we don't have massive "boxes"----

fits<-c(WRF,WRA)-1
Fit = matrix(fits,1,2)

#Calculate error bars----
sem<-c(sd(wRF)/sqrt(length(wRF)),sd(wRA)/sqrt(length(wRA)))
eblb = matrix(sem,1,2) 


#Plot barplot----
textsize<-1.6
numsize<-1.5
par(mar=c(7.1,5.7,4.1,2.1))
Fit.Plot <- barplot(fits, col="gray", 
    names.arg=c("",""),
    ylim=c(0,0.085), 
    ylab="relative fitness", 
    cex.lab=textsize, xpd=FALSE, axes=FALSE) 
superpose.eb(Fit.Plot, Fit, eblb, col="black", lwd=2, length=0.03)

mtext("Restricted\nCommunity", side=1, at=Fit.Plot[1], line=2.85, cex=textsize)
mtext("Restricted\nAlone", side=1, at=Fit.Plot[2], line=2.85, cex=textsize)

act<-seq(0,0.08,0.02)
axis(2,  at = act, labels=act+1,cex.axis=numsize)

lines(c(0.05,Fit.Plot[2]+Fit.Plot[1]),c(0,0),col="white")

#Statistical analysis on fitnesses----
t.test(wRA,wRF)

#Add appropriate letters on barplot----
lift<-0.003
text((Fit.Plot[1]+Fit.Plot[2])/2,max(fits[1]+sem[1]+lift,fits[2]+sem[2]+lift),"*",cex=3)

#------------------------------------------------SI FIGURE 3



#SI FIGURE 4------------------------------------------------
#-----------------------------------------------------------

#Read in 12x16 simulation at time point 36----
RPS.SIM.df<-read.table("RPS_SIM_t_100.txt", header=TRUE)

summary(RPS.SIM.df)

#Extract fitnesses from RF, UF, and RA----
wRF<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & !RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world 
& RPS.SIM.df$number_of_transfers==36 & RPS.SIM.df$Sen>0 & RPS.SIM.df$Pro>0 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30
+RPS.SIM.df$M31+RPS.SIM.df$M32+RPS.SIM.df$M33)>0,]$mean_ending_fitness

wUF<-RPS.SIM.df[!RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & !RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world
& RPS.SIM.df$number_of_transfers==36 & RPS.SIM.df$Sen>0 & RPS.SIM.df$Pro>0 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30
+RPS.SIM.df$M31+RPS.SIM.df$M32+RPS.SIM.df$M33)>0,]$mean_ending_fitness

wRA<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world
& RPS.SIM.df$number_of_transfers==36 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30+RPS.SIM.df$M31+RPS.SIM.df$M32
+RPS.SIM.df$M33)>0,]$mean_ending_fitness


#Compute means of RF, UF, and RA----
WRF<-mean(wRF)
WRA<-mean(wRA)
WUF<-mean(wUF)

#We want to plot from 1 to the mean fitnesses on our barplot.
#Here we shift our fitnesses down by one so that when the ps file is 
#imported into Illustrator, we don't have massive "boxes"----

fits<-c(WRF,WUF,WRA)-1
Fit = matrix(fits,1,3)

#Calculate error bars----
sem<-c(sd(wRF)/sqrt(length(wRF)),sd(wUF)/sqrt(length(wUF)),sd(wRA)/sqrt(length(wRA)))
eblb = matrix(sem,1,3) 


#Plot barplot----
textsize<-1.6
numsize<-1.5
par(mar=c(7.1,5.7,4.1,2.1))
Fit.Plot <- barplot(fits, col="gray", 
    names.arg=c("","",""),
    ylim=c(0,0.055), 
    ylab="relative fitness", 
    cex.lab=textsize, xpd=FALSE, axes=FALSE) 
superpose.eb(Fit.Plot, Fit, eblb, col="black", lwd=2, length=0.03)

mtext("Restricted\nCommunity", side=1, at=Fit.Plot[1], line=2.85, cex=textsize)
mtext("Unrestricted\nCommunity", side=1, at=Fit.Plot[2], line=2.85, cex=textsize)
mtext("Restricted\nAlone", side=1, at=Fit.Plot[3], line=2.85, cex=textsize)

act<-seq(0,0.055,0.01)
axis(2,  at = act, labels=act+1,cex.axis=numsize)

lines(c(0.05,Fit.Plot[3]+Fit.Plot[1]),c(0,0),col="white")


#Statistical analysis on fitnesses----
x<-c(rep("a", length(wRF)),rep("b", length(wUF)),rep("c", length(wRA)))
y<-c(wRF,wUF,wRA)
fit.df<-data.frame(ind=x,fit=y)

mod <- aov(fit ~ ind, data=fit.df)
summary(mod)
TukeyHSD(mod)

#Add appropriate letters on barplot----
lift<-0.003
text(Fit.Plot[1],c(fits[1]+sem[1]+lift),"a",cex=textsize)
text(Fit.Plot[2],c(fits[2]+sem[2]+lift),"b",cex=textsize)
text(Fit.Plot[3],c(fits[3]+sem[3]+lift),"c",cex=textsize)

#------------------------------------------------SI FIGURE 4


#SI FIGURE 5------------------------------------------------
#------------------------------------------------------------

#Read in 12x16 simulation at time point 36----
RPS.SIM.df<-read.table("RPS_SIM_t_100.txt", header=TRUE)

#Extract fitnesses from RF, Perm, and Shad----
wRF<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & !RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world 
& RPS.SIM.df$number_of_transfers==36 & RPS.SIM.df$Sen>0 & RPS.SIM.df$Pro>0 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30
+RPS.SIM.df$M31+RPS.SIM.df$M32+RPS.SIM.df$M33)>0,]$mean_ending_fitness

wPerm<-RPS.SIM.df[RPS.SIM.df$is_restricted & RPS.SIM.df$permute_world & !RPS.SIM.df$r_alone & !RPS.SIM.df$shadow_world
& RPS.SIM.df$number_of_transfers==36 & RPS.SIM.df$Sen>0 & RPS.SIM.df$Pro>0 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30
+RPS.SIM.df$M31+RPS.SIM.df$M32+RPS.SIM.df$M33)>0,]$mean_ending_fitness

wShad<-RPS.SIM.df[RPS.SIM.df$is_restricted & !RPS.SIM.df$permute_world & !RPS.SIM.df$r_alone & RPS.SIM.df$shadow_world
& RPS.SIM.df$number_of_transfers==36 & (RPS.SIM.df$Res+RPS.SIM.df$M27+RPS.SIM.df$M28+RPS.SIM.df$M29+RPS.SIM.df$M30+RPS.SIM.df$M31+RPS.SIM.df$M32
+RPS.SIM.df$M33)>0,]$mean_ending_fitness

#Compute means of RF, Perm, and Shad----
WRF<-mean(wRF)
WPerm<-mean(wPerm)
WShad<-mean(wShad)

#We want to plot from 1 to the mean fitnesses on our barplot.
#Here we shift our fitnesses down by one so that when the ps file is 
#imported into Illustrator, we don't have massive "boxes"----

fits<-c(WRF,WPerm,WShad)-1
Fit = matrix(fits,1,3)

summary(wPerm)

#Calculate error bars----
sem<-c(sd(wRF)/sqrt(length(wRF)),sd(wPerm)/sqrt(length(wPerm)),sd(wShad)/sqrt(length(wShad)))
eblb = matrix(sem,1,3) 

#Plot barplot----
textsize<-1.6
numsize<-1.5
par(mar=c(7.1,5.7,4.1,2.1))
Fit.Plot <- barplot(fits, col="gray", 
    names.arg=c("","",""),
    ylim=c(0,0.055), 
    ylab="relative fitness", 
    cex.lab=textsize, xpd=FALSE, axes=FALSE) 
superpose.eb(Fit.Plot, Fit, eblb, col="black", lwd=2, length=0.03)

mtext("Restricted\nCommunity", side=1, at=Fit.Plot[1], line=2.85, cex=textsize)
mtext("Restricted\nCommunity with\nPermutation", side=1, at=Fit.Plot[2], line=4.5, cex=textsize)
mtext("Restricted\nAlone with\nShadowing", side=1, at=Fit.Plot[3], line=4.5, cex=textsize)

act<-seq(0,0.055,0.01)
axis(2,  at = act, labels=act+1,cex.axis=numsize)

lines(c(0.05,Fit.Plot[3]+Fit.Plot[1]),c(0,0),col="white")


#Statistical analysis on fitnesses----
x<-c(rep("a", length(wRF)),rep("b", length(wPerm)),rep("c", length(wShad)))
y<-c(wRF,wPerm,wShad)
fit.df<-data.frame(ind=x,fit=y)

mod <- aov(fit ~ ind, data=fit.df)
summary(mod)
TukeyHSD(mod)

#Add appropriate letters on barplot----
lift<-0.0018
text(Fit.Plot[1],c(fits[1]+sem[1]+lift),"a",cex=textsize)
text(Fit.Plot[2],c(fits[2]+sem[2]+lift),"b",cex=textsize)
text(Fit.Plot[3],c(fits[3]+sem[3]+lift),"c",cex=textsize)

#------------------------------------------------SI FIGURE 5



#SI FIGURE 6------------------------------------------------
#[run Figure 2 first]---------------------------------------

treat<-c("S","S","S","S","S","U","U","U","U","U","A","A","A","A","A")

#Plotting and analysis on number of generations----
R.st.gen<-c(0,0,0,0,0)
R.un.gen<-c(0,0,0,0,0)
R.alone.gen<-c(0,0,0,0,0)
tps<-c(0,6,12,18,24,30)

for(tp in tps) {
R.st.gen<-R.st.gen + log2((40^6)*((R.st[R.st$Transfer==(tp+6),]$Num_cells)/(R.st[R.st$Transfer==tp,]$Num_cells)))
R.un.gen<-R.un.gen + log2((40^6)*((R.un[R.un$Transfer==(tp+6),]$Num_cells)/(R.un[R.un$Transfer==tp,]$Num_cells)))
R.alone.gen<-R.alone.gen + log2((40^6)*((R.alone[R.alone$Transfer==(tp+6),]$Num_cells)/(R.alone[R.alone$Transfer==tp,]$Num_cells)))
}

gens<-c(R.st.gen,R.un.gen,R.alone.gen)
gen.df<-data.frame(treat,gens)
model <- aov(gens ~ treat, data=gen.df)
fligner.test(gens~treat,gen.df)
summary(model)
TukeyHSD(model)
gens<-c(mean(R.st.gen),mean(R.un.gen),mean(R.alone.gen))
sem<-c(sd(R.st.gen)/sqrt(length(R.st.gen)),sd(R.un.gen)/sqrt(length(R.un.gen)),sd(R.alone.gen)/sqrt(length(R.alone.gen)))
Gen = matrix(gens,1,3)
boxplot(R.st.gen,R.un.gen,R.alone.gen,las=1,
names=c("Restricted/nFull","Unrestricted/nFull","Restricted/nAlone"),
xlab="treatment",ylab="number of generations",cex.lab=1.5,
col="gray")
text(.7,197.5,"a",cex=1.75)
text(1.9,197.5,"a",cex=1.75)
text(3.1,195.5,"b",cex=1.75)

#Plotting and analysis on number of divisions----
R.st.div<-c(0,0,0,0,0)
R.un.div<-c(0,0,0,0,0)
R.alone.div<-c(0,0,0,0,0)
tps<-c(0,6,12,18,24,30)
for(tp in tps) {
R.st.div<-R.st.div+(((7-5*(1/40))*R.st[R.st$Transfer==(tp+6),]$Num_cells) + ((5-7*(1/40))*R.st[R.st$Transfer==(tp),]$Num_cells))
R.un.div<-R.un.div+(((7-5*(1/40))*R.un[R.un$Transfer==(tp+6),]$Num_cells) + ((5-7*(1/40))*R.un[R.un$Transfer==(tp),]$Num_cells))
R.alone.div<-R.alone.div+(((7-5*(1/40))*R.alone[R.alone$Transfer==(tp+6),]$Num_cells) + ((5-7*(1/40))*R.alone[R.alone$Transfer==(tp),]$Num_cells))
}

divs<-c(R.st.div,R.un.div,R.alone.div)
div.df<-data.frame(treat,divs)
model2 <- aov(divs ~ treat, data=div.df)
fligner.test(divs~treat,div.df)
summary(model2)
TukeyHSD(model2)

divs.m<-c(mean(R.st.div),mean(R.un.div),mean(R.alone.div))
sem<-c(sd(R.st.div)/sqrt(length(R.st.div)),sd(R.un.div)/sqrt(length(R.un.div)),sd(R.alone.div)/sqrt(length(R.alone.div)))

Div = matrix(divs.m,1,3)   
eblb = matrix(sem,1,3) 

textsize<-1.6
numsize<-1.4
par(mar=c(7.1,5.7,4.1,2.1))
Div.Plot <- barplot(divs.m, col="gray", 
    names.arg=c("","",""),
    log="y", ylim=c(1,5*10^12),
    ylab="",
    cex.lab=textsize, cex.names=textsize, cex.axis=numsize, axes=FALSE) 
superpose.eb(Div.Plot, Div, eblb, col="black", lwd=2, length=0.03)

mtext("Restricted\nCommunity", side=1, at=Fit.Plot[1], line=2.85, cex=textsize)
mtext("Unrestricted\nCommunity", side=1, at=Fit.Plot[2], line=2.85, cex=textsize)
mtext("Restricted\nAlone", side=1, at=Fit.Plot[3], line=2.85, cex=textsize)

mtext("number of divisions", side=2, at=2500000, line=3.6, cex=textsize)

axis(2, las=1, at = c(1,10000,100000000,1000000000000), 
labels=c(expression(10^0),expression(10^4),expression(10^8),expression(10^12)),cex.axis=1.5)

#------------------------------------------------SI FIGURE 6