# A Consistent Discrete Version of a Nonautonomous SIRVS Model

Authors: Joaquim Mateus; César M. Silva; S. Vaz
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

A family of discrete nonautonomous SIRVS models with general incidence is obtained from a continuous family of models by applying Mickens nonstandard discretization method. Conditions for the permanence and extinction of the disease and the stability of disease-free solutions are determined. Concerning extinction and persistence, the consistency of those discrete models with the corresponding continuous model is discussed: if the time step is sufficiently small, when we have extinction (permanence) for the continuous model, we also have extinction (permanence) for the corresponding discrete model. Some numerical simulations are carried out to compare the different possible discretizations of our continuous model using real data.


---

## Body

## 1. Introduction

Most of the epidemiological models in the literature are continuous models. In spite of this, recently, there has been a growing interest in discrete-time models [1–7]. In this work, we will use Mickens nonstandard difference (NSFD) scheme to achieve a discretization of a family of continuous epidemiological models with vaccination and general incidence function considered in [8]. We have multiple objectives: firstly, we want to obtain conditions for extinction and permanence of the disease for the discrete family; next, having a continuous and a corresponding discrete family of models, we wish to discuss the problem of consistency of the discrete models with the corresponding continuous ones; finally, we intend to present some simulation results.The dynamical consistency of a numerical scheme with the associated continuous system is not a precise definition. By the expression “dynamical consistency,” it is meant that the numerical solutions replicate some of the properties of the continuous systems solutions.We will consider the dynamical consistency regarding the permanence and extinction of the disease: whenever there is extinction (permanence) of the disease for the continuous-time model, the same holds for the discrete-time one. Several papers [9–15] discuss the dynamical consistency with respect to some particular properties of discrete epidemiological models obtained from continuous models by some NSFD scheme [16]. We note that while the papers cited above consider autonomous models, in the present work we discuss dynamical consistency for a nonautonomous model. To the best of our knowledge, this is the first work where the consistency of a discretized epidemiological model with the original continuous model is discussed in the nonautonomous context.We remark that, even in the very particular case of autonomous models with mass action incidence and assuming that there is no disease related death, it follows from results in [17] that, when R0>1, the continuous autonomous model has one or two endemic equilibrium points coexisting with the disease-free equilibrium. The existence and stability of the equilibriums are governed by R0 and three additional thresholds: Rq, Rv, and Rqv. These thresholds determine the qualitative behaviour of the system. This very particular situation shows that, even in the autonomous context, the qualitative behaviour can be difficult to determine.One of the motivations for our work was that the difficulties increase considerably when we deal with a general nonautonomous situation. Thus, we decided to discuss an aspect of consistency that, nevertheless, is very important from the point of view of biomathematics: can we obtain consistency between continuous and discrete-time models from the point of view of persistence and extinction of the disease?Regarding our simulation results, we considered two different types of computational experiments. Our first set of simulation results are designed to compare different possible discretizations of our continuous models. After this discussion, we apply our model to a real situation, considering data from the incidence of measles in France in the period 2012-2016. To the possible extent, this data is used to estimate our model parameters and the computational results obtained are compared with real data.The law of mass action states that the rate of change in the disease incidence is directly proportional to the product of the number of susceptible and infective individuals and was the paradigm in the classic models in epidemiology. This is why classical models usually consider a bilinear incidence rateβSI, where S and I denote, respectively, the number of susceptible and infective individuals, to model the disease transmission. In spite of this, it is sometimes important to consider other forms of incidence functions. Another usual assumption is the time independence of the parameters model parameters: in fact, the majority of the epidemiological models in the literature are given by a system of autonomous differential or difference equations. Nevertheless, the assumption that the parameters are independent of time is not very realistic in many situations and it is useful to consider nonautonomous models that, for instance, allow the discussion of environmental and demographic effects that change with time [18, 19]. In this work the family of models considered is nonautonomous and the incidence rates are taken from a large class of functions.Our model generalizes one obtained by Mickens nonstandard finite difference method from the continuous model [8] (see Section 2). In [20], a discrete nonautonomous epidemic model with vaccination and mass action incidence was obtained by Mickens method. We emphasize that, in the particular mass action case, our model is not exactly similar to the model in [20], although Mickens rules were considered in both. We briefly compare computationally these two slightly different models in Section 5.The model we will consider is the following:(1)Sn+1-Sn=Λn-βnφSn+1,In-μn+pnSn+1+ηnVn+1In+1-In=βnφSn+1,In+σnψVn+1,In-μn+αn+γnIn+1Rn+1-Rn=γnIn+1-μnRn+1Vn+1-Vn=pnSn+1-μn+ηnVn+1-σnψVn+1,In,n∈N, where the classes S, I, R, and V correspond, respectively, to susceptible, infective, recovered, and vaccinated individuals and the parameter functions have the following meanings: Λn denotes the inflow of newborns in the susceptible class; the function βnφ is the incidence (into the infective class) from the susceptible individuals; the function σnψ is the incidence (into the infective class) from the vaccinated individuals; μn are the natural deaths; pn represents the susceptibles vaccination; ηn represents the immunity loss and consequence influx in the susceptible class; αn are the deaths occurring in the infective class; γn is the recovery. We will assume that (Λn),(μn),(pn),(ηn),(αn), (βn), (σn), and (γn) are bounded and nonnegative sequences and that there are positive constants wμ, wΛ, wp, kφ, and kψ such that (H1)
the functionsφ:R2→R and ψ:R2→R are nonnegative and differentiable in (R0+)2 and the functions R0+∋x→∂2φ(x,0) and R0+∋x→∂2ψ(x,0) are nondecreasing and Lipschitz, with Lipschitz constants kφ and kψ,(H2)
we haveφ(x,0)=ψ(x,0)=φ(0,y)=ψ(0,y)=0 for all x,y∈R0+,(H3)
limsupn→+∞⁡∏k=nn+ωμ1/1+μk<1,(H4)
liminfn→+∞⁡∑k=n+1n+ωΛΛk>0 and liminfn→+∞⁡∑k=n+1n+ωppk>0,(H5)
functionsR+∋y↦ϕ(x,y)/y and R+∋y↦ψ(x,y)/y are nonincreasing.In this work, we prove, when our conditions prescribe extinction (permanence) for the continuous model we also have extinction (permanence) for the corresponding discrete model as long as the time step is smaller than some constant (that depends on some model parameters and on the threshold condition). We also consider a family of examples of the periodic system of period1 such that the continuous and the discrete-time system with time step h=1/L is not consistent, highlighting the importance of knowing that for time steps smaller than some explicit value we have consistency.


## 2. Discretization of the Continuous Model

We start with a nonautonomous SIRVS model that is slightly less general than the one considered in [8] and generalizes the one in [21]. Namely, we consider the model:(2)S′=Λt-βtφSI-μt+ptS+ηtVI′=βtφS+σtψV-μt-αt-γtIR′=γtI-μtRV′=ptS-μt+ηtV-σtψVI.We assume that the functionsΛ,μ,p,η,α, β, σ, and γ belong to the class C1(R0+) and are nonnegative and bounded. We also require that(C1)
the functionsφ:R→R and ψ:R→R are nonnegative, nondecreasing, differentiable, and Lipschitz with Lipschitz constants kφ and kψ, respectively;(C2)
φ(0)=ψ(0)=0;(C3)
there isω>0 such that liminft→+∞⁡∫tt+ωμ(s)ds>0.In order to obtain threshold conditions for model (2), it was considered in [8] the following auxiliary system:(3)x′=Λt-μt+ptx+ηtyy′=ptx-μt+ηty.And for each solution (x∗(t),y∗(t)) of (3) with positive initial conditions, it was shown that the numbers(4)RClλ=liminft→∞⁡∫tt+λβsφx∗s+σsψy∗s-μs-αs-γsds,RCuλ=limsupt→∞⁡∫tt+λβsφx∗s+σsψy∗s-μs-αs-γsds are independent of the particular solution.Using the above numbers, the following results are contained in results obtained in [8].Theorem 1 (Theorem 1 of [8]).
Assume that conditions (C1), (C2), and (C3) hold. Then, if there is a constantλ>0 such that RCl(λ)>0, then the infectives I are permanent in system (2).Theorem 2 (Theorem 2 of [8]).
Assume that conditions (C1), (C2), and (C3) hold. Then if there is a constantλ>0 such that RCu(λ)<0, then the infectives I go to extinction in system (2).In the literature, several models were discretized using Mickens NSFD schemes [22–35]. Next, we will apply Micken’s nonstandard method to obtain a discrete version of system (2).Letϕ:R0+→R be a positive continuous function such that(5)limh→0⁡ϕh=0.Given h∈R+, we let t=nh, with n∈N, and identify S′(t) with (6)Snh+h-Snhϕh. After deciding a nonlocal representation for the incidence function and that terms that do not correspond to an interaction will be considered in the n+1 time, the first equation in (2) becomes(7)Sn+1h-Snh=ϕhΛnh-βnhφSn+1hInh-μnh+pnhSn+1h+ηnhVnh+h.Writing Sn=S(nh), In=I(nh), Vn=V(nh), Λn=ϕ(h)Λ(nh), βn=ϕ(h)β(nh), μn=ϕ(h)μ(nh), pn=ϕ(h)p(nh), and ηn=ϕ(h)η(nh), we have (8)Sn+1-Sn=Λn-βnφSn+1In-μn+pnSn+1+ηnVn+1. Proceeding similarly for the other equations, we obtain the following discrete model: (9)Sn+1-Sn=Λn-βnφSn+1In-μn+pnSn+1+ηnVn+1In+1-In=βnφSn+1In+σnψVn+1In-μn+αn+γnIn+1Rn+1-Rn=γnIn+1-μnRn+1Vn+1-Vn=pnSn+1-μn+ηnVn+1-σnψVn+1In,n∈N0. We will consider a model that contains this one to obtain some of our results. Namely, based on model (9), in Sections 3 and 4 we will study model (1) that has a more general form for the incidence function.Now, we need to make some definitions. We say that(i)
the infectives(In) arepermanent if for any solution (Sn,In,Rn,Vn) of (1) with initial conditions S0,I0,R0,V0>0 there are constants 0<m<M such that(10)m<liminfn→∞⁡In⩽limsupn→∞⁡In<M;(ii)
the infectives(In) go toextinction if for any solution (Sn,In,Rn,Vn) of (1) with initial conditions S0,I0,R0,V0⩾0 we have limn→∞⁡In=0.Similar definitions can be made for the other compartments. For instance, if there exist constants0<m<M such that for any solution (Sn,In,Rn,Vn) of (1) with initial conditions S0,I0,R0,V0>0 we have(11)m<liminfn→∞⁡Sn⩽limsupn→∞⁡Sn<M we say that the susceptibles are permanent.


## 3. Permanence and Extinction in the Discrete Model

In this section, we will extend the results obtained for the model with the usual mass action incidence in [20] to our generalized family of models. Namely, suitable thresholds are defined and conditions for persistence and extinction of the disease are obtained. As a corollary of our results, we consider the periodic case where we have a unique number that establishes the boundary between the regions of permanence and extinction. Although the proofs of our results are inspired in [20], some difficulties must be dealt with. In particular, it was necessary to understand the right conditions to impose to the incidence functions in order to overcome the technical difficulties.To lighten the reading, the proofs of our results are presented in the appendix.


### 3.1. Auxiliary Results

Consider the auxiliary system,(12)xn+1=Λn+ηnyn+1+xn1+μn+pnyn+1=pnxn+1+yn1+μn+ηn.Note that the auxiliary system describes the behaviour of the system in the absence of infection. If (Λn), (μn), (pn), (ηn), (αn), (μn), (σn), and (βn) are constant sequences, then the linear system (12) becomes autonomous and corresponds to the linearization of the equations for (Sn) and (Vn) in the classical (autonomous) SIRVS model.In order to proceed we need to recall some notions. A solution(un) of some system of difference equations un+1=fn(un) is said to beattractive if for all n0∈N and all ε>0 there is σ(n0)>0 and T(ε,n0,u0)∈N such that if (u¯n) is a solution with u0-u¯0<σ(n0) then un-u¯n<ε, for all n⩾n0+T(ε,n0,u0). Additionally, if some solution is attractive and we can take T to be only dependent on ε, we say that it isuniformly attractive.The following theorem furnishes some simple properties of system (12).Lemma 3 (lemma 2.2 of [20]).
Assume that conditions (H3) and (H4) hold. Then(i)
all solutions(xn,yn) of system (12) with initial condition x0⩾0 and y0⩾0 are nonnegative for all n∈N0;(ii)
each fixed solution(xn,yn) of (12) is bounded and globally uniformly attractive for all n∈N0;(iii)
if(xn,yn) is a solution of (12) and (x~n,y~n) is a solution of the system(13)xn+1=Λn+ηnyn+1+xn+fn1+μn+pnyn+1=pnxn+1+yn+gn1+μn+ηn,with(x~0,y~0)=(x0,y0) then there is a constant L>0, only depending on μn, satisfying (14)supn∈N0⁡x~n-xn+y~n-yn⩽Lsupn∈N0⁡fn+gn;(iv)
there exist constantsm,M>0 such that, for each solution (xn,yn) of (12), we have(15)m⩽liminfn→∞⁡xn⩽limsupn→∞⁡xn⩽M,m⩽liminfn→∞⁡yn⩽limsupn→∞⁡yn⩽M;(v)
when system (12) is ω–periodic, it has a unique positive ω–periodic solution which is globally uniformly attractive.We have the following lemma.Lemma 4.
Assume that condition (H5) holds. Then we have the following:(i)
all solutions(Sn,In,Rn,Vn) of (1) with nonnegative initial conditions are nonnegative for all n∈N0;(ii)
all solutions(Sn,In,Rn,Vn) of (1) with positive initial conditions are positive for all n∈N0;(iii)
there is a constantM>0 such that, if (Sn,In,Rn,Vn) is a solution of (1) with nonnegative initial conditions, then(16)limsupn→+∞⁡Sn+In+Rn+Vn<M.Proof.
S e e the appendix.For eachλ and each particular solution ξn∗=(xn∗,yn∗) of (12) with x0∗>0 and y0∗>0 we define the numbers(17)RDlξ∗,λ=liminfn→∞⁡∏k=nn+λ1+βk∂2φxk+1∗,0+σk∂2ψyk+1∗,01+μk+αk+γk,(18)RDuξ∗,λ=limsupn→∞⁡∏k=nn+λ1+βk∂2φxk+1∗,0+σk∂2ψyk+1∗,01+μk+αk+γk,where ∂if denotes the partial derivative of f with respect to the i-th variable. Contrarily to what one could expect, the next lemma shows that the numbers above do not depend on the particular solution ξn=(xn,yn) of (12) with xn(0)>0 and yn(0)>0.Lemma 5.
Assume that (H1), (H3), and (H4) hold. If(ξ1∗)n=((x1)n∗,(y1)n∗) and (ξ2∗)n=((x2)n∗,(y2)n∗) are two solutions of (12) with xi∗(0)>0 and yi∗(0)>0, i=1,2, then(19)RDlξ1∗,λ=RDlξ2∗,λ,RDuξ1∗,λ=RDuξ2∗,λ.Proof.
Seethe Appendix.By Lemma5 we can drop the dependence of the particular solution and simply write RDl(λ) and RDu(λ) instead of RDl(ξ∗,λ) an RDu(ξ∗,λ), respectively.


### 3.2. Extinction and Permanence

We have the following result about the extinction of the disease.Theorem 6 (extinction of the disease).
Assume that conditions (H1) to (H5) hold. Then(a)
if there is a constantλ>0 such that RDu(λ)<1, then the infectives (In) go to extinction;(b)
any solution(xn∗,0,0,yn∗), where (xn∗,yn∗) is a particular solution of system (12), is globally uniformly attractive.Proof.
S e e   t h e   A p p e n d i x .We have the following result about the permanence of the disease.Theorem 7 (permanence of the disease).
Assume that conditions (H1) to (H5) hold. If there is a constantλ>0 such that RDl(λ)>1 then the infectives (In) are permanent in system (1).Proof.
See the Appendix.We consider now the particular periodic case: assume that all parameters of system (1) are periodic with period ω∈N. By (v) in Lemma 3, there is an ω-periodic disease-free solution of (12), ξ∗=(xn∗,yn∗)n∈N. Thus, in the periodic setting, (17) and (18) become both equal to(20)RDperξ∗=∏k=0ω-11+βk∂2φxk+1∗,0+σk∂2ψyk+1∗,01+μk+αk+γk.Therefore, we obtain the corollary.Corollary 8 (periodic case).
Assume that all coefficients areω-periodic in (1) and that conditions (H1) to (H5) hold. Then (a)
ifRDper(ξ∗)<1 then the infectives (In) go to extinction;(b)
the disease-free solution(xn∗,0,0,yn∗), where (xn∗,yn∗)n∈N is an disease-free ω-periodic solution of (12), is globally attractive;(c)
ifRDper(ξ∗)>1, then the infectives (In) are permanent.Proof.
S e e   t h e Appendix.


## 3.1. Auxiliary Results

Consider the auxiliary system,(12)xn+1=Λn+ηnyn+1+xn1+μn+pnyn+1=pnxn+1+yn1+μn+ηn.Note that the auxiliary system describes the behaviour of the system in the absence of infection. If (Λn), (μn), (pn), (ηn), (αn), (μn), (σn), and (βn) are constant sequences, then the linear system (12) becomes autonomous and corresponds to the linearization of the equations for (Sn) and (Vn) in the classical (autonomous) SIRVS model.In order to proceed we need to recall some notions. A solution(un) of some system of difference equations un+1=fn(un) is said to beattractive if for all n0∈N and all ε>0 there is σ(n0)>0 and T(ε,n0,u0)∈N such that if (u¯n) is a solution with u0-u¯0<σ(n0) then un-u¯n<ε, for all n⩾n0+T(ε,n0,u0). Additionally, if some solution is attractive and we can take T to be only dependent on ε, we say that it isuniformly attractive.The following theorem furnishes some simple properties of system (12).Lemma 3 (lemma 2.2 of [20]).
Assume that conditions (H3) and (H4) hold. Then(i)
all solutions(xn,yn) of system (12) with initial condition x0⩾0 and y0⩾0 are nonnegative for all n∈N0;(ii)
each fixed solution(xn,yn) of (12) is bounded and globally uniformly attractive for all n∈N0;(iii)
if(xn,yn) is a solution of (12) and (x~n,y~n) is a solution of the system(13)xn+1=Λn+ηnyn+1+xn+fn1+μn+pnyn+1=pnxn+1+yn+gn1+μn+ηn,with(x~0,y~0)=(x0,y0) then there is a constant L>0, only depending on μn, satisfying (14)supn∈N0⁡x~n-xn+y~n-yn⩽Lsupn∈N0⁡fn+gn;(iv)
there exist constantsm,M>0 such that, for each solution (xn,yn) of (12), we have(15)m⩽liminfn→∞⁡xn⩽limsupn→∞⁡xn⩽M,m⩽liminfn→∞⁡yn⩽limsupn→∞⁡yn⩽M;(v)
when system (12) is ω–periodic, it has a unique positive ω–periodic solution which is globally uniformly attractive.We have the following lemma.Lemma 4.
Assume that condition (H5) holds. Then we have the following:(i)
all solutions(Sn,In,Rn,Vn) of (1) with nonnegative initial conditions are nonnegative for all n∈N0;(ii)
all solutions(Sn,In,Rn,Vn) of (1) with positive initial conditions are positive for all n∈N0;(iii)
there is a constantM>0 such that, if (Sn,In,Rn,Vn) is a solution of (1) with nonnegative initial conditions, then(16)limsupn→+∞⁡Sn+In+Rn+Vn<M.Proof.
S e e the appendix.For eachλ and each particular solution ξn∗=(xn∗,yn∗) of (12) with x0∗>0 and y0∗>0 we define the numbers(17)RDlξ∗,λ=liminfn→∞⁡∏k=nn+λ1+βk∂2φxk+1∗,0+σk∂2ψyk+1∗,01+μk+αk+γk,(18)RDuξ∗,λ=limsupn→∞⁡∏k=nn+λ1+βk∂2φxk+1∗,0+σk∂2ψyk+1∗,01+μk+αk+γk,where ∂if denotes the partial derivative of f with respect to the i-th variable. Contrarily to what one could expect, the next lemma shows that the numbers above do not depend on the particular solution ξn=(xn,yn) of (12) with xn(0)>0 and yn(0)>0.Lemma 5.
Assume that (H1), (H3), and (H4) hold. If(ξ1∗)n=((x1)n∗,(y1)n∗) and (ξ2∗)n=((x2)n∗,(y2)n∗) are two solutions of (12) with xi∗(0)>0 and yi∗(0)>0, i=1,2, then(19)RDlξ1∗,λ=RDlξ2∗,λ,RDuξ1∗,λ=RDuξ2∗,λ.Proof.
Seethe Appendix.By Lemma5 we can drop the dependence of the particular solution and simply write RDl(λ) and RDu(λ) instead of RDl(ξ∗,λ) an RDu(ξ∗,λ), respectively.


## 3.2. Extinction and Permanence

We have the following result about the extinction of the disease.Theorem 6 (extinction of the disease).
Assume that conditions (H1) to (H5) hold. Then(a)
if there is a constantλ>0 such that RDu(λ)<1, then the infectives (In) go to extinction;(b)
any solution(xn∗,0,0,yn∗), where (xn∗,yn∗) is a particular solution of system (12), is globally uniformly attractive.Proof.
S e e   t h e   A p p e n d i x .We have the following result about the permanence of the disease.Theorem 7 (permanence of the disease).
Assume that conditions (H1) to (H5) hold. If there is a constantλ>0 such that RDl(λ)>1 then the infectives (In) are permanent in system (1).Proof.
See the Appendix.We consider now the particular periodic case: assume that all parameters of system (1) are periodic with period ω∈N. By (v) in Lemma 3, there is an ω-periodic disease-free solution of (12), ξ∗=(xn∗,yn∗)n∈N. Thus, in the periodic setting, (17) and (18) become both equal to(20)RDperξ∗=∏k=0ω-11+βk∂2φxk+1∗,0+σk∂2ψyk+1∗,01+μk+αk+γk.Therefore, we obtain the corollary.Corollary 8 (periodic case).
Assume that all coefficients areω-periodic in (1) and that conditions (H1) to (H5) hold. Then (a)
ifRDper(ξ∗)<1 then the infectives (In) go to extinction;(b)
the disease-free solution(xn∗,0,0,yn∗), where (xn∗,yn∗)n∈N is an disease-free ω-periodic solution of (12), is globally attractive;(c)
ifRDper(ξ∗)>1, then the infectives (In) are permanent.Proof.
S e e   t h e Appendix.


## 4. Consistency

In this section, under the additional assumption that the parameter functionsΛ, μ, η, and p are constant, we will get a result stating that when our integral conditions prescribe extinction (persistence) for the continuous-time model, then the discrete-time conditions prescribe extinction (persistence) for the corresponding discrete-time models, as long as the time step is less than some constant. Throughout this section, we assume that the parameter functions Λ, μ, η, and p are constant functions and ϕ(h) will be the function used in the discretization of the derivative.We consider the continuous-time model (2) and, for a given time step h, the corresponding discrete-time model, that is, the discrete-time model with parameters βkh=ϕ(h)β(kh), σkh=ϕ(h)σ(kh), Λkh=ϕ(h)Λ, μkh=ϕ(h)μ, pkh=ϕ(h)p, ηkh=ϕ(h)η, αkh=ϕ(h)α(kh), and γkh=ϕ(h)γ(kh).For a given time steph>0, the expressions RDl(λ) and RDu(λ) in (17) and (18) become, in our context,(21)RDlλ,h=liminfn→∞∏k=nn+λ1+βkhφxk+1∗+σkhψyk+1∗1+μkh+αkh+γkh,RDuλ,h=limsupn→∞∏k=nn+λ1+βkhφxk+1∗+σkhψyk+1∗1+μkh+αkh+γkh, where (xk∗,yk∗) is the solution of the (in our context autonomous) system (12).We have the following result.Theorem 9.
For system (2), assume that Λ(t)=Λ, μ(t)=μ, η(t)=η, and p(t)=p for all t⩾0 and that the functions α(t), γ(t), β(t), and σ(t) are differentiable, nonnegative, and bounded and have bounded derivative. Assume also that conditions (C1) to (C3) hold and let (22)hmaxu=-RCuλsupt⩾0⁡f′tλ+1,hmaxl=RClλsupt⩾0⁡f′tλ+1,where (23)ft=βtφΛμ+ημμ+η+p+σtψpΛμμ+η+p-μ-αt-γt. Then,(a)
ifRCu(λ)<0, then RDu(⌊λ/h⌋,h)<1 for all h∈]0,hmaxu[;(b)
ifRCl(λ)>0, then RDl(⌊λ/h⌋,h)>1 for all for all h∈]0,hmaxl[.Proof.
Observe that(xn,yn)=(a,b), n∈N and (x(t),y(t))=(a,b), t∈R, where (24)a,b=Λμ+ημμ+η+p,pΛμμ+η+p, are, respectively, solutions of system (12) and system (3). Thus,(25)RCuλ=limsupt→∞⁡∫tt+λβsφa+σsψb-μs+αs+γsds,RDuλh,h=limsupn→∞∏k=nn+λ/h1+βkhφa+σkhψb1+μkh+αkh+γkh.
By contradiction, assume that(26)RCuλ<0,and that there is a sequence (hm)m∈N such that hm→0 as m→+∞ and(27)RDuλhm,hm=limsupn→∞∏k=nn+λ/hm1+βkhmφa+σkhmψb1+μkhm+αkhm+γkhm⩾1,for all m∈N. By (27), we conclude that, for each m∈N, there are sequences (hm)m∈N and (nm,r)r∈N such that hm→0 as m→+∞, nm,r→+∞ as r→+∞ and(28)∏k=nm,rnm,r+λ/hm1+βkhmφa+σkhmψb>1-hm∏k=nm,rnm,r+λ/hm1+μkhm+αkhm+γkhm.By (28), we have(29)∑k=nm,rnm,r+λ/hmβkhmφa+σkhmψb-μkhm-αkhm-γkhm>Bnm,r,λ,hm-Anm,r,λ,hm-Cnm,r,λ,hmhm,where(30)An,L,h≔-h+h∏k=nn+L/h1+βkhφa+σkhψb-h∑k=nn+L/hβkhφa+σkhψb,(31)Bn,L,h≔-h+h∏k=nn+L/h1+μkh+αkh+γkh-h∑k=nn+L/hμkh+αkh+γkh,(32)Cn,L,h=h∏k=nn+L/h1+μkh+αkh+γkh.and, multiplying both sides by hm, we get (33)ϕhm∑k=nm,rnm,r+λ/hmhmβkhmφa+σkhmψb-μkhm-αkhm-γkhm>Bnm,r,λ,hm-Anm,r,λ,hm-Cnm,r,λ,hm.
We also have(34)Anm,r,λ,hm⩽hm∑k=2λ/hmλhmkβuφa+σuψbkϕhmk⩽hm∑k=0λ/hmλhmkβuφa+σuψbkϕhmk=hm1+βuφa+σuψbϕhmλ/hm.Noting that, by (5), we have(35)limm→+∞⁡1+βuφa+σuψbϕhmλ/hm=limm→+∞⁡1+βuφa+σuψb1/ϕhm1/ϕhmϕhmλ/hm=eβuφa+σuψbλ and that a convergent sequence is bounded; by (34) there is C1>0 such that(36)Anm,r,λ,hm⩽C1hm.Similarly, we have(37)Bnm,r,λ,hm⩽∑k=2λ/hmλhmkμu+αu+γukϕhmk⩽∑k=0λ/hmλhmkμu+αu+γukϕhmk=hm1+μu+αu+γuϕhmλ/hm.Using (5) again, we get(38)limm→+∞⁡1+μu+αu+γuϕhmλ/hm⩽eμu+αu+γuλ. There is C2>0 such that(39)Bnm,r,λ,hm⩽C2hm.Finally, we have(40)Cnm,r,λ,hm=hm∏k=nm,rnm,r+λ/hm1+μkhm+αkhm+γkhm=hm1+3ϕhmmax⁡μu,αu,γuλ/hm+1. According to (5), we obtain(41)limm→+∞⁡1+3ϕhmmax⁡μu,αu,γuλ/hm+1=limm→+∞1+3max⁡μu,αu,γu1/ϕhm1/ϕhmϕhmλ/hm+1=e3max⁡μu,αu,γuλ. There is C3>0 such that(42)Cnm,r,λ,hm⩽C3hm.
Thus,(43)Bnm,r,λ,hm-Anm,r,λ,hm-Cnm,r,λ,hm⩽Anm,r,λ,hm+Bnm,r,λ,hm+Cnm,r,λ,hm⩽C1+C2+C3hm,for all m⩾M. Since the right hand side of (43) is independent of nm,r, we conclude that(44)Bnm,r,λ,hm-Anm,r,λ,hm-Cnm,r,λ,hm⟶0,as m→+∞, uniformly in r.
On the other hand, we note that theC1 function f:R0+→R given by(45)ft=βtφa+σtψb-μt-αt-γtis Riemann-integrable on any bounded interval I⊂R0+.
We have that(46)∑k=nm,rnm,r+λ/hmhmfkhm+λ-λhmhmfnm,rhm+λhmhm is a Riemann sum of (47)∫nm,rhmnm,rhm+λβsφa+σsψb-μs+αs+γsds with respect to the partition (48)nm,rhm,nm,rhm+hm,…,nm,rhm+λhmhm,nm,rhm+λ of size hm of the interval [nm,r,nm,r+λ]. Note that(49)sm,r≔λ-λhmhmfnm,rhm+λhmhm⩽hmfu≔sm and sm→0 as m→+∞, uniformly in r.
Sincef is C1 with bounded derivative, for any h>0 we have (50)fx-fx+h⩽Ch, where C=supt⩾0⁡|f′(t)|. We conclude that(51)∑k=nm,rnm,r+λ/hmhmfkhm+sm,r-∫nm,rhmnm,rhm+λfsds<Chm2λhm+Chm2<Cλ+1hm,thus(52)∑k=nm,rnm,r+λ/hmhmfkhm<∫nm,rhmnm,rhm+λfsds-sm,r+Cλ+1hm, and therefore (53)∑k=nm,rnm,r+λ/hmϕhmhmfkhm<ϕhm∫nm,rhmnm,rhm+λfsds+Cλ+1hm.By (53) we conclude that, given δ>0, there is rm∈N such that, for all r⩾rm,(54)ϕhm∑k=nm,rmnm,rm+λ/hmhmfkhm<ϕhmRCuλ+δ+Cλ+1hm.
Finally, recalling thatRCu(λ)<0, by assumption, by the arbitrariness of δ>0 and the fact that hm→0 as m→+∞, we obtain for sufficiently large m∈N, (55)0⩽ϕhm∑k=nm,rmnm,rm+λ/hmhmfkhm<0, which is a contradiction. We obtain (a).
A similar argument allows us to prove (b). In fact, assuming by contradiction that(56)RClλ>0,and that there is a sequence (hm)m∈N such that hm→0 as m→+∞ and(57)RDlλhm,hm=liminfn→∞∏k=nn+λ/hm1+βkhmφa+σkhmψb1+μkhm+αkhm+γkhm⩽1, it is possible to conclude that(58)ϕhm∑k=nm,rnm,r+λ/hmhmβkhmφa+σkhmψb-μkhm-αkhm-γkhm<Bnm,r,λ,hm-Anm,r,λ,hm+Cnm,r,λ,hm,where Anm,r,λ,hm, Bnm,r,λ,hm, and Cnm,r,λ,hm are given, respectively, by (30), (31), and (42) and still satisfy (36), (39), and (42). Consequently, given δ>0, there is rm∈N such that, for all r⩾rm,(59)ϕhm∑k=nm,rnm,r+λ/hmhmfkhm>ϕhm∫nm,rhmnm,rhm+λfsds-δ+Cλ+1hm.Recalling that RCl(λ)>0, by assumption and since δ>0 is arbitrary, we obtain for sufficiently large m∈N, (60)0⩾ϕhm∑k=nm,rmnm,rm+λ/hmfkhm>0, which is a contradiction. We obtain (b) and the theorem follows.Next, for eachL∈N, we give an example of a periodic system of period 1 such that the continuous and the discrete-time system with time step h=1/L are not consistent; namely, we will have persistence for the continuous-time model and extinction for the discrete-time model with time step h=1/L.Example 10.
LetL∈N. Consider in system (2) that ϕ(x)=ψ(x)=x, that, with the exception of σ and β, all parameters are constant, that Λ=μ, and that (61)σt=βt=d1+csin22πLt1+cos⁡2πt. We obtain a periodic system of period 1.
In this context,(xn,yn)=(a,b), n∈N, and (x(t),y(t))=(a,b), t∈R, where (62)a,b=μ+ημ+η+p,pμ+η+p, are, respectively, solutions of system (12) and system (3). It is now possible to compute the number RCl(1). In fact, noting that x∗(t)+y∗(t)=1, we get(63)RCl1=∫01βsx∗s+σsy∗s-μ-α-γds=∫01d1+csin22πLt1+cos⁡2πtds-μ-α-γ=d1+c2-μ-α-γ. We can also compute RDl(1,1/L). Namely, we have (64)RDl1,1L=1+d/L1+μ+α+γ.
If we letd be sufficiently small so that d<(μ+α+γ)L, or in other words, d<(μ+α+γ) and c be sufficiently large so that c>2/d(μ+γ+α-d), we obtain (65)RCl1>1⟺1+d1+c/21+μ+α+γ>1,RDl1,1L<1⟺1+d/L1+μ+α+γ<1. So we conclude that we do not have consistency for time step 1/L.
LetL=6 and consider the continuous model with the following parameters μ=Λ=0.25, γ=0.3, α=0.05, η=0.05, p=2/3, d=0.6, and c=1.5. In Figure 1, we plot function β (or similarly σ) and the component I(t) of the solution of system (2) given by the solver of Mathematica® (that we take to represent the solution of the continuous-time model) and the solution of the discrete-time model (9) with time step 1/6. As can be seen, the infectives are persistent in the continuous-time model but go to extinction in the discrete-time model. We have inconsistency in this case.Figure 1
Left: functionβ; right: inconsistency (time step=1/6).Note that, changingβ(t) and σ(t) slightly, we can construct an example of a periodic system with period 1 where the infectives in the continuous-time model go to extinction but, in the discrete-time model with time step h=1/L, the infectives are persistent.Furthermore, we emphasize that this lack of consistency is not a result of the discretization method used but simply a result of the fact that the time steps lead to a situation where the pointsn/L where the functions β and σ are evaluated (in order to obtain the discrete-time parameters) correspond to minimums of β and σ.


## 5. Simulation

Our objective in this section is twofold. On the one hand, we want to consider different incidence functionsφ, corresponding to different discretizations of our continuous model, and compare the several discrete models obtained. We do this in the first subsection. On the other hand, we want to use our model to describe a real situation. We do this in the second subsection where we consider data from the incidence of measles in France in the period 2012-2016.


### 5.1. Simulation with Several NSFD Schemes

In this subsection we do some simulation to illustrate our results. To begin, we compare our model (1) with mass action incidence (φ(S,I)=SI and ψ(V,I)=VI) with Zhang’s model [20]. We use the following set of parameters: ϕ(h)=h+0.2h2, Λ=0.5, μ(t)=γ(t)=δ(t)=0.3, α(t)=0.05, η=0.05, p=2/3 and (66)βt=σt=b1+0.3cos⁡tπ2. Setting b=0.3 we obtain RCu(4)=-0.6<0 and thus we conclude that we have extinction for the continuous model. Taking time steps equal to 4, 1, and 0.5, we get RDu(0,4)=RDl(0,4)=1, RDu(3,1)=0.644<1, and RD(7,0.5)=0.601<1 and we conclude that we have extinction for time steps 1 and 0.5. For these parameters, we have consistency in the sense of Theorem 9 as long as the time step is less than 0.05. Clearly, there is numerical evidence that there is consistency even for higher time steps. Figure 2 illustrates this situation.Figure 2
S I; ϕ(h)=h+0.2h2; step-size: 4, 1, 0.5.Changingb to 0.9 we obtain RCl(4)=3.4>0 and thus we conclude that we have persistence. Taking time steps equal to 2, 1, and 0.5, we get RDl(1,2)=3.201>1, RDl(3,1)=5.9>1, and RDl(7,0.5)=10.2>1 and we conclude that we have persistence for all these time steps. Figure 3 illustrates this situation. Figures 2 and 3 suggest that numerically our model is slightly better than Zhang’s model, at least for large time steps.Figure 3
S I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Next, we compare our model with the discretized model obtained by Euler method and the output of the Mathematica solver ODE (that uses a Runge-Kutta method). Consideringb=0.3, we get extinction for the continuous-time model, as we already saw. Taking time steps equal to 2, 1, and 0.5, we can see in Figure 4 that for all methods considered and all time steps we have extinction, although the behaviour of our model shadows better the behaviour given by Mathematica’s solver, at least for these time steps.Figure 4
S I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Changingb to 0.9 we already saw that we get persistence for the continuous model. Figure 5 illustrates this situation.Figure 5
S I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Next, we change our incidence function and considerϕ(S,I)=SI/(1+0.7I), maintaining the set of parameters. Letting b=0.3 we have extinction for the continuous model and letting b=0.9 we have persistence for the continuous-time model. Note that the thresholds RCu, RCl, RDu, and RDl are similar to the mass action case. Figures 6 and 7 illustrate this situation.Figure 6
S I / 1 + 0.7 I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Figure 7
S I / 1 + 0.7 I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Doing corresponding simulations and comparisons for our model withϕ(h)=(1-e-0.002h)/(0.002) instead of ϕ(h)=h+0.2h2, we can draw the same conclusions regarding extinction/persistence, relation to Zhang’s model and the model obtained by Euler method.


### 5.2. Simulation with Real Data

In this subsection, we present some simulation regarding measles. This disease is endemic in some countries such as France. In that country, with the measles outbreak in 2011, a vaccination policy that lowered the number of reported cases was introduced. We will focus on measles in France, between 2012 and 2016. For a study concerning the period before 2012, see [36]. For our parameters estimation, we gathered information from several websites. We considered standard incidence functions ψ(Vn+1,In)=Vn+1In/Pn and φ(Sn+1,In)=Sn+1In/Pn, where Pn is the total population. Inspired in the time series for the infectives (https://ecdc.europa.eu), we considered σn=0.03 and βn given by(67)βn=3.8+10sinn+1π6,ifn12⩽52.7,otherwise. The remaining parameters were considered time independent and were inspired in data contained in the websites http://www.worldbank.org, https://data.oecd.org, and http://www.geoba.se. Namely, we took the mortality rate μn=0.0007, the newborns Λn=50000, the disease induced mortality αn=0.000375, the immunity loss ηn=0.001, the vaccination rate pn=0.001, and the recovery rate γn=0.957. We used the initial conditions S0=7.20428×106, I0=106, V0=5.84372×107, and R0=1.81918×104. In Figure 8, we plot the real data for the infectives and the output given by our model.Figure 8
Measles (2012-2016), simulation.Can be seen, in a general way, that our model behaves in the same manner as the real data. It seems that if the vaccination policy in France continues to be very strict, it may decrease the number of cases.


## 5.1. Simulation with Several NSFD Schemes

In this subsection we do some simulation to illustrate our results. To begin, we compare our model (1) with mass action incidence (φ(S,I)=SI and ψ(V,I)=VI) with Zhang’s model [20]. We use the following set of parameters: ϕ(h)=h+0.2h2, Λ=0.5, μ(t)=γ(t)=δ(t)=0.3, α(t)=0.05, η=0.05, p=2/3 and (66)βt=σt=b1+0.3cos⁡tπ2. Setting b=0.3 we obtain RCu(4)=-0.6<0 and thus we conclude that we have extinction for the continuous model. Taking time steps equal to 4, 1, and 0.5, we get RDu(0,4)=RDl(0,4)=1, RDu(3,1)=0.644<1, and RD(7,0.5)=0.601<1 and we conclude that we have extinction for time steps 1 and 0.5. For these parameters, we have consistency in the sense of Theorem 9 as long as the time step is less than 0.05. Clearly, there is numerical evidence that there is consistency even for higher time steps. Figure 2 illustrates this situation.Figure 2
S I; ϕ(h)=h+0.2h2; step-size: 4, 1, 0.5.Changingb to 0.9 we obtain RCl(4)=3.4>0 and thus we conclude that we have persistence. Taking time steps equal to 2, 1, and 0.5, we get RDl(1,2)=3.201>1, RDl(3,1)=5.9>1, and RDl(7,0.5)=10.2>1 and we conclude that we have persistence for all these time steps. Figure 3 illustrates this situation. Figures 2 and 3 suggest that numerically our model is slightly better than Zhang’s model, at least for large time steps.Figure 3
S I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Next, we compare our model with the discretized model obtained by Euler method and the output of the Mathematica solver ODE (that uses a Runge-Kutta method). Consideringb=0.3, we get extinction for the continuous-time model, as we already saw. Taking time steps equal to 2, 1, and 0.5, we can see in Figure 4 that for all methods considered and all time steps we have extinction, although the behaviour of our model shadows better the behaviour given by Mathematica’s solver, at least for these time steps.Figure 4
S I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Changingb to 0.9 we already saw that we get persistence for the continuous model. Figure 5 illustrates this situation.Figure 5
S I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Next, we change our incidence function and considerϕ(S,I)=SI/(1+0.7I), maintaining the set of parameters. Letting b=0.3 we have extinction for the continuous model and letting b=0.9 we have persistence for the continuous-time model. Note that the thresholds RCu, RCl, RDu, and RDl are similar to the mass action case. Figures 6 and 7 illustrate this situation.Figure 6
S I / 1 + 0.7 I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Figure 7
S I / 1 + 0.7 I; ϕ(h)=h+0.2h2; step-size: 2, 1, 0.5.Doing corresponding simulations and comparisons for our model withϕ(h)=(1-e-0.002h)/(0.002) instead of ϕ(h)=h+0.2h2, we can draw the same conclusions regarding extinction/persistence, relation to Zhang’s model and the model obtained by Euler method.


## 5.2. Simulation with Real Data

In this subsection, we present some simulation regarding measles. This disease is endemic in some countries such as France. In that country, with the measles outbreak in 2011, a vaccination policy that lowered the number of reported cases was introduced. We will focus on measles in France, between 2012 and 2016. For a study concerning the period before 2012, see [36]. For our parameters estimation, we gathered information from several websites. We considered standard incidence functions ψ(Vn+1,In)=Vn+1In/Pn and φ(Sn+1,In)=Sn+1In/Pn, where Pn is the total population. Inspired in the time series for the infectives (https://ecdc.europa.eu), we considered σn=0.03 and βn given by(67)βn=3.8+10sinn+1π6,ifn12⩽52.7,otherwise. The remaining parameters were considered time independent and were inspired in data contained in the websites http://www.worldbank.org, https://data.oecd.org, and http://www.geoba.se. Namely, we took the mortality rate μn=0.0007, the newborns Λn=50000, the disease induced mortality αn=0.000375, the immunity loss ηn=0.001, the vaccination rate pn=0.001, and the recovery rate γn=0.957. We used the initial conditions S0=7.20428×106, I0=106, V0=5.84372×107, and R0=1.81918×104. In Figure 8, we plot the real data for the infectives and the output given by our model.Figure 8
Measles (2012-2016), simulation.Can be seen, in a general way, that our model behaves in the same manner as the real data. It seems that if the vaccination policy in France continues to be very strict, it may decrease the number of cases.


## 6. Conclusions

We considered a discretization procedure, based on Mickens NSFD scheme, to get a discrete-time model from a continuous time with vaccination and incidence given by a general function. For a family of models containing the previous discrete-time model, we achieved results on the persistence and the extinction of the disease (Theorems6 and 7). They contain the results of Zhang [20] as a particular case. Our threshold conditions depend on the parameters of the model and of the incidence function derivative, with respect to the infectives, computed on some disease-free solution. This agrees with the continuous counterparts of these results [8].We also considered the problem of establishing the consistency of the continuous-time model and the discrete-time model for small time steps, in the sense that if the time step is small enough when we have persistence (extinction) for the continuous-time model we also have persistence (extinction) for the discrete-time model (at least for situations where Theorems1 and 2 allow us to conclude that we have persistence or extinction). Assuming the differentiability of parameters, our result on this direction, Theorem 9, furnishes an interval [0,a], where a depends only on the parameters of the model and their derivatives, where there is consistency.We present an example of a periodic system of period1 where the continuous and the discrete-time system with time step h=1/L are not consistent. Namely, for that time step, we will have persistence for the continuous-time model and extinction for the discrete-time model. These examples show the importance of knowing that for time steps smaller than some explicit value we have consistency, a type of result like the one in Theorem 9.Finally, we carried out some simulations to illustrate our results. As one might expect our simulations furnish evidence that we may have consistency in intervals whose lengths are several times bigger than the length of the given interval in Theorem9. Additionally, we used our model to describe a real situation, namely, the case of measles incidence in France in the period 2012-2016, and compared our results with the real-time series for the infectives. We found in general that the predictive behaviour of our model is very similar to the real data.We remark that this work is far from giving a complete answer to the problem of consistency between our discrete and continuous-time models. Moreover, in addition to a further discussion of the qualitative behaviour of the continuous and discrete systems, in future work it would be very interesting to discuss the convergence of the methods and to carry out some numerical analysis in our general nonautonomous setting.



---

*Source: 8152032-2018-06-24.xml*