# Optimal Antiviral Treatment Strategies of HBV Infection Model with Logistic Hepatocyte Growth

Authors: Hassan Laarabi; Abdelhadi Abta; Mostafa Rachik; Jamal Bouyaghroumni
Journal: ISRN Biomathematics
 (2013)

---

## Abstract

This study considers an optimal therapy strategy for HBV infection by incorporating two controls laws into a previous hepatitis B viral infection model with logistic hepatocyte growth. Our goal is to
maximize the number of healthy cells and to minimize the cost of the therapy. In this context, the existence of an optimal control is proved. The optimal control is obtained by solving the optimality system which was composed of three nonlinear ODEs with initial conditions and three nonlinear adjoint ODEs with transversality conditions. The results were analysed and interpreted numerically using MATLAB.


---

## Body

## 1. Introduction

Hepatitis B virus (HBV) infection is a major global public-health problem. Of the approximately 2 billion people who have been infected worldwide, more than 350 million are chronic carriers of HBV. Approximately 15–40% of infected patients will develop cirrhosis, liver failure, or hepatocellular carcinoma (HCC) which is the fifth most frequent cancer, killing 300 000–500 000 people each year. The treatment of chronic infection may be necessary to reduce the risk of cirrhosis and liver cancer. Chronically, infected individuals with persistently elevated serum alanine aminotransferase, a marker of liver damage, and HBV DNA levels are candidates for therapy. Although none of the available drugs can clear the infection, they can stop the virus from replicating and prevent liver damage such as cirrhosis and liver cancer. Treatments include antiviral drugs such as lamivudine, adefovir, tenofovir, telbivudine, and entecavir and immune system modulators such as interferon alpha. However, some individuals are much more likely to respond than others; this might be because of the genotype of the infecting virus or the patient’s heredity. The treatment works by reducing the viral load (the amount of virus particles as measured in the blood), which in turn reduces viral replication in the liver.Mathematical modelling and analysis of hepatitis B infections have been explored extensively over the last decade. One of the earliest models was proposed by Nowak et al. [1] who modeled the infection of healthy hepatocytes by free virions as a mass action process. This makes the viral basic reproduction number dependent on the homeostatic liver size. Gourley et al. [2], Min et al. [3], Eikenberry et al. [4], Ciupe et al. [5, 6], and Hews et al. [7] all replace this mass-action process with a standard incidence function to increase the richness of the dynamics.Optimal control theory has found wide-ranging applications in biological and ecological problems [8]. In biomedical problems, techniques from control theory are of great use in developing optimal therapeutic strategies. The treatment regimen is usually taken to be the control variable, with the aim of minimizing the detrimental effects of the medical condition.We consider the mathematical model for HBV developed by Hews et al. [7], who modeled the constant infusion of healthy hepatocytes with a logistic growth term and the mass action infection term by a standard incidence function. The dynamics of this model are governed by the following equations:
(1)dxdt=rx(t)(1-T(t)K)-βv(t)x(t)T(t),x(0)=x0≥0,dydt=βv(t)x(t)T(t)-ay(t),y(0)=y0≥0,dvdt=γy(t)-μv(t),v(0)=v0≥0,
with
(2)T(t)=x(t)+y(t),
where x(t) is the mass of healthy hepatocytes, y(t) is the mass of infected hepatocytes, and v(t) is the mass of free virions. Healthy hepatocytes grow at a rate that depends on the homeostatic liver size, K, at a maximum per capita proliferation rate r. Hepatocytes become infected at maximum rate β and die at rate a. Free virions are created by infected hepatocytes at per capita rate γ, and virions either disintegrate or are destroyed by the immune system at rate μ. All parameters in the mathematical model are strictly positive.


## 2. The Optimal Control Problem

We introduce a controlu(t)=(u1(t),u2(t)) to the previously mentioned model simulating the antiviral therapy. Hence, the previous model becomes
(3)dxdt=rx(t)(1-T(t)K)-(1-u1(t))βv(t)x(t)T(t),x(0)=x0≥0,dydt=(1-u1(t))βv(t)x(t)T(t)-ay(t),y(0)=y0≥0,dvdt=(1-u2(t))γy(t)-μv(t),v(0)=v0≥0.
The control functions u1(t) and u2(t) are bounded Lebesgue integrable functions. The control u1(t) represents the efficiency of drug therapy in blocking new infection, and the control u2(t) represents the efficiency of drug therapy in inhibiting viral production.Ifu2=1, the inhibition is 100% effective, whereas if u2=0, there is no inhibition.The problem is to maximize the objective functional(4)J(u1(t),u2(t))=∫0tend‍[x(t)-(A12u12(t)+A22u22(t))]dt,
where the parameters A1≥0 and A2≥0 are based on the benefits and costs of the treatment. Our target is to maximize the objective functional defined in (4) by increasing the number of the uninfected cells, decreasing the viral load, and minimizing the cost of treatment. In other words, we are seeking an optimal control pair (u1*,u2*) such that
(5)J(u1*,u2*)=max{J(u1,u2)(u1,u2)∈𝒰},
where 𝒰 is the control set defined by
(6)𝒰={[0,tf]u=(u1,u2):uimeasurable,0≤ui≤1,t∈[0,tf],i=1,2}.


### 2.1. Existence of an Optimal Control Pair

The existence of the optimal control pair can be obtained using a result by Fleming and Rishel [9] and by Lukes [10].Theorem 1.
Consider the control problem with system equations (3). There exists an optimal control u*(t)=(u1*,u2*)∈𝒰 such that
(7)J(u1*,u2*)=max(u1,u2)∈𝒰J(u1,u2).Proof.
To use an existence result, Theorem III.4.1 from [9], we must check the following properties.(1)
The set of controls and corresponding state variables is nonempty.(2)
The control𝒰 set is convex and closed.(3)
The RHS of the state system is bounded by a linear function in the state and control variables.(4)
The integrand of the objective functional is concave on𝒰.(5)
There exist constantsc1, c2>0, and β>1 such that the integrand L(x;u1;u2) of the objective functional satisfies
(8)L(x,u1,u2)≤c2-c1(|u1|2+|u2|2)β/2.In order to verify these conditions, we use a result by Lukes [10, Theorem 9.2.1, page 182] to give the existence of solutions of ODEs (3) with bounded coefficients, which gives condition 1. We note that the solutions are bounded. Our control set satisfies condition 2. Since our state system is bilinear in u1 and u2, the RHS of (3) satisfies condition 3, using the boundedness of the solutions. Note that the integrand of our objective functional is concave. Also, we have the last condition needed
(9)x(t)-(A12u12(t)+A22u22(t))≤c2-c1(|u1|2+|u2|2)β/2,
where c2 depends on the upper bound on x and c1>0 since A1, A2>0. We conclude that there exists an optimal control pair.


### 2.2. Characterization of the Optimal Control

Pontryagin’s maximum principle [11] provides the necessary conditions for an optimal control problem. This principle converts (3), (4), and (5) into a problem of maximizing a Hamiltonian, H, pointwise with respect to u1 and u2:
(10)H(x,y,v,u,λi,t)=L(x,u,t)+λ1dx(t)dt+λ2dy(t)dt+λ3dv(t)dt,
where
(11)L(x,u,t)=x(t)-(A12u12(t)+A22u22(t)),λi, i=1,2,3, are the adjoint functions to be determined suitably.Theorem 2.
Given optimal controlsu1*, u2* and solutions x*(t), y*(t), and v*(t) of the corresponding state system (3), there exist adjoint variables λ1, λ2, and λ3 satisfying
(12)dλ1(t)dt=-1-λ1(t)(r(1-TK)-rxK-(1-u1)βvyT2)-λ2(t)(1-u1)βvyT2,dλ2(t)dt=-λ1(t)(-rxK+(1-u1)βvxT2)+λ2(t)(a+(1-u1)βvxT2)-λ3(t)γ(1-u2),dλ3(t)dt=(1-u1)βvxT(λ1(t)-λ2(t))+λ3(t)μwithtransversalityconditionsλi(tend)=0,i=1,2,3.
Furthermore, the optimal control u*(t) is given by
(13)u1*(t)=min(max((λ1(t)-λ2(t))βv*(t)x*(t)TA1,0),1),u2*(t)=min(max(-λ3(t)γy*(t)A2,0),1).Proof.
We use the Hamiltonian (10) in order to determine the adjoint equations and the transversality conditions. By putting x(t)=x*(t),y(t)=y*(t),  and  v(t)=v*(t) and differentiating the Hamiltonian with respect to x, y, and v, we obtain
(14)dλ1(t)dt=-∂H∂x=-1-λ1(t)(r(1-TK)-rxK-(1-u1)βvyT2)-λ2(t)(1-u1)βvyT2,dλ2(t)dt=-∂H∂y=-λ1(t)(-rxK+(1-u1)βvxT2)+λ2(t)(a+(1-u1)βvxT2)-λ3(t)γ(1-u2),dλ3(t)dt=-∂H∂v=(1-u1)βvxT(λ1(t)-λ2(t))+λ3(t)μ,
and by using the optimality conditions, we find
(15)∂H∂u1=-A1u1+λ1(t)βv*(t)x*(t)T-λ2(t)βv*(t)x*(t)T=0,atu1=u1*(t)whichgivesu1*(t)=(λ1(t)-λ2(t))βv*(t)x*(t)TA1,∂H∂u2=-A2u2-λ3(t)γy*(t)=0,atu2=u2*(t)whichgivesu2*(t)=-λ3(t)γy*(t)A2(t).
By analyzing the bounds on the controls and the representation in the interior of the set, we obtain u1*(t) and u2*(t) in the form of (13).


### 2.3. The Optimality System

The optimality system consists of the state system coupled with the adjoint system with the initial conditions, the transversality conditions, and the characterization of the optimal control.Consider(16)dx*dt=rx*(1-T*K)-(1-u1*)βv*x*T*,dy*dt=(1-u1*)βv*x*T*-ay*,dv*dt=(1-u2*)γy*-μv*,dλ1dt=-1-λ1(r(1-T*K)-rx*K-(1-u1*)βv*y*T*2)-λ2(1-u1*)βv*y*T*2,dλ2dt=-λ1(-rx*K+(1-u1*)βv*x*T*2)+λ2(a+(1-u1*)βv*x*T*2)-λ3γ(1-u2*),dλ3dt=(1-u1*)βv*x*T*(λ1-λ2)+λ3μ,withλ1(tend)=0,λ2(tend)=0,λ3(tend)=0,x(0)=x0,y(0)=y0,v(0)=v0.


## 2.1. Existence of an Optimal Control Pair

The existence of the optimal control pair can be obtained using a result by Fleming and Rishel [9] and by Lukes [10].Theorem 1.
Consider the control problem with system equations (3). There exists an optimal control u*(t)=(u1*,u2*)∈𝒰 such that
(7)J(u1*,u2*)=max(u1,u2)∈𝒰J(u1,u2).Proof.
To use an existence result, Theorem III.4.1 from [9], we must check the following properties.(1)
The set of controls and corresponding state variables is nonempty.(2)
The control𝒰 set is convex and closed.(3)
The RHS of the state system is bounded by a linear function in the state and control variables.(4)
The integrand of the objective functional is concave on𝒰.(5)
There exist constantsc1, c2>0, and β>1 such that the integrand L(x;u1;u2) of the objective functional satisfies
(8)L(x,u1,u2)≤c2-c1(|u1|2+|u2|2)β/2.In order to verify these conditions, we use a result by Lukes [10, Theorem 9.2.1, page 182] to give the existence of solutions of ODEs (3) with bounded coefficients, which gives condition 1. We note that the solutions are bounded. Our control set satisfies condition 2. Since our state system is bilinear in u1 and u2, the RHS of (3) satisfies condition 3, using the boundedness of the solutions. Note that the integrand of our objective functional is concave. Also, we have the last condition needed
(9)x(t)-(A12u12(t)+A22u22(t))≤c2-c1(|u1|2+|u2|2)β/2,
where c2 depends on the upper bound on x and c1>0 since A1, A2>0. We conclude that there exists an optimal control pair.


## 2.2. Characterization of the Optimal Control

Pontryagin’s maximum principle [11] provides the necessary conditions for an optimal control problem. This principle converts (3), (4), and (5) into a problem of maximizing a Hamiltonian, H, pointwise with respect to u1 and u2:
(10)H(x,y,v,u,λi,t)=L(x,u,t)+λ1dx(t)dt+λ2dy(t)dt+λ3dv(t)dt,
where
(11)L(x,u,t)=x(t)-(A12u12(t)+A22u22(t)),λi, i=1,2,3, are the adjoint functions to be determined suitably.Theorem 2.
Given optimal controlsu1*, u2* and solutions x*(t), y*(t), and v*(t) of the corresponding state system (3), there exist adjoint variables λ1, λ2, and λ3 satisfying
(12)dλ1(t)dt=-1-λ1(t)(r(1-TK)-rxK-(1-u1)βvyT2)-λ2(t)(1-u1)βvyT2,dλ2(t)dt=-λ1(t)(-rxK+(1-u1)βvxT2)+λ2(t)(a+(1-u1)βvxT2)-λ3(t)γ(1-u2),dλ3(t)dt=(1-u1)βvxT(λ1(t)-λ2(t))+λ3(t)μwithtransversalityconditionsλi(tend)=0,i=1,2,3.
Furthermore, the optimal control u*(t) is given by
(13)u1*(t)=min(max((λ1(t)-λ2(t))βv*(t)x*(t)TA1,0),1),u2*(t)=min(max(-λ3(t)γy*(t)A2,0),1).Proof.
We use the Hamiltonian (10) in order to determine the adjoint equations and the transversality conditions. By putting x(t)=x*(t),y(t)=y*(t),  and  v(t)=v*(t) and differentiating the Hamiltonian with respect to x, y, and v, we obtain
(14)dλ1(t)dt=-∂H∂x=-1-λ1(t)(r(1-TK)-rxK-(1-u1)βvyT2)-λ2(t)(1-u1)βvyT2,dλ2(t)dt=-∂H∂y=-λ1(t)(-rxK+(1-u1)βvxT2)+λ2(t)(a+(1-u1)βvxT2)-λ3(t)γ(1-u2),dλ3(t)dt=-∂H∂v=(1-u1)βvxT(λ1(t)-λ2(t))+λ3(t)μ,
and by using the optimality conditions, we find
(15)∂H∂u1=-A1u1+λ1(t)βv*(t)x*(t)T-λ2(t)βv*(t)x*(t)T=0,atu1=u1*(t)whichgivesu1*(t)=(λ1(t)-λ2(t))βv*(t)x*(t)TA1,∂H∂u2=-A2u2-λ3(t)γy*(t)=0,atu2=u2*(t)whichgivesu2*(t)=-λ3(t)γy*(t)A2(t).
By analyzing the bounds on the controls and the representation in the interior of the set, we obtain u1*(t) and u2*(t) in the form of (13).


## 2.3. The Optimality System

The optimality system consists of the state system coupled with the adjoint system with the initial conditions, the transversality conditions, and the characterization of the optimal control.Consider(16)dx*dt=rx*(1-T*K)-(1-u1*)βv*x*T*,dy*dt=(1-u1*)βv*x*T*-ay*,dv*dt=(1-u2*)γy*-μv*,dλ1dt=-1-λ1(r(1-T*K)-rx*K-(1-u1*)βv*y*T*2)-λ2(1-u1*)βv*y*T*2,dλ2dt=-λ1(-rx*K+(1-u1*)βv*x*T*2)+λ2(a+(1-u1*)βv*x*T*2)-λ3γ(1-u2*),dλ3dt=(1-u1*)βv*x*T*(λ1-λ2)+λ3μ,withλ1(tend)=0,λ2(tend)=0,λ3(tend)=0,x(0)=x0,y(0)=y0,v(0)=v0.


## 3. Numerical Simulations

### 3.1. The Improved GSS1 Method

The resolution of the optimality system is created by improving the Gauss-Seidel-like implicit finite-difference method developed by Gumel et al. [12], presented in [13–15], and denoted by the GSS1 method. It consists of discretizing the interval [0,tend] by the points tk=kl(k=0,1,…,n), where l is the time step. Next, we define the state and adjoint variables x(t), y(t), v(t), λ1(t), λ2(t), and λ3(t) and the controls u1(t) and u2(t) in terms of nodal points xk, yk, vk, λ1k, λ2k, λ3k, u1k, and u2k with x0, y0, v0, λ10, λ20, λ30, u10 and u20 as the state and adjoint variables and the controls at initial time t0. xn, yn, vn, λ1n, λ2n, λ3n, u1n, and u2n are the state and adjoint variables and the controls at final time tend.As it is well known, the approximation of the time derivative by its first-order forward difference is given, for the first state variablex(t), by
(17)dx(t)dt=liml→0x(t+l)-x(t)l.
We use the GSS1 method to adapt it to our case as follows: we examine the variables one by one by blocking all other values to the most recently calculated
(18)xk+1-xkl=rxk+1(1-xk+1+ykK)-βvkxk+1xk+1+yk,yk+1-ykl=(1-u1k)βvkxk+1xk+1+yk+1-ayk+1,vk+1-vkl=(1-u2k)γyk+1-μvk+1.
By applying an analogous technology, we approximate the time derivative of the adjoint variables by their first-order backward difference, and we use the appropriated scheme as follows:
(19)λ1n-k-λ1n-k-1l=-1-λ1n-k-1×(r(1-Tk+1K)-rxk+1K-(1-u1k)βvk+1yk+1(Tk+1)2)-λ2n-k(1-u1k)βvk+1yk+1(Tk+1)2,λ2n-k-λ2n-k-1l=λ1n-k-1(rxk+1K-(1-u1k)βvk+1xk+1(Tk+1)2)+λ2n-k-1((1-u1k)βvk+1xk+1(Tk+1)2+a)-λ3n-kγ(1-u2k),λ3n-k-λ3n-k-1l=(1-u1k)βxk+1Tk+1(λ1n-k-1-λ2n-k-1)+λ3n-k-1μ.


### 3.2. Numerical Results

By making some parameter value choices, computer simulation experiments are done to verify the effectiveness of our new model by comparing the disease progression before and after introducing the two optimal control variablesu1*(t) and u2*(t). We use the following parameters and initial values:
(20)x0=140,y0=145,v0=18,μ=0.73,γ=250,β=0.0014,r=0.01,a=0.0693,K=340,A1=900,A2=1000.
In Figures 1, 2, 3, 4, and 5, the control individuals are marked by solid blue curves, while the individual without control is marked by dashed black curves.Figure 1
Shows a significant difference in the mass of healthy hepatocytes with and without control from the very beginning of treatment (in the first three days), and after that, it begins go to the stable state.Figure 2
We see that the mass of infected hepatocytes, with control, decreases rapidly from the beginning days of treatment.Figure 3
We see that the mass of free virions decreases rapidly from the first week of the treatment.Figure 4
Represents the controlsu1* and u2*.Figure 5
Represents the controlsu1* and u2*.


## 3.1. The Improved GSS1 Method

The resolution of the optimality system is created by improving the Gauss-Seidel-like implicit finite-difference method developed by Gumel et al. [12], presented in [13–15], and denoted by the GSS1 method. It consists of discretizing the interval [0,tend] by the points tk=kl(k=0,1,…,n), where l is the time step. Next, we define the state and adjoint variables x(t), y(t), v(t), λ1(t), λ2(t), and λ3(t) and the controls u1(t) and u2(t) in terms of nodal points xk, yk, vk, λ1k, λ2k, λ3k, u1k, and u2k with x0, y0, v0, λ10, λ20, λ30, u10 and u20 as the state and adjoint variables and the controls at initial time t0. xn, yn, vn, λ1n, λ2n, λ3n, u1n, and u2n are the state and adjoint variables and the controls at final time tend.As it is well known, the approximation of the time derivative by its first-order forward difference is given, for the first state variablex(t), by
(17)dx(t)dt=liml→0x(t+l)-x(t)l.
We use the GSS1 method to adapt it to our case as follows: we examine the variables one by one by blocking all other values to the most recently calculated
(18)xk+1-xkl=rxk+1(1-xk+1+ykK)-βvkxk+1xk+1+yk,yk+1-ykl=(1-u1k)βvkxk+1xk+1+yk+1-ayk+1,vk+1-vkl=(1-u2k)γyk+1-μvk+1.
By applying an analogous technology, we approximate the time derivative of the adjoint variables by their first-order backward difference, and we use the appropriated scheme as follows:
(19)λ1n-k-λ1n-k-1l=-1-λ1n-k-1×(r(1-Tk+1K)-rxk+1K-(1-u1k)βvk+1yk+1(Tk+1)2)-λ2n-k(1-u1k)βvk+1yk+1(Tk+1)2,λ2n-k-λ2n-k-1l=λ1n-k-1(rxk+1K-(1-u1k)βvk+1xk+1(Tk+1)2)+λ2n-k-1((1-u1k)βvk+1xk+1(Tk+1)2+a)-λ3n-kγ(1-u2k),λ3n-k-λ3n-k-1l=(1-u1k)βxk+1Tk+1(λ1n-k-1-λ2n-k-1)+λ3n-k-1μ.


## 3.2. Numerical Results

By making some parameter value choices, computer simulation experiments are done to verify the effectiveness of our new model by comparing the disease progression before and after introducing the two optimal control variablesu1*(t) and u2*(t). We use the following parameters and initial values:
(20)x0=140,y0=145,v0=18,μ=0.73,γ=250,β=0.0014,r=0.01,a=0.0693,K=340,A1=900,A2=1000.
In Figures 1, 2, 3, 4, and 5, the control individuals are marked by solid blue curves, while the individual without control is marked by dashed black curves.Figure 1
Shows a significant difference in the mass of healthy hepatocytes with and without control from the very beginning of treatment (in the first three days), and after that, it begins go to the stable state.Figure 2
We see that the mass of infected hepatocytes, with control, decreases rapidly from the beginning days of treatment.Figure 3
We see that the mass of free virions decreases rapidly from the first week of the treatment.Figure 4
Represents the controlsu1* and u2*.Figure 5
Represents the controlsu1* and u2*.


## 4. Conclusion

In this paper, we have proposed an HBV model with control variables.Our goal was to maximize the mass of healthy hepatocytes and to minimize the cost of therapy, and effectively we can clearly notice the benefit of putting a control on the decreasing kinetic of infected hepatocytes. We believe that the analysis presented in this paper, combined with pharmacokinetics studies, could play an important role in developing improved HBV treatment regimen.



---

*Source: 912835-2013-06-26.xml*