# Discrete compartmental models in a nutshell

## From continuous to discrete time

In its simplest form, a SIR model is typically written in continuous time as:

$\frac{dS}{dt} = - \beta \frac{S_t I_t}{N_t}$

$\frac{dI}{dt} = \beta \frac{S_t I_t}{N_t} - \gamma I_t$

$\frac{dR}{dt} = \gamma I_t$

Where $$\beta$$ is an infection rate and $$\gamma$$ a removal rate, assuming ‘R’ stands for ‘recovered’, which can mean recovery or death.

For discrete time equivalent, we take a small time step $$t$$ (typically a day), and write the changes of individuals in each compartment as:

$S_{t+1} = S_t - \beta \frac{S_t I_t}{N_t}$

$I_{t+1} = I_t + \beta \frac{S_t I_t}{N_t} - \gamma I_t$

$R_{t+1} = R_t + \gamma I_t$

## Stochastic processes

The discrete model above remains deterministic: for given values of the rates $$\beta$$ and $$\gamma$$, dynamics will be fixed. It is fairly straightforward to convert this discrete model into a stochastic one: one merely needs to uses appropriate probability distributions to model the transfer of individuals across compartments. There are at least 3 types of such distributions which will be useful to consider.

### Binomial distribution

This distribution will be used to determine numbers of individuals leaving a given compartment. While we may be tempted to use a Poisson distribution with the rates specified in the equations above, this could lead to over-shooting, i.e. more individuals leaving a compartment than there actually are. To avoid infecting more people than there are susceptibles, we use a binomial distribution, with one draw for each individual in the compartment of interest. The workflow will be:

1. determine a per-capita probability of leaving the compartment, based on the original rates specified in the equations; if the rate at which each individual leaves a compartment is $$\lambda$$, then the corresponding probability of this individual leaving the compartment in one time unit is $$p = 1 - e^{- \lambda}$$.

2. determine the number of individuals leaving the compartment using a Binomial distribution, with one draw per individual and a probability $$p$$

As an example, let us consider transition $$S \rightarrow I$$ in the SIR model. The overall rate at which this change happens is $$\beta \frac{S_t I_t}{N_t}$$. The corresponding per susceptible rate is $$\beta \frac{I_t}{N_t}$$. Therefore, the probability for an individual to move from S to I at time $$t$$ is $$p_{(S \rightarrow I), t} = 1 - e^{- \beta \frac{I_t}{N_t}}$$.

### Poisson distribution

Poisson distributions will be useful when individuals enter a compartment at a given rate, from ‘the outside’. This could be birth or migration (for $$S$$), or introduction of infections from an external reservoir (for $$I$$), etc. The essential distinction with the previous process is that individuals are not leaving an existing compartment.

This case is simple to handle: one just needs to draw new individuals entering the compartment from a Poisson distribution with the rate directly coming from the equations.

For instance, let us now consider a variant of the SIR model where new infectious cases are imported at a constant rate $$\epsilon$$. The only change to the equation is for the infected compartment:

$I_{t+1} = I_t + \beta \frac{S_t I_t}{N_t} + \epsilon - \gamma I_t$

where:

• individuals move from $$S$$ to $$I$$ according to a Binomial distribution $$\mathcal{B}(S_t, 1 - e^{- \beta \frac{I_t}{N_t}})$$

• new infected individuals are imported according to a Poisson distribution $$\mathcal{P}(\epsilon)$$

• individual move from $$I$$ to $$R$$ according to a Binomial distribution $$\mathcal{B}(I_t, 1 - e^{- \gamma})$$

### Multinomial distribution

This distribution will be useful when individuals leaving a compartment are distributed over several compartments. The Multinomial distribution will be used to determine how many individuals end up in each compartment. Let us assume that individuals move from a compartment $$X$$ to compartments $$A$$, $$B$$, and $$C$$, at rates $$\lambda_A$$, $$\lambda_B$$ and $$\lambda_C$$. The workflow to handle these transitions will be:

1. determine the total number of individuals leaving $$X$$; this is done by summing the rates ($$\lambda = \lambda_A + \lambda_B + \lambda_C$$) to compute the per capita probability of leaving $$X$$ $$(p_(X \rightarrow ...) = 1 - e^{- \lambda})$$, and drawing the number of individuals leaving $$X$$ ($$n_{_(X \rightarrow ...)}$$) from a binomial distribution $$n_{(X \rightarrow ...)} \sim B(X, p_(X \rightarrow ...))$$

2. compute relative probabilities of moving to the different compartments (using $$i$$ as a placeholder for $$A$$, $$B$$, $$C$$): $$p_i = \frac{\lambda_i}{\sum_i \lambda_i}$$

3. determine the numbers of individuals moving to $$A$$, $$B$$ and $$C$$ using a Multinomial distribution: $$n_{(X \rightarrow A, B, C)} \sim \mathcal{M}(n_{(X \rightarrow ...)}, p_A, p_B, p_C)$$

# Implementation using odin

## Deterministic SIR model

We start by loading the odin code for a discrete, stochastic SIR model:

As said in the previous vignette, remember this looks and parses like R code, but is not actually R code. Copy-pasting this in a R session will trigger errors.

We then use odin to compile this model:

## <odin_model> object generator
##   Public:
##     initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
##     ir: function ()
##     set_user: function (..., user = list(...), unused_user_action = NULL)
##     initial: function (step)
##     rhs: function (step, y)
##     update: function (step, y)
##     contents: function ()
##     transform_variables: function (y)
##     run: function (step, y = NULL, ..., use_names = TRUE)
##   Private:
##     ptr: NULL
##     use_dde: NULL
##     odin: NULL
##     variable_order: NULL
##     output_order: NULL
##     n_out: NULL
##     ynames: NULL
##     interpolate_t: NULL
##     cfuns: list
##     dll: discrete.deterministic.sir7184bd4d
##     user: I_ini S_ini beta gamma
##     registration: function ()
##     set_initial: function (step, y, use_dde)
##   Parent env: <environment: namespace:discrete.deterministic.sir7184bd4d>
##   Locked objects: TRUE
##   Locked class: FALSE
##   Portable: TRUE

Note: this is the slow part (generation and then compilation of C code)! Which means for computer-intensive work, the number of times this is done should be minimized.

The returned object sir_generatoris an R6 generator that can be used to create an instance of the model: generate an instance of the model:

## <odin_model>
##   Public:
##     contents: function ()
##     initial: function (step)
##     initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
##     ir: function ()
##     rhs: function (step, y)
##     run: function (step, y = NULL, ..., use_names = TRUE)
##     set_user: function (..., user = list(...), unused_user_action = NULL)
##     transform_variables: function (y)
##     update: function (step, y)
##   Private:
##     cfuns: list
##     dll: discrete.deterministic.sir7184bd4d
##     interpolate_t: NULL
##     n_out: 0
##     odin: environment
##     output_order: NULL
##     ptr: externalptr
##     registration: function ()
##     set_initial: function (step, y, use_dde)
##     use_dde: FALSE
##     user: I_ini S_ini beta gamma
##     variable_order: list
##     ynames: step S I R

x is an ode_model object which can be used to generate dynamics of a discrete-time, deterministic SIR model. This is achieved using the function x\$run(), providing time steps as single argument, e.g.:

##       step         S        I         R
##  [1,]    0 1000.0000 1.000000 0.0000000
##  [2,]    1  999.8002 1.099800 0.1000000
##  [3,]    2  999.5805 1.209517 0.2099800
##  [4,]    3  999.3389 1.330125 0.3309317
##  [5,]    4  999.0734 1.462696 0.4639442
##  [6,]    5  998.7814 1.608403 0.6102138
##  [7,]    6  998.4604 1.768530 0.7710541
##  [8,]    7  998.1076 1.944486 0.9479071
##  [9,]    8  997.7198 2.137811 1.1423557
## [10,]    9  997.2937 2.350191 1.3561368
## [11,]   10  996.8254 2.583469 1.5911558

## Stochastic SIR model

The stochastic equivalent of the previous model can be formulated in odin as follows:

We can use the same workflow as before to run the model, using 10 initial infected individuals (I_ini = 10):

## <odin_model> object generator
##   Public:
##     initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
##     ir: function ()
##     set_user: function (..., user = list(...), unused_user_action = NULL)
##     initial: function (step)
##     rhs: function (step, y)
##     update: function (step, y)
##     contents: function ()
##     transform_variables: function (y)
##     run: function (step, y = NULL, ..., use_names = TRUE)
##   Private:
##     ptr: NULL
##     use_dde: NULL
##     odin: NULL
##     variable_order: NULL
##     output_order: NULL
##     n_out: NULL
##     ynames: NULL
##     interpolate_t: NULL
##     cfuns: list
##     dll: discrete.stochastic.sir8ccbfa6f
##     user: I_ini S_ini beta gamma
##     registration: function ()
##     set_initial: function (step, y, use_dde)
##   Parent env: <environment: namespace:discrete.stochastic.sir8ccbfa6f>
##   Locked objects: TRUE
##   Locked class: FALSE
##   Portable: TRUE

This gives us a single stochastic realisation of the model, which is of limited interest. As an alternative, we can generate a large number of replicates using arrays for each compartment:

## <odin_model> object generator
##   Public:
##     initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
##     ir: function ()
##     set_user: function (..., user = list(...), unused_user_action = NULL)
##     initial: function (step)
##     rhs: function (step, y)
##     update: function (step, y)
##     contents: function ()
##     transform_variables: function (y)
##     run: function (step, y = NULL, ..., use_names = TRUE)
##   Private:
##     ptr: NULL
##     use_dde: NULL
##     odin: NULL
##     variable_order: NULL
##     output_order: NULL
##     n_out: NULL
##     ynames: NULL
##     interpolate_t: NULL
##     cfuns: list
##     dll: discrete.stochastic.sir.arrays4ec523f1
##     user: I_ini S_ini beta gamma nsim
##     registration: function ()
##     set_initial: function (step, y, use_dde)
##   Parent env: <environment: namespace:discrete.stochastic.sir.arrays4ec523f1>
##   Locked objects: TRUE
##   Locked class: FALSE
##   Portable: TRUE

## A stochastic SEIRDS model

This model is a more complex version of the previous one, which we will use to illustrate the use of all distributions mentioned in the first part: Binomial, Poisson and Multinomial.

The model is contains the following compartments:

• $$S$$: susceptibles
• $$E$$: exposed, i.e. infected but not yet contagious
• $$I_R$$: infectious who will survive
• $$I_D$$: infectious who will die
• $$R$$: recovered
• $$D$$: dead

There are no birth of natural death processes in this model. Parameters are:

• $$\beta$$: rate of infection
• $$\delta$$: rate at which symptoms appear (i.e inverse of mean incubation period)
• $$\gamma_R$$: recovery rate
• $$\gamma_D$$: death rate
• $$\mu$$: case fatality ratio (proportion of cases who die)
• $$\epsilon$$: import rate of infected individuals (applies to $$E$$ and $$I$$)
• $$\omega$$: rate waning immunity

The model will be written as:

$S_{t+1} = S_t - \beta \frac{S_t (I_{R,t} + I_{D,t})}{N_t} + \omega R_t$

$E_{t+1} = E_t + \beta \frac{S_t (I_{R,t} + I_{D,t})}{N_t} - \delta E_t + \epsilon$

$I_{R,t+1} = I_{R,t} + \delta (1 - \mu) E_t - \gamma_R I_{R,t} + \epsilon$

$I_{D,t+1} = I_{D,t} + \delta \mu E_t - \gamma_D I_{D,t} + \epsilon$

$R_{t+1} = R_t + \gamma_R I_{R,t} - \omega R_t$

$D_{t+1} = D_t + \gamma_D I_{D,t}$

The formulation of the model in odin is:

## <odin_model> object generator
##   Public:
##     initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
##     ir: function ()
##     set_user: function (..., user = list(...), unused_user_action = NULL)
##     initial: function (step)
##     rhs: function (step, y)
##     update: function (step, y)
##     contents: function ()
##     transform_variables: function (y)
##     run: function (step, y = NULL, ..., use_names = TRUE)
##   Private:
##     ptr: NULL
##     use_dde: NULL
##     odin: NULL
##     variable_order: NULL
##     output_order: NULL
##     n_out: NULL
##     ynames: NULL
##     interpolate_t: NULL
##     cfuns: list
##     dll: discrete.stochastic.seirdsf8855d06
##     user: E_ini S_ini beta delta epsilon gamma_D gamma_R mu omega
##     registration: function ()
##     set_initial: function (step, y, use_dde)
##   Parent env: <environment: namespace:discrete.stochastic.seirdsf8855d06>
##   Locked objects: TRUE
##   Locked class: FALSE
##   Portable: TRUE

Several runs can be obtained without rewriting the model, for instance, to get 100 replicates:

## [1] 366 600
##    S.1 E.1 Id.1 Ir.1 R.1 D.1  S.2 E.2 Id.2 Ir.2
## 1 1000   1    0    0   0   0 1000   1    0    0
## 2 1000   1    0    0   0   0 1000   1    0    0
## 3 1000   0    0    1   0   0 1000   2    0    0
## 4 1000   0    0    1   0   0 1000   1    1    0
## 5 1000   0    0    1   0   0 1000   1    1    0
## 6 1000   0    0    1   0   0 1000   0    2    0

It is then possible to explore the behaviour of the model using a simple function:

This is a sanity check with a null infection rate and no imported case:

Another easy case: no importation, no waning immunity:

A more nuanced case: persistence of the disease with limited import, waning immunity, low severity, larger population: