Objectives

• Design your favourite model with pen & paper.
• Write down a table with all transitions and jump intensities for this model.
• Implement the model in R using the package adaptivetau
• Simulate and plot stochastic trajectories to investigate on the dynamics of your model.

Design your favourite model

Using only pen & paper, design a compartmental model with boxes and arrows. You can design either:

• A model you would like to work with for a particular disease in your research. In this case, try to choose one not too complicated, say with less than 10 compartments, otherwise it will be tricky to debug.
• The SEITL model, which is one of the 5 models used to analyse the two-wave flu outbreak in Tristan da Cunha. You can use the information provided in this page to design the SEITL model.

Include labels for the rates of flow between compartments and define all the symbols and parameters you use. You should also give a brief description of what the model does as well as some guess values for the parameters and initial conditions.

Write the transition table

Once you have designed your model, write down the associated transition table. The table should contain 3 columns:

• Event: describe the type of transition (e.g. infection)
• Transition: change in the state vector corresponding to the transition (e.g. \((s, i, r)\to(s-1, i+1, r)\))
• Jump intensity: intensity at which the transition occur (e.g. \(\beta si/N\))

If you chose the SEITL model you check your table with our solution.

Implement the model in R

The adaptivetau package can be installed with

install.packages("adaptivetau")

Once installed, it is loaded with

library(adaptivetau)

The adaptivetau package uses a different syntax from the deSolve package. Instead of providing a function to calculate the rates of change at each time point, one specifies a list of transitions and their rates (jump intensities). Examples for how this is done can be found in the adaptivetau vignette.

For the SIR model, we could write

SIR_transitions <- list(
    c(S = -1, I = 1), # infection
    c(I = -1, R = 1) # recovery
    )

SIR_rateFunc <- function(x, parameters, t) {
    
    # define model parameters in term of the natural parameters
    beta <- parameters["R0"]/parameters["D_inf"] 
    nu <- 1/parameters["D_inf"]

    # create temporary variables for states to simplify the writing of the rates below
    S <- x["S"] 
    I <- x["I"] 
    R <- x["R"]
    N <- S + I + R
    
    # return rates
    return(c(
        beta * S * I / N, # infection 
        nu * I # recovery
        )) 
}

You can now re-use and modify the code of the SIR model to implement your model in R. If you chose the SEITL model and you need more help, you can use and complete this partial solution.

Simulate & plot stochastic epidemics

To run the stochastic model, we then use the ssa.adaptivetau function, which takes a vector of initial conditions, the list of transitions and rate function, a named vector of parameters, and the final time (with simulations starting at time 0).

parameters <- c(R0 = 5, D_inf = 1)
run <- ssa.adaptivetau(init.values = c(S = 999, I = 1, R = 0), transitions = SIR_transitions, 
    rateFunc = SIR_rateFunc, params = parameters, tf = 10)
head(run)

##           time   S I R
## [1,] 0.0000000 999 1 0
## [2,] 0.4173075 998 2 0
## [3,] 0.4178569 997 3 0
## [4,] 0.4363569 996 4 0
## [5,] 0.4525435 995 5 0
## [6,] 0.4801278 995 4 1

The simplest way to plot the trajectories is using plot. To plot the output of the stochastic SIR run above, we first convert it to a data frame (ssa.adaptivetau returns a matrix) using data.frame

run_df <- data.frame(run)

We can then plot the number of infected against time using

plot(x = run_df$time, y = run_df$I, type = "line")

Spend some time playing with your stochastic model by changing the parameter values and see how it affects the epidemic. In particular, try a value of \(R_0\) just above 1 to see the early stochastic extinction.

Going further

# get output at days 1, 2, ... , 10
run_I_times <- approx(x = run_df$time, y = run_df$I, xout = 1:10, method = "constant")

print(run_I_times)

## $x
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $y
##  [1]  16 325 363 170  77  32  16   7   3   1