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

# ─────────────────────────────────────────────
# mrgsolve model: Linezolid — 1-CMT IV Infusion
# PopPK in COVID-19-ARDS Patients on vv-ECMO
# Milaković et al. (2024) Pharmaceutics 16:253
# doi: 10.3390/pharmaceutics16020253
# ─────────────────────────────────────────────
model_code <- "
$PARAM @annotated
CL    : 5.9   : Clearance (L/h)
V     : 41.1  : Volume of distribution (L)
FUB   : 0.69  : Fraction unbound (unitless)

$CMT @annotated
CENT  : Central compartment (mg)

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

$TABLE
double Cp    = CENT / V;          // total plasma conc (mg/L = ug/mL)
double Cfree = Cp * FUB;          // free plasma concentration

$CAPTURE Cp Cfree
"

mod <- mcode("linezolid_popPK", model_code, quiet = TRUE)

# ─────────────────────────────────────────────
# 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; margin-top: 2px; }
  .metric-cmax  .metric-value { color: #8b5cf6; }
  .metric-ctrough .metric-value { color: #e74c3c; }
  .metric-auc   .metric-value { color: #10b981; }
  .metric-half  .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; }
  .footer-brand {
    text-align: center; padding: 12px; font-size: 11px;
    color: #adb5bd; border-top: 1px solid #dee2e6; margin-top: 16px;
  }
  .footer-brand a { color: #8b5cf6; text-decoration: none; }
  .pkpd-badge {
    display: inline-block; background: #8b5cf6; color: white;
    font-size: 10px; font-weight: 700; padding: 2px 7px;
    border-radius: 10px; letter-spacing: 0.5px; margin-left: 4px;
  }
")

# Helper — build one metric card div
metric_card <- function(id, css_class, label) {
  div(class = paste("metric-card", css_class),
    div(class = "metric-value", textOutput(id, inline = TRUE)),
    div(class = "metric-label", label)
  )
}

# ─────────────────────────────────────────────
# UI
# ─────────────────────────────────────────────
ui <- page_sidebar(
  title = "Linezolid PopPK Simulator — vv-ECMO / COVID-19",
  theme = app_theme,

  sidebar = sidebar(
    width = 290,
    h6("Dosing Regimen", class = "text-muted fw-bold"),
    radioButtons("regimen", "Dosing interval",
      choices = list("600 mg q8h (study)" = 8, "600 mg q12h (standard)" = 12),
      selected = 8
    ),
    sliderInput("dose",   "Dose (mg)",        min = 200, max = 1200, value = 600, step = 100),
    numericInput("infusion_dur", "Infusion duration (h)", value = 0.5, min = 0.1, max = 2, step = 0.1),
    hr(),
    h6("Population PK Parameters", class = "text-muted fw-bold"),
    sliderInput("CL", "Clearance — CL (L/h)", min = 1, max = 15, value = 5.9, step = 0.1),
    sliderInput("V",  "Volume of distribution — V (L)", min = 10, max = 100, value = 41.1, step = 0.5),
    hr(),
    h6("PK/PD Target", class = "text-muted fw-bold"),
    numericInput("mic", "MIC (mg/L)",     value = 2,    min = 0.03, max = 32,  step = 0.01),
    checkboxInput("log_scale", "Log-scale Y axis", value = FALSE),
    hr(),
    div(class = "footer-brand",
      "Powered by",
      tags$a(href = "https://pkpdbuilder.com", target = "_blank", "PKPDBuilder"),
      span(class = "pkpd-badge", "AI")
    )
  ),

  # ── Metric cards ──────────────────────────────
  layout_column_wrap(
    width = 1/4, fill = FALSE,
    metric_card("cmax_val",    "metric-card metric-cmax",    "Cmax (µg/mL)"),
    metric_card("ctrough_val", "metric-card metric-ctrough", "Ctrough (µg/mL)"),
    metric_card("auc24_val",   "metric-card metric-auc",     "fAUC24/MIC"),
    metric_card("half_val",    "metric-card metric-half",    "t½ (h)")
  ),

  # ── PK Plot ───────────────────────────────────
  card(
    full_screen = TRUE,
    card_header("Linezolid Plasma Concentration–Time Profile"),
    plotOutput("pk_plot", height = "480px")
  ),

  # ── PTA summary ──────────────────────────────
  card(
    card_header("PK/PD Target Assessment"),
    tableOutput("pta_table")
  ),

  # ── Tabs: Model Info / References ────────────
  navset_card_underline(
    nav_panel("Model Information",
      p(strong("Drug:"), "Linezolid (oxazolidinone antibiotic — Gram-positive coverage)"),
      p(strong("Structural model:"), "One-compartment model with first-order elimination and zero-order IV infusion input"),
      p(strong("Population:"), "Adult patients with COVID-19-associated ARDS (CARDS) on veno-venous ECMO (vv-ECMO)"),
      p(strong("Study:"), "Single-centre, prospective PK study (n = 9 patients; 53 steady-state plasma samples)"),
      p(strong("Data analysis:"), "NONMEM® v7.5, FOCE-I; bootstrap 95% CI (n = 1000 runs)"),
      tags$table(class = "table table-sm table-bordered mt-2",
        tags$thead(tags$tr(
          tags$th("Parameter"), tags$th("Estimate"), tags$th("Bootstrap 95% CI"), tags$th("IIV (%CV)")
        )),
        tags$tbody(
          tags$tr(tags$td("V (L)"),    tags$td("41.1"), tags$td("30.57–51.77"), tags$td("36.3%")),
          tags$tr(tags$td("CL (L/h)"), tags$td("5.9"),  tags$td("4.98–7.10"),  tags$td("24.8%")),
          tags$tr(tags$td("Proportional error"), tags$td("11.4%"), tags$td("7.9–14.8%"), tags$td("—")),
          tags$tr(tags$td("Covariance CL–V"), tags$td("−0.0901"), tags$td("−0.1353–−0.0212"), tags$td("—"))
        )
      ),
      p(strong("PK/PD targets (Linezolid):"),
        tags$ul(
          tags$li("85%Tf>MIC: time that free concentration exceeds MIC ≥ 85%"),
          tags$li("fAUC24/MIC ≥ 80: free-drug AUC over 24 h / MIC"),
          tags$li("Fraction unbound (fub) = 0.69 (literature value)")
        )
      ),
      p(strong("Therapeutic range:"), "Trough 2–8 µg/mL; trough > 8 µg/mL → elevated thrombocytopenia risk"),
      p(strong("Key finding:"), "600 mg q8h achieves ≥ 85%Tf>MIC in 91.5% of simulated ECMO patients vs 67.9% with q12h dosing (at MIC = 2 mg/L).")
    ),
    nav_panel("References",
      div(class = "ref-box",
        p(tags$b("Primary reference:")),
        p("Milaković D, Kovačević T, Kovačević P, et al. Population Pharmacokinetic Model of Linezolid and Probability of Target Attainment in Patients with COVID-19-Associated Acute Respiratory Distress Syndrome on Veno-Venous Extracorporeal Membrane Oxygenation. ",
          em("Pharmaceutics"), " 2024;16(2):253.",
          tags$a(href = "https://doi.org/10.3390/pharmaceutics16020253", target = "_blank", " DOI: 10.3390/pharmaceutics16020253")
        ),
        hr(),
        p(tags$b("Supporting references:")),
        tags$ul(
          tags$li("Bandín-Vilaret E, et al. A Review of Population Pharmacokinetic Analyses of Linezolid. Clin Pharmacokinet. 2022;61:789–817."),
          tags$li("Lin B, et al. Expert consensus statement on therapeutic drug monitoring and individualization of linezolid. Front Public Health 2022;10:967311."),
          tags$li("Roger C, et al. Clinical Pharmacokinetics and Pharmacodynamics of Oxazolidinones. Clin Pharmacokinet. 2018;57:559–575.")
        )
      )
    )
  )
)

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

  # ── Reactive: run mrgsolve simulation ────────
  sim_data <- reactive({
    shiny::req(input$CL, input$V, input$dose, input$infusion_dur, input$regimen, input$mic)

    ii  <- as.numeric(input$regimen)   # dosing interval (h)
    dur <- input$infusion_dur
    dose <- input$dose
    cl   <- input$CL
    v    <- input$V
    fub  <- 0.69

    # Simulate 5 dosing intervals (steady-state window), show last 2
    n_doses <- 5
    sim_end  <- n_doses * ii

    # Build event table — multiple IV infusions
    evs <- lapply(seq(0, (n_doses - 1) * ii, by = ii), function(t) {
      ev(amt = dose, cmt = 1, rate = dose / dur, time = t, evid = 1)
    })
    ev_combined <- do.call(c, evs)

    p <- param(mod, CL = cl, V = v, FUB = fub)

    out <- mrgsim(p, ev_combined,
                  tgrid = seq(0, sim_end, by = 0.1),
                  output = "df") |>
      dplyr::filter(time >= (n_doses - 2) * ii)

    # Normalise time to last 2 intervals
    out$time_norm <- out$time - (n_doses - 2) * ii
    out$Cp_free   <- out$Cfree
    out$AUC24     <- out$Cp_free   # placeholder — computed below
    out
  })

  # ── Derived PK metrics ───────────────────────
  pk_metrics <- reactive({
    shiny::req(sim_data())
    df  <- sim_data()
    ii  <- as.numeric(input$regimen)
    mic <- input$mic
    cl  <- input$CL
    v   <- input$V
    fub <- 0.69

    # Use the last dosing interval
    last_int <- df |> dplyr::filter(time_norm >= (max(df$time_norm) - ii) &
                                      time_norm <= max(df$time_norm))

    cmax   <- max(last_int$Cp, na.rm = TRUE)
    ctrough <- last_int$Cp[which.min(abs(last_int$time_norm - max(last_int$time_norm)))]

    # fAUC24 = free dose administered per day / CL  (analytical)
    doses_per_day <- 24 / ii
    fAUC24 <- (doses_per_day * input$dose * fub) / cl
    fAUC_MIC <- fAUC24 / mic

    # t½
    t_half <- 0.693 * v / cl

    # %Tf>MIC (free conc in last interval > MIC)
    free_conc <- last_int$Cp * fub
    tf_above  <- sum(free_conc > mic) / length(free_conc) * 100

    list(
      cmax    = round(cmax, 2),
      ctrough = round(ctrough, 2),
      fauc24  = round(fAUC24, 1),
      fauc_mic = round(fAUC_MIC, 1),
      t_half  = round(t_half, 2),
      tf_mic  = round(tf_above, 1)
    )
  })

  # ── Metric outputs ───────────────────────────
  output$cmax_val    <- renderText({ paste0(pk_metrics()$cmax) })
  output$ctrough_val <- renderText({ paste0(pk_metrics()$ctrough) })
  output$auc24_val   <- renderText({ paste0(pk_metrics()$fauc_mic) })
  output$half_val    <- renderText({ paste0(pk_metrics()$t_half) })

  # ── PK Plot ──────────────────────────────────
  output$pk_plot <- renderPlot({
    shiny::req(sim_data())
    df   <- sim_data()
    ii   <- as.numeric(input$regimen)
    mic  <- input$mic
    fub  <- 0.69

    # Add free conc column
    df$Cp_free <- df$Cp * fub

    # Reshape to long for ggplot
    df_long <- df |>
      dplyr::select(time_norm, Cp, Cp_free) |>
      tidyr::pivot_longer(cols = c(Cp, Cp_free),
                          names_to = "type",
                          values_to = "conc") |>
      dplyr::mutate(type = dplyr::case_when(
        type == "Cp"      ~ "Total plasma",
        type == "Cp_free" ~ "Free (unbound)"
      ))

    # Therapeutic window bands
    trough_low  <- 2
    trough_high <- 8

    p <- ggplot(df_long, aes(x = time_norm, y = conc, color = type, linetype = type)) +
      # Therapeutic window (shaded)
      annotate("rect", xmin = -Inf, xmax = Inf,
               ymin = trough_low, ymax = trough_high,
               fill = "#10b981", alpha = 0.08) +
      # Toxicity threshold
      geom_hline(yintercept = trough_high, linetype = "dotted", color = "#e74c3c", linewidth = 0.7) +
      # Lower bound
      geom_hline(yintercept = trough_low, linetype = "dotted", color = "#10b981", linewidth = 0.7) +
      # MIC line
      geom_hline(yintercept = mic * fub, linetype = "dashed", color = "#f59e0b", linewidth = 0.8) +
      # Concentration curves
      geom_line(linewidth = 1.1) +
      scale_color_manual(values = c("Total plasma" = "#8b5cf6", "Free (unbound)" = "#06b6d4")) +
      scale_linetype_manual(values = c("Total plasma" = "solid", "Free (unbound)" = "dashed")) +
      # Annotations
      annotate("text", x = max(df_long$time_norm) * 0.98, y = trough_high + 0.3,
               label = "Toxicity threshold (8 µg/mL)", hjust = 1, size = 3, color = "#e74c3c") +
      annotate("text", x = max(df_long$time_norm) * 0.98, y = mic * fub + 0.2,
               label = paste0("Free MIC (", round(mic * fub, 3), " µg/mL)"),
               hjust = 1, size = 3, color = "#f59e0b") +
      labs(
        x     = "Time (h)",
        y     = "Linezolid Concentration (µg/mL)",
        color = "Concentration type",
        linetype = "Concentration type",
        caption = paste0(
          input$dose, " mg IV q", ii, "h | CL = ", input$CL, " L/h | V = ", input$V, " L | ",
          "t½ = ", pk_metrics()$t_half, " h | %Tf>MIC = ", pk_metrics()$tf_mic, "%"
        )
      ) +
      theme_minimal(base_size = 13) +
      theme(
        legend.position  = "bottom",
        plot.caption     = element_text(size = 10, color = "gray50"),
        panel.grid.minor = element_blank()
      )

    # Log scale if requested
    if (isTRUE(input$log_scale)) {
      df_pos <- df_long |> dplyr::filter(conc > 0)
      if (nrow(df_pos) > 0) {
        p <- p %+% df_pos + scale_y_log10()
      }
    }

    p
  }, res = 110)

  # ── PTA Table ────────────────────────────────
  output$pta_table <- renderTable({
    m  <- pk_metrics()
    ii <- as.numeric(input$regimen)
    mic <- input$mic

    pass_icon <- function(x) if (x) "✅ PASS" else "❌ FAIL"

    data.frame(
      `PK/PD Target`         = c(
        "fAUC₂₄/MIC ≥ 80",
        "fAUC₂₄/MIC ≥ 100",
        "% Time free > MIC ≥ 85%",
        "Trough 2–8 µg/mL",
        "Trough < 8 µg/mL (safety)"
      ),
      `Simulated Value`      = c(
        round(m$fauc_mic, 1),
        round(m$fauc_mic, 1),
        paste0(m$tf_mic, "%"),
        round(m$ctrough, 2),
        round(m$ctrough, 2)
      ),
      `Target`               = c("≥ 80", "≥ 100", "≥ 85%", "2–8 µg/mL", "< 8 µg/mL"),
      `Assessment`           = c(
        pass_icon(m$fauc_mic >= 80),
        pass_icon(m$fauc_mic >= 100),
        pass_icon(m$tf_mic   >= 85),
        pass_icon(m$ctrough  >= 2 & m$ctrough <= 8),
        pass_icon(m$ctrough  <  8)
      ),
      check.names = FALSE
    )
  }, striped = TRUE, hover = TRUE, bordered = TRUE)
}

shinyApp(ui, server)