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

# ── Vancomycin PopPK Model ──
# Goti V, Chaturvedula A, Fossler MJ, Mok S, Jacob JT.
# Ther Drug Monit. 2018;40:212-221.
# Two-compartment IV infusion model with creatinine clearance and dialysis covariates.

model_code <- '
$PARAM @annotated
CL    : 4.5   : Clearance for non-dialysis subject at CrCL 120 mL/min (L/h)
VC    : 58.4  : Central volume for 70-kg non-dialysis subject (L)
VP    : 38.4  : Peripheral volume for 70-kg subject (L)
Q     : 6.5   : Intercompartmental clearance (L/h)
CRCL  : 62    : Creatinine clearance (mL/min)
WT    : 79    : Body weight (kg)
DIAL  : 0     : Hemodialysis status (0 = no, 1 = yes)

$CMT @annotated
CENT   : Central compartment (mg)
PERIPH : Peripheral compartment (mg)

$MAIN
// Final model from Table 2:
// TVCL = theta1 * (CrCL/120)^theta2 * theta3^DIAL
// TVVc = theta4 * (WT/70) * theta5^DIAL

double CLi = CL * pow(CRCL / 120.0, 0.8) * pow(0.7, DIAL);
double VCi = VC * (WT / 70.0) * pow(0.5, DIAL);
double VPi = VP * (WT / 70.0);
double Qi  = Q;

// Safety floors
CLi = CLi > 0.05 ? CLi : 0.05;
VCi = VCi > 1.0 ? VCi : 1.0;
VPi = VPi > 1.0 ? VPi : 1.0;

$ODE
dxdt_CENT   = - (CLi / VCi + Qi / VCi) * CENT + (Qi / VPi) * PERIPH;
dxdt_PERIPH =   (Qi / VCi) * CENT - (Qi / VPi) * PERIPH;

$TABLE
double CP = CENT / VCi;
CP = CP > 0 ? CP : 0;

$CAPTURE @annotated
CP  : Plasma vancomycin concentration (mg/L)
CLi : Individual clearance (L/h)
VCi : Individual central volume (L)
VPi : Individual peripheral volume (L)
Qi  : Individual intercompartmental clearance (L/h)
'

mod <- mcode("vancomycin_goti_2018", model_code)

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; }
  .ref-box {
    background: #f0f4ff;
    border-left: 4px solid #8b5cf6;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
  .nomogram-box {
    background: #f8fafc;
    border: 1px solid #e2e8f0;
    border-radius: 8px;
    padding: 12px;
    margin-top: 10px;
    font-size: 13px;
  }
  .note-box {
    background: #fff7ed;
    border-left: 4px solid #f59e0b;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
")

recommend_nomogram <- function(crcl, dialysis) {
  if (isTRUE(dialysis)) {
    return(list(
      loading_mgkg = 15,
      maint_mgkg_day = 6.7,
      crcl_band = "Hemodialysis",
      source = "Paper optimized nomogram"
    ))
  }

  if (crcl <= 30) {
    list(loading_mgkg = 15, maint_mgkg_day = 6.7, crcl_band = "≤30 mL/min", source = "Paper optimized nomogram")
  } else if (crcl <= 60) {
    list(loading_mgkg = 25, maint_mgkg_day = 12.9, crcl_band = "31-60 mL/min", source = "Paper optimized nomogram")
  } else if (crcl <= 90) {
    list(loading_mgkg = 25, maint_mgkg_day = 23.9, crcl_band = "60-90 mL/min", source = "Paper optimized nomogram")
  } else if (crcl <= 120) {
    list(loading_mgkg = 25, maint_mgkg_day = 38.2, crcl_band = "91-120 mL/min", source = "Paper optimized nomogram")
  } else if (crcl <= 150) {
    list(loading_mgkg = 25, maint_mgkg_day = 37.8, crcl_band = "121-150 mL/min", source = "Paper optimized nomogram")
  } else {
    list(loading_mgkg = 25, maint_mgkg_day = 37.7, crcl_band = ">150 mL/min", source = "Paper optimized nomogram")
  }
}

build_regimen <- function(weight, loading_mgkg, maint_mgkg_day, interval_h, infusion_h, n_days) {
  maint_amt <- maint_mgkg_day * weight * interval_h / 24
  load_amt <- loading_mgkg * weight
  n_maint <- max(ceiling(n_days * 24 / interval_h) - 1, 0)

  events <- list()

  if (load_amt > 0) {
    events[[length(events) + 1]] <- ev(
      time = 0,
      amt = load_amt,
      cmt = 1,
      rate = load_amt / infusion_h
    )
  }

  if (maint_amt > 0 && n_maint > 0) {
    events[[length(events) + 1]] <- ev(
      time = interval_h,
      amt = maint_amt,
      cmt = 1,
      rate = maint_amt / infusion_h,
      ii = interval_h,
      addl = n_maint - 1
    )
  }

  if (length(events) == 0) {
    return(NULL)
  }

  do.call(c, events)
}

calc_auc <- function(time, conc) {
  if (length(time) < 2 || length(conc) < 2) {
    return(NA_real_)
  }
  sum(diff(time) * (head(conc, -1) + tail(conc, -1)) / 2)
}

calc_terminal_half_life <- function(cl, vc, vp, q) {
  k10 <- cl / vc
  k12 <- q / vc
  k21 <- q / vp
  trace_term <- k10 + k12 + k21
  disc <- trace_term^2 - 4 * k10 * k21
  disc <- ifelse(disc < 0, 0, disc)
  beta <- (trace_term - sqrt(disc)) / 2
  if (beta <= 0) {
    return(NA_real_)
  }
  log(2) / beta
}

paper_performance <- data.frame(
  dose_number = factor(
    c("After first dose", "After second dose", "After third dose", "After fourth dose", "After fifth dose", "After sixth dose"),
    levels = c("After first dose", "After second dose", "After third dose", "After fourth dose", "After fifth dose", "After sixth dose")
  ),
  Current = c(27.5, 44.5, 43.6, 39.2, 35.0, 31.4),
  Optimized = c(46.7, 48.3, 48.3, 47.1, 45.4, 43.7)
)

ui <- page_sidebar(
  title = "Vancomycin PopPK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,

    h6("Patient Characteristics"),
    sliderInput("weight", "Body Weight (kg)", min = 40, max = 180, value = 79, step = 1),
    sliderInput("crcl", "Creatinine Clearance (mL/min)", min = 5, max = 180, value = 62, step = 1),
    checkboxInput("dialysis", "Intermittent Hemodialysis", value = FALSE),

    hr(),
    h6("Dosing Regimen"),
    numericInput("loading_mgkg", "Loading Dose (mg/kg)", value = 25, min = 0, max = 30, step = 1),
    numericInput("maint_mgkg_day", "Maintenance Dose (mg/kg/day)", value = 23.9, min = 0, max = 80, step = 0.1),
    sliderInput("interval", "Dosing Interval (h)", min = 8, max = 48, value = 12, step = 4),
    sliderInput("infusion_h", "Infusion Duration (h)", min = 0.5, max = 4, value = 1.5, step = 0.5),
    numericInput("n_days", "Simulation Duration (days)", value = 5, min = 2, max = 14),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE),

    hr(),
    actionButton("apply_nomogram", "Apply Optimized Nomogram", class = "btn-primary"),
    uiOutput("recommended_regimen"),

    div(class = "note-box",
      strong("Important limitation:"),
      tags$br(),
      "The paper modeled dialysis as a binary covariate on CL and Vc. Actual dialysis session timing and extracorporeal removal were not parameterized, so this simulator should not be used for patient-specific dialysis dosing."
    )
  ),

  layout_column_wrap(
    width = 1 / 4,
    fill = FALSE,
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("cmax")),
      div(class = "metric-label", "Cmax (mg/L)")
    ),
    div(class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("ctrough")),
      div(class = "metric-label", "Ctrough (mg/L)")
    ),
    div(class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("auc")),
      div(class = "metric-label", "AUCτ (mg·h/L)")
    ),
    div(class = "metric-card metric-info",
      div(class = "metric-value", textOutput("thalf")),
      div(class = "metric-label", "Terminal t½ (h)")
    )
  ),

  card(
    card_header("Vancomycin Concentration-Time Profile"),
    full_screen = TRUE,
    plotOutput("pkPlot", height = "500px")
  ),

  navset_card_underline(
    title = "Drug Information",
    nav_panel(
      "Model Information",
      markdown(
        paste0(
          "## Vancomycin Population PK Model\n\n",
          "**Source:** Goti V, Chaturvedula A, Fossler MJ, Mok S, Jacob JT. *Hospitalized Patients With and Without Hemodialysis Have Markedly Different Vancomycin Pharmacokinetics: A Population Pharmacokinetic Model-Based Analysis.* Ther Drug Monit. 2018;40:212-221.\n\n",
          "### Model structure\n",
          "- **Structural model:** 2-compartment IV infusion model\n",
          "- **Estimation:** NONMEM FOCEI\n",
          "- **Population:** 1812 hospitalized adults, including 336 on intermittent hemodialysis\n",
          "- **Residual error:** combined additive + proportional\n\n",
          "### Final parameter estimates\n",
          "- **CL:** 4.5 L/h\n",
          "- **Vc:** 58.4 L\n",
          "- **Vp:** 38.4 L\n",
          "- **Q:** 6.5 L/h\n",
          "- **CrCL effect on CL:** exponent 0.8\n",
          "- **Dialysis effect on CL:** 0.7-fold\n",
          "- **Dialysis effect on Vc:** 0.5-fold\n\n",
          "### Covariate relationships\n",
          "- **TVCL = 4.5 × (CrCL/120)^0.8 × 0.7^DIAL**\n",
          "- **TVVc = 58.4 × (WT/70) × 0.5^DIAL**\n",
          "- **Vp** is weight-normalized here as **38.4 × WT/70** to preserve volume scaling used in the manuscript context\n\n",
          "### Clinical context\n",
          "- Therapeutic trough window discussed in the paper: **10-20 mg/L**\n",
          "- Hemodialysis status materially lowers both clearance and central volume\n",
          "- The paper's optimized nomogram improved therapeutic-range attainment after the first six doses\n"
        )
      )
    ),
    nav_panel(
      "Nomogram Performance",
      plotOutput("nomogramPlot", height = "350px"),
      div(class = "nomogram-box",
        tags$strong("Published optimized nomogram (Table 4):"),
        tags$ul(
          tags$li("CrCL ≤30 mL/min: 6.7 mg/kg/day, loading 15 mg/kg"),
          tags$li("CrCL 31-60 mL/min: 12.9 mg/kg/day, loading 25 mg/kg"),
          tags$li("CrCL 60-90 mL/min: 23.9 mg/kg/day, loading 25 mg/kg"),
          tags$li("CrCL 91-120 mL/min: 38.2 mg/kg/day, loading 25 mg/kg"),
          tags$li("CrCL 121-150 mL/min: 37.8 mg/kg/day, loading 25 mg/kg"),
          tags$li("CrCL >150 mL/min: 37.7 mg/kg/day, loading 25 mg/kg"),
          tags$li("Hemodialysis patients: mean maintenance dose 6.7 mg/kg/day, loading 15 mg/kg")
        )
      )
    ),
    nav_panel(
      "References",
      div(class = "ref-box",
        tags$h5("Key references"),
        tags$ol(
          tags$li("Goti V, Chaturvedula A, Fossler MJ, Mok S, Jacob JT. Ther Drug Monit. 2018;40:212-221."),
          tags$li("Rybak MJ et al. Vancomycin therapeutic guidelines. Clin Infect Dis. 2009;49:325-327."),
          tags$li("Marsot A et al. Vancomycin: a review of population pharmacokinetic analyses. Clin Pharmacokinet. 2012;51:1-13.")
        ),
        tags$h5("Drug context"),
        tags$ul(
          tags$li(tags$strong("Drug class:"), " Glycopeptide antibiotic"),
          tags$li(tags$strong("Route modeled:"), " Intravenous infusion"),
          tags$li(tags$strong("Primary elimination:"), " Renal"),
          tags$li(tags$strong("Typical monitoring target in the paper:"), " Trough 10-20 mg/L")
        )
      )
    )
  ),

  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.")
  )
)

server <- function(input, output, session) {
  observeEvent(input$apply_nomogram, {
    reco <- recommend_nomogram(input$crcl, input$dialysis)
    updateNumericInput(session, "loading_mgkg", value = reco$loading_mgkg)
    updateNumericInput(session, "maint_mgkg_day", value = reco$maint_mgkg_day)
  })

  output$recommended_regimen <- renderUI({
    reco <- recommend_nomogram(input$crcl, input$dialysis)
    dose_amt <- reco$maint_mgkg_day * input$weight * input$interval / 24
    div(
      class = "nomogram-box",
      strong("Optimized paper-based recommendation"),
      tags$br(),
      paste0("Band: ", reco$crcl_band),
      tags$br(),
      paste0("Loading dose: ", sprintf("%.1f", reco$loading_mgkg), " mg/kg"),
      tags$br(),
      paste0("Maintenance: ", sprintf("%.1f", reco$maint_mgkg_day), " mg/kg/day"),
      tags$br(),
      paste0("Equivalent dose at current interval: ", sprintf("%.0f", dose_amt), " mg every ", input$interval, " h")
    )
  })

  sim_data <- reactive({
    shiny::req(input$weight, input$crcl, input$interval, input$infusion_h, input$n_days)

    reg <- build_regimen(
      weight = input$weight,
      loading_mgkg = input$loading_mgkg,
      maint_mgkg_day = input$maint_mgkg_day,
      interval_h = input$interval,
      infusion_h = input$infusion_h,
      n_days = input$n_days
    )

    shiny::req(!is.null(reg))

    out <- mod %>%
      param(
        WT = input$weight,
        CRCL = input$crcl,
        DIAL = as.numeric(input$dialysis)
      ) %>%
      ev(reg) %>%
      mrgsim(end = input$n_days * 24, delta = 0.1) %>%
      as.data.frame()

    out <- out |>
      mutate(
        time_h = time,
        dose_interval = input$interval,
        last_start = max(time_h) - input$interval
      )

    out
  })

  last_interval <- reactive({
    d <- sim_data()
    end_time <- max(d$time_h)
    start_time <- max(end_time - input$interval, 0)
    d |>
      dplyr::filter(time_h >= start_time, time_h <= end_time)
  })

  output$cmax <- renderText({
    d <- last_interval()
    sprintf("%.1f", max(d$CP, na.rm = TRUE))
  })

  output$ctrough <- renderText({
    d <- last_interval()
    sprintf("%.1f", d$CP[nrow(d)])
  })

  output$auc <- renderText({
    d <- last_interval()
    auc_val <- calc_auc(d$time_h, d$CP)
    sprintf("%.1f", auc_val)
  })

  output$thalf <- renderText({
    d <- sim_data()
    th <- calc_terminal_half_life(
      cl = d$CLi[1],
      vc = d$VCi[1],
      vp = d$VPi[1],
      q = d$Qi[1]
    )
    sprintf("%.1f", th)
  })

  output$pkPlot <- renderPlot({
    d <- sim_data()
    d_plot <- if (input$log_scale) d |> dplyr::filter(CP > 0.001) else d

    p <- ggplot(d_plot, aes(x = time_h, y = CP)) +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 10, ymax = 20, fill = "#10b981", alpha = 0.10) +
      geom_hline(yintercept = c(10, 20), linetype = "dashed", color = "#10b981", alpha = 0.7) +
      geom_line(color = "#8b5cf6", linewidth = 0.9) +
      labs(
        x = "Time (hours)",
        y = "Vancomycin Concentration (mg/L)",
        title = paste0(
          "Vancomycin ",
          sprintf("%.0f", input$loading_mgkg * input$weight),
          " mg load, then ",
          sprintf("%.0f", input$maint_mgkg_day * input$weight * input$interval / 24),
          " mg q", input$interval, "h"
        ),
        subtitle = "Green band shows the 10-20 mg/L trough window discussed in the paper"
      ) +
      theme_minimal(base_size = 14)

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

    p
  })

  output$nomogramPlot <- renderPlot({
    d <- data.frame(
      dose_number = rep(paper_performance$dose_number, 2),
      Percent = c(paper_performance$Current, paper_performance$Optimized),
      Regimen = rep(c("Current nomogram", "Optimized nomogram"), each = nrow(paper_performance))
    )

    ggplot(d, aes(x = dose_number, y = Percent, color = Regimen, group = Regimen)) +
      geom_line(linewidth = 1) +
      geom_point(size = 2.5) +
      scale_color_manual(values = c("Current nomogram" = "#94a3b8", "Optimized nomogram" = "#8b5cf6")) +
      labs(
        x = NULL,
        y = "% in therapeutic range",
        color = NULL,
        title = "Published therapeutic-range attainment across the first six troughs"
      ) +
      theme_minimal(base_size = 13) +
      theme(
        legend.position = "top",
        axis.text.x = element_text(angle = 20, hjust = 1)
      )
  })
}

shinyApp(ui = ui, server = server)