An Optimal Control Approach to HIV Immunology

This paper introduces a mathematical model which describes the dynamics of the spread of HIV in the human body. This model is comprised of a system of ordinary differential equations that involve susceptible cells, infected cells, HIV, immune cells and immune active cells. The distinguishing feature in the proposed model with respect to other models in the literature is that it takes into account cells that represent two distinct mechanisms of the immune system in the defense against HIV: the non-HIV-activated cells and the HIV-activated cells. With a view at minimizing the side effects of a treatment that employs a drug combination designed to attack the HIV at various stages of its life cycle, we introduce control variables that represent the infected patient’s medication. The optimal control rule that prescribes the medication for a given time period is obtained by means of Pontryagin’s Maximum Principle.

HIV infection results in a chronic, progressive disease that can lead to the destruction of the immune system. The disease is characterized by a high rate of viral replication, which results in the emergence of more virulent variants. HIV infection is currently characterized by the count of CD4+ T cells, by the amount of viral particles in the blood (viral load) and also by the clinical symptoms. Not all patients develop every stage of the disease, and the time elapsed between the infection and the manifestation of different clinical symptoms is highly variable, even though the causes of such a variation remain partly unknown.
To reproduce, HIV joins the membrane of the T4 cell, which is vital to the immune response. The virus releases its RNA and an enzyme, which produces the DNA of the virus. Then, the DNA of the virus enters the nucleus and joins to the DNA of the cell, taking full control. The result of this union is the pro-viral DNA, that produces the messenger RNA, which contains the genetic code of the virus. The messenger RNA then reaches the cytoplasm and produces virions, which leave the host cell as newly formed HIV's. Thus, when joined to a T4 cell, a single virus produces many potential threats to other cells.
By making quantification possible and unfolding non-trivial equilibrium points, the analysis of viral load in HIV infection has facilitated the management of the disease. It turns out that an exponential decrease in viral levels in plasma can be attained by the reverse transcriptase inhibitors and protease that are included in Antiretroviral Therapy (ART) [1] [2].
When HIV viruses invade the human body, they attack the CD4+ T cells in their way. When attacked, these auxiliary cells signal the presence of an invader to other immune cells (CD8+ T cells). The CD8+ T cells then respond to this signal and become Cytotoxic T Lymphocytes (CTL) by attempting to destroy the infected cells [3]- [6]. This process, which is not exploited in typical HIV models, plays an important role in the proposed approach. Indeed, a novel feature of the present work is the introduction of a variable to represent the CD8+ T cells. Such a variable is extremely important to the model, since it enables the decision maker to evaluate the interaction between the CD8+ T cells, CTL and the other variables in the model, such as the virus load.
This work proposes a simple mathematical model to describe the dynamics of the HIV taking into account the immune response. The proposed model introduces changes to several existing models in the literature [7]- [10]. As mentioned above, one such change is the introduction of a new variable to provide a more detailed description of the defense of the immune system. This variable represents the number of unactivated CD8+ T defense cells, which can become activated (i.e. HIV-specific, or CTL) after being warned by some CD4+ T cell. In contrast to the formulation proposed by [4], which accounts for the activated cells but does not model the activation mechanism, the proposed model keeps track of both the unactivated CD8+ T cells and the activated cells, thus taking into account the dynamics of the activation process.
Even though ART has produced undisputed advances in the treatment of HIV infection, it has been argued that the inhibitors that comprise ART may cause adverse effects; see for example [11]- [13]. Hence, one can argue that a compromise should be reached between the benefits obtained from ART and the adverse effects that it may cause. The ideal treatment should keep the benefits to a maximum degree while also minimizing the adverse effects. It is the development of such a treatment that we address in this paper, making use of the optimal control theory framework. We propose a dynamic model to represent the dynamics of the HIV infection, which takes into account the effects of the ART therapy and introduces control parameters that determine the intensity of the medication. An optimal control problem is then proposed to determine the optimal medication levels in such a way to maximize the benefits of the therapy, while keeping the medication to a minimum efficient level. The side effects assessment is another novelty of the paper, which also strives to provide insights into the effects in the optimal treatment by changing the parameters of the optimal control problem. In order to achieve our goal and minimize the side effects, optimal control theory is applied to the HIV infection model encompassing drug treatments. That, on the other hand, produces a problem whose solution is numerically obtained by standard algorithms, such as the gradient descent method; for more details on these algorithms we refer to [14]. The gradient algorithm makes use of both the solution to the system's dynamics under a given control sequence, i.e., the medication levels for the complete treatment, and the solution of the co-state equations under the same medication, to generate an improved control sequence, and so on, until convergence is attained. The solutions to both the system's dynamics and the co-state equations are obtained by means of finite difference algorithms.
This work is organized as follows. Section 2 introduces the model, while Section 3 derives the trivial equilibrium point and analyzes the stability of the model. Section 4 introduces the optimal control problem, which derives the optimal medication levels. In Section 5, numerical experiments are proposed to illustrate the proposed approach and shed light into the behavior of the controlled system. Finally, Section 6 concludes the paper.

Model Formulation
Let x represent the susceptible cells, i.e. the cells that can be infected with HIV, and define y as the cells that are already infected. In addition, assume that v represents the free viruses in the body, z denotes the defense cells (CD8+ T and B) and a z corresponds to the activated defense cells, i.e., cytotoxic T cells and plasma (activated CD8+ T and B cells). In the proposed model, the interaction of these variables is represented by the following set of differential equations: The parameters in the above equations are described in Table 3. Note that uninfected cells x are produced at a constant rate x λ , and decay at rate x µ . In addition, these cells are infected by the free viruses at rate v β . As for the infected cells y, they are produced from the interaction of their uninfected counterparts and the viruses, at rate v β , decay at rate y µ and are eliminated by the activated defense cells at rate v p . Free viruses v are generated from infected cells at rate v k , decay at rate v µ and are eliminated by means of the activated defense cells a z , at rate v p . The defense cells, in turn, are generated at a constant rate z λ , and decay at rate z µ .
Furthermore, they become activated by the viruses at rate z β . The activated defense cells are generated from the defense cells in the presence of the virus, at rate z β , and decay at rate z µ . It is worth stressing that these cells fight against the infected cells y and the free viruses v, as previously mentioned in the description of the model dynamics.
The virus replication mechanism is depicted in Figure 1. Observe that free viruses (v) and uninfected cells (x, CD4+ T) produce infected cells (y) that, in turn, produce new free viruses. Figure 2 illustrates the activation of the defense cells (z, CD8+ T). Note that they are activated in the presence of free viruses.
In the model, one can notice the introduction of two variables z and a z to represent the defense cells (CD8+ T) and the HIV activated defense cells (CTL), respectively. Figure 3 illustrates the action of the activated defense cells, which eliminate infected cells y and free viruses v.
It is worth point out that, in the proposed model, the immune CTL cells are activated directly by the free viruses. In contrast, some previous works in the literature assume they are activated by infected cells, healthy cells and also CTL cells [15]. Observe that this is a simplification of the proposed model, which takes advantage of the fact that all elements that act to activate the immune CTL cells are mediated directly or indirectly by the virus. This hypothesis in our model simplifies the complex phenomena of activation and inactivation of the immune system by cytokines and dendritic cells (antigen presenting cells) [15]. Such a simplification is very useful in the optimal control formulation, which results in a larger and more complex set of equations.     (1) and solving for the remaining unknowns, one finds the trivial equilibrium point of the dynamic system in Equation (1), described below:

Stability of Trivial Equilibrium Point
By Hartman-Grobman's theorem [16], one can say that an equilibrium point is stable if the sign of the real part of the eigenvalues of the Jacobian matrix, evaluated at this point, is negative. Making use of this theorem, it is possible to show that if the basic reproduction number of the virus 0 R is less than one unit, the trivial equilibrium point o P in (2) is stable. In such a scenario the infection does not propagate in the human body.
The Jacobian matrix of the system (1), ( ) J P , is given by: evaluated at values corresponding to the equilibrium points. The Jacobian matrix evaluated at the point of trivial equilibrium (2) is given by: The eigenvalues of the matrix ( ) o J P are the roots r of the characteristic polynomial ( ) p r given by: It can be noted that Applying Routh-Hurwitz criterion [17], one can see that the system is stable if the last term in Equation (6) is positive. Hence, the system is stable if Lets assume that one virus infects one person. During its average survival time ( ) 1 v µ , the possibility of encountering and infecting one target cell v β in a completely susceptible population of CD4+ T cells This infected cell, releases on average v k viruses. Hence, 0 R is the average number of secondary viruses produced by one virus infecting a cell in a completely susceptible population of CD4+ T cells.
So, if 0 1 R < (basic reproduction number of the virus), as it can be seen in (7), the trivial equilibrium point o P given in (2) is stable. That is, in this case, the infection does not propagate in the body.

Non-Trivial Equilibrium Points
Equating to zero all the differential equations given in (1), one comes up to the following general formulation for the equilibrium points of the dynamical system: where: It is possible to observe that when 0 v = in (9), and (10) or (11), we obtain the trivial equilibrium point given in (2). However, if 0 v ≠ , i.e, for an HIV infected individual, equaling Equation (10) and (11) it is possible to obtain, after some numerical manipulation, the following polynomial equation of degree 3: where the coefficients are: µ µ µ µ β µ µ µ β β µ µ β β λ µ µ β β λ µ µ µ β λ µ µ µ β λ µ β λ The coefficients of Equation (12) are positive when 0 1 R < , implying that there is no positive solution; hence 0 v = is a solution. When 0 1 R > , 0 0 a < , and as 0 R increases, the coefficient 1 a becomes negative, followed by 2 a . Hence, for 0 1 R > , a unique change of sign occurs between successive coefficients, and according to Descartes rule of signs, there is a unique positive solution.

Optimal Control Formulation
In order to study the effect of treating HIV infection with the use of a cocktail of drugs, we introduce two control variables ( 1 u and 2 u ) in the dynamic system described in Equation (1). The control variable 1 u represents the effect of inhibitors of reverse transcriptase, integrase and input. These drugs protect the uninfected cells x, preventing their infection, i.e., keeping them from turning into infected cells y. To account for this effect, we introduce in the basic model of Equation (1) a new state variable p x , which represents the cells that are protected by the action of the inhibitors. The control variable 2 u simulates the effect of the protease inhibitor, which blocks infected cells y, preventing them from releasing new viruses v in the body. To account for this action, we introduce a new variable b y in our model. This variable describes the infected cells that are blocked by the action of the protease inhibitor. With the modifications described above, the studied dynamic system (1) can be reformulated as: The initial conditions were obtained by running the uncontrolled model from Equation (1) Alternatively, doctors could prescribe average values of the sequences 1 u and 2 u over predetermined intervals, such as monthly intervals for example, re-evaluating the treatment between successive intervals.
Observe that, if the dosages 1 u and 2 u were to be kept constant, one could replace x µ for ( ) 1 x u µ + and y µ for ( ) 2 y u µ + in terms involving 1 u and 2 u , and all the results regarding the equilibrium points and their stability obtained in the foregoing section would be retrieved here. In that case, the number of primary virus replication after treatment (control), represented by c R , is given by: Note that the terms multiplying 0 R in the above equation arise due to the (now constant) control variables ( ) In the remainder of this paper, we analyze the system's dynamics with 1 : u t + →  and 2 : u t + →  varying over time, and strive to obtain an optimal treatment, prescribing the values of these variables at each time in such a way that an optimal compromise between the efficiency of the treatment and its side effects is obtained.
It is worth reinforcing that, in relation to typical HIV models that simulate the production rates and infectability of the virus, the proposed model includes two new variables in the dynamics: protected uninfected cells p x and blocked infected cells b y , which do not release new viruses. This modeling choices are explained partly by the fact that the antiretroviral treatment stops the cell transcription process. The cell becomes healthy and the virus is destroyed, but that happens only after the contact between the cell and the virus. Bearing that in mind, instead of having the control be set in the infection force by a mass action law, we assume that the cells became resistant to infection (or protected). We also assume that the effect of the protease inhibitor is to prevent the release of viable viruses. Consequently, only infected but untreated cells release viable viruses. Upon being treated, the infected cell becomes blocked and does not release viable viruses thereafter.
The introduction of the new variables, namely protected uninfected cells p x and blocked infected cells b y , allows a better understanding of the dynamics. Moreover, taking into account the number of such cells is important to assess the effects of the prescribed treatment. Another innovation is the introduction of the control variables, which enables us to model the compromise between medication and side-effects.
Observing the model we note that the defense cells CD8+ T (z) produce the activated cells CTL ( ) a z , while the latter contribute to the reduction of infected cells (y) and free viruses (v). Moreover, the control variable 1 u , corresponding to the reverse transcriptase inhibitor, protects the cells CD4+ T (x) from being infected, thus producing protected cells p x . In addition, the control variable 2 u corresponds to the integrase inhibitor and keeps the infected cells y from producing viruses v. The infected cells under the effect of this inhibitor are called blocked cells, and are denoted by b y . Table 1 provides a description of the variables that were added to the model. Figure 4 illustrates the dynamics of the system described by Equation (17).

Objective Function
System (17) is now optimized with respect to control parameters 1 u and 2 u . For this reason, both are allowed to vary with time. To reach a compromise between medication and side effects, we introduce an optimal control problem aimed at maximizing the number of protected cells while also mitigating the side effects. The optimal control problem is defined as follows:  while also trying to minimize the drug administrations ( 1 u and 2 u ). We are assuming that the higher 1 u and 2 u , the higher the side effects.

The Hamiltonian
To find the optimal control variables * 1 u and * 2 u that solve Problem (20), we make use of Pontryagin's maximum Table 1. New variables after control.

Control variables Label
Reverse transcriptase reverse, integrase and entrance inhibitors where , 1, , 7 j w j =  , are the co-state variables, that determine the adjoint systems. It is well known that the optimal solution must satisfy both the original and the adjoint system of equations.The variables 1 v and 2 v are penalty multipliers (slack variables), added to the model to ensure that the constraints 1 0 u ≥ and 2 0 u ≥ are satisfied. For the optimal controls * 1 u and * 2 u , it holds that Hence, the optimal control for Problem (20) is characterized by (25). 1 u and * 2 u that simultaneously solve the initial value problem in Equation (17) and the final value problem in Equation (26). In this paper, we find the optimal control trajectories iteratively, employing a classical specialized gradient algorithm [14] and solve the differential equation systems by means of finite difference methods.

Numerical Results
For the numerical simulations we have used the dataset described in Table 2 and Table 3, which was also studied in [8] [9] [20]. For this data set, the reproduction number 0 3.6 R = indicates an infection. - [23]. The system rapidly approaches the non-trivial equilibrium point.

The Uncontrolled System
In Figure 5, one can see a steep increase in the number of viruses and infected cells in the first month. In contrast, the number of susceptible cells presents a substantial decrease over the same period. In the subsequent  periods, these variables stabilize and gradually approach a non-trivial equilibrium point. As previously mentioned, the proposed model describes the typical behavior of the immune system in the presence of HIV. The difference is that it also allows a separate analysis of the activated defense cells. In Figure 5, one can notice a migration of immune cells into the compartment of activated defense cells. We believe that this could lead, over the years, to the saturation of the immune system. Furthermore, the results seem to suggest a trend for a regular non-trivial dynamical system balance.

Adding Medication: The Optimal Control Problem
The optimal control problem in Equation (20) was solved for a one year treatment period. The initial condition for the state variables was obtained by running the uncontrolled system in Equation (1) for a one year period. The last value of each variable (after 365 days) was taken as the initial value for the same variable in the controlled problem. In this work, we first consider Case 1, with 3 1 10 c − = , 2 1 c = and 3 1 c = , which will be used as our benchmark for model evaluation. Note that, under such parameters, the objective function in (20) evaluates the number of protected cells versus the side effects attained from the medication to be taken to accomplish such a number. The small value of 1 c is due to the large number of immune cells, which render the first term in Equation (20) quite significant. The optimal trajectories of the system for the selected parameters are depicted in Figure 6.
Note in Figure 6 that the first week of treatment causes a substantial decrease in the unprotected (susceptible) CD4+ T cells, with an increase in the number of protected cells. That is because the treatment makes most of these cells become protected cells. Note that, since the medication efficiently fights the infection, the number of HIV-specific CD8+ T cells can be reduced with respect to the uncontrolled dynamics in Figure 5. That leaves the imune system free to combat other infections. With regards to the infected cells, blocked and unblocked, they rapidly vanish. Note that the unblocked cells vanish more rapidly than their blocked counterparts. Note also that the number of remaining viruses rapidly decays to zero, its trivial equilibrium point. It is also noteworthy that the optimal dosage of the transcriptase inhibitor ( ) 1 u is significantly higher than that of the integrase inhibitor ( ) 2 u . In the next section we verify how the variation in the parameters of the optimal control problem affects the dynamics. It is also noteworthy that the medication levels are much higher in the early treatment period, and are stabilized to lower levels in a very short time span.

Sensitivity Analysis
In this section we examine the effects of the parameters of the optimal control problem in Equation (20). We vary the parameters 2 c and 3 c and verify the effects of the variations in the controlled system dynamics. The E. F. Arruda et al. 1127 objective is to verify the implications of different compromises between immune response and side effects into the medication levels as time elapses.
In Case 2, we examine the influence of parameter 2 c , making 1 0.001 c = , 2 100 c = and 3 1 c = . In this case, we assume that the reverse transcriptase inhibitor presents more side effects than the protease inhibitor, thus causing an augmented value of 2 c with respect to Case 1. The results are conveyed in Figure 7. It can be noted that, with the increase of 2 c there is a decrease in both the number of protected cells ( ) p x and unprotected cells (x). The evolution of the virus shows that, in this case, the viruses get extinct more slowly than in Case 1. The same applies to the infected cells y and b y . Apparently, the changes do not cause significant changes in the values of z and a z . With regards to the control variables, i.e. the medication levels, the levels of 1 u are reduced in this example.
In Case 3 we make 1 0.001 c = , 2 1 c = and 3 100 c = , assuming that 2 u produces more side effects than 1 u . The results are depicted in Figure 8. One notes that the viruses get extinct slightly slower than in Case 1, with the same applying to infected cells y. Once again, the changes do not cause significant changes in the values of z and a z . The medication levels 1 u and 2 u are similar to those in Case 1. The increase in 3 c does not seem to affect the levels of 2 u significantly, for these levels were already very reduced in Case 1. In Case 4 we make 1 2 0.001, 100 c c = = and 3 0.01 c = . In this case, in addition to 1 u being very harmful, we assume that 2 u has reduced side effects. The results are presented in Figure 9. One notices that the dynamics of the system remains similar to that in the other cases. One difference is that the viruses and the unblocked infected cells y get extinct much faster. Another significant difference is that, in this case, the levels of 2 u exceed those of 1 u .

Concluding Remarks
This paper introduces a novel model of HIV dynamics. In contrast to typical HIV dynamics models, the proposed model explicitly describes the protected CD4+ T cells and the HIV-specific CD8+ T cells. That allows us to explicitly understand and quantify these effects, thus providing a better understanding of the system's dynamics.
The dynamics of the proposed model can be influenced by the use of Anti-retroviral Therapy (ART), which has produced significant advances in the treatment of HIV infection. Such a therapy, however, is associated to side effects, which should be avoided whenever possible. To take account of the side effects and produce a desirable compromise between treatment effectiveness and side effects, we propose an optimal control approach  which prescribes an optimal treatment aimed at maximizing the benefits of the treatment, while also minimizing the side effects. The optimal control problem is solved by means of Pontryagin's maximum principle and the optimal medication levels are numerically found by a standard gradient algorithm. The proposed optimal control formulation is analyzed in the light of selected numerical examples, which provide insight into the behavior of the system under different compromises between medication effectiveness and side effects. We evaluate the sensitivity of the method with respect to the parameters in the optimal control functional by means of a series of examples, which shed light on the optimal trajectories of the variables with respect to different compromises between efficiency and side-effects. We analyze four different cases with distinct levels of side effect introduced by drugs 1 u and 2 u , which give rise to different treatment policies that take into account the prescribed compromises between these side effects, with a view at keeping the drug benefits while keeping the side effects to the minimum.