Hepatitis B virus (HBV) is a partially double-stranded DNA virus and can cause liver diseases such as hepatitis, cirrhosis and liver cancer, and eventually death. When the virus invades the human body, it enters the liver cells by binding to the receptor. Then, HBV-DNA is transported to the nucleus of liver cells and converted into cccDNA under the action of host enzymes. Since cccDNA is not easily degraded, chronic hepatitis B is difficult to be cured [1, 2]. Despite the highly effective recombinant vaccine has been available from 1981 [3], and more and more population are actively vaccinated, the number of the HBV infected is still increasing. More than two billion people have been infected worldwide, and 350 million people are suffering from chronic hepatitis B [4]. HBV has become a major threat to global public health.

Mathematical models can be used to simulate and predict the process of biodynamics. Researchers analyzed the dynamics of the model from different aspects. An immune model was given by Su et al. [5] and the existence, uniqueness and global stability of the model were also analyzed. Anwarud [6] proposed a HBV model with the effect of vaccine and treatment and gave the optimal control strategy. Li et al. [7] analyzed the network vulnerability based on Markov criticality and applied it in lots of biosystem networks. Kalyan [8] considered the intracellular delay in the production of the infected hepatocytes and incorporated a time lag in his model. The spatial diffusion was also used by Khalid [9]. However, most of the papers are limited in the classical integer order.

Fractional calculus is an arbitrary order generalization of classical integral calculus which has the characters of memory, history and nonlocality. There are many definitions of fractional calculus. Among them, Caputo order is one of the widely used types. Kai Diethel [10] summarized the application of Caputo in various field and the numerical methods of fractional order. Devendra [11] derived analysis and numerical solutions of nonlinear fractional problems in nanofluid dynamics, nanoscale heat conduction and carbon nanotube current. Ivan [12] gave a fractional order Ebola epidemic model to study its outbreak. Min Ma [13] considered the diffusion issue in fractional order models and gave the numerical method. Hatıra [14] proposed HIV model with Caputo and constant proportional Caputo operators. Fatmawati [15] gave the Caputo order model of HIV transmission with awareness effect. To reduce the transmission of HIV, Caputo order model for the transmission of HIV epidemic with optimal control was proposed by Parvaiz [16]. The theory of Caputo fractional order in dynamic system has also been derived. The local stability was similar to that in integer order and derived by E. Ahmed [17]. The global stability of fractional order was also given by Boukhouima [18]. And other theorems [19] of the Caputo order systems have been derived and widely used in papers. However, the singular kernel of the Caputo fractional derivative causes a lot of difficult problems both in integral calculation and discretization.

To avoid the shortage of the singular kernel, some researchers have proposed new definitions of fractional order. These definitions include the Caputo-Fabrizio derivative defined by an exponential decay law [20] and Atangana–Baleanu derivative defined by a Mittag–Leffler law [21]. Although Caputo-Fabrizio order does not conform to the traditional fractional calculus [22], it is also used in lots of biodynamic system. Esmehan [23] gave a cancer-immune system and discussed the model by adding IL-12 cytokine and anti-PD-L1 inhibitor. Elvin [24] proposed and analyzed a Caputo–Fabrizio fractional model for the HIV/AIDS epidemic which includes an antiretroviral treatment compartment. Human liver model involving Caputo–Fabrizio derivative was also proposed by Dumitru [25]. Hakimeh [26] proposed a Caputo–Fabrizio model for hearing loss due to Mumps virus with optimal control. The basic theorems of this order have also been proved, including the existence and uniqueness of the solution and the local stability of the equilibrium point [27].

Many papers study HBV with mathematic models. The earliest study can be traced back to the three-state model proposed by Novak et al. [28] and the analysis of the model and basic reproduction number were given. With the development of HBV, researchers have found that the bilinear incidence can not describe the real infection within body well. So Min et al. [29] modified βxv to β x v x + y which can be seen as the standard incidence. Similarly, Y. J. et al. [30] gave a model by using the standard incidence rate and the logistic function. They also gave the optimal therapy in its corresponding control model. The factional order is also used in HBV models. Salman proposed a within-host fractional order HBV model with Caputo order [31]. Compared with models above, this paper considered not only Caputo fractional order, but also the cure of infected hepatocytes. At the same time, Caputo-Fabrizio order can be also used in HBV models. Saeed Ahmad [32] gave a Caputo-Fabrizio order model describing the dynamics of hepatitis B infectious disease and used the Banach fixed point theory approach. Fei Gao [33] proposed the Caputo-Fabrizio with delay, analyzed and simulated the stability of its equilibrium points. However, the infection incidence is bilinear and the real data of HBV infection is not used in his model.

So far, almost no paper considered the Caputo-Fabrizio order model with standard incidence which is more reasonable. And Caputo-Fabrizio order also has the character of memory and can effectively avoid the singularity of Caputo order. Besides, there are seldom papers considered parameter estimation based on the real data and optimal control in Caputo Fabrizio order HBV infection models. Based on the above discussion, we first propose a Caputo-Fabrizio fractional order model with standard incidence. Then, we use the clinical data to estimate the parameters of our model. Finally, we propose an optimal control model based on the treatment for chronic hepatitis B and give the optimal therapy strategy. Using the actual data to estimate the parameters of the HBV model can make our model have medical significance. And the optimal treatment strategy can reduce the cost and side effects while ensuring the therapeutic effect.

In general, the process of hepatitis B infection [34] mainly includes susceptible hepatocytes, infected hepatocytes, and hepatitis B virus. Uninfected hepatocytes are produced at a certain rate. HBV particles can invade hepatocytes through interaction with the receptors on the liver surface. They reproduce new viruses by various functions and release them into the bloodstream with the lysis of infected cells. The hepatitis B virus genes can naturally be lost which lead to a certain proportion of the infected cells transforming into uninfected hepatocytes. At the same time, liver cells and viruses will gradually die out. The simplified process is shown in Figure 1. And based on this process, we first give the Caputo-Fabrizio model based on the standard incidence as follows: D t α 0 C F x = λ α − d α x − β α x v x + y + δ α y D t α 0 C F y = β α x v x + y − a α y − δ α y D t α 0 C F v = k α y − μ α v (1) Here x t , y t , v t are the healthy hepatocytes, infected hepatocytes and hepatitis B virus respectively. α is the fractional order. Uninfected cells produce at a constant rate of λ ^α . Healthy and infected hepatocyte die at the rate of d ^α x and a ^α y respectively. The infected rate is β α x v x + y . Infected hepatocytes can be cured by noncytolytic processes at the rate of δ per cell. Free virus is produced from infected cells at the rate of k ^α y and attenuated at the rate of μ ^α v. The initial values are given as x t 0 = x 0 , y t 0 = y 0 , v t 0 = v 0 .

Assume that λ = λ ^α , d = d ^α , β = β ^α , δ = δ ^α , a = a ^α , k = k ^α , μ = μ ^α , model 1) will be transformed to model (2). And the meaning of parameters in our framework is given in Table 1. D t α 0 C F x = λ − d x − β x v x + y + δ y D t α 0 C F y = β x v x + y − a y − δ y D t α 0 C F v = k y − μ v (2) Since our model is a fractional order model, CF represents the Caputo-Fabrizio order, and here we will give some basic definition of Caputo-Fabrizio fractional calculus in Definition 1 and 2.

Let H 1 a , b = f | f ∈ L 2 a , b and f ′ ∈ L 2 a , b , where L 2 a , b is the space of square integrable functions on the interval a , b .

Definition 1

([21]) The definition of Caputo-Fabrizio order differential was later modified by Losada and Nieto as D t α a C F f ( t ) = 2 1 − α 2 − α M α ∫ a t f ′ ( x ) exp − α t − x 1 − α d x . (3)

Definition 2.

([21]) Let α ∈ 0,1 , then fractional integral of α for function f(x) can be expressed as I t α a C F f t = 2 1 − α 2 − α M α f t + 2 α 2 − α M α ∫ a t f x d x . (4)

Remark 1.

In Definition 1, the definition of Eq. 4 is the mean of f t and ∫ a t f x d x . Therefore 2 1 − α 2 − α M α + 2 α 2 − α M α = 1 , which implies M α = 2 2 − α , α ∈ 0,1 .

 The rest of this paper is organized as follows. Section 3 gives the analysis and simulation of the local stability of equilibrium points. The parameter estimation are given in Section 4, and the last two sections give the theorems and simulation of optimal control model and our conclusion respectively.

In biology models, the stable analysis of the system can help us to understand the process of the biology phenomenon. The following two basic lemmas [35] can be used for stability analysis. Lemma 1 give the characteristic equation form of Caputo-Fabrizio order linear system. Lemma two give the criterion of local stability.

Lemma 1.

Let F t ∈ R n and M ∈ R ^n×n , then the characteristic equation related to the linear system D t α t 0 C F F t = M F t supplemented with the order α ∈ (0, 1) is given by det s I n × n − 1 − α M − α M = 0 .

Lemma 2.

Suppose that the matrix s I n × n − 1 − α M − α M is invertible. Then, if and only if all roots of the characteristic equation have negative real parts, the fractional linear Caputo-Fabrizio system is asymptotically stable at the free equilibrium point.

 There are two equilibria in model (2): The disease-free equilibrium E 0 = λ d , 0,0 and the endemic-free equilibrium E 1 = λ R 0 a − a + d , ( R 0 − 1 ) λ R 0 a − a + d , k μ ( R 0 − 1 ) λ R 0 a − a + d . And the basic reproduction number [36] R ₀ of system 2) is β k μ ( a + δ ) . Here we will give the analysis of the two equilibria.

Theorem 1.

The disease-free equilibrium E ₀ is locally asymptotically stable if R ₀ < 1. If R ₀ > 1, it is unstable.

Proof. The Jacobian matrix at E ₀ is given: J 1 = − d δ − β 0 − a − δ β o k − μ (5) From Lemma 1, the characteristic equation of system (2) is s I − 1 − α J 1 − α J 1 = 0 . We can denote three roots as s 1 = α − 1 + α s 2 = 2 α 1 − α μ a + δ R 0 − 1 − c 1 + c 1 2 + 4 μ a + δ R 0 − 1 2 1 + 1 − α a + δ + μ + 1 − R 0 a + δ 1 − α 2 s 3 = 2 α 1 − α μ a + δ R 0 − 1 − c 1 − c 1 2 + 4 μ a + δ R 0 − 1 2 1 + 1 − α a + δ + μ + 1 − R 0 a + δ 1 − α 2 (6) Here c 1 = a + δ + μ α .

 It is obvious to see that s ₁ < 0 as 0 < α < 1. If R ₀ < 1, s ₂ < 0, s ₃ < 0. Otherwise s ₂ > 0 and s ₃ > 0. From Lemma 2, it can be proved that when R ₀ < 1, the disease free equilibrium E ₀ is locally asymptotically stable. When R ₀ > 1, E ₀ is unstable.

 The local stability of disease-free equilibrium indicates that HBV infection will not persist. Basic reproduction number R ₀ is an important index to measure whether the disease will die out. When R ₀ < 1, the disease will die out, and when R ₀ > 1, we will have the endemic equilibrium E ₁. We give the analysis of the its stability in Theorem 2.

Theorem 2.

If R ₀ > 1, the endemic equilibrium E ₁ is asymptotically stable.

Proof. The Jacobian matrix at E ₁ is given as J 2 = − d − β v 1 y 1 x 1 + y 1 2 δ + β x 1 v 1 x 1 + y 1 2 − β x 1 x 1 + y 1 β v 1 y 1 x 1 + y 1 2 − β x 1 v 1 x 1 + y 1 2 − a − δ β x 1 x 1 + y 1 o k − μ (7) From Lemma 1, the characteristic equation of system (2) is (I − (1 − α)J ₂) −α J ₂ = 0 which can be rewritten as a₁s³ + a₂s² + a₃s¹ + a₄ = 0.

 Let l 1 = β v 1 y 1 x 1 + y 1 2 , l 2 = β x 1 v 1 x 1 + y 1 2 , l 3 = β x 1 x 1 + y 1 , we will get that a 4 = μ α 3 a l 1 + d l 2 a 3 = α 2 a d + d δ + a l 1 + d l 2 + d μ + l 1 μ + l 2 μ + α 2 1 − α 3 a l 1 μ + 3 d l 2 μ a 2 = α a + d + δ + l 1 + l 2 + μ + 1 − α 2 α 3 a l 1 μ + 3 d l 2 μ + 2 1 − α a d + d δ + a l 1 + d l 2 + d μ + l 1 μ + l 2 μ a 1 = 1 + 1 − α a + d + δ + l 1 + l 2 + μ + 1 − α 3 3 a l 1 μ + 3 d l 2 μ + 1 − α 2 d a + μ + δ + a l 1 + d l 2 + l 1 μ + l 2 μ (8) As 0 < α < 1, thus a₁, a₂, a₃, a₄ > 0.

 Set d 1 = a 2 a 1 , d 2 = a 3 a 1 , d 3 = a 4 a 1 , q = ad + dδ + al ₁ + dl ₂ + dμ + l ₁ μ + l ₂ μ, then a 2 a 3 − a 1 a 4 = q a + d + δ + l 1 + l 2 + μ ( a d + d δ + d μ + l 1 μ + l 2 μ ) + 1 − α 2 α 2 2 q 3 a l 1 μ + 4 d l 2 μ + α 3 2 q a l 1 + d l 2 + 1 − α 3 α 3 μ 2 a l 1 + d l 2 8 a l 1 + 11 d l 2 (9) As 0 < α < 1, a ₂ a ₃ − a ₁ a ₄ > 0 and d 1 d 2 − d 3 = a 2 a 3 − a 1 a 4 a 1 2 > 0 . Thus we can prove the stability of E ₁.

 The local asymptotic stability of the endemic equilibrium indicates the possible physical state of the patients without drug intervention. When R ₀ > 1, the disease will always exist.

First we use the predictor-corrector Adams-Bashforth-Moulton (ABM) technique of the Caputo-Fabrizio equation to carry out the simulation. we will give the simulation results of equilibria E ₀ and E ₁. And the initial values are given in Table. 2.

For the disease-free equilibrium E ₀, we select the parameters in Table. 3 which is satisfied the condition in Theorem 1. The simulation results are shown in Figure 2. As we have proved the local stability of the equilibrium in Section 3 and Figure 2 can also validates this proof. Here x can approach to 500 which is λ d , y and v are 0. However, due to the lack of theory, it is difficult to prove the global stability. As we choose the initial values which is quite different from each other, the simulation results always tend to the same value. Therefore, we can think E ₀ might be globally asymptotically stable. Similarly, we use the parameters in Table 4 to the simulate the endemic equilibrium E ₁, the simulation results is shown in Figure 3 which shows the simulation results is consistent with Theorem 2.

Then, we use the real data in Table 5 [37] to estimate the parameters of our model. The estimation range of parameters are 1 × 1 0 5 , 1 × 1 0 7 , 0,100 for λ and k and 0,1 for others which are in line with the actual human physiological structure. Our estimated results are shown in Table 6 which make R ₀ < 1 and simulation is shown in Figure 4. The fractional order is estimated as 0, 7,519, this shows that fractional order is applicable in the model. k is estimated as 0.0574 and μ is 0.799, this shows the drugs have been used and have obvious effect. The estimated results of a and d are 0.01 and 0.0038, which is similar to the data that is used in papers [30]. From Figure 4, we can see that our model can fit the data well. The hepatitis B virus can reduced rapidly in the first few weeks. And then it tends to a more gentle downward trend. However, the infected hepatocytes reduce exponentially, but decreasing rate is slower than that of hepatitis B virus. Due to the reduce of infected hepatocytes and hepatitis B virus, healthy hepatocytes are less likely infected, therefore they are in a state of first increasing and then regional stability. This phenomenon is consistent with Theorem 1.

In paper [38], a combined treatment of HBV infection was designed with Entecavir (ETV) and Chinese herbal Tiaogan Jianpi Jiedu granules (TGJPJD). ETV can be used to reduce the replication of hepatitis B virus [39]. TGJPJD can kill hepatitis B virus by activating immunity. And based on the combined therapy of HBV, we propose our control model. D t α 0 C F x = λ − d x − β x v x + y + δ y D t α 0 C F y = β x v x + y − a y − δ y D t α 0 C F v = 1 − u 1 ( t ) k y − 1 + u 2 ( t ) μ v (10)

The parameters have the same meaning as given in Table. 1 u ₁(t) (0 < u ₁(t) < 1) and u ₂(t) (0 < u ₂(t) < 1) are the control variables and represent the drug effect of ETV and TGJPJD respectively. In order to achieve the therapeutic effect with lower the cost and side effects, we give the following objective function. Here a ₂, a ₃, b ₂, b ₃, c ₁, c ₂ represent the corresponding weight of each variable. J u = 1 2 a 2 y 2 + a 3 v 2 + 1 2 ∫ 0 t f ( b 2 y 2 + b 3 v 2 + c 1 u 1 2 + c 2 u 2 2 ) d t . (11)

The Hamilton function of J u can be written as H = p 1 λ − d x − β x v x + y + δ y + p 2 β x v x + y − a y − δ y + p 3 ( 1 − u 1 ) k y − ( 1 + u 2 ) μ v + b 2 y 2 + b 3 v 2 + c 1 u 1 2 + c 2 u 2 2 (12) Optimization process has many applications in dynamic systems [40–42]. In our paper, the optimal control problem is a process in which we find the best u(t) to minimize the performance index functional J(t). By using Pontryagin’s minimum principlex [43], we can transform the original problem into a constrained functional and find its theoretical optimization. And the optimal control therapy our model based the above theorem which can make J(u) minimum are given as follows: u 1 = min max k y p 3 c 1 , 0 , 1 u 2 = min max μ p 3 v c 2 , 0 , 1 (13) Using the parameters as we estimated in Section 5, we give the simulation results in Figure 5 and Figure 6. From the simulation in Figure 5 and Figure 6, we can see that the dosage of the two drugs should be kept at the maximum at the initial stage, but the dosage of TGJPJD can be gradually reduced at about 138 days and will last 12 days. The dosage of ETV could be reduced at about 150 days, but the reduction of ETV would result in the rebound of virus, it is necessary to further strengthen the dosage of TGJPJD. So the dose of TGJPJD will increase at about 150 days and then reduce. Finally, it will be the lowest dose at about 174 days. With the enhancement of the immunity by TGJPJD, the two drugs can be maintained at a lower level to achieve better therapeutic effect.

In Figure 6, the healthy hepatocytes will always keep increasing and then begin to stabilize. And the infected cells will reduce due to the drugs and the reduce of virus. hepatitis B virus will first reduce at about 150 days. And then due to the decrease of the ETV, the number of viruses showed a rebound trend. However, by adjusting the dosage of TGJPJD granule, the final number of virus can be reduced to less than 500, reaching the standard of cure.

In order to show the advantages of optimal control of therapeutic strategies, we also compare the objective function values of optimal as we give above with the maximizing drug dosage (u ₁ = u ₂ = 0.9) and minimizing drug dosage (u ₁ = u ₂ = 0.1), and show them in Table. 7. The results show that with the optimal therapy, the cost can be minimized, and the corresponding treatment effect can also be satisfying.

In this paper, a Caputo-Fabrizio fractional HBV model with standard incidence were proposed. The local stability of the equilibria were analyzed and global stability were simulated. At the same time, parameters were estimated based on the real data. With the estimated parameters, our model can simulate the data well. Finally, based on clinic combination therapy of HBV with ETV and TGJPJD, we proposed an optimal control model and gave the optimal therapy, which showed therapeutic effect can be achieved with lower drug dosage and side effects. However, cellular immunity and humoral immunity are not carefully considered in our model. In the optimal control model, the metabolism of drugs in the liver is not analyzed, and HBV replication in cells is not reflected. For our future work, we will apply the latest therapeutic regimens for chronic hepatitis B in our model and analyze it according to the corresponding clinical data. The metabolism of drugs should also be considered. Besides, we will consider the effect of cellular immunity and humoral immunity which is essential.