Statistical determination of synergy based on Bliss definition of drugs independence

Summary

Although synergy is a pillar of modern pharmacology, toxicology, and medicine, there is no consensus on its definition despite its nearly one hundred-year history. Moreover, methods for statistical determination of synergy that account for variation of response to treatment are underdeveloped and if exist are reduced to the traditional t-test, but do not comply with the normal distribution assumption. We offer statistical models for estimation of synergy using an established definition of Bliss drugs’ independence. Although Bliss definition is well-known, it remains a theoretical concept and has never been applied for statistical determination of synergy with various forms of treatment outcome. We rigorously and consistently extend the Bliss definition to detect statistically significant synergy under various designs: (1) in vitro, when the outcome of a cell culture experiment with replicates is the proportion of surviving cells for a single dose or multiple doses, (2) dose-response methodology, (3) in vivo studies in organisms, when the outcome is a longitudinal measurement such as tumor volume, and (4) clinical studies, when the outcome of treatment is measured by survival. For each design, we developed a specific statistical model and demonstrated how to test for independence, synergy, and antagonism, and compute the associated p-value.

Keywords:

Drug therapy – Cancer treatment – Normal distribution – Breast tumors – Drug administration – Drug screening – Synergy testing – Drug interactions

Introduction

The definition of synergy is one of the most controversial concepts in biology and medicine. Despite its fundamental importance in pharmacology, toxicology and experimental medicine, including cancer research, many, sometimes contradictory, definitions of synergy and drugs interaction exist in the literature [1], [2].

There are two main competing definitions of drugs independence: based on (a) Loewe [3] definition of additivity (isobologram approach) and implied combination index (CI) and (b) Bliss definition of independence [4] or, equivalently, Webb [5] fractional product. Herein, we discuss the advantages and limitations of these two major definitions, and then apply Bliss definition of independence for determination of statistically significant synergy in popular experimental and clinical trial settings in biology and medicine. Although several statistical packages in R for determination of synergy exist, such as hbim [6], mixlow [7], COMBIA [8], and CImbinator [9], they lack rigorous statistical testing and computation of the p-value for various treatments/drugs interaction designs under one methodological umbrella, as we propose in the present work.

Advantages and limitations of Loewe additivity and combination index

According to Loewe, in the case of no interaction between drugs A and B, the drug combination with doses d_A and d_B leading to the same mortality M must satisfy the so called “median-effect” equation

This equations defines a segment on the plane with dose concentrations (d_A, d_B) splitting the positive quadrant into two parts corresponding to synergy and antagonism, and as such is visually attractive. The CI can be computed using the popular CalcuSyn software available at http://www.biosoft.com/w/calcusyn.htm [10] using the following steps: (a) conduct n experiments with drugs A and B alone at a set of doses {d_Ai, i = 1, …, n} and {d_Bi, i = 1, …, n}, and the same doses in combination when the drugs are given simultaneously {(d_Ai, d_Bi), i = 1, …, n}, (b) construct individual dose-response curves for each drug, (c) for each mortality value M with the drug combination find equivalent single-agent doses D_1A and D_1B by inverting the individual dose-response relationships from the previous step, and finally (d) compute n values of the left-hand side of Eq (1) and report the mean and standard deviation of the CI.

Although visually attractive, this approach has several limitations:

• Statistical testing of drugs independence based on CI or Loewe definition (2) is troublesome because the uncertainty estimation of D_1A and D_1B is not taken into account. Moreover, the suggested grading of synergy based on the CI value, as suggested in [11], ignores the uncertainty of the CI estimation. A more comprehensive and scientifically accepted way to test for drug interaction is to use methods of statistical hypothesis testing through computation of the p- or d-values [12], [13]. Since these statistics are reported, synergy may be claimed when the drugs interaction is not statistically different from being independent. In short, CI <1 is not sufficient to claim synergy due to inevitable stochastic nature of experimental data.

• The criterion for independence (1) requires knowledge of single-drug dose-response relationships, and particularly in the case of (2) knowledge of EC₅₀. This requirement is often not satisfied when basic preliminary cell culture experiments are performed using just a few drug concentrations. Simply put, individual single-drug dose-response relationships are often not available or derived just using few data points and as such prone to large uncertainty not taken into account when computing the CI.

• Why is the killing benchmark set at 50%? Although EC₅₀ is an established parameter in toxicology it may be not appropriate in other applications where the point of interest is the tail of the dose-response curve. For example, epidemiology is concerned with pollutant concentration with relatively small mortality effect. On the other hand, cancer researchers mostly look at the response to the drug where cancer cell mortality is close to 1. Unlike Loewe additivity, Bliss independence does not use EC₅₀ as the benchmark of the drug effect and therefore applicable on the entire range of doses.

• It is well-known that the response to treatment is better to model on the log scale. Indeed, the authors of the CI approach use the log-transformed data for the single-agent dose-response relationships [14], [15], [16], [17], [18], yet drugs in Eq (1) are expressed on the linear scale. The Loewe additivity model assumes that drugs interplay in the linear fashion and as such has neither theoretical nor empirical foundation. We agree with Geary [19] who noted that “This metric is ad hoc based neither on physiology nor on formal axioms.”

• CI is difficult to generalize to many other drug interaction settings such as cell culture in the presence of an affected control group, longitudinal organism experiments such as tumor growth in animals, or survival in clinical studies.

Most authors who discuss drug interaction provide definitions, but few publications address the topic of statistical determination of drug interaction. Data on response to treatment derived from experiments/trials in biology and medicine are subject to considerable variation. Simply reporting a CI is insufficient to claim synergy because CI < 1 may happen due to natural biological variation and stochastic nature of experiment outcomes. We assert that determination of synergy is incomplete without rigorous statistical testing against the null hypothesis of drugs independence and reporting the respective p-value.

The basis for our statistical developments is Bliss definition of drugs/treatments independence. We adjust this definition to various designs in biology and medicine and show how to test for synergy with maximum statistical power.

Probabilistic interpretation of Bliss independence

The goal of this section is to illustrate Bliss independence by the formula for the probability of survival in independent events where the proportion of survived cells in assay is treated as the probability of surviving of a single cell.

According to Webb [5], drugs A and B act independently if the surviving fraction (SF) of cells or organisms upon simultaneous administration, denoted S_AB, is equal the product of SFs, S_A and S_B, when drugs are given separately, or in mathematical language,

We say that there is synergy if S_AB < S_AS_B and antagonism if S_AB > S_AS_B. This definition of independence has a clear interpretation: after treatment with drug A the proportion of cells in the population is S_A. When drug B is subsequently given after drug A, drug B affects only living cells; therefore if drugs act independently, the proportion of living cells upon simultaneous drug administration is S_AS_B.

This definition of synergy has straightforward probabilistic interpretation [20]: associate surviving fractions S_A and S_B in the population of cells with probabilities of death of a single cell due to single-agent administration of either drug A or B, respectively. If the events of survival are denoted as A and B, then their probabilities are Pr(A) = S_A and Pr(B) = S_B. If the drugs act independently, the probability of survival (as follows from the probability theory, [21]) upon simultaneous administration of both drugs is the product of individual probabilities:

The independent action of drugs is illustrated in Fig 1 using the Venn diagram.

The definition of independence (3) expressed in terms of survival is equivalent to Bliss [4] independence expressed in terms of mortality,

The definition of drugs independence expressed in Eq (4) can be explained as follows: the sum M_A + M_B is the the total number of cells killed when two drugs act independently, but it contains double-counted cells M_AM_B that must be subtracted, as cell can only die once.

Below, we systematically apply the definition of drugs independence either in the form of survival (3) or equivalently in the form of mortality (4) to various settings of drug interaction in vitro and in vivo for statistical determination of synergy or antagonism.

Treatment interaction with single dose

This section deals with a simple yet fundamental study design for investigation of treatment (e.g. drug) interaction at specific dose. As above, the effect of a drug is associated with a surviving fraction (SF) of cells or organisms. It is critical to recognize the importance of replicates in each drug group to account for natural variation of response to treatment to make a statistically-informed statement about possible treatment independence, synergy, or antagonism. In many previous synergy analyses, such as those implemented in CompuSyn software mentioned above, replicates were reduced to averages by which the statistical variation of the treatment effect has been down-played. In contrast, our statistical approach to synergy is more realistic because it addresses the inevitable variation of response to treatment in cells or organisms, both and between treatment groups [20].

We underscore the importance of the single-dose assay as a preliminary step to further employment of a full-range dose investigation because it requires fewer samples. More fundamentally, synergy may exists for a specific pair of treatment doses and not an entire range of doses, therefore, a single-dose synergy may be overlooked in the traditional whole-curve dose-response approach. We emphasize that CI cannot be used with single-dose experiments because EC₅₀ is not available.

Treatment interactions with replicates

Let there be n₁ replicate experiments in drug group A, n₂ replicates in drug group B, and n₃ replicates in drug group D when drugs with the single dose from groups A and B are given simultaneously. Since the observed SFs are positive it is convenient to model the variation of SF on the log scale, which can be expressed through the exponential function as follows:

The true unknown μs are expected to be negative and connected to unknown true surviving fractions from the previous section as S A = e μ 1, S B = e μ 2 and S A B = e μ 3 . According to Bliss, drugs act independently if

We will test this hypothesis by statistical means using replicates in each group.

An advantage of representation (5) is that after log transformation, the model turns into an analysis of variances (ANOVA) model. Specifically, let the log fractions in each group be y_1i = ln A_i, y_2i = ln B_i and y_3i = ln D_i. Then the system (5) can be rewritten as an ANOVA model [13], [22], [23] with three groups,

This hypothesis is tested using the test statistic T which has a t-distribution with ∑ i = 1 3 n i - 3 degrees of freedom, namely,

Other than log transformation can be applied to surviving fraction to comply with the normal distribution assumption. For example, one may use the logit transformation, ln (S/(1 − S)), or inverse Weibull cumulative distribution function transformation [27]. The advantage of the log transformation is that the Bliss definition of independence turns into a linear equation easy to test by standard linear statistical theory.

Several authors used ANOVA models to test Loewe additivity condition based on the CI, or test the Bliss independence null hypothesis H₀: S_AB − S_AS_B = 0 using the traditional t-test. Particularly, the CI-based ANOVA models have been criticized in [2] and [28]. We have to comment that the normal distribution assumption is not adequate here because it makes the possibility of CI or surviving fraction taking negative values, which makes no sense. We, however, test the Bliss independence on the log scale, as follows from Eq (6), and therefore normal distribution and negative values are permissible.

Example: Daphnia acute test with CuSO₄ and NiCl

We illustrate treatment interaction with a single dose of each treatment using 48-hour acute tests with Daphnia [29], [30]. In each experiment, organisms in water were treated with two stressors, a single dose of NiCl (1.8 mg/L), CuSO₄ (7.0 μg/L) or the combination, and the numbers of surviving organisms were counted after 48 hours. There were three repeated experiments (replicates) with NiCl, four replicates with CuSO₄ and four replicates with the combination. The results are portrayed in Fig 2. The right-most group “Independence” was derived from the individual groups by computing SF under the assumption that treatments act independently, i.e. on the log scale as pairwise sums y_1j + y_2k, where j = 1, 2, 3 and k = 1, 2, 3, 4. Note that although all calculations with SF were done on the log scale we used the original scale in % on the y-axis for easy interpretation (the scale is not uniform). For these data, y ¯ 1 + y ¯ 2 - y ¯ 3 = 0 . 156 and the ratio of SF under independence to SF when both treatments are applied simultaneously is exp(y₁ + y₂ − y₃) = 1.168 which means that synergy was 17%, which means that 17% more cells would survive if stressors acted independently. However, according to the t-test (7) this result is not statistically significant because p = 0.172—there is no sufficient evidence to reject the null hypothesis that NiCl and CuSO₄ act independently.

It is worthwhile to note that determination of statistical significance of drug interaction is impossible using the isobologram and Combination Index (CI) approaches because (a) replicates are substituted with average values so that variation between replicates is ignored and (b) the EC₅₀ values are not available because the data exist only for single dose.

Treatment interaction in the presence of control group

In many biological experiments treatments, elicit delayed but not immediate effects on cell/organism. During such time, cells may die or proliferate by natural means, which can be captured through inclusion of a control/sham group C and the associated surviving fraction S_C. Following our probabilistic approach we equate the proportion of cells that survived in group C with the probability of dying of a single cell, S_C = Pr(C). Now we turn our attention to the group treated with drug A. Appealing to probability theory, we treat fraction S_A as the proportion of cells killed by drug A but exclude those cells that naturally died as in the control group, S_A = Pr(A ∩ C). Thus, the complimentary probability 1 − S_A is the probability of death either due to drug A or natural causes. The probability of survival associated with the effect of drug A can be viewed as the conditional probability [21]

Note that this formula connects the treatment effect of drug A with the observed surviving fraction in groups A and C. We do not observe the conditional event A|C, but treat this as imaginary occurrence to account for cell loss due to natural death, or reversely, to cell proliferation observed in control group. A similar probability of survival can be obtained for the population of cells treated with drug B alone, and for simultaneous treatment with both drugs A and B,

Following the rule of independence defined in Eq (3) we arrive at a more general definition of synergy when the population in the control group reduces due to natural cell death,

The statistical model for testing hypothesis (8) is similar to what was described earlier with the addition of replicates from the control group C j = e μ 0 + ε 0 j . Denoting log SF as y_0i = ln C_i the system can be rewritten as an ANOVA model with four groups as follows:

Under this hypothesis the ratio of the mean difference to the standard error

Slinker [31] defines synergy in terms close to (9) and urges use of ANOVA to test the independence of treatment interaction. However, Slinker was not specific regarding what biological quantity/effect to use for μ and simply referred to “experimental outcome.” Instead, (a) we apply treatment interaction to the data when the outcome of an experiment is survival fraction, (b) we derived the null hypothesis using the conditional probability, and (c) as a part of our formal derivation we use the log scale i.e. μ is the log SF and therefore comply with the normal distribution assumption.

Example. Testing for drug synergy for cancer cells

These cell culture data come in the form of the surviving fraction of cells per dish. These experiments were intended to test for synergy between the drugs BYL and GSK in ZR75-1 breast cancer cells in vitro [32]. There were four treatment groups and three replicates per group: control, 1μM BYL, 1μM GSK, and the combination at 1μM each. For each treatment group, we first computed 9 pairwise differences from the control group on the log scale as ln A_j − ln C_k, ln B_j − ln C_k, and ln D_j − ln C_k; see Fig 3. Secondly, data were plotted with the y-axis as % surviving cells with respect to the control group. The analysis and graphics representation are similar to the Daphnia case with no control group. The right-most box depicts the imaginable surviving fraction (SF) if drugs acted independently according to Bliss. Those values were computed based on 27 triple combinations from groups A, B, and C as follows: (ln A_j − ln C_k) + (ln B_l − ln C_k). If BYL and GSK acted independently the two right-most boxes would be the same height. However, the SF under the additive assumption (right-most box) is higher than that from the actual data (combination group); therefore, we may expect treatment synergy between BYL and GSK. Indeed, the ratio of SF under Bliss independence to SF from the combination group is 1.65, which can be interpreted as 65% synergy. The problem is that the SFs exhibit variation and the T-statistic computed by formula (10) yields 1.7 producing p = 0.129 that indicates that observed synergy is not statistically significant.

Thus, despite a fairly high value of synergy, the p-value is not below 0.05 because of the variation in the group surviving fraction.

Dose-response relationship

In this section, we discuss statistical determination of treatment interaction when cell/organism survival experiments are conducted using a range of doses. Two types of analyses are presented: (a) model-free statistical determination, when survival experiments with the individual and combination treatments are dose-escalated in parallel, and (b) estimation of a bivariate dose-response relationship with any set of doses.

Model-free statistical determination of synergy

This method was originally suggested by Webb [5] and is commonly referred to as the fractional product method [1]. The idea is compare the SF from the combination group with the SF computed as the product of the SFs from the two single-agent groups. Although the idea is straightforward, no rigorous statistical model or statistical test for synergy has been developed. For example, [33] suggest estimating the beta-slope coefficient in the linear regression of the SF in the combination group on the product of the SFs from single-agent treatment groups. However, this approach has a number of limitations: (i) it does not comply with the normal distribution assumption because the SF is in the range from 0 to 1, (ii) it assumes that the single-agent SFs do not have variation (fixed) but the SF from the combination group does, (iii) no statistical test for the Bliss independence hypothesis H₀: β = 1 has been developed. In contrast, statistical methods described below comply with the normal distribution because the analysis of the SF is carried out on the log scale and rigorous statistical tests of Bliss independence are developed for different assay designs.

We propose to test Bliss independence by employing the theory of a linear model similar but not equivalent to the ANOVA model discussed above. Two types of design are studied here: incomplete and complete pairwise design. In incomplete design, the same number of doses are used in individual and combination treatments. In complete design, all pairwise combination of treatments are required. An advantage of the model-free statistical determination of synergy is that no dose-response relationship is required, but experiments with combined treatments must be included with the same concentrations as individual treatment. Although the models described below do not account for an affected control group, their generalization is straightforward following the idea in the preceding section.

Incomplete pairwise design

Let S_Ai and S_Bi be SFs of cells in the ith experiment with dose d_Ai of drug A and dose d_Bi of drug B for i = 1, 2, .., n. Let S_Di be the SF in the ith experiment when drugs A and B are combined with doses (d_Ai, d_Bi). If drugs are acting independently ln S_Ai + ln S_Bi and ln S_Di must be close. To test the difference, the following statistical model under multiplicative error scheme is suggested. Denote x_i = ln S_Ai, y_i = ln S_Bi, and z_i = ln S_Di and let

As an example, we illustrate this test using the data from [34]. For these data, the estimate of the delta-parameter is 1.64 which indicates that there is antagonism between the drugs with two-sided p = 6.43 × 10⁻⁷.

Complete pairwise design

Under this design simultaneous application of drugs is conducted at all pairs of the doses. Let S_Ai be the SF of cells in the ith experiment with dose d_Ai of drug A where i = 1, 2, …, n_A. Similarly, let S_Bj be the SF in the jth experiment of drug B with dose d_Bj, where j = 1, 2, …, n_B. It is assumed that in n_An_B experiments (group D), drugs are given in combination (d_Ai, d_Bj) and result in the SF S_Dij. According to Bliss independence we need to compare values ln S_Ai + ln S_Bj with ln S_Dij. As in the previous model, we assume a multiplicative error scheme with x_i = ln S_Ai, y_j = ln S_Bj, and z_ij = ln S_Dij, where

Choice of drug concentration

It is important to chose the right set of doses to conduct experiments when dose-response relationship is estimated or when synergy is being tested by statistical means. A simple rule for estimation of individual dose-response curves has been devised in [29] to optimally choose drug concentrations that induce mortality (or reduce SF) to 0.122 and 0.878, the pivotal points on the logistic response curve.

In choosing optimal drug concentrations to detect synergy it should be remembered that under the null hypothesis of drug interaction the SF is the product of individual SFs. The optimal statistical identification of the product (e.g., when p-value is minimal) is when variance is maximal. If SFs are computed from Bernoulli random variables (dead/alive) for two drugs X and Y, we want to maximize var(XY) where Pr(X = 1) = S_A and Pr(Y = 1) = S_B. The optimal choice is when S_AS_B = 0.5. As an example, we suggest the combinations of drugs that leads to products 0.4, 0.5, and 0.6, which can be achieved using the following individual SFs: 0.6, 0.7, and 0.77, or their combination, depending on the size of the study.

The two-drug copula mortality function

In the analysis above, we were concerned with statistical testing for synergy. But if synergy is detected, how do we predict the mortality probability given two synergistic treatments? The answer relies on the two-drug dose-response function. Several two-drug dose-response relationships, as a part of the response surface methodology, have been suggested [1], but none satisfy the properties formulated below. Most relationships been developed in connection with CI, not Bliss independence, such as [14]; see reviews [2] and [35]. Several authors associated the interaction effect with the product of drug concentrations as customarily used in the linear model framework [36], but justification is lacking. A typical example of an ad hoc dose-response model [37] is a quadratic polynomial and as such can be negative and/or not an increasing function of dose concentration.

The goal of this work is to present a novel rigorous two-drug copula mortality function M(d_A, d_B) that satisfies the following properties:

• Log scale: Model M is expressed on the log scale (i.e., drug concentrations enter the model as ln d_A and ln d_B).

• Singe-drug inheritance in the absence of the other drug, the two-agent copula model collapses to a single-drug model, M_A(d_A), or M(d_A, 0) = M_A(d_A), the mortality function of drug A applied alone. Similarly, we require that M(0, d_B) = M_B(d_B), the mortality function of drug B.

• Independence/synergy parameter: the model depends on parameter ρ, which determines independence, synergy, or antagonism, or more rigorously, (a) when ρ = 0 the model collapses to the Bliss independence model M(d_A, d_B) = M_A(d_A) + M_B(d_B) − M_A(d_A)M_B(d_B), (b) when ρ < 0 we have synergy M(d_A, d_B) > M_A(d_A) + M_B(d_B) − M_A(d_A)M_B(d_B), and (c) when ρ > 0 we have antagonism, i.e. M(d_A, d_B) < M_A(d_A) + M_B(d_B) − M_A(d_A)M_B(d_B).

A two-drug model simplifies statistical determination of independence, synergy, or antagonism: estimate ρ; if the null hypothesis H₀: ρ = 0 is not rejected, we claim Bliss independence; otherwise, synergy or antagonism depending on the sign of ρ. None of the two-agent dose-response models suggested previously satisfy the above properties. The goal of this section is to introduce a family of novel dose-response functions that describe the mortality effect in the presence of two drugs (d_A, d_B) for any given single-drug models M_A(d_A) and M_B(d_B), and use it for statistically determination of synergy.

The basis for our creation is (a) recognition that the mortality function with one or more drugs can be viewed as a cumulative distribution function (cdf) of one or multiple random variables, respectively, routinely used in probability theory, and (b) apply copula to construct a two-agent mortality function using individual mortality functions treated as marginal cdfs in the probability theory [38], [39], [40], [41]. As shown in Appendix (Supporting Information), the two-drug copula mortality function is expressed through a double integral over the bivariate standard normal density with the correlation coefficient ρ as

Correlation coefficient ρ receives a new meaning in our copula mortality function. When ρ < 0, drugs complement each other and therefore the killing effect increases. Conversely, when ρ > 0 the killing effect diminishes when drugs are applied simultaneously. The following theorem defines the limits of mortality when two drugs are completely antagonistic or synergistic as extreme cases when ρ approaches 1 or −1.

Theorem. The properties of the two-drug copula mortality function. (a) Complete antagonism: when two drugs are the same, i.e. M_A = M_B and ρ → 1, then M(d_A, d_B) → M_A(max(d_A, d_B)). (b) Complete synergy: when ρ → −1 then M(d_A, d_B) → M_A(d_A) + M_B(d_B) if M_A(d_A) + M_B(d_B) < 1 and M(d_A, d_B) → 1 otherwise.

The proof is in Appendix (Supporting Information). Complete antagonism indicates that the drugs have completely overlapping effect, e.g. when the drugs affect the same set of cellular receptors, so the mortality is defined by the maximum dose. Complete synergy occurs when the affected cellular receptors do not overlap; then the joint mortality is the sum of single mortalities.

The two-drug copula function is illustrated in Fig 4 by contours of the two-drug function (11) with mortality held at the 50% level. To comply with the Loewe additivity approach, the single-drug mortality functions for this example share the same slope m and therefore give rise to the two-drug function. Specifically, drugs A and B have the same m = 1, 2 but different E C 50 A = 0 . 4 and EC₅₀ = 0.5 and therefore the Loewe independence is depicted as a segment (black), which connects the two EC₅₀s. The fact that Loewe and Bliss independence are not equivalent is well-known and can be seen from this plot (the black and green lines are different). When m = 1 synergy will be overestimated and antagonism will be underestimated using Loewe independence, compared with that derived from our two-drug copula model, because the green contour line (ρ = 0) is below the black Loewe segment. The bottom-left corner depicts synergy and the top-right corner depict antagonism. If drugs complement each other (ρ = −0.5), a smaller dose can lead to the same 50% kill, and on the other hand, increased doses are required to get 50% kill. On the other hand, if drugs have overlapping effect (ρ = 0.5) and are therefore antagonistic. Conclusions change when drugs have a stronger effect on mortality (m = 2): synergy may be claimed according to the copula model while the Loewe approach concludes antagonism (the green contour line is above the black Loewe segment).

Example: Testing lethal effects of insecticides

To illustrate our two-agent copula mortality model, we used classic data on lethal effects of two insecticides, rotenone and pyrethrins, on fruit flies from experiments described by Finney [42] and later analyzed by [15] for determination of synergy/antagonism. The original data contain two series of experiments with the ratio of rotenone:pyrethrins as 1:5 and 1:15 mixtures; we used only the first data series for illustration. Using the CI, the authors declared the 1:5 mixture to be “mildly synergistic.” We apply our statistical model to determine if the synergy is statistically significant.

The estimated probit-based two-drug copula mortality function depicted as a surface in 3D is shown in Fig 5 and the output of the nonlinear least squares algorithm (the nls function in R) is listed in Fig 6.

All parameters are statistically significant. Note that the two drugs have slopes greater than one. They termed by Chou and Talalay [15] as of “high order,” which explains why the 50% contour line points away from the origin as in the right plot in Fig 4. The rho-parameter = −0.6678, and this negative value confirms synergy. The null hypothesis that drugs act independently is rejected with p = 0.0172, and the one-sided p-value for the synergy is 0.0086 = 0.0172/2. Details on the algorithmic implementation and the R code are found in Supporting Information.

Tumor growth experiments in vivo

Before moving to clinical studies, drugs are often tested in animals in vivo. Several authors apply CI for tumor volume at a specific time point [43], [44], [45], [46]. Here we apply Bliss definition to tumor volume data measured over time t in four groups of animals: control (C), group with drug A, group with drug B, and group D with a combination of drugs A and B (all drugs are given in a single dose). Our analysis of drug interaction relies on the following assumptions: (1) tumor volume is reflective of number of cancer cells, (2) at the time when treatments begins (baseline), tumor volume is the same in all four groups, (3) drugs are given over the course of treatment and result in exponential growth (or decay) of tumors with a group-specific rate. The combination of these assumptions leads to group-specific exponential growth/decay of tumors

Below we demonstrate that the hypothesis of drug interaction in tumor growth experiments can be reduced to drug interaction in the presence of a control group as discussed above. Indeed, fix time t and consider the proportion of survived cells in three treatment groups as S_A = A(t)/C(t), S_B = B(t)/C(t), and S_AB = D(t)/C(t). Following Bliss definition, drugs act independently if S_AB = S_AS_B or on the log scale if ln S_AB = ln S_A + ln S_B. However,

This condition is similar to condition (9), but is expressed in terms of rate of tumor growth. We say that there is synergy if β₁ + β₂ − β₃ − β₀ > 0 is statistically supported by the p-value computed from the estimates of the slope and their variances obtained from the output of the R function lme from the library nlme (see details in Supporting Information). In other words, synergy means that rate of tumor growth in the drug combination group yields a slower rate than if drugs acted independently, compared to the control group. Typically, each group of animals is estimated separately using the theory of mixed models [20] to account for differences in individual trajectory and baseline tumor volume. The following example illustrates this approach.

Example: Combination of two drugs (EHT1864 and fulvestrant) for breast tumor treatment

The data [48] and fit are depicted in Fig 7. The R code is listed in Supporting Information. Since exponential function on the log scale turns into a linear function it is advantageous to plot tumor volume data on the log scale but show tickmarks on the y-axis using actual tumor volume. To address inter-animal heterogeneity the data are fitted by a linear function with random effect that reflects the differences of tumor volume at the time of treatment [20]. For example, the rate of breast tumor growth in the control group is estimated as β ^ 0 = 0 . 0432 , which means that tumors grow by 4.32% every week if treated with vehicle. The combination of drugs (EHT1864 plus fulvestrant) induces tumor shrinkage at the rate of 0.97% per week. The above formula for p-value yields 0.006, indicating a statistically significant drug interaction between EHT1864 and fulvestrant.

Survival curves

An advantage of Bliss independence over Loewe additivity, and the related CI, is that drug interaction can be applied to clinical data when the outcome of treatment is survival, e.g. recurrence-free survival, time to recurrence/progression, progression-free survival, or overall survival. We then treat proportions of patients that survive by time t as the probability of an individual to stay alive.

Let treatments A and B, when applied separately, lead to survival curves S_A(t) and S_B(t), respectively, and let S_AB(t) be the survival curve if drugs are applied simultaneously. If drugs were acting independently, as follows from Bliss definition we expect the survival curve to be

If drugs/treatments are synergistic, the observed survival S_AB(t) must be greater than (12), or geometrically, the curve S_AB(t) is above the curve under independence assumption (12). To test for drugs independence, we compare the survival curve (12) with S_AB(t) using standard statistical tests for survival curves, such as the logrank test.

We illustrate this approach using a study of the interaction between nivolumab and ipilimumab for the treatment of patients with metastatic melanoma ([49], [50]). Although the authors hint at possible drug synergy, it was not tested. According to our analysis, the drugs acted independently in two groups of patients, “Intent-to-treat population” and “Patients with PD-L1-negative tumors” (Figs 8 and 9). The survival curve denoted as “drugs independence” and computed as the right-hand side of Eq (12) is expected under the assumption that the two drugs acted independently. This curve is very close to the patient-derived survival curve S_AB denoted “Nivolumab + ipilimumab” as evidence that the two drugs act independently. A formal statistical logrank test produced the p-values 0.58 and 0.37 (using R package survival) suggesting that the hypothesis that these two drugs acted independently cannot be rejected despite the fact that the survival curve with the drug combination is above the single-drug survival curves.

Discussion

Bliss definition of drugs independence has a solid probabilistic justification that relies upon the interpretation of the drug effect as the probability that an individual cell or organism dies. Sometimes, Bliss definition of independence is criticized for possibility of claiming synergy when the same drugs are applied. Although theoretically this is a limitation of Bliss definition, fortunately in practice we never split the same drug into two doses. In the end, the endpoint of dose-response experiments is surviving fraction, not our perception of the biological mechanism of synergy. On the other hand, we argue that Loewe additivity definition may claim independence for biologically synergetic drugs. Indeed, imagine two different cancer drugs that act on different set of receptors but equivalent in terms of their individual effect, i.e. they have the same dose-response relationship. In addition, let us assume that cancer cell die if only two set of receptors are damaged, that is, drugs are biologically synergistic. However, if they have the same dose-response relationship, according to Loewe definition, drugs will act independently.

We offer a rigorous and systematic implementation of Bliss independence applied to common design of preclinical (in vitro and animal) and clinical (e.g., survival) studies. Statistical evaluation and computation of a p- and d-value should be a required component when reporting on drug synergy or antagonism because of considerable variation in response to treatment across organisms and replicates. Moreover, the existing shortage of statistical tests for synergy leads to frequent overstatements and false positive findings. Herein, we demonstrated that Bliss independence can be tested using comprehensive statistical models tailored to major biomedical study designs. Unlike prior reports describing statistically testing for drug independence using the traditional t-test we developed statistical models that comply with the normal distribution required for correct application of many statistical tests and computation of the p-value.

Supporting information

S1 File [pdf]
R Codes and Appendix.

S2 File [csv]
Raw data on Daphnia experiments with replicates.

S3 File [csv]
Raw data on cancer cell experiments in the presence of control group.

S4 File [csv]
Raw data on tumor growth in mice analyzed by a linear mixed model on the log scale.