EVOLUTION OF UNCONTROLLED PROLIFERATION AND THE ANGIOGENIC SWITCH IN CANCER

The major goal of evolutionary oncology is to explain how malignant traits evolve to become cancer “hallmarks.” One such hallmark—the angiogenic switch—is difficult to explain for the same reason altruism is difficult to explain. An angiogenic clone is vulnerable to “cheater” lineages that shunt energy from angiogenesis to proliferation, allowing the cheater to outcompete cooperative phenotypes in the environment built by the cooperators. Here we show that cellor clone-level selection is sufficient to explain the angiogenic switch, but not because of direct selection on angiogenesis factor secretion— angiogenic potential evolves only as a pleiotropic afterthought. We study a multiscale mathematical model that includes an energy management system in an evolving angiogenic tumor. The energy management model makes the counterintuitive prediction that ATP concentration in resting cells increases with increasing ATP hydrolysis, as seen in other theoretical and empirical studies. As a result, increasing ATP hydrolysis for angiogenesis can increase proliferative potential, which is the trait directly under selection. Intriguingly, this energy dynamic allows an evolutionary stable angiogenesis strategy, but this strategy is an evolutionary repeller, leading to runaway selection for extreme vascular hypoor hyperplasia. The former case yields a tumor-on-atumor, or hypertumor, as predicted in other studies, and the latter case may explain vascular hyperplasia evident in certain tumor types.

1. Introduction.The notion that cancer is essentially an evolutionary phenomenon has a deep history in oncology (see [18,19,36,39,43] for reviews).Evolutionary forces shape the disease on at least three levels.At the individual organism level, natural selection favors tumor suppressor mechanisms that keep cancer at bay at least through the reproductive years.At the clinical level, natural selection explains treatment failure, as first recognized by Law in 1952 [35].But most important to our understanding of cancer pathobiology is evolution at the individual cell level.Here, genetic and genomic alterations generate the raw material of natural selection-alternative genotypes that produce differential reproductive success among neighboring cells.Mutations that confer a proliferative phenotype tend to expand in the population.Within proliferative lineages further mutations occur, typically at a highly accelerated rate due to genetic instability.Selection then sifts among new mutants, favoring those with a proliferative or survival advantage.The result is neoplasia.In the proper environment, selective pressures can drive the neoplasm to malignancy.Such evolutionary trajectories have been mapped using chromosomal/genomic markers from the 1950s [12,24,37,47] to the present (e.g., [34,40,45,49,52]).
Phenotypic traits predicted to be differentially favored in the evolution of malignancy have been laid out in a classic review paper by Hanahan and Weinberg [21], which has been updated recently [22].The evolutionary benefits of most of these traits are essentially obvious-for example, a clone inheriting gene amplification of a receptor tyrosine kinase will acquire a greater proliferative potential.Less obvious are the traits' costs.Using the same example, a proliferative lineage must allocate significant resources to protein production to support proliferation.Such a drain, coupled with deranged energetic controls, could decrease resources for maintenance, like the Na + /K + ATPase pump, making highly proliferative cells less robust.
Of all the traits identified by Hanahan and Weinberg, the angiogenic switch is one of the more perplexing from an evolutionary viewpoint.Superficially, angiogenesis seems obviously beneficial.Tumor growth in most cases requires neoangiogenesis [20].Tumors failing to generate an angiogenic signal tend to remain tiny, insignificant proliferative foci.But the evolutionary picture is less clear.Selection does not act on the tumor per se.It acts on clones, individual cells or their genomes.Tumorlevel characteristics, like vascularization, arise as a consequence of this selection; they are not traits.Analyzed using the proper unit of selection, the angiogenic phenotype is difficult to explain for the same reason altruism is difficult to explain.
The evolutionary problem, however, is why angiogenic tumors are not invaded by mutant lineages with their angiogenesis signaling apparatus knocked out.By avoiding secretion, such a mutant would seemingly gain an energetic advantage over its angiogenic competitors.In this view, one sees angiogenesis as an altruistic trait susceptible to a cheater phenotype.Such cheater lineages, or "hypertumors," have been predicted by mathematical models of tumor angiogenesis [42,44].In these models, selection always favored more proliferative over less proliferative clones.Angiogenesis, on the other hand, was always evolutionarily neutral.However, cell energetics were not included in these models, so the could not properly evaluate metabolic tradeoffs between cell proliferation and angiogenesis signaling.Physiological capabilities of cells are limited by their energy states.This limitation largely defines the costs associated with the malignant traits favored during carcinogenesis.
Here we extend the previous models by introducing intracellular adenylate metabolism to more adequately address the evolution of angiogenesis and proliferative potential within tumors under the hypotehsis that costs and benefits of these "life history strategies" are a consequence of their effect on cell energetics.We introduce the modeling details in section 2 and parameterizations in section 3.In section 4 we present results of the model at intracellular (section 4.2), tissue (section 4.3) and evolutionary (section 4.4) scales.The model predicts that the angiogenic switch can evolve from selection at the level of cell clones, but not directly because of its ability to grow new blood vessels.Instead, angiogenesis arises either as an evolutionary rider on the back of another trait favored by selection, including possibly proliferative potential, or via runaway selection for ever-increasing angiogenesis signaling.In both cases, the angiogenic switch is an indirect consequence of a complex interaction between natural selection and the cell's management of energy charge.
2. The Models.Since metabolic effectors are intracellular while tumor growth is a tissue-level phenomenon, we take a multiscale approach and model metabolism on a spatial scale of micrometers and seconds to minutes and tumor dynamics in terms of millimeters and hours.These scales are connected by the chemical potential energy available to tumor cells.At the tissue level, tumor growth and angiogenesis factor secretion depends on available chemical energy, primarily in the form of ATP.In turn, the available energy is determined by oxygen and nutrient supply, themselves functions of tumor vascularization and under the control of angiogenesis factors secreted by tumor cells.
2.1.Cell Energetics Model.Physiologists typically measure a cell's "energy" status by a sort of weighted average of "high-energy" phosphates-β and γ phosphoryl groups-per adenylate molecule.(For convenience I abuse the term "adenylate" slightly by restricting its meaning to adenosine 5 mono-, di-and triphosphate-AMP, ADP and ATP, respectively.)This measure, called the energy charge and denoted φ, is usually defined as In cells, the enzyme adenylate kinase tightly controls the relative amounts of the various adenylate compounds via the reaction,
Forward and backward reaction rates are nearly equal in most cells and large relative to other metabolic processes involving adenylate [23], so typically, All other important reactions involving adenylate can be classified as either (i) synthesis de novo, (ii) irreversible destruction, (iii) conversion of one adenylate form to another; e.g., ATP → ADP by hydrolysis, or ADP → ATP by glycolysis or oxidative phosphorylation.Adenylate is constructed de novo in two main pathways: salvage of adenosine from nucleic acid breakdown, and construction of AMP from inosine monophosphate (IMP) which itself is made from phosphoribosyl pyrophosphate (PRPP) and an appropriate amino acid-glutamine, glycine or aspartate-via a complex pathway [8,9,17].No distinction among various sources of adenylate will be made in the model.
The second class of reactions-adenylate destruction-includes reactions mediated by 5 nucleotidases, referred to as AMP phosphatase by Ataullakhanov and Vitvitsky [4].For example, ecto-5 nucleotidase removes phosphate from AMP, creating the nucleoside adenosine.Also, AMP can be broken down by various enzymes in the AMP deaminase family, which convert AMP to IMP.In addition to these mechanisms, adenylate can also "disappear" as ATP is incorporated as structural elements in nucleic acids.Although these adenosine phosphates can be recovered via nucleic acid breakdown, in the model this reaction acts as a sink for ATP.
Interconversions among various forms of adenosine phosphate present a much more diverse array of reactions.In addition to the adenylate kinase reaction described above, ATP is converted to ADP by a variety of hydrolytic reactions, most notably the membrane-bound Na + /K + ATPase, which controls cell volume and accounts for up to 40% of ATP hydrolysis in resting eukaryotic cells [1, pg. 210].ATP is also converted to ADP in the adenosine kinase reaction, which serves as another source of AMP synthesis de novo.ATP can be converted to AMP via a number of reactions, most notably the following: (i) adenylation of amino acids for protein synthesis, (ii) the phosphoribosyl pyrophosphate synthetase reaction, Ribose 5-P + ATP → PRPP + AMP, and (iii) a two-step reaction in which adenylyl cyclase first converts ATP to cyclic-AMP (cAMP) and pyrophosphate, followed by the conversion of cAMP to AMP by cAMP phosphodiesterase.Finally, ATP can be regenerated from ADP by glycolysis and oxidative phosphorylation.
Our immediate goal is to model all these reactions in a realistic yet simple way.The energetics model operates on a much faster time and smaller spatial scale than the tissue-level tumor model to come.Also, intracellular ATP concentrations are about 3 mM (M = mol/L), which equates to about a billion ATP per cell [1,13].The ATP concentration in resting cells is about 10 times that of ADP and 100 times that of AMP, implying tens of millions of AMP per cell.Cells hydrolyze ATP at a rate nearly 10 8 per second [33].Therefore, in the long run, stochastic fluctuations both in time and space are assumed to have negligible effects at the level of a whole tumor.So, for simplicity we choose to model adenylate dynamics with a system of ODEs, using the pioneering work of Ataullakhanov, Vitvitsky and colleagues [2,3,4,5,38] as our starting point.
Let A i , i ∈ {1, 2, 3}, represent adenosine with i 5 phosphates.Then our adenylate model takes the following form: Relative Glycolysis Rate q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 1.Proposed model of the unscaled glycolysis rate, g, as a function of cell energy charge (curve) with erythrocyte glycolysis data from [4] superimposed (points; see text).
where a, b and c are constants that can be decomposed as follows: The rate of AMP synthesis de novo is denoted α a .The basic rate of the adenylate kinase reaction, both forward and back, is k.Parameter γ a is the basic rate of the adenosine kinase reaction (the product of a mass action constant and the concentration of adenosine, assumed to be constant).ATP is incorporated into nucleic acid or otherwise irreversibly destroyed at per-molecule rate µ, and AMP destruction by 5 nucleotidases and AMP deaminases is modeled by the function f (A 1 , A 3 ).The parameters β 1 and β 2 represent the the cell's constant demand for ATP to maintain normal physiological function.Parameter β 1 can be interpreted as the energy demanded by general biosynthesis, whereas β 2 is energy demand for maintenance physiology, especially that required to power the Na + /K + -ATPase to maintain cell volume.In addition to these purely maintenance functions, cells can allocate ATP to protein synthesis in support of either proliferation or TAF secretion.The allocation to proliferation, denoted η p , is assumed to be under genetic control and invariant over changes in local vascularization, v. On the other hand, ATP allocation to TAF secretion, η s (v), varies with local vascular density.In general, η s (v) should approach 0 as v → 0 or as v → ∞, reaching a unique maximum at some relatively small v [15,25,42].Here we assume that TAF secretion is proportional to ATP allocated for that function, and follow [42] in assuming the following form for such allocation: ) Parameter ηs is the basic commitment to TAF secretion, which like η p is assumed to be under genetic control.In this model TAF secretion peaks at v = ξ −1 , which is assumed to be constant across tumor cell strains.
The function G(φ, v) represents ATP regeneration from ADP by both glycolysis and oxidative phosphorylation.However, both for simplicity and because many malignant tumors rely on glycolysis to regenerate ATP, in this model we will assume that oxidative phosphorylation is negligible.
The rate of glycolysis appears to be a very simple function of intracellular ATP concentration, at least in human erythrocytes [4].Ataullakhanov and Vitvitsky's [4] data suggest that glycolysis rate in erythrocytes can be well modeled as a quadratic function of ATP concentration, although they observed enormous inter-individual variation.Following their lead, we combine data among individuals in their data set and calculate the best-fit quadratic function: 0.06595 + 1.62033A 3 − 0.70186A 2 3 which yields r 2 = 0.8467.However, there is good reason to believe that ATP glycolysis rate is tied to energy charge, not ATP concentration per se [4,6].Assuming that the data from Ataullakhanov and Vitvitsky varies because of variations in total adenylate among cells in their samples, one can transform the quadratic model above to a function of energy charge with a zero at φ = 1.Arbitrarily setting the peak of this function to unity yields, thorough a straightforward transformation, the curve in Fig. 1, which also shows the data from [4] similarly transformed.This graph suggests the following model of glycolysis: where g(φ) is the rate of glycolysis in normalized units of amount per time.Introducing a scaling facter s, which most naturally has units fmol/min (see section 3), converts the normalized model into Biologically, s measures the maximum glycolysis rate, which depends on the amount of reduced carbon available from the blood; therefore, s is an increasing, saturating function of vascularization, v, since the plasma concentration of glucose is finite and glycolysis rate has a biochemically imposed upper limit.Let this physiological maximum be s max .Then for the scaling factor in equation (6) we propose the following model: The sensitivity parameter of 0.1 was chosen to match the model in [42].
2.2.Tumor Growth Model.To model tumor growth at the tissue level we adapt model (3) to the model of Nagy [42].In particular, for a tumor comprising only one parenchyma cell phenotype, where x(t), y(t) and z(t) represent tumor parenchyma mass, immature vascular endothelial cell (VEC) mass and total microvessel length, respectively.Functions Φ(A 3 ) and Ψ(v, A 3 ) model per capita growth rates for parenchyma and VEC populations, respectively, γ is the rate at which endothelial cells become incorporated into functional microvessels, and δ represents the microvessel remodeling rate.We assume that parenchyma cell proliferation exhibits Michalis-Menten dynamics as a function of ATP allocated to proliferation, with maximum p and sensitivity k s .Also, as a first approximation we assume that mortality is inversely proportional to ATP concentration, with constant m.Therefore, where Ā3 is the equilibrium [ATP] from model (3).Unlike parenchyma cells, which proliferate on their own, immature vascular endothelial cells require an angiogenic signal to stimulate proliferation.The strength of that signal is assumed to be proportional to ATP allocated to producing and secreting the molecules involved, specifically η s (v)A 3 .However, the response of immature VECs to this signal is nonlinear because biological limitations bar cell proliferation from running off to an arbitrarily large rate.Here we assume that VEC response to this signal obeys a Michaelis-Menten function of angiogenic signal strength.To this we add constant VEC maturation and death rates to produce where β is the sum of VEC death and maturation rates.Biologically, α represents the maximum per-capita proliferation rate of VECs, and k v measures VEC sensitivity to the angiogenic signal.
We view Φ and Ψ as essentially placeholder functions with the correct general properties, not detailed hypotheses about proliferation and angiogenic dynamics.As in [42], the precise forms of these functions have little impact on the qualitative results as long as they are generally saturating or unimodal functions of A 3 .

Initial Conditions and Parameterization.
The notation used in these models is summarized in Tables 1 and 2. The model's complexity is reflected in the length of the parameter list.However, reasonable estimates for most can be obtained.In some cases, these estimates can be argued from biological first principles.Tumor microvessel length density (= z/x) VEC proliferation rate Others are founded on entirely empirical grounds.Nevertheless, some are still a priori unobtainable.For parameters in this last category, we attempt to define their magnitudes by identifying values that produce reasonable biological behavior.Such estimates amount to testable biological hypotheses on parameter measures that are not yet established empirically.Below we develop our estimates for all parameters and identify their sources or derivations.

3.1.
Estimates from biological measurements.Initial conditions for the adenylate model are derived directly from measures of ATP concentrations in cells.Standard textbooks-Alberts et al. [1, pg. 67] for example-tend to claim that the generic resting cell maintains an equilibrium of about 1 billion ATP molecules.This value agrees well with Du et al.'s [13] measurements in anesthetized rat brain, in which they find an ATP concentration of 3.0 mM, which equates to almost exactly 1 billion ATP per cell assuming a cell volume equivalent to a 10 µm diameter sphere, or a concentration of about 1.5 fmol/cell.Therefore we initially take A 3 (0) = 1.5 fmol (per cell).Since resting cells have about 10 times more ATP than ADP and about 100 times more ATP than AMP [23], we initially set A 1 (0) = 0.1A 2 (0) and A 2 (0) = 0.1A 3 (0).These values will be modified slightly based on model output (see section 3.3).
Next we attack the functions G(v, φ) and f (A 1 , A 3 ).Beginning with the former, imagine a cell in a abundantly vascularized tissue.Data from erythrocytes [4] suggest that the mean maximum rate of glycolysis for such a cell is approximately 2 mmol/L/hr, or about 1.7×10 −2 fmol/cell/min, again assuming a cell with a volume equal to a 10 µm diameter sphere.However, these data probably represent a lower bound for maximum glycolysis rate in generic cells, since erythrocytes are relatively metabolically inactive.Indeed, Kilburn et al. [33] found that mouse LS cells (in vitro) hydrolyze ATP at 11.7 fmol/cell/min.Assuming a resting cell maintains an ADP concentration of 0.15 fmol, a constant energy charge, φ, and has a constant nutrient supply, then Kilburn et al.'s data imply G( φ, •) = 78 min −1 .From the basic form of G(v, φ) in Fig. 1 and the observation that maximum glycolysis rates are approximately 5 times resting rates [4], then maximum per ADP glycolysis rate should be around 390 min −1 , which we take as our estimate of s max .The model for AMP destruction, f (A 1 , A 3 ), is taken directly from the work of Ataullakhanov, Vitvitsky and colleages [4,38].In essence, their models assume that activities of the two AMP-destroying enzyme classes are essentially negligible at physiological AMP concentrations in resting cells, approximately 0.015 fmol/cell.However, they suggest that the activity of 5 nucleotidases dominates at relatively low AMP concentration, whereas AMP deaminases dominate at high concentrations.In particular, Martinov et al. [38] use to model action of AMP deaminases, whereas their model of 5 nucleotidases takes the following form: where M 1 = 0. we take f (A 1 , A 3 ) to be the sum of expressions ( 11) and ( 12) with these parameter values (Fig. 2).Parameters β 1 and β 2 can also be estimated from data.If Kilburn et al.'s [33] estimate of ATP hydrolysis rate represents generic cells at rest, then (β 1 + β 2 ) Â3 = 11.7 fmol/cell/min.With Â3 = 1.5 fmol/cell, then β 1 + β 2 ≈ 8 min −1 .For simplicity, we set β 1 = β 2 = 4 as defaults.
The cost of TAF secretion is very difficult to measure, largely because the precise nature of the TAF signal varies among tissues.In most tissues, however, vascular endothelial growth factor (VEGF) and various angiopoietins play significant roles.Recently, Karoubi et al. [31] measured the cost of VEGF secretion in an assay in which eukaryotic cells were transfected with a VEGF-carrying plasmid.In one assay they measured the secretion rate of VEGF 165 to be 8 fg/cell/hr, which equates to about 3.6 × 10 −4 fmol/cell/hr.Given that VEGF 165 's molecular mass is 22.32 kDa, synthesis of this compound requires adenylation of 191 residues per molecule.So, secretion at this rate burns a minimum of 6.9 × 10 −2 fmol of ATP per cell per hour.In other assays, Karoubi et al. boosted VEGF secretion up to 40 × this value by increasing the number of plasmids per cell.How this level of synthesis relates to natural VEGF production is hard to evaluate.But certainly these calculations only account for synthesis of VEGF, not secretion.Nor does it measure synthesis and secretion of other compounds in the TAF signal.So, here we assume that TAF secretion requires no more than 0.7 fmol/cell/hr = 1.17 × 10 −2 fmol/cell/min, 10 times the basic value measured by Karoubi et al.Since max η s (v) = η s (ξ −1 ) = ηs ξ −1 e −1 is assumed to be 1.17 × 10 −2 , our estimate of ηs ≈ 0.11 min −1 , assuming ξ −1 = 0.3 (see below).

3.2.
Estimates from biological principles and hypotheses.As already noted, the adenylate kinase reaction is very fast relative to other metabolic reactions involving adenylate.This observation justifies the assumption that and also In numerical investigations we set k = 10 6 .
In resting cells we assume that AMP construction via adenosine kinase (γ a ) and ATP loss through nucleic acid construction (µ) are much smaller than β 1 + β 2 .Therefore, as a first approximation we set γ a = µ = 0.01.
To determine η p , we note that, with few very atypical exceptions-eggs of the frog Xenopus for example-mammalian cells double no faster than once every 12 hours.To reproduce itself, a eukaryotic cell requires approximately 1.6 × 10 4 fmol of ATP [33].Therefore, a cell growing at its maximum rate is burning ATP at a rate of approximately 22 fmol min −1 .If such a cell were able to maintain its ATP level at 1.5 fmol, then we expect (β 1 + β 2 + η p )1.5 ≈ 22 fmol/min, assuming µ and γ a (maintenance amino acid adenylation and the background adenosine kinase reaction rates, respectively) are minor drains on ATP.With the default values for β 1 and β 2 , this expectation implies an upper bound for η p of 7 min −1 .However, we leave open the possibility that pathology can cause this upper bound to increase significantly.

3.3.
Estimates from model fitting.Obtaining a priori values for the remaining parameters presents considerable technological challenges, but these parameters can be estimated from the model by considering a cell either at rest or in equilibrium within a healthy, but dynamic, tissue.Resting cells burn negligible ATP for both proliferation and TAF secretion (η p = η s (v) ≡ 0), and their adenylate pools are in equilibrium.With all other parameters set to their defaults, we find that "tuning" α a = 5.725 × 10 −5 maintains a resting cell's adenylate in equilibrium near their initial concentrations.To be precise, with α a and all other parameters set to their defaults, adenylate concentrations equilibrate at Â1 = 0.0250, Â2 = 0.1938 and Â3 = 1.5027 (Fig. 2).Therefore, in all subsequent numerical investigations we set the default adenylate concentrations to these values.
To obtain an estimate of p, the proliferation response parameter, we consider a cell dividing at its maximum rate.In section 3.  η p ≈ 7 min −1 .The adenylate model predicts that maintaining such a proliferation commitment incurs a drain on ATP (Fig. 3) which the adenylate model predicts will eventually equilibrate at about 1.4301 fmol.Therefore, to maintain a maximal proliferation rate of 0.06 hr −1 , equivalent to a 12 hour doubling time assuming exponential growth, requires Note that the base time scale for adenylate dynamics is min −1 while that of the tissue-level tumor model is hr −1 .The maximum VEC response to TAF, α, the VEC sensitivity parameter, k v , and the parenchyma mortality parameter, m, are tuned in the default model to produce a homeostatic tissue when there is a small commitment to proliferation-that is, η p is somewhat arbitrarily set to 1, representing proliferation required to replace natural cell death-in a healthy tissue in which v = 1.In such a tissue, there should be no net tissue growth ⇒ Φ( Ā3 ) = 0, or net vascular change ⇒ Ψ(1, Ā3 ) = 0.When η p = 1, then Ā3 = 1.9589 (obtained by solving model 3 to equilibrium under default parameters), so Φ = 0 is satisfied when m = 0.0698 fmol/hr.Tissue homeostasis otherwise arises when α = 0.1 hr −1 , k v = 0.0115 fmol/min (Fig. 4).This function nicely matches the Gammack et al. tumor growth model [15] as adapted in [42].
The remaining parameters, ξ (reciprocal of vascular density at which TAF secretion reaches it maximum), β (immature VEC death/maturation rate), γ (rate at which immature VECs construct functional microvessels) and δ (microvessel remodeling rate) take the same values as in [42]   .Red surface is the graph of the manifold on which X = 0 (adenylate kinase in equilibrium) and the translucent grey manifold is the set of points satisfying Y = constant (conservation of adenylate).Black curve is a numerical solution of the full model (equations [3], [4], [6], [7], [11] and [12]) with default parameters and no simplifying assumptions.I) Dynamics on time scale τ 1 , representing equilibration of adenylate kinase (note that this portion of the trajectory lays above the red manifold); II) dynamics on time scale τ 2 .Dot is the initial condition.Compare Fig. 6. function by replacing f 2 (A1, A3) in Equation ( 12) by f2 (A1, A3) = M2A3 A1 , since the relative error f2− f2 f2 for A 1 > 10 −3 is less than 5%.From Figure 2 we assume that f (A 1 , A 3 ) has the same order of magnitude as α a .Hence the scales in the adenylate model (3) are s 1 = 10 6 , s 2 = 1, s 3 = 10 −5 , corresponding to time scales τ 1 = 10 6 t, τ 2 = t, and τ 3 = 10 −5 t.
We expand the variables A i into these time scales, define two new variables and expand them as follows: Figure 6.Adenylate dynamics, with selected initial conditions.Red manifold is defined in Figure 5, and the two translucent manifolds satisfy A 3 = mA 2 and A 3 = m 2 A 1 , m = Ĝ/(2c − b) (see equations ( 21) and ( 22)).The heavy black trajectory is the same solution as shown in Figure 5, continued to show dynamics on time scale τ 3 , labeled III.Blue curves are various solutions with different starting concentrations of ATP.The heavy dot embedded in all three manifolds is an apparently asymptotically stable equilibrium attracting all these solutions.
The time evolution on the ultrafast time scale, τ 1 , is therefore described by the following system: , We find that on this time scale the initial amount of adenylate stays constant at Y (0), and the ultrafast adenylate kinase reaction equilibrates in finite time to X = 0 and a finite A 2,1 .These dynamics are shown in Fig. 5 (labeled "phase I") for a particular solution; the red surface and grey plane satisfy X = 0 and Y = Y (0) for that solution, respectively.
The precise point at which this ultrafast system equilibrates depends on the initial conditions; however, we are interested in the qualitative dynamics and so continue by assuming that at this equilibrium we have Y = 1 and consider two possible states for the evolution on the next time scale: i) A i,1 = 1 3 corresponding to an energy charge φ 1 = 1 2 and a glycolysis term G 1 = s, and ii) A 1,1 = A 2,1 = 0, A 3,1 = 1 corresponding to an energy charge of φ 1 = 1 and a glycolysis term G 1 = 0. On the time scale τ 2 , case (i) yields and in case (ii) we have Qualitatively, both cases are essentially the same.Notice that a + b − c = γ a − µ; biologically, this quantity represents the difference between rates of adenosine creation via the adenosine kinase reaction and removal for biosynthesis.For the remainder of the paper we make the simplifying assumption that γ a = µ.As a result, the total number of adenylates is still conserved, X 2 approaches an equilibrium exponentially fast, and A 2,2 (t), after an exponential transient, grows or decays linearly on the time scale τ 2 .For case (i) the energy charge, φ 2 (τ 2 ), grows linearly while for case (ii) it decays linearly.This linear behavior continues until the term f (A 1 , A 3 ) has the same order of magnitude as α a .From Figure 2 we can estimate that A 1 should increase from zero to O(10 −2 ) or O(10 −3 ) on time scale τ 2 = t, which would imply fast equilibration within about 0.01 minute; hence, energy balance among the three species of adenylate is restored to a very good approximation via enzymatic reactions and glycolysis within a few seconds.This equilibration on time scale τ 2 for the default model is evident in Figure ( 5) (the solution limb labeled "II").
Finally, total adenylate is regulated on time scale τ 3 .Note that Therefore, with our assumption that γ a = µ, an equilibrium for the entire system requires α a = f (A 1 , A 3 ).This balance between de novo AMP synthesis and AMP destruction is achieved with a very fine adjustment on the ultra-slow time scale τ 3 , corresponding to weeks according to the parameterization used here.This somewhat implausibly slow adjustment in total adenylate probably arises from a slight error in the parameterization of f , which in turn leads to an error in α a .Nevertheless, the qualitative behavior of the model is plausible (see Discussion), so we assume in the tissue-level model to come that the adenylate dynamics are locked in the (quasi-) equilibrium elucidated below.3) (no simplifying assumptions) onto the A 1 , A 3 plane for values of η p ∈ {0, 1, . . ., 12} with v = 1 and all other parameters set to default values.Solid blue lines are the initial conditions, which intersect at the equilibrium associated with η p = 0. Dashed blue line is the graph of equation (22) for η p = 12.Fixed points are black circles, blue arrows show progression of fixed points with increasing η p and black arrows show direction of solution trajectories.Evident are dynamics on time scales τ 2 (rapid movement from initial condition to cusp, left to right) and τ 3 (slow adjustment along the line defined by equation ( 22)).Dynamics on time scale τ 1 are irrelevant here since the initial conditions satisfy X = 0. Compare Fig. 8A.
Numerical results from the default parameterization (Figs. 3 and 6) and all others investigated (not shown) suggest that energy charge equilibrates on timescale τ 2 .From the definition of φ, since dY /dt = α a − f (A 1 , A 3 ).At a fixed point for the entire system, then, α a = f ( Â1 , Â3 ); therefore, at equilibrium.Substituting the explicit expressions for the derivatives followed by a bit of algebra shows that Ĝ where Ĝ is the equilibrium glycolysis rate.At equilibrium we also have k(A 2 2 − A 1 A 3 ) + aA 3 = 0, Returning to the explicit expression for f (A 1 , A 3 ) (equations ( 11) and ( 12)) and given α a = f (A 1 , A 3 ), Ĝ (equation ( 21)) and approximation (22), the equilibrium value for A 1 , say Â1 , can be approximated by solving Since Ĝ equilibrates so rapidly, identifying the (approximate) equilibrium values for all species of adenylate becomes a straightforward numerical problem, the solution to which is used to determine the quasi-equilibrium ATP concentration in the tissue-level model.
Although relation (23) was constructed to identify the location of the equilibrium adenylate concentrations, the argument leading to it generates insight into the dynamics on timescale τ 3 .Speaking roughly, the τ 2 dynamics transition to the τ 3 dynamics as energy charge, and therefore glycolysis, equilibrates and α a approaches f (A 1 , A 3 ).This is evident in Figures 6 and 7.In the former, the limb labeled "II" represents the τ 2 dynamics.As solutions approach the planes defined by equations ( 21) and ( 22) (translucent grey in the figure), solutions abruptly switch to the final equilibration phase on timescale τ 3 (labeled "III" in the figure) creating an apparent cusp.In fact, the solution remains smooth.Figure 7 shows a series of solutions in which only one parameter was varied (η p ).In each case solutions shoot to the right on timescale τ 2 from their initial condition (crossing of solid horizontal and vertical lines; note that here the dynamics on τ 1 are not evident as these initial conditions already satisfy the equilibrium conditions for system (16)).As they approach the plane satisfying condition (22), they transition to the creep towards equilibrium as they adjust total adenylate to balance de novo construction, α a , with destruction, f (A 1 , A 3 ).

Adenylate Dynamics.
With the basic qualitative dynamics understood, we now move on to a biological assessment of adenylate dynamics and the effects of variations in the two evolutionarily important parameters-proliferation and angiogenesis factor secretion.A cell that increases commitment to proliferation (increases η p ) may degrade, but may even bolster, its equilibrium ATP and total adenylate concentrations.Numerical investigation (Figs. 3 and 7) suggests that a resting cell (η p , ηs = 0) at normal tissue vascularization (v = 1) which modestly ramps up its proliferative machinery tends to increase its ATP reserves.This occurs in the model because the cell overcompensates its increasing ATP demands by a disproportionate decrease in total adenylate destruction.However, if demand for ATP becomes severe, starting just above the break-even commitment of η p = 1, equilibrium ATP concentrations begin to decline with increasing η p (Fig. 8A), reflecting  1, with v = 1 (normal tissue vascularization).the physiological stress of proliferation that can no longer be accommodated by manipulations of total adenylate.Increasing η p to 14 min −1 , approximately twice our estimated physiological maximum, causes total adenylate and ATP to decline to about 65% of normal (Fig. 8A).These observations are consistent with the model of Ataullakhanov and Vitvitsky and real cells [4].
Energy charge, on the other hand, is reasonably well buffered from changes in ATP demanded for proliferation.In a physiologically normal cell, this model predicts that energy charge will not drop beyond 92% of its normal value when η p rises to 14 min −1 (Fig. 8B).This result is largely explained by the homeostatic mechanism represented by f (A 1 , A 3 ), as suggested by Ataullakhanov and Vitvitski [4].
A cell's energy status has a similar, somewhat counterintuitive relationship with tissue vascularization in this model.In particular, as vascular density increases, ATP concentration may increase or decrease.Since blood supplies nutrientsmainly glucose and, for aerobic cells, O 2 -required for immediate regeneration of ATP, one might predict that equilibrium ATP concentration is a monotonically increasing function of vascularization, v.This model contradicts that suggestion; in highly vascularized tissues, ATP concentration may decline with increasing v (Fig. 9).As vascular density increases beyond the normal range, the model predicts that total adenylate begins to decline, decreasing the ATP concentration.Whether this prediction represents actual biological behavior or is simply model artefact is an open question that can be addressed experimentally.

4.3.
Tumor dynamics.This model predicts that cells suffer diminishing returns as they ramp up ATP hydrolysis to support proliferation.The tumor growth rate as a function of proliferative commitment (Φ(A 3 ) viewed as a function of η p ) is concave, with the mode at about η p = 3.5 (Fig. 10C).At this point the rate at which proliferation increases with ATP "spent" on proliferation is matched by rate at which necrosis increases due to ATP deficiency.As more and more ATP is drained for proliferation, necrosis outpaces reproduction, so growth rate declines.However, vascularization starts to increase as blood vessel growth can more easily keep pace with tumor growth (Fig. 10A).When η p exceeds about 10.5, ATP reserves become so degraded that the tumor begins to regress (Fig. 10B and C). Figure 10."Equilibrium" tumor properties as a function of tumor cell proliferation commitment (η p ), in a tumor with phenotypically identical cells; that is, all cells use the same η p .Although in general the tumor size is changing, these are equilibrium values in the sense that vascularization and [ATP] have stabilized.Panels represent stable (A) tissue vascularization, (B) mean intracellular ATP concentration, and (C) tumor growth rate.Left-hand limit is a bifurcation point between situations when a tumor with stable vasculature is possible (to the right) and not (to the left; see Fig. 11).Note that panel (B) differs from Fig. 8A; that latter assumes a constant v = 1.
Since equilibrium ATP concentration Ā3 is ultimately a function of vascularization v, both Φ and Ψ in model ( 8) can be viewed as functions of v. Following [42], one can therefore define w = y/z and differentiate it and v to obtain the system, as an alternative to model (8).A fixed point in model ( 24) represents a particular solution of model (8) in which the tumor grows or decays exponentially with an invariant density of microvessels and immature vascular endothelial cells.The biologically important equilibria of model ( 24) require v, w > 0, and therefore at the equilibrium [42].Since δ, v ≥ 0, no such fixed point can exist for an inviable tumor (v at which Φ(v) < 0).As shown in [42], model (24) can allow a saddle-node bifurcation if Φ and Ψ cross as variations in parameters move the two functions relative to each other (Fig. 11).Biologically this bifurcation matters only if relations (25) are met, so assume this is true.The saddle, at say (v, ŵ), has the property that Φ (v) < Ψ (v), and the (locally asymptotically stable) node has the inequality reversed.If lim inf v→∞ Φ(v) > 0, and lim sup v→∞ then by Proposition 1 in [42], the "right-hand" fixed point, i.e., that associated with the largest v, is always a node.In the model considered here, conditions ( 26) and ( 27) are not guaranteed (see e.g., Fig. 11), leading to an interesting biological possibility.If the bifurcation arises as in Figure 11, then at the the largest v satisfying Φ(v) = Ψ(v) we observe Φ < Ψ and Φ(v) < Ψ(v) for all v > v. Perturbations from v in the direction of increasing v will therefore lead to runaway (ever-increasing) vascular hyperplasia as blood vessels grow faster than the tumor for ever.We term this dynamic hypervascularization to distinguish it from another, entirely different mechanism of hyperplasia predicted by this model (see section 4.4).
Dynamic hypervascularization would be characterized by a "weedy" growth of blood vessels that the model predicts will choke the tumor to death.Typically, ATP concentration must remain finite.(Biologically this can be taken as axiomatic.)Therefore, if we let since lim v→∞ η s (v) = 0 and by assumption β > 0. Since v is destined to remain greater than v, and Φ(v) < Ψ(v) for all v > v, we have Φ(v) < 0; the tumor must regress.
4.4.Evolutionary Dynamics.We begin our analysis of the evolutionary dynamics by considering only 1-dimensional evolutionary strategies.In particular, only a single model parameter is, for now, assumed to be affected by natural selection.Let s ∈ R represent this parameter, and define S = {s; 0 ≤ s < ∞} as the strategy space.Our first goal is to identify evolutionary stable (or unbeatable [16]) strategies (ESSs) within this strategy space.As originally defined by Maynard Smith and Price in the 1970s [27,26], an ESS is any evolutionary strategy that is invulnerable to invasion by any other strategy that arises as a rare mutant in a population in which all other individuals are using the ESS.
Here we employ standard generalizations of Maynard Smith and Price's concept as developed in (phenotypic) adaptive dynamics theory (see Geritz et al. [16] for a review).First we extend the model to include two strains: a resident with mass x r , and a rare mutant with mass x m .Following [42], the 2-strain model takes the following form: where A 3i and Ψ i are the ATP concentration and angiogenesis signal output in cells of strain i ∈ {r, m}, respectively.
Following Geritz et al. [16], let s r , s m ∈ S be the strategies employed by residents and mutants, respectively.In practice, s represents a parameter in model ( 29) that takes different values in the different strains.We view s m as a strategy used by a vanishingly rare mutant in a resident population of cells otherwise using strategy s r .For now we ignore cases in which more than one mutant strain competes with the resident.We also assume that the mutant challenger arises in a tumor growing along an attracting solution for the resident population; equivalently, we assume the mutant arises in a population already sitting on one of the locally stable equilibria of model (24) representing a monomorphic tumor composed of the resident strain.
Let f (s m , s r ) be the fitness, sensu Metz et al. [41], of a mutant using strategy s m in a resident population using s r .(Maynard Smith and Price [27] and Geritz et al. [16] use a notation of the form f sr (s m ).)Given the results in section 4.3, the fitness of the mutant immediately follows from Proposition 2 in [42]; specifically, f (s m , s r ) = Φ(A 3m ; s m ) − Φ(A 3r ; s r ), (30) assuming the two strains have otherwise identical parameters.Evolution of proliferation potential.Set s = η p ; that is, consider the ATP commitment to proliferation as the parameter under selection.The fitness gradientthe derivative of the fitness in equation (30) in the direction of the mutant strategy at the point where mutant and resident strategies coincide (see [16] for technical definitions and theory)-is plotted across potential resident strategies for default parameters in Figure 12.The singular strategy (resident strategy for which the fitness gradient is zero) is marked by a vertical dashed line and is evidently locally convergence-stable as defined by [16] (roughly, it is an evolutionarily attracting strategy).The pairwise invasibility plot (PIP)-a plot of the zero contour of the fitness function (30)-for the same default tumor (Fig. 13) shows that the singular strategy is an ESS; that is, no mutant strategy in the physiologically relevant strategy set can invade a population of residents employing the singular strategy.Therefore, the singular strategy is also an evolutionary endpoint [16].So, in this model, tumors characterized by the default parameters will evolve to an ATP proliferation investment η p ≈ 3.67 if mutations cause only small changes in η p [16].4.4.2.Evolution of angiogenesis potential.In contrast to the relatively simple evolutionary dynamics of proliferative potential, angiogenesis presents a much less straightforward evolutionary picture.Interestingly, this model predicts the possibility that an angiogenic phenotype can evolve based on indirect benefits at the individual level of selection, but the result is typically runaway selection causing either vascular hypo-or hyperplasia.In either case, the model predicts the eventual production of a hypertumor.We begin by studying a tumor with "normal" cell turnover; that is, we fix η p = 1 and study evolution of ATP earmarked for TAF secretion (s = ηs ).In Figure 14 we see a singular strategy at an ATP-TAF investment rate of ηs ≈ 0.0524, well within the region of tumor viability.The PIP (Fig. 15) shows that this strategy is an ESS.However, it is not locally convergence-stable.Apparently for any resident strategy greater than (less than, respectively) the singular strategy, there exist mutant strategies with a higher (respectively, lower) TAF investment that can invade the resident population (Fig. 14), and the fitness gradient as a function of resident strategy has a positive slope at the singular strategy (in fact, throughout the range of tumor viability; Fig. 14).Therefore, the ESS is an evolutionary repeller [16].Tumors with initial investment in TAF secretion to the right of the ESS will experience runaway selection always favoring increasing TAF secretion, eventually pushing the tumor into the right-hand inviability region.The result is evolutionary suicide [48] leading to a hypertumor [42].This picture is mirrored on the other side of the singular strategy, with runaway selection eventually causing evolutionary suicide this time by favoring decreasing TAF secretion.So here we have a single mechanism, selection for TAF secretion, generating hypertumors with apparently two distinct natural histories-one would observe either hyper-or hypoplasia of the vascular net before tumor collapse.The former case contrasts with dynamic hypervascularization introduced in section 4.3.Here evolutionary forces generate the phenotype, and so we refer to this situation as evolutionary hypervascularization.
Variations in ATP investment for proliferation (η p ) changes the location of the TAF secretion ( ηs ) ESS, but it remains an evolutionary repeller.As η p increases, both boundaries for the zone of viability rapidly retreat from each other, and the graph of the fitness gradient declines while retaining the same basic shape, at least for moderate increases in η p (Fig. 16).The ESS accelerates to the right (Fig. 16) with no qualitative changes in the PIP (data not shown).If η p has reached its evolutionary endpoint at default (= 3.67; Figs. 12 and 13), the ESS for ηs has moved well beyond what is physiologically reasonable ( ηs > 90) but remains a repeller.Therefore, if a tumor has reached its evolutionary endpoint for η p before ηs has undergone any major evolutionary change, we can expect that ηs will be to the left of its ESS and the fitness gradient at this point will be negative.As a result, we would expect the evolutionary destruction of the tumor's vascular infrastructure as selection continually favors mutants with less investment in TAF secretion (Fig. 17).This situation represents the original hypertumor mechanism predicted in [42]. 5. Discussion.In this paper we study the evolution of proliferative and angiogenic potential in tumors by extending a previous evolutionary model of angiogenic tumor growth [42].The extension adds evolutionary costs in the form of energy requirements associated with particular proliferation and angiogenesis strategies.The result is a multiscale model, with cell energetics modeled on space and time scales of micrometers and seconds to minutes, respectively, and tumor dynamics modeled at millimeters and hours.The energetics model is adapted from the work of Ataullakhanov, Vitvitsky and colleagues [2,3,4,5,38].
As noted in the introduction, we focus on angiogenesis because, of the various hallmarks of cancer [21,22], angiogenic potential relies on cooperative production of the angiogenic signal.As with other altruistic traits, angiogenic "behavior" should be susceptible to mutant cheater clones that conserve energy by refusing to cooperate.In the case of cancer, the cheater clone would be identified histologically as a local (potentially extensive) hypoxic necrosis caused by microvascular hypoplasia.Given that all angiogenic tumors are potentially susceptible to cheaters, one expects either that hypertumor necrosis is fairly common or that some general mechanisms exist to prevent it.
This model suggests a simple, but nontrivial, hypothesis for the latter possibility.In the previous model, which represented small tumors limited only by access to microvessels [42], natural selection always favored clones with the highest proliferative potential (difference between division and cell death rates), but was neutral towards angiogenic growth factor secretion.If ATP allocation were zero-sum-that is, ATP allocated to secretion necessarily decreases ATP available for proliferation-then natural selection would secondarily disfavor angiogenic clones since these clones have less energy to support proliferation.The current model, which includes a more realistic energy management system, implies that this simple picture is incorrect.ATP allocation is not zero-sum.Reallocation of ATP away from proliferation in favor of angiogenic secretion can at the same time increase ATP available for proliferation due to a nonlinear interaction between ATP hydrolysis rate and total cell adenylate.Cells in a "resting" state, represented here by the default parameters, respond to increasing ATP hydrolysis by slowing adenylate destruction and thereby elevating total cell adenylate (Fig. 7).Homeostatic mechanisms, primarily glycolysis, then buoy up energy charge by transferring AMP to ATP.As a result, the cell has more ATP than before, and by default allocates more to proliferation even though no mutations acted directly on proliferative ability.As in the previous model, natural selection acts exclusively on proliferative potential.However, this model predicts that individual selection acting on clonal lineages can favor mutations that increase angiogenesis, but only via a pleiotropic effect on the energetic, and therefore proliferative, phenotype expressed by the clone.Cheaters cannot benefit by halting secretion because by doing so they decrease total adenylate and therefore ATP available for proliferation.However, the picture just painted applies only in relatively permissive environments.In more challenging environments, the homeostatic mechanisms switch into a mode in which ATP allocation becomes zero-sum or worse.As originally suggested by Atkinson [6], when the cell has difficulty adjusting total adenylate upwards to compensate for ATP hydrolysis-for example from ischemia (Fig. 9) or some increase in ATP hydrolysis rate (Fig. 7)-it maintains energy charge by increasing the AMP destruction rate (see equations (11), (12) and Fig. 2).Therefore, in such a situation, any increase in ATP hydrolysis for angiogenesis results in loss of total adenylate, total ATP and decreased proliferative potential.Therefore angiogenesis is selected against in a highly proliferative clone (high η p ).In such clones, the fitness gradient for angiogenesis potential (η s ) is negative for all physiologically reasonable values (Fig. 16).Therefore, only clones that down-regulate angiogenesis secretion are favored, resulting in runaway selection for ever-decreasing angiogenic potential.
Here we have a classic conflict between individual and group (tissue) level benefits [11,50,56,57].The group benefits from angiogenesis, but individuals benefit from proliferation.Because cells in nutrient-poor tumors cannot adjust total adenylate to compensate for increased ATP hydrolysis, proliferation and angiogenesis become antagonistic traits.Since selection acts exclusively on individual (clonal) traits, proliferation evolves at the expense of the group and its angiogenic potential.Given enough time, individual-level selection will eventually produce a favored lineage incapable of generating sufficient angiogenesis to support the tumor, resulting in evolutionary suicide [48] by a "hypertumor," or tumor-on-a-tumor [42,43,44].Such a situation applies to a tumor proliferating at its default ESS (Fig. 17), suggesting that a non-angiogenic clone, and therefore a hypertumor, is the ultimate evolutionary endpoint when both proliferative (η p ) and angiogenic (η s ) potential are under selection, although the techniques employed here cannot demonstrate this conjecture.
Selection for angiogenic potential's pleiotropic benefits generates another interesting prediction.Although this model typically admits an ESS for an angiogenic phenotype in the zone of tumor viability, the ESS is as an evolutionary repellersmall perturbations on either side cause runaway selection in the direction of the perturbation.Perturbations "to the left" of the ESS result in evolutionary suicide, as just described.Perturbations to the right spark runaway selection towards ever-increasing angiogenic signalling, resulting in massive microvessel hyperplasia ("evolutionary hypervascularization").In this case, evolutionary suicide is again inevitable in this model because extreme secretion rates degrade cellular energy stores beyond that required to maintain osmotic balance in the cell.In reality, this result is unlikely since other biochemical and biophysical limits besides cell death exist on angiogenic potential.
Nevertheless, such evolutionary hypervascularization could be responsible for microvascular hyperplasia characteristic of high-grade gliomas.Traditionally the explanation for microvascular hyperplasia in glioma is hypoxia, perhaps caused by thrombosis [10,32,53,54].According to this model, evolutionary hypervascularization should be limited to tumors with relatively low proliferative potential.Highly proliferative clones hydrolyze ATP a rates significantly above resting cells.Such clones therefore tend to erode total adenylate as decribed above, removing ATP for angiogenesis signaling.Note that this prediction does not imply that the phenomenon occurs only in slowly growing tumors.Tumor growth rate is a function of vascularization and nutritional support in addition to proliferative abilities conferred by genotype.To be precise, the prediction is that, for clones compared in the same environments, less proliferative genotypes would be more susceptible to evolutionary hypervascularization than would more proliferative clones.A relatively low-proliferative tumor hypervascularized by this evolutionary mechanism may nevertheless grow rapidly because of its dense nutritional infrastructure.
The hypotheses distilled from this model depend critically on one key behavior of the energetics model, namely the counterintuitive increase in total adenylate resulting from increased ATP hydrolysis.This physiological behavior has empirical support-indeed, it was precisely this behavior in erythrocytes that led Ataulakhanov, Vitvitsky, Martinov and their colleagues to formulate their adenylate model [2,3,4,5,38], which largely inspired model (3).This compensatory response creates the conditions making the angiogenic singular strategy an evolutionary repeller.Therefore, the model's key prediction is largely robust to changes in the adenylate model as long as the compensatory mechanism remains.For example, there is good empirical support for the assumption that total adenylate is regulated primarily via AMP [38], at least in erythrocytes where adenylate dynamics are relatively easily studied.However, total adenylate in model ( 3) is regulated exclusively by feedback on AMP destruction.Alternative models in which adenylate is regulated by its production rate have been around for some time [28,29] and would make an interesting alternative to the formulation used here.However, as long as resting cells increase adenylate concentration in response to increased ATP hydrolysis, the model's evolutionary predictions will remain unchanged.

Figure 2 .
Figure 2. AMP destruction function f (A 1 , A 3 ) for A 3 = 1.5027, the equilibrium value under the default parameters.AMP destruction is dominated by 5 nucleotidases when [AMP] falls below 10 −1 .Above this threshold, AMP deaminases dominate.Adapted from[4,38].Dashed horizontal line represents f (A 1 , A 2 ) = α a for default value of α a and the default equilibrium value of A 3 .Dashed vertical line represents default equilibrium value of A 1 .

Figure 3 .
Figure 3. Adenylate and energy charge dynamics of a cell reproducing at its maximum rate.Parameters and initial conditions were set to default values listed in table 1 except that η p = 7, its theoretical maximum, and ηs = 0. (A) Total adenylate (top line), ATP, ADP and AMP concentrations.(B) Energy charge.

Figure 4 .
Figure 4. Properties of a healthy tissue.All parameters are set to their default values listed in Table 1.(A) Parenchyma (Φ) and VEC (Ψ) growth functions as functions of tumor vascularization.(B) Total adenylate and ATP concentration per cell at equilibrium as a function of vascularization. .

4 . 4 . 1 .
Results.In this section we begin with a treatment of the adenylate dynamics predicted by model (3) before moving to the tissue-level model and its evolutionary dynamics.Time scales.From Table3we can estimate the orders of magnitude for the parameters in the adenylate dynamics model: a, b, c = O(1), G = O(10 2 ), α a = O(10 −5 ) and k = O(10 6 ).In addition, we can simplify the adenylate destruction

Figure 5 .
Figure5.Adenylate dynamics on time scales τ 1 and τ 2 (see text).Red surface is the graph of the manifold on which X = 0 (adenylate kinase in equilibrium) and the translucent grey manifold is the set of points satisfying Y = constant (conservation of adenylate).Black curve is a numerical solution of the full model (equations[3],[4],[6],[7],[11] and[12]) with default parameters and no simplifying assumptions.I) Dynamics on time scale τ 1 , representing equilibration of adenylate kinase (note that this portion of the trajectory lays above the red manifold); II) dynamics on time scale τ 2 .Dot is the initial condition.Compare Fig.6.

Figure 7 .
Figure 7. Equilibration of adenylate with changes in proliferation committment.Red curve is the graph of the solution to α a = f (A 1 , A 3 ).Black trajectories are projections of solutions of model (3) (no simplifying assumptions) onto the A 1 , A 3 plane for values of η p ∈ {0, 1, . . ., 12} with v = 1 and all other parameters set to default values.Solid blue lines are the initial conditions, which intersect at the equilibrium associated with η p = 0. Dashed blue line is the graph of equation(22) for η p = 12.Fixed points are black circles, blue arrows show progression of fixed points with increasing η p and black arrows show direction of solution trajectories.Evident are dynamics on time scales τ 2 (rapid movement from initial condition to cusp, left to right) and τ 3 (slow adjustment along the line defined by equation (22)).Dynamics on time scale τ 1 are irrelevant here since the initial conditions satisfy X = 0. Compare Fig.8A.

Figure 8 .
Figure 8. Cellular energy status as a function of ATP allocated to proliferation (η p ) over twice its estimated physiological range.(A) Equilibrium ATP and total adenylate concentrations.(B) Energy charge.All other parameters were set to their values indicated in Table1, with v = 1 (normal tissue vascularization).

Figure 9 .
Figure 9. Equilibration of adenylate with changes in vascularization.Curves as in figure7with the following exceptions: here η p = 1 and v now varies on the set {0.1, 0.2, . . ., 1.5}.Concentration of AMP declines with v in this range.

Figure 11 .
Figure 11.Saddle-node bifurcation with increasing commitment to proliferation.Here we view η p as the bifurcation parameter.Red and blue curves represent Ψ and Φ, respectively, for three different values of η p : 0.2, 0.255 and 0.3.The proliferation function Φ increases uniformly with η p .For η p < 0.255, Φ < Ψ and there are no biologically relevant fixed points.At about η p ≈ 0.255, we see a saddle-node bifurcation at the point where Φ ≈ Ψ.As η p continues to increase, saddle and node separate with the saddle to the right.

Figure 12 .
Figure 12.Fitness gradient for η p (ATP investment to proliferation) throughout its range of viability.Grey regions on either side represent regions in which the tumor is inviable.Dashed vertical line marks the ESS.All other parameters are set to default values.

Figure 14 .
Figure 14.Fitness gradient for ATP investment in TAF secretion ( ηs ) throughout the zone of viability.Grey areas represent regions of tumor inviability.Vertical dashed line indicates the singular strategy.All other parameters are fixed at their defaults; in particular, η p = 1.

Figure 15 .
Figure 15.Close-up view of pairwise invasibility plot near ESS for for ATP allocation to angiogenesis factor secretion (η s ) under default parameters.The PIP remains qualitatively the same throughout the range of biologically reasonable values of ηs .

Figure 16 .
Figure 16.Fitness gradients for ATP investment in TAF secretion ( ηs ) for various proliferation investments (η p ∈ {1.2, 1.3, 1.4, 1.5, 1.6}).Vertical dashed lines indicate singular strategies (all of which are ESSs), and arrows indicate the direction of increasing η p .All other parameters are fixed at their defaults.

Figure 17 .
Figure 17.Fitness gradient for ATP investment in TAF secretion ( ηs ) at the default ESS proliferation investments (η p ≈ 3.67).Gray region on the left is the tumor inviability region.

Table 1 .
Dependent variables and functions used in this research.

Table 2 .
Parameters and their default values, representing a "resting" cell.Sources are given in section 3.
2 we established that, for such a cell,