## Abstract

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.

## Introduction

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.

### 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 ( ) and/or death rate ( );

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 , 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 ( ) of proliferation, death, and rare phenotype alteration events:

**Proliferation and Death:**

**Alteration:**

where

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:**
is the number of cells in subpopulation
at time
,
is the intermediate number of cells in subpopulation
at time
before cells undergo phenotype alterations, ^{prolif}
is the effective proliferation rate, ^{death}
is the effective death rate,
is the random number of cells of subpopulation
that transforms into cell type
at time
,
is the proliferation rate of all subpopulations when untreated,
is the proliferation probability of subpopulation
,
is the basal (natural) apoptosis rate,
is the drug-induced apoptosis rate of subpopulation
, and
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
, we define
as the equivalent dose (17) of drug 1 with the same magnitude of effect as the dose combination (
,
):
is the median effective dose of drug
(the dose of drug
that affects 50% of the population). Resistance is defined as the dose sensitivity (denoted by
, in %) of subpopulation
to the dose of drug
. The interaction parameter
in the Greco model describes synergism if
> 0, additivity if
= 0, and antagonism if
< 0. As different subpopulations sense different fractions of doses
and
, the equivalent dose of drug 1 that is effectively sensed by each subpopulation is denoted by
, where
can be *SS, RS, SR*, or *RR*. From the equivalent dose calculated in Eq. B, we can calculate drug effect from:
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,
, and drug-induced apoptosis rate,
, can be calculated by multiplying the normalized effect by modulating parameters for proliferation,
, and apoptosis,
, as follows:
**Parameters:**
is the dose of drug
,
is the 50% effective dose of drug
,
is the interaction parameter of drug combination in the Greco model,
is the potency ratio constant,
is the equivalent dose of drug 1 sensed by subpopulation
that has the same magnitude of effect as the combined effect of the (
,
) dose pair,
is the normalized drug effect on subpopulation
from treatment,
is the scaling factor for proliferation effect, and
is the scaling factor for apoptotic effect.

#### C. The time until double resistance arises, *T*_{RR}

_{RR}

We derive an analytical approximation for the probability that *RR* cells will arise at time
(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:

**Parameters:**
is the time for doubly-resistant *RR* cells to appear,
is the net growth rate of singly-resistant subpopulations *RS* and *SR* collectively,
is the net growth rate of sensitive *SS* subpopulation,
is the number of *SS* cells at time 0, and
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.

## Materials and Methods

### 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
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
, denoted
, signifies the proliferative potential of cells in the subpopulation under the given therapeutic condition. When exposed to treatment,
decreases from 1 (full proliferative capacity, e.g., untreated cells) to any fraction or 0 (nonproliferative). Meanwhile, the apoptosis potential of subpopulation
was represented by the sum of the cell's basal apoptosis rate,
, and the drug-induced apoptosis rate,
. 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
and basal apoptosis rate
without treatment. During treatment, the proliferation rate would be scaled down by each subpopulation's proliferation probability
, conferring differential fitness across the subpopulations, summarized in the effective proliferation rate ^{prolif}
and effective apoptosis rate ^{death}
(Eq. A). Our main simulations started with 10^{9}*SS* cells (representing a minimum size for detection; ref. 45) and terminated when tumor size reached zero (representing eradication) or 10^{12} 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 , and drug-induced apoptosis rate .

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
. Under monotherapy with drug 1, the proliferation probabilities of cells resistant to drug 1 (
and
) remain approximately unchanged, whereas
and
decline to *x <* 1. Conversely, under drug 2,
and
are approximately unchanged, whereas
and
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 to denote the strength of interaction between a combination, where a positive, zero, or negative signifies synergism, additivity, or antagonism, respectively. For flexibility in defining whether anticancer drugs target proliferation or survival, our model introduced two modulating parameters: for describing the maximum effect on proliferation, and for describing the maximum effect on apoptosis (Eqs. C and D). A larger or indicates a greater effect on proliferation or apoptosis, respectively. For our main simulations, we simulated both effects at once, assuming and . (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., 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 (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
(time to reach 10^{12} cells) and
(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}).

### Approximation of (the time until double resistance arises)

To generalize how synergism and antagonism affect evolution, we derived an analytical approximation of (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 is

Here,
is the net growth rate of singly-resistant subpopulations *RS* and *SR* collectively,
is the net growth rate of subpopulation *SS*,
is the alteration rate, and
is the number of sensitive *SS* cells at time 0. In interpreting the approximation, we focus on
> 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
or
approach 1 (see the denominator of the exponent).

## Results

### 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 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
of 0.69 (the same as in the additive control). Interestingly, these results showed that as treatment became more antagonistic,
and
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
.

### 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
(see Materials and Methods). Figure 3A–E validates the applicability of our analytical equation, showing that the analytical approximation of
agrees closely with the results from Monte Carlo simulations.

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

### 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* (
) 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).

## Discussion

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,
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
and
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
, the definition of success in the Constant-Dose and Constant-Efficacy Methods becomes straightforward: by substituting the appropriate
and
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.

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

## Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

## Authors' Contributions

**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

## Acknowledgments

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.

## Footnotes

**Note:**Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).E.C. Saputra and L. Huang are co-first authors of this article.

- Received May 12, 2017.
- Revision received September 29, 2017.
- Accepted February 12, 2018.
- Published first April 23, 2018.

- ©2018 American Association for Cancer Research.