Abstract
Background
Antitumor therapies aim at reducing to zero the number of tumor cells in a host within their end or, at least, aim at leaving the patient with a sufficiently small number of tumor cells so that the residual tumor can be eradicated by the immune system. Besides severe sideeffects, a key problem of such therapies is finding a suitable scheduling of their administration to the patients. In this paper we study the effect of varying therapyrelated parameters on the final outcome of the interplay between a tumor and the immune system.
Results
This work generalizes our previous study on hybrid models of such an interplay where interleukins are modeled as a continuous variable, and the tumor and the immune system as a discretestate continuoustime stochastic process. The hybrid model we use is obtained by modifying the corresponding deterministic model, originally proposed by Kirschner and Panetta. We consider Adoptive Cellular Immunotherapies and Interleukinbased therapies, as well as their combination. By asymptotic and transitory analyses of the corresponding deterministic model we find conditions guaranteeing tumor eradication, and we tune the parameters of the hybrid model accordingly. We then perform stochastic simulations of the hybrid model under various therapeutic settings: constant, piecewise constant or impulsive infusion and daily or weekly delivery schedules.
Conclusions
Results suggest that, in some cases, the delivery schedule may deeply impact on the therapyinduced tumor eradication time. Indeed, our model suggests that Interleukinbased therapies may not be effective for every patient, and that the piecewise constant is the most effective delivery to stimulate the immuneresponse. For Adoptive Cellular Immunotherapies a metronomic delivery seems more effective, as it happens for other antiangiogenesis therapies and chemotherapies, and the impulsive delivery seems more effective than the piecewise constant. The expected synergistic effects have been observed when the therapies are combined.
Introduction
A key problem of antitumor therapies is finding a suitable scheduling of their administration to the patients. Of course a major problem in medical oncology is avoiding severe therapyrelated side effects which, unfortunately, may cause the death of the patient. However, also in the ideal case of no sideeffects, a therapy aims at reducing to zero the number of tumor cells in the host, within its end. Indeed, also if a single tumor cell remains the patient has a tumor. Actually, the requirement might theoretically be milder by accepting to leave the patient with a sufficiently small number of tumor cells so that the residual tumor can be eradicated by the immune system. In any case, both the duration and the scheduling of a therapy becomes of great relevance, as experimentally shown and theoretically studied [1]. In this paper we shall focus on a computational study of some kinds of immunotherapies, whose underlying key idea is to modify the natural interplay between tumor cells and immune system, by boosting the latter.
Tumor cells are characterized by a vast number of genetic and epigenetic events eventually leading to the appearance of specific tumor antigens, called neoantigens. Such antigens trigger antitumor actions by the immune system [2], thus resulting in the tumorimmune system interaction taking place. These observations provided a theoretical basis to the hypothesis of immune surveillance, i.e. the immune system may act to control and, in some case, to eliminate tumors [3]. Only recently studies in molecular oncology and epidemiology accumulated evidences of this [4]. The competitive interaction between tumor cells and the immune system is extremely complex. As such, a neoplasm may very often escape from immune control. A dynamic equilibrium may also be established with the tumor surviving in a microscopic undetectable "dormant" steady state [5]. If this is the case, on the one hand a dormant tumor may induce metastases, on the other hand over a long period of time, a significant fraction of the average host life span, the neoplasm may develop strategies to circumvent the action of the immune system, thus restarting to grow [2,4,6,7]. We stress that this evolutionary adaptation, termed "immunoediting" [4], may negatively impact on the effectiveness of immunotherapies, as shown in [8].
Regarding immunotherapies, although the basic idea of immunotherapy is simple and promising [9], the clinical results are controversial since a huge intersubjects variability is observed [1012]. Immunotherapies are divided into two broad classes: passive and active therapies [13]. Among the passive ones, the Adoptive Cellular Immunotherapy (ACI) consisting in injecting cultured activated immune effectors in the diseased host [14,15] is probably the most important. Active immunotherapies aim at stimulating the immune response by expanding, for instance, the proliferation of cytotoxic T cells. Among these, a prominent role is played by Interleukinbased therapies [16,17].
Regarding the mathematical modeling of tumorimmune system interactions and and related therapies, many papers have appeared using various approaches. For instance, ordinary differential equations (ODEs) are used in [5,8,1316,1828], the theory of kinetic active particles is used by Bellomo and Forni in [29,30] and hybrid agentbased models have been introduced by Motta and Lollini [31,32]. In [14] Kirschner and Panetta proposed a largely influential ODEbased model of TumorImmune system (TIS) interplay, whose variables are tumor cells, effector cells and the concentration of interleukins2 (IL2). This model is able to explain various kinds of experimentally observed tumor size oscillations [3338] as well as both macroscopic and microscopic constant equilibria. Although a vast array of behaviors is mimicked by the solutions of the KirschnerPanetta (KP) model, the tumorfree equilibrium is unstable for all biologically meaningful values of the parameters. However, in [39] we have shown that resetting the model in a hybrid setting where the interleukins are modeled with a continuous variable and the tumor and the immune system are modeled with discretestate continuoustime stochastic process, the eradication via immune surveillance can be correctly reproduced. Since eradication is a fundamental topic in the study of immunotherapies, here we extend our hybrid version of the KP model to investigate the effects of both interleukinbased therapies and ACIs. Although our hybrid version of KP model is highly idealized, we think that it can provide useful information on the design of the above mentioned therapies.
Methods
In the next sections we recall the KP model [14], its hybrid definition [39] and we extend the hybrid model with general immunotherapies.
The deterministic KirschnerPanetta model
In [14] the following model of the dynamics of tumorimmune system interaction was proposed
where T_{* }(t), E_{* }(t) and I (t) denote, respectively, the densities of tumor cells, effectors of the immune system and interleukins. The tumor induces the recruitment of the effectors at a linear rate cT thus c may be seen as a measure of the immunogenicity of the tumor. In other words, according to [14]c is "a measure of how different the tumor is from self". The proliferation of effectors is stimulated by the interleukins. The average lifespan of effectors is and the average degradation time for interleukin is . The source of interleukin is modeled as linearly depending on effectors, and it also depends on the tumor burden. Finally, continuous infusion immunotherapy may be delivered when effectors and interleukins are injected at constant rates σ_{E}, σ_{I }≥ 0.
In the case of no therapy, i.e. σ_{E}, = σ_{I }= 0, the main results obtained in [14] are that (i) the tumorfree equilibrium point (0, 0, 0) is always unstable, (ii) it exists positive c_{m }≪ 1 such that for c ∈ (0, c_{m}) there is only one is locally stable equilibrium, whose size is very large due to the low value of c, (iii) it exists a c_{M }> 0 such that if c ∈ (c_{m}, c_{M }) there is a unique periodic solution whose period and amplitude decrease if c increases and, finally, (iv) when c >c_{M }there is a unique globally stable equilibrium, whose size is a decreasing function of c. Thus, we note that this model explicitly precludes the possibility of tumor suppression in absence of immunotherapies. If idealized infinitely long constant continuous infusion therapies are considered the behavior of the system is complex, but in all the possible meaningful combinations of the parameters it is possible to find regions where globally stable limit cycles exist, as well as regions where there is cancer suppression. Moreover, there is a threshold value such that, for higher values of σ_{I}, there is an unbounded growth of effectors, leading to severe sideeffects.
A hybrid model with constant therapies
As discussed in [39], the lowlevel oscillations predicted by model (1) make stochastic effects on the cell populations worth investigating. Unfortunately, it is possible to see that a purely stochastic model with discrete populations becomes computationally too hard to analyze, being the number of IL2 huge. In this case a hybrid approach, despite being more costly than the deterministic counterpart, still permits a feasible analysis.
We now recall the hybrid model of TIS interplay that we defined in [39] and that extends the deterministic model (1). Variables T and E of the hybrid model are obtained from the densities T_{* }and E_{* }in (1) converted into total number of cells by means of the volume V (e.g. the blood and bone marrow volumes for leukemia). We have T_{* }= TV ^{1 }and E_{* }= EV ^{1}. This leads to the ODE system
Note that the modified deterministic model (2) is obtained by the original KirschenerPanetta deterministic model (1) by means of a nonsingular linear transformation. As it is very wellknown linear transformations of the state variables do not change the topological properties of the solutions, and as a consequence these transformations do not change all the stability properties of equilibria [40,41].
From (2) a bidimensional stochastic process ruling the dynamics of T (t) and E (t) is linked to a scalar differential equation ruling the dynamics of I (t). In [39] it is discussed that the dynamics of I (t) can be assumed to be approximated by a linear ODE with randomly varying coefficients, which are constant in the intervals between two consecutive stochastic events. We briefly discuss how this model is simulated, for a detailed description of the underlying algorithm we refer to [39]. The algorithm adopted is an extension of the Gillespie Stochastic Simulation Algorithm (SSA) [42,43]. The SSA simulates a trajectory of the continuoustime discretstate Markov process underlying the system. To use this algorithm we write the stochastic events as reactions modeling birth and death, which for this model are described in Table 1. The set of all the events is denoted as . At each step the SSA solves equations to discover the putative time for the next reaction to fire, and probabilistically decides which reaction fires. Intuitively, given to each reaction a propensity function a_{i }(x) is associated. With the system state x at time t the value a_{i }(x) dt gives the probability of the next reaction to fire in the infinitesimal time [t, t + dt). All the propensity functions in are timeindependent, meaning that they depend on a state which is constant in between two stochastic events. Differently, R_{4 }depends on the continuous part of the hybrid model, i.e. a_{4 }(t) depends on I (t), and is hence timedependent.
Table 1. Birth and death reactions for the hybrid model
We remark that the original SSA is assumed to simulate only timeindependent chemical reactions. The algorithm used to simulate this model takes inspiration from algorithms simulating hybrid systems with timeindependent propensity functions [44,45] and algorithms simulating purely stochastic systems with timedependent propensity functions [46,47]. In [39] the original SSA equations (i.e. see [43]) for the putative time for the next reaction are rephrased in this setting. In particular, with the last stochastic event fired at time t_{n }the putative time τ for the next reaction to fire is determined by solving
with χ a random number with distribution Exp(1). The solution of equation (3) has no analytical form, thus requiring iterative methods to find its solution.
We remark that the same equation in the case of only timeindependent propensity functions yields the well known SSA strategy to generate exponential jumps. As far as model analysis is concerned, differently from the deterministic model (1), the stochastic simulations performed in [39] show that, at least in some cases, suppression of the neoplasm might be reached, thanks to the conjunction of the intrinsic tendency of the TumorImmune System to oscillate with the stochastic dynamics.
A hybrid model with general therapies
In this paper we consider more general immunotherapies than the constant ones considered in [14,39]. Thus, here we allow that the therapy related influxes σ_{E }and σ_{I }are functions of time, and in next sections we shall focus on two common periodic scheduling of a therapy. In the deterministic setting, model (2) can trivially be modified by adding timedependent immunotherapies, that is
It is important to notice that in the realistic case of finite duration therapies the deterministic system always predicts tumor regrowth, being the tumorfree equilibrium (0, 0, 0) unstable. Practically, if at the end of the therapy the solution is very close to the tumorfree equilibrium, immediately after the end the tumor restarts growing. We remark that in the oncological context it is important the state in which the tumor is at the end of a therapy, e.g. the tumor shrinkage up to an undetectable size, but it is far more important what is observed during the followup visits, i.e. that the tumor has not regrown. In this framework the deterministic model is unable to reproduce the reality and it structurally gives negative answers on the effectiveness of the therapy. These reasons, combined with the possible lowlevel oscillations of model (1) makes again stochastic effects worth investigating. Along the line of [39], this shall permit to have an estimation of the average value of the random therapyinduced eradication time t_{er }better than the one obtainable in the deterministic framework, that is better than min{t  T (t) <1} which would be actually very rough. As in [39], switching from the deterministic to the hybrid version of the model is furthermore justified by the fact that, at the best of our knowledge, there is no proof that the average of the state variables of a nonlinear birth and death stochastic processes follows the dynamics of the corresponding ODE system, when variables assume low values. In fact, the ODE system would be a good approximation of the purely stochastic model if the state variables were "sufficiently large". In any case, by using the above deterministic approximation of t_{er }its standard deviation could never be computed.
We now present an SSAbased algorithm to simulate the hybrid version of the above model. This new algorithm extends the algorithm which simulates model (2) in a natural way. The reactions in Table 1 are left unchanged provided that the propensity function for R_{7 }is modified as a_{7}(t) = V σ_{E }(t) to reflect the timedependency of σ_{E }(t). As a consequence in this case we have two timedependent propensity functions a_{4 }(t) and a_{7 }(t), hence we define . Let us model the eventinduced changes in the values of T and E by the vectors
Here the ith components of ν_{T }and ν_{E}, i.e. ν_{T,i }and ν_{E,i}, describe how the reaction R_{i }affects the population T and E, respectively. With the last stochastic event fired at time t_{n }we define
and
The exact simulation algorithm for the hybrid system with generic immunotherapies is Algorithm 1 and is defined in Table 2. With the current state x and the propensity functions in evaluated, i.e. the putative time for the next stochastic event follows
where χ is a random number with distribution Exp(1). Notice that equation (5) is a natural extension of equation (3) and, as such, it does not admit a general analytical solution. Consequently numerical methods to find its solution are again required. In the following, when looking for a solution of equation f(x) = 0 in [a, b] with f(a)f(b) <0 and f continuous we always adopt the bisection method with double stopping criteria b  a <10^{8 }∧ f(x) <10^{6}. Once the jump is determined, given , the next event to fire find is R_{j }if
where r is a random number U[0, 1] and a_{i}(t + τ ) = a_{i}(x) if . It is important to notice that equations (5) and (6) would reduce to the standard equations used in the SSA if the system was entirely stochastic.
Some considerations about the ODEs constituting the model (4) are worth discussing. If the last stochastic event happened at time t_{n }and the next will happen at time t_{n+1 }the equation for I' reads as
which, since in such intervals the other two state variables are constant, is a linear ODE with constant input and constant coefficients. Thus, given I(t_{n}) = I_{n }its analytical solution in (t_{n}, t_{n+1}) is
where
Notice that equation (8) is necessary to evaluate, at each step of Algorithm 1, the value of a_{4}(t). In the next sections we will consider specific types of immunotherapies and, according to their mathematical modeling, we may discuss on equation (8).
Table 2. The Hybrid Simulation Algorithm (Algorithm 1)
Values of the parameters
We discuss now the values of the parameters used to simulate the model. In Table 3 we recall the parameters used in [39] to perform the simulations. Those values and ranges are originally given in [14,48]. Note that those values pertain to mice and they were taken from [20,49], where accurate fitting of real data concerning laboratory animals was performed. As in [39] volume V is estimated to be 3.2 ml. This follows by the body weight of a femal chimeric mouse ranging from 20 to 40 grams, and by considering also that their blood volume ranges from 5.8 to 8 ml per 100 grams. It follows that V reasonable ranges from 1.16 up to 3.2 ml. We remark that we will give the parameters for the functions σ_{I}(t) and σ_{E}(t) in the next sections when the therapies will be discussed.
Table 3. Values of the parameters
Interleukinbased immunotherapies
In the next paragraphs we discuss how to specialize model (4) with either a piecewise constant or an impulsive interleukinbased immunotherapy. In both cases, we consider a therapy starting at time t_{s}, ending at time t_{e }and consisting of the injection of molecules of IL2. As already said a really uninterrupted continuous infusion therapy in [t_{s}, t_{e}], such as those in [14,39], is not realistic. Here we consider a continuous infusion therapy delivered at preset times
where θ_{i }∈ [t_{s}, t_{e}] for i = 0, . . . , k. In the following, with a slight abuse of terminology, we refer to each θ_{i }∈ Θ as a therapy session.
Piecewise constant therapy
In the case of a continuous infusion each therapy session has a duration of A timeunits and a constant influx rate d_{i }in the ith infusion. We model the therapy as
where A <min{θ_{i+1 } θ_{i } i = 0, . . . , k  1} (i.e. no overlap between two therapy sessions), and g(t, A) is the "window function" of amplitude A, i.e.
where H(·) the wellknown Heavyside function. Although this scenario can be simulated by Algorithm 1, some considerations about the equations that need to be solved are worth discussing. In particular, when solving for t >t_{n }equation (8) by using equation (10) the key point is to evaluate the integration part of Equation (8), i.e.
where t = t_{n }+ φ and η (x) = σ_{I}(x)e^{μIx}. Since Θ = {θ_{0}, θ_{1}, . . . , θ_{k}}, we start by considering the set of therapy sessions which start and complete in the time window (t_{n}, t_{n }+ φ), such a set is
Now, since Θ contains ordered timepoints then also θ(t_{n}, φ) does. We want to get the index of the minimum and the maximum θ_{j}'s from θ(t_{n}, φ); we define
so that the indexes min and max denote such values. To solve R(t_{n}, t_{n }+ φ) we split the integral in three time intervals [t_{n}, θ_{min}), [θ_{min}, θ_{max }+ A] and (θ_{max }+ A, t_{n }+ φ] so that
The integral in [θ_{min}, θ_{max }+ A] is the summation of θ(t_{n}, φ) nonzero basic integrals over the subintervals [θ_{j}, θ_{j }+ A] for any θ_{j }∈ θ(t_{n}, φ), that is
Notice that for any θ_{j }∈ θ(t_{n}, φ) it holds that σ_{I}(x) = 0 for any x ∈ (θ_{j }+ A, θ_{j+1}) (i.e. when the therapy is not delivered) yielding . Furthermore, since
then the overall integral evaluates as
Notice that this quantity can be easily computed in a iterative fashion. The cases of the rightmost and leftmost integration intervals are similarly accounted. Let us consider the leftmost interval [t_{n}, θ_{min}), we consider whether t_{n }is included in in the (min  1)th therapysession. We have this definition
if θ_{min1 }< t_{n }≤ θ_{min1 }+ A and 0 otherwise. This holds since in the uppermost case . Similarly, in the interval (θ_{max}, t_{n }+ φ] by cases on the relation between t_{n }+ φ and θ_{max+1}
if θ_{max+1 }≤ t_{n }+ φ < θ_{max+1 }+ A and 0 otherwise. This holds since in the uppermost case . We remark that all these combinations of cases are necessary because of all the possible combinations of the parameters t and t_{n }with the set Θ.
Impulsive therapy
In many cases the infusions are very short implying that the in flux rate reaches very large values, so that one may approximate σ_{I}(t) as a train of pulses, i.e.
Here, u_{i }is the ith injected dose of molecules of IL2 and δ(t  θ_{i}) is the Dirac's delta function centered at t = θ_{i}. By simple algebraic manipulations it is possible to see that here function R(t_{n}, t) is given by
where . Finally, we stress that an alternative to the use of the Dirac's generalized functions is to represent the drug deliveries as impulsive kicks given to I(t), thus becoming an impulsive differential equation [50,51]. Indeed, we have that at time θ_{i }∈ Θ the value of I(θ_{i}) increases by u_{i }units thus yielding
for i = 0, . . . , k.
Adoptive cellular immunotherapies
In the next paragraphs we discuss how to specialize model (4) with either a piecewise constant or an impulsive ACI. As in the previous section we consider a set of k therapy sessions Θ = {θ_{i } i = 0, . . . , k}.
Piecewise constant therapy
As before we assume
This scenario can be simulated by Algorithm 1 where, in this case, the function A_{7}(τ) is given by
Impulsive therapy
We consider the case of an impulsive ACI where at each θ_{i }∈ Θ there is the rapid infusion of w_{i }cultured effectors. Thus, we could proceed similarly to the case of impulsive ILbased therapy. However, here the introduction of generalized functions does not lead to simplifications, and as a consequence we shall model the therapy as
Thus to the "natural" stochastic events external deterministic events are superimposed. As a consequence, at such times (i) the differential equation ruling the dynamics of I(t) must be updated and (ii) all the propensity functions involving E change value. This means that the integral terms of equation (5) in Algorithm 1 should be modified to consider such changes. However, results in [52] allow us for a modification of the algorithm that we discuss now. One of the fundamental properties of the SSA is that, with the system at time t, for any possible value of τ the propensity functions are constant in the time interval [t, t + τ). The same holds in some configurations of the hybrid model studied in [39]. However, this does not hold in the case of impulsive ACI since, if t_{n }is the time of the nth stochastic event and θ_{j }>t_{n }is the closest injection after t_{n}, the propensity functions a_{3 }and a_{5 }change after θ_{j}, invalidating the property. Of course, the timedependent propensity functions a_{4 }and a_{7 }are, by definition, potentially nonconstant. We argue that there is a similarity between schedulingbased SSAs for systems with delays and this hybrid system (i.e. consider this ACI as a set of events scheduled at the preset times in Θ). Although this differs from such algorithms where the scheduling times are stochastically chosen during the simulation, this allows to modify Algorithm 1 into Algorithm 2 presented in Table 4, which is indeed inspired by the algorithms discussed in [5357]. As for such algorithms, the correctness of Algorithm 2 relies on the existence of a general schema of SSAbased algorithm with piecewise constant timedependent propensity functions [52] where is stated tha, and this algorithm respects that schema.
Table 4. The Hybrid Simulation Algorithm (Algorithm 2)
Combining IL2 therapies and ACI
Combining therapies requires combining the results of the previous sections. To shorten the presentation we briefly discuss how to perform the simulations of the hybrid model with combined immunotherapies.
Whenever an impulsive ACI is considered, independently of the IL2 therapy, the model can be simulated by Algorithm 2. In the other cases the model can be simulated by Algorithm 1.
Results
In the next sections we perform both asymptotic and transitory analysis of the solutions of deterministic system (4). We show the existence of deterministic conditions that guarantee the eradication of the tumor, and that will be used to tune the parameters of our hybrid model when performing its simulations, which are discussed in the forthcoming section.
Deterministic asymptotic analysis
With the aid of elementary dynamical systems theory [41] and by using some mathematical properties summarized in the additional file 1, here we briefly investigate how the therapies may influence the asymptotic behavior of the solutions of the deterministic system (4). We shall focus on the mathematically idealized case of infinite length of the therapy and on the deterministic conditions guaranteeing the eradication of the tumor.
Additional file 1. Supplementary Materials (text). Results about scalar differential inequalities, scalar linear ODEs with periodic coefficients, impulsive ODEs and combined impulsive immunotherapies are given.
Format: PDF Size: 200KB Download file
This file can be viewed with: Adobe Acrobat Reader
IL2 immunotherapy
We start from the case of the delivery of a ILbased monotherapy (i.e. σ_{E}(t) = 0). By setting T = E = 0 it follows I(t) = J_{∞}(t), where J_{∞}(t) is the asymptotic solution of
Thus, we found a tumorfree state
which is, however, unstable. In fact, this follows by setting
and linearizing since the equation for tumor cells reads as . Even though this means that the tumorfree state (0, 0, J_{∞}(t)) is unstable, this does not mean that under an ILbased monotherapy a tumor cannot be eradicated. In fact, as we show, there exist conditions implying that (T, E, I) → (0, +∞, +∞), namely the tumour is eradicated. This can be verified by using basic differential inequalities recalled in the additional file 1. Indeed, from the differential inequality
it follows that I(t) >J(t) and, for large times, I(t) ≥ J_{∞}(t). In turn, the inequality
yields the inequality
Finally, from the properties of periodical linear differential equations recalled the additional file 1 we can observe that
implies that (T(t), E(t)) → (0, +∞). In other words, the deterministic model predicts that eradication is possible only at the price of killing the patient for the excess of stimulation of the immune system, that is E(t) → +∞. We remark that the above inference does not consider finite length therapies. In fact, if at the end of the therapy the tumor has been eradicated, the number of effectors and the interleukin density will both decay. Inequality (17) is the condition which we obtain by the deterministic model. We now refine such a condition to account for the three therapies that we mentioned so far: (i) constant, (ii) piecewise constant and (iii) impulsive. It holds that:
(i) when a constant therapy where σ_{I}(t) = σ_{I }is considered J_{∞ }= σ_{I }/μ_{I }and condition (17) is
(ii) if we consider the piecewise constant σ_{I}(t) of equation (10) where the therapy sessions are periodically scheduled with period P, duration A and the rate of each session is the same, that is ∀θ_{i }∈ Θ. d_{i }= d, we have that
where
In this case condition (17) reads as
where
and γ_{1 }<1.
(iii) if we consider the impulsive σ_{I}(t) of equation (11) where all the therapy sessions are periodically scheduled with period P and the rate of each session is ∀θ_{i }∈ Θ.u_{i }= U we have
where t mod P is the standard modulo operation. Then the eradication condition (17) is
ACI
When only ACI is delivered (i.e. σ_{I}(t) = 0) there is the following tumorfree asymptotic solution
where ε_{E}(t) is the asymptotic solution of
which is the equation for E when T = I = 0. In the case of periodically delivered therapy with period P, the linearized equation for tumor cells is
and hence the local eradication condition (17) is
However, by averaging both the sides of (20) yields
which implies that
Finally, in this case the local stability condition corresponding to equation (17) becomes
As for the case of IL2 monotherapy we now focus on the therapies considered so far, we have that:
(i) if we consider constant σ_{E}(t) then condition (22) is σ_{E }>rg_{T }V/a.
(ii) if we consider a piecewise constant σ_{E}(t) with infinite length t_{e }= +∞, period P where θ_{i }= iP, duration A and injection rate b_{i }= b for any θ_{i }∈ Θ, then similarly to the case of the ILbased monotherapy one can show that
where
if 0 ≤ t < A and
if 0 ≤ t < P and the eradication condition can be consequently computed.
(iii) if we consider the impulsive therapy σ_{E}(t) where the common injection rate is w_{i }= w for i = 1, . . . , +∞, by using the results reported in the additional file 1 we have that
and hence the eradication condition (22) is
Combined therapies
Finally, we shortly consider when both the therapies are delivered. In this case there is the following tumorfree asymptotic solution
where ε_{∞}(t) is the asymptotic solution of
In the case of synchronous delivery with common period P the local eradication condition becomes
As for the monothrapeutic case in some of the scenarios that we mentioned it is possible to infer analytical local eradication conditions. In the additional file 1 we show the derivation of the eradication condition for combined impulsive therapies with synchronous delivery.
Global stability of the eradication
We conclude by investing the global stability of the eradication. Since for sufficiently large t it is I(t) ≥ J_{∞}(t) then from the differential inequality
follows that, for large times, E(t) ≥ ε_{∞}(t). Then by
it follows that T(t) < X(t), where
Finally, it follows that if ∀x(0) > 0. X(t) → 0 then also T(t) → 0, and the tumor eradication is globally asymptotically stable [41].
Deterministic transitory analysis
Results of the previous section refer to the highly idealized case of a infinite horizon therapy. However, real therapies have a finite duration and, more important, the host organism has a finite lifespan. Thus, in this as well as in other applications of computational biology and medicine it is natural to wonder whether such results can be used at all [8,28]. This is critically related to the velocity at which the solutions of the equations studied in the previous section tend to their asymptotic solutions.
As far as the IL2 monotherapy is concerned, the velocity of growth of E(t)  which in turn determines the velocity of reduction of T(t)  is ruled by the difference p_{E}〈J_{∞}(t)/(g_{E }+ J_{∞}(t))〉  μ_{E}. Moreover, independently of the initial conditions the function J(t) converges to J_{∞}(t) in some multiple of average degradation time of the interleukin, i.e. 1/μ_{I}, which is small. Thus, this means that very soon the asymptotic solution is reached. Observe now that since J_{∞}(t)/(g_{E }+ J_{∞}(t)) <1 it follows that if p_{E }< μ_{E }then the constraint (17) is never fulfilled. We stress that this is the case for the values listed in Table 3. In practice, since in general 〈J_{∞}(t)/(g_{E }+ J_{∞}(t))〉 <<1 and since μ_{E }is small, it follows that p_{E }must be far larger than μ_{E}. This requires that, for some of the settings that we will simulate in the next sections, the value of p_{E }could be different from the one given in Table 3. Biologically, this might substantially reduce the number of patients to whom the IL2 monotherapy might be effective.
Further discussions are worth. In the case of impulsive therapy unless IL2 is injected every few hours J_{∞}(t) → 0 rapidly, so that the eradication is unlikely unless large doses are delivered. This is also mirrored in the corresponding local eradication condition (19). Furthermore, in equation (16) it is also very important parameter g_{E}, whose value used in our simulations is very large. Thus, if μ_{I}P and g_{E }are large we may roughly say that  unless huge doses are delivered or p_{E }is particularly large  the coefficient of E in inequality (16) is almost always comparable to μ_{E}, so that a large rate of injection is required to fulfill eradication condition (19).
In the case of piecewise continuous delivery of IL2 J_{∞}(t) takes few hours to get sufficiently closer its maximum plateau value, as given by equation (18). This suggests that for this kind of drug delivery the duration of each therapy session, i.e. A, should be a quite large fraction of the unity. This, of course, poses some practical problems since the patients should receive very long daily infusions. However, in some recent clinical trials on cyrcadian rythmstuned delivery of chemotherapy some special 24hours infusors have been experimented [58]. Roughly speaking, the above fact might be related to the "indirect effect" of ILbased therapies. In fact, they aim at triggering the expansion of the number of effector cells which, in turn, kill tumor cells.
Differently, as far as the ACI monotherapy is concerned, the velocity of convergence of ε_{E}(t) is some multiple of the average lifespan of the effector, i.e. 1/μ_{E}. Such a value is generally quite big and, for instance, it is 33 days about in our simulations. This means that the convergence is very slow and that, unless the duration of the therapy is exceptionally long, the results of the asymptotic analysis cannot be used as a basis for the stochastic simulations. Similar considerations can be done for the combined therapy.
Finally, it is important to recall that the conditions that derived in the previous section are of local nature. In order to guarantee the eradication for generic nonsmall initial conditions (T(0), E(0), I(0)) the constraint specific to the simulated therapy has to be largely fulfilled.
Stochastic simulations
We performed stochastic simulations of the model under various therapeutic settings, whose results are now reported. We have mostly considered a singlemonth daily therapy, i.e. Θ = {1, . . . , 30}. When different schedules are considered the parameters are explicitly reported. In all the figures representing simulation days and number of cells are given on the xaxis and the yaxis, respectively. To perform simulations a JAVA implementation of model (4) and Algorithms 1 and 2 has been developed.
IL2 immunotherapy
In Figure 1 a single stochastic run of the hybrid model with only IL2 based immunotherapies is shown. We simulated both piecewise constant (left) and impulsive (right) IL2 daily immunotherapies with initial configuration T(0) = 10^{5 }and E(0) = I(0) = 0. In the former case each therapy session lasts A = 0.2 days, i.e. 4.8 hours, and d_{i }= 4 · 10^{7}/A for i = 1, . . . , 30. In the latter case at each therapy session the value u_{i }for the injection is u_{i }= 4·10^{7 }for i = 1, . . . , 30. Thus the total injected drug per day is the same in both cases. Moreover, according to the transitory analysis, in both cases we twentyfold increased the value of p_{E }reported in Table 3, i.e. we used p_{E }= 20 · 0.1245. As we already said, this requirement shall reduce the number of patients to whom this therapy might be effective. All the other parameters are as in Table 3.
Figure 1. Single run, IL2 monotherapy. Singlerun of piecewise constant (left) and impulsive (right) IL2 daily immunotherapy. In (left) A = 0.2 and d_{i }= 4 · 10^{7}/A for i = 1, . . . , 30. In (right) u_{i }= 4 · 10^{7 }for i = 1, . . . , 30. In both cases T(0) = 10^{5}, E(0) = I(0) = 0 and p_{E }= 20 · 0.1245. All the other parameters are as in Table 3.
Notice that (i.e. compare with the transitory analysis) the piecewise constant immunotherapy seems more efficient than the impulsive one. Indeed (i) in the piecewise case the eradication of the tumor is reached at t_{e }≈ 20 days and the maximum tumor size is around 10^{6 }= 10 · T(0). Differently, (ii) in the impulsive case the tumor is eradicated a few days later, i.e. t_{e }≈ 23 days, and the maximum tumor size is almost 20·T(0). At the eradication day the number of effector cells is around 4 · 10^{6 }in both cases, whereas the density of IL2 is of the order of 10^{7 }in (left) and 10^{4 }in (right).
Adoptive Cellular Immunotherapy
In Figure 2 a single stochastic run of the hybrid model with only ACI is shown. As for ILbased therapy, we simulated both piecewise constant (left) and impulsive (right) daily ACIs with initial configuration T(0) = 10^{5 }and E(0) = I(0) = 0. In the former case each therapy session lasts A = 0.2 days and b_{i }= 25 · 10^{4 }for i = 1, . . . , 30. In the latter case at each therapy session 5·10^{4 }effector cells are injected, i.e. w_{i }= 5 · 10^{4 }for i = 1, . . . , 30. Thus the number of injected effectors is equal in both cases. All the other parameters are as in Table 3.
Figure 2. Single run, daily ACI monotherapy. Singlerun of piecewise constant (left) and impulsive (right) daily ACI. In (left) A = 0.2 and b_{i }= 5 · 10^{4}/(V A) for i = 1, . . . , 30. In (right) w_{i }= 5 · 10^{4 }for i = 1, . . . , 30. In both cases T(0) = 10^{5 }and E(0) = I(0) = 0. All the other parameters are as in Table 3.
Note that the figures show no remarkable difference in the tumor response. In particular, in both the simulations the eradication is obtained at around day 15. In both cases, at the eradication day the number of effector cells is around 6 · 10^{5}, and the density of IL2 is of the order of 10^{2}.
Finally, to discover the relation between the frequency of the therapy sessions and the dosage of each session, in Figure 3 a single stochastic run of the hybrid model with impulsive weekly ACI is shown. In those simulations we used an impulsive ACI with a weekly schedule (i.e. θ_{i }= 7i + 1 for i = 0, . . . , 3) with dosage w_{i }= 35 · 10^{4 }for i = 1, . . . , 30. Note that this means that once a week a number of effectors is injected equivalent to the number of effectors given in an entire week of Figure 2 (right). The other parameters are as in Table 3.
Figure 3. Single run, weekly ACI monotherapy. Singlerun of piecewise constant (left) and impulsive (right) weekly ACI, i.e. in both panels θ_{i }= 7i + 1 for i = 0, . . . , 3. In (left) A = 1.4 days and b_{i }= 35 · 10^{4}/V A i = 0, . . . , 3. In (right) w_{i }= 35 · 10^{4 }for i = 1, . . . , 3. All the other parameters are as in Table 3.
Notice that it seems that the immune response is slightly better stimulated with this therapy setting than the one in Figure 2 (right). In fact, in this case the eradication day is around 12. The number of effector cells and the density of IL2 are similar to those in Figure 2 (right).
Combined therapies
In Figure 4 a single stochastic run of the hybrid model with combined impulsive IL2 and ACI daily immunotherapies is shown. In (left) both the therapies are given at the same day (i.e. with the same Θ). In (right) the therapies are asynchronous with a shift of 0.5 days, i.e. and for i = 1, . . . , 30. Again, the initial configuration is T(0) = 10^{5}, E(0) = I(0) = 0. The parameters for the therapies dosage and duration are the same for both cases. As far as IL2 therapy is concerned, at each therapy session the value u_{i }for the injection is u_{i }= 10^{6 }for i = 1, . . . , 30. Differently, in each ACI session 10^{4 }effector cells are injected, i.e. w_{i }= 10^{4 }for i = 1, . . . , 30. Notice that both the values are smaller than those used in the corresponding monotherapy scenarios of Figures 1 (right) and 2 (right). Moreover, as in the case of single IL2 therapies the value of p_{E }is set to p_{E }= 20 · 0.1245. All the other parameters are as in Table 3.
Figure 4. Single run, combined immunotherapies. Singlerun of synchronous (left) and asynchronous (right) combined impulsive IL2 and ACI daily immunotherapies. The asynchronous delivery is a shift of 0.5 days. In (left) u_{i }= 10^{6 }for i = 1, . . . , 30 and in (right) w_{i }= 10^{4 }for i = 1, . . . , 30. In both cases T(0) = 10^{5}, E(0) = I(0) = 0 and p_{E }= 20 · 0.1245. All the other parameters are as in Table 3.
As expected, this combined therapy eradicates even though the parameters are lower than those used in the scenarios where the single therapies are used (i.e. Figure 1 (right) and 2 (right)). In both cases the eradication is observed a few days after the end of the therapy. In the case of synchronous therapy t_{er }≈ 39, in the asynchronous case t_{er }≈ 37. In both cases at the day of eradication the number of effector cells is around 2 · 10^{5 }and the density of IL2 is around 10. In both therapies the maximum size of the tumor is almost equal, i.e. it is around 2.5·10^{5}. Finally, since in both cases the proliferation of effector cells is almost equal, it seems that no remarkable differences are observed with these therapy schedules.
Larger initial tumor and ACI
We analyzed the effect of varying the initial number of tumor cells in a scenario with impulsive ACI. Again, we analyzed this scenario because it was the one which permitted a computationally easier analysis. In Figure 5 a single stochastic run of the hybrid model with T(0) = 10^{6 }(left) and T(0) = 10^{7 }(right) is shown. In both cases E(0) = I(0) = 0. Notice that in this case T(0) is either 10 or 100 times larger than the value used in Figure 2. In left panel of Figure 5 at each therapy session 5 · 10^{4 }effector cells are injected, i.e. w_{i }= 5 · 10^{4 }for i = 1, . . . , 30, and the eradication is reached at ≈ 24.5 days. On the contrary, in the case where T(0) = 10^{7 }and the same schedule is applied the eradication was not observed and the tumor size reached, quite rapidly, a size of the order of 10^{9}. In order to obtain eradication also for T(0) = 10^{7}, as shown in the right panel of Figure 5, we increased the number of effectors injected at each therapy session, i.e. w_{i }= 20 · 10^{4 }for i = 1, . . . , 30 (i.e. each injection is 4 times bigger than the one shown in the left panel). In both cases all the other parameters are as in Table 3.
Figure 5. Single run, ACI monotherapy and larger tumor. Singlerun of impulsive daily ACI with T(0) = 10^{6 }(left) and T(0) = 10^{7 }(right). In both cases E(0) = I(0) = 0. In (left) w_{i }= 5 · 10^{4 }for i = 1, . . . , 30. In (right) w_{i }= 20 · 10^{5 }for i = 1, . . . , 30. All the other parameters are as in Table 3.
In both cases the eradication is observed close to the end of the therapy (i.e. day 25) with the number of effector cells being around 10^{6 }(left) and 10^{8 }(right), and the density of IL2 of the order of 10^{2 }(left) and 10^{7 }(right). Notice that in (left) the line of the effectors is interrupted at the eradication time since the simulation is interrupted when T = 0. Moreover, Differently form the other figures here the scales on the yaxes are different since the maximum in (right) is 100 times bigger than the maximum in (left).
Probabilistic analysis of IL2 therapy
As we already said some of the simulations we performed are timeconsuming, especially when the value of T or E become huge. This happens, for instance, in the scenarios of Figure 6 where the simulation time spans from 1 to 30 hours. However, singlerun of stochastic simulations are not much informative and, when possible, any conclusion should be supported by a big number of averaged simulations. To this extent, we performed multiple runs of both the piecewise constant and the impulsive IL2 immunotherapies of Figure 1. This was possible since their running time is of the order of some minutes.
Figure 6. Probability density function, IL2 monotherapy. Empirical evaluation of ϱ(t_{er}) for different piecewise constant (left) and impulsive (right) daily IL2 immunotherapies. In (left) A = 0.2 and d_{i }= 4 · 10^{7}/A for i = 1, . . . , 30. In (right) u_{i }= 4 · 10^{7 }for i = 1, . . . , 30. In both cases T(0) = 10^{5}, E(0) = I(0) = 0 and p_{E }= 20 · 0.1245. All the other parameters are as in Table 3. The densities are obtained by performing 10^{2 }simulations for each configuration.
We defined the following timedependent property over a single simulation: we want to evaluate
meaning that t_{er }is the eradication time for tumor cells before 70 days, more than twice the duration of the simulated therapies. We evaluated the empirical probability density function of t_{er}, denoted ϱ(t_{er}), by performing 10^{2 }simulations for both the scenarios in Figure 1. In all simulations we used the same configuration used in such a figure, that is in (left) A = 0.2 and d_{i }= 4 · 10^{7}/A for i = 1, . . . , 30 whereas in (right) u_{i }= 4 · 10^{7 }for i = 1, . . . , 30. In both cases T(0) = 10^{5}, E(0) = I(0) = 0 and p_{E }= 20 · 0.1245. All the other parameters are as in Table 3.
In case of daily delivery of the IL2 monotherapy, Figure 6 shows the evaluation of ϱ(t_{er}) for piecewise constant (left) and impulsive therapies (right). It is remarkable that for all the simulations the eradication is always reached (i.e. 10^{2 }times out of 10^{2 }simulations).
In Table 5 the average of t_{er}, denoted 〈t_{er}〉 and the standard deviation σ of t_{er }are evaluated for the piecewise constant (left) and impulsive (right) IL2 immunotherapies of Figure 6. By means of the Wilcoxon statistical test, the nonparametric equivalent of the Ttest, we compared the observed realizations of t_{er}. We obtained that the differences in the values of Table 5 are statistically meaningful since p <2.2 · 10^{16 }for the Wilcoxon statistical test.
Table 5. Averages and standard deviation, daily IL2 monotherapy
Probabilistic analysis of ACI
As for the case of IL2 monotherapy the results on ACIs in Figure 2 motivated us to to investigat the relationship between impulsive and piecewise constant ACIs, as well as the influence of the period P between two consecutive therapy sessions. Again, we considered the same timedependent property over a single simulation, i.e. t_{er }= min{t ≤ 70  T (t) = 0}. Also in this case 70 days is more than twice the duration of the therapies that we simulated. We evaluated the empirical probability density function ϱ(t_{er}) by performing 10^{2 }simulations for each of a set of parameter configurations. In all simulations we used the initial configuration T(0) = 10^{5 }and E(0) = I(0) = 0.
Figure 7 shows the evaluation of ϱ(t_{er}) for a daily ACIs with piecewise constant (left) and impulsive (right) infusions. We studied the effect of varying the dosage of the therapy on t_{er}. Given w_{* }∈ {5, 2.5, 2, 1.5, 1}, in (left) each therapy session lasts A = 0.2 days and b_{i }= w_{* }· 10^{4}/(V A) for i = 1, . . . , 30 and in (right) w_{i }= w_{* }· 10^{4 }for i = 1, . . . , 30. The other parameters are in Table 3.
Figure 7. Probability density function, daily ACI monotherapy. Empirical evaluation of ϱ(t_{er}) for different piecewise constant (left) and impulsive (right) daily ACIs. In (left) b_{i }= w_{* }· 10^{4}/(V A) for i = 1, . . . , 30 and in (right) w_{i }= w_{* }· 10^{4 }for i = 1, . . . , 30 where w_{* }∈ {5, 2.5, 2, 1.5, 1}. The value of w_{* }is given in the figure, all the other parameters are as in Table 3. The densities are obtained by performing 10^{2 }simulations for each value of w_{*}.
For any simulation (i.e. 10^{2 }times out of 10^{2 }simulations) the eradication was found if w_{* }≠ 1. In the case of w_{* }= 1 only half of the simulations predicted eradication before day 70. However, at around day 70 the tumor was small (i.e. always less than 100 cells) meaning that the eradication could have been reached immediately in the days after the 70th. Interestingly, in both cases it seems that [1.5; 2] is a range for w_{* }to have eradication before the end of the therapy, as often desired. In Table 6 the average of t_{er}, i.e. 〈t_{er}〉, and its standard deviation σ are evaluated for the piecewise constant (left) and impulsive (right) ACIs of Figure 7.
Table 6. Averages and standard deviation, daily ACI monotherapy
In case of weekly delivered therapy, we considered the simulated therapies where the quantity of effectors are injected each week is equal to the one injected per week in the daily therapies above described. This implies θ_{k }= 1+7 _{* }k, k = 0, . . . , 4 and b_{k }= w_{* }· 10^{4}/(V A) for piecewise constant therapy (left) and w_{k }= w_{* }· 10^{4 }(right). The densities are plotted in Figure 8, and the corresponding mean and standard deviation are shown in Table 7.
Figure 8. Probability density function, weekly ACI monotherapy. Empirical evaluation of ϱ(t_{er}) for different piecewise constant (left) and impulsive (right) weekly ACIs. In (left) b_{i }= w_{* }· 7 · 10^{4}/(V A) for i = 1, . . . , 30 and in (right) w_{i }= w_{* }· 7 · 10^{4 }for i = 1, . . . , 5 where w_{* }∈ {5, 2.5, 2, 1.5, 1}. The value of w_{* }is given in the figure, all the other parameters are as in Table 3. The densities are obtained by performing 10^{2 }simulations for each value of w_{*}.
Table 7. Averages and standard deviation, weekly ACI monotherapy
By means of the Wilcoxon statistical test we compared the observed realizations of t_{er}, grouped by kind of delivery (i.e. impulsive or piecewise constant) and by frequency (i.e. daily or weekly). We obtained that (i) if delivered daily the eradication times of the piecewise constant and the impulsive ACIs are not statistically different (i.e. p > 0.05 for all doses). Also, (ii) if delivered daily the eradication time of the impulsive therapy is significantly smaller than the one with piecewise constant therapy (i.e. p <10^{12}) and (iii) for impulsive ACI the eradication time of the weekly delivered therapy is significantly smaller than the one of the daily delivered therapy (i.e. p <10^{12}). Finally, (iv) for impulsive ACI the eradication time in the weekly delivery is not statistically different from the one of the daily delivered therapy (i.e. p > 0.05).
Conclusions
In this work we extended our hybrid model [39] with ILbased immunotherapies and Adoptive Cellular Immunotherapies (ACIs), both modeled as piecewise constant or impulsive functions. We performed analytical analysis of the corresponding deterministic model, inspired by earlier work by Panetta and Kirschner [14]. We analyzed our hybrid model via stochastic simulations which seem to suggest results of some interest, which we briefly summarize:
(i) by the transitory analysis it turns out that ILbased immunotherapies require very large values of the parameter p_{E}, which might substantially reduce the number of patients to whom it may be used as monotherapy;
(ii) in ILbased immunotherapies the piecewise constant delivery seems more effective for tumor eradication than the impulsive one although at the price of very long infusion sessions;
(iii) in a daily delivered ACI the piecewise constant delivery seems more or less equivalent to the impulsive one;
(iv) in a ACI the impulsive delivery seems slightly more effective than the daily delivery: less frequent deliveries of larger doses ensure a slightly more rapid eradication than frequent deliveries of smaller doses. Note that the latter type of delivery is called metronomic delivery, and it is of great relevance for other antitumor therapies such as antiangiogenesis therapies and chemotherapies [1,59,60]. Furthermore, for those therapies the metronomic delivery is often more effective;
(v) in a ACI the weekly impulsive delivery seems slightly more effective than the weekly piecewise constant delivery;
(vi) when combined impulsive therapies are considered both the synchronous and the asynchronous delivery seem to be effective and no remarkable differences are observable.
Other more predictable effects were observed such as the synergistic effects of combined therapies, or the dependence of the eradication on the initial values. Of course, these results are strongly linked to the specific model, to its ability in describing the dynamics of real tumors and to the chosen parameters.
As far as the model is concerned, we have previously stressed that maybe the hypothesis that the linear antigenic effect cT due to the tumor size should be corrected by assuming a saturating stimulation cT/(1+dT); here we also add that the assumption that E' linearly depend on E could be corrected, as there are cases where this dependence might be nonlinear (see [26] and references therein). Note also that, although computationally useful, representing the piecewise constant delivery of ACI by means of a continuous input σ_{E}(t) is only an approximation. Indeed, in reality the infusion should be more realistically represented as a series of injections of a group of cells each Δt ≪ 1 time units. The time interval Δt should be modeled as a Poisson random variable.
As far as the parameters are concerned, in order to obtain more general biological inferences an extensive and systematic exploration of the space of parameters is mandatory. Of course this will require the exploitation of intelligent algorithms (e.g. approximated stochastic simulations [61,62]) to tackle the computational hardness of model analysis.
Finally, here we have only explored the effects of the intrinsic stochasticity on the dynamics of tumorimmune system interplay under therapy. However, it has been shown that without therapy the extrinsic stochasticity may play a significant role in shaping tumor evasion from the immune control [28]. Moreover, it has also been proposed that realistic bounded stochastic fluctuations affecting chemotherapy may deeply influence the outcome of chemotherapies of solid vascularized tumors [63].
Note that the inclusion of realistic extrinsic noise would require minor changes in the proposed hybrid simulation algorithms besides the inclusion of the stochastic nonlinear equations for correlated bounded noises [28,63]. However, that would require extensive numerical simulations (e.g. a higher number of samples of the stochastic process underlying the hybrid system) when inferring heuristic probability densities of eradication times, for instance.
List of abbreviations used
ACI: Adoptive Cellular Immunotherapy. IL or IL2: Interleukin2. ODE: Ordinary Differential Equation. TIS: TumorImmune system (interplay). KP: KirschnerPanetta (model). SSA: Stochastic Simulation Algorithm.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
AdO conceived of the study and defined the model, GC implemented the model and performed the simulations, AdO, GC and RB analysed the model and wrote the manuscript. All authors read and approved the final version of the manuscript.
Acknowledgements
This work was conducted within the framework of the official agreement on "Computational Medicine" between the European Institute of Oncology and the University of Pisa. The research of AdO has been done in the framework of the Integrated Project "Pmedicine  from data sharing and integration via VPH models to personalized medicine" (project ID: 270089), partially funded by the European Commission under the 7th framework program. The authors would like to thank the reviewers for their comments that helped to improve the manuscript.
This article has been published as part of BMC Bioinformatics Volume 13 Supplement 4, 2012: Italian Society of Bioinformatics (BITS): Annual Meeting 2011. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/13/S4.
References

d'Onofrio A, Gandolfi A, Rocca A: The dynamics of tumorvasculature interaction suggests lowdose, timedense antiangiogenic schedulings.
Cell Prolif 2009, 42:317329. PubMed Abstract  Publisher Full Text

Pardoll D: Does the immune system see tumors as foreign or self?
Annu Rev Immunol 2003, 21:807839. PubMed Abstract  Publisher Full Text

Dunn GP, Old LJ, Schreiber RD: The three Es of cancer immunoediting.
Annu Rev Immunol 2004, 22:329360. PubMed Abstract  Publisher Full Text

d'Onofrio A: Tumor evasion from immune system control: strategies of a MISS to become a MASS.
Chaos Solitons and Fractals 2007, 31:261268. Publisher Full Text

Whiteside TL: Tumorinduced death of immune cells: its mechanisms and consequences.
Semin Cancer Biol 2002, 12:4350. PubMed Abstract  Publisher Full Text

Vicari AP, Caux C, Trinchieri G: Tumor escape from immune surveillance through dendritic cell inactivation.
Semin Cancer Biol 2002, 12:3342. PubMed Abstract  Publisher Full Text

d'Onofrio A, Ciancio A: Simple biophysical model of tumor evasion from immune system control.
Phys Rev E Stat Nonlin Soft Matter Phys 2011, 84:031910. PubMed Abstract  Publisher Full Text

De Vito VT Jr, Hellman J, Rosenberg SA: Cancer: Principles and Practice of Oncology. 9th North American Edition edition. Philadelphia: JP Lippincott; 2005.

Agarwala SS: New applications of cancer immunotherapy.
Semin Oncol 2002, 29(3 Suppl 7):14. PubMed Abstract

Bleumer I, Oosterwijk E, De Mulder P, Mulders PF: Immunotherapy for renal cell carcinoma.
Eur Urol 2003, 44:6575,. PubMed Abstract  Publisher Full Text

Kaminski JM, Summers JB, Ward MB, Huber MR, Minev B: Immunotherapy and prostate cancer.
Cancer Treat Rev 2003, 29:199209. PubMed Abstract  Publisher Full Text

De Pillis LG, Radunskaya AE, Wiseman CL: A Validated Mathematical Model of CellMediated Immune Response to Tumor Growth.
Cancer Res 2005, 65:79507958. PubMed Abstract  Publisher Full Text

Kirschner D, Panetta JC: Modeling immunotherapy of the tumor  immune interaction.
J Math Biol 1998, 37:235252. PubMed Abstract  Publisher Full Text

Kronik N, Kogan Y, Vainstein V, Agur Z: Improving alloreactive CTL immunotherapy for malignant gliomas using a simulation model of their interactive dynamics.
Cancer Immunol Immunother 2008, 57:425439. PubMed Abstract  Publisher Full Text

Cappuccio A, Elishmereni M, Agur Z: Cancer immunotherapy by interleukin21: potential treatment strategies evaluated in a mathematical model.
Cancer Res 2006, 66:72937300. PubMed Abstract  Publisher Full Text

Meloni G, et al.: How long can we give interleukin2? Clinical and immunological evaluation of AML patients after 10 or more years of IL2 administration.
Leukemia 2002, 16:20162018. PubMed Abstract  Publisher Full Text

De Lisi C, Rescigno A: Immune surveillance and neoplasia: a minimal mathematical model.
Bull Math Biol 1977, 39(2):201221. PubMed Abstract

Stepanova NV : Course of the immune reaction during the development of a malignant tumor.

Kuznetsov VA, Makalkin IA, Taylor MA, Perelson AS: Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis.
Bull Math Biol 1994, 56:295321. PubMed Abstract

Nani F, Freedman HI: A mathematical model of cancer treatment by immunotherapy.
Math Biosci 2000, 163:159199. PubMed Abstract  Publisher Full Text

Kuznetsov VA, Knott GD: Modeling tumor regrowth and immunotherapy.
Math Comp Mod 2001, 33:12751287. Publisher Full Text

de Vladar HP, González JA: Dynamic response of cancer under the influence of immunological activity and therapy.
J Theor Biol 2004, 227:335348. PubMed Abstract  Publisher Full Text

d'Onofrio A: A general framework for modeling tumorimmune system competition and immunotherapy: Mathematical analysis and biomedical inferences.
Phys D 2005, 208:220235. Publisher Full Text

Kirschner D, Tsygvintsev A: On the global dynamics of a model for tumor immunotherapy.
Math Biosci Eng 2009, 6(3):573583. PubMed Abstract

d'Onofrio A: Tumorimmune system interaction: modeling the tumorstimulated proliferation of effectors and immunotherapy.
Math Mod and Meth in App Sci 2006, 16:13751401. Publisher Full Text

d'Onofrio A, Gatti F, Cerrai P: Delayinduced oscillatory dynamics of tumorimmune system interaction.
Math and Comp Mod 2010, 51:572591. Publisher Full Text

d'Onofrio A: Boundednoiseinduced transitions in a tumorimmune system interplay.
Phys Rev E Stat Nonlin Soft Matter Phys 2010, 81:021923. PubMed Abstract  Publisher Full Text

Bellomo N, Delitala M: . From the mathematical kinetic, and stochastic game theory to modelling mutations, onset, progression and immune competition of cancer cells.
Phys of Life Rev 2008, 5:183206. Publisher Full Text

Bellomo N, Forni G: Complex multicellular systems and immune competition: New paradigms looking for a mathematical theory.
Current Topics In Developmental Biology 2008, 81:485502. PubMed Abstract  Publisher Full Text

Pappalardo F, Lollini PL, Castiglione F, Motta S: Modeling and simulation of cancer immunoprevention vaccine.
Bioinformatics 2005, 21(12):28912897. PubMed Abstract  Publisher Full Text

Pennisi M, Pappalardo F, Palladini A, Nicoletti G, Nanni P, Lollini PL, Motta S: Modeling the competition between lung metastases and the immune system using agents.
BMC Bioinformatics 2010, 11(Suppl 7):S13. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Kennedy BJ: . Cyclic leukocyte oscillations in chronic myelogenous leukemia during hydroxyurea therapy.
Blood 1970, 35:751760. PubMed Abstract  Publisher Full Text

Vodopick H, Rupp EM, Edwards CL, Goswitz FA, Beauchamp JJ: Spontaneous cyclic leukocytosis and thrombocytosis in chronic granulocytic leukemia.
N Engl J Med 1972, 286:284290. PubMed Abstract  Publisher Full Text

Gatti R, et al.: Cyclic Leukocytosis in Chronic Myelogenous Leukemia: New Perspectives on Pathogenesis and Therapy.
Blood 1973, 41:771783. PubMed Abstract  Publisher Full Text

Mehta BC, Agarwal MB: Cyclic oscillations in leukocyte count in chronic myeloid leukemia.
Acta Haematol 1980, 63:6870. PubMed Abstract  Publisher Full Text

Sohrabi A, Sandoz J, Spratt JS, Polk HC: Recurrence of breast cancer: Obesity, tumor size, and axillary lymph node metastases.
JAMA 1980, 244(3):264265. PubMed Abstract  Publisher Full Text

Tsao H, Cosimi AB, Sober AJ: Ultralate recurrence (15 years or longer) of cutaneous melanoma.
Cancer 1997, 79(12):23612370. PubMed Abstract  Publisher Full Text

Caravagna G, d'Onofrio A, Milazzo P, Barbuti R: Tumour suppression by immune system through stochastic oscillations.
J Theor Biol 2010, 265(3):336345. PubMed Abstract  Publisher Full Text

Guckenheimer J, Holmes PJ: Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields. Springer, New York; 1983.

Hale JK, Kocak H: Dynamics and Bifurcation. Springer, New York; 1991.

Gillespie DT: A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions.

Gillespie DT: Exact Stochastic Simulation of Coupled Chemical Reactions.
J of Phys Chem 1977, 81:23402361. Publisher Full Text

Salis H, Kaznessis Y: Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions.
J Chem Phys 2005, 122:054103. Publisher Full Text

Alfonsi A, Cances E, Turinici G, Di Ventura B, Huisinga W: Exact simulation of hybrid stochastic and deterministic models for biochemical systems.

Lecca P: A timedependent extension of Gillespie algorithm for biochemical stochastic πcalculus. In Proceedings of the 2006 ACM symposium on Applied computing 2327 April 2006 Dijon. Edited by Hisham M. Haddad:ACM; 2006:137144.

Anderson DF : A modified next reaction method for simulating chemical systems with time dependent propensities and delays.
J Chem Phys 2007, 127:214107. PubMed Abstract  Publisher Full Text

Kirschner D, Arciero JC, Jackson TL: A Mathematical Model of TumorImmune Evasion and siRNA Treatment.

DeBoer RJ, Hogeweg P, Hub F, Dullens J, DeWeger RA, DenOtter W: Macrophage T Lymphocyte interactions in the antitumor immune response: A mathematical model.
J Immunol 1985, 134:27482758. PubMed Abstract  Publisher Full Text

d'Onofrio A: Stability property of Pulse Vaccination Strategy in SEIR epidemic model.
Math Biosci 2002, 179(1):5772. PubMed Abstract  Publisher Full Text

d'Onofrio A: Pulse Vaccination Strategy in SIR epidemic model: global stability, vaccine failures and double vaccinations.
Math and Comp Model 2002, 36(45):461478. Publisher Full Text

Shahrezaei V, Ollivier JF, Swain PS: Colored extrinsic fluctuations and stochastic gene expression.
Mol Syst Biol 2008, 4:196. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Caravagna G: Formal Modeling and Simulation of Biological Systems With Delays. PhD thesis. Università di Pisa, Dipartimento di Informatica; 2011.

Barrio M, Burrage K, Leier A, Tian T: Oscillatory regulation of Hes1: discrete stochastic delay modelling and simulation.
PLoS Comput Biol 2006, 2(9):e117. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bratsun D, Volfson D, Tsimring LS, Hasty J: Delayinduced Stochastic Oscillations in Gene Regulation.
PNAS 2005, 102(41):1459314598. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Barbuti R, Caravagna G, MaggioloSchettini A, Milazzo P: On the Interpretation of Delays in Delay Stochastic Simulation of Biological Systems.
In Proceedings of the 2nd Int Workshop on Computational Models for Cell Processes (CompMod'09) EPTCS Edited by Back R, Petre I, de Vink EP. 2009, 6:1729.

Barbuti R, Caravagna G, MaggioloSchettini A, Milazzo P: Delay Stochastic Simulation of Biological Systems: A Purely Delayed Approach.

Levi F, Schibler U: Circadian Rhythms: Mechanisms and Therapeutic Implications.
Annu Rev Pharmacol Toxicol 2007, 47:593628. PubMed Abstract  Publisher Full Text

Kerbel RS, Kamen BA: The antiangiogenic basis of metronomic chemotherapy.
Nat Rev Cancer 2004, 4:423436. PubMed Abstract  Publisher Full Text

Orlando L, Cardillo A, Rocca A, Balduzzi A, Ghisini R, Peruzzotti G, Goldhirsch A, D'Alessandro C, Cinieri S, Preda L, Colleoni M: Prolonged clinical benefit with metronomic chemotherapy in patients with metastatic breast cancer.
Anticancer Drugs 2006, 17:961967. PubMed Abstract  Publisher Full Text

Gillespie DT: Approximated Accelerated Stochastic Simulation of Chemically Reacting Systems.
J of Chem Phys 2001, 115(4):17161733. Publisher Full Text

Gillespie DT, Petzold LR: Numerical Simulation for Biochemical Kinetics. In System modeling in cell biology: from concepts to nuts and bolts. MIT Press; 2001.

d'Onofrio A, Gandolfi A: Resistance to antitumor chemotherapy due to boundednoiseinduced transitions.
Phys Rev E Stat Nonlin Soft Matter Phys 2010, 82:061901. PubMed Abstract  Publisher Full Text

Daalhuis ABO: Hypergeometric function. [http://dlmf.nist.gov/15] webcite
In NIST Handbook of Mathematical Functions Edited by Olver FWJ, et al. Cambridge University Press; 2010.