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

Authors: Yongmei Su; Sinuo Liu; Shurui Song; Xiaoke Li; Yongan Ye
Journal: Complexity
 (2020)

---

## Abstract

In this paper, a fractional-order HBV model was set up based on standard mass action incidences and quasisteady assumption. The basic reproductive numberR0 and the cytotoxic T lymphocytes’ immune-response reproductive number R1 were derived. There were three equilibrium points of the model, 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. The simulation also showed that the clinical individual differences should be reflected by the fractional-order model.


---

## Body

## 1. Introduction

Viral infection is a major global problem, and mathematical models are an important tool for the study of biological phenomenon [1–5] and viral infectious disease [6–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–12]. These 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:(1)Dαx=λ−dx+ρx1−x+yxmax−βxv,Dαy=βxv−ay,Dαv=ky−γv.Within-host HBV models of fractional order were also discussed [17, 18], and the model in paper [17] was as follows:(2)Dαx=λ−dx−βxv+δy,Dαy=βxv−a+δy,Dαv=ky−γv,where x,y,andv represent the concentration of uninfected hepatocytes cells, infected cells, and viruses, the death rate of them is d,a,andγ,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 βandk 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–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.Though 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:(3)Dαx=λ−dx−βxyx+y+δy,Dαy=βxyx+y−ay−pyz−δy,Dαz=cyz−bz,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.This paper is organized as follows. In Section2, some definitions and lemmas of fractional-order differential equation are cited. In Section 3, we mainly discuss the existence and uniqueness of positive solutions. The stability analysis is given in Section 4. The simulation is given in Section 5. This paper ends with a conclusion in Section 6 and discussion in Section 7.


## 2. 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.Definition 1.
(see [14]). The Caputo fractional derivatives of order α>0,n−1<α<n,n∈N, are defined as(4)Dαft=1Γn−α∫0tfnst−sα+1−nds,whereΓn−α=∫0∞tn−αe−tdt.Definition 2.
(see [20]). The discriminant Df of a polynomial(5)fx=xn+a1xn−1+a2xn−2+⋯+an,is defined by Df=−1nn−1/2Rf,f′, where f′ is the derivative of f. Rf,f′ is the determinant of the corresponding 2n−1⊗2n−1 Sylvester matrix. The Sylvester matrix is formed by filling the matrix beginning with the upper left corner with the coefficients of fx 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. The process is then repeated for the coefficients of f′x.
The following lemmas were useful to our arguments.Lemma 1 (see [23]).
For a fractional-order system:(6)Dαxt=ft,xt,xt0=x0,with 0<α<1andx∈Rn, the equilibrium point of the system is locally asymptotically stable if all eigenvalues λi of Jacobian matrix J=∂f/∂x evaluated at the equilibrium point satisfy(7)argλi>απ2.Lemma 2 (see [24]).
For the polynomial equation,(8)Pλ=λn+d1λn−1+d2λn−2+⋯+dn=0.The conditions which make all the roots of (8) satisfy (7) are displayed as follows:(i)
Forn = 1, the condition is d1>0(ii)
Forn = 2, the conditions are either Routh–Hurwitz conditions or(9)d1<0,4d2>d12,tan−14d2−d12d1>απ2.(iii)
Forn = 3,(a)
If the discriminantDP of Pλ is positive, then Routh–Hurwitz conditions are the necessary and sufficient conditions, that is, d1>0,d3>0,d1d2>d3(b)
IfD<0,d1≥0,d2≥0,andd3>0, then (7) holds when α<2/3(c)
IfDP<0,d1<0, and d2<0, then (7) holds when α>2/3(d)
IfDP<0,d1>0,d2>0, and d1d2=d3, then (7) holds for all αϵ0,1Lemma 3 (see [25]).
Assume that the vector functionf: R+×R3⟶R3 satisfies the following conditions:(1)
Functionft,xt is Lebesgue measurable with respect to t on R+(2)
Functionft,xt is continuous with respect to xt on R3(3)
∂ft,xt/∂x is continuous with respect to xt on R3(4)
ft,X≤ω+λx,∀t∈R+,X∈R3
Hereω and λ are two positive constants; then, initial value problems DαXt=ft,xt, Xt0=x0, α∈0,1 have a unique solution.Lemma 4 (see [26]).
Letmt∈Cqt0,T,R and assume that ∀t1∈t0,T,mt1=0andmt>0, t∈t0,t1; then, Dqmt1<0.


## 3. 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(10)Nt=x+y+pcz,and we have(11)DαNt=λ−dx−ay−pbcz≤λ−hx+y+pcz,in which h=mind,a,b, so we have(12)Nt≤−λh+N0Eα−htα+λh.LetD=x+y+p/cz≤λ/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 inD, system (3) has a unique solution Xt=xt,yt,ztT, and the solution will remain nonnegative for all t≥0.Proof 1.
Firstly, we prove the existence and uniqueness of the solution. We denote(13)ft,X=λ−dx−βxyx+y+δyβxyx+y−ay−pyz−δycyz−bz.
Obviouslyft,X satisfies conditions (1)–(3) of Lemma 3, and we only prove that system (3) satisfies the last condition (4) of Lemma 3. Let(14)η=λ00,A1=−dδ00−a−δ000−b,A2=−β00β00000,A3=0000−p0000,A4=00000000c,Xt=xtytzt.
Hence,(15)ft,X=A1Xt+ytxt+ytA2Xt+ztA3Xt+ytA4Xt+η,≤A1xt+A2Xt+cpmA3Xt+mA4Xt+η,≤A1+A2+cpmA3+mA4Xt+η,in which m=λ/h. By Lemma 3, system (3) has a unique solution.
We now prove the solution is nonnegative for allt≥0. From model (3), we know(16)Dαxx−0=λ≥0,Dαyy=0=0,Dαyz=0=0.
By Lemma4, the solution is nonnegative. In summary, the system has a unique nonnegative solution.


## 4. Stable Analysis

In this section, we will discuss the stability of model (3). This system always has an infection-free equilibrium E0=x0,0,0, where x0=λ/d.The basic reproduction number is(17)R0=βa+δ.WhenR0>1, system (3) will have an immune-absence equilibrium E1=x1,y1,0, where(18)x1=λa+δa+δd+aβ−a−δ,y1=λβ−a−δa+δd+aβ−a−δ,which means the uninfected cells and infected cells coexist but the immune response is not activated yet, that is, cy1<b. Further, we will give the immune-response reproductive number R1=cy1/b when R1>1, that is, cy1>b, which means immune response is activated. So, whenR1>1, there is another immune-response equilibrium E2=x2,y2,z2, where(19)x2=Δ+Δ2+4dδy22+λy22dc,y2=bc,z2=βcx2pcx2+b−a+δp,in which Δ=λ−dy2−βy2+δy2. Note that x2, y2, and z2 must be positive, and we can prove z2>0 holds only when R1>1. Hence, E2 exists if and only if R1>1.Now, we introduce the main theorem.Theorem 2.
For system (3), if R0<1, the equilibrium E0 is globally asymptotically stable. If R0>1, the equilibrium E0 is unstable.Proof 2.
Jacobian matrixJ at E0 is(20)J=λ+dβ−δ00λ+a+δ−β000λ+b.
The characteristic equation for the infection-free equilibriumE0 is given as follows:(21)λ+dλ+a+δ−βλ+b=0.
We can see that the characteristic rootsλ1=−d<0,λ2=−b<0, and when R0<1, λ3=β−a−δ<0, all the characteristic roots satisfied argλi=π>απ/2, which shows that E0 is locally asymptotically stable by Lemma 1. If R0>1,argλ3=0<απ/2,E0 is unstable.
Let(22)Vx=y+pcz,and we have(23)DαV=Dαy+pcDαz,=βxyx+y−ay−pyz−δy+pccyz−bz,=xx+y∗βy−a+δy−pbcz,≤βy−a+δy=a+δyR0−1.
SinceR0<1, we have DαV≤0. Let M=x,y,z∈D,DαV=0; obviously,M⊂x,y,z∈D,y=0. 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 x0,y0,z0∈E; then, yt=0. By the first and the third equations of (3), we have(24)xt=−λd+x0Eα−dtα+λd,zt=z0Eα−btα,and when, xt⟶λ/d, zt⟶0, by the invariance of E, xt=λ/d, zt=0. By t⟶∞the Lyapunov–LaSalle theorem [26], when t⟶∞, all solutions in the set D approach E0. Noting that E0 is locally asymptotically stable, the infection-free equilibrium E0 is globally asymptotically stable.Theorem 3.
For system (3), when R0>1,(1)
The equilibriumE1 is locally asymptotically stable if R1<1(2)
IfR1>1, the equilibrium E1 is unstableProof 3.
The characteristic equation for the equilibriumE1 is given as(25)detλ+d+βy12x1+y12βx12x1+y12−δ0−βy12x1+y12λ−βx12x1+y12+a+δpy100λ−cy1+b=0,that is,(26)λ−cy1+bλ2+a+δ+d+βy12−x12x1+y12λ+a+δd+βay12−dx12x1+y12=0.
We can seeλ1=cy1−b, and when R1<1, λ1 is negative and argλ1=π>απ/2.
Next, we consider the following equation:(27)λ2+a+δ+d+βy12−x12x1+y12λ+a+δd+βay12−dx12x1+y12=0.
We let(28)B=a+δ+d+βy12−x12x1+y12=a+δ+d+βy1−x1x1+y1=a+δR0−1+d,C=a+δd+βay12−dx12x1+y12=da+σa+δR0−1+aa+δ2R0−12β.
SinceR0>1, we can see B>0,C>0; hence, (27) has two negative real roots, and we denote them by λ2 and λ3 (argλ2=π>απ/2 and argλ3=π>απ/2), so when R1<1, the equilibrium E1 is locally asymptotically stable. When R1>1, λ1 is positive, and argλ1=0<απ/2, E1 is unstable.
Finally, we discuss the local stability of the immune-response equilibrium stateE2, and the Jacobian matrix at E2 is given by(29)J=−d−βy22x2+y22−βx22x2+y22+δ0βy22x2+y22βx22x2+y22−a−δ−pz2−py20−cz2cy2−b.
The corresponding characteristic equation is(30)Pλ=λ3+a1λ2+a2λ+a3=0,where(31)a1=d+a+δ+pz2+βy22−x22x2+y22,=d+βx2x2+y2+βy2−x2x2+y2,=d+βy2x2+y2>0,a2=βy22x2+y22a+pz2+pbz2+da+δ+pz2−βx22dx2+y22,=βy22x2+y22a+pz2+pbz2+dβx2x2+y2−βx22dx2+y22,=βy22x2+y22a+pz2+pbz2+βdx2y2x2+y22>0,a3=pbz2d+βy22x2+y22>0.
It is easy to verify that(32)a1a2>pbz2d+βy2x2+y2>pbz2d+βy22x2+y22=a3.
Based on Definition2, we obtain(33)DP=−1a1a2a3001a1a2a332a1a200032a1a200032a1a2=18a1a2a3+a1a22−4a12−4a22−27a32.
By Lemma2, we have the following theorem.Theorem 4.
For system (3), when R1>1,(1)
If the discriminantDP of Pλ is positive, the equilibrium E2 is locally asymptotically stable for 0<α<1(2)
IfDP<0, then the equilibrium E2 is locally asymptotically stable for 0<α<2/3


## 5. Numerical Simulation

### 5.1. Simulation of the Immune-Response Equilibrium

We first simulate the stability ofE2, and the locally stable condition of E2 is given by Theorem 4. The parameters we used are listed in Table 1 which can make R1=578.5076>1, DP=−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. The 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). The simulation shows that the immune-response equilibrium E2 is locally asymptotically stable for 0<α<2/3, which is consistent with the conclusion of Theorem 4.Table 1
Parameters in Figure1.
ParameterValueParameterValueλ10c0.051β0.01d0.18a0.003p0.006b0.011δ0.000000018Figure 1
Simulation results of immune-response equilibriumE2 with different α.
(a)(b)(c)


### 5.2. 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.Table 2
HBeAg data of successful-treatment group.
01224364860728496108 (week)I819.5134.5512.444.231.120.610.440.380.380.44II874.3112.151.881.651.150.820.360.40.370.31Table 3
HBeAg data of unsuccessful-treatment group.
01224364860728496108 (week)I36.79138.6221.36314.88393.78391.75378.63411.61407.91342.63II44.1185.05140.27141.82125.86142.63158.75194.480.6243.72By Theorem2, the infection-free equilibrium E0=x0,0,0 is asymptotically stable if R0<1. When simulating the successful-treatment group, the parameters we choose should make R0<1. Considering the biological significance of the parameters, we choose the parameters in Table 4 as follows.Table 4
Parameters in Figure2.
ParameterValueParameterValueλ2.53e + 005c0.31e − 007β6.9e − 004a4∗3.79e − 003d3.79e − 003p1.5e − 004b0.9e − 002δ1.9e − 004A human liver contains approximately2×1011 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×1011/3000. Since the half-life of a hepatocyte is about half a year, we choose d=0.00379 similar to [28]. The half-life of CTL is about 77 days [29], so we choose b=0.009. Other parameters were chosen by simulation. The parameters in Table 4 can make R0=0.045<1 hold. We give three simulation results with α=1, α=0.9, and α=0.8; the simulations are shown in Figure 2. The dynamic routs of uninfected cells and CTL cells are shown in Figures 2(a) and 2(b), and the HBeAg and infected cells are shown in Figure 2(c). The simulation shows that the infection-free equilibrium E0 is asymptotically stable with different 0<α<1 which is consistent with the conclusion of Theorem 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.Figure 2
Simulation results of successful-treatment group.
(a)(b)(c)When we simulate the unsuccessful-treatment group in Table3, 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 R0>1, there exists the equilibrium E1 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. The death rate of infected cells should decrease without immune activities, so we choose parameter a which is smaller than that in Figure 2. Since E1 is asymptotically stable, if R1<1, the parameter we choose in Table 5 can make R0=2.4253>1 and R1=0.1351<1 hold. We also used α=0.9, α=0.85, and α=0.8 to simulate, and the simulation is shown in Figure 3. The 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 E1 is asymptotically stable with different α which is consistent with the conclusion of Theorem 3. From the simulation, there still exists the same change trend but different change rate with different α, and the positive correlation between log-infected cells and log-HBeAgis reflected roughly. The 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.Table 5
Parameters in Figure3.
ParameterValueParameterValueλ2.53e + 005c0.31e − 009β(6.9e − 003)∗2a3.79e − 003d3.79e − 003p1.5e − 004b0.9δ1.9e − 004Figure 3
Simulation results of successful-treatment group.
(a)(b)(c)


## 5.1. Simulation of the Immune-Response Equilibrium

We first simulate the stability ofE2, and the locally stable condition of E2 is given by Theorem 4. The parameters we used are listed in Table 1 which can make R1=578.5076>1, DP=−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. The 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). The simulation shows that the immune-response equilibrium E2 is locally asymptotically stable for 0<α<2/3, which is consistent with the conclusion of Theorem 4.Table 1
Parameters in Figure1.
ParameterValueParameterValueλ10c0.051β0.01d0.18a0.003p0.006b0.011δ0.000000018Figure 1
Simulation results of immune-response equilibriumE2 with different α.
(a)(b)(c)


## 5.2. 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.Table 2
HBeAg data of successful-treatment group.
01224364860728496108 (week)I819.5134.5512.444.231.120.610.440.380.380.44II874.3112.151.881.651.150.820.360.40.370.31Table 3
HBeAg data of unsuccessful-treatment group.
01224364860728496108 (week)I36.79138.6221.36314.88393.78391.75378.63411.61407.91342.63II44.1185.05140.27141.82125.86142.63158.75194.480.6243.72By Theorem2, the infection-free equilibrium E0=x0,0,0 is asymptotically stable if R0<1. When simulating the successful-treatment group, the parameters we choose should make R0<1. Considering the biological significance of the parameters, we choose the parameters in Table 4 as follows.Table 4
Parameters in Figure2.
ParameterValueParameterValueλ2.53e + 005c0.31e − 007β6.9e − 004a4∗3.79e − 003d3.79e − 003p1.5e − 004b0.9e − 002δ1.9e − 004A human liver contains approximately2×1011 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×1011/3000. Since the half-life of a hepatocyte is about half a year, we choose d=0.00379 similar to [28]. The half-life of CTL is about 77 days [29], so we choose b=0.009. Other parameters were chosen by simulation. The parameters in Table 4 can make R0=0.045<1 hold. We give three simulation results with α=1, α=0.9, and α=0.8; the simulations are shown in Figure 2. The dynamic routs of uninfected cells and CTL cells are shown in Figures 2(a) and 2(b), and the HBeAg and infected cells are shown in Figure 2(c). The simulation shows that the infection-free equilibrium E0 is asymptotically stable with different 0<α<1 which is consistent with the conclusion of Theorem 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.Figure 2
Simulation results of successful-treatment group.
(a)(b)(c)When we simulate the unsuccessful-treatment group in Table3, 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 R0>1, there exists the equilibrium E1 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. The death rate of infected cells should decrease without immune activities, so we choose parameter a which is smaller than that in Figure 2. Since E1 is asymptotically stable, if R1<1, the parameter we choose in Table 5 can make R0=2.4253>1 and R1=0.1351<1 hold. We also used α=0.9, α=0.85, and α=0.8 to simulate, and the simulation is shown in Figure 3. The 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 E1 is asymptotically stable with different α which is consistent with the conclusion of Theorem 3. From the simulation, there still exists the same change trend but different change rate with different α, and the positive correlation between log-infected cells and log-HBeAgis reflected roughly. The 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.Table 5
Parameters in Figure3.
ParameterValueParameterValueλ2.53e + 005c0.31e − 009β(6.9e − 003)∗2a3.79e − 003d3.79e − 003p1.5e − 004b0.9δ1.9e − 004Figure 3
Simulation results of successful-treatment group.
(a)(b)(c)


## 6. 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 numbersR0 and the cytotoxic T lymphocytes’ immune-response reproductive number R1. When R0<1, we have proved that E0 is globally asymptotically stable which meant the infected person can eventually recover automatically even when infected with a large number of HBV. When R0>1 and R1<1, E1 was locally asymptotically stable which meant the one infected by HBV with no immune response would be infected persistently. When R1>1, we also gave the local stable condition of E2, 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. The simulation showed that our model can simulate the positive correlation between log-infected 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.


## 7. 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.



---

*Source: 8864403-2020-09-19.xml*