##Make figures for revised Evolution paper setwd("C:/Users/Erik Osnas/Documents/House Finch/Virulence Evol MS/Two Spp Model/Manusc Figs/Fig 1") #Figure 1 dat = read.table(file="contraint.figure.dat",header=TRUE) plot(1/dat$parvalues,dat$ESP.1,col="black",ylim=c(0,20),xlim=c(0,0.12), ylab="Evolutionarily Stable Virulence ",xlab="Constraint (1/c)", type="l",lwd=2) lines(1/dat$parvalues,dat$ESP.2,col=gray(0.5),lwd=2) lines(1/dat$parvalues,dat$ESP.3,col="black",lwd=2) segments(x0=1/30,y0=2.5,x1=1/30,y1=10,lty=2) #Alt. Fig 1 plot(dat$parvalues,dat$ESP.1,col="black",ylim=c(0,20),xlim=c(0,100), ylab="Evolutionarily Stable Virulence ",xlab="Constraint (c)", type="l",lwd=2) lines(dat$parvalues,dat$ESP.2,col=gray(0.5),lwd=2) lines(dat$parvalues,dat$ESP.3,col="black",lwd=2) segments(x0=30,y0=2.5,x1=30,y1=10,lty=2) x2=matrix(NA,nrow=1000,ncol=10) c=c(1:10)*10 x1=seq(0.001,100,length=1000) for(i in 1:10){x2[,i]=c[i]/x1} matplot(x=x1,y=x2,type="l",lty=1,xlim=c(0,20),ylim=c(0,20),col=c(1,2,4)) points(dat$ESP.1,dat$parvalues/dat$ESP.1,pch=20) points(dat$ESP.2,dat$parvalues/dat$ESP.2,pch=20,col=gray(0.5)) points(dat$ESP.3,dat$parvalues/dat$ESP.3,pch=20) points(dat$parvalues/dat$ESP.1,dat$parvalues,pch=20) points(dat$parvalues/dat$ESP.2,dat$parvalues,pch=20,col=gray(0.5)) points(dat$parvalues/dat$ESP.3,dat$parvalues,pch=20) plot(dat$ESP.1,dat$parvalues/dat$ESP.1,pch=20,xlim=c(0,10),ylim=c(0,30)) points(dat$ESP.2,dat$parvalues/dat$ESP.2,pch=20,col=gray(0.5)) points(dat$ESP.3,dat$parvalues/dat$ESP.3,pch=20) plot(dat$parvalues/dat$ESP.1,dat$ESP.1,pch=20,xlim=c(0,10),ylim=c(0,30)) points(dat$parvalues/dat$ESP.2,dat$ESP.2,pch=20,col=gray(0.5)) points(dat$parvalues/dat$ESP.3,dat$ESP.3,pch=20) #Fgiure 2 #Original ESdata = read.table(file="figure2.sij is 0.1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=16, ylab="Evolutionarily Stable Virulence (Species 1)", xlab=expression(paste("Background Mortality ","(",mu,")"," in Species 1"))) points(ESdata$parvalues,ESdata$ESP.2,col=gray(0.5),pch=16) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=16) ESdata=read.table(file="figure2.sij is 1.dat",header=TRUE) points(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=16,cex=0.1) points(ESdata$parvalues,ESdata$ESP.2,col=gray(0.5),pch=16,cex=0.1) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=16,cex=0.1) ESdata=read.table(file="figure2.sij is 0.5.dat",header=TRUE) points(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=16,cex=0.5) points(ESdata$parvalues,ESdata$ESP.2,col=gray(0.5),pch=16,cex=0.5) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=16,cex=0.5) #Revised par(mfrow=c(3,1),oma=c(5,15,5,15),mar=c(0.5,0.5,0.5,0.5)) ESdata = read.table(file="figure2.sij is 0.1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),type="l",lwd=3, ann=FALSE,xaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,lwd=1,col="black") lines(ESdata$parvalues,ESdata$ESP.3,lwd=3,col="black") mtext(expression(paste(S[ij]," = 0.1")), side=2, font=2, line=5) ESdata = read.table(file="figure2.sij is 0.5.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),type="l",lwd=3, ann=FALSE,xaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,lwd=1,col="black") lines(ESdata$parvalues,ESdata$ESP.3,lwd=3,col="black") mtext(expression(paste(S[ij]," = 0.5")), side=2, font=2, line=5) ESdata = read.table(file="figure2.sij is 1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,cex=1,ann=FALSE) lines(ESdata$parvalues,ESdata$ESP.2,lwd=1,col="black") lines(ESdata$parvalues,ESdata$ESP.3,lwd=3,col="black") mtext(expression(paste(S[ij]," = 1")), side=2, font=2, line=5) mtext(expression(paste("Background Mortality ","(",mu,")"," in Species 1")),side=1, line=3) par(mfrow=c(1,1)) mtext("Evolutionarily Stable Virulence (Species 1)", side=2, line=-6.25) #Figure 3, used to be 3a par(mfrow=c(3,1),oma=c(5,15,5,15),mar=c(0.5,0.5,0.5,0.5)) ESdata=read.table(file="figure3a.sij is 0.1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,xaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,col="black") #,pch=16,cex=0.1) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"a",font=2) mtext(expression(paste(S[ij]," = 0.1")), side=2, font=2, line=5) ESdata=read.table(file="figure3a.sij is 0.5.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,xaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,col="black") points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"b",font=2) mtext(expression(paste(S[ij]," = 0.5")), side=2, font=2, line=5) ESdata=read.table(file="figure3a.sij is 1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20) lines(ESdata$parvalues,ESdata$ESP.2,col="black", lwd=1) #,pch=20,cex=0.5) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"c",font=2) mtext(expression(paste(S[ij]," = 1")), side=2, font=2, line=5) #Fill in above! ESdata=read.table(file="figure3a.sij is 1.fill.dat",header=TRUE) points(ESdata$parvalues,ESdata$ESP.1,col="black",pch=20) lines(ESdata$parvalues,ESdata$ESP.2,col="black", lwd=1) #cex=0.5,pch=20) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) mtext(expression(paste("Relative Mortality (",rho,") in Species 2")),side=1, line=3) par(mfrow=c(1,1)) mtext("Evolutionarily Stable Virulence (Species 1)", side=2, line=-6.25) #Figure 4, used to be 3b par(mfrow=c(3,1),oma=c(5,15,5,15),mar=c(0.5,0.5,0.5,0.5)) ESdata=read.table(file="figure3b.sij is 0.1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,xaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,col="black") #,pch=16,cex=0.1) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"a",font=2) mtext(expression(paste(S[ij]," = 0.1")), side=2, font=2, line=5) ESdata=read.table(file="figure3b.sij is 0.5.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,xaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,col="black") points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"b",font=2) mtext(expression(paste(S[ij]," = 0.5")), side=2, font=2, line=5) ESdata=read.table(file="figure3b.sij is 1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20) lines(ESdata$parvalues,ESdata$ESP.2,col="black", lwd=1) #,pch=20,cex=0.5) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"c",font=2) mtext(expression(paste(S[ij]," = 1")), side=2, font=2, line=5) #Fill in above! ESdata=read.table(file="figure3b.sij is 1.fill.dat",header=TRUE) points(ESdata$parvalues,ESdata$ESP.1,col="black",pch=20) lines(ESdata$parvalues,ESdata$ESP.2,col="black", lwd=1) #cex=0.5,pch=20) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) mtext(expression(paste("Relative Transmission (",phi,") in Species 2")),side=1, line=3) par(mfrow=c(1,1)) mtext("Evolutionarily Stable Virulence (Species 1)", side=2, line=-6.25) ##Figure 3 alternate 6 pannel figure par(mfrow=c(3,2),oma=c(6,7.5,4,7.5),mar=rep(0,4),pty="s") ESdata=read.table(file="figure3a.sij is 0.1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,xaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,col="black") #,pch=16,cex=0.1) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"a",font=2) mtext(expression(paste(S[ij]," = 0.1")), side=2, font=2, line=5) ESdata=read.table(file="figure3b.sij is 0.1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,xaxt="n",yaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,col="black") #,pch=16,cex=0.1) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"d",font=2) ESdata=read.table(file="figure3a.sij is 0.5.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,xaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,col="black") points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"b",font=2) mtext(expression(paste(S[ij]," = 0.5")), side=2, font=2, line=5) ESdata=read.table(file="figure3b.sij is 0.5.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,xaxt="n",yaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,col="black") points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"e",font=2) ESdata=read.table(file="figure3a.sij is 1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20) lines(ESdata$parvalues,ESdata$ESP.2,col="black", lwd=1) #,pch=20,cex=0.5) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"c",font=2) mtext(expression(paste(S[ij]," = 1")), side=2, font=2, line=5) mtext(expression(paste("Relative Mortality (",rho,")")),side=1, line=3) #Fill in above! ESdata=read.table(file="figure3a.sij is 1.fill.dat",header=TRUE) points(ESdata$parvalues,ESdata$ESP.1,col="black",pch=20) lines(ESdata$parvalues,ESdata$ESP.2,col="black", lwd=1) #cex=0.5,pch=20) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) ESdata=read.table(file="figure3b.sij is 1.dat",header=TRUE) plot(ESdata$parvalues,ESdata$ESP.1,col="black",ylim=c(0,25),pch=20,yaxt="n") lines(ESdata$parvalues,ESdata$ESP.2,col="black", lwd=1) #,pch=20,cex=0.5) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) text(0,25,"f",font=2) #Fill in above! ESdata=read.table(file="figure3b.sij is 1.fill.dat",header=TRUE) points(ESdata$parvalues,ESdata$ESP.1,col="black",pch=20) lines(ESdata$parvalues,ESdata$ESP.2,col="black", lwd=1) #cex=0.5,pch=20) points(ESdata$parvalues,ESdata$ESP.3,col="black",pch=20) mtext(expression(paste("Relative Transmission (",phi,")")),side=1, line=3) par(mfrow=c(1,1)) mtext("Evolutionarily Stable Virulence in Species 1", side=2, line=3.5) mtext("in Species 2", side=1, line=7) ##Figure 4 source("C:/Users/Erik Osnas/Documents/House Finch/Virulence Evol MS/Two Spp Model/AD.sim.function.r") source("C:/Users/Erik Osnas/Documents/House Finch/Virulence Evol MS/Two Spp Model/two.species.ess.r") par.1=list( psi1=0.1,psi2=0.1,constraint=30, B1=NA,B2=NA,s12=1,s21=1,phi1=1,phi2=1, d1=1,d2=1,rho1=1,rho2=1, alpha1=NA,alpha2=NA, sig1=1,sig2=1,gamma1=50,gamma2=50) AD.sim.1=ADsim(par=par.1,sdmut=0.1) matplot(x=AD.sim.1$num,y=AD.sim.1$vir.1,type="p",col=1,pch=".") matpoints(x=AD.sim.1$num,AD.sim.1$vir.2,col=1,pch=".",) AD.sim.2=ADsim(par=par.1,sdmut=1) matplot(x=AD.sim.2$num,y=AD.sim.2$vir.1,type="p",col=1,pch=".") matpoints(x=AD.sim.2$num,AD.sim.2$vir.2,col=1,pch=".",) par.2=list( psi1=0.1,psi2=0.1,constraint=30, B1=NA,B2=NA,s12=0.5,s21=0.5,phi1=1,phi2=1, d1=1,d2=1,rho1=1,rho2=1, alpha1=NA,alpha2=NA, sig1=1,sig2=1,gamma1=50,gamma2=50) AD.sim.3=ADsim(par=par.2,sdmut=0.1) matplot(x=AD.sim.3$num,y=AD.sim.3$vir.1,type="p",col=1,pch=".") matpoints(x=AD.sim.3$num,AD.sim.3$vir.2,col=1,pch=".") AD.sim.4=ADsim(par=par.2,sdmut=1) matplot(x=AD.sim.4$num,y=AD.sim.4$vir.1,type="p",col=1,pch=".") matpoints(x=AD.sim.4$num,y=AD.sim.4$vir.2,col=1,pch=".") setwd("C:/Users/Erik Osnas/Documents/House Finch/Virulence Evol MS/Two Spp Model/Manusc Figs/Fig 1") write.table(AD.sim.1,file="AD.sim.sdmut 0.1 sij 1.dat") write.table(AD.sim.2,file="AD.sim.sdmut 1 sij 1.dat") write.table(AD.sim.3,file="AD.sim.sdmut 0.1 sij 0.5.dat") write.table(AD.sim.4,file="AD.sim.sdmut 1 sij 0.5.dat") write.table(PIP,file="PIP.sij 0.5.dat") write.table(PIP,file="PIP.sij 1.dat") ##Define function to make PIP plots plotPIP=function(PIP=NA,a=NA){ image(a,a,PIP,col=gray(c(1,0.5)),breaks=c(-1000,0,1000),ann=FALSE, xlab="Resident Virulence",ylab="Mutant Virulence") ##Try to find region of ploymorphism size=length(a) coexit=matrix(NA,nrow=size,ncol=size) for(i in 1:size){ for(j in 1:size){ if(is.na(PIP[i,j])){ # points(a[i],a[j],pch=19,col="blue",cex=0.1) } else{if(PIP[i,j]>0){ if(PIP[j,i]>0){ if(i!=j){ points(a[i],a[j],pch=19,col=gray(0.7),cex=0.1) points(a[j],a[i],pch=19,col=gray(0.7),cex=0.1) } } } } #if(is.na(PIP[j,i])){points(a[i],a[j],pch=19,col="blue",cex=0.1)} } } } ############################# ##Read data and make figures AD.sim.1=list(num=numeric(100),vir.1=matrix(NA,100,100),vir.2=matrix(NA,100,100)) temp=read.table(file="AD.sim.sdmut 0.1 sij 1.dat") AD.sim.1$num=temp$num AD.sim.1$vir.1=as.matrix(temp[,2:101]) AD.sim.1$vir.2=as.matrix(temp[,102:201]) AD.sim.2=list(num=numeric(100),vir.1=matrix(NA,100,100),vir.2=matrix(NA,100,100)) temp=read.table(file="AD.sim.sdmut 1 sij 1.dat") AD.sim.2$num=temp$num AD.sim.2$vir.1=as.matrix(temp[,2:101]) AD.sim.2$vir.2=as.matrix(temp[,102:201]) AD.sim.3=list(num=numeric(100),vir.1=matrix(NA,100,100),vir.2=matrix(NA,100,100)) temp=read.table(file="AD.sim.sdmut 0.1 sij 0.5.dat") AD.sim.3$num=temp$num AD.sim.3$vir.1=as.matrix(temp[,2:101]) AD.sim.3$vir.2=as.matrix(temp[,102:201]) AD.sim.4=list(num=numeric(100),vir.1=matrix(NA,100,100),vir.2=matrix(NA,100,100)) temp=read.table(file="AD.sim.sdmut 1 sij 0.5.dat") AD.sim.4$num=temp$num AD.sim.4$vir.1=as.matrix(temp[,2:101]) AD.sim.4$vir.2=as.matrix(temp[,102:201]) par(mfrow=c(2,3),oma=c(5,5,5,0),mar=rep(0,4),pty="s") PIP.1=as.matrix(read.table(file="PIP.sij 0.5.dat")) PIP.2=as.matrix(read.table(file="PIP.sij 1.dat")) a=seq(0.001,20,length=200) par(mar=c(0,0,0,2)) plotPIP(PIP.1,a) box() mtext("Mutant Virulence",side=2, line=2) mtext("Resident Virulence",side=1, line=2.5) mtext(expression(paste(S[ij]," = 0.5")),side=2, line=3.5) mtext("a",font=2,adj=0,line=0.5) mtext("PIP",side=3, line=2) par(mar=c(0,2,0,0)) matplot(x=AD.sim.3$num,y=AD.sim.3$vir.1,type="p",col=1,pch=16,cex=0.1,ylim=c(0,20),yaxt="n",ann=FALSE) mtext("Time",side=1, line=2.5) mtext("Virulence",side=2, line=1) mtext("b",font=2,adj=0,line=0.5) mtext("v = 0.1",side=3, line=2) par(mar=c(0,1,0,1)) matplot(x=AD.sim.4$num,y=AD.sim.4$vir.1,type="p",col=1,pch=16,cex=0.1,ylim=c(0,20),yaxt="n",ann=FALSE) mtext("Time",side=1, line=2.5) mtext("c",font=2,adj=0,line=0.5) mtext("v = 1.0",side=3, line=2) a=seq(0.001,10,length=200) par(mar=c(0,0,0,2)) plotPIP(PIP.2,a) box() mtext("Mutant Virulence",side=2, line=2) mtext("Resident Virulence",side=1, line=2.5) mtext(expression(paste(S[ij]," = 1.0")),side=2, line=3.5) mtext("d",font=2,adj=0,line=0.5) par(mar=c(0,2,0,0)) matplot(x=AD.sim.1$num,y=AD.sim.1$vir.1,type="p",col=1,pch=16,cex=0.1,yaxt="n",ann=FALSE) matpoints(x=AD.sim.1$num,AD.sim.1$vir.2,col=1,pch=16, cex=0.1) mtext("Time",side=1, line=2.5) mtext("Virulence",side=2, line=1) mtext("e",font=2,adj=0,line=0.5) par(mar=c(0,1,0,1)) matplot(x=AD.sim.2$num,y=AD.sim.2$vir.1,type="p",col=1,pch=16,cex=0.1,yaxt="n",ann=FALSE) matpoints(x=AD.sim.2$num,AD.sim.2$vir.2,col=1,pch=16,cex=0.1) points(x=AD.sim.2$num,y=AD.sim.2$vir.1[,60], pch=16, cex=0.5) points(x=AD.sim.2$num,y=AD.sim.2$vir.2[,60], pch=16, cex=0.5) mtext("Time",side=1, line=2.5) mtext("f",font=2,adj=0,line=0.5) ##Illustrate single runs for(i in 1:length(AD.sim.2$num)){ matplot(x=AD.sim.2$num,y=AD.sim.2$vir.1,type="p",col=1,pch=".",yaxt="n",ann=FALSE) matpoints(x=AD.sim.2$num,AD.sim.2$vir.2,col=1,pch=".") points(x=AD.sim.2$num,y=AD.sim.2$vir.1[,i]) points(x=AD.sim.2$num,y=AD.sim.2$vir.2[,i]) text(x=80,y=9, labels=i) Sys.sleep(0.5) }