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

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.

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−𝐾−𝛽𝑣(𝑡)𝑥(𝑡)𝑇(𝑡),𝑥(0)=𝑥₀≥0,𝑑𝑦=(𝑡)𝑥(𝑡)𝑑𝑡𝛽𝑣𝑇(𝑡)−𝑎𝑦(𝑡),𝑦(0)=𝑦₀≥0,𝑑𝑣(𝑡)(𝑡)(0)𝑑𝑡=𝛾𝑦−𝜇𝑣,𝑣=𝑣₀≥0, (1) with 𝑇(𝑡)=𝑥(𝑡)+𝑦(𝑡), (2) where 𝑥(𝑡) is the mass of healthy hepatocytes, 𝑦(𝑡) is the mass of infected hepatocytes, and 𝑣(𝑡) is the mass of free virions. Healthy hepatocytes grow at a rate that depends on the homeostatic liver size, 𝐾, at a maximum per capita proliferation rate 𝑟. Hepatocytes become infected at maximum rate 𝛽 and die at rate 𝑎. 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 control 𝑢(𝑡)=(𝑢₁(𝑡),𝑢₂(𝑡)) to the previously mentioned model simulating the antiviral therapy. Hence, the previous model becomes 𝑑𝑥(𝑡)𝑇(𝑡)𝑑𝑡=𝑟𝑥1−𝐾−1−𝑢₁(𝑡)𝛽𝑣(𝑡)𝑥(𝑡)𝑇(𝑡),𝑥(0)=𝑥₀≥0,𝑑𝑦=𝑑𝑡1−𝑢₁(𝑡)𝛽𝑣(𝑡)𝑥(𝑡)𝑇(𝑡)(𝑡),𝑦(0)−𝑎𝑦=𝑦₀≥0,𝑑𝑣=𝑑𝑡1−𝑢₂(𝑡)𝛾𝑦(𝑡)−𝜇𝑣(𝑡),𝑣(0)=𝑣₀≥0. (3) The control functions 𝑢₁(𝑡) and 𝑢₂(𝑡) are bounded Lebesgue integrable functions. The control 𝑢₁(𝑡) represents the efficiency of drug therapy in blocking new infection, and the control 𝑢₂(𝑡) represents the efficiency of drug therapy in inhibiting viral production.

If 𝑢₂=1, the inhibition is 100% effective, whereas if 𝑢₂=0, there is no inhibition.

The problem is to maximize the objective functional 𝐽𝑢₁(𝑡),𝑢₂(𝑡)=^𝑡end₀𝑥(𝑡)−𝐴₁2𝑢²₁(𝑡)+𝐴₂2𝑢²₂(𝑡)𝑑𝑡, (4) where the parameters 𝐴₁≥0 and 𝐴₂≥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 (𝑢^∗₁,𝑢^∗₂) such that 𝐽𝑢^∗₁,𝑢^∗₂𝐽𝑢=max¹,𝑢₂𝑢₁,𝑢₂∈𝒰, (5) where 𝒰 is the control set defined by 𝒰=𝑢=𝑢₁,𝑢₂∶𝑢_𝑖measurable,0≤𝑢_𝑖≤1,𝑡∈0,𝑡_𝑓,𝑖=1,2. (6)

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 𝑢^∗(𝑡)=(𝑢^∗₁,𝑢^∗₂)∈𝒰 such that 𝐽𝑢^∗₁,𝑢^∗₂=max_{(𝑢1,𝑢2)∈𝒰}𝐽𝑢₁,𝑢₂. (7)

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 constants 𝑐₁, 𝑐₂>0, and 𝛽>1 such that the integrand 𝐿(𝑥;𝑢₁;𝑢₂) of the objective functional satisfies 𝐿𝑥,𝑢₁,𝑢₂≤𝑐₂−𝑐₁𝑢₁²+𝑢₂^2𝛽/2. (8)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 𝑢₁ and 𝑢₂, 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 𝑥(𝑡)−𝐴₁2𝑢²₁(𝑡)+𝐴₂2𝑢²₂(𝑡)≤𝑐₂−𝑐₁𝑢₁²+𝑢₂^2𝛽/2, (9) where 𝑐₂ depends on the upper bound on 𝑥 and 𝑐₁>0 since 𝐴₁, 𝐴₂>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, 𝐻, pointwise with respect to 𝑢₁ and 𝑢₂: 𝐻𝑥,𝑦,𝑣,𝑢,𝜆_𝑖,𝑡=𝐿(𝑥,𝑢,𝑡)+𝜆₁𝑑𝑥(𝑡)𝑑𝑡+𝜆₂𝑑𝑦(𝑡)𝑑𝑡+𝜆₃𝑑𝑣(𝑡),𝑑𝑡 (10) where 𝐿()(𝑡)−𝐴𝑥,𝑢,𝑡=𝑥¹2𝑢²₁(𝑡)+𝐴₂2𝑢²₂(𝑡), (11)𝜆_𝑖, 𝑖=1,2,3, are the adjoint functions to be determined suitably.

Theorem 2. Given optimal controls 𝑢^∗₁, 𝑢^∗₂ and solutions 𝑥^∗(𝑡), 𝑦^∗(𝑡), and 𝑣^∗(𝑡) of the corresponding state system (3), there exist adjoint variables 𝜆₁, 𝜆₂, and 𝜆₃ satisfying 𝑑𝜆₁(𝑡)𝑑𝑡=−1−𝜆₁(𝑡)𝑟𝑇1−𝐾−𝑟𝑥𝐾−1−𝑢₁𝑦𝛽𝑣𝑇²−𝜆₂(𝑡)1−𝑢₁𝑦𝛽𝑣𝑇²,𝑑𝜆₂(𝑡)𝑑𝑡=−𝜆₁(𝑡)−𝑟𝑥𝐾+1−𝑢₁𝑥𝛽𝑣𝑇²+𝜆₂(𝑡)𝑎+1−𝑢₁𝑥𝛽𝑣𝑇²−𝜆₃(𝑡)𝛾1−𝑢₂,𝑑𝜆₃(𝑡)=𝑑𝑡1−𝑢₁𝑥𝛽𝑣𝑇𝜆₁(𝑡)−𝜆₂(𝑡)+𝜆₃(𝑡)𝜇𝑤𝑖𝑡ℎ𝑡𝑟𝑎𝑛𝑠𝑣𝑒𝑟𝑠𝑎𝑙𝑖𝑡𝑦𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛𝑠𝜆_𝑖𝑡_𝑒𝑛𝑑=0,𝑖=1,2,3. (12) Furthermore, the optimal control 𝑢^∗(𝑡) is given by 𝑢^∗₁(𝑡)𝜆=minmax¹(𝑡)−𝜆₂(𝑡)𝛽𝑣^∗(𝑡)𝑥^∗(𝑡)𝑇𝐴₁,𝑢,0,1_∗2(𝑡)−𝜆=minmax³(𝑡)𝛾𝑦^∗(𝑡)𝐴₂,0,1. (13)

Proof. We use the Hamiltonian (10) in order to determine the adjoint equations and the transversality conditions. By putting 𝑥(𝑡)=𝑥^∗(𝑡),𝑦(𝑡)=𝑦^∗(𝑡),  and  𝑣(𝑡)=𝑣^∗(𝑡) and differentiating the Hamiltonian with respect to 𝑥, 𝑦, and 𝑣, we obtain 𝑑𝜆₁(𝑡)𝑑𝑡=−𝜕𝐻𝜕𝑥=−1−𝜆₁(𝑡)𝑟𝑇1−𝐾−𝑟𝑥𝐾−1−𝑢₁𝑦𝛽𝑣𝑇²−𝜆₂(𝑡)1−𝑢₁𝑦𝛽𝑣𝑇²,𝑑𝜆₂(𝑡)𝑑𝑡=−𝜕𝐻𝜕𝑦=−𝜆₁(𝑡)−𝑟𝑥𝐾+1−𝑢₁𝑥𝛽𝑣𝑇²+𝜆₂(𝑡)𝑎+1−𝑢₁𝑥𝛽𝑣𝑇²−𝜆₃(𝑡)𝛾1−𝑢₂,𝑑𝜆₃(𝑡)𝑑𝑡=−𝜕𝐻=𝜕𝑣1−𝑢₁𝑥𝛽𝑣𝑇𝜆₁(𝑡)−𝜆₂(𝑡)+𝜆₃(𝑡)𝜇, (14) and by using the optimality conditions, we find 𝜕𝐻𝜕𝑢₁=−𝐴₁𝑢₁+𝜆₁(𝑡)𝛽𝑣^∗(𝑡)𝑥^∗(𝑡)𝑇−𝜆₂(𝑡)𝛽𝑣^∗(𝑡)𝑥^∗(𝑡)𝑇=0,at𝑢₁=𝑢^∗₁(𝑡)𝑢whichgives_∗1(𝑡)=𝜆₁(𝑡)−𝜆₂(𝑡)𝛽𝑣^∗(𝑡)𝑥^∗(𝑡)𝑇𝐴₁,𝜕𝐻𝜕𝑢₂=−𝐴₂𝑢₂−𝜆₃(𝑡)𝛾𝑦^∗(𝑡)=0,at𝑢₂=𝑢^∗₂(𝑡)𝑢whichgives_∗2(𝑡)𝜆=−³(𝑡)𝛾𝑦^∗(𝑡)𝐴₂(𝑡). (15) By analyzing the bounds on the controls and the representation in the interior of the set, we obtain 𝑢^∗₁(𝑡) and 𝑢^∗₂(𝑡) 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 𝑑𝑥^∗𝑑𝑡=𝑟𝑥^∗𝑇1−^∗𝐾−1−𝑢^∗₁𝛽𝑣^∗𝑥^∗𝑇^∗,𝑑𝑦^∗=𝑑𝑡1−𝑢^∗₁𝛽𝑣^∗𝑥^∗𝑇^∗−𝑎𝑦^∗,𝑑𝑣^∗=𝑑𝑡1−𝑢^∗₂𝛾𝑦^∗−𝜇𝑣^∗,𝑑𝜆₁𝑑𝑡=−1−𝜆₁𝑟𝑇1−^∗𝐾−𝑟𝑥^∗𝐾−1−𝑢^∗₁𝛽𝑣^∗𝑦^∗𝑇^∗2−𝜆₂1−𝑢^∗₁𝛽𝑣^∗𝑦^∗𝑇^∗2,𝑑𝜆₂𝑑𝑡=−𝜆₁−𝑟𝑥^∗𝐾+1−𝑢^∗₁𝛽𝑣^∗𝑥^∗𝑇^∗2+𝜆₂𝑎+1−𝑢^∗₁𝛽𝑣^∗𝑥^∗𝑇^∗2−𝜆₃𝛾1−𝑢^∗₂,𝑑𝜆₃=𝑑𝑡1−𝑢^∗₁𝛽𝑣^∗𝑥^∗𝑇^∗𝜆₁−𝜆₂+𝜆₃𝜇,with𝜆₁𝑡_end=0,𝜆₂𝑡_end=0,𝜆₃𝑡_end𝑥(0)=0,=𝑥₀,𝑦(0)=𝑦₀,𝑣(0)=𝑣₀. (16)

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,𝑡_end] by the points 𝑡_𝑘=𝑘𝑙(𝑘=0,1,…,𝑛), where 𝑙 is the time step. Next, we define the state and adjoint variables 𝑥(𝑡), 𝑦(𝑡), 𝑣(𝑡), 𝜆₁(𝑡), 𝜆₂(𝑡), and 𝜆₃(𝑡) and the controls 𝑢₁(𝑡) and 𝑢₂(𝑡) in terms of nodal points 𝑥^𝑘, 𝑦^𝑘, 𝑣^𝑘, 𝜆^𝑘₁, 𝜆^𝑘₂, 𝜆^𝑘₃, 𝑢^𝑘₁, and 𝑢^𝑘₂ with 𝑥⁰, 𝑦⁰, 𝑣⁰, 𝜆⁰₁, 𝜆⁰₂, 𝜆⁰₃, 𝑢⁰₁ and 𝑢⁰₂ as the state and adjoint variables and the controls at initial time 𝑡₀. 𝑥^𝑛, 𝑦^𝑛, 𝑣^𝑛, 𝜆^𝑛₁, 𝜆^𝑛₂, 𝜆^𝑛₃, 𝑢^𝑛₁, and 𝑢^𝑛₂ are the state and adjoint variables and the controls at final time 𝑡_end.

As it is well known, the approximation of the time derivative by its first-order forward difference is given, for the first state variable 𝑥(𝑡), by 𝑑𝑥(𝑡)𝑑𝑡=lim_𝑙→0𝑥(𝑡+𝑙)−𝑥(𝑡)𝑙. (17) 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 𝑥^𝑘+1−𝑥^𝑘𝑙=𝑟𝑥^𝑘+1𝑥1−^𝑘+1+𝑦^𝑘𝐾−𝛽𝑣^𝑘𝑥^𝑘+1𝑥^𝑘+1+𝑦^𝑘,𝑦^𝑘+1−𝑦^𝑘𝑙=1−𝑢^𝑘₁𝛽𝑣^𝑘𝑥^𝑘+1𝑥^𝑘+1+𝑦^𝑘+1−𝑎𝑦^𝑘+1,𝑣^𝑘+1−𝑣^𝑘𝑙=1−𝑢^𝑘₂𝛾𝑦^𝑘+1−𝜇𝑣^𝑘+1. (18) 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: 𝜆₁^𝑛−𝑘−𝜆₁^{𝑛−𝑘−1}𝑙=−1−𝜆₁^{𝑛−𝑘−1}×𝑟𝑇1−^𝑘+1𝐾−𝑟𝑥^𝑘+1𝐾−1−𝑢^𝑘₁𝛽𝑣^𝑘+1𝑦^𝑘+1𝑇^𝑘+12−𝜆₂^𝑛−𝑘1−𝑢^𝑘₁𝛽𝑣^𝑘+1𝑦^𝑘+1𝑇^𝑘+12,𝜆₂^𝑛−𝑘−𝜆₂^{𝑛−𝑘−1}𝑙=𝜆₁^{𝑛−𝑘−1}𝑟𝑥^𝑘+1𝐾−1−𝑢^𝑘₁𝛽𝑣^𝑘+1𝑥^𝑘+1𝑇^𝑘+12+𝜆₂^{𝑛−𝑘−1}1−𝑢^𝑘₁𝛽𝑣^𝑘+1𝑥^𝑘+1𝑇^𝑘+12+𝑎−𝜆₃^𝑛−𝑘𝛾1−𝑢^𝑘₂,𝜆₃^𝑛−𝑘−𝜆₃^{𝑛−𝑘−1}𝑙=1−𝑢^𝑘₁𝛽𝑥^𝑘+1𝑇^𝑘+1𝜆₁^{𝑛−𝑘−1}−𝜆₂^{𝑛−𝑘−1}+𝜆₃^{𝑛−𝑘−1}𝜇. (19)

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 variables 𝑢^∗₁(𝑡) and 𝑢^∗₂(𝑡). We use the following parameters and initial values: 𝑥₀=140,𝑦₀=145,𝑣₀𝐴=18,𝜇=0.73,𝛾=250,𝛽=0.0014,𝑟=0.01,𝑎=0.0693,𝐾=340,₁=900,𝐴₂=1000. (20) 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 controls 𝑢∗1 and 𝑢∗2.

Figure 5 

Represents the controls 𝑢∗1 and 𝑢∗2.

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.

Acknowledgments

The authors would like to thank the anonymous referees for their valuable comments. Some of this work was funded by the network of systems theory affiliated to Hassan II Academy of Science and Technology.