In maomlab/BayesPharma: Tools for Bayesian Analysis of Non-Linear Pharmacology Models

MuSyC Epistasis Model {#sec-MuSyC}

When two different treatments are made in an assay, their combined effect may be stronger or weaker than what would be expected with an additive model, the treatments are said to be epistatic For sigmoidal dose-response models, one treatment may effect of the other in two different ways; by either shifting the maximal response (efficacy) or by shifting the dose needed to cause the response (potency). A range of statistical models have been proposed that capture different aspects of synergy, notably Bliss independence[@Bliss1956-hf] and Loewe additivity[@Loewe1926-sr] models can be used to test for significant efficacy or potency interactions, respectively. The SynergyFinder R package[@Ianevski2022-cb] and the synergy python package[@Wooten2021-cr] can be used to visualize treatment interactions, compute a range of synergy scores, and test if the interactions are significant.

Recently Meyer et al.[@Meyer2019-zr,@Wooten2021-lg] derived an integrated functional synergistic sigmoidal dose-response called the Multi-dimensional Synergy of Combinations (MuSyC) method, which has the Loewe and Bliss models as special cases. They implemented a Bayesian model-fitting strategy in Matlab, and a maximum likelihood model-fitting into the synergy python package. To make the model more accessible to the pharmacology community, in this section, we briefly review the MuSyC functional form, describe a Bayesian implementation in Stan/BRMS, and illustrate using the model to re-analyze how drugs and voltage may interact to modulate the current through a potassium channel.

MuSyC Functional Form: The functional form for the MuSyC model gives an equation for the response $\color{brown}{E_d}$ at doses of $\color{teal}{d_1}$ and $\color{teal}{d_2}$ of the two treatments and has $9$ free parameters $\theta = \left({\color{purple}{C_1}}, {\color{purple}{C_2}}, {\color{brown}{E_0}}, {\color{brown}{E_1}}, {\color{brown}{E_2}}, {\color{brown}{E_3}}, {\color{purple}{h_1}}, {\color{purple}{h_2}}, {\color{purple}{\alpha}}\right)$

\begin{align} {\color{brown}{E_d}} = \mbox{MuSyC}({\color{teal}{d_1}}, {\color{teal}{d_2}}; \theta) &= \frac{ {\color{purple}{C_1}}^{\color{purple}{h_1}}{\color{purple}{C_2}}^{\color{purple}{h_2}}{\color{brown}{E_0}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}{\color{purple}{C_2}}^{\color{purple}{h_2}}{\color{brown}{E_1}} + {\color{purple}{C_1}}^{\color{purple}{h_1}}{\color{teal}{d_2}}^{\color{purple}{h_2}}{\color{brown}{E_2}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}{\color{teal}{d_2}}^{\color{purple}{h_2}}{\color{purple}{\alpha}} {\color{brown}{E_3}} }{ {\color{purple}{C_1}}^{\color{purple}{h_1}}{\color{purple}{C_2}}^{\color{purple}{h_2}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}{\color{purple}{C_2}}^{\color{purple}{h_2}} + {\color{purple}{C_1}}^{\color{purple}{h_1}}{\color{teal}{d_2}}^{\color{purple}{h_2}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}{\color{teal}{d_2}}^{\color{purple}{h_2}}{\color{purple}{\alpha}}} \end{align} The parameters $\color{brown}{E_0}$, $\color{brown}{E_3}$, give response values for the extreme values for the doses: $\mbox{MuSyC}({\color{teal}{0}}, {\color{teal}{0}}; \theta) = \color{brown}{E_0}$ and $\mbox{MuSyC}({\color{teal}{\infty}}, {\color{teal}{\infty}}; \theta) = \color{brown}{E_3}$. Setting one of the doses to zero e.g., ${\color{teal}{d_2}} = 0$, the MuSyC functional form reduces to a sigmoid function of the other, where $\mbox{MuSyC}({\color{teal}{d_1}}, {\color{teal}{0}}; \theta) = Sigmoid({\color{teal}{d_1}}; \phi)$, where the half maximal activity is $\mbox{AC}_{50} = \color{purple}{C_1}$ and the slope at the half maximal activity is ${\color{purple}{\mbox{hill}}} = {\color{purple}{h_1}}$ and if $h_1 > 0$, then ${\color{brown}{\mbox{top}}} = {\color{brown}{E_1}}$ and ${\color{brown}{\mbox{bottom}}} = {\color{brown}{E_0}}$, otherwise the assignment is reversed. See Appendix XXX for a derivation.

To interpret these parameters if we set $\color{teal}{d_2}=0$, then \begin{align} \color{brown}{E_d} &= \frac{ {\color{purple}{C_1}}^{\color{purple}{h_1}}{\color{brown}{E_0}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}{\color{brown}{E_1}} }{ {\color{purple}{C_1}}^{\color{purple}{h_1}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}} \end{align} which is the Hill equation, which we modeled above \ref{sec:hill}. If we then additionally set $\color{teal}{d_1}=0$ then $\color{brown}{E_d}=\color{brown}{E_0}$, in the limit as ${\color{teal}{d_1}}\rightarrow \infty$ then ${\color{brown}{E_d}}\rightarrow {\color{brown}{E_1}}$, and if ${\color{teal}{d_1}}=\color{purple}{C_1}$ then ${\color{brown}{E_d}} = ({\color{brown}{E_0}} + {\color{brown}{E_2}})/2$, which is the half maximal response (either the $\color{brown}{\mbox{IC}{50}}$ if treatment $1$ is an inhibitor or $\color{brown}{\mbox{EC}{50}}$ if treatment $1$ is agonist). The slope at ${\color{teal}{d_1}}={\color{purple}{C_1}}$ is \begin{align} \frac{\mathrm{d}\;\color{brown}{E_d}}{\mathrm{d}\color{teal}{d_1}} &= {\color{purple}{C_1}}^{v}{\color{brown}{E_0}} \frac{\mathrm{d}}{\mathrm{d}\color{teal}{d_1}} \frac{1}{{\color{purple}{C_1}}^{\color{purple}{h_1}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}} + {\color{brown}{E_1}} \frac{\mathrm{d}}{\mathrm{d}\color{teal}{d_1}} \frac{{\color{teal}{d_1}}^{h_1}}{{\color{purple}{C_1}}^{\color{purple}{h_1}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}}\ &= {\color{purple}{C_1}}^{h_1}{\color{brown}{E_0}} \frac{ h_1{\color{teal}{d_1}}^{{\color{purple}{h_1}}-1}}{\left({\color{purple}{C_1}}^{\color{purple}{h_1}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}\right)^2} + {\color{brown}{E_1}} \frac{{\color{purple}{C_1}}^{\color{purple}{h_1}}h_1{\color{teal}{d_1}}^{{\color{purple}{h_1}}-1}}{\left({\color{purple}{C_1}}^{\color{purple}{h_1}} + {\color{teal}{d_1}}^{\color{purple}{h_1}}\right)^2}\ &= ({\color{brown}{E_0}} + {\color{brown}{E_1}}) \end{align}

The evaluation of the functional form for ${\color{brown}{E_d}}$ is numerically unstable due to the exponentiation. To transform using the $\texttt{log-sum-exp}$ trick, let

\begin{align} \texttt{numerator_parts} = [\ &{\color{purple}{h_1}}\log({\color{purple}{C_1}}) + {\color{purple}{h_2}}\log({\color{purple}{C_2}}) + \log({\color{brown}{E_0}}),\ &{\color{purple}{h_1}}\log({\color{teal}{d_1}}) + {\color{purple}{h_2}}\log({\color{purple}{C_2}}) + \log({\color{brown}{E_1}}),\ &{\color{purple}{h_1}}\log({\color{purple}{C_1}}) + {\color{purple}{h_2}}\log({\color{teal}{d_2}}) + \log({\color{brown}{E_2}}),\ &{\color{purple}{h_1}}\log({\color{teal}{d_1}}) + {\color{purple}{h_2}}\log({\color{teal}{d_2}}) + \log({\color{brown}{E_3}}) + \log({\color{purple}{\alpha}}) ]\ \texttt{denominator_parts} = [\ &{\color{purple}{h_1}}\log({\color{purple}{C_1}}) + {\color{purple}{h_2}}\log({\color{purple}{C_2}}),\ &{\color{purple}{h_1}}\log({\color{teal}{d_1}}) + {\color{purple}{h_2}}\log({\color{purple}{C_2}}),\ &{\color{purple}{h_1}}\log({\color{purple}{C_1}}) + {\color{purple}{h_2}}\log({\color{teal}{d_2}}),\ &{\color{purple}{h_1}}\log({\color{teal}{d_1}}) + {\color{purple}{h_2}}\log({\color{teal}{d_2}})]\ \end{align} Then for a vector $x = [x_1, x_2, \dots, x_n]$, let $\texttt{log_sum_exp}(x) = \mbox{log}\left(\mbox{exp}(x_1) + \mbox{exp}(x_2) + \dots + \mbox{exp}(x_n)\right)$. Then

$$ {\color{brown}{E_d}} = \mbox{exp}!\left(\texttt{log-sum-exp}(\texttt{numerator_parts}) - \texttt{log-sum-exp}(\texttt{denominator_parts})\right). $$ To implement the $\mbox{MuSyC}$ model in Stan we use the following parameterization.

\scriptsize

# the brms MuSyC formula with given covariates
synergy_formula <- MuSyC_formula(predictors = covariates)

# will generate a formula like this
synergy_formula_alt <- brms::brmsformula(
  # The Stan MuSyC function is defined in BayesPharma::MuSyC_stanvar()
  response ~ MuSyC(
    logd1 - logd1scale,
    logd2 - logd2scale,
    logE0,
    logC1, logE1, h1,
    logC2, logE2, h2,
    logE3, logalpha),
  nl = TRUE) +
  # The free parameters are regressed against the given covariates
  brms::lf(logE0              ~ covariates) +
  brms::lf(logC1 + logE1 + h1 ~ covariates) +
  brms::lf(logC2 + logE2 + h2 ~ covariates) +
  brms::lf(logE3 + logalpha   ~ covariates)

\normalsize

Note that if the logd1scale and logd1scale values are not provided in the the data, when the model is run, they are automatically computed as the mean value of the doses and is used to make the model easier to fit.

MuSyC model BayesPharma

#' Drug Synergy #' MuSyC Drug Synergy model #' #' Assume that the response metric decreases with more effective drugs #' Let E3 be the effect at the maximum concentration of both drugs #' #' #' Special cases: #' * dose additive model: alpha1 = alpha2 = 0 #' * loewe: h1 = h2 = 1 #' * CI: E0 = 1, E1 = E2 = E3 = 0 #' the drug effect is equated with percent inhibition #' * bliss drug independence model: #' E0 = 1, E1 = E2 = E3 = 0, alpha1 = alpha2 = 1 #' @param d1 Dose of drug 1 #' @param d2 Dose of drug 2 #' #' @param E0 effect with no drug treatment #' #' # params for drug 1 by it self #' @param s1 drug 1 hill slope #' @param C1 drug 1 EC50 #' @param E1 drug 1 maximum effect #' #' # params for drug 2 by it self #' @param s2 drug 2 hill slope #' @param C2 drug 2 EC50 #' @param E2 drug 2 maximum effect #' #' @param beta synergistic efficacy #' percent increase in a drug combination's effect #' beyond the most efficacious single drug. #' #' beta > 0 => synergistic efficacy #' the effect at the maximum concentration of both drugs (E3) exceeds the #' maximum effect of either drug alone (E1 or E2) #' #' beta \< 0 => antagonistic efficacy #' at least one or both drugs are more efficacious as #' single agents than in combination #' #' @param alpha1 synergistic potency #' how the effective dose of drug 1 #' is altered by the presence of drug 2 #' @param alpha2 synergistic potency #' how the effective dose of drug 2 #' is altered by the presence of drug 1 #' #' alpha > 1 => synergistic potency #' the EC50 decreases because of the addition of the other drug, #' corresponding to an increase in potency #' #' 0 \<= alpha \< 1 => antagonistic potency #' the EC50 of the drug increases as a result of the other drug, #' corresponding to a decrease in potency #' #' alpha1 == alpha2 if detailed balance #' @export generate_MuSyC_effects \<- function( d1, d2, E0, s1, C1, E1, s2, C2, E2, alpha, E3) { h1 \<- MuSyC_si_to_hi(s1, C1, E0, E1) h2 \<- MuSyC_si_to_hi(s2, C2, E0, E2) numerator \<- C1\^h1 * C2\^h2 * E0 + d1\^h1 * C2\^h2 * E1 + C1\^h1 * d2\^h2 * E2 + d1\^h1 * d2\^h2 * E3 * alpha denominator \<- C1\^h1 * C2\^h2 + d1\^h1 * C2\^h2 + C1\^h1 * d2\^h2 + d1\^h1 * d2\^h2 * alpha numerator / denominator }

#' Create a formula for the MuSyC synergy model #' #' @description setup a defaulMuSyC synergy model formula to predict #' the E0, C1, E1, s1, C2, E2, s2, log10alpha, and E3alpha #' parameters. #' #' @param predictors Additional formula objects to specify predictors of #' non-linear parameters. i.e. what perturbations/experimental differences #' should be modeled separately? (Default: 1) should a random effect be taken #' into consideration? i.e. cell number, plate number, etc. #' @return brmsformula #' #' @examples #'\dontrun{

' Consider observations made using 4 different drugs and the column header

' containing the labels for the 4 different drugs is predictors.

'

' MuSyC_formula(predictors = 0 + predictors)

'

' Consider that the cell_ID was recorded and the noise from the cell_ID is to

' be accounted for.

'

' MuSyC_formula(predictors = 0 + predictors + (1|cell_ID))

'} #' #' @export MuSyC_formula \<- function( predictors = 1, ...) {

predictor_eq <- rlang::new_formula(
  lhs = quote(E0 + C1 + E1 + s1 + C2 + E2 + s2 + log10alpha + E3alpha),
  rhs = rlang::enexpr(predictors))

brms::brmsformula(
  response ~ (C1^h1 * C2^h2 * E0 +
      d1^h1 * C2^h2 * E1 +
      C1^h1 * d2^h2 * E2 +
      d1^h1 * d2^h2 * E3alpha
    ) / (
      C1^h1 * C2^h2 +
      d1^h1 * C2^h2 +
      C1^h1 * d2^h2 +
      d1^h1 * d2^h2 * 10^log10alpha),
  brms::nlf(d1 ~ dose1 / d1_scale_factor),
  brms::nlf(d2 ~ dose2 / d2_scale_factor),
  brms::nlf(h1 ~ s1 * (4 * C1) / (E0 + E1)),
  brms::nlf(h2 ~ s2 * (4 * C2) / (E0 + E2)),
  predictors_eq,
  nl = TRUE,
  ...)

}

#' Fit the MuSyC synergy model by dose #' #' @param data data.frame of experimental data #' with columns: dose1, dose2, n_positive, count, [] #' @param group_vars quosures list #' dplyr::vars(...) columns to over when fitting synergy model #' @param C1_prior prior distribution for Ed when d1=d1_IC50, d2=0 #' @param C2_prior prior distribution for Ed when d1=0, d2=d2_IC50 #' @param s1_prior prior distribution for d(Ed)/d(d1) when d1=d1_IC50, d2=0 #' @param s2_prior prior distribution for d(Ed)/d(d2) when d1=0, d2=d2_IC50 #' @param log10alpha_prior prior distribution for alpha synergy parameter #' @param E0_prior prior distribution for Ed when d1=0, d2=0 #' @param E1_prior prior distribution for Ed when d1=Inf, d2=0 #' @param E2_prior prior distribution for Ed when d1=0, d2=Inf #' @param E3_alpha_prior prior distribution for Ed scaled by alpha when d1=Inf, #' d2=Inf #' @param C1_init initial sampling distribution for the C1 parameter #' @param C2_init initial sampling distribution for the C2 parameter #' @param s1_init initial sampling distribution for the s1 parameter #' @param s2_init initial sampling distribution for the s2 parameter #' @param log10alpha_init initial sampling distribution for the alpha parameter #' @param E0_init initial sampling distribution for the E0 parameter #' @param E1_init initial sampling distribution for the E1 parameter #' @param E2_init initial sampling distribution for the E2 parameter #' @param E3_alpha_init initial sampling distribution for the E3 parameter #' @param combine combine the grouped models into a single brms model #' @param verbose verbose output #' #' @param iter number of stan NUTS sampling steps #' @param cores number of cores used for sampling #' @param stan_model_args stan model arguments #' @param control stan control arguments #' #' The #' #' bernoulli_inf(n_positive / count) = #' Ed \~ MuSyC(d1, d2, C_params, E_params, s_params, alpha) #' #' To improve numeric stability, the d1 and d2 and C1 and C2 variables #' are scaled to improve numeric stability: #' #' d1 = dose1/max(dose1) #' d2 = dose2/max(dose2) #' drug1_IC50 = C1 * max(dose1) #' drug2_IC50 = C2 * max(dose2) #' #' Functional form: #' Ed \~ ( #' C1\^h1 * C2\^h2 * E0 + #' d1\^h1 * C2\^h2 * E1 + #' C1\^h1 * d2\^h2 * E2 + #' d1\^h1 * d2\^h2 * E3 * alpha #' ) / ( #' C1\^h1 * C2\^h2 + #' d1\^h1 * C2\^h2 + #' C1\^h1 * d2\^h2 + #' d1\^h1 * d2\^h2 * alpha #' ) #' #' #' #' #' #' ############################################## #' # Proof of the definitions of the parameters # #' ############################################## #' #' Claim: When d1=0 and d2=0 then Ed = E0 #' Ed = (C1\^h1 * C2\^h2 * E0) / (C1\^h1 * C2\^h2) #' = E0 #' #' Claim: When d1=0 and d2 -> Inf then Ed = E2 #' Ed = (C2\^h2 * E0 + d2\^h2 * E2) / (C2\^h2 + d2\^h2) #' = (d2\^h2 * E2) / (d2\^h2) #' = E2 #' #; Claim: When d1=0 and d2=C2 then Ed = (E0 + E2) / 2 #' When d1>0 and d2 -> Inf then Ed #' Ed = (C1\^h1 * C2\^h2 * E0 + C1\^h1 * C2\^h2 * E2) / #' (C1\^h1 * C2\^h2 + C1\^h1 * C2\^h2) #' = (E0 + E2) / 2 #' #'@export MuSyC_model \<- function( data, group_vars = vars(compound), formula = MuSyC_formula(), prior = MuSyC_prior(), init = MuSyC_init(), combine = FALSE, verbose = FALSE, iter = 8000, cores = 4, stan_model_args = list(verbose = FALSE), control = list( adapt_delta = .99, max_treedepth = 12), model_evaluation_criteria = c("loo", "bayes_R2"), ...) {

if (is.data.frame(well_scores)) { grouped_data \<- well_scores \|> dplyr::group_by(!!!group_vars) \|> dplyr::mutate( d1_scale_factor = max(dose1), d2_scale_factor = max(dose2)) \|> tidyr::nest() \|> dplyr::ungroup() }

if (verbose) { cat("Fitting MuSyC model\n") }

model \<- brms::brm_multiple( formula = formula, data = grouped_data\$data, family = binomial("identity"), prior = prior, init = init, # stanvars = c( # brms::stanvar( # scode = " real d1_scale_factor = max(dose1));", # block ="tdata", # position = "end"), # brms::stanvar( # scode = " real d2_scale_factor = max(dose2));", # block ="tdata", # position = "end"), # brms::stanvar( # scode = " real drug1_IC50 = b_C1 * d1_scale_factor);", # block ="genquant", # position = "end"), # brms::stanvar( # scode = " real drug2_IC50 = b_C2 * d2_scale_factor;", # block ="genquant", # position = "end")), combine = FALSE, data2 = NULL, iter = iter, cores = cores, stan_model_args = stan_model_args, control = control, ...)

if (!is.null(model_evaluation_criteria)) { # evaluate fits model \<- model \|> purrr::imap(function(model, i) { group_index \<- grouped_data[i, ] \|> dplyr::select(-data) group_index_label \<- paste0( names(group_index), ":", group_index, collapse = ",") cat("Evaluating model fit for", group_index_label, "...\n", sep = "") model \<- model \|> brms::add_criterion( criterion = model_evaluation_criteria, model_name = paste0("MuSyC:", group_index_label), reloo = TRUE) model }) } grouped_data \|> dplyr::mutate( model = model) }

maomlab/BayesPharma documentation built on Aug. 24, 2024, 8:45 a.m.