Modeling the Adaptive Immune Response in an HBV Infection Model with Virus to Cell Transmission in Both Liver with CTL Immune Response and the Extrahepatic Tissue

: The objective of this paper is to investigate a mathematical model describing the infection of hepatitis B virus (HBV) in intrahepatic and extrahepatic tissues. Additionally, the model includes the effect of the cytotoxic T cell (CTL) immunity, which is described by a linear activation rate by infected cells. The positivity and boundedness of solutions for non-negative initial data are proven, which is consistent with the biological studies. The local stability of the equilibrium is established. In addition to this, the global stability of the disease-free equilibrium and the endemic equilibrium is fulﬁlled by using appropriate Lyapanov functions. Finally, numerical simulations are performed to support our theoretical ﬁndings. It has been revealed that the fractional-order derivatives have no inﬂuence on the stability but only on the speed of convergence toward the equilibria.


Introduction
The hepatitis B virus (HBV) attacks liver cells (hepatocytes) and kills almost a million people every year. The HBV is known as a global public health issue [1,2], and it infects around two billion people, with more than 360 million chronic carriers. This serious infection can be easily transferred by any infected bodily fluid contact [3][4][5], and it can cause acute or chronic disease following transmission. Many mathematical models have been constructed to represent and better comprehend the dynamics of this deadly viral illness [6][7][8][9]. All of these models take into account how the HBV virus interacts with both healthy and diseased liver cells. However, HBV infection models with two compartments are rare. Recently, the model for describing the adaptive immune response in an HBV infection model with virus to cell transmission in both liver with CTL immune response and the extrahepatic tissue is studied in [10]. They have used the following system of differential equation: dV(t) dt = k 1 I 1 + k 2 I 2 − σV, where we denote the liver and the second extrahepatic compartment as compartments C 1 and C 2 , respectively. In this model, H 1 , I 1 , H 2 , I 2 , V and W denote the concentrations of uninfected hepatocytes in C 1 , infected hepatocytes in C 1 , uninfected hepatocytes in C 2 , infected hepatocytes in C 2 , free virus and the CTL immune response (in C 1 ) at time t, respectively. In addition, s i is the recruitment rate of healthy cells and 1 d i is the average lifespan of uninfected cells in compartment C i (i = 1, 2). The healthy cells become infected by free virus at a rate b 1 H 1 V, infected cells in compartment C i (i = 1, 2) die at a rate a i I i , and infected cells in compartment C 1 are cleared by the CTL immune response at a rate pI 1 W. Free virus (V) grow in blood at a rate k 1 I 1 + k 2 I 2 , decay at a rate σV. Finally, CTLs (W) expand in response to viral antigen derived from infected cells at a rate qI 1 W and decay in the absence of antigenic stimulation at a rate rW.
In this paper, we will study the same above problem but using the fractional derivative equations that will be considered is as follows: where α is the order of the fractional derivative and H 1 (0) = H 10 , are the initial data, such as our dynamical model with two HBV proliferative compartments: One is the liver, whereas the other is the extrahepatic compartment, which includes serum, peripheral blood mononuclear cells and perihepatic lymph nodes. According to the investigations, the immune response has an influence, or it has no impact on the extrahepatic compartment. It is considered that the CTL response plays a half function in clearing infected hepatocytes. Furthemore, it has no role in the second proliferative compartment of extrahepatic tissue. We notice that the fractional derivative order can be considered as an index of memory in many biological and physical problems [11]. For instance, for evolution problems in biology, we can find this kind of derivative to model HBV [12]. The fractional derivative equations were used study other models such as Kawahar and KdV equations [13][14][15][16][17]. Different papers studied the behavior of many viral infections by using the fractional derivative equations [18][19][20][21][22]. We notice that the Laplace operator has been used to study the stability of a fractional-order delayed predator-prey system [23]. In addition, fractional derivative have shown its efficiency in studying many biological systems [24][25][26]. Recently, optimal control problems have been studied using fractional derivative equations [27]. These same techniques were applied to study the interaction between tumor cells and the immune system [28]. Therefore, our motivation in this paper is to take into account the memory effect in our biological dynamical system by incorporating the fractional derivative instead of the ordinary one.
The paper is organized as follows. The next section is devoted to the existence, positivity and boundedness of solutions, which is followed in Section 4 by the local stability analysis. In Section 5, we study the global stability of the equilibrium. We construct an appropriate numerical algorithm and give some numerical simulations in Section 6, and the last section concludes the work.

Preliminary Results
Is this section, we recall some preliminary definitions of the fractional order integral, Caputo fractional derivative and Mittag-Leffler function [29]. Definition 1. The fractional integral of order α > 0 of a function ψ : R + → R is defined by where Γ (.) stands for Gamma function.
Definition 2. The Caputo fractional derivative of order α > 0 of a function ψ : R + → R is given as follows D α ψ(t) = I n−a D n ψ(t) where D = d d t and n − 1 ≤ α ≤ n, n ∈ N. In addition, if 0 < α ≤ 1, we have is called the Mittag-Leffler function of parameter α.
Let f : R n → R n where n ≥ 1. Consider the fractional order system with 0 < α ≤ 1, and X 0 ∈ R n . For the global existence of solution of system (2), we have the following lemma.

Lemma 1. Assume that f satisfies the following conditions
• f (X) and (∂ f /∂X)(X) are continuous on R n . • f (X) ≤ c 1 + c 2 X for all X ∈ R n , with c 1 and c 2 are two positive constants. Then, the system (2) has a unique solution defind on R + .
The proof of this lemma follows immediately from [20].

Positivity and Boundedness of Solutions
It is commonly known that any solutions reflecting cell densities in issues involving cell population evolution should be non-negative and bounded. As a result, demonstrating the model's positivity and boundedness will be important. First of all, for biological reasons, the initial data H 10 , I 10 , H 20 , I 20 , V 0 and W 0 must be larger than or equal to 0. The main result of this section is given as follows: For any non-negative initial conditions, there exists an unique solution of system (1) defined on R + . In addition, this solution is non-negative and bounded t ≥ 0. Proof. The model (1) can be described as follows: Note that if α = 1, (3) will be a system with ODEs. By using the results in [19], we establish the existence of solutions. In the case of FDEs, let This implies that Hence, the proprieties of the previous Lemma are satisfied. Then, system (1) has a unique solution on R + . Now, we will show that the non-negative solutions R 6 + , we have: Therefore, all solutions initiating in R 6 + are positive. Next, we will prove that these solutions remain bounded.
Therefore, W is bounded. From the fifth equation of (1), we obtain t 0 e σ(ξ−t) (I 1 (ξ) + I 2 (ξ))dξ We conclude that all the solutions are bounded and also, each local solution can be prolonged up to any positive time, which means that the solution of the problem exists globally.

Local Stability of Equilibrium
At any equilibrium system, we have by the first and third equation of this system, we obtain Substituting them into the second and fourth equations of the same system, they yield, respectively, When V = 0 and W = 0, substituting I 1 and I 2 of (3) into the fifth equation of (4) We can see that this function f (V) is decreasing with respect to V. We have Thus, (1) has a boundary equilibrium E 1 H 2 , V (1) , 0 when f (0) > 1. When W = 0, from the last equation of system (4), we have I 1 = r q := I  (4) gives Then, a necessary condition on the existence of the positive equilibrium is s 1 > a 1 r q , and for the positive equilibrium E 2 H In the case s 1 > a 1 r q , we have Then, according to the monotonicity of function g(V), we have lim V→0 + g(V) = +∞; then, g(V) = 1 has a unique root in the interval V , 1 σ k 1 s 1 For all this, we will determine the steady states of the model (1). The basic infection reproductive number of the system (1) is given by From a biological point of view, R 0 denotes the average number of secondary infections generated by one infected cell when all cells are uninfected. Moreover, the same system has the following disease-free equilibrium In addition, model (1) admits two endemic steady states , where V (1) is the positive root of (4). This endemic equilibrium exists when R 0 ≥ 1, and Here r and V (2) is the positive root of (4). This endemic equilibriam exists when R 1 ≥ 1.

Global Stability of the Equilibrium
In this section, we will study the global stability of each equilibrium of system (1) by using some suitable Lyapunov functional and by using LaSalle's invariant principle [4]. For the infection-free equilibrium E f , we have the following result: , 0, 0, 0 is globally asymptotically stable.

Proof. When
We can choose a positive number m 1 satisfying the inequality Furthermore, for the given m 1 , we choose a positive number m 2 satisfying the inequality 1 g When m 1 and m 2 are given, we have a Lyapunov function the derivative of L f along solutions of our model is given by  [8] that the infectionfree equilibrium E f is globally asymptotically stable whenever R f ≤ 1.
Proof. First, we define a Lyapunov function L 0 as follows: calculating the derivative of L f along positive solutions of system (4), it follows that Thus, it follows from LaSalle's invariance principle [8] that the infection-free equilibrium E 1 is globally asymptotically stable whenever R 0 ≤ 1.
For the second endemic equilibrium E 2 , we have the following result: 2 , V (2) , W (2) of system (1) is globally asymptotically stable.

Numerical Simulations
In the previous sections, we have studied the theoretical part of our fractional problem (1). In this present section, we will present some numerical simulations to the same problem. We notice that the simulation parameters were inspired from [10]. Figure 1 shows the infection dynamics for the following parameters as in the first numerical result of [10]: s 1 = 1, s 2 = 0.8, b 1 = 0.08, b 2 = 0.1, d 1 = 0.2, d 2 = 0.2, a 1 = 0.2, a 2 = 0.2, k 1 = 1, k 2 = 1.2, p = 0.3, r = 0.2, σ = 5 and q = 0.1. Within these chosen parameters, we have that the basic reproduction number is less than unity R 0 = 0.88 < 1, and we observe the convergence of the curves, which corresponds to the stability of the disease-free equilibrium E f = (5, 0, 4, 0, 0, 0). This result is completely is in good agreement with [10] and also confirm our theoretical result given in Theorem 1. Figure 2 shows the infection dynamics for the following parameters: s 1 = 1, s 2 = 0.8, b 1 = 0.08, b 2 = 0.1, d 1 = 0.4, d 2 = 0.2, a 1 = 0.2, a 2 = 0.2, k 1 = 1, k 2 = 1.2, p = 0.3, r = 0.5, σ = 2.3 and q = 0.56. Within these parameters, we can easily compute the reproduction numbers R 0 = 1.4783 > 1 and R 1 = 0.8929 < 1, which means that the first one is greater than unity while the second is less than one. This predicts the numerical stability of the first endemic equilibrium E 1 . Indeed, we observe that the curves converge toward the first endemic equilibrium E 1 = (2.009, 0.9807, 2.4844, 1.5155, 1.2224, 0), which confirms our theoretical finding concerning the stability of E 1 .

Conclusions
In this paper, we have studied a mathematical model describing HBV infection in intrahepatic and extrahepatic tissues. In our suggested model, we have taken into consideration the effect of CTL immune response. The positivity and boundedness of solutions for non-negative initial data were proven, which is consistent with the biological background. The local stability of the equilibrium was established. Moreover, the global stability of the disease free equilibrium and the endemic equilibrium was also fulfilled by using some appropriate Lyapanov functional. Numerical tests were performed to support our theoretical findings. In the end of this paper, we have studied the effect of the fractional derivative on the numerical stability of each equilibrium. It has been revealed the fractional order derivative has no influence on the stability but only on the speed of convergence toward the equilibria.