Modifying SIR model to include stochasticity

Question

I am trying to establish a method of estimating infectious disease parameters by comparing real epidemic curves with simulations of a stochastic SIR model. To construct the stochastic SIR model, I am using the deSolve package and instead of using fixed parameter values I would like to draw the parameter value used in the equations at each time point from a Poisson distribution centered on the original parameter values.

Using the parameter beta as an example, beta represents the average number of transmission events per capita and is the product of the average number of contacts and the probability that transmission occurs upon contact. Realistically, there is variation in the number of contacts a person will have and since transmission is also a probabilistic event there is variation surrounding this too. So even if the average transmission rate were to be 2.4 (for example), an individual can go on to infect 0, 1, 2 or 3 ... etc. people with varying probabilities.

I have tried to incorporate this into my code below using the rpois function and reassigning the parameters used in the equations to the outputs of the rpois.

I have run my code with the same initial values and parameters multiple times and all the curves are different indicating that SOMETHING "stochastic" is going on, but I am unsure whether the code is sampling using the rpois at each time point or just once at the beginning. I have only started coding very recently so do not have much experience.

I would be grateful if anyone more experienced than myself could verify what my code is ACTUALLY doing and whether it is sampling using rpois at each time point or not. If not I would be grateful for any suggestions for achieving this. Perhaps a loop is needed?

library('deSolve')
library('reshape2')
library('ggplot2')

#MODEL INPUTS
 initial_state_values <- c(S = 10000,
                      I = 1,
                      R = 0)

 #PARAMETERS
  parameters <- c(beta = 2.4,
            gamma = 0.1)


 #POISSON MODELLING OF PARAMETERS
 #BETA
  beta_p <- rpois(1, parameters[1])

 #GAMMA
  infectious_period_p <- rpois(1, 1/(parameters[2]))
 
  gamma_p <- 1/infectious_period_p

 #TIMESTEPS
  times <- seq(from = 0, to = 50,by = 1)

 #SIR MODEL FUNCTION 
  sir_model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
  N <- S + I + R   
  lambda <- beta_p * I/N 
  dS <- -lambda * S
  dI <- lambda*S - gamma_p*I
 dR <- gamma_p*I

return(list(c(dS, dI, dR)))
 })
 }

 output<- as.data.frame(ode(y= initial_state_values,
                       times = times,
                       func = sir_model,
                       parms = parameters))

Do you want to run different simulations with different gamma and beta parameters, or do you want the value of these parameters to change during one single simulation? — Javier, May 31 '21 at 14:25
As @Javier says, this is an important point. In either case: if you use a solver with automatic time step, you can NOT draw the random parameters within the model function. This can work only with Euler method. In all other cases, you should define time variant parameters externally, as a so-called forcing variable. — tpetzoldt, May 31 '21 at 14:41
I would like the parameters to be drawn from a Poisson distribution at each time step, so they can take different values within the same simulation — fil0607, May 31 '21 at 17:41

tpetzoldt · Accepted Answer · 2021-06-24T05:37:05.890

The code given in the question runs the model with constant parameters over time. Here an example with parameters varying over time. However, this setting assumes that for a given time step, the parameters are equal for all indidividuals of the population. If you want to have individual variability, one can either use a matrix formulation for different sub-populations or use an individual model instead.

Model with fluctuating population parameters:

library('deSolve')

initial_state_values <- c(S = 10000,
                          I = 1,
                          R = 0)

parameters <- c(beta = 2.4, gamma = 0.1)

times <- seq(from = 0, to = 50, by = 1) # note time step = 1!

# +1 to add one for time = zero
beta_p <- rpois(max(times) + 1, parameters[1])

infectious_period_p <- rpois(max(times) + 1, 1/(parameters[2]))

gamma_p <- 1/infectious_period_p

sir_model <- function(time, state, parameters) {
  # cat(time, "\n") # show time steps for debugging
  with(as.list(c(state, parameters)), {
    
    # this overwrites the parms passed via parameters
    beta  <- beta_p[floor(time)  + 1]
    gamma <- gamma_p[floor(time) + 1]
    
    N <- S + I + R   
    lambda <- beta * I/N 
    
    dS <- -lambda * S
    dI <- lambda * S - gamma * I
    dR <- gamma * I
    
    list(c(dS, dI, dR))
  })
}

output <- ode(y = initial_state_values,
              times = times,
              func = sir_model,
              parms = parameters)

plot(output)

this is great thank you ! That is exactly what I was looking for... I think my method would be improved by considering variability in parameters at the individual level as you mention. How would I use a matrix formulation to do this? — fil0607, May 31 '21 at 17:44
A SIR model in matrix formulation can be found in [this post](https://stackoverflow.com/questions/65113230/r-desolve-package-ode-function-change-a-matrix-of-parameters-in-sir-model-a). Some addional ideas can be found [here](https://github.com/tpetzoldt/covid/tree/master/models) and also some [links](https://github.com/tpetzoldt/covid). — tpetzoldt, May 31 '21 at 18:05
Yes the matrix formulation will make it more realistic, thank you :) — fil0607, Jun 01 '21 at 11:51
But a matrix formulation can also get somewhat complex. After some thinking, I would personally prefer a stochatic individual-based model instead. Here the question is if your true aim is to simulate a SIR-like model, or if this was only an example. — tpetzoldt, Jun 01 '21 at 13:03

score 2 · Answer 2 · answered May 31 '21 at 18:25

Here another, slightly more generalized version. It is added as a second answer, to keep the original version compact and simple. The new version differs with respect to the following:

generalized, so that it can work with fixed parameters and stochastic forcing
pass parameters as list
run a basic Monte-Carlo simulation

library('deSolve')

sir_model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {

    # this overwrites the parms passed via parameters
    if (time_dependent) {
      beta  <- beta_p[floor(time)  + 1]
      gamma <- gamma_p[floor(time) + 1]
    }
    N <- S + I + R
    lambda <- beta * I/N

    dS <- -lambda * S
    dI <- lambda * S - gamma * I
    dR <- gamma * I

    list(c(dS, dI, dR))
  })
}

initial_state_values <- c(S = 10000, I = 1, R = 0)
times <- seq(from = 0, to = 50, by = 1) # note time step = 1!

## (1) standard simulation with constant parameters
parameters <- c(beta = 2.4, gamma = 0.1)
out0 <- ode(y= initial_state_values,
            times = times,
            func  = sir_model,
            parms = c(parameters, time_dependent = FALSE))
plot(out0)

## (2) single simulation with time varying parameters
beta_p <- rpois(max(times) + 1, parameters[1])
infectious_period_p <- rpois(times + 1, 1/(parameters[2]))
gamma_p <- 1/infectious_period_p

## here we need pass the vectorized parameters globally
## for simplicity, it can also be done as list
out1 <- ode(y = initial_state_values, times = times,
          func = sir_model, parms = c(time_dependent = TRUE))
plot(out0, out1)

## (3) a sample of simulations
monte_carlo <- function(i) {
  #parameters <- c(beta = 2.4, gamma = 0.1)
  beta_p <- rpois(max(times) + 1, parameters[1])
  infectious_period_p <- rpois(max(times) + 1, 1 / (parameters[2]))
  gamma_p <- 1/infectious_period_p
  ode(y = initial_state_values, times = times,
      func = sir_model, parms = list(beta_p  = beta_p,
                                  gamma_p = gamma_p,
                                  time_dependent = TRUE))
}

## run 10 simulations
out_mc <- lapply(1:10, monte_carlo)
plot(out0, out_mc, mfrow=c(1, 3))

Monte Carlo Simulation

Thank you for taking the time to give me such a detailed response, it's extremely useful! — fil0607, Jun 01 '21 at 11:50

Modifying SIR model to include stochasticity

2 Answers2

Linked