QUALITATIVE ANALYSIS OF KINETIC-BASED MODELS FOR TUMOR-IMMUNE SYSTEM INTERACTION

. A mathematical model, based on a mesoscopic approach, describing the competition between tumor cells and immune system in terms of ki- netic integro-diﬀerential equations is presented. Four interacting components are considered, representing, respectively, tumors cells, cells of the host en- vironment, cells of the immune system, and interleukins, which are capable to modify the tumor-immune system interaction and to contribute to destroy tumor cells. The internal state variable (activity) measures the capability of a cell of prevailing in a binary interaction. Under suitable assumptions, a closed set of autonomous ordinary diﬀerential equations is then derived by a moment procedure and two three-dimensional reduced systems are obtained in some partial quasi-steady state approximations. Their qualitative analysis is ﬁnally performed, with particular attention to equilibria and their stability, bifurca- tions, and their meaning. Results are obtained on asymptotically autonomous dynamical systems, and also on the occurrence of a particular backward bifur- cation.


1.
Introduction. In the last several years many papers have dealt with the problem of devising reliable dynamical models of tumor development [1,3,11,12,13,14,16,17,29]. A large number of such papers make use of systems of ODEs with Lotka-Volterra or Verhulst (logistic) terms for describing the interactions between malignant and immune cells. In spite of their simplicity, these models of tumor growth and possible remission can reasonably describe the different dynamics of cancer development [28]. The most assessed tools in the literature for modeling tumor dynamics are analysis of the equilibrium points, bifurcation diagrams, inspection of the phase plane, or of the phase space, when the competition process is mediated by the presence of additional participating populations [13]. This allows to identify conditions which are critical for tumor growth. Often it can be shown that by changing the values of some control parameter the domain of attraction of the tumor-free equilibrium can be enlarged, and such domain represents a safety region, since, for any initial condition in that region, tumor is annihilated by the immune response. Indeed, when dealing with this kind of problems, it is of primary interest to determine what happens when a (voluntary or undesired) perturbation is injected into the system, moving the trajectory away from the equilibrium point. Interaction terms are typically quadratic, to mimic simple binary interactions between cells of the various populations.
Nonlinear models of quadratic type for tumor dynamics are also obtained within the frame of more detailed and sophisticated theories, and, among them, kinetic approaches have become quite popular, and have proved to be significantly effective and reliable in their predictions. Just to quote few additional more recent contributions we can mention, without pretending to be exhaustive, the papers [3,5,11,20,23,24]. This research line follows the stream of evolution of dominance in population dynamics, where dominance represents any possible internal state, attribute, activity, or performance capability possessed by single individuals. Here binary individual encounters at microscopic level are described by stochastic models leading to Chapman-Kolmogoroff equations in the frame of the theory of Markov processes [21]. The approach is essentially the same as for the derivation of the nonlinear Boltzmann equation of gas kinetic theory [10]. The dependent variables to be investigated are the "dominance" distribution functions, whose first few state moments provide the macroscopic observables. As typical of kinetic theory, exact evolution equations at macroscopic level may be derived for the above physical quantities by taking moments of the microscopic nonlinear integro-differential equations, but the resulting set of differential equations turns out not to be closed. A kinetic approach to immunology problems is motivated not only by the better insight allowed by a deeper description, but also by the fact that the stage of the early growth of a tumor belongs to the so-called free cells regime, in which tumor cells are not yet condensed in a macroscopically observable spatial structure, and interactions between tumor and immune system occur at a cellular level. This stage is particularly important since the competition between tumor cells and immune system can still lead to the depletion of the tumor. At the same time, spatial effects are of minor importance, to leading order, in the balance equations, which implies considerable simplifications in the analytical investigation.
In the present paper a four species model, proposed and validated already in the literature, will be considered [1], in which the competition tumor-immune system is mediated by the presence of the host environment (other cells of the body) and by an additional chemical species of interleukins [4], considered, like the others, as a population of individuals, capable to enhance the immune response without destroying tumor themselves. The main properties and features of the adopted kinetic model are recalled in the next Section for the readers' convenience. In this paper the model is generalized by the addition of extra mechanisms, not of mutual interaction nature, which may affect the specific dynamics of the species, like birth and death processes. Possible different types of immune cells are here represented, as typical in the pertinent literature [3,4,11], by a single population, aiming at a simple description capable to qualitatively reproduce the overall basic behaviour of the immune defense. They could be included at the price of additional technical, but not essential, difficulties. Analogously, it is well known that there exists a great variety of interleukins, of different origin and with quite different behaviour and pecularities (see for instance the seminal paper [22] and the recent work [26]), but in our mathematical model we are mainly interested in capturing their overall effect as immune system supporters, which will be described in a statistical average manner by single interaction rates. For a numerical solution of the resulting integro-differential system with quadratic nonlinearities, a suitable discretization technique is needed anyhow. In this way, the description of the evolution of the various components involves a finite number of key macroscopic parameters, deduced appropriately from the actual collision frequencies and probability distributions characterizing the microscopic interactions, which are instead functions of a continuous kinetic variable, and would be quite hard to determine by comparison with experiments. All those microscopic functions will be kept arbitrary in the general presentation of the model, and will try to cover all binary interactions of any type that might occur among different cells, without having in mind any specific biological problem. In this work a discretization is achieved by integration over partial ranges with respect to the kinetic (state) variable, and grouping together all individual cells with the value of state in the same range to form single separated populations. Such achievement is made possible by technical simplifying assumptions on the microscopic interaction parameters, taken to be constant or piecewise constant with respect to their state variables, which are certainly crude, but account for the interaction mechanisms at least in an average way, and allow a much deeper analytical investigation. In this way, in fact, a closed set of autonomous ODEs is derived, representing a sort of macroscopic continuity equations in the sense of kinetic theory. The qualitative analysis of the evolution problem can then be performed in the well established framework of the theory of dynamical systems [18], as well as the study of the bifurcations occurring on the whole range of variation of parameters [25]. Extensive numerical simulations (very partially shown) have been performed in order to test and improve analytical predictions, by using random selected values of the dimensionless parameters, aiming mainly at analyzing in depth the essential features of the model rather than at focusing on the numerical ranges of major immunological interest. In the spirit of the present work, we keep thus generic the numerical values for our numerical simulations, so that experimentalists can single out the ranges of their interest and see the effects predicted there by our model. The paper is organized as follows. In the next Section, after discussing, at a formal level, kinetic equations for a population of different cells with conservative, destructive, and proliferative events, we proceed next to their specialization to the considered four population model of tumor-immune system competition. In the same Section 3 we reduce the problem to a four dimensional dynamical system and single out the dimensionless physical parameters that are crucial in the evolution. In the following two Sections we further specialize such a dynamical system to two limiting situations of most practical interest, in which the dimension of the phase space reduces to three. These are the cases on which our attention is mainly focused, and where we obtain the main original results of this work. However, in order to control reliability of such reductions, comparative results are plotted for the time evolution of the various species in the 3D and 4D environments. Section 4 is concerned with the case where the host environment is a sort of infinite background whose state is not affected by the process going on. The resulting reduced three-dimensional system of ODEs for tumor cells, immune system and interleukins can be investigated in the framework of the theory of the asymptotically autonomous differential systems [30], and we will show that its asymptotic behaviours can be deduced from an easier two dimensional limit system. In Section 5 the role of background is played by interleukins, that have reached a quasi-steady state condition. A remarkable feature of this latter reduced system for the interactions between tumor, immune system and host environment is the presence of a backward bifurcation, usually related to epidemic models [19], with reversed stability of the colliding equilibria.
2. Balance equations at cellular level. As anticipated in the Introduction, we consider a system of N = 4 different populations, labeled by an index i, each individual (cell) being endowed with an internal state variable (activity) u, ranging in the real interval (−1, 1). Such a state variable is not necessarily related to any biological property, but, in our mathematical model, denotes simply an attribute or competing capability labelling the rank of an individual inside its species, namely its effectiveness in performing the job that the species is expected to do, and in prevailing when interacting with other species. We assume that only binary interactions are effective in the evolution, and restrict ourselves to space homogeneous conditions. A detailed knowledge of the state of the whole system is provided by the four distribution functions (densities in phase space) f i (u, t), smooth non-negative functions, from which one can deduce the cellular densities n i (the actual observable macroscopic quantities of practical interest) simply by Balance equations in phase space may be derived by equating the rate of change at time t of the i-th populations in the elementary activity interval du to the corresponding net production rate (gain minus loss) due to all events of any nature taking place in the system, including any effect coming from external sources, treatments, spontaneous mortality, and so on. An important contribution to the exchange rates comes of course from the mutual interactions among cells, and these can be described by suitable pairwise "collision" operators Q ij , representing the effects on cells of type i of binary encounters with those of type j, in a way that closely resembles models and methods of gas kinetic theory [10], where individuals are molecules, state variable is velocity, and interactions are actual mechanical collisions. We may write the kinetic equations in the form where J i collects all contributions of events different from cellular interactions. As a first step, we build up the interactive operators Q ij . A significant difference with respect to gas dynamics is that encounters are not conservative, but, as typical of other disciplines, like transport theory [15], they may lead to disappearance or proliferation of a participating populations. We shall obtain however integral operators of Boltzmann type by resorting to an equivalent probabilistic formulation [6] in terms of suitable interaction probabilities per unit time (collision frequencies) and creation kernels. More precisely, let η ij (u, v) = η ji (v, u) ≥ 0 denote the collision frequency for a conservative encounter between an i-th cell with activity u and a jth cell in state v, and let ψ ij (u, v; w) ≥ 0 represent the probability density that, after this interaction, the i-th cell ends up in the state w, with the obvious normalization Similarly, we shall denote by d ij (u, v) the collision frequency of an (i, u) cell with a (j, v) cell in an encounter which is not conservative for the population i, and by µ ij (u, v) ≤ d ij (u, v) the reduced collision frequency which is relevant to proliferative interactions only. For the latter events, the expected density of i cells which end up in the state w will be labeled by ε ij (u, v; w), and the integral provides the average number of i cells generated in the proliferative encounter (i, u)-(j, v). Such a number is in general greater than unity. At this point, the count of the number of gains and losses leads to the explicit expression for the kinetic "collision" operator in its usual nonlocal form of integral type. Kinetic equations (2) are explicit once all previous probabilities, as well as all physical specific parameters making up the non-interactive operators J i , are known or appropriately modeled. For practical purposes, one is mainly interested in the evolution of the macroscopic quantities n i , and hopefully a set of ordinary differential equations could be obtained from the integro-differential equations (2) by integration with respect to the state variable u. The result, however, is not closed in general, since cell densities do not factor out directly from the integrals. For instance, integration of (3) over u ∈ (−1, 1) yields, for i = 1, . . . , 4, where of course conservative interactions are not influential, and the positive or negative contribution to n i of the general cells of type j depends on the sign of µ ij m ij − d ij . If such parameters were constant, the collision contribution (4) would reduce to a quadratic form in the densities n i .
We proceed now to the formulation of a kinetic model for the considered problem by specifying, in a very simple but yet realistic manner, inspired by their own physical meaning, the probabilistic quantities in (3), as well as the additional operator J i in (2). The present biological model is in the frame of a research strategy established several years ago [4] and further developed by many authors, and represents a significant generalization to a much more complicated scenario of a similar approach proposed in [20].

3.
A kinetic model for tumor-immune system competition. The four populations making up the physical system are assumed to represent, respectively, tumor cells, cells of the host environment, cells of the immune system, and interleukins, labeled by an index i increasing from 1 to 4. As per the pertinent literature, the latter population plays a role in the overall interaction and contributes to the destruction of tumor by strengthening the action of the immune system. The value of the state u of each cell is a measure of its capability of prevailing in binary interaction, and we will call active all cells with positive state, and passive those with a negative one. Collision frequencies and transition probabilities are specialized as follows, where for simplicity the former will be also taken as positive constant in the domain where they do not vanish (resulting thus piecewise constant), which resembles the popular Maxwell molecule assumption of rarefied gas dynamics [10].
A tumor cell is destroyed by interaction with an active cell of the immune system, but proliferates in interactions with passive immune cells. These interactions are conservative for the immune system, whose activity however always decreases, and is changed from positive to negative in the former event. This is quantified by (5) where U denotes Heaviside function.
An interaction between a tumor cell and a cell of the host environment always ends up with tumor proliferation, namely In addition, the host environment is supposed to be endowed with a self-consistent control mechanism which tends to establish, for an optimal functioning, a given distribution f * 2 (u) of the host cells, with a strength depending linearly on the instantaneous deviation through a rate parameter ν 2 (u). In other words, An encounter between a cell of the immune system and an interleukine is conservative for both populations, and increases the state of the immune system in such a way that a passive cell always undergoes a transition to a positive state. Explicitly In addition, interleukins are subject to decay in time at a given rate α 4 (u), but, as well known, there are possible mechanisms by which they can be replaced, and they may be produced even by immune system and tumor. Here we are mostly interested in the case when the production term of interleukins is mainly due to an external medical treatment, thus we neglect the dependence of the interleukin production on immune system effectors and tumor cells and assume that the entire production rate is well represented, in the simplest possible way, by an external positive source γ 4 (u) acting on the body. In other words This completes the list of possible processes that are considered significant for the evolution of our four components system. All other interaction parameters appearing in (3) are then equal to zero, as well as the remaining non-interactive operators J 1 and J 3 in (2). Balance equations for our simple kinetic model read then as (10) An expected feature of the model is that integration over u of the third equation leads to the conclusion that the immune system population n 3 is constant in time, since these cells undergo only conservative interactions, and simply change their state without any birth nor death. Another remark may concern interleukins, whose activity variable is actually effective only when the probability distribution ψ 34 (v, w; u) is sensitive enough to variations of its argument w, otherwise all interleukins would behave in the same way when interacting with immune system, and the internal state variable would be useless for them, at least as long as also the dependences of γ 4 and α 4 on u are weak enough. However, for the sake of mathematical generality, we shall keep interleukins depend on their activity variable, though such a dependence can not be "seen" at macroscopic level, as shown below.
Equations (10) belong, apart from the addition of stabilizing damping terms, to a class of kinetic equations for which mathematical well posedness is well established (see for instance [2,20]) on the basis of the theory of approximate solutions in the sense of [27]. However, we are mainly interested here, as a first preliminary approach to the kinetic description of the microscopic process, in the derivation and analysis of reliable macroscopic equations for the observable moments n i . In this respect we notice that equations (10) lend themselves to an integration over the state variable that would single out only macroscopic quantities (moments of the distribution functions), provided all m ij , ν i and α i were taken to be constant. Therefore, we shall stick in the sequel to this strong, but still reasonable simplifying assumption, and leave investigation of other experience-based shapes for those parameters to a future work. However, we realize that one of the cutoffs present in the collision frequencies introduces a further hindrance towards the derivation of a self-consistent set of ODEs for macroscopic densities, namely the appearance in the first equation of a further unknown, the partial density of passive immune cells, say n − 3 , with and n 3 = n + 3 + n − 3 . Since n 3 is determined by the initial conditions, only one between active and passive cells density may be considered as an effective timevarying unknown, and an equation for it is simply obtained by integration of the third equation in (10) over the relevant partial interval. Upon defining we end up with the four dimensional dynamical system in the four non-negative unknowns n 1 , which, in spite of its drastic simplifications, incorporates the expected essential features that may affect tumor evolution in the body. The set (13), along with its own interest, carries a non negligible meaning also at kinetic level, since knowledge of the densities and estimates on the u-dependence of the ψ and ε functions allows, under the present assumptions, also the calculation via (10) of the actual distribution functions. In any case, it is convenient, as usual, to cast the set (13) in dimensionless form, by measuring all densities in units of a typical density, such as n * 2 , and time in units of a characteristic time, for which a proper choice seems to be the inverse tumor proliferation rate in the absence of immune system, τ = (d 12 (m 12 − 1)n * 2 ) −1 . If X 3 = n + 3 /n * 2 , and X i = n i /n * 2 for all other populations, denote dimensionless densities, and the same symbol as before is kept for the new time variable, this procedure yields where all dimensionless parameters are positive, and X = n 3 /n * 2 represents the (constant) size of the overall immune system. The physical meaning of the other seven parameters can be sketched as follows: is the rate at which active immune cells become passive by interaction with tumor; is the proliferation rate of tumor due to interaction with passive immune system; • C =η represents the rate at which interleukins are injected into the system by external sources; is the spontaneous death rate of interleukins; • F = 1 m 12 − 1 measures the destruction rate of host environment due to the action of tumor; represents the spontaneous convergence rate of the host environment towards its saturation (equilibrium) value. Of course, in our scaling, the proliferation rate of tumor by interaction with the host environment is unity. All unknowns are non-negative, and X 3 can not exceed the upper bound X.
As it can be seen from (14), cells of type 2 and 4 are little affected by binary interactions, especially if one considers that host environment is typically much denser than all other populations, and very often is well approximated by a given background in equilibrium [15], with negligible effects of binary encounters on its population. In addition, the last equation, relevant to interleukins, could be solved independently from the others to yield a non autonomous three dimensional dynamical system. For these reasons, we shall investigate in detail in the next sections, analytically as far as possible, two important subcases of the evolution problem (14), in order to emphasize the role played by either of these two auxiliary (but essential) populations in the process, where the actual competing cells are indeed tumor and immune system. Analysis in four dimensions will be hopefully resumed in future work. 4. Qualitative analysis of the reduced model: Tumor-immune systeminterleukins. In this section we shall be concerned with the physical situation in which the host environment has a very prompt and effective reaction to any perturbation of its natural equilibrium state, and is able to re-establish it in an exceedingly small time. This fact can be quantified in a limiting procedure by letting the parameter ν 2 (and then G) tend to ∞, in a sort of zero-order Chapman Enskog expansion, leading, in the language of kinetic theory, to hydrodynamic macroscopic equations in the asymptotic limit [10]. In practice, the second equation in (14) is replaced, to leading order, by X 2 − 1 = 0, and the host environment becomes a sort of huge background, essentially unaffected by the interactive process going on, as conceivable in an initial stage of tumor development. In this partial quasi-steady state approximation, it is convenient to rename variables X 1 , X 3 , and X 4 as Y 1 , Y 2 , and Y 3 , respectively, and to rewrite the set of ODEs as     Ẏ a three-dimensional dynamical system depending on 6 scalar parameters. It can be easily proved that the first octant is a positively invariant set, thus the positivity of solutions, starting from positive initial conditions, is guaranteed. Moreover, the planes Y 1 = 0 and Y 3 = D/E are invariant sets for the trajectories, and then cannot be crossed. The system (15) admits the equilibrium points where the first represents the optimal working conditions of the organism (no tumoral cells and immune system fully active) whereas the second, which makes sense only when it belongs to the phase space, i.e. A ≥ 1/X, represents a scenario of coexistence of tumor and immune system. The local stability properties of equilibrium states E 1 and E 2 can be easily determined by the analysis of the eigenvalues of the Jacobian matrix associated to the system (15). The Jacobian J(E 1 ) is a lower triangular matrix with diagonal elements λ 1 = 1 − AX, λ 2 = −C D E and λ 3 = −E. Therefore, E 1 is locally asymptotically stable if and only if A > 1/X, namely only in presence of the coexistence equilibrium state E 2 . As regards the stability of E 2 , the Jacobian matrix J(E 2 ) is given by with eigenvalue λ 3 = −E < 0 and real eigenvalues λ 1 and λ 2 of opposite sign since the first minor J 33 has trace and determinant both negative in the admissibility domain of E 2 (A > 1/X). Therefore, E 2 is a saddle point when it exists, with a two-dimensional stable and a one-dimensional unstable manifolds, respectively. The phase portrait of (15) for A > 1/X is presented in Fig. 1 (parameter values A = 5, B = 2, C = 1, D = 1.5, E = 1, X = 1/3). Only initial conditions with Y 3 (0) < D/E have been chosen, since an initial level of interleukins above the saturation value D/E is not realistic. The optimal working condition AX > 1, that allows tumor depletion, links a measure of the intensity of the immune system reaction to tumor (AX) to the tumor proliferation rate by interaction with the host environment, which is 1 in the present scaling; thus, such condition quantifies how the reaction of the immune system should be stronger than the proliferation rate of the tumor in order to be able to deplete it. Time evolution of the various components are shown in Fig. 2 for parameter values selected as in Fig. 1. In order to numerically test accuracy of the present steady state approximation, the trend of each component is compared to the one resulting from the numerical solution of system (14), and for this reason all species have been re-labelled according to the notation of Sec. 3. Fig. 2 (left) is relevant to the option G = 5, while in Fig. 2 (right) we take G = 50, with F = 1 for both. Initial conditions are given in all cases by X 1 (0) = 0.115, X 2 (0) = 0.5, X 3 (0) = 0.332, X 4 (0) = 0.02, which determine a point in the 4D phase space belonging to the basin of attraction of the recovery equilibrium state E 1 = 0, 1, X, D E for system (14). Differences are marked only for the "crucial" components X 1 and X 3 ; X 2 would be unity for (15) and X 4 is the same for both systems. It is clear how discrepancies are quite small, and limited to an initial stage of the evolution, even for relatively small values of the parameter G. As expected, the agreement becomes better and better for increasing G, as comparison in Fig. 2 shows, and differences are hardly visible on almost all the time range in the right panel. As a general comment, we observe in Fig. 2 a typical behavior for tumor recovery, in which the fight produces initially a decrease of both tumor and immune system, but the latter prevails in the long run, reverses its trend and gets restored to its equilibrium value, whereas tumor is depleted.
System (15) and its dynamics can be investigated in the framework of the theory of asymptotically autonomous differential systems (see [30,9] and the references therein). Given the differential equationṡ with f, g continuous functions, locally Lipschitz in x, y ∈ R n , respectively, equation (16) is called asymptotically autonomous -with limit equation (17) -if , t → ∞, locally uniformly in x ∈ R n .
The autonomous system (15) can be rewritten as an asymptotically autonomous planar system: by integrating the last equation in (15) we get and substituting it into the second equations in (15) we obtain the following equivalent 2D non-autonomous differential system   having the planar limit system   The limit system has been already investigated in [20]; it admits the equilibria the latter being admissible if and only if A ≥ 1/X, with second coordinate ranging from X (when A = 1/X) to 0 (when A → +∞). The border equilibrium state (0, X) turned out to be a stable node for A > 1/X and a saddle for A < 1/X; the other stationary point is a saddle when it exists. A transcritical bifurcation occurs between the two equilibria when A = 1/X. Moreover, it has been shown in [20] that, for A > 1/X, the basin of attraction of the "optimal" equilibrium (0, X) is given by the domain bounded by the lines Y 1 = 0, Y 2 = 0, Y 2 = X and by the stable manifold of the coexistence saddle point. Such results allow to prove the following theorem Theorem 4.1. Let A > 1/X; there exists a bounded invariant region R in the phase plane [0, +∞) × [0, X] such that every forward solution of (18) starting in R converges towards the equilibrium (0, X) of (19) as t → ∞.
Proof. We will apply Corollary 2.2 of [9] (see Appendix). First, let us recall that equilibrium E 2 of system (15) is a saddle with a two-dimensional stable manifold; such a manifold has a trace γ on the invariant plane Y 3 = D/E given by the stable manifold of the coexistence equilibrium of the limit system (19). Let us consider the closed and bounded two-dimensional region R delimited by the lines Y 1 = 0, Y 2 = 0, Y 2 = X and by the curve obtained by the intersection of the stable manifold of equilibrium E 2 with the plane Y 3 = 0. Such a curve must necessary lie on the left of γ, in accordance with the nullclines (surfaces) of the 3D autonomous system (15) and the resulting sign of the components of its vector field (illustrated in fig. 3, same parameter values as in fig. 1); therefore, the region R is strictly contained in the basin of attraction of equilibrium (0, X) of system (19), and then the first two hypotheses are satisfied, taking as D the interior of R. Moreover, if we choose ρ(Y 1 ,   Figure 3. Nullcline surfaces of system (15) The situation is illustrated in Fig. 4: trajectories of the equivalent non-autonomous system (18) have been compared with those of the limit system (19), with the same initial points; the dotted line represents the intersection on Y 3 = 0 of the tangent plane at E 2 (in the phase space) to its stable manifold, and it is taken as an approximation of the right boundary of the domain R, that is the intersection of the stable manifold of E 2 with Y 3 = 0. Parameter values are A = 5, B = 2, C = 1, D = 2, E = 1, X = 1/3 and it has been chosen Y 30 = 1 for all trajectories. It can be noticed that trajectories of both systems originated inside the region R converge towards E 1 in different ways. In the region between the border of R and the curve γ, trajectories from the same initial point have different destiny, in agreement with the fact that for the non-autonomous system the interleukine population is below its saturation value, which will be reached only asymptotically in time.
For a given initial state, tumor depletion and recovery could be obtained by suitably strengthening the interleukine population. This possibility can be qualitatively examined by choosing an initial point outside of the basin of attraction of the "safety" equilibrium E 1 of system (15) for A > 1/X, which is delimited by the planes Y 3 = 0, Y 3 = D/E, Y 2 = 0, Y 2 = X, Y 1 = 0 and by the two-dimentional stable manifold of E 2 . Then, by varying the parameter D representing the supply rate of interleukins, we try to find a positive threshold value that may lead to tumor depletion even in this case. An example is provided in fig. 5, where we represent the trajectories of the equivalent non-autonomous system (18)  values for A = 5, B = 2, C = 1, Y 30 = 1, X = 1/3, E = 0.5 . When D overcomes the threshold D * 1.43, the trajectories that escaped to infinity for smaller D get reversed and tends asymptotically to the point (0, X), giving tumor depletion. Of course the threshold D * is a function of parameters, and in particular of the initial data. In Table 1 we show the values of D * versus Y 10 , obtained by simulating trajectories starting from (Y 10 , X).  Table 1. Threshold values D * versus initial data Y 10 .
It is worth noticing that, as expected, D * is increasing with Y 10 , and moreover the possibility of interleukins degradation (E = 0) implies greater values of D * with respect to the case considered in [20], in which the same example has been presented in absence of degradation. which the role of a fixed constant background is played by interleukins, that are supposed to have reached, after a short transient, the equilibrium saturation density determined by their supply and decay rates. Again, from a mathematical point of view, this can be justified in an asymptotic procedure in which Γ 4 and α 4 (thus D and E) are of comparable magnitude and large enough. The fourth equation in (14) is replaced by X 4 = D/E, so that density of this population remains constant in the evolution, and for the other three we have the set of ODEs where we have set C * = CD/E, and again we are left with a three-dimensional dynamical system depending on 6 scalar parameters.
It can be easily proved that the first octant turns out to be a positive invariant set, implying positivity of solutions starting from positive initial data; moreover, the plane X 1 = 0 is invariant for trajectories. The system (21) admits, for all positive values of parameters, the equilibrium representing the optimal condition for the organism, with extinction of tumor cells and immune system fully active. Other possible equilibrium points E = (X 1 , X 2 , X 3 ) are characterized by and X 1 positive solutions to the quadratic algebraic equation The discriminant is non negative if and only if A ≥ A * , where A * is a positive quantity depending on the parameters values and always less than 1/X, given by From the Descartes' rule of signs, when A ≥ A * , it follows that: -if C * ≤ G(1 + BX)/(F X) (namely β ≥ 0) then equation (24) has no positive root for A ≤ 1/X (when γ ≥ 0) and 1 positive root for A > 1/X; -if C * > G(1 + BX)/(F X) (namely β < 0) then equation (24) has 2 positive roots for A * ≤ A < 1/X (when γ > 0) and 1 positive root for A > 1/X; the two positive roots coincide when A = A * , and when A = 1/X then a positive root becomes zero. Once a positive root X 1 of eq. (24) is found, it always yields a positive equilibrium state for system (21), which components X 2 and X 3 are given in (23). The above discussion emphasizes the critical value A = 1/X as a transcritical bifurcation value. It is remarkable that system (21) can admit three equilibrium states for suitable parameters values, contrary to the scenario occurring for the interaction between tumor cells, immune system and interleukins investigated in the previous section.
The linear stability of the 'optimal' equilibrium state E 1 is determined by the Jacobian matrix and then E 1 is locally asymptotically stable when A > 1/X, otherwise is unstable, as for the model (15). The equilibrium state is then a nonhyperbolic point for A = 1/X. To determine its local stability for such a critical value and to settle the question of the existence and stability of another equilibrium bifurcated by the nonhyperbolic point, as found above, we will make use of Theorem 4.1 of [8] (summarized in the Appendix), which is based on the use of the center manifold theory [18]. That theorem prescribes the role of the coefficients a and b of the normal form representing the system dynamics on the central manifold, in deciding the direction of the transcritical bifurcation occurring at φ = 0 (see Appendix and the notation defined therein). In particular, if a > 0 and b < 0, then the bifurcation is forward; if a < 0 and b < 0 then the bifurcation is backward (see also [7]).
Proof. We apply Theorem 4.1 of [8] to system (21) to investigate the bifurcation occurring when A = 1/X. Assumption A1 follows from the Jacobian J(E 1 ) given in (25) evaluated at A = 1/X, as discussed above. Let w = (w 1 , w 2 , w 3 ) T be a right eigenvector of J(E 1 )| A=1/X associated to λ 1 = 0; it follows w = 1, having negative components in correspondence of positive components of the equilibrium E 0 , as allowed by the Remark in Appendix. Furthermore, the left eigenvector v = (v 1 , v 2 , v 3 ) satisfying v · w = 1 is given by v = (1, 0, 0) T . The coefficients a and b defined in Theorem 4.1 of [8] can be now explicitly computed; it follows that: where f 1 denotes the first component of the vector field associated to system (21). The coefficient b is always negative so that, according to Theorem 4.1 of [8], it is the sign of the coefficient a which decides the local dynamics around the equilibrium E 1 for A = 1/X. The coefficient a has the same sign of β in the quadratic equation (24) and thus -if C * < G(1 + BX)/(F X) (namely β > 0) then a > 0 and, according to point 3 of Theorem 4.1 of [8], the positive equilibrium (E 2 ) appearing for A > 1/X is unstable and coexists with E 1 which is locally asymptotically stable; then a transcritical bifurcation of forward type occurs at A = 1/X; -if C * > G(1 + BX)/(F X) (namely β < 0) then a < 0 and, according to point 2 of Theorem 4.1 of [8], in a left neighborhood of A = 1/X there exists a positive and locally asymptotically stable equilibrium (E 3 ) coexisting with E 1 which is unstable, and coincides with it when A = 1/X; therefore, in this case a transcritical bifurcation of backward type occurs at A = 1/X.
In presence of a backward bifurcation (namely, when C * > G(1 + BX)/(F X)) the critical value A = A * , where the discriminant of equation (24) vanishes, plays the role of a saddle-node bifurcation value. In fact, the two positive equilibrium states E 2 and E 3 , which are admissible for A ≥ A * , coincide when A = A * ; by some algebra we find that J(E 2 )| A=A * has a simple zero eigenvalue and there are no other eigenvalues on the imaginary axis, thus rank J(E 2 ) | A=A * =2, while rank J(E 2 ) |∂f /∂A | A=A * = 3. Then the stability properties of the equilibrium state E 2 (relevant to the maximum root X 1 of equation (24)) follow from the results in ( [25], p. 253): the bifurcation can be only of saddle-nodes type, since E 3 is stable as proved in Theorem 5.1, and then E 2 is unstable.
The situations about equilibria and their stability are summarized in the bifurcation diagrams reported in Figs. 6 (forward bifurcation) and 7 (backward bifurcation).
The most important feature of this reduced system describing the interactions between tumor, immune system and host environment is the occurrence, for suitable parameter values, of a backward bifurcation [8]; such a bifurcation, which is usually related to epidemic models [19] but with reversed stability properties, can be then found also in this context. However, it is not present in the reduced system investigated in the previous section, describing interactions between tumor cells, immune system and interleukins.
In case of backward bifurcation, it is worth noticing that the system, for proper initial states and parameter values, can evolve towards a scenario characterized by the presence of tumor cells coexisting at equilibrium with immune system and host environment; even if this situation is not the optimal one, it can represent the tumor latency observed in many clinical cases. However, it is remarkable that the level  of the adimensionalized cellular density of the tumor at the stable equilibrium E 3 prescribed by this mathematical model is relatively low (for all values of A) with respect to the corresponding value of the other positive unstable steady state E 2 (see fig. 7). Under such conditions, a locally attractive steady state thus exists even below the threshold 1/X for the crucial parameter A, when the optimal equilibrium (22) is unstable. In order to illustrate the above situation, we show in Fig. 8 the phase portrait relevant to parameter values A = 0.9, B = 1, C * = 6.1, F = 1, G = 1, X = 1 (backward bifurcation). The positive equilibrium state E 3 is locally asymptotically stable, and coexists with the 'optimal' equilibrium E 1 and the positive equilibrium E 2 , which are both unstable. Phase trajectories either diverge or tend to E 3 , and one can easily see the typical saddle behaviour of E 1 and E 2 .
Examples of time evolution of the different components in the safety region attracted by the stable tumor free equilibrium E 1 are given in Fig. 9. Like in the previous Section, comparison is made of the present results to the ones following from the complete four dimensional system (14). Chosen parameter values are related to those of Fig. 1, namely A = 5, B = 2, C = 1, X = 1/3, F = 1, G = 1, with D/E = 1.5, and initial conditions are the same of Fig. 2. The difference is in the parameter E, which takes values 1 in the left panel and 10 in the right panel of Fig. 9. Notice that here the dummy variable X 4 would be equal to D/E for system (21), while the host component X 2 would be different in the two approaches, but the difference can not be seen on the scale of the figure. Comments are essentially the same as in Sec. 4. The present "steady state approximation" seems to be  reasonable already relatively far from the asymptotic limit regime; discrepancies tend to disappear for increasing values of the large parameter, as well as, for fixed parameter values, when time increases. Appendix.
Theorem. ( [9], Corollary 2.2). Let R be a subset of R 2 such that any equilibrium of the limit system (17) in R is the only equilibrium in a sufficiently small neighborhood. Further assume that exist a subset Y of R 2 and an open simply connected subset D of R 2 whit the following properties: • Every bounded forward orbit of the differential system (16) in R has its ω-limit set in Y . • All possible periodic orbits of the limit system (17) in Y and the closures of all possible orbits of (17) that chain equilibria of (17) cyclically in Y are contained in D. • g is continuously differentiable on D and there is a real-valued continuously differentiable function ρ on D such that div(ρg) is either strictly positive almost everywhere on D or strictly negative almost everywhere on D.
Then every bounded forward solution of the limit system (17) in R and every bounded forward solution of the system (16) in R converges towards an equilibrium of the limit system (17) as time tends to infinity.  Theorem. ( [8], Theorem 4.1). Let us consider the system of ODEs with parameter φ dx dt = f (x, φ), f : R n × R → R n and f ∈ C 2 (R n × R); Let x = 0 be an equilibrium of (28). Assume: A1. A = D x f (0, 0) is the linearization matrix of system (28) around the equilibrium x = 0 with φ evaluated at 0. Zero is a simple eigenvalue of A and all other eigenvalues of A have negative real parts; A2. Matrix A has a (nonnegative) right eigenvector w and a left eigenvector v corresponding to zero eigenvalue. Let f k denotes the kth component of f , and Then the local dynamics of system (28) around x = 0 are totally determined by a and b.
(The proof can be found in [8], and in the same paper Table 3 well illustrates these results).
Remark. Taking into account Remark 1 in [8], if the equilibrium of interest in the above theorem is a non negative equilibrium x 0 , then the requirement that w is non negative is not necessary. When some components in w are negative, one can still apply the theorem provided that w(j) > 0 whenever x 0 (j) = 0; instead, if x 0 (j) > 0, then w(j) need not to be positive. Here w(j) and x 0 (j) denote the j-th component of w and x 0 , respectively.