Research article Special Issues

A theta-scheme approximation of basic reproduction number for an age-structured epidemic system in a finite horizon

  • 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,nR0 as n+. Some numerical simulations are provided to certify the theoretical results.

    Citation: Wenjuan Guo, Ming Ye, Xining Li, Anke Meyer-Baese, Qimin Zhang. A theta-scheme approximation of basic reproduction number for an age-structured epidemic system in a finite horizon[J]. Mathematical Biosciences and Engineering, 2019, 16(5): 4107-4121. doi: 10.3934/mbe.2019204

    Related Papers:

    [1] Junyuan Yang, Rui Xu, Xiaofeng Luo . Dynamical analysis of an age-structured multi-group SIVS epidemic model. Mathematical Biosciences and Engineering, 2019, 16(2): 636-666. doi: 10.3934/mbe.2019031
    [2] Simone De Reggi, Francesca Scarabel, Rossana Vermiglio . Approximating reproduction numbers: a general numerical method for age-structured models. Mathematical Biosciences and Engineering, 2024, 21(4): 5360-5393. doi: 10.3934/mbe.2024236
    [3] Toshikazu Kuniya, Mimmo Iannelli . R0 and the global behavior of an age-structured SIS epidemic model with periodicity and vertical transmission. Mathematical Biosciences and Engineering, 2014, 11(4): 929-945. doi: 10.3934/mbe.2014.11.929
    [4] Zhiping Liu, Zhen Jin, Junyuan Yang, Juan Zhang . The backward bifurcation of an age-structured cholera transmission model with saturation incidence. Mathematical Biosciences and Engineering, 2022, 19(12): 12427-12447. doi: 10.3934/mbe.2022580
    [5] Yicang Zhou, Zhien Ma . Global stability of a class of discrete age-structured SIS models with immigration. Mathematical Biosciences and Engineering, 2009, 6(2): 409-425. doi: 10.3934/mbe.2009.6.409
    [6] Xi-Chao Duan, Xue-Zhi Li, Maia Martcheva . Dynamics of an age-structured heroin transmission model with vaccination and treatment. Mathematical Biosciences and Engineering, 2019, 16(1): 397-420. doi: 10.3934/mbe.2019019
    [7] Abdennasser Chekroun, Mohammed Nor Frioui, Toshikazu Kuniya, Tarik Mohammed Touaoula . Global stability of an age-structured epidemic model with general Lyapunov functional. Mathematical Biosciences and Engineering, 2019, 16(3): 1525-1553. doi: 10.3934/mbe.2019073
    [8] Christine K. Yang, Fred Brauer . Calculation of R0 for age-of-infection models. Mathematical Biosciences and Engineering, 2008, 5(3): 585-599. doi: 10.3934/mbe.2008.5.585
    [9] Chayu Yang, Jin Wang . Computation of the basic reproduction numbers for reaction-diffusion epidemic models. Mathematical Biosciences and Engineering, 2023, 20(8): 15201-15218. doi: 10.3934/mbe.2023680
    [10] Mostafa Adimy, Abdennasser Chekroun, Claudia Pio Ferreira . Global dynamics of a differential-difference system: a case of Kermack-McKendrick SIR model with age-structured protection phase. Mathematical Biosciences and Engineering, 2020, 17(2): 1329-1354. doi: 10.3934/mbe.2020067
  • 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,nR0 as n+. Some numerical simulations are provided to certify the theoretical results.


    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=a0k(σ)eσ0μ(η)dη1γ(1eγσ)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.

    In this section, we first present the age-structured SIRS epidemic model developed by [11],

    {(t+a)S(t,a)=μ(a)S(t,a)λ(a,t)S(t,a)+γ(a)R(t,a),(t+a)I(t,a)=λ(a,t)S(t,a)(μ(a)+ν(a)+δ(a))I(t,a),(t+a)R(t,a)=ν(a)I(t,a)(μ(a)+γ(a))R(t,a),S(t,0)=Λ,t[0,+),S(0,a)=S0(a),a(0,A)I(t,0)=0,t[0,+),I(0,a)=I0(a),a(0,A)R(t,0)=0,t[0,+),R(0,a)=R0(a),a(0,A) (2.1)

    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

    λ(a,t)=A0β(a,ϱ)I(ϱ,t)dϱ.

    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.

    Table 1.  Meanings of all parameters.
    Parameters Meanings
    μ(a) the natural mortality of the population
    β(a,ϱ) the age-dependent transmission coefficient
    γ(a) the rate of removed individuals who lose immunity returning to the susceptible class
    A the maximum age
    ν(a) the natural recovery rate of the infective individuals
    δ(a) the disease inducing death rate

     | Show Table
    DownLoad: CSV

    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

    Tφ(a):=[T1φ1(a)T2φ2(a)T3φ3(a)]=[dφ1(a)daμ(a)φ1(a)λ(a,t)φ1(a)dφ2(a)da(μ(a)+ν(a)+δ(a))φ2(a)dφ3(a)da(μ(a)+γ(a))φ3(a)], (2.2)

    φ(a)=(φ1(a),φ2(a),φ3(a))D(T), where the domain D(T) is given as

    D(T):={φX:φiis absolutely continuous on[0,A],ddaφiXandφ(0)=(0,0,0)}.

    The disease-free equilibrium of model (2.1) is E=(E0(a),0,Er(a)), where Er(a)=ea0(μ(η)+γ(η))dη, and E0(a)=γ(a)Er(a)a0eaϱμ(η)dηdϱ is the density of the susceptible population at age a in the disease-free state. Then we define a nonlinear operator F:XX by

    Fφ(a):=[F1φ1(a)F2φ2(a)F3φ3(a)]=[γ(a)φ3(a)E0(a)A0β(a,ϱ)φ2(ϱ)dϱν(a)φ2(a)]. (2.3)

    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

    ddtu(t)=Tu(t)+Fu(t),u(0)=u0X. (2.4)

    Next, we mainly consider the second equation of (2.1). By simple calculation, the positive inverse (T2)1 is defined as follows

    (T2)1φ2(a):=a0eaϱ(μ(η)+ν(η)+δ(η))dηφ2(ϱ)dϱ,φ2Y:=L1(0,A).

    Then, according to [11], we can give the next generation operator K by

    Kφ2(a):=F2(T2)1φ2(a)=E0(a)A0β(a,ϱ)ϱ0eϱρ(μ(η)+ν(η)+δ(η))dηφ2(ρ)dρdϱ.

    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),

    D(B):={Y:is absolutely continuous on[0,A],ddaYand(0)=0}.

    Hence, we discretize the following system

    ddtI(t)=BI(t)+GI(t),I(0)=I0Y (2.5)

    into a system of ordinary differential equations in Yn:=Rn, nN. 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

    ddtI(t)=BnI(t)+GnI(t),I(0)=I0Yn, (2.6)

    where I(t) and I0 are ncolumn vectors, Bn and Gn are nsquare matrices with the following form

    Bn:=[θM11Δa001Δa(1θ)M1θM21Δa001Δa(1θ)Mn1θMn1Δa]n×n,
    Gn:=[N0[(1θ)β01+θβ11]ΔaN0[(1θ)β02+θβ12]ΔaN0[(1θ)β0n+θβ1n]ΔaN1[(1θ)β11+θβ21]ΔaN1[(1θ)β12+θβ22]ΔaN0[(1θ)β1n+θβ2n]ΔaNn1[(1θ)βn1,1+θβn1]ΔaNn1[(1θ)βn1,2+θβn2]ΔaNn1[(1θ)βn1,n+θβnn]Δa],

    where Mi=μi+νi+δi(i=1,,n), Ni=(1θ)E0i+θE0i+1(i=0,,n1). 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:YYn and J:YnY as follows

    {(Pn)k:=1Δaak+1ak(a)da,k=0,1,,n1,Y,(Jnψ)(a):=n1k=0ψkχ(ak,ak+1](a),ψ=(ψ1,ψ2,,ψn)Yn, (2.7)

    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

    χ(ak,ak+1](a)={1,a(ak,ak+1],0,a(ak,ak+1].

    From Section 4.1 in [19], we know that for all nN, Pn1 and Jn1. We denote Yn is the norm in Yn, and

    ψYn:=Δan1k=0|ψk|,ψ=(ψ1,ψ2,,ψn)Yn. (2.8)

    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+JnKnPnKY=0, then R0,nR0 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

    (Bn)1=[1θM1+1Δa00(1)3((1θ)M11Δa)(θM1+1Δa)(θM2+1Δa)1θM2+1Δa0(1)n+1n1i=1((1θ)Mi1Δa)nk=1(θMk+1Δa)n1i=2(1Δa(1θ)Mi)nk=2(θMk+1Δa)1θMn+1Δa], (2.9)

    then we have

    KnψYn=Gn(Bn)1ψYnΔan1k=0ˉE0ˉβΔaθ(μ_+ν_+δ_)n1k=0|ψk|=AˉE0ˉβθ(μ_+ν_+δ_)ψYn,θ[12,1],

    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,

    limh0A0|E0(a+h)β(a+h,ϱ)E0(a)β(a,ϱ)|da=0uniformly forϱR, (2.10)

    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+JnKnPnKY=0.

    Proof. For any Y, we obtain

    JnKnPnKY=JnGn(Bn)1PnG(B)1YJnGn(Bn)1PnJnGnPn(B)1Y+JnGnPn(B)1G(B)1YJnGn(Bn)1PnPn(B)1Yn+JnGnPn(B)1G(B)1YL(Bn)1PnPn(B)1Yn+JnGnPn(B)1G(B)1Y. (2.11)

    Since Jn1, and for any nN, GnA¯E0ˉβ, we have L=JnGn=A¯E0ˉβ. Next we estimate the first term in the right-hand of (2.11), then

    (Bn)1Pn(B)1PnXn=(Bn)1Pn(B)(B)1(Bn)1(Bn)Pn(B)1Yn(Bn)1Pn(B)(B)1(Bn)Pn(B)1YnAPn(B)ϕ(Bn)PnϕYn,

    where ϕ:=(B)1D(B), and for any ψ=(ψ1,ψ2,,ψn)Yn,

    (Bn)1ψYnΔank=11θ(μ_+ν_+δ_)+1Δan1k=0|ψk|AψYn,

    namely, (Bn)1A. From (2.7), we obtain

    (Bn)1Pn(B)1PnYnAPn(B)ϕ(Bn)PnϕYnAΔan1k=0|Pn(B)ϕ(Bn)Pnϕ|AΔan1k=0|1Δaak+1ak(ddaϕ(a)+(μ(a)+ν(a)+δ(a))ϕ(a))da(1θ)(μ(k)+ν(k)+δ(k))Δaakak1ϕ(a)da1Δaak+1akϕ(a)da1Δaakak1ϕ(a)daΔaθ(μ(k+1)+ν(k+1)+δ(k+1))Δaak+1akϕ(a)da|,

    where a0=a1=0. By the mean value theorem, we have

    (Bn)1Pn(B)1PnYnAΔan1k=0|ddaϕ(ηk+1)+(μ(ηk+1)+ν(ηk+1)+δ(ηk+1))ϕ(ηk+1)(1θ)(μ(k)+ν(k)+δ(k))ϕ(ρk)1Δa(ϕ(ξk+1)ϕ(ξk))θ(μ(k+1)+ν(k+1)+δ(k+1))ϕ(ζk+1)|AΔan1k=0(|ddaϕ(ηk+1)ddaϕ(εk+1)|+|(μ(ηk+1)+ν(ηk+1)+δ(ηk+1))ϕ(ηk+1)(μ(k)+ν(k)+δ(k))ϕ(ϱk)|)+|θ(μ(k)+ν(k)+δ(k))ϕ(ρk)θ(μ(k+1)+ν(k+1)+δ(k+1))ϕ(ζk+1)|AΔan1k=0[ω(ϕ,2Δa)+ω(μ+ν+δ,Δa)ω(ϕ,Δa)+ω(θ(μ+ν+δ),2Δa)ω(ϕ,2Δa)],

    where ω(f,r) denotes the modulus of continuity. We know that ω(f,r) is defined by sup|xy|r|f(x)f(y)| with the following property

    ω(f,r)0,asr0.

    Hence, (Bn)1Pn(B)1PnYn0 holds. Then we consider the second term of (2.11) as follows

    JnGnPn(B)1G(B)1Y=JnGnPnϕGϕY=n1k=0ak+1ak|nj=1((1θ)E0k+θE0k+1)((1θ)βkj+θβk+1,j)jj1ϕ(ϱ)dϱA0E0(a)β(a,ϱ)ϕ(ϱ)dϱ|da=n1k=0ak+1ak|nj=1((1θ)E0k+θE0k+1)((1θ)βkj+θβk+1,j)jj1ϕ(ϱ)dϱnj=1jj1[(1θ)E0+θE0][(1θ)β+θβ]ϕ(ϱ)dϱ|dan1k=0ak+1aknj=1jj1|(1θ)2E0kβkj+(1θ)θE0kβk+1,j+θ(1θ)E0k+1βkj+θ2E0k+1βk+1,j(1θ)2E0β+(1θ)θE0β+θ(1θ)E0β+θ2E0β||ϕ(ϱ)|dϱdan1k=0ak+1aknj=1jj1|(1θ)2ω(E0,Δa)ω(β,Δa)+(1θ)θω(E0,Δa)ω(β,Δa)+θ(1θ)ω(E0,Δa)ω(β,Δa)+θ2ω(E0,Δa)ω(β,Δa)||ϕ(ϱ)|dϱdaAω(E0,Δa)ω(β,Δa)ϕY0asn+, (2.12)

    where ω(E0,Δa)0(Δa0) and ω(β,Δa)0(Δa0), respectively. Hence,

    JnGnPn(B)1G(B)1Y0.

    Combine the above discussion, we have limn+JnKnPnKY=0.

    By virtue of Assumption 2.1 and Lemma 2.1, we know that Theorem 2.1 holds. Namely, R0,nR0 as n+, preserving algebraic multiplicity 1.

    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

    {(t+a)S(t,a)=μ(a)S(t,a)λ(a,t)S(t,a)+γ(a)R(t,a)+σS(t,a)dBtdt,(t+a)I(t,a)=λ(a,t)S(t,a)(μ(a)+ν(a)+δ(a))I(t,a)+σI(t,a)dBtdt,(t+a)R(t,a)=ν(a)I(t,a)(μ(a)+γ(a))R(t,a),S(t,0)=Λ,t[0,+),S(0,a)=S0(a),a(0,A)I(t,0)=0,t[0,+),I(0,a)=I0(a),a(0,A)R(t,0)=0,t[0,+),R(0,a)=R0(a),a(0,A). (2.13)

    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

    {(t+a)I(t,a)=E0(a)A0β(a,ϱ)I(t,a)dϱ(μ(a)+ν(a)+δ(a))I(t,a)+σ(a)I(t,a)dBtdt,I(t,0)=0,t[0,+),I(0,a)=I0(a),a(0,A). (2.14)

    According to the general definition of the stochastic basic reproduction number, the following two operators are defined on Y:=L1(0,A)

    {A(a)=dda(a)(μ(a)+ν(a)+δ(a))(a),F(a)=E0(a)(1σ22)A0β(a,ϱ)(ϱ)dϱ,1σ22>0, (2.15)

    and

    D(A):={Y:is absolutely continuous on[0,A],ddaYand(0)=0}.

    Using A and F to rewrite (2.14) as

    ddtI(t)=AI(t)+FI(t),I(0)=I0. (2.16)

    Then we have

    (A)1(a):=a0eaϱ(μ(η)+ν(η)+δ(η))dη(ϱ)dϱ,Y.

    The next generation operator T is shown by

    T(a):=F(A)1(a)=E0(a)(1σ22)A0β(a,ϱ)ϱ0eϱρ(μ(η)+ν(η)+δ(η))dη(ρ)dρdϱ,Y.

    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, nN. Then the system (2.16) is discretized into the following equation

    ddtI(t)=AnI(t)+FnI(t),I(0)=I0Yn, (2.17)

    where An is defined as the same as Bn(An:=Bn), and

    Fn:=[Q0[(1θ)β01+θβ11]ΔaQ0[(1θ)β0n+θβ1n]ΔaQn1[(1θ)βn1,1+θβn1]ΔaQn1[(1θ)βn1,n+θβnn]Δa]

    where θ[12,1], Qi=(1σ22)[(1θ)E0i+θE0i+1](i=0,1,,n1).

    Theorem 2.2. From Theorem 2.1, we know that T is irreducible, compact and strictly positive. If

    limΔ0EJnTnPnTY=0

    for any Y, then

    Rs0,nRs0asΔ0,preserving algebraic multiplicity 1,

    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

    TnψYn=Fn(An)1ψYnΔan1k=0ˉE0(1σ22)ˉβΔaθ(μ_+ν_+δ_)n1k=0|ψk|=AˉE0(1σ22)ˉβθ(μ_+ν_+δ_)ψYn,θ[12,1],1σ22>0,

    where E_0 is the lower bound of E0.

    Next, we verify that limΔ0JnTnPnTY=0. For any Y, we have

    JnTnPnTY=JnFn(An)1PnF(A)1YJnFn(An)1PnJnFnPn(A)1Y+JnFnPn(A)1F(A)1YJnFn(An)1PnPn(A)1Yn+JnFnPn(A)1F(A)1YAˉE0ˉβ(An)1PnPn(A)1Yn+JnFnPn(A)1F(A)1Y, (2.18)

    where the first term of (2.18)

    (An)1PnPn(A)1Yn=(Bn)1PnPn(B)1Yn

    is similar to the first term in the right-hand of (2.11), so it is easy to see that

    (An)1PnPn(A)1Yn0.

    Next we estimated the second term of (2.18). Let ϖ:=(A)1D(A), we obtain

    JnFnPn(A)1F(A)1Y=JnFnPnϖFϖY=n1k=0ak+1ak|nj=1[(1θ)E0i+θE0i+1](1σ22)jj1ϖ(ϱ)dϱA0E0(a)(1σ22)β(a,ϱ)ϖ(ϱ)dϱ|da=n1k=0ak+1ak|nj=1((1θ)E0k+θE0k+1)(1σ22)((1θ)βkj+θβk+1,j)jj1ϖ(ϱ)dϱnj=1jj1((1θ)E0+θE0)(1σ22)((1θ)β+θβ)ϖ(ϱ)dϱ|da(1σ22)n1k=0ak+1aknj=1jj1|(1θ)2E0kβkj+(1θ)θE0kβk+1,j+θ(1θ)E0k+1βkj+θ2E0k+1βk+1,j(1θ)2E0β+(1θ)θE0β+θ(1θ)E0β+θ2E0β||ϕ(ϱ)|dϱda(1σ22)n1k=0ak+1aknj=1jj1|(1θ)2ω(E0,Δa)ω(β,Δa)+(1θ)θω(E0,Δa)ω(β,Δa)+θ(1θ)ω(E0,Δa)ω(β,Δa)+θ2ω(E0,Δa)ω(β,Δa)||ϕ(ϱ)|dϱdaA(1σ22)ω(E0,Δa)ω(β,Δa)ϕY0asn+. (2.19)

    Thus, JnFnPn(A)1F(A)1Y0 holds. Hence, we obtain the desired assertion.

    In conclusion, Theorem 2.2 holds, which implies that Rs0,nRs0 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 Rn0R0 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.

    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)a0eaϱμ(η)dηdϱ=0.25e(0.45aa42×104)a0eϱ(ϱ32×104+0.2)a(a32×104+0.2)dϱ. Based on numerical integration for E0(a), we obtain Fig. 1(b).

    Figure 1.  Parameters used in the numerical example.

    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)×106+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,nR0(Rs0,nRs0) as n+.

    Let θ=0.5, and choose R0,10002.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 RR0,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].

    Figure 2.  Logarithmic plots of the threshold R0,n (a) and the error RR0,n with respect to the reference value R=2.57673470573749 (b).
    Figure 3.  Computer simulations of the threshold R0,n with different values of θ.

    In this example, let σ=0.1, and we also choose Rs0,10002.56385103220880=:Rs as a reference for Rs0. Similarly, the threshold value Rs0,n for the discretized system (2.17) respect to the reference value Rs (see Fig. 4(a)) and the error RsRs0,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.

    Figure 4.  Logarithmic plots of the threshold Rs0,n (a) and the error RsRs0,n with respect to the reference value Rs=2.56385103220880 (b).
    Figure 5.  Computer simulations of the threshold Rs0,n with different values of σ (a) and the 3D simulation of Rs0,n corresponding to θ[0.5,1] and σ[0,1] (b).

    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.

    The research is supported by the Natural Science Foundation of China (Grant number 11661064).

    The authors declare there is no conflict of interest.



    [1] X. Zhang, D. Jiang, A. Alsaedi, et al., Stationary distribution of stochastic SIS epidemic model with vaccination under regime switching, Appl. Math. Lett., 59 (2016), 87–93.
    [2] P. Driessche and J. Watmough, A simple SIS epidemic model with a backward bifurcation, J. Math. Biol., 40 (2000), 525–540.
    [3] W. Guo, Y. Cai, Q. Zhang, et al., Stochastic persistence and stationary distribution in an SIS epidemic model with media coverage, Physica A, 492 (2018), 2220–2236.
    [4] J. Pan, A. Gray, D. Greenhalgh, et al., Parameter estimation for the stochastic SIS epidemic model, J. Stat. Inference Stoch. Process, 17 (2014), 75–98.
    [5] Y. Cai, Y. Kang and W. Wang, A stochastic SIRS epidemic model with nonlinear incidence rate, Appl. Math. Comput., 305 (2017), 221–240.
    [6] S. Busenberg, M. Iannelli and H. Thieme, Global behavior of an age-structured epidemic model, Siam J. Math. Anal., 22 (1991), 1065–1080.
    [7] B. Cao, H. Huo and H. Xiang, Global stability of an age-structure epidemic model with imperfect vaccination and relapse, Physica A, 486 (2017), 638–655.
    [8] H. Inaba, Age-structured population dynamics in demography and epidemiology, Springer, Sin- gapore, 2017.
    [9] T.Kuniya, Globalstabilityanalysiswithadiscretizationapproachforanage-structuredmultigroup SIR epidemic model, Nonlinear Anal. Real World Appl., 12 (2011), 2640–2655.
    [10] K. Toshikazu, Numerical approximation of the basic reproduction number for a class of age- structured epidemic models, Appl. Math. Lett., 73 (2017), 106–112.
    [11] H. Inaba, Threshold and stability results for an age-structured epidemic model, J. Math. Biol., 28 (1990), 411–434.
    [12] O. Diekmann, J. Heesterbeek and J. Metz, On the definition and the computation of the basic reproduction ratio R 0 , in models for infectious diseases in heterogeneous populations, J. Math. Biol., 28 (1990), 365–382.
    [13] T. Kuniya and R. Oizumi, Existence result for an age-structured SIS epidemic model with spatial diffusion, Nonlinear Anal. Real World Appl., 23 (2015), 196–208.
    [14] N. Bacaër, Approximation of the basic reproduction number R 0 for Vector-Borne diseases with a periodic vector population, B. Math. Biol., 69 (2007), 1067–1091.
    [15] Z. Xu, F. Wu and C. Huang, Theta schemes for SDDEs with non-globally Lipschitz continuous coefficients, J. Comput. Appl. Math., 278 (2015), 258–277.
    [16] X. Mao and L. Szpruch, Strong convergence and stability of implicit numerical methods for stochastic differential equations with non-globally Lipschitz continuous coefficients, J. Comput. Appl. Math., 238 (2013), 14–28.
    [17] F. Chatelin, The spectral approximation of linear operators with applications to the computation of eigenelements of differential and integral operators, Siam Rev., 23 (1981), 495–522.
    [18] A. Berman and R. Plemmons, Nonnegative matrices in the mathematical sciences, Academic press, New York, 1979.
    [19] K. Ito and F. Kappel, The Trotter-Kato theorem and approximation of PDEs, Math. Comput., 67 (1998), 21–44.
    [20] B. Pagter, Irreducible compact operators, Math. Z., 192 (1986), 149–153.
    [21] M. Krein, Linear operators leaving invariant acone in a Banach space, Amer. Math. Soc. Transl.,10 (1962), 3–95.
    [22] W. Guo, Q. Zhang, X. Li, et al., Dynamic behavior of a stochastic SIRS epidemic model with media coverage, Math. Method. Appl. Sci., 41 (2018), 5506-5525.
    [23] C. Mills, J. Robins and M. Lipsitch, Transmissibility of 1918 pandemic influenza, Nature, 432 (2004), 904–906.
  • This article has been cited by:

    1. Dimitri Breda, Francesco Florian, Jordi Ripoll, Rossana Vermiglio, Efficient numerical computation of the basic reproduction number for structured populations, 2021, 384, 03770427, 113165, 10.1016/j.cam.2020.113165
    2. Dimitri Breda, Toshikazu Kuniya, Jordi Ripoll, Rossana Vermiglio, Collocation of Next-Generation Operators for Computing the Basic Reproduction Number of Structured Populations, 2020, 85, 0885-7474, 10.1007/s10915-020-01339-1
    3. Dimitri Breda, Simone De Reggi, Francesca Scarabel, Rossana Vermiglio, Jianhong Wu, Bivariate collocation for computing R0 in epidemic models with two structures, 2022, 116, 08981221, 15, 10.1016/j.camwa.2021.10.026
    4. Huizi Yang, Zhanwen Yang, Shengqiang Liu, Numerical threshold of linearly implicit Euler method for nonlinear infection-age SIR models, 2023, 28, 1531-3492, 70, 10.3934/dcdsb.2022067
    5. Kangkang Chang, Qimin Zhang, Numerical approximation of basic reproduction number for an age‐structured HIV infection model with both virus‐to‐cell and cell‐to‐cell transmissions, 2021, 44, 0170-4214, 12851, 10.1002/mma.7586
    6. X. Liu, M. Zhang, Z.W. Yang, Numerical threshold stability of a nonlinear age-structured reaction diffusion heroin transmission model, 2024, 204, 01689274, 291, 10.1016/j.apnum.2024.06.016
    7. Simone De Reggi, Francesca Scarabel, Rossana Vermiglio, Approximating reproduction numbers: a general numerical method for age-structured models, 2024, 21, 1551-0018, 5360, 10.3934/mbe.2024236
  • Reader Comments
  • © 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4587) PDF downloads(527) Cited by(7)

Figures and Tables

Figures(5)  /  Tables(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog