Evolutionary dynamics of cancer in response to targeted combination therapy

In solid tumors, targeted treatments can lead to dramatic regressions, but responses are often short-lived because resistant cancer cells arise. The major strategy proposed for overcoming resistance is combination therapy. We present a mathematical model describing the evolutionary dynamics of lesions in response to treatment. We first studied 20 melanoma patients receiving vemurafenib. We then applied our model to an independent set of pancreatic, colorectal, and melanoma cancer patients with metastatic disease. We find that dual therapy results in long-term disease control for most patients, if there are no single mutations that cause cross-resistance to both drugs; in patients with large disease burden, triple therapy is needed. We also find that simultaneous therapy with two drugs is much more effective than sequential therapy. Our results provide realistic expectations for the efficacy of new drug combinations and inform the design of trials for new cancer therapeutics. DOI: http://dx.doi.org/10.7554/eLife.00747.001


Branching process model
As described in the main text, we model the evolution of resistance as a stochastic branching process, in which each resistance profile is identified as a distinct type. Prior to treatment, all cell types divide at rate b and die at rate d. Thus the tumor expands at rate r = b − d prior to treatment. Treatment is initiated when the tumor has reached a detection size of M cells. During treatment, we suppose that cell types sensitive to at least one drug (i.e. those with resistance profiles other than 1 . . . 1) divide and die at rates b and d , respectively, with r = b − d < 0. Fully resistant cells (those with profile 1 . . . 1) continue to divide and die at rates b and d, respectively. (More generally, one could suppose that each resistance profile has distinct birth and death rates during treatment, but we do not consider such generality here.) With each reproduction, one of the daughter cells acquires each potential point mutation with probability equal to the point mutation rate u. (Our model assumes that only one of the daughter cells can acquire mutations, which amounts to rescaling the mutation rate by a factor of two.) Thus, for each combination of drugs i 1 , . . . , i m , the probability that one of the n i 1 ...i d resistance mutations occurs in a daughter cell is (1 − u) n i 1 ...im ≈ 1 − n i 1 ...im u. (This approximation assumes that un i 1 ...im 1-that is, resistance mutations are rare-so that the possibility of multiple resistance mutations in a single reproduction event can be disregarded.)

Paths to full resistance
We are most interested in the emergence of simultaneous resistance to all D drugs, leading to treatment failure. Such multi-drug resistance may occur via multiple paths, where a path is a sequence of resistance mutations leading from 0 . . . 0 (no resistance) to 1 . . . 1 (resistance to all drugs). An example for three drugs is the path 000 → 010 → 111 (resistance to drug 2 first, then 1 and 3 simultaneously). Each path involves some number m of mutations, 1 ≤ m ≤ D. Along a path there are m + 1 resistance profiles, indexed j = 0, . . . , m. The initial resistance profile 0 . . . 0 is indexed j = 0, while the final, 1 . . . 1, is indexed j = m. The number of potential point mutations that would lead from profile j − 1 to profile j, along a particular path, is denoted ν j . The number of potential point mutations that would lead from profile j −1 to a profile not on the path in question is denoted η j . In short, ν j and η j are the numbers of "on-path" and "off-path" mutations, respectively, from profile j − 1. So, for example, for the path 000 → 010 → 111 we have ν 1 = n 2 η 1 = n 1 + n 3 + n 12 + n 13 + n 23 + n 123 ν 2 = n 13 + n 123 η 2 = n 1 + n 3 + n 12 + n 23 .
Thus, with each division of a cell of profile j − 1, we have the following probabilities of events (assuming (ν j + η j )u 1, i.e. resistance mutations are rare): 1 − (ν j + η j )u two daughters of profile j − 1 ν j u one daughter of profile j − 1 and one of profile j η j u one daughter of profile j − 1 and one off-path daughter.
1.4 The rare-mutation, large-tumor-size limit In human cancers, the point mutation rate is very small (u ∼ 10 −9 ) and the number of cells in a detectable tumor is very large (M ∼ 10 9 ). Therefore, we concentrate on results that are asymptotically exact under the following limits: Above, k will be 1 or 2, depending on the result being presented. We will represent such a limit by the arrow .

Number of resistant mutants at detection
We are first interested in following question: how many cells in the tumor are resistant to all D drugs at the time of detection?
We assume for the moment that ν k + η k = ν + η for all pairs k, = 1, . . . , m with k = . Under this assumption, the solution to the system (2) is given by We note that the above expressions for x j (t) are expectations over all possible trajectories, including those in which the tumor becomes extinct. Since we are only interested in tumors that grow to detectable size, we divide these expressions by the tumor survival probability, which is r/b. This yields an exact expression for the expected number of cells of each type at time t, conditioned on the survival of the tumor.
However, we wish to know the expected number of resistant cells not at a fixed time from when the tumor was initiated, but at the moment the total number of cells reaches M . The time T for the tumor to reach M cells is a random quantity; however, we can approximate T by using the deterministic growth law x(t) = b/r e rt for the total population of cells. (In other words, we set the total cell population equal to its expected value conditioned on non-extinction.) We then solve x(T ) = M , yielding the approximation Simulation results (not shown) suggest that this approximation is exact in the large tumor size, small mutation rate limit (1). We approximate the expected number x det m of fully resistant cells at detection (arising via this particular path) by substituting approximation (4) for T into the expression for x m (t) given in (3). This yields (5) To simplify the above expression, we define µ = (b/r) u log(M r/b). We note that µ → 0 in the rare mutation, large tumor size limit u → 0, M → ∞, M u = constant. Substituting µ into (5) and replacing e −(ν k +η k )µ by its Taylor expansion, we obtain For any collection of m distinct nonzero real numbers α 1 , . . . , α m -in our case, we are interested in α j = ν j + η j -the following combinatorial identity holds: We save the proof of this identity for Section 7. Using this identity to simplify (6), we obtain Interestingly, this formula does not involve the numbers η j of off-path mutations. As a consequence, we see that this result does not depend on the assumption that ν k + η k = ν + η for k = , and holds regardless of whether this condition is satisfied.

One drug
For single-drug therapy (D = 1), there is only one path to resistance, and the expected number of resistant cells at detection is given by (5), which simplifies to (We recall from above that µ = (b/r) u log(M r/b).)

Two drugs
For two-drug therapy, there are three paths to consider: Path 1: 00 → 10 → 11 In this case we have ν 1 = n 1 and ν 2 = n 2 + n 12 .

Arbitrary number of drugs with no cross-resistance
Suppose a therapy consists of D ≥ 1 drugs, and there are no mutations that simultaneously confer resistance to multiple drugs. In this case, using (8), we obtain a remarkably simple formula for the expected number of resistant cells at detection: For example, in the case of triple therapy with no cross-resistance, we expect approximately M n 1 n 2 n 3 µ 3 resistant cells at detection (accurate to order µ 4 ). Notice that the factorial in (8) is cancelled by the number, D!, of possible paths to full resistance.

Generating functions for branching processes
Generating functions are a powerful tool for analyzing stochastic processes. In a generating function, the probabilities of different events are recorded as coefficients in a power series, allowing these probabilities to be easily manipulated. Here we introduce the generating functions for the one-and two-type branching processes, which we will later use to derive probabilities of resistance and of treatment success.

One-type branching process
We first consider the one-type branching process. We introduce the random variable Y (t) to represent the number of cells at time t, given that there was one such cell at time t = 0. The generating function for this process is then defined as In words, the generating function is a time-dependent power series in which the coefficient of z y equals the probability that there were y cells at time t. There is a well-known closed-form expression for this generating function (Athreya and Ney, 2004): In particular, the probability that a lineage of type 1 cells survives for time t is given by 3.2 Two-type branching process We now turn to a branching process involving two types, labeled 1 and 2, with one-way mutation of rate uν 2 from type 1 to type 2. This process can be described by the following rates: We introduce the random variables Y 1 (t) and Y 2 (t) to represent the numbers of type 1 and type 2 cells at time t, given that the process was initiated with a single type 1 cell at time t = 0. The generating function for this process is defined as In words, the generating function is a time-dependent power series in the variables z 1 and z 2 , with the coefficient of z y 1 1 z y 2 2 equal to the probability that there are y 1 cells of type 1 and y 2 cells of type 2 at time t.
A closed-form expression for this branching process was discovered by Antal and Krapivsky (2011). To state this solution, we first define the following constants: Next, we let F denote the hypergeometric function 2 F 1 , and we define the following functions of a real number x: Third, we define the following quantities, which depend on the arguments z 1 , z 2 , and t, of the generating function φ m=2 (z 1 , z 2 ; t): Finally, we state Antal and Krapivsky's (2011) formula for the generating function of the two-type branching process in terms of the above quantities and functions: . (14) 4 Probability of resistance at time of detection We now turn to the question of whether at least one resistant cell exists at the time the tumor reaches detectable size. Again we consider each path to resistance separately.

One-step paths
A number of works (Coldman and Goldie, 1983;Dewanji et al., 2005;Komarova and Wodarz, 2005;Iwasa et al., 2006) have investigated the probability that resistance exists at the start of treatment, in the case that this resistance can be achieved through a single mutation. Here we follow the approach of Dewanji et al. (2005), in which we suppose that type 0 (sensitive) cells grow deterministically, and that type 1 (resistant) cells arise as a Poisson process, with rate depending on the current number of type 0 cells. Specifically, we approximate the growth of type 0 cells (conditioned on non-extinction) by Noting that type 0 cells divide at rate b, and each division produces a type 1 mutant with probability uν 1 , we suppose that type 1 (resistant) cells arise as a Poisson process with rate buν 1 x 0 (t) at time t. Each line of type 1 cells thus created is described by a one-type branching process initiated at the time of mutation. We let the random variable X 1 (t) denote the number of type 1 cells at time t. The probability that there are no resistant cells at the time of detection (T = 1/r log(M r/b)) is equal to the probability that none of the type 1 mutations that arise in the interval [0, T ] survive to time T . The probability that a single type 1 mutation arising at time 0 ≤ s ≤ T survives until time T is given by 1 − φ m=1 (0; T − s). Recalling that these mutations arise at rate buν 1 x 0 (s) for 0 ≤ s ≤ T , we can write this probability as Substituting from (12) and simplifying, we arrive at This result was also obtained by Iwasa et al. (2006).

Two-step paths
We can apply similar methods to two-step paths (m = 2). The probability that, time t after a type 1 cell arises, this cell's lineage includes at least one surviving type 2 cell, can be written as 1 − φ m=2 (1, 0, t). Thus a type 1 mutation that arises at time s, 0 ≤ s ≤ T , will give rise to at least one living type 2 cell at time T , with probability 1 − φ m=2 (1, 0, T − s). Following the arguments in the m = 1 case, the probability of no type 2 cells at time T can be written as

Low-mutation expansion for two-step paths
For low mutation rates, we can expand (17) as follows. First, we work through the construction of φ m=2 (z 1 , z 2 ; t), substituting z 1 = 1 and z 2 = 0, and expanding all quantities and functions to low orders in u. This yields: Combining the above expansions yields the following expansion for φ m=2 (1, 0; t)− 1: Substituting (18) into formula (17) for the probability of no resistance yields Above, we have also used the substitution M = b/r e rT . We simplify the integral in (19) as follows: The remaining integral above is positive and bounded above by where Li 2 is the dilogarithm function. In particular, this bound is constant with respect to T , and therefore becomes negligible in comparison to log(b/r)T as T becomes large. This implies the following asymptotic formula for the integral in (19): Finally, substituting the above result into (19), and additionally substituting T = log(M r/b)/r, we obtain Our formulas (17) and (20) improve on results obtained by Haeno et al. (2007), as we discuss in Section 6.2.

Overall probability of resistance
One drug For a single drug, the probability that resistance exists at the time of detection can be obtained directly from (16) with ν 1 = n 1 , yielding Two drugs For two drugs, we must consider the one-step path 00 → 11 and the two two-step paths 00 → 10 → 11 and 00 → 01 → 11. The probability that resistance exists at the time of detection equals one minus the probability that no cells of profile 11 are generated via any of these paths. We write this probability as p res = 1 − p 1 p 2 , where is obtained from (16) with ν 1 = n 12 , and is obtained from (20) with ν 1 = n 1 , ν 2 = n 2 +n 12 for the path 00 → 10 → 11, and ν 1 = n 2 , ν 2 = n 1 + n 12 for the path 00 → 01 → 11. Rewriting in terms of s = 1 − d/b yields the formulas for p 1 and p 2 presented in the main text.

Probability of tumor eradication
We now turn to the ultimate success or failure of multi-drug therapy. We define the therapy as successful if the tumor ultimately becomes extinct; otherwise, the tumor grows exponentially and there is a relapse. This differs from the question of analyzed in Section 4-the probability that resistance is present at the start of treatment-for two reasons. First, resistance that exists at the start of treatment may disappear due to stochastic drift. Second, new resistance mutations may appear during therapy. We recall that, for all cell types sensitive to at least one drug (that is, all cells with resistance profiles other than 11 . . . 1), the birth and death rates during treatment are denoted b and d , respectively, with r = b − d < 0. Cells with resistance to all drugs are unaffected.
We separate the question of ultimate treatment outcome into resistance arising during tumor expansion and resistance arising during treatment. We let p ↑ and p ↓ denote the probability that no resistance mutations leading to relapse arise during expansion and treatment, respectively. The overall probability of eradication can then be written

Resistance arising during expansion
We further separate into paths leading toward resistance.

One-step paths
As in Section 4.1, we suppose that type 1 (resistant) cells arise as a Poisson process with rate buν 1 x 0 (t) at time t, with x 0 (t) = b/r e rt . We also recall that each type 1 (resistant) mutation that arises has probability r/b of ultimately escaping stochastic drift (leading to unchecked tumor growth and patient relapse). Using similar reasoning to Section 4.1, the probability that no resistance mutations leading to treatment failure arise by this path, during tumor expansion, can be written as This result was previously obtained by Komarova (2006, Appendix A).

Two-step paths
During the treatment phase, the lineage of each type 2 (fully resistant) cell will ultimately disappear with probability χ 2 = d/b. For a type 1 cell that is present at the start of treatment, the probability χ 1 that its lineage will ultimately disappear can be obtained from the results of Antal and Krapivsky (2011): Above, b and d are the birth and death rates, respectively, of type 1 cells during treatment, and r = b − d . Suppose a type 1 mutation gives rise to y 1 type 1 cells and y 2 type 2 cells at the time of detection. Then the probability that all lineages of these cells disappear during treatment (and thus none of them cause eventual relapse) is χ y 1 1 χ y 2 2 . Overall, for a type 1 mutation that arises at time s, 0 ≤ s ≤ T , the probability that the lineage of this cell eventually disappears can be written as Following the arguments of previous sections, the probability that no type 1 mutation leading to treatment failure arises during expansion can be written as This formula involves the expression φ m=2 (χ 1 , χ 2 ; T − s). This expression cannot be evaluated directly using formula (14) for φ m=2 , obtained by Antal and Krapivsky (2011), because this formula is undefined (has a removable singularity) at z 2 = χ 2 = d/b. We therefore derive an expression for φ m=2 (z 1 , d/b; t) from first principles. To begin, we note that dynamics of the two-type branching process (13) satisfy the backward Kolmogorov equations, which can be written as: Above, φ ≡ φ m=2 (z 1 , z 2 ; t) is the generating function for the two-type process starting with a cell of type 1, whileφ(z 1 , z 2 ; t) is the analogous generating function starting with a cell of type 2. From the second equation in (24), we can see thatφ for all z 1 and all t > 0. Thus, for the fixed value z 2 = χ 2 = d/b, the backward Kolmogorov equations (24) reduce to This is a Ricatti equation with constant coefficients, which can be solved using standard techniques. We find the solution Above,

Low-mutation expansion for two-step paths
We now derive an expression for p ↑ m=2 , the probability that no resistance arises via the path in question during tumor expansion, that is asymptotically exact in the limit of small mutation rate u and large detection size M . We know χ 2 = d/b, and from (22) we have The formula (26) for θ ± admits the following low-mutation expansions: We substitute these expansions into (25) and note that, for values of t that contribute significantly to the integral in (23) (with t identified as T − s), the terms containing e rt eclipse those that are constant in t. As u → 0 and e rt → ∞, we have To obtain p ↑ m=2 we substitute (27)  (28) The integral on the right-hand side above simplifies as follows: Finally, by expressing everything in terms of the tumor size at detection, M = b/r e rT , we arrive at 5.2 Resistance arising during treatment

One-step paths
For one-step paths, the probability that no resistance arises during treatment was obtained by Michor et al. (2006): and p ↓ 2 = exp −M u 2 2n 1 n 2 + n 12 (n 1 + n 2 ) with a = b/r − b /r as above. Substituting r/b = s and r /b = s yields the formulas presented in the main text. We can also combine the above expressions to obtain Formula (32) improves on prior results of Komarova (2006), as we discuss in Section 6.3.

Comparison to previous results
Aspects of the dynamics of combination cancer therapy and resistance have been investigated in a number of previous works-notably Komarova and Wodarz (2005), Komarova (2006), and Haeno et al. (2007). Our results improve on the results previously obtained in these works and provide a closer match to simulations. This improvement is due in part to our use of recent advances in the theory of branching processes (Antal and Krapivsky, 2011), which allow us to obtain results that are asymptotically exact in the rare-mutation, large-tumor-size limit.

Expected number of resistant cells at detection
The question of the expected number of resistant cells in a tumor of detectable size has also been investigated by Iwasa et al. (2006), in the case of one drug, and Haeno et al. (2007), in the case of two drugs (with no mutations conferring resistance to both simultaneously). For one drug, in the case that resistant cells have the same division and death rates as sensitive cells in the absence of treatment, formula (10) of Iwasa et al. (2006) gives the following expression for the expected number of resistant cells:    Table S3: Probability p res of resistance at detection for dual therapy with no cross-resistance, as calculated using simulation, using our formulas (17) and (20), and using the formulas (35) and (36)  (see Appendix A). In contrast, our formulas (17) and (20) for P X 2 (T ) = 0 utilize the exact solution obtained by Antal and Krapivsky (2011). The two formulas derived by Haeno et al. (2007) both describe the probability that, at the time of detection, at least one cell contains two mutations arising in a specified order. (In other words, their formulas apply to the case of a particular two-step path.) Their main formula can be expressed in our notation as They also present an alternative formula: .
(36) In Table S3 we compare values of the overall probability of no resistance at detection p res (along all paths) as obtained from simulation, as derived from our formulas (17) and (20), and as derived from formulas (35) and (36) of Haeno et al. (2007).

Probability of tumor eradication
Turning to the question of whether therapy will successfully eradicate a tumor, our results in the case of one-step paths to resistance agree with those of Komarova (2006) and Michor et al. (2006), as noted where these results are presented.
For two-step paths, Komarova (2006) obtained closed-form approximations for the probability of tumor eradication, in the special case that an equal number of mutations are required for each step (ν 1 = ν 2 ), and that treatment affects only the death rate of tumor cells, not the division rate (b = b). Komarova's (2006) expression for the probability p ↑ -that no mutations leading to eventual relapse arise during tumor expansion-can be expressed in our notation as Above, ν = ν 1 = ν 2 is the number of resistance mutations at each step, and M 0 is an extra parameter representing tumor's initial size. Our approach of considering the entire history of the tumor, conditioned on its survival, amounts to setting M 0 = b/r. (We also note that Komarova uses P ↑ and P ↓ to denote the probabilities that mutations leading to relapse do arise during the two respective phases; thus P ↑ in Komarova's notation corresponds to 1 − p ↑ in ours, and similarly for P ↓ and p ↓ .) The approximation (37) for p ↑ is based on first determining the probability that at least one type 2 mutation arises, then multiplying that quantity by the survival probability of each such mutation. This method is accurate if few type 2 mutants are likely to be generated (M u 2 ν 2 1), but loses accuracy if many type 2 mutations may arise, because it does not take into account the individual fate of each mutation. Our expression (29), based on the exact formula for φ m=2 (z 1 , z 2 ; t) obtained by Antal and Krapivsky (2011), does not have this limitation.
For the probability p ↓ that no mutations leading to relapse arising during treatment, Komarova (2006) obtained This expression coincides with our formula (31) in the special case ν 1 = ν 2 = ν and b = b. Our result (31) can therefore be seen as a generalization of Komarova's result (38) to the case that mutation numbers may be unequal and treatment may affect tumor cell division rates. Table S4 compares our formula (32) for the overall probability p erad of tumor eradication for two-drug therapy to the results of Komarova (2006).
Finally, we note that Komarova (2006) also provides computational recipes to obtain probabilities for treatment success to arbitrary numerical precision. This methodology has been applied to a number of specific questions regarding treatment of chronic myeloid leukemia (Komarova and Wodarz, 2005;Komarova et al., 2009;Katouli and Komarova, 2010;Komarova, 2011). 7 Proof of Identity (7) In this section we prove the mathematical identity used in Section 2.1. We state this identity in an equivalent form: . Now consider a closed curve γ in the complex plane, whose interior contains the points α 1 , . . . , α m , as well as zero. Applying the residue theorem and Cauchy's integral formula yields On the other hand, since F has no poles outside of γ other than possibly at ∞ (in the case s ≥ m) we can also express this integral in terms of the residue of F at ∞: In particular, for s = m we have Res [F (z), ∞] = (−1) m . Combining with (40) verifies the desired result (39), and moreover, provides a method to calculate similar expressions with s > m.