SEIR model equations

Our system of equations for an SEIR model is therefore:

\[\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*}\].

Population dynamics - the open epidemic

One simple change is to include vital dynamics - birth and death rates.

• Represented here by \(b\) and \(m\)
• People are born into the \(S\) compartments and die at an equal rate from each compartment.
• This is known as an open epidemic.

This is especially important over long time periods (years)

Vital dynamics equations

Our SEIR 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*}\]

Where \(b\) and \(m\) are the birth and mortality rates, and \(0<{b,m}<1\)

Note that \(m\) represents the general underlying mortality rate, not mortality specifically related to this epidemic.

A stable population

It is common practice to assume a stable population by setting \(m=b\)

• Births and deaths cancel out
• The population is stable at \(N\)

Note: In an open epidemic it is more common to work with proportions of people in compartments, rather than total numbers.

SEIR with vital dynamics - code

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
  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)
}

Example: the first 9 months (no vital dynamics)

The first nine months of this model without population dynamics (ie. with \(b=m=0\)).

Example: the first 9 months (with vital dynamics)

The same system with \(m=b = \frac{1}{7*365}\). This is for exaggeration purposes - it assumes an average life expectancy of 7 years.

Example: vital dynamics extended

The same SIR model plotted for five years. \(S\) gradually increases, a threshold is reached and there is another small peak (around day 1500).

Often used to model diseases such as Measles (pre-vaccine) where there can be a low level of the disease present in the population, with occasional outbreaks.

Example: vital dynamics extended even longer

The \(I\) compartment for a measles-like disease plotted for 70 years.

Here, \(\beta = \frac{13}{7},\,\sigma=\frac{1}{8},\,\gamma = \frac{1}{7}\) and \(m = b = \frac{1}{365\times 70}\)

Waning immunity

So far we have focussed on infectious diseases that lead to full immunity following infection.

In the SIR model, once people are in the \(R\) (Removed) compartment, they stay there forever.

In fact, with many infectious diseases:

• there is temporary immunity following infection
• after some period of time people are susceptible once again to the same infection.

This is an SEIRS model.

Waning immunity: equations

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).

Waning immunity: R code

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)
}

Example: Waning immunity

If we set \(w=0\) then people don’t leave the \(R\) compartment and immunity doesn’t wane:

Example: Waning immunity

Now we set \(w=0.01\) (so that people retain immunity for an average of 100 days).

Example: Waning immunity - long term

Plotting over 5 years, we see that an equilibrium is reached

Unlike the models we have seen so far, people cycle through the compartments.

Example: respiratory syncytial virus (RSV)

From White et al. (2007).

Human respiratory syncytial virus (RSV):

• One of the leading [respiratory] causes of hospitalisation for young children
• Once a child has been infected, they can be reinfected (even into adulthood).
• After being infected once, the infection can be milder
• For some time after being infected, people are less likely to be infected

The point of this example is to show you a real compartmental infectious disease model, and some of its features - don’t worry about retaining the detail!

Example: RSV model

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)

The birth and mortality rates are assumed equal (\(\mu\)) so that \[ S + I_S + I_R + R = N.\]

Example: RSV diagram

Example: RSV equations

There is a term for every transfer in the model:

\[\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*}\]

Notice that infection involves a function \(\Lambda\left(t\right)\)

Example: RSV varying transmissibility

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)}. \end{equation}\]

In White et al. (2007), \(Q\left(t\right)\) is used to represent seasonal variation (more on that later).

For now we set \(Q\left(t\right)=\frac{b}{\tau+\mu}\), so that

\[ \Lambda\left(t\right) = b\left(I_S + \eta I_R\right). \] Here,

• \(b\) represents ‘overall infectiousness’
• \(\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\).
  • If \(\eta = 0\) then people in \(I_R\) aren’t infectious.
  • If \(\eta = 1\) then people in \(I_R\) and \(I_S\) are equally infectious.

Example: running the RSV model

To run this model, we need

• initial values for the four compartments,
• values for each of the parameters \(\mu, \alpha, \tau, \rho, \sigma, b\) and \(\eta\).

We set initial proportions to:

\[\begin{align*} S & = 0.999\\ I_S& = 0.001\\ I_R& = 0,\\ R & = 0. \end{align*}\]

and parameter values to:

Example: RSV model outputs

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.

Example: Plotting Q

Exercise

Remember that the seasonality function \(Q\left(t\right)\) is involved in the model via the transmission function \(\Lambda\left(t\right)\),

In groups, discuss what you expect to see as we include seasonality in the model. How do you think it will differ from the version with no seasonality?

Example: RSV model output

After the initial outbreak, the system settles into a cycle.

For this plot, we set a=0.5, phi=0.