From continuous to discrete time

In its simplest form, a SIR model is typically written in continuoustime 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 rateand\(\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 thechanges 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 ofthe rates\(\beta\) and\(\gamma\), dynamics will be fixed. It isfairly straightforward to convert this discrete model into a stochasticone: one merely needs to uses appropriate probability distributions tomodel the transfer of individuals across compartments. There are atleast 3 types of such distributions which will be useful toconsider.

Deterministic SIR model

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

path_sir_model<-system.file("examples/discrete_deterministic_sir.R",package ="odin")

## Core equations for transitions between compartments:update(S)<- S- beta* S* I/ Nupdate(I)<- I+ beta* S* I/ N- gamma* Iupdate(R)<- R+ gamma* I## Total population size (odin will recompute this at each timestep:## automatically)N<- S+ I+ R## Initial states:initial(S)<- S_ini# will be user-definedinitial(I)<- I_ini# will be user-definedinitial(R)<-0## User defined parameters - default in parentheses:S_ini<-user(1000)I_ini<-user(1)beta<-user(0.2)gamma<-user(0.1)

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

We then useodin to compile this model:

sir_generator<- odin::odin(path_sir_model)

## Loading required namespace: pkgbuild

sir_generator

## <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.sir4de329d6##     user: I_ini S_ini beta gamma##     registration: function () ##     set_initial: function (step, y, use_dde) ##     update_metadata: function () ##   Parent env: <environment: namespace:discrete.deterministic.sir4de329d6>##   Locked objects: TRUE##   Locked class: FALSE##   Portable: TRUE

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

The returned objectsir_generatoris an R6 generator thatcan be used to create an instance of the model: generate an instance ofthe model:

x<- sir_generator$new()x

## <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.sir4de329d6##     interpolate_t: NULL##     n_out: 0##     odin: environment##     output_order: NULL##     ptr: externalptr##     registration: function () ##     set_initial: function (step, y, use_dde) ##     update_metadata: function () ##     use_dde: FALSE##     user: I_ini S_ini beta gamma##     variable_order: list##     ynames: step S I R

x is anode_model object which can be usedto generate dynamics of a discrete-time, deterministic SIR model. Thisis achieved using the functionx$run(), providing timesteps as single argument, e.g.:

sir_col<-c("#8c8cd9","#cc0044","#999966")x$run(0:10)

##       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

x_res<- x$run(0:200)par(mar =c(4.1,5.1,0.5,0.5),las =1)matplot(x_res[,1], x_res[,-1],xlab ="Time",ylab ="Number of individuals",type ="l",col = sir_col,lty =1)legend("topright",lwd =1,col = sir_col,legend =c("S","I","R"),bty ="n")

An example of deterministic, discrete-timeSIR model

Stochastic SIR model

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

path_sir_model_s<-system.file("examples/discrete_stochastic_sir.R",package ="odin")

## Core equations for transitions between compartments:update(S)<- S- n_SIupdate(I)<- I+ n_SI- n_IRupdate(R)<- R+ n_IR## Individual probabilities of transition:p_SI<-1-exp(-beta* I/ N)# S to Ip_IR<-1-exp(-gamma)# I to R## Draws from binomial distributions for numbers changing between## compartments:n_SI<-rbinom(S, p_SI)n_IR<-rbinom(I, p_IR)## Total population sizeN<- S+ I+ R## Initial states:initial(S)<- S_iniinitial(I)<- I_iniinitial(R)<-0## User defined parameters - default in parentheses:S_ini<-user(1000)I_ini<-user(1)beta<-user(0.2)gamma<-user(0.1)

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

sir_s_generator<- odin::odin(path_sir_model_s)sir_s_generator

## <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.sir09143884##     user: I_ini S_ini beta gamma##     registration: function () ##     set_initial: function (step, y, use_dde) ##     update_metadata: function () ##   Parent env: <environment: namespace:discrete.stochastic.sir09143884>##   Locked objects: TRUE##   Locked class: FALSE##   Portable: TRUE

x<- sir_s_generator$new(I_ini =10)

set.seed(1)x_res<- x$run(0:100)par(mar =c(4.1,5.1,0.5,0.5),las =1)matplot(x_res[,1], x_res[,-1],xlab ="Time",ylab ="Number of individuals",type ="l",col = sir_col,lty =1)legend("topright",lwd =1,col = sir_col,legend =c("S","I","R"),bty ="n")

An example of stochastic, discrete-time SIRmodel

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

path_sir_model_s_a<-system.file("examples/discrete_stochastic_sir_arrays.R",package ="odin")

## Core equations for transitions between compartments:update(S[])<- S[i]- n_SI[i]update(I[])<- I[i]+ n_SI[i]- n_IR[i]update(R[])<- R[i]+ n_IR[i]## Individual probabilities of transition:p_SI[]<-1-exp(-beta* I[i]/ N[i])p_IR<-1-exp(-gamma)## Draws from binomial distributions for numbers changing between## compartments:n_SI[]<-rbinom(S[i], p_SI[i])n_IR[]<-rbinom(I[i], p_IR)## Total population sizeN[]<- S[i]+ I[i]+ R[i]## Initial states:initial(S[])<- S_iniinitial(I[])<- I_iniinitial(R[])<-0## User defined parameters - default in parentheses:S_ini<-user(1000)I_ini<-user(1)beta<-user(0.2)gamma<-user(0.1)## Number of replicatesnsim<-user(100)dim(N)<- nsimdim(S)<- nsimdim(I)<- nsimdim(R)<- nsimdim(p_SI)<- nsimdim(n_SI)<- nsimdim(n_IR)<- nsim

sir_s_a_generator<- odin::odin(path_sir_model_s_a)sir_s_a_generator

## <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.arraysb818bd3a##     user: I_ini S_ini beta gamma nsim##     registration: function () ##     set_initial: function (step, y, use_dde) ##     update_metadata: function () ##   Parent env: <environment: namespace:discrete.stochastic.sir.arraysb818bd3a>##   Locked objects: TRUE##   Locked class: FALSE##   Portable: TRUE

x<- sir_s_a_generator$new()

set.seed(1)sir_col_transp<-paste0(sir_col,"66")x_res<- x$run(0:100)par(mar =c(4.1,5.1,0.5,0.5),las =1)matplot(x_res[,1], x_res[,-1],xlab ="Time",ylab ="Number of individuals",type ="l",col =rep(sir_col_transp,each =100),lty =1)legend("left",lwd =1,col = sir_col,legend =c("S","I","R"),bty ="n")

100 replicates of a stochastic, discrete-timeSIR model

A stochastic SEIRDS model

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

The model is contains the following compartments:

• \(S\): susceptibles
• \(E\): exposed, i.e. infected butnot yet contagious
• \(I_R\): infectious who willsurvive
• \(I_D\): infectious who willdie
• \(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 symptomsappear (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 ofinfected individuals (applies to\(E\)and\(I\))
• \(\omega\): rate waningimmunity

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 inodin is:

path_seirds_model<-system.file("examples/discrete_stochastic_seirds.R",package ="odin")

## Core equations for transitions between compartments:update(S)<- S- n_SE+ n_RSupdate(E)<- E+ n_SE- n_EI+ n_import_Eupdate(Ir)<- Ir+ n_EIr- n_IrRupdate(Id)<- Id+ n_EId- n_IdDupdate(R)<- R+ n_IrR- n_RSupdate(D)<- D+ n_IdD## Individual probabilities of transition:p_SE<-1-exp(-beta* I/ N)p_EI<-1-exp(-delta)p_IrR<-1-exp(-gamma_R)# Ir to Rp_IdD<-1-exp(-gamma_D)# Id to dp_RS<-1-exp(-omega)# R to S## Draws from binomial distributions for numbers changing between## compartments:n_SE<-rbinom(S, p_SE)n_EI<-rbinom(E, p_EI)n_EIrId[]<-rmultinom(n_EI, p)p[1]<-1- mup[2]<- mudim(p)<-2dim(n_EIrId)<-2n_EIr<- n_EIrId[1]n_EId<- n_EIrId[2]n_IrR<-rbinom(Ir, p_IrR)n_IdD<-rbinom(Id, p_IdD)n_RS<-rbinom(R, p_RS)n_import_E<-rpois(epsilon)## Total population size, and number of infectedsI<- Ir+ IdN<- S+ E+ I+ R+ D## Initial statesinitial(S)<- S_iniinitial(E)<- E_iniinitial(Id)<-0initial(Ir)<-0initial(R)<-0initial(D)<-0## User defined parameters - default in parentheses:S_ini<-user(1000)# susceptiblesE_ini<-user(1)# infectedbeta<-user(0.3)# infection ratedelta<-user(0.3)# inverse incubationgamma_R<-user(0.08)# recovery rategamma_D<-user(0.12)# death ratemu<-user(0.7)# CFRomega<-user(0.01)# rate of waning immunityepsilon<-user(0.1)# import case rate

seirds_generator<- odin::odin(path_seirds_model)seirds_generator

## <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.seirds7393a829##     user: E_ini S_ini beta delta epsilon gamma_D gamma_R mu omega##     registration: function () ##     set_initial: function (step, y, use_dde) ##     update_metadata: function () ##   Parent env: <environment: namespace:discrete.stochastic.seirds7393a829>##   Locked objects: TRUE##   Locked class: FALSE##   Portable: TRUE

x<- seirds_generator$new()seirds_col<-c("#8c8cd9","#e67300","#d279a6","#ff4d4d","#999966","#660000")set.seed(1)x_res<- x$run(0:365)par(mar =c(4.1,5.1,0.5,0.5),las =1)matplot(x_res[,1], x_res[,-1],xlab ="Time",ylab ="Number of individuals",type ="l",col = seirds_col,lty =1)legend("left",lwd =1,col = seirds_col,legend =c("S","E","Ir","Id","R","D"),bty ="n")

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

x_res<-as.data.frame(replicate(100, x$run(0:365)[,-1]))dim(x_res)

## [1] 366 600

x_res[1:6,1:10]

##    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

seirds_col_transp<-paste0(seirds_col,"1A")par(mar =c(4.1,5.1,0.5,0.5),las =1)matplot(0:365, x_res,xlab ="Time",ylab ="Number of individuals",type ="l",col =rep(seirds_col_transp,100),lty =1)legend("right",lwd =1,col = seirds_col,legend =c("S","E","Ir","Id","R","D"),bty ="n")

100 replicates of a stochastic, discrete-timeSEIRDS model

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

check_model<-function(n =50,t =0:365,alpha =0.2, ...,legend_pos ="topright") {  model<- seirds_generator$new(...)  col<-paste0(seirds_col,"33")  res<-as.data.frame(replicate(n, model$run(t)[,-1]))  opar<-par(no.readonly =TRUE)on.exit(par(opar))par(mar =c(4.1,5.1,0.5,0.5),las =1)matplot(t, res,xlab ="Time",ylab ="",type ="l",col =rep(col, n),lty =1)mtext("Number of individuals",side =2,line =3.5,las =3,cex =1.2)legend(legend_pos,lwd =1,col = seirds_col,legend =c("S","E","Ir","Id","R","D"),bty ="n")}

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

check_model(beta =0,epsilon =0)

Stochastic SEIRDS model: sanity check with noinfections

Another easy case: no importation, no waning immunity:

check_model(epsilon =0,omega =0)

Stochastic SEIRDS model: no importation orwaning immunity

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

check_model(t =0:(365*3),epsilon =0.1,beta = .2,omega = .01,mu =0.005,S_ini =1e5)

Stochastic SEIRDS model: endemic state in alarger population