Deterministic and Stochastic Study for a Microscopic Angiogenesis Model: Applications to the Lewis Lung Carcinoma

Angiogenesis modelling is an important tool to understand the underlying mechanisms yielding tumour growth. Nevertheless, there is usually a gap between models and experimental data. We propose a model based on the intrinsic microscopic reactions defining the angiogenesis process to link experimental data with previous macroscopic models. The microscopic characterisation can describe the macroscopic behaviour of the tumour, which stability analysis reveals a set of predicted tumour states involving different morphologies. Additionally, the microscopic description also gives a framework to study the intrinsic stochasticity of the reactive system through the resulting Langevin equation. To follow the goal of the paper, we use available experimental information on the Lewis lung carcinoma to infer meaningful parameters for the model that are able to describe the different stages of the tumour growth. Finally we explore the predictive capabilities of the fitted model by showing that fluctuations are determinant for the survival of the tumour during the first week and that available treatments can give raise to new stable tumour dormant states with a reduced vascular network.


Introduction
Angiogenesis, the formation of new vessels from pre-existing ones, is a normal and vital process in growth and development of animal organisms. It is required during the repair mechanism of damaged tissues like wound healing processes. However, angiogenesis is also essential in the transition of avascular forms of solid tumours into cancers that are able to metastasise and cause the lethal outcome of the disease. Hence, angiogenesis promotes cancer growth.
In general, for tumours reaching size of 1-2 mm 3 the necrotic core inside the tumour is formed [1,2]. At this stage, cancer cells start to secrete angiogenic factors e.g. FGF, VEGF, VEGFR, Ang1 and Ang2, which promote proliferation and differentiation of endothelial cells, smooth muscle cells, and fibroblasts, initiating the formation of new blood vessels. Those vessels provide the nutrients for growing cancer mass and help to remove the metabolism waste products. On the other hand, achieving a better control on angiogenesis would also improve the cancer treatment allowing anti-cancer drugs to penetrate efficiently the tumour structure through a good functioning vessel network, and hence reduce the tumour mass.
One of the most well known mathematical models describing the influence of angiogenesis processes on the tumour dynamics was proposed by Hahnfeldt et al. [3]. This model and its extensions were later studied, among others, by d'Onofrio & Gandolfi [4][5][6] and by Piotrowska & Foryś [7,8]. Additionally, within last years, different treatment protocols were introduced to the Hahnfeldt et al. type models, and optimal control theory was applied to these modified models to find a proper timing and dosing of drug administration (see e.g. [9][10][11][12][13]), settling those models as the paradigm of angiogenesis modelling. However, global stability of a positive steady state is a weak point of that type of models without delays, since newly formed vessels usually have highly unstable structure. Moreover, feedback loops observed in the biological systems can lead to oscillatory dynamics, [14]. To reflect the complex nature of the vessels formation process, Arakelyan et al. [15] proposed an intricate computational model, which was compared with implanted human ovarian carcinoma reported in [16]. Next, that complex model was simplified to the system of three equations with two time delays. The delays presented in this latter model reflect the length of feedback loops considered as detailed reactions in the original computational model [15]. Later in [17], it was showed that independently of the magnitude of delays, the positive steady state is always unstable and the model can not reflect a stable behaviour of newly formed vessels observed for less aggressive tumours, [18]. In [19], combining the ideas presented by Hahnfeldt et al. [3] and Agur et al. [20] a model of three differential equations with discrete delays describing the process of formation of new vessels that could reflect both, stable and unstable, structures of vessels observed in reality was proposed. Next, in [21] a detailed stability analysis of this model was performed showing that for some parameters a hysteresis loop can be observed and multiple stability switches can occur with increasing time delay.
The models described above are mean-field, in the sense that they deal with the dynamics of the tumour without taking into account its spatial organisation. On the other hand, there are several models studying this spatial dependence focusing on the patterning of the vascular network without delving into the microscopic details [22][23][24][25].
In the current work we aim to bring together the microscopic reactions occurring at cellular level with the mean-field models looking for the essential features of the non-linear interaction between cells, vessels and endothelial growth factors, with the goal to build a solid understanding that can be expanded into the spatial models.
Such a microscopic description starts by identifying the system of reactions that governs the tumour growth. The advantage of having the microscopic description is twofold. On one hand, a formal derivation of the macroscopic system of ODE and the meaning of the parameters become available, giving a justification for the type of models considered in [19] and [21].
In contrast, the microscopic description is intrinsically stochastic, formally described by a multivariate birth-death process [26], that can be decisive in the macroscopic tumour dynamics giving place to effects not present in a pure deterministic description [27][28][29], and feature an upcoming effort to integrate them into experimental data [30,31]. For this reason we also analyse the stochastic system through exact simulations of the kinetic reactions and by using the corresponding Chemical Langevin Equation [32] that gives analytical insight on the effect of intrinsic noise in the macroscopic evolution of the tumour.
Additionally, to give soundness to the prediction of the model we use experimental data to give meaningful values to each parameter of the model. The result is a working model able to predict the dynamics of the Lewis lung carcinoma.

Microscopic description
The dynamical evolution of the tumour is determined by the non-linear interactions between the tumour cells, N, and the effective number of blood vessel cells, V, or to be precise the number of endothelial cells that lines the interior surface of blood vessels. This interaction takes place by means of VEGF, the regulating proteins secreted by the tumour cells that promote the angiogenesis process, P. It is also useful to define the relative quantity E V/N, which measures the number of vessel cells that are available per tumour cell i.e. gives a measure of the amount of nutrient each tumour cell receives.
The different events characterising the dynamics of the system take place at a microscopic level as a set of events changing individually the population of each species. Following the molecular analogue, such events will be referred here on as reactions. Some of these reactions depend on the absolute number of tumour cells/vessels/proteins (denoted by N/V/P, respectively) while others depend on the concentration of each of those species in the system where O is the parameter relating concentrations with the number of cells (or molecules), giving a measure of the size of the system. The larger is O for a certain concentration, the greater will be the corresponding number of cells (proteins). The parameters used to describe the particular reactions are summarised in Table 1.
Cell reactions. The first considered reaction is the cell proliferation in which a cell can divide giving place to two daughter cells, where α > 0 is the maximal proliferation rate and g b (Á) is a decreasing function capturing intracellular competition for resources. Such a competition will increase with tumour cell concentration N Ã , causing a reduction in the proliferation rate. On the other hand, the proliferation capability increases with the amount of nutrients that the tumour cells receive through the vessels f 1 (E). These assumptions capture a switching behaviour in which the cells proliferate if the environmental conditions are favourable (N Ã < K + f 1 (E)) ( Fig 1A). Thus, for a large concentration of available nutrients, the cells duplicate with a maximum rate α, (g b (0) = 1), while the division decay for a lower concentration of nutrients (lim x ! 1 g b (x) = 0). In particular, we assume that g b (Á) has the Hill form with n 1 > 0. The actual concentration threshold (K + f 1 (E)) depends non-linearly on the concentration of endothelial cells through E and it is related to the maximal tumour size that can be reached assuming the effective vessel density is equal to E. The parameter K corresponds to the maximal tumour size if no vasculature is present. Thus, following the ideas of [19][20][21], it is assumed that f 1 (E) is an increasing Hill function of E such that f 1 (0) = 0. Additionally, in a nutrient saturation environment there is a maximum response b 1 .
On the other hand, under hostile conditions cells can also die following similar arguments of competition and saturation of resources as the growth rate, where δ N is the maximum death rate and function g d (Á) grows with the concentration of tumour cells and decreases with the nutrient supply, and analogously to [19][20][21], is assumed a Hill increasing function, (see Fig 1C). Note that this description assumes the same characteristic tumour size K + f 1 (E) for both, the birth and death rate, as a first approximation in order to reduce the complexity of the model. Protein reactions. The regulation proteins, that stimulate vessel growth, are secreted by tumour cells, The rate of stimulatory protein secretion f 2 (E) also depends on the supply of nutrients to the tumour cells (E) and it is assumed to be a non-negative decreasing function of E that satisfies (lim E ! 1 f 2 (E) = 0), see [20]. On the other hand, if no vascular network is present, the VEGF proteins are secreted at a maximum rate a 2 (f 2 (0) = a 2 ) with a characteristic concentration c 2 , giving place to the used expression for (see Fig 1B). Additionally, the proteins degrade with a constant rate δ P , Vessel reactions. Finally, the regulatory proteins secreted by the tumour cells control the angiogenesis process, thus Similarly, to the tumour cell proliferation, the vascular generation/death is controlled by two switching mechanisms that follow a rate proportional to the number of vessel cells and a factor dependent on the protein concentration reaching maximum generation/death rates b 3 , a 3 > 0, and, following the ideas of [20], can be described as (see Fig 1B), where n 3 , n 4 > 0 and m 3 is a characteristic protein concentration under which the production and degradation of vessels reach half of their maximum values. Specifically, in the case b 3 = a 3 , m 3 is the protein concentration for which vessel cell birth and death are in balance. More detailed vessel birth and death reactions (9) and (10) would include competition between vessels for the available protein. This is not included in this model as a first approximation to vessel cell dynamics. Table 1 gathers the introduced parameters for the three species and their interpretation.

Macroscopic system
Once the reaction scheme is set, a coarse grained (macroscopic) approximation can be performed by expressing the resulting temporal dynamics as a system-size expansion. The leading order term gives the macroscopic behaviour, while the subsequent orders capture the stochastic nature of the system (for details, see [33]). The resulting macroscopic dynamics follow Systems (13)-(15) describes the evolution in time of the three magnitudes defining the state of the tumour (tumour cells, proteins and vessel cells) in such a way that any tumour (finite N, P, V ! 0) will evolve to a steady state (see Theorem A in S1 Text). Thus, to determine and classify the tumour growth dynamics is mandatory to find the steady states of the system . Different steady states can be associated with different tumour morphologies with differentiated biological consequences.
Macroscopic tumour morphologies. The steady states of Eqs (13)-(15) will be labelled alphabetically: A, B and C i . The trivial steady state always exists, and corresponds to the case in which the tumour disappears completely. This state is always unstable (see Lemma D in S1 Text). In the rest of steady states (N Ã 6 ¼ 0), i.e. the final state of the tumour will always have a finite number of tumour cells so the relative measure E = V/N is mathematically well defined. To study the existence and stability of these solutions comes in handy to rewrite Eq (15) in terms of E without any loss of generality In order to have a finite tumour (N Ã 6 ¼ 0), a balance of production and death of tumour cells is required (see Eq (13)) where the solution γ is defined uniquely because g b is a decreasing positive function, g d is increasing non-negative, g d (0) = 0, and g b (1) = 0. In particular, if n 1 = n 2 , then γ is given by the following analytical expression Similarly, in order to reach the steady state condition dE/dt = 0 in Eq (16) requires E = 0 or a balance in the vessel production and degradation through f 3b (P Ã ) = f 3d (P Ã ). If E = 0, from Eq (17) we obtain N Ã = Kγ, since f 1 (0) = 0. Using Eq (14) we get P Ã = a 2 N Ã /δ P , obtaining the steady state This solution corresponds with the case in which the vessel network disappear from the tumour, while the number of tumour cells reach a constant amount, i.e. we have an avascular tumour. The steady state B is locally asymptotically stable if the condition f 3b (a 2 Kγ/δ P ) < f 3d (a 2 Kγ/δ P ) holds; for details see Proposition E in S1 Text.
Finally, it only remains to consider the steady state corresponding with E 6 ¼ 0, and the vessel balance f 3d (P Ã ) = f 3b (P Ã ). This last equation has exactly one positive solution denoted byP Ã (see of Fig 1D). Moreover, in the particular case a 3 = b 3 we haveP Ã ¼ m 3 . On the other hand, solution to Eq (17) gives N Ã = γ(K + f 1 (E)). Introducing this expression of N Ã into Eq (14) we obtain This equation can have multiple solutions, determining the steady states C i . The existence of these steady states correspond with the positive zeros of the auxiliary function For the functions f 1 and f 2 given by Eqs (3) and (7), respectively, there can be at most three solutions to Eq (18) (see Fig 2 and Proposition B in S1 Text). The stability of a steady state C i ¼ ðN Ã C i ;P Ã ; V Ã C i Þ depends on the behaviour of the function h given by Eq (19) at the point Fig 2). Otherwise, it is locally asymptotically stable (black circles in Fig 2) (See Proposition G and Remark H in S1 Text for details).
The number of stable states depends on the sign of h(0) and the local maximum and minimum values of h(E), in a generic case, having four different outcomes gathered in Fig 2: No steady state (D); one stable steady state (A); two steady states, one stable and one unstable (C); three steady states, two stable and one unstable (B). In all of them the stability alternates with increasing E, being the largest E C i always stable. Changing the parameters of the model one can change the shape of h(E) or shift h(E) vertically, thus changing the number and positions of steady states C i . This is illustrated in Fig 3 where it is indicated the existence and stability of the steady states C and B in dependence on the parameter δ P .
In Table 2 the existence and stability results of the steady states of the systems (13)-(15) are summarised

Stochastic analysis
Dynamical systems (13)-(15) correspond to the macroscopic limit in which the population number is big enough to neglect the intrinsic fluctuations of the system. Nevertheless, there is always a finite population number in a real tumour and the study of stochastic effects is necessary to understand the biological scenario. Fluctuations appear naturally in the definition of the microscopic system (Eqs 1-10).
Stochastic description of the model. The exact stochastic dynamics describe a multivariate birth-death process [26] that can be tackled analytically by means of the corresponding Master Equation [33,34]. Unfortunately, the analytical solution of the Master Equation exists only for a few simple situations, e.g. when birth and death are constants and the corresponding Master Equation is linear. This is not true for the current case and numerical simulations of the birth-death process are required. This can be achieved by means of the Gillespie algorithm [35], which generates exact stochastic realisations of the microscopic reactions.  (13)-(15) for the tumour cell density N* (A) and the protein density P* (B) as a function of δ P . In the interval δ Pcr,1 < δ P < δ Pcr,2 (for more details see Proposition B in S1 Text) a coexistence of the stable steady states is observed. Stability is indicated with solid lines (stable solutions) and dashed lines (unstable solutions).   (13), (14) and (16) Even though the system can be simulated exactly using Gillespie algorithm, when the number of reactions is high, the simulations can be computationally heavy. Additionally, such simulations do not provide analytical insight into the magnitude of the fluctuations. For this reason it is often useful to approximate the stochastic process, for instance, through state-space truncating procedures [36] or by approximating the fluctuations as a superimposed random diffusion [32,37]. This latter case includes the Chemical Langevin Equation (CLE), which introduces the fluctuations as a multiplicative noise to the macroscopic evolution of each one of the species in Eqs (13)- (15). The intensity of the noise is determined by each reaction channel rate and its stoichiometry [32][33][34]. This leads to the following system of stochastic differential equations By ξ i (t) we denote white Gaussian noises with the following properties where δ(t−t 0 ) denotes the Dirac delta distribution and δ ij is Kronecker's delta. The different noise terms in Eq (20) are uncorrelated between themselves because each reaction only involves changes in one of the species. The CLE is valid as long as a time scale exists such that all the reactions take place without a relevant change in their reaction probability. This is similar to other expansions of volume O in which the stochastic effect is reduced to an order O 1/2 , see e.g. [32,33]. Note that the deterministic limit (Eq 13) is recovered for O ! 1.
In the current analysis we use both, the Gillespie algorithm and the Langevin description. The first returns exact trajectories while the second is used to determine analytical properties or obtain stochastic trajectories when the populations of the species is large enough. Actually, later in the text, the Langevin system will be tested for small populations to assess the agreement of both descriptions obtaining that the CLE approximation is enough in most biological relevant tumour scenarios. Table 2. Steady states of the macroscopic systems (13)- (15).

Steady state Stability
A = (0, 0, 0) unstable B = (Kγ, a 2 N*/δ P , 0) stable for f 3b (a 2 Kγ/δ P ) < f 3d (a 2 Kγ/δ P ) unstable for f 3b (a 2 Kγ/δ P ) > f 3d (a 2 Kγ/δ P ) Behaviour of the stochastic system. In contrast to the deterministic model, the resulting stochastic trajectories have the ability to explore the whole dynamical landscape. Concretely, one of the effects of the intrinsic noise is the ability to jump between different co-existing stable steady states. This is presented in Fig 4 where the used parameters correspond with previous studies [19,21] for which the bi-stability behaviour (see Fig 5 and Remark H in S1 Text) is observed, (see S1 Fig). This switching is not available for large O (O > 1000), where the trajectories follow probability distributions around the corresponding deterministic solution. For smaller volume sizes these transitions become more frequent and the noise plays a more important role determining the final state of the tumour (see Figs 4,6,7 and S1 Fig). This effect is also seen as a broadening of the probability distributions around each stable state, providing variability to the tumour characteristic even when it can be described within a single morphology. Additionally, as expected, for very small sizes (O * 10) the Langevin approximation can not recapitulate results predicted by the Gillespie algorithm. Concretely, the CLE overestimates the fluctuations giving a higher switching probability. Nevertheless, as will be discussed in the next section, in a biologically relevant scenario O ) 10. Finally, it is interesting to point out that the transition between states for a growing tumour happens usually at early times when the populations of species is small (Fig 4), even though there is a continuous flux over time that predicts a complete depletion of the macroscopic state for a long enough time (see S1 Fig). A precise quantification of this effect in which jump rate between states changes in time is analysed in the Results section.
Fluctuations also give the system the possibility to visit unstable states such as state B for which V Ã = 0. Nevertheless, once the number of vessel cells reaches zero, the tumour is unable to create more vessels Eq (9). Therefore, the number of vessel cells remains zero during the rest of the evolution of the tumour, describing an avascular state. This can also be predicted from Eq (20), where for V = 0 the deterministic and stochastic terms become null. Thus, even though the state B can be unstable, it will always capture a stochastic trajectory that reaches V = 0 in the absorbent set of states (N, P, 0).
Likewise, the state A = (0, 0, 0) (all the tumour cells disappear) is an absorbent state. Even though state A is linearly unstable, once a fluctuations drives the system in state A cells cannot proliferate and the system is trapped i.e. fluctuations can make the tumour disappear. Again, the absorbent state A, which is always unstable, has a very different dynamics when the stochastic nature of the tumour is considered.
Such behaviour is shown in Fig 8, where a stable tumour becomes avascular and finally dies due to the fluctuations. Since this is an stochastic effect, the rate of absorbency will decrease with increasing O and disappear in the deterministic limit (Figs 6 and 7).

Estimation of the parameters for the Lewis lung carcinoma
The macroscopic models (13)-(15) gives a general description for the evolution of different tumours compatible with the microscopic framework proposed. Nevertheless, the model predictions will vary for different tumour physiologies that come determined by the different parameters that describe the microscopic reactions.     In particular we apply our model to the lung carcinoma (Lewis) mouse tumour for the tumour cell line from [38]. An optimal estimation of model parameters should make use of experiments under identical conditions e.g. run by the same lab. Unfortunately, this information is not available, and parameters had to be extracted from different experimental set-ups and sometimes even different mouse lines. Quite often, such experiments only return qualitative data, such as the case of luminescence measurements, and quantitative readouts are only relative.
Under these circumstances, the parameter estimation is driven in two parts. First, by exhausting all the parameters that can be obtained directly from experimental measurements, this first estimation consists of measurements on the VEGF and the tumour cells physiology. This will give a working range of parameters where the final parameter values can be inferred by fitting experimental data on tumour growth to the dynamics predicted by the model.
Since the experimental data provided report the progression of tumour volume (e.g. see [38]), it comes in handy to rewrite the populations in systems (13)- (15) in terms of absolute numbers instead of densities, by rescaling the magnitudes with the parameter O without any loss of generality, Estimation of VEGF parameters. The VEGF dynamics comes mainly determined by two parameters, namely, its maximal production rate a 2 and its degradation rate δ P . Nevertheless, VEGF interacts in general through three major isoforms of mouse VEGF (VEGF 120 , VEGF 164 and VEGF 188 ), which have different expression levels depending on the species and tissue type [39] and [40]. For mouse skeletal muscle, the VEGF isoform has been reported to have mRNA expression ratio of (VEGF 164 +VEGF 188 ) to VEGF 120 has been measured to be 92:8, see Results in [40], while for the lung tissue this ratio is 82:18. The VEGF isoform distribution and dynamics has already been studied in a thorough in silico model for the mouse skeletal muscle [41], where the secretion rate of VEGF was estimated to be equal to 5875.2 molecules per cell per day; see Results in [41]. Additionally, the lower and upper bounds of the secretion rate of VEGF 164 +VEGF 188 has been found to be 864 and 17000 molecules per cell per day, respectively. Thus, using the findings on the skeletal muscle mouse with the ratio for lung tissue we can set bounds for total VEGF production rates to be 939 and 19000 molecules per cells per  day. Since the estimations made in [41] where done extrapolating information from skeletal muscle tissue to lung tissue we will assume a looser range that will be refined later in the parameter fitting. On the other hand the VEGF degradation rate δ P of the three isoforms can be considered comparable and determined from in vitro half-life measurements [42,43], obtaining a degradation rate of 16.6 day −1 . Again, this datum correspond to a different tissue, in this case human cells. Thus, we treat that estimation as a reference for the mouse tumour.
Estimation of tumour cells growth parameters. The maximal tumour growth rate i.e. the tumour cell proliferation rate with nutrient profusion has been reported for the lung carcinoma (Lewis) C57BL/6 [44] showing a doubling time of 1-1.7 days. This gives a range for the maximal growth rate parameter α in the range 0.41-0.69 day −1 .
On the other hand, in the absence of nutrient supply from blood vessels the growth of the tumour is hindered reaching a maximal number of tumour cells that in our model comes described by the parameter product g K . This can be related to the first stages of cancer disease when cells are supplied only by diffusion and no angiogenesis takes place. The typical diameter of the tumour that can be achieved in this avascular phase of growth is 2-3 mm, [1,2]. Since the average diameter of the Lewis lung cancer (LLC) cell equal to 7 μm [38], the average number of cells at the end of avascular growth can be estimated to be in the range 2.33 × 10 7 -7.87 × 10 7 (cells).
All the estimated parameters above mentioned, their ranges and their bibliographical sources are summarised in Table 3.
Note, that due to a lack of experimental data, none of our estimations take into account the number of vessel cells V. For instance, in [38] the final vessel cells density for LLC tumours was determined via counting the number of the microvessels per high-power field (HPF) within hot spot area. Similarly, in [45] for the same type of tumour, the number of microvessels were counted for each section and the mean number of microvessels/HPF were calculated as a marker for microvessel density in tumour tissue. In both cases, such data can not be used to interpret the actual number of cell vessels V. Similar problems are found in measures of VEGF such as Fig 5B in [45], where the concentration of plasma level of VEGF (expressed in pg/mL) at single time point was measured. Such blood samples were drawn by a terminal cardiac puncture after opening of the chest wall and collected in heparin anti-coagulated tubes preventing to use those data as a measure of VEGF concentration in the tumour tissue.
Fitting of the parameters using tumour volume growth data. With the parameter ranges obtained, a further inference can be carried out by using volume data growth for the Lewis lung cancer (LLC) published in [38]. Again, assuming spherical tumour cells of radius r = 3.5 μm the volume measurements of LLC cultivated in vivo can be directly converted to number of cells. The resulting profile of number of cells in time can be compared with realisations of the model for different parameter values. The optimal values of the parameters can be found by optimising a distance function between the experimental data and the model prediction. In this case we minimise the mean square error (MSE) between both where M is the number of experimental measurements, x data i are the measured values and x sym i the values of the solution of ODE for respective time points. Due to the number of parameters, the optimisation has been done using Particle Swarm Optimization (PSO), a stochastic optimisation technique that explores efficiently the parameter space [46,47]. The PSO was carried out using the Matlab 1 implementation of the algorithm.
The resulting optimal set of parameters shows a good agreement with the experimental data (MSE = 0.07%) predicting LLC volume growth during the first months (Fig 9B). The corresponding values of the estimated data are given in Table 3. In order to give soundness to the parameter estimation, a sensitivity study was performed on the inferred parameters. The results of the sensitivity analysis are discussed below in the Results section, while computational details and specific results are gathered in S2 Text.

Results
The fitted model predicts a single stable state C coexisting with unstable states A and B. Thus, the model suggests that the natural growth of the tumour results as the evolution towards the state C, which describes a large tumour with well developed vasculature. This dynamic shows a very rapid secretion of the VEGF proteins within few first hours, that is followed by the fast increase of the vessels network (number of endothelial cells) within few first days of the experiment, see Fig 9. After the first day the active VEGF returns to lower levels. This active VEGF stabilises during the first week causing the future tumour growth. The fast dynamics of VEGF within the first week can be explained by the injection of a large tumour that need to develop vascular network fast.
The performed sensitivity analyses (see S2 Text) shows that after homogeneous perturbations on the parameter set (±2% over each parameter), the resulting perturbed tumour morphology is more sensitive to certain specific parameters. Parameters δ P and b 3 dominate on determining the number of tumour cells N and effective vessel E, whereas m 3 and a 3 are also important for determining the number of regulating proteins P. It is interesting to point out that the sensitivity of P to different parameters strongly changes during the evolution of the tumour. On the other hand, the analysis shows that K and tumour cell death rate δ N are the parameters with less influence i.e. the ones with correlation coefficient in the interval (-0.2, 0.2). Including the initial conditions P(0) and E(0) in the sensitivity analysis (N(0) is given by experimental conditions) did not change the sensitivity results (results not shown), suggesting that the results are robust to small perturbations of the initial conditions.

Intrinsic fluctuations determine the survival of the tumour during the first week
The stochastic dynamics of the tumour reveals that the characteristic evolution predicted by the deterministic equations is preserved. Nevertheless, the fluctuation in the population of species can induce on occasions the death of the tumour as described in Materials and Methods where the tumour dynamics is absorbed by the trivial state A. For the estimated parameters, the survival of the tumour is determined during the first week when the population of the different species is small (see Fig 10A, 10B and 10C). On the other hand, all the tumours that survive this first week, grow considerably unaffected by the stochastic effects (see Fig 10D, 10E  and 10F).
A close analysis of the first days of the tumour growth shows that the noise has almost no influence on the dynamics of the model during the first two days (see Fig 11). At the beginning of the second day the noise increases dramatically with two possible outcomes. Either the tumour keeps growing overcoming the noise effects, or a fluctuation drives the tumour to the absorbent state A. An analysis of the noise intensity of the CLE Eq (20) can be evaluated analytically along the deterministic trajectory (see Fig 11B). This result points out how the relative intensity of the noise can be neglected during the first two days and grows suddenly until day four were it decays again for the rest of the evolution of the tumour, setting a time window that determines the survival of the tumour. During this time window, around 11% of the growing  [38] (red circles) and solutions of model (23) for δ P = 8.9112 day −1 (blue lines) for the best fitted parameters, with MSE error equal to 0.07%. (C), (D), (E) and (F), comparison of solutions P*, E and V* of system (23) with δ P = 8.9112 day −1 (blue lines). Solid green lines show the behaviour of the model for a high suppression of VEGF (δ P = 2 × 10 4 day −1 ). For δ P > δ P cr,2 the model has only one stable steady state with a small amount of tumour cells. Rest of parameters are given in Table 3. doi:10.1371/journal.pone.0155553.g009 Angiogenesis Model for Lewis Lung Carcinoma tumours die spontaneously. This quantity remains practically constant along the rest of the life of the tumour that shows less than a 13% of spontaneous death at day 50 (see S2 Fig). A more detailed analysis (see Remark J in S1 Text) shows that the noise intensity of the tumour cells increases with the effective amount of nutrient f 1 (E), revealing that the increase in the intensity of the noise between days 2 and 4 is the aftermath of the fast vascularisation predicted during the first days of the tumour growth.

VEGF degradation rate controls the appearance of dormant tumours
For the parameters inferred, there is only one stable state (C). Nevertheless, a further analysis on the dependence of the steady states with biological meaningful changes in the parameters shows that new steady states C i can arise, for details see S1 Text. This is the case of VEGF degradation δ P that can be modified through the administration of a drug (e.g. bevacizumab commercial name Avastin), a humanised monoclonal antibody, that binds the vascular endothelial growth factor A (VEGF-A) and thus inhibits VEGF-A binding to receptors on the surface of endothelial cells, effectively increasing δ P .
A stability analysis on δ P reveals that for higher values of δ P new stable states C i arise. Concretely, over a threshold δ Pcr,1 ' 160 day −1 two stable steady C 1 , C 3 coexist (see Figs 3 and 12). The new stable tumour morphology C 1 is considerably smaller. Nevertheless, the evolution to this state requires to increase δ P beyond δ Pcr,2 ' 1.2 × 10 4 day −1 where the original state C 3 disappears resulting in a single stable steady state with a comparably smaller amount of tumour cells (see Fig 12 and green lines in Fig 9).
The size of the new tumour morphology C 1 is 8 × 10 7 , corresponding with the number of tumour cells for the maximal tumour size ( K g) that can be achieved without direct nutrient supply from the blood vessels. This is due the fact that for sufficiently large δ P , the number of cells N of steady states C 1 and B merge (see Fig 12). This result suggests that large enough inhibition of VEGF activity, induces dormant tumours with diameter around 2 mm. Moreover, the effective vessel network (E) is of order 10 −4 , that is around nine orders of magnitude lower than the one for the untreated value of δ P = 8.9 day −1 , preventing active support of the cancer cells in nutrients and in turn also preventing the fast tumour growth (see Fig 12). Additionally, the stability analysis (which results are summarised in Table 2) shows that the number of proteins in the steady state does not depend on δ P for the positive steady states C i . For the biologically inferred parameters the number of proteins is low and close to 100 not causing strong stimulation of the endothelial cells, as in the steady state efficient vessel network has been already created and intensive stimulation of the endothelial cells is not needed any more.

Summary and Discussion
In this paper we have investigated the microscopic foundations of the types of model of tumour angiogenesis proposed in [19,20] and later studied in [21]. Our model predicts the same  Table 3. The stability of existing steady states in indicated by solid lines, while instability by dashed. Colours indicate the steady states B (green) C i (blue). Additionally, E = 0 for state B, therefore it is not indicated on the righthand side panel. qualitative behaviour (regarding the existence and stability of multiple steady states that show hysteresis), being able to reflect the possible instability of the blood vessel formation and structure that is observed in the experiments, see [14,18], but without the requirement of temporal delays. This was done by describing the tumour dynamics from a set of reactions from which the deterministic macroscopic ODE system was derived. This macroscopic solution predicts different behaviours of the tumour, namely, the tumour disappears (state A), the tumour reaches an avascular state (vicinity of state B) or the tumour evolves towards a discrete set of tumour morphologies (states C i ). This kind of behaviour is not observed in a family of classical Hahnfeldt et al. [3] angiogenesis models without delay where all solutions stabilise at a unique positive steady state. On the other hand, in a modified Hahnfeld et al. model, proposed in [48], there can exist multiple positive steady states, but due to a lack of the hysteresis loop, the solution cannot jump between stable steady states when parameters change.
The microscopic description opens the door to the study of the fluctuations intrinsic to reactive systems. This stochasticity allows the transition between the different tumour steady states, providing new behaviours such as reaching the avascular absorbing set of states (N, P, 0) or the spontaneous death of the tumour.
In order to explore further the predictions of the model we applied the model to the Lewis lung cancer. The parameters were inferred using the available experimental data in the literature obtaining a good agreement with experimental observations on the tumour growth dynamics in time, during the different phases of tumour formation. A study on the intrinsic noise effects on the tumour dynamics shows that the noise effects are of key importance during the first stage of the tumour formation when fluctuations can destroy the tumour. This stochastic time window of extinction is dictated by the fast vascularisation predicted during the first days of the tumour growth that could be readily tested by a careful tracking of the tumour cell population during this time period. Alternatively, the fast vessel growth could be a consequence of the lack of intervessel competence in Eqs (11) and (12), posing the necessity of such a mechanism for a proper tumour description. In both cases, a deep understanding of this effect and the dependence of this time window with the tumour physiology could provide insight into to tumour prevention and treatment.
We also studied how variations in the degradation of VEGF regulate the tumour progression. The biological in vivo degradation rate only predicts one possible tumour state. Nevertheless, the increase of the degradation rate introduces a new stable dormant tumour state that is orders of magnitude smaller. This is biologically relevant since the degradation rate can be easily modified with existing drugs. These two mechanisms of tumour control, spontaneous dead and dormant tumour, are just two cases in which this model helps to determine changes in the tumour morphology that could be addressed experimentally. Nevertheless a more exhaustive analysis may reveal novel strategies in oncology. We hope that this model inspires new experiments and measurements that push forward the knowledge of the field.
Overall, the presented model introduces a simplified reaction scheme being able to exhibit measurable dynamics of tumour growth and vessel formation. Such a reduced model has proved to be complex enough to accommodate experimental data. On the other hand, one can think of a more deeper description of the angiogenesis process described, considering, for instance, vessel cells competition or more complex regulatory functions. However, such approach would require more quantitative data than the one available for the Lewis lung carcinoma. Alternatively, extensions of the model could also address new biological questions such as a spatial description of the angiogenesis process together with the influence of additional vessel growth stimuli factor (e.g. some kinds of interleukins). Moreover, we hope that the general description of the model makes the analysis easily transferable to other types of cancer motivating the use of experimental data in, otherwise, purely theoretical models.
Supporting Information S1 Text. Mathematical derivation details. Theorems and proofs concerning the existence, uniqueness, boundedness of the solutions; as well as the existence of the positive steady states and conditions for the stability of the steady states for the macroscopic model. (PDF) S2 Text. Sensitivity analysis of the estimated parameters. Detailed information and results of the sensitivity analysis of the estimated parameters in Table 3