Chapter 14 The Next Generation Method

Goals

Learn to complete all the steps in the Next Generation Method for relating the epidemiological quantity \(\mathcal{R}_0\) to parameters and other quantities used in a compartmental model
Build understanding of compartmental models as a way to map disease spread, or other kinds of spread, as a generation-to-generation phenomenon
Increase awareness of how eigenvalues show us specific information about the long-term behavior of a system of differential equations

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.

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.
Create two different \(2 \times 2\) matrices and demonstrate how to multiply them.
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.

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

estimated \(\beta\) visually by comparing the model’s solution curves with data, or
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.

Consider the SI model below:

\begin{align*} \frac{dS}{dt} \amp = - \beta S I \\ \frac{dI}{dt} \amp = \beta S I \text{.} \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?
Next: determine the possible values of \(S\) and \(I\) for which the SI model is at equilibrium.
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 \text{.} \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.

For the system to be at equilibrium, each of the derivatives must equal \(0\text{.}\)
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.
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.3. 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.

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{.}\)
🔗
🔗
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. The concept of a next-generation matrix comes from Diekmann and Heesterbeek.

A good source to read on this is the book Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation by O. Diekmann and J.A.P. Heesterbeek, published in the Wiley Series in Mathematical and Computational Biology in 2000.

The steps below appear in the paper “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission” by van den Driessche and Watmough.

The authors are P. van den Driessche and James Watmough, and the article was published by the journal Mathematical Biosciences in 2002.

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 given at the bottom of the step. There will then be additional chances to practice the Next Generation Method, so that you will have the chance to try this process 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 \text{.} \end{align*}

(a)

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.

Answer.

For the SEIR model, the E and I compartments are stages of infection. The S and R compartments are not stages of infection.

(b)

For the equation(s) selected in Step (a), 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.

Answer.

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

(c)

Separate the right-hand side obtained in Step (b) 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.

Answer.

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

(d)

Compute the disease-free equilibrium, as described within and just after Activity 14.2, for your compartmental model.

Answer.

The disease-free equilibrium is \((\overline{S}, \overline{E}, \overline{I}, \overline{R})=(N, 0, 0, 0)\text{,}\) where \(N=S+E+I+R\) is the entire population being studied.

(e)

Compute matrices \(F\) and \(V\text{,}\) noticing these are different from \(\widehat{F}\) and \(\widehat{V}\text{,}\) 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*}

Answer.

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

(f)

Compute the Next Generation Matrix, which is the matrix

\begin{equation*} K = F \, V^{-1} \end{equation*}

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 positive eigenvalue of \(F \,V^{-1}\) is the value of \(\mathcal{R}_0\text{.}\)

Answer.

Compute matrix

\begin{equation*} V^{-1}=\begin{pmatrix} \frac{1}{\kappa} \amp 0 \\ \frac{1}{\gamma} \amp \frac{1}{\gamma} \end{pmatrix}\text{.} \end{equation*}

Then compute matrix

\begin{equation*} FV^{-1}=\begin{pmatrix} \frac{\beta N}{\gamma} \amp \frac{\beta N}{\gamma} \\ 0 \amp 0 \end{pmatrix}\text{.} \end{equation*}

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

\begin{equation*} \mathcal{R}_0 = \frac{\beta N}{\gamma}\text{.} \end{equation*}

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

described in detail following the image — Figure 14.6. 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

\begin{equation*} \widehat{F}= \begin{pmatrix} \beta SI \end{pmatrix} \end{equation*}

and

\begin{equation*} \widehat{V}= \begin{pmatrix} \gamma I \end{pmatrix}\text{.} \end{equation*}

Computing partial derivatives with respect to \(I\text{,}\) and evaluating at the disease-free equilibrium \((S, I, R)=(N, 0, 0)\text{,}\) results in

\begin{equation*} F = \begin{pmatrix} \beta N \end{pmatrix} \end{equation*}

and

\begin{equation*} V = \begin{pmatrix} \gamma \end{pmatrix}\text{.} \end{equation*}

Computing \(FV^{-1}\) results in

\begin{equation*} FV^{-1} = \begin{pmatrix} \frac{\beta N}{\gamma} \end{pmatrix} \end{equation*}

and the one eigenvalue of \(FV^{-1}\) tells us that

\begin{equation*} \mathcal{R}_0 = \frac{\beta N}{\gamma}\text{.} \end{equation*}

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 (a) through (d) of the Next Generation Method cannot be automated, as they require human judgment. In particular, Step (a) requires determining which compartments involve some stage of infection, and in Step (c), we must distinguish new infections from all other stages of infection. In Step (d), 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 (e) of the Next Generation Method, and then we use Python for the grand finale in Step (f).

Activity 14.7 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.7.

The code block below produces the eigenvalues for the SEIR model. In other words, this code does the computation in Step (f) 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.

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 Step 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.
Use lines 10 and 14 of the Python code to describe the syntax we use for entering matrices into Python.
From line 18, state how we tell Python to compute the inverse of matrix \(V\text{,}\) also written as \(V^{-1}\text{.}\) From line 20, state how we tell Python to 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{,}\) and give the name we assign to that set of eigenvalues.
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.
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.

In Step 5 of Activity 14.4, we produced matrices for \(F\) and \(V\text{.}\) These matrices used letters, including Greek letters, to represent some parts of the model. We now designate each of those letters to be a Symbol so that Python knows how to work with it. As an example, we write gamma = Symbol('gamma') to give the name “gamma” to the Greek letter \(\gamma\text{.}\) Python then knows to treat the entire text string gamma as the single parameter value \(\gamma\) whenever \(\gamma\) appears in a mathematical expression, such as in matrix \(V\text{.}\)
To enter a matrix, we type Matrix([]). Inside the square brackets, we enter each row of the matrix, using its own set of square brackets, with commas separating the matrix entries. This system extends to matrices larger than \(2\times 2\text{.}\)
These commands appear on the right-hand side of the equals \((=)\) sign in each of the given lines. Notice that in each case we name the output. In particular, we give the name R0eval to the output of the eigenvals command so that we can print this output, using the command in line 26.
The print commands in this code block all include text, in quotation marks, followed by the name of a mathematical expression we have either entered into Python or asked Python to compute. As a result, Python prints this exact text first, telling us what we are seeing in the mathematical expression appearing just after the text.

For instance, the code in line 19 tells Python to print the text “V inverse =” and then to print the actual matrix \(V^{-1}\) that has been computed based on the matrix \(V\) we entered.
For matrix \(FV^{-1}\) resulting from the Next Generation Method shown in Activity 14.4, the eigenvalues are 0, having multiplicity \(1\text{,}\) and N*beta/gamma, which equals \(\frac{\beta N}{\gamma}\text{,}\) having multiplicity 1.

It is not unusual for the eigenvalue equaling \(\mathcal{R}_0\) to look different in Python than in mathematical writing. Rearranging and rewriting the Python version of \(\mathcal{R}_0\) often makes the result easier to work with, both mathematically and for understanding its meaning.

In the case of this activity’s output: the eigenvalue N*beta/gamma is a positive real number, and is larger than 0, because all three values \(N\text{,}\) \(\beta\text{,}\) and \(\gamma\) are positive real numbers in the SEIR model that led to this result.

Activity 14.8 below provides a somewhat larger model on which to practice the Next Generation Method. It also introduces a Python code block for producing eigenvalues when the Next Generation Method involves a \(3\times 3\) matrix.

Activity 14.8.

The book Mathematical Epidemiology

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

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.8?
🔗
🔗
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?
This model separates Recovered people from those who die from influenza. To emphasize, the R compartment here stands for Recovered, which is different than what we have seen in some models where the R stood for Removed, which could mean recovery, death, or other departure from the population. 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.
The code block below produces the eigenvalues for the SEIAR model. This code can also be used as a template for producing eigenvalues when the Next Generation Method produces a \(3\times 3\) matrix. Compare this code with the code in Activity 14.7. Name at least two ways the code blocks are similar, indicating ways the Next Generation Method computation is similar, and name at least two ways the code blocks are different, indicating things that change for the Next Generation Method computation in these two different examples.
One way to write the value of \(\mathcal{R}_0\) for this model is the following:

\begin{align*} \mathcal{R}_0 \amp = \beta N \left( \frac{p}{\alpha} + \frac{\delta (1-p)}{\eta} \right) . \end{align*}

Write the value of \(\mathcal{R}_0\) given by the Python code block below, and show the algebra steps needed to demonstrate that Python’s result matches the formula shown here. Then explain, mathematically, what the fractions \(\frac{p}{\alpha} \) and \(\frac{\delta (1-p)}{\eta}\) indicate about movement through the SEIAR model.

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.
🔗
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.
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 (unless we let \(f=1\text{,}\) sending everyone from the \(I\) compartment into the \(R\) compartment).
The answers here are open-ended.
Python writes the eigenvalue that is \(\mathcal{R}_0\) in the form

\begin{gather*} -N \beta (\alpha \delta p - \alpha \delta - \eta p)/(\alpha \eta ) . \end{gather*}

First, we write this in more standard fraction notation. In the same step, we move the minus sign (\(-\)) from before \(N \beta\text{,}\) choosing instead to distribute this minus sign to the terms inside the parentheses:

\begin{gather*} \frac{N \beta (-\alpha \delta p + \alpha \delta + \eta p)}{\alpha \eta}. \end{gather*}

Next, we group terms within the parentheses, splitting the terms up so they form two separate fractions, and rearrange the order of some terms to better match the answer we eventually want to reach:

\begin{gather*} \beta N \left( \frac{\eta p}{\alpha \eta} + \frac{\alpha \delta -\alpha \delta p}{\alpha \eta} \right). \end{gather*}

Finally, we divide \(\eta\) from the top and bottom of the first fraction, and we divide \(\alpha\) from the top and bottom of the second fraction. In doing so, we note that these are parameters that should always have strictly positive values, which means we can divide by them and not be dividing by \(0\text{:}\)

\begin{gather*} \beta N \left( \frac{p}{\alpha} + \frac{\delta (1-p)}{\eta} \right). \end{gather*}

For mathematical explanation, we compare this form of \(\mathcal{R}_0\) with \(\mathcal{R}_0=\frac{\beta N}{\gamma}\) in some of the simpler models we have studied. In those, the value of \(\gamma\) in the denominator indicates that the more rapidly people depart the I compartment, the less able a disease is to spread, making \(\mathcal{R}_0\) lower. Here, since we split into compartments I and A, we need to separately evaluate the rapidity with which people depart each of these compartments. A fraction \(p\) of people in the I compartment depart at rate \(\alpha\text{.}\) This is the first fraction in parentheses for \(\mathcal{R}_0\) in this exercise. The remaining fraction of people, \(1-p\text{,}\) depart compartment A at rate \(\eta\text{.}\) However, these people also have less impact on new infections because of the parameter \(\delta\text{.}\) This lesser impact is measured by the value of \(\delta\) in the second fraction within \(\mathcal{R}_0\text{.}\)

When using the Next Generation Method to work with real-world values of \(\mathcal{R}_0\text{,}\) be sure to focus your attention on models representing disease spread without human interventions such as vaccines and isolation. The quantity \(\mathcal{R}_0\) does not apply in these cases because the definition of \(\mathcal{R}_0\) involves a single infectious case in an otherwise completely Susceptible population. By design, public health measures such as vaccines and isolation cause a population to not be fully Susceptible.

In book chapters where we use \(\mathcal{R}_0\) to compute the value of \(\beta\) in a model, we first work with a model having no intervention. In Chapter 6, in Exercise 1 of the For Further Thought 6.4 section, we start with a standard SIR model to compute \(\mathcal{R}_0\text{.}\) The exercise says to use (5.1), which is the formula that results from the Next Generation Method. Once \(\beta\) is computed, we can then experiment with different vaccine levels to change the model. This topic returns in Chapter 8, with a multi-step explanation near the end of Section 8.2. The explanation walks through steps for first estimating \(\beta\) using a model without isolation, using the appropriate formula generated by the Next Generation Method, and only then building isolation into the model.

In summary, knowing how to use the Next Generation Method provides the freedom to develop new models and to then correctly determine the formula relating a model’s parameters with the epidemiological constant \(\mathcal{R}_0\text{.}\) When using the Next Generation Method, take care to keep in mind the complete definition of \(\mathcal{R}_0\text{,}\) both for biological accuracy and so that the model’s parameters will be reasonable when running simulations.

For Further Thought 14.4 For Further Thought

1.

Consider the SEIR model with births and deaths, shown in Figure 14.10, which includes a constant term \(\Pi\) representing births per time unit, as well as a death rate \(\mu\) from each compartment. This model comes from the book Mathematical Epidemiology.

⁴

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 6, which was written by James Watmough and Pauline van den Driessche.

Complete all steps of the Next Generation Method to show how to determine the value of \(\mathcal{R}_0\) for this model, numbering the steps (a) through (f), as shown in Activity 14.4. As a hint: you should get

\begin{align*} \mathcal{R}_0 \amp = \frac{\kappa \beta N}{(\mu + \kappa)(\mu + \gamma)} \end{align*}

as your result. Demonstrate fully how you complete each step of the Next Generation Method in order to achieve this result.

2.

Consider the SEIR model shown in Figure 14.10. In this exercise, analyze the ways each parameter within the \(\mathcal{R}_0\) formula affects the size of \(\mathcal{R}_0\text{.}\) (Note: the formula for \(\mathcal{R}_0\) for this model is given above in Exercise 1.) In particular, consider the following.

How do \(\beta\) and \(N\) affect \(\mathcal{R}_0\text{?}\) That is, as \(\beta\) increases, what happens to \(\mathcal{R}_0\text{?}\) As \(\beta\) decreases, what happens to \(\mathcal{R}_0\text{?}\) Similarly, how do increases or decreases in \(N\) affect \(\mathcal{R}_0\text{?}\)
How does \(\kappa\) affect \(\mathcal{R}_0\text{,}\) given that \(\kappa\) appears in both the numerator and denominator of \(\mathcal{R}_0\text{?}\)
Consider the fraction \(\frac{\kappa}{\mu + \kappa }\text{.}\) How large or small can this fraction become? How does this fraction affect \(\mathcal{R}_0\) as \(\kappa\) becomes larger or smaller, and how does this fraction affect \(\mathcal{R}_0\) as \(\mu\) becomes larger or smaller?
Consider the fraction \(\frac{1}{\mu + \gamma }\text{.}\) How does this fraction affect \(\mathcal{R}_0\) as \(\gamma\) becomes larger or smaller, and how does this fraction affect \(\mathcal{R}_0\) as \(\mu\) becomes larger or smaller?
If we set \(\mu = 0\text{,}\) we can simplify the value of \(\mathcal{R}_0 =\frac{\kappa \beta N}{(\mu + \kappa)(\mu + \gamma)}\text{.}\) What result do we get when we simplify in this way, and how can we interpret this result in relation to the SEIR model in this activity, as well as in relation to an SEIR model without birth and death terms, such as the model in Activity 7.8?

3.

Revisit the SIRS model from Activity 14.5 and introduce birth and death terms in a similar way as is shown in the model in Figure 14.10. Then use the steps of the Next Generation Method to compute \(\mathcal{R}_0\) for the resulting SIRS model with births and deaths, numbering the steps (a) through (f), as shown in Activity 14.4. Demonstrate fully how you complete each step of the Next Generation Method.

Prev Top Next