Skip to main content

Mathematical Epidemiology and More: A Course on Compartmental Models

Meredith L. Greer

Chapter 14 * The Next Generation Method

***This chapter is under construction.***
This chapter uses many linear algebra concepts. Before beginning, we pause to review.

Exploration 14.1. Linear Algebra Review.

To do the work of this chapter, it will help to have certain linear algebra concepts freshly in mind. This Exploration prepares us for what comes next.

(a)

First review matrices. Use textbooks, online information, or other means to refresh your about matrices in general. Then respond to the following.
  1. Create a \(2 \times 2\) matrix. Create another matrix, this time \(3 \times 3\text{.}\) These matrices are called square matrices. Why can we use the name “square”? And what could a matrix look like if it is not a square matrix? To demonstrate, create at least two differently sized matrices that are not square, and explain how you know that these matrices are not square.
  2. Create two different \(2 \times 2\) matrices and demonstrate how to multiply them.
  3. Describe what an inverse matrix is. Then write the formula for computing the inverse of a matrix of the form:
    \begin{equation*} \begin{pmatrix} a \amp b \\ c \amp d \end{pmatrix}. \end{equation*}
    For each of the two \(2 \times 2\) matrices you created in part (2) above, show how to compute the inverse of the matrix, using the matrix inverse formula you wrote. If it happens that either of your matrices has no inverse, explain how you know this. If neither of your matrices has an inverse, then create a new \(2 \times 2\) matrix that does have an inverse, and compute its inverse.

(b)

Next review eigenvalues.
  1. How do we describe an eigenvalue mathematically? Your response should include at least one mathematical formula, along with a description of what the formula means.
  2. Research eigenvalues to find ways that eigenvalues are used in applications of mathematics. Write descriptions of at least two examples of applications of eigenvalues.

Section 14.1 Mathematical Background for the Next Generation Method

So far, to determine the value of \(\beta\) in a compartmental model, we have done one of two things:
  1. estimated \(\beta\) visually by comparing the model’s solution curves with data, or
  2. used a formula for \(\mathcal{R}_0\text{,}\) such as \(\mathcal{R}_0 = (\beta N)/\gamma\text{,}\) to substitute in values we know for \(\mathcal{R}_0\text{,}\) \(N\text{,}\) and \(\gamma\text{,}\) rearrange, and solve for \(\beta\text{.}\)
However, the formula for \(\mathcal{R}_0\) is not always \(\mathcal{R}_0 = (\beta N)/\gamma\text{.}\) Systems of differential equations having different compartments, or systems in which the flow between compartments is different than what we have worked with before, can have different formulas for \(\mathcal{R}_0\text{.}\) These formulas are usually not obvious. However, a method exists that allows us to determine the formula. This is the Next Generation Method.
To use the Next Generation Method, we will need to be familiar with matrices, eigenvalues, equilibrium values for a system of differential equations, and partial derivatives. Exploration 14.1 was a reminder about matrices and eigenvalues, and we will build more knowledge about these linear algebra concepts. But first, use Activity 14.2 and Activity 14.3 to learn about equilibrium values and partial derivatives.

Activity 14.2.

This activity identifies equilibrium values for two systems of differential equations that we have worked with as compartmental models.
  1. Consider the SI model below:
    \begin{align*} \frac{dS}{dt} \amp = - \beta S I \\ \frac{dI}{dt} \amp = \beta S I. \end{align*}
    For a system of differential equations, being at equilibrium means that the dependent variables in the system are not changing as the independent variables change. For the SI model, this means that the populations \(S\) and \(I\) are not changing as time \(t\) changes. Using knowledge from calculus: what can we say about each of the derivatives \(dS/dt\) and \(dI/dt\) when the system is at equilibrium?
  2. Next: determine the possible values of \(S\) and \(I\) for which the SI model is at equilibrium.
  3. Next consider the SEIR model below:
    \begin{align*} \frac{dS}{dt} \amp = - \beta S I \\ \frac{dE}{dt} \amp = \beta S I - \kappa E \\ \frac{dI}{dt} \amp = \kappa E - \gamma I \\ \frac{dR}{dt} \amp = \gamma I. \end{align*}
    Can you determine possible values of \(S\text{,}\) \(E\text{,}\) \(I\text{,}\) and \(R\) for which the SEIR model is at equilibrium? In doing this, be sure the overall population is nonzero, that is, \(S+E+I+R \neq 0\text{.}\)
Answer.
  1. For the system to be at equilibrium, each of the derivatives must equal \(0\text{.}\)
  2. If \(S=0\text{,}\) the SI model is at equilibrium. This corresponds to a situation where the entire population is in the \(I\) compartment.
    Alternatively, if \(I=0\text{,}\) the SI model is at equilibrium. This corresponds to a situation where the entire population is in the \(S\) compartment.
  3. The SEIR model is at equilibrium when \(E=0\) and \(I=0\text{.}\) The rest of the population can be divided in any way between the \(S\) and \(R\) compartments.
You may have noticed in Activity 14.2, part 3, that there are many ways to arrive at equilbrium. It is true that the \(E\) and \(I\) compartments should be empty, but the population can be distributed in any way across compartments \(S\) and \(R\text{.}\)
For the work we will do in this chapter, we decide the population distribution by seeking the disease-free equilibrium. This is the particular equilibrium value that occurs in a situation where disease has not yet appeared in a population. In the SEIR model, for instance, equilibrium requires \(E=I=0\text{,}\) with the rest of the population split somehow across the \(S\) and \(R\) compartments. Disease-free equilibrium adds the extra requirement that the population all be in the \(S\) compartment. This means \(S=N\text{,}\) where \(N\) represents the entire population, as first defined in Activity 5.2. Then the other three populations, \(E\text{,}\) \(I\text{,}\) and \(R\text{,}\) all equal \(0\text{.}\)
One way to write the disease-free equilibrium is \((\overline{S}, \overline{E}, \overline{I}, \overline{R})=(N, 0, 0, 0)\text{.}\) The bar over each of the population letters \(S\text{,}\) \(E\text{,}\) \(I\text{,}\) and \(R\) represents an equilibrium value for that population.
For the models we have studied so far, disease-free equilibrium requires that the entire population start in the \(S\) compartment, but there exist compartmental models where the initial population can be split across multiple compartments.
In addition to the concept of the disease-free equilibrium value, we also work with partial derivatives when applying the Next Generation Method. Gain practice with partial derivatives in Activity 14.3.

Activity 14.3.

This course requires mainly the calculus concepts from a Calculus I course. These concepts include working with derivatives, understanding that a derivative is a rate of change, knowing that the slope on a graph relates directly to the derivative of a function, and recognizing that positive slopes (or, equivalently, positive derivatives) correspond to increasing functions and, similarly, negative slopes (or, equivalently, negative derivatives) correspond to decreasing functions.
The Next Generation Method also requires us to compute partial derivatives. To see how this works, progress through the steps below.
  1. Start by computing derivatives as you would during Calculus I.
    1. For the function \(f(x)=3x^2 + 5x + 7\text{,}\) compute its derivative \(f'(x)\text{.}\) This can also be written \(\frac{d}{dx}(3x^2 + 5x + 7)\text{.}\)
    2. Then, for the related function \(g(x)=ax^2 + bx + c\text{,}\) compute its derivative. The values \(a\text{,}\) \(b\text{,}\) and \(c\) are all constants. This derivative can also be written \(\frac{d}{dx}(ax^2 + bx + c)\text{.}\)
  2. Next, we introduce partial derivatives. Whereas in part (1) we computed derivatives with respect to \(x\text{,}\) we now allow ourselves to compute derivatives with respect to other variables in a formula.
    1. Consider this differential equation from the SEIR model: \(\frac{dE}{dt} = \beta S I - \kappa E\text{.}\) Compute the derivative of the right-hand side, with respect to \(E\text{.}\) That is, compute \(\frac{\partial}{\partial E}(\beta S I - \kappa E)\text{.}\) When computing the derivative with respect to \(E\text{:}\) think of \(S\text{,}\) \(I\text{,}\) \(\beta\text{,}\) and \(\kappa\) as constants. Think of \(E\) as the only variable.
      You may notice that we have switched from the notation \(\frac{d}{dx}\) to the notation \(\frac{\partial}{\partial E}\text{.}\) We use \(\partial\) rather than \(d\) when there are multiple choices for our variable of interest: we are focusing on \(E\text{,}\) but we could have chosen any other variable in our system of equations. The idea of partial derivatives is that there are multiple variables involved, and as each of them changes, the entire system is affected, but we focus on one variable at a time, therefore looking at part of the change of the system rather than the total change of the system.
      This emphasis on one variable at a time benefits our understanding by showing how change in each variable—in our case, change in each compartment—affects the outcome of the entire system.
    2. Then compute these two derivatives: \(\frac{\partial}{\partial S}(\beta S I - \kappa E)\) and \(\frac{\partial}{\partial I}(\beta S I - \kappa E)\text{.}\) The first of these thinks of only \(S\) as a variable. The second thinks of only \(I\) as a variable.
Answer.
    1. The derivative is \(\frac{d}{dx}(3x^2 + 5x + 7) = 6x+5\text{.}\)
    2. The derivative is \(\frac{d}{dx}(ax^2 + bx + 7c) = 2ax+b\text{.}\)
    1. The derivative is \(\frac{\partial}{\partial E}(\beta S I - \kappa E) = - \kappa \text{.}\)
    2. The derivative with respect to \(S\) is \(\frac{\partial}{\partial S}(\beta S I - \kappa E) = \beta I \) and the derivative with respect to \(I\) is \(\frac{\partial}{\partial I}(\beta S I - \kappa E) = \beta S \text{.}\)

Section 14.2 The Next Generation Method: Steps to Follow

We now go through the steps of the Next Generation Method.

Activity 14.4.

Begin with the SEIR model, shown below, and follow the steps described. While in many activities it is good to try all the steps before checking the answer, in this activity you may wish to compare your results after each step with the “Answer” link at the bottom of the activity. There will then be additional chances to practice the Next Generation Method, so that you will have the chance to try the Next Generation Method without checking answers at each step.
\begin{align*} \frac{dS}{dt} \amp = - \beta S I \\ \frac{dE}{dt} \amp = \beta S I - \kappa E \\ \frac{dI}{dt} \amp = \kappa E - \gamma I \\ \frac{dR}{dt} \amp = \gamma I. \end{align*}
  1. Decide which compartments of the SEIR model involve individuals experiencing some stage of infection. If a person has contracted an illness, they are in some stage of infection throughout the time they are incubating, whether contagious or not, and whether showing symptoms or not. They stop being in a stage of infection when they are completely done with the illness.
  2. For the equation(s) selected in Step 1, write the equation(s) in matrix form. To do so, place the derivatives on the left-hand side into a vector. On the right-hand side, write the equations in a vector.
  3. Separate the right-hand side into two vectors, \(\widehat{F}\) and \(\widehat{V}\text{,}\) so that all of the following are true:
    • the right-hand side can be written as \(\widehat{F}-\widehat{V}\text{;}\)
    • vector \(\widehat{F}\) includes all new infection terms, that is, the terms within the model that indicate when individuals who were not at any stage of infection newly enter a compartment in which they are at some stage of infection; and
    • vector \(\widehat{V}\) is the vector of all other movement between compartments.
  4. Compute the disease-free equilibrium for your compartmental model.
  5. Compute matrices \(F\) and \(V\) (notice these do not have “hats”: they are different from \(\widehat{F}\) and \(\widehat{V}\)) by computing partial derivatives of each of the terms in \(\widehat{F}\) and \(\widehat{V}\text{,}\) with respect to each of the populations of compartments experiencing infection, following the pattern shown below in the case of the SEIR model, and then evaluating your results at the disease-free equilibrium:
    \begin{align*} F = \left. \begin{pmatrix} \frac{\partial (\beta SI)}{\partial E} \amp \frac{\partial (\beta SI)}{\partial I} \\ \frac{\partial (0)}{\partial E} \amp \frac{\partial (0)}{\partial I} \end{pmatrix} \right| _{(S, E, I, R) = (N, 0, 0, 0)} \end{align*}
    \begin{align*} V = \left. \begin{pmatrix} \frac{\partial (\kappa E)}{\partial E} \amp \frac{\partial (\kappa E)}{\partial I} \\ \frac{\partial (-\kappa E + \gamma I)}{\partial E} \amp \frac{\partial (-\kappa E + \gamma I)}{\partial I} \end{pmatrix} \right| _{(S, E, I, R) = (N, 0, 0, 0)} . \end{align*}
  6. Compute the Next Generation Matrix, which is the matrix
    \begin{gather*} K = F \, V^{-1} \end{gather*}
    and determine the largest positive eigenvalue of \(K\text{.}\)
    To do this: first compute \(V^{-1}\text{.}\) Then compute \(F \,V^{-1}\text{.}\) Then determine the eigenvalues of \(F \,V^{-1}\text{.}\)
    The largest possible eigenvalue of \(F \,V^{-1}\) is the value of \(\mathcal{R}_0\text{.}\)
Answer.
  1. For the SEIR model, the E and I compartments are stages of infection. The S and R compartments are not stages of infection.
  2. \begin{gather*} \begin{pmatrix} \frac{dE}{dt} \\ \frac{dI}{dt} \end{pmatrix} = \begin{pmatrix} \beta S I - \kappa E \\ \kappa E - \gamma I \end{pmatrix} \end{gather*}
  3. \begin{gather*} \begin{pmatrix} \frac{dE}{dt} \\ \frac{dI}{dt} \end{pmatrix} = \begin{pmatrix} \beta S I \\ 0 \end{pmatrix} - \begin{pmatrix} \kappa E \\ -\kappa E + \gamma I \end{pmatrix} \end{gather*}
    That is:
    \begin{gather*} \widehat{F} = \begin{pmatrix} \beta S I \\ 0 \end{pmatrix} \end{gather*}
    and
    \begin{gather*} \widehat{V} = \begin{pmatrix} \kappa E \\ -\kappa E + \gamma I \end{pmatrix} \end{gather*}
  4. We did this for the SEIR model in Activity 14.2. The disease-free equilibrium is \((S, E, I, R)=(N, 0, 0, 0)\text{,}\) where \(N\) is the total number of people in the entire population being studied.
  5. \begin{align*} F = \left. \begin{pmatrix} 0 \amp \beta S \\ 0 \amp 0 \end{pmatrix} \right| _{(S, E, I, R) = (N, 0, 0, 0)} = \begin{pmatrix} 0 \amp \beta N \\ 0 \amp 0 \end{pmatrix} \end{align*}
    \begin{align*} V = \left. \begin{pmatrix} \kappa \amp 0 \\ -\kappa \amp \gamma \end{pmatrix} \right| _{(S, E, I, R) = (N, 0, 0, 0)} = \begin{pmatrix} \kappa \amp 0 \\ -\kappa \amp \gamma \end{pmatrix} \end{align*}
  6. Compute matrix \(V^{-1}=\begin{pmatrix} \frac{1}{\kappa} \amp 0 \\ \frac{1}{\gamma} \amp \frac{1}{\gamma} \end{pmatrix}.\) Then compute matrix \(FV^{-1}=\begin{pmatrix} \frac{\beta N}{\gamma} \amp \frac{\beta N}{\gamma} \\ 0 \amp 0 \end{pmatrix}.\) The eigenvalues of \(FV^{-1}\) are \(0\) and \(\frac{\beta N}{\gamma}\text{.}\) Of these, the largest eigenvalue is \(\frac{\beta N}{\gamma}\text{,}\) and therefore \(\mathcal{R}_0 = \frac{\beta N}{\gamma}\text{.}\)

Activity 14.5.

Now that we have shown all the steps for computing a formula for \(\mathcal{R}_0\) in the SEIR model, try again for the SIRS model, shown in Figure 14.1.
(Write this)
Figure 14.1. The compartmental diagram for one possible SIRS model, where “S”stands for “Susceptible”, “I” stands for “Infectious”, and “R” stands for “Recovered”. In this model, individuals lose immunity some time after recovering. They then return to the S compartment.
Hint.
Only one compartment of the SIR model involves individuals experiencing some stage of infection, namely, the I compartment. This means that when we write the equations in matrix form, we use a \(1 \times 1\) matrix. Then \(\widehat{F}= \begin{pmatrix} \beta SI \end{pmatrix}\) and \(\widehat{V}= \begin{pmatrix} \gamma I \end{pmatrix}.\) Computing partial derivatives with respect to \(I\text{,}\) and evaluating at the disease-free equilibrium \((S, I, R)=(N, 0, 0)\text{,}\) results in \(F = \begin{pmatrix} \beta N \end{pmatrix}\) and \(V = \begin{pmatrix} \gamma \end{pmatrix}.\) Computing \(FV^{-1}\) results in \(FV^{-1} = \begin{pmatrix} \frac{\beta N}{\gamma} \end{pmatrix}\text{,}\) and the one eigenvalue of \(FV^{-1}\) tells us that \(\mathcal{R}_0 = \frac{\beta N}{\gamma}\text{.}\)
We have now seen the steps of the Next Generation Method for both SEIR and SIR models. The final step of the Next Generation Method requires a significant amount of linear algebra, especially as models become more complicated. In Section 14.3, we introduce Python code to help with this.

Section 14.3 Python Code for the Next Generation Method

For the most part, steps 1 through 4 of the Next Generation Method cannot be automated, as they require human judgment. In particular, Step 1 requires determining which compartments involve some stage of infection, and in Step 3, we must distinguish new infections from all other stages of infection. In Step 4, we use logical thinking to help us write down the population distribution in a situation where there has not been an outbreak previously.
For the later steps of the Next Generation Method, especially computing matrix inverses and eigenvalues, technology can be extremely helpful. These are doable using linear algebra, yet they are not the focus of this course. For this textbook, students are asked to complete by hand Step 5 of the Next Generation Method, and then we use Python for the grand finale in Step 6.
Activity 14.6 provides one possible set of Python commands for computing the value of \(\mathcal{R}_0\) in the SEIR model shown in Activity 14.4.

Activity 14.6.

The code block below produces the eigenvalues for the SEIR model, doing the computation in Step 6 of the Next Generation Method as described in Activity 14.4. This code can also be used as a template for producing eigenvalues when the Next Generation Method produces a \(2\times 2\) matrix.
Use the code to answer the questions below.
  1. Four expressions are listed as Symbol in lines 4 through 7 of the Python code. Compare these expressions with the entries in matrices \(F\) and \(V\) that we solved for in part 5 of Activity 14.4. Show where each Symbol occurs in matrices \(F\) and \(V\text{,}\) and state why these are all the terms we need to define as Symbol so that Python will recognize them.
  2. Use lines 10 and 14 of the Python code to describe the syntax we use for entering matrices into Python.
  3. From line 18, state how we compute the inverse of matrix \(V\text{,}\) also written as \(V^{-1}\text{.}\) From line 20, state how we compute the matrix \(FV^{-1}\text{,}\) which multiples matrix \(F\) with matrix \(V^{-1}\text{.}\) From line 25, state the command used to calculate the eigenvalues of matrix \(FV^{-1}\text{.}\)
  4. The lines 11, 15, 19, 21, and 26 all include the print command. Write out exactly what Python prints for each of these lines when we evaluate the code, and explain what each of these printed results means, in the context of the Next Generation Method.
  5. When we compute eigenvalues for a matrix, it is possible that an eigenvalue may occur more than once. The number of times an eigenvalue occurs is called the multiplicity of the eigenvalue. Python tells us the eigenvalues of a matrix, along with the multiplicity of each eigenvalue. Write out the eigenvalues of \(FV^{-1}\) and the multiplicity of each eigenvalue. Explain how you know that each eigenvalue is a real number (and not a complex number, meaning that no part of the eigenvalue includes \(i=\sqrt{-1}\)). Then state which eigenvalue is the largest, and confirm that this eigenvalue must represent a positive real number.
Hint.
Only
Activity 14.7 below introduces a Python code block for producing eigenvalues when the Next Generation Method involves a \(3\times 3\) matrix.

Activity 14.7.

The book Mathematical Epidemiology
 1 
Mathematical Epidemiology was edited by Fred Brauer, Pauline van den Driessche, and Jianhong Wu, and published by Springer in 2008. The model we describe comes from Chapter 12, which was written by Fred Brauer, and is labeled SLIAR. Our version is labeled SEIAR to stay consistent with notation we have used in other chapters.
introduces a model of influenza having five compartments. This model has a Susceptibles (S) compartment that is similar to what we have seen before. There is then an Exposed-and-incubating (E) compartment for people incubating influenza who are not yet able to spread influenza to others. After that, people may be Infectious (I), meaning they have symptoms and can spread influenza, or they may be Asymptomatic (A), meaning they do not have symptoms but can still spread influenza to others. In this model, Asymptomatic (A) people are less able to spread influenza than Infectious (I) people. Many people then move into a Recovered (R) compartment, but those who die depart the model rather than moving into the R compartment. The compartmental diagram for this model appears in Figure 14.2. In this activity, we explore the SEIAR model and investigate Python code that can help us determine \(\mathcal{R}_0\) for this model using the Next Generation Method.
(Write this)
Figure 14.2. The compartmental diagram for one possible SEIAR model for influenza, where “S”stands for “Susceptible”, “E” stands for “Exposed-and-incubating”, “I” stands for “Infectious”, “A” stands for “Asymptomatic” and “R” stands for “Recovered”. In this model, Asymptomatic individuals have no symptoms but can spread influenza, though they are not as contagious as individuals in the Infectious compartment.
  1. The parameter \(p\) may take on values \(0 \leq p \leq 1\text{.}\) Use the questions below to think through the role of \(p\) in the SEIAR model.
    1. What does it mean that the arrow from E to I includes the term \(p \kappa E\text{,}\) and the arrow from E to A includes the term \((1-p) \kappa E\text{?}\) What would happen if we added the terms on these two arrows and simplified the result?
    2. How would the result in part (a) compare with the formula for moving from the E compartment to the I compartment in a simpler model with an Exposed-and-incubating compartment, such as the model in Activity 7.4?
  2. Asymptomatic people in the SEIAR model are able to spread influenza but may not be able to spread influenza as easily as people in the Infectious compartment. Given this information, what possible values for \(\delta\) make sense?
  3. This model separates Recovered people from those who die from influenza. We have seen that the R compartment stands for Recovered. As further information, the possible values for the parameter \(f\) are \(0 \leq f \leq 1\text{.}\) Explain how \(f\) distinguishes people who have died from people who have recovered, and compare this approach with how we have handled deaths in models where the R compartment stood for Removed.
Answer.
    1. The parameter \(p\) splits up the term \(\kappa E\text{,}\) sending a fraction \(p\) from the Exposed-and-incubating compartment to the Infectious compartment, and the other fraction \(1-p\) to the Asymptomatic compartment.
      Notice that since \(0 \leq p \leq 1\text{,}\) then \(0 \leq 1-p \leq 1\text{.}\) This means that \(0 \leq p\kappa E \leq \kappa E\) and \(0 \leq (1-p)\kappa E \leq \kappa E\text{.}\)
      Adding the two terms: \(p\kappa E + (1-p)\kappa E = p \kappa E + \kappa E - p \kappa E = \kappa E\text{.}\)
    2. In a simpler model, the expression \(\kappa E\) would represent the full number of people transitioning from the Exposed-and-incubating compartment to the Infectious compartment within a single unit of time. In the SEIAR model, \(\kappa E\) is split into two pieces, \(p\kappa E\) and \((1-p)\kappa E\text{,}\) so that the total number of people departing the E compartment is the same, but those people go to two possible compartments, the Infectious compartment and the Asymptomatic compartment.
  2. It makes sense to have \(0 \lt \delta \lt 1\text{.}\) Notice that letting \(\delta = 0\) would means there is no spread of influenza from Asymptomatic people, and letting \(\delta = 1\) would means Asymptomatic people spread influenza just as much as people in the Infectious compartment. For values of \(\delta\) strictly between \(0\) and \(1\text{,}\) Asymptomatic people spread influenza to some extent, but strictly less than people in the Infectious compartment
    For this model, we wrote that Asymptomatic people spread influenza less than Infectious people, so we should not allow \(\delta \gt 1\text{.}\) However, the idea of \(\delta \gt 1\) could be useful in some other model with different assumptions.
  3. In a single unit of time, \(\alpha I\) people depart compartment I. A fraction \(f\) of these people move into compartment R, meaning they have recovered. A fraction \(1-f\) of these people do not move into a named compartment. These people die, and in this model, people who have died are no longer represented as being in a compartment. This means that in this model, the total population of people in compartments does not stay the same across time.
The code block below produces the eigenvalues for the SLIAR model. This code can also be used as a template for producing eigenvalues when the Next Generation Method produces a \(3\times 3\) matrix. Comparing this code with the previous code may help show how similar the code is when computing the Next Generation matrix for different models, as well as showing what to change for a different model.
PrevTopNext