Additional Approaches

In braidrm: Fitting Combined Action with the BRAID Response Surface Model

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(braidrm)

Introduction

This will come as no surprise to anyone who has performed even a cursory review of the literature, but BRAID it not the only approach available for combination analysis. The debate over the best methods and standards for quantifying drug combinations goes back nearly 100 years, and has yet to be truly resolved. And while we (unsurprisingly) feel that the BRAID method is the best method for comprehending combined action overall, we can certainly acknowledge that many circumstances call for a different approach. So we have tried, with the braidrm package, to include a set of tools for quantifying drug combinations using many of the most popular approaches available. These are certainly not the only implementations of these metrics and models, but it is our hope that they are still sufficiently flexible to be of use.

Response Surface Methods

Though it is the model that we have found to be the most effective as an analytical tool, BRAID is not the only response surface method in the combination analysis space. The braidrm package includes implementations of two other parametric response surface models, one which predates BRAID by about 25 years, and one that followed it.

The URSA Model

The universal response surface approach (URSA) was proposed by Greco, Park, and Rustum [-@Greco1990] 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)

ufit2 <- fitUrsaModel(measure ~ concA+concB, synergisticExample)
coef(ufit2)

ufit3 <- fitUrsaModel(measure ~ concA+concB, antagonisticExample)
coef(ufit3)

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. [-@Wooten2021] 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)

mfit2 <- fitMusycModel(measure ~ concA+concB, synergisticExample)
coef(mfit2)

mfit3 <- fitMusycModel(measure ~ concA+concB, antagonisticExample)
coef(mfit3)

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.

The Deviation Methods

In its most basic form, pharmacological "synergy" is defined as an observed effect in combination that is greater than what would be expected based on the effect of either drug in isolation. It makes sense then that many approaches to combination analysis would tackle the question directly, and calculate how a measured surface deviates from the expected, putatively non-interacting surface. We refer to these methods as the "deviation" methods, and they quantify synergy and antagonism (either at particular dose pairs or in the combination overall) by estimating such deviations. The primary differences between them are what model they use as the "expected" non-interacting surface. The braidrm package includes functions for evaluating four such methods: Bliss deviation, HSA deviation, Loewe deviation, and ZIP $\delta$. These are described briefly here:

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 [@Bliss1939], 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 [@Loewe1926]. 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 [@Yadav2015]. 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)))
mean(deviationSurface(concs1,act1,method="HSA",increasing=TRUE))
mean(deviationSurface(concs1,act1,method="Loewe"))
mean(deviationSurface(concs1,act1,method="ZIP",range=c(0,1)))

concs2 <- cbind(synergisticExample$concA,additiveExample$concB)
act2 <- synergisticExample$measure
mean(deviationSurface(concs2,act2,method="Bliss",range=c(0,1)))
mean(deviationSurface(concs2,act2,method="HSA",increasing=TRUE))
mean(deviationSurface(concs2,act2,method="Loewe"))
mean(deviationSurface(concs2,act2,method="ZIP",range=c(0,1)))

concs3 <- cbind(antagonisticExample$concA,additiveExample$concB)
act3 <- antagonisticExample$measure
mean(deviationSurface(concs3,act3,method="Bliss",range=c(0,1)))
mean(deviationSurface(concs3,act3,method="HSA",increasing=TRUE))
mean(deviationSurface(concs3,act3,method="Loewe"))
mean(deviationSurface(concs3,act3,method="ZIP",range=c(0,1)))

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

The Combination Index

Of course no discussion of combination analysis would be complete without the combination index [@Chou1984]. One of the most widely cited and widely used methods for combination analysis, the combination index (and equivalent methods such as the sum of FICs [@Berenbaum1989], observed over expected, and interaction index [@Berenbaum1978]) quantifies the degree of interaction in a combination by the extent to which constant ratio combinations are potentiated relative to Loewe additivity. It can be a bit difficult to wrap one's head around, but briefly, the combination index for a given combination, for a particular effect level and a particular dose ratio, is the "dose" of that constant ratio combination required to produce that particular effect, divided by the dose that would be required were the combination purely additive. As a result, additive surfaces give combination indices near 1, synergistic surfaces, being more potent, produce combination indices below 1, and antagonistic surfaces produce combination indices above 1.

The nature of the combination index makes specifying its inputs rather unwieldy, so braidrm defaults to an input format similar to that of the deviation methods. The result is a set of combination index values for all dose ratios at which a sufficient number of dose pairs were measured:

estimateCombinationIndices(concs1,act1,level=c(0.5))
estimateCombinationIndices(concs2,act2,level=c(0.5))
estimateCombinationIndices(concs3,act3,level=c(0.5))

We can see that for many dose ratios (particularly those near an equal mixture of drug A and drug B), the additive surface exhibits a combination index near 1, the synergistic surface exhibits a combination index around 0.5, and the antagonistic surface gives values up around 2. But these values are not consistent across dose ratios, and unfortunately, the combination index method does not expect them to be. The combination index, therefore, is not a robust metric describing the overall behavior of the surface, making it difficult to compare between combinations of very different dosages, Hill slopes, and potencies.

braidrm documentation built on Sept. 30, 2024, 9:40 a.m.