Skip to contents

Derive: MuSyC Model -- Synergy Analysis

derive_MuSyC_model.Rmd

MuSyC Epistasis Model

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(Bliss 1956) and Loewe additivity(Loewe 1926) models can be used to test for significant efficacy or potency interactions, 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.Wooten et al. (2021) 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 99 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 h1>0h_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 . 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 11 is an inhibitor or $\color{brown}{\mbox{EC}_{50}}$ if treatment 11 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=[x1,x2,…,xn]x = [x_1, x_2, \dots, x_n], let πš•πš˜πš_πšœπšžπš–_πšŽπš‘πš™(x)=log(exp(x1)+exp(x2)+…+exp(xn))\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 MuSyC\mbox{MuSyC} model in Stan we use the following parameterization.

# 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)

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?) #β€˜ #’ #’ (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”) }

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) }

Bliss, C I. 1956. β€œThe Calculation of Microbial Assays.” Bacteriol. Rev. 20 (4): 243–58. https://doi.org/10.1128/br.20.4.243-258.1956.
Ianevski, Aleksandr, Anil K Giri, and Tero Aittokallio. 2022. β€œSynergyFinder 3.0: An Interactive Analysis and Consensus Interpretation of Multi-Drug Synergies Across Multiple Samples.” Nucleic Acids Research. https://doi.org/10.1093/nar/gkac382.
Loewe, S. 1926. β€œEffect of Combinations: Mathematical Basis of Problem.” Arch. Exp. Pathol. Pharmakol. 114: 313–26.
Meyer, Christian T, David J Wooten, B Bishal Paudel, Joshua Bauer, Keisha N Hardeman, David Westover, Christine M Lovly, Leonard A Harris, Darren R Tyson, and Vito Quaranta. 2019. β€œQuantifying Drug Combination Synergy Along Potency and Efficacy Axes.” Cell Syst 8 (2): 97–108.e16. https://doi.org/10.1016/j.cels.2019.01.003.
Wooten, David J, and RΓ©ka Albert. 2021. β€œSynergy: A Python Library for Calculating, Analyzing and Visualizing Drug Combination Synergy.” Bioinformatics 37 (10): 1473–74. https://doi.org/10.1093/bioinformatics/btaa826.
Wooten, David J, Christian T Meyer, Alexander L R Lubbock, Vito Quaranta, and Carlos F Lopez. 2021. β€œMuSyC Is a Consensus Framework That Unifies Multi-Drug Synergy Metrics for Combinatorial Drug Discovery.” Nat. Commun. 12 (1): 4607. https://doi.org/10.1038/s41467-021-24789-z.