Hepatitis B Virus Dynamics : Modeling , Analysis , and Optimal Treatment Scheduling

Modeling, analysis, and control of hepatitis B virus (HBV) infection have attracted the interests ofmathematicians during the recent years. Several mathematical models exist and adequately explain the HBV dynamics as well as the effect of antiviral drug therapies. However, none of thesemodels can completely exhibit all that is observed clinically and account the full course of infection. Besides model inaccuracies thatHBVdynamicsmodels suffer from, some disturbances/uncertainties fromdifferent sourcesmay arise in the modeling. In this paper, the HBV dynamics is described by a system of nonlinear ordinary differential equations. The disturbances or uncertainties are modeled in the HBV dynamics model as additive bounded disturbances. The model is incorporated with two types of drug therapies which are used to inhibit viral production and prevent new infections. The model can be considered as nonlinear control system with control input is defined to be dependent on the drug dose and drug efficiency. We developed treatment schedules for HBV infected patients by using multirate model predictive control (MPC). The MPC is applied to the stabilization of the uninfected steady state of the HBV dynamics model. The inherent robustness properties of the MPC against additive disturbances are also shown.

Several drug therapies have been proposed for treating persons with chronic HBV including adefovir dipivoxil, alpha-interferon, lamivudine, pegylated interferon, entecavir, telbivudine, and tenofovir [33].Hepatitis antiviral drugs prevent replication of HBVs and save the liver from cirrhosis and cancer.During the treatment, the viral load is reduced and consequently the viral replication in liver is decreased [34].
Optimal control theory can be used to optimize the drug doses required in the treatment of HBV infected patients [33,35,36].In this paper, the optimal treatment schedules will be designed using model predictive control (MPC) method.MPC obtains a feedback control by solving a finite horizon optimal control problem at each time instant using the current state of the system as the initial state for the optimization and applying "the first part" of the optimal control [37].MPC has many advantages which include its facility of handling constraints, its closed-loop stability, and inherent robustness.Moreover, MPC solves optimal control problem on-line for the current state of the plant which is a mathematical programming problem and is much simpler than determining the feedback solution by dynamic programming [37].
Recently, MPC has been applied for determining the optimal treatment schedules for HIV infected patients in [38][39][40][41].As far as the authors know, the application of MPC to the HBV dynamics model has never been studied before.The aim of the present paper is to apply a multirate MPC approach to an HBV dynamics model.The MPC is constructed based on the approximate discrete-time model.The closed-loop stability of the multirate MPC will be shown.The inherent robustness properties of the MPC against disturbances will also be shown.The disturbances may arise from modeling errors, or disturbances arise from immune system fluctuating or immune effect of a coinfection.

HBV Infection Model
We will use the mathematical model proposed in [22], incorporating to allow some additive disturbances and taking into account the effect of two types of antiviral drugs Here , , and  denote the concentration of uninfected hepatocytes, infected hepatocytes, and free virions, respectively.Uninfected hepatocytes are assumed to be produced at the constant rate  and die at the rate of .Uninfected hepatocytes can also be created by proliferation of existing  cells.Here we represent the proliferation by a logistic function in which  is the maximum proliferation rate of target cells, and  max is the  concentration at which proliferation shuts off.The rate of infection is given by saturation functional response  − 1 /(1 + ), where  is the infection rate constant characteristic of the infection efficiency and  > 0 is constant.The death rate of infected hepatocytes is .Free virions are assumed to be produced from infected hepatocytes at the rate of  − 2  and  is the clearance rate of viral particles.Functions   ,  = 1, 2, 3, describe the model uncertainties/disturbances.Two types of drugs are considered; the first is used to prevent the virus from infecting the cell and is represented by the chemotherapy function  − 1 , and the second is used to prevent the infected cells from producing new viruses and is represented by the function  − 2 .Here,   () =     (),  = 1, 2 are the control inputs, where the parameters  1 and  2 are the efficiencies of the drugs and  1 () and  2 () are the drug doses (see [42]).
All the parameters of the model are supposed to be positive.System (1)-( 3) can be considered as a nonlinear control system with (, , )  as the state vector and ( 1 ,  2 )  is the control input vector.We are now ready to present a study on the basic mathematical properties of system (1)- (3).
We assume that the model uncertainty/disturbances are given by where  is the external disturbances and Ψ  satisfies the following bound     Ψ  (, , ,  1 ,  2 , )     ≤   ,   ≥ 0,  = 1, 2, 3. ( Now we show that under which conditions the nonnegative orthant R 3 + is positively invariant for (1)-( 3): This means that the nonnegative orthant R 3 + is positively invariant, namely, if a trajectory starts in the nonnegative orthant, it remains there.Proposition 1. Assume that the disturbances satisfy the boundedness condition (5) and  max > + 1 .Then there exist positive numbers  1 and  2 for which the compact set is positively invariant.
Proof.From (1) we have This means that, even if  reaches  max , it should decrease.
Note that Ω contains all the biologically relevant states, thus we can restrict the state space of the system to the compact set Ω. Since the drug doses cannot be arbitrarily increased, we then may consider only a compact control constraint set.

Multirate MPC for Nonlinear Systems
In this section, we outline the multirate MPC design for nonlinear systems in the presence of bounded disturbances and give a review on the results obtained in [43,44].
The set of real and natural numbers (including zero) are denoted, respectively, by R and N. The symbol R ≥0 denotes the set of nonnegative real numbers.A continuous function for every  ≥ 0, it is strictly decreasing in  for every  > 0 and (, ) → 0 when  → ∞.Let B Δ = { ∈ R  : ‖‖ ≤ Δ} will be used in R  .
Consider a continuous-time nonlinear control system with additive disturbances given by where () ∈ R  , () ∈  ⊂ R  , () ∈  ⊂ R  are the state, control input, and disturbances, respectively,  : R  ×  → R  is continuous and Lipschitz continuous with respect to  in any compact set and (0, 0) = 0,  is compact and 0 ∈ ,  is compact and 0 ∈ .
The control is taken to be a piecewise constant signal where  > 0 is the control sampling period which is fixed.
In this paper we address the problem of state feedback stabilization of (15) under "low measurement rate." More precisely, we will assume that state measurements can be performed at the time instants   ,  = 0, 1, . ..: where   is the measurement sampling period.We assume that   = ℓ for a fixed integer ℓ > 0.
Assumption 11.Let X 1 = X + B 1 , (i)  : X 1 → R is continuous, positive definite radially unbounded, and Lipschitz continuous in any compact set, (ii)  ℎ (, ) is continuous with respect to  and , uniformly continuous in small ℎ, and Lipschitz continuous in any compact set, (iii) there exist an ℎ * > 0 and two class-K ∞ functions  1 and  2 such that the inequality holds for all  ∈ X 1 ,  ∈  and ℎ ∈ (0, ℎ * ]. Assumption 12.There exist ℎ * > 0 and  > 0 such that for all  ∈ G  = { : () ≤ } there exists a () ∈  such that inequality holds for all ℎ ∈ (0, ℎ * ]. We define the value function, which represents the optimal value of ( 26) for a given initial condition as follows Let u * = { * 0 , . . .,  * −1 } be the solution of this optimization problem, then the first ℓ elements of u * are applied at the state , that is, Let ℎ * 0 denote the minimum of the values ℎ * generated by Assumptions 9-12.Let Δ  and Δ  be such numbers that, for each  ∈ X,  ∈ , ‖‖ ≤ Δ  , ‖‖ ≤ Δ  .

MPC for the HBV Model
In this section we apply the MPC method proposed in Section 3. We will show that, with a suitable choice of  and functions  and  ℎ , the assumptions of the previous section can be satisfied.Introduce new variables as  1 =  −  0 ,  2 = ,  3 = .In these new variables the model ( 1)-( 3) takes the form of (15) with and  = ( 1 ,  2 ,  3 )  .Let the compact set X be defined as where  1 and  2 are given in Proposition 1.With this definition,  satisfies all regularity assumptions, and according to Proposition 1, X is positively invariant.
To verify Assumptions 11 and 12, we linearized the nominal system (35) around the origin in case of constant controllers, that is,  1 () =  1 >  (1)   ,  2 () =  2 >  (2)    with  (1)   +  (2)    =   , where   is given in Proposition 3. Let   be the coefficient matrix of the linearized system and  = ( −  0 , , )  .Then the discrete-time model for the linearized system is given by Let the sampling period be chosen to be  = 0.5 and  1 =  2 = 1.The running cost and the terminal cost can be chosen as where  is a positive definite diagonal matrix and  is a positive definite symmetric matrix satisfying the Lyapunov equation for the discrete-time system (37) Using the parameters given in Table 1, we have verified Assumption 5 numerically by fixing one controller, for example,  1 and solve the inequality (13) with respect to  2 .We have found that for  1 = 1, then inequality ( 13) is satisfied when  2 > 0.20397.
From ( 38)-( 39), Assumption 11 is satisfied.Assumption 10 follows from Proposition 6. Assumption 9 is fulfilled as well, if we choose a suitable numerical integration scheme (e.g., the Runge-Kutta formula).To verify Assumption 12, the matrix  has been chosen through a series of numerical experiments  = diag (1,20,3).It has been verified numerically by solving a constrained minimization problem with several starting points, for which Assumption 12 is satisfied over the whole set X. Thus all assumptions of the proposed method can be satisfied with suitable choice of the parameters of the MPC method.

Numerical Results
We performed simulation studies using the following parameter values which are listed in Table 1.
We assumed that the state measurements were performed at the instants ℓ,  = 0, 1, . . .with ℓ = 4 and the horizon length  was chosen to be  = 8.All computations were carried out by MATLAB.Especially, the optimal control sequence was computed by the fmincon code of the Optimization toolbox.Simulations for the continuous-time system were carried out using ODE45 program in MATLAB for two cases, untreated case and treated case.
(i) Untreated Case.In this case, we simulate the nominal system (1)-(3) (i.e.,   () = 0,  = 1, 2, 3) without treatment (i.e.,  1 =  2 = 0).Using the parameters given in Table 1 and choose  = 0.005, we have obtained (/) 0 = 3.33333 > 1 and  0 (0.0) = 1.60987 > 1, moreover, all the conditions of Theorem 8 are satisfied.Therefore the uninfected steady state  0 is unstable and the infected steady state  1 exists and is globally asymptotically stable.To show the simulation results for the untreated case, we assume that the infection occurs with a certain amount of virus particles  = 0.001.Thus the initial conditions for the untreated case are (0) =  0 = 1000, (0) = 0, and (0) = 0.001.From Figures 1, 2, and 3, we can see that, the concentration of uninfected hepatocytes is decaying while the concentrations of the infected hepatocytes, and free viruses are increasing.Also we note that the trajectory tends to the stable infected steady state  1 .Figures 1-3 show also the effect of the saturation parameter  on the evolution of uninfected hepatocytes, infected hepatocytes and free viruses, respectively.It can be seen that, as  is increased, the evolution of the disease is postponed.
(ii) Treated Case.In this case, the treatment is designed via MPC strategy for the nominal and disturbed HBV infection models.In this case, we assume that the system is in the where the parameters () are uniformly distributed random numbers on [0, 1], and   = −  when the system states lie  + , the lower bound   has to be chosen as the following: to guarantee that the positive orthant R 3 + is positively invariant.Figure 4 shows that, when the MPC is applied, the concentration of uninfected hepatocytes is increasing.From Figures 5 and 6, we can see that the concentrations of infected hepatocytes and free virions are decaying during  the treatment.From Figures 4-6, we observe that when the MPC strategy is applied in the presence of bounded disturbances, the trajectory of the system tends to a ball around the uninfected steady state  0 and remains there.The drug doses  1 =  1 / 1 and  2 =  2 / 2 as functions of the time for MPC-(I), MPC-(II), and MPC-(III) cases are shown in Figures 7 and 8, respectively.

Conclusion
In this paper we have studied a mathematical model describing the HBV dynamics.The model has been incorporated to allow some additive disturbances.Under the effect of two types of drug therapies the HBV dynamics model is considered as a nonlinear control system, where the control input is defined to be dependent on the drug dose and drug efficiency.We have applied MPC method for determining the optimal treatment schedules and for stabilizing the HBV infection system around the uninfected steady state.The inherent robustness properties of the MPC have been established.

Figure 1 :
Figure 1: The evolution of uninfected hepatocytes for untreated case.

Figure 2 :
Figure 2: The evolution of infected hepatocytes for untreated case.

Figure 3 :
Figure 3: The evolution of free virus for untreated case.

Figure 4 :
Figure 4: The evolution of uninfected hepatocytes under MPC.

Figure 5 :Figure 6 :
Figure 5: The evolution of infected hepatocytes under MPC.

Figure 7 :Figure 8 :
Figure 7: The drug dose of the first type designed by MPC.

Table 1 :
The values of the parameters in the HBV dynamics model.