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

Authors: A. M. Elaiw; M. A. Alghamdi; Shaban Aly
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

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


---

## Body

## 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)T˙=s-dT+aT[1-TTmax]-βe-u1TV1+bV+w1,(2)I˙=βe-u1TV1+bV-δI+w2,(3)V˙=pe-u2I-cV+w3.
Here T, I, and V denote the concentration of uninfected hepatocytes, infected hepatocytes, and free virions, respectively. Uninfected hepatocytes are assumed to be produced at the constant rate sand die at the rate of dT. Uninfected hepatocytes can also be created by proliferation of existing T cells. Here we represent the proliferation by a logistic function in which a is the maximum proliferation rate of target cells, and Tmax is the T concentration at which proliferation shuts off. The rate of infection is given by saturation functional response βe-u1TV/(1+bV), where β is the infection rate constant characteristic of the infection efficiency and b>0 is constant. The death rate of infected hepatocytes is δI. Free virions are assumed to be produced from infected hepatocytes at the rate of pe-u2I and cV is the clearance rate of viral particles. Functions wi, i=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 e-u1, and the second is used to prevent the infected cells from producing new viruses and is represented by the function e-u2. Here, ui(t)=ψimi(t),i=1,2 are the control inputs, where the parameters ψ1 and ψ2 are the efficiencies of the drugs and m1(t) and m2(t) 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 (T,I,V)′ as the state vector and (u1,u2)′ 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(4)wi(t)=Ψi(T,I,V,u1,u2,ϑ),i=1,2,3,
where ϑ is the external disturbances and Ψi satisfies the following bound
(5)∥Ψi(T,I,V,u1,u2,ϑ)∥≤ϵi,ϵi≥0,i=1,2,3.
Now we show that under which conditions the nonnegative orthant ℝ+3 is positively invariant for (1)–(3):
(6)T˙|(T=0)=s+w1≥0,ifw1≥-s,I˙|(I=0)=βe-u1TV1+bV+w2≥0,ifw2≥-βe-u1TV1+bV,(T,V)≥0,V˙|(V=0)=pe-u2I+w3≥0,ifw3≥-pe-u2I,I≥0.
This means that the nonnegative orthant ℝ+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 dTmax>s+ϵ1. Then there exist positive numbers M1 and M2 for which the compact set
(7)Ω={(T,I,V):0≤T,I≤M1,0≤V≤M2}
is positively invariant.Proof.
From (1) we have
(8)T˙≤s-dT+aT[1-TTmax]+w1≤s-dT+aT[1-TTmax]+ϵ1,T˙|T=Tmax≤s-dTmax+ϵ1≤0.
This means that, even if T reaches Tmax, it should decrease.
LetTtot=T+I, then
(9)T˙tot=s-dT+aT[1-TTmax]-δI+w1+w2≤s-dT-δI+aT-aT2Tmax+ϵ1+ϵ2=s-dT-δI-aTmax[T-Tmax2]2+aTmax4+ϵ1+ϵ2≤s-σTtot+aTmax4+ϵ1+ϵ2,
where σ=min{d,δ}. Hence 0≤Ttot(t)≤(aTmax+4(s+ϵ1+ϵ2))/4σ for all t≥0 if Ttot(0)≤(aTmax+4(s+ϵ1+ϵ2))/4σ. It follows that 0≤T(t), I(t)≤M1 for all t≥0 if T(0)+I(0)≤M1, where M1=(aTmax+4(s+ϵ1+ϵ2))/4σ. On the other hand, V˙≤pM1-cV+ϵ3. Let M2=(pM1+ϵ3)/c, then we get that V(t)≤M2 for all t>0, whenever this inequality holds for t=0 and I(0)≤M1.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 uj(t)=u¯j,j=1,2, wi(t)=0, i=1,2,3, t≥0. System (1)–(3) has two steady states, uninfected steady state E0=(T0,0,0) and infected steady state E1=(T1,I1,V1), where
(10)T0=Tmax2a[a-d+(a-d)2+4asTmax],T1=Tmax2a[+(a-d-βe-u¯1b)2+4asTmax+4acδpe-u¯2bTmaxa-d-βe-u¯1b+(a-d-βe-u¯1b)2+4asTmax+4acδpe-u¯2bTmax],I1=cpe-u¯2V1,V1=1bcδ[pβe-(u¯1+u¯2)T1-cδ].
It has been shown in [22] that the infected steady state E1 exists if (pβe-(u¯1+u¯2)/δc)T1>1.The basic reproductive ratio for system (1)–(3) is given by
(11)R0(u¯1,u¯2)=pβe-(u¯1+u¯2)T1δc.Theorem 2 (see [22]).
For the nominal system (1)–(3) (i.e., wi(t)=0, i=1,2,3), if (pβe-(u¯1+u¯2)/δc)T0<1, then E0 is locally asymptotically stable.Proposition 3.
The nominal system (1)–(3) is locally asymptotically controllable to E0 with piecewise constant controllers.Proof.
Letu1(t)=u^1 and u2(t)=u^2 with u^1+u^2>uc, where
(12)uc=max{0,ln(pβT0δc)}.
Then (pβe-(u¯1+u¯2)/δc)T0<1; therefore, the corresponding trajectory will tend to E0 as t→∞.Theorem 4 (see [22]).
For the nominal system (1)–(3), if R0(u¯1,u¯2)<1, then E0 is globally asymptotically stable.Assumption 5.
There existu¯1≥0 and u¯2≥0 such that
(13)Tmaxpβe-(u¯1+u¯2)2aδc×[(a-d-βe-u¯1b)2+4asTmax+4acδpe-u¯2bTmaxa-d-βe-u¯1b+(a-d-βe-u¯1b)2+4asTmax+4acδpe-u¯2bTmax]<1.Proposition 6.
Suppose that Assumption5 is valid. Then the nominal system (1)–(3) is globally asymptotically controllable to E0 with piecewise constant controllers.Proof.
Letu1(t)=u¯1 and u2(t)=u¯2 such that u¯1 and u¯2 satisfy inequality (13), then R0(u¯1,u¯2)<1, and the corresponding trajectory will tend to E0 as t→∞.Theorem 7 (see [22]).
Consider the nominal system (1)–(3). Suppose that(i)
R0(u¯1,u¯2)>1,(ii)
δcR0(u¯1,u¯2)(R0(u¯1,u¯2)-1)+R02(u¯1,u¯2)((s/T1)+a(T1/Tmax))((s/T1)+a(T1/Tmax)+δ+c)>βe-u¯1δcV1/(δ+c).Then the positive steady state E1 is locally asymptotically stable.
Define(14)m=Tmax2a[a-d-βe-u¯1b+(a-d-βe-u¯1b)2+4asTmax].Theorem 8 (see [22]).
Consider the nominal system (1)–(3). Suppose that(i)
R0(u¯1,u¯2)>1,(ii)
δcR0(u¯1,u¯2)(R0(u¯1,u¯2)-1)+R02(u¯1,u¯2)((s/T1)+a(T1/Tmax))((s/T1)+a(T1/Tmax)+δ+c) > βe-u¯1δcV1/(δ+c).Then the positive steady state E1 is globally asymptotically stable provided that one of the following two assumptions holds:(i)
(Tmax/2a)(a-d)<m<T0<(Tmax/2a)(a-d+δ),(ii)
m>(Tmax/2a)(a-d+δ).


## 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, σ(s)>0 for all s>0 and it is strictly increasing. It is of class-𝒦∞ if it is of class-𝒦 and σ(s)→∞ when s→∞. A continuous function β:ℝ≥0×ℝ≥0→ℝ≥0 is of class-𝒦ℒ if β(s,τ) is of class-𝒦 in s for every τ≥0, it is strictly decreasing in τ for every s>0 and β(s,τ)→0 when τ→∞. Let ℬΔ={z∈ℝn:∥z∥≤Δ} will be used in ℝn.Consider a continuous-time nonlinear control system with additive disturbances given by(15)z˙(t)=f(z(t),u(t))+w(t),z(0)=z0,
where z(t)∈ℝn, u(t)∈U⊂ℝm, w(t)∈W⊂ℝp are the state, control input, and disturbances, respectively, f:ℝn×U→ℝn is continuous and Lipschitz continuous with respect to z in any compact set and f(0,0)=0,U is compact and 0∈U, W is compact and 0∈W.The control is taken to be a piecewise constant signal(16)u(t)=u(iτ)=:ui,fort∈[iτ,(i+1)τ),i∈ℕ,
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 jτm, j=0,1,…:
(17)yj:=z(jτm),j=0,1,…,
where τm is the measurement sampling period. We assume that τm=ℓτ for a fixed integer ℓ>0.For a given functionw:ℝ≥0→ℝn, we use the following notation: wτ[i]:={w(t),t∈[iτ,(i+1)τ]}, where i∈ℕ. We denote the norm ∥w∥∞:=ess.sups≥0∥w(s)∥. We assume that there exists μ>0 such that W⊂ℬμ. Let us define
(18)𝒲μ={w∈L[0,∞)∞:w(t)∈Wa.e.t∈[0,∞)with∥w∥∞≤μ},𝒲ℓρ={w(i)={wτ[iℓ],…,wτ[(i+1)ℓ-1]},w∈𝒲ρ,i=0,1,…}.
We will assume that there is a compact set 𝒳⊂ℝn containing the origin, that is positively invariant with respect to system (15) for any w(·)∈W and any piecewise constant controller u¯∈U. Let t↦ΦE(t,z,u¯,w¯(·)) denote the solution of (15) with given u¯, w¯, and z=z(0). Then the exact discrete-time model can be defined as
(19)zi+1=F~τE(zi,ui,wτ[i]),
where F~τE(z,u,wτ):=ΦE(τ;z,u,wτ).Letu(i)={u0(i),…,uℓ-1(i)}, w(i)={wτ[iℓ],…,wτ[(i+1)ℓ-1]}, and ℱℓE(ξ,u,w):=ΦE(τℓ,ξ,u,w), then the exact ℓ-step discrete-time model is given by
(20)ξi+1E=ℱℓE(ξiE,u(i),w(i)),ξ0E=z0.
We note that the exact discrete-time models (19) and (20) describe, respectively, the behavior of the system at the time instants kτ and kℓτ,k=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
(21)x˙(t)=f(x(t),u(t)),x(0)=z(0),
and its exact discrete-time model is given by
(22)xi+1E=FτE(xiE,ui).
We note that, since f is typically nonlinear, FτE 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
(23)xi+1A=Fτ,hA(xiA,ui),
where h 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 anh*>0 such that(i)
Fτ,hA(0,0)=0, Fτ,hA is continuous in both variables, uniformly continuous in h∈(0,h*], and Lipschitz continuous with respect to x in any compact set, uniformly in small h,(ii)
there exists aγ∈𝒦 such that
(24)∥FτE(x,u)-Fτ,hA(x,u)∥≤τγ(h),for all x∈𝒳, all u∈U, and h∈(0,h*].Assumption 10.
There exists anh*>0 such that the nominal exact discrete-time model (22) is practically asymptotically controllable from 𝒳 to the origin with piecewise constant controllers for all h∈(0,h*] (see e.g., [43] for the definition).For the solutions of (19), (22), and (23) with u={u0,u1,…},w={wτ[0],wτ[1],…} and x0 we will use the notations ΦiE(x0,u,w),ϕiE(x0,u), and ϕiA(x0,u), respectively.The following problem is to be solved: for givenτ and ℓ find a control strategy
(25)vh:𝒳⟶U×U×⋯×U︸ℓtimes,vh(x)={u0(x),…,uℓ-1(x)},
using the nominal approximate discrete-time model (23), to practically stabilize the exact discrete-time system (19).LetN∈ℕ with N≥ℓ be given. Let (23) be subject to the cost function
(26)Jτ,h(N,x,u)=∑i=0N-1‍τlh(xiA,ui)+g(xNA),
where u={u0,…,uN-1},xiA=ϕiA(x,u),i=0,1,…,N, lh and g are given functions satisfying the following assumptions.Assumption 11.
Let𝒳1=𝒳+ℬ1, (i) g:𝒳1→ℝ is continuous, positive definite radially unbounded, and Lipschitz continuous in any compact set,
(ii)lh(x,u) is continuous with respect to x and u, uniformly continuous in small h, and Lipschitz continuous in any compact set,
(iii) there exist anh*>0 and two class-𝒦∞ functions φ1 and φ2 such that the inequality
(27)φ1(∥x∥)≤lh(x,u)≤φ2(∥x∥)+φ2(∥u∥),
holds for all x∈𝒳1,u∈U and h∈(0,h*].Assumption 12.
There existh*>0 and η>0 such that for all x∈𝒢η={x:g(x)≤η} there exists a κ(x)∈U such that inequality
(28)τlh(x,κ(x))+g(Fτ,hA(x,κ(x)))≤g(x),
holds for all h∈(0,h*].We define the value function, which represents the optimal value of (26) for a given initial condition as follows
(29)VN(x)=inf{Jτ,h(N,x,u):ui∈U}.
Let u*={u0*,…,uN-1*} be the solution of this optimization problem, then the first ℓ elements of u* are applied at the state x, that is,
(30)vh(x)={u0*(x),…,uℓ-1*(x)}.
Let h0* denote the minimum of the values h* generated by Assumptions 9–12. Let Δx and Δu be such numbers that, for each x∈𝒳,u∈U, ∥x∥≤Δx, ∥u∥≤Δu.Theorem 13 (see [43]).
If Assumptions9–12 hold, then(i)
there exist anh1* with 0<h1*≤h0*, and a constant VmaxA independent of N, such that VN(x)≤VmaxA for all x∈𝒳,h∈(0,h1*] and N∈ℕ,(ii)
there exist constantsN*,LVand δV and functions σ1,σ2∈𝒦∞ such that for all x∈𝒳, N>N*, h∈(0,h1*] and i=1,2,…,ℓ,
(31)σ1(∥x∥)≤VN(x)≤σ2(∥x∥),VN(ϕiA(x,vh(x)))-VN(x)≤-τφ1(∥x∥).Moreover, for all x,y∈𝒳1 with ∥x-y∥≤δV,
(32)|VN(x)-VN(y)|≤LV∥x-y∥,
for all h∈(0,h1*].
Clearly𝒳⊂{x:VN(x)≤VmaxA}.Theorem 14 (see [38, 44]).
Suppose that Assumptions9–12 are valid and N is chosen such that N≥N*. Then, there exist β∈𝒦ℒ,θ∈𝒦∞, μ*>0 and for any δ>0 there exists an h*>0 such that for any x0∈𝒳, and h∈(0,h*] the trajectory of the ℓ-step exact discrete-time system
(33)ξi+1E=ℱℓE(ξiE,vh(ξiE),w(i)),ξ0E=x0,
with the ℓ-step MPC vh and w(i)∈𝒲ℓμ* satisfies
(34)∥ξiE∥≤β(∥x0∥,iℓτ)+θ(μ*)+δ,∀i≥0.


## 4. MPC for the HBV Model

In this section we apply the MPC method proposed in Section3. We will show that, with a suitable choice of N and functions g and lh, the assumptions of the previous section can be satisfied. Introduce new variables as z1=T-T0,z2=I,z3=V. In these new variables the model (1)–(3) takes the form of (15) with
(35)f(z,u)=(s-d(z1+T0)+a(z1+T0)[1-z1+T0Tmax]-βe-u1(z1+T0)z31+bz3βe-u1(z1+T0)z31+bz3-δz2pe-u2z2-cz3),
and w=(w1,w2,w3)′.Let the compact set𝒳 be defined as
(36)𝒳={(z1,z2,z3):-T0≤z1≤M1-T0,0≤z2≤M1,0≤z3≤M2},
where M1 and M2 are given in Proposition 1.With this definition,f satisfies all regularity assumptions, and according to Proposition 1, 𝒳 is positively invariant.To verify Assumptions11 and 12, we linearized the nominal system (35) around the origin in case of constant controllers, that is, u1(t)=u¯1>uc(1), u2(t)=u¯2>uc(2) with uc(1)+uc(2)=uc, where uc is given in Proposition 3. Let AC be the coefficient matrix of the linearized system and x=(T-T0,I,V)′. Then the discrete-time model for the linearized system is given by
(37)x(k+1)=eACτx(k).
Let the sampling period be chosen to be τ=0.5 and u¯1=u¯2=1. The running cost and the terminal cost can be chosen as
(38)lh(x,u)=0.1x′Qx+(u1-uc(1))2+(u2-uc(2))2,(39)g(x)=x′Px,
where P is a positive definite diagonal matrix and Q is a positive definite symmetric matrix satisfying the Lyapunov equation for the discrete-time system (37)
(40)Q=-(Aτ′PAτ-P),Aτ=eACτ.
Using the parameters given in Table 1, we have verified Assumption 5 numerically by fixing one controller, for example, u¯1 and solve the inequality (13) with respect to u¯2. We have found that for u¯1=1, then inequality (13) is satisfied when u¯2>0.20397.Table 1
The values of the parameters in the HBV dynamics model.
Parameters
d
s
δ
p
c
β
a
T
max
ψ
1
ψ
2
Values
0.072
36
0.5
5
3
0.001
0.108
1500
0.01
0.01From (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 P has been chosen through a series of numerical experiments P=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 Table1.We assumed that the state measurements were performed at the instantsjℓτ,j=0,1,… with ℓ=4 and the horizon length N was chosen to be N=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., wi(t)=0, i=1,2,3) without treatment (i.e., u1=u2=0). Using the parameters given in Table 1 and choose b=0.005, we have obtained (pβ/δc)T0=3.33333>1 and R0(0.0)=1.60987>1, moreover, all the conditions of Theorem 8 are satisfied. Therefore the uninfected steady state E0 is unstable and the infected steady state E1 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 V=0.001. Thus the initial conditions for the untreated case are T(0)=T0=1000, I(0)=0, and V(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 E1. Figures 1–3 show also the effect of the saturation parameter b on the evolution of uninfected hepatocytes, infected hepatocytes and free viruses, respectively. It can be seen that, as b 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, E1|(u1=u2=0). In this case, we choose the saturation parameter b=0.005. Thus the initial conditions for the MPC are given by T(0)=482.96, I(0)=73.19, and V(0)=121.98.The disturbances are simulated by(41)wi(t)∈[ηi,ϵi],wi(t)=wi(j)=ηi+(ϵi-ηi)r(j),t∈[jτ,(j+1)τ),i=1,2,3,j=0,1,…,
where the parameters r(j) are uniformly distributed random numbers on [0,1], and ηi=-ϵi when the system states lie in the interior of the positive orthant ℝ+3. At the boundary of ℝ+3, the lower bound ηi has to be chosen as the following:
(42)η1=max{-s,-ϵ1},η2=max{-e-u1β(z1+T0)z31+bz3,-ϵ2},η3=max{-pe-u2z2,-ϵ3},
to guarantee that the positive orthant ℝ+3 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):wi(t)=0, i=1,2,3,(2)
MPC-(II)wi(t)∈[ηi,ϵi], with ϵ1=20, ϵ2=5, ϵ3=10,(3)
MPC-(III)wi(t)∈[ηi,ϵi], with ϵ1=40, ϵ2=10, ϵ3=20.Figure4 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 E0 and remains there. The drug doses m1=u1/ψ1 and m2=u2/ψ2 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.



---

*Source: 712829-2013-04-18.xml*