Our model is based on a multitype branching process (see ‘Materials and methods’). Similar mathematical modeling has successfully predicted the dynamics of acquired resistance, including the timing of treatment failure, in colorectal cancer patients treated with the EGFR inhibitor panitumumab (Diaz et al., 2012), and has led to specific recommendations for combination therapies to treat CML (Komarova et al., 2009; Katouli and Komarova, 2010). Our current work builds on these previous studies by using recent advances in the mathematical theory of branching processes (Antal and Krapivsky, 2011), which enable us to obtain results that are exact in the biologically relevant limit of many tumor cells and small mutation rate.

To obtain key parameters for our model, we have studied the dynamics of 68 index lesions in 20 melanoma patients receiving the BRAF inhibitor vemurafenib. The data from six patients that represented distinct patterns of responses are shown in Figure 1. Patients P1 and P2 achieved complete responses, and their lesions became undetectable. Patient P3 had stable disease, with tumors remaining approximately the same size throughout treatment. Patients P4 to P6 all had partial remissions, with some lesions shrinking and others unchanging or regrowing during treatment. As expected, the smallest lesions were the ones most likely to become undetectable when the agent was effective.

Figure 1 Download asset Open asset Variability in treatment response to monotherapy among six patients. Patients were treated with the BRAF inhibitor vemurafenib. Patients P1 and P2 achieved a complete response. Patient P3 had stable disease. Patients P4, P5, and P6 had partial responses. The minimal detection size (indicated by discontinuous red line) was assumed to be ≈63 × 106 cells. https://doi.org/10.7554/eLife.00747.003

For 21 lesions in our vemurafenib dataset, two pretreatment measurements were available. Using these data, we calculate the average net growth rate of these lesions to be 0.01 per day, which is consistent with previous reports (Friberg and Mattson, 1997; Eskelin et al., 2000). The estimated average time between cell divisions in the absence of cell death in melanoma cells is 7 days (Rew and Wilson, 2000), implying a birth (cell division) rate of b = 0.14 per day. We set this as the typical birth rate, and additionally explore birth rates that correspond to a wide range of 1–14 days between cell divisions (Supplementary file 3). To achieve the observed net growth rate, we set the cell death rate to d = b − 0.01 (typical d = 0.13). We assume that these birth and death rates are valid for all cell types prior to treatment. For simplicity, we assume that these birth and death rates remain constant for all cell types prior to treatment, and neglect variations in the growth rate due to spatial and metabolic constraints in solid tumors (Bozic et al., 2012).

A given cancer therapy will reduce the birth rate and/or increase the death rate of tumor cells. A cell type is defined as sensitive if the treatment in question would cause its death rate to exceed its birth rate; otherwise, it is resistant. The key parameters describing a particular combination treatment are its effects on the birth and death rates of cells and the number of point mutations that have the potential to confer resistance. Consider a treatment with two drugs, 1 and 2. We denote by n 1 (respectively, n 2 ) the number of point mutations that have the potential to confer resistance to drug 1 alone (respectively, drug 2 alone). We denote by n 12 the number of point mutations that have the potential to confer resistance to both drug 1 and drug 2 (cross-resistance mutations). We assume that drugs in a combination treatment are given at concentrations tolerable by patients, and define the numbers of resistance mutations (n 1 , n 2 , n 12 ) relative to these concentrations (Katouli and Komarova, 2010).

A crucial quantity for the effects of combination therapy is the expected number, X, of resistant cells at the start of treatment in a lesion containing M cells. From the dynamics of our branching process model (see Supplementary file 1), we obtain

X ≈ M [ n 12 μ + ( n 1 n 2 + n 12 2 ( n 1 + n 2 − n 12 ) ) μ 2 ] .

Here μ = u s log ( Ms ) (log denotes the natural logarithm), where s = 1 − d/b is the survival probability of the branching process initiated with a single cell and u is the point mutation rate, ∼10−9 for most cancers. As µ is small, this formula can be further simplified. If there is at least one possible mutation that could in principle confer resistance to both drugs, n 12 ≥ 1, then X ≈ M n 12 µ. In this case, the expected number of cells resistant to both drugs is independent of the numbers of mutations, n 1 and n 2 , that have the potential to confer resistance to each individual drug. Intuitively, this means that tumor cells are much more likely to become resistant to dual therapy through the occurrence of one mutation conferring resistance to both drugs simultaneously than through sequential mutations conferring resistance to each drug separately. If there is no mutation that could confer resistance to both drugs simultaneously (no cross-resistance), then n 12 = 0 and we obtain X ≈ M n 1 n 2 µ2. This quantity scales with the square of the point mutation rate, so the number of resistant cells in a tumor will be much smaller than for the case n 12 > 0. In general, the expected number of cells resistant to combination therapy with k drugs, with no cross-resistance, is X ≈ M n 1 n 2 … n k µk (proof in Supplementary file 1).

We emphasize, however, that resistance is the outcome of random mutation, division, and death events, and consequently may arise in one lesion but not in another, even if these lesions are otherwise identical. We therefore also obtain formulas for the probability that resistance to combination therapy is present at the time of detection. This probability can be computed as p res = 1 − p 1 p 2 . Here, p 1 is the probability that there is no resistance at detection that arose in a single mutational step, due to one of the n 12 possible cross-resistance mutations. p 2 is the probability that no such resistance arises in two mutational steps. These probabilities can be expressed as follows (proofs in Supplementary file 1, Section 4):

p 1 = exp ( Mun 12 log ( s ) 1 − s )

p 2 ≈ exp [ Mu 2 ( 2 n 1 n 2 + n 12 ( n 1 + n 2 ) ) log ( s ) log ( Ms ) s ( 1 − s ) ] .

As above, s = 1 − d/b is the survival probability of the branching process initiated with a single cell. The quantity 2n 1 n 2 + n 12 (n 1 + n 2 ) in the expression for p 2 represents the number of possible two-step mutational paths to dual resistance.

We turn now to the dynamics of the treatment response. Once treatment starts, sensitive cells decline, but resistant cells continue to grow. We assume that resistant cells maintain the pretreatment birth and death rates, b and d, respectively, during treatment. To obtain estimates for the birth rate b′ and death rate d′ of sensitive cells during treatment, we calculate that the 68 lesions in our dataset declined at median rate b′ − d′ = −0.03 per day (−0.01 and −0.07 being 10th and 90th percentile, respectively). Thus, we set the typical death rate of sensitive cells during treatment to d′ = b′ + 0.03, and additionally explore cases when treatment is less (d′ = b′ + 0.01) or more effective (d′ = b′ + 0.07). As a default in our simulations, we suppose that treatment affects only the death rate (b′ = b), but our mathematical analysis applies also to the case that treatment affects the birth rate.

Figure 2 shows computer simulations of single lesions in response to targeted therapies. Previous studies (Engelman et al., 2007; Corcoran et al., 2010; Diaz et al., 2012; Ellis et al., 2012; Misale et al., 2012; Straussman et al., 2012; Wilson et al., 2012) suggest that about 50 different mutations can confer resistance to a typical targeted therapeutic agent. Assuming that there are 50 or more potential resistance mutations, monotherapy will eventually fail in all lesions that can be detected by conventional imaging (Figure 2A,B) even when the death rate d’ conferred by the therapy is far higher than usually observed in practice (Figure 2C). Small lesions, however, can decrease below the detection limit and appear to be eradicated for years before re-emerging (Figure 2B,C). This result is important, as it explains why tumors can recur after long periods of remission without the need to invoke processes involving cancer stem cells, angiogenesis, or immune escape (Hensel et al., 2012). Note that results similar to those obtained by simulation are observed in several of the individual lesions from actual patients graphed in Figure 1.

Figure 2 Download asset Open asset Tumor response to mono and dual therapy. The tumor grows exponentially until a certain detection size, M, is reached, at which point treatment is initiated. The number of point mutations that could in principle confer resistance to monotherapy is n = 50. For dual therapy, the number of point mutations that could confer resistance to drugs 1 and 2 separately is given by n 1 = 50 and n 2 = 50. The number of point mutations that could confer resistance to both drugs simultaneously is given by n 12 . The point mutation rate was assumed to be u = 10−9 and the rate of cell division b = 0.14 per day and is unaffected by treatment. The rate of cell death before treatment is d = 0.13 per day; it is increased to d’ for sensitive cells during treatment. (A)–(C) For clinically detectable sizes (M = 1010, 109, 108, depending on the location of the tumors and the detection methods used), monotherapy leads to a temporary shrinkage of the tumor but is always followed by tumor regrowth. (D) Due to stochastic fluctuations the few resistant cells present at the start of treatment go extinct and the lesion is eradicated. (E) Treatment leads to a temporary shrinkage of the tumor followed by regrowth. (F) The tumor decreases slowly in response to dual therapy, but resistant cells eventually evolve and cause treatment failure. https://doi.org/10.7554/eLife.00747.005

The results predicted to occur with dual therapy are shown in Figure 2D–F. Here, we also assume that there are 50 mutations that have the potential to confer resistance to either drug alone, but also that there is at least one mutation that can confer resistance to both drugs simultaneously. Intuitively, one might imagine that the existence of even a single cell resistant to both drugs at the start of therapy will automatically result in treatment failure. However, our results show that this is not necessarily true, and that the response depends on the size of the lesion, the number of cross-resistant cells, and the effects of the therapy on the balance between cell birth and cell death. Three examples illustrate these points. In Figure 2D, there is a small number of cells resistant to both drugs at the initiation of dual therapy, but these cells are lost by stochastic drift and the lesion is eradicated. In Figure 2E, there is a greater, but still relatively small number (∼100), of cells resistant to both drugs. The lesion shrinks at first, but eventually progresses due to preexisting cross-resistance mutations within it. In the third lesion, the few cells resistant to both drugs at the initiation of therapy are lost to stochastic drift, but the cytolytic effects of the drug combination are less pronounced than in the other two cases (d’ = 0.15 instead of 0.17 or 0.21). The relatively slow decrease in lesion size enables the generation of de novo resistance mutations during treatment and the lesion eventually recurred (Figure 2F).

In summary, treatment failure can be caused either by the preexistence of resistance to both drugs in a small number of tumor cells (Figure 2E) or the emergence of resistant cells during treatment (Figure 2F). Taking both of these possibilities into account, the probability, p erad , that dual therapy eradicates a lesion containing M cells at the start of treatment is given by

(1) p erad = p 1 ↑ p 1 ↓ p 2 ↑ p 2 ↓ .

p 1 ↑ is the probability that no 1-step resistant lineage arises (and survives) prior to treatment. p 1 ↓ is the probability that no 1-step resistant lineage arises (and survives) during treatment. p 2 ↑ is the probability that no 2-step resistant lineage arises (and survives) prior to treatment. p 2 ↓ is the probability that no 2-step resistant lineage arises (and survives) during treatment. Here, ‘steps’ refers to the number of mutations (one or two) needed to achieve dual resistance, and ‘lineage’ refers to the descendants of a single cell that has achieved dual resistance via a particular mutational path. The therapy is successful if there is no resistant lineage arising in any of these four scenarios; since these are independent events, the overall success probability is obtained by multiplying the corresponding probabilities as shown in equation (1). The probabilities that no 1-step resistant lineages arise before ( p 1 ↑ ) or during treatment ( p 1 ↓ ) and survive are given by Komarova and Wodarz (2005)

p 1 ↑ = exp ( − Mu n 12 )

and Michor et al. (2006)

p 1 ↓ = exp ( Mu n 12 s s ' ) .

Here s = 1 − d/b as above, and s′ = 1 − d′ /b′, where b′ and d′ are birth and death rates of cells sensitive to at least one drug during treatment (note that s’<0). The probabilities that no 2-step resistant lineages arise before ( p 2 ↑ ) or during ( p 2 ↓ ) treatment and survive can be calculated as:

p 2 ↑ = exp [ Mu 2 s ' − s s s ' ( n 1 ( n 2 + n 12 ) log ( 1 s M + u ( n 2 + n 12 ) s ' − s s s ' ) + n 2 ( n 1 + n 12 ) log ( 1 s M + u ( n 1 + n 12 ) s ' − s s s ' ) ) ]

and

p 2 ↓ = exp ( − M u 2 ( 2 n 1 n 2 + n 12 ( n 1 + n 2 ) ) s s ' 2 ) .

The proofs of these results are provided in Supplementary file 1, Section 5. Excellent agreement between equation (1) and simulation results is shown in Figure 3.

Figure 3 Download asset Open asset Probability of tumor eradication for two-drug combination therapy. A single mutation conferring cross-resistance to both drugs (n 12 = 1) can prohibit any hope for a successful dual therapy. Solid curves show analytical results for dual therapy and dashed curve shows analytical results for a typical monotherapy, both are calculated using equation (1). Markers (square, triangle, circle, diamond) indicate simulation results (averages of 106 runs). Parameter values: birth rate b = 0.14, death rate d = 0.13, death rate for sensitive cells during treatment d’ = 0.17, point mutation rate u = 10−9. https://doi.org/10.7554/eLife.00747.006

Although modeling of single neoplastic lesions is the norm in theoretical studies, most patients with advanced cancers have multiple lesions and curing a patient requires eradication of all lesions. Equation (1) can be used to evaluate which combination treatments will be successful in typical patients with multiple metastatic lesions.

To determine the total extent of disease in typical patients who enroll for clinical trials, we quantified all radiographically detectable metastases in 22 such patients: 7 with pancreatic ductal adenocarcinomas, 11 with colorectal carcinomas, and 6 with melanomas—a different cohort than that depicted in Figure 1, in which only index lesions (those easiest to measure) were evaluated. The number of metastatic lesions in the 22 patients described in Table 1 ranged from 1 to 30, and their total tumor burden ranged from 9 × 108 to 3 × 1011 cells (see Supplementary file 2).

Table 1 Probability of treatment failure for combination therapy in patients https://doi.org/10.7554/eLife.00747.007 Patient Primary tumor type Number of metastases Total tumor burden (number of cells) Probability of treatment failure Monotherapy Dual therapy: n 12 = 1 Dual therapy: n 12 = 0 N1 Pancreas 18 2.6 × 1011 1 1 0.283 N2 Colon 25 2.3 × 1011 1 1 0.26 N3 Melanoma 26 1.7 × 1011 1 1 0.203 N4 Melanoma 30 1.4 × 1011 1 1 0.172 N5 Colon 21 1.0 × 1011 1 1 0.128 N6 Melanoma 8 9.8 × 1010 1 1 0.12 N7 Colon 25 9.1 × 1010 1 1 0.112 N8 Pancreas 8 7.4 × 1010 1 1 0.092 N9 Pancreas 23 6.4 × 1010 1 1 0.08 N10 Pancreas 5 5.5 × 1010 1 1 0.069 N11 Colon 14 5.4 × 1010 1 1 0.068 N12 Rectal 23 4.8 × 1010 1 1 0.061 N13 Melanoma 9 4.1 × 1010 1 1 0.052 N14 Pancreas 13 4.1 × 1010 1 1 0.051 N15 Pancreas 8 3.3 × 1010 1 1 0.042 N16 Melanoma 7 2.2 × 1010 1 1 0.028 N17 Melanoma 10 2.1 × 1010 1 1 0.027 N18 Colon 4 2.0 × 1010 1 1 0.026 N19 Melanoma 9 1.8 × 1010 1 1 0.023 N20 Colon 3 1.6 × 109 1 0.881 0.002 N21 Melanoma 21 1.3 × 109 1 0.828 0.002 N22 Pancreas 1 8.5 × 108 1 0.677 0.001

For each of these 22 patients, we used equation (1) to calculate the probability that monotherapy or dual therapy would eradicate all the patients’ lesions. We find that monotherapy will fail in all 22 patients (Table 1 and Supplementary file 3), as expected from the simulations in Figure 2A–C and from clinical experience. If there is even one possible mutation that can in principle confer resistance to both drugs, then our model shows that dual therapy has also only a small chance of curing patients, even those with the smallest tumor burden. In our cohort of 22 patients, none are expected to be cured under these circumstances (Table 1). Only if there are no potential mutations that can confer cross-resistance will dual therapy be successful in eradicating all lesions. In the cohort described in Table 1, we calculate that eight patients (those with the smallest tumor burden) would have >95% probability of cure. Those with the largest tumor burden would still have a >20% probability of tumor recurrence. Additional simulations show that therapy with three agents will also not cure patients if there is even one mutation that can confer resistance to all three agents. Similar conclusions hold if we vary parameter values within a reasonable range (Supplementary file 3). We note that in patients whose tumors have high cell turnover (time between cell divisions of 1 day, corresponding to b = 1), even dual therapy with no cross-resistance mutations would be expected to fail in 37% of patients described in Table 1 (Supplementary file 3).

Graphical representations of the simulated responses of two patients with multiple metastatic lesions are shown in Figure 4. With monotherapy in patient N1 (Figure 4A), all lesions are predicted to regress, but then recur within a year or so after the initiation of therapy (Figure 4B, left panel). Treatment failure in most lesions would also occur after dual therapy when there is at least one mutation that could confer resistance to both agents, although the length of remission will be longer than with monotherapy (Figure 4B, middle panel). In patient N11, with less disease burden, dual therapy will fail to eradicate several of the lesions when there is a possibility of a single cross-resistance mutation, but there is hope of cure if no such cross-resistance mutations are possible (Figure 4C,D).

Figure 4 Download asset Open asset Treatment response dynamics to monotherapy and dual therapy in two patients. (A) Depiction of all 18 detectable metastases in patient N1, who had a particularly heavy tumor burden (scale 1:4). (B) Simulated treatment of patient N1, comparing monotherapy with n = 50 resistance mutations and dual therapy with n 1 = n 2 = 50 resistance mutations to the individual drugs and one (n 12 = 1) or no (n 12 = 0) cross-resistance mutations to both drugs. (C) Depiction of all 14 detectable metastases in patient N11, who had a more typical tumor burden (scale 1:4). (D) Simulated treatment of patient N11. Parameter values for simulations in (B) and (D): birth rate b = 0.14; death rate d = 0.13; death rate for sensitive cells during treatment d′ = 0.17; point mutation rate u = 10−9. https://doi.org/10.7554/eLife.00747.008

In current clinical practice, it is common to administer targeted agents sequentially: once relapse occurs, a second, often experimental, agent is administered. The model described above can also be used to predict the relative effectiveness of sequential vs simultaneous therapies of a single lesion with two drugs. When there is a possibility of a single mutation conferring resistance to both drugs, sequential combination therapy will ‘always’ fail. In ∼74% of lesions, the failure is due to mutations that were present prior to the treatment with the first drug, whereas in ∼26% of the lesions, failure is due to the development of cells resistant to drug 2 during treatment with drug 1 (Figure 5A and Figure 5—figure supplement 1). With simultaneous treatment, it is possible to eradicate ∼26% of the lesions even when cross-resistance mutations are possible (Figure 5B). When there is no possibility of a mutation conferring cross-resistance to both drugs, the differences are even more striking: sequential therapy fails in 100% of cases (Figure 5C), whereas simultaneous therapy succeeds in >99% of lesions of the identical size (Figure 5D).

Figure 5 with 1 supplement with 1 supplement see all Download asset Open asset Sequential vs simultaneous therapy with two drugs. (A) If there is even a single mutation that confers cross-resistance to both drugs (n 12 = 1), then sequential therapy will fail in all cases. In 73.7% of the cases, this failure is due to the exponential growth of fully resistant cells that were present at the start of treatment. In the remaining 26.3% of cases, the failure is due to resistance mutations that developed during therapy with the first drug. (B) With simultaneous therapy, 26.3% of patients can be cured under the same circumstances. In the remaining patients (73.7%), cross-resistant mutations existed prior to the therapy and their expansive growth will ensure treatment failure whether treatment is simultaneous or sequential (see Figure 5—figure supplement 1 for further details). (C) and (D) If the two drugs have no resistance mutations in common (n 12 = 0), then simultaneous therapy is successful with a probability of 99.9% while sequential therapy still fails in all cases. https://doi.org/10.7554/eLife.00747.009

One of the most important aspects of the cancer stem cell hypothesis revolves around therapeutic resistance. Evidence to date does not indicate that cancer stem cells are innately resistant to either single drugs or drug combinations. However, the precise proportion of cancer stem cells (among all cancer cells) has a dramatic effect on the development of resistance. This effect can be studied using our model if we use the number of cancer stem cells as an effective population size in our formulas and adjust other parameters to account for the stem cell dynamics (Tomasetti and Levy, 2010) (i.e., the birth rate should correspond to the rate of symmetric renewal, the rate of symmetric differentiation should be added to the death rate, and an effective mutation rate for stem cells should be introduced to account for mutations that occur during asymmetric division). For example, if cancer stem cells represent only 0.1% of cancer cells, then the development of resistance to single agents or combinations is roughly 0.1% as likely as if 100% of the cancer cells have the capacity to repopulate the tumor. The fraction of cancer stem cells appears to be this low in CML, perhaps explaining the remarkable success of imatinib (Michor et al., 2005). In solid tumors, however, the fraction of cancer stem cells seems much higher, usually higher than 5% and in some cases close to 100% (Shackleton et al., 2009). This issue is further complicated by the fact that the situation is plastic, with non-stem cells converting to cancer stem cells under certain conditions (Gupta et al., 2009). As better approaches to quantify cancer stem cells in solid tumors become available, our estimates of the likelihood of therapeutic success will be improved.

If resistance has a fitness cost, then we expect a smaller number of resistant cells at the start of treatment and correspondingly a higher chance of treatment success. We used computer simulations to verify our results in the case when there is a cost for resistance, by assuming that each resistance mutation decreases the net growth rate of the cell by up to 10%. The results are shown in Table 2. For combination therapies with drugs that have resistance mutations in common, the probability of eradicating a lesion is only marginally affected by costly resistance. For dual therapies with no cross-resistance mutations, treatment has a high chance of eradicating all but the largest lesions, whether or not resistance is costly. In the case of large lesions with high cell turnover rates (the case in which even dual therapies with no cross-resistance might fail), costly resistance increases the chance of treatment success. For example, if each resistance mutation decreases the net growth rate of cells by 10%, the probability that dual therapy with no cross-resistance mutations will eradicate a lesion of size 1011 in which cells divide on average every day is 68% (compared with 47% in the case of neutral resistance).

Table 2 Simulation results for the probability of treatment failure when resistance is costly https://doi.org/10.7554/eLife.00747.011 Dual therapy: Number of cells Birth rate Probability of treatment failure n 1 = n 2 n 12 c = 0% c = 1% c = 5% c = 10% 50 0 109 0.14 0.0 0.0 0.0 0.0 50 0 109 1 0.01 0.01 0.01 0.0 50 1 109 0.14 0.74 0.73 0.72 0.7 50 1 109 1 0.74 0.74 0.72 0.7 50 0 1011 0.14 0.12 0.11 0.08 0.06 50 0 1011 1 0.53 0.51 0.42 0.32 50 1 1011 0.14 1.0 1.0 1.0 1.0 50 1 1011 1 1.0 1.0 1.0 1.0

Some therapies may directly eliminate tumor cells (d′ > d), whereas others may impede their division (b′ < b). Our formulas account for both of these possibilities. Overall, the rate b′ − d′ of tumor decline is of primary importance, and whether this is achieved by eliminating cells or suppressing division has only a minor effect on treatment outcomes. For example, consider a dual therapy with n 1 = n 2 = 50, n 12 = 1, applied to a lesion of size M = 109, with other parameters as inferred from our dataset. If this therapy shrinks the tumor at rate −0.03 per day by increasing cell death, the eradication probability is 26%. If the therapy instead suppresses division, this probability increases to 29%, because there are fewer chances for resistance mutations during treatment.

While our typical parameter values are derived from the melanoma dataset, our analytical results can accommodate parameter values from any other type of cancer, once they become available. Furthermore, our results are qualitatively robust across a wide range of birth and death rates (Supplementary file 3). The parameters with the strongest effects on the success of combination treatments—apart from the number of cross-resistance mutations—are lesion sizes and point mutation rate. Thus, we expect that combination treatments will be more effective in cancers with small fractions of tumor stem cells (small effective population size of lesions) and less effective in cancers with significantly increased point mutation rates.