Optimal control of a delayed HIV model

We propose a model for the human immunodeficiency virus type 1 (HIV-1) infection with intracellular delay and prove the local asymptotical stability of the equilibrium points. Then we introduce a control function representing the efficiency of reverse transcriptase inhibitors and consider the pharmacological delay associated to the control. Finally, we propose and analyze an optimal control problem with state and control delays. Through numerical simulations, extremal solutions are proposed for minimization of the virus concentration and treatment costs.

1. Introduction. Infection by human immunodeficiency virus type 1 (HIV-1) has many quantitative features [26]. Mathematical models for HIV infection can provide insights into the dynamics of viral load in vivo and may play a significant role in the development of a better understanding of HIV/AIDS and drug therapies [37]. Cytotoxic T lymphocytes (CTLs) play a critical role in antiviral defense by attacking virus-infected cells. It is believed that CTLs are the main host immune factor that determine virus load [22]. When HIV invades the body, it targets the CD4+ T cells. These cells can be considered the command centers of the immune system. The CTLs are cells that set out to eliminate infection by killing infected cells [6]. Several mathematical models have been proposed for HIV-1 infection with CTLs response: see, e.g., [1,6,22,36] and references cited therein.
Time delay plays an important role in the dynamics of HIV infection. Intracellular delay, that is, the delay between initial infection of a cell by HIV and the release of new virions, was considered in the models proposed by [5,13,18,19,20,21,34,37]. Here, we enrich the undelayed mathematical model proposed by [22], which considers the action of CTLs in the immune system, by introducing a discrete time delay that represents an intracellular delay. State delays for such type of models have been already introduced, e.g., in [12]. However, in our case we also model the important pharmacological delay that occurs between the administration of drug and its appearance within cells, due to the time required for drug absorption, distribution, and penetration into the target cells [27]. In the context of anticancer therapy, the idea to represent delay effects in drug kinetics and dynamics was presented in [33] and developed in [15].
Optimal control is a branch of mathematics developed to find optimal ways to control a dynamic system [4,8,28]. Optimal control theory has been applied with success to HIV models: see, e.g., [6,12,14,30,32] and references cited therein. Here, we introduce a control function, which represents the efficiency of reverse transcriptase inhibitors, and consider a delay in the control function representing the pharmacological delay. Our aim is to determine the control function that minimizes the concentration of virus and the treatment costs. To the best of our knowledge, this is the first time an optimal control HIV problem with delay in state and control variables is investigated.
The paper is organized as follows. The model with intracellular delay is formulated in Section 2 and local stability is proved for any time delay. In Section 3, we introduce a control function in the delayed model of Section 2 and analyze an optimal control problem with intracellular and pharmacological delays. Section 4 is devoted to numerical simulations for the stability of the equilibrium points and the computation of extremals for the optimal control problem with state and control delays. We compare the extremal of our optimal control problem with state and control delays with the solutions of the uncontrolled problem and the control problem with delay in the state variable only. We end with Section 5, where we discuss the established results.
2. Intracellular delayed mathematical model. In this section, we propose a delayed mathematical model for HIV-1 infection. We consider the undelayed model proposed by [22] and introduce a discrete intracellular time delay. The model considers four state variables: Z(t) represents the concentration of uninfected cells, I(t) represents the concentration of infected cells, V (t) represents the concentration of free virus particles, and T (t) represents the concentration of CTLs at time t. The following assumptions are made to describe the cell dynamics [22]: uninfected cells are produced at a constant rate λ, and die at a rate mZ. Infected cells are produced from uninfected cells and free viruses at a rate rV Z and die at rate uI (the average lifetime of an infected cell is 1/u). Free viruses are produced from infected cells at rate kI and declines at rate vV (the average lifetime of a free virus particle is 1/v). The rate of CTLs proliferation in response to antigen is given by aIT . In the absence of stimulation, CTLs decay at rate nT . Infected cells are killed by CTLs at rate sIT . The intracellular delay, τ , represents the time needed for infected cells to produce virions after viral entry [12,37], called the eclipse phase [25]. The model we propose is given by the following system of ordinary differential equations: The initial conditions for system (1) are of continuous functions mapping the interval [−τ, 0] into R 4 . The usual local existence, uniqueness and continuation results apply [11,16]. Moreover, from biological OPTIMAL CONTROL OF A DELAYED HIV MODEL   3 meaning, we further assume that the initial functions are nonnegative: From [37, Theorem 2.1], it follows that all solutions of (1) satisfying (2) and (3) are bounded for all time t ≥ 0, which ensures not only local existence but the existence of a unique solution (Z(t), I(t), V (t), T (t)) of (1) with initial conditions (2)-(3) for all time t ≥ 0. The equilibrium points are independent of the delays. Their stability depends, however, on the delays. The equilibrium points of (1) are studied in [29,37]. System (1) has an infection-free equilibrium E 0 = λ m , 0, 0, 0 , which is the only biologically meaningful equilibrium, if R 0 = kλr muv < 1. Let R 1 = knr mav . If 1 < R 0 < 1 + R 1 , then system (1) has a unique CTL-inactivated infection equilibrium E 1 given by Whenever R 0 > 1 + R 1 , system (1) has also a CTL-activated infection equilibrium E 2 given by The proofs of these facts are found in [29,37]. Here we prove the local asymptotic stability of the equilibrium points E 0 , E 1 and E 2 for any time delay τ .
Theorem 2.1 (Local stability of the equilibrium points of (1)). If R 0 > 1, then the infection-free equilibrium E 0 is unstable for any time-delay τ ≥ 0. If R 0 < 1, then E 0 is locally asymptotically stable for any time-delay τ ≥ 0. If R 0 = 1, then we have a critical case. If R 0 > 1 + R 1 , then the CTL-inactivated infection equilibrium E 1 is unstable for any time-delay τ ≥ 0. If R 0 < 1 + R 1 , then E 1 is locally asymptotically stable for any time-delay τ ≥ 0. If R 0 > 1 + R 1 , then the CTL-activated infection equilibrium E 2 is locally asymptotically stable for any time-delay τ ≥ 0.
Proof. Consider the following coordinate transformation: where (Z,Ī,V ,T ) denotes any equilibrium point of system (1). The linearized system of (1) is of form We can express system (4) in matrix form as follows: where A 1 and A 2 are 4 × 4 matrices given by  The characteristic equation of system (4) for any equilibrium point is given by (see, e.g., [16]), where I d denotes the identity matrix of dimension 4, that is, (i) Stability of the infection-free equilibrium E 0 . The characteristic equation at E 0 is given by Assume that τ = 0. In this case, the equation (6) becomes We need to prove that all the roots of the characteristic equation have negative real parts. It is easy to see that y 1 = −n and y 2 = −m are roots of equation (7) and both are real negative roots. Thus, we just need to consider the third term of the above equation. Let Using the Routh-Hurwitz criterion, we know that all roots of p(y) have negative real parts if and only if the coefficients a i of p(y) are strictly positive. In our case, Hence, if R 0 < 1, then all roots of the characteristic equation (7) have negative real parts. Therefore, E 0 is locally asymptotically stable for τ = 0. Suppose now that τ > 0. To prove the stability of E 0 we use Rouché theorem, so we need to prove that all the roots of the characteristic equation (6) cannot intersect the imaginary axis, i.e., the characteristic equation cannot have pure imaginary roots. Suppose the reverse, i.e., that there exists w ∈ R such that y = wi is a solution of (6). Replacing y in the third term of (6), we get Note that we do not need to consider the full equation (6) because we already know that the remaining part of this equation has just two real negative solutions. By using the Euler formula and separating the real and imaginary parts of the above equation, we obtain that By adding up the squares of both equations and using the fundamental trigonometric formula, we obtain that which is the same as and equivalent to Hence, we have w 2 < 0, which is a contradiction. Therefore, we proved that if R 0 < 1, then the characteristic equation (6) cannot have pure imaginary roots and the infection-free equilibrium E 0 is locally asymptotically stable for any strictly positive time-delay. Suppose now that R 0 > 1. We know that the characteristic equation (6) has two real negative roots y = −n and y = −m. Thus, we need to check if the remaining roots of have negative real parts. It is easy to see that q(0) = uv − kλr m < 0, because we are assuming R 0 > 1. On the other hand, lim y→+∞ q(y) = +∞. Therefore, by continuity of q(y), there is at least one positive root of the characteristic equation (6). Hence, we conclude that E 0 is unstable. Finally, we need to analyse the case R 0 = 1, i.e., muv = kλr. In this case the characteristic equation (6) becomes To prove the stability, we need to check again if all the roots of the above equation have negative real parts. Note that y = 0, y = −n and y = −m are solutions of this equation, so we just need to prove that the remaining roots cannot have nonnegative real parts. Assuming that y = a + bi with a ≥ 0 is a solution of the above equation, then By using the Euler formula and by separating the real and imaginary parts, we get Adding up the squares of both equations and using the fundamental trigonometric formula, we obtain which is a contradiction because This proves that 0 is the unique root of (8) that does not have negative real part.
After some basic simplifications, we haveĪa−n < 0. Hence, if R 0 > 1+R 1 , then the characteristic equation (9) has a positive root and, consequently, the equilibrium E 1 is not locally asymptotically stable. On the other hand, if 1 < R 0 < 1 + R 1 , then y =Īa − n is a real negative root of the characteristic equation (9) and we just need to analyze the equation Consider τ = 0. From equation (10) we have where A = m+u+v+V r > 0, B−D =V ru+V rv+mu+mv > 0, C−E =V ruv > 0 and A(B − D) > (C − E). Therefore, from the Routh-Hurwitz criterion, it follows that all roots of (11) have negative real parts. Hence, E 1 is locally asymptotically stable for τ = 0. Let τ > 0. Suppose that (10) has pure imaginary roots, wi. By replacing y in (10) by wi, we get If we separate the real and imaginary parts, then we obtain −Aw 2 + C = E cos(wτ ) + Dw sin(wτ ), −w 3 + Bw = Dw cos(wτ ) − E sin(wτ ).
By adding up the squares of both equations, and using the fundamental trigonometric formula, we obtain that which is equivalent to we have that the left hand-side of equation (11) is strictly positive, which implies that this equation is not possible. Therefore, (9) does not have imaginary roots, which implies that E 1 is locally asymptotically stable for any time delay τ ≥ 0.
(iii) Stability of CTL-activated infection equilibrium E 2 . Assume R 0 > 1 + R 1 . The characteristic equation (5) at Suppose that there is a wi, w ∈ R, such that y = wi is root of equation (12). Then, (13) Since wi + m +V r 2 > |wi + m| 2 and |v + wi| 2 > |v| 2 , it follows from (13) that We conclude that the left hand-side of (13) is always strictly greater than the right hand-side, which implies that this equation is impossible. Hence, the solutions of the characteristic equation (12) cannot be pure imaginary. Therefore, by Rouchè theorem, E 2 is locally asymptotically stable for any time-delay τ ≥ 0.
3. Optimal control of the HIV model with intracellular and pharmacological delays. In the human system, RNA molecules are produced from DNA. Nevertheless, there are enzymes that make the reverse process, i.e., they can obtain DNA molecules from RNA. Such an enzyme is called a reverse transcriptase. One kind of such enzymes are found in HIV-1. As a result, when a virus particle infects a T-cell, it comes into the kernel of the cell and makes the reverse transcriptase process converting the RNA viral molecules into DNA viral molecules, which are then combined with DNA molecules of the CTLs. Hence, CTLs work to create new viruses instead of doing the defense job they are supposed to do in the immune system. Nowadays, there are drugs that can inhibit the reverse transcriptase, which allow the CTLs to keep their natural work. In this section, we formulate an optimal control problem for HIV-1 infection, with time delay in state and control variables, and derive extremals for the minimization of virus by the use of drugs that inhibit the reverse transcriptase of CTLs. We introduce a control function c(t) in model (1), t ∈ [0, t f ], that represents the efficiency of the reverse transcriptase inhibitors, which block a new infection. Due to the importance of the pharmacological delay in the HIV treatment, we consider a discrete time delay in the control variable c(t), denoted by ξ, which represents the delay that occurs between the administration of a drug and its appearance within the cells, due to the time required for drug absorption, distribution, and penetration into the target cells [27]. We propose the following control system with discrete time delay in the state and control variables: The initial conditions for the state variables I and T and, due to the delays, initial functions for the state variables Z and V and control c, are given by We note that values (15) are the only ones that are considered in our numerical simulations (Section 4). The control function c(t) is bounded between 0 and 1. If it takes the value 0, then the drug therapy for the transcriptase reversion has no efficacy. If the control takes the value 1, then it will be 100% effective. Precisely, we consider following set of admissible control functions: We consider the L 1 objective functional which measures the concentration of virus and the treatment costs for the period of time under study. The optimal control problem consists in determining a control function c(·) ∈ L 1 ([0, t f ], R) that minimizes the cost functional (17) subject to the control system (14), initial conditions (15) and control constraints (16). In Section 4.2, we present numerical solutions for three cases of delays τ and ξ and weights w = 1 and w = 5. To apply the optimality conditions given by the Minimum Principle for Multiple Delayed Optimal Control Problems of [10, Theorem 3.1], we introduce the delayed state variables ζ(t) = Z(t−τ ), η(t) = V (t−τ ) and the control variable ω(t) = c(t − ξ). Using the adjoint variable λ = (λ Z , λ I , λ V , λ T ) ∈ R 4 , the Hamiltonian for the cost functional (17) and the control system (14) is given by The adjoint equations are given by where the subscripts denote partial derivatives and χ [0,t f −τ ] is the characteristic function in the interval [0, t f − τ ] (see [10]). Since the terminal state is free, i.e., (Z(t f ), I(t f ), V (t f ), T (t f )) ∈ R 4 , the transversality conditions are

OPTIMAL CONTROL OF A DELAYED HIV MODEL 9
To characterize the optimal control c, we introduce the following switching function: The minimality condition of the Minimum Principle [10, Theorem 3.1] gives the control law Similar arguments can also be used to solve related optimal control problems, e.g., one may consider an additional constraint on the final virus concentration or inclusion of the final values of this concentration in the cost functional.

Numerical simulations.
In this section, we study numerically the stability of the delayed model (1) proposed in Section 2 and the solution of the optimal control problem proposed in Section 3. We consider the initial conditions (15) and the parameter values as given in Table 1, which are based on [12]. of system (1) is locally asymptotically stable for any time delay τ ≥ 0. In Figure 1, we observe the stability of system (1) in a time interval of 500 days and a time delay of 0.5 days (τ = 0.5). In Figure 2, we compare the behavior of system (1) for τ = 0 (no delay) and delay τ = 0.5. The first local maximum of concentration of infected cells, virus and CTLs is smaller in the delayed case (τ = 0.5). The local maxima are  Table 1 and time delay τ = 0.5.  similar, although they are attained at latter in the delayed case, when compared to the nondelayed situation. At the end of 50 days, the values of the variables Z(t), I(t), V (t) and T (t) are similar in delayed and nondelayed cases.

4.2.
Optimal control problem with state and control delays. In this section, we present numerical solutions to the delayed optimal control problem (14)- (17) in the time interval [0, 50] days and consider three cases: Case 1: τ = ξ = 0 (no delays); Case 2: τ = 0.5, ξ = 0 (intracellular delay τ only); Case 3: τ = 0.5, ξ = 0.2 (intracellular delay τ and pharmacological delay ξ). As before, we consider the parameter values from Table 1 and the weight parameters w = 1 and w = 5 in the cost functional (17). To solve the delayed optimal control problem (14)-(17), we discretize the control problem on a sufficiently fine grid [10] and obtain a nonlinear optimization problem (NLP). The NLP is implemented using the Applied Modeling Programming Language AMPL [9], which can be interfaced with several large-scale nonlinear optimization solvers like the interior-point solver Ipopt; see [35]. We mostly use N = 2500 grid nodes and the trapezoidal rule as integration method to compute the solution with an error tolerance of eps = 10 −9 . In all three cases, the computed controls are bang-bang with only one switch at t s : For the weight w = 1, we obtain the following numerical results: A zoom into the controls and switching functions, in a neighborhood of the switching time t s , is displayed in Figure 3. The state trajectories in the three cases are very similar on the terminal time interval [15,50], while the concentration of uninfected cells Z(t) is nearly identical on the whole time interval [0,50]. To display the effect of the delays on the state variables I, V, T , Figure 4 shows a comparison of the state trajectories in Case 1 (no delays) and Case 3 (state and control delays). We see that the delay in the control c implies an increase of the concentration of infected cells I(t) in the first two days (the delay on the drug effect), which is also responsible for an increase on the concentration of the free virus particles V and CTL cells T .  The bang-bang controls and the switching functions in Figure 3 do not only match the switching condition (19) but satisfy also the so-called strict bang-bang property [24] with respect to the Minimum Principle: The strict bang-bang property enables us to check second-order sufficient conditions (SSC) for the bang-bang control in the non-delayed Case 1. In the delayed Cases 2 and 3, no sufficient conditions are available in the literature. In Case 1, we consider the so-called Induced Optimization Problem (IOP), where the switching time t s in (20) is the only optimization variable. Hence, we optimize the function J(t s ) = J(c) with respect to t s . The IOP can be solved using the arc-parametrization method [17,24] and its implementation in the optimal control package NUDOCCCS [2]. We obtain the highly accurate numerical results J(c) = 475.1854, t s = 47.0903, J (t s ) = 5.0962 > 0.
In view of the strict bang-bang property (21) and the positive second derivative J (t s ) = 5.0962, we conclude from Theorem 7.10 in [24] that the bang-bang control in Case 1 provides a strict strong minimum.
Since SSC hold, it follows from the standard sensitivity result in finite-dimensional optimization [7] (cf. also [3]) that the switching time t 1 is locally a C 1 -function with respect to all parameters p in the system. The state trajectories are locally C 1 -functions except at the switching time t s . The code NUDOCCCS [2] allows to compute the sensitivity derivatives dt s /dp and dy(50)/dp for y ∈ {J(c), Z, I, V } at a nominal parameter value p 0 . The sensitivities dT (50)/dp are very small so that we do not list them here. Choosing, e.g., the parameter p ∈ {w, r, v}, where w is the weight parameter in the functional (17) and the parameters r and v are as in Table 1, we get the following sensitivity derivatives at their nominal values w 0 = 1, r 0 = 0.0014 and v 0 = 1: parameter dt 1 /dp dJ(c)/dp dZ(50)/dp dI(50)/dp dV (50)/dp The sensitivity derivatives quantify our more intuitive feeling on how the switching time changes under parameter perturbations. As an example, let us increase the weight parameter w = 1 + ∆ for the control in the objective J(c) (17). Then the switching time t s (w) decreases and has the approximative value Similar Taylor expansion approximations hold for the other quantities. It is an interesting exercise to show that the sensitivity derivative dJ(c)/dp agrees with t s = 47.09. Finally, Figure 5 displays a comparison of the controlled state variables with the uncontrolled ones in Case 2.  The computed switching functions (18) match the control law (19) and satisfy the strict bang-bang property (21). SSC can only be verified in the non-delayed Case 1. NUDOCCCS computes the second derivative J (t s ) = 25.43 > 0. Hence, the bang-bang control provides a strict strong minimum in view of [24, Theorem 7.10].

5.
Conclusion and discussion. In this paper we have considered not only intracellular delay (delay τ in the state variables) as done in the literature [12], but also a pharmacological delay (delay ξ in the control function). The pharmacological delay causes an increase of the concentration of the infected cells in an initial interval of time. However, after this increase, related to the delay in the action of the drugs in the cells, the concentrations of infected cells, virus and CTL cells associated to the extremal solution of the optimal control problem with both delays in state and control variables, decrease significantly. The extremal control is bang-bang and switches from its maximal value one to zero. This type of control is easier to implement, from a medical point of view, when compared to controls found in [12] for L 2 cost functionals. We observe that the extremal control derived from the application of Pontryagin's necessary optimality condition [10, Theorem 3.1] to our multiple delayed optimal control problem, is associated to a marked reduction of the concentration of infected cells, virus and CTLs, as well as treatment costs, and to an increase of the uninfected target cells. Sufficient optimality conditions could only be checked for the non-delayed solution in Case 1. In this case, we could also perform a local sensitivity analysis by computing the sensitivity derivatives. It remains an open and challenging question to prove and verify sufficient optimality conditions for delayed bang-bang controls.