|
|
|
# Simulating dynamics under a SIR model
|
|
|
|
|
|
|
|
The SIR epidemiological model describes a population divided into three compartments: the S compartment containing individuals susceptible to a disease, the I compartment containing individuals infected with the disease, and the R compartment containing individuals who have recovered from the disease. Individuals in these compartments are subject to two distinct processes. In this model, there are two events/reactions :
|
|
|
|
1. the infection at a transmission rate $\beta$ that requires contact between a susceptible individual and an infected one and that generates a new infected individual,
|
|
|
|
2. the recovery at a rate $\gamma$ where an infected individual is moved to the R compartment.
|
|
|
|
|
|
|
|
This SIR model can be described by a system of reactions such as
|
|
|
|
|
|
|
|
```math
|
|
|
|
S+I \xrightarrow[]{\beta S I} I+I \\
|
|
|
|
I \xrightarrow[]{\gamma I} R
|
|
|
|
```
|
|
|
|
|
|
|
|
or by a system of differential equations such as
|
|
|
|
|
|
|
|
```math
|
|
|
|
\begin{aligned}
|
|
|
|
\frac{dS}{dt} &= -\beta S I \\
|
|
|
|
\frac{dI}{dt} &= \beta S I - \gamma I
|
|
|
|
\end{aligned}
|
|
|
|
```
|
|
|
|
|
|
|
|
## Defining the SIR model in TiPS
|
|
|
|
|
|
|
|
In R, with `TiPS`, this model can be written as
|
|
|
|
|
|
|
|
``` r
|
|
|
|
reactions <- c("S [beta*S*I] -> I",
|
|
|
|
"I [gamma*I] -> R")
|
|
|
|
```
|
|
|
|
|
|
|
|
Let's now build the simulator:
|
|
|
|
|
|
|
|
We then build the simulator that will allow us to run multiple trajectories:
|
|
|
|
|
|
|
|
``` r
|
|
|
|
sir_simu <- build_simulator(reactions)
|
|
|
|
```
|
|
|
|
|
|
|
|
This typically takes 10-15' as it involves compilation.
|
|
|
|
|
|
|
|
## Defining simulations parameters
|
|
|
|
-------------------------------
|
|
|
|
|
|
|
|
To run numerical simulations, we define the initial values of the state variables,
|
|
|
|
|
|
|
|
``` r
|
|
|
|
initialStates <- c(I = 1, S = 9999, R = 0)
|
|
|
|
```
|
|
|
|
|
|
|
|
the time range of the simulations,
|
|
|
|
|
|
|
|
``` r
|
|
|
|
time <- c(0, 20)
|
|
|
|
```
|
|
|
|
|
|
|
|
and the parameters values
|
|
|
|
|
|
|
|
``` r
|
|
|
|
theta <- list(gamma = 1, beta = 2e-4)
|
|
|
|
```
|
|
|
|
|
|
|
|
For the *τ*-leap and mixed algorithms, a time step is also required:
|
|
|
|
|
|
|
|
``` r
|
|
|
|
dT <- 0.001
|
|
|
|
```
|
|
|
|
|
|
|
|
## Running simulations
|
|
|
|
|
|
|
|
In some simulations, the total rate may be zero before the upper time limit is reached, because of stochasticity or parameter values. In this case, the simulation is considered to have failed and is halted.
|
|
|
|
|
|
|
|
To bypass these failures, we can define the following wrapper:
|
|
|
|
|
|
|
|
``` r
|
|
|
|
safe_run <- function(f, ...) {
|
|
|
|
out <- list()
|
|
|
|
while(! length(out)) {out <- f(...)}
|
|
|
|
out
|
|
|
|
}
|
|
|
|
```
|
|
|
|
|
|
|
|
A safe version of our simulator `sir_simu()` is then:
|
|
|
|
|
|
|
|
``` r
|
|
|
|
safe_sir_simu <- function(...) safe_run(sir_simu, ...)
|
|
|
|
```
|
|
|
|
|
|
|
|
### Direct method
|
|
|
|
|
|
|
|
A trajectory using Gillespie's direct method is obtained by
|
|
|
|
|
|
|
|
``` r
|
|
|
|
traj_dm <- safe_sir_simu(
|
|
|
|
paramValues = theta,
|
|
|
|
initialStates = initialStates,
|
|
|
|
times = time,
|
|
|
|
method = "exact")
|
|
|
|
```
|
|
|
|
|
|
|
|
The output consists of a named list containing the reactions of the model (with `$reactions`), the parameter values (with `$values`), the time range (with `$times`), the algorithm used to simulate (with `$algo`), the time-step in case the algorithm is *τ*-leap or the mixed algorithm (with `$dT`) and finally the simulated trajectory (with `$traj`) :
|
|
|
|
|
|
|
|
``` r
|
|
|
|
names(traj_dm)
|
|
|
|
#> [1] "reactions" "values" "times" "method" "tau" "seed"
|
|
|
|
#> [7] "traj"
|
|
|
|
```
|
|
|
|
|
|
|
|
The simulated trajectory is also a named list, where each simulated reaction event is recorded `$Reaction`, along with the time at which it occured `$Time`, the number of times it occured `$Nrep` (if *τ*-leap or mixed algorithm chosen), and the size of each compartment througt time, here `$I` `$R` `$S`.
|
|
|
|
|
|
|
|
``` r
|
|
|
|
head(traj_dm$traj)
|
|
|
|
#> Time Reaction Nrep S I R
|
|
|
|
#> 1 0.00000000 init 1 9999 1 0
|
|
|
|
#> 2 0.09966524 S [beta*S*I] -> I 1 9998 2 0
|
|
|
|
#> 3 0.13240687 I [gamma*I] -> R 1 9998 1 1
|
|
|
|
#> 4 0.63636802 S [beta*S*I] -> I 1 9997 2 1
|
|
|
|
#> 5 0.76161520 S [beta*S*I] -> I 1 9996 3 1
|
|
|
|
#> 6 0.76990828 S [beta*S*I] -> I 1 9995 4 1
|
|
|
|
```
|
|
|
|
|
|
|
|
The trajectory can readily be plotted using the `plot()` function:
|
|
|
|
|
|
|
|
``` r
|
|
|
|
plot(traj_dm)
|
|
|
|
```
|
|
|
|
|
|
|
|
<img src="vignettes/TiPS_files/figure-markdown_github/unnamed-chunk-13-1.png" width="576" style="display: block; margin: auto;" />
|
|
|
|
|