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

model_code <- '
$PARAM @annotated
V1      : 2.9     : Plasma volume (L)
VI      : 12.0    : Interstitial volume (L)
CL_lin  : 0.00938 : Linear clearance (L/h)
Q       : 0.05    : Intercompartmental CL (L/h)
kon_m   : 2520    : On-rate membrane HER2 (1/umol/L/h)
koff_m  : 1.26    : Off-rate membrane HER2 (1/h)
kint_m  : 0.10    : Internalization rate membrane complex (1/h)
kdeg_m  : 0.10    : Degradation rate membrane HER2 (1/h)
ksyn_m  : 0.000357 : Synthesis rate membrane HER2 (umol/L/h)
ksh_p   : 0.00202 : Shedding rate constant (1/h)
lambda  : 0.01    : Lymphatic transport ECD IS to plasma (1/h)
kdeg_p  : 0.20    : Degradation plasma ECD HER2 (1/h)
kon_s   : 2520    : On-rate soluble IS HER2 (1/umol/L/h)
koff_s  : 1.26    : Off-rate soluble IS HER2 (1/h)
kelm_s  : 0.01    : Elimination drug-soluble IS complex (1/h)
kon_p   : 2520    : On-rate soluble plasma HER2 (1/umol/L/h)
koff_p  : 1.26    : Off-rate soluble plasma HER2 (1/h)
kelm_p  : 0.20    : Elimination drug-soluble plasma complex (1/h)

$CMT @annotated
AP    : Drug amount plasma (umol)
AI    : Drug amount interstitial (umol)
RM    : Membrane-bound HER2 amount in IS (umol)
DRM   : Drug-membrane HER2 complex amount (umol)
RS    : Shed HER2 amount in IS (umol)
DRS   : Drug-soluble IS HER2 complex amount (umol)
RP    : Shed HER2 amount in plasma (umol)
DRP   : Drug-soluble plasma HER2 complex amount (umol)

$ODE
double CP  = AP / V1;
double CI  = AI / VI;
double CRM = RM / VI;
double CRS = RS / VI;
double CRP = RP / V1;

dxdt_AP = -(CL_lin / V1) * AP - Q * (CP - CI) - kon_p * CP * CRP * V1 + koff_p * DRP;

dxdt_AI = Q * (CP - CI) - kon_m * CI * CRM * VI + koff_m * DRM - kon_s * CI * CRS * VI + koff_s * DRS;

dxdt_RM = (ksyn_m - kdeg_m * CRM - ksh_p * CRM) * VI - kon_m * CI * CRM * VI + koff_m * DRM;

dxdt_DRM = kon_m * CI * CRM * VI - (koff_m + kint_m) * DRM;

dxdt_RS = ksh_p * RM - lambda * RS - kon_s * CI * CRS * VI + koff_s * DRS;

dxdt_DRS = kon_s * CI * CRS * VI - (koff_s + kelm_s) * DRS;

dxdt_RP = lambda * RS - kdeg_p * RP - kon_p * CP * CRP * V1 + koff_p * DRP;

dxdt_DRP = kon_p * CP * CRP * V1 - (koff_p + kelm_p) * DRP;

$TABLE
double CP_ugmL   = (AP / V1) * 146.0;
double CI_ugmL   = (AI / VI) * 146.0;
double ECD_ngmL  = (RP / V1) * 100000.0;

// Membrane-only receptor occupancy (original)
double RO_pct    = (DRM / (DRM + RM + 1e-20)) * 100.0;

// Total HER2 occupancy: all bound drug / all HER2 (membrane + shed IS + shed plasma)
double RO_total  = ((DRM + DRS + DRP) / (RM + DRM + RS + DRS + RP + DRP + 1e-20)) * 100.0;

$CAPTURE CP_ugmL CI_ugmL ECD_ngmL RO_pct RO_total
'

mod <- mcode("trastuzumab_pbpk_tmdd_v2", 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; }
  .metric-danger .metric-value { color: #ef4444; }
  .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; }
")

# ── Helper: run one simulation for a given regimen ────────────────────────────
run_sim <- function(mod, dp, n_cycles, wt, ecd_val, affinity_ratio) {
  ksh_p_val   <- ecd_val / 7241
  RM_ss       <- 0.042
  RS_ss       <- ksh_p_val * RM_ss / 0.01
  RP_ss       <- 0.01 * RS_ss / 0.2
  ksyn_m_val  <- (0.1 + ksh_p_val) * 0.0035
  koff_s_val  <- 1.26 * affinity_ratio

  dose_load_mg   <- dp$load_mgkg  * wt
  dose_maint_mg  <- dp$maint_mgkg * wt
  dose_load_umol  <- dose_load_mg  / 146
  dose_maint_umol <- dose_maint_mg / 146
  inf_rate_load  <- dose_load_umol  / dp$inf_load_h
  inf_rate_maint <- dose_maint_umol / dp$inf_maint_h
  interval_h     <- dp$interval_days * 24

  evs <- lapply(seq_len(n_cycles), function(i) {
    if (i == 1) {
      ev(amt = dose_load_umol,  cmt = 1, rate = inf_rate_load,  time = 0)
    } else {
      ev(amt = dose_maint_umol, cmt = 1, rate = inf_rate_maint, time = (i - 1) * interval_h)
    }
  })
  ev_all   <- do.call(c, evs)
  end_time <- n_cycles * interval_h

  mod %>%
    param(ksh_p = ksh_p_val, ksyn_m = ksyn_m_val, koff_s = koff_s_val) %>%
    init(AP = 0, AI = 0, RM = RM_ss, DRM = 0, RS = RS_ss, DRS = 0, RP = RP_ss, DRP = 0) %>%
    ev(ev_all) %>%
    mrgsim(end = end_time, delta = 1) %>%
    as.data.frame() %>%
    mutate(time_days = time / 24)
}

# ── All standard regimen definitions ─────────────────────────────────────────
all_regimens <- list(
  weekly = list(label = "Weekly (4→2 mg/kg)",  load_mgkg = 4, maint_mgkg = 2, interval_days = 7,  inf_load_h = 1.5, inf_maint_h = 0.5),
  q3w    = list(label = "Q3W (8→6 mg/kg)",     load_mgkg = 8, maint_mgkg = 6, interval_days = 21, inf_load_h = 1.5, inf_maint_h = 0.5),
  q4w    = list(label = "Q4W (8→6 mg/kg)",     load_mgkg = 8, maint_mgkg = 6, interval_days = 28, inf_load_h = 1.5, inf_maint_h = 0.5),
  q6w    = list(label = "Q6W (8→8 mg/kg)",     load_mgkg = 8, maint_mgkg = 8, interval_days = 42, inf_load_h = 1.5, inf_maint_h = 0.5)
)
comparison_colors <- c(
  weekly = "#8b5cf6",
  q3w    = "#10b981",
  q4w    = "#f59e0b",
  q6w    = "#ef4444"
)

# ── UI ────────────────────────────────────────────────────────────────────────
ui <- page_sidebar(
  title = "Trastuzumab Minimal PBPK-TMDD Target Shedding Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings", width = 340,
    h6("\U0001F489 Dosing Regimen"),
    selectInput("regimen", "Dosing Regimen",
      choices = c(
        "Standard weekly (4\u21922 mg/kg)"  = "weekly",
        "Every 3 weeks (8\u21926 mg/kg)"    = "q3w",
        "Every 4 weeks (8\u21926 mg/kg)"    = "q4w",
        "Every 6 weeks (8\u21928 mg/kg)"    = "q6w",
        "Custom"                             = "custom"
      )
    ),
    conditionalPanel(
      condition = "input.regimen == 'custom'",
      sliderInput("loading_dose", "Loading Dose (mg/kg)", min = 2, max = 8, value = 4, step = 1),
      sliderInput("maint_dose",   "Maintenance Dose (mg/kg)", min = 1, max = 8, value = 2, step = 1),
      sliderInput("interval",     "Dosing Interval (days)", min = 7, max = 42, value = 7, step = 7)
    ),
    numericInput("n_cycles", "Number of Cycles", value = 8, min = 1, max = 24),
    hr(), h6("\U0001F9EC Patient Parameters"),
    sliderInput("wt", "Body Weight (kg)", min = 40, max = 120, value = 70, step = 5),
    sliderInput("ecd_her2", "Plasma ECD HER2 (ng/mL)", min = 0, max = 2200, value = 15, step = 5),
    hr(), h6("\U0001F52C Antibody Engineering (Advanced)"),
    sliderInput("affinity_ratio", "Affinity Ratio (soluble/membrane)", min = 1, max = 1000, value = 1, step = 1),
    hr(), h6("\U0001F4CA Display Options"),
    radioButtons("ro_type", "Receptor Occupancy Method",
      choices = c(
        "Membrane HER2 only"             = "ro_membrane",
        "Total HER2 (membrane + soluble)" = "ro_total"
      ),
      selected = "ro_total"
    ),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE),
    hr(), h6("\U0001F504 Dose Comparison"),
    numericInput("comp_weeks", "Comparison Duration (weeks)", value = 24, min = 4, max = 52),
    checkboxGroupInput("comp_regimens", "Regimens to Compare",
      choices  = c("Weekly (4→2 mg/kg)" = "weekly",
                   "Q3W (8→6 mg/kg)"   = "q3w",
                   "Q4W (8→6 mg/kg)"   = "q4w",
                   "Q6W (8→8 mg/kg)"   = "q6w"),
      selected = c("weekly", "q3w", "q4w", "q6w")
    )
  ),
  navset_card_tab(
    title = "Trastuzumab PBPK-TMDD Simulator",
    full_screen = TRUE,

    # ── Tab 1: PK/PD Simulation ──────────────────────────────────────────────
    nav_panel("PK/PD Simulation",
      layout_columns(
        col_widths = c(3, 3, 3, 3),
        div(class = "metric-card metric-success",
          div(class = "metric-value", textOutput("trough_drug")),
          div(class = "metric-label", "Trough Drug (\u00b5g/mL)")
        ),
        div(class = "metric-card metric-primary",
          div(class = "metric-value", textOutput("target_occ")),
          div(class = "metric-label", textOutput("ro_label_metric"))
        ),
        div(class = "metric-card metric-info",
          div(class = "metric-value", textOutput("thalf")),
          div(class = "metric-label", "t\u00bd (days)")
        ),
        div(class = "metric-card metric-warning",
          div(class = "metric-value", textOutput("ecd_impact")),
          div(class = "metric-label", "ECD-HER2 Impact")
        )
      ),
      layout_columns(
        col_widths = c(6, 6),
        plotOutput("pkPlot", height = "420px"),
        plotOutput("roPlot", height = "420px")
      )
    ),

    # ── Tab 2: Dose Comparison ───────────────────────────────────────────────
    nav_panel("Dose Comparison",
      h5("All selected regimens simulated for the same patient over the same time horizon."),
      layout_columns(
        col_widths = c(6, 6),
        plotOutput("compPKPlot", height = "420px"),
        plotOutput("compROPlot", height = "420px")
      )
    ),

    # ── Tab 3: Model Information ─────────────────────────────────────────────
    nav_panel("Model Information", markdown("
## Trastuzumab (Herceptin\u00ae) \u2014 Minimal PBPK-TMDD with HER2 Target Shedding

**Model Type:** Minimal PBPK-TMDD with target shedding  
**Drug:** Trastuzumab (Herceptin\u00ae), anti-HER2 monoclonal antibody, MW = 146 kDa  
**Indication:** HER2-positive metastatic breast cancer  

### Mechanism
HER2 exists in two forms:  
1. **Membrane-bound (RM)** \u2014 the therapeutic target on tumor cells  
2. **Shed soluble ECD HER2** \u2014 cleaved extracellular domain circulating in interstitial fluid and plasma  

Trastuzumab binds both forms. Soluble ECD HER2 acts as a **drug sink**, competing with membrane-bound HER2 for drug binding and reducing target occupancy.

### Receptor Occupancy Methods
| Method | Formula | Interpretation |
|--------|---------|----------------|
| **Membrane-only** | DRM / (RM + DRM) | Fraction of membrane HER2 occupied |
| **Total HER2** | (DRM + DRS + DRP) / (RM + DRM + RS + DRS + RP + DRP) | Fraction of all HER2 forms (membrane + shed IS + shed plasma) occupied |

The **Total HER2 occupancy** provides a more complete picture when shed ECD HER2 is elevated, reflecting the drug's engagement across all HER2 pools.

### Dosing Regimens
| Regimen | Loading | Maintenance | Interval |
|---------|---------|-------------|----------|
| Standard Weekly | 4 mg/kg IV (90 min) | 2 mg/kg IV (30 min) | 7 days |
| Every 3 Weeks | 8 mg/kg IV (90 min) | 6 mg/kg IV (30 min) | 21 days |
| Every 4 Weeks | 8 mg/kg IV (90 min) | 6 mg/kg IV (30 min) | 28 days |
| Every 6 Weeks | 8 mg/kg IV (90 min) | 8 mg/kg IV (30 min) | 42 days |

*Q4W and Q6W are exploratory regimens; Q6W uses a higher maintenance dose to compensate for the extended dosing interval.*

### Key Clinical Insight
Patients with **ECD HER2 > 500 ng/mL** have drastically reduced trough concentrations and lower receptor occupancy, potentially explaining reduced clinical efficacy in this subgroup.

### Model Parameters
| Parameter | Value | Units | Description |
|-----------|-------|-------|-------------|
| V1 | 2.9 | L | Plasma volume |
| VI | 12.0 | L | Interstitial volume |
| CL_lin | 0.00938 | L/h | Linear clearance |
| Q | 0.05 | L/h | Intercompartmental CL |
| kon_m | 2520 | (\u00b5mol/L)\u207b\u00b9\u00b7h\u207b\u00b9 | Association rate (membrane HER2) |
| koff_m | 1.26 | h\u207b\u00b9 | Dissociation rate (membrane HER2) |
| kint_m | 0.10 | h\u207b\u00b9 | Internalization rate (drug-HER2 complex) |
| kdeg_m | 0.10 | h\u207b\u00b9 | Membrane HER2 degradation |
| [RM]ss | 0.0035 | \u00b5mol/L | Baseline membrane HER2 (IS) |
| \u03bb | 0.01 | h\u207b\u00b9 | Lymphatic transport (IS \u2192 plasma) |
| kdeg_p | 0.20 | h\u207b\u00b9 | Plasma ECD degradation |
| kelm_p | 0.20 | h\u207b\u00b9 | Plasma complex elimination |
    ")),

    # ── Tab 4: References ────────────────────────────────────────────────────
    nav_panel("References", div(class = "ref-box",
      tags$h5("\U0001F4DA Key References"),
      tags$ol(
        tags$li("Betts AM et al. (2013) Incorporating Target Shedding Into a Minimal PBPK-TMDD Model for Monoclonal Antibodies. ",
                tags$em("CPT Pharmacometrics Syst Pharmacol"), " 2:e30. doi:10.1038/psp.2013.6"),
        tags$li("Baselga J et al. (2010) Relationship between tumor biomarkers and efficacy in EMILIA. ",
                tags$em("J Clin Oncol"), " 28(16)."),
        tags$li("Mager DE, Jusko WJ (2001) General pharmacokinetic model for drugs exhibiting target-mediated drug disposition. ",
                tags$em("J Pharmacokinet Pharmacodyn"), " 28:507\u2013532."),
        tags$li("Bruno R et al. (2005) Population pharmacokinetics of trastuzumab in patients with HER2+ breast cancer. ",
                tags$em("J Clin Pharmacol"), " 45:286\u2013292."),
        tags$li("Gibiansky L et al. (2008) Approximations of the target-mediated drug disposition model and identifiability of model parameters. ",
                tags$em("J Pharmacokinet Pharmacodyn"), " 35:573\u2013591."),
        tags$li("Shah DK, Betts AM (2012) Towards a platform PBPK model to characterize the plasma and tissue disposition of monoclonal antibodies in preclinical species and human. ",
                tags$em("J Pharmacokinet Pharmacodyn"), " 39:67\u201386.")
      ),
      tags$h5("\U0001F48A Therapeutic Context"),
      tags$ul(
        tags$li(tags$strong("Class:"), " Humanized monoclonal antibody (IgG1)"),
        tags$li(tags$strong("Target:"), " HER2 (ErbB2) receptor"),
        tags$li(tags$strong("Indication:"), " HER2-positive metastatic breast cancer"),
        tags$li(tags$strong("Mechanism:"), " Inhibits HER2 signaling, induces ADCC, blocks receptor shedding"),
        tags$li(tags$strong("Key finding:"), " High plasma ECD HER2 (>500 ng/mL) reduces drug exposure and target occupancy")
      )
    ))
  ),
  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"),
    " \u2022 Built by Sunny \u2600\uFE0F (Husain Attarwala's AI Assistant)", br(),
    tags$span(style = "font-size: 10px;", "For research and educational purposes only. Not for clinical decision-making."))
)

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

  get_dosing_params <- reactive({
    switch(input$regimen,
      weekly = all_regimens$weekly,
      q3w    = all_regimens$q3w,
      q4w    = all_regimens$q4w,
      q6w    = all_regimens$q6w,
      custom = list(
        load_mgkg    = input$loading_dose,
        maint_mgkg   = input$maint_dose,
        interval_days= input$interval,
        inf_load_h   = 1.5,
        inf_maint_h  = 0.5
      )
    )
  })

  sim_data <- reactive({
    shiny::req(input$wt, input$n_cycles)
    run_sim(mod, get_dosing_params(), input$n_cycles, input$wt, input$ecd_her2, input$affinity_ratio)
  })

  sim_data_zero_ecd <- reactive({
    shiny::req(input$wt, input$n_cycles)
    run_sim(mod, get_dosing_params(), input$n_cycles, input$wt, 0, input$affinity_ratio)
  })

  # ── Comparison simulations ─────────────────────────────────────────────────
  comp_data <- reactive({
    shiny::req(input$comp_weeks, input$comp_regimens, input$wt)
    comp_weeks <- max(4, input$comp_weeks)
    comp_days  <- comp_weeks * 7

    bind_rows(lapply(input$comp_regimens, function(rname) {
      dp <- all_regimens[[rname]]
      n  <- ceiling(comp_days / dp$interval_days)
      n  <- max(1, n)
      d  <- run_sim(mod, dp, n, input$wt, input$ecd_her2, input$affinity_ratio)
      d  <- d[d$time_days <= comp_days, ]
      d$regimen <- all_regimens[[rname]]$label
      d$regimen_key <- rname
      d
    }))
  })

  # ── Metrics ────────────────────────────────────────────────────────────────
  output$trough_drug <- renderText({
    d  <- sim_data()
    dp <- get_dosing_params()
    interval_h    <- dp$interval_days * 24
    last_dose_time <- (input$n_cycles - 1) * interval_h
    tw <- d[d$time >= last_dose_time + 2 & d$time <= last_dose_time + interval_h, ]
    if (nrow(tw) == 0) return("--")
    sprintf("%.1f", min(tw$CP_ugmL, na.rm = TRUE))
  })

  output$ro_label_metric <- renderText({
    if (input$ro_type == "ro_total") "% Total HER2 Occupancy" else "% Membrane HER2 Occupancy"
  })

  output$target_occ <- renderText({
    d  <- sim_data()
    dp <- get_dosing_params()
    interval_h    <- dp$interval_days * 24
    last_dose_time <- (input$n_cycles - 1) * interval_h
    lc <- d[d$time >= last_dose_time & d$time <= last_dose_time + interval_h, ]
    if (nrow(lc) == 0) return("--")
    col_name <- if (input$ro_type == "ro_total") "RO_total" else "RO_pct"
    sprintf("%.1f%%", mean(lc[[col_name]], na.rm = TRUE))
  })

  output$thalf <- renderText({
    cl_lin    <- 0.00938
    vss       <- 2.9 + 12.0
    thalf_days <- log(2) / (cl_lin / vss) / 24
    sprintf("%.1f", thalf_days)
  })

  output$ecd_impact <- renderText({
    d  <- sim_data()
    d0 <- sim_data_zero_ecd()
    dp <- get_dosing_params()
    interval_h    <- dp$interval_days * 24
    last_dose_time <- (input$n_cycles - 1) * interval_h
    tw  <- d[ d$time  >= last_dose_time + 2 & d$time  <= last_dose_time + interval_h, ]
    tw0 <- d0[d0$time >= last_dose_time + 2 & d0$time <= last_dose_time + interval_h, ]
    if (nrow(tw) == 0 || nrow(tw0) == 0) return("--")
    trough_ecd  <- min(tw$CP_ugmL,  na.rm = TRUE)
    trough_zero <- min(tw0$CP_ugmL, na.rm = TRUE)
    if (trough_zero < 0.01) return("--")
    red <- max(0, (1 - trough_ecd / trough_zero) * 100)
    sprintf("-%.0f%%", red)
  })

  # ── PK plot ────────────────────────────────────────────────────────────────
  output$pkPlot <- renderPlot({
    d <- sim_data()
    p <- ggplot(d, aes(x = time_days, y = CP_ugmL)) +
      geom_hline(yintercept = 10, linetype = "dashed", color = "#10b981", alpha = 0.7) +
      annotate("text", x = max(d$time_days) * 0.02, y = 12,
               label = "Therapeutic target: 10 \u00b5g/mL", hjust = 0, size = 3.5, color = "#10b981") +
      geom_line(color = "#8b5cf6", linewidth = 0.8) +
      labs(x = "Time (days)", y = "Plasma Trastuzumab (\u00b5g/mL)",
           title = paste0("Trastuzumab PK | ECD HER2 = ", input$ecd_her2, " ng/mL")) +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(face = "bold", size = 14))
    if (input$log_scale) p <- p + scale_y_log10()
    p
  })

  # ── RO plot ────────────────────────────────────────────────────────────────
  output$roPlot <- renderPlot({
    d        <- sim_data()
    col_name <- if (input$ro_type == "ro_total") "RO_total" else "RO_pct"
    y_label  <- if (input$ro_type == "ro_total") "Total HER2 Receptor Occupancy (%)"
                else "Membrane HER2 Receptor Occupancy (%)"
    plot_title <- if (input$ro_type == "ro_total")
      "Total HER2 Receptor Occupancy (Membrane + Shed ECD)"
    else
      "Membrane HER2 Receptor Occupancy"

    ggplot(d, aes(x = time_days, y = .data[[col_name]])) +
      geom_line(color = "#f59e0b", linewidth = 0.8) +
      ylim(0, 100) +
      labs(x = "Time (days)", y = y_label, title = plot_title) +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(face = "bold", size = 14))
  })

  # ── Comparison PK plot ─────────────────────────────────────────────────────
  output$compPKPlot <- renderPlot({
    d <- comp_data()
    shiny::req(nrow(d) > 0)

    reg_levels  <- unique(d$regimen)
    reg_keys    <- unique(d$regimen_key)
    color_map   <- setNames(comparison_colors[reg_keys], reg_levels)

    p <- ggplot(d, aes(x = time_days, y = CP_ugmL, color = regimen)) +
      geom_hline(yintercept = 10, linetype = "dashed", color = "grey60", alpha = 0.8) +
      annotate("text", x = 0, y = 11.5, label = "10 \u00b5g/mL target",
               hjust = 0, size = 3.2, color = "grey50") +
      geom_line(linewidth = 0.9) +
      scale_color_manual(values = color_map) +
      labs(
        x = "Time (days)",
        y = "Plasma Trastuzumab (\u00b5g/mL)",
        color = "Regimen",
        title = paste0("PK Comparison — All Regimens | ECD HER2 = ", input$ecd_her2,
                       " ng/mL | BW = ", input$wt, " kg")
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title   = element_text(face = "bold", size = 13),
        legend.position = "bottom"
      )
    if (input$log_scale) p <- p + scale_y_log10()
    p
  })

  # ── Comparison RO plot ─────────────────────────────────────────────────────
  output$compROPlot <- renderPlot({
    d        <- comp_data()
    shiny::req(nrow(d) > 0)
    col_name  <- if (input$ro_type == "ro_total") "RO_total" else "RO_pct"
    y_label   <- if (input$ro_type == "ro_total") "Total HER2 Occupancy (%)"
                 else "Membrane HER2 Occupancy (%)"

    reg_levels <- unique(d$regimen)
    reg_keys   <- unique(d$regimen_key)
    color_map  <- setNames(comparison_colors[reg_keys], reg_levels)

    ggplot(d, aes(x = time_days, y = .data[[col_name]], color = regimen)) +
      geom_line(linewidth = 0.9) +
      scale_color_manual(values = color_map) +
      ylim(0, 100) +
      labs(
        x = "Time (days)",
        y = y_label,
        color = "Regimen",
        title = paste0("Receptor Occupancy Comparison — ", y_label)
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title      = element_text(face = "bold", size = 13),
        legend.position = "bottom"
      )
  })
}

shinyApp(ui = ui, server = server)