BRAID: A Unifying Paradigm for the Analysis of Combined Drug Action

With combination therapies becoming increasingly vital to understanding and combatting disease, a reliable method for analyzing combined dose response is essential. The importance of combination studies both in basic and translational research necessitates a method that can be applied to a wide range of experimental and analytical conditions. However, despite increasing demand, no such unified method has materialized. Here we introduce the Bivariate Response to Additive Interacting Doses (BRAID) model, a response surface model that combines the simplicity and intuitiveness needed for basic interaction classifications with the versatility and depth needed to analyze a combined response in the context of pharmacological and toxicological constraints. We evaluate the model in a series of simulated combination experiments, a public combination dataset, and several experiments on Ewing’s Sarcoma. The resulting interaction classifications are more consistent than those produced by traditional index methods, and show a strong relationship between compound mechanisms and nature of interaction. Furthermore, analysis of fitted response surfaces in the context of pharmacological constraints yields a more concrete prediction of combination efficacy that better agrees with in vivo evaluations.


Derivation of the BRAID Model
To describe the combined effect of two doses of compounds A and B, the BRAID model assumes that the individual effects of both compounds can be modeled by a Hill or log-logistic equation 1 : A Loewe additive combination is one that everywhere satisfies the equation where X is the combined effect of D A and D B , and ID X,A and ID X,B are the concentrations of drug A and B alone required to produce the effect X. Thus, any effect produced by any combination of drugs A and B must be also produced at some concentration by both drugs A and B in isolation; one consequence of this is that both drugs must produce the same range of effects, and have the same maximal effects. We begin with the equation for a Loewe additive surface when Hill slopes and maximal effects are equal: In order for this equation to give an additive surface when E f = E f,A = E f,B and n a = n b , it must be that n = n a = n b . Thus, more generally, n should be a symmetric function of n a and n b that simplifies to either value when they are equal. Given the multiplicative nature of exponents, the geometric mean seems most appropriate.

If we assume (without loss of arbitrariness) that E f,A is the larger effect of E f,A and E f,B
, and that E f = E f,A , then we get the following equation, the BRAID model of additive combined action: Note that this equation assumes that That is, the effects of drug A and B are in the same direction, and the maximal effect of drug A is at least as large as that of drug B. It is easy to show that this equation simplifies to the behavior of each drug alone when the other drug is not present; that is, .
Substituting E f and n back into the above equation and expressing ̃ and ̃ in terms of E A (D A ) and E B (D B ) allows us to write the following implicit form of the additive BRAID equation: Given that this equation describes the BRAID model of additivity, a reasonable method for introducing interaction would be to add a multiplicative interaction term to the equation. However, if this interaction term involves the products of the two terms on the right hand side of the equation, than any negative coefficient on the interaction term will create areas of dose-pairspace that produce undefined effects. We therefore introduce an interaction term that is the geometric mean of the two terms on the right hand side, inspired by Greco, Park and Rustum 2 : Interestingly, an alternative model of interaction can be achieved by altering the exponent in the denominator of the additive BRAID model, rather than the base. By adjusting the exponent √ with a multiplicative factor δ, one gets what we call the delta-BRAID model: In this model, δ > 1 implies synergy, δ < 1 implies antagonism, and δ = 1 implies additivity. The synergistic surfaces produced by the delta-BRAID model are very similar to those produced by the standard BRAID model, but the antagonistic surfaces differ drastically, behaving in a similar fashion to the antagonistic surfaces produced by the model of Greco, Park, and Rustum 2 . Though this model is in some ways more elegant than the standard BRAID model, we find its results to be unsatisfactory; in particular, even extreme δ-antagonism cannot create surfaces in which the potency of one drug is reduced by the presence of another. Though some have argued that such response surfaces are invalid 2 , we have observed such effects in several experiments, and have thus found the delta-BRAID model to be insufficient.
An additional complication arises if one expects the maximal effect of a combination to differ from the maximal effect of either drug alone. Though we have not observed this in practice, it is theoretically conceivable, and thus worth considering. Fortunately, it is not difficult to adjust the standard BRAID model to account for such a circumstance. To incorporate the widest range of possible response surfaces, including both kappa-and delta-interaction, and a differing combined maximal effect, we developed the 10 parameter extended-BRAID or eBRAID model: In practice, we have found the eight parameter standard BRAID model to be entirely sufficient; nevertheless, all functions in the braidrm R package support the full 10-parameter eBRAID model if desired.

Derivation of Bayesian Concentration Correction
If we assume that the true underlying relationship between concentration and effect for a single compound is modeled by a function f(•), and that the measured value of an effect Ê is normally distributed around the true underlying effect E, then the probability of a measured effect Ê given an actual dose ̂ is where G(•) is a standard normal distribution and ϵ is the standard deviation of the experimental noise. Furthermore, if we assume that the actual concentration ̂ underlying a measurement is log-normally distributed around the intended concentration D, then the probability of a particular actual concentration is where η is the standard log-deviation of the concentration error. According to Bayes' rule,

Rationale for Fitting Logarithmic Transform of Cell Survival
Recent approaches to pharmacodynamic modeling have extended use of the Hill equation from a model of static equilibria of reactions to a component of dynamic models in which the Hill equation describes the rate at which a particular process or event (e.g. cell growth/killing) occurs [3][4][5] . For example, the growth or death of a particular cell population might be modeled as where N is the population size at time t, λ is a baseline growth rate, and ε is a maximal kill rate for a particular drug. Following this approach, we can model the growth or death of a cell population in the presence of a drug combination as It is easy to see that the right-hand side of this equation is another instance of the BRAID model, this time with E 0 = logN 0 + λt e . Hence, fitting the BRAID model to the logarithmic transform of the measured cell population is appropriate both from a statistical and a mechanistic viewpoint.

Potentiation and IAE Calculations
To determine the potentiation of a compound A by the presence of a dose D B of a second compound B at a particular effect level E, one must know the dose of A alone required to produce the effect E and the dose of A required to produce the effect E in combination with the dose D B of B that produces the same effect. In the case of the BRAID surface model, this requires a partial inversion of the BRAID equation. Fortunately, this task, though cumbersome, is relatively straightforward. On the assumption that E 0 < E < E f,A and The process for determining a dose of compound B which produces an effect E in combination with a dose D A of compound A proceeds similarly. It is worth noting that when κ is negative, some inputs produce two solutions; the equation above gives the larger, more conservative solution. All of these calculations are performed by the R function invertBRAIDrsm in the braidrm package. The potentiation can then be described by taking the ratio of the amount of compound A needed to produce the effect E alone ( −1 ( ,~,0)) with amount needed in the presence of dose D B ( −1 ( ,~, )).
The traditional therapeutic index consists of the ratio of the maximum achievable dose subject to certain constraints (pharmacokinetics, toxicity, metabolism, etc.) to the minimum dose required to produce a desired effect 6 . This is equivalent to calculating the ratio of the range of achievable doses to the range of doses too low to produce the desired effect. Performing the analogous calculation for dose-pairs gives the index of achievable efficacy, or IAE. The IAE for a given combination at an effect level E is given by: where H(•) is the Heaviside step function equal to 1 for values greater than 0, and 0 for values less than 0, and AC is the space of achievable dose pairs given pharmacological constraints such as solubility and toxicity. The inclusion of E f,A -E simply ensures that the input to the step function has the appropriate sign. Taking the square root keeps the IAE at the same scale as the therapeutic index, so that, for example, the IAE of a drug combined with itself will be the same as the therapeutic index of that drug alone.
In the simplest case AC can be considered a rectangle of dose-pairs bounded by the maximum achievable concentration of either drug; this is assumed in the IAE calculations reported in our paper. The form of AC, however, may be much more complicated: if both drugs are identical, for example, AC is bounded above by a diagonal of slope -1 representing dose pairs whose concentrations sum to the maximum achievable concentration of the drug. AC may even be determined by a threshold placed on another BRAID surface representing overlapping toxicity. In its current form, the function calculateIAE in the braidReports package supports only the use of a rectangle of achievable doses if the drugs are different, or the lower triangle of achievable doses if the two drugs are the same.
In practice, it is not possible to calculate this quantity analytically; but it can be estimated to any desired level of accuracy by calculating the predicted BRAID effect at N points in the space of achievable dose pairs. We use a 300-by-300 grid of dose pairs within the bounding box of achievable dose pairs to generate estimates of the IAE and corresponding 95% confidence intervals.