GLOBAL STABILITY OF A DIFFUSIVE AND DELAYED HBV INFECTION MODEL WITH HBV DNA-CONTAINING CAPSIDS AND GENERAL INCIDENCE RATE

. The aim of this paper is to study the dynamics of a new chronic HBV infection model that includes spatial diﬀusion, three time delays and a general incidence function. First, we analyze the well-posedness of the initial value problem of the model in the bounded domain. Then, we deﬁne a thresh- old parameter R 0 called the basic reproduction number and show that our model admits two possible equilibria, namely the infection-free equilibrium E 1 as well as the chronic infection equilibrium E 2 . Further, by constructing two appropriate Lyapunov functionals, we prove that E 1 is globally asymptotically stable when R 0 < 1, corresponding to the viruses are cleared and the disease dies out; if R 0 > 1, then E 1 becomes unstable and the equilibrium point E 2 appears and is globally asymptotically stable, which means that the viruses persist in the host and the infection becomes chronic. An application is pro- vided to conﬁrm the main theoretical results. Additionally, it is worth saying that, our results suggest theoretically useful method to control HBV infection and these results can be applied to a variety of possible incidence functions presented in a series of other papers.


1.
Introduction. As one of the top three infectious diseases in China, Hepatitis B can result in acute or chronic liver diseases and put people at high risk of death from cirrhosis of the liver and liver cancer. Nowadays, it is a major global health problem and is carrying an enormous economic and social burdens. According to World Health Organization (WHO), roughly 0.24 billion people are chronically infected with hepatitis B. Besides, more than 0.78 million people die every year due to complications of hepatitis B, including fatty liver, liver cirrhosis and hepatocellular carcinoma (HCC) [1]. In the last decades, a number of dynamic models with respect to HBV infection have been introduced [21,8,3,34,16,4,29], and these models have been proved to be of use in understanding pathogenesis and designing better treatments protocols. A basic HBV infection model, containing uninfected target cells (presumably hepatocytes, the same as below), infected cells and free viruses, has been analyzed quantificationally by Nowak et al. [27]. To accurately determine the dynamics of infection and clearance in three acutely infected chimpanzees, Murray et al. [26] proposed the mathematical model with all intracellular components of HBV infection. In the literature [25], Murray et al. researched a simplified version of the HBV infection model described in [26] and found that the half-life of HBV virions is approximately 4 hours. These models however do not incorporate the uninfected cells and the intracellular HBV DNA-containing capsids simultaneously. For this purpose, Manna and Chakrabarty [22] was the first to contain both uninfected cells and intracellular HBV DNA-containing capsids aside from infected cells and free viruses. The corresponding mathematical model is as follows: dD(t) dt = aI(t) − βD(t) − δD(t),

GLOBAL STABILITY OF A DIFFUSIVE AND DELAYED HBV INFECTION MODEL 4225
in which state variables H, I, D and V are the concentrations of uninfected (susceptible) hepatocytes, infected hepatocytes, intracellular HBV DNA-containing capsids and free HBVs, respectively. In the first state equation in (1), healthy hepatocytes are created at a constant rate s, either from differentiation of progenitor cells or by direct proliferation of mature hepatocytes [29]. Furthermore, µ is the death rate of healthy hepatocytes and k is a constant infection rate characterized by the infection efficiency. In the second state equation in (1), δ represents the natural death rate of infected hepatocytes. In the third state equation in (1), the parameter a indicates the production rate of intracellular HBV DNA-containing capsids per infected hepatocyte, β denotes the rate at which these capsids are transmitted to blood which then lead to the growth of free virions, and δ is per capita death rate for infected hepatocytes. The parameter c that appears in the last state equation in (1) stands for the clearance rate of virions in plasma. It is assumed that all parameters in (1) are strictly positive and µ ≤ δ [22]. Note that the forenamed model (1) assumed that the uninfected hepatocytes which are exposed to free HBVs immediately become infected and the maturation process of the capsids is instantaneous despite the fact that this two delays actually exist [21,22,36,13]. To make model (1) much closer to reality, Manna and Chakrabarty [21] incorporated the eclipse phase and capsids maturation period, which are denoted as τ 1 and τ 2 respectively. More specifically, they assumed that the uninfected hepatocytes become actively infected at time t but are infected by free HBVs at time t − τ 1 , and the infected cells create new intracellular HBV DNAcontaining capsids at time t but are penetrated by viruses at time t − τ 2 . Therefore, in the document [21], they considered the following refined model with two time delays: obtained the equilibrium solutions of system (2) and proved their global stability by constructing Liapunov functionals.
In model (2), the key assumptions are that cells and viruses are well mixed in space at all times, and the mobility of cells and viruses are ignored. Several models [37,23,11,35,40] have corrected this problem by including a spatial component and adding Fickian diffusion for the virus while assuming that hepatocytes cannot move under normal conditions. Another important feature is the fact that there is a lag between the maturation time of the intracellular HBV DNA-containing capsids and the time of the mature capsids producing free viruses, and some researchers have incorporated the exponentially decay functions during the delays into the viral dynamic models [13,23,11,40,12]. Besides, the viral infection rate in model (2) is assumed to be bilinear in the free viruses and uninfected cells, which is not very appropriate. As pointed out by Min et al. [24] and Chen et al. [2], the bilinear incidence rate indicates that a person with a smaller liver may be less prone to infection than an individual with a larger one. Because of this, the bilinear incidence rate is replaced by standard incidence in [2], by Holling type II functional response in [20], by Beddington-DeAngelis incidence in [5], by Crowley-Martin incidence in [38], by general functional response in [13]. Motivated by above three biological facts, we establish the following diffused HBV model with three time delays and a general incidence rate (excluding bilinear incidence and saturation incidence): for t > 0, x ∈ Ω, with initial conditions and homogeneous Neumann boundary conditions where Ω is a connected, bounded open set in R n with smooth boundary ∂Ω and ∂ ∂ n denotes the outward normal derivative on ∂Ω. The boundary conditions in (5) imply that the virus populations do not move across the boundary ∂Ω. In biological terms, H(x, t), I(x, t), D(x, t) and V (x, t) denote the densities of uninfected cells, infected cells, intracellular HBV DNA-containing capsids and free viruses at position x and time t, respectively. The diffusion coefficient is indicated by d v and ∆ = ∂ 2 ∂x 2 is the Laplacian operator. The third delay τ 3 represents the time necessary for the newly produced HBV DNA-containing capsids to become free viruses and we assume the surviving rates for the infected cells from time t−τ 1 to time t, the immature capsids from time t − τ 2 to time t, as well as the immature free viruses from time t − τ 3 to time t obey exponentially decay functions. The biological meanings of all other parameters are the same as those given in (2). As in [23,11,14], the general incidence function f (H, I, V ) is assumed to be continuously differentiable in the interior of R 3 + and satisfies the following conditions: In the light of the above description, we know that the novelty of this model is that it concerns the general incidence function, the diffusion for free HBVs, survival rates and three delays. More precisely, firstly, in our model, the general incidence function (eliminating bilinear incidence and saturation incidence), instead of mass action response, is used to describe the growth rate of infected cells. In this case, this new model not only generalizes HBV models containing a variety of famous incidence forms [24,2,20,5,38] but also removes the disadvantage of the bilinear rate in (2). Secondly, on account of spatially uniform models are not sufficient to give a realistic picture of disease diffusion, the spatial effects and Fickian diffusion for the virus are introduced. Meanwhile, we assume that the movement of cells is slow enough that it can be ignored. Thirdly, our research is more complicated, because it consists of three probabilities of surviving and one delay except two delays appear in model (2). Moreover, the effect of multi-time delays is considered as an important factor, which will be close to reality.
The outline of this paper is as follows. In Section 2, firstly, we prove that our model is well posed by showing the existence, non-negativity and boundedness of solutions for model (3)-(5). After, we determine the basic reproduction number and existence of all equilibria. In Section 3, we discuss the local and global stability of the infection-free equilibrium E 1 . In Section 4, by constructing appropriate Lyapunov functional, we prove the global stability of the chronic infection equilibrium E 2 . In Section 5, an application is given to confirm our analytical theory. Finally, a brief conclusion and discussion are given in Section 6.
2. Basic results. In this section, we establish the existence, non-negativity and boundedness of solutions of model (3)-(5) because this model describes the evolution of liver cells and free viruses. Furthermore, we give the basic reproduction number R 0 and existence conditions of the positive homogeneous equilibria of the model.

2.1.
Global existence and well-posedness of solutions. Before proceeding we shall put forward some notations that will be used in the following analysis. Let X = BC U (Ω, R 4 ) be the set of all bounded and uniformly continuous functions fromΩ to R 4 equipped with a norm · X defined by φ X = sup x∈Ω |φ(x)|, where | · | denotes the Euclidean norm on R 4 . In this case, we introduce symbol Additionally, for any given a continuous function Proof. For any x ∈Ω and φ = (φ 1 , φ 2 , φ 3 , φ 4 ) T ∈ C, we define F = (F 1 , F 2 , F 3 , F 4 ) : C → X as follow:

TING GUO, HAIHONG LIU, CHENGLIN XU AND FANG YAN
It is easy to see that F is locally Lipschitz in C. Then we reformulate model (3)-(5) as the following abstract functional differential equation: According to standard existence theory [33,6], it is not hard to deduce that there exists a unique local solution for model (6) on [0, T max ), where T max is the maximal existence time for solutions of model (6). In addition, the inequalities H(x, t) ≥ 0, Next, we show that the solutions of (3)-(5) are bounded. Firstly, we define a new variable: From the first three equations of (3), we get 0)}}, and thus G(x, t) ≤ M , which implies that G(x, t) is bounded and so are H, I and D.
Then we can prove that V (x, t) is bounded. Using the boundedness of D and model (3)-(5), we obtain the following system IfṼ (t) is a solution to the ordinary differential equation By the comparison principle [28], we get V (x, t) ≤Ṽ (t). Hence, The above analysis supports the conclusion that H(x, t), I(x, t), D(x, t) and V (x, t) are bounded onΩ × [0, T max ). Therefore, we deduce that T max = +∞ from the standard theory for semilinear parabolic systems [15]. This completes the proof.

2.2.
Basic reproduction number and existence of positive equilibria. In this subsection, we show that model (3)-(5) has two possible equilibria. Moreover, existence of these equilibria is determined by a threshold parameter R 0 , which is given by Here, R 0 is called the basic reproduction number and describes the average number of new infected cells derived from one infected cell when the free HBVs have just entered the body. In the expression (7), 1 δ is the average survival time of an infectious cell. During this period a virus-producing cell generates a intracellular HBV DNA-containing capsids per unit time. β β+δ gives the amount of free viruses produced from an intracellular HBV DNA-containing capsid during its survival period. e −α1τ1 , e −α2τ2 and e −α3τ3 denote the probabilities of surviving the infected cells from time t − τ 1 to time t, the immature capsids from time t − τ 2 to time t, as well as the immature free viruses from time t − τ 3 to time t, respectively. 1 c represents the average life expectancy of a virus and f s µ , 0, 0 represents the value of the function f at the beginning of the infection process. By multiplying the above quantities together, we get the expected number of newly infected cells generated by one infected cell, that is R 0 . Furthermore, it is important to note the following remark.
Remark 1. If no delays are considered (τ 1 = τ 2 = τ 3 = 0) or the mortalities during the three delays are ignored (α 1 = α 2 = α 3 = 0), our R 0 coincides with the basic reproduction number of model (2). Whereas, in presence of the three delays, R 0 is a decreasing function of the mortalities. This implies that, any one of the mortalities during the three delays can decrease R 0 . Hence, ignoring the mortality during the delay in a viral model will overestimate R 0 . In other words, our R 0 is biologically well defined.
To simplify the analysis, we show only the existence conditions of positive homogeneous equilibria of model (3)-(5) and have the following theorem: Proof. Obviously, model (3)-(5) always has an infection-free equilibrium E 1 , which represents the state that the cellular infection initiated by HBVs will ultimately die out. To find the positive equilibrium, we study the following system: A short calculation gives V = βe −α3τ3 c D, D = ae −α2τ2 β + δ I and s − µH − δe α1τ1 I = 0.

TING GUO, HAIHONG LIU, CHENGLIN XU AND FANG YAN
This means that we can get the following equation: In order to guarantee I = s−µH δe α 1 τ 1 ≥ 0, we must have H ≤ s µ . Thus there is not equilibrium point if H > s µ . Now, we redefine the function g on interval 0, s µ as follows: From (9), it is easy to show that and Thus, when R 0 > 1, there exists at least one positive equilibrium point E 2 = (H 2 , I 2 , D 2 , V 2 ) ∈ R 4 > 0. Next, we show that E 2 is the unique chronic infection equilibrium of system (3)-(5). Using hypothesis (T 1), we have This imply g is strictly increasing in the interior of the feasible region. Combining (10) and (11), it follows that if R 0 > 1, there exists a unique chronic infection equilibrium E 2 with H 2 ∈ 0, s µ and I 2 , D 2 , V 2 > 0.
3. Stability of the infection-free equilibrium E 1 . The purpose of this section is to give a rigorous investigation for the local and global stability of the infectionfree equilibrium E 1 . First, we analyze the local asymptotic stability of E 1 . To do so, we need to determine the characteristic equation about this point. Let 0 = η 1 < η 2 < · · · < η n < · · · be the eigenvalues of the operator −∆ on Ω with homogeneous Neumann boundary conditions, and denote the eigenfunction space corresponding to η i in C 1 (Ω) by E(η i ) for i = 1, 2, · · · . Let X = [C 1 (Ω)] 4 , φ ij : j = 1, · · · , dim E(η i ) be an orthonormal basis of E(η i ) and X ij = cφ ij : c ∈ R 4 . Then Suppose E * = (H * , I * , D * , V * ) be any feasible steady state of system (3)-(5), E 1 or E 2 , and consider the following change of variables: Then by substituting Z 1 (x, t), Z 2 (x, t), Z 3 (x, t) and Z 4 (x, t) into model (3) and linearizing, we obtain a new system of the form This new system can be equivalently expressed by Thus, λ is an eigenvalue if and only if the matrix −λI−Qη i +B +Ce −λτ1 +Le −λτ2 + N e −λτ3 has determinant zero for all i ≥ 1, here I denotes 4 × 4 identity matrix. In addition, we get the characteristic equation at equilibrium point E * as follow: In the following theorem, we will obtain the local stability of the steady state E 1 .

GLOBAL STABILITY OF A DIFFUSIVE AND DELAYED HBV INFECTION MODEL 4233
Squaring and adding the two equations of (16) gives Letting ψ = ω 2 yields Clearly, all real roots of (17) are negative provided R 0 < 1. Therefore, we conclude that E 1 is locally asymptotically stable for any time delay τ 1 , τ 2 , τ 3 ≥ 0 when R 0 < 1, completing the proof.
Theorem 3.1 only establishes local stability of infection-free equilibrium E 1 . However, the global stability of equilibrium is very useful in researching the fundamental question of whether this equilibrium be induced ultimately. So, in the next content, we focus on the mathematical analysis of the global dynamics of E 1 . Moreover, for the global stability of E 1 , we have the following theorem: is globally asymptotically stable.
Proof. Define the following Lyapunov functional: dξ + e α1τ1 I(x, t) + δ a e α1τ1+α2τ2 D(x, t) where H 1 = s µ . Obviously, the summation of the first three terms in the right-hand side of U is a non-negative number. Indeed, if H(x, t) ≥ H 1 , then with the aid of hypothesis (T 1). If H(x, t) < H 1 holds, apply similar reasoning to the above formula, we have the same conclusion that the functional U is nonnegative. For convenience, we will use the following notations: w = w(x, t) and w τi = w(x, t − τ i ), i = 1, 2, 3 for any w ∈ H, I, D, V . Calculating the time derivative of U along solutions of system (3)-(5), we obtain According to the Divergence theorem and the homogeneous Neumann boundary condition (5), we have Besides, we use hypothesis (T 1) again to note that since the function f (H, I, V ) is strictly monotonically increasing with respect to H, it follows that Thus, we have that if R 0 < 1, dU dt ≤ 0 for all H, I, D, V ≥ 0. Furthermore, it is easy to see that for dU dt = 0, then V = 0 and H = s µ hold. When the conditions V = 0 and H = s µ are satisfied, combined with system (3)-(5), we have D = I = 0. That is to say, the largest compact invariant set in (H, I, D, V ) ∈ R 4 + : dU dt = 0 is the singleton E 1 . From LaSalle invariance principle [39], we conclude that the infection-free equilibrium E 1 of system (3)-(5) is globally asymptotically stable when R 0 < 1.

4.
Global stability of the chronic infection equilibrium E 2 . According to the above analysis, we know E 1 becomes unstable and the chronic infection equilibrium E 2 emerges when R 0 > 1. Thus, in this section, we discuss the global stability of E 2 by constructing a Lyapunov functional based on the Volterra function It is obvious that the function R attains its strict global minimum at 1 and satisfies R(1) = 0. And in the next theorem, we will make use of the following further assumption about f :  Proof. Consider the following Lyapunov functional:

GLOBAL STABILITY OF A DIFFUSIVE AND DELAYED HBV INFECTION MODEL 4237
On the basis of the assumption (T 2), we have Using the Neumann boundary condition (5) and the Divergence theorem, we obtain Ω ∆V dx = 0 and Further R(z) ≥ 0 for any z > 0 implies that dL dt ≤ 0 with equality if and only if H = H 2 , I = I 2 , D = D 2 and V = V 2 . Thus, the largest compact invariant set in (H, I, D, V ) dL dt = 0 is the singleton E 2 . By LaSalle invariance principle [10], we show that the chronic infection equilibrium E 2 is globally asymptotically stable when R 0 > 1.
Based upon the above analysis, we know that the extinction and persistence of viral infections crucially depend on the basic reproduction number R 0 . What is more, it is evident that R 0 is a decreasing function on death rates α 1 , α 2 and α 3 . Hence, the neglect of death rates must result in increase in size of R 0 , as shown in model (2). In Fig. 5, we remark that R 0 becomes large enough when death rates α 1 , α 2 and α 3 tend to 0, which confirms the result in Remark 1.
On the other hand, when the basic reproduction number R 0 is less than one, the viruses are cleared and the infection dies out (the case when E 1 is globally asymptotically stable). Moreover, by the precise expression of R 0 , we find that we can reduce R 0 to below one by increasing delays τ 1 , τ 2 and τ 3 . Thus, a good strategy to control the HBV should focus on any drugs that can prolong the values of three delays. Numerical simulation for the relationship between R 0 and three delays are shown in Fig. 6.  6. Conclusion and discussion. In this paper, based on the fact of HBV infection [19,25,22] and motivated by the works [13,23,24,2], we assume that only the viruses move freely in liver, and then establish a 4-dimensional diffusion HBV model with three delays and a general incidence rate. For this mathematical model, we define the basic reproduction number R 0 that acts as a threshold to predict whether the disease persist in the host. When the general incidence function f is assumed to meet biologically reasonable conditions (T 1) and (T 2), we discussed the global stability of the infection-free equilibrium E 1 and the chronic infection equilibrium E 2 by constructing suitable Lyapunov functionals and using LaSalle invariance principle. More precisely, we have shown that E 1 is globally asymptotically stable whenever R 0 < 1. In this case, all positive solutions converge to E 1 Figure 6. The plots of the basic reproduction number R 0 as a function of three delays τ 1 , τ 2 and τ 3 . Here, s = 2.6×10 7 , µ = 0.01, δ = 0.053, a = 150, β = 0.87, c = 3.8, b 1 = 0.01, α 1 = 0.2, α 2 = 0.28, α 3 = 0.1 and k = 2.4 × 10 −3 . (a) τ 3 = 4, (b) τ 2 = 6, and (c) τ 1 = 5.8. and the disease can be eradicated ultimately. When R 0 > 1, E 1 becomes unstable and there appears a chronic infection equilibrium E 2 which is globally asymptotically stable. In this case, all positive solutions converge to E 2 and the disease will be persistent in the host. Compared with model (2), these results imply that the diffusion of free viruses and time delays have no effect on the global stability of the HBV infection model under homogeneous Neumann boundary conditions. On the other hand, notice that R 0 is a decreasing function of the three delays. Hence, we can reduce the value of R 0 to a level lower than one by increasing three delays in an effort to prevent the viruses. Moreover, ignoring the mortalities during the three delays and the third delay in model (3) will overestimate R 0 . Thus, this study is more realistic and may help to analyze the dynamical behavior of other models including commonly used incidence rates [2,20,5,38].