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

We consider a class of k-dimensional reaction-diffusion epidemic models (k=1,2,⋯) that are developed from autonomous ODE systems. We present a computational approach for the calculation and analysis of their basic reproduction numbers. Particularly, we apply matrix theory to study the relationship between the basic reproduction numbers of the PDE models and those of their underlying ODE models. We show that the basic reproduction numbers are the same for these PDE models and their associated ODE models in several important scenarios. We additionally provide two numerical examples to verify our analytical results.


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 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, ℛ 0 PDE , and that of its underlying ODE model, ℛ 0 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 ℛ 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 ≥ 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 ℛ 0 PDE and the ODE-based ℛ 0 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 ℛ 0 PDE and then analyze the relationship between ℛ 0 PDE and ℛ 0

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.

Reaction-diffusion epidemic model
Let n be a positive integer and U t, x = u 1 t, x , …, u n t, x 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 ≤ i ≤ n at time t and location x.We are concerned with the k-dimensional spatial domain 0, 1 k , where k ≥ 1 is an integer.We consider the following reaction-diffusion epidemic system with appropriate initial conditions.In this model, d i x 1 ≤ i ≤ n denotes the diffusion rate at location x and is assumed to be continuously differentiable on 0, 1 k .ℱ i U denotes the rate of generation for newly infected individuals in compartment i, and with V i + denoting the transfer rate of individuals into compartment i and V i − the transfer rate of individuals out of compartment i.Note that ℱ i U and V i U , i = 1, 2, ⋯, 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, v is the unit normal vector on the boundary 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 ≤ m < n.Consequently, the set of all disease-free steady states is defined as 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 , ℝ m associated with the following linear reaction-diffusion equation: (2.4) 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 F T t U m x and the distribution of the total new infections is represented by ∫ 0 +∞ F T 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 Γ : ℝ m ℝ m by where, for T , with c j x = c j1 x , c j2 x , …, c jk x , 1 ≤ j ≤ 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 ℛ 0 PDE , based on Eqs (2.6) and (2.8), and the analysis of its relationship with the basic reproduction number ℛ 0 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.

Numerical formulation and matrix analysis
Let λ be an eigenvalue of the operator L such that L ϕ x = λϕ x for an eigenvector T .Then from Eq (2.8) we have . Based on the condition (B3), this equation can be written as Pick a sufficiently large integer N > 0 and denote for any integers 0 ≤ j 1 , j 2 , …, j k ≤ N and 1 ≤ r ≤ k.Apply the standard centered difference scheme to Eq (3.2) on the spatial domain [0, 1] k .Then for any 0 ≤ j 1 , …, j k ≤ N, we obtain and ψ p , 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 ≤ r ≤ 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 = a ij p is a (N + 1) k × (N + 1) k matrix, 1 ≤ i, j ≤ (N + 1) k , and Note that for any 0 ≤ j 1 , …, j k ≤ N, the coefficient of ψ p j 1 …j k in Eq (3. 3) is a diagonal entry of A p , which is equal to a positive number 2kN 2 d 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 Nc pr j 1 …j k are nonpositive for all 1 ≤ r ≤ k.Hence for any eigenvalue λ of A p , by the Gershgorin Circle Theorem, there Thus, Re λ ≥ v p > 0 and A p is invertible.Moreover, we have This leads to the following lemma.
Lemma 3.1.Let N > N 0 and λ A p be an eigenvalue of matrix A p .Then, for 1 ≤ p ≤ m, the real part of λ A p satisfies Re λ A p ≥ v p , and, consequently, A p is invertible.
In addition, for any 0 ≤ j 1 , …, j k ≤ N, if we fix ψ p j 1 … j r + 1 …j k = ψ p j 1 …j k = ψ p j 1 … j r − 1 …j k = 1 for all 1 ≤ r ≤ k in Eq (3.3), then the left-hand side of Eq (3. 3) is the sum of the 1 + ∑ r = 1 k j r (N + 1) k − r -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 ρ A p −1 = 1/v p .We obtain the following result.
Lemma 3.2.For all N > N 0 , we have ρ A p T T , and Then A is invertible and Ψ ≈ A −1 Φ 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 × (N + 1) k identity matrix and ⊗ denotes the Kronecker product that is defined as follows: for any r × s matrix M = m ij and p × q matrix Q, With the substitution of Ψ ≈ A −1 Φ 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 Hence, for any ε > 0, we can pick N sufficiently large such that Letting ε 0, we obtain our central result for the computation of ℛ 0 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 ℛ 0 PDE .
In what follows, we apply matrix theory to conduct an analysis of ℛ 0 PDE and its connection to , based on Eq (3.8).We first introduce the following lemma.
Lemma 3.3.Assume that X = x ij is an m × m matrix and Y ij 1 ≤ i, j ≤ m are n × 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 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 ℛ 0 PDE and ℛ 0 ODE in the following three theorems.Proof.( 1) Let e = (1, 1…, 1) T be a vector with all the (N + 1) k − 1 entries being 1 's, and let . Since the sum of each row of A i is v i , we have where α i is a (N + 1) k − 1 -dimensional vector and S i is a (N + 1) k − 1 × (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 (2) This statement directly follows from Lemma 3.2, since

Theorem 3.2. If the matrix set
2) is a commuting family where each pair of matrices commute with each other, then ℛ 0 Proof.By Theorem 3.1 (1), it suffices to show that ℛ 0 is a commuting family.Hence, there exists a nonsingular matrix Q such that where B i is an upper triangular matrix with diagonal elements α i1 , …, α i(N + 1) k and α ij ≤ 1/v i for i = 1, …, m; j = 1, …, (N + 1) k .Consequently, by Lemma 3.3, we have where O j = F 11 α 1j ⋯ F 1m α mj ⋮ ⋮ ⋮ F m1 α 1j ⋯ F mm α mj for 1 ≤ j ≤ (N + 1) k .Thus, Eq (3.9) indicates that Since F is nonnegative by assumptions (A1) and (A4), then O j ≤ F V −1 , where Next, we provide sufficient and necessary conditions to characterize the scenarios where 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 δ p , σ pr , and continuous functions d x , g r x , 1 ≤ p ≤ m, 1 ≤ r ≤ k, such that d p x = δ p d x , c pr x = σ pr g r x and δ p σ qr = δ q σ pr , σ pi σ qj = σ qi σ pj for any 1 ≤ p, q ≤ m and 1 ≤ r, i, j ≤ k.
Proof.For any 0 ≤ j 1 , …, j k ≤ N, we denote the 1 + ∑ r = 1 k j r (N + 1) k − r -th entry of a and G p satisfies G p Ψ p j 1 …j k = ∑ r = 1 k c pr j 1 …j k ψ p j 1 … j r + 1 …j k + ψ p j 1 … j r − 1 …j k .
Thus, for any 1 ≤ p, q ≤ m, the equality Since Eq (3.11) holds for any N > N 0 , we can conclude that for any 1 ≤ p, q ≤ m.Following similar algebraic manipulations as those in [24], we can conclude the following: (i) H p H q = H q H p implies that there exist constants δ p , 1 ≤ p ≤ m, and a continuous function d x such that (ii) G p G q = G q G p implies that there exist constants σ pr , 1 ≤ p ≤ m, 1 ≤ r ≤ k, and continuous functions g r x , 1 ≤ r ≤ k, such that c pr x = σ pr g r x , σ pi σ qj = σ qi σ pj for any 1 ≤ p, q ≤ m, 1 ≤ r, i, j ≤ k.Substitute functions d p x = δ p d x , c pr x = σ pr g r x into matrices H p , H q , G p and G q .Then for any 1 ≤ p, q ≤ m, we can obtain that H p G q + G p H q = H q G p + G q H p holds if and only if δ p σ qr = δ q σ pr for all 1 ≤ r ≤ 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

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.

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 ∈ [0, 1] 2 .We study the following 2D reaction-diffusion SIR system, which is extended from the model presented in [4]: The constant parameters Λ, α, μ, and γ 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 ∂[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 = αΛ μ and V = μ + γ.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 ℛ 0 PDE and ℛ 0 ODE for this model.
PDE is computed by our numerical method based on Eq (3.8).The values of 1, where A 1 is the matrix obtained in Eq (3.4) from the single infectious compartment I.Meanwhile, since 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 ℛ 0 PDE based on ρ F ⊗ I (N + 1) 2 A −1 almost perfectly match ℛ 0 ODE , and this pattern continues for all N ≥ 40.
Next, we verify that ℛ 0 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 ℛ 0 PDE < 1 and ℛ 0 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 ℛ 0 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 ℛ 0 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.

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 ∈ [0, 1] 3 , with Neumann boundary conditions and appropriate initial conditions.The parameters α and β represent, respectively, the direct and indirect transmission rates, ξ 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 τ 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.From Corollary 3.4, we obtain To provide numerical evidence for this analytical relationship, we compare ℛ 0 PDE and ℛ 0 ODE in Figure 3, where we plot ρ F ⊗ I (N + 1) 3 A −1 versus N N = 1, 2, ⋯ for the model (4.2).Note 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 PDE , we set the pathogen diffusion rate as d B x = 10 −3 sin 100 x 1 + x 2 + x 3 + 1.01   in this test.We observe a pattern similar to that in Figure 1.Specifically, the numerical approximations of ℛ 0 PDE based on ρ F ⊗ I (N + 1) 3 A −1 coincide with ℛ 0 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 ℛ 0 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 ℛ 0 PDE < 1 and ℛ 0 PDE > 1, and set the initial concentration of the pathogen as B 0, x 1 , x 2 , x 3 = 10 5 3 − x 1 − 0.5 2 − x 2 − 0.5 2 − x 3 − 0.5 2 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.

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, ℛ 0 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 ℛ 0 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 ℛ 0 PDE = ℛ 0 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 ℛ 0 PDE can be conveniently handled by using ℛ 0 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 reactiondiffusion models.We hope to extend the methodology to reaction-convectiondiffusion systems and possibly other PDE epidemic models in our future research.Two typical scenarios of B t, x 1 , 0.5, 0.5 versus t and

(N + 1
) k − dimensional vector β by (β) j 1 …j k and write A p = N 2 H p + N 2 G p + v p I (N + 1) k for p = 1, …, m, where H p satisfies H p Ψ p j 1 …j k = − d

ℛ 0 PDECorollary 3 . 1 .Corollary 3 . 2 .PCorollary 3 . 3 .PCorollary 3 . 4 .
for the reaction-diffusion system (2.1) and the basic reproduction number ℛ 0 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).If there exist constants δ p , 1 ≤ p ≤ m, and a continuous function d x such that d p x = δ p d x for all p = 1, …, m in system (2.1), then ℛ 0 P DE = ℛ 0 ODE .(A scenario with constant diffusion rates.)If the diffusion rates of all the infected compartments are positive constants in system (2.1), then ℛ 0 (A scenario with uniform diffusion patterns.)If d p x = d q x for all the infected compartments 1 ≤ p, q ≤ m in system (2.1), then ℛ 0 (A scenario with partial diffusion.)If d p x = 0 for p = 1, …, m − 1 and d m x ≥ d 0 > 0 in system (2.1), then ℛ 0 P DE = ℛ 0 ODE .