The search for effective combination therapies for cancer has focused heavily on synergistic combinations because they exhibit enhanced therapeutic efficacy at lower doses. Although synergism is intuitively attractive, therapeutic success often depends on whether drug resistance develops. The impact of synergistic combinations (vs. antagonistic or additive combinations) on the process of drug-resistance evolution has not been investigated. In this study, we use a simplified computational model of cancer cell numbers in a population of drug-sensitive, singly-resistant, and fully-resistant cells to simulate the dynamics of resistance evolution in the presence of two-drug combinations. When we compared combination therapies administered at the same combination of effective doses, simulations showed synergistic combinations most effective at delaying onset of resistance. Paradoxically, when the therapies were compared using dose combinations with equal initial efficacy, antagonistic combinations were most successful at suppressing expansion of resistant subclones. These findings suggest that, although synergistic combinations could suppress resistance through early decimation of cell numbers (making them “proefficacy” strategies), they are inherently fragile toward the development of single resistance. In contrast, antagonistic combinations suppressed the clonal expansion of singly-resistant cells, making them “antiresistance” strategies. The distinction between synergism and antagonism was intrinsically connected to the distinction between offensive and defensive strategies, where offensive strategies inflicted early casualties and defensive strategies established protection against anticipated future threats. Our findings question the exclusive focus on synergistic combinations and motivate further consideration of nonsynergistic combinations for cancer therapy.

Significance: Computational simulations show that if different combination therapies have similar initial efficacy in cancers, then nonsynergistic drug combinations are more likely than synergistic drug combinations to provide a long-term defense against the evolution of therapeutic resistance. Cancer Res; 78(9); 2419–31. ©2018 AACR.

Major Findings

We characterize the theoretical circumstances under which additive, synergistic, or antagonistic combination therapies would suppress evolution of resistance. Antagonistic therapies minimize the competitive advantage of cells that develop single-drug resistance, and thus offer superior performance in delaying resistance, for cases where the dosing limits allow sufficient efficacy. Conversely, synergistic therapies can delay the evolution of resistance if they are sufficiently effective to decimate naïve cancer cells faster than resistant cells can arise. Finally, we use the inherent symmetry between therapeutic efficacy and evolutionary fitness to explain why synergism is a “proefficacy” strategy, whereas antagonism is an “antiresistance” strategy.

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.

Figure 1.

Modeling framework. A, Isobologram showing lines of equal efficacy (isoboles) for additive, synergistic, and antagonistic two-drug combinations, based on Loewe's additivity model. Isobole is defined as the curve of dose combination points having equal efficacy. The x- and y-intercepts represent the single dose of drug 1 only, or drug 2 only, that will generate equal efficacy, and all points on the isobole connecting the two intercepts also have the same efficacy. Drugs are noninteracting if their combined effect lies on the linear additive isobole (solid line). Greco's interaction parameter α = 0 signifies additivity, α < 0 signifies antagonism, and α > 0 signifies synergism. B, Flowchart of the evolutionary tumor cell population model, including “Proliferation and Death phase” to change the cell numbers in each subgroup, “Alteration phase” to allow cells to change subgroup membership, as well as a test for termination. The “Alteration phase” box shows the state transition model for phenotype alterations, where the four states represent four subpopulations: fully-sensitive (SS), resistant to drug 1 (RS), resistant to drug 2 (SR), and doubly-resistant (RR). State transitions are bidirectional and independent, with uniform probability μ = 10−6. Simulation begins with 109SS cells and terminates when cell numbers reach 1012 or zero.

Figure 1.

Modeling framework. A, Isobologram showing lines of equal efficacy (isoboles) for additive, synergistic, and antagonistic two-drug combinations, based on Loewe's additivity model. Isobole is defined as the curve of dose combination points having equal efficacy. The x- and y-intercepts represent the single dose of drug 1 only, or drug 2 only, that will generate equal efficacy, and all points on the isobole connecting the two intercepts also have the same efficacy. Drugs are noninteracting if their combined effect lies on the linear additive isobole (solid line). Greco's interaction parameter α = 0 signifies additivity, α < 0 signifies antagonism, and α > 0 signifies synergism. B, Flowchart of the evolutionary tumor cell population model, including “Proliferation and Death phase” to change the cell numbers in each subgroup, “Alteration phase” to allow cells to change subgroup membership, as well as a test for termination. The “Alteration phase” box shows the state transition model for phenotype alterations, where the four states represent four subpopulations: fully-sensitive (SS), resistant to drug 1 (RS), resistant to drug 2 (SR), and doubly-resistant (RR). State transitions are bidirectional and independent, with uniform probability μ = 10−6. Simulation begins with 109SS cells and terminates when cell numbers reach 1012 or zero.

Close modal
Quick Guide to Equations and Assumptions
A. Simulation of cancer evolution in the presence of combination therapy

Major assumptions:

  • Tumor cells are treated with synergistic, additive, or antagonistic drug combinations such that the difference in therapy affects the proliferation probability (⁠|P{P_i}$|⁠) and/or death rate (⁠|{a_i}$|⁠);

  • The population heterogeneity is categorized into four subpopulations: fully-sensitive (SS), resistant to drug 1 (RS), resistant to drug 2 (SR), or doubly-resistant (RR); we neglect heterogeneity within each subpopulation;

  • Mechanisms of drug resistance are independent for each drug; a fully-sensitive cell cannot become resistant to both drugs without undergoing at least two discrete alteration events;

  • The number of tumor cells in each subpopulation grows or shrinks over time according to an exponential model.

The model tracks |{N_{i,t}}$|⁠, the number of cells over time in each of the sensitive or resistant subpopulations, by simulating evolutionary growth dynamics. The evolutionary dynamics are described by the following stochastic Poisson process (⁠|{\rm{Po}}$|⁠) of proliferation, death, and rare phenotype alteration events:

Proliferation and Death:|$N_{i,t}^* = {\rm{Po}}( {{N_{i,t}}( {1{ +\, ^{{\rm{prolif}}}}{k_i}} )( {1{ -\, ^{{\rm{death}}}}{k_i}} )} )$|

Alteration:|{N_{i,t +\! 1}} = N_{i,t}^* - \sum\nolimits_{j \ne i} {{M_{i \to j,t}}} + \sum\nolimits_{j \ne i} {{M_{j \to i,t}}}, $|

where

formula

Each term in Eq. A is a Poisson random number with mean specified in the parentheses. Phenotype alterations can occur in both forward and reverse directions, and direct transitions between SS and RR are disallowed.

Parameters:|{N_{i,t}}$| is the number of cells in subpopulation |i$| at time |t$|⁠, |N_{i,t}^*$| is the intermediate number of cells in subpopulation |i$| at time |t$| before cells undergo phenotype alterations, prolif|{k_i}$| is the effective proliferation rate, death|{k_i}$| is the effective death rate, |{M_{i \to j,t}}$| is the random number of cells of subpopulation |i$| that transforms into cell type |j$| at time |t$|⁠, |p$| is the proliferation rate of all subpopulations when untreated, |P{P_i}$| is the proliferation probability of subpopulation |i$|⁠, |{a_b}$| is the basal (natural) apoptosis rate, |{a_i}$| is the drug-induced apoptosis rate of subpopulation |i$|⁠, and |\mu $| is the phenotype alteration rate.

B. The effect of combination therapy

Major assumptions:

  • Drug interaction is defined following Loewe's additivity model (14);

  • Each drug combination has a constant potency ratio (15).

We use Greco and colleagues' (16) response surface model for two-drug combinations, as a basis to quantify the effect of combination therapy. Assuming a constant potency ratio |R$|⁠, we define |{D_1}$| as the equivalent dose (17) of drug 1 with the same magnitude of effect as the dose combination (⁠|{d_1}$|⁠, |{d_2}$|⁠):

formula

|{\rm{E}}{{\rm{D}}_{50,y}}$| is the median effective dose of drug |y$| (the dose of drug |y$| that affects 50% of the population). Resistance is defined as the dose sensitivity (denoted by |{S_{y,i}}$|⁠, in %) of subpopulation |i$| to the dose of drug |y$|⁠. The interaction parameter |\alpha $| in the Greco model describes synergism if |\alpha $| > 0, additivity if |\alpha $| = 0, and antagonism if |\alpha $| < 0. As different subpopulations sense different fractions of doses |{d_1}$| and |{d_2}$|⁠, the equivalent dose of drug 1 that is effectively sensed by each subpopulation is denoted by |{D_{1,i}}$|⁠, where |i$| can be SS, RS, SR, or RR. From the equivalent dose calculated in Eq. B, we can calculate drug effect from:

formula

|{E_i}$| can be implemented as either the normalized reduction in proliferation probability, or the normalized drug-induced apoptosis rate, and can be used to parameterize the strength of drug effect in proliferation and/or apoptosis.

The absolute parameter values for proliferation probability, |P{P_i}$|⁠, and drug-induced apoptosis rate, |{a_i}$|⁠, can be calculated by multiplying the normalized effect by modulating parameters for proliferation, |P{P_c}$|⁠, and apoptosis, |{a_c}$|⁠, as follows:

formula

Parameters:|{d_y}$| is the dose of drug |y$|⁠, |{\rm{E}}{{\rm{D}}_{50,y}}$| is the 50% effective dose of drug |y$|⁠, |\alpha $| is the interaction parameter of drug combination in the Greco model, |R$| is the potency ratio constant, |{D_{1,i}}$| is the equivalent dose of drug 1 sensed by subpopulation |i$| that has the same magnitude of effect as the combined effect of the (⁠|{d_1}$|⁠, |{d_2}$|⁠) dose pair, |{E_i}$| is the normalized drug effect on subpopulation |i$| from treatment, |P{P_c}$| is the scaling factor for proliferation effect, and |{a_c}$| is the scaling factor for apoptotic effect.

C. The time until double resistance arises, TRR

We derive an analytical approximation for the probability that RR cells will arise at time |t$| (see Supplementary Text S4.1, Supplementary Equation S16 for the full equation). Assuming the two drugs are delivered symmetrically at the same effective dose, the probability is given by:

formula

Parameters:|{T_{RR}}$| is the time for doubly-resistant RR cells to appear, |{k_{RS + SR}}$| is the net growth rate of singly-resistant subpopulations RS and SR collectively, |{k_{SS}}$| is the net growth rate of sensitive SS subpopulation, |{n_{SS,0}}$| is the number of SS cells at time 0, and |{n_{RS + SR,0}}$| is the number of RS and SR cells collectively at time 0.

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.

Evolutionary tumor cell population model

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.

Evolutionary changes.

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−6 (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).

Cell fitness.

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.

Simulation model.

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 109SS cells (representing a minimum size for detection; ref. 45) and terminated when tumor size reached zero (representing eradication) or 1012 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.

Drug-effect parameters

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.

Table 1.

Table of fitness parameters, illustrated here by proliferation probability |P{P_i}$|⁠, where |i$| references the four subpopulations: fully-sensitive (SS), resistant to drug 1 (RS), resistant to drug 2 (SR), and doubly-resistant (RR)

|P{P_i}$|UntreatedDrug 1Drug 2Drugs 1 and 2Synergistic exampleAntagonistic example
|P{P_{SS}}$| x y ?? 0.35 0.56 
|P{P_{RS}}$| y y y = 0.70 y = 0.70 
|P{P_{SR}}$| x x x = 0.73 x = 0.73 
|P{P_{RR}}$| 1.00 1.00 
|P{P_i}$|UntreatedDrug 1Drug 2Drugs 1 and 2Synergistic exampleAntagonistic example
|P{P_{SS}}$| x y ?? 0.35 0.56 
|P{P_{RS}}$| y y y = 0.70 y = 0.70 
|P{P_{SR}}$| x x x = 0.73 x = 0.73 
|P{P_{RR}}$| 1.00 1.00 

NOTE: x and y are arbitrary variables between 0 and 1. Indicated values illustrate conceptual trends and not exact fitness.

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.

Constant-Dose Method and Constant-Efficacy Method

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 (EDX) 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 ED50,X = 1 mg/kg, while drug Y is delivered at ED50,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 EDX 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 1012 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 EDX).

Approximation of |{{\bm{T}}_{{\bm{RR}}}}$| (the time until double resistance arises)

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

formula

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).

Dynamics of tumor evolution under the Constant-Dose Method and the Constant-Efficacy Method

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 EDX. 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.

Figure 2.

The dynamics of resistance evolution under the Constant-Dose Method and the Constant-Efficacy Method. Simulated growth dynamics of the subpopulations under two-drug combination therapies, from Monte Carlo simulations (n = 10,000). Indigo, SS (sensitive to both drugs); magenta, RS (resistant to drug 1); red, SR (resistant to drug 2, eclipsing the RS curves); cyan, RR (resistant to both drugs); and yellow, the total population. Red numbers in the “Time” axis indicate the time until double resistance emerges, |{T_{RR}}$|⁠. In the simulation parameters at right, an asterisk indicates which parameter is used for defining constant dose or constant efficacy. All simulations can be compared against the additive case, as a control simulation (middle row on the right). Each pair of drugs is delivered at symmetric effective doses (⁠|{\rm{E}}{{\rm{D}}_{\rm{X}}}$|⁠), meaning |{\rm{E}}{{\rm{D}}_{{\rm{X}},1}} = {\rm{E}}{{\rm{D}}_{{\rm{X}},2}}$|⁠. A, Cell numbers over time, simulated according to the parameters at right. Antagonistic and synergistic combination therapies are compared according to the Constant-Dose Method, where all drug pairs have the same combinations of effective doses. B, Repeated simulations except using the Constant-Efficacy Method, where all drug pairs have the same combined effect on sensitive cells, |{E_{SS}} = 0.69$|⁠, while allowing effective doses to vary (see Supplementary Fig. S2C; parameters: |{\rm{E}}{{\rm{D}}_{50,1}}$| = 10, |{\rm{E}}{{\rm{D}}_{50,2}}$| = 10, |p$| = 0.95, |{a_b}$| = 0.12, |P{P_c}$| = 1, |{a_c}$| = 0.254.)

Figure 2.

The dynamics of resistance evolution under the Constant-Dose Method and the Constant-Efficacy Method. Simulated growth dynamics of the subpopulations under two-drug combination therapies, from Monte Carlo simulations (n = 10,000). Indigo, SS (sensitive to both drugs); magenta, RS (resistant to drug 1); red, SR (resistant to drug 2, eclipsing the RS curves); cyan, RR (resistant to both drugs); and yellow, the total population. Red numbers in the “Time” axis indicate the time until double resistance emerges, |{T_{RR}}$|⁠. In the simulation parameters at right, an asterisk indicates which parameter is used for defining constant dose or constant efficacy. All simulations can be compared against the additive case, as a control simulation (middle row on the right). Each pair of drugs is delivered at symmetric effective doses (⁠|{\rm{E}}{{\rm{D}}_{\rm{X}}}$|⁠), meaning |{\rm{E}}{{\rm{D}}_{{\rm{X}},1}} = {\rm{E}}{{\rm{D}}_{{\rm{X}},2}}$|⁠. A, Cell numbers over time, simulated according to the parameters at right. Antagonistic and synergistic combination therapies are compared according to the Constant-Dose Method, where all drug pairs have the same combinations of effective doses. B, Repeated simulations except using the Constant-Efficacy Method, where all drug pairs have the same combined effect on sensitive cells, |{E_{SS}} = 0.69$|⁠, while allowing effective doses to vary (see Supplementary Fig. S2C; parameters: |{\rm{E}}{{\rm{D}}_{50,1}}$| = 10, |{\rm{E}}{{\rm{D}}_{50,2}}$| = 10, |p$| = 0.95, |{a_b}$| = 0.12, |P{P_c}$| = 1, |{a_c}$| = 0.254.)

Close modal

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 EDX 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}}}}$|⁠.

Stochastic simulations and analytical approximations show the merits of synergism under the Constant-Dose Method, and antagonism under the Constant-Efficacy Method

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.

Figure 3.

Histograms showing |{T_{RR}},$| the time until double resistance emerges. For each frame, the histogram of 10,000 stochastic simulations (gray bars) is superimposed with the probability distribution curve of |{T_{RR}}$| generated from the analytical approximation. The parameters are specified in the tables in Fig. 2 for the additive control (A), the antagonistic example (B), and the synergistic example (C) of the Constant-Dose Method in Fig. 2A, and the antagonistic example (D) and the synergistic example (E) of the Constant-Efficacy Method in Fig. 2B.

Figure 3.

Histograms showing |{T_{RR}},$| the time until double resistance emerges. For each frame, the histogram of 10,000 stochastic simulations (gray bars) is superimposed with the probability distribution curve of |{T_{RR}}$| generated from the analytical approximation. The parameters are specified in the tables in Fig. 2 for the additive control (A), the antagonistic example (B), and the synergistic example (C) of the Constant-Dose Method in Fig. 2A, and the antagonistic example (D) and the synergistic example (E) of the Constant-Efficacy Method in Fig. 2B.

Close modal

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).

Figure 4.

|{T_{{\rm{lethal}}}}$| and |{T_{RR}}$| from treatment with different combinations, according to the Constant-Dose Method and the Constant-Efficacy Method. A and C were obtained by running Monte Carlo simulations with 10,000 experiments for each treatment scenario. A, Under the Constant-Dose Method, increasing α from antagonistic (⁠|\alpha < 0$|⁠) to synergistic (⁠|\alpha > 0$|⁠) prolongs the time to reach a lethal tumor burden, |{T_{{\rm{lethal}}}}$| (solid line), and the time for RR cells to emerge, |{T_{RR}}$| (dashed line). The gray shaded areas indicate the respective SDs of |{T_{{\rm{lethal}}}}$| and |{T_{RR}}$|⁠. B, The analytical approximation of |{T_{RR}}$| under the Constant-Dose Method agrees with the results in A, because the probability distribution of |{T_{RR}}$| from the approximation shifts to higher values as treatment becomes more synergistic (i.e., from blue to red), indicating more time before doubly-resistant cells emerge. C, In contrast, under the Constant-Efficacy Method, decreasing |\alpha $| from synergistic (⁠|\alpha \,> \,0$|⁠) to antagonistic (⁠|\alpha < 0$|⁠) prolongs |{T_{{\rm{lethal}}}}$| and |{T_{RR}}$|⁠. D, The analytical approximation of |{T_{RR}}$| under the Constant-Efficacy Method shows that the probability distribution of |{T_{RR}}$| shifts to higher values and becomes more widely spread as treatment becomes more antagonistic (i.e., from red to blue; parameters: the same as in Fig. 2, except for the |\alpha $| values.)

Figure 4.

|{T_{{\rm{lethal}}}}$| and |{T_{RR}}$| from treatment with different combinations, according to the Constant-Dose Method and the Constant-Efficacy Method. A and C were obtained by running Monte Carlo simulations with 10,000 experiments for each treatment scenario. A, Under the Constant-Dose Method, increasing α from antagonistic (⁠|\alpha < 0$|⁠) to synergistic (⁠|\alpha > 0$|⁠) prolongs the time to reach a lethal tumor burden, |{T_{{\rm{lethal}}}}$| (solid line), and the time for RR cells to emerge, |{T_{RR}}$| (dashed line). The gray shaded areas indicate the respective SDs of |{T_{{\rm{lethal}}}}$| and |{T_{RR}}$|⁠. B, The analytical approximation of |{T_{RR}}$| under the Constant-Dose Method agrees with the results in A, because the probability distribution of |{T_{RR}}$| from the approximation shifts to higher values as treatment becomes more synergistic (i.e., from blue to red), indicating more time before doubly-resistant cells emerge. C, In contrast, under the Constant-Efficacy Method, decreasing |\alpha $| from synergistic (⁠|\alpha \,> \,0$|⁠) to antagonistic (⁠|\alpha < 0$|⁠) prolongs |{T_{{\rm{lethal}}}}$| and |{T_{RR}}$|⁠. D, The analytical approximation of |{T_{RR}}$| under the Constant-Efficacy Method shows that the probability distribution of |{T_{RR}}$| shifts to higher values and becomes more widely spread as treatment becomes more antagonistic (i.e., from red to blue; parameters: the same as in Fig. 2, except for the |\alpha $| values.)

Close modal

Pre-existing resistance

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).

Figure 5.

The relative advantages of synergistic or antagonistic therapies in simulations with pre-existing single-drug resistance. Stochastic simulations (10,000 experiments for each scenario) of cancer populations with increasing numbers of pre-existing singly-resistant (RS and SR) cells at the start of treatment. We define |{( {{T_{RR}}} )_{{\rm{strategy}}}}\ $|to be the number of generations until double resistance arises under a given strategy (synergistic, additive, or antagonistic), and ΔTRR to be the difference |{( {{T_{RR}}} )_{{\rm{syn}}}} - {( {{T_{RR}}} )_{{\rm{ant}}}}$|⁠. Using the Constant-Dose Method, synergism consistently results in longer |{T_{RR}}$| than antagonism (i.e., ΔTRR > 0), and when using the Constant-Efficacy Method, antagonism consistently results in longer |{T_{RR}}$| than synergism (i.e., ΔΔTRR < 0). However, the differences between synergism and antagonism diminish as the number of pre-existing singly-resistant cells increases. In this figure, RS and SR cells are present in equal numbers at the beginning of the simulation. The synergistic simulation (⁠|\alpha $| = 0.4) is compared with |\alpha $| = −0.4 antagonism for the Constant-Dose Method and |\alpha $| = −0.3 for the Constant-Efficacy Method. The change in |\alpha $| is to avoid simulations where the cancer is eradicated, without changing the other key parameters. All other parameters are the same as in Fig. 2.

Figure 5.

The relative advantages of synergistic or antagonistic therapies in simulations with pre-existing single-drug resistance. Stochastic simulations (10,000 experiments for each scenario) of cancer populations with increasing numbers of pre-existing singly-resistant (RS and SR) cells at the start of treatment. We define |{( {{T_{RR}}} )_{{\rm{strategy}}}}\ $|to be the number of generations until double resistance arises under a given strategy (synergistic, additive, or antagonistic), and ΔTRR to be the difference |{( {{T_{RR}}} )_{{\rm{syn}}}} - {( {{T_{RR}}} )_{{\rm{ant}}}}$|⁠. Using the Constant-Dose Method, synergism consistently results in longer |{T_{RR}}$| than antagonism (i.e., ΔTRR > 0), and when using the Constant-Efficacy Method, antagonism consistently results in longer |{T_{RR}}$| than synergism (i.e., ΔΔTRR < 0). However, the differences between synergism and antagonism diminish as the number of pre-existing singly-resistant cells increases. In this figure, RS and SR cells are present in equal numbers at the beginning of the simulation. The synergistic simulation (⁠|\alpha $| = 0.4) is compared with |\alpha $| = −0.4 antagonism for the Constant-Dose Method and |\alpha $| = −0.3 for the Constant-Efficacy Method. The change in |\alpha $| is to avoid simulations where the cancer is eradicated, without changing the other key parameters. All other parameters are the same as in Fig. 2.

Close modal

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.

Constant-dose and constant-efficacy comparisons

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.

Comparative fitness rewards from synergism and antagonism

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.

Figure 6.

Synergistic therapy as a “proefficacy” strategy and antagonistic therapy as an “antiresistance” strategy. All plots were computed using the parameters in Supplementary Table S1, with fitness defined as |1 - {E_i}$| (efficacy) in equations B and C, and resistance defined by symmetrically increasing |{S_{1,i}}$| and |{S_{2,i}}$| from 0% to 100%. Drug effect plots are bar charts showing the efficacy of treatments delivered as single agents or in combination, and evolution plots show the fitness of the cancer cells over a spectrum of drug-resistance levels. The drug effect plot (A) and the evolution plot (B) for the Constant-Dose Method, and the drug effect plot (C) and the evolution plot (D) for the Constant-Efficacy Method. By this analysis, synergism is a “proefficacy” strategy because it has greater-than-additive efficacy, but it suffers from large Δfitness/Δresistance in the evolution plot, meaning that synergism rewards cells that start to develop resistance. Conversely, antagonism is an “anti-resistance” strategy because it has smaller Δfitness/Δresistance, but it suffers less-than-additive efficacy.

Figure 6.

Synergistic therapy as a “proefficacy” strategy and antagonistic therapy as an “antiresistance” strategy. All plots were computed using the parameters in Supplementary Table S1, with fitness defined as |1 - {E_i}$| (efficacy) in equations B and C, and resistance defined by symmetrically increasing |{S_{1,i}}$| and |{S_{2,i}}$| from 0% to 100%. Drug effect plots are bar charts showing the efficacy of treatments delivered as single agents or in combination, and evolution plots show the fitness of the cancer cells over a spectrum of drug-resistance levels. The drug effect plot (A) and the evolution plot (B) for the Constant-Dose Method, and the drug effect plot (C) and the evolution plot (D) for the Constant-Efficacy Method. By this analysis, synergism is a “proefficacy” strategy because it has greater-than-additive efficacy, but it suffers from large Δfitness/Δresistance in the evolution plot, meaning that synergism rewards cells that start to develop resistance. Conversely, antagonism is an “anti-resistance” strategy because it has smaller Δfitness/Δresistance, but it suffers less-than-additive efficacy.

Close modal

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.

Synergism as an offensive strategy and antagonism as a defensive strategy

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.

Caveats

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−6, 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.

Toxicity and dosing

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.

1.
Flaherty
KT
,
Infante
JR
,
Daud
A
,
Gonzalez
R
,
Kefford
RF
,
Sosman
J
, et al
Combined BRAF and MEK inhibition in melanoma with BRAF V600 mutations
.
N Engl J Med
2012
;
367
:
1694
703
.
2.
Sharma
P
,
Allison
JP
. 
Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential
.
Cell
2015
;
161
:
205
14
.
3.
Bae
SY
,
Hong
J-Y
,
Lee
H-J
,
Park
HJ
,
Lee
SK
. 
Targeting the degradation of AXL receptor tyrosine kinase to overcome resistance in gefitinib-resistant non-small cell lung cancer
.
Oncotarget
2015
;
6
:
10146
60
.
4.
Choi
YJ
,
Kim
SY
,
So
KS
,
Baek
I-J
,
Kim
WS
,
Choi
SH
, et al
AUY922 effectively overcomes MET- and AXL-mediated resistance to EGFR-TKI in lung cancer cells
.
PLoS One
2015
;
10
:
e0119832
.
5.
Komarova
N
. 
Stochastic modeling of drug resistance in cancer
.
J Theor Biol
2006
;
239
:
351
66
.
6.
Beckman
RA
,
Schemmann
GS
,
Yeang
C-H
. 
Impact of genetic dynamics and single-cell heterogeneity on development of nonstandard personalized medicine strategies for cancer
.
Proc Natl Acad Sci U S A
2012
;
109
:
14586
91
.
7.
Huang
L
.
Simulation and network analysis of MAPK pathway regulation
.
Singapore
:
Singapore-MIT Alliance National University of Singapore
; 
2012
.
8.
Bozic
I
,
Reiter
JG
,
Allen
B
,
Antal
T
,
Chatterjee
K
,
Shah
P
, et al
Evolutionary dynamics of cancer in response to targeted combination therapy
.
eLife
2013
;
2
:
e00747
.
9.
Coldman
AJ
,
Goldie
JH
. 
A model for the resistance of tumor cells to cancer chemotherapeutic agents
.
Math Biosci
1983
;
65
:
291
307
.
10.
Chou
T-C
. 
Theoretical basis, experimental design, and computerized simulation of synergism and antagonism in drug combination studies
.
Pharmacol Rev
2006
;
58
:
621
81
.
11.
Cavalli
F
,
Hansen
HH
,
Kaye
SB
,
European Society for Medical Oncology
.
Textbook of medical oncology
.
London; New York
:
Taylor & Francis
; 
2004
.
12.
Nussbaum
JC
,
Jackson
A
,
Namarika
D
,
Phulusa
J
,
Kenala
J
,
Kanyemba
C
, et al
Combination flucytosine and highdose fluconazole compared with fluconazole monotherapy for the treatment of cryptococcal meningitis: a randomized trial in Malawi
.
Clin Infect Dis
2010
;
50
:
338
44
.
13.
Lehár
J
,
Krueger
AS
,
Avery
W
,
Heilbut
AM
,
Johansen
LM
,
Price
ER
, et al
Synergistic drug combinations tend to improve therapeutically relevant selectivity
.
Nat Biotechnol
2009
;
27
:
659
66
.
14.
Loewe
S
. 
Die quantitativen Probleme der Pharmakologie
.
Ergeb Physiol Biol Chem Exp Pharmakol
1928
;
27
:
47
187
.
15.
Loewe
S
. 
The problem of synergism and antagonism of combined drugs
.
Arzneimittelforschung
1953
;
3
:
285
90
.
16.
Greco
WR
,
Bravo
G
,
Parsons
JC
. 
The search for synergy: a critical review from a response surface perspective
.
Pharmacol Rev
1995
;
47
:
331
85
.
17.
Tallarida
RJ
. 
An overview of drug combination analysis with isobolograms
.
J Pharmacol Exp Ther
2006
;
319
:
1
7
.
18.
Chait
R
,
Craney
A
,
Kishony
R
. 
Antibiotic interactions that select against resistance
.
Nature
2007
;
446
:
668
71
.
19.
Hegreness
M
,
Shoresh
N
,
Damian
D
,
Hartl
D
,
Kishony
R
. 
Accelerated evolution of resistance in multidrug environments
.
Proc Natl Acad Sci U S A
2008
;
105
:
13977
81
.
20.
Michel
J-B
,
Yeh
PJ
,
Chait
R
,
Moellering
RC
,
Kishony
R
. 
Drug interactions modulate the potential for evolution of resistance
.
Proc Natl Acad Sci U S A
2008
;
105
:
14918
23
.
21.
Yeh
PJ
,
Hegreness
MJ
,
Aiden
AP
,
Kishony
R
. 
Drug interactions and the evolution of antibiotic resistance
.
Nat Rev Microbiol
2009
;
7
:
460
6
.
22.
Foo
J
,
Michor
F
. 
Evolution of acquired resistance to anti-cancer therapy
.
J Theor Biol
2014
;
355
:
10
20
.
23.
Nowell
PC
. 
The clonal evolution of tumor cell populations
.
Science
1976
;
194
:
23
8
.
24.
Gatenby
RA
,
Brown
J
,
Vincent
T
. 
Lessons from applied ecology: cancer control using an evolutionary double bind
.
Cancer Res
2009
;
69
:
7499
502
.
25.
Sun
D
,
Dalin
S
,
Hemann
MT
,
Lauffenburger
DA
,
Zhao
B
. 
Differential selective pressure alters rate of drug resistance acquisition in heterogeneous tumor populations
.
Sci Rep
2016
;
6
:
36198
.
26.
Brown
JS
,
Cunningham
JJ
,
Gatenby
RA
. 
Aggregation effects and population-based dynamics as a source of therapy resistance in cancer
.
IEEE Trans Biomed Eng
2017
;
64
:
512
8
.
27.
Greene
J
,
Lavi
O
,
Gottesman
MM
,
Levy
D
. 
The impact of cell density and mutations in a model of multidrug resistance in solid tumors
.
Bull Math Biol
2014
;
76
:
627
53
.
28.
Waclaw
B
,
Bozic
I
,
Pittman
ME
,
Hruban
RH
,
Vogelstein
B
,
Nowak
MA
. 
A spatial model predicts that dispersal and cell turnover limit intratumour heterogeneity
.
Nature
2015
;
525
:
261
4
.
29.
Komarova
NL
. 
Spatial interactions and cooperation can change the speed of evolution of complex phenotypes
.
Proc Natl Acad Sci U S A
2014
;
111
(
Suppl_3
):
10789
95
.
30.
Yeang
C-H
,
Beckman
RA
. 
Long range personalized cancer treatment strategies incorporating evolutionary dynamics
.
Biol Direct
2016
;
11
:
56
.
31.
Cunningham
JJ
,
Gatenby
RA
,
Brown
JS
. 
Evolutionary dynamics in cancer therapy
.
Mol Pharm
2011
;
8
:
2094
100
.
32.
Zhao
B
,
Sedlak
JC
,
Srinivas
R
,
Creixell
P
,
Pritchard
JR
,
Tidor
B
, et al
Exploiting temporal collateral sensitivity in tumor clonal evolution
.
Cell
2016
;
165
:
234
46
.
33.
Yu
HA
,
Sima
C
,
Feldman
D
,
Liu
LL
,
Vaitheesvaran
B
,
Cross
J
, et al
Phase 1 study of twice weekly pulse dose and daily low-dose erlotinib as initial treatment for patients with EGFR -mutant lung cancers
.
Ann Oncol
2017
;
28
:
278
84
.
34.
Jonsson
VD
,
Blakely
CM
,
Lin
L
,
Asthana
S
,
Matni
N
,
Olivas
V
, et al
Novel computational method for predicting polytherapy switching strategies to overcome tumor heterogeneity and evolution
.
Sci Rep
2017
;
7
:
44206
.
35.
Chmielecki
J
,
Foo
J
,
Oxnard
GR
,
Hutchinson
K
,
Ohashi
K
,
Somwar
R
, et al
Optimization of dosing for EGFR-mutant non-small cell lung cancer with evolutionary cancer modeling
.
Sci Transl Med
2011
;
3
:
90ra59
.
36.
Pritchard
JR
,
Lauffenburger
DA
,
Hemann
MT
. 
Understanding resistance to combination chemotherapy
.
Drug Resist Updat
2012
;
15
:
249
57
.
37.
Komarova
NL
,
Wodarz
D
. 
Combination therapies against chronic myeloid leukemia: short-term versus long-term strategies
.
Cancer Res
2009
;
69
:
4904
10
.
38.
Maley
CC
,
Reid
BJ
,
Forrest
S
. 
Cancer prevention strategies that address the evolutionary dynamics of neoplastic cells: simulating benign cell boosters and selection for chemosensitivity
.
Cancer Epidemiol Biomark Prev
2004
;
13
:
1375
84
.
39.
Rodriguez-Brenes
IA
,
Komarova
NL
,
Wodarz
D
. 
Tumor growth dynamics: insights into evolutionary processes
.
Trends Ecol Evol
2013
;
28
:
597
604
.
40.
Obe
G
.
Mutations in man [Internet]
.
Berlin, Heidelberg
:
Springer Berlin Heidelberg
; 
1984
[cited 2017 Jul 17]. Available from
: http://public.eblib.com/choice/publicfullrecord.aspx?p=3091185
41.
Duesberg
P
,
Stindl
R
,
Hehlmann
R
. 
Explaining the high mutation rates of cancer cells to drug and multidrug resistance by chromosome reassortments that are catalyzed by aneuploidy
.
Proc Natl Acad Sci U S A
2000
;
97
:
14295
300
.
42.
Nim
TH
,
Luo
L
,
Clément
M-V
,
White
JK
,
Tucker-Kellogg
L
. 
Systematic parameter estimation in data-rich environments for cell signalling dynamics
.
Bioinforma Oxf Engl
2013
;
29
:
1044
51
.
43.
Shi
Y
,
Mellier
G
,
Huang
S
,
White
J
,
Pervaiz
S
,
Tucker-Kellogg
L
. 
Computational modelling of LY303511 and TRAIL-induced apoptosis suggests dynamic regulation of cFLIP
.
Bioinforma Oxf Engl
2013
;
29
:
347
54
.
44.
Singh
A
,
Settleman
J
. 
EMT, cancer stem cells and drug resistance: an emerging axis of evil in the war on cancer
.
Oncogene
2010
;
29
:
4741
51
.
45.
Friberg
S
,
Mattson
S
. 
On the growth rates of human malignant tumors: implications for medical decision making
.
J Surg Oncol
1997
;
65
:
284
97
.
46.
Baskar
R
,
Dai
J
,
Wenlong
N
,
Yeo
R
,
Yeoh
K-W
. 
Biological response of cancer cells to radiation treatment
.
Front Mol Biosci
2014
;
1
:
24
.
47.
Archer
CD
,
Parton
M
,
Smith
IE
,
Ellis
PA
,
Salter
J
,
Ashley
S
, et al
Early changes in apoptosis and proliferation following primary chemotherapy for breast cancer
.
Br J Cancer
2003
;
89
:
1035
41
.
48.
Illmer
T
,
Schaich
M
,
Platzbecker
U
,
Freiberg-Richter
J
,
Oelschlägel
U
,
von Bonin
M
, et al
P-glycoprotein-mediated drug efflux is a resistance mechanism of chronic myelogenous leukemia cells to treatment with imatinib mesylate
.
Leukemia
2004
;
18
:
401
8
.
49.
Zardavas
D
,
Irrthum
A
,
Swanton
C
,
Piccart
M
. 
Clinical management of breast cancer heterogeneity
.
Nat Rev Clin Oncol
2015
;
12
:
381
94
.
50.
Kahen
E
,
Yu
D
,
Harrison
DJ
,
Clark
J
,
Hingorani
P
,
Cubitt
CL
, et al
Identification of clinically achievable combination therapies in childhood rhabdomyosarcoma
.
Cancer Chemother Pharmacol
2016
;
78
:
313
23
.
51.
Meacham
CE
,
Morrison
SJ
. 
Tumour heterogeneity and cancer cell plasticity
.
Nature
2013
;
501
:
328
37
.
52.
Fu
F
,
Nowak
MA
,
Bonhoeffer
S
. 
Spatial heterogeneity in drug concentrations can facilitate the emergence of resistance to cancer therapy
.
PLOS Comput Biol
2015
;
11
:
e1004142
.
53.
Coldman
AJ
,
Murray
JM
. 
Optimal control for a stochastic model of cancer chemotherapy
.
Math Biosci
2000
;
168
:
187
200
.
54.
Foo
J
,
Michor
F
. 
Evolution of resistance to targeted anti-cancer therapies during continuous and pulsed administration strategies
.
PLoS Comput Biol
2009
;
5
:
e1000557
.
55.
Foo
J
,
Choi
N
,
Leder
K
,
Mumenthaler
S
,
Pao
W
,
Michor
F
, et al
The impact of microenvironmental heterogeneity on the evolution of drug resistance in cancer cells
.
Cancer Inform
2015
;
19
.

Supplementary data