An Age-structured Model with Immune Response of Hiv Infection: Modeling and Optimal Control Approach

This paper develops and analyzes an age-structured model of HIV infection with various compartments, including target cells, infected cells, viral loads and immune effector cells, to provide a better understanding of the interaction between HIV and the immune system. We show that the proposed model has one uninfected steady state and several infected steady states. We conduct a local stability analysis of these steady states by using a generalized Jacobian matrix method in conjunction with the Laplace transform. In addition, we consider various techniques and ideas from optimal control theory to derive optimal therapy protocols by using two types of dynamic treatment methods representing reverse transcriptase inhibitors and protease inhibitors. We derive the necessary conditions (an optimality system) for optimal control functions by considering the first variations of the Lagrangian. Further, we obtain optimal therapy protocols by solving a large optimality system of equations through the use of a difference scheme based on the Runge-Kutta method. The results of numerical simulations indicate that the optimal therapy protocols can facilitate long-term control of HIV through a strong immune response after the discontinuation of the therapy.


Introduction
Substantial progress has been made in the mathematical modeling of HIV infection and treatment strategies [2,20,21,24,26,29,30].Existing models explore mainly the dynamics of CD4+ T helper cells, virus production/clearance, and the effects of various anti-retroviral drug treatments.In general, during primary HIV infection, the viral load in plasma increases, peaks, and then declines within several weeks.Stafford et al. [30] address this decline by developing various models of primary HIV infection and comparing predictions from the models with data from patients.Callaway and Perelson [8] investigate the ability of several models of HIV infection to explain how low viral loads can be sustained.Perelson and Nelson [24] explain how mathematical modeling can facilitate the discovery of some important features of HIV pathogenesis and impact the way in which HIV patients are treated with efficient anti-retroviral drugs.Nowak et al. [23] provide some analytic solutions for the emergence of a resistant virus under a single-drug therapy.Herz et al. [13] introduce a delay model for considering the time between the infection of target cells and the production of virions.Subsequently, Culshaw and Ruan [9] examine the effects of a delay on the stability of the endemical equilibrium.
Age-structured models have received considerable attention from researchers interested in exploring the epidemiology of HIV.Thieme and Castillo-Chavez [31] examine the effects of infectionage-dependent infectivity on the dynamics of HIV transmission in a homogeneously mixing population.Kirschner and Webb [16] propose an age-structured model to account for the mechanism underlying AIDS chemotherapy.Recent studies have developed age-structured models for analyzing the within-host dynamics of HIV to provide a better understanding of the interaction between the immune system and HIV [22,27].Nelson et al. [22] propose and analyze an age-structured model that allows for variations in the death rate for infected T cells and the production rate for virions.Rong et al. [27] extend the age-structured model by incorporating combination therapies and considering the effects of anti-retroviral therapy on viral dynamics.
Previous studies have demonstrated that the long-term use of a single anti-retroviral drug can lead to the emergence of a resistant virus [18,23].Thus, the most widely used strategy for HIV patients is highly active anti-retroviral therapy (HAART), which uses one or more reverse transcriptase inhibitors (RTIs) and a protease inhibitor (PI).However, although HAART can successfully reduce and maintain low viral loads in many patients, its prolonged use entails severe complications.In addition, HAART can be extremely costly in developing countries.In this regard, a number of researchers have searched for optimal treatment strategies that can reduce drug side effects, virus mutations, and complex/expensive medication burdens.There is no general consensus on optimal treatment strategies or interruption schemes.One way to consider optimal treatment schedules is to use a mathematical model for HIV infection in conjunction with control theory.Some studies have suggested continuous optimal treatment schedules based on open-loop control [2,3,10,17], and others have used feedback control laws chosen on a real-time basis through the observation of the full or partial state of the system [4,7,28].More recently, Jang et al. [14] consider the free terminal time optimal control problem to derive the minimum treatment duration as well as the optimal multidrug treatment method.Banks et al. [5] explore some problems associated with optimal feedback control and state estimators for HIV infection by considering the state-dependent Riccati equation (SDRE) approach.Banks et al. [6] apply a receding horizon observer to an HIV feedback control problem with a long measurement time to derive optimal treatment methods for HIV progression.Kwon et al. [19] investigate an optimal control problem for an age-structured model with three compartments including target cells, infected cells, and viral loads.
This paper develops an age-structured model with the immune response of HIV infection to derive an optimal multidrug therapy.The proposed model extends the age-structured models in [22,27] by incorporating immune effector cells (cytotoxic T lymphocytes or CTLs) and combining two types of therapies (RTIs and PIs).In addition, we investigate the steady states of the model with no drug therapy and their local stability.The analytical results indicate that the proposed model has one uninfected steady state that is unstable and two infected steady states that are stable.In addition, the numerical results demonstrate that we can determine an optimal treatment scheme in which a patient moves to an immune-dominant state after discontinuing his or her therapy.
The rest of this paper is organized as follows: Section 2 formulates and analyzes an agestructured HIV model that generalizes the model in [22] by incorporating immune effector cells.Section 3 presents the formulation for the optimal control problem and derives the corresponding optimality system by considering the first variations of the Lagrangian.Section 4 introduces an algorithm based on the gradient method to solve an optimality system and presents the results of numerical simulations, demonstrating the effectiveness of the continuous optimal therapy, and Section 5 concludes.The Appendix provides detailed calculations.

Age-Structured Model
We now introduce an age-structured model of HIV infection.The model has four state variables: uninfected CD4+ T cells T (t); infected CD4+ T cells structured by the age a of their infection, T * (a, t); virions V (t); and immune effector cells E(t).A system of three ordinary differential equations and one first order hyperbolic equation describing HIV dynamics is given by This model assumes that uninfected CD4+ T cells are produced at the constant rate λ and die at the rate d.Here kV T , where k denotes the infection rate, indicates the infection process in which infected cells T * are produced through interactions between uninfected target cells T and virions V .The death rate δ(a) and the virion production rate P (a) for T * are assumed to be functions of the age of cellular infection, a, and virions V are assumed to be cleared at the constant rate c.In addition, we assume that da dt = 1; that is, the time unit for the age of infection is the same as that for the clock time.Because the model includes a first order hyperbolic equation, we need to introduce the boundary and initial conditions.
Infected CD4+ T cells of age zero are created by infection, that is, In addition, we impose specific initial conditions T (0) = T 0 , T * (a, 0) = T * 0 (a), V (0) = V 0 and E(0) = E 0 .The control term ǫ 1 (t) represents the effectiveness of RTIs, which block new infection.Thus, the infection rate k is reduced to (1 − ǫ 1 (t))k, where 0 ≤ a 1 ≤ ǫ 1 (t) ≤ b 1 < 1.The control term ǫ 2 (t) represents the effectiveness of PIs, which reduce the number of infectious virions.Thus, the production rate P (a) is reduced to (1 − ǫ 2 (t))P (a), where 0 ≤ a 2 ≤ ǫ 2 (t) ≤ b 2 < 1.Here a i and b i , (i = 1, 2) represent the minimum and maximum drug efficacy, respectively.Infected CD4+ T cells are cleared by the action of immune effector cells E. This action can be expressed as m(a)T * (a, t)E(t).The articles [2,3]  With the above boundary and initial conditions and a smooth enough control function, there exists a unique solution to system (2.1) that remains bounded and nonnegative for t > 0 (see [19,22,32]).
Model (2.1) contains several constant parameters and function parameters that must be assigned for numerical simulations.Table 1 provides a summary of the descriptions and numerical values for these parameters, which are extracted mainly from [2,3,5].
It is reasonable to expect that a ≤ a max so that the integral in (2.1) is not necessarily an infinite integral.The virus production kernel P (a) has the maximum production rate P max because cellular resources ultimately limit how rapidly virions can be produced.

Analysis of the Age-Structured Model
with the initial condition Here (T e , T * e (a), V e , E e ) is the steady state of the age-structured model.We can easily find a trivial or uninfected steady state: To consider other steady states, we need the following new notations: Solving the second equation in (2.3) with the initial condition, we have The third equation in (2.3) implies that Thus, we get T e = c kN P . (2.7) The first equation in (2.3) implies (2.9) In addition, substituting (2.9) into the fourth equation in (2.3) and rearranging yields which implies that (2.10) Finally, we get the non-trivial or infected steady states and we determine E e by solving equation (2.10).
A general analysis of the proposed model's steady states and their local stability is challenging because of the high nonlinearities in equation (2.10).However, we are still interested in the ability of the model to demonstrate multiple locally asymptotically stable steady states, and thus, we calculate steady states numerically and investigate their stability by using a generalized Jacobian matrix method with given numerical values for the parameters.This process requires the method of characteristics and the Laplace transform of each age-dependent term (see [22] and the references therein).Now we need an explicit functional form of the virus production kernel P (a).According to previous research [22], we may consider one functional form of kernels that may capture the features of the biology, where the maximum production rate P max is provided because cellular resources ultimately limit how rapidly virions can be produced.We can define a delayed or non-delayed exponential function as follows: where β controls how rapidly the saturation level P max is reached and d 1 denotes the delay in virus production (that is, it takes some time d 1 after the initial infection for the first virions to be produced).For simplicity, we assume that the functions δ(a) and m(a) are constants for deriving steady states and considering their stability: δ(a) = δ = 0.7 and m(a) = m = 0.01.

Stability of the Steady States
In this section, we demonstrate that the model has multiple locally asymptotically stable steady states by calculating its steady states and conducting a stability analysis.The standard Jacobian matrix method is generally used to determine the stability of a system of ordinary differential equations.For an age-structured model with immune response, we employ a generalized Jacobian matrix method in conjunction with the Laplace transform for the stability analysis as follows: Note that we can derive a general solution to the second equation in (2.1): where T * 0 (a) is the initial value and B(t) = kV (t)T (t) is the boundary value of T * (a, t) because we assume that ǫ 1 = ǫ 2 = 0 in this section.Substituting (2.13) into the system (2.1), we can rewrite the system (2.1) as (2.14) Since the functions P (a) , T * 0 , and E are bounded above, we have Once the stability of (2.14) is found, it then gives the stability of (2.1).Thus, we now consider where the notations N P and N 1 are defined in (2.5) by assuming δ(a) = δ and m(a We assume that the perturbations ( T (t), V (t), B(t), E(t)) are sufficiently small for t ≥ 0 and that they are zero for t < 0. By linearizing the exponential function (e x ≈ 1 + x), we obtain Similarly, we have where For sufficiently small perturbations, that is, B and E ≪ 1, we may assume that | 1.By (2.18) and linearizing the function ( (2.20) We now consider the Laplace transform (L) of (2.20) by denoting transforms with hats.By Fubini's theorem and changing variables, we obtain In addition, we have and The Laplace transform of (2.20) now form the following linear system: and We can solve the linear system by using Cramer's rule to determine

State
for some functions f T , f V , f B and f E .Here A(s) is a matrix that depends on the Laplace variable s in system (2.21).Taking the inverse Laplace transform, we determine the growth rate s of the solutions ( T (t), V (t), B(t), E(t)) at the poles of the solutions (2.22).After checking that there are no common factors in f T , f V , f B or f E , we find that these poles are at the zeroes of det(A(s)).Therefore, we conclude that if the zeroes of det(A(s)) all have negative real parts, then the equilibrium (T e , T * e , V e , E e ) is locally asymptotically stable.Given the specified parameters in Table 1 and by assuming P max = 80, β = 5, and d 1 = 0 in equation (2.12), we determine that model (2.1) has four physical steady states and several nonphysical steady states of which at least one state variable is negative.Here we omit the nonphysical steady states.Table 2 provides a summary of the numerical values and stability results for these steady states.By the above argument, we can easily check the trivial or uninfected steady state EQ 0 is locally unstable.In addition, we can show that there are one locally unstable equilibrium and two locally asymptotically stable equilibria for an infected patient with no treatment (ǫ 1 = ǫ 2 = 0).In fact, the states of the system converge to an "unhealthy" steady state EQ 1 in which target cells are substantially depleted and the viral load is extremely high when a small amount of virus is introduced in the uninfected steady state EQ 0 (see Figure 1).In addition, the age-structured model has a "healthy" stable steady state EQ 3 in which a strong immune response develops for successful control of viral infection.One of the objectives of this paper is to derive optimal strategies for moving the model system to EQ 3 , not to EQ 1 , as shown in Figure 1.

A B
Figure 1: EQ 0 indicates the uninfected locally unstable equilibrium; EQ 1 , the "unhealthy" locally asymptotically stable equilibrium with attraction region A; and EQ 3 , the "healthy" locally asymptotically stable equilibrium with attraction region B. The dashed line indicates an uncontrolled trajectory, and the solid line, a controlled trajectory.

An Optimal Control Problem
In this section, we formulate an optimal control problem based on model (2.1) to derive optimal treatment strategies in which HIV patients reach the "healthy" steady state EQ 3 (i.e., immune control of viral infection).We first define the objective functional as follows: where R 1 , R 2 , Q, and S are the weight constants for controls, viral loads, and immune effector cells, respectively.The first two terms represent the systemic cost of drug treatment.This cost is based on the actual treatment cost as well as the severity of unintended side effects of drugs.
Here the objective is to minimize both the systemic cost of drug treatment and the viral load while maximizing the number of immune effector cells to achieve a "healthy" steady state EQ 3 .The two control functions ǫ 1 (t) and ǫ 2 (t) represent the efficacy of RTIs and PIs satisfying 0 , respectively.The case ǫ i (t) = b i , (i = 1, 2) represents maximum efficacy of drugs.We consider the following optimal control problem: Problem 3.1.We seek optimal controls ǫ * 1 and ǫ * 2 such that subject to a system of the partial and ordinary differential equations (2.1) with the boundary condition T * (0, t) = (1 − ǫ 1 (t))kV (t)T (t) and the initial conditions The basic framework of an optimal control problem is to prove the existence of an optimal control and then characterizes the optimal control by using the optimality system.Remark 3.2.We note that the existence of optimal control functions to Problem 3.1 can be shown by using the same arguments as in Theorem 3.2 in [19].
Theorem 3.3.Given optimal controls and solutions to the corresponding constraint equations (2.1) that minimize the objective functional (3.1), there exists an adjoint variable Y = (ξ 1 (t), ξ 2 (a, t), ξ 3 (t), ξ 4 (t)) with the terminal conditions ξ 1 (T f ) = 0, ξ 2 (a, T f ) = 0, ξ 3 (T f ) = 0, and ξ 4 (T f ) = 0 and the boundary condition lim a→∞ ξ 2 (a, t) = 0.In addition, the optimal control functions ǫ * 1 and ǫ * 2 are given by (3.3) Proof.We apply the necessary conditions for finding the stationary points of L. Setting to zero the first variations of the Lagrangian with respect to Y yields the constraints (2.1).Setting to zero the first variations with respect to X yields the adjoint equations (3.2):The first adjoint equation can be obtained by taking ∂L ∂T = 0.
Since variations T in the state are arbitrary, we get with the terminal condition ξ 1 (T f ) = 0. Using similar arguments, the first-order necessary condition ∂L ∂T * = 0 yields the equation for ξ 2 in the adjoint system (3.2).

∂L ∂T
Since variations T * in the state are arbitrary, we have the terminal condition ξ 2 (a, T f ) = 0, and the boundary condition lim a→∞ ξ 2 (a, t) = 0 .
The adjoint equations for ξ 3 and ξ 4 in (3.2) may also be derived from for any arbitrary variations V and E in the state.Finally, to derive the optimality conditions, we must consider the controls that adhere to the bounds.Because of the bounds on the controls, the derivatives of the objective functional may not be zero at the optimal control.However, the idea is essentially the same as that in the case without bounds, and we can derive explicit expressions for the optimal control functions (3.3) by taking the approach in [2,3,19].
To get the explicit formulation of the optimal controls ǫ * 1 and ǫ * 2 , we explore the necessary optimality conditions These numerical results illustrate that it is possible to move a patient from a state of early infection to a "healthy" stable steady state (EQ 3 ) in which a strong immune response successfully controls the viral load and consequently restores uninfected target cells without requiring drugs.This suggests one possible scenario in which optimal treatment strategies facilitate long-term control of HIV after the discontinuation of the therapy.

Conclusions
This paper formulates and analyzes an age-structured model of HIV infection with various compartments (including target cells, infected cells, viral loads, and immune effector cells) that is subject to multiple drug treatments as control functions.The proposed model possesses an uninfected equilibrium and several infected equilibria.We conduct a detailed stability analysis of the steady states by using a generalized Jacobian matrix method in conjunction with the Laplace transform.We then apply techniques and ideas from optimal control theory to determine optimal therapy protocols.We derive an optimality system from which optimal solutions can be determined.We define finite difference approximations based on the Runge-Kutta method and implement them in conjunction with a gradient method to solve the optimality system.We design optimal treatment  strategies to move the state of the system to a "healthy" steady state with low viral loads and high levels of immune effector cells.In addition, the results of the numerical simulations indicate that the optimal therapy protocols can facilitate long-term control of HIV by inducing a strong immune response after the discontinuation of the therapy.

0 P
With no drug treatment (ǫ 1 = ǫ 2 = 0), model (2.1) has several steady states.To determine the steady states, we solve the equations 0 = λ − dT e − kT e V e dT * e da = −δ(a)T * e (a) − m(a)E e T * e (a)T * e (a)da − cV e 0

(2. 15 )
Let (T e , V e , B e , E e ) be the solutions to 0 = λ − dT e − kV e T e 0 = B e N P − cV e 0 = kV e T e − B e 0
d 0.5 (cells/µl) Saturation constant for the death of immune effector cells δ E 0.1 (1/day) Natural death rate of immune effector cells Table 1: Descriptions and numerical values of the parameters in the HIV model Remark 2.1.

Table 2 :
Numerical values and stability results for steady states of the age-structured model with no treatment