Stochastic epidemic models

6 April 2022

Julien Arino

Department of Mathematics & Data Science Nexus
University of Manitoba*

Canadian Centre for Disease Modelling
Canadian COVID-19 Mathematical Modelling Task Force
NSERC-PHAC EID Modelling Consortium (CANMOD, MfPH, OMNI/RÉUNIS)

* The University of Manitoba campuses are located on original lands of Anishinaabeg, Cree, Oji-Cree, Dakota and Dene peoples, and on the homeland of the Métis Nation.

Outline

  • Why stochasticity matters
  • Discrete time Markov chains
  • Continuous time Markov chains

Remarks / Resources

This is a user-oriented course: I barely, if at all, touch on the algorithms; instead, I focus on how to use them

Code is available in this Github repo in the CODE directory

Some of the slides are inspired from slides given to me by Linda Allen (Texas Tech). I recommend books and articles by Linda for more detail

Why stochasticity matters

Running example - SIS model without demography

Constant total population

Basic reproduction number:

In the deterministic world, rules the world

  • If , the disease dies out (disease-free equilibrium)
  • If , it becomes established at an endemic equilibrium

Next slides: K, , (and )

In stochastic world, make that '' rules-ish'' ()

center

When , extinctions happen quite frequently

Types of stochastic systems discussed today

  • Discrete-time Markov chains (DTMC)
  • Continuous-time Markov chains (CTMC)

But there are many others. Of note:

  • Branching processes (BP)
  • Stochastic differential equations (SDE)

Discrete time Markov chains

Discrete-time Markov chains

: probability vector, with describing the probability that at time , the system is in state ,

for al# l , of course

State evolution governed by

where is a stochastic matrix (row sums all equal 1), the transition matrix, with entry , where sequence of random variables describing the state

If constant, homogeneous DTMC

Time often ''recast'' so that

Important remark

The DTMC world lives at the interface between probabilists, who like to think of as a row vector, as a column-stochastic matrix and thus write

and linear algebraists, who prefer column vectors and row-stochastic matrices,

So check the direction to understand whether you are using or

Advantages of DTMC

As a teacher of modelling: base theory of DTMC uses a lot of linear algebra and graph theory; usually really appreciated by students

Regular DTMC (with primitive transition matrices) allow to consider equilibrium distribution of probability

Absorbing DTMC (with reducible transition matrices) allow the consideration of time to absorption, mean first passage time, etc.

See for instance book of Kemeny and Snell

DTMC for example SIS system

Since , consider only the infected. To simulate as DTMC, consider a random walk on ( Gambler's ruin problem)

Denote , and

Transition matrix

# Make the transition matrix
T = mat.or.vec(nr = (Pop+1), nc = (Pop+1))
for (row in 2:Pop) {
  I = row-1
  mv_right = gamma*I*Delta_t # Recoveries
  mv_left = beta*I*(Pop-I)*Delta_t # Infections
  T[row,(row-1)] = mv_right
  T[row,(row+1)] = mv_left
}
# Last row only has move left
T[(Pop+1),Pop] = gamma*(Pop)*Delta_t
# Check that we don't have too large values
if (max(rowSums(T))>1) {
  T = T/max(rowSums(T))
}
diag(T) = 1-rowSums(T)

Simulating a DTMC

library(DTMCPack)
sol = MultDTMC(nchains = 20, tmat = T, io = IC, n = t_f)

gives 20 realisations of a DTMC with transition matrix T, initial conditions IC (a vector of initial probabilities of being in the different states) and final time t_f

See code on Github

Going a bit further

DTMCPack is great for obtaining realisations of a DTMC, but to study it in more detail, markovchain is much more comprehensive

library(markovchain)
mcSIS <- new("markovchain", 
             states = sprintf("I_%d", 0:Pop),
             transitionMatrix = T,
             name = "SIS")

Note that interestingly, markovchain overrides the weird default "* is Hadamard, %*% is usual" R matrix product rule, so mcSIS*mcSIS does , not 

> summary(mcSIS)
SIS  Markov chain that is composed by: 
Closed classes: 
I_0 
Recurrent classes: 
{I_0}
Transient classes: 
{I_1,I_2,I_3,I_4,I_5,I_6,I_7,I_8,I_9,I_10,I_11,I_12,I_13,I_14,I_15,
I_16,I_17,I_18,I_19,I_20,I_21,I_22,I_23,I_24,I_25,I_26,I_27,I_28,
I_29,I_30,I_31,I_32,I_33,I_34,I_35,I_36,I_37,I_38,I_39,I_40,I_41,
I_42,I_43,I_44,I_45,I_46,I_47,I_48,I_49,I_50,I_51,I_52,I_53,I_54,
I_55,I_56,I_57,I_58,I_59,I_60,I_61,I_62,I_63,I_64,I_65,I_66,I_67,
I_68,I_69,I_70,I_71,I_72,I_73,I_74,I_75,I_76,I_77,I_78,I_79,I_80,
I_81,I_82,I_83,I_84,I_85,I_86,I_87,I_88,I_89,I_90,I_91,I_92,I_93,
I_94,I_95,I_96,I_97,I_98,I_99,I_100}
The Markov chain is not irreducible 
The absorbing states are: I_0
Function Role
absorbingStates absorbing states of the transition matrix, if any
steadyStates the vector(s) of steady state(s) in matrix form
meanFirstPassageTime matrix or vector of mean first passage times
meanRecurrenceTime vector of mean number of steps to return to each recurrent state
hittingProbabilities matrix of hitting probabilities for a Markov chain
meanAbsorptionTime expected number of steps for a transient state to be absorbed by any recurrent class
absorptionProbabilities transient states being absorbed by each recurrent state
canonicForm canonic form of transition matrix
period the period of an irreducible DTMC
summary DTMC summary

Continuous time Markov chains

Continuous-time Markov chains

CTMC similar to DTMC except in way they handle time between events (transitions)

DTMC: transitions occur each

CTMC: and transition times follow an exponential distribution parametrised by the state of the system

CTMC are roughly equivalent to ODE

ODE to CTMC : focus on different components

Weight Transition Effect
, new infection of a susceptible
, recovery of an infectious

Will use and omit dynamics

Gillespie's algorithm (SIS model with only I eq.)

set and
while {}

  • Draw from
  • Draw from
  • Find such that
  • switch {}
    • 1: New infection,
    • 2: End of infectious period,

You can also use a package

  • You can implement Gillespie's algorithm yourself, it is a good exercise..
  • But in R there also exists a few packages allowing you to do that easily
  • They have the advantage of implementing tau-leaping (more on this later)

Simulating a CTMC

library(GillespieSSA2)
IC <- c(S = (Pop-I_0), I = I_0)
params <- c(gamma = gamma, beta = beta)
reactions <- list(
  reaction("beta*S*I", c(S=-1,I=+1), "new_infection"),
  reaction("gamma*I", c(S=+1,I=-1), "recovery")
)
set.seed(NULL)
sol <- ssa(
    initial_state = IC,
    reactions = reactions,
    params = params,
    method = ssa_exact(),
    final_time = t_f,
)
plot(sol$time, sol$state[,"I"], type = "l",
     xlab = "Time (days)", ylab = "Number infectious")

Sometimes in a CTMC things go bad

  • Recall that the inter-event time is exponentially distributed
  • Critical step of the Gillespie algorithm:
    • weight of all possible events (propensity)
    • Draw from
  • So the inter-event time if becomes very large for some
  • This can cause the simulation to grind to a halt

Example: a birth and death process

  • Individuals born at per capita rate
  • Individuals die at per capita rate
  • Let's implement this using classic Gillespie

Gillespie's algorithm (birth-death model)

set and
while {}

  • Draw from
  • Draw from
  • Find such that
  • switch {}
    • 1: Birth,
    • 2: Death, 
b = 0.01   # Birth rate
d = 0.01   # Death rate
t_0 = 0    # Initial time
N_0 = 100  # Initial population

# Vectors to store time and state. Initialise with initial condition.
t = t_0
N = N_0

t_f = 1000  # Final time

# We'll track the current time and state (could also just check last entry in t
# and N, but will take more operations)
t_curr = t_0
N_curr = N_0

while (t_curr<=t_f) {
  xi_t = (b+d)*N_curr
  # The exponential number generator does not like a rate of 0 (when the 
  # population crashes), so we check if we need to quit
  if (N_curr == 0) {
    break
  }
  tau_t = rexp(1, rate = xi_t)
  t_curr = t_curr+tau_t
  v = c(b*N_curr, xi_t)/xi_t
  zeta_t = runif(n = 1)
  pos = findInterval(zeta_t, v)+1
  switch(pos,
         { 
           # Birth
           N_curr = N_curr+1
         },
         {
           # Death
           N_curr = N_curr-1
         })
  N = c(N, N_curr)
  t = c(t, t_curr)
}

Drawing at random from an exponential distribution

If you do not have an exponential distribution random number generator.. We want from , i.e., has probability density function

Use cumulative distribution function

which has values in . So draw from and solve for

Last one did not go well

  • Wanted 1000 time units (days?)
  • Interrupted around 500 () because I lost patience
  • (Slide before: the sim stopped because the population went extinct, I did not stop it!)
  • At stop time
    • (and as well, of course!)
    • time was moving slowly
> tail(diff(t))
[1] 1.357242e-04 1.291839e-04 5.926044e-05 7.344020e-05 1.401148e-04 4.423529e-04

Tau-leaping to the rescue!

  • Will not go into details
  • Approximation method (compared to classic Gillespie, which is exact)
  • Roughly: consider "groups" of events instead of individual events
  • Good news: GillespieSSA2 (which we saw earlier) and adaptivetau