Martin Zumsande
email: zumsande@mpipks-dresden.mpg.de

Thilo Gross

Bifurcations and Chaos in the MAPK signaling cascade

Keywords: Modelling, Signalling networks, Phosphorylation, Bistability

The mitogen-activated protein kinase (MAPK) cascade is an important signaling cascade in eukaryotes. We use the approach of generalized modeling to analyze the dynamics of the MAPK cascade and identify key mechanisms of instability. Furthermore, we report sustained multi-mode oscillations and potentially chaotic behavior caused by a sequestration-based feedback mechanism. Finally, we investigate the interplay between sequestration and external feedback loops. Our analysis thereby confirms, extends and generalizes previous results obtained by conventional modeling and points out the diversity of dynamics that sequestration can bring about.

Introduction

The mitogen-activated protein kinase (MAPK) pathway is involved in the regulation of numerous cell functions, among them proliferation, apoptosis, differentiation and cell motility [START_REF] Seger | The MAPK signaling cascade[END_REF][START_REF] Chang | Mammalian MAP kinase signalling cascades[END_REF][START_REF] Johnson | Mitogen-Activated Protein Kinase Pathways Mediated by ERK, JNK, and p38 Protein Kinases[END_REF]. Moreover, in many types of cancer, dysfunctions in the MAPK pathway have been detected [START_REF] Downward | Targeting RAS signalling pathways in cancer therapy[END_REF][START_REF] Roberts | Targeting the Raf-MEK-ERK mitogenactivated protein kinase cascade for the treatment of cancer[END_REF]. Therefore, a heightened understanding of the pathway may potentially lead to the discovery of new drug targets for cancer treatment.

While there are several variants of MAPK cascades depending on organism and signal type, the topological motif shown in Fig. 1 is conserved in MAPK pathways throughout all eukaryotic cells [START_REF] Seger | The MAPK signaling cascade[END_REF]. In the figure, we show a prominent example of a MAPK cascade, the growth-factor induced Ras-Raf-MEK-ERK cascade. When a cell is exposed to certain growth factors, a series of intermediate reactions is triggered that ultimately lead to the activation of the enzyme Ras, which itself phosphorylates Raf, a protein also called MAPKKK (MAP-Kinase-Kinase-Kinase). In its activated state, Raf can phosphorylate the protein MEK (MAPKK) at two different sites. The double-phosphorylated MEK can then act as an enzyme in the double phosphorylation of ERK (MAPK), which activates several transcription factors and downstream kinases in its phosphorylated state. Metabolic Control Analysis (MCA) suggests that these reactions downstream from Ras, the actual cascade, are most important for the dynamics [START_REF] Hornberg | Control of MAPK signalling: from complexity to what really matters[END_REF].

The first mathematical model of the MAPK cascade was formulated by [START_REF] Huang | Ultrasensitivity in the mitogenactivated protein kinase cascade[END_REF], showing that the system can act as an ultrasensitive switch. Later, it was found that adding explicite negative feedback leads to oscillations [START_REF] Kholodenko | Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades[END_REF] and that the combination of double phosphorylation and enzyme sequestration introduces an implicit feedback leading to bistability [START_REF] Markevich | Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades[END_REF] and oscillations [START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF]. Most investigations adopted the ranges of model parameters (concentrations and rate constants) from [START_REF] Huang | Ultrasensitivity in the mitogenactivated protein kinase cascade[END_REF], based on measurements in Xenopus oocyte extracts. However, little is known about the possible parameter ranges in which eukaryotic cells can operate in different species and at different developmental stages. It has been argued before that the complicated structure of a three-layered double-phosphorylation cascade could have evolved as a result of an optimization of the response time to signals [START_REF] Frey | How quantitative measures unravel design principles in multi-stage phosphorylation cascades[END_REF]. However, there may also be adaptive advantages for the realization chosen by nature that are based on other dynamical properties of the underlying models.

In this paper we apply the approach of generalized modeling [START_REF] Gross | Generalized models as an universal approach to the analysis of nonlinear dynamical systems[END_REF][START_REF] Steuer | Structural kinetic modeling of metabolic networks[END_REF][START_REF] Gross | Generalized models reveal stabilizing factors in food webs[END_REF] to the MAPK cascade. The main advantage of this approach is that it avoids simulation,

A c c e p t e d m a n u s c r i p t

which is numerically demanding and can be unreliable for some systems. While generalized models reveal only certain types of information, this information is gained very efficiently, allowing us to explore a large parameter space.

Our analysis identifies the key parameters that govern the stability of stationary states of the MAPK cascade.

We then analyze the transitions in which the stability is lost. This reveals regions of oscillatory behavior and provides evidence for bistability, generalizing previous results. Then we apply the generalized analysis to identify a parameter region in which complex dynamics can be observed. Finally, we investigate the effect of external feedback mechanisms on the dynamics of the cascade.

Formulation of the model

Generalized modeling

The example of the MAPK signaling pathway presents a typical challenge: While the structure of the underlying network of protein interactions is relatively well known, the functioning of the system is closely linked to its dynamics and therefore difficult to observe in experiments. Conventional modeling approaches aiming at exploring the dynamics of the system are forced to restrict every process to a specific rate law. It is therefore often not clear whether the observed dynamics represents robust behavior of the system or arises only because of the specific choice of functional forms. Moreover, the exploration of the long-term behavior of conventional models by direct simulation is numerically inefficient. Computational restraints therefore impose strong limits on the size of the parameter space that can be explored.

The approach of generalized modeling is based on the insight that certain information can be extracted more efficiently by considering the whole class of plausible models [START_REF] Gross | Generalized models as an universal approach to the analysis of nonlinear dynamical systems[END_REF][START_REF] Steuer | Structural kinetic modeling of metabolic networks[END_REF][START_REF] Gross | Generalized models reveal stabilizing factors in food webs[END_REF]. Instead of parameterizing the specific rate laws in the model, we directly parameterize the Jacobian matrix, governing the dynamics close to all possible steady states in all models that are consistent with the observed interaction network. Only subsequently is the class of models under consideration narrowed down by fixing certain parameters of the generalized model.

Before we focus specifically on the MAPK cascade, let us introduce the approach of generalized modeling as such. For this purpose consider a given network of biochemical interactions regulating the concentrations of a number of proteins. If the concentrations are sufficiently high to neglect stochastic effects, the dynamics of the system can be captured by a system of differential equations of the form

d dt S(t) = N f (S), ( 1
)
where S is a vector of protein concentrations, N is the stoichiometric matrix and f is a flux vector depending on the S i .

Assuming that a positive, but not necessarily stable, steady state S * exists, this equation can be normalized with respect to the steady state concentrations and written in the form

d dt x = Λμ(x) ( 2 )
with

x i = S i /S * i , Λ i j = N i j f j (S * )/S * i and μ i (x) = f i (S)/f i (S * ).
The local dynamics close to the steady state S * is governed by the corresponding Jacobian matrix defined as J i j = ∂ S i (t)/∂ S j S * . Since Λ does not depend on x the Jacobian of our system can be written as

J = Λθ μ x (3) with θ μ j x i = ∂μ j (x) ∂x i S * . ( 4
)
The steady state under consideration is stable if all eigenvalues of the Jacobian have negative real parts and unstable if at least one eigenvalue has a positive real part.

Changing the parameters can bring about a loss of stability when it causes a single real eigenvalue or a complex conjugate pair of eigenvalues to cross the imaginary axis of the complex plane, acquiring a positive real part. The former case corresponds to a saddlenode bifurcation in which the number of steady states changes, while the latter corresponds to a Hopf bifurcation in which a limit cycle emerges.

The entries of Λ and θ μ

x constitute a complete parameterization for the entire range of possible Jacobians that are consistent with the underlying topology of biochemical interactions. Since we have not specified the righthand-side of the equations of motion, Eq. 2, we cannot compute the steady states of the system with the chosen degree of generality. Therefore, also the entries of the Jacobian matrix, which depend on the steady state, remain unknown. However, the quantities appearing in the Jacobian are constants and can therefore be interpreted as unknown parameters with the same right as In general, the generalized parameters can be related to parameters appearing in conventional models [START_REF] Van Voorn | Stabilization due to predator interference: comparison of different analysis approaches[END_REF][START_REF] Stiefs | Evidence of chaos in ecoepidemic models[END_REF]. Given a steady state from a conventional model, it is straightforward to compute the corresponding generalized parameters. Conversely, a single set of generalized parameters can be translated to a family of conventional models, exhibiting a corresponding steady state. While the difficulty of the latter depends on the functional forms that are assumed in the conventional model, it is no major obstacle in the models with mass action used in this paper.

The parameters identified by the generalized modeling procedure have an intuitive biological interpretation:

The entries of Λ have the dimension of an inverse time and represent characteristic time scales of the model. They are therefore called scale parameters in generalized modeling. The entries of θ μ

x are logarithmic derivatives of the fluxes. The corresponding parameters are called exponent parameters in the context of generalized models and elasticities in the context of metabolic control theory [START_REF] Fell | Metabolic control analysis: a survey of its theoretical and experimental developement[END_REF]). Note, for example that if one of the processes in the model is linear f i ∼ S j the corresponding exponent parameter is θ i j = 1, irrespective of the constant of proportionality. More generally, for any power-law relationship f i ∼ S j p the corresponding parameter is the exponent θ i j = p. Even for functions that are not power laws, the exponent parameters can be interpreted as a measure of the local nonlinearity.

In the steady state under consideration all fluxes in the model have to balance N f (S * ) = 0, which constitutes a constraint for the permissible values of scale parameters. To satisfy this constraint and reduce the number of free parameters it can be advantageous to use an alternative parameterization based on the flux modes of the model [START_REF] Steuer | From structure to dynamics of metabolic pathways: Application to the plant mitochondrial TCA cycle[END_REF]. In this parameterization the scale parameters denote the strength of the independent flux modes that are consistent with the stationarity condition.

In summary, the generalized model provides a parameterization of the Jacobian in all steady states that can be found in a class of models defined by a given interaction topology. While the analysis of the Jacobian can only reveal certain types of information it does so with a high degree of robustness and numerical efficiency. Although numerical simulation of conventional models is useful to supplement insights gained from generalized modeling, the generalized model allows us to explore dynamics close to steady states without resorting to numerical simulation. The analysis of the Jacobian matrix benefits from the availability of highly efficient algorithms for the computation of eigenvalues. Therefore, in models of intermediate complexity, say less than 50 dynamical variables, several tens of billions parameter sets can be explored in reasonable computational time [START_REF] Gross | Generalized models reveal stabilizing factors in food webs[END_REF]. By comparison, for the present model already the analysis of just tens of millions of parameter sets yields sound statistical results.

Model of the MAPK cascade

In our model of the MAPK cascade we use the following naming convention: The symbol S i j is used to denote the concentrations of proteins and protein complexes which appear as substrates in the model. The superscript index i indicates the layer of the cascade in which the specific protein is involved, while the subscript index j enumerates all proteins in the respective layer. Note that, although we avoid the use of asterisks, in the following all symbols S i j denote concentrations in the steady state.

The concentrations of the additional free kinases and phosphatases that do not appear as substrates themselves are labeled K i R and P i R , respectively. The total concentrations of the substrates, kinases and phosphatases in a given layer of the cascade are denoted S i T , K i T and P i T . As these concentrations are conserved in the context of the model, they constitute constant external parameters. These naming conventions are illustrated in Fig. 2.

In the following we consider three different models describing subsystems of the MAPK cascade of different complexity. The first of these models contains only a double phosphorylation step as it occurs in the second layer of the cascade whereas the first and third layer are neglected. The second model contains the first and the second layer of the cascade, but neglects the third. Finally, the third model contains all three layers. The scope of the different models is indicated by boxes in Fig. 2. The kinase that catalyzes the phosphorylation in the second layer is called K 1 R in the first model in order to indicate that, in the scope of this model, it is not a substrate by itself. In the two larger models it is called S 0 1 as it appears as a substrate in the first layer.

For the dynamics of the cascade, enzyme sequestration is suspected to be of central importance [START_REF] Blüthgen | Effects of sequestration on signal transduction cascades[END_REF][START_REF] Legewie | Competing docking interactions can bring about bistability in the MAPK cascade[END_REF][START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF]. Sequestration can affect the dynamics because an enzyme that is bound to a substrate cannot at the same time participate in other reactions, including its own dephosphorylation. Without external feedback, sequestration effects provide the only means by which the concentrations on the lower layers of the cascade can affect the

A c c e p t e d m a n u s c r i p t

dynamics on the layers above. To capture the sequestration effects explicitly we follow [START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF] and break each phosphorylation step into two separate processes, for example

S 0 0 + K 0 R k 0 k 1 S 0 2 k 2 → S 0 1 + K 0 R , (5)
The first step represents the reversible binding of the kinase K to the unphosphorylated substrate S 0 0 , while the second step represents the action of the kinase resulting in the release of the phosphorylated substrate S 0 1 . Other approaches such as Michaelis-Menten based enzyme kinetics can miss sequestration effects. Although a sequestration-based approximation scheme has been proposed recently [START_REF] Ventura | A hidden feedback in signaling cascades is revealed[END_REF], the efficiency of generalized modeling does not rely on the number of reactions or variables being small. We therefore avoid this approximation and represent the complexes explicitly in the model.

Assuming that the concentrations of proteins are sufficiently high to be considered as continuous variables, we capture the dynamics of the cascade by a system of coupled ODEs. For the phosphorylation step in the top layer of the cascade we obtain

d dt S 0 0 = -f 1 (S 0 0 , K 0 R ) + f 2 (S 0 2 ) d dt S 0 2 = f 1 (S 0 0 , K 0 R ) -f 2 (S 0 2 ) d dt S 0 1 = f 3 (S 0 2 ) (6)
where f 1 , f 2 , and f 3 are general functions, which are not yet restricted to specific functional forms.

Starting from the set of differential equations, the Jacobian of the generalized model is obtained by following the steps outlined in the previous section. While we omit this derivation here the resulting equations are shown in the supplementary material. Also, the complete derivation of the generalized model for the onelayer model is shown for the purpose of illustration in the appendix.

As the model explicitly describes the complexes formed by enzymes with their substrates, it is reasonable to assume that the remaining processes are governed by mass action. We therefore utilize the ability of generalized models to deal with general functions only when considering the effect of external feedbacks, while the ability to explore large parameter spaces is used throughout the whole study. With this choice our model corresponds directly to the conventional models considered in [START_REF] Huang | Ultrasensitivity in the mitogenactivated protein kinase cascade[END_REF][START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF].

Note that even if it is assumed that the processes are linear, the system of equations is nonlinear as it includes bilinear terms. Using mass conservation to express the concentrations of certain proteins as a function of others (see appendix) changes some of the bilinear terms into explicit quadratic nonlinearities.

Investigation of the dynamics

Stability analysis of the single-layer model

Following [START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF], we started by investigating the second layer of the cascade in isolation. This comparably simple system is still analytically tractable and has been investigated in earlier work [START_REF] Markevich | Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades[END_REF][START_REF] Ortega | Bistability from double phosphorylation in signal transduction. Kinetic and structural requirements[END_REF].

To identify the parameters that have a strong impact on the stability of steady states we created a large ensemble of 10 6 parameter sets where the parameters in each set were generated as follows: We first drew the concentrations of S 1 3 , S 1 5 and K 1 R from a uniform distribution. These concentrations were then normalized such that

K 1 T = S 1 3 + S 1 5 + K 1 R = 1. Then, S 1 4 , S 1 6
and P 1 R were drawn randomly and normalized such that

P 1 T = S 1 4 + S 1 6 + P 1 R = 1.
To ensure that sequestration effects can have an impact on the dynamics we set the total concentration S 1 T of the substrate to 10, so that there is considerably more substrate than kinases and phosphatases. A fraction of the substrate is already bound in the complexes S 1 3 , S 1 5 , S 1 4 , and S 1 6 . We distribute the remainder r = S 1 T -S 1 3 -S 1 5 -S 1 4 -S 1 6 by drawing random concentrations for S 1 0 , S 1 1 , and S 1 2 and normalizing so that S 1 0 + S 1 1 + S 1 2 = r. Finally we drew the strength of all flux modes in the model randomly from a uniform distribution.

In order to assess the impact of the individual parameters on stability we computed the stability of the steady state in each sample parameter set. Stable steady states were assigned the stability value 1 whereas unstable states were assigned the stability value 0. For each parameter we then computed the Pearson correlation coefficient of the focal parameter with the stability value across the whole ensemble.

For the single layer the correlation coefficients are shown in Fig. 3A. The concentrations of the unbound kinase K 1 R and phosphatase P 1 R show a correlation coefficient of 0.21 with stability. Consequently the system is likely to be stable if these concentrations are high and unstable if they are low. This shows that instability occurs most likely if a large fraction of the kinase and phosphatase is sequestered in complexes. Particularly, the complexes S 1 3 and S 1 6 , which are involved in the production of the single-phosphorylated S 1 1 are negatively

A c c e p t e d m a n u s c r i p t

correlated with stability, showing that a high concentration of these complexes has a strong destabilizing effect. By contrast, the complexes S 1 4 and S 1 5 which appear in reactions decreasing S 1 1 are positively correlated with stability.

More detailed insights in the effect of important parameters can be gained by plotting the fraction p of randomly drawn steady states that are stable. For this purpose we generated another ensemble of 10 7 parameter sets as described above. Figure 4 shows p as a function of selected combinations of two parameter values. In Fig. 4A, p is shown as a function of the free enzyme concentrations K 1 R and P 1 R . The figure reveals that low concentrations of either K 1 R and P 1 R promote instability. This confirms that sequestration plays a role in the destabilization of steady states. The dependence of p on the parameters S 1 3 and S 1 6 , which were found to have a strongly destabilizing impact, is plotted in Fig. 4B. The figure shows that the fraction of stable states changes very sharply as a threshold is crossed. By comparison, even if the stabilizing parameters S 1 4 and S 1 5 are very low p only drops to 0.5 (Fig. 4C). The comparison confirms the lesser impact of the two stabilizing parameters, already seen in the correlation analysis. This conclusion is further supported by the direct comparison of the destabilizing parameter S 1 3 with the stabilizing parameter S 1 4 (Fig. 4D). In summary we conclude that the main source of instability in the single-layer system lies in high concentrations of the complexes S 1 3 and S 1 6 .

Mechanisms of instability

To understand the mechanism leading to instability in more detail we consider the eigenvectors of the system's Jacobian. In an unstable steady state the eigenvector corresponding to largest (positive) eigenvalue represents the direction of perturbation in the state space of the system that grows most rapidly in time. Therefore, the eigenvector also marks the direction in which the system leaves the steady state after destabilization occurs. We computed the typical direction of escape from the unstable steady state by averaging over the eigenvectors of the largest eigenvalue in 10 6 unstable parameter sets generated randomly. The components of this averaged direction of escape are shown symbolically in Fig. 5. The figure illustrates that the instability causes the proteins to shift into either one of the phosphorylation cycles, while depleting the other.

The nature of the instability described above can be linked to the destabilizing effect of S 1 3 and S 1 6 . If there is a large concentration of S 1 3 (or, S 1 6 ) then the kinase (phosphatase) is not available for the formation of the complex S 1 5 (S 1 4 ). This implies the primary reaction involving S 1 1 as a substrate is binding in the complex S 1 4 (S 1 5 ). This leads to accumulation in the left (right) phosphorylation loop and consequently increases S 1 3 (S 1 6 ) further. Thus a positive feedback loop is formed that causes the instability.

Relation to bistability

Unstable steady states are often located on a separatrix that divides the basins of attraction of different attractors. Also in the MAPK cascade previous works have shown that a single layer of the cascade can support bistable dynamics [START_REF] Markevich | Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades[END_REF][START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF], in which the system can approach either one of two stable steady states. To confirm that the unstable steady states observed in the generalized model mark the separatrix between two stable steady states, we simulated the dynamics of the conventional model used in [START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF] using a stiff numerical integrator based on the modified extended backward differentiation formula (MEBDF) [START_REF] Hairer | Solving ordinary differential equations II: Stiff and differential-algebraic problems[END_REF]. We randomly selected 100 starting points corresponding to unstable states computed in the generalized model. For all of these points we found that multiple trajectories starting close to the steady state eventually converged to two different final states. An example time series is shown in the supplementary material. This result suggests that, in a class of reasonable conventional models of the second layer of the MAPK cascade, the longterm dynamics is bistable whenever we find an unstable steady state in the generalized model.

Stability in larger models

As a next step we considered the stability of two larger subsystems of the cascade. We proceeded analogously to the investigation of the single-layer model by first generating a large ensemble of randomly drawn parameter sets and then correlating the stability of the corresponding steady states with the parameter values.

In the investigation of the larger models, describing two and three layers of the cascade, we slightly altered our procedure for generating random parameter sets in order to generate sets that are close to empirical observations, using data from [START_REF] Huang | Ultrasensitivity in the mitogenactivated protein kinase cascade[END_REF]. The details of this procedure are explained in the appendix.

The results of a correlation analysis are shown in Fig. 3B,C. Figure 3B shows the stability correlations in our second model, which consists of the first and second layer of the cascade. While there are now more parameters than in model 1 (Fig. 3A), the additional parameters, describing the top layer of the cascade have

A c c e p t e d m a n u s c r i p t

a weaker impact on the stability of steady states, as evidenced by the low values of the corresponding correlation coefficients. Note furthermore that the pattern of correlation coefficients in the second model is very similar to that found in the first model. In particular, there is a repeated motif of strong positive correlations in P 1 R , S 1 4 , and S 1 5 and strong negative correlations in S 1 3 and S 1 6 .

In the third model, consisting of all three layers of the cascade, the pattern of correlations observed above reappears twice. As shown in Fig. 3C, the correlation coefficients corresponding to both the first and second layer of the cascade are now quite small, indicating a much reduced impact on the dynamics. Nevertheless the pattern of correlations is still visible in the second layer. More importantly, the pattern reappears in a stronger form in the third layer. Thus, the correlation analysis shows that the basic mechanism of instability, the positive feedback induced by sequestration of substrates into the complexes S 3 and S 6 , remains the same in all three models.

Oscillatory instabilities

In the next step, we focused on the dynamics that can be observed in the two and three layer systems after the stability of a steady state has be lost. For this purpose we searched for the bifurcation points, marking a change in the stability of a steady state. Since the probability of finding a bifurcation point by random sampling is of measure zero, we generated a ensemble of 10 7 random parameter sets as described above, but selected only those parameter sets in which the real part of the leading eigenvalue of the Jacobian in the steady state was in the range [0,0.01], indicating proximity to a bifurcation point. Those parameter sets were then settled to the bifurcation manifold by a Newton fixed-point algorithm [START_REF] Kelley | Solving nonlinear equations with Newton's method[END_REF].

In the single-layer model this procedure yields only bifurcation points of saddle-node type, which in this model mark the onset of bistability. The absence of Hopf bifurcations, which typically correspond to the onset of oscillatory dynamics, supports the previously observed absence of oscillations in the single-layer model. In the two and three-layer models we find both Hopf and saddle-node bifurcations, which shows that both bistability and oscillations are generically possible in these models.

Further insights can be gained by investigating the dynamics close to unstable steady states. For all unstable steady states found in 10 7 samples we checked whether the leading eigenvalue, which governs the departure from the steady state, is real or complex valued.

In the one-layer model, we found no unstable steady states with complex leading eigenvalues, providing further evidence for the absence of oscillatory dynamics. By contrast, leading complex eigenvalues are found in 4% of the unstable steady states of the two-layer system and in 34% of the unstable steady states of the threelayer system. This suggests that oscillatory dynamics appear more prominently in the three-layer cascade.

Complex global dynamics

Since we saw that oscillatory dynamics is absent in the one-layer system, but appears when a second layer is added, one can ask whether even more complex dynamics such as quasiperiodicity or chaos can occur in the three-layer system. Generalized modeling, being essentially based on a local analysis, cannot detect such global dynamics directly. However, generalized models can be used to search for local bifurcations of higher codimension that can serve as proxies for complex dynamics.

In order to find a bifurcation of codimension two, one must adjust two parameters of the system simultaneously to their exact bifurcation values. Therefore, codimension-2 bifurcations are almost never observed directly in experiments and have, for the lack of examples, received relatively little attention in the mathematical literature (Guckenheimer and Holmes, 1997)(p. 420).

Nevertheless the identification of codimension-2 bifurcations can reveal many insights in the dynamics of a system because of their role of organizing centers for different types of dynamics.

An example of a codimension-2 bifurcation is shown in Fig. 6 in a three-parameter bifurcation diagram, which was created by a triangulation procedure described in [START_REF] Stiefs | Computation and visualization of bifurcation surfaces[END_REF]. The parameters are given in Table S1 in the supplementary material. In the threedimensional parameter space a surface of Hopf bifurcation points (light surface) and a surface of saddle-node bifurcation points (dark surface) intersect. On the intersection line of the surfaces the Jacobian has both a single zero eigenvalue and a purely imaginary eigenvalue pair, marking the points on the intersection line as codimension-2 Gavrilov-Guckenheimer bifurcation points. It is well known from normal form analysis that close to this type of bifurcation quasiperiodic dynamics should generally occur (Kuznetsov, 1995)(p. 345). Another example of a codimension-2 bifurcation is the double-Hopf bifurcation, which is characterized by the presence of two purely imaginary eigenvalue pairs. This bifurcation generically involves the creation of chaotic and quasiperiodic regions (Kuznetsov, 1995)(p. 369).

A c c e p t e d m a n u s c r i p t

While the presence of Gavrilov-Guckenheimer and double Hopf bifurcations cannot guarantee complex long-term behavior, they serve as proxies indicating parameter regions in which such complex dynamics is likely. We therefore used the generalized model to search specifically for these bifurcations. Explicit numerical simulation of the conventional model was then used to explore the dynamics close to these bifurcation points. In this way the greater efficiency of generalized models is combined with the higher predictive power of conventional models. In the model of the MAPK cascade we observed interesting dynamics only close to double-Hopf bifurcations, which appear exclusively in the three-layer system. A two-parameter bifurcation diagram of a double-Hopf bifurcation and two example time series from different regions are shown in Fig. 7. The corresponding parameters are listed in Table S1 and Table S2 in the supplementary material. Both timeseries show evidence of complex mixed-mode oscillations. Moreover the spike-like time-series shown in Fig. 7C exhibits an irregularity that is indicative of Šilnikov chaos. Similar dynamics were also observed at multiple points in the neighborhood (not shown). We are therefore confident that the system shows chaotic long-term behavior in a finite parameter region. However, to confirm the chaotic nature of the dynamics and to determine the parameter ranges in which they occur will require more extensive numerical investigations, including the computation of Lyapunov exponents, which exceed the scope of this paper.

Explicit Feedback

In addition to the sequestration-mediated negative feedback inherent in the cascade, we also considered the implications of adding explicit feedback loops. While the structure of the MAPK cascade is widely conserved throughout eukaryotic cells, this is not necessarily true for the types of feedback that may depend on species and tissue [START_REF] Kolch | When kinases meet mathematics: the systems biology of MAPK signalling[END_REF]. Here we consider two types of mechanisms: a) In fibroblasts, inhibitory phosphorylation of SOS, a protein upstream from Ras (K 0 R ), by ERK (S 2

2 ) effectively inhibits the formation of the Ras-Raf complex (S 0 2 ) [START_REF] Langlois | Negative Feedback Regulation and Desensitization of Insulin-and Epidermal Growth Factor-stimulated p21ras Activation[END_REF][START_REF] Kholodenko | Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades[END_REF]. b) In COS1-cells, activated ERK (S 2

2 ) feeds back to MEK (S 1 0 ) during cell adhesion by phosphorylating it at the site T292, which inhibits phosphorylation by PAK at the adjacent site S298. The suppressed phosphorylation facilitates the formation ERK-Raf complexes (S 1 3 ,S 1 5 ) [START_REF] Eblen | Mitogen-activated protein kinase feedback phosphorylation regulates MEK1 complex formation and activation during cellular adhesion[END_REF], resulting in a negative feedback of S 2 2 on S 1 0 . As it is widely believed that the MAPK cascade is utilized mainly as an ultrasensitive switch, we ask whether the feedback loops act to suppress other dynamics, such as oscillations. We therefore investigated whether the oscillatory regime is reduced by introducing feedback in the system. To study this question in conventional models is not trivial for two reasons: First, in a conventional model the effect of the feedback has to be restricted to a specific functional form, which can in practice be difficult to derive. Second, and perhaps more importantly, to study the effect of feedbacks in isolation in conventional models is not trivial. The question is typically how the stability of a given steady state observed in nature depends on the nonlinearity of the feedback. However, if feedback parameters in conventional models are changed one generally observes the combined effect of both the altered nonlinearity and the corresponding shift in the steady state. Both of these difficulties are easily avoided in generalized models as they are designed to deal with unknown functions and parameterize the position of the steady state independently from the nonlinearities in the system.

We studied the effect of the feedback of type a) and b) by adding an inhibitory link from the final product of the cascade, S 2 2 , back to the phosphorylation steps in the first (a) and second (b) layer, respectively. In our model inhibitory feedback implies that the first derivative of the feedback function is negative. The effect of the feedback can therefore be modeled by introducing

θ f = θ μ 1 S 2 2 < 0 (feedback a) and θ f = θ μ 7 S 2 2 = θ μ 13 S 2 2 < 0
(feedback b). In order to determine the effect of the feedback we generated an ensemble consisting of 10 7 parameter sets for each of the two feedback mechanisms. The parameter sets were drawn randomly, as described above, except for the additional feedback parameter θ f , which were drawn randomly from a uniform distribution in the interval [-2, 0]. This range is for instance consistent with inhibition modeled by a Hill function with exponent two. For each parameter set, we then computed the spectrum of the Jacobian with and without feedback. While the stability of steady states corresponding to individual data sets changes when the feedback is switched on, the net effect on stability on the whole ensemble is small for both types of feedback (not shown). A stronger effect of the feedback becomes apparent if one considers the dynamics subsequent to the loss of stability. Figure 8 shows the change in the number of unstable states in which the leading eigenvalue of the Jacobian is complex or real, respectively. For both feedback mechanisms we observe that as the feedback parameter increases the number of unstable steady states with complex leading eigenvalues increases, while the number of unstable steady states

A c c e p t e d m a n u s c r i p t

with real leading eigenvalues decreases. Although the observation of an unstable steady state with complex leading eigenvalues does not imply that the long-term dynamics of the system is oscillatory, the increase of such states over a wide parameter range strongly suggests that also oscillatory long-term dynamics is promoted rather than suppressed by the two feedback loops under consideration.

Discussion

In this paper we used generalized modeling to analyze the dynamics of a class of models of different subsystems describing the MAPK cascade. By a correlation analysis, building on tens of millions of parameter sets, we determined the impact of the individual model parameters on the stability of the steady states. This analysis confirmed that sequestration of enzymes has a strong impact on the dynamics of the cascade. In particular, we showed that instability is likely if large portions of the kinase and phosphatase acting on the lowest level of the cascade are sequestered into complexes (S 2 3 and S 2 6 in our notation). In this case a positive feedback loop is formed that destabilizes stationary states.

In a second step we used bifurcation analysis, spectral analysis of unstable states, and explicit simulation to investigate the dynamics subsequent to the loss of stability. In a subsystem, consisting of only one layer of the cascade, bifurcation analysis of the generalized model revealed that steady states lose their stability only in saddle-node bifurcations. By contrast, in the twoand three-layer subsystems Hopf bifurcations as well as saddle-node bifurcations occur. In the present model it was confirmed by numerical simulations that the saddlenode bifurcations mark the onset of bistability, while the Hopf bifurcations mark the onset of oscillations. Our results are therefore in agreement with earlier work [START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF] showing that, without external feedback, at least two layers of the cascade are required to observe oscillatory dynamics, while bistability already occurs in a single layer of the cascade.

Using a combination of generalized and conventional modeling we then identified a parameter region in which complex and potentially chaotic mixed-mode oscillations occur. To determine the exact nature of these oscillations and the parameter ranges in which they can be observed is a promising question for future studies, but will probably require more extensive numerical simulations and the computation of Lyapunov exponents.

Finally, we have investigated the impact of two known feedback loops acting on the MAPK cascade. Our analysis indicates that, under general conditions, the effect of these loops is to reduce bistability and promote oscillatory dynamics in the cascade.

Sustained oscillations in cells can be generated by several different mechanisms [START_REF] Tyson | Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell[END_REF]. The most prominent cause of oscillations is negative feedback, sometimes in connection with a time delay. Explicit negative feedback led to sustained oscillations in various models of the MAPK cascade [START_REF] Kholodenko | Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades[END_REF][START_REF] Shankaran | Rapid and sustained nuclearcytoplasmic ERK oscillations induced by epidermal growth factor[END_REF] and indeed enlarged the oscillatory parameter regime in our model, too. However, negative feedback cannot explain the oscillations without explicit feedback, which must therefore rely on a different mechanism. It is known that in a system with positive feedback that exhibits bistability, a slow process can cause the system to jump between the two stable states. As a result, relaxation oscillations occur that usually have a pulselike shape [START_REF] Kholodenko | Cell signalling dynamics in time and space[END_REF]. Indeed, it was found that in the MAPK cascade, the second layer is in the regime of bistability when the whole system oscillates [START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF]. The bistable double phosphorylation loop is controlled by a slow process in the top level, thus alternating between the two stable states. It is conceivable that a third process in the full cascade could interact with an already oscillating twolayer system and thereby cause chaotic behavior.

An important question is whether oscillatory dynamics in the MAPK cascade play a role in vivo. Based on the early modeling results, it is often assumed that the cascade is utilized as an ultrasensitive switch, and other dynamics appear only as nonfunctional byproducts. Our results illustrate that the oscillatory parameter space is quite large and therefore likely to be encountered in vivo. This view is supported by the recent experimental observation of MAPK oscillations [START_REF] Hilioti | Oscillatory phosphorylation of yeast Fus3 MAP kinase controls periodic gene expression and morphogenesis[END_REF][START_REF] Shankaran | Rapid and sustained nuclearcytoplasmic ERK oscillations induced by epidermal growth factor[END_REF].

Further support for the functional role of oscillations comes from an evolutionary argument. If oscillations were indeed non-functional one could suspect that specific mechanisms suppressing oscillations should have evolved. Yet, the basic topology of the MAPK cascade that can sustain oscillations is conserved through evolution. By contrast, external feedback loops, that could potentially suppress oscillations, differ among organisms and cell types. Moreover, the two examples of feedback loops studied here were found to promote oscillations instead of suppressing them. If the MAPK cascade is indeed used both as a switch and as an oscillator under different physiological conditions, it is intuitive that the system should be very well conserved, as any further mutations are unlikely to maintain both of the cascade's function.

To ascribe biological function also to the irregular os-

A c c e p t e d m a n u s c r i p t

cillations reported here is highly speculative. Nevertheless, it is conceivable that the observed spike-like dynamics could play a useful role. One could alternatively suspect that the complex oscillations are a byproduct of a mechanism that enhances differences between individual cells.

Such a mechanism may be useful because if in a population of cells, a substantial fraction of the populations changed their behavior at a similar level of an external stimulus, then the response would become switch-like also on the population level. This could have detrimental consequences because it could seriously decrease the dynamical stability on the population level and produce unnecessarily strong responses. It may therefore be advantageous to sustain a high sensitivity not only to the stimulus but also to enzyme concentrations that may differ between individual cells. Thereby, signals in different cells could be triggered at different levels of the stimulus, so that the response of the population is smooth, while the response of the individual cell stays sharp. Weakly chaotic dynamics could possibly arise as a byproduct of this internal sensitivity.

The question whether the detected forms of dynamics can appear in vivo is closely connected to the question how realistic the space of parameters assessed in our study is. The ranges that we assigned to the general parameters were motivated by experimental data. While the total concentrations of the involved enzymes have been measured in one cell type [START_REF] Huang | Ultrasensitivity in the mitogenactivated protein kinase cascade[END_REF], less is known about the remaining parameters. Note, however, that large differences in the fluxes f i (see supplementary material) are to be expected since also the total concentration listed in Table 5 vary strongly between different layers. Therefore the parameter ranges in which oscillations and complex dynamics were observed in this study are not unreasonable. Nevertheless, future experimental work covering further variables and different cell types and developmental stages would be highly desirable to clarify the dynamics and function of the MAPK cascade in vivo.

= f 13 (S 1 1 , K 1 R ) -f 14 (S 1 5 ) -f 15 (S 1 5 ) d dt S 1 6 = f 16 (S 1 2 , P 1 R ) -f 17 (S 1 6 ) -f 18 (S 1 6 )
From these equations, the stoichiometric matrix N is constructed according to the rule that N i, j = 1 if f j+6 appears in the equation for S 1 i with a positive sign, N i j = -1 if f j+6 appears in the equation with a negative sign and N i j = 0 else. Note that we use f j+6 for the sake of a consistent notation. Since the first layer of the cascade is not included in this model, the fluxes start with f 7 here.

Next, Λ i j = N i j f j+6 /S 1 i is constructed from the stoichiometric matrix.

Assuming mass action, the matrix of derivatives is given by

θ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ - S 1 1 S 1 0 - S 1 2 S 1 0 - S 1 3 S 1 0 - S 1 3 K 1 R - S 1 4 S 1 0 - S 1 5 S 1 0 - S 1 5 K 1 R - S 1 6 S 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 - S 1 4 P 1 R 0 - S 1 6 P 1 R 0 0 0 1 0 0 0 0 0 1 0 0 1 0 - S 1 3 K 1 R 0 - S 1 5 K 1 R 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 - S 1 4 P 1 R 0 - S 1 6 P 1 R 0 0 0 0 0 1 0 0 0 0 0 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
The entries of θ that are equal to 1 indicate linear dependencies.

The remaining nonzero entries arise due to indirect effects because of mass conservation. Here we show the calculation for θ μ 1

x 1 = ∂μ 1 ∂x 1 S=S * (7)

A c c e p t e d m a n u s c r i p t

Assuming mass action, f 7 = α 7 S 1 0 K 1 R with a rate constant α 7 . It follows that

θ μ 1 x 1 = f 7 f * 7 = S 1 0 K 1 R S 1, * 0 K 1, * R (8)
Because of mass conservation of the substrate and the kinase,

S 1 0 = S 1 T -S 1 1 -S 1 2 -S 1 3 -S 1 4 -S 1 5 -S 1 6 and K 1 R = K 1 T -S 1 3 -S 1 5 .
Therefore, taking the derivative with respect to x 1 = S 1 1 /S 1, * 1 gives

θ μ 1 x 1 = ∂μ 1 ∂x 1 S=S * = S 1, * 1 S 1, * 0 K 1, * R ∂K 1 R S 1 0 ∂S 1 1 S=S * = - S 1, * 1 S 1, * 0 (9)
The Jacobian J = Λθ μ x can then be obtained by matrix multiplication. The resulting matrix is shown in the supplementary material.

Parameter Ranges

In the two-layer and the three-layer model we use parameters inferred from biologically reasonable ranges. We choose the parameters based on data from Xenopus oocyte extracts [START_REF] Huang | Ultrasensitivity in the mitogenactivated protein kinase cascade[END_REF].

Parameter Value S 0 T 2nM S 1 T 1.2μM S 2 T 1.2μM K 0 T 0.3nM P 0 T 0.3nM P 1 T 0.3nM P 2 T 120nM
We vary all total concentrations parameters by factors of 5 in the random sampling. The steady-state concentrations of all dynamic variables are then determined by the algorithm outlined in the main part of the paper.

Since we do not know of experimental data for typical steady-state fluxes f, we tried to infer ranges for these quantities using information on rate constants. For each steady-state flux f i , there is a rate constant α i of the conventional model. From [START_REF] Qiao | Bistability and oscillations in the Huang-Ferrell model of MAPK signaling[END_REF] we adopt ranges of 30 -750min -1 for α 3n+2 and α 3n with n = 1 . . . 10. Combining this with the steady-state concentrations, we can assign ranges to most fluxes. The remaining fluxes are determined by the flux modes. For this reason, and also because continuation of the generalized parameters in the bifurcation analysis alters also the parameters of the conventional models, we cannot guarantee that all α i of conventional models remain restricted to these ranges. The α i of the samples for which we show simulations are given in the supplementary material. Figure 5: Instability in a single layer of the cascade. The components of the eigenvector corresponding to the largest eigenvalue of the Jacobian are analyzed. For each dynamic variable, the color indicates the averaged value of the corresponding entry of the eigenvector. All eigenvectors were normalized so that S 1 5 = 1 before taking the average over 10 6 samples. The diagram shows that the systems leaves the region of a typical unstable state by accumulating mass in either one of the two phosphorylation cycles, while depleting the other.

A c c e p t e d m a n u s c r i p t

Figure 6: Three-parameter bifurcation diagram of a codimension-2 Gavrilov-Guckenheimer bifurcation. The codimension-2 bifurcation is formed at the intersection line of a Hopf bifurcation (light gray) with a saddle-node bifurcation (dark gray). Bifurcations of higher codimension such as Gavrilov-Guckenheimer and double Hopf bifurcations can serve as proxies that indicate parameter regions in which complex dynamics can potentially be observed. See the supplementary material for additional parameters. T and S 2 T . Hopf bifurcations are drawn in red, saddlenode bifurcations in blue. At the green line (AP), a Hopf pair loses its imaginary part, leading to two distinct real eigenvalues. The numbers denote the number of complex and real eigenvalues with a positive real part, in this order. This means that the stable region is marked by 00. The points B and C at which we show numerical integration results are drawn in. B and C: Complex oscillations and chaos close to the double-Hopf bifurcation.

Figure 8: The effect of external feedback. Color-coded is the change in relative abundance of unstable stationary states as feedback is switched on. This change is shown separately for unstable states with leading complex eigenvalues (left) and leading real eigenvalues (right). The top row corresponds to feedback to the first layer (type a) while the bottom row corresponds to feedback to the second layer (type b). In both cases stronger feedback (large values of θ) means that more unstable steady states with a complex and less states with a real leading eigenvalue occur. This suggests that the feedbacks promote rather than inhibit oscillatory dynamics.

  are usually introduced in conventional modeling.

Figure 1 :

 1 Figure 1: Diagrammatric representation of the MAPK signaling cascade. The pathway consists of three layers, a single phosphorylation loop at the top layer and two double-phosphorylation loops at the other layers. Phosphorylated forms are indicated by P's in the figure. The fully phosphorylated product in each level acts as a kinase phosphorylating the proteins on the level below.

Figure 2 :

 2 Figure 2: Schematic overview of the MAPK cascade. The solid arrows correspond to biochemical reactions while the dashed arrows denote the binding of an enzyme in the formation of an substrateenzyme complex. The dotted kinase K 1 R (1-layer model) is replaced by S 0 1 in the larger models. The three different subsystems that we analyzed separately are indicated by boxes. The gray dotted lines represent two types of feedback, which are only included in the model in Sec. 3.7.

Figure 3 :

 3 Figure3: Correlations of parameters with stability. We consider three different subsystems of the cascade, a single layer (A), two connected layers (B) and three layers (C). The bars show the correlation of a given parameter with the stability of the steady states in a large ensemble of randomly generated parameter sets. The statistical errors, which are determined by bootstrap resampling, are on the order of the line width. The plots show a conserved pattern of correlations and indicate the parameters that are important for stability.

Figure 4 :

 4 Figure 4: Connection between steady-state concentrations and stability. Contour lines computed from 10 7 random parameter sets show the fraction of stable systems depending on the two concentrations. The top left panel shows that instability is promoted by low concentrations of free kinase and phosphatase, highlighting the importance of sequestration effects. The model parameters with the strongest impact on stability are the concentrations of substrates S 1 3 and S 1 6 (top right).

Figure 7 :

 7 Figure 7: Complex dynamics close to a double-Hopf bifurcation. A: Two-dimensional bifurcation diagram with a double Hopf bifurcation. Parameters are S 0T and S 2 T . Hopf bifurcations are drawn in red, saddlenode bifurcations in blue. At the green line (AP), a Hopf pair loses its imaginary part, leading to two distinct real eigenvalues. The numbers denote the number of complex and real eigenvalues with a positive real part, in this order. This means that the stable region is marked by 00. The points B and C at which we show numerical integration results are drawn in. B and C: Complex oscillations and chaos close to the double-Hopf bifurcation.

  Figure 1

  Figure 2

Acknowledgments

The authors would like to thank Mark Fricker for insightful discussion of the results, especially regarding the potentially beneficial effects of weak chaos.

Appendix

Construction of the Jacobian

We demonstrate the construction of the Jacobian for the one-layer system, consisting of two phosphorylation steps. In a general form, the associated system of ODEs is


Giulio Caravagna
email: caravagn@di.unipi.it

Alberto D'onofrio
email: alberto.donofrio@ifom-ieo-campus.it

Paolo Milazzo
email: milazzo@di.unipi.it

Roberto Barbuti
email: barbuti@di.unipi.it

Tumour suppression by immune system through stochastic oscillations

Keywords:

The well-known Kirschner-Panetta model for Tumour-Immune System interplay [Kirschner and Panetta, J. Math. Biol 37 (3), 1998] reproduces a number of features of this essential interaction, but it excludes the possibility of tumour suppression by the immune system in the absence of therapy. Here we present a hybrid-stochastic version of that model. In this new framework, we show that in reality the model is also able to reproduce the suppression, through stochastic extinction after the first spike of an oscillation.

Introduction

The tumour-immune system interaction takes place because tumour cells are characterised by a vast number of genetic and epigenetic events leading to the appearance of specific antigens, called neoantigens, that trigger antitumoural actions by the immune system [START_REF] Pardoll | Does the Immune System See Tumours as Foreign or Self?[END_REF]. These observations provided a theoretical basis [START_REF] Burnet | Cancer -a biological approach[END_REF] to the empirical hypothesis of immune surveillance, i.e. that the immune system may act to control and also, in some case, to eliminate tumours [START_REF] Ehrlich | Ueber den jetzigen Stand der Karzinomforschung[END_REF]. Only in recent years has the study of cancer immunobiology accumulated a sufficient amount of evidence to show that tumours may be suppressed by immune system effectors. This has been made possible by the application of new molecular techniques and by performing a large number of epidemiological studies [START_REF] Dunn | The three ES of Cancer Immunoediting[END_REF].

The competitive interaction between tumour cells and the immune system involves a considerable number of events and molecules in a extremely complex way. As a consequence, the immune system is not able to eliminate a neoplasm in all cases, since it may escape from its control. Of course, a dynamic equilibrium may also be established, such that the tumour may survive in a microscopic steady state, which is undetectable by diagnostic equipment (d'Onofrio 2007). However, let us consider a tumour which is constrained by the immune system in a microscopic steady state. On the one hand, the presence of a small colony of tumour cells may result in metastases. On the other hand, over a long period of time, a significant fraction of the mean human life span, according to [START_REF] Dunn | The three ES of Cancer Immunoediting[END_REF]), the neoplasm may develop multiple strategies to circumvent the action of the immune system [START_REF] Whiteside | Tumour-induced death of immune cells: its mechanisms and consequences[END_REF][START_REF] Pardoll | Does the Immune System See Tumours as Foreign or Self?[END_REF][START_REF] Dunn | The three ES of Cancer Immunoediting[END_REF][START_REF] Vicari | Tumour escape from immune surveillance through dendritic cell inactivation[END_REF]. This, in the long term,

A c c e p t e d m a n u s c r i p t

may allow it to evade immune surveillance and to re-commence growing to its carrying capacity.

The tumour has itself adapted to survive in a hostile environment, in which the antitumour immune response is activated (d'Onofrio 2007[START_REF] Cattani | Metamodeling the learning-hiding competition between tumours and the immune system: A kinematic approach[END_REF]). However, the tumour-immune system interaction may also show attractors that are different from a constant equilibrium, namely oscillations. In fact, one may observe both 'short term-small amplitude' oscillations [START_REF] Kennedy | Cyclic leukocyte oscillations in chronic myelogenous leukemia during hydroxyurea therapy[END_REF][START_REF] Vodopick | Spontaneous cyclic leukocytosis and thrombocytosis in chronic granulocytic leukemia[END_REF][START_REF] Gatti | Cyclic Leukocytosis in Chronic Myelogenous Leukemia: New Perspectives on Pathogenesis and Therapy[END_REF][START_REF] Mehta | Cyclic oscillations in leukocyte count in chronic myeloid leukemia[END_REF] and patterns of remission-recurrence such as alternation of long phases of dormancy (where a fallacious apparent immune surveillance seems reached) followed by phases where the tumour seems to escape [START_REF] Tsao | Ultra-late recurrence (15 years or longer) of cutaneous melanoma[END_REF][START_REF] Blumberg | The relationship of blook transfustion, tumor staging, and cancer recurrence[END_REF][START_REF] Sohrabi | Recurrence of breast cancer: Obesity, tumor size, and axillary lymph node metastases[END_REF]. Finally, the study of the tumour-immune system interaction led to the proposal and implementation of an interesting therapeutic approach: immunotherapy [START_REF] Vito | Cancer: principles and practice of Oncology[END_REF], consisting in stimulating the immune system in order to better fight, and hopefully eradicate, a cancer. The basic idea of immunotherapy is simple and promising, but the results obtained in medical investigations are globally controversial [START_REF] Agarwala | New Applications of Cancer Immunotherapy[END_REF], Bleumer et al 1003[START_REF] Kaminski | Immunotherapy and prostate cancer[END_REF].

Regarding the mathematical modelling of the above interactions, many papers have appeared using either a finite dimensional approach [START_REF] De Lisi | Immune surveillance and neoplasia: a minimal mathematical model[END_REF][START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF][START_REF] Stepanova | Course of the immune reaction during the development of a malignant tumour[END_REF], Kuznetzov et al. 1994, Kuznetzov and Knott 2001[START_REF] Galach | Dynamics of the tumour-immune system competition: the effect of time delay[END_REF][START_REF] Szymanska | Analysis of the immunotherapy models in the context of cancer dynamics[END_REF][START_REF] De | Dynamic response of cancer under the influence of immunological activity and therapy[END_REF][START_REF] De Pillis | A Validated Mathematical Model of Cell-Mediated Immune Response to Tumour Growth[END_REF], d'Onofrio 2005, d'Onofrio 2007[START_REF] Cappuccio | Cancer immunotherapy by interleukin-21 treatment strategies evaluated in a mathematical model[END_REF][START_REF] Kronik | Improving alloreactive CTL immunotherapy for malignant gliomas using a simulation model of their interactive dynamics[END_REF][START_REF] Kirschner | On the global dynamics of a model for tumor immunotherapy[END_REF] or the theory of kinetic active particles (Bellomo and[START_REF] Bellomo | From the mathematical kinetic, and stochastic game theory to modelling mutations, onset, progression and immune competition of cancer cells[END_REF]. Some of the above works focus on sustained oscillations that are related to delays in the process of detecting the presence of the tumour by the immune system [START_REF] Galach | Dynamics of the tumour-immune system competition: the effect of time delay[END_REF][START_REF] Szymanska | Analysis of the immunotherapy models in the context of cancer dynamics[END_REF][START_REF] Cattani | Metamodeling the learning-hiding competition between tumours and the immune system: A kinematic approach[END_REF], whereas others show that oscillations are also possible in the case of absence of such delays (d'Onofrio 2006, Kirschner and[START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF]. For example, in (d'Onofrio 2006) an analytical study is conducted and conditions are stated for the existence and uniqueness of a limit cycle for a large class of models.

In [START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF], Kirschner and Panetta proposed a largely influential tridimensional deterministic model, whose variables are tumour cells, effector cells and the concentration of interleukins. This model is able to explain both the above type of tumour size oscillations, as well as, of course, constant equilibria, both macroscopic and microscopic. The global behaviour of the model was recently investigated in [START_REF] Kirschner | On the global dynamics of a model for tumor immunotherapy[END_REF]. Although a vast array of behaviours is mimicked by the solutions of the Kirschner-Panetta model, the tumour-free equilibrium is unstable for all biologically meaningful values of the parameters of the model.

However, we noticed that in the transitory oscillations before reaching a stable equilibrium, and also in sustained stable large oscillations, very small values of the tumour size are predicted by this deterministic model. These values could be considered well below the level of detection but they cannot be considered as eradication. Thus, we propose here a stochastic version of the Kirschner-Panetta model, which may better describe the dynamics of the tumour-immune system interplay in the case of low levels of both tumour and immune cells. In the bio-mathematical literature it has been stressed how deterministic models may exhibit oscillating solutions having no counterpart in their corresponding stochastic version, where population extinction is, instead, predicted [START_REF] Nasell | Measles outbreaks are not chaotic, in Mathematical approaches for Emerging and Reemerging Infectious Diseases: Models, Methods, and Theory[END_REF].

In order to assess whether this also happens in the Kirschner-Panetta model, we propose here an extended model where the dynamics of tumour cells and immune effectors are described by a suitable stochastic process, whereas the dynamics of the concentration of interleukins is modelled by an ordinary differential equation. Since the stochastic processes are interlinked with the dynamics of the interleukins (and vice versa), the main features of this model are that it contains time-dependent propensity functions and combines both the features of deterministic and stochastic models.

This provides a useful tool through which we aim to offer a theoretical contribution towards a better understanding of the dynamical aspects of a key phenomenon: the suppression of a tumour by the immune system. In particular, we focused on: i) assessing the role of antigenicity of tumour cells [START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF] in shaping those dynamics; ii) investigating the interplay between the so-called intrinsic noise and the macroscopic tendency of tumour-immune system interplay to oscillate (d'Onofrio 2006, Kirschner and[START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF].

A c c e p t e d m a n u s c r i p t

The present work is organised as follows: in Section 2 we present the Kirschner-Panetta model and we summarise its main features, with particular regard to those linked to the onset of sustained oscillations. In Section 3 we propose a hybrid stochastic version of the model. In Section 4.1 we use the stochastic Kirschner-Panetta model to investigate the role of stochasticity in tumour suppression. The work ends with some concluding remarks.

Overview of the Kirschner-Panetta model

In [START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF] the following model of the dynamics of tumour-immune system interaction was proposed:

T * = rT * (1 -bT * ) - aT * g 2 + T * E * E * = p 1 I g 1 + I E * -μ 2 E * + cT * + s 1 (1)
I = p 2 T * E * g 3 + T * -μ 3 I + s 2
where T * (t), E * (t) and I(t) denote, respectively, the densities of tumour cells, effectors of the immune system and interleukins. The tumour induces the recruitment of the effectors at a linear rate cT , thus, the parameter c may be seen as a measure of the immunogenicity of the tumour, i.e. 'a measure of how different the tumour is from self' [START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF]. The proliferation of effectors is stimulated by the interleukins. A continuous infusion immunotherapy is delivered; in the general case, both effectors and interleukins are injected at rates s 1 ≥ 0 and s 2 ≥ 0. The source of interleukin is modelled as linearly depending on effectors, and it also depends on the tumour burden. Finally, the average lifespan of effectors is μ -1 2 and the average degradation time for interleukin is μ -1 3 . In the absence of therapy (s 1 = s 2 = 0), the main results obtained in [START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF] are the following: the unique tumour suppression equilibrium point (0, 0, 0) is always unstable, there is a value 0 < c 0 << 1 such that for 0 < c < c 0 there are three equilibria of which two are unstable and one is locally stable. As it is intuitive, the antigenicity of the tumour being very small, this locally stable equilibrium is very close to the carrying capacity in the absence of immune response. Furthermore, there is a value c max > 0 such that if c ∈ (c 0 , c max ) there is a unique periodic solution; the period and amplitude of the limit cycle decreases as far as c approaches c max . Finally, when c > c max there is a unique globally stable equilibrium, which is reached through damped oscillations. Again, this equilibrium is a decreasing function of c, and its size is small due to the high immunogenicity of the tumour. The presence of a very small equilibrium state that is not detectable might seem to be a sub-optimal outcome, so that one may consider it as practically equivalent to tumour suppression. However, it must be remembered that, as we mentioned in the introduction, the Volterrian interaction between tumour cells and the immune system is also evolutionary: tumour cells evolve their ability to evade the control [START_REF] Dunn | The three ES of Cancer Immunoediting[END_REF], d'Onofrio 2007) and, as a consequence, no microscopic steady state may be considered permanent and safe for the patient (d'Onofrio 2007).

Thus, we note that this model explicitly precludes the possibility of tumour suppression, which follows from the absence of a baseline influx of effectors. We remark that this holds in the absence of therapy, and that is the case which we will consider in the rest of the paper. In the presence of therapies, the behaviour of the system becomes more complex, but in all the possible 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 the injection rate of the interleukins greater than such a value, there is an unbounded growth of effectors, leading to capillary leak syndrome.

A c c e p t e d m a n u s c r i p t 3 A Hybrid Stochastic version of the Kirschner-Panetta model

In this section we define a hybrid stochastic version of the model (1). In principle, indeed, before defining a stochastic equivalent of the model (1), we should switch from variables representing densities to variables representing the total number of cells. Thus, we would obtain an ODE model from which it would be straightforward to obtain a three dimensional stochastic process. However, since the molecular weight of interleukins is 15000 Dalton, the average number of proteins is huge. As a consequence, we may assume that the dynamics of I(t) is well approximated by a linear differential equation with randomly varying coefficients, which, however, in the intervals between two consecutive stochastic events are evidently constant.

As a consequence, we introduce as a new parameter the volume V of interest (e.g. the blood and bone marrow volumes for leukaemia) and we define:

T * = T V , E * = E V (2)
thus obtaining:

T = rT 1 - b V T - aT g 2 V + T E E = p 1 I g 1 + I E -μ 2 E + cT + V s 1 (3) I = p 2 V T E g 3 V + T -μ 3 I + s 2
From the above model, we define a hybrid stochastic process, where a bi-dimensional stochastic process is linked to a scalar differential equation ruling the dynamics of the concentration of interleukins.

Namely, in the interval between the n-th event (at time denoted as t n ) and the (n + 1)-th event, the above mentioned ODE is:

I = p 2 V T (t n )E(t n ) g 3 V + T (t n ) + s 2 -μ 3 I (4)
which, since in (t n , t n+1 ) the other two state variables are constant, is a linear ordinary differential equation with constant input and constant coefficients. Thus, the solution of (4) in (t n , t n+1 ) has the following analytical form:

I(t) = B n + (I n -B n ) Exp (-μ 3 (t -t n )) (5)
where B n , for the sake of notation simplicity, is the following function of (T n , E n ):

B n := 1 μ 3 p 2 V T (t n )E(t n ) g 3 V + T (t n ) + s 2 . ( 6
)
The stochastic model of the evolution of the tumour cells and the effector cells is defined by the following events and the corresponding propensity functions:

(R 1 ) T → T + 1 a 1 = r 2 T (R 2 ) T → T -1 a 2 = r 2 bV -1 T 2 (R 3 ) T → T -1 a 3 = (aT E)/(g 2 V + T ) (R 4 ) E → E + 1 a 4 (t) = p 1 E(I(t)/(g1 + I(t))) (R 5 ) E → E -1 a 5 = μ 2 E (R 6 ) E → E + 1 a 6 = cT (R 7 ) E → E + 1 a 7 = V s 1
Notice that all the propensity functions except for a 4 (t) are time-constant. In order to analyse this model, we need to define an algorithm for computing its time evolution. Although in the literature

A c c e p t e d m a n u s c r i p t

algorithms for simulating hybrid systems have been presented [START_REF] Salis | Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions[END_REF]Kaznessis 2005, Alfonsi et al. 2004), in this model we combine the hybrid approach with time-dependent propensity functions [START_REF] Lecca | A time-dependent extension of Gillespie algorithm for biochemical stochastic πcalculus[END_REF][START_REF] Anderson | A modified next reaction method for simulating chemical systems with time dependent propensities and delays[END_REF]. Consequently, we define, based on the main ideas of Gillespie's algorithm [START_REF] Gillespie | A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions[END_REF][START_REF] Gillespie | Exact Stochastic Simulation of Coupled Chemical Reactions[END_REF], a stochastic evolution algorithm for this hybrid model. The event-induced changes in the values of T and E can be represented by the following sevendimensional vectors ν T = (1, -1, -1, 0, 0, 0, 0) T r and ν E = (0, 0, 0, 1, -1, 1, 1) T r .

Namely, ν T [i] and ν E [i] denote how the event R i affects the population T and E, respectively. Furthermore, we denote with S the set of time-constant propensity functions, namely the set {a 1 , a 2 , a 3 , a 5 , a 6 , a 7 }.

The unique time varying propensity is given by the formula:

a 4 (t) = E n p 1 I(t) g 1 + I(t)
.

We define the function:

A(τ ) = tn+τ tn a 4 (t)dt.
Given the analytical form of I(t) in formula (5), A(τ ) reads as follows:

A(τ ) = p 2 E n (B n + g 1 ) μ 3 τμ 3 (B n + g 1 ) + g 1 Ln g 1 + I n -B n + e τ μ3 (B n + g 1 ) + I n . ( 7
)
The exact simulation algorithm for the hybrid system under study can consequently be defined as follows:

1. Input the initial state: (T 0 , E 0 , I 0 ) and the maximum simulation time Θ M 2. Initialise the time t = t 0 and the system state (T, E, I) = (T 0 , E 0 , I 0 ).

3. If t ≥ Θ M end the simulation, else 4. Evaluate all the a j with j ∈ S and their sum a 0 = j∈S a j .

5. The putative time for the next reaction τ is determined by solving the following transcendental equation:

A(τ ) + a 0 τ = χ,
where χ is a random number with distribution Exp(1). For the numerical solution of this equation we used Newton-Raphson method.

6. Given a random number r uniformly distributed in the interval [0, 1], the next event is determined in accordance with

j-1 i=1 a i < r • a 0 ≤ j i=1 a i 7. Update the system state to (T + ν T [j], E + ν E [j], I(t + τ ))
and the clock to t + τ .

8. Go to step 3.

Finally, in order to perform the simulations of both the deterministic and stochastic model, we have to define numerical values for the parameters.

Concerning the parameters, we used the values and ranges given in the references (Kirschner andPanetta 1998, Kirschner et al 2004), as listed in Table 1. Note that those values pertain to mice, since they were taken from the influential papers [START_REF] Deboer | Macrophage T Lymphocyte interactions in the anti-tumor immune response: A mathematical model[END_REF], Kuznetzov et al. 1994), where accurate fitting of real data concerning laboratory animals were performed. As far as parameter V is concerned, by taking into account that in a chimeric mouse the body weight ranges from 20 grams for female mouse up to 40 grams, and that their blood volume ranges from 5.8 ml per 100 grams up to 8 ml per 100 grams, it follows that a reasonable range for V is from 1.16 ml up to 3.2 ml, in our simulations we assumed: V = 3.2ml.

In our simulations we focused on two key parameters: c, which reflects the antigenic potential of the tumour, and s 1 which represents the injection rate of effectors for adoptive cellular immunotherapy.

Stochastic behaviour

Stochastic suppression of tumour in absence of immunotherapy

In [START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF]) it has been shown that the antigenicity parameter c is a key parameter for model (1) since, for very low values of c, the tumour burden reaches values near its carrying capacity, for large values of c the tumour reaches a small steady state, whereas, for intermediate values, lying in a quite narrow range R = (c min , c max ), there is the onset of oscillations through Hopf bifurcations at c = c min and at c = c max . This is biologically sound since, of course, the immunogenicity of a tumour is an essential factor so that the immune system may attack the tumour with some effectiveness.

However, simulating (3) we noticed that, although for c slightly higher than c min the maximum value attained during the oscillations may be of the order of the carrying capacity, the minima of T (t) seem very small. Indeed, if we plot time courses of tumour size in a Log10 scale, as in Figure 1, we may easily see that there are temporal ranges where the deterministic setting predicts that T (t) < 1 cell.

This suggested to us that a range Ω ⊂ R might exist where, by using the stochastic version of the Kirschner-Panetta model, the oscillating interplay tumour-effectors-interleukins might lead to the eradication of the disease.

This hypothesis was confirmed by the stochastic simulations. All the simulations of the hybrid model weredone by using a Java implementation of this model and of the algorithm presented in Section 3. The Java source code of the implementation is freely available at the URL: http://www.di.unipi.it/msvbio.

We present, in Figure 2, two realisations (for c = 0.01 and for c = 0.02) of the stochastic process corresponding to the deterministic time-courses shown in Figure 1. We set the worst possible initial conditions for effectors and interleukins, i.e. E(0) = 0, I(0) = 0.

It is immediately noticeable that the stochastic behaviour is deeply different from that predicted by the deterministic model since, in the stochastic simulation, only the first spike of the oscillation is achieved and then the tumour eradication is reached. In particular, with c = 0.01 eradication is observed around day 160 and with c = 0.02 around day 120. The other parameter are as in Table 1.

Thus, denoting as t er the random variable time of occurrence of the eradication we calculated the associated empirical probability density (t er ; c) with c = 0.02. The result is shown in Figure 3. Such a density, computed over 10 3 stochastic runs of the model, has a peak at around day 130; all the simulations have shown eradication with a simulation time bounded to 200 days.

Note that the stochastic simulations for some values of c become computationally intractable for sample hundreds of runs. For example, at c = 0.01 a single stochastic run took almost 70 hours on a Apple MAC1 .

Consequently, we decided to perform a deeper analysis on the deterministic model, which is much less computationally expensive, considering as eradication the first time instant in which T (t) < 1. We calculated, for the values of Log 10 (c) in [-4; -1] the average eradication time t er (c) in function of c, as reported in Figure 4. We note that for c approaching c min the putative eradication time becomes very large. However, notice that in reality the patient dies before that the carrying capacity is reached. Moreover, after decreasing for nearly all the range Ω, the curve has a minimum and the eradication time restarts growing. Notice that for c = 0.02 (∼ Log 10 (c) = -2.6) the analysis of the deterministic model gives roughly similar predictions of the stochastic behaviour observed in Figure 3.

In order to assess the influence of non-null initial value for effectors, in Figure 5, we show the effect of increasing E 0 on the empirical probability density of t er with c = 0.02 in the stochastic model. These densities are computed over 10 2 stochastic runs. All the initial configurations we considered are of the form (T 0 , E 0 , I 0 ) = (1, 10 i , 0) with i = 1, . . . , 5. In particular, for i = 1, 2, 3 the densities have the peaks in the interval [120; 140] while, for i = 4, 5 the densities become bimodal with peaks in the intervals [120; 140] and [0; 10]. For the particular initial configurations, which are not shown in the figure, having 10 6 or more initial effectors, all the eradications are observed almost instantaneously.

The same kind of analysis was performed by varying the initial numbers of interleukins. In particular, in Figure 6, we show the effect of increasing I 0 on the empirical probability density of t er with 1. In both cases there is tumour suppression, and further spikes of the original oscillations disappear. 1. c = 0.02 in the stochastic model. Also these densities are computed over 10 2 stochastic runs. All the initial configurations we considered are of the form (T 0 , E 0 , I 0 ) = (1, 0, 10 i ) with i = 1, 2, 4, 6. In this case, unlike the previous ones, we observed single modal probability densities with peaks in the interval [110; 135].

Stochastic oscillations in absence of immunotherapy

The stochastic version of the Kirschner-Panetta model exhibits stochastic oscillations for a small range of values of the antigenicity 0.03 ≈< c <≈ 0.032, whereas for higher values damped oscillations appears (i.e. with c = 0.033 the system reaches an equilibrium with around 60000 tumour cells). In Figure 7 we show two examples of oscillations; for c = 0.03 and for c = 0.032. Note that in Figure 7 we also plotted the minimum after the first spike of the oscillation. Figure 5: Empirical probability density of t er with c = 0.02 and variable initial size of effectors: E(0) = 10 i with i = 1 (solid line), i = 2 (dashed line), i = 3 (semi-dashed line) on the left and i = 4 (solid line) and i = 5 (dashed line) on the right. For each configuration the density has been computed over 10 2 stochastic runs. The other parameters are as in Table 1.

Tumour suppression in the case of large c

If the value of c is further increased the average value of the stochastic oscillations of tumour cells progressively reduces. The reduction is such that the suppression of the tumour gradually reappears.

For example, we performed a series of simulations by varying c (2000 simulations 2000 days long for each value of c), and we observed that at least one observed eradication occurred for c greater than a threshold value 13, the value for which two eradications over 2000 simulations of the stochastic process were observed. At c = 30 tumour eradication was observed in 20% of simulations, and for c = 50 almost half of the simulations had the tumour eradication as outcome (956 out of 2000).

In Figure 8 one simulation of the stochastic process for c = 71.1 is shown, whereas in Figure 9 (left) we plotted the percentage of observed eradications in function of c. As far as the times when the first eradication is observed, we note that in many cases, the minimal eradication time is very small. Thus, we plotted in Figure 9 (right) the modal value of t erad for c ∈ (13, 71.5), i.e. the value of t er for which the maximum number of eradications was observed.

Immunotherapy

Although, as we stressed in the introduction, here we are essentially interested in investigating the effect of stochasticity in the unperturbed time course of a tumour, for the sake of completeness, we briefly report some simulations of an adoptive cellular immunotherapy (s 1 > 0). IL(0) = 10 i with i = 1 (solid line), i = 2 (dashed line), i = 4 (semi-dashed line) and i = 6 (dotted line). For each configuration the density has been computed over 10 2 stochastic runs. The other parameter are as in Table 1. The other parameters are as in Table 1. In both subfigures an insert is plotted zooming in on the minimal values reached at the end of the first spike of the oscillation.

The deterministic setting is quite pessimistic since a therapy only eradicates for t → ∞, so that after a therapy having a realistic finite time span the tumour restarts growing. Thus the maximum beneficial effect of an immunotherapy would be moving the state of tumour from the region of convergence of a macroscopic (not compatible with the life of the host organism) to the basin of attraction of a small steady state.

The stochastic simulations, instead, suggests that an immunotherapy may eradicate the disease in a finite time, as shown, for the case of an adoptive cellular immunotherapy (Figure 10).

We focused on the antigenicity value c = 10 -4 which is, for both the deterministic and stochastic models, a value such that no eradication is observed in the absence of therapy. Regarding the therapy, we modelled an adoptive cellular immunotherapy from day 70 to day 130 of simulation. We used s 1 = 0.0 for the first 70 days of simulation and, for the following days, we used s 1 = 10s1 cr where s1 cr is the threshold value for the local stability of the suppression in the deterministic setting [START_REF] Kirschner | Modeling immunotherapy of the tumour -immune interaction[END_REF].

Concluding Remarks

In order to investigate the phenomenon of spontaneous onco-suppression by the immune system and at the same time to assess whether the Kirschener-Panetta approach may reproduce this important The other parameters are as in Table 1.

phenomenon, we proposed here a hybrid model that extends the Kirschner-Panetta model. In our extended model the dynamics of tumour cells and immune effectors are described by a suitable stochastic process, whereas the dynamics of the concentration of interleukins are modelled by an ordinary differential equation. Since the stochastic processes are interlinked with the dynamics of the interleukins (and vice versa), the main features of this model are that it contains time-dependent propensity functions and combines both the features of deterministic and stochastic models.

In order to simulate the model, it was necessary to define a proper stochastic simulation algorithm. Indeed, although in the literature we find algorithms for simulating time-dependent propensity functions [START_REF] Lecca | A time-dependent extension of Gillespie algorithm for biochemical stochastic πcalculus[END_REF][START_REF] Anderson | A modified next reaction method for simulating chemical systems with time dependent propensities and delays[END_REF] as well as algorithms for simulating hybrid models [START_REF] Salis | Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions[END_REF]Kaznessis 2005, Alfonsi et al. 2004), the combination of both is missing. Consequently, for the sake of analysing this model we defined an ad-hoc algorithm based on the ideas of Gillespie's Stochastic Simulation Algorithm [START_REF] Gillespie | A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions[END_REF][START_REF] Gillespie | Exact Stochastic Simulation of Coupled Chemical Reactions[END_REF] which permits an efficient implementation.

We would point out that the model we present is not a completely new one, but is a stochastic model built on top of a well-known and biologically backgrounded model. The simulations and their biological interpretation, instead, are new.

Our hybrid stochastic version of the Kirschner-Panetta model suggests that, at least in some cases, complete immune surveillance with suppression of the neoplasm might be reached thanks to the conjunction of the intrinsic tendency of the tumour-immune system to oscillate, which is significantly evidenced by the deterministic model, with the stochastic dynamics. We note that this may empirically be inferred by the deterministic Kirschner-Panetta model by considering the tumour eradicated whenever T (t) < 1. In fact, from a quantitative analysis of the models, we observe significant differences in the tumour eradication times. Moreover, by following the stochastic approach we have 1.

been able to give a probability distribution of the eradication time, and the discussion of Section 4.1 regarding the different behaviours of the system by varying its initial conditions would not have been possible by considering only the average eradication time provided by the deterministic model. The relevant role of oscillations showed an intriguing role for the immunogenicity of the disease, which is represented by the parameter c in the model. Indeed, the deterministic immunogenicity path can be summarised as follows: for low values of c the tumour size is large and the oscillations have very large periods, for intermediate values of c we observe intermediate oscillating tumour size. Finally, for large values of c we observe small tumour size reached after damped oscillations. Conversely,, the stochastic path forecasts immune suppression for intermediate c, as well as, of course, for quite large c. This is quite a surprising result which deserves further biological investigations and which may have interesting implications for some immunotherapies, since for intermediate c the suppression range is small. Moreover, since for many cases of tumour suppression in the case of low c tumour sizes do reach very large values, we must stress that in these cases the eradication is only theoretical since the host unfortunately dies. However, for larger values of c the tumour size is large but compatible with life, so that the model may contribute to a better understanding some sporadic cases of spontaneous permanent regression of macroscopic tumours [START_REF] Krikorian | Spontaneous regression of non/hodgkin lymphoma: a report of nine cases[END_REF]. Of course for large values of c the probability of having a remission is large (e.g. for c = 50 we observed that the suppression is reached in 50% of the simulated paths). Finally, we observe that also stochastic sustained oscillations are possible.

When considering adoptive cellular immunotherapies, our simulations suggest that a stochastic model may qualitatively be used in tuning the average length of the therapies, as briefly illustrated in Figure 10.

In summary, our work shows that the Kirschner and Panetta model is even more powerful than it seems, since our stochastic version of it may also describe the most important phenomenon of the tumour-immune system interplay, namely the total suppression of the tumour by the immune system in the absence of human intervention.

Finally, as far as the phenomena of late recurrence of tumours are concerned, since they did not appear in our stochastic simulations, we think that they may have their main root in the phenomenon of immunoediting [START_REF] Dunn | The three ES of Cancer Immunoediting[END_REF], d'Onofrio 2007) and we intend to study them, following the approach used in (d'Onofrio 2007), by allowing, in the stochastic model, a slow decrease of some key parameters.

A c c e p t e d m a n u s c r i p t

and the authors are grateful for this support. 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. We wish to thank the anonymous referees for their suggestions that contributed to significantly improving the paper.

  -1 baseline strength of the killing rate by immune effectors g2 10 5 ml -1 50% reduction factor of the killing rate by immune effectors p1 0.1245 ml × days -1 baseline strength of the interleukin-stimulated growth rate of effectors g1 2 × 10 7 ml -1 50% reduction factor of interleukin-stimulated growth rate of effectors μ2 0ml -1 × days -1 injection rate of effectors for adoptive cellular immunotherapy s2 0.0 ml -1 × days -1 injection rate of interleukins for immunotherapy

Figure 1 :

 1 Figure1: Deterministic dynamics of T (t) for c = 0.01 (upper panels), and c = 0.02 (lower panels). In left panels the curves are plotted in natural scale, whereas in right panels the plot is in Log10 scale. The other parameter are as in Table1.

Figure 2 :

 2 Figure 2: Stochastic dynamics of T (t) for c = 0.01 (left panel), and c = 0.02 (right panel). The other parameter are as in Table1. In both cases there is tumour suppression, and further spikes of the original oscillations disappear.

Figure 3 :

 3 Figure3: Empirical probability density of t er with c = 0.02 computed over 10 3 stochastic runs. The other parameters are as in Table1.

Figure 4 :

 4 Figure 4: Eradication time vs. Log 10 (c) for the deterministic model. The other parameters are as in Table1. Left panel: complete plot, right panel: zoom to final part.

Figure 6 :

 6 Figure6: Empirical probability density of t er with c = 0.02 and variable initial size of interleukins: IL(0) = 10 i with i = 1 (solid line), i = 2 (dashed line), i = 4 (semi-dashed line) and i = 6 (dotted line). For each configuration the density has been computed over 10 2 stochastic runs. The other parameter are as in Table1.

Figure 7 :

 7 Figure 7: Stochastic oscillations for c = 0.03 (left) and c = 0.032 (right). The other parameters are as in Table1. In both subfigures an insert is plotted zooming in on the minimal values reached at the end of the first spike of the oscillation.

Figure 8 :Figure 9 :

 89 Figure 8: Stochastic suppression for c = 1.1c * . The other parameters are as in Table1.

Figure 10 :

 10 Figure10: Deterministic vs stochastic models of adoptive cellular immunotherapy with c = 10 -4 and s 1 = 5400 delivered for 60 days after day 70. The other parameter are as in Table1.

Table 1 :

 1 Parameters used in most of the simulations of the model. Tumour antigenicity c is the only variable parameter that is reported in the description of each simulation.

Apple Macbook

2.4 Ghz Intel Core 2 Duo, 4GB 1067 Mhz running Mac OS X 10.5.8.

Acknowledgements

The work of A. d'Onofrio was conducted within the framework of the EU Integrated Project "Advancing Clinico-Genomic Trials on Cancer -ACGT". The ACGT project is partly funded by the EU