Extinction and permanence of a general non-autonomous discrete-time SIRS epidemic model

1.   Introduction

In recent years, mathematical epidemic models have been extensively investigated in order to understand the dynamics and control the spreading of various diseases, see for instance ^{[4,5,6,7,11,13,14,15,19,20,23,24,25,26,27]} and the references therein. Related to this work are the following studies. Zhang et al. ^[26] investigated a non-autonomous SIRS epidemic model with a standard incidence rate and distributed delays of the form

where $S(t)$, $I(t)$, and $R(t)$ are the numbers of susceptible, infectious, and recovered individuals at time $t$, respectively. The distributed delays are used to model the infection mechanisms of some diseases, where infected individuals may not be infectious until some time after becoming infected, and that the infectivity is a function of the duration since infection, up to some maximum duration. Various (time-dependent) parameters in (1.1) are defined as follows: $\varLambda$ is the growth (or recruitment) rate of the population, $\beta$ is the daily contact rate, that is the average number of contacts per day, $\mu_i$, for $i = 1, 2, 3$, are the instantaneous per capita mortality rates of $S$-, $I$-, and $R$-classes, respectively, $\gamma$ is the rate that the recovered individual loses immunity and returns to be susceptible, $\xi$ is a time taken for an infected individual to become infectious and $p(\xi)$ is the distributed proportion of the population taking time $\xi$ after being infected to become infectious, $\tau$ is the infected period, and $k$ is the recovery rate that can incorporate basic medical treatment and prevention of the disease. The main result in ^[26] is that they obtain two threshold values $R_\star$ and $R^\star$, which depend on the parameters of the model, so that the disease is permanent if $R_\star > 1$ while if $R^\star < 1$ then the disease extincts. Global behavior of the model is also studied using the Lyapunov functional method.

Enatsu et al. ^[3] investigated an autonomous version of (1.1), i.e. $\varLambda, \beta, \mu_1$, $\gamma, \mu_2, k, \mu_3$ are constants, but with a general nonlinear incidence rate and distributed delays of the form

where the nonlinear function $g(I)$ satisfies

(A1) $g(I)$ is continuous and monotone increasing on $[0, \infty)$ with $g(0) = 0$, and

(A2) $g(I)/I$ is monotone decreasing on $(0, \infty)$ with $\lim_{I\to0^+}g(I)/I = 1$.

Using the Lyapunov functional method, the authors derived the basic reproduction number $R_0$ and established sufficient conditions of the rate of immunity loss for the global asymptotic stability of an endemic equilibrium for the model.

To mitigate the impact of the disease, some control and prevention should be included into modeling the disease. ^[16,22] investigated SIR and SIRS epidemic models where the impact of health care resources especially hospital beds is included, so that the recovery rate is a nonlinear function of $I$ of the form

Here, according to ^[22], $k_1$ is the maximum per capita recovery rate when health care resources are adequate and the number of infected people is low, $k_0$ is the minimum per capita recovery rate due to the lack of basic clinical resources, and $b$ is a parameter measuring the availability of hospital beds. Complex dynamics of the models are then derived in those papers.

Based on the above discussions, it is interesting to consider the following general non-autonomous SIRS epidemic model

where $S(t)$, $I(t)$, and $R(t)$ and $\varLambda(t), \beta(t), \mu_1(t), \gamma(t), \mu_2(t), \mu_3(t), p(\xi)$ are the same as explained in (1.1) and (1.2). In the model (1.4), however, we relax the assumptions (A1), (A2) and also we do not require the nonlinear recovery rate $k(I, t)$ to have the form (1.3). Namely, for $k(I, t)$, we assume that there are positive constants $k_0, k_1$ such that

for all $I\geq0$ and $t\geq0$. For $g(I)$, it is assumed to be continuous, $g(I)\geq0$ on $[0, \infty)$, $g(0) = 0$ and $g(I) > 0$ for $I > 0$, and also that $\lim_{x\to0^+}(g(x)/x)$ exists and is positive. Now we define the function

In epidemiology, $f$ is the transmission function of the disease. Observe that $f$ is continuous and $f(I) > 0$ on $[0, \infty)$. Clearly, $g(I) = If(I)$, so the incidence rate appearing in (1.4) can be expressed as

Let us mention some specific examples of incidence rates. If $f(I) = 1$ is constant, then $g(I)S = IS$ is simply the standard (or bilinear) incidence rate; if $f(I) = \alpha_3/(1+\alpha_1 I)$, then $g(I)S = \alpha_3IS/(1+\alpha_1 I)$ is the saturated incidence rate ^[2]; if $f(I) = \alpha_3/(1+\alpha_2 I+\alpha_1I^2)$, then $g(I)S = \alpha_3IS/(1+\alpha_2 I+\alpha_1I^2)$ is the non-monotone incidence rate ^[8,19,20].

An interesting and important question is whether the results obtained for the continuous-time SIRS model (1.4) can be extended to the corresponding discrete-time epidemic model. In this work, we investigate the following discrete-time non-autonomous SIRS epidemic model: For a step size $h > 0$ and the maximum infectious period $m\in\mathbb{N}$, consider

where $S_n, I_n$, and $R_n$ are the numbers of susceptible, infectious, and recovered individuals at time $n$, respectively, and $p_j\ge 0\, (j = 0, 1, \dots, m)$ is the distributed proportion of the population taking time $j$ to become infectious, and $\varLambda(n)$, $\beta(n)$, $\mu_1(n)$, $\gamma(n)$, $\mu_2(n)$, $k(I_n, n)$, $\mu_3(n)$ are the discretized values of the corresponding continuous-time functions explained in (1.1), (1.2), and (1.5). We have obtained (1.8) from (1.4) by applying Mickens's nonstandard finite difference discretization (see ^[10,17,27]; see also Appendix for a brief introduction). Let us list some known results for (1.8). A global attractivity for the autonomous version of (1.8) with $k = k(n)$ and a bit more general $f$ is investigated in ^[17]. In ^[27], the threshold conditions are obtained for (1.8) with the non-delay standard incidence rate.

This paper is organized as follows. In Section 2 we introduce notations and basic results used in this work. Then we derive the positivity and boundedness of solutions $(S_n, I_n, R_n)$ of (1.8) in Section 3. In Section 4, a threshold condition is obtained so that the permanence of disease for the model (1.8) is guaranteed. Then, we derive in Section 5 a threshold condition so that the extinction of the disease for (1.8) holds. In Section 6, some numerical simulations are presented. Finally, a discussion is given in Section 7.

2.   Notations and preliminaries

Notation. In this work, we denote

● $G_n = \beta(n)\sum\limits_{j = 0}^mp_jf(I_{n-j})I_{n-j}$ and $K_n = k(I_n, n)$.

● $\varPi_1(n) = 1+h\mu_1(n)$, $\varPi_2(n) = 1+h\mu_2(n)$, and $\varPi_3(n) = 1+h(\mu_3(n)+\gamma(n))$.

● For a sequence $\{a(n)\}$, denote $a^u = \limsup\limits_{n\to\infty}a(n)$, $a^l = \liminf\limits_{n\to\infty}a(n)$, $a^M = \sup\limits_{n\in \mathbb{N}_0}a(n)$, $a^L = \inf\limits_{n\in \mathbb{N}_0}a(n)$. Same notations are applied to functions as well, e.g. $f^M = \sup_{I\geq0}f(I)$ and $G^M = \sup_{n\in \mathbb{N}_0}G_n$.

The system (1.8) is said to be permanent if there exist constants $l_i, L_i > 0$ ($i = 1, 2, 3$) such that any solution $(S_n, I_n, R_n)$ satisfies

On the other hand, (1.8) is said to exhibit extinction provided $\lim_{n\to \infty}I_n = 0$ for any solution.

We assume the following conditions for (1.8). First of all, $\varLambda(n), \beta(n), \mu_1(n), \mu_2(n), \mu_3(n), \gamma(n)$ are bounded and positive for all $n\in \mathbb{N}_0$, $\sum_{j = 0}^{m}p_j = 1$, and there are constants $k_0 > 0, k_1 > 0$ such that $k_0\leq K(I, n)\leq k_1$ for all $I\geq0$ and $n\in \mathbb{N}_0$. We list all other assumptions that will be used in this work.

(H1) $f\geq0$ is a bounded continuous function on $[0, \infty)$.

(H2) $\mu_1(n)\leq\min\{\mu_2(n), \mu_3(n)\}$ for all $n\in \mathbb{N}_0$.

(H3) There exist $n_0\in \mathbb{N}$, $\chi\geq1$, and $0 < q < 1$ such that

(H4) There exist $\lambda > 0$ and $r\in \mathbb{N}\cup\{0\}$ such that

The following key result will be used in the study of disease-free states of (1.8) and also employed repeatedly in the investigation of the permanence and extinction phenomena. The proof is inspired by Lemma 2 in ^[21] and Lemma 2.2 in ^[24].

Lemma 1. Let $\{a(n)\}, \{b(n)\}$ be sequences such that $a(n) > 0$ and $0 < b(n) < 1$ for all $n\in \mathbb{N}_0$. Assume that $a^u < \infty$, $b^l > 0$, and that there exist $n_0$, $\chi > 0$, and $q\in(0, 1)$ such that

Let $x_n$ be a positive solution to the equation

with initial value $x_0\geq0$. Then the following results hold:

(1) There exist constants $x^\star\geq0, x_\star\geq0$ depending only upon $\{a(n)\}$ and $\{b(n)\}$ such that

independent of the initial condition $x_0$. In fact, we have

for any $n_1\geq n_0$, and the following estimate holds

Moreover, each fixed solution of (2.3) is globally uniformly attractive, i.e. if $x_n'$ is also a solution of (2.3) with initial value $x_0'\geq0$, then

(2) Suppose there exist constants $a_0 > 0$ and $r\in \mathbb{N}\cup\{0\}$ such that

Then $x^\star > 0$ and $x_\star > 0$.

(3) Suppose $y_n$ satisfies

Then

If in addition $y_0\leq x_0$ then $y_n\leq x_n$ for all $n$.

(4) Suppose $z_n$ satisfies

Then

If in addition $z_0\geq x_0$ then $z_n\geq x_n$ for all $n$.

(5) Let $\{a(n)\}, \{b(n)\}$ be as above and $\{ \tilde{b}(n)\}$ be another sequence satisfying $\tilde{b}(n)\in(0, 1)$ for all $n\in \mathbb{N}_0$, $\tilde{b}^l > 0$, and (2.2) holds with $b(n), \chi, q, n_0$ replaced by $\tilde{b}(n), \chi, q, \tilde{n}_0$, respectively. Assume also that there exist constants $\varepsilon > 0$ and $N_0\in \mathbb{N}$ such that

Let $\tilde{x}_n$ be a positive solution of the equation

Then there is a constant $C = C(a^u, x^\star, \chi, q) > 0$ such that

where according to (1), $\tilde{x}^\star = \limsup_{n\to \infty} \tilde{x}_n$ and $\tilde{x}_\star = \liminf_{n\to \infty} \tilde{x}_n$.

(6) For each $\varepsilon > 0$ small, there exist positive integers $N_\varepsilon$ and $D_\varepsilon$ such that if $n_1\geq N_\varepsilon$ and $n_2-n_1\geq D_\varepsilon$, then

Proof. (1) For $N > n_1$, it is directly to show that

where $A(N, n_1) = \sum_{n = n_1}^{N-1}a(n)b(n)b(n+1)\cdots b(N-1)$ and $B(N, n_1) = b(n_1)b(n_1+1)\cdots b(N-1)$. For fixed $n_1\geq n_0$, we have by assumption (2.2) that

and

so $\{A(N, n_1)\}_{N = n_1+1}^ \infty$ is bounded. Now $\lim_{N\to \infty}(x_{N}-A(N, n_1)) = \lim_{N\to \infty}B(N, n_1)x_{n_1} = 0$, it follows by the boundedness of $\{A(N, n_1)\}_{N = n_1+1}^ \infty$ that

and the two limits do not depend on $n_1$. This proves the first part and (2.4), (2.5).

For each $\varepsilon > 0$, by taking $n_1$ large enough, we get $\sup_{n\geq n_1}a(n)\leq a^u+\varepsilon$, hence

proving (2.6).

The last part of (1) is true because

(3) and (4) immediately follow from the proof of (1). We prove (2). By assumption $b^l > 0$, so there exist $b_0\in(0, 1)$ and $N_0\in \mathbb{N}$ such that $b(n)\geq b_0$ for all $n\geq N_0$. Fix $n_1\geq \max\{n_0, N_0\}$. Then we have for all $N\geq n_1+r+1$ that

Since $a_0b_0^{r+1} > 0$, we have $x^\star\geq x_\star = \liminf_{N\to\infty}A(N, n_1) > 0$.

(5) Recall the well-known fact^* that, for bounded sequences,

^* Let $c = \limsup| \tilde{x}_n-x_n|$, $a = \limsup \tilde{x}_n, b = \limsup x_n$. By the triangle inequality $\tilde{x}_n\leq| \tilde{x}_n-x_n|+x_n$, we get $a\leq c+b$, similarly, we also have $b\leq c+a$. Thus $|a-b|\leq c$. For the other fact, we simply use that $\liminf x_n = -\limsup(-x_n)$.

So it suffices to show that $\limsup_{n\to \infty}| \tilde{x}_n-x_n|\leq C\varepsilon$. For each $n$, we have

By the assumption (2.8), it follows that

Note that $\limsup_{n\to \infty}\varepsilon(a(n)+x_n)\leq\varepsilon(a^u+x^\star)$. Using part (3) and (2.6), we have

so $\limsup_{n\to\infty}| \tilde{x}_n-x_n|\le C\epsilon$, where $C = (a^u+x^\star)\chi q/(1-q) > 0$.

(6) We have $x_n\leq x^M: = \max x_n$ for all $n$. For $\varepsilon > 0$, there exists $N_\varepsilon\geq n_0$ such that if $n\geq N_\varepsilon$ then $x_{n}\geq x_\star-\frac{\varepsilon}{2}$. By (2.2), there exists $D_\varepsilon$ such that if $N-n_1\geq D_\varepsilon$ then $b(n_1)b(n_1+1)\cdots b(N)x^M\leq\frac{\varepsilon}{2}$. Now let $n_1\geq N_\varepsilon$ and $N\geq n_1+D_\varepsilon$. Then we have

as desired.

The following elementary result is employed throughout this work.

Lemma 2. The system (1.8) can be expressed as

where $\varPi_1(n) = 1+h\mu_1(n), \varPi_2(n) = 1+h\mu_2(n)$, $K_n = k(I_n, n)$, and $\varPi_3(n) = 1+h(\mu_3(n)+\gamma(n))$.

Proof. By (1.8), we have $S_{n+1} = h\varLambda(n)-hG_nS_{n+1}-h\mu_1(n)S_{n+1}+h\gamma(n) R_{n+1}+S_n$, which directly leads to the first expression. Similarly, we have $I_{n+1} = hG_nS_{n+1}-h(\mu_2(n)+k(I_n, n))I_{n+1}+I_n$ and $R_{n+1} = hk(I_n, n)I_{n+1}-h(\mu_3(n)+\gamma(n))R_{n+1}+R_n$, so the other two expressions follow.

3.   Positivity and boundedness

We prove the positivity of solutions of (1.8).

Proposition 3. Assume $f(I)\geq0$ for all $I\geq0$. Then every solution of the problem (1.8) is positive, that is $S_n > 0, I_n > 0, R_n > 0$ for all $n > 0$.

Proof. For simplicity, we omit the dependence of $\varPi_1, \varPi_2, \varPi_3, \gamma$ on $n$ in the following calculations. By (2.9), one can directly manipulate the identities to obtain

thus

From the initial condition in (1.8), it is easily seen that $G_0\geq0$ so $S_1 > 0$. Observe that $\varPi_i > 1$ for $i = 1, 2, 3$ and $(\varPi_2+K_n)\varPi_3-h\gamma K_n > 0$. Using (2.9), we also have $I_1 > 0$ and $R_1 > 0$. Applying the same argument, we obtain $G_{n}, K_n\geq0$ and so $S_{n+1}, I_{n+1}, R_{n+1} > 0$ for all $n$.

In this work, we are interested in the extinction and permanence of disease of the model (1.8), so the asymptotic behaviors of disease-free states of the system is crucial. By definition, a disease-free state of (1.8) is a solution $(S_n, I_n, R_n)$ where $I_n = 0$ for all $n\geq0$. From Lemma 2, the system is reduced in this case to

Let $N > n_1\geq n_0$. By (H2), (H3), and Proposition 3, we have that

as $N\to \infty$. Similarly, $R_{N}\to0$. Thus the asymptotic behaviors of the disease-free states of (1.8) are closely related to the properties of solutions $S_n^0$ of the following equation

The following proposition gives the uniform upper and lower bounds for any positive solution of (3.1).

Proposition 4. Suppose that (H3) and (H4) hold. Then there exist constants $S^{0, \star}, S^{0}_\star > 0$ such that

for any positive solution $S_n^0$ of (3.1) regardless of the initial condition. In fact, we have

for all $n_1\geq n_0$

Proof. Set $a(n) = h\Lambda(n)$ and $b(n) = \frac{1}{1+\mu_1(n)}$. It is obvious that $a^u < \infty$ and $b^l > 0$. Also, the hypotheses (H3) and (H4) are equivalent to (2.2) and (2.7) in Lemma 1, respectively. Hence, the proof immediately follows from Lemmas 1(1) and (2).

The following proposition is a generalized version of the preceding one. It will be used in the proof of the permanence of disease (Theorem 8).

Proposition 5. Suppose that (H3) and (H4) hold. Let $\nu\geq0$ and define $\mu_1^\nu$ by

Then there exist constants $S^{\nu, \star}, S^{\nu}_\star > 0$ such that if $S^\nu_n$ is a positive solution of the equation

then

Proof. As in the proof of Proposition 4, choosing $a(n) = h\Lambda(n)$ and $b(n) = \frac{1}{1+\mu^\nu_1(n)}$, and applying (H3) and (H4), we get the conditions (2.2) and (2.7). By Lemmas 1(1) and (2), the result is true and the quantities $S^{\nu, \star}, S^{\nu}_\star$ are given by replacing $\mu_1$ in the formulas (3.2) and (3.3) with $\mu_1^\nu$.

Notice that

The following result shows that the "total population" of the model (1.8) is bounded above.

Proposition 6. Assume $f(I)\geq0$ for all $I\geq0$ and that (H2), (H3) hold. Then the total population $T_n: = S_n+I_n+R_n$ for (1.8) satisfies

where $S^{0, \star}$ is given by (3.2).

Proof. Adding the equations in (1.8), applying the hypothesis (H2), and Proposition 3, we get

which implies

The desired conclusion now follows from Lemma 1(3) and (H3).

As a corollary, we obtain the boundedness of solutions $(S_n, I_n, R_n)$ of (1.8).

Corollary 7. Assume $f(I)\geq0$ for all $I\geq0$ and that (H2), (H3) hold. Then

Proof. Since $f\geq0$, we have by Proposition 3 that $S_n > 0, I_n > 0, R_n > 0$ for all $n > 0$. Then, we have $S_n\leq T_n, I_n\leq T_n$, and $R_n\leq T_n$, hence the desired conclusion follows from Proposition 6.

4.   Permanence of disease

Now we prove our first main result.

Theorem 8. Suppose that (H1)-(H4) hold, and that

where $S^{0}_\star$ is given by (3.3), $\beta^l = \liminf_{n\to\infty}\beta(n)$, and $\mu_2^u = \limsup_{n\to\infty}\mu_2(n)$. Then the model (1.8) is permanent. Moreover, we will show that if

where $K^u = \limsup_{n\to \infty}K_n\leq k_1$, then (1.8) is permanent.

Proof. It suffices to prove the second statement because $\frac{\beta^lf(0)S^{0}_\star}{\mu_2^u+K^u}\geq \mathcal{R}_\infty$. Our proof is inspired by ^[18] and ^[14]. By the positivity (Proposition 3) and boundedness of solutions (Corollary 7), it suffices to prove uniform lower bounds (2.1) for $S_n, I_n,$ and $R_n$ of (1.8).

By the assumption (4.2) and (H1), there exist $\theta_0 > 0$ and $\varepsilon_0 > 0$ sufficiently small such that

This also implies $\beta^l-\theta_0 > 0$ and $S_\star^0-\theta_0 > 0$. We also assume $K^l-\theta_0 > 0$.

Choose $Q_0\in \mathbb{N}$ independent of solutions to (1.8) such that

Estimate for $S_n$. There is a constant $l_S > 0$ such that $\liminf\limits_{n\to\infty}S_n\ge l_S$.

Proof of Estimate for $S_n$. By the second estimate in Corollary 7, we can take $P_0\in \mathbb{N}$ which depend on solution $(S_n, I_n, R_n)$ so that

Let $n > n_1\geq n_0+\max\{P_0, Q_0\}$. By Proposition 3, we have $R_{n+1}\geq0$. Also, $\beta(n)\sum_{j = 0}^mp_jf(I_{n-j})I_{n-j}\leq(\beta^u+\theta_0)f^M(S^{0, \star}+\theta_0) < \infty$. Using the first equation in (1.8), we get

that is

where $\tilde{\mu}_1(n) = \mu_1(n)+(\beta^u+\theta_0)f^M(S^{0, \star}+\theta_0)$. Applying Lemma 1(4) and Proposition 5, we find that $\liminf_{n\to\infty}S_n\geq \tilde S^{0}_\star > 0$ where

Estimate for $I_n$. There is a constant $l_I > 0$ such that $\liminf\limits_{n\to\infty}I_n\ge l_I$.

Proof of Estimate for $I_n$. We begin by proving a continuity result for the sum in (3.5) as $\nu\to0$.

Claim 1. Let $\theta_0 > 0$ be as specified in (4.3). Then there exist $\nu_0, n_0', \rho$, depending only on parameters in (1.8) and $\theta_0$ and satisfying

such that the following statement holds:

where $\mu_1^{\nu_0}: = \mu_1+\beta^Mf^M\nu_0$ and $S^0_\star$ is given by (3.3).

Proof of Claim 1. First we specify $\nu_0$. By (3.6), we have

and the limit does not depend on $n_1$. We apply Lemma 1(5) with $a(n) = h\varLambda(n), b(n) = 1/(1+h\mu_1(n))$, and $\tilde{b}(n) = 1/(1+h\mu_1^\nu(n))$ and $\nu > 0$. Note that $\tilde{b}(n)\geq\frac{1}{1+h(\mu_1^M+\beta^Mf^M\nu)} > 0$, where $\mu_1^M = \max_n\mu_1(n)$, and $\tilde{b}(n)\leq b(n)$ for all $n$, so $\tilde{b}^l > 0$ and $\tilde{b}$ satisfies (2.2). Also, $|b(n)- \tilde{b}(n)|\leq h\beta^Mf^M\nu$, so applying Lemma 1(5), there is a constant $C > 0$ independent of $\nu$ such that

Now we choose $\nu_0$ small enough so that $\nu_0\leq\varepsilon_0$ and

Then we have

Now we choose $n_0', \rho$. By Lemma 1(6), there exist $n_0', \rho$ depending only on the $\varLambda, \mu_1, f, m$, and $\theta_0$, so that

Combining the above two estimates, we obtain the desired result.

Fix $\theta_0 > 0$ and $\nu_0, n_0', \rho$ as above. Notice that $n_0'\geq n_0+Q_0$, so according to (4.4) we get

Next, we prove that if $I_n\leq\nu_0$ on an interval of length at least $\rho(m+1)$, then $S_{N+1} > S^0_\star-\theta_0$, where $N$ is the right endpoint.

Claim 2. Let $n_1\geq n_0'$, $N\geq n_1+\rho m$, and $(S_n, I_n, R_n)$ be a solution of (1.8). Then the following statement holds. If $0\leq I_n\leq\nu_0$ for all $n\in[n_1-m, N]$, then

Proof of Claim 2. Since $I_p\leq\nu_0$ on $[n_1-m, N]$, we have for all $n\in[n_1, N]$ that

This gives using (1.8) that $S_{n+1}\geq h(\varLambda(n)-\beta^Mf^M\nu_0 S_{n+1}-\mu_1(n)S_{n+1})+S_{n}$, i.e.

By induction, it is directly to see that

by the positivity of $S_{n_1}$. Applying Claim 1, the last term on the right hand side is greater than $S^{0}_\star-\theta_0$, which implies $S_{N+1} > S^0_\star-\theta_0$. So the claim is proved.

In addition to the previous estimate for $S_{N+1}$, at a right endpoint $N$, we also have the following boosting estimate for $I_{N+1}$.

Claim 3. Let $n_1\geq n_0'$, $N\geq n_1+\rho m$, and $(S_n, I_n, R_n)$ be a solution of (1.8). Then the following statement holds. If $0\leq I_n\leq\nu_0$ for all $n\in[n_1-m, N]$, then

and $\kappa > 1$ is the constant given by

Proof of Claim 3. Observe that $\kappa > 1$ because $\xi > 1$. We can apply Claim 2 to get $S_{N+1} > S^0_\star-\theta_0$ and by (4.7) with that $N\geq n_0'$, we also have $\beta(N)\geq\beta^l-\theta_0$. Using (4.3) and that $\nu_0\leq\varepsilon_0$, we obtain

for $j = 0, 1, \ldots, m$. Also, noting that $\mu_2(N)+K_N\leq\mu_2^u+K^u+\theta_0$ by (4.7). Applying Lemma 2, we have

where we have used that $\sum_{j = 0}^mp_j = 1$.

Claim 4. It is impossible that $I_n \leq\nu_0$ for all sufficiently large $n$.

Proof of Claim 4. Suppose on the contrary that there is $N_0\geq n_0'$ such that $I_n\leq\nu_0$ for all $n\geq N_0$. Denote $N_1 = N_0+(\rho+1)m$, and define $\mathcal{I}_n$ as in (4.9) above. Since $I_n\leq\nu_0$ for all $n\in[N_0, N_1]$, we have by Claim 3 that

This implies in particular that $\mathcal{I}_{N_1+1} = \min_{n\in[N_1+1-m, N_1+1]}I_n \geq \mathcal{I}_{N_1}$ because $\kappa > 1$. Repeating the preceding argument with that $\mathcal{I}_{N_1+1}\geq \mathcal{I}_{N_1}$, we get $I_{N_1+2}\geq\kappa \mathcal{I}_{N_1+1}\geq\kappa \mathcal{I}_{N_1}$. Continuing the argument, we finally obtain

This implies

Let $N_2 = N_1+1+\rho m$. Since $I_n\leq\nu_0$ for all $n\in[N_1+1-m, N_2]$ and observe that $\mathcal{I}_{N_2} \geq\kappa \mathcal{I}_{N_1}$, it follows by the same argument as above that $I_n\geq\kappa \mathcal{I}_{N_2}\geq\kappa^2 \mathcal{I}_{N_1}$ for all $n\geq N_2+1$. Repeating the process, we obtain, for $N_{l}: = N_{l-1}+1+\rho m$, that

Recalling $\kappa > 1$, so $\kappa^l\to\infty$ as $l\to\infty$. Now by selecting $l$ large enough, we get from the above that there is $n > N_0$ such that $I_n\geq\kappa^{l} \mathcal{I}_{N_1} > \nu_0$. This contradicts that $I_n\leq\nu_0$ for all $n\geq N_0$. Therefore the claim is proved.

After the above preparations, now we can prove the lower estimate for $I_n$. According to Claim 4, there are two possibilities: either

(i) $I_n > \nu_0$ for all $n$ sufficiently large, or

(ii) $I_n$ oscillates about $\nu_0$ for large $n$.

Obviously, if (i) occurs then the desired estimate follows, namely

So assume (ii).

It suffices to prove the following statement. Suppose $n_1, N$ are such that $N > n_1\geq n_0'$ and

Then

Denote

Let $n > n_1$, so $n-1\geq n_1\geq n_0'$. By (4.7), we have $\mu_2(n-1)+K_{n-1}\leq\mu_2^u+K^u+\theta_0$. Employing the second identity in Lemma 2 and the positivity of $G_{n-1}, S_n$, we get

By induction, we have for all $p\in \mathbb{N}$ satisfying $n-p\geq n_1$ that

Case $N\leq n_2$. For any $n\in(n_1, N)$, it follows by taking $p = n-n_1$, which satisfies $p\leq\Delta$, that

where we have employed the fact that $I_{n_1}\geq\nu_0$. Thus in this case (4.11) is true.

Case $n_2 < N\leq n_3$. We use the preceding case on $(n_1, n_2)$ to conclude that $I_n\geq l_I$ for all $n\in (n_1, n_2]$. This implies $\mathcal{I}_{n_2}\geq l_I$. We show that $I_n\geq l_I$ for $n\in (n_2, N]$ as well. Since $I_n\leq\nu_0$ on $[n_1+1, n_2]$ and $n_2\geq(n_1+1)+(\rho+1)m$, it follows by Claim 3 that

If $n_2+1 = N$ then we are done. Otherwise, we employ that $I_n\leq\nu_0$ on $[n_1+2, n_2+1]$ to continue the preceding argument and get

By induction, we finally conclude that $I_n\geq l_I$ for all $n\in(n_2, N]$ as desired.

Cases that $n_3 < N\leq n_4$ and so on can be proved by the same argument and is omitted.

Now (4.11) is proved hence establishing the desired estimate for $I_n$.

Estimate for $R_n$. There exists a positive constant $l_R$ such that $\liminf_{n\rightarrow \infty}R_n\ge l_R$.

Proof of Estimate for $R_n$. First observe that, from the estimate for $I_n$, we get that

where $l_I$ is given by (4.11). Choose $Q_0' > 0$ so that

Note that $K^l-\theta_0 > 0$ according to (4.3). Let $n\geq n_0'+Q_0'$. Using the third equation of Lemma 2, we get

We can apply parts (2) and (4) of Lemma 1 to find that there is a constant $l_R > 0$, independent of solutions to (1.8) such that

This establish the estimate for $R_n$, therefore it completes the proof of Theorem 8.

5.   Extinction of disease

We prove that under a certain condition the disease of the model (1.8) always extincts.

Theorem 9. Suppose that (H1)-(H4) hold and that

where $S^{0, \star}$ is given by (3.2), $\beta^u = \limsup\limits_{n\to\infty}\beta(n)$, and $\mu_2^l = \liminf\limits_{n\to\infty}\mu_2(n)$. Then (1.8) exhibits the extinction of the disease. Moreover, we will show that if

where $K^l: = \liminf_{n\to\infty}K_n\geq k_0$, then (1.8) exhibits the extinction of the disease.

Proof of Theorem 9. It suffices to prove the second statement because $\frac{\beta^u(\sup_{I\in[0, S^{0, \star}]}f(I))S^{0, \star}}{\mu_2^l+K^l}\leq \mathcal{R}_0$. We shall repeatedly employ the results on positivity (Proposition 3) and boundedness of solutions (Proposition 6 and Corollary 7).

By the assumption (5.2) and the continuity of $f$ (H1), we choose $\theta_0$ and $\varepsilon_0 > 0$ such that

Warning! Here and below, for simplicity, the parameters $\theta_0, \varepsilon_0, \xi, n_0', \kappa, \ldots$ have been reused and are in no connection to corresponding ones appeared in the proof of Theorem 8.

Let $(S_n, I_n, R_n)$ be a solution of (1.8). We also choose $n_0' > 0$ such that if $n\geq n_0'$ then

and

Now let $n\geq n_0'+m$. We observe that

So $G_n = \beta(n)\sum\limits_{j = 0}^mp_jf(I_{n-j})I_{n-j}\leq(\beta^u+\theta_0)(\sup_{I\in[0, S^{0, \star}+\varepsilon_0]}f(I)) \tilde{ \mathcal{I}}_n$, where

Applying the second equation in (2.9) and that $S_{n+1}\leq S^{0, \star}+\theta_0$, we get

In other words, we now obtain

where $0 < \kappa < 1$ is the constant given by

Employing (5.5), we are going to show that

for all $n\geq n_0'+m$.

By (5.5) we have $I_{n+1}\leq\kappa \tilde{ \mathcal{I}}_n$, so $\tilde{ \mathcal{I}}_{n+1} = \max_{p\in[n+1-m, n+1]}I_p\leq \tilde{ \mathcal{I}}_n$. Then by (5.5) again $I_{n+2}\leq\kappa \tilde{ \mathcal{I}}_{n+1}$, hence $I_{n+2}\leq\kappa \tilde{ \mathcal{I}}_n$. Similarly, (5.5) gives $I_{n+3}\leq\kappa \tilde{ \mathcal{I}}_{n+2}$ and $\tilde{ \mathcal{I}}_{n+2} = \max_{[n+2-m, n+2]}I_p\leq \tilde{ \mathcal{I}}_n$, so $I_{n+3}\leq\kappa \tilde{ \mathcal{I}}_n$. By induction, we get $I_{n+p}\leq\kappa \tilde{ \mathcal{I}}_n$ for all $p = 1, 2, \ldots$. Now

proving (5.7).

We apply the conclusion from (5.7). Setting

we directly get

Since $\kappa < 1$, the sequence $\{ \tilde{ \mathcal{I}}_{N_j}\}$ is monotonically decreasing. It is easy to see that $\tilde{ \mathcal{I}}_{N_j}\leq\kappa^j \tilde{ \mathcal{I}}_{N_0}$, hence

Now we prove the extinction of the disease, i.e. $\lim_{n\to\infty}I_n = 0$.

Let $\varepsilon > 0$. By (5.8) we can find $j$ such that $\tilde{ \mathcal{I}}_{N_{j}}\leq\varepsilon$. Now let $n\geq N_{j}$. We have $n\in[N_{r-1}, N_r)$ for some $r\geq j+1$. Then $I_n\leq \tilde{ \mathcal{I}}_{N_r}\leq \tilde{ \mathcal{I}}_{N_j}\leq\varepsilon$. Therefore $I_n\to0$ as desired.

6.   Numerical simulations

In this section, the model (1.8) is considered with non-monotone incidence rate

For simplicity, some parameters are fixed as follows: $h = 1$, $m = 2$, $p_0 = 0.2$, $p_1 = 0.3$, $p_2 = 0.5$, $\alpha_3 = 2$, $\alpha_1 = 1$, $\{\varLambda(n)\} = (1, 1.2, 1, 1.2, \dots)$, $\{\gamma(n)\} = (0.3, 0.3, 0.3, \dots)$, $\{\mu_2(n)\} = (0.5, 0.8, 0.5, 0.8, \dots)$, $\{\mu_1(n)\} = \{\mu_3(n)\} = (0.2, 0.4, 0.2, 0.4, \dots)$, $\{\beta(n)\} = (1, 1, 1, \ldots)$. We also impose the following initial conditions: $S_j = 3$, $I_j = 1$, $R_j = 0.5$ ($j = 0, -1, -2$). From the formulas (3.2) and (3.3), we have $S^{0, \star} = S^{0}_\star = 3.2353$.

Now, we present the examples and numerical simulations for different $\alpha_2, k_0, k_1$.

Example 1. Choose $\alpha_2 = 0.3$. We have $f(I) = \frac{g(I)}{I} = \frac{2}{1+0.3I+I^2}$. Since $f$ is a decreasing function, we have $f(0) = \sup_{I\in[0, S^{0, \star}]}f(I) = 2$. Also $\beta^u = \beta^l = 1$. We use $K_n = k_0+(k_1-k_0)\frac{10}{I_n+10}$.

(a) If $k_0 = 8$ and $k_1 = 9$, then we get $K^u = K^l = 9$ and $\frac{\beta^u(\sup_{I\in[0, S^{0, \star}]}f(I))S^{0, \star}}{\mu_2^l+K^l} = 0.6811$, $\mathcal{R}_0 = 0.7612 < 1$. Figure 1(a) indicates that the disease exhibits extinction.

(b) If $k_0 = 1$ and $k_1 = 3$, we get $K^u = 2.9556$, $K^l = 2.9525$, and $\frac{\beta^lf(0)S^{0}_\star}{\mu_2^u+K^u} = 1.7229$, and $\mathcal{R}_\infty = 1.7029 > 1$. Figure 1(b) indicates that the disease is permanent.

Example 2. Choose $\alpha_2 = -1$. We have $f(I) = \frac{g(I)}{I} = \frac{2}{1-I+I^2}$. It is directly to calculate that $f(0) = 2$ and $\sup_{I\in[0, S^{0, \star}]}f(I) = 2.6667$. Again $\beta^u = \beta^l = 1$. We use $K_n = k_0+(k_1-k_0)\frac{10}{I_n+10}$.

(a) If $k_0 = 9$ and $k_1 = 10$, then we get $K^u = K^l = 10$, $\frac{\beta^u(\sup_{I\in[0, S^{0, ~~ \star}~~ ]}~~~ f(I))S^{0, \star}}{\mu_2^l+K^l} = 0.8217$, and $\mathcal{R}_0 = 0.9082 < 1$. Figure 2(a) indicates that the disease exhibits extinction.

(b) If $k_0 = 2$ and $k_1 = 4$, we get $K^u = 3.9618, K^l = 3.9595$, $\frac{\beta^lf(0)S^{0}_\star}{\mu_2^u+K^u} = 1.3589$, and $\mathcal{R}_\infty = 1.3480 > 1$. Figure 2(b) indicates that the disease is permanent.

7.   Conclusions

In this paper, we investigate the discrete-time non-autonomous SIRS epidemic model (1.8), which is a discretization by the nonstandard finite difference method of the continuous-time model (1.4). In the model, a general nonlinear incidence rate with distributed delays is included together with a nonlinear recovery rate which takes into account the effect of health care resources such as the hospital beds (1.3). Two threshold parameters $\mathcal{R}_0$ and $\mathcal{R}_\infty$ are obtained so that if $\mathcal{R}_0 < 1$ then the disease dies out while if $\mathcal{R}_\infty > 1$ then the disease is permanent.

In the special case of autonomous (1.8), $f$ is a decreasing function, and $k$ a constant, our results imply that $\mathcal{R}_0 = \mathcal{R}_\infty = \frac{\beta f(0)\varLambda}{\mu_1(\mu_2+k)}$, which gives the basic reproduction number of the model. The same number is reported in Corollary 5.4 of ^[27] when $f = 1$. This result is also in line with the known results for the continuous-time models ^[12] when $f = 1$, and ^[3] when $f$ is decreasing and $f(0) = 1$.

For the non-autonomous model (1.8), we obtain the following corollary. Suppose that, as $n\to\infty$, we have $\beta(n)\to\beta'$, $\varLambda(n)\to\varLambda'$, $\mu_1(n)\to\mu_1'$, $\mu_2(n)\to\mu_2'$, $k(n)\to k'$ ($k$ does not depend on $I_n$), and $f = 1$. Then we obtain the reproduction number of the model to be $\mathcal{R}_0 = \mathcal{R}_\infty = \frac{\beta'S'}{\mu_2'+k'}$, which is the same number implied by Theorem 4.1 and Theorem 5.1 in ^[27]. Comparing this result to that from the continuous-time model in ^[25,26], we find an improvement because the threshold parameters obtained in those two papers do not lead to the reproduction number of the model.

A novelty of this paper is that our results are mathematically more tractable (compared e.g. to ^[27]) so they can be effective tools in practice to help the policy-makers fight the spread of the disease. To employ the condition (5.2), for example, one can first set the "ultimate" contact rate $\beta^\infty$ and health care resources $K^\infty$ so that the condition is met with $\beta^u, K^l$ replaced by $\beta^\infty, K^\infty$ (assuming all other parameters are available and fixed). Then set a certain starting time $N_0$ and control the transmission rate and the hospital beds to satisfy $\beta(n)\leq\beta^\infty$ and $K_n\geq K^\infty$ for all $n\geq N_0$ onward. It then follows from Theorem 9 that the spreading of the disease will be suppressed eventually.

For further investigations, it is interesting to extend the results of this paper to the model (1.8) where the transmission function $f$ depends not only on $I$ but also on $S$, such as the Beddington-DeAngelis function $f(S, I) = S/(1+m_1S+m_2I)$ and the saturated incidence $f(S, I) = S/(1+m_1S)(1+m_2I)$. It is also interesting to apply the technique in this paper to explore the threshold dynamics for other non-autonomous models such as a model with vaccination, a model with age structures, multi-strain diseases, etc.

Another interesting question is the chaotic dynamics of epidemic models with seasonality in the transmission rate ^[1,9]. It was shown in ^[1] for the classical SIR model that the disease dies out when $R_0 < 1$, where $R_0$ is the reproduction number, while if $R_0 > 1$ the model admits periodic and aperiodic patterns together with sensitive dependence on the initial conditions of the solution. The SIR model with logistic growth rate was explored in ^[9] and it was shown that the condition $R_0 < 1$ is not sufficient to guarantee the extinction of the disease due to backward bifurcation and the model exhibits persistent strange attractors. For our model (1.8), the condition $\mathcal{R}_0 < 1$ always implies the elimination of the disease regardless of the initial infectives. Moreover, for the periodic forced model with $f = 1$, it was shown in Corollary 5.3 ^[27] that if $R_0 < 1$ the disease dies out while it is permanent when $R_0 > 1$, where $R_0$ is given explicitly by $R_0 = \prod_{k = 0}^{\omega-1}(\frac{1+\beta(k)z_k^\ast}{1+\mu(k)+\gamma(k)+\alpha(k)})^{1/\omega}$ and $z_k^\ast$ is a unique $\omega$-periodic solution of (3.1). In our work, $f$ is a nonlinear function, so it is an interesting open problem to see whether a chaotic behavior can happen when $\mathcal{R}_0 > 1$ due to the seasonally forced contact rate $\beta$.

Appendix

We briefly discuss the basic idea of the nonstandard finite difference (NSFD) method that enables to get the discrete model (1.8) from the continuous model (1.4). NSFD is a discretization technique consisting of some rules with the aim to preserve significant properties of the related continuous problem (such as positivity, boundedness, stability, etc.) and avoid numerical instabilities. The following rules were proposed by Mickens ^[10] for constructing a NSFD scheme from a continuous problem:

Rule 1. The orders of the discrete derivatives of the scheme should be equal to the orders of the corresponding derivatives of the differential equation.

Rule 2. The denominator function for each discrete derivative should, in general, be expressed as a function of step-size which is more complicated than the conventional one. This rule is not strict and in this work, to get (1.8), we choose $\phi(h) = h$.

Example. Consider a first-order differential equation of the form

Conventionally, the discrete derivative for $dx/dt$ in (8.1) is given by

but, in the NSFD scheme, the denominator $\Delta t$ is replaced by a denominator function $\phi(h)$ with the property

Rule 3. Nonlinear terms should, in general, be replaced by nonlocal discrete representations using more than one mesh point. For example, the nonlinear term $x^2$ can be replaced by a nonlocal representation evaluated at two mesh points such as $x_{k+1}x_k$ or $2x^2_{k}-x_{k+1}x_{k}$.

To get (1.8), the incidence rate (1.7) is discretized by $\beta(n)\sum_{j = 0}^mp_jf(I_{n-j})I_{n-j}S_{n+1}$.

Rule 4. Special conditions that hold for the solutions of the differential equations should also hold for the solutions of the finite difference scheme. An important example is the positivity of solutions when the solutions represent some physical positive quantities. If the discrete equations allow their solutions to become negative, then numerical instabilities will occur.

We have shown the basic properties including Propositions 3 and 6.

Rule 5. The finite difference scheme should not introduce extraneous or spurious solutions that do not correspond to any solution of the corresponding differential equations.

Conflict of interest

The authors declare that there are no conflicts of interest.

Acknowledgments

The authors would like to thank the reviewers for valuable comments and suggestions which are tremendously helpful to improve the manuscript. We are thankful for financial support by Science Achievement Scholarship of Thailand (SAST). We are also supported by the 90th Anniversary of Chulalongkorn University Fund, Thailand (Ratchadaphiseksomphot Endowment Fund).