suppressPackageStartupMessages({
  library(shiny)
  library(bslib)
  library(mrgsolve)
  library(dplyr)
  library(ggplot2)
})

model_code <- '
$PARAM @annotated
CL      : 0.00649   : Linear clearance at reference weight (L/h)
VC      : 3.04      : Central volume at reference weight (L)
Q       : 0.0207    : Intercompartmental clearance (L/h)
VP      : 1.74      : Peripheral volume at reference weight (L)
KA      : 0.0126    : SC absorption rate constant (1/h)
F1      : 0.566     : SC bioavailability
R0      : 0.00184   : Baseline total target concentration (umol/L; 1.84 nM)
KSS     : 0.0000296 : Quasi-steady-state constant (umol/L; 0.0296 nM)
KINT    : 0.000141  : Complex internalization/elimination rate (1/h)
KDEG    : 0.097     : Target degradation rate constant (1/h)
WTREF   : 82.1      : Reference weight (kg)
WT      : 82.1      : Body weight (kg)
EXPCL   : 0.943     : Weight exponent on CL
EXPVC   : 0.762     : Weight exponent on Vc
EXPVP   : 0.646     : Weight exponent on Vp
MW_KDA  : 150       : Antibody molecular weight assumption (kDa)

$CMT @annotated
DEPOT  : SC depot (umol)
CENT   : Total drug in central compartment (umol)
PERIPH : Drug in peripheral compartment (umol)
RTOT   : Total target concentration in central compartment (umol/L)

$MAIN
F_DEPOT = F1;
RTOT_0  = R0;

$ODE
double cli_ode = CL * pow(WT / WTREF, EXPCL);
double vci_ode = VC * pow(WT / WTREF, EXPVC);
double vpi_ode = VP * pow(WT / WTREF, EXPVP);
double ksyn_ode = KDEG * R0;
double ctot_ode = CENT / vci_ode;
double disc_ode = pow(ctot_ode - RTOT - KSS, 2.0) + 4.0 * KSS * ctot_ode;
double rad_ode = disc_ode > 0.0 ? disc_ode : 0.0;
double cfree_ode_raw = 0.5 * ((ctot_ode - RTOT - KSS) + sqrt(rad_ode));
double cfree_ode = cfree_ode_raw > 0.0 ? cfree_ode_raw : 0.0;
double cbound_ode = RTOT * cfree_ode / (KSS + cfree_ode + 1e-12);

// Drug amounts are tracked in umol; RTOT/cfree_ode/cbound_ode are concentrations in umol/L.
dxdt_DEPOT  = -KA * DEPOT;
dxdt_CENT   = KA * DEPOT - cli_ode * ctot_ode - Q * cfree_ode + (Q / vpi_ode) * PERIPH - KINT * cbound_ode * vci_ode;
dxdt_PERIPH = Q * cfree_ode - (Q / vpi_ode) * PERIPH;
dxdt_RTOT   = ksyn_ode - KDEG * RTOT - (KINT - KDEG) * cbound_ode;

$TABLE
double cli_tbl = CL * pow(WT / WTREF, EXPCL);
double vci_tbl = VC * pow(WT / WTREF, EXPVC);
double vpi_tbl = VP * pow(WT / WTREF, EXPVP);
double ctot_tbl = CENT / vci_tbl;
double disc_tbl = pow(ctot_tbl - RTOT - KSS, 2.0) + 4.0 * KSS * ctot_tbl;
double rad_tbl = disc_tbl > 0.0 ? disc_tbl : 0.0;
double cfree_tbl_raw = 0.5 * ((ctot_tbl - RTOT - KSS) + sqrt(rad_tbl));
double cfree_tbl = cfree_tbl_raw > 0.0 ? cfree_tbl_raw : 0.0;
double rfree_tbl = RTOT * KSS / (KSS + cfree_tbl + 1e-12);
double cbound_tbl = RTOT - rfree_tbl;
double CP_total = ctot_tbl * MW_KDA;
double CP_free  = cfree_tbl * MW_KDA;
double RO_pct   = RTOT > 1e-12 ? 100.0 * cbound_tbl / RTOT : 0.0;
double RTOTout  = RTOT;
double Rfreeout = rfree_tbl;
double CLi = cli_tbl;
double VCi = vci_tbl;
double VPi = vpi_tbl;
double Ctot = ctot_tbl;
double Cfree = cfree_tbl;

$CAPTURE @annotated
CP_total : Total plasma concentration (ug/mL)
CP_free  : Free plasma concentration (ug/mL)
RO_pct   : Central target occupancy (%)
Ctot     : Total drug concentration in central compartment (umol/L)
Cfree    : Free drug concentration in central compartment (umol/L)
RTOTout  : Total target concentration in central compartment (umol/L)
Rfreeout : Free target concentration in central compartment (umol/L)
CLi      : Clearance adjusted for weight (L/h)
VCi      : Central volume adjusted for weight (L)
VPi      : Peripheral volume adjusted for weight (L)
'

mod <- mcode("vixarelimab_tmdd_qss", 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; }
  .small-note {
    font-size: 12px;
    color: #6c757d;
    line-height: 1.4;
  }
")

mg_to_umol <- function(mg, mw_kda = 150) {
  mg / mw_kda
}

interval_hours <- function(interval_choice) {
  if (identical(interval_choice, "single")) {
    return(0)
  }
  as.numeric(interval_choice)
}

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

build_vixarelimab_events <- function(
  route = c("SC", "IV"),
  sc_dose_mg = 360,
  iv_dose_mgkg = 5,
  wt_kg = 82.1,
  interval_choice = "168",
  n_doses = 8,
  use_loading = TRUE,
  loading_mg = 720,
  infusion_hours = 1,
  followup_weeks = 8,
  mw_kda = 150
) {
  route <- match.arg(route)
  n_doses <- max(1L, as.integer(round(n_doses)))
  ii <- interval_hours(interval_choice)
  followup_h <- max(24, followup_weeks * 7 * 24)

  if (route == "SC") {
    maint_dose_mg <- as.numeric(sc_dose_mg)
    maint_amt <- mg_to_umol(maint_dose_mg, mw_kda = mw_kda)
    loading_amt <- mg_to_umol(as.numeric(loading_mg), mw_kda = mw_kda)
    cmt <- 1

    if (ii <= 0 || n_doses == 1) {
      events <- ev(time = 0, amt = maint_amt, cmt = cmt)
      dose_times_h <- 0
      regimen_label <- sprintf("SC %.0f mg single dose", maint_dose_mg)
    } else if (isTRUE(use_loading)) {
      events <- ev(time = 0, amt = loading_amt, cmt = cmt)
      if (n_doses > 1) {
        maint_ev <- ev(time = ii, amt = maint_amt, cmt = cmt, ii = ii, addl = n_doses - 2)
        events <- c(events, maint_ev)
      }
      dose_times_h <- c(0, if (n_doses > 1) seq(ii, by = ii, length.out = n_doses - 1) else numeric(0))
      regimen_label <- sprintf(
        "SC 720 mg load → %.0f mg every %s",
        maint_dose_mg,
        dplyr::case_when(
          ii == 168 ~ "week",
          ii == 336 ~ "2 weeks",
          ii == 672 ~ "4 weeks",
          TRUE ~ paste(ii, "h")
        )
      )
    } else {
      events <- ev(time = 0, amt = maint_amt, cmt = cmt, ii = ii, addl = n_doses - 1)
      dose_times_h <- seq(0, by = ii, length.out = n_doses)
      regimen_label <- sprintf(
        "SC %.0f mg every %s",
        maint_dose_mg,
        dplyr::case_when(
          ii == 168 ~ "week",
          ii == 336 ~ "2 weeks",
          ii == 672 ~ "4 weeks",
          TRUE ~ paste(ii, "h")
        )
      )
    }
  } else {
    dose_mg <- as.numeric(iv_dose_mgkg) * wt_kg
    dose_amt <- mg_to_umol(dose_mg, mw_kda = mw_kda)
    cmt <- 2

    if (ii <= 0 || n_doses == 1) {
      events <- ev(time = 0, amt = dose_amt, rate = dose_amt / infusion_hours, cmt = cmt)
      dose_times_h <- 0
      regimen_label <- sprintf("IV %.1f mg/kg single dose", as.numeric(iv_dose_mgkg))
    } else {
      events <- ev(
        time = 0,
        amt = dose_amt,
        rate = dose_amt / infusion_hours,
        cmt = cmt,
        ii = ii,
        addl = n_doses - 1
      )
      dose_times_h <- seq(0, by = ii, length.out = n_doses)
      regimen_label <- sprintf(
        "IV %.1f mg/kg every %s",
        as.numeric(iv_dose_mgkg),
        dplyr::case_when(
          ii == 168 ~ "week",
          ii == 336 ~ "2 weeks",
          ii == 672 ~ "4 weeks",
          TRUE ~ paste(ii, "h")
        )
      )
    }
  }

  last_dose_time <- max(dose_times_h)
  sim_end <- last_dose_time + max(followup_h, if (ii > 0) ii else 0)

  list(
    events = events,
    dose_times_h = dose_times_h,
    ii = ii,
    n_doses = n_doses,
    last_dose_time = last_dose_time,
    sim_end = sim_end,
    regimen_label = regimen_label
  )
}

simulate_vixarelimab <- function(
  route = c("SC", "IV"),
  sc_dose_mg = 360,
  iv_dose_mgkg = 5,
  wt_kg = 82.1,
  interval_choice = "168",
  n_doses = 8,
  use_loading = TRUE,
  loading_mg = 720,
  followup_weeks = 8,
  delta = 6
) {
  route <- match.arg(route)

  schedule <- build_vixarelimab_events(
    route = route,
    sc_dose_mg = sc_dose_mg,
    iv_dose_mgkg = iv_dose_mgkg,
    wt_kg = wt_kg,
    interval_choice = interval_choice,
    n_doses = n_doses,
    use_loading = use_loading,
    loading_mg = loading_mg,
    followup_weeks = followup_weeks
  )

  out <- mod %>%
    param(WT = wt_kg) %>%
    ev(schedule$events) %>%
    mrgsim(end = schedule$sim_end, delta = delta, rtol = 1e-8, atol = 1e-10) %>%
    as.data.frame() |>
    mutate(
      time_days = time / 24,
      time_weeks = time / (24 * 7)
    )

  list(data = out, schedule = schedule)
}

compute_linear_half_life_days <- function(wt_kg = 82.1) {
  cli <- 0.00649 * (wt_kg / 82.1)^0.943
  vci <- 3.04 * (wt_kg / 82.1)^0.762
  vpi <- 1.74 * (wt_kg / 82.1)^0.646
  q <- 0.0207

  k10 <- cli / vci
  k12 <- q / vci
  k21 <- q / vpi

  disc <- (k10 + k12 + k21)^2 - 4 * k10 * k21
  beta <- 0.5 * ((k10 + k12 + k21) - sqrt(pmax(disc, 0)))

  if (!is.finite(beta) || beta <= 0) {
    return(NA_real_)
  }

  log(2) / beta / 24
}

compute_exposure_metrics <- function(sim_df, schedule) {
  if (is.null(sim_df) || nrow(sim_df) == 0) {
    return(list(
      cmax = NA_real_,
      ctrough = NA_real_,
      auc = NA_real_,
      ro_trough = NA_real_
    ))
  }

  if (schedule$ii > 0 && schedule$n_doses > 1) {
    interval_end <- schedule$last_dose_time + schedule$ii
    interval_df <- sim_df |>
      dplyr::filter(time >= schedule$last_dose_time, time <= interval_end)

    trough_df <- sim_df |>
      dplyr::filter(time <= interval_end) |>
      dplyr::filter(time == max(time, na.rm = TRUE))
  } else {
    interval_df <- sim_df
    trough_df <- sim_df |>
      dplyr::filter(time == max(time, na.rm = TRUE))
  }

  if (nrow(interval_df) < 2) {
    interval_df <- sim_df
  }

  list(
    cmax = max(interval_df$CP_total, na.rm = TRUE),
    ctrough = trough_df$CP_total[[1]],
    auc = trap_auc(interval_df$time, interval_df$CP_total),
    ro_trough = trough_df$RO_pct[[1]]
  )
}

model_info_md <- paste(
  "## Vixarelimab TMDD-QSS PK Simulator",
  "",
  "**Model structure**  ",
  "Two-compartment PK with a quasi-steady-state (QSS) TMDD approximation in the **central compartment**. Drug is tracked as total amount in the depot, central, and peripheral compartments, while the target is tracked as **total central target concentration**.",
  "",
  "**Core equations**  ",
  "- dA_DEPOT/dt = -Ka × A_DEPOT  ",
  "- dA_CENT/dt = Ka × A_DEPOT − CL × C_free − Q × C_free + Q/Vp × A_PERIPH − Kint × C_bound × Vc  ",
  "- dA_PERIPH/dt = Q × C_free − Q/Vp × A_PERIPH  ",
  "- dR_tot/dt = Kdeg × R0 − Kdeg × R_tot − (Kint − Kdeg) × C_bound",
  "",
  "with the standard QSS free-drug solution in the central compartment:",
  "",
  "C_free = 0.5 × ((C_tot − R_tot − Kss) + sqrt((C_tot − R_tot − Kss)^2 + 4 × Kss × C_tot))",
  "",
  "**Published parameterization used for this build**  ",
  "- CL = 0.00649 L/h  ",
  "- Vc = 3.04 L  ",
  "- Q = 0.0207 L/h  ",
  "- Vp = 1.74 L  ",
  "- Ka = 0.0126 h^-1  ",
  "- F1 = 0.566  ",
  "- R0 = 0.00184 umol/L (1.84 nM)  ",
  "- Kss = 0.0000296 umol/L (0.0296 nM)  ",
  "- Kint = 0.000141 h^-1  ",
  "- Kdeg = 0.097 h^-1",
  "",
  "**Body-weight scaling**  ",
  "Reference weight is 82.1 kg. Weight affects CL (0.943), Vc (0.762), and Vp (0.646).",
  "",
  "**Important approximation / caveat**  ",
  "This is a **single-target central TMDD-QSS approximation**. Receptor occupancy is reported for the central target pool and does **not** distinguish soluble versus membrane-bound OSMR in tissue. That matches the paper's own caveat that central occupancy is only a proxy for tissue target engagement.",
  sep = "\n"
)

ui <- page_sidebar(
  title = "Vixarelimab PK Simulator",
  theme = app_theme,

  sidebar = sidebar(
    title = "Simulation Settings",
    width = 360,

    accordion(
      open = c("Regimen", "Patient"),
      accordion_panel(
        "Regimen",
        icon = bsicons::bs_icon("capsule"),
        radioButtons(
          "route",
          "Administration route",
          choices = c("Subcutaneous (SC)" = "SC", "Intravenous (IV)" = "IV"),
          selected = "SC"
        ),
        conditionalPanel(
          condition = "input.route == 'SC'",
          selectInput(
            "sc_dose_mg",
            "Maintenance SC dose (mg)",
            choices = c("120 mg" = 120, "360 mg" = 360, "540 mg" = 540, "720 mg" = 720),
            selected = 360
          ),
          checkboxInput("use_loading", "Use 720 mg loading dose", value = TRUE)
        ),
        conditionalPanel(
          condition = "input.route == 'IV'",
          selectInput(
            "iv_dose_mgkg",
            "IV dose (mg/kg)",
            choices = c("0.3" = 0.3, "1.5" = 1.5, "5" = 5, "7.5" = 7.5, "10" = 10, "20" = 20),
            selected = 5
          )
        ),
        selectInput(
          "interval_choice",
          "Dosing interval",
          choices = c(
            "Single dose" = "single",
            "Weekly (Q1W)" = "168",
            "Every 2 weeks (Q2W)" = "336",
            "Every 4 weeks (Q4W)" = "672"
          ),
          selected = "168"
        ),
        numericInput("n_doses", "Number of doses", value = 8, min = 1, max = 24),
        sliderInput("followup_weeks", "Follow-up after last dose (weeks)", min = 4, max = 24, value = 8, step = 1),
        div(
          class = "small-note",
          "Defaults reflect the clinically relevant regimen: 720 mg SC loading dose followed by 360 mg weekly maintenance."
        )
      ),
      accordion_panel(
        "Patient",
        icon = bsicons::bs_icon("person"),
        sliderInput("wt_kg", "Body weight (kg)", min = 40, max = 140, value = 82.1, step = 0.5),
        card(
          card_header("Reference covariates", class = "bg-light"),
          p(strong("Reference weight:"), "82.1 kg"),
          p(strong("Weight effects:"), "CL 0.943, Vc 0.762, Vp 0.646"),
          p(class = "small-note", "Model-based simulations in the paper indicated only a minor clinical effect of body weight.")
        )
      ),
      accordion_panel(
        "Display",
        icon = bsicons::bs_icon("gear"),
        checkboxInput("show_free", "Overlay free concentration", value = TRUE),
        checkboxInput("log_scale", "Log scale (Y-axis)", value = FALSE),
        checkboxInput("show_dose_lines", "Show dose-time markers", value = TRUE),
        div(
          class = "small-note",
          "Receptor occupancy is a central-compartment proxy, not a tissue membrane occupancy estimate."
        )
      )
    )
  ),

  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 (ug/mL)")
    ),
    div(
      class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("ctrough")),
      div(class = "metric-label", "Ctrough / final (ug/mL)")
    ),
    div(
      class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("auc")),
      div(class = "metric-label", "AUC (last interval, ug·h/mL)")
    ),
    div(
      class = "metric-card metric-info",
      div(class = "metric-value", textOutput("thalf")),
      div(class = "metric-label", "Approx. terminal t1/2 (days)")
    ),
    div(
      class = "metric-card metric-danger",
      div(class = "metric-value", textOutput("ro_trough")),
      div(class = "metric-label", "Target occupancy at trough (%)")
    )
  ),

  card(
    card_header(textOutput("plot_title")),
    full_screen = TRUE,
    plotOutput("pkPlot", height = "500px")
  ),

  navset_card_underline(
    title = "Drug Information",
    nav_panel(
      "Target Engagement",
      plotOutput("roPlot", height = "500px"),
      div(
        class = "small-note",
        style = "padding: 0 15px 15px 15px;",
        "This occupancy trace is calculated from the central QSS target pool and is best interpreted as a mechanistic proxy rather than direct tissue occupancy."
      )
    ),
    nav_panel("Model Information", markdown(model_info_md)),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Key References"),
        tags$ol(
          tags$li(
            "Tang F, Dang S, Arjomandi A, et al. (2026) ",
            tags$em("Population Pharmacokinetic Analysis of Vixarelimab in Healthy Volunteers and Patients With Chronic Pruritic Conditions"),
            ". CPT: Pharmacometrics & Systems Pharmacology 15:e70230. ",
            tags$a(href = "https://doi.org/10.1002/psp4.70230", target = "_blank", "doi:10.1002/psp4.70230")
          ),
          tags$li(
            "Gibiansky L, Gibiansky E, Kakkar T, Ma P. (2008) ",
            tags$em("Approximations of the target-mediated drug disposition model and identifiability of model parameters"),
            ". J Pharmacokinet Pharmacodyn 35:573-591."
          ),
          tags$li(
            "Sofen H, Bissonnette R, Yosipovitch G, et al. (2023) ",
            tags$em("Efficacy and Safety of Vixarelimab in Moderate-to-Severe Prurigo Nodularis"),
            ". EClinicalMedicine 57:101826."
          )
        ),
        tags$h5("Clinical context captured in this app"),
        tags$ul(
          tags$li("Single IV doses: 0.3, 1.5, 5, 7.5, 10, 20 mg/kg"),
          tags$li("Single SC: 1.5 mg/kg"),
          tags$li("Repeated SC: 360 mg QW"),
          tags$li("Repeated SC: 720 mg load + 360 mg QW"),
          tags$li("Repeated SC: 120 mg Q4W, 360 mg Q4W, 540 mg Q4W"),
          tags$li("Open-label / continuation context: 360 mg Q2W")
        ),
        tags$h5("Interpretation notes"),
        tags$ul(
          tags$li("Total and free plasma concentrations are reported in ug/mL."),
          tags$li("A molecular weight of 150 kDa was assumed for dose-to-amount conversion from mg to umol."),
          tags$li("TMDD parameters are implemented directly from Table 3, with receptor terms represented in umol/L for numerical consistency."),
          tags$li("120 mg SC Q4W is expected to underperform for sustained target engagement.")
        )
      )
    )
  ),

  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) {
  sim_bundle <- reactive({
    simulate_vixarelimab(
      route = input$route,
      sc_dose_mg = as.numeric(input$sc_dose_mg),
      iv_dose_mgkg = as.numeric(input$iv_dose_mgkg),
      wt_kg = input$wt_kg,
      interval_choice = input$interval_choice,
      n_doses = input$n_doses,
      use_loading = isTRUE(input$route == "SC") && isTRUE(input$use_loading) && input$interval_choice != "single",
      followup_weeks = input$followup_weeks
    )
  })

  sim_data <- reactive({
    sim_bundle()$data
  })

  schedule <- reactive({
    sim_bundle()$schedule
  })

  exposure_metrics <- reactive({
    compute_exposure_metrics(sim_data(), schedule())
  })

  output$plot_title <- renderText({
    schedule()$regimen_label
  })

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

  output$ctrough <- renderText({
    metrics <- exposure_metrics()
    sprintf("%.1f", metrics$ctrough)
  })

  output$auc <- renderText({
    metrics <- exposure_metrics()
    if (!is.finite(metrics$auc)) {
      return("—")
    }
    sprintf("%.0f", metrics$auc)
  })

  output$thalf <- renderText({
    thalf_days <- compute_linear_half_life_days(input$wt_kg)
    if (!is.finite(thalf_days)) {
      return("—")
    }
    sprintf("%.1f", thalf_days)
  })

  output$ro_trough <- renderText({
    metrics <- exposure_metrics()
    sprintf("%.1f", metrics$ro_trough)
  })

  output$pkPlot <- renderPlot({
    d <- sim_data()
    shiny::req(nrow(d) > 0)

    d_plot <- if (input$log_scale) {
      d |>
        dplyr::filter(CP_total > 0.01) |>
        dplyr::filter(!is.na(CP_total))
    } else {
      d
    }

    p <- ggplot(d_plot, aes(x = time_days))

    if (isTRUE(input$show_dose_lines)) {
      dose_df <- data.frame(time_days = schedule()$dose_times_h / 24)
      p <- p + geom_vline(
        data = dose_df,
        aes(xintercept = time_days),
        inherit.aes = FALSE,
        color = "#d1d5db",
        linewidth = 0.5,
        linetype = "dashed"
      )
    }

    p <- p +
      geom_line(aes(y = CP_total, color = "Total"), linewidth = 1.0)

    if (isTRUE(input$show_free)) {
      p <- p +
        geom_line(aes(y = CP_free, color = "Free"), linewidth = 0.9, linetype = "22")
    }

    p <- p +
      scale_color_manual(values = c("Total" = "#8b5cf6", "Free" = "#f59e0b")) +
      labs(
        x = "Time (days)",
        y = "Vixarelimab concentration (ug/mL)",
        color = NULL,
        title = schedule()$regimen_label,
        subtitle = sprintf("Weight %.1f kg • Central TMDD-QSS approximation", input$wt_kg)
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(face = "bold"),
        legend.position = "top"
      )

    if (input$log_scale) {
      p <- p + scale_y_log10(labels = scales::label_number())
    }

    p
  })

  output$roPlot <- renderPlot({
    d <- sim_data()
    shiny::req(nrow(d) > 0)

    ggplot(d, aes(x = time_days, y = RO_pct)) +
      geom_line(color = "#10b981", linewidth = 1.0) +
      scale_y_continuous(limits = c(0, 100)) +
      labs(
        x = "Time (days)",
        y = "Central target occupancy (%)",
        title = "Predicted target occupancy over time",
        subtitle = "Central target pool only; soluble vs membrane target not separated"
      ) +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(face = "bold"))
  })
}

shinyApp(ui = ui, server = server)