EFTA00603737

LETTER doE10.1038/nature14971 A spatial model predicts that dispersal and cell turnover limit intratumour heterogeneity Bartlomiej WaclawI, Ivana Bozic13, Meredith E. Pittman'', Ralph H. Hruban4, Bert Vogelstein“ & Martin A. NowaIc2-1.6 Most cancers in humans are large, measuring centimetres in We first model the expansion of a metastatic lesion derived from a diameter, and composed of many billions of cells'. An equivalent cancer cell that has escaped its primary site (for example, breast or mass of normal cells would be highly heterogeneous as a result of colorectal epithelium) and travelled through the circulation until it the mutations that occur during each cell division. What is remark- lodged at a distant site (for example, lung or liver). The cell initiating able about cancers is that virtually every neoplastic cell within a the metastatic lesion is assumed to have all the driver gene mutations large tumour often contains the same core set of genetic altera- needed to expand. Motivated by histopathological images (Fig. la), we tions, with heterogeneity confined to mutations that emerge late model the lesion as a conglomerate of balls of cells (see Methods and during tumour growth'''. How such alterations expand within the Extended Data Fig. 1). Cells occupy sites in a regular three-dimen- spatially constrained three-dimensional architecture of a tumour, sional lattice (Extended Data Fig. 2a, b). Cells replicate stochastically and come to dominate a large, pre-existing lesion, has been with rates proportional to the number of surrounding empty sites 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 d 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. Figure I I Structure of solid neoplasms. a, Hepatocellular carcinoma Most attempts at explaining the genetic make-up of tumours composed of balls of cells (circled in green) separated by non-neoplastic tissue assume well-mixed populations of cells and do not incorporate spatial (asterisk). b, Adjacent section of the bottom tumour in a immunolabelled constraintr'°. Several models of the genetic evolution of expanding with the proliferation marker Ki67. The edge ofthe tumour is delineated in red; tumours have been developed in the past"-' 4, but they assume either the centre is marked with a green circle. Proliferation is decreased in the centre very few mutations' or one- or two-dimensional growth".". when compared to the edge of the neoplasm. c, d, Higher magnification of Conversely, models that incorporate spatial limitations have been the centre (c) and the edge (d) with each proliferating neoplastic cell marked developed to help to understand processes such as tumour metabolism1s, by a green dot. The blue nuclei without green dots are non-proliferating. The angiogenesism" and cell migration", but these models ignore gen- 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 etics. Here, we formulate a model that combines spatial growth and d is above the red line; non-neoplastic (normal liver) is below the red genetic evolution, and use the model to describe the growth of primary line. Comparison of c with d shows that proliferation of neoplastic cells is tumours and metastases, as well as the development of resistance to decreased in the centre as compared to the edge of the lesion (quantified in therapeutic agents. 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 Non-spatial models point to the size of a tumour as a crucial deter- faster at the edge of the tumour. This is supported by experimental minant of chemotherapeutic drug resistance'''. To determine data (Fig. I b-d and Extended Data Table I). A cell with no cancer cell whether a spatial model would similarly predict this dependency in neighbours replicates at the maximal rate of b = In(2) = 0.69 days-I, a clinically relevant time frame, we calculated tumour regrowth prob- in which b denotes the initial birth rate, equivalent to 24h cell- abilities after targeted therapies. We assume that the cell that initiates doubling time, and a cell that is completely surrounded by other cancer the lesion is susceptible to treatment, otherwise the treatment would cells does not replicate. Cells can also mutate, but we assume all muta- have no effect on the mass, and that the probability of a resistant tions are passengers (they do not confer fitness advantages). After mutation is 10- ' (Methods); only one such mutation is needed for a replication, a cell moves with a small probability (M) to a nearby place regrowth. close to the surface of the lesion and creates a new lesion. This 'sprout- Figure 3a shows snapshots from a simulation (Supplementary ing of initial lesions could be due to short-range migration after an Video I) performed before and after the administration of a typical epithelial-to-mesenchymal transition" and consecutive reversion to a targeted therapy at time T= 0. At first, the size of the lesion (-3 mm at non-motile phenotype. Alternatively, it could be the result of another T= 0) rapidly decreases, but I month later resistant clones begin to process such as angiogenesis (Methods), through which the tumour proliferate and form tumours of microscopic size. Such resistant sub- gains better access to nutrients. The same model governs the evolution clones are predicted to be nearly always present in lesions of sizes that of larger metastatic lesions that have already developed extensive vas- can be visualized by clinical imaging techniquesS12. By 6 months after culature. Cells die with a death rate (d) independent of the number of treatment, the lesions have regrown to their original size. The evolu- neighbours, and are replaced by empty sites (non-neoplastic cells tion of resistance is a stochastic process—some lesions shrink to zero within the local tumour environment). and some regrow (Extended Data Fig. 4a). Figure 3b, c shows the If there is little dispersal (M .6 0), the shape of the tumour becomes probability of regrowth versus the time from the initiation of the lesion roughly spherical as it grows to a large size (Fig. 2a and Supplementary to the onset of treatment upon varying net growth rates b-d and Video 2). However, even a very small amount of dispersal markedly dispersal probabilities. Regardless of growth rate, the capacity to affects the predicted shape. For M> 0, the tumour forms a conglom- migrate makes it more likely that regrowth will occur sooner, particu- erate of `balls' (Fig. 2b, Extented Data Fig. 2c and Supplementary larly for more aggressive cancers, that is, those which have higher net Video 3), much like those observed in actual metastatic lesions, with growth rates (Fig. 3b). This conclusion is in line with recent theoretical the balls separated by islands of non- neoplastic stromal cells mixed work on evolving populations of migrating cells''. If resistant muta- with extracellular matrix. In addition to this remarkable change in tions additionally increase the dispersal probability before or during topology, dispersal strongly affects the growth rate and doubling time treatment, regrowth is faster (Extended Data Fig. 4b, c). of the tumour. Although the size (N) of the tumour increases with time Having shown that the predictions of the spatial model are consist- (7) from initiation as without dispersal (Extended Data Fig. 3a, b), ent with metastatic lesion growth and regrowth times, we turn to it grows much faster (—exp(const X 7) for large 7) when M> primary tumours. In contrast to metastatic lesions, here the situation (Fig. 2c). This also remains true for long-range dispersal in which M is considerably more complex because the tumour cells are continually affects the probability of escape from the primary tumour into the acquiring new driver gene mutations that can endow them with fitness circulation to create new lesions in distant organs (metastasis). advantages over adjacent cells within the same tumour. Our model of a Using plausible estimates for the rates of cell birth, death and dispersal primary tumour assumes that it is initiated via a single driver gene probability, we calculate that it takes 8 years for a lesion to grow from mutation that provides a selective growth advantage over normal one cell to one billion cells in the absence of dispersal (Al = 0), but less neighbouring cells. Each subsequent driver gene mutation reduces than 2 years with dispersal (Fig. 2c). The latter estimate is consistent the death rate as d = b(l — s)k, in which k is the number of driver with experimentally determined rates of metastasis growth as well mutations in the cell (k a 0, and s is the average fitness advantage as clinical experience, while the conventional model (without dis- per driver. Almost identical results are obtained if driver gene muta- persal) is not. tions increase cell birth rather than decrease cell death, or affect both a 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 oflesions (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) areIn. Continuous lines (extrapolation) arc 6,000 X 100.43T (gee,) 1,000 X 100 77. (blue). 1010 M=10-5 102 0 10 20 30 40 50 60 T (months) 2 I NATURE I VOL 000 I 00 MONTH 2015 COM Macmillan Publishers Limited. All rights reserved EFTA00603738 LETTER RESEARCH s-0% b s-1% 7 r.0 .6 44444 -2 MAIM OS• b 0 d Net growth rate 0.56 0.34 days" Net growth rate 0.16 0.07 days', lor E 3. r.. is 2.5 100 ' 81" 'r ' 100 N• I d [I I I _ Ir. I. .6-2.0 X 80 00 s-1% 60 1 I .IAI I . 44.1.4 GAs 2c 40 S 13 $13 11 t' 40 , I IA .0 IA - 10. Ur m. o IA -104 !! A 2o Z 1. °: 20 laic!! M - 10 ' Cl` 20 I n ' IA .10 ' 00 0.6 1.0 1.5 2.0 0 1st.. s-0% s-1% 0 2 4 6 8 10 12 10 20 30 40 Figure 4 Genetic diversity is strongly reduced by the emergence of driver T (month) T(month) mutations. a-f. For all, M = 0 and the initial net growth rate = 0.007 clays-1 Figure 3 I Treatment success rates depend on the net growth rate of (d = 0.996). The three most abundant genetic alterations (GAs) have been tumours. a, Time snapshots before and during therapy (M = Resistant colour-coded using red (R), green (G) and blue (8) (c). Each section is 80 cells subpopulations that cause the tumour to regrow after treatment can be seen thick. Combinations of the three basic colours correspond to cells having two or at T = I month. b, c, Probability of tumour regrowth (P,g,„,$) as a function three of these genetic alterations. a, No drivers—separated, conical sectors of time after treatment initiation, for different dispersal probabilities (M) and emerge in different parts of the lesion, each corresponding to a different clone. net growth rates of the resistant cells. A higher net growth rate (b) leads to b, Drivers with selective advantage s = 1% lead to clonal expansions and a high regrowth probability, so that 50% of tumours regrow 6 months after many cells have all three genetic alterations (white area). d Genetic diversity treatment is initiated when M = I 0-5. c., Tumours with lower net growth rates can be determined quantitatively by randomly sampling pairs ofcells separated require >20 months to achieve the same probability of regrowth. Number of by distance r and counting the number of shared genetic alterations. e, The samples = Ito 800 per point (282 on average). Error bars areM. Sec number of shared genetic alterations versus the normalized distance il<r> Methods for details. 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 cell birth and cell death (Extended Data Fig. 514; the most important samples = 50 per data point. Error bars arc M. parameter is the fitness gain, s, conferred by each driver mutation. Figure 4a shows that in the absence of any new driver mutations (as historically been viewed as a feature of cancer associated with late for a perfectly normal cell growing in utero), donal subpopulations events in tumorigenesis, such as invasion through basement mem- would be restricted to small, localized areas. Each of these areas has at branes or vascular walls, this classical view of migration pertains to least one new genetic alteration, but none of them confers a fitness the ability of cancer cells to migrate over large distance?'. Instead, our advantage (they are 'passengers'). In an early tumour, in which the analysis reveals that even small amounts of localized cellular move- centre cell contains the initiating driver gene mutation, the same struc- ment are able to markedly reshape a tumour. Moreover, we predict ture would be observed—as long as no new driver gene mutations have that the rate of tumour growth can be substantially altered by a change yet appeared. The occurrence of a new driver gene mutation, however, in dispersal rate of the cancer cells, even in the absence of any changes markedly alters the spatial distribution of cells. In particular, the het- in doubling times or net growth rates of the cells within the tumour. erogeneity observed in normal cells (Fig. 4a) is substantially reduced Some of our predictions could be experimentally tested using new cell (Fig. 4b and Supplementary Video 5). The degree of heterogeneity can labelling techniques'''. Our results could also greatly inform the be quantified by calculating the number of genetic alterations (passen- interpretation of mutations in genes whose main functions seem to gers plus drivers) shared between two cells separated by various dis- be related to the cytoskeleton or to cell adhesion rather than to cell tances (Fig. 4d-f). The genetic diversity is markedly decreased (Fig. 4e), birth, death, or differentiationn". For example, cells that have lost the even with relatively small fitness advantages (s = 1%). This also has expression of E-cadherin (a cell adhesion protein) are more migratory implications for the number of genetic alterations that will be present than normal cells with intact E-cadherin expression", and loss of in a macroscopic fraction (for example, >50%) of all cells. Figure 4f E-cadherin in pancreatic cancer has been associated with poorer pro- shows that this number is many times larger for s = I% than s = 0%. gnosis", in line with our predictions. Furthermore, our model predicts that virtually all cells within a large Online Content Methods along with any additional Extended Data display items tumour will have at least one new driver gene mutation after 5 years of and Source Data. are available in the online version of the paper: references unique growth (Extended Data Fig. 5a).The faster the clonal expansion occurs to these sections appear only in the online paper. (the larger s is), the smaller the number of passenger mutations (Extended Data Fig. 5d, e). Our results are also robust to changes to Received 1September 2014: accepted 23 July 2015. the model (Methods and Extended Data Figs 5 and 6). We stress that Published online 26 August 2015. an important prerequisite for limiting heterogeneity is cell turnover in 1. Vogelstein.B. eta?. Cancergenome landscapes. Science 339,1546-1558(2013). the tumour, because in the spatial setting cells with driver mutations 2. Yachida. S. et aL Distant metastasis mass late during the genetic evolution of can 'percolate through the tumour only if they replace other cells. In pancreatic cancer. Nature 467,1114-1117 (2010) the absence of cell turnover, tumours are much more heterogeneous 3. Sottoriva.A.et at Intiatumor heterogeneity in human glicblastoma reflects cancer evolutionanidynamia, Proc. NatI Acad. Sci USA 110.4009-4014 (2013) (Extended Data Fig. 6d). 4. Navin. N. eta?. Tumour evolution inferred by single-cell sequencing. NaMe 472, In summary, our model accounts for many facts observed clinically 90-94 (2011) and experimentally. Our results are robust and many assumptions can 5. Gerlinger. M. et at Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. N. Engl.I Med 366,883-892 (2012). be relaxed without qualitatively affecting the outcome (Methods and & Gatenby.R.A & Vincent.T.L.An evolutionary modelof carcinogenesis.CancerRes Supplementary Information). Although tumour cell migration has 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 24. Talmadge.J. E. & I.J.AACR Centennial Series: the !hom of cancer progression. Proc. NatI Acad. Sci USA 107.18545-18550 (2010). metastasis: historical perspective. Cancer Res. 70, 5649-5669 (2010). 9. Beerenwinkel. N.et Genetic progression and the waiting time to cancer. PLOS 25. Alcolea....etaL Differentiation imbalance in singJe oesophageal progenitor cells Camput aid 3, e225 (2007). causes clonal immortalization and field change. Nature Cell 8101.16.615-622 10. Durrett. R. & Moseley. S.Evolution of resistance and progression to disease during (2014). clonal expansion of cancer. Theor. Poput. Bid 77, 42-48 (2010). 2& Weber. K. eta?. RGB marking facilitates multicolorclonal cell tracking, Nature Med 11. Gonzalez-Garcia. L Sole. RV. & Costa. J. Metapopulation dynamics and spatial 17, 504-509 (2011). ----' • •-- ' 1A-c A - `4-c-• " cA cricA " cc" c/AcIn` 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 13. Martens. E. A. Kostadmov. R. Maley. C. C. & Hailatschek O. Spatial structure integrin-mediated adhesion and cell migration.Smaff GTPases 5.e27958(2014). increases the waiting time for cancer. Newt Phys. 13.115014 (2011} 29. Gall, T. M. H.& Frampton.A EGene of the month: E-cadherin (CDH1).J. Dia 14. Andersco.A. R. A. Weaver. A. M.. Cummings. P. T.& Quaranta. V. Tumor Pathot. 66.928-932 (2013). morphology and phenotypic evolution driven by selective pressure from the 30. Winter. J. M. et af. Absence of E-cadherin expression distinguishes noncohesive microenvironment. CM' 127, 905-915 (2006). from cohesive pancreatic cancer. Cfin. Cancer Res. 14, 412-418(2008). 15. Kim.Y..Magdalena.A.S&Othrner.H.G.A hybridmodelkr tumor.spheroid growth Supplementary Information is available in the online version of the paper. in vitro I: theoretical development and early results.Math.Models Methods Ap#.Sci. 17.1773-1798(2007). Acknowledgements Support from The John Templeton Foundation is gratefully 16. McDougall. S R., Anderson, A. R. & Chaplain. M.A Mathematical modeling of acknowledged B.W. was supported by the Leverhulme Trust Eery-Career Fellowship. dynamic adaptive tumour-induced angicgenesis: clinical implications and and the Royal Society of Edinburgh Personal Research Fellowship. 1.8. was supported therapeutic targeting strategies. J. Meer. Blosi. 241, 564-589 (2006). by FoundationalQuestions in Evolutionary Bio Grant 17. Hawkins-Deemer, A.. Rockne. R. C..Anderson, A. RA & Swanson. K FL Modeting S.M.acknowledge support from The Virginia and M. Ludwig Fund for Cancer tumor-associated edema in gliomas during anti-angiogenic therapy and its Research. The Lustgarten Foundation for Pancreatic Cancer Research. The Sd imnact nn imao4:444tornnr Pm& Chun( 2 SA ontm, 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 pr ams 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 ISOZIC. Liven. U.& NOWalt.IA A LlytlatOKS Ot targeted cancer tnerapy. MOOS I1101. contributions from= and R.H.H. Med.18. 311-316 (2012). 21. Bozk, L& Nowak M. A.Timirg and heterogeneity of mutations associated with Author Information Re .nts and perrnissions information is available at drug resistance in metastatic cancers.Proc. NatfAcad. Sci. USA 111. The authors declare no competing financial interests. 15964-15968 (2014). Readers are welcome to comment on the online version of the paper. 22 Turke, A. B. et aL Preexistence and clonal selection of MET amplification in EGFR Cares ndence and requests for materials should be addressed to= mutant NSCLC. Cancer CO 17, 77-88 (2010). 4 I NATURE I VOL 000 100 MONTH 4015 INOIS Macmillan Publishers Limited. All rights reserved EFTA00603740 LETTER METHODS bacterial colonies"' show that the structure of the lattice (or the lack thereof in No statistical methods were used to predetermine sample size. Experiments were off-lattice models) has a marginal effect on genetic heterogeneity. not randomized and investigators were not blinded to allocation during experi- Asynchronous cell division. Division times ofrelated cells remain correlated for a ments and outcome assessment few generations. However, stochastic cell division implemented in our model is a Spatial model for tumour evolution. Tumour modelling has a long tradition". good approximation for a large mass of cells and is much less computationally Many models of spatially expanding tumours were proposed in the past"'"', expensive than modelling a full cell cycle. but they either assume very fewn-m-s'-"."."-" or no new mutations at allumat".". Replication faster at the boundary than in the Interior. Several studies have or one- or two-dimensional growthwwn'". On the other hand, well-mixed described a higher proliferation rate at the leading edge of tumours, and this has models with several mutations'" do not often include space, and computa- been associated with a more aggressive clinical course". To estimate the range of tional models aimed at being more biologically realistic"' require too much values of death rate d for our model. we used the proliferation marker Ki67. computing resources (time and memory) to simulate realistically large tumours Representative formalin-fixed, paraffin-embedded tissue blocks were selected (N., 10' cells). Our model builds on theEden latticemodel' and combines spatial from four small chromophobe renal cell carcinomas and six small hepatocellular growth and accumulation ofmultiple mutations. Since we focus on the interplay of carcinomas by the pathologist (=.). A section of each block was immunola- genetics. spatial expansion and short-range dispersal of cells. for simplicity we do belled for IC67 using the Ventana Benchmark XT system. Around 8-12 images. not explicitly model metabolism". tissue mechanics, spatial heterogeneity of tis- depending on the size of the lesion, were acquired from each tumour. Fields were sues, different types of cells present or angiogenesis". chosen at random from the leading edge and the middle of the tumour and were A tumour is made ofnon-overlapping balls (microlesions) ofcells. Tumour cells not necessarily 'hot spots' of proliferative activity. Using an Image) macro, each occupy sites of a regular 3D square lattice(Moore neighbourhood, 26 neighbours). Ki67-positive tumour nucleus was labelled green by the pathologist. and each Empty lattice sites are assumed to be either normal cells or filled with extracellular Ki67-negative tumour nucleus was labelled red. Other cell types (endothelium, matrix and are not modelled explicitly. Each cell in the model is described by its fibroblasts and inflammatory cells) were not labelled. The proliferation rate was position and a list of genetic alterations that have occurred since the initial neo- then calculated using previously descrthed methods". Statistical significance ofthe plastic cell, and the information about whether a given mutation is a passenger. results was determined using a Kolmogorov-Smimov nap-sample test (signifi- driver. or resistance-carrying mutation. A passenger mutation does not affect the cance level 0.05). The study was approved by the Institutional Review Board of the net growth rate whereas a driver mutation increases it by disrupting tight regu- Johns Hopkins University School of Medicine. In all ten tumours, the proliferation lation ofcellular divisions and shifts the balance towards increased proliferation or rate at the leading edge of the tumour was grmter than that at the centre by a factor decreased apoptosis. The changes can also be epigenetic and we do not distinguish of 1.25 to 6 (Extended Data Table I). Comparing the density of proliferating cells between different types of alterations. We assume that each genetic alteration to our model gives cf.., 0.56 (range: d = 0.176...0.86), which is what we assume in occurs only once ('infinite allele model"). The average numbers of all genetic the simulations of aggressive lesions. alterations, driver and resistant genetic alterations produced in a single replication Equal fitness of all cells in metastatic lesions. We assume that cells in a meta- event are denoted by y, 7d. and yr. respectively. When a cell replicates, each of the static lesion are already very fit since they contain multiple driven. Indeed. studies daughter cells receives n new genetic alterations of each type (n being generally ofprimary tumours and their matched metastases usually fail to find driver muta- different in both cells)drawn at random from the Poisson probabilitydistribution: 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. e-712(11.12)" Experimental evidence in microbes" and (to a lesser extent) in eukaryotes" sug- POO — (I) n! gests that fitness gains due to individual mutations are largest at the beginning of in which x denotes the type of genetic alteration. an evolutionary process and that the effects oflater mutations are much smaller.It In model A shown in Figs 2-1. replication occurs stochastically. with rate remains to be seen how well these results apply to late genetic alterations in proportional to the number of empty sites surrounding the replicating cell. and cancer" but if true. new driven occurring in the lesion are unlikely to spread death occurs with constant rate depending only on the number of drivers. We also through the population before the lesion reaches a clinically relevant size. simulated other scenarios (models B, C and D. see below). Driver mutations Dispersant,our model. cells detach from the lesion and attach again at a different increase the net growth rate (the difference between proliferation and death) either location in the tissue. This can be viewed either as cells migrating from one place to by increasing the birth rate or decreasing the death rate by a constant factor 1 + s. another one. or as a more generic mechanism that allows tumour cells to get better in which s> 0. access to nutrients by dispersing within the tissue. hence providing a growth Dispersal is modelled by moving an offspring cell to a nearby position where it advantage over cells that did not disperse. Some mechanisms that do not involve starts a new microksion (Extended Data Fig. la). Klicroksions repel each other; a active motion (that is. cells becoming motile) are discussed below. 'shoving algorithm's" (Extended Data Fig. Ib) ensures they do not merge. Migration. Cancer cells are known to undergo epithelial-to-mesenchymal trans- Code availability. The computer code (available at httpifwww2.ph.ed.ac.uki ition. the origin of which is thought to be epigenetic". This involves a cell becom- —bwadawfcancer-code) can handle up to 1 X 109 cells. which corresponds to ingmotik and miming some distance. If the cell finds the right environment,it can tumours that are clinically meaningful and can be observed by conventional switch back to the non-motile phenotype and start a new lesion. Motility can be medical imaging (diameter >1cm). The algorithm is discussed in details in the enhanced by tissue fluidization due to replication and death". Instead of mod- Supplementary Information. It is not an exact kinetic Monte Carlo algorithm elling the entire cycle (epithelial-mesenchymal-epithelial). we only model the because such an algorithm would be too slow to simulate large tumours. A com- final outcome (a cell has moved some distance). parison with kinetic Monte Carlo for smaller tumours (Supplementary Tumour buds. Many tumours exhibit focally invasive cell clusters, also known as Information) shows that both algorithms produce consistent results. tumour buds. Their proliferation rate is less than that of cells in the main turnout'. Model parameters. The initial birth rate 6= In(2) 0.69 days "'. which corre- We propose that tumour buds contain cells that have not yet completed epithdial- sponds to a 21h minimum doubling time. The initial death rate d = 0...0.995b to.meserichymal transition and therefore they proliferate slower. depends on the aggressiveness of the tumour (larger values = less aggressive lesion). Single versus cluster migration. Ref. 82 found that circulating cancer cells can In simulations of targeted therapy. we assume that. before treatment, = 0.69 travel in clusters of 2-50 cells, and that such clusters can initiate metastatic foci. days"' and d= 0.56 = 0.35 dayss '. whereas during treatment b= 035 days- ' They report that approximately one-halfof the metastatic foci they examined were and d = 0.69 days- ', that is, birth and death rates swap places This rather arbitrary initiated by single circulating cancer cells, and that circulating cancer cell dusters choice leads to the regrowth time of about 6 months which agrees well with clinical initiated the other half. The authors also note that the cells forming a duster are evidence. Mutation probabilities are 7 = 0.02,7d = 1 X 10-5.7, = I X 10-7. in line probably neighbouring cancer cells from the primary tumour. This means that the with experimental evidence and theoretical work"'". Since there are no reliable genetic make-up of cells within a newly established lesion will be very similar. data on the dispersal probability M. we have explored a range of values between regardless of its origin (single cell versus a small cluster of cells). Therefore, the M = 1 X 10-7 and I X 10-1. An parameters are summarized in Extended Data Fig. ability to travel in clusters should not affect the genetic heterogeneity or regrowth Ic, see also further discussion in Supplementary Information. probability as compared to single-cell dispersal from our modeL Validity of the assumptions of the modeL Our model is deliberately oversim- Angiogenesis. We do not explicitly model angiogenesis for two reasons. First. plified. However, many of the assumptions we make can be experimentally jus- most genetic alterations that can either change the growth rate or be detected tified or shown not to qualitatively affect the model. experimentally must occur at early stages of tumour growth as explained before. Three-dimensional regular lattice of cells. The 3D Moore neighbourhood was Hence, the genetic make-up of the tumour is determined primarily by what chosen because it is computationally fast and introduces relatively fever artefacts happens before angiogenesis. Second, local dispersal from the model mimics related to lattice symmetries. Real tissues are much less regular and the number of tumour cells interspersing with the vascularized tissue and getting better access nearest neighbours is different". However. recent simulations of similar models of 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 regression, and some do not. If only resistant cells can migrate, regrowth is faster cells and the extracellular matrix. We do not explicitly include biomechanics (see, (Extended Data Fig. 4b, c). Extended Data Fig. 4d-g shows regrowth probabilities however, below).in contrast to more realistic models"". as this would not allow us P,,,,,„„th for different treatment scenarios not mentioned in the main text. depend- to simulate lesions larger than about 1 X 106 cells. Instead, we take experimentally ing on whether the drug is cytostatic (bite„,tr,„„, = 0) o cytocidal = b). determined values for birth and death rates. values that are affected by biomecha- and whether d = 0 or d>0 before treatment. In Extended Data Fig. 4d, cells nics. as the parameters of our model. replicate and die only on the surface. and the core is 'quiescent'—cells are still Isolated balls of cells. In our simulations. balls of cells are thought to be separated alive there but cannot replicate unless outer layers are removed by treatment by normal, vascularized tissue which delivers nutrients to the tumour. The envir- (Supplementary Videos 6 and 7). Pitp.0.0 does not depend on the dispersal prob- onment of each ball is the same. and there are no interactions between the balls ability Af at all. and is close to 100% for N> 108 cells. a size that is larger than for other than mechanical repulsion. This represents a convenient mathematical con- d > 0 (Extended Data Fig. 40. It can be shown that P,,,,„,„h = 1 — exp(-7,4 trivance and qualitatively recapitulates what is observed in stained sections of Extended Data Fig. 4e is for the cytostatic drug (Oh...um= = 0): this is actual tumours (Fig. la). We investigated under which circumstances the balls also equivalent to the eytocidal dnigif the tumour has a necrotic core (cells are dead of cancer cells would mechanically repel each other. see Extended Data Fig. 7 for a but still occupy physical volume). In this case, increases with Af because graphical summary of the results. We simulated a biomechanical, off-lattice model more resistant cells are on the surface for larger M (cells can replicate only on the of normal tissue composed of 'ducts' lined with epithelial cells and separated by surface in this scenario). Extended Data Fig. 4f. g shows models with cell death stroma (Supplementary Information. section 8). Mechanical interactions between present even in the absence of treatment (d = 0.9b) but occurring only at the cells were modelled using an approach similar to that described previously'". "s surface. unlike in Fig. 3 where cells also die inside the tumour. Death increases with model parameters taken from refs 59, 60.85-88. We assumed cancer cells to owing to a larger number of cellular division necessary to obtain the same be of epithelial origin. as are most cancers". Cancer cells that invaded different size and hence more opportunities to mutate. areas ofepithelium grew into balls that remained separated by thin slices of stroma Relaxing the assumptions of the model. Figure 4 shows that even a small fitness (Supplementary Videos 8-11). This 'encapsulation' of tumour microlesions was advantage substantially reduces genetic diversity through the process of donal possible owing to the supportive nature of stroma that is able to mechanically resist expansion, see also Supplementary Videos 4 and S. We now demonstrate that this expansion of balls of cancer cells. Encapsulation is essential if the balls are to repel also applies to modified versions of the model, proving its robustness. each other. If the tissue is 'fluidized' by random replication and death. the balk Exact values ofMend s has no qualltathv effect. Extended Data Fig. 8b, e shows quickly merge (Supplementary Video 12). Another important factor are differ- that the average number of shared genetic alterations is larger in the presence of ences in mechanical properties of tumour and normal cells": it is known that drivers also in the case of non-zero dispersal (elf >0), and its numerical value is differences in cdlular adhesion and stiffness promote segregation ofdifferent types almost the same as for M = 0 (Fig.4). Extended Data Fig. 8c. f shows that as long as of cells",. s> 0 and regardless of its exact value, driver mutations reduce genetic diversity in In reality. mkrolesions within the primary tumour are less symmetric and some the tumour compared to the case s = 0. Extended Data Fig. Sa-c shows how many of them are better described as 'protrusions bulging out from the main tumour driver mutations are expected to be present in a randomly chosen cell from a tissue, owing to biomechanical instabilities: see. for example, refs 93.94. However, tumour that is Tyears old. Neither dispersal nor the way drivers affect growth (via stroma may still provide enough spatial separation. and the capillary network of birth or death rate) has a significant effect on the number of drivers per cell blood vessels—either due to tumour angiogenesis or preexisting in the invaded (Extended Data Fig. 56. c). A small discrepancy visible in Extended Data Fig. 5b tissue—may provide enough nutrients to the lesions so that our assumption of is caused by a slightly asymmetric way death and birth is treated in our model, see independently growing balls of cells remains valid. Therefore, we believe that the Supplementary Information. modelling the tumour as a collection of non- or weakly-interacting microksions Model B. Cells replicate with constant rate if there is at least one empty neighbour. is essentially correct. We also note that the existence of isolated balls is not neces- In the absence of drivers, genetic alterations are distributed evenly throughout the sary to explain our qualitative results: reduced heterogeneity and increased growth lesion (Extended Data Fig. 6b) but they often occur independently and the number in the presence ofmigration. Supplementary Video 13 shows that even if the tissue of frequent genetic alterations is low (Extended Data Fig. 6e). Drivers cause clonal is homogeneous and highly dynamical and there are no isolated balls of cells, expansion as in model A. migration leads to a considerable speedup of growth as compared to the case with Model C. Cells replicate regardless of whether there are empty sites surrounding no migration (Supplementary Video 14). them or not. When a cell replicates, it pushes away other cells towards the surface Tumour geometry and heterogeneity in the absence of driver mutations. (Supplementary Information). Extended Data Fig. bc,e shows that this again leads Supplementary Videos 2 and 3 illustrate the process of growth of a tumour with to clonal expansion which decreases diversity. maximally N= 107 cells, for Af = 0 and Af = l0"". respectively, and for d = 0.5. Model D. Replication/death occurs only on the surface and the core of the tumour Extended Data Fig. 2 shows snapshots from a single simulation for M = 0, is static Extended Data Fig. 6d shows that driver mutations cannot spread to the N— 10'. and d = 0 (no death. Extended Data Fig. 2a) and d = 0.9 (Extended other side of the lesion and conical clonal sectors can be seen even for s> 0. The Data Fig. 2b). In the latter case. cells are separated by empty sites (normal cells/ number of frequent genetic alterations is the same fors = 0 and s = 1%, indicating extracellular matrix). Extended Data Fig. 2c shows that the tumour is almost that genetic heterogeneity is not lowered by clonal expansion. This demonstrates spherically symmetric for Af = 0. The symmetry is lost for small but non-zero that cell turnover inside the tumour is very important for reducing heterogeneity. Af, and restored for larger M when the balk become smaller and their number To obtain the same (low) heterogeneity as for models a-c, the probability of driver increases. Extended Data Fig. 2c also shows that metastatic tumours contain many mutations must be much larger in model D (Extended Data Fig. 60. donal sectors with passenger mutations. Extended Data Fig. 8a shows that the Drivers affecting M. We investigated three scenarios in which drivers affect fraction G(r) of genetic alterations that are the same in two randomly sampled cells (1) only the dispersal probability Al—, (1 + q)M. in which q> 0 is the 'migration (Fig. 4) separated by distance r quickly decreases with r. indicating increased genetic heterogeneity owing to passenger mutations. 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. Targeted therapy of metastatic lesions. Models of cancer treatments" "° often Extended Data Fig. 3c shows that growth is unaffected in cases (1, 3) compared assume either no spatial structure or do not model the emergence ofresistance. We to the neutral case. For (2) the tumour growth rate increases significantly when the 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 tumour is larger than N= I X 106 cells. This shows that migration increases the same birth and death rates. but after the treatment resistant cells continue to 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 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 constant (motility) of organisms. death rate of sensitive cells during treatment is greater than the birth rate. or the Six-site (von Neumann) neighbourhood. We simulated a model in which each tumours would not be sensitive to the drug. For example. in Fig. 3 treatment cell has only six neighbours (von Neumann neighbourhood) instead of 26(Moore increases the death rate and decreases the growth rate of susceptible cells, the neighbourhood). Extended Data Fig. 9 compares models A and C for the two growth rate of resistant cells after therapy is identical to that of the sensitive cells neighbourhoods and show that there is only a small quantitative difference in the before treatment. d = 0.5b in the absence of treatment. elf = I0 6, and treatment growth curves for model A (model C is unaffected), but that the shape of the ball of begins when the tumour has N= 10' cells. cells deviates more from the spherical one for the six-site neighbourhood. see also Note that our model assumes the drug is uniformly distributed in the tumour": section 7 in the Supplementary Information. it is known that drug gradients can speed up the onset of resistance'. 31. Supplementary Video 1 and Extended Data Fig. 4a show that. since the process of resistance acquisition is stochastic. some tumours regrow after an initial 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 68. Diaz.L A.JretaL The molecularevolution ofacquireclresistance to targeted EGFR imnartanra nt roll aelhoninn Math MM iaint 22 1F1-1PA rinnci blockade in colorectal cancers.Natare 486.537-540 (2012). n. 69. Honda. H.. MoritA T.& Tanabe. A. Establishment of epidermal cell columns in mammalian skin: computer simulation.). Them. Biol. 81, 745-759 (1979). Lavrentoem. M. U. as Nelson U. K. survival prObablatieS at sonencaiirontiers 70. Ali.A. Somfai. E.& Grosskinsky. S. Reproduction-time statistics and segregation patterns in growing populations. Pros. Rev. £ 85, 021923 (2012). 71. Korolev. IC S.. Xavier. J. B.. Nelson. D. R. & Foster, K. R. Data from: a quantitative test of population genetics using a:eta-genetic patterns in bacterial colonies. Dryad Digital Repository. http://dxdolorg/10.5061/dryad.3147q (2011). 72. Gong P..Wang.Y.. Liu, G..Zhang J.&Wang Z.New insight intoKi67expression at the invasive front in breast cancer.PLoS ONES, e54912 (2013} 73. Ellison. T. kat at A single institution's 26-year experience with nonfunctional spatialy structured tissue.). Math. Blot http://dxdoicrg/10.1007/s00285. pancreatic neuroendocrinetumas:a validation of current stalingsystemsand a 015-0912-1(2015). new prognostic nomogram.Mn Su% 259, 204-212 (2014). 38. Gerlee. P.& Nelander. S. The impact of phenotypic switching on glioblastoma 74. Jones. S, Oat Comparative lesion sequencing provides insights into tumor growth and invasion. PLOS Comput Biot. 8, e1002556 (2012). evolution. Proc. Natl Aced So. USA 105,4283-4288 (2008). 39. Gonzalez-Garda,I, Solo. R V.& Costa.). Metapopulation dynamics and spatial 75. Lindstrom. L S. eta/. Ctinically used breast cancer markers such as estrogen heterogeneity in cancer. Prat Nan Acad. Sot USA 99,13085-13089 (2002). receptor, progesterone receptor. and human epidermal growth factor receptor 2 40. Sehyo.C.C.efaL Model forin vivo progression of tumors based on co-evolving cell are unstable throughout tumor progression.). Ctn. Oncot. 30, 2601-2608 population and vasculature. Sd. Rep. 1, 31(2011). (2012). 41. Torquato.S. Toward an 'sing model of cancer and beyond. Piro. Bid. 8, 015017 76. Voss. M. H. et at Tumor genetic analyses of patients with metastatic renal cell (2011). carcinoma and extended benefit from mTOR inhibitor therapy. C& CancevRes 42. Reiter.J. G.. Bozic. LAllen. B, Chattelee, K. & Nowak. M.A. The effect of one 20, 1955-1964 (2014). additional driver mutation on tumor progression. Eva 6, 34-45 (2013). 77. Wiser. M.J, Ribeck N. & Lensld, R. E Long-term dynamics of adaptation in 43. Kansal,A.R., Torquato.S.Harsh.G.R.Chixca.E A & Deisboecit. T. S. Simulated asexual populations. Science 342, 1364-1367 (2013). brain tumor growth dynamics using a three-dimensional cellular automaton. 78. White. T. C. Increased mRNA Levels of ERGI6. CDR, and MORI correlate with J. meat. Biol. 203, 367-382 (2000). increases in azole resistance in Candda afbicans isolates from a patient infected 44. Kemal. A. R. Torquato.S..Chiacca. EA & Deisboeck T. S. Emergence of a with human immunodeficiency virus. Antimicrob. Agents Chemother. 41. subpopulation in a corn putational model of tumor growth.). Thew. Biol. 207, 1482-1487 (1997). 431-441(2000). 79. McGranaten, N. etaL Clonal status of actionable driver events and the timing of 45. Antal.T. Knpivsky, P. L. & Nowak M.A. Spatial evolution of tumors with mutational processes in cancer evolution. So. Trend. Med 7, 283 (2015). successive driver mutations. Phys. Rev. £92.022705 (2015} 80. Ranft,J.etaL Fluidization of tissues bycell division and apoptosis.Proc NattAcad. 46. Endeding H, HLatky. L & tlahnfeldt. P. Migration rules: tumours are Sd. USA 107, 20863-20868 (2010). conglomerates of self-metastases. Br. J. Cancer 100, 1917-1925 (2009). 81. LeBleu, V. S. et al. PGC-12 mediates mitochondrial biogenesis and oxidative 47. Sottoriva etaL Cancer stem cell tumor modelreveals invasive morphology and phosphorylation in cancer cells to promote metastasis. Nature Cell Biol. 16. increased phenotypical heterogeneity. Cancer Res. 70. 46-56 (2010). 992-1003(2014} 48. Schaller. G. & Meyer-Hermann. M. Multicellular tumor spheroid in an off-lattice 82. Aceto, N. eta,. Circulating tumor cell dusters are oligodonal precursors of breast voronai-delaunay cell modelPhys. Rev. E71, 051910 (2005} cancer metastasis. Cell 158,1110-1122 (2014). 49. Radszuweit. M. Block. M.Hengstler. J.G.Sctibll, E & Drasdo, D. Comparing the 83. Sciume, G. eta?. A multiphase model for three-dimensional tumor growth. growth kinetics of cell populations in two and three dimensions.Phys. Rev. £79. Newt Phys. 15. 015005 (2013} 051907 (2009). 84. Charras,a&Sahai.E Physical influences of the extracellutarenvironment on cell 50. Moglia.R.Guisoni. N.& Albano, E.V.Interfacial properties in a discrete model for migration. Nature Rev. Mot. CellBiol. 15, 813-824 (2014). tumor. growth.Phys, Rev. E87, 032713 (2013} 85. Jiao.Y.& Torquato. S. Diversity of dynamics and morphologies of invasive solid 51. Foo. J. Leder, K. & Ryser. M. Multifocality and recurrence risk: a quantitative tumors. AIPAdvances 2.011003 (2012). model of field cancerization. J. moor. Biol. 355, 170-184 (2014). 52. PoleszczukJ.Hahnfeldt,P.&Enderling H.Evolution and phenotypic selection of 8,6. Galle.J.. LoeMer.M.& Drasdo, D. Modeling the effect of deregulated proliferation and apoptosis on the growth dynamics of epithelial cell populations in wIto. cancer stem cells. PLOS Cornput Bid. 11, e1004025 (2015). Biophys 188, 62-75 (2005} 53. Durrett. R. Schmidt. D.& Schweinsberg. J. A waiting time problem arisirg from the study of multi-stage carcinogenesisAnn.Appt hob& 19, 676-718 (2009). 87. Chen. E. J.. Novakofski. J.. Jenkins. W. K. & O'Brien. W. D. Jr. Young's modulus measurements of soft tissues with application to elasticity imaging. Ultrasonics. 54. Spencer, S. L.Eletryman. M. J.. Garda, J. A. &Abbott, D.M ordinary differential ferroelectrics and frequency control. IEEE Transactions 43.191-194 (199e. equation model for the multistep transformation to cancer.). Them Blot 231. 515-524(2004). 88. Samani, Bishop.), Luginbuhl. C.& Pieties. D. B. Measuring the elastic 55. Kim, Y.& Othrner. H. G.A hybrid model of tumor-stromal interactions in breast modulus of ex vivo small tissue samples/Wry& Med. Biol. At 2183 (2033). 89. Weinbe .R.A The ' of Cancer (Garland Science. 2007). cancer. &AL Math. Blot 75, 1304-1350 (2013). 56. Ramis-Conde,I, Chaplain. M.A. J., Anderson, A. RA & Drasdo. D. Multi-scale 90. Lekka. cancerous humanbladdercellsstudied by modelling of cancer cell intravasation: the role of cad herins in metastasis. Phys. scanning force microscopy. Ew. Biophys.) 28. 312-316 (1999} Blot 6.016008 (2009} 91. Gonzalez-Rodriguez. D., Guevoikian, K, Douezan. S.& Brochard-WyarL F. Soft 57. Swanson. K. R. erg Quantifying the role of angicgenesis in malignant matter models oideveioping tissues and tumors.Science 338.910-917(2012). progression of gliomas: in silica modeling integrates imaging and histology. 92. Stirba T.V. eta?. Fine tuning of tissues' viscosity and surface tension through Cancer Res. 71.7366-7375 (2011). contractility suggests a new role for l-catenin. PLoS ONE 8. e52554 (2013). 58. Taloni.A. eta?. Mechanical properties of growing melanocytic nevi and the 93. Drasdo. D. Bucking instabilities of one-layered growing tissues. Phys. Rev. Lett progression to melanoma. PLoS ONE 9, e94229 (2014). 84, 4244-4247 (2000). 59. Drasdo. D. & Hbhme, S.A single-cell-based model of tumor growth in vitro: 94. Base n. M., Joanny.J.-F., Prost, J. & Risler. T. Undulation instability of epithelial monclayers and spheroids. Phys. Biol. 2.133 (2005). tissues. Phys. Rev. Left 106, 158101(2011). 60. Drasdo. D.. Hoehme. S.& Block. M. On the role of physics in the growth and 95. Bozic. I. Mat Evolutionary dynamics of cancer in response to targeted pattem formation of multi-cellular systems: what can we learn from individual- combination therapy. Bite 2, e00747 (2013} cell based models?). Stet Phys. 128, 287-345 (2007). 96. Stylianopoulce,T. & Jain. R K. Combining two strategies to improve perfusion 61. Jiang, V., Pjesivac-Grbovic.J.CantrellC & Freyer. J. P. A multiscale model for and drug delivery in solid tumors. Proc. Natl Acad. Sci. USA 110,18632-18637 avascular tumor growth.Thophys. J. 89.3884-3894 (2005). (2013). 62. Eden. M. in A Two-Dimensional Growth Process (eels Family. F. & Woe* T.) 97. Foo.J.& Mictior. F. Evolution of acquired resistance to anti-cancer therapy. 265-283 (World Scientific.1961} J. Them. Blot 355,10-20 (2014} 63. Hard. D. L & Clark. A. G. Principles of Population Genetics (Sinauer 98. Goldie, J. H.& Coldman.A.J.The genetic origin of drug resistance in neoplasms: Associates. 1997). implications for systemic therapy. Cancer Res. 44.3643-3653 (1984} 64. Kreft J. U. Booth. G. & Wimpenny.J.W.T. BacSim. a simulator for individual- 99. Cokiman, A.J. & Goktia. J. H.A stochastk model for the origin and treatment of based modelling of bacterial colony growth. Microbiology 144, 3275-3287 tumors containing drug-resistant cells. BeN. Math Biel 48, 279-292 (1986). (1998). 100. Cokiman, A. J. & Goldie. J. H. A model for the resistance of tumor cells to cancer 65. Lardon, LA et at iCynoMICS: Next-generation individual-based modelling of chemotherapeutic agents. Math Mosel 65.291-307 (1983). biofilms. Enviton. Microbiel 13, 2416-2434 (2011} 101. Greulich, P.. Waclaw. B. &Allen. R J. Mutational pathway determines whether 66. Jones. &et at Comparative lesion sequencing provides insights into tumor drug gradients accelerate evolution of drug-resistant cells. Phys. Rev. Lett 109, evolution. Arc NallAcad. Sci. USA 105.4283-4288 (2008). 088101 (2012} 67. Wang. T. L ei af. Prevalence of somatic alterations in the colorectal cancer cell 102. Korolev. K S. eta?. Selective sweeps in growing microbial colonies. Phys. Biot. 9, genome.Prot Nate Acad. Sd. USA 99, 307E-3080 (2002). 026008 (2012} 4.2015 Macmillan Publishers Limited. All rights reserved EFTA00603743 RESEARCH LETTER a 2 1 ....... II- ........- . .......... Xi a X1 + overlap overlap b growth reduction reduction C Parameter Meaning Value References b Birth rate In 2 a 0.69 days., Thiswork d Death rate Thiswork Selective advantage 0 ... 0.1(0 ... 10%) This week 7 Mutation probability (all GAS) 0.02 C) 18.66.67.681 7e Mutation probab.11ty (driversonly) 4.10-50 (8.66.67.681 7, Mutation peobablety (resistance- 104 (8,66,67,68] carrying mutations) M Dispersal probability 0... 10-4 This work Extended Data Figure 1I Details of the model. a. A sketch showing how arc both moved apart along the line connecting their centres of mass (green dispersal is implemented: (1) A ball of cells of radius R„ in which the centre line) by as much as necessary to reduce the overlap to zero. The process is is at X,. is composed of tumour cells and normal cells (blue and empty squares repeated for all overlapping balls as many times as needed until there is no in the zoomed-in rectangle (2)). A cell at position x, with respect to the centre overlap. c. Summary of all parameters used in the model. If. for a given of the ball attempts to replicate (3). If replication is successful, the cell parameter, many different values have been used in different plots, a range of migrates with probability M and creates a new microlesion (4). The position X., values used is shown. Birth and death rates can also depend on the number of this new ball of cells is determined as the endpoint of the vector that starts of driver mutations, see Methods. Asterisk, parameter estimated from other at X, and has direction x, and length R,. b. Overlap reduction between the quantities available in the literature. balls ofcells. When a growingball begins to overlap with another ban (red), they (NOM Macmillan Publishers Limited. All rights reserved EFTA00603744 LETTER a ita wito time T b Otte. M=0 M=1O' Extended Data Figure 2 1Simulation snapshots. a, b, A few snapshots of tumour growth for no dispersal, and d = 0 (a) and d = 0.96 (b). To visua0ze clonal sectors, cells have been colour-coded by making the colour a M=106 M=104 arc non-cancer cells mixed with extracellular matrix. Note that images arc not to scale. c, Tumour shapes for N= I X 10', d = 0.96. and different dispersal probability M. Images not to scale, the tumour for M = I X 10-5 is larger heritabk trait and changing each of its RGB components by a small random than the one for M = 0. fraction whenever a cell mutates. The initial cell is grey. Empty space (white) .C.2015 Macmillan Publishers Limited. All rights reserved EFTA00603745 RESEARCH LETTER a b 1000 10' 800 600 400 200 1000 0 5 10 50 100 20 40 60 80 100 T [months) T [months] 10' N 10 - 8 2 1000 10 5 10 15 20 25 30 35 T [motel Extended Data Figure 31 Tumour size as a function of time. a, Growth of drivers increase both the net growth rate (s = 10%) and At (3. green) drivers a tumour without dispersal (M = 0), for d = 0.85. For large times (7), the either (with probability 1/2) increase Al tenfold (q = 9) or increase the net number of cells grows approximately as const X T 3. The tumour reaches size growth rate by s = 10%; (4. blue) drivers increase only the net growth rate by N = 1 X 10' cells (horizontal line) after about 100 months (8 years) of growth. s = 1094 and (5, black dashed line) neutral case with M = 1 X 10-7. which is b, The same data are plotted in the linear scale, with N replaced by 'linear indistinguishable from (1). In all cases (1-3) the initial dispersal probability extension' N1fd. c, Tumour size versus time when drivers affect the dispersal M= I X 10-7. PointsLeEesent average value over 40-100 simulations per data probability. In all cases, d= 0.95, and (1. black) drivers increase the dispersal point, error bars are =. rate tenfold (q = 9) but have no effect on the net growth rate (2, red) g..2015 Macmillan Publishers Limited. All rights reserved EFTA00603746 LETTER RESEARCH a b 10' 10' 10'' a 1000 d' woo 1000 100 10D 100 10 10 10 -150-100 -50 0 50 100 150 -150 -100 -50 0 50 100 -150-100 -50 0 50 100 150 T (days) T 'days) T (days) d e f 9 100 100 100 d=0.9b eo 80 80 60 • ' 60 60 I II g 40 '2 40 40 O 0: it 20 20 20 I I 0 •i I l 104 10' 10' 107 los la vs 106 10' ioe 10 70' los 106 101 10 N (cells) N (cells) N (cells) N Extended Data Figure 4 I Simulation of targeted therapy. a-c, The total The probability is plotted as a function of tumour size N just before the therapy number of cells in the tumour (black) and the number or resistant cells (red) commences. d, Before treatment, cells replicate only on the surface. Cells in versus time, during growth (T< 0) and treatment (T> 0), for —100 the core arequiescent and do not replicate. Therapy kills cells on the surface and independent simulations, for d = 0.56 for T <0. Therapy begins when cells in the core resume proliferation when liberated by treatment. e, As in = I X 106 cells. After treatment, many tumours die out (N decreases to zero) d, but drug is cytostatk and does not kill cells but inhibits their growth. The but those with resistant cells will regrow sooner or later. a.M = 0 for all cells at results for Pwp,,,,a, arc identical if the drug is cytotoxic and the tumour has a all times. b, M = 0 for all cells for T < 0 and M = 10-4 for resistant cells for necrotic core (cells die inside the tumour and cannot replicate even if the T> 0. c,M = 0 for non-resistant and M = 10-5 for resistant cells at all times.In surface is removed). f. Before treatment, cells replicate and die on the surface. all three cases, P II, is very similar. 36 t 5% (mean ± M.) (a), 25 ± 4% The core is quiescent Therapy kills cells on the surface (cytotoxic drug). (b),and 27 ± 6%for (c).d-g. Regrowth probability for four treatment scenarios g„ As in f, but therapy only inhibits growth (c)tostatic drug). In all cases not discussed in the main text. Data points correspond to three dispersal (d-g) error bars represent = from 8-1,000 simulations per point. probabilities: M = 0 (red), M = I X 10-5 (green), and M= I X 10-4 (blue). (g2015 Macmillan Publishers Limited. All rights reserved EFTA00603747 RESEARCH LETTER a b 2O 4 2.5- drivers I; .5=1% e if I ll 4 2.0 affect d, '§ 1.5 • 13 I 11 I drivers affect b, ,Iird k 1.5- k 1.0 if i II 1% ka I I li s 2% aw I 1.o II il tP5% e 0.5 I I II 11 s=1% I. 0.0 kilt.i. aerie 2 4 6 8 2 4 6 8 10 0 2 4 6 8 T (years] T (years) T (years] d e 40 40 .8 3° 130 20 a20 10 if: 10 0 0 0 1 2 3 4 0 1 2 3 4 drivers per cell drivers per eel Extended Data Figure 5 I Accumulation of driver and passenger genetic they affect death or birth rate. c, Dispersal does not affect the rate of driver alterations. a-c, The number of drivers per cell in the primary tumour lotted accumulation. d e, The number of passenger mutations (PMs) per cell versus as a function of time (10-100 simulations per point, error bars denote =.). the number of driver mutations per cell. More passenger mutations arc present a, M = 0 and three different driver selective advantages. For s = 1%, cells for smaller driver selective advantage (d). and this is independent of the accumulate on average one driver mutation within 5 years. The time can be dispersal probability M (e) in the regime of small M. Data points correspond to significantly lower for very strong drivers (s> 1%). b,The rateat which drivers independent simulations. accumulate depends mainly on their selective advantage and not on whether O2015 Macmillan Publishers Limited. All rights reserved EFTA00603748 LETTil= a C d Roclicaban terrpty space Repliczem'O it spaces0 Conslent cep .:.v.pon tale Suglace proAthrdearn e • f A r 0.5 a b c d Extended Data Figure 61 Genetic diversity in a single lesion for different ('necrotic core') is static. In all cases, N = I X 10', d = 0.996. e. Number of models. a-d. Representative simulation snapshots, with genetic alterations genetic alterations present in at least 50% of cells for identical parameters as colour-coded as in Fig. 4. Tops = O. bottoms = I%.a, Model A from the main in a-d. In all cases except surface growth (d). drivers inaease genetic text in which cells replicate with rates proportional to the number of empty homogeneity, as measured by the number of most frequent genetic alterations. nearby sites. b, Model B, the replication rate is constant and non-zero if there is Results averaged over 50-100 simulations error bars denote-. f, Model D. at least one empty site nearby. and zero otherwise. c, Model C. cells replicate at a with yd = 2 X I0-4 instead of 4 X 10-5, that is. drivers occur five times constant rate and push away other cells to make space for their progeny. more often. In this case, driver mutations arise earlier than in d. and the d. Model D. cells replicatehlk only on the surface, the interior of the tumour tumour becomes more homogeneous. O2015 Macmillan Publishers Limited. AEI rights reserved EFTA00603749 LETTER a b Parameter Meaning Value References Length of the simulation box 400-800 µm This work Cell diameter 10 um() [59.60.861 9 Elongation speed 0.208 [p.m/h1(1) This work Dynamic viscosity of the intracellular fluid 2.108 Pas [86) Ese Young modulus of epithelialkancer cells 1 kPa () [86) Es Young modulus of stroma 1 kPa (4) [87, 881 Poisson's ratio 1/3 () [69,60,861 Cell-cell adhesion energy 200 µl/ma () [861 d e . • • • .• t • •.• • • • • 4 c . f. I ' r • • . • p• , T.3600 days • •- • • , P. • I 4 • , . ••• • • n • $ A.C>0 4•••?... - • I ' i 1. S,. i. Extended Data Figure 71 The off-lattice model. a, Summary of all two nearby ducts repel each other as they grow as a consequence ofmechanical parameters used in the modeL Asterisk typical value, varies between different forces exerted on each other. d, The balls coalesce if growth is able to break types of tissues; dagger symbol, equivalent to 24 h minimal doubling time; the separating extracellular matrix. e, If the balls arc not encapsulated, double dagger symbol, based on the assumption that macroscopic elastic they quickly merge. 1, Isolated balls of cells are not required to speed up growth; properties of tissues such as liver, pancreases or mammary glands are primarily migration (left) can cause the tumour to expand much faster even if determined by the elastic properties of stroma. b, Simulation snapshot of a individual microlesions merge together. as opposed to the case with no normal tissue before the invasion of cancer cells. c, Two balls of cancer cells in migration (right). O2015 Macmillan Publishers Limited. All rights reserved EFTA00603750 LETTER RESEARCH a d —1 0.4 v.- 12 0.3 Po rig 0.2 2 g 0.1 tic'> b 100 M=10-7 3.51 6 3.0 40 s=1% s=0%. d=0.99b 1 g 2.5 2. 2.0 g 1.5 g 1.0 20 0.5 at. 0.0 0.5 1.0 1.5 2.0 2.5 30 rkr> C f 3.0 g 0 2.5 2.0 1.5 c 1.0 0.5 0.5 1.0 1.5 2.0 recta Extended Data Figure 8I Genetic diversity quantified. a:rumours are much N= 1 X 107, and Al = I X 10-7.1n the presence of drivers. G(r) decays slower. more genetically heterogeneous in the absence of driver mutations (s = 0) (see indicating more homogeneous tumours. c, The exact value of the selective Fig. 4). The plot shows the fraction G(r) of genetic alterations (GAs) shared advantage of driver mutations is not important (all curves G(r) look the same, between the cells as function of their separation (distance r) in the tumour. The except for s = 0) as long as s> 0. d-f, Number of genetic alterations present fraction quickly decreases with increasing r. The distance in the figure is in at least 50% of cells for identical parameters as in a-c. correspondingly. normalized by the average distance <r> between any two cells in the tumour. Drivers substantially increase the level of genetic homogeneity. In all panels For a spherical tumour. <r> is approximately equal to half of the tumour the results have been averaged over 30-100 simulations, with error bars diameter. b, Fraction of shared genetic alterations for s = 1% and s = 0%. as= 4.2015 Macmillan Publishers Limited. All rights reserved EFTA00603751 RESEARCH LETTER a Model A. d=0.9b to a 10 C 100 100 10 5 10 T (days) 50 100 Model A. d=0 95b. s=514 d J1 tl 10 100 T [days) Model C. no death 1000 to. 10 100 1000 10 15 20 25 T [days) T [days) Extended Data Figure 91 Growth curves for the 26-nearest neighbours very little difference in the growth curves between the 6n and 26n models. A (26n, red curves) and the 6-nearest neighbours (6n, green curves) models. small asymmetry in the shape is caused by faster-growing cells with driver a, Model A (as in the main text), no death. The tumour grows about twice as mutations.d.ModetC(exponentialgrowth).Growth is thc same in both 6nand slow in the 6ta model Pictures show tumour snapshots for both models: 26n models, but the shape is more asphcric for the Si model. This is there is no visible difference in the shape. b, Model A, death d = 0.8b. The probably caused by shifting cells along the coordinate axes and not along the additional blue curve is for the 6n model, with modified replication probability shortest path to the surface when making space for new cells. All plots show to account for missing neighbours as explained in the Supplementary the mean (average over 50-100 simulations) and M. Information. c, Model A. with death d= 0.95b. and drivers s = 5%. There is 42015 Macmillan Publishers Limited. All rights reserved EFTA00603752 LE1011=3 Extended Data Table 1 I Experimental results for the percentage of proliferating cells in the centre versus the edge of solid tumours Edge Center Ratio center:edge p-value Case Tumor type Images KI % Total no. of Images Kr % Total no. of cells cells 1 Chromophobe 8 3.46 1013 3 1.11 2561 0.32 0.05 RCC 2 Chromophobe 5 3.07 508 3 1.05 938 0.34 0.28 RCC 3 Chromophobe 5 2.63 524 3 0.44 697 0.17 0.03 ROC 4 Chromophobe 7 1.58 581 2 0.53 958 0.34 0.17 RCC 5 HCC 7 17.14 892 2 9.74 1637 0.57 0.05 6 HCC 7 51.71 1079 4 32.84 2562 0.64 0.03 7 HCC 6 47.37 435 3 19.97 1397 0.42 0.09 6 HCC 7 19.02 895 4 13.78 1191 0.72 0.35 9 HCC 6 15.09 1074 3 11.98 1094 0.79 0.33 10 HCC 9 29.84 1305 2 20.87 2457 0.70 022 Summary 1.4 Chromophobe 25 2.69 11 0.81 0.30 0.00002 RCC 5-10 HCC 42 30.0 18 19.1 0.64 0.007 A reprteentalive section leech tunas was labelled to the wale/atm/ marter ise OM end mage or the temou et the leering edge end the centre were acquired es descroed (klatnOCIS) Psettestion is merseSlyincreesee et the larding ear arwl this is stalistically similar/ (Summar% gremegaar-Sm mem two-same* note. /35). The .101Oft reboot the number d preareretng cells in me realest the edge 40.50 fronts0.17-07% HCC.Ms/alocelluler carcinoma RCC. renelan artisan. O2015 Macmillan Publishers Limited. All rights reserved EFTA00603753