## Introduction

In this exercise I demonstrate the use of Gillespie’s direct method for compartmental epidemic models.

My approach will be to wrap a function SIR.onestep (that generates both the next event time and event type) within another function that calls SIR.onestep iteratively until stopping criterion is reached.

Here is illustration for $$SIR$$ model without demography.

SIR.onestep <- function(x, params){
X <- x[2]
Y <- x[3]
Z <- x[4]
N <- X+Y+Z
with(as.list(params),{
R <- c(beta*X*Y/N, gamma*Y)
E <- matrix(c(-1,1,0,
0,-1,1),
ncol=3, byrow=TRUE)
S <- -log(runif(1))/sum(R)   # Step 1 of Gillespie
U <- runif(1)
rescaledU <- U*sum(R)
m <- min(which(cumsum(R) >= rescaledU))
x <- x[2:4] + E[m,]
return(out <- c(S, x))
}
)
}

Here we write SIR.model to iteratively call SIR.onestep to simulate one epidemic. Note that because there is no demography, I had to put a break in to terminate the loop if the epidemic dies out. I did this with the code if (x[3] == 0) break.

SIR.model <- function(x, params, nstep=1000){
output <- array(dim=c(nstep+1, 4))
colnames(output) <- c('Time', 'Susceptibles', 'Infected', 'Recovered')
output[1,] <- x
for(i in 1:nstep){
output[i+1,] <- x <- SIR.onestep(x, params)
if (x[3] == 0) break
}
return(output)
}

Now, I can set up some initial conditions and simulate. At the end of this code chunk I remove any rows that are NA and convert to a data frame.

xstart <- c(time=0, X=100, Y=15, Z=0)
params <- list(beta=60, gamma=365/13)
result <- SIR.model(xstart, params, nstep=1000)
result <- as.data.frame(na.omit(result))

The variable Time is just the time increments. I calculate the cumulative sum of this to get chronological time.

t <- c(cumsum(result$Time)) Now I can plot the number of infected over time. plot(t, result$Infected, col='red', lwd=2, type='l',
xlab='Time',
ylab='Infected')