library(shiny)
library(bslib)
library(mrgsolve)
library(dplyr)
library(ggplot2)

# ── PD Model: Meropenem/Gentamicin/Ciprofloxacin Bacterial Killing ──
# Source: Sadouki et al. (2025) Sci Rep 15:45244
# E. coli NCTC 12241 static time-kill model (nlmixr2)

model_code <- '
$PARAM @annotated
// Growth parameters
knet   : 1.35    : Net bacterial growth rate (1/h)
BMAX   : 10.0    : Log10 max carrying capacity (log10 CFU/mL)
B0     : 5.29    : Log10 initial inoculum (log10 CFU/mL)

// Meropenem PD parameters
EMAX_MER  : 4.18   : Meropenem Emax
IC50_MER  : 0.0781 : Meropenem IC50 (mg/L)
HILL_MER  : 2.76   : Meropenem Hill coefficient
BETA_MER  : 0.922  : Meropenem regrowth BETA
TAU_MER   : 0.570  : Meropenem regrowth TAU
TAU_MER_COEFF : 0.00155 : MER conc effect on TAU

// Gentamicin PD parameters
EMAX_GEN  : 5.47   : Gentamicin Emax
IC50_GEN  : 1.12   : Gentamicin IC50 (mg/L)
HILL_GEN  : 3.63   : Gentamicin Hill coefficient
BETA_GEN  : 0.829  : Gentamicin regrowth BETA
TAU_GEN   : 0.517  : Gentamicin regrowth TAU
// Emax model on GEN BETA
BETA_GEN_IC50 : 5.72  : IC50 for GEN conc on BETA
BETA_GEN_EMAX : -2.97 : Emax for GEN conc on BETA

// Ciprofloxacin PD parameters
EMAX_CIP  : 4.55   : Ciprofloxacin Emax
IC50_CIP  : 0.0106 : Ciprofloxacin IC50 (mg/L)
HILL_CIP  : 3.58   : Ciprofloxacin Hill coefficient
BETA_CIP  : 0.674  : Ciprofloxacin regrowth BETA
TAU_CIP   : 0.359  : Ciprofloxacin regrowth TAU
// Emax model on CIP BETA
BETA_CIP_IC50 : 0.017 : IC50 for CIP conc on BETA
BETA_CIP_EMAX : -4.0  : Emax for CIP conc on BETA

// Drug concentrations (static, set by user)
C_MER : 0.0  : Meropenem concentration (mg/L)
C_GEN : 0.0  : Gentamicin concentration (mg/L)
C_CIP : 0.0  : Ciprofloxacin concentration (mg/L)

// Interaction parameters
COMBO_BETA_COEFF : -1.0 : Categorical combo effect on BETA
IC50_CIP_MER_COEFF : -0.353 : MER presence effect on CIP IC50
IC50_CIP_GEN_COEFF : -0.576 : GEN presence effect on CIP IC50

// Meropenem degradation
KDEG_MER : 0.01316 : Meropenem degradation rate (1/h)

// Inoculum effect
LOW_INOC : 0 : Low inoculum flag (0=10^5, 1=10^3)
INOC_B0  : -0.326 : Low inoculum effect on B0
INOC_BMAX : -0.11  : Low inoculum effect on BMAX

$CMT @annotated
BACT : Bacterial load (log10 CFU/mL)
CMER : Meropenem concentration (mg/L)

$MAIN
double B0_adj = B0 + LOW_INOC * INOC_B0;
double BMAX_adj = BMAX + LOW_INOC * INOC_BMAX;
BACT_0 = pow(10.0, B0_adj);
CMER_0 = C_MER;

$ODE
// Current meropenem concentration (degrades over time)
double cmer_now = CMER > 1e-12 ? CMER : 0.0;
double cgen_now = C_GEN;
double ccip_now = C_CIP;

// Determine if combination (2+ drugs present)
int n_drugs = 0;
if (cmer_now > 1e-12) n_drugs++;
if (cgen_now > 1e-12) n_drugs++;
if (ccip_now > 1e-12) n_drugs++;
int is_combo = n_drugs >= 2 ? 1 : 0;

// Meropenem effect
double eff_mer = 0.0;
if (cmer_now > 1e-12) {
  double conc_ratio_mer = pow(cmer_now, HILL_MER) / (pow(IC50_MER, HILL_MER) + pow(cmer_now, HILL_MER));
  double tau_mer_eff = TAU_MER + TAU_MER_COEFF * cmer_now;
  double beta_mer_eff = BETA_MER;
  if (is_combo) beta_mer_eff = beta_mer_eff + COMBO_BETA_COEFF;
  double wane_mer = 1.0 - beta_mer_eff * (1.0 - exp(-SOLVERTIME * tau_mer_eff));
  if (wane_mer < 0.0) wane_mer = 0.0;
  eff_mer = EMAX_MER * conc_ratio_mer * wane_mer;
}

// Gentamicin effect
double eff_gen = 0.0;
if (cgen_now > 1e-12) {
  double conc_ratio_gen = pow(cgen_now, HILL_GEN) / (pow(IC50_GEN, HILL_GEN) + pow(cgen_now, HILL_GEN));
  // Emax model on BETA for gentamicin
  double beta_gen_drug = BETA_GEN_EMAX * pow(cgen_now, 20.0) / (pow(BETA_GEN_IC50, 20.0) + pow(cgen_now, 20.0));
  double beta_gen_eff = BETA_GEN + beta_gen_drug;
  if (is_combo) beta_gen_eff = beta_gen_eff + COMBO_BETA_COEFF;
  // TAU with inverse concentration relationship
  double tau_gen_eff = TAU_GEN;
  if (cgen_now > 0.01) tau_gen_eff = TAU_GEN / cgen_now;
  double wane_gen = 1.0 - beta_gen_eff * (1.0 - exp(-SOLVERTIME * tau_gen_eff));
  if (wane_gen < 0.0) wane_gen = 0.0;
  eff_gen = EMAX_GEN * conc_ratio_gen * wane_gen;
}

// Ciprofloxacin effect
double eff_cip = 0.0;
if (ccip_now > 1e-12) {
  // Interaction: MER/GEN presence reduces CIP IC50
  double ic50_cip_eff = IC50_CIP;
  if (cmer_now > 1e-12) ic50_cip_eff = ic50_cip_eff * (1.0 + IC50_CIP_MER_COEFF);
  if (cgen_now > 1e-12) ic50_cip_eff = ic50_cip_eff * (1.0 + IC50_CIP_GEN_COEFF);
  double conc_ratio_cip = pow(ccip_now, HILL_CIP) / (pow(ic50_cip_eff, HILL_CIP) + pow(ccip_now, HILL_CIP));
  // Emax model on BETA for ciprofloxacin
  double beta_cip_drug = BETA_CIP_EMAX * pow(ccip_now, 20.0) / (pow(BETA_CIP_IC50, 20.0) + pow(ccip_now, 20.0));
  double beta_cip_eff = BETA_CIP + beta_cip_drug;
  if (is_combo) beta_cip_eff = beta_cip_eff + COMBO_BETA_COEFF;
  double wane_cip = 1.0 - beta_cip_eff * (1.0 - exp(-SOLVERTIME * TAU_CIP));
  if (wane_cip < 0.0) wane_cip = 0.0;
  eff_cip = EMAX_CIP * conc_ratio_cip * wane_cip;
}

// Total drug effect (additive)
double total_effect = eff_mer + eff_gen + eff_cip;

// Bacterial dynamics: logistic growth - killing
double bact_now = BACT > 1.0 ? BACT : 1.0;  // Floor at 1 CFU/mL
double growth = knet * (1.0 - bact_now / pow(10.0, BMAX_adj));
double net_rate = growth - total_effect;
dxdt_BACT = net_rate * bact_now;

// Meropenem degradation
dxdt_CMER = -KDEG_MER * CMER;

$TABLE
double LOG10_CFU = log10(BACT > 1.0 ? BACT : 1.0);
double CMER_OUT = CMER > 0.0 ? CMER : 0.0;

$CAPTURE @annotated
LOG10_CFU : Bacterial load (log10 CFU/mL)
CMER_OUT  : Meropenem concentration (mg/L)
'

mod <- mcode("pd_interaction", model_code)

# ── MIC values ──
MIC_MER <- 0.03   # mg/L
MIC_GEN <- 1.0    # mg/L
MIC_CIP <- 0.015  # mg/L

# ── Theme ──
app_theme <- bs_theme(
  version = 5,
  bootswatch = "flatly",
  primary = "#8b5cf6"
) |> bs_add_rules("
  .metric-card {
    background: #f8f9fa;
    border-radius: 8px;
    padding: 15px;
    margin: 5px;
    text-align: center;
    border: 1px solid #dee2e6;
  }
  .metric-value {
    font-size: 24px;
    font-weight: bold;
    color: #2c3e50;
  }
  .metric-label {
    font-size: 12px;
    color: #7f8c8d;
  }
  .metric-success .metric-value { color: #10b981; }
  .metric-warning .metric-value { color: #f59e0b; }
  .metric-primary .metric-value { color: #8b5cf6; }
  .metric-info .metric-value { color: #0dcaf0; }
  .metric-danger .metric-value { color: #ef4444; }
  .drug-badge {
    display: inline-block;
    padding: 3px 10px;
    border-radius: 12px;
    font-size: 11px;
    font-weight: 600;
    margin: 2px;
  }
  .badge-mer { background: #dbeafe; color: #1e40af; }
  .badge-gen { background: #dcfce7; color: #166534; }
  .badge-cip { background: #fef3c7; color: #92400e; }
")

# ── UI ──
ui <- page_sidebar(
 title = "Antibiotic PD Interaction Simulator",
  theme = app_theme,

  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,

    h6("Antibiotic Selection"),
    checkboxInput("use_mer", span("Meropenem", class = "drug-badge badge-mer"), value = TRUE),
    conditionalPanel(
      condition = "input.use_mer",
      sliderInput("mic_mult_mer", "Meropenem (× MIC)", min = 0.25, max = 16, value = 4, step = 0.25,
                  post = paste0(" × ", MIC_MER, " mg/L"))
    ),

    checkboxInput("use_gen", span("Gentamicin", class = "drug-badge badge-gen"), value = FALSE),
    conditionalPanel(
      condition = "input.use_gen",
      sliderInput("mic_mult_gen", "Gentamicin (× MIC)", min = 0.25, max = 16, value = 4, step = 0.25,
                  post = paste0(" × ", MIC_GEN, " mg/L"))
    ),

    checkboxInput("use_cip", span("Ciprofloxacin", class = "drug-badge badge-cip"), value = FALSE),
    conditionalPanel(
      condition = "input.use_cip",
      sliderInput("mic_mult_cip", "Ciprofloxacin (× MIC)", min = 0.25, max = 16, value = 4, step = 0.25,
                  post = paste0(" × ", MIC_CIP, " mg/L"))
    ),

    hr(),
    h6("Experimental Conditions"),
    radioButtons("inoculum", "Inoculum Size",
                 choices = c("10^5 CFU/mL (standard)" = "high", "10^3 CFU/mL (low)" = "low"),
                 selected = "high"),

    hr(),
    checkboxInput("show_growth", "Show growth control", value = TRUE),
    checkboxInput("show_combo_compare", "Compare all 2-way combos", value = FALSE),

    hr(),
    card(
      card_header("MIC Values", class = "bg-light"),
      p(strong("Meropenem:"), "0.03 mg/L"),
      p(strong("Gentamicin:"), "1.0 mg/L"),
      p(strong("Ciprofloxacin:"), "0.015 mg/L"),
      p(em("E. coli NCTC 12241"))
    )
  ),

  # Metrics row
  layout_column_wrap(
        width = 1/4,
        fill = FALSE,
    div(class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("max_kill")),
      div(class = "metric-label", "Max Kill (Δlog10)")
    ),
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("nadir")),
      div(class = "metric-label", "Nadir (log10 CFU/mL)")
    ),
    div(class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("cfu_24h")),
      div(class = "metric-label", "24h Load (log10 CFU/mL)")
    ),
    div(class = "metric-card metric-danger",
      div(class = "metric-value", textOutput("regrowth_status")),
      div(class = "metric-label", "Regrowth at 24h")
    )
  ),

  # Main time-kill plot
  card(
    card_header("Time-Kill Curve"),
    full_screen = TRUE,
    plotOutput("tkPlot", height = "480px")
  ),

  # Tabs
  navset_card_tab(
    title = "Analysis",
    nav_panel("Concentration-Effect",
      plotOutput("ceplot", height = "400px")
    ),
    nav_panel("Combo Comparison",
      plotOutput("comboPlot", height = "450px")
    ),
    nav_panel("Model Information",
      markdown(paste0(
        "## Pharmacodynamic Interaction Model\n\n",
        "**Source:** Sadouki Z, Wey EQ, Read L, et al. *Pharmacodynamic interactions among meropenem, ",
        "ciprofloxacin and gentamicin in an in-vitro model.* Scientific Reports. 2025;15:45244.\n\n",
        "**DOI:** [10.1038/s41598-025-29354-y](https://doi.org/10.1038/s41598-025-29354-y)\n\n",
        "### Model Structure\n\n",
        "- **Bacterial Growth:** Logistic growth model (knet = 1.35 h⁻¹, Bmax = 10¹⁰ CFU/mL)\n",
        "- **Drug Effect:** Sigmoidal Emax model for each drug\n",
        "- **Regrowth:** Time-dependent waning of drug effect (β, τ parameters)\n",
        "- **Interactions:** Additive drug effects; synergy via reduced ciprofloxacin IC50\n",
        "- **Meropenem Degradation:** First-order (≈10% loss in 8h at 37°C)\n\n",
        "### Key Findings\n\n",
        "- **Gentamicin + Ciprofloxacin:** Synergistic — 57.6% reduction in CIP IC50\n",
        "- **Meropenem + Ciprofloxacin:** Synergistic — 35.3% reduction in CIP IC50\n",
        "- **Meropenem + Gentamicin:** Indifferent (no synergy on IC50)\n",
        "- **All 2/3-way combos:** Prevent bacterial regrowth (vs. monotherapy)\n",
        "- **Three-way vs Two-way:** No additional benefit of adding meropenem to GEN+CIP\n\n",
        "### References\n\n",
        "- [Original Paper (Nature Scientific Reports)](https://doi.org/10.1038/s41598-025-29354-y)\n",
        "- [Meropenem FDA Label](https://www.accessdata.fda.gov/drugsatfda_docs/label/2017/209776lbl.pdf)\n",
        "- [CLSI M26-A: Time-Kill Methods](https://clsi.org/)\n\n",
        "### Therapeutic Context\n\n",
        "This combination is used as empiric broad-spectrum therapy in severe sepsis, ",
        "particularly for MDR gram-negative infections. Recent MHRA regulations restrict ",
        "quinolone use, making understanding of PD interactions critical for antibiotic stewardship.\n"
      ))
    )
  ),

  # Footer
  div(style = "text-align: center; padding: 20px; margin-top: 30px;
      border-top: 1px solid #e9ecef; color: #6c757d; font-size: 12px;",
    "Powered by ",
    tags$a(href = "https://www.pkpdbuilder.com", target = "_blank",
           style = "color: #8b5cf6; font-weight: 500;", "PKPDBuilder.com"),
    " \u2022 Built by Sunny \u2600\ufe0f (Husain Attarwala's AI Assistant)",
    br(),
    tags$span(style = "font-size: 10px;",
      "For research and educational purposes only. Not for clinical decision-making.")
  )
)

# ── Helper: run simulation ──
run_sim <- function(mod, c_mer, c_gen, c_cip, low_inoc = 0, end_time = 24) {
  out <- mod %>%
    param(C_MER = c_mer, C_GEN = c_gen, C_CIP = c_cip, LOW_INOC = low_inoc) %>%
    mrgsim(end = end_time, delta = 0.1) %>%
    as.data.frame()
  out
}

# ── Server ──
server <- function(input, output, session) {

  low_inoc <- reactive({ if (input$inoculum == "low") 1 else 0 })

  sim_data <- reactive({
    c_mer <- if (input$use_mer) input$mic_mult_mer * MIC_MER else 0
    c_gen <- if (input$use_gen) input$mic_mult_gen * MIC_GEN else 0
    c_cip <- if (input$use_cip) input$mic_mult_cip * MIC_CIP else 0
    run_sim(mod, c_mer, c_gen, c_cip, low_inoc())
  })

  growth_data <- reactive({
    run_sim(mod, 0, 0, 0, low_inoc())
  })

  # Metrics
  output$max_kill <- renderText({
    d <- sim_data()
    b0 <- d$LOG10_CFU[1]
    nadir <- min(d$LOG10_CFU)
    sprintf("%.1f", b0 - nadir)
  })

  output$nadir <- renderText({
    sprintf("%.1f", min(sim_data()$LOG10_CFU))
  })

  output$cfu_24h <- renderText({
    d <- sim_data()
    sprintf("%.1f", tail(d$LOG10_CFU, 1))
  })

  output$regrowth_status <- renderText({
    d <- sim_data()
    nadir <- min(d$LOG10_CFU)
    final <- tail(d$LOG10_CFU, 1)
    if (final > nadir + 1) "Yes" else "No"
  })

  # Time-kill plot
  output$tkPlot <- renderPlot({
    d <- sim_data()

    # Build label
    drugs <- c()
    if (input$use_mer) drugs <- c(drugs, paste0("MER ", input$mic_mult_mer, "×MIC"))
    if (input$use_gen) drugs <- c(drugs, paste0("GEN ", input$mic_mult_gen, "×MIC"))
    if (input$use_cip) drugs <- c(drugs, paste0("CIP ", input$mic_mult_cip, "×MIC"))
    label <- if (length(drugs) > 0) paste(drugs, collapse = " + ") else "No Drug"

    d$Regimen <- label

    p <- ggplot(d, aes(x = time, y = LOG10_CFU)) +
      geom_line(aes(color = Regimen), linewidth = 1.2)

    if (input$show_growth) {
      gd <- growth_data()
      gd$Regimen <- "Growth Control"
      p <- p + geom_line(data = gd, aes(x = time, y = LOG10_CFU, color = Regimen),
                         linewidth = 1, linetype = "dashed")
    }

    # LOD line
    p <- p +
      geom_hline(yintercept = 1, linetype = "dotted", color = "grey50") +
      annotate("text", x = 0.5, y = 1.3, label = "LOD", color = "grey50", size = 3) +
      labs(x = "Time (hours)", y = expression(log[10]~"CFU/mL"),
           title = "Static Time-Kill Curve",
           subtitle = "E. coli NCTC 12241") +
      scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 12, 16, 20, 24)) +
      scale_y_continuous(limits = c(0, 12)) +
      scale_color_manual(values = c("Growth Control" = "grey40",
                                    setNames("#8b5cf6", label))) +
      theme_minimal(base_size = 14) +
      theme(legend.position = "top",
            legend.title = element_blank(),
            panel.grid.minor = element_blank())
    p
  })

  # Concentration-effect plot
  output$ceplot <- renderPlot({
    mic_mults <- c(0.25, 0.5, 1, 2, 4, 8, 16)
    results <- data.frame()

    for (m in mic_mults) {
      # Monotherapies
      if (input$use_mer || (!input$use_mer && !input$use_gen && !input$use_cip)) {
        d <- run_sim(mod, m * MIC_MER, 0, 0, low_inoc())
        results <- rbind(results, data.frame(
          MIC_mult = m, Drug = "Meropenem",
          Nadir = min(d$LOG10_CFU), CFU_24h = tail(d$LOG10_CFU, 1)
        ))
      }
      if (input$use_gen || (!input$use_mer && !input$use_gen && !input$use_cip)) {
        d <- run_sim(mod, 0, m * MIC_GEN, 0, low_inoc())
        results <- rbind(results, data.frame(
          MIC_mult = m, Drug = "Gentamicin",
          Nadir = min(d$LOG10_CFU), CFU_24h = tail(d$LOG10_CFU, 1)
        ))
      }
      if (input$use_cip || (!input$use_mer && !input$use_gen && !input$use_cip)) {
        d <- run_sim(mod, 0, 0, m * MIC_CIP, low_inoc())
        results <- rbind(results, data.frame(
          MIC_mult = m, Drug = "Ciprofloxacin",
          Nadir = min(d$LOG10_CFU), CFU_24h = tail(d$LOG10_CFU, 1)
        ))
      }
    }

    shiny::req(nrow(results) > 0)

    ggplot(results, aes(x = MIC_mult, y = Nadir, color = Drug)) +
      geom_line(linewidth = 1) +
      geom_point(size = 2.5) +
      scale_x_log10(breaks = mic_mults, labels = paste0(mic_mults, "×")) +
      labs(x = "Concentration (× MIC)", y = expression("Nadir " * log[10] ~ "CFU/mL"),
           title = "Concentration-Effect Relationship",
           subtitle = "Monotherapy nadir bacterial load across MIC multiples") +
      scale_color_manual(values = c("Meropenem" = "#3b82f6", "Gentamicin" = "#22c55e",
                                    "Ciprofloxacin" = "#f59e0b")) +
      theme_minimal(base_size = 14) +
      theme(legend.position = "top", legend.title = element_blank(),
            panel.grid.minor = element_blank())
  })

  # Combo comparison plot
  output$comboPlot <- renderPlot({
    mic_mult <- if (input$use_mer) input$mic_mult_mer else 4
    if (input$use_gen) mic_mult <- input$mic_mult_gen
    if (input$use_cip) mic_mult <- input$mic_mult_cip

    combos <- list(
      "Growth Control" = c(0, 0, 0),
      "MER alone" = c(mic_mult * MIC_MER, 0, 0),
      "GEN alone" = c(0, mic_mult * MIC_GEN, 0),
      "CIP alone" = c(0, 0, mic_mult * MIC_CIP),
      "MER + GEN" = c(mic_mult * MIC_MER, mic_mult * MIC_GEN, 0),
      "MER + CIP" = c(mic_mult * MIC_MER, 0, mic_mult * MIC_CIP),
      "GEN + CIP" = c(0, mic_mult * MIC_GEN, mic_mult * MIC_CIP),
      "MER + GEN + CIP" = c(mic_mult * MIC_MER, mic_mult * MIC_GEN, mic_mult * MIC_CIP)
    )

    all_data <- data.frame()
    for (nm in names(combos)) {
      cc <- combos[[nm]]
      d <- run_sim(mod, cc[1], cc[2], cc[3], low_inoc())
      d$Regimen <- nm
      all_data <- rbind(all_data, d[, c("time", "LOG10_CFU", "Regimen")])
    }

    all_data$Regimen <- factor(all_data$Regimen, levels = names(combos))

    colors <- c("Growth Control" = "grey40",
                "MER alone" = "#93c5fd", "GEN alone" = "#86efac", "CIP alone" = "#fcd34d",
                "MER + GEN" = "#3b82f6", "MER + CIP" = "#f97316",
                "GEN + CIP" = "#22c55e", "MER + GEN + CIP" = "#8b5cf6")

    ggplot(all_data, aes(x = time, y = LOG10_CFU, color = Regimen)) +
      geom_line(linewidth = 0.9) +
      geom_hline(yintercept = 1, linetype = "dotted", color = "grey50") +
      scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 12, 16, 20, 24)) +
      scale_y_continuous(limits = c(0, 12)) +
      scale_color_manual(values = colors) +
      labs(x = "Time (hours)", y = expression(log[10]~"CFU/mL"),
           title = paste0("All Regimen Comparison at ", mic_mult, "× MIC"),
           subtitle = "E. coli NCTC 12241 — Monotherapy vs Combinations") +
      theme_minimal(base_size = 14) +
      theme(legend.position = "right", legend.title = element_blank(),
            panel.grid.minor = element_blank())
  })
}

shinyApp(ui = ui, server = server)