Additional Approaches

The URSA Model

The universal response surface approach (URSA) was proposed by Greco, Park, and Rustum (1990) as an early attempt to numerically fit a quantitative model of synergy and antagonism to measured data. Like BRAID, it used pharmacological additivity as its basis, but unlike BRAID it hewed more closely to the precise concept of Loewe additivity. They started with the observation that in a Loewe additive surface the following relation always holds:

\[ 1 = \frac{D_A}{{ID}_{X,A}} + \frac{D_B}{{ID}_{X,B}} \]

where \({ID}_{X,A}\) and \({{ID}_{X,B}}\) are the doses of drug A and drug B that produce the same effect in isolation that the drug combination of \(D_A\) and \(D_B\) produce together. When the individual combinations are assumed to follow a standard Hill model, this can be written more explicilty (if more cumbersomely) as:

\[ 1 = \frac{D_A}{{ID}_{M,A}\left(\frac{E-E_0}{E_f-E}\right)^{1/n_a}} + \frac{D_B}{{ID}_{M,B}\left(\frac{E-E_0}{E_f-E}\right)^{1/n_b}} \]

The URSA model generalized this relationship to allow for synergistic and antagonistic interactions by adding an interaction term that is very nearly a classic product term:

\[ 1 = \frac{D_A}{{ID}_{M,A}\left(\frac{E-E_0}{E_f-E}\right)^{1/n_a}} + \frac{D_B}{{ID}_{M,B}\left(\frac{E-E_0}{E_f-E}\right)^{1/n_b}} + \alpha\frac{D_AD_B}{{ID}_{M,A}{ID}_{M,B}\left(\frac{E-E_0}{E_f-E}\right)^{1/{2n_a}}\left(\frac{E-E_0}{E_f-E}\right)^{1/{2n_b}}} \]

When \(alpha\) is positive, less of both doses is needed to sum to a full 1, so lower combined doses produce the same effect, resulting in synergy. When \(alpha\) is negative, the doses must be increased to reach the full value of 1, resulting in antagonism.

The URSA model can be fit to combination data using the fitUrsaModel() function:

ufit1 <- fitUrsaModel(measure ~ concA+concB, additiveExample)
coef(ufit1)
#>         IDMA         IDMB           na           nb        alpha           E0 
#>  0.847395180  0.990519924  2.378979964  2.690630700 -0.153006078  0.001130186 
#>           Ef 
#>  1.002739329

ufit2 <- fitUrsaModel(measure ~ concA+concB, synergisticExample)
coef(ufit2)
#>        IDMA        IDMB          na          nb       alpha          E0 
#>  0.99993303  1.01984863  2.21919208  2.16306644  8.25168963 -0.03841135 
#>          Ef 
#>  0.99986881

ufit3 <- fitUrsaModel(measure ~ concA+concB, antagonisticExample)
coef(ufit3)
#>         IDMA         IDMB           na           nb        alpha           E0 
#>  1.450038016  1.547877012  3.208442983  2.654243308 -0.292355825  0.004973926 
#>           Ef 
#>  0.992736539

Encouraginly, the URSA model fits the additive surface with an \(alpha\) near 0, the synergistic surface with a clearly positive \(alpha\) and the antagonistic surface with a negative \(alpha\). However, the modle does suffer some significiant limitations.

Because it is an extension of Loewe additivity, it is hindered by the same constraints, most notably that the maximal effects of both drugs must be identical for surface to be definable at all. Further, because URSA, like Loewe additivity, is defined by an insoluble implicit equation, evaluation of the surface is necessarily a numerical approximation, which is both slow and difficult to extend to other tasks that require derivatives. Finally, while we have included antagonistic URSA surfaces (with \(alpha\) ranging down to \(-1\)), the product term in the URSA equation actually produces undefined values at high doses for any negative \(alpha\) value, so predictions of the model at high doses should be treated with caution.

The MuSyC Model

Proposed much more recently, the Multidimensional Synergy of Combinations (or MuSyC) model of Wooten et al. (2021) is a beautifully versatile response surface model that can be fit with up to twelve free parameters in its most flexible form. Though it is described as a synthesis of all competing methods, it is best understood as an extension of the principles of Bliss independence (albeit a very significant generalization of Bliss surfaces). The model treats a combined therapy as a four-state system, in which the measured target (cells, organisms, bound proteins, etc.) transitions between the four states as a function of the two drug concentrations. The four states are:

• \(U\): Targets unaffected by either drug
• \(A_1\) Targets affected by drug 1, but not drug 2
• \(A_2\): Targets affected by drug 2, but not drug 1
• \(A_{12}\): Targets affected by both drugs

Each of these states produces a distinct level of the measured effect, and the overal measured effect is a weight average of those four governed by the relative occupancy of the four states. The simplest form would simply assume that any affected state would produce a 100% effect, so the measured effect would simply be the total proportion of targest that are affected in any way; but the strength of the MuSyC models is that each of these values can be numerically varied to produce a much wider range of observed behaviors, inlcuding partial and observed behaviors.

Interaction in the combination is modeled in several ways. First, because the effect level measured in state \(A_{12}\) can differ from the maximal effects resulting from either drug, the combinatoin can have a higher (or even lower) maximal effect than the individual agents. Wooten et al. refer to this as “efficacy synergy”, and while it can be included in the full BRAID model, it is much less a component of traditional BRAID analysis. Second, because the pharmacological transitions induced by each drug can differ depending on whether they are already affected by the other drug, the effect of one drug can potentiate the effect of the other; a phenomenon Wooten et al. refer to as “potency synergy”. (It is also possible for the effect of one drug to alter the Hill slope of the other, but this is much more difficult to analyze, and is usually not included.)

Similar to URSA, the MuSyC model can be fit to combination data using the fitMusycModel() function:

mfit1 <- fitMusycModel(measure ~ concA+concB, additiveExample)
coef(mfit1)
#>         IDMA         IDMB           na           nb      alpha12      alpha21 
#>  0.800139368  0.949387585  1.943812580  2.261862297  2.060378090  1.621066302 
#>           E0          EfA          EfB           Ef 
#> -0.008523318  1.027721216  1.037863582  0.972305644

mfit2 <- fitMusycModel(measure ~ concA+concB, synergisticExample)
coef(mfit2)
#>        IDMA        IDMB          na          nb     alpha12     alpha21 
#>  1.19668859  0.96655322  1.52782988  2.26339095 18.64303977 18.85244067 
#>          E0         EfA         EfB          Ef 
#> -0.05573616  1.11679988  0.99075362  0.99841452

mfit3 <- fitMusycModel(measure ~ concA+concB, antagonisticExample)
coef(mfit3)
#>        IDMA        IDMB          na          nb     alpha12     alpha21 
#> 1.265611199 1.240076136 2.965970304 2.977843271 0.019376190 0.048681256 
#>          E0         EfA         EfB          Ef 
#> 0.001749911 0.997186076 0.957373529 1.106652844

Encouragingly, the MuSyC fits posit no significant efficacy synergy, as the underlying surfaces all share the same maximal effect. The potency synergy parameters \(\alpha_{12}\) and \(\alpha_{21}\), however, show more distinction. These parameters represent the degree to which the effect of one drug increase the potency of the other; with \(\alpha_{12}\) representing the degree to which being affected by drug 1 increases the potency of drug2, and \(\alpha_{21}\) representing the reverse. (Unlike BRAID, MuSyC can represent synergies operating in both directions.) These values are well above the baseline 1 for the synergistic surface, and well below 1 for the antagonistic surface. The model does posit some slight synergy even for the additive surface, but this is unsurprising as additive surfaces with higher Hill slopes are often marked as synergistic by Bliss-based and Bliss-inspired methods.

Though it is quite flexible and very powerful, we have still found BRAID to be a more reliable method overall. The high-dimensionality of the interaction in MuSyC (with three or even five different parameters governing interaction) makes an intuitive understanding of the interactions difficult, and comparison across surfaces nearly impossible. While we agree that potency synergy and efficacy synergy are distinct, and deserve to be treated as such, we have found efficacy synergy to be more rare, and in most cases very difficult to distinguish from more basic potentiation. Additionally, the basis of the model around Bliss independence, however implicit, is out of sync with our overall experience that additivity is generally a better, ess biased zero point. Nevertheless, we frequently make use of the MuSyC model in our work, particularly when quantifying efficacy synergy directly is desired, and hope that users of the braidrm package will make use of MuSyC as well.

Bliss Deviation

Bliss deviation, perhaps the most commonly used deviation approach, quantifies the interaction in a response surface by modeling the effect of the two drugs as probabilistically independent events. The measured effect reflects the proportion of targets that are affected or unaffected by either drug, and the probability that a target will be unaffected by two doses in combination is equal to the product of the probabilities that it would be unaffected by each dose individually. This principle, Bliss independence (Bliss 1939), is expressed mathematically by the following equation:

\[ A_{12} = A_1 + A_2 - A_1A_2 \] where \(A_{12}\) is the proportion of targets affected by the combination, and \(A_1\) and \(A_2\) are the proportions affected by either dose individually. Bliss independence is extremely popular, due the intuitive nature of combining independent events and the mathematical simplicity of the non-interacting model.

Highest Single Agent

Even simpler than Bliss, the highest single agent (HSA) model of non-interaction simply assumes that the effect of two combined doses will simply be the maximum of the effects of the two doses individually:

\[ A_{12} = \max\left(A_1,A_2\right) \]

Despite its extreme simplicity, the model is still fairly widely used, often in parallel with Bliss independence as it will always necessarily be lower than the Bliss independent prediction.

Loewe Additive

Though Loewe-based methods are more common as response surface models, deviations from a Loewe additive response surface can still be measured directly (Loewe and Muischnek 1926). The only constraint is that the maximal effects of both drugs must be identical so that a Loewe additive surface is well defined for all dose pairs. Mathematically, Loewe additivity is defined as:

\[ 1 = \frac{D_A}{{ID}_{X,A}} + \frac{D_B}{{ID}_{X,B}} \]

where \({ID}_{X,A}\) and \({{ID}_{X,B}}\) are the doses of drug A and drug B that produce the same effect in isolation that the drug combination of \(D_A\) and \(D_B\) produce together.

Zero-Interaction Potency

A more recent entry to the deviation method landscape, the zero-interaction potency (or ZIP) model of Yadav et al. produces a more robust measure of response surface deviation than the Bliss deviation method on which it is based (Yadav et al. 2015). The ZIP reference surface is a Bliss independent surface, but before the combined effects are calculated, the dose-response behaviors of both individual drugs are fit with a Hill model to produce stabler, smoother effects.
Furthermore, before the reference surface is subtracted from them, individual combined measurements are replaced with smoothed values resulting from fitting each set of measurements with a shared dose of either drug with a partial dose response curve. The result, though mathematically identical to Bliss deviation at its center, is a smoother, more robust deviation surface.

Deviation Functions in braidrm

All four deviation methods can be accessed using the deviationSurface() function in the braidrm package. The function produces a full deviation surface, but can be summed or averaged to produce an estimated deviation metric:

concs1 <- cbind(additiveExample$concA,additiveExample$concB)
act1 <- additiveExample$measure
mean(deviationSurface(concs1,act1,method="Bliss",range=c(0,1)))
#> [1] 0.03101272

mean(deviationSurface(concs1,act1,method="HSA",increasing=TRUE))
#> [1] 0.033698

mean(deviationSurface(concs1,act1,method="Loewe"))
#> [1] -0.00859843

mean(deviationSurface(concs1,act1,method="ZIP",range=c(0,1)))
#> [1] 0.02199841

concs2 <- cbind(synergisticExample$concA,additiveExample$concB)
act2 <- synergisticExample$measure
mean(deviationSurface(concs2,act2,method="Bliss",range=c(0,1)))
#> [1] 0.1348243

mean(deviationSurface(concs2,act2,method="HSA",increasing=TRUE))
#> [1] 0.1400639

mean(deviationSurface(concs2,act2,method="Loewe"))
#> [1] 0.09550931

mean(deviationSurface(concs2,act2,method="ZIP",range=c(0,1)))
#> [1] 0.1250206

concs3 <- cbind(antagonisticExample$concA,additiveExample$concB)
act3 <- antagonisticExample$measure
mean(deviationSurface(concs3,act3,method="Bliss",range=c(0,1)))
#> [1] -0.06080502

mean(deviationSurface(concs3,act3,method="HSA",increasing=TRUE))
#> [1] -0.04623165

mean(deviationSurface(concs3,act3,method="Loewe"))
#> [1] -0.09061148

mean(deviationSurface(concs3,act3,method="ZIP",range=c(0,1)))
#> [1] -0.03461665

Details about the additional parameters for each of the different deviation methods can be found in the documentation for deviationSurface().