model <- function(t, state, parms) { with(as.list(c(state,parms)), { f <- v*R/(R+k) tlag <- t - lambda if (tlag < 0) lags <- rep(0,5) # no initial infection else lags <- lagvalue(tlag) # returns 5 lags dtR <- - f*e*(B0+B1) dtB0 <- f*B0 - delta*B0*P dtB1 <- f*B1 dtM <- delta*B0*P - delta*lags[2]*lags[5] dtP <- b*delta*lags[2]*lags[5] - 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,m=0) s <- c(R=350,B0=5e6,B1=500,M=0,P=5e6) tweak <- "nsol$B<-nsol[,\"B0\"]+nsol[,\"B1\"]+nsol[,\"M\"]" run(7,0.1,ymin=1e3,ymax=1e10,log="y",delay=T,tweak=tweak) 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,delay=T) p[free] <- f0$par s <- c(R=350,B0=4382322,B1=1000,M=0,P=5.58e6) free <- c("B1","b","delta","lambda") f1 <- fit(fig2B,tstep=0.1,ymin=1e3,ymax=1e10,log="y",fun=log1p,free=free,tweak=tweak,lower=0,delay=T) p[free[2:4]] <- f1$par[2:4] # 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-7*tstep times <- seq(0,7,by=tstep) run(times=times,ymin=1,ymax=1e10,log="y",delay=T,tweak=tweak,events=list(func=fun,time=times))