Dynamics of a stochastic HBV infection model with drug therapy and immune response

1.   Introduction

Hepatitis B virus (HBV) infection has become a worldwide problem and has been widely concerned ^{[1,2,3,4,5,6]}. WHO estimated that in 2019,296 million people were living with chronic HBV infection that may develop into cirrhosis and hepatocellular carcinoma. WHO indicated that treatment with medicines (such as interferon, nucleoside and nucleotide analogs) can slow the progression of cirrhosis, reduce incidence of liver cancer and improve long term survival. In addition, WHO estimated that in 2021, 12 to 25% of chronic HBV infection needed treatment ^[7]. However, in view of the side effects and cost factors, some hepatitis B patients often terminate antiviral treatment ^[8]. Therefore, in order to make full use of drugs to treat hepatitis B, it is very necessary to develop more effective treatment strategies ^[9].

In recent years, some mathematical models have been formulated to analyze the inhibitory effect of drug treatment on the epidemic ^{[10,11,12,13]}. Tsiang et al. proposed a HBV model based on the classical three-dimensional model established by Nowak in ^[14], in which they assumed that drugs could inhibit the production of the virus ^[15]. Furthermore, infected hepatocytes have the ability to recover, which is caused by a non cytolytic process ^[16,17]. Salman et al. considered a fractional-order model for HBV infection with cure of infected cells and discussed locally asymptotic stability of fixed point ^[18]. Lewin et al. studied a new HBV model and showed that combination therapy is more effective than single therapy for the treatment of HBV infection ^[19]. Generally speaking, T cells can proliferate and differentiate into cytotoxic effector T cells after receiving antigen stimulation, so as to rupture the target cells ^[20,21]. Yang et al. proposed a HBV fractional order model with cellular immune response ^[22]. Considering the effect of cytotoxic T lymphocytes (CTLs), Yosyingyong et al. ^[23] established the following HBV model with some treatments and immune response:

where $x$, $y$, $v$ and $z$ represent the total number of uninfected hepatocytes, infected hepatocytes, free virus and CTLs, respectively. The uninfected hepatocytes are produced by the constant rate $\Lambda$. $\beta$ is the effective contact rate between uninfected hepatocytes and virus. $u_1$ is the efficacy of the drug in blocking new infection such that the infection rate under this drug is $(1-u_1)\beta xv$. The infected hepatocytes may recover to uninfected hepatocytes at a constant rate $p$ by non-cytolytic cure process. The natural death rate of uninfected and infected hepatocytes is $\sigma$. Infected hepatocytes are eliminated by CTLs at rate $q$. Free virus are produced from infected cells at rate $m$ and die at rate $\mu$. $u_2$ is the efficacy of the drug in inhibiting viral production such that the virions production rate under this drug is $(1-u_2)my.$ CTLs are produced at a constant rate $k$ from the thymus, and at rate $s$ due to the stimulation of infected cells, and die at rate $\varepsilon$.

For model (1), the basic reproduction number of the virus is

and the authors obtained the following results ^[23]:

(A1) If $R_0 < 1$, $E_0 = \left(\Lambda/\sigma, 0, 0, k/\varepsilon\right)$ is globally asymptotically stable;

(A2) If $R_0 > 1$, $E^* = \left(x^*, y^*, v^*, z^*\right)$ is globally asymptotically stable, where

and $y^*$ is a positive root of equation

On the other hand, HBV can inevitably be disturbed by uncertain factors in the process of transmission; such as white noise, which is induced by the disturbance of parameters (e.g., the infection rate, recovery rate and so on) ^{[24,25,26,27]}. Therefore, stochastic differential equation (SDE) models are widely used by researchers in the field of biomathematics ^{[28,29,30,31,32]}. Dalal et al. established a stochastic HIV infection model with highly active antiretroviral therapy and carried out an analysis on the asymptotic behaviour of solution of the model ^[33]. Mahrouf et al. studied two stochastic viral infection models with general incidence rate: one is the three-dimensional virus model with cure rate ^[34]; the other is the four-dimensional virus model with immune response represented by CTL cells ^[35]. Xie et al. presented a stochastic HBV infection model by assuming that the infection rate fluctuates around some average value due to continuous fluctuation in the environment ^[36]. Rihan et al. investigated the impact of high-order stochastic perturbations on the dynamics of delay differential model of HBV infection with immune system ^[37]. Ma et al. formulated a stochastic viral infection model with the two modes of transmission and immune impairment and gave sufficient conditions for the extinction and persistence of the disease ^[25]. Fatini et al. considered a stochastic viral infection model to describe the role of lytic and nonlytic immune responses and studied the near optimal problem of drug therapy in blocking new infection to reduce the infection of susceptible host cells ^[38]. Unfortunately, no attempt has been made to investigate the dynamics of a stochastic HBV model with both drug therapy in blocking new infection and in inhibiting viral production, and CTL immune response.

As we already know, there are few studies on the effects of the white noise on an HBV model with antiviral therapies and CTL immune response ^{[39,40,41,42]}. The aim of this paper is to bring together stochastic differential equation and the HBV model to investigate how the environmental fluctuations affect extinction and persistence of the HBV infection, so as to promote the study of the dynamic behavior of the stochastic HBV model.

In this paper, we suppose that the random disturbances mainly affect the infection rate $\beta$ in model (1). Let $\beta\rightarrow \beta+\delta \mathrm{d}B(t)$, then corresponding to deterministic model (1), we can derive the stochastic HBV model by the following stochastic differential equations:

where $\delta$ represents the intensity of environmental forcing, $B(t)$ the standard Brownian motion defined on a complete filtered probability space $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq0}, \mathbb{P})$ with the filtration $\{\mathcal{F}_t\}_{t\geq0}$ satisfying the usual condition (i.e., it is increasing and right continuous while $\mathcal{F}_0$ contains all $\mathbb{P}$-null sets).

For the existence and uniqueness of the solution of SDE model (2), we can obtain the following results:

Theorem 1. For any initial value $(x(0), y(0), v(0), z(0))\in \mathbb{R}^{4}_{+} = \{(x, y, v, z):x(0) > 0, y(0) > 0, v(0) > 0, z(0) > 0\}$, SDE model (2) has a unique global positive solution $X(t) = (x(t), y(t), v(t), z(t)), t\geq0$, and the solution will remain in $\mathbb{R}^{4}_{+}$ with probability one.

The proof of Theorem 1 can be shown by using the similar method given in ^[36] and hence it is left to readers.

The paper is organized as follows. In Section 2, we study the stability of the disease-free equilibrium of model (2). And in the next Section, the dynamics of stochastic model (2) around the endemic equilibrium of model (1) is obtained. Section 4 is devoted to illustrating our analytical results by using some numerical examples. In Section 5, we provide a brief discussion and the summary of the main results.

2.   Stochastic extinction of the HBV infection

Consider a $d$-dimensional stochastic differential equation

Assume $f(0, t) = 0$ and $g(0, t) = 0$ for all $t\geq t_0$, and Eq (3) has the solution $x(t)\equiv0$ corresponding to the initial value $x(t_0) = 0$, which is called the trivial solution or equilibrium position.

It is clear that $E_0 = \left(\frac{\Lambda}{\sigma}, 0, 0, \frac{k}{\varepsilon}\right)$ is the disease-free equilibrium of models (1) and (2). As shown in (A1), if $R_0 < 1$, the disease-free equilibrium $E_0$ of model (1) is globally stable. In this section, we will show that the trivial solution $E_0$ of model (2) is stochastically asymptotically stable in the large.

Lemma 1. ^[43] If there exists a positive-definite decrescent radially unbounded function $V(x, t)\in C^{2, 1}(\mathbb{R}^d\times [t_0, \infty); \mathbb{R}_+)$ such that $LV(x, t)$ is negative-definite, then the trivial solution of Eq (3) is stochastically asymptotically stable in the large.

Theorem 2. Let $X(t)$ be the solution of model (2) with any initial value $X(0)\in \Gamma$. If $R_0 < 1,$ then the disease-free equilibrium $E_0 = \left(\frac{\Lambda}{\sigma}, 0, 0, \frac{k}{\varepsilon}\right)$ of model (2) is stochastically asymptotically stable in the large. Here, $\Gamma$ is the invariant set of model (2):

Proof. Set $T = x-\dfrac{\Lambda}{\sigma}, Q = z-\dfrac{k}{\varepsilon},$ and model (2) can be rewritten as:

Let

where $c$ is a positive constant to be determined.

By the Itô's formula, we get

If $R_0 < 1$, then

Choosing $c = \dfrac{(1-u_1)\beta\Lambda}{(p+\sigma+\frac{qk}{\varepsilon})\mu}$, we have

Consider a new Lyapunov function

By the Itô's formula, we get

Let $V(T, y, v, Q)\in C^{2, 1}(\mathbb{R}^4; \mathbb{R}_+)$ such that $V(T, y, v, Q) = V_{1}(T, y, v, Q)+V_{2}(T, y, v, Q)$, which is positive-definite decrescent radially-unbounded function, then

which is negative-definite. Hence, it follows from Lemma 1 that the disease-free equilibrium $E_0$ of model (2) is stochastically asymptotically stable in the large.

Remark 1. In Theorem 2, if $R_0 < 1$, no matter how strong the white noise is, $E_0$ of model (2) is stochastically asymptotically stable in the large, which means that the HBV infection will be eliminated with probability one. Moreover, for extremely large noise intensity, eradication of HBV infection is faster in stochastic HBV model (2) than in deterministic HBV model (1).

3.   Stochastic persistence of the HBV infection

As mentioned in (A2), when $R_0 > 1$, model (1) has a unique endemic equilibrium $E^* = (x^*, y^*, v^*, z^*),$ which is globally stable. However, since SDE model (2) has no endemic equilibrium, we will study the asymptotic behavior of the solution of model (2) around $E^*$ to reflect whether HBV infection is persistent to some extent.

Let

and it follows from (4) that $x(t)\leq W, y(t)\leq W, v(t)\leq W, z(t)\leq W.$ Then, we provide the following theorem about the stochastic persistence of the HBV infection.

Theorem 3. For any given initial value $X(0) = (x(0), y(0), v(0), z(0))\in \Gamma$, if $R_0 > 1$ and

hold, the solution of model (2) has the following property:

where,

Proof. Let

Using Itô's formula, we obtain

Consider a new Lyapunov function

By Itô's formula, we obtain

Let

By Itô's formula, we have

Let

By Itô's formula, we get

Define $V = V_{1}+n_1V_{2}+n_2V_3+n_3V_4,$ and

Choose $n_1 = \dfrac{2(1-u_1)^2\beta^2\Lambda^2}{\sigma^3\mu}, n_2 = \dfrac{p\mu+n_1(1-u_2)^2m^2}{\sigma\mu}, n_3 = \dfrac{2n_2Wz^*q^2(\varepsilon+\sigma)^2}{\sigma ks^2},$ and we can obtain

where $\eta: = \min\left\{\dfrac{3\sigma}{4}-\dfrac{\delta^2(1-u_1)^2(1-u_2)^2m^2\Lambda^2}{\sigma^2\mu^2}, \dfrac{3n_2\sigma}{4}, \dfrac{n_1\mu}{4}, \dfrac{n_2\varepsilon q^2}{s^2}\right\}.$

As a result,

By integrating both sides from 0 to $t$ and then taking the expectation on both sides, it yields

where $M: = \dfrac{\delta^2(1-u_1)^2(1-u_2)^2m^2\Lambda^2(x^*)^2}{\sigma^2\mu^2}+n_3s(W^2+y^*z^*).$

Hence,

Remark 2. In Theorem 3, in the case of $R_0 > 1$, if the white noise intensity is small enough to satisfy $\delta^2 < \Delta,$ the solution of stochastic model (2) is fluctuating around the endemic equilibrium $E^*$ of model (1), which reveals that the HBV infection will persist with probability one.

4.   Numerical simulations

In this section, with the help of MATLAB, we provide numerical simulations to show the effects of environmental fluctuation and drug therapies on the HBV dynamics from the analytic results. Through employing the Milstein's method in ^[44], we obtain the discretization equations of model (2):

where $\xi_i, i = 1, 2, \ldots, n$ are independent Gaussian random variables $N(0, 1)$. And we take the parameter values as in ^[23]:

4.1.   The impact of environmental noise on the HBV dynamics

From Theorem 2 and Remark 1, we can see that when $R_0 < 1$, the white noise has no effect on the dynamic behavior of model (2). Therefore, we will mainly focus on the impact of noise intensity on the resulting dynamics of model (2) in the case of $R_0 > 1$.

Fix $u_1 = u_2 = 0.5$, and we can easily calculate that

We adopt $\delta = 1\times10^{-10}$, then simple computations show that $\delta^2 = 1\times10^{-20} < \Delta$. According to Theorem 3, we can conclude that the HBV infection persists with probability one in this case. The numerical simulations are shown in Figure 1, which suggest that the solutions $(x(t), y(t), v(t), z(t))$ of stochastic HBV model (2) will fluctuate around endemic equilibrium $E^*$ of deterministic HBV model (1). If we adopt $\delta = 5\times10^{-10}, 7\times10^{-10}, 9\times10^{-10}, 11\times10^{-10}, 13\times10^{-10}$ and $15\times10^{-10}$, then $\delta^2 > 3.5504\times10^{-20}$; obviously, the conditions of Theorem 3 are not satisfied. The numerical results are shown in Figure 2, which suggests that the solution of model (2) no longer fluctuates around the equilibrium $E^*$. In addition, though Figure 2(a)–(c) suggests that the HBV infection still persists with probability one, the abundance of infected hepatitis cells decreases gradually. As the noise intensity increases, the HBV infection will die out with probability one (see Figure 2(d)–(f)).

4.2.   The effects of drug therapies on the HBV dynamics

We adopt the values of $(u_1, u_2)$ in Table 1. Clearly, when we take $\delta = 2\times10^{-9}$, the conditions of Theorem 3 are not satisfied. We have repeated the simulation 10,000 times keeping all parameters fixed and given the average extinction time of $v(t)$ when $\delta = 0$ and $\delta = 2\times10^{-9}$. See Figure 3.

On the one hand, for stochastic model (2) with $\delta = 2\times10^{-9}$, we find that the effects of drug treatment on the average extinction time of virus are complex. In more detail, when only one drug therapy is adopted, i.e., $u_1 = 0$ or $u_2 = 0$, the numerical simulations shown in Figure 3(a), (b) suggests that as the efficacy of the drug increases, the average extinction time of $v(t)$ increases at first and then decreases. In addition, if drug therapies both in blocking new infection and in inhibiting viral production are adopted, there are two cases: one is that when one of them is not very effective, i.e., $u_1 = 0.4$ or $u_2 = 0.4$, the average extinction time of $v(t)$ is still increasing at first and then decreasing as the efficacy of the another drug increases (see Figure 3(c), (d)); the other is that when one of them is significantly effective, i.e., $u_1 = 0.8$ or $u_2 = 0.8$, the average extinction time of $v(t)$ is decreasing as the efficacy of the another drug increases (see Figure 3(e), (f)). On the other, for stochastic model (2) with $\delta = 0$, i.e., deterministic model (1), as the efficacy of a drug increases, the average extinction time of $v(t)$ always decreases regardless of the efficacy of another drug. Moreover, the examples shed the interesting light that environmental fluctuation shortens the average extinction time of virus to some extent.

Next, we take $\delta = 2\times10^{-9}$ as an example to discuss the effects of two drug therapies on the clearance of HBV infection for stochastic model (2). From Figure 3, when $(u_1, u_2) = (0.8, 0.8)$, the average extinction time for $v(t)$ is 62.8657; when $(u_1, u_2) = (0.8, 0)$, the average extinction time for $v(t)$ is 65.7663; when $(u_1, u_2) = (0, 0.8)$, the average extinction time for $v(t)$ is 63.4311. We can assert that combination therapy is more effective than single therapy to some extent. Moreover, when $(u_1, u_2) = (0.4, 0.8)$, the average extinction time for $v(t)$ is 63.1429; when $(u_1, u_2) = (0.8, 0.4)$, the average extinction time for $v(t)$ is 64.7367. We can conclude that inhibition of virus production is partly more effective than blocking new infection.

5.   Conclusions

In this paper, we proposed a stochastic HBV infection model with drug therapies and CTL immune response to study the effect of environmental randomness on the virus dynamics, in which we assume that infected hepatocytes can be cured by non-cytolytic process and the infection rate $\beta$ is disturbed by random factors. And the value of this study lies in two aspects.

Mathematically, on the one hand, we prove that when $R_0 < 1$, the solution $E_0$ of model (2) is stochastically asymptotically stable in the large, which reveals that HBV infection dies out with probability one no matter how strong the noise is (Theorem 2 and Remark 1). That is to say, if HBV infection is eradicated for deterministic HBV model (1), environmental fluctuation will not cause the outbreak of HBV infection. On the other, we show that when the noise intensity is small enough and $R_0 > 1$, the solution of random model (2) fluctuates around endemic equilibrium $E^*$ of deterministic model (1), which reveals that HBV infection persists with probability one (Theorem 3 and Remark 2).

Epidemiologically, HBV infection is eliminated for both stochastic model (2) and deterministic model (1) when $R_0 < 1$, so we partially provide the effects of the environment fluctuations on the dynamics of stochastic HBV model (2) when $R_0 > 1$. If the noise intensity is relatively small, the HBV infection persists and the dynamics of random model (2) is consistent with that of deterministic model (1) (Figures 1 and 2(a)–(c)), which reveal that we can ignore the noise and use deterministic model (1) to describe the HBV infection dynamics to a great extent. It is noteworthy that in the case of $R_0 > 1$, bigger intensity of noise $\delta$ may lead to elimination of HBV infection almost surely (Figure 2(d)–(f)). Therefore, it is necessary to take stronger random fluctuations in order to effectively control the outbreak of HBV infection. This indicates that we cannot ignore the influence of random factors; that is, it is better to use SDE HBV model (2) rather than deterministic model (1) to describe the HBV infection dynamics.

In addition, unlike the method of simulating paths of the solutions of the proposed models used in ^[23,33,45], we study the effects of drug therapies on the dynamics of HBV infection through calculating the average extinction time of the solution $v(t)$ for model (2) (see Figure 3). Our numerical examples demonstrate rather strikingly that environmental noise can effectively shorten the extinction time of virus. What's more, the therapeutic effect of blocking new infection is not as significant as that of inhibiting virus production to a certain extent. Furthermore, in some sense, the combination of drugs that inhibit the production of virus and drugs that block new infection is more conducive to the elimination of HBV infection than single drugs. These findings may have some potential clinical value during treatment.

It should be noted that this work is based on the assumption that the infection of hepatocytes, the production of virus and the activation of immune response occur instantaneously. However, it usually takes a certain period of time for the virus to infect healthy hepatocytes until the new virus is produced and released from infected hepatocytes, or for the body to receive antigen stimulation and produce CTLs. Therefore, it is necessary to introduce time delay into stochastic HBV model, which is also our future work.

Acknowledgements

The authors would like to thank the anonymous referees for very helpful suggestions and comments which led to improvements of our original manuscript. The authors would like to thank the editor and the referees for their helpful comments. Y. Tan, Y. Cai and W. Wang were supported by the National Natural Science Foundation of China (Grant numbers 12071173 and 1217011789), the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (20KJB110025) and the Huaian Key Laboratory for Infectious Diseases Control and Prevention (HAP201704). Z. Peng was supported by the National Natural Science Foundation of China (Grant number 82073673), the National S & T Major Project Foundation of China (2018ZX10715002) and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). K. Wang was supported by the National Natural Science Foundation of China (Grant number 12171396). R. Yao was supported by the Natural Science Basic Research Program of Shaanxi Province, China (2021JZ-21).

Conflict of interest

The authors declare that they have no competing interests.