In ICI3D/ici3d-pkg: Materials for the International Clinics on Infectious Disease Dynamics and Data (ICI3D) Workshops

This tutorial simulates a simple population model where mortality is the only process. The differential equation version of the model is:

$$ {dN \over dt} = - \mu N $$

where $N$ is the population size as a function of time $t$ and $\mu$ is the mortality rate.

This R code defines a function that calculates the population over a time range, using the analytical solution to the differential equation:

analytic <- function(
  initial = initialPopulationSize,
  rate = mortalityRate,
  stepSize = timeStep,
  maxT = maxTime
) {
  times = seq(0, maxT, stepSize)
  ts <- data.frame(
    time = times,
    populationSize = initial * exp(-rate*times)
  )
  return(ts)
}

This R code defines a function that calculates the population over a time range, using a discrete time approximation to the differential equation:

discreteTime <- function(
  initial = initialPopulationSize,
  rate = mortalityRate,
  stepSize = timeStep,
  maxT = maxTime
) {
  tt <- 0
  pop <- initial
  ts <- data.frame(time = tt, populationSize = pop)
  while (tt<maxT){
    nextt <- tt + stepSize
    nextpop <- pop - rate*pop*stepSize
    ts <- rbind(
      ts, c(time = nextt, populationSize = nextpop)
    )
    tt <- nextt
    pop <- nextpop
  }
  return(ts)
}

This R code defines a function that calculates the population over a time range, using a stochastic simulation of the process represented by the differential equation model:

individual <- function(
  initial = initialPopulationSize,
  rate = mortalityRate,
  stepSize = timeStep,
  maxT = maxTime
) {
  deathTimes <- rexp(initial, rate)
  ts <- data.frame(time = 0, populationSize = initial)
  for(tt in seq(stepSize, maxT, stepSize)){
    ts <- rbind(
      ts,
      c(time = tt, populationSize = sum(deathTimes >= tt)))
  }
  return(ts)
}

Set the parameter values

initialPopulationSize <- 10 # number of individuals
mortalityRate <- 0.05       # per capital deaths per day
timeStep <- 1                               # days
maxTime <- 30                               # days

Setup plotting

par(bty='L',lwd=3,mar=c(4,4,1,1))
plot(NA,NA,ylim=c(0,initialPopulationSize),xlim=c(0,maxTime),
     ylab='Population size',xlab='Time')
plotter <- function(ts, col) with(ts, lines(time, populationSize, col=col))

Now run the three functions to see what you get

plotter(analytic(), col="black")
plotter(discreteTime(), col="green")
plotter(individual(), col="red")

Try running these three lines multiple times without resetting the parameters or re-making the plot. Which functions give you different outcomes each time?

Now try changing the parameter values above and re-running the functions (you may want to re-make the plot as well). How does changing each of the values change the output? Can you get the greencurve to diverge from the black curve? How?

Stochastic SIR simulation

Compartments: (S,I,R) = (susceptible, infectious, removed)

Transitions:

| Event | Change | Rate | |:---------------------:|:-------------------------:|:-----------:| | Spillover (S) | (S, I, R)->(S-1, I+1, R) | $\lambda$ | | Infection (S) | (S, I, R)->(S-1, I+1, R) | $\beta I S$ | | Recovery/Removal (I) | (S, I, R)->(S, I-1, R+1) | $\gamma I$ |

transitions <- list(
  spillover = c(S = -1, I =  1, R = 0),
  infection = c(S = -1, I =  1, R = 0),
  recovery  = c(S =  0, I = -1, R = 1),
  end = c(S = 0, I = 0, R = 0)
)

N=50   # population size

parms=list(
  lambda=.02,      # spillover rate
  beta=.2,           # contact rate
  gamma=.1        # recovery rate
)

event <- function(prevevent, params){
  with(c(prevevent, params),{
    # update rates
    rates <- c(
      spillover = ifelse(S>0, lambda, 0), # no spillover infections if S depleted
      infect = beta*I*S/N,
      recover = gamma*I
    )

    totRate <- sum(rates)

    # if the event rate has gone to 0, skip to end of simulation
    if (totRate==0) {
      eventTime <- final.time
      eventType <- "end"
    } else {
      # calculate time of the next event
      eventTime <- time + rexp(1, totRate)

      # choose type of event
      eventType <- sample(
        c("spillover","infection","recovery"),
        size=1, replace=F, prob=rates/totRate
      )
    }

    # determine changes based on the event type
    update <- transitions[[eventType]]

    return(data.frame(
      time = eventTime,
      as.list(update + c(S, I, R)),
      outcome = eventType
    ))
  })
}

simulateSIR <- function(t,y,params){
  with(as.list(y),{
    ts <- data.frame(
      time=0,
      S=round(S), I=round(I), R=round(R),
      outcome = "spillover"
    )
    nextEvent <- ts
    while(nextEvent$time < final.time){
      nextEvent <- event(nextEvent,params)
      ts <- rbind(ts, nextEvent)
    }
    return(ts)
  })
}

final.time=400
parms=parms

# run the simulation
tsTest <- simulateSIR(final.time, c(S=N-1,I=1,R=0), parms)

# plot the infectious individuals over time
plot(
  tsTest$time, tsTest$I,
  type='s', main="Number Infected",
  ylim=c(0,N), xlim=c(0,final.time), xlab="Time", ylab="",
  bty="n", cex.main=2, cex.lab=1.5, cex.axis=1.25, lwd=2
)

# plot all compartments over time
cols <- c(S="black",I="red2",R="purple")
plot(
  tsTest$time, tsTest$S,
  type='s', ylim=c(0,N+1), bty="n",
  ylab="Number of individuals", xlab="time", lwd=2,
  col = cols["S"]
)
lines(tsTest$time, tsTest$I, type='s', col = cols["I"], lwd = 2)
lines(tsTest$time, tsTest$R, type='s', col = cols["R"], lwd = 2)
legend(
  x = "right", names(cols),
  lty = 1, col = cols,
  bty = "n", lwd = 2
)

# print the total number of infections
sum(tsTest$outcome == "infection")

# print the number of spillover events
sum(tsTest$outcome == "spillover") - 1

ICI3D/ici3d-pkg documentation built on July 2, 2023, 1:59 p.m.