Global Dynamics of a Virus-Immune System with Virus-Guided Therapy and Saturation Growth of Virus

We considered a piecewise virus-immune dynamic model to investigate the effectiveness of the HIV virus loads-guided structured treatment interruptions (STIs). To better describe the biological reality, we extended the existing models by taking the carrying capacity of the virus loads into consideration to indicate the saturated growth of virus loads. We initially investigated the sliding dynamics of the proposed model and then obtained the global dynamics of the proposed model. Our main results showed that the system can exhibit very complex and diverse dynamic behaviors including a globally asymptotically stable equilibrium, bistable equilibra, and tristable equilibria, depending on the dynamics of the subsystems and the threshold level. In particular, an interesting result indicated that, with a proper threshold condition, the virus-guided therapy policy can successfully control the virus loads far below its carrying capacity and maintain the activity of the immune system for the case that the effector cells always go to zero without therapy or with continuous therapy. The finding suggested that the optimal strategy should be individual-based due to coexistence of multiple stable steady states, depending on the threshold conditions and the initial levels of viral loads and effector cells of the patients.


Introduction
Highly active antiretroviral therapy (HAART) plays significant roles in improving survival and reducing morbidity and mortality in HIV patients [1,2].Nevertheless, lifelong HAART confronts many complicated problems such as adherence difficulties and the evolution of drug resistance [3][4][5][6].Structured treatment interruptions (STIs), the alternative strategies, have been designed to effectively achieve sustained specific immunity for early therapy in HIV infection.Indeed, STIs are beneficial for some chronically infected individuals who may need to take drugs throughout their lives as well as the reconstruction of patients' immune system when they are not taking the drugs [7].
Since the first HIV patient received the structured antiretroviral therapy in 1999 [8], STIs have been widely investigated by many researchers.Basically, it includes two kinds of treatment regimes: carrying out the therapies at fixed moments or depending on the CD4 cell counts and/or viral loads.For example, some studies designed a scheduled treatment that therapy was initiated only when the CD4 cell counts decreased below 350 cells/ [9][10][11][12].Aiming to maintain CD4 cell counts higher than 350 cells/ and HIV virus loads less than 100,000 copies/, Ruiz et al. designed an experiment to evaluate the safety of CD4 cell and HIV virusguided STIs [13].Many works have been done to compare the STIs with the continuous antiretroviral therapy (ART), but conflicting results have been reported [7,11,[13][14][15][16][17].In particular, Maggiolo et al. made a conclusion that the structured treatment can result in the similar benefits for the patients with the continuous therapies [7], while the SMART study group found that STIs may also endanger the patients' life if their immune responses are not sufficiently stimulated [18].STI is still an important topic since lifelong ART brings great challenges to all HIV positive individuals.Moreover, as immune therapy arises, combining the ART with immune therapy (such as interleukin-2 (IL-2) treatment) may be a good choice for treating HIV patients [19].Modelling and quantifying the efficacy of the combination of STIs and immune therapy remains unclear and falls within the scope of this study.
Mathematical models are useful tools for investigating the dynamics of HIV-1 infection [20,21].Although many studies modelled the continuous therapy for HIV-infected patients [22][23][24], few attempts studied the STIs through mathematical modelling.In 2000, using a standard population model, Bonhoeffer et al. argued the benefits and risk generated from STIs [25].In the study [26], Adams et al. analyzed the efficacy of STIs with fixed moments through the optimal control theory.Further, Rosenberg et al. pointed out that it is needed to include drug resistance and CD4 cell counts when making the structured treatment [27].In 2006, Park et al. developed a mathematical model to investigate the impact of finite times of structured treatment guided by CD4 cell counts on HIV patients [28].Then, Tang et al. [29] formulated a piecewise system to describe the CD4 cell-guided STIs, quantitatively explored STI strategies, and explained some controversial conclusions from different clinical studies.To further model immune therapy, Tang et al. [30] proposed a piecewise virus-immune dynamic model considering the combined antiretroviral therapy with interleukin-2 (IL-2) treatment.The model is given as follows: with where  and  represent the HIV virus loads and the effector cells, respectively.Note that a critical value of virus loads, denoted by   , is set to trigger the combined therapy.That is, only when the virus loads exceed the threshold, we carry out the antiretroviral therapy and immune therapy (such as the interleukin-2 (IL-2) treatment) simultaneously for the patients, with  1 being the elimination rate of HIV virus due to antiretroviral therapy and  2 representing the growth rate of the effector cells due to immune therapy.Here,  is the growth rate of HIV virus which incorporates both their multiplication and death,  is the death rate of the effector cells,  denotes the rate of binding of the effector cells to the virus loads,  represents the rate of inactivation of the effector cells, and /(1 + ) denotes that effector cell multiplication due to immune response has a maximum value as the HIV virus loads become large sufficiently.In [30], the authors mainly discussed the global dynamics of system (1), and they further considered the efficacy of the effector cell-guided or both the effector cell and the HIV virus-guided piecewise therapies [31,32].Their results indicated that HIV virus can be controlled under some conditions with virus or effector cell-guided therapy, but still there are some parameter regions that the threshold policy can not prevent HIV from growing infinitely.As we can see from model (1), they simply assumed the linear growth of HIV virus.However, this is not what the thing looks like.Some clinical facts indicate that there should be "saturation effects" for the growth of HIV virus [33,34].The linear growth of HIV virus is unreasonable.Note that including the saturated growth of the HIV virus makes the model more natural but significantly changes the dynamics of the basic immune-virus system.Little is known about whether HIV virus-guided therapy based on this more natural description of its growth can completely control the growth of HIV and falls within the scope of this study.
The main purpose of this paper is to propose a piecewise model with the saturated growth for virus loads to describe HIV virus-guided therapy and then evaluate the effectiveness of the threshold policy through analyzing the global dynamics of the proposed model.The remaining parts of this paper are organized as follows.In Section 2, we initially propose the model, give some basic definitions of the Filippov system and also briefly conclude the dynamic behaviors of two subsystems.In Section 3, we mainly discuss the sliding mode and sliding dynamics, including the existence of the sliding domain and the pseudoequilibrium.In Section 4, we investigate the global dynamics of the proposed model.Finally, we provide the discussion and some biological conclusions.

Filippov Model and Preliminaries
In this paper, we extend model (1) by including the saturated growth of HIV virus, and hence we propose the following model: with Here, 1/ represents the carrying capacity of virus loads.
In order to consider the trajectory of the Filippov vector field (  1 (),   2 ()), through a point  ∈ Σ, we distinguish the following regions on Σ according to whether or not the vector field points towards it: (a) Crossing region: (b) Sliding region: (c) Escaping region: Further, we can define a scalar differential equation on the sliding domain by employing the Filippov method.Let with 0 <  < 1.Then, we have that describes the dynamics of system (6) on the sliding segment Σ  .In what follows, we give the definitions of two types of equilibria of Filippov system (6).
Definition 1.A point  * is referred to as a real equilibrium of system (6) if Both the real equilibrium and virtual equilibrium are called regular equilibria.
Definition 2. A point  * is referred to as a pseudoequilibrium if it is an equilibrium of the sliding mode of system (6); i.e., the equilibrium of the equation   = () = 0 satisfies ( * ) = 0.
Before exploring the global dynamics of the Filippov system (3), it is essential to have a clear view of the dynamics of two basic subsystems.Thus, we first briefly make some conclusions on the dynamics of subsystem  1 and subsystem  2 .
Next, we briefly investigate the stability of the equilibria of system  1 . 1 is locally stable if (/)/(1 + /) − / −  < 0, which is equivalent to Δ 1 < 0 or Δ 1 > 0,  11 < 0 or Δ 1 > 0, and  12 > 0. Define a function (, ) = 1/, and calculate where / and / are defined in system  1 .Since all the parameters are positive and solutions to system  1 initiating from  2 + are nonnegative,  does not change its sign in the first quadrant.Hence according to Dulac-Bendixson criterion, we get that there are no limit cycles or homoclinic connections for system  1 .Therefore, if there is no positive equilibrium, then the boundary equilibrium  1 is globally asymptotically stable because system  1 is bounded.Further, if Δ 1 > 0 and  11 < 1/ <  12 hold true, it is easy to verify that (1/) > 0, and hence the boundary equilibrium  1 is unstable.Therefore, under this scenario, positive equilibrium  11 is globally asymptotically stable.When there are two positive equilibria of system  1 , the dynamics of system  1 have been shown in [38] that the equilibrium  12 is always an unstable saddle and  11 and  1 are bistable.Then the dynamics of system  1 can be concluded as follows.

Sliding Dynamics
According to the definition of () in Section 2, we have Let () < 0, and one yields It is reasonable to assume that the therapy should be done before virus growing and approaching its carrying capacity; therefore, in this study, we always assume that   < 1/ holds true.Then, if Next, we employ Utkin's equivalent control method introduced in [39] to examine the sliding dynamics in the region Σ  .It follows from () = 0 and the first equation of system (3) that Solving ( 23) with respect to Φ yields Then substituting Φ into the second equation of system (3) gives Then, the sliding dynamics of system (3) can be described as (25).Let   (  , ) = 0, and we obtain two roots This indicates that there are two possible pseudoequilibria  0 = (  , 0) and   = (  ,   ).In order to verify whether they are located at the sliding segment, we consider two situations 0 <   <  2 and  2 <   <  1 .
It is easy to see that  0 can not be a pseudoequilibrium when 0 <   <  2 , while it is always a pseudoequilibrium for  2 <   <  1 .Thus, we only need to analyze the existence of   .First, we prove the following lemmas.
(ii) There is a unique positive equilibrium for subsystem  1 (denoted by Case A22); that is, For such case, it follows from Lemmas 6 and 7 that   is feasible for  21 <   <  11 .Especially, if  11 >  2 , then there are two pseudoequilibria  0 and   for  2 <   <  11 , and there is no pseudoequilibrium for 0 <   <  21 .
(iii) There are two positive equilibria for subsystem  1 (denoted by Case A23); that is, Here, we consider three subcases:  Case A3.There are two positive equilibria for subsystem  2 .Under this case, there can not be a unique positive equilibrium for subsystem  1 according to Lemma 5.Then, consider the following: (i) There is no positive equilibrium for subsystem  1 (denoted by Case A31); that is, It follows from Lemma 6 that   is feasible for  21 <   <  22 .
Having discussed the existence of the pseudoequilibria, we take the derivative of   with respect to  to analyze the stability of the pseudoequilibria.Note that  Therefore, if ( 1 ) holds true, then the pseudoequilibrium   is locally stable on the sliding segment when it is feasible, and  0 is unstable; if ( 2 ) holds true, then   becomes infeasible and the pseudoequilibrium  0 is locally stable on the sliding segment when it exists.

Global Dynamics of System (3)
In this section, we will investigate the global dynamics of Filippov system (3).Initially, we shall rule out the existence of limit cycle for system (3).Here, we take Case A11 as an example to illustrate that there is not any limit cycles.The detailed proof is given in the appendix.Note that, for the other cases, we can prove the nonexistence of limit cycle for system (3) using the similar methods, and here we omit the detailed proof.

Lemma 8. For Case A11, there is no limit cycle for Filippov system (3).
Corresponding to the discussion of the sliding dynamics, we then consider the following cases to examine the global dynamics of system (3).Case A1.There is no positive equilibrium for subsystem  2 .
(i) There is no positive equilibrium for subsystem  1 (i.e., Case A11).Then, we discuss the global dynamics of the Filippov system by considering the following subsituations: (a) If 0 <   <  2 , then we know that  2 is a real equilibrium which is stable while  1 is a virtual equilibrium which is unstable.As we mentioned above, there is no limit cycle for system (3), and moreover there is not any other stable equilibrium, thus the local stability of  2 indicates that  2 is globally asymptotically stable (shown in Figure 2(a)).(b) If  2 <   <  1 and ( 2 ) hold true, then we have that  2 becomes a virtual equilibrium.It follows from the discussion in Section 3 that the pseudoequilibrium  0 is stable on the sliding segment Σ 1  .Similarly, system (3) does not have any limit cycle.Therefore, the local stability of  0 indicates that  0 is globally asymptotically stable (shown in Figure 2(b)).
(c) If  2 <   <  1 and ( 1 ) hold true, then there are two pseudoequilibria  0 and   on the sliding segment Σ 1  .As we proved in Section 3,   is locally stable while  0 is unstable.Similar to subcase (b), we have that   is globally asymptotically stable (shown in Figure 2(c)).
Based on the above discussions, the global dynamics of system (3) can be concluded as follows.
(ii) There is a unique positive equilibrium for subsystem  1 (i.e., Case A12).If 0 <   <  2 ,  1 and  11 are virtual equilibria, and  2 is a real equilibrium which is stable.As mentioned above, under this scenario, there is no limit cycle or no pseudoequilibrium for Filippov system (3), and thus  2 is globally stable.If  2 <   <  21 ,  2 becomes a virtual equilibrium, and there emerges the stable pseudoequilibrium  0 on the sliding segment Σ 1  .If  21 <   <  11 , the same results hold compared with the case  2 <   <  21 if ( 2 ) holds true, while there are two pseudoequilibria  0 (unstable) and   (stable) on Σ 1  if ( 1 ) holds true.Furthermore, If  11 <   <  1 ,  11 becomes a real equilibrium which is locally stable.There is a pseudoequilibrium  0 on Σ 1  which is unstable.Then, the global dynamics of Filippov system (3) are concluded as follows.

)).
(iii) There are two positive equilibria for subsystem  1 (i.e., Case A13).If there is only one steady state (boundary equilibrium or positive equilibrium or pseudoequilibrium), then it is globally stable.We can analyze the dynamics of system (3) under the same scenario using the similar methods.However, there is a different situation for such case.That is, if  12 <   < min{ 1 ,  22 } holds true, then  11 is a real equilibrium; hence it is locally stable.Furthermore,  0 is a locally stable pseudoequilibrium if ( 2 ) holds true, but if ( 1 ) holds true,  0 lose its stability and   becomes a locally stable pseudoequilibrium.In conclusion, if  12 <   < min{ 1 ,  22 },  11 is bistable with  0 (or   ) for ( 2 ) (or ( 1 )) being satisfied.Based on the discussion above, the global dynamics of system (3) for Case 13 can be concluded as follows.
Cases A2 (A21, A22, and A23) and A31.All the dynamics of system (3) for four cases are similar to the former three cases, and thus we just omit the proof here.The global dynamics of these cases are summarized in Table 1 with the dynamical behaviors shown in Figures 5-8, respectively.Case A32.There are two positive equilibria for both subsystems  1 and  2 , respectively.And the horizontal compartments of these equilibria satisfy  21 <  11 <  12 <  22 .In terms of the existence and stability of the positive equilibria, the dynamics of the subsystems are similar to the focused situation in study [30].Thus, based on the relationships of   with respect to  21 ,  11 ,  12 , and  22 , we can analyze the dynamics of system (3) through the following subcases.
Thus, based on the above discussions, we have the following conclusions.

Conclusion and Discussion
Based on the basic immune-virus model, we considered the threshold policy depending on virus loads for the HIV-infected patients and obtained a nonsmooth system.Compared with the model in [30], we extend it by taking the saturated growth of virus into consideration, and hence our model is better in line with the biological reality.Meanwhile, our model also includes a piecewise elimination rate of HIV virus and growth rate of the effector cells to indicate that the antiretroviral therapy and interleukin-2 (IL-2) treatment are carried out simultaneously once the HIV virus loads increase above a threshold level.
We first briefly concluded the dynamics of two subsystems.Based on the properties of the subsystems, we discussed the sliding mode and sliding dynamics with different threshold conditions.In more detail, the sliding segments are   there is always a pseudoequilibrium  0 on the sliding segment Σ 1  .And we further discussed the existence of the other pseudoequilibrium   by considering eight different cases.Correspondingly, we then investigated the global dynamics of system (3) for various cases.Our main results show that the system can exhibit very rich dynamic behaviors: (a) one of the equilibria is globally asymptotically stable, which can be the pseudoequilibrium  0 or the pseudoequilibrium   or the boundary equilibrium  2 of subsystem  2 or the positive equilibrium  11 or  21 of subsystem  1 or  2 ; (b) two equilibria are bistable, particularly, which can be that the positive equilibrium  11 is bistable with  0 (or   or  2 ), or the boundary equilibrium  2 is bistable with the pseudoequilibrium   (or the positive equilibrium  21 ); (c) three equilibria are locally tristable; that is, the pseudoequilibrium   , the boundary equilibrium  2 , and the positive equilibrium  11 are stable for  12 <   <  22 .
Global dynamics of our proposed model means that there is the case that the boundary equilibrium is globally stable for both the subsystems.Note that the immune-virus system without any therapy is actually system  1 (the free system) while the immune-virus system with continuous therapy can be described by system  2 (the control system).Comparing system  2 with system  1 , we can find that the continuous therapy can shrink the final density of the virus (i.e.,  2 <  1 ).It follows that the continuous therapy failed to maintain the activity of the immune system since the effector cell goes to zero finally.However, threshold policy indicated that, with a proper threshold condition, the STIs can not only successfully control virus loads to a previously given level which can be far below its natural carrying capacity (i.e.,  1 ) but also maintain the activity of the immune system (corresponding to the case that the pseudoequilibrium   is globally stable), as Theorem 9 suggested.Moreover, according to the dynamics of subsystem  1 for Case 21 we can see that virus loads will finally approach its carrying capacity if there is no therapy for the patients.With the continuous therapy, the virus loads of the patients can be successfully controlled far below its carrying capacity and their immune system can always maintain the activity; that is, equilibrium  21 is globally asymptotically stable for system  2 .For the virus-guided threshold policy, if the threshold level is set to satisfy  21 <   <  2 , then the pseudoequilibrium   is globally stable, and   will tend to  21 as   →  21 .This means that the threshold policy can obtain the similar effects for the patients compared with the continuous therapy in terms of virus loads being controlled, which is in agreement with the clinical study [7].
It is worth mentioning that if we also assume that the two subsystems have two positive equilibria (i.e., Case A32), our model results, due to the saturated growth in model equations, show that the virus can not grow to infinity.Under this situation, we also found that the dynamics of system (3) are similar to those of the control system (system  2 ) when the threshold condition is relatively small (i.e., 0 <   <  21 ); that is, the positive equilibrium  21 and the boundary equilibrium  2 are bistable.And if the threshold is relatively high (i.e.,   >  22 ), the dynamics of system (3) are similar to the free system  1 .Moreover, for middle values of the threshold (i.e.,  12 <   <  22 ), then equilibria   ,  11 , and  2 are locally tristable.Therefore, the bistability and the tristability suggest that the optimal strategy should also be individual-based by taking the initial conditions of the patients and the threshold conditions into consideration.which contradicts (A.4).Thus, we preclude the existence of the limit cycle C surrounding the sliding segment  1  2 .The proof is completed.
Professor Robert Cheke for his valuable comments and for his generous help in polishing the manuscript.

Figure 1 :
Figure 1: Phase plane x-y for Filippov epidemic model (3), showing the switching line ( =   ), sliding segment  1  2 , and the diagrams  1 and  2 split from the limit cycle C by the switching line.

Figure 8 :∫
Figure 8: Phase plane of Filippov system (3) for Case 31 with different threshold conditions.(a) The regular equilibrium  21 and the boundary equilibrium  2 are bistable, here   = 3; (b) the pseudoequilibrium   and the boundary equilibrium  2 are bistable, here   = 5.5; (c) the boundary equilibrium  2 is globally asymptotically stable, here   = 7; (d) the pseudoequilibrium  0 is globally asymptotically stable.Parameter values are chosen as follows:  = 2.55,  = 0.1, and  2 = 0.1, and other parameter values are shown in Figure 5.