Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
CC-BY_NC

Produced with R version 3.2.3 and pomp version 1.3.1.2.

Objectives

This tutorial develops a general class of dynamic models with particular relevance for biological systems. We have the following goals:

  1. Dynamic systems can often be represented in terms of flows between compartments. We will develop the concept of a compartment model for which we specify rates for the flows between compartments.

  2. We develop deterministic and stochastic interpretations of a compartment model.

  3. We introduce Euler’s method to simulate from dynamic models, and we apply it to both deterministic and stochastic compartment models.


12.1 Compartment models


12.2 Compartment models in epidemiology: the SIR model and its generalizations

We will develop deterministic and stochastic representations of a susceptible-infected-recovered (SIR) system, a fundamental class of models for disease transmission dynamics. We will do this using notation which generalizes to more complex systems (Bretó et al. 2009).

\(S = \text{susceptible}\)
\(I = \text{infected and infectious}\)
\(R = \text{recovered and/or removed}\)


12.3 The ordinary differential equation (ODE) interpretation of the SIR model

Together with initial conditions specifying \(S(0)\), \(I(0)\) and \(R(0)\), we just need to write down ODEs for the flow counting processes. These are, \[ dN_{SI}/dt = \mu_{SI}(t) \, S(t),\] \[ dN_{IR}/dt = \mu_{IR}\, I(t).\]


12.4 The simple continuous-time Markov chain interpretation of the SIR model

\[\begin{array}{lcl} {\mathbb{P}}\big[N_{SI}(t+\delta){{=\,}}N_{SI}(t)+1\big] &=& \mu_{SI}(t) \, S(t) \delta + o(\delta) \\ {\mathbb{P}}\big[N_{SI}(t+\delta){{=\,}}N_{SI}(t)\big] &=& 1-\mu_{SI}(t) \, S(t) \delta + o(\delta) \\ {\mathbb{P}}\big[N_{IR}(t+\delta){{=\,}}N_{IR}(t)+1\big] &=& \mu_{IR} \, I(t) \delta + o(\delta) \\ {\mathbb{P}}\big[N_{IR}(t+\delta){{=\,}}N_{IR}(t)\big] &=& 1-\mu_{IR}(t) \, I(t) \delta + o(\delta) \end{array}\]


12.5 Exercises

12.5.1 From Markov chain to ODE

Find the expected value of \(N_{SI}(t+\delta)-N_{SI}(t)\) and \(N_{IR}(t+\delta)-N_{IR}(t)\) given the current state, \(S(t)\), \(I(t)\) and \(R(t)\). Take the limit as \(\delta\to 0\) and show that this gives the ODE model.

  • A simple counting process is one which cannot count more than one event at a time (Wikipedia: Point_process). Thus, in a technical sense, the SIR Markov chain model we have written is simple. One may want to model the extra randomness resulting from multiple simultaneous events: someone sneezing in a bus; large gatherings at football matches; etc. This extra randomness may even be critical to match the variability in data.

  • Later in the course, we may see situations where this extra randomness plays an important role. Setting up the model using counting processes, as we have done here, turns out to be useful for this.


12.6 Euler’s method for ordinary differential equations (ODEs)

  1. Euler’s method is the simplest (following the KISS principle).

  2. Euler’s method extends naturally to stochastic models, both continuous-time Markov chains models and stochastic differential equation (SDE) models.

  3. Close approximation of the numerical solutions to a continuous-time model is less important than it may at first appear, a topic worth further discussion…


12.7 Some comments on using continuous-time models and discretized approximations


12.8 Euler’s method for a discrete SIR model

\[\begin{array}{lcl} {\mathbb{P}}\big[N_{SI}(t+\delta){{=\,}}N_{SI}(t)+1\big] &=& \mu_{SI}(t) \, S(t) \delta + o(\delta), \\ {\mathbb{P}}\big[N_{IR}(t+\delta){{=\,}}N_{IR}(t)+1\big] &=& \mu_{IR} \, I(t) \delta + o(\delta). \end{array}\]


  1. A Poisson approximation. \[\tilde N_{SI}(t+\delta)= \tilde N_{SI}(t) + \mathrm{Poisson}\big[\mu_{SI}\big(\tilde I(t)\big) \, \tilde S(t) \,\delta\big],\] where \(\mathrm{Poisson}(\mu)\) is a Poisson random variable with mean \(\mu\) and \[\mu_{SI}\big(\tilde I(t)\big) = \beta\, \tilde I(t).\]


  1. A binomial approximation with transition probabilities approximated by rate times time. \[\tilde N_{SI}(t+\delta) = \tilde N_{SI}(t) + \mathrm{Binomial}\big[\tilde S(t),\mu_{SI}\big(\tilde I(t)\big) \, \delta\big),\] where \(\mathrm{Binomial}(n,p)\) is a binomial random variable with mean \(np\) and variance \(np(1-p)\).


  1. A binomial approximation with exponential transition probabilities.

\[ \tilde N_{SI}(t+\delta)= \tilde N_{SI}(t) + \mathrm{Binomial}\big[\tilde S(t),1-\exp\big\{-\mu_{SI}\big(\tilde I(t)\big) \delta\big\}\big].\]



12.9 Exercises

12.9.1 Theoretical exercise: Compartment models via stochastic differential equations (SDEs)

The Euler method extends naturally to stochastic differential equations. A natural way to add stochastic variation to an ODE \(dx/dt=h(x)\) is \[ dX/dt = h(X) + \sigma \, dB/dt\] where \(\{B(t)\}\) is Brownian motion and so \(dB/dt\) is Brownian noise. Then, an Euler approximation \(\tilde X(t)\) is generated by \[ \tilde X\big( \,(k+1)\delta\,\big) = \tilde X( k\delta) + \delta\, h\big(\, \tilde X(k\delta)\,\big) + \sigma \sqrt{\delta} \, Z_k\] where \(Z_1,Z_2,\dots\) is a sequence of independent standard normal random variables, i.e., \(Z_k\sim N[0,1]\). Although SDEs are often considered an advanced topic in probability, the Euler approximation doesn’t demand much more than familiarity with the normal distribution.

Write down a stochastic Euler method for an SDE representation of the SIR model. Consider some difficulties that might arise with non-negativity constraints, and propose some practical way one might deal with that issue.

  • A useful method to deal with positivity constraints is to use Gamma noise rather than Brownian noise (Bhadra et al. 2011,Laneri et al. (2010)). SDEs driven by Gamma noise can be investigated by Euler solutions simply by replacing the Gaussian noise by an appropriate Gamma distribution.


12.9.2 Conceptual exercise: Euler’s method vs Gillspie’s algorithm

  • A widely used, exact simulation method for continuous time Markov chains is Gillspie’s algorithm. We do not put much emphasis on Gillespie’s algorithm here. Why? When would you prefer an implementation of Gillespie’s algorithm to an Euler solution?

  • Numerically, Gillespie’s algorithm is often approximated using so-called tau-leaping methods. These are closely related to Euler’s approach. Is it reasonable to call a suitable Euler approach a tau-leaping method?


12.10 Compartmental models in pomp.

As an example that we can probe in some depth, let’s look at an isolated outbreak of influenza that occurred in a boarding school for boys in England (Anonymous 1978). Let’s examine the data:

bsflu <- read.table("bsflu_data.txt")
head(bsflu)
##              B  C day
## 1978-01-22   1  0   1
## 1978-01-23   6  0   2
## 1978-01-24  26  0   3
## 1978-01-25  73  1   4
## 1978-01-26 222  8   5
## 1978-01-27 293 16   6

The variable B refers to boys confined to bed and C to boys in convalescence. Let’s restrict our attention for the moment to the B variable.

bsflu <- subset(bsflu,select=c(day,B))
ggplot(data=bsflu,aes(x=day,y=B))+geom_line()+geom_point()

Let’s assume that \(B\) indicates the number of boys confined to bed the preceding day and that the disease follows the simple SIR model. Our tasks will be, first, to estimate the parameters of the SIR and, second, to decide whether or not the SIR model is an adequate description of these data.

Below is a diagram of the SIR model. The host population is divided into three classes according to their infection status: S, susceptible hosts; I, infected (and infectious) hosts; R, recovered and immune hosts. The rate at which individuals move from S to I is the force of infection, \(\lambda=\beta\,I/N\), while that at which individuals move into the R class is \(\gamma\).

Let’s look at how we can view the SIR as a POMP model. The unobserved state variables, in this case, are the numbers of individuals, \(S\), \(I\), \(R\) in the S, I, and R compartments, respectively. It’s reasonable in this case to view the population size \(N=S+I+R\), as fixed. The numbers that actually move from one compartment to another over any particular time interval are modeled as stochastic processes. In this case, we’ll assume that the stochasticity is purely demographic, i.e., that each individual in a compartment at any given time faces the same risk of exiting the compartment.

To implement the model in pomp, the first thing we need is a stochastic simulator for the unobserved state process. We’ve seen that there are several ways of approximating the process just described for numerical purposes. An attractive option here is to model the number moving from one compartment to the next over a very short time interval as a binomial random variable. In particular, we model the number, \({\Delta}{N_{SI}}\), moving from S to I over interval \({\Delta}{t}\) as \[{\Delta}{N_{SI}} \sim {\mathrm{Binomial}\left(S,1-e^{-\lambda{\Delta}{t}}\right)},\] and the number moving from I to R as \[{\Delta}{N_{IR}} \sim {\mathrm{Binomial}\left(I,1-e^{-\gamma{\Delta}{t}}\right)}.\]

A Csnippet that encodes such a simulator is as follows:

sir_step <- Csnippet("
  double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
  double dN_IR = rbinom(I,1-exp(-gamma*dt));
  S -= dN_SI;
  I += dN_SI - dN_IR;
  R += dN_IR;
")

At day zero, we’ll assume that \(I=1\) and \(R=0\), but we don’t know how big the school is, so we treat \(N\) as a parameter to be estimated and let \(S(0)=N-1\). Thus an initializer Csnippet is

sir_init <- Csnippet("
  S = N-1;
  I = 1;
  R = 0;
")

We fold these Csnippets, with the data, into a pomp object thus:

pomp(bsflu,time="day",t0=0,rprocess=euler.sim(sir_step,delta.t=1/6),
     initializer=sir_init,paramnames=c("N","Beta","gamma"),
     statenames=c("S","I","R")) -> sir

Now let’s assume that the case reports, \(B\), result from a process by which new infections result in confinement with probability \(\rho\), which we can think of as the probability that an infection is severe enough to be noticed by the school authorities. Since confined cases have, presumably, a much lower transmission rate, let’s treat \(B\) as being a count of the number of boys who have moved from I to R over the course of the past day. We need a variable to track this. Let’s modify our Csnippet above, adding a variable \(H\) to track the incidence. We’ll then replace the rprocess with the new one.

sir_step <- Csnippet("
  double dN_SI = rbinom(S,1-exp(-Beta*I/N*dt));
  double dN_IR = rbinom(I,1-exp(-gamma*dt));
  S -= dN_SI;
  I += dN_SI - dN_IR;
  R += dN_IR;
  H += dN_IR;
")

sir_init <- Csnippet("
  S = N-1;
  I = 1;
  R = 0;
  H = 0;
")

pomp(sir,rprocess=euler.sim(sir_step,delta.t=1/6),initializer=sir_init,
     paramnames=c("Beta","gamma","N"),statenames=c("S","I","R","H")) -> sir

Now, we’ll model the data, \(B\), as a binomial process, \[B_t \sim {\mathrm{Binomial}\left(H(t)-H(t-1),\rho\right)}.\] But we have a problem, since at time \(t\), the variable H we’ve defined will contain \(H(t)\), not \(H(t)-H(t-1)\). We can overcome this by telling pomp that we want H to be set to zero immediately following each observation. We do this by setting the zeronames argument to pomp:

pomp(sir,zeronames="H") -> sir

Now, to include the observations in the model, we must write both a dmeasure and an rmeasure component:

dmeas <- Csnippet("lik = dbinom(B,H,rho,give_log);")
rmeas <- Csnippet("B = rbinom(H,rho);")

and put these into our pomp object:

sir <- pomp(sir,rmeasure=rmeas,dmeasure=dmeas,statenames="H",paramnames="rho")

Let’s perform some simulations to verify that things are working. To do so, we’ll need some parameters. A little thought will get us some ballpark estimates. In the data, it looks like there were a total of 1540 infections, so the population size, \(N\), must be somewhat in excess of this number. In fact, we can use the final-size equation \[R_0 = -\frac{\log{(1-f)}}{f},\] where \(f=R(\infty)/N\) is the final size of the epidemic, together with the idea that \(R_0\) must be, say, around 1.5, to estimate that \(f\approx 0.6\), whence \(N\approx 2600\). If the infectious period is roughly 1 da, then \(1/\gamma \approx 1~\text{da}\) and \(\beta = \gamma\,R_0 \approx 1.5~\text{da}^{-1}\).

sims <- simulate(sir,params=c(Beta=1.5,gamma=1,rho=0.9,N=2600),
                 nsim=20,as=TRUE,include=TRUE)

ggplot(sims,mapping=aes(x=time,y=B,group=sim,color=sim=="data"))+
  geom_line()+guides(color=FALSE)


12.11 Exercises

12.11.1 Explore the SIR model

Fiddle with the parameters to see if you can’t find parameters for which the data are a more plausible realization.


12.11.2 The SEIR model

Below is a diagram of the so-called SEIR model. This differs from the SIR model in that infected individuals must pass a period of latency before becoming infectious.

Modify the codes above to construct a pomp object containing the flu data and an SEIR model. Perform simulations as above and adjust parameters to get a sense of whether improvement is possible by including a latent period.


12.11.3 Rethinking the boarding-school flu data

In the preceding, we’ve been assuming that \(B_t\) represents the number of boys sent to bed on day \(t\). Actually, this isn’t correct at all. As described in the report (Anonymous 1978), \(B_t\) represents the total number of boys in bed on day \(t\). Since boys were potentially confined for more than one day, the data count each infection multiple times. On the other hand, we have information about the total number of boys at risk and the total number who were infected. In fact, we know from Anonymous (1978) that \(N=763\) boys were at risk and \(512\) boys in total spent between 3 and 7 days away from class (either in bed or convalescent). Moreover, there is information in the data on the number of boys, \(C_t\), convalescent at day \(t\). Since \(1540~\text{boy-da}/512~\text{boy} \approx 3~\text{da}\), we know that the average duration spent in bed was 3 da and, since \(\sum_t\!C_t=0\), we can infer that the average time spent convalescing was \(0~\text{boy-da}/512~\text{boy} \approx 0~\text{da}\).

require(reshape2)
ggplot(data=melt(bsflu,id="day"),mapping=aes(x=day,y=value,color=variable))+
  geom_line()+geom_point()

Formulate a model with both confinement and convalescent stages. Implement it in pomp using a compartmental model like that diagrammed below.

You will have to give some thought to just how to model the relationship between the data (\(B\) and \(C\)) and the state variables. How many parameters can reasonably be fixed? How many must be estimated? Obtain some ballpark estimates of the parameters and simulate to see if you can plausibly explain the data as a realization of this model.


Acknowledgment

These notes draw on material developed for a short course on Simulation-based Inference for Epidemiological Dynamics by Aaron King and Edward Ionides, taught at the University of Washington Summer Institute in Statistics and Modeling in Infectious Diseases, 2015.

References

Anonymous. 1978. Influenza in a boarding school. British Medical Journal 1:587.

Bhadra, A., E. L. Ionides, K. Laneri, M. Pascual, M. Bouma, and R. Dhiman. 2011. Malaria in Northwest India: Data analysis via partially observed stochastic differential equation models driven by Lévy noise. Journal of the Americal Statistical Association 106:440–451.

Bretó, C., D. He, E. L. Ionides, and A. A. King. 2009. Time series analysis via mechanistic models. Annals of Applied Statistics 3:319–348.

Keeling, M., and P. Rohani. 2007. Modeling infectious diseases in humans and animals. Princeton University Press, Princeton.

Laneri, K., A. Bhadra, E. L. Ionides, M. Bouma, R. Yadav, R. Dhiman, and M. Pascual. 2010. Forcing versus feedback: Epidemic malaria and monsoon rains in NW India. PLoS Comput. Biol. 6:e1000898.