An Efficient Therapy Strategy under a Novel HIV Model

Abstract

By incorporating the chemotherapy into a previous model describing the interaction of the immune system with the human immunodeficiency virus (HIV), this paper proposes a novel HIV virus spread model with control variables. Our goal is to maximize the number of healthy cells and, meanwhile, to minimize the cost of chemotherapy. In this context, the existence of an optimal control is proved. Experimental results show that, under this model, the spread of HIV virus can be controlled effectively.

1. Introduction

Numerous studies have been devoted to the description and understanding of the spread of infectious diseases (especially, the acquired immunodeficiency syndrome (AIDS)) [1–18]. Mathematical modeling of the human immunodeficiency virus (HIV) viral dynamics has offered many insights into the pathogenesis and treatment of HIV [1, 2, 4–10, 12–16, 18]. Consequently, many mathematical models have been developed to depict the relationships among HIV, etiological agent for AIDS and CD4⁺T lymphoblasts, which are the targets for the virus [13]. Some of these models investigate how to avoid an excessive use of drugs because it might be toxic to human body and, hence, cause damages [1, 4–6, 8–11, 14, 15, 17, 18].

Recently, Sedaghat et al. [13] proposed a model, which describes the law governing the transition of two populations of target cells, the T cells (the abbreviation of the CD4⁺T lymphoblasts) and the M cells (say, macrophages, T cells in a lower state of activation, or another cell type), in the effect of free virus (see Figure 1). The T cells produce most of the plasma virus and are responsible for the first-phase decay, while the M cells are responsible for the second-phase decay. T cells are classified into three categories: 𝑇_𝑈 cells (uninfected T cells), 𝑇₁ cells (early-stage infected T cells), and 𝑇₂ cells (late-stage infected T cells). Let 𝑇_𝑈, 𝑇₁ and 𝑇₂ denote the numbers of 𝑇_𝑈 cells, 𝑇₁ cells, and 𝑇₂ cells, respectively. Likewise, M cells are classified into three categories: 𝑀_𝑢 cells (uninfected M cells), 𝑀₁ cells (early-stage infected M cells), and 𝑀₂ cells (late-stage infected M cells). Let 𝑀_𝑢, 𝑀₁, and 𝑀₂ denote the numbers of 𝑀_𝑈 cells, 𝑀₁ cells and 𝑀₂ cells, respectively. Besides, let 𝑉 denote the number of free viruses. Sedaghat et al. [13] made the following reasonable assumptions.(A₁)𝑇_𝑈 cells are produced with constant rate 𝜃_𝑇. 𝑀_𝑈 cells are produced with constant rate 𝜃_𝑀. (A₂)𝑇_𝑈 cells become 𝑇₁ cells with constant rate 𝛽_𝑇. 𝑀_𝑈 cells become 𝑀₁ cells with constant rate 𝛽_𝑀. (A₃)𝑇₁ cells become 𝑇₂ cells with constant rate 𝑘_𝑇. 𝑀₁ cells become 𝑀₂ cells with constant rate 𝑘_𝑀. (A₄)These cells die with constant rates 𝛿_𝑇𝑈, 𝛿_𝑇1, 𝛿_𝑇2, 𝛿_𝑀𝑈, 𝛿_𝑀1, and 𝛿_𝑀2 respectively. (A₅)Free viruses (𝑉) are cleared at a rate 𝑐, produced by 𝑇₂ cells with a burst size of 𝑁_𝑇, and produced by 𝑀₂ cells with a burst size of 𝑁_𝑀, respectively.

Figure 1 

The HIV model.

Under these assumptions, Sedaghat et al. [13] deduced the following system of ordinary differential equations:𝑑𝑇_𝑈𝑑𝑡=𝜃_𝑇−𝛿_𝑇𝑈𝑇_𝑈−𝛽_𝑇𝑇_𝑈𝑉,𝑑𝑇₁𝑑𝑡=𝛽_𝑇𝑇_𝑈𝑉−𝛿_𝑇1+𝑘_𝑇𝑇₁,𝑑𝑇₂𝑑𝑡=𝑘_𝑇𝑇₁−𝛿_𝑇2𝑇₂,𝑑𝑀_𝑈𝑑𝑡=𝜃_𝑀−𝛿_𝑀𝑈𝑀_𝑈−𝛽_𝑀𝑀_𝑈𝑉,𝑑𝑀₁𝑑𝑡=𝛽_𝑀𝑀_𝑈𝑉−𝛿_𝑀1+𝑘_𝑀𝑀₁,𝑑𝑀₂𝑑𝑡=𝑘_𝑀𝑀₁−𝛿_𝑀2𝑀₂,𝑑𝑉𝑑𝑡=𝑁_𝑇𝑇₂+𝑁_𝑀𝑀₂−𝑐𝑉. (1.1) For a highly simplified version of this system, Sedaghat et al. [13] derived its analytic solution.

It is well known [5, 6, 8–11, 13, 15, 17] that there are mainly two categories of anti-HIV drugs: the reverse transcriptase inhibitors (RTIs), which prevent new HIV infection by disrupting the conversion of viral RNA into DNA inside of T cells, and the protease inhibitors (PIs), which reduce the number of virus particles produced by actively-infected T cells.

In consideration of this, this paper introduces a novel HIV model by incorporating the drug dosage into the above-mentioned model. Our goal is to maximize the number of healthy cells and, meanwhile, to minimize the cost of chemotherapy. In this context, the existence of an optimal control strategy is proved. Experimental results show that, under this model, the spread of HIV virus can be controlled effectively.

2. Presentation of a New Model

For our purpose, let us introduce the following notations (see Figure 2): 𝑢₁(𝑡):the dosage of RTI at time 𝑡, which is assumed to take values in the interval [0,1];𝑢₂(𝑡):the dosage of PI at time 𝑡, which is assumed to take values in [0,1];𝛾₁:the capability of preventing 𝑇_𝑈 cells from becoming 𝑇₁ cells with per unit dosage of RTI;𝛾₂:the capability of preventing 𝑀_𝑈 cells from becoming 𝑀₁ cells with per unit dosage of RTI;𝛼₁:the capability of preventing 𝑇₂ cells from producing viruses with per unit dosage of PI;𝛼₂:the capability of preventing 𝑀₂ cells from producing viruses with per unit dosage of PI.

Figure 2 

The HIV model with therapy strategy.

Next, let us consider the following assumptions.(A₆)Due to the effect of RTIs, 𝑇_𝑈 cells become 𝑇₁ cells with rate 𝛽_𝑇(1−𝑢₁(𝑡))𝛾₁, and 𝑀_𝑈 cells become 𝑀₁ cells with rate 𝛽_𝑀(1−𝑢₁(𝑡))𝛾₂, where 𝛾₁ and 𝛾₂ are constants. (A₇)Due to the effect of PIs, Free viruses (𝑉) are produced by 𝑇₂ and 𝑀₂ cells with a burst size of 𝛼₁(1−𝑢₂(𝑡))𝑁_𝑇 and 𝛼₂(1−𝑢₂(𝑡))𝑁_𝑀, respectively, where 𝛼₁ and 𝛼₂ are constants.

Under assumptions (A₁)–(A₇), we can derive the following system of ordinary differential equations:𝑑𝑇_𝑈𝑑𝑡=𝜃_𝑇−𝛿_𝑇𝑈𝑇_𝑈−𝛽_𝑇𝑉𝑇_𝑈1−𝑢₁𝛾₁,𝑑𝑇₁𝑑𝑡=𝛽_𝑇𝑉𝑇_𝑈1−𝑢₁𝛾₁−𝛿_𝑇1+𝑘_𝑇𝑇₁,𝑑𝑇₂𝑑𝑡=𝑘_𝑇𝑇₁−𝛿_𝑇2𝑇₂,𝑑𝑀_𝑈𝑑𝑡=𝜃_𝑀−𝛿_𝑀𝑈𝑀_𝑈−𝛽_𝑀𝑉𝑀_𝑈1−𝑢₁𝛾₂,𝑑𝑀₁𝑑𝑡=𝛽_𝑀𝑉𝑀_𝑈1−𝑢₁𝛾₂−𝛿_𝑀1+𝑘_𝑀𝑀₁,𝑑𝑀₂𝑑𝑡=𝑘_𝑀𝑀₁−𝛿_𝑀2𝑀₂,𝑑𝑉𝑑𝑡=𝛼₁𝑁_𝑇𝑇₂1−𝑢₂+𝛼₂𝑁_𝑀𝑀₂1−𝑢₂−𝑐𝑉. (2.1)

Our target is to maximize the objective functional by increasing the number of healthy T and M cells and minimizing the cost based on the percentage effect of the chemotherapy given. For that purpose, we introduce the following objective functional 𝒥𝑢₁(𝑡),𝑢₂(𝑡)=^𝑡1_𝑡0𝐵₁𝑇_𝑈+𝐵₂𝑀_𝑈−𝐴₁𝑢²₁+𝐴₂𝑢²₂𝑑𝑡, (2.2) where 𝐵₁,𝐵₂ represent the benefit of per 𝑇_𝑈 cell and per 𝑀_𝑈 cell, respectively, and 𝐴₁,𝐴₂ represent the cost of per unit RTI and per unit PI, respectively. Our goal is to obtain an optimal control pair (𝑢^∗₁,𝑢^∗₂) such that𝒥𝑢^∗₁,𝑢^∗₂=max𝒥𝑢₁,𝑢₂∶𝑢₁,𝑢₂∈𝒰, (2.3) where 𝒰 is the admissible control set defined by𝒰=𝑈₁×𝑈₂,𝑈₁=𝑈₂=𝑢(𝑡)∶𝑢measurable,0≤𝑢(𝑡)≤1,𝑡∈𝑡₀,𝑡₁. (2.4)

3. Existence of an Optimal Control Pair

For our purpose, let us introduce the following four assumptions. (A₈)The set of control and corresponding state variables is nonempty.(A₉)The admissible control set 𝒰 is closed and convex.(A₁₀)All the right hand sides of equations of system (2.1) are continuous, bounded above by a sum of bounded control and state, and can be written as a linear function of 𝑢 with coefficients depending on time and state.(A₁₁)There exist positive constants 𝑐₁,𝑐₂ and 𝛽>1 such that the integrand (denoted by 𝐿(𝑦,𝑢,𝑡)) of the objective functional (2.2) is concave and satisfies the condition 𝐿(𝑦,𝑢,𝑡)≤𝑐₁−𝑐₂(𝑢²₁+𝑢²₂)^𝛽/2.

In what follows, it is always assumed that assumptions (A₁)–(A₇) hold.

Theorem 3.1. Consider system (2.1) with initial conditions, and the objective functional (2.2). There exists 𝑢^∗=(𝑢^∗₁,𝑢^∗₂) such that 𝒥𝑢^∗₁,𝑢^∗₂=max_𝑢∈𝑈𝒥𝑢₁,𝑢₂. (3.1)

Proof. It suffices to verify the assumptions (A₈)–(A₁₁) with respect to the seven ODEs of system (2.1).
Since the coefficients involved in the system are bounded, and each state variable of the system is bounded on the finite time interval, it follows by a result (see Appendix A) from [19] we can obtain the existence to the solution of the system (2.1).
The control set 𝒰=𝑈₁×𝑈₂ is obviously closed and convex, because both 𝑈₁ and 𝑈₂ are closed and convex sets.
By definition, each right hand side of the ODEs of system (2.1) is continuous and can be written as a linear function of 𝑢 with coefficients depending on time and states. The fact that all state variables 𝑇_𝑈, 𝑇₁, 𝑇₂, 𝑀_𝑈, 𝑀₁, 𝑀₂, 𝑉, and 𝒰 are bounded on [𝑡₀,𝑡₁], implies the rest of assumption (A₁₀).
It is easy to see that 𝐿(𝑦,𝑢,𝑡) is concave in 𝒰. By setting 𝑐₁=max{𝐵₁𝑇_𝑈+𝐵₂𝑀_𝑈}, 𝑐₂=inf(𝐴₁,𝐴₂) and 𝛽=2, we can derive 𝐿𝑦,𝑢,𝑡=𝐵₁𝑇_𝑈+𝐵₂𝑀_𝑈−𝐴₁𝑢²₁+𝐴₂𝑢²₂≤𝑐₁−𝑐₂𝑢²₁+𝑢²₂. (3.2)
The proof is complete.

4. Optimally Controlling Chemotherapy

In this section, we discuss the theorem related to the characterization of the optimal control. This result depends on the Pontryagin’s Maximum Principle, which gives necessary conditions for the optimal control. First, we rewrite the system (2.1) in the following vector notation: 𝑑𝑦(𝑡)𝑑𝑡=𝐴𝑦,𝑢,𝑡;∀𝑡>𝑡₀𝑦𝑡,∀𝑢∈𝑈,₀=𝑦₀, (4.1) where 𝑦(𝑡) and 𝐴(𝑦,𝑢,𝑡) are given by 𝑦(𝑡)=𝑇_𝑈(𝑡),𝑇₁(𝑡),𝑇₂(𝑡),𝑀_𝑈(𝑡),𝑀₁(𝑡),𝑀₂(𝑡),𝑉(𝑡)^𝑇,𝐴𝑦,𝑢,𝑡=𝑔₁𝑦,𝑢,𝑡,𝑔₂𝑦,𝑢,𝑡,…,𝑔₆𝑦,𝑢,𝑡,𝑔₇𝑦,𝑢,𝑡^𝑇. (4.2)

The Hamiltonian associated with our problem is 𝐻𝑦,𝑢,𝑝,𝑡=𝐿𝑦,𝑢,𝑡+𝜆(𝑡)^𝑇𝐴𝑦,𝑢,𝑡, (4.3) where the adjoint vector 𝜆(𝑡) is defined by the adjoint equation 𝑑𝜆(𝑡)𝑑𝑡=−𝐴_𝑦𝜆(𝑡)−𝐿_𝑦,𝜆𝑡₁=0. (4.4)

Here 𝐴_𝑦=𝜕𝑔₁𝜕𝑇_𝑈𝜕𝑔₂𝜕𝑇_𝑈𝜕𝑔₃𝜕𝑇_𝑈𝜕𝑔₄𝜕𝑇_𝑈𝜕𝑔₅𝜕𝑇_𝑈𝜕𝑔₆𝜕𝑇_𝑈𝜕𝑔₇𝜕𝑇_𝑈𝜕𝑔₁𝜕𝑇₁𝜕𝑔₂𝜕𝑇₁𝜕𝑔₃𝜕𝑇₁𝜕𝑔₄𝜕𝑇₁𝜕𝑔₅𝜕𝑇₁𝜕𝑔₆𝜕𝑇₁𝜕𝑔₇𝜕𝑇₁𝜕𝑔₁𝜕𝑇₂𝜕𝑔₂𝜕𝑇₂𝜕𝑔₃𝜕𝑇₂𝜕𝑔₄𝜕𝑇₂𝜕𝑔₅𝜕𝑇₂𝜕𝑔₆𝜕𝑇₂𝜕𝑔₇𝜕𝑇₂𝜕𝑔₁𝜕𝑀_𝑈𝜕𝑔₂𝜕𝑀_𝑈𝜕𝑔₃𝜕𝑀_𝑈𝜕𝑔₄𝜕𝑀_𝑈𝜕𝑔₅𝜕𝑀_𝑈𝜕𝑔₆𝜕𝑀_𝑈𝜕𝑔₇𝜕𝑀_𝑈𝜕𝑔₁𝜕𝑀₁𝜕𝑔₂𝜕𝑀₁𝜕𝑔₃𝜕𝑀₁𝜕𝑔₄𝜕𝑀₁𝜕𝑔₅𝜕𝑀₁𝜕𝑔₆𝜕𝑀₁𝜕𝑔₇𝜕𝑀₁𝜕𝑔₁𝜕𝑀₂𝜕𝑔₂𝜕𝑀₂𝜕𝑔₃𝜕𝑀₂𝜕𝑔₄𝜕𝑀₂𝜕𝑔₅𝜕𝑀₂𝜕𝑔₆𝜕𝑀₂𝜕𝑔₇𝜕𝑀₂𝜕𝑔₁𝜕𝑉𝜕𝑔₂𝜕𝑉𝜕𝑔₃𝜕𝑉𝜕𝑔₄𝜕𝑉𝜕𝑔₅𝜕𝑉𝜕𝑔₆𝜕𝑉𝜕𝑔₇=()𝜕𝑉𝐸,𝐹, (4.5)

where 𝐸=−𝛿_𝑇𝑈−𝛽_𝑇𝑉1−𝑢₁𝛾₁𝛽_𝑇𝑉1−𝑢₁𝛾₁00−𝛿_𝑇1+𝑘_𝑇𝑘_𝑇00−𝛿_𝑇2000,000000000𝐹=00000000000𝛼₁𝑁_𝑇1−𝑢₂−𝛿_𝑀𝑈−𝛽_𝑀𝑉1−𝑢₁𝛾₂𝛽_𝑀𝑉1−𝑢₁𝛾₂𝛿000−_𝑀1+𝑘_𝑀𝑘_𝑀000−𝛿_𝑀2𝛼₂𝑁_𝑀1−𝑢₂.000−𝑐 (4.6) In addition, the 𝐿_𝑦 in system (4.3) is 𝐿_𝑦=𝜕𝐿𝜕𝑇_𝑈,𝜕𝐿𝜕𝑇₁,𝜕𝐿𝜕𝑇₂,𝜕𝐿𝜕𝑀_𝑈,𝜕𝐿𝜕𝑀₁,𝜕𝐿𝜕𝑀₂,𝜕𝐿𝜕𝑉^𝑇,=𝐵₁,0,0,𝐵₂,0,0,0^𝑇. (4.7)

Next, adding the penalty term will give us the optimality condition 𝜉𝑦,𝑢,𝜆,𝑡=𝐻𝑦,𝑢,𝜆,𝑡+Γ(𝑢(𝑡))𝜔(𝑡), (4.8) where Γ is an operator from ℝ² to ℝ⁴ defined by Γ(𝑢(𝑡))=1−𝑢₁(𝑡),𝑢₁(𝑡),1−𝑢₂(𝑡),𝑢₂(𝑡),𝜔(𝑡)=𝜔₁₁(𝑡)𝜔₁₂(𝑡)𝜔₂₁(𝑡)𝜔₂₂(𝑡), (4.9) where all 𝜔_𝑖𝑗,𝑖,𝑗=1,2 are nonnegative penalty multipliers satisfying the following conditions: 1−𝑢^∗₁(𝑡)𝜔₁₁(𝑡)=𝑢^∗₁(𝑡)𝜔₁₂(𝑡)=1−𝑢^∗₂(𝑡)𝜔₂₁(𝑡)=𝑢^∗₂(𝑡)𝜔₂₂(𝑡)=0. (4.10)

According to the Pontryagin’s Maximum Principle, if the control 𝑢^∗(𝑡) and the corresponding state 𝑦^∗(𝑡) constitute an optimal pair, there exists an adjoint vector 𝜆(𝑡) defined system (4.4) such that the function 𝜉(𝑦,𝑢,𝜆,𝑡) defined by (4.8) reaches its maximum on the set 𝒰 at the point 𝑢^∗. This gives the following result.

Theorem 4.1. Given an optimal control pair 𝑢^∗(𝑡)=(𝑢^∗₁(𝑡),𝑢^∗₂(𝑡)) and a solution 𝑦^∗(𝑡)=(𝑇^∗_𝑈(𝑡),𝑇^∗₁(𝑡),𝑇^∗₂(𝑡),𝑀^∗_𝑈(𝑡),𝑀^∗₁(𝑡),𝑀^∗₂(𝑡),𝑉^∗(𝑡)) of the corresponding system, then there exist seven adjoint variables 𝜆₁(𝑡),𝜆₂(𝑡),…,𝜆₇(𝑡) satisfying 𝑑𝜆₁=𝛿𝑑𝑡^𝑇𝑈+𝛽_𝑇𝑉1−𝑢₁𝛾₁𝜆₁−𝛽_𝑇𝑉1−𝑢₁𝛾₁𝜆₂−𝐵₁,𝑑𝜆₂=𝛿𝑑𝑡^𝑇1+𝑘_𝑇𝜆₂−𝑘_𝑇𝜆₃,𝑑𝜆₃𝑑𝑡=𝛿_𝑇2𝜆₃−𝛼₁𝑁_𝑇1−𝑢₂𝜆₇,𝑑𝜆₄=𝛿𝑑𝑡^𝑀𝑈+𝛽_𝑀𝑉1−𝑢₁𝛾₂𝜆₄−𝛽_𝑀𝑉1−𝑢₁𝛾₂𝜆₅−𝐵₂,𝑑𝜆₅=𝛿𝑑𝑡^𝑀1+𝑘_𝑀𝜆₅−𝑘_𝑀𝜆₆,𝑑𝜆₆𝑑𝑡=𝛿_𝑀2𝜆₆−𝛼₂𝑁_𝑀1−𝑢₂𝜆₇,𝑑𝜆₇𝑑𝑡=𝑐𝜆₇+𝛽_𝑇𝑇_𝑈1−𝑢₁𝛾₁𝜆₁−𝛽_𝑇𝑇_𝑈1−𝑢₁𝛾₁𝜆₂+𝛽_𝑀𝑀_𝑈1−𝑢₁𝛾₂𝜆₄−𝛽_𝑀𝑀_𝑈1−𝑢₁𝛾₂𝜆₅, (4.11) with the final conditions 𝜆₁𝑡₁=𝜆₂𝑡₁=⋯=𝜆₇𝑡₁=0. (4.12) Furthermore, 𝑢^∗₁(𝑡)=min{max{0,𝑅₁(𝑡)},1},𝑢^∗₂(𝑡)=min{max{0,𝑅₂(𝑡)},1}, where 𝑅₁(𝑡)=𝑉^∗2𝐴₁𝛽_𝑇𝑇^∗_𝑈𝛾₁𝜆₁−𝜆₂+𝛽_𝑀𝑀^∗_𝑈𝛾₂𝜆₄,𝑅−𝜆5₂(𝑡)𝜆=−⁷2𝐴₂𝛼₁𝑁_𝑇𝑇₂+𝛼₂𝑁_𝑀𝑀₂. (4.13)

Proof. According to the previous section, an optimal couple (𝑦^∗(𝑡),𝑢^∗(𝑡)) exists for maximizing the objective functional (2.2) subject to the system (2.1). Therefore, by Pontryagin’s Maximum Principle, there exists a vector 𝜆(𝑡)=(𝜆₁(𝑡),…,𝜆₇(𝑡))^𝑇 satisfying 𝜆(𝑡)𝑑𝑡=−𝜕𝐻𝜕𝑦=−𝐿_𝑦−𝐴_𝑦𝜆(𝑡). (4.14) That yields 𝜆₁(𝑡)𝑑𝑡=−𝜕𝑔₁𝑦^∗,𝑢^∗,𝑡𝜕𝑇_𝑈,…,𝜕𝑔₇𝑦^∗,𝑢^∗,𝑡𝜕𝑇_𝑈𝜆(𝑡)−𝜕𝐿𝑦^∗,𝑢^∗,𝑡𝜕𝑇_𝑈,𝜆₂(𝑡)𝑑𝑡=−𝜕𝑔₁𝑦^∗,𝑢^∗,𝑡𝜕𝑇₁,…,𝜕𝑔₇𝑦^∗,𝑢^∗,𝑡𝜕𝑇₁𝜆(𝑡)−𝜕𝐿𝑦^∗,𝑢^∗,𝑡𝜕𝑇₁,𝜆₃(𝑡)𝑑𝑡=−𝜕𝑔₁𝑦^∗,𝑢^∗,𝑡𝜕𝑇₂,…,𝜕𝑔₇𝑦^∗,𝑢^∗,𝑡𝜕𝑇₂𝜆(𝑡)−𝜕𝐿𝑦^∗,𝑢^∗,𝑡𝜕𝑇₂,𝜆₄(𝑡)𝑑𝑡=−𝜕𝑔₁𝑦^∗,𝑢^∗,𝑡𝜕𝑀_𝑈,…,𝜕𝑔₇𝑦^∗,𝑢^∗,𝑡𝜕𝑀_𝑈𝜆(𝑡)−𝜕𝐿𝑦^∗,𝑢^∗,𝑡𝜕𝑀_𝑈,𝜆₅(𝑡)𝑑𝑡=−𝜕𝑔₁𝑦^∗,𝑢^∗,𝑡𝜕𝑀₁,…,𝜕𝑔₇𝑦^∗,𝑢^∗,𝑡𝜕𝑀₁𝜆(𝑡)−𝜕𝐿𝑦^∗,𝑢^∗,𝑡𝜕𝑀₂,𝜆₆(𝑡)𝑑𝑡=−𝜕𝑔₁𝑦^∗,𝑢^∗,𝑡𝜕𝑀₂,…,𝜕𝑔₇𝑦^∗,𝑢^∗,𝑡𝜕𝑀₂𝜆(𝑡)−𝜕𝐿𝑦^∗,𝑢^∗,𝑡𝜕𝑀₂,𝜆₇(𝑡)𝑑𝑡=−𝜕𝑔₁𝑦^∗,𝑢^∗,𝑡𝜕𝑉,…,𝜕𝑔₇𝑦^∗,𝑢^∗,𝑡𝜆(𝑡)−𝑦𝜕𝑉𝜕𝐿^∗,𝑢^∗,𝑡.𝜕𝑉 (4.15)
Through simple calculations, we derive system (4.11). The Pontryagin’s Maximum Principle gives the following necessary conditions to obtain the optimal pair (𝑦^∗,𝑢^∗): 𝜕𝜉𝑦^∗,𝑢^∗,𝜆,𝑡𝜕𝑢₁𝑦=0,𝜕𝜉^∗,𝑢^∗,𝜆,𝑡𝜕𝑢₂=0, (4.16) where 𝜉(𝑦^∗,𝑢^∗,𝜆,𝑡)=𝐻(𝑦^∗,𝑢^∗,𝜆,𝑡)+Γ(𝑢^∗(𝑡))𝜔(𝑡). From (4.10) and (4.16), we have 𝜕𝜉𝑦^∗,𝑢^∗,𝜆,𝑡𝜕𝑢₁⟹𝑦=0𝜕𝐿^∗,𝑢^∗,𝑡𝜕𝑢₁+𝜕𝜆(𝑡)^𝑇𝐴𝑦^∗,𝑢^∗,𝑡𝜕𝑢₁+𝜕Γ𝑢^∗(𝑡)𝜔(𝑡)𝜕𝑢₁=0, (4.17) which implies 𝑢^∗₁(𝑡)=𝑉^∗2𝐴₁𝛽_𝑇𝑇^∗_𝑈𝛾₁𝜆₁−𝜆₂+𝛽_𝑀𝑀^∗_𝑈𝛾₂𝜆₄−𝜆5. (4.18) On the other hand, 𝜕𝜉𝑦,𝑢^∗,𝜆,𝑡𝜕𝑢₂⟹𝑦=0𝜕𝐿^∗,𝑢^∗,𝑡𝜕𝑢₂+𝜕𝜆(𝑡)^𝑇𝐴𝑦,𝑢^∗,𝑡𝜕𝑢₂+𝜕Γ𝑢^∗(𝑡)𝜔(𝑡)𝜕𝑢₂=0, (4.19) which indicates 𝑢^∗₂(𝑡)𝜆=−⁷2𝐴₂𝛼₁𝑁_𝑇𝑇₂+𝛼₂𝑁_𝑀𝑀₂. (4.20)
Now from the constraint condition, the following three cases arise. Case 1. 𝑡∈{𝑡∶0<𝑢^∗₁(𝑡)<1} and 𝜔₁₁(𝑡)=𝜔₁₂(𝑡)=0. Then 𝑢^∗₁(𝑡)=𝑅₁(𝑡).Case 2. 𝑡∈{𝑡∶𝑢^∗₁(𝑡)=0} and 𝜔₁₁(𝑡)=0. Then 0=𝑢^∗₁(𝑡)=𝑅₁(𝑡)+𝜔₁₂(𝑡), which implies 𝑢^∗₁(𝑡)≥𝑅₁(𝑡) because 𝜔₁₂(𝑡)≥0.Case 3. 𝑡∈{𝑡∶𝑢^∗₁(𝑡)=1} and 𝜔₁₂(𝑡)=0. Then 𝑢^∗₁(𝑡)=𝑅₁(𝑡)−𝜔₁₁(𝑡), which leads to 1=𝑢^∗₁(𝑡)≤𝑅₁(𝑡), owing to 𝜔₁₁≥0.Hence, we have 𝑢^∗₁(𝑡)=min{max{0,𝑅₁(𝑡)},1}. Similarly, we can get that 𝑢^∗₂(𝑡)=min{max{0,𝑅₂(𝑡)},1}.
The proof is complete.

Now, the optimality system is given by incorporating the optimal control pair in the state system coupled with the adjoint system. Thus, we have𝑑𝑦^∗(𝑡)𝑦𝑑𝑡=𝐴^∗,𝑢^∗,𝑡;∀𝑡>𝑡₀,𝑑𝜆(𝑡)𝑑𝑡=−𝐴_𝑦∗𝜆(𝑡)−𝐿_𝑦∗,𝑦^∗𝑡₀=𝑦^∗₀,𝜆𝑡₁=0. (4.21) We substitute the expressions 𝑢^∗=(𝑢^∗₁,𝑢^∗₂) in the above system. The uniqueness of the solution of the optimality system can be derived by a standard method (refer to [6] for more details on the proof).

5. Numerical Algorithm and Results

The resolution of the optimal system is created improving the Gauss-Seidel-like implicit finite-difference method developed by [7] and denoted by GSS1 method. It consists on discretizing the interval [𝑡₀,𝑡₁] at the points 𝑡_𝑘=𝑘𝑙+𝑡₀(𝑘=0,1,…,𝑛), where 𝑙 is the time step.

In the following, we define the state and adjoint variables 𝑇_𝑈(𝑡), 𝑇₁(𝑡), 𝑇₂(𝑡), 𝑀_𝑈(𝑡), 𝑀₁(𝑡), 𝑀₂(𝑡), 𝑉(𝑡), 𝜆₁(𝑡)∼𝜆₇(𝑡) and the controls 𝑢₁(𝑡) and 𝑢₂(𝑡) in terms of nodal points 𝑇^𝑘_𝑈, 𝑇^𝑘₁, 𝑇^𝑘₂, 𝑀^𝑘_𝑈, 𝑀^𝑘₁, 𝑀^𝑘₂, 𝑉^𝑘, 𝜆^𝑘₁∼𝜆^𝑘₇, 𝑢^𝑘₁, 𝑢^𝑘₂ as the state and adjoint variables and the controls at initial time 𝑡₀, while 𝑇^𝑛_𝑈, 𝑇^𝑛₁, 𝑇^𝑛₂, 𝑀^𝑛_𝑈, 𝑀^𝑛₁, 𝑀^𝑛₂, 𝑉^𝑛, 𝜆^𝑛₁∼𝜆^𝑛₇, 𝑢^𝑛₁, 𝑢^𝑛₂ as the state and adjoint variables and the controls at final time 𝑡₁. As it is well known that the approximation of the time derivative by its first-order forward-difference is given, for the first state variable 𝑇_𝑈, by 𝑑𝑇_𝑈(𝑡)𝑑𝑡=lim_𝑙→0𝑇_𝑈(𝑡+𝑙)−𝑇_𝑈(𝑡)𝑙. (5.1)

We use the scheme developed by Gumel et al. [7] in the following way: 𝑇_𝑈^𝑘+1−𝑇^𝑘_𝑈𝑙=𝜃_𝑇−𝛿_𝑇𝑈𝑇_𝑈^𝑘+1−𝛽_𝑇𝛾₁𝑉^𝑘1−𝑢^𝑘₁𝑇_𝑈^𝑘+1. (5.2) Analogously, we have 𝑇₁^𝑘+1−𝑇^𝑘₁𝑙=𝛽_𝑇𝛾₁𝑉^𝑘1−𝑢^𝑘₁𝑇_𝑈^𝑘+1−𝛿_𝑇1+𝑘_𝑇𝑇₁^𝑘+1,𝑇₂^𝑘+1−𝑇^𝑘₂𝑙=𝑘_𝑇𝑇₁^𝑘+1−𝛿_𝑇2𝑇₂^𝑘+1,𝑀_𝑈^𝑘+1−𝑀^𝑘_𝑈𝑙=𝜃_𝑀−𝛿_𝑀𝑈𝑀_𝑈^𝑘+1−𝛽_𝑀𝛾₂𝑉^𝑘1−𝑢^𝑘₁𝑀_𝑈^𝑘+1,𝑀₁^𝑘+1−𝑀^𝑘₁𝑙=𝛽_𝑀𝛾₂𝑉^𝑘1−𝑢^𝑘₁𝑀_𝑈^𝑘+1−𝛿_𝑀1+𝑘_𝑀𝑀₁^𝑘+1,𝑀₂^𝑘+1−𝑀^𝑘₂𝑙=𝑘_𝑀𝑀₁^𝑘+1−𝛿_𝑀2𝑀₂^𝑘+1,𝑉^𝑘+1−𝑉^𝑘𝑙=𝛼₁𝑁_𝑇𝑇₂^𝑘+11−𝑢^𝑘₂+𝛼₂𝑁_𝑀𝑀₂^𝑘+11−𝑢^𝑘₂−𝑐𝑉^𝑘+1. (5.3)

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}𝑙=𝛿_𝑇𝑈+𝛽_𝑇𝑉^𝑘+11−𝑢^𝑘₁𝛾₁𝜆₁^{𝑛−𝑘−1}−𝛽_𝑇𝑉^𝑘+11−𝑢^𝑘₁𝛾₁𝜆₂^𝑛−𝑘−𝐵₁,𝜆₂^𝑛−𝑘−𝜆₂^{𝑛−𝑘−1}𝑙=𝛿_𝑇1+𝑘_𝑇𝜆₂^{𝑛−𝑘−1}−𝑘_𝑇𝜆₃^𝑛−𝑘,𝜆₃^𝑛−𝑘−𝜆₃^{𝑛−𝑘−1}𝑙=𝛿_𝑇2𝜆₃^{𝑛−𝑘−1}−𝛼₁𝑁_𝑇1−𝑢^𝑘₂𝜆₇^𝑛−𝑘,𝜆₄^𝑛−𝑘−𝜆₄^{𝑛−𝑘−1}𝑙=𝛿_𝑀𝑈+𝛽_𝑀𝑉^𝑘+11−𝑢^𝑘₁𝛾₂𝜆₄^{𝑛−𝑘−1}−𝛽_𝑀𝑉^𝑘+11−𝑢^𝑘₁𝛾₂𝜆₅^𝑛−𝑘−𝐵₂,𝜆₅^𝑛−𝑘−𝜆₅^{𝑛−𝑘−1}𝑙=𝛿_𝑀1+𝑘_𝑀𝜆₅^{𝑛−𝑘−1}−𝑘_𝑀𝜆₆^𝑛−𝑘,𝜆₆^𝑛−𝑘−𝜆₆^{𝑛−𝑘−1}𝑙=𝛿_𝑀2𝜆₆^{𝑛−𝑘−1}−𝛼₂𝑁_𝑀1−𝑢^𝑘₂𝜆₇^𝑛−𝑘,𝜆₇^𝑛−𝑘−𝜆₇^{𝑛−𝑘−1}𝑙=𝑐𝜆₇^{𝑛−𝑘−1}+𝛽_𝑇𝑇_𝑈^𝑘+11−𝑢^𝑘₁𝛾₁𝜆₁^{𝑛−𝑘−1}−𝜆₂^{𝑛−𝑘−1}+𝛽_𝑀𝑀_𝑈^𝑘+11−𝑢^𝑘₁𝛾₂𝜆₄^{𝑛−𝑘−1}−𝜆₅^{𝑛−𝑘−1}. (5.4)

Hence, we can establish an algorithm to solve the optimality system and then to compute the optimal control pair by employing the GSS1 method (5.2)–(5.4) that we denote by IGSS1 method here (see Appendix B).

5.1. 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 𝑢^∗₁(𝑡), 𝑢^∗₂(𝑡). For the following parameters and initial values:

𝜃_𝑇=10, 𝜃_𝑀=10, 𝛿_𝑇𝑈=0.02, 𝛿_𝑇1=0.5, 𝛿_𝑇2=1, 𝛿_𝑀𝑈=0.0495, 𝛿_𝑀1=0.0495, 𝛿_𝑀2=0.0495, 𝛽_𝑇=0.00008, 𝛽_𝑀=0.00008, 𝑘_𝑇=0.1, 𝑘_𝑀=0.1, 𝑁_𝑇=100, 𝑁_𝑀=100, 𝑇⁰_𝑈=1000, 𝑇⁰₁=0, 𝑇⁰₂=0, 𝑀⁰_𝑈=1000, 𝑀⁰₁=0, 𝑀⁰₂=0, 𝑉⁰=1, 𝑐=0.03.

The experimental results obtained are listed in Table 1 (in which “before treatment’’ and “after treatment’’ are denoted by BT and AT, resp.).

Table 1 

For more clearness, it is better to present these comparative results by the following graphs. When the viruses attack the human body, uninfected T and M cells decrease (see Figures 3 and 4).

Figure 3 

Uninfected T cells.

Figure 4 

Uninfected M cells.

The viruses do not stop to proliferate and so its abundance dramatically increases (see Figure 5). However, after introducing the optimal controls, the situation changes. A few days later, the effect of chemotherapy starts to appear; which explains the growth of uninfected T and M cells and the diminishing of viruses (see Figure 6).

Figure 5 

Virus population before optimal controls.

Figure 6 

Virus population after optimal controls.

Finally, the optimal controls 𝑢^∗₁(𝑡), 𝑢^∗₂(𝑡) for drug administration are presented through Figures 7 and 8.

Figure 7 

Optimal control variable 𝑢∗1(𝑡).

Figure 8 

Optimal control variable 𝑢∗2(𝑡).

6. Conclusions

By incorporating the chemotherapy into a previous model describing the interaction of the immune system with the human immunodeficiency virus (HIV), this paper has proposed a novel HIV virus spread model with control variables. Our goal is to maximize the number of healthy cells and, meanwhile, to minimize the cost of chemotherapy. In this context, the existence of an optimal control has been proved. Experimental results show that, under this model, the spread of HIV virus can be controlled effectively.

Our next work is to study other kinds of models, especially those with impulsive drug effect.

Appendices

A. The Theorem Used in Theorem 3.1

The equations 𝑥̇(𝑡)=𝑓(𝑡,𝑥(𝑡)),𝑥|_𝑡=𝜏=𝜉, (A.1) where (𝜏,𝜉)∈𝐷, with 𝐷 a nonempty open subset of ℛ×ℛ^𝑛 and 𝑓∶𝐷→ℛ^𝑛, are called a Cauchy problem or initial-value problem.

A solution to the Cauchy Problem is defined to be any pair (𝐼,𝜙) in which 𝐼 is an open subinterval of ℛ containing 𝜏,𝜙∶𝐼→ℛ^𝑛 is absolutely continuous, (𝑡,𝜙(𝑡))∈𝐷 for all 𝑡∈𝐼, and 𝜙 satisfies the above two equations at 𝑎.𝑒.𝑡∈𝐼.

For 𝑥∈ℛ^𝑛 with coordinates 𝑥_𝑖, define a norm on ℛ^𝑛 by |𝑥|=max_{1≤𝑖≤𝑛}𝑥_𝑖. (A.2)

The following theorem applies the Lebesgue integral and the hypothesis is stated in terms of the rectangular subset of ℛ×ℛ^𝑛 centered about (𝜏,𝜉), 𝑅_𝑎,𝑏=(𝑡,𝑥)∶|𝑡−𝜏|≤𝑎,𝑥−𝜉≤𝑏,𝑎>0,𝑏>0. (A.3)

Theorem A.1 (see [19, p.182]). The Cauchy problem has a solution if for some 𝑅_𝑎,𝑏⊂𝐷 centered about (𝜏,𝜉) the restriction of 𝑓 to 𝑅_𝑎,𝑏 is continuous in 𝑥 for fixed 𝑡, measurable in 𝑡 for fixed 𝑥, and satisfies 𝑓(𝑡,𝑥)≤𝑚(𝑡),(𝑡,𝑥)∈𝑅_𝑎,𝑏, (A.4) for some 𝑚 integrable over the interval [𝜏−𝑎,𝜏+𝑎].

B. An Algorithm Using the GSS1 Method

Algorithm B.1. Step 1. 𝑇_𝑈𝑡₀⟵𝑇⁰_𝑈,𝑇₁𝑡₀⟵𝑇⁰₁,𝑇₂𝑡₀⟵𝑇⁰₂,𝑀_𝑈𝑡₀⟵𝑀⁰_𝑈,𝑀₁𝑡₀⟵𝑀⁰₁,𝑀₂𝑡₀⟵𝑀⁰₂,𝑉𝑡₀⟵𝑉⁰,𝜆₁𝑡_𝑛⟵0,𝜆₂𝑡_𝑛⟵0,𝜆₃𝑡_𝑛⟵0,𝜆₄𝑡_𝑛𝜆⟵0,₅𝑡_𝑛⟵0,𝜆₆𝑡_𝑛⟵0,𝜆₇𝑡_𝑛𝑢⟵0,₁𝑡₀⟵0,𝑢₂𝑡₀⟵0. (B.1)Step 2. for 𝑘=1,…,𝑛 do 𝑇^𝑘_𝑈⟵𝑙𝜃_𝑇+𝑇_𝑈^𝑘−11+𝑙𝛿_𝑇𝑈+𝑙𝛽_𝑇𝑉^𝑘−11−𝑢₁^𝑘−1𝛾₁,𝑇^𝑘₁⟵𝑙𝛽_𝑇𝛾₁𝑉^𝑘−11−𝑢₁^𝑘−1𝑇^𝑘_𝑈+𝑇₁^𝑘−11+𝛿_𝑇1+𝑘_𝑇𝑙,𝑇^𝑘₂⟵𝑙𝑘_𝑇𝑇^𝑘₁+𝑇₂^𝑘−11+𝑙𝛿_𝑇2,𝑀^𝑘_𝑈⟵𝑙𝜃_𝑀+𝑀_𝑈^𝑘−11+𝑙𝛿_𝑀𝑈+𝑙𝛽_𝑀𝛾₂𝑉^𝑘−11−𝑢₁^𝑘−1,𝑀^𝑘₁⟵𝑙𝛽_𝑀𝛾₂𝑉^𝑘−11−𝑢₁^𝑘−1𝑀^𝑘_𝑈+𝑀₁^𝑘−11+𝛿_𝑀1+𝑘_𝑀𝑙,𝑀^𝑘₂⟵𝑙𝑘_𝑀𝑀^𝑘₁+𝑀₂^𝑘−11+𝑙𝛿_𝑀2,𝑉^𝑘⟵𝑙𝛼₁𝑁_𝑇𝑇^𝑘₂1−𝑢₂^𝑘−1+𝑙𝛼₂𝑁_𝑀𝑀^𝑘₂1−𝑢₂^𝑘−1+𝑉^𝑘−1,𝜆1+𝑐𝑙_1𝑛−𝑘⟵𝑙𝐵₁+𝑙𝛽_𝑇𝑉^𝑘1−𝑢₁^𝑘−1𝛾₁𝜆₂^{𝑛−𝑘+1}+𝜆₁^{𝑛−𝑘+1}1+𝑙𝛿_𝑇𝑈+𝑙𝛽_𝑇𝑉^𝑘1−𝑢₁^𝑘−1𝛾₁,𝜆₂^𝑛−𝑘⟵𝜆₂^{𝑛−𝑘+1}+𝑙𝑘_𝑇𝜆₃^{𝑛−𝑘+1}1+𝛿_𝑇1+𝑘_𝑇𝑙,𝜆₃^𝑛−𝑘⟵𝜆₃^{𝑛−𝑘+1}+𝑙𝛼₁𝑁_𝑇1−𝑢₂^𝑘−1𝜆₇^{𝑛−𝑘+1}1+𝑙𝛿_𝑇2,𝜆₄^𝑛−𝑘⟵𝑙𝐵₂+𝑙𝛽_𝑀𝑉^𝑘1−𝑢₁^𝑘−1𝛾₂𝜆₅^{𝑛−𝑘+1}+𝜆₄^{𝑛−𝑘+1}1+𝑙𝛿_𝑀𝑈+𝑙𝛽_𝑀𝑉^𝑘1−𝑢₁^𝑘−1𝛾₂,𝜆₅^𝑛−𝑘⟵𝜆₅^{𝑛−𝑘+1}+𝑙𝑘_𝑀𝜆₆^{𝑛−𝑘+1}1+𝛿_𝑀1+𝑘_𝑀𝑙,𝜆₆^𝑛−𝑘⟵𝜆₆^{𝑛−𝑘+1}+𝑙𝛼₂𝑁_𝑀1−𝑢₂^𝑘−1𝜆₇^{𝑛−𝑘+1}1+𝑙𝛿_𝑀2,𝜆₇^𝑛−𝑘⟵𝜆₇^{𝑛−𝑘+1}+𝑙1−𝑢₁^𝑘−1𝛽_𝑇𝑇^𝑘_𝑈𝛾₁𝜆₂^𝑛−𝑘−𝜆₁^𝑛−𝑘+𝛽_𝑀𝑀^𝑘_𝑈𝛾₂𝜆₅^𝑛−𝑘−𝜆₄^𝑛−𝑘,𝑅1+𝑙𝑐_𝑘1⟵𝑉^𝑘2𝐴₁𝛽_𝑇𝑇^𝑘_𝑈𝛾₁𝜆₁^𝑛−𝑘−𝜆₂^𝑛−𝑘+𝛽_𝑀𝑀^𝑘_𝑈𝛾₂𝜆₄^𝑛−𝑘−𝜆₅^𝑛−𝑘,𝑅^𝑘₂𝜆⟵−^7𝑛−𝑘2𝐴₂𝛼₁𝑁_𝑇𝑇^𝑘₂+𝛼₂𝑁_𝑀𝑀^𝑘₂,𝑢^𝑘₁⟵minmax0,𝑅^𝑘₁,1,𝑢^𝑘₂⟵minmax0,𝑅^𝑘₂,1. (B.2)Step 3. for 𝑘=1,…,𝑛, denote 𝑇^∗_𝑈𝑡_𝑘=𝑇^𝑘_𝑈,𝑀^∗_𝑈𝑡_𝑘=𝑀^𝑘_𝑈,𝑢^∗₁𝑡_𝑘=𝑢^𝑘₁,𝑢^∗₂𝑡_𝑘=𝑢^𝑘₂. (B.3)

It is easy to conclude that this algorithm takes 𝑂(𝑛) execution time.

Acknowledgments

This work was financially supported by the Fundamental Research Funds for the Central Universities (no. CDJXS10181130), the New Century Excellent Talent Project of China (no. NCET-05-0759), and the National Natural Science Foundation of China (no. 10771227).