SPATIOTEMPORAL DYNAMICS OF A FRACTIONAL MODEL FOR HEPATITIS B VIRUS INFECTION WITH CELLULAR IMMUNITY

In this paper, we propose and investigate a fractional diffusive model for hepatitis B virus (HBV) infection with capsids and immune response presented by cytotoxic T lymphocyte (CTL) cells. We derive the conditions for global asymptotic stability of the steady states of the model in terms of the basic reproduction number R0 and the immune response reproduction number R1. By constructing appropriate Lyapunov functionals, it is shown that the infection-free equilibrium is globally asymptotically stable when R0 ≤ 1, the immune-free infection equilibrium is globally asymptotically stable when R1 ≤ 1 < R0 and the infection equilibrium with CTL immune response is globally asymptotically stable when R1 > 1. Numerical simulations are performed to illustrate the analytical results. Mathematics Subject Classification. 26A33, 37B25, 37N25, 92D30. Received May 25, 2020. Accepted December 19, 2020.


Introduction
Many mathematical models used fractional differential equations (FDEs) have been proposed and analyzed to gain insights into the dynamics of HBV infection in vivo. In [32], Zhou and Sun proposed and investigated the dynamics of a fractional-order HBV infection model with the standard incidence function instead of the mass action term for the infection process. Salman and Youssef [27] also considered a fractional-order model for HBV infection with cure of infected cells and discussed the local asymptotic stability of equilibria. Shi et al. [28] in their paper constructed and analyzed the global dynamics of a fractional-order model with Holling II fuctional response to describe the transmission of HBV. Further, Bachraoui et al. [2] proposed a fractional order-model for HBV infection with capsids and CTL immune response that improved and generalized the mathematical models formulated by ordinary differential equations (ODEs) in [19,20] and also the FDE models introduced in [6,27,32] by considering the Hattaf's incidence rate [11] that includes the common types such as the bilinear incidence rate, the saturated incidence rate and the Beddington-De Anglis functional response [4,9].
Mathematical modeling of HBV infection using the partial differential equations (PDEs) has attracted the attention of many researchers. Wang and Wang [30] proposed an HBV model with spatial dependence. They considered the random mobility of virions and extended the basic virus infection model proposed by Nowak et al. [23]. They also neglected the spatial mobility for the susceptible and infected hepatocytes. The movement of virus was adopted via Fickian diffusion, where the population flux of virus is proportional to the concentration gradient with negative proportionality constant. Geng et al. [10] considered the mobility of capsids and viruses and applied the nonstandard finite difference (NSFD) scheme to discretize a continuous HBV infection model with capsids. Their work was an extension to the work of Manna and Chakrabarty [17]. Manna and Hattaf [16] studied an HBV infection model which contains two arms of immunity, three time delays, capsids, general incidence rate, and allow the movement of capsids and viruses by diffusion. Manna [18] investigated the role of the CTL immune response in a reaction-diffusion model of HBV with capsids. Moreover, Hattaf and Yousfi [12] developed a mathematical model HBV infection with two modes of transmission, spatial diffusion for both HBV DNA-containing capsids and viruses, and three distributed delays.
All the models mentioned above are formulated either by FDEs or by PDEs. In this study, we will develop a mathematical model governed by both FDEs and PDEs in order to describe the dynamics of HBV infection under the effects of diffusion and memory. To do this, Section 2 deals with model formulation and some preliminary results. Section 3 is devoted to equilibria and their global stability. It is followed by numerical simulations in Section 4 to illustrate the analytical results. The study ends with a conclusion in Section 5.

Model formulation and preliminary results
The memory is an important characteristic of adaptive immunity which means that the immune system can remember the antigens that previously activated it and launch a more intense immune reaction when encountering the same antigen a second time. The classical integer derivative does not reflect this characteristic because it is a local operator unlike the fractional derivative. Indeed, the computation of the fractional derivative at time t = t 1 take into account the entire history from the starting point t = t 0 up to the point t = t 1 . Recently, Baleanu et al. [3] analyzed the interactions between various tumor cell populations and immune system via a system of FDEs. A study in [1] was devoted to the modelling immune systems based on Atangana-Baleanu fractional derivative. Other fractional models with one or two arms adaptive immunity have been proposed in [2,5,29]. However, these above FDE models assumed that cells and viruses are well mixed and ignored their mobility.
Motivated by the above biological and mathematical considerations, we propose the following nonlinear system of FDEs and PDEs: where H(x, t), I(x, t), C(x, t) ,V (x, t) and Z(x, t) are the concentrations of uninfected hepatocytes, infected hepatocytes, HBV DNA-containing capsids, virions and CTL cells at location x and time t, respectively. The uninfected hepatocytes are produced at a constant rate s with a natural death rate µ and become infected at rate f (H, V )V . As in [21], the parameter δ is the death rate for infected hepatocytes and capsids. The parameters a, β and c are, respectively, the production rate of capsids from infected hepatocytes, the rate of export of the capsids to blood, and the clearance rate of virions. Here, the capsids are exported from the cell at the rate β to produce the virions [22]. The capsids are killed by CTL cells at rate p while q and σ denote CTL responsiveness rate and decay rate of CTL cells in absence of antigenic stimulation, respectively. The positive constants d C and d V denote respectively, the diffusion coefficients of capsids and virus. ∆ = n i=1 ∂ 2 ∂x 2 i is the Laplacian operator.
Further, ∂ α is the fractional derivative in the Caputo sense and α is a parameter that describes the order of the fractional time-derivative with α ∈ (0, 1]. Here, the infection transmission is modeled by Hattaf-Yousfi functional response [11] of the form f (H, , where α 0 , α 1 , α 2 , α 3 ≥ 0 are the saturation factors measuring the inhibitory or psychological effect, and k is a positive constant rate describing the infection process. This incidence function generalizes several forms existing in the literature such as the bilinear incidence when α 0 = 1 and α 1 = α 2 = α 3 = 0, the saturation incidence when α 0 = 1 and α 1 = α 3 = 0 or α 2 = α 3 = 0, the Beddington-DeAngelis functional response when α 0 = 1 and α 3 = 0, the Crowley-Martin functional response [7] when α 0 = 1 and α 3 = α 1 α 2 , the specific functional response introduced by Hattaf et al. [13] when α 0 = 1, and the incidence function used in [31] when α 0 = α 3 = 0 and α 1 = α 2 = 1. Additionally, the Hattaf-Yousfi functional response was used by Riad et al. [25] to describe the dynamics of labour market and by Mahrouf et al. [15] to study the dynamical behavior of a viral infection model with immune response and stochastic perturbations. Throughout this paper, we consider system (2.1) with initial conditions and zero-flux boundary conditions where Ω is a bounded domain in R n with smooth boundary ∂Ω, and ∂ ∂ν denotes the outward normal derivative on ∂Ω. From the biological point of view, these conditions mean that the capsids and free virus particles do not move across the boundary ∂Ω.
It is obvious that system (2.1) has one infection-free equilibrium E 0 (H 0 , 0, 0, 0, 0), where H 0 = s µ . Then we define the basic reproduction number of (2.1) as follows The remaining equilibria of system (2.1) satisfy the following algebraic equations The last equation (2.7) implies that either Z = 0 or I = σ q . Each of these cases will lead to one of the other equilibria.
First, consider the case Z = 0. Then by We have g 1 (0) = − δc(β+δ) βa < 0 and g 1 s µ = δc(β+δ) βa (R 0 − 1) and When R 0 > 1, we deduce that system (2.1) admits a unique immune-free infection equilibrium For the case when . Thus, there does not exist any biologically feasible steady state whenever Then we can easily obtain that g 2 (0) = − cq (β + δ) (s − µH) βσa < 0 and In addition to the threshold parameter R 0 , we define the CTL immune response reproduction number R 1 by which describes the average number of CTL immune cells activated by infected hepatocytes in case of successful HBV infection. Here, q denotes the rate of CTL response activation, 1 b represents the average life expectancy for CTL cells and I 1 is the number of infected hepatocytes at the immune-free equilibrium E 1 .
If R 1 > 1, then H 1 < s µ − σδ qµ and g 2 s µ − σδ qµ > 0. Therefore, there exists a unique infection equilibrium with The above discussions gives to the following results.
(i) If R 0 ≤ 1, then system (2.1) has a unique infection-free equilibrium of the form E 0 (H 0 , 0, 0, 0, 0), where H 0 = s µ . (ii) If R 0 > 1, then system (2.1) has an immune-free infection equilibrium of the form E 1 (H 1 , (iii) If R 1 > 1, then system (2.1) has a unique infection equilibrium with CTL immune response of the and By a simple computation, we have the following remark.

Global stability
First, we discuss the global stability of the infection-free equilibrium E 0 and the immune-free infection equilibrium E 1 . (i) The infection-free equilibrium E 0 is globally asymptotically stable when R 0 ≤ 1.
Proof. In order to show (i), we consider the following Lyapunov functional where Φ(x) = x − 1 − ln(x) for x > 0. By using the property of fractional derivatives given in [8], we have Since R 0 ≤ 1, we have that D α W 0 (t) ≤ 0. In addition, the largest invariant set in {(H, I, C, V, Z) | ∂ α t W 0 = 0} is the singleton {E 0 }. It follows from LaSalle's invariance principale [14] that E 0 is globally asymptotically stable when R 0 ≤ 1. This completes the proof.
Next, we construct the Lyapunov functional for system (2.1) at E 1 : The derivative of W 1 (t) along the positive solutions of (2.1) satisfies: By applying s = µH 1 + f (H 1 , V 1 )V 1 , we get and the equality holds only for H = H 1 , I = I 1 , C = C 1 , V = V 1 and Z = 0. Therefore, D α W 1 ≤ 0 if R 1 < 1. Further, the largest compact invariant set in {(H, I, C, V, Z) | ∂ α t W 1 = 0} is the singleton {E 1 }. By the LaSalle's invariance principale, we conclude that E 1 is globally asymptotically stable.
Finally, we investigate the global stability of the third equilibrium E 2 .
Theorem 3.2. The infection equilibrium with CTL immune response E 2 is globally asymptotically stable if R 1 > 1.
Proof. Consider the following Lyapunov functional Calculating the derivative of W 2 (t) along the positive solutions of (2.1), we have By applying the equality From (3.1), we deduce that D α L 2 (t) ≤ 0. Further, the equality D α L 2 (t) = 0 holds only for H = H 2 , I = I 2 , C = C 2 , V = V 2 and Z = Z 2 . Hence, the largest compact invariant set in {(H, I, C, V, Z) | D α L 2 (t) = 0} is the singleton {E 2 }. Then by the LaSalle's invariance principale, we conclude that E 2 is globally asymptotically stable.

Numerical simulations
In this section, we validate our theoretical results by numerical simulations. Let ∆t be the time step size, Ω = [x min , x max ] and ∆x = (x max − x min ) /N be the space step size with N is a positive integer. The grid points for the space are x i = x min + i∆x for i ∈ {0, . . . , N }, and for time are t m = m∆t for m ∈ I N. Based on Grünwald-Letnikov method [24] which recently used in [26], the Caputo fractional derivative can be approximated by For simplicity, we denote the approximations of (H, I, C, V, Z) solution of system (2.1) at the discretized point . By applying (4.1), we find Thus, For numerical illustrations, we take in the whole section Ω = [0, 4] and d C = d V = 0.1. The initial conditions used are: H(x, 0) = 40000, I(x, 0) = 500, C(x, 0) = 300, V (x, 0) = 5 × 10 5 (1 + 0.2 cos( π 3 x)) and Z(x, 0) = 100.
If we take k = 3 × 10 −5 and we keep the values of the other parameters, we find that R 0 = 7.1978 > 1 and R 1 = 0.0328 < 1. Therefore, the solution of system (2.1) converges to the immune-free infection equilibrium E 1 (1.2916 × 10 8 , 3.7275 × 10 3 , 2.0762 × 10 6 , 6.1977 × 10 5 , 0) which is globally asymptotically stable (see Fig. 2).  To confirm numerically the third point of Theorem 3.2, we choose q = 4.4 × 10 −6 and we keep the same values of the other parameters as the second case. By calculation, we have R 1 = 3.2802 > 1. Then model (2.1) has an infection equilibrium with CTL immune response E 2 (1.2923 × 10 8 , 1.1364 × 10 3 , 6.3295 × 10 5 , 1.9894 × 10 5 , 158.7) which is globally asymptotically stable. Figure 3 demonstrates this result.  Here, we observe that the time for arriving at the steady state E 0 with large fractional derivative order value is faster than that with small fractional derivative order value.
On the other hand, Figure 5 shows the impact of the parameter k when it depends on x on the dynamical behavior of the model (2.1). In this case, the solution of (2.1) converges to E 0 even if the value of the basic reproduction number R 0 is sometimes inferior or superior to 1. Then the basic reproduction number R 0 given by the formula (2) is unable to determine the dynamics of the model when the parameters depends on the variable space. We will study this problem in our future work.

Conclusion
In this work, we have studied a time-fractional diffusion model for HBV infection with capsids, CTL immune response and spatial diffusion in capsids and virus. We have shown that the global dynamics of the model is completely determined by two threshold parameters that are the basic reproduction number R 0 and the reproduction numbers for cellular immunity R 1 . We have proved that the infection-free equilibrium E 0 is globally asymptotically stable when R 0 ≤ 1, which means that the virions are cleared and eventually the HBV infection dies out. Whereas, when R 0 > 1 the HBV persists in the liver and two steady states occur, the immune-free infection equilibrium E 1 which is globally asymptotically stable if R 1 ≤ 1, and the infection equilibrium with CTL immune response E 2 which is globally asymptotically stable if R 1 > 1.
From the above analytical and numerical results, we conclude that the diffusion and the order of fractional derivative in sense of Caputo have no effects on the stability of the three equilibria, but it can affect the time for arriving to these equilibria. For instance, when the order of fractional derivative increases the solutions of the model converge rapidly to the equilibrium points (see, Fig. 4).