Abstract

This present study proposes an optimal control problem, with the final goal of implementing an optimal treatment protocol which could maximize the survival time of patients and minimize the cost of drug utilizing a system of ordinary differential equations which describes the interaction of the immune system with the human immunodeficiency virus (HIV). Optimal control problem transfers into a modified problem in measure space using an embedding method in which the existence of optimal solution is guaranteed by compactness of the space. Then the metamorphosed problem is approximated by a linear programming (LP) problem, and by solving this LP problem a suboptimal piecewise constant control function, which is more practical from the clinical viewpoint, is achieved. The comparison between the immune system dynamics in treated and untreated patients is introduced. Finally, the relationships between the healthy cells and virus are shown.

1. Introduction

Human immunodeficiency virus infects CD4+ T-cells, which are an important part of the human immune system, and other target cells. The infected cells produce a large number of viruses. Medical treatments for HIV have greatly improved during the last two decades. Highly active antiretroviral therapy (HAART) allows for the effective suppression of HIV-infected individuals and prolongs the time before the onset of acquired immune deficiency syndrome (AIDS) for years or even decades and increases life expectancy and quality to the patient. But antiretroviral therapy cannot eradicate HIV from infected patients because of the long-lived infected cells and sites within the body where drugs may not achieve effective levels [13]. HAART contains two major types of anti-HIV drugs: reverse transcriptase inhibitors (RTI), and protease inhibitors (PI). Reverse transcriptase inhibitors prevent HIV from infecting cells by blocking the integration of the HIV viral code into the host cell genome while protease inhibitors prevent infected cells from replication of infectious virus particles, and can reduce and maintain viral load below the limit of detection in many patients. Moreover, treatment with either type of drugs can also increase the CD4+ T-cell counts that are target cells for HIV.

Many of the host-pathogen interaction mechanisms during HIV infection and progression to AIDS are still unknown. Mathematical modeling of HIV infection is of interest to the medical community as no adequate animal models exist to test the efficacy of drug regimes. These models can test different assumptions and provide new insights into questions that are difficult to answer by clinical or experimental studies. A number of mathematical models have been formulated to describe various aspects of the interaction of HIV with healthy cells. Some of these models are addressed in [4]. The basic model of HIV infection is presented by Wodarz and Nowak [5], which contains three state variables: healthy CD4+ T-cells, infected CD4+ T-cells, and concentration of free virus. Their model has been modified to offer important theoretical insights into immune control of the virus, based on treatment strategies, while maintaining a simple structure [6]. Furthermore, this modified model has been developed to guess the natural evolution of HIV infection, as qualitatively described in several clinical studies [7].

Some authors have used mathematical models for HIV infection in conjunction with control theory to achieve appropriate goals. For example, these goals may include maximizing the level of healthy CD4+ T-cells and minimizing the cost of treatment [811], maximizing the level of healthy CD4+ T-cells while minimizing both the cost of treatment and viral load [12], minimizing both the HIV population and systemic costs to body while maximizing immune response [13, 14], and maximizing both the healthy CD4+ T-cell counts and immune response while minimizing the cost of treatment [15], maximizing the healthy CD4+ T-cell counts and minimizing both the side effects and drug resistance [16].

The papers [1721] consider only RTI medication while the papers [22, 23] consider only PIs. In [2427], all effects of a HAART medication are combined to one control variable in the model. In [2832], dynamical multidrug therapies based on RTIs and PIs are designed.

In this paper, a mathematical model of HIV dynamic is considered that includes the effect of antiretroviral therapy, and an analysis of optimal control is performed regarding appropriate goals.

The paper is organized as follows: in Section 2, the underlying HIV mathematical model is described. Our formulation of the control problem, which attempts to prolong the survival time of patient as long as possible, is described in Section 3. Approximating the obtained optimal control problem by an LP problem is the subject of Section 4. Numerical results obtained from solving the LP problem are presented in Section 5. Finally, Section 6 is assigned to concluding remarks.

2. Presentation of a Working Model

In this paper, the pathological behavior of HIV is considered which is modeled with the simplified version of a system of ordinary differential equations (ODEs) as described in [17]. This model, which is consistent with clinical data, is given as follows:𝑑𝑃(𝑡)𝑑𝑡=𝐼𝑃𝑃+𝛽0𝑃(𝑡)𝜏𝑝𝑃(𝑡)𝑐𝑃𝑎(𝑡)𝐶(𝑡)𝑃(𝑡),(1)𝑑𝑃(𝑡)=𝑑𝑡𝜏𝑃𝑃(𝑡)𝜏𝑃𝑃(𝑡)𝑐𝑃𝑎(𝑡)𝐶(𝑡)𝑃(𝑡),(2)𝑑𝑎(𝑡)𝑑𝑡=𝑎(𝑡)(𝜅𝜃𝛾𝐶(𝑡)),(3)𝑑𝐶(𝑡)𝑑𝑡=𝑎(𝑡)𝜀𝐼𝐶+𝛼𝐶(𝑡)𝑃(𝑡)𝑃0𝜐𝜏𝐶𝐶(𝑡).(4) Most of the terms in this model have straightforward interpretations. 𝑃() and 𝑃() denote the amounts of immature CD4+ T-cells and mature ones, respectively. The term 𝑎() indicates HIV particles, and 𝐶() designates cytotoxic T-cells specific for HIV (CTLs) as a function of time. Here, 𝐼𝑃 is the constant rate that 𝑃 cells are produced, 𝜏𝑝 is the rate of maturation of 𝑃 cells into 𝑃 cells, and 𝜏𝑝 is the rate of natural death of 𝑃 cells. Furthermore, 𝛽 is the amplifying coefficient of the linear feedback effect of 𝑃 cells decrease on the influx of 𝑃 cells at time 𝑡. Free virus particles 𝑎(𝑡) eliminate 𝑃() cells at a rate proportional to 𝑐𝑃𝑎(𝑡)𝐶(𝑡)𝑃(𝑡) at time 𝑡. Similarly, 𝑐𝑃𝑎(𝑡)𝐶(𝑡)𝑃(𝑡) is the rate of elimination of 𝑃 cells. The term 𝜃 characterizes the growth rate of HIV particles, and 𝛾 is the rate of inactivation of HIV products mediated by cytotoxic 𝐶 cells. 𝐼𝐶 is the influx of 𝐶 cell precursors, 𝜀 is their maturation rate, 𝛼 is the proliferation rate of 𝐶 cells under the antigenic stimulation by HIV products, and 𝜏𝐶 is their natural death rate. Helper T-cells effect on maturation and proliferation of 𝐶 cells is expressed by the ratio 𝑃(𝑡)/𝑃0, and 𝜈 is introduced to characterize the intensity of this helper effect. Chemotherapeutic agent was simulated by decreasing the value 𝜅, that is, the HIV proliferation rate. Lower value for 𝜅 corresponds to higher RTI-drug doses.

3. Optimal Control Formulation

In this section, we formulate an optimal control problem that identifies the inhibition parameter 𝜅 in (3), with a function of the control variable. In particular, we will replace the parameter 𝜅 with the function 1𝑢(𝑡). This choice then identifies the control variable 𝑢(𝑡) with the rate of inhibition of virus reproduction, which is modeled as a simple function of drug dosage.

In clinical practice, the following guidelines are used typically.(i)Antiretroviral therapy is initiated at 𝑡0, the time at which the CD4+ T-cell count falls below 350 cells/μL.(ii)The transition from HIV to AIDS is marked by a CD4+ T-cell count below 200 cells/μL.(iii)A person is said to have full-blown AIDS when his/her CD4+ T-cell count falls below CD4+crit, typically around of 50 cells/μL.

This paper aims to propose a drug regimen that delays the onset of full-blown AIDS and prolongs survival as much as possible, while one is going to minimize the drug costs. This can be modeled as follows.

Assume that the onset of full-blown AIDS occurs after time 𝑡𝑓. Hence, we should have𝑃(𝑡)CD4+crit𝑡,𝑡0,𝑡𝑓𝑡,𝑃𝑓=CD4+crit.(5) A problem arising from the use of most chemotherapies is the multiple and sometimes harmful side effects, as well as the ineffectiveness of treatment after a certain time due to the capability of the virus to mutate and become resistant to the treatment. Global effects of these phenomena can be considered by imposing limited treatment interval [22], that is, treatment lasting for a given period from time 𝑡0 to 𝑡0+𝜂. Therefore, the support of the control function 𝑢() must be in the treatment interval 𝑡supp𝑢0,𝑡0+𝜂.(6) Here, we follow [8, 22] in assuming that the costs of the treatment is proportional to 𝑢2(𝑡) at time 𝑡. Therefore, the overall cost of the treatment is 𝑡𝑓𝑡0𝑢2(𝑡)𝑑𝑡. So, the following functional should be maximized:𝜎𝑡𝑓,𝑢=𝑡𝑓𝜆𝑡𝑓𝑡0𝑢2(𝑡)𝑑𝑡.(7) Parameter 𝜆 is used to set the relative importance between maximizing the survival time 𝑡𝑓 and minimizing the systemic cost to the body. Setting 𝑃=𝑥1, 𝑃=𝑥2, 𝑎=𝑥3, and 𝐶=𝑥4, the system of differential equations (1)–(4) can be represented in a generalized form as=𝐼̇𝑥(𝑡)=𝑔(𝑡,𝑥(𝑡),𝑢(𝑡))𝑃𝑥+𝛽02𝑥2𝜏𝑃𝑥1𝑐𝑃𝑥3𝑥4𝑥1𝜏𝑃𝑥1𝜏𝑃𝑥2𝑐𝑃𝑥3𝑥4𝑥2𝑥3(1𝑢)𝜃𝛾𝑥4𝜀𝐼𝐶+𝛼𝑥4𝑥3𝑥2𝑥02𝜐𝜏𝐶𝑥4,𝑥(0)=𝑥0.(8) Assume that 𝐾 denotes the set of all measurable control functions 𝑢()[0,1], where 𝑢() satisfies (6), and the corresponding solution of (8) at final time 𝑡𝑓 satisfies (5). Therefore, we are seeking for 𝑢()𝐾 such that 𝜎𝑡𝑓𝑡,𝑢𝜎𝑓,𝑢,𝑢𝐾.(9) Setting 𝑓0(𝑡,𝑥(𝑡),𝑢(𝑡))=1𝜆𝑢2(𝑡), then the optimal drug regimen problem, while ignoring 𝑡0, can be represented as:max𝑡𝑓,𝑢𝐾𝑡𝑓𝑡0𝑓0(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡(10) subject to𝑥𝑡̇𝑥=𝑔(𝑡,𝑥(𝑡),𝑢(𝑡)),(11)0=𝑥𝑡0,𝑥2𝑡𝑓=CD4+crit𝑥,(12)2(𝑡)CD4+crit𝑡,𝑡0,𝑡𝑓.(13)

This optimal control problem is referred to as OCP. Some problems may arise in the quest of solving OCP. The set 𝐾 may be empty. If 𝐾 is not empty, the functional measuring the performance of the system may not achieve its maximum in the set 𝐾. In order to overcome these difficulties, in the next section we transfer the OCP into a modified problem in measure space.

4. Approximation of OCP by Linear Programming Problem

Using measure theory for solving optimal control problems based on the idea of Young [33], which was applied for the first time by Wilson and Rubio [34], has been theoretically established by Rubio in [35]. Then, the method has been extended for approximating the time optimal problems by an LP model [36]. Here, this approach is used.

4.1. Functional Space

We assume that the state variables 𝑥() and the control input 𝑢(), respectively, get their values in the compact sets 𝐴=𝐴1×𝐴2×𝐴3×𝐴44 and 𝑈. Set 𝐽=[𝑡0,𝑡𝑓]. Here, we are going to derive weak forms for (11)–(13).

Definition 1. A triple 𝑝=[𝑡𝑓,𝑥,𝑢] is said to be admissible if the following conditions hold.(i)The vector function 𝑥() is absolutely continuous and belongs to 𝐴 for all 𝑡𝐽.(ii)The function 𝑢() takes its values in the set 𝑈 and is Lebesgue measurable on 𝐽.(iii)𝑝 satisfies in (11)–(13), on 𝐽0, that is, the interior set of 𝐽.It is assumed that the set of all admissible triples is nonempty and denotes it by 𝑊. Let 𝑝 be an admissible triple, 𝐵 be an open ball in 5 containing 𝐽×𝐴, and let 𝐶(𝐵) be the space of all real-valued continuous differentiable functions on it. Let 𝜑𝐶(𝐵), and define 𝜑𝑔 as follows: 𝜑𝑔(𝑡,𝑥(𝑡),𝑢(𝑡))=𝜑𝑥(𝑡,𝑥(𝑡))𝑔(𝑡,𝑥(𝑡),𝑢(𝑡))+𝜑𝑡(𝑡,𝑥(𝑡))(14) for each [𝑡,𝑥(𝑡),𝑢(𝑡)]Ω, where Ω=𝐽×𝐴×𝑈. The function 𝜑𝑔 is in the space 𝐶(Ω), the set of all continuous functions on the compact set Ω. Since 𝑝=[𝑡𝑓,𝑥,𝑢] is an admissible triple, we have 𝑡𝑓𝑡0𝜑𝑔(𝑡,𝜉(𝑡),𝑢(𝑡))𝑑𝑡=𝑡𝑓𝑡0𝜑𝑥(𝑡,𝑥(𝑡))̇𝑥(𝑡)+𝜑𝑡(𝑡𝑡,𝑥(𝑡))𝑑𝑡=𝜑𝑓𝑡,𝑥𝑓𝑡𝜑0𝑡,𝑥0=Δ𝜑,(15) for all 𝜑𝐶(𝐵). Let 𝐷(𝐽0) be the space of all infinitely differentiable real-valued functions with compact support in 𝐽0. Define 𝜓𝑛(𝑡,𝑥(𝑡),𝑢(𝑡))=𝑥𝑛(𝑡)𝜓(𝑡)+𝑔𝑛𝐽(𝑡,𝑥(𝑡),𝑢(𝑡))𝜓(𝑡),𝑛=1,2,3,4,𝜓𝐷0.(16) Assume 𝑝=[𝑡𝑓,𝑥,𝑢] be an admissible triple. Since the function 𝜓() has compact support in 𝐽0, 𝜓(𝑡0)=𝜓(𝑡𝑓)=0. Thus, for 𝑛=1,2,3,4, and for all 𝜓𝐷(𝐽0), from (16) and using integration by parts, we have 𝑡𝑓𝑡0𝜓𝑛(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡=𝑡𝑓𝑡0𝑥𝑛(𝑡)𝜓(+𝑡)𝑑𝑡𝑡𝑓𝑡0𝑔𝑛(𝑡,𝑥(𝑡),𝑢(𝑡))𝜓(𝑡)𝑑𝑡=0.(17) Also, by choosing the functions which are dependent only on time, we have 𝑡𝑓𝑡0𝜗(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡=𝑎𝜗,𝜗𝐶1(Ω),(18) where 𝐶1(Ω) is the space of all functions in 𝐶(Ω) that depend only on time and 𝑎𝜗 is the integral of 𝜗() on 𝐽.
Equations (15), (17), and (18) are the weak forms of (11)–(13). Note that the constraints (12) are considered on the right-hand side of (15) by choosing suitable functions 𝜑𝐶(𝐵) which are monomials of 𝑥2. Furthermore, the constraint (13) is considered, by choosing an appropriate set 𝐴. Now, we consider the following positive linear functional on 𝐶(Ω): Γ𝑝𝐹𝐽𝐹(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡,𝐹𝐶(Ω).(19)

Proposition 1. Transformation 𝑝Γ𝑝 of admissible triples in 𝑊 into the linear mappings Γ𝑝 defined in (19) is an injection.

Proof. We must show that if 𝑝1𝑝2, then Γ𝑝1Γ𝑝2. Let 𝑝𝑗=[𝑡𝑓𝑗,𝑥𝑗,𝑢𝑗],𝑗=1,2 be different admissible triples. If 𝑡𝑓1=𝑡𝑓2, then there is a subinterval of [𝑡0,𝑡𝑓1], say 𝐽1, where 𝑥1(𝑡)𝑥2(𝑡) for each 𝑡𝐽1. A continuous function 𝐹 can be constructed on Ω so that the right-hand sides of (19) corresponding to 𝑝1 and 𝑝1 are not equal. For instance, one can make 𝐹 independent of 𝑢, equal zero for all 𝑡 outside 𝐽1, and such that it is positive on the appropriate portion of 𝑥1(), and zero on the 𝑥2(), then the linear functionals are not equal. In other words, if 𝑡𝑓1𝑡𝑓2, then Γ𝑝1 and Γ𝑝2 have different domains and are not equal.

Thus, from (15), (17), and (18), one can conclude that maximizing the functional (10) over admissible space 𝑊, changes to the following optimization problem in functional space:max𝑝𝑊Γ𝑝𝑓0(20) subject toΓ𝑝(𝜑𝑔)=Δ𝜑,𝜑𝐶Γ(𝐵),(21)𝑝(𝜓𝑛𝐽)=0,𝑛=1,2,3,4,𝜓𝐷0Γ,(22)𝑝(𝜗)=𝑎𝜗,𝜗𝐶1(Ω).(23)

4.2. Measure Space

Let 𝑀+(Ω) denote the space of all positive Radon measures on Ω. By the Riesz representation theorem [35], there exists a unique positive Radon measure 𝜇 on Ω such thatΓ𝑝(𝐹)=𝐽=𝐹(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡Ω𝐹(𝑡,𝑥,𝑢)𝑑𝜇𝜇(𝐹),𝐹𝐶(Ω).(24) So, we may change the functional space of the optimization problem to measure space. In other words, the optimization problem (20)–(23) can be converted to the following optimization problem in measure space:Maximize𝜇𝑀+(Ω)𝜇𝑓0(25) subject to𝜇(𝜑𝑔)=Δ𝜑,𝜑𝐶(𝐵),(26)𝜇(𝜓𝑛𝐽)=0,𝑛=1,2,3,4,𝜓𝐷0,(27)𝜇(𝜗)=𝑎𝜗,𝜗𝐶1(Ω).(28)

We will consider maximization of (25) over the set 𝑄 of all positive Radon measures on Ω, satisfying (26)–(28). The main advantages of considering this measure theoretic form of the problem is the existence of optimal measure in the set 𝑄 where this point can be studied in a straightforward manner without having to impose conditions such as convexity which may be artificial.

Define function 𝐼𝑄𝑅 as 𝐼(𝜇)=𝜇(𝑓0). The following theorem guarantees the existence of an optimal solution.

Theorem 1. The measure theoretical problem of maximizing (25) with constraints (26)–(28) has an optimal solution, say 𝜇, where 𝜇𝑄.

Proof. The so-called constraints (27) and (28) are special cases of (26) [35]. So, the set Q can be written as 𝑄=𝜑𝐶(𝐵)𝜇𝑀+(Ω)𝜇(𝜑𝑔)=Δ𝜑.(29) Assume that 𝑝=[𝑡𝑓,𝑥,𝑢] is an admissible triple. It is well known that the set {𝜇𝑀+(Ω)𝜇(1)=𝑡𝑓𝑡0} is compact in the weak* topology. Furthermore, the set 𝑄 as intersection of inverse image of closed singleton sets {Δ𝜑} under the continuous functions 𝜇𝜇(𝜑𝑔) is also closed. Thus, 𝑄 is a closed subset of a compact set. This proves the compactness of the set 𝑄. Since the functional 𝐼, mapping the compact set 𝑄 on the real line, is continuous and thus takes its maximum on the compact set 𝑄.

Next, based on analysis in [35], the problem (25)–(28) is approximated by an LP problem, and a triple 𝑝 which approximates the action of 𝜇𝑄 is achieved.

4.3. Approximation

The problem (25)–(28) is an infinite-dimensional linear programming problem, and we are mainly interested in approximating it. First, the maximization of 𝐼 is considered not over the set 𝑄, but over a subset of it denoted by requiring that only a finite number of constraints (26)–(28) be satisfied. Let {𝜑𝑖𝑖=1,2,}, {𝜓𝑗𝑗=1,2,}, and {𝜗𝑠𝑠=1,2,} be the sets of total functions, respectively, in 𝐶(𝐵), 𝐷(𝐽0), and 𝐶1(Ω). The first approximation is completed by choosing finite number of functions 𝜑𝑖s, 𝜓𝑗s, and 𝜗𝑠s. Now we have the following propositions.

Proposition 2. Consider the linear program problem consisting of maximizing the function 𝐼 over the set 𝑄𝑀 of measures in 𝑀+(Ω) satisfying: 𝜇𝜑𝑔𝑖=Δ𝜑𝑖,𝑖=1,,𝑀.(30) Then, 𝐽𝑀max𝑄𝑀𝐼 tends to 𝐽=max𝑄𝐼 as 𝑀.

Proof. We have 𝑄1𝑄2𝑄𝑀𝑄 and hence, 𝐽1𝐽2𝐽𝑀𝐽. The sequence {𝐽𝑗}𝑗=1 is nonincreasing and bounded, so, it converges to a number 𝜁 such that 𝜁𝐽. We show that 𝜁=𝐽. Set 𝑅𝑀=1𝑄𝑀. Then, 𝑅𝑄 and 𝜁max𝑅𝐼. It is sufficient to show 𝑅𝑄. Assume 𝜇𝑅 and 𝜑𝐶(𝐵). Since the linear combinations of the functions {𝜑𝑗,𝑗=1,2,} are uniformly dense in 𝐶(𝐵), there is a sequence {𝜑𝑘}span{𝜑𝑗,𝑗=1,2,}, such that 𝜑𝑘 tends to 𝜑 uniformly as 𝑘. Hence, 𝑆1, 𝑆2, and 𝑆3 tend to zero as 𝑘 where 𝑆1=sup|𝜑𝑥𝜑𝑘𝑥|, 𝑆2=sup|𝜑𝑡𝜑𝑘𝑡|, and 𝑆3=sup|𝜑𝜑𝑘|. Since 𝜇𝑅 and the functional 𝑓𝜇(𝑓) is linear, 𝜇(𝜑𝑘𝑔)=Δ𝜑𝑘 and ||𝜇(𝜑𝑔||=||)Δ𝜑𝜇(𝜑𝑔)Δ𝜑𝜇𝜑𝑘𝑔+Δ𝜑𝑘||=||||Ω𝜑𝑥(𝑡,𝑥)𝜑𝑘𝑥𝑔+𝜑(𝑡,𝑥)(𝑡,𝑥,𝑢)𝑡(𝑡,𝑥)𝜑𝑘𝑡(𝑡,𝑥)𝑑𝜇Δ𝜑Δ𝜑𝑘||||𝑆1Ω||||𝑔(𝑡,𝑥,𝑢)𝑑𝜇+𝑆2Ω𝑑𝜇+2𝑆3.(31) The right-hand side of the above inequality tends to zero as 𝑘, and the left-hand side is independent of 𝑘; therefore 𝜇(𝜑𝑔)=Δ𝜑. Thus, 𝑅𝑄 and 𝜁𝐽, which implies 𝜁=𝐽.

Proposition 3. The measure 𝜇 in the set 𝑄𝑀 at which the functional 𝐼 attains its maximum has the form 𝜇=𝑀𝑗=1𝛼𝑗𝛿𝑧𝑗,(32) where 𝛼𝑗0, 𝑧𝑗Ω, and 𝛿(𝑧) is unitary atomic measure with the support being the singleton set {𝑧𝑗}, characterized by 𝛿(𝑧)(𝐹)=𝐹(𝑧), 𝑧Ω.

Proof . See [35].

Therefore, our attention is restricted to finding a measure in the form of (32), which maximizes the functional 𝐼 and satisfies in 𝑀 number of the constraints (26)–(28). Thus, by choosing the functions 𝜑𝑖, 𝑖=1,2,,𝑀1, 𝜓𝑘, 𝑘=1,,𝑀2, and 𝜗𝑠, 𝑠=1,,𝑆, the infinite dimensional problem (25)–(28) is approximated by the following finite dimensional nonlinear programming (NLP) problem: Maximize𝛼𝑗0,𝑧𝑗Ω𝑀𝑗=1𝛼𝑗𝑓0𝑧𝑗(33) subject to𝑀𝑗=1𝛼𝑗𝜑𝑔𝑖𝑧𝑗=Δ𝜑𝑖,𝑖=1,,𝑀1,(34)𝑀𝑗=1𝛼𝑗𝜓𝑛𝑘𝑧𝑗=0,𝑘=1,,𝑀2,𝑛=1,2,3,4,(35)𝑀𝑗=1𝛼𝑗𝜗𝑠𝑧𝑗=𝑎𝜗𝑠,𝑠=1,,𝑆,(36) where 𝑀=𝑀1+4𝑀2+𝑆. Clearly, (33)–(36) is an NLP problem with 2 𝑀 unknowns: 𝛼𝑗 and 𝑧𝑗, 𝑗=1,,𝑀. One is interested in LP problem. The following proposition enables us to approximate the NLP problem (33)–(36) by a finite dimensional LP problem.

Proposition 4. Let Ω𝑁={𝑦1,𝑦2,,𝑦𝑁} be a countable dense subset of Ω. Given 𝜀>0, a measure 𝑣𝑀+(Ω) can be found such that: ||𝑣𝑓0𝜇𝑓0||||𝑣𝜑𝜀,𝑔𝑖𝜇𝜑𝑔𝑖||𝜀,𝑖=1,,𝑀1,||𝑣𝜓𝑛𝑘𝜇𝜓𝑛𝑘||𝜀,𝑘=1,,𝑀2||𝑣𝜗,𝑛=1,2,3,4,𝑠𝜇𝜗𝑠||𝜀,𝑠=1,,𝑆,(37) where the measure 𝑣 has the form 𝑣=𝑀𝑗=1𝛼𝑗𝛿𝑦𝑗,(38) and the coefficients 𝛼𝑗, 𝑗=1,,𝑀, are the same as optimal measure (32), and 𝑦𝑗Ω𝑁, 𝑗=1,,𝑀.

Proof. We rename the functions 𝑓0, 𝜑𝑔𝑖’s, 𝜓𝑛𝑘’s, and 𝜗𝑠’s sequentially as 𝑗, 𝑗=1,2,,𝑀+1. Then, for 𝑗=1,,𝑀+1, ||𝜇𝑣𝑗||=|||𝑀𝑖=1𝛼𝑖𝑗𝑧𝑖𝑗𝑦𝑖|||𝑀𝑖=1𝛼𝑖max𝑖,𝑗||𝑗𝑧𝑖𝑗𝑦𝑖||.(39)𝑗s are continuous. Therefore, max𝑖,𝑗 can be made less than 𝜀/𝑀𝑗=1𝛼𝑗 by choosing 𝑦𝑖, 𝑖=1,2,,𝑀, sufficiently near 𝑧𝑖.

For constructing a suitable set Ω𝑁, which preserves the relation (6), 𝐽 is divided to 𝑆 subintervals as follows:𝐽𝑠=𝑡0+(𝑠1)Δ𝑇𝑆1,𝑡0+𝑠Δ𝑇,𝑆1𝑠=1,2,,𝑆1,𝐽𝑆=𝑡𝑙,𝑡𝑓,(40) where 𝑡𝑙 is a lower bound for optimal time 𝑡𝑓, which can be obtained by using a search algorithm based on golden section [36] or Fibonnaci search method [37]. Let 𝑆 be the largest number such that 𝐽𝑆[𝑡0,𝑡0+𝜂]. Set 𝐽1=𝑆𝑠=1𝐽𝑠, 𝐽2=𝑆𝑠=𝑆+1𝐽𝑠, Ω1=𝐽1×𝐴×𝑈, and Ω2=𝐽2×𝐴×{0}. Moreover, the intervals 𝐴𝑖 (𝑖=1,2,3,4) and 𝑈 are divided, respectively, into 𝑛𝑖 and 𝑚 subintervals. So, the sets Ω𝑖, 𝑖=1,2, are partitioned into 𝑁1=𝑆𝑛1𝑛2𝑛3𝑛4𝑚 and 𝑁2=(𝑆𝑆)𝑛1𝑛2𝑛3𝑛4 cells, respectively. One point is chosen from each cell. In this way, we will have a grid of points, which are numbered sequentially as 𝑦𝑗=(𝑡𝑗,𝑥1𝑗,,𝑥4𝑗,𝑢𝑗), 𝑗=1,,𝑁, where 𝑁=𝑁1+𝑁2.

Therefore, according to (38), the NLP problem (33)–(36) is converted to the following LP problem:Maximize𝛼𝑗0𝑀𝑗=1𝛼𝑗𝑓0𝑦𝑗(41) subject to𝑁𝑗=1𝛼𝑗𝜑𝑔𝑖𝑦𝑗=Δ𝜑𝑖,𝑖=1,,𝑀1,(42)𝑁𝑗=1𝛼𝑗𝜓𝑛𝑘𝑦𝑗=0,𝑘=1,,𝑀2,𝑛=1,2,3,4,(43)𝑁𝑗=1𝛼𝑗𝜗𝑠𝑦𝑗=𝑎𝜗𝑠,𝑠=1,,𝑆.(44)

Here, we discuss suitable total functions 𝜑𝑖s, 𝜓𝑘s, and 𝜗𝑠s. The functions 𝜑𝑖s can be taken to be monomials of 𝑡 and the components of the vector 𝑥 as follows:𝑡𝑖𝑥𝑗2,𝑥𝑗2𝑥𝑖,𝑖{0,1},𝑗{1,2,},{1,3,4}.(45) In addition, we choose some functions with compact support in the following form [36, 37]: 𝜓2𝑟1(𝑡)=sin2𝜋𝑟𝑡𝑡0Δ𝑇𝑡𝑡𝑙𝜓0otherwise,2𝑟(𝑡)=1cos2𝜋𝑟𝑡𝑡0Δ𝑇𝑡𝑡𝑙0otherwise,(46) where 𝑟=1,2, and Δ𝑇=𝑡𝑙𝑡0. Finally, the following functions are considered that are dependent on 𝑡 only:𝜗𝑠(𝑡)=1𝑡𝐽𝑠0otherwise,(47) where 𝐽𝑠, 𝑠=1,,𝑆, are given by (40). These functions are used to construct the approximate piecewise constant control [3537]. By the above definition of 𝜗𝑠, we consider 𝑡𝑓 as an unknown variable in the constraints (44) which can be written as𝑗=1𝛼𝑗=Δ𝑇𝑆1(𝑆1)𝑗=(𝑆2)+1𝛼𝑗=Δ𝑇𝑆1𝑆𝑗=(𝑆1)+1𝛼𝑗=𝑡𝑓𝑡𝑙,(48) where =𝑁/𝑆. Of course, we need only to construct the control function 𝑢(), since 𝑥() can be obtained by solving the ODEs (8). By using simplex method, a nonzero optimal solution 𝛼𝑖1,𝛼𝑖2,,𝛼𝑖𝑘, 𝑖1<𝑖2<<𝑖𝑘 of the LP problem (41)–(44) can be found where 𝑘 cannot exceed the number of constraints, that is, 𝑘𝑀1+𝑀2+𝑆. Setting 𝛼𝑖0=𝑡0, a piecewise control function 𝑢() approximating the optimal control is constructed based on these nonzero coefficients as follows [35, 36]:𝑢𝑢(𝑡)=𝑖𝑗𝑡𝑗1=0𝛼𝑖,𝑗=0𝛼𝑖0otherwise,𝑗=1,2,,𝑘,(49) where 𝑢𝑖𝑗is the 6th component of 𝑦𝑖𝑗.

To start the proposed method, one needs to have 𝑡𝑙. Here, a bisection method is proposed to find the desired lower bound 𝑡𝑙 for optimal time 𝑡𝑓. This algorithm has a simple structure and is started with a given upper bound 𝑡𝑢, where it is assumed that the lower bound starts with 𝑡𝑙=0. Assuming that 𝑡𝑓(𝑡𝑙) denotes the solution of LP problem (41)–(44) corresponding to the given lower bound 𝑡𝑙, the bisection method is outlined as follows.

Algorithm 1 (estimation of the lower bound 𝑡𝑙). First, let 𝜏=[𝑡𝑙,𝑡𝑢], where 𝑡𝑙=0 and 𝑡𝑢 is an upper bound for 𝑡𝑓.
Step 1. Let 𝑎=(𝑡𝑙+𝑡𝑢)/2 and solve the corresponding LP problem to find 𝑡𝑓(𝑎). If no feasible solution is found for the corresponding LP problem or 𝑡𝑓(𝑎)=𝑎, set 𝑡𝑢=𝑎; else set 𝑡𝑙=𝑎.
Step 2. If the length of the interval 𝜏=[𝑡𝑙,𝑡𝑢] is small enough, then choose 𝑡𝑙 as a good estimation for lower bound 𝑡𝑓 else, go to Step 1.

5. Numerical Results

In this implementation, we set 𝑀1=10 and choose functions 𝜑𝑖, 𝑖=1,2,,𝑀1, from 𝐶(𝐵) as follows:𝑥1,𝑥2,𝑥3,𝑥4,𝑥22,𝑥32,𝑥1𝑥2,𝑥3𝑥2,𝑥4𝑥2,𝑡𝑥2.(50) Furthermore, we set 𝑆=13 and 𝑀2=6. Setting 𝑢=0 in (8), we find that at 𝑡0=1642, 𝑥(𝑡0)=(8.19,35,0.05,0.04). Model parameters are chosen as follows [17]:𝐼𝑃=1.0,𝛽=0.01,𝜏𝑝=0.2,𝑐𝑝=0,𝜏𝑝𝑐=0.001,𝑝=20,𝜃=0.02,𝛾=0.8,𝜀=0.154,𝐼𝐶=0.2,𝛼=0.3,𝑣=3,(51) and the following initial condition is used:𝑥(0)=(5,100,0.0005,0).(52) Besides, 𝜆 is set to 𝜆=10, and the length of treatment is set to 𝜂=500 (days). By using controllability on the dynamical control system, one can assume 𝐴1=[7,10], 𝐴2=[5,85], 𝐴3=[0,5], 𝐴4=[0,0.1], and 𝑈=[0,1]. Furthermore, the number of partitions in the construction of the set Ω𝑁 are 𝑛1=4, 𝑛2=10, 𝑛3=4, 𝑛4=4, and 𝑚1=4. Initially, an upper bound for optimal time 𝑡𝑓 is set to 𝑡𝑢=5297.45. The results of implementing Algorithm 1 are summarized in Table 1. Setting 𝑡𝑙=3642, we have an LP problem with 16000 unknowns and 47 constraints which is solved by the linprog code of the optimization toolbox in MATLAB. The total CPU time required on a laptop with CPU 2.20 GHz and 0.99 GB of RAM was 17.23 minutes. The suboptimal time has been found 𝑡𝑓=3943.2. The resulting suboptimal control and the response of the system to the obtained control function are depicted in Figures 1 and 2, respectively. Moreover, we found 𝑃(𝑡𝑓)=4.9360, which is close to the exact value, that is, 5 (𝐶𝐷4+crit%). Note that the normal level of mature CD4+ T-cells is about 1000 cells/μL. The relationships between the CD4+ T-cells, CTLs, and virus during the different stages of the disease are shown in Figure 3 as a phase space diagram.

6. Conclusion

In this paper, we considered a dynamical system which describes the various aspects of the interaction of HIV with the immune system, to construct an optimal control problem which maximizes survival time of patients. A measure theoretical method is used to solve such kind of problems. The method is not iterative, and it does not need any initial guess of the solution, and numerical results confirmed the effectiveness of this approach.

Numerical results show that in presence of treatment, the survival time of patients can be considerably prolonged. From Figures 2(b) and 2(c), it is concluded that in presence of treatment (solid lines), the virus is controlled to very low levels and CD4+ T-cells are maintained at high levels for relatively long time. From Figure 2(d), an increase in CTL’s occurs in response to therapy.

Figure 3(a) shows an inverse correlation between CD4+ T-cells and virus particles. Furthermore, Figure 3(b) shows a clear correlation between the level of CTLs in the blood and HIV progression. As the virus increases upon initial infection, CTLs increase in order to decrease the virus. But this situation changes after about 1000th day due to destruction of CD4+ T-cells. Because these cells play an essential role in stimulation of immune response and signal other immune cells to eliminate infection by killing infected cells. After the 1642nd day, an increase in immune response can be observed which is due to recovery of CD4+ T-cells in response to treatment. Immune response increases for a while after discontinuation of therapy but ultimately becomes extinct.