Asymptotic dynamics of some t-periodic one-dimensional model with application to prostate cancer immunotherapy

In the case of some specific cancers, immunotherapy is one of the possible treatments that can be considered. Our study is based on a mathematical model of patient-specific immunotherapy proposed in Kronik et al. (PLoS One 5(12):e15,482, 2010). This model was validated for clinical trials presented in Michael et al. (Clin Cancer Res 11(12):4469–4478, 2005). It consists of seven ordinary differential equations and its asymptotic dynamics can be described by some t-periodic one-dimensional dynamical system. In this paper we propose a generalised version of this t-periodic system and study the dynamics of the proposed model. We show that there are three possible types of the model behaviour: the solution either converges to zero, or diverges to infinity, or it is periodic. Moreover, the periodic solution is unique, and it divides the phase space into two sub-regions. The general results are applied to the PC specific case, which allow to derive conditions guaranteeing successful as well as unsuccessful treatment. The results indicate that a single vaccination is not sufficient to cure the cancer.


Introduction
Over the last few decades, cancer immunotherapy has become a growing area of active research, with a wide variety of approaches developed and clinically implemented in different cancer indications. The rational basis for immunotherapy is the assumption that cancer evolves to evade the control of the host's immune system, and that a proper boost, e.g. by restoration of an impaired function or by countering immunosuppressive mechanisms employed by the cancer, would allow immune system to regain control over the disease (de Visser et al. 2006;Hanahan and Weinberg 2011). Thus, the common denominator of the different modes of cancer immunotherapy is the attempt to manipulate one or several components of patient's own immune system that interact with the disease (in contrast to the chemical or biological therapies that directly target the cancer or cancer-supporting cells). This adds further levels of complexity to the system, making understanding and rationalisation of the treatment effects a non-trivial task.
In order to effectively respond to a developing malignancy, the immune system has to properly recognise the threat, induce massive production of several types of immune cells that perform different complementary tasks, create the suitable microenvironment and successfully deliver the relevant cells and substances to the disease location. Cancer evades the immune control by acquiring the ability to impede these processes at one or more stages, e.g., by reducing the efficacy of antigen recognition, or suppressing the immune response by secreting anti-inflammatory or pro-regulatory cytokines. Many different ways in which cancer may protect itself from the immune system have been discovered (Dunn et al. 2002;Schreiber et al. 2011). When these evasion strategies are understood, the immunotherapy approaches can be designed in order to overcome them by enhancing the impaired immune processes. The variety of suggested immune treatments includes therapeutic vaccines, cytokines, and monoclonal antibodies enhancing tumor-killing mechanisms of the adaptive immune response (Gulley and Drake 2011;Melero et al. 2014). Each treatment type targets a specific component of the complicated network of interactions between the immune system and cancer. Yet, the suggested treatments encounter significant problems when translated into clinic, eventually showing limited efficacy and low response rate, even though sporadic cases of remarkable effectiveness are observed in individual patients. The latter observation raises the possibility that effective design of immune intervention is undermined by our limited understanding of the complicated dynamics of interactions between the two major players, the immune system and the cancer, and how it varies from patient to patient .
Mathematical analysis is a powerful way to gain understanding of behaviour of complex systems, in particular, in biology. Indeed, cancer-immune interactions have been extensively investigated using various mathematical models; see recent reviews in de Pillis et al. (2014) and Eftimie et al. (2011) and references therein. One popular approach is to represent the studied system by a set of differential equations, whereas the time-dependent variables represent the quantities of cells or substances of interest.
The equations' structure reflects the understanding of the most important processes and interactions, whereas the parameters represent the rates or relations that characterise these processes. Such a model can be then investigated analytically or numerically (by simulations) in order to determine the spectrum of plausible behaviours and, if possible, gain insight into optimal design of therapeutic intervention. Hitherto, many different models have been formulated and studied, most of them emphasising the theoretical aspects of the immunotherapy, with limited application to actual experimental and clinical data. Yet, if the model's parameters are retrieved from experiments, or by matching model output versus real clinical data, one can obtain not only a general insight into the systems' dynamics, but rather actually predict the quantitative response of the disease to the treatment. Further, this allows designing and testing various alternative treatment strategies that would be most beneficial, on a population or individual level. Following this line of thought, a model of progression of prostate cancer (PCa) under vaccination treatment was developed and calibrated in Kronik et al. (2010), using data published in the literature, as well as individual profiles of prostate-specific antigen (PSA) in response to vaccination by allogeneic whole-cell vaccination in a phase 2 clinical study (Michael et al. 2005). The model, comprising a system of 7 ordinary differential equations (ODEs), was calibrated individually for each patient, by fitting the individual response data, given as PSA levels measured before, during and after the treatment. It was able to accurately reproduce individual PSA dynamics in response to immunotherapy, for most of the patients. In a subsequent work, it was also demonstrated that individual parameters can be identified at early stages of the treatment, and thus, the model can be used to predict further dynamics of response and eventually optimise the individual treatment schedule, given a welldefined end-point; e.g., stabilisation of PSA levels (Kogan et al. 2012). Thus, such a model may be relevant to a clinical application of immunotherapy, with a potential to dramatically change the approach to treatment design. All the above investigation was carried by numerical simulations, which has the obvious drawback of exhibiting model behaviour only under the specific sets of parameters values, at which the model is numerically solved. Here, we wish to expand this investigation by an analytical study of the asymptotic behaviour of this specific model. This will allow further insight into its predictive ability, scope and limits of applications, which can be important to further utilisation of the model, and the above approach as a whole.
This article is built as follows. Section 2 presents a general result concerning asymptotic dynamics of one ODE with the right hand-side F(t, x) being t-periodic and monotonic in x. Section 3 briefly presents and explains the model to be studied. Section 4 shows that the system dynamics is asymptotically one dimensional and particularly simple when only one boost of immunotherapy is given. Finally, Sect. 5 presents the analysis for the case of periodic impulsive vaccination. In that case giving the treatment periodically, we asymptotically obtain a t-periodic right-hand side of the equation, and the general theorem from Sect. 2 applies.

General theorems
We consider a general Cauchy problem of the forṁ where the function F fulfils the following conditions.
(A1) F is continuous and locally Lipschitz continuous in x for all (t, Remark 1 Notice that due to t-periodicity of F it is enough to consider initial time t 0 ∈ [0, 1). Moreover, Assumptions (A1)-(A4) imply that unique solutions of Eq. (1) exist globally in time (that is for every t ≥ 0) independently of t 0 ∈ [0, 1) and Then any solution of Eq. (1) with x 0 > 0 tends to +∞ as t → +∞.
Proof Let t ∈ [0, 1) and n ∈ N be arbitrary. Integrating Eq. (1) on the interval [t + n, t + n + 1] we obtain ln x(t + n + 1) and therefore monotonicity of F in x implies Under the assumption F A > 0 the sequence (x(t + n)) n∈N satisfies the inequality Theorem 2 Let F fulfil (A1)-(A4) and Moreover, assume that Then there exists a curve γ : , then x(t 0 + 1) = γ (t 0 ) and the curve γ prolonged in a periodic way on [1, +∞) is a periodic solution of Eq. (1).
Proof According to Remark 1 solutions of Eq. (1) are unique, and thus trajectories cannot cross each other. We show that the behaviour of the solution is thus determined by the inequality between values at t = t 0 and t = t 0 + 1, cf. Fig. 1 for which yields For the casex > x(t 0 ) > x(t 0 + 1) analogous arguments leads to the following conclusion It remains to prove that a periodic solution exists. We show that if x 0 is small enough, then x(t) → 0, while if x 0 is large enough, then the solution diverges to +∞. The continuous dependance of solutions on initial data would complete the proof.
To this end, first we give some estimates of the solutions. Due to (A4) there exists r > 0 such that |F(t, x)| < r for all (t, x) ∈ D. Thus, for any t ∈ [t 0 , t 0 + 1) we have and, Now, consider x 0 sufficiently small. Then Similarly, we can show that for x 0 sufficiently large x(t 0 + 1) > x(t 0 ). Clearly, As f (t) = lim x→∞ F(t, x) > 0 and due to continuity and t-periodicity of F the function f is continuous periodic and thus there exists ε > 0 such that f (t) > ε.
The uniform convergence of F(t, ·) to f (t) implies that there existsx such that for all x ≥x the inequality F(t, x) > ε/2 holds. Thus, for x 0 >x we havė

Now, continuous dependance and uniqueness of solutions of Eq.
(1) finishes the proof.
Remark 2 Notice, that if F(t, x) is independent of t, that is F(t, x) = G(x) for some locally Lipschitz continuous function G, then 1. either there existsx > 0 such that G(x) = 0, and then x(t) → 0 for 0 < x 0 <x, while for x 0 >x we have x(t) → +∞ as t → +∞; 2. or G(x) > 0 for x > 0, and then x(t) → +∞ for any x 0 > 0.
In the next section we use our results to analyse the PC immunotherapy model proposed in Kronik et al. (2010).

A model for PCa immunotherapy
In this section we present and briefly explain the model for vaccination treatment of PCa which was proposed in Kronik et al. (2010). This model describes interactions among the tumour antigen, either endogenous (V p ) or injected as a vaccine extracted from cell lines (V ), several types of immune cells and prostate cancer cells, P. The injected vaccine is assumed to be taken up by the naive dendritic cells. It is injected into the dermis where it stimulates maturation of dendritic cells (we assume that there is a large pool of immature dendritic cells) into mature antigen-presenting cells (D m ) at the rate k i ; cf. Banchereau and Palucka (2005) and Banchereau and Steinman (1998). Each dendritic cell takes up an amount of vaccine, n V , during maturation. The mature antigen-presenting dendritic cells migrate from skin into lymph nodes at the rate k m with the probability α l to join the pool of functional antigen-presenting dendritic cells (D C ). The functional dendritic cells become exhausted (Kajino et al. 2007) and turn into regulatory dendritic cells (D R ) (Langenkamp et al. 2000) at the rate k C R, and the latter die with rate μ D . Moreover, functional antigen-presenting dendritic cells recruit and activate tumor-specific cytotoxic T lymphocytes (CTLs; C) (Janeway et al. 2005) that can kill tumor cells, at the rate a C . In parallel, regulatory dendritic cells recruit regulatory lymphocytes, known as Tregs (R) (Hollenbaugh and Dutton 2006), at the rate a R . The Tregs inactivate CTLs (George et al. 2003) with rate proportional to both cell type concentrations, with the coefficient k C R. Both types of lymphocytes, CTLs and Tregs, are also dying with rates μ C and μ R , respectively. The existence of regulatory cells cannot be ignored as increased number of regulatory CD4 + CD25 high T cells were detected in the blood and tumour tissue of early stage PCa patients (Miller et al. 2006). Finally, the tumour cells, P, are assumed to grow exponentially with rate g in absence of CTLs, while the latter kill the tumour cells with rate proportional to both populations, with coefficient a P and additional factor that decrease tumour killing efficacy with increasing tumor burden (cf. Kronik et al. 2008;Kogan et al. 2010;Piotrowska et al. 2013 for additional information about this saturation). All the variables related to the specific cell populations, reflect amounts of the corresponding cells at time t. The above-described cascade of immune reactions that eventually leads to the suppression of cancer cells, is represented by the following system of ordinary differential equationsV = −k l n V V, with an initial condition reflecting one boost of V 0 to the system without immunity, that is (V 0 , 0, 0, 0, 0, 0, P 0 ), V 0 , P 0 > 0. Later, we study addition of multiple boosts, which can be represented by adding appropriate δ-functions (i.e. V i · δ(t − t i ) for each boost V i at time t i ) to the first equation in (2). For the comparison with previous analysis, note, that in Kogan et al. (2012) and Kronik et al. (2010) it was assumed that V p = 0, i.e. the endogenous immune response induced by tumour antigens is either fully suppressed, or negligible. In this case, without treatment, asymptotically, the tumour always will grow without limit, while other system components will stay at zero. All other parameters were estimated (averages and biologically plausible ranges) in Kronik et al. (2010). For facilitation of our analysis, we non-dimensionalise this system, in order to reduce the number of parameters. New parameter values can be easily computed from the formulae below and the estimations reported in Kronik et al. (2010). We define the following variables and parameters: These substitutions transform the system into the following form: The non-dimensional formulation reduces the number of parameters from 14 to 9, without changing the qualitative dynamics. In fact, the original variables values are simple linear functions of those in (3), and can be easily recomputed, once the parameters are given. In the following, we study the behaviour of system (3).

Dynamics of Eqs. (3) after one boost
Studying the dynamics of Eqs.
(3) we can integrate subsequent equations one by one obtaining as k m = k v for the parameter values estimated in Kronik et al. (2010). Therefore, d m → V p k m , and moreover if V p = 0, then d m → 0 exponentially. Next, we can also obtain an explicit formula for d c , and now as k C R = k m in the original parameter estimation, we have not only exponential functions, but also combinations of exponential functions with a linear term t. However, the full expression for d c is not important for the asymptotic behaviour of the solutions. The first five equations of (3) form a linear subsystem, and it is easily seen that d c → V p k m k C R and again d c → 0 for V p = 0, and the convergence is of the order te −k C R t . Similarly, Eventually, for every ε > 0 there existst > 0 such that Let us consider one-dimensional dynamical system governed bẏ According to Remark 2, the dynamics of Eq. (5) depends on the magnitude of c ∞ . Clearly, -If c ∞ is small, such that Eq. (5) has no positive equilibrium, then x → ∞.
-If c ∞ is sufficiently large, then there exists a positive equilibriumx = ac ∞ g − 1 of Eq. (5) and for 0 < x 0 <x the solution tends to 0, while for x 0 >x the solution tends to ∞.
As ε is arbitrary, to have a positive equilibriump one needs This condition is equivalent to Corollary 1 To achieve a cure of the disease one needs g < a μ D μ R k := g max and V P > gk m μ C μ D μ R k C R k(g max − g) .
Corollary 1 means that to achieve the cure after one boost the tumour cannot be highly reproductive and moreover natural influx V P of mature dendritic cells must be sufficiently large. It is obvious that in this case, Moreover, if c ∞ < g a , then there is no positive equilibrium and p also tends to ∞. Notice, that if V p = 0 we have c ∞ = 0, and therefore for any positive initial data p(0) = p 0 we have p(t) → ∞. Kronik et al. (2010) the cure cannot be achieved after one boost.

Impulsive vaccination
Following Corollary 2, when V p = 0 (or when it is very small) we need to apply more vaccination boosts in order to reduce tumour load. In the following we assume V p = 0 and consider the sequence of boosts of the same magnitude, V 0 , given each time interval Δt, starting at t = 0. We obtain impulsive equation for V which can be easily solved on intervals of the form [nΔt, (n + 1)Δt], n ∈ N. Clearly, for t ∈ (nΔt, (n + 1)Δt) we obtain where V (nΔt) + is the level of vaccine just after (n + 1)th boost, as the first boost is given at t = 0. It is obvious that V tends to a periodic function as t → ∞. More precisely, Then the asymptotic dynamics of D m is governed bẏ and as V is periodic, in the limit we also obtain a periodic expression for d m : where we assume Δt = 1, for simplicity (one can always choose appropriate time units and scale all the parameters and variables accordingly). Similarly, we can show that the functions d c (t), d r (t), r (t) and c(t) are asymptotically periodic with the period equal to 1. Clearly, asymptotic equation for any of these variables can be expressed in the following waẏ where F and G are periodic with the period 1, G(t) > 0 (in fact, for all the variables except c the function G is constant, Lemma 1 For any solution x(t) of Eq. (7) there exists a periodic function x * with period 1 such that |x(t) − x * (t)| → 0 as t → ∞.
Proof The solution of Eq. (7) reads Therefore, Next, let us take t ∈ (0, 1) and calculate This proves that x(t) converges to the periodic function Now, we need to studẏ where x = p and H (t) = ac(t) is smooth and periodic. We can show that the dynamics of Eq. (9) depends on an average of H . Let us define H A = 1 0 H (s)ds. Theorem 3 If g > H A , then x(t) → ∞, as t → +∞, for any t 0 , and x 0 > 0. If g < H A , then for any t 0 ∈ [0, 1) there exists x * (t 0 ) > 0 such that Proof It is a simple corollary from Theorems 1 and 2.

Conditions for cure and unsuccessful treatment
Now, we come back to Eq. (6), and consider together the magnitude of one boost V 0 with the frequency of vaccine application. For every ε > 0 there existst such that for any t >t we have is sufficient to cure the disease for p 0 <p min = ac min g − 1. On the other hand, if p 0 >p max = ac max g − 1 or c max < g a , then p(t) → ∞ as t → ∞.
At the end we should notice, that the result of Theorem 3 can be used for better approximation of the conditions for cure, but then the average value of c has to be calculated explicitly, which is possible but tedious task.

Conclusion
The model proposed in Kronik et al. (2010) and analysed in the present work, describes the interaction between advanced PCa and immune system under vaccination treatment. The immune response is governed by a one-way signalling cascade starting from the tumour-specific antigen, which can be endogenous (represented by V p ) or externally administered (represented by V ). The signal proceeds through antigenpresenting cells to the activation of cytotoxic effector cells (c). In parallel a built-in delayed negative feedback (carried out by r ) inhibits results from the same stimulation, and inhibits the immune attack on the cancer. Without treatment, this system stabilises at certain level of effector cells (c = c ∞ ), which contributes to the reduction of the tumour load. Corollary 1 gives the condition for a possibility of tumour control. In brief, these conditions demand that the pro-inflammatory arm of the system will be stronger than pro-regulatory arm and also be strong in comparison to the growth rate of the tumour. When this is the case, the outcome depends on the tumour load at the beginning of the process: if it is higher than a thresholdp, computed in Sect. 4, the disease will progress, while if it is lower thanp, disease will be eliminated.
One boost of vaccination has a transient effect on the system. A signal is sent from V to c, as described above, but eventually it wanes off, due to elimination in all the cell populations along the chain of reaction. Therefore, c temporarily increases above it's steady level, c ∞ , but eventually goes back down to that level. It is reasonable to assume that at the beginning of the treatment the disease was not under control (i.e., p >p, (otherwise the treatment would not be needed). Then, if the conditions of Corollary 1 hold, the boost V 0 should be high enough, so that the transient increase in c will reduce the tumour load p below the thresholdp. If this happens, the immune system may restore the control and eliminate the disease without additional treatments. If on the other hand, such substantial boost is unfeasible, or the parameters make it is impossible (violating assumptions of Corollary 1), a single vaccination cannot be effective against cancer. In particular this applies when, as assumed in Kogan et al. (2012) and Kronik et al. (2010), V p is zero or negligible.
In the latter case of ineffective endogenous antigen response, the hope is that periodic treatment may help. Indeed, as shown in Sect. 5, when the vaccination is applied periodically, the signal sent to the effector cells reaches higher levels, and can be practically evaluated to be higher than some minimal level c > c min > c ∞ . Thus, in principle a prolonged may generate an on-going immune reaction at levels higher than endogenous. As shown by Corollary 3, there is a possibility for elimination of tumour of any initial size, given that the treatment frequency (Δt) and doses (V 0 ) can be controlled. Importantly, using this Corollary, or, if needed, more explicit numerical calculations, the minimal requirements of treatment intensity parameters (i.e., Δt and V 0 ) can be computed, given the values of the model parameters. This opens the path for optimisation of the immune treatment regimens, both on population and individual levels. We hope that our analysis will facilitate and motivate further advancement of mathematical and computational methods in the field of oncoimmunotherapy.