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

# ── mrgsolve model ─────────────────────────────────────────────────────────────
# One-compartment model: IV infusion, first-order elimination
# Population PK of piperacillin in neonates with LOS / severe infections
# CL: allometric BW + Rhodin sigmoidal PMA maturation + SCr power covariate
# V:  allometric BW scaling only
# Ref: Rhodin MM et al. Pediatr Nephrol 2009;24(1):67-76 (maturation params)

model_code <- '
$PARAM @annotated
TVCL   : 0.748  : Mature CL at BW=1.76kg SCr=0.4 mg/dL (L/h)
TVV    : 0.866  : V at BW=1.76 kg (L)
BW     : 1.76   : Body weight (kg)
PMA    : 32     : Postmenstrual age (weeks)
SCR    : 0.40   : Serum creatinine (mg/dL)
BWREF  : 1.76   : Reference body weight (kg)
SCRREF : 0.40   : Reference SCr (mg/dL)
TM50   : 47.7   : TM50 for GFR maturation Rhodin 2009 (weeks)
HLCOEF : 3.4    : Hill coefficient maturation Rhodin 2009
THSCR  : -0.635 : SCr exponent on CL
FU     : 0.70   : Unbound fraction of piperacillin

$CMT @annotated
CENT : Central compartment (mg)

$MAIN
double mat  = pow(PMA, HLCOEF) / (pow(PMA, HLCOEF) + pow(TM50, HLCOEF));
double CLi  = TVCL * pow(BW / BWREF, 0.75) * mat * pow(SCR / SCRREF, THSCR);
double Vi   = TVV  * (BW / BWREF);

$ODE
dxdt_CENT = -(CLi / Vi) * CENT;

$TABLE
double CP  = CENT / Vi;    // total plasma concentration (mg/L)
double CPF = CP * FU;      // free (unbound) concentration (mg/L)

$CAPTURE CP CPF CLi Vi mat
'

mod <- mcode("piperacillin_neonatal", model_code)

# ── Theme ───────────────────────────────────────────────────────────────────────
app_theme <- bs_theme(
  version = 5, bootswatch = "flatly", primary = "#8b5cf6"
) |> bs_add_rules("
  .metric-card {
    background: #f8f9fa; border-radius: 8px; padding: 14px 10px;
    text-align: center; border: 1px solid #dee2e6; height: 100%;
  }
  .metric-value { font-size: 22px; font-weight: bold; color: #2c3e50; }
  .metric-label { font-size: 11px; color: #7f8c8d; margin-top: 4px; }
  .metric-sub   { font-size: 10px; margin-top: 2px; }
  .clr-green  .metric-value { color: #10b981; }
  .clr-red    .metric-value { color: #dc3545; }
  .clr-amber  .metric-value { color: #f59e0b; }
  .clr-purple .metric-value { color: #8b5cf6; }
  .clr-teal   .metric-value { color: #0d9488; }
  .pk-note { font-size: 11px; color: #6c757d; background: #f1f3f5;
             border-radius: 4px; padding: 6px 8px; margin-top: 4px; font-family: monospace; }
  .pkpd-footer { text-align: center; font-size: 11px; color: #aaa; padding: 10px 0 4px; }
  .ref-box { background: #f0f4ff; border-left: 4px solid #8b5cf6;
             padding: 12px 16px; border-radius: 4px; margin-top: 8px; font-size: 13px; }
")

# ── UI ──────────────────────────────────────────────────────────────────────────
ui <- page_sidebar(
  title = "Piperacillin Neonatal PopPK Simulator",
  theme = app_theme,

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

    h6("Dosing", style = "font-weight:600; margin-bottom:4px;"),
    sliderInput("dose_pkg", "Dose (mg/kg per dose)",
                min = 30, max = 150, value = 100, step = 10),
    selectInput("interval", "Dosing Interval",
                choices = list("q8h (every 8 h)" = 8, "q12h (every 12 h)" = 12),
                selected = 12),
    selectInput("inf_dur", "Infusion Duration",
                choices = list("30 min (0.5 h)" = 0.5, "1 h" = 1,
                               "2 h" = 2, "3 h" = 3, "4 h" = 4),
                selected = 0.5),
    numericInput("n_days", "Simulation Duration (days)", value = 2, min = 1, max = 3),

    hr(),
    h6("Patient Characteristics", style = "font-weight:600; margin-bottom:4px;"),
    sliderInput("bw",  "Body Weight (kg)", min = 0.5, max = 5.0, value = 1.76, step = 0.1),
    sliderInput("pma", "Postmenstrual Age (weeks)", min = 23, max = 44, value = 32, step = 1),
    sliderInput("scr", "Serum Creatinine (mg/dL)", min = 0.1, max = 1.2, value = 0.40, step = 0.05),
    uiOutput("pk_params_ui"),

    hr(),
    h6("PK/PD Target", style = "font-weight:600; margin-bottom:4px;"),
    selectInput("mic", "MIC Breakpoint",
                choices = list(
                  "4 mg/L"  = 4,
                  "8 mg/L"  = 8,
                  "16 mg/L — Enterobacterales (EUCAST)" = 16,
                  "32 mg/L — Pseudomonas aeruginosa"   = 32
                ),
                selected = 16),
    selectInput("target_pct", "Target %fT>MIC",
                choices = list("≥ 50% (FDA label)" = 50,
                               "≥ 75% (conservative NICU)" = 75),
                selected = 75),

    hr(),
    checkboxInput("show_free", "Show free (unbound) concentration", value = TRUE),
    checkboxInput("log_scale", "Log scale (Y-axis)", value = FALSE)
  ),

  # ── Metric cards ────────────────────────────────────────────────────────────
  layout_columns(
    col_widths = c(3, 3, 3, 3),
    uiOutput("card_cmax"),
    uiOutput("card_ctrough"),
    uiOutput("card_ftmic"),
    uiOutput("card_target")
  ),

  # ── PK plot ─────────────────────────────────────────────────────────────────
  card(
    full_screen = TRUE,
    style = "height: 500px; margin-top: 10px;",
    card_body(padding = 0,
      plotOutput("pkPlot", height = "100%")
    )
  ),

  # ── Info / References tabs ───────────────────────────────────────────────────
  div(style = "margin-top: 10px;",
  navset_card_underline(

    nav_panel("Model Information",
      markdown("
## Piperacillin Neonatal Population PK Model

**Drug:** Piperacillin (administered as piperacillin-tazobactam 8:1 fixed ratio)  
**Population:** Neonates (PNA ≤ 28 days) with late-onset sepsis or severe infection  
**Model structure:** One-compartment, IV infusion, first-order elimination  
**Software:** NONMEM v7.4 (FOCE-I); ADVAN1 TRANS2  

### Parameter Estimates (Table 2)

| Parameter | Estimate | RSE | Description |
|-----------|----------|-----|-------------|
| θ_CL | 0.748 L/h | 8% | CL at full GFR maturation, BW = 1.76 kg, SCr = 0.4 mg/dL |
| θ_SCr | −0.635 | 36% | SCr power-model effect on CL |
| θ_V | 0.866 L | 8% | V at BW = 1.76 kg |
| ω²_CL | 38.3% CV | 12% | Interindividual variability in CL |
| ω²_V | 37.7% CV | 17% | Interindividual variability in V |
| σ_prop | 11.4% | 36% | Proportional residual error |

### Covariate Model

**Clearance:**  
CL_i = 0.748 × (BW_i / 1.76)^0.75 × mat(PMA_i) × (SCr_i / 0.4)^(−0.635)

**Sigmoidal maturation function (Rhodin et al., 2009):**  
mat(PMA) = PMA^3.4 / (PMA^3.4 + 47.7^3.4)  
→ TM50 = 47.7 weeks, Hill = 3.4

**Volume:**  
V_i = 0.866 × (BW_i / 1.76)

### PK/PD Targets

| Target | Threshold | Basis |
|--------|-----------|-------|
| ≥ 50% fT>MIC | FDA label | Standard β-lactam time-dependent killing |
| ≥ 75% fT>MIC | Conservative | Recommended for immunocompromised neonates |
| Free Ctrough < 50 mg/L | Safety threshold | Exploratory nephrotoxicity cutoff |

**MIC breakpoints used in this study:**  
- 16 mg/L → Enterobacterales  
- 32 mg/L → *Pseudomonas aeruginosa*  

**Protein binding:** fu = 0.70 (piperacillin ~30% protein bound)

### Model Validation

- Development: 25 neonates, 65 plasma samples  
- External validation: 8 neonates — MAPE 25.89%, MPE 2.29 mg/L (low bias)  
- Bootstrap (n = 1,000), pcVPC, and NPDE confirmed predictive accuracy
"     )
    ),

    nav_panel("References",
      div(class = "ref-box",
        tags$p(tags$strong("Primary citation:")),
        tags$p("Population pharmacokinetics and dose optimization of piperacillin-tazobactam in neonates with late-onset sepsis and severe infections. 2026."),
        tags$p(tags$strong("Maturation model:")),
        tags$p("Rhodin MM, et al. Human renal function maturation: a quantitative description using weight and postmenstrual age. ",
               tags$em("Pediatric Nephrology"), " 2009;24(1):67-76."),
        tags$p(tags$strong("PK comparison:")),
        tags$p("Cohen-Wolkowiez M, et al. Population pharmacokinetics of piperacillin using scavenged samples from premature infants. ",
               tags$em("Pharmacotherapy"), " 2012;32(8):717-719."),
        tags$p(tags$strong("Kong et al. (2020) pooled piperacillin neonates–adults model.")),
        tags$p(tags$strong("Built with:")),
        tags$p(
          tags$a("PKPDBuilder", href = "https://pkpdbuilder.com", target = "_blank"),
          " — AI-powered pharmacokinetic simulation platform."
        )
      )
    )
  )),

  # ── Footer ──────────────────────────────────────────────────────────────────
  div(class = "pkpd-footer",
      "Simulations based on published neonatal piperacillin PopPK · ",
      tags$a("PKPDBuilder.com", href = "https://pkpdbuilder.com", target = "_blank"),
      " | For research and educational use only")
)

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

  # Compute individual PK parameters ─────────────────────────────────────────
  pk_params <- reactive({
    shiny::req(input$bw, input$pma, input$scr)
    bw  <- input$bw
    pma <- input$pma
    scr <- input$scr
    mat   <- pma^3.4 / (pma^3.4 + 47.7^3.4)
    CLi   <- 0.748 * (bw / 1.76)^0.75 * mat * (scr / 0.4)^(-0.635)
    Vi    <- 0.866 * (bw / 1.76)
    thalf <- Vi / CLi * log(2)
    list(CLi = CLi, Vi = Vi, thalf = thalf, mat = mat)
  })

  # PK parameter display in sidebar ──────────────────────────────────────────
  output$pk_params_ui <- renderUI({
    p <- pk_params()
    div(class = "pk-note",
        sprintf("CL %.3f L/h  |  V %.3f L  |  t½ %.1f h  |  mat %.0f%%",
                p$CLi, p$Vi, p$thalf, p$mat * 100))
  })

  # Run mrgsolve simulation ───────────────────────────────────────────────────
  sim_data <- reactive({
    shiny::req(input$dose_pkg, input$bw, input$pma, input$scr,
               input$interval, input$inf_dur, input$n_days)

    bw       <- input$bw
    dose_mg  <- input$dose_pkg * bw
    interval <- as.numeric(input$interval)
    inf_dur  <- as.numeric(input$inf_dur)
    # Guard: infusion must be shorter than interval
    inf_dur  <- min(inf_dur, interval * 0.9)
    n_days   <- input$n_days
    n_doses  <- ceiling(24 * n_days / interval)

    dose_ev <- ev(
      amt  = dose_mg,
      rate = dose_mg / inf_dur,
      ii   = interval,
      addl = n_doses - 1,
      cmt  = 1
    )

    sim_end <- interval * n_doses

    mod |>
      param(BW = bw, PMA = input$pma, SCR = input$scr) |>
      ev(dose_ev) |>
      mrgsim(end = sim_end, delta = 0.1) |>
      as.data.frame()
  })

  # Compute PK/PD metrics from last dosing interval ──────────────────────────
  metrics <- reactive({
    df       <- sim_data()
    shiny::req(nrow(df) > 0)
    interval <- as.numeric(input$interval)
    mic_val  <- as.numeric(input$mic)
    tgt      <- as.numeric(input$target_pct)

    total_time  <- max(df$time)
    last_start  <- total_time - interval
    last_int    <- dplyr::filter(df, time >= last_start, time <= total_time)

    cmax_free    <- max(last_int$CPF, na.rm = TRUE)
    cmax_total   <- max(last_int$CP,  na.rm = TRUE)
    trough_row   <- dplyr::filter(last_int, time == max(last_int$time))
    ctrough_free <- if (nrow(trough_row) > 0) trough_row$CPF[1] else 0
    ctrough_tot  <- if (nrow(trough_row) > 0) trough_row$CP[1]  else 0

    # %fT>MIC: fraction of interval where free conc exceeds MIC
    n_above    <- sum(last_int$CPF > mic_val, na.rm = TRUE)
    pct_ft_mic <- min(100, n_above * 0.1 / interval * 100)

    target_met  <- pct_ft_mic >= tgt
    safety_risk <- ctrough_free > 50    # exploratory nephrotoxicity threshold

    list(
      cmax_free = cmax_free, cmax_total = cmax_total,
      ctrough_free = ctrough_free, ctrough_tot = ctrough_tot,
      pct_ft_mic = pct_ft_mic,
      target_met = target_met,
      safety_risk = safety_risk
    )
  })

  # ── Metric cards ─────────────────────────────────────────────────────────────
  output$card_cmax <- renderUI({
    m   <- metrics()
    val <- if (input$show_free) m$cmax_free else m$cmax_total
    lbl <- if (input$show_free) "Free Cmax (mg/L)" else "Total Cmax (mg/L)"
    div(class = "metric-card clr-purple",
        div(class = "metric-value", sprintf("%.1f", val)),
        div(class = "metric-label", lbl))
  })

  output$card_ctrough <- renderUI({
    m   <- metrics()
    val <- if (input$show_free) m$ctrough_free else m$ctrough_tot
    lbl <- if (input$show_free) "Free Ctrough (mg/L)" else "Total Ctrough (mg/L)"
    warn <- input$show_free && m$safety_risk
    cls <- if (warn) "metric-card clr-amber" else "metric-card clr-teal"
    div(class = cls,
        div(class = "metric-value", sprintf("%.1f", val)),
        if (warn) div(class = "metric-sub", style = "color:#f59e0b;", "⚠ >50 mg/L (safety limit)"),
        div(class = "metric-label", lbl))
  })

  output$card_ftmic <- renderUI({
    m   <- metrics()
    tgt <- as.numeric(input$target_pct)
    cls <- if (m$target_met) "metric-card clr-green" else "metric-card clr-red"
    div(class = cls,
        div(class = "metric-value", sprintf("%.0f%%", m$pct_ft_mic)),
        div(class = "metric-label", paste0("%fT>MIC (target ≥", tgt, "%)")))
  })

  output$card_target <- renderUI({
    m   <- metrics()
    tgt <- as.numeric(input$target_pct)
    mic <- as.numeric(input$mic)
    if (m$target_met) {
      div(class = "metric-card clr-green",
          div(class = "metric-value", "TARGET MET ✓"),
          div(class = "metric-label", sprintf("≥%d%% fT>MIC at %d mg/L", tgt, mic)))
    } else {
      div(class = "metric-card clr-red",
          div(class = "metric-value", "NOT MET ✗"),
          div(class = "metric-label", sprintf("<%d%% fT>MIC at %d mg/L", tgt, mic)))
    }
  })

  # ── PK Plot ──────────────────────────────────────────────────────────────────
  output$pkPlot <- renderPlot({
    df <- sim_data()
    shiny::req(nrow(df) > 0)

    mic_val <- as.numeric(input$mic)
    fu      <- 0.70

    if (input$show_free) {
      df$CONC     <- df$CPF
      y_lab       <- "Free Piperacillin Concentration (mg/L)"
      mic_line    <- mic_val
      safety_line <- 50
      mic_lbl     <- paste0("MIC = ", mic_val, " mg/L (free drug)")
    } else {
      df$CONC     <- df$CP
      y_lab       <- "Total Piperacillin Concentration (mg/L)"
      mic_line    <- mic_val / fu
      safety_line <- 50 / fu
      mic_lbl     <- paste0("MIC equiv = ", round(mic_val / fu, 1), " mg/L (total)")
    }

    if (input$log_scale) df <- dplyr::filter(df, .data$CONC > 0.001)

    m   <- metrics()
    tgt <- as.numeric(input$target_pct)
    title_col <- if (m$target_met) "#10b981" else "#dc3545"
    title_txt <- sprintf(
      "%s %d%% fT>MIC = %.0f%%  (%s)",
      if (m$target_met) "✓" else "✗",
      tgt, m$pct_ft_mic,
      if (m$target_met) "target met" else "target NOT met"
    )
    subtitle_txt <- sprintf(
      "Dose: %d mg/kg q%sh × %sh infusion  |  BW %.1f kg  PMA %d wk  SCr %.2f mg/dL",
      input$dose_pkg, input$interval, input$inf_dur,
      input$bw, input$pma, input$scr
    )

    p <- ggplot(df, aes(x = .data$time, y = .data$CONC)) +
      # Shaded region below MIC (sub-therapeutic zone)
      annotate("rect",
               xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = mic_line,
               fill = "#dc354515", color = NA) +
      geom_line(color = "#8b5cf6", linewidth = 1.3) +
      # MIC line
      geom_hline(yintercept = mic_line, linetype = "dashed",
                 color = "#dc3545", linewidth = 0.9) +
      annotate("text",
               x = max(df$time) * 0.97, y = mic_line,
               label = mic_lbl, color = "#dc3545",
               size = 3.5, vjust = -0.6, hjust = 1) +
      # Safety threshold line
      geom_hline(yintercept = safety_line, linetype = "dotted",
                 color = "#f59e0b", linewidth = 0.8) +
      annotate("text",
               x = max(df$time) * 0.97, y = safety_line,
               label = "Safety threshold\n(free Ctrough = 50 mg/L)",
               color = "#f59e0b", size = 3.0, vjust = -0.2, hjust = 1) +
      labs(
        x = "Time (h)",
        y = y_lab,
        title    = title_txt,
        subtitle = subtitle_txt
      ) +
      theme_minimal(base_size = 13) +
      theme(
        plot.title    = element_text(face = "bold", size = 14, color = title_col),
        plot.subtitle = element_text(color = "#555", size = 11),
        panel.grid.minor = element_blank(),
        plot.margin = margin(12, 16, 8, 12)
      )

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

    p
  })
}

shinyApp(ui, server)