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

Alberto Gandolfi
email: gandolfi@iasi.cnr.it

Chemotherapy of vascularised tumours: Role of vessel density and the effect of vascular "pruning"

Keywords: Tumour, angiogenesis, chemotherapy, pruning, multistability

In this work we propose to model chemotherapy taking into account the mutual interaction between tumour growth and the development of tumour vasculature. By adopting a simple model for this interaction, and assuming that the efficacy of a drug can be modulated by the vessel density, we study the constant continuous therapy, the periodic bolus-based therapy, and combined therapy in which a chemoterapic drug is associated with an anti-angiogenic agent. The model allows to represent the vessel-disrupting activity of some standard chemotherapic drugs, and shows, in the case of constant continuous drug administration, the possibility of multiple stable equilibria. The multistability suggests an explanation for some sudden losses of control observed during therapy, and for the beneficial effect of vascular "pruning" exherted by anti-angiogenic agents in combined therapy. Moreover, in case of periodic therapies in which the drug amount administered per unit time is constant ("metronomic" delivery), the model predicts a response, as function of the bolus frequency, significantly influenced by the extent of the anti-angiogenic activity of the chemotherapic drug and by the dependence of the drug efficacy on the vessel density.

Introduction

The use of chemotherapies in oncology has been one of the major steps forward in the so called war against cancer [START_REF] Peckham | Oxford Textbook of Oncology[END_REF][START_REF] Boyle | Measuring progress against cancer in Europe: has the 15% decline targeted for 2000 come about?[END_REF]. Chemotherapy has a clinical role so important that in the common usage of many languages the word chemotherapy nowadays uniquely denotes anti-tumour chemotherapy. A huge corpus of experimental and clinical literature has been produced in the last 60 years, but chemotherapy also triggered the production of a large amount of theoretical researches due to its apparently simple translation into mathematical models (Skipper 1986, Wheldon 1988[START_REF] Agur | A theoretical analysis of interval drug dosing for cell-cycle-phasespecific drugs[END_REF][START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF], Swierniak et al. 1996[START_REF] Panetta | A mathematical model of breast and ovarian cancer treated with paclitaxel[END_REF][START_REF] Ledzewicz | Drug resistance in cancer chemotherapy as an optimal control problem[END_REF], Simeoni et al. 2004, 1 A c c e p t e d m a n u s c r i p t [START_REF] Ribba | A mathematical model of Doxorubicin treatment efficacy for non-Hodgkin's lymphoma: investigation of the current protocol through theoretical modelling results[END_REF]. Quite interestingly, and unlike from other fields of biomedicine, here there has been a limited but important interplay between theoretical and experimental-clinical scientists [START_REF] Goldie | A mathematic model for relating the drug sensitivity of tumors to their spontaneous mutation rate[END_REF][START_REF] Norton | The Norton-Simon hypothesis revisited[END_REF][START_REF] Panetta | Using pharmacokinetic and pharmacodynamic modeling and simulation to evaluate importance of schedule in topotecan therapy for pediatric neuroblastoma[END_REF], Ubezio and Cameron 2008).

Since their first applications, it has been plain that a number of serious side effects are related to the use of cytotoxic chemicals to cure tumors, for the simple reason that these agents, which never are perfectly selective, kill a more or less wide range of proliferating cells. However, even when the amount of side effects was relatively limited, a number of failures related to the emergence of resistance of tumour cells to the drug were noticed. Progresses in molecular biology allowed to explain this phenomenon: tumour cells are strongly characterized by a mutation rate significantly higher than that of healthy cells [START_REF] Schimke | Gene amplification, drug resistance, and cancer[END_REF][START_REF] Goldie | Drug resistance in cancer: a perspective[END_REF]). Thus, the response of tumour cells to chemotherapy is characterized by a considerable evolutionary ability to enhance the cell survival in an environment becoming hostile. Moreover, small subpopulations of cells intrinsically unsensitive to the treatment may be present ab initio in the tumour and, after some time, these cells may become the dominant population leading to the failure of therapy. Considerable research efforts have been and are devoted in finding means to overcome drug resistance [START_REF] Frame | New strategies in controlling drug resistance[END_REF].

Chemotherapy aims at targeting the main peculiar characteristic of tumour cells, i.e. their proliferative derangement. However, tumour cells show a vast array of genetic and epigenetic events, and of microscopic and macroscopic interactions with other cellular populations. As a consequence, the study of these phenomena may open the way for the creation of new therapies. J. Folkman [START_REF] Folkman | Anti-angiogenesis: new concept for therapy of solid tumors[END_REF][START_REF] Folkman | Tumor angiogenesis: a possible control point in tumor growth[END_REF]) stressed in early seventies that during the progression of tumours the development of a vascular network inside the tumour mass becomes necessary to support the tumour growth, and he named this process neo-angiogenesis. The tumour angiogenesis is a very complex phenomenon, driven by pro-angiogenic factors released by the tumour cells lacking a full level of nutrients. Interestingly enough, tumour cells also release anti-angiogenic chemicals to modulate the growth of the vessel network. In such a way, a solid tumour deploys a sophisticated strategy to control its own growth. However, Folkman did not only study the neo-angiogenesis. He also had the idea that controlling the development of the tumoral vessel network could be a powerful way to control, in turn, the neoplastic growth via the reduction of nutrients supply, and he termed this new kind of therapy 'anti-angiogenic therapy'. Per se, this new way of controlling the tumor burden appears intriguing. Moreover, it has another fascinating advantage respect to conventional chemotherapy: endothelial cells forming blood vessels are far more genetically stable than tumour cells, so that the anti-angiogenic agents should be less subject to resistance phenomena.

It is worth noting that tumours differ from normal tissues also in density, topology and functionality of their vessel network, which is characterized by a remarkable degree of intricacy as well as by a variety of disfunctionalities. The neovessel network that brings nutrients to the tumour also is the route to deliver chemotherapic drugs. As a consequence, the poor functionality of the network itself may be expected to affect the efficacy of chemotherapy negatively. Moreover, since in the angiogenesis process endothelial cells proliferate, chemotherapic drugs have been found to exert an anti-angiogenic action [START_REF] Browder | Antiangiogenic scheduling of chemotherapy improves efficacy against experimental drug-resistant cancer[END_REF]. Summarizing these observations, one may infer that:

• the theoretical characterization (e.g. by mathematical models) of chemotherapy must take into account the existence of neovessels and their dynamics;

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

• neovessels can be destroyed by chemotherapy;

• conversely, the excessive reduction of the vessel network by an anti-angiogenic agent might be detrimental for chemotherapy.

However, it has been observed that an anti-angiogenic agent has not necessarily a detrimental effect on the efficacy of a chemotherapic drug: on the contrary, better results have been obtained by combined therapy [START_REF] Browder | Antiangiogenic scheduling of chemotherapy improves efficacy against experimental drug-resistant cancer[END_REF][START_REF] Klement | Continuous low-dose therapy with vinblastine and VEGF receptor 2 antibody induces sustained tumor regression without overt toxicity[END_REF], Siemann and Rojiani 2002). To explain these findings, R.K. Jain hypothesized that the preliminary delivering of such an agent, by 'pruning' the vessel network, may regularize it with beneficial consequences for a succesive therapy with cytotoxic anti-tumour drugs [START_REF] Jain | Normalizing tumor vasculature with anti-angiogenic therapy: A new paradigm for combination therapy[END_REF][START_REF] Jain | Vascular normalization as a rationale for combining chemotherapy with antiangiogenesis agents[END_REF].

A full description of all these phenomena to create a comprehensive model of chemotherapy and of the combination of chemotherapy and anti-angiogenic therapy could require a multi-scale approach. Instead, we aim at assessing whether a relevant part of the interplay between a given chemotherapy, the neovessels and a possible anti-angiogenic therapy may simply be explained by ecological interaction between the two involved cellular populations: tumour celles and endothelial cells. Elementary models of this interaction, indeed, have been proved able to reproduce the main features of tumour growth and of anti-angiogenic therapies [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF], d'Onofrio et al. 2009a, Swierniak 2009).

The present paper is organized as follows: after a brief introduction of some basic concepts (Section 2), we summarize in Section 3 recent results on modelling the growth of a tumour and of its vasculature. In Section 4, we introduce our model of chemotherapy that assumes the dependence of the drug efficacy on the vessel density, and includes the possible vessel disrupting action of the chemotherapeutic drug. In Section 5, we study the effects of a constant uninterrupted chemotherapy, by stressing a number of phenomena linked to the possible multistability of the model. In the next Section 6, we study by simulations effects arising in the case of periodic boli-based therapy. In Section 7, we estend our model by allowing the co-presence of an anti-angiogenic therapy, and after studying the case of continuous constant therapy, we focus on the role of vascular "pruning". Final remarks conclude the work.

Preliminary concepts

The most simple model of chemotherapy can be written by coupling an empirical law for the unperturbed growth of the tumour, with a loss term depending on the drug concentration according to the so called "log-kill" hypothesis (Skipper 1986). If the insurgence of the drug resistance is disregarded, and a logistic law is assumed for the unperturbed growth, we can write:

N = α 1 - N (t) N ∞ N (t) -γc(t)N (t), (1)
where N (t) denotes the tumour size, N ∞ the (fixed) carrying capacity of the logistic growth, and α the proliferation rate constant. The coefficient γ represents the sensitivity of the cells to the drug, and c(t) is the drug concentration in blood.

In case of c(t) = C constant, which is the case of a continuous constant infusion treatment, it is easy to see from (1) that if C > α γ ,

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

N will tend to zero for any positive initial condition, i.e. asymptotic tumour eradication is achieved. Conversely, if C < α/γ there will exist a positive equilibrium smaller than N ∞ , precisely N e = (1 -γC/α)N ∞ . Note that the eradication would not be possible, whatever the magnitude of the drug concentration may be, if the growth law in (1) were assumed Gompertzian (d'Onofrio 2005). This fact is due to the unboundedness of the relative growth rate of the Gompertz law, N /N ∼ ln(N ∞ /N ), for N going to zero. Such an unboundedness clearly contradicts the physical basis of cell kinetics, since in any cell populations the relative growth rate has a finite upper bound because the cell cycle duration has a finite lower bound. This pathology was recognized by Wheldon (Wheldon 1988), who proposed the so-called Gomp-exp growth law, in which the growth is exponential for N < N , N given, and for N ≥ N it obeys to the Gompertz law. In the following, therefore, we will consider only tumour growth laws with bounded relative growth rate.

From the point of view of the cure, the really relevant measure of the tumour size is the number of viable cells, whereas when measurements on solid tumours are performed, usually the tumour volume is recorded. The two quantities have different dynamics [START_REF] Bertuzzi | Regression and regrowth of tumour cords following single-dose anticancer treatment[END_REF], Simeoni et al. 2004, Ubezio and Cameron 2008), since in the evolution of the tumour volume the delay in the dead cells degradation and in the reabsorption of waste fluids play an important role. In this paper, however, we will view for simplicity the number of viable tumour cells and the tumour volume as proportional quantities, and identify the tumour size with the tumour volume, V (t).

The appreciation of the role of angiogenesis in tumour development has led to the concept of a varying carrying capacity, defined as the tumour size potentially sustainable by the vascular network existing at a given time [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF]). The carrying capacity may be assumed proportional to the extent of the actual tumour vasculature [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF]. By introducing in Eq. (1) the variable carrying capacity, K(t), and by describing its dynamics according to [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF], we may obtain the model (see (d'Onofrio 2007, d'Onofrio et al. 2009b)):

V = α 1 - V (t) K(t) V (t) -γc(t)V (t), (2)
K = bV (t) -(dV (t) 2/3 + μ)K(t), (3)
where b, d, and μ (with b > μ) are constants related to the stimulation, inibition and natural loss of the vasculature respectively. In the case of a costant infusion therapy, the behaviour predicted by model ( 2),(3) is not substantially different from that of model (1). It can be easily derived that there exists a critical value for C above which there is eradication,

C > α γ b -μ b ,
and below which a positive equilibrium is established. Note that now the critical value of the drug concentration depends also on the parameters characterizing the dynamics of the vasculature. However, as we will see in the following, model (2),(3) can be a useful starting point towards extensions aimed at describing the possible dependence of the efficacy of the cytotoxic drug on the "quality" of the vascular network, the anti-angiogenic effect of some chemotherapeutic drugs, and the synergy between vessel targeting and tumour-cell targeting drugs.

A c c e p t e d m a n u s c r i p t 3 A family of models for the tumour growth

To describe the interplay between the tumour and its vasculature, we adopt a family of models previously proposed (d'Onofrio and Gandolfi 2009), which includes as particular cases the models in [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF][START_REF] Sachs | Simple ODE models of tumour growth and anti-angiogenic or radiation treatment[END_REF], d'Onofrio and Gandolfi 2004). In this family of models, we assume that (i) the carrying capacity of the tumour vasculature is simply proportional to the amount of vessels, and (ii) the specific growth rate of the tumour, V /V , and the specific "birth" rate of vessels depend on the ratio between the carrying capacity and the tumour size. Since the ratio K/V may be interpreted as proportional to the tumour vessel density, assumption (ii) agrees with the model proposed by [START_REF] Agur | Hopf point analysis for angiogenesis models[END_REF](see also [START_REF] Forys | Critical-point analysis for three-variable cancer angiogenesis model[END_REF]). Following [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF], the growth of the neo-vasculature will be antagonized by endogenous factors and we can write:

V = V F ( K V ) ( 4 ) K = K β( K V ) -ψ(V ) -μ (5)
where:

• F (u) > 0, F (1) = 0, 0 < lim u→+∞ F (u) < +∞, and 0 > lim u→0 + F (u) ≥ -∞; • β(+∞) = 0, β (u) < 0 and lim u→0 + β(u) = β 0 ≤ +∞;
• ψ(0) = 0, ψ (u) > 0 and lim u→+∞ ψ(u) = +∞;

• β(1) > μ.

In Eq. ( 5), ψ(V ) represents the vessels loss due to the possible accumulation into the tumour of an endogenous inhibitory factor secreted by the tumour cells, and μ represents the natural loss of vessels.

The assumption that F (+∞) < +∞ is the only restriction we set on the family of models proposed in (d'Onofrio and Gandolfi 2009), and it means that V /V will have a finite upper bound. Thus the Gompertz law is not included in F . As an example of possible expressions for F (u) we may consider the generalised-logistic:

F (u) = α(1 -u -ν ), ν >0. (6)
We mention here that a generalised-logistic growth with ν = 2/3 can be predicted for spherical tumours supplied by diffusible nutrients at the periphery [START_REF] Bodnar | Three types of simple DDE's describing tumour growth[END_REF], assuming the cell proliferation rate proportional to nutrient concentration and a uniform cell death. The function β(u) includes power laws

β(u) = bu -δ , δ >0, (7)
so that it may be lim u→0 + β(u) = +∞, functions such as

β(u) = β M 1 1 + ku n , n ≥ 1,

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

which have lim u→0 + β(u) finite and correspond to Hill functions in the variable u -1 , and combinations of the above two expressions as:

β(u) = β 1 1 1 + ku n + β 2 1 u . (8)
The expression ( 7) with δ = 1 yields Kβ(K/V ) = bV , as proposed by [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF], whereas (8) yields

Kβ( K V ) = β 1 1 + k( K V ) n -1 K + β 2 V.
The above expression for the growth term of the tumour vasculature distinguishes the contribution of the endothelial cell proliferation, which is assumed depending on the vessel density, and the input of new endothelial cells due to their migration from the peritumoral regions or to the influx of circulating endothelial progenitors [START_REF] Rafii | Vascular and haematopoietic stem cells: novel targets for anti-angiogenesis therapy?[END_REF]. The latter contribution is simply taken proportional to the tumour size. On the other hand, the prescribed properties of the function β exclude those modifications of the Hahnfeldt's model proposed in [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF]) and (d'Onofrio and Gandolfi 2004) that disregard the role of the tumour mass in the proantiangiogenic activity.

Concerning the function ψ, we recall that ψ(V ) = dV 2/3 has been assumed in [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF].

It is easy to see that model ( 4)-( 5) admits a unique positive equilibrium V e = K e , and it can be shown that this equilibrium is globally asymptotically stable in R 2 + (d'Onofrio and Gandolfi 2009).

As an example, Fig. 1 reports the tumour volume growth, predicted by the following model

V = α 1 -( V K ) 0.5 V, (9)
K = bV -dV 2/3 K, (10)
for different initial values of K. Note that the equation for K is that proposed in [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF] with μ = 0. The growth curve of mouse xenografts is taken as reference for this simulation. As the panels A and B show, different sets of parameters α and (b, d) can give very similar growth curves. When the dynamics of the vasculature is fast ((b, d) high, panel A), the system tends to reach quickly the K-nullcline, the evolution is thus scarcely influenced by the initial value of K, and the growth velocity is mainly under the control of the intrinsic proliferation rate α. Conversely, when the dynamics of the vasculature is rather low ((b, d) low, panel B), this dynamics starts to be limiting and thus to have a comparable growth velocity the proliferation rate of tumour cells must be higher, and the influence of the K(0) value is increased. Panels C and D show the time course of the ratio ρ(t) = K(t)/V (t). Since all the growths tend to an equilibrium state with V e = K e , ρ(t) always will tend to 1. In the case of higher α and lower (b, d) (Panel D), ρ goes to the unity more rapidly and monotonically.

The chemotherapy model

To model the effect of a cytotoxic chemotherapeutic agent of concentration c(t), we have to modify model ( 4)-( 5) by adding a log-kill term in equation ( 4), as well as a term representing the possible cytotoxic actions of the chemotherapeutic drug on endothelial cells.

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

The killing efficacy of a blood-born agent on the tumour cells will depend on its actual concentration at the cell site, and thus it will be influenced by the geometry of the vascular network and by the extent of blood flow. The efficacy of a drug will be higher if vessels are close each other and sufficiently regular to permit a fast blood flow; it will be lower if vessels are distanced, or irregular and turtuous so to hamper the flow. To represent simply these phenomena, we assume that the log-kill term to be added in Eq. ( 4) be dependent on the vessel density, i.e. in our model on the ratio ρ = K/V , by writing it as γ(ρ)c(t)V (t). We suppose γ(0) = 0, and γ increasing for ρ small. For larger values of the vessel density we make two hypotheses: either γ(ρ) continues to increase tending to a saturation value, or it starts to decrease after having reached a unique absolute maximum. According to these assumptions we propose the following model

V = V F ( K V ) -γ( K V )c(t) (11) K = K β( K V ) -ψ(V ) -μ -χc(t) (12)
where:

• c(t) ≥ 0 is the blood concentration of the cytotoxic drug;

• χ ≥ 0. The case χ = 0 corresponds to chemoterapic agents having no vessel disrupting effects;

• γ(ρ), with γ(0) = 0, is continuous and either increasing with γ(+∞) < +∞, or it exists ρ M > 0 such that at ρ = ρ M the function γ has an absolute maximum and γ(+∞) ≥ 0.

We shall call γ(ρ) the efficacy curve of the drug. A simple sufficient condition, although quite sharp, for the asymptotic eradication of the tumour is stated in the following:

Proposition 4.1 Let us define C o = sup ρ∈[1,+∞) F (ρ) γ(ρ) If C o < +∞ and min t∈[t 0 ,+∞) c(t) > C o for some t 0 ≥ 0, ( 13
)
then lim t→+∞ (V (t), K(t)) = (0, 0),
for any positive initial condition, i.e. global eradication of the tumor is achieved.

It is interesting to note that, as far as the role of γ(ρ) in the asymptotic behaviour is concerned, only the values of this function for ρ > 1 matter.

Constant infusion therapy: eradication and equilibrium points

Eradication

In the case of c(t) = C, i.e. in the case of constant continuous infusion therapy, the previous Proposition guarantees that if C > C o asymptotic eradication will occur. However, even if

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

C < C o , with C o ≤ +∞, there may be tumour eradication. In such a case, the anti-angiogenic effect of the drug as well as of the spontaneous loss of vessels play a necessary role. If 0 < C < C o , by setting V = 0, equilibrium values for the ratio ρ are determined by solving the equation

F (ρ) = Cγ(ρ), (14)
which, since (13) does not hold, has n ≥ 1 solutions:

1 < ρ 1 (C) < ρ 2 (C) < • • • < ρ n (C).
In the case of γ constant, obviously we would have n = 1. The corresponding equilibrium values V i (C) may be obtained by solving

ψ(V ) = β (ρ i (C)) -μ -χC, (15)
which, provided that

β (ρ i (C)) -μ -χC > 0,
has the unique positive solution:

V i (C) = ψ -1 (β (ρ i (C)) -μ -χC)
Of course, this suggests that if there exists

C * < C o such that β (ρ 1 (C * )) -μ -χC * = 0,
then there is no positive solution of (15), since β(ρ) is decreasing, and there is eradication. This fact is rigorously shown in the following proposition:

Proposition 5.1 Let a positive C * < C o exist such that β (ρ 1 (C * )) -μ -χC * = 0. ( 16
) If C ≥ C * , (17)
then the tumour is globally eradicated.

We remark that an eradication threshold C * < C o can exist only if χ > 0 and/or μ > 0, that is only if the chemotherapeutic drug exerts a vessel disrupting effect and/or a spontaneous loss of vessel is present. The converse is true if lim C→C o ρ 1 (C)) = +∞. In this case, in fact, if χ > 0 and/or μ > 0 a value C * < C o always exists, because β(ρ) tends to zero for ρ increasing.

In the case of delivering of a purely chemotherapeutic agent, i.e. χ = 0, the eradication condition becomes

β (ρ 1 (C)) -μ ≤ 0, (18)
and then the possibility of eradication with doses smaller than C o relies on the spontaneous loss of vessels. Since the rate of this loss is likely to be small [START_REF] Hahnfeldt | Tumour development under angiogenic signalling: A dynamic theory of tumour growth, treatment response, and postvascular dormancy[END_REF]), it appears unlikely that an eradication threshold significantly smaller than C o can exist in this case. Of course, if μ = 0 and C o = +∞ the eradication would be impossible.

A c c e p t e d

m a n u s c r i p t

Non eradication: multiple equilibria and their local stability

In the case in which the eradication condition (17) cannot be fulfilled, indicating with m ≤ n the maximum index i such that

β (ρ i (C)) -μ -χC > 0
we have m ≥ 1 coexisting equilibria:

E i = (V i (C), K i (C)), i= 1, . . . , m, with K i = ρ i V i .
Remark Note that the integer m is a function of C. In fact, in the case of nonmonotone γ(ρ), for small C there is typically a single equilibrium with ρ slightly greater than 1, for large C no or at most one equilibrium point is present (in this case the equation F (ρ) = γ(ρ)C can have one solution but there may be no equilibria because the eradication condition may be satisfied), whereas for intermediate C values there may be multiple equilibria.

From ( 13) we see that the tumour size at the equilibrium points under therapy, V i (C), i = 1, . . . , m, will be smaller than the equilibrium tumour size in the absence of therapy also if μ = χ = 0. We have in fact

β(1) -μ > β(ρ i (C)) -μ -χC,
and, since β(ρ) is decreasing, the above inequality is true also if μ = χ = 0. Whereas the role of the chemoterapic action appears to be, in the framework of the present model, the determination of a steady-state vessel density greater than 1 through the equation

F (ρ) = γ(ρ)C,
the eventual value of the tumour volume is prescribed by the vessel dynamics. The role of the anti-angiogenic action of the drug is that of setting the tumour size to a value lower than the value which would be set in the absence of any anti-angiogenic action.

It is interesting to note that the local stability of the co-existing equilibria depends on the dynamics of V and on the shape of the function γ(u) through a nice geometrical property, as shown in the following proposition:

Proposition 5.2 Let E i = (V i , K i ) be a critical point of (11)-(12). If γ (ρ i ) < F (ρ i ), ( 19
)
then E i is locally asymptotically stable, whereas if

γ (ρ i ) > F (ρ i ), ( 20
)
then E i is unstable.
In case of unique equilibrium point (m = 1), it is easy matter to show that it must be locally asymptotically stable. As far as its global stability is concerned, it holds the following:

C < min ρ∈(0,ρ M ) F (ρ) -β (ρ) γ (ρ) (21)
then E is globally asymptotically stable.

Figure 2 illustrates the response to a constant continuous infusion therapy of model ( 11)-( 12) with F , β and ψ as in ( 9)-( 10), and μ = 0. In the simulations we used the following simple form for γ(ρ):

γ(ρ) = max γ 4 9 ρ(3 -ρ), 0 , (22)
so that γ attains its maximum, equal to γ, for ρ = 3/2. With such a choice for the efficacy function of the drug, the treatment cannot be eradicative if the chemotherapic drug does not exert also an anti-angiogenic effect. However, for each value of C, a unique equilibrium exists, with a reduced tumour volume, and this equilibrium point can asymptotically be reached. In the absence of any vessel disrupting effect (χ = 0, solid lines), the final volume is smaller in the case of model 1, as expected because this model is characterized by a smaller value of the proliferation rate α of tumour cells. The response pattern can be changed if the drug also exert a vessel disrupting action (dashed lines): the final volumes are smaller, and the smaller volume is predicted now by model 2, since the slow vessel dynamics which is characteristic of such a model is impaired at a larger extent by the anti-angiogenic action. This simulation show how the underlying vessel dynamics can affect the response to a chemotherapeutic treatment of different tumours that appear very similar on the basis of their growth curves. We note that this behaviour is not necessarily related to the non-monotonic shape of γ(ρ), it can be achieved also for γ increasing or constant.

Bifurcations

In this section we shall deal with the analysis of bifurcations by taking as bifurcation parameter the drug concentration C. Preliminarly we note that, by applying the implicit function theorem, it follows that:

dρ i dC = γ(ρ i ) F (ρ i ) -γ (ρ i )C .
As a consequence, it holds that if i is odd then ρ i (C) > 0, else if i is even then ρ i (C) < 0. This implies that if m ≥ 3 then there is at least a hysteresis bifurcation [START_REF] Hale | Dynamics and bifurcations[END_REF].

In fact, let us consider for the sake of simplicity a C such that m( C) = 3. Note that this is possible both in the case of non-monotone γ functions, and for γ increasing. We have:

ρ 1 ( C) > 0 , ρ 2 ( C) < 0 , ρ 3 ( C) > 0,
(see middle panel of figure 3). It follows that by increasing C there is a threshold value C 2 > C such that at C = C 2 the equilibria E 1 ad E 2 will collide and only E 3 will last (see upper panel of figure 3 and figure 4), which implies that the asymptotic tumour size -which is given by eq. ( 15) -may undergo a 'jump' towards a far lower value. On the contrary, by decreasing C one has that it exists a C 1 < C < C 2 such that at C = C 1 E 2 and E 3 will collide and only E 1 will last (see lower panel of figure 3 and figure 4), with a possible jump of tumour size towards a higher value.

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

We note that the latter transition could occur even if C were constant and the efficacy function γ(ρ) were multiplied by a parameter continuously decreasing. A time-varying, decreasing γ might be the result of the development of drug resistance, and thus we might expect in such a case that the tumour control exerted by the chemotherapy can be abruptly lost.

Figure 5 illustrates this phenomenon. Let us consider the model ( 11)-( 12) with F , β and ψ as in ( 9)-( 10), μ = 0, and

γ(ρ, t) = γ 1 + ρ-ρm σ 2 e -(t-t 0 )/τ Heav(t -t 0 ), c(t) = CHeav(t -t 0 ),
where t 0 is the time at which the treatment starts, and Heav(•) is the Heaviside function. In the simulation, the tumour is treated from t = 15 day with a constant continuous infusion of drug. Its volume starts moving towards the equilibrium value characterized by the higher ρ value, and for a relatively long period tumour control is achieved with only a mild regrowth.

As the value of γ decreases, however, an abrupt transition occurs towards high values of the tumour volume, when the root of γ(ρ, t)C = F (ρ) jumps to a low ρ value. The tumour behaviour is qualitatively the same even in the presence of an anti-angiogenic effect of the drug.

Periodic therapy

In this section we shall assess the effects of a therapy periodically delivered by means of boluses. In the case of a chemotherapic agent with monoexponential pharmacokinetics, the drug concentration profile in blood is asymptotically given by

c(t) = D c e -qM od(t,T ) 1 -e -qT , ( 23
)
where T is the period of the delivering, q is the clearance rate of the drug, and D c is the amount of drug delivered as bolus over the distribution volume. Such a profile of c(t) is characterized by the following average value:

< c(t) >= D c qT . (24)
For purely anti-angiogenic therapies, experimental studies suggested that intensifying the delivering frequency by maintaining the total administered dose per unit time constant, gives a better response [START_REF] Drixler | Continuous administration of angiostatin inhibits accelerated growth of colorectal liver metastases after partial hepatectomy[END_REF][START_REF] Kisker | Continuous administration of endostatin by intraperitoneally implanted osmotic pump improves the efficacy and potency of therapy in a mouse xenograft tumor model[END_REF]. Such kind of drug scheduling has been called (with a slight abuse) metronomic, and the improving of the response with the increase of the delivering frequency is termed "metronomic effect". It has also been suggested that metronomic schedulings might also improve the results of chemotherapies [START_REF] Kerbel | The anti-angiogenic basis of metronomic chemotherapy[END_REF][START_REF] Orlando | Prolonged clinical benefit with metronomic chemotherapy in patients with metastatic breast cancer[END_REF]). However, Browder et al. (2000) showed that, at least in some cases, the metronomic effect is not the strict rule and that there is a somewhat optimal T for the reduction of the tumour volume. Thus, here, we aim at verifying whether our model is able to reproduce these effects.

In our simulations we used, as for Fig 2, model ( 11)-( 12) with F , β and ψ as in ( 9)-( 10), μ = 0, and γ(ρ) firstly given by ( 22). Moreover, since in metronomic therapy the dose delivered per time unit is constant, we impose

D c = qT M 0 , (25)
A c c e p t e d m a n u s c r i p t so that: < c(t) >= M 0 .

After setting α = ln 2/1.5 day -1 , b = 4.64, d = 0.01, and γM 0 = 0.2, we simulated the model response changing T and the extent of the anti-angiogenic action, i.e the value of χ. The model, driven by a periodic c(t), produces an asymptotically periodic V (t) and we computed the average tumour volume, V m =< V (t) >, when this regimen is in practice reached. Note that with the γ function given by ( 22), the continuous infusion therapy does not guarantee the eradication of the tumor whatever be the value of the drug concentration.

Figure 6 shows how the average asymptotic tumour volume depends on T for different values of χ, both in the case of constant γ, and in the case of drug efficacy changing with ρ. In the case of γ constant, V m is an increasing function of T for all the values of χ tested, i.e. there is full metronomic effect (upper panel of figure 6). Quite different is the behaviour in case of non-constant γ(ρ) (lower panel). For χ = 0, namely in the case of delivering of an agent with pure chemotherapic action, the average volume V m decreases as T decreases, i.e. there is full metronomic effect. By setting χ to small values, V m is again increasing with T (not shown) and, of course, it sets to lower values because of the synergy between the chemotherapic and the anti-angiogenic action. However, for χ/γ = 0.5 or greater, V m becomes non-monotone and it exists a T min such that the average volume has a minimum. Thus our model was able to mimick the metronomic effect in the absence of antiangiogenic action of drug, but also the presence of a partial metronomic effect characterized by an optimal delivering period, T min , when a significant anti-angiogenic effect is present. However, the difference between the values of V m for T = 1 and for T = T min results small. Thus, for χ/γ = 0.5 and for χ/γ = 0.65 it would be more correct to say that there is a more or less large (depending on χ) interval for T where V m is approximately constant and delivering with the highest frequencies is not advantageous. The effect of higher values of χ is reported in Fig. 7. If for χ/γ = 1 a (slight) anti-metronomic effect extends up to about T = 9 days, for higher χ values the therapy tends to be more and more characterized by a metronomic effect. It must be observed, however, that when χ/γ attains the high values of 5 and 7 the drug is exherting mainly an anti-angiogenic, vessel disrupting action, and a better efficacy of time-dense schedulings of pure anti-angiogenic drugs is indeed expected (d'Onofrio et al. 2009a, d'Onofrio andGandolfi 2009).

The effect of varying the period of a periodic therapy (keeping the average concentration constant) can be subtly modulated by the shape of the function γ(ρ). Assuming the form

γ(ρ) = γMin(1, ρ/A),
which has a monotone profile, Figure 8 shows how different behaviours can be achieved by varying the value of A. The brisk change going from A = 3 to A = 3.5 corresponds, in the constant continuous therapy, to a bifurcation due to the occurrence of three equilibrium points.

Finally, we note that the metronomic effect appears to be related basically to the nonlinear nature of the tumour growth law. A slight metronomic effect (not shown) was indeed found using the model of Eq. (1). However, modelling the vessel dynamics and the dependence of the drug cytotoxicity on the vessel density makes it possible to expand the variety of the response patterns.

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

In the case in which the eradication condition is not fulfilled, indicating with m ≤ n the maximum index i such that

β (ρ i (C)) θ(H) -μ -χC -ηG > 0,
we have m coexisting equilibria

E 1 = (V 1 (C, G, H), ρ 1 (C)), . . . , E m = (V m (C, G, H), ρ m (C)),
whose local stability properties are established by Proposition (5.2). Note that the number m itself will depend on C, G, H. As we have seen before, if a unique equilibrium point there exists, it will be locally asymptotically stable. Concerning its global stability, it holds the following:

Proposition 7.2 Let E be a unique equilibrium point for Denoting as ρ M ≤ +∞ the value such that γ(ρ) is maximum, if

C < min ρ∈(0,ρ M ) F (ρ) -β (ρ)θ(H) γ (ρ) , ( 31
)
then E is globally asymptotically stable.

We note that in the case of a cytostatic anti-angiogenic agent and in the absence of any vessel disrupting activity by both the drugs (θ(H) < 1 and χ = η = 0), the eradication condition for

C < C o becomes θ(H)β (ρ 1 (C)) -μ ≤ 0, that is H > θ -1 μ β (ρ 1 (C)) . (32)
Thus, the tumour eradication might be difficult if μ is small, unless θ(H) has a quite rapid decrease. Of course, if μ = 0 the eradication would be impossible. In the case of an antiangiogenic drug which is purely vessel disrupting (θ(H) = 1 and η > 0), the eradication for C < C o is possible even if the chemotherapeutic drug does not impair the vasculature, provided that:

G > G * = β (ρ 1 (C)) -μ χ .
As far as the asymptotic outcome of a constant infusion therapy is concerned, the administration of an anti-angiogenic drug, therefore, always produces a synergistic effect. If chemotherapy alone were not eradicative, the concomitant administration of an anti-angiogenic drug could induce tumour eradication, or could set the asymptotic tumour size to a value lower than the value reachable by using the chemotherapeutic drug alone.

Complex dynamics can arise, in the case in which multiple equilibria are present. In such a case, as the following simulations will show, the asymptotic outcome of a constant continuous infusion therapy can be dramatically altered if the administration of the chemoterapic drug follows a preliminary transient vessel disrupting treatment. This behaviour matches, at least qualitatively, experimental observations on the positive effect of preliminary anti-angiogenic treatments on chemotherapy [START_REF] Jain | Normalizing tumor vasculature with anti-angiogenic therapy: A new paradigm for combination therapy[END_REF]. For these observations, the causal role of the

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

vascular "pruning", in order to obtain a "normalized" more efficient vascular network deprived of turtous and immature vessels, has been hypothesized [START_REF] Jain | Vascular normalization as a rationale for combining chemotherapy with antiangiogenesis agents[END_REF].

The simulations of Figures 9 and 10 are performed assuming

γ(ρ) = γ 1 + ρ-ρm σ 2 ,
and

c(t) = Heav(t -t 2 )C, g(t) = -t 1 ) -Heav(t -t 2 ))G,
where t 1 and t 2 are the times at which the anti-angiogenic treatment and, respectively, the chemotherapic treatment start. Figure 9 show how a mild chemotherapic treatment, capable if administered alone of only stopping the growth, become able, if proceeded by a transient (although rather vigorous) anti-angiogenic treatment, of achieving a remarkable tumour reduction. As the phase-portraits illustrate, the vessel disruption, that means the reduction of ρ, allows to exit from the basin of attraction of the equilibium point with large tumour volume, and to enter the basin of the second locally stable equilibrium, characterized by a smaller tumour volume. It may be interesting to note that this shift can be achieved also by a very short anti-angiogenic treatment (pulse treatment), provided that its intensity be sufficiently high (see Figure 10). Finally, we want to observe that:

• If in the case of three roots ρ 1 < ρ 2 < ρ 3 , the root ρ 3 corresponds to tumour eradication, the preliminary anti-angiogenic treatment can allow to exit from the basin of attraction of the equilibium point with larger tumour volume, and to reach eradication.

• The shift of the state of the system to a more favourable basin of attraction can be achieved, in principle, also by a purely antiproliferative anti-angiogenic treatment (H > 0, G = 0).

Concluding remarks

Despite its inherent simplicity, our model appears able to explain a wide array of biomedically significant behaviours. Most traditional models of chemotherapy are characterized by monostability, namely they suggest that the effect of delivering a chemotherapic drug is either to eradicate the tumor or to drive the tumour into a unique equilibrium state, so that varying slightly the delivering of the chemotherapic agent would in any case have the effect of slightly moving the equilibrium value. Although allowing some nice biological inferences, this appears an excessively simple scenario, largely different from the clinical reality where more complex and puzzling phenomena are observed (see e.g. [START_REF] Peckham | Oxford Textbook of Oncology[END_REF]). Thanks to the possibility of being multistable, our model may represent an improvement in the understanding of relevant biological phenomena and also, at least conceptually, in suggesting therapeutic plans.

In particular, the multistability of our model may give a contribution in understanding some cases of failures of therapy, in which the tumour suddenly restarts growing. These occurrences might simply be explained as hysteresis bifurcations, driven by the gradual onset of drug resistance, that cause jumps in V (t).

Moreover, we have shown how multistability and the beneficial effect of vascular "pruning" by anti-angiogenic agents can be strictly interlinked. In fact, if the tumour under a given (constant continuous) chemotherapy admits two stable equilibria, it appears in principle possible by the delivery of an anti-angiogenic treatment before the chemotherapy, to move the state in the basin of attraction of the equilibrium with the smaller (or null) tumour volume. In general (i.e also in the case of monostable behaviour), the model predicted synergy between constant continuous chemotherapy and anti-angiogenic treatment, and we were not able to find combinations in which the anti-angiogenic action was detrimental for the chemotherapy. However, we cannot exclude a priori that particular cases may exist where competition can arise. In fact, although the synergy between antiangiogenic therapy and chemotherapy may appear as a trivial intrinsic feature of our family of models, this property cannot strictly be considered as embedded in our dynamical systems because they are not cooperative. Furthermore, our model is able to describe how the effect of a periodic delivery of a cytotoxic drug can change by varying the time intervals between the administrations, keeping constant the average drug concentration. Our simulations show indeed that the vessel-disrupting action of some chemotherapic agents may imply the existence of an optimal delivery frequency that minimize the average tumour volume. This is a result of some interest also when the minimum is quite flat, since in such a case our simulations suggest that there is a quite large range of inter-boli time intervals where there is no sensible metronomic effect, i.e. an advantage in increasing the bolus frequency. In addition, the behaviour under a periodic therapy that keep constant the average drug concentration appears to be modulated, in a rather complex way, by the dependence of the drug efficacy on the vessel density.

We stress that we have considered here only tumour growth laws with bounded relative growth rate, so excluding the Gompertz law since its potentially unbounded relative growth rate is unphysical. An interesting property of models of anti-angiogenesis therapy incorporating a Gompertzian tumour growth law, namely the impossibility, in case of periodic therapy, of setting a sufficient condition for the tumour eradication only in terms of the average value of the drug concentration, is fully retained by models with generalised logistic growth law and ν < 1, for which the function F (K/V ) is bounded (d'Onofrio et al. 2009a, d'Onofrio andGandolfi 2009). From a strict mathematical point of view, however, the results of this paper on the tumour eradication in the presence of vessel disrupting activity and those on the location, local stability and bifurcation of multiple equilibria continue to hold also in case of F unbounded.

The model assumes that the carrying capacity of the vasculature is simply proportional to the vessel amount. A more realistic view might consider the carrying capacity dependent on the vessel amount but also on the vessel density. The same volume of vessels can have, in fact, a different capability of supplying oxygen and nutrients to the tumour cells because of different vessel functionality. We intend to analyse the consequence of this hypothesis in a future investigation. Preliminary results suggest that if the carrying capacity of vessels, for a fixed tumour volume, is an increasing function of the vessel amount, the qualitative behaviour of the present model should be substantially conserved.

Of course our model is oversimplified in many other important aspects. We must mention the distinction between mature vessels and immature vessels, and the spatial heterogeneity of the vessels network, with the consequent heterogeneity of the distribution of cells between the proliferating and quiescent compartments. This last phenomenon could be significant in the case of drugs characterized by a marked cycle-specificity. Moreover, the onset of drug resistance is of paramount importance in all chemotherapies. To account for the above Proof Preliminarly, for the sake of notation simplicity, let us define the following auxiliary function

H(ρ) = F (ρ) γ(ρ) . Since it is ∂V ∂V E i = -ρ i γ(ρ i )H (ρ i ), ∂V ∂K E i = γ(ρ i )H (ρ i ), ∂K ∂V E i = -ρ 2 i β (ρ i ) + K i Ψ (V i ) , ∂K ∂K E i = ρ i β (ρ i ),
we get the following characteristic polynomial:

λ 2 + ρ i γ(ρ i )H (ρ i ) -β (ρ i ) λ + K i Ψ (V i )γ(ρ i )H (ρ i ).
Thus if H (ρ i ) < 0 then E i is unstable, whereas if H (ρ i ) > 0 then E i is locally asymtptically stable (remember that β < 0). Finally, since

H (ρ i ) = γ(ρ i )F (ρ i ) -γ (ρ i )F (ρ i ) γ 2 (ρ i ) = F (ρ i ) -γ (ρ i ) γ(ρ i ) , our claim immediately follows. ♦ Proposition 5.3 Proof From div 1 V K V , K = 1 V 2 -F (ρ) + β (ρ) + γ (ρ)C .

  Let E be a unique equilibrium point for (11)-(12). Denoting as ρ M ≤ +∞ the value such that γ(ρ) is maximum. If

  , a far more complex modelling would be required. As a final remark, we observe that in this initial work we focused on general properties and so we framed our study in a general setting, showing that non-trivial results -in particular the possible presence of multistability -are substantially independent of the choice of specific functions. A comprehensive analysis of the dependence of the tumour response on the specific functional form of the rates F , β, γ and Ψ, might reveal finer details of the problem under study.

  we have, from the Poincare's tricothomy, that condition (21) implies the global asymptotic stability of E. ♦ Simeoni M.,Magni P., Cammia C., De Nicolao G., Croci V., Pesenti E., Germani M., Poggesi I., Rocchetti M., 2004

Figure 1 :

 1 Figure 1: Tumour volume (panels A and B), and ρ = K/V (panels C and D) vs. time, according to the model of Eqs. (9)-(10). Panels A and C: α = ln 2/1.5 day -1 , b = 4.64, d = 0.01. Panels B and D: α = ln 2/0.9 day -1 , b = 0.464, d = 0.001. In all panels: V (0) = 200 mm 3 , ρ(0) = 2 (dashed), ρ(0) = 4 (solid), ρ(0) = 8 (dotted line).

Figure 2 :Figure 3 :

 23 Figure 2: Tumour volume (panels A and B), and ρ = K/V (panels C and D) vs. time, before and after a treatment with constant drug concentration starting at t = 30 days. Panels A and C: α = ln 2/1.5 day -1 , b = 4.64, d = 0.01, γC = 0.2 day -1 , χ = 0 (solid lines). Panels B and D: α = ln 2/0.9 day -1 , b = 0.464, d = 0.001, γC = 0.2 day -1 , χ/γ = 0.65 (dashed lines). In all panels: V (0) = 200 mm 3 , ρ(0) = 4.

Figure 4 :

 4 Figure 4: Bifurcation diagram. In red the central unstable equilibrium, in black the two locally asymptotically stable equilibria. Functions F and γ as in Fig. 3.

Figure 5 :

 5 Figure 5: Tumour response to a constant continuous therapy, when Eq. (12) has 3 positive roots at the start of treatment, and drug efficacy decays with time (mimicking the rise of drug resistance).Model of Eqs. (9),(10), with γ(ρ, t)c(t) = 0.5Heav (t -15) 1/(1 + ((ρ -2)/0.35) 2 ) Exp (-0.01(t -15)). Other parameters as in the text.

Figure 6 :

 6 Figure 6: Periodic therapy. Average tumour volume V m vs. the period T for different values of χ. Upper panel: γ(ρ) constant; lower panel: γ(ρ) as in formula (22). In both panels: χ = 0 (solid line), χ/γ = 0.5 (dashed line), χ/γ = 0.65 (dotted line). Other parameters given in the text.

Figure 7 :

 7 Figure 7: Periodic therapy. Average tumour volume V m vs. the period T for different values of χ. γ(ρ) as in formula (22). χ/γ = 1 (solid line), χ/γ = 5 (dashed line), χ/γ = 7 (dotted line). Other parameters as Fig.6.

Figure 8 :

 8 Figure 8: Periodic therapy. Average tumour volume V m vs. the period T . Pure chemotherapy (χ = 0) with γ(ρ) = γMin(1, ρ/A). A = 1 (solid line), A = 3 (dashed line), A = 3.25 (dotted line), A = 3.5 (thick line). Other parameters as in Fig. 6.

  . Predictive pharmacokinetic-pharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents. Cancer Res. 64, 1094-1101. Swierniak A., 2009. Comparison of six models of antiangiogenic therapy. Appl. Math. 36,333-348. Ubezio P., Cameron D., 2008. Cell killing and resistance in pre-operative breast cancer chemotherapy. BMC Cancer 8, 201. Wheldon T., 1989. Mathematical Models in Cancer Research. Hilger Publishing.

	Skipper H.E., 1986. On mathematical modeling of critical variables in cancer treatment (goals:
	better understanding of the past and better planning in the future). Bull. Math. Biol.
	48,253-278.
	Swierniak A., Polanski A., Kimmel M., 1996. Optimal control problems arising in cell-cycle-
	specific cancer chemotherapy. Cell Prolif. 29, 117-139.

Acknowledgements

The authors want to thank Professors Urszula Ledzewicz and Heinz Schättler for useful discussions and comments, as well as two anonymous referees that helped us to improve this work. The work of one of the authors (A. d'O.) has been done in the framework of the Integrated Project "Advancing Clinical-Genomic Trials on Cancer -ACGT". The ACGT project (FP6-2005-IST-026996) is partly funded by the EC and the authors are grateful for this support.

A c c e p t e d m a n u s c r i p t 7 Combined chemotherapy and anti-angiogenic therapy. Role of vascular "pruning"

In this Section we will consider the case of a combined therapy in which a true anti-angiogenic agent is administered together with a chemotherapeutic drug. The most general setting is the following:

• a chemotherapeutic agent, with concentration c(t);

• an anti-angiogenic agent exerting vessel disruption, with concentration g(t);

• an anti-angiogenic agent, with concentration h(t), exerting (directly or indirectly) a cytostatic action on the endothelial cells.

Of course, the same drug can exert both vessel disruption and inhibition of endothelial cell proliferation, and in this case it shall be g(t) = h(t).

The model ( 11),( 12) can be easily extended, obtaining

where θ(u) is a decreasing function such that 0 ≤ θ ∞ < θ(u) ≤ 1. θ(u) = 1 corresponds to the case of absence of cytostatic action on the endothelial cells. The coefficient η ≥ 0 represents the vessel disrupting efficacy of the anti-angiogenic agent.

As it easy to recognize, almost all the results of the previous sections can be translated for the model ( 26),(27). We state the main results in the case of constant continuous administration of all the drugs, i.e when c(t) = C, g(t) = G, and h(t) = H.

If C o < +∞ and C > C o , asymptotic eradication of tumour is achieved thanks to Proposition 4.1. If 0 < C < C o ≤ +∞, the equation

where ρ i (C) is a generic solution of Eq. ( 14), has the unique positive solution

then there is no positive equilibria and there is asymptotic eradication. This property is stated in the following proposition:

then the tumour is globally eradicated.

10 Appendix: Proofs Proposition 4.1

Proof Rewriting (11) as follows

and taking into account that F (ρ) < 0 for ρ < 1 and F (ρ) > 0 for ρ > 1, if (13) holds we have V < 0 for any positive V . Thus it follows that lim t→+∞ V (t) = 0.

Thus the equation for K asymptotically becomes: 17) is satisfied, we have

This means that in the set

it is K > 0, implying that either V (t) → 0 + and minlim t→+∞ ρ(t) > ρ 1 , or that the orbit enters in the set:

which is positively invariant and where V < 0, so that V (t) → 0 + . ♦ . At t = 30 days, the tumour is treated with a vessel disrupting agent for 10 days (solid lines) or 30 days (small dashed lines). Thereafter, the anti-angiogenic treatment ceases and the tumour is treated with a pure chemoterapic drug. The effect of the chemotherapic drug alone, given at t = 30, is also shown (dashed line). In panels A and B the equilibrium points are reported as large dots (the unstable equilibrium is red). Model of Eqs. ( 9)-( 10


Masanori Takezawa
email: m.takezawa@uvt.nl

Michael E Price
email: michael.price@brunel.ac.uk

Revisiting "The evolution of reciprocity in sizable groups": Continuous reciprocity in the repeated Nperson prisoner's dilemma

Keywords: evolution of cooperation, reciprocal altruism, reciprocity, collective action, assortative interaction

For many years in evolutionary science, the consensus view has been that while reciprocal altruism can evolve in dyadic interactions, it is unlikely to evolve in sizable groups. This view had been based on studies which have assumed cooperation to be discrete rather than continuous (i.e., individuals can either fully cooperate or else fully defect, but they cannot continuously vary their level of cooperation). In real world cooperation, however, cooperation is often continuous. In this paper, we re-examine the evolution of reciprocity in sizable groups by presenting a model of the n-person prisoner's dilemma that assumes continuous rather than discrete cooperation. This model shows that continuous reciprocity has a dramatically wider basin of attraction than discrete reciprocity, and that this basin's size increases with efficiency of cooperation (marginal per capita return). Further, we find that assortative interaction interacts synergistically with continuous reciprocity to a much greater extent than it does with discrete reciprocity. These results suggest that previous models may have underestimated reciprocity's adaptiveness in groups. However, we also find that the invasion of continuous reciprocators into a population of unconditional defectors becomes realistic only within a narrow parameter space in which the efficiency of cooperation is close to its maximum bound. Therefore our model suggests that continuous reciprocity can evolve in large groups more easily than discrete reciprocity only under unusual circumstances.

Introduction

Direct reciprocity has long been regarded as a plausible route for the evolution of cooperation in genetically unrelated dyads [START_REF] Axelrod | The evolution of cooperation[END_REF][START_REF] Trivers | The evolution of reciprocal altruism[END_REF], particularly in humans. However, human sociality extends far beyond dyadic cooperation, and is characterized by collective action [START_REF] Olson | The logic of collective action: Public goods and the theory of groups[END_REF], i.e., a large group of unrelated members involved in the joint production of a shared resource. Based on the pessimistic results of models of the evolution of reciprocity in large groups (Bender & Mookherjee, 1987;Boyd & Richerson, 1988;[START_REF] Joshi | Evolution of cooperation by reciprocation within structured demes[END_REF][START_REF] Taylor | Anarchy and cooperation[END_REF], the standard view in biology has for years been that reciprocity is unlikely to evolve in groups much larger than dyads (e.g. [START_REF] Boyd | The evolution of altruistic punishment[END_REF][START_REF] Fehr | Don't lose your reputation[END_REF][START_REF] Fehr | The nature of human altruism[END_REF][START_REF] Gächter | Reciprocity, culture and human cooperation: previous insights and a new cross-cultural experiment[END_REF][START_REF] Hagen | Game theory and human evolution: A critique of some recent interpretations of experimental games[END_REF][START_REF] Hauert | Volunteering as red queen mechanism for cooperation in public goods games[END_REF][START_REF] Henrich | Cultural group selection, coevolutionary processes and large-scale cooperation[END_REF][START_REF] Henrich | Cooperation, punishment, and the evolution of human institutions[END_REF][START_REF] Kurzban | Experiments investigating cooperative types in humans: A complement to evolutionary theory and simulations[END_REF][START_REF] Sigmund | A tale of two selves[END_REF][START_REF] Suzuki | Reputation and the evolution of cooperation in sizable groups[END_REF]. This view has led many theorists to suggest alternative routes for the evolution of group cooperation, for example genetic or cultural group selection [START_REF] Boyd | The evolution of altruistic punishment[END_REF][START_REF] Gintis | Strong reciprocity and human sociality[END_REF][START_REF] Henrich | Cultural group selection, coevolutionary processes and large-scale cooperation[END_REF][START_REF] Wilson | Reintroducing group selection to the human behavioral[END_REF], costly signalling [START_REF] Gintis | Costly signaling and cooperation[END_REF], opting out of a group [START_REF] Hauert | Via freedom to coercion: The emergence of costly punishment[END_REF], and indirect reciprocity [START_REF] Panchanathan | Indirect reciprocity can stabilize cooperation without the second-order free rider problem[END_REF].

Reciprocity in dyads is traditionally modelled as the tit-for-tat strategy [START_REF] Axelrod | The evolution of cooperation[END_REF] that continues cooperating as long as the partner cooperates, but stops cooperating if the partner did not cooperate in the previous round. Tit-for-tat is a discrete strategy that contributes either fully or not at all. Models of reciprocity in groups have defined reciprocity as a discrete strategy that continues cooperating as long as all other A c c e p t e d m a n u s c r i p t 3 group members cooperate, but stops cooperating if one or more co-members defect (Bender & Mookherjee, 1987;Boyd & Richerson, 1988;[START_REF] Joshi | Evolution of cooperation by reciprocation within structured demes[END_REF][START_REF] Taylor | Anarchy and cooperation[END_REF].

However, the evidence that people actually engage in a discrete all-or-nothing strategyalso called a trigger strategy -is weak [START_REF] Ostrom | Rules, games, and common-pool resources[END_REF][START_REF] Watabe | Choice of strategies in a social dilemma[END_REF][START_REF] Watabe | Choice of strategies in social dilemma supergames[END_REF]. Furthermore, many instances of real-world cooperation seem to be better modelled as continuous rather than discrete, because individuals continuously vary their degree of cooperation from full defection to full cooperation. Bshary & Bronstein (2004), for example, review real world examples of interspecific mutualism and find many instances of continuous cooperation. Continuous cooperation is also common in human groups: members tend to modulate their contribution levels in order to approximately match the mean co-member contribution [START_REF] Croson | Theories of commitment, altruism and reciprocity: evidence from linear public goods games[END_REF][START_REF] Croson | Reciprocity, matching and conditional cooperation in two public goods games[END_REF][START_REF] Fischbacher | Are people conditionally cooperative? Evidence from a public good experiment[END_REF][START_REF] Kurzban | Experiments investigating cooperative types in humans: A complement to evolutionary theory and simulations[END_REF][START_REF] Kurzban | Incremental commitment and reciprocity in a real time public goods game[END_REF].

The representation of reciprocity as a discrete strategy could be regarded as a trivial or inevitable simplification of a complex reality, of the kind that is often required in formal models. However, we found that replacing discrete reciprocity with continuous mean-matching reciprocity effected significant changes in model results. In Section 2 below, we first review a model of discrete reciprocity in the n-person prisoner's dilemma (Boyd & Richerson, 1988) that is widely cited as having shown the difficulty of the evolution of (discrete) reciprocity in large groups. We then compare this model with a new model that replaces discrete reciprocity with continuous reciprocity. We find that like discrete reciprocity, continuous reciprocity can be evolutionarily stable against an unconditional defection strategy, while unconditional defection can also be evolutionarily

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

4 stable against both reciprocity strategies. The most important finding is that the basin of attraction for continuous reciprocity gets wider as cooperation gets more efficient, while the basin of attraction for discrete reciprocity remains quite small even when cooperation is extremely efficient. We also investigate the influence of assortative interaction and find it to have a strong synergic effect with continuous reciprocity but not with discrete reciprocity. These results hold even when an unconditional cooperation strategy exists in the population.

However, we find that the advantage of continuous reciprocity is strictly constrained by the efficiency of cooperation; the invasion of continuous reciprocity into a population of defectors becomes realistic only when cooperation is extremely efficient, i.e. when a contribution from a single individual produces a very large benefit for the entire group. Because conditions of such efficiency are probably relatively rare, this finding suggests that in a population divided up into large groups, continuous reciprocity could invade only under unusual circumstances. In Section 3, we demonstrate the robustness of continuous reciprocity in highly efficient public goods situations by extending the strategy space from one in which continuous reciprocators can only match the mean partner contribution to one in which they can contribute above or below this mean.

Discrete versus Continuous Reciprocity

Repeated n-person Prisoner's Dilemma Game

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

To model the evolution of cooperation in sizable groups, we considered a large population subdivided into randomly-formed groups of size n which play repeated nperson prisoner's dilemma games. Every round, individual members decide whether to contribute c to their group, in order to create the benefit Bc which is divided equally among all n group members. B is assumed to be smaller than n, so a member's private return is always less than that member's contribution (Bc/n < c), and groups are thus public good-producing collective actions characterized by social dilemmas. The game is repeated with the probability w.

The evolution of discrete reciprocity

First, let us consider two classes of strategies: unconditional defector (D) and discrete reciprocator (T a ). D never contributes, while T a contributes fully (c = 1) in the first round and continues contributing fully as long as a or more group co-members also contribute fully. Let V(x|y) be the payoff to an individual with a strategy x in a group with y reciprocators (in this section, T a ). When D is common, the expected payoff to D is V(D|0) = 0 while the payoff to a rare T a is V(T a |1) = B/n -1 when a = n -1, and (B/n -1)

/ (1 -w) when a < n -1. In both cases, the payoff to T a is smaller than 0 and thus a rare = B (n -1) / n and this is smaller than the expected payoff to T n-1 when w is sufficiently large (Boyd & Richerson, 1988). Thus, T n-1 is an evolutionary stable strategy that prevents the intrusion of rare D when interactions are sufficiently iterated. In the following, we will consider only T n-1 (hereafter referred to simply as T).

It is proved that, when T is evolutionarily stable, there is a unique unstable internal equilibrium (Boyd & Richerson, 1988). However, the basin of attraction for a pure T equilibrium gets smaller quickly as group size increases. The proportion of T at an

internal equilibrium is p 1 B /n w(B 1) /(1 w) w w 1/ n
. Figure 1a illustrates how this proportion of T changes as a function of the efficiency of cooperation B/n (marginal per capita return) and w when group size is 100. i The upper two lines in Figure 1a show that T cannot increase its share in a population unless it already composes a large majority of that population; if the proportion of T falls below 85-90%, D increases its share.

The evolution of continuous reciprocity

We now replace the discrete strategy T with a simple continuous reciprocal strategy R.

Like T, R contributes fully in the first round. R then matches the average contribution made by the other n-1 players in the previous round. When D is common, the expected payoff to D is V(D|0) = 0. The expected payoff to rare R is V(R|1) = B/n -1 because R contributes fully only in the first round and contributes nothing from the second round.

As V(R|1) < V(D|0), D is evolutionarily stable. On the other hand, when R is common, it is proved that the expected payoff to R, V(R|n) = (B -1) / (1 -w), is larger than D's expected payoff, V(D|n -1), and R is evolutionarily stable when The basin of attraction for R is drastically larger than that for T, however. Figure 1a shows the minimum proportion of R and T that must exist in a population in order to prevent the intrusion of D (i.e., the proportion of R and T at an unstable internal equilibrium). In contrast to T, the proportion of R is generally small, especially when the efficiency of cooperation, B/n, is high. For instance, when n = 100 and B/n is 0.95 (so 1 unit of cooperation [c = 1] produces 95 units of benefit [B = 95]), R can proliferate even when its initial proportion is only 10% or less. This is markedly different from T, which must compose 85-95% of the population in order to evolve. We increased group size in increments of 100 from 100 to 500, and confirmed that this advantage of R over T holds across all these group sizes. Figure 1c shows how the basin of attraction for T and R changes in response to the size of benefit of cooperation (B) and group size (n); the basin of attraction for R changes almost linearly as benefit of cooperation approaches the maximum size (i.e., as B approaches n). For all three group sizes (n=20, 50 and 100), the proportion of R at internal equilibrium gets lower than 10% when the efficiency of cooperation (B/n) is larger than about 0.9 (i.e., B ≈ 18, 45 and 90 for n = 25, 50 and 100, respectively) which indicates that the invasion of R becomes realistic only in a narrow range of B that is close to B's maximum value. On the other hand, the basin of attraction for T is less sensitive to the size of B as group size increases.

Synergic effect between continuous reciprocity and assortative interaction

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

In dyadic cooperation, a strong synergic effect between reciprocity and assortative interaction is known to exist [START_REF] Axelrod | The evolution of cooperation[END_REF]; rare reciprocators can invade a population of unconditional defectors much more easily if reciprocators can preferentially interact with other reciprocators and thus decrease the risk of being exploited by defectors. On the other hand, a model of discrete group reciprocity (Boyd & Richerson, 1988) found that assortative group formation does not facilitate the invasion of rare reciprocators very much.

Figure 2 shows the threshold number of interactions that must be exceeded in order for a rare reciprocator strategy, T or R, to invade a population of unconditional defectors given a specific group size and level of assortative interaction (r) (see Appendix B). The figure illustrates that both T and R can invade when group size is very small (n = 3 or 5) and when groups continue to interact for a very long time. The threshold value for number of interactions is slightly smaller for R than T, implying that R can more readily evolve, although this difference is minor. As group size gets larger, it quickly gets much tougher or even impossible for both T and R to invade (see also Boyd & Richerson, 1988).

We also investigated the synergic effect of reciprocity and assortative interaction by considering the size of the basin of attraction for reciprocators. We did find a strong synergic effect, but only for R: with assortative interaction, the basin of attraction for R remains large even when cooperation is relatively inefficient. Figure 1b shows the minimum proportion of discrete reciprocators, T, necessary for T to increase its share when assortative interaction is fairly low (r = 1/16) and very low (r = 1/64), and suggests

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

that assortation does not help T very much: in both cases, T cannot prevent the invasion of D if its proportion in a population is smaller than around 85-95%. In contrast, a small degree of assortative interaction drastically decreases the value of efficiency of cooperation (B/n) that is necessary for R to evolve. When r = 1/16, a small proportion of R can proliferate even when cooperation is relatively inefficient (B/n = 0.13). Again, similar results are obtained even when group size is increased from 100 to 500, in increments of 100. Figure 1d shows how the basin of attraction for T and R changes in response to the size of benefit of cooperation (B) and the group size (n) when r = 1/16. As when there is no assortative interaction existing (i.e., Figure 1c), the basin of attraction for R changes almost linearly as B increases to its maximum value (i.e., n/{(n-1)r+1}; if B is larger than this point, even a rare unconditional cooperator can increase its share in the population of unconditional defectors).

Where does the synergic effect between assortative interaction and continuous reciprocity come from? First, T gains the benefits of mutual cooperation only when all of its co-members are also T; the introduction of just one D co-member induces the collapse of cooperation. The introduction of one or more D also results in the collapse of cooperation when D and R coexist in the same group. However, R gains some benefit from mutual cooperation because of its gradual decrease of cooperation. As the number of D co-members increases, R stops cooperating more quickly and thus minimizes the extent of its own exploitation. R is thus able to reduce cooperation when interacting with too many Ds, and to maintain a moderately high level of cooperation with fewer Ds. As a result, payoffs to D and R marginally increase as the number of reciprocators increases, with C and T, but only a small minority of this area when it coexists with C and R.

Again, similar results were obtained when group size was increased incrementally to 500.

Extending the strategy space to allow below-and above-mean continuous reciprocity

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

We next investigated the stability of the continuous mean matching strategy against other variants of continuous reciprocity, such as those who give above or below the mean. Following [START_REF] Le | Evolutionary dynamics of the continuous iterated Prisoner's dilemma[END_REF], we modelled the strategy space of continuous reciprocal strategies using a parameter, ρ, which ranges from 0 to 2. When ρ = 1, it is identical with the R mean matching strategy. As ρ gets larger, a player gets more generous than R and contributes more than the mean contribution. When ρ = 2, it is identical with unconditional cooperation. As ρ gets smaller than 1 it gets less generous, and when ρ = 0 it is identical with unconditional defection. Contribution at the round t (≠ 1) is ρu t-1 when 0 ≤ ρ ≤ 1 and u t-1 + (ρ-1)(1-u t-1 ) when 1 < ρ ≤ 2 where u t-1 is the mean contribution made by the other players at round t-1. For simplicity, we assumed that all strategies other than D (ρ = 0) contribute fully in the first round. As we could not find a closed form expression of the payoffs, we conducted a numerical simulation by fixing a number of interactions instead of using a continuation probability of interactions.

Numerical analysis

We first assumed that the population consists of a single strategy, ρ k , and checked if a rare single strategy, ρ r , could invade this population. Figure 5 is the payoff map of various combinations of ρ k and ρ r that shows local dynamics of the system. This system has one non-cooperative attractor, ρ k = 0: when the common strategy ρ k is below the threshold value q, a more generous strategy (ρ r > ρ k ) cannot invade the population.

When the common strategy is more generous than this threshold line q but less generous than the mean matching strategy (i.e., q < ρ k < 1), a more generous strategy (ρ r > ρ k ) can invade the population. Remember that strategies 1 ≤ ρ k < q' are neutrally stable against population becomes too generous (ρ k > q'), any strategy that is less generous than the mean matching strategy (ρ k < 1) can invade the population. We found that the basin of attraction for ρ k = 0 is rather small as long as efficiency of cooperation remains high.

Agent-based simulations

We then conducted agent-based simulations to further investigate global dynamics of the system and the stability of cooperative states in a heterogeneous population where agents' strategy ρ can take any value ranging from 0 to 2. At the first generation, the entire population was occupied by ρ = 1. At the end of each generation, each individual produced a number of offspring that was proportional to the payoff that the individual received in that generation, and Gaussian error mutation ~N(0,SD) was added to each individual's strategy. Figure 6 shows the results of simulations of 5,000 agents that were randomly subdivided into 50 groups of 100 agents at the beginning of each generation and that played 20 rounds ( w = 0.95) of the repeated n-person prisoner's dilemma game for 5,000 generations under several mutation sizes (SD of Gaussian mutation: 0.002, 0.02, 0.05 and 0.07).

As the mutation size increased, the system fluctuated more frequently. In general, however, the system exhibited cooperative polymorphism: the average of ρ fluctuated between cooperative (ρ > 1) and less cooperative values (ρ 0.2~1) across generations but never arrived or stayed at a non-cooperative equilibrium (ρ = 0).

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

Discussion

This model investigated reciprocity under rather specific conditions, where new groups were formed and dismissed in each generation, in the absence of errors of perception and implementation. We intentionally chose those settings so that our model would be directly comparable to a very influential study [START_REF] Boyd | The evolution of reciprocity in sizable groups[END_REF]) that is often cited for suggesting that (discrete) reciprocity is unlikely to have played a very important role in the evolution of group cooperation (e.g. [START_REF] Boyd | The evolution of altruistic punishment[END_REF][START_REF] Fehr | Don't lose your reputation[END_REF][START_REF] Fehr | The nature of human altruism[END_REF][START_REF] Gächter | Reciprocity, culture and human cooperation: previous insights and a new cross-cultural experiment[END_REF][START_REF] Hagen | Game theory and human evolution: A critique of some recent interpretations of experimental games[END_REF][START_REF] Hauert | Volunteering as red queen mechanism for cooperation in public goods games[END_REF][START_REF] Henrich | Cultural group selection, coevolutionary processes and large-scale cooperation[END_REF][START_REF] Henrich | Cooperation, punishment, and the evolution of human institutions[END_REF][START_REF] Kurzban | Experiments investigating cooperative types in humans: A complement to evolutionary theory and simulations[END_REF][START_REF] Sigmund | A tale of two selves[END_REF][START_REF] Suzuki | Reputation and the evolution of cooperation in sizable groups[END_REF]. Our results suggest that models which represent reciprocity as discrete rather than continuous may underestimate the adaptiveness of reciprocity in groups.

Although some real world group cooperation surely involves all-or-nothing discrete decisions, continuous cooperation has been observed in numerous species (Bshary & Bronstein, 2004). The disregard of reciprocity's role in the evolution of collective action is surprising, given that (continuous) reciprocity is the most commonlyobserved strategy pursued by subjects in n-person cooperation experiments [START_REF] Croson | Theories of commitment, altruism and reciprocity: evidence from linear public goods games[END_REF][START_REF] Croson | Reciprocity, matching and conditional cooperation in two public goods games[END_REF][START_REF] Fischbacher | Are people conditionally cooperative? Evidence from a public good experiment[END_REF][START_REF] Kurzban | Experiments investigating cooperative types in humans: A complement to evolutionary theory and simulations[END_REF][START_REF] Kurzban | Incremental commitment and reciprocity in a real time public goods game[END_REF]Yamagishi, 1986;Yamagisi & Sato, 1986). In this study, we pointed out the difference between the discrete reciprocity that has prevailed in past modelling efforts, and the continuous reciprocity that has been observed in these experiments. While the In the absence of assortative interaction, when the efficiency of cooperation (B/n) approached around 0.9 (i.e., B = 18, 45 and 90 for groups of 20, 50 and 100 individuals, respectively), the proportion of continuous reciprocators necessary for invading a population of defectors dropped to below 10% (Figure 1c). When groups were assortatively formed, the value of B required for continuous reciprocity to invade the population fell significantly (Figure 1d). For instance, when degree of assortative group formation was fairly weak (r = 1/16), the proportion of continuous reciprocators necessary for invading the population dropped below 10% when the values of B were around 8, 11, and 12.5 for groups of size 20, 50 and 100, respectively. In both cases with and without assortative interaction, invasion of continuous reciprocity seems to be realistic only in narrow range where the size of the benefit from cooperation comes close to its maximum bound. When the value of the benefit from cooperation was sufficiently large, continuous reciprocity was stable even in a heterogeneous population where strategy space was extended to allow for below-and above-mean reciprocity.

It is important to remember that the value of B indicates the extent to which one unit of contribution by one member benefits the member's entire group. For example, if B

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

= 20, this means that when one member expends one unit of effort, it produces 20 units of public good for the group. Cooperation may under some extreme circumstances lead to extraordinarily efficient public good provisioning-for example, contributing to the discovery of a new food source during a time of starvation, or helping to build a higher wall in order to defend a village against a large-scale attack. Further, with some public goods, for example public buildings or roads, B tends to increase as n increases, because the total benefit produced often increases with the number of people who utilize the good.

However, it is not clear that such efficient public goods production would have been a

realistic aspect of the societies in which human adaptations for reciprocity evolved.

The model presented suggests that continuous reciprocity can evolve in large groups only under a narrow range of circumstances involving unusually high productive efficiency. In this sense, the puzzle of the evolution of large-scale cooperation still remains to be solved. On the other hand, the overall effects of replacing discrete with continuous reciprocity were significant. Future models of the evolution of cooperation should incorporate real-world behavioural regularities that are observed in the context of both cooperation and punishment (e.g., Carpenter 2007), while still keeping models as simple as possible.

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

i Throughout this paper, we manipulate the value of B/n as a model parameter instead of independently varying the values of B and n. B/n is frequently used as an index of efficiency of cooperation; in order for the game to be a public goods dilemma, its value needs to lie between 0 and c (= 1 in our model). Differences between the two reciprocal strategies were exhibited mainly when this index's value was close to its upper limit.

Note that the value of B, the amount of benefit produced from one unit of contribution, varies with the increase of the group size, n, even when the efficiency of cooperation index remains constant. population is evolutionarily stable, and an all-R or all-T population is neutrally stable, as C can receive the same payoff as R or T when D is absent from the population. In both 4a and 4b, the line QC separates two regions. Within DQC, the system evolves towards a non-cooperative equilibrium (point D). Within QCR or QCT, the system evolves towards a point on the line CR or CT (i.e., a mix of strategies C and R or T). On the line QR or QT, the system evolves towards an all-reciprocator equilibrium (point R or point T).

When R is present in the population (4a), the basin of attraction of cooperative strategies C and R is much larger than that of D; however when T is present (4b), the basin of attraction of cooperative strategies C and T is much smaller than that of D. For both 4a and 4b, the parameter values are n = 100, B/n = 0.13, w = 0.95, and r = 1/16. and maximum number of rounds = 10 ( w = 0.9). The parameters q and q' are strongly influenced by the efficiency of cooperation. For instance, when B/n is decreased to 0.65, q increases to 0.36~0.41 and q' decreases to 1.61.

p(x) n 1 x n n r (1 r)p r r r p x (1 r)(1 p) (1 ) n 1 x .
In this equation, r is a parameter, ranging from zero to one, determining the degree of assortative interaction; groups are randomly formed when r = 0. Here, r plays the same role as Hamilton's kinship coefficient [START_REF] Hamilton | The genetical evolution of social behavior[END_REF] and even unconditional cooperation can evolve when {(n-1)r+1}B/n > c, which is called the inclusive fitness effect [START_REF] Hamilton | Innate social aptitudes of man: An approach from evolutionary genetics[END_REF]. When r > 0, of particular interest is the situation in which this equation does not hold, and unconditional cooperation can therefore not evolve.

We analyzed the influence of the assortative interactions in two different ways:

invasion analysis and basin of attraction analysis. First, we investigated a situation where an invading strategy is very rare. Assume that p is the proportion of a rare invading strategy in a population. As p approaches 0, the above mentioned probability p(x) where V(R| m x 1) is (x 1)B/n 1 when x < n-1 and (x 1)B /n 1 1 w when x = n-1.

Second, we investigated a situation in which the proportion p of an invading strategy takes any value other than zero. Figure 1b in the paper shows the results of this analysis.

C. When an unconditional cooperator (C) exists:

When the three strategies R, C and D coexist in a group, the contribution made by any R player in round i is z i v i 1 s v j 1 v i 1 s(1 v i 1 ) /(1 r) , where p(x,y) is the probability of finding x R and y = n-1-x C in a group given that the proportion of R and C in a population is p m and p k , respectively. In general, the probability that a focal player of a strategy s 1 will find himself in a group of j s 1 players, k s 2 players and n-x-1 s 3 players is p( j,k) j,k,n 1 j ! r (1 r) p 1 r j (1 r) p 2 (1 k (1 r)(1 p 1 p 2 ) ( ( 1)

n 1 j k
given that the proportion of s 1 and s 2 in a population are p 1 and p 2 and that the degree of assortative interaction is r.

  T a cannot invade a population when D is common. When T a is common, the expectedpayoff to T a is V(T a |n) = (B -1) / (1 -w). The expected payoff to a rare D is V(D|n-1) = B(n 1) / n1 w when a < n -1. As V(T a |n) < V(D|n-1) when n > 2, T a cannot prevent the intrusion of D. On the other hand, when a = n -1, the expected payoff to D is V(D |n-1)

  1) / n and 1-w(n-2)/(n-1) is always smaller than 1, R can avoid an invasion of D in a narrower range of conditions than can T.

  to D and T mostly linearly increase. This non-linearity advantages R, especially when assortative interaction exists (see Figure3).2.5. The third strategy: unconditional cooperatorsSo far we have investigated interactions between only two strategies (D and either R or T) at one time. Do results change when a third strategy -unconditional cooperation (C) -is added? When D is completely absent, C receives the same benefit as T or R, and C can drift into the populations composed of reciprocators (T or R). However, C is easy prey to D, and after drifting onto a population, C may even help D to beat a reciprocal strategy. In order to examine the influence of unconditional cooperators, we extended the model by adding C and examined the evolutionary dynamics of the system. Appendix C discusses the payoffs and Figure 4 illustrates the dynamics of three strategies, C, D and either T (Figure 4a) or R (Figure 4b), when n = 100. In both figures, any points on the line QC are unstable internal equilibria and the area above the line QC (including the point C but excluding the point Q) is the basin of attraction for D. Any points on the line TC or RC other than the point C are neutrally stable fixed points, and the area below the line QC is the basin of attraction for a mixture of cooperative strategies. We found that D's basin of attraction composes a large majority of the triangular area when it coexists

  ρ r ≥ 1, as neither perception nor implementation error exists in the model and both the common and invading strategies receive the same fitness value. If a

  and continuous reciprocity may at first glance seem trivial, it affected model results significantly: compared to discrete reciprocity, the continuous mean matching strategy greatly enlarged the basin of attraction where reciprocity can evolve. However, this advantage of continuous reciprocity was heavily constrained by the size of the benefit from cooperation; the basin of attraction for continuous reciprocators increased almost linearly as the benefit from cooperation (B) increased.

  and by a grant to the Santa Fe Institute from the James S. McDonnell Foundation-21 st Century Collaborative Award Studying Complex Systems.

Figure 1 :

 1 FIGURE CAPTIONS

Figure 2 :

 2 Figure 2: Minimum number of interactions that must be exceeded for rare

Figure 3 :

 3 Figure 3: Payoff to defectors (D) and discrete (T) and continuous (R) reciprocity

Figure 4 :

 4 Figure 4: Evolutionary dynamics of unconditional cooperators (C), unconditional

Figure 5 :

 5 Figure 5: Comparisons of payoffs of common strategy ρ k and rare invading strategy ρ r . Vertical axis represents a common strategy in a population and horizontal axis

Figure 6 :

 6 Figure 6: Evolution of generosity (ρ) in the agent-based computer simulations. For

  1 r) n 1 x . When D is common in the population, the fitness of a rare invading R is obtained by recalculating W(R) in the section (i) after replacing p(x) with the new value derived in this section. As the average fitness of the common strategy D is equal to the baseline population fitness W 0 , R can invade the population when W(R) > W 0 . The threshold value of the number of interactions that must be exceeded by R in Figure 2 in the paper was derived by numerically solving this equation given that the other parameters were fixed. The threshold value for T was calculated in the same manner given that W (T) W 0 p(x) V(T | m x 1)

  s = k/(n-1), and m and k are the number of R and C players in the group, respectively. Thus, the average payoff to R in a group of m R and k C is

  average payoffs to C and D in the same group are V(C|m,k) = w i 1 B z i m /n Hence, the average payoffs to R, C, and D in a population are W (R) W 0 p(x, y) V (R | m x 1,k y)

Acknowledgements

Masanori Takezawa was financially supported by the Max Planck Institute for Human Development and the Tilburg Institute for Behavioral Economic Research. Michael Price was funded by the Indiana University Workshop in Political Theory and Policy Analysis,

When there are m continuous reciprocators and n -m defectors in a group, the contribution made by any R in round i is z i = v i-1 where v = (m-1)/(n-1). Thus in this group, in round i the payoff to R is B z i m / n z i and the payoff to

strictly larger than V(D|n-1) when the following condition is satisfied:

The expected payoff to R in a group of m R players and n -m D players is V(R|m) =

w . Assuming random group formation, the average payoff to R at the population level is

, where W 0 is the baseline fitness, p(x) is the probability for a focal player to find himself

, when n>m+k, and

w i , respectively. The average payoffs to the three