A Comparison of Methods for Calculating the Basic Reproductive Number for Periodic Epidemic Systems

Abstract

When using mathematics to study epidemics, oftentimes the goal is to determine when an infection can invade and persist within a population. The most common way to do so uses threshold quantities called reproductive numbers. An infection’s basic reproductive number (BRN), typically denoted $R_0$, measures the infection’s initial ability to reproduce in a naive population and is tied mathematically to the stability of the disease-free equilibrium. Next-generation methods have long been used to derive $R_0$ for autonomous continuous-time systems; however, many diseases exhibit seasonal behavior. Incorporating seasonality into models may improve accuracy in important ways, but the resulting non-autonomous systems are much more difficult to analyze. In the literature, two principal methods have been used to derive BRNs for periodic epidemic models. One, based on time-averages, is simple to apply but does not always describe the correct threshold behavior. The other, based on linear operator theory, is more general but also more complicated, and no detailed explanations of the necessary computations have yet been laid out. This paper reconciles the two methods by laying out an explicit procedure for the second and then identifying conditions (and some important classes of models) under which the two methods agree. This allows the use of the simpler method, which yields interpretable closed-form expressions, when appropriate, and illustrates in detail the simplest possible case where they disagree. Results show that seasonality alone cannot affect disease persistence, but must act in conjunction with non-hierarchical heterogeneity in the infected population, in order to do so.

Appendix A: Development of the Linear Operator Approach

1.1 A.1 Floquet Theory

Before developing the linear operator method for non-autonomous systems, some background is needed for periodic dynamical systems. Consider the system:

$$\begin{aligned} \dot{x}=A(t)x,\quad x \in \mathbb {R}^n \end{aligned}$$

(3)

where $t \rightarrow A(t)$ is a T-periodic continuous matrix-valued function. A fundamental matrix of a system of n homogeneous ordinary differential equations is a matrix-valued function, $\Psi (t)$, whose columns are linearly independent solutions of the system. Floquet’s theorem gives a canonical form for fundamental matrix solutions (Chicone 1999).

Theorem A.1

(Floquet’s Theorem) If $\Psi (t)$ is a fundamental matrix solution of the T-periodic system (3), then, for all $t \in \mathbb {R}$,

$$\begin{aligned} \Psi (t + T) = \Psi (t)\Psi ^{-1}(0)\Psi (T). \end{aligned}$$

In addition, there is a matrix B (which may be complex) such that

$$\begin{aligned} e^{TB} = \Psi ^{-1}(0)\Psi (T) \end{aligned}$$

and a T-periodic matrix function $t \rightarrow P(t)$ (which may be complex-valued) such that $\Psi (t) = P(t)e^{tB}$ for all $t \in \mathbb {R}$. Also, there is a real matrix R and a real 2T-periodic matrix function $t \rightarrow Q(t)$ such that $\Psi (t) = Q(t)e^{tR}$ for all $t \in \mathbb {R}$.

This theorem gives a representation called a Floquet normal form for the fundamental matrix $\Psi (t)$. For an $\omega$-periodic system define a monodromy matrix, $\Phi$, in terms of a fundamental matrix $\Psi (t)$ for the system, as $\Phi =\Psi (\omega )\Psi ^{-1}(t_0)$. A matrix $\Psi (t)$ is a principal fundamental matrix if it is a fundamental matrix and there exists a $t_0$ such that $\Psi ^{-1}(t_0)=I$. So if the solution to the $\omega$-periodic system is a principal fundamental matrix, then the monodromy matrix is just $\Phi =\Psi (\omega )$.

The problem becomes how these matrices are calculated. Most systems cannot be solved analytically, so numerical work must be done. The fundamental matrix is computed in terms of Taylor series, using the original matrix A(t) to various powers. This can be done since the solution is of the form $e^{TB} = \Psi ^{-1}(0)\Psi (T)$, by using the matrix exponential definition:

$$\begin{aligned} e^{tA} = \sum _{0}^{\infty }\frac{t^k}{k!}A^k =I+tA+\frac{t^2}{2!}A^2 +\frac{t^3}{3!}A^3 +\cdots . \end{aligned}$$

1.2 A. 2 Linear Operator Method

Consider a setting similar to van den Driessche and Watmough (2002) where a heterogeneous population is grouped into n homogeneous compartments. In this model, certain parameters will be assumed to be $\omega$-periodic, thus giving a non-autonomous system. As before, sort the compartments so that the first m compartments represent infected classes. Consider again the set $\mathbf X _s$ to be all disease-free states, where now values in $\mathbf X _s$ can be disease-free periodic solutions and not just equilibria. Let $\mathscr {F}_i(t,x)$ be the input rate of newly infected individuals in the ith compartment, $\mathscr {V}_i^+(t,x)$ be the input rate of individuals by other means, and $\mathscr {V}_i^-(t,x)$ be the rate of transfer out of compartment i. The model is then given by:

$$\begin{aligned} \frac{{\hbox {d}}x_i}{{\hbox {d}}t}=\mathscr {F}_i(t,x)-\mathscr {V}_i(t,x)=f_i(t,x),\quad i=1,\ldots ,n \end{aligned}$$

where $\mathscr {V}_i=\mathscr {V}_i^--\mathscr {V}_i^+$. As with the procedure for autonomous models, some assumptions must be made on the system, to ensure that the model is well posed and makes biological sense. Five of them are analogous to those imposed by van den Driessche and Watmough; two more are specific to periodic models.

• (A1) For each $1\le i\le n$, the functions $\mathscr {F}_i(t,x), \mathscr {V}_i^+(t,x),$ and $\mathscr {V}_i^-(t,x)$ are nonnegative and continuous on $\mathbb {R} \times \mathbb {R}_+^n$ and continuously differentiable with respect to x.

• (A2) There is a real number $\omega >0$ such that for each $1\le i\le n$, the functions $\mathscr {F}_i(t,x), \mathscr {V}_i^+(t,x),$ and $\mathscr {V}_i^-(t,x)$ are $\omega$-periodic in t. (This is new and establishes $\omega$-periodicity of the system.)

• (A3) If $x_i=0$, then $\mathscr {V}_i^-=0$. If $x \in \mathbf X _s$, then $\mathscr {V}_i^-=0$ for $i=1,\ldots m$.

• (A4) $\mathscr {F}_i=0$ if $i>m$.

• (A5) If $x \in \mathbf X _s$, then $\mathscr {F}_i=0$ and $\mathscr {V}_i^+=0$ for $i=1,\ldots m$.

The remaining two assumptions are more lengthy to state. The next one (A6) is analogous to (A5) from the autonomous case. Assume the model has a disease-free periodic solution $x_0(t)$. It must be verified that $x_0(t)$ is linearly asymptotically stable in the disease-free subspace, $\mathbf X _s$. Define an $(n-m)\times (n-m)$ matrix

$$\begin{aligned} M(t):=\left( \frac{\partial f_i(t,x_0(t))}{\partial x_j}\right) _{m+1\le i,j\le n}. \end{aligned}$$

Let $\Phi _M$ be the monodromy matrix of the linear $\omega$-periodic system $\frac{{\hbox {d}}z}{{\hbox {d}}t}=M(t)z$. To verify local stability, one must show that the spectral radius of the monodromy matrix is less than one, or (A6) $\rho (\Phi _M)<1$.

For the last assumption, (A7), we require that the infection eventually die out if no new infections arise. Following the notation of van den Driessche and Watmough (2002), define two $m\times m$ matrices

$$\begin{aligned} F(t)=\left[ \frac{\partial \mathscr {F}_i(t,x_0(t))}{\partial x_j}\right] _{1\le i,j\le m}, \quad V(t)=\left[ \frac{\partial \mathscr {V}_i(t,x_0(t))}{\partial x_j}(x_0)\right] _{1\le i,j\le m}. \end{aligned}$$

Let Y(t, s), $t\ge s$ be the evolution operator of the linear $\omega$-periodic system $\frac{{\hbox {d}}y}{{\hbox {d}}t}=-V(t)y$. That is, for each $s\in \mathbb {R}$ the $m\times m$ matrix Y(t, s) satisfies

$$\begin{aligned} \frac{{\hbox {d}}}{{\hbox {d}}t}Y(t,s)=-V(t)Y(t,s),\quad \forall t\ge s, \, Y(s,s)=I \end{aligned}$$

where I is the $m\times m$ identity matrix. The monodromy matrix $\Phi _{-V}$ of the system then equals $Y(\omega ,0)$. The last assumption imposed is that the internal evolution of individuals in the infectious compartments is dissipative and exponentially decays in many cases, i.e., that (A7) $\rho (\Phi _{-V})<1$.

Stability theory for equilibria of nonlinear autonomous systems has been well established. In 1960, the Hartman–Grobman theorem was proved, showing that the local behavior of the system around a hyperbolic equilibrium point can be studied using the linearized system around the same point (Aulbach and Wanner 2000). In 1975, this theorem was extended by Kenneth Palmer to non-autonomous systems (Palmer 1975). This allows the use of the linearized system in (A6) to be used to study the stability of the disease-free solution.

Taking assumptions (A1)–(A7) as given, the basic reproductive number for the epidemic model can then be calculated, considering perturbations about the disease-free periodic state $x_0(t)$. By standard theory of linear periodic systems, Hale (1980), there exist a $K>0$ and $\alpha >0$ such that

$$\begin{aligned} \Vert Y(t,s)\Vert \le Ke^{-\alpha (t-s)},\quad \forall t\ge s, \, s\in \mathbb {R}. \end{aligned}$$

It then follows that

$$\begin{aligned} \Vert Y(t,t-a)F(t-a)\Vert \le K\Vert F(t-a)\Vert e^{-\alpha a}, \quad \forall t\in \mathbb {R}, \, a\in [0,\infty ). \end{aligned}$$

Given the $\omega$-periodic system, let $\phi (s)$ denote the ($\omega$-periodic in s) distribution vector of infectious individuals. Then $F(s)\phi (s)$ is the distribution vector of individuals who were newly infected at time s. Given $t\ge s$, then $Y(t,s)F(s)\phi (s)$ gives the distribution of those infected individuals who were newly infected at time s and remain infected at time t. Then

$$\begin{aligned} \psi (t):=\int _{-\infty }^{t}Y(t,s)F(s)\phi (s){\hbox {d}}s=\int _{0}^{\infty }Y(t,t-a)F(t-a)\phi (t-a){\hbox {d}}a \end{aligned}$$

is the distribution of accumulative new infections at time t produced by all those infected individuals introduced at times prior to t.

Let $C_{\omega }$ be the ordered Banach space of all $\omega$-periodic functions from $\mathbb {R}$ to $\mathbb {R}^m$, with the max norm and the positive cone $C_{\omega }^+:=\{\phi \in C_{\omega }: \phi (t)\ge 0, \forall t\in \mathbb {R}\}$. Now define a new linear operator $L:C_{\omega } \rightarrow C_{\omega }$ by

$$\begin{aligned} (L\phi )(t)=\int _{0}^{\infty }Y(t,t-a)F(t-a)\phi (t-a){\hbox {d}}a, \quad \forall t\in \mathbb {R},\, \phi \in C_{\omega }. \end{aligned}$$

Call L the next infection operator, following the motivation of van den Driessche and Watmough, and then the spectral radius of L is given by:

$$\begin{aligned} R_{LO}:=\rho (L) \end{aligned}$$

for the periodic epidemic model.

Following the approach from Bacaër and Guernaoui (2006), one can obtain another linear operator on $C_{\omega }$, using the same notation as Wang and Zhao (2008):

$$\begin{aligned} (\bar{L}\phi )(t)= & {} \int _{0}^{\infty }\!\!\!\!F(t)Y(t,t-a)\phi (t-a){\hbox {d}}a\\= & {} F(t)\!\int _{0}^{\infty }\!\!\!Y(t,t-a)\phi (t-a){\hbox {d}}a, \quad \forall t\in \mathbb {R},\, \phi \in C_{\omega }. \end{aligned}$$

The spectral radius of $\bar{L}$, $\rho (\bar{L})$, was defined in Bacaër and Guernaoui (2006) as the basic reproductive number. Wang and Zhao showed that the basic reproductive number for L and the basic reproductive number for $\bar{L}$ coincide, but that the kernels have different biological interpretations.

As before with autonomous systems, the question becomes whether the quantity thus defined identifies the usual threshold behavior of disease invasion, where if $R_0<1$ then the disease-free periodic solution is stable, while it is unstable if $R_0>1$. Wang and Zhao proved that indeed this is the case, but in order to characterize and calculate $R_0$ for periodic systems they needed to take a different approach to actually calculate the BRN. The reason for this is that linear operators are difficult to work with, so they decided to approach it using Floquet theory instead.

Consider the linear $\omega$-periodic equation

$$\begin{aligned} \frac{{\hbox {d}}w}{{\hbox {d}}t}=\left[ -V(t)+\frac{F(t)}{\lambda }\right] w, \quad t\in \mathbb {R} \end{aligned}$$

(4)

with parameter $\lambda \in (0,\infty )$. Let $W(t,s,\lambda ), t\ge s, s\in \mathbb {R}$ be the evolution operator of the system (4) on $\mathbb {R}^m$. Wang and Zhao showed that the linear operator $W(t,s,\lambda )$ is positive in $\mathbb {R}^m$ for each $t\ge s, s\in \mathbb {R}$. The Perron–Frobenius theorem implies that $\rho (W(\omega , 0, \lambda ))$ is an eigenvalue of $W(\omega , 0, \lambda )$ with a nonnegative eigenvector. The following theorems from Wang and Zhao (2008) connect $R_0$ first to W and then to the monodromy matrix $\Phi _{F-V}$ for (4):

Theorem

(2.1)

• i. If $\rho (W(\omega , 0, \lambda ))=1$ has a positive solution $\lambda _0$, then $\lambda _0$ is an eigenvalue of L, and hence $R_0>0$.

• ii. If $R_0>0$, then $\lambda =R_0$ is the unique solution of $\rho (W(\omega , 0, \lambda ))=1$.

• iii. $R_0=0$ if and only if $\rho (W(\omega , 0, \lambda ))<1$ $\forall \lambda >0$.

Theorem

(2.2)

• i. $R_0=1$ if and only if $\rho (\Phi _{F-V}(\omega ))=1$

• ii. $R_0>1$ if and only if $\rho (\Phi _{F-V}(\omega ))>1$

• iii. $R_0<1$ if and only if $\rho (\Phi _{F-V}(\omega ))<1$

Thus the disease-free solution, $x_0(t)$, is asymptotically stable if $R_0<1$ and unstable if $R_0>1$.

These results show that to find the basic reproductive number, one needs to find the monodromy matrix $\Phi _{F-V}(\lambda )$ of the system (4) as a function of the parameter $\lambda$. Then find the spectral radius of $\Phi _{F-V}(\lambda )$ and solve the equation $\rho (\Phi _{F-V}(\lambda ))=1$ for $\lambda$. This value of $\lambda$ is the basic reproductive number $R_0$. The threshold behavior for the disease-free solution can then be studied.