Dynamics of a stochastic HBV infection model with cell-to-cell transmission and immune response

: In this paper, considering the proven role of exosomes and the inevitable randomization within-host, we establish a hepatitis B virus (HBV) model with cell-to-cell transmission and CTL immune response from a deterministic framework to a stochastic di ﬀ erential equation (SDE). By introducing the reproduction number R 0 , we prove that R 0 can be used to govern the stochastic dynamics of the SDE HBV model. Under certain assumptions, if R 0 ≤ 1, the solution of the SDE model always ﬂuctuates around the infection-free equilibrium of the deterministic model, which indicates that the HBV will eventually disappear almost surely; if R 0 > 1, under extra conditions, the solution of the SDE model ﬂuctuates around endemic equilibrium of the corresponding deterministic model, which leads to the stochastic persistence of the HBV with probability one. One of the most interesting ﬁndings is that the ﬂuctuation amplitude is positively related to the intensity of the white noise, which can provide us some useful control strategies to regulate HBV infection dynamics.


Introduction
Hepatitis B, including acute and chronic disease, is caused by the hepatitis B virus (HBV) infection and has become a major global health problem. World Health Organization (WHO) estimates that in 2019, 325 million people globally live with hepatitis B and/or C, while in 2015, 257 million people were living with chronic hepatitis B infection (defined as hepatitis B surface antigen positive), of them, 887,000 deaths, mostly from cirrhosis and hepatocellular carcinoma (i.e., primary liver cancer) [1].
Generally, cellular immunopathological reaction plays an important role in controlling HBV infection. For example, cytotoxic T lymphocytes (CTLs) can specifically attack the target infected host cells [2], which will induce hepatocyte damage. And 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 [3][4][5][6][7] and the references therein.
Studies have shown that human immunodeficiency virus (HIV) and hepatitis C virus (HCV) can spread by two fundamental modes within-host, one by virus-to-cell infection through the extracellular space and the other by cell-to-cell transfer involving direct cell-to-cell contact [4,[8][9][10][11][12][13][14][15][16]. Especially, Sourissea et al. [17] had reported that viral transfer in vitro via cell-to-cell contact is much more rapid and efficient than infection by free virus because it avoids several biophysical and kinetic barriers. Monel et al. [18] showed that direct cell-to-cell transmission in vivo is also more potent and more efficient. And Yang and co-authors [19] showed a previously unappreciated role of exosomes in HBV transmission and natural killer cells dysfunction during chronic hepatitis B (CHB) infection. These results showed that cell-to-cell transmission is reasonable apart from virus-to-cell mode in HBV infection because exosomes can transfer genetic material between cells. Therefore, during HBV infection, uninfected hepatocytes can be infected not only by newly released free virus, but also by contacting with infected hepatocytes. Combining with the CTLs population, to describe HBV from a single patient's point of view, we first propose the following HBV infection model: where x, y, v and z denote the total numbers of uninfected hepatocytes, infected hepatocytes, free virus and CTLs, respectively. It is assumed that normal hepatocytes are produced by bone marrow and other organs at a constant rate of λ 1 and the natural death rate is d 1 x. β 1 is the effective contact rate between uninfected hepatocytes and virus, β 2 stands for the effective contact rate between uninfected and infected hepatocytes, a represents the natural death rate of the infected hepatocytes. Infected hepatocytes are eliminated by CTLs at rate pyz, free virus are produced from infected cells at rate ky and die at rate γv. CTLs are produced at a constant λ 2 from the thymus, and at rate qyz due to the stimulation of infected cells [20], and die at rate d 2 z. Note that 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. Thus, in addition to the traditional ordinary differential equations (ODEs), such as model (1), more attention has been paid to the stochastic differential equations (SDEs) which takes Brownian motion into consideration [21][22][23][24][25][26][27][28][29][30][31] and the references therein). Concretely, [31] was about HIV infection and did not consider cell-to-cell transmission. Especially, references [5,26,32] investigated the stochastic dynamics of HBV infection. After assuming the total number of hepatocytes is constant, Bertacchi et al. [5] proposed a stochastic HBV model for the infection within a patient which was treated with two drugs, either sequentially or contemporaneously, developed a two-step mutation which was resistant to both drugs, and studied the deterministic approximation of the stochastic model and gave a biological interpretation of its asymptotic behaviour. Luzyanina and Bocharov [26] considered a reduced ODEs model which only describes the interaction between the HBV and the CTL response, and Xie et al. [32] investigated the stochastic HBV dynamics. Unfortunately, these studies had not considered the CTL response and cell-to-cell transmission in their SDE HBV model.
Thanks to the insightful work of Luzyanina and Bocharov [26], we know that the variations can affect either the replication of the virus or its elimination kinetics, e.g., via parameters β 2 , γ and d 2 , respectively. And we can explore their impacts on the dynamics of HBV infection could be the extension of the deterministic description of the virus-CTL interaction to include the stochastic forcing in a multiplicative way.
To investigate the influence of fluctuating (vs. constant) rates β 2 , γ and d 2 on the model solutions, following [26], we randomize these parameters as follows.
Let p ∈ [p 1 , p 2 ] be a parameter being randomized and p 1 , p 2 are its low and upper bounds. We assume that p varies randomly according to p(t) =p + σpξ(t), wherep is the value of p around which we randomize, ξ(t) is a standard Gaussian random variable for each t and σp > 0 is the intensity of the noise. We adopt σp as σp = min(p − p 1 , p 2 −p) 3 to ensure that p remains in the interval [p 1 , p 2 ] with probability 0.997. This implies that about 99.7% of values drawn from a normal distribution are within three standard deviations σp away from the mean.
More precisely, let then, corresponding to the deterministic model (1), we can further derive the following SDE HBV model: where σ i (i = 1, 2, 3) represent the intensities of the white noises, and B i (t) (i = 1, 2, 3) are independent standard Brownian motions. It is worthy to note that, if the number of individuals in the population gets very large, Kurtz [33,34] proved that if there is a density-dependent stochastic jump process (such as the number of viral particles, the infected/uninfected hepatocytes etc..), then the process that we obtain when rescaling by a constant and large quantity can be approximated by the solution of the corresponding ODE.
The purpose of the present study is to elaborate the virus infection dynamics to random perturbations in the virus replication and immune responses parameters. In the next section, we give the global dynamics of the deterministic HBV model (1), including the basic reproduction number and global stability of infection-free equilibrium and infection equilibrium. In Section 3, we show the existence and uniqueness of the global positive solutions of the SDE HBV model (2). In Section 4, we give the asymptotic property of the positive solution of model (2) around infection-free equilibrium. And in Section 5, the dynamics of the stochastic model around infection equilibrium is obtained. In Section 6, numerical simulations of model (2) illustrate the correctness of the theoretical results and the correlation between fluctuation amplitude and intensity of the white noise is further studied. In Section 7, we discuss our new findings in the view of epidemiological implications.

Global threshold dynamics of the deterministic model (1)
We are interested in the dynamics of viral infection rather than the initial processes of infection. The initial condition of model (1) is assumed as the form x(0) > 0, y(0) > 0, v(0) > 0 and z(0) > 0. Based on the initial conditions, it is clear that the solutions of model (1) are non-negative and ultimately bounded.

Model (1) always has an infection-free equilibrium
According to the definition and algorithm of the basic reproduction number of virus in [35], we can obtain the basic reproduction number of model (1) as follows: Obviously, R 0 increases as β 1 and β 2 increasing.
After some algebraic operations, it is easy to know that, if R 0 > 1, model (1) has an infection equilibrium E * = (x * , y * , v * , z * ), where and y * is the positive root of Eq (4) in 0, (R 0 − 1)d 1 β : On the global stability of the equilibria of model (1), we have the following results.
Proof. We first prove the global stability of the disease-free equilibrium E 0 .
Define a Lyapunov function where m will be determined below.
Taking the time derivative of L 1 (x, y, v, z) with respect to the solution of model (1), we can obtain: Then Hence, M 1 = {E 0 }. Therefore it follows from Lyapunov-Lasalle invariance principle that E 0 is globally asymptotically stable when R 0 ≤ 1.
Next, we prove the global stability of the positive equilibrium E * . Consider the following Lyapunov function: At the infection equilibrium E * , we can get: Taking the time derivation of L 2 (x, y, v, z) along with the solution of model (1), and considering Eq (5), we can obtain: .
Since the arithmetical mean m+n 2 of m and n is greater than or equal to their geometrical mean So the largest invariant set of model (1) on the set M * is the singleton {E * }. By the Lasalle principle, when R 0 > 1, the infection equilibrium E * is globally asymptotically stable.
The proof is completed.
Remark 1. From Theorem 2.1, we know that R 0 can be used to govern the threshold dynamics of deterministic HBV model (1). That is, if R 0 ≤ 1, the HBV infection will go to extinction; while if R 0 > 1, the HBV infection will persist.

The existence and uniqueness of the solution of the SDE model (2)
Theorem 3.1. For any initial value X(0) = (x(0), y(0), v(0), z(0)) ∈ R 4 + , there exist positive constants N 1 and N 2 such that the set Ω is almost surely positively invariant of the SDE HBV model (2), where Proof. It is clear that all solutions of model (2) starting from R 4 From the third equation of model (2), we have where and Clearly, A 1 (t) and C 1 (t) are continuous adapted increasing process on t ≥ 0 with A 1 (0) = C 1 (0) = 0. By [36, Theorem 1.3.9, p.14], we can obtain From the first two and the fourth equations, we have It follows that where and Clearly, A 2 (t) and C 2 (t) are continuous adapted increasing process on t ≥ 0 with A 2 (0) = C 2 (0) = 0. By [36, Theorem 1.3.9, p.14], we can obtain s., which means that there is a positive constant N 2 such that Summarize the discussions above, we can conclude that the invariant set of model (2) is given as follows: The proof is completed. (2) has a unique globally positive solution X(t) = (x(t), y(t), v(t), z(t)), t ≥ 0, and the solution will remain in R 4 + with probability one.
Proof. It is obvious that the coefficients of model (2) satisfy the local Lipschitz conditions. Then, there is a unique local solution X(t) on (t ∈ [0, τ e )) for any given initial value X(0) ∈ R 4 + , where τ e is the explosion time. In order to show this solution is global, we next prove τ e = ∞ a.s.

Stochastic extinction of the HBV infection
In this section, in order to investigate whether the HBV infection is stochastically extinct, we study the dynamics around E 0 .
and R 0 ≤ 1, then the solution of model (2) has the following property: where (2) can be rewritten as: Let By the Itô's formula, we get: Noting that then we can obtain:

Consider a new Lyapunov function
where n is positive constant, we compute When R 0 ≤ 1, we get We choose n satisfying Then, Let So we obtain lim sup Considering the definition of ξ in Eq (8), we obtain the conclusion (7). This ends the proof.
Remark 2. From Theorem 4.1, we know that, when R 0 ≤ 1, if the white noise intensity is small and satisfies the solution of the stochastic model (2) is fluctuating around the disease-free equilibrium point E 0 of model (1). Epidemiologically speaking, the HBV infection dies out with probability one.

Stochastic persistence of the HBV infection
In this section, we discuss the stochastic dynamics around the positive equilibrium E * of the ODE HBV model (1) when R 0 > 1.
and R 0 > 1, then the solution of model (2) has the following property: lim sup Proof. Let By using Itô's formula, we can obtain: Consider the following Lyapunov function By Itô's formula, when σ 2 2 < γ, we can obtain:

Consider a new Lyapunov function
Because a 2 ≤ 2(a − b) 2 + 2b 2 , we can get: Let and we can derive that As a result, By integrating both sides from 0 to t and then taking the expectation on both sides, it yields Then, we can get And hence, there is lim sup Remark 3. From Theorem 5.1, we know that, when R 0 > 1, if the noise intensity is small and satisfies then the solution of the stochastic model (2) fluctuates around the infection equilibrium E * of the deterministic model (1). Epidemiologically, the HBV infection persist with probability one.

Numerical simulations
In this section, we continue to study the SDE model (2) through a numerical approach. We aim to investigate how environmental fluctuations affect disease spreading on the HBV dynamics.
As an example, the parameters are taken as follows: Then we obtain and the corresponding deterministic model (1) has a disease-free equilibrium E 0 = (1.2, 0, 0, 2.0) and a unique epidemic equilibrium E * = (1.140, 0.025, 0.015, 2.077), which is globally asymptotically stable according to Theorem 2.1.
Next, we use different values of σ i (i = 1, 2, 3) to understand the role of the noise strength in the resulting HBV dynamics for the SDE model (2).

Example 1:
We adopt (σ 1 , σ 2 , σ 3 ) = (0.15, 0.08, 0.1) and (0.2, 0.1, 0.15), respectively, then simple computations show that σ 2 1 < 0.056, σ 2 2 < 0.100, σ 2 3 < 0.050, hence the conditions in Theorem 5.1 hold. From Theorem 5.1, we know that the HBV persist with probability one. The numerical results are shown in Figures 1 and 2, we know that, after some initial transients the solutions x(t), y(t), v(t), z(t) of the SDE model (2) fluctuation around the endemic equilibrium E * = (1.140, 0.025, 0.015, 2.077) of the deterministic model (1). Comparing Figures 1 and 2, we can see that, if the noise is lower, the amplitude of fluctuation is small ( Figure 1); while if the noise is higher, the amplitude of fluctuation is big (Figure 2). For the sake of learning the impact of random perturbation on the HBV disease dynamics, we have repeated the simulation 10000 times keeping all parameters fixed and never observed any extinction scenario up to t = 100. And we give the histogram of the approximate stationary distributions of x(t) (a), y(t) (b), v(t) (c) and z(t) (d) at t = 100 for the stochastic model (2) with σ 1 = 0.15, σ 2 = 0.08 and σ 3 = 0.1 ( Figure 3) and σ 1 = 0.2, σ 2 = 0.1 and σ 3 = 0.15 (Figure 4), respectively (the numerical method can be found in [23]). One can find that x(t), y(t), v(t) and z(t) distributes normally, and also the solution to the SDE model (2) suggests for lower σ i (i = 1, 2, 3), the amplitude of fluctuation is slight and the oscillations are more symmetrically distributed ( Figure 3); while for higher σ i that the amplitude of fluctuation is remarkable and the distribution of the solution is skewed (Figure 4).
Example 2: If we adopt (σ 1 , σ 2 , σ 3 ) = (0.8, 0.65, 0.6), then σ 2 1 > 0.056, σ 2 2 > 0.100, σ 2 3 > 0.050, obviously, the conditions of Theorem 5.1 are not satisfied. The numerical results are shown in Figure 5, which suggest that the HBV infection dies out with probability one. That is, in the case of R 0 > 1, if the environmental forcing intensity σ i (i = 1, 2, 3) is so large that the conditions of Theorem 5.1 are not satisfied, the HBV could die out with probability one. Obviously, the extinction of the HBV is induced by the noise.
Form Example 2, we know that in the case of R 0 > 1, there exists extinction of the HBV infection in the SDE model (2). Next, we will focus on the extinction time of the HBV infection in the stochastic model (2).    parameters as in Figure 5, we can calculate the average extinction time of y(t) with σ 1 , and the results are shown in Figure 6. For example, when σ 1 = 0.4, the average extinction time for y(t) is 196.6; when σ 1 = 0.8, it is 79.6. We can conclude that the average extinction time decreases with the increase of noise intensity σ 1 . Finally, we will investigate the effect of the effective contact rate β 2 on the HBV dynamics of the SDE HBV model (2) (11). In this case, we can obtain and the values of R 0 are given in Table 1

Conclusions
For the study of an HBV infection model, much attention is usually paid to the persistence and extinction of HBV infection. In the deterministic HBV model, the value of the basic reproductive number is usually the key to determining the persistent and extinction of the HBV. When R 0 ≤ 1, the HBV will die out; when R 0 > 1, the HBV will persist (see Theorem 2.1). But the SDE HBV model, such as model (2), does not have any equilibrium, so we have to study the asymptotic behavior of the solution of the SDE model around the equilibrium of the deterministic HBV model to reflect the trend of the HBV infection to some extent.
In the present paper, we investigate the stochastic dynamics of an HBV infection model with cell-to-cell transmission and CTL immune response, and analyze the stochastic fluctuation of positive solutions of SDE model (2) around the infection-free equilibrium E 0 (Figures 7(a)-7(c)) and the infection equilibrium E * (Figures 1, 2 and Figure 7(d)-7(f)) of the deterministic model (1), respectively, which reveals that the influence of the white noise on HBV dynamics. More precisely, when the deterministic model (1) is disturbed by small white noise σ i (i = 1, 2, 3), the HBV infection is extinct (see Theorem 4.1 and Remark 2) or persist (see Theorem 5.1 and Remark 3), and the solutions of model (2) will fluctuate around the infection equilibrium E * of the deterministic model (1). Furthermore, via numerical simulations, we find that the fluctuation amplitude is positively related to the intensity of the white noise, which means, epidemiologically, as the intensity of white noise decreases, the fluctuation of the solution of model (2) around the equilibrium of the corresponding deterministic model (1) becomes smaller and smaller (Figures 1-4). In these cases, to a great extent, we can ignore the noise and use the deterministic model (2) to describe the HBV infection dynamics. However, when the noise σ i (i = 1, 2, 3) is sufficiently large, the infected hepatocytes y(t) on relatively large disturbance of white noise goes to extinction, and the uninfected hepatocytes x(t) generates persistence ( Figure 5). In these cases, we cannot ignore the effect of noise, that is, we should use the stochastic HBV model (2) rather than the deterministic model (1) to describe the HBV infection dynamics. In addition, it should be pointed out that the increase in β 2 results in the corresponding increase in R 0 . From the biological point of view, one can find that under certain forcing intensity, reducing the value of the effective contact rate β 2 so that the value of R 0 is less than unity could lead to the stochastically extinction of the HBV infection. These results can provide us with some useful control strategies to regulate HBV infection dynamics.