A mathematical model of systemic inhibition of angiogenesis in metastatic development
S. Benzekry^{*}^{*}*Corresponding author. Email: sebastien.benzekry@inria.fr, A. Gandolfi, P. Hahnfeldt
Inria team MC2, Institut de Mathématiques de Bordeaux, Bordeaux, France
Istituto di Analisi dei Sistemi ed Informatica “Antonio Ruberti”  CNR, Roma, Italy
Center of Cancer Systems Biology, GRI, Tufts University School of Medicine, Boston, 02142, USA
Abstract. We present a mathematical model describing the time development of a population of tumors subject to mutual angiogenic inhibitory signaling. Based on biophysical derivations, it describes organismscale population dynamics under the influence of three processes: birth (dissemination of secondary tumors), growth and inhibition (through angiogenesis). The resulting model is a nonlinear partial differential transport equation with nonlocal boundary condition. The nonlinearity stands in the velocity through a nonlocal quantity of the model (the total metastatic volume). The asymptotic behavior of the model is numerically investigated and reveals interesting dynamics ranging from convergence to a steady state to bounded nonperiodic or periodic behaviors, possibly with complex repeated patterns. Numerical simulations are performed with the intent to theoretically study the relative impact of potentiation or impairment of each process of the birth/growth/inhibition balance. Biological insights on possible implications for the phenomenon of “cancer without disease” are also discussed.
Key words: Nonlinear renewal equation; Cancer modeling; Metastasis; Angiogenesis inhibition;
AMS subject classification: 35B40, 35Q92, 92C50
Contents
1. Introduction
Metastasis, the process by which secondary tumors are shed from a primary lesion to colonize local or distant sites, is a complex process that is responsible for a very large proportion (90%) of deaths by a cancer disease [15]. Despite significant efforts in the last decade to strengthen our understanding of metastatic cancer biology, global mechanisms are still poorly understood [23].
Among the wide variety of topics to be addressed, we focus here on signaling interactions within a population of tumors that impact on the global dynamics of the system, at the organism scale. Indeed, it has been long known that several tumors simultaneously growing within the same host influence each other, mostly in an inhibitory fashion. Most of the experimental studies were conducted in twotumors experimental systems where impaired growth was observed for a second inoculum in the presence of a preexisting tumor, a phenomenon that has been termed concomitant resistance [30] (see [11] for a review). This distant inhibition not only occurs between two artificially implanted tumors but also between a primary tumor and its metastases [13, 19]. Such an observation is of fundamental importance for cancer biology as it impacts on the temporal development of the disease, but also has clinical implications in terms of metastatic dormancy [1] and surgery [31].
Several hypotheses have been proposed to explain this phenomenon, among which athrepsia (deprivation of nutriments for the second tumor due to monopolization of nutrients by the first one) or immune enhancement. Indeed, it is known [23] that some tumors are immunogenic, i.e. that they provoke a hostile reaction from the immune system. This reaction could be triggered by the presence of a first tumor and suppress the growth of a subsequent implant. However, the importance of this last hypothesis was tempered in the 1980’s when concomitant resistance was shown to occur in immunodeficient mice [20].
In the 1990’s, Judah Folkman and his team put forward a novel hypothesis, based on their discovery of endogenous inhibitors of tumor neoangiogenesis [27, 28]. Angiogenesis is the process of creation of new blood vessels from an existing vasculature and was shown to be of fundamental importance for a tumor’s development [18] in order to ensure access to nutrients. Indeed, in the absence of angiogenesis, a neoplastic lesion cannot grow beyond 23 mm in diameter [18]. Observation that a tumor not only stimulates this process (through production of growth factors) but also regulates it by producing angiogenesis inhibitors opened the way to a new explanation of concomitant resistance. A primary tumor could distantly inhibit a secondary tumor through inhibitory angiogenic signaling, an hypothesis that is strengthened by the observation that angiogenic inhibitors have significant longer halflives than angiogenic stimulators (of the order of hours for inhibitors such as angiostatin [27] or thrombospondin1 [32] and of the order of minutes for stimulating agents such as VEGF [16]). The formers are thus transferred to the central circulation and from there systemically distributed to distant sites. Considering the large and unequivocal body of support for the role angiogenesis inhibition plays in the maintenance of tumor dormancy [24, 27, 28, 32, 36, 33] we will focus on systemic inhibition of angiogenesis (SIA) as the major process for tumortumor interactions at the organism scale.
2. Model
The main idea, originally introduced by Iwata et al. [25], and further studied in [3, 12] is to represent the population of metastatic colonies as a density structured in tumor macroscopic traits (such as the volume in [25]) and to derive a transport partial differential equation on that reflects mass conservation during the growth, endowed with a nonlocal boundary condition for dissemination (birth) of new metastases.
A major limitation of this initial model for our purpose is that it does not take angiogenesis into account, although this is a fundamental process of tumor development. By combining the approach of [25] with the model of [22] for tumor growth under angiogenic control, we developed in previous work a new model taking into account the angiogenic process in the growth of each tumor [4, 5, 7]. Combined to the fact that it is written at the level of the organism makes it an adapted framework for modeling SIA viewed as inhibiting interactions across a population of tumors.
The model we propose here adds an inhibitory component to the two main processes of the previous model (angiogenic growth and metastatic dissemination). We add a new variable representing the circulating concentration of an endogenous angiogenesis inhibitor standing for all possible inhibitory molecules (examples being endostatin, angiostatin and thrombospondin1). It impacts on the growth of each tumor. As a general modeling principle, we want to be parsimonious and describe the major dynamics of the system with as few parameters as possible.
Tumors are seen as individuals whose state is described by two traits: volume and carrying capacity , the latter representing environmental limitations that constrain growth of the tumor. This carrying capacity is here assumed to be equivalent to the vascular support. The primary tumor state is denoted while the modelâs main variable for representation of the population of secondary tumors is denoted by and stands for the physiologically structured density of metastases with volume and carrying capacity at time . We assume that growth of the tumors is governed by a growth rate denoted , that spread of new metastases is driven by an volumedependent dissemination rate denoted and that the repartition of metastases at birth is given by a measure . The precise expressions of these coefficients shall be described below. We consider some fixed final time and a physiological domain where the distribution of metastases has its support, which means that only metastases having volume larger than some minimal volume and nonnegative carrying capacity are considered biologically relevant. In the formula below, stands for the external normal to the boundary of the domain . The notation stands for the subset of the boundary where the flux is pointing inward, i.e. where . The function denotes the initial distribution of the metastatic colonies. Overall, the model writes
(2.1) 
2.1. Tumor growth and systemic angiogenesis inhibition
We assume that all the tumors (primary and secondaries) share the same growth model and parameters. We do so in order to reduce the number of parameters and focus on the dynamical properties of the system. Our population of tumors should thus be viewed as a population of identical entities in mutual interactions parallel to global development. The growth velocity of each tumor is given by a vector field . Following the approach of [22] we assume
In the previous expression, the first line is the rate of change of the tumor volume (where is a constant parameter driving the proliferation kinetics) and the second line is the rate of change of the carrying capacity . The main idea of this tumor growth model is to start from a gompertzian growth of the tumor volume (that could be replaced by any carrying capacitylike growth model, see [14]) and to assume that is a dynamical variable representing the tumor environment limitations (here limited to the vascular support) changing over time. The balance between a stimulation term and an inhibition term denoted (assumed here to depend on the global state of the system represented by the density ) governs the dynamics of . For the stimulation term we follow [22] and assume
where the parameter is related to the concentration of angiogenic stimulating factors such as VEGF or bFGF. This last quantity was derived to be constant in [22] from the consideration of very fast clearance of angiogenic stimulators [16]. This fact also strengthens our assumption of a local stimulatory term, as stimulating agents don’t spread through the organism.
For the inhibition term, Hahnfeldt et al. [22] only considered a local inhibition coming from the tumor itself. Our main modeling novelty is to consider in addition a global inhibition coming from the release in the circulation of angiogenic inhibitors by the total (primary + secondary) population of tumors. The following is an extension of the biophysical analysis performed in [22]. Let us consider a spherical tumor of radius R inside the host body. The host is represented, for simplicity, by a single compartment of volume in which concentrations are assumed spatially uniform. Let be the inhibitor concentration inside the tumor at radial distance . Let the intratumor clearance be zero [22]. At quasi steady state, solves the following diffusion equation:
where is the inhibitor production rate and is the inhibitor diffusion coefficient. This equation is endowed with the boundary condition , being the systemic concentration of the inhibitor resulting from a primary tumor volume and secondary tumors density at time . Solving this equation (using that ) we obtain
From this expression we compute the mean inhibitor concentration in the tumor to obtain
where is a sensitivity coefficient. For , considering that the total flux of inhibitors produced by a tumor with volume is and assuming that the inhibitor production rate is the same in all the tumors, we have
where is an elimination constant from the blood circulation. Defining , we get
(2.2) 
endowed with zero initial condition ().
Overall, the explicit expression of the metastases growth rate is given by
(2.3) 
where and
(2.4) 
Note that we retrieve here the local term from the analysis of [22]. Our analysis results in addition of a global term for the effect of systemic inhibition of angiogenesis.
For the primary tumor, we assume the growth velocity, hence the dynamics of , to be given by
(2.5) 
where is the initial volume of a tumor and its initial carrying capacity.
2.2. Metastatic spreading
There is no clear consensus in the biological literature about metastases being able to metastasize or not [35, 9, 34]. However, we argue here that cancer cells that acquired the ability to metastasize should conserve it when establishing in a new site. Moreover, since metastasis is a long process before being detectable [35, 9, 34] (in particular because tumors could remain dormant during large time periods), the absence of clear proof in favor of metastases from metastases could be due to the short duration of the experiments compared to the time required for a secondary generation of tumors to reach a visible size. Here we are interested in long time behaviors and, although metastases from metastases could be neglected in first approximation, we think this second order term is relevant in our setting and chose to include it in our modeling, following clinical evidences of secondgeneration metastases [2].
Successful metastatic seeding results from one malignant cell being able to overcome various adverse events including: detachment from the tumor, migration in the local microenvironment, intravasation in local blood (or lymphatic) vessels, survival in the circulation, escape from immune surveillance, extravasation, survival in a new environment and eventually establishment of a new colony at the distant site (see [21] for more details). We regroup here all these events into one emission rate quantifying the number of successfully created metastases per unit of time, and neglect intricate description of all these processes. We assume that very small metastases do not metastasize arguing that they do not have access to the blood circulation and hence include a threshold below which tumors do not spread new individuals. Apart from addition of this threshold, the expression of is the one from [25] and is given by
(2.6) 
where and are coefficients quantifying the overall metastatic aggressiveness of the cancer disease. The parameter represents an intrinsic metastatic potential of the cancer cells. On the other hand, parameter represents the microenvironment dimension of metastatic dissemination. It lies between 0 and 1 and is the third of the fractal dimension of the vasculature of the tumor under consideration. For instance, if vasculature develops superficially then , whereas a fully penetrating vasculature would give a value of . We chose here to take a dissemination rate only depending on the volume because simulations showed that adding a monotonous dependence on did not improve the flexibility of the model in simulations while adding at least one parameter, in opposition to our parsimony principle.
Stating a balance law for the number of metastases when growing in size gives the first equation of (2.1). The boundary condition, i.e. the second equation of (2.1), states that the entering flux of tumors equals the newly disseminated ones. These result from two sources: spreading from the primary tumor, represented by the term and second generation tumors coming from the metastases themselves, described by the term . The map stands for the volume and carrying capacity distribution of metastases at birth. We assume that newly created tumors have all the size and some initial carrying capacity and thus take
(2.7) 
i.e. the Dirac distribution centered in and refer to [6] for more detailed considerations on going from an absolutely continuous density to the Dirac mass considered here. We allow metastases to exit the domain by imposing the boundary condition only where the flux points inward. In view of expression (2.3), the growth velocity pushes tumors out of the domain when is less than . This can happen when global inhibition is strong enough such that tumors can cross the line , i.e. when . These tumors are then removed from the population, corresponding to the death of metastases caused by nutrient deprivation.
From the solution of the model, relevant macroscopic quantities can be defined such as the total number of metastases: or the total metastatic burden: .
2.3. Model nondimensionalization
We are not interested here in calibrating the model parameters or outputs to relevant biological values (see [8] for such a purpose) but will rather focus our interest on theoretical exploration of possible dynamics emerging from the model. Consequently, we rescale parameters and variables in order to essentialize characteristic properties of the system. Performing the following transformations  where is the maximal reachable volume under the Hahnfeldt model [14]  on variables
and on coefficients
gives equation (2.1) in the new variables, with new velocity given by
(2.8) 
and new differential equation on
(2.9) 
where and . Note that this nondimensionalization made parameters and disappear from the equations. We will consider the resulting model in the following but drop the tildes to simplify the notations.
3. Simulation results
We present now simulations of the nondimensionalized model (2.1, 2.52.9). For simulations based on comparison to experimental data we refer to [8] and focus here on the theoretical dynamical behavior of the model. It results from the balance between three processes: dissemination (governed by the function ) growth (controlled by the growth rate ) and systemic inhibition of angiogenesis (represented by the inhibitor quantity ). For the sake of simplicity we will reduce our analysis to one parameter per process: dissemination will be represented by parameter , growth by parameter and inhibition by parameter . Throughout the study we fixed most parameters to (see Table 1), on the notable exception of , fixed to hence assuming superficial development of the vasculature. We use this parameter set as a basis point for exploration of the parameter space. Numerical simulations were performed adapting a previously developed discretization scheme based on the straightening of the characteristics [4].
Growth 



SIA 


Dissemination 

3.1. Linear versus nonlinear dynamics
When SIA is not considered in the model (), the velocity of the transport equation does not depend on and the model is linear. It falls then in the range of classical renewal equations from structured population dynamics (see [29] for general theory) that exhibit asynchronous exponential growth governed by the first eigenvalue of the underlying operator. This means that asymptotically
where convergence occurs in norm with weight , the dual normalized first eigenvector, and is a particular first eigenvector of the linear operator (these eigenspaces have dimension 1). The eigenvalue is the unique solution of the following spectral equation
where is the flow associated to the vector field . We refer to [5] for detailed statements and proofs of results on the asymptotic behavior of the linear model and only illustrate here in Figure 1 the result and depict an example of the shape of the first eigenvector (projected on the axis).
Metastatic burden  Final volume distribution 
When considering the nonlinear model that we introduce here, nontrivial asymptotic behavior is observed with oscillations of the system (Figure 2). First, under the impulsion of the primary tumor dissemination, metastases appear and start to grow until reaching a substantial metastatic burden that in turn yields growth suppression of both primary and secondary tumors. At some point, inhibition is strong enough to push metastases out of the domain, inducing a decrease in total number of metastases and metastatic burden that translates then into lower inhibitory pressure (due to clearance of the inhibitor) allowing the metastases to regrow and restart the process.
Metastatic burden  Number of metastases 
Primary tumor volume  Systemic inhibitor 
3.2. Exploration of the growth/dissemination/inhibition balance
We explore now the parameter space for possible different qualitative dynamics, by varying 10 fold above and below the base values of the parameters respectively controlling growth (), dissemination () and inhibition (). The oscillations that were observed with the base parameter, although recovered in most of the situations, are not the only possible situation as more complex dynamics are found.
Simulation results of individual increase of each parameter are reported in Figure 3. As appears, disruption of the base regime of parameters from Table 1 (where all the forces in presence are in relative equilibrium) towards more pronounced impact of either of the constitutive processes of our model generates more complex dynamics. Moreover, different parameters have different impact on the global behavior.
Potentiation of the growth velocity (as well as resistance to the inhibition pressure) through increase of parameter results in an asymptotic behavior of the global metastatic burden which, while still being periodic, repeats a much more complex pattern, revealing interesting underlying dynamics. In particular, observation that same value of metastatic burden does not always yield same future evolution implies that no autonomous ordinary differential equation can be derived for the dynamics of , since same initial condition potentially leads to different future evolution. Indeed, when rereaches a previous value, the state of the global dynamical system is different because the composition of the tumors population (represented by ) is. Interestingly, this happens despite the fact that the growth rate depends on only through (equations (2.8, 2.9)).
Increase of the metastatic aggressiveness of the system (parameter ) results in densification of the oscillations and amplitude increase, yielding sharp repeated peaks of metastatic growth. Violent increases of the total metastatic burden are followed by similarly violent decreases that make the system reach almostzero values.
Stronger inhibition pressure delays the stabilization of the system to an oscillatory regimen, intensifies the oscillations frequency and lowers their amplitude. Note that in this situation, as well as in the base situation, the total metastatic volume remains away from zero, suggesting a nonnegligible amount of longlasting residual disease.
Turning our interest to the opposite situation, i.e. 10 fold decrease of the individual parameters, shows yet other interesting behaviors. Small value of generates an oscillatory asymptotic behavior with very low amplitude of the oscillations. From a biological point of view, this suggests the possibility of an homeostatic state of the system where all the forces in presence equilibrate to give a stable state where metastases don’t grow while still remaining present in the organism, possibly with small volumes that would make the micrometastases undetectable. Here, this homeostatic state stabilizes around the third of the maximal reachable size with a stable underlying distribution of metastases where volume of the largest metastasis is lower than 0.15. This result could suggest a possible explanation of reported cases of population of asymptomatic occult metastases with no evidence of primary lesion [37, 17, 26, 10] as resulting from mutual inhibitory interactions between tumors. This observation is substantiated by the simulation of the model with lower initial volume and carrying capacity (see Figure 5).
Small value of coefficient significantly delays emergence of the oscillations, since it takes more time to reach a sufficient amount of total tumor burden to trigger the dynamics. Lowering the value of has no clear impact on the frequency of oscillations but consistently increases their amplitude. Indeed, lower power of inhibitory pressure allows the metastases to grow to larger volumes.
Base
Large  Large  Large 
Base
Small  Small  Small 
In Figure 5 we report a few other examples of interesting dynamics arising from simulation of the model. In the case of large intrinsic metastatic potential () and inhibitor clearance set to , after an initial sharp peak of metastatic burden, we observe firing episodes of metastatic disease of increasing intensity but delayed appearance, while rest periods are characterized by almostzero amount of cancer mass. Although not completely relevant because of the nonbiological values of the parameters, biological analogy would suggest possible violent bursts of metastatic development separated by possibly long periods of minimal (and occult) residual disease that could even lead to endogenous elimination of the cancer. Indeed, under the hypothesis of strong emission of metastases and nonnegligible individual inhibitory power, it is to be expected fast exponential increase of the metastatic number of individuals, which in turn generates strong inhibition of the total population growth. However, what happens in this situation is not just acceleration of the dynamics and it remains intriguing that subsequent burst relapse with higher intensity than previous ones. A heuristic description of what happens is the following. The first burst is particular as it results from the initial condition of the primary tumor that concentrates most of the metastatic mass when it is small. This mass is eliminated all at once during the evacuation phase of the first burst that ends with smaller metastatic burden than at initiation (about in our simulation). By the time that metastatic mass recovers significative value (say for instance), production of inhibitors has occurred that changes the condition compared to the starting point and results in lower amplitude of the second burst. At the end of this second burst, much deeper metastatic burden has been reached (of the order of ) despite smaller zenith as in the first burst, because now metastases continuously outflew. With this smaller initial metastatic burden, the system had time to eliminate the inhibitor and when crosses again, it does so with smaller value of , hence producing a higher amplitude of the relapse, which in turn provokes a smaller postrelapse burden. Repeating the same mechanism explains the following bursts. It should be noted that simulating this same situation for larger times and plotting the result in logscale (Figure 5) reveals globally bounded behavior with seemingly nonperiodic orbit, thus adding another feature to the diversity of the system’s dynamics.
Yet another interesting dynamics is observed for and . Tormented patterns occur while still generating a periodic behavior, underlining the complexity of the dynamics of the density. On the opposite to this widely varying behavior, the model numerically exhibits convergence to a steady state for the total metastatic burden when initial conditions (for both primary tumor and metastases) are set to . Same apparent convergence also occurs for the number of metastases and the amount of inhibitor (simulations not shown). Looking closer to the volume distribution of metastases at the end of the simulation reveals concentration of the density to the smallest possible volume, suggesting convergence to a Dirac mass located in .
, (logscale)  
Increasing sharp relapses  Bounded and nonperiodic 
, ,  
Complex periodic behavior  Possible convergence to a stable equilibrium 
4. Conclusion
We have presented a mathematical model for systemic inhibition of angiogenesis designed to theoretically study the relative balance of three important biological processes happening in the development of a population of metastases, namely dissemination, local growth and global (systemic) inhibition.
Simulations of this nonlinear model revealed interesting dynamics that underline the complexity. These results illustrate the interest of mathematical and computational models as useful tools for simulating the global result of a complex biology. Indeed, despite the apparent simplicity of our model that reduced here metastatic development to only three essential dynamical processes (birth, growth and inhibition), it would not be possible to guess a priori the result of the combination of the three since, in some situations, this combination even appears as very tormented and unpredictable (see Figure 5) while in others convergence to a steady state is numerically observed.
From a mathematical point of view, the numerical study we performed suggests a wide array of possible asymptotical behavior of the nonlinear model. During the exploration of the parameter space, we observed periodical behaviors with various underlying patterns, from simple oscillations to complex shapes of repeated patterns. In this context, an interesting problem would be to quantify amplitude and frequency of the oscillations in terms of the parameters of the model.
Periodicity was not always the rule since we also observed nonperiodic patterns or convergence to a steady state. These observations suggest possible bifurcations in the parameter space of the infinitedimensional dynamical system. Mathematical study of these dynamical properties, in particular how to go from a stable steadystate to a limit cycle could be interesting perspectives of our work.
From a biological viewpoint, the existence of bounded solutions of the system suggests the possibility of a stable homeostatic burden of metastases that remain in equilibrium due to mutual inhibitory interactions. In our simulations, this happens in situations where growth is substantially altered, either by a reduced growth velocity (parameter ) or smaller initial volume and carrying capacity of a newborn tumor with respect to what was considered here as base values (it should be noted however that within this base parameter set, is the tenth of the maximal reachable volume, which represents an unrealistically large initial tumor volume, making thus smaller values more biologically relevant). This modeling result could shed light on the reported observations (from forensic autopsy studies) of multiple small metastatic foci present in the organism of healthy individuals [37, 26], with a prevalence rate of up to 99% (in the case of thyroid cancer) of the population, that yielded J. Folkman to use the term of “cancer without disease” [17]. However, as shown in [8] from calibration of the model to experimental data, elucidation of this phenomenon from systemic inhibition of angiogenesis only would require a very high production rate (or, equivalently, efficacy parameter ) of the systemic inhibitor. It is more likely that SIA contributes to generate such a situation in combination with other biological players (such as the immune system for instance).
References
 [1] J. AguirreGhiso. Models, mechanisms and clinical evidence for cancer dormancy. Nat Rev Cancer, 7 (11):834â846, 2007.
 [2] D. A August, P. H Sugarbaker, and P. D Schneider. Lymphatic dissemination of hepatic metastases. Implications for the followup and treatment of patients with colorectal cancer. Cancer, 55(7):1490â1494, April 1985.
 [3] D. Barbolosi, A. Benabdallah, F. Hubert, and F. Verga. Mathematical and numerical analysis for a model of growing metastatic tumors. Math Biosci, 218(1):1â14, March 2009.
 [4] S. Benzekry. Mathematical and numerical analysis of a model for antiangiogenic therapy in metastatic cancers. ESAIM, Math. Model. Numer. Anal., 46(2):207â237, 2012.
 [5] S. Benzekry. Mathematical analysis of a twodimensional population model of metastatic growth including angiogenesis. J Evol Equ, 11(1):187â213, December 2011.
 [6] S. Benzekry. Passing to the limit 2D1D in a model for metastatic growth. J Biol Dynam, 6(sup1):19â30, January 2012.
 [7] S. Benzekry, N. André, A. Benabdallah, J. Ciccolini, C. Faivre, F. Hubert and D. Barbolosi. Modelling the impact of anticancer agents on metastatic spreading. Math Model Nat Phenom, 7(1):306â336, 2012.
 [8] S. Benzekry, A. Gandolfi and P. Hahnfeldt. Global Dormancy of Metastases due to Systemic Inhibition of Angiogenesis. submitted, 2013.
 [9] A. Bethge, U. Schumacher, A. Wree and G. Wedemann. Are metastases from metastases clinical relevant? Computer modelling of cancer spread in a case of hepatocellular carcinoma. PloS One, 7(4):e35689, January 2012.
 [10] W. C Black and H. G Welch. Advances in diagnostic imaging and overestimation of disease prevalence and the benefits of therapy. N Engl J Med, 328(17):1237â1243, April 1993.
 [11] P. Chiarella, J. Bruzzo, R. P Meiss, and R. A Ruggiero. Concomitant tumor resistance. Cancer Lett, 324(2):133â141, May 2012.
 [12] A. Devys, T. Goudon and P. Lafitte. A model describing the growth and the size distribution of multiple metastatic tumors. Discrete Contin Dyn SystB, 12(4):731â767, August 2009.
 [13] W. D Dewys. Studies correlating the growth rate of a tumor and its metastases and providing evidence for tumorrelated systemic growthretarding factors. Cancer Res, 32(2):374â379, 1972.
 [14] A. dâOnofrio and A. Gandolfi. Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by Hahnfeldt et al. (1999). Math Biosci, 191(2):159â184, October 2004.
 [15] I. J Fidler and S. Paget. The pathogenesis of cancer metastasis: the âseed and soilâ hypothesis revisited. Nat Rev Cancer, 3(6):453â458, 2003.
 [16] J. Folkman. Angiogenesis inhibitors generated by tumors. Mol Med, 1(2):120â2, January 1995.
 [17] J. Folkman and R. Kalluri. Cancer without disease. Nature, 427(6977):787, 2004.
 [18] J. Folkman. Tumor angiogenesis: therapeutic implications. N Engl J Med, 285:1182â1186, 1971.
 [19] E. Gorelik, S. Segal and M. Feldman. Growth of a local tumor exerts a specific inhibitory effect on progression of lung metastases. Int J Cancer, 21(5):617â625, 1978.
 [20] E. Gorelik. Resistance of tumorbearing mice to a second tumor challenge. Cancer Res, 43(1):138â145, 1983.
 [21] G. P Gupta and J. Massagu Ìe. Cancer metastasis: Building a framework. Cell, 127(4):679â695, November 2006.
 [22] P. Hahnfeldt, D. Panigraphy, J. Folkman and L. Hlatky. Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment, response and postvascular dormancy. Cancer Res, 59(19):4770â4775, October 1999.
 [23] D. Hanahan and R. A Weinberg. Hallmarks of cancer: the next generation. Cell, 144(5):646â674, March 2011.
 [24] L. Holmgren, M. S OâReilly and J. Folkman. Dormancy of micrometastases: balanced proliferation and apoptosis in the presence of angiogenesis suppression. Nat Med, 1(2):149â153, 1995.
 [25] K Iwata, K Kawasaki and N. Shigesada. A dynamical model for the growth and size distribution of multiple metastatic tumors. J Theor Biol, 203(2):177â186, March 2000.
 [26] M. Nielsen, J.L. L Thomsen, S. Primdahl, U. Dyreborg, and J. A Andersen. Breast cancer and atypia among young and middleaged women: a study of 110 medicolegal autopsies. Br J Cancer, 56(6):814â819, 1987.
 [27] M. S OâReilly, L. Holmgren, Y. Shing, C. Chen, R. A Rosenthal, M. Moses, W. S Lane, Y. Cao, E. H Sage and J. Folkman. Angiostatin: a novel angiogenesis inhibitor that mediates the suppression of metastases by a Lewis lung carcinoma. Cell, 79(2):315â28, October 1994.
 [28] M. S OâReilly, T. Boehm, Y. Shing, N. Fukai, G. Vasios, W. S Lane, E. Flynn, J. R Birkhead, B. R Olsen and J. Folkman. Endostatin: an endogenous inhibitor of angiogenesis and tumor growth. Cell, 88(2):277â285, 1997.
 [29] B. Perthame. Transport equations in biology. Birkhauser Verlag, 2007.
 [30] R. T Prehn. Two competing influences that may explain concomitant tumor resistance. Cancer Res, 53(14):3266â9, July 1993.
 [31] M. Retsky, R. Demicheli, W. Hrushesky, M. Baum and I. Gukas. Surgery triggers outgrowth of latent distant disease in breast cancer: an inconvenient truth? Cancers, 2(2):305â337, March 2010.
 [32] E. Rofstad and B. Graff. Thrombospondin1mediated metastasis suppression by the primary tumor in human melanoma xenografts. J Invest Dermatol, 117(5):1042â9, November 2001.
 [33] A. Sckell, N. Safabakhsh, M. Dellian and R. K Jain. Primary tumor sizedependent inhibition of angiogenesis at a secondary site: an intravital microscopic study in mice. Cancer Res, 58(24):5866â5869, 1998.
 [34] E. V Sugarbaker, A. M Cohen and A. S Ketcham. Do metastases metastasize? Annals of surgery, 174(2):161â6, August 1971.
 [35] C. R Tait, D. Dodwell and K Horgan. Do metastases metastasize? J Pathol, 203(1):515â8, May 2004.
 [36] O. V Volpert, J. Lawler and N. P Bouck. A human fibrosarcoma inhibits systemic angiogenesis and the growth of experimental metastases via thrombospondin1. Proc Natl Acad Sci USA, 95(11):6343â8, May 1998.
 [37] H. G Welch and W. C Black. Overdiagnosis in cancer. J Natl Cancer Inst, 102(9):605â13, May 2010.