library(deSolve) system("R CMD SHLIB sir_fig2.c") dyn.load("sir_fig2.dll") parms=c(p=0.5,d=0.75,m=73,K=25000,r_val=seq(0.01,19.91,le=200),beta_val=seq(0.01,19.91,le=200),df=17.4,rf=175,Kf=2,a=0.004) times=seq(0,600,by=0.1) Y=c(S=500,I=10,R =490,N=3,F=1000) r_val=seq(0.01,19.91,le=200) beta_val=seq(0.01,19.91,le=200) out=NULL val <- expand.grid(r_val,beta_val) out <- sapply(1:nrow(val), function(x){ out = ode(y=Y,times=times,func="derivs", parms=c(p=0.5,d=0.75,m=73,K=25000,beta=val[x,2],r=val[x,1],df=17.4,rf=175,Kf=2, a=0.004),dllname="sir_fig2",initfunc="initmod",nout=1,outnames="Sum") out = tail(out,200) minN = c(min(out[,2]),min(out[,3]),min(out[,4])) return(minN) }) out = t(out) write.table(out,"figure2.txt") # FIGURE r_val=seq(0.01,19.91,le=200) beta_val=seq(0.01,19.91,le=200) out3=out>1 out3[which(is.na(out3))]="FALSE" z=NULL for (i in 1:nrow(out3)){ if (out3[i,1]=="FALSE" & out3[i,2]=="FALSE" & out3[i,3]=="FALSE") {z[i]=0} if (out3[i,1]=="TRUE" & out3[i,2]=="FALSE" & out3[i,3]=="FALSE") {z[i]=1} if (out3[i,1]=="TRUE" & out3[i,2]=="TRUE" & out3[i,3]=="FALSE") {z[i]=2} if (out3[i,1]=="TRUE" & out3[i,2]=="TRUE" & out3[i,3]=="TRUE") {z[i]=3} if (out3[i,1]=="TRUE" & out3[i,2]=="FALSE" & out3[i,3]=="TRUE") {z[i]=1} } color=c("black",grey(0.4),grey(0.8),"white") pdf("figure2.pdf",width=3.1,height=3.5, pointsize=9, family="Times") par(mgp=c(1.7,0.5,0)) image(r_val,beta_val,matrix(z,nrow=length(r_val)), xlim=c(0.01,19.91), ylim=c(0.01,19.91), xlab="Rat intrinsic birth rate, r", ylab=expression(paste("Transmission rate, ",beta)), col=color) y=NULL w=NULL for (i in 1:200) {y=rbind(y,(17.4+1-exp(-0.004*25000))/(2*(1-exp(-0.004*25000))))} points(r_val,y,t="l",cex=1.5) text(18.5,8.6,labels="Ro=1") for (i in 1:200) {w=rbind(w,0.75)} points(w,beta_val,t="l",cex=1.5) text(1.7,0.6,labels="r=d") legend("topright",legend=c("0","S","S,I","S,I,R"),fill=color,bg="white") dev.off()