Analysis of a stochastic HBV infection model with delayed immune response

1.   Introduction

Hepatitis B is a liver infection, caused by the hepatitis B virus (HBV), is responsible for more than 400 million chronic infections worldwide ^[1]. HBV is a leading cause of broad-spectrum liver diseases such as hepatitis, cirrhosis and liver cancer ^[2]. Some people with HBV are sick for only a few weeks (known as 'acute' infection), but for others, the disease progresses to a serious, lifelong illness known as 'chronic' hepatitis B. Majority of infected adults successfully clear the virus and acquires lifelong immunity ^[3]. The immune response to HBV-encoded antigens is responsible both for viral clearance and for disease pathogenesis during this infection. While the humoral antibody response to viral envelope antigens contributes to the clearance of circulating virus particles, the cellular immune response to the envelope, nucleocapsid, and polymerase antigens eliminates infected cells. During acute HBV infection, cytotoxic T lymphocytes (CTLs) can directly attack infected hepatocytes and participate in the pathogenesis of liver disease by orchestrating diverse components of the immune system; see ^[4,5].

In fact, the process intercellular transmission and virus-to-cell infection is not instantaneous but needs to be completed over a period of time, so it is necessary to consider the effect of time-delays on the HBV system ^[6,7,8,9]. Moreover, the parameters for growth and interactions depend on the state and nature of the virus, the condition of the immune system, and the environment in which the interaction takes place the body ^[10]. The environment of the body is determined by the overall health of the individual. One way to explore the impact of body environmental factors on the dynamics of HBV infection could be the extension of the deterministic description of the virus-CTL interaction to include the stochastic forcing either in an additive or multiplicative way. Mathematical models to investigate the dynamics of HBV transmission within environmental noise have been studied by many researchers, among them ^{[11,12,13,14]}.

Many other mathematical models have been designed to evaluate the effect of public health programs and provided long-term predictions regarding the disease prevalence and control ^{[15,16,17,18,19]}. More and more attentions have been paid to the study of virus dynamics within-host, which can provide insights into virus infection and dynamics, as well as to how an infection can be reduced or even eradicated, see ^{[20,21,22,23,24]}. However, most of these approaches are based on deterministic models and do not consider the randomness in cell transmission and effect of environmental variability.

Motivated by the mentioned biological and mathematical considerations, in the present paper, we investigate the dynamics of stochastic delay differential equations (SDDEs) of HBV model with cell-to-cell transmission and CTLs immune response. For more realistic situation of the development process of the disease, we incorporate the effect of multiple time-delays and randomization within a host. The organization of the rest of this paper is as follows: In Section 2, we propose a stochastic delay differential model for HBV infection. In Section 3, we investigate existence and uniqueness of the global positive solution, and study existence of stationary distribution in Section 4. Possible extinction of the disease is studied in Section 5. Some numerical results and simulations are provided in Section 6 to show the effectiveness of the theoretical results. Concluding remarks are discussed in the last Section.

2.   Time-delay dstochastic model of HBV

Although HBV replication at the cellular level is not fully understood, in this paper, we propose a new model of stochastic delay differential equations (SDDEs) for HBV replications in one host. We assume that during HBV infection, uninfected (healthy) hepatocytes can be infected not only by newly released free virus, but also by contacting with infected hepatocytes. We also assume that the cytotoxic T lymphocytes (CTLs) can specifically attack the target infected host cells. Of course there is an intracellular time-delay (time-lag) between the infection of a cell and the viral particles emission, and virus production. Time-delay is also required to represent incubation period, the time required for the production of new virus particles. Herein, based on the basic model of Nowak et al. ^[25], we introduce a delay differential model to combine the CTLs population with HBV infection. The model takes the form

$H(t)$, $I(t)$, $V(t)$, and $D(t)$, respectively, denote the healthy hepatocytes that are not infected by the viruses, the infected hepatocytes which are infected by viruses, hepatitis B viruses and CTLs. Time-delay $\tau_1$ is considered, in the first term of the second equation, to justify the required time between initial infection of a cell by HBV and the release of new virions. It is also incorporated in the second term to consider the reaction time that healthy hepatocytes become infected by the infected cell contacts and then transformed into the infected hepatocytes; While $\tau_2$ stands for the time necessary for the newly produced particles to become mature then infectious particles. The healthy cells become infected either by free viruses at rate $\beta_1HV$ (virus-to-cell infection mode), or by direct contact with an infected cell at rate $\beta_2HI$ (cell-to-cell transmission mode). Hence, the term $\beta_1HV+\beta_2HI$ represents the total infection rate of uninfected cells. Infected cells are eliminated by CTLs at rate $\beta_3ID$, $a$ is the production rate of free viruses by infected cells; While CTLs are produced at a constant $\eta_2$ from the thymus and at the rate $\beta_4ID$ as a result of stimulation of infected cells (see Figure 1). (Existence of the equilibrium points and basic reproduction number ${\cal R}_0$ for the deterministic system (2.1) are given in the Appendix.) The description of the model parameters is presented in Table 1.

As a matter of fact, there are inevitably random disturbances in the process of HBV infection within-host, such as temperature fluctuation, mood fluctuation and other physiological rhythm changes, which may affect the dynamics of HBV infection. Taking this into consideration enables a lot of authors to introduce randomness into deterministic model of biological systems to reveal the effect of environmental variability, see ^[26,27,28].

For more realistic situation of the development process of the disease, we incorporate the effect of randomization within-host by introducing nonlinear perturbation on the natural death rate with white noise into each equation of system (2.1). 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. Therefore, we propose a delayed stochastic model of the form

subject to the initial conditions

Here, ${\cal C}$ is the family of Lebesgue integrable functions ${\cal C}([-\tau, 0], \mathbb{R}^4_+)$. Such that, $\nu_{11}, \nu_{12}, \nu_{21}, \nu_{22}, \nu_{31}$, $\nu_{32}, \nu_{41}, \nu_{42}$, represent the intensities of the white noise and $W_i$, $(i = 1, 2, 3, 4)$ is a real-valued standard Brownian motion defined on a complete probability space $(\Omega, {\cal A}, \mathbb{P})$ satisfying the usual conditions ^[29]. We assume that the parameters $\alpha_i$, $i = 1, 2, 3, 4$, are distributed by some non linear stochastic noise ^[30]. The random perturbation may be dependent on square of the state variables $H$, $I$, $V$ and $D$ of system (2.2), respectively, that is to say $\alpha_1\rightarrow \alpha_1-(\nu_{11}H+\nu_{12})dW_1$, $\alpha_2\rightarrow \alpha_2-(\nu_{21}+\nu_{22}I)dW_2$, $\alpha_3\rightarrow \alpha_3-(\nu_{31}+\nu_{32}V)dW_3$, $\alpha_4\rightarrow \alpha_4-(\nu_{41}+\nu_{42}D)dW_4$.

3.   Existence and uniqueness of the global positive solution

In this section, we provide some conditions that guarantee a unique global positive solution of the SDDEs system (2.2). This can be achieved that if the coefficients of the system realize the growth and Lipschitzian conditions, then there will be a positive solution.

Theorem 1. For any initial value system (2.3), there is a unique positive solution $(H(t),$$I(t),$$V(t),$$D(t))$ of system (2.2), on $t\geq-\tau$ and the solution will remain in $\mathbb{R}^4_+$ almost surely (a.s.).

Proof. Since all the coefficients of system (2.2) are Lipschitz continuous, therefore, there is a unique local solution $(H(t), I(t), V(t), D(t))$ on $[-\tau, \tau_e)$, where $\tau_e$ is an explosion time. To show this solution is global, one may need to show $\tau_e = \infty$ a.s. (almost surely). Let $l_0 > 0$ be sufficiently large so that $(H(t), I(t), V(t), D(t)) = \{(\varphi_1(t), \varphi_2(t), \varphi_3(t), \varphi_4(t)):-\tau \leq t \leq 0\} \in {\cal C}([-\tau, 0];\mathbb{R}_+^4)$ all lie within the interval $[\frac{1}{l_0}, l_0]$. For each integer $l\geq l_0$, define the stopping time

let $\inf \phi = \infty$. $\tau_l$ is increasing with $l$ and let $\tau_\infty = \lim_{l\rightarrow \infty}\tau_l$, then $\tau_\infty \leq \tau_e$ and by showing $\tau_\infty = \infty$ a.s., the aim is to conclude that $\tau_e = \infty$ a.s. If this assertion is erroneous, then there exists a pair of constants $T > 0$ and $\epsilon \in (0, 1)$ such that $\mathbb{P}\{\tau_\infty\leq T\} > \epsilon$. Therefore, there is an integer $l_1\geq l_0$ such that

Define a ${\cal C}^2$-function ${\cal G}:\mathbb{R}_+^4\rightarrow \mathbb{R}_+$ as follows:

where $l_2$ is a positive constants to be determined later. By It$\hat{o}$'s formula, we have

Here,

Choosing $l_2 = \beta_3/\beta_4$, yields

where ${\cal B}$ is a positive constant. It follows that ${\cal L}{\cal G}(.)$ is bounded. Since the following proof is standard and it is similar to the method in the literature ^[31], so it is omitted. Therefore, the proof is completed.

4.   Existence of ergodic stationary distribution

Herein, we construct a suitable stochastic Lyapunov function to study existence of a unique ergodic stationary distribution of the positive solutions to system (2.2). Ergodic stationary distribution of a stochastic model is one of the most important and significant characteristics. Ergodic property of a stochastic HBV epidemic model means that the stochastic model has a unique stationary distribution which predicts the persistence of the disease in the future under some restrictions on the intensity of white noise, that is the stochastic model fluctuate in a neighborhood of the infected equilibrium, ${\cal E}^*$ (defined in the Appendix) of the corresponding undisturbed model for all time regardless of the initial conditions.

First, assume that $X(t)$ is a regular time-homogenous Markov process in $\mathbb{R}^d$, illustrated by the SDDE

The diffusion matrix of the process $X(t)$ is

Lemma 1. ^[32]. The Markov process $X(t)$ has a unique ergodic stationary distribution $\pi(.)$ if there exist a bounded domain ${\cal U}\subset \mathbb{R}^d$ with regular boundary $\Gamma$ and

(i): there is a positive number ${\cal M}$ such that $\sum_{i, j = 1}^d \lambda_{ij}(x)\xi_i\xi_j\geq {\cal M} |\xi|^2, x\in {\cal U}, \xi \in \mathbb{R}^d.$

(ii): there exists a nonnegative ${\cal C}^2$-function $V$ such that ${\cal L}V$ is negative for any $\mathbb{R}^d \setminus {\cal U}.$

Theorem 2. Assume that

where

then system (2.2) admits a unique stationary distribution $\pi(.)$ and it has the ergodic property.

Proof. In order to prove Theorem 2, it is enough to validate conditions $(i)$ and $(ii)$ of Lemma 1.

We first prove condition $(i)$. The diffusion matrix of system (2.2) is given by

Let ${\cal U}$ be any bounded domain in $\mathbb{R}_+^4$, then there exists a positive constant

such that

for any $(H, I, V, D)\in \bar{{\cal U}}_\sigma, \xi = (\xi_1, \xi_2, \xi_3, \xi_4)\in\mathbb{R}_+^4$. Thus, we have verified that condition $(i)$ of Lemma 1 is satisfied. We then prove condition $(ii)$ of Lemma 1. For any relatively small $\theta\in(0, 1)$, we define

where,

such that $f_1 = \frac{8}{3(1-\theta)f_2^\theta}, \quad f_2 = 2\sqrt[3]{\frac{\eta_1}{(1-\theta)\nu_{42}^2}}.$ Clearly, $\lim\limits_{\theta\rightarrow0^+}{\cal R}_0^s(\theta) = {\cal R}_0^s$. Since ${\cal R}_0^s(\theta)$ is continuous and ${\cal R}_0^s > 1$, we can choose relatively small $\theta$ such that ${\cal R}_0^s(\theta) > 1$. By system (2.2), we have

and

Define

where $c_1$, $c_2$, $h_1$, $h_2$, $m_1$, $m_2$, $m_3$ and $k_1, k_2, k_3$ are positive constants which will be determined later. By It$\hat{o}$ formula to ${\cal V}_1$, we obtain

Let

Therefore,

Thus, from systems (4.4) and (4.8), we have

By It$\hat{o}$'s formula to ${\cal V}_2$, we get

such that,

Let

Thus, we get

From systems (4.5), (4.7) and (4.10), we have

Applying It$\hat{o}$ formula to ${\cal V}_3$, one gets

Hence,

Let $m_1 = \frac{1}{(1-\theta)\nu_{31}^\theta}$. Therefore,

By It$\hat{o}$'s formula to ${\cal V}_4$, we have

Applying It$\hat{o}$'s formula to ${\cal V}_5$, one obtains

Additionally, we have

Let

Therefore, we obtain

Define a ${\cal C}^2$-function $\tilde{{\cal V}}:\mathbb{R}^4_+\rightarrow \mathbb{R}$ in the following form

where ${\cal M} > 0$ is a sufficiently large number satisfying the following condition

and

Noting that $\tilde{{\cal V}}(H, I, V, D)$ is not only continuous, but also tends to $+\infty$ as $(H, I, V, D)$ approches the boundary of $\mathbb{R}_+^4$, and $\Vert(H, I, V, D)\Vert\rightarrow \infty$. Hence, $\tilde{{\cal V}}$ must have a minimum point $(H_0, I_0, V_0, D_0)$ in the interior of $\mathbb{R}_+^4$. We define a ${\cal C}^2$-function ${\cal V}:\mathbb{R}_+^4\rightarrow \mathbb{R}_+$ as follows:

By It$\hat{o}$'s formula, it follows that

Now, one can construct a bounded open domain ${\cal U}_\epsilon$ such that condition $(ii)$ of Lemma 1 satisfies. Define a bounded open set, for arbitrary $\epsilon > 0$, as follows

Therefore, we need to prove ${\cal L}{\cal V}\leq-1$ for $(H, I, V, D)\in \mathbb{R}_+^4\setminus{\cal U}_\epsilon$. Clearly, $\mathbb{R}_+^4\setminus{\cal U}_\epsilon = \cup_{i = 1}^8{\cal U}_i$, such that

Choosing

Case I: For any $(H, I, V, D)\in {\cal U}_1$, by system (4.23), one obtains

Let $-\frac{\eta_1}{\epsilon}+\lambda_2\leq-1$, yields ${\cal L}{\cal V}\leq-1$.

Case II: For any $(H, I, V, D)\in{\cal U}_2$ from system (4.23), one may have

choosing $-\frac{(1-\theta)\nu_{11}^{\theta+2}}{4\epsilon^{\theta+2}}+\lambda_2\leq-1$, yields ${\cal L}{\cal V}\leq-1$.

Case III: For any $(H, I, V, D)\in{\cal U}_3$ from system (4.23), we have

we select $\epsilon$ small enough such that the following condition holds

Therefore, we have

Case V: For any $(H, I, V, D)\in{\cal U}_4$ from system (4.23), we have

Let $-\frac{a}{\epsilon}+\lambda_2\leq-1$, then we have ${\cal L}{\cal V}\leq -1$.

Case VI: For any $(H, I, V, D)\in{\cal U}_5$ from system (4.23), we have

choose $-\frac{\eta_2}{\epsilon}+\lambda_2\leq-1$, yields ${\cal L}{\cal V}\leq-1$.

Case VII: For any $(H, I, V, D)\in{\cal U}_6$ from system (4.23), we have

let $-\frac{(1-\theta)\nu_{22}^{\theta+2}}{4\epsilon^{2(\theta+2)}}+\lambda_2\leq-1$, then we obtain ${\cal L}{\cal V}\leq-1$.

Case IV: For any $(H, I, V, D)\in{\cal U}_7$ from system (4.23), we have

choose $-\frac{(1-\theta)\nu_{32}^{\theta+2}}{4\epsilon^{3(\theta+2)}}+\lambda_2\leq-1$, we obtain that ${\cal L}{\cal V}\leq-1$.

Case VI: For any $(H, I, V, D)\in{\cal U}_8$ from system (4.23), we have

let $-\frac{(1-\theta)\nu_{32}^{\theta+2}}{4\epsilon^{\theta+2}}+\lambda_2\leq-1$, then ${\cal L}{\cal V}\leq-1$. Therefore, condition $(ii)$ of Lemma 1 satisfies, such that system (2.2) identifies a unique stationary distribution $\pi(.)$.

Remark 1. If ${\cal R}_0^s > 1$ the solution of system (2.2) fluctuates around the endemic equilibrium of the undisturbed system (2.1) under certain conditions. This means that the disease will be persistent, provided that the intensities of white noise are adequately small; See Figure 2.

5.   Extinction of the disease

For the undisturbed system (2.1), if ${\cal R}_0 \le 1$, the disease-free equilibrium ${\cal E}_0$ is globally asymptotically stability and HBV infection will die out. However, for ${\cal R}_0 > 1$, ${\cal E}^*$ is globally asymptotically stable and ${\cal E}_0$ is unstable; see ^[11]. Now, we investigate the possible extinction of the disease $I(t)$ for the stochastic system (2.2). Define

Here,

Theorem 3. For any given initial value system (2.3) on $t\geq-\tau$, if ${\cal R}_0^e < 0$, the solution $(H(t), I(t), V(t), D(t))$ of system (2.2) satisfies the following:

Such that $\lim\limits_{t\rightarrow \infty}I(t) = 0$ a.s., i.e., the disease, $I(t)$, will die out exponentially with probability one and $\lim\limits_{t\rightarrow \infty}V(t) = 0$ a.s. Additionally, the distributions $H(t)$ and $D(t)$ converge weakly to the measures which have the densities $\pi(x)$ and $\pi(y)$, which will be determined later.

Proof. Assume the following auxiliary logistic equations with nonlinear stochastic perturbation

with initial value $X_0 = H_0 > 0$.

with initial value $Y_0 = D_0 > 0$. Therefore, systems (5.3) and (5.4) have the ergodic property ^[30], such that the invariant densities are given by

where ${\cal N}_1$ and ${\cal N}_2$ are constants such that

Let $X(t)$ be the solution of system (5.3) with initial value $X_0 = H_0 > 0$, by the comparison theorem ^[33], one obtains $H(t)\leq X(t)$ for any $t\geq0$ a.s. Similarly, assume that $Y(t)$ be the solution of system (5.4) with initial value $Y_0 = D_0 > 0$, one may have $D(t)\leq Y(t)$ such that $I(t) < 1 /\beta_4$.

On the other hand, by ^[34], consider the vector $(\chi_1, \chi_2) = \Big(\frac{a}{\alpha_3}, \sqrt{\frac{a\beta_1\eta_1}{\alpha_1\alpha_3(\alpha_2-\frac{\beta_2\eta_1}{\alpha_1})}}\Big)$. Therefore, we can derive that there exists a left eigenvector of

corresponding to the spectral radius of ${\cal A}_0$; $\rho({\cal A}_0) = \sqrt{\frac{a\beta_1\eta_1}{\alpha_1\alpha_3(\alpha_2-\frac{\beta_2\eta_1}{\alpha_1})}}$, which can be denoted as $\sqrt{\frac{a\beta_1\eta_1}{\alpha_1\alpha_3(\alpha_2-\frac{\beta_2\eta_1}{\alpha_1})}}\Big)(\chi_1, \chi_2) = (\chi_1, \chi_2){\cal A}_0$.

Consider a ${\cal C}^2$-function $\hat{{\cal V}}:\mathbb{R}^2_+\rightarrow \mathbb{R}_+$ by

where $\varphi_1 = \frac{\chi_1}{(\alpha_2-\frac{\beta_2\eta_1}{\alpha_1})}$ and $\varphi_2 = \frac{\chi_2}{\alpha_3}$, by It$\hat{o}$ formula to $\hat{{\cal V}}$, one may have

such that

Hence, we have

by substituting the values of $\chi_1$ and $\chi_2$ from system (5.1), we obtain

Additionally, squaring both sides of system (5.5) and by Cauchy inequality, one may have

Hence,

From systems (5.6), (5.7) and (5.8), we have

Integrating both sides of system (5.9) then dividing it by $t$, we obtain

Assume that $M_i(t)$, $i = 1, 2, 3, 4$, are real valued continuous local martingale vanishing at $t = 0$; Such that

In addition, their quadratic form are given by

By strong law of large numbers ^[35], one obtaines

Hence, we have

Let $\epsilon$, $\zeta_1$ and $T_1$ are positive constants, by the exponential martingale inequality ^[36], for each $T_1\geq 1$, such that $j = 3, 4$, we have

Let $T_1 = n$, $\epsilon = 1$ and $\zeta_1 = 2\ln n$, we get

By using Borel-Cantelli Lemma ^[36], there is $\Omega_0\subset\Omega$ with $\mathbb{P}(\Omega_0) = 1$ such that for $\rho \in \Omega_0$ there exists an integer $n_0 = n_0(\rho)$, such that

For all $0\leq t\leq n \wedge n\geq n_0(\rho)$ a.s. That is, for $0\leq n-1\leq t\leq n$, one obtains

Noting that $X(t)$ is ergodic and $\int_0^\infty x\pi(x)dx < \infty$ a.s., therefore, we have

In view of systems (5.13) and (5.18) by taking superior limit on both sides of system (5.17), one gets

If ${\cal R}_0^e < 0$, then $\lim\sup\limits_{t\rightarrow \infty} \frac{\ln I(t)}{t} < 0$, and $\lim\sup\limits_{t\rightarrow \infty} \frac{\ln V(t)}{t} < 0$, a.s. which implies that $\lim_{t\rightarrow \infty}I(t) = 0$ and $\lim\limits_{t\rightarrow \infty}V(t) = 0$ a.s. Therefore, the disease will die out exponentially with probability one.

We arrive at the following Remark.

Remark 2. From $0 < \Phi_0 < 1$, when the stochastic perturbations $\nu_{21}$ and $\nu_{31}$ are sufficiently large so that ${\cal R}_0^e < 0$, the stochastic system (2.2) displays disease extinction with probability one. Also, stochastic perturbations lead to the eradication of the disease in the stochastic system faster than the undisturbed system; See Figure 3.

6.   Numerical simulations

In this section, some numerical simulations are given to validate the theoretical results, through Milstein's Higher order method ^[41,42], to numerically solve SDDEs system (2.2). The discretization transformation takes the form

The independent Gaussian random variables denoted as $\xi_{i, k}$, $(i = 1, 2, 3, 4)$, which follow the distribution $N(0, 1)$, $m_1$ and $m_2$ are integers such that the time-delays can be expressed in terms of the step-size $\Delta t$ as $\tau_1 = m_1\Delta t$ and $\tau_2 = m_2\Delta t$. Initial values are taken fixed $(0.4, 0.2, 0.7, 0.5)$.

Example 1. Herein, we use parameter values of Table 2 to analyze the dynamics of systems (2.1) and (2.2), with time-delays $\tau_1 = 1$ and $\tau_2 = 2$, we choose $\nu_{11} = 0.003$, $\nu_{12} = 0.001$, $\nu_{21} = 0.004$, $\nu_{22} = 0.001$, $\nu_{31} = 0.001$, $\nu_{32} = 0.006$, $\nu_{41} = 0.003$, $\nu_{42} = 0.004$. Direct calculations leads to ${\cal R}_0^s > 1$, the condition of Theorem 2 satisfies. The stationary distribution illustrates that the solution of the stochastic system (2.2) fluctuate in a neighborhood of the endemic equilibrium ${\cal E}^*$ of the corresponding undisturbed system (2.1), which means that the disease is persistent for all time regardless of the initial conditions if the scale of random perturbations is relatively small. Therefore, system (2.2) admits a unique ergodic stationary distribution $\pi(.)$. The simulation also indicates that the endemic equilibrium is asymptotically stable for the undisturbed system; See Figure 2.

Example 2. Now, we compare the solution of the SDDEs system (2.2) with the undisturbed system (2.1) around the disease free equilibrium point, using parameter values of Table 2. Figure 3 shows that the disease free equilibrium point is stable for the undisturbed system. However, for the SDDEs system (2.2), we choose $\nu_{11} = 0.1$, $\nu_{12} = 0.1$, $\nu_{21} = 0.6$, $\nu_{22} = 0.6$, $\nu_{31} = 0.7$, $\nu_{32} = 0.3$, $\nu_{41} = 0.2$, $\nu_{42} = 0.1$, by a simple calculation we have $\Phi_0\approx 0.067 < 1$ and ${\cal R}_0^e < 0$, so that conditions of Theorem 3 hold, Therefore, the disease will be extinct with probability one. It is shown that the larger intensity of white noise set may help to eliminate the disease faster than the model without noise.

7.   Concluding remarks

In the present work, we investigated the impact of high-order stochastic perturbations on the dynamics of delay differential model of HBV infection with both virus-to-cell and cell-to-cell transmissions, intracellular delay, and immunity. The effect of stochastic perturbations on the persistence and possible extinction of the disease have been studied in detail. By utilizing Lyapunov functional, we proved the existence and uniqueness of an ergodic stationary distribution of positive solutions to the system, where the solution fluctuates around endemic equilibrium of the corresponding deterministic model and leads to the stochastic persistence of the disease with probability one. The model has a unique stationary distribution which is ergodic if ${\cal R}_0^s > 1$. In addition, we formulate sufficient conditions for complete extinction of the disease by constructing a suitable stochastic Lyapunov function. Under some conditions, the disease can die out exponentially with probability one. Some numerical simulations, by using Milstein's scheme, are carried out to show the effectiveness of the obtained results. The intensity of white noise plays an important role in the treatment of HBV and other infectious diseases.

The incorporation of intracellular time-delays and stochastic perturbations (noise) in the model is assumed to give a clearer view in the interpretation of the analytical result and this has important implications on therapeutic options and drug development. Some other interesting topics deserve further investigation. One may take other kinds of environmental noise into account, such as the Levy noise ^[43]. In addition, motivated by the work in ^[44,45] the deterministic system (2.1) can be extended to include fractional derivatives in the model in order to consider long-run memory of the dynamic of the disease.

Acknowledgement

This work has been generously funded by the UAE university (UAE), fund # 31S435-UPAR -5-2020. The authors believe that the editor and reviewers' suggestions have been very helpful in improving the manuscript.

Conflict of interest

The authors declare no conflict of interest in this paper.

Appendix: Equilibrium points and ${\cal R}_0$ of the deterministic model

The deterministic system (2.1) admits two equilibrium points, namely; The disease-free equilibrium, ${\cal E}_0 = (H_0, 0, 0, D_0)$, where $H_0 = \frac{\eta_1}{\alpha_1}$ and $D_0 = \frac{\eta_2}{\alpha_4}$; The infected equilibrium, ${\cal E}^* = (H^*, I^*, V^*, D^*)$, where

with $\alpha_4 > \beta_4 I^*$ such that $I^*$ is the positive root of

To determine the expression of the basic reproduction number, we utilize the next generation matrix approach ^[46]. Therefore, we have

The basic reproduction number is the spectral radius of $({\cal F}\mathbb{V}^{-1})$, i.e. ${\cal R}_0 = \rho({\cal F}\mathbb{V}^{-1})$. Hence,

where ${\cal R}_{01} = \frac{a\alpha_4\beta_1H_0}{\alpha_3(\beta_3\eta_2+\alpha_2\alpha_4)}$ and ${\cal R}_{02} = \frac{\alpha_4\beta_2H_0}{\beta_3\eta_2+\alpha_2\alpha_4}$. From biological point of view, ${\cal R}_{01}$ stands for the average number of secondary infected cells produced by an infectious virion and ${\cal R}_{02}$ represents the average number of secondary infected cells produced by an infected cell.