Volume 13 Supplement 4
Italian Society of Bioinformatics (BITS): Annual Meeting 2011
Finetuning antitumor immunotherapies via stochastic simulations
 Giulio Caravagna^{1},
 Roberto Barbuti^{2} and
 Alberto d'Onofrio^{3}Email author
DOI: 10.1186/1471210513S4S8
© Caravagna et al.; licensee BioMed Central Ltd. 2012
Published: 28 March 2012
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 [10–12]. 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, 13–16, 18–28], 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 [33–38] 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
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 ${\mu}_{E}^{1}$ and the average degradation time for interleukin is ${\mu}_{I}^{1}$. 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.
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].
Birth and death reactions for the hybrid model
reaction  Propensity  reaction  propensity 

R_{1} : T ↦ T + 1  a_{1} = r_{2}T  R_{2} : T ↦ T  1  ${a}_{2}=\frac{{r}_{2}b}{V}{T}^{2}$ 
R_{3} : T ↦ T  1  ${a}_{3}=\frac{{p}_{T}TE}{{g}_{T}V+T}$  R_{4} : E ↦ E + 1  ${a}_{4}\left(t\right)=\frac{{p}_{E}EI\left(t\right)}{{g}_{E}+I\left(t\right)}$ 
R_{5} : E ↦ E  1  a_{5} = μ_{ E }E  R_{6} : E ↦ E + 1  a_{6} = cT 
R_{7} : E ↦ E + 1  a_{7} = Vσ_{ E } 
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
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.
where r is a random number U[0, 1] and a_{ i }(t + τ ) = a_{ i }(x) if $i\in \mathcal{S}$. 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.
The Hybrid Simulation Algorithm (Algorithm 1)
Require: (T_{0}, E_{0}, I_{0}), t_{0}, t_{ stop }. 

1: set the initial state to (T_{0}, E_{0}, I(t_{0})) and the initial time t to t_{0}; 
2: while t < t_{ stop } do 
3: let x be the current state, for $j\in \mathcal{S}$ evaluate a_{ j }(x), define ${a}_{0}^{\mathcal{S}}={\sum}_{j\in \mathcal{S}}{a}_{j}\left(\mathbf{x}\right)$; 
4: let χ be a random number with distribution Exp(1), solve the transcendental equation 
${A}_{4}\left(\tau \right)+{a}_{0}^{S}\tau +{A}_{7}\left(\tau \right)=\chi $ 
and then define ${a}_{0}={a}_{0}^{S}+{a}_{4}\left(t+\tau \right)+{a}_{7}\left(t+\tau \right)$; 
5: let r be a random number U[0, 1], for the next event to fire find j by solving 
$\sum _{i=1}^{j1}{a}_{i}\left(t+\tau \right)<r\cdot {a}_{0}\le \sum _{i=1}^{j}{a}_{i}\left(t+\tau \right)$ 
where if $i\in \mathcal{S}$ then a_{ i }(t + τ ) = a_{ i }(x); 
6: update (T, E, I(t)) to (T + ν_{T,j}, E + ν_{E,j}, I(t + τ )) and change clock to t + τ ; 
7: end while 
Values of the parameters
Values of the parameters
Par.  Value  Unit  Description 

r  0.18  days ^{1}  baseline growth rate of the tumor 
b  10^{9}  ml ^{1}  carrying capacity of the tumor 
a  1  ml/days  baseline strength of the killing rate by immune effectors 
c  10^{4}  days ^{1}  tumor antigenicity 
V  3.2  Ml  blood and bone marrow volumes for leukemia 
g _{ T }  10^{5}  ml ^{1}  50% reduction factor of the killing rate by immune effectors 
g _{ E }  2 · 10^{7}  pg/l  50% reduction factor of ILstimulated growth rate of effectors 
g _{ I }  10^{3}  ml ^{1}  50% reduction factor of production rate of interleukins 
p _{ E }  0.1245  days ^{1}  baseline strength of the ILstimulated growth rate of effectors 
p _{ I }  5  pg/days  baseline strength of production rate of interleukins 
μ _{ E }  0.03  days ^{1}  inverse of average lifespan of effectors 
μ _{ I }  10  days ^{1}  loss/degradation rate of IL_{2} 
Interleukinbased immunotherapies
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
if θ_{max+1}≤ t_{ n } + φ < θ_{max+1}+ A and 0 otherwise. This holds since in the uppermost case ${\int}_{{\theta}_{max}+A}^{{t}_{n}+\phi}\eta \left(x\right)dx={\int}_{{\theta}_{max+1}}^{{t}_{n}+\phi}\eta \left(x\right)dx$. 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
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
Impulsive therapy
The Hybrid Simulation Algorithm (Algorithm 2)
Require: (T_{0}, E_{0}, I_{0}), t_{0}, t_{ stop }. 

1: initialize the simulation as for Algorithm 1; 
2: while t < t_{ stop } do 
3: pick a value for τ as in Algorithm 1; 
4: get θ_{ next } = min{θ_{ i } ∈ Θ  θ_{ i } >t} 
5: if τ < θ_{ next } then 
6: fire a reaction as in Algorithm 1; 
7: else 
8: update (T, E, I(t)) to (T, E + w_{ next }, I(t + θ_{ next })) and change clock to t + θ_{ next }; 
9: end if 
10: end while 
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.
IL2 immunotherapy
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:
and γ_{1} < 1.
ACI
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.
if 0 ≤ t < P and the eradication condition can be consequently computed.
Combined therapies
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
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
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
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}.
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
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
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
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).
Averages and standard deviation, daily IL2 monotherapy
〈t_{ er }〉  σ  〈t_{ er }〉  σ 

19.34  0.33  22.71  0.41 
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.
Averages and standard deviation, daily ACI monotherapy
w _{*}  〈t_{ er }〉  σ  〈t_{ er }〉  σ 

5  15.5  1.1180  15.2  1.7204 
2.5  24.0  2.0  23.5  1.7078 
2  27.6  1.9720  27  2 
1.5  35.6  3.0397  34.7  3.1638 
1  62.0  4.3204  61.0  4.3204 
Averages and standard deviation, weekly ACI monotherapy
w _{*}  〈t_{ er }〉  σ  〈t_{ er }〉  σ 

5  14.98  0.80  11.39  0.93 
2.5  23.12  1.22  19.27  1.30 
2  26.81  1.28  22.95  1.23 
1.5  33.66  2  28.9  2.42 
1  61.18  3.75  53.19  5.70 
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.
Declarations
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.
Authors’ Affiliations
References
 d'Onofrio A, Gandolfi A, Rocca A: The dynamics of tumorvasculature interaction suggests lowdose, timedense antiangiogenic schedulings. Cell Prolif 2009, 42: 317–329. 10.1111/j.13652184.2009.00595.xView ArticlePubMedGoogle Scholar
 Pardoll D: Does the immune system see tumors as foreign or self? Annu Rev Immunol 2003, 21: 807–839. 10.1146/annurev.immunol.21.120601.141135View ArticlePubMedGoogle Scholar
 Ehrlich P: Ueber den jetzigen Stand der Karzinomforschung. Ned Tijdschr Geneeskd 1909, 5: 273–290.Google Scholar
 Dunn GP, Old LJ, Schreiber RD: The three Es of cancer immunoediting. Annu Rev Immunol 2004, 22: 329–360. 10.1146/annurev.immunol.22.012703.104803View ArticlePubMedGoogle Scholar
 d'Onofrio A: Tumor evasion from immune system control: strategies of a MISS to become a MASS. Chaos Solitons and Fractals 2007, 31: 261–268. 10.1016/j.chaos.2005.10.006View ArticleGoogle Scholar
 Whiteside TL: Tumorinduced death of immune cells: its mechanisms and consequences. Semin Cancer Biol 2002, 12: 43–50. 10.1006/scbi.2001.0402View ArticlePubMedGoogle Scholar
 Vicari AP, Caux C, Trinchieri G: Tumor escape from immune surveillance through dendritic cell inactivation. Semin Cancer Biol 2002, 12: 33–42. 10.1006/scbi.2001.0400View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 De Vito VT Jr, Hellman J, Rosenberg SA: Cancer: Principles and Practice of Oncology. 9th North American Edition edition. Philadelphia: JP Lippincott; 2005.Google Scholar
 Agarwala SS: New applications of cancer immunotherapy. Semin Oncol 2002, 29(3 Suppl 7):1–4.View ArticlePubMedGoogle Scholar
 Bleumer I, Oosterwijk E, De Mulder P, Mulders PF: Immunotherapy for renal cell carcinoma. Eur Urol 2003, 44: 65–75,. 10.1016/S03022838(03)00191XView ArticlePubMedGoogle Scholar
 Kaminski JM, Summers JB, Ward MB, Huber MR, Minev B: Immunotherapy and prostate cancer. Cancer Treat Rev 2003, 29: 199–209. 10.1016/S03057372(03)000057View ArticlePubMedGoogle Scholar
 De Pillis LG, Radunskaya AE, Wiseman CL: A Validated Mathematical Model of CellMediated Immune Response to Tumor Growth. Cancer Res 2005, 65: 7950–7958.PubMedGoogle Scholar
 Kirschner D, Panetta JC: Modeling immunotherapy of the tumor  immune interaction. J Math Biol 1998, 37: 235–252. 10.1007/s002850050127View ArticlePubMedGoogle Scholar
 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: 425–439. 10.1007/s002620070387zView ArticlePubMedGoogle Scholar
 Cappuccio A, Elishmereni M, Agur Z: Cancer immunotherapy by interleukin21: potential treatment strategies evaluated in a mathematical model. Cancer Res 2006, 66: 7293–7300. 10.1158/00085472.CAN060241View ArticlePubMedGoogle Scholar
 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: 2016–2018. 10.1038/sj.leu.2402566View ArticlePubMedGoogle Scholar
 De Lisi C, Rescigno A: Immune surveillance and neoplasia: a minimal mathematical model. Bull Math Biol 1977, 39(2):201–221.Google Scholar
 Stepanova NV : Course of the immune reaction during the development of a malignant tumor. Bioph 1980, 24: 917–923.Google Scholar
 Kuznetsov VA, Makalkin IA, Taylor MA, Perelson AS: Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis. Bull Math Biol 1994, 56: 295–321.View ArticlePubMedGoogle Scholar
 Nani F, Freedman HI: A mathematical model of cancer treatment by immunotherapy. Math Biosci 2000, 163: 159–199. 10.1016/S00255564(99)000589View ArticlePubMedGoogle Scholar
 Kuznetsov VA, Knott GD: Modeling tumor regrowth and immunotherapy. Math Comp Mod 2001, 33: 1275–1287. 10.1016/S08957177(00)003149View ArticleGoogle Scholar
 de Vladar HP, González JA: Dynamic response of cancer under the influence of immunological activity and therapy. J Theor Biol 2004, 227: 335–348. 10.1016/j.jtbi.2003.11.012View ArticlePubMedGoogle Scholar
 d'Onofrio A: A general framework for modeling tumorimmune system competition and immunotherapy: Mathematical analysis and biomedical inferences. Phys D 2005, 208: 220–235. 10.1016/j.physd.2005.06.032View ArticleGoogle Scholar
 Kirschner D, Tsygvintsev A: On the global dynamics of a model for tumor immunotherapy. Math Biosci Eng 2009, 6(3):573–583.View ArticlePubMedGoogle Scholar
 d'Onofrio A: Tumorimmune system interaction: modeling the tumorstimulated proliferation of effectors and immunotherapy. Math Mod and Meth in App Sci 2006, 16: 1375–1401. 10.1142/S0218202506001571View ArticleGoogle Scholar
 d'Onofrio A, Gatti F, Cerrai P: Delayinduced oscillatory dynamics of tumorimmune system interaction. Math and Comp Mod 2010, 51: 572–591. 10.1016/j.mcm.2009.11.005View ArticleGoogle Scholar
 d'Onofrio A: Boundednoiseinduced transitions in a tumorimmune system interplay. Phys Rev E Stat Nonlin Soft Matter Phys 2010, 81: 021923.View ArticlePubMedGoogle Scholar
 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: 183–206. 10.1016/j.plrev.2008.07.001View ArticleGoogle Scholar
 Bellomo N, Forni G: Complex multicellular systems and immune competition: New paradigms looking for a mathematical theory. Current Topics In Developmental Biology 2008, 81: 485–502.View ArticlePubMedGoogle Scholar
 Pappalardo F, Lollini PL, Castiglione F, Motta S: Modeling and simulation of cancer immunoprevention vaccine. Bioinformatics 2005, 21(12):2891–2897. 10.1093/bioinformatics/bti426View ArticlePubMedGoogle Scholar
 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. 10.1186/1471210511S7S13PubMed CentralView ArticlePubMedGoogle Scholar
 Kennedy BJ: . Cyclic leukocyte oscillations in chronic myelogenous leukemia during hydroxyurea therapy. Blood 1970, 35: 751–760.PubMedGoogle Scholar
 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: 284–290. 10.1056/NEJM197202102860603View ArticlePubMedGoogle Scholar
 Gatti R, et al.: Cyclic Leukocytosis in Chronic Myelogenous Leukemia: New Perspectives on Pathogenesis and Therapy. Blood 1973, 41: 771–783.PubMedGoogle Scholar
 Mehta BC, Agarwal MB: Cyclic oscillations in leukocyte count in chronic myeloid leukemia. Acta Haematol 1980, 63: 68–70. 10.1159/000207373View ArticlePubMedGoogle Scholar
 Sohrabi A, Sandoz J, Spratt JS, Polk HC: Recurrence of breast cancer: Obesity, tumor size, and axillary lymph node metastases. JAMA 1980, 244(3):264–265. 10.1001/jama.1980.03310030040023View ArticlePubMedGoogle Scholar
 Tsao H, Cosimi AB, Sober AJ: Ultralate recurrence (15 years or longer) of cutaneous melanoma. Cancer 1997, 79(12):2361–2370. 10.1002/(SICI)10970142(19970615)79:12<2361::AIDCNCR10>3.0.CO;2PView ArticlePubMedGoogle Scholar
 Caravagna G, d'Onofrio A, Milazzo P, Barbuti R: Tumour suppression by immune system through stochastic oscillations. J Theor Biol 2010, 265(3):336–345. 10.1016/j.jtbi.2010.05.013View ArticlePubMedGoogle Scholar
 Guckenheimer J, Holmes PJ: Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields. Springer, New York; 1983.View ArticleGoogle Scholar
 Hale JK, Kocak H: Dynamics and Bifurcation. Springer, New York; 1991.View ArticleGoogle Scholar
 Gillespie DT: A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. J of Comp Phys 19t6, 22(4):403–434.View ArticleGoogle Scholar
 Gillespie DT: Exact Stochastic Simulation of Coupled Chemical Reactions. J of Phys Chem 1977, 81: 2340–2361. 10.1021/j100540a008View ArticleGoogle Scholar
 Salis H, Kaznessis Y: Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions. J Chem Phys 2005, 122: 054103. 10.1063/1.1835951View ArticleGoogle Scholar
 Alfonsi A, Cances E, Turinici G, Di Ventura B, Huisinga W: Exact simulation of hybrid stochastic and deterministic models for biochemical systems. INRIA Tech Report 2004, 5435.Google Scholar
 Lecca P: A timedependent extension of Gillespie algorithm for biochemical stochastic π calculus. In Proceedings of the 2006 ACM symposium on Applied computing 23–27 April 2006 Dijon. Edited by: Hisham M. Haddad:ACM; 2006:137–144.Google Scholar
 Anderson DF : A modified next reaction method for simulating chemical systems with time dependent propensities and delays. J Chem Phys 2007, 127: 214107. 10.1063/1.2799998View ArticlePubMedGoogle Scholar
 Kirschner D, Arciero JC, Jackson TL: A Mathematical Model of TumorImmune Evasion and siRNA Treatment. Discr and Cont Dyn Systems 2004, 4: 39–58.View ArticleGoogle Scholar
 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: 2748–2758.Google Scholar
 d'Onofrio A: Stability property of Pulse Vaccination Strategy in SEIR epidemic model. Math Biosci 2002, 179(1):57–72. 10.1016/S00255564(02)000950View ArticlePubMedGoogle Scholar
 d'Onofrio A: Pulse Vaccination Strategy in SIR epidemic model: global stability, vaccine failures and double vaccinations. Math and Comp Model 2002, 36(4–5):461–478. 10.1016/S08957177(02)001760View ArticleGoogle Scholar
 Shahrezaei V, Ollivier JF, Swain PS: Colored extrinsic fluctuations and stochastic gene expression. Mol Syst Biol 2008, 4: 196.PubMed CentralView ArticlePubMedGoogle Scholar
 Caravagna G: Formal Modeling and Simulation of Biological Systems With Delays. PhD thesis. Università di Pisa, Dipartimento di Informatica; 2011.Google Scholar
 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. 10.1371/journal.pcbi.0020117PubMed CentralView ArticlePubMedGoogle Scholar
 Bratsun D, Volfson D, Tsimring LS, Hasty J: Delayinduced Stochastic Oscillations in Gene Regulation. PNAS 2005, 102(41):14593–14598. 10.1073/pnas.0503858102PubMed CentralView ArticlePubMedGoogle Scholar
 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: 17–29.Google Scholar
 Barbuti R, Caravagna G, MaggioloSchettini A, Milazzo P: Delay Stochastic Simulation of Biological Systems: A Purely Delayed Approach. Trans on Comp Sys Bio 2011, XIII 6575: 61–84.View ArticleGoogle Scholar
 Levi F, Schibler U: Circadian Rhythms: Mechanisms and Therapeutic Implications. Annu Rev Pharmacol Toxicol 2007, 47: 593–628. 10.1146/annurev.pharmtox.47.120505.105208View ArticlePubMedGoogle Scholar
 Kerbel RS, Kamen BA: The antiangiogenic basis of metronomic chemotherapy. Nat Rev Cancer 2004, 4: 423–436. 10.1038/nrc1369View ArticlePubMedGoogle Scholar
 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: 961–967. 10.1097/01.cad.0000224454.46824.fcView ArticlePubMedGoogle Scholar
 Gillespie DT: Approximated Accelerated Stochastic Simulation of Chemically Reacting Systems. J of Chem Phys 2001, 115(4):1716–1733. 10.1063/1.1378322View ArticleGoogle Scholar
 Gillespie DT, Petzold LR: Numerical Simulation for Biochemical Kinetics. In System modeling in cell biology: from concepts to nuts and bolts. MIT Press; 2001.Google Scholar
 d'Onofrio A, Gandolfi A: Resistance to antitumor chemotherapy due to boundednoiseinduced transitions. Phys Rev E Stat Nonlin Soft Matter Phys 2010, 82: 061901.View ArticlePubMedGoogle Scholar
 Daalhuis ABO: Hypergeometric function.In NIST Handbook of Mathematical Functions Edited by: Olver FWJ, et al. Cambridge University Press; 2010. [http://dlmf.nist.gov/15]Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.