Introduction

This tutorial provides a mathematical and theoretical background for deterministic, compartmental models (DCMs), with instructions on how to run the built-in models designed for learning in EpiModel. For information on how to extend these models to simulate novel epidemiological processes, see the related tutorial New DCMs with EpiModel.

Deterministic compartmental models solve differential equations representing analytic epidemic systems in continuous time. The models are deterministic because their solutions are fixed mathematical functions of the input parameters and initial conditions, with no stochastic variability in the disease and demographic transition processes. The models are compartmental because they divide the population into groups representing discrete disease states (e.g., susceptible and infected), and further on demographic, biological, and behavioral traits that influence disease transmission. In contrast to the stochastic models presented below, individuals within the population are not discretely represented.

A Basic SI Model

Starting with a basic Susceptible-Infected (SI) disease model in which there is random mixing within the population, the size of each compartment over time is represented by the equations:

\[ \begin{align} \frac{dS}{dt} &= -\lambda S \notag \\[4pt] \frac{dI}{dt} &= \lambda S \end{align} \]

where \(\lambda\) is the force of infection and represents \(\frac{\beta c I}{N}\). \(\beta\) is the probability of transmission per contact, \(c\) is the rate of contact per person per unit time, \(I\) is the number infected at time \(t\) and and \(N\) is the population size at time \(t\) (we drop the \(t\) subscript as it is implicit throughout).

Because ``contact’’ has been defined several ways in the modeling literature, we use the word act to represent the action, such as face-to-face discussion or sexual intercourse, by which disease may be transmitted. The force of infection is multiplied by the current state sizes of the \(S\) to solve the differential equation yielding the rate of change for the compartment.

To simulate a deterministic model in EpiModel, use the dcm function. Prior to running the model, it is necessary to parameterize it. Model parameters, initial conditions, and controls settings are input into three helper functions that organize the model framework. In param.dcm, the epidemic model parameters are entered. The inf.prob argument sets the transmission probability per act, and act.rate sets the acts per person per unit time. The init.dcm function collects the initial conditions for the model, and since this is an SI model it is necessary to specify the initial number susceptible and infected at \(t_1\). The control.dcm finally collects other structural model controls like the model type and number of time steps for the simulation.

param <- param.dcm(inf.prob = 0.2, act.rate = 0.25)
init <- init.dcm(s.num = 500, i.num = 1)
control <- control.dcm(type = "SI", nsteps = 500)

The model parameters, initial conditions, and controls are then entered as inputs into the dcm function to run the model and save the output in our object mod.

mod <- dcm(param, init, control)

Options for analyzing the results are consistent across all the model classes and types. Printing the model object provides basic information on model input and output, including model parameters and data variables for plotting and analysis.

mod
EpiModel Simulation
=======================
Model class: dcm

Simulation Summary
-----------------------
Model type: SI
No. runs: 1
No. time steps: 500
No. groups: 1

Model Parameters
-----------------------
inf.prob = 0.2
act.rate = 0.25

Model Output
-----------------------
Variables: s.num i.num si.flow num

The output shows that two compartments and one flow are available. In EpiModel, regardless of the model class, compartments are discrete disease states and flows are the transitions between states; when demographic processes are introduced, flows represent transitions in and out of the population. The i.num compartment is the size of the infected population at each of the solved time steps. The endogenous disease flow names represent the starting and ending state: si.flow is the number of people moving from \(S\) to \(I\) at each time step. In epidemiological terms, i.num and si.flow are disease prevalence and incidence.

To plot the results of the model we use the generic plot function, which by default plots all the model compartment state sizes over time.

plot(mod)

After examining the plot, one can investigate the the size of each compartment at a specific time step. That is available with the summary function. At \(t_{150}\), 22.5% of the population have become infected, and disease incidence is 4.37 new infections.

summary(mod, at = 150)
EpiModel Summary
=======================
Model class: dcm

Simulation Summary
-----------------------
Model type: SI
No. runs: 1
No. time steps:
No. groups: 1

Model Statistics
------------------------------
Time: 150    Run: 1 
------------------------------ 
                n    pct
Suscept.  112.845  0.225
Infect.   388.155  0.775
Total     501.000  1.000
S -> I      4.311     NA
------------------------------

SIR Model with Demography

In a Susceptible-Infectious-Recovered (SIR) model, infected individuals transition from disease into a life-long recovered state in which they are never again susceptible. Here we model an SIR disease by adding to our basic SI model a recovery process. We also introduce demographic processes so that persons enter and exit the population through births and deaths. The model is represented by the following system of differential equations:

\[ \begin{align} \frac{dS}{dt} &= -\lambda S + bN - \mu_s S \notag \\[4pt] \frac{dI}{dt} &= \lambda S - \nu I - \mu_i I \notag \\[4pt] \frac{dR}{dt} &= \nu I - \mu_r R \end{align} \]

where \(b\) is the birth rate, \(\mu\) are the mortality rates specific for each compartment, and \(\nu\) is the recovery rate; note that across EpiModel, birth and mortality rates are more generally referred to as arrival and departure rates, respectively. In an SIR model, the recovery rate is the reciprocal of the average duration of disease; likewise, the reciprocal of the death rates are the average life expectancy for persons in those compartments.

Simulation

In EpiModel, introducing new transition processes into the model is straightforward. In param.dcm, parameters for the recovery rate, birth rate, and state-specific death rates are entered. These parameters imply that the birth rate is slightly higher than the underlying death rate among susceptibles, and that there is disease-induced mortality because the di.rate is larger than the other two death rates. In init.dcm, it is necessary to specify the number of initially recovered, even if that is 0. In control.dcm, the dt argument may be used to obtain model results in fractional time units (i.e., results are available for \(t_1\), \(t_{1.5}\), \(\dots\), \(t_{499.5}\), \(t_{500}\)).

param <- param.dcm(inf.prob = 0.2, act.rate = 1, rec.rate = 1/20,
                   a.rate = 1/95, ds.rate = 1/100, di.rate = 1/80, dr.rate = 1/100)
init <- init.dcm(s.num = 1000, i.num = 1, r.num = 0)
control <- control.dcm(type = "SIR", nsteps = 500, dt = 0.5)
mod <- dcm(param, init, control)

Plotting

Next we plot the results of the model to demonstrate several plot arguments. First, the par function is used to change some default graphical options. In the left plot, the popfrac=FALSE argument plots the compartment size (rather than prevalence) and alpha increases the transparency of the lines for better visualization. By default, the plot function will plot the prevalences for all compartments in the model, but in the right plot we override that using the y argument to specify that disease incidence (the si.flow element of the model object) should be plotted.

par(mar = c(3.2, 3, 2, 1), mgp = c(2, 1, 0), mfrow = c(1, 2))
plot(mod, popfrac = FALSE, alpha = 0.5,
     lwd = 4, main = "Compartment Sizes")
plot(mod, y = "si.flow", lwd = 4, col = "firebrick",
     main = "Disease Incidence", legend = "n")