Eradication Conditions of Infected Cell Populations in the 7-Order HIV Model with Viral Mutations and Related Results

: In this paper, we study possibilities of eradication of populations at an early stage of a patient’s infection in the framework of the seven-order Stengel model with 11 model parameters and four treatment parameters describing the interactions of wild-type and mutant HIV particles with various immune cells. We compute ultimate upper bounds for all model variables that deﬁne a polytope containing the attracting set. The theoretical possibility of eradicating HIV-infected populations has been investigated in the case of a therapy aimed only at eliminating wild-type HIV particles. Eradication conditions are expressed via algebraic inequalities imposed on parameters. Under these conditions, the concentrations of wild-type HIV particles, mutant HIV particles, and infected cells asymptotically tend to zero with increasing time. Our study covers the scope of acceptable therapies with constant concentrations and values of model parameters where eradication of infected particles/cells populations is observed. Sets of parameter values for which Stengel performed his research do not satisfy our local asymptotic stability conditions. Therefore, our exploration develops the Stengel results where he investigated using the optimal control theory and numerical dynamics of his model and came to a negative health prognosis for a patient. The biological interpretation of these results is that after a sufﬁciently long time, the concentrations of wild-type and mutant HIV particles, as well as infected cells will be maintained at a sufﬁciently low level, which means that the viral load and the concentration of infected cells will be minimized. Thus, our study theoretically conﬁrms the possibility of efﬁcient treatment beginning at the earliest stage of infection. Our approach is based on a combination of the localization method of compact invariant sets and the LaSalle theorem.


Introduction
Human immunodeficiency virus (HIV) infection as the cause of AIDS has attracted the attention of many researchers from various fields around the world since the 1980s. In particular, the interest of many scientists has been focused on the elaboration and studies of mathematical models which describe the immunological response to infection with HIV. There are different types of such dynamical models that characterize interactions of HIV with CD4-expressing cells including helper T cells, macrophages, and natural killer cells. The basic studies in this area are contained in seminal works of [1][2][3][4][5]. The mentioned researches were continued in papers [6][7][8][9][10][11][12][13][14] where various dynamical issues related to HIV models are explored, such as positive invariance properties, boundedness of positive half trajectories of this model, stability analysis of equilibrium points, the existence of an orbitally asymptotically stable periodic solution, and bifurcations. These articles did not address the propensity of viral mutations to replicate HIV. Taking HIV mutations into account gives a more realistic picture of the human infectious process. This leads to more complete and complex models. However, as far as the authors know, HIV models with mutant infection have not been sufficiently studied yet [15][16][17][18][19][20][21][22].
In [16], Stengel constructed the seven-order model of HIV mutant infection and discussed the feasibility and effectiveness of the optimal therapies for his model in which virtual therapies for wild-type infections are incorporated: x 1 x 2 − (a 6 + a 9 )x 3 ; x 4 = a 9 x 3 − a 4 x 4 ; x 5 = a 3 a 4 a 10 x 4 − a 1 x 5 − a 2 a 11 x 2 x 5 + a 3 a 4 x 7 ; x 6 = a 2 a 11 x 2 x 5 − (a 6 + a 9 )x 6 ; where the notation is utilized. The model (5) has been obtained from 4D equations in [3] by adding equations describing the dynamics of the mutant HIV strain. The four state variables of the initial model of [3] represent concentrations of free, wild-type HIV particles (x 1 ), uninfected Th cells (x 2 ), latently infected Th cells (x 3 ), and productively infected Th cells (x 4 ) in both the periphery and lymphoid organs. Other state variables of (1) represent concentrations of the mutant HIV strain (x 5 ), proviral Th cells infected by the mutant strain (x 6 ), and Th cells productively infected by the new strain (x 7 ). Parameters have the following biological meanings: a 1 is the death rate of free virions; a 2 is the rate at which CD4 + T cells become infected by free virions; a 3 is the number of free virions produced by x 4 cells; a 4 is the death rate of the actively infected CD4 + T cell population; a 5 is the source term for uninfected CD4 + T cells; a 6 is the death rate of the uninfected CD4 + T cell population; a 7 is the growth rate for the CD4 + T cell population; a 8 is the maximum CD4 + T cell population level; a 9 is the rate at which x 3 cells convert to actively infected cells; a 10 is the mutation rate; a 11 is the fitness of the mutant strain.
In [23], it is indicated that the main damage to the immune system occurs in the first weeks of infection, when the diversity of virions is low. The model in [16] describes the interaction of the HIV-immune system only at an early stage of infection, when wild-type virions and virions of the first mutant strain attack the patient, but new mutant strains have not yet appeared.
Treatment or control parameters u i are supposed to be constant, u i ∈ [0, 1), and are defined as follows: they are concentrations of protease inhibitor (u 1 ), fusion inhibitor (u 2 ), Th cell enhancer (u 3 ), and reverse transcription inhibitor (u 4 ). According to Equation (1), applied therapy can affect wild HIV particles, but do not possess direct effects on the mutant HIV strain. Stengel raised the question of whether the complex interactions described in the upper three equations of (1), as well as the presence of the mutant HIV strain variable in the second equation, may have a "good effect" on the patient's health. Based on clinical practice, according to which "HIV infection is never cured, and requires continuous treatment to maintain a state in remission", Stengel studied the possibility of optimal therapy for treating HIV infection at certain parameter values. Mathematically, his research is based on the steepest descent algorithm and Pontryagin's maximum principle.
The purpose of our work is to investigate some qualitative features of the model (1). Our research provides a positive answer to Stengel's question for certain ranges of parameter values within the broad framework of the rigorous dynamic analysis (1) carried out in this article. With this goal, we find equilibrium points and provide local asymptotic stability (LAS) conditions for the infection-free equilibrium point, prove the existence of the attracting set, and calculate ultimate upper bounds for the polytope containing the attracting set. Further, we show that the dynamics of the model (1) theoretically makes it possible to eradicate the infection by an appropriate choice of treatment parameters if the model parameters satisfy a number of algebraic inequalities. Namely, we find the curious dynamic property of (1), which is that the LAS conditions of the infection-free equilibrium point imply its global asymptotic stability (GAS) conditions, that is, these GAS conditions cannot be improved. Another interesting issue found here is that these conditions do not depend on controls u 3 and u 4 . Moreover, we describe the case when these conditions do not depend on rest controls u 1 and u 2 as well.
The ranges of parameters at which the cure is achieved in the studied model can be considered as target ranges for real biomedical problems, as well as when choosing parameter control. The biological feasibility issues of numerical values are not the focus of this study.
Biologically, the global asymptotic eradication of infection (GAS) property means that after a sufficiently long observation period, concentrations of wild and mutant HIV particles and infected cells are maintained at a fairly low level. Thus, the results of our study may be applicable in subsequent studies in the case of early initiation of therapy for the patient, when the time period after infection is short.
Our research is based on the localization method of compact invariant sets (LMCIS) [24], and the LaSalle theorem. It should be noted that earlier, the LMCIS was effectively used in the study of many models taken from chaos theory, see, for example [25], cosmology [26], mathematical oncology [27,28], mathematical inclusions [29], and others.
The structure of this paper is the following. In Section 2 we describe the LMCIS which is utilized in combination with the LaSalle theorem for obtaining conditions for the locations of ω-limit sets in coordinate planes. Section 3 contains preliminary remarks. In Section 4, formulas for equilibrium points and LAS conditions for the infection-free equilibrium point E 0 are provided. In Section 5 we derive ultimate upper bounds for all state variables; these bounds define the localization polytope containing the attracting set. Section 6 contains main results of this paper: we present the GAS conditions, and, besides, we concern two other issues of ultimate dynamics of (1). Finally, concluding remarks are given in Section 7.

On the Localization Problem of Compact Invariant Sets
For the reader's convenience, we give reminders on a few helpful notions. We consider a nonlinear systemẋ This function is used in the solution of the localization problem of compact invariant sets, and is called a localizing function. Suppose that the function h is not the first integral of the system (2). By h| U , we denote the restriction of h on a set U ⊂ R n .
Assertion 1 (see [24]). For any h(x) ∈ C 1 (R n ), all compact invariant sets of the system (2) located in U are contained in the set

Preliminary Remarks
The system (1) is defined in which is the positively invariant domain.
In order to simplify computations, we introduce the following notations: Utilizing these notations, we come to the system We note that (4) is the extension of the four-dimensional system from [3], in which the resistant virus modification is included. The latter system has the form or after using notations introduced above, Below by f , we denote the vector field corresponding to (4). Next, we point to several properties of system dynamics of (4).
(2) This system possesses the invariant plane x 1 = x 3 = x 4 = 0; this is the case when there are no wild-type HIV particles and no Th cells infected by them. The subsystem defined on this plane is explored in Section 7. (3) Suppose that a 5 is a source term for uninfected CD4 + T cells that is zero, the concentration of uninfected Th cells (x 2 ) is zero as well, and between the death rate and growth rate of the uninfected CD4 + T cell population, the following inequalities are fulfilled: a 6 > a 7 and u 3 < a 6 a −1 7 − 1. Then the system (4) becomes a linear asymptotically stable system for which free, wild-type HIV particles and all other cell populations vanish after a sufficiently long observation time.

Equilibrium Points
Firstly, in order to find equilibrium points, we eliminate variables x 3 , x 4 , x 6 , x 7 using equations f i (x) = 0, i = 3, 4, 6, 7. As a result, we get the system: Then we come to two cases: (1) x 1 = 0; (2) x 1 = 0. In the first case, we obtain the system It follows from the second equation that either and we obtain the equilibrium point and we get the equation respecting x 5 : Its maximal root is given by We notice that x 50 ≥ 0, provided a 5 ≥ d 5 .
In that way, we come to the system of equations x 1 x 22 = 0 with respect to x 1 , x 5 . Eliminating x 5 we obtain the first equation in the form We come to the equation This equation has the unique positive root of the form This solution gives us the equilibrium point , r 3 x 10 , c 1 r 3 x 10 a 6 + a 9 , c 1 a 9 r 3 x 10 a 4 (a 6 + a 9 ) .
The equilibrium point E 2 ∈ R 7 +,0 if conditions are fulfilled.

Example 1.
Let us select the following set of parameters: Then the system (6)  Now we find the stability condition of E 0 . The Jacobian matrix taken at E 0 has the following form: where by * we denote non-essential elements of J 0 .
The spectrum of the matrix J 0 contains one evident eigenvalue a 8 x 20 < 0 and eigenvalues of the matrixĴ 0 which has the block-triangular form: Therefore, the spectrum ofĴ 0 is the union of spectra of matrices Both matrices A 1 and A 2 are matrices of the following special type where all α i and β i are positive. The characteristic polynomial (with an opposite sign) of the matrix A has the form One can show (for instance, using the Hurwitz criterion) that roots of ξ A (t) are negative if, and only if α 1 α 2 α 3 > β 1 β 2 β 3 .
Hence, the condition (8) is the stability condition of the matrix A.
Coming back to matrices A 1 , A 2 we write the stability condition of E 0 in the form: a 4 (a 6 + a 9 )(a 1 + b 1 x 20 ) > b 2 b 4 a 9 x 20 , a 4 (a 6 + a 9 )(a 1 + c 1 x 20 ) > c 1 c 3 a 9 x 20 , or in notations (7) we obtain the condition In particular, E 0 is LAS with d 1 < 1, d 2 < 1, which means that E 1 , E 2 are contained in the half-space x 2 < 0.

Ultimate Upper Bounds
In this section, we derive upper bounds for the dynamics of the system (4) in nonnegative orthant. Upper bounds give us ultimate maximal values for all cell populations involved into the model. We construct a polytope containing all compact invariant sets in R 7 +,0 . This polytope is a positively invariant set and, therefore, limits biologically significant bounded dynamics. Moreover, the polytope turns out to be globally attractive, so all semi-trajectories are bounded as t → +∞, and their ω-limit sets are contained in the polytope. Now we demonstrate that this polytope can be constructed by using linear localizing functions.
1. Firstly, we apply the function h 1 (x) = x 2 . Then (we recall that L f h is the Lie derivative of h with respect to f ) and the set S(h 1 ) is defined by The last formula can be rewritten as where The Equation (10) is quadratic, respecting x 2 , and its larger root increases while q 0 decreases and the left summand increases. Thus, The value of h 1,inf is reached at the maximum values of the variables x 1 , x 3 , x 4 , x 5 , x 6 , x 7 . We see that h 1,inf = 0. Thus, we get the localizing set {x ∈ R n | 0 ≤ x 2 ≤ x 2,max }.
2. Let us take the next function h 2 (x) = x 2 + x 3 . Then Applying the change x 3 = h 2 − x 2 in the equation L f h 2 = 0, describing the set S(h 1 ), we get that Solving this linear equation respecting h 2 , we obtain that with It follows from (12) that h 2,inf = 0 and The fractional linear function (13) is monotonic within (0, ∞) and we derive that This estimate can be improved if we use the value x 2 max in (13) instead of x 2 → +∞. Thus, we get the localization set {x ∈ R n | x 2 + x 3 ≤ h 2,sup }. Therefore, we come to the estimate 3. Let us apply the localization function h 3 (x) = x 4 . In this case, the set S(h 3 ) is given by Hence, we obtain the localization set defined by 4. Let us employ the localization function h 4 (x) = x 1 . The set S(h 4 ) is defined by Thus, we obtain the localization set defined by 5. Next, let us utilize the localization function h 5 (x) = x 2 + x 6 . Then The set S(h 5 ) is defined by Consequently, we have on S(h 5 ) that As a result, we derive that Therefore, we have the localization set {x ∈ R n | x 2 + x 6 ≤ h 5,sup } which provides the bound for x 6 : 6. Now we take the localization function h 6 (x) = x 7 . Then the set S(h 6 ) is given by Taking into account the bound for x 6 , we get that 7. Now we take the localization function h 7 (x) = x 5 . The set S(h 7 ) is defined by the equation Using bounds x 4,max ; x 7,max , we get the localizing set x 7,max = a 9 c 2 a 4 a 1 h 2,sup + a 9 c 3 a 4 a 1 h 5,sup .

Remark 1.
The formula x 2 max = x 20 provides the best upper bound for the concentration of uninfected Th cells. Indeed, it was established that x 20 = h 1,sup and, at the same time, x 20 is the coordinate of an equilibrium point. We see that the bound h 1,sup is not refinable.
To summarize all these results, we arrive at: Theorem 1. All compact invariant sets of the system (4) located in R 7 +,0 are contained in the polytope as well. Here, x 4,max = a 9 a 4 x 3,max , x 5,max = c 2 a 9 a 4 a 1 x 3,max + c 3 a 9 a 4 a 1 x 6,max , x 7,max = a 9 a 4 x 6,max .

Remark 2.
Based on this theorem, we obtain the lower bound for x 2 : x 3,max + x 4,max + x 6,max + x 7,max .

On the Location of ω-Limit Sets
It follows from Theorem 1 that all ω-limit sets are located in Π. Generally speaking, this polytope is not a positive invariant domain. Let us take a smaller polytope Π 0 defined by We recall that this polytope is the localization set obtained as a result of applying the iterative procedure using localizing functions h 1 , . . . , h 7 . Herewith, we note that L f h i (x) < 0, with x taken from the domain h i (x) > h i,sup and other inequalities are satisfied. Indeed, L f h i (x) in this domain keeps its sign, and one can verify that for large values of h i (x), the Lie derivative L f h i (x) is negative. This conclusion means that the vector field f is directed inward Π 0 on the boundary of Π 0 , and this polytope is positively invariant.
Actually, one can prove a stronger assertion. Arguing as above, we conclude that for any C > 0, the polytope Π 0 (C) defined by inequalities is positively invariant. It follows from this fact that a trajectory exiting a point M 0 ∈ R 7 +,0 , remains in the polytope Π 0 (C) for sufficiently large C. However, Π 0 (C) is a compact set. Therefore, the ω-limit set of the trajectory is a compact set belonging to Π 0 . We conclude that Π 0 is a globally attracting set. Theorem 2. If conditions (9) hold for the system (4), then the equilibrium point E 0 attracts all trajectories in R 7 +,0 .
Proof. Let us take the function with positive parameters η i , i = 1, . . . , 5. Then we compute that The condition L f h 8 (x) ≤ 0 holds if the system of inequalities , η 2 (a 6 + a 9 ) > a 9 , holds in this polytope, that is, under the condition 0 ≤ x 2 ≤ x 2,max = x 20 . Excluding parameters η 2 , η 5 , we come to the system Next, we exclude η 4 : a 9 c 3 a 4 (a 6 + a 9 ) Finally, we come to the inequalities , a 9 c 3 a 4 (a 6 + a 9 ) which are equivalent to inequalities (9). Thus, if inequalities (9) hold, then the system (14) has a solution η * i , i = 1, . . . , 5. Taking these values, we have L f h 8 (x) ≤ 0 in polytope Π 1 and h 8 (x) = 0 if that is, on the axis Ox 2 . The unique compact invariant contained in the half axis Ox 2 is the equilibrium point E 0 . Now our assertion is followed from the LaSalle theorem.

1.
Conditions of Theorem 2 do not depend on controls u 3 and u 4 .

2.
If the condition d 2 < 1 holds, then Theorem 2 is true. Indeed, taking into account (3), we get that this condition implies the condition d 1 < 1 and, consequently, conditions (9). Note that the condition d 2 < 1 does not depend on controls, that is, it is satisfied with zero values of controls.
Theorem 3. Suppose that Then all ω-limit sets in R 7 +,0 are located in the invariant plane Proof. We take h 10 = η 1 x 1 + η 2 x 3 + x 4 . Then We obtain L f h 10 ≤ 0 in Π 1 if the following conditions hold: These inequalities have a solution with respect to η 1 , η 2 if the following condition holds: Since the last inequality is satisfied in virtue of (17), then for the corresponding choice of parameters η 1 , η 2 , we have that L f h 10 ≤ 0 in Π 1 and Similarly, arguing as in the previous result, we get the desirable conclusion.
Conditions (20) are met with the appropriate choice of η i if the condition (19) is true. Therefore, choosing the proper parameters η i , we have that L f h 9 ≥ 0 within Π 1 and L f h 6 = 0 if x 3 = x 4 = x 5 = x 6 = x 7 = 0 therein. Next, all trajectories eventually go into Π 1 and remain there. Utilizing the LaSalle theorem with respect to the domain Π 1 , we get that all ω-limit sets of the system trajectories lie in the plane x 1 x 2 . However, there is only one compact invariant set in the plane x 1 x 2 -the equilibrium point E 0 . We conclude that this point is GAS. However, this condition contradicts the condition (9) of LAS of E 0 . This means that the assumption (19) is not true, and we have the condition (18).

Concluding Remarks
In this work, we carry out an analysis of ultimate dynamics of the seven-order Stengel model as mentioned below: • Calculate equilibrium points; • Present local stability conditions; • Find ultimate upper bounds for all variables of this model that define the polytope containing all ω-limit sets.
Our principal contribution is the exploration of wild-type and the possibility of mutant HIV particle eradication, as well as infected cells in the model (1) at the early stage of HIV infection provided the treatment affects only wild-type HIV viruses.
In particular, we establish that the LAS condition to the infection-free equilibrium point E 0 implies the GAS condition to E 0 . These conditions do not depend on control parameters u 3 , u 4 and under the additional assumption that d 2 < 1, the GAS conditions to E 0 do not depend on controls u i , i = 1, . . . , 4, thus, they can be chosen as equal to zero.
To summarize, the system (4) may possess complex behaviour only when conditions (9) are violated. In this case, the location of the attracting set is described in Theorem 1. However, analysis of other features of ultimate dynamics of (4) remains a very difficult task.