THE FOUR-DIMENSIONAL KIRSCHNER-PANETTA TYPE CANCER MODEL: HOW TO OBTAIN TUMOR ERADICATION?

. In this paper we examine ultimate dynamics of the four-dimen-sional model describing interactions between tumor cells, eﬀector immune cells, interleukin -2 and transforming growth factor-beta. This model was elaborated by Arciero et al. and is obtained from the Kirschner-Panetta type model by introducing two various treatments. We provide ultimate upper bounds for all variables of this model and two lower bounds and, besides, study when dynam- ics of this model possesses a global attracting set. The nonexistence conditions of compact invariant sets are derived. We obtain bounds for treatment pa- rameters s 1 , 2 under which all trajectories in the positive orthant tend to the tumor-free equilibrium point. Conditions imposed on s 1 , 2 under which the tumor population persists are presented as well. Finally, we compare tumor eradication/ persistence bounds and discuss our results.


1.
Introduction. The Kirschner-Panetta equations [2] have a great influence on the modelling tumor dynamics under immunotherapy. These equations describe dynamics of interactions between tumor cells, effector immune cells and interleukin-2 (IL-2). One of the promising generalizations of this model is obtained by incorporating in these equations yet another differential equation characterizing the dynamics of the suppressor cytokine, transforming growth factor-β (TGF-β), [1]. It is wellknown [16] that TGF-β may display both inhibitory activity and stimulating activity on the growth of most of cells depending on type of cells, their differentiation and activation state. The production of TGF-β by tumor cells greatly challenges the immune system through the promotion of angiogenesis, enhancing tumor growth and metastasis. Tumors can evade immune surveillance by secreting various immunosuppressive factors including (interleukin-10) IL-10 and TGF-β, [8,17]. It was indicated in [17] that activation of the immunosuppressed immune system by cytokine IL-2 therapy is a possible strategy to limit the malignant immuno-modulatory activities of TGF-β. This type of therapy may be applied alone or in combination with other immunotherapeutic approaches and is used in the clinical practice, [16]. This possible approach to cancer treatment is formalized in this paper by introducing into the model from [1] two treatment parameters s 1,2 . Parameters s 1,2 are included in the equations for the same cells populations as in the Kirschner-Panetta (KP)-model. We recall that s 1 is the treatment term that represents an external source of IL-2 injected into the system; s 2 is the treatment term that represents an external source of effector cells such as LAK (lymphocytes activated killer cells) or TIL (tumor-infiltrating lymphocytes). Following to [1] we consider that s 1,2 are constants. So we come to the following model in the non-dimensional forṁ In equation (1) x(t) describes the number of tumor cells at the moment t; y(t) describes the concentration of effector molecules at the moment t; z(t) describes TGF-β's immuno-suppressive and growth stimulatory effects in the single tumorsite compartment; w(t) describes the number of immune cells at the moment t.
All parameters are supposed to be positive excepting s 1,2 which are nonnegative. In the first equation r is the cancer growth rate; a is the cancer clearance term. The proliferation of tumor cells due to the response to TGF-β is denoted by p 2 and is modeled by Michaelis-Menten kinetics. In the second equation p 3 is the rate of IL-2 production in the presence of effector cells; g 4 is half-saturation constant; α is a measure of inhibition. In the third equation p 4 is maximal rate of TGF-β production; τ is the critical tumor cells population at which angiogenesis switch occurs. In the fourth equation c is known as the antigenicity of the tumor which measures the ability of the immune system to recognize tumor cells; γ is inhibitory parameter; µ 1 is the death rate of immune cells; p 1 is the proliferation rate of immune cells; q 1 is the rate of anti-proliferative effect of TGF-β; q 2 is half-saturation constant. More details concerning these parameters are contained in [1].
We notice that if we put z = 0, p 4 = 0 in (1) we get the KP-model which was created for studying the immune response to tumors under special types of immunotherapy.
Dynamics of the KP-model has been studied in [1,2,3,14] and some others. In particular, ultimate upper and lower bounds for state variables of the KP-model have been derived in different cases. The main result of [14] consists in global asymptotic tumor clearance conditions obtained under various assumptions imposed on the ratio between the proliferation rate of the immune cells and their mortality rate. To the best of the authors' knowledge, up to now there have not been published any results concerning rigorous dynamical analysis of (1). We notice that the system (1) introduced for the case s 1 = s 2 = 0 was explored in [1] only by means of numerical simulations of its dynamics.
In this paper we consider the tumor growth system pertained to the broad class of life sciences models which possess the following characteristic feature: there is a tumor-free equilibrium point, which is the most preferable state of the system. From the biological point of view some deviations from this equilibrium point can be dangerous and cause fatal outcomes.
Therefore, the control goal is to return the system to the indicated equilibrium point and keep it in a sufficiently small neighborhood of this state. It can be done by different treatment injections.
This article establishes the existence conditions for the positively invariant polytope that has a biological meaning. Further, two types of conditions are found: the first one is the tumor persistence (impossibility to achieve the control goal), another one is the tumor eradication (possibility of global asymptotical stabilization at the tumor-free equilibrium point).
The goal and the novelty of this work consists in studies of ultimate dynamics (1) in case of applied treatments s i , i = 1, 2. Namely, we find upper and lower ultimate bounds for all variables of the system (1) and establish conditions under which (1) is dissipative in the Levinson sense; 2) we propose the nonexistence conditions of compact invariant sets in the positive orthant; 3) we deduce the global asymptotic tumor eradication conditions; 4) we describe the tumor persistence conditions.
In other words, the three most important results for understanding of any dynamical system's behavior are actually proved: the existence conditions of the global system's attractor; the coincidence of this attractor with the tumor-free equilibrium point; the presence of a local attractor of the system that does not contain the tumor-free equilibrium point, but attracts almost all ("perturbed") trajectories of the system.
Our approach is based on the localization method of compact invariant sets in which the first order extremum conditions are utilized, see [4,5,6]. We also mention that earlier this method has been successfully utilized in studies of various cancer tumor growth models, see e.g. [7,9,10,11,12,13,14,15] and references therein.
The remainder of the paper is organized as follows. In Section 2 we briefly present useful results. In Section 3 under some condition we obtain formulae for a polytope containing all compact invariant sets. This polytope provides us ultimate upper and lower bounds for all variables of the system (1), see Theorem 3.4. In Section 4 under the same condition as in Theorem 3.4 we show in Theorem 4.1 that this polytope contains the attracting set of the system (1). In Section 5 we present the nonexistence conditions of compact invariant sets in R 4 +,0 ∩ {x > 0}, see Theorem 5.1. Further, in Section 6 using results of Theorems 3.4; 4.1; 5.1 we derive conditions of global asymptotic stability with respect to the tumor-free equilibrium point (TFEP), see Theorem 6.1. In the latter theorem we provide bounds for treatments parameters s i , i = 1, 2, for which the global asymptotic tumor eradication process may be observed. In Section 7 we describe the persistence tumor conditions which are compared with tumor eradication bounds of Theorem 7.1 in Section 8. Section 9 contains the concluding remarks.
In what follows, we examine dynamics of (1) in the positive orthant 2. Some useful results. We consider a nonlinear systeṁ where v is a C 1 −differentiable vector field; x ∈ R n is the state vector. Let h(x) be a C 1 −differentiable function such that h is not the first integral of (2). By Assume that we are interested in the localization of all compact invariant sets contained in the set U . Further, we define Assertion 1 [5,6] For any h(x) ∈ C 1 (R n ) all compact invariant sets of the system (2) located in U are contained in the localization set K(h; U ) defined by the formula Then for any τ 0 > 0 the extended localization set then every trajectory of the system (2) goes into the setK(h; U ; τ 0 ) in finite time.
3. Formulae for a polytope containing all compact invariant sets. Here we find localization sets for the system (1). Let f be a vector field of this system. Lemma 3.1. All compact invariant sets in R 4 +,0 are located in the set Proof. We apply the function h 1 = x and get that and h 1,inf (R 4 +,0 ) = 0. On the set S(h 1 ) ∩ R 4 +,0 ∩ {x > 0} the inequality holds. Using the last inequality to calculate the supremum we obtain Let η = min(g 4 ; 1)ap −1 3 . Lemma 3.2. All compact invariant sets in R 4 +,0 are located in the set Proof. We apply the function h 2 = x + ηy and compute In R 4 +,0 we get the inequality because z 1+z ≤ 1 and η was chosen earlier such that After calculating the supremum we obtain that Now let us use the function h 3 = y and compute Therefore, and we come to the desirable conclusion.
Lemma 3.3. All compact invariant sets in R 4 +,0 are located in the set Proof. We apply the function h 4 = z and get that which leads to the desirable conclusion.
Theorem 3.4. If then all compact invariant sets are located in the polytope .

ALEXANDER P. KRISHCHENKO AND KONSTANTIN E. STARKOV
Proof. We apply the function h 5 = w and transform equation L f h 5 = 0 in the equality w(µ 1 + y 1 + y ( q 1 z z + q 2 − p 1 )) = cx 1 + γz + s 1 . Therefore, and all compact invariant sets are located in the set {w min ≤ w ≤ w max } ∩ M = Π.

Corollary 1.
If then all compact invariant sets are located in the set

4.
On the dissipativity in the sense of Levinson. Below we shall establish conditions under which the system (1) is dissipative in the sense of Levinson. Here we recall that the system (2) is called dissipative in the sense of Levinson if there exists r > 0 such that for any x ∈ R n we have that lim t→∞ sup |ϕ(x, t)| < r; here |ϕ(x, t)| is the Euclidean norm of the solution ϕ(x, t) of the system (2) starting in time t = 0 at the point x ∈ R n . In this case there exists a bounded set which attracts any trajectory in R n .
Theorem 4.1. If condition (7) is fulfilled then the system (1) is dissipative in sense of Levinson in R 4 Proof. Firstly, we note that extended localization setŝ where τ i > 0, i = 1, 2, 3, have the form of the set (3) and, by Assertions 2,3 (see remark below), are positively invariant and every trajectory goes into these sets in finite time. Therefore, the intersection of these setsM =K 1 ∩K 2 ∩K 3 is a positively invariant set and every trajectory goes into this set in finite time.
Next, if condition (7) is fulfilled then for some sufficiently small τ 2 > 0 we get We fix such value of τ 2 and find the localization set in the set {h 1 ≥x max } ∩ R 4 +,0 the localizing function h 1 = x is equal tox max + ∆ 1 > 0, where ∆ 1 ≥ 0, and therefore, in the set {h 4 ≥ẑ max }∩R 4 +,0 the localizing function h 4 = z is equal toẑ max +∆ 3 , where ∆ 3 ≥ 0, and therefore, in the set {h 5 ≥ŵ max }∩M the localizing function h 5 = w is equal toŵ max +∆ 4 > 0, where ∆ 4 ≥ 0, and therefore, 5. The nonexistence conditions of compact invariant sets in R 4 +,0 ∩{x > 0}. Under condition (8) all compact invariant sets lying in the set R 4 +,0 ∩ {x > 0} are contained in the set O 1 := M 1 ∩ {x > 0} (see Corollary 1). Below we apply localizing function h 1 = x and show that its derivative is negative in the set O 1 if some inequality holds. Therefore, in the case (8) this inequality is a nonexistence condition of compact invariant sets in the set R 4 +,0 ∩ {x > 0}. This condition may hold both in case of the existence of the TFEP and its nonexistence. It means the nonexistence of bounded tumor persistence dynamics, for example, the tumor dormancy. As a corollary, we describe the property of the nonexistence of periodic orbits and tumor persistence equilibrium points (TPEPs) in some range of model and treatment parameters.

ALEXANDER P. KRISHCHENKO AND KONSTANTIN E. STARKOV
Theorem 5.1. Suppose that (8) and where Then there are no compact invariant sets in the set R 4 Proof. Let us apply the function h 1 = x and find that In order to find w 0 we consider and getη Therefore ifη (0) ≤ 0 then w 0 =η(0) and ifη (0) > 0 then w 0 =η(x * ).
6. Global asymptotic stability respecting the TFEP. If the system (1) has the TFEP The TFEP is asymptotically stable if r < aw 1 i.e.
Proof. If conditions of this theorem are fulfilled then all trajectories of the system (1) go into bounded positively invariant setK 4 ; the TFEP exists and is asymptotically stable. Therefore, in order to prove Theorem 6.1 it is sufficient to show that the TFEP is the unique compact invariant set of the system.
The system (1) has no compact invariant sets in R 4 +,0 ∩ {x > 0} (see Theorem 5.1). The TFEP is the only compact invariant set of the system in the invariant plane {x = 0}. Indeed, let us consider the system restricted on this planė Next, applying localizing functions x; y; w we obtain localization sets for compact invariant set of the system (12) Now let us prove that the system (1) has no compact invariant set C for which C ∩ R 4 +,0 ∩ {x > 0} = ∅ and C ∩ {x = 0} = ∅. Indeed, otherwise C ∩ {x = 0} = E 1 and there exists a point P ∈ C ∩ R 4 +,0 ∩ {x > 0}. In this case the α-limit set of the trajectory starting at the point P is a nonempty compact invariant set D and E 1 / ∈ D because otherwise the point E 1 is not stable. Therefore, the nonempty compact invariant set D is a subset of R 4 +,0 ∩ {x > 0} and the statement of the theorem follows from this contradiction with Theorem 5.1.
2. Further, we provide conditions (8) and (9) under which there are no compact invariant sets in the set R 4 +,0 ∩{x > 0}. As a result, there is no conditions for tumor dormancy. In particular, the system (1) has neither TPEPs nor periodic orbits.
3. We find conditions (7); (9) and (11) under which the TFEP attracts all trajectories in R 4 +,0 . The biological sense of this behavior consists in asymptotic eradication of tumor cells which means that after a while the tumor cells population will be under control.
4. Tumor eradication and tumor persistence bounds are compared in Section 8. One can point to the following essential difference of dynamics of (1) in cases s 1,2 = 0, [1], and in case s 1,2 > 0 under varying antigenicity c. Namely, it was noticed in [1] that there are many negative scenarios including uncontrolled tumor growth and damped oscillations around the TPEP, which corresponds to tumor dormancy. In our work, because of the proper assignment of treatment parameters s 1,2 satisfying Theorem 6.1 for given model parameters a, k, r, µ 1 , g 4 , p 1 , p 2 , p 3 , p 4 one may achieve tumor eradication regardless the value of c, where c > 0. However, the tumor persistence bound s per 1 depends on parameter c. All assertions are formulated in terms of simple algebraic inequalities imposed on parameters of the model and treatments. These inequalities are stable for sufficiently small perturbations caused by imprecise knowledge of parameters' values which is convenient in applications.