Skip to main content

Fine-tuning anti-tumor immunotherapies via stochastic simulations

Abstract

Background

Anti-tumor 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 side-effects, 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 therapy-related 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 discrete-state continuous-time 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 Interleukin-based 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, piece-wise 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 therapy-induced tumor eradication time. Indeed, our model suggests that Interleukin-based therapies may not be effective for every patient, and that the piece-wise constant is the most effective delivery to stimulate the immune-response. For Adoptive Cellular Immunotherapies a metronomic delivery seems more effective, as it happens for other anti-angiogenesis therapies and chemotherapies, and the impulsive delivery seems more effective than the piece-wise constant. The expected synergistic effects have been observed when the therapies are combined.

Introduction

A key problem of anti-tumor therapies is finding a suitable scheduling of their administration to the patients. Of course a major problem in medical oncology is avoiding severe therapy-related side effects which, unfortunately, may cause the death of the patient. However, also in the ideal case of no side-effects, 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 neo-antigens. Such antigens trigger anti-tumor actions by the immune system [2], thus resulting in the tumor-immune 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 inter-subjects 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 Interleukin-based 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 agent-based models have been introduced by Motta and Lollini [31, 32]. In [14] Kirschner and Panetta proposed a largely influential ODE-based model of Tumor-Immune system (T-IS) interplay, whose variables are tumor cells, effector cells and the concentration of interleukins-2 (IL-2). 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 Kirschner-Panetta (KP) model, the tumor-free 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 discrete-state continuous-time 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 interleukin-based 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 Kirschner-Panetta model

In [14] the following model of the dynamics of tumor-immune system interaction was proposed

T ′ * = r T * ( 1 - b T * ) - p T T * g T + T * E * E ′ * = p E I g E + I E * - μ E E * + c T * + σ E I ′ = p I T * E * g I + T * - μ I I + σ I
(1)

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 μ E - 1 and the average degradation time for interleukin is μ 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 tumor-free 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 side-effects.

A hybrid model with constant therapies

As discussed in [39], the low-level 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 IL-2 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 T-IS 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

T ′ = r T 1 - b V T - p T T g T V + T E E ′ = p E I g E + I E - μ E E + c T + V σ E I ′ = p I V T E g I V + T - μ I I + σ I .
(2)

Note that the modified deterministic model (2) is obtained by the original Kirschener-Panetta deterministic model (1) by means of a nonsingular linear transformation. As it is very well-known 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 bi-dimensional 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 continuous-time 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 R. 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 S = R \ R 4 to each reaction R i ∈S 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 S are time-independent, meaning that they depend on a state which is constant in between two stochastic events. Differently, R4 depends on the continuous part of the hybrid model, i.e. a4 (t) depends on I (t), and is hence time-dependent.

Table 1 Birth and death reactions for the hybrid model

We remark that the original SSA is assumed to simulate only time-independent chemical reactions. The algorithm used to simulate this model takes inspiration from algorithms simulating hybrid systems with time-independent propensity functions [44, 45] and algorithms simulating purely stochastic systems with time-dependent 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

τ ∑ j ∈ S a j ( x ) + ∫ t n t n + τ a 4 ( t ) dt=χ
(3)

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 time-independent 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 Tumor-Immune 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 time-dependent immunotherapies, that is

T ′ = r T 1 - b V T - p T T g T V + T E E ′ = p E I g E + I E - μ E E + c T + V σ E ( t ) I ′ = p I V T E g I V + T - μ I I + σ I ( t ) .
(4)

It is important to notice that in the realistic case of finite duration therapies the deterministic system always predicts tumor re-growth, being the tumor-free equilibrium (0, 0, 0) unstable. Practically, if at the end of the therapy the solution is very close to the tumor-free 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 follow-up visits, i.e. that the tumor has not re-grown. 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 low-level 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 therapy-induced 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 SSA-based 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 R7 is modified as a7(t) = V σ E (t) to reflect the time-dependency of σ E (t). As a consequence in this case we have two time-dependent propensity functions a4 (t) and a7 (t), hence we define S=R\ R 4 , R 7 . Let us model the event-induced changes in the values of T and E by the vectors

ν T ν E = 1 - 1 - 1 0 0 0 0 0 0 0 1 - 1 1 1 .

Here the i-th components of ν T and ν E , i.e. νT,iand ν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

A 4 ( τ ) = ∫ t n t n + τ a 4 ( t ) d t

and

A 7 ( τ ) = ∫ t n t n + τ a 7 ( t ) d t = V ∫ t n t n + τ σ E ( t ) d t .

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 S evaluated, i.e. a 0 S = ∑ j ∈ S a j ( x ) the putative time for the next stochastic event follows

A 4 ( τ ) + a 0 S τ + A 7 ( τ ) = χ
(5)

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 a 0 = a 0 S + a 4 ( t + τ ) + a 7 ( t + τ ) , the next event to fire find is R j if

∑ i = 1 j - 1 a i ( t + τ ) < r ⋅ a 0 ≤ ∑ i = 1 j a i ( t + τ )
(6)

where r is a random number U[0, 1] and a i (t + τ ) = a i (x) if i∈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.

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 tn+1the equation for I' reads as

I ′ = p I V T ( t n ) E ( t n ) g I V + T ( t n ) + σ I ( t ) - μ I I
(7)

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 , tn+1) is

I ( t ) = B n + ( I n - B n ) e - μ I ( t - t n ) + e - μ I t ∫ t n t σ I ( z ) e μ I Z d z
(8)

where

B n = 1 μ I p I V T ( t n ) E ( t n ) g I V + T ( t n ) .
(9)

Notice that equation (8) is necessary to evaluate, at each step of Algorithm 1, the value of a4(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

Interleukin-based immunotherapies

In the next paragraphs we discuss how to specialize model (4) with either a piece-wise constant or an impulsive interleukin-based 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 IL-2. 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 pre-set times

Θ = { θ i | i = 0 , … , k }

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.

Piece-wise constant therapy

In the case of a continuous infusion each therapy session has a duration of A time-units and a constant influx rate d i in the i-th infusion. We model the therapy as

σ I ( t ) = 0 if t s < t < t e ∑ i = 0 k d i g ( t - θ i , A ) otherwise
(10)

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.

g ( t , A ) = H ( t ) - H ( t - A )

where H(·) the well-known 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.

R ( t n , t ) = ∫ t n t σ I ( x ) e μ I x d x = ∫ t n t n + φ η ( x ) d x

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

θ ( t n , φ ) = { θ j ∈ Θ | t n < θ j ∧ θ j + A < t n + φ } .

Now, since Θ contains ordered time-points then also θ(t n , φ) does. We want to get the index of the minimum and the maximum θ j 's from θ(t n , φ); we define

θ m i n = min ( θ ( t n , φ ) ) θ m a x = max ( θ ( t n , φ ) )

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

R ( t n , t n + φ ) = ∫ t n θ m i n η ( x ) d x + ∫ θ m i n θ m a x + A η ( x ) d x + ∫ θ m a x + A t n + φ η ( x ) d x .

The integral in [θ min , θ max + A] is the summation of |θ(t n , φ)| non-zero basic integrals over the sub-intervals [θ j , θ j + A] for any θ j ∈ θ(t n , φ), that is

∫ θ m i n θ m a x + A η ( x ) d x = ∑ θ j ∈ θ ( t n , φ ) ∫ θ j θ j + A g ( x ) d x .

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 ∫ θ j + A θ j + 1 η ( x ) dx=0. Furthermore, since

∫ θ j θ j + A σ I ( x ) e μ I x d x = d j μ I e μ I x θ j θ j + A = d j e μ I θ j ( e μ I A - 1 ) μ I

then the overall integral evaluates as

∫ θ m i n θ m a x + A η ( x ) d x = e μ I A - 1 μ I ∑ θ j ∈ θ ( t n , φ ) d j e μ I θ j .

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 therapy-session. We have this definition

∫ t n θ m i n η ( x ) d x = d m i n - 1 μ I e μ I ( θ m i n - 1 + A ) - e μ I t n

if θmin-1< t n ≤ θmin-1+ A and 0 otherwise. This holds since in the uppermost case ∫ t n θ m i n η ( x ) dx= ∫ t n θ m i n - 1 + A η ( x ) dx. Similarly, in the interval (θ max , t n + φ] by cases on the relation between t n + φ and θmax+1

∫ θ m a x + A t n + φ η ( x ) d x = d m a x + 1 μ I ( e - μ I t - e μ I θ m a x + 1 )

if θmax+1≤ t n + φ < θmax+1+ A and 0 otherwise. This holds since in the uppermost case ∫ θ m a x + A t n + φ η ( x ) dx= ∫ θ m a x + 1 t n + φ η ( x ) 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

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.

σ I ( t ) = 0 if t s < t < t e ∑ i = 0 N u i δ ( t - θ i ) , otherwise .
(11)

Here, u i is the i-th injected dose of molecules of IL-2 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

R ( t n , t ) = ∑ i = i s i e u i e μ I θ i H ( t - θ i )

where { θ i s , … , θ i e } =Θ∩ ( t n , t ) . 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

I ( θ i + ) = I ( θ i - ) + u i
(12)

for i = 0, . . . , k.

Adoptive cellular immunotherapies

In the next paragraphs we discuss how to specialize model (4) with either a piece-wise constant or an impulsive ACI. As in the previous section we consider a set of k therapy sessions Θ = {θ i | i = 0, . . . , k}.

Piece-wise constant therapy

As before we assume

σ E ( t ) = 0 if t s < t < t e ∑ i = 0 k d i g ( t - θ i , A ) otherwise .
(13)

This scenario can be simulated by Algorithm 1 where, in this case, the function A7(Ï„) is given by

A 7 ( τ ) = V Ψ ( t ) + V ∑ i = i s i e d i ( t - θ i ) H ( t - θ i ) - V ∑ i = i s i e d i ( t - θ i - A ) H ( t - θ i - A )

where { θ i s , … , θ i e } =Θ∩ ( t n , t ) and

Ψ ( t ) = H ( θ i s - 1 + A - t n ) d i s - 1 ( t - t n ) H ( t - θ i s - 1 ) - ( t - θ i s - 1 - A ) H ( t - θ i s - 1 - A ) .

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 IL-based therapy. However, here the introduction of generalized functions does not lead to simplifications, and as a consequence we shall model the therapy as

E ( θ i + ) = E ( θ i - ) + w i .
(14)

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 n-th stochastic event and θ j >t n is the closest injection after t n , the propensity functions a3 and a5 change after θ j , invalidating the property. Of course, the time-dependent propensity functions a4 and a7 are, by definition, potentially non-constant. We argue that there is a similarity between scheduling-based 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 [53–57]. As for such algorithms, the correctness of Algorithm 2 relies on the existence of a general schema of SSA-based algorithm with piece-wise constant time-dependent propensity functions [52] where is stated tha, and this algorithm respects that schema.

Table 4 The Hybrid Simulation Algorithm (Algorithm 2)

Combining IL-2 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 IL-2 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.

IL-2 immunotherapy

We start from the case of the delivery of a IL-based mono-therapy (i.e. σ E (t) = 0). By setting T = E = 0 it follows I(t) = J∞(t), where J∞(t) is the asymptotic solution of

J ′ = - μ I J + σ I ( t ) .
(15)

Thus, we found a tumor-free state

T F I = ( 0 , 0 , J ∞ ( t ) )

which is, however, unstable. In fact, this follows by setting

( T , E , I ) = ( 0 , 0 , J ∞ ( t ) ) + ( T l , E l , I l ) ,

and linearizing since the equation for tumor cells reads as T l ′ =r T l . Even though this means that the tumor-free state (0, 0, J∞(t)) is unstable, this does not mean that under an IL-based mono-therapy 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

I ′ ≥ - μ I I + σ I ( t )

it follows that I(t) >J(t) and, for large times, I(t) ≥ J∞(t). In turn, the inequality

p E I ( t ) g E + I ( t ) > p E J ∞ ( t ) g E + J ∞ ( t )

yields the inequality

E ′ ≥ p E J ∞ ( t ) g E + J ∞ ( t ) - μ E E .
(16)

Finally, from the properties of periodical linear differential equations recalled the additional file 1 we can observe that

p E J ∞ ( t ) g E + J ∞ ( t ) ≥ μ E
(17)

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) piece-wise 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

p E σ I g E μ I + σ I ≥ μ E .

(ii) if we consider the piece-wise 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

J ∞ ( t ) = d μ I - β 1 e - μ I t if 0 ≤ t < A β 2 e - μ I ( t - A ) if 0 ≤ t < P
(18)

where

β 1 = d μ I e μ I P - e μ I A e μ I P - 1 β 2 = d μ I - β 1 e - μ I A .

In this case condition (17) reads as

p E P ( 1 - ω 1 - ω 2 ) ≥ μ E

where

ω 1 = g E μ I ( g E + α 1 ) log e μ I A - γ 1 1 - γ 1 ω 1 = 1 μ I log e μ I ( P - A ) + γ 2 1 + γ 2 γ 1 = β 1 g E + α 1 γ 2 = β 2 g E

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

J ∞ ( t ) = U 1 - e - μ I P e - μ I ( t mod P )

where t mod P is the standard modulo operation. Then the eradication condition (17) is

p E P μ I log u e μ I P + g E ( e μ I P - 1 ) u + g E ( e μ I P - 1 ) > μ E .
(19)

ACI

When only ACI is delivered (i.e. σ I (t) = 0) there is the following tumor-free asymptotic solution

( 0 , ε E ( t ) , 0 )

where ε E (t) is the asymptotic solution of

ε ′ ( t ) = - μ E ε + σ E ( t )
(20)

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

T l ′ = r - a g T V ε E ( t ) τ ι

and hence the local eradication condition (17) is

a g T V ⟨ ε E ( t ) ⟩ ≥ r .
(21)

However, by averaging both the sides of (20) yields

ε E ( P ) - ε E ( 0 ) = 0 = - μ E ⟨ ε E ( t ) ⟩ + σ E ( t )

which implies that

⟨ ε E ( t ) ⟩ = ⟨ σ E ( t ) ⟩ μ E .

Finally, in this case the local stability condition corresponding to equation (17) becomes

a g T V ⟨ σ E ( t ) ⟩ μ E ≥ r .
(22)

As for the case of IL-2 mono-therapy 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 piece-wise 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 IL-based mono-therapy one can show that

ε ∞ ( t ) = b μ E - b μ E ξ e - μ E t

where

ξ = e μ E P - e μ E A e μ E P - 1

if 0 ≤ t < A and

ε ∞ ( t ) = b μ E - b μ E ξ e - μ E A e - μ E ( t - A )

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

ε E ( t ) = w 1 - e - μ E P e - μ E P

and hence the eradication condition (22) is

a g T V w P μ E ≥ r .
(23)

Combined therapies

Finally, we shortly consider when both the therapies are delivered. In this case there is the following tumor-free asymptotic solution

( 0 , ε ∞ ( t ) , J ∞ ( t ) )

where ε∞(t) is the asymptotic solution of

ε ′ ( t ) = p E J ∞ ( t ) g E + J ∞ ( t ) - μ E ε + σ E ( t ) .
(24)

In the case of synchronous delivery with common period P the local eradication condition becomes

a g T V ⟨ ε ∞ ( t ) ⟩ ≥r.
(25)

As for the mono-thrapeutic 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

E ′ ( t ) > p E J ∞ ( t ) g E + J ∞ ( t ) - μ E E + σ E ( t )

follows that, for large times, E(t) ≥ ε∞(t). Then by

T ′ < r T 1 - b V T - a T g T V + T ε ∞ ( t )

it follows that T(t) < X(t), where

X ′ = r X 1 - b V X - a X g T V + x ε ∞ ( t ) .

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 IL-2 mono-therapy 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 IL-2 mono-therapy might be effective.

Further discussions are worth. In the case of impulsive therapy unless IL-2 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 piece-wise continuous delivery of IL-2 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 rythms-tuned delivery of chemotherapy some special 24-hours infusors have been experimented [58]. Roughly speaking, the above fact might be related to the "indirect effect" of IL-based 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 mono-therapy 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 non-small 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 single-month 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 x-axis and the y-axis, respectively. To perform simulations a JAVA implementation of model (4) and Algorithms 1 and 2 has been developed.

IL-2 immunotherapy

In Figure 1 a single stochastic run of the hybrid model with only IL-2 based immunotherapies is shown. We simulated both piece-wise constant (left) and impulsive (right) IL-2 daily immunotherapies with initial configuration T(0) = 105 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 · 107/A for i = 1, . . . , 30. In the latter case at each therapy session the value u i for the injection is u i = 4·107 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
figure 1

Single run, IL-2 mono-therapy. Single-run of piece-wise constant (left) and impulsive (right) IL-2 daily immunotherapy. In (left) A = 0.2 and d i = 4 · 107/A for i = 1, . . . , 30. In (right) u i = 4 · 107 for i = 1, . . . , 30. In both cases T(0) = 105, 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 piece-wise constant immunotherapy seems more efficient than the impulsive one. Indeed (i) in the piece-wise case the eradication of the tumor is reached at t e ≈ 20 days and the maximum tumor size is around 106 = 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 · 106 in both cases, whereas the density of IL-2 is of the order of 107 in (left) and 104 in (right).

Adoptive Cellular Immunotherapy

In Figure 2 a single stochastic run of the hybrid model with only ACI is shown. As for IL-based therapy, we simulated both piece-wise constant (left) and impulsive (right) daily ACIs with initial configuration T(0) = 105 and E(0) = I(0) = 0. In the former case each therapy session lasts A = 0.2 days and b i = 25 · 104 for i = 1, . . . , 30. In the latter case at each therapy session 5·104 effector cells are injected, i.e. w i = 5 · 104 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
figure 2

Single run, daily ACI mono-therapy. Single-run of piece-wise constant (left) and impulsive (right) daily ACI. In (left) A = 0.2 and b i = 5 · 104/(V A) for i = 1, . . . , 30. In (right) w i = 5 · 104 for i = 1, . . . , 30. In both cases T(0) = 105 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 · 105, and the density of IL-2 is of the order of 102.

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 · 104 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
figure 3

Single run, weekly ACI mono-therapy. Single-run of piece-wise 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 · 104/V A i = 0, . . . , 3. In (right) w i = 35 · 104 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 IL-2 are similar to those in Figure 2 (right).

Combined therapies

In Figure 4 a single stochastic run of the hybrid model with combined impulsive IL-2 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. Θ i I L =i+ 0 .5 and Θ i E =i for i = 1, . . . , 30. Again, the initial configuration is T(0) = 105, E(0) = I(0) = 0. The parameters for the therapies dosage and duration are the same for both cases. As far as IL-2 therapy is concerned, at each therapy session the value u i for the injection is u i = 106 for i = 1, . . . , 30. Differently, in each ACI session 104 effector cells are injected, i.e. w i = 104 for i = 1, . . . , 30. Notice that both the values are smaller than those used in the corresponding mono-therapy scenarios of Figures 1 (right) and 2 (right). Moreover, as in the case of single IL-2 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
figure 4

Single run, combined immunotherapies. Single-run of synchronous (left) and asynchronous (right) combined impulsive IL-2 and ACI daily immunotherapies. The asynchronous delivery is a shift of 0.5 days. In (left) u i = 106 for i = 1, . . . , 30 and in (right) w i = 104 for i = 1, . . . , 30. In both cases T(0) = 105, 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 · 105 and the density of IL-2 is around 10. In both therapies the maximum size of the tumor is almost equal, i.e. it is around 2.5·105. 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) = 106 (left) and T(0) = 107 (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 · 104 effector cells are injected, i.e. w i = 5 · 104 for i = 1, . . . , 30, and the eradication is reached at ≈ 24.5 days. On the contrary, in the case where T(0) = 107 and the same schedule is applied the eradication was not observed and the tumor size reached, quite rapidly, a size of the order of 109. In order to obtain eradication also for T(0) = 107, 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 · 104 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
figure 5

Single run, ACI mono-therapy and larger tumor. Single-run of impulsive daily ACI with T(0) = 106 (left) and T(0) = 107 (right). In both cases E(0) = I(0) = 0. In (left) w i = 5 · 104 for i = 1, . . . , 30. In (right) w i = 20 · 105 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 106 (left) and 108 (right), and the density of IL-2 of the order of 102 (left) and 107 (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 y-axes are different since the maximum in (right) is 100 times bigger than the maximum in (left).

Probabilistic analysis of IL-2 therapy

As we already said some of the simulations we performed are time-consuming, 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, single-run 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 piece-wise constant and the impulsive IL-2 immunotherapies of Figure 1. This was possible since their running time is of the order of some minutes.

Figure 6
figure 6

Probability density function, IL-2 mono-therapy. Empirical evaluation of ϱ(t er ) for different piece-wise constant (left) and impulsive (right) daily IL-2 immunotherapies. In (left) A = 0.2 and d i = 4 · 107/A for i = 1, . . . , 30. In (right) u i = 4 · 107 for i = 1, . . . , 30. In both cases T(0) = 105, 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 102 simulations for each configuration.

We defined the following time-dependent property over a single simulation: we want to evaluate

t e r = min { t ≤ 70 | T ( t ) = 0 }

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 102 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 · 107/A for i = 1, . . . , 30 whereas in (right) u i = 4 · 107 for i = 1, . . . , 30. In both cases T(0) = 105, 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 IL-2 mono-therapy, Figure 6 shows the evaluation of ϱ(t er ) for piece-wise constant (left) and impulsive therapies (right). It is remarkable that for all the simulations the eradication is always reached (i.e. 102 times out of 102 simulations).

In Table 5 the average of t er , denoted 〈t er 〉 and the standard deviation σ of t er are evaluated for the piece-wise constant (left) and impulsive (right) IL-2 immunotherapies of Figure 6. By means of the Wilcoxon statistical test, the non-parametric equivalent of the T-test, 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 IL-2 mono-therapy

Probabilistic analysis of ACI

As for the case of IL-2 mono-therapy the results on ACIs in Figure 2 motivated us to to investigat the relationship between impulsive and piece-wise constant ACIs, as well as the influence of the period P between two consecutive therapy sessions. Again, we considered the same time-dependent 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 102 simulations for each of a set of parameter configurations. In all simulations we used the initial configuration T(0) = 105 and E(0) = I(0) = 0.

Figure 7 shows the evaluation of ϱ(t er ) for a daily ACIs with piece-wise 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* · 104/(V A) for i = 1, . . . , 30 and in (right) w i = w* · 104 for i = 1, . . . , 30. The other parameters are in Table 3.

Figure 7
figure 7

Probability density function, daily ACI mono-therapy. Empirical evaluation of ϱ(t er ) for different piece-wise constant (left) and impulsive (right) daily ACIs. In (left) b i = w* · 104/(V A) for i = 1, . . . , 30 and in (right) w i = w* · 104 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 102 simulations for each value of w*.

For any simulation (i.e. 102 times out of 102 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 70-th. 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 piece-wise constant (left) and impulsive (right) ACIs of Figure 7.

Table 6 Averages and standard deviation, daily ACI mono-therapy

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* · 104/(V A) for piece-wise constant therapy (left) and w k = w* · 104 (right). The densities are plotted in Figure 8, and the corresponding mean and standard deviation are shown in Table 7.

Figure 8
figure 8

Probability density function, weekly ACI mono-therapy. Empirical evaluation of ϱ(t er ) for different piece-wise constant (left) and impulsive (right) weekly ACIs. In (left) b i = w* · 7 · 104/(V A) for i = 1, . . . , 30 and in (right) w i = w* · 7 · 104 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 102 simulations for each value of w*.

Table 7 Averages and standard deviation, weekly ACI mono-therapy

By means of the Wilcoxon statistical test we compared the observed realizations of t er , grouped by kind of delivery (i.e. impulsive or piece-wise constant) and by frequency (i.e. daily or weekly). We obtained that (i) if delivered daily the eradication times of the piece-wise 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 piece-wise 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 IL-based 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 IL-based 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 IL-based immunotherapies the piece-wise 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 piece-wise 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 anti-tumor therapies such as anti-angiogenesis 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 piece-wise 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 tumor-immune 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.

Abbreviations

ACI:

Adoptive Cellular Immunotherapy. IL or IL-2: Interleukin-2. ODE: Ordinary Differential Equation. T-IS: Tumor-Immune system (interplay). KP: Kirschner-Panetta (model). SSA: Stochastic Simulation Algorithm.

References

  1. d'Onofrio A, Gandolfi A, Rocca A: The dynamics of tumor-vasculature interaction suggests low-dose, time-dense antiangiogenic schedulings. Cell Prolif 2009, 42: 317–329. 10.1111/j.1365-2184.2009.00595.x

    Article  PubMed  Google Scholar 

  2. 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.141135

    Article  CAS  PubMed  Google Scholar 

  3. Ehrlich P: Ueber den jetzigen Stand der Karzinomforschung. Ned Tijdschr Geneeskd 1909, 5: 273–290.

    Google Scholar 

  4. 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.104803

    Article  CAS  PubMed  Google Scholar 

  5. 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.006

    Article  Google Scholar 

  6. Whiteside TL: Tumor-induced death of immune cells: its mechanisms and consequences. Semin Cancer Biol 2002, 12: 43–50. 10.1006/scbi.2001.0402

    Article  CAS  PubMed  Google Scholar 

  7. 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.0400

    Article  CAS  PubMed  Google Scholar 

  8. 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.

    Article  PubMed  Google Scholar 

  9. 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 

  10. Agarwala SS: New applications of cancer immunotherapy. Semin Oncol 2002, 29(3 Suppl 7):1–4.

    Article  PubMed  Google Scholar 

  11. Bleumer I, Oosterwijk E, De Mulder P, Mulders PF: Immunotherapy for renal cell carcinoma. Eur Urol 2003, 44: 65–75,. 10.1016/S0302-2838(03)00191-X

    Article  CAS  PubMed  Google Scholar 

  12. Kaminski JM, Summers JB, Ward MB, Huber MR, Minev B: Immunotherapy and prostate cancer. Cancer Treat Rev 2003, 29: 199–209. 10.1016/S0305-7372(03)00005-7

    Article  CAS  PubMed  Google Scholar 

  13. De Pillis LG, Radunskaya AE, Wiseman CL: A Validated Mathematical Model of Cell-Mediated Immune Response to Tumor Growth. Cancer Res 2005, 65: 7950–7958.

    CAS  PubMed  Google Scholar 

  14. Kirschner D, Panetta JC: Modeling immunotherapy of the tumor - immune interaction. J Math Biol 1998, 37: 235–252. 10.1007/s002850050127

    Article  CAS  PubMed  Google Scholar 

  15. 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/s00262-007-0387-z

    Article  PubMed  Google Scholar 

  16. Cappuccio A, Elishmereni M, Agur Z: Cancer immunotherapy by interleukin-21: potential treatment strategies evaluated in a mathematical model. Cancer Res 2006, 66: 7293–7300. 10.1158/0008-5472.CAN-06-0241

    Article  CAS  PubMed  Google Scholar 

  17. Meloni G, et al.: How long can we give interleukin-2? Clinical and immunological evaluation of AML patients after 10 or more years of IL2 administration. Leukemia 2002, 16: 2016–2018. 10.1038/sj.leu.2402566

    Article  CAS  PubMed  Google Scholar 

  18. De Lisi C, Rescigno A: Immune surveillance and neoplasia: a minimal mathematical model. Bull Math Biol 1977, 39(2):201–221.

    CAS  Google Scholar 

  19. Stepanova NV : Course of the immune reaction during the development of a malignant tumor. Bioph 1980, 24: 917–923.

    Google Scholar 

  20. 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.

    Article  CAS  PubMed  Google Scholar 

  21. Nani F, Freedman HI: A mathematical model of cancer treatment by immunotherapy. Math Biosci 2000, 163: 159–199. 10.1016/S0025-5564(99)00058-9

    Article  CAS  PubMed  Google Scholar 

  22. Kuznetsov VA, Knott GD: Modeling tumor re-growth and immunotherapy. Math Comp Mod 2001, 33: 1275–1287. 10.1016/S0895-7177(00)00314-9

    Article  Google Scholar 

  23. 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.012

    Article  PubMed  Google Scholar 

  24. d'Onofrio A: A general framework for modeling tumor-immune system competition and immunotherapy: Mathematical analysis and biomedical inferences. Phys D 2005, 208: 220–235. 10.1016/j.physd.2005.06.032

    Article  Google Scholar 

  25. Kirschner D, Tsygvintsev A: On the global dynamics of a model for tumor immunotherapy. Math Biosci Eng 2009, 6(3):573–583.

    Article  PubMed  Google Scholar 

  26. d'Onofrio A: Tumor-immune system interaction: modeling the tumor-stimulated proliferation of effectors and immunotherapy. Math Mod and Meth in App Sci 2006, 16: 1375–1401. 10.1142/S0218202506001571

    Article  Google Scholar 

  27. d'Onofrio A, Gatti F, Cerrai P: Delay-induced oscillatory dynamics of tumor-immune system interaction. Math and Comp Mod 2010, 51: 572–591. 10.1016/j.mcm.2009.11.005

    Article  Google Scholar 

  28. d'Onofrio A: Bounded-noise-induced transitions in a tumor-immune system interplay. Phys Rev E Stat Nonlin Soft Matter Phys 2010, 81: 021923.

    Article  PubMed  Google Scholar 

  29. 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.001

    Article  Google Scholar 

  30. 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.

    Article  PubMed  Google Scholar 

  31. Pappalardo F, Lollini PL, Castiglione F, Motta S: Modeling and simulation of cancer immunoprevention vaccine. Bioinformatics 2005, 21(12):2891–2897. 10.1093/bioinformatics/bti426

    Article  CAS  PubMed  Google Scholar 

  32. 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/1471-2105-11-S7-S13

    Article  PubMed Central  PubMed  Google Scholar 

  33. Kennedy BJ: . Cyclic leukocyte oscillations in chronic myelogenous leukemia during hydroxyurea therapy. Blood 1970, 35: 751–760.

    CAS  PubMed  Google Scholar 

  34. 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/NEJM197202102860603

    Article  CAS  PubMed  Google Scholar 

  35. Gatti R, et al.: Cyclic Leukocytosis in Chronic Myelogenous Leukemia: New Perspectives on Pathogenesis and Therapy. Blood 1973, 41: 771–783.

    CAS  PubMed  Google Scholar 

  36. Mehta BC, Agarwal MB: Cyclic oscillations in leukocyte count in chronic myeloid leukemia. Acta Haematol 1980, 63: 68–70. 10.1159/000207373

    Article  CAS  PubMed  Google Scholar 

  37. 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.03310030040023

    Article  CAS  PubMed  Google Scholar 

  38. Tsao H, Cosimi AB, Sober AJ: Ultra-late recurrence (15 years or longer) of cutaneous melanoma. Cancer 1997, 79(12):2361–2370. 10.1002/(SICI)1097-0142(19970615)79:12<2361::AID-CNCR10>3.0.CO;2-P

    Article  CAS  PubMed  Google Scholar 

  39. 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.013

    Article  CAS  PubMed  Google Scholar 

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

    Chapter  Google Scholar 

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

    Chapter  Google Scholar 

  42. 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.

    Article  Google Scholar 

  43. Gillespie DT: Exact Stochastic Simulation of Coupled Chemical Reactions. J of Phys Chem 1977, 81: 2340–2361. 10.1021/j100540a008

    Article  CAS  Google Scholar 

  44. 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.1835951

    Article  Google Scholar 

  45. 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 

  46. Lecca P: A time-dependent 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 

  47. 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.2799998

    Article  PubMed  Google Scholar 

  48. Kirschner D, Arciero JC, Jackson TL: A Mathematical Model of Tumor-Immune Evasion and siRNA Treatment. Discr and Cont Dyn Systems 2004, 4: 39–58.

    Article  Google Scholar 

  49. DeBoer RJ, Hogeweg P, Hub F, Dullens J, DeWeger RA, DenOtter W: Macrophage T Lymphocyte interactions in the anti-tumor immune response: A mathematical model. J Immunol 1985, 134: 2748–2758.

    CAS  Google Scholar 

  50. d'Onofrio A: Stability property of Pulse Vaccination Strategy in SEIR epidemic model. Math Biosci 2002, 179(1):57–72. 10.1016/S0025-5564(02)00095-0

    Article  PubMed  Google Scholar 

  51. 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/S0895-7177(02)00176-0

    Article  Google Scholar 

  52. Shahrezaei V, Ollivier JF, Swain PS: Colored extrinsic fluctuations and stochastic gene expression. Mol Syst Biol 2008, 4: 196.

    Article  PubMed Central  PubMed  Google Scholar 

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

    Google Scholar 

  54. 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.0020117

    Article  PubMed Central  PubMed  Google Scholar 

  55. Bratsun D, Volfson D, Tsimring LS, Hasty J: Delay-induced Stochastic Oscillations in Gene Regulation. PNAS 2005, 102(41):14593–14598. 10.1073/pnas.0503858102

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  56. Barbuti R, Caravagna G, Maggiolo-Schettini A, Mi-lazzo 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 

  57. Barbuti R, Caravagna G, Maggiolo-Schettini A, Milazzo P: Delay Stochastic Simulation of Biological Systems: A Purely Delayed Approach. Trans on Comp Sys Bio 2011, XIII 6575: 61–84.

    Article  Google Scholar 

  58. Levi F, Schibler U: Circadian Rhythms: Mechanisms and Therapeutic Implications. Annu Rev Pharmacol Toxicol 2007, 47: 593–628. 10.1146/annurev.pharmtox.47.120505.105208

    Article  CAS  PubMed  Google Scholar 

  59. Kerbel RS, Kamen BA: The anti-angiogenic basis of metronomic chemotherapy. Nat Rev Cancer 2004, 4: 423–436. 10.1038/nrc1369

    Article  CAS  PubMed  Google Scholar 

  60. 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.fc

    Article  CAS  PubMed  Google Scholar 

  61. Gillespie DT: Approximated Accelerated Stochastic Simulation of Chemically Reacting Systems. J of Chem Phys 2001, 115(4):1716–1733. 10.1063/1.1378322

    Article  CAS  Google Scholar 

  62. 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 

  63. d'Onofrio A, Gandolfi A: Resistance to antitumor chemotherapy due to bounded-noise-induced transitions. Phys Rev E Stat Nonlin Soft Matter Phys 2010, 82: 061901.

    Article  PubMed  Google Scholar 

  64. 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 

Download references

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 "P-medicine - 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/1471-2105/13/S4.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Alberto d'Onofrio.

Additional information

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.

Electronic supplementary material

12859_2012_5105_MOESM1_ESM.pdf

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. (PDF 200 KB)

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Caravagna, G., Barbuti, R. & d'Onofrio, A. Fine-tuning anti-tumor immunotherapies via stochastic simulations. BMC Bioinformatics 13 (Suppl 4), S8 (2012). https://doi.org/10.1186/1471-2105-13-S4-S8

Download citation

  • Published:

  • DOI: https://doi.org/10.1186/1471-2105-13-S4-S8

Keywords