On the global dynamics of a model for tumor immunotherapy

Understanding the dynamics of human hosts and tumors is of critical importance. A mathematical model was developed that explored the immune response to tumors that was used to study a special type of treatment [3]. This treatment approach uses elements of the host to boost its immune response in the hopes that the host can clear the tumor. This model was extensively studied using numerical simulation, however no global analytical results were originally presented. In this work we explore the global dynamics to show under what conditions tumor clearance can be achieved.

1. Introduction.Cancer cells are host cells that begin to proliferate in an uncontrolled and non-specific way that leads to the development of a tumor.Mathematical models of tumor-immune dynamics have added to our understanding of how host immune cells and cancerous cells evolve and interact (see [1] for a review).Most of these models, do not address spatial aspects of tumor dynamics, but focus on ordinary differential equations (ODE) over time.Others have looked specifically at models of tumors as they grow spatially, but these tend not to focus on any aspect of the host's response (c.f.[9]).Still other models have explored aspects of treatment and how that can reduce or eliminate the tumor (c.f.[10]).Two ODE models of host-tumor dynamics looked specifically at the idea of using host immune factors as treatment strategies to augment or abrogate host-tumor dynamics [3] and [11].What is different about this work is that the treatment that was being considered as not drugs, but products derived specifically from the host immune response.This was considered in an effort to boost the host's own immune response.The idea for this type of treatment arose from clinical studies that attempted this in patients, yielding mixed results [5,6,7,8].An additional paper explored optimal control methods for determining the best treatment strategy [12].They designed a control functional to maximize the effector cells and interleukin-2 concentration while minimizing the tumor cells.The two modeling studies examined tumor-host systems using numerical bifurcation techniques and determined local stability criteria; however no global analysis was performed.Here, we extend the work of the 574 DENISE KIRSCHNER AND ALEXEI TSYGVINTSEV first model [3] by performing analysis on the system that allows us to make specific recommendations based on the results.To this end, Let Ẋ = V (X), X ∈ M be a system of n nonlinear ordinary differential equations where M ⊂ R n is a certain open set.
We take an arbitrary differentiable function I : R n → R. Using the vector field V (X) one can calculate the directional derivative K(X) = D V .I(X) =< V (X), gradI(X) >.The surface K(X) = 0 divides M into domains D i such that M = ∪ i D i and the sign of K is constant in each D i .That has a clear dynamical interpretation known as Lyapunov: in the region D i where K > 0 (K < 0) the solutions X(t) of the differential equation intersect the constant level surfaces I(X) = c, c ∈ R in the direction of increasing (decreasing) c.Thus, once the geometry of the level surfaces of I is well understood, some useful information about the behavior of solutions X(t) in D i can be obtained.It is natural to call I the quasi-Lyapunov function.
We apply this method in Section 2 of our paper where the functions I are chosen in a particular way to guarantee that its directional derivatives have a simple algebraic form.Our analytical results predict global dynamics for this nonlinear ODE system and allow us to predict under what conditions the tumor can be cleared.
2. The model.To begin we present the mathematical model that is studied.A full explanation is given in [3], but we summarize briefly as follows.Tumor cells are tracked as a continuous variable as they are large and generally homogeneous; they are defined as y(t).Immune cells are also large in number and are those cells that have been stimulated and are ready to respond to the foreign matter (known as antigen) (known as effector cells); they are define as x(t).Finally, effector molecules are represented as concentrations z(t).These are self-stimulating proteins for effector cells which produce them.The equations that describe the interactions of these three state variables are given in [3]: In equation (1a), the first term represents stimulation by the tumor to generate effector immune cells.The parameter c is known as the "antigenicity" or the strength of this characteristic.Since tumor cells begin as self, c represents how different the tumor cells are from self cells (i.e., how foreign).The second term in (1a) represents natural death and the third is the proliferative enhancement effect of IL-2.s 1 represents a treatment term where by a physician administers effector cells that have been taken from a patient, stimulated to a large degree, and then subsequently infused back into the patient.In equation (1b), the first term is a logistic growth term for tumor growth, and the second is a clearance term by the immune effector cells.In the final equation, (1c), IL-2 is produced by effector cells (in a Michaelis-Menton fashion, i.e. dose response) and decays via a known half-life (third term).The second term, s 2 , is a treatment term that represents administration of IL-2 (manufactured) by a physician to a patient, to again stimulate effector cell growth and proliferation.
To help with interpretation of the mathematical results, we present a Now, we move to perform analysis on this system.Previously, it was shown in [3] that the system of 3 ordinary differential equations (1) possesses a number of fixed points whose stable (unstable) properties depend on parameter values.The fixed points with y = 0 play a key role since those states are "tumor free" states.The stability conditions of these points are of medical importance and was the main object of the study performed in [3].Nevertheless, the Lyapunov linear stability is rather a local phenomena: it guarantees the convergence of neighborhood trajectories to a given fixed point only for sufficiently close initial conditions.We consider the paper [3] as an important step in investigation of the dynamics of (1) in the vicinity of its fixed points.The present paper contributes to understanding of the global dynamics of (1) and is based on quasi-Lyapunov functions techniques.
3. Description of quasi-Lyapunov functions I(y, z), J(y, z) and its constant levels.In this section we are going to study some aspects of global dynamics of (1).Being non-linear, these equations do not have obvious solution and one has to apply a qualitative approach.We perform this analysis with the help of a quasi-Lypunov function method described briefly in Introduction.
First, we write the equations (1b) and (1c) in the form where V 1,2 are 2-dimensional vector fields given by: (3) Proposition 1.The function is a real analytic first integral of the vector field V 2 defined for any y > 0, z > 0.
The function is a real analytic first integral of the vector field V 1 defined in the region Proof.The systems of differential equations ( ẏ, ż) = V 1,2 can be easy solved by quadratures since the equations for ẏ in both cases depend only on y.The expressions for the first integrals I and J can be derived straightforwardly from these solutions.
The goal is to study the dynamics of system (1) using the evolution of constant levels of I, J under the flow.The restriction to D of constant levels I(y, z) = const is a family of curves starting from the vertical line y = b −1 and ending at the zaxis.They intersect both axes transversally.One can check that ∂I/∂y > 0 for y ∈ (0, b −1 ), hence all curves I(y, z) = const are strictly decreasing ones.The constant levels J(y, z) = const, restricted to D, all start from the point (b −1 , 0) approaching the z-axis asymptotically so that z → ±∞ as y → +0.For our purpose it is sufficient to represent the constant levels of I, J in D , shown schematically in Fig. 1.
The constant levels of J(y, z) and I(y, z) together with directions of increasing of its constant levels Remark 1.The elementary analysis of (1) shows that the following domains (7) are invariant under the flow.In this paper, the only solutions of (1) we study are those starting from U 1 .
The directional derivatives of I(y, z) and J(y, z) along V i , i = 1, 2 are given by the following expressions We note that according to Proposition 1 the following identities hold: D V1 J = D V2 I = 0. Thus, the derivative ( 8) is the full one with respect to the vector field, V , of (1) i.e.
It follows from ( 8) and ( 9), since x > 0, that the zero-velocity conditions The previous equation defines a certain curve Γ ⊂ D. As seen from ( 11

DENISE KIRSCHNER AND ALEXEI TSYGVINTSEV
The curve Γ separates D into two regions D 1 and D 2 so that and Let us assume that the solution curve (x(t), y(t), z(t)) of ( 1) intersects the cylindric surface defined by i.e. the projection X yz = (y(t), z(t)) intersects the curve Γ in yz-plane.Let < , > denotes the scalar product in R 2 .We introduce the normal vector N to Γ pointed from D 1 to D 2 .It can be easily calculated with help of ( 11): N = grad(z − L(y)).
The expression W =< N, V 1 > +x < N, V 2 >, evaluated at the point x > 0, (y, z) ∈ Γ is positive or negative depending on whether the trajectory enters or leaves the cylindric surface C described above.The simple calculation shows with The algebraic equation P (y) = 0 can have N = 0, 1, 2, 3 different roots y i in the interval Y = (0, b −1 ) depending on parameter values.That defines a division of Y into the union of disjoint intervals Y = i Y i such that the sign of P considered in Y i is constant.
As follows from the previous discussion, at the moments of intersection of X yz (t) with Γ, the corresponding direction of the trajectory (from D 1 to D 2 or inverse) can be calculated based on division Y = i Y i and values of x and y.For example, the following two conditions indicate that (x(t), y(t)) transverses Γ from (19) Let X(t) = (x(t), y(t), z(t)) be an arbitrary solution of (1) starting from the point P = (x(0), y(0), z(0)).We denote X yz (t) = (y(t), z(t)), P yz = (y(0), z(0)) their projections on yz-plane.The analysis given above allows to describe possible types of behavior of X(t).Theorem 3.1.For every solution of (1) starting from U 1 there exists T > 0 and z max > 0 such that z(t) < z max for all t > T i.e. the z-component will always reach a regime of bounded behaviour.
For every solution of (1) starting below the surface z = z min = s 2 /µ 3 we have two following possibilities A) ∀ t > 0, z(t) < z min and lim t→+∞ (y(t), z(t)) = (0, s 2 /µ 3 ) B) There exists T > 0 such that z(t) > z min for all t > T i.e. starting from a certain moment z(t) stays all the time above the plane z = z min .
Proof.That follows from the analysis of evolution of constant levels I = const under the action of flow of system (1) (see Fig. 2).There are several cases to consider.Case I. Let P yz ∈ D 1 .In this region both functions I(Y (t)), J(Y (t)) increase, so the curve X yz (t), t > 0 behaves in D 1 in such a way it transverses the levels of I, J in the direction of its increasing as shown in Fig. 2. At a certain moment it meets the curve Γ, on which the derivatives İ(Y (t)), J(Y (t)) change sign simultaniously.After, once X yz belongs to D 2 , it starts to move in the direction of decreasing of I, J until it meets Γ again etc.We can have oscillations with X yz jumping from one side of Γ to another and never going far from it.We note that some examples of such oscillatory behaviour were established in [3].Case II.If P yz ∈ D 2 then X yz (t) goes in the direction of decreasing of I, J until it meets Γ and enters D 1 .The value z max is defined by intersection of the curve Z : I(y, z) = γ = const tangent to Γ with the z-axis.We note, that z max depends only on parameters of the system and not on the particular solution.The explicit formula for z max requires solving of a quartic algebraic equation.Nevertheless one easy estimates it with help of (12).Indeed, the line I(y, z) = c * = I(b −1 , M ) will lie above Z.Therefore, it is sufficient to solve I(0, z) = c * with respect to z.That gives the result Case III.It is possible that X yz (t), starting from a certain moment, will stay in D 1 or D 2 and tends asymptotically to Γ.
We note that the above conclusions can be derived based uniquely on the analysis of constant levels of I and one could avoid the use of J. Nevertheless we believe that J is worth describing as it will contribute to future work.

Dynamics including x(t).
The quasi-Lyapunov function K(x, y).In this section we will proceed to further analyze of (1) with help of a new quasi-Lyapunov function constracted in another way.The idea is to introduce the variable x in the study which was missed in quasi-Lypunov functions I, J from the previous section which depended on y and z only.We will use a new quasi-Lyapunov function K(x, y) which will be a key tool in the proof of the following theorem.Theorem 4.1.Let parameters of the system (1) satisfy the inequalities or Then every solution with initial conditions of the form x 0 > 0, 0 < y 0 < b −1 , z 0 > 0 will have a tumor free limit behavior lim t→+∞ y(t) = 0.
In the case the limit lim t→+∞ y(t) = 0 holds for solutions whose initial conditions satisfy Proof.We write the system (1) in the following separated form where U 1,2 are 3-dimensional vector fields given by The idea is to look for a solution of D U1 K = 0 functionally independent from J(y, z).
The system of 3 differential equations ( ẋ, ẏ, ż) = U 1 can be integrated by quadratures since the equation for ẏ ż depend only on y and z respectively Using this solution it is easy to derive the corresponding first integral K K(x, y) = br 2 x + (c + s 1 b) log(1 − by) − bs 1 log y .
(27) One sees that K is analytic in D since the logarithmic functions it contains depend on 1 − by and y which are positive in this domain.
As in the previous section, we calculate the derivative where E(y, z) = P (y)z + g 1 Q(y) , (29) with P, Q obtained by collecting the terms containing the factor z.
We do not provide explicit formulas for P, Q (which are large) but give the following relations defining them: It is important to study the evolution of K(x, y) under the action of flow of the vectot field U 2 .Here one follows the same ideas presented in Section 3 where the evolution of level lines of I and J was analysed.The levels K(x, y) = const in the plane xy are shown on Fig. 3.
We are interested how the curve L = {E(y, z) = 0} crosses the domain D. Solving E(y, z) = 0 with respect to z one obtains where Depending on parameters, different scenarios of crossing L∩D can happen.Below we indicate some interesting cases and give their dynamical interpretations.
As easily follows from the analysis of levels of K (see Fig. 3) together with (1b), once the inequalities (21) or (22) hold, every solution of (1) starting from U 1 possesses the following property lim t→+∞ y(t) = 0 , i.e. it converges to the tumor free state.At the same time, the previous analysis does not allow us to prescribe the limit behavior of x(t) or z(t) as t → +∞.Now, using again (31) and (32) one checks that if (33) is satisfied together with then D U2 K > 0 will hold in U 2 .After all substitutions, the conditions (33) and (37) can be written as (23).
The theorem is proved.

Conclusion and biological implications.
Here we study a host-tumor model with two different types of treatments: one that enhances the type of effector immune cells that can fight a tumor, and the other is enhancment of effector molecules that stimulates effector immune cells to proliferate.In a previous numerical study of this model [3], they showed that in the system with combined treatments (i.e. both s 1 > 0 and s 2 > 0) that there are regions in phase space where the tumor can be cleared.Here, we study the same system rigorously, and find, in Theorem 3, that under specific conditions we can define exactly what conditions allow for tumor clearance.Sections 3 and 4 of the paper are devoted to the global dynamics of (1) with help of quasi-Lyapunov functions.Interestingly, the inequalities (21), ( 22), (23) depend on c, the antigenicity of the tumor, and thus there are still likely regions where c is so small (or even 0) that the tumor will not be cleared.This means that for tumors that are very similar to self (i.e.c very small) there will not be enough of a signal to the effector cells to clear the tumor.Future work could consider what conditions could be derived for cases where tumors are minimally antigenic and how the immune response could be boosted to be effective in those cases.We would like to comment on the connection between studies performed in Section 3 and Section 4. Indeed, they are closely related.In Section 3, with help of quasi-Lyapunov functions I(y, z) and J(y, z) we could study the dynamics in the projection yz-space only i.e. the information including x variable was missing.In Section 4 this gap was compensated for introducing the quasi-Lyapunov function K(x, y) allowing the capture of the time behaviour of x variable.So, together, I, J, K functions constitute a full set giving information about the dynamics in xyz-space.We believe more can be derived from the study of these quasi-Lyapunov functions and there may be others.That will be reported in our future work.

O D 1 D 2 Γ:
levels of I(y, z) : levels of J(y, z) : the curve Γ : the solution projection
), the only real solutions of L(y) = s 2 /µ 3 are y = 0 and y = b −1 .So, the graph Γ of z = L(y) lies above the line z = s 2 µ 3 for y ∈ (0, b −1 ) with L(y) reaching a maximum M > s 2 /µ 3 in this interval since L(0) = L(b −1 ).To determine explicitly the critical points of L one has to solve a cubic equation.The elementary calculation provides the following upper bound on M M = max y∈(0,b −1 )

Figure 3 .
Figure 3.The level lines K = const with direction of increasing. z

Table 1 .
Table of Parameters for ease of parameter interpretation: Parameter Values.