Multiobjective Optimal Control of HIV Dynamics

Abstract

Various aspects of the interaction of HIV with the human immune system can be modeled by a system of ordinary differential equations. This model is utilized, and a multiobjective optimal control problem (MOOCP) is proposed to maximize the CD4+ T cells population and minimize both the viral load and drug costs. The weighted sum method is used, and continuous Pareto optimal solutions are derived by solving the corresponding optimality system. Moreover, a model predictive control (MPC) strategy is applied, with the final goal of implementing Pareto optimal structured treatment interruptions (STI) protocol. In particular, by using a fuzzy approach, the MOOCP is converted to a single-objective optimization problem to derive a Pareto optimal solution which among other Pareto optimal solutions has the best satisfaction performance. Then, by using an embedding method, the problem is transferred into a modified problem in an appropriate space in which the existence of solution is guaranteed by compactness of the space. The metamorphosed problem is approximated by a linear programming (LP) model, and a piecewise constant solution which shows the desired combinations of reverse transcriptase inhibitor (RTI) and protease inhibitor (PI) drug efficacies is achieved.

1. Introduction

Human Immunodeficiency Virus (HIV) 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 for the patient but eradication of HIV infection does not seem possible with currently available antiretroviral drugs. This is due primarily to the establishment of a pool of latently infected T-cells and sites within the body where drugs may not achieve effective levels [1–3]. 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. 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 drug can also increase the CD4+ T-cell count 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 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 [4]. The basic model of HIV infection is presented by Perelson et al. [5] that contains three state variables: healthy CD4+ T-cells, infected CD4+ T-cells, and concentration of free virus. Another model is presented in [6] that although maintaining a simple structure, it offers important theoretical insights into immune control of the virus based on treatment strategies. Furthermore, this modified model is developed to describe the natural evolution of HIV infection, as qualitatively described in several clinical studies [7].

The problem of designing drug administration in HIV infected patients using mathematical models can be considered as multi-objective optimal control problems. For example, these objectives may include maximizing the level of healthy CD4+ T-cells and minimizing the cost of treatment [8–15], maximizing immune response and minimizing both the cost of treatment and viral load [16–21], maximizing both the level of healthy CD4+ T-cells and immune response and minimizing the cost of treatment [22], maximizing the level of healthy CD4+ T-cells while minimizing both the side effects and drug resistance [23], and maximizing survival time of patient subject to drug cost [24].

When a multi-objective problem is treated, each objective conflicts with one another and, unlike a single-objective optimization, the solution to this problem is not a single point, but a family of solutions known as the Pareto-optimal set. Among these solutions, the designer should find the best compromise taking into proper account the attributes and characteristics of the handled problem.

There exists a wide variety of methods that can be used to compute Pareto optimal points. A widely used technique consists of reducing the multi-objective problem to a single-objective one by means of “scalarization” procedure. The weighted sum (WS) method is a commonly used scalarization technique which consists of assigning each objective function a weight coefficient and then optimizing the function obtained by summing up all the objective functions scaled by their weight coefficients. Nevertheless, it has intrinsic drawbacks. The WS method is highly scale dependent, an equal distribution of weights does not necessarily lead to an even spread along the Pareto front, and that points in a nonconvex part of the Pareto front cannot be obtained [25].

Recent scalar multiple objective optimization techniques such as Normal Boundary Intersection (NBI) [26] and Normalized Normal Constraint (NNC) [27] have been found to mitigate the disadvantages of the WS method. Recently, these methods have been successfully combined with direct optimal control approaches for the efficient solution of multi-objective optimal control problems. For example, in [28], a successful application of NBI and NNC for the multiple objective optimal control of (bio) chemical processes has been reported, and in [29] several scalarization techniques for multi-objective optimization, for example, WS, NNC, and NBI have been integrated with fast deterministic direct optimal control approaches.

The papers [8, 13, 22, 24, 30, 31] consider only RTI medication while the papers [14, 32] consider only PIs. In [6, 19–21], all effects of a HAART medication are combined to one control variable in the model. In [7, 10–12, 15, 17, 18, 33–35], dynamical multidrug therapies based on RTIs and PIs are designed.

In the considered control approaches, the amount of medications can be either continuous or on-off type. This treatment is also known as structured treatment interruption (STI). STI has received considerable attentions as it might reduce the risk of HIV mutating to strains which are resistant to current medication regimens. STI approach might also reduce possible long-term toxicity of the drugs. A concise summary of clinical STI studies, including protocols and results, is presented in [35].

In this paper, we consider a mathematical model of HIV dynamics that includes the effect of antiretroviral therapy, and we perform analysis of optimal control regarding maximizing the CD4+ T-cell counts and minimizing both the viral load and cost of drugs.

The paper is organized as follows: in Section 2, the underlying HIV mathematical model is presented. Our formulation of the MOOCP is described in Section 3. In Section 4, the weighted sum method is used, and the MOOCP is converted to a single-objective optimal control problem. First, we assume that the cost of the drug regime varies quadratically with the amount of drug administered, and then we characterize the continuous Pareto optimal solutions using Pontryagin’s Maximum Principle. Secondly, it is shown that the presence of linear cost functionals leads to STI-type Pareto optimal solutions, and an MPC-based technique is proposed to find this kind of solutions. In Section 5, a fuzzy approach is utilized and the MOOCP is converted to a single-objective optimization problem, with the final goal of finding a Pareto optimal solution which has the best satisfaction performance among other Pareto optimal solutions. Moreover, the converted problem is approximated by an LP model applying a measure-theoretical approach, and a piecewise constant solution is achieved. Some numerical experiments are provided in Section 6. The last section is the conclusion.

2. Presentation of a Working Model

In this paper, the pathogenesis of HIV is modeled with a system of ordinary differential equations (ODEs) described in [7]. This model can be viewed as an extension of basic HIV Model of Perelson et al. [5]̇̇̇̇̇𝑥=𝑠−𝑑𝑥−𝑟𝑥𝑣,𝑦=𝑟𝑥𝑣−𝑎𝑦−𝑓𝑦𝑧,𝑤=𝑐𝑥𝑦𝑤−𝑞𝑦𝑤−𝑏𝑤,𝑧=𝑞𝑦𝑤−ℎ𝑧,𝑣=𝑘1−𝑢_𝑃̇𝑦−𝜏𝑣,𝑟=𝑟₀−𝑢_𝑅. (2.1)(2.2)(2.3)(2.4)(2.5)(2.6) Most of the terms in the model have straightforward interpretations as following.

The first equation represents the dynamics of the concentration of healthy CD4+ T-cells (𝑥). The healthy CD4+ T-cells are produced from a source, such as the thymus, at a constant rate 𝑠, and die at a rate 𝑑𝑥. The cells are infected by the virus at a rate 𝑟𝑥𝑣. The second equation describes the dynamics of the concentration of infected CD4+ T-cells (𝑦). The infected CD4+ T-cells result from the infection of healthy CD4+ T-cells and die at a rate 𝑎𝑦 and are killed by cytotoxic T-lymphocyte effectors CTLe (𝑧) at a rate 𝑓𝑦𝑧. The population of cytotoxic T-cells is divided into precursors or CTLp (𝑤), and effectors or CTLe (𝑧). Equations (2.3) and (2.4) describe the dynamics of these compartments. In accordance with experimental, findings establishment of a lasting CTL response depends on CD4+ T-cell help, and that HIV impairs T-helper cell function. Thus, proliferation of the CTLp population is given by 𝑐𝑥𝑦𝑤 and is proportional to both virus load (𝑦) and the number of uninfected T-helper cells (𝑥). CTLp differentiation into effectors occurs at a rate 𝑞𝑦𝑤. Finally, CTLe dies at a rate ℎ𝑧. Equation (2.5) describes the dynamics of the free virus particles (𝑣). These free virus particles are produced from infected CD4+ T-cells at a rate 𝑘𝑦 and are cleared at a rate 𝜏𝑣. The model also contains an index of the intrinsic virulence or aggressiveness of the virus (𝑟). This index increases linearly in the case of an untreated HIV-infected individual, with a growth rate that depends on the constant 𝑟₀. Finally, (2.6) describes the dynamic of this index. In model variables 𝑢 and 𝑢 denote protease inhibitors (PI) and reverse transcriptase inhibitors (RTI), respectively. The effect of PI drugs is modeled by reducing the proliferation rate of viruses from infected cells, while the effect of RTI drugs is modeled by reducing the infection rate, and in this way, blocking the infection of CD4+ T-cells by free virus. Hence, in this model the RTI drugs have an effect on virulence because their main role is halting cellular infection and preventing virus production by reducing the production rate from infected CD4+ T cells.

The model has several parameters that must be assigned for numerical simulations. The descriptions, numerical values, and units of the parameters are summarized in Table 1. These descriptions and values are taken from [7]. We note that (2.1)–(2.6) with these parameters model dynamics of fast progressive patients (FPP).

Table 1 

Parameter values for the HIV model.

3. Multiobjective Optimal Control Formulation

A problem arising from the use of most chemotherapy 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. Although we do not intend to model effects of resistance or side effects, we impose a condition called a limited treatment window, that monitors the global effects of these phenomena. The treatment lasts for a given period from time 𝑡₀ to 𝑡_𝑓, say. In clinical practice, antiretroviral therapy is initiated at 𝑡₀, the time at which CD4+ T-cell counts reach 350 cells/μL. We would like to maximize levels of healthy CD4+ T-cells, minimize level of viral load, and also we want to minimize the systemic costs of treatments. Thus, the following objective functionals are to be maximized simultaneously𝐼₁(𝑢)=^𝑡𝑓_𝑡0𝑥(𝑡)𝑑𝑡,𝐼₂(𝑢)=−^𝑡𝑓_𝑡0𝑣(𝑡)𝑑𝑡,𝐼₃(𝑢)=−^𝑡𝑓_𝑡0𝑢^𝑚_𝑃(𝑡)𝑑𝑡,𝐼₄(𝑢)=−^𝑡𝑓_𝑡0𝑢^𝑚_𝑃(𝑡)𝑑𝑡, (3.1) where 𝑥(𝑡) and 𝑣(𝑡) are the solution of ODEs (2.1)–(2.6) corresponding to control pair 𝑢=(𝑢_𝑃,𝑢_𝑅) and 𝑚 is a positive integer. Let 𝑈 be the set of all measurable control pairs (𝑢_𝑃,𝑢_𝑅) that 𝑎₁≤𝑢_𝑃(𝑡)≤𝑏₁,𝑎₂≤𝑢_𝑅(𝑡)≤𝑏₂,∀𝑡∈𝑡₀,𝑡_𝑓. (3.2) Therefore, the optimal drug regimen problem can be represented asMaximize_𝑢∈𝑈𝐼₁(𝑢),𝐼₂(𝑢),𝐼₃(𝑢),𝐼₄(𝑢). (3.3) In general, there does not exist a pair of control functions (𝑢_𝑃,𝑢_𝑅)∈𝑈 that renders the maximum value to each functional 𝐼_𝑖,𝑖=1,2,3,4, simultaneously, and one uses the concept of Pareto optimality in the sense of following definition.

Definition 3.1. A pair 𝑢^∗=(𝑢^∗_𝑃,𝑢^∗_𝑅)∈𝑈 is said to be an Pareto optimal solution of the problem (3.3) if, and only if, there exist no (𝑢_𝑃,𝑢_𝑅)∈𝑈 such that 𝐼_𝑖(𝑢^∗)≤𝐼_𝑖(𝑢) for all 𝑖∈{1,2,3,4}, and 𝐼_𝑖(𝑢^∗)<𝐼_𝑖(𝑢) for some 𝑖∈{1,2,3,4}.

As seen from Definition 3.1, in general there exist an infinite number of Pareto optimal solutions. Actually, we should select one control function among the set of Pareto optimal solutions. We are going to find a Pareto optimal solution of problem (3.3), combining the techniques and methods from the classic multi-objective optimization and control theory.

4. The Weighted Sum Method

The weighted sum method scalarizes set of objectives into a single-objective by multiplying each objective with user supplied weights. The weight of an objective is usually chosen in proportion to the objective’s relative importance in the problem. After, a composite objective functional 𝐼(𝑢) can be formed by summing the weighted objectives and the MOOCP given in (3.3) is converted to a single-objective optimization problem as follows:Maximize_𝑢∈𝑈𝐼(𝑢)=⁴_𝑗=1𝑤_𝑗𝐼_𝑗(𝑢). (4.1)

Theorem 4.1. The solution of the weighted sum problem (4.1) is Pareto optimal if the weighting coefficients are positive, that is, 𝑤_𝑖>0 for all 𝑖=1,2,3,4.

Proof. Assume 𝑢^∗∈𝑈 is the optimal solution of problem (4.1). If 𝑢^∗ is not a Pareto optimal solution, then there is 𝑢∈𝑈 such that 𝐼_𝑖(𝑢^∗)≤𝐼_𝑖(𝑢) for all 𝑖∈{1,2,3,4}, and 𝐼_𝑗(𝑢^∗)<𝐼_𝑗(𝑢) for some 𝑗∈{1,2,3,4}. Since 𝑤_𝑖>0 for all 𝑖=1,2,3,4, thus, 𝐼(𝑢^∗)<𝐼(𝑢), and this contradicts the optimality of 𝑢^∗.

4.1. Continuous Solutions

Here, we followed [8, 32] in assuming that systemic costs of the PI and RTI drugs treatment are proportional to 𝑢²_𝑃(𝑡) and 𝑢²_𝑅(𝑡) at time 𝑡, respectively, (𝑚=2). We note that the existence of the solution can be obtained using a result by Fleming and Rishel [36]. That is rather straightforward to show that the right-hand sides of (2.1)–(2.6) are bounded by a linear function of the state and control variables and that the integrand of the composed objective functional 𝐼(𝑢) is convex on 𝑈 and is bounded. We now proceed to compute candidates for a Pareto optimal solution. To this end, we apply the Pontryagin’s Maximum Principle [37] and begin by defining the Lagrangian (which is the Hamiltonian augmented with penalty terms for the constraints) to be𝐿𝑥,𝑦,𝑤,𝑧,𝑣,𝑟,𝜆₁,…,𝜆₆,𝑢_𝑃,𝑢_𝑅=𝑤₁𝑥−𝑤₂𝑣−𝑤₃𝑢²_𝑃−𝑤₄𝑢²_𝑅+𝜆₁(𝑠−𝑑𝑥−𝑟𝑥𝑣)+𝜆₂𝑟𝑥𝑣−𝑎𝑦−𝑓𝑦𝑧+𝜆₃𝑐𝑥𝑦𝑤−𝑞𝑦𝑤−𝑏𝑤+𝜆₄𝑞𝑦𝑤−ℎ𝑧+𝜆₅𝑘1−𝑢_𝑃𝑦−𝜏𝑣+𝜆₆𝑟₀−𝑢_𝑅+𝜔₁₁𝑢_𝑃−𝑎₁+𝜔₁₂𝑏₁−𝑢_𝑃+𝜔₂₁𝑢_𝑅−𝑎₂+𝜔₂₂𝑏₂−𝑢_𝑅, (4.2) where 𝜔_𝑖𝑗(𝑡)≥0 are the penalty multipliers satisfying𝜔₁₁𝑢_𝑃−𝑎₁=𝜔₁₂𝑏₁−𝑢_𝑃=0at𝑢^∗_𝑃,𝜔₂₁𝑢_𝑅−𝑎₂=𝜔₂₂𝑏₂−𝑢_𝑅=0at𝑢^∗_𝑅. (4.3) Here, (𝑢^∗_𝑃,𝑢^∗_𝑅) is the optimal control pair yet to be found. Thus, the Maximum Principle gives the existence of adjoint variables satisfying 𝜆̇₁=−𝜕𝐿𝜕𝑥=−𝑤₁+(𝑑+𝑣𝑟)𝜆₁−𝑣𝑟𝜆₂−𝑐𝑦𝑤𝜆₃,𝜆̇₂=−𝜕𝐿=𝜆𝜕𝑦𝑎+𝑓𝑧₂−𝑐𝑥𝑤−𝑞𝑤𝜆₃−𝑞𝑤𝜆₄−𝑘1−𝑢_𝑃𝜆₅,𝜆̇₃=−𝜕𝐿=𝜆𝜕𝑤−𝑐𝑥𝑦+𝑞𝑦+𝑏₃−𝑞𝑦𝜆₄,𝜆̇₄=−𝜕𝐿𝜕𝑧=𝑓𝑦𝜆₂+ℎ𝜆₄,𝜆̇₅=−𝜕𝐿𝜕𝑣=𝑤₂+𝑥𝑟𝜆₁−𝑥𝑟𝜆₂+𝜏𝜆₅,𝜆̇₆=−𝜕𝐿𝜕𝑟=𝑥𝑣𝜆₁−𝑥𝑣𝜆₂, (4.4) where 𝜆_𝑖(𝑡_𝑓)=0,𝑖=1,…,6, are the transversality conditions. The Lagrangian is maximized with respect to 𝑢_𝑃 and 𝑢_𝑅 at the optimal pair (𝑢^∗_𝑃,𝑢^∗_𝑅), so the partial derivatives of the Lagrangian with respect to 𝑢_𝑃 and 𝑢_𝑅 at 𝑢^∗_𝑃 and 𝑢^∗_𝑅 are zero. That is, 𝜕𝐿/𝜕𝑢_𝑃=0 and 𝜕𝐿/𝜕𝑢_𝑅=0.

Since 𝜕𝐿𝜕𝑢_𝑃=−2𝑤₃𝑢_𝑃−𝑘𝑦𝜆₅+𝜔₁₁−𝜔₁₂=0, (4.5)

then we have 𝑢^∗_𝑃(𝑡)=−𝑘𝑦(𝑡)𝜆₅(𝑡)+𝜔₁₁(𝑡)−𝜔₁₂(𝑡)2𝑤₃. (4.6)

To determine an explicit expression for the optimal control without 𝜔₁₁,𝜔₁₂, we consider the following three cases.(i)On the set {𝑡∣𝑎₁<𝑢^∗_𝑃(𝑡)<𝑏₁}, we set 𝜔₁₁(𝑡)=𝜔₁₂(𝑡)=0. Hence, the optimal control is 𝑢^∗_𝑃(𝑡)=−𝑘𝑦(𝑡)𝜆₅(𝑡)2𝑤₃. (4.7)(ii)On the set {𝑡∣𝑢^∗_𝑃(𝑡)=𝑏₁}, we set 𝜔₁₁(𝑡)=0. Hence, the optimal control is 𝑢^∗_𝑃(𝑡)=𝑏₁=−𝑘𝑦(𝑡)𝜆₅(𝑡)−𝜔₁₂(𝑡)2𝑤₃ (4.8)which implies that −𝑘𝑦(𝑡)𝜆₅(𝑡)2𝑤₃≥𝑏₁since𝜔₁₂(𝑡)≥0. (4.9)(iii)On the set {𝑡∣𝑢^∗_𝑃(𝑡)=𝑎₁}, we set 𝜔₁₂(𝑡)=0. Hence, the optimal control is 𝑎₁=𝑢^∗₁(𝑡)=−𝑘𝑦(𝑡)𝜆₅(𝑡)+𝜔₁₁(𝑡)2𝑤₃ (4.10) which implies that−𝑘𝑦(𝑡)𝜆₅(𝑡)2𝑤₃≤𝑎₁since𝜔₁₁(𝑡)≥0. (4.11)

Combining all the three cases in compact form gives𝑢^∗_𝑃(𝑡)=max𝑎₁,min𝑏₁,−𝑘𝑦(𝑡)𝜆₅(𝑡)2𝑤₃. (4.12) Using similar arguments, we also obtain the following expression for the second optimal control function𝑢^∗_𝑅(𝑡)=max𝑎₂,min𝑏₂,−𝜆₆(𝑡)2𝑤₄. (4.13) We point out that the optimality system consists of the state system (2.1)–(2.6) with the initial conditions, the adjoint (or co-state) system (4.4) with the terminal conditions, together with the expressions (4.12) and (4.13) for the control functions. The results of implementing the optimal control policy obtained by this procedure are shown in Section 6.

4.2. STI Solution

From a therapeutic point of view, it may be unsafe to administrate drugs at a dose less than maximum because virus mutations may occur (see, e.g., [21] and references therein for a more exhaustive discussion on this point). Therefore, standard HAART protocols require persistent drugs uptake at maximum value. However, a number of clinical and theoretical studies attempted STI protocols in which periods of therapy at maximum dosages are alternated with periods of treatment suspension [21, 38]. The reasons for these attempts can be found in several side effects of HAART, such as serious hepatic damages and the high cost of the therapy, but also in the evidence that appropriate suspension periods may enhance the immune response of the patient.

In this section, we follow [24] in assuming that the cost of the drug regime varies linearly with the amount of drug administered. Thus, 𝐼₃ and 𝐼₄ are considered as follow:𝐼₃(𝑢)=−^𝑡𝑓_𝑡0𝑢_𝑃(𝑡)𝑑𝑡,𝐼₄(𝑢)=−^𝑡𝑓_𝑡0𝑢_𝑅(𝑡)𝑑𝑡. (4.14)

We now proceed to compute candidates for a Pareto STI solution. To this end, we apply the Pontryagin’s Maximum Principle again and begin by defining the Hamiltonian to be:𝐻𝑥,𝑦,𝑤,𝑧,𝑣,𝑟,𝜆₁,…,𝜆₆,𝑢_𝑃,𝑢_𝑅=𝑤₁𝑥−𝑤₂𝑣−𝑤₃𝑢_𝑃−𝑤₄𝑢_𝑅+𝜆₁(𝑠−𝑑𝑥−𝑟𝑥𝑣)+𝜆₂𝑟𝑥𝑣−𝑎𝑦−𝑓𝑦𝑧+𝜆₃𝑐𝑥𝑦𝑤−𝑞𝑦𝑤−𝑏𝑤+𝜆₄𝑞𝑦𝑤−ℎ𝑧+𝜆₅𝑘1−𝑢_𝑃𝑦−𝜏𝑣+𝜆₆𝑟₀−𝑢_𝑅 (4.15)

and the maximum principle requires that𝐻𝑥^∗,𝑦^∗,𝑤^∗,𝑧^∗,𝑣^∗,𝑟^∗,𝜆^∗₁,…,𝜆^∗₆,𝑢_𝑃,𝑢_𝑅≤𝐻𝑥^∗,𝑦^∗,𝑤^∗,𝑧^∗,𝑣^∗,𝑟^∗,𝜆^∗₁,…,𝜆^∗₆,𝑢^∗_𝑃,𝑢^∗_𝑅,∀𝑢_𝑃,𝑢_𝑅∈𝑈 (4.16) in which the optimal controls and corresponding states and costates satisfying (2.1)–(2.6) and (4.4) are indicated by a star superscript *. From (4.16) and by performing straightforward calculations it is concluded that the optimal pair (𝑢^∗_𝑃,𝑢^∗_𝑅) should maximize the following expression: −𝑤₃𝑢_𝑃−𝑤₄𝑢_𝑅−𝑘𝜆^∗₅𝑢_𝑃𝑦^∗−𝜆^∗₆𝑢_𝑅. (4.17) Therefore,𝑢^∗_𝑅(𝑡)=𝑏₂if𝑆₁𝑢<0,_𝑅,sin𝑔if𝑆₁𝑎=0,₂if𝑆₁𝑢>0,^∗𝑃(𝑡)=𝑏₁if𝑆₂𝑢<0,_𝑃,sin𝑔if𝑆₂𝑎=0,₁if𝑆₂>0. (4.18) where 𝑆₁=𝑤₄+𝜆^∗₆(𝑡) and 𝑆₂=𝑤₃+𝑘𝜆^∗₅(𝑡)𝑦^∗(𝑡) are the switching functions which determine the type of optimal control pair (𝑢^∗_𝑃,𝑢^∗_𝑅). Now we consider an opportunity for the control pair (𝑢^∗_𝑃,𝑢^∗_𝑅) to contain singular arcs. Let us assume that the switching function 𝑆₁ is zero on the singular interval 𝐽₁=[𝑡₁,𝑡₂]. To find a singular control 𝑢_𝑅,sin𝑔, we use the fact that 𝑆₁=0, and 𝑆̇₁=0, on 𝐽₁; hence, 𝜆̇^∗₆(𝑡)=0. Therefore, from the last two equations (4.4), it is concluded that 𝜆^∗₁(𝑡)=𝜆^∗₂(𝑡),𝜆^∗₅(𝑡)𝑤=−²𝜏+𝑤𝜍+²𝜏𝑒^{𝜏(𝑡−𝑡1)}∀𝑡∈𝐽₁, (4.19) where 𝜍=𝜆^∗₅(𝑡₁). Since 𝜆̇^∗₁(𝑡)=𝜆̇^∗₂(𝑡), equating right-hand sides of first two equations (4.4) and then substituting (4.19) into the resulting expression yield𝑥^∗=𝐾₁𝑡,𝑦^∗,𝑤^∗,𝑧^∗,𝜆^∗₁,𝜆^∗₃,𝜆^∗₄,𝑢^∗_𝑃, (4.20) where, 𝐾₁𝑡,𝑦,𝑤,𝑧,𝜆₁,𝜆₃,𝜆₄,𝑢_𝑃𝑤=𝑦+¹𝑐𝑤𝜆₃+𝑞𝑐𝜆1−⁴𝜆₃+𝜆𝑎+𝑓𝑧−𝑑¹𝑐𝑤𝜆₃−𝑘1−𝑢_𝑃𝑐𝑤𝜆₃−𝑤₂𝜏+𝑤𝜍+²𝜏𝑒^{𝜏(𝑡−𝑡1)}. (4.21) Differentiating both sides of (2.1) and then substituting (4.20) into the resulting expression, lead to𝑢_𝑅,sin𝑔(𝑡)=𝐾̈₁𝑡,𝑦^∗,𝑤^∗,𝑧^∗,𝜆^∗₁,𝜆^∗₃,𝜆^∗₄,𝑢^∗_𝑃−𝐾₂𝑟^∗,𝑥^∗,𝑦^∗,𝑣^∗,𝑢^∗_𝑃𝑥^∗𝑣^∗,𝑡∈𝐽₁, (4.22)

where𝐾₂𝑟,𝑥,𝑦,𝑣,𝑢_𝑃=(𝑠−𝑑𝑥−𝑟𝑥𝑣)(−𝑑−𝑟𝑣)−𝑟₀−𝜏𝑟𝑥𝑣−𝑘𝑟𝑥𝑦1−𝑢_𝑃. (4.23) Obviously, 𝐾̈₁ can be calculated using the Chain rule. Now, we find the value of singular control 𝑢_𝑃,sin𝑔. With respect to this fact that 𝑑^ℓ𝑆₂/𝑑𝑡^ℓ=0,ℓ=0,1,2, and from second last equation (4.4) and (2.2) it is easy to see that during the corresponding singular interval, say 𝐽₂, 𝑦̈^∗=2𝑤₃𝑟^∗𝑥^∗𝑣^∗−𝑎𝑦^∗−𝑓𝑦^∗𝑧^∗2−2𝑘𝑤₂+𝑥^∗𝑟^∗𝜆^∗₁−𝑥^∗𝑟^∗𝜆^∗₂+𝜏𝜆^∗₅𝑟^∗𝑥^∗𝑣^∗−𝑎𝑦^∗−𝑓𝑦^∗𝑧^∗𝑦^∗2𝑤₃𝑦^∗+𝑘𝑦^∗2𝜆^∗₅. (4.24) Moreover, differentiating both sides of (2.2) and then substituting (4.24) into the resulting expression yield 𝑢_𝑃,sin𝑔(𝑡)=𝐾₃𝑥^∗,𝑦^∗,𝑤^∗,𝑧^∗,𝑣^∗,𝑟^∗,𝑢^∗_𝑅𝑘𝑟^∗𝑥^∗𝑦^∗−2𝑟^∗𝑥^∗𝑣^∗−𝑎𝑦^∗−𝑓𝑦^∗𝑧^∗𝑟^∗𝑥^∗𝑤₃+𝑘𝜆^∗₅𝑤₃𝑟^∗𝑥^∗𝑣^∗−𝑎𝑦^∗−𝑓𝑦^∗𝑧^∗𝑘𝑦^∗2−𝜆^∗₅,𝑡∈𝐽₂, (4.25) where𝐾₃𝑥,𝑦,𝑤,𝑧,𝑣,𝑟,𝑢_𝑅=𝑟₀−𝑢_𝑅𝑥𝑣+𝑠𝑟𝑣−(𝑎+𝑑+𝜏)𝑥𝑟𝑣−𝑟²𝑣²𝑥+𝑘𝑟𝑥𝑦+𝑎²𝑦+2𝑎𝑓+𝑓ℎ𝑦𝑧−𝑓𝑞𝑦²𝑤−𝑓𝑟𝑥𝑣𝑧+𝑓²𝑦𝑧². (4.26)

Note that the sets 𝐽₁ and 𝐽₂ are considered such that (𝑢_𝑃,sin𝑔,𝑢_𝑅,sin𝑔)∈𝑈. Moreover, for the control pair (𝑢^∗_𝑃,𝑢^∗_𝑅) to be optimal, the Kelley’s condition [39] must be satisfied (−1)^𝛾𝜕𝜕𝑢_𝑅𝑑^2𝛾𝑑𝑡^2𝛾𝑆₁≥0,(−1)^𝜁𝜕𝜕𝑢_𝑃𝑑^2𝜁𝑑𝑡^2𝜁𝑆₂≥0, (4.27) where 𝛾 and 𝜁 are known as the order of singularity.

The above discussion shows that using linear cost functionals, while ignoring singular arcs, leads to an STI solution. Unfortunately, solving the corresponding optimality system is not analytically possible and the gradient method [40] fails to converge; therefore, we resort to model predictive control- (MPC-) based approach. MPC is a control technology widely used in many areas, especially in the process industries [41], for systems with a large number of controlled and manipulated variables, which interact significantly. In general, feedback control technologies, and MPC in particular, have started to gain significant attention in the biomedical area [12, 19–21]. A thorough overview of the history of MPC and its various incarnations can be found in [42]. In order to give a formal description of the proposed MPC algorithm, some preliminary definitions are necessary. Let Δ be the sampling time, we denote with 𝑋(𝑖) the state of the system at time 𝑡_𝑖=𝑖Δ, that is,𝑋(𝑖)=𝑥𝑡_𝑖𝑦𝑡_𝑖𝑤𝑡_𝑖𝑧𝑡_𝑖𝑣𝑡_𝑖𝑟𝑡_𝑖=𝑥(𝑖)𝑦(𝑖)𝑤(𝑖)𝑧(𝑖)𝑣(𝑖)𝑟(𝑖). (4.28) Next, let 𝑈(𝑖)=[𝑢_𝑃(𝑡_𝑖)𝑢_𝑅(𝑡_𝑖)]=[𝑢_𝑃(𝑖)𝑢_𝑅(𝑖)] denote the vector of PI and RTI drugs taken in the time interval [𝑖Δ,(𝑖+1)Δ], and rewrite the system model (2.1)–(2.6) in the integrated form 𝑋(𝑖+1)=𝐹(𝑋(𝑖),𝑈(𝑖)), where 𝐹(⋅) is a vector-value function obtained by numerical integration of (2.1)–(2.6) over the sampling Δ, assuming a constant drug uptake during the sampling time. With these definitions, it is now possible to state the following finite-horizon optimal control problem (FHOCP) to be solved at each discrete time 𝑘:Maximize_𝑈𝑘^{𝑁+𝑘−1}_𝑖=𝑘𝑤₁𝑥(𝑖)−𝑤₂𝑣(𝑖)−𝑤₃𝑢_𝑃(𝑖)−𝑤₄𝑢_𝑅(𝑖)(𝑘)𝑋𝑗𝑗𝑈𝑗∈𝑎subjectto∶𝑋isgiven,𝑋𝑗+1=𝐹,𝑈,𝑗=𝑘,…,𝑁+𝑘−1,₁,𝑏₁×𝑎₂,𝑏₂,𝑗=𝑘,…,𝑁+𝑘−1, (4.29) in which 𝑁 is an integer(control horizon), and 𝑈_𝑘=[𝑈(𝑘)⋯𝑈(𝑁+𝑘−1)] is the decision variable. Let 𝑈^∗_𝑘 be a maximizing control input sequence. The first input in the optimal sequence, that is, 𝑈^∗(𝑘), is injected into the system, and the entire optimization is repeated at subsequent control intervals.

5. Piecewise Constant Solution Using Fuzzy Aggregation and Embedding Method

In this section, as in continuous solutions, we set 𝑚=2. With respect to the nature of controls 𝑢_𝑃 and 𝑢_𝑅, that indicate the PI and RTI drug efficacies as a function of time, piecewise constant control is more practical from the clinical point of view. In this section, we are going to propose such kind of solutions for problem (3.3).

Setting, 𝜉=(𝑥,𝑦,𝑤,𝑧,𝑣,𝑟) and 𝑢=(𝑢_𝑃,𝑢_𝑅), the system of differential equations (2.1)–(2.6) can be represented in a generalized form as𝜉̇(𝑡)=𝑔(𝑡,𝜉(𝑡),𝑢(𝑡))=𝑠−𝑑𝜉₁−𝜉₁𝜉₅𝜉₆𝜉₁𝜉₅𝜉₆−𝑎𝜉₂−𝑓𝜉₂𝜉₄𝑐𝜉₁𝜉₂𝜉₃−𝑞𝜉₂𝜉₃−𝑏𝜉₃𝑞𝜉₂𝜉₃−ℎ𝜉₄𝑘1−𝑢_𝑃𝜉₂−𝜏𝜉₅𝑟₀−𝑢_𝑅,𝜉(0)=𝜉₀. (5.1)In the absence of importance on the objectives and without knowledge of the possible level of attainment for those objectives, finding a Pareto optimal solution that the best satisfies the decision maker can be viewed as a fuzzy problem. In this situation, for each of the objective functionals, 𝐼_𝑖,𝑖=1,2,3,4, of problem (3.3) we associate a linear membership function 𝜇_𝑖(𝐼_𝑖(𝑢)) as 𝜇_𝑖𝐼_𝑖(𝑢)=0,𝐼_𝑖(𝑢)<𝑚_𝑖,𝐼_𝑖(𝑢)−𝑚_𝑖𝑀_𝑖−𝑚_𝑖,𝑚_𝑖≤𝐼_𝑖(𝑢)≤𝑀_𝑖,1,𝐼_𝑖(𝑢)>𝑀_𝑖, (5.2)where 𝑚_𝑖 or 𝑀_𝑖 denotes the value of the objective functional such that the degree of membership function is 0 or 1, respectively, and it is depicted in Figure 1. Following the principle of the fuzzy decision of Bellman and Zadeh [43], the MOOCP (3.3) can be interpreted as the following problem: Maximize_𝑢∈𝑈min_𝑖=1,…,4𝜇_𝑖𝐼_𝑖(𝑢). (5.3)Let 𝜌 be a given small positive number. Introducing an auxiliary variable 𝜆 and by adding the augmented term 𝜌∑⁴_𝑘=1𝐼_𝑖(𝑢) to it, a solution to problem (3.3), can be obtained by solving the following problem:max_{𝜆,𝑢∈𝑈}𝜆+𝜌⁴_𝑖=1𝐼_𝑖(𝑢)subjectto𝑚_𝑖≤𝐼_𝑖(𝑢)−𝜆𝑀_𝑖−𝑚_𝑖̇()𝑡,𝑖=1,2,3,4𝜉=𝑔𝑡,𝜉,𝑢,𝜉₀=𝜉₀. (5.4)(5.5)(5.6)

Figure 1 

Linear membership function for the ith fuzzy objective.

Theorem 5.1. If (𝜆^∗,𝑢^∗) is the optimal solution of the problem (5.4)–(5.6), then 𝑢^∗ is a Pareto optimal solution of (3.3).

Proof. If 𝑢^∗ is not a Pareto optimal solution for the problem (3.3), then there exists 𝑢∈𝑈 such that 𝐼_𝑖(𝑢^∗)≤𝐼_𝑖(𝑢) for all 𝑖∈{1,2,3,4}, and 𝐼_𝑗(𝑢^∗)<𝐼_𝑗(𝑢) for some 𝑗∈{1,2,3,4}. Obviously, 𝜆^∗=min_𝑖=1,…,4{𝜇_𝑖(𝐼_𝑖(𝑢^∗))}≤min_𝑖=1,…,4{𝜇_𝑖(𝐼_𝑖(𝑢))}=𝜆 and ∑⁴_𝑖=1𝐼_𝑖(𝑢^∗)<∑⁴_𝑖=1𝐼_𝑖(𝑢). Therefore, 𝜆^∗+𝜌∑⁴_𝑖=1𝐼_𝑖(𝑢^∗)<𝜆+𝜌∑⁴_𝑖=1𝐼_𝑖(𝑢). This contradicts the assumption that (𝜆^∗,𝑢^∗) is the optimal solution for problem (5.4)–(5.6).

Using the measure theory for solving optimal control problems based on the idea of Young [44], which was applied for the first time by Wilson and Rubio [45], has been theoretically established by Rubio in [46]. In the rest of this section, we follow their approach and we transfer the problem (5.4)-(5.6) into a modified problem in an appropriate space that the existence of the optimal solution (𝜆^∗,𝑢^∗) is guaranteed.

5.1. Transformation to Functional Space

We assume that state variable 𝜉(⋅) and control input 𝑢(⋅) get their values in the compact sets 𝐴=𝐴₁×⋯×𝐴₆⊂ℜ⁶ and 𝑈=𝑈₁×𝑈₂⊂ℜ², respectively. Set 𝐽=[𝑡₀,𝑡_𝑓].

Definition 5.2. A pair 𝑝=[𝜉,𝑢] 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 (5.5) and (5.6) on 𝐽⁰, that is, the interior set of 𝐽.We assume that the set of all admissible pairs is nonempty and denote it by 𝑊. Let 𝑝 be an admissible pair and 𝐵 an open ball in ℜ⁶ containing 𝐽×𝐴, and 𝐶′(𝐵) the space of all real-valued continuous differentiable functions on it. Let 𝜑∈𝐶′(𝐵) and define 𝜑^𝑔 as follows: 𝜑^𝑔((𝑡)(𝑡=((𝑡𝑡,𝜉,𝑢))𝑑𝜑𝑡,𝜉))𝑑𝑡=⁶_𝑗=1𝜕𝜑(𝑡,𝜉(𝑡))𝜕𝜉_𝑗𝑔_𝑗((𝑡)(𝑡+((𝑡𝑡,𝜉,𝑢))𝜕𝜑𝑡,𝜉))𝜕𝑡 (5.7) for each [𝑡,𝜉(𝑡),𝑢(𝑡)]∈Ω, where Ω=𝐽×𝐴×𝑈. The function 𝜑^𝑔 is in the space 𝐶(Ω), the set of all continuous functions on the compact set Ω. Since 𝑝=[𝜉,𝑢] is an admissible pair, we have ^𝑡𝑓_𝑡0𝜑^𝑔(𝑡,𝜉(𝑡),𝑢(𝑡))𝑑𝑡=𝜑𝑡_𝑓,𝜉𝑡_𝑓−𝜑𝑡₀,𝜉𝑡₀=Δ𝜑, (5.8) for all 𝜑∈𝐶′(𝐵). Let 𝐷(𝐽⁰) be the space of infinitely differentiable all real-valued functions with compact support in 𝐽⁰. Define 𝜓^𝑛(𝑡,𝜉(𝑡),𝑢(𝑡))=𝜉_𝑛(𝑡)𝜓^′(𝑡)+𝑔_𝑛(𝑡,𝜉(𝑡),𝑢(𝑡))𝜓(𝑡),𝑛=1,…,6,∀𝜓∈𝐷𝐽⁰. (5.9) Assume 𝑝=[𝜉,𝑢] be an admissible pair. Since the function 𝜓(⋅) has compact support in 𝐽⁰, so, 𝜓(𝑡₀)=𝜓(𝑡_𝑓)=0. Thus, for 𝑛=1,…,6, and for all 𝜓∈𝐷(𝐽⁰), by integrating both sides of (5.9) and using integration by parts, we have ^𝑡𝑓_𝑡0𝜓^𝑛(𝑡,𝜉(𝑡),𝑢(𝑡))𝑑𝑡=^𝑡𝑓_𝑡0𝜉_𝑛(𝑡)𝜓^′(𝑡)𝑑𝑡+^𝑡𝑓_𝑡0𝑔_𝑛(𝑡,𝜉(𝑡),𝑢(𝑡))𝜓(𝑡)𝑑𝑡=𝜉_𝑛(𝑡)𝜓(𝑡)∣^𝑡𝑓_𝑡0−^𝑡𝑓_𝑡0𝜉̇_𝑛(𝑡)−𝑔_𝑛(𝑡,𝜉(𝑡),𝑢(𝑡))𝜓(𝑡)𝑑𝑡=0. (5.10) Also by choosing the functions which are dependent only on time, we have ^𝑡𝑓_𝑡0𝜗(𝑡,𝜉(𝑡),𝑢(𝑡))𝑑𝑡=𝑎_𝜗,∀𝜗∈𝐶¹(Ω), (5.11) where 𝐶¹(Ω) is the space of all functions in 𝐶(Ω) that depend only on time and 𝑎_𝜗 is the integral of 𝜗 on 𝐽. Equations (5.8), (5.10), and (5.11) are weak forms of (5.6). Now, we consider the following positive linear functional on 𝐶(Ω)Γ_𝑝∶𝐹⟶_𝐽𝐹(𝑡,𝑥(𝑡),𝑢(𝑡))𝑑𝑡,∀𝐹∈𝐶(Ω). (5.12)

Proposition 5.3. Transformation 𝑝→Γ_𝑝 of admissible pairs in 𝑊 into the linear mappings Γ_𝑝 defined in (5.12) is an injection.

Proof. See [46].

Thus, the problem (5.4)–(5.6) is converted to the following optimization problem in functional space: maximize_{𝜆,𝑝∈𝑊}𝜆+Γ_𝑝𝑓₀subjectto𝑚_𝑖≤Γ_𝑝𝑓_𝑖−𝜆𝑀_𝑖−𝑚_𝑖Γ𝑖=1,2,3,4,_𝑝𝜑^𝑔(𝐵)Γ=Δ𝜑,𝜑∈𝐶′,_𝑝𝜓^𝑛=0,𝑛=1,…,6,𝜓∈𝐷𝐽⁰Γ_𝑝,(𝜗)=𝑎_𝜗,𝜗∈𝐶¹(Ω), (5.13)(5.14)(5.15)(5.16)(5.17) where 𝑓₀(𝑡,𝜉,𝑢)=𝜌(𝜉₁−𝜉₅−𝑢²_𝑃−𝑢²_𝑅), 𝑓₁(𝑡,𝜉,𝑢)=𝜉₁, 𝑓₂(𝑡,𝜉,𝑢)=−𝜉₅, 𝑓₃(𝑡,𝜉,𝑢)=−𝑢²_𝑃, and 𝑓₄(𝑡,𝜉,𝑢)=−𝑢²_𝑅.

5.2. Transformation to Measure Space

Let 𝑀⁺(Ω) denote the space of all positive Radon measures on Ω. By the Riesz representation theorem [46], there exists a unique positive Radon measure𝜇 on Ω such that Γ_𝑝(𝐹)=_𝐽𝐹(𝑡,𝜉(𝑡),𝑢(𝑡))𝑑𝑡=_Ω𝐹(𝑡,𝜉,𝑢)𝑑𝜇≡𝜇(𝐹),𝐹∈𝐶(Ω). (5.18)So, we may change the space of optimization problem to measure space. In other words, the optimization problem in functional space (5.13)–(5.17) can be replaced by the following new problem in measure space: maximize_{𝜆,𝜇∈𝑀+(Ω)}𝜆+𝜇𝑓₀subjectto𝑚_𝑖≤𝜇𝑓_𝑖−𝜆𝑀_𝑖−𝑚_𝑖𝜇𝜑,𝑖=1,2,3,4,_𝑔(𝐵)𝜇𝜓=Δ𝜑,𝜑∈𝐶′,_𝑛=0,𝑛=1,…,6,𝜓∈𝐷𝐽⁰𝜇(𝜗),=𝑎_𝜗,𝜗∈𝐶¹(Ω). (5.19)(5.20)(5.21)(5.22)(5.23) Obviously, 𝜆 takes its values in the compact set [0,1]⊆ℜ. We shall consider the maximization of (5.21) over the set 𝑄 of all pairs (𝜇,𝜆)∈𝑀⁺(Ω)×[0,1] satisfying (5.20)–(5.23). The main advantages of considering this measure theoretic form of the problem is the existence of optimal pair (𝜇^∗,𝜆^∗)∈𝑄, 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 𝜎(𝜇,𝜆)=𝜆+𝜇(𝑓₀). The following theorem guarantees the existence of an optimal solution.

Theorem 5.4. Function 𝜎 attains its maximum (𝜇^∗,𝜆^∗) on the set𝑄.

Proof. The set 𝑄can be written as 𝑄=𝑄₁∩𝑄₂, (5.24) where 𝑄₁=𝜇,𝜆∈𝑀⁺(Ω)×[0,1]∶𝑚_𝑖≤𝜇𝑓_𝑖−𝜆𝑀_𝑖−𝑚_𝑖,𝑄,𝑖=1,2,3,4₂=𝜇,𝜆∈𝑀⁺(Ω)×[0,1]∶𝜇satisfiesin(5.21)-(5.23). (5.25) We will show in the next section that (5.22) and (5.23) are special cases of (5.21). Therefore, the set 𝑄₂ can be written as 𝑄₂=Π×[0,1], where Π=_{𝜑∈𝐶′(𝐵)}𝜇∈𝑀⁺(Ω)∶𝜇𝜑^𝑔=Δ𝜑. (5.26) It is well known that the set {𝜇∈𝑀⁺(Ω)∶𝜇(1)=𝑡_𝑓−𝑡₀} is compact in weak*-topology [46]. Furthermore, for 𝜑∈𝐶′(𝐵), the set Π as intersection of inverse image of the 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 Π. Hence, the compactness of [0,1] implies the compactness of 𝑄₂ with the product topology. Besides, for 𝑖=1,2,3,4, the set 𝑄₁ as intersection of the inverse image of the closed sets [𝑚_𝑖,∞) under the continuous functions (𝜇,𝜆)→𝜇(𝑓_𝑖)−𝜆(𝑀_𝑖−𝑚_𝑖) is also closed. Therefore, the set 𝑄 as a closed subset of the compact set 𝑄₂ is compact also. Since the function 𝜎 mapping the compact set 𝑄 on the real line is continuous so it has a maximum on the compact set 𝑄.

5.3. Approximation

Since 𝑀⁺(Ω) is an infinite-dimensional space, the problem (5.19)–(5.23) 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 (5.21)–(5.23) be satisfied. Let {𝜑_𝑖∶𝑖=1,2,…}, {𝜓_𝑗∶𝑗=1,2,…}, and {𝜗_𝑠∶𝑠=1,2,…} be the sets of total functions, respectively, in 𝐶′(𝐵), 𝐷(𝐽⁰) and 𝐶¹(Ω).

Remark 5.5. From (5.7) and (5.9), it can be seen that (5.22) and (5.23) are also achieved from (5.21) by setting 𝜑(𝑡,𝜉(𝑡))=𝜉_𝑛(𝑡)𝜓(𝑡) and 𝜑(𝑡,𝜉(𝑡))=∫^𝑡₀𝜗(𝜏)𝑑𝜏, respectively.

Now, the first approximation will be completed by choosing finite number of total functions and using following propositions.

Proposition 5.6. Consider the linear program consisting of maximizing the function 𝜎 over the set 𝑄(𝐾) of pairs (𝜇,𝜆) in 𝑀⁺(Ω)×[0,1] satisfying 𝑚_𝑖≤𝜇𝑓_𝑖−𝜆𝑀_𝑖−𝑚_𝑖𝜇𝜑i=1,2,3,4,_𝑔𝑗=Δ𝜑_𝑗,𝑗=1,…,𝐾. (5.27) Then, 𝜂_𝐾≡max_𝑄(𝐾)𝜎 tends to 𝜂^∗≡max_𝑄𝜎 as 𝐾→∞.

Proof. We have 𝑄(1)⊇𝑄(2)⊇⋯⊇𝑄(𝑀)⊇⋯⊇𝑄 and hence, 𝜂₁≥𝜂₂≥⋯≥𝜂_𝑀≥⋯≥𝜂^∗. Therefore, {𝜂_𝑗} is nonincreasing and bounded sequence then converges to a number 𝜁 such that 𝜁≥𝜂^∗. We show that, 𝜁=𝜂^∗. Set 𝑅≡⋂^∞_𝐾=1𝑄(𝐾). Then, 𝑅⊇𝑄 and 𝜁≡max_𝑅𝜎. It is sufficient to show that 𝑅⊆𝑄. Assume (𝜇,𝜆)∈𝑅 and 𝜑∈𝐶^′(𝐵). Since the Linear combinations of the functions {𝜑_𝑗,𝑗=1,2,…} are uniformly dense in 𝐶^′(𝐵), there is the sequence {𝜑_𝑘}∈span{𝜑_𝑗,𝑗=1,2,…} such that 𝜑_𝑘 tends to 𝜑 uniformly as 𝑘→∞. Hence, 𝑆₁, 𝑆₂, and 𝑆₃ tend to zero as 𝑘→∞, where, 𝑆₁=sup𝜑_𝜉(𝑡,𝜉)−𝜑_𝑘𝜉(𝑡,𝜉),𝑆₂=sup𝜑_𝑡(𝑡,𝜉)−𝜑_𝑘𝑡(𝑡,𝜉),𝑆₃=sup𝜑(𝑡,𝜉)−𝜑_𝑘(𝑡,𝜉). (5.28) We have 𝜇∈𝑅 and functional 𝑓→𝜇(𝑓) is linear. Therefore, 𝜇(𝜑_𝑘^𝑔)=Δ𝜑_𝑘 and 𝜇𝜑^𝑔−Δ𝜑=𝜇𝜑^𝑔−Δ𝜑−𝜇𝜑_𝑘^𝑔+Δ𝜑_𝑘=_Ω𝜑_𝜉(𝑡,𝜉)−𝜑_𝑘𝜉(𝑡,𝜉)𝑔(𝑡,𝜉,𝑢)+𝜑_𝑡(𝑡,𝜉)−𝜑_𝑘𝑡(𝑡,𝜉)𝑑𝜇−Δ𝜑−Δ𝜑_𝑘≤𝑆_1Ω𝑔(𝑡,𝜉,𝑢)𝑑𝜇+𝑆_2Ω𝑑𝜇+2𝑆₃. (5.29) Since the right-hand side of the above inequality tends to zero as 𝑘→∞, while left-hand side is independent of 𝑘, therefore, 𝜇(𝜑^𝑔)=Δ𝜑. Thus 𝑅⊆𝑄 and 𝜂^∗≥𝜁 which implies 𝜁=𝜂^∗.

Proposition 5.7. An optimal pair (𝜇^∗,𝜆^∗) in the set 𝑄(𝐾) at which the function 𝜎 attains its maximum can be found where the measure 𝜇^∗ has the form 𝜇^∗=^𝐾+4_𝑗=1𝛼^∗_𝑗𝛿𝑧^∗_𝑗, (5.30) where 𝛼^∗_𝑗≥0,𝑧^∗_𝑗∈Ω, and 𝛿(𝑧) is unitary atomic measure with the support being the singleton set {𝑧^∗_𝑗}, characterized by 𝛿(𝑧)(𝐹)=𝐹(𝑧),𝑧∈Ω.

Proof. Assume that (𝜇^∗,𝜆^∗) is an optimal solution of maximizing the function 𝜎 over the set 𝑄(𝐾). Denote the set of all measures 𝜇∈𝑀⁺(Ω) satisfying in (5.27) for a fixed 𝜆=𝜆^∗, by 𝑄_𝜆^∗. The set 𝑄_𝜆^∗ is a weak*-compact subset of 𝑀⁺(Ω) [46]. Therefore, from Theorem A.5 of [46], 𝑄_𝜆^∗ has at least one extreme point in the form of (5.30).

Therefore, our attention is restricted to finding a measure in the form of (5.30) which maximizes the function 𝜎 and satisfies in (5.20) and K number of the constraints (5.21)–(5.23). Thus, by choosing the functions 𝜑_𝑖, 𝑖=1,2,…,𝐾₁, 𝜓_𝑘,𝑘=1,…,𝐾₂, and 𝜗_𝑠,𝑠=1,…,𝑆, the infinite-dimensional problem (5.19)–(5.23) is approximated by the following finite-dimensional nonlinear programming (NLP) problem: maximize_{𝜆,𝛼𝑗≥0,𝑧𝑗∈Ω}𝜆+^𝐾+4_𝑗=1𝛼_𝑗𝑓₀𝑧_𝑗subjectto𝑚_𝑖≤^𝐾+4_𝑗=1𝛼_𝑗𝑓_𝑖𝑧_𝑗−𝜆𝑀_𝑖−𝑚_𝑖,𝑖=1,2,3,4,_𝐾+4𝑗=1𝛼_𝑗𝜑^𝑔_𝑖𝑧_𝑗=Δ𝜑_𝑖,𝑖=1,…,𝐾₁,_𝐾+4𝑗=1𝛼_𝑗𝜓^𝑛_𝑘𝑧_𝑗=0,𝑘=1,…,𝐾₂,𝑛=1,…,6,_𝐾+4𝑗=1𝛼_𝑗𝜗_𝑠𝑧_𝑗=𝑎_𝜗𝑠,𝑠=1,…,𝑆, (5.31) where 𝐾=𝐾₁+6𝐾₂+𝑆. Clearly, (5.31) is an NLP problem with 2(𝐾+4)+1 unknowns 𝜆, 𝛼_𝑗,𝑧_𝑗,𝑗=1,…,𝐾+4. One is interested in LP problem. The following proposition enables us to approximate the NLP problem (5.31) by a finite-dimensional LP problem.

Proposition 5.8. Let Ω_𝑁={𝑦₁,𝑦₂,…,𝑦_𝑁} be a countable dense subset of Ω, where 𝑁 is a sufficiently large number. Given 𝜀>0, a measure 𝑣∈𝑀⁺(Ω) can be found such that 𝑣𝑓_𝑖−𝜇^∗𝑓_𝑖𝑣𝜑≤𝜀,𝑖=0,1,2,3,4,_𝑔𝑗−𝜇^∗𝜑^𝑔_𝑗≤𝜀,𝑗=1,…,𝐾₁,𝑣𝜓^𝑛_𝑘−𝜇^∗𝜓^𝑛_𝑘≤𝜀,𝑘=1,…,𝐾₂𝑣𝜗,𝑛=1,…,6,_𝑠−𝜇^∗𝜗_𝑠≤𝜀,𝑠=1,…,𝑆, (5.32) where the measure 𝑣 has the form 𝑣=^𝐾+4_𝑗=1𝛼^∗_𝑗𝛿𝑦_𝑗, (5.33) and the coefficients 𝛼^∗_𝑗,𝑗=1,…,𝐾+4, are the same as optimal measure (5.30), and 𝑦_𝑗∈Ω_𝑁,𝑗=1,…,𝐾+4.

Proof. We rename the functions 𝑓_𝑖s, 𝜑^𝑔_𝑗s, 𝜓^𝑛_𝑘s, and 𝜗_𝑠s sequentially as ℎ_𝑗,𝑗=1,2,…,𝐾+5. Then, for 𝑗=1,…,𝐾+5, 𝜇^∗−𝑣ℎ_𝑗=^𝐾+4_𝑖=1𝛼^∗_𝑖ℎ_𝑗𝑧^∗_𝑖−ℎ_𝑗𝑦_𝑖≤^𝐾+4_𝑖=1𝛼^∗_𝑖max_𝑖,𝑗ℎ_𝑗𝑧^∗_𝑖−ℎ_𝑗𝑦_𝑖. (5.34) The functions ℎ_𝑗,𝑗=1,2,…,𝐾+5, are continues. Therefore, max_𝑖,𝑗|ℎ_𝑗(𝑧^∗_𝑖)−ℎ_𝑗(𝑦_𝑖)| can be made less than 𝜀/∑^𝐾+4_𝑗=1𝛼^∗_𝑗 by choosing 𝑦_𝑖,𝑖=1,2,…,𝐾+4, sufficiently near 𝑧^∗_𝑖.

For constructing a suitable set Ω_𝑁, 𝐽 is divided to 𝑆 subintervals as follows: 𝐽_𝑠=𝑡₀+(𝑠−1)Δ,𝑡₀+𝑠Δ,𝑠=1,2,…,𝑆, (5.35) where Δ=(𝑡_𝑓−𝑡₀)/𝑆. Moreover, the intervals 𝐴_𝑖(𝑖=1,…,6) and 𝑈_𝑗(𝑗=1,2) are divided, respectively, into 𝑛_𝑖 and 𝑠_𝑗 subintervals. So, the set Ω is divided into 𝑁=𝑆𝑛₁𝑛₂𝑛₃𝑛₄𝑛₅𝑛₆𝑠₁𝑠₂ cells. One point is chosen from each cell. In this way, we will have a grid of points, which are numbered sequentially as 𝑦_𝑗=(𝑡_𝑗,𝜉_1𝑗,…,𝜉_6𝑗,𝑢_𝑃𝑗,𝑢_𝑅𝑗), 𝑗=1,…,𝑁.

The desired numbers 𝑚_𝑖 and 𝑀_𝑖 in definition of the membership functions 𝜇_𝑖(𝐼_𝑖(𝑢)),𝑖=1,2,3,4 can be set𝑚_𝑖=min_𝑢∈𝑈^𝑡𝑓_𝑡0𝑓_𝑖(𝑡,𝜉,𝑢)𝑑𝑡,𝑀_𝑖=max_𝑢∈𝑈^𝑡𝑓_𝑡0𝑓_𝑖(𝑡,𝜉,𝑢)𝑑𝑡. (5.36)

Remark 5.9. With respect to this fact that the maximum level of healthy cells is attained with full HAART, we have 𝑀₁=𝐼₁(𝑏₁,𝑏₂). Similarly, the other numbers 𝑚_𝑖 and 𝑀_𝑖 can be calculated without explicitly solving auxiliary optimal control problems. The values of these variables are shown in Table 3.

Therefore, according to (5.33) the NLP problem (5.31) is converted to the following LP problem:maximize_{𝜆,𝛼𝑗≥0}𝜆+^𝑁_𝑗=1𝛼_𝑗𝑓₀𝑦_𝑗subjectto𝑚_𝑖≤^𝑁_𝑗=1𝛼_𝑗𝑓_𝑖𝑦_𝑗−𝜆𝑀_𝑖−𝑚_𝑖𝑖=1,2,3,4,_𝑁𝑗=1𝛼_𝑗𝜑^𝑔_𝑖𝑦_𝑗=Δ𝜑_𝑖,𝑖=1,…,𝐾₁,_𝑁𝑗=1𝛼_𝑗𝜓^𝑛_𝑘𝑦_𝑗=0,𝑘=1,…,𝐾₂,𝑛=1,…,6,_𝑁𝑗=1𝛼_𝑗𝜗_𝑠𝑦_𝑗=𝑎_𝜗𝑠,𝑠=1,…,𝑆. (5.37)

Here, we discuss suitable total functions 𝜑_𝑖, 𝑖=1,…,𝐾₁,𝜓_𝑘,𝑘=1,…,𝐾₂, and 𝜗_𝑠,𝑠=1,…,𝑆. The functions 𝜑_𝑖s can be taken to be monomials of 𝑡 and the components of the vector 𝜉 as follows:𝑡^𝑖𝜉_𝑗,𝑖∈{0,1,2,…},𝑗∈{1,…,6}. (5.38)In addition, we choose some functions with compact support in the following form [46]: 𝜓_2𝑟−1(𝑡)=sin2𝜋𝑟𝑡−𝑡₀𝑡_𝑓−𝑡₀,𝜓_2𝑟(𝑡)=1−cos2𝜋𝑟𝑡−𝑡₀𝑡_𝑓−𝑡₀,𝑟=1,2,…. (5.39)

Finally, the following functions are considered that are dependent on 𝑡 only:𝜗_𝑠(𝑡)=1𝑡∈𝐽_𝑠,0otherwise, (5.40)

where 𝐽_𝑠,𝑠=1,…,𝑆, are given by (5.35). These functions are used to construct the approximate piecewise constant controls [46]. Of course, we need only to construct the control function 𝑢(⋅), since 𝜉(⋅) can be obtained by solving the ODEs (2.1)–(2.6). By using simplex method, a nonzero optimal solution 𝛼^∗_𝑖1,𝛼^∗_𝑖2,…,𝛼^∗_𝑖𝑘, 𝑖₁<𝑖₂<⋯<𝑖_𝑘, of the LP problem (5.37) can be found where 𝑘 cannot exceed the number of constraints, that is, 𝑘≤𝐾+4. Setting 𝛼^∗_𝑖0=𝑡₀, a piecewise constant control function 𝑢(⋅) approximating the optimal solution is constructed based on these nonzero coefficients as follows [46]:𝑢(𝑡)=𝑢_𝑃𝑖𝑗,𝑢_𝑅𝑖𝑗𝑡∈^𝑗−1_ℎ=0𝛼^∗_𝑖ℎ,^𝑗_ℎ=0𝛼^∗_𝑖ℎ0otherwise,𝑗=1,2,…,𝑘, (5.41)

where 𝑢_𝑃𝑖𝑗and 𝑢_𝑅𝑖𝑗are, respectively, the 8th and the 9th components of grid point 𝑦_𝑖𝑗.

6. Numerical Results

In this section we show the response of the HIV model to the continuous, STI and piecewise constant controls presented in Sections 4 and 5.

6.1. Continuous Solutions

As we discussed in Section 4.1, in the case of continuous solutions the optimality systems, is a two-point boundary value problem because of the forward-in-time nature of the state system (2.1)–(2.6) with the initial conditions and the backward-in-time nature of the adjoint system (4.4) with the terminal conditions. We consider a gradient method for the solution of the optimal control problem. The algorithm proceeds as follows.(i)Choose an initial guess of control.(ii)Solve the state system forward in time with initial conditions using an initial guess of control.(iii)Solve the adjoint system backward in time with terminal conditions.(iv)Update the controls in each iteration using the optimality conditions (4.12) and (4.13).(v)Continue the iterations until convergence is achieved.

For further discussion of the gradient method, we refer the interested reader to [40]. We use the parameters for solving the optimality system from Table 1. Antiretroviral therapy is initiated at 𝑡₀, the time at which the CD4+ T-cell count falls below 350 cells/μL, and treatment was simulated for 750 days. Furthermore, control variables 𝑢_𝑃 and 𝑢_𝑅 are bounded by 𝑎₁=0, 𝑏₁=0.7, 𝑎₂=0, and 𝑏₂=9×10⁻¹⁰, and the following initial condition is used (see [7]):𝑥(0)=10³cells𝜇L⁻¹,𝑣(0)=104copiesmL⁻¹,𝑦(0)=0cells𝜇L⁻¹,𝑤(0)=10⁻³cellsmL⁻¹,𝑧(0)=10⁻⁷cells𝜇L⁻¹,𝑟(0)=2×10⁻⁷mLcopies⁻¹day⁻¹. (6.1)We ran simulations with five different values of weight factor 𝑤. The corresponding control functions are presented in Figure 2. As we increase 𝑤₃ and 𝑤₄, thereby increasing the cost of therapy, the control functions are decreasing. In Figure 8, we plot only the response of the HIV model to the obtained control function corresponding to weights 𝑤₁=1, 𝑤₂=0.1, 𝑤₃=21000, and 𝑤₄=10²¹ for brevity. Since the magnitude of CD4+ T cell, and virus population is much larger than the magnitude of the cost of the drugs treatment in the objective functional 𝐼(𝑢), these differences in magnitudes are balanced by this choice of 𝑤. Moreover, the performance of the gradient method with this choice of 𝑤 is shown in Figure 3. The number of required iterations to achieve the convergence is 307. Note that in Figure 3, ‖⋅‖ denotes the infinity norm, and 𝑥^(𝑖) denotes the state 𝑥(𝑡) in the 𝑖th iteration.

(a)

(b)

Figure 2 

Continuous solutions with 𝑤1=1 and 𝑤2=0.1 and different weights 𝑤3 and 𝑤4, on control cost: 𝑤3=25000 and 𝑤4=2×1021 (the x shape), 𝑤3=24000 and 𝑤4=1.5×1021 (the dashed dotted line), 𝑤3=19000 and 𝑤4=1021 (the solid line), 𝑤3=21000 and 𝑤4=1021 (the dotted line), 𝑤3=18000 and 𝑤4=8×1020 (the dashed line).

Figure 3 

Convergence of the gradient method.

6.2. MPC-Based STI Solutions

We implement therapy protocols based on the MPC algorithm (4.29) described in Section 4. All computations are carried out by MATLAB. We choose a sampling time of five days, Δ=5. Consequently, we do not worry about creating an explicit discretization of our differential equation; we simply use a numerical simulator to approximate our discretization. Moreover, we chose horizon length 𝑁=8. In particular, a finite horizon and a finite control space mean that we have, for each horizon, a finite number of possible control sequences. We implement MPC by using the Simulated Annealing (SA) metaheuristic combined with the Steepest Descent Explorer [47] for searching this space.

We ran simulations with three different values of weight factor 𝑤. The corresponding control functions are plotted in Figures 4, 5, and 6. In Figure 8, we present only the response of the HIV model to the STI solution corresponding to weights 𝑤₁=1, 𝑤₂=0.1, 𝑤₃=1800, and 𝑤₄=5×10⁵ for brevity.

Figure 4 

STI solution with 𝑤1=1, 𝑤2=0.1, 𝑤3=2000, and 𝑤4=105.

Figure 5 

STI solution with 𝑤1=1, 𝑤2=0.1, 𝑤3=1800, and 𝑤4=5×105.

Figure 6 

STI solution with 𝑤1=1, 𝑤2=0.1, 𝑤3=1500 and 𝑤4=5×105.

6.3. Piecewise Constant Solutions

In our implementation, we choose 𝐾₁=12 number of functions from the set 𝐶′(𝐵) as 𝑡^𝑖𝜉_𝑗,𝑖∈{0,1,},𝑗∈{1,…,6}. (6.2)

Furthermore, we set 𝐾₂=2 and 𝑆=6. By using controllability on the dynamical control system, the considered ranges for states and controls and the corresponding number of partitions in the construction of the set Ω_𝑁 are summarized in Table 2. Note that the selected values from the set 𝑈₁ in the construction of 𝑦_𝑗s are 0, 0.4, and 0.7. These values indicate off, moderate, and strong PI-therapy, respectively. Similarly, the corresponding values for RTI control are 0, 5×10⁻¹⁰, and 9×10⁻¹⁰ [7]. The desired values for 𝑚_𝑖s and 𝑀_𝑖s are summarized in Table 3. Implementing the corresponding LP model (5.37), we achieve 𝜆^∗=0.3902, which is close to the obtained number from Table 3, that is, 𝜆=min_𝑖=1,…,4{𝜇_𝑖(𝐼_𝑖(𝑢))}=0.3867. Figure 7 shows the resulting piecewise constant control pair. We solved the HIV model (2.1)–(2.6) with no treatment (𝑢_𝑃=𝑢_𝑅≡0) and with fully efficacious treatment (𝑢_𝑃≡0.7,𝑢_𝑅≡9×10⁻¹⁰). The numerical results for these two cases as well as the response of the systems to piecewise constant control, STI control, and continuous control are shown in Figure 8 for comparison. Except for STI solution, which is derived by considering linear cost functionals, the values of the objective functionals as well as corresponding values of the membership functions for other cases are given in Table 3.

Table 2 

Considered ranges for states and controls and corresponding number of partitions.

Table 3 

The values of the objective functionals and corresponding values of the membership functions.

Figure 7 

The piecewise constant control pair 𝑢∗𝑃 and 𝑢∗𝑅.

(a)

(b)

(c)

(d)

Figure 8 

Dynamic behavior of the state variables versus time in the presence of continues controls (the dashed line), piecewise constant controls (the solid line), STI controls (the dotted line), with no treatment (𝑢𝑃=𝑢𝑅≡0) (the × shape), and with fully efficacious treatment (𝑢𝑃≡0.7,𝑢𝑅≡9×10−10) (the dashed dotted line).

From Figures 8(a) and 8(b), we see drop in the number of CD4+ T cells, and a rise in viral load following the initial infection until about the third month. After this time, CD4+ T cells start recovering and virus starts decreasing due to the immune response, but can never eradicate virus completely. Then CD4+ T cells level decreases and viral load increases due to destruction of immune system in absence of treatment. Figures 8(b) and 8(c) show a clear correlation between the CTLe and virus population. As the virus increases upon initial infection, CTLe increases in order to decrease the virus. Once this is accomplished, virus decreases. Then virus grows back slowly, and this triggers an increase in the CTLe population. CTLe further increases in an attempt to keep the virus at constant levels but lose the battle because of virus-induced impairment of CD4+ T cell function, in absence of treatment. Memory CTL responses depend on the presence of CD4+ T cell help. Figures 8(a) and 8(b) show that, in presence of treatment, the virus is controlled to very low levels and CD4+ T cells are maintained above the critical levels. Therefore, immune response expands for relatively long time successfully. Furthermore, these Figures indicate inverse coloration between viral load and CD4+ T cells level. From Figures 8(c) and 8(d) interestingly, a decrease in CTLs occur in response to therapy can be observed. The extent of the decrease is directly correlated with the increase in treatment effectiveness which is consistent with experimental findings [48]. From Figure 8(a), we conclude that in the presence of treatment, the level of CD4+ T cells decreases very slowly as compared with untreated patient. Moreover, this figure shows that the average CD4+ T cells concentration with piecewise constant control is higher than continuous control. Interestingly, from Figure 8(d) we see that the intensity of immune response in the presence of continuous control is more than full treated patient until about the 42th month, and the intensity of immune response with piecewise constant control is higher than other controls after about the 31th month.

7. Conclusion

In this paper, a six-dimensional HIV model is considered, and a multi-objective optimal control problem is proposed to design therapy protocols to treat the HIV infection. The model permits drug therapy (RTIs and PIs) as controllers. We derived continuous treatment strategies by solving the corresponding optimality systems with the gradient method. In addition, structured treatment interruption (STI) protocols are found by using ideas from model predictive control. The results of simulations show that the proposed MPC-based method can effectively be applied to find STI-type solutions. In particular, a novel idea based on fuzzy aggregation and measure theory is proposed to find a compromise Pareto optimal solution, and numerical results confirmed the effectiveness of this approach. The method is not iterative and it does not need any initial guess of the solution and it can be utilized to solve multi-objective optimal control problems.