Global Dynamics of a Discrete-Time MERS-Cov Model

Abstract

In this paper, we have investigated the global dynamics of a discrete-time middle east respiratory syndrome (MERS-Cov) model. The proposed discrete model was analyzed and the threshold conditions for the global attractivity of the disease-free equilibrium (DFE) and the endemic equilibrium are established. We proved that the DFE is globally asymptotically stable when ${\mathcal{R}}_0\le 1$. Whenever ${\tilde{\mathcal{R}}}_0>1$, the proposed model has a unique endemic equilibrium that is globally asymptotically stable. The theoretical results are illustrated by a numerical simulation.

1. Introduction

Mankind has been suffering from several epidemic diseases during every age of human being history. In the past century, many disease outbreaks have occurred, such as influenza, malaria, dengue fever, SARS, H1N1, H7N1, pestilence, AIDS, MERS and Covid-19. Recently, most of the viruses that affect the mankind life are the coronaviruses, which are a single standard RNA virus that were firstly reported in 1960s [1].

Severe acute respiratory syndrome (SARS) is an airborne virus which was firstly found in Foshan, China on 16 November 2002 [2]. The spread of the SARS virus was through small droplets of saliva, similar to influenza with an incubation period of 2–7 days. The SARS virus transmission was from bats to humans. The number of reported cases was 8098 and the number of death cases was 774 [3]. The middle east respiratory syndrome (MERS) is a viral respiratory disease caused by novel coronavirus MERS-Cov, that was first reported in Saudi Arabia in 2012 [4]. The infection of MERS may occur through both animals, such as camels and bats, and humans by a direct or indirect contact [5]. The incubation period of MERS is 2 to 14 days [6]. Since 2012, the reported infected cases of MERS are 2494 with a death rate of 34% [7]. In early December 2019, the first case of Covid-19 was reported in Wuhan, China [8]. The symptoms of Covid-19 are similar to the other coronaviruses SARS and MERS, such as colds and high temperature, which lead to pneumonia and hence respiratory system failure and death [9].

A vast number of mathematical models have been proposed to understand the dynamics of the coronaviruses. The SIS, SIR and SEIR models are studied by many researchers to investigate the dynamics of the coronaviruses [10,11,12,13,14,15,16,17,18,19,20,21,21]. The model we have proposed in this paper is investigating the dynamics of the middle east respiratory syndrome coronavirus. We have studied the asymptomatic, symptomatic and hospitalized individuals’ effects on the spread of the virus.

Yong and Owen [22] have discussed the dynamical transmission model of MERS-Cov in two areas. They proved that the model has two equilibrium points, a disease-free equilibrium point and an endemic equilibrium point. Moreover, the disease dies out whenever the basic reproductive number is less than unity. Usaini et al. [17] proposed a new deterministic mathematical model for the transmission dynamics of MERS-Cov. They found that there is a unique endemic equilibrium point and there is no infection free equilibrium point due to the constant influx of latent immigrants. Lee et al. [16] have designed a dynamic transmission model to analyze the MERS-Cov outbreak in the Republic of Korea. This model incorporates the time-dependent parameters and the pulse of infections. Moreover, they estimated the basic reproductive number, ${\mathcal{R}}_e$, and showed that it is decreased which indicates that the MERS-Cov outbreak in the Republic of Korea had a low transmissibility.

Scientists have studied several types of dynamical epidemic models such as continuous-time type models that are described by differential equations and discrete-time type models that are described by difference equations. Currently, the scientists have paid more attention to the investigation of bigdata, which have made the discrete-time type epidemic models more interesting. In addition, discrete-time models have more dynamical behaviors [23].

Numerous studies have been done to explore the discrete-time epidemic models. L. Wang et al. [23] have studied a class of discrete SIRS epidemic models with disease courses. The authors have computed the basic reproduction number ${\mathcal{R}}_0$. In addition, they have proved that the disease-free equilibrium is globally attractive when ${\mathcal{R}}_0<1$ and the disease is permanent whenever ${\mathcal{R}}_0>1$. Y. Wang et al. [24] have introduced Lyapunov functions for a class of discrete SIRS epidemic models with nonlinear incidence rate and varying population size. The authors have established the sufficient and necessary conditions on the global asymptotic stability of the disease-free equilibrium and endemic equilibrium with general nonlinear incidence rate $\beta Sg\left(I\right)$ and different death rates. X. Fan et al. [25] have investigated a class of SEIRS epidemic models with a general nonlinear incidence function. In addition, they considered a discrete SEIRS model with standard incidence function. Moreover, the authors have shown that the model has a disease-free equilibrium, is globally attractive when the basic reproduction number ${\mathcal{R}}_0\le 1$ and the disease is permanent whenever ${\mathcal{R}}_0>1$. Batarfi et al. [18] have proposed a nonlinear mathematical model for MERS-Cov with two discrete-time delays. They computed the reproduction number, ${\mathcal{R}}_0$, and proved that there exists a disease-free equilibrium point when ${\mathcal{R}}_0\le 1$ and there is an endemic equilibrium point whenever ${\mathcal{R}}_0>1$. M. Khan et al. [26] have investigated a discrete-time TB model which is parameterized by the cases in the Pakistani of Khyber Pakhtunkhwa between 2002 and 2017. The authors have computed the reproduction number ${\mathcal{R}}_0$ which showed that the discrete-time TB model is stable at the disease-free equilibrium point when ${\mathcal{R}}_0<1$ and the model is globally asymptotically stable for the endemic equilibrium point whenever ${\mathcal{R}}_0>1$. Moreover, the authors have compared the discrete-time model with the continuous-time model. M. Safi et al. [27] have considered a discrete-time mathematical model that is obtained from the continuous-time model in [28]. The authors investigated the stability of the model and they proved that the model is globally asymptotically stable when ${\mathcal{R}}_0<1$.

This paper is organized as follows: A discrete MERS-Cov epidemic model is considered in Section 2. In Section 3, we present the fundamental properties of the discrete model. The stability analysis of the disease-free equilibrium is carried out in Section 4. The existence and stability analysis of endemic equilibrium point is conducted in Section 5. In Section 6, numerical simulations are provided to illustrate the obtained results and the results are concluded.

2. Model Formulation

The backward difference scheme is applied to propose the following discrete epidemic MERS-Cov coronavirus model. This model is called SEAIHR model. The population $N\left(n\right)$ is divided into the following compartments: $S\left(n\right)$ is the susceptible individuals, $E\left(n\right)$ is the exposed individuals, $A\left(n\right)$ is the asymptomatic individuals, $I\left(n\right)$ is the symptomatic individuals, $H\left(n\right)$ is the hospitalized individuals and $R\left(n\right)$ is the recovered individuals. Hence, $N\left(n\right)=S\left(n\right)+E\left(n\right)+A\left(n\right)+I\left(n\right)+H\left(n\right)+R\left(n\right)$. Susceptible people (S) is increased by the recruitment of individuals into the population, at a rate of $\mathsf{\Pi }$. This class is decreased by infection (with the rate of $\lambda$). Furthermore, this population is decreased by natural death (at a rate of $\mu$; populations in all classes are assumed to have the same natural death rate). Exposed individuals (E) are generated with the rate of $\lambda$ and reduced by progression to the asymptomatic individuals (A) at rate of $\sigma$ and to the symptomatic individuals (I) at rate of (r). The class of asymptomatic individuals (A) increased to the exposed individuals at a rate of $\sigma$ and it is reduced by progression to the exposed at rate of (r) and to the recovered individuals (R) at rate of ${\gamma }_1$. The symptomatic class (I) is increased by the exposed people at a rate of ($r\sigma$) and decreased to the hospitalized individuals (isolation) (H) at rate of $\varphi$, to the recovered people at rate of ${\gamma }_2$ and disease-induced death (at a rate ${\delta }_1$). The hospitalized (isolated) class is increased by the symptomatic individuals (I) at rate of $\varphi$ and decreased to recovery at rate of ${\gamma }_3$ and disease-induced death (at a rate ${\delta }_2$). Figure 1 illustrates the model (1) by a schematic diagram. Therefore, the SEAIHR model is governed by the following difference equations:

$$\left\{\begin{array}{cc}S(n+1)-S\left(n\right)\hfill & =\mathsf{\Pi }-\lambda (n+1)S(n+1)-\mu S(n+1),\hfill \\ E(n+1)-E\left(n\right)\hfill & =\lambda (n+1)S(n+1)-k_{1}E(n+1),\hfill \\ A(n+1)-A\left(n\right)\hfill & =(1-r)\sigma E(n+1)-k_{2}A(n+1),\hfill \\ I(n+1)-I\left(n\right)\hfill & =r\sigma E(n+1)-k_{3}I(n+1),\hfill \\ H(n+1)-H\left(n\right)\hfill & =\varphi I(n+1)-k_{4}H(n+1),\hfill \\ R(n+1)-R\left(n\right)\hfill & ={\gamma }_{1}A(n+1)+{\gamma }_{2}I(n+1)+{\gamma }_{3}H(n+1)-\mu R(n+1),\hfill \end{array}\right.$$

(1)

where $\lambda (n+1)={\displaystyle \frac{\beta [A(n+1)+{\eta }_{1}I(n+1)+{\eta }_{2}H(n+1)]}{N(n+1)}}$, $k_1=\sigma +\mu$, $k_2=\mu +{\gamma }_1$, $k_3=\mu +\varphi +{\delta }_1+{\gamma }_2$, $k_4=\mu +{\delta }_2+{\gamma }_3$, $\mathsf{\Pi }$ is the recruitment rate of susceptible people corresponding to births and immigration, $\mu$ is the natural death rate, $\beta$ is the contact rate, ${\displaystyle \frac{1}{\sigma }}$ is the mean time of incubation period, ${\displaystyle \frac{1}{\varphi }}$ is the mean time from data of symptoms onset to data of hospitalization, r is the clinical outbreak rate in all the infected cases, ${\displaystyle \frac{1}{{\gamma }_1}}$ is the mean infections period of asymptomatic infected person for survivors, ${\displaystyle \frac{1}{{\gamma }_2}}$ is the mean duration of infected person for survivors, ${\displaystyle \frac{1}{{\gamma }_3}}$ is the mean duration for hospitalized cases for survivors and ${\delta }_i(i=1,2)$ is the disease-induced death rate of infectious. Subject to the following initial conditions:

$$S\left(0\right)>0,E\left(0\right)>0,A\left(0\right)>0,I\left(0\right)>0,H\left(0\right)>0,R\left(0\right)>0.$$

(2)

It is worthy to mention that such models are potentially analytically solvable by using the Lie algebra method through matrix exponentials [29]. Based on the Kolmogorov equation and the Wei–Norman method, the analytical solution for the proposed model can be obtained in terms of matrix exponentials (for more details about Lie algebra method, see [30,31]).

3. Fundamental Properties

Lemma 1.

The model (1) with initial conditions (2) has a unique positive solution $\left(S\left(n\right),E\left(n\right),A\left(n\right),I\left(n\right),H\left(n\right),R\left(n\right)\right)$, $\forall n=1,2,3,\cdots$

To prove Lemma 1, we will apply the mathematical induction as follows:

Proof. 

Rewrite the model (1) as follows:

$$\left\{\begin{array}{cc}A(n+1)\hfill & =a_{1}E(n+1)+b_1\hfill \\ I(n+1)\hfill & =a_{2}E(n+1)+b_2\hfill \\ H(n+1)\hfill & =a_{3}E(n+1)+b_3\hfill \\ R(n+1)\hfill & =a_{4}E(n+1)+b_4\hfill \\ S(n+1)\hfill & =a_{5}E(n+1)+b_5\hfill \end{array}\right.$$

(3)

where

$$a_1=\frac{(1-r)\sigma }{1+k_2},b_1=\frac{A\left(n\right)}{1+k_2}$$

$$a_2=\frac{r\sigma }{(1+k_3)},b_2=\frac{I\left(n\right)}{(1+k_3)}$$

$$a_3=\frac{r\varphi \sigma }{(1+k_3)(1+k_4)},b_3=\frac{\varphi I\left(n\right)}{(1+k_3)(1+k_4)}+\frac{H\left(n\right)}{1+k_4}$$

$$a_4=\frac{\gamma a_1+{\gamma }_2{a}_2+{\gamma }_3{a}_3}{1+\mu },b_4=\frac{{\gamma }_1{b}_1+{\gamma }_2{b}_2+{\gamma }_3{b}_3+R\left(n\right)}{1+\mu }$$

$$a_5=-\frac{1+k_1}{1+\mu },b_5=\frac{\mathsf{\Pi }+S\left(n\right)+E\left(n\right)}{1+\mu }$$

Substituting the expression (3) in the second equation of model (1) yields -4.6cm0cm

$${\displaystyle E(n+1)=\frac{E\left(n\right)}{1+k_1}+\frac{\beta \left[a_{5}E(n+1)+b_5\right]\left[(a_1+{\eta }_1{a}_2+{\eta }_2{a}_3)E(n+1)+b_1+{\eta }_1{b}_2+{\eta }_2{b}_3\right]}{(1+k_1)[(1+{\sum }_{i=1}^5{a}_i)E(n+1)+{\sum }_{i=1}^5{b}_i]}}$$

(4)

Let $u=E(n+1)$. Hence, we define

$${\displaystyle f\left(u\right)=u-\frac{E\left(n\right)}{1+k_1}-\frac{\beta \left[a_{5}u+b_5\right]\left[(a_1+{\eta }_1{a}_2+{\eta }_2{a}_3)u+b_1+{\eta }_1{b}_2+{\eta }_2{b}_3\right]}{(1+k_1)[(1+{\sum }_{i=1}^5{a}_i)u+{\sum }_{i=1}^5{b}_i]}=0}$$

(5)

Let $A=1+{\sum }_{i=1}^5{a}_i$, $B={\sum }_{i=1}^5{b}_i$, $C=\beta a_5(a_1+{\eta }_1{a}_2+{\eta }_2{a}_3)$, $D=\beta a_5(b_1+{\eta }_1{b}_2+{\eta }_2{b}_3)$, $E=b_5(a_1+{\eta }_1{a}_2+{\eta }_2{a}_3)$ and $F=b_5(b_1+{\eta }_1{b}_2+{\eta }_2{b}_3)$. Then, Equation (5) can be written as follows:

$$f\left(u\right)=\frac{\left(A(1+k_1)C\right)u^2-(AE\left(n\right)+D+E)u-(BE\left(n\right)+F)}{(1+k_1)(Au+B)}$$

Which is simplified to

$$f\left(u\right)={\mathsf{\Gamma }}_{1}u+{\mathsf{\Gamma }}_2-\frac{{\mathsf{\Gamma }}_3}{u}$$

where ${\mathsf{\Gamma }}_1=1-\frac{C}{A(1+k_1)}$, ${\mathsf{\Gamma }}_2={\displaystyle \frac{BC-AB(1+k_1)-A^{2}E\left(n\right)-A(D+E)}{A^2(1+k_1)}}$ and ${\mathsf{\Gamma }}_3={\displaystyle \frac{A^{2}F-A{B}^2(1+k_1)+B^{2}C-AB(D+E)}{A^2}}$. Since ${\mathsf{\Gamma }}_1>0$, see Figure 2, then ${lim}_{u\to \infty }f\left(u\right)>0$. Moreover, $f^{\prime }\left(u\right)={\mathsf{\Gamma }}_1+{\displaystyle \frac{{\mathsf{\Gamma }}_3}{u^2}}$ and hence, $f^{\prime }\left(u\right)>0$, then $f\left(u\right)$ is an increasing function and ${\displaystyle f\left(0\right)=-\frac{\beta b_5(b_1+{\eta }_1{b}_2+{\eta }_2{b}_3)}{(1+k_1){\sum }_{i=1}^5{b}_i}-\frac{E\left(n\right)}{1+k_1}<0}$, then $f\left(u\right)=0$ has a unique positive solution and therefore, Equation (4) has a unique solution $E(n+1)>0$ which implies that a unique $A(n+1)$, $I(n+1)$, $H(n+1)$ and $R(n+1)$ exist with $A(n+1)>0$, $I(n+1)>0$, $H(n+1)>0$ and $R(n+1)>0$.

Let $w=S(n+1)$. Thus, from the first equation in model (1), we define

$$g\left(w\right)=(1+\mu )w+\frac{\beta \left[(a_1+{\eta }_1{a}_2+{\eta }_2{a}_3)E(n+1)+b_1+{\eta }_1{b}_2+{\eta }_2{b}_3\right]w}{w+M}-S\left(n\right)-\mathsf{\Pi }$$

where $M=E(n+1)+A(n+1)+I(n+1)+H(n+1)+R(n+1)$. Since $g\left(w\right)$ is continuous and $g^{\prime }\left(w\right)>0$, then $g\left(w\right)$ is an increasing function with $g\left(0\right)=-S\left(n\right)-\mathsf{\Pi }<0$. Thus, by the intermediate value theorem $g\left(w\right)=0$ has a unique positive solution and therefore, there exists a unique $S(n+1)>0$ which completes the proof. □

Lemma 2.

Any solution $\left(S\right(n),E(n),A(n),I(n),H(n),R(n\left)\right)$ of model (1) with initial conditions (2) satisfies ${\displaystyle \underset{n\to \infty }{lim\; sup}N\left(n\right)\le \frac{\mathsf{\Pi }}{\mu }}$.

Proof. 

Since,

$$\begin{array}{cc}\hfill N(n+1)& =S(n+1)+E(n+1)+A(n+1)+I(n+1)+H(n+1)+R(n+1)\hfill \\ \hfill & =\mathsf{\Pi }-\mu \left[S\right(n+1)+E(n+1)+A(n+1)+I(n+1)+H(n+1)\hfill \\ \hfill & +R(n+1)]-{\delta }_{1}I(n+1)-{\delta }_{2}H(n+1)+S\left(n\right)+E\left(n\right)+A\left(n\right)+I\left(n\right)\hfill \\ \hfill & +H\left(n\right)+R\left(n\right)\hfill \\ \hfill & =\mathsf{\Pi }-\mu \left[S\right(n+1)+E(n+1)+A(n+1)+I(n+1)+H(n+1)\hfill \\ \hfill & +R(n+1)]-{\delta }_{1}I(n+1)-{\delta }_{2}H(n+1)+N\left(n\right)\hfill \\ \hfill & =\mathsf{\Pi }-\mu N(n+1)-{\delta }_{1}I(n+1)-{\delta }_{2}H(n+1)+N\left(n\right)\hfill \\ \hfill N(n+1)+\mu N(n+1)& =\mathsf{\Pi }-{\delta }_{1}I(n+1)-{\delta }_{2}H(n+1)+N\left(n\right)\hfill \\ \hfill N(n+1)& =\frac{\mathsf{\Pi }+N\left(n\right)-{\delta }_{1}I(n+1)-{\delta }_{2}H(n+1)}{1+\mu }\hfill \\ \hfill & \le \frac{\mathsf{\Pi }+N\left(n\right)}{1+\mu },n=0,1,2,\cdots \hfill \\ \hfill & =\frac{\mathsf{\Pi }}{1+\mu }+\frac{N\left(n\right)}{1+\mu }.\mathrm{Upon}\mathrm{using}\mathrm{the}\mathrm{iteration}\mathrm{method}\mathrm{we}\mathrm{get}\hfill \\ \hfill N(n+1)& \le \frac{\mathsf{\Pi }}{1+\mu }+\frac{\mathsf{\Pi }}{{(1+\mu )}^2}+\frac{\mathsf{\Pi }}{{(1+\mu )}^3}+\cdots +\frac{\mathsf{\Pi }}{{(1+\mu )}^{(n+1)}}+\frac{N\left(0\right)}{{(1+\mu )}^{(n+1)}}\hfill \\ \hfill & =\frac{\mathsf{\Pi }}{\mu }\left[1-\frac{1}{{(1+\mu )}^{(n+1)}}\right]+\frac{N\left(0\right)}{{(1+\mu )}^{(n+1)}}.\hfill \end{array}$$

Hence,

$$\underset{n\to \infty }{lim\; sup}N(n+1)\le \underset{n\to \infty }{lim\; sup}N\left(n\right)\le \frac{\mathsf{\Pi }}{\mu }$$

□

Therefore, the region

$${\displaystyle \mathcal{D}=\left\{(S\left(n\right),E\left(n\right),A\left(n\right),I\left(n\right),H\left(n\right),R\left(n\right))\in {\mathcal{R}}_{+}^6|N\left(n\right)\le \frac{\mathsf{\Pi }}{\mu }\right\}}$$

is positively invariant.

4. Disease-Free Equilibrium (DFE)

4.1. Local Stability of DFE

The unique disease-free equilibrium (DFE) of the system (1) is given by

$${\mathcal{E}}_0=\left(\frac{\mathsf{\Pi }}{\mu },0,0,0,0,0\right)$$

To compute the basic reproduction number of model (1), we will apply the next generation operator method [32,33,34,35,36]. The matrix of the new infection terms, F, and the matrix of the transition terms, V, that are associated with the model (1) are given by:

$$F=\left[\begin{array}{cccc}0& \beta & \beta {\eta }_1& \beta {\eta }_2\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\mathrm{and}V=\left[\begin{array}{cccc}{k}_1& 0& 0& 0\\ -(1-r)\sigma & k_2& 0& 0\\ -r\sigma & 0& k_3& 0\\ 0& 0& -\varphi & k_4\end{array}\right]$$

Following [34], the basic reproduction number is denoted by ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)$ and is given by

$${\mathcal{R}}_0=\frac{\beta \sigma [(1-r)k_3{k}_4+r{\eta }_1{k}_2{k}_4+r\varphi {\eta }_2{k}_2]}{k_1{k}_2{k}_3{k}_4}$$

The proof of following lemma can be deduced from the proof of Theorem 2 in [34], which is the following “Consider the disease transmission model given by (1). If $x_0$ is a disease-free equilibrium of the model, then $x_0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$, but unstable if ${\mathcal{R}}_0>1$, where ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)$.”

Lemma 3.

The DFE point ${\displaystyle {\mathcal{E}}_0=\left(\frac{\mathsf{\Pi }}{\mu },0,0,0,0,0\right)}$ of model (1) is locally asymptotically stable (LAS) when ${\mathcal{R}}_0<1$ and unstable when ${\mathcal{R}}_0>1$.

4.2. Global Stability of DFE

In this section, the global attractivity of the disease-free equilibrium of model (1) is investigated and we can obtain the following result.

Theorem 1.

The DFE of the model (1) is globally-asymptotically stable (GAS) in $\mathcal{D}$ whenever ${\mathcal{R}}_0\le 1.$

Proof. 

Consider the following Lyapunov function

$$F_1\left(n\right)=p_1{k}_{1}E\left(n\right)+p_{2}A\left(n\right)+p_{3}I\left(n\right)+p_{4}H\left(n\right)$$

where

$$p_1=\frac{\sigma [(1-r)k_3{k}_4+r{\eta }_1{k}_2{k}_4+r\varphi {\eta }_2{k}_2]}{k_1{k}_2{k}_3},p_2=\frac{k_4}{k_2},p_3=\frac{{\eta }_1{k}_4+{\eta }_2\varphi }{k_3}\mathrm{and}{p}_4={\eta }_2$$

The backward difference of $F_1$ is denoted by $\Delta F_1$ and is given by

$$\begin{array}{cc}\hfill \Delta F_1& =F_1(n+1)-F_1\left(n\right)\hfill \\ \hfill & =p_1[E(n+1)-E\left(n\right)]+p_2[A(n+1)-A\left(n\right)]+p_3[I(n+1)-I\left(n\right)]\hfill \\ \hfill & +p_4[H(n+1)-H\left(n\right)]\hfill \\ \hfill & =p_1\left[\lambda (n+1)S(n+1)-k_{1}E(n+1)\right]+p_2[(1-r)\sigma E(n+1)\hfill \\ \hfill & -k_{2}A(n+1)]+p_3[r\sigma E(n+1)-k_{3}I(n+1)]+p4[\varphi I(n+1)\hfill \\ \hfill & -k_{4}H(n+1)]\hfill \\ \hfill \mathrm{Sin}\mathrm{ce}& \frac{S(n+1)}{N(n+1)}\le 1\mathrm{in}\mathcal{D},\mathrm{then}\hfill \\ \hfill F_1(n+1)-F_1\left(n\right)& \le p_1\left[\beta [A(n+1)+{\eta }_{1}I(n+1)+{\eta }_{2}H(n+1)]-k_{1}E(n+1)\right]\hfill \\ \hfill & +p_2[(1-r)\sigma E(n+1)-k_{2}A(n+1)]+p_3[r\sigma E(n+1)\hfill \\ \hfill & -k_{3}I(n+1)]+p4[\varphi I(n+1)-k_{4}H(n+1)]\hfill \\ \hfill & =k_4({\mathcal{R}}_0-1)[A(n+1)+{\eta }_{1}I(n+1)+{\eta }_{2}H(n+1)]\hfill \end{array}$$

This implies that $\Delta F_1=F_1(n+1)-F_1\left(n\right)\le 0\mathrm{whenever}{\mathcal{R}}_0\le 1$ and $\Delta F_1=0$ if and only if $E(n+1)=A(n+1)=I(n+1)=H(n+1)=0$. Hence, $(E,A,I,H)\to (0,0,0,0)$ as $n\to \infty$. Upon setting $E=A=I=H=0$ in the first and last equations in model (1), we get ${\displaystyle S\to \frac{\mathsf{\Pi }}{\mu }}$ and $R\to 0$ as $n\to \infty$. Thus, the maximum invariable set in $\{(S,E,A,I,H,R):F_1\left(n\right)=0\}$ is a disease-free equilibrium ${\mathcal{E}}_0$. Following the theorems of stability of difference equations (Theorem 6.3 in [36]), every solution of the equations in model (1) with the initial conditions in $\mathcal{D}$ approaches ${\mathcal{E}}_0$ as $n\to \infty$. Thus, the disease-free equilibrium ${\mathcal{E}}_0$ of model (1) is globally attractive. Hence, the proof is completed. □

It is worthy to remark that similar techniques have been used in the proof of stability with feedback (for more details see [37]).

5. Endemic Equilibria

5.1. Existence of the Endemic Equilibrium Point EEP

Let ${\mathcal{E}}_1=(S_1,E_1,A_1,I_1,H_1,R_1)$ be and endemic equilibrium point for the model (1). Hence, we can conclude the following lemma.

Lemma 4.

The model (1) has a unique endemic equilibrium point ${\mathcal{E}}_1\in {\mathbb{R}}_{+}^6$, whenever ${\mathcal{R}}_0>1$.

Proof. 

By solving the equations of the model (1) at steady-state, we get:

$$\begin{array}{cc}\hfill S_1& =\frac{\mathsf{\Pi }}{{\lambda }_1+\mu },E_1=\frac{{\lambda }_1{S}_1}{k_1},A_1=\frac{(1-r)\sigma {\lambda }_1{S}_1}{k_1{k}_2}\hfill \\ \hfill I_1& =\frac{r\sigma {\lambda }_1{S}_1}{k_1{k}_3},H_1=\frac{r\varphi \sigma {\lambda }_1{S}_1}{k_1{k}_3{k}_4}\hfill \\ \hfill R_1& =\frac{r\sigma {\lambda }_1{S}_1}{\mu k_1{k}_2{k}_3{k}_4}\left[(1-r){\gamma }_1{k}_3{k}_4+r{\gamma }_2{k}_2{k}_4+r{\gamma }_3\varphi k_2\right]\hfill \end{array}$$

(6)

where

$${\lambda }_1=\frac{\beta [A_1+{\eta }_1{I}_1+{\eta }_2{H}_1]}{N_1}\mathrm{and}{N}_1=S_1+E_1+A_1+I_1+H_1+R_1$$

Substitute Equation (6) in the expression of ${\lambda }_1$ to get

$${\mathcal{R}}_0=1+\mathsf{\Gamma }{\lambda }_1$$

where

$$\mathsf{\Gamma }=\frac{\mu k_2{k}_3{k}_4+(1-r)\sigma (\mu +{\gamma }_1)k_3{k}_4+r\sigma (\mu +{\gamma }_2)k_2{k}_4+r\varphi \sigma (\mu +{\gamma }_3)k_2}{\mu k_1{k}_2{k}_3{k}_4}$$

Since $\mathsf{\Gamma }>0$, then ${\displaystyle \frac{{\mathcal{R}}_0-1}{\mathsf{\Gamma }}>0}$ if and only if ${\mathcal{R}}_0>1$. This implies that $S_1,E_1,A_1,I_1,H_1,R_1>0$ if and only if ${\mathcal{R}}_0>1$. □

5.2. Stability of the Endemic Equilibrium

Although no global asymptotic stability result is given here for the endemic equilibrium point ${\mathcal{E}}_1$, extensive numerical simulations suggest that ${\mathcal{E}}_1$ is a GAS in $\mathcal{D}/{\mathcal{D}}_0$, whenever ${\mathcal{R}}_0>1$. Hence, we add the following conjecture.

Conjecture 1.

The unique endemic equilibrium point (EEP) of the model (1) is globally asymptotically stable (GAS) in $\mathcal{D}$ if ${\mathcal{R}}_0>1$.

For mathematical convenience, we provide the proof of the conjecture (1) for the special case when the associated disease-induced mortality is neglected, i.e, ${\delta }_1={\delta }_2=0$. Therefore, we define

$${\mathcal{D}}_0=\left\{(S,E,A,I,H,R)\in \mathcal{D}|E=A=I=H=R=0\right\}$$

which is the basin of attraction of the disease-free equilibrium point ${\mathcal{E}}_0$, to conduct the global stability of the unique endemic equilibrium point for this case. Hence, the following result is obtained.

Theorem 2.

The unique endemic equilibrium point (EEP) of the model (1) with ${\delta }_1={\delta }_2=0$ is globally asymptotically stable (GAS) in $\mathcal{D}/{\mathcal{D}}_0$ if

$${\tilde{\mathcal{R}}}_0={\mathcal{R}}_0{|}_{{\delta }_1={\delta }_2=0}>1$$

Proof. 

Set ${\delta }_1={\delta }_2=0$ in model (1) and consider the following Lyapunov function

$$\begin{array}{cc}\hfill F_2\left(n\right)& =\frac{1}{2}\left[(S\left(n\right)-S_1)+(E\left(n\right)-E_1)+(A\left(n\right)-A_1)+(I\left(n\right)-I_1)+(H\left(n\right)-H_1)\right.\hfill \\ \hfill & {\left.+(R\left(n\right)-R_1)\right]}^2\hfill \\ \hfill & =\frac{1}{2}{\left[N\left(n\right)-N_1\right]}^2\hfill \end{array}$$

Then, the backward difference of $F_2$ is given by

$$\begin{array}{cc}\hfill \Delta F_2& =F_2(n+1)-F_2\left(n\right)\hfill \\ \hfill & =\frac{1}{2}{\left[N(n+1)-N_1\right]}^2-\frac{1}{2}{\left[N\left(n\right)-N_1\right]}^2\hfill \\ \hfill & =\frac{1}{2}\left[N(n+1)-N\left(n\right)\right]\left[N(n+1)+N\left(n\right)-2N_1\right]\hfill \\ \hfill & =-\frac{1}{2}{\left[N(n+1)-N\left(n\right)\right]}^2+\left[N(n+1)-N_1\right]\left[N(n+1)-N\left(n\right)\right]\hfill \\ \hfill & \le \left[N(n+1)-N_1\right]\left[N(n+1)-N\left(n\right)\right].\hfill \end{array}$$

Upon setting ${\delta }_1={\delta }_2=0$ and adding the equations of model (1), we get

$N(n+1)-N\left(n\right)=\mathsf{\Pi }-\mu N(n+1)$ and ${\displaystyle N_1=\frac{\mathsf{\Pi }}{\mu }}$. Thus,

$$\begin{array}{cc}\hfill \Delta F_2& \le \left[N(n+1)-N_1\right]\left[\mathsf{\Pi }-\mu N(n+1)\right].\hfill \\ \hfill & =\left[N(n+1)-N_1\right]\left[\mu N_1-\mu N(n+1)\right].\hfill \\ \hfill & =-\mu {\left[N(n+1)-N_1\right]}^2\hfill \\ \hfill & \le 0\hfill \end{array}$$

Therefore, $F_2$ is a Lyapunov function on $\mathcal{D}/{\mathcal{D}}_0$ and hence, by applying the theorem of stability of difference equations (Theorem 6.3 in [36]), we obtain that every solution of the equations in model (1) with ${\delta }_1={\delta }_2=0$ approaching the unique endemic equilibrium point as $n\to \infty$ whenever ${\tilde{\mathcal{R}}}_0>1.$ □

6. Discussion and Conclusions

In this section, we will investigate the numerical simulation of the proposed model (1). The values of the model (1) parameters are listed in

Table 1. The values of the parameters of model (1) when ${\mathcal{R}}_0\le 1$ and when ${\mathcal{R}}_0>1$.

Parameter	${\mathcal{R}}_0\le 1$	${\mathcal{R}}_0>1$	Parameter	${\mathcal{R}}_0\le 1$	${\mathcal{R}}_0>1$
$\mathsf{\Pi }$	136	136	${\gamma }_1$	0.0337	0.0337
$\mu$	0.000035	0.000035	${\gamma }_2$	0.0486	0.0486
r	0.5	0.5	${\gamma }_3$	0.0535	0.0535
$\beta$	0.05	0.10	${\delta }_1$	0.1	0.1
$\sigma$	0.157	0.157	${\delta }_2$	0.1	0.1
$\varphi$	0.2	0.2	${\eta }_1$	0.03	0.03
${\eta }_2$	0.04	0.04

below.

Figure 3 (Up) shows that the disease dies out which means that the disease-free equilibrium point of model (1) is globally attractive and Figure 3 (Down) shows that the disease is permanent.

Figure 4 displays that the number of cumulative cases of infection with treatment is larger than cumulative cases of infection without treatment.

Figure 5 depicts the relation between the basic reproduction number ${\mathcal{R}}_0$ and the natural death rate $\mu$ for several values of the contact rate $\beta$. It shows that ${\mathcal{R}}_0$ decreases as $\mu$ increases and ${\mathcal{R}}_0$ increases as $\beta$ increases.