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')
```