LETTER doE10.1038/nature14971

LETTER doE10.1038/nature14971 A spatial model predicts that dispersal and cell turnover limit intratumour heterogeneity Bartlomiej WaclawI, Ivana Bozic 13, Meredith E. Pittman'', Ralph H. Hruban4, Bert Vogelstein“ & Martin A. NowaIc2-1.6 Most cancers in humans are large, measuring centimetres in diameter, and composed of many billions of cells'. An equivalent mass of normal cells would be highly heterogeneous as a result of the mutations that occur during each cell division. What is remark- able about cancers is that virtually every neoplastic cell within a large tumour often contains the same core set of genetic altera- tions, with heterogeneity confined to mutations that emerge late during tumour growth'''. How such alterations expand within the spatially constrained three-dimensional architecture of a tumour, and come to dominate a large, pre-existing lesion, has been unclear. Here we describe a model for tumour evolution that shows how short-range dispersal and cell turnover can account for rapid cell mixing inside the tumour. We show that even a small selective advantage of a single cell within a large tumour allows the descendants of that cell to replace the precursor mass in a clinically relevant time frame. We also demonstrate that the same median- isms can be responsible for the rapid onset of resistance to chemo- therapy. Our model not only provides insights into spatial and temporal aspects of tumour growth, but also suggests that target- ing short-range cellular migratory activity could have marked effects on tumour growth rates. Tumour growth is initiated when a single cell acquires genetic or epigenetic alterations that change the net growth rate of the cell (birth minus death), and enable its progeny to outgrow surrounding cells. As these small lesions grow, the cells acquire additional alterations that cause them to multiply even faster and to change their metabolism to survive better the harsh conditions and nutrient deprivation. This progression eventually leads to a malignant tumour that can invade surrounding tissues and spread to other organs. Typical solid tumours contain about 30-70 clonal amino-acid-changing mutations that have accumulated during this multi-stage progression'. Most of these muta- tions are believed to be passengers that do not affect growth, and only —5-10% are drivers that provide cells with a small selective growth advantage. Nevertheless, a major fraction of the mutations, particu- larly the drivers, are present in 30-100% of neoplastic cells in the primary tumour, as well as in metastatic lesions derived from it's. Most attempts at explaining the genetic make-up of tumours assume well-mixed populations of cells and do not incorporate spatial constraintr'°. Several models of the genetic evolution of expanding tumours have been developed in the past"-'4, but they assume either very few mutations' or one- or two-dimensional growth".". Conversely, models that incorporate spatial limitations have been developed to help to understand processes such as tumour metabolism 1s, angiogenesism" and cell migration", but these models ignore gen- etics. Here, we formulate a model that combines spatial growth and genetic evolution, and use the model to describe the growth of primary tumours and metastases, as well as the development of resistance to therapeutic agents. We first model the expansion of a metastatic lesion derived from a cancer cell that has escaped its primary site (for example, breast or colorectal epithelium) and travelled through the circulation until it lodged at a distant site (for example, lung or liver). The cell initiating the metastatic lesion is assumed to have all the driver gene mutations needed to expand. Motivated by histopathological images (Fig. la), we model the lesion as a conglomerate of balls of cells (see Methods and Extended Data Fig. 1). Cells occupy sites in a regular three-dimen- sional lattice (Extended Data Fig. 2a, b). Cells replicate stochastically with rates proportional to the number of surrounding empty sites d Figure I I Structure of solid neoplasms. a, Hepatocellular carcinoma composed of balls of cells (circled in green) separated by non-neoplastic tissue (asterisk). b, Adjacent section of the bottom tumour in a immunolabelled with the proliferation marker Ki67. The edge ofthe tumour is delineated in red; the centre is marked with a green circle. Proliferation is decreased in the centre when compared to the edge of the neoplasm. c, d, Higher magnification of the centre (c) and the edge (d) with each proliferating neoplastic cell marked by a green dot. The blue nuclei without green dots are non-proliferating. The red circle in c demonstrates an example of cells (inflammatory cells) that were not included in the count of neoplastic cells. The neoplastic tissue in d is above the red line; non-neoplastic (normal liver) is below the red line. Comparison of c with d shows that proliferation of neoplastic cells is decreased in the centre as compared to the edge of the lesion (quantified in Extended Data Table I). 'School of Physics end Astronomy. Univesity of Edinburgh Olt Peter Guthrie Tad Reed. Echnburei CH9 ND. UK.2Prowaen for Evolutonery Dynamics, Harvwd University. One Battle Square. Cambridge. Massachusetts 02138. USA. aDepartment et Mathematics. lianrard thwersity. One Oxford Street. Cambridee. Massachusetts 02133. USA' he Sol Goldman Pancreatic Cancer Research Center. Department et Palk...sly. Johns Flocains University Sdicee of Medicine. 401 North Broadway. Weinberg 2202. Baltimore. Maryland 21231.USA SLucleig Center and Howard Hushes medial InSlittite, Johns Hopkins Kimmel Cancer Center. 1650 Oriws Street Bellmore. Maryland 21287. USA. `Department oe Organismicand Evolubonary Ulm/. Renard Unmake. 26 Word Street Cambridge. Massashusetts 02138. USA. 00 MONTH 2015 I VOL 000 I NATURE I 1 0.2015 Macmillan Publishers Limited. Al rights reserved EFTA00603737 RESEARCH LETTER (non- neoplastic cells or extracellular matrix), hence replication is faster at the edge of the tumour. This is supported by experimental data (Fig. I b-d and Extended Data Table I). A cell with no cancer cell neighbours replicates at the maximal rate of b = In(2) = 0.69 days-I, in which b denotes the initial birth rate, equivalent to 24h cell- doubling time, and a cell that is completely surrounded by other cancer cells does not replicate. Cells can also mutate, but we assume all muta- tions are passengers (they do not confer fitness advantages). After replication, a cell moves with a small probability (M) to a nearby place close to the surface of the lesion and creates a new lesion. This 'sprout- ing of initial lesions could be due to short-range migration after an epithelial-to-mesenchymal transition" and consecutive reversion to a non-motile phenotype. Alternatively, it could be the result of another process such as angiogenesis (Methods), through which the tumour gains better access to nutrients. The same model governs the evolution of larger metastatic lesions that have already developed extensive vas- culature. Cells die with a death rate (d) independent of the number of neighbours, and are replaced by empty sites (non-neoplastic cells within the local tumour environment). If there is little dispersal (M .6 0), the shape of the tumour becomes roughly spherical as it grows to a large size (Fig. 2a and Supplementary Video 2). However, even a very small amount of dispersal markedly affects the predicted shape. For M> 0, the tumour forms a conglom- erate of `balls' (Fig. 2b, Extented Data Fig. 2c and Supplementary Video 3), much like those observed in actual metastatic lesions, with the balls separated by islands of non- neoplastic stromal cells mixed with extracellular matrix. In addition to this remarkable change in topology, dispersal strongly affects the growth rate and doubling time of the tumour. Although the size (N) of the tumour increases with time (7) from initiation as without dispersal (Extended Data Fig. 3a, b), it grows much faster (—exp(const X 7) for large 7) when M> (Fig. 2c). This also remains true for long-range dispersal in which M affects the probability of escape from the primary tumour into the circulation to create new lesions in distant organs (metastasis). Using plausible estimates for the rates of cell birth, death and dispersal probability, we calculate that it takes 8 years for a lesion to grow from one cell to one billion cells in the absence of dispersal (Al = 0), but less than 2 years with dispersal (Fig. 2c). The latter estimate is consistent with experimentally determined rates of metastasis growth as well as clinical experience, while the conventional model (without dis- persal) is not. a 1010 102 0 10 20 30 40 50 60 T (months) M=10-5Non-spatial models point to the size of a tumour as a crucial deter- minant of chemotherapeutic drug resistance'''. To determine whether a spatial model would similarly predict this dependency in a clinically relevant time frame, we calculated tumour regrowth prob- abilities after targeted therapies. We assume that the cell that initiates the lesion is susceptible to treatment, otherwise the treatment would have no effect on the mass, and that the probability of a resistant mutation is 10-' (Methods); only one such mutation is needed for a regrowth. Figure 3a shows snapshots from a simulation (Supplementary Video I) performed before and after the administration of a typical targeted therapy at time T= 0. At first, the size of the lesion (-3 mm at T= 0) rapidly decreases, but I month later resistant clones begin to proliferate and form tumours of microscopic size. Such resistant sub- clones are predicted to be nearly always present in lesions of sizes that can be visualized by clinical imaging techniquesS 12. By 6 months after treatment, the lesions have regrown to their original size. The evolu- tion of resistance is a stochastic process—some lesions shrink to zero and some regrow (Extended Data Fig. 4a). Figure 3b, c shows the probability of regrowth versus the time from the initiation of the lesion to the onset of treatment upon varying net growth rates b-d and dispersal probabilities. Regardless of growth rate, the capacity to migrate makes it more likely that regrowth will occur sooner, particu- larly for more aggressive cancers, that is, those which have higher net growth rates (Fig. 3b). This conclusion is in line with recent theoretical work on evolving populations of migrating cells''. If resistant muta- tions additionally increase the dispersal probability before or during treatment, regrowth is faster (Extended Data Fig. 4b, c). Having shown that the predictions of the spatial model are consist- ent with metastatic lesion growth and regrowth times, we turn to primary tumours. In contrast to metastatic lesions, here the situation is considerably more complex because the tumour cells are continually acquiring new driver gene mutations that can endow them with fitness advantages over adjacent cells within the same tumour. Our model of a primary tumour assumes that it is initiated via a single driver gene mutation that provides a selective growth advantage over normal neighbouring cells. Each subsequent driver gene mutation reduces the death rate as d = b(l — s)k, in which k is the number of driver mutations in the cell (k a 0, and s is the average fitness advantage per driver. Almost identical results are obtained if driver gene muta- tions increase cell birth rather than decrease cell death, or affect both Figure 21 Short-range dispersal affects size, shape and growth rate of tumours. a. b. A spherical lesion in the absence of dispersal (M = 0) (a) and a conglomerate of lesions (b), each initiated by a cell that has migrated from a previous lesion, for low but non-zero migration (M = 10-6). Colours reflect the degree of genetic similarity cells with similar colours have similar genetic alterations. The death rate is d = 0.8b, which corresponds to a net growth rate of 0.26 = 0.14 days- I, and N = 10' cells. c, Dispersal (M > 0) causes the tumour to grow faster in time. Each point = 100 sanples, error bars (too small to be visible) are In. Continuous lines (extrapolation) arc 6,000 X 1 00.43T (gee,) 1,000 X 100 77. (blue). 2 I NATURE I VOL 000 I 00 MONTH 2015 COM Macmillan Publishers Limited. All rights reserved EFTA00603738 LETTER RESEARCH 7 r.0 .6 44444 -2 MAIM OS• b 100 X 80 60 t' 40 °: 20 , I 11 1 laic!! 0 1st.. 2 4 6 8 10 12 T (month) 0 Net growth rate 0.56 0.34 days" Net growth rate 0.16 0.07 days', ' 81" 'r ' 100 Ir.Id [II I_ 00 . I. IA .0 IA - 10. r U I . r m . o IAI I .. IA-104 M - 10 ' Cl` 20 IA .10 ' 0 ' In 10 20 30 40 T(month) Figure 3 I Treatment success rates depend on the net growth rate of tumours. a, Time snapshots before and during therapy (M = Resistant subpopulations that cause the tumour to regrow after treatment can be seen at T = I month. b, c, Probability of tumour regrowth (P,g,„,$) as a function of time after treatment initiation, for different dispersal probabilities (M) and net growth rates of the resistant cells. A higher net growth rate (b) leads to a high regrowth probability, so that 50% of tumours regrow 6 months after treatment is initiated when M = I 0-5. c., Tumours with lower net growth rates require >20 months to achieve the same probability of regrowth. Number of samples = Ito 800 per point (282 on average). Error bars are M. Sec Methods for details. cell birth and cell death (Extended Data Fig. 514; the most important parameter is the fitness gain, s, conferred by each driver mutation. Figure 4a shows that in the absence of any new driver mutations (as for a perfectly normal cell growing in utero), donal subpopulations would be restricted to small, localized areas. Each of these areas has at least one new genetic alteration, but none of them confers a fitness advantage (they are 'passengers'). In an early tumour, in which the centre cell contains the initiating driver gene mutation, the same struc- ture would be observed—as long as no new driver gene mutations have yet appeared. The occurrence of a new driver gene mutation, however, markedly alters the spatial distribution of cells. In particular, the het- erogeneity observed in normal cells (Fig. 4a) is substantially reduced (Fig. 4b and Supplementary Video 5). The degree of heterogeneity can be quantified by calculating the number of genetic alterations (passen- gers plus drivers) shared between two cells separated by various dis- tances (Fig. 4d-f). The genetic diversity is markedly decreased (Fig. 4e), even with relatively small fitness advantages (s = 1%). This also has implications for the number of genetic alterations that will be present in a macroscopic fraction (for example, >50%) of all cells. Figure 4f shows that this number is many times larger for s = I% than s = 0%. Furthermore, our model predicts that virtually all cells within a large tumour will have at least one new driver gene mutation after 5 years of growth (Extended Data Fig. 5a).The faster the clonal expansion occurs (the larger s is), the smaller the number of passenger mutations (Extended Data Fig. 5d, e). Our results are also robust to changes to the model (Methods and Extended Data Figs 5 and 6). We stress that an important prerequisite for limiting heterogeneity is cell turnover in the tumour, because in the spatial setting cells with driver mutations can 'percolate through the tumour only if they replace other cells. In the absence of cell turnover, tumours are much more heterogeneous (Extended Data Fig. 6d). In summary, our model accounts for many facts observed clinically and experimentally. Our results are robust and many assumptions can be relaxed without qualitatively affecting the outcome (Methods and Supplementary Information). Although tumour cell migration has s-0% b s-1% d 2 40 S ! 44.1.4 G A As c ! 2o 00 0.6 1.0 1.5 2.0 s-0% s-1% E 3. N • lor is 2.5 s-1% .6-2.0 13$13 Z 1. Figure 4 Genetic diversity is strongly reduced by the emergence of driver mutations. a-f. For all, M = 0 and the initial net growth rate = 0.007 clays-1 (d = 0.996). The three most abundant genetic alterations (GAs) have been colour-coded using red (R), green (G) and blue (8) (c). Each section is 80 cells thick. Combinations of the three basic colours correspond to cells having two or three of these genetic alterations. a, No drivers—separated, conical sectors emerge in different parts of the lesion, each corresponding to a different clone. b, Drivers with selective advantage s = 1% lead to clonal expansions and many cells have all three genetic alterations (white area). d Genetic diversity can be determined quantitatively by randomly sampling pairs of cells separated by distance r and counting the number of shared genetic alterations. e, The number of shared genetic alterations versus the normalized distance il<r> decreases much more slowly for the case with (red) than without (blue) driver mutations. f, The total number of genetic alterations present in at least 50% of all cells is much larger for s = 1% than for s = 0%. Number of samples = 50 per data point. Error bars arc M. historically been viewed as a feature of cancer associated with late events in tumorigenesis, such as invasion through basement mem- branes or vascular walls, this classical view of migration pertains to the ability of cancer cells to migrate over large distance?'. Instead, our analysis reveals that even small amounts of localized cellular move- ment are able to markedly reshape a tumour. Moreover, we predict that the rate of tumour growth can be substantially altered by a change in dispersal rate of the cancer cells, even in the absence of any changes in doubling times or net growth rates of the cells within the tumour. Some of our predictions could be experimentally tested using new cell labelling techniques'''. Our results could also greatly inform the interpretation of mutations in genes whose main functions seem to be related to the cytoskeleton or to cell adhesion rather than to cell birth, death, or differentiationn". For example, cells that have lost the expression of E-cadherin (a cell adhesion protein) are more migratory than normal cells with intact E-cadherin expression", and loss of E-cadherin in pancreatic cancer has been associated with poorer pro- gnosis", in line with our predictions. Online Content Methods along with any additional Extended Data display items and Source Data. are available in the online version of the paper: references unique to these sections appear only in the online paper. Received 1 September 2014: accepted 23 July 2015. Published online 26 August 2015. 1. Vogelstein.B. eta?. Cancergenome landscapes. Science 339,1546-1558(2013). 2. Yachida. S. et aL Distant metastasis mass late during the genetic evolution of pancreatic cancer. Nature 467,1114-1117 (2010) 3. Sottoriva.A.et at Intiatumor heterogeneity in human glicblastoma reflects cancer evolutionanidynamia, Proc. Nat I Acad. Sci USA 110.4009-4014 (2013) 4. Navin. N. eta?. Tumour evolution inferred by single-cell sequencing. NaMe 472, 90-94 (2011) 5. Gerlinger. M. et at Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. N. Engl. I Med 366,883-892 (2012). & Gatenby. R.A & Vincent.T.L.An evolutionary modelof carcinogenesis.CancerRes 63.6212-6220 (2003). 00 MONTH 2015 I VOL 000 I NATURE 13 VMS Macmillan Publishers Limited. All rights reserved EFTA00603739 RESEARCH LETTER 7. Johnston. M. D.. Edwards. C. M. Balmer. W. F. Maini. P. K. & Chapman, S.J. Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer. Proc. Nan Aced Sri. USA 104.4008-4013 (2007). 8. Bozic. Let Accumulation of driver and passenger mutations during tumor progression. Proc. Nat I Acad. Sci USA 107.18545-18550 (2010). 9. Beerenwinkel. N.et Genetic progression and the waiting time to cancer. PLOS Cam put aid 3, e225 (2007). 10. Durrett. R. & Moseley. S.Evolution of resistance and progression to disease during clonal expansion of cancer. Theor. Poput. Bid 77, 42-48 (2010). 11. Gonzalez-Garcia. L Sole. RV. & Costa. J. Metapopulation dynamics and spatial ----' • •-- ' 1A-c A -`4-c-• "cA cricA "cc" c/AcIn` 13. Martens. E. A. Kostadmov. R. Maley. C. C. & Hailatschek O. Spatial structure increases the waiting time for cancer. Newt Phys. 13.115014 (2011} 14. Andersco.A. R. A. Weaver. A. M.. Cummings. P. T.& Quaranta. V. Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. CM' 127, 905-915 (2006). 15. Kim.Y..Magdalena.A.S&Othrner.H.G.A hybrid model kr tumor. spheroid growth in vitro I: theoretical development and early results. Math.Models Methods Ap#.Sci. 17.1773-1798(2007). 16. McDougall. S R., Anderson, A. R. & Chaplain. M.A Mathematical modeling of dynamic adaptive tumour-induced angicgenesis: clinical implications and therapeutic targeting strategies. J. Meer. Blosi. 241, 564-589 (2006). 17. Hawkins-Deemer, A.. Rockne. R. C..Anderson, A. RA & Swanson. K FL Modeting tumor-associated edema in gliomas during anti-angiogenic therapy and its imnact nn imao4:444tornnr Pm& Chun( 2 SA ontm, ISOZIC. Liven. U.& NOWalt. IA A LlytlatOKS Ot targeted cancer tnerapy. MOOS I1101. Med.18. 311-316 (2012). 21. Bozk, L& Nowak M. A.Timirg and heterogeneity of mutations associated with drug resistance in metastatic cancers.Proc. Natf Acad. Sci. USA 111. 15964-15968 (2014). 22 Turke, A. B. et aL Preexistence and clonal selection of MET amplification in EGFR mutant NSCLC. Cancer CO 17, 77-88 (2010). 24. Talmadge.J. E. & I.J.AACR Centennial Series: the !hom of cancer metastasis: historical perspective. Cancer Res. 70, 5649-5669 (2010). 25. Alcolea....etaL Differentiation imbalance in singJe oesophageal progenitor cells causes clonal immortalization and field change. Nature Cell 8101.16.615-622 (2014). 2& Weber. K. eta?. RGB marking facilitates multicolorclonal cell tracking, Nature Med 17, 504-509 (2011). 27. Bordeleau.F.. Akoser. T. A. & Reinharl-King C.A. Physical biology in cancer. 5. The rocky road of metastasis: the role of cytoskeletal mechanics in cell migratory res- ponseto3Dmatrixtopcgraphy.AmJ.Physia CettPhysiol 306,C110 -C120(2014). 28 Lawson. C. D. & Burridge. K The on-off relationship of Rho and Rac during integrin-mediated adhesion and cell migration.Smaff GTPases 5.e27958(2014). 29. Gall, T. M. H.& Frampton.A EGene of the month: E-cadherin (CDH1).J. Dia Pathot. 66.928-932 (2013). 30. Winter. J. M. et af. Absence of E-cadherin expression distinguishes noncohesive from cohesive pancreatic cancer. Cfin. Cancer Res. 14, 412-418(2008). Supplementary Information is available in the online version of the paper. Acknowledgements Support from The John Templeton Foundation is gratefully acknowledged B.W. was supported by the Leverhulme Trust Eery-Career Fellowship. and the Royal Society of Edinburgh Personal Research Fellowship. 1.8. was supported by FoundationalQuestions in Evolutionary Bio Grant S.M.acknowledge support from The Virginia and M. Ludwig Fund for Cancer Research. The Lustgarten Foundation for Pancreatic Cancer Research. The Sd Goldman Center for Pancreatic Cancer Research. and Nil grants CA43460 and CA62924. Author Contributions I.B. and B.Y. designed the study. B.W.wrote the computer prams and made simulations. B.W.. I.B. andfl made analytic calculations.= and R.H.H.canied out experimental work All authorsdiscussed the results. The man uscr_S_was written primarily by B.W.,=. I.B. and BV.. with contributions from= and R.H.H. Author Information Re .nts and perrnissions information is available at The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper. Cares ndence and requests for materials should be addressed to= 4 I NATURE I VOL 000 100 MONTH 4015 INOIS Macmillan Publishers Limited. All rights reserved EFTA00603740 LETTER METHODS No statistical methods were used to predetermine sample size. Experiments were not randomized and investigators were not blinded to allocation during experi- ments and outcome assessment Spatial model for tumour evolution. Tumour modelling has a long tradition". Many models of spatially expanding tumours were proposed in the past"'"', but they either assume very fewn-m-s'-"."."-" or no new mutations at allumat".". or one- or two-dimensional growthwwn'". On the other hand, well-mixed models with several mutations'" do not often include space, and computa- tional models aimed at being more biologically realistic"' require too much computing resources (time and memory) to simulate realistically large tumours (N., 10' cells). Our model builds on the Eden lattice model' and combines spatial growth and accumulation of multiple mutations. Since we focus on the interplay of genetics. spatial expansion and short-range dispersal of cells. for simplicity we do not explicitly model metabolism". tissue mechanics, spatial heterogeneity of tis- sues, different types of cells present or angiogenesis". A tumour is made of non-overlapping balls (microlesions) of cells. Tumour cells occupy sites of a regular 3D square lattice (Moore neighbourhood, 26 neighbours). Empty lattice sites are assumed to be either normal cells or filled with extracellular matrix and are not modelled explicitly. Each cell in the model is described by its position and a list of genetic alterations that have occurred since the initial neo- plastic cell, and the information about whether a given mutation is a passenger. driver. or resistance-carrying mutation. A passenger mutation does not affect the net growth rate whereas a driver mutation increases it by disrupting tight regu- lation of cellular divisions and shifts the balance towards increased proliferation or decreased apoptosis. The changes can also be epigenetic and we do not distinguish between different types of alterations. We assume that each genetic alteration occurs only once ('infinite allele model"). The average numbers of all genetic alterations, driver and resistant genetic alterations produced in a single replication event are denoted by y, 7d. and y r. respectively. When a cell replicates, each of the daughter cells receives n new genetic alterations of each type (n being generally different in both cells)drawn at random from the Poisson probabilitydistribution: e-712(11.12)" POO— (I) n! in which x denotes the type of genetic alteration. In model A shown in Figs 2-1. replication occurs stochastically. with rate proportional to the number of empty sites surrounding the replicating cell. and death occurs with constant rate depending only on the number of drivers. We also simulated other scenarios (models B, C and D. see below). Driver mutations increase the net growth rate (the difference between proliferation and death) either by increasing the birth rate or decreasing the death rate by a constant factor 1 + s. in which s> 0. Dispersal is modelled by moving an offspring cell to a nearby position where it starts a new microksion (Extended Data Fig. la). Klicroksions repel each other; a 'shoving algorithm's" (Extended Data Fig. Ib) ensures they do not merge. Code availability. The computer code (available at httpifwww2.ph.ed.ac.uki —bwadawfcancer-code) can handle up to 1 X 109 cells. which corresponds to tumours that are clinically meaningful and can be observed by conventional medical imaging (diameter >1 cm). The algorithm is discussed in details in the Supplementary Information. It is not an exact kinetic Monte Carlo algorithm because such an algorithm would be too slow to simulate large tumours. A com- parison with kinetic Monte Carlo for smaller tumours (Supplementary Information) shows that both algorithms produce consistent results. Model parameters. The initial birth rate 6= In(2) 0.69 days "'. which corre- sponds to a 21 h minimum doubling time. The initial death rate d = 0...0.995b depends on the aggressiveness of the tumour (larger values = less aggressive lesion). In simulations of targeted therapy. we assume that. before treatment, = 0.69 days"' and d= 0.56 = 0.35 dayss '. whereas during treatment b= 035 days-' and d = 0.69 days-', that is, birth and death rates swap places This rather arbitrary choice leads to the regrowth time of about 6 months which agrees well with clinical evidence. Mutation probabilities are 7 = 0.02,7d = 1 X 10-5.7, = I X 10-7. in line with experimental evidence and theoretical work"'". Since there are no reliable data on the dispersal probability M. we have explored a range of values between M = 1 X 10-7 and I X 10-1. An parameters are summarized in Extended Data Fig. Ic, see also further discussion in Supplementary Information. Validity of the assumptions of the modeL Our model is deliberately oversim- plified. However, many of the assumptions we make can be experimentally jus- tified or shown not to qualitatively affect the model. Three-dimensional regular lattice of cells. The 3D Moore neighbourhood was chosen because it is computationally fast and introduces relatively fever artefacts related to lattice symmetries. Real tissues are much less regular and the number of nearest neighbours is different". However. recent simulations of similar models of bacterial colonies"' show that the structure of the lattice (or the lack thereof in off-lattice models) has a marginal effect on genetic heterogeneity. Asynchronous cell division. Division times of related cells remain correlated for a few generations. However, stochastic cell division implemented in our model is a good approximation for a large mass of cells and is much less computationally expensive than modelling a full cell cycle. Replication faster at the boundary than in the Interior. Several studies have described a higher proliferation rate at the leading edge of tumours, and this has been associated with a more aggressive clinical course". To estimate the range of values of death rate d for our model. we used the proliferation marker Ki67. Representative formalin-fixed, paraffin-embedded tissue blocks were selected from four small chromophobe renal cell carcinomas and six small hepatocellular carcinomas by the pathologist (=.). A section of each block was immunola- belled for IC67 using the Ventana Benchmark XT system. Around 8-12 images. depending on the size of the lesion, were acquired from each tumour. Fields were chosen at random from the leading edge and the middle of the tumour and were not necessarily 'hot spots' of proliferative activity. Using an Image) macro, each Ki67-positive tumour nucleus was labelled green by the pathologist. and each Ki67-negative tumour nucleus was labelled red. Other cell types (endothelium, fibroblasts and inflammatory cells) were not labelled. The proliferation rate was then calculated using previously descrthed methods". Statistical significance of the results was determined using a Kolmogorov-Smimov nap-sample test (signifi- cance level 0.05). The study was approved by the Institutional Review Board of the Johns Hopkins University School of Medicine. In all ten tumours, the proliferation rate at the leading edge of the tumour was grmter than that at the centre by a factor of 1.25 to 6 (Extended Data Table I). Comparing the density of proliferating cells to our model gives cf.., 0.56 (range: d = 0.176...0.86), which is what we assume in the simulations of aggressive lesions. Equal fitness of all cells in metastatic lesions. We assume that cells in a meta- static lesion are already very fit since they contain multiple driven. Indeed. studies of primary tumours and their matched metastases usually fail to find driver muta- tions present in the metastases that were not present in the primary lesions'', although there are notable exceptions, see, for example. refs 75 and 76. Experimental evidence in microbes" and (to a lesser extent) in eukaryotes" sug- gests that fitness gains due to individual mutations are largest at the beginning of an evolutionary process and that the effects of later mutations are much smaller. It remains to be seen how well these results apply to late genetic alterations in cancer" but if true. new driven occurring in the lesion are unlikely to spread through the population before the lesion reaches a clinically relevant size. Dispersant, our model. cells detach from the lesion and attach again at a different location in the tissue. This can be viewed either as cells migrating from one place to another one. or as a more generic mechanism that allows tumour cells to get better access to nutrients by dispersing within the tissue. hence providing a growth advantage over cells that did not disperse. Some mechanisms that do not involve active motion (that is. cells becoming motile) are discussed below. Migration. Cancer cells are known to undergo epithelial-to-mesenchymal trans- ition. the origin of which is thought to be epigenetic". This involves a cell becom- ing motik and miming some distance. If the cell finds the right environment, it can switch back to the non-motile phenotype and start a new lesion. Motility can be enhanced by tissue fluidization due to replication and death". Instead of mod- elling the entire cycle (epithelial-mesenchymal-epithelial). we only model the final outcome (a cell has moved some distance). Tumour buds. Many tumours exhibit focally invasive cell clusters, also known as tumour buds. Their proliferation rate is less than that of cells in the main turnout'. We propose that tumour buds contain cells that have not yet completed epithdial- to.meserichymal transition and therefore they proliferate slower. Single versus cluster migration. Ref. 82 found that circulating cancer cells can travel in clusters of 2-50 cells, and that such clusters can initiate metastatic foci. They report that approximately one-half of the metastatic foci they examined were initiated by single circulating cancer cells, and that circulating cancer cell dusters initiated the other half. The authors also note that the cells forming a duster are probably neighbouring cancer cells from the primary tumour. This means that the genetic make-up of cells within a newly established lesion will be very similar. regardless of its origin (single cell versus a small cluster of cells). Therefore, the ability to travel in clusters should not affect the genetic heterogeneity or regrowth probability as compared to single-cell dispersal from our modeL Angiogenesis. We do not explicitly model angiogenesis for two reasons. First. most genetic alterations that can either change the growth rate or be detected experimentally must occur at early stages of tumour growth as explained before. Hence, the genetic make-up of the tumour is determined primarily by what happens before angiogenesis. Second, local dispersal from the model mimics tumour cells interspersing with the vascularized tissue and getting better access to nutrients, which is one of the outcomes of angiogenesis. @2015 Macmillan Publishers Limited. AU rights reserved EFTA00603741 RESEARCH LETTER Biomechanics of tumours. Growth is affected by the mechanical properties of cells and the extracellular matrix. We do not explicitly include biomechanics (see, however, below). in contrast to more realistic models"". as this would not allow us to simulate lesions larger than about 1 X 106 cells. Instead, we take experimentally determined values for birth and death rates. values that are affected by biomecha- nics. as the parameters of our model. Isolated balls of cells. In our simulations. balls of cells are thought to be separated by normal, vascularized tissue which delivers nutrients to the tumour. The envir- onment of each ball is the same. and there are no interactions between the balls other than mechanical repulsion. This represents a convenient mathematical con- trivance and qualitatively recapitulates what is observed in stained sections of actual tumours (Fig. la). We investigated under which circumstances the balls of cancer cells would mechanically repel each other. see Extended Data Fig. 7 for a graphical summary of the results. We simulated a biomechanical, off-lattice model of normal tissue composed of 'ducts' lined with epithelial cells and separated by stroma (Supplementary Information. section 8). Mechanical interactions between cells were modelled using an approach similar to that described previously'". "s with model parameters taken from refs 59, 60.85-88. We assumed cancer cells to be of epithelial origin. as are most cancers". Cancer cells that invaded different areas of epithelium grew into balls that remained separated by thin slices of stroma (Supplementary Videos 8-11). This 'encapsulation' of tumour microlesions was possible owing to the supportive nature of stroma that is able to mechanically resist expansion of balls of cancer cells. Encapsulation is essential if the balls are to repel each other. If the tissue is 'fluidized' by random replication and death. the balk quickly merge (Supplementary Video 12). Another important factor are differ- ences in mechanical properties of tumour and normal cells": it is known that differences in cdlular adhesion and stiffness promote segregation of different types of cells",. In reality. mkrolesions within the primary tumour are less symmetric and some of them are better described as 'protrusions bulging out from the main tumour tissue, owing to biomechanical instabilities: see. for example, refs 93.94. However, stroma may still provide enough spatial separation. and the capillary network of blood vessels—either due to tumour angiogenesis or preexisting in the invaded tissue—may provide enough nutrients to the lesions so that our assumption of independently growing balls of cells remains valid. Therefore, we believe that modelling the tumour as a collection of non- or weakly-interacting microksions is essentially correct. We also note that the existence of isolated balls is not neces- sary to explain our qualitative results: reduced heterogeneity and increased growth in the presence of migration. Supplementary Video 13 shows that even if the tissue is homogeneous and highly dynamical and there are no isolated balls of cells, migration leads to a considerable speedup of growth as compared to the case with no migration (Supplementary Video 14). Tumour geometry and heterogeneity in the absence of driver mutations. Supplementary Videos 2 and 3 illustrate the process of growth of a tumour with maximally N= 107 cells, for Af = 0 and Af = l0"". respectively, and for d = 0.5. Extended Data Fig. 2 shows snapshots from a single simulation for M = 0, N— 10'. and d = 0 (no death. Extended Data Fig. 2a) and d = 0.9 (Extended Data Fig. 2b). In the latter case. cells are separated by empty sites (normal cells/ extracellular matrix). Extended Data Fig. 2c shows that the tumour is almost spherically symmetric for Af = 0. The symmetry is lost for small but non-zero Af, and restored for larger M when the balk become smaller and their number increases. Extended Data Fig. 2c also shows that metastatic tumours contain many donal sectors with passenger mutations. Extended Data Fig. 8a shows that the fraction G(r) of genetic alterations that are the same in two randomly sampled cells (Fig. 4) separated by distance r quickly decreases with r. indicating increased genetic heterogeneity owing to passenger mutations. Targeted therapy of metastatic lesions. Models of cancer treatments" "° often assume either no spatial structure or do not model the emergence of resistance. We assume that the cell that initiated the lesion was sensitive to treatment but its progeny may become resistant. Before the therapy commences. all cells have the same birth and death rates. but after the treatment resistant cells continue to proliferate with the same rate, whereas susceptible cells are assigned different rates as described above. Resistant cells can emerge before and during the therapy. The death rate of sensitive cells during treatment is greater than the birth rate. or the tumours would not be sensitive to the drug. For example. in Fig. 3 treatment increases the death rate and decreases the growth rate of susceptible cells, the growth rate of resistant cells after therapy is identical to that of the sensitive cells before treatment. d = 0.5b in the absence of treatment. elf = I0 6, and treatment begins when the tumour has N= 10' cells. Note that our model assumes the drug is uniformly distributed in the tumour": it is known that drug gradients can speed up the onset of resistance'. Supplementary Video 1 and Extended Data Fig. 4a show that. since the process of resistance acquisition is stochastic. some tumours regrow after an initial regression, and some do not. If only resistant cells can migrate, regrowth is faster (Extended Data Fig. 4b, c). Extended Data Fig. 4d-g shows regrowth probabilities P,,,,,„„th for different treatment scenarios not mentioned in the main text. depend- ing on whether the drug is cytostatic (bite„,tr,„„, = 0) o cytocidal = b). and whether d = 0 or d>0 before treatment. In Extended Data Fig. 4d, cells replicate and die only on the surface. and the core is 'quiescent' —cells are still alive there but cannot replicate unless outer layers are removed by treatment (Supplementary Videos 6 and 7). Pitp.0.0 does not depend on the dispersal prob- ability Af at all. and is close to 100% for N> 108 cells. a size that is larger than for d > 0 (Extended Data Fig. 40. It can be shown that P,,,,„,„h = 1 — exp(-7,4 Extended Data Fig. 4e is for the cytostatic drug (Oh...um= = 0): this is also equivalent to the eytocidal dnig if the tumour has a necrotic core (cells are dead but still occupy physical volume). In this case, increases with Af because more resistant cells are on the surface for larger M (cells can replicate only on the surface in this scenario). Extended Data Fig. 4f. g shows models with cell death present even in the absence of treatment (d = 0.9b) but occurring only at the surface. unlike in Fig. 3 where cells also die inside the tumour. Death increases owing to a larger number of cellular division necessary to obtain the same size and hence more opportunities to mutate. Relaxing the assumptions of the model. Figure 4 shows that even a small fitness advantage substantially reduces genetic diversity through the process of donal expansion, see also Supplementary Videos 4 and S. We now demonstrate that this also applies to modified versions of the model, proving its robustness. Exact values of Mend s has no qualltathv effect. Extended Data Fig. 8b, e shows that the average number of shared genetic alterations is larger in the presence of drivers also in the case of non-zero dispersal (elf >0), and its numerical value is almost the same as for M = 0 (Fig.4). Extended Data Fig. 8c. f shows that as long as s> 0 and regardless of its exact value, driver mutations reduce genetic diversity in the tumour compared to the case s = 0. Extended Data Fig. Sa-c shows how many driver mutations are expected to be present in a randomly chosen cell from a tumour that is Tyears old. Neither dispersal nor the way drivers affect growth (via birth or death rate) has a significant effect on the number of drivers per cell (Extended Data Fig. 56. c). A small discrepancy visible in Extended Data Fig. 5b is caused by a slightly asymmetric way death and birth is treated in our model, see the Supplementary Information. Model B. Cells replicate with constant rate if there is at least one empty neighbour. In the absence of drivers, genetic alterations are distributed evenly throughout the lesion (Extended Data Fig. 6b) but they often occur independently and the number of frequent genetic alterations is low (Extended Data Fig. 6e). Drivers cause clonal expansion as in model A. Model C. Cells replicate regardless of whether there are empty sites surrounding them or not. When a cell replicates, it pushes away other cells towards the surface (Supplementary Information). Extended Data Fig. bc,e shows that this again leads to clonal expansion which decreases diversity. Model D. Replication/death occurs only on the surface and the core of the tumour is static Extended Data Fig. 6d shows that driver mutations cannot spread to the other side of the lesion and conical clonal sectors can be seen even for s> 0. The number of frequent genetic alterations is the same fors = 0 and s = 1%, indicating that genetic heterogeneity is not lowered by clonal expansion. This demonstrates that cell turnover inside the tumour is very important for reducing heterogeneity. To obtain the same (low) heterogeneity as for models a-c, the probability of driver mutations must be much larger in model D (Extended Data Fig. 60. Drivers affecting M. We investigated three scenarios in which drivers affect (1) only the dispersal probability Al—, (1 + q)M. in which q> 0 is the 'migration fitness advantage (no change in b, d), (2) both Af and 4 that is, (d.AO— (d(1—s).(1+q)115) with s, q > 0, (3) either Al or d. with probability 1/2. Extended Data Fig. 3c shows that growth is unaffected in cases (1, 3) compared to the neutral case. For (2) the tumour growth rate increases significantly when the tumour is larger than N= I X 106 cells. This shows that migration increases the overall fitness advantage. in line with ref. 102. which shows that fixation probabil- ity is determined by the product of the exponential growth rate and diffusion constant (motility) of organisms. Six-site (von Neumann) neighbourhood. We simulated a model in which each cell has only six neighbours (von Neumann neighbourhood) instead of 26(Moore neighbourhood). Extended Data Fig. 9 compares models A and C for the two neighbourhoods and show that there is only a small quantitative difference in the growth curves for model A (model C is unaffected), but that the shape of the ball of cells deviates more from the spherical one for the six-site neighbourhood. see also section 7 in the Supplementary Information. 31. 4.2015 Macmillan Publishers Limited. All rights reserved EFTA00603742 LETTER RESEARCH 32. Anderson. A R. A. A hybrid mathematical model of solid tumour invasion: the imnartanra nt roll aelhoninn Math MM iaint 22 1F1-1 PA rinnci n. Lavrentoem. M. U. as Nelson U. K. survival prObablatieS at sonencaiirontiers spatialy structured tissue.). Math. Blot http://dxdoicrg/10.1007/s00285. 015-0912-1 (2015). 38. Gerlee. P.& Nelander. S. The impact of phenotypic switching on glioblastoma growth and invasion. PLOS Comput Biot. 8, e1002556 (2012). 39. Gonzalez-Garda, I, Solo. R V.& Costa.). Metapopulation dynamics and spatial heterogeneity in cancer. Prat Nan Acad. Sot USA 99,13085-13089 (2002). 40. Sehyo.C.C.efaL Model forin vivo progression of tumors based on co-evolving cell population and vasculature. Sd. Rep. 1, 31 (2011). 41. Torquato. S. Toward an 'sing model of cancer and beyond. Piro. Bid. 8, 015017 (2011). 42. Reiter.J. G.. Bozic. LAllen. B, Chattelee, K. & Nowak. M.A. The effect of one additional driver mutation on tumor progression. Eva 6, 34-45 (2013). 43. Kansal,A.R., Torquato.S.Harsh.G.R.Chixca.E A & Deisboecit. T. S. Simulated brain tumor growth dynamics using a three-dimensional cellular automaton. J. meat. Biol. 203, 367-382 (2000). 44. Kemal. A. R. Torquato.S..Chiacca. EA & Deisboeck T. S. Emergence of a subpopulation in a corn putational model of tumor growth.). Thew. Biol. 207, 431-441(2000). 45. Antal.T. Knpivsky, P. L. & Nowak M.A. Spatial evolution of tumors with successive driver mutations. Phys. Rev. £92.022705 (2015} 46. Endeding H, HLatky. L & tlahnfeldt. P. Migration rules: tumours are conglomerates of self-metastases. Br. J. Cancer 100, 1917-1925 (2009). 47. Sottoriva etaL Cancer stem cell tumor model reveals invasive morphology and increased phenotypical heterogeneity. Cancer Res. 70. 46-56 (2010). 48. Schaller. G. & Meyer-Hermann. M. Multicellular tumor spheroid in an off-lattice voronai-delaunay cell model Phys. Rev. E71, 051910 (2005} 49. Radszuweit. M. Block. M. Hengstler. J.G.Sctibll, E & Drasdo, D. Comparing the growth kinetics of cell populations in two and three dimensions. Phys. Rev. £79. 051907 (2009). 50. Moglia.R.Guisoni. N.& Albano, E.V. Interfacial properties in a discrete model for tumor. growth.Phys, Rev. E87, 032713 (2013} 51. Foo. J. Leder, K. & Ryser. M. Multifocality and recurrence risk: a quantitative model of field cancerization. J. moor. Biol. 355, 170-184 (2014). 52. PoleszczukJ.Hahnfeldt,P.&Enderling H. Evolution and phenotypic selection of cancer stem cells. PLOS Cornput Bid. 11, e1004025 (2015). 53. Durrett. R. Schmidt. D.& Schweinsberg. J. A waiting time problem arisi rg from the study of multi-stage carcinogenesisAnn.Appt hob& 19, 676-718 (2009). 54. Spencer, S. L. Eletryman. M. J.. Garda, J. A. &Abbott, D.M ordinary differential equation model for the multistep transformation to cancer.). Them Blot 231. 515-524(2004). 55. Kim, Y.& Othrner. H. G.A hybrid model of tumor-stromal interactions in breast cancer. &AL Math. Blot 75, 1304-1350 (2013). 56. Ramis-Conde, I, Chaplain. M.A. J., Anderson, A. RA & Drasdo. D. Multi-scale modelling of cancer cell intravasation: the role of cad herins in metastasis. Phys. Blot 6.016008 (2009} 57. Swanson. K. R. erg Quantifying the role of angicgenesis in malignant progr