Lecture 4 - Vital Dynamics and extending models
1 Overview
So far we have explored the relatively simple closed SIR and SEIR models. However, real world infectious diseases are often much more complex.
In this lecture, we look at some examples of extending the SIR model to incorporate more realistic features:
- Adding an ‘E’ compartment (recap from workshop)
- Vital dynamics (births and deaths)
- Waning immunity
and go on to look at an example from the literature where a more complex infectious disease scenario is modelled.
1.1 Adding an \(E\) compartment
Let’s now investigate how we can add a pre-infectious (or exposed) \(E\) compartment, to create an SEIR model.
Figure 1.1: The trajectory in an SEIR model.
Remember that the exposed compartment contains people in the latent period. These people have been infected, but are not yet able to infect others (they are not yet infectious).
This additional compartment means that we need to rethink how people progress through the model, and our system of equations. Rather than going straight from \(S\) to \(I\), people now spend some time in \(E\).
1.1.1 Rate of change in susceptibles (\(S\))
When people leave the \(S\) compartment, they are now going to the \(E\) compartment. However, they are still being infected by people in the \(I\) compartment, so the equation for \(S\) remains the same:
\[\frac{dS}{dt} = -\frac{\beta I S}{N}.\]
1.1.2 Rate of change in pre-infectious (\(E\))
The people leaving \(S\) above are joining \(E\), so the first term will be \(\frac{\beta I S}{N}\).
As people’s latent period comes to an end, they move to the \(I\) compartment. We will model this with a rate, as with recovery. We will use \(\sigma\) for this rate parameter as the proportion of pre-infectious (exposed) people who become infectious at each time step.
Therefore we have:
\[\frac{dE}{dt} = \frac{\beta I S}{N} - \sigma E.\]
1.1.3 Rate of change in infectious (\(I\))
People are now joining the \(I\) compartment from \(E\) with rate \(\sigma\). This leads to a \(\sigma E\) term. We also still have people leaving \(I\) with recovery rate \(\gamma\), to go to \(R\). Therefore:
\[\frac{dI}{dt} = \sigma E - \gamma I. \]
The rate of change in recovered/removed (\(R\)) is unchanged, since people still recover from the \(I\) compartment.
Our system of equations for an SEIR model is therefore
$$ \[\begin{align*} \frac{dS}{dt} & = -\frac{\beta I S}{N}\\ \frac{dE}{dt} & = \frac{\beta I S}{N} - \sigma E\\ \frac{dI}{dt} & = \sigma E - \gamma I.\\ \frac{dR}{dt} & = \gamma I. \end{align*}\] $$
2 Population dynamics - the open epidemic
[Fig 4.17, p84 - rewrite a bit!]
One simple change we can make, particularly if we are studying the behaviour of a disease long term, is to include vital dynamics - birth and death rates, represented here by \(b\) and \(m\) respectively. Because people are now entering and leaving the system, this is known as an open epidemic.
Figure 2.1: The trajectory in an SEIR model with vital dynamics.
Our equations now become
\[ \begin{align*} \frac{dS}{dt} &= \color{green}{b N}- \color{purple}{m S} -\frac{\beta I S}{N}\\ \frac{dE}{dt} &= \frac{\beta I S}{N} - \sigma E - \color{purple}{m E} \\ \frac{dI}{dt} &= \sigma E - \gamma I - \color{purple}{m I}\\ \frac{dR}{dt} &= \gamma I - \color{purple}{m R} \end{align*}. \]
People are born into the Susceptible compartment at a rate of \(b\), represented by the green \(b N\). People from within each category die at rate \(m\), represented by the purple terms. Note that this term \(m\) represents the general underlying mortality rate, not mortality specifically related to this epidemic.
It is common practice (though perhaps not entirely realistic) to assume a stable population, so that the births and deaths cancel out, and therefore \(m=b\).
In the long term, as more people are born into the \(S\) compartment, new waves of the disease can occur.
Note In an open epidemic, where the size of the population may fluctuate, it is more common to work with proportions of people in compartments, rather than total numbers.
2.1 Example: the open epidemic
To illustrate this, let’s start with parameter values \[ \begin{align*} \beta &= \frac{2}{7} \\ \sigma & = \frac{1}{2} \\ \gamma & = \frac{1}{14} \end{align*} \] and initial values (now in terms of proportions) \[ \begin{align*} S & = 0.999\\ E & = 0 \\ I& = 0.001\\ R & = 0 \end{align*} \] The values for \(\beta\) and \(\gamma\) suggest a long infectious period, but that not much contact is being made by infectious people (or that the probability of transmission is low).
We code this model as follows
library(deSolve)
sir_vital = function (t, y, params)
{
# label state variables from within the y vector
S = y[1] # susceptibles
E = y[2] # Exposed
I = y[3] # infectious
R = y[4] # recovered
N = S + E + I + R
# Pull parameter values from the params vector (including new m and b parameters)
beta = params["beta"]
gamma = params["gamma"]
sigma = params["sigma"]
m = params["m"]
b=params["b"]
# Define equations, now with m and b terms
dS = b*N -m*S-(beta*S*I)/N
dE = beta*S*I/N - sigma*E - m*E
dI = sigma*E - gamma*I -m*I
dR = gamma * I - m*R
# return list of gradients
res = c(dS, dE, dI, dR)
list(res)
}Figure 2.2 shows the first nine months of this model without population dynamics (ie. with \(b=m=0\)).
Figure 2.2: An SEIR model without vital dynamics shown for 9 months.
Figure 2.3 shows the same system with \(m=b = \frac{1}{7*365}\). This is for exaggeration purposes - it assumes an average life expectancy of 7 years.
# set values for parameters, initial values and model output timepoints
params_list = c(
beta = 2/7,
gamma = 1 / 14,
sigma = 1/2,
m = 1/(365*7),
b=1/(365*7)
)
# Initial values for differential equations
initial_values = c (S = 0.999, E=0, I = 0.001, R = 0)
# Time points for model output (39 weeks in increments of 0.1 days)
timepoints = seq (0, 39*7, by=0.1)
# solve ODEs
output_vital = ode (y=initial_values, times = timepoints, func = sir_vital, parms= params_list)
# coerce to a data.frame object
outputv_df = as.data.frame(output_vital)
# and then to longer format
# Here, all population values will be collected into a new column 'value',
# and there will be a corresponding 'Compartment' column
outputv_tidy = pivot_longer(outputv_df, cols = c("S", "E", "I", "R"), names_to = "Compartment")
# Plot compartent populations against time, with a separate colour for each compartment
ggplot(data=outputv_tidy, aes(x=time, y=value, col=Compartment)) + geom_line()+
ylab("People") + xlab("Time (days)")
Figure 2.3: An SIR model with vital dynamics shown for 9 months.
Notice that in Figure 2.2, the compartment populations remain constant once the epidemic has died down; almost everyone has had the disease and recovered, and there will be no further change. By contrast, in Figure 2.3, the population within each compartment continues to change, as people are born into the susceptible compartment, and die from each compartment. This is most noticeable by the \(S\) line gradually increasing, and the \(R\) line gradually decreasing.
Figure 2.4 shows the same SEIR model as in Figure 2.3, but now plotted for five years. As the population of susceptibles gradually increases, a threshold is reached and there is another small wave of infections (around day 1500).
Figure 2.4: An SEIR model with vital dynamics shown for 5 years.
This sort of model is often used to model diseases such as Measles (pre-vaccine) where there can be a low level of the disease present in the population. Figure 2.5 shows the infectious compartment for a measles-like disease over 70 years, this time with a population life expectancy of 70 years.
Figure 2.5: The infectious compartment for a measles-like disease over 70 years.
Here, \(\beta = \frac{13}{7},\,\sigma=\frac{1}{8},\,\gamma = \frac{1}{7}\) and \(m = b = \frac{1}{365\times 70}\). We see that at first there are large outbreaks with very little occurrence between, but after a number of years there is much less fluctuation, and always some infection present in the population.
3 Waning immunity
So far we have focussed on infectious diseases that lead to full immunity following infection (at least enough to assume that in a model). In SIR model terms, once people are in the R (Removed) compartment, they stay there forever. In fact, with many infectious diseases, there is a temporary immunity following infection, but then after some period of time people are susceptible once again to the same infection.
Figure 3.1: The trajectory in an SEIR model with waning immunity.
We use the parameter \(w\) to denote the rate at which immunity is lost, and our equations become
\[ \begin{align*} \frac{dS}{dt} &= \color{blue}{wR} -\frac{\beta I S}{N}\\ \frac{dE}{dt} &= \frac{\beta I S}{N} - \sigma E \\ \frac{dI}{dt} &= \sigma E - \gamma I \\ \frac{dR}{dt} &= \gamma I - \color{blue}{wR} \end{align*}. \]
The average time for which someone retains their immunity will therefore be \(\frac{1}{w}\) (in terms of the model’s time units).
To solve this model using ode we will need to define a new function
library(deSolve)
sir_wane = function (t, y, params)
{
# label state variables from within the y vector
S = y[1] # susceptibles
E = y[2] # Exposed
I = y[3] # infectious
R = y[4] # recovered
N = S + E + I + R
# Pull parameter values from the params vector, including new 'w' parameter
beta = params["beta"]
gamma = params["gamma"]
sigma = params["sigma"]
w = params["w"]
# Define equations, now with w terms
dS = w*R - (beta*S*I)/N
dE = beta*S*I/N - sigma*E
dI = sigma*E - gamma*I
dR = gamma * I - w*R
# return list of gradients
res = c(dS, dE, dI, dR)
list(res)
}3.1 Example: waning immunity
To illustrate the effect of waning immunity, we will use the same parameter and initial values as in our open epidemic example.
\[ \begin{align*} \beta &= \frac{2}{7} \\ \sigma & = \frac{1}{2} \\ \gamma & = \frac{1}{14} \end{align*} \] and initial values (now in terms of proportions) \[ \begin{align*} S & = 0.999\\ E & = 0 \\ I& = 0.001\\ R & = 0 \end{align*} \] If we set \(w=0\) then people don’t leave the \(R\) compartment and immunity doesn’t wane, as in Figure 3.2 (and Figure 2.2).
Figure 3.2: An SEIR model without waning immunity shown for 9 months.
In Figures 3.3 and 3.4 we have set \(w=0.01\) (meaning that people retain immunity for an average of 100 days).
Figure 3.3: An SIR model with waning immunity shown for 9 months.
Figure 3.4: An SEIR model with waning immunity shown for 5 years.
Notice that an equilibrium is still reached, but with the number of people in the \(R\) compartment much fewer, and a constant proportion of people in the exposed and infectious compartments. Note that although the proportions remain in equilibrium, people are (notionally) cycling through the compartments, becoming infectious, recovering, returning to the susceptible compartment and becoming re-infected, and so on. This is in constrast to the example with no waning immunity in Figure 3.2, where the same people remain in the susceptible and removed compartments forever. It is also in constrast to the open epidemic in Figure 2.3, where any particular person is only infected once.
4 Example: respiratory syncytial virus (RSV)
The following example comes from White et al. (2007), and demonstrates some of the complexity that can be represented with a compartmental model structure.
Human respiratory syncytial virus (RSV) is one of the leading causes of hospitalisation for young children who have a respiratory infection. Once a child has been infected, a child (or subsequently adult) can be reinfected.
A simplified version of the model developed by White et al. (2007) is shown in the equations below and in Figure 4.1. In this model, there are four compartments:
- \(S\) - susceptible people, including those who have never yet been infected
- \(I_S\) - those infected for the first time (or from the \(S\) compartment)
- \(R\) - those recovered, but still susceptible to reinfection (but not as susceptible as those in \(S\))
- \(I_R\) - those infected from the \(R\) compartment (must have been infected at least once before)
In this model the birth and mortality rates are assumed equal (\(\mu\)) so that \[ S + I_S + I_R + R = N.\]
The equations governing the dynamics are as follows:
\[\begin{align*} \frac{dS}{dt} & = \mu N-\frac{\Lambda\left(t\right) S}{N} + \frac{\alpha \tau}{\rho}\left(I_S + I_R + R\right) - \mu S\\ \frac{dI_S}{dt} & = \frac{\Lambda\left(t\right) S}{N} - \left(\tau + \mu\right) I_S \\ \frac{dI_R}{dt} & = \frac{\sigma\Lambda\left(t\right) R}{N} - \left(\frac{\tau}{\rho} + \mu\right)I_R\\ \frac{dR}{dt} & = \left(1-\frac{\alpha}{\rho}\right)\tau I_S + \left(1-\alpha\right)\frac{\tau}{\rho}I_R - R\left(\frac{\sigma\Lambda\left(t\right)}{N}+\frac{\alpha\tau}{\rho}+\mu\right). \end{align*}\] These correspond to the diagram in Figure 4.1.
Figure 4.1: The RSV model from White et al. (2007).
This model contains many more parameters than we are used to, but here is an idea of what is happening:
People begin in compartment \(S\). At rate \(\Lambda\left(t\right)\) they are infected, and join the \(I_S\) compartment. Their infection is lost at rate \(\tau\) (so the average duration of infectiousness in \(I_S\) is \(\frac{1}{\tau}\)). Proportion \(\frac{\alpha}{\rho}\) return to the \(S\) (fully susceptible) compartment, while proportion \(1-\frac{\alpha}{\rho}\) join the \(R\) (recovered) compartment. Note that to make sense of these proportions, we should have \(0\leq{\alpha}\leq{\rho}\).
From the \(S\) compartment, people become reinfected at rate \(\Lambda\left(t\right)\), but in the \(R\) compartment the rate of reinfection is damped by the ‘partial susceptibility factor’ \(\sigma\) (where we assume \(0\leq{\sigma}\leq{1}\)), so that the rate of infection from \(R\) is \(\sigma\Lambda\left(t\right)\). When people in \(R\) are reinfected, they join compartment \(I_R\) (Re-infected),
The rate of people leaving \(I_R\) (to either \(R\) or \(S\)) is \(\frac{\tau}{\rho}\), where \(0<\rho\leq{1}\). This represents a reduced duration of infectiousness for those who have been infected from the \(R\) compartment, of \(\frac{\rho}{\tau}\).
To represent waning immunity, people move from \(R\) to \(S\) at rate \(\frac{\alpha\tau}{\rho}\). This means that setting \(\alpha=0\) would create a model where immunity did not wane.
The function \(\Lambda\left(t\right)\) (which governs infectiousness) is given by
\[\begin{equation} \Lambda\left(t\right) = Q\left(t\right)\times{}\left(I_S + \eta I_R\right)\times{\left(\tau+\mu\right)}. \tag{4.1} \end{equation}\]
In White et al. (2007), \(Q\left(t\right)\) is used to represent seasonal variation. We will ignore this for now and set \(Q\left(t\right)=\frac{b}{\tau+\mu}\), which leads to
\[ \Lambda\left(t\right) = b\left(I_S + \eta I_R\right). \] Here, \(b\) represents ‘overall infectiousness’ and \(\eta\) is a damping parameter (with \(0\leq{\eta}\leq{1}\)) so that people in \(I_R\) are less infectious than those in \(I_S\).
To run this model, we therefore need initial values for the four compartments, and values for each of the parameters \(\mu, \alpha, \tau, \rho, \sigma, b\) and \(\eta\). In Figure 4.2 we have set the parameter values to: \[ \begin{align*} \mu& = \frac{1}{365\times{70}}\\ \alpha& = 0.16\\ \tau& = \frac{1}{6}\\ \rho & = 0.8\\ \sigma & = 0.3\\ b & = 0.5\\ \eta & = 0.8 \end{align*} \] and the initial populations (as proportions) to: \[\begin{align*} S & = 0.999\\ I_S& = 0.001\\ I_R& = 0,\\ R & = 0. \end{align*}\] As Figure 4.2 shows, this system reaches an equilibrium by around day 60. Note that people are still moving between compartments, but that the proportions in each remain steady.
Figure 4.2: Compartment populations in the RSV model described above, with no seasonality.
4.1 Introducing seasonality
In the work above, we set \(Q\left(t\right)\) to be constant, removing the seasonal variation from the model. However, seasonal variability is an important component of many infectious disease models, and so we will now include it. First, let’s explore the function \(Q\left(t\right)\).
White et al. (2007) specify \[ Q\left(t\right) = \frac{b\left(a\cos\left[2\pi\left(t-\phi\right)\right]+1\right)}{\tau+\mu}. \] Here \(b, \tau\) and \(\mu\) represent the overall infectiousness, rate of loss of primary infection and birth/mortality rate respectively, as before, and \(t\) represents time (in years). The peak in transmission is given by \(0\leq{\phi}\leq{1}\) as a fraction of the year, so that \(\phi=0\) is \(1^{st}\) January and \(\phi=1\) is \(31^{st}\) December. The relative amplitude of \(Q\left(t\right)\) is given by \(0\leq{a}\leq{1}\), where \(a=0\) gives no seasonal variation, and \(a=1\) gives puts transmission at 0 at the minumum and twice the mean at the peak.
Figure 4.3: The function Q, with various values for phi and a.
Remember that the seasonality function \(Q\left(t\right)\) is involved in the model via the function \(\Lambda\left(t\right)\) from Equation (4.1).
Exercise 4.1 In groups, discuss what you expect to see as we include seasonality in the model.
Figure 4.4: Compartment populations in the RSV model described above, with seasonality. For this plot, we set a=0.5, phi=0.
We can see that after the initial wave, when \(I_S\) reaches its highest, the system settles into a cycle. The seasonal change in infectiousness stops the sytem from reaching a steady state as in Figure 4.2. You can read more about this model and how the parameter values were chosen to fit datasets from many different locations in White et al. (2007).
5 Summary
In this lecture we have looked at some developments to SIR-type models that enable them to model more complex behaviour.
- By including vital dynamics we can model diseases over very long periods of time
- By including waning immunity we can model diseases for which immunity is not permanent
- We can include nuances such as differing recovery rates and seasonal variations in transmission, as shown by our RSV example.