# Ergodic Stationary Distribution of a Stochastic Hepatitis B Epidemic Model with Interval-Valued Parameters and Compensated Poisson Process

Authors: Driss Kiouach; Yassine Sabbar
Journal: Computational and Mathematical Methods in Medicine
 (2020)

---

## Abstract

Hepatitis B epidemic was and is still a rich subject that sparks the interest of epidemiological researchers. The dynamics of this epidemic is often modeled by a system with constant parameters. In reality, the parameters associated with the Hepatitis B model are not certain, but the interval in which it belongs to can readily be determined. Our paper focuses on an imprecise Hepatitis B model perturbed by Lévy noise due to unexpected environmental disturbances. This model has a global positive solution. Under an appropriate assumption, we prove the existence of a unique ergodic stationary distribution by using the mutually exclusive possibilities lemma demonstrated by Stettner in 1986. Our main effort is to establish an almost perfect condition for the existence of the stationary distribution. Numerical simulations are introduced to illustrate the analytical results.


---

## Body

## 1. Introduction

Hepatitis B is an enormous defiance and a great global health issue caused by the Hepatitis B virus (HBV) [2]. Chronic HBV can be transmitted by sexual contact, through the touch, by impregnation with polluted blood, or by the direct transmission of Hepatitis B from the mother to a fetus during pregnancy (vertical transmission) [3]. According to the recent statistics of world health organization (WHO) [4], about 350 million people worldwide have been infected and carrying HBV. This serious infection is responsible for approximately 600,000 deaths each year [5]. Because of the high severity of HBV infection and a large number of deaths associated with it, it is compulsory to improve our control of this virus. Mathematical models are a vigorous tool to simulate and control the spread of the HBV infection. There exist many previous interesting works committed to studying Hepatitis B transmission. For example, Anderson and May [6] analyzed a straightforward mathematical model for illustrating the role of carrier individuals on the spread of HBV. In [7, 8], the authors developed the impact of vaccination and other controlling measures of HBV outbreak. They showed that the booster vaccine of Hepatitis B is very necessary and useful. Khan et al. [9] formulated the characteristics of HBV disease transmission and proposed the following deterministic Susceptible (S)—Infected (I)—Recovered (R) model:(1)S˙t=A−βStIt−μ+θSt,I˙t=βStIt−μ+δ+rIt,R˙t=δIt+θSt−μRt,with initial data S0=S0>0, I0=I0>0, and R0=R0>0. The positive parameters of the deterministic model (1) are given in the following list. The deterministic model constructed above can be improved by taking into account the unpredictable biological conditions [10–14]. Also, environmental fluctuations have important effects on the growth and propagation of an epidemic disease [15, 16]. Khan et al. [9] discussed the dynamics of a stochastic Hepatitis B epidemic model with varying population size. They supposed that the effect of the random fluctuations is manifested as a perturbation in the Hepatitis B transmission rate. To confer the realistic aspect to our study and make it biologically reasonable, in this study, we extend the work of Khan et al. [9] to the case of Lévy noise perturbation. We take into consideration the effects due to some unexpected and severe environmental disturbances (tsunami, floods, earthquakes, hurricanes, whirlwinds, etc.) on the disease outbreak [17, 18]. Thus, we consider the following model:(2)dSt=A−βStIt−μ+θStdt−σSt−It−dWt−∫ZηuSt−It−N˜dt,du,dIt=βStIt−μ+δ+rItdt+σSt−It−dWt+∫ZηuSt−It−N˜dt,du,dRt=δIt+θSt−μRtdt,where St− and It− are the left limits of St and It, respectively. Wt is a real-valued Brownian motion with intensity σ>0 defined on a complet probablity space Ω,ℱ,ℙ with a filtration ℱtt≥0 satisfying the usual conditions. N is a Poisson counting measure with the compensator N˜ and the characteristic measure ν on a measurable subset Z of 0,∞ satisfying νZ<∞. Wt is independent of N. We assume that ν is a Lévy measure such that N˜dt,du=Ndt,du−νdudt. The bounded function η:Z×Ω⟶ℝ is BZ×ℱt-measurable and continuous with respect to ν.In system (2), we assume that model parameters (see Table 1) are precisely known and constant. However, this hypothesis may not be validated due to the lack of data and errors of measurements. It is more realistic to study Hepatitis B dynamics with interval-valued parameters. Recently, Pal et al. [19] used interval-valued parameters to analyze the prey-predator model due to the lack of precise biological data such as prey and predator population growth rates. The same logic was applied for epidemic models. In [20], the authors treated a cholera epidemic model with uncertain parameters. They investigated the stability condition of equilibrium points. Bao et al. [21] studied a stochastic SIRS model that includes Lévy jumps and interval parameters. They established the stochastic threshold which determines the extinction and persistence in the mean of disease. In [22], the authors studied an imprecise SIR epidemic model. They solved the optimal control problem.Table 1
Biological meanings of the parameters in model (1).
Parameters
Interpretation
A
The recruitment rate corresponding to births and immigration.
μ
The natural mortality rate.
β
The transmission rate.
δ
The rate of individuals leavingI to R.
r
The disease-related death rate.
θ
The successful vaccination rate.In this paper, we consider the Hepatitis B epidemic model with stochastic transmissions and Lévy noise. To make our model more realistic, we consider imprecise biological parameters. To the best of our knowledge, the existence of a stationary distribution of system (2) with imprecise parameters remains not proved. In the next section, we propose a solution to the mentioned problem by considering an original method different from the Lyapunov approach described in [23]. Before proving the existence of unique a stationary distribution in Subsection 2.3, we demonstrate the well-posedness of the model (2) with interval-valued parameters in Subsection 2.2. Simulation examples are proposed in Subsection 2.4 to illustrate our theoretical study.


## 2. Main Results

### 2.1. Imprecise Stochastic Hepatitis B Model

Before showing the main result of this paper, we first present some definitions of interval numbers and interval-valued functions which are used in our study. Then, we construct the imprecise stochastic Hepatitis B model.Definition 1 (see [19]).
An interval numberZ is defined as Z=zˇ,z^=xzˇ≤x≤z^,x∈ℝ where ℝ is the set of all real numbers and zˇ and z^ are the lower and upper limits of the interval numbers, respectively. Furthermore, any real number z can be represented in terms of interval number as z,z.Definition 2 (see [19]).
An interval-valued function for the intervalx,y can be represented by the following function:(3)ψp=x1−pyp,forp∈0,1.Theorem 1.
The following stochastic differential equation with interval-valued parameters(4)dSt=A¯−β¯StIt−μ¯+θ¯Stdt−σ¯St−It−dWt−∫ZηuSt−It−N˜dt,du,dIt=β¯StIt−μ¯+δ¯+r¯Itdt+σ¯St−It−dWt+∫ZηuSt−It−N˜dt,du,dRt=δ¯It+θ¯St−μ¯Rtdt,where A¯∈Aˇ,A^, β¯∈βˇ,β^, μ¯∈μˇ,μ^, θ¯∈θˇ,θ^, δ¯∈δˇ,δ^, r¯∈rˇ,r^ and σ¯∈σˇ,σ^, is provided an interval-valued functional form of parameters by the following stochastic differential equation (SDE):(5)dSt=A^1−pAˇp−β^1−pβˇpStIt−μ^1−pμˇp+θ^1−pθˇpStdt−σ^1−pσˇpSt−It−dWt−∫ZηuSt−It−N˜dt,du,dIt=β^1−pβˇpStIt−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇpItdt+σ^1−pσˇpSt−It−dWt+∫ZηuSt−It−N˜dt,du,dRt=δ^1−pδˇpIt+θ^1−pθˇpSt−μ^1−pμˇpRtdt,for p∈0,1.
The proof is similar to that in [19] and hence is omitted.


### 2.2. Well-Posedness of the Stochastic Model (5)

To investigate the dynamical behavior of the model (5), the first concerning thing is whether the model is well-posed. Theorem 2 is a prerequisite for analyzing the long-run behavior of the model (5). From epidemiological considerations, it is reasonable to suppose that the intensity of Lévy jumps cannot exceed environmental carrying capacity. Hence, we impose the following standard assumption:Assumption 1.
The functionηu is bounded and A^1−pAˇp/μ^1−pμˇpηu≤Γ<1, u∈Z.
By using the Lyapunov analysis method (as mentioned in [21]), we shall verify that the solution of system (5) is global and positive.Theorem 2.
For any initial valueS0,I0,R0∈ℝ+3, there exists a unique positive solution St,It,Rt of system (5) on t≥0, and the solution will remain in ℝ+3 with probability one. That is to say, the solution St,It,Rt∈ℝ+3 for all t≥0 almost surely.Proof.
Since the coefficients of system (5) satisfy the local Lipschitz condition, then for any initial value S0,I0,R0∈ℝ+3, there is a unique local solution St,It,Rt on 0,τe, where τe is the explosion time. To show that the solution is global, we only need to prove that τe=∞ a.s. Let ϵ0>0 be sufficiently large such that S0, I0, R0 lying within the interval 1/ϵ0,ϵ0. For each integer ϵ≥ϵ0, we define the following stopping time:(6)τϵ=inft∈0,τe:minSt,It,Rt≤1ϵor maxSt,It,Rt≥ϵ,where, throughout this paper, we set inf∅=∞ (as usual, ∅ denotes the empty set). Clearly, τϵ is increasing as ϵ⟶∞. Set τ∞=limϵ⟶∞τϵ whence τ∞≤τe. If we can prove that τ∞=∞ a.s., then τe=∞ and the solution St,It,Rt∈ℝ+3 for all t≥0 almost surely. Specifically, to complete the proof, all we need is only to prove that τ∞=∞ a.s. If this statement is false, then there exists a pair of positive constants T>0 and k∈0,1 such that(7)ℙτ∞≤T>k.
Hence, there is an integerϵ1≥ϵ0 such that(8)ℙτϵ≤Tkfor allϵ≥ϵ1.
Fort≥τϵ and each ϵ,(9)dS+I+R=A^1−pAˇp−μ^1−pμˇpS+I+R−r^1−prˇpIdt≤A^1−pAˇp−μ^1−pμˇpS+I+Rdt.
Then,(10)St+It+Rt≤A^1−pAˇpμ^1−pμˇp+e−μ^1−pμˇptS0+I0+R0−A^1−pAˇpμ^1−pμˇp≤A^1−pAˇpμ^1−pμˇp,ifS0+I0+R0≤A^1−pAˇpμ^1−pμˇpS0+I0+R0,ifS0+I0+R0>A^1−pAˇpμ^1−pμˇp≔C.
Define the following LyapunovC2 function V:ℝ+3⟶ℝ+ by(11)VS,I,R=S−1−lnS+I−1−lnI+R−1−lnR.
Obviously, this function is nonnegative which can be seen fromx−1−lnx>0 for x>0.
For0≤t≤τϵ∧T, using Itô’s formula, we obtain that(12)dVS,I,R=ℒVS,I,Rdt−σ^1−pσˇpSdWt+σ^1−pσˇpIdWt−∫ZηuSI+ln1−ηuIN˜dt,du+∫ZηuSI−ln1+ηuSN˜dt,du,where ℒ is the differential operator, and(13)ℒVS,I,R=1−1SA^1−pAˇp−β^1−pβˇpSI−μ^1−pμˇp+θ^1−pθˇpS+1−1Iβ^1−pβˇpSI−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇpI+1−1Rδ^1−pδˇpI+θ^1−pθˇpS−μ^1−pμˇpR+12σ^2−2pσˇ2pI2+12σ^2−2pσˇ2pS2−∫Zln1−ηuI+ηuIνdu−∫Zln1+ηuS−ηuSνdu≤A^1−pAˇp+β^1−pβˇpC+θ^1−pθˇp+σ^2−2pσˇ2pC2+r^1−prˇp+3μ^1−pμˇp+δ^1−pδˇp+∫ZH1νdu+∫ZH2νdu,where (14)H1=−ln1−ηuI−ηuI,H2=−ln1+ηuS+ηuS.
By assumption1, we have 1−ηuI>0. In addition, by Taylor–Lagrange’s formula, we show that(15)H1=ηuI−ηuI+η2uI221−κηuI2≤Γ221−Γ2,κ∈0,1.
Similarly, we get(16)H2=−ηuS+ηuS+η2uS221+κηuS2≤Γ221−Γ2,κ∈0,1.
Therefore,(17)ℒVS,I,R≤A^1−pAˇp+β^1−pβˇpC+θ^1−pθˇp+σ^2−2pσˇ2pC2+r^1−prˇp+3μ^1−pμˇp+β^1−pβ^p+δ^1−pδ^p+θ^1−pθˇp+Γ21−Γ2νZ≔C˜,where C˜ is a positive constant. Integrating both sides of (12) from 0 to τϵ∧T, and taking expectation, we get(18)EVSτϵ∧T,Iτϵ∧T,Rτϵ∧T≤VS0,I0,R0+C˜T.
SettingΩϵ=τϵ≤T for ϵ≥ϵ0 and by (8), we have ℙΩϵ≥k. For ω∈Ωϵ, there is some component of Sτϵ, Iτϵ, and Rτϵ equals either ϵ or 1/ϵ. Hence, VSτϵ,Iτϵ,Rτϵ is not less than ϵ−1−lnϵ or 1/ϵ−1−ln1/ϵ. Consequently,(19)VS0,I0,R0+C˜T≥E1ΩϵVSτϵ,ω,Iτϵ,ω,Rτϵ,ω≥kϵ−1−lnϵ∧1ϵ−1−ln1ϵ.
Extendingϵ to ∞ leads to the contradiction. Thus, τ∞=∞ a.s. which completes the proof of the theorem.
Remark 2.5. From mathematical and biological considerations, we can study the disease dynamics of the model (5) in the following bounded set:(20)Δ=S,I,R∈ℝ+3:S+I+R≤A^1−pAˇpμ^1−pμˇpa.s..
Therefore, the regionΔ is almost surely positively invariant set by system (5).


### 2.3. Existence and Uniqueness of a Stationary Distribution to System (5)

Our aim in this subsection is to give the appropriate condition for the SDE model (5) which has a unique ergodic stationary distribution. To this end, we introduce the following lemma known as mutually exclusive possibilities. It was proved by Stettner [1].Lemma 1 (see [1]).
LetXt∈ℝn be a stochastic Feller process, then either an ergodic probability measure exists, or(21)limt⟶∞supν1t∫0tℙ∫u,X0,Σνdxdu=0,for any compact setΣ∈ℝn,where the supremum is taken over all initial distributions ν on Rd and ℙt,X0,Σ is the probability for Xt∈Σ with X0=X0∈ℝn.
For convenience, we introduce the following notation. Let(22)ℛ0s=1μ^1−pμˇp+δ^1−pδˇp+r^1−prˇpβ^A^1−pβˇAˇpμ^1−pμˇp+θ^1−pθˇp−σ˜2A^2−2pAˇ2p2μ^2−2pμˇ2p,where σ˜2=σ^2−2pσˇ2p+∫Zη2u/1−Γ2νdu.
For the ergodicity of system (5), we have the following result.Theorem 3.
Ifℛ0s>1, the stochastic system (5) admits a unique stationary distribution and it has the ergodic property for any initial value S0,I0,R0∈Δ.Proof.
The following proof is divided into three steps:Step I. Similar to the proof of Lemma 3.2 in [24] or Theorem 2.5 in [25], we briefly verify the Feller property of the SDE model (5). The main purpose of the next steps is to prove that (21) is impossible.Step II. Define(23)W=lnI+β^1−pβˇpμ^1−pμˇp+θ^1−pθˇpS.Applying Itô’s formula gives(24)dWt=β^1−pβˇpSt−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇp−12σ^2−2pσˇ2pS2t−+∫Zln1+ηuSt−−ηuSt−νdudt+σ^1−pσˇpSt−dWt+∫Zln1+ηuSt−N˜dt,du+β^A^1−pβˇAˇpμ^1−pμˇp+θ^1−pθˇpdt−β^2−2pβˇ2pμ^1−pμˇp+θ^1−pθˇpStItdt−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇpμ^1−pμˇp+θ^1−pθˇpStdt−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇpσ^1−pσˇpSt−It−dWt−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇp∫ZηuSt−It−N˜dt,du.Noting that0<S<A^1−pAˇp/μ^1−pμˇp, the equality (24) can be rewritten as follows:(25)dWt≥β^A^1−pβˇA^pμ^1−pμˇp+θ^1−pθˇp−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇp−σ^A^2−2pσˇA^2p2μ^2−2pμˇ2p−β^2−2pβˇ2pμ^1−pμˇp+θ^1−pθˇpStItdt+∫Zln1+ηuSt−−ηuSt−νdudt+σ^1−pσˇpSt−dWt+∫Zln1+ηuSt−N˜dt,du−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇpσ^1−pσˇpSt−It−dWt−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇp∫ZηuSt−It−N˜dt,du.Integrating the inequality (25) from 0 to t leads to(26)Wt−W0≥β^A^1−pβˇAˇpμ^1−pμˇp+θ^1−pθˇp−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇp−σ^A^2−2pσˇA^2p2μ^2−2pμˇ2p−β^2−2pβˇ2pμ^1−pμˇp+θ^1−pθˇp∫0tSsIsds+∫0t∫Zln1+ηuSs−−ηuSs−νduds+K1t+K2t+K3t+K4t,where(27)K1t=∫0tσ^1−pσˇpSs−dWs,K2t=−β^1−pβˇpσ^1−pσˇpμ^1−pμˇp+θ^1−pθˇp∫0tSs−Is−dWs,K3t=∫0t∫Zln1+ηuSs−N˜ds,du,K4t=−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇp∫0t∫ZηuSs−Is−N˜ds,du.The quadratic variation ofK1 is defined by K1,K1t=∫0tσ^2−2pσˇ2pS2sds. Therefore, we get(28)lim supt⟶∞K1,K1tt=σ^2−2pσˇ2plim supt⟶∞1t∫0tS2s−ds≤σ^2−2pσˇ2pA^2−2pAˇ2pμ^2−2pμˇ2p<∞a.s.Similarly, we have(29)lim supt⟶∞K2,K2tt=β^2−2pβˇ2pσ^2−2pσˇ2pμ^1−pμˇp+θ^1−pθˇp2lim supt⟶∞1t∫0tS2s−I2s−ds<∞a.s.By the assumption1, we deduce that(30)ln1−Γ≤ln1+ηuSs−≤ln1+Γ.Then(31)lim supt⟶∞K3,K3tt=lim supt⟶∞13∫0t∫Zln1+ηuSs−2νduds≤maxln1+Γ2,ln1−Γ2νZ<∞a.s.,lim supt⟶∞K4,K4tt≤β^A2−2pβAˇ2pμ^1−pμˇp+θ^1−pθˇpμ^1−pμˇp2νZ<∞a.s.According to the strong law of large numbers for local martingales [26], one can conclude that(32)limt⟶∞1tKit=0,a.s.,i=1,2,3,4.By using (16) and assumption 1, we get(33)1t∫0t∫Zln1+ηuSs−−ηuSs−νduds≥−12A^2−2pAˇ2pμ^2−2pμˇ2p∫Zη2u1−Γ2νdu.Let(34)σ˜2=σ^2−2pσˇ2p+∫Zη2u1−Γ2νdu.Therefore,(35)lim inft⟶∞1t∫0tβ^1−pβˇpSsIsds≥μ^1−pμˇp+θ^1−pθˇpβ^1−pβˇpβ^A^1−pβˇAˇpμ^1−pμˇp+θ^1−pθˇp−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇp−σ˜2A^2−2pA^2p2μ^2−2pμˇ2p.Thus, we can derive that(36)lim inft⟶∞1t∫0tβ^1−pβˇpSsIsds≥μ^1−pμˇp+θ^1−pθˇpβ^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.Step III. To continue our analysis, we need to set the following subsets:(37)Ω1=S,I,R∈ℝ+3S≥ϵ,and,I≥ϵ,Ω2=S,I,R∈ℝ+3S≤ϵ,Ω3=S,I,R∈ℝ+3I≤ϵ,whereϵ>0 is a positive constant to be determined later. It then follows from (36) that(38)liminft⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω1du≥liminft⟶+∞1t∫0tEβ^1−pβˇpSuIudu−limsupt⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω2du−limsupt⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω3du≥μ^1−pμˇp+θ^1−pθˇpβ^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1−2β^A^1−pβˇAˇpϵμ^1−pμˇp.
We can choose(39)ϵ≤μ^1−pμˇp+θ^1−pθˇpμ^1−pμˇp4β^2−2pβˇ2pA^1−pAˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1,then, we obtain(40)liminft⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω1du≥μ^1−pμˇp+θ^1−pθˇp2β^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.
Leta and b two real numbers greater than 1 such that 1/a+1/b=1. By utilizing Young inequality xy≤xa/a+yb/b for all x, y>0, we get(41)liminft⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω1du≤liminft⟶+∞1t∫0tEa−1ϖβ^1−pβˇpSuIua+b−1ϖ−b1Ω1du≤a−1ϖβ^1−pβˇpaA^1−pAˇpμ^1−pμˇp2a+liminft⟶+∞1t∫0tEb−1ϖ−b1Ω1du,where ϖ is a positive constant satisfying(42)ϖa≤a4β^1−pβˇp−a+1A^1−pAˇpμ^1−pμˇp−2aμ^1−pμˇp+θ^1−pθˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1.
From (41), we deduce that(43)liminft⟶+∞1t∫0tE1Ω1du≥μ^1−pμˇp+θ^1−pθˇpbϖb4β^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.
Setting(44)Ω4=S,I,R∈ℝ+3S≥ζ,or,I≥ζ,Σ=S,I,R∈ℝ+3ϵ≤S≤ζ,and,ϵ≤I≤ζ,where ζ>0 is a positive constant to be explained in the following. By using the Tchebychev inequality, we can observe that(45)E1Ω4≤ℙSt≥ζ+ℙIt≥ζ≤1ζESt+It≤1ζA^1−pAˇpμ^1−pμˇp.
Choosing(46)1ζ≤μ^1−pμˇp+θ^1−pθˇpbϖb8β^A^1−pβˇAˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1.
We thus obtain(47)limsupt⟶+∞1t∫0tE1Ω4du≤μ^1−pμˇp+θ^1−pθˇpbϖb8β^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1.
According to (43), one can derive that(48)liminft⟶+∞1tE∫0t1Σdu≥liminft⟶+∞1t∫0tE1Ω1du−limsupt⟶+∞1t∫0tE1Ω4du≥μ^1−pμˇp+θ^1−pθˇpbϖb8^β1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.
Based on the above analysis, we have determined a compact domainΣ⊂ℝ+3 such that(49)liminft⟶+∞1t∫0tℙu,S0,I0,R0,Σdu≥μ^1−pμˇp+θ^1−pθˇpbϖb8β^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.
Applying similar arguments to those in [24], we show the uniqueness of the ergodic stationary distribution of our model (5), denoted by π⋅. This completes the proof.


### 2.4. Numerical Simulations

In this subsection, in order to show different dynamical results of the stochastic model (2) under imprecise parameter values, we present some numerical simulations. We use Milstein’s method to simulate the trajectories of the stochastic model (5). The parameters values are given in the following list. For the purpose of showing the effects of imprecise parameters and Lévy noise on Hepatitis B dynamics, we have realized the simulation 10000 times. We assume that ηu=0.03, Z=0,∞, and νZ=1. Then, we obtain the following results: noticing that the assumption 1 is always held with parameters’ value in Table 2. From Figures 1–3, we show the existence of the unique stationary distributions for St, It, and Rt of model (5) at t=300, where the smooth curves are the probability density functions of St, It, and Rt, respectively. It can be obviously observed that the solution of the SDE model (5) persists in the mean. Furthermore, different values of the parameter imprecision p can also crucially affect the persistence of Hepatitis B (see Table 3).Table 2
Parameters’ value used in numerical simulations.
Notation
Value
Source
Notation
Value
A
^
0.4
[9]
A
ˇ
0.6
μ
^
0.09
[9]
μ
ˇ
0.2
β
^
0.1
[9]
β
ˇ
0.2
δ
^
0.3
[9]
δ
ˇ
0.5
r
^
0.1
[9]
r
ˇ
0.3
θ
^
0.2
Assumed
θ
ˇ
0.3
σ
^
0.08
Assumed
σ
ˇ
0.1Figure 1
The trajectories and histogram of solution of model (5) with initial value S0,I0,R0=0.5,0.3,0.2 and p=1.
(a)
(b)
(c)
(d)Figure 2
The trajectories and histogram of solution of model (5) with initial value S0,I0,R0=0.5,0.3,0.2 and p=0.5.
(a)
(b)
(c)
(d)Figure 3
The trajectories and histogram of solution of model (5) with initial value S0,I0,R0=0.5,0.3,0.2 and p=0.
(a)
(b)
(c)
(d)Table 3
Examples of some value ofp and their numerical illustrations.
Parameter imprecision
Value ofℛ0s
Simulation result
p
=
1
1.1192
Figure1
p
=
0.5
1.5812
Figure2
p
=
0
2.1816
Figure3


## 2.1. Imprecise Stochastic Hepatitis B Model

Before showing the main result of this paper, we first present some definitions of interval numbers and interval-valued functions which are used in our study. Then, we construct the imprecise stochastic Hepatitis B model.Definition 1 (see [19]).
An interval numberZ is defined as Z=zˇ,z^=xzˇ≤x≤z^,x∈ℝ where ℝ is the set of all real numbers and zˇ and z^ are the lower and upper limits of the interval numbers, respectively. Furthermore, any real number z can be represented in terms of interval number as z,z.Definition 2 (see [19]).
An interval-valued function for the intervalx,y can be represented by the following function:(3)ψp=x1−pyp,forp∈0,1.Theorem 1.
The following stochastic differential equation with interval-valued parameters(4)dSt=A¯−β¯StIt−μ¯+θ¯Stdt−σ¯St−It−dWt−∫ZηuSt−It−N˜dt,du,dIt=β¯StIt−μ¯+δ¯+r¯Itdt+σ¯St−It−dWt+∫ZηuSt−It−N˜dt,du,dRt=δ¯It+θ¯St−μ¯Rtdt,where A¯∈Aˇ,A^, β¯∈βˇ,β^, μ¯∈μˇ,μ^, θ¯∈θˇ,θ^, δ¯∈δˇ,δ^, r¯∈rˇ,r^ and σ¯∈σˇ,σ^, is provided an interval-valued functional form of parameters by the following stochastic differential equation (SDE):(5)dSt=A^1−pAˇp−β^1−pβˇpStIt−μ^1−pμˇp+θ^1−pθˇpStdt−σ^1−pσˇpSt−It−dWt−∫ZηuSt−It−N˜dt,du,dIt=β^1−pβˇpStIt−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇpItdt+σ^1−pσˇpSt−It−dWt+∫ZηuSt−It−N˜dt,du,dRt=δ^1−pδˇpIt+θ^1−pθˇpSt−μ^1−pμˇpRtdt,for p∈0,1.
The proof is similar to that in [19] and hence is omitted.


## 2.2. Well-Posedness of the Stochastic Model (5)

To investigate the dynamical behavior of the model (5), the first concerning thing is whether the model is well-posed. Theorem 2 is a prerequisite for analyzing the long-run behavior of the model (5). From epidemiological considerations, it is reasonable to suppose that the intensity of Lévy jumps cannot exceed environmental carrying capacity. Hence, we impose the following standard assumption:Assumption 1.
The functionηu is bounded and A^1−pAˇp/μ^1−pμˇpηu≤Γ<1, u∈Z.
By using the Lyapunov analysis method (as mentioned in [21]), we shall verify that the solution of system (5) is global and positive.Theorem 2.
For any initial valueS0,I0,R0∈ℝ+3, there exists a unique positive solution St,It,Rt of system (5) on t≥0, and the solution will remain in ℝ+3 with probability one. That is to say, the solution St,It,Rt∈ℝ+3 for all t≥0 almost surely.Proof.
Since the coefficients of system (5) satisfy the local Lipschitz condition, then for any initial value S0,I0,R0∈ℝ+3, there is a unique local solution St,It,Rt on 0,τe, where τe is the explosion time. To show that the solution is global, we only need to prove that τe=∞ a.s. Let ϵ0>0 be sufficiently large such that S0, I0, R0 lying within the interval 1/ϵ0,ϵ0. For each integer ϵ≥ϵ0, we define the following stopping time:(6)τϵ=inft∈0,τe:minSt,It,Rt≤1ϵor maxSt,It,Rt≥ϵ,where, throughout this paper, we set inf∅=∞ (as usual, ∅ denotes the empty set). Clearly, τϵ is increasing as ϵ⟶∞. Set τ∞=limϵ⟶∞τϵ whence τ∞≤τe. If we can prove that τ∞=∞ a.s., then τe=∞ and the solution St,It,Rt∈ℝ+3 for all t≥0 almost surely. Specifically, to complete the proof, all we need is only to prove that τ∞=∞ a.s. If this statement is false, then there exists a pair of positive constants T>0 and k∈0,1 such that(7)ℙτ∞≤T>k.
Hence, there is an integerϵ1≥ϵ0 such that(8)ℙτϵ≤Tkfor allϵ≥ϵ1.
Fort≥τϵ and each ϵ,(9)dS+I+R=A^1−pAˇp−μ^1−pμˇpS+I+R−r^1−prˇpIdt≤A^1−pAˇp−μ^1−pμˇpS+I+Rdt.
Then,(10)St+It+Rt≤A^1−pAˇpμ^1−pμˇp+e−μ^1−pμˇptS0+I0+R0−A^1−pAˇpμ^1−pμˇp≤A^1−pAˇpμ^1−pμˇp,ifS0+I0+R0≤A^1−pAˇpμ^1−pμˇpS0+I0+R0,ifS0+I0+R0>A^1−pAˇpμ^1−pμˇp≔C.
Define the following LyapunovC2 function V:ℝ+3⟶ℝ+ by(11)VS,I,R=S−1−lnS+I−1−lnI+R−1−lnR.
Obviously, this function is nonnegative which can be seen fromx−1−lnx>0 for x>0.
For0≤t≤τϵ∧T, using Itô’s formula, we obtain that(12)dVS,I,R=ℒVS,I,Rdt−σ^1−pσˇpSdWt+σ^1−pσˇpIdWt−∫ZηuSI+ln1−ηuIN˜dt,du+∫ZηuSI−ln1+ηuSN˜dt,du,where ℒ is the differential operator, and(13)ℒVS,I,R=1−1SA^1−pAˇp−β^1−pβˇpSI−μ^1−pμˇp+θ^1−pθˇpS+1−1Iβ^1−pβˇpSI−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇpI+1−1Rδ^1−pδˇpI+θ^1−pθˇpS−μ^1−pμˇpR+12σ^2−2pσˇ2pI2+12σ^2−2pσˇ2pS2−∫Zln1−ηuI+ηuIνdu−∫Zln1+ηuS−ηuSνdu≤A^1−pAˇp+β^1−pβˇpC+θ^1−pθˇp+σ^2−2pσˇ2pC2+r^1−prˇp+3μ^1−pμˇp+δ^1−pδˇp+∫ZH1νdu+∫ZH2νdu,where (14)H1=−ln1−ηuI−ηuI,H2=−ln1+ηuS+ηuS.
By assumption1, we have 1−ηuI>0. In addition, by Taylor–Lagrange’s formula, we show that(15)H1=ηuI−ηuI+η2uI221−κηuI2≤Γ221−Γ2,κ∈0,1.
Similarly, we get(16)H2=−ηuS+ηuS+η2uS221+κηuS2≤Γ221−Γ2,κ∈0,1.
Therefore,(17)ℒVS,I,R≤A^1−pAˇp+β^1−pβˇpC+θ^1−pθˇp+σ^2−2pσˇ2pC2+r^1−prˇp+3μ^1−pμˇp+β^1−pβ^p+δ^1−pδ^p+θ^1−pθˇp+Γ21−Γ2νZ≔C˜,where C˜ is a positive constant. Integrating both sides of (12) from 0 to τϵ∧T, and taking expectation, we get(18)EVSτϵ∧T,Iτϵ∧T,Rτϵ∧T≤VS0,I0,R0+C˜T.
SettingΩϵ=τϵ≤T for ϵ≥ϵ0 and by (8), we have ℙΩϵ≥k. For ω∈Ωϵ, there is some component of Sτϵ, Iτϵ, and Rτϵ equals either ϵ or 1/ϵ. Hence, VSτϵ,Iτϵ,Rτϵ is not less than ϵ−1−lnϵ or 1/ϵ−1−ln1/ϵ. Consequently,(19)VS0,I0,R0+C˜T≥E1ΩϵVSτϵ,ω,Iτϵ,ω,Rτϵ,ω≥kϵ−1−lnϵ∧1ϵ−1−ln1ϵ.
Extendingϵ to ∞ leads to the contradiction. Thus, τ∞=∞ a.s. which completes the proof of the theorem.
Remark 2.5. From mathematical and biological considerations, we can study the disease dynamics of the model (5) in the following bounded set:(20)Δ=S,I,R∈ℝ+3:S+I+R≤A^1−pAˇpμ^1−pμˇpa.s..
Therefore, the regionΔ is almost surely positively invariant set by system (5).


## 2.3. Existence and Uniqueness of a Stationary Distribution to System (5)

Our aim in this subsection is to give the appropriate condition for the SDE model (5) which has a unique ergodic stationary distribution. To this end, we introduce the following lemma known as mutually exclusive possibilities. It was proved by Stettner [1].Lemma 1 (see [1]).
LetXt∈ℝn be a stochastic Feller process, then either an ergodic probability measure exists, or(21)limt⟶∞supν1t∫0tℙ∫u,X0,Σνdxdu=0,for any compact setΣ∈ℝn,where the supremum is taken over all initial distributions ν on Rd and ℙt,X0,Σ is the probability for Xt∈Σ with X0=X0∈ℝn.
For convenience, we introduce the following notation. Let(22)ℛ0s=1μ^1−pμˇp+δ^1−pδˇp+r^1−prˇpβ^A^1−pβˇAˇpμ^1−pμˇp+θ^1−pθˇp−σ˜2A^2−2pAˇ2p2μ^2−2pμˇ2p,where σ˜2=σ^2−2pσˇ2p+∫Zη2u/1−Γ2νdu.
For the ergodicity of system (5), we have the following result.Theorem 3.
Ifℛ0s>1, the stochastic system (5) admits a unique stationary distribution and it has the ergodic property for any initial value S0,I0,R0∈Δ.Proof.
The following proof is divided into three steps:Step I. Similar to the proof of Lemma 3.2 in [24] or Theorem 2.5 in [25], we briefly verify the Feller property of the SDE model (5). The main purpose of the next steps is to prove that (21) is impossible.Step II. Define(23)W=lnI+β^1−pβˇpμ^1−pμˇp+θ^1−pθˇpS.Applying Itô’s formula gives(24)dWt=β^1−pβˇpSt−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇp−12σ^2−2pσˇ2pS2t−+∫Zln1+ηuSt−−ηuSt−νdudt+σ^1−pσˇpSt−dWt+∫Zln1+ηuSt−N˜dt,du+β^A^1−pβˇAˇpμ^1−pμˇp+θ^1−pθˇpdt−β^2−2pβˇ2pμ^1−pμˇp+θ^1−pθˇpStItdt−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇpμ^1−pμˇp+θ^1−pθˇpStdt−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇpσ^1−pσˇpSt−It−dWt−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇp∫ZηuSt−It−N˜dt,du.Noting that0<S<A^1−pAˇp/μ^1−pμˇp, the equality (24) can be rewritten as follows:(25)dWt≥β^A^1−pβˇA^pμ^1−pμˇp+θ^1−pθˇp−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇp−σ^A^2−2pσˇA^2p2μ^2−2pμˇ2p−β^2−2pβˇ2pμ^1−pμˇp+θ^1−pθˇpStItdt+∫Zln1+ηuSt−−ηuSt−νdudt+σ^1−pσˇpSt−dWt+∫Zln1+ηuSt−N˜dt,du−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇpσ^1−pσˇpSt−It−dWt−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇp∫ZηuSt−It−N˜dt,du.Integrating the inequality (25) from 0 to t leads to(26)Wt−W0≥β^A^1−pβˇAˇpμ^1−pμˇp+θ^1−pθˇp−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇp−σ^A^2−2pσˇA^2p2μ^2−2pμˇ2p−β^2−2pβˇ2pμ^1−pμˇp+θ^1−pθˇp∫0tSsIsds+∫0t∫Zln1+ηuSs−−ηuSs−νduds+K1t+K2t+K3t+K4t,where(27)K1t=∫0tσ^1−pσˇpSs−dWs,K2t=−β^1−pβˇpσ^1−pσˇpμ^1−pμˇp+θ^1−pθˇp∫0tSs−Is−dWs,K3t=∫0t∫Zln1+ηuSs−N˜ds,du,K4t=−β^1−pβˇpμ^1−pμˇp+θ^1−pθˇp∫0t∫ZηuSs−Is−N˜ds,du.The quadratic variation ofK1 is defined by K1,K1t=∫0tσ^2−2pσˇ2pS2sds. Therefore, we get(28)lim supt⟶∞K1,K1tt=σ^2−2pσˇ2plim supt⟶∞1t∫0tS2s−ds≤σ^2−2pσˇ2pA^2−2pAˇ2pμ^2−2pμˇ2p<∞a.s.Similarly, we have(29)lim supt⟶∞K2,K2tt=β^2−2pβˇ2pσ^2−2pσˇ2pμ^1−pμˇp+θ^1−pθˇp2lim supt⟶∞1t∫0tS2s−I2s−ds<∞a.s.By the assumption1, we deduce that(30)ln1−Γ≤ln1+ηuSs−≤ln1+Γ.Then(31)lim supt⟶∞K3,K3tt=lim supt⟶∞13∫0t∫Zln1+ηuSs−2νduds≤maxln1+Γ2,ln1−Γ2νZ<∞a.s.,lim supt⟶∞K4,K4tt≤β^A2−2pβAˇ2pμ^1−pμˇp+θ^1−pθˇpμ^1−pμˇp2νZ<∞a.s.According to the strong law of large numbers for local martingales [26], one can conclude that(32)limt⟶∞1tKit=0,a.s.,i=1,2,3,4.By using (16) and assumption 1, we get(33)1t∫0t∫Zln1+ηuSs−−ηuSs−νduds≥−12A^2−2pAˇ2pμ^2−2pμˇ2p∫Zη2u1−Γ2νdu.Let(34)σ˜2=σ^2−2pσˇ2p+∫Zη2u1−Γ2νdu.Therefore,(35)lim inft⟶∞1t∫0tβ^1−pβˇpSsIsds≥μ^1−pμˇp+θ^1−pθˇpβ^1−pβˇpβ^A^1−pβˇAˇpμ^1−pμˇp+θ^1−pθˇp−μ^1−pμˇp+δ^1−pδˇp+r^1−prˇp−σ˜2A^2−2pA^2p2μ^2−2pμˇ2p.Thus, we can derive that(36)lim inft⟶∞1t∫0tβ^1−pβˇpSsIsds≥μ^1−pμˇp+θ^1−pθˇpβ^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.Step III. To continue our analysis, we need to set the following subsets:(37)Ω1=S,I,R∈ℝ+3S≥ϵ,and,I≥ϵ,Ω2=S,I,R∈ℝ+3S≤ϵ,Ω3=S,I,R∈ℝ+3I≤ϵ,whereϵ>0 is a positive constant to be determined later. It then follows from (36) that(38)liminft⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω1du≥liminft⟶+∞1t∫0tEβ^1−pβˇpSuIudu−limsupt⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω2du−limsupt⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω3du≥μ^1−pμˇp+θ^1−pθˇpβ^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1−2β^A^1−pβˇAˇpϵμ^1−pμˇp.
We can choose(39)ϵ≤μ^1−pμˇp+θ^1−pθˇpμ^1−pμˇp4β^2−2pβˇ2pA^1−pAˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1,then, we obtain(40)liminft⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω1du≥μ^1−pμˇp+θ^1−pθˇp2β^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.
Leta and b two real numbers greater than 1 such that 1/a+1/b=1. By utilizing Young inequality xy≤xa/a+yb/b for all x, y>0, we get(41)liminft⟶+∞1t∫0tEβ^1−pβˇpSuIu1Ω1du≤liminft⟶+∞1t∫0tEa−1ϖβ^1−pβˇpSuIua+b−1ϖ−b1Ω1du≤a−1ϖβ^1−pβˇpaA^1−pAˇpμ^1−pμˇp2a+liminft⟶+∞1t∫0tEb−1ϖ−b1Ω1du,where ϖ is a positive constant satisfying(42)ϖa≤a4β^1−pβˇp−a+1A^1−pAˇpμ^1−pμˇp−2aμ^1−pμˇp+θ^1−pθˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1.
From (41), we deduce that(43)liminft⟶+∞1t∫0tE1Ω1du≥μ^1−pμˇp+θ^1−pθˇpbϖb4β^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.
Setting(44)Ω4=S,I,R∈ℝ+3S≥ζ,or,I≥ζ,Σ=S,I,R∈ℝ+3ϵ≤S≤ζ,and,ϵ≤I≤ζ,where ζ>0 is a positive constant to be explained in the following. By using the Tchebychev inequality, we can observe that(45)E1Ω4≤ℙSt≥ζ+ℙIt≥ζ≤1ζESt+It≤1ζA^1−pAˇpμ^1−pμˇp.
Choosing(46)1ζ≤μ^1−pμˇp+θ^1−pθˇpbϖb8β^A^1−pβˇAˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1.
We thus obtain(47)limsupt⟶+∞1t∫0tE1Ω4du≤μ^1−pμˇp+θ^1−pθˇpbϖb8β^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1.
According to (43), one can derive that(48)liminft⟶+∞1tE∫0t1Σdu≥liminft⟶+∞1t∫0tE1Ω1du−limsupt⟶+∞1t∫0tE1Ω4du≥μ^1−pμˇp+θ^1−pθˇpbϖb8^β1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.
Based on the above analysis, we have determined a compact domainΣ⊂ℝ+3 such that(49)liminft⟶+∞1t∫0tℙu,S0,I0,R0,Σdu≥μ^1−pμˇp+θ^1−pθˇpbϖb8β^1−pβˇpμ^1−pμˇp+δ^1−pδˇp+r^1−prˇpℛ0s−1>0a.s.
Applying similar arguments to those in [24], we show the uniqueness of the ergodic stationary distribution of our model (5), denoted by π⋅. This completes the proof.


## 2.4. Numerical Simulations

In this subsection, in order to show different dynamical results of the stochastic model (2) under imprecise parameter values, we present some numerical simulations. We use Milstein’s method to simulate the trajectories of the stochastic model (5). The parameters values are given in the following list. For the purpose of showing the effects of imprecise parameters and Lévy noise on Hepatitis B dynamics, we have realized the simulation 10000 times. We assume that ηu=0.03, Z=0,∞, and νZ=1. Then, we obtain the following results: noticing that the assumption 1 is always held with parameters’ value in Table 2. From Figures 1–3, we show the existence of the unique stationary distributions for St, It, and Rt of model (5) at t=300, where the smooth curves are the probability density functions of St, It, and Rt, respectively. It can be obviously observed that the solution of the SDE model (5) persists in the mean. Furthermore, different values of the parameter imprecision p can also crucially affect the persistence of Hepatitis B (see Table 3).Table 2
Parameters’ value used in numerical simulations.
Notation
Value
Source
Notation
Value
A
^
0.4
[9]
A
ˇ
0.6
μ
^
0.09
[9]
μ
ˇ
0.2
β
^
0.1
[9]
β
ˇ
0.2
δ
^
0.3
[9]
δ
ˇ
0.5
r
^
0.1
[9]
r
ˇ
0.3
θ
^
0.2
Assumed
θ
ˇ
0.3
σ
^
0.08
Assumed
σ
ˇ
0.1Figure 1
The trajectories and histogram of solution of model (5) with initial value S0,I0,R0=0.5,0.3,0.2 and p=1.
(a)
(b)
(c)
(d)Figure 2
The trajectories and histogram of solution of model (5) with initial value S0,I0,R0=0.5,0.3,0.2 and p=0.5.
(a)
(b)
(c)
(d)Figure 3
The trajectories and histogram of solution of model (5) with initial value S0,I0,R0=0.5,0.3,0.2 and p=0.
(a)
(b)
(c)
(d)Table 3
Examples of some value ofp and their numerical illustrations.
Parameter imprecision
Value ofℛ0s
Simulation result
p
=
1
1.1192
Figure1
p
=
0.5
1.5812
Figure2
p
=
0
2.1816
Figure3


## 3. Discussion

In the study of the dynamics of stochastic systems, the existence of an ergodic stationary distribution is one of the most important and significant characteristics. For this purpose, we have used the Feller property and mutually exclusive possibilities lemma to establish the sharp and optimal condition for the existence of the stationary distribution without employing the classical Lyapunov method. To ensure the realistic aspect of our model, we replaced constant parameters in the model (2) by imprecise ones. Based on Theorem 4.2 in [23], for any π-integrable function g:ℝ+⟶ℝ,(50)ℙlimt⟶+∞1t∫0tgXsds=∫ℝ+gxπxdx=1.The ergodic property for HBV means that the stochastic model has a unique stationary distribution which predicts the survival of the infected population in the future. That means the HBV persists for all time regardless of the initial conditions [27]. Furthermore, the ergodic property grants a reason why the integral average of a solution of system (5) converges to a fixed point whilst the system may fluctuate around as time goes by.



---

*Source: 9676501-2020-01-22.xml*