Stability Analysis and Clinic Phenomenon Simulation of a Fractional-Order HBV Infection Model

In this paper, a fractional-order HBVmodel was set up based on standardmass action incidences and quasisteady assumption.'e basic reproductive number R0 and the cytotoxic T lymphocytes’ immune-response reproductive number R1 were derived. 'ere were three equilibrium points of themodel, and stable analysis of each equilibrium point was given with corresponding hypothesis about R0 or R1. Some numerical simulations were also given based on HBeAg clinical data, and the simulation showed that there existed positive logarithmic correlation between the number of infected cells and HBeAg, which was consistent with the clinical facts. 'e simulation also showed that the clinical individual differences should be reflected by the fractional-order model.


Introduction
Viral infection is a major global problem, and mathematical models are an important tool for the study of biological phenomenon [1][2][3][4][5] and viral infectious disease [6][7][8][9] because they can help us to understand the dynamics of some infectious diseases and some chronic viral infections.
Mathematical models were also used to interpret experimental and clinical results in the fields of (anti-) HIV, HBV, and HCV infections [10][11][12]. ese models were all set up with ordinary differential equation. In recent years, fractional differential equation models were often used in biology because the researchers found that the biological cell membranes have electron conductivity, which could be classified as a fractional-order model [13,14]. In addition, some biological models established by fractional differential equations have proved to be more advantageous than integers [14]. In particular, the biggest difference between the fractional-order model and the integer-order model is that the fractional-order model has the memory, while the characteristic of the immune response contains the memory [14].
So, when we discuss virus immune models, fractional mathematical models have become important tools.
Arafa et al. [14] proposed a fractional-order HIV infection model, Wang et al. [15] proposed a fractional-order HIV infection model, considering the logistic growth of the healthy CD4 cells, and Yan and Kou [16] had further proposed the following HIV model: them is d, a, and c, reprectively, the uninfected hepatocytes are supposed to be produced at rate λ, β is the infection rate at which the uninfected cell becomes infected, and the infected hepatocytes are cured by noncytolytic processes at a constant rate δ per cell. Based on model (2), Cardoso et al. [18] also disused a fractional model of hepatitis B with drug therapy by representing β and k in model (2) with (1 − ρ)β and (1 − μ)k, in which ρ and μ represent the drug efficacy.
It should be pointed out that the bilinear incidences βxv were used to describe the infection between uninfected cells and virus in [14][15][16][17][18], but Min et al. changed it to standard incidence function (βxv/(x + y)) to describe the HBV infection model, which seemed more reasonable because it is independent of the number of total cells of liver [19]. On the other hand, Bartholdy et al. and Wodarz et al. [20] found that the turnover of free virus was much faster than that of infected cells, which means quasisteady assumption could be used, that is, the amount of free virus is simply proportional to the number of infected cells, so the number of infected cells y can also be considered as a measure of virus load. Based on quasisteady assumption, Guo and Cai [21] discussed a HBV infection model by (βxy/(x + y)) instead of (βxv/(x + y)), but the model is integral order. On the other hand, model (2) does not include the immune cell, while the cytotoxic T lymphocyte (CTL) immune response after viral infection is universal and necessary to eliminate or control the disease.
ough fractional differential equations have proved to be a good choice to describe biological phenomenon, most discussions only focus on mathematic analysis and numerical simulation, and so far, almost no papers have used fractional-order models to explain clinic phenomenon about HBV. In this paper, based on the above discussion and quasisteady assumption, a HBV model was set up as follows: in which 0 < α < 1, the meanings of x and y are the same as those in model (2), z represents the number of CTL, where the immune response is assumed to get stronger at a rate cyz and decays exponentially at a rate bz, which is proportional to their current concentration, and the parameter p expresses the efficacy of the nonlytic component. is paper is organized as follows. In Section 2, some definitions and lemmas of fractional-order differential equation are cited. In Section 3, we mainly discuss the existence and uniqueness of positive solutions. e stability analysis is given in Section 4. e simulation is given in Section 5. is paper ends with a conclusion in Section 6 and discussion in Section 7.

Primary Concept and Lemma
It is known that fractional derivative has a variety of definitions [17,22]. In this paper, we used the Caputo fractional derivatives which are defined as follows.
e discriminant D(f) of a polynomial f(x) � x n + a 1 x n− 1 + a 2 x n− 2 + · · · + a n , is the determinant of the corresponding (2n − 1) ⊗ (2n − 1) Sylvester matrix. e Sylvester matrix is formed by filling the matrix beginning with the upper left corner with the coefficients of f(x) and then shifting down one row and one column to the right and filling in the coefficients starting there until they hit the right side. e process is then repeated for the coefficients of f ′ (x). e following lemmas were useful to our arguments.
Lemma 1 (see [23]). For a fractional-order system: with 0 < α < 1 and x ∈ R n , the equilibrium point of the system is locally asymptotically stable if all eigenvalues λ i of Jacobian matrix J � (zf/zx) evaluated at the equilibrium point satisfy Lemma 2 (see [24]). For the polynomial equation, e conditions which make all the roots of (8) satisfy (7) are displayed as follows: (i) For n � 1, the condition is d 1 > 0 (ii) For n � 2, the conditions are either Routh-Hurwitz conditions or (iii) For n � 3, necessary and sufficient conditions, that is, Lemma 3 (see [25]). Assume that the vector function f: R + × R 3 ⟶ R 3 satisfies the following conditions: Here ω and λ are two positive constants; then, initial value Lemma 4 (see [26]

The Existence and Uniqueness of Positive Solutions
For the proof of the existence and uniqueness about the positive solution, we firstly prove that there exists a positively invariant region for system (3). Let and we have Let D � x + y + (p/c)z ≤ (λ/h), x, y, z ≥ 0 ; it is easy to see that D is a positively invariant region for model (3).

Theorem 1. For any initial condition in D, system (3) has a unique solution X(t) � (x(t), y(t), z(t)) T , and the solution will remain nonnegative for all t ≥ 0.
Proof 1. Firstly, we prove the existence and uniqueness of the solution. We denote Obviously f(t, X) satisfies conditions (1)-(3) of Lemma 3, and we only prove that system (3) satisfies the last condition (4) of Lemma 3. Let Hence, Complexity 3 in which m � (λ/h). By Lemma 3, system (3) has a unique solution.
We now prove the solution is nonnegative for all t ≥ 0. From model (3), we know By Lemma 4, the solution is nonnegative. In summary, the system has a unique nonnegative solution.

Stable Analysis
In this section, we will discuss the stability of model (3). is system always has an infection-free equilibrium e basic reproduction number is , which means the uninfected cells and infected cells coexist but the immune response is not activated yet, that is, Further, we will give the immune-response reproductive number , where , in which Δ � λ − dy 2 − βy 2 + δy 2 . Note that x 2 , y 2 , and z 2 must be positive, and we can prove z 2 > 0 holds only when R 1 > 1. Hence, E 2 exists if and only if R 1 > 1. Now, we introduce the main theorem.
e characteristic equation for the infection-free equilibrium E 0 is given as follows: We can see that the characteristic roots all the characteristic roots satisfied |arg λ i | � π > α(π/2), which shows that E 0 is locally asymptotically stable by Lemma 1.
and we have Since Let E be the largest positively invariant subset of M; obviously, E is nonempty since ((λ/d), 0, 0) ∈ E. Let (x(t), y(t), z(t)) be the solution of system (3) with initial value (x 0 , y 0 , z 0 ) ∈ E; then, y(t) � 0. By the first and the third equations of (3), we have and when, [26], when t ⟶ ∞, all solutions in the set D approach E 0 . Noting that E 0 is locally asymptotically stable, the infection-free equilibrium E 0 is globally asymptotically stable.
Finally, we discuss the local stability of the immuneresponse equilibrium state E 2 , and the Jacobian matrix at E 2 is given by e corresponding characteristic equation is where It is easy to verify that Complexity 5 (32) Based on Definition 2, we obtain By Lemma 2, we have the following theorem.

Simulation of the Immune-Response Equilibrium.
We first simulate the stability of E 2 , and the locally stable condition of E 2 is given by eorem 4. e parameters we used are listed in Table 1 which can make R 1 � 578.5076 > 1, D(P) � − 0.1296 < 0 hold. We choose α � 0.55, α � 0.6 and α � 0.65, which satisfy 0 < α < (2/3), and the simulation is shown in Figure 1. e dynamic routes of uninfected and infected cells are shown in Figures 1(a) and 1(b), and the CTL cells are shown in Figure 1(c). e simulation shows that the immune-response equilibrium E 2 is locally asymptotically stable for 0 < α < (2/3), which is consistent with the conclusion of eorem 4.

Simulation of Correlation between HBeAg and Infected Cells.
A positive correlation was found between log-HBV DNA and log-HBeAg in [27], that is, the higher the HBeAg, the higher the HBV DNA load. Since by quasisteady assumption, the amount of free virus is simply proportional to the number of infected cells, the number of infected cell should also be positive logarithmic correlative with HBeAg. Since we could not get the specific number of infected cells from clinic, HBeAg can be quantified. In the following part, we will testify the positive correlation between the infected cell and the HBeAg based on the clinical data of HBeAg by numerical simulation. We use four CHB patients' HBeAg clinical data of 108 weeks with the treatment of entecavir from Dongzhimen Hospital; two patients were from successful-treatment group, and the other two patients were from unsuccessful-treatment group; the data are shown in Tables 2 and 3.
By eorem 2, the infection-free equilibrium E 0 � (x 0 , 0, 0) is asymptotically stable if R 0 < 1. When simulating the successful-treatment group, the parameters we choose should make R 0 < 1. Considering the biological significance of the parameters, we choose the parameters in Table 4 as follows.
A human liver contains approximately 2 × 10 11 hepatocytes [28], and a patient has a total of about 3000 ml plasma. Usually, clinical testing quantification is based on one millilitre. Consequently, we can assume that (λ/d) ≈ (2 × 10 11 /3000). Since the half-life of a hepatocyte is about half a year, we choose d � 0.00379 similar to [28]. e half-life of CTL is about 77 days [29], so we choose b � 0.009. Other parameters were chosen by simulation. e parameters in Table 4 can make R 0 � 0.045 < 1 hold. We give three simulation results with α � 1, α � 0.9, and α � 0.8; the simulations are shown in Figure 2 e simulation shows that the infection-free equilibrium E 0 is asymptotically stable with different 0 < α < 1 which is consistent with the conclusion of eorem 2, but when α � 1, the dynamic route is obviously different. From the simulation in Figure 2(c), we can also observe the positive correlation between log-infected cells and log-HBeAg for different α < 1, which is consistent with clinical fact. On the other hand, from the clinical data of HBeAg, we can see that the change rate was different between different patient, though they have same change trend, which shows that individual differences do exist in clinical therapy even with the same medicine and the same dose; the clinical phenomenon may be explained by our model to some extent. From the simulation, we can see that the larger the α, the slower the decrease rate of infected cells and CTL in the former stage, but the faster the decrease rate in the later stage; the individual difference may be reflected by different α, and the fractional-order model should be a good tool to describe HBV infection.
When we simulate the unsuccessful-treatment group in Table 3, the patients' ALT level was lower than normal standard, so we can think that the patient has not yet developed an immune response to the hepatitis B virus. When R 0 > 1, there exists the equilibrium E 1 without immune activities, so we choose the immune parameter b larger than that in Table 4 and the immune parameter c smaller than that in Table 4. e death rate of infected cells should decrease without immune activities, so we choose parameter a which is smaller than that in Figure 2. Since E 1 is asymptotically stable, if R 1 < 1, the parameter we choose in Table 5 can make R 0 � 2.4253 > 1 and R 1 � 0.1351 < 1 hold. We also used α � 0.9, α � 0.85, and α � 0.8 to simulate, and the simulation is shown in Figure 3. e dynamic routes of uninfected cells and CTL cells are shown in Figures 3(a) and 3(b), and the HBeAg and infected cells are shown in Figure 3(c). From the simulation, the dynamic routes of uninfected cells, infected cells, and CTL cells show that the immune-absence equilibrium E 1 is asymptotically stable with different α which is consistent with the conclusion of eorem 3. From the simulation, there still exists the same change trend but 6 Complexity          Complexity different change rate with different α, and the positive correlation between log-infected cells and log-HBeAgis reflected roughly. e positive correlation is not so good for patient 2's HBeAg data, but the trend is consistent between the infected cells and HBeAg. We can see there also exists some difference with different α which may reflect the individual difference.

Conclusions
In this paper, based on the fact that immune response has memory, we discussed a fractional-order HBV model with standard mass action incidences, and we obtained the basic reproductive numbers R 0 and the cytotoxic T lymphocytes' immune-response reproductive number R 1 . When R 0 < 1, we have proved that E 0 is globally asymptotically stable which meant the infected person can eventually recover automatically even when infected with a large number of HBV. When R 0 > 1 and R 1 < 1, E 1 was locally asymptotically stable which meant the one infected by HBV with no immune response would be infected persistently. When R 1 > 1, we also gave the local stable condition of E 2 , which meant the infected person has immune response to HBV, but the persistent infection still existed. Furthermore, we gave some simulations with different α to test our theoretical results and some clinical phenomena. e simulation showed that our model can simulate the positive correlation between loginfected cells and log-HBeAg to some extent. On the other hand, even with the same initial conditions and the same parameters, we also found that there existed some difference in the dynamic routes with different order α, which may reflect the individual difference. From the analysis and simulation, we can see the fractional-order model maybe more reasonable to describe HBV immune course.

Discussion
In our paper, we found that the fractional-order model can reflect the characteristics of immune memory and clinical individual differences, but it should be pointed out that there are only three variables in our model; the model only described HBV infection roughly. Furthermore, we could not obtain the number of uninfected cells, infected cells, and CTL; when we performed the simulation, we only chose the parameters from biological meaning, not from real clinical data, so there are still many works to be done in the future.

Data Availability
e HBeAg data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.