# Timing and heterogeneity of mutations associated with drug resistance in metastatic cancers

See allHide authors and affiliations

Edited by Herbert Levine, Rice University, Houston, TX, and approved October 8, 2014 (received for review June 28, 2014)

## 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 radiographically detectable metastatic lesions harbor multiple resistant 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.

## Abstract

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 consequence 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 subpopulation. Using the novel approach of ordering the resistant subclones 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 clone divided by the median of the first is

Acquired 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 resistance 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 colorectal cancer patients is a *fait 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). Successful 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 mutations in patients with metastatic cancers.

First mathematical investigations of the evolution of resistance to cancer therapy were concerned with calculating the probability 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–Delbrück model for accumulation of resistant cells in exponentially growing bacterial populations (23). Probability distribution for the number of resistant cells in a population of a certain size in the fully stochastic formulation of the Luria–Delbrück 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 *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

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. 1*A*). A reasonable assumption for the number of point mutations that can provide resistance to a targeted 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 *M* cells. Typical radiographically detectable lesions are *u*, leading to resistance is the product of the point mutation rate μ, which is on the order of *M* and small *u* limit and mostly focus on the case when

Tumor sizes at which successful resistant mutants are produced can be viewed as a Poisson process on *u* (*SI Text*) (10, 17). The number of successful mutant lineages is thus Poisson distributed with mean *k*th mutant appeared, which survived stochastic drift (Fig. 1*A*), then *k*th clone appeared when the total population size was *k* times the size of the *k*th clone. The probability that exactly *k* clones are present in the population of size *M* is

Counting new successful resistant clones in the order of appearance, we calculate the probability distribution for the number of cells in the *k*th resistant clone. In particular, if *k*th clone simplifies to**1** and exact computer simulations of the stochastic process is shown in Fig. 1*B*.

The mean number of cells in the *k*th resistant clone is *k*th subclone is given by*k* and *j* is

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 (

In *SI Text*, we calculate the probability distribution for the ratio of resistant clone sizes

Fig. 2 shows different realizations of the stochastic process of evolution of resistance in metastatic lesions containing

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 mutation 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 fragments 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 mutant 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 excellent agreement between the data and our model predictions. For example, the median ratio of the sizes of the first two resistant 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 that the ratio of resistant clone sizes is independent of parameters.

Our mathematical results describe the relative sizes of resistant clones ordered by age, whereas the experimental data in Table 1 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, *M*, and small mutation rate, *u*, with

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* is

## Discussion

In this paper we describe the heterogeneity of mutations providing resistance to cancer therapy that can be found in any one 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 process and that the number of lesions is much smaller than *k*th 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 depend on the parameters of the process, their ratios are generally parameter free. The universality of the clone ratio statistics follows 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 *k*th largest clone is *k*th and *j*th largest clones is

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–Delbrück model in which sensitive cells grow deterministically and calculated the number of individual resistant clones and the probability distribution for the number of cells in a single resistant clone after time *t*. In another study (28), mathematical modeling along with in vitro growth rates of cells harboring 12 point mutations providing resistance to BCR-ABL (fusion of breakpoint cluster region gene and Abelson murine leukemia viral oncogene homolog 1) 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

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 mutations is further amplified when taking into account multiple 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 continuous time multitype branching process (34). The growth of a lesion is initiated by a single cell sensitive to the drug. Sensitive cells produce a resistant cell at each division with probability *u* and each resistant cell produced by sensitive cells starts a new resistant type.

### Analysis.

In our analysis, we use the approximation that resistant cells produced by sensitive cells appear as a Poisson process on the number of sensitive cells (17). For more details and derivations of our results, please see *SI Text*.

### Simulations.

We perform Monte Carlo simulations of the multitype branching process using the Gillespie algorithm (35). Between 5,000 and 10,000 surviving 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.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: ibozic{at}math.harvard.edu or martin_nowak{at}harvard.edu.

Author contributions: I.B. designed research; I.B. and M.A.N. performed research; I.B. and M.A.N. analyzed data; and I.B. and M.A.N. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1412075111/-/DCSupplemental.

## References

- ↵.
- Vogelstein B, et al.

- ↵
- ↵
- ↵
- ↵.
- Komarova NL,
- Wodarz D

- ↵
- ↵
- ↵.
- Katayama R, et al.

- ↵
- ↵
- ↵
- ↵.
- Antonescu CR, et al.

- ↵.
- O’Hare T,
- Eide CA,
- Deininger MW

- ↵
- ↵.
- Bozic I, et al.

- ↵
- ↵.
- Iwasa Y,
- Nowak MA,
- Michor F

- ↵
- ↵.
- Komarova NL,
- Wodarz D

- ↵
- ↵.
- Haeno H,
- Iwasa Y,
- Michor F

- ↵
- ↵.
- Luria SE,
- Delbrück M

- ↵.
- Kessler DA,
- Levine H

- ↵.
- Kessler DA,
- Levine H

- ↵
- ↵.
- Bailey NTJ

- ↵
- ↵.
- Bettegowda C, et al.

- ↵
- ↵
- ↵
- ↵
- ↵.
- Athreya KB,
- Ney PE

- ↵

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Biophysics and Computational Biology

- Physical Sciences
- Applied Mathematics