# Global Dynamics of a Periodic SEIRS Model with General Incidence Rate

Authors: Eric Ávila-Vales; Erika Rivero-Esquivel; Gerardo Emilio García-Almeida
Journal: International Journal of Differential Equations
 (2017)

---

## Abstract

We consider a family of periodic SEIRS epidemic models with a fairly general incidence rate of the formSf(I), and it is shown that the basic reproduction number determines the global dynamics of the models and it is a threshold parameter for persistence of disease. Numerical simulations are performed using a nonlinear incidence rate to estimate the basic reproduction number and illustrate our analytical findings.


---

## Body

## 1. Introduction

Epidemiological models in mathematics have been recognized as valuable tools in analyzing the dynamics of an infectious disease nowadays. They are used to describe the spread of disease and also to make control measures known to avoid its persistence, for example, via vaccination terms or treatment terms. These models consider the total population divided into compartments, given by the biological assumptions on the model and represented by functions depending on timet. The most common categories used are susceptible (S), infected (I), recovered (R), exposed (E), quarantined (Q), and vaccinated (V), and the dynamics of model is given by transmission rates from a compartment to another. We have then indicated that the models could be of type SIR, SIRS, SEIR, SEIRS, SEIVR, SEIQV, and so forth.To ensure that the model can give a justified qualitative description of the disease, the choice of the incidence rate plays an important role. An incidence rate is defined as the number of new health related events or cases of a disease in a population exposed to the risk in a given time period. Some examples are the bilinear incidence rate, the saturated incidence rate, or a general incidence rate. The bilinear incidence rate has been repeatedly used by several authors. It is given byβSI, where β is the transmission rate and the product SI represents the contact between infected and susceptible individuals (based on the law of mass action). It was introduced by Kermack and McKendrick [1] in 1927, and even when it is mathematically simple to use, it faces multiple problems and challenges when it is used to describe disease propagation among gregarious animals or persons [2], because it goes to infinity when I becomes larger. In order to improve the modelling process to study the dynamics of infection among a large population, Capasso and Serio [3] in 1978 introduced a saturated incidence rate by studying the Cholera epidemic spread in Bari, given by βSI/1+kI, where β is the transmission rate and k the saturation constant. Unlike the bilinear incidence, saturated incidence does not grow up without a limit, but it goes to a saturation limit as I goes to infinity. Multiple types of saturated incidence have been used in the literature; see, for example, [2] for a list of them. To avoid the use of a single incidence function, the use of a general incidence rate that includes a family of particular functions with similar properties has become a topic of interest by several authors (see, e.g., [4–8]).The basic reproduction (represented byR0) is defined as “the average number of secondary cases produced by a single infected case when it is introduced in a susceptible population” and it has an important role in the study of disease transmission. In biological terms, usually when this number is less than one, the disease is eradicated from population, but when it is greater than one, the infection persists. Mathematically, it is of interest to compute a threshold parameter with the properties of the basic reproduction number. A method to compute this number for certain compartmental disease models is via the next-generation matrix method developed in [9]; however, it is not useful when the model presents time periodic seasonal terms. Authors like [10, 11] have defined its basic reproduction number for periodic models as an average, to give some results about extinction or persistence of infection. However Bacaër and Guernaouni in [12] introduced the definition of basic reproduction number for periodic environments, and, later, Wang and Zhao [13] made a formal definition of it, via the monodromy matrix.In the present work, we focus on a family of SEIRS epidemic models with a time periodic seasonal term, improving the model of Moneim and Greenhalgh in [14], by introducing an incidence rate with a general function taken from [4] and the references therein.We propose the following SEIRS model:(1)dSdt=μN1-p-βtSfI-μ+rtS+δRdEdt=βtSfI-μ+σEdIdt=σE-μ+γIdRdt=μNp+rtS+γI-μ+δR,where N=S+E+I+R is the total population size, with S,E,I,R denoting the fractions of population that are susceptible, exposed, infected, and recovered, respectively. β(t) is the transmission rate and it is a continuous, positive T-periodic function. p (0≤p≤1) is the vaccination rate of all newborn children. r(t) is the vaccination rate of all susceptibles in the population and it is a continuous, positive periodic function with period LT, where L is an integer. μ is the common per capita birth and death rate. σ,γ, and δ are the per capita rates of leaving the latent stage, infected stage, and recovered stage, respectively. It is assumed that all parameters are positive constants.Bai and Zhou in [5] answered some open problems stated in [14], they also showed that their condition is a threshold between persistence and extinction of the disease via the framework established in [13]. They assumed that the incidence was bilinear. In our study, the nonlinear assumptions on function f are listed below (see [4]).(A1)
f:R+→R+ is continuously differentiable.(A2)
f(0)=0, f′(0)>0 and f(I)>0 for all I>0.(A3)
f(I)-If′(I)≥0.Under these assumptions, functionf(I) includes various types of incidence rate; in particular, when f(I)=I, we are on the bilinear case considered in [14].In addition, we assume the following extra conditions (see [15]).(A4)
f′′(0)≤0.(A5)
There existsϵ∗>0 such that when 0<I<ϵ∗, f(I)≥f(0)+If′(0)+1/2I2f′′(0).This set of assumptions on the functionf allows for more general incidence functions than the bilinear one, like saturated incidence functions and functions of the form βSI/(1+kIq); in particular, in the case when q>1, they represent psychological or media effects depending on the infected population. In this last case the incidence function is nonmonotone on I. (A3) regulates the value of f(I) comparing it with the value at I of a line containing the origin of slope f′(I) (note that this line varies as I increases), (A4) requires a concave f(I) at the origin, and (A5) imposes the geometrical condition that in a small neighborhood of the origin f(I) must lie between the tangent line of f at I and a concave parabola tangent to f at I.We consider a family ofSEIRS epidemic models with periodic coefficients and general incidence rate in epidemiology. Then we show that the global dynamics of solutions is determined by the basic reproduction number R0, generalizing the results in [5]. The layout of this paper is as follows: In Section 2, we prove the existence of a disease-free periodic solution and we introduce the basic reproduction number via the theory developed in [12, 13]. In Section 3, we adapt the arguments given in [5] to prove that the disease-free periodic solution of system (1) is globally asymptotically stable if R0<1 and it is persistent when R0>1. Finally, in Section 4, we give some numerical simulations of our results, making a comparison between our basic reproduction number R0 and the average reproduction number R0T used by several authors (see, e.g., [10, 11]).


## 2. The Basic Reproduction Number

First of all, we prove nonnegativity of the solutions under nonnegative initial conditions.Theorem 1.
LetS0,E0,I0,R0≥0. The solution (S(t),E(t),I(t),R(t)) of (1) with (2)S0,E0,I0,R0=S0,E0,I0,R0is nonnegative in the sense that S(t),E(t),I(t),R(t)≥0, ∀t>0, and satisfies S(t)+E(t)+I(t)+R(t)=N, with N constant.Proof.
LetN(t)=S(t)+E(t)+I(t)+R(t); then, adding all equations of system (1), we can see that dN/dt=0, so the value of N is constant. Now, set x(t)=(S(t),E(t),I(t),R(t)) as the solution of system (1) under initial conditions x0=(S(0),E(0),I(0),R(0))=(S0,E0,I0,R0)≥0. By the continuity of solutions, for all of S(t),E(t),I(t) and R(t) that have a positive initial value at t=0, we have the existence of an interval (0,t0) such that S(t),E(t),I(t),R(t)≥0 for 0<t<t0. We will prove that t0=∞.
IfS(t1)=0 for a t1≥0 and other components of x(t) remain nonnegative at t=t1, then (3)dSdtt1=μN1-p+δRt1≥0,implying that whenever the solution x(t) touches the S-axis, the derivative of S is nondecreasing and the function S(t) does not cross to negative values. Similarly, when E(t1)=0 for a t1≥0 and other components remain nonnegative, we have (4)dEdtt1=βt1St1fIt1≥0.
WhenI(t1)=0 for a t1≥0 and other components remain nonnegative, (5)dIdtt1=σEt1≥0.
Finally, whenR(t1)=0 for a t1≥0 and other components remain nonnegative, (6)dRdtt1=μNp+rt1St1+γIt1≥0.
Therefore, wheneverx(t) touches any of the axes S=0, E=0, I=0, or R=0, it never crosses them.In order to make the analysis of the model in a simpler way from now on, we make a reduction of dimension in system (1) making R=N-S-E-I, obtaining the following:(7)dSdt=μN1-p-βtSfI-μ+rtS+δN-S-E-I,dEdt=βtSfI-μ+σE,dIdt=σE-μ+γI.The dynamics of system (1) is equivalent to that of (7); moreover, due to positivity of solutions, we have S+E+I≤N, so we study the dynamics of system (7) in the region(8)X=S,E,I∈R+3:S+E+I≤N.A disease-free periodic solution can be found for (7). To find it, set E=I=0; then, from the first equation of (7) we can obtain the following initial value problem:(9)dSdt=μN1-p-μ+rtS+δN-S,S0=S0∈R+.From [5, 14], the equation above admits a unique positive LT-periodic solution given by(10)S^t=e-∫0tμ+rs+δdsS^0+Nμ1-p+δ∫0te∫0sμ+rξ+δdξds,where (11)S^0=Nμ1-p+δ∫0LTe∫0sμ+rξ+δdξdse∫0LTμ+rs+δds-1.Therefore,(S^(t),0,0) is a disease-free periodic solution of (7); moreover, from [5] we have that S^(t)≤N; therefore, (S^(t),0,0) lives in X.Using the notation of [9], we sort the compartments so that the first two compartments correspond to infected individuals. Let x=(E,I,S) and define(i)
Fi: the rate of new infection in compartment i,(ii)
Vi+: the rate of individuals into compartment i by other means,(iii)
Vi-: the rate of individuals transfer out of compartment i.System can be written as(12)x′t=βtSfI-μ+σEσE-μ+γIμN1-p-βtSfI-μ+rtS+δN-S-E-I=F-V,where V=V--V+,(13)F=βtSfI00,V+=0σEμN1-p+δN,V-=μ+σEμ+γIβtSfI+δS+E+I+μ+rtS.Linearizing system (12) around the disease-free solution, we obtain the matrix of partial derivatives J(0,0,S^)=DF(0,0,S^)-DV(0,0,S^), where(14)DF0,0,S^=0βtS^f′00000000DV0,0,S^=μ+σ00-σμ+γ0δβtS^f′0+δδ+μ+rt.Using Lemma1 of [9], we part DF and DV and set(15)Ft=0βtS^f′000,Vt=μ+σ0-σμ+γ.For a compartmental epidemiological model based on an autonomous system, the basic reproduction number is determined by the spectral radius of the next-generation matrixFV-1 (which is independent of time) [9]. The definition of basic reproduction number for nonautonomous systems has been studied for multiple authors; see, for example, [12, 13]. Particularly, Wang and Zhao in [13] extended the work of [9] to include epidemiological models in periodic environments. They introduced the next infection operator L:CLT→CLT given by(16)Lϕt=∫0∞Yt,t-aFt-aϕt-ada,∀t∈R,ϕ∈CLT,where CLT is the ordered Banach space of all LT periodic functions from R to R2, which is equipped with the maximum norm. ϕ(s)∈CLT is the initial distribution of infectious individuals in this periodic environment, and Y(t,s),t≥s is the evolution operator of the linear periodic system:(17)dydt=-Vty,meaning that, for each s∈R, the 2×2 matrix Y satisfies(18)dYt,sdt=-VtYt,s,∀t≥s,Ys,s=I2×2.L
ϕ is the distribution of accumulative new infections at time t produced by all those infected individuals ϕ(s) introduced before t, with kernel K(t,a)=Y(t,t-a)F(t-a). The coefficient Ki,j(t,a) in row i and column j represents the expected number of individuals in compartment Ii that one individual in compartment Ij generates at the beginning of an epidemic per unit time at time t if it has been in compartment Ij for a units of time, with I1=E,I2=I [16].Letr0>0,r0 is an eigenvalue of L if there is a nonnegative eigenfunction v(t)∈CLT such that(19)Lv=r0v.Therefore, the basic reproduction number is defined as(20)R0≔ρL,the spectral radius of L. The basic reproduction number can be evaluated by several numerical methods and approximations [15–17]; in Section 4 we discuss this topic.


## 3. The Threshold Dynamics ofR0

### 3.1. Disease Extinction

Theorem 2.
LetR0 be defined as (20); then the disease-free periodic solution (S^(t),0,0) is asymptotically stable if R0<1 and unstable if R0>1.Proof.
We use Theorem2.2 of [13] and check conditions (A1)–(A7). Conditions (A1)–(A5) are clearly satisfied from the definitions of F and V given in Section 2. We prove only conditions (A6) and (A7). Define (21)Mt≔-μ+rt+δ,and let ΦM(t) be the monodromy matrix of system(22)dzdt=Mtz.(A6)
ρ(ΦM(LT))<1. Let ΨM be a fundamental matrix for system dz/dt=M(t)z, with M defined as before and LT periodic; the monodromy matrix ΦM(LT) is given by ΦM(LT)=ΨM-1(0)ΨM(LT). The general solution of (22) is (23)zt=Kexp⁡-∫0tμ+rs+δds,soΨM=exp(-∫0t(μ+r(s)+δ)ds) and ΨM-1=exp(∫0t(μ+r(s)+δ)ds). Note that ΨM-1(0)=1, so ΦM(LT)=ΨM(LT) and(24)ΦMLT=exp-∫0LTμ+rs+δds.Due to the fact thatΦM(LT) is a constant, its eigenvalue is itself and ρ(ΦM(LT))<1 for μ,δ,r(s)>0.(A7)
ρ(Φ-V(LT))<1. Solving the system dz/dt=-V(t)z, we arrive at the general solution (25)zt=c1γ-σσ1e-μ+σt+c201e-μ+γt,so(26)Ψ-Vt=γ-σσe-μ+σt0e-μ+σte-μ+γt.ComputingΦ-V(LT)=Ψ-V-1(0)Ψ-V(LT), we have(27)Φ-VLT=e-μ+σLT00e-μ+γLT.Clearly,ρ(Φ-V(LT))=max{e-(μ+σ)LT,e-(μ+γ)LT}<1 for μ,γ,σ>0.Note 1.
Due to the fact thatΨA is a fundamental solution of a periodic system, we can always choose it such that Ψ(0)=I, so the monodromy matrix satisfies ΦA(LT)=ΨA(LT). This property is used in further analysis.In order to prove the global stability of the disease-free periodic solution, we enunciate some useful definitions and some lemmas.LetA(t) be continuous, cooperative, irreducible, and ω-periodic k×k matrix function, and ΨA(t) the fundamental matrix of system x′(t)=A(t)x(t). Denote by ρ(ΨA(ω)) the spectral radius of ΨA(ω).Lemma 3.
Letp=1/ωln⁡ρ(ΨA(ω)). Then there exists a positive, ω-periodic function v(t) such that eptv(t) is a solution of x′(t)=A(t)x(t) (see proof in Lemma  2.1 of [18]).Lemma 4.
Functionf(I) of model (1) satisfies f(I)≤f′(0)I, ∀I≥0.Proof.
Using assumptions on functionf, we have(28)ddIfII=If′I-fII2≤0,so function f(I)/I decreases ∀I>0 and then f(I)/I≤limI→0+fI/I=f′(0).Lemma 5.
Let(S(t),E(t),I(t)) be a solution of system (7) with initial conditions (S0,E0,I0)≥0, and (S^(t),0,0) the disease-free periodic solution of (7); then(29)limsupt→∞St-S^t≤0.Proof.
Proof is similar to Lemma4.1 of [14]. S(t) satisfies the first equation of system (7); then (30)dSdt=μN1-p-βtSfI-μ+rtS+δN-S-E-I≤Nμ1-p+δ-μ+rt+δS.Let X(t)=S(t)-S^(t); then (31)dXdt=μ+rt+δS^-S-βtSfI-δE+I≤-μ+rt+δX.Using Gronwall’s inequality X(t)≤X(0)e-∫0t(μ+r(s)+δ)ds, (32)St-S^t≤S0-S^0e-∫0tμ+rs+δds=S0-S^0e-μ+δte-∫0trsds.Taking limits in both sides, we obtain that limsupt→∞S(t)-S^(t)≤0.Now, we are able to enunciate our theorem for global stability of disease-free periodic solution.Theorem 6.
The disease-free periodic solution(S^(t),0,0) of system (7) is globally asymptotically stable if R0<1.Proof.
From Theorem2 we have that (S^(t),0,0) is unstable for R0>1 and asymptotically stable for R0<1, so it is sufficient to prove that any solution (S(t),E(t),I(t)) with nonnegative initial conditions (S0,E0,I0) approaches (S^,0,0) as t tends to infinity.
Letϵ>0; from Lemma 5 we have (33)limsupt→∞St-S^t=limt→∞supτ≥tSτ-S^τ=L≤0,so there exists a N>0 such that for all t1>N(34)-ϵ<supt≥t1St-S^t-L<ϵ,which implies that supt≥t1(S(t)-S^(t))<ϵ+L≤ϵ. Then, from the definition of supremum, we have that for all t>t1(35)St-S^t≤supt≥t1St-S^t<ϵ.
Then, we have proved that for allϵ>0 we can find a t1>0 such that S(t)<ϵ+S^(t) for all t>t1.
Now, using Lemma4, for ϵ>0 we can find a t1>0 such that for t>t1(36)dEdt=βtSfI-μ+σE≤βtStf′0I-μ+σEt(37)<βtf′0S^t+ϵIt-μ+σEt.We consider the following perturbed subsystem:(38)dE¯dt=βtf′0S^+ϵI¯-μ+σE¯,dI¯dt=σE¯-μ+γI¯,which can be rewritten as (39)dE¯dt,dI¯dtT=Ft-VtE¯,I¯T+ϵHtE¯,I¯T,with F(t),V(t) defined in (15) and(40)Ht=0βtf′000.
Matrix(F-V+ϵH)(t) is LT-periodic, cooperative, irreducible, and continuous. Using Lemma 3, if q=1/LTln⁡ρ(ΨF-V+ϵH(LT)), then there exists a positive and LT-periodic function v(t)=(v1(t),v2(t))T such that eqtv(t) is solution of system (38). Note that for all k>0, function keq(t-ti)v(t-ti) is also a solution of system (38) with initial condition kv(0) at t=ti.
Choose at¯>t1 and α1>0 such that (E(t¯),I(t¯))T≤α1v(0); then from (37) we have that (41)dEdt,dIdtT≤F-VE,IT+ϵHE,IT,and using a comparison principle (see, e.g., [19] Theorem B.1), we have (E(t),I(t))T≤α1eq(t-t¯)v(t-t¯) for all t>t¯.
From Theorem2.2 of [13], R0<1 iff ρ(ΦF-V(LT))<1. By the continuity of the spectrum for matrices (see [20], Section  II.5.8), we can choose ϵ>0 small enough so that ρ(ΦF-V+ϵH(LT))<1 and then q<0 (see Note 1). Thus, using positivity of solutions and comparison, (42)0≤limt→∞Et≤limt→∞α1eqt-t¯v1t-t¯=0.And similarly for I, we obtain that(43)limt→∞Et=0limt→∞It=0.
We need only to prove thatS(t) approaches S^. At disease-free periodic solution R^(t)=N-S^(t), where R^ satisfies (44)dR^dt=μNp+rtS^-μ+δR^.Thus, R(t)=N-S(t)-E(t)-I(t) satisfies(45)dR-R^dt=rtS-S^+γI-μ+δR-R^.
Letϵ1>0 be arbitrary and rmax=maxu∈[0,LT]r(u). Due to (43) we can find a t2>0 such that I(t)<ϵ1 for t>t2; moreover, we can find a t3>0 such that S(t)≤S^(t)+ϵ1 for t>t3. Then, let t4=max{t2,t3}; we have for t>t4(46)dR-R^dt≤rmax+γϵ1-μ+δR-R^.Multiplying in both sides by e(μ+δ)t and integrating from t4 to t, we obtain(47)R-R^≤R-R^t4e-μ+δt-t4+ϵ1rmax+γμ+δ1-e-μ+δt-t4.
So,limsupt→∞(R-R^)(t)≤ϵ1(rmax+γ)/μ+δ, where ϵ1(rmax+γ)/μ+δ is arbitrarily small. Then limsupt→∞(R-R^)(t)≤0, and using similar arguments for S and ϵ2>0, we can find a t5>0 with R(t)≤R^(t)+ϵ2/2 for t>t5. Also, from (43), we can find t6>0 with E(t)+I(t)<ϵ2/2 for t>t6, so, for t>max{t5,t6}, we have(48)St=N-Et-It-Rt≥N-R^t-ϵ2=S^t-ϵ2.
Or, equivalently,S(t)-S^(t)≥-ϵ2, with ϵ2 being arbitrarily small, and this implies that liminft→∞(S-S^)(t)≥0. We conclude by comparison and using Lemma 5 that limt→∞S(t)=S^(t), completing the proof.Theorem6 shows that disease will completely disappear as long as R0<1. Thus, reducing and keeping R0 below the unity would be sufficient to eradicate infection, even in a periodic environment and a general incidence rate.


### 3.2. Disease Persistence

Uniform persistence is an important concept in population dynamics, since it characterizes the long-term survival of some or all interacting species in an ecosystem [21].In this section we consider the dynamics of the periodic model whenR0>1. We will show that actually R0 is a threshold parameter for the extinction and the uniform persistence of the disease. Our results are inspired by [5, 15, 18, 22].LetP:X→X be the Poincaré map associated with system (7); that is, (49)Px0=ϕLT,x0,∀x0∈X,where X is defined in (8) and ϕ(t,x0) is the unique solution of system (7) with ϕ(0,x0)=x0. We define the following sets: (50)X0≔S,E,I∈X:E>0,I>0,∂X0≔X∖X0.Note that∂X0 is not the boundary of X0, but it is a standard notation of persistence theory.Lemma 7.
SetX0 is positively invariant under system (7).Proof.
Letx0=(S0,E0,I0)∈X0, that is, E0>0,I0>0, and let(51)ϕt,x0=St,Et,Itbe the solution of (7) with ϕ(0,x0)=x0. Due to nonnegativity of solutions and assumptions on function β(t) and f(I), we have (52)dEdt=βtSfI-μ+σE≥-μ+σE,∀t>0.Using a comparison theorem (see, e.g., [19] Appendix B.1), we have for all t>0(53)Et≥Ke-μ+σt>0,withK=E0>0.Similarly, (54)dIdt=σE-μ+γI≥-μ+γI,so, (55)It≥I0e-μ+γt>0,∀t>0.Therefore, ϕ(t,x0) remains on X0 for all t>0.To use persistence theory developed in [21], we show that(56)M∂=S,0,0:S≥0,where (57)M∂≔S0,E0,I0∈∂X0:PmS0,E0,I0∈∂X0,∀m≥0.Letx0=(S0,0,0)∈X and (S(t),E(t),I(t)) be the solution that passes through that initial condition. We have that ϕ(t,x0)=(S1(t),0,0), with S1(t) being a solution of (9) and S1(0)=S0 being a solution that satisfies the initial condition. By uniqueness of solutions we have E(t)=0=I(t)∀t≥0, so x0 lives on M∂.Now, ifx0∈M∂, we want x0=(S0,0,0). We prove an equivalent sentence: if x0∈∂X0∖{(S,0,0):S≥0}, then it does not belong to M∂. Consider an initial point x0=(S0,E0,I0)∈∂X0∖{(S,0,0):S≥0}; then E0>0,I0=0, or E0=0,I0>0. Suppose that E>0 and I0=0; then ϕ(t,x0) holds (58)dIdt0=σE0>0.By continuity ofE(t) and sign of derivative of I, we have that, for small 0<t≪1, E(t)>0,I(t)>0, so, for 0<t≪1, ϕ(t,x0)∈X0. Using invariance of X0 (Lemma 7) we have ϕ(t,x0)∈X0 for all t>1. Finally, for a m>0 such that mLT>1, we have Pm(x0)=ϕ(mLT,x0)∈X0 and this implies (56). By the existence of a disease-free periodic solution (proved in Section 2), it is clear that there is one fixed point of P in M∂ given by M0=(S^(0),0,0) ([23]).Now, we are in a position to introduce the following result of uniform persistence of the disease.Theorem 8.
LetR0>1; then there exists an ϵ>0 such that any solution (S(t),E(t)I(t)) of (7) with initial values (S(0),E(0),I(0))∈X0 satisfies(59)liminft→∞Et≥ϵ,liminft→∞It≥ϵ.Proof.
We first prove thatP is uniformly persistent (see Definition 1.3.2 from [21]) with respect to (X0,∂X0), because this implies that the solution of (7) is uniformly persistent with respect to (X0,∂X0) (see [21], Theorem 3.1.1). Clearly, X0 is relatively open in X, so ∂X0 is relatively closed.
Define(60)Ws≔x0∈X0:limm→∞Pmx0-M0=0;we show that Ws(M0)∩X0=∅.
By Theorem2.2 of [13], R0>1 if and only if r(ΨF-V(LT))>1. Choose an η>0 small enough with the property S^(t)-η>0,∀t>0 (see Appendix A). For α>0, let us consider the following perturbed equation:(61)dS¯dt=Nμ1-p+δ-2δα-βtf′0α+μ+rt+δS¯.System above admits a unique positive LT-periodic solution of the form(62)S^t,α=e-∫0tβsf′0α+μ+rs+δdsS^0,α+Nμ1-p+Nδ-2δα∫0te∫0sβξf′0α+μ+rξ+δdξdswhit S^(t,0)=S^(t), which is globally attractive for all solutions of (61) (see Appendix B for proof), and with(63)S^0,α=Nμ1-p+Nδ-2δα∫0LTe∫0sβξf′0α+μ+rξ+δdξdse∫0LTβsf′0α+μ+rs+δds-1.
SinceS^(0,α) is continuous in α, for all ϵ>0 there is a δ>0 such that for |α|<δ we have |S^(0,α)-S^(0,0)|<ϵ. Moreover, by continuity of solutions with respect to initial values we can find for all η¯>0 an ϵ¯>0 such that if |S^(0,α)-S^(0,0)|<ϵ¯, then (64)S^t,α-S^t,0<η¯.
Therefore, forη established before, we can find α small enough such that S^(t,α)>S^(t)-η,  ∀t>0.
Again, by continuity of solutions with respect to initial values, for this smallα>0, there exists a δ>0 such that for all (S0,E0,I0)∈X0 with (S0,E0,I0)-M0≤δ we have ϕ(t,(S0,E0,I0))-ϕ(t,M0)<α,  ∀t∈[0,LT].
We now claim that(65)limsupm→∞PmS0,E0,I0-M0≥δ,∀S0,E0,I0∈X0.By contradiction, suppose that(66)limsupm→∞PmS0,E0,I0-M0<δ,forsomeS0,E0,I0∈X0.
Without loss of generality, we can assume thatPm(S0,E0,I0)-M0<δ for all m≥0 (see Appendix C). From the discussion above, ϕ(t,Pm(S0,E0,I0))-ϕ(t,M0)<α, ∀m≥0 and t∈[0,LT].
For anyt≥0, let t=mLT+t1, where t1∈[0,LT) and m=[t/LT] is the greatest integer less than or equal to t/LT. Then, we get (67)ϕt,S0,E0,I0-ϕt,M0=ϕt1,PmS0,E0,I0-ϕt,M0<α.
If we setϕ(t,(S0,E0,I0))=(S(t),E(t),I(t)), then we have E(t)≤α, I(t)≤α, ∀t≥0, and from the first equation of (7) and Lemma 4 we arrive at(68)dSdt≥Nμ1-p+δ-2δα-βtf′0α+μ+rt+δS¯,which is exactly the equation in (61). Since the unique periodic solution of (61) is globally attractive, we have for S¯(t,α) solution of (61) that limt→∞S¯(t,α)=S^(t,α). So for η given before, there exists T>0 such that for all t≥T(69)S¯t,α-S^t,α<η,or equivalently S¯(t,α)>S^(t,α)-η. Moreover, from previous analysis, S^(t,α)-η>S^(t)-η; therefore, using comparison principle on (68) we arrive at(70)St≥S^t-ηfor t>T.
We haveE(t),I(t)≤α, and α is fixed small, so we can take α<ϵ∗ and use assumption (A5) in Introduction (see Appendix D) to obtain(71)dEdtdIdt≥F-V-ηH-αKE,IT,where F,V are defined in (15), H is defined in (40), and (72)K=0-12βtf′′0S^-η00.
By Theorem2.2 of [13], we have R0>1 iff ρ(ΦF-V(LT))>1. By continuity of spectrum (see [20] Section II), we can find α,ϵ such that (73)ρΦF-V-ηH-αK>1.Consider the auxiliary system (74)dE2dtdI2dt=F-V-ηH-αKE2,I2T;then, using Lemma 3 there exists a solution of (71) with the form ep2tv2(t), with p2=1/LTln⁡(ρ(ΦF-V-ηH-αK(LT)))>0. Choose a t2>T and a small number α2>0 such that (E2(t2),I2(t2))T≥α2v2(0). Using comparison principle we get (E(t),I(t))≥α2v2(t-t2)ep2(t-t2), which implies E(t)→∞ and I(t)→∞. This leads to a contradiction.
The claim above shows thatP is weakly uniformly persistent with respect to (X0,∂X0). Note that P has a global attractor S^(0) (see Lemma 5). It follows that M0 is an isolated invariant set in X, Ws(M0)∩X0=∅. Every orbit in M∂ converges to M0 and M0 is acyclic. By the acyclicity theorem on uniform persistence for maps ([21] Theorem 1.3.1 and Remark 1.3.1), it follows that P is uniformly persistent with respect to (X0,∂X0); that is, there exists ϵ>0 such that any solution of (7) satisfies limt→∞E(t)≥ϵ, limt→∞I(t)≥ϵ.


## 3.1. Disease Extinction

Theorem 2.
LetR0 be defined as (20); then the disease-free periodic solution (S^(t),0,0) is asymptotically stable if R0<1 and unstable if R0>1.Proof.
We use Theorem2.2 of [13] and check conditions (A1)–(A7). Conditions (A1)–(A5) are clearly satisfied from the definitions of F and V given in Section 2. We prove only conditions (A6) and (A7). Define (21)Mt≔-μ+rt+δ,and let ΦM(t) be the monodromy matrix of system(22)dzdt=Mtz.(A6)
ρ(ΦM(LT))<1. Let ΨM be a fundamental matrix for system dz/dt=M(t)z, with M defined as before and LT periodic; the monodromy matrix ΦM(LT) is given by ΦM(LT)=ΨM-1(0)ΨM(LT). The general solution of (22) is (23)zt=Kexp⁡-∫0tμ+rs+δds,soΨM=exp(-∫0t(μ+r(s)+δ)ds) and ΨM-1=exp(∫0t(μ+r(s)+δ)ds). Note that ΨM-1(0)=1, so ΦM(LT)=ΨM(LT) and(24)ΦMLT=exp-∫0LTμ+rs+δds.Due to the fact thatΦM(LT) is a constant, its eigenvalue is itself and ρ(ΦM(LT))<1 for μ,δ,r(s)>0.(A7)
ρ(Φ-V(LT))<1. Solving the system dz/dt=-V(t)z, we arrive at the general solution (25)zt=c1γ-σσ1e-μ+σt+c201e-μ+γt,so(26)Ψ-Vt=γ-σσe-μ+σt0e-μ+σte-μ+γt.ComputingΦ-V(LT)=Ψ-V-1(0)Ψ-V(LT), we have(27)Φ-VLT=e-μ+σLT00e-μ+γLT.Clearly,ρ(Φ-V(LT))=max{e-(μ+σ)LT,e-(μ+γ)LT}<1 for μ,γ,σ>0.Note 1.
Due to the fact thatΨA is a fundamental solution of a periodic system, we can always choose it such that Ψ(0)=I, so the monodromy matrix satisfies ΦA(LT)=ΨA(LT). This property is used in further analysis.In order to prove the global stability of the disease-free periodic solution, we enunciate some useful definitions and some lemmas.LetA(t) be continuous, cooperative, irreducible, and ω-periodic k×k matrix function, and ΨA(t) the fundamental matrix of system x′(t)=A(t)x(t). Denote by ρ(ΨA(ω)) the spectral radius of ΨA(ω).Lemma 3.
Letp=1/ωln⁡ρ(ΨA(ω)). Then there exists a positive, ω-periodic function v(t) such that eptv(t) is a solution of x′(t)=A(t)x(t) (see proof in Lemma  2.1 of [18]).Lemma 4.
Functionf(I) of model (1) satisfies f(I)≤f′(0)I, ∀I≥0.Proof.
Using assumptions on functionf, we have(28)ddIfII=If′I-fII2≤0,so function f(I)/I decreases ∀I>0 and then f(I)/I≤limI→0+fI/I=f′(0).Lemma 5.
Let(S(t),E(t),I(t)) be a solution of system (7) with initial conditions (S0,E0,I0)≥0, and (S^(t),0,0) the disease-free periodic solution of (7); then(29)limsupt→∞St-S^t≤0.Proof.
Proof is similar to Lemma4.1 of [14]. S(t) satisfies the first equation of system (7); then (30)dSdt=μN1-p-βtSfI-μ+rtS+δN-S-E-I≤Nμ1-p+δ-μ+rt+δS.Let X(t)=S(t)-S^(t); then (31)dXdt=μ+rt+δS^-S-βtSfI-δE+I≤-μ+rt+δX.Using Gronwall’s inequality X(t)≤X(0)e-∫0t(μ+r(s)+δ)ds, (32)St-S^t≤S0-S^0e-∫0tμ+rs+δds=S0-S^0e-μ+δte-∫0trsds.Taking limits in both sides, we obtain that limsupt→∞S(t)-S^(t)≤0.Now, we are able to enunciate our theorem for global stability of disease-free periodic solution.Theorem 6.
The disease-free periodic solution(S^(t),0,0) of system (7) is globally asymptotically stable if R0<1.Proof.
From Theorem2 we have that (S^(t),0,0) is unstable for R0>1 and asymptotically stable for R0<1, so it is sufficient to prove that any solution (S(t),E(t),I(t)) with nonnegative initial conditions (S0,E0,I0) approaches (S^,0,0) as t tends to infinity.
Letϵ>0; from Lemma 5 we have (33)limsupt→∞St-S^t=limt→∞supτ≥tSτ-S^τ=L≤0,so there exists a N>0 such that for all t1>N(34)-ϵ<supt≥t1St-S^t-L<ϵ,which implies that supt≥t1(S(t)-S^(t))<ϵ+L≤ϵ. Then, from the definition of supremum, we have that for all t>t1(35)St-S^t≤supt≥t1St-S^t<ϵ.
Then, we have proved that for allϵ>0 we can find a t1>0 such that S(t)<ϵ+S^(t) for all t>t1.
Now, using Lemma4, for ϵ>0 we can find a t1>0 such that for t>t1(36)dEdt=βtSfI-μ+σE≤βtStf′0I-μ+σEt(37)<βtf′0S^t+ϵIt-μ+σEt.We consider the following perturbed subsystem:(38)dE¯dt=βtf′0S^+ϵI¯-μ+σE¯,dI¯dt=σE¯-μ+γI¯,which can be rewritten as (39)dE¯dt,dI¯dtT=Ft-VtE¯,I¯T+ϵHtE¯,I¯T,with F(t),V(t) defined in (15) and(40)Ht=0βtf′000.
Matrix(F-V+ϵH)(t) is LT-periodic, cooperative, irreducible, and continuous. Using Lemma 3, if q=1/LTln⁡ρ(ΨF-V+ϵH(LT)), then there exists a positive and LT-periodic function v(t)=(v1(t),v2(t))T such that eqtv(t) is solution of system (38). Note that for all k>0, function keq(t-ti)v(t-ti) is also a solution of system (38) with initial condition kv(0) at t=ti.
Choose at¯>t1 and α1>0 such that (E(t¯),I(t¯))T≤α1v(0); then from (37) we have that (41)dEdt,dIdtT≤F-VE,IT+ϵHE,IT,and using a comparison principle (see, e.g., [19] Theorem B.1), we have (E(t),I(t))T≤α1eq(t-t¯)v(t-t¯) for all t>t¯.
From Theorem2.2 of [13], R0<1 iff ρ(ΦF-V(LT))<1. By the continuity of the spectrum for matrices (see [20], Section  II.5.8), we can choose ϵ>0 small enough so that ρ(ΦF-V+ϵH(LT))<1 and then q<0 (see Note 1). Thus, using positivity of solutions and comparison, (42)0≤limt→∞Et≤limt→∞α1eqt-t¯v1t-t¯=0.And similarly for I, we obtain that(43)limt→∞Et=0limt→∞It=0.
We need only to prove thatS(t) approaches S^. At disease-free periodic solution R^(t)=N-S^(t), where R^ satisfies (44)dR^dt=μNp+rtS^-μ+δR^.Thus, R(t)=N-S(t)-E(t)-I(t) satisfies(45)dR-R^dt=rtS-S^+γI-μ+δR-R^.
Letϵ1>0 be arbitrary and rmax=maxu∈[0,LT]r(u). Due to (43) we can find a t2>0 such that I(t)<ϵ1 for t>t2; moreover, we can find a t3>0 such that S(t)≤S^(t)+ϵ1 for t>t3. Then, let t4=max{t2,t3}; we have for t>t4(46)dR-R^dt≤rmax+γϵ1-μ+δR-R^.Multiplying in both sides by e(μ+δ)t and integrating from t4 to t, we obtain(47)R-R^≤R-R^t4e-μ+δt-t4+ϵ1rmax+γμ+δ1-e-μ+δt-t4.
So,limsupt→∞(R-R^)(t)≤ϵ1(rmax+γ)/μ+δ, where ϵ1(rmax+γ)/μ+δ is arbitrarily small. Then limsupt→∞(R-R^)(t)≤0, and using similar arguments for S and ϵ2>0, we can find a t5>0 with R(t)≤R^(t)+ϵ2/2 for t>t5. Also, from (43), we can find t6>0 with E(t)+I(t)<ϵ2/2 for t>t6, so, for t>max{t5,t6}, we have(48)St=N-Et-It-Rt≥N-R^t-ϵ2=S^t-ϵ2.
Or, equivalently,S(t)-S^(t)≥-ϵ2, with ϵ2 being arbitrarily small, and this implies that liminft→∞(S-S^)(t)≥0. We conclude by comparison and using Lemma 5 that limt→∞S(t)=S^(t), completing the proof.Theorem6 shows that disease will completely disappear as long as R0<1. Thus, reducing and keeping R0 below the unity would be sufficient to eradicate infection, even in a periodic environment and a general incidence rate.


## 3.2. Disease Persistence

Uniform persistence is an important concept in population dynamics, since it characterizes the long-term survival of some or all interacting species in an ecosystem [21].In this section we consider the dynamics of the periodic model whenR0>1. We will show that actually R0 is a threshold parameter for the extinction and the uniform persistence of the disease. Our results are inspired by [5, 15, 18, 22].LetP:X→X be the Poincaré map associated with system (7); that is, (49)Px0=ϕLT,x0,∀x0∈X,where X is defined in (8) and ϕ(t,x0) is the unique solution of system (7) with ϕ(0,x0)=x0. We define the following sets: (50)X0≔S,E,I∈X:E>0,I>0,∂X0≔X∖X0.Note that∂X0 is not the boundary of X0, but it is a standard notation of persistence theory.Lemma 7.
SetX0 is positively invariant under system (7).Proof.
Letx0=(S0,E0,I0)∈X0, that is, E0>0,I0>0, and let(51)ϕt,x0=St,Et,Itbe the solution of (7) with ϕ(0,x0)=x0. Due to nonnegativity of solutions and assumptions on function β(t) and f(I), we have (52)dEdt=βtSfI-μ+σE≥-μ+σE,∀t>0.Using a comparison theorem (see, e.g., [19] Appendix B.1), we have for all t>0(53)Et≥Ke-μ+σt>0,withK=E0>0.Similarly, (54)dIdt=σE-μ+γI≥-μ+γI,so, (55)It≥I0e-μ+γt>0,∀t>0.Therefore, ϕ(t,x0) remains on X0 for all t>0.To use persistence theory developed in [21], we show that(56)M∂=S,0,0:S≥0,where (57)M∂≔S0,E0,I0∈∂X0:PmS0,E0,I0∈∂X0,∀m≥0.Letx0=(S0,0,0)∈X and (S(t),E(t),I(t)) be the solution that passes through that initial condition. We have that ϕ(t,x0)=(S1(t),0,0), with S1(t) being a solution of (9) and S1(0)=S0 being a solution that satisfies the initial condition. By uniqueness of solutions we have E(t)=0=I(t)∀t≥0, so x0 lives on M∂.Now, ifx0∈M∂, we want x0=(S0,0,0). We prove an equivalent sentence: if x0∈∂X0∖{(S,0,0):S≥0}, then it does not belong to M∂. Consider an initial point x0=(S0,E0,I0)∈∂X0∖{(S,0,0):S≥0}; then E0>0,I0=0, or E0=0,I0>0. Suppose that E>0 and I0=0; then ϕ(t,x0) holds (58)dIdt0=σE0>0.By continuity ofE(t) and sign of derivative of I, we have that, for small 0<t≪1, E(t)>0,I(t)>0, so, for 0<t≪1, ϕ(t,x0)∈X0. Using invariance of X0 (Lemma 7) we have ϕ(t,x0)∈X0 for all t>1. Finally, for a m>0 such that mLT>1, we have Pm(x0)=ϕ(mLT,x0)∈X0 and this implies (56). By the existence of a disease-free periodic solution (proved in Section 2), it is clear that there is one fixed point of P in M∂ given by M0=(S^(0),0,0) ([23]).Now, we are in a position to introduce the following result of uniform persistence of the disease.Theorem 8.
LetR0>1; then there exists an ϵ>0 such that any solution (S(t),E(t)I(t)) of (7) with initial values (S(0),E(0),I(0))∈X0 satisfies(59)liminft→∞Et≥ϵ,liminft→∞It≥ϵ.Proof.
We first prove thatP is uniformly persistent (see Definition 1.3.2 from [21]) with respect to (X0,∂X0), because this implies that the solution of (7) is uniformly persistent with respect to (X0,∂X0) (see [21], Theorem 3.1.1). Clearly, X0 is relatively open in X, so ∂X0 is relatively closed.
Define(60)Ws≔x0∈X0:limm→∞Pmx0-M0=0;we show that Ws(M0)∩X0=∅.
By Theorem2.2 of [13], R0>1 if and only if r(ΨF-V(LT))>1. Choose an η>0 small enough with the property S^(t)-η>0,∀t>0 (see Appendix A). For α>0, let us consider the following perturbed equation:(61)dS¯dt=Nμ1-p+δ-2δα-βtf′0α+μ+rt+δS¯.System above admits a unique positive LT-periodic solution of the form(62)S^t,α=e-∫0tβsf′0α+μ+rs+δdsS^0,α+Nμ1-p+Nδ-2δα∫0te∫0sβξf′0α+μ+rξ+δdξdswhit S^(t,0)=S^(t), which is globally attractive for all solutions of (61) (see Appendix B for proof), and with(63)S^0,α=Nμ1-p+Nδ-2δα∫0LTe∫0sβξf′0α+μ+rξ+δdξdse∫0LTβsf′0α+μ+rs+δds-1.
SinceS^(0,α) is continuous in α, for all ϵ>0 there is a δ>0 such that for |α|<δ we have |S^(0,α)-S^(0,0)|<ϵ. Moreover, by continuity of solutions with respect to initial values we can find for all η¯>0 an ϵ¯>0 such that if |S^(0,α)-S^(0,0)|<ϵ¯, then (64)S^t,α-S^t,0<η¯.
Therefore, forη established before, we can find α small enough such that S^(t,α)>S^(t)-η,  ∀t>0.
Again, by continuity of solutions with respect to initial values, for this smallα>0, there exists a δ>0 such that for all (S0,E0,I0)∈X0 with (S0,E0,I0)-M0≤δ we have ϕ(t,(S0,E0,I0))-ϕ(t,M0)<α,  ∀t∈[0,LT].
We now claim that(65)limsupm→∞PmS0,E0,I0-M0≥δ,∀S0,E0,I0∈X0.By contradiction, suppose that(66)limsupm→∞PmS0,E0,I0-M0<δ,forsomeS0,E0,I0∈X0.
Without loss of generality, we can assume thatPm(S0,E0,I0)-M0<δ for all m≥0 (see Appendix C). From the discussion above, ϕ(t,Pm(S0,E0,I0))-ϕ(t,M0)<α, ∀m≥0 and t∈[0,LT].
For anyt≥0, let t=mLT+t1, where t1∈[0,LT) and m=[t/LT] is the greatest integer less than or equal to t/LT. Then, we get (67)ϕt,S0,E0,I0-ϕt,M0=ϕt1,PmS0,E0,I0-ϕt,M0<α.
If we setϕ(t,(S0,E0,I0))=(S(t),E(t),I(t)), then we have E(t)≤α, I(t)≤α, ∀t≥0, and from the first equation of (7) and Lemma 4 we arrive at(68)dSdt≥Nμ1-p+δ-2δα-βtf′0α+μ+rt+δS¯,which is exactly the equation in (61). Since the unique periodic solution of (61) is globally attractive, we have for S¯(t,α) solution of (61) that limt→∞S¯(t,α)=S^(t,α). So for η given before, there exists T>0 such that for all t≥T(69)S¯t,α-S^t,α<η,or equivalently S¯(t,α)>S^(t,α)-η. Moreover, from previous analysis, S^(t,α)-η>S^(t)-η; therefore, using comparison principle on (68) we arrive at(70)St≥S^t-ηfor t>T.
We haveE(t),I(t)≤α, and α is fixed small, so we can take α<ϵ∗ and use assumption (A5) in Introduction (see Appendix D) to obtain(71)dEdtdIdt≥F-V-ηH-αKE,IT,where F,V are defined in (15), H is defined in (40), and (72)K=0-12βtf′′0S^-η00.
By Theorem2.2 of [13], we have R0>1 iff ρ(ΦF-V(LT))>1. By continuity of spectrum (see [20] Section II), we can find α,ϵ such that (73)ρΦF-V-ηH-αK>1.Consider the auxiliary system (74)dE2dtdI2dt=F-V-ηH-αKE2,I2T;then, using Lemma 3 there exists a solution of (71) with the form ep2tv2(t), with p2=1/LTln⁡(ρ(ΦF-V-ηH-αK(LT)))>0. Choose a t2>T and a small number α2>0 such that (E2(t2),I2(t2))T≥α2v2(0). Using comparison principle we get (E(t),I(t))≥α2v2(t-t2)ep2(t-t2), which implies E(t)→∞ and I(t)→∞. This leads to a contradiction.
The claim above shows thatP is weakly uniformly persistent with respect to (X0,∂X0). Note that P has a global attractor S^(0) (see Lemma 5). It follows that M0 is an isolated invariant set in X, Ws(M0)∩X0=∅. Every orbit in M∂ converges to M0 and M0 is acyclic. By the acyclicity theorem on uniform persistence for maps ([21] Theorem 1.3.1 and Remark 1.3.1), it follows that P is uniformly persistent with respect to (X0,∂X0); that is, there exists ϵ>0 such that any solution of (7) satisfies limt→∞E(t)≥ϵ, limt→∞I(t)≥ϵ.


## 4. Numerical Simulations

In this section we provide some numerical simulations to illustrate the results obtained in our theorems and compare them with previous results.To improve previous models used in references, we use a particular function(75)fI=I1+aI,a≥0,which includes the case f(I)=I used in [5]. One can check that function (75) satisfies conditions (A1)–(A5). Using this function, system (7) is rewritten as(76)dSdt=μN1-p-βtSI1+aI-μ+rtS+δN-S-E-I,dEdt=βtSI1+aI-μ+σE,dIdt=σE-μ+γI.Set an initial populationN=2,200,000 and take time t in years. Suppose μ=0.02 per year, corresponding to an average human life time of 50 years. Following [5] take the parameters as follows: σ=38.5 per year, γ=100 per year, p=0.85,δ=0, and a=1. Choose the periodic transmission as β(t)=β0+0.0002cos⁡(2πt), with β0 being the transmission parameter, and the periodic vaccination rate r(t)=0.1+0.004cos⁡(2πt). Both functions have period LT=1.There exists multiple methods for computing the basic reproduction number, via numerical approximations, or finding a positive solution of the equationρ(W(LT,0,λ))=1 (see Theorem 2.1 of [13]). In order to compare our work with previous works, we approximate the basic reproduction number with its average value R0T, used by several authors as the reproduction number (for example [10, 11]), so define(77)R0T=ρFV-1,where V is given by (15) and (78)F=0βS^f′000,with β,[S^] being the average of functions β,S^ defined as [β]=1/LT∫0LTβ(t)dt,S^=1/LT∫0LTS^(t)dt. Computing each average, we obtain (79)R0T=549.6702634β0,so R0T>1 for β0∈(0.001819272510,∞).Following Theorem2.1 of [13], to compute R0, let W(t,s,λ),t≥s, be the evolution operator of the system(80)dwdt=-Vt+Ftλwthat is, for each λ∈(0,∞), dW(t,s,λ)/dt=-V(t)+F(t)/λW(t,s,λ), ∀t≥s, and W(s,s,λ)=I2×2. With this operator, R0>0 is the unique solution of ρ(W(LT,0,λ))=1.Example 1.
To illustrate our results, fixβ0=0.0018. Computing R0T, we have R0T=0.9894064741, which is a first approximation of R0. To solve system (80) numerically, we substitute the terms of expression of S^(t) in (10): (81)S^t=e-0.1200000000t-0.0006366197724sin⁡6.283185307t54999.33689+6600.0∫0.0te0.1200000000s+0.0006366197724sin⁡6.283185307sds
The previous integral cannot be computed analytically, so we approachS^(t) using Taylor expansion around 0 (remember that we want so solve ρ(W(LT,0,λ))=1, where LT=1), so even when we cannot find an explicit expression for S^(t), the Taylor expansion is a good way to estimate it in (0,1). It could be of interest to also use an approach of S^(t) around t=1 and compare the results with those obtained in the present work (see Section 5 for a discussion about this topic).
Setting an initial valueλ0=0.98 and letting λi=λ0+i(0.0001), we solve system (80) numerically for each λi (using initial conditions w(0)=(1,0) and w(0)=(0,1), to satisfy W(0,0)=I2×2), and compute ρ1=ρ(W(LT,0,λi)) until ρ1~1. With previous process we arrive at ρ1=1.00120166209265 for λ=0.9872 and ρ1=0.997826338969630 for λ=0.9873; therefore R0∈(0.9872,0.9873). Using a finer step size 0.0000001 to have more accuracy, we arrive at R0~0.9872355<1.
Set initial values asS(0)=1,500,000, E(0)=400,000, I(0)=40,000, and R(0)=N-(S(0)+E(0)+I(0)).
There exist multiple numerical methods to compute and plot the solutions of nonautonomous differential equations; see, for example, the Adomian method, the homotopy analysis method, or the modified homotopy methods (see, e.g., [24, 25]). For this work we use Matlab algorithms (ODE 45) to graph the solution of system (76) with these initial conditions. Figures 2 and 3 shows the results. We can see that I(t),E(t) goes to zero, while S(t),R(t) tend to stabilize; also S(t) is tending to S^(t) with values between 54,000 and 56,000 (see Figure 1); this shows the results obtained in Theorem 6.Figure 1
S component of infection-free periodic solution. Time t is given in years. (a) S^(t), because it is a periodic function of period 1, is plotted only in [0,1]. (b) Taylor expansion of S^ around t=0 of order t10.
(a)
(b)Figure 2
Solution of exposed and infected populations of SEIRS system whenR0<1. We can see that both approach zero when time goes to infinity. Time t is given in years.Figure 3
Solution of susceptible and recovered populations of SEIRS system whenR0<1. We can see that S approaches S^(t) (see Figure 1), and R approaches R^(t)=N-S^(t). Time t is given in years.Example 2.
Now, chooseβ0=0.005. As we can see in Figures 4 and 5, the solutions of system (1) remain persistent when t tends to infinity; this fact suggests that R0>1 from Theorem 8. In fact, if we compute the basic reproduction number and its average (using the process described in example 1), R0T=3.298021580 and R0∈(2.7456,2.7457); therefore it is bigger than one. In fact, this shows the results of persistence obtained in Theorem 8.Figure 4
Solution of exposed and infected individuals of SEIRS system whenR0>1. Both E and I remain persistent when time goes to infinity. Time t is given in years.Figure 5
Solution of susceptible and recovered populations of SEIRS system whenR0>1. Time t is given in years.


## 5. Conclusion

In this paper we presented a model with seasonal fluctuation with a general incidence functionSf(I) that includes the bilinear case βSI (studied by [5]) and a family of saturated incidence rate of the form βSI/(1+kIq). We proved the existence of a disease-free periodic solution (S^(t),0,0) and defined the basic reproduction number R0, proving that it is a threshold parameter for disease, in the sense that when R0<1, the disease-free periodic solution is globally asymptotically stable, and when R0>1, the disease is persistent. A next step of this work is to consider a family of incidence rates more generally, changing Sf(I) by f(S,I) and trying to obtain results of persistence and stability similar to the ones obtained in this work. Another interesting topic is to ask what the behavior of system at R0=1 is, in order to complete the analysis that we have made.Several authors (e.g., [10, 11]) define R0 as an average, which we denoted as R0T to distinguish between it and the basic reproduction number defined by [13], via the monodromy matrix (which is a real threshold parameter for extinction and persistence of disease). We compute R0T, approximate R0 (with the help of Taylor theorem), and compare these values, obtaining that R0T is not equal to R0; moreover R0T>R0 in both examples (similar comparisons can be observed also in the works made by [13, 17]). This fact suggests that the use of R0T for persistence overestimates the threshold. To emphasize this conclusion, it would be helpful to find an example where R0<1 but R0T>1 and then compute the solutions to observe the behavior (we affirm that the disease will go extinct due to Theorem 6).To obtain the estimation ofR0 we used a code in Maple, which is based on numerical computing of ρ1=ρ(W(LT,0,λi)) until ρ1~1, where λi=λ0+Δλi, Δλ is the step size, and the initial estimation λ0 is taken as R0T-ϵ. For this approximation we have used a Taylor expansion of the periodic solution S^(t); another interesting possibility could be varying the approximation used for R0, for example, changing the Taylor approach of S^(t) around t=1 instead of t=0. The graphs of the solutions were obtained with ODE 45 from Matlab, but other methods can be used to improve them, for example, Adomian methods or homotopy methods [24, 25]. The Maple code used to estimate R0 is available for anyone who wants to use it.



---

*Source: 5796958-2017-11-09.xml*