Computation of the basic reproduction numbers for reaction-diffusion epidemic models

1.   Introduction

Partial differential equations (PDEs) of the reaction-diffusion type are extensively used in the modeling of infectious diseases ^{[1,2,3,4,5,6,7,8]}. The focus of the present paper is to study the basic reproduction numbers for a class of reaction-diffusion epidemic models constructed from autonomous systems of ordinary differential equations (ODEs). The underlying ODE models represent the dynamics of disease transmission and spread that are spatially homogeneous, whereas the PDE models, with diffusion terms added, emphasize the movement and dispersal of the hosts and pathogens over a (typically heterogeneous) spatial domain.

The basic reproduction number, commonly denoted by $\mathcal{R}_0$, is a critical quantity to quantify the transmission risk of an infectious disease. It measures the expected number of secondary infections produced by one infective individual in a completely susceptible population. $\mathcal{R}_0$ is often used to characterize the threshold behavior of an epidemic, with disease eradication if $\mathcal{R}_0 < 1$ and disease persistence if $\mathcal{R}_0 > 1$. The theory of the basic reproduction numbers for ODE-based autonomous epidemic models is well developed, and the calculation of $\mathcal{R}_0$ follows a standard procedure based on the next-generation matrix technique ^[9,10].

Many efforts have been devoted to the definition, calculation, and analysis of the basic reproduction numbers for PDE models, such as reaction-diffusion epidemic systems, which are more complex than their ODE counterparts. Thieme ^[11] introduced a theoretical framework to study $\mathcal{R}_0$, defined as the spectral radius of a resolvent-positive operator, of such reaction-diffusion equations. Using the theory of principal eigenvalues, Wang and Zhao ^[12] defined $\mathcal{R}_0$ as the spectral radius of a next-infection operator for a general class of reaction-diffusion models. They showed that the $\mathcal{R}_0$ of such a PDE system is the same as that of the underlying ODE system when the diffusion rates are positive constants and the next-generation matrices are independent of the spatial location. Asymptotic profiles of $\mathcal{R}_0$ for reaction-diffusion epidemic systems with constant diffusion rates were investigated in ^{[13,14,15,16]}. Particularly, Chen and Shi ^[14] showed that $\mathcal{R}_0$ approaches the spectral radius of a spatially averaged next-generation matrix when the diffusion rates tend to infinity, and $\mathcal{R}_0$ approaches the maximum value of a local reproduction number when the diffusion rates tend to zero. Moreover, reproduction numbers for time-periodic reaction-diffusion systems were discussed in ^[17,18]. Other theoretical studies related to the reaction-diffusion epidemic models and their basic reproduction numbers include ^{[1,6,8,19,20,21,22,23]}.

The goal of the present work is twofold: (1) to produce practically useful means to compute the basic reproduction numbers of reaction-diffusion epidemic models, since the classical next-generation matrix technique for autonomous ODE systems is no longer applicable; and (2) to gain a deeper understanding of the relationship between the basic reproduction number of a reaction-diffusion PDE model, $\mathcal{R}_0^{{\text{PDE}}}$, and that of its underlying ODE model, $\mathcal{R}_0^{{\text{ODE}}}$. To that end, we will focus on a class of reaction-diffusion epidemic systems that are developed by adding diffusion terms to autonomous ODE systems, where the diffusion rates generally are functions of the location variables {representing the spatial heterogeneity.} In a prior study ^[24], we proposed a numerical method to compute the value of $\mathcal{R}_0$ for such reaction-diffusion epidemic models on one-dimensional (1D) spatial domains. The essential idea is the reduction of the infinite-dimensional operator eigenvalue problem for a PDE system to a finite-dimensional matrix eigenvalue problem.

In the present study, we will make a nontrivial extension of the methodology in ^[24] to reaction-diffusion epidemic systems on $k$-dimensional spatial domains, where $k \geq 1$ can be any positive integer. Such an extension, {in addition to making the theory and methodology more complete, would facilitate the study of more practical applications} where these PDE models are utilized to investigate the transmission and spread of infectious diseases {in the real world.} Based on the numerical formulation, we will use matrix theory to analyze the relationship between the PDE-based $\mathcal{R}_0^{{\text{PDE}}}$ and the ODE-based $\mathcal{R}_0^{{\text{ODE}}}$. The matrix analysis involved in the current work for $k$-dimensional models is significantly harder than that for one-dimensional models. We will show that under several important scenarios, such as the presence of a single infected compartment, constant diffusion rates, uniform diffusion of the infected compartments, and partial diffusion in a system, the two basic reproduction numbers equal each other.

We organize the remainder of this paper as follows. In Section 2, we present the $k$-dimensional reaction-diffusion epidemic system and the definition of its basic reproduction number. In Section 3, we describe the details of our computational method for $\mathcal{R}_0^{{\text{PDE}}}$ and then analyze the relationship between $\mathcal{R}_0^{{\text{PDE}}}$ and $\mathcal{R}_0^{{\text{ODE}}}$. In Section 4, we provide specific numerical examples to verify our analytical findings. Finally, we conclude the paper with some discussion in Section 5.

2.   Reaction-diffusion epidemic model

Let $n$ be a positive integer and $U(t, x) = \big(u_1(t, x), ..., u_n(t, x) \big)^T$ be a vector-valued function that represents the hosts and pathogens related to an infectious disease, with each $u_i(t, x)$ denoting the density of the (host or pathogen) population in compartment $i$ ($1\leq i \leq n$) at time $t$ and location $x$. We are concerned with the $k$-dimensional spatial domain $[0, 1]^k$, where $k\ge1$ is an integer. We consider the following reaction-diffusion epidemic system

with appropriate initial conditions. In this model, $d_i(x)$ ($1\le i\le n$) denotes the diffusion rate at location $x$ and is assumed to be continuously differentiable on $[0, 1]^k$. $\mathcal{F}_i(U)$ denotes the rate of generation for newly infected individuals in compartment $i$, and $\mathcal{V}_i(U) = \mathcal{V}^-_i(U)-\mathcal{V}^+_i(U)$, with $\mathcal{V}_i^+$ denoting the transfer rate of individuals into compartment $i$ and $\mathcal{V}_i^-$ the transfer rate of individuals out of compartment $i$. Note that $\mathcal{F}_i(U)$ and $\mathcal{V}_i(U)$, $i = 1, 2, \cdots, n$, are functions of $U$ only, so that the PDE model (2.1) is associated with an underlying ODE model, discussed in Appendix A. In addition, $\nu$ is the unit normal vector on the boundary $\partial[0, 1]^k$.

Without loss of generality, we assume that $U_I = U_I(t, x) = (u_1, ..., u_m)^T$ denotes all the infected compartments in the vector $U$, where $1 \leq m < n$. Consequently, the set of all disease-free steady states is defined as $U_s = \{U\ge0 : u_i = 0, \; i = 1, ..., m\}.$

System (2.1) can be re-written as

where

In what follows, we investigate the PDE system (2.2), where our results can be easily applied to the original system (2.1) under the condition (2.3).

Following the theoretical framework in ^[12], we let $T(t)$ be the solution semigroup on $C([0, 1]^k, \mathbb{R}^m)$ associated with the following linear reaction-diffusion equation:

Let the distribution of the initial infections, i.e., $U_I(0, x)$, be $U_m(x) = (u_1(x), ..., u_m(x))^T$. Then the distribution of these infections after time $t > 0$ is given by $T(t)(U_m(x))$. Let $F$ be the generation matrix of new infections (see Appendix A). Then the distribution of new infections at any time $t > 0$ is $FT(t)(U_m(x))$ and the distribution of the total new infections is represented by $\int_0^{+\infty}FT(t)(U_m(x))dt.$ Therefore, the next-generation operator $L$, which maps the distribution of initial infections to the distribution of the total infective individuals generated during the infectious period, is defined by

Consequently, the basic reproduction number for the PDE system (2.2) is the spectral radius of the operator $L$:

Meanwhile, we introduce the operator $\Gamma: \; \mathbb{R}^m\longmapsto\mathbb{R}^m$ by

where, for $U = (u_1, \cdots, u_m)^T \in \mathbb{R}^m$, $(D_m(x)\Delta)U = \big(d_1(x)\Delta u_1, \, ..., \, d_m(x)\Delta u_m \big)^T$ and $(C_m(x)\cdot\nabla)U = \left(c_1(x)\cdot\nabla u_1, \, ..., \, c_m(x)\cdot\nabla u_m\right)^T$, with $c_j(x) = (c_{j1}(x), c_{j2}(x), ..., c_{jk}(x))$, $1\le j \le m$. Then we can obtain an essential characterization of the next-infection operator,

with details provided in Appendix B.

Below we focus our attention on the numerical calculation of $\mathcal{R}_0^{{\text{PDE}}}$, based on Eqs (2.6) and (2.8), and the analysis of its relationship with the basic reproduction number $\mathcal{R}_0^{{\text{ODE}}}$ of the underlying ODE model (A.1). Throughout our discussion, we assume that the conditions (A1)–(A4) in Appendix A and (B1)–(B3) in Appendix B hold.

3.   Numerical formulation and matrix analysis

Let $\lambda$ be an eigenvalue of the operator $L$ such that $L(\phi(x)) = \lambda\phi(x)$ for an eigenvector $\phi(x) = (\phi_1(x), ..., \phi_m(x))^T$. Then from Eq (2.8) we have

Let $\psi(x) = -\Gamma^{-1}(\phi(x))$, where $\psi(x) = (\psi_1(x), ..., \psi_m(x))^T$. Then $-\Gamma(\psi(x)) = \phi(x)$. Based on the condition (B3), this equation can be written as

Pick a sufficiently large integer $N > 0$ and denote

for any integers $0\le j_1, j_2, ..., j_k\le N$ and $1\le r\le k$. Apply the standard centered difference scheme to Eq (3.2) on the spatial domain $[0, \, 1]^k$. Then for any $0\le j_1, ..., j_{k}\le N$, we obtain

and $\psi_p^{j_1...(-1)...j_{k}} = \psi_p^{j_1...(1)...j_{k}}, \; \psi_p^{j_1...(N+1)...j_{k}} = \psi_p^{j_1...(N-1)...j_{k}}$, where $(-1), (1), (N+1)$, and $(N-1)$ are all in the $(j_r)$ position inside the permutation $j_1...j_k$ for any $1\le r\le k$ by the Neumann boundary condition. Next, we can write the above $(N+1)^k$ approximate equations of (3.3) in the following matrix form

where $A_p = \left(a^{(p)}_{ij}\right)$ is a $(N+1)^k\times(N+1)^k$ matrix, $1\le i, j\le(N+1)^k$, and

Note that for any $0 \le j_1, ..., j_k\le N$, the coefficient of $\psi_p^{j_1...j_{k}}$ in Eq (3.3) is a diagonal entry of $A_p$, which is equal to a positive number $2kN^2d_p^{j_1...j_k}+v_p$. Define

and where $d_0$ is a positive lower bound for the diffusion rates (see condition (B2)). Then for $N > N_0$, the off-diagonal entries $-N^2d_p^{j_1...j_{k}}+\frac12Nc_{pr}^{j_1...j_{k}}$ and $-N^2d_p^{j_1...j_{k}}-\frac12Nc_{pr}^{j_1...j_{k}}$ are nonpositive for all $1\le r\le k$. Hence for any eigenvalue $\lambda$ of $A_p$, by {the Gershgorin Circle Theorem}, there exists $1\le i\le (N+1)^k$ such that

Thus, Re$(\lambda)\ge v_p > 0$ and $A_p$ is invertible. Moreover, we have $\rho({A_p^{-1}})\le1/$Re$(\lambda)\le1/v_p, \; 1\le p\le m.$ This leads to the following lemma.

Lemma 3.1. Let $N > N_0$ and $\lambda_{A_p}$ be an eigenvalue of matrix $A_p$. Then, for $1\le p\le m$, the real part of $\lambda_{A_p}$ satisfies Re$(\lambda_{A_p})\ge v_p$, and, consequently, $A_p$ is invertible.

In addition, for any $0 \le j_1, ..., j_k\le N$, if we fix $\psi_p^{j_1...(j_r+1)...j_{k}} = \psi_p^{j_1...j_{k}} = \psi_p^{j_1...(j_r-1)...j_{k}} = 1$ for all $1\le r\le k$ in Eq (3.3), then the left-hand side of Eq (3.3) is the sum of the $\left(1+ \sum\limits_{r = 1}^kj_r(N+1)^{k-r}\right)$-th row of the matrix $A_p$, which is obviously equal to $v_p$. Hence the sum of each row of $A_p$ is $v_p$, which implies $v_p$ is an eigenvalue of $A_p$, and consequently, $1/v_p$ is an eigenvalue of $A_p^{-1}$. Therefore $\rho({A_p^{-1}}) = 1/v_p$. We obtain the following result.

Lemma 3.2. For all $N > N_0$, we have $\rho({A_p^{-1}}) = 1/v_p, \; 1\le p\le m$.

Denote $\Psi = (\Psi_1^T, ..., \Psi_m^T)^T$, $\Phi = (\Phi_1^T, ..., \Phi_m^T)^T$, and

Then $A$ is invertible and $\Psi\approx A^{-1}\Phi$ by Eq (3.4). It follows from Eq (3.1) that

which yields

for any integer $N > 0$, where $I_{(N+1)^k}$ is the $(N+1)^k\times(N+1)^k$ identity matrix and $\otimes$ denotes the Kronecker product that is defined as follows: for any $r\times s$ matrix $M = (m_{ij})$ and $p\times q$ matrix $Q$,

With the substitution of $\Psi\approx A^{-1}\Phi$ into Eq (3.6), our numerical formulation leads to

From the basic theory of finite difference schemes ^[25,26], the solution of Eq (3.7) (or, equivalently, Eq (3.3)) converges to the solution of Eq (3.1) (or, equivalently, Eq (3.2)) when $N \to \infty$. Hence, for any $\varepsilon > 0$, we can pick $N$ sufficiently large such that

Letting $\varepsilon \to 0$, we obtain our central result for the computation of $\mathcal{R}_0^{{\text{PDE}}}$:

We have reduced the original operator eigenvalue problem (3.1) to a matrix eigenvalue problem (3.7). Since there are many efficient numerical techniques available for computing eigenvalues of matrices ^[27,28], our method facilitates practical evaluation of the basic reproduction number for such a reaction-diffusion epidemic model.

Additionally, our numerical formulation provides important insight into the property of $\mathcal{R}_0^{{\text{PDE}}}$. In what follows, we apply matrix theory to conduct an analysis of $\mathcal{R}_0^{{\text{PDE}}}$ and its connection to $\mathcal{R}_0^{{\text{ODE}}}$, based on Eq (3.8). We first introduce the following lemma.

Lemma 3.3. Assume that $X = (x_{ij})$ is an $m\times m$ matrix and $Y_{ij} \; (1\le i, j\le m)$ are $n\times n$ matrices. If there exists a nonsingular matrix $P$ such that $P^{-1} Y_{ij}P = U_{ij}$ for all $i, j = 1, ..., m$, where $U_{ij}$ is an upper triangular matrix with diagonal entries $y^{(1)}_{ij}, ..., y^{(n)}_{ij}$, then

The proof of Lemma 3.3 is similar to that of Lemma 4.2 in ^[24]. Now we state our main results regarding the relationship between $\mathcal{R}_0^{{\text{PDE}}}$ and $\mathcal{R}_0^{{\text{ODE}}}$ in the following three theorems.

Theorem 3.1. (1) In general, we have $\mathcal{R}_0^{PDE}\ge\mathcal{R}_0^{ODE}$.

(2) If $F$ is a triangular matrix, then $\mathcal{R}_0^{PDE} = \mathcal{R}_0^{ODE}$. Particularly, if $m = 1$, then $\mathcal{R}_0^{PDE} = \mathcal{R}_0^{ODE}$.

Proof. (1) Let $e = (1, 1..., 1)^T$ be a vector with all the $(N+1)^k-1$ entries being 1's, and let $P = \begin{bmatrix}1 & 0\\ e &I_{(N+1)^k-1}\end{bmatrix}$. Since the sum of each row of $A_i$ is $v_i$, we have

where $\alpha_i$ is a $[(N+1)^k-1]$-dimensional vector and $S_i$ is a $[(N+1)^k-1]\times[(N+1)^k-1]$ matrix with $(N+1)^k-1$ rows and $(N+1)^k-1$ columns. Similar to the proof of Lemma 4.2 in ^[24], we obtain that $\det(\lambda I_m-FV^{-1})$ is a factor of $\det\big[\lambda I_{m(N+1)^k}-(F\otimes I_{(N+1)^k})A^{-1}\big]$; that is, each eigenvalue of $FV^{-1}$ is an eigenvalue of $(F\otimes I_{(N+1)^k})A^{-1}$. Thus, $\rho((F\otimes I_{(N+1)^k})A^{-1})\ge\rho(FV^{-1})$. Letting $N\to\infty$, we obtain $\mathcal{R}_0^{{\text{PDE}}}\ge\mathcal{R}_0^{{\text{ODE}}}$.

(2) This statement directly follows from Lemma 3.2, since

□

Theorem 3.2. If the matrix set $\{A_i\}_{i = 1}^{m}$ for system (2.2) is a commuting family where each pair of matrices commute with each other, then $\mathcal{R}_0^{PDE} = \mathcal{R}_0^{ODE}$.

Proof. By Theorem 3.1(1), it suffices to show that $\mathcal{R}_0^{{\text{PDE}}}\le\mathcal{R}_0^{{\text{ODE}}}$. Since $\{A_i\}_{i = 1}^{m}$ is a commuting family of matrices, then $\{A_i^{-1}\}_{i = 1}^{m}$ is a commuting family. Hence, there exists a nonsingular matrix $Q$ such that

where $B_i$ is an upper triangular matrix with diagonal elements $\alpha_{i1}, \; ..., \; \alpha_{i (N+1)^k}$ and $|\alpha_{ij}|\le 1/v_i$ for $i = 1, ..., m$; $j = 1, ..., (N+1)^k$. Consequently, by Lemma 3.3, we have

where $O_j = \begin{bmatrix} F_{11}\alpha_{1j} & \cdots & F_{1m}\alpha_{mj} \\ \vdots & \vdots & \vdots \\ F_{m1}\alpha_{1j} & \cdots & F_{mm}\alpha_{mj} \end{bmatrix}$ for $1\le j\le (N+1)^k$. Thus, Eq (3.9) indicates that

Since $F$ is nonnegative by assumptions (A1) and (A4), then $|O_j|\le FV^{-1}$, where $|O_j| = \begin{bmatrix} |F_{11}\alpha_{1j}| & \cdots & |F_{1m}\alpha_{mj}| \\ \vdots & \vdots & \vdots \\ |F_{m1}\alpha_{1j}| & \cdots & |F_{mm}\alpha_{mj}| \end{bmatrix}$. We thus obtain $\; \rho(O_j)\le\rho(|O_j|)\le\rho(FV^{-1}).$ Therefore, $\rho((F\otimes I_{(N+1)^k})A^{-1})\le\rho(FV^{-1})$, which implies $\mathcal{R}_0^{{\text{PDE}}}\le\mathcal{R}_0^{{\text{ODE}}}$. □

Next, we provide sufficient and necessary conditions to characterize the scenarios where $\{A_p\}_{p = 1}^{m}$ is a commuting family.

Theorem 3.3. For any integer $N > N_0$, the matrix set $\{A_p\}_{p = 1}^{m}$ associated with system (2.2) is a commuting family if and only if there exist constants $\delta_p, \, \sigma_{pr}$, and continuous functions $d(x), \, g_r(x), \, 1\le p\le m, \, 1\le r\le k$, such that

and

for any $1\le p, q\le m$ and $1\le r, i, j\le k$.

Proof. For any $0\le j_1, ..., j_k\le N$, we denote the $\left(1+ \sum_{r = 1}^kj_r(N+1)^{k-r}\right)$-th entry of a $(N+1)^k$-dimensional vector $\beta$ by $(\beta)_{j_1...j_k}$ and write $A_p = N^2H_p+ \frac{N}2G_p+v_pI_{(N+1)^k}$ for $p = 1, ..., m$, where $H_p$ satisfies

and $G_p$ satisfies

Thus, for any $1\le p, q\le m$, the equality

is equivalent to

Since Eq (3.11) holds for any $N > N_0$, we can conclude that

for any $1\le p, q\le m$. Following similar algebraic manipulations as those in ^[24], we can conclude the following: (i) $H_pH_q = H_qH_p$ implies that there exist constants $\delta_p$, $1\le p\le m$, and a continuous function $d(x)$ such that

(ii) $G_pG_q = G_qG_p$ implies that there exist constants $\sigma_{pr}, \; 1\le p\le m, \; 1\le r\le k$, and continuous functions $g_r(x), \; 1\le r\le k$, such that

for any $1\le p, q\le m, \; 1\le r, i, j\le k$. Substitute functions $d_p(x) = \delta_pd(x), \; c_{pr}(x) = \sigma_{pr}g_r(x)$ into matrices $H_p, \, H_q, \, G_p$ and $G_q$. Then for any $1\le p, q\le m$, we can obtain that $H_pG_q+G_pH_q = H_qG_p+G_qH_p$ holds if and only if

for all $1 \le r \le k$. We thus complete the proof. □

The conclusions in Theorems 3.2 and 3.3 can easily apply to the original PDE system (2.1) under the constraint (2.3). In each of the following scenarios, the basic reproduction number $\mathcal{R}_0^{{\text{PDE}}}$ for the reaction-diffusion system (2.1) and the basic reproduction number $\mathcal{R}_0^{{\text{ODE}}}$ for its ODE counterpart (A.1) are the same. These results cover several special, but important, cases associated with the PDE model (2.1).

Corollary 3.1. If there exist constants $\delta_p$, $1\leq p \leq m$, and a continuous function $d(x)$ such that $d_p(x) = \delta_pd(x)$ for all $p = 1, \, ..., \, m$ in system (2.1), then $\mathcal{R}_0^{PDE} = \mathcal{R}_0^{ODE}$.

Corollary 3.2. (A scenario with constant diffusion rates.) If the diffusion rates of all the infected compartments are positive constants in system (2.1), then $\mathcal{R}_0^{PDE} = \mathcal{R}_0^{ODE}$.

Corollary 3.3. (A scenario with uniform diffusion patterns.) If $d_p(x) = d_q(x)$ for all the infected compartments ($1\le p, \, q \le m$) in system (2.1), then $\mathcal{R}_0^{PDE} = \mathcal{R}_0^{ODE}$.

Corollary 3.4. (A scenario with partial diffusion.) If $d_p(x) = 0$ for $p = 1,\, ...,\, m-1$ and $d_m(x) \ge d_0 > 0$ in system (2.1), then $\mathcal{R}_0^{PDE} = \mathcal{R}_0^{ODE}$.

4.   Two examples

Several numerical examples concerned with one-dimensional (1D) reaction-diffusion epidemic models were presented in ^[24] to demonstrate that they have the same basic reproduction numbers as those of their ODE counterparts. Now we extend the numerical studies to two-dimensional (2D) and three-dimensional (3D) spatial domains to verify some of our analytical predictions in Section 3.

4.1.   A 2D SIR model

We consider a host population that moves on a 2D spatial domain represented by $[0, 1]^2$, where the motion can be described by a diffusion process. Let $S$, $I$ and $R$ be the density of the susceptible, infected, and recovered individuals, respectively, and $d_S(x)$, $d_I(x)$ and $d_R(x)$ be their associated diffusion rates with $x = (x_1, \, x_2) \in[0, 1]^2$. We study the following 2D reaction-diffusion SIR system, which is extended from the model presented in ^[4]:

The constant parameters $\Lambda$, $\alpha$, $\mu$, and $\gamma$ denote the recruitment rate, transmission rate, natural death rate, and disease recovery rate, respectively. Disease-induce mortality is not included here. Neumann boundary conditions are imposed on the boundary $\partial[0, 1]^2$ and appropriate initial conditions are provided at $t = 0$.

Obviously, $I$ is the only infected compartment in system (4.1); i.e., $m = 1$, so that Theorem 3.1(2) applies. From the underlying ODE system, we obtain $F = \frac{\alpha\Lambda}{\mu}$ and $V = \mu+\gamma$. Then the basic reproduction number of the PDE system (4.1) is the same as that of its corresponding ODE system, based on Theorem 3.1(2):

To verify this relationship, Figure 1 compares $\mathcal{R}_0^{{\text{PDE}}}$ and $\mathcal{R}_0^{{\text{ODE}}}$ for this model. $\mathcal{R}_0^{{\text{PDE}}}$ is computed by our numerical method based on Eq (3.8). The values of $\rho((F\otimes I_{(N+1)^2})A^{-1}) = \rho(\frac{\alpha\Lambda}{\mu}A_1^{-1})$ versus $N$ ($N = 1, 2, \cdots$) are plotted in Figure 1, where $A_1$ is the matrix obtained in Eq (3.4) from the single infectious compartment $I$. Meanwhile, since $\mathcal{R}_0^{{\text{ODE}}} = \rho(FV^{-1})$ does not depend on $N$, it is represented by a horizontal line in the graph. We set the diffusion rate of the infected individuals as $d_I(x) = \sin(100(x_1+x_2))+2$ in this test. We observe that when $N$ is sufficiently large, the numerical values of $\mathcal{R}_0^{{\text{PDE}}}$ based on $\rho((F\otimes I_{(N+1)^2})A^{-1})$ almost perfectly match $\mathcal{R}_0^{{\text{ODE}}}$, and this pattern continues for all $N \ge 40$.

Next, we verify that $\mathcal{R}_0^{{\text{PDE}}} = 1$ can be used as a threshold to distinguish the two dynamical behaviors between disease eradication and disease persistence for the model (4.1). To that end, we consider two typical scenarios of $\mathcal{R}_0^{{\text{PDE}}} < 1$ and $\mathcal{R}_0^{{\text{PDE}}} > 1$ by selecting appropriate parameter values, and run the simulation for the model (4.1) with an initial infection density $I(0, x_1, x_2) = 200$ in each scenario. Figure 2(a) displays the surface plot of $I(t, \, x_1, \, x_2 = 0.5)$ versus $t$ and $x_1$ with $\mathcal{R}_0^{{\text{PDE}}} = 0.88 < 1$, where $I$ approaches $0$ over time for all $x_1$. In contrast, Figure 2(b) displays the surface plot of $I(t, \, x_1, \, x_2 = 0.5)$ versus $t$ and $x_1$ with $\mathcal{R}_0^{{\text{PDE}}} = 3.20 > 1$, where $I$ increases from its initial value and remains positive for all the time. Surface plots at other fixed values of $x_2$ and those with $x_1$ fixed (not shown here) are qualitatively similar.

4.2.   A 3D model for environmentally transmitted diseases

Environmentally transmitted diseases continue to impose a significant public health burden throughout the world ^[29]. The transmission dynamics of many such diseases can be studied by SIR-B models ^[30,31,32], where the B compartment typically represents the concentration of the pathogen in the contaminated environment. Here, we focus on airborne infections, where the pathogenic particles (especially those tiny particles such as aerosols) can float and move in the air for an extended period of time. These airborne pathogens include bacteria such as Mycobacterium, Staphylococcus, and Legionella ^[33], viruses such as Varicella, Hantavirus, and SARS-CoV-2 ^[34], and other microorganisms such as fungi ^[35].

We consider {a reaction-diffusion SIR-B model} that incorporates both the indirect (i.e., airborne) and direct (i.e., human-to-human) transmission routes for an airborne infection. We assume that the pathogen undergoes a diffusion process in the air in a 3D domain represented by $[0, 1]^3$. We also assume that, compared to the pathogen diffusion and dispersal, the average spatial movement of human hosts is slow and can be disregarded in our model. We thus obtain {the following} PDE system

for $t > 0$ and $x = (x_1, \, x_2, \, x_3) \in [0, 1]^3$, with Neumann boundary conditions and appropriate initial conditions. The parameters $\alpha$ and $\beta$ represent, respectively, the direct and indirect transmission rates, $\xi$ is the rate of contribution from infected individuals (through coughing, sneezing, etc.) to the pathogen in the air, $r$ is the intrinsic growth rate of the pathogen (for viruses, we may set $r = 0$), $K$ is the carrying capacity of the pathogen growth, and $\tau$ is the removal rate of the pathogen from the air. A bilinear incidence form is used to represent both the direct and indirect transmission routes.

The model (4.2) is a partially diffusive PDE system since the diffusion component is only incorporated into the pathogen equation. Therefore, Corollary 3.4 predicts that $\mathcal{R}_0^{{\text{PDE}}}$ of system (4.2) equals $\mathcal{R}_0^{{\text{ODE}}}$ of the underlying ODE system. The infectious compartments are $I$ and $B$. From the associated ODE system, we easily obtain

From Corollary 3.4, we obtain

To provide numerical evidence for this analytical relationship, we compare $\mathcal{R}_0^{{\text{PDE}}}$ and $\mathcal{R}_0^{{\text{ODE}}}$ in Figure 3, where we plot $\rho((F\otimes I_{(N+1)^3})A^{-1})$ versus $N$ ($N = 1, 2, \cdots$) for the model (4.2). Note that $(F\otimes I_{(N+1)^3})A^{-1} = \begin{pmatrix} \frac{\alpha\Lambda}{\mu(\mu+\gamma)}I_{(N+1)^3} & \frac{\beta\Lambda}{\mu}A_2^{-1}\\ \frac{\xi}{\mu+\gamma}I_{(N+1)^3} & rA_2^{-1} \end{pmatrix}$ and $A_2$ is the matrix obtained in Eq (3.4) from the infectious compartment $B$. In order to clearly show the convergence of the numerical values of $R_0^{{\text{PDE}}}$, we set the pathogen diffusion rate as $d_B(x) = 10^{-3}\big(\sin(100(x_1+x_2+x_3))+1.01 \big)$ in this test. We observe a pattern similar to that in Figure 1. Specifically, the numerical approximations of $\mathcal{R}_0^{{\text{PDE}}}$ based on $\rho((F\otimes I_{(N+1)^3})A^{-1})$ coincide with $\mathcal{R}_0^{{\text{ODE}}}$ for all $N > 10$.

The stability properties for a more general partially diffusive SIR-B model were analyzed in ^[36] and it was shown that $\mathcal{R}_0^{{\text{PDE}}} = 1$ provided a threshold for the transition between disease eradication and disease persistence. We now provide some numerical verification for the PDE model (4.2). We again consider two typical scenarios with $\mathcal{R}_0^{{\text{PDE}}} < 1$ and $\mathcal{R}_0^{{\text{PDE}}} > 1$, and set the initial concentration of the pathogen as $B(0, x_1, x_2, x_3) = 10^5 \big(3 - (x_1 - 0.5)^2 - (x_2 - 0.5)^2 - (x_3 - 0.5)^2 \big)$ for the model (4.2). We then run the simulation in each scenario and plot the pathogen concentration $B(t, \, x_1, \, x_2 = 0.5, \, x_3 = 0.5)$ versus $t$ and $x_1$ in Figure 4, where we clearly see the eradication of the pathogen in panel (a) and the persistence of the infection in panel (b). Other surface plots with various values of $x_1$, $x_2$, and $x_3$ have similar behaviors and are not shown here.

5.   Conclusions

In this paper, we are concerned with a class of reaction-diffusion epidemic models whose underlying ODE systems are autonomous. We have proposed a computational approach to efficiently calculate and analyze the basic reproduction numbers, $\mathcal{R}_0^{{\text{PDE}}}$, for such PDE models. The present work is an extension of our previous study in ^[24] from one-dimensional spatial domains to $k$-dimensional domains for any positive integer $k$. This extension contributes to a more holistic understanding for the basic reproduction numbers of reaction-diffusion epidemic systems, and allows broader applications of the methodology in epidemiological studies.

Our numerical method transfers the computation of $\mathcal{R}_0^{{\text{\text{PDE}}}}$, defined as the spectral radius of an operator that is infinite-dimensional, to the calculation of the principal eigenvalue associated with a finite-dimensional matrix. Such a formulation enables us to apply the matrix theory to analyze and compare the basic reproduction numbers for the PDE system and its underlying ODE system. We have found that $\mathcal{R}_0^{{\text{PDE}}} = \mathcal{R}_0^{{\text{ODE}}}$ in several special but important cases, such as (i) a single infected compartment in the system; (ii) constant diffusion rates; (iii) uniform diffusion patterns in the infected compartments; and (iv) the presence of partial diffusion. For these scenarios, the computation of $\mathcal{R}_0^{{\text{PDE}}}$ can be conveniently handled by using $\mathcal{R}_0^{{\text{ODE}}}$, saving unnecessary efforts in the operator and eigenvalue analysis associated with the PDE models.

Prior studies on the basic reproduction numbers of reaction-diffusion epidemic models are often focused on the analysis of the asymptotic profiles when the constant diffusion rates tend to zero or infinity (e.g., ^{[13,14,15,16]}). In contrast, our work aims to explore a more general relationship between the basic reproduction numbers of the PDE models with variable diffusion rates and those of their underlying autonomous ODE models. Another unique feature of our study is that the analytical work is inspired by the numerical formulation and involves only elementary numerical techniques and matrix theory. The findings in this study help to efficiently quantify the risk of disease transmission for a large class of reaction-diffusion models. We hope to extend the methodology to reaction-convection-diffusion systems and possibly other PDE epidemic models in our future research.

Use of AI tools declaration

The authors declare that they have not used any Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was partially supported by the National Institutes of Health under grant number 1R15GM131315. We thank the two anonymous reviewers for their helpful comments that have improved the quality of the original manuscript.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix

A.   Underlying ODE model

If system (2.1) is homogeneous in space; i.e., $U = (u_1(t), ..., u_n(t))^T$, then the PDE model (2.1) is reduced to the following ODE model

where $\mathcal{F}(U) = (\mathcal{F}_1(U), ..., \mathcal{F}_n(U))^T$ and $\mathcal{V}(U) = \mathcal{V}^-(U) - \mathcal{V}^+(U) = (\mathcal{V}_1(U), ..., \mathcal{V}_n(U))^T$. Based on the framework in ^[10] for ODE epidemic models, we make the following standard assumptions:

(A1) $\mathcal{F}_i(U), \, \mathcal{V}_i^+(U), \, \mathcal{V}_i^-(U)$ are non-negative and continuously differentiable, $1\le i\le n$.

(A2) If $u_i = 0$, then $\mathcal{V}_i^- = 0, \; 1\le i\le m$.

(A3) $\mathcal{F}_i = 0$ for $i > m$.

(A4) If $U\in U_s$, then $\mathcal{F}_i = \mathcal{V}_i^+ = 0, \; i = 1, ..., m$.

Suppose that $U^0 = (0, ..., 0, u_{m+1}^0, ..., u_n^0)$ is a disease-free steady state of model (A.1) which is spatially independent. Then the basic reproduction number for the ODE system (A.1) is defined as the spectral radius of the next-generation matrix ^[10]; i.e.,

where $F$ and $V$ are $m\times m$ constant matrices with $(i, j)$ entry $F_{ij} = \frac{\partial\mathcal{F}_i}{\partial u_j}(U^0)$ and $V_{ij} = \frac{\partial\mathcal{V}_i}{\partial u_j}(U^0)$, respectively.

B.   Next-generation operator

For $t > 0$ and any solution $\phi(t, x)$ of Eq (2.4),

Hence, $\Gamma$ is the generator of the $C_0$-semigroup $T(t)$ on $C([0, 1]^k, \mathbb{R}^m)$. Note that $T(t)$ is a positive semigroup since $T(t)C([0, 1]^k, \mathbb{R}_+^m)\subset C([0, 1]^k, \mathbb{R}_+^m)$ for all $t\ge0$. Let $\sigma(\Gamma)$ denote the spectrum of the operator $\Gamma$. It then follows from Theorem 3.12 in ^[11] that

where $s(\Gamma) = \sup\{{\text{Re}}\lambda : \lambda \in \sigma(\Gamma)\}$ is the spectral bound of $\Gamma$. Due to the loss of infected individuals from natural and disease-induced mortalities, the internal evolution of individuals in the infectious compartments is usually dissipative and exponentially decaying. We thus assume

(B1) $-V$ is cooperative and $s(\Gamma) < 0$.

Using assumption (B1) and setting $\lambda = 0$ in Eq (B.2), we obtain

or $L = -F\Gamma^{-1}$. As done in ^[24], we make two additional assumptions regarding the PDE system (2.1) and its underlying ODE system (A.1):

(B2) There exists a constant $d_0$ such that $d_i(x)\ge d_0 > 0$ for any $x\in[0, 1]^k$ and $1\le i\le m$.

(B3) $V = {\text{diag}}(v_1, ..., v_m)$ with $v_i > 0, \; i = 1, ..., m$.

Condition (B2) sets a minimal diffusion rate at all spatial locations for the PDE model. Condition (B3) is satisfied by many common epidemic models. In case $V$ is not a diagonal matrix, it is often possible to introduce a new infection vector, $\widetilde{\mathcal{F}}(U)$, and a new transfer vector, $\widetilde{\mathcal{V}}(U)$, in the ODE system (A.1) such that $\frac{dU}{dt} = \mathcal{F} - \mathcal{V} = \widetilde{\mathcal{F}} - \widetilde{\mathcal{V}}$, where the vectors $\widetilde{\mathcal{F}}$ and $\widetilde{\mathcal{V}}$ produce a new generation matrix $\widetilde{F}$ and a diagonal transitive matrix $\widetilde{V}$, respectively. From ^[10], we know that $\rho(FV^{-1})$ and $\rho \big(\widetilde{F} \widetilde{V}^{-1} \big)$ are equivalent in characterizing the disease threshold in the sense that they are simultaneously higher than (or lower than) unity.