
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.
Citation: Chayu Yang, Jin Wang. Computation of the basic reproduction numbers for reaction-diffusion epidemic models[J]. Mathematical Biosciences and Engineering, 2023, 20(8): 15201-15218. doi: 10.3934/mbe.2023680
[1] | Liqiong Pu, Zhigui Lin . A diffusive SIS epidemic model in a heterogeneous and periodically evolvingenvironment. Mathematical Biosciences and Engineering, 2019, 16(4): 3094-3110. doi: 10.3934/mbe.2019153 |
[2] | Kazuo Yamazaki, Xueying Wang . Global stability and uniform persistence of the reaction-convection-diffusion cholera epidemic model. Mathematical Biosciences and Engineering, 2017, 14(2): 559-579. doi: 10.3934/mbe.2017033 |
[3] | Chichia Chiu, Jui-Ling Yu . An optimal adaptive time-stepping scheme for solving reaction-diffusion-chemotaxis systems. Mathematical Biosciences and Engineering, 2007, 4(2): 187-203. doi: 10.3934/mbe.2007.4.187 |
[4] | Lin Zhao, Zhi-Cheng Wang, Liang Zhang . Threshold dynamics of a time periodic and two–group epidemic model with distributed delay. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1535-1563. doi: 10.3934/mbe.2017080 |
[5] | Aníbal Coronel, Fernando Huancas, Ian Hess, Alex Tello . The diffusion identification in a SIS reaction-diffusion system. Mathematical Biosciences and Engineering, 2024, 21(1): 562-581. doi: 10.3934/mbe.2024024 |
[6] | Nicolas Bacaër, Cheikh Sokhna . A reaction-diffusion system modeling the spread of resistance to an antimalarial drug. Mathematical Biosciences and Engineering, 2005, 2(2): 227-238. doi: 10.3934/mbe.2005.2.227 |
[7] | Azmy S. Ackleh, Keng Deng, Yixiang Wu . Competitive exclusion and coexistence in a two-strain pathogen model with diffusion. Mathematical Biosciences and Engineering, 2016, 13(1): 1-18. doi: 10.3934/mbe.2016.13.1 |
[8] | Wenzhang Huang, Maoan Han, Kaiyu Liu . Dynamics of an SIS reaction-diffusion epidemic model for disease transmission. Mathematical Biosciences and Engineering, 2010, 7(1): 51-66. doi: 10.3934/mbe.2010.7.51 |
[9] | Danfeng Pang, Yanni Xiao . The SIS model with diffusion of virus in the environment. Mathematical Biosciences and Engineering, 2019, 16(4): 2852-2874. doi: 10.3934/mbe.2019141 |
[10] | Blessing O. Emerenini, Stefanie Sonner, Hermann J. Eberl . Mathematical analysis of a quorum sensing induced biofilm dispersal model and numerical simulation of hollowing effects. Mathematical Biosciences and Engineering, 2017, 14(3): 625-653. doi: 10.3934/mbe.2017036 |
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.
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 R0, 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. R0 is often used to characterize the threshold behavior of an epidemic, with disease eradication if R0<1 and disease persistence if R0>1. The theory of the basic reproduction numbers for ODE-based autonomous epidemic models is well developed, and the calculation of R0 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 R0, 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 R0 as the spectral radius of a next-infection operator for a general class of reaction-diffusion models. They showed that the R0 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 R0 for reaction-diffusion epidemic systems with constant diffusion rates were investigated in [13,14,15,16]. Particularly, Chen and Shi [14] showed that R0 approaches the spectral radius of a spatially averaged next-generation matrix when the diffusion rates tend to infinity, and R0 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, RPDE0, and that of its underlying ODE model, RODE0. 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 R0 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 RPDE0 and the ODE-based RODE0. 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 RPDE0 and then analyze the relationship between RPDE0 and RODE0. In Section 4, we provide specific numerical examples to verify our analytical findings. Finally, we conclude the paper with some discussion in Section 5.
Let n be a positive integer and U(t,x)=(u1(t,x),...,un(t,x))T be a vector-valued function that represents the hosts and pathogens related to an infectious disease, with each ui(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
{∂ui∂t=∇⋅(di(x)∇ui)+Fi(U)−Vi(U),1≤i≤n,t>0,x∈[0,1]k;∂ui∂ν=0,1≤i≤n,t>0,x∈∂[0,1]k, | (2.1) |
with appropriate initial conditions. In this model, di(x) (1≤i≤n) denotes the diffusion rate at location x and is assumed to be continuously differentiable on [0,1]k. Fi(U) denotes the rate of generation for newly infected individuals in compartment i, and Vi(U)=V−i(U)−V+i(U), 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 Fi(U) and Vi(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, ν is the unit normal vector on the boundary ∂[0,1]k.
Without loss of generality, we assume that UI=UI(t,x)=(u1,...,um)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 Us={U≥0:ui=0,i=1,...,m}.
System (2.1) can be re-written as
{∂ui∂t=di(x)Δui−ci(x)⋅∇ui+Fi(U)−Vi(U),1≤i≤n,t>0,x∈[0,1]k;∂ui∂ν=0,1≤i≤n,t>0,x∈∂[0,1]k, | (2.2) |
where
ci(x)=(ci1(x),ci2(x),...,cik(x))=−∇di(x),1≤i≤n. | (2.3) |
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,Rm) associated with the following linear reaction-diffusion equation:
{∂ui∂t=di(x)Δui−ci(x)⋅∇ui−Vi(U),1≤i≤m,t>0,x∈[0,1]k;∂ui∂ν=0,1≤i≤m,t>0,x∈∂[0,1]k. | (2.4) |
Let the distribution of the initial infections, i.e., UI(0,x), be Um(x)=(u1(x),...,um(x))T. Then the distribution of these infections after time t>0 is given by T(t)(Um(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)(Um(x)) and the distribution of the total new infections is represented by ∫+∞0FT(t)(Um(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
L:Um(x)⟼∫+∞0FT(t)(Um(x))dt. | (2.5) |
Consequently, the basic reproduction number for the PDE system (2.2) is the spectral radius of the operator L:
RPDE0=ρ(L). | (2.6) |
Meanwhile, we introduce the operator Γ:Rm⟼Rm by
Γ=Dm(x)Δ−Cm(x)⋅∇−V, | (2.7) |
where, for U=(u1,⋯,um)T∈Rm, (Dm(x)Δ)U=(d1(x)Δu1,...,dm(x)Δum)T and (Cm(x)⋅∇)U=(c1(x)⋅∇u1,...,cm(x)⋅∇um)T, with cj(x)=(cj1(x),cj2(x),...,cjk(x)), 1≤j≤m. Then we can obtain an essential characterization of the next-infection operator,
L=−FΓ−1, | (2.8) |
with details provided in Appendix B.
Below we focus our attention on the numerical calculation of RPDE0, based on Eqs (2.6) and (2.8), and the analysis of its relationship with the basic reproduction number RODE0 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.
Let λ be an eigenvalue of the operator L such that L(ϕ(x))=λϕ(x) for an eigenvector ϕ(x)=(ϕ1(x),...,ϕm(x))T. Then from Eq (2.8) we have
−FΓ−1(ϕ(x))=λϕ(x). | (3.1) |
Let ψ(x)=−Γ−1(ϕ(x)), where ψ(x)=(ψ1(x),...,ψm(x))T. Then −Γ(ψ(x))=ϕ(x). Based on the condition (B3), this equation can be written as
−(dp(x)Δψp(x)−cp(x)⋅∇ψp(x)−vpψp(x))=ϕp(x),1≤p≤m. | (3.2) |
Pick a sufficiently large integer N>0 and denote
dj1...jkp=dp(j1N,j2N,...,jkN),cj1...jkpj=cpr(j1N,j2N,...,jkN),ψj1...jkp=ψp(j1N,j2N,...,jkN),ϕj1...jkp=ϕp(j1N,j2N,...,jkN), |
for any integers 0≤j1,j2,...,jk≤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≤j1,...,jk≤N, we obtain
−dj1...jkpk∑r=1ψj1...(jr+1)...jkp−2ψj1...jkp+ψj1...(jr−1)...jkp1/N2+k∑r=1cj1...jkprψj1...(jr+1)...jkp−ψj1...(jr−1)...jkp2/N+vpψj1...jkp≈ϕj1...jki, | (3.3) |
and ψj1...(−1)...jkp=ψj1...(1)...jkp,ψj1...(N+1)...jkp=ψj1...(N−1)...jkp, where (−1),(1),(N+1), and (N−1) are all in the (jr) position inside the permutation j1...jk 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
ApΨp≈Φp, | (3.4) |
where Ap=(a(p)ij) is a (N+1)k×(N+1)k matrix, 1≤i,j≤(N+1)k, and
Ψp=(ψ0...00p,...,ψ0...0Np,...,ψ0...N...0p,...,ψ0...N...Np,...,ψN...N0p,...,ψN...NNp)T, |
Φp=(ϕ0...00p,...,ϕ0...0Np,...,ϕ0...N...0p,...,ϕ0...N...Np,...,ϕN...N0p,...,ϕN...NNp)T. |
Note that for any 0≤j1,...,jk≤N, the coefficient of ψj1...jkp in Eq (3.3) is a diagonal entry of Ap, which is equal to a positive number 2kN2dj1...jkp+vp. Define
N0=max1≤p≤m{Np},whereNp=max{||cpr(x)||∞2d0,1≤r≤k}, |
and where d0 is a positive lower bound for the diffusion rates (see condition (B2)). Then for N>N0, the off-diagonal entries −N2dj1...jkp+12Ncj1...jkpr and −N2dj1...jkp−12Ncj1...jkpr are nonpositive for all 1≤r≤k. Hence for any eigenvalue λ of Ap, by {the Gershgorin Circle Theorem}, there exists 1≤i≤(N+1)k such that
|λ−a(p)ii|≤∑j≠i|a(p)ij|=|a(p)ii−vp|. |
Thus, Re(λ)≥vp>0 and Ap is invertible. Moreover, we have ρ(A−1p)≤1/Re(λ)≤1/vp,1≤p≤m. This leads to the following lemma.
Lemma 3.1. Let N>N0 and λAp be an eigenvalue of matrix Ap. Then, for 1≤p≤m, the real part of λAp satisfies Re(λAp)≥vp, and, consequently, Ap is invertible.
In addition, for any 0≤j1,...,jk≤N, if we fix ψj1...(jr+1)...jkp=ψj1...jkp=ψj1...(jr−1)...jkp=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+k∑r=1jr(N+1)k−r)-th row of the matrix Ap, which is obviously equal to vp. Hence the sum of each row of Ap is vp, which implies vp is an eigenvalue of Ap, and consequently, 1/vp is an eigenvalue of A−1p. Therefore ρ(A−1p)=1/vp. We obtain the following result.
Lemma 3.2. For all N>N0, we have ρ(A−1p)=1/vp,1≤p≤m.
Denote Ψ=(ΨT1,...,ΨTm)T, Φ=(ΦT1,...,ΦTm)T, and
A=diag(A1,...,Am). |
Then A is invertible and Ψ≈A−1Φ by Eq (3.4). It follows from Eq (3.1) that
Fψ(x)=−FΓ−1(ϕ(x))=λϕ(x), | (3.5) |
which yields
(F⊗I(N+1)k)Ψ=λΦ | (3.6) |
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=(mij) and p×q matrix Q,
M⊗Q=[m11Q⋯m1sQ⋮⋱⋮mr1Q⋯mrsQ]. |
With the substitution of Ψ≈A−1Φ into Eq (3.6), our numerical formulation leads to
(F⊗I(N+1)k)A−1Φ≈λΦ. | (3.7) |
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→∞. Hence, for any ε>0, we can pick N sufficiently large such that
|ρ((F⊗I(N+1)k)A−1)−ρ(L)|<ε. |
Letting ε→0, we obtain our central result for the computation of RPDE0:
RPDE0=limN→∞ρ((F⊗I(N+1)k)A−1). | (3.8) |
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 RPDE0. In what follows, we apply matrix theory to conduct an analysis of RPDE0 and its connection to RODE0, based on Eq (3.8). We first introduce the following lemma.
Lemma 3.3. Assume that X=(xij) is an m×m matrix and Yij(1≤i,j≤m) are n×n matrices. If there exists a nonsingular matrix P such that P−1YijP=Uij for all i,j=1,...,m, where Uij is an upper triangular matrix with diagonal entries y(1)ij,...,y(n)ij, then
det[x11Y11⋯x1mY1m⋮⋮⋮xm1Ym1⋯xmmYmm]=n∏k=1det[x11y(k)11⋯x1my(k)1m⋮⋮⋮xm1y(k)m1⋯xmmy(k)mm]. |
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 RPDE0 and RODE0 in the following three theorems.
Theorem 3.1. (1) In general, we have RPDE0≥RODE0.
(2) If F is a triangular matrix, then RPDE0=RODE0. Particularly, if m=1, then RPDE0=RODE0.
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=[10eI(N+1)k−1]. Since the sum of each row of Ai is vi, we have
P−1A−1iP=[1/viαTi0Si], |
where αi is a [(N+1)k−1]-dimensional vector and Si 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 det(λIm−FV−1) is a factor of det[λIm(N+1)k−(F⊗I(N+1)k)A−1]; that is, each eigenvalue of FV−1 is an eigenvalue of (F⊗I(N+1)k)A−1. Thus, ρ((F⊗I(N+1)k)A−1)≥ρ(FV−1). Letting N→∞, we obtain RPDE0≥RODE0.
(2) This statement directly follows from Lemma 3.2, since
ρ((F⊗I(N+1)k)A−1)=max1≤i≤m{ρ(FiiA−1i)}=max1≤i≤m{Fii/vi}=ρ(FV−1). |
Theorem 3.2. If the matrix set {Ai}mi=1 for system (2.2) is a commuting family where each pair of matrices commute with each other, then RPDE0=RODE0.
Proof. By Theorem 3.1(1), it suffices to show that RPDE0≤RODE0. Since {Ai}mi=1 is a commuting family of matrices, then {A−1i}mi=1 is a commuting family. Hence, there exists a nonsingular matrix Q such that
QA−1iQ−1=Bi, |
where Bi is an upper triangular matrix with diagonal elements αi1,...,αi(N+1)k and |αij|≤1/vi for i=1,...,m; j=1,...,(N+1)k. Consequently, by Lemma 3.3, we have
det(λIm(N+1)k−(F⊗I(N+1)k)A−1)=(N+1)k∏j=1det(λIm−Oj), | (3.9) |
where Oj=[F11α1j⋯F1mαmj⋮⋮⋮Fm1α1j⋯Fmmαmj] for 1≤j≤(N+1)k. Thus, Eq (3.9) indicates that
ρ((F⊗I(N+1)k)A−1)=max1≤j≤(N+1)k{ρ(Oj)}. | (3.10) |
Since F is nonnegative by assumptions (A1) and (A4), then |Oj|≤FV−1, where |Oj|=[|F11α1j|⋯|F1mαmj|⋮⋮⋮|Fm1α1j|⋯|Fmmαmj|]. We thus obtain ρ(Oj)≤ρ(|Oj|)≤ρ(FV−1). Therefore, ρ((F⊗I(N+1)k)A−1)≤ρ(FV−1), which implies RPDE0≤RODE0.
Next, we provide sufficient and necessary conditions to characterize the scenarios where {Ap}mp=1 is a commuting family.
Theorem 3.3. For any integer N>N0, the matrix set {Ap}mp=1 associated with system (2.2) is a commuting family if and only if there exist constants δp,σpr, and continuous functions d(x),gr(x),1≤p≤m,1≤r≤k, such that
dp(x)=δpd(x),cpr(x)=σprgr(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≤j1,...,jk≤N, we denote the (1+∑kr=1jr(N+1)k−r)-th entry of a (N+1)k-dimensional vector β by (β)j1...jk and write Ap=N2Hp+N2Gp+vpI(N+1)k for p=1,...,m, where Hp satisfies
(HpΨp)j1...jk=−dj1...jkpk∑r=1(ψj1...(jr+1)...jkp−2ψj1...jkp+ψj1...(jr−1)...jkp) |
and Gp satisfies
(GpΨp)j1...jk=k∑r=1cj1...jkpr(ψj1...(jr+1)...jkp+ψj1...(jr−1)...jkp). |
Thus, for any 1≤p,q≤m, the equality
ApAq=AqAp |
is equivalent to
N2(HpHq−HqHp)+N2(HpGq+GpHq−HqGp−GqHp)+14(GpGq−GqGp)=0. | (3.11) |
Since Eq (3.11) holds for any N>N0, we can conclude that
HpHq=HqHp,HpGq+GpHq=HqGp+GqHp,GpGq=GqGp, |
for any 1≤p,q≤m. Following similar algebraic manipulations as those in [24], we can conclude the following: (i) HpHq=HqHp implies that there exist constants δp, 1≤p≤m, and a continuous function d(x) such that
dp(x)=δpd(x); |
(ii) GpGq=GqGp implies that there exist constants σpr,1≤p≤m,1≤r≤k, and continuous functions gr(x),1≤r≤k, such that
cpr(x)=σprgr(x),σpiσqj=σqiσpj |
for any 1≤p,q≤m,1≤r,i,j≤k. Substitute functions dp(x)=δpd(x),cpr(x)=σprgr(x) into matrices Hp,Hq,Gp and Gq. Then for any 1≤p,q≤m, we can obtain that HpGq+GpHq=HqGp+GqHp 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 RPDE0 for the reaction-diffusion system (2.1) and the basic reproduction number RODE0 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 δp, 1≤p≤m, and a continuous function d(x) such that dp(x)=δpd(x) for all p=1,...,m in system (2.1), then RPDE0=RODE0.
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 RPDE0=RODE0.
Corollary 3.3. (A scenario with uniform diffusion patterns.) If dp(x)=dq(x) for all the infected compartments (1≤p,q≤m) in system (2.1), then RPDE0=RODE0.
Corollary 3.4. (A scenario with partial diffusion.) If dp(x)=0 for p=1,...,m−1 and dm(x)≥d0>0 in system (2.1), then RPDE0=RODE0.
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.
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 dS(x), dI(x) and dR(x) be their associated diffusion rates with x=(x1,x2)∈[0,1]2. We study the following 2D reaction-diffusion SIR system, which is extended from the model presented in [4]:
∂S∂t=∇⋅(dS(x)∇S)+Λ−αSI−μS,x∈[0,1]2,t>0;∂I∂t=∇⋅(dI(x)∇I)+αSI−(μ+γ)I,x∈[0,1]2,t>0;∂R∂t=∇⋅(dR(x)∇R)+γI−μR,x∈[0,1]2,t>0. | (4.1) |
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):
RPDE0=ρ(FV−1)=αΛμ(μ+γ). |
To verify this relationship, Figure 1 compares RPDE0 and RODE0 for this model. RPDE0 is computed by our numerical method based on Eq (3.8). The values of ρ((F⊗I(N+1)2)A−1)=ρ(αΛμA−11) versus N (N=1,2,⋯) are plotted in Figure 1, where A1 is the matrix obtained in Eq (3.4) from the single infectious compartment I. Meanwhile, since RODE0=ρ(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 dI(x)=sin(100(x1+x2))+2 in this test. We observe that when N is sufficiently large, the numerical values of RPDE0 based on ρ((F⊗I(N+1)2)A−1) almost perfectly match RODE0, and this pattern continues for all N≥40.
Next, we verify that RPDE0=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 RPDE0<1 and RPDE0>1 by selecting appropriate parameter values, and run the simulation for the model (4.1) with an initial infection density I(0,x1,x2)=200 in each scenario. Figure 2(a) displays the surface plot of I(t,x1,x2=0.5) versus t and x1 with RPDE0=0.88<1, where I approaches 0 over time for all x1. In contrast, Figure 2(b) displays the surface plot of I(t,x1,x2=0.5) versus t and x1 with RPDE0=3.20>1, where I increases from its initial value and remains positive for all the time. Surface plots at other fixed values of x2 and those with x1 fixed (not shown here) are qualitatively similar.
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
∂S∂t=Λ−(αI+βB)S−μS;∂I∂t=(αI+βB)S−(μ+γ)I;∂R∂t=γI−μR;∂B∂t=∇⋅(dB(x)∇B)+ξI+rB(1−BK)−τB, | (4.2) |
for t>0 and x=(x1,x2,x3)∈[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.
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 RPDE0 of system (4.2) equals RODE0 of the underlying ODE system. The infectious compartments are I and B. From the associated ODE system, we easily obtain
F=[αΛμβΛμξr]andV=[μ+γ00τ]. |
From Corollary 3.4, we obtain
RPDE0=ρ(FV−1)=12(αΛμ(μ+γ)+rτ+√(αΛμ(μ+γ)−rτ)2+4ξβΛμτ(μ+γ)). | (4.3) |
To provide numerical evidence for this analytical relationship, we compare RPDE0 and RODE0 in Figure 3, where we plot ρ((F⊗I(N+1)3)A−1) versus N (N=1,2,⋯) for the model (4.2). Note that (F⊗I(N+1)3)A−1=(αΛμ(μ+γ)I(N+1)3βΛμA−12ξμ+γI(N+1)3rA−12) and A2 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 RPDE0, we set the pathogen diffusion rate as dB(x)=10−3(sin(100(x1+x2+x3))+1.01) in this test. We observe a pattern similar to that in Figure 1. Specifically, the numerical approximations of RPDE0 based on ρ((F⊗I(N+1)3)A−1) coincide with RODE0 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 RPDE0=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 RPDE0<1 and RPDE0>1, and set the initial concentration of the pathogen as B(0,x1,x2,x3)=105(3−(x1−0.5)2−(x2−0.5)2−(x3−0.5)2) for the model (4.2). We then run the simulation in each scenario and plot the pathogen concentration B(t,x1,x2=0.5,x3=0.5) versus t and x1 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 x1, x2, and x3 have similar behaviors and are not shown here.
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, RPDE0, 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 R\text{PDE}0, 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 RPDE0=RODE0 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 RPDE0 can be conveniently handled by using RODE0, 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.
The authors declare that they have not used any Artificial Intelligence (AI) tools in the creation of this article.
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.
The authors declare that there is no conflict of interest.
If system (2.1) is homogeneous in space; i.e., U=(u1(t),...,un(t))T, then the PDE model (2.1) is reduced to the following ODE model
dUdt=F(U)−V(U),t>0, | (A.1) |
where F(U)=(F1(U),...,Fn(U))T and V(U)=V−(U)−V+(U)=(V1(U),...,Vn(U))T. Based on the framework in [10] for ODE epidemic models, we make the following standard assumptions:
(A1) Fi(U),V+i(U),V−i(U) are non-negative and continuously differentiable, 1≤i≤n.
(A2) If ui=0, then V−i=0,1≤i≤m.
(A3) Fi=0 for i>m.
(A4) If U∈Us, then Fi=V+i=0,i=1,...,m.
Suppose that U0=(0,...,0,u0m+1,...,u0n) 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.,
RODE0=ρ(FV−1), | (A.2) |
where F and V are m×m constant matrices with (i,j) entry Fij=∂Fi∂uj(U0) and Vij=∂Vi∂uj(U0), respectively.
For t>0 and any solution ϕ(t,x) of Eq (2.4),
lims→0+T(s)ϕ(t,x)−ϕ(t,x)s=lims→0+ϕ(t+s,x)−ϕ(t,x)s=∂ϕ∂t(t,x)=Γ(ϕ). | (B.1) |
Hence, Γ is the generator of the C0-semigroup T(t) on C([0,1]k,Rm). Note that T(t) is a positive semigroup since T(t)C([0,1]k,Rm+)⊂C([0,1]k,Rm+) for all t≥0. Let σ(Γ) denote the spectrum of the operator Γ. It then follows from Theorem 3.12 in [11] that
(λI−Γ)−1(ϕ)=∫+∞0eλtT(t)(ϕ)dt,∀λ>s(Γ),ϕ∈C([0,1]k,Rm), | (B.2) |
where s(Γ)=sup{Reλ:λ∈σ(Γ)} is the spectral bound of Γ. 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(Γ)<0.
Using assumption (B1) and setting λ=0 in Eq (B.2), we obtain
L(ϕ)=F∫+∞0T(t)(ϕ)dt=−FΓ−1(ϕ), | (B.3) |
or L=−FΓ−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 d0 such that di(x)≥d0>0 for any x∈[0,1]k and 1≤i≤m.
(B3) V=diag(v1,...,vm) with vi>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, ˜F(U), and a new transfer vector, ˜V(U), in the ODE system (A.1) such that dUdt=F−V=˜F−˜V, where the vectors ˜F and ˜V produce a new generation matrix ˜F and a diagonal transitive matrix ˜V, respectively. From [10], we know that ρ(FV−1) and ρ(˜F˜V−1) are equivalent in characterizing the disease threshold in the sense that they are simultaneously higher than (or lower than) unity.
[1] |
E. Bertuzzo, R. Casagrandi, M. Gatto, I. Rodriguez-Iturbe, A. Rinaldo, On spatially explicit models of cholera epidemics, J. R. Soc. Interface, 7 (2010), 321–333. https://doi.org/10.1098/rsif.2009.0204 doi: 10.1098/rsif.2009.0204
![]() |
[2] |
R. S. Cantrell, C. Cosner, The effects of spatial heterogeneity in population dynamics, J. Math. Biol., 29 (1991), 315–338. https://doi.org/10.1007/BF00167155 doi: 10.1007/BF00167155
![]() |
[3] | R. S. Cantrell, C. Cosner, Spatial Ecology via Reaction-Diffusion Equations, Wiley, 2003. https://doi.org/10.1002/0470871296 |
[4] |
K. I. Kim, Z. Lin, Q. Zhang, An SIR epidemic model with free boundary, Nonlinear Anal. Real World Appl., 14 (2013), 1992–2001. https://doi.org/10.1016/j.nonrwa.2013.02.003 doi: 10.1016/j.nonrwa.2013.02.003
![]() |
[5] |
A. Rinaldo, E. Bertuzzo, L. Mari, L. Righetto, M. Blokesch, M. Gatto, et al., Reassessment of the 2010-2011 Haiti cholera outbreak and rainfall-driven multiseason projections, PNAS, 109 (2012), 6602–6607. https://doi.org/10.1073/pnas.1203333109 doi: 10.1073/pnas.1203333109
![]() |
[6] |
F. B. Wang, J. Shi, X. Zou, Dynamics of a host-pathogen system on a bounded spatial domain, Commun. Pure Appl. Anal., 14 (2015), 2535–2560. https://doi.org/10.3934/cpaa.2015.14.2535 doi: 10.3934/cpaa.2015.14.2535
![]() |
[7] |
X. Wang, D. Gao, J. Wang, Influence of human behavior on cholera dynamics, Math. Biosci., 267 (2015), 41–52. https://doi.org/10.1016/j.mbs.2015.06.009 doi: 10.1016/j.mbs.2015.06.009
![]() |
[8] |
X. Wang, D. Posny, J. Wang, A reaction-convection-diffusion model for cholera spatial dynamics, Discrete Contin. Dyn. Syst. - Ser. B, 21 (2016), 2785–2809. https://doi.org/10.3934/dcdsb.2016073 doi: 10.3934/dcdsb.2016073
![]() |
[9] |
O. Diekmann, J. A. P. Heesterbeek, A. J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous population, J. Math. Biol., 28 (1990), 365–382. https://doi.org/10.1007/BF00178324 doi: 10.1007/BF00178324
![]() |
[10] |
P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48. https://doi.org/10.1016/S0025-5564(02)00108-6 doi: 10.1016/S0025-5564(02)00108-6
![]() |
[11] |
H. R. Thieme, Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity, SIAM J. Appl. Math., 70 (2009), 188–211. https://doi.org/10.1137/080732870 doi: 10.1137/080732870
![]() |
[12] |
W. Wang, X. Q. Zhao, Basic reproduction numbers for reaction-diffusion epidemic models, SIAM J. Appl. Dyn. Syst., 11 (2012), 1652–1673. https://doi.org/10.1137/120872942 doi: 10.1137/120872942
![]() |
[13] |
L. J. S. Allen, B. M. Bolker, Y. Lou, A. L. Nevai, Asymptotic profiles of the steady states for an SIS epidemic reaction-diffusion model, Discrete Contin. Dyn. Syst., 21 (2008), 1–20. https://doi.org/10.3934/dcds.2008.21.1 doi: 10.3934/dcds.2008.21.1
![]() |
[14] |
S. Chen, J. Shi, Asymptotic profiles of basic reproduction number for epidemic spreading in heterogeneous environment, SIAM J. Appl. Math., 80 (2020), 1247–1271. https://doi.org/10.1137/19M128907 doi: 10.1137/19M128907
![]() |
[15] |
P. Magal, G. F. Webb, Y. Wu, On the basic reproduction number of reaction-diffusion epidemic models, SIAM J. Appl. Math., 79 (2019), 284–304. https://doi.org/10.1137/18M1182243 doi: 10.1137/18M1182243
![]() |
[16] |
P. Song, Y. Lou, Y. Xiao, A spatial SEIRS reaction-diffusion model in heterogeneous environment, J. Differ. Equations, 267 (2019), 5084–5114. https://doi.org/10.1016/j.jde.2019.05.022 doi: 10.1016/j.jde.2019.05.022
![]() |
[17] |
X. Lin, Q. Wang, Asymptotic behavior of the principal eigenvalue and basic reproduction ratio for time-periodic reaction-diffusion systems with time delay, Discrete Contin. Dyn. Syst. - Ser. B, 28 (2023), 3955–3984. https://doi.org/0.3934/dcdsb.2022250 doi: 10.3934/dcdsb.2022250
![]() |
[18] |
L. Zhang, X. Q. Zhao, Asymptotic behavior of the basic reproduction ratio for periodic reaction-diffusion systems, SIAM J. Math. Anal., 53 (2021), 6873–6909. https://doi.org/10.1137/20M1366344 doi: 10.1137/20M1366344
![]() |
[19] |
C. Barril, A. Calsina, J. Ripoll, A practical approach to R0 in continuous-time ecological models, Math. Methods Appl. Sci., 41 (2018), 8432–8445. https://doi.org/10.1002/mma.4673 doi: 10.1002/mma.4673
![]() |
[20] |
D. Breda, F. Florian, J. Ripoll, R. Vermiglio, Efficient numerical computation of the basic reproduction number for structured populations, J. Comput. Appl. Math., 384 (2021), 113165. https://doi.org/10.1016/j.cam.2020.113165 doi: 10.1016/j.cam.2020.113165
![]() |
[21] |
J. Ge, C. Lei, Z. Lin, Reproduction numbers and the expanding fronts for a diffusion-advection SIS model in heterogeneous time-periodic environment, Nonlinear Anal. Real World Appl., 33 (2017), 100–120. https://doi.org/10.1016/j.nonrwa.2016.06.005 doi: 10.1016/j.nonrwa.2016.06.005
![]() |
[22] |
Y. Lou, X. Q. Zhao, A reaction-diffusion malaria model with incubation period in the vector population, J. Math. Biol., 62 (2011), 543–568. https://doi.org/10.1007/s00285-010-0346-8 doi: 10.1007/s00285-010-0346-8
![]() |
[23] |
H. Zhao, K. Wang, H. Wang, Basic reproduction ratio of a mosquito-borne disease in heterogeneous environment, J. Math. Biol., 86 (2023), 32. https://doi.org/10.1007/s00285-023-01867-y doi: 10.1007/s00285-023-01867-y
![]() |
[24] |
C. Yang, J. Wang, Basic reproduction numbers for a class of reaction-diffusion epidemic models, Bull. Math. Biol., 82 (2020), 111. https://doi.org/10.1007/s11538-020-00788-x doi: 10.1007/s11538-020-00788-x
![]() |
[25] | R. D. Richtmyer, K. W. Morton, Difference Methods for Initial-Value Problems, Second Edition, Krieger Publication Company, 1994. |
[26] | J. W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, Springer-Verlag New York, 1995. https://doi.org/10.1007/978-1-4899-7278-1 |
[27] | G. H. Golub, C. F. van Loan, Matrix Computations, Third Edition, Johns Hopkins University Press, 1996. |
[28] | Y. Saad, Numerical Methods for Large Eigenvalue Problems, Revised Edition, SIAM, 2011. |
[29] |
C. P. Gerba, Environmentally transmitted pathogens, Environ. Microbiol., 2015 (2015), 509–550. https://doi.org/10.1016/B978-0-12-394626-3.00022-3 doi: 10.1016/B978-0-12-394626-3.00022-3
![]() |
[30] |
Z. Mukandavire, S. Liao, J. Wang, H. Gaff, D. L. Smith, J. G. Morris, Estimating the reproductive numbers for the 2008-2009 cholera outbreaks in Zimbabwe, PNAS, 108 (2011), 8767–8772. https://doi.org/10.1073/pnas.1019712108 doi: 10.1073/pnas.1019712108
![]() |
[31] |
D. Posny, J. Wang, Modeling cholera in periodic environments, J. Biol. Dyn., 8 (2014), 1–19. https://doi.org/10.1080/17513758.2014.896482 doi: 10.1080/17513758.2014.896482
![]() |
[32] |
J. H. Tien, D. J. Earn, Multiple transmission pathways and disease dynamics in a waterborne pathogen model, Bull. Math. Biol., 72 (2010), 1506–1533. https://doi.org/10.1007/s11538-010-9507-6 doi: 10.1007/s11538-010-9507-6
![]() |
[33] |
T. M. Nguyen, D. Ilef, S. Jarraud, L. Rouil, C. Campese, D. Che, et al., A community-wide outbreak of legionnaires disease linked to industrial cooling towers–How far can contaminated aerosols spread? J. Infect. Dis., 193 (2006), 102–111. https://doi.org/10.1086/498575 doi: 10.1086/498575
![]() |
[34] |
T. Greenhalgh, J. L. Jimenez, K. A. Prather, Z. Tufekci, D. Fisman, R. Schooley, Ten scientific reasons in support of airborne transmission of SARS-CoV-2, Lancet, 397 (2021), 1603–1605. https://doi.org/10.1016/S0140-6736(21)00869-2 doi: 10.1016/S0140-6736(21)00869-2
![]() |
[35] |
L. D. Stetzenbach, Airborne infectious microorganisms, Encycl. Microbiol., 2009 (2009), 175–182. https://doi.org/10.1016/B978-012373944-5.00177-2 doi: 10.1016/B978-012373944-5.00177-2
![]() |
[36] |
K. Yamazaki, C. Yang, J. Wang, A partially diffusive cholera model based on a general second-order differential operator, J. Math. Anal. Appl., 501 (2021), 125181. https://doi.org/10.1016/j.jmaa.2021.125181 doi: 10.1016/j.jmaa.2021.125181
![]() |
1. | Anna M. Langmüller, Joachim Hermisson, Courtney C. Murdock, Philipp W. Messer, Catching a wave: On the suitability of traveling-wave solutions in epidemiological modeling, 2024, 00405809, 10.1016/j.tpb.2024.12.004 |