OPTIMAL CONTROL OF AN HIV MODEL WITH CTL CELLS AND LATENTLY INFECTED CELLS

. This paper deals with an optimal control problem for an human immunodeﬁciency virus (HIV) infection model with cytotoxic T-lymphocytes (CTL) immune response and latently infected cells. The model under con-sideration describes the interaction between the uninfected cells, the latently infected cells, the productively infected cells, the free viruses and the CTL cells. The two treatments represent the eﬃciency of drug treatment in inhibiting vi- ral production and preventing new infections. Existence of the optimal control pair is established and the Pontryagin’s minimum principle is used to charac- terize these two optimal controls. The optimality system is derived and solved numerically using the forward and backward diﬀerence approximation. Final- ly, numerical simulations are performed in order to show the role of optimal therapy in controlling the infection severity.


1.
Introduction. The human immunodeficiency virus (HIV) is a virus that gradually weakens the immune system. It is considered as the main cause for several deadly diseases after the resulting acquired immunodeficiency syndrome (AIDS) is reached. After this stage, the immune system fails to play a crucial role, which is to protect the whole body against harmful intruders. This failure is due to destruction of the vast majority of CD4 + T cells by the HIV virus, reducing them to an account below 200 cells per µl [3]. On the other hand, according to the world health organization [22], 36.7 million people living with HIV, 1.8 million people becoming newly infected with HIV and more than 1 million deaths annually. Based on these alarming statistics, HIV becomes a major global public health issue. Hence, mathematical modeling of human immunodeficiency virus (HIV) infection is one of the great current interest [2,4,7,8,[14][15][16][17]20,21,23]. The basic viral infection model with four dynamics components including the uninfected cells, the infected cells, the free viruses and the cytotoxic T-lymphocytes (CTL) immune response was first studied in [14]. However, upon the virus entry to the uninfected cell, there exists a latent period after that the cell becomes infected. In a recent work, the model describing HIV viral dynamics taking into account this fifth compartment representing the infected cells in latent stage is formulated and studied in [2]. The authors study the global stability of the endemic states and illustrate the numerical simulations in order to show the numerical stability for each problem steady state. This paper will be focused on studying an optimal control for the HIV infection model given in [2]. For this purpose, we will consider the following nonlinear differential equations: With the initial conditions x(0) = x 0 , s(0) = s 0 , y(0) = y 0 , v(0) = v 0 and z(0) = z 0 . In this model, x, s, y, v and z denote the concentration of uninfected cells, latently infected cells, infected cells, free virus, and CTL cells, respectively. Susceptible host cells CD4 + T cells are produced at a rate λ, die at a rate d 1 x and become infected by virus at a rate k 1 xv x + v . The exposed cells die at a rate (d 2 + k 2 )s. The infected cells increase at rate k 2 s, decay at rate d 3 y and are killed by the CTL response at a rate pyz. Free virus is produced by infected cells at a rate ay and decay at a rate d 4 v. Finally, CTLs expand in response to viral antigen derived from infected cells at a rate cyz and decay in the absence of antigenic stimulation at a rate bz. Note that this model (1), includes a saturated rate, called the saturated mass action [2,17] which describes better the rate of viral infection. The two constants u 1 and u 2 stand for the efficiency of treatment in blocking new infection and in inhibiting viral production, respectively. So far, there is no effective treatment that eradicates HIV virus; However there are some therapies that reduce HIV viral replication including reverse transcriptase inhibitors (RTIs) and protease inhibitors (PIs) [11]. We note that 54% of infected adults and 43% of infected children are receiving lifelong antiretroviral therapy (ART) [22]. Hence, the need of studying the effect of optimal therapy control on HIV viral dynamics. The paper is organized as follows. The next section is dedicated to the qualitative analysis of the model, followed in Section 3 by the optimization analysis of the model. In Section 4, we construct an appropriate numerical algorithm and give some numerical simulations. We conclude in the last section.
2. Qualitative analysis of the model.

Positivity and boundedness.
It is straightforward to show that the solutions of the system (1) are positive for any positive initial conditions. We examine the boundedness of the solutions to the system (1) with non-negative initial conditions. For this purpose, we introduce the variable X(t) = x(t)+s(t)+y(t)+ d 3 2a v(t)+ p c z(t). Using the system (1), we havė this proves, by virtue of Gronwalls Lemma, that X(t) is bounded, and so are the functions x(t), s(t), y(t), v(t) and z(t). Therefore, every local solution can be prolonged up to any time t m > 0, which means that the solution exists globally.
2.1. Steady states. First, the basic reproduction number of our problem is given as follows: The model (1) admits three steady states that are typical in viral dynamics models i) The disease-free steady state E f = ( λ d1 , 0, 0, 0, 0). ii) The immune-free steady state E 1 = (x 1 , s 1 , y 1 , v 1 , z 1 ), where iii) The infected steady state E 2 = (x 2 , s 2 , y 2 , v 2 , z 2 ), where From the components of E 1 , it is clear that when R 0 > 1 this endemic point exists. In what follows, we will use the following CTL immune response reproduction number R CT L to classify the model system dynamics .
Depending on the value of this CTL immune reproduction number R CT L ; in other words depending on R 0 and the parameters describing the effectiveness of CTL like c and b, we will study the stability of the endemic equilibria. First, observe that the second endemic state E 2 = (x 2 , s 2 , y 2 , v 2 , z 2 ) exists when R CT L > 1. This endemic state is also called an interior equilibrium.
We explain the existence of this endemic equilibrium E 2 as follows. We recall first, that in this state both the free viruses and CTL cells are present. Assume that R 0 > 1, in the total absence of CTL immune response, the infected cell load per unit time is ) . Via the fifth equation of model (1), CTL cells are reproduced, due to infected cells stimulating, with the rate we will have the existence of the endemic equilibrium E 2 .
2.2. Critical drug efficacy. In the system (1), the efficiency of drugs RTIs and PIs are denoted by the two terms 1 − u 1 and 1 − u 2 , respectively. Therefore the values, u i = 0 and u i = 1 (i = 1, 2), reflect the completely ineffectiveness and the perfectly effectiveness of the chosen drug. For simplicity, the efficiency of the two drugs RTIs and PIs will be combined to a new term that reflects the overall efficacy 1 − θ = (1 − u 1 )(1 − u 2 ). The basic reproduction number with this new term is . In order to find the critical drug θ crit , we put R 0 = 1 On the other hand, if θ ≤ θ crit , then implies that the disease persists.
3. The optimal control problem. Now, to state the optimization problem, we first consider that u 1 and u 2 vary with time. The problem (1) becomes The objective functional that we seek to minimize is defined by the following quadratic costs: where t f is the period of treatment, Q 1 , Q 2 , Q 3 are the cost coefficients [12] for s, y and v respectively. Also, A 1 and A 2 represent the financial and physiological costs associated with the administration of RTIs and PIs respectively. A 1 > 0, A 2 > 0 reflect the extent of toxicity of the administered drugs. The two control functions, u 1 (t) and u 2 (t) are assumed to be bounded and Lebesgue integrable. The quadratic terms u 2 1 (t) and u 2 2 (t) reflect the expectation that the side-effects of the drugs are nonlinear [1,18]. Our target is to minimize the objective functional defined in equation (4) by decreasing the number of the exposed cells, infected cells and the viral load and minimizing the cost of treatment. In other words, we are seeking optimal control pair (u * 1 , u * 2 ) such that: Here a i and b i act as the lower and upper bound on the therapeutic efficacy u i (i = 1, 2).

3.1.
Existence of an optimal control pair. In order to prove the existence of an optimal control pair and the uniqueness of the optimal system, we need the upper bounds of the state variables [5]. Note that the super-solutions x(t), s(t), y(t), v(t) and z(t) are uniformly bounded. Now, we prove the existence of an optimal control pair. We illustrate the sufficient condition using the result of Fleming and Rishel (see Theorem 4.1, pp 68−69 in [6]).
Theorem 3.1. There exists an optimal control pair (u * 1 , u * 2 ) ∈ U such that Proof. To use the existence result in [6], we must check the following properties: (P 1 ) The set of controls and corresponding state variables is nonempty.
The control U set is convex and closed. (P 3 ) The right-hand side of each equation of (3) is continuous, is bounded above by a sum of the bounded control and the state, and can be written as a linear function of a control pair (u 1 , u 2 ) with coefficients depending on time and state. (P 4 ) The integrand of the objective functional L( X, u, t) is convex on U . (P 5 ) There exists constants c 1 and c 2 > 0, such that the integrand L( X, u, t) of the objective functional satisfies Where with X = (x, s, y, v, z) and u = (u 1 , u 2 ).
(P 1 ) The boundedness of the system of state equations under the two controls (3) ensures the existence of a solution. We can therefore conclude that the set of controls and the corresponding state variables is non-empty, which gives condition (P 1 ). (P 2 ) The control set is convex and closed by definition, which gives condition (P 2 ). (P 3 ) It can be easily seen that the right-hand side of the each equation of (3) is continuous and can be written as: Since the system is bilinear in u 1 , u 2 and the solutions of the system (3) are bounded, we have where X = (x, s, y, v, z) and u = (u 1 , u 2 ) and c 0 depends on the coefficients of the state system (3). (P 4 ) Let X = (x, s, y, v, z) be the state vector and u = (u 1 , u 2 ) and u = (u 1 , u 2 ) be two control vectors. Then, for 0 < < 1, we have So, L is convex on U . (P 5 ) For the last part, we have where c 1 is any non-negative real number and c 2 = min(A 1 , A 2 ) . We conclude that there exists an optimal control pair u 3.2. Optimality system. We have already proved the existence of an optimal control pair which minimizes the functional (4) subject to (3). We now use Pontryagin's minimum principle to derive the necessary conditions for the optimal control [13].
However, it will be worthy to notice that there exist other effecient optimal control methods for some problems dealing with switching or time-delay optimal control [9,10,19]. Accordingly the Lagrangian for our control problem will be defined as follows: with w 11 (t), w 12 (t), w 21 (t), w 22 (t) ≥ 0 are penalty multipliers satisfying and By applying Pontryagin's minimum principle in state, we obtain the following theorem Theorem 3.2. For any optimal controls u * 1 , u * 2 , and solutions x * , y * , s * , v * and z * of the corresponding state system (3), there exist adjoint variables, λ 1 , λ 2 , λ 3 , λ 4 and λ 5 satisfying the equations with the transversality conditions λ i (t f ) = 0, i = 1, ..., 5.
Moreover, the optimal control is given by Proof. The adjoint equations and transversality conditions can be obtained by using Pontryagin's minimum principle [13], such that The optimal control u * 1 and u * 2 can be solved from the optimality conditions, ∂L ∂u 1 (t) = 0, ∂L ∂u 2 (t) = 0.
Then, by calculating the derivative of L, we obtain using (7), (8) and the boundedness of the two controls in U , we obtain the following optimal control: .

this fact, implies that
We can derive the integral equations and their estimates for the remaining seven states and adjoint variables. Combining all the ten of these estimates, we get whereC 1 andC 2 are dependent on the system coefficients as well as the bounds on state and adjoint variables. Then, If we choose such that >C 1 +C 2 and t f < 1 2 ln( −C 1 C 2 ), then p i =p i and q i =q i for i = 1, ..., 5. Hence, the solution of the optimality system is unique for t f sufficiently small. Then, the optimal controls, u * 1 and u * 2 are characterized in terms of the unique solution of the optimality system. The above optimal controls give an optimal treatment strategy for the HIV infected patient under the scenario of these two types of drug treatments.
4. Numerical Simulations. In this section, we illustrate the numerical simulation of the various results for the pharmacokinetic model, as well as the consequence of application of optimal combination therapy of RTIs and PIs. Now, we study the influence of efficiency u 1 and u 2 (RTIs and PIs respectively) on the basic reproduction number R 0 . Recall that the infection die out (persists) whenever R 0 < 1(> 1), which is equivalent to θ > θ crit (< θ crit ). The parameters of our numerical simulations are inspired from [2]; i.e. λ = 10, d 1 = 0.0139, k 1 = 0.04, d 2 = 0.0495, k 2 = 1.1, d 3 = 0.5776 , a = 40, d 4 = 0.6, p = 0.0024, c = 0.15 and b = 0.5. Prior to the initiation of treatment R 0 = 4.41 > 1, and u 1 = 0, u 2 = 0. The goal is to keep increasing u 1 and u 2 such that R 0 is driven to a value less than unity. In this case, u 1 and u 2 must be chosen such that θ > θ crit = 0.77. We illustrate this by a surface plot and a contour plot in Fig. 1. One can easily observe that for u 1 = 0 and u 2 = 0 the value of R 0 reaches its maximum value of 4.41. By increasing u 1 and u 2 from 0 to 1, we observe that the value of R 0 gradually decreases and tends towards 0 (corresponding to u 1 = 1, u 2 = 1). This clearly reflects the impact of the efficiency in terms of clearance of the infection.
For the purpose of this work, we will assume that the minimum and maximum of the controls to be a i = 0.01 and b i = 0.99, i = 1, 2. The optimal control pair u * = (u * 1 , u * 2 ) is obtained for a period of 120 days of observation. The cost coefficients that were introduced in the definition of the objective functional (4) are Q 1 = 10 −5 , Q 2 = 10 −5 , Q 3 = 10 −5 , A 1 = 5000 and A 2 = 5000 [12].
In order to solve our optimization system, we will use a numerical scheme based on forward and backward finite difference approximation. Hence, we will have the following numerical algorithm Step 1: Step 2: for i = 0, ... , n-1, do: 2 , a 2 )), end for Step 3: for i = 1, ... , n, write The numerical algorithm. Figure 2 shows that with control the amount of the uninfected cells is higher than those observed for without control case.     The role of therapy control is also observed in Fig. 5. It was shown that with control, the number of HIV virus dies out after the first days of therapy, while without control it stays equal to 138.88. This indicates the impact of the administrated therapy in controlling viral replication.
The CTL cells are clearly affected by the control. This is shown in Fig. 6; indeed the curve of CTL cells converges towards zero with control, while without any   The two optimal controls u 1 and u 2 corresponding to blocking new infections and inhibiting viral production are represented in Fig. 7. The two curves represent the drug administration schedule during the period of treatment. Both controls start from zero and oscillate between zero and one. When the first immune boosting drug is administered at full scale, the second drug is at its lowest and vice versa. In this case, the new infection is totally blocked.
In Fig. 8, we test the hypothesis that a patient stops treatment at 120 th day of the infection, we notice that the curves of u 1 and u 2 drop toward zero. This sudden discontinuation of treatment after 120 th day results in the rebound of the levels of HIV virions (see Fig. 9). From this figure, it can be seen that the uninfected cells fall dramatically to a very low level; while the infected cells suddenly jump to high  Figure 7. The optimal control u 1 (left) and the optimal control u 2 (right) versus time. amount after stopping the therapy. One can conclude that therapeutic protocol should be administered over the entire period of infection.  Conclusion and discussion. In this paper, we have presented a model of the human immunodeficiency virus (HIV) infection where the patient is subjected to the combination therapy of reverse transcriptase inhibitors (RTIs) and protease inhibitors (PIs). We have studied the global stability for the disease-free steady state, the immune-free steady state and the infected steady states and obtained a critical drug efficacy in terms of model parameters. A combination efficacy of both RTIs and PIs was defined in such a way that the persistence or clearance of infection depend on whether this combination efficacy is less or greater than the critical drug efficacy.
We have formulated an optimal control problem with the objective is minimizing the levels of latently infected cells, infected cells and HIV virus as well as the therapeutic cost of the combination treatment. The existence and uniqueness of the optimal controls using Pontryagin's minimum principle is established. The problem was solved numerically using backward and forward finite difference scheme. It was shown that with the two optimal treatments, the number of the healthy CD4 + T cells increases remarkably while the number of latently infected CD4 + T cells and infected CD4 + T cells decreases significantly. In addition, it was also observed that, with the control strategy, the viral load decreases considerably compared with the model without control which can improve the patient's life quality. Also, it was observed, that the drug therapy should be administrated without any stops during all the period of the infection. Moreover, the determination of the critical drug efficacy and the optimal efficacy would be of great use for the biomedical community in designing individual treatment protocols.