# Analysis and Simulation of a Fractional Order Optimal Control Model for HBV

Authors: Liangli Yang; Yongmei Su; Xiaoke Li; Yongan Ye
Journal: Journal of Function Spaces
 (2020)

---

## Abstract

In this paper, we propose a fractional optimal control model of anti-HBV infection based on saturation incidence and logistic proliferation of uninfected cells for the first time. We derive the basic reproduction numberR0 and the cytotoxic T lymphocyte immune response reproductive number R1 and give the stable analysis based on R0 and R1. We analyse the optimal control condition and give two optimal control strategies about entecavir monotherapy and combination therapy of traditional Chinese medicine and entecavir with different fractional orders by simulation. The simulation shows that combination therapy may be a good choice in anti-HBV infection therapy. We also compare the objective function values of the optimal control strategies with other constant control strategies; the comparison shows that the optimal therapy can get similar or better treatment effect with less drug dose and side effects.


---

## Body

## 1. Introduction

Hepatitis B virus infection is a major global public health problem. Based on the information published by the World Health Organization on World Hepatitis Day 2019, hepatitis B is the second deadly epidemic, and the number of people infected with the hepatitis virus is 9 times that of HIV [1]. The prevalence of chronic hepatitis B has brought a huge medical burden on society. At present, the most widely used drugs for chronic hepatitis B are nucleot(s)ide analogues (NUCs), such as tenofovir and entecavir, which are used to inhibit viral DNA polymerase and reverse transcriptase activity [2]. The combined therapy of NUCs plus Chinese herbal medicine (CHM) is widely accepted in China, which has been recognized as a prospective alternative approach [3, 4].In order to better understand the transmission mechanism of various infectious diseases, many mathematical models were set up to enhance our understanding of the dynamics of infectious diseases and chronic viral infections [5–8]. Mathematic models are also used in the study of anti-HBV infection treatment. The initial model of HBV infection was proposed in 1996 by Zeuzem et al. [9]:
(1)dxdt=λ−dx−βxv,dydt=βxv−ay,dvdt=ky−μv,where xt,yt,andvt represent the concentration of uninfected cells, infected cells, and viruses at time t, respectively.Uninfected cells are assumed to be produced at the constant rate λ, die at the rate of dx, and become infected at the rate of βxv. Infected cells are thus produced at the rate of βxv and die at the rate ay. Free virions are generated from infected cells at the rate of ay and decay at the rate of μv.The infection rate in model (1) is assumed to be linear with respect to the virus. However, the basic reproduction number of model (1) is given by R0=βk/aμx^, where x^ represents the number of total liver cells, which implies that an individual with a larger liver will be more difficult to be cured than a person with a smaller one. Paper [10] changed βxv in (1) to βxv/x+v and gave the following model:
(2)dxdt=λ−dx−βxvx+v,dydt=βxvx+v−ay,dvdt=ky−μv.The reproductive numberR0 of (2) is βk/aμ which seems more reasonable because it is no more dependent on the total number of liver cells. On the other hand, when modelling virus infection, the host’s immune response should not be ignored; Nowak and Bangham [11] proposed the following immune model:
(3)dxdt=λ−dx−βxv,dydt=βxv−ay−pyz,dvdt=ky−μv,dzdt=cyz−bz,where z represents the number of CTL cells, cyz means the activated rate of the immune cells, and pyz indicates the rate at which infected liver cells are eliminated by CTL immune cells. Then, Su et al. changed βxv into βxv/x+v and gave the model as follows [12]:
(4)dxdt=λ−dx−βxvx+v,dydt=βxvx+v−ay−pyz,dvdt=ky−μv,dzdt=cyz−bz.Perelson and Nelson [13] added logistic function in his HIV model, that is dx/dt=λ−dx−px1−x/xmax−βxv. Based on paper [13], Song and Neumann proposed a HIV viral infection model with logistic function and saturated mass action incidence [14]. Then, Yu et al. gave a HBV model with logistic function and standard mass action incidence [15]:
(5)dxdt=λ−dx+px1−xxmax−βxvx+y,dydt=βxvx+y−ay,dvdt=ky−μv.In recent years, more and more fractional order models were used in the biological immune system, because the fractional order model has the memory, while the characteristic of the immune response contains the memory [16, 17]. So, when we set the virus immune model, fractional mathematical models have become an important choice. Many fractional order HIV infection models were set up [18, 19].So far, some fractional order models about HBV have been set up [20–25]. Papers [20–22] considered the epidemic transmission SEIR model of HBV. In paper [23], a within-host fractional order HBV model was proposed:
(6)Dαx=λ−dx−βxv+δy,Dαy=βxv−a+δy,Dαv=ky−μv.In this model,xt,yt, and vt have the same meaning of (1); infected hepatocytes are cured by noncytolytic processes at a constant rate δ per cell. But the model still used the bilinear mass action incidence βxv. Paper [24] also considered the HBV models based on the bilinear mass action incidence βxv. Paper [25] added δy to the first equation in (4) and changed it to fractional order. The model is given as follows:
(7)Dαx=λ−dx−βxvx+v+δy,Dαy=βxvx+v−ay−δy−pyz,Dαv=ky−μv,Dαz=cyz−bz.On the other hand, in the process of anti-HBV treatment, we hope not only to lower the levels of HBV and the infected hepatocytes during and at the end of therapy but also to minimize the therapeutic side effects and the cost of drugs, so it is very important to use optimal control theory to study the anti-HBV treatment model.Paper [26] provided a control model, which uses the linear incidence. Here, ut0≤ut≤1 means the control and μ0≤μ≤1 is the efficacy of the antiviral drug:
(8)dxdt=λ−dx−βxv,dydt=βxv−ay,dvdt=1−μuty−cv.For the fractional order model, paper [27] gave the analysis of the control model. Here, xt, yt, and vt have the same meaning of (1) and u1andu2 represent the drug effect on HBV by interferon (IFN) and lamivudine (LAM):
(9)Dαx=Π−δx−1−u1βxv1+ax+py,Dαy=1−u1βxv1+ax−b+py,Dαv=cy−γv−u2v.Based on the above discussion, according to the clinical anti-HBV combination therapy of traditional Chinese and western medicine [4], we consider a HBV fractional order model as follows:
(10)Dαx=λ−dx−βxvx+v+δy+qx1−xxmax,Dαy=βxvx+v−ay−δy−1+u3pyz,Dαv=1−u1ky−1+u2μv,Dαz=1+u3cyz−bz.Here,x,y,v,andz have the same meaning as those in model (3). 0<α<1.u10≤u1≤1 represents the treatment effect of entecavir (ETV), which can block the replication of virus. u20≤u2≤1 represents the treatment effect of Tiaogan Jianpi Jiedu Granules (TGJPJD), which can accelerate the death of virus. u30≤u3≤1 represents the treatment effect of Tiaogan-Yipi Granule (TGYP) which can enhance immunity.The structure of this paper is given as follows. In Section2, we show some definitions of fractional order and related lemmas. In Section 3, the existence and uniqueness of the positive solution are discussed. In Section 4, we give the stable analysis of our system. The necessary conditions for an optimal control are derived in Section 5. The numerical simulation and the conclusion are given, respectively, in Sections 6 and 7.


## 2. Basic Concepts and Lemmas

First of all, we give some basic definitions and inferences about fractional order. There are several definitions of fractional derivatives, but we only consider the Caputo derivation used in this paper.Definition 1.
The Caputo fractional derivative of orderαα>0 of a function ff:0,∞⟶R is given:
(11)Dαf=1Γn−α∫atfnt−sα−n+1ds.
Here,n=α+1, and Γ· is the Euler gamma function. In this paper, we only discuss the situation that 0<α<1. Function (11) will become
(12)Dαf=1Γ1−α∫atf′t−sαds.Lemma 2 (see [28] (Laplace transform).
The Laplace transform of formula (12) is
(13)LDαft=sαFs−∑k=0n−1fk0sα−k−1.Lemma 3 (see [29] (Generalized mean value theorem).
Suppose thatfx∈Ca,b and Dαfx∈Ca,b, for 0<α<1, we have
(14)fx=fa+1ΓαDαfξx−aα.
Notice. Suppose thatft∈Ca,b and Dαft∈Ca,b, for 0<α<1, we have the following:
(1)
IfDαft>0,t0∈a,b, then neighbourhood N satisfies ft>fa,∀t∈N(2)
IfDαft<0,t0∈a,b, then neighbourhood N satisfies ft<fa,∀t∈NLemma 4 (see [30]).
For a fractional order system (15) of which 0<α<1 and x∈Rn,
(15)Dαxt=ft,xt,xt0=x0.
Assumeft,x:R+×Rn⟶Rn satisfies the following conditions:
(1)
ft,xt is Lebesgue measurable with respect tot on R+(2)
ft,xt and ∂ft,xt/∂x are continuous for all x∈Rn(3)
ft,xt≤ω+λx,∀t∈R+,X∈Rn; here, ω,λ>0 are two positive constants
Then, the initial problem (15) has unique and positive solution on t0,∞.Lemma 5 (see [31]).
For the system (15), the equilibrium point is locally asymptotically stable, if all eigenvalues λi of Jacobian matrix J=∂f/∂x evaluated at the equilibrium point satisfies
(16)argλi>απ2.Lemma 6.
For the discriminantDf of polynomial equation
(17)fx=xn+a1xn−1+a2xn−2+⋯+anis defined by Df=−1nn−1/2Rf,f′. Here, Rf·f′ represents the determinant for the corresponding 2n−1⊗2n−1 Sylvester matrix.Lemma 7.
For the polynomial equation(18)Pλ=λn+d1λn−1+d2λn−2+⋯+dn,the conditions which make all the roots of (18) satisfy (16) are displayed as follows:
(1)
Forn=1, the condition is d>1(2)
Forn=2, the conditions are either Routh-Hurwitz conditions or
(19)d1<0,4d2>d12,tan−14d2−d12d1>απ2Forn=3(a)
If the discriminantDP of (18) is positive, then Routh-Hurwitz conditions are the necessary and sufficient conditions, that is d1>0,d3>0,d1d2>d3(b)
IfDP<0,d1>0,d2>0, and d3>0, then (16) holds when α<2/3(c)
IfDP<0,d1<0,d2<0, and d3≥0, then (16) holds when α>2/3(d)
IfDP<0,d1>0,d2>0, and d1d2=d3, then (16) holds for all α∈0,1(4)
Whenn>3, suppose Δ1,Δ2,⋯,Δn is the Routh-Hurwitz discriminant, that is,
(20)Δ1=a1>0,Δ2=a11a3a2>0,Δ3=a110a3a2a1a5a4a3>0.
Forα∈0,1, (16) is sufficient but not necessary.


## 3. The Existence and Uniqueness of Positive Solutions

We first analyse the system (10) without control, that is, u1=u2=u3=0.Theorem 8.
The solution of system (10) is always nonnegative.Proof.
For system (10), we can get that
(21)Dαxx=0=λ+δy≥0,Dαyy=0=βxvx+v≥0,Dαvv=0=ky≥0,Dαzz=0=0.
Based on Lemma2, we can prove that xt,yt,vt≥0 for t∈t0,∞. For zt, assume that there exists t1 satisfying that
(22)zt1=0,zt>0,∀t∈t0,t1.
We can find thatDαz≥−bz,t∈t0,t1 and zt≥z0Eα−btα,t∈t0,t1, which implies zt1>0. This result is contradicted with the assumption, so zt>0 for any t≥t0.Theorem 9.
There exists anM>0, such that xt,yt,vt,zt<M, for any t∈0,∞.Proof.
Let(23)Nt=xt+yt+pczt,N0=N0=x0+y0+pcz0.
Calculate theα-order derivatives on both sides, respectively,
(24)DαNt=λ−dx−ay−pbcz+qx1−xxmax≤λ+qxmax4−dx−ay−pbcz≤h1−h2x+y+pcz.
Here,h1=λ+qxmax/4 and h2=mind,a,b. So, we can get
(25)Nt≤h1tαEα,α+1−h2tα+N0Eα,1−h2tα=−h1h2+N0Eα−h2tα+h1h2.
Eα−x is completely monotonic for x>0 [32], limt→∞Eα−μtα=0 and Eα0=1. Here, we assume M1=maxh1/h2,N0. So, Nt≤M1, xt≤M1,yt≤M1,andzt≤M1. Similarly, we can have the following inequation:
(26)vt≤−M1kμ+v0Eα−μtα+M1kμ.
Assume thatM2=maxM1k/μ,v0, vt≤M2.
LetM=maxM1,M2, so xt,yt,vt,zt<M.Theorem 10.
System (10) has a unique positive solution on 0,∞.Proof.
LetR4+=x,y,v,z,x,y,v,z≥0. Suppose Xt=xt,yt,vt,zt, system (10) can be transformed into the following form:
(27)ft,X=λ−dx−βxvx+v+δy+qx1−xxmaxβxvx+v−ay−δy−pyzky−μvcyz−bz.
Obviously, (27) satisfies conditions (1) and (2) of Lemma 4. We only prove that system (33) satisfies condition (3) of Lemma 4. Let
(28)η=λ+qx1−xxmax000,A1=−d0000−a0000−μ0000−b,C=−β000β00000000000,A3=00000−p0000000000,A4=000000000000000c,A5=0δ000−δ000k000000,Xt=xtytvtzt.
So, we can get(29)ft,X=A1Xt+vtxt+vtA2Xt+ztA3Xt+ytA4Xt+A5yt+η≤A1Xt+A2Xt+MA3Xt+cpMA4Xt+A5Xt+η≤A1+A2+MA3+cpMA4+A5Xt+h1.
Here,h1=λ+qxmax/4, and M is a bound value mentioned in Theorem 9. By Lemma 4, system (10) has a unique solution on 0,∞.


## 4. Stable Analysis

In this section, we will discuss the stability of the system (10). The system always has an infection-free equilibrium E0=x0,0,0,0, where
(30)x0=xmax2qq−d+q−d2+4qλxmax.Here, we have the basic reproduction number as(31)R0=βkμa+δ.WhenR0>1, the system (10) will have an immune-absence equilibrium E1=x1,y1,v1,0; here, N1=q−d/R0−1−μq/k and
(32)x1=1R0−1v1,y1=μkv1,v1=R0−1xmax2q−N1+N12+4qλR0−12xmax,which means that the infected cells and virus coexist but the immune response is not activated yet, that is, cy1<b. Further, we will have the cytotoxic T lymphocyte immune response reproductive number
(33)R1=μckbR0−1xmax2q−N1+N12+4qλR0−12xmax.WhenR1>1, it means that immune response is activated; the system (10) will have an immune-response equilibrium E2=x2,y2,v2,z2, where
(34)x2=m+n−a23a1,y2=bc,v2=kbμc,z2=λ−dx2−ay2+qx21−x2/xmaxpy2.Herem,n are given as follows:
(35)a1x23+a2x22+a3x2+a4=0,a1=−qxmax,a2=q−d+v2qxmax,a3=q−d−βv2+λ+δy2,a4=v2λ+δy2,uu=9a1a2a3−27a12a4−2a2354a13,vv=34a1a33−a22a32−18a1a2a3a4+27a12a42+4a4a2318a12,m=uu+vv3,n=a22−3a1a39a1m.Theorem 11.
For the model of (10): when R0<1, the infection-free equilibrium point E0 is locally asymptotically stable; when R0>1, E0 is unstable.Proof.
The Jacobi matrix for equilibriumE0=x0,0,0,0 is shown as
(36)JE0=−q−d2+4qλxmaxδ−β00−a+δβ00k−μ0000−b.
So, the characteristic equation forE0 is
(37)lE−JE0=l+q−d2+4qλxmaxl+bl2+a+δ+μl−βk1−1R0=0.
Suppose thatB=a+δ+μ,C=−βk1−1/R0, we can have the four characteristic roots: l1=−q−d2+4qλ/xmax<0, l2=−b<0, l3,4=−B/2±B2−4C/4. When R0<1, −B/2+B2−4C/4<0, from Lemma 5, E0 is locally asymptotically stable. When R0>1, −B/2+B2−4C/4>0, from Lemma 5, E0 is unstable.Theorem 12.
WhenR0>1andR1<1, if DP>0, E1 is locally asymptotically stable for 0<α<1.If DP<0, then E1 is locally asymptotically stable for 0<α<2/3. If R1>1, E1 is unstable.Proof.
The characteristic equation forE1 is given as follows:
(38)lE−JE1=l−cy1+bl3+a1l2+a2l+a3=0.
Let(39)N1=q−dR0−1−μqk,N2=N12+4qλR0−12xmax,a1=a+δ+μ+βv12x1+v12+R0−1−N1+N2>0,a2=μa+δ1−1R0+βv12x1+v12μ+k+a+δ+μR0−1−N1+N2>0,a3=μa+δR0−1−N1+N21−1R0+aβv12μx1+v12>0.
WhenR1<1, the characteristic root l1=cy1−b<0.
LetN3=R0−1−N1+N2, we have
(40)a1a2−a3=βkx12x1+v12+μ2N3+βv12x1+v12+a+δ+βv12x1+v12+N3·μN3+βv12x1+v12+βv12x1+v12+N3a+δ+μa+δ1−1R0>0.
From Lemma6, the discriminant DP of the polynomial equation l3+a1l2+a2l+a3 is as follows:
(41)DP=18a1a2a3−a1a22−4a12−4a22−27a13.According to Lemma7, when R1<1, if DP>0,E1 is locally asymptotically stable for 0<α<1.If DP<0, then E1 is locally asymptotically stable for 0<α<2/3. If R1>1, the equilibrium E1 is unstable.Theorem 13.
WhenR0>1andR1>1, E2 is locally asymptotically stable for 0<α<1.Proof.
The characteristic equation forE2 is given as follows:
(42)lE−JE2=l4+a1l3+a2l2+a3l+a3=0.
Here, letd−q+2qx2/xm=N4≥0,
(43)a1=a+δ+μ+N4+pz2+βv22x2+v22,a2=N4δ+N4+βv22x2+v22μ+a+pz2+μ+bpz2+βk1R0−x22x2+v22,a3=bpz2μ+N4+N4μpz2+βv22x2+v22bpz2+aμ+μpz2+N4βk1R0−x22x2+v22,a4=μbpz2N4+βv22x2+v22.
WhenR0>1,R1>1, we have x22/x2+v22<1/R0, which ensures an≥0,n=1,2,3. So,
(44)Δ2=a2N4+βv22x2+v22+δ+μβk1R0−x22x2+v22+N4δ+μ+N4+pz2N42+βk1R0−x22x2+v22+2N4+μδ+μ+bδ+p2z22N4+μ+b+βv22x2+v22βk1R0−x22x2+v22+μ+pz2δ+μ+pz2+N4δ+2μ+pz2+aN42+βk1R0−x22x2+v22+βv22x2+v22βv22x2+v22+βv22x2+v22+pz22βv22x2+v22+μ+c+2N4βv22x2+v22+δ+μ+pz2+βv22x2+v222μ+pz2≥0,Δ3>a3βv22x2+v22N4+βv22x2+v22+a2βv22μx2+v22N42+2N4βv22x2+v22+δ+μ+bN4+βv22x2+v222+pz2μbβv22x2+v22+N4+βv22x2+v22N4+3βv22x2+v22+βv22μ2x2+v22+aβv22μx2+v22N4+βv12x2+v22+δ+μβv22μx2+v22+N4δ+μ≥0.According to Lemma7, Theorem 13 can be proofed.


## 5. Optimal Control Problem for Combination of Traditional Chinese and Western Medicines

In this section, we analyse the optimal control problem for system (10). Considering the high cost and side effects of long-term use of the highest dose drugs and in order to achieve the treatment effect, we need to analyse the optimal control model to find optimal treatment strategy. Here, u1,u2,andu3 represent the effects of the drugs ETV, TGJPJD, and TGYP, respectively. Lettf be the endpoint of treatment, we choose the objective function of the optimal control as follows [33]:
(45)Ju=12a2y2+a3v2+12∫0tfb2y2+b3v2+c1u12+c2u22+c3u32dt.Here,a2,a3,b2,b3,c1,c2,c3 represent the corresponding weight of each variable. LetXt=xt,yt,vt,ztT,ut=u1t,u2t,u3tT,X0=x0,y0,v0,z0T, then (9) can be rewritten as follows:
(46)DαXt=fXt,ut,t,X0=X0.Then, its Hamilton function ofJu can be expressed as
(47)H=p1λ−dx−βxvx+v+δy+qx1−xxmax+p2βxvx+v−ay−δy−1+u3pyz+p41+u3cyz−bz+p31−u1ky−1+u2μv+b2y2+b3v2+c1u12+c2u22+c3u32.Using Pontryagin’s minimum principle, the necessary conditions for the existence of the optimal solution of system (10) are shown as follows:
(48)Dαx=λ−dx−βxvx+v+δy+qx1−xxmax,Dαy=βxvx+v−ay−δy−1+u3pyz,Dαv=1−u1ky−1+u2μv,Dαz=1+u3cyz−bz,X0=X0,Dαp1=p1q−d−βv2x+v2−2qxxmax+βv2x+v2p2,Dαp2=δp1−a+δ+1+u3pzp2+1−u1kp3+1+u3czp4,Dαp3=−βx2x+v2p1+βx2x+v2p2−1+u2μp3,Dαp4=1+u3pyp2−cy−bp4,−kyp3+c1u1=0,−μp3v+c2u2=0,c3u3−p2pyz+p4cyz=0,p1tf=0,p2tf=a2ytf,p3tf=a3vtf,p4tf=0.So, we get the optimal control as(49)u1=minmaxkyp3c1,0,1,u2=minmaxμp3vc2,0,1,u3=minmaxp2pyz−p4cyzc3,0,1.


## 6. Simulation

In this section, we use the numerical method provided in [34] to simulate our system (10). The parameters are given in Table 1.Table 1
The numerical value of the parameters for system.
ParametersValueReferenceParametersValueReferenceλ6.67×106×d[35]p1×10−4[35]d3.785×10−3[37]k1[35]β1×10−4[10]μ0.35[35]δ3.8×10−3[23]c0.95×10−7[35]a3.38×d[35]b4.5×10−3[35]q0.3[15]xm1×108[36]According to the clinical trials [4], we choose 108 weeks as a treatment period, that is, the initial time t0=0 and the end time tf=7×108 (days). u1,u2,u3∈0,1, by which 0 means no medicine treatment and 1 means the medicine is fully functional. Due to the fact that the human body cannot fully absorb all of the medicine, we assume a maximum drug effect of 0.98. In this simulation, we assume that there is a positive correlation between the efficacy and the dose in a certain range. As α⟶1, the influence of memory decreases [37]. Here, we choose α=0.95,0.9,and0.85 to see the difference between different fractional orders with the lower memory effect.In the following part, we will give optimal control strategies according to different treatment protocols.Strategy 1.
ETV monotherapyi.e.,u1≠0,u2=u3=0.
In this case, the objective function (45) is transformed to
(50)Ju=12a2y2+a3v2+12∫0tfb2y2+b3v2+c1u12dt.The initial condition isx0=6.6×106,y0=2.3×106,v0=2.3×107,z0=200. We give the simulation with different orders α=0.95,0.9,0.85; the corresponding dynamic route simulation is shown in Figure 1, and the optimal control strategy is shown in Figure 2.Figure 1
State dynamics of uninfected and infected cells, free virus, and CTL cells with Strategy1.
(a)(b)(c)(d)Figure 2
Optimal control Strategy1 with ETV monotherapy and α=0.95,0.9,0.85.From the simulation in Figure1, we can see that the change rate and the terminal value of infected cells, virus, and CTL are obviously different with different order α even with the same initial condition and the same change trend, which is consistent with the clinical fact that there do exist individual difference, so individual difference may be reflected by different α, and the fractional order model should be a good tool to describe HBV infection.From the simulation in Figure2, we find that the maximum dose should only be used at the earlier stage of treatment. And the drug dose can be lowered after a period of high-dose treatment. But the time point of lowering the drug dose is different with different order α. If we use constant control u1=0.98 (maximum effect) and u1=0.5 (minimum effect), we give the objective function value of the optimal control and constant control with different orders in Table 2. We can see that the objective function values are lower than those of the maximum and minimum constant control, which shows that the optimal therapy can get similar or better treatment effect with less drug dose and side effects.Table 2
The corresponding number of objective function of Strategy1 with different order and control strategies.
αOptimal controlMaximum effectMinimum effect0.953.27×1083.59×1084.76×1080.92.06×1082.36×1082.87×1080.851.22×1081.55×1081.52×108Strategy 2.
Combination of ETV and CHMi.e.,u1≠0,u2≠0,u3≠0.In this section, we will give the optimal control strategies of combination therapy of ETV and Chinese medicine with different orders asα=0.95,0.9,0.85. We choose the same initial condition and the same parameters as those of Strategy 1 and the objective function (45). As mentioned before, ETV can block the replication of virus; we use u1 to represent ETV treatment effect, and u2 represents the treatment effect of TGJPJD, which can accelerate the death of virus. u3 represents the treatment effect of TGYP which can enhance immunity. The simulations are shown in Figures 3 and 4.Figure 3
State dynamics of uninfected and infected cells, free virus, and CTL cells with Strategy2.
(a)(b)(c)(d)Figure 4
Optimal control strategy with ETV, TGJPJD, and TGYP andα=0.95,0.9,0.85.From the simulation in Figure3, we can see that dynamic routes are similar as that in Figure 1. But the terminal values of infected cells and virus are obviously lower than those in Figure 1, which show that the Chinese medicine can enhance the death of infected cells and virus. The simulation can also show that the individual difference may be reflected by different α.From the simulation in Figure4, we also find that the maximum dose should only be used at the earlier stage of treatment. The drug dose can be lowered after a period of high-dose treatment with different order α. And the time point of lowering the drug dose is also different with different order α. We also find that the time point of lowering the drug dose is different with different medicine; that is, ETV can be lowered at the earliest time, and then, TGJPJD and TGYP need to be taken at the maximum dose for the longest time which shows the importance of enhancing the immune in anti-HBV therapy. Moreover, for Strategy 1, the time point of lowering the ETV drug dose is about 100 days, but for strategy 2, the time point is only about 60 days.We also compared the objective function values of constant controlu1=u2=u3=0.98 (maximum effect) and u1=u2=u2=0.5 (minimum effect) with the optimal control. The results are shown in Table 3. We can see that the optimal objective function values are always the lowest for different orders. By comparing the optimal objective function values in Table 2 with those in Table 3, we find that combination of ETV and CHM can greatly reduce objective function values, which show that combination of ETV and CHM may be a good choice in anti-HBV infection therapy.Table 3
The corresponding number of objective function of Strategy2 with different order and control strategies.
αOptimal controlMaximum effectMinimum effect0.951.81×1082.12×1085.78×1080.91.28×1081.59×1083.44×1080.851.01×1081.34×1081.75×108


## 7. Conclusion and Discussion

In this paper, based on the combination therapy of traditional Chinese medicine and Western medicine, we proposed a fractional optimal control model of anti-HBV infection based on saturation incidence and logistic proliferation of uninfected cells for the first time. After giving the stable analysis of our system, we discussed the necessary conditions of the optimal control problem. Two optimal control strategies about ETV monotherapy and combination therapy of ETV and CHM with different fractional orders were given by simulation.In the simulation, we suppose that there is a positive correlation between drug use and effectiveness. From the simulations, we know that the dynamic routes were different with different orders even with the same parameters, which may show that the individual difference could be reflected by differentα.By comparing the dynamic routes between ETV monotherapy and combination therapy of ETV and CHM, we found that combination therapy can not only obtain better treatment effect but also reduce the taking time and dose of ETV. The simulation shows that combination therapy may be a good choice in anti-HBV infection therapy. We also compared the objective function values of the optimal control strategies with other constant control strategies; the comparison showed the optimal therapy can get similar or better treatment effect with less drug dose and side effects.



---

*Source: 8175217-2020-06-23.xml*