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

Abstract

Modeling, analysis, and control of hepatitis B virus (HBV) infection have attracted the interests of mathematicians 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 these models can completely exhibit all that is observed clinically and account the full course of infection. Besides model inaccuracies that HBV dynamics models suffer from, some disturbances/uncertainties from different sources may 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.

1. Introduction

Hepatitis B virus (HBV) infection is among the most common causes of hepatitis and can result in serious liver diseases such as chronic hepatic insufficiency, hepatocellular carcinoma, and cirrhosis. Over the last decade, much collaborated effort involving biologists and mathematicians has been devoted towards designing mathematical models of HBV dynamics [1–11], human immunodeficiency virus (HIV) [12–20], and hepatitis C virus (HCV) [21]. Mathematical modeling and model analysis of the HBV dynamics are important for exploring possible mechanisms and dynamical behaviors of the viral infection process, estimating key parameter values, and guiding development of efficient antiviral drug therapies. Stability analysis of HBV dynamics models has been studied by many authors (see e.g., [22–32]).

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

2. 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 𝑇̇=𝑠−𝑑𝑇+𝑎𝑇1−𝑇𝑇_max−𝛽𝑒^−𝑢1𝑇𝑉1+𝑏𝑉+𝑤₁,𝐼̇=𝛽𝑒^−𝑢1𝑇𝑉1+𝑏𝑉−𝛿𝐼+𝑤₂,𝑉̇=𝑝𝑒^−𝑢2𝐼−𝑐𝑉+𝑤₃. (1)(2)(3) 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 𝜓₁ and 𝜓₂ are the efficiencies of the drugs and 𝑚₁(𝑡) and 𝑚₂(𝑡) 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 (𝑢₁,𝑢₂)^′ 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 𝑤_𝑖(𝑡)=Ψ_𝑖𝑇,𝐼,𝑉,𝑢₁,𝑢₂,𝜗,𝑖=1,2,3, (4) where 𝜗 is the external disturbances and Ψ_𝑖 satisfies the following bound Ψ_𝑖𝑇,𝐼,𝑉,𝑢₁,𝑢₂,𝜗≤𝜖_𝑖,𝜖_𝑖≥0,𝑖=1,2,3. (5) Now we show that under which conditions the nonnegative orthant ℝ³₊ is positively invariant for (1)–(3): 𝑇̇|_(𝑇=0)=𝑠+𝑤₁≥0,if𝑤₁𝐼|̇≥−𝑠,_(𝐼=0)=𝛽𝑒^−𝑢1𝑇𝑉1+𝑏𝑉+𝑤₂≥0,if𝑤₂≥−𝛽𝑒^−𝑢1𝑇𝑉,()𝑉̇|1+𝑏𝑉𝑇,𝑉≥0,_(𝑉=0)=𝑝𝑒^−𝑢2𝐼+𝑤₃≥0,if𝑤₃≥−𝑝𝑒^−𝑢2𝐼,𝐼≥0. (6) This means that the nonnegative orthant ℝ³₊ 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>𝑠+𝜖₁. Then there exist positive numbers 𝑀₁ and 𝑀₂ for which the compact set Ω=(𝑇,𝐼,𝑉)∶0≤𝑇,𝐼≤𝑀₁,0≤𝑉≤𝑀₂ (7) is positively invariant.

Proof. From (1) we have 𝑇̇≤𝑠−𝑑𝑇+𝑎𝑇1−𝑇𝑇_max+𝑤₁𝑇≤𝑠−𝑑𝑇+𝑎𝑇1−𝑇_max+𝜖₁,𝑇̇|_𝑇=𝑇max≤𝑠−𝑑𝑇_max+𝜖₁≤0. (8) This means that, even if 𝑇 reaches 𝑇_max, it should decrease.
Let 𝑇_tot=𝑇+𝐼, then 𝑇̇_tot=𝑠−𝑑𝑇+𝑎𝑇1−𝑇𝑇_max−𝛿𝐼+𝑤₁+𝑤₂≤𝑠−𝑑𝑇−𝛿𝐼+𝑎𝑇−𝑎𝑇²𝑇_max+𝜖₁+𝜖₂𝑎=𝑠−𝑑𝑇−𝛿𝐼−𝑇_max𝑇𝑇−^max2²+𝑎𝑇_max4+𝜖₁+𝜖₂≤𝑠−𝜎𝑇_tot+𝑎𝑇_max4+𝜖₁+𝜖₂, (9) where 𝜎=min{𝑑,𝛿}. Hence 0≤𝑇_tot(𝑡)≤(𝑎𝑇_max+4(𝑠+𝜖₁+𝜖₂))/4𝜎 for all 𝑡≥0 if 𝑇_tot(0)≤(𝑎𝑇_max+4(𝑠+𝜖₁+𝜖₂))/4𝜎. It follows that 0≤𝑇(𝑡), 𝐼(𝑡)≤𝑀₁ for all 𝑡≥0 if 𝑇(0)+𝐼(0)≤𝑀₁, where 𝑀₁=(𝑎𝑇_max+4(𝑠+𝜖₁+𝜖₂))/4𝜎. On the other hand, 𝑉̇≤𝑝𝑀₁−𝑐𝑉+𝜖₃. Let 𝑀₂=(𝑝𝑀₁+𝜖₃)/𝑐, then we get that 𝑉(𝑡)≤𝑀₂ for all 𝑡>0, whenever this inequality holds for 𝑡=0 and 𝐼(0)≤𝑀₁.

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.

Let us compute the steady states of system (1)–(3) under constant controller in the absence of the disturbances that is for 𝑢_𝑗(𝑡)=𝑢_𝑗,𝑗=1,2, 𝑤_𝑖(𝑡)=0, 𝑖=1,2,3, 𝑡≥0. System (1)–(3) has two steady states, uninfected steady state 𝐸₀=(𝑇₀,0,0) and infected steady state 𝐸₁=(𝑇₁,𝐼₁,𝑉₁), where 𝑇₀=𝑇_max2𝑎𝑎−𝑑+(𝑎−𝑑)²+4𝑎𝑠𝑇_max,𝑇₁=𝑇_max2𝑎𝑎−𝑑−𝛽𝑒^−𝑢1𝑏+𝑎−𝑑−𝛽𝑒^−𝑢1𝑏²+4𝑎𝑠𝑇_max+4𝑎𝑐𝛿𝑝𝑒^−𝑢2𝑏𝑇_max,𝐼₁=𝑐𝑝𝑒^−𝑢2𝑉₁,𝑉₁=1𝑏𝑐𝛿𝑝𝛽𝑒^{−(𝑢1+𝑢2)}𝑇₁−𝑐𝛿. (10) It has been shown in [22] that the infected steady state 𝐸₁ exists if (𝑝𝛽𝑒^{−(𝑢1+𝑢2)}/𝛿𝑐)𝑇₁>1.

The basic reproductive ratio for system (1)–(3) is given by 𝑅₀𝑢₁,𝑢₂=𝑝𝛽𝑒^{−(𝑢1+𝑢2)}𝑇₁.𝛿𝑐 (11)

Theorem 2 (see [22]). For the nominal system (1)–(3) (i.e., 𝑤_𝑖(𝑡)=0, 𝑖=1,2,3), if (𝑝𝛽𝑒^{−(𝑢1+𝑢2)}/𝛿𝑐)𝑇₀<1, then 𝐸₀ is locally asymptotically stable.

Proposition 3. The nominal system (1)–(3) is locally asymptotically controllable to 𝐸₀ with piecewise constant controllers.

Proof. Let 𝑢₁(𝑡)=𝑢₁ and 𝑢₂(𝑡)=𝑢₂ with 𝑢₁+𝑢₂>𝑢_𝑐, where 𝑢_𝑐=max0,ln𝑝𝛽𝑇₀.𝛿𝑐 (12) Then (𝑝𝛽𝑒^{−(𝑢1+𝑢2)}/𝛿𝑐)𝑇₀<1; therefore, the corresponding trajectory will tend to 𝐸₀ as 𝑡→∞.

Theorem 4 (see [22]). For the nominal system (1)–(3), if 𝑅₀(𝑢₁,𝑢₂)<1, then 𝐸₀ is globally asymptotically stable.

Assumption 5. There exist 𝑢₁≥0 and 𝑢₂≥0 such that 𝑇_max𝑝𝛽𝑒^{−(𝑢1+𝑢2)}×2𝑎𝛿𝑐𝑎−𝑑−𝛽𝑒^−𝑢1𝑏+𝑎−𝑑−𝛽𝑒^−𝑢1𝑏²+4𝑎𝑠𝑇_max+4𝑎𝑐𝛿𝑝𝑒^−𝑢2𝑏𝑇_max<1. (13)

Proposition 6. Suppose that Assumption 5 is valid. Then the nominal system (1)–(3) is globally asymptotically controllable to 𝐸₀ with piecewise constant controllers.

Proof. Let 𝑢₁(𝑡)=𝑢₁ and 𝑢₂(𝑡)=𝑢₂ such that 𝑢₁ and 𝑢₂ satisfy inequality (13), then 𝑅₀(𝑢₁,𝑢₂)<1, and the corresponding trajectory will tend to 𝐸₀ as 𝑡→∞.

Theorem 7 (see [22]). Consider the nominal system (1)–(3). Suppose that(i)𝑅₀(𝑢₁,𝑢₂)>1, (ii)𝛿𝑐𝑅₀(𝑢₁,𝑢₂)(𝑅₀(𝑢₁,𝑢₂)−1)+𝑅²₀(𝑢₁,𝑢₂)((𝑠/𝑇₁)+𝑎(𝑇₁/𝑇_max))((𝑠/𝑇₁)+𝑎(𝑇₁/𝑇_max)+𝛿+𝑐)>𝛽𝑒^−𝑢1𝛿𝑐𝑉₁/(𝛿+𝑐). Then the positive steady state 𝐸₁ is locally asymptotically stable.
Define 𝑇𝑚=^max2𝑎𝑎−𝑑−𝛽𝑒^−𝑢1𝑏+𝑎−𝑑−𝛽𝑒^−𝑢1𝑏²+4𝑎𝑠𝑇_max. (14)

Theorem 8 (see [22]). Consider the nominal system (1)–(3). Suppose that(i)𝑅₀(𝑢₁,𝑢₂)>1, (ii)𝛿𝑐𝑅₀(𝑢₁,𝑢₂)(𝑅₀(𝑢₁,𝑢₂)−1)+𝑅²₀(𝑢₁,𝑢₂)((𝑠/𝑇₁)+𝑎(𝑇₁/𝑇_max))((𝑠/𝑇₁)+𝑎(𝑇₁/𝑇_max)+𝛿+𝑐) > 𝛽𝑒^−𝑢1𝛿𝑐𝑉₁/(𝛿+𝑐).Then the positive steady state 𝐸₁ is globally asymptotically stable provided that one of the following two assumptions holds:(i)(𝑇_max/2𝑎)(𝑎−𝑑)<𝑚<𝑇₀<(𝑇_max/2𝑎)(𝑎−𝑑+𝛿), (ii) 𝑚>(𝑇_max/2𝑎)(𝑎−𝑑+𝛿).

3. 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 ℝ and ℕ. The symbol ℝ_≥0 denotes the set of nonnegative real numbers. A continuous function 𝜎∶ℝ_≥0→ℝ_≥0 is of class-𝒦 if 𝜎(0)=0, 𝜎(𝑠)>0 for all 𝑠>0 and it is strictly increasing. It is of class-𝒦_∞ if it is of class-𝒦 and 𝜎(𝑠)→∞ when 𝑠→∞. A continuous function 𝛽∶ℝ_≥0×ℝ_≥0→ℝ_≥0 is of class-𝒦ℒ if 𝛽(𝑠,𝜏) is of class-𝒦 in 𝑠 for every 𝜏≥0, it is strictly decreasing in 𝜏 for every 𝑠>0 and 𝛽(𝑠,𝜏)→0 when 𝜏→∞. Let ℬ_Δ={𝑧∈ℝ^𝑛∶‖𝑧‖≤Δ} will be used in ℝ^𝑛.

Consider a continuous-time nonlinear control system with additive disturbances given by 𝑧̇(𝑡)=𝑓(𝑧(𝑡),𝑢(𝑡))+𝑤(𝑡),𝑧(0)=𝑧₀, (15) where 𝑧(𝑡)∈ℝ^𝑛, 𝑢(𝑡)∈𝑈⊂ℝ^𝑚, 𝑤(𝑡)∈𝑊⊂ℝ^𝑝 are the state, control input, and disturbances, respectively, 𝑓∶ℝ^𝑛×𝑈→ℝ^𝑛 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 𝑢(𝑡)=𝑢(𝑖𝜏)=∶𝑢_𝑖,for𝑡∈[𝑖𝜏,(𝑖+1)𝜏),𝑖∈ℕ, (16) 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,…: 𝑦_𝑗∶=𝑧𝑗𝜏^𝑚,𝑗=0,1,…, (17) where 𝜏^𝑚 is the measurement sampling period. We assume that 𝜏^𝑚=ℓ𝜏 for a fixed integer ℓ>0.

For a given function 𝑤∶ℝ_≥0→ℝ^𝑛, we use the following notation: 𝑤_𝜏[𝑖]∶={𝑤(𝑡),𝑡∈[𝑖𝜏,(𝑖+1)𝜏]}, where 𝑖∈ℕ. We denote the norm ‖𝑤‖_∞∶=ess.sup_𝑠≥0‖𝑤(𝑠)‖. We assume that there exists 𝜇>0 such that 𝑊⊂ℬ_𝜇. Let us define 𝒲_𝜇=𝑤∈𝐿^∞_[)0,∞(𝑡)[)‖𝑤‖∶𝑤∈𝑊a.e.𝑡∈0,∞with_∞,𝒲≤𝜇_𝜌ℓ=𝐰^(𝑖)=𝑤_𝜏[𝑖ℓ],…,𝑤_𝜏[(𝑖+1)ℓ−1],𝑤∈𝒲_𝜌,𝑖=0,1,…. (18) We will assume that there is a compact set 𝒳⊂ℝ^𝑛 containing the origin, that is positively invariant with respect to system (15) for any 𝑤(⋅)∈𝑊 and any piecewise constant controller 𝑢∈𝑈. Let 𝑡↦Φ^𝐸(𝑡,𝑧,𝑢,𝑤(⋅)) denote the solution of (15) with given 𝑢, 𝑤, and 𝑧=𝑧(0). Then the exact discrete-time model can be defined as 𝑧_𝑖+1=𝐹^∼𝐸_𝜏𝑧_𝑖,𝑢_𝑖,𝑤_𝜏[𝑖], (19) where 𝐹^∼𝐸_𝜏(𝑧,𝑢,𝑤_𝜏)∶=Φ^𝐸(𝜏;𝑧,𝑢,𝑤_𝜏).

Let 𝐮^(𝑖)={𝑢₀^(𝑖),…,𝑢^(𝑖)_ℓ−1}, 𝐰^(𝑖)={𝑤_𝜏[𝑖ℓ],…,𝑤_𝜏[(𝑖+1)ℓ−1]}, and ℱ^𝐸_ℓ(𝜉,𝐮,𝐰)∶=Φ^𝐸(𝜏ℓ,𝜉,𝐮,𝐰), then the exact ℓ-step discrete-time model is given by 𝜉^𝐸_𝑖+1=ℱ^𝐸_ℓ𝜉^𝐸_𝑖,𝐮^(𝑖),𝐰^(𝑖),𝜉^𝐸₀=𝑧₀. (20) We note that the exact discrete-time models (19) and (20) describe, respectively, the behavior of the system at the time instants 𝑘𝜏 and 𝑘ℓ𝜏,𝑘=0,1,….

In this work, the construction of multirate MPC is based on the nominal prediction and only small disturbances are allowed. The nominal system of (15) is given by 𝑥̇(𝑡)=𝑓(𝑥(𝑡),𝑢(𝑡)),𝑥(0)=𝑧(0), (21) and its exact discrete-time model is given by 𝑥^𝐸_𝑖+1=𝐹^𝐸_𝜏𝑥^𝐸_𝑖,𝑢_𝑖. (22) We note that, since 𝑓 is typically nonlinear, 𝐹^𝐸_𝜏 in (22) is not known in most cases, therefore the controller design can be carried out by means of the nominal approximate discrete-time model 𝑥^𝐴_𝑖+1=𝐹^𝐴_𝜏,ℎ𝑥^𝐴_𝑖,𝑢_𝑖, (23) where ℎ is a modeling parameter, which is typically the step size of the underlying numerical method. The applied numerical scheme approximation has to ensure the closeness of the exact models in the following sense.

Assumption 9. There exists an ℎ^∗>0 such that(i)𝐹^𝐴_𝜏,ℎ(0,0)=0, 𝐹^𝐴_𝜏,ℎ is continuous in both variables, uniformly continuous in ℎ∈(0,ℎ^∗], and Lipschitz continuous with respect to 𝑥 in any compact set, uniformly in small ℎ,(ii)there exists a 𝛾∈𝒦 such that 𝐹^𝐸_𝜏(𝑥,𝑢)−𝐹^𝐴_𝜏,ℎ(𝑥,𝑢)≤𝜏𝛾(ℎ), (24)for all 𝑥∈𝒳, all 𝑢∈𝑈, and ℎ∈(0,ℎ^∗].

Assumption 10. There exists an ℎ^∗>0 such that the nominal exact discrete-time model (22) is practically asymptotically controllable from 𝒳 to the origin with piecewise constant controllers for all ℎ∈(0,ℎ^∗] (see e.g., [43] for the definition).

For the solutions of (19), (22), and (23) with 𝐮={𝑢₀,𝑢₁,…},𝐰={𝑤_𝜏[0],𝑤_𝜏[1],…} and 𝑥₀ we will use the notations Φ^𝐸_𝑖(𝑥₀,𝐮,𝐰),𝜙^𝐸_𝑖(𝑥₀,𝐮), and 𝜙^𝐴_𝑖(𝑥₀,𝐮), respectively.

The following problem is to be solved: for given 𝜏 and ℓ find a control strategy 𝐯_ℎ∶𝒳⟶𝑈×𝑈×⋯×𝑈_ℓtimes,𝐯_ℎ(𝑥)=𝑢₀(𝑥),…,𝑢_ℓ−1(𝑥), (25) using the nominal approximate discrete-time model (23), to practically stabilize the exact discrete-time system (19).

Let 𝑁∈ℕ with 𝑁≥ℓ be given. Let (23) be subject to the cost function 𝐽_𝜏,ℎ(𝑁,𝑥,𝐮)=^𝑁−1_𝑖=0𝜏𝑙_ℎ𝑥^𝐴_𝑖,𝑢_𝑖+𝑔𝑥^𝐴_𝑁, (26) where 𝐮={𝑢₀,…,𝑢_𝑁−1},𝑥^𝐴_𝑖=𝜙^𝐴_𝑖(𝑥,𝐮),𝑖=0,1,…,𝑁, 𝑙_ℎ and 𝑔 are given functions satisfying the following assumptions.

Assumption 11. Let 𝒳₁=𝒳+ℬ₁, (i) 𝑔∶𝒳₁→ℝ 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-𝒦_∞ functions 𝜑₁ and 𝜑₂ such that the inequality 𝜑₁(‖𝑥‖)≤𝑙_ℎ(𝑥,𝑢)≤𝜑₂(‖𝑥‖)+𝜑₂(‖𝑢‖), (27) holds for all 𝑥∈𝒳₁,𝑢∈𝑈 and ℎ∈(0,ℎ^∗].

Assumption 12. There exist ℎ^∗>0 and 𝜂>0 such that for all 𝑥∈𝒢_𝜂={𝑥∶𝑔(𝑥)≤𝜂} there exists a 𝜅(𝑥)∈𝑈 such that inequality 𝜏𝑙_ℎ(𝑥,𝜅(𝑥))+𝑔𝐹^𝐴_𝜏,ℎ(𝑥,𝜅(𝑥))≤𝑔(𝑥), (28) holds for all ℎ∈(0,ℎ^∗].

We define the value function, which represents the optimal value of (26) for a given initial condition as follows 𝑉_𝑁(𝑥)=inf𝐽_𝜏,ℎ(𝑁,𝑥,𝐮)∶𝑢_𝑖∈𝑈. (29) Let 𝐮^∗={𝑢^∗₀,…,𝑢^∗_𝑁−1} be the solution of this optimization problem, then the first ℓ elements of 𝐮^∗ are applied at the state 𝑥, that is, 𝐯_ℎ(𝑥)=𝑢^∗₀(𝑥),…,𝑢^∗_ℓ−1(𝑥). (30) Let ℎ^∗₀ denote the minimum of the values ℎ^∗ generated by Assumptions 9–12. Let Δ_𝑥 and Δ_𝑢 be such numbers that, for each 𝑥∈𝒳,𝑢∈𝑈, ‖𝑥‖≤Δ_𝑥, ‖𝑢‖≤Δ_𝑢.

Theorem 13 (see [43]). If Assumptions 9–12 hold, then(i)there exist an ℎ^∗₁ with 0<ℎ^∗₁≤ℎ^∗₀, and a constant 𝑉^𝐴_max independent of 𝑁, such that 𝑉_𝑁(𝑥)≤𝑉^𝐴_max for all 𝑥∈𝒳,ℎ∈(0,ℎ^∗₁] and 𝑁∈ℕ,(ii)there exist constants 𝑁^∗,𝐿_𝑉and 𝛿_𝑉 and functions 𝜎₁,𝜎₂∈𝒦_∞ such that for all 𝑥∈𝒳, 𝑁>𝑁^∗, ℎ∈(0,ℎ^∗₁] and 𝑖=1,2,…,ℓ, 𝜎₁(‖𝑥‖)≤𝑉_𝑁(𝑥)≤𝜎₂𝑥,𝑉(‖‖)_𝑁𝜙^𝐴_𝑖𝑥,𝐯_ℎ(𝑥)−𝑉_𝑁(𝑥)≤−𝜏𝜑₁(‖𝑥‖). (31)Moreover, for all 𝑥,𝑦∈𝒳₁ with ‖𝑥−𝑦‖≤𝛿_𝑉, 𝑉_𝑁(𝑥)−𝑉_𝑁𝑦≤𝐿_𝑉𝑥−𝑦, (32) for all ℎ∈(0,ℎ^∗₁].
Clearly 𝒳⊂{𝑥∶𝑉_𝑁(𝑥)≤𝑉^𝐴_max}.

Theorem 14 (see [38, 44]). Suppose that Assumptions 9–12 are valid and 𝑁 is chosen such that 𝑁≥𝑁^∗. Then, there exist 𝛽∈𝒦ℒ,𝜃∈𝒦_∞, 𝜇^∗>0 and for any 𝛿>0 there exists an ℎ^∗>0 such that for any 𝑥₀∈𝒳, and ℎ∈(0,ℎ^∗] the trajectory of the ℓ-step exact discrete-time system 𝜉^𝐸_𝑖+1=ℱ^𝐸_ℓ𝜉^𝐸_𝑖,𝐯_ℎ𝜉^𝐸_𝑖,𝐰^(𝑖),𝜉^𝐸₀=𝑥₀, (33) with the ℓ-step MPC 𝐯_ℎ and 𝐰^(𝑖)∈𝒲^𝜇∗_ℓ satisfies 𝜉^𝐸_𝑖≤𝛽𝑥₀,𝑖ℓ𝜏+𝜃𝜇^∗+𝛿,∀𝑖≥0. (34)

4. 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 𝑧₁=𝑇−𝑇₀,𝑧₂=𝐼,𝑧₃=𝑉. In these new variables the model (1)–(3) takes the form of (15) with 𝑓(𝑧,𝑢)=_{𝑧𝑠−𝑑1+𝑇0𝑧+𝑎1+𝑇0𝑧1−1+𝑇0𝑇max−𝛽𝑒−𝑢1𝑧1+𝑇0𝑧31+𝑏𝑧3𝛽𝑒−𝑢1𝑧1+𝑇0𝑧31+𝑏𝑧3−𝛿𝑧2𝑝𝑒−𝑢2𝑧2−𝑐𝑧3}, (35) and 𝑤=(𝑤₁,𝑤₂,𝑤₃)^′.

Let the compact set 𝒳 be defined as 𝒳=𝑧₁,𝑧₂,𝑧₃∶−𝑇₀≤𝑧₁≤𝑀₁−𝑇₀,0≤𝑧₂≤𝑀₁,0≤𝑧₃≤𝑀₂, (36) where 𝑀₁ and 𝑀₂ are given in Proposition 1.

With this definition, 𝑓 satisfies all regularity assumptions, and according to Proposition 1, 𝒳 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, 𝑢₁(𝑡)=𝑢₁>𝑢_𝑐⁽¹⁾, 𝑢₂(𝑡)=𝑢₂>𝑢_𝑐⁽²⁾ with 𝑢_𝑐⁽¹⁾+𝑢_𝑐⁽²⁾=𝑢_𝑐, where 𝑢_𝑐 is given in Proposition 3. Let 𝐴_𝐶 be the coefficient matrix of the linearized system and 𝑥=(𝑇−𝑇₀,𝐼,𝑉)^′. Then the discrete-time model for the linearized system is given by 𝑥(𝑘+1)=𝑒^𝐴𝐶𝜏𝑥(𝑘). (37) Let the sampling period be chosen to be 𝜏=0.5 and 𝑢₁=𝑢₂=1. The running cost and the terminal cost can be chosen as 𝑙_ℎ(𝑥,𝑢)=0.1𝑥^′𝑄𝑥+𝑢₁−𝑢_𝑐⁽¹⁾²+𝑢₂−𝑢_𝑐⁽²⁾²,𝑔(𝑥)=𝑥^′𝑃𝑥, (38)(39) where 𝑃 is a positive definite diagonal matrix and 𝑄 is a positive definite symmetric matrix satisfying the Lyapunov equation for the discrete-time system (37) 𝑄=−𝐴^′_𝜏𝑃𝐴_𝜏−𝑃,𝐴_𝜏=𝑒^𝐴𝐶𝜏. (40) Using the parameters given in Table 1, we have verified Assumption 5 numerically by fixing one controller, for example, 𝑢₁ and solve the inequality (13) with respect to 𝑢₂. We have found that for 𝑢₁=1, then inequality (13) is satisfied when 𝑢₂>0.20397.

Table 1 

The values of the parameters in the HBV dynamics model.

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 𝒳. Thus all assumptions of the proposed method can be satisfied with suitable choice of the parameters of the MPC method.

5. 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., 𝑢₁=𝑢₂=0). Using the parameters given in Table 1 and choose 𝑏=0.005, we have obtained (𝑝𝛽/𝛿𝑐)𝑇₀=3.33333>1 and 𝑅₀(0.0)=1.60987>1, moreover, all the conditions of Theorem 8 are satisfied. Therefore the uninfected steady state 𝐸₀ is unstable and the infected steady state 𝐸₁ 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)=𝑇₀=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 𝐸₁. 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.

Figure 1 

The evolution of uninfected hepatocytes for untreated case.

Figure 2 

The evolution of infected hepatocytes for untreated case.

Figure 3 

The evolution of free virus for untreated case.

(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 infected steady state before initiating the treatment, that is, 𝐸₁|_{(𝑢1=𝑢2=0)}. In this case, we choose the saturation parameter 𝑏=0.005. Thus the initial conditions for the MPC are given by 𝑇(0)=482.96, 𝐼(0)=73.19, and 𝑉(0)=121.98.

The disturbances are simulated by 𝑤_𝑖(𝑡)∈𝜂_𝑖,𝜖_𝑖,𝑤_𝑖(𝑡)=𝑤_𝑖𝑗=𝜂_𝑖+𝜖_𝑖−𝜂_𝑖𝑟𝑗,𝑡∈𝑗𝜏,𝑗+1𝜏,𝑖=1,2,3,𝑗=0,1,…, (41) where the parameters 𝑟(𝑗) are uniformly distributed random numbers on [0,1], and 𝜂_𝑖=−𝜖_𝑖 when the system states lie in the interior of the positive orthant ℝ³₊. At the boundary of ℝ³₊, the lower bound 𝜂_𝑖 has to be chosen as the following: 𝜂₁=max−𝑠,−𝜖₁,𝜂₂=max−𝑒^−𝑢1𝛽𝑧₁+𝑇₀𝑧₃1+𝑏𝑧₃,−𝜖₂,𝜂₃=max−𝑝𝑒^−𝑢2𝑧₂,−𝜖₃, (42) to guarantee that the positive orthant ℝ³₊ is positively invariant.

We apply the MPC to the HBV model with and without disturbances. The following cases are shown in the simulation results: (1)MPC-(I): 𝑤_𝑖(𝑡)=0, 𝑖=1,2,3,(2)MPC-(II) 𝑤_𝑖(𝑡)∈[𝜂_𝑖,𝜖_𝑖], with 𝜖₁=20, 𝜖₂=5, 𝜖₃=10,(3)MPC-(III) 𝑤_𝑖(𝑡)∈[𝜂_𝑖,𝜖_𝑖], with 𝜖₁=40, 𝜖₂=10, 𝜖₃=20.

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 𝐸₀ and remains there. The drug doses 𝑚₁=𝑢₁/𝜓₁ and 𝑚₂=𝑢₂/𝜓₂ as functions of the time for MPC-(I), MPC-(II), and MPC-(III) cases are shown in Figures 7 and 8, respectively.

Figure 4 

The evolution of uninfected hepatocytes under MPC.

Figure 5 

The evolution of infected hepatocytes under MPC.

Figure 6 

The evolution of free virus under MPC.

Figure 7 

The drug dose of the first type designed by MPC.

Figure 8 

The drug dose of the second type designed by MPC.

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

Acknowledgments

This project was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, under Grant no. 408/130/1431. The authors, therefore, acknowledge with thanks DSR technical and financial support. The authors are grateful to the anonymous reviewers for constructive suggestions and valuable comments, which improved the quality of the paper.