Stability and Feedback Stabilization of HIV Infection Model with Two Classes of Target Cells

We study the stability and feedback stabilization of the uninfected steady state of a human immunodeficiency virus HIV infection model. The model is a 6-dimensional nonlinear ODEs that describes the interaction of the HIV with two classes of target cells, CD4 T cells and macrophages, and takes into account the Cytotoxic T Lymphocytes CTLs immune response. Lyapunov function is constructed to establish the global asymptotic stability of the uninfected steady state of the model. In a control system framework, the HIV infection model incorporating the effect of Highly Active AntiRetroviral Therapy HAART is considered as a nonlinear control system with drug dose as control input. We developed treatment schedules for HIV-infected patients by using Model Predictive Control MPCbased method. The MPC is constructed on the basis of an approximate discrete-timemodel of the HIV infection model. TheMPC is applied to the stabilization of the uninfected steady state of the HIV infection model. Besides model inaccuracies that HIV infection model suffers from, some disturbances/uncertainties from different sources may arise in the modelling. In this work the disturbances are modelled in the HIV infection model as additive bounded disturbances. The robustness of the MPC against small model uncertainties or disturbances is also shown.


Introduction
In the last decade many mathematical models have been proposed to describe the immunological response to infection with the human immunodeficiency virus HIV .HIV is responsible of acquired immunodeficiency syndrome AIDS .Some of these models mainly represent the interaction of the HIV with the CD4 T cells, others take into account the Cytotoxic T Lymphocytes CTLs immune response 1 .A tremendous effort has been made in studying the basic properties of these models such as positive invariance properties, boundedness of the model solutions, and stability analysis of the steady states see e.g., 2-14 .
Studying such properties is important for understanding the associated characteristics of the HIV dynamics.The treatment of HIV-infected patients is of major importance in today's social medicine.Currently, the most important categories of anti-HIV drugs are reverse transcriptase inhibitors RTIs drugs and protease inhibitors PIs drugs.Reverse transcriptase inhibitors prevent the HIV from infecting cells by blocking the integration of the HIV viral code into the host cell genome.Protease inhibitors prevent already infected host cells from producing infectious virus particles.Recently, Highly Active AntiRetroviral Therapies HAARTs which consist of one or more RTI and a PI can suppress viral load below detectable levels and consequently prolong time to the onset of AIDS.Perleson et al., observed that after the rapid first phase of decay during the initial 1-2 weeks of antiretroviral treatment, plasma virus levels declined at a considerably slower rate 15 .This second phase of viral decay was attributed to the turnover of a longer-lived virus reservoir of infected cells.These cells are called macrophages and considered as the second-target cell for the HIV.Therefore, the HIV infection model with two classes of target cells is more accurate than the model with one class of target cells see 8, 16 .Some HIV infection models exist to describe the interaction process of the HIV not only with the CD4 T cells but also with the macrophages which are the crucial immune responses and play important roles in phagocytosis see e.g., 8, 16-19 .The basic properties of the two target cells models are recently studied in 3, 20-22 .In 23, 24 we have studied a class of virus infection models assuming that the virus attacks multiple classes of target cells.However, in 3, 20-24 , the CTL immune response is neglected.The role of CTL cells is to attack the infected cells.The first purpose of the present paper is to study the basic properties of an HIV infection model which describes the interaction of the HIV with two target cells, CD4 T cells and macrophages and takes into consideration the CTL immune response.The global stability of the uninfected steady state of the model is established using a Lyapunov approach.
Optimal treatment scheduling of HIV infection using a control theoretic approach is the subject of substantial research activity.In 17,18,[25][26][27][28] , an open-loop optimal controller is designed via Pontryagin's Maximum Principle.A major drawback of such controller is its lack of robustness against disturbances or model uncertainties.In fact, the interaction of HIV with immune system is not very clear and complicated which leads to model inaccuracies and parameter uncertainties.Moreover, some disturbances may arise from immune system fluctuation or immune system reactions due to coinfections.Therefore, the design of optimal treatment schedules based on an open-loop optimal controller may lead to undesired results.To overcome this drawback, a feedback controller which has inherent robustness property against disturbances has to be designed.Recently, model predictive control MPC method is developed for determining optimal treatment schedules for HIV patients 20, 29-33 .The MPC method obtains the feedback control by solving a finite horizon optimal control problem at each time instant using the current state of the system as the initial state for the optimization and applying "the first part" of the optimal control.The study of stability and robustness properties of MPC schemes has been the subject of intensive research in recent years see e.g., 34,35 . In 20,30,31,33 , the HIV infection model is discretized and the MPC is constructed on the basis of an approximate discrete-time model.Sufficient conditions are established which guarantee that the MPC designed via approximate discrete-time model stabilizes the original continuous-time model.In 29, 32 , the effect of the discretization of the differential equations on the stability analysis was not considered.The importance of approximate discrete-time design is supported by a series of counter examples see e.g., 36 , which show that the controller that stabilizes the approximate discrete-time model, may fail to stabilize the original continuous-time system.In 29-33 , the MPC is applied to the one target cell i.e., CD4 T cell model.The HIV infection model with two target cells is considered in 20 , but the CTL immune response is neglected.
The second purpose of the present paper is to develop treatment schedules for HIVinfected patients by using MPC.The construction of MPC is based on the approximate discrete-time model of the model.The HIV infection model with additive disturbances is also considered, and the inherent robustness properties of the MPC are established.
The layout of the paper is as follows: in Section 2, we introduce the HIV infection model and study its basic properties.In Section 3, the HIV infection model with additive disturbance is outlined.In Section 4 we outline the MPC design for nonlinear control systems and summarize the main results obtained in 36, 37 .Application of MPC to the HIV infection model is given in Section 5. Section 6 presents the simulation results.The last section is the conclusion.

HIV Infection Model
We will study the mathematical model of HIV infection which describes the interaction of HIV with two cocirculation populations of target cells, potentially representing CD4 T cells and macrophages and takes into consideration the CTL immune response.The model is a modification of the HIV infection model presented in 3, 8, 16, 20 , which include an additional state variable for the CTL cells: Here x t , x 1 t , y t , y 1 t , v t , and z t represent the concentrations of uninfected CD4 T cells, infected CD4 T cells, uninfected macrophages, infected macrophages, free virus particles, and CTL cells, respectively, at time t.The populations of the uninfected CD4 T cells and macrophages are described by 2.1 and 2.3 , respectively, where λ 1 and λ 2 represent, respectively, the rates of which new CD4 T cell and macrophages are generated from sources within the body, d 1 , d 2 are the death rate constants, and β 1 , β 2 are the infection rate constants.Equations 2.2 and 2.4 describe the population dynamics of the infected CD4 T cells and macrophages and show that they die with rate constant a and killed at rate px 1 z and py 1 z.The virus particles are produced by the infected CD4 T cells and infected macrophages with rate constant k and are cleared from plasma with rate constant r.The CTL cells are produced at a rate c x 1 y 1 z and are decayed at a rate bz.The effect of the RTI and PI drugs is represented by the chemotherapy functions e −ψ 1 m 1 t and e −ψ 2 m 2 t , where the parameters ψ 1 and ψ 2 are the efficiencies of RTI and PI drugs, respectively.The variables m 1 t and m 2 t are the doses from drugs administrated to RTI and PI, respectively see 38 .All the parameters of the model are supposed to be positive.
We introduce control functions as u i t ψ i m i t , i 1, 2. Then system 2.1 -2.6 can be considered as a nonlinear control system with x, x 1 , y, y 1 , v, z as the state vector and u 1 , u 2 is the control input vector.We are now ready to present a study on the basic mathematical properties of the model.

Positive Invariance
Now we show that model 2.1 -2.6 is biologically acceptable in the sense that no population goes negative:

2.7
This means that the nonnegative orthant R 6 is positively invariant, namely, if a trajectory starts in the nonnegative orthant, it remains there.The boundedness of the solutions of model 2.1 -2.6 will be given in the following proposition.

Proposition 2.1.
There exist positive numbers L 1 , L 2 , and L 3 such that the compact set, where

Steady States
We will compute the steady states of system 2.1 -2.6 under constant controllers, that is, for u j t u j , j 1, 2, t ≥ 0. It will be explained in the following that the behavior of model 2.1 -2.6 crucially depends on the basic reproduction number given by where We note that R 0 can be written as where 12 are the basic reproduction numbers of each CD4 T and macrophages dynamics separately.
The immune strength of CD4 T and macrophages are given, respectively, by I 1 cλ 1 /ab and then there exist three steady states E 0 , E 1 , and E 2 .
Proof.The proof can follow the same lines as that of Proposition 2 in 3 and Theorem 2.2 in 39 .
In the present paper, we are interested in studying the stability of the uninfected steady state of the HIV infection model.

Global Stability of E 0
In this section, we prove the global stability of the uninfected steady state by using a Lyapunov approach.
Proof.By the method of 3, 5 , we consider a Lyapunov function:

2.13
We note that W is defined, continuous, and positive definite for all x, x 1 , y, y 1 , v, z > 0. Also, the global minimum W 0 occurs at the uninfected steady state E 0 .Further, function W along the trajectories of 2.1 -2.6 satisfies

2.14
Since I 1 > 0, then the last term of 2.14 is negative.Also, since the arithmetical mean is greater than or equal to the geometrical mean, then the first two terms of 2.14 are less than or equal to zero.Therefore, if then R 0 u 1 , u 2 < 1, therefore the corresponding trajectory will tend to E 0 as t → ∞.
Remark 2.5.We observe that, if R 0 < 1, then it is sure that R 1 < 1 and R 2 < 1.But if one considers the three-dimensional model 2.1 -2.2 , and 2.5 and designs a controller such that R 1 < 1, then the whole system may be unstable around E 0 , because R 0 > 1.This shows the importance of considering the effect of the macrophages in the HIV dynamics.We note that the HIV infection model suffers from model inaccuracies or disturbances that may arise from different sources such as, modeling errors, immune system fluctuation, immune effect of a coinfection, measurement noise, and estimation errors.Therefore, in the next section we will consider the presence of the disturbances in the HIV infection model.

HIV Infection Model with Additive Disturbances
The model inaccuracies or disturbances can be modeled in the HIV infection model as additive disturbances see 20 .Therefore, model 2.1 -2.6 can be written as: In 3.1 -3.6 , w i t describe model uncertainties/disturbances which are assumed to satisfy the following bound: It is important to show that in the presence of the disturbances, the nonnegative orthant R 6 is positively invariant and the solutions of model 3.1 -3.6 are bounded.

Positive Invariance
For the model 3.1 -3.6 , the following conditions guarantee that the nonnegative orthant R 6 is positively invariant see 20 :

3.8
We note that the conditions in 3.8 give a lower bound of the disturbances only at the the boundary of R 6 .Proposition 3.1.Suppose that the disturbances satisfy the bound 3.7 then there exist positive numbers M 1 , M 2 , and M 3 such that the compact set, p/c 6 .

MPC for Nonlinear Systems
In this section, we outline the MPC design for nonlinear systems and give a review on the results obtained in 36, 37 .We have shown in the preceding sections that the HIV system states for the nominal model 2.1 -2.6 and disturbed model 3.1 -3.6 can be taken from compact sets.Moreover, since the drug dosage of HAART can not arbitrarily increased, thus the controller can also be taken from a compact set.
The sets of real and natural numbers including zero are denoted, respectively, by R and N. The notation R denote the set of real numbers in the interval 0, ∞ .A continuous function σ : R → R is of class-K if σ 0 0, σ s > 0 for all s > 0 and it is strictly τ is of class-K in s for every τ ≥ 0, it is strictly decreasing in τ for every s > 0 and β s, τ → 0 when τ → ∞.The Euclidean norm of a vector x is denoted as x .Given ρ > 0 we define B ρ to be a ball of radius ρ centered at the origin.We introduce the following notation Y ρ Y ∩ B ρ for any set Y.
Consider a continuous-time nonlinear control system given by ẋ t f x t , u t , x 0 x 0 , 4.1 where x t ∈ R n , u t ∈ U ⊂ R m are the state and control input, respectively, f : R n ×U → R n is continuous and Lipschitz continuous with respect to x in any compact set and f 0, 0 0, and U is compact and 0 ∈ U.
The control is taken to be a piecewise constant function where T > 0 is the control sampling period which is fixed.We will assume that there is a compact set X ⊂ R n containing the origin that is positively invariant with respect to system 4.1 for any piecewise constant controller u ∈ U. Let t → φ E t; x, u denote the solution of 4.1 with given u, and x x 0 .The exact discretetime model which describes the behavior of the system at the sampling instants iT , i ∈ N is given by where F E T x, u : φ E T ; x, u .We note that, since f is typically nonlinear, F E T in 4.3 is not known in most cases, therefore the controller design can be carried out by means of an approximate discrete-time model where h is a modelling parameter, which is typically the step size of the underlying numerical method.The applied numerical scheme approximation has to ensure the closeness of the exact models in the following sense.Let X be a compact set such that X ⊃ X.
Assumption A1.There exists an h * > 0 such that i F A T,h 0, 0 0, F A T,h is continuous in both variables uniformly in h ∈ 0, h * , and Lipschitz continuous w.r.t x in any compact set, uniformly in small h, ii there exists a γ ∈ K such that for all x ∈ X, all u ∈ U, and h ∈ 0, h * .We note that, under the conditions on the function f, then Assumption A1 can be proven for many one-step numerical methods.
The problem is to define a state-feedback controller: using the approximate discrete-time model 4.4 , to practically stabilize the exact discretetime system 4.3 .
Since we want to find a state-feedback controller, it seems to be reasonable to investigate when it does exist.The next assumption formulates, roughly speaking, a necessary condition for the existence of a stabilizing feedback.Assumption A2.There exists an h * > 0 such that the exact discrete-time model 4.3 is practically asymptotically controllable from X to the origin with piecewise constant controllers for all h ∈ 0, h * see e.g., 36, 40 for the definition .
For the solutions of 4.3 and 4.4 with u {u 0 , u 1 , . ..}, and x 0 we will use the notations φ E i x 0 , u and φ A i x 0 , u , respectively.Let N ∈ N be given and let 4.4 be subject to the cost function where u {u 0 , u 1 , . . ., u N−1 }, x A i φ A i x, u , i 0, 1, . . ., N, l h and g are given functions, satisfying the following assumptions.
To ensure the existence and the stabilizing property of the proposed controller, the following assumptions are needed.Assumption A3.Let X 1 X B 1 , i g : X 1 → R is continuous, positive definite radially unbounded, and Lipschitz continuous in any compact set, ii l h x, u is continuous with respect to x and u, uniformly in small h, and Lipschitz continuous in any compact set, iii there exist an h * > 0 and two class-K ∞ functions ϕ 1 and ϕ 2 such that the inequality, The terminal cost function g and/or a terminal constraint set given explicitly or implicitly play a crucial role in establishing the desired stabilizing property.We will assume that g has to be a local control Lyapunov function.
Assumption A4.There exist h * > 0 and η > 0 such that for all x ∈ G η {x : g x ≤ η} there exists a κ x ∈ U such that inequality holds true for all h ∈ 0, h * .
Consider the optimization problem:

4.10
If this optimization problem has a solution denoted by u * x {u * 0 x , u * 1 x , . . ., u * N−1 x }, then the first element of u * is applied at the state x, that is, The value function for the optimal control problem is

4.12
Let h * 0 denote the minimum of the values h * generated by Assumptions A1-A4.

Theorem 4.1 see 36 . If Assumptions A1-A4 hold true, then
i there exist an h * 1 with 0 < h * 1 ≤ h * 0 , and a constant V A max independent of N, such that ii there exist constants N * , L V , and δ V and functions σ 1 , σ 2 ∈ K ∞ such that for all x ∈ X, N > N * , h ∈ 0, h * 1 and i 1, 2, . ..,

4.13
Moreover, for all x, y ∈ X 14 Clearly X ⊂ {x : V N x ≤ V A max }.Theorem 4.1 shows that undersuitable conditions; the state feedback MPC renders the origin to be asymptotically stable for the approximate discrete-time model.These conditions concern directly with the data of the problem and the design parameters the horizon length N, the stage cost l δ , the terminal cost g, and the terminal constraint set G η of the method but not the results of the design procedure.Theorem 4.2 see 36 .Suppose that Assumptions A1-A4 are valid and N is chosen such that N ≥ N * , then, there exists β ∈ KL, and for any δ > 0 there exists an h * > 0 such that for any x 0 ∈ X and h ∈ 0, h * the trajectory of the exact discrete-time system with the MPC, v h x E i 1 Theorem 4.2 lays the foundation for the design of a state feedback MPC via an approximate discrete-time model to achieve practical stability of the exact discrete-time model.Achieving practical stability of the exact model requires that the approximation error can be made sufficiently small.By application of some of one-step numerical approximation formula with possibly variable step size e.g., a Runge-Kutta formula , the approximation error can be made sufficiently small.

Robustness Properties of the MPC
In this section we show the inherent robustness properties of the MPC against to bounded disturbances.We consider the continuous-time model 4.1 with additive disturbances ż t f z t , u t w t , z 0 x 0 , 4.17 where z t ∈ X and w t ∈ W are the state and disturbance, respectively, W is compact, 0 ∈ W. We assume that there exists μ > 0 such that W ⊂ β μ .For a given w : R → R we use the following notation w T i : {w t , t ∈ iT, i 1 T , i ∈ N}.Let us define

4.18
The exact discrete-time model for model 4.17 can be given as: In this paper, the MPC is constructed on the basis of the approximate discrete-time model 4.4 for the nominal model 4.1 .The MPC algorithm consists of performing the following steps at certain instants t i iT : 1 measure the current state of the system z E i ; 2 compute the open-loop optimal control u * to the problem 0 is applied to the system in the interval iT, i 1 T the remaining {u * 1 , . . ., u * N−1 } is discarded ; 4 the procedure is repeated from 1 for the next sampling instants t i 1 i 1 T .
Theorem 4.3 see 37 .Suppose that A1-A4 are valid and N is chosen such as N ≥ N * .Then there exist β ∈ KL, θ ∈ K ∞ , μ * > 0 and for any δ > 0 there exists an h * > 0 such that for any x 0 ∈ X and h ∈ 0, h * , the trajectory of the exact discrete-time model with the MPC, Theorems 4.3 show that under suitable conditions the state feedback MPC practically stabilizes the exact discrete-time model for sufficiently small integration parameter and disturbances.When the full state of the system is not available for feedback, an observer can be designed for estimating the unknown states.In 41 , it is shown that under a set of conditions the output feedback MPC practically stabilize the exact discrete-time model of the plant for sufficiently small approximation and estimation errors.

MPC for the HIV Infection Model
In this section we apply the MPC method proposed in Section 4 to the nominal model 2.1 -2.6 as well as the disturbed model 3.1 -3.6 .We will show that, with a suitable choice of N and functions g and l h , the assumptions of Section 4 can be satisfied.We introduce new variables by the definition ξ ., ξ 6 , then in these new variables the model 3.1 -3.6 takes the form of 4.17 with and w w 1 , w 2 , . . ., w 6 .The disturbance vector w is assumed to be bounded in a compact set containing the origin.Let the compact set X be defined as

5.2
where M 1 , M 2 , and M 3 are as in Proposition 3.1.For the nominal model, the compact set X is defined as X by replacing M i with L i , i 1, 2, 3 where L i are given in Proposition 2.1.
To verify Assumptions A3 and A4, we linearized the nominal system 5.1 i.e., w i 0 around the origin in case of constant controllers, that is, , where u c is given in Proposition 3.1.Let A C be the coefficient matrix of the linearized system.Then the discrete-time model for the linearized system is given by: ξ i 1 e A C T ξ i , i ≥ 0.

5.3
Let the sampling period be chosen to be T 1 and, u 1 u 2 2. The running cost and the terminal cost can be chosen as: where α i are positive weighting constants, P is a positive definite diagonal matrix, and Q is a positive definite symmetric matrix satisfying the Lyapunov equation for the discrete-time system 5.3 : From 5.4 , Assumption A3 is satisfied.Assumption A2 follows from Corollary 2.4 and Assumption A1 holds also true if we choose a suitable numerical integration scheme e.g., the Runge-Kutta formula .To verify Assumption A4, the weights α i and the matrix P have been chosen through a series of numerical experiments as α 1 0.1, α 2 0.05, α 3 0.01, and P diag 0.1, 1, 0.5, 1, 0.1, 1 .It has been verified numerically by solving a constrained  minimization problem with several starting points that Assumption A4 is satisfied over the whole set X. Thus all assumptions of the proposed method can be satisfied with suitable choice of the parameters of the MPC method.

Numerical Results
We perform simulation studies using the parameter values which are listed in Table 1.
All computations are carried out by MATLAB, in particular, the optimal control sequence is computed by the fmincon code of the Optimization Toolbox.To reduce the computational complexity we chose horizon length N to be N 8. Simulations for the continuous-time system are carried out using ode45 program in MATLAB.
In Figures 1, 2, 3, 4, 5, and 6, we show the evolution of the HIV infection model variables for two cases.

Untreated Case
In this case no treatment is used i.e., u 1 u 2 0 .We can see that, for the parameters given in Table 1, R 0 0, 0 6.04 > 1, and . Therefore E 0 is unstable, E 1 exists and stable, while E 2 does not exist.To show the simulation results for this case, we assume that the infection occurs with a certain amount of virus particles v 0.001 and CTL cells z 0.03.Thus the initial conditions for the untreated case are x 0 x 0 , y 0 y 0 , x 1 0 y 1 0 0, v 0 0.001, and z 0 0.03.From Figures 1, 2 can be seen that the concentrations of uninfected CD4 T cells, macrophages, and CTL cells are decaying, while the concentrations of infected CD4 T cells, infected macrophages, and free viruses are increasing.Also we note that the trajectory tends to the stable infected steady state E 1 0.488, 0.5026, 0.1319, 0.0136, 5.1626, 0 .

Treated Case
In this case, the treatment is designed via MPC strategy for the nominal and disturbed HIV infection models.In this case, we assume that the treatment is initiated with 50 time units after the onset of infection.Thus the initial conditions for the MPC are given by x 0 0.4845, x 1 0 0.5024, y 0 0.132, y 1 0 0.0137, v 0 5.1874, and z 0 0.0139.The disturbances are simulated by w i t ∈ η i , i , w i t w i j η i i − η i r j , t ∈ jT, j 1 T , i 1, . . ., 6, j 0, 1, . . ., 6.1 where the parameters r j 's are uniformly distributed random numbers on 0, 1 and η i − i when the system states lie in the interior of the positive orthant R 6 .At the boundary of R 6 , the lower bound η i has to be chosen as the following:  1 and 3 show that when the MPC is applied, the number of uninfected CD4 T cells is increasing as well as the macrophages.This means that the HAART helps the immune system to recover with some fluctuations due to the presence of disturbances.From Figures 2 and 4 we can see that the number of infected CD4 T cells and infected macrophages is decaying during the treatment.Figure 5 shows that after initiation of HAART the viral load drops quickly and it can be kept under a suitable level, with a small controller, corresponding to rather mild dosage of HAART. Figure 6 shows that the CTL cells are decaying for untreated as well as treated cases but with a faster rate than untreated case.The model predictive controller as a function of the time for cases I and II is shown in Figure 7.It is observed that the treatment is initiated with a stronger dosage of HAART and sequentially decreasing over time.Thus we can say that when the MPC strategy is applied in the presence of bounded disturbances the trajectory of the system tends to a ball around the uninfected steady state E 0 and remains there i.e., practical stability .We observe that, for the disturbance-free Case I , the size of the ball is very small due to small numerical errors.For cases II and III , the size of the ball becomes larger and larger by increasing the bounds of the disturbances.

Conclusion
The basic properties of the 6-dimensional model that describes the interaction of the HIV with two target cells, CD4 T cells, and macrophages and takes into account the Cytotoxic T Lymphocytes CTL immune response were studied.The HIV infection model was incorporated to allow some additive disturbances.The HIV infection model incorporates the effect of HAART is considered as a nonlinear control system, where the control input is  defined to be dependent on the drug dose and drug efficiency.The proposed MPC method is applied for determining HAART schedules and stabilizing the HIV infection system around the uninfected steady state.The inherent robustness properties of the MPC were established.

2 Figure 7 :
Figure 7: MPC for cases I and II .
The global stability of E 0 follows from LaSalle's Invariance Principle.
The proof can follow the same as Proposition 2.1.It can be easily to show that M 1 λ 1 λ 2 /σ, M 2 cM 1 /p, M 3 9 is positively invariant.Proof.

Table 1 :
The values of the parameters in the HIV infection model.