Vlastimil Křivan

On the Gause predator-prey model with a refuge: A fresh look at the history

Keywords: discontinuous differential equations, Filippov solution, Holling functional response, limit cycle, Lotka-Volterra model, population dynamics

This article re-analyses a prey-predator model with a refuge introduced by one of the founders of population ecology G. F. Gause and his co-workers to explain discrepancies between their observations and predictions of the Lotka-Volterra prey-predator model. They replaced the linear functional response used by Lotka and Volterra by a saturating functional response with a discontinuity at a critical prey density. At concentrations below this critical density prey were effectively in a refuge while at a higher densities they were available to predators. Thus, their functional response was of the Holling type III. They analyzed this model and predicted existence of a limit cycle in predator-prey dynamics. In this article I show that their model is ill posed, because trajectories are not well defined. Using the Filippov method, I define and analyze solutions of the Gause model. I show that depending on parameter values, there are three possibilities: (1) trajectories converge to a limit cycle, as predicted by F. G. Gause, (2) trajectories converge to an equilibrium, or (3) the prey population escapes predator control and grows to infinity.

three experimental predator-prey systems consisting of prey (Aleuroglyphus agilis) and predatory (Cheyletus eruditus) mites, prey (Paramecium caudatum) and predatory (Didinium nasatum) protists, and yeast (Saccharomyces exiguus) and protists (Paramecium bursaria). In none of these experiments population dynamics were consistent with the Lotka-Volterra neutrally stable limit cycles. In experiments with protists and mites a prolonged coexistence of both prey and predators was obtained only when both species were regularly added to the system. In a completely homogeneous environment Didinium destructed all prey and it collapsed subsequently.

When the environment was not homogeneous and there was a refuge for prey, prey survived in the refuge but the predator population collapsed. The situation was different when protists fed on the yeast. There was strong experimental evidence that population dynamics tended to a limit cycle that was independent from initial population numbers. These observations lead [START_REF] Gause | Further studies of interaction between predators and prey[END_REF] to search for discrepancies in assumptions of the Lotka-Volterra predator-prey model when applied to their experimental systems. First, they observed that protists were not able to feed on yeast at low densities, because at low yeast densities the prey formed into a sediment at the bottom which was not accessible to predators inhabiting the water column. Thus, prey were effectively in a refuge when at low concentrations. When prey reached above the critical density, they re-appeared in the water column and became accessible to predators. Second, they postulated that consumption of prey was a saturating function of prey density. Third, they observed that the predator population did not start to decrease until most of the predators were in the refuge.

To describe their observations mathematically, [START_REF] Gause | Further studies of interaction between predators and prey[END_REF] substituted the linear consumption rate used in the Lotka-Volterra model by a saturating function that was zero below the critical prey density and had a jump (discontinuity) at the critical prey density. Hereafter I will refer to this model as the Gause model. Models with a saturating functional response were introduced to theoretical ecological literature later on by e. g., [START_REF] Rosenzweig | Graphical representation and stability conditions of predator-prey interactions[END_REF] who replaced the linear functional response in the Lotka-Volterra model by the Holling type II [START_REF] Holling | Some characteristics of simple types of predation and parasitism[END_REF] functional response. Besides the fact that the Rosenzweig-MacAthur model assumes negative density dependent prey growth, another difference between these two models is the jump in the functional response in the Gause model. The Rosenzweig-MacArthur model became one of the key models of prey-predator interactions because it documents that prey-predator coexistence is not limited to an equilibrium. So did the Gause model. However, the mechanism that leads to the limit cycle in the Rosenzweig-MacArthur model is entirely different from the mechanism that causes fluctuations in the Gause model. In the Rosenzweig-MacArthur model the limit cycle is caused by interaction between bottom-up and top-down prey regulation. When prey are limited by resources (i. e., when environmental carrying capacity is low), bottom-up regulation is strong and keeps predator and prey population dynamics at an equilibrium. As environmental carrying capacity increases due to enrichment, bottom-up regulation gets weaker and the population equilibrium becomes unstable due to the destabilizing Holling type II functional response. In the case of the Gause model there is no bottom-up regulation of prey growth so the nature of the limit cycle is solely due to top-down regulation and refuge presence. In fact, the Gause functional response is of the Holling type III, so their analysis clearly showed that such a functional response can lead to a limit cycle in predator-prey population dynamics.

In general, proving existence and uniqueness of a limit cycle in predator-prey population dynamics is not trivial. For example, the limit cycle in the Rosenzweig-McArthur model appears when the stable equilibrium undergoes the Hopf bifurcation and its uniqueness was proved only much later (e. g., [START_REF] Huang | Conditions for uniqueness of limit cycles in general predator-prey systems[END_REF].

Certainly, these concepts were unknown when [START_REF] Gause | Further studies of interaction between predators and prey[END_REF] analyzed their system. Therefore, it is quite remarkable that using a simple geometrical argument, they were able to predict that trajectories of their model converge to a limit cycle (see their Figure 5). However, there is one crucial problem with their analysis. The Gause model was described by a differential equation with a discontinuous right-hand side (due to the jump in the functional response). Such differential equations may not have solutions in the usual sense 1 and the Gause model is an example. A concept of a "solution" for such models was introduced later on by A. F. [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF] (see also [START_REF] Aubin | Differential inclusions[END_REF][START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF]. No such mathematical concept existed at the time when Gause with his co-workers analysed their model to achieve a better fit with experimental data. This is really a remarkable example of biological research that used a mathematical methodology that was not yet developed at that time.

In this article I will explain how solutions of the Gause model can be defined and will analyze this model with respect to parameter values and initial population densities. I will show that trajectories of the Gause model can converge to a limit cycle as correctly predicted by [START_REF] Gause | Further studies of interaction between predators and prey[END_REF], but they can also converge to an equilibrium, or prey can escape completely predator regulation. (

The Gause model

) 1
Here R is prey density, C is predator density, r is the per capita prey population growth rate, f is the Gause functional response specified below, e is the efficiency rate with which captured prey are converted to new predators, and m is the predator mortality rate. There are three distinctive features of this functional response: 1. below a critical prey population threshold (denoted by R c ) prey are not consumed, 2. the functional response has a discontinuity at the threshold, 3. consumption saturates with increasing prey density [START_REF] Gause | Further studies of interaction between predators and prey[END_REF]. The first feature corresponds to the refuge of a fixed size R c . The jump at the critical prey density was motivated Unless we set dx dt = 0 when x = 0, this differential equation has no solution starting at x(0) = 0.

by their observation that "... predators ... do not seriously decrease in concentration until the destruction of the prey down to this threshold..." [START_REF] Gause | Further studies of interaction between predators and prey[END_REF]. This suggest that the functional response in the vicinity of R c is quite steep, and can be approximated by a functional response with a jump at R c . In fact, as we will see below, this assumption allowed authors to analyze the predator-prey population dynamics using a simple geometric argument. Moreover these authors also asked that "the tangent to f ... crosses the ordinate and not the abscissa". This latter condition excludes coexistence of prey and predators in a locally asymptotically stable equilibrium and it is interesting to note that this seems to be the first occurrence of the now well known (in)stability condition for prey-predator interactions [START_REF] Murdoch | Predation and population stability[END_REF]. A prototype of such a functional response is (see the solid curve in Figure 1A)

f (R) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 0 R < R c λR 1 + hλR R ≥ R c .
(2)

Here parameter λ describes the search rate of a predator and h is the handling time a predator needs to process one unit of prey. R c is the critical prey density below which prey are not accessible to predators. Thus, for above-critical prey density, f is the Holling type II functional response. The Gause functional response (2) is a limiting case of the Holling type III functional response (dashed curve in Figure 1A)

f III (R) = λR μ /(R μ + R μ s )R 1 + hλR μ /(R μ + R μ s )R (3)
when the exponent μ tends to infinity. I will show now that model ( 1) is not well defined in the sense that solutions cannot be continued in forward time once prey density reaches the critical threshold R c , and predator density is high enough.

For prey population densities below the threshold (R < R c ), prey not being eaten by predators grow exponentially while predators die exponentially, i. (4) Thus, at each point of the prey-predator density phase space to the left of the critical prey density R = R c , trajectories move in the south-east direction.

For prey population densities above the threshold (R > R c ), population dynamics (1) are given by the Lotka-Volterra model with the Holling type II functional

response dR dt = rR -C λR 1 + hλR dC dt = eλR 1 + hλR -m C.
(5)

First, I will recall properties of model ( 5) when no refuge exists (R c = 0, e. g., [START_REF] Murdoch | Predation and population stability[END_REF][START_REF] Svirezhev | Stability of biological communites[END_REF][START_REF] Křivan | Prey-predator models[END_REF]. For small handling times that satisfy h < e/m this model has an interior, unstable equilibrium

E = {R , C } = { m λ(e-mh)
, re λ(e-mh) } (Appendix A). Depending on parameters and initial conditions, model ( 5) has two types of solutions (Appendix A). For small handling times that satisfy h < e/(r + m) all trajectories of model ( 5) spiral outward from the interior equilibrium with ever increasing amplitude. For larger handling times (e/m > h > e/(r + m)) there exist trajectories that tend to infinity without spiraling around the equilibrium (see the bottom trajectory in Figure 2B,F). In such a case, the prey population completely escapes predator regulation.

Second, I will assume a refuge exists (R c > 0) and I will consider model (1) with the Gause functional response (2). For R > R c the prey isocline coincides with the prey isocline of model ( 5) and it is given by

C = r λ (1 + λhR)
(thick dashed line in Figure 2). For small prey densities satisfying R < R c the prey isocline is not defined. Let us consider the line R = R c in the prey-predator density phase space. The isocline analysis shows that in the vicinity and to the right of this line, the prey population decreases provided predator density is high enough, i. e., 5). This new (multivalued) vector field is given by

C > C c = r λ (1 + λhR c ). ( 6
F = αf 1 + (1 -α)f 2
where 0 ≤ α ≤ 1. In other words, this vector field associates to every point along the vertical part of prey isocline a whole set of possible directions given by F . This definition of the Filippov field re-defines the functional response (2) at the critical prey density to

f (R) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 R < R c [0, λRc 1+λhRc ] R = R c λR 1 + hλR R > R c . ( 7
)
It follows that under this new definition f (R c ) is the line segment [0, λRc 1+λhRc ] that fills the gap in the Gause functional response (Figure 1B). This definition is very natural as it reflects the fact that at the critical prey density the functional response does not specify exactly the prey consumption by predators. Thus, functional response ( 7) is a very reasonable approximation of a continuous Holling type III functional response that is steep enough at prey densities close to the critical prey density (Figure 1B).

Appendix B shows that model (1) with Filippov vector field at the critical prey density has uniquely defined trajectories for every initial condition.

To analyze the Gause model we need to know the dynamics along the vertical part of the prey isocline. When trajectory of model ( 1 is horizontal there). It is interesting to observe that there exists a new equilibrium 8), provided this point is on the vertical part of the prey isocline, i. e., when C F > C c . This holds whenever the prey critical density satisfies R c > m/(λ(emh)), i. e., when the predator isocline of model ( 5) is to the left of the prey critical line R = R c , or, in other words, if the prey refuge is large enough. Now I will analyze qualitative behavior of the solutions of the Gause model with respect to parameters and initial population densities. I will consider two cases, depending whether the critical prey population density R c is smaller or larger than is the equilibrium prey density R = m/(λ(emh)) of model ( 5). First, I will assume that the critical prey density is smaller (R c < R , Figure 2A-D). Figure 2A shows the case where the handling time is so small that all trajectories of model ( 5) spiral around the equilibrium. In this case, starting from any initial population densities, the corresponding trajectory must fall at certain time onto the vertical part of the prey isocline. The trajectory will then follow dynamics described by ( 8) and because in this case there is no equilibrium on the vertical part of the prey isocline, the trajectory must reach at a positive time the lower end-point (R c , C c ) (i. e., the corner of the prey isocline). At this point it will leave the isocline and a limit cycle is formed (the heavy solid cycle in Figure 2A). All trajectories will reach this limit cycle in a finite time. This is the limit cycle predicted by [START_REF] Gause | Further studies of interaction between predators and prey[END_REF].

E F = (R F , C F ) = (R c , erR c /m) of model (
For larger handling times (i. e., e/m > h > e/(r + m)), model ( 5) also has trajectories that do not spiral around the equilibrium. There are two possibilities depending on the position of the critical point (R c , C c ). Either the trajectory that starts at this point returns to the critical line R = R c (Figure 2B) in which case a limit cycle exists. However, this limit cycle is not globally stable, because there are trajectories that tend to infinity (see the bottom trajectory in Figure 2B). Numerical simulations show that for yet higher handling times the trajectory that starts at the point (R c , C c ) tends to infinity (Figure 2C) and predators lose control over their prey.

Second, I consider the case where R c > R (Figure 2E,F). In this case, equilibrium E of model ( 5) is located to the left of the prey critical density line in the prey-predator phase space where dynamics are described by ( 4). Thus, it is no longer an equilibrium of the Gause model. However, the equilibrium E F is on the prey isocline and therefore is the only equilibrium of the Gause model. Similar to the case where the prey critical density line was to the left of the predator isocline, for small handling times all trajectories must fall onto the vertical part of prey isocline and then move along it to reach the equilibrium E F (Figure 2E). For higher handling times, some trajectories will escape predator regulation (Figure 2F) so that the equilibrium is then only locally asymptotically stable.

These results are summarized in Figure 3. This is an example of a bifurcation diagram in the (h, R c ) parameter space. The solid curve in Figure 3 corresponds to the transition from the locally stable limit cycle in Figure 2A to the locally stable equilibrium (panel E). This curve is given by equation R c = R * , i. e., R c = m/(λ(emh)). Such bifurcation was called the boundary-focus bifurcation (see the case BF 3 in Figure 5 in [START_REF] Kuznetsov | One-parameter bifurcations in planar Filippov systems[END_REF]. The dashed line in Figure 3 corresponds to transition between Figure 2B and Figure 2C, i. e., disappearance of the limit cycle in Figure 2B, and was obtained by numerical simulations.

Discussion

In this article I have analyzed a prey-predator model introduced by [START_REF] Gause | Further studies of interaction between predators and prey[END_REF] to explain cycles in prey-predator population dynamics which they observed in some experiments. They replaced the linear functional response used by Lotka and

Volterra by a functional response that was zero below some critical prey threshold density, had a "jump" at this threshold, and was saturating at high prey densities.

Such a functional response is an extreme form of the Holling type III sigmoid functional response that was introduced later on [START_REF] Holling | Some characteristics of simple types of predation and parasitism[END_REF]. Using this functional response [START_REF] Gause | Further studies of interaction between predators and prey[END_REF] were able to show that predators and prey can coexist along a limit cycle. This is perhaps the very first models in ecology that shows species coexistence is possible not only at an equilibrium. They also showed that a prey refuge can lead to predator-prey coexistence as the same model without a refuge does not predict such a coexistence. Moreover, they derived a general stability condition for predator-prey models. Thus, this article was really fundamental for further development of population ecology. Although their analysis is very elegant, there is one crucial problem because their model is not well defined as there are initial population densities for which no solutions exist. The question then is if and in which sense the predictions these authors obtained from the model are correct or not. In this article, using the approach developed by [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF], I showed how trajectories in the Gause model can be defined. I also analyzed this model with respect to parameters and initial population densities. [START_REF] Gause | Further studies of interaction between predators and prey[END_REF] predicted that the model has a globally stable limit cycle (see their Figure 5). I showed that predators and prey can coexist along a limit cycle provided that the refuge size and handling times are not too high (see the region LC in Figure 3). In fact, for small handling times this limit cycle is globally stable, but for slightly higher handling times there are also trajectories along which both populations tend to infinity (Fig- ure 2B). Along such trajectories prey "escape"predator control due to low predator efficiency to handle prey items. For yet higher handling times the limit cycle disap-pears and predators cannot control prey growth (see the region in Figure 3 denoted as ∞ ).

A different situation occurs when prey critical density is high enough, i. e., when the prey refuge is large. Then, for small handling times a globally asymptotically stable equilibrium point exists. As handling time increases, there will be again some trajectories that will tend to infinity.

As prey below the critical density are effectively in a refuge, this analysis nicely demonstrates the stabilizing role of refugia that protect a fixed number of prey (e. g., [START_REF] Rosenzweig | Graphical representation and stability conditions of predator-prey interactions[END_REF][START_REF] Hassell | Stability in insect host-parasite models[END_REF][START_REF] Smith | Models in ecology[END_REF][START_REF] Murdoch | Predation and population stability[END_REF][START_REF] Hassell | The dynamics of arthropod predator-prey systems[END_REF][START_REF] Sih | Prey refuges and predator-prey stability[END_REF][START_REF] Ives | Antipredator behaviour and the population dynamics of simple predator-prey systems[END_REF][START_REF] Ruxton | Short term refuge use and stability of predator-prey models[END_REF][START_REF] Hochberg | Refuge evolution and the population dynamics of coupled host-parasitoid associations[END_REF][START_REF] Křivan | Effects of optimal antipredator behavior of prey on predator-prey dynamics: role of refuges[END_REF]. When the refuge is large enough, i. e., when the critical prey density below which prey are protected is high, provided predators are able to regulate the prey population (i. e., when handling times are small) then the two species coexist at an equilibrium. As refuge size decreases, less prey are protected and prey-predator coexistence occurs along a limit cycle.

Because the Gause functional response is of the Holling type III, analysis of the Gause model shows that such a functional response can lead to predator and prey coexistence without any prey density dependence. In fact, the effect of a refuge on predator-prey stability is somewhat parallel to the effect of enrichment on predatorprey stability in the Rosenzweig-MacArthur model. While in the Rosenzweig-MacArthur model a higher environmental capacity weakens the bottom-up control and leads to population oscillations, the effect of refuge is just opposite. A larger refuge stabilizes otherwise unstable population equilibrium.

A similar analysis for the Lotka-Volterra model where functional response is linear with a "jump" at the critical prey density shows the same pattern (Figure 2D).

Here all trajectories converge to a global attractor which is bounded by the largest Lotka-Volterra cycle that is to the right of the vertical part of the prey isocline. For larger R c the population dynamics converge on an equilibrium (results not shown) exactly as in the case with the Holling type II functional response.

The ultimate reason for species coexistence at an equilibrium or along a limit cycle in the Gause model is the fact that the prey isocline has a vertical segment.

This vertical segment then limits the amplitude of species fluctuations as suggested by [START_REF] Rosenzweig | Graphical representation and stability conditions of predator-prey interactions[END_REF] in their seminal work. Since then it was shown that such isoclines arise when prey-predator models consider optimal foraging of predators [START_REF] Křivan | Optimal foraging and predator-prey dynamics[END_REF][START_REF] Van Baalen | Alternative food, switching predators, and the persistence of predator-prey systems[END_REF], optimal activity level of prey and/or predators [START_REF] Křivan | The Lotka-Volterra predator-prey model with foraging-predation risk trade-offs[END_REF] or optimal use of refuges [START_REF] Křivan | Effects of optimal antipredator behavior of prey on predator-prey dynamics: role of refuges[END_REF]. All these models show that prey isoclines with vertical segments (or equivalently predator isoclines with horizontal segment) can arise naturally in real populations. In fact, using models with discontinuous right-hand sides is a very natural methodology for analyzing such systems and can lead to a much simpler analysis when compared to models with strong non-linearities. Such an approach was used for example by [START_REF] Crowley | Dispersal and the stability of predator-prey interactions[END_REF] to demonstrate an existence of a limit cycle in [START_REF] Rosenzweig | Graphical representation and stability conditions of predator-prey interactions[END_REF] model where prey isocline was assumed to be vertical at the prey carrying capacity. This corresponds to strong density dependence that operates only when prey density reaches the carrying capacity.

It is likely that the Gause predator-prey model was one of the first models in biology where discontinuous differential equations were used. It is remarkable that this happened even before mathematicians provided a formal definition of a solution for such systems. Although it does not seem that Gause work motivated some research in this field, this example clearly shows that biology can be a source of interesting problems for mathematicians. Since that time, similar models were used in ecology (e. g., [START_REF] Křivan | Optimal foraging and predator-prey dynamics[END_REF][START_REF] Meza | Threshold policies control for predator-prey systems using a control Liapunov function approach[END_REF][START_REF] Dercole | Bifurcation analysis of piecewise smooth ecological models[END_REF][START_REF] Křivan | The Lotka-Volterra predator-prey model with foraging-predation risk trade-offs[END_REF], and in gene networks (e. g., [START_REF] Edwards | Analysis of continuous-time switching networks[END_REF][START_REF] De Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF][START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF][START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF]. Moreover, a complete qualitative theory was developed for two dimensional models with a discontinuity [START_REF] Kuznetsov | One-parameter bifurcations in planar Filippov systems[END_REF], and a systematic approach (called "the puzzle method") for analyzing such models was developed [START_REF] Dercole | Bifurcation analysis of piecewise smooth ecological models[END_REF]. Moreover, software for numerical analysis of such models is also available (Dercole andKuznetsov, 2004, 2005;[START_REF] Piiroinen | An event-driven method to simulate Filippov systems with accurate computing of sliding motions[END_REF][START_REF] Thota | Tc-hat ( T C): A novel toolbox for the continuation of periodic trajectories in hybrid dynamical systems[END_REF]. These methodological advances open further possibilities to study models with discontinuities that were not available at times when [START_REF] Gause | Further studies of interaction between predators and prey[END_REF] wrote their pioneering treatise. and the refuge size. The solid curve corresponds to the transition from the limit cycle (Figure 2A) to the stable equilibrium (Figure 2E) and is given by R c = m λ(e-mh) . The dashed curve corresponds to the set of parameters where the limit cycle disappears.

LC stands for the limit cycle, Eq stands for the interior equilibrium, ∞ means that no bounded attractor exists, and the right region (Cons. ext.) denotes the set of parameters for which consumers extinct and prey grow exponentially.

  [START_REF] Gause | Further studies of interaction between predators and prey[END_REF] considered the following adaptation of the Lotka-Volterra predator-prey model dR dt = rR -Cf (R) dC dt = (ef (R)m)C.

1

  A trivial example is the differential equation dx dt = 1 when x < 0 and dx dt = -1 when x > 0.

  )Here the point (R c , C c ) corresponds to the intersection of the prey isocline of model (5) with the critical prey line R = R c (the corner of the thick dashed line in Figure2). Similarly, to the left of this line, the prey population increases. It is clear that above the critical predator density C c trajectories of the Gause model are pushed from both sides to the line R = R c and they cannot leave this line if they exist. Thus, trajectories of model (1) cannot cross the critical line R = R c above the point C c and, by definition (the prey isocline consists of points in the prey-predator density phase space where prey population neither increases nor decreases,[START_REF] Rosenzweig | Graphical representation and stability conditions of predator-prey interactions[END_REF], this half-line must be a part of the prey isocline of the Gause model. Thus, the prey isocline of the Gause model is L-shaped (the thick dashed line in Figure2). However, the most interesting feature of the Gause model is the fact that once a trajectory falls on the vertical part of the prey isocline, it cannot be continued any further. Indeed, as the "trajectory" cannot leave the line R = R c above the point C c it follows that dR dt = 0. But the right-hand side of model (1) evaluated at R = R c gives dR c dt = 0. In other words, the Gause model is not well posed because its trajectories are not defined when they fall on the vertical part of the prey isocline. This is a consequence of the fact that the Gause functional response has a "jump" at the critical prey density, because such models may not have solutions.In other words it is not clear how[START_REF] Gause | Further studies of interaction between predators and prey[END_REF] defined the vector field in their Figure 5, panel 1' for R = R c . It took more than twenty years after Gause et al. (1936) published their article before Filippov (1960) (see also Filippov, 1988) introduced a new solution concept for such models. The crucial step is provided through suitable definition of the vector field at the critical prey density, which I briefly describe now. The Filippov solution concept applied to the Gause model defines a new vector field at the critical prey density R c as the line segment with end points given by the two adjacent vector fields f 1 and f 2 . Here f 1 = (rR c , -mC) stands for the vector field defined by the right-hand side of model (4) and f 2 = (rR c -CλR c /(1 + λhR c ), eλR c C/(1 + λhR c ) -mC) for the vector field defined by model (

  ) falls on the vertical part of the prey isocline and the predator density satisfies C > C c , the trajectory cannot leave the isocline and it must move vertically, i. e., dR/dt = 0. This implies that along such a trajectory α = 1 -r(1 + λhR c ) λC , which corresponds to the vector f F = (0, erR c -mC) ∈ F, describes population dynamics along the line R = R c as long as C > C c . At the point (R c , C c ) the trajectory leaves the line R = R c and because the consumer density is low, it enters the region where R > R c (this also follows from considerations in Appendix B). Thus, model (8) provides a formal definition of the vector field along the horizontal line in Figure 5, panel 1' of Gause et al. (1936) (I remark that axes are swapped in Gause et al. (1936) so the line R = R c

Figure legends : Figure 1 :

 :1 Figure legends: Figure 1: Panel A (solid line) shows the Gause functional response (2). Panel B shows the Filippov regularization (7) (solid line) of the Gause functional response. The dashed line in both panels is the Holling type III functional response given by (3). Parameters: R c = 1, λ = 10, h = 0.1, μ = 10.

Figure 2 :

 2 Figure 2: Panels A-D assume that the critical prey density (R c = 0.2) is smaller than the equilibrium prey density R of model (5), while panels E and F assume the opposite case (R c = 0.6). Panel A shows the trajectories of the Gause model when handling time is small (h = 0.1) and all trajectories converge to a limit cycle.For intermediate handling times (h = 0.5, panel B) a locally asymptotically stable limit cycle exists, but there are also trajectories that tend to infinity. For yet larger handling times (h = 0.7, panel C) all trajectories diverge from the equilibrium. Panel D shows the case where handling time is zero. Panel E assumes low handling time (h = 0.1) in which case all trajectories converge to an equilibrium point that is located at the vertical part of the prey isocline. Panel F assumes a larger handling time (h = 0.7) where the equilibrium is only locally asymptotically stable because some trajectories tend to infinity. Other parameters: r = 1, λ = 1, m = 0.2, e = 0.5.

Figure 3 :

 3 Figure 3: A bifurcation diagram for the Gause model with respect to handling times

  Figure 1:

Acknowledgments I thank people at the NIMBioS who made my sabbatical fellowship there very pleasant and fruitful. A part of this work was conducted while a Sabbatical Fellow at the National Institute for Mathematical and Biological Synthesis, an Institute sponsored by the National Science Foundation, the U.S. Department of Homeland Security, and the U.S. Department of Agriculture through NSF Award #EF-0832858, with additional support from The University of Tennessee, Knoxville. The Institute of Entomology is funded by the Academy of Sciences of the Czech Republic (Z50070508).

Appendix A: Analysis of model ( 5) At the interior equilibrium E = {R , C } = { m λ(e-mh) , re λ(e-mh) } of model ( 5), the Gause criterion

the [START_REF] Dulac | Recherche des cycles limites[END_REF] criterion (see also [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF] with the Dulac function 1 RC excludes existence of a limit cycle for model (5). When prey density is high, model ( 5) asymptotically tends to

The prey isocline is the line C = hrR. Below this isocline prey population increases while above it decreases. Let us consider a line given by C = kR with 0 < k < hr (i. e., this line is below the prey isocline in the prey-predator density phase space).

Trajectories of model (9) cross this line upwards provided k > h(m + r)e because under this condition dC/dR > k at every point of this line. It follows that when

h < e/(m + r) then k > h(m + r)e for every 0 < k < hr and trajectories of model (9) will reach the prey isocline, i. e., trajectories of model (5) will spiral around its equilibrium. However, when h > e/(m + r) there will be trajectories of model ( 9) that never reach the prey isocline (see bottom trajectories in Figure 2B,F). Along such trajectories prey population escapes predator control.

Appendix B: Existence and uniqueness of trajectories of the Gause model

Let n = (n R , n C ) = (1, 0) be the vector perpendicular to the line R = R c in the prey-predator density phase space. Projection of the two vector fields given by the right-hand sides of (4) (denoted as f 1 ) and ( 5) (denoted as f 2 ) are n, f 1 = rR c and n, f 2 = rR c -λRcC 1+λhRc . If n, f 1 > 0 and n, f 2 < 0 trajectories are pushed from both below and above to the line R = R c . These conditions hold when predator population density is large enough, i. e., C > C c .

The existence of trajectories for the Gause model follows from general existence theorems that can be found in [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF] (see also [START_REF] Colombo | Selective strategies in food webs[END_REF].

Uniqueness of trajectories for the Gause model follows from the fact that n, f 1 = n, f 2 + λRcC 1+λhRc . Thus, it is impossible that at some points of discontinuity n, f 1 < 0 and n, f 2 > 0 which would imply non-uniqueness of trajectories [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF][START_REF] Colombo | Selective strategies in food webs[END_REF].


J W Schofield
email: james.schofield@st-hughs.ox.ac.uk

E A Gaffney

R A Gatenby

P K Maini

Tumour Angiogenesis: The Gap between Theory and Experiments

Keywords: Avascular, Spheroid, Modelling, Cornea

A common experimental technique for viewing in vivo angiogenesis utilises tumours implanted into a test animal cornea. The cornea is avascular but the tumour promotes vascularisation from the limbus and the new blood vessels can be readily observed through the transparent cornea. Many of the early mathematical models for tumour angiogenesis used this scenario as their experimental template and as such assumed that there is a large gap, of the order of 2 mm, between the tumour and neighbouring vasculature at the onset of angiogenesis. In this work we consider whether the assumption that there is a significant gap between the tumour and neighbouring vasculature is unique to intra-cornea tumour implants, or whether this characterises avascular tumour growth more generally. To do this we utilise a simple scaling argument, derive a multi-compartment model for tumour growth, and consider in vivo images. This analysis demonstrates that the corneal implant experiments and the corresponding mathematical models cannot generally be applied to a clinical setting.

Introduction 1

Folkman [START_REF] Folkman | The vascularization of tumours[END_REF] proposed that tumours un-2 dergo three distinct phases of growth: an initial avascu-3 lar phase of slow limited growth, then an angiogenesis 4 phase where the tumour cultivates its own blood supply, 5 and finally a phase of rapid vascularised growth. Folk-6 man also proposed that inhibiting angiogenesis could 7 provide an effective anti-tumour therapy. 8

A common technique for observing in vivo angiogen-9 models of angiogenesis have been developed to recreate the results of intra-corneal tumour implant experiments by assuming there is a significant gap, of the magnitude of the width of the cornea, between the tumour and the neighbouring vasculature [START_REF] Balding | A mathematical model of tumour-induced capillary growth[END_REF][START_REF] Chaplain | A model mechanism for the chemotactic response of endothelial cells to tumour angiogenesis factor[END_REF][START_REF] Byrne | Mathematical models for tumour angiogenesis: Numerical simulations and nonlinear wave solutions[END_REF][START_REF] Chaplain | Mathematical modelling of angiogenesis[END_REF][START_REF] Plank | Lattice and non-lattice models of angiogenesis[END_REF]). These investigations have been very successful at modelling the cornea implant experiments and have been well verified in context; however it is not obvious that a significant gap should exist between a tumour and the neighbouring capillaries when an in vivo tumour grows in "normal" vascular tissue. Furthermore, while many modelling explorations have previously invoked analogous assumptions when considering normally vascularised tissue, recent imaging suggests that there is no significant gap between an in vivo avascular tumour and its neighbouring capillaries [START_REF] Shubik | Vascularization of tumors: A review[END_REF], as illustrated by Figure 5(a).

Thus there is a need to explore this potential discrepancy between a modelling framework for vascular tumour angiogenesis and recent observations, and we firstly explore this by utilising a simple scaling argument. For the intra-cornea implant experiments to directly apply to a typical in vivo tumour the in vivo tumour must develop and/or maintain a significant gap between itself and the neighbouring capillaries. Consider a spherical tumour of radius r tum (t), a function of time t, growing in vascular tissue. Then the volume of the tumour, V tum , is

V tum = 4πr 3 tum 3 .
Now suppose that following the initiation of cancer in a single cell there is a region of non-cancerous cells between the tumour and the surrounding vascu-lature which remains of constant width as the tumour 57 grows. To be consistent with the intercapillary dis-58 tance in vascular tissue the distance between the tumour 59 and neighbouring vasculature should be approximately 60 50 μm [START_REF] Fait | Microvascular patterns of the human large intestine: Morphometric studies of vascular parameters in corrosion casts[END_REF] and the volume of tissue in such 61 a gap between the tumour and surrounding vasculature 62

(in μm 3 ), V gap , is then 63

V gap = 4π (r tum + 50) 3 -r 3 tum 3 .
It can then be shown that the maintenance of this gap 64 of constant width requires the non-cancerous tissue to 65

increase in volume at the relative rate compared to can-66 cer tissue of 67

1 V gap dV gap dt 1 V tum dV tum dt = r tum (r tum + 50) 2 -r 3 tum (r tum + 50) 3 -r 3 tum .
To use this result to compare the cellular proliferation 68 rates of cancerous and non-cancerous cells we recall 69 that tumours of radius three cell widths, which corre-70 sponds to a r tum of approximately 80 μm (Casciari et al., 71 1992;[START_REF] Buckley | MRI measurement of cell volume fraction in the perfused rat hippocampal slice[END_REF] will typically consist of only 72 proliferating cells [START_REF] Larue | Microenvironmental regulation of proliferation in multicellular spheroids is mediated through differential expression of cyclin-dependent[END_REF]. For this value of 73 r tum , the tumour radius, we find that the non-cancerous 74 cells must initially proliferate approximately 50% as 75 quickly as the cancerous cells. Having non-cancerous 76 cells proliferate half as quickly as cancerous cells may 77 be faster proliferation than would be anticipated, but is 78 not completely unfeasible for particularly slow grow-79 ing tumours or tumours growing in tissue with a par-80 ticularly high cell turnover rate. However in general no 81 obvious mechanism exists which would upregulate non-82 cancerous cell proliferation to this extent. This suggests 83 that any initial gap between the tumour and the neigh-84 bouring vasculature will begin to close as the tumour 85 growth rate of the tumour.

For the mathematical models based upon the intracornea experiments to be more widely applicable the assumption that there is a significant gap between the tu- We now consider how the distance between a typical in vivo avascular tumour and the neighbouring capillaries evolves as the tumour grows. In particular we consider whether a significant gap can be maintained between the tumour and capillaries, and if not then why does the gap close and how quickly does it do so. The results of this study will then imply whether or not the mathematical models based upon the intra-cornea implant experiments can be applied in a general context.

The Model

To consider how the distance between an in vivo tumour and the neighbouring capillaries evolves as the tumour grows we model the growth of an avascular tumour beginning with a small collection of cells and finishing with an avascular tumour which has reached diffusion limited saturation. We consider populations of healthy non-cancerous cells, proliferating cancerous cells, quiescent cancerous cells and necrotic cells. We assume that the system can be modelled as a continuum of cells and exhibits radial symmetry. These assumptions form an idealised model for avascular tumour growth, but one which constitutes a best case scenario for maintaining a significant gap between the tumour 119 and neighbouring vasculature. If a significant gap be-120 tween the tumour and neighbouring vasculature cannot 121 be maintained by our idealised model, then we would 122 not expect it to be maintained with the introduction of 123 asymmetries which effectively reduce the gap in at least 124 one direction. 125

The domain of our model consists of the tumour 126 spheroid and a region of the non-cancerous cells out-127 side the tumour up to the nearest capillaries, as shown 128

in Figure 1. 129 rate of oxygen consumption by cells [START_REF] Lin | Oxygen diffusion in a spheroid cell with non-linear uptake kinetics[END_REF]. 139

We assume that proliferating cells can become qui-140 escent and vice versa. Whilst the local oxygen con-141 centration can be a factor in determining this change of 142 state, other factors such as cell-cell interaction and the 143 availability of growth factors are also likely to be influ-144 ential [START_REF] Mueller-Klieser | Tumour biology and experimental therapeutics[END_REF]. The effect of all these 145 potential factors is too complex to be considered in de-146 tail. Instead we note that proliferating cells are generally 147 found on the outer rim of an avascular tumour and qui-148 escent cells are generally observed further towards the 149 tumour's centre [START_REF] Larue | Microenvironmental regulation of proliferation in multicellular spheroids is mediated through differential expression of cyclin-dependent[END_REF]Mueller-Klieser, 150 2000). To model this we postulate that proliferating 151 cells tend to become quiescent if too far from the edge 152 of the tumour and quiescent cells tend to become prolif-153 erative if too close to the edge of the tumour. We also 154 assume that becoming quiescent does not affect a cell's 155 death rate. The resulting equations for the proliferating 156 and quiescent cell populations are 157

∂p ∂t + ∇ • vp = [k m (c) -k d (c)] p -αpH ([S tum -R] -r) + βqH (r -[S tum -R]) , ∂q ∂t + ∇ • vq = -k d (c)q + αpH ([S tum -R] -r) -βqH (r -[S tum -R]) ,
where p(r, t) and q(r, t) are the densities of proliferat-158 ing and quiescent cells respectively, v(r, t) is the local 159 cell velocity, S tum (t) is the radius of the tumour, R is 160 the width of the proliferating rim, α is the rate at which 161 proliferating cells become quiescent outside the prolif-162 erating rim, β is the rate at which quiescent cells become 163 proliferative inside the proliferating rim, and H(x) is the 164 Heaviside function (which takes the value 0 if x is neg-165 ative, and 1 if x is positive).

To investigate whether it is feasible for a significant gap to be maintained between the tumour and the vasculature as the tumour grows we also consider the dynamics of the non-cancerous cells immediately outside the tumour. We assume that non-cancerous cells proliferate more slowly than cancerous cells but otherwise exhibit similar dynamics; in particular we assume that their proliferation depends upon the local oxygen concentration in the same manner for both cell types. We also assume that necrotic cells undergo linear degradation. Then the evolution of non-cancerous cells and necrotic cells is given by

∂h ∂t + ∇ • (vh) = (λk m (c) -k d (c)) h, ∂n ∂t + ∇ • vn = k d (c)(h + p + q) -νn,
where h(r, t), and n(r, t) are the densities of noncancerous cells and necrotic cells respectively, λ < 1 is the relative proliferation rate of non-cancerous cells, and ν is the linear degradation rate of necrotic cells.

Note that the kinetics of healthy and cancerous cells can be chosen to allow the canceous cells to proliferate at lower oxygen concentrations than the non-cancerous cells; this does not significantly alter the modelling predictions and the model presented above constitutes the best case scenario for exploring gap maintenance. We estimate the relative proliferation rate of non-cancerous cells, λ, using the fact that in the absence of any tumour the non-cancerous cell population should be in equilibrium. Hence λ is chosen so that in the absence of any tumour, oxygen feedback regulates cell proliferation and necrosis in exactly the way required to maintain the observed intercapillary distance.

We close our model by assuming that the cell volume fraction of tissue is constant, which implies the follow-

∇ • v = V L (λh + p) k m (c) V + (V N -V L ) (h + p + q) k d (c) V - V N νn V .
where V L is the mean volume of a live cell, V N is the 199 mean volume of a necrotic cell, and V is the cell volume 200 fraction, and that oxygen diffuses out of the capillaries 201 and is taken up by living cells. If we assume that the rate 202 of oxygen uptake for the normal processes which keep 203 cells alive also follows Michaelis-Menten kinetics and 204 the rate of additional oxygen uptake required for cell 205 proliferation is proportional to the rate of proliferation, 206 then the evolution of the local oxygen concentration is 207 given by 208

∂c ∂t +∇•(vc) = ∇•(D∇c)-γk m (c)(λh+ p)-k n (c)(h+ p+q),
where c(r, t) is the local oxygen concentration, k n (c) is 209 the Michaelis-Menten kinetics for the oxygen uptake for 210 normal processes, γ is the constant of proportionality 211 for additional oxygen uptake for mitosis, and D is the 212 diffusion coefficient of oxygen. 213

We take the initial state of our tumour to be the largest 214 collection of cells which can be assumed to be entirely 215 proliferating, which is a tumour of radius equal to three 216 cell widths [START_REF] Larue | Microenvironmental regulation of proliferation in multicellular spheroids is mediated through differential expression of cyclin-dependent[END_REF]. We also assume that 217 the tumour front moves with the local cell velocity, so 218

that 219 dS tum dt = v(S tum , t),
as do the neighbouring capillaries. Note that if the capil-220 laries do not move with the local velocity vector field for 221 structural reasons the separation between tumour and 222 capillary would reduce even more quickly so that, once 223 more, we are considering a best scenario for gap main-224 tenance. Similar comments apply for the initial position 225 of the vasculature, which is assumed to be 50μm from 226 the tumour, consistent with the intercapillary distance [START_REF] Fait | Microvascular patterns of the human large intestine: Morphometric studies of vascular parameters in corrosion casts[END_REF]. For our final boundary condition we assume the capillary network is a constant source of oxygen.

Where possible we estimate parameters directly from the literature, but it is necessary to estimate some parameters by performing an a posteriori analysis of parameters. To do this we numerically simulate the evolution of the tumour for different values of the unknown parameters and refine the estimates until the modelled tumour growth is consistent with observed tumour growth. In particular we ensure that our modelled tumour's growth curve, size at saturation, relationship between necrotic core radius and tumour radius, and oxygen diffusion distance are all consistent with observed behaviour. A summary of all the parameters used in our model and their values is presented in Table 1.

(Note that the oxygen concentration on the edge of our domain can be non-dimensionalised out of our model and hence need not be estimated explicitly).

Given this model we now simulate the evolution of our tumour whilst paying particular attention to the distance between the tumour and the neighbouring capillaries.

Results

Figure 2 shows the growth dynamics for the modelled tumour. Figure 2(a) shows that our modelled tumour undergoes an initial period of exponential growth, then the growth rate slows and the tumour saturates at a radius of approximately 1 mm, which is consistent with observed tumour growth. Table 1: The parameters in our model. given the parameters in Table 1.

tumour and the neighbouring vasculature as the tumour grows assuming that there is no obstruction between the tumour and vasculature. This is likely to be valid when the tumour is a secondary metastasis, in contrast to a carcinogenesis of a primary carcinoma with a basement tercapillary distance, to evolve into a gap of 2000 μm, 312 which is observed following an intra-cornea tumour im-313 plant, requires non-cancerous cells to proliferate on av-314 erage approximately 160 times faster than this. 315

Discussion and Conclusions 316

Our modelled tumour showed that simple biologi-317 cally motivated dynamics for the proliferation of non-318 cancerous cells caused the distance between the tumour 319 and neighbouring vasculature to decrease quickly com-320 pared to the growth rate of the tumour, so that at satu-321 ration the neighbouring capillaries effectively lie on the 322 surface of the tumour, unless there is some obstruction 323 which prevents the tumour from reaching the surround-324 ing vasculature, such as a basement membrane. This 325 was despite taking care to ensure that any idealisations 326 in our model contributed to the model being a best case 327 scenario for maintaining the distance between the tu-328 mour and vasculature. 329

Non-cancerous cell proliferation can be forced to 330 maintain a constant gap between the tumour and vas-331 culature, but this requires non-cancerous cells to ini-332 tially proliferate approximately 55% as quickly as tu-333 mour cells, which are typically taken to have lost growth 334 control. This result is consistent with our initial scaling 335 argument, which motivated the problem, which showed 336 that maintaining a uniform gap between a tumour and 337 the vasculature requires the non-cancerous cells to ini-338 tially proliferate approximately 50% as quickly as the 339 cancerous cells. Whilst the rate of proliferation required 340 to maintain such a gap is high, it is not unfeasible. How-341 ever the absence of any biologically motivated mecha-342 nism which would up-regulate non-cancerous cell pro-343 liferation in this way, and having shown that oxygen feedback is not sufficient, suggests that non-cancerous cell proliferation is not generally up-regulated to this extent.

We have found that the rate of proliferation required to maintain a gap of constant width between the tumour and vasculature is high and demands extreme parameter values. However, for an in vivo tumour growing in vascular tissue to evolve into a state which is consistent with that seen following an intra-cornea tumour implant would require non-cancerous cells to proliferate on average approximately 160 times faster that was required to maintain a gap of constant width. This would imply that some non-cancerous cells would be proliferating at least two orders of magnitude faster than cancerous cells which is in conflict with our general understanding of tumour systems. This is clearly impossible. Furthermore, a 2 mm gap between the tumour and neighbouring vasculature would greatly exceed the diffusion distance of oxygen and it is not clear how an in vivo tumour in such a configuration would be supported. From this we conclude that the 2 mm gap following intra-cornea tumour implants is unique to that setting and cannot be included in a general model of tumour blood flow dynamics.

Our model assumes that capillaries are carried with the movement of local cells and neglects any angiogenic stimuli or mechanical properties of the capillaries which might resist their movement. We have also neglected any contact inhibition of normal, non-cancerous, cell growth. If any of these effects are included then the non-cancerous cell growth rates required for gap maintenance will be even higher. Hence our model is a best case scenario for gap maintenance and in an in vivo tumour we would expect the gap to close even quicker than it does in our model. Given this and that the gap closes quickly in our modelled tumour we conclude that in a well developed in vivo tumour growing in vascular tissue there will not be a significant gap between the tumour and neighbouring capillaries, except possibly for the presence of an intact basement membrane preventing the tumour from reaching the surrounding vasculature, as with the carcinogenesis of a primary carcinoma.

Unless the tumour is geometrically unconstrained our simple scaling argument still applies and the gap between the tumour and surrounding vasculature should only be the width of the basement membrane, which will be relatively thin. This conclusion implies that mathematical models of angiogenesis typically should not consider an initial gap between the tumour and the vasculature, especially for micrometastases. Furthermore, the studies based on corneal implants serve as building blocks for developing models for angiogenesis in that setting, but do not immediately translate to a clinical setting and are likely to overestimate the importance of angiogenesis and the effectiveness of anti-angiogens without suitable modification. Furthermore, the new vasculature cultivated by tumour angiogenesis tends to be distorted and leaky compared to normal vasculature though this will only exist within the tumour interior if angiogenesis initiates from the surrounding capillaries encroaching the tumour surface. This could mean that predictions of healthy tis- images of intra-corneal tumour angiogenesis.

  mour and neighbouring capillaries requires further validation, especially since it appears to demand extreme parameter choices in the context of a simple scaling argument. Alternatively, if no such gap exists at the onset of angiogenesis then this poses the question of what causes the gap to close.

Figure 1 :

 1 Figure 1: A schematic cross-section of the domain being considered.

Figure 2

 2 Figure 2(b) shows the width of the gap between the

  (a) The growth curve of our modelled tumour. (b) The width of the gap between the tumour and neighbouring vasculature as a function of tumour radius for our modelled tumour.

Figure 2 :

 2 Figure 2: The growth dynamics of our modelled tumour,

  membrane separating the surrounding vasculature from the microtumour. It can be seen that the width of the gap decreases significantly as the tumour grows and the capillaries effectively reach the surface of the tumour well before the tumour growth saturates.To further highlight this we compare the configurations of the initial state of our modelled tumour, the final state of our modelled tumour, and a typical tumour following an intra-cornea implant in Figure3. Each of the schematics in Figure3is on the same scale with a 1 mm scale bar shown. In Figure3(a) there is initially a small gap of width 50 μm between our modelled tumour and neighbouring capillaries. In Figure3(b) it can be seen that when our modelled tumour has grown to its saturated size, as given by diffusion limited control, there is effectively no gap between the tumour and the vasculature. However, in Figure3(c) it can be seen that following tumour implantation into the cornea there is a large gap between the tumour and the vasculature. The difference between Figures 3(b) and 3(c) is marked and highlights that the results of our model conflict with observed angiogenesis in the corneal implant experiments where capillaries are seen to sprout towards and into an avascular tumour. To investigate whether the lack of a gap between our modelled tumour and the neighbouring vasculature is due to a failing of our model, or whether this properly captures in vivo behaviour, we now suppose that noncancerous cell proliferation is somehow upregulated so as to maintain a constant gap between the tumour and (a) The initial configuration of our model. (b) The final configuration of our model (at diffusion limited saturation). (c) The configuration corresponding to a typical cornea implant experiment.

Figure 3 :

 3 Figure 3: A schematic comparison of the configurations of the tumour corresponding to (a) the initial conditions taken in our model, (b) the final state of our modelled tumour, and (c) a typical tumour following implantation into the cornea, each to scale with a 1 mm scale bar shown. The central shaded regions correspond to the tumour, the red lines outside the tumour represent the capillaries and any unshaded regions correspond to non-cancerous cells.

  sue exposure to cytotoxic chemotherapies during therapeutic interventions are overestimated within corneal implant motivated models because of differing detailed predictions of vasculature integrity. These conclusions are supported by experimental observations. The image of the avascular tumour in Figure 5(a) shows small blood vessels growing directly into the periphery of a tumour that was initially avascular. It 414 demonstrates that blood vessels penetrate directly into 415 the tumour-host interface rather than via an elaborate 416 network of vessels adjacent to the tumour edge as seen 417 in the cornea experiments, such as that shown in Fig-418 ure 5(b). In Figure 5(b) it can be seen that the corneal 419 implant is surrounded by an extremely large region of 420 avascular tissue and the resulting angiogenesis across a 421 large gap between the tumour and vasculature does not 422 represent angiogenesis for an in vivo tumour growing in 423 vascular tissue. 424 5. Summary 425 We have considered the evolution of the gap between 426 an avascular tumour and the neighbouring vasculature 427 as the tumour grows. A simple mathematical model 428 illustrates that a well established yet avascular tumour 429 will not exhibit a gap between the tumour and the vas-430 culature and that the neighbouring capillaries effectively 431 lie on the surface of the tumour. These conclusions are 432 reinforced by a simple scaling argument and images of 433 in vivo tumours. This illustrates that mathematical and 434 computational models based on the cornea implant ex-435 periments, which assume a significant gap between the 436 tumour and neighbouring vasculature, cannot be immeby a DTA award from the EP-440 SRC. PKM was partially supported by a Royal Society 441 Wolfson Research Merit Award. PKM and RAG were 442 partially supported by NIH grant 1U54CA143970-01. 443 The authors would like to thank Dr. T. Velpandian of 444 the all India institute of medical sciences for providing 445 (a) An in vivo avascular tumour (the white region) growing in a mouse mound chamber. The tumour has been added to the wound chamber in a matrigel disc, and on the outer ring of the disc endothelial buds which were obtained from scrotal fat were placed. A 200 μm length bar is shown. (b) Angiogenesis following tumour implantation into a test animal's cornea (Velpandian).

Figure 5 :

 5 Figure 5: A comparison of an in vivo tumour and a tumour following implantation into a test animal's cornea.