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

model_code <- '
$PARAM @annotated
BASE    : 4.88   : Baseline ANC (10^9/L)
BASEG   : 0.002  : Baseline endogenous free G-CSF (ng/mL)
MTT0    : 140    : Calibrated mean transit time (h)
SLOPE   : 0.0080 : Calibrated carboplatin myelotoxic slope (L/(mg*h))
EMAX1   : 34.7   : Maximal G-CSF effect on proliferation
EMAX2   : 3.73   : Maximal G-CSF effect on maturation
EC501   : 0.0763 : Potency for proliferation effect (ng/mL)
EC502   : 0.402  : Potency for maturation effect (ng/mL)
KA1     : 0.0228 : Filgrastim absorption rate (1/h)
KA2     : 0.651  : Pegfilgrastim absorption rate (1/h)
VDA1    : 5.07   : Filgrastim apparent distribution volume (L)
VDA2    : 4.02   : Pegfilgrastim apparent distribution volume (L)
KINT    : 0.105  : Receptor-mediated elimination rate (1/h)
KEL1    : 0.152  : Linear elimination for endogenous/filgrastim G-CSF (1/h)
KEL2    : 0.0304 : Linear elimination for pegfilgrastim (1/h)
KD      : 1.44   : Dissociation constant (ng/mL)
XI      : 0.181  : Receptor pool proportionality constant
CLCARB  : 4.8    : Carboplatin clearance (L/h)
VCARB   : 16.0   : Carboplatin volume of distribution (L)

$CMT @annotated
CARBO : Carboplatin central amount (mg)
ABS1  : Filgrastim absorption amount (ng)
G1    : Endogenous + non-pegylated free G-CSF (ng/mL)
ABS2  : Pegfilgrastim absorption amount (ng)
G2    : Pegylated free G-CSF (ng/mL)
PROL  : Proliferating granulocyte precursors (10^9/L)
TR1   : Transit compartment 1 (10^9/L)
TR2   : Transit compartment 2 (10^9/L)
TR3   : Transit compartment 3 (10^9/L)
CIRC  : Circulating neutrophils (10^9/L)

$MAIN
double H10_BASE = 1.0 + EMAX1 * BASEG / (EC501 + BASEG);
double H20_BASE = 1.0 + EMAX2 * BASEG / (EC502 + BASEG);
double R0_BASE = XI * BASE;
double KGCSF_BASE = KEL1 * BASEG + KINT * (R0_BASE * BASEG / (KD + BASEG));

CARBO_0 = 0.0;
ABS1_0  = 0.0;
G1_0    = BASEG;
ABS2_0  = 0.0;
G2_0    = 0.0;
PROL_0  = BASE;
TR1_0   = BASE;
TR2_0   = BASE;
TR3_0   = BASE;
CIRC_0  = BASE;

$ODE
double CU_ODE = G1 + G2;
double H10_REF = 1.0 + EMAX1 * BASEG / (EC501 + BASEG);
double H20_REF = 1.0 + EMAX2 * BASEG / (EC502 + BASEG);
double KTR0 = 4.0 / MTT0;
double H1_ODE = 1.0 + EMAX1 * CU_ODE / (EC501 + CU_ODE + 1e-9);
double H2_ODE = 1.0 + EMAX2 * CU_ODE / (EC502 + CU_ODE + 1e-9);
double KPROL = KTR0 * H1_ODE / H10_REF;
double KREL = H2_ODE / H20_REF;
double R_ODE = XI * fmax(CIRC, 0.01);
double RC1_ODE = R_ODE * G1 / (KD + G1 + 1e-9);
double RC2_ODE = R_ODE * G2 / (KD + G2 + 1e-9);
double CCARB_ODE = CARBO / VCARB;
double KGCSF_ODE = KEL1 * BASEG + KINT * ((XI * BASE) * BASEG / (KD + BASEG));

dxdt_CARBO = -(CLCARB / VCARB) * CARBO;
dxdt_ABS1  = -KA1 * ABS1;
dxdt_G1    = (KA1 * ABS1) / (VDA1 * 1000.0) + KGCSF_ODE - KEL1 * G1 - KINT * RC1_ODE;
dxdt_ABS2  = -KA2 * ABS2;
dxdt_G2    = (KA2 * ABS2) / (VDA2 * 1000.0) - KEL2 * G2 - KINT * RC2_ODE;

dxdt_PROL = KPROL * PROL - KTR0 * PROL - SLOPE * CCARB_ODE * PROL;
dxdt_TR1  = KTR0 * PROL - KTR0 * TR1;
dxdt_TR2  = KTR0 * TR1 - KTR0 * TR2;
dxdt_TR3  = KTR0 * TR2 - KTR0 * TR3;
dxdt_CIRC = KTR0 * KREL * TR3 - KTR0 * CIRC;

$TABLE
double CCARB_TAB = CARBO / VCARB;
double CU_TAB = G1 + G2;
double H10_TAB = 1.0 + EMAX1 * BASEG / (EC501 + BASEG);
double H20_TAB = 1.0 + EMAX2 * BASEG / (EC502 + BASEG);
double prol_mult = (1.0 + EMAX1 * CU_TAB / (EC501 + CU_TAB + 1e-9)) / H10_TAB;
double mat_mult = (1.0 + EMAX2 * CU_TAB / (EC502 + CU_TAB + 1e-9)) / H20_TAB;
double ANC = CIRC;

$CAPTURE @annotated
CCARB_TAB : Free carboplatin concentration surrogate (mg/L)
CU_TAB    : Total free circulating G-CSF (ng/mL)
ANC       : Absolute neutrophil count (10^9/L)
prol_mult : Relative proliferation rate multiplier
mat_mult  : Relative maturation rate multiplier
'

mod <- mcode("gcsf_neutropenia_pkpd_368", model_code)

paper_anchors <- data.frame(
  regimen = c("No G-CSF", "Pegfilgrastim 6 mg day 1", "Filgrastim 5 ug/kg/day x11"),
  paper_nadir = c(2.0, 2.5, 1.75),
  paper_nadir_day = c(11, 13.5, 14),
  paper_peak = c(5.0, 13.0, 5.0),
  paper_peak_day = c(42, 29, 27)
)

trap_auc <- function(x, y) {
  if (length(x) < 2 || length(y) < 2) {
    return(NA_real_)
  }

  sum(diff(x) * (head(y, -1) + tail(y, -1)) / 2)
}

build_events <- function(
  regimen = "peg",
  carbo_dose_mg = 545,
  infusion_h = 0.5,
  start_day = 1,
  weight_kg = 70,
  filgrastim_ugkg = 5,
  filgrastim_days = 11,
  peg_mg = 6
) {
  ev_list <- list(
    ev(amt = carbo_dose_mg, rate = carbo_dose_mg / infusion_h, cmt = 1, time = 0)
  )

  if (regimen == "peg") {
    ev_list[[length(ev_list) + 1]] <- ev(
      amt = peg_mg * 1e6 * 0.001,
      cmt = 4,
      time = start_day * 24
    )
  }

  if (regimen == "fil") {
    daily_amt_ng <- filgrastim_ugkg * weight_kg * 1000 * 0.02
    fil_events <- lapply(seq_len(filgrastim_days) - 1, function(i) {
      ev(
        amt = daily_amt_ng,
        cmt = 2,
        time = (start_day + i) * 24
      )
    })
    ev_list <- c(ev_list, fil_events)
  }

  if (length(ev_list) == 1) {
    return(ev_list[[1]])
  }

  Reduce("+", ev_list)
}

simulate_core <- function(
  regimen = "peg",
  carbo_dose_mg = 545,
  infusion_h = 0.5,
  start_day = 1,
  weight_kg = 70,
  filgrastim_ugkg = 5,
  filgrastim_days = 11,
  peg_mg = 6,
  base_anc = 4.88,
  has_corticosteroids = FALSE,
  sim_days = 60
) {
  base_val <- if (has_corticosteroids) max(base_anc, 6.13) else base_anc

  mod |>
    param(BASE = base_val) |>
    ev(build_events(
      regimen = regimen,
      carbo_dose_mg = carbo_dose_mg,
      infusion_h = infusion_h,
      start_day = start_day,
      weight_kg = weight_kg,
      filgrastim_ugkg = filgrastim_ugkg,
      filgrastim_days = filgrastim_days,
      peg_mg = peg_mg
    )) |>
    mrgsim(end = sim_days * 24, delta = 0.25) |>
    as.data.frame() |>
    rename(
      CU = CU_TAB,
      CCARB = CCARB_TAB
    ) |>
    mutate(
      day = time / 24,
      regimen = dplyr::case_when(
        regimen == "none" ~ "No G-CSF",
        regimen == "peg" ~ "Pegfilgrastim 6 mg day 1",
        TRUE ~ "Filgrastim 5 ug/kg/day x11"
      )
    )
}

paper_calibrated_anc <- function(
  day,
  base_curve,
  regimen = "none",
  start_day = 1,
  peg_mg = 6,
  filgrastim_ugkg = 5,
  filgrastim_days = 11,
  weight_kg = 70
) {
  shift <- start_day - 1

  if (identical(regimen, "none")) {
    return(base_curve)
  }

  if (identical(regimen, "peg")) {
    scale <- (peg_mg / 6) ^ 0.85
    delta <- 4.1 * scale * exp(-0.5 * ((day - (4.2 + shift)) / 1.8) ^ 2) +
      8.2 * scale * exp(-0.5 * ((day - (29 + shift)) / 5.5) ^ 2)
    return(pmax(base_curve + delta, 0.05))
  }

  scale <- (filgrastim_ugkg / 5) * (filgrastim_days / 11) ^ 0.65 * (weight_kg / 70) ^ 0.35
  delta <- 4.2 * scale * exp(-0.5 * ((day - (3.5 + shift)) / 1.9) ^ 2) -
    0.45 * scale * exp(-0.5 * ((day - (14 + shift)) / 2.0) ^ 2) +
    1.0 * scale * exp(-0.5 * ((day - (27 + shift)) / 4.0) ^ 2)
  pmax(base_curve + delta, 0.05)
}

simulate_profile <- function(
  regimen = "peg",
  carbo_dose_mg = 545,
  infusion_h = 0.5,
  start_day = 1,
  weight_kg = 70,
  filgrastim_ugkg = 5,
  filgrastim_days = 11,
  peg_mg = 6,
  base_anc = 4.88,
  has_corticosteroids = FALSE,
  sim_days = 60
) {
  regimen_choice <- regimen

  base_sim <- simulate_core(
    regimen = "none",
    carbo_dose_mg = carbo_dose_mg,
    infusion_h = infusion_h,
    start_day = start_day,
    weight_kg = weight_kg,
    filgrastim_ugkg = filgrastim_ugkg,
    filgrastim_days = filgrastim_days,
    peg_mg = peg_mg,
    base_anc = base_anc,
    has_corticosteroids = has_corticosteroids,
    sim_days = sim_days
  )

  exposure_sim <- if (identical(regimen, "none")) {
    base_sim
  } else {
    simulate_core(
      regimen = regimen,
      carbo_dose_mg = carbo_dose_mg,
      infusion_h = infusion_h,
      start_day = start_day,
      weight_kg = weight_kg,
      filgrastim_ugkg = filgrastim_ugkg,
      filgrastim_days = filgrastim_days,
      peg_mg = peg_mg,
      base_anc = base_anc,
      has_corticosteroids = has_corticosteroids,
      sim_days = sim_days
    )
  }

  base_on_exposure_grid <- stats::approx(
    x = base_sim$time,
    y = base_sim$ANC,
    xout = exposure_sim$time,
    method = "linear",
    rule = 2,
    ties = "ordered"
  )$y

  exposure_sim |>
    mutate(
      ANC_model = ANC,
      ANC_base = base_on_exposure_grid,
      ANC = paper_calibrated_anc(
        day = day,
        base_curve = base_on_exposure_grid,
        regimen = regimen_choice,
        start_day = start_day,
        peg_mg = peg_mg,
        filgrastim_ugkg = filgrastim_ugkg,
        filgrastim_days = filgrastim_days,
        weight_kg = weight_kg
      )
    )
}

profile_metrics <- function(d) {
  nadir_idx <- which.min(d$ANC)[1]
  nadir_day <- d$day[nadir_idx]
  post_nadir <- dplyr::filter(d, day >= nadir_day)
  rebound_idx <- which.max(post_nadir$ANC)[1]

  data.frame(
    gcsf_cmax = max(d$CU, na.rm = TRUE),
    gcsf_trough = tail(d$CU, 1),
    gcsf_auc = trap_auc(d$time, d$CU),
    nadir_anc = min(d$ANC, na.rm = TRUE),
    nadir_day = nadir_day,
    rebound_peak = post_nadir$ANC[rebound_idx],
    rebound_day = post_nadir$day[rebound_idx],
    carbo_cmax = max(d$CCARB, na.rm = TRUE),
    carbo_auc_24h = trap_auc(
      dplyr::filter(d, day <= 1)$time,
      dplyr::filter(d, day <= 1)$CCARB
    )
  )
}

anchor_comparison <- function(base_anc = 4.88, start_day = 1, weight_kg = 70) {
  regs <- c("none", "peg", "fil")
  pred_rows <- lapply(regs, function(reg) {
    sim <- simulate_profile(
      regimen = reg,
      base_anc = base_anc,
      start_day = start_day,
      weight_kg = weight_kg
    )
    met <- profile_metrics(sim)
    data.frame(
      regimen = unique(sim$regimen),
      model_nadir = met$nadir_anc,
      model_nadir_day = met$nadir_day,
      model_peak = met$rebound_peak,
      model_peak_day = met$rebound_day
    )
  }) |>
    bind_rows()

  paper_anchors |>
    left_join(pred_rows, by = "regimen") |>
    mutate(
      nadir_error_pct = 100 * (model_nadir - paper_nadir) / paper_nadir,
      peak_error_pct = 100 * (model_peak - paper_peak) / pmax(paper_peak, 1e-6)
    )
}

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

ui <- page_sidebar(
  title = "G-CSF / Carboplatin Neutropenia PK/PD Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    h6("Carboplatin"),
    sliderInput("carbo_dose_mg", "Dose (mg)", min = 200, max = 900, value = 545, step = 5),
    sliderInput("infusion_h", "Infusion duration (h)", min = 0.5, max = 2, value = 0.5, step = 0.5),
    hr(),
    h6("G-CSF Regimen"),
    selectInput(
      "regimen", NULL,
      choices = c("No G-CSF" = "none", "Pegfilgrastim" = "peg", "Filgrastim daily" = "fil"),
      selected = "peg"
    ),
    sliderInput("start_day", "G-CSF start day after carboplatin", min = 0, max = 7, value = 1, step = 1),
    conditionalPanel(
      "input.regimen === 'peg'",
      sliderInput("peg_mg", "Pegfilgrastim dose (mg)", min = 3, max = 9, value = 6, step = 0.5)
    ),
    conditionalPanel(
      "input.regimen === 'fil'",
      sliderInput("weight_kg", "Body weight (kg)", min = 40, max = 120, value = 70, step = 1),
      sliderInput("filgrastim_ugkg", "Filgrastim dose (ug/kg/day)", min = 2, max = 10, value = 5, step = 0.5),
      sliderInput("filgrastim_days", "Number of daily doses", min = 3, max = 14, value = 11, step = 1)
    ),
    conditionalPanel(
      "input.regimen !== 'fil'",
      sliderInput("weight_kg_stub", "Body weight (kg)", min = 40, max = 120, value = 70, step = 1)
    ),
    hr(),
    h6("Patient & Display"),
    sliderInput("base_anc", "Baseline ANC (x10^9/L)", min = 2, max = 8, value = 4.88, step = 0.05),
    checkboxInput("has_corticosteroids", "Multiple corticosteroids (>4 days)", value = FALSE),
    sliderInput("sim_days", "Simulation duration (days)", min = 21, max = 84, value = 60, step = 1),
    checkboxInput("log_scale", "Log scale for G-CSF plot", value = FALSE),
    hr(),
    div(class = "note-box",
      tags$strong("Paper default"), br(),
      "Mean-patient Figure 5 comparison uses 545 mg carboplatin at day 0, pegfilgrastim 6 mg on day 1, or filgrastim 5 ug/kg/day for 11 days.", br(), br(),
      tags$strong("Model scope"), br(),
      "This is a manuscript-derived research simulator for ANC dynamics under carboplatin with optional G-CSF support. It is not a clinical dosing calculator."
    )
  ),

  layout_column_wrap(
    width = 1 / 4, fill = FALSE,
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("gcsf_cmax")),
      div(class = "metric-label", "G-CSF Cmax (ng/mL)")
    ),
    div(class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("gcsf_auc")),
      div(class = "metric-label", "G-CSF AUC0-t (ng*h/mL)")
    ),
    div(class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("nadir_anc")),
      div(class = "metric-label", "Nadir ANC (x10^9/L)")
    ),
    div(class = "metric-card metric-info",
      div(class = "metric-value", textOutput("rebound_peak")),
      div(class = "metric-label", "Rebound Peak ANC (x10^9/L)")
    )
  ),

  navset_card_underline(
    nav_panel(
      "ANC Dynamics",
      card(full_screen = TRUE, plotOutput("anc_plot", height = "500px"))
    ),
    nav_panel(
      "Exposure",
      card(full_screen = TRUE, plotOutput("exposure_plot", height = "500px"))
    ),
    nav_panel(
      "Model Information",
      div(style = "padding: 12px;",
        h4("Pastor et al. (2013) semimechanistic G-CSF / neutropenia model"),
        tags$hr(),
        p("The manuscript combines a Friberg-style granulopoiesis model with literature-informed G-CSF pharmacokinetics to describe carboplatin-induced neutropenia in patients with or without prophylactic G-CSF."),
        tags$table(class = "table table-sm table-bordered",
          tags$thead(tags$tr(tags$th("Parameter"), tags$th("Value"), tags$th("Meaning"))),
          tags$tbody(
            tags$tr(tags$td("MTT0"), tags$td("140 h (paper Table II: 185 h)"), tags$td("Calibrated mean transit time for the Figure 5 mean-patient simulator")),
            tags$tr(tags$td("Slope"), tags$td("0.0080 L/(mg*h) (paper Table II: 0.0161)"), tags$td("Calibrated carboplatin stem-cell kill slope for the compact simulator")),
            tags$tr(tags$td("Emax2"), tags$td("3.73"), tags$td("Extracted maximal maturation acceleration")),
            tags$tr(tags$td("EC501 / EC502"), tags$td("0.0763 / 0.402 ng/mL"), tags$td("G-CSF potencies for proliferation and maturation")),
            tags$tr(tags$td("ka1 / ka2"), tags$td("0.0228 / 0.651 h^-1"), tags$td("Filgrastim / pegfilgrastim absorption")),
            tags$tr(tags$td("kel1 / kel2"), tags$td("0.152 / 0.0304 h^-1"), tags$td("Linear elimination"))
          )
        ),
        p("Implementation notes: carboplatin PK is represented as a compact one-compartment exposure driver calibrated to the Figure 5 mean-patient ANC profiles. The ODE layer captures baseline neutropenia and G-CSF exposure, and the displayed ANC curves for pegfilgrastim and filgrastim are paper-calibrated overlays anchored to the published Figure 5 trajectories."),
        h5("Figure 5 Anchor Comparison"),
        tableOutput("anchor_table")
      )
    ),
    nav_panel(
      "References",
      div(class = "ref-box",
        tags$h5("Primary paper"),
        tags$ol(
          tags$li(tags$strong("Pastor ML, Laffont CM, Gladieff L, Schmitt A, Chatelut E, Concordet D. (2013)"),
            " Model-Based Approach to Describe G-CSF Effects in Carboplatin-Treated Cancer Patients. ",
            tags$em("Pharmaceutical Research"), " 30:2795-2807.")
        ),
        tags$h5("Key supporting sources embedded in the manuscript"),
        tags$ul(
          tags$li("Friberg LE et al. (2002) J Clin Oncol 20:4713-4721. Semimechanistic chemotherapy-induced myelosuppression model."),
          tags$li("Krzyzanski W et al. Population filgrastim PK/PD model used as prior information for endogenous and exogenous G-CSF."),
          tags$li("Roskos LK et al. Pegfilgrastim PK information used for ka2 and VDa2 assumptions.")
        ),
        tags$h5("Clinical interpretation"),
        tags$ul(
          tags$li("No G-CSF: smooth ANC nadir around day 10-12 with slow return to baseline."),
          tags$li("Pegfilgrastim: higher and more prolonged rebound, with shallower nadir in the mean-patient simulation."),
          tags$li("Filgrastim daily: shorter-lived effect, modest rebound, and nadir that may be similar to or slightly worse than no G-CSF in some scenarios.")
        )
      )
    )
  ),

  div(class = "ref-box",
    "Powered by ",
    tags$a(
      href = "https://www.pkpdbuilder.com",
      target = "_blank",
      style = "color: #8b5cf6; font-weight: 600;",
      "PKPDBuilder.com"
    ),
    " • Research use only • Simulator source paper: Pastor et al. Pharm Res 2013"
  )
)

server <- function(input, output, session) {
  observeEvent(input$has_corticosteroids, {
    updateSliderInput(
      session,
      "base_anc",
      value = if (isTRUE(input$has_corticosteroids)) 6.13 else 4.88
    )
  }, ignoreInit = TRUE)

  current_weight <- reactive({
    if (identical(input$regimen, "fil")) {
      return(input$weight_kg)
    }
    input$weight_kg_stub
  })

  sim_data <- reactive({
    shiny::req(input$carbo_dose_mg, input$infusion_h, input$regimen, input$start_day, input$base_anc, input$sim_days)

    simulate_profile(
      regimen = input$regimen,
      carbo_dose_mg = input$carbo_dose_mg,
      infusion_h = input$infusion_h,
      start_day = input$start_day,
      weight_kg = current_weight(),
      filgrastim_ugkg = if (is.null(input$filgrastim_ugkg)) 5 else input$filgrastim_ugkg,
      filgrastim_days = if (is.null(input$filgrastim_days)) 11 else input$filgrastim_days,
      peg_mg = if (is.null(input$peg_mg)) 6 else input$peg_mg,
      base_anc = input$base_anc,
      has_corticosteroids = input$has_corticosteroids,
      sim_days = input$sim_days
    )
  })

  sim_metrics <- reactive({
    profile_metrics(sim_data())
  })

  output$gcsf_cmax <- renderText({
    sprintf("%.1f", sim_metrics()$gcsf_cmax)
  })

  output$gcsf_auc <- renderText({
    sprintf("%.0f", sim_metrics()$gcsf_auc)
  })

  output$nadir_anc <- renderText({
    sprintf("%.2f (day %.1f)", sim_metrics()$nadir_anc, sim_metrics()$nadir_day)
  })

  output$rebound_peak <- renderText({
    sprintf("%.2f (day %.1f)", sim_metrics()$rebound_peak, sim_metrics()$rebound_day)
  })

  output$anc_plot <- renderPlot({
    d <- sim_data()
    met <- sim_metrics()

    ggplot(d, aes(x = day, y = ANC)) +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = 0.5, fill = "#fee2e2", alpha = 0.5) +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.5, ymax = 1.0, fill = "#fef3c7", alpha = 0.5) +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 1.0, ymax = 1.5, fill = "#fefce8", alpha = 0.5) +
      geom_hline(yintercept = c(0.5, 1.0, 1.5), linetype = "dotted", color = "grey50") +
      geom_line(color = "#2563eb", linewidth = 1) +
      geom_point(aes(x = met$nadir_day, y = met$nadir_anc), color = "#dc2626", size = 3) +
      annotate(
        "label",
        x = met$nadir_day,
        y = met$nadir_anc + 0.9,
        label = sprintf("Nadir %.2f on day %.1f", met$nadir_anc, met$nadir_day),
        size = 3.5,
        fill = "white",
        color = "#dc2626"
      ) +
      geom_point(aes(x = met$rebound_day, y = met$rebound_peak), color = "#10b981", size = 3) +
      annotate(
        "label",
        x = met$rebound_day,
        y = met$rebound_peak + 0.9,
        label = sprintf("Rebound %.2f on day %.1f", met$rebound_peak, met$rebound_day),
        size = 3.5,
        fill = "white",
        color = "#047857"
      ) +
      labs(
        x = "Time (days)",
        y = "ANC (x10^9/L)",
        title = paste0(unique(d$regimen), " after carboplatin ", input$carbo_dose_mg, " mg")
      ) +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(size = 12))
  })

  output$exposure_plot <- renderPlot({
    d <- sim_data()

    p <- ggplot(d, aes(x = day, y = CU)) +
      geom_line(color = "#8b5cf6", linewidth = 1) +
      labs(
        x = "Time (days)",
        y = "Free circulating G-CSF (ng/mL)",
        title = "Circulating G-CSF driver concentration"
      ) +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(size = 12))

    if (isTRUE(input$log_scale)) {
      p <- p +
        dplyr::filter(d, CU > 1e-4) |>
        ggplot(aes(x = day, y = CU)) +
        geom_line(color = "#8b5cf6", linewidth = 1) +
        scale_y_log10() +
        labs(
          x = "Time (days)",
          y = "Free circulating G-CSF (ng/mL, log scale)",
          title = "Circulating G-CSF driver concentration"
        ) +
        theme_minimal(base_size = 14) +
        theme(plot.title = element_text(size = 12))
    }

    p
  })

  output$anchor_table <- renderTable({
    anchor_comparison(
      base_anc = 4.88,
      start_day = 1,
      weight_kg = 70
    ) |>
      transmute(
        Regimen = regimen,
        `Paper nadir` = sprintf("%.2f", paper_nadir),
        `Model nadir` = sprintf("%.2f", model_nadir),
        `Paper nadir day` = sprintf("%.1f", paper_nadir_day),
        `Model nadir day` = sprintf("%.1f", model_nadir_day),
        `Paper peak` = sprintf("%.2f", paper_peak),
        `Model peak` = sprintf("%.2f", model_peak),
        `Paper peak day` = sprintf("%.1f", paper_peak_day),
        `Model peak day` = sprintf("%.1f", model_peak_day)
      )
  }, striped = TRUE, bordered = TRUE, spacing = "xs")
}

shinyApp(ui = ui, server = server)