This paper focuses on numerical approximation of the basic reproduction number R0, which is the threshold defined by the spectral radius of the next-generation operator in epidemiology. Generally speaking, R0 cannot be explicitly calculated for most age-structured epidemic systems. In this paper, for a deterministic age-structured epidemic system and its stochastic version, we discretize a linear operator produced by the infective population with a theta scheme in a finite horizon, which transforms the abstract problem into the problem of solving the positive dominant eigenvalue of the next-generation matrix. This leads to a corresponding threshold R0,n. Using the spectral approximation theory, we obtain that R0,n→R0 as n→+∞. Some numerical simulations are provided to certify the theoretical results.
1.
Introduction
Many classical SIS (Susceptible-Infective-Susceptible) and SIRS (Susceptible-Infective-Recovered-Susceptible) models have been developed to study disease outbreaks [1,2,3,4,5]. Since certain diseases (e.g., childhood diseases) are age dependent, age-structured epidemic models have attracted the attention of many scholars [6,7,8,9,10,11]. In [6], Busenberg found that a sharp threshold (defined by the spectral radius of a positive linear operator) exists and can determine the global behavior of an age-structured epidemic model. In [7], Cao investigated the existence and global stability of all equilibria for an age-structured epidemic model with imperfect vaccination and relapse. It was found that, if the threshold is less than 1, the disease-free equilibrium is globally and asymptotically stable; if the threshold is greater than 1, the endemic equilibrium is globally stable. In reference [9], by discretizing the multigroup model, the authors transformed a PDE (Partial Differential Equations) system into an ODE (Ordinary Differential Equations) system, and proved that the global asymptotic stability of each equilibrium of the discretized system is completely determined by threshold R0. The threshold is defined as the basic reproduction number, which denotes the expected value of secondary cases produced by infective individuals during the entire infectious period when the entire population are susceptible [12].
As the threshold that controls disease outbreaks, R0 plays an extremely important role in assessing disease transmission trend and in reducing disease burden. However, for most age-structured epidemic equations such as the system in [13], the basic reproduction number, R0=∫a†0k(σ)e−∫σ0μ(η)dη1γ(1−e−γσ)dσ∫R˜P(ω)dω, is merely a theoretical expression of the next generation operator, and is always difficult to calculate. It is a common practice to use numerical approaches to approximate the threshold value [10,14]. Since many widely used epidemic models do not satisfy the global Lipschitz coefficients required for using the explicit Euler-Maruyama (EM) scheme, we propose the semi-implicit theta-scheme [15,16], which is known as the backward EM when θ=1, to approximate the exact basic reproduction number. We also estimate the approximate error of the exact basic reproduction number and the numerical threshold.
The novelty of this paper is that we use the theta scheme to discrete the linear operator produced by the infective population in a finite dimensional horizon, so that we can find out the spectral radius, which is the positive dominant eigenvalue of a nonnegative irreducible matrix defined by the next-generation operator. Subsequently, based on the spectral approximation theory [17], we obtain the threshold that converges to the exact basic reproduction number under a relatively weak condition (i.e., the compactness of the next-generation operator needs to be satisfied). These results are expected to be useful for studying infectious diseases.
The rest of this paper is organized as follows: in Section 2, the theta scheme is constructed based on the operator theory, and the scheme yields the numerical approximation of the basic reproduction number for a deterministic and a stochastic age-structure epidemic system. Section 3 presented several numerical simulations to illustrate the theoretical results. Concluding remarks are given in Section 4.
2.
Numerical approximation for the basic reproduction number
2.1. Theta scheme approximation for the deterministic age-structured SIRS system
In this section, we first present the age-structured SIRS epidemic model developed by [11],
where S(t,a), I(t,a) and R(t,a) denote the density of susceptible, infective and recovered individuals of age a at time t, respectively. Define the force of infectious λ(a,t) by
The condition S(t,0)=Λ means that the newborns are all susceptible, Λ is the recruitment rate of the population. S0(a), I0(a) and R0(a)∈L1(0,A) for ∀a∈[0,A]. All parameters are positive and their meanings are shown in Table 1.
Let us consider system (2.1) on the Banach space X:=L1(0,A)×L1(0,A)×L1(0,A). Let T be a linear operator defined by
φ(a)=(φ1(a),φ2(a),φ3(a))⊤∈D(T), where the domain D(T) is given as
The disease-free equilibrium of model (2.1) is E=(E0(a),0,Er(a)), where Er(a)=e−∫a0(μ(η)+γ(η))dη, and E0(a)=γ(a)Er(a)∫a0e−∫aϱμ(η)dηdϱ is the density of the susceptible population at age a in the disease-free state. Then we define a nonlinear operator F:X→X by
Let u(t)=(S(t,⋅),I(t,⋅).R(t,⋅))⊤, together with (2.2) and (2.3), system (2.1) has been rewritten as the following abstract Cauchy problem
Next, we mainly consider the second equation of (2.1). By simple calculation, the positive inverse (−T2)−1 is defined as follows
Then, according to [11], we can give the next generation operator K by
Based on the definition in [12], the basic reproduction number R0 is defined as r(K), where r(K) is the spectral radius of the operator K.
Since the form of r(K) is abstract, we can not calculate R0 explicitly. To avoid misunderstanding, we let B=T2,G=F2, φ2=ℏ∈D(B),
Hence, we discretize the following system
into a system of ordinary differential equations in Yn:=Rn, n∈N. Let Δa=A/n, ak:=kΔa,βkj:=β(ak,aj),μk:=μ(ak),νk:=ν(ak) and δk:=δ(ak),k=0,1,…,n,j=1,2,…,n. Then the abstract Cauchy system (2.5) is discretized as
where I(t) and I0 are n−column vectors, Bn and Gn are n−square matrices with the following form
where Mi=μi+νi+δi(i=1,⋯,n), Ni=(1−θ)E0i+θE0i+1(i=0,⋯,n−1). The additional parameter θ∈[0,1] allows us to control the implicitness of the numerical scheme [16], for technical reasons we always require θ≥12. Here we denote the next generation matrix Kn:=Gn(−Bn)−1, R0,n:=r(Kn) is the threshold corresponding to R0, and R0,n can be analyzed in a finite horizon. Since −Bn is a nonsingular M-matrix, and (−Bn)−1 is positive. Hence, from the Perron-Frobenius theorem [18], we know that r(Kn) is the positive dominant eigenvalue with algebraic multiplicity 1.
We give two bounded linear operators P:Y→Yn and J:Yn→Y as follows
where k is the kth entry of a vector, ⊤ is the transpose of matrix ψ, and χ(ak,ak+1](a) is the indicator function which implies that
From Section 4.1 in [19], we know that for all n∈N, ‖Pn‖≤1 and ‖Jn‖≤1. We denote ‖⋅‖Yn is the norm in Yn, and
Next, we apply the spectral approximation theory to present the convergence theorem of the basic reproduction number.
Theorem 2.1. Assuming that K is compact, if for any ℏ∈Y, limn→+∞‖JnKnPnℏ−Kℏ‖Y=0, then R0,n→R0 as n→+∞, preserving algebraic multiplicity 1.
Proof. It is easy to see K is strictly positive and irreducible, then from Theorem 3 in [20] and the Krein-Rutman theorem in [21], yield that R0=r(K)>0 is the maximum eigenvalue of operator K. By a simple calculation, the inverse matrix of −Bn is shown as follows
then we have
where ˉE0 and ˉβ denote the upper bounds of E0 and β, respectively. μ_, ν_ and δ_ denote the lower bounds of μ, ν and δ, respectively. They are both finite positive.
In addition, we give the following assumption to make that K is compact.
Assumption 2.1. For any h>0,
where E0β is extended by E0(a)β(a,ϱ)=0 for any a,ϱ∈(−∞,0)∪(A,∞).
The above assumption implies that the operator K keep the compactness [11, Assumption 4.4]. In order to prove JnKnPn converges to K point by point, we provide the following lemma.
Lemma 2.1. For all ℏ∈Y, limn→+∞‖JnKnPnℏ−Kℏ‖Y=0.
Proof. For any ℏ∈Y, we obtain
Since ‖Jn‖≤1, and for any n∈N, ‖Gn‖≤A¯E0ˉβ, we have L=‖Jn‖‖Gn‖=A¯E0ˉβ. Next we estimate the first term in the right-hand of (2.11), then
where ϕ:=(−B)−1ℏ∈D(B), and for any ψ=(ψ1,ψ2,⋯,ψn)⊤∈Yn,
namely, ‖(−Bn)−1‖≤A. From (2.7), we obtain
where a0=a−1=0. By the mean value theorem, we have
where ω(f,r) denotes the modulus of continuity. We know that ω(f,r) is defined by sup|x−y|≤r|f(x)−f(y)| with the following property
Hence, ‖(−Bn)−1Pnℏ−(−B)−1Pnℏ‖Yn→0 holds. Then we consider the second term of (2.11) as follows
where ω(E0,Δa)→0(Δa→0) and ω(β,Δa)→0(Δa→0), respectively. Hence,
Combine the above discussion, we have limn→+∞‖JnKnPnℏ−Kℏ‖Y=0.
By virtue of Assumption 2.1 and Lemma 2.1, we know that Theorem 2.1 holds. Namely, R0,n→R0 as n→+∞, preserving algebraic multiplicity 1.
2.2. Theta scheme approximation for the stochastic age-structured SIRS system
In this section, we seem the natural mortality μ(a) as a random variable μ(a)−σdBtdt, where Bt is a standard Brownian motion, σ is the intensity of noise perturbation. Then, replace μ(a) with μ(a)−σdBtdt in system (2.1), we can obtain a stochastic age-structured SIRS model
Next, we analysis the stochastic basic reproduction number. In the same way, we take the infective population of system (2.13) into account, and substitute S(t,a)=E0(a) into it, we derive
According to the general definition of the stochastic basic reproduction number, the following two operators are defined on Y:=L1(0,A)
and
Using A and F to rewrite (2.14) as
Then we have
The next generation operator T is shown by
Samely, we define r(T) as the basic reproduction number Rs0 of the stochastic system (2.13), and Rs0,n:=r(T) is the threshold corresponding to Rs0.
Next, we discretize (2.16) in Yn:=Rn, n∈N. Then the system (2.16) is discretized into the following equation
where An is defined as the same as Bn(An:=Bn), and
where θ∈[12,1], Qi=(1−σ22)[(1−θ)E0i+θE0i+1](i=0,1,⋯,n−1).
Theorem 2.2. From Theorem 2.1, we know that T is irreducible, compact and strictly positive. If
for any ℏ∈Y, then
where Jn and Pn are defined as (2.7).
Proof. Obviously, Rs0=r(T)>0, and r(T) is the spectral radius of operator T. We know that (−An)−1=(−Bn)−1, and (−Bn)−1 is given by (2.9). Then we have
where E_0 is the lower bound of E0.
Next, we verify that limΔ→0‖JnTnPnℏ−Tℏ‖Y=0. For any ℏ∈Y, we have
where the first term of (2.18)
is similar to the first term in the right-hand of (2.11), so it is easy to see that
Next we estimated the second term of (2.18). Let ϖ:=(−A)−1ℏ∈D(A), we obtain
Thus, ‖JnFnPn(−A)−1ℏ−F(−A)−1ℏ‖Y→0 holds. Hence, we obtain the desired assertion.
In conclusion, Theorem 2.2 holds, which implies that Rs0,n→Rs0 as Δ→0, preserving algebraic multiplicity 1.
Remark 2.1. Compared with [10], our paper has two advantages:
● [10] employed a backward Euler method to approximate R0, and obtain the numerical threshold Rn0→R0 as n→∞. In present paper, we propose a θ method is know as the backward EM when θ=1, and the explicit Euler-Maruyama(EM) scheme when θ=0. The θ scheme has the parameter θ, and different θ values give different convergence rates. Therefore, we can use the θ method to find the optimal convergence rate. And the backward Euler method is a special case when θ=1 of our method. Our work provides an extension of [10].
● A deterministic age-structured epidemic model is discussed in [10], but in present paper, we studied not only the deterministic system but also the stochastic age-structured epidemic model, and the stochastic system is more practical.
3.
Numerical simulations
In this section, numerical examples are shown to verify our Theorems. In what follows, let A=100, μ(a)=0.2(1+a3103) ([13], see Fig. 1(a)), γ(a)=γ=0.25, ν(a)=ν=0.1 and δ(a)=δ=0.05 (see [22]). Thus, E0(a)=γ(a)Er(a)∫a0e−∫aϱμ(η)dηdϱ=0.25e(−0.45a−a42×104)∫a0eϱ(ϱ32×104+0.2)−a(a32×104+0.2)dϱ. Based on numerical integration for E0(a), we obtain Fig. 1(b).
In this example, we do not specify what kinds of influenza-like disease it is, and the value of R0 is in the range of 2-3 [23]. We assumption that the disease is more likely to transmission between individuals with similar ages [10], then we let β(a,ϱ)=kJ(a−ϱ), where k=200 and J(x)=0.6(−x2+1002)×10−6+0.001 is a normalized distance function. Thus, we can easily verify that Assumption 2.1 is true. Hence, Theorem 2.1 and 2.2 hold, which implies that R0,n→R0(Rs0,n→Rs0) as n→+∞.
3.1. Numerical approximation of R0,n for the deterministic system
Let θ=0.5, and choose R0,1000≈2.57673470573749=:R∗ as a reference value for R0. From Fig. 2(a), we see that the threshold R0,n for the discretized system (2.6) respect to the reference value R∗ as n increases. Further more, the error R∗−R0,n converges to zero as n increases (see Fig. 2(b)). In Fig. 3, we show the numerical simulations of R0,n at θ=0.5, θ=0.7 and θ=0.9, respectively. It is obvious to see that the value of θ has a certain impact on the convergence rate of R0,n. The bigger value of θ, the faster rate of convergence. This implies that the backward EM method would make the convergence faster. Our paper verified the work of [10].
3.2. Numerical approximation of Rs0,n for the stochastic system
In this example, let σ=0.1, and we also choose Rs0,1000≈2.56385103220880=:R∗s as a reference for Rs0. Similarly, the threshold value Rs0,n for the discretized system (2.17) respect to the reference value R∗s (see Fig. 4(a)) and the error R∗s−Rs0,n converges to zero as n increases (see Fig. 4(b)). Fig. 5(a) give a comparation for Rs0,n at σ=0.1, σ=0.5 and σ=0.8, respectively. We can see that the intensity of environmental disturbance has a great influence on the threshold Rs0,n. The higher value of σ, the smaller value of Rs0,n. This means that the intensity of environmental fluctuation can reduce the threshold of disease outbreak, which may be a better measure to control disease outbreak. We show a 3D simulation of Rs0,n corresponding to θ∈[0.5,1] and σ∈[0,1] in Fig. 5(b), the effect of σ on the threshold Rs0,n with the change of θ is further explained.
4.
Concluding remarks
For the age-structure epidemic model, the basic reproduction number is defined as an integral and difficult to be estimated. Hence, it is necessary to approximate it using numerical methods. This paper investigates the numerical approximation of two basic reproduction numbers for deterministic and stochastic age-structured epidemic systems, respectively. We use the theta scheme to discrete the infective population in a finite space, so that the two abstract basic reproduction numbers can be calculated explicitly. Afterward, using the spectral approximation theory, we obtain the numerical threshold that converges to the exact value as n increases. We also estimate the approximation error between the exact basic reproduction number and its numerical approximation. Finally, several numerical simulations are shown to illustrate our theoretical results. The numerical results show that, for the deterministic system, the convergence rate of R0,n is faster when θ is bigger under the condition of θ∈[12,1]. For θ∈[0,12], the proof of the pointwise convergence in Lemma 2.2 remains challenging, and is warranted to be investigated in a future study. For the stochastic system, the appropriate noise intensity can reduce the threshold of disease outbreak.
Acknowledgments
The research is supported by the Natural Science Foundation of China (Grant number 11661064).
Conflict of interest
The authors declare there is no conflict of interest.