## Abstract

Acquired drug resistance is a major limitation for cancer therapy. Often, one genetic alteration suffices to confer resistance to an otherwise successful therapy. However, little is known about the dynamics of the emergence of resistant tumor cells. In this article, we consider an exponentially growing population starting from one cancer cell that is sensitive to therapy. Sensitive cancer cells can mutate into resistant ones, which have relative fitness α prior to therapy. In the special case of no cell death, our model converges to the one investigated by Luria and Delbrück. We calculate the probability of resistance and the mean number of resistant cells once the cancer has reached detection size *M*. The probability of resistance is an increasing function of the detection size *M* times the mutation rate *u*. If *Mu* ≪ 1, then the expected number of resistant cells in cancers with resistance is independent of the mutation rate *u* and increases with *M* in proportion to for advantageous mutants with relative fitness , to for neutral mutants (α = 1), but converges to an upper limit for deleterious mutants (). Further, the probability of resistance and the average number of resistant cells increase with the number of cell divisions in the history of the tumor. Hence a tumor subject to high rates of apoptosis will show a higher incidence of resistance than expected on its detection size only.

CANCER therapy often fails because acquired resistance enables cancer cells to grow despite the continuous administration of therapy. In some cases, a single genetic alteration is sufficient to cause resistance to cancer treatment. For example, treatment of chronic myeloid leukemia with the chemotherapeutic agent imatinib fails due to resistance acquired through a single point mutation in the kinase domain of BCR-ABL, the oncogene driving this disease (Gorre *et al.* 2001). This also holds true for mutations in the Bcl-2 (Schmitt *et al.* 2000), PI3K and RAS (Gupta *et al.* 2001), or NF-κB (Bentires-Alj *et al.* 2003) genes, respectively.

Consider the clonal expansion initiated by a single cancer cell that is sensitive to cancer therapy. The offspring of this cancer cell increase exponentially in abundance until the tumor is detected and cancer therapy is initiated. During clonal expansion, mutant cells can be generated that are resistant to therapy. In this article, we address the following questions: What is the probability that a tumor of a certain size contains resistant mutants? And what is the distribution of these mutants in the tumor?

These questions arise in many different contexts. Resistance mutations emerging during the exponential growth phase of a bacterial infection cause failure of antibiotics treatment (Lenski 1998; Levin *et al.* 2000; Levy 2001). Viral infections are treated with antiviral drugs only after an initial phase of exponential growth during which resistance mutations arise (Nowak *et al.* 1991; Bonhoeffer and Nowak 1997; Ribeiro *et al.* 1998; Nowak and May 2000; Richman 2001). Inherited childhood retinoblastoma ensues if retinal cells receive the second hit in the retinoblastoma tumor suppressor gene during the development of the retina (Knudson 1971). Chronic myeloid leukemia is diagnosed only after an extensive period of clonal growth and resistance is often detected at diagnosis (Michor *et al.* 2005). Of those chronic myeloid leukemia patients commencing imatinib (Gleevec) in early and late chronic phase, 6 and 14% develop resistance within the first year of therapy (Branford *et al.* 2003).

In this article, we develop a mathematical model to calculate the probability of resistance of an exponentially growing cell population that starts from one sensitive cell and reaches detection size *M*. We determine the probability distribution of resistant cells in the tumor of size *M*. Our model is related to the question asked by Luria and Delbrück (1943) concerning the mutations of bacteria that confer resistance to phages. The distribution of the number of mutants in an exponentially growing population is called Luria–Delbrück distribution. This distribution has been studied for more than half a century (Tlsty *et al.* 1989; Zheng 1999; Frank 2003). Several different models have been proposed that are all based on pure birth processes and do not include the possibility of cell death (Zheng 1999). However, in most situations of cancer growth, bacterial and viral infections, death of the cancer cells or infecting agents cannot be neglected. In this article, we consider the case in which each cancer cell, virus, or bacterium experiences random birth or death. The total population size is growing exponentially with random fluctuations. Furthermore, most theoretical studies of the Luria–Delbrück distribution consider the situation in which the time period of exponential expansion is known. Conversely, in many medical situations the time between the initial start of growth and the detection of a cancer or infection is unknown, whereas the population size of the cells or infectious agents can be estimated. Therefore, we are interested in the probability of resistance until a clonally expanding population has reached a certain size *M*. We introduce a birth-and-death process model and calculate the number of resistant mutants in the population once it has reached size *M*. The probability distribution depends on the growth rate of the exponentially expanding agents, the mutation rate per replication, the relative fitness of mutants, and the detection size *M.* We derive formulas for the probability of resistance, for the mean number of mutants, and for the distribution of the number of mutants. These results are confirmed by exact computer simulations of the stochastic process.

## THE MODEL

Consider an exponentially expanding cancer cell population starting from one cell that is sensitive to therapy. The cell population follows a continuous-time branching process. In a short time interval of length , a sensitive cancer cell divides with probability and dies with probability . A resistant cancer cell arises by mutation with probability *u* per cell division. Resistant cancer cells also follow a continuous-time branching process. In a time interval of length , a resistant cancer cell divides with probability and dies with probability . The relative fitness of resistant cancer cells as compared to sensitive cancer cells is given by . Resistant cancer cells can be advantageous (), neutral (), or deleterious (). Both sensitive and resistant cancer cells grow on average exponentially, which means and . Since the cancer cell population is initiated by a single cell, it might go extinct due to random events. Otherwise, the cancer cell population reaches detection size *M*. For the sake of mathematical simplicity, we assume that clinical detection occurs when the population size of sensitive cancer cells reaches *M*.

Denote by the expected number of resistant cancer cells produced by a tumor with exactly *x* sensitive cancer cells. We assume that the number of resistant cancer cells follows a Poisson distribution with mean . Consider the clones each initiated by one resistant cancer cell. Their population sizes follow continuous-time branching processes. Let *Z* be the number of resistant cancer cells at the time of diagnosis. The generating function of *Z* is given by . The generating function is useful for the calculation of the probability of resistance, of the average number of resistant cells, and also of higher moments of *Z.* The total number of resistant cells, *Z*, is the sum of the number of resistant cancer cells in each clone. The generating function for the number of resistant cancer cells at diagnosis in the clone whose initiating resistant cell was produced when there were *x* sensitive cancer cells is denoted by . The generating function of the total number of resistant cells, *Z*, is given by(1)See appendix a for derivation of Equation 1.

Let us consider the Markov chain of the number of sensitive cancer cells, *x*. Their number follows a birth-and-death process (Feller 1971). The cell division rate is denoted by *r*, the cell death rate is denoted by *d*, and resistance mutants arise with probability *u* per cell division. We consider the process starting from *x* = 1; it can end either with extinction, *x =* 0, or with exponential increase until reaching the detection size, *x* = *M*. The transition rate from *x* to *x* + 1 is given by , which can be approximated as because *u* ≪ 1. The transition rate from *x* to *x* − 1 is given by *dx*. Both 0 and *M* are absorbing states of the process. Mutation does not change the state.

Let be the probability that there are *x* sensitive cells at time *t.* Then we have(2a)where *x* = 2, 3, …, *M* − 2. Furthermore, we have(2b)(2c)The initial condition is and . The production rate of resistant cancer cells is given by the product of the cell division rate, *r*, the mutation rate per cell division, *u*, and the number of sensitive cancer cells, *x*. The expected time until the number of sensitive cancer cells is *x* is given by . Then the expected number of resistant cancer cells produced when there are *x* sensitive cancer cells is given by(3)The factor 1 − *d*/*r* arises because we consider the conditional mean number of resistant cancer cells and exclude evolutionary trajectories in which sensitive cancer cells go extinct.

In appendix b, we can derive . Therefore, equal numbers of resistant cancer cells are produced by any number of sensitive cancer cells, *x* = 1, 2, … , *M* − 1. Then the generating function, Equation 1, becomes(4)

After a single resistant cancer cell has been created at time *t*, the number of resistant cancer cells follows a continuous-time branching process with the generating function(5)The derivation is explained in appendix c. Equation 5 is the generating function for the number of resistant cancer cells at time *t* after the occurrence of the resistance mutation. Here we are interested in the number of resistant cancer cells once the number of sensitive cancer cells has reached detection size *M*. The time period during which the number of sensitive cancer cells increases from *x* to *M* is not a fixed constant, but rather a stochastic variable. Assuming deterministic population growth of the sensitive cancer cell population once a resistant cancer cell has been produced, we have . With and , the generating function becomes(6)

## THE PROBABILITY OF RESISTANCE

Let us first calculate the probability of resistance. This probability is given by 1 minus the probability of extinction of resistant cancer cells (Equation 6 with ). The probability of nonextinction of a clone of resistant cancer cells is given byThen the probability of resistance is given by(7)Here and . The sum over *x* was replaced by an integral because *M* is very large. Equation 7 can be used for any fitness of resistant cancer cells as compared to sensitive cancer cells. Resistant cancer cells can be advantageous (), neutral (), or deleterious (). If *Mu ≪* 1, then we can approximate Equation 7 as .

When resistant cancer cells are neutral as compared to sensitive cells, we have and . Then Equation 7 can be simplified as(8)

The traditional Luria–Delbrück distribution assumes that there is no cell death, . In that case, Equation 7 becomes(9)This result can also be derived from the fact that it takes *M* − 1 cell divisions until one cell has generated a population of *M* cells; if resistant mutants are produced at rate *u* per cell division and the produced mutants survive until detection time, then the probability of resistance is given by Equation 9.

The analytic formula for the probability of resistance is confirmed by the exact computer simulation of the stochastic process. We define two integer variables for the numbers of sensitive and resistant cancer cells, and , subject to the constraint . Each process is initiated with one sensitive cancer cell, *x* = 1. Let . The transition probabilities between states are given byOnce the total population size reaches the detection size, , the process is stopped.

We perform many independent runs to account for random fluctuations. Then we determine the fraction of runs in which resistant cells were produced. Figure 1 shows the very good agreement between the analytic formula for the probability of resistance, *P*, given by Equations 7 and 8 (lines) and the exact computer simulation of the stochastic process for the probability of resistance (circles).

Figure 1 also shows that the probability of resistance, *P*, depends on the product *Mu* rather than on *M* and *u* separately. If *Mu* ≪ 1, then *P* is proportional to *Mu* (Figure 1a). If the relative fitness of resistant cells, α, is constant, then *P* decreases with increasing growth rate of sensitive cells, *r* (Figure 1b). Therefore, a slowly growing cancer is more likely to develop resistance than a quickly growing one, because the former undergoes a larger number of cell divisions than the latter until reaching a certain size. If the growth rate of sensitive cancer cells, *r*, is constant, then *P* increases with the relative fitness of resistant cells, α (Figure 1c). This implies that advantageous mutations are more dangerous than deleterious mutations under the assumption of constant mutation rates.

## THE AVERAGE NUMBER OF RESISTANT CELLS

Let us now calculate the average number of resistant cancer cells from the generating function, Equation 4. We consider the conditional average, *i.e.*, the average number of resistant cells in cancers that have developed resistance before reaching size *M*. This average number is related to the time period during which cancer therapy is effective. The conditional average number of resistant cancer cells is given by . Note that . Then we have(10)Here and .

If *Mu* ≪ 1, we can further simplify as(11)

If resistant cells are neutral, because *a = r* and *b* = *d*, Equation 10 becomes(12)Figure 2 shows the agreement between the analytic formula for the conditional average number of resistant cancer cells given by Equations 10 and 12 (line) and the exact computer simulation of the stochastic process (circles). If resistant cancer cells are deleterious or neutral, then there is very good agreement between the analytic prediction and the exact computer simulation. If resistant cells are advantageous, however, the equations above greatly overestimate the number of resistant cells. Therefore, we consider a different approach for calculating the average number of advantageous resistant cancer cells in the next section.

If resistant cells are deleterious or neutral, then the dependence of the conditional average number of resistant cancer cells, , on the mutation rate, *u*, is weak. If *Mu* is small, is even independent of the mutation rate. Note that is the conditional average of resistant cells—the probability of resistance, however, increases with the mutation rate, as discussed above. The dependence of on the detection size, *M*, is less than linear. If resistant cancer cells are neutral, then increases in proportion to ln *M*. If resistant cancer cells are deleterious, then does not increase strongly with *M* and converges to an upper limit (Figure 2a). Furthermore, decreases with increasing growth rate of sensitive cancer cells, *r*, but the decrease is faster for deleterious than for neutral resistant cancer cells (Figure 2b). Finally, increases with the relative fitness of resistant cancer cells, α (Figure 2c).

## HEURISTIC FORMULA FOR THE NUMBER OF ADVANTAGEOUS RESISTANT CANCER CELLS

An unexpected result of our analysis is that the number of resistant cancer cells, , produced by *x* sensitive cancer cells is almost independent of *x*. If we classify new resistant cancer cells according to the number of sensitive cancer cells present at the time of the production of resistant cancer cells, then the latter are evenly distributed from 1 to *M* − 1. Once one resistant cancer cell has been produced, its offspring increase in number until reaching . This is the expected number of resistant cancer cells once the sensitive cancer cells, *x*, reach detection size *M*, and it is the average over exponentially growing lineages and extinct lineages. The mean number of resistant cancer cells conditional to nonextinction is given by , where . Equation 10 gives the average number of resistant cancer cells when *x* is distributed uniformly from 1 to *M* − 1. However, the conditional number of resistant cancer cells, , can increase beyond *M*. Hence Equation 10 or 11 cannot be used for advantageous resistant cells, .

We introduce a truncation as follows. The stochastic process stops when the sum of sensitive and resistant cancer cells reaches the detection size, = *M*. Let us approximate this process in the model by stopping when the number of sensitive cancer cells reaches half of the threshold size, *M*/2. We define a threshold value of *x* as , leading to . Then we replace in Equation 11 by for . Hence, if *Mu* ≪ 1, we have(13)Figure 2c shows the agreement of Equation 13 (lines) with the direct computer simulation (circles). Equation 13 shows the parameter dependence of the conditional average number of resistant cancer cells, , very clearly. First, is independent of the mutation rate, *u*, and of the growth rate of sensitive cancer cells, *r.* Note that only the conditional mean number of resistant cancer cells is independent of *u* and *r*; the probability of resistance, *P*, increases with *u* and decreases with *r*. Second, increases as a power function of *M*. Third, increases with increasing fitness of resistant cells, α.

When there is no cell death, , and when *Mu* ≪ 1, we have(14a)(14b)(14c)This result corresponds to the traditional Luria–Delbrück distribution (Zheng 1999).

## THE DISTRIBUTION OF THE NUMBER OF RESISTANT CANCER CELLS

On the basis of the approximation for the deterministic growth of resistant cancer cells, we can derive the distribution of the number of resistant cancer cells at detection time. We use the result of the uniform distribution of the number of sensitive cancer cells when a resistant cancer cell arises. The number of resistant cancer cells is given by , where and . Let . The probability that the number of resistant cancer cells is greater than is given by(15a)Note that this and all the following expressions give the probabilities conditional to the existence of some resistant cancer cells at diagnosis. These formulas apply if *y* is clearly smaller than *M*. Similarly, the probability that the number of resistant cells is between and is given by(15b)Figure 3 shows the agreement between the computer simulation and Equation 15. If ≤ 10% of the threshold *M*, then Equation 15 fits very well to the results of the computer simulation.

The same argument applies to the case of deleterious or neutral resistant cancer cells (). Here the conditional average number of resistant cancer cells reaches less than the maximum, *M*, but Equation 15 holds for . Figure 3 illustrates the results of the computer simulations and Equation 15.

Comparing these results with direct computer simulations, we find that the formulas above are very accurate for advantageous resistant cancer cells with *y* < 10% of *M.* However, they are not accurate for deleterious resistant cancer cells, especially for small *y*, because the logic of Equations 15 relies on the exponential growth of the resistant cancer cell lineages. Many deleterious resistant cancer cells eventually go extinct, but might not yet have disappeared at the time of detection. This effect causes an underestimation for small *y.* In appendix d, we derive a formula for the number of resistant cancer cells that follows from the generating function (Equations 1 and 6) by Taylor expansion with respect to . When *Mu* ≪ 1, we have(16)Equation 16 is especially useful for small *y* and is accurate for deleterious resistant cancer cells. If there is no cell death, we have explicit solutions (see appendix d). Equation 15 is accurate for advantageous and neutral resistant cells, and Equation 16 is accurate for deleterious resistant cells (Figure 3).

## DISCUSSION

In this article, we have analyzed the dynamics of resistance mutations emerging from an exponentially growing population of cancer cells. Resistant cancer cells can be advantageous, neutral, or deleterious as compared to sensitive cancer cells in the absence of therapy. We use a continuous-time branching process to calculate the probability of resistance once the tumor has reached a certain size. We also calculate the probability distribution of resistant cancer cells conditional to the emergence of resistance in the tumor until clinical diagnosis. The accuracy of our model is confirmed by exact computer simulations of the underlying stochastic process. Tables 1 and 2 show the formulas for the probability of resistance and for the average number of resistant cancer cells with and without cell death, respectively.

Consider the number of resistant cancer cells that arise when the number of sensitive cancer cells is *x*. One might expect that more resistant cancer cells are produced when *x* is large, because more sensitive cancer cells divide more often. On the other hand, if *x* is large, events changing *x* (*i.e.*, cell death and cell division) occur faster, and the number of sensitive cancer cells, *x*, stays at any given value for a shorter time. Our model shows that exactly the same number of resistant cancer cells is produced by different numbers of sensitive cancer cells, *x*. This uniform distribution allows us to derive a formula for the probability that the number of resistant cells at diagnosis exceeds a certain number.

The probability of resistance, *P*, is given by Equation 7. This equation is accurate for both deleterious and advantageous resistant cancer cells. If resistant cancer cells are neutral as compared to sensitive cancer cells, then the probability of resistance is given by an even simpler formula (Equation 8). The probability of resistance, *P*, increases with the detection size, *M*, the mutation rate, *u*, and the relative fitness of resistant cells, α. Conversely, *P* decreases with increasing growth rate of sensitive cancer cells, *r.* This implies that the probability of resistance increases with the number of cell divisions that occurred until clinical detection of the cancer. The number of cell divisions increases if the growth rate is small or the death rate is large. Therefore, cancer cell clones with genetic instability (Lengauer *et al.* 1998; Loeb 2001; Nowak *et al.* 2002) might have a larger probability of generating resistant cells not only because of increased mutation rates, but also due to increased apoptosis rates, as such cell clones are more likely to receive deleterious or lethal mutations.

We derive equations for the mean number of resistant cancer cells, , conditional that the cancer has produced at least one resistant cell until diagnosis. Interestingly, is independent of the mutation rate, *u*. However, increases with the relative fitness of resistant cells, α. The dependence of on the detection size *M* differs among neutral, deleterious, and advantageous resistant cancer cells: increases as if they are advantageous and increases as if they are neutral, but it increases and saturates to an upper limit for large *M* if resistant cells are deleterious. Note that this result holds for the number of resistant cancer cells averaged over tumors with one or more resistant cells and only when *Mu* ≪ 1.

One intriguing result is that the conditional average number of resistant cancer cells, , does not increase with *M* if *M* is large and if resistant cancer cells are deleterious (see Equations 11 and 14a). The conditional number of deleterious resistant cancer cells increases exponentially, but their growth rate is slower than the growth rate of sensitive cancer cells. However, the number of tumors with newly emerging resistance increases. Therefore, the number of resistant cancer cells averaged over tumors with preexisting resistance and tumors with newly emerging resistance is constant for very large *M*.

The unconditional mean number of resistant cancer cells at diagnosis is given by . This number always increases with the detection size, *M*. It is proportional to *M* for deleterious resistant cells, but is faster than linear for advantageous resistant cells. Therefore, a decrease in the detection size reduces the risk of resistance, especially if resistance is caused by advantageous mutations. In this case, however, resistance caused by deleterious mutations remains unchanged and the relative importance of deleterious mutations increases.

The total risk of resistance, , increases with the mutation rate, *u*, and decreases with the growth rate of sensitive cancer cells, *r.* Therefore, cancer cell clones with genetic instability have a much greater risk of resistance than cells with normal mutation rates.

Here we assumed constant rates of cell division and cell death. However, the tumor growth rate is likely to slow down once its size becomes very large due to interference between cancer cells. Our results suggest that the reduction of *r − d* would enhance the risk of resistance. A detailed analysis of the effect of cell–cell interactions is an important subject of future theoretical investigations.

## APPENDIX A: DERIVATION OF EQUATION 1

Let be the number of resistant cancer cells produced when there are *x* = 1, 2, … , *M* − 1 sensitive cancer cells. The clone produced by each resistant cancer cell contributes to the number of resistant cells at diagnosis. The number of resistant cells at diagnosis has the generating function . From the assumption that follow a Poisson distribution with mean and that they are independent of each other, their probability distribution is given by . The generating function of the total number of resistant cancer cells is given by(A1)This is Equation 1.

## APPENDIX B: SOLUTION OF EQUATION 3

We consider the integral of both sides of Equation 2 from *t* = 0–. We note and With , we have(B1a)(B1b)(B1c)From Equations 2a–2c, Equation B1b becomes , and then Equation B1a becomes . Hence . From Equation B1c, we haveWe further have(B2)except for those *x* that are very close to *M.* Using this result, we have(B3)

## APPENDIX C: DERIVATION OF THE BRANCHING PROCESS GENERATING FUNCTION

Let . Considering the events happening in a short time interval of length , we haveWe have used , because two lineages starting from different cells at a given time behave independently (Iwasa *et al*. 2004).

In the limit of , this leads to(C1)Using the initial condition , we have Equation 5 in the text. It satisfies (because the sum of all probabilities is 1), , which is the mean number of mutants, and , which is the extinction probability.

## APPENDIX D: THE NUMBER OF RESISTANT CELLS

We can derive a formula for the number of resistant cancer cells, *y*, directly from the generating function (Equation 1). Considering the case *Mu* ≪ 1, we have(D1)From Equation 6, we have(D2)Here and . Hence we have(D3)The coefficient of is(D4)This is the unconditional probability that the number of resistant cancer cells is *k.* The mean number conditional to the existence of resistance is given byThis is Equation 16 in the text.

If there is no cell death, we have *d* = *b =* 0, and hence(D5)

If in addition the resistant cancer cells are neutral (), we have(D6)

## Acknowledgments

This work was done during Y.I.'s visit to the Program for Evolutionary Dynamics at Harvard University in 2004 and 2005. A Japan Society for the Promotion of Science grant-in-aid for scientific research to Y.I. is gratefully acknowledged. The Program for Evolutionary Dynamics is supported by Jeffery Epstein.

## Footnotes

Communicating editor: P. Oefner

- Received August 18, 2005.
- Accepted October 3, 2005.

- Copyright © 2006 by the Genetics Society of America