# The Continuous and Discrete Stability Characterization of Hepatitis B Deterministic Model

Authors: Shuo Li; Amjid Hussain; Ihsan Ullah Khan; Amine El Koufi; Arif Mehmood
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The hepatitis B infection is a global epidemic disease which is a huge risk to the public health. In this paper, the transmission dynamics of hepatitis B deterministic model are presented and studied. The basic reproduction number is attained and by applying it, the local as well as global stability of disease-free and endemic equilibria of continuous hepatitis B deterministic model are discussed. To better understand the dynamics of the disease, the discrete nonstandard finite difference (NSFD) scheme is produced for the continuous model. Different criteria are employed to check the local and global stability of disease-free and endemic equilibria for the NSFD scheme. Our findings demonstrate that the NSFD scheme is convergent for all step sizes and consequently reasonable in all respect for the continuous deterministic epidemic model. All the aforementioned properties and their effects are also proved numerically at each stage to show their mathematical as well as biological feasibility. The theoretical and numerical findings used in this paper can be employed as a helpful tool for predicting the transmission of other infectious diseases.


---

## Body

## 1. Introduction

Hepatitis is a general term that means inflammation of the liver. This disease can cause both acute and chronic infections. The acute stage is usually defined as the first six months of virus infection. During this phase, the immune system is capable to manipulate the infection of human body. The two primary indications of the acute stage are feeling sick and having a high temperature, which subsides after few weeks due to the immune system. Chronic disease affects the liver ability to perform life-sustaining processes such as removing dangerous transmitted substances from the blood, collecting sugar levels, and converting it to useable energy forms [1]. Hepatitis B virus (HBV) is one of the world most serious health problems [2]. HBV has a large rate of deaths, both from acute and chronic infection [3]. HBV is spread by blood transfusion and gets transmitted to the newly born child during pregnancy from affected mother. Vaccination is the most enchanting and efficient process in newly born children to decrease the occurrence of HBV [4].HBV can induce chronic infection which can lead to death from cirrhosis and liver cancer if not treated properly. The HBV fatality rates are among the higher causes of universal deaths [5]. Some pharmacological therapies for chronic HBV have been proposed, including alpha interferon, lamivudine, pegylated interferon, tenofovir disoproxil, entecavir, and telbivudine [6]. During treatment, the viral load is decreased, which reduces virus-related reproduction in the liver [7]. The vaccine against HBV is available since 1982, but still its transmission continues to rise [8, 9]. According to Sheikhan and Ghoreishi [6], HBV can also live beyond the mortal body. HBV can survive on the outer part of the body for at least seven days, and it can be transmitted to any unimmunized human body during this time.The use of mathematical modeling helps us to concentrate on the procedure by which an infectious disease spreads throughout an area. Many mathematical models are constructed by researchers from all around the world to understand different types of infectious diseases and their dynamic characteristics. In [10–13], the mathematical models of fractional order derivatives have been employed to investigate and evaluate the transmission of various infectious diseases. The authors not only examined the precise qualitative characteristics of the formulated models but also offered numerical simulations to verify the obtained theoretical findings. In [14, 15], the authors presented vaccination effects on HBV transmission with control strategies by using different age structures in the population. In [16–20], several specialized models of HBV transmission dynamics have been focused on the impact of commitment and control measures like vaccination and antiviral therapy. Din et al. [21] performed a detailed analysis of stability, showing that the reproduction number determines the entire dynamic activities of the system. Recently, in [22], the author discussed and analyzed the stochastic SACR model for HBV transmission and left the deterministic model unsolved. The author investigated the analytical results, including the stability of disease-free and endemic equilibria only for the continuous stochastic model. The purpose of the present work is continuous and discrete characterization of the hepatitis B deterministic model. Different criteria are used to discuss the local as well as global stability of disease-free and endemic equilibria for the continuous deterministic model. The discrete NSFD scheme is constructed for the continuous model to display its sustainability and biological suitability. The NSFD scheme constructed for the model is dynamically consistent with the original system for any step size. Our theoretical and numerical findings indicate that the NSFD scheme retains the essential qualitative characteristics of the continuous model. Consequently, this scheme is not only realistic but also verifies various features of the continuous model. The results acquired through this scheme are very precise and accurate.The paper is structured as follows: In Section2, the HBV epidemic model is presented and associated parameters are explained. The existing equilibria and reproduction number are established for the deterministic model in Section 3. By using the reproduction number, the local and global stability of disease-free and endemic equilibria for the continuous model are discussed in Section 4. The discrete NSFD scheme is constructed in Section 5 to analyze the convergence and divergence of disease-free and endemic equilibria for the proposed model. Our calculations show that the NSFD scheme is an effective and powerful technique that presents a clear portrait of the continuous model. The numerical simulations are also provided which strengthen our theoretical results. Finally, a brief conclusion is presented in the last section.


## 2. Mathematical Model for HBV

In order to define the stochastic HBV disease model with variable population environment, it is required to put some conditions on the epidemic model. It is assumed that the total populationNt at time t is divided into four classes, i.e., susceptible St, acutely infected At, chronically infected Ct, and recovered Rt where Nt=St+At+Ct+Rt. The second supposition is that all state variables and parameters of the proposed model are non-negative. The function Ω ⟶ Ω+ΩB with Ω2>0 denotes the concentration of white noise, where Bt shows the normal which satisfy B0=0. Based on all above information, the stochastic hepatitis B epidemic model [22] illustrated by the system of four stochastic differential equations is defined as follows:(1)dSt=p−βStCt−τ+μ1Stdt−ΩStAtdBt,dAt=βStCt−γ+μ1+ωAtdt+ΩStAtdBt,dCt=ωAt−γ2+μ1+εCtdt,dRt=τSt+γ2Ct+ωAt−μ1Rtdt.By puttingΩ  = 0, the model (1) deduces into the following deterministic model:(2)dStdt=p−βStCt−τ+μ1St,dAtdt=βStCt−γ+μ1+ωAt,dCtdt=ωAt−γ2+μ1+εCt,dRtdt=τSt+γ2Ct+ωAt−μ1Rt,dCtdt=ωAt−γ2+μ1+εCt.Our main aim is continuous and discrete characterization of model (2). It is assumed that all the parameters are positive constants where the parameters and their explanations are provided in Table 1.Table 1
Parameters included in model (2) and their explanations.
ParameterParameter descriptionpThe constant birth rateβThe transmission rateτThe vaccination rateμ1Natural death rateγThe constant recovery rate for acutely infected individualsωThe moving rate of acutely infected individuals to chronic stageγ2The disease induced death rateεThe constant recovery rate for chronically infected individualsAs the total populationN is denoted by(3)N=S+A+C+R.So, by employing model (2), we attain(4)dNdt=p−μ1S−γA−μ1A−μ1C+ωA−μ1R.From previous equation, we can write that(5)dNdt≤p−μ1S. (6)limt⟶∞supN≤pμ1.Therefore, the feasible region for model (2) becomes(7)K=S,A,C,R∈R4,N≤pμ1.


## 3. Equilibria of the Model and Basic Reproduction Number

### 3.1. Equilibria of the Model

The following two equilibria exist for model (2):


#### 3.1.1. Disease Free Equilibrium (DFE) Point

To find DFE point, we take all other classes equal to zero except the susceptible class, i.e., ifAt=0,Ct=0,Rt=0, then we get St=p/τ+μ1. Therefore, the DFE point denoted by E0S0,A0,C0,R0 becomes E0 = p/t+μ1,0,0,0.


#### 3.1.2. Disease Endemic Equilibrium (DEE) Point

To find DEE point, we simultaneously solve the proposed model (2) for S, A, C and R. If the DEE point is denoted by E∗S∗,A∗,C∗,R∗, then from model (2), we get(8)S∗t=γ+μ1+ωA∗tβγ2+μ1+εC∗t,A∗t=γ2+μ1+εC∗tω,C∗t=ωA∗tγ2+μ1+ε,R∗t=pτ+μ1μ1.


### 3.2. The Basic Reproduction NumberR0

The quantityR0 is the most crucial threshold related to any infectious disease. It assists to find out whether an infectious disease will transmit through population or not [23]. If R0<1 throughout its infectious period, then infection does not grow. On the other hand, if R0>1 then infection grows and disease remains in the population. To obtain R0, we employ transmission and translation matrices Fx and Vx, respectively. The previously discussed matrices can be demonstrated as(9)Fx=βStCt0andVx=−γ+μ1+ωAtωAt−γ2+μ1+εCt.From the abovementioned matrices, we get(10)F=0βpτ+μ100,V=γ+μ1+ω0−ωγ2+μ1+ε.As we know that(11)R0=ρFV−1.Therefore, usingF and V, we obtain(12)R0=βpωγ+μ1+ωγ2+μ1+ετ+μ1.


## 3.1. Equilibria of the Model

The following two equilibria exist for model (2):


### 3.1.1. Disease Free Equilibrium (DFE) Point

To find DFE point, we take all other classes equal to zero except the susceptible class, i.e., ifAt=0,Ct=0,Rt=0, then we get St=p/τ+μ1. Therefore, the DFE point denoted by E0S0,A0,C0,R0 becomes E0 = p/t+μ1,0,0,0.


### 3.1.2. Disease Endemic Equilibrium (DEE) Point

To find DEE point, we simultaneously solve the proposed model (2) for S, A, C and R. If the DEE point is denoted by E∗S∗,A∗,C∗,R∗, then from model (2), we get(8)S∗t=γ+μ1+ωA∗tβγ2+μ1+εC∗t,A∗t=γ2+μ1+εC∗tω,C∗t=ωA∗tγ2+μ1+ε,R∗t=pτ+μ1μ1.


## 3.1.1. Disease Free Equilibrium (DFE) Point

To find DFE point, we take all other classes equal to zero except the susceptible class, i.e., ifAt=0,Ct=0,Rt=0, then we get St=p/τ+μ1. Therefore, the DFE point denoted by E0S0,A0,C0,R0 becomes E0 = p/t+μ1,0,0,0.


## 3.1.2. Disease Endemic Equilibrium (DEE) Point

To find DEE point, we simultaneously solve the proposed model (2) for S, A, C and R. If the DEE point is denoted by E∗S∗,A∗,C∗,R∗, then from model (2), we get(8)S∗t=γ+μ1+ωA∗tβγ2+μ1+εC∗t,A∗t=γ2+μ1+εC∗tω,C∗t=ωA∗tγ2+μ1+ε,R∗t=pτ+μ1μ1.


## 3.2. The Basic Reproduction NumberR0

The quantityR0 is the most crucial threshold related to any infectious disease. It assists to find out whether an infectious disease will transmit through population or not [23]. If R0<1 throughout its infectious period, then infection does not grow. On the other hand, if R0>1 then infection grows and disease remains in the population. To obtain R0, we employ transmission and translation matrices Fx and Vx, respectively. The previously discussed matrices can be demonstrated as(9)Fx=βStCt0andVx=−γ+μ1+ωAtωAt−γ2+μ1+εCt.From the abovementioned matrices, we get(10)F=0βpτ+μ100,V=γ+μ1+ω0−ωγ2+μ1+ε.As we know that(11)R0=ρFV−1.Therefore, usingF and V, we obtain(12)R0=βpωγ+μ1+ωγ2+μ1+ετ+μ1.


## 4. Local and Global Stability of DFE and DEE Points for the Deterministic Model

In the following section, we first discuss the local and global stability of DFE point for the deterministic HBV disease model (2):


### 4.1. Local and Global Stability of DFE Point

To discuss the local stability, we assume(13)F=p−βStCt−τ+μ1St,G=βStCt−γ+μ1+ωAt,H=ωAt−γ2+μ1+εCt,I=τSt+γ2Ct+ωAt−μ1Rt.In the following theorem, we first discuss the local stability of DFE point by using Routh−Hurwitz criterion [24, 25]:Theorem 1.
The DFE pointE0 of model (2) is locally asymptotically stable whenever R0<1.Proof.
Let us take the Jacobian matrix as follows:(14)J=∂F∂S∂F∂A∂F∂C∂F∂R∂G∂S∂G∂A∂G∂C∂G∂R∂H∂S∂H∂A∂H∂C∂H∂R∂I∂S∂I∂A∂I∂C∂I∂R.
We first find all the derivatives included in (14) as follows:(15)∂F∂S=−βCt−τ+μ1,∂F∂A=0,∂F∂C=−βSt,∂F∂R=0,∂G∂S=βCt,∂G∂A=−γ+μ1+ω,∂G∂C=βSt,∂G∂R=0,∂H∂S=0,∂H∂A=ω,∂H∂C=−γ2+μ1+ε,∂H∂R=0,∂I∂S=τ,∂I∂A=ω,∂I∂C=γ2,∂I∂R=−μ1.
By replacing all the derivatives in (14), we get(16)J=−βCt−τ+μ10−βSt0βCt−γ+μ1+ωβSt00ω−γ2+μ1+ε0τωγ2−μ1.
By putting DFE pointE0 = p/t+μ1,0,0,0, we get(17)JE0=−τ+μ10−βpτ+μ100−γ+μ1+ωβpτ+μ100ω−γ2+μ1+ε0τωγ2−μ1.
In order to find eigenvalues, we consider(18)−τ+μ1−λ0−βpτ+μ100−γ+μ1+ω−λβpτ+μ100ω−γ2+μ1+ε−λ0τωγ2−μ1−λ=0.
The characteristic equation for the abovementioned equation becomes(19)−τ+μ1−λ−μ1−λλ2+λP+Q=0,where(20)P=γ+μ1+ω+γ2+μ1+ε,Q=γ+μ1+ωγ2+μ1+ε1−R0.
The two negative roots of (19) are λ1=−μ1 and λ2=−τ+μ1. Also, it is clear that P>0 and Q>0 whenever R0<1. So, by using Routh−Hurwitz criterion, the other two roots of λ2+λP+Q=0 must have negative real parts. Therefore, we deduce that E0 is locally asymptotically stable for R0<1. This completes the proof.Theorem 2.
The DFE pointE0 of model (2) is globally asymptotically stable whenever R0≤1.Proof.
In order to demonstrate the global stability of DFE pointE0 of model (2), we construct the Lyapunov function as(21)L=S−S0+A+C+R.
From (21), the following can easily be obtained:(22)dLdt=dSdt+dAdt+dCdt+dRdt.
After simple calculations, we get(23)dLdt=p−τ+μ1S−γ+μ1A−γ2+μ1+εC−μ1Rt,dLdt=−τ+μ1S−S0−γ+μ1A−γ2+μ1+εC−μ1Rt,dLdt=−τ+μ1S−S0+γ+μ1A+γ2+μ1+εC+μ1Rt≤0.
Thus,dL/dt≤0 for R0≤1. Also note that dL/dt=0 if and only if S=S0andA=C=R=0. Hence, Labzai et al. [26] imply that E0 is globally asymptotically stable, as shown in Figure 1(a).Figure 1
The solutions of HBV model (2) obtained through ODE-45. (a) Stable DFE point whenever R0≤1 and p=0.05, (b) stable DEE point whenever R0≥1 and p=2. Other parameters remain fixed, as β=0.60,γ=0.2,ω=0.5,τ=0.001,μ1=0.2,γ2=0.4,ε=0.5.
(a)(b)


### 4.2. Local and Globally Stability of DEE Point

Theorem 3.
The DEE pointE∗ of model (2) is locally asymptotically stable whenever R0>1.Proof.
In the similar way as in Theorem1, the Jacobian matrix can be written as(24)J=−βCt−τ+μ10−βSt0βCt−γ+μ1+ωβSt00ω−γ2+μ1+ε0τωγ2−μ1.
By putting DEE pointE∗, we get(25)JE∗=−βωA∗tγ2+μ1+ε−τ+μ10−βγ+μ1+ωA∗tγ2+μ1+εC∗t0βωA∗tγ2+μ1+ε−γ+μ1+ωβγ+μ1+ωA∗tγ2+μ1+εC∗t00ω−γ2+μ1+ε0τωγ2−μ1.
In order to find the eigenvalues, we consider(26)−βωA∗tγ2+μ1+ε−τ+μ1−λ0−βγ+μ1+ωA∗tγ2+μ1+εC∗t0βωA∗tγ2+μ1+ε−γ+μ1+ω−λβγ+μ1+ωA∗tγ2+μ1+εC∗t00ω−γ2+μ1+ε−λ0τωγ2−μ1−λ=0.
The characteristic equation of abovementioned equation becomes−μ1−λλ3+Bλ2+Dλ+E=0.
The abovementioned equation gives one negative eigenvalueλ1=−μ1. The other eigenvalues can be obtained from(27)λ3+Bλ2+Dλ+E=0,where(28)B=βωA∗tγ2+μ1+ε+τ+μ1+γ+μ1+ω+γ2+μ1+ε,D=βωA∗tγ+μ1+ωγ2+μ1+ε+βωA∗t+τ+μ1γ+μ1+ω+τ+μ1γ2+μ1+ε+γ+μ1+ωγ2+μ1+ε+βωγ+μ1+ωA∗tτ+μ1γ2+μ1+εC∗t. (29)E=βωA∗tτ+μ1γ+μ1+ωA∗tγ2+μ1+εC∗t1−γ+μ1+ωA∗tpR0>0.
It is clear thatB,D,E>0 whenever R0>1. Also,(30)BD−E=βωA∗tγ2+μ1+ε+τ+μ1+γ+μ1+ω+γ2+μ1+εβωA∗tγ+μ1+ωγ2+μ1+ε+βωA∗t+τ+μ1γ+μ1+ω+τ+μ1γ2+μ1+ε+γ+μ1+ωγ2+μ1+ε+βωγ+μ1+ωA∗tτ+μ1γ2+μ1+εC∗t−βωA∗tτ+μ1γ+μ1+ωA∗tγ2+μ1+εC∗t1−γ+μ1+ωA∗tpR0>0.
Hence, by applying Routh−Hurwitz criterion, all the solutions of (27) must have negative real parts if and only if R0>1. Therefore, E∗ is locally asymptotically stable whenever R0>1.Theorem 4.
The DEE pointE∗ of model (2) is globally asymptotically stable whenever R0≥1.Proof.
In order to demonstrate the global stability of DEE pointE∗ of model (2), we construct the Lyapunov function as(31)Ψ=12S−S∗+A−A∗+C−C∗+R−R∗2.
Now, we calculate the derivative with respect to the time of (31) and then using model (2), we get(32)dΨdt=S−S∗+A−A∗+C−C∗∗p−τ+μ1S−γ+μ1A−γ2+μ1+εC−μ1R.
If we putΦ1=τ+μ1,Φ2=γ+μ1,Φ3=γ2+μ1+ε,Φ4=μ1. Then, after simple arrangement, we get(33)dΨdt=S−S∗+A−A∗+C−C∗+R−R∗∗−S∗R0Φ1−SΦ1−Φ2A−Φ3C−Φ4R,dΨdt=−S−S∗+A−A∗+C−C∗+R−R∗∗S∗R0Φ1+SΦ1+Φ2A+Φ3C+Φ4R,dΨdt=−S−S∗+A−A∗+C−C∗+R−R∗∗S∗R0+SΦ1+Φ2A+Φ3C+Φ4R.
Since the right-hand side of (33) has a negative sign, so the derivative on right hand side is less than or equal to zero, i.e., dΨ/dt≤0. Substituting S=S∗,A=A∗,C=C∗,R=R∗ in (33), dΨ/dt yields zero, i.e., dΨ/dt=0. Therefore, the largest invariant set in S,A,C,RϵR4:dΨ/dt=0 is the singleton invariant set E∗, where E∗ is the DEE point. Then, by applying invariant principle of LaSalle et al. [27], it implies that E∗ is globally asymptotically stable, as shown in Figure 1(b).


## 4.1. Local and Global Stability of DFE Point

To discuss the local stability, we assume(13)F=p−βStCt−τ+μ1St,G=βStCt−γ+μ1+ωAt,H=ωAt−γ2+μ1+εCt,I=τSt+γ2Ct+ωAt−μ1Rt.In the following theorem, we first discuss the local stability of DFE point by using Routh−Hurwitz criterion [24, 25]:Theorem 1.
The DFE pointE0 of model (2) is locally asymptotically stable whenever R0<1.Proof.
Let us take the Jacobian matrix as follows:(14)J=∂F∂S∂F∂A∂F∂C∂F∂R∂G∂S∂G∂A∂G∂C∂G∂R∂H∂S∂H∂A∂H∂C∂H∂R∂I∂S∂I∂A∂I∂C∂I∂R.
We first find all the derivatives included in (14) as follows:(15)∂F∂S=−βCt−τ+μ1,∂F∂A=0,∂F∂C=−βSt,∂F∂R=0,∂G∂S=βCt,∂G∂A=−γ+μ1+ω,∂G∂C=βSt,∂G∂R=0,∂H∂S=0,∂H∂A=ω,∂H∂C=−γ2+μ1+ε,∂H∂R=0,∂I∂S=τ,∂I∂A=ω,∂I∂C=γ2,∂I∂R=−μ1.
By replacing all the derivatives in (14), we get(16)J=−βCt−τ+μ10−βSt0βCt−γ+μ1+ωβSt00ω−γ2+μ1+ε0τωγ2−μ1.
By putting DFE pointE0 = p/t+μ1,0,0,0, we get(17)JE0=−τ+μ10−βpτ+μ100−γ+μ1+ωβpτ+μ100ω−γ2+μ1+ε0τωγ2−μ1.
In order to find eigenvalues, we consider(18)−τ+μ1−λ0−βpτ+μ100−γ+μ1+ω−λβpτ+μ100ω−γ2+μ1+ε−λ0τωγ2−μ1−λ=0.
The characteristic equation for the abovementioned equation becomes(19)−τ+μ1−λ−μ1−λλ2+λP+Q=0,where(20)P=γ+μ1+ω+γ2+μ1+ε,Q=γ+μ1+ωγ2+μ1+ε1−R0.
The two negative roots of (19) are λ1=−μ1 and λ2=−τ+μ1. Also, it is clear that P>0 and Q>0 whenever R0<1. So, by using Routh−Hurwitz criterion, the other two roots of λ2+λP+Q=0 must have negative real parts. Therefore, we deduce that E0 is locally asymptotically stable for R0<1. This completes the proof.Theorem 2.
The DFE pointE0 of model (2) is globally asymptotically stable whenever R0≤1.Proof.
In order to demonstrate the global stability of DFE pointE0 of model (2), we construct the Lyapunov function as(21)L=S−S0+A+C+R.
From (21), the following can easily be obtained:(22)dLdt=dSdt+dAdt+dCdt+dRdt.
After simple calculations, we get(23)dLdt=p−τ+μ1S−γ+μ1A−γ2+μ1+εC−μ1Rt,dLdt=−τ+μ1S−S0−γ+μ1A−γ2+μ1+εC−μ1Rt,dLdt=−τ+μ1S−S0+γ+μ1A+γ2+μ1+εC+μ1Rt≤0.
Thus,dL/dt≤0 for R0≤1. Also note that dL/dt=0 if and only if S=S0andA=C=R=0. Hence, Labzai et al. [26] imply that E0 is globally asymptotically stable, as shown in Figure 1(a).Figure 1
The solutions of HBV model (2) obtained through ODE-45. (a) Stable DFE point whenever R0≤1 and p=0.05, (b) stable DEE point whenever R0≥1 and p=2. Other parameters remain fixed, as β=0.60,γ=0.2,ω=0.5,τ=0.001,μ1=0.2,γ2=0.4,ε=0.5.
(a)(b)


## 4.2. Local and Globally Stability of DEE Point

Theorem 3.
The DEE pointE∗ of model (2) is locally asymptotically stable whenever R0>1.Proof.
In the similar way as in Theorem1, the Jacobian matrix can be written as(24)J=−βCt−τ+μ10−βSt0βCt−γ+μ1+ωβSt00ω−γ2+μ1+ε0τωγ2−μ1.
By putting DEE pointE∗, we get(25)JE∗=−βωA∗tγ2+μ1+ε−τ+μ10−βγ+μ1+ωA∗tγ2+μ1+εC∗t0βωA∗tγ2+μ1+ε−γ+μ1+ωβγ+μ1+ωA∗tγ2+μ1+εC∗t00ω−γ2+μ1+ε0τωγ2−μ1.
In order to find the eigenvalues, we consider(26)−βωA∗tγ2+μ1+ε−τ+μ1−λ0−βγ+μ1+ωA∗tγ2+μ1+εC∗t0βωA∗tγ2+μ1+ε−γ+μ1+ω−λβγ+μ1+ωA∗tγ2+μ1+εC∗t00ω−γ2+μ1+ε−λ0τωγ2−μ1−λ=0.
The characteristic equation of abovementioned equation becomes−μ1−λλ3+Bλ2+Dλ+E=0.
The abovementioned equation gives one negative eigenvalueλ1=−μ1. The other eigenvalues can be obtained from(27)λ3+Bλ2+Dλ+E=0,where(28)B=βωA∗tγ2+μ1+ε+τ+μ1+γ+μ1+ω+γ2+μ1+ε,D=βωA∗tγ+μ1+ωγ2+μ1+ε+βωA∗t+τ+μ1γ+μ1+ω+τ+μ1γ2+μ1+ε+γ+μ1+ωγ2+μ1+ε+βωγ+μ1+ωA∗tτ+μ1γ2+μ1+εC∗t. (29)E=βωA∗tτ+μ1γ+μ1+ωA∗tγ2+μ1+εC∗t1−γ+μ1+ωA∗tpR0>0.
It is clear thatB,D,E>0 whenever R0>1. Also,(30)BD−E=βωA∗tγ2+μ1+ε+τ+μ1+γ+μ1+ω+γ2+μ1+εβωA∗tγ+μ1+ωγ2+μ1+ε+βωA∗t+τ+μ1γ+μ1+ω+τ+μ1γ2+μ1+ε+γ+μ1+ωγ2+μ1+ε+βωγ+μ1+ωA∗tτ+μ1γ2+μ1+εC∗t−βωA∗tτ+μ1γ+μ1+ωA∗tγ2+μ1+εC∗t1−γ+μ1+ωA∗tpR0>0.
Hence, by applying Routh−Hurwitz criterion, all the solutions of (27) must have negative real parts if and only if R0>1. Therefore, E∗ is locally asymptotically stable whenever R0>1.Theorem 4.
The DEE pointE∗ of model (2) is globally asymptotically stable whenever R0≥1.Proof.
In order to demonstrate the global stability of DEE pointE∗ of model (2), we construct the Lyapunov function as(31)Ψ=12S−S∗+A−A∗+C−C∗+R−R∗2.
Now, we calculate the derivative with respect to the time of (31) and then using model (2), we get(32)dΨdt=S−S∗+A−A∗+C−C∗∗p−τ+μ1S−γ+μ1A−γ2+μ1+εC−μ1R.
If we putΦ1=τ+μ1,Φ2=γ+μ1,Φ3=γ2+μ1+ε,Φ4=μ1. Then, after simple arrangement, we get(33)dΨdt=S−S∗+A−A∗+C−C∗+R−R∗∗−S∗R0Φ1−SΦ1−Φ2A−Φ3C−Φ4R,dΨdt=−S−S∗+A−A∗+C−C∗+R−R∗∗S∗R0Φ1+SΦ1+Φ2A+Φ3C+Φ4R,dΨdt=−S−S∗+A−A∗+C−C∗+R−R∗∗S∗R0+SΦ1+Φ2A+Φ3C+Φ4R.
Since the right-hand side of (33) has a negative sign, so the derivative on right hand side is less than or equal to zero, i.e., dΨ/dt≤0. Substituting S=S∗,A=A∗,C=C∗,R=R∗ in (33), dΨ/dt yields zero, i.e., dΨ/dt=0. Therefore, the largest invariant set in S,A,C,RϵR4:dΨ/dt=0 is the singleton invariant set E∗, where E∗ is the DEE point. Then, by applying invariant principle of LaSalle et al. [27], it implies that E∗ is globally asymptotically stable, as shown in Figure 1(b).


## 5. The NSFD Scheme

The main objective of this subsection is to develop a dynamically reliable discrete NSFD scheme for system (2). The NSFD scheme has been taken successfully to a variety of challenges, including ecology [28, 29], epidemiology [30, 31], and population models [32, 33]. To develop the NSFD scheme for system (2), we use Sn, An, Cn, and Rn as numerical approximations of St, AT,Ct,and Rt at t=nh, where n=0,1,2…., and h denotes the time-step size. By applying the concept of Mickens [34], we can discretize model (2) as follows:(34)Sn+1−Snh=p−βSn+1tCnt−τ+μ1Sn+1t,An+1−Anh=βSn+1tCnt−γ+μ1+ωAn+1t,Cn+1−Cnh=ωAnt−γ2+μ1+εCn+1t,Rn+1−Rnh=τSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t.The discrete NSFD model (10) can be rearranged to get explicit form as(35)Sn+1=Sn+hp1+hβCnt+hτ+μ1,An+1=An+hβCntSn+1t1+hγ+μ1+ω,Cn+1=Cn+hωAnt1+hγ2+μ1+ε,Rn+1=Rn+hτSn+1t+γ2Cn+1t+ωAn+1t1+hμ1.


### 5.1. Local and Global Stability of DFE Point for NSFD Scheme

To obtain the local stability of the DFE point, we assume that(36)F1=Sn+1=Sn+hp1+hβCnt+hτ+μ1,F2=An+1=An+hβCntSn+1t1+hγ+μ1+ω,F3=Cn+1=Cn+hωAnt1+hγ2+μ1+ε,F4=Rn+1=Rn+hτSn+1t+γ2Cn+1t+ωAn+1t1+hμ1.Theorem 5.
For allh>0, the DFE point E0 of the NSFD model (11) is locally asymptotically stable whenever R0<1.Proof.
Let us take the Jacobian matrix as follows:(37)J=∂F1∂S∂F1∂A∂F1∂C∂F1∂R∂F2∂S∂F2∂A∂F2∂C∂F2∂R∂F3∂S∂F3∂A∂F3∂C∂F3∂R∂F4∂S∂F4∂A∂F4∂C∂F4∂R.
First, we find all the derivatives of matrix (13) as follows:(38)∂F1∂S=11+hβCnt+τ+μ1,∂F1∂A=0,∂F1∂C=−Sn+hphβ1+hβCnt+τ+μ12,∂F1∂R=0,∂F2∂S=hβCnt1+hγ+μ1+ω,∂F2∂A=11+hγ+μ1+ω,∂F2∂C=hβSn+1t1+hγ+μ1+ω,∂F2∂R=0,∂F3∂S=0,∂F3∂A=hω1+hγ2+μ1+ε,∂F3∂C=11+hγ2+μ1+ε,∂F3∂R=0,∂F4∂S=hτ1+hμ1,∂F4∂A=hω1+hμ1,∂F4∂C=hγ21+hμ1,∂F4∂R=11+hμ1.
By replacing all the derivatives in (37), we get(39)J=11+hβCnt+τ+μ10−Sn+hphβ1+hβCnt+τ+μ120hβCnt1+hγ+μ1+ω11+hγ+μ1+ωhβSn+1t1+hγ+μ1+ω00hω1+hγ2+μ1+ε11+hγ2+μ1+ε0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1.
After putting the DFE pointE0, we get(40)JE0=11+hτ+μ10−p/t+μ1+hphβ1+hτ+μ120011+hγ+μ1+ωhβp1+hγ+μ1+ωτ+μ100hω1+hγ2+μ1+ε11+hγ2+μ1+ε0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1.
In order to find the eigenvalues, we consider(41)11+hτ+μ1−λ0−p/t+μ1+hphβ1+hτ+μ120011+hγ+μ1+ω−λhβp1+hγ+μ1+ωτ+μ100hω1+hγ2+μ1+ε11+hγ2+μ1+ε−λ0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1−λ=0.
The abovementioned equation gives the following characteristic equation:(42)11+hμ1−λ11+hτ+μ1−λλ2+Lλ+M=0,where(43)L=11+hγ+μ1+ω+11+hγ2+μ1+ε,M=hβpω1+hγ2+μ1+ε1+hγ+μ1+ωτ+μ1.
The two roots of the (42) are λ1=1/1+hμ1<1 and λ2=1/1+hτ+μ1<1. Also, it is clear that L>0 and M>0 whenever R0<1. So, by using Routh−Hurwitz criterion, the other two roots of λ2+λL+M=0 must have negative real parts. Therefore, we conclude that the DFE point E0 of the discrete NSFD model (11) is locally asymptotically stable whenever R0<1. This completes the proof.
In the following theorem, we now prove the global stability of the DFE pointE0. To prove it, we use the criterion employed by Vaz and Torres [35].Theorem 6.
For allh>0, the DFE point E0 of the NSFD model (11) is globally asymptotically stable whenever R0≤1, as shown in Figures 2(a)–2(d).Figure 2
The solutions of the HBV model (2) obtained through the NSFD scheme whenever R0≤1 and (a) h=0.01, (b) h=1, (c) h=10 and (d) h=20. Other parameters remain fixed as p=0.5,β=0.60,γ=0.2,ω=0.5,τ=0.001,μ1=0.2,γ2=0.4,ε=0.5.
(a)(b)(c)(d)Proof.
If we chooseε>0, then there exists an integer n0, such that for any n≥n0,Sn+1<p/τ+μ1+ε. We consider the sequence Tn defined by(44)Tn=hβSn+1Cn+An+ωC3Cn+μ1C2Rn,where C2=γ2+μ1+ε,C3=γ+μ1+ω. From abovementioned equation, we can write that(45)Tn+1−Tn=hβSn+2Cn+An+1+ωC3Cn+1+μ1C2Rn+1−hβSn+1Cn−An−ωC3Cn−μ1C2Rn,=hβSn+2Cn+An+1−An+ωC3Cn+1−Cn+μ1C2Rn+1−Rn.
After simple calculations, we obtain(46)=hβSn+2Cn+hβSn+1tCnt−γ+μ1+ωAn+1t+ωC3hωAnt−γ2+μ1+εCn+1t+μ1C2hτSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t−hβSn+1Cn.(47)=hβSn+2Cn−hγ+μ1+ωAn+1t+ωC3hωAnt−γ2+μ1+εCn+1t+μ1C2hτSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t
LetC1=γ+μ1+ω, then(48)=hβSn+2Cn−hC1An+1t+ωC3hωAnt−γ2+μ1+εCn+1t+μ1C2hτSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t=hβSn+2Cn−C1An+1t+ωC3ωAnt−γ2+μ1+εCn+1t+μ1C2τSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t=hβSn+2Cn+ω−C1An+1t+μ1C2γ2−ωC3γ2+μ1+εCn+1t+μ1C2τSn+1t−μ1Rn+1t.
IfU=μ1γ2C3−C2ωγ2+μ1+ε, then(49)=hβSn+2Cn−μ1γ2C3−C2ωγ2+μ1+εCn+1t+ω−C1An+1t+μ1C2τSn+1t−UCn+1t−μ1Rn+1t.
We can selectβ, a very small positive number such that(50)βSn+2Cn≤qAn+1+C1Cn+1+C2Rn+1.
After simple rearrangement, we can write that(51)βSn+2Cn≤qAn+1+C1An+1+C2An+1,≤qAn+1+C1ωC2An+1−C2Uωγ2C2+ωγ2C2An+1−UAn+1C3,≤qAn+11+C1ωC2+C2ωC2γ2−UC3,≤qAn+11+βpτ+μ1C1R0βp−ωγ2U+1−UC3.
IfR0≤1, and because it is imprecise, we reach the conclusion that Tn+1−Tn≤0 and limn⟶∞In=0 for any n≥0. The sequence Tnn=0∞ is a monotonic decreasing and limn⟶∞Sn=p/τ+μ1. Hence, the DFE point E0 is globally asymptotically stable.


### 5.2. Local and Global Stability of DEE Point for NSFD Scheme

Theorem 7.
For allh>0, the DEE point E∗ of the NSFD model (11) is locally asymptotically stable whenever R0>1.Proof.
In the similar way as in Theorem5, the Jacobian matrix can be obtained as(52)J=11+hβCnt+τ+μ10−Sn+hphβ1+hβCnt+τ+μ120hβCt1+hγ+μ1+ω11+hγ+μ1+ωhβSt1+hγ+μ1+ω00hω1+hγ2+μ1+ε11+hγ2+μ1+ε0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1.
By putting DEE pointE∗, we get(53)JE∗=11+hβC∗+τ+μ10−S∗+hphβ1+hβCnt+τ+μ120hβC∗1+hγ+μ1+ω11+hγ+μ1+ωhβS∗1+hγ+μ1+ω00hω1+hγ2+μ1+ε11+hγ2+μ1+ε0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1.
To find the eigenvalues of (53), we consider(54)11+hβC∗+τ+μ1−λ0−S∗+hphβ1+hβCnt+τ+μ120hβC∗1+hγ+μ1+ω11+hγ+μ1+ω−λhβS∗1+hγ+μ1+ω00hω1+hγ2+μ1+ε11+hγ2+μ1+ε−λ0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1−λ=0.
The abovementioned equation gives the following deterministic equation:(55)11+hμ1−λλ3+λ2K+λL+M=0.
The abovementioned equation provides one eigenvalueλ1=1/1+hμ1<1. The remaining eigenvalues can be obtained from(56)λ3+λ2K+λL+M=0,where(57)K=11+hγ+μ1+ω+11+hγ2+μ1+ε−11+hβC∗+τ+μ1,L=11+hγ+μ1+ω1+hβC∗+τ+μ1+11+hγ2+μ1+ε1+hβC∗+τ+μ1+11+hγ+μ1+ω1+hγ2+μ1+ε+h2βS∗ω1+hγ2+μ1+ε,M=h3β2C∗S∗+hpω1+hγ+μ1+ω1+hγ2+μ1+ε1+hβC∗t+τ+μ12−11+hγ+μ1+ω1+hγ2+μ1+ε1+hβC∗+τ+μ1.
It is clear thatK,L,M>0 whenever R0>1. Also,(58)KL−M=11+hγ+μ1+ω+11+hγ2+μ1+ε−11+hβC∗+τ+μ111+hγ+μ1+ω1+hβC∗+τ+μ1+11+hγ2+μ1+ε1+hβC∗+τ+μ1+11+hγ+μ1+ω1+hγ2+μ1+ε+h2βS∗ω1+hγ2+μ1+ε−h3β2C∗S∗+hpω1+hγ+μ1+ω1+hγ2+μ1+ε1+hβC∗t+τ+μ12+11+hγ+μ1+ω1+hγ2+μ1+ε1+hβC∗+τ+μ1>0.
Therefore, by applying Routh−Hurwitz criterion, all the solutions of the (56) must have negative real parts whenever R0>1. Hence, the DEE point E∗ of the NSFD scheme is locally asymptotically stable whenever R0>1.Theorem 8.
For allh>0, the DEE point E∗ of the NSFD model (11) is globally asymptotically stable whenever R0≥1, as shown in Figures 3(a)–3(d)).Figure 3
The solutions of HBV model (2) obtained through NSFD scheme whenever R0≥1 and (a) h=0.01, (b) h=1, (c) h=10 and (d) h=20. Other parameters remain fixed as p=0.5,β=1.60,γ=0.2,ω=0.5,τ=0.001,μ1=0.1,γ2=0.4,ε=0.5.
(a)(b)(c)(d)Proof.
We construct a sequenceVn such that(59)Vn=1hA∗fSnS∗+1hS∗fAnA∗+ωC∗hC3S∗C∗fCnC∗+βR∗hC2S∗A∗fRnR∗.where C2=γ+μ1+ω and C3 = γ2+μ1+ε. Let fx=x−1−lnx,x∈R+. Clearly fx≥0, and the equality holds true if x=1. We have(60)fSn+1S∗−fSnS∗=Sn+1−SnS∗−lnSn+1Sn≤Sn+1−S∗Sn+1−SnSn+1S∗=Sn+1−S∗Sn+1S∗hp−βSn+1tCnt−τ+μ1Sn+1t=Sn+1−S∗Sn+1S∗hβS∗ttC∗+τ+μ1S∗t−βSn+1tCnt−τ+μ1Sn+1t=hSn+1−S∗−τ+μ1Sn+1−S∗−βC∗Sn+1−S∗=−τ+μ1Sn+1−S∗2Sn+1S∗+βhC∗1−S∗Sn+1AnA∗Sn+1S∗−1.
Similarly.(61)fAn+1A∗−fAnA∗=An+1−AnA∗−lnAn+1An≤An+1−A∗An+1A∗≤An+1−A∗An+1A∗βSn+1tCnt−γ+μ1+ωAn+1t≤An+1−A∗An+1A∗βSn+1tCnt−Cn+1S∗C∗−ωCn+1A∗C∗+ωAn+1t=1−A∗An+1ωhA∗C∗An+1A∗−Cn+1C∗+1−A∗An+1hS∗C∗C∗βCnSn+1S∗−Cn+1S∗C∗fCn+1C∗−fCnC∗=Cn+1−CnC∗−lnCn+1Cn≤Cn+1−C∗Cn+1−CnCn+1C∗≤Cn+1−C∗Cn+1C∗ωAnt−γ2+μ1+εCnt≤Cn+1−C∗Cn+1C∗ωAnt−C3Cn+1=C3h1−C∗Cn+1An+1A∗−Cn+1C∗. (62)fRn+1R∗−fRnR∗=Rn+1−RnR∗−lnRn+1Rn≤Rn+1−R∗Rn+1−RnRn+1R∗≤Rn+1−R∗Rn+1R∗τSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t=h1−R∗Rn+1An+1A∗−Rn+1R∗.
The difference ofTn satisfies(63)Tn+1−Tn=1hA∗fSn+1S∗−fSnS∗+1hS∗fAn+1A∗−fAnA∗,+ωI∗hC3S∗A∗fCn+1C∗−fCnC∗+βR∗hC2S∗A∗fRn+1R∗−fRnR∗,≤−τ+μ1Sn+1−S∗2Sn+1S∗+βhC∗1−S∗Sn+1AnA∗Sn+1S∗−1−1−A∗An+1ωhA∗C∗An+1A∗−Cn+1C∗+1−A∗An+1hS∗C∗C∗βCnSn+1S∗−Cn+1S∗C∗−ωI∗hC3S∗A∗Cn+1−C∗Cn+1C∗ωAnt−C3Cn+1−βR∗hC2S∗A∗1−R∗Rn+1An+1A∗−Rn+1R∗,≤−βhSn+1−S∗2Sn+1S∗−ωhA∗C∗C∗Cn+1AnA∗Sn+1S∗−2−AnA∗+S∗Sn+1+An+1A∗,−1hS∗C∗An+1Cn+1C∗+A∗Cn+1An+1C∗−2−βR∗hC2S∗A∗C∗Rn+1Cn+1R∗+R∗Cn+1Rn+1C∗−2,≤−βhSn+1−S∗2Sn+1S∗−1hA∗fS∗Sn+1+fAn+1A∗+fC∗Cn+1AnA∗Sn+1S∗−fAnA∗−1hS∗fA∗Cn+1An+1C∗+fC∗An+1Cn+1A∗−βR∗hC2S∗A∗fA∗Rn+1An+1R∗+fR∗An+1Rn+1A∗.
Therefore,Tn is a monotonic decreasing sequence for any n≥0. Since Tn≥0 and limn⟶∞Tn+1−Tn=0, we obtain limn⟶∞Sn+1=S∗,limn⟶∞An+1=A∗, limn⟶∞Cn+1=C∗ and limn⟶∞Rn+1=R∗. Hence, according to Vaz et al. [35], we conclude that when R0≥1, the DEE point E∗ is globally asymptotically stable.


## 5.1. Local and Global Stability of DFE Point for NSFD Scheme

To obtain the local stability of the DFE point, we assume that(36)F1=Sn+1=Sn+hp1+hβCnt+hτ+μ1,F2=An+1=An+hβCntSn+1t1+hγ+μ1+ω,F3=Cn+1=Cn+hωAnt1+hγ2+μ1+ε,F4=Rn+1=Rn+hτSn+1t+γ2Cn+1t+ωAn+1t1+hμ1.Theorem 5.
For allh>0, the DFE point E0 of the NSFD model (11) is locally asymptotically stable whenever R0<1.Proof.
Let us take the Jacobian matrix as follows:(37)J=∂F1∂S∂F1∂A∂F1∂C∂F1∂R∂F2∂S∂F2∂A∂F2∂C∂F2∂R∂F3∂S∂F3∂A∂F3∂C∂F3∂R∂F4∂S∂F4∂A∂F4∂C∂F4∂R.
First, we find all the derivatives of matrix (13) as follows:(38)∂F1∂S=11+hβCnt+τ+μ1,∂F1∂A=0,∂F1∂C=−Sn+hphβ1+hβCnt+τ+μ12,∂F1∂R=0,∂F2∂S=hβCnt1+hγ+μ1+ω,∂F2∂A=11+hγ+μ1+ω,∂F2∂C=hβSn+1t1+hγ+μ1+ω,∂F2∂R=0,∂F3∂S=0,∂F3∂A=hω1+hγ2+μ1+ε,∂F3∂C=11+hγ2+μ1+ε,∂F3∂R=0,∂F4∂S=hτ1+hμ1,∂F4∂A=hω1+hμ1,∂F4∂C=hγ21+hμ1,∂F4∂R=11+hμ1.
By replacing all the derivatives in (37), we get(39)J=11+hβCnt+τ+μ10−Sn+hphβ1+hβCnt+τ+μ120hβCnt1+hγ+μ1+ω11+hγ+μ1+ωhβSn+1t1+hγ+μ1+ω00hω1+hγ2+μ1+ε11+hγ2+μ1+ε0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1.
After putting the DFE pointE0, we get(40)JE0=11+hτ+μ10−p/t+μ1+hphβ1+hτ+μ120011+hγ+μ1+ωhβp1+hγ+μ1+ωτ+μ100hω1+hγ2+μ1+ε11+hγ2+μ1+ε0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1.
In order to find the eigenvalues, we consider(41)11+hτ+μ1−λ0−p/t+μ1+hphβ1+hτ+μ120011+hγ+μ1+ω−λhβp1+hγ+μ1+ωτ+μ100hω1+hγ2+μ1+ε11+hγ2+μ1+ε−λ0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1−λ=0.
The abovementioned equation gives the following characteristic equation:(42)11+hμ1−λ11+hτ+μ1−λλ2+Lλ+M=0,where(43)L=11+hγ+μ1+ω+11+hγ2+μ1+ε,M=hβpω1+hγ2+μ1+ε1+hγ+μ1+ωτ+μ1.
The two roots of the (42) are λ1=1/1+hμ1<1 and λ2=1/1+hτ+μ1<1. Also, it is clear that L>0 and M>0 whenever R0<1. So, by using Routh−Hurwitz criterion, the other two roots of λ2+λL+M=0 must have negative real parts. Therefore, we conclude that the DFE point E0 of the discrete NSFD model (11) is locally asymptotically stable whenever R0<1. This completes the proof.
In the following theorem, we now prove the global stability of the DFE pointE0. To prove it, we use the criterion employed by Vaz and Torres [35].Theorem 6.
For allh>0, the DFE point E0 of the NSFD model (11) is globally asymptotically stable whenever R0≤1, as shown in Figures 2(a)–2(d).Figure 2
The solutions of the HBV model (2) obtained through the NSFD scheme whenever R0≤1 and (a) h=0.01, (b) h=1, (c) h=10 and (d) h=20. Other parameters remain fixed as p=0.5,β=0.60,γ=0.2,ω=0.5,τ=0.001,μ1=0.2,γ2=0.4,ε=0.5.
(a)(b)(c)(d)Proof.
If we chooseε>0, then there exists an integer n0, such that for any n≥n0,Sn+1<p/τ+μ1+ε. We consider the sequence Tn defined by(44)Tn=hβSn+1Cn+An+ωC3Cn+μ1C2Rn,where C2=γ2+μ1+ε,C3=γ+μ1+ω. From abovementioned equation, we can write that(45)Tn+1−Tn=hβSn+2Cn+An+1+ωC3Cn+1+μ1C2Rn+1−hβSn+1Cn−An−ωC3Cn−μ1C2Rn,=hβSn+2Cn+An+1−An+ωC3Cn+1−Cn+μ1C2Rn+1−Rn.
After simple calculations, we obtain(46)=hβSn+2Cn+hβSn+1tCnt−γ+μ1+ωAn+1t+ωC3hωAnt−γ2+μ1+εCn+1t+μ1C2hτSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t−hβSn+1Cn.(47)=hβSn+2Cn−hγ+μ1+ωAn+1t+ωC3hωAnt−γ2+μ1+εCn+1t+μ1C2hτSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t
LetC1=γ+μ1+ω, then(48)=hβSn+2Cn−hC1An+1t+ωC3hωAnt−γ2+μ1+εCn+1t+μ1C2hτSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t=hβSn+2Cn−C1An+1t+ωC3ωAnt−γ2+μ1+εCn+1t+μ1C2τSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t=hβSn+2Cn+ω−C1An+1t+μ1C2γ2−ωC3γ2+μ1+εCn+1t+μ1C2τSn+1t−μ1Rn+1t.
IfU=μ1γ2C3−C2ωγ2+μ1+ε, then(49)=hβSn+2Cn−μ1γ2C3−C2ωγ2+μ1+εCn+1t+ω−C1An+1t+μ1C2τSn+1t−UCn+1t−μ1Rn+1t.
We can selectβ, a very small positive number such that(50)βSn+2Cn≤qAn+1+C1Cn+1+C2Rn+1.
After simple rearrangement, we can write that(51)βSn+2Cn≤qAn+1+C1An+1+C2An+1,≤qAn+1+C1ωC2An+1−C2Uωγ2C2+ωγ2C2An+1−UAn+1C3,≤qAn+11+C1ωC2+C2ωC2γ2−UC3,≤qAn+11+βpτ+μ1C1R0βp−ωγ2U+1−UC3.
IfR0≤1, and because it is imprecise, we reach the conclusion that Tn+1−Tn≤0 and limn⟶∞In=0 for any n≥0. The sequence Tnn=0∞ is a monotonic decreasing and limn⟶∞Sn=p/τ+μ1. Hence, the DFE point E0 is globally asymptotically stable.


## 5.2. Local and Global Stability of DEE Point for NSFD Scheme

Theorem 7.
For allh>0, the DEE point E∗ of the NSFD model (11) is locally asymptotically stable whenever R0>1.Proof.
In the similar way as in Theorem5, the Jacobian matrix can be obtained as(52)J=11+hβCnt+τ+μ10−Sn+hphβ1+hβCnt+τ+μ120hβCt1+hγ+μ1+ω11+hγ+μ1+ωhβSt1+hγ+μ1+ω00hω1+hγ2+μ1+ε11+hγ2+μ1+ε0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1.
By putting DEE pointE∗, we get(53)JE∗=11+hβC∗+τ+μ10−S∗+hphβ1+hβCnt+τ+μ120hβC∗1+hγ+μ1+ω11+hγ+μ1+ωhβS∗1+hγ+μ1+ω00hω1+hγ2+μ1+ε11+hγ2+μ1+ε0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1.
To find the eigenvalues of (53), we consider(54)11+hβC∗+τ+μ1−λ0−S∗+hphβ1+hβCnt+τ+μ120hβC∗1+hγ+μ1+ω11+hγ+μ1+ω−λhβS∗1+hγ+μ1+ω00hω1+hγ2+μ1+ε11+hγ2+μ1+ε−λ0hτ1+hμ1hω1+hμ1hγ21+hμ111+hμ1−λ=0.
The abovementioned equation gives the following deterministic equation:(55)11+hμ1−λλ3+λ2K+λL+M=0.
The abovementioned equation provides one eigenvalueλ1=1/1+hμ1<1. The remaining eigenvalues can be obtained from(56)λ3+λ2K+λL+M=0,where(57)K=11+hγ+μ1+ω+11+hγ2+μ1+ε−11+hβC∗+τ+μ1,L=11+hγ+μ1+ω1+hβC∗+τ+μ1+11+hγ2+μ1+ε1+hβC∗+τ+μ1+11+hγ+μ1+ω1+hγ2+μ1+ε+h2βS∗ω1+hγ2+μ1+ε,M=h3β2C∗S∗+hpω1+hγ+μ1+ω1+hγ2+μ1+ε1+hβC∗t+τ+μ12−11+hγ+μ1+ω1+hγ2+μ1+ε1+hβC∗+τ+μ1.
It is clear thatK,L,M>0 whenever R0>1. Also,(58)KL−M=11+hγ+μ1+ω+11+hγ2+μ1+ε−11+hβC∗+τ+μ111+hγ+μ1+ω1+hβC∗+τ+μ1+11+hγ2+μ1+ε1+hβC∗+τ+μ1+11+hγ+μ1+ω1+hγ2+μ1+ε+h2βS∗ω1+hγ2+μ1+ε−h3β2C∗S∗+hpω1+hγ+μ1+ω1+hγ2+μ1+ε1+hβC∗t+τ+μ12+11+hγ+μ1+ω1+hγ2+μ1+ε1+hβC∗+τ+μ1>0.
Therefore, by applying Routh−Hurwitz criterion, all the solutions of the (56) must have negative real parts whenever R0>1. Hence, the DEE point E∗ of the NSFD scheme is locally asymptotically stable whenever R0>1.Theorem 8.
For allh>0, the DEE point E∗ of the NSFD model (11) is globally asymptotically stable whenever R0≥1, as shown in Figures 3(a)–3(d)).Figure 3
The solutions of HBV model (2) obtained through NSFD scheme whenever R0≥1 and (a) h=0.01, (b) h=1, (c) h=10 and (d) h=20. Other parameters remain fixed as p=0.5,β=1.60,γ=0.2,ω=0.5,τ=0.001,μ1=0.1,γ2=0.4,ε=0.5.
(a)(b)(c)(d)Proof.
We construct a sequenceVn such that(59)Vn=1hA∗fSnS∗+1hS∗fAnA∗+ωC∗hC3S∗C∗fCnC∗+βR∗hC2S∗A∗fRnR∗.where C2=γ+μ1+ω and C3 = γ2+μ1+ε. Let fx=x−1−lnx,x∈R+. Clearly fx≥0, and the equality holds true if x=1. We have(60)fSn+1S∗−fSnS∗=Sn+1−SnS∗−lnSn+1Sn≤Sn+1−S∗Sn+1−SnSn+1S∗=Sn+1−S∗Sn+1S∗hp−βSn+1tCnt−τ+μ1Sn+1t=Sn+1−S∗Sn+1S∗hβS∗ttC∗+τ+μ1S∗t−βSn+1tCnt−τ+μ1Sn+1t=hSn+1−S∗−τ+μ1Sn+1−S∗−βC∗Sn+1−S∗=−τ+μ1Sn+1−S∗2Sn+1S∗+βhC∗1−S∗Sn+1AnA∗Sn+1S∗−1.
Similarly.(61)fAn+1A∗−fAnA∗=An+1−AnA∗−lnAn+1An≤An+1−A∗An+1A∗≤An+1−A∗An+1A∗βSn+1tCnt−γ+μ1+ωAn+1t≤An+1−A∗An+1A∗βSn+1tCnt−Cn+1S∗C∗−ωCn+1A∗C∗+ωAn+1t=1−A∗An+1ωhA∗C∗An+1A∗−Cn+1C∗+1−A∗An+1hS∗C∗C∗βCnSn+1S∗−Cn+1S∗C∗fCn+1C∗−fCnC∗=Cn+1−CnC∗−lnCn+1Cn≤Cn+1−C∗Cn+1−CnCn+1C∗≤Cn+1−C∗Cn+1C∗ωAnt−γ2+μ1+εCnt≤Cn+1−C∗Cn+1C∗ωAnt−C3Cn+1=C3h1−C∗Cn+1An+1A∗−Cn+1C∗. (62)fRn+1R∗−fRnR∗=Rn+1−RnR∗−lnRn+1Rn≤Rn+1−R∗Rn+1−RnRn+1R∗≤Rn+1−R∗Rn+1R∗τSn+1t+γ2Cn+1t+ωAn+1t−μ1Rn+1t=h1−R∗Rn+1An+1A∗−Rn+1R∗.
The difference ofTn satisfies(63)Tn+1−Tn=1hA∗fSn+1S∗−fSnS∗+1hS∗fAn+1A∗−fAnA∗,+ωI∗hC3S∗A∗fCn+1C∗−fCnC∗+βR∗hC2S∗A∗fRn+1R∗−fRnR∗,≤−τ+μ1Sn+1−S∗2Sn+1S∗+βhC∗1−S∗Sn+1AnA∗Sn+1S∗−1−1−A∗An+1ωhA∗C∗An+1A∗−Cn+1C∗+1−A∗An+1hS∗C∗C∗βCnSn+1S∗−Cn+1S∗C∗−ωI∗hC3S∗A∗Cn+1−C∗Cn+1C∗ωAnt−C3Cn+1−βR∗hC2S∗A∗1−R∗Rn+1An+1A∗−Rn+1R∗,≤−βhSn+1−S∗2Sn+1S∗−ωhA∗C∗C∗Cn+1AnA∗Sn+1S∗−2−AnA∗+S∗Sn+1+An+1A∗,−1hS∗C∗An+1Cn+1C∗+A∗Cn+1An+1C∗−2−βR∗hC2S∗A∗C∗Rn+1Cn+1R∗+R∗Cn+1Rn+1C∗−2,≤−βhSn+1−S∗2Sn+1S∗−1hA∗fS∗Sn+1+fAn+1A∗+fC∗Cn+1AnA∗Sn+1S∗−fAnA∗−1hS∗fA∗Cn+1An+1C∗+fC∗An+1Cn+1A∗−βR∗hC2S∗A∗fA∗Rn+1An+1R∗+fR∗An+1Rn+1A∗.
Therefore,Tn is a monotonic decreasing sequence for any n≥0. Since Tn≥0 and limn⟶∞Tn+1−Tn=0, we obtain limn⟶∞Sn+1=S∗,limn⟶∞An+1=A∗, limn⟶∞Cn+1=C∗ and limn⟶∞Rn+1=R∗. Hence, according to Vaz et al. [35], we conclude that when R0≥1, the DEE point E∗ is globally asymptotically stable.


## 6. Conclusions

In the present paper, we proposed and investigated the deterministic HBV disease model. The aim is to understand the disease dynamics and generate methods for managing the disease transmission among the people. The continuous as well as discrete description of the model is discussed, and different essential properties are developed. Both continuous and discrete outcomes can be used as efficient tool for monitoring the transmission of HBV epidemic diseases. The basic reproduction numberR0 is constructed for the model and by applying it, the local as well as global stability of DFE and DEE points are discussed for the continuous model. The discrete NSFD scheme is developed for the model which is not only unconditionally convergent but also produces precise findings that are mathematically and biologically consistent with the associated continuous model. The local and global stability of both DFE and DEE points are demonstrated for the NSFD scheme by utilizing various criteria and conditions. The qualities of the NSFD scheme are presented, which shows that the results of the NSFD scheme are qualitatively exact and proficient. The NSFD scheme is a simple technique which shows that the continuous and discrete models behave identically, providing mathematically stable and high-quality response. By using the NSFD scheme, we are able to better explain the results that would be valuable to society and to the area of medicine as well. The outcomes provided in this paper can be used as a helpful apparatus for forecasting the development of HBV epidemic diseases. To support our theoretical results, numerical simulations are offered at every stage.



---

*Source: 1893665-2022-11-16.xml*