model <- function(t, state, parms) { with(as.list(c(state,parms)), { f <- v*R/(R+k) dtR <- - f*e*(B0+(1-s)*B1) dtB0 <- f*B0 - delta*B0*P dtB1 <- (1-s)*f*B1 dtM <- delta*B0*P - M/lambda dtP <- b*M - delta*P*(B0+B1) return(list(c(dtR, dtB0, dtB1, dtM, dtP))) }) } p <- c(b=80,delta=5e-8,e=5e-7,lambda=0.4,v=1.4,k=1,s=0,m=0) s <- c(R=350,B0=5e6,B1=500,M=0,P=5e6) after <- "state<-ifelse(state<0,0,state)" tweak <- "nsol$B<-nsol[,\"B0\"]+nsol[,\"B1\"]+nsol[,\"M\"]" run(7,0.1,ymin=1e3,ymax=1e10,log="y",tweak=tweak,after=after) fig2B0 <- read.csv("Levin_pg13_fig2_B0.csv",header=TRUE) fig2B <- read.csv("Levin_pg13_fig2.csv",header=TRUE) timePlot(fig2B0,ymin=1e3,ymax=1e10,log="y") timePlot(fig2B,ymin=1e3,ymax=1e10,log="y") free <- c("v") s <- c(R=350,B0=5189702,B1=0,M=0,P=0) f0 <- fit(fig2B0,tstep=0.1,ymin=1e3,ymax=1e10,log="y",fun=log1p,lower=0,free=free,tweak=tweak) p[free] <- f0$par p["s"] <- 0.001 s <- c(R=350,B0=4382322,B1=1000,M=0,P=5.58e6) free <- c("B1","b","delta","lambda","s") f1 <- fit(fig2B,tstep=0.1,ymin=1e3,ymax=1e10,log="y",fun=log1p,free=free,tweak=tweak,after=after,lower=0) p[free[2:5]] <- f1$par[2:5] # now with stochastic generation of B1: fun <- function(t, y, parms){ y <- ifelse(y<0, 0, y) # Set negatives to zero x <- parms["m"]*y["M"] # Expected number of BIMs if (x > 1) { y["M"] <- y["M"] - x; y["B1"] <- y["B1"] + x } else if (runif(1,0,1) < x) { # Create one BIM y["M"] <- y["M"] - 1; y["B1"] <- y["B1"] + 1 } return(y) } s["B1"] <- 0 tstep <- 0.01 p["m"] <- 1e-6*tstep times <- seq(0,7,by=tstep) run(times=times,ymin=1,ymax=1e10,log="y",tweak=tweak,events=list(func=fun,time=times))