# MAGGIE DATA PROCESSING
# 

library(plyr)
library(lme4)
library(emmeans)

realized<-read.table("GrowthRates.csv",header=TRUE,sep=",")

realized$Genus<-relevel(realized$Genus,ref="Synechococcus")

#Model genera first, then individual strains, and do pairwise comparisons

growthmodel <- lmer(RealizedGR~CO2*Genus+(1|Strain),REML=FALSE,data=realized)
summary(growthmodel)
emmeans(growthmodel,pairwise~CO2|Genus)

growthmodel.strain <- lm(RealizedGR~CO2*Strain,data=realized)
summary(growthmodel.strain)
emmeans(growthmodel.strain,pairwise~CO2|Strain) 

growthmodel.exponential.genus<-lmer(ExponentialGR~CO2*Genus+(1|Strain),REML=FALSE,data=realized)
growthmodel.exponential.strain<-lm(ExponentialGR~CO2*Strain,data=realized)
summary(growthmodel.exponential.genus)
summary(growthmodel.exponential.strain)
emmeans(growthmodel.exponential.genus,pairwise~CO2|Genus)
emmeans(growthmodel.exponential.strain,pairwise~CO2|Strain)

growthmodel.lag.genus<-lmer(Lag~CO2*Genus+(1|Strain),REML=FALSE,data=realized)
growthmodel.lag.strain<-lm(Lag~CO2*Strain,data=realized)
summary(growthmodel.lag.genus)
summary(growthmodel.lag.strain)
emmeans(growthmodel.lag.genus,pairwise~CO2|Genus)
emmeans(growthmodel.lag.strain,pairwise~CO2|Strain)

#Aggregate data
table1<-ddply(realized,~Strain*CO2,summarise,mean=mean(RealizedGR),ex=mean(ExponentialGR),lag=mean(Lag),length=length(RealizedGR),pCO2=mean(CO2))

# Make plots

real.g.plot<-summary(emmeans(growthmodel,c("CO2","Genus")))
bar<-real.g.plot[c("emmean","lower.CL","upper.CL")]
bar<-as.data.frame(bar,row.names=paste(real.g.plot$CO2,real.g.plot$Genus,sep=" "))

tiff("Fig.Genus.Realized.tiff",compression="lzw",units="in",width=5,height=5,res=600)
coords<-barplot(bar$emmean,names.arg=row.names(bar),yaxt="n",xaxt="n",ylim=c(0,.6),col=c("White","Gray"))
title(ylab=expression(Realized~Growth~Rate~(d^{-1})),mgp=c(2.3,1,0))
axis(1,at=c((coords[1]+coords[2])/2,(coords[3]+coords[4])/2),labels=c(expression(italic("Synechococcus")),expression(italic("Prochlorococcus"))),lty=1)
axis(2,at=c(0,.1,.2,.3,.4,.5,.6))
box()
for(i in 1:4){segments(coords[i],bar$lower.CL[i],coords[i],bar$upper.CL[i])}
legend(3,0.58,legend=c("400 ppm","800 ppm"),fill=c("white","gray"))
text((coords[3]+coords[4])/2,.425,"***",cex=2)
dev.off()

ex.g.plot<-summary(emmeans(growthmodel.exponential.genus,c("CO2","Genus")))
bar<-ex.g.plot[c("emmean","lower.CL","upper.CL")]
bar<-as.data.frame(bar,row.names=paste(ex.g.plot$CO2,ex.g.plot$Genus,sep=" "))

tiff("Fig.Genus.Exponential.tiff",compression="lzw",units="in",width=5,height=5,res=600)
coords<-barplot(bar$emmean,names.arg=row.names(bar),yaxt="n",xaxt="n",ylim=c(0,.6),col=c("White","Gray"))
title(ylab=expression(Exponential~Growth~Rate~(d^{-1})),mgp=c(2.3,1,0))
axis(1,at=c((coords[1]+coords[2])/2,(coords[3]+coords[4])/2),labels=c(expression(italic("Synechococcus")),expression(italic("Prochlorococcus"))),lty=1)
axis(2,at=c(0,.1,.2,.3,.4,.5,.6))
box()
for(i in 1:4){segments(coords[i],bar$lower.CL[i],coords[i],bar$upper.CL[i])}
legend(3,0.58,legend=c("400 ppm","800 ppm"),fill=c("white","gray"))
dev.off()

lag.g.plot<-summary(emmeans(growthmodel.lag.genus,c("CO2","Genus")))
bar<-lag.g.plot[c("emmean","lower.CL","upper.CL")]
bar<-as.data.frame(bar,row.names=paste(lag.g.plot$CO2,lag.g.plot$Genus,sep=" "))

tiff("Fig.Genus.Lag.tiff",compression="lzw",units="in",width=5,height=5,res=600)
coords<-barplot(bar$emmean,names.arg=row.names(bar),ylab="Lag Time (d)",yaxt="n",xaxt="n",ylim=c(0,10),col=c("White","Gray"))
axis(1,at=c((coords[1]+coords[2])/2,(coords[3]+coords[4])/2),labels=c(expression(italic("Synechococcus")),expression(italic("Prochlorococcus"))),lty=1)
axis(2,at=c(0,2,4,6,8,10))
box()
for(i in 1:4){segments(coords[i],bar$lower.CL[i],coords[i],bar$upper.CL[i])}
legend(.5,9,legend=c("400 ppm","800 ppm"),fill=c("white","gray"))
text((coords[3]+coords[4])/2,9,"***",cex=2)
dev.off()

# Use emmeans output to produce strain-level tables.

#GRR

grr<-read.table("GRR.csv",header=TRUE,sep=",")

grr$Genus<-factor(grr$Genus,levels=c("Diazotroph","Synechococcus","Prochlorococcus"))
grr.nonatl<-subset(grr,Strain!="VOL3")
grr.noweirdos<-subset(grr.nonatl,Strain!="WH8102"&xgrr<1.4)

#compare pro vs. syn

rgrrmodel<-lm(rgrr~Genus,data=grr)
rgrrmodel.nonatl<-lm(rgrr~Genus,data=grr.nonatl)
rgrrmodel.noweirdos<-lm(rgrr~Genus,data=grr.noweirdos)
xgrrmodel<-lm(xgrr~Genus,data=grr)
xgrrmodel.nonatl<-lm(xgrr~Genus,data=grr.nonatl)
xgrrmodel.noweirdos<-lm(xgrr~Genus,data=grr.noweirdos)
xgrrmodel.proveverybody<-lm(xgrr~Flag,data=grr.noweirdos)

summary(rgrrmodel)
summary(rgrrmodel.nonatl)
summary(rgrrmodel.noweirdos)
summary(xgrrmodel)
summary(xgrrmodel.nonatl)
summary(xgrrmodel.noweirdos)
summary(xgrrmodel.proveverybody)

emmeans(rgrrmodel.noweirdos,pairwise~Genus)
emmeans(xgrrmodel.noweirdos,pairwise~Genus)
emmeans(xgrrmodel.proveverybody,pairwise~Flag)

wilcox.test(xgrr~Flag,data=grr)

morris<-subset(grr,Dataset=="Morris"&Strain!="VOL3")
morris$Genus<-factor(morris$Genus,c("Synechococcus","Prochlorococcus"))

tiff("Fig.RGRR.tiff",compression="lzw",units="in",width=5,height=5,res=600)
boxplot(rgrr~Genus,data=morris,ylim=c(0.4,1.6),xaxt="n",xlab="",ylab="Realized Growth Rate Response")
plot(jitter(c(1,1,1,1,1,1,2,2,2,2),factor=.2),c(subset(morris,Genus=="Synechococcus")$rgrr,subset(morris,Genus=="Prochlorococcus")$rgrr),xaxt="n",xlab="",ylab="Realized Growth Rate Response",xlim=c(0.5,2.5),ylim=c(0.4,1.6))
axis(side=1,at=c(1,2),tick=TRUE,labels=c(expression(italic("Synechococcus")),expression(italic("Prochlorococcus"))))
abline(h=1,lty=1)
segments(1,1.5,2,1.5)
text(1.5,1.6,"*",cex=2)
points(2,1.22,pch=16)
dev.off()

tiff("Fig.XGRR.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(jitter(as.numeric(subset(grr,Genus=="Synechococcus"&Dataset=="Morris")$Flag),factor=.5),subset(grr,Genus=="Synechococcus"&Dataset=="Morris")$xgrr,ylim=c(0.8,2.2),axes=FALSE,xlab="",ylab="",xlim=c(.5,2.5),pch=16)
par(new=TRUE)
plot(jitter(as.numeric(subset(grr,Genus=="Synechococcus"&Dataset=="Dutkiewicz")$Flag),factor=.5),subset(grr,Genus=="Synechococcus"&Dataset=="Dutkiewicz")$xgrr,ylim=c(0.8,2.2),axes=FALSE,xlab="",ylab="",xlim=c(.5,2.5),pch=17)
par(new=TRUE)
plot(jitter(as.numeric(subset(grr,Genus!="Synechococcus"&Dataset=="Morris")$Flag),factor=.5),subset(grr,Genus!="Synechococcus"&Dataset=="Morris")$xgrr,ylim=c(0.8,2.2),axes=FALSE,xlab="",ylab="",xlim=c(.5,2.5))
par(new=TRUE)
plot(jitter(as.numeric(subset(grr,Genus!="Synechococcus"&Dataset=="Dutkiewicz")$Flag),factor=.5),subset(grr,Genus!="Synechococcus"&Dataset=="Dutkiewicz")$xgrr,ylim=c(0.8,2.2),xaxt="n",xlab="",ylab="Exponential Growth Rate Response",xlim=c(.5,2.5),pch=2)
axis(side=1,at=c(1,2),tick=TRUE,labels=c("Cyanobacteria",expression(italic("Prochlorococcus"))))
abline(h=1,lty=1)
segments(1,2.15,2,2.15)
text(1.5,2.2,"*",cex=2)
dev.off()

# COMPETITION FIGURES

data<-read.table('CompDataProcessed.csv',header=TRUE,sep=",")
data$CO2<-as.factor(data$CO2)

model <- lm(Spro~CO2*InitProFreq,data=data)
summary(model)

# Rerunning analysis with the one outlier removed; only Maggie data

data.out<-read.table('CompDataProcessed.outlierremoved.csv',header=TRUE,sep=",")
data.out$CO2<-as.factor(data.out$CO2)
data.maggie<-subset(data.out,Experiment=="Maggie")

model.naive<-lm(Spro~CO2,data=data.out)
summary(model.naive)

emmeans(model.naive,pairwise~CO2)

# In a naive model, absolutely no effect of CO2 on treatment (p = 0.573)

# Let's plot naive Pro fitness vs. CO2 
tiff("Fig.Fitness.tiff",compression="lzw",units="in",width=5,height=5,res=600)
labels<-c("400 ppm","800 ppm")
boxplot(Spro~CO2,data=data.out,ylab="",xaxt="n",ylim=c(-1.5,1),col=c("white","gray"),xlab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[realized]~(d^{-1})),mgp=c(2.3,1,0))
axis(side=1,at=c(1,2),tick=TRUE,labels=labels)
abline(h=0)
abline(h=0,lty=1)
dev.off()

tiff("Fig.FitnessEx.tiff",compression="lzw",units="in",width=5,height=5,res=600)
labels<-c("400 ppm","800 ppm")
boxplot(SproEx~CO2,data=data.out,ylab="",xaxt="n",ylim=c(-1.5,1),col=c("white","gray"),xlab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[exponential]~(d^{-1})),mgp=c(2.3,1,0))
axis(side=1,at=c(1,2),tick=TRUE,labels=labels)
abline(h=0)
abline(h=0,lty=1)
dev.off()

#model with pro frequency included

model.maggie <- lm(Spro~CO2*InitProFreq,data=data.maggie)

summary(model.maggie)

# Main effect of CO2 is weak (0.268 +/- 0.156, p = 0.0887) and initial Pro frequency is  very strongly significant (-0.473 +/- 0.109, p = 3.47 x 10^-5), but interaction is not significant (p = 0.338).

#Now bring in ALL data

model.out <- lm(Spro~CO2*InitProFreq,data=data.out)
summary(model.out)

# With additional experimental data, the effect of CO2 and its interaction is completely gone.  Effect of initial Pro frequency is -0.444 +/- 0.093, p = 4.99 e -6.

#rerun without CO2

model.out.noco2<-lm(Spro~InitProFreq,data=data.out)
summary(model.out.noco2)

#main effect -0.437, p = 2.56e-8

#What about density dependence? Checking total density, plus initial pro and syn density.

model.total.dens <- lm(Spro~CO2*log10(InitDens), data=data.out)
summary(model.total.dens)

#Significant negative impact of initial density, borderline effects of CO2 and interaction with density

model.total.dens.noco2 <- lm(Spro~log10(InitDens), data=data.out)
summary(model.total.dens.noco2)

#NOT significant.

model.pro.dens <- lm(Spro~CO2*log10(InitPro), data=data.out)
summary(model.pro.dens)

#Strong negative impact of Pro initial density

model.pro.dens.noco2<-lm(Spro~log10(InitPro),data=data.out)
summary(model.pro.dens.noco2)

#Strong negative impact of Pro initial density

model.syn.dens <- lm(Spro~CO2*log10(InitSyn),data=data.out)
summary(model.syn.dens)

#Negative impact of CO2, counteracted by Initial Syn Density

model.syn.dens.noco2 <- lm(Spro~log10(InitSyn),data=data.out)
summary(model.syn.dens.noco2)

#Positive impact of initial syn density

model.full.dens <- lm(Spro~CO2*log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.full.dens)

model.full.dens.noco2 <- lm(Spro~log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.full.dens.noco2)

#Both full models; positive impact of initial Syn density and negative impact of initial Pro density

BIC(model.out,model.out.noco2,model.total.dens,model.total.dens.noco2,model.pro.dens,model.pro.dens.noco2,model.syn.dens,model.syn.dens.noco2,model.full.dens,model.full.dens.noco2)

#Best model is model.full.dens.noco2, then the no CO2 frequency dependent model. Worst models are the total density models and Syn density models.

# Predicted equilibrium isocline from final full.dens model:
#spro=logSyn*b+logPro*c+a
#0=0.3053logSyn-.186logPro-.6421
#.3053 logSyn = .186 logPro + .642
#logSyn = .609 logPro + 2.103

isocline.dens<-c(-coef(model.full.dens.noco2)[1]/coef(model.full.dens.noco2)[2],-coef(model.full.dens.noco2)[3]/coef(model.full.dens.noco2)[2])

model.out.ProGR<-lm(ProGR~CO2*InitProFreq,data=data.out)
summary(model.out.ProGR)
slope.grfreq.low<-coef(model.out.ProGR)[3]
slope.grfreq.high<-coef(model.out.ProGR)[3]+coef(model.out.ProGR)[4]
int.grfreq.low<-coef(model.out.ProGR)[1]
int.grfreq.high<-coef(model.out.ProGR)[1]+coef(model.out.ProGR)[2]

model.out.noco2.ProGR<-lm(ProGR~InitProFreq,data=data.out)
summary(model.out.noco2.ProGR)

model.total.dens.ProGR<- lm(ProGR~CO2*log10(InitDens), data=data.out)
summary(model.total.dens.ProGR)

model.total.dens.noco2.ProGR<- lm(ProGR~log10(InitDens), data=data.out)
summary(model.total.dens.noco2.ProGR)

model.pro.dens.ProGR<- lm(ProGR~CO2*log10(InitPro), data=data.out)
summary(model.pro.dens.ProGR)
slope.dens.low<-coef(model.pro.dens.ProGR)[3]
slope.dens.high<-coef(model.pro.dens.ProGR)[3]+coef(model.pro.dens.ProGR)[4]
int.dens.low<-coef(model.pro.dens.ProGR)[1]
int.dens.high<-coef(model.pro.dens.ProGR)[1]+coef(model.pro.dens.ProGR)[2]

model.pro.dens.noco2.ProGR<-lm(ProGR~log10(InitPro),data=data.out)
summary(model.pro.dens.noco2.ProGR)

model.syn.dens.ProGR<- lm(ProGR~CO2*log10(InitSyn),data=data.out)
summary(model.syn.dens.ProGR)

model.syn.dens.noco2.ProGR<- lm(ProGR~log10(InitSyn),data=data.out)
summary(model.syn.dens.noco2.ProGR)

model.full.dens.ProGR<- lm(ProGR~CO2*log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.full.dens.ProGR)

model.full.dens.noco2.ProGR<- lm(ProGR~log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.full.dens.noco2.ProGR)

BIC(model.out.ProGR,model.out.noco2.ProGR,model.total.dens.ProGR,model.total.dens.noco2.ProGR,model.pro.dens.ProGR,model.pro.dens.noco2.ProGR,model.syn.dens.ProGR,model.syn.dens.noco2.ProGR,model.full.dens.ProGR,model.full.dens.noco2.ProGR)

#Best model is the Pro Density model with or without CO2; second best are the freq. dependence models. All are DRAMATICALLY better than relative fitness models. CO2 is always significantly negative, and counteracted by the density term (only marginal in the case of Syn density) -- except in "full dens" model.

model.out.SynGR<-lm(SynGR~CO2*InitProFreq,data=data.out)
summary(model.out.SynGR)

model.out.noco2.SynGR<-lm(SynGR~InitProFreq,data=data.out)
summary(model.out.noco2.SynGR)

model.total.dens.SynGR<- lm(SynGR~CO2*log10(InitDens), data=data.out)
summary(model.total.dens.SynGR)

model.total.dens.noco2.SynGR<- lm(SynGR~log10(InitDens), data=data.out)
summary(model.total.dens.noco2.SynGR)

model.pro.dens.SynGR<- lm(SynGR~CO2*log10(InitPro), data=data.out)
summary(model.pro.dens.SynGR)

model.pro.dens.noco2.SynGR<-lm(SynGR~log10(InitPro),data=data.out)
summary(model.pro.dens.noco2.SynGR)

model.syn.dens.SynGR<- lm(SynGR~CO2*log10(InitSyn),data=data.out)
summary(model.syn.dens.SynGR)

model.syn.dens.noco2.SynGR<- lm(SynGR~log10(InitSyn),data=data.out)
summary(model.syn.dens.noco2.SynGR)
syngr.syndens.int=coef(model.syn.dens.noco2.SynGR)[1]
syngr.syndens.slope=coef(model.syn.dens.noco2.SynGR)[2]

model.full.dens.SynGR<- lm(SynGR~CO2*log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.full.dens.SynGR)

model.full.dens.noco2.SynGR<- lm(SynGR~log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.full.dens.noco2.SynGR)

BIC(model.out.SynGR,model.out.noco2.SynGR,model.total.dens.SynGR,model.total.dens.noco2.SynGR,model.pro.dens.SynGR,model.pro.dens.noco2.SynGR,model.syn.dens.SynGR,model.syn.dens.noco2.SynGR,model.full.dens.SynGR,model.full.dens.noco2.SynGR)

#Best model is Syn density, no CO2. CO2 is negative and counteracted by frequency term in freq dependence model, but not in any other ones.  Syn density is negative, Pro density has no impact.

***** USE EXPONENTIAL GR FITNESSES; ARE ANY DIFFERENT?
model.exgr.maggie <- lm(SproEx~CO2*InitProFreq,data=data.maggie)

summary(model.exgr.maggie)

model.exgr.out <- lm(SproEx~CO2*InitProFreq,data=data.out)
summary(model.exgr.out)

model.exgr.out.noco2<-lm(SproEx~InitProFreq,data=data.out)
summary(model.exgr.out.noco2)

model.exgr.total.dens <- lm(SproEx~CO2*log10(InitDens), data=data.out)
summary(model.exgr.total.dens)

model.exgr.total.dens.noco2 <- lm(SproEx~log10(InitDens), data=data.out)
summary(model.exgr.total.dens.noco2)

model.exgr.pro.dens <- lm(SproEx~CO2*log10(InitPro), data=data.out)
summary(model.exgr.pro.dens)

model.exgr.pro.dens.noco2<-lm(SproEx~log10(InitPro),data=data.out)
summary(model.exgr.pro.dens.noco2)

model.exgr.syn.dens <- lm(SproEx~CO2*log10(InitSyn),data=data.out)
summary(model.exgr.syn.dens)

model.exgr.syn.dens.noco2 <- lm(SproEx~log10(InitSyn),data=data.out)
summary(model.exgr.syn.dens.noco2)

model.exgr.full.dens <- lm(SproEx~CO2*log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.exgr.full.dens)

model.exgr.full.dens.noco2 <- lm(SproEx~log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.exgr.full.dens.noco2)

BIC(model.exgr.out,model.exgr.out.noco2,model.exgr.total.dens,model.exgr.total.dens.noco2,model.exgr.pro.dens,model.exgr.pro.dens.noco2,model.exgr.syn.dens,model.exgr.syn.dens.noco2,model.exgr.full.dens,model.exgr.full.dens.noco2)

#Best model.exgr is the no CO2 frequency dependent model. Worst models are the total density models and Syn density models.

model.exgr.out.ProExGR<-lm(ProExGR~CO2*InitProFreq,data=data.out)
summary(model.exgr.out.ProExGR)

model.exgr.out.noco2.ProExGR<-lm(ProExGR~InitProFreq,data=data.out)
summary(model.exgr.out.noco2.ProExGR)

model.exgr.total.dens.ProExGR<- lm(ProExGR~CO2*log10(InitDens), data=data.out)
summary(model.exgr.total.dens.ProExGR)

model.exgr.total.dens.noco2.ProExGR<- lm(ProExGR~log10(InitDens), data=data.out)
summary(model.exgr.total.dens.noco2.ProExGR)

model.exgr.pro.dens.ProExGR<- lm(ProExGR~CO2*log10(InitPro), data=data.out)
summary(model.exgr.pro.dens.ProExGR)

model.exgr.pro.dens.noco2.ProExGR<-lm(ProExGR~log10(InitPro),data=data.out)
summary(model.exgr.pro.dens.noco2.ProExGR)

model.exgr.syn.dens.ProExGR<- lm(ProExGR~CO2*log10(InitSyn),data=data.out)
summary(model.exgr.syn.dens.ProExGR)

model.exgr.syn.dens.noco2.ProExGR<- lm(ProExGR~log10(InitSyn),data=data.out)
summary(model.exgr.syn.dens.noco2.ProExGR)

model.exgr.full.dens.ProExGR<- lm(ProExGR~CO2*log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.exgr.full.dens.ProExGR)

model.exgr.full.dens.noco2.ProExGR<- lm(ProExGR~log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.exgr.full.dens.noco2.ProExGR)

BIC(model.exgr.out.ProExGR,model.exgr.out.noco2.ProExGR,model.exgr.total.dens.ProExGR,model.exgr.total.dens.noco2.ProExGR,model.exgr.pro.dens.ProExGR,model.exgr.pro.dens.noco2.ProExGR,model.exgr.syn.dens.ProExGR,model.exgr.syn.dens.noco2.ProExGR,model.exgr.full.dens.ProExGR,model.exgr.full.dens.noco2.ProExGR)

#Best model.exgr is the Pro Density model.exgr without CO2. CO2 is always significantly negative, and counteracted by the density term -- except in "full dens" and "syn dens" models.

model.exgr.out.SynExGR<-lm(SynExGR~CO2*InitProFreq,data=data.out)
summary(model.exgr.out.SynExGR)

model.exgr.out.noco2.SynExGR<-lm(SynExGR~InitProFreq,data=data.out)
summary(model.exgr.out.noco2.SynExGR)

model.exgr.total.dens.SynExGR<- lm(SynExGR~CO2*log10(InitDens), data=data.out)
summary(model.exgr.total.dens.SynExGR)

model.exgr.total.dens.noco2.SynExGR<- lm(SynExGR~log10(InitDens), data=data.out)
summary(model.exgr.total.dens.noco2.SynExGR)

model.exgr.pro.dens.SynExGR<- lm(SynExGR~CO2*log10(InitPro), data=data.out)
summary(model.exgr.pro.dens.SynExGR)

model.exgr.pro.dens.noco2.SynExGR<-lm(SynExGR~log10(InitPro),data=data.out)
summary(model.exgr.pro.dens.noco2.SynExGR)

model.exgr.syn.dens.SynExGR<- lm(SynExGR~CO2*log10(InitSyn),data=data.out)
summary(model.exgr.syn.dens.SynExGR)

model.exgr.syn.dens.noco2.SynExGR<- lm(SynExGR~log10(InitSyn),data=data.out)
summary(model.exgr.syn.dens.noco2.SynExGR)
syngr.syndens.int=coef(model.exgr.syn.dens.noco2.SynExGR)[1]

model.exgr.full.dens.SynExGR<- lm(SynExGR~CO2*log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.exgr.full.dens.SynExGR)

model.exgr.full.dens.noco2.SynExGR<- lm(SynExGR~log10(InitSyn)+log10(InitPro),data=data.out)
summary(model.exgr.full.dens.noco2.SynExGR)

BIC(model.exgr.out.SynExGR,model.exgr.out.noco2.SynExGR,model.exgr.total.dens.SynExGR,model.exgr.total.dens.noco2.SynExGR,model.exgr.pro.dens.SynExGR,model.exgr.pro.dens.noco2.SynExGR,model.exgr.syn.dens.SynExGR,model.exgr.syn.dens.noco2.SynExGR,model.exgr.full.dens.SynExGR,model.exgr.full.dens.noco2.SynExGR)

#best model is syn dens no co2 model.  Same conclusions as before.

#*******

#Okay, so what do I want to plot? 1. Def frequency dependence. 2. Also density dependence of Pro growth rate, with 3. density dependence of Pro fitness as a supp figure.  Also convergence on equilibrium 4. frequency and 5. concentration and 6. fitness

# 1. Let's plot the regression of fitness on Initial Pro Frequency.

data.jeff<-subset(data.out,Experiment!="Maggie")

tiff("Fig.Spro.FreqDep.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$InitProFreq,subset(data.maggie,CO2=="0")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(0,1),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$InitProFreq,subset(data.maggie,CO2=="1")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(0,1),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$InitProFreq,subset(data.jeff,CO2=="0")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(0,1),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$InitProFreq,subset(data.jeff,CO2=="1")$Spro,ylim=c(-1.5,1),xlim=c(0,1),pch=17,xlab=expression(Initial~Frequency~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[realized]~(d^{-1})),mgp=c(2.3,1,0))
legend(.55,.95,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)

#Estimate equilibrium points:

equil<--coef(model.out.noco2)[1]/coef(model.out.noco2)[2]
abline(model.out.noco2,lty=2)
abline(v=equil,lty=2)

dev.off()

#Based on the model, they should coexist at 17% under both CO2 conditions.

#2. Init Pro Dens vs. Pro growth rate
tiff("Fig.ProGR.ProDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitPro),subset(data.maggie,CO2=="0")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitPro),subset(data.maggie,CO2=="1")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitPro),subset(data.jeff,CO2=="0")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitPro),subset(data.jeff,CO2=="1")$ProGR,ylim=c(-.5,1.5),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(Realized~Growth~Rate~of~italic("Prochlorococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(5,1.5,c("400 ppm","800 ppm"),pch=c(1,2),lty=c(3,2),cex=.8)
abline(h=0,lty=1)

abline(int.dens.low,slope.dens.low,lty=3)
abline(int.dens.high,slope.dens.high,lty=2)

dev.off()

#3. Init Pro Dens vs. Pro fitness
tiff("Fig.Spro.ProDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitPro),subset(data.maggie,CO2=="0")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1.5),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitPro),subset(data.maggie,CO2=="1")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1.5),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitPro),subset(data.jeff,CO2=="0")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1.5),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitPro),subset(data.jeff,CO2=="1")$Spro,ylim=c(-1.5,1.5),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[realized]~(d^{-1})),mgp=c(2.3,1,0))
legend(5,1.5,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)

#Estimate equilibrium points:

pro.slope<-coef(model.pro.dens.noco2)[2]
pro.int<-coef(model.pro.dens.noco2)[1]
pro.equil<--pro.int/pro.slope

abline(pro.int,pro.slope,lty=2)
abline(v=pro.equil,lty=2)

dev.off()

#4. Plot convergence of Frequency on transfer #
tiff("Fig.Spro.Freq.Convergence.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$Transfer-0.1,subset(data.maggie,CO2=="0")$InitProFreq,axes=FALSE,ann=FALSE,ylim=c(0,1),xlim=c(1,7),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$Transfer+0.1,subset(data.maggie,CO2=="1")$InitProFreq,axes=FALSE,ann=FALSE,ylim=c(0,1),xlim=c(1,7),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$Transfer-0.1,subset(data.jeff,CO2=="0")$InitProFreq,axes=FALSE,ann=FALSE,ylim=c(0,1),xlim=c(1,7),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$Transfer+0.1,subset(data.jeff,CO2=="1")$InitProFreq,ylim=c(0,1),xlim=c(1,7),pch=17,xlab="Transfer",ylab=expression(Initial~Frequency~of~italic("Prochlorococcus")))
legend(4.5,.9,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=equil,lty=2)
dev.off()

#Over time, fitness converges on 0 and freq approaches that predicted by the regression analysis, suggesting equilibrium.

#5. make bubble plot with bubble size and color representing transfer #
tiff("Fig.BubblePro.tiff",compression="lzw",units="in",width=5,height=5,res=600)
symbols(log10(data.out$InitPro),log10(data.out$InitSyn),circles=data.out$Transfer,inches=.2,fg=data.out$Transfer,xlab=expression(Log~Initial~Density~of~italic("Prochlorococcus")),ylab=expression(Log~Initial~Density~of~italic("Synechococcus")),xlim=c(2,6.5),ylim=c(2,6))

legend(5.8,4.5,c(1:6),col=c(1:6),pch=1,pt.cex=c(1:6)/1.5,bty="n",y.intersp=c(1,1,1.1,1.2,1.3,1.42),cex=.8,title="Transfer",x.intersp=2)

# Determine isocline where Pro is 16.6% of the overall population:
#Pro/.166=total
#total = syn + pro
#Pro/.166=syn+pro
#Pro/.166-pro=syn
#pro/.166-.166pro/.166=syn
#(pro-.166pro)/.166=syn
#log(pro-.166pro)-log(.166)=log(syn)
#log(pro(1-.166))-log(.166)=log(syn)
#log(pro)+log(1-.166)-log(.166)=log(syn)
#log(pro)+.701=log(syn)
#intercept = .701, slope = 1

abline(.701,1,lty=2)
#abline(isocline.dens,lty=1)

dev.off()

# 6. Let's plot the Pro fitness as a function of transfer number.
tiff("Fig.Spro.Convergence.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$Transfer-0.1,subset(data.maggie,CO2=="0")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(1,7),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$Transfer+0.1,subset(data.maggie,CO2=="1")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(1,7),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$Transfer-0.1,subset(data.jeff,CO2=="0")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(1,7),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$Transfer+0.1,subset(data.jeff,CO2=="1")$Spro,ylim=c(-1.5,1),xlim=c(1,7),pch=17,xlab="Transfer",ylab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[realized]~(d^{-1})),mgp=c(2.3,1,0))
legend(4,-1,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

#Supp: Now let's plot high and low CO2 vs. Init Syn Density
tiff("Fig.Spro.SynDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
par("mar"=c(5,4,4,1)+.1)
plot(log10(subset(data.maggie,CO2=="0")$InitSyn),subset(data.maggie,CO2=="0")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitSyn),subset(data.maggie,CO2=="1")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitSyn),subset(data.jeff,CO2=="0")$Spro,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitSyn),subset(data.jeff,CO2=="1")$Spro,ylim=c(-1.5,1),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Synechococcus")),ylab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[realized]~(d^{-1})),mgp=c(2.3,1,0))
legend(2,.95,c("400 ppm","800 ppm"),pch=c(1,2),lty=c(3,2),cex=.8)
abline(h=0,lty=1)

#Estimate equilibrium points:

slope.low<-coef(model.syn.dens)[3]
slope.high<-coef(model.syn.dens)[3]+coef(model.syn.dens)[4]
int.low<-coef(model.syn.dens)[1]
int.high<-coef(model.syn.dens)[1]+coef(model.syn.dens)[2]
equil.high<--int.high/slope.high

abline(int.low,slope.low,lty=3)
abline(int.high,slope.high,lty=2)
abline(v=equil.high,lty=2)

dev.off()

#Supp. Init Freq vs. Pro GR
tiff("Fig.ProGR.FreqDep.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot((subset(data.maggie,CO2=="0")$InitProFreq),subset(data.maggie,CO2=="0")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(0,1),pch=1)
par(new=TRUE)
plot((subset(data.maggie,CO2=="1")$InitProFreq),subset(data.maggie,CO2=="1")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(0,1),pch=2)
par(new=TRUE)
plot((subset(data.jeff,CO2=="0")$InitProFreq),subset(data.jeff,CO2=="0")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(0,1),pch=16)
par(new=TRUE)
plot((subset(data.jeff,CO2=="1")$InitProFreq),subset(data.jeff,CO2=="1")$ProGR,ylim=c(-.5,1.5),xlim=c(0,1),pch=17,xlab=expression(Initial~Frequency~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(Realized~Growth~Rate~of~italic("Prochlorococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(.7,1.5,c("400 ppm","800 ppm"),pch=c(1,2),lty=c(3,2),cex=.8)
abline(h=0,lty=1)

abline(int.grfreq.low,slope.grfreq.low,lty=3)
abline(int.grfreq.high,slope.grfreq.high,lty=2)
abline(v=pro.equil,lty=2)

dev.off()

#Supp. Init Syn vs. Pro GR
tiff("Fig.ProGR.SynDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitSyn),subset(data.maggie,CO2=="0")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitSyn),subset(data.maggie,CO2=="1")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitSyn),subset(data.jeff,CO2=="0")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitSyn),subset(data.jeff,CO2=="1")$ProGR,ylim=c(-.5,1.5),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Synechococcus")),ylab="")
title(ylab=expression(Realized~Growth~Rate~of~italic("Prochlorococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(5,1.5,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)

dev.off()

#Supp. Init Freq vs. Syn GR
tiff("Fig.SynGR.FreqDep.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot((subset(data.maggie,CO2=="0")$InitProFreq),subset(data.maggie,CO2=="0")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(0,1),pch=1)
par(new=TRUE)
plot((subset(data.maggie,CO2=="1")$InitProFreq),subset(data.maggie,CO2=="1")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(0,1),pch=2)
par(new=TRUE)
plot((subset(data.jeff,CO2=="0")$InitProFreq),subset(data.jeff,CO2=="0")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(0,1),pch=16)
par(new=TRUE)
plot((subset(data.jeff,CO2=="1")$InitProFreq),subset(data.jeff,CO2=="1")$SynGR,ylim=c(-.5,2),xlim=c(0,1),pch=17,xlab=expression(Initial~Frequency~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(Realized~Growth~Rate~of~italic("Synechococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(0,2,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)

dev.off()

#Supp. Init Syn vs. Syn GR
tiff("Fig.SynGR.SynDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitSyn),subset(data.maggie,CO2=="0")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitSyn),subset(data.maggie,CO2=="1")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitSyn),subset(data.jeff,CO2=="0")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitSyn),subset(data.jeff,CO2=="1")$SynGR,ylim=c(-.5,2),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Synechococcus")),ylab="")
title(ylab=expression(Realized~Growth~Rate~of~italic("Synechococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(2,2,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
abline(syngr.syndens.int,syngr.syndens.slope,lty=2)
dev.off()

#Supp. Init Pro vs. Syn GR
tiff("Fig.SynGR.ProDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitPro),subset(data.maggie,CO2=="0")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitPro),subset(data.maggie,CO2=="1")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitPro),subset(data.jeff,CO2=="0")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitPro),subset(data.jeff,CO2=="1")$SynGR,ylim=c(-.5,2),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(Realized~Growth~Rate~of~italic("Synechococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(2,2,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

# Supp. ProGR vs. transfer number.
tiff("Fig.ProGR.Convergence.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$Transfer-0.1,subset(data.maggie,CO2=="0")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(1,7),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$Transfer+0.1,subset(data.maggie,CO2=="1")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(1,7),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$Transfer-0.1,subset(data.jeff,CO2=="0")$ProGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(1,7),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$Transfer+0.1,subset(data.jeff,CO2=="1")$ProGR,ylim=c(-.5,1.5),xlim=c(1,7),pch=17,xlab="Transfer",ylab="")
title(ylab=expression(Realized~Growth~Rate~of~italic("Prochlorococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(5,1.5,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

# Supp. SynGR vs. transfer number.
tiff("Fig.SynGR.Convergence.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$Transfer-0.1,subset(data.maggie,CO2=="0")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(1,7),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$Transfer+0.1,subset(data.maggie,CO2=="1")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(1,7),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$Transfer-0.1,subset(data.jeff,CO2=="0")$SynGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(1,7),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$Transfer+0.1,subset(data.jeff,CO2=="1")$SynGR,ylim=c(-.5,2),xlim=c(1,7),pch=17,xlab="Transfer",ylab="")
title(ylab=expression(Realized~Growth~Rate~of~italic("Synechococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(5,2,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

# Supp. InitPro vs. transfer number.
tiff("Fig.InitPro.Convergence.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$Transfer-0.1,log10(subset(data.maggie,CO2=="0")$InitPro),axes=FALSE,ann=FALSE,ylim=c(2,6),xlim=c(1,7),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$Transfer+0.1,log10(subset(data.maggie,CO2=="1")$InitPro),axes=FALSE,ann=FALSE,ylim=c(2,6),xlim=c(1,7),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$Transfer-0.1,log10(subset(data.jeff,CO2=="0")$InitPro),axes=FALSE,ann=FALSE,ylim=c(2,6),xlim=c(1,7),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$Transfer+0.1,log10(subset(data.jeff,CO2=="1")$InitPro),ylim=c(2,6),xlim=c(1,7),pch=17,xlab="Transfer",ylab="")
title(ylab=expression(Log~Initial~Density~of~italic("Prochlorococcus")~(cells~mL^{-1})),mgp=c(2.3,1,0))
legend(5,6,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

# Supp. InitSyn vs. transfer number.
tiff("Fig.InitSyn.Convergence.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$Transfer-0.1,log10(subset(data.maggie,CO2=="0")$InitSyn),axes=FALSE,ann=FALSE,ylim=c(2,6),xlim=c(1,7),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$Transfer+0.1,log10(subset(data.maggie,CO2=="1")$InitSyn),axes=FALSE,ann=FALSE,ylim=c(2,6),xlim=c(1,7),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$Transfer-0.1,log10(subset(data.jeff,CO2=="0")$InitSyn),axes=FALSE,ann=FALSE,ylim=c(2,6),xlim=c(1,7),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$Transfer+0.1,log10(subset(data.jeff,CO2=="1")$InitSyn),ylim=c(2,6),xlim=c(1,7),pch=17,xlab="Transfer",ylab="")
title(ylab=expression(Log~Initial~Density~of~italic("Synechococcus")~(cells~mL^{-1})),mgp=c(2.3,1,0))
legend(5,6,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

#*********

# 1ex. Let's plot the regression of Ex fitness on Initial Pro Frequency.

tiff("Fig.SproEx.FreqDep.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$InitProFreq,subset(data.maggie,CO2=="0")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(0,1),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$InitProFreq,subset(data.maggie,CO2=="1")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(0,1),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$InitProFreq,subset(data.jeff,CO2=="0")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(0,1),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$InitProFreq,subset(data.jeff,CO2=="1")$SproEx,ylim=c(-1.5,1),xlim=c(0,1),pch=17,xlab=expression(Initial~Frequency~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[exponential]~(d^{-1})),mgp=c(2.3,1,0))
legend(.55,.95,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)

#Estimate equilibrium points:

equil<--coef(model.exgr.out.noco2)[1]/coef(model.exgr.out.noco2)[2]
abline(model.exgr.out.noco2,lty=2)
abline(v=equil,lty=2)

dev.off()

#Based on the no interaction assumption, they should coexist at 0.4% under both CO2 conditions.

#2ex. Init Pro Dens vs. Pro Ex growth rate

tiff("Fig.ProExGR.ProDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitPro),subset(data.maggie,CO2=="0")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitPro),subset(data.maggie,CO2=="1")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitPro),subset(data.jeff,CO2=="0")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitPro),subset(data.jeff,CO2=="1")$ProExGR,ylim=c(-.5,1.5),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(Exponential~Growth~Rate~of~italic("Prochlorococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(5,1.5,c("400 ppm","800 ppm"),pch=c(1,2),lty=c(3,2),cex=.8)
abline(h=0,lty=1)

int.dens.ex.low<-coef(model.exgr.pro.dens.ProExGR)[1]
int.dens.ex.high<-coef(model.exgr.pro.dens.ProExGR)[1]+coef(model.exgr.pro.dens.ProExGR)[2]
slope.dens.ex.low<-coef(model.exgr.pro.dens.ProExGR)[3]
slope.dens.ex.high<-coef(model.exgr.pro.dens.ProExGR)[3]+coef(model.exgr.pro.dens.ProExGR)[4]

abline(int.dens.ex.low,slope.dens.ex.low,lty=3)
abline(int.dens.ex.high,slope.dens.ex.high,lty=2)

dev.off()

#3ex. Init Pro Dens vs. Pro Ex fitness
tiff("Fig.SproEx.ProDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitPro),subset(data.maggie,CO2=="0")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1.5),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitPro),subset(data.maggie,CO2=="1")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1.5),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitPro),subset(data.jeff,CO2=="0")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1.5),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitPro),subset(data.jeff,CO2=="1")$SproEx,ylim=c(-1.5,1.5),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[exponential]~(d^{-1})),mgp=c(2.3,1,0))
legend(4,1.5,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)

#Estimate equilibrium points:

pro.slope<-coef(model.exgr.pro.dens.noco2)[2]
pro.int<-coef(model.exgr.pro.dens.noco2)[1]
pro.equil<--pro.int/pro.slope

abline(pro.int,pro.slope,lty=2)
abline(v=pro.equil,lty=2)

dev.off()

# 6. Let's plot the Pro fitness Ex as a function of transfer number.
tiff("Fig.SproEx.Convergence.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$Transfer-0.1,subset(data.maggie,CO2=="0")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(1,7),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$Transfer+0.1,subset(data.maggie,CO2=="1")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(1,7),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$Transfer-0.1,subset(data.jeff,CO2=="0")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(1,7),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$Transfer+0.1,subset(data.jeff,CO2=="1")$SproEx,ylim=c(-1.5,1),xlim=c(1,7),pch=17,xlab="Transfer",ylab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[exponential]~(d^{-1})),mgp=c(2.3,1,0))
legend(4,-1,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

#Supp: Now let's plot high and low CO2 Exponential Fitness vs. Init Syn Density
tiff("Fig.SproEx.SynDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
par("mar"=c(5,4,4,1)+.1)
plot(log10(subset(data.maggie,CO2=="0")$InitSyn),subset(data.maggie,CO2=="0")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitSyn),subset(data.maggie,CO2=="1")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitSyn),subset(data.jeff,CO2=="0")$SproEx,axes=FALSE,ann=FALSE,ylim=c(-1.5,1),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitSyn),subset(data.jeff,CO2=="1")$SproEx,ylim=c(-1.5,1),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Synechococcus")),ylab="")
title(ylab=expression(italic("Prochlorococcus")~italic("s")[exponential]~(d^{-1})),mgp=c(2.3,1,0))
legend(2,.95,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)

dev.off()

#Supp. Init Freq vs. Pro GREx
tiff("Fig.ProExGR.FreqDep.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot((subset(data.maggie,CO2=="0")$InitProFreq),subset(data.maggie,CO2=="0")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(0,1),pch=1)
par(new=TRUE)
plot((subset(data.maggie,CO2=="1")$InitProFreq),subset(data.maggie,CO2=="1")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(0,1),pch=2)
par(new=TRUE)
plot((subset(data.jeff,CO2=="0")$InitProFreq),subset(data.jeff,CO2=="0")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(0,1),pch=16)
par(new=TRUE)
plot((subset(data.jeff,CO2=="1")$InitProFreq),subset(data.jeff,CO2=="1")$ProExGR,ylim=c(-.5,1.5),xlim=c(0,1),pch=17,xlab=expression(Initial~Frequency~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(Exponential~Growth~Rate~of~italic("Prochlorococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(.7,1.5,c("400 ppm","800 ppm"),pch=c(1,2),lty=c(3,2),cex=.8)
abline(h=0,lty=1)

int.grfreq.ex.low<-coef(model.exgr.out.ProExGR)[1]
int.grfreq.ex.high<-coef(model.exgr.out.ProExGR)[1]+coef(model.exgr.out.ProExGR)[2]
slope.grfreq.ex.low<-coef(model.exgr.out.ProExGR)[3]
slope.grfreq.ex.high<-coef(model.exgr.out.ProExGR)[3]+coef(model.exgr.out.ProExGR)[4]

abline(int.grfreq.ex.low,slope.grfreq.ex.low,lty=3)
abline(int.grfreq.ex.high,slope.grfreq.ex.high,lty=2)

dev.off()

#Supp. Init Syn vs. Pro GREx
tiff("Fig.ProExGR.SynDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitSyn),subset(data.maggie,CO2=="0")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitSyn),subset(data.maggie,CO2=="1")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitSyn),subset(data.jeff,CO2=="0")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitSyn),subset(data.jeff,CO2=="1")$ProExGR,ylim=c(-.5,1.5),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Synechococcus")),ylab="")
title(ylab=expression(Exponential~Growth~Rate~of~italic("Prochlorococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(5,1.5,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)

int.syn.ex<-coef(model.exgr.syn.dens.noco2.ProExGR)[1]
slope.syn.ex<-coef(model.exgr.syn.dens.noco2.ProExGR)[2]
abline(int.syn.ex,slope.syn.ex,lty=2)

dev.off()

#Supp. Init Freq vs. Syn GREx
tiff("Fig.SynExGR.FreqDep.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot((subset(data.maggie,CO2=="0")$InitProFreq),subset(data.maggie,CO2=="0")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(0,1),pch=1)
par(new=TRUE)
plot((subset(data.maggie,CO2=="1")$InitProFreq),subset(data.maggie,CO2=="1")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(0,1),pch=2)
par(new=TRUE)
plot((subset(data.jeff,CO2=="0")$InitProFreq),subset(data.jeff,CO2=="0")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(0,1),pch=16)
par(new=TRUE)
plot((subset(data.jeff,CO2=="1")$InitProFreq),subset(data.jeff,CO2=="1")$SynExGR,ylim=c(-.5,2),xlim=c(0,1),pch=17,xlab=expression(Initial~Frequency~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(Exponential~Growth~Rate~of~italic("Synechococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(0,2,c("400 ppm","800 ppm"),pch=c(1,2),lty=c(3,2),cex=.8)
abline(h=0,lty=1)

int.synfreq.ex.low<-coef(model.exgr.out.SynExGR)[1]
int.synfreq.ex.high<-coef(model.exgr.out.SynExGR)[1]+coef(model.exgr.out.SynExGR)[2]
slope.synfreq.ex.low<-coef(model.exgr.out.SynExGR)[3]
slope.synfreq.ex.high<-coef(model.exgr.out.SynExGR)[3]+coef(model.exgr.out.SynExGR)[4]

abline(int.synfreq.ex.high,slope.synfreq.ex.high,lty=2)
abline(int.synfreq.ex.low,slope.synfreq.ex.low,lty=3)

dev.off()

#Supp. Init Syn vs. Syn GREx
tiff("Fig.SynExGR.SynDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitSyn),subset(data.maggie,CO2=="0")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitSyn),subset(data.maggie,CO2=="1")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitSyn),subset(data.jeff,CO2=="0")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitSyn),subset(data.jeff,CO2=="1")$SynExGR,ylim=c(-.5,2),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Synechococcus")),ylab="")
title(ylab=expression(Exponential~Growth~Rate~of~italic("Synechococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(2,2,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
abline(syngr.syndens.int,syngr.syndens.slope,lty=2)

int.synsyn.ex<-coef(model.exgr.syn.dens.noco2.SynExGR)[1]
slope.synsyn.ex<-coef(model.exgr.syn.dens.noco2.SynExGR)[2]

abline(int.synsyn.ex,slope.synsyn.ex,lty=2)

dev.off()

#Supp. Init Pro vs. Syn GREx
tiff("Fig.SynExGR.ProDens.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(log10(subset(data.maggie,CO2=="0")$InitPro),subset(data.maggie,CO2=="0")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=1)
par(new=TRUE)
plot(log10(subset(data.maggie,CO2=="1")$InitPro),subset(data.maggie,CO2=="1")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=2)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="0")$InitPro),subset(data.jeff,CO2=="0")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(2,6),pch=16)
par(new=TRUE)
plot(log10(subset(data.jeff,CO2=="1")$InitPro),subset(data.jeff,CO2=="1")$SynExGR,ylim=c(-.5,2),xlim=c(2,6),pch=17,xlab=expression(Log~Initial~Density~of~italic("Prochlorococcus")),ylab="")
title(ylab=expression(Exponential~Growth~Rate~of~italic("Synechococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(2,2,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

# Supp. ProExGR vs. transfer number.
tiff("Fig.ProExGR.Convergence.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$Transfer-0.1,subset(data.maggie,CO2=="0")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(1,7),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$Transfer+0.1,subset(data.maggie,CO2=="1")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(1,7),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$Transfer-0.1,subset(data.jeff,CO2=="0")$ProExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,1.5),xlim=c(1,7),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$Transfer+0.1,subset(data.jeff,CO2=="1")$ProExGR,ylim=c(-.5,1.5),xlim=c(1,7),pch=17,xlab="Transfer",ylab="")
title(ylab=expression(Exponential~Growth~Rate~of~italic("Prochlorococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(5,1.5,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

# Supp. SynExGR vs. transfer number.
tiff("Fig.SynExGR.Convergence.tiff",compression="lzw",units="in",width=5,height=5,res=600)
plot(subset(data.maggie,CO2=="0")$Transfer-0.1,subset(data.maggie,CO2=="0")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(1,7),pch=1)
par(new=TRUE)
plot(subset(data.maggie,CO2=="1")$Transfer+0.1,subset(data.maggie,CO2=="1")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(1,7),pch=2)
par(new=TRUE)
plot(subset(data.jeff,CO2=="0")$Transfer-0.1,subset(data.jeff,CO2=="0")$SynExGR,axes=FALSE,ann=FALSE,ylim=c(-.5,2),xlim=c(1,7),pch=16)
par(new=TRUE)
plot(subset(data.jeff,CO2=="1")$Transfer+0.1,subset(data.jeff,CO2=="1")$SynExGR,ylim=c(-.5,2),xlim=c(1,7),pch=17,xlab="Transfer",ylab="")
title(ylab=expression(Exponential~Growth~Rate~of~italic("Synechococcus")~(d^{-1})),mgp=c(2.3,1,0))
legend(5,2,c("400 ppm","800 ppm"),pch=c(1,2),cex=.8)
abline(h=0,lty=1)
dev.off()

#***************

#Master BIC table:

BIC.master<-BIC(model.out,model.out.noco2,model.total.dens,model.total.dens.noco2,model.pro.dens,model.pro.dens.noco2,model.syn.dens,model.syn.dens.noco2,model.full.dens,model.full.dens.noco2,model.out.ProGR,model.out.noco2.ProGR,model.total.dens.ProGR,model.total.dens.noco2.ProGR,model.pro.dens.ProGR,model.pro.dens.noco2.ProGR,model.syn.dens.ProGR,model.syn.dens.noco2.ProGR,model.full.dens.ProGR,model.full.dens.noco2.ProGR,model.out.SynGR,model.out.noco2.SynGR,model.total.dens.SynGR,model.total.dens.noco2.SynGR,model.pro.dens.SynGR,model.pro.dens.noco2.SynGR,model.syn.dens.SynGR,model.syn.dens.noco2.SynGR,model.full.dens.SynGR,model.full.dens.noco2.SynGR,model.exgr.out,model.exgr.out.noco2,model.exgr.total.dens,model.exgr.total.dens.noco2,model.exgr.pro.dens,model.exgr.pro.dens.noco2,model.exgr.syn.dens,model.exgr.syn.dens.noco2,model.exgr.full.dens,model.exgr.full.dens.noco2,model.exgr.out.ProExGR,model.exgr.out.noco2.ProExGR,model.exgr.total.dens.ProExGR,model.exgr.total.dens.noco2.ProExGR,model.exgr.pro.dens.ProExGR,model.exgr.pro.dens.noco2.ProExGR,model.exgr.syn.dens.ProExGR,model.exgr.syn.dens.noco2.ProExGR,model.exgr.full.dens.ProExGR,model.exgr.full.dens.noco2.ProExGR,model.exgr.out.SynExGR,model.exgr.out.noco2.SynExGR,model.exgr.total.dens.SynExGR,model.exgr.total.dens.noco2.SynExGR,model.exgr.pro.dens.SynExGR,model.exgr.pro.dens.noco2.SynExGR,model.exgr.syn.dens.SynExGR,model.exgr.syn.dens.noco2.SynExGR,model.exgr.full.dens.SynExGR,model.exgr.full.dens.noco2.SynExGR)

order.BIC<-order(BIC.master$BIC)
BIC.master<-BIC.master[order.BIC, ]

#Aggregate Syn and Pro growth rates, together and alone

together<-read.table("Together.csv",header=TRUE,sep=",")

table4<-ddply(together,~CO2*Together*Strain,summarise,meanGR=mean(RealizedGR),sdGR=sd(RealizedGR),meanXGR=mean(ExponentialGR),sdXGR=sd(ExponentialGR),meanLag=mean(Lag),sdLag=sd(Lag),length=length(RealizedGR))
write.csv(table4,"Table4.CoCultureGR.csv")

# Test difference of growth rate in single vs. co-culture. Outlier removed.

modelco<-lm(RealizedGR~CO2*Strain*Together,data=together)

summary(modelco)

emmeans(modelco,pairwise~CO2|(Together+Strain))
emmeans(modelco,pairwise~Together|(CO2+Strain))

#Significant effects of strain and togetherness. Pro grows faster together than separate under both CO2 conditions (-0.18, p=0.007 CO2-; -.35, p < 0.0001 CO2+), but there is no effect of togetherness on the Syn growth rate.  Moreover, there is a significant decrease in Pro growth rate at high CO2 alone, but not together. (-0.21, p = 0.02, vs 0.03, p = 0.5)

modelcoex<-lm(ExponentialGR~CO2*Strain*Together,data=together)
summary(modelcoex)
emmeans(modelcoex,pairwise~CO2|(Together+Strain))
emmeans(modelcoex,pairwise~Together|(CO2+Strain))

#Marginal effect of togetherness on Pro exponential growth rate under both circumstances (p = 0.06). 

modelcolag<-lm(Lag~CO2*Strain*Together,data=together)
summary(modelcolag)
emmeans(modelcolag,pairwise~CO2|(Together+Strain))
emmeans(modelcolag,pairwise~Together|(CO2+Strain))

#For Pro, significantly greater lag at high CO2 alone, but not when together.  Significant reduction in lag time together under BOTH conditions, but much greater at high CO2.  No effect for Syn, and in fact lags were not significantly different for 0 for any strains except for Pro alone.

# Plot Growth rate alone vs. together for each genus
tiff("Fig.TogetherPro.tiff",compression="lzw",units="in",width=5,height=5,res=600)
labels<-c("Alone","Mixed")
boxplot(RealizedGR~CO2*Together,data=subset(together,Strain=="Prochlorococcus"),ylab="",xaxt="n",ylim=c(-0.3,2),col=c("white","gray","white","gray"))
title(ylab=expression(Realized~Growth~Rate~(d^{-1})),mgp=c(2.3,1,0))
axis(side=1,at=c(1.5,3.5),tick=TRUE,labels=labels)
legend(3.5,2,c("400 ppm","800 ppm"),fill=c("white","gray"),cex=.8)
abline(h=0,lty=2)
segments(c(1,2,1),c(1.6,1.4,1.2),c(3,4,2),c(1.6,1.4,1.2))
text(c(2,3,1.5),c(1.7,1.5,1.3),c("**","***","*"),cex=2)
dev.off()

tiff("Fig.TogetherSyn.tiff",compression="lzw",units="in",width=5,height=5,res=600)
labels<-c("Alone","Mixed")
boxplot(RealizedGR~CO2*Together,data=subset(together,Strain=="Synechococcus"),ylab="",xaxt="n",ylim=c(-0.3,2),col=c("white","gray","white","gray"))
title(ylab=expression(Realized~Growth~Rate~(d^{-1})),mgp=c(2.3,1,0))
axis(side=1,at=c(1.5,3.5),tick=TRUE,labels=labels)
legend(.5,2,c("400 ppm","800 ppm"),fill=c("white","gray"),cex=.8)
abline(h=0,lty=2)
dev.off()

tiff("Fig.TogetherProEx.tiff",compression="lzw",units="in",width=5,height=5,res=600)
labels<-c("Alone","Mixed")
boxplot(ExponentialGR~CO2*Together,data=subset(together,Strain=="Prochlorococcus"),ylab="",xaxt="n",ylim=c(-0.3,2),col=c("white","gray","white","gray"))
title(ylab=expression(Exponential~Growth~Rate~(d^{-1})),mgp=c(2.3,1,0))
axis(side=1,at=c(1.5,3.5),tick=TRUE,labels=labels)
legend(3.5,2,c("400 ppm","800 ppm"),fill=c("white","gray"),cex=.8)
abline(h=0,lty=2)
segments(c(1,2),c(1.6,1.4),c(3,4),c(1.6,1.4))
text(c(2,3),c(1.7,1.5),c("p = 0.05","p = 0.07"))
dev.off()

tiff("Fig.TogetherSynEx.tiff",compression="lzw",units="in",width=5,height=5,res=600)
labels<-c("Alone","Mixed")
boxplot(ExponentialGR~CO2*Together,data=subset(together,Strain=="Synechococcus"),ylab="",xaxt="n",ylim=c(-0.3,2),col=c("white","gray","white","gray"))
title(ylab=expression(Exponential~Growth~Rate~(d^{-1})),mgp=c(2.3,1,0))
axis(side=1,at=c(1.5,3.5),tick=TRUE,labels=labels)
legend(.5,2,c("400 ppm","800 ppm"),fill=c("white","gray"),cex=.8)
abline(h=0,lty=2)
segments(3,1.8,4,1.8)
text(3.5,1.9,"*",cex=2)
dev.off()

tiff("Fig.TogetherProLag.tiff",compression="lzw",units="in",width=5,height=5,res=600)
labels<-c("Alone","Mixed")
boxplot(Lag~CO2*Together,data=subset(together,Strain=="Prochlorococcus"),ylab="",xaxt="n",ylim=c(-2,20),col=c("white","gray","white","gray"))
title(ylab="Lag (d)",mgp=c(2.3,1,0))
axis(side=1,at=c(1.5,3.5),tick=TRUE,labels=labels)
legend(3.5,20,c("400 ppm","800 ppm"),fill=c("white","gray"),cex=.8)
abline(h=0,lty=2)
segments(c(1,2,1),c(17.5,16,19),c(3,4,2),c(17.5,16,19))
text(c(2,3,1.5),c(18,16.5,19.5),c("***","***","***"),cex=2)
dev.off()

tiff("Fig.TogetherSynLag.tiff",compression="lzw",units="in",width=5,height=5,res=600)
labels<-c("Alone","Mixed")
boxplot(Lag~CO2*Together,data=subset(together,Strain=="Synechococcus"),ylab="",xaxt="n",ylim=c(-2,20),col=c("white","gray","white","gray"))
title(ylab="Lag (d)",mgp=c(2.3,1,0))
axis(side=1,at=c(1.5,3.5),tick=TRUE,labels=labels)
legend(.5,20,c("400 ppm","800 ppm"),fill=c("white","gray"),cex=.8)
abline(h=0,lty=2)
dev.off()