library(EpidSim)
library(ape)
library(magrittr)

Simulating trajectories of an SIR model

Building the simulator

In epidsim, the reactions of the epidemiological model are described by the following formalism:

reactions <- list("S+I+=[beta*S*I]I+I", "I=[gamma*I]R")

One can alternatively adopt the formalism of the deSolve package and then translate to the formalism of the epidsim package thanks to the read.reactions() function:

reactions <- read.reactions(c("dI =  beta * S * I - gamma * I",
                              "dS = -beta * S * I",
                              "dR =  gamma * I"))

which gives:

reactions
#> [[1]]
#> [1] "I=[gamma*I]R"
#> 
#> [[2]]
#> [1] "S+I+=[beta*S*I]I+I"

Let’s now define the names of the parameters:

paramNames <- list("beta", "gamma")

and of the compartments:

indivNames <- list("S", "*I", "R")

Note the * in front of the deme compartment. Next step is to build the simulator from reactions, paramNames and indivNames:

sir.simu <- build.simulator(reactions, paramNames, indivNames)

This typically takes 10-15’ as it involves compilation.

Defining simulations parameters

Let’s now define the initial values of the state variables:

initialStates <- c(I = 1, S = 9999, R = 0)

the time range of the simulations:

time <- c(0, 20)

and the parameters values:

theta <- list(c(gamma = 1, beta = 2e-4))

as well as the time step required by the \(\tau\)-leap and mixed algorithms:

dT <- .08

Running simulations

Sometimes the simulations fail. In order to bypass these failures, let’s define the following wrapper:

safe_run <- function(f, ...) {
  out <- list()
  while(! length(out)) {out <- f(...)}
  out
}

And a safe version of sir.simu() is:

safe_sir.simu <- function(...) safe_run(sir.simu, ...)

Direct method

traj_dm <- safe_sir.simu(
  paramValues = theta,
  initialStates = initialStates,
  times = time)

which gives

head(traj_dm)
#>        Time           Reaction Nrep    S I R
#> 1 0.0000000               init    1 9999 1 0
#> 2 0.6293988 S+I+=[beta*S*I]I+I    1 9998 2 0
#> 3 0.7088931       I=[gamma*I]R    1 9998 1 1
#> 4 0.9973502 S+I+=[beta*S*I]I+I    1 9997 2 1
#> 5 1.0454721       I=[gamma*I]R    1 9997 1 2
#> 6 1.0603469 S+I+=[beta*S*I]I+I    1 9996 2 2

and

plot(I ~ Time, traj_dm, type = "l", col = "red")

\(\tau\)-leap method

traj_tl <- safe_sir.simu(
  paramValues = theta,
    initialStates = initialStates,
    times = time,
    dT = dT)

which gives

head(traj_tl)
#>   Time           Reaction Nrep    S I R
#> 1 0.00               init    1 9999 1 0
#> 2 0.40 S+I+=[beta*S*I]I+I    1 9998 2 0
#> 3 0.80 S+I+=[beta*S*I]I+I    1 9997 3 0
#> 4 0.88 S+I+=[beta*S*I]I+I    1 9996 4 0
#> 5 0.96 S+I+=[beta*S*I]I+I    2 9994 6 0
#> 6 1.04 S+I+=[beta*S*I]I+I    1 9993 7 0

and

plot(I ~ Time, traj_tl, type = "l", col = "blue")

Mixed method

traj_mm <- safe_sir.simu(
    paramValues = theta,
    initialStates = initialStates,
    times = time,
    dT = dT,
    switchingMode = TRUE,
    nTrials = 2)

which gives:

head(traj_mm)
#>         Time           Reaction Nrep    S I R
#> 1 0.00000000               init    1 9999 1 0
#> 2 0.09781542 S+I+=[beta*S*I]I+I    1 9998 2 0
#> 3 0.23492386 S+I+=[beta*S*I]I+I    1 9997 3 0
#> 4 0.24771052 S+I+=[beta*S*I]I+I    1 9996 4 0
#> 5 0.47468083       I=[gamma*I]R    1 9996 3 1
#> 6 0.63350217 S+I+=[beta*S*I]I+I    1 9995 4 1

and

plot(I ~ Time, traj_mm, type = "l", col = "green")

Comparing the three methods

nb <- 100
traj_dm100 <- purrr::rerun(nb,
                           safe_sir.simu(paramValues = theta,
                                         initialStates = initialStates,
                                         times = time))
traj_tl100 <- purrr::rerun(nb,
                           safe_sir.simu(paramValues = theta,
                                         initialStates = initialStates,
                                         times = time,
                                         dT = dT))
traj_mm100 <- purrr::rerun(nb,
                           safe_sir.simu(paramValues = theta,
                                         initialStates = initialStates,
                                         times = time,
                                         dT = dT,
                                         switchingMode = TRUE,
                                         nTrials = 2))
frame <- function(...) {
  plot(NA, xlab = "Time", ylab = "I", xlim = c(0, 20), ...)
}
ylim <- c(0, max(unlist(purrr::map(c(traj_dm100, traj_tl100, traj_mm100), ~ .[["I"]]))))
alpha <- .5
frame(ylim = ylim)
purrr::walk(traj_dm100[1:10], with, lines(Time, I, col = adjustcolor("red", alpha)))
purrr::walk(traj_tl100[1:10], with, lines(Time, I, col = adjustcolor("blue", alpha)))
purrr::walk(traj_mm100[1:10], with, lines(Time, I, col = adjustcolor("green", alpha)))
legend("topright", c("direct method", "tau-leap", "mixed method"),
       col = c("red", "blue", "green"), bty = "n", lty = 1, lwd = 2)

frame(ylim = ylim)
purrr::walk(traj_dm100, with, lines(Time, I, col = adjustcolor("red", alpha)))

frame(ylim = ylim)
purrr::walk(traj_tl100, with, lines(Time, I, col = adjustcolor("blue", alpha)))

frame(ylim = ylim)
purrr::walk(traj_mm100, with, lines(Time, I, col = adjustcolor("green", alpha)))

Comparing with the adaptivetau package

sir <- function(f = adaptivetau::ssa.adaptivetau) {
  require(magrittr)
  
  transitions <- list(c(S = -1, I = +1),
                      c(I = -1, R = +1))
  
  lvrates <- function(x, params, t) {
    with(c(x, params),
         c(beta * S * I,
           gamma * I)
        )
  }
  
  f(initialStates, transitions, lvrates, as.list(unlist(theta)), tf = max(time)) %>% 
    data.frame() %>% 
    dplyr::rename(Time = time)
}
sir_dm <- purrr::rerun(nb, sir(adaptivetau::ssa.exact))
sir_tl <- purrr::rerun(nb, sir(adaptivetau::ssa.adaptivetau))
frame(ylim = ylim)
purrr::walk(sir_dm, with, lines(Time, I, col = adjustcolor("red", alpha)))

frame(ylim = ylim)
purrr::walk(sir_tl, with, lines(Time, I, col = adjustcolor("blue", alpha)))

The following function calculates quantiles from a list of simulation:

quantiles <- function(x) {
  require(magrittr)
  x %>% 
    dplyr::bind_rows() %>% 
    dplyr::select(Time, I) %>%
    dplyr::mutate(Time = cut(Time, 0:ceiling(max(Time)), include.lowest = TRUE)) %>% 
    dplyr::group_by(Time) %>% 
    dplyr::summarise(list(tibble::enframe(quantile(I, c(.25, .5, .75)), "probs", "quantile"))) %>% 
    tidyr::unnest() %>% 
    tidyr::spread(probs, quantile) %>% 
    setNames(c("time", "lwr", "median", "upr"))
}

Let’s try it:

(traj_dm100q <- quantiles(traj_dm100))
#> # A tibble: 22 x 4
#>    time     lwr median   upr
#>    <fct>  <dbl>  <dbl> <dbl>
#>  1 [0,1]      2      4     6
#>  2 (1,2]      7     13    21
#>  3 (2,3]     20     36    62
#>  4 (3,4]     53    100   164
#>  5 (4,5]    136    252   402
#>  6 (5,6]    301    525   791
#>  7 (6,7]    595    905  1227
#>  8 (7,8]    954   1295  1483
#>  9 (8,9]   1246   1425  1509
#> 10 (9,10]  1181   1372  1483
#> # … with 12 more rows

And for the other ones:

traj_tl100q <- quantiles(traj_tl100)
traj_mm100q <- quantiles(traj_mm100)
sir_dm_q <- quantiles(sir_dm)
sir_tl_q <- quantiles(sir_tl)

The following function the quantiles to a plot:

add_distribution <- function(x, col = "red", alpha = .5) {
  with(x, {
    polygon(c(time, rev(time)), c(lwr, rev(upr)),
            col = adjustcolor(col, alpha), border = NA)
    lines(time, median, col = col)
  })
}

Let’s try it: comparing the three methods of the EpidSim package:

alpha2 <- .3
frame(ylim = ylim)
add_distribution(traj_dm100q, "red", alpha2)
add_distribution(traj_tl100q, "blue", alpha2)
add_distribution(traj_mm100q, "green", alpha2)
legend("topright", c("direct method", "tau-leap", "mixed method"),
       fill = adjustcolor(c("red", "blue", "green"), alpha2), bty = "n")

Comparing EpidSim and adaptivetau packages, for the direct method:

frame(ylim = ylim)
add_distribution(traj_dm100q, "red", alpha2)
add_distribution(sir_dm_q, "blue", alpha2)
legend("topright", c("EpidSim", "adaptivetau"),
       fill = adjustcolor(c("red", "blue"), alpha2), bty = "n")

and for the tau-leap method:

frame(ylim = ylim)
add_distribution(traj_tl100q, "red", alpha2)
add_distribution(sir_tl_q, "blue", alpha2)
legend("topright", c("EpidSim", "adaptivetau"),
       fill = adjustcolor(c("red", "blue"), alpha2), bty = "n")

TO DO: comparing execution times. But hard to do when EpidSim sometimes fails.

Simulating phylogenies with the coalescent

Here we consider a simple BD process with 10 time intervals.

Simulating a trajectory

Specifying the reactions:

reactions <- list("I+=[beta*I]I+I", "I-=[gamma*I]")

Building the simulator:

gil.simu <- build.simulator(
  Reactions = reactions,
  ParameterNames = list("R0", "gamma", "beta"),
  CompartmentNames = list("*I"))

A safe version of it:

safe_gil.simu <- function(...) safe_run(gil.simu, ...)

Initial states:

initialStates <- c(I = 1)

A list of parameters values:

theta <- tibble::tribble(
   ~R0,    ~beta,
  2.15, .8827628,
  1.74, .714422,
  1.80, .7390572,
  1.43, .5871399,
  1.35, .5542929,
  1.54, .6323045,
  2.03, .8334923,
  1.39, .5707164,
  1.26, .51734,
  1.78, .7308455) %$%
  purrr::map2(R0, beta, ~ c(R0 = .x, gamma = .410587333412841, beta = .y))

Time values:

times <- seq(1926.17, by = 8.683, le = 11)

dT:

dT <- .08

\(\tau\)-leap method:

trjctr_tl <- safe_gil.simu(
  paramValues = theta,
  initialStates = initialStates,
  times = times,
  dT = dT,
  info = TRUE,
  seed = 119700)

Mixed method:

trjctr_mm <- safe_gil.simu(
  paramValues = theta,
  initialStates = initialStates,
  times = times,
  dT = dT,
  info = FALSE,
  switchingMode = TRUE,
  nTrials = 10)

Let’s compare the simulations:

plot(I ~ Time, trjctr_tl, type = "l", col = "red")
with(trjctr_mm, lines(Time, I, col = "blue"))
legend("center", c("tau-leap", "mixed method"),
       col = c("red", "blue"), bty = "n", lty = 1, lwd = 2)

Simulating a phylogeny

The sampling dates are in the following text file:

dates <- system.file("extdata", "BD-10_dates.txt", package = "EpidSim")

Here we used the simulated trajectory to simulate a phylogeny using the coalescent:

simulate.tree(
  reactions = reactions,
  trajectory = trjctr_tl,
  dates = dates,
  sampled = list("I" = 1), # the type of individuals that are sampled and their proportion (??)
  root = "I", # type of individual at the root of the tree
  outputTreeFile = "BD_Backward_Tree.newick",
  isFullTrajectory = FALSE, # deads do not generate leaves
  seed = 0,
  nTrials = 1,
  addInfos = FALSE, # additional info for each node
  nexus = FALSE) # output format is newick

Reading the simulated tree and visualizing it:

treeFig <- read.tree("BD_Backward_Tree.newick")
ape::plot.phylo(treeFig, cex = .5)

Forward simulation of phylogenies

Here we simulate phylogenies at the same time as the trajectories.

Reactions du modèle SIS

reactions <- list("S+I+=[beta*S*I]I+I", "I=[gamma*I]S")

Building the simulator.

This produces an ERROR:

forward.sis <- build.forward.tree.simulation(reactions,
                                             list("gamma", "beta"),
                                             list("*I", "S"))

Simulating the trajectory:

traj <- forward.sis(
  paramValues = list(c(beta = .2, gamma = .1)),
  initialStates = c(I = 1, S = 99999),
  times = c(1980, 2020),
  treeFile = "BD_Backward_Tree.newick",
  samplingDates = system.file("extdata", "SIS-6.dates", package = "EpidSim"),
  sampled = list("I" = 1))