## Abstract

Interactions among drugs play a critical role in the killing efficacy of multi-drug treatments. Recent advances in theory and experiment for three-drug interactions enable the search for emergent interactions—ones not predictable from pairwise interactions. Previous work has shown it is easier to detect synergies and antagonisms among pairwise interactions when a rescaling method is applied to the interaction metric. However, no study has carefully examined whether new types of normalization might be needed for emergence. Here, we propose several rescaling methods for enhancing the classification of the higher order drug interactions based on our conceptual framework. To choose the rescaling that best separates synergism, antagonism and additivity, we conducted bacterial growth experiments in the presence of single, pairwise and triple-drug combinations among 14 antibiotics. We found one of our rescaling methods is far better at distinguishing synergistic and antagonistic emergent interactions than any of the other methods. Using our new method, we find around 50% of emergent interactions are additive, much less than previous reports of greater than 90% additivity. We conclude that higher order emergent interactions are much more common than previously believed, and we argue these findings for drugs suggest that appropriate rescaling is crucial to infer higher order interactions.

## 1. Introduction

Multi-drug treatments are an important tool [1–5], in particular, for combatting bacteria that are highly resistant to the individual use of traditional antibiotics [6–11]. The efficacy and efficiency of these combination therapies are substantially affected by how the specific drugs interact. Thus, a useful categorization scheme for interactions is needed that uses the additive case [6,12,13]—drugs do not interact at all—as a baseline. Along with this is the concept of emergence—effects of drug combinations that cannot be predicted from lower order interactions among subsets of the drugs [14]. Relative to these baselines, interactions are generally categorized as a type of synergy if the combination kills more efficiently than is expected from the additive case or from lower order interactions (figure 1). Conversely, when the drug interaction reduces the effect of each drug, the interaction is called antagonistic (figure 1), which itself contains special cases: (i) buffering, in which one drug completely masks the effect of the other drug and (ii) suppression, in which the effectiveness of the drugs in combination is weaker than at least one drug by itself [8,15–17].

It is a challenge to quantify interactions with a metric that has clear boundaries between these cases and effectively identifies and distinguishes between interaction types. Overcoming this challenge often requires a rescaling or normalization of basic metrics, and for pairwise interactions, some effective methods have been discovered [8,18–20]. Rescaling for combinations of more than two drugs has previously been done using the most straightforward generalization of the pairwise method [14,21], but unlike the pairwise case, there are several possibilities for rescaling metrics for higher order interactions, even among three drugs. Here, we explore these possible rescalings and identify the one that is best at categorizing and thus distinguishing among three-drug interactions. Importantly, this particular rescaling method is different from the rescaling previously used in the literature by ourselves and others. Consequently, this could reveal new insights for interactions among three or more drugs, as well as other objects such as proteins [22] and predators [23–27], and some existing results may need to be revisited and revised [14,28].

For the case of two drugs, there are two common methods to analyse and categorize interactions. First, Loewe additivity categorizes the interaction based on the strength of inhibition on bacterial growth when drug concentrations are varied simultaneously [12]. When the same effective concentration (relative to each drug's minimum inhibitory concentration) of the two drugs kill bacteria at a rate that is independent of the relative fractions of each drug, the interaction is regarded as additive. This Loewe measure is motivated by the simplest case in which a drug does not interact with itself. Depending on the direction of the divergence from this additive case, the relationship is considered as either synergy or antagonism.

Second, Bliss independence (BI) defines additivity to be the case when the presence of one drug does not affect another drug's per cent reduction of bacterial growth rate [6]. This definition breaks down if one of the drugs and the pairwise combination are both lethal such that bacteria cannot grow. Interactions are antagonistic according to BI when the deviation from additivity (DA) is positive and are synergistic when DA is negative [8,18]. BI offers a simple measurement of the epistatic interactions because it relies on less data and its results are more easily calculable and interpretable [17].

The DA measures yield a unimodal distribution around the additive case for drug interactions. Because of the unimodal shape, it is challenging to delineate boundaries and tease apart synergistic, additive and antagonistic cases. To overcome this hurdle, a rescaling method of the DA measure was proposed in Segre *et al*. [18] for pairwise interactions. The rescaling normalizes DA with respect to the pairwise drug fitness for two reference cases: lethal synergy and complete buffering (i.e. the combined two drugs have the same effect as the strongest single-drug alone). This rescaled form of DA leads to an interaction distribution that exhibits three dominant peaks, with clear spacing between them, hereafter referred to as a trimodal distribution. For both simulated and empirical results, these peaks enable a straightforward separation among synergistic, antagonistic and additive interactions [8,18–20]. These peaks are observed at the exact location expected theoretically.

Given that higher order (greater than 2) drug combinations are increasingly being used to combat drug-resistant pathogens, it is important to have similarly effective, though not necessarily similar in form, rescaled measures for higher order interactions within a complex environment. This is a complicated task, because effects at all levels—single drug, pairwise combination, triple combination, etc.—may need to be taken into account (figure 1). For example, a three-drug combination could have interactions arising from the three different pairwise combinations as well as an interaction that only emerges when all three drugs are present. Recently, a novel method to characterize and quantify emergent interactions in three-drug combinations has been introduced. Beppler *et al*. [14] present a framework that compares the higher order interaction with expectations based on its lower order component interactions. First, the direct extension of the DA metric allows identification of three-drug interaction compared solely to the single-drug effects [14,21,29]. Next, the new emergent three-way interaction (E3) incorporates the pairwise interactions in the model and determines the effect of three-drug interaction that is beyond the effects from all two-drug combinations. The generalization of the interaction formulae for combinations of more than three drugs is followed by the conceptual derivation of each metric, and the special case with four drugs is provided in Beppler *et al*. [14]. Moreover, a recent numerical model by Wood *et al*. [28] used maximum entropy estimation to predict the higher order effects relative to single and pairwise drug effects for any number of drug combinations. They showed that their numerical estimation is consistent with an algebraic expression that is equivalent to E3 for three-drug combinations. However, our emergent *N*-way interaction metric differs from Wood *et al*.'s model in the sense that it quantifies the deviations from the expectations from all lower order effects, not just pairwise effects, so for four drugs or more, our model and Wood *et al*.'s model will differ.

In this paper, we show that there are several choices for rescaling emergent interactions and that the specific choice of rescaling plays a crucial role in identifying interactions among drugs. We analyse the previously defined rescaling method—direct extension of the two-drug rescaling—of the emergent three-way measure and establish new rescaling methods that greatly improve the characterization of emergent properties. These methods are defined by exploring possible reference cases of synergy and antagonism that arise with higher order drug combinations. Lethality always serves the reference case for synergies as it offers the most extreme case of synergism regardless of number of drugs in the environment. Therefore, there is only one simple extension of the rescaling method for two-drug combinations for synergies. However, for antagonistically interacting drugs, the definition for the reference case of complete buffering could vary based on the drug system and also the interaction metric (DA versus E3). Because DA quantifies the overall interaction with respect to individual drug effects, the strongest single-drug effect defines the complete buffering case as in the two-drug system. However, for the emergent interaction, it matters whether the buffering is defined relative to some subset or to all of the lower order drug combinations (single drug, two drug or some combination). Our new rescaling methods cover the possibilities for the definition of the antagonistic buffering, and our evaluation of these possibilities eventually leads to enhanced classification of the emergent interaction.

To empirically study whether these different scaling methods effectively separate the histogram for the metric into a trimodal distribution, making it straightforward to separate synergy, additivity and antagonism, we selected 14 antibiotics and systematically investigated the effects of three-drug combinations on the growth rate of a bacterium, wild-type *Escherichia coli*. These data allow us to apply several rescaling methods and hence to identify which method is best at distinguishing among interactions for three-drug combinations and emergent interactions. We further apply our rescaling analysis to the three-antibiotic combination data presented in Wood *et al*. [28]. We compare the original analysis in the Wood *et al*. paper [28] with the new analysis from our rescaling method. Finally, we present a straightforward generalization of our new rescaling methods to emergent interactions that involve more than three drugs.

## 2. Material and methods

### 2.1. Experimental details

#### 2.1.1. Bacteria

The bacteria used in these experiments was the *E. coli* strain BW25113, the wild-type strain (*lacI*^{q} *rrnB _{T14}* Δ

*lacZ*Δ

_{WJ16}hsdR514*araBAD*Δ

_{AH33}*rhaBAD*) [30] derived from the strain W1485 background [31]. A single colony was used to inoculate cultures for glycerol stocks stored at −80°C. A single colony from this glycerol culture was used to inoculate cultures in LB media (10 g l

_{LD78}^{−1}tryptone, 5 g l

^{−1}yeast extract and 10 g l

^{−1}NaCl). These cultures were resuspended in MC buffer and stored at 4°C. Bacteria for experiments were grown by inoculating 20 µl of the MC stock into 2 ml LB daily and growing for 5 h at 37°C; 25 µl of a 10

^{−4}dilution of this culture in LB was used to inoculate into 975 µl cultures for experiments.

#### 2.1.2. Antibiotics

Antibiotics used include clindamycin hydrochloride (Sigma C-5269), chloramphenicol succinate sodium salt (Sigma C3787), fusidic acid sodium salt (Sigma F0881), erythromycin (Sigma Aldrich E-6376), ciprofloxacin hydrochloride (MP Biomedicals 199020), cefoxitin sodium salt (Fluka C4786), ampicillin (Sigma A9518), nitrofurantoin (Sigma N7878), trimethoprim (Sigma T7883), tobramycin sulfate (Sigma T-1783), streptomycin sulfate (Sigma Aldrich S-6501), gentamicin sulfate salt (Sigma G1264), vancomycin hydrochloride (Sigma V2002) and doxycycline hyclate (Sigma D9891). The dosage, mechanism of action and abbreviations of these antibiotics are given in table 1.

#### 2.1.3. Growth experiments

Antibiotic concentrations were chosen to reduce growth by 15–35% as compared to the no-drug control (LB). These sub-lethal dosages were first determined by testing a range of concentrations for each antibiotic and were then used in all single-drug, two-drug and three-drug conditions. Each triple-drug experiment included a no-drug control, three single-drug conditions, three two-drug conditions and the three-drug combination. In all cases, antibiotics were added to the previously determined sub-lethal concentrations in 1 ml LB and inoculated as described above; 100 µl was aliquoted into 4–6 wells for each condition in a 96-well plate. These cultures were grown overnight for 18 h at 215 r.p.m. and 37°C. Optical density readings at 600 nm (OD_{600}) at 18 h were used to calculate growth rates as compared to the no-drug control at 18 h by taking their ratios at this time point. This procedure has been previously used by us [32] and others [33,34]. Each three-drug experiment was repeated at least three times. Data are represented as median, maximum and minimum.

Growth experiment data were also obtained from Wood *et al.* [28], and methods for those empirical measurements are detailed in that paper. Here, the Wood *et al*. [28] triple-drug combination data (with non-zero concentration for each drug in the combination) are analysed using the exact same analytical DA and E3 measures, calculations and cut-offs for significance (all these methods are explained in the subsequent sections) as for our empirical data.

#### 2.1.4. Definition of lethal

Lethal antibiotic concentrations are all concentrations above which no bacterial growth can be measured. Experiments were performed to determine the maximum OD_{600} measurements with error that represent the lethal case. Bacteria for these experiments were grown and inoculated in 96-well plates in the same manner as described above. Three conditions were tested: LB only (no cells), LB + cells and LB + streptomycin (STR) 9.2 µg ml^{−1} + cells. STR 9.2 µg ml^{−1} was chosen as an extremely high antibiotic concentration that would ensure no bacterial growth, even when the populations were re-inoculated in no-drug environments. Thus, the LB only (no cells) and the LB + STR 9.2 µg ml^{−1} + cells conditions could be used to determine the error in OD_{600} measurements that represent no bacterial growth. The LB + cells condition was used as the positive control and 100% growth reference point. The LB only and LB + cells conditions were each replicated in 16 wells. The LB + STR 9.2 µg ml^{−1} + cells condition was replicated in 64 wells. After 18 h of growth at 37°C and 215 r.p.m., the OD_{600} measurements were gathered. From these extremely high drug concentrations as well as populations with no bacterial growth at all, we obtained mean OD_{600} measurements of 0.044 with an error of 0.003. In addition, we further tested bacteria populations at low OD, but above 0.047, to confirm that those populations continued to grow (electronic supplementary material, figure S1). Thus, all growth measurements below 0.047 represent lethal cases.

Based on the analysis of the E3 measure, we find that when the three-drug combination and one of the pairwise combinations are both lethal, the rescaled E3 measure identifies the emergent interaction as either antagonistic buffering or lethal synergy. This situation is consistent with both definitions. However, unless the three-drug combination is chosen to be non-lethal, the effect of the third drug on the outcome is not obvious as the pairwise combination already represents lethality. For this reason, we identified all such cases as inconclusive, which is consistent with previous work on large drug interaction networks (e.g. Yeh *et al*. [8]).

### 2.2. Cut-off values for the rescaled emergent interaction measure

Based on previous work [8,14] and the analysis of the resulting distributions of [E3]_{R2} over all the drug triples (see Model Framework and Results), the interaction is identified as synergistic when rescaled E3 is less than −0.5, antagonistic if it is greater than 0.5, and additive otherwise. For antagonistically identified triples, we further choose 1.3 as the cut-off between the special cases of buffering and suppression.

### 2.3. Model framework

#### 2.3.1. Three-way interaction measures

In general, a measure of the efficacy of a treatment is how much it inhibits bacterial growth rate relative to growth in the absence of drugs. This measure is equivalent to a relative fitness for the bacteria that is typically denoted by *w*_{D}, where D stands for a single drug or a mixture of several drugs. Here, the fitness measures are symmetric in the ordering of drug indices, for example, *w*_{XY} = *w*_{YX} for drugs X and Y. As discussed in the Introduction, BI is when the per cent reduction in growth rate by a single drug is independent of the presence of other drugs and is expressed in an equation as *w*_{XY} = *w*_{X}*w*_{Y} [6]. Accordingly, the DA is defined by with the general interpretation that a large-enough negative value of DA_{X,Y} implies synergy between drugs, such that the combined effect is greater than would be predicted based on the single effects. Conversely, a large-enough positive DA_{X,Y} means that the drugs are acting antagonistically.

Identifying the existence of some type of interaction among more than two drugs can be defined analogously (i.e. via generalization of DA). For three drugs, the DA measure becomes: , quantifying interactions at any level that contribute to the overall interaction [14,21,28,29]. With more sophisticated modelling and measures, it is also possible to identify true emergence—the overall interaction is not just a result of interactions among subsets of the drugs. That is, it is important to distinguish between effects that arise from lower order interactions (such as pairwise interactions that yield an apparent three-way effect) and those that arise from emergent interactions that require all of the drugs to be present to manifest their unified effect in killing bacteria. A recent model that is capable of making this distinction was introduced [14] and was termed the E3 measure.

The logic of the E3 measure is that all possible pairwise contributions are correctly weighted and subtracted from the overall interaction; hence it quantifies any triple-drug interaction that does not originate from the pairwise interactions. The weighting is determined based on the expected three-way effect when only two drugs (pairwise interaction X and Y) interact and the third drug (Z) is additive with them. Based on the three-way DA measure, this would give an effect of *w*_{z}DA_{X,Y} because *w*_{Z} would factor from all terms and the remaining terms are just the definition of the pairwise DA measure. For three drugs, there are three pairwise combinations, so subtracting all possible pairwise combinations yields
2.1

Rewriting this equation solely in terms of the relative fitness gives 2.2

Notably, this E3 measure includes every possible relative fitness in the three-drug system. Both the E3 and DA measures are symmetric in drugs X, Y and Z and can be easily adapted to higher order interactions that involve more than three drugs. More detailed discussion of the derivation of higher order interaction metrics can be found in Beppler *et al*. [14].

#### 2.3.2. Rescaling three-way interaction measures

To easily quantify the interaction strength and the separation of interaction classes, the interaction measure must be rescaled. For pairwise interactions, Segre *et al*. [18] established a rescaling method that greatly enhances the discovery of antagonism. Their normalization is based on two limiting reference cases for the synergistically and antagonistically interacting drug pairs. The normalization factor for synergy is defined by substituting the lethal case into the DA measure (i.e.), whereas the antagonistic interaction is rescaled by the complete buffering case in which the two-drug effect is the same as the fitness of the single drug with the stronger effect (i.e. ). With this rescaling the pairwise interactions yield a trimodal distribution centred on these reference cases (modes at −1, 0 and 1 as expected theoretically) and provides a clear cut-off between synergistic, additive and antagonistic interactions [8,18,20,35].

To classify interactions among three-drug combinations, we also need to establish an appropriate rescaling or normalization. In particular, the reference cases for negative and positive measures must be properly defined for each triple-drug interaction measure. It is expected that there is a biologically and empirically grounded choice of rescaling that will provide a clear distinction between the interaction classes, as in the case with two drugs. Moreover, as our rescaling measures are based on special cases of the DA or E3 measures, they are symmetric by construction.

The case of lethality is uniquely defined for each three-way interaction measure. On the other hand, buffering can be defined in several different ways, and the choice of definition may depend on the type of interaction measure (emergent versus overall) being considered.

Because DA quantifies the deviations from the case that all single-drug effects are combined additively, it is practical to define buffering with respect to the effect of the strongest individual drug. However, the definition of buffering for E3 is not unique because E3 captures the effect of the three-drug combination relative to all lower order effects. Hence, buffering could be defined relative to the strongest individual drug, the strongest pairwise effect or the strongest of all of these. As a result of this ambiguity, the choice of rescaling might cause misidentification of antagonism. In the subsequent sections, we first present the extension of the two-drug rescaling method to three-way interaction measures, and then construct new rescaling factors that will help to more comprehensively evaluate and characterize the emergent properties of drug combinations.

#### 2.3.3. Extension of two-drug rescaling method (Rescale 0)

Beppler *et al*. [14] and Sanjuán *et al*. [21] introduced a rescaling that directly extends the method of Segre *et al*. [18]. Negative measures are scaled relative to extreme lethal synergy such that the triple-drug combination kills off the bacteria, i.e. *w*_{XYZ} = 0. When the interaction measure is positive, it is scaled with respect to the buffering case defined to be when the fitness of bacteria exposed to the triple-drug combination is the same as when exposed to just the single drug with the strongest efficacy, i.e. . Hence, the normalization factor (Scale 0, denoted by S0) for the positive E3 measure is defined by
Then, the rescaled E3 measure with S0 (Rescale 0) takes the form
Based on the discussion above, the current published rescaling method [14,18,21] is the most appropriate way of characterizing interaction type based on the DA measure and the synergistic emergent interactions.

#### 2.3.4. New rescaling methods for emergent three-way interaction (E3)

To explore alternative rescaling methods that may lead to a clearer identification of antagonistic interactions, we revisit the possible definitions of buffering. Based on the buffering definition, the antagonistic buffering cases would be mapped to 1 via the corresponding rescaled metric. Note that the synergistic interactions are characterized by the previously defined rescaling method; hence rescaling methods presented here map the lethal synergies to −1. All of the different rescaling methods are summarized in table 2.

#### 2.3.5. Rescale 1: buffering relative to pairwise drug effects

Higher order combination therapies of *N* drugs are analysed relative to the lower order combinations (i.e. *N* − 1, *N* − 2, etc.) of subsets of the drugs. Therefore, another type of buffering is when the effect of all three drugs combined is exactly the same as that of the most powerful pairwise combination (i.e. the lowest relative fitness in the presence of any pairwise combination of drugs). Consequently, for antagonistically interacting drugs (i.e. E3 > 0), we introduce the normalization factor S1 by substituting the minimum of all pairwise drug fitnesses for the triple-drug fitness in E3 (table 2):

Hence, the rescaled E3 measure (Rescale 1) is

In this way, we can assess if any of the pairwise drug combinations disguise the effect of the cooperation and antagonism of the remaining drug in the triple-drug combination.

#### 2.3.6. Rescale 2: buffering relative to single and pairwise drug effects

In this section, we propose a generalization of the above S1-rescaling scheme that accounts for *all* lower order drug components. This generalization allows for the possibility that a single-drug therapy might be more powerful and offer more effective treatment than any pairwise combination (i.e. the lowest relative fitness in the presence of a single drug is lower than the relative fitness in the presence of any pairwise combination). This extension defines the E3 measure relative to the minimum of the fitness values attained by any single or pairwise components as (table 2)

Thus, the new rescaled E3 measure (Rescale 2) takes the form

Rescale 1 and 2 yield the same results if the minimum of the pairwise drug fitnesses is less than the minimum of the single-drug fitnesses, which is frequently the case. Notably, Rescale 2 spans a broader set of circumstances than Rescale 1.

We further note that it is straightforward to generalize Rescale 2 when emergent properties with more than three drugs are analysed. For the case of *N* drugs, the scaling factor S2 is given by calculating the emergent measure by replacing the *N*-drug fitness by the minimum of the bacterial fitness among all possible subsets of the drugs. For instance, four-drug combinations are rescaled relative to the case when

#### 2.3.7. Rescale 3: buffering relative to effects from pairwise interactions

Because the E3 measure captures the emergent triple-drug interaction that cannot be explained by the pairwise effects, we consider another definition of buffering to be relative to the expected pairwise contribution to the overall three-way interaction (DA). Note that based on the DA and E3 measures, this is not the same as the relative fitness in the presence of a pairwise drug combination because each of those relative fitnesses is multiplied by the proper weighting (i.e. the relative fitness in the presence of the remaining single drug). Replacing the DA term (DA_{X,Y,Z}) in the E3 expression (2.1) by the minimum of pairwise interaction contributions defines our last scaling method, S3, via

Substituting the pairwise fitness measures into the minimized quantity in S3, we obtain which is a symmetric expression in terms of the pairwise contributions. S3 can be expressed purely in terms of the fitness parameters as

So when the unscaled E3 measure is positive, we rescale the E3 measure (Rescale 3) as (table 2)

Because experimental concentrations of individual drugs are chosen such that each reduces the growth about 15–35%, we expect somewhat similar results as in the Rescale 1 because we are simply multiplying by a number between 0.65 and 0.85 (relative fitness of a single drug).

## 3. Results

To distinguish different types of interactions manifested in triple-drug combination therapies, we use a newly proposed emergent three-way measure [14]. To include as much information as possible, we show all data in the histogram results (figure 2). We observe that the unscaled version of the emergent measure yields a unimodal distribution around the additive case (figure 2*a*). We find that Rescale 0 which defines buffering based on the comparison of the triple-drug combination with the single drug effects (previously used method in [14,21]) results in a distribution with modes at −1 and 0, hence identifying emergent synergies. Yet, Rescale 0 does not result in a clear distinction for the antagonistic triples (figure 2*c*).

Importantly, using the new rescaling methods we propose here, we see that Rescale 2, which defines antagonistic buffering relative to the single or pairwise drug effects, successfully maps (figure 2*b*) lethal synergies to −1, additive interactions to 0 and the antagonistic buffering interactions to 1. Consequently, the resulting distribution from Rescale 2 leads to a multimodal distribution with three dominant peaks (see Silverman test results for multimodality in electronic supplementary material, table S1 and figure S2). The distribution of Rescale 2 without the inconclusive cases—when one of the pairwise combinations and triple combination are lethal—is given in electronic supplementary material, figure S3. Although the peak at 1 is diminished when the inconclusive cases are excluded, our main result still holds as the identification of synergistic and antagonistic interactions is enhanced. Note that the same result is obtained by Rescale 1 where the triple combination is only compared with the most effective pairwise combination (figure 2*d*). This is because the pairwise components in the analysed data are always stronger than individual drugs. Moreover, in agreement with theoretical predictions, Rescale 3, which defines antagonistic buffering relative to the most powerful pairwise interaction, yields similar results as in both Rescale 1 and Rescale 2. Notably, Rescale 3 leads to cleaner separation between the different modes than the other choices of rescaling (figure 2*e*).

Via Rescale 2, we find 38 synergistic, 78 antagonistic buffering and 47 antagonistic suppressive emergent interactions from 364 different triple-drug combination experiments. Of the remaining combinations, 165 are additive and 36 are inconclusive due to lethality to the pairwise and triple combination. Interaction classifications found via each of the rescaling methods are summarized in the electronic supplementary material, table S2. In figure 3, we provide a comparison between the existing rescaling method (Rescale 0) and the newly proposed method Rescale 2 in terms of the identification by each of the synergistic or antagonistic triples. We find that the majority of emergent antagonisms (90%) are either classified as additive or underestimated (buffering versus suppression) according to Rescale 0. A comparison of emergent interaction types obtained via the other rescaling methods is given in electronic supplementary material, table S3.

We additionally give two examples of highly synergistic and highly suppressive emergent interactions. First, a triple-drug combination of ciprofloxacin + clindamycin + streptomycin yields emergent lethal synergy (figure 4*a*) according to all rescale methods studied here, whereas such an interaction is not apparent with the unscaled version of E3. Next, the drug combination with erythromycin + cefoxitin + tobramycin shows highly suppressive emergent interaction (figure 4*b*) according to rescaling methods 1–3, while it is barely identified as antagonistic according to Rescale 0 (previously used method) and is not even identified as antagonistic with the unscaled version.

To see if we obtain similar results (e.g. unimodal distribution from unscaled data, multimodal distribution using scaled methods) with other datasets, we also analysed the data from Wood *et al*. [28]. We find that the unscaled measure for this dataset also yields a unimodal distribution (figure 5, inset), seeming to imply that almost no emergent interactions exist. However, when we apply rescaling methods to this dataset, we observe multimodal distributions (figure 5; electronic supplementary material, figure S4) that are similar to the results from the analysis of our own data. The peak at –1 (corresponding to lethal synergies by definition) is not visible for these data. This is due to the fact that the vast majority of the data (93%) are coming from just six combinations of three drugs that are varied in their concentrations to create a larger set of data. These particular combinations and the range of concentrations used result in either emergent antagonism or emergent additivity for almost all experiments. These results are consistent with previous reports that different drug concentrations do not change a drug interaction from antagonistic to synergistic [8]. Therefore, the data in this study mainly represent the effect of varying drug dosages for a small set of drug combinations, whereas our study measures a much larger set of drug combinations and thus analyses a much larger set of potential interactions and classifications. Even for the limited set of drug combinations studied in Wood *et al*. [28], we find our rescaling methods greatly enhance the identification of interactions. Specifically, we find only 53% of the measurements are additive with the remaining measurements being synergistic (3%) or antagonistic (44%) emergent interactions, as opposed to more than 97% additivity reported by Wood *et al*. [28]. That is, the choice of rescaling leads to dramatically different conclusions about the presence and overall prevalence of emergent interactions. Moreover, our new analysis reveals that the Wood *et al*. [28] data also frequently exhibit emergent interactions and that these interactions tend to be much more antagonistic, and using our new rescale, different datasets yield the same conclusions. Consequently, applying our new rescaling methods to different datasets (ours and Wood *et al*.'s [28]) yields similar conclusions about an increase in antagonism for higher order emergent interactions, suggesting the generality of our findings and the need for rescaling.

## 4. Discussion

Through a systematic analysis of relative growth rates of bacteria exposed to single, pairwise, and triple-drug combinations, we introduced rescaling methods for clearly delineating and categorizing types of interactions among multiple drugs. Owing to its unimodal shape, we concluded that the unscaled emergent measure is not useful for identifying and distinguishing among drug interactions. We overcame this problem by constructing new rescaling methods relative to natural reference frames of synergistic and antagonistic cases. As a result, our rescaling methods of the emergent measure allow us to clearly distinguish interaction types via multimodal characteristics of the emergent E3 measure distribution.

We have shown that defining antagonistic buffering as the masking effect by the pairwise combinations (Rescale 1) or by all possible lower order combinations (Rescale 2) lead to nearly identical results. Intriguingly, Rescale 3, where buffering occurs when the overall interaction effect is masked by the strongest pairwise interaction, is somewhat similar to Rescale 1 and 2 but yields important differences for categorizing interactions. The consistency of all three new rescaling results gives confidence that the identification of the emergent properties is well established. Among these methods, we identify Rescale 2 and Rescale 3 as the best at categorizing emergent interactions based on the distributions they produce. Rescale 2 may be more widely adopted because of ease of interpreting from the fitness measurements, hence comparing through standard bar graphs. On the other hand, Rescale 3 may be desirable because it is more naturally connected to the underlying formulation of E3 itself.

There are three baselines for rescaling: −1 (lethal synergy), 0 (additivity) and 1 (antagonistic buffering). Performing a test of multimodality for the distribution of E3 with Rescale 2, we find at least four modes, with three dominant ones occurring at −1, 0 and 1 (electronic supplementary material, table S1). These modes suggest that our empirical data do segregate into the three extreme cases upon which our rescaling metrics are designed. Intriguingly, there may be a fourth peak that could be a weaker form of antagonism (see multimodal Silverman test in electronic supplementary material, figure S2). This mode is consistent with features of the pairwise distribution of interactions although that limited dataset makes it difficult to see the peak [8]. As we obtain more empirical data, this fourth peak and other characteristic features of this distribution will be explored in more detail.

Applying the rescaling methods we proposed here substantially changes the existing results for the multi-drug interactions from Wood *et al*. [28]. They suggest that most of the interactions in the multi-drug environment can be explained by pairwise and single-drug effects. However, via the use of our rescaling methods, we identify substantial numbers of non-additive interactions from the same dataset. Hence, we conclude that for many combinations the higher order effects are not simply the result of lower order effects. These conclusions further strengthen the importance of the choice of an appropriate rescaling method when identifying higher order emergent interactions.

Our rescaling method enhances the mode for emergent synergy (at −1) even more than the Segre *et al*. [18] rescaling method enhances the mode for synergy for pairwise interactions in a yeast epistasis study. This can be explained by considering the differences of the fitness ranges for yeast knockouts and drug combinations. The fitness of the single knockout is typically very close to the wild-type fitness (=1) [18]. As interaction type is determined based on the comparison between single and double mutants, lethal synergies (when fitness of double knockout is close to 0) are obvious even with the unscaled measure. However, drugs at most useful concentrations lead to lower fitnesses than those that result from single-gene knockouts. Because of this reduced difference between lethality and the single effects, lethal drug synergies are not as obvious with the unscaled measure as they were in the gene epistasis studies. Intriguingly, the rescaling of the emergent interaction measure helps uncover the strength of the synergistic interaction. Moreover, in addition to the modes expected for synergistic, additive and antagonistic buffering interactions, the rescaling for emergent drug measures yields a substantial number of suppressive triples, thus yielding even more useful information for classifying drug interactions. Taken together, our new rescaling method offers a strong and robust identification and categorization scheme for three-drug interactions and very possibly higher order interactions in general.

The identification of highly synergistic and suppressive emergent interactions that we established in this study could be especially important for antibiotic research. This importance and clinical relevance is because highly synergistic triple-drug combinations are of utmost importance due to the high efficacy of the treatment compared with pairwise interactions [1,36–40]. On the other hand, identifying suppressive interactions may be especially valuable because it has been shown that, counterintuitively, these interactions may slow and thus suppress the evolution of antibiotic resistance [8,15,17,41–43]. That is, there may be a trade-off between killing efficiency and the evolution of resistance. Synergistic combinations may be especially good at killing bacteria but may also increase the likelihood and rate of evolution of resistance [17,42–44]. Thus, this trade-off can also be seen as balancing the good of an individual with the good of the public.

We note that it can be very important clinically whether three-drug interactions are different from the expectations based on all two-drug interactions. This is because an emergent synergistic interaction would be a potentially novel therapeutic option, whereas a three-way synergy that just resulted from a pairwise synergy would not gain much benefit by adding the third drug. Indeed, in such a case, the third drug could be irrelevant in terms of killing bacteria, and thus likely should not be used in order to decrease the toxic side effects to the patient. This rescaling method will allow us to determine if an additional third drug gives a real benefit. As mentioned above, drug interactions can change the topography of fitness landscapes, and this is important because ultimately, it would be very useful to design landscapes such that the bacteria would have a difficult time reaching maximal fitness peaks with multi-drug resistance. Emergent interactions, as opposed to DA interactions, could change the shape of the fitness landscape more dramatically.

One caveat is that the empirical data are for fixed drug concentrations, so that BI is used to determine interaction type. BI is a much simpler and straightforward way of measuring interactions than Loewe additivity, which requires examination of interactions across a gradient of concentrations for each of the drugs. Although Loewe integrates more information about interactions, the ease of measurement and calculation of BI has led to its use in a huge number of studies [8,14,21,29,45–49], including our present one.

Methods developed here for capturing the higher order interactions are applicable to other complex systems. These systems include, for example, protein and gene interaction networks (e.g. [22]), food webs (e.g. [23–25]) and transportation networks (e.g. [50]). Hence, our rescaled emergent measure offers a systematic and straightforward method for uncovering the complex interactions that occur in a wide range of systems.

In summary, we have shown that our new rescaling methods lead to a clear distinction between different categories—synergistic, additive, antagonistic buffering and antagonistic suppressive—of interactions, and that these distinctions do not exist when looking at results for the unscaled measure or the previously used rescaling method (i.e. Rescale 0). Therefore, our new rescaling methods are required to accurately characterize cases of emergent interactions among multiple drugs. In particular, Rescale 2 and Rescale 3—synergy defined relative to lethality for both cases and antagonistic buffering relative to the strongest lower order effect (Rescale 2) or interaction (Rescale 3)—offer a direct and more distinct measurement of epistasis that is straightforward to generalize to higher order interactions. Therefore, we propose that Rescale 2 and 3 could be a useful tool in future studies that examine complex drug interactions or other complex systems with interacting components.

## Data accessibility

The datasets for this article have been uploaded as part of the electronic supplementary material.

## Authors' contributions

E.T., V.M.S. and P.J.Y. conceived of the project and designed the study; E.T. and V.M.S. developed new analytical tools; C.B., C.W. and Z.M. conducted the experiments; E.T. and V.M.S. conducted the analyses; and E.T., C.B., V.M.S. and P.J.Y. wrote the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

This work was funded by a James F. McDonnell Complex Systems Scholar Award to V.M.S., NSF DBI Career Award to V.M.S. and UCLA Faculty Career Development Award to P.J.Y.

## Acknowledgements

We thank Elliot M. Landaw and Tom Chou for discussion, and Tina Manzhu Kang, Emily Vargas, Borna Shirvani and Gabriel Ting for assistance in the laboratory.

- Received April 27, 2016.
- Accepted May 17, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.