DS9 Document EFTA01201684

1 C Timing and heterogeneity of mutations associated with drug resistance in metastatic cancers Ivana Bozic'' and Martin A. Nowa ka'bi 'Program for Evolutionary Dynamics, Department of Mathematics, and °Department of Organismic and Evolutionary Biology, Harvard University, Cambridge, MA 02138 Edited by Herbert Levine, Rice University, Houston, 1X and approved October 8, 2010 (received for review June 28, 2010) Targeted therapies provide an exciting new approach to combat human cancer. The immediate effect is a dramatic reduction in disease burden, but in most cases, the tumor returns as a conse- quence of resistance. Various mechanisms for the evolution of resistance have been implicated, including mutation of target genes and activation of other drivers. There is increasing evidence that the reason for failure of many targeted treatments is a small preexisting subpopulation of resistant cells; however, little is known about the genetic composition of this resistant subpopu• lation. Using the novel approach of ordering the resistant sub• dones according to their time of appearance, here we describe the full spectrum of resistance mutations present in a metastatic lesion. We calculate the expected and median number of cells in each resistant subclone. Surprisingly, the ratio of the medians of successive resistant clones is independent of any parameter in our model; for example, the median of the second done divided by the median of the first is —1. We find that most radiographically detectable lesions harbor at least 10 resistant subclones. Our pre- dictions are in agreement with clinical data on the relative sizes of resistant subclones obtained from liquid biopsies of colorectal can• cer patients treated with epidermal growth factor receptor (EGFR) blockade. Our theory quantifies the genetic heterogeneity of re- sistance that exists before treatment and provides information to design treatment strategies that aim to control resistance. cancer I drug resistance I heterogeneity I mathematical biology A d resistance to treatment is a major impediment to successful eradication of cancer. Patients presenting with early-stage cancers can often be cured surgically, but patients with metastatic disease must be treated with systemic therapies (1). Traditional treatments such as chemotherapy and radiation that exploit the enhanced sensitivity of cancer cells to DNA damage have serious side effects and, although curative in some cases, often fail due to intrinsic or resistance acquired during treatment Targeted therapies, a new class of drugs, inhibit specific molecules implicated in tumor development and are typically less harmful to normal cells compared with chemotherapy and radiation (2-5). In the case of many targeted treatments, patients initially have a dramatic response (6, 7), only to be followed by a regrowth of most of their lesions several months later (8-10). Acquired re- sistance is often a consequence of genetic alterations (usually point mutations) in the drug target itself or in other genes (10-14). Recently, mathematical modeling and clinical data were used to show that acquired resistance to an epidermal growth factor receptor (EGFR) inhibitor panitumumab in metastatic co- lorectal cancer patients is a fair accompli, because typical detectable metastatic lesions are expected to contain hundreds of cells resistant to the drug before the start of treatment (10). These cells would then expand during treatment, repopulate the tumor, and cause treatment failure. Similar conclusions should hold for targeted treatments of other solid cancers (15). Suc- cessful treatment requires drugs that are effective against the preexisting resistant subpopulation and must take into account the (possible) heterogeneity of resistance mutations present in the patient's lesions. In this article we use mathematical modeling to investigate the heterogeneity of drug-resistant mu- tations in patients with metastatic cancers. First mathematical investigations of the evolution of resistance to cancer therapy were concerned with calculating the proba- bility that cells resistant to chemotherapy are present in a tumor of a certain size (16). Later studies expanded these results to include the effects of a fitness advantage or disadvantage provided by resistance mutations (17, 18), multiple mutations needed to achieve resistance to several drugs (15, 19-21), and density limitations caused by geometric constraints (22). These studies used generalizations of the famous Luria-Delbrfick model for accumulation of resistant cells in exponentially growing bac- terial populations (23). Probability distribution for the number of resistant cells in a population of a certain size in the fully sto- chastic formulation of the Luria—Delbrfick model was recently calculated in the large population size limit (24, 25). The focus of above studies was describing the total number of all resistant cells, rather than the composition of the resistant population (26). Results We model the growth of a metastatic lesion as a branching process (27) that starts from a single cell (the founder cell of the metastasis) that is sensitive to treatment. Sensitive cells divide with rate b and die with rate d. The net growth rate of sensitive cells is r=b —d. During division, one of the daughter cells receives a resistance mutation with probability u. Resistant mutations can be neutral in the absence of treatment, which means they have the same birth and death rates as sensitive cells, and we initially focus on this case. We also expand our theory to the more general case where resistant cells are nonneutral, which means they have birth and death rates bR and de. respectively. If Significance Metastatic dissemination to surgically inaccessible sites is the major cause of death in cancer patients. Targeted therapies, often initially effective against metastatic disease, invariably fail due to resistance. We use mathematical modeling to study heterogeneity of resistance to treatment and describe for the first time, to our knowledge, the entire ensemble of resistant subclones present in metastatic lesions. We show that radio- graphically detectable metastatic lesions harbor multiple re- sistant subclones of comparable size and compare our predictions to clinical data on resistance-associated mutations in colorectal cancer patients. Our model provides important information for the development of second-line treatments that aim to inhibit known resistance mutations. Author contributions: IS. designed research; I.B. and M.A.N. performed research; I.B. and M.A.N. analyzed data; and IS. and M.A.N. mote the paper. The authors declare no conflict of Interest. This article is a NOS Direct Submissico. 'To whom correspondence may be addressed. Email: lbozIcOmath.harvard.edv or martinnowakeharvard.edu. This article contains supporting information online at www.pnas.orgilookup/supplidoi:10. 1073/pnal.141207511INDCSUPplemantal. www.paes.orgegweovio.tozneuistetzozsiti PNAS Early Edition I 1 of 5 EFTA01201684 kg VA a c = (bR —dR)1(b —d)> 1, then resistance mutations are advanta- geous before treatment; if c < 1, they are deleterious. A resistant cell may appear in the population and be lost due to stochastic drift or it can establish a resistant subclone. We number the resistant subclones that survive stochastic drift by the order of appearance (Fig. IA). A reasonable assumption for the number of point mutations that can provide resistance to a tar- geted drug is on the order of 100 (10, 28). Thus, the different resistant subclones will typically contain different resistance mutations, especially if we only focus on the largest ones. We calculate the number and sizes of resistant subclones in a metastatic lesion containing Al cells. Typical radiographically detectable lesions are — I cm in diameter and contain —109 cells. The mutation rate, u, leading to resistance is the product of the point mutation rate p, which is on the order of —10-9 per base pair per cell division, and the number of point mutations that can confer resistance, which is —100. In our analysis we will Cumulative distribution 1 0.9 0.8 0.7 0.6 0.5 O4 0.3 0.2 0.1 0 1 10 100 1000 Clone size 10000 100000 Ng. 1. Evolution of resistance in a metastatic lesion. (A) As the lesion (green) grows from one cell to detectable size, new resistant subclones ap- pear. Some of them are lost to stochastic drift (yellow and pink), while others survive (purple, red and orange triangle). Instead of looking at the time of appearance of new clones, our approach takes into account the total size of the lesion when the resistance mutation first occurred. (8) Agreement be- tween computer simulations and formula (I) for the cumulative distribution function for the number of cells in the first four resistant clones. The first subclone contains 10 or fewer cells with probability 0.06, between 10 and 100 cells with probability 0.30, between 100 and 1000 cells with probability 0.07 and more than 1000 cells with probability 0.13. The second subclone contains more than 100 cells with probability 0.36. Parameters b= 0 25, d=0.181, M=109, u=42.104 . assume a large M and small u limit and mostly focus on the case when Mu >> I. Tumor sizes at which successful resistant mutants are pro- duced can be viewed as a Poisson process on [0.,M] with rate u (SI Tat) (10, 17). The number of successful mutant lineages is thus Poisson distributed with mean 2 =Mu. If Mk is the number of cancer cells in the lesion when the kth mutant appeared, which survived stochastic drift (Fig. 1A), then Mkti —Mk is exponen- tially distributed with mean 1/u. Therefore, we expect that the kth clone appeared when the total population size was Mk —klu and that roughly the size of the first clone is k times the size of the kth clone. The probability that exactly k clones are present in the population of size Al is Ake-A/k1 Counting new successful resistant clones in the order of ap- pearance, we calculate the probability distribution for the num- ber of cells in the kth resistant clone. In particular, if k 4Z Mu, the cumulative distribution function for the number of resistant cells in the kth clone simplifies to Fk(y) As I (mu MU dy . Ill The excellent agreement between Formula 1 and exact computer simulations of the stochastic process is shown in Fig. lB. The mean number of cells in the kth resistant clone is E(Yi ) Ps [bMulr][log(rIbu)— I] and E(Yk) s /Wu l[r(k — 1)] for k> 2. The median for the number of cells in the kth subclone is given by bMu Med(Yk) —r (211* —1). [2] Interestingly, the ratio of the means of the two subclones k and/ is (i — 1)/(k — I) for k. j> 1. The ratio of their medians is Med(Yk) 21/k — I MeA i'l 2 3 — I These ratios are independent of any parameters of the process. In particular, the ratio of the medians of the first and second clone is in —1, which implies that they have comparable size (same order of magnitude). Liquid biopsy data were used to obtain estimates for the birth and death rates of cells in metastatic lesions and the number of point mutations providing resistance to the EGFR inhibitor panitumumab in colorectal cancer (10). The resulting parameter values (b=0.25 and d = 0.181 per day, point mutation rate p =10-9 per base pair per replication, and 42 point mutations conferring resistance) can be used to calculate the mean and median sizes of the resistant subclones in a metastatic lesion containing M = 109 cells. The mean numbers of cells in the first, second, and third appearing resistant clone are E(Y1)x2237, E(Y2) PS 152, and E(Y3) PS76, respectively. However, the mean for Yi, the size of the first resistant clone, is heavily influenced by the realizations of the stochastic process in which the first resistance mutation appeared early and is not a good summary of the probability distribution for Yi. Namely, the realizations in which the number of cells in the first clone is greater than the mean (2,237) account for less than 7% of all cases. The median number of cells in the first resistant clone [Med(Yi )] for the above parameters is 152, whereas the medians for Y2 and Y3 are 63 and 40, respectively. In SI Tar, we calculate the probability distribution for the ratio of resistant clone sizes YI/Yk and show that it is also in- dependent of the parameters of the process. Even though the first appearing clone is expected to be the largest, followed by the second clone and so on, we show that this ordering is often I31 2 of 5 I www.poes.cregfCgilda/10.107342nat.1412075111 Bozic and Nowak EFTA01201685 1 violated. In 31% of lesions, the first successful subclone is smaller than the second one: on the other hand, in 24% of lesions the first subclone is at least 10 times larger than the second one. Fig. 2 shows different realizations of the stochastic process of evolution of resistance in metastatic lesions containing 109 and 109 cancer cells. The same parameters were used to generate all lesions. The size of each subclone is shown (in number of cells), and the subclones are ordered by their time of appearance. In lesion LI, the first three subclones are the largest, and each have around 100 cells. Lesion L.5 contains only two subclones, whereas L6 contains seven subclones, but none has more than 10 cells. In each lesion of total size 109 cells, there are more than 10 resistant subclones. In L7, the two largest subclones contain 1,500 and 460 cells. In IS, there are five subclones of about 100 cells. In Table 1, we show clinical data for the number of circulating tumor DNA (ctDNA) fragments harboring mutations in five genes associated with resistance to anti-EGFR treatment in 18 colorectal cancer patients who developed more than one muta- tion in those genes (29). These mutations were not detectable in patients' serum before therapy, but became detectable during the course of anti-EGFR treatment. The number of ctDNA frag- ments correlates with the number of tumor cells harboring that mutation: it was previously estimated (using the tumor burdens and pretreatment ctDNA levels measured in patients who had KRAS mutations in their tumors before therapy) that one mu- tant DNA fragment per milliliter of serum corresponds to 44 million mutant cells in the patient's tumor (10). Thus, the ratios of the resistant clone sizes can be obtained from the ratios of the numbers of ctDNA fragments harboring resistance-associated mutations. These data provide a unique opportunity to test our theory and compare the relative sizes of resistant clones inferred from the data with those predicted using our model. Assuming that resistance-associated mutations with higher ctDNA counts appeared before those with lower ctDNA counts, we find ex- cellent agreement between the data and our model predictions. For example, the median ratio of the sizes of the first two re- sistant clones inferred from clinical data (29) is 2.21, whereas our model predicts 2.51. The median ratio of the sizes of the first and third clones from clinical data are 4.3, and our model predicts 4.12 (Table 1). This comparison is parameter free, as we showed A Lesion size M-108 cells LI L2 in two • too to II 0 II moo too to 'coo Z • I 1 2 3 4 5 6 7 8 9 ID B Lesion size M109 cells L7 U I 1 1000 4 t000 (.... I 0 100 100 I/ 10 10, Z 1 1 2 3 4 3 6 7 8 9 10 1 2345678910 100 2 3 4 5 6 7 8 9 10 10 that the ratio of resistant clone sizes is independent of parameters. Our mathematical results describe the relative sizes of re- sistant clones ordered by age, whereas the experimental data in Table I are ordered by size, which serves as a proxy for age, because exact clonal age is unknown. We quantify the extent to which this difference in clonal ordering by size vs. age influences our statistics using exact computer simulations (Table 1). In the relevant parameter regime of large lesion size, At, and small mutation rate, a, with Mu >, 1, the results are largely in- dependent of parameters (median ratios of clone sizes vary by <10% for different parameter combinations). We show simula- tion results for median ratios of clone sizes when clones are ordered by size for typical parameter values (10). As we see in Table 1, the ordering of experimental data by size does not sig- nificantly change the results of our analysis. We can generalize our approach to the case when resistance mutations are not neutral, but provide a fitness effect already before treatment (formulas shown in Si Text). In Table 2, we compare the predicted medians for the first five resistant clones in a metastatic lesion containing M= 109 cells when resistance is deleterious, neutral, or advantageous. We see from Table 2 that even if resistant cells are only 10% as fit as sensitive cells, they will still be present in typical lesions. The average number of resistant cells produced until the lesion reaches size M is Mu/s. Here s= 1 —d/b is the survival probability of sensitive cells, which is the probability that the lineage of a single sensitive cell will not die out. For typical parameter values (i.e., those used in Table 2), the number of resistant cells produced by sensitive cells in a single lesion is — 150. Resistant cells that are 10% as fit as sensitive cells have a survival probability of 4%; so on average, six of them will form surviving clones. The effect that mutations can cause treatment failure, although they have high fitness cost is a consequence of the high number of resistant mutants pro- duced by billion(s) of sensitive cells in a lesion and the specific properties of the branching process, namely the independence of lineages. Discussion In this paper we describe the heterogeneity of mutations pro- viding resistance to cancer therapy that can be found in any one L3 woo too 10 VIII_ woo too to 1.3 moo too I0 L6 1 2 1 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 910 1 2 3 4 3 6 7 8 910 1 2 3 4 5 6 7 8 9 10 L9 LI0 1000 100 to io 1 2 3 4 5 6 7 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 9 10 LII 1000 100 10 L12 hid 2 3 4 3 6 7 8 9 IS Resistant clones (order of appearance) Fig. 2. Resistant subclones in metastatic lesions. Different realizations of the same stochastic process are shown in each panel. (A) Six lesions of size 10° and (8) six lesions of size le cells. The first ten resistant clones are shown, which survived until time of detection. They are ordered according to their time of appearance. Parameter values for all simulations: b =0.25, d=0.181, u=42.104. 17 e BR Bozic and Nowak PNAS Batty Edition I 3 of 5 EFTA01201686 Table 1. Comparison of predicted ratios of resistant clone sizes and ratios obtained from clinical data FA a 1 Patient Yi* Y2 112 114 YI/Y2 YI/Y3 YI/Y4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Median from patients Predicted median Predicted median (order by size) 168 90 129 120 82 80 948 120 28 15 114 40 6,760 4,940 220 30 848 374 61 25 244 83 429 400 394 13 308 265 130 13 28 13 131 45 250 173 1.87 1.08 30 1.03 2.73 104 100 7.9 9.12 9.48 1.87 2.85 4,100 3,900 1.37 1.65 1.73 7.33 135 133 2.27 6.28 6.38 2.44 57 2.94 4.28 100 1.07 4.29 4 30.31 98.5 208 139 1.16 1.48 2.22 10 2.15 12 11 2.91 10.92 11.91 58 3t 1.45 431 8.06 2.21 43 7.22 2.51 4.12 5.74 2.05 3.63 5.25 •Number of circulating tumor DNA (ctDNA) fragments per milliliter (VI to Y4) harboring different mutations associated with resistance to anti-EGFR agents in colorectal cancer patients treated with EGFR blockade (29). Ratio of resistant clone sizes is given by the ratio of the ctDNA counts for any two resistance-associated muta- tions. We assumed that mutations with higher ctDNA counts in the patient data appeared before mutations with smaller ctDNA counts. We also report predicted median ratios obtained from computer simulations when clones are ordered by size (parameters: b =0.25, d=0.181, M=109, u =42 x 10-9). metastatic lesion. Our results can be generalized to take into account all of the patient's lesions, assuming that they evolve according to the same branching proetts and that the number of lesions is much smaller than l/u. In that case, the probability distribution for the size of the kth appearing resistant clone in the patient's cancer is given by Formula 1 if we let M be the number of cancer cells in all of the patient's lesions. All our results generalize similarly. Although the mean and median clone sizes in our model de- pend on the parameters of the process, their ratios are generally parameter free. The universality of the clone ratio statistics fol- lows from the fact that the skeleton of our branching process, which includes only cells with infinite line of descent, can be approximated by a Yule (pure birth) process (30). It has been shown that in the limit of large lesion size M and small mutation rate u, the statistics of the relevant clones in a branching process with death remain approximately Yule (31). Similarly, it can be shown that in the Yule process, in the above limits, the mean size of the kth largest clone is —Mu/(k-1), and the ratio of the mean sizes of the kth and jth largest clones is —(f —1)/(k — 1) (31, 32). This formula is exactly the result we obtain for the ratio of mean clone sizes even though we order clones by age. A few recent investigations studied the dynamics of single clones resistant to therapy (28, 33). In one of the studies (33), the authors used a generalization of the Luria-Delbriick model in which sensitive cells grow deterministically and calculated the number of individual resistant clones and the probability distri- bution for the number of cells in a single resistant clone after time I. In another study (28), mathematical modeling along with in vitro growth rates of cells harboring 12 point mutations Table 2. Sizes of resistant clones when resistance Is deleterious, neutral, or advantageous c.(bR -dR)/(b-d) First done* Second clone Third done Fourth clone Fifth done 0.01 0 0 0 0 0 0.1 10 6 4 2 0.5 27 17 13 11 10 0.7 50 26 19 15 13 0.9 103 46 30 23 18 0.95 125 54 35 26 20 1 152 63 40 29 23 1.05 186 74 45 33 25 t.1 229 87 52 37 28 •Median number of cells in the first five successful resistant clones in a metastatic lesion with M=109 cells when resistant cells are less fit than sensitive cells (c < 1), neutral (c =1), and more fit than sensitive cells (c> 1). We fix the birth and death rate of sensitive cells, b= a 25 and d=0.181, and the death rate of resistant cells dR =d. We vary the relative fitness of resistant cells, G and let the birth rate of resistant cells be bit =dR +c(b—d). Mutation rate u =42 x 104. For c =0.1 we report simulation results, and for c> 0 1, we use Eq. S13; see SI Text for details. 4 of 5 I wwwonas.orgfcgildoi/10.10734ffias.14I2075111 Bozic and Nowak EFTA01201687 providing resistance to BCR-ABL (fusion of breakpoint cluster re- gion gene and Abelson murine leukemia viral oncogene homolog inhibitor imatinib were used to calculate the number of resistant clones and the expected number of resistant cells with a particular resistance mutation at the time of diagnosis of chronic myeloid leukemia. The authors found that at most one resistant clone is expected to be present, as the total number of CML stem cells at diagnosis is estimated to be approximately M— 100.000 cells and is much smaller than the billions of cells typically present in a single detectable lesion of a solid tumor. In this paper, we use a different mathematical technique and the novel approach of ordering the resistant clones according to their time of appearance, which allows us for the first time, to our knowledge, to describe the full spectrum of resistance mutations present in a lesion. Our study is challenging the conventional view of the evolution of resistance in cancer. For every therapy that is opposed by multiple potential resistance mutations, which is the case for every targeted drug developed thus far, we can expect multiple resistant clones of comparable size in every lesion. Our theory provides a precise quantification of the relative sizes of those resistant subclones. The heterogeneity of resistance muta- tions is further amplified when taking into account multiple 1. Vogelstein B. et al (2013) Cancer gencrne landscapes. Science 331(61271:1546-1593. 2. Sawyers C (2004) Targeted cancer therapy. Nature 432(7015)194-297. 3. Midor F. et al. (2005) Dynamics of chronic Roelold leukaemia. mature 43547046): 1267-127D 4. Gerber DE. Minna .ID (2010) ALK inhibition for non-small cell lung cancer. From dis- covery to therapy in record time. Cancer CO 18(6):548-55I. S. Komerova NI.. Wcdarx D (2013) Targeted Cancer Treatment In Seiko: Small Molecule Mhibitors and Oncolytic Viruses (Springer. New York). 6. Chapman PH, et al; BRIM-3 Study Group (2011) improved survival with vemurafenib in melanoma with BRAE V600E mutation. N Engl f Med 364(26):2507-2516. 7. Maensondo M. et al.; Nonh.Fast Japan Stud/Group0010) Gef it In lb Or ChemOtheraPY for nerrunatacell lung cancer with mutated EGFR. N Eng: Med 362(25):2380-2368. 8. Katayama R. et al. (2011) Therapeutic strategies to overcome aizotist resistance n non-small cell lung cancers harboring the fusion oncogene EMIA-ALK. Prot Nati AM Sci USA 108(1).753S-7SM. 9. Sosman lA, et at. (2012) Survival in BRAE V600-mutant advanced melanoma treated with vemurafenta. N Eng/ / Med 366(83:707-714 10. Dias LA, Jr, et al. (2012) The molecular evolution of acquired resistance to targeted EGFR blockade in colorectal cancers. Nature 466(7404)537-540. FaCi W, et at. (2005) Acquired resistance of lung adenocarcinomas to gelitinib or er- lotinib is associated with a second mutation in the EGFR kinase domain. PloS Med 2(3):173. 12. Antonescu CR, et aL (200S)Acquired resistance to imatinib in gastrointestinal stromal tumor occurs through secondary gene mutation. Clin Cancer Res 11(1)4182-4190. 13. Mare T. Elde CA. DelnInger MW (2007) Ba-AN kinase domain mutations. chug re- sistance, and the road to a cure for chronic myeloid leukemia. Wood 110PP-2242-2299. 19. Misak S. et al. (2012) Emergence of KRAS mutations and acquired resistance to anti- EGFR therapy in colorectal cancer. Nature 986(7409).532-536. IS. Boa I, et al (2013) Evolutionary dynamics of cancer in response to targeted corn. dilation therapy. 'Life 2:e00747. 16. Coldman Al, Goitre NI (1983) A model for the resistance of tumor cells to cancer chemotherapeutic agents. Math Shad 65(2)291-307. 17. Iwasa Y, Nowak MA. Midler F (2006) Evolution of resistance during clonal expansion. Genetics 172(4):2557-2566. It Durrett R, Moseley $ (2010) Evolution of resistance and progression to disease during clonal expansion of cancer. Theor Popo( Viol 7701:42-4a metastatic lesions in a patient. This information is pertinent to the development of second line treatments that aim to inhibit known resistance mutations. Materials and Methods Model. We model the growth and evolution of a metastatic lesion as a con- tinuous time multitype branching process (34). The growth of a lesion is initiated by a single cell sensitive to the drug. Sensitive cells produce a re- sistant cell at each division with probability ta and each resistant cell pro- duced by sensitive cells starts a new resistant type. Analysis. In ow analysis, we use the approximation that resistant cells pro- duced bysensitive cells appear as a Poisson process on the number of sensitive cells (17). For more details and derivations of our results, please see SI Ten. Simulations. We perform Monte Carlo simulations of the multitype branching process using the Gillespie algorithm (35). Between 5,000 and 10,000 sur- viving runs are used for each parameter combination. ACKNOWLEDGMENTS. We thank Bert Vogelstein for critical reading of the manuscript and Rick Durrett for discussion during the conception of this work. We are grateful for the support from Foundational Questions in Evolutionary Biology Grant RFP-12-17 and the John Templeton Foundation. 19. Komarova N4 Wodar2 0 (2005) Drug resistance in cancer: Principles of emergence and prevention. Prot Nat) Acad Sci USA 102(27):9714-9719. 20. Komarova N (2006) Stochastic modeling of drug resistance in cancer. 1 Timor Riot 239(3):351-366. 21. Hoene H, Iwasa Y, Medico F (2007) The evOlubOn of two mutations during clonal expansion. Genetics 17700:2209-2221. 22. Bozic I, Allen B, Nowak MA (2012) Dynamics of targeted cancer therapy. Trench Mot Med 18(6):311-316. 23. tuna SE, DelbrOck M (1943) Mutations of bacteria from virus sensitivity to virus re- sistance. Genetics 24(6)491-511. 29. Kessler DA Levine N (2013) Large population solution of the stochastic Luriattelauck evolution model. hot Nat? Aced Sci USA 110(29):11682-11687. 25. Kessler DA. Levine H (2014) Scaling solution in the large population limit of the general asyrnmetdc stochastic Luria-Oelduck evolution process. arXtv:1404.2407. 24. Fool, Midler F (2019) Evolution of acquired resistance toanticancer therapy.) Theor Riot 355:10-20. 27. Bailey NT) (1964) The Elements of Stochastic PrOCeUeS With Appecations to the maven Sciences (Wiley. New Yor). 28. Leder K, et at. (2011) Fitness conferred by BCR.ABL kinase domain mutations dc- terrnines the risk of pre-existing resistance in chronic myeloid leukemia. PloS ONE 601)427682. 29. Bettegowda C. et al. (2014) Detection of circulating tumor DNA In early- and late- stage human malignancies. Sc! Trans: Med 6(224):224ra24. 30. O'Connell N (1993) Yule process approximation for the skeleton of a branching process. I ADA, Frobab 30(3):725-729. 31. Ivianrubla SC. Zanette OH (2002) At the boundary between biological and cuttural evolution: The origin of surname distributions. f Thor 85o1216(4)A61-477. 32. Maruvka YE, %nab NM, Kessler DA (2010) Universal features of surname distribution in a subsample of a growing population. I Theor Rid 262(2):245-256. 33. Denary A, Luebedt EG, Moolgavkar SH (2005) A generalized Luna-Delbrkk model. Math alC4C1197(2k140-152. 39. Athreya KB, Ney PE (1972) Breathing Processes (Springer-Verlag. Berlin). 35. Gillespie DT (1977) Exact stochastic gmulation of coupled dwmical reactions. I PhYl Chem 61R512340-2361. mg Bozic and Nowak PNAS Early Edition I 5 of 5 EFTA01201688 Supporting Information IA VA a Bozic and Nowak 10.1073/pnas.1412075111 Si Text The Model. We model the growth of a metastatic lesion as a branching process (1) that starts from a single cell sensitive to treatment. Sensitive cells divide with rate b and die with rate d. The net growth rate of sensitive cells is r=b— d. During division. one of the daughter cells receives a resistance mutation with probability U. Resistant cells have birth and death rates bR and dR. Resistant mu- tations can be neutral in the absence of treatment, which means they have the same birth and death rates as sensitive cells, and we initially focus on this case. Alternatively, if c = (bR —4)/(b —d)). 1, then resistance mutations are advantageous before treatment; if c < 1, they are deleterious. We assume that mutation rate a is small final lesion size M S large, Mu > 1, and we are mostly interested in the behavior of the early surviving clones. Furthermore, since mutation rate u is small, we assume that the size of the resistant population is much smaller than the size of the sensitive population, and approximate the size of the sensitive population with the size of the lesion. Rate of Production of Mutants. We use the result (2) that the col- lection of tumor sizes at which resistance mutations are produced can be viewed (approximated) as a homogeneous Poisson process on "I. hl] with intensity u/(1 —d/b). The reasoning follows from the fact that the average total number of resistance mutations produced when there are exactly x sensitive cells in the population is given by 00 Rx bar ! kW& —41b [SI] where jx(t) is the probability that there are exactly x sensitive cells after time t and the factor I —d lb comes from only looking at lineages in which the tumor population did not go extinct. In the small mutation rate a limit, one can neglect the production of mutants to calculateMt) as pertaining to a single type branching process on the sensitive cells. Even though it seems that Iwasa et al. (2) were not aware of it, fr(t) was derived by Bailey (1) f2(1) = (1 —0(1 —fi)tpx-i. [S21 with a= (de —d)/(be" —d) and /I= (be —b)/(be" —d). Plugging in the expression for L into Eq. 1 leads to the rate at which mutants are produced when there are x sensitive cells R.T =u/(1 —d/b). (S3] A more intuitive way to prove this result is as follows: when the population contains exactly x sensitive cells, the probability that they will produce a mutant before going to x —1 or x + 1 sensitive cells is bit 1(b + d). The (average) number of occurrences of ex- actly x cells in the process is one plus the (average) number of returns of a biased random walk with p = bl(b +d). Multiplying the probability of producing a mutant cell while at state x with the number of occurrences of that state leads to fix =u/(1—d/b). Each mutant cell survives stochastic drift with probability I —d/b, so the tumor sizes at which mutations that survive stochastic drift are produced can be viewed as a Poisson prooess on . MI with intensity u. Because Al is large and a is small, we can replace the interval [l. M] by [0. Al without losing much accuracy (3). Size of the kth Resistant Clone. Let MR be the size of the sensitive population when the kth successful resistant mutant appears. Furthermore, let Yk denote the size of the kth resistant sub- population that survives stochastic drift when there are Al- sensitive cells, conditioned on A4 <M. By the time the sensitive population reached size Al, n can be approximated by MV/MR, where V is an exponentially distributed random variable with mean b/(b —d) (4). If FR(y)= Pr[YR sy] is the cumulative distri- bution of Yk, then, expanding on the reasoning in ref. 3, we have Fk(Ofts 1 — POW Age Ilk% SMI if =1 I Prob.Density[MR =zIMR SM] x Pr [V Yiriz]elz o 1 =1 i m (zu)k-le-=u ( 1 &E l (Mu)t rya\ Mb P ) exP(— as )dt (k — I)! 160 0 Evaluating the integral above leads to &Ws's I ( Mu t r(k)—r(k.mu +y — dylb) Afu .e.y —4/0 — nk.mtr) where r(k).(k— I)! and r(a,z)= Jr 9-le'dt is the incomplete Gamma function. The probability density function for Yk, is given by lyri _knit) — r(k. mu +lb) issi My) a rk(bMu)k(bMu + r(k)—r(k.mu) ' In particular, fork C Mu, we have F&&) # I chi Mil_dy/b) , 1S6) and fk(y) a k(I -d/b)(Mu)k(Mu +y-dylb)-14. [S7] Comparison of Formula S6 and the cumulative distribution func- tion for n obtained from 5.000 runs of the exact computer simulation of the branching process is shown in Fig. 18. Calculating the expected number of cells in the kth clone using Formula S7 (integrating from 0 to Al) and expanding it in the small u limit, we obtain bMu 2 E(Yk) Pli—I) de0(11 )! for k >2 and E(Yt) a b—Afr u (logL— 1) +0(u2). 1591 We can also obtain the median for the number of cells in the kth clone, Yki2 , from the cumulative distribution function (S6) v I/2 = MIN (2tik _ 0. k 15101 Ratio of Reactant Clone Sizes. To more precisely determine the relationship between the sizes of different subclones, we next Bock and Nowak wanv.pnas.orgkgikefdent/shortila1207SII I 1 of 2 EFTA01201689 IA a calculate the probability distribution of }.11 /Yr,: the ratio of sizes of first and kth clone. We again use the fact that the random variable describing the size of the kth clone, Yk, can be ap- proximated by VkM/Mk, where Ilk —Exp[bl(b -d)), and that Mk, size of the population on the arrival of the kth clone, is the sum of k exponential random variables with mean l/u. We have Pr[Ydrk Sri =Pr[VaM/Mr SxVkM/Mk] = Pz[VIIK SxAdi I A fk]. We note that = Vi/Vk is the ratio of two independent, identically distributed exponential random variables and thus its probability density function is fz(z)= II(z + 1)2. Similarly, W =MI/Mk — rii„tymi..1]+ - 1.2])— k - I] is a p-distributed ran- dom variable with probability density function fw(w)= (k- 1)(1 -w)k-2. It follows that Sx] = Plc Sx1V] a 1 = (k-I)(1-w) k-2 J —dzdw (z + 1)2 = — - k - k — 2F1[1,1.1 + c. -xi) I +x [S11] where 2F1 is the hypergeometric function. Notably, this distri- bution depends only on k and not on any parameters of the process. For example, the probability that the ratio of sizes of the first and the second clone is smaller than x is Pr[Yi /Y2 Lx] — I log(I +x) x (S12] In particular, the first successful subclone that appears is smaller than the second appearing subclone in 31% of cases. The prob- ability that the first appearing subclone is twice as large or larger than the second appearing clone is 55%, and the probability that it is 10 or more times larger is 24%. Nonneutral Resistance. In this section, we will obtain similar results for the probability distributions of clone sizes and their ratios in the case in which resistant cells have birth and death rates bR and 4, respectively. Tumor sizes at which resistance mutations ap- pear can still be viewed as a homogeneous Poisson process on [0,M] with intensity u/(I — d/b). However, the lineages of newly produced resistance mutations will escape extinction with prob- ability 1 —dR/bR, so the successful resistant subclones in this scenario will arrive with rate :IR = u(1 - dR IbR)I(1 - d lb). An- other change from the neutral case is that the size of the kth clone when the total population size is M can be approximated with r k wimkrU, where Mk is the population size when the kth successful resistance mutation appeared and U is an expo- nentially distributed random variable with mean bR/(bR -4). We recall that c = (bR —4)/(b—d). 1. flaky NTI (1964) The Elements of Stochestk Processes *ills AppEcations to the Natural Sciences (IMley, New York), 2. Masa Y, Nowak MA Micher F (2006) Evolution of resistance during clonal expansion. Gentles 172(4)2557-2566. With the above caveats, we can write the derivation of the cumulative distribution function for the size of the kth appearing resistant subclone, yk, similarly as before F*(y) = 1- Pli(Al IMO ` 11 ayImk =I—/ Prob.Density[Afk =21Alk M] x Pr [U dz JJJ i m uR(zuR)" r n or — (1 rffi I! re -m") (k — I)! fra a rRyt ) x exp Rb M• az. (S131 The difference is that in the nonneutral case the integral has to be evaluated numerically. For the ratio of clone sizes Yk/Yi in the nonneutral case, we have prErlinsx)=Pr[utwimocsxUdm/mkr] =Pr[Ul/Uk sy(Mi")1 = pk t1c_1 c J (241)2 dz dm, 0 0 = k— dw. I +wx c ‘1 t 0 1S141 We see that even when resistance is not neutral, the ratio of clone sizes depends only on the order of appearance and the relative fitness c and not on M, is, and the specific birth and death rates of cells. Ow formulas rely on the approximation Y —e3p[(bR -dR)t]Ll for the size of a resistant clone Y. where U is an exponentially distrib- uted random variable with mean bR/(bR -4) and timer is mea- sured from the appearance of the founder cell of the clone. This approximation assumes large r and loses accuracy for r close to 0 [e.g., it predicts that the average size of a clone at time 0 is bR/(bR —4) rather than I]. Thus, our formulas lose accuracy when successful resistant clones are produced shortly before reaching sae Al. Successful resistant clones appear as a Poisson process on the number of sensitive cells with rate 15R= u(1— dR/bR)/(I — d lb) when resistance is not neutral, and the first such clone will ap- pear when there are —1/uR sensitive cells. For our approxima- tion to hold, Al must be significantly larger than 1/uR, which is equivalent to Muc(b/bR) I. In other words, we assume that the relative fitness c is such that the expected number of suc- cessful resistant clones, MuR, is much larger than I. 3. Diaz IA Jr, et at. (2012) The molecular evolution of acquired resistance to targeted EGFR blockade in colorectal cancers. 'More 486(7404):537-540. 4. Durrett A Moseley S (2010) Evolution of resatance and progression to disease during clonal expansion of cancer. Theo, POput 6,01 77(1 )42-03 Bo& and Nowak inmeiv.pnas.Org/Cgikentent/shOrl/14120751I I 2 aft EFTA01201690