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

# ─────────────────────────────────────────────────────────────────────────────
# Antidepressant Tapering PK/PD Simulator
# Paper: "Alternate-day dosing to taper antidepressants risks severe withdrawal
#         effects: an in silico analysis" — O'Neill J et al., J Affect Disord 2025
#         DOI: 10.1016/j.jad.2025.015265
# 5 drugs: citalopram, escitalopram, sertraline, duloxetine, mirtazapine
# Model: 1-compartment oral PK + Michaelis-Menten receptor occupancy
# Ka values computed from Bateman equation (exact tmax solution)
# ─────────────────────────────────────────────────────────────────────────────

# PK parameters from Table 1; PD from Table 2 (paper)
# Ka solved numerically from Bateman: tmax = ln(Ka/Ke)/(Ka-Ke)
drug_params <- list(
  citalopram = list(
    label        = "Citalopram",
    V_L          = 840,      # 12 L/kg × 70 kg
    KA           = 1.034,    # Bateman exact: tmax=3.9h, t1/2=35h
    F            = 0.80,
    tmax_h       = 3.9,
    t12_h        = 35,
    V_Lkg        = 12,
    rec_type     = "SERT",
    Emax         = 83.98,
    EC50_ngml    = 4.412,
    default_dose = 20,
    min_dose     = 1,
    max_dose     = 80
  ),
  escitalopram = list(
    label        = "Escitalopram",
    V_L          = 1330,     # 19 L/kg × 70 kg
    KA           = 0.953,    # Bateman exact: tmax=4h, t1/2=30h
    F            = 0.80,
    tmax_h       = 4,
    t12_h        = 30,
    V_Lkg        = 19,
    rec_type     = "SERT",
    Emax         = 69.32,
    EC50_ngml    = 2.329,
    default_dose = 10,
    min_dose     = 1,
    max_dose     = 40
  ),
  sertraline = list(
    label        = "Sertraline",
    V_L          = 1400,     # 20 L/kg × 70 kg
    KA           = 0.4665,   # Bateman exact: tmax=6.5h, t1/2=25.9h
    F            = 0.756,
    tmax_h       = 6.5,
    t12_h        = 25.9,
    V_Lkg        = 20,
    rec_type     = "SERT",
    Emax         = 92.01,
    EC50_ngml    = 4.457,
    default_dose = 50,
    min_dose     = 1,
    max_dose     = 400
  ),
  duloxetine = list(
    label        = "Duloxetine",
    V_L          = 1640,     # 23.43 L/kg × 70 kg
    KA           = 0.365,    # Bateman exact: tmax=6h, t1/2=12h
    F            = 0.50,
    tmax_h       = 6,
    t12_h        = 12,
    V_Lkg        = 23.43,
    rec_type     = "SERT",
    Emax         = 90.75,
    EC50_ngml    = 2.696,
    default_dose = 60,
    min_dose     = 5,
    max_dose     = 240
  ),
  mirtazapine = list(
    label        = "Mirtazapine",
    V_L          = 264,      # 3.77 L/kg × 70 kg
    KA           = 2.864,    # Bateman exact: tmax=1.5h, t1/2=16.7h
    F            = 0.50,
    tmax_h       = 1.5,
    t12_h        = 16.7,
    V_Lkg        = 3.77,
    rec_type     = "alpha2",
    Emax         = 89.9,
    EC50_ngml    = 0.5961,
    default_dose = 15,
    min_dose     = 0.5,
    max_dose     = 90
  ),
  venlafaxine = list(
    label        = "Venlafaxine XR",
    V_L          = 525,      # 7.5 L/kg × 70 kg (Effexor XR PI)
    KA           = 0.4666,   # Bateman exact: tmax=5.5h, t1/2=15h (XR effective t½)
    F            = 0.45,     # 45% bioavailability (XR, Ellingrod & Perry 1994)
    tmax_h       = 5.5,
    t12_h        = 15,
    V_Lkg        = 7.5,
    rec_type     = "SERT",
    Emax         = 80.0,     # Sørensen et al. 2022 / Meyer et al. 2001 (PET)
    EC50_ngml    = 30.0,     # Lower SERT affinity vs SSRIs; ~30 ng/mL (Meyer 2001)
    default_dose = 75,
    min_dose     = 10,
    max_dose     = 375
  )
)

# Pre-compute CL for each drug (CL = V × ln2 / t½)
for (nm in names(drug_params)) {
  p <- drug_params[[nm]]
  p$CL_Lh <- p$V_L * log(2) / p$t12_h
  drug_params[[nm]] <- p
}

# ─── mrgsolve model: 1-CMT oral with Michaelis-Menten receptor occupancy ────
model_code <- "
$PARAM @annotated
CL      : 37.5    : Clearance (L/h)
V       : 1400.0  : Volume of distribution (L)
KA      : 0.4665  : First-order absorption rate (1/h)
F1      : 0.756   : Oral bioavailability
EMAX    : 92.01   : Maximum receptor occupancy (%)
EC50    : 4.457   : EC50 for receptor occupancy (ng/mL)

$CMT @annotated
DEPOT   : Oral depot (mg)
CENT    : Central compartment (mg)

$MAIN
F_DEPOT = F1;

$ODE
dxdt_DEPOT = -KA * DEPOT;
dxdt_CENT  = KA * DEPOT - (CL / V) * CENT;

$TABLE
double CP_ngml = (CENT / V) * 1000.0;
double CP_pos  = CP_ngml > 0.0 ? CP_ngml : 0.0;
double RO      = (EMAX * CP_pos) / (EC50 + CP_pos);

$CAPTURE @annotated
CP_pos : Plasma concentration (ng/mL)
RO     : Receptor occupancy (%)
"

mod <- mcode("antidepressant_tapering", model_code)

# ─── App 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: 22px; font-weight: bold; color: #2c3e50; }
  .metric-label  { font-size: 11px; 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; }
  .ref-box {
    background: #f0f4ff;
    border-left: 4px solid #8b5cf6;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
  .ref-box a { color: #8b5cf6; }
")

# ─── UI ──────────────────────────────────────────────────────────────────────
ui <- page_sidebar(
  title = "Antidepressant Tapering PK/PD Simulator",
  theme = app_theme,

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

    h6("Drug Selection"),
    selectInput(
      "drug", "Antidepressant",
      choices = c(
        "Citalopram"      = "citalopram",
        "Escitalopram"    = "escitalopram",
        "Sertraline"      = "sertraline",
        "Duloxetine"      = "duloxetine",
        "Mirtazapine"     = "mirtazapine",
        "Venlafaxine XR"  = "venlafaxine"
      ),
      selected = "venlafaxine"
    ),

    hr(),
    h6("Dosing Regimen"),
    numericInput("dose",       "Dose (mg)",            value = 50,  min = 0.1, max = 500, step = 0.5),
    numericInput("interval_h", "Dosing interval (h)",  value = 24,  min = 4,   max = 168, step = 4),
    numericInput("n_days",     "Simulation (days)",    value = 14,  min = 3,   max = 60),

    hr(),
    h6("Tapering Comparisons"),
    checkboxInput("compare_qod", "Overlay every-other-day (Q48H) dosing", value = FALSE),
    checkboxInput("compare_bd",  "Overlay twice-daily (Q12H) dosing",     value = FALSE),

    hr(),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE)
  ),

  # ── Metric Cards ────────────────────────────────────────────────────────────
  layout_column_wrap(
    width = 1/4,
    fill  = FALSE,
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("cmax_txt")),
      div(class = "metric-label", "Cmax,ss (ng/mL)")
    ),
    div(class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("ctrough_txt")),
      div(class = "metric-label", "Ctrough,ss (ng/mL)")
    ),
    div(class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("ro_max_txt")),
      div(class = "metric-label", "Max RO,ss (%)")
    ),
    div(class = "metric-card metric-info",
      div(class = "metric-value", textOutput("ro_var_txt")),
      div(class = "metric-label", "ΔRO within interval (pp)")
    )
  ),

  # ── Main Tabset ──────────────────────────────────────────────────────────────
  navset_card_underline(
    title = NULL,
    height = "550px",

    nav_panel(
      title = "PK Profile",
      plotOutput("pkPlot", height = "470px")
    ),

    nav_panel(
      title = "Receptor Occupancy",
      plotOutput("roPlot", height = "470px")
    ),

    nav_panel(
      title = "ΔRO Summary",
      div(style = "padding: 10px;",
        p(style = "font-size: 13px; color: #6c757d; margin-bottom: 10px;",
          "Receptor occupancy variation (ΔRO) across dosing intervals at the selected dose.
           Higher ΔRO = greater risk of withdrawal symptoms."),
        tableOutput("roTable")
      )
    ),

    nav_panel(
      title = "Model Info",
      div(style = "padding: 15px; overflow-y: auto; max-height: 470px;",
        markdown("
## Antidepressant Tapering PK/PD Simulator

This simulator models plasma concentration and receptor occupancy for five commonly
prescribed antidepressants to evaluate pharmacological tapering strategies.

### Model Structure
- **PK:** 1-compartment oral absorption (depot → central, first-order Ka)
- **PD:** Michaelis-Menten receptor occupancy equation

`RO(%) = Emax × Cp / (EC50 + Cp)` where Cp is in ng/mL

- SERT occupancy for: citalopram, escitalopram, sertraline, duloxetine
- α₂-adrenergic receptor occupancy for: mirtazapine

### Withdrawal Risk Thresholds (Maudsley Deprescribing Guidelines)

| ΔRO within dosing interval | Withdrawal Risk |
|---|---|
| < 5 percentage points | Mild |
| 5–10 pp | Moderate |
| > 10 pp | Severe |

### PK Parameters (Table 1)

| Drug | t½ (h) | V (L/kg) | Bioavailability | Tmax (h) |
|---|---|---|---|---|
| Citalopram | 35 | 12 | 80% | 3.9 |
| Escitalopram | 30 | 19 | 80% | 4.0 |
| Sertraline | 25.9 | 20 | 75.6% | 6.5 |
| Duloxetine | 12 | 23.43 | 50% | 6.0 |
| Mirtazapine | 16.7 | 3.77 | 50% | 1.5 |
| Venlafaxine XR | 15 | 7.5 | 45% | 5.5 |

### PD Parameters (Table 2)

| Drug | Emax (%) | EC50 (ng/mL) | Receptor |
|---|---|---|---|
| Citalopram | 83.98 | 4.412 | SERT |
| Escitalopram | 69.32 | 2.329 | SERT |
| Sertraline | 92.01 | 4.457 | SERT |
| Duloxetine | 90.75 | 2.696 | SERT |
| Mirtazapine | 89.9 | 0.5961 | α₂-adrenergic |
| Venlafaxine XR | 80.0 | 30.0 | SERT |

### Key Finding
Every-other-day (Q48H) dosing causes a **2–5× increase** in receptor occupancy
variation compared to daily (Q24H) dosing, substantially raising risk of
antidepressant withdrawal symptoms at all standard doses.
        ")
      )
    ),

    nav_panel(
      title = "References",
      div(style = "padding: 15px; overflow-y: auto; max-height: 470px;",
        div(class = "ref-box",
          tags$h5("Primary Publication"),
          tags$p(
            tags$strong(
              "O'Neill J, Horowitz MA, Sørensen A, et al. Alternate-day dosing to taper ",
              "antidepressants risks severe withdrawal effects: an in silico analysis."
            ),
            " J Affect Disord, 2025. ",
            tags$a(
              href = "https://doi.org/10.1016/j.jad.2025.015265",
              target = "_blank",
              "doi:10.1016/j.jad.2025.015265"
            )
          ),
          tags$h5("PK Data Sources"),
          tags$ul(
            tags$li(tags$strong("Citalopram:"), " Forest Pharmaceuticals Inc. (2009). Celexa PI."),
            tags$li(tags$strong("Escitalopram:"), " Aurobindo Pharma (2024). Lexapro PI."),
            tags$li(tags$strong("Sertraline:"), " Lupin Healthcare (2023); Alhadab & Brundage (2020) CPT; Kristensen et al. (1998)."),
            tags$li(tags$strong("Duloxetine:"), " Westanmo et al. (2005) Ann Pharmacother; Krka UK Ltd. (2023)."),
            tags$li(tags$strong("Mirtazapine:"), " Voortman & Paanakker (1995) J Clin Psychiatry.")
          ),
          tags$h5("Receptor Occupancy Data"),
          tags$ul(
            tags$li("Sørensen A et al. (2022). SERT Emax. J Psychopharmacol."),
            tags$li("Meyer JH et al. (2004). [11C]-DASB PET imaging."),
            tags$li("Smith DF et al. (2007). α₂ occupancy with mirtazapine. Eur Neuropsychopharmacol."),
            tags$li("Amann U et al. (2024). Duloxetine EC50. J Clin Psychopharmacol.")
          ),
          tags$h5("Clinical Guideline"),
          tags$ul(
            tags$li("Horowitz MA, Taylor D (2024). Maudsley Deprescribing Guidelines.")
          )
        )
      )
    )
  ),

  # ── 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"
    ),
    " • Built by Sunny ☀️ (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 for a given drug + regimen ───────────────────────
run_sim_drug <- function(p, dose_mg, interval_h, n_days) {
  n_doses <- max(1L, as.integer(floor(n_days * 24 / interval_h)))
  addl    <- n_doses - 1L
  ev1 <- ev(amt = dose_mg, cmt = 1, ii = interval_h, addl = addl)
  mod %>%
    param(
      CL   = p$CL_Lh,
      V    = p$V_L,
      KA   = p$KA,
      F1   = p$F,
      EMAX = p$Emax,
      EC50 = p$EC50_ngml
    ) %>%
    ev(ev1) %>%
    mrgsim(end = n_days * 24, delta = 0.25) %>%
    as.data.frame()
}

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

  dp <- reactive({
    shiny::req(input$drug)
    drug_params[[input$drug]]
  })

  # Update dose range/default when drug changes
  observeEvent(input$drug, {
    p <- drug_params[[input$drug]]
    updateNumericInput(session, "dose",
      value = p$default_dose,
      min   = p$min_dose,
      max   = p$max_dose
    )
  })

  # ── Primary simulation ────────────────────────────────────────────────────
  sim_primary <- reactive({
    shiny::req(input$dose, input$interval_h, input$n_days)
    d <- run_sim_drug(dp(), input$dose, input$interval_h, input$n_days)
    d$regimen <- paste0("Q", input$interval_h, "H — ", input$dose, " mg")
    d
  })

  # ── QOD (Q48H) overlay ───────────────────────────────────────────────────
  sim_qod <- reactive({
    shiny::req(input$dose, input$n_days)
    d <- run_sim_drug(dp(), input$dose, 48, input$n_days)
    d$regimen <- paste0("Q48H (every-other-day) — ", input$dose, " mg")
    d
  })

  # ── Q12H (twice-daily) overlay ───────────────────────────────────────────
  sim_bd <- reactive({
    shiny::req(input$dose, input$n_days)
    half_dose <- input$dose / 2
    d <- run_sim_drug(dp(), half_dose, 12, input$n_days)
    d$regimen <- paste0("Q12H — ", half_dose, " mg/dose (same daily total)")
    d
  })

  # Combined plot data
  plot_data <- reactive({
    d <- sim_primary()
    if (input$compare_qod) d <- bind_rows(d, sim_qod())
    if (input$compare_bd)  d <- bind_rows(d, sim_bd())
    d
  })

  # Steady-state window (last dosing interval of primary sim)
  ss_data <- reactive({
    d <- sim_primary()
    # Use the actual last dose time based on chosen interval (not always n_days-1 days)
    last_start <- floor((input$n_days * 24 / input$interval_h) - 1) * input$interval_h
    d |> dplyr::filter(time >= last_start)
  })

  # ── Metric outputs ────────────────────────────────────────────────────────
  output$cmax_txt <- renderText({
    d <- ss_data()
    shiny::req(nrow(d) > 0)
    sprintf("%.2f", max(d$CP_pos, na.rm = TRUE))
  })

  output$ctrough_txt <- renderText({
    d <- ss_data()
    shiny::req(nrow(d) > 0)
    sprintf("%.2f", min(d$CP_pos, na.rm = TRUE))
  })

  output$ro_max_txt <- renderText({
    d <- ss_data()
    shiny::req(nrow(d) > 0)
    sprintf("%.1f%%", max(d$RO, na.rm = TRUE))
  })

  output$ro_var_txt <- renderText({
    d <- ss_data()
    shiny::req(nrow(d) > 0)
    delta_ro <- max(d$RO, na.rm = TRUE) - min(d$RO, na.rm = TRUE)
    sprintf("%.1f pp", delta_ro)
  })

  # ── PK Concentration Plot ─────────────────────────────────────────────────
  output$pkPlot <- renderPlot({
    d <- plot_data()
    shiny::req(nrow(d) > 0)

    d_plot <- if (input$log_scale) {
      d |> dplyr::filter(CP_pos > 0.001)
    } else {
      d
    }

    p_info <- dp()
    p <- ggplot(d_plot, aes(x = time / 24, y = CP_pos, colour = regimen)) +
      geom_line(linewidth = 0.9) +
      scale_colour_manual(
        values = c("#8b5cf6", "#f59e0b", "#10b981"),
        name   = "Regimen"
      ) +
      labs(
        x        = "Time (days)",
        y        = "Plasma Concentration (ng/mL)",
        title    = paste0(p_info$label, " — Concentration-Time Profile"),
        subtitle = paste0(
          "t\u00bd = ", p_info$t12_h, " h  |  Vd = ", p_info$V_Lkg,
          " L/kg  |  F = ", round(p_info$F * 100), "%"
        )
      ) +
      theme_minimal(base_size = 13) +
      theme(legend.position = "bottom")

    if (input$log_scale) p <- p + scale_y_log10()
    p
  })

  # ── Receptor Occupancy Plot ───────────────────────────────────────────────
  output$roPlot <- renderPlot({
    d <- plot_data()
    shiny::req(nrow(d) > 0)

    p_info   <- dp()
    rec_lbl  <- if (p_info$rec_type == "SERT") "SERT Occupancy (%)" else "\u03b1\u2082-Adrenergic Receptor Occupancy (%)"

    ggplot(d, aes(x = time / 24, y = RO, colour = regimen)) +
      # Risk band shading
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0,  ymax = 100, fill = "#fee2e2", alpha = 0.06) +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0,  ymax = 75,  fill = "#fef3c7", alpha = 0.06) +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0,  ymax = 50,  fill = "#d1fae5", alpha = 0.08) +
      geom_line(linewidth = 0.9) +
      scale_colour_manual(
        values = c("#8b5cf6", "#f59e0b", "#10b981"),
        name   = "Regimen"
      ) +
      scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 20)) +
      labs(
        x        = "Time (days)",
        y        = rec_lbl,
        title    = paste0(p_info$label, " — Receptor Occupancy"),
        subtitle = paste0(
          "Emax = ", p_info$Emax, "%  |  EC\u2085\u2080 = ", p_info$EC50_ngml, " ng/mL"
        )
      ) +
      theme_minimal(base_size = 13) +
      theme(legend.position = "bottom")
  })

  # ── ΔRO Table: all intervals ─────────────────────────────────────────────
  output$roTable <- renderTable({
    shiny::req(input$dose, input$n_days)
    p <- dp()

    intervals <- c(12, 24, 48, 72, 96, 168)
    int_labels <- c("Q12H", "Q24H (daily)", "Q48H (QOD)", "Q72H", "Q96H", "Q168H (weekly)")

    rows <- mapply(function(ii_h, lbl) {
      d <- run_sim_drug(p, input$dose, ii_h, input$n_days)
      last_start <- floor((input$n_days * 24 / ii_h) - 1) * ii_h
      ss <- d |> dplyr::filter(time >= last_start)
      if (nrow(ss) < 2) return(NULL)

      ro_max <- max(ss$RO, na.rm = TRUE)
      ro_min <- min(ss$RO, na.rm = TRUE)
      delta  <- ro_max - ro_min

      risk <- dplyr::case_when(
        delta <  5  ~ "Mild",
        delta <= 10 ~ "Moderate",
        TRUE        ~ "Severe"
      )

      data.frame(
        "Interval"         = lbl,
        "RO max (%)"       = sprintf("%.1f", ro_max),
        "RO min (%)"       = sprintf("%.1f", ro_min),
        "\u0394RO (pp)"    = sprintf("%.1f", delta),
        "Withdrawal Risk"  = risk,
        check.names = FALSE
      )
    }, intervals, int_labels, SIMPLIFY = FALSE)

    do.call(rbind, Filter(Negate(is.null), rows))
  }, striped = TRUE, hover = TRUE, bordered = TRUE, width = "100%")
}

shinyApp(ui = ui, server = server)