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

# inspired by
# https://www.jumpingrivers.com/blog/knitr-default-options-settings-hooks/
knitr::opts_chunk$set(
  echo = FALSE,
  cache = TRUE,
  cache.path = "cache/apply_MuSyC_KCNQ/",
  fig.path = "apply_MuSyC_KCNQ_files/",
  fig.retina = 2, # Control using dpi
  fig.width = 6,  # generated images
  fig.height = 5, # generated images
  fig.pos = "t",  # pdf mode
  fig.align = "center",
  dpi = if (knitr::is_latex_output()) 72 else 300, 
  out.width = "100%")
options(width = 300)

#| echo=FALSE
suppressWarnings(suppressMessages(library(tidyverse)))
library(BayesPharma)

Re-analysis Figure 1 of (Li, et al., 2021) Molecular Basis for Ligand Activation of the Human KCNQ2 Channel

KCNQ2 is a voltage gated potassium channel important for re-establishing neuronal homeostasis following an action potential. Dis-regulation of channels in the KCNQ family can cause epilepsy, tinnitus, and depression. While there has been sustained interest in developing small-molecule based therapeutics targeting KCNQ2, the only FDA approved drug that targets KCNQ2, retigabine, was given a black-box warning due to unwanted side-effects, even though it was effective at treating epilepsy. Excitingly there has been recent progress in using Cryo-EM to structurally characterize ion channels including KCNQ2, which promise to support structure-based drug design.

Recently, [@Li2021-ef] used Cryo-EM to characterize the structure of KCNQ2 in complex with retigabine, which interacts in the membrane on the outside of the pore, and ztz240, which interacts with the voltage sensor domain. To relate the structure and function, [@Li2021-ef] measured conductance as a function of voltage (G-V curves) in the presence of varying doses of retigabine and ztz240 using whole cell patch-clamp electrophysiology in Chinese hamster ovary (CHO)-K1 cells overexpressing KCNQ2. Both compounds are thought to be agonists, that is with increasing concentration of drug, less negative voltages are required to open the channel.

In this case study, the aim is to re-analyze the effects of retigabine and ztz240 using the data collected in [@Li2021-ef] and presented in panels B and D of figure 1. A key idea is to recognize that voltage and drug treatment can be thought of independent perturbations and the drug effect can be framed as characterizing the interaction between these perturbations.

To begin, we will load and plot the data,

#| echo=FALSE,
#| message=FALSE

read_range <- function(
    sheet,
    range,
    treatment,
    dose_uM,
    is_control){
  readxl::read_excel(
    path = here::here(
      "inst", "extdata", "conductance_KCNQ", "Li_CellRes_2021_fig1.xlsx"),
    sheet = {{sheet}},
    range = range,
    col_names = FALSE,
    .name_repair = ~ paste0("replica_", 1:length(.))) |>
    dplyr::mutate(
      treatment = {{treatment}},
      dose_uM = {{dose_uM}},
      is_control = {{is_control}},
      voltage = seq(-90, 60, by = 10),
      .before = 1) |>
    tidyr::pivot_longer(
      cols = c(-treatment, -dose_uM, -is_control, -voltage),
      names_to = "replica",
      values_to = "conductance")
}

data <- dplyr::bind_rows(
  read_range(
    sheet = "fig1b", range = "E3:I18",
    treatment = "ztz240", dose_uM = 10, is_control = TRUE),
  read_range(
    sheet = "fig1b", range = "L3:P18",
    treatment = "ztz240", dose_uM = 10, is_control = FALSE),
  read_range(
    sheet = "fig1b", range = "E22:I37",
    treatment = "ztz240", dose_uM = 5, is_control = TRUE),
  read_range(
    sheet = "fig1b", range = "L22:P37",
    treatment = "ztz240", dose_uM = 5, is_control = FALSE),  
  read_range(
    sheet = "fig1b", range = "E40:I55",
    treatment = "ztz240", dose_uM = 3, is_control = TRUE),
  read_range(
    sheet = "fig1b", range = "L40:P55",
    treatment = "ztz240", dose_uM = 3, is_control = FALSE), 
  read_range(
    sheet = "fig1b", range = "E58:J73",
    treatment = "ztz240", dose_uM = 1, is_control = TRUE),
  read_range(
    sheet = "fig1b", range = "L58:Q73",
    treatment = "ztz240", dose_uM = 1, is_control = FALSE),   
  read_range(
    sheet = "fig1b", range = "E77:I92",
    treatment = "ztz240", dose_uM = 0.1, is_control = TRUE),
  read_range(
    sheet = "fig1b", range = "L77:P92",
    treatment = "ztz240", dose_uM = 0.1, is_control = FALSE),

  read_range(
    sheet = "fig 1d", range = "D2:I17",
    treatment = "retigabine", dose_uM = 30, is_control = TRUE),
  read_range(
    sheet = "fig 1d", range = "M2:R17",
    treatment = "retigabine", dose_uM = 30, is_control = FALSE),
  read_range(
    sheet = "fig 1d", range = "D20:I35",
    treatment = "retigabine", dose_uM = 10, is_control = TRUE),
  read_range(
    sheet = "fig 1d", range = "M20:R35",
    treatment = "retigabine", dose_uM = 10, is_control = FALSE),  
  read_range(
    sheet = "fig 1d", range = "D39:H54",
    treatment = "retigabine", dose_uM = 5, is_control = TRUE),
  read_range(
    sheet = "fig 1d", range = "M39:Q54",
    treatment = "retigabine", dose_uM = 5, is_control = FALSE), 
  read_range(
    sheet = "fig 1d", range = "D57:H72",
    treatment = "retigabine", dose_uM = 3, is_control = TRUE),
  read_range(
    sheet = "fig 1d", range = "M57:Q72",
    treatment = "retigabine", dose_uM = 3, is_control = FALSE),   
  read_range(
    sheet = "fig 1d", range = "D76:H91",
    treatment = "retigabine", dose_uM = 1, is_control = TRUE),
  read_range(
    sheet = "fig 1d", range = "M76:Q91",
    treatment = "retigabine", dose_uM = 1, is_control = FALSE),
  read_range(
    sheet = "fig 1d", range = "D95:H110",
    treatment = "retigabine", dose_uM = 0.3, is_control = TRUE),
  read_range(
    sheet = "fig 1d", range = "M95:Q110",
    treatment = "retigabine", dose_uM = 0.3, is_control = FALSE))
data

\normalsize

#| echo=FALSE,
#| fig.width=4,
#| fig.heigh=4,
#| fig.cap="Reproduction of Figure 1b and 1d from (Li, et al., 2021)",
#| dependson=c("load-data")
plot_data <- data |>
  dplyr::mutate(dose_label = ifelse(is_control, 0, dose_uM)) |>
  dplyr::arrange(dose_label) |>
  dplyr::mutate(
    dose_label = ifelse(dose_label == 0, "Control", paste0(dose_uM, " uM")) |>
      forcats::fct_inorder()) |>
  dplyr::arrange(desc(treatment)) |>
  dplyr::mutate(
    treatment = treatment |> forcats::fct_inorder())

data_normalize <- plot_data |>
  dplyr::filter(is_control) |>
  dplyr::group_by(treatment, dose_uM, voltage) |>
  dplyr::summarize(
    conductance_mean = mean(conductance),
    .groups = "drop") |>
  dplyr::group_by(treatment, dose_uM) |>
  dplyr::summarize(
    conductance_max = max(conductance_mean),
    .groups = "drop")

plot_data <- plot_data |>
  dplyr::left_join(
    data_normalize,
    by = c("treatment", "dose_uM")) |>
  dplyr::mutate(
    conductance_normalized = conductance / conductance_max)

plot_data_replica_summary <- plot_data |>
  dplyr::group_by(treatment, dose_label, voltage) |>
  dplyr::summarize(
    conductance_mean = mean(conductance_normalized),
    conductance_se = sd(conductance_normalized) / sqrt(dplyr::n()),
    .groups = "drop")

ggplot2::ggplot() +
  ggplot2::theme_bw() +
  ggplot2::geom_line(
    data = plot_data_replica_summary,
    mapping = ggplot2::aes(
      x = voltage,
      y = conductance_mean,
      color = dose_label,
      group = dose_label)) +
  ggplot2::geom_errorbar(
    data = plot_data_replica_summary,
    mapping = ggplot2::aes(
      x = voltage,
      ymin = conductance_mean - conductance_se,
      ymax = conductance_mean + conductance_se,
      color = dose_label,
      group = dose_label),
    width = 2) +
  ggplot2::geom_point(
    data = plot_data_replica_summary,
    mapping = ggplot2::aes(
      x = voltage,
      y = conductance_mean,
      color = dose_label,
      group = dose_label)) +
  ggplot2::facet_wrap(
    facets = dplyr::vars(treatment), ncol = 1) +
  ggplot2::scale_color_discrete("Drug Dose") +
  ggplot2::scale_x_continuous("Voltage") +
  ggplot2::scale_y_continuous("G/Gmax")

Then, for each treatment and dose, we will fit a sigmoid curve, which has an comparable functional form to the what they call the Boltzmann equation for G/Gmax as a function of voltage. To make the fitting more stable, we will transform the treatment and response scales.

#| dependson=c("load-data")
model_data <- data |>
  dplyr::filter(
    is_control == FALSE) |>
  dplyr::transmute(
    conductance,
    voltage = voltage/100,
    treatment = treatment,
    doselabel = paste(dose_uM, "uM") |> as.factor(),
    logdose = log10(dose_uM) - 6)

#| echo=TRUE,
#| results='hide',
#| message=FALSE,
#| dependson=c("prepare-model-data")

model_conductance <- BayesPharma::sigmoid_model(
  data = model_data,
  formula = BayesPharma::sigmoid_agonist_formula(
    treatment_variable = "voltage",
    treatment_units = "mV",
    response_variable = "conductance",
    response_units = "% Gmax",
    predictors = 0 + treatment:doselabel),
  prior = BayesPharma::sigmoid_agonist_prior(
    ec50 = brms::prior(normal(-0.2, 1), nlpar = "ec50"),
    hill = brms::prior(normal(3, 2), nlpar = "hill", lb = 0)),
  init = BayesPharma::sigmoid_agonist_init(
    ec50 = \() stats::runif(1, min = -7, max = -5),
    hill = \() stats::runif(1, min = 0.8, max = 1.2),
    bottom = \() stats::runif(1, min = -.1, max = 0.1),
    top = \() stats::runif(1, min = 0.8, max = 1.2)),
  cores = 4)

#| echo=FALSE,
#| dependson=c("fit-model-conductance")
model_conductance

#| echo=TRUE,
#| message=FALSE,
#| dependson=c("model-conductance")

model_conductance |>
  BayesPharma::plot_prior_posterior_densities(refresh = 0)

#| echo=TRUE,
#| message=FALSE,
#| dependson=c("model-conductance")

model_conductance |>
  BayesPharma::plot_posterior_draws()

Fit MuSyC Model

#| echo=FALSE,
#| message=FALSE,
#| dependson=c("prepare-model-data")

model_MuSyC <- BayesPharma::MuSyC_model(
  data = model_data,
  formula = BayesPharma::MuSyC_formula(
    treatment_1_variable = "voltage",
    treatment_2_variable = "logdose",
    response_variable = "conductance"),
  prior = BayesPharma::MuSyC_prior(
    logE0 = log(0),
    logE1 = log(1),
    logC1 = brms::prior(normal(-0.2, 1), nlpar = "logC1"),
    h1 = brms::prior(normal(1, 1), nlpar = "h1"),
    logE2 = log(0),
    logC2 = brms::prior(normal(-6, 2), nlpar = "logC2"),
    h2 = brms::prior(normal(-1, 1), nlpar = "h2"),    
    logE3 = log(1),
    logalpha = brms::prior(normal(0, 1), nlpar = "logalpha")),
  init = BayesPharma::MuSyC_init(
    logE0 = NULL,
    logE1 = NULL,
    logC1 = \() stats::runif(n = 1, min = -0.4, max = 0),
    h1 = \() stats::runif(n = 1, min = 0.8, max = 1.2),
    logE2 = NULL,
    logC2 = \() stats::runif(n = 1, min = -7, max = -5),
    h2 = \() stats::runif(n = 1, min = -1.2, max = -0.8),
    logE3 = NULL,
    logalpha = \() stats::runif(n = 1, min = -0.2, max = 0.2)),
  control = list(
      adapt_delta = .99,
      max_treedepth = 12),
  stanvars = BayesPharma::MuSyC_stanvar(),
  cores = 4)

Join the conductance model with the dose information to model how the voltage dependence depends on the drug dose

#| echo=FALSE,
#| message=FALSE,
#| dependson=c("prepare-model-data", "model-conductance")

data_conductance_by_voltage <- model_conductance |>
  posterior::summarize_draws() |>
  tidyr::separate_wider_delim(
    cols = variable,
    delim = "_",
    names = c("variable_type", "variable", "predictors_label"),
    too_few = "align_start") |>
  dplyr::filter(variable_type == "b") |>
  dplyr::inner_join(
    model_data |>
      dplyr::distinct(treatment, logdose, doselabel) |>
      dplyr::mutate(
        predictors_label = paste0(
          "treatment", treatment, ":",
          "doselabel", doselabel |> stringr::str_replace(" ", ""))),
    by = "predictors_label")
data_conductance_by_voltage

#| echo=FALSE,
#| message=FALSE,
#| dependson=c("extract-voltage-by-dose")
plot <- ggplot2::ggplot(
  data = data_conductance_by_voltage) +
  ggplot2::theme_bw() +
  ggplot2::geom_ribbon(
    mapping = ggplot2::aes(
      x = logdose,
      ymin = q5,
      ymax = q95,
      fill = treatment),
    alpha = .2) +
  ggplot2::geom_line(
    mapping = ggplot2::aes(
      x = logdose,
      y = mean,
      color = treatment),
    linewidth = 1.8) +
  ggplot2::facet_wrap(
    facets = dplyr::vars(variable),
    scales = "free_y") +
  ggplot2::theme(legend.position = "bottom") +
  ggplot2::scale_color_viridis_d("Treatment", begin=.2, end=0.8) +
  ggplot2::scale_fill_viridis_d("Treatment", begin=0.2, end=0.8) +
  ggplot2::scale_x_continuous("Log[Molar]") +
  ggplot2::scale_y_continuous("Parameter Value")
plot

#| echo=FALSE,
#| message=FALSE,
#| dependson=c("prepare-model-data")


formula_bilevel <- brms::brmsformula(
  conductance ~ sigmoid(
    sigmoid(Eec50, Ehill, Etop, Ebottom, logdose),
    sigmoid(Hec50, Hhill, Htop, Hbottom, logdose),
    vtop,
    vbottom,
    voltage),
  vtop + vbottom ~ 1,
  Eec50 + Ehill + Etop + Ebottom ~ 0 + treatment,
  Hec50 + Hhill + Htop + Hbottom ~ 0 + treatment,
  nl = TRUE)

prior_bilevel <- c(
  brms::prior(prior = normal(-6, 2.5), nlpar = "Eec50"),
  brms::prior(prior = normal(-1.5, 1), nlpar = "Ehill", ub = 0.01),
  brms::prior(prior = normal(-.1, 0.05), nlpar = "Etop"),
  brms::prior(prior = normal(-.4, 0.05), nlpar = "Ebottom"),

  brms::prior(prior = normal(-6, 2.5), nlpar = "Hec50"),
  brms::prior(prior = normal(1, 1), nlpar = "Hhill", lb = 0.01),
  brms::prior(prior = normal(5, 2), nlpar = "Htop"),
  brms::prior(prior = normal(3, 0.1), nlpar = "Hbottom"),    

  brms::prior(prior = normal(1, 0.3), nlpar = "vtop"),
  brms::prior(prior = normal(0, 0.3), nlpar = "vbottom"))

init_bilevel <- list(
  b_Eec50 = \() stats::runif(n = 1, min = -7, max = -5),
  b_Ehill = \() stats::runif(n = 1, min = -1.6, max = -1.4),
  b_Etop = \() stats::runif(n = 1, min = -.12, max = -.08),
  b_Ebottom = \() stats::runif(n = 1, min = -.42, max = -0.38),

  b_Hec50 = \() stats::runif(n = 1, min = -7, max = -5),
  b_Hhill = \() stats::runif(n = 1, min = 0.8, max = 1.2),
  b_Htop = \() stats::runif(n = 1, min = 4, max = 6),
  b_Hbottom = \() stats::runif(n = 1, min = 2.9, max = 3.1),

  b_vtop = \() stats::runif(n = 1, min = 0.7, max = 1.3),
  b_vbottom = \() stats::runif(n = 1, min = -0.3, max = .3)) |>
  `class<-`(c("bpinit")) |>
  BayesPharma::eval_init(
    sdata = brms::make_standata(
      formula = formula_bilevel,
      data = model_data,
      prior = prior_bilevel),
    chains = 4)

model_bilevel <- brms::brm(
  data = model_data,
  formula = formula_bilevel,
  prior = prior_bilevel,
  init = init_bilevel,
  control = list(
      adapt_delta = .99,
      max_treedepth = 12),
  stanvars = BayesPharma::sigmoid_stanvar(),
  cores = 4)

model_bilevel |> brms::expose_functions(vectorize = TRUE)

#| echo=FALSE,
#| message=FALSE,
#| dependson=c("fit-bilevel-model")
model_bilevel

#| echo=FALSE,
#| message=FALSE,
#| dependson=c("fit-model-conductance", "fit-bilevel-model")

model_MuSyC <- model_MuSyC |> brms::add_criterion("loo")
model_bilevel <- model_bilevel |> brms::add_criterion("loo")

brms::loo_compare(model_MuSyC, model_bilevel)

#| echo=FALSE,
#| message=FALSE,
#| dependson=c("fit-MuSyC-model", "fit-bilevel-model")

model_subsets01 <- data.frame(
  fraction = seq(0.01, 0.09, by=0.01)) |>
  dplyr::rowwise() |>
  dplyr::do({
    fraction <- .$fraction[1]
    cat(
      "Re-fitting the models with only ", fraction, " fraction of the data\n",
      sep = "")
    model_data_subset <- model_data |>
      dplyr::sample_frac(size = fraction)
    tibble::tibble(
      fraction = fraction,
      model_MuSyC = model_MuSyC |>
        stats::update(
          newdata = model_data_subset |>
            dplyr::mutate(
              logd1scale = mean(voltage),
              logd2scale = mean(logdose)),
          control = list(
            adapt_delta = .99,
            max_treedepth = 12),
          cores = 4,
          silent = TRUE) |>
        brms::add_criterion(
          "loo",
          moment_match = TRUE) |>
        list(),
      model_bilevel = model_bilevel |>
        stats::update(
          newdata = model_data_subset,
          cores = 4,
          control = list(
            adapt_delta = .99,
            max_treedepth = 12),
          silent = TRUE) |>
        brms::add_criterion(
          "loo",
          moment_match = TRUE) |>
        list())
  }) |>
  dplyr::ungroup()


model_subsets_loo_compare <- model_subsets01 |>
  dplyr::rowwise() |>
  dplyr::do({
    z <- .
    dplyr::bind_rows(
      data.frame(
        model_id = "model_MuSyC",
        elpd_loo = z$model_MuSyC$criteria$loo$estimates[1, 1],
        elpd_loo_se = z$model_MuSyC$criteria$loo$estimates[1, 2]),
      data.frame(
        model_id = "model_bilevel",
        elpd_loo = z$model_bilevel$criteria$loo$estimates[1, 1],
        elpd_loo_se = z$model_bilevel$criteria$loo$estimates[1, 2])) |>
      dplyr::mutate(
        fraction = z$fraction,
        .before = 1)
  })


ggplot2::ggplot(data = model_subsets_loo_compare) +
  ggplot2::theme_bw() +
  ggplot2::geom_ribbon(
    mapping = ggplot2::aes(
      x = fraction,
      ymin = elpd_loo - elpd_loo_se,
      ymax = elpd_loo + elpd_loo_se,
      fill = model_id,
      group = model_id),
    alpha = 0.3) +
  ggplot2::geom_line(
    mapping = ggplot2::aes(
      x = fraction,
      y = elpd_loo,
      color = model_id,
      group = model_id),
    linewidth = 1.2) +
  ggplot2::scale_color_discrete(
    "Model ID") +
  ggplot2::scale_fill_discrete(
    "Model ID") +
  ggplot2::scale_x_continuous(
    "Fraction of dataset") +
  ggplot2::scale_y_continuous(
    "Expected Log Probability of the Data (Higher is Better)")

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