Skip to main content

Introducing Mathematical Biology: The SIR model with demographics

Introducing Mathematical Biology
The SIR model with demographics
    • Notifications
    • Privacy
  • Project HomeNatural Sciences Collection: Anatomy, Biology, and Chemistry
  • Projects
  • Learn more about Manifold

Notes

Show the following:

  • Annotations
  • Resources
Search within:

Adjust appearance:

  • font
    Font style
  • color scheme
  • Margins
table of contents
  1. Cover
  2. Title Page
  3. Copyright
  4. Table Of Contents
  5. Introduction
  6. Population ecology
    1. Single population models
    2. Interacting populations 1: competition
    3. Interacting populations 2: predator-prey
  7. Infectious disease
    1. Epidemics in human populations
    2. The SIR model with demographics
    3. Diseases of ecological populations
    4. Evolution and adaptive dynamics
  8. Immune and cell dynamics
    1. A within-host Covid-19 model
    2. A within-host HIV-I model
    3. Introducing models of cancer dynamics
    4. A model of cancer volume dynamics
  9. Gene networks
    1. Introducing gene networks
    2. Autoregulation 1: auto-repression
    3. Autoregulation 2: auto-activation
    4. Longer negative feedback networks
    5. A two-gene toggle switch
  10. Pharmacokinetics
    1. Single intravenous bolus dose
    2. Repeated intravenous bolus doses
    3. Single and repeated oral doses
    4. A two compartment bolus model
  11. Background reviews
    1. Phase portraits
    2. Linear stability analysis
    3. Bifurcations
  12. Final thoughts and acknowledgements
  13. References

5

The SIR model with demographics

Including Births and Deaths

In our initial epidemic model in the previous chapter we had only two mechanisms involved – infection and recovery. Let us now add some increased detail into the model by including birth and death processes. Given the timescales of epidemic processes, this likely means looking at dynamics over a much longer time period, at least in human populations. Let’s begin by drawing a schematic of the system again:

Schematic of the SIR model with births and deaths now included. Three boxes are shown labelled susceptible, infected and recovered respectively. Starting from the left, an arrow points from blank space into the susceptible box, from the susceptible into the infected box and from the infected into the recovered box. In addition, an arrow leaves each box into blank space.
Schematic of the SIR model with births and deaths now included.

We will make the simplifying assumption that the birth rate and death rate are equal, so that the population would be at equilibrium in the absence of disease (this is a questionable assumption for many ecological populations, but perhaps not too unreasonable for modern human populations). If we assume all individuals produce offspring at rate \mu, that everyone is born susceptible, and that every individual has the same death rate, also \mu, this leads us to the equations,

\begin{align} &\frac{dS}{dt} = \mu N - \beta SI - \mu S\\ &\frac{dI}{dt} = \beta SI - (\gamma +\mu)I\\ &\frac{dR}{dt} = \gamma I-\mu R \end{align}

Once more, dN/dt=0 so we can eliminate R=N-(S+I) (again, this is somewhat by design – if the rates of birth and death were not equal this simplification could not be made). One thing to note here is that the infectious period has changed to be 1/(\gamma+\mu). When we define R_0, therefore, in our updated model we will have R_0=\beta N/(\gamma+\mu).

Endemic disease

At this point we could non-dimensionalise our system as we have seen previously. This would allow us to reduce the number of parameters in our model to make life easier, as well as revealing potentially useful information about the scales involved. You will see some studies do this and others don’t. For now, let us continue with our analysis without doing so, since it means we retain the clear biological meanings of all of our parameters and variables.

First we should find the equilibria of our system, where dS/dt=0 and dI/dt=0 simultaneously. Recall in the previous model this only happened when I=0, with this condition satisfying both ODEs and leaving us with a line of equilibria. Now, we find more ‘standard’ unique equilibria. There are two cases where dI/dt=0:

  • I^*=0, giving dS/dt=0\implies S^*=N;

This means the population is disease-free.

  • S^*=\dfrac{\gamma+\mu}{\beta}, giving dS/dt=0\implies I^*=\dfrac{\mu(N-S^*)}{\beta S^*}=\dfrac{\mu}{\beta}\left(\dfrac{N}{S^*}-1\right);

This means the disease is endemic.

We can do a bit more manipulation of that last equilibrium,

\begin{equation} I^*=\frac{\mu}{\beta}\left(\dfrac{N}{S^*}-1\right)=\frac{\mu}{\beta}\left(\frac{\beta N}{\gamma+\mu}-1\right)=\frac{\mu}{\beta}\left(R_0-1\right). \end{equation}

While both expressions are equally correct, the latter one is nice in that R_0 is included in the expression. We now have two important terms from epidemiology. An epidemic occurs when a disease initially spreads through a population. A disease is endemic when it remains at a steady level within the population.

Having found these equilibria we now want to classify their stability. Let’s look at the Jacobian for our system,

\begin{align} J=&\left( \begin{array}{cc} -\beta I^*-\mu & -\beta S^*\\ \beta I^* & \beta S^*-(\gamma+\mu) \end{array} \right) \end{align}

We shall now deal with each of our equilibria in turn.

Disease-free equilibrium

\begin{equation*} J=\left( \begin{array}{cc} -\mu & -\beta N\\ 0 & \beta N-(\gamma+\mu) \end{array} \right) \end{equation*}

Having a zero in an off-diagonal (for a 2×2 matrix) makes our life much easier, as we can just read off the eigenvalues as the two diagonal entries. In this case we therefore have \lambda_1=-\mu and \lambda_2=\beta N-(\gamma+\mu)=(\gamma+\mu)(R_0-1). Therefore, if R_0\lt1, the disease-free equilibrium is unstable.

Endemic equilibrium

\begin{equation*} J=\left( \begin{array}{cc} - \mu R_0& -(\gamma+\mu)\\ \mu(R_0-1)& 0 \end{array} \right) \end{equation*}

In this case we cannot read off our eigenvalues, so we shall instead look at the signs of the trace and determinant:

  • tr(J) = -\mu R_0\lt 0
  • \det(J) =\mu(\gamma+\mu)(R_0-1)

Since the trace is always negative, from the determinant we see that if R_0\gt 1, the endemic equilibrium is stable.

The stable equilibrium can be a node or a spiral depending on the precise parameter values. Generally, high \gamma, intermediate R_0 and low \mu make a spiral more likely.

Putting these results together we can see that when R_0\gt 1 the disease will become endemic.

Phase portraits

We can explore these two scenarios further with phase portraits. First we remember that again we must have S+I\le N, so we only have a triangular region in our phase portraits to really worry about. To find the nullclines:

\begin{align*} &dS/dt=0\implies I=\frac{\mu}{\beta}\left(\frac{N}{S}-1\right)\\ &dI/dt=0\implies I=0,S=\frac{\gamma+\mu}{\beta}. \end{align*}

Exercises

Draw the two qualitatively different phase portraits for this system.
Click for solution

The two phase portraits are sketched below. If you start sketching a phase portrait you should find that the only way you can produce two qualitatively different phase portraits is again about the placement of the vertical nullcline, now given by S=(\gamma+\mu)/\beta, which in turn relates to whether or not R_0\gt1 or not. In the first diagram, all trajectories will eventually approach the bottom-right corner at the disease-free equilibrium. In the second diagram, you should see the new endemic equilibrium in the middle where the two nullclines cross. Trajectories approach this equilibrium, moving around it anti-clockwise (as I once read, if you can’t remember clockwise/anti-clockwise directions, clockwise is the way you dial a phone).

Two hand-drawn phase portraits. In both cases population densities of I are on the vertical axis with no numeric values, and population densities of S are on the horizontal axis with no numeric values. In each case a line runs from I=N to S=N, creating a triangular region with the axes. In both cases there is a vertical line at S=gamma+mu/beta and a curve coming from high values of I at S=0, entering the triangular region then saturating slightly and crossing through I=0 at S=N.1. The vertical line is outside the triangular region to the right. Trajectories above the blue line move south-west and those below it move south-east. 2. The vertical line is inside the traingular region. Two trajectories are shown both spiralling anti-clockwise into the central equilibrium.
Two phase portraits for the SIR model. Blue lines give the susceptible nullclines, green lines the infected nullclines, and the black line the boundary of biological feasibility, with trajectories in pink.

We see nice demonstrations of the behaviour we predicted from our stability analysis here. In particular, when R_0\gt 1 the disease will persist in the long-term.

R_0 is a quantity that is often estimated for diseases when they emerge, and gives a good indication of how contagious they are and (as we shall see below) how easy they will be to eradicate. We can also use R_0 as a useful bifurcation parameter to draw a bifurcation diagram of this system below. Note again the key value of R_0=1 where the transcritical bifurcation occurs.

image0 (and a small area of the dashed line is shown below I=0)." width="694" height="421">
Bifurcation diagram of the SIR model, plotting the equilibria of infected densities in terms of R_0.

Vaccination

Suppose we want to prevent a disease spreading through a vaccination programme. What proportion, p, of the population do we need to vaccinate to succeed? If our vaccine works perfectly then the pre-infection density of susceptibles reduces from N to N(1-p). Recall that for the disease to die out we want R_0=\beta N/(\gamma+\mu)\lt 1. After vaccination this becomes,

\begin{align*} \frac{\beta N(1-p)}{\gamma+\mu} \implies p\gt 1-\frac{1}{R_0} \end{align*}

Note that not all of the population needs to be vaccinated – just enough to prevent the disease from spreading freely. This is important as there may be certain groups, those undergoing medical treatment, for example, who cannot be vaccinated. This population-wide protection created by significant vaccination is called herd immunity. The higher the R_0 of the disease, the greater proportion of the population that needs to be vaccinated.

Exercises

A novel infectious disease has been detected in a small town with a population of 10,000 people. The yearly birth and death rates of the population are estimated to be 0.02. From previous outbreaks in other towns the estimated transmission coefficient is \beta=7.5\times10^{-4} and the (yearly) recovery rate is 4.48. A vaccine exists that will perfectly prevent infections from occurring. What proportion of the population must be vaccinated to achieve herd immunity?
Click for solution

Using the values given we can find the value for R_0 as,

\begin{equation} R_0=\frac{\beta N}{\gamma+\mu}=\frac{7.5}{4.5}=\frac{5}{3}. \end{equation}

Since the vaccine prevents all infections we can use our herd immunity threshold as above to give,

\begin{align*} p\gt &1-\frac{1}{R_0}\\ \gt& 1 - \frac{3}{5}. \end{align*}

We therefore need to vaccinate at least 40% of the population to achieve herd immunity.

A point about herd immunity

This method to calculate the herd immunity threshold relies on using R_0 as our key parameter, which is the relative speed of spread if the initial population is disease-free. During the Covid-19 pandemic, proposals were made for a “herd immunity strategy” that would essentially rely on the disease dying away once enough of the population – the herd immunity threshold – had immunity, but immunity that was acquired after being infected, not through vaccination. There are two big drawbacks to this approach (see chapter references).

Firstly, as we have seen here, if the disease persists for any length of time – such that birth and death processes start to introduce population turnover – we’d be unwise to assume the disease would naturally burn out. There are a whole host of other factors that would also allow for a disease to persist in the long-term, including waning immunity and lack of immunity to mutations to name just two.

Secondly, and related to R_0, if we use the herd immunity threshold as calculated above, the idea is that we have protected (almost) the whole population until a vaccine was developed, and then we can move enough individuals out of the susceptible compartment to stop the disease from ever spreading. If instead we allow the disease to spread the key parameter is no longer R_0 but R_t, as discussed in the last chapter. We know R_t equals 1 exactly at the peak of an epidemic, and at this point the proportion that has been infected will be the same as the herd immunity threshold. However, in this approach we now have a potentially large proportion of the population who are still infected, as well as a large proportion who remain susceptible. So while the number of infections will start to fall (because we are in the region of R_t\lt1) we will still have a considerable number of onward infections. Say that we calculate that a disease has R_0=3, meaning the herd immunity threshold is 66.7%. As we have just said, that is now the proportion that has been infected at the peak of the infection. By the time the infection has actually died out, the proportion infected will be over 90%.

The lesson here is to be careful about how you interpret models!

Key Takeaways

  • We can extend the SIR model to be more realistic for long-term predictions by including births and deaths.
  • We now get an endemic equilibrium, where the disease will remain in the population over the long term.
  • R_0 remains a key quantity, and is helpful in determining the threshold for herd immunity.

Chapter references

  • The discussion on herd immunity is based on the paper by Best & Ashby (2022).
  • The content in the Infectious diseases section was influenced by the textbooks, Mathematical Biology 1 by Murray and Modelling Infectious Diseases in Humans and Animals by Keeling & Rohani.

Annotate

Next Chapter
Diseases of ecological populations
PreviousNext
Biology

Copyright © 2023

            Introducing Mathematical Biology by Alex Best is licensed under a Creative Commons Attribution 4.0 International License, except where otherwise noted.
Powered by Manifold Scholarship. Learn more at
Opens in new tab or windowmanifoldapp.org