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

model_code <- '
$PARAM @annotated
V1      : 19.7 : Central volume of distribution (L)
V2      : 15.0 : Peripheral volume of distribution (L)
Q       : 52.4 : Intercompartmental clearance (L/h)
KM      : 8.04 : Michaelis-Menten constant (mg/L)
TVVMAX  : 1900 : Typical Vmax at BSA 1.94 m2 (mg/h)
EXPO    : 1.14 : BSA exponent on Vmax
BSA     : 1.94 : Body surface area (m2)
ETA_VM  : 0    : Inter-individual variability on Vmax

$CMT @annotated
CENT   : Central compartment (mg)
PERIPH : Peripheral compartment (mg)
AUC    : Area under concentration-time curve (mg*h/L)

$MAIN
double VMAXi = TVVMAX * pow(BSA / 1.94, EXPO) * exp(ETA_VM);
double K12 = Q / V1;
double K21 = Q / V2;

$ODE
double CP = CENT / V1;
double ELIM = (CP > 0) ? (CP * VMAXi / (KM + CP)) : 0;
dxdt_CENT = -ELIM + K21 * PERIPH - K12 * CENT;
dxdt_PERIPH = K12 * CENT - K21 * PERIPH;
dxdt_AUC = CP;

$TABLE
double CP_MG_L = CENT / V1;
double AUCOUT = AUC;

$CAPTURE CP_MG_L AUCOUT VMAXi
'

mod <- mcode("pkpd_app_347_5fu", 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;
    min-height: 108px;
  }
  .metric-value {
    font-size: 24px;
    font-weight: bold;
    color: #2c3e50;
  }
  .metric-label {
    font-size: 12px;
    color: #7f8c8d;
    margin-top: 4px;
  }
  .metric-subtext {
    font-size: 11px;
    color: #6c757d;
    margin-top: 6px;
  }
  .metric-success .metric-value { color: #10b981; }
  .metric-warning .metric-value { color: #f59e0b; }
  .metric-primary .metric-value { color: #8b5cf6; }
  .metric-info .metric-value { color: #0dcaf0; }
  .info-box {
    background: #f8f9ff;
    border: 1px solid #d8ddff;
    border-radius: 8px;
    padding: 12px 14px;
    margin-top: 10px;
    font-size: 13px;
  }
  .ref-box {
    background: #f0f4ff;
    border-left: 4px solid #8b5cf6;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
  .table-note {
    color: #6c757d;
    font-size: 12px;
    margin-top: 8px;
  }
")

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

calc_terminal_half_life_min <- function(bsa) {
  vmax <- 1900 * (bsa / 1.94)^1.14
  cl_lin <- vmax / 8.04
  k10 <- cl_lin / 19.7
  k12 <- 52.4 / 19.7
  k21 <- 52.4 / 15.0
  beta <- 0.5 * ((k10 + k12 + k21) - sqrt((k10 + k12 + k21)^2 - 4 * k21 * k10))
  log(2) / beta * 60
}

regimen_defaults <- function(preset) {
  switch(
    preset,
    "Figure 5A: 400 + 2400 mg/m² over 46 h" = list(
      regimen = "bolus_plus_infusion",
      bolus_m2 = 400,
      infusion_m2 = 2400,
      infusion_h = 46
    ),
    "Figure 5B: 200 + 2400 mg/m² over 46 h" = list(
      regimen = "bolus_plus_infusion",
      bolus_m2 = 200,
      infusion_m2 = 2400,
      infusion_h = 46
    ),
    "Figure 5C: 2400 mg/m² over 24 h" = list(
      regimen = "infusion_only",
      bolus_m2 = 0,
      infusion_m2 = 2400,
      infusion_h = 24
    ),
    "Figure 5D: 2800 mg/m² over 24 h" = list(
      regimen = "infusion_only",
      bolus_m2 = 0,
      infusion_m2 = 2800,
      infusion_h = 24
    ),
    list(
      regimen = "bolus_plus_infusion",
      bolus_m2 = 400,
      infusion_m2 = 2400,
      infusion_h = 46
    )
  )
}

build_events <- function(regimen, bsa, bolus_m2, infusion_m2, infusion_h, bolus_h = 0.1) {
  event_rows <- list()

  if (regimen == "bolus_plus_infusion" && bolus_m2 > 0) {
    bolus_amt <- bolus_m2 * bsa
    event_rows[[length(event_rows) + 1]] <- data.frame(
      time = 0,
      amt = bolus_amt,
      cmt = 1,
      evid = 1,
      rate = bolus_amt / bolus_h
    )
  }

  if (infusion_m2 > 0 && infusion_h > 0) {
    infusion_amt <- infusion_m2 * bsa
    infusion_start <- if (regimen == "bolus_plus_infusion" && bolus_m2 > 0) bolus_h else 0
    event_rows[[length(event_rows) + 1]] <- data.frame(
      time = infusion_start,
      amt = infusion_amt,
      cmt = 1,
      evid = 1,
      rate = infusion_amt / infusion_h
    )
  }

  bind_rows(event_rows) |>
    arrange(time)
}

simulate_population <- function(
  regimen,
  bsa,
  bolus_m2,
  infusion_m2,
  infusion_h,
  sim_h,
  n_subjects = 200,
  delta = 0.05
) {
  eta_sd <- sqrt(0.10)
  idata <- data.frame(
    ID = seq_len(n_subjects),
    ETA_VM = rnorm(n_subjects, mean = 0, sd = eta_sd),
    BSA = bsa
  )
  events <- build_events(regimen, bsa, bolus_m2, infusion_m2, infusion_h)
  events_all <- merge(data.frame(ID = idata$ID), events, by = NULL) |>
    arrange(ID, time)

  sim <- mod %>%
    idata_set(idata) %>%
    data_set(events_all) %>%
    mrgsim(end = sim_h, delta = delta) %>%
    as.data.frame() |>
    mutate(CP_MG_L = pmax(CP_MG_L, 1e-6))

  plot_data <- sim |>
    summarise(
      median_cp = median(CP_MG_L, na.rm = TRUE),
      p05_cp = quantile(CP_MG_L, probs = 0.05, na.rm = TRUE),
      p95_cp = quantile(CP_MG_L, probs = 0.95, na.rm = TRUE),
      .by = time
    )

  auc_data <- sim |>
    summarise(
      auc = trap_auc(time, CP_MG_L),
      cmax = max(CP_MG_L, na.rm = TRUE),
      trough = dplyr::last(CP_MG_L),
      .by = ID
    )

  list(
    plot_data = plot_data,
    auc_data = auc_data
  )
}

simulate_typical <- function(regimen, bsa, bolus_m2, infusion_m2, infusion_h, sim_h, delta = 0.05) {
  events <- build_events(regimen, bsa, bolus_m2, infusion_m2, infusion_h)

  mod %>%
    param(BSA = bsa, ETA_VM = 0) %>%
    ev(events) %>%
    mrgsim(end = sim_h, delta = delta) %>%
    as.data.frame() |>
    mutate(CP_MG_L = pmax(CP_MG_L, 1e-6))
}

ui <- page_sidebar(
  title = "5-Fluorouracil Population PK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    selectInput(
      "preset",
      "Paper Preset",
      choices = c(
        "Custom",
        "Figure 5A: 400 + 2400 mg/m² over 46 h",
        "Figure 5B: 200 + 2400 mg/m² over 46 h",
        "Figure 5C: 2400 mg/m² over 24 h",
        "Figure 5D: 2800 mg/m² over 24 h"
      ),
      selected = "Figure 5A: 400 + 2400 mg/m² over 46 h"
    ),
    radioButtons(
      "regimen",
      "Regimen Type",
      choices = c(
        "Bolus + continuous infusion" = "bolus_plus_infusion",
        "Continuous infusion only" = "infusion_only"
      ),
      selected = "bolus_plus_infusion"
    ),
    numericInput("bsa", "Body surface area (m²)", value = 1.94, min = 1.2, max = 2.8, step = 0.01),
    numericInput("bolus_m2", "Bolus dose (mg/m²)", value = 400, min = 0, max = 800, step = 50),
    numericInput("infusion_m2", "Infusion dose (mg/m²)", value = 2400, min = 200, max = 4000, step = 100),
    numericInput("infusion_h", "Infusion duration (h)", value = 46, min = 4, max = 72, step = 1),
    numericInput("sim_h", "Simulation horizon (h)", value = 60, min = 8, max = 120, step = 1),
    checkboxInput("log_y", "Log concentration scale", value = FALSE),
    div(
      class = "info-box",
      tags$strong("Paper targets"),
      tags$ul(
        tags$li("Target AUC0-inf: 20-30 mg*h/L"),
        tags$li("Best paper-supported starting regimens: 200 + 2400 mg/m² over 46 h, or 2800 mg/m² over 24 h"),
        tags$li("Limited sampling: 0.1, 1, 3 h for bolus + infusion; 0.8, 2, 5 h for continuous infusion")
      )
    )
  ),
  navset_card_underline(
    title = "5-Fluorouracil PopPK Explorer",
    nav_panel(
      "Simulation",
      layout_column_wrap(
        width = 1 / 4,
        fill = FALSE,
        div(
          class = "metric-card metric-success",
          div(class = "metric-value", textOutput("cmax")),
          div(class = "metric-label", "Median Cmax (mg/L)"),
          div(class = "metric-subtext", textOutput("cmax_pi"))
        ),
        div(
          class = "metric-card metric-warning",
          div(class = "metric-value", textOutput("ctrough")),
          div(class = "metric-label", "Median end-of-cycle C (mg/L)"),
          div(class = "metric-subtext", "At the selected simulation horizon")
        ),
        div(
          class = "metric-card metric-primary",
          div(class = "metric-value", textOutput("auc")),
          div(class = "metric-label", "Median AUC0-inf surrogate (mg*h/L)"),
          div(class = "metric-subtext", textOutput("auc_pi"))
        ),
        div(
          class = "metric-card metric-info",
          div(class = "metric-value", textOutput("thalf")),
          div(class = "metric-label", "Approx terminal t1/2 (min)"),
          div(class = "metric-subtext", "Low-concentration linearized approximation")
        )
      ),
      card(
        full_screen = TRUE,
        card_header("Population prediction interval"),
        plotOutput("pk_plot", height = "500px")
      ),
      card(
        card_header("Paper-aligned regimen summary"),
        tableOutput("regimen_summary"),
        div(
          class = "table-note",
          "The ribbon and summary table are based on 200 virtual subjects using the reported IIV on Vmax (omega^2 = 0.10)."
        )
      )
    ),
    nav_panel(
      "Model Information",
      markdown(
        "
## Final 5-FU model from Tan et al. 2026

- Structure: two-compartment IV model with Michaelis-Menten elimination
- Covariate: body surface area on Vmax only
- Central volume (V1): 19.7 L
- Peripheral volume (V2): 15.0 L
- Q: 52.4 L/h
- Km: 8.04 mg/L
- Typical Vmax: 1900 mg/h at BSA 1.94 m²
- BSA exponent on Vmax: 1.14
- Reported proportional residual error: 34%

### Clinical context
- Standard FOLFOX-like schedule in the paper: 400 mg/m² bolus followed by 2400 mg/m² over 46 h
- Paper conclusion: 200 + 2400 mg/m² over 46 h and 2800 mg/m² over 24 h improved median AUC toward the 20-30 mg*h/L target window
- Limited sampling strategies:
  - Bolus + infusion: 0.1, 1, and 3 h
  - Continuous infusion only: 0.8, 2, and 5 h

### Notes
- This app is a population PK simulator, not a certified clinical dosing tool.
- The displayed AUC uses a long enough single-cycle simulation horizon to approximate AUC0-inf for the selected regimen.
        "
      )
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Primary paper"),
        tags$ol(
          tags$li(
            "Tan Z, Sancho-Araiz A, Voller S, et al. Advancing Precision Dosing of 5-FU: Population PK Model Development, Limited Sampling Strategies, and Fit-for-Use Application. Clin Pharmacokinet. 2026. doi:10.1007/s40262-026-01645-1."
          )
        ),
        tags$h5("Drug context"),
        tags$ul(
          tags$li(tags$strong("Drug class:"), " Antimetabolite fluoropyrimidine"),
          tags$li(tags$strong("Primary route here:"), " Intravenous bolus plus prolonged infusion, or infusion only"),
          tags$li(tags$strong("Typical therapeutic target:"), " AUC0-inf 20-30 mg*h/L"),
          tags$li(tags$strong("Common regimens in the paper:"), " FOLFOX, FOLFIRI, FOLFIRINOX, FLOT")
        )
      )
    )
  ),
  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", 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$preset, {
    if (input$preset == "Custom") {
      return()
    }
    defaults <- regimen_defaults(input$preset)
    updateRadioButtons(session, "regimen", selected = defaults$regimen)
    updateNumericInput(session, "bolus_m2", value = defaults$bolus_m2)
    updateNumericInput(session, "infusion_m2", value = defaults$infusion_m2)
    updateNumericInput(session, "infusion_h", value = defaults$infusion_h)
    updateNumericInput(
      session,
      "sim_h",
      value = max(defaults$infusion_h + 12, 36)
    )
  }, ignoreInit = FALSE)

  sim_bundle <- reactive({
    shiny::req(input$bsa, input$infusion_m2, input$infusion_h, input$sim_h)

    population <- simulate_population(
      regimen = input$regimen,
      bsa = input$bsa,
      bolus_m2 = input$bolus_m2,
      infusion_m2 = input$infusion_m2,
      infusion_h = input$infusion_h,
      sim_h = input$sim_h
    )

    typical <- simulate_typical(
      regimen = input$regimen,
      bsa = input$bsa,
      bolus_m2 = input$bolus_m2,
      infusion_m2 = input$infusion_m2,
      infusion_h = input$infusion_h,
      sim_h = input$sim_h
    )

    list(population = population, typical = typical)
  })

  output$cmax <- renderText({
    stats <- sim_bundle()$population$auc_data
    sprintf("%.2f", median(stats$cmax, na.rm = TRUE))
  })

  output$cmax_pi <- renderText({
    stats <- sim_bundle()$population$auc_data
    sprintf(
      "90%% PI %.2f-%.2f",
      quantile(stats$cmax, probs = 0.05, na.rm = TRUE),
      quantile(stats$cmax, probs = 0.95, na.rm = TRUE)
    )
  })

  output$ctrough <- renderText({
    stats <- sim_bundle()$population$auc_data
    sprintf("%.4f", median(stats$trough, na.rm = TRUE))
  })

  output$auc <- renderText({
    stats <- sim_bundle()$population$auc_data
    sprintf("%.1f", median(stats$auc, na.rm = TRUE))
  })

  output$auc_pi <- renderText({
    stats <- sim_bundle()$population$auc_data
    sprintf(
      "90%% PI %.1f-%.1f",
      quantile(stats$auc, probs = 0.05, na.rm = TRUE),
      quantile(stats$auc, probs = 0.95, na.rm = TRUE)
    )
  })

  output$thalf <- renderText({
    sprintf("%.1f", calc_terminal_half_life_min(input$bsa))
  })

  output$pk_plot <- renderPlot({
    d <- sim_bundle()$population$plot_data
    shiny::req(nrow(d) > 0)

    if (isTRUE(input$log_y)) {
      d <- dplyr::filter(d, median_cp > 0, p05_cp > 0, p95_cp > 0)
    }

    p <- ggplot(d, aes(x = time, y = median_cp)) +
      geom_ribbon(aes(ymin = p05_cp, ymax = p95_cp), fill = "#8b5cf6", alpha = 0.18) +
      geom_line(color = "#8b5cf6", linewidth = 1) +
      geom_vline(
        xintercept = if (input$regimen == "bolus_plus_infusion" && input$bolus_m2 > 0) 0.1 else 0,
        linetype = "dotted",
        color = "#6c757d"
      ) +
      geom_vline(
        xintercept = if (input$regimen == "bolus_plus_infusion" && input$bolus_m2 > 0) 0.1 + input$infusion_h else input$infusion_h,
        linetype = "dashed",
        color = "#10b981"
      ) +
      labs(
        x = "Time (hours)",
        y = "5-FU concentration (mg/L)",
        title = "Population concentration-time profile",
        subtitle = "Solid line = median; ribbon = 90% prediction interval"
      ) +
      theme_minimal(base_size = 14)

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

    p
  })

  output$regimen_summary <- renderTable({
    stats <- sim_bundle()$population$auc_data
    data.frame(
      regimen = if (input$regimen == "bolus_plus_infusion") "Bolus + infusion" else "Infusion only",
      bsa_m2 = sprintf("%.2f", input$bsa),
      bolus_mg_m2 = sprintf("%.0f", input$bolus_m2),
      infusion_mg_m2 = sprintf("%.0f", input$infusion_m2),
      infusion_h = sprintf("%.0f", input$infusion_h),
      median_auc_mg_h_L = sprintf("%.1f", median(stats$auc, na.rm = TRUE)),
      auc_90_pi = sprintf(
        "%.1f-%.1f",
        quantile(stats$auc, probs = 0.05, na.rm = TRUE),
        quantile(stats$auc, probs = 0.95, na.rm = TRUE)
      ),
      check.names = FALSE
    )
  }, striped = TRUE, spacing = "s")
}

shinyApp(ui = ui, server = server)