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

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

library(braidrm)

The BRAID Model

The bivariate response to additive interacting doses (or BRAID) model is a parametric response surface model of the effect of combined doses of two agents. A full specification of the BRAID equation can be found in vignette("derivation"), but briefly, the BRAID model is specified by 9 parameters: a baseline minimal effect observed when no agent is present, three dose response parameters for each agent, an overall maximal effect parameter, and an interaction paramter $\kappa$ indicating the presence of a synergistic or antagonistic interaction.

In the functions of the braidrm package, BRAID surfaces are specified by a numeric parameter vector. For a full, length-9 parameter vector, the values are the following

IDMA: The dose of median effect of the first agent
IDMB: The dose of median effect of the second agent
na: The Hill slope (or sigmoidicity) of the first agent
nb: The Hill slope (or sigmoidicity) of the second agent
kappa: The BRAID interaction parameter, which is between -2 and 0 for antagonistic surfaces, 0 for BRAID additive surfaces, and greater than 0 for synergistic surfaces
E0: The minimal effect observed when neither agent is present
EfA: The maximal effect of the first agent, theoretically observed at infinite concentration
EfB: The maximal effect of the second agent, theoretically observed at infinite concentration
Ef: The overall maximal effect of the combination

The order of these parameters is meant to mirror the order of single-agent dose-response parameters in the basicdrm package: "potency" parameters, "shape" parameters (including interaction), minimal effect, maximal effects.

In many cases, the response surface being modeled or studied will only need a subset of these parameters. Nearly all combination analyses assume (implicitly or explicitly) that the overall maximal effect of a combination (the parameter $E_f$) will be equal to the "larger" of the two individual maximal effects. (We place "larger" in quotes here is it refers to the effect value that is further from the minimal effect, but not necssarily the numerically greater value.) In some cases, it is even assumed that all maximal effects are simply equal to each other. For simplicity, many braidrm functions support shortened BRAID parameter vectors that carry these assumptions. If a BRAID parameter vector with eight values is passed to a braidrm function, it is assumed that the ninth parameter Ef has been left out, and the it should be equal to the "larger" of the two individual maximal effects. If a BRAID paramter vector with seven values is passed, it is assumed that parameters EfA and EfB have been leftout, and that they are both equal to the seventh value (assumed to be Ef).

Evaluate BRAID Surfaces

Creating a BRAID parameter vector is as simple as creating a numeric vector:

basicParameter <- c(1, 1, 3, 3, 0, 0, 100)

This vector specifies, in order, that:

The dose of median effect of drug A is 1 (the units are not specified, but for this vignette we'll assume they are micromolar)
The dose of median effect of drug B is 1
The Hill slope of drug A is 3
The Hill slope of drug B is 3
The interaction parameter $\kappa$ is 0 (additivity)
The minimal effect is 0
The overall maximal effect is 100

Because this vector is length seven, parameters EfA and EfB are implied and assumed to be equal to Ef. We could specify exactly the same surface with a full-length parameter vector:

fullParameter <- c(1, 1, 3, 3, 0, 0, 100, 100, 100)

Evaluating a BRAID surface requires the concentration or concentrations of the first drug, the concentration or concentrations of the second, and the BRAID parameter vector. For example:

evalBraidModel(1, 0, basicParameter)
evalBraidModel(0, 1, basicParameter)
evalBraidModel(1, 1, basicParameter)

The first two outputs demonstrate that, as expected, when either drug is present at 1 micromolar, we observe an effect of 50, halfway between the minimal effect and the maximal effect. The third output is noticeably higher, as when both drugs are present at 1 micromolar, the effect of the individual doses is compounded. We can however produce the same effect as 1 micromolar of either drug alone by halving their doses:

evalBraidModel(0.5, 0.5, basicParameter)

This is because the parameter vector represents an additive combination of two drugs with identical pharmacological properties. (Note that this does not work for all response surfaces, or even all BRAID additive surfaces. It is only in the case of identical Hill slopes and maximal effects that BRAID additivity aligns perfectly with the more classical Loewe additivity.) [@Loewe1926]

We can also see that using the full, length-9 parameter vector produces exactly the same results:

evalBraidModel(1, 0, fullParameter)
evalBraidModel(0, 1, fullParameter)
evalBraidModel(0.5, 0.5, fullParameter)

Note what happens, however, when we introduce an interaction to the response surface (in this case, as $\kappa$ is positive, a synergistic interaction):

synergyParameter <- c(1, 1, 3, 3, 1.5, 0, 100)

evalBraidModel(1, 0, synergyParameter)
evalBraidModel(0, 1, synergyParameter)
evalBraidModel(0.5, 0.5, synergyParameter)

While the effect of the individual drugs is unchanged, their combined effect is increased. This is the classic pharmacological definition of synergy: an effect in combination which is larger than what would be expected. What precisely is "expected" for any given pair of drugs is one of the primary debates of combination analysis; BRAID additivity is only one answer, albeit the one around which we have built the BRAID system.

Fitting BRAID Surfaces

set.seed(20240804)

Of course, the most common usage of the BRAID model is fitting it to experimental data to summarize and quantify that data. The main workhorse function for this task is braidrm, which produces a fit object of class braidrm. We can see it in action with the example data-set additiveExample:

head(additiveExample)
additiveFit <- braidrm(measure ~ concA + concB, additiveExample)
summary(additiveFit)

The data frame additiveExample is a synthetically generated set of response surface measurements, and contains four columns: concA, containing the dose of drug A for each measurement; concB, containing the dose of drug B for each measurement; truth containing the ground truth response generated by an additive response surface with effect ranging from 0 to 1; and measure, a noisy measurement sampled from a normal distribution centered on the ground truth with a standard deviation of 0.07. By default, braidrm() fits an eight-paramter BRAID surface (one which assumes the overall maximal effect is equal to the larger of the two individual maximal effects) to the given data with a moderate prior on the interaction parameter $\kappa$. (See vignette("modelsAndConstraints") for more details on specifying the BRAID model to be fit, and vignette("bayesianKappa") for a more in-depth explanation of the Bayesian stabilization of $\kappa$). By default braidrm() also estimates bootstrapped confidence intervals on all fit parameters, as can be seen the printed summary; note in particular that 0 lies within the confidence interval on $\kappa$, indicating no statistically significant deviation from BRAID additivity. (Estimating confidence intervals can take several seconds; so if you are running many fits and do not need the confidence intervals, you can forgo them by setting the parameter getCIs to FALSE.)

We get a very different result, however, when we fit the example data-set synergisticExample, which has the same format as additiveExample but was generated with a synergistic parameter vector with a $\kappa$ of 2:

synergisticFit <- braidrm(measure ~ concA + concB, synergisticExample)
summary(synergisticFit)

Though most of the parameter estimates are very similar (indeed the generating parameter vectors are identical except for $\kappa$), the confidence interval on $\kappa$ lies well above zero, centered quite close the true value of 2.

Another useful fitting function is findBestBraid(), which uses the Bayesian or Akaike information criterion to select among several candidate models (again, see vignette("modelsAndConstraints") for more details) [@Schwarz1978, @Akaike1974]. This allows the user to stabilize underdetermined values to reasonable defaults and reduces the frequency of wildly implausible over-fitting. We have run the function on antagonisticExample, which, as the name implies, contains a synthetically generated set of response surface measurements taken on an antagonistic response surface with a true $\kappa$ of $-1$. The usage of findBestBraid() is very similar to that of braidrm():

bestFit <- findBestBraid(measure ~ concA + concB, antagonisticExample,
                         defaults = c(0,1))
summary(bestFit)

The defaults parameter here indicates that, absent any other evidence, we expect the minimal effect to 0 and the maximal effects to be 1, so simpler models that assume these fixed values and explain the data equally well should be preferred. This time the confidence interval on $\kappa$ lies well below zero and comfortably encloses the true value of $-1$. No confidence intervals on either minimal or maximal effects are included, indicating that the most parsimonious model was one which fixed E0 at 0 and the three maximal effects at the default 1.