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

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 number R0 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.


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][6][7][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]: where xðtÞ, yðtÞ, and vðtÞ 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 R 0 = ðβk/aμÞx, wherex 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: The reproductive number R 0 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: 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]: Perelson and Nelson [13] added logistic function in his 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]: 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][21][22][23][24][25]. Papers [20][21][22] considered the epidemic transmission SEIR model of HBV. In paper [23], a withinhost fractional order HBV model was proposed: In this model, xðtÞ, yðtÞ, and vðtÞ 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: 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, uðtÞð0 ≤ uðtÞ ≤ 1Þ means the control and μð0 ≤ μ ≤ 1Þ is the efficacy of the antiviral drug: For the fractional order model, paper [27] gave the analysis of the control model. Here, xðtÞ, yðtÞ, and vðtÞ have the 2 Journal of Function Spaces same meaning of (1) and u 1 and u 2 represent the drug effect on HBV by interferon (IFN) and lamivudine (LAM): 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: Here, x, y, v, and z have the same meaning as those in model (3). 0 < α < 1:u 1 ð0 ≤ u 1 ≤ 1Þ represents the treatment effect of entecavir (ETV), which can block the replication of virus. u 2 ð0 ≤ u 2 ≤ 1Þ represents the treatment effect of Tiaogan Jianpi Jiedu Granules (TGJPJD), which can accelerate the death of virus. u 3 ð0 ≤ u 3 ≤ 1Þ represents the treatment effect of Tiaogan-Yipi Granule (TGYP) which can enhance immunity.
The structure of this paper is given as follows. In Section 2, 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.
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 Lemma 6. For the discriminant Dð f Þ of polynomial equation is defined by Here, Rð f · f ′ Þ represents the determinant for the corresponding ð2n − 1Þ ⊗ ð2n − 1Þ Sylvester matrix.

The Existence and Uniqueness of Positive Solutions
We first analyse the system (10) without control, that is, Theorem 8. The solution of system (10) is always nonnegative.

Stable Analysis
In this section, we will discuss the stability of the system (10). The system always has an infection-free equilibrium E 0 = ðx 0 , 0, 0, 0Þ, where Here, we have the basic reproduction number as When R 0 > 1, the system (10) will have an immuneabsence equilibrium E 1 = ðx 1 , y 1 , v 1 , 0Þ; here, N 1 = ððq − dÞ/ R 0 − 1 ð ÞÞ − ðμq/kÞ and which means that the infected cells and virus coexist but the immune response is not activated yet, that is, cy 1 < b. Further, we will have the cytotoxic T lymphocyte immune response reproductive number 5 Journal of Function Spaces When R 1 > 1, it means that immune response is activated; the system (10) will have an immune-response equilibrium E 2 = ðx 2 , y 2 , v 2 , z 2 Þ, where Here m, n are given as follows: Theorem 11. For the model of (10): when R 0 < 1, the infectionfree equilibrium point E 0 is locally asymptotically stable; when R 0 > 1, E 0 is unstable.
Proof. The characteristic equation for E 1 is given as follows: Let When R 1 < 1, the characteristic root l 1 = cy 1 − b < 0.

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 7 Journal of Function Spaces control model to find optimal treatment strategy. Here, u 1 , u 2 , and u 3 represent the effects of the drugs ETV, TGJPJD, and TGYP, respectively. Let t f be the endpoint of treatment, we choose the objective function of the optimal control as follows [33]: Here, a 2 , a 3 , b 2 , b 3 , c 1 , c 2 , c 3 represent the corresponding weight of each variable. Let XðtÞ = ðxðtÞ, yðtÞ, vðtÞ, zðtÞÞ T , uðtÞ = ðu 1 ðtÞ, u 2 ðtÞ, u 3 ðtÞÞ T , Xð0Þ = ðx 0 , y 0 , v 0 , z 0 Þ T , then (9) can be rewritten as follows: Then, its Hamilton function of JðuÞ can be expressed as Using Pontryagin's minimum principle, the necessary conditions for the existence of the optimal solution of system (10) are shown as follows: : So, we get the optimal control as

Simulation
In this section, we use the numerical method provided in [34] to simulate our system (10). The parameters are given in Table 1.
According to the clinical trials [4], we choose 108 weeks as a treatment period, that is, the initial time t 0 = 0 and the end time t f = 7 × 108 (days). u 1 , u 2 , u 3 ∈ ½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  [35] q 0.3 [15] x m 1 × 10 8 [36] 8 Journal of Function Spaces 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, and 0: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 monotherapy ði:e:,u 1 ≠ 0, u 2 = u 3 = 0Þ.
From the simulation in Figure 1, 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 Figure 2, 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 u 1 = 0:98 (maximum effect) and u 1 = 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.
Strategy 2. Combination of ETV and CHM ði:e:,u 1 ≠ 0, u 2 ≠ 0, u 3 ≠ 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 u 1 to represent ETV treatment effect, and u 2 represents the treatment effect of TGJPJD, which can accelerate the death of virus. u 3 represents the treatment effect of TGYP which can enhance immunity. The simulations are shown in Figures 3 and 4.
From the simulation in Figure 3, 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 Figure 4, 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 control u 1 = u 2 = u 3 = 0:98 (maximum effect) and u 1 = u 2 = u 2 = 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.

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

10
Journal of Function Spaces 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 11 Journal of Function Spaces comparison showed the optimal therapy can get similar or better treatment effect with less drug dose and side effects.

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors declare that they have no conflicts of interest.