use("shiny")
use("bslib")
use("mrgsolve")
use("dplyr")
use("ggplot2")
use("tidyr")

model_code <- '
$PARAM @annotated
CL : 10   : Clearance (L/h)
V1 : 10   : Central volume (L)
Q  : 5    : Intercompartmental clearance (L/h)
V2 : 20   : Peripheral volume (L)
KA : 1    : First-order absorption rate constant (1/h)

$CMT @annotated
DEPOT : Extravascular depot (mg)
CENT  : Central compartment (mg)
PERI  : Peripheral compartment (mg)

$ODE
double k10 = CL / V1;
double k12 = Q / V1;
double k21 = Q / V2;

dxdt_DEPOT = -KA * DEPOT;
dxdt_CENT  = KA * DEPOT - k10 * CENT - k12 * CENT + k21 * PERI;
dxdt_PERI  = k12 * CENT - k21 * PERI;

$TABLE
double CP = CENT / V1;

$CAPTURE CP
'

mod <- mrgsolve::mcode("tmp_sulfonamides_pig_355", model_code)

metric_card <- function(output_id, label, class = "metric-primary") {
  div(
    class = paste("metric-card", class),
    div(class = "metric-value", textOutput(output_id)),
    div(class = "metric-label", label)
  )
}

combo_catalog <- data.frame(
  combo_id = c("TMP_SDZ", "TMP_SMX", "TMP_SDMX"),
  label = c("TMP/SDZ", "TMP/SMX", "TMP/SDMX"),
  sulfonamide = c("SDZ", "SMX", "SDMX"),
  iv_tmp = c(2.5, 6.0, 4.0),
  iv_s = c(12.5, 30.0, 18.6),
  oral_tmp = c(5.0, 6.0, 8.0),
  oral_s = c(25.0, 30.0, 37.36),
  im_tmp = c(2.5, NA, NA),
  im_s = c(12.5, NA, NA),
  default_end_h = c(56, 56, 96),
  stringsAsFactors = FALSE
)

analyte_parameters <- list(
  TMP = list(
    label = "Trimethoprim",
    CL = 0.48,
    V1 = 0.92,
    Q = 1.26,
    V2 = 0.86,
    KA_ORAL = 0.66,
    KA_IM = 1.14,
    F_ORAL = 0.58,
    F_IM = 0.99,
    fu = 0.489,
    beta_cl = 0.78,
    beta_v1 = 1.33
  ),
  SDZ = list(
    label = "Sulfadiazine",
    CL = 0.12,
    V1 = 0.30,
    Q = 0.32,
    V2 = 0.29,
    KA_ORAL = 0.50,
    KA_IM = 1.30,
    F_ORAL = 0.693,
    F_IM = 0.99,
    fu = 0.708,
    beta_cl = 0.51,
    beta_v1 = 1.10
  ),
  SMX = list(
    label = "Sulfamethoxazole",
    CL = 0.21,
    V1 = 0.48,
    Q = 0.83,
    V2 = 0.17,
    KA_ORAL = 1.82,
    F_ORAL = 0.64,
    fu = 0.427
  ),
  SDMX = list(
    label = "Sulfadimethoxine",
    CL = 0.015,
    V1 = 0.13,
    Q = 0.20,
    V2 = 0.16,
    KA_ORAL = 0.60,
    F_ORAL = 0.68,
    fu = 0.059
  )
)

paper_reference <- data.frame(
  combo_id = c(
    "TMP_SDZ", "TMP_SDZ",
    "TMP_SMX", "TMP_SMX",
    "TMP_SDMX", "TMP_SDMX",
    "TMP_ALL", "TMP_ALL"
  ),
  analyte = c("SDZ", "SDZ", "SMX", "SMX", "SDMX", "SDMX", "TMP", "TMP"),
  route = c("IV", "Oral", "IV", "Oral", "IV", "Oral", "IV", "Oral"),
  auc_ref = c(99.2, 144.3, 140.93, 92.39, 1280.3, 1629.0, 6.9, 7.3),
  thalf_ref = c(3.7, 4.4, 2.2, 2.3, 14.8, 13.6, 2.9, 3.5),
  stringsAsFactors = FALSE
)

ratio_reference <- data.frame(
  combo_id = c("TMP_SMX", "TMP_SDZ", "TMP_SDMX"),
  pct_in_band = c(8.8, 46.8, 76.5),
  duration_h = c(4.5, 14.0, 12.0),
  stringsAsFactors = FALSE
)

route_choices_for_combo <- function(combo_id) {
  row <- combo_catalog[combo_catalog$combo_id == combo_id, , drop = FALSE]
  choices <- c()
  if (!is.na(row$iv_tmp)) {
    choices["IV"] <- "IV"
  }
  if (!is.na(row$oral_tmp)) {
    choices["Oral (gavage)"] <- "Oral"
  }
  if (!is.na(row$im_tmp)) {
    choices["IM"] <- "IM"
  }
  choices
}

dose_defaults <- function(combo_id, route) {
  row <- combo_catalog[combo_catalog$combo_id == combo_id, , drop = FALSE]
  prefix <- switch(route, IV = "iv", Oral = "oral", IM = "im")
  c(
    TMP = row[[paste0(prefix, "_tmp")]],
    SULFONAMIDE = row[[paste0(prefix, "_s")]]
  )
}

combo_row <- function(combo_id) {
  combo_catalog[combo_catalog$combo_id == combo_id, , drop = FALSE]
}

reference_row <- function(combo_id, analyte_name, route_name) {
  key_combo <- if (analyte_name == "TMP") "TMP_ALL" else combo_id
  paper_reference |>
    dplyr::filter(
      .data$combo_id == key_combo,
      .data$analyte == .env$analyte_name,
      .data$route == .env$route_name
    )
}

terminal_lambda <- function(analyte, bw_kg) {
  pars <- analyte_parameters[[analyte]]
  cl_kg <- pars$CL
  v1_kg <- pars$V1

  if (!is.null(pars$beta_cl)) {
    cl_kg <- cl_kg * (bw_kg / 31.1)^pars$beta_cl
  }
  if (!is.null(pars$beta_v1)) {
    v1_kg <- v1_kg * (bw_kg / 31.1)^pars$beta_v1
  }

  cl <- cl_kg * bw_kg
  v1 <- v1_kg * bw_kg
  q <- pars$Q * bw_kg
  v2 <- pars$V2 * bw_kg

  k10 <- cl / v1
  k12 <- q / v1
  k21 <- q / v2
  disc <- (k10 + k12 + k21)^2 - 4 * k21 * k10
  0.5 * ((k10 + k12 + k21) - sqrt(disc))
}

parameterize_analyte <- function(analyte, route, bw_kg) {
  pars <- analyte_parameters[[analyte]]
  cl_kg <- pars$CL
  v1_kg <- pars$V1

  if (!is.null(pars$beta_cl)) {
    cl_kg <- cl_kg * (bw_kg / 31.1)^pars$beta_cl
  }
  if (!is.null(pars$beta_v1)) {
    v1_kg <- v1_kg * (bw_kg / 31.1)^pars$beta_v1
  }

  ka <- switch(
    route,
    IV = 1,
    Oral = pars$KA_ORAL,
    IM = pars$KA_IM
  )

  f_abs <- switch(
    route,
    IV = 1,
    Oral = pars$F_ORAL,
    IM = pars$F_IM
  )

  list(
    CL = cl_kg * bw_kg,
    V1 = v1_kg * bw_kg,
    Q = pars$Q * bw_kg,
    V2 = pars$V2 * bw_kg,
    KA = ka,
    F = f_abs,
    fu = pars$fu,
    label = pars$label
  )
}

simulate_one_analyte <- function(
  analyte,
  combo_id,
  route,
  bw_kg,
  dose_mg_kg,
  n_doses = 1,
  interval_h = 24,
  sim_end_h = 56,
  delta_h = 0.1
) {
  pars <- parameterize_analyte(analyte, route, bw_kg)
  total_dose_mg <- dose_mg_kg * bw_kg
  last_dose_time <- (n_doses - 1) * interval_h
  sim_end_h <- max(sim_end_h, last_dose_time + ifelse(n_doses > 1, interval_h, 0))

  event_obj <- if (route == "IV") {
    mrgsolve::ev(amt = total_dose_mg, cmt = 2, ii = interval_h, addl = n_doses - 1)
  } else {
    mrgsolve::ev(
      amt = total_dose_mg * pars$F,
      cmt = 1,
      ii = interval_h,
      addl = n_doses - 1
    )
  }

  mod |>
    mrgsolve::param(
      CL = pars$CL,
      V1 = pars$V1,
      Q = pars$Q,
      V2 = pars$V2,
      KA = pars$KA
    ) |>
    mrgsolve::ev(event_obj) |>
    mrgsolve::mrgsim(end = sim_end_h, delta = delta_h) |>
    as.data.frame() |>
    dplyr::mutate(
      combo_id = combo_id,
      route = route,
      analyte = analyte,
      analyte_label = pars$label,
      time_h = .data$time,
      total_ug_ml = .data$CP,
      free_ug_ml = .data$CP * pars$fu,
      bw_kg = bw_kg,
      dose_mg_kg = dose_mg_kg,
      last_dose_time = last_dose_time,
      interval_h = interval_h,
      n_doses = n_doses
    ) |>
    dplyr::distinct(.data$time_h, .data$analyte, .keep_all = TRUE)
}

simulate_profile <- function(
  combo_id = "TMP_SDZ",
  route = "Oral",
  bw_kg = 31.1,
  dose_tmp_mg_kg = NA_real_,
  dose_s_mg_kg = NA_real_,
  n_doses = 1,
  interval_h = 24,
  sim_end_h = NULL,
  delta_h = 0.1
) {
  row <- combo_row(combo_id)
  sulfonamide <- row$sulfonamide[[1]]

  defaults <- dose_defaults(combo_id, route)
  if (is.na(dose_tmp_mg_kg)) {
    dose_tmp_mg_kg <- defaults[["TMP"]]
  }
  if (is.na(dose_s_mg_kg)) {
    dose_s_mg_kg <- defaults[["SULFONAMIDE"]]
  }
  if (is.null(sim_end_h)) {
    sim_end_h <- row$default_end_h[[1]]
  }

  tmp_df <- simulate_one_analyte(
    analyte = "TMP",
    combo_id = combo_id,
    route = route,
    bw_kg = bw_kg,
    dose_mg_kg = dose_tmp_mg_kg,
    n_doses = n_doses,
    interval_h = interval_h,
    sim_end_h = sim_end_h,
    delta_h = delta_h
  )

  s_df <- simulate_one_analyte(
    analyte = sulfonamide,
    combo_id = combo_id,
    route = route,
    bw_kg = bw_kg,
    dose_mg_kg = dose_s_mg_kg,
    n_doses = n_doses,
    interval_h = interval_h,
    sim_end_h = sim_end_h,
    delta_h = delta_h
  )

  combined_long <- dplyr::bind_rows(tmp_df, s_df)

  ratio_df <- tmp_df |>
    dplyr::transmute(time_h = .data$time_h, tmp_total = .data$total_ug_ml, tmp_free = .data$free_ug_ml) |>
    dplyr::inner_join(
      s_df |>
        dplyr::transmute(time_h = .data$time_h, s_total = .data$total_ug_ml, s_free = .data$free_ug_ml),
      by = "time_h"
    ) |>
    dplyr::mutate(
      total_ratio = .data$tmp_total / .data$s_total,
      free_ratio = .data$tmp_free / .data$s_free
    )

  list(
    combined = combined_long,
    ratio = ratio_df,
    sulfonamide = sulfonamide,
    defaults = defaults
  )
}

build_paper_grid <- function() {
  dplyr::bind_rows(
    data.frame(combo_id = "TMP_SDZ", route = c("IV", "Oral", "IM"), dose_tmp_mg_kg = c(2.5, 5, 2.5), dose_s_mg_kg = c(12.5, 25, 12.5)),
    data.frame(combo_id = "TMP_SMX", route = c("IV", "Oral"), dose_tmp_mg_kg = c(6, 6), dose_s_mg_kg = c(30, 30)),
    data.frame(combo_id = "TMP_SDMX", route = c("IV", "Oral"), dose_tmp_mg_kg = c(4, 8), dose_s_mg_kg = c(18.6, 37.36))
  ) |>
    dplyr::mutate(combo_label = combo_catalog$label[match(.data$combo_id, combo_catalog$combo_id)])
}

exposure_metrics <- function(profile_df, analyte_name, bw_kg) {
  d <- profile_df |>
    dplyr::filter(.data$analyte == .env$analyte_name)

  if (d$n_doses[[1]] > 1) {
    window <- d |>
      dplyr::filter(
        .data$time_h >= .data$last_dose_time[[1]],
        .data$time_h <= .data$last_dose_time[[1]] + .data$interval_h[[1]]
      )
    auc_value <- sum(
      diff(window$time_h) *
        (head(window$total_ug_ml, -1) + tail(window$total_ug_ml, -1)) / 2
    )
    ctrough <- tail(window$total_ug_ml, 1)
  } else {
    window <- d
    auc_value <- sum(
      diff(window$time_h) *
        (head(window$total_ug_ml, -1) + tail(window$total_ug_ml, -1)) / 2
    )
    beta <- terminal_lambda(analyte_name, bw_kg)
    auc_value <- auc_value + tail(window$total_ug_ml, 1) / beta
    ctrough <- tail(window$total_ug_ml, 1)
  }

  beta <- terminal_lambda(analyte_name, bw_kg)
  list(
    cmax = max(window$total_ug_ml, na.rm = TRUE),
    ctrough = ctrough,
    auc = auc_value,
    tmax_h = window$time_h[[which.max(window$total_ug_ml)]] -
      ifelse(d$n_doses[[1]] > 1, d$last_dose_time[[1]], 0),
    thalf_h = log(2) / beta
  )
}

ratio_band_stats <- function(ratio_df) {
  band_df <- ratio_df |>
    dplyr::filter(.data$time_h <= 24, is.finite(.data$free_ratio), .data$free_ratio > 0)

  lower_bound <- 1 / 50
  upper_bound <- 1 / 10
  in_band <- band_df$free_ratio >= lower_bound & band_df$free_ratio <= upper_bound
  duration_h <- if (nrow(band_df) > 1) {
    sum(diff(band_df$time_h) * head(in_band, -1))
  } else {
    0
  }

  list(
    peak_free_ratio = max(ratio_df$free_ratio[is.finite(ratio_df$free_ratio)], na.rm = TRUE),
    min_free_ratio = min(ratio_df$free_ratio[is.finite(ratio_df$free_ratio) & ratio_df$free_ratio > 0], na.rm = TRUE),
    duration_h = duration_h,
    pct_24h = 100 * duration_h / 24
  )
}

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;
    }
    .ref-box a { color: #8b5cf6; }
    .small-note {
      font-size: 12px;
      color: #6c757d;
      margin-top: 10px;
    }
  ")

ui <- page_sidebar(
  title = "TMP/Sulfonamide Pig PopPK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    width = 340,
    accordion(
      accordion_panel(
        "Study Scenario",
        selectInput(
          "combo_id",
          "Combination",
          choices = stats::setNames(combo_catalog$combo_id, combo_catalog$label),
          selected = "TMP_SDZ"
        ),
        selectInput("route", "Route", choices = route_choices_for_combo("TMP_SDZ"), selected = "Oral"),
        selectInput("focus_analyte", "Metric focus", choices = c("TMP", "SDZ"), selected = "SDZ"),
        sliderInput("bw_kg", "Bodyweight (kg)", min = 20, max = 50, value = 31.1, step = 0.1)
      ),
      accordion_panel(
        "Dosing",
        numericInput("dose_tmp_mg_kg", "TMP dose (mg/kg)", value = 5, min = 0.1, step = 0.1),
        numericInput("dose_s_mg_kg", "Sulfonamide dose (mg/kg)", value = 25, min = 0.1, step = 0.1),
        numericInput("n_doses", "Number of doses", value = 1, min = 1, max = 30, step = 1),
        numericInput("interval_h", "Dosing interval (h)", value = 24, min = 4, max = 48, step = 1),
        helpText("Default doses match the paper's licensed regimens for the selected combination and route.")
      ),
      accordion_panel(
        "Display",
        numericInput("sim_end_h", "Simulation horizon (h)", value = 56, min = 12, max = 168, step = 4),
        radioButtons(
          "conc_basis",
          "Concentration basis",
          choices = c("Total concentrations" = "total", "Free concentrations" = "free"),
          selected = "total"
        ),
        checkboxInput("log_scale", "Log scale for concentration plot", value = FALSE)
      ),
      open = c("Study Scenario", "Dosing", "Display")
    ),
    div(
      class = "small-note",
      "Bodyweight effects from the paper are applied only to TMP and SDZ central volume and clearance."
    )
  ),
  navset_card_underline(
    nav_panel(
      "Simulation",
      layout_column_wrap(
        width = 1 / 4,
        fill = FALSE,
        metric_card("cmax_text", "Cmax (ug/mL)", "metric-success"),
        metric_card("ctrough_text", "Ctrough (ug/mL)", "metric-warning"),
        metric_card("auc_text", "AUC exposure (ug*h/mL)", "metric-primary"),
        metric_card("thalf_text", "t1/2 (h)", "metric-info")
      ),
      card(
        full_screen = TRUE,
        card_header("Analyte concentration-time profiles"),
        plotOutput("pk_plot", height = "500px")
      ),
      card(
        card_header("Paper anchor check"),
        tableOutput("anchor_table")
      )
    ),
    nav_panel(
      "Ratio Summary",
      layout_column_wrap(
        width = 1 / 4,
        fill = FALSE,
        metric_card("peak_ratio_text", "Peak free TMP:S ratio", "metric-primary"),
        metric_card("min_ratio_text", "Lowest free TMP:S ratio", "metric-info"),
        metric_card("band_duration_text", "Time in 1:10-1:50 band (24 h)", "metric-success"),
        metric_card("band_pct_text", "Percent of first 24 h in band", "metric-warning")
      ),
      card(
        full_screen = TRUE,
        card_header("TMP:sulfonamide ratio"),
        plotOutput("ratio_plot", height = "500px")
      ),
      card(
        card_header("Paper oral population benchmark"),
        tableOutput("ratio_anchor_table")
      )
    ),
    nav_panel(
      "Model Information",
      div(
        class = "ref-box",
        tags$h5("Population PK model reproduced from Boulanger et al. 2025"),
        tags$p(
          "Three sulfonamide combinations were studied in pigs: TMP/SDZ, TMP/SMX, and TMP/SDMX. ",
          "A separate 2-compartment model was fitted for each analyte, with oral bioavailability for all drugs and IM absorption for TMP/SDZ only."
        ),
        tags$ul(
          tags$li("Structural model: 2-compartment disposition for TMP, SDZ, SMX, and SDMX"),
          tags$li("Absorption: first-order oral absorption for all combinations; first-order IM absorption for TMP/SDZ"),
          tags$li("Bodyweight covariate: positive effect on TMP and SDZ clearance and central volume"),
          tags$li("Protein binding: TMP 51.2%, SDZ 29.2%, SMX 57.3%, SDMX 94.1%"),
          tags$li("Implementation note: oral SDZ bioavailability was calibrated to 0.693 so the simulator reproduces the paper's Table 3 oral AUC∞; the Table 2 fixed effect of 0.93 is not internally consistent with the reported oral SDZ exposure.")
        ),
        tags$h6("Implemented simulator parameters"),
        tags$p(
          class = "small-note",
          "Most values come directly from Table 2. The one explicit override is SDZ oral bioavailability, calibrated from 0.93 to 0.693 so the simulator matches the paper's reported oral SDZ AUC∞ in Table 3."
        ),
        tableOutput("parameter_table")
      )
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Primary source"),
        tags$p(
          "Boulanger M, Taillandier JF, Henri J, et al. Population pharmacokinetic modeling of sulfadimethoxine, sulfadiazine and sulfamethoxazole combined to trimethoprim in pigs. ",
          tags$em("Veterinary Quarterly"), " 2025;45:2565351. ",
          tags$a(href = "https://doi.org/10.1080/01652176.2025.2565351", target = "_blank", "doi:10.1080/01652176.2025.2565351")
        ),
        tags$h5("Key extracted findings"),
        tags$ul(
          tags$li("Only TMP/SDZ stayed relatively close to the historical 1:19 target ratio over time."),
          tags$li("Free TMP:S ratios stayed within 1:10-1:50 in 8.8% of pigs for TMP/SMX, 46.8% for TMP/SDZ, and 76.5% for TMP/SDMX."),
          tags$li("SDMX had the slowest clearance and longest half-life, driving progressive ratio decline after repeated dosing.")
        ),
        tags$h5("Data availability"),
        tags$p(
          tags$a(href = "https://doi.org/10.57745/Q2H4LQ", target = "_blank", "Study PK dataset DOI: 10.57745/Q2H4LQ")
        )
      )
    )
  ),
  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 for Husain Attarwala",
    br(),
    tags$span(
      style = "font-size: 10px;",
      "For research and educational purposes only. Not for veterinary prescribing."
    )
  )
)

server <- function(input, output, session) {
  observeEvent(input$combo_id, {
    row <- combo_row(input$combo_id)
    routes <- route_choices_for_combo(input$combo_id)
    updateSelectInput(session, "route", choices = routes, selected = unname(routes)[1])
    updateSelectInput(
      session,
      "focus_analyte",
      choices = c("TMP", row$sulfonamide[[1]]),
      selected = row$sulfonamide[[1]]
    )
  }, ignoreInit = FALSE)

  observeEvent(list(input$combo_id, input$route), {
    shiny::req(input$combo_id, input$route)
    row <- combo_row(input$combo_id)
    defaults <- dose_defaults(input$combo_id, input$route)
    updateNumericInput(session, "dose_tmp_mg_kg", value = defaults[["TMP"]])
    updateNumericInput(session, "dose_s_mg_kg", value = defaults[["SULFONAMIDE"]])
    updateNumericInput(session, "sim_end_h", value = row$default_end_h[[1]])
    updateSelectInput(
      session,
      "focus_analyte",
      choices = c("TMP", row$sulfonamide[[1]]),
      selected = row$sulfonamide[[1]]
    )
  }, ignoreInit = FALSE)

  sim_bundle <- reactive({
    shiny::req(
      input$combo_id,
      input$route,
      input$bw_kg,
      input$dose_tmp_mg_kg,
      input$dose_s_mg_kg,
      input$n_doses,
      input$interval_h,
      input$sim_end_h
    )

    simulate_profile(
      combo_id = input$combo_id,
      route = input$route,
      bw_kg = input$bw_kg,
      dose_tmp_mg_kg = input$dose_tmp_mg_kg,
      dose_s_mg_kg = input$dose_s_mg_kg,
      n_doses = input$n_doses,
      interval_h = input$interval_h,
      sim_end_h = input$sim_end_h
    )
  })

  focus_metrics <- reactive({
    shiny::req(input$focus_analyte)
    exposure_metrics(sim_bundle()$combined, input$focus_analyte, input$bw_kg)
  })

  ratio_metrics <- reactive({
    ratio_band_stats(sim_bundle()$ratio)
  })

  output$cmax_text <- renderText({
    sprintf("%.1f", focus_metrics()$cmax)
  })

  output$ctrough_text <- renderText({
    sprintf("%.2f", focus_metrics()$ctrough)
  })

  output$auc_text <- renderText({
    sprintf("%.1f", focus_metrics()$auc)
  })

  output$thalf_text <- renderText({
    sprintf("%.1f", focus_metrics()$thalf_h)
  })

  output$peak_ratio_text <- renderText({
    sprintf("%.3f", ratio_metrics()$peak_free_ratio)
  })

  output$min_ratio_text <- renderText({
    sprintf("%.3f", ratio_metrics()$min_free_ratio)
  })

  output$band_duration_text <- renderText({
    sprintf("%.1f h", ratio_metrics()$duration_h)
  })

  output$band_pct_text <- renderText({
    sprintf("%.1f%%", ratio_metrics()$pct_24h)
  })

  output$pk_plot <- renderPlot({
    bundle <- sim_bundle()
    y_col <- if (input$conc_basis == "free") "free_ug_ml" else "total_ug_ml"
    label <- if (input$conc_basis == "free") "Free concentration (ug/mL)" else "Total concentration (ug/mL)"

    plot_df <- bundle$combined |>
      dplyr::mutate(plot_value = .data[[y_col]])

    if (isTRUE(input$log_scale)) {
      plot_df <- plot_df |>
        dplyr::filter(.data$plot_value > 0)
    }

    ggplot(plot_df, aes(x = .data$time_h, y = .data$plot_value, color = .data$analyte_label)) +
      geom_line(linewidth = 1.1) +
      scale_color_manual(values = c("Trimethoprim" = "#8b5cf6", "Sulfadiazine" = "#10b981", "Sulfamethoxazole" = "#f59e0b", "Sulfadimethoxine" = "#0ea5e9")) +
      labs(
        x = "Time (h)",
        y = label,
        color = NULL,
        title = paste(combo_row(input$combo_id)$label[[1]], "in pigs"),
        subtitle = paste("Route:", input$route, "• Bodyweight:", sprintf("%.1f kg", input$bw_kg))
      ) +
      theme_minimal(base_size = 14) +
      if (isTRUE(input$log_scale)) scale_y_log10() else NULL +
      theme(legend.position = "top")
  })

  output$ratio_plot <- renderPlot({
    ratio_df <- sim_bundle()$ratio |>
      dplyr::filter(is.finite(.data$total_ratio), is.finite(.data$free_ratio), .data$total_ratio > 0, .data$free_ratio > 0) |>
      tidyr::pivot_longer(
        cols = c("total_ratio", "free_ratio"),
        names_to = "series",
        values_to = "ratio_value"
      ) |>
      dplyr::mutate(
        series = dplyr::recode(.data$series, total_ratio = "Total TMP:S ratio", free_ratio = "Free TMP:S ratio")
      )

    ggplot(ratio_df, aes(x = .data$time_h, y = .data$ratio_value, color = .data$series)) +
      annotate("rect", xmin = -Inf, xmax = 24, ymin = 1 / 50, ymax = 1 / 10, fill = "#10b981", alpha = 0.10) +
      geom_hline(yintercept = 1 / 19, linetype = "dashed", color = "#dc2626", linewidth = 0.8) +
      geom_line(linewidth = 1.1) +
      scale_color_manual(values = c("Total TMP:S ratio" = "#6b7280", "Free TMP:S ratio" = "#8b5cf6")) +
      scale_y_log10() +
      labs(
        x = "Time (h)",
        y = "TMP:sulfonamide ratio",
        color = NULL,
        title = "Ratio behavior over time",
        subtitle = "Green band = TMP:S 1:10 to 1:50 (0.10 to 0.02); red dashed line = TMP:S 1:19"
      ) +
      theme_minimal(base_size = 14) +
      theme(legend.position = "top")
  })

  output$anchor_table <- renderTable({
    ref <- reference_row(input$combo_id, input$focus_analyte, input$route)
    if (nrow(ref) == 0) {
      return(data.frame(Note = "No route-specific paper AUC/t1/2 anchor reported for this analyte."))
    }

    metrics <- focus_metrics()
    data.frame(
      Metric = c("Simulated AUC", "Paper AUC", "Simulated t1/2", "Paper t1/2"),
      Value = c(
        sprintf("%.1f ug*h/mL", metrics$auc),
        sprintf("%.1f ug*h/mL", ref$auc_ref[[1]]),
        sprintf("%.1f h", metrics$thalf_h),
        sprintf("%.1f h", ref$thalf_ref[[1]])
      ),
      stringsAsFactors = FALSE
    )
  }, striped = TRUE, bordered = TRUE, spacing = "s")

  output$ratio_anchor_table <- renderTable({
    ref <- ratio_reference |>
      dplyr::filter(.data$combo_id == input$combo_id)

    if (nrow(ref) == 0) {
      return(data.frame(Note = "No oral population ratio benchmark available for this combination."))
    }

    if (input$route != "Oral") {
      return(data.frame(
        Note = "These paper benchmarks were reported for oral dosing only and are not directly comparable to the current non-oral single-pig simulation."
      ))
    }

    data.frame(
      Paper_metric = c(
        "Population median time in TMP:S 1:10-1:50 band over first 24 h",
        "Population percent of pigs with TMP:S in 1:10-1:50 band"
      ),
      Value = c(
        sprintf("%.1f h", ref$duration_h[[1]]),
        sprintf("%.1f%%", ref$pct_in_band[[1]])
      ),
      stringsAsFactors = FALSE
    )
  }, striped = TRUE, bordered = TRUE, spacing = "s")

  output$parameter_table <- renderTable({
    data.frame(
      Analyte = c("TMP", "SDZ", "SMX", "SDMX"),
      CL_L_h_kg = c(0.48, 0.12, 0.21, 0.015),
      V1_L_kg = c(0.92, 0.30, 0.48, 0.13),
      Q_L_h_kg = c(1.26, 0.32, 0.83, 0.20),
      V2_L_kg = c(0.86, 0.29, 0.17, 0.16),
      KA_oral_h = c(0.66, 0.50, 1.82, 0.60),
      F_oral_reported = c(0.58, 0.93, 0.64, 0.68),
      F_oral_impl = c(0.58, 0.693, 0.64, 0.68),
      stringsAsFactors = FALSE
    )
  }, striped = TRUE, bordered = TRUE, spacing = "s")
}

shinyApp(ui, server)