Combination Therapy and the Evolution of Resistance: The Theoretical Merits of Synergism and Antagonism in Cancer

Modern cancer therapeutics have excellent initial efficacy, but resistance often develops in a span of months. Investigating combination therapy for combating cancer resistance is currently of great interest in the clinical setting (1, 2), at the bench (3, 4), and in computational modeling (5–9). The combined effect of two drugs can be categorized as synergistic, additive, or antagonistic, depending on whether it is greater than, equal to, or less than the sum of the drugs' individual effects (Fig. 1A; ref. 10). As synergism by definition has the greatest potency relative to total dose, and as toxicity often increases monotonically with dose, much focus has been given to finding potent synergistic combinations. However, combination toxicity is complex (11–13), and the side effects of multidrug treatments remain speculative.

Antagonism is sometimes misconstrued as an “antidote” effect, where one drug cancels out the efficacy of another (a phenomenon called super-antagonism). Super-antagonism is counter-productive, but antagonism is not; “less-than-additive” merely implies that the second drug gives smaller additional benefit than in the additive case, but is still beneficial. In the field of infectious diseases, work by Kishony and colleagues (18–21) discovered that antagonistic antibiotic combinations could delay the development of bacterial resistance. Can this concept be applied to cancer (22) as a strategy for delaying resistance evolution?

Early work by Nowell (23) framed the development of cancer drug resistance as a process of mutation and evolutionary selection. Theoretical simulations of the Darwinian dynamics of drug-sensitive and -resistant subclones in heterogeneous cancers have described how evolutionary trade-offs (24, 25), aggregation effects (26), variable cell densities (26–28), or spatial interactions (29) can create differential selective pressures among subclones and influence tumor evolution in response to therapy. These insights from mathematical oncology have inspired the design of “evolutionarily-enlightened therapies” (26), which consider factors such as future states of resistance (6, 30), evolutionary trade-offs (31), and temporal subclone vulnerabilities (32), and predict optimal scheduling to guide clinical studies (33–35). Such models have also been used to study combination therapies for cancer (5–9), but more theoretical evolution work is needed to understand the long-term impact of synergistic (vs. nonsynergistic) therapy. A combination has more efficacy than monotherapy due to the simultaneous actions of both drugs. However, if some cells develop resistance to one drug, they will escape not only the effect of the one drug, but also its enhancement/masking of the second drug. Hence, we ask how the efficacies of different combinations would be vulnerable to the development of single-drug resistance.

In this work, we investigate the theoretical effect of two-drug combination therapies on the evolutionary dynamics of resistance in a tumor cell population. We ask how Kishony's discoveries about the risk of synergistic combinations during antimicrobial therapy might apply to cancer. As there are many types of cancer and anticancer treatments, we abstracted the broad landscape using binary resistance (5–9) to define four subpopulations of cells: fully-sensitive (SS), resistant to drug 1 (RS), resistant to drug 2 (SR), and doubly-resistant (RR). We simulated the number of cells over time using a simple nonspatial model of cellular alteration and proliferation. The fitness of each phenotype is quantified by adapting the concept of dose equivalence (14, 17) to Greco and colleagues' (16) response surface model for two-drug combinations. We establish two comparison methods to construct a fair comparison between different classes of treatments: the Constant-Dose Method uses dose as its basis of comparison, whereas the Constant-Efficacy Method uses efficacy on fully-sensitive cells as its basis of comparison. Our findings may provide a conceptual framework to guide future experiments in specific cancer systems.

We developed a stochastic computational model with the following features: a simple exponential process representing tumor growth; drug-dependent cell fitness parameters representing cellular effects of treatment; and stochastic introduction of resistance phenotypes, based on the preliminary model in ref. 7. These features were implemented using a time-step simulation of proliferation, death, and phenotype alteration.

Between the extremes of full sensitivity and double resistance, certain types of partial resistance during combination therapy are anticipated to occur, especially single-drug resistance (36). Hence, we simulated the tumor dynamics of four subpopulations changing independently over time: cells that are sensitive to both drugs (SS), resistant to drug 1 (RS), resistant to drug 2 (SR), and resistant to both drugs (RR). This categorization followed the approach of Coldman and Goldie (9), which has also been adapted by others (6, 37). Each subpopulation represents a phenotype class of cells that may be genetically or metabolically diverse, but categorized for the presence/absence of resistance. Instead of tracking the birth and death of individual cells (5, 8, 9, 38), our model tracked the growth and death of the subpopulations (6), each following exponential dynamics. Exponential growth and exponential decay provide a first-order approximation of observed population sizes according to empirical studies of solid and liquid tumors (39).

At each time step, each cell has a small probability of switching into a different subpopulation, according to a Markov transition model (Fig. 1B, “Alteration phase”). Sensitive cells could acquire resistance to any one drug; subsequently, singly-resistant cells could acquire resistance to the second drug. Hence, a fully-sensitive cell would have to acquire two independent alterations to become doubly-resistant, making this model incapable of describing cases where resistance opportunities are nonindependent. We permitted state transitions to be bidirectional (details in Supplementary Text S1), assuming all alterations to be independent. We set the alteration rate |\mu $| to be a classically cited rate of human gene mutation, 10⁻⁶ (40), bearing in mind the complexities not covered, such as nonmutational alterations (e.g., epigenetics), cancers with orders-of-magnitude higher/lower mutation rates (41), and mutations in loci with dominant, recessive, or mutator genes (38).

We abstracted cell fitness as two holistic parameters, proliferation and apoptosis, rather than parameterizing molecular mechanisms of fitness (42, 43). The proliferation probability of subpopulation |i$|⁠, denoted |P{P_i}$|⁠, signifies the proliferative potential of cells in the subpopulation under the given therapeutic condition. When exposed to treatment, |P{P_i}$| decreases from 1 (full proliferative capacity, e.g., untreated cells) to any fraction or 0 (nonproliferative). Meanwhile, the apoptosis potential of subpopulation |i$| was represented by the sum of the cell's basal apoptosis rate, |{a_b}$|⁠, and the drug-induced apoptosis rate, |{a_i}$|⁠. We excluded trivial trials where cancer always decreased or always increased regardless of treatment, by focusing on cases where SS shrank and RR grew under treatment.

Growth and death were simulated using a cycle of proliferation, death, and phenotype alteration (Fig. 1B; Eq. A). Although resistance can sometimes be accompanied by a phenotypic cost (24), it can also be accompanied by fitness advantages such as dedifferentiation and increased aggressiveness (44). Our model assumed that resistance conferred neither fitness cost nor advantage, meaning that the subpopulations had uniform proliferation rate |p$| and basal apoptosis rate |{a_b}$| without treatment. During treatment, the proliferation rate would be scaled down by each subpopulation's proliferation probability |P{P_i}$|⁠, conferring differential fitness across the subpopulations, summarized in the effective proliferation rate ^prolif|{k_i}$| and effective apoptosis rate ^death|{k_i}$| (Eq. A). Our main simulations started with 10⁹SS cells (representing a minimum size for detection; ref. 45) and terminated when tumor size reached zero (representing eradication) or 10¹² cells (representing a lethal tumor burden; ref. 45). Simulations with pre-existing resistance in the starting population are also considered. The Supplementary MATLAB file provides the simulation code.

Anticancer therapies can be antiproliferative, proapoptotic, or both (46, 47). Our model employed two drug-effect parameters to quantify these effects: reduction in proliferation probability |P{P_i}$|⁠, and drug-induced apoptosis rate |{a_i}$|⁠.

A central data structure of the method is the table of fitness parameters for the different subpopulations, illustrated in Table 1 using the fitness parameter |P{P_i}$|⁠. Under monotherapy with drug 1, the proliferation probabilities of cells resistant to drug 1 (⁠|P{P_{RS}}$| and |P{P_{RR}}$|⁠) remain approximately unchanged, whereas |P{P_{SS}}$| and |P{P_{SR}}$| decline to x < 1. Conversely, under drug 2, |P{P_{SR}}$| and |P{P_{RR}}$| are approximately unchanged, whereas |P{P_{SS}}$| and |P{P_{RS}}$| decrease to y < 1. When both drugs are given in combination, the combined effect is less straightforward to determine, contingent upon the type of drug interaction.

To generalize the parameterization of proliferation and apoptosis beyond specific cases, instead of using experimental data (6, 8, 30), we calculated the theoretical combination effect using a method based on the response surface model by Greco and colleagues (16). We applied the concept of “dose equivalence” (ref. 17; assumed in Loewe's additivity model, ref. 14) to define the dose of one drug that generates an effect equivalent to a combination (Eq. B, derivation in Supplementary Text S2). Greco's model designates an interaction parameter |\alpha $| to denote the strength of interaction between a combination, where a positive, zero, or negative |\alpha $| signifies synergism, additivity, or antagonism, respectively. For flexibility in defining whether anticancer drugs target proliferation or survival, our model introduced two modulating parameters: |P{P_c}$| for describing the maximum effect on proliferation, and |{a_c}$| for describing the maximum effect on apoptosis (Eqs. C and D). A larger |P{P_c}$| or |{a_c}$| indicates a greater effect on proliferation or apoptosis, respectively. For our main simulations, we simulated both effects at once, assuming |P{P_c} = 1$| and |{a_c}\ \, \sim\ \!\!\!0.254$|⁠. (Proliferation-specific or apoptosis-specific effects are shown in Supplementary Text S3 and Supplementary Fig. S1.)

Our model assumed resistance to be a binary effect, defined as the ability to “ignore” 90% of the dose (e.g., |{S_{1,RS}} = 10\% $| in Eq. B). This definition is an abstract simplification for any molecular mechanism of resistance, simply lowering the dose–response curve. Binary resistance is a simplified discretization of actual resistance mechanisms that are often graded (e.g., expression of drug-efflux pumps; ref. 48), a promising topic for future modeling.

In combination therapy, determining maximum tolerable doses (MTD) can be complex. Evidence suggests that combination toxicity can be “nonadditive” relative to the individual toxicities; drugs with overlapping toxicity may cause side effects at doses lower than the MTD of either drug (11), whereas antagonistic combinations sometimes generate little increase in adverse events compared with monotherapy (12). The complexity of combination toxicities (and MTD) creates uncertainty for establishing a fair dosing method for comparing alternative combination therapies.

Therefore, we defined two complementary methods of comparison: the Constant-Dose Method and the Constant-Efficacy Method. The Constant-Dose Method defines fairness using dose as the equalizer—all combinations were dosed using the same combination of effective doses, regardless of synergism or antagonism. Effective dose (ED_X) refers to the dose of a single drug that affects X% of the population. Hence, two drugs may be delivered at different absolute concentrations yet the same effective dose (e.g., if drug X is delivered at ED_50,X = 1 mg/kg, while drug Y is delivered at ED_50,Y = 0.5 mg/kg). Under this method, efficacy increases from antagonistic, additive, to synergistic.

The Constant-Efficacy Method defined fairness by using combined efficacy as the organizing equalizer. Doses were set so that all combination treatments had equal efficacy toward the fully-sensitive tumor bulk (SS cells), including the possibility that the administered ED_X differs across combinations. This method assumes that antagonistic drugs can be delivered at doses higher than synergistic drugs without violating safety limits (for example, if toxicity mirrors efficacy).

As comparative measures for the effectiveness of resistance suppression, we evaluated how different combination therapies affected |{T_{{\rm{lethal}}}}$| (time to reach 10¹² cells) and |{T_{RR}}$| (time until doubly-resistant RR cells arise). The simulation parameters were calculated with the described approach for the constant-dose and constant-efficacy comparisons (Supplementary Fig. S2), where each combination was dosed symmetrically (i.e., with the same ED_X).

To generalize how synergism and antagonism affect evolution, we derived an analytical approximation of |{T_{RR}}$| (Eq. S16 of Supplementary Text S4.1). The full derivation allows arbitrary level of pre-existing resistance and arbitrary dosing. If we set dosing to be symmetric (Eq. E) and set pre-existing resistance to zero, then the full analytical approximation (Equation S16) can be simplified as follows. The probability that double resistance will arise (for the first time) at generation |t$| is

Here, |{k_{RS \,+ \,SR\ }}$| is the net growth rate of singly-resistant subpopulations RS and SR collectively, |{k_{SS}}$| is the net growth rate of subpopulation SS, |\mu $| is the alteration rate, and |{n_{SS,0\ }}$| is the number of sensitive SS cells at time 0. In interpreting the approximation, we focus on |{k_{RS\, + \,SR\ }}$|> 1, because then some singly-resistant cells will continue to exist, permitting the question of how long until at least one transforms into RR. Note that there are critical points when |{k_{SS}}$| or |{k_{RS \,+ \,SR\ }}$|approach 1 (see the denominator of the exponent).

Monte Carlo simulations (n = 10,000) were run for synergistic, additive, and antagonistic treatments (using additive treatment as control) according to the Constant-Dose Method, dosing the combinations at a constant combination of ED_X. Figure 2 shows five individual simulations. Looking at the additive control, tumor response exhibits three stages. Firstly, the SS subpopulation (indigo) goes down, causing the total tumor mass (yellow) to decline (“tumor regression stage”). Next, the singly-resistant subpopulations SR (red) and RS (magenta, eclipsed by red) arise and proliferate (“resistance evolution stage”). Finally, doubly-resistant cells (cyan) arise at approximately 27.3 generations, eventually causing catastrophically fast growth (“tumor relapse stage”). This three-stage tumor response is common to many of our outcomes. The same high-level conclusion was also achieved in ref. 6.

In the constant-dose simulations, synergism caused a steeper initial decrease of the SS curve compared with antagonism (Fig. 2A). Because SS cells comprise the majority of the population, the total tumor mass shrank more rapidly in response to treatment. Meanwhile, antagonistic combinations had weaker efficacy against the SS subpopulation and decreased its size slowly. Consequently, there were more cells per time-step with opportunity to transition toward resistance (compared with in the synergistic case), making RR cells emerge sooner. Therefore, not only are antagonistic combinations less effective against the naïve cells, but by letting a large tumor population persist, they also allow resistance to arise more quickly. Once the untreatable RR cells appeared, the clonal expansion of these highly fit cells followed, as seen in the steep slope of the RR curve. RR cells soon made up most of the total population and expanded the population size.

Simulations were again run for synergistic and antagonistic treatments according to the Constant-Efficacy Method, dosing each drug pair with symmetric ED_X to achieve |{E_{SS}}$| of 0.69 (the same as in the additive control). Interestingly, these results showed that as treatment became more antagonistic, |{T_{RR}}$| and |{T_{{\rm{lethal}}}}$|were prolonged, giving better outcome. Because all combinations have equal efficacy toward SS cells, the simulations showed similar initial responsiveness, meaning the same downward slopes of the SS curves in Fig. 2B. However, the singly-resistant subpopulations, SR and RS, quickly went up in the synergistic case, whereas this increase was not as strong in the antagonistic case. Because the singly-resistant subpopulations expanded quickly in the synergistic case, tumor relapse occurred earlier than in the antagonistic case (after ∼75 generations in the synergistic plot vs. after ∼90 generations in the antagonistic plot). The earlier expansion of singly-resistant subpopulations, which are one step away from RR, promoted the emergence of doubly-resistant cells. The fast increase of the RR curve then accelerated |{T_{{\rm{lethal}}}}$|⁠.

To test whether the above observations were consistent, we repeated the stochastic simulations (n = 10,000) for a range of α values, according to the Constant-Dose and the Constant-Efficacy Methods. To confirm the general nature of the findings, we also used the analytical equation for |{T_{RR}}$| (see Materials and Methods). Figure 3A–E validates the applicability of our analytical equation, showing that the analytical approximation of |{T_{RR}}$| agrees closely with the results from Monte Carlo simulations.

Stochastic simulations showed that under the Constant-Dose Method, increasing synergism prolongs |{T_{RR}}$| and |{T_{{\rm{lethal}}}}$| (Fig. 4A). This advantage of synergism was also shown by the analytical approximation of |{T_{RR}}$|⁠, by varying |{k_{SS}}$| while fixing |{k_{RS + SR}}$|⁠. Lowering |{k_{SS}}$| means increasing synergism, which shifts the probability distribution toward higher |{T_{RR}}$| (i.e., better anticancer outcomes; Fig. 4B). In contrast, stochastic simulations showed that under the Constant-Efficacy Method, |{T_{RR}}$| and |{T_{{\rm{lethal}}}}$| increased as the treatment became more antagonistic (Fig. 4C). This finding was confirmed by the analytical equation for |{T_{RR}}$|⁠, by varying |{k_{RS + SR}}$| while holding |{k_{SS}}\ $|fixed. Lowering |{k_{RS + SR}}$| means increasing the antagonism of the combination therapy, which broadens the probability distribution of |{T_{RR}}$|⁠, peaking at significantly higher |{T_{RR}}$| (Fig. 4D).

Given the prevalence of pre-existing resistance in many cancers, we asked if the merits of synergism in constant-dose scenarios and antagonism in constant-efficacy scenarios would change if some singly-resistant cells were present in the starting population. When we added some pre-existing singly-resistant cells, the relative merits of synergism or antagonism were unchanged (Fig. 5), although the absolute times to RR (⁠|{T_{RR}}$|⁠) were globally shorter. As the number of pre-existing singly-resistant cells increased, there was a threshold above which the choice of synergism or antagonism made little difference, because RR cells would quickly evolve and drive tumor relapse. This result was not limited to cases where RS and SR cells were present in equal numbers, as shown by randomly varying the levels of pre-existing singly-resistant RS and SR cells (Supplementary Fig. S3).

In this work, we simulated a simplified cancer cell population to understand how different drug interaction types (i.e., synergistic, additive, antagonistic) affected the long-term evolution of resistance. We modeled the growth trends of four cell subpopulations (fully-sensitive, two types of singly-resistant, and doubly-resistant) responding to two-drug combination treatments, and we performed side-by-side comparisons of synergistic versus nonsynergistic therapies, according to the Constant-Dose or the Constant-Efficacy Method. Most simulations displayed three stages of tumor response: regression, resistance evolution, and relapse.

The two comparison methods gave opposite results, reflecting that these paradigms are opposite extremes of the plausible dosing spectrum. Under the Constant-Dose Method, synergism was more effective than antagonism at reducing cell numbers in the short term, and suppressing the growth of resistant cells in the long term (Fig. 2A). In contrast, under the Constant-Efficacy Method, antagonism was better at suppressing the expansion of resistant cells (Fig. 2B). As previously implied in refs. 5, 6, 8, 9, |{T_{RR}}$| is a surrogate measure for the evolutionary success of a combination treatment, because delaying the emergence of doubly-resistant cells prolongs the time until acquiring a lethal tumor burden (Fig. 4A and C). Our model, by construction, gave nearly equal fitness to RR cells regardless of dosing method (RR cells ignored 90% of any dose). Because the RR subpopulation would grow exponentially after arising, the number of generations between |{T_{RR}}$| and |{T_{{\rm{lethal}}}}$| would be approximately constant. (Possible exceptions include therapy-failure cases where both drugs are weak against SS cells, allowing the expansion of non-RR cells to drive relapse.) Looking at the analytical approximation for |{T_{RR}}$|⁠, the definition of success in the Constant-Dose and Constant-Efficacy Methods becomes straightforward: by substituting the appropriate |{k_{SS}}$| and |{k_{RS + SR}}$| values, synergism always wins in the constant-dose scenario, and antagonism always wins in the constant-efficacy scenario. What this means in practice depends on MTDs, which are more complex than the one-dimensional spectrum between Constant-Dose and Constant-Efficacy Methods.

In a constant-dose comparison, synergistic combinations are superior because they strongly impair the fitness of sensitive cells. The high efficacy of synergistic combinations enables rapid immediate killing of sensitive cells, shrinking the tumor cell population during early-stage treatment. In theory, this initial success decreases the number of opportunities (and slows the speed) for resistance to evolve. This result implies that when a treatment has a significantly better initial efficacy than its alternative, it may be able to generate a better long-term outcome, even if this efficacy is achieved via synergism. Meanwhile, the results of the constant-efficacy comparison may seem counterintuitive, but they can be explained (below) by studying how synergism and antagonism give differential fitness rewards to singly-resistant cells, which can be understood by considering the inverse symmetry between combination efficacy and evolutionary fitness.

The superiority of antagonism under the constant-efficacy comparison is intrinsic to the definitions of synergism and antagonism. When synergistic and antagonistic combinations are compared by dose (e.g., Constant-Dose Method, Fig. 6A), the increase in efficacy from single-drug to two-drug treatment, by definition, is larger for synergistic drugs (red bars) than antagonistic (blue bars). Meanwhile, treatment efficacy, regardless of drug mechanisms, is determined by the impairment of the cancer, which is measured in our model as a decreased evolutionary fitness. This inverse symmetry between efficacy and fitness means that the difference in fitness as a result of a synergistic or antagonistic treatment is also inherent to the definition of synergism or antagonism—in this constant-dose case, synergistic treatments decrease fitness more significantly than antagonistic treatments.

However, the impact of resistance can be abstractly described as a reduced sensitivity to drug dose, making full sensitivity resemble a two-drug treatment and single resistance resemble a single-drug treatment. This means that, as cells develop single resistance and dose sensitivity decreases, the cells will gain more fitness under synergistic treatments than under antagonistic treatments, implying that synergism intrinsically gives fitness rewards to singly-resistant cells. This is also evident when we view resistance as a spectrum. We defined a fitness function (Fig. 6B) that summarizes the cancer phenotype as a function of fractional resistance. Plotting fitness versus resistance (an “evolution plot”) shows the relative reward from synergism or antagonism (Δfitness/Δresistance) as cells develop resistance. Δfitness/Δresistance is intrinsically greater for synergistic than for antagonistic treatments.

Applying this concept to the Constant-Efficacy Method, which considers an alternative dosing scenario where the combinations have equal efficacies toward sensitive wild-type cells (equal efficacies for Drugs 1+2, Fig. 6C), the evolution plot (Fig. 6D) shows that although both treatments give equal fitness to fully-sensitive cells, synergistic treatments will reward partial resistance with higher fitness than antagonistic treatments, because of the greater Δfitness/Δresistance of synergism. To explain this observation, a synergistic combination needs both drugs working simultaneously for their benefits to ensue. If either drug is evaded, or if the sensitivity of the cells to the full drug doses lessens, not only will the single-drug effect diminish, but the enhancement of efficacy from the combination effect will also be lost. Because the enhancement of efficacy from synergism is greater than that from antagonism, synergism will cause a bigger reduction of efficacy when cells develop resistance, thereby rewarding partial resistance with a much higher fitness. This makes synergistic strategies more fragile to the development of partial resistance. Conversely, by not rewarding partially-resistant cells with much fitness advantage, antagonism will slow down the expansion of partially-resistant subclones. The cautionary lesson is that if two combinations have the same initial efficacy, and if one combination achieves its efficacy through synergism, then the less synergistic option would be preferable, because its efficacy would exhibit less deterioration against a subpopulation with single-drug resistance. Real-world dosing is complex and may lie along a spectrum between the absolute ideals of constant dose and constant efficacy. Supplementary Fig. S4 illustrates a potential hybrid.

The parallel between efficacy and fitness suggests an insight into the contrast between offensive and defensive strategies to combat resistance evolution. Defense refers to preparing to counter anticipated future attacks, such as foreseeable forms of single-drug resistance or partial resistance. We assert that synergism is an offensive (proefficacy) strategy, because resistance is suppressed by decimating cancer cells before resistance can emerge. Meanwhile, antagonism is a defensive (antiresistance) strategy; although it may not be optimally effective at reducing cell numbers initially, it suppresses singly-resistant subpopulations by not giving them as much fitness advantage as synergism might. Our findings illustrate that, in some cases, treatments that do not yield better initial response can give a better long-term outcome, because of their effectiveness in suppressing the evolution of resistance. Therefore, the idea of a defensive strategy is aligned with multiple lines of research (7, 25, 32, 49) that advocate a shift in therapeutic approach, toward “treat-to-stabilize” versus “treat-to-eradicate” (49).

The qualitative difference in how combinations reward different resistance phenotypes is independent of the mathematical formalism for defining synergism and antagonism, but it is affected strongly by the formalism for MTD. Dosing moves the evolution curves up and down along the fitness axis, while not changing the curve shape significantly. Therefore, the advisability of the proefficacy or antiresistance strategy depends on the dosing of the comparisons. Nonetheless, the allowable doses of different combinations are highly variable because dosing is limited by the toxicity of the combination (not the sum of the individual toxicities). In addition, the effect of a combination may switch between synergism and antagonism depending on dose (16, 50), context (50), or tumor progression. Although analyses over a variety of parameters showed that our conclusions were robust (Fig. 5; Supplementary Figs. S1, S3, S5, and S6), the variability of dosing strategies and the heterogeneity of cancers introduce uncertainty, and it remains to be determined what combination will be superior for treating a particular instance. One might be tempted to interpret our results as recommending synergistic drugs for naïve tumors and nonsynergistic drugs for tumors with partial resistance or fast mutation, but empirical testing would be needed to determine which strategy could be applied for each specific case.

There are several caveats for interpreting our results. Firstly, the state-transition probability is independent for each cell, so the time required to create the first doubly-resistant cell depends powerfully on the number of singly-resistant cells. In reality, some forms of resistance might not arise in proportion to the number of precursor cells involved (e.g., cancer-associated macrophages or fibroblasts). Secondly, the assumption that there are no direct transitions from fully-sensitive (SS) to doubly-resistant (RR) would be unreasonable for combinations where resistance to one drug provides some intrinsic resistance to the other drug (e.g., if the drugs target the same pathway). Such cases violate our model's independence assumption, and regardless how a biological system might achieve a direct leap from SS to RR, such a transition would have dire effects on |{T_{RR}}$| (Supplementary Fig. S7). Our model also assumes a uniform phenotype-transition rate of 10⁻⁶, but mutation rates may vary in vivo (41), and this could alter the outcome. Future work should model specific resistance mechanisms. In the absence of that, we simulated varying rates of alteration (Supplementary Fig. S5), which yielded the following observations: the relative benefits of synergism and antagonism were consistent across a range of alteration rates; the benefits of antagonism were asymmetrically better than the benefits of synergism for slow-evolving tumors; and the benefits of either strategy were lost in fast-evolving tumors.

Our model simplifies the diversity of resistance mechanisms by using phenotypic categories for single- or multidrug resistance. This ignores several issues: the presence of hierarchical heterogeneity (51); the fitness heterogeneity arising from cooperation and cheating (29); and the fitness penalties associated with drug resistance (24). We simulated cases where drug resistance was associated with a fitness cost (or advantage). When drug resistance incurred a fitness cost, the relative benefits of synergism under constant dose, and of antagonism under constant efficacy, were magnified (Supplementary Fig. S6) in keeping with the slower speed of evolution. Conversely, when drug resistance brought a fitness advantage, evolution was faster and the relative benefits became smaller. Future work can represent resistance mechanisms for each subpopulation more explicitly. Our model also neglects spatial considerations, which can affect the speed of tumor growth, relapse (28), and accumulation of mutations (29). In addition, assuming that drugs are uniformly distributed within the tumor ignores an important spatial effect in which regional variations in perfusion and drug penetration may create heterogeneity in drug concentrations, potentially forming single-drug “sanctuaries.” Single-drug sanctuaries can accelerate the evolution of multidrug resistance by providing escape sites for singly-resistant cells (52). Thus, there may be significant variations in the dynamics of evolution in different tumor regions.

Our study has not explicitly accounted for toxicity, but toxicity is fundamental to our dosing questions. Toxicity determines whether antagonistic drugs can be administered with equal efficacy to other combinations, and whether synergistic drugs can be administered with equal dose to other combinations, thus defining whether situations resembling the constant-dose and constant-efficacy simulations are valid. Note that for combination therapies, the individual doses are not the sole determinant of toxicity, because combinations often have “nonadditive” toxicities that are either more or less severe than the individual toxicities (11, 12). Because toxicity affects tolerable dose limits, understanding combination toxicity is the barrier to determining the preferable treatment strategy.

In prior research, combination toxicity has been addressed through several approaches. Coldman and Murray (53) used Bayes' theorem to calculate the probability of a toxic event for combinations of chemotherapeutics, assuming that toxicity depends only on the total killing of normal cells. Other studies (54, 55) have used toxicity constraints in designing optimal dosing schedules, defining toxicity constraints as the maximum tolerable duration of a treatment pulse delivered at a particular dose, based on clinical data. These studies, as well as (6, 30), found that short pulses of high mono-therapeutic doses could minimize the probability of resistance without crossing the toxicity limit. The studies (6, 30) also found that adaptive switching of treatment yielded better outcomes than prolonged treatment. Future work should study how pulsed or variable dosing strategies would affect combination toxicity and resistance evolution, as a means for addressing toxicity limits. Another unmet need is determining how pharmacokinetic clearance could affect toxicity.

In summary, our work demonstrates the potential for different combination therapies to combat resistance evolution in cancer. Overall, synergism provides a proefficacy approach that can suppress resistance if it is sufficiently efficacious to decimate total cell numbers. Otherwise, synergism becomes a double-edged sword; the reason synergistic combinations are attractive (i.e., a steep increment in efficacy when adding a second drug) creates a dangerously high increment in evolutionary fitness, if some cells develop single-drug resistance. Meanwhile, antagonism provides an antiresistance approach; the reason antagonistic treatments have been avoided (i.e., a poor increase in efficacy when adding a second drug) becomes a protective measure against future resistance, because partially-resistant cells would have little fitness advantage over sensitive cells. The constant-dose experiment suggests that, given the option of therapies where one has significantly higher initial efficacy than others, the more efficacious one would give a better long-term outcome. On the other hand, the constant-efficacy experiment cautions that if two combinations have similar efficacy, then the more antagonistic one would provide better long-term defense against resistance. Our study urges an open-minded consideration of combinations according to their empirical long-term impact, rather than presuming that synergism in the initial efficacy would necessarily produce better long-term outcomes.

No potential conflicts of interest were disclosed.

Conception and design: E.C. Saputra, L. Huang, L. Tucker-Kellogg

Development of methodology: E.C. Saputra, L. Huang, Y. Chen, L. Tucker-Kellogg

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): E.C. Saputra

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): E.C. Saputra, L. Huang, L. Tucker-Kellogg

Writing, review, and/or revision of the manuscript: E.C. Saputra, L. Huang, L. Tucker-Kellogg

Study supervision: L. Tucker-Kellogg

We are grateful to Daniel Tan Shao Weng and Dawn Lau Pingxi for empirical studies of resistance evolution, to David Virshup for helpful comments on the article, to Roy Welsch for discussions about the analytical approximation, to Narendra Suhas Jagannathan and Marie-Véronique Clément for discussions of cancer state transitions, and to Low Boon Chuan and Chen Yu Zong for excellent support of L. Huang. This research was supported by the St. Baldrick's Foundation and the Singapore Ministry of Health's National Medical Research Council (NMRC) under its Open Fund Individual Research Grant scheme (OFIRG15nov062). L. Huang was supported by a doctoral fellowship from the Singapore–MIT Alliance.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.