Global dynamics of a periodic SEIRS model with general incidence rate

TL;DR

This paper analyzes a family of periodic SEIRS epidemic models with general incidence rates, demonstrating that the basic reproduction number governs the global dynamics and persistence of the disease, supported by numerical simulations.

Contribution

It introduces a general periodic SEIRS model with nonlinear incidence rates and establishes the basic reproduction number as a threshold for disease persistence and global behavior.

Findings

Basic reproduction number determines global dynamics.

Numerical simulations validate analytical results.

Model accommodates general nonlinear incidence rates.

Abstract

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

Equations

$$\frac{dS}{dt}$$

$$\frac{dE}{dt}$$

$$\frac{dI}{dt}$$

$$\frac{dR}{dt}$$

$$(S(0),E(0),I(0),R(0))=(S_0,E_0,I_0,R_0)$$

$$\frac{dS}{dt}(t=t_1)=\mu N(1-p)+\delta R\ge 0,$$

$$\frac{dE}{dt}(t=t_1)=\beta (t)Sf(I)\ge 0.$$

$$\frac{dI}{dt}(t=t_1)=\sigma E\ge 0.$$

$$\frac{dR}{dt}(t=t_1)=\mu Np+r(t)S+\gamma I\ge 0.$$

$$\frac{dS}{dt}$$

$$\frac{dE}{dt}$$

$$\frac{dI}{dt}$$

$$X=\{(S,E,I)\in {\mathbb{R}}_{+}^3:S+E+I\le N\}.$$

$$\frac{dS}{dt}=\mu N(1-p)-(\mu +r(t))S+\delta (N-S),S(0)=S^0\in {\mathbb{R}}_{+}.$$

$$\hat{S}(t)=e^{-{\int }_0^t(\mu +r(s)+\delta )ds}\left(\hat{S}(0)+N(\mu (1-p)+\delta ){\int }_0^t{e}^{{\int }_0^s(\mu +r(\xi )+\delta )d\xi }ds\right),$$

$$\hat{S}(0)=\frac{N(\mu (1-p)+\delta ){\int }_0^{LT}{e}^{{\int }_0^s(\mu +r(\xi )+\delta )d\xi }ds}{e^{{\int }_0^{LT}(\mu +r(s)+\delta )ds}-1}.$$

$$x^{\prime }(t)$$

$$=\mathcal{F}-\mathcal{V},$$

$$\mathcal{F}$$

$${\mathcal{V}}^{-}$$

$$D\mathcal{F}(0,0,\hat{S})$$

$$D\mathcal{V}(0,0,\hat{S})$$

$$F(t)=\left(\begin{array}{cc}0& \beta (t)\hat{S}{f}^{\prime }(0)\\ 0& 0\end{array}\right),V(t)=\left(\begin{array}{cc}\mu +\sigma & 0\\ -\sigma & \mu +\gamma \end{array}\right).$$

$$(\mathcal{L}\varphi )(t)={\int }_0^{\infty }Y(t,t-a)F(t-a)\varphi (t-a)da,\forall t\in \mathbb{R},\varphi \in C_{LT},$$

$$\frac{dy}{dt}=-V(t)y,$$

$$\frac{dY(t,s)}{dt}=-V(t)Y(t,s),\forall t\ge s,Y(s,s)=I_{2\times 2}.$$

$$\mathcal{L}v=r_{0}v.$$

$${\mathcal{R}}_0:=\rho (\mathcal{L}),$$

$$M(t):=-(\mu +r(t)+\delta ),$$

$$\frac{dz}{dt}=M(t)z.$$

$$z(t)=K\exp (-{\int }_0^t(\mu +r(s)+\delta )ds),$$

$${\Phi }_M(LT)=\exp (-{\int }_0^{LT}(\mu +r(s)+\delta )ds).$$

$$z(t)=c_1\left(\begin{array}{c}\frac{\gamma -\sigma }{\sigma }\\ 1\end{array}\right)e^{-(\mu +\sigma )t}+c_2\left(\begin{array}{c}0\\ 1\end{array}\right)e^{-(\mu +\gamma )t},$$

$${\Psi }_{-V}(t)=\left(\begin{array}{cc}\frac{\gamma -\sigma }{\sigma }{e}^{-(\mu +\sigma )t}& 0\\ e^{-(\mu +\sigma )t}& e^{-(\mu +\gamma )t}\end{array}\right).$$

Full text

1 Introduction

Epidemiological models have been recognized as valuable tools in analyzing the spread and control of infectious diseases. In the study of epidemiological models, 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. Incidence rate has been developed by many authors. In order to model this disease transmission process several authors employ the incidence functions: The earliest one is the bilinear incidence rate $\beta SI$ used by Kermack and Mckendrick [8] in 1927, where $\beta,S$ and $I$ denote the transmission rate, the number of susceptible population and the infectious population respectively. It is based on the law of mass action which is not realistic. So there is a need to modify the classical linear incidence rate to study the dynamics of infection among large population. In 1978, Capasso and Serio [6] introduced a saturated incidence rate by research of the Cholera epidemic spread in Bari. Also in 1978, May and Anderson [1] proposed the saturated incidence rate.

In the present work we focus on SEIRS epidemic models. We improve the model of Moneim and Greenhalgh in [11], introducing an incidence rate with a general function taken from [4] and the references therein.

We propose the following SEIRS model:

[TABLE]

Where $N=S+E+I+R$ is the total population size, with $S,E,I,R$ denoting the fractions of population that are suceptible, exposed, infected and revovered, respectively. $\beta(t)$ is the transmission rate and it is a continuous, positive T -periodic function. $p$ ($0\leq p\leq 1$) is the vaccination rate of all new-born 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. $\mu$ is the common per capita birth and death rate. $\sigma,\gamma$ and $\delta$ are the per capita rates of leaving the latent stage, infected stage and recovered stage, respectively. It is assumed that parameters are positive constants.

Bai and Zhou in [5] answered some open problems stated in [11] , they also shown that their condition is a treshold between persistence and extinction of the disease via the framework established in [16]. They assumed that the incidence was bilinear. In our study, the nonlinear assumptions on function $f$ are listed below (see [4] ):

• A1)

$f:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is continuously differentiable.

• A2)

$f(0)=0,f^{\prime}(0)>0$ and $f(I)>0$ for all $I>0$.

• A3)

$f(I)-If^{\prime}(I)\geq 0$.

Under these assuptions, function $f(I)$ includes various types of incidence rate, particularly, when $f(I)=I$, we are on the bilinear case considered by Moneim.

In addition, we assume following extra conditions (see [13]):

• A4)

$f^{\prime\prime}(0)\leq 0$.

• A5)

There exists $\epsilon^{*}>0$ such that when $0<I<\epsilon^{*}$, $f(I)\geq f(0)+If^{\prime}(0)+\frac{1}{2}I^{2}f^{\prime\prime}(0).$

This set of assumptions on the function f allows for more general incidence functions than the bilinear one, like saturated incidence functions and functions of the form $\beta SI/(1+kI^{q})$, 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 non monotone 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^{\prime}(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.

The layout of this paper is as follows.In section 2,we introduce the basic reproduction number via the theory developed in [3], [16]. In section 3, we adapt the arguments given in [5] to prove that the disease free periodic solution is globally asymptotically stable if $\mathcal{R}_{0}<1$ and that if $\mathcal{R}_{0}>1$ system (1) is persistent.

We consider a family of $SEIR$ models with periodic coefficients with general incidence rate in epidemiology. We show that the global dynamics is determined by the basic reproduction number $\mathcal{R}_{0}$ Our results generalize the ones in [5].

2 The basic reproduction number

First of all, we prove non-negativity of solutions under non-negative initial conditions.

Theorem 1**.**

Let $S_{0},E_{0},I_{0},R_{0}\geq 0$, the solution of (1) with

[TABLE]

is non negative in sense that $S(t),E(t),I(t),R(t)\geq 0$ , $\forall t>0$, and satisfies $S(t)+E(t)+I(t)+R(t)=N$, with $N$ constant.

Proof.

Let $x(t)=(S(t),E(t),I(t),R(t))$ be the solution of system (2) under initial conditions $x_{0}=(S(0),E(0),I(0),R(0))=(S_{0},E_{0},I_{0},R_{0})\geq 0$, by continuity of solution, 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,t_{0})$ such that $S(t),E(t),I(t),R(t)\geq 0$ for $0<t<t_{0}$. We will prove that $t_{0}=\infty$.

If $S(t_{1})=0$ for a $t_{1}\geq 0$ and other components remain non-negative at $t=t_{1}$, then

[TABLE]

this implies that whenever the solution $x(t)$ touches the $S$-axis, the derivative of $S$ is non decreasing and the function $S(t)$ does not cross to negative values. Similarly: When $E(t_{1})=0$ for a $t_{1}>0$ and other components remain non-negative:

[TABLE]

When $I(t_{1})=0$ for a $t_{1}>0$ and other components remain non-negative:

[TABLE]

Finally, when $R(t_{1})=0$ for a $t_{1}>0$ and other components remain non-negative:

[TABLE]

Therefore, whenever $x(t)$ touches any of the axis $S=0,E=0,I=0,R=0$, it never crosses them.

Now, let $N(t)=S(t)+E(t)+I(t)+R(t)$, then adding all equations of system we can see that $\frac{dN}{dt}=0$, so the value of $N$ is constant. ∎

To reduce the system (1), let With $R=N-S-E-I$, system (1) is reduced to

[TABLE]

The dynamics of system (1) is equivalent to that of (2). Moreover, due to positivity of solutions $S+E+I\leq N$, so we study the dynamic of system (2) in the region

[TABLE]

A disease free periodic solution can be found for (2). To find it, set $E=0=I$, then from first equation of (2):

[TABLE]

From [5] and [11], the equation above admits a unique positive LT-periodic solution given by:

[TABLE]

where

[TABLE]

Then, $(\hat{S}(t),0,0)$ is a periodic infection free solution of (2), moreover, from [5] we have that $\hat{S}(t)<N$, therefore, $(\hat{S}(t),0,0)$ lives in $X$.

Using the notation of [15], we sort the compartments so that the first 2 compartments correspond to infected individuals. Let $x=(E,I,S)$ and define

• •

$\mathcal{F}_{i}$: the rate of new infection in compartment i.

• •

$\mathcal{V}_{i}^{+}:$ the rate of individuals into compartment i by other means.

• •

$\mathcal{V}_{i}^{-}:$ the rate of transfer individuals out of compartment i.

System can be written as

[TABLE]

where $\mathcal{V}=\mathcal{V}^{-}-\mathcal{V}^{+}$,

[TABLE]

Linearizing system (6) around the disease free solution, we obtain the matrix of partial derivatives $J(0,0,\hat{S})=D\mathcal{F}(0,0,\hat{S})-D\mathcal{V}(0,0,\hat{S})$, where

[TABLE]

Using lemma 1 of [15], we part $D\mathcal{F}$ and $D\mathcal{V}$ and set

[TABLE]

For a compartmental epidemiological model based on an autonomous system, the basic reproduction number is determined by the spectral radius of the next-generation matrix $FV^{-1}$(which is independent of time) [15]. The definition of basic reproduction number for non autonomous systems has been studied for multiple authors, see for example [3] and [16]. Particularly, Wang and Zhao in [16] extended the work of [15] to include epidemiological models in periodic environments. They introduced the next infection operator $\mathcal{L}:C_{LT}\rightarrow C_{LT}$ given by

[TABLE]

where $C_{LT}$ is the ordered Banach space of all $LT$ periodic functions form $\mathbb{R}$ to $\mathbb{R}^{2}$, which is equipped with the maximum norm. $\phi(s)\in C_{LT}$ is the initial distribution of infectious individuals in this periodic environment, and $Y(t,s)$, $t\geq s$ is the evolution operator of the linear periodic system:

[TABLE]

that means, for each $s\in\mathbb{R}$, the 2x2 matrix $Y$ satisfies

[TABLE]

$\mathcal{L}\phi$ is the distribution of accumulative new infections at time $t$ produced by all those infected individuals $\phi(s)$ introduced before $t$, with kernel $K(t,a)=Y(t,t-a)F(t-a)$. The coefficient $K_{i,j}(t,a)$ in row $i$ and column $j$ represents the expected number of individuals in compartment $I_{i}$ that one individual in compartment $I_{j}$ generates at the beginning of an epidemic per unit time at time $t$ if it has been in compartment $I_{j}$ for $a$ units of time, with $I_{1}=E,I_{2}=I$ [2].

Let $r_{0}>0,$ $r_{0}$ is an eigenvalue of $\mathcal{L}$ if there is a nonnegative eigenfunction $v(t)\in C_{LT}$ such that

[TABLE]

Therefore, the basic reproduction number is defined as

[TABLE]

the spectral radius of $\mathcal{L}$. The basic reproduction number can be evaluated by several numerical methods and approximations ([2], [10],[13]).

3 The threshold dynamics of $R_{0}$

3.1 Disease extinction

Theorem 2**.**

Let $\mathcal{R}_{0}$ be defined as (15), then the disease free periodic solution $(\hat{S}(t),0,0)$ is asymptotically stable if $\mathcal{R}_{0}<1$ and unstable if $\mathcal{R}_{0}>1$ .

Proof.

We use theorem 2.2 of [16], and check conditions (A1)-(A7). Conditions (A1)-(A5) are clearly satisfied from definition of $\mathcal{F}$ and $\mathcal{V}$ given in section 2. We prove only condition (A6) and (A7). Define

[TABLE]

and let $\Phi_{M}(t)$ be the monodromy matrix of system

[TABLE]

• (A6)

$\rho(\Phi_{M}(LT))<1$.

Let $\Psi_{M}$ be a fundamental matrix for system $\frac{dz}{dt}=M(t)z$, with $M$ defined as before and $LT$ periodic, the monodromy matrix $\Phi_{M}(LT)$ is given by $\Phi_{M}(LT)=\Psi_{M}^{-1}(0)\Psi_{M}(LT)$. The general solution of (16) is

[TABLE]

so $\Psi_{M}=\exp(-\int_{0}^{t}(\mu+r(s)+\delta)ds)$ and $\Psi_{M}^{-1}=\exp(\int_{0}^{t}(\mu+r(s)+\delta)ds)$. Note that $\Psi_{M}^{-1}(0)=1$, so $\Phi_{M}(LT)=\Psi_{M}(LT)$

[TABLE]

Due to the fact that $\Phi_{M}(LT)$ is a constant, its eigenvalue is itself and $\rho(\Phi_{M}(LT))<1$ for $\mu,\delta,r(s)>0$.

• (A7)

$\rho(\Phi_{-V}(LT))<1$.

Solving the system $\frac{dz}{dt}=-V(t)z$, we arrive to the general solution

[TABLE]

so

[TABLE]

Computing $\Phi_{-V}(LT)=\Psi_{-V}^{-1}(0)\Psi_{-V}(LT)$ we have

[TABLE]

Clearly, $\rho(\Phi_{-V}(LT))=\max\{e^{-(\mu+\sigma)LT},e^{-(\mu+\gamma)LT}\}<1$ for $\mu,\gamma,\sigma>0$.

∎

Note 1**.**

Due to $\Psi_{A}$ is a fundamental solution of a periodic system, we can always choose it such that $\Psi(0)=I$, so the monodromy matrix satisfies $\Phi_{A}(LT)=\Psi_{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.

Let $A(t)$ continuous, cooperative, irreducible and $\omega-$periodic $k\times k$ matrix function, and $\Psi_{A}(t)$ the fundamental matrix of system $x^{\prime}(t)=A(t)x(t)$. Denote by $\rho(\Psi_{A}(\omega))$ the spectral radius of $\Psi_{A}(\omega)$ .

Lemma 1**.**

Let $p=\frac{1}{\omega}\ln\rho(\Psi_{A}(\omega))$. Then there exists a positive, $\omega$-periodic function $v(t)$ such that $e^{pt}v(t)$ is a solution of $x^{\prime}(t)=A(t)x(t)$ (see proof in Lemma 2.1 of [19]).

Lemma 2**.**

Function $f(I)$ of model (1) satisfy that $f(I)\leq f^{\prime}(0)I$, $\forall I\geq 0$.

Proof.

Using assumptions on function $f$ we have

[TABLE]

so function $f(I)/I$ decreases $\forall I>0$ and then $\frac{f(I)}{I}\leq\lim_{I\rightarrow 0^{+}}\frac{f(I)}{I}=f^{\prime}(0)$ ∎

Lemma 3**.**

Let $(S(t),E(t),I(t))$ a solution of system (2) with initial conditions $(S_{0},E_{0},I_{0})\geq 0$, and $(\hat{S}(t),0,0)$ the infection free periodic solution of (2), then

[TABLE]

Proof.

Proof is similar to Lemma 4.1 of [11]. $S(t)$ satisfies first equation of system (2),then

[TABLE]

Let $X(t)=S(t)-\hat{S}(t)$, then

[TABLE]

Using Gronwall’s inequality $X(t)\leq X(0)e^{-\int_{0}^{t}(\mu+r(s)+\delta)ds},$ so

[TABLE]

Applying limit in both sides, we obtain $\limsup_{t\rightarrow\infty}S(t)-\hat{S}(t)\leq 0$. ∎

Now, we are able to enunciate our theorem for global stability of infection free periodic solution.

Theorem 3**.**

The infection free periodic solution $(\hat{S}(t),0,0)$ of system (2) is globally asymptotically stable if $\mathcal{R}_{0}<1$.

Proof.

From theorem (2) we have $(\hat{S}(t),0,0)$ is unstable for $\mathcal{R}_{0}>1$ and asymptotically stable for $\mathcal{R}_{0}<1$, so it is sufficient to prove that any solution $(S(t),E(t),I(t))$ with non-negative initial conditions $(S_{0},E_{0},I_{0})$ approaches to $(\hat{S},0,0)$.

Let $\epsilon>0$, from Lemma (20) we have

[TABLE]

so there exist a $N>0$ such that for all $t_{1}>N$

[TABLE]

this implies that $\sup_{t\geq t_{1}}(S(t)-\hat{S}(t))<\epsilon+L\leq\epsilon$. Then, from definition of supremum we have for all $t>t_{1}$

[TABLE]

Then, we have proved that for all $\epsilon>0$ we can find a $t_{1}>0$ such that $S(t)<\epsilon+\hat{S}(t)$ for all $t>t_{1}$.

Now, using lemma (2) , for $\epsilon>0$ we can find a $t_{1}>0$ such that for $t>t_{1}$

[TABLE]

We consider the following perturbated sub-system:

[TABLE]

which can be rewritten as

[TABLE]

with $F(t),V(t)$ defined in (10) and

[TABLE]

Matrix $(F-V+\epsilon H)(t)$ is LT-periodic, cooperative, irreducible and continuous. Using lemma (1), if $q=\frac{1}{LT}\ln\rho(\Psi_{F-V+\epsilon H}(LT))$ then there exist a positive and LT-periodic function $v(t)=(v_{1}(t),v_{2}(t))^{T}$ such that $e^{qt}v(t)$ is solution of system (23). Note that for all $k>0$, function $ke^{q(t-t_{i})}v(t-t_{i})$ is also a solution of system (23) with initial condition $kv(0)$ at $t=t_{i}$ .

Choose a $\bar{t}>t_{1}$ and $\alpha_{1}>0$ such that $(E(\bar{t}),I(\bar{t}))^{T}\leq\alpha_{1}v(0)$, then from (22)

[TABLE]

and using comparison principle (see for instance [14] theorem B.1), we have $(E(t),I(t))^{T}\leq\alpha_{1}e^{q(t-\bar{t})}v(t-\bar{t})$ for all $t>\bar{t}$.

From theorem 2.2 of [16], $\mathcal{R}_{0}<1$ iff $\rho(\Phi_{F-V}(LT)<1$. By the continuity of the spectrum for matrices (see [7], Section II.5.8 ) we can choose $\epsilon>0$ small enough that $\rho(\Phi_{F-V+\epsilon H}(LT)<1$ and then $q<0$ (see note (1) ). So, using positivity of solutions and comparison:

[TABLE]

And similarly for I. We obtain

[TABLE]

We need only prove that $S(t)$ approaches to $\hat{S}$. At infection free solution $\hat{R}(t)=N-\hat{S}(t),$ where $\hat{R}$ satisfies equation

[TABLE]

So $R(t)=N-S(t)-E(t)-I(t)$ satisfies

[TABLE]

Let $\epsilon_{1}>0$ arbitrary and $r_{max}=\max_{u\in[0,LT]}r(u)$. Due to (25) we can find a $t_{2}>0$ such that $I(t)<\epsilon_{1}$ for $t>t_{2}$, moreover we can find a $t_{3}>0$ such that $S(t)\leq\hat{S}(t)+\epsilon_{1}$ for $t>t_{3}$. Then, let $t_{4}=\max\{t_{2},t_{3}\}$, we have for $t>t_{4}$

[TABLE]

Multiplying in both sides by $e^{(\mu+\delta)}t$ and integrating from $t_{4}$ to $t$ we obtain

[TABLE]

So, $\limsup_{t\rightarrow\infty}(R-\hat{R})(t)\leq\frac{\epsilon_{1}(r_{max}+\gamma)}{\mu+\delta}$, where $\frac{\epsilon_{1}(r_{max}+\gamma)}{\mu+\delta}$ is arbitrarily small. Then $\limsup_{t\rightarrow\infty}(R-\hat{R})(t)\leq 0$ and using similar arguments to used for $S$ for $\epsilon_{3}>0$ we can find a $t_{5}>0$ with $R(t)\leq\hat{R}(t)+\frac{\epsilon_{3}}{2}$ for $t>t_{5}$. Also, from (25) we can find $t_{6}>0$ with $E(t)+I(t)<\frac{\epsilon_{3}}{2}$ for $t>t_{6}$, so for $t>\max\{t_{5},t_{6}\}$ we have

[TABLE]

Or equivalently, $S(t)-\hat{S}(t)\geq-\epsilon_{2}$, with $\epsilon$ arbitrarily small and this implies that $\liminf_{t\rightarrow\infty}(S-\hat{S})(t)\geq 0$. We conclude by comparison and using lemma (20) that $\lim_{t\rightarrow\infty}S(t)=\hat{S}(t)$ completing the proof.

∎

Theorem (3) shows that disease will completely die as long as $\mathcal{R}_{0}<1$. So, reducing and keeping $\mathcal{R}_{0}$ 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 [20].

In this section we consider the dynamics of the periodic model when $\mathcal{R}_{0}>1$. We will show that actually, $\mathcal{R}_{0}$ is a threshold parameter for the extinction and the uniform persistence of the disease. The results are inspired by [5], [13], [18] and [19].

Let $P:X\rightarrow X$ be the Poincaré map associated with system (2), that is

[TABLE]

where $X$ is defined in (3) and $\phi(t,x_{0})$ is the unique solution of system (2) with $\phi(0,x_{0})=x_{0}$. We define the following sets:

[TABLE]

Note that $\partial X_{0}$ is not the boundary of $X_{0}$, but it is a standard notation of persistence theory.

Lemma 4**.**

The set $X_{0}$ is positively invariant under system (2).

Proof.

Let $x_{0}=(S_{0},E_{0},I_{0})\in X_{0}$, ie $E_{0}>0,I_{0}>0$ and $\phi(t,x_{0})=(S(t),E(t),I(t))$ the solution of (2) with $\phi(0,x_{0})=x_{0}$. Due to non negativity of solutions and assumptions on function $\beta(t)$ and $f(I)$, we have:

[TABLE]

Using comparison theorem (see for instance [14] Appendix B.1) we have for all $t>0$:

[TABLE]

Similarly,

[TABLE]

so,

[TABLE]

Therefore, $\phi(t,x_{0})$ remains on $X_{0}$ for all $t>0$. ∎

Set:

[TABLE]

To use persistence theory developed in [20], we show that

[TABLE]

Let $x_{0}=(S_{0},0,0)\in X$ and $(S(t),E(t),I(t))$ the solution that passes trough that initial condition. We have that $\phi(t,x_{0})=(S_{1}(t),0,0)$, with $S_{1}(t)$ solution of (4) and $S_{1}(0)=S_{0}$ is a solution that satisfies the initial condition. By uniqueness of solutions we have $E(t)=0=I(t)$ $\forall t\geq 0$, so $x_{0}$ lives on $M_{\partial}$.

Now, if $x_{0}\in M_{\partial}$ we want $x_{0}=(S_{0},0,0)$. We prove an equivalent sentence: if $x_{0}\in\partial X_{0}\backslash\{(S,0,0):S\geq 0\}$ then it does not belong to $M_{\partial}$. Consider an initial point $x_{0}=(S_{0},E_{0},I_{0})\in\partial X_{0}\backslash\{(S,0,0):S\geq 0\}$, then $E_{0}>0,I_{0}=0$ or $E_{0}=0,I_{0}>0$. Suppose $E>0$ and $I_{0}=0$, then $\phi(t,x_{0})$ holds

[TABLE]

By continuity of $E(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$, $\phi(t,x_{0})\in X_{0}$. Using invariance of $X_{0}$ ( Lemma (4) ) we have $\phi(t,x_{0})\in X_{0}$ for all $t>1$. Finally, for an $m>0$ such that $mLT>1$ we have $P^{m}(x_{0})=\phi(mLT,x_{0})\in X_{0}$ and this implies (29). By discussion in section 2 is clear that there is one fixed point of $P$ in $M_{\partial}$: $M_{0}(\hat{S}(0),0,0)$ ( [12] ).

Now, we are in position to introduce the following result of uniform persistence of the disease.

Theorem 4**.**

Let $\mathcal{R}_{0}>1,$ then there exists an $\epsilon>0$ such that any solution $(S(t),E(t)I(t))$ of (2) with initial values $(S(0),E(0),I(0))\in X_{0}$ satisfies

[TABLE]

Proof.

We first prove that $P$ is uniformly persistent (see definition 1.3.2 from [20] ) with respect to $(X_{0},\partial X_{0})$, because this implies that the solution of (2) is uniformly persistent with respect to $(X_{0},\partial X_{0})$ ( [20], theorem 3.1.1 ). Clearly, $X_{0}$ is relatively open in $X$, so $\partial X_{0}$ is relatively closed.

Define

[TABLE]

we show that $W^{s}(M_{0})\cap X_{0}=\emptyset.$

By theorem 2.2 of [16], $\mathcal{R}_{0}>1$ iff $r(\Psi_{F-V}(LT))>1$. Choose an $\eta>0$ small enough with the property $\hat{S}(t)-\eta>0,\forall t>0$ (see appendix (B) ). For $\alpha>0$, let us consider the following perturbed equation:

[TABLE]

System above admits a unique positive $LT-$periodic solution of the form:

[TABLE]

whith $\hat{S}(t,0)=\hat{S}(t)$ and which is globally attractive in $\mathcal{R}_{+}$ (see appendix (D) ) with

[TABLE]

Since $\hat{S}(0,\alpha)$ is continuous in $\alpha$, then for all $\epsilon>0$ there is a $\delta>0$ such that for $|\alpha|<\delta$ we have $|\hat{S}(0,\alpha)-\hat{S}(0,0)|<\epsilon$. Moreover, by continuity of solutions with respect to initial values we can find for all $\bar{\eta}>0$ an $\bar{\epsilon}>0$ such that if $|\hat{S}(0,\alpha)-\hat{S}(0,0)|<\bar{\epsilon}$ then

[TABLE]

Therefore, for $\eta$ established before, we can find $\alpha$ small enough such that $\hat{S}(t,\alpha)>\hat{S}(t)-\eta$, $\forall t>0$.

Again, by continuity of solutions with respect to initial values, for this small $\alpha>0$, there exists a $\delta>0$ such that for all $(S_{0},E_{0},I_{0})\in X_{0}$ with $\|(S_{0},E_{0},I_{0})-M_{0}\|\leq\delta$ then $\|\phi(t,(S_{0},E_{0},I_{0}))-\phi(t,M_{i})\|<\alpha$, $\forall t\in[0,LT]$.

We now claim that

[TABLE]

By contradiction, suppose that

[TABLE]

Without loss of generality, we can assume that $\|P^{m}(S_{0},E_{0},I_{0})-M_{0}\|<\delta$ for all $m\geq 0$ (see appendix (A)). From above discussion, $\|\phi(t,P^{m}(S_{0},E_{0},I_{0}))-\phi(t,M_{0})\|<\alpha$, $\forall m\geq 0$ and $t\in[0,LT]$.

For any $t\geq 0,$ let $t=mLT+t_{1}$, where $t_{1}\in[0,LT)$ and $m=[\frac{t}{LT}]$ is the greatest integer less than or equal to $\frac{t}{LT}$. Then, we get

[TABLE]

If we set $\phi(t,(S_{0},E_{0},I_{0}))=(S(t),E(t),I(t))$ , then we have $E(t)\leq\alpha,I(t)\leq\alpha$, $\forall t\geq 0$, and from first equation of (2) and lemma (2) we arrive to:

[TABLE]

Which is exactly the equation in (31). Since the unique periodic solution of (31) is globally attractive in $\mathcal{R}_{+}$, we have for $\bar{S}(t,\alpha)$ solution of (31) that $\lim_{t\rightarrow\infty}\bar{S}(t,\alpha)=\hat{S}(t,\alpha)$. So for $\eta$ given before there exists $T>0$ such that for all $t\geq T$

[TABLE]

or equivalently $\bar{S}(t,\alpha)>\hat{S}(t,\alpha)-\eta$. Moreover, from previous analysis, $\hat{S}(t,\alpha)-\eta>\hat{S}(t)-\eta$, therefore, using comparison principle on (37) we arrive to

[TABLE]

for $t>T$.

Due to $E(t),I(t)\leq\alpha$, and $\alpha$ is fixed small, we can take $\alpha<\epsilon^{*}$ and use assumption (A5) in introduction and hence (see appendix (C) )

[TABLE]

where $F,V,$ are defined in (10), $H$ is defined as (24) and

[TABLE]

By theorem 2.2 of [16], we have $\mathcal{R}_{0}>1$ iff $\rho(\Phi_{F-V}(LT))>1$. By continuity of spectrum (see [7] Section II ) we can find $\alpha,\epsilon$ such that

[TABLE]

Consider the auxiliar system

[TABLE]

then, using lemma (1) there exist a solution of (39) with the form $e^{p_{2}t}v_{2}(t)$, with $p_{2}=\frac{1}{LT}\ln(\rho(\Phi_{F-V-\eta H-\alpha K}(LT)))>0$. Choose a $t_{2}>T$, and a small number $\alpha_{2}>0$ such that $(E_{2}(t_{2}),I_{2}(t_{2}))^{T}\geq\alpha_{2}v_{2}(0)$. Using comparison principle we get $(E(t),I(t))\geq\alpha_{2}v_{2}(t-t_{2})e^{p_{2}(t-t_{2})},$ which implies $E(t)\rightarrow\infty$ and $I(t)\rightarrow\infty$. This leads a contradiction.

The above claim shows that $P$ is weakly uniformly persistent with respect to $(X_{0},\partial X_{0}).$ Note that $P$ has a global attractor $\hat{S}(0)$ ( (20) ). It follows that $M_{0}$ is an isolated invariant set in $X$, $W^{s}(M_{0})\cap X_{0}=\emptyset$. Every orbit in $M_{\partial}$ converges to $M_{0}$ and $M_{0}$ is acyclic. By the aciclity theorem on uniform persistence for maps( [20] Theorem 1.3.1 and Remark 1.3.1 )is follows that $P$ is uniformly persistent with respect to $(X_{0},\partial X_{0})$, that is, there exists $\epsilon>0$ such that any solution of (2) satisfies $\lim_{t\rightarrow\infty}E(t)\geq\epsilon,$ $\lim_{t\rightarrow\infty}I(t)\geq\epsilon.$

∎

4 Numerical simulations

In this section we provide some numerical simulations to illustrate the results obtained in our theorems and compare with previous results.

To improve previous models used in references, we use a particular function

[TABLE]

which includes the case $f(I)=I$ used in [5]. One can check that function (40) satisfies conditions A1)-A5). Using this function, system (2) is rewritten as:

[TABLE]

Set an initial population $N=2,200,000$ and take time $t$ in years. Suppose $\mu=0.02$ per year, corresponding to an average human life time of 50 years. Following [5] take the parameters as follows: $\sigma=38.5$ per year, $\gamma=100$ per year, $p=0.85,\delta=0$. Choose the periodic transmission as $\beta(t)=\beta_{0}+0.0002\cos(2\pi t)$, with $\beta_{0}$ the transmission parameter, and the periodic vaccination rate $r(t)=0.1+0.004\cos(2\pi 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 $\rho(W(LT,0,\lambda))=1$ (see theorem 2.1 of [16] ). In order to compare with previous works, we firstly approximate the basic reproduction number with its average value $\mathcal{R}_{0}^{T}$ (used by several authors as [9] and [17]), so define:

[TABLE]

where $V$ is given by (10) and

[TABLE]

with $[\beta],[\hat{S}]$ the average of functions $\beta,\hat{S}$ defined as $[\beta]=\frac{1}{LT}\int_{0}^{LT}\beta(t)dt,$ $\hat{S}=\frac{1}{LT}\int_{0}^{LT}\hat{S}(t)dt$. Computing each average we obtain $R_{0}^{T}=549.6702634\beta_{0}$, so $R_{0}^{T}>1$ for $\beta_{0}\in(0.001819272510,\infty).$

Following theorem 2.1 of [16], to compute $\mathcal{R}_{0}$, let $W(t,s,\lambda),t\geq s,$ the evolution operator of the system

[TABLE]

ie, for each $\lambda\in(0,\infty)$, $\frac{dW(t,s,\lambda)}{dt}=\left(-V(t)+\frac{F(t)}{\lambda}\right)W(t,s,\lambda),\forall t\geq s$ and $W(s,s,\lambda)=I_{2\times 2}$, then $\mathcal{R}_{0}>0$ is the unique solution of $\rho(W(LT,0,\lambda))=1.$

Example 1**.**

To illustrate our results, first fix $\beta_{0}=0.0018$. Computing $R_{0}^{T}$ we have $R_{0}^{T}=0.9894064741$, which is a first approximation of $R_{0}$. To solve system (43) numerically, we substitute the terms of expression of $\hat{S}(t)$ in (5):

[TABLE]

Due to we can not compute analytically the term

[TABLE]

we approach $\hat{S}(t)$ using Taylor expansion around 0 (remember that we want so solve $\rho(W(LT,0,\lambda))=1,$ where $LT=1$), so even when we can not find an explicit expression for $\hat{S}(t)$, the Taylor expansion is a good way to estimate it in $(0,1)$ .

Setting an initial value $\lambda_{0}=0.98$ and letting $\lambda_{i}=\lambda_{0}+i(0.0001)$, we solve system (43) numerically for each $\lambda_{i}$ (using initial conditions $w(0)=(1,0)$ and $w(0)=(0,1)$, to satisfy $W(0,0)=I_{2\times 2}$ ), and compute $\rho_{1}=\rho(W(LT,0,\lambda_{i}))$ until $\rho_{1}\sim 1$  . With previous process we arrive to $\rho_{1}=1.00120166209265$ for $\lambda=0.9872$ and $\rho_{1}=0.997826338969630$ for $\lambda=0.9873$, therefore $\mathcal{R}_{0}\in(0.9872,0.9873)$. Using a finer step size $0.0000001$ to have more accuracy, we arrive to $\mathcal{R}_{0}\sim 0.9872355<1$.

Set initial values as $S(0)=1,500,000,E(0)=400,000,I(0)=40,000,R(0)=N-(S(0)+E(0)+I(0))$.

We use Matlab algorithms to graph the solution of system (41) with these initial conditions. Figure (2) shows the results. We can see that $I(t),E(t)$ goes to cero, while $S(t),R(t)$ tend to stabilize, also $S(t)$ is tending to $\hat{S}(t)$ with values between 54,000 and 56,000 (see figure (1) ), this shows the results obtained in theorem (3) .

Example 2**.**

Now, choose $\beta_{0}=0.005$. As we can see in figure (3), the solutions of system remain persistent when $t$ tends to infty, this fact suggest that $\mathcal{R}_{0}>1$ from theorem (30). In fact, if we compute the basic reproduction number and its average (using the process described in example 1), $\mathcal{R}_{0}^{T}=3.298021580$ and $\mathcal{R}_{0}\in(2.7456,2.7457)$, therefore it is bigger than one. In fact, this shows the results of persistence obtained in theorem (30).

5 Conclusion

In this paper we presented a model with seasonal fluctuation with a general incidence function that includes the bilinear case studied by [5]. We proved that $\mathcal{R}_{0}$ is a threshold parameter for stability and persistence of system, giving also some numerical simulations that show these results.

Several authors (for example [9] and [17]) define $\mathcal{R}_{0}^{T}$ as the basic reproduction number, but we can see in numerical simulations that $\mathcal{R}_{0}^{T}$ is not equal to $\mathcal{R}_{0}$ defined by [16], which is a real threshold parameter for extinction and persistence of disease.

To obtain the estimation of $R_{0}$ we used a code in Maple, which is based on numerical computing $\rho_{1}=\rho(W(LT,0,\lambda_{i}))$ until $\rho_{1}\sim 1$, where $\lambda_{i}=\lambda_{0}+\Delta_{\lambda}i$, $\Delta_{\lambda}$ is the step size and the initial estimation $\lambda_{0}$ is taken as $R_{0}^{T}-\epsilon$. The Maple code is available to anyone who wants to use it.

6 Acknowledgements

This article was supported in part by Mexican SNI under grant 15284 and 33365.

Appendix A Appendix A: Assumption used in (30)

Let $f(m):=\|P^{m}(S_{0},E_{0},I_{0})-M_{i}\|$. If

[TABLE]

then we have $L=\lim_{m\rightarrow\infty}\left(\sup_{n\geq m}f(n)\right)<\delta$. For all $\epsilon>0$ there exists a $M_{\epsilon}>0$ such that if $m\geq M_{\epsilon}$ then $-\epsilon<\sup_{n\geq m}f(n)-L<\epsilon$. Particularly, for $\epsilon=\frac{\delta-L}{2}>0$ we have

[TABLE]

or equivalently, $\sup_{n\geq m}f(n)<\delta$ for $m\geq M_{\delta-L}$. Moreover, for all $n\geq m$ with $m\geq M_{\delta-L}$ we have $f(n)<\sup_{n\geq m}f(n)<\delta$. Therefore, $\|P^{n}(S_{0},E_{0},I_{0})-M_{i}\|<\delta$, $\forall n\geq M_{\delta_{L}}$.

We can take $(S^{1}_{0},E^{1}_{0},I^{1}_{0})=P^{M_{\delta-L}}(S_{0},E_{0},I_{0})$ as initial condition and therefore,

[TABLE]

making our assumption valid.

So, we can assume without loss of generality that $\|P^{m}(S_{0},E_{0},I_{0})-M_{i}\|<\delta$ for all $m\geq 0$ .

Appendix B Appendix B

Note that $\hat{S}(t)$ has a positive minimum $min$(is periodic, positive and continuous, so it is bounded for $t\in[0,LT]$ and then for all $t>0$) and we can choose a $min>\eta>0$, sufficiently small such that $\hat{S}(t)-\eta>0$.

Appendix C Appendix C: expression 36

From system (2) $\frac{dE}{dt}=\beta(t)Sf(I)-(\mu+\sigma)E,$ with $S(t)>\hat{S}(t)-\eta$ for $t>T$, so

[TABLE]

Using assumption (A5) for $f(I)$ and positivity of $\hat{S}(t)-\eta$, we have also

[TABLE]

So,

[TABLE]

Due to $0<I<\alpha$ and $f^{\prime\prime}(0)\leq 0$, then $I^{2}<\alpha I$ and $f^{\prime\prime}(0)I^{2}\geq f^{\prime\prime}(0)\alpha I$, applying this we arrive to

[TABLE]

This expression can be written as (39).

Appendix D Appendix D: Auxiliar from theorem 4

The system used in the proof of theorem (30) is

[TABLE]

Solving the equation above, we arrive to the general solution

[TABLE]

where $p(s)=\beta(s)f^{\prime}(0)\alpha+\mu+r(s)+\delta$. We shall examine the behaviour of an arbitrary solution $S$. For each $n=0,1,...$ we can use an initial time $\bar{t}_{0}=t_{0}+nLT$ with initial point $\bar{S}(\bar{t}_{0})$ and see that:

[TABLE]

Due to $p(s)$ is a periodic function, then

[TABLE]

where $s-nLT\geq t_{0}$. Then

[TABLE]

And using the change of variable $u=s-LT$, then

[TABLE]

Equation (46) gives a recursive relationship between the solution at $t_{0}+nLT$ and after $LT$ times. If we set $S_{n}=S(t_{0}+nLT)$,then for each solution $S$ this relationship is described by:

[TABLE]

with $F$ the right side of (46). If we take $S_{i}$ and $S_{j}$, two different values of ${S_{n}}$, then

[TABLE]

Then, $F(S)$ is a contracting map and by Banach fixed point theorem $F$ has a unique fixed point $S_{i}$ such that $S_{i+1}=F(S_{i})=S_{i}$, or equivalently, $S(t_{0}+iLT)=S(t_{0}+(i+1)LT)$. This fixed point can be found for any $S$ that is solution of differential equation with arbitrary initial condition $S(t_{0})$ at any time $t_{0}$. The fixed point has the form:

[TABLE]

So, define the function

[TABLE]

$S^{*}$ is a periodic function with period $LT$ and is continuously differentiable with respect to $t$. One can check (by computing the derivative) that $S^{*}(t)$ is a solution of differential equation, so by existence and uniqueness of solutions it can be rewritten as (33) with initial condition (34).