Derive: MuSyC Model -- Synergy Analysis
derive_MuSyC_model.Rmd## Error in library(tidyverse): there is no package called 'tidyverse'MuSyC synergy model
When two different treatments are combined they may interact. For
end-point assays, if the response is stronger or weaker than what would
be expected with an additive model, the treatments are said to be
epistatic. For sigmoidal dose-response models, however, the analysis may
be more complicated. One drug may not only may shift the maximal
response (efficacy) of the other, but it may also shift the effective
dose and shape of the response (potency). Historically a range of models
have been proposed that capture different aspects of synergy, for
example the Bliss independence(Bliss 1956)
and Loewe additivity(Loewe 1926) are
null-models for no synergistic efficacy or potency, respectively. The
SynergyFinder R package(Ianevski,
Giri, and Aittokallio 2022) and the synergy python
package(Wooten and Albert 2021) can be
used to visualize treatment interactions, compute a range of synergy
scores, and test if the interactions are significant.
Recently Meyer et al.(Meyer et al. 2019,
Wooten2021–lg) derived an integrated functional synergistic
sigmoidal dose-response, 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 \(\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}\):
\[\begin{align} \color{brown}{E_d} &= \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}\]
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. To transform using the \(\mbox{log\_sum\_exp}\) trick, let
\[\begin{align*} \mbox{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}}) ]\\ \mbox{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 \[ E_d = \mbox{exp}\!\left(\mbox{log\_sum\_exp}(\mbox{numerator\_parts}) - \mbox{log\_sum\_exp}(\mbox{denominator\_parts})\right). \]
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?) #‘ #’ #’
(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,
[
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”) }
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, “…”, 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) }