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

Authors: A. M. Elaiw; A. M. Shehata
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

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 (MPC-)based method. The MPC is constructed on the basis of an approximate discrete-time model of the HIV infection model. The MPC 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.


---

## Body

## 1. 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–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 HIV-infected 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 Section2, 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.


## 2. 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:
(2.1)x˙=λ1-d1x-β1e-ψ1m1xv,(2.2)x˙1=β1e-ψ1m1xv-ax1-px1z,(2.3)y˙=λ2-d2y-e-ψ1m1β2yv,(2.4)y˙1=e-ψ1m1β2yv-ay1-py1z,(2.5)v˙=e-ψ2m2k(x1+y1)-rv,(2.6)z˙=c(x1+y1)z-bz.Herex(t), x1(t), y(t), y1(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, d1, d2 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 px1z and py1z. 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(x1+y1)z and are decayed at a rate bz. The effect of the RTI and PI drugs is represented by the chemotherapy functions e-ψ1m1(t) and e-ψ2m2(t), where the parameters ψ1 and ψ2 are the efficiencies of RTI and PI drugs, respectively. The variables m1(t) and m2(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 asui(t)=ψimi(t), i=1,2. Then system (2.1)–(2.6) can be considered as a nonlinear control system with (x,x1,y,y1,v,z)′ as the state vector and (u1,u2)′ is the control input vector. We are now ready to present a study on the basic mathematical properties of the model.


### 2.1. Positive Invariance

Now we show that model (2.1)–(2.6) is biologically acceptable in the sense that no population goes negative:
(2.7)x˙|(x=0)=λ1≥0,x˙1|(x1=0)=β1e-u1xv≥0,(x,v≥0),y˙|(y=0)=λ2≥0,y˙1|(y1=0)=β2e-u1yv≥0,(y,v)≥0,v˙|(v=0)=ke-u2(x1+y1)≥0,(x1,y1≥0),z˙|(z=0)=0.This means that the nonnegative orthantℝ+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 numbersL1, L2, and L3 such that the compact set,
(2.8)Γ1={(x,x1,y,y1,v,z)∈R+6:0≤x,x1,y,y1≤L1,0≤z≤L2,0≤v≤L3},
is positively invariant.Proof.
LetX=x+x1+y+y1+(p/c)z, then
(2.9)X˙≤λ1+λ2-σX,
where σ=min{d1,d2,a,b}. Hence 0≤X(t)≤(λ1+λ2)/σ for all t≥0 if X(0)≤(λ1+λ2)/σ. It follows that 0≤x(t),x1(t),y(t),y1(t)≤L1,0≤z(t)≤L2 for all t≥0 if x(0)+x1(0)+y(0)+y1(0)+(p/c)z(0)≤L1, where L1=(λ1+λ2)/σ,L2=cL1/p. On the other hand,
(2.10)v˙(t)≤ke-u2L1-rv≤kL1-rv,
then 0≤v(t)≤L3, for all t≥0 if v(0)≤L3, where L3=kL1/r.


### 2.2. Steady States

We will compute the steady states of system (2.1)–(2.6) under constant controllers, that is, for uj(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
(2.11)R0(u-1,u-2)=e-(u-1+u-2)k(β1x0+β2y0)ar,
where x0=λ1/d1,y0=λ2/d2. We note that R0 can be written as R0=R1+R2, where
(2.12)R1(u-1,u-2)=e-(u-1+u-2)kβ1λ1ard1,R2(u-1,u-2)=e-(u-1+u-2)kβ2λ2ard2
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 I1=cλ1/ab and I2=cλ2/ab (see [12]).Lemma 2.2.
(i) IfR0≤1, then model (2.1)–(2.6) has only one steady state E0.(ii) IfR0>1 and (I1R1/I1+R1)+(I2R2/I2+R2)<1, then there exist two steady states E0 and E1.(iii) IfR0>1 and (I1R1/I1+R1)+(I2R2/I2+R2)>1, then there exist three steady states E0, E1, and E2.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.


### 2.3. Global Stability ofE0

In this section, we prove the global stability of the uninfected steady state by using a Lyapunov approach.Theorem 2.3.
IfR0≤1, then E0 is globally asymptotically stable.Proof.
By the method of [3, 5], we consider a Lyapunov function:
(2.13)W(x,x1,y,y1,v,z)=x0[xx0-ln⁡(xx0)-1]+y0[yy0-ln⁡(yy0)-1]+x1+y1+aeu-2kv+pcz.
We note thatW is defined, continuous, and positive definite for all (x,x1,y,y1,v,z)>0. Also, the global minimum W=0 occurs at the uninfected steady state E0. Further, function W along the trajectories of (2.1)–(2.6) satisfies
(2.14)dWdt=(1-x0x)(λ1-d1x-β1e-u-1xv)+(1-y0y)(λ2-d2y-β2e-u-1yv)+β1e-u-1xv-ax1-px1z+β2e-u-1yv-ay1-py1z+aeu-2k(ke-u-2(x1+y1)-rv)+pc(c(x1+y1)z-bz)=λ1[2-xx0-x0x]+λ2[2-yy0-y0y]+areu-2k[R0-1]v-pλ1aI1z.
SinceI1>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 R0≤1 then dW/dt≤0 for all x,y,v,z>0. The maximal compact invariant set in {(x,x1,y,y1,v,z)∈Γ1:dW/dt=0} is the singleton {E0} when R0≤1. The global stability of E0 follows from LaSalle’s Invariance Principle.Corollary 2.4.
System (2.1)–(2.6) is globally asymptotically controllable to E0 with piecewise constant controllers.Proof.
Letu1(t)=u-1 and u2(t)=u-2 with u-1+u-2>uc, where
(2.15)uc=max⁡{0,ln⁡(kβ1x0+kβ2y0ar)},
then R0(u-1,u-2)<1, therefore the corresponding trajectory will tend to E0 as t→∞.Remark 2.5.
We observe that, ifR0<1, then it is sure that R1<1 and R2<1. But if one considers the three-dimensional model (2.1)-(2.2), and (2.5) and designs a controller such that R1<1, then the whole system may be unstable around E0, because R0>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.


## 2.1. Positive Invariance

Now we show that model (2.1)–(2.6) is biologically acceptable in the sense that no population goes negative:
(2.7)x˙|(x=0)=λ1≥0,x˙1|(x1=0)=β1e-u1xv≥0,(x,v≥0),y˙|(y=0)=λ2≥0,y˙1|(y1=0)=β2e-u1yv≥0,(y,v)≥0,v˙|(v=0)=ke-u2(x1+y1)≥0,(x1,y1≥0),z˙|(z=0)=0.This means that the nonnegative orthantℝ+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 numbersL1, L2, and L3 such that the compact set,
(2.8)Γ1={(x,x1,y,y1,v,z)∈R+6:0≤x,x1,y,y1≤L1,0≤z≤L2,0≤v≤L3},
is positively invariant.Proof.
LetX=x+x1+y+y1+(p/c)z, then
(2.9)X˙≤λ1+λ2-σX,
where σ=min{d1,d2,a,b}. Hence 0≤X(t)≤(λ1+λ2)/σ for all t≥0 if X(0)≤(λ1+λ2)/σ. It follows that 0≤x(t),x1(t),y(t),y1(t)≤L1,0≤z(t)≤L2 for all t≥0 if x(0)+x1(0)+y(0)+y1(0)+(p/c)z(0)≤L1, where L1=(λ1+λ2)/σ,L2=cL1/p. On the other hand,
(2.10)v˙(t)≤ke-u2L1-rv≤kL1-rv,
then 0≤v(t)≤L3, for all t≥0 if v(0)≤L3, where L3=kL1/r.


## 2.2. Steady States

We will compute the steady states of system (2.1)–(2.6) under constant controllers, that is, for uj(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
(2.11)R0(u-1,u-2)=e-(u-1+u-2)k(β1x0+β2y0)ar,
where x0=λ1/d1,y0=λ2/d2. We note that R0 can be written as R0=R1+R2, where
(2.12)R1(u-1,u-2)=e-(u-1+u-2)kβ1λ1ard1,R2(u-1,u-2)=e-(u-1+u-2)kβ2λ2ard2
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 I1=cλ1/ab and I2=cλ2/ab (see [12]).Lemma 2.2.
(i) IfR0≤1, then model (2.1)–(2.6) has only one steady state E0.(ii) IfR0>1 and (I1R1/I1+R1)+(I2R2/I2+R2)<1, then there exist two steady states E0 and E1.(iii) IfR0>1 and (I1R1/I1+R1)+(I2R2/I2+R2)>1, then there exist three steady states E0, E1, and E2.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.


## 2.3. Global Stability ofE0

In this section, we prove the global stability of the uninfected steady state by using a Lyapunov approach.Theorem 2.3.
IfR0≤1, then E0 is globally asymptotically stable.Proof.
By the method of [3, 5], we consider a Lyapunov function:
(2.13)W(x,x1,y,y1,v,z)=x0[xx0-ln⁡(xx0)-1]+y0[yy0-ln⁡(yy0)-1]+x1+y1+aeu-2kv+pcz.
We note thatW is defined, continuous, and positive definite for all (x,x1,y,y1,v,z)>0. Also, the global minimum W=0 occurs at the uninfected steady state E0. Further, function W along the trajectories of (2.1)–(2.6) satisfies
(2.14)dWdt=(1-x0x)(λ1-d1x-β1e-u-1xv)+(1-y0y)(λ2-d2y-β2e-u-1yv)+β1e-u-1xv-ax1-px1z+β2e-u-1yv-ay1-py1z+aeu-2k(ke-u-2(x1+y1)-rv)+pc(c(x1+y1)z-bz)=λ1[2-xx0-x0x]+λ2[2-yy0-y0y]+areu-2k[R0-1]v-pλ1aI1z.
SinceI1>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 R0≤1 then dW/dt≤0 for all x,y,v,z>0. The maximal compact invariant set in {(x,x1,y,y1,v,z)∈Γ1:dW/dt=0} is the singleton {E0} when R0≤1. The global stability of E0 follows from LaSalle’s Invariance Principle.Corollary 2.4.
System (2.1)–(2.6) is globally asymptotically controllable to E0 with piecewise constant controllers.Proof.
Letu1(t)=u-1 and u2(t)=u-2 with u-1+u-2>uc, where
(2.15)uc=max⁡{0,ln⁡(kβ1x0+kβ2y0ar)},
then R0(u-1,u-2)<1, therefore the corresponding trajectory will tend to E0 as t→∞.Remark 2.5.
We observe that, ifR0<1, then it is sure that R1<1 and R2<1. But if one considers the three-dimensional model (2.1)-(2.2), and (2.5) and designs a controller such that R1<1, then the whole system may be unstable around E0, because R0>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.


## 3. 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:
(3.1)x˙=λ1-d1x-β1e-u1xv+w1,(3.2)x˙1=β1e-u1xv-ax1-px1z+w2,(3.3)y˙=λ2-d2y-β2e-u1yv+w3,(3.4)y˙1=β2e-u1yv-ay1-py1z+w4,(3.5)v˙=ke-u2(x1+y1)-rv+w5,(3.6)z˙=c(x1+y1)z-bz+w6.In (3.1)–(3.6), wi(t) describe model uncertainties/disturbances which are assumed to satisfy the following bound:
(3.7)‖wi(t)‖≤ϵi,ϵi≥0,i=1,…,6.It is important to show that in the presence of the disturbances, the nonnegative orthantℝ+6 is positively invariant and the solutions of model (3.1)–(3.6) are bounded.


### 3.1. Positive Invariance

For the model (3.1)–(3.6), the following conditions guarantee that the nonnegative orthant ℝ+6 is positively invariant (see [20]):
(3.8)x˙|(x=0)=λ1+w1≥0,ifw1≥-λ1,x˙1|(x1=0)=β1e-u-1xv+w2≥0,ifw2≥-β1e-u-1xv,y˙|(y=0)=λ2+w3≥0,ifw3≥-λ2,y˙1|(y1=0)=β2e-u-1yv+w4≥0,ifw4≥-β2e-u-1yv,v˙|(v=0)=ke-u-2(x1+y1)+w5≥0,ifw5≥-ke-u-2(x1+y1),z˙|(z=0)=w6,ifw6≥0.We note that the conditions in (3.8) give a lower bound of the disturbances only at the the boundary of ℝ+6.Proposition 3.1.
Suppose that the disturbances satisfy the bound (3.7) then there exist positive numbers M1,M2, and M3 such that the compact set,
(3.9)Γ2={(x,x1,y,y1,v,z)∈R+6:0≤x,x1,y,y1≤M1,0≤z≤M2,0≤v≤M3},
is positively invariant.Proof.
The proof can follow the same as Proposition2.1. It can be easily to show that M1=(λ1+λ2+ϵ)/σ, M2=cM1/p, M3=(kM1+ϵ5)/r where ϵ=(ϵ1+ϵ2+ϵ3+ϵ4+(p/c)ϵ6).


## 3.1. Positive Invariance

For the model (3.1)–(3.6), the following conditions guarantee that the nonnegative orthant ℝ+6 is positively invariant (see [20]):
(3.8)x˙|(x=0)=λ1+w1≥0,ifw1≥-λ1,x˙1|(x1=0)=β1e-u-1xv+w2≥0,ifw2≥-β1e-u-1xv,y˙|(y=0)=λ2+w3≥0,ifw3≥-λ2,y˙1|(y1=0)=β2e-u-1yv+w4≥0,ifw4≥-β2e-u-1yv,v˙|(v=0)=ke-u-2(x1+y1)+w5≥0,ifw5≥-ke-u-2(x1+y1),z˙|(z=0)=w6,ifw6≥0.We note that the conditions in (3.8) give a lower bound of the disturbances only at the the boundary of ℝ+6.Proposition 3.1.
Suppose that the disturbances satisfy the bound (3.7) then there exist positive numbers M1,M2, and M3 such that the compact set,
(3.9)Γ2={(x,x1,y,y1,v,z)∈R+6:0≤x,x1,y,y1≤M1,0≤z≤M2,0≤v≤M3},
is positively invariant.Proof.
The proof can follow the same as Proposition2.1. It can be easily to show that M1=(λ1+λ2+ϵ)/σ, M2=cM1/p, M3=(kM1+ϵ5)/r where ϵ=(ϵ1+ϵ2+ϵ3+ϵ4+(p/c)ϵ6).


## 4. 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ℝ and ℕ. The notation ℝ+ denote the set of real numbers in the interval [0,∞). A continuous function σ:ℝ+→ℝ+ is of class-𝒦 if σ(0)=0, σ(s)>0 for all s>0 and it is strictly increasing. It is of class-𝒦∞ if it is of class-𝒦 and σ(s)→∞ when s→∞. A continuous function β:ℝ+×ℝ+→ℝ+ is of class-𝒦ℒ if β(s,τ) is of class-𝒦 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 ℬρ to be a ball of radius ρ centered at the origin. We introduce the following notation Yρ=Y∩ℬρ for any set Y.Consider a continuous-time nonlinear control system given by(4.1)x˙(t)=f(x(t),u(t)),x(0)=x0,
where x(t)∈ℝn, u(t)∈U⊂ℝm are the state and control input, respectively, f:ℝn×U→ℝ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(4.2)u(t)=u(iT)=:ui,fort∈[iT,(i+1)T),i∈N,
where T>0 is the control sampling period which is fixed.We will assume that there is a compact set𝒳¯⊂ℝ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 discrete-time model which describes the behavior of the system at the sampling instants iT, i∈ℕ is given by
(4.3)xi+1E=FTE(xiE,ui),x0E=x0,
where FTE(x,u)∶=ϕE(T;x,u). We note that, since f is typically nonlinear, FTE 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
(4.4)xi+1A=FT,hA(xiA,ui),x0A=x0,
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 𝒳 be a compact set such that 𝒳⊃𝒳¯.Assumption A1.
There exists anh*>0 such that(i)
FT,hA(0,0)=0, FT,hA 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γ∈𝒦 such that
(4.5)‖FTE(x,u)-FT,hA(x,u)‖≤Tγ(h),for all x∈𝒳, all u∈U, and h∈(0,h*].
We note that, under the conditions on the functionf, then Assumption A1 can be proven for many one-step numerical methods.The problem is to define a state-feedback controller:(4.6)vh:X⟶U
using the approximate discrete-time model (4.4), to practically stabilize the exact discrete-time 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 anh*>0 such that the exact discrete-time model (4.3) is practically asymptotically controllable from 𝒳 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={u0,u1,…}, and x0 we will use the notations ϕiE(x0,u) and ϕiA(x0,u), respectively.
LetN∈ℕ be given and let (4.4) be subject to the cost function
(4.7)JT,h(N,x,u)=∑i=0N-1‍Tlh(xiA,ui)+g(xNA),
where u={u0,u1,…,uN-1}, xiA=ϕiA(x,u), i=0,1,…,N, lh 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𝒳1=𝒳+ℬ1, (i) g:𝒳1→ℝ is continuous, positive definite radially unbounded, and Lipschitz continuous in any compact set,(ii)lh(x,u) is continuous with respect to x and u, uniformly in small h, and Lipschitz continuous in any compact set,(iii) there exist anh*>0 and two class-𝒦∞ functions φ1 and φ2 such that the inequality,
(4.8)φ1(‖x‖)≤lh(x,u)≤φ2(‖x‖)+φ2(‖u‖),holds for all x∈𝒳1,u∈U and h∈(0,h*].
The terminal cost functiong 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 existh*>0 and η>0 such that for all x∈𝒢η={x:g(x)≤η} there exists a κ(x)∈U such that inequality
(4.9)g(FT,hA(x,κ(x)))-g(x)≤-Tlh(x,κ(x))
holds true for all h∈(0,h*].Consider the optimization problem:(4.10)PT,hA(N,x):min⁡{JT,δ(N,x,u):ui∈U}.If this optimization problem has a solution denoted byu*(x)={u0*(x),u1*(x),…,uN-1*(x)}, then the first element of u* is applied at the state x, that is,
(4.11)vh(x)=u0*(x).The value function for the optimal control problem is(4.12)VN(x)=JT,h(N,x,u*(x)).Leth0* 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 anh1* with 0<h1*≤h0*, and a constant VmaxA independent of N, such that VN(x)≤VmaxA   for all x∈𝒳, h∈(0,h1*] and N∈ℕ,(ii)
there exist constantsN*,LV, and δV and functions σ1,σ2∈𝒦∞ such that for all x∈𝒳, N>N*, h∈(0,h1*] and i=1,2,…,
(4.13)σ1(‖x‖)≤VN(x)≤σ2(‖x‖),VN(FT,hA(x,vh(x)))-VN(x)≤-Tφ1(‖x‖).
Moreover, for allx,y∈𝒳1 with ∥x-y∥≤δV(4.14)|VN(x)-VN(y)|≤LV‖x-y‖
for all h∈(0,h1*].Clearly𝒳⊂{x:VN(x)≤VmaxA}.Theorem4.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 𝒢η) of the method but not the results of the design procedure.Theorem 4.2 (see [36]).
Suppose that Assumptions A1–A4 are valid andN is chosen such that N≥N*, then, there exists β∈𝒦ℒ, and for any δ>0 there exists an h*>0 such that for any x0∈𝒳 and h∈(0,h*] the trajectory of the exact discrete-time system with the MPC, vh(4.15)xi+1E=FTE(xiE,vh(xiE)),x0E=x0,
satisfies
(4.16)‖xiE‖≤β(‖x0‖,iT)+δ,foralli≥0.Theorem4.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.


### 4.1. 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
(4.17)z˙(t)=f(z(t),u(t))+w(t),z(0)=x0,
where z(t)∈𝒳 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:ℝ+→ℝ we use the following notation wT[i]∶={w(t),t∈[iT,(i+1)T),i∈N}. Let us define
(4.18)Wμ={w∈L[0,∞)∞:w(t)∈W,a.e.t∈[0,∞)with‖w‖∞≤μ},Wμ={wT[i],w∈Wμ,i=0,1,2,…}.The exact discrete-time model for model (4.17) can be given as:
(4.19)zi+1E=FTE(ziE,ui,wT[i]),z0E=x0.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 ti=iT:(1)
measure the current state of the systemziE;(2)
compute the open-loop optimal controlu* to the problem PT,hA(N,ziE);(3)
the controlvh(ziE)∶=u0* is applied to the system in the interval [iT,(i+1)T) (the remaining {u1*,…,uN-1*} is discarded);(4)
the procedure is repeated from (1) for the next sampling instantsti+1=(i+1)T.Theorem 4.3 (see [37]).
Suppose that A1–A4 are valid andN is chosen such as N≥N*. Then there exist β∈𝒦ℒ,θ∈𝒦∞,μ*>0 and for any δ>0 there exists an h*>0 such that for any x0∈𝒳 and h∈(0,h*], the trajectory of the exact discrete-time model with the MPC, vh(4.20)zi+1E=FTE(ziE,vh(ziE),wT[i]),z0E=x0,
and wT[i]∈Wμ*satisfies
(4.21)‖ziE‖≤β(‖x0‖,iT)+θ(μ*)+δ,foralli≥0.Theorems4.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.


## 4.1. 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
(4.17)z˙(t)=f(z(t),u(t))+w(t),z(0)=x0,
where z(t)∈𝒳 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:ℝ+→ℝ we use the following notation wT[i]∶={w(t),t∈[iT,(i+1)T),i∈N}. Let us define
(4.18)Wμ={w∈L[0,∞)∞:w(t)∈W,a.e.t∈[0,∞)with‖w‖∞≤μ},Wμ={wT[i],w∈Wμ,i=0,1,2,…}.The exact discrete-time model for model (4.17) can be given as:
(4.19)zi+1E=FTE(ziE,ui,wT[i]),z0E=x0.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 ti=iT:(1)
measure the current state of the systemziE;(2)
compute the open-loop optimal controlu* to the problem PT,hA(N,ziE);(3)
the controlvh(ziE)∶=u0* is applied to the system in the interval [iT,(i+1)T) (the remaining {u1*,…,uN-1*} is discarded);(4)
the procedure is repeated from (1) for the next sampling instantsti+1=(i+1)T.Theorem 4.3 (see [37]).
Suppose that A1–A4 are valid andN is chosen such as N≥N*. Then there exist β∈𝒦ℒ,θ∈𝒦∞,μ*>0 and for any δ>0 there exists an h*>0 such that for any x0∈𝒳 and h∈(0,h*], the trajectory of the exact discrete-time model with the MPC, vh(4.20)zi+1E=FTE(ziE,vh(ziE),wT[i]),z0E=x0,
and wT[i]∈Wμ*satisfies
(4.21)‖ziE‖≤β(‖x0‖,iT)+θ(μ*)+δ,foralli≥0.Theorems4.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.


## 5. MPC for the HIV Infection Model

In this section we apply the MPC method proposed in Section4 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 lh, the assumptions of Section 4 can be satisfied. We introduce new variables by the definition ξ1=x-x0, ξ2=x1, ξ3=y-y0, ξ4=y1, ξ5=v, ξ6=z. Let ξ=(ξ1,ξ2,..,ξ6)′, then in these new variables the model (3.1)–(3.6) takes the form of (4.17) with
(5.1)f(ξ,u)=(λ1-d1(ξ1+x0)-β1e-u1(ξ1+x0)ξ5β1e-u1(ξ1+x0)ξ5-aξ2-pξ2ξ6λ2-d2(ξ3+y0)-β2e-u1(ξ3+y0)ξ5β2e-u1(ξ3+y0)ξ5-aξ4-pξ4ξ6ke-u2(ξ2+ξ4)-rξ5c(ξ2+ξ4)ξ6-bξ6),
and w=(w1,w2,…,w6)′. The disturbance vector w is assumed to be bounded in a compact set containing the origin. Let the compact set 𝒳 be defined as
(5.2)X={ξ∈R6:-x0≤ξ1≤M1-x0,-y0≤ξ3≤M1-y0,0≤ξ2,ξ4≤M1,0≤ξ5≤M3,0≤ξ6≤M2},
where M1, M2, and M3 are as in Proposition 3.1. For the nominal model, the compact set 𝒳¯ is defined as 𝒳 by replacing Mi with Li, i=1,2,3 where Li are given in Proposition 2.1.To verify Assumptions A3 and A4, we linearized the nominal system (5.1) (i.e., wi=0) around the origin in case of constant controllers, that is, u1(t)=u¯1>u1c,u2(t)=u¯2>u2c with u1c+u2c=uc, where uc is given in Proposition 3.1. Let AC be the coefficient matrix of the linearized system. Then the discrete-time model for the linearized system is given by:
(5.3)ξi+1=eACTξi,i≥0.Let the sampling period be chosen to beT=1 and, u¯1=u¯2=2. The running cost and the terminal cost can be chosen as:
(5.4)lh(ξ,u)=α1ξ′Qξ+α2(u1-u1c)2+α3(u2-u2c)2,g(ξ)=ξ′Pξ,
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):
(5.5)Q=-(AT′PAT-P),AT=eACT.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 𝒳¯. Thus all assumptions of the proposed method can be satisfied with suitable choice of the parameters of the MPC method.


## 6. Numerical Results

We perform simulation studies using the parameter values which are listed in Table1.Table 1
The values of the parameters in the HIV infection model.
Parameters
λ
1
d
1
β
1
λ
2
d
2
β
2
a
k
r
c
p
b
Values
0.3
0.1
0.1
0.02
0.1
0.01
0.5
1.0
0.1
0.01
0.03
0.02All 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 lengthN to be N=8. Simulations for the continuous-time system are carried out using ode45 program in MATLAB.In Figures1, 2, 3, 4, 5, and 6, we show the evolution of the HIV infection model variables for two cases.Figure 1
The evolution of uninfected CD4+ T cells.Figure 2
The evolution of infected CD4+ T cells.Figure 3
The evolution of uninfected macrophages.Figure 4
The evolution of infected macrophages.Figure 5
The evolution of free viruses.Figure 6
The evolution of CTL cells.Untreated Case
In this case no treatment is used (i.e.,u1=u2=0). We can see that, for the parameters given in Table 1, R0(0,0)=6.04>1, and (I1R1/(I1+R1))+(I2R2/(I2+R2))=0.299<1. Therefore E0 is unstable, E1 exists and stable, while E2does 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)=x0, y(0)=y0,x1(0)=y1(0)=0, v(0)=0.001, and z(0)=0.03. From Figures 1, 2, 3, 4, 5, and 6, it 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 E1=(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 with50 time units after the onset of infection. Thus the initial conditions for the MPC are given by x(0)=0.4845, x1(0)=0.5024, y(0)=0.132, y1(0)=0.0137, v(0)=5.1874, and z(0)=0.0139. The disturbances are simulated by wi(t)∈[ηi,ϵi],
(6.1)wi(t)=wi(j)=ηi+(ϵi-ηi)r(j),t∈[jT,(j+1)T),i=1,…,6,j=0,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 ℝ+6. At the boundary of ℝ+6, the lower bound ηi has to be chosen as the following:
(6.2)η1=max⁡{-λ1,-ϵ1},η2=max⁡{-β1e-u1(ξ1+x0)ξ5,-ϵ2},η3=max⁡{-λ2,-ϵ3},η4=max⁡{-β2e-u1(ξ3+y0)ξ5,-ϵ4},η5=max⁡{-ke-u2(ξ2+ξ4),-ϵ5},η6=0,
to guarantee that the positive orthant ℝ+6 is positively invariant.
Figures1, 2, 3, 4, 5, and 6 show also the application of the MPC to the HIV infection models in following cases.
Case (I): nominal HIV infection model that is, wi(t)=0.
Case (II): disturbed HIV infection model with ϵ1=0.03, ϵ2=0.005, ϵ3=0.0001, ϵ4=0.0005, ϵ5=0.04, ϵ6=0.0001.
Case (III): disturbed HIV infection model with ϵ1=0.1, ϵ2=0.01, ϵ3=0.005, ϵ4=0.001, ϵ5=0.08, ϵ6=0.0002.Figures1 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 E0 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.Figure 7
MPC for cases (I) and (II).


## 7. 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.



---

*Source: 963864-2012-08-15.xml*