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

model_code <- '
$PARAM @annotated
CL   : 6.56  : Apparent clearance CL/F (L/h)
V1   : 49.3  : Apparent central volume V1/F (L)
V2   : 30.7  : Apparent peripheral volume 1 V2/F (L)
V3   : 37.1  : Apparent peripheral volume 2 V3/F (L)
Q1   : 3.52  : Intercompartmental clearance Q1/F (L/h)
Q2   : 0.769 : Intercompartmental clearance Q2/F (L/h)
KA   : 1.02  : Absorption rate constant (1/h)
D1   : 1.22  : Duration of zero-order absorption (h)
FR_ZO : 0.5  : Fraction absorbed via zero-order
FREL : 0.896 : Relative bioavailability capsule vs tablet
WT   : 75    : Body weight (kg)
POP  : 0     : Population (0=cancer, 1=healthy)
IND  : 0     : CYP3A4 inducer coadministration
INH  : 0     : CYP3A4 inhibitor coadministration
ALBL : 0     : Low albumin <30 g/L
ALPH : 0     : Elevated ALP > ULN
ETA_CL : 0   : IIV on CL
ETA_V1 : 0   : IIV on V1
ETA_V2 : 0   : IIV on V2
ETA_V3 : 0   : IIV on V3
ETA_KA : 0   : IIV on Ka
ETA_D1 : 0   : IIV on D1
ETA_F1 : 0   : IIV on F1
PROP_ERR : 0 : Proportional residual error

$CMT @annotated
DEPOT  : Oral depot for first-order absorption (mg)
CENT   : Central (mg)
PERIPH1 : Peripheral 1 (mg)
PERIPH2 : Peripheral 2 (mg)

$MAIN
double CLi  = CL * exp(ETA_CL) * pow(WT/75.0, 0.75) * pow(1.30, IND) * pow(0.922, INH) * pow(0.837, ALBL) * pow(0.883, ALPH) * pow(1.15, POP);
double V1i  = V1 * exp(ETA_V1) * (WT / 75.0);
double V2i  = V2 * exp(ETA_V2) * (WT / 75.0);
double V3i  = V3 * exp(ETA_V3) * (WT / 75.0);
double Q1i  = Q1 * pow(WT/75.0, 0.75);
double Q2i  = Q2 * pow(WT/75.0, 0.75);

// Simultaneous first- and zero-order absorption (Gupta et al. 2016)
// Dose split: FR_ZO fraction via zero-order into CENT, (1-FR_ZO) via first-order through DEPOT
F_DEPOT = FREL * (1.0 - FR_ZO) * exp(ETA_F1);
F_CENT  = FREL * FR_ZO * exp(ETA_F1);
D_CENT  = D1 * exp(ETA_D1);

$ODE
dxdt_DEPOT   = -KA * exp(ETA_KA) * DEPOT;
dxdt_CENT    = KA * exp(ETA_KA) * DEPOT - (CLi/V1i)*CENT - (Q1i/V1i)*CENT + (Q1i/V2i)*PERIPH1 - (Q2i/V1i)*CENT + (Q2i/V3i)*PERIPH2;
dxdt_PERIPH1 = (Q1i/V1i)*CENT - (Q1i/V2i)*PERIPH1;
dxdt_PERIPH2 = (Q2i/V1i)*CENT - (Q2i/V3i)*PERIPH2;

$TABLE
double CP = CENT / V1i * 1000.0 * (1.0 + PROP_ERR);

$CAPTURE CP
'

mod <- mcode("lenvatinib", model_code)

# IIV omega^2 values from Gupta et al. 2016 Table 2
# %CV -> omega^2: omega^2 = log(1 + (CV/100)^2)
omega2_CL <- log(1 + 0.255^2)
omega2_V1 <- log(1 + 0.228^2)
omega2_V2 <- log(1 + 0.390^2)
omega2_V3 <- log(1 + 0.303^2)
omega2_KA <- log(1 + 0.548^2)
omega2_D1 <- log(1 + 0.767^2)
omega2_F1 <- log(1 + 0.302^2)

sigma_prop_patient <- 0.333
sigma_prop_clinpharm <- 0.169

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; }
")

ui <- page_sidebar(
  title = "Lenvatinib (Lenvima\u00AE) Population PK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings", width = 340,
    h6("\U0001F48A Dosing"),
    radioButtons("dose_mode", "Dosing Mode",
                 choices = c("Standard Regimen" = "standard", "Custom Schedule" = "custom"),
                 selected = "standard", inline = TRUE),
    conditionalPanel(
      condition = "input.dose_mode == 'standard'",
      sliderInput("dose", "Dose (mg)", min = 4, max = 24, value = 24, step = 2),
      sliderInput("interval", "Dosing Interval (h)", min = 12, max = 48, value = 24, step = 12),
      numericInput("n_days", "Duration (days)", value = 14, min = 1, max = 28)
    ),
    conditionalPanel(
      condition = "input.dose_mode == 'custom'",
      textAreaInput("custom_doses", "Dosing Schedule (Day, Dose mg — one per line)",
                    value = "0, 24\n7, 20\n14, 14",
                    rows = 6, placeholder = "Day, Dose (mg)\n0, 24\n7, 20"),
      div(style = "font-size: 11px; color: #6c757d; margin-top: -8px; margin-bottom: 8px;",
          "Enter one dose event per line: Day, Dose (mg).",
          br(), "Each line starts a new regimen from that day.",
          br(), "Example: dose reduce from 24→20 mg on day 7."),
      numericInput("custom_interval", "Dosing Interval (h)", value = 24, min = 8, max = 72, step = 4),
      numericInput("custom_end_day", "Simulation End (days)", value = 28, min = 1, max = 56)
    ),
    hr(), h6("\U0001F9EC Patient Characteristics"),
    sliderInput("wt", "Body Weight (kg)", min = 40, max = 130, value = 75, step = 5),
    radioButtons("pop", "Population", choices = c("Cancer Patient" = 0, "Healthy Volunteer" = 1), selected = 0, inline = TRUE),
    hr(), h6("\u2699\uFE0F Covariates"),
    checkboxInput("cyp_ind", "CYP3A4 Inducer (e.g., rifampin)", value = FALSE),
    checkboxInput("cyp_inh", "CYP3A4 Inhibitor (e.g., ketoconazole)", value = FALSE),
    checkboxInput("alb_low", "Low Albumin (<30 g/L)", value = FALSE),
    checkboxInput("alp_high", "Elevated ALP (> ULN)", value = FALSE),
    hr(), checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE),
    hr(), h6("\U0001F4CA Population Simulation"),
    sliderInput("pi_pct", "Prediction Interval (%)", min = 50, max = 99, value = 90, step = 5),
    numericInput("n_subj", "Virtual Subjects (N)", value = 500, min = 50, max = 2000, step = 50),
    checkboxInput("include_ruv", "Include Residual Variability", value = TRUE),
    radioButtons("pop_ruv", "RUV Type",
                 choices = c("Patient (33.3% CV)" = "patient",
                             "Clinical Pharm (16.9% CV)" = "clinpharm"),
                 selected = "patient", inline = FALSE)
  ),
  navset_card_tab(
    title = "Lenvatinib PK Simulator",
    full_screen = TRUE,
    nav_panel("Individual Simulation",
      layout_column_wrap(
        width = 1/4,
        fill = FALSE,
        div(class = "metric-card metric-success", div(class = "metric-value", textOutput("cmax")), div(class = "metric-label", "Cmax,ss (ng/mL)")),
        div(class = "metric-card metric-warning", div(class = "metric-value", textOutput("ctrough")), div(class = "metric-label", "Ctrough,ss (ng/mL)")),
        div(class = "metric-card metric-primary", div(class = "metric-value", textOutput("auc")), div(class = "metric-label", "AUC\u03C4,ss (ng\u00B7h/mL)")),
        div(class = "metric-card metric-info", div(class = "metric-value", textOutput("thalf")), div(class = "metric-label", "Effective t\u00BD (h)"))
      ),
      plotOutput("pkPlot", height = "500px"),
      div(style = "font-size: 12px; color: #6c757d; padding: 8px;",
        "Typical patient simulation (population mean parameters). No inter-individual variability."
      )
    ),
    nav_panel("Population Simulation",
      layout_column_wrap(
        width = 1/4,
        fill = FALSE,
        div(class = "metric-card metric-success", div(class = "metric-value", textOutput("pop_cmax")), div(class = "metric-label", "Median Cmax,ss (ng/mL)")),
        div(class = "metric-card metric-warning", div(class = "metric-value", textOutput("pop_ctrough")), div(class = "metric-label", "Median Ctrough,ss (ng/mL)")),
        div(class = "metric-card metric-primary", div(class = "metric-value", textOutput("pop_cmax_range")), div(class = "metric-label", "Cmax,ss PI Range")),
        div(class = "metric-card metric-info", div(class = "metric-value", textOutput("pop_ctrough_range")), div(class = "metric-label", "Ctrough,ss PI Range"))
      ),
      plotOutput("popPlot", height = "500px"),
      div(style = "font-size: 12px; color: #6c757d; padding: 8px;",
        "Monte Carlo simulation with inter-individual variability (IIV) on CL, V1, V2, V3, Ka, and F1. ",
        "Shaded region = prediction interval. Solid line = typical patient. Dashed line = population median."
      )
    ),
    nav_panel("Dosing Schedule",
      conditionalPanel(
        condition = "input.dose_mode == 'standard'",
        div(style = "padding: 20px; color: #6c757d;",
          h5("Standard Regimen Active"),
          p("Switch to 'Custom Schedule' in the sidebar to define variable dosing."),
          tableOutput("standard_sched_table")
        )
      ),
      conditionalPanel(
        condition = "input.dose_mode == 'custom'",
        div(style = "padding: 20px;",
          h5("Custom Dosing Schedule"),
          tableOutput("custom_sched_table"),
          plotOutput("dose_timeline", height = "200px")
        )
      )
    ),
    nav_panel("Model Information", markdown("
## Lenvatinib (Lenvima\u00AE) \u2014 3-Compartment Population PK Model

**Source:** Gupta A et al. *Br J Clin Pharmacol* (2016) 81:1124\u20131133

### Model Structure
- **Type:** Three-compartment with simultaneous first- and zero-order absorption and linear elimination
- **Software:** NONMEM 7.2.0 (original); mrgsolve (this simulator)
- **Data:** Pooled from 15 clinical studies (n=779 subjects), 10,265 plasma concentrations

### Population PK Parameters (Typical Values)
| Parameter | Value | Description |
|-----------|-------|-------------|
| CL/F | 6.56 L/h | Apparent clearance (%CV 25.5) |
| V1/F | 49.3 L | Central volume |
| V2/F | 30.7 L | Peripheral volume 1 |
| V3/F | 37.1 L | Peripheral volume 2 |
| Q1/F | 3.52 L/h | Intercompartmental CL (V1\u2194V2) |
| Q2/F | 0.769 L/h | Intercompartmental CL (V1\u2194V3) |
| Ka | 1.02 h\u207B\u00B9 | First-order absorption rate constant |
| D1 | 1.22 h | Duration of zero-order absorption |
| F1 | 0.896 | Relative bioavailability (capsule vs tablet) |
| t\u00BD | ~28 h | Terminal half-life |

### Inter-Individual Variability (IIV)
| Parameter | %CV | Omega\u00B2 |
|-----------|-----|---------|
| CL/F | 25.5% | 0.063 |
| V1/F | 22.8% | 0.051 |
| V2/F | 39.0% | 0.141 |
| V3/F | 30.3% | 0.088 |
| Ka | 54.8% | 0.261 |
| D1 | 76.7% | 0.453 |
| F1 | 30.2% | 0.087 |

### Residual Variability
| Population | Proportional %CV |
|------------|-----------------|
| Patient Studies | 33.3% |
| Clinical Pharmacology | 16.9% |

### Covariate Effects on CL/F
- **Body weight:** Allometric (power 0.75), reference 75 kg
- **CYP3A4 inducers:** +30% CL/F (e.g., rifampin)
- **CYP3A4 inhibitors:** \u22127.8% CL/F (e.g., ketoconazole)
- **Low albumin (<30 g/L):** \u221216.3% CL/F
- **Elevated ALP (> ULN):** \u221211.7% CL/F
- **Population:** Healthy subjects +15% CL/F vs cancer patients

### Clinical Context
- **Approved dose (RR-DTC):** 24 mg once daily (capsules)
- **Dose reductions:** 20 mg \u2192 14 mg \u2192 10 mg (for toxicity management)
- **Steady-state AUC:** ~3710 ng\u00B7h/mL at 24 mg QD (SELECT trial)
- **Linear PK** confirmed across 3.2\u201332 mg dose range
    ")),
    nav_panel("References", div(class = "ref-box",
      tags$h5("\U0001F4DA Key References"),
      tags$ol(
        tags$li("Gupta A, Jarzab B, Capdevila J, Shumaker R, Hussein Z. Population pharmacokinetic analysis of lenvatinib in healthy subjects and patients with cancer. Br J Clin Pharmacol. 2016;81(6):1124\u20131133."),
        tags$li("Schlumberger M, Tahara M, Wirth LJ, et al. Lenvatinib versus placebo in radioiodine-refractory thyroid cancer (SELECT). N Engl J Med. 2015;372(7):621\u2013630."),
        tags$li("Shumaker R, et al. Effect of lenvatinib (E7080) on the QTc interval: results from a thorough QT study in healthy volunteers. Cancer Chemother Pharmacol. 2014;73:1109\u20131117."),
        tags$li(tags$a(href = "https://www.accessdata.fda.gov/drugsatfda_docs/label/2021/206947s021lbl.pdf", target = "_blank", "Lenvima\u00AE (lenvatinib) FDA Prescribing Information"))
      ),
      tags$h5("\U0001F48A Therapeutic Context"),
      tags$ul(
        tags$li(tags$strong("Drug:"), " Lenvatinib mesylate (Lenvima\u00AE, E7080)"),
        tags$li(tags$strong("Class:"), " Oral multikinase inhibitor (VEGFR1-3, FGFR1-4, PDGFR\u03B1, RET, KIT)"),
        tags$li(tags$strong("Indications:"), " RR-DTC (24 mg QD), RCC (18 mg + everolimus), HCC (8/12 mg), endometrial cancer (20 mg + pembrolizumab)"),
        tags$li(tags$strong("Formulation:"), " 4 mg and 10 mg capsules"),
        tags$li(tags$strong("Metabolism:"), " CYP3A4 and aldehyde oxidase; biliary excretion (~64% feces, ~25% urine)"),
        tags$li(tags$strong("Protein binding:"), " 97.9\u201398.6% (primarily albumin)")
      )
    ))
  ),
  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 <- function(input, output, session) {

  # ── Helper: parse custom dosing text into event objects ──
  parse_custom_doses <- function(text, interval_h, end_day) {
    lines <- strsplit(trimws(text), "\n")[[1]]
    lines <- lines[nchar(trimws(lines)) > 0]
    events <- list()
    for (ln in lines) {
      parts <- as.numeric(strsplit(trimws(ln), "[,\\s]+")[[1]])
      if (length(parts) >= 2 && !any(is.na(parts))) {
        day_start <- parts[1]
        dose_mg <- parts[2]
        events <- c(events, list(data.frame(day = day_start, dose = dose_mg)))
      }
    }
    if (length(events) == 0) return(NULL)
    sched <- do.call(rbind, events)
    sched <- sched[order(sched$day), ]

    # Build mrgsolve events: each row defines a regimen segment
    ev_list <- list()
    for (i in seq_len(nrow(sched))) {
      start_h <- sched$day[i] * 24
      if (i < nrow(sched)) {
        next_start_h <- sched$day[i + 1] * 24
      } else {
        next_start_h <- end_day * 24
      }
      n_doses_seg <- max(0, floor((next_start_h - start_h) / interval_h))
      addl_seg <- max(n_doses_seg - 1, 0)
      # Depot (first-order) and central (zero-order) events
      ev_dep <- ev(amt = sched$dose[i], cmt = 1, time = start_h, ii = interval_h, addl = addl_seg)
      ev_cen <- ev(amt = sched$dose[i], cmt = 2, time = start_h, ii = interval_h, addl = addl_seg)
      ev_list <- c(ev_list, list(ev_dep), list(ev_cen))
    }
    do.call(c, ev_list)
  }

  # ── Reactive: simulation end time (hours) ──
  sim_end_h <- reactive({
    if (input$dose_mode == "custom") input$custom_end_day * 24 else input$n_days * 24
  })

  # ── Reactive: dosing interval for metrics ──
  last_interval_h <- reactive({
    if (input$dose_mode == "custom") input$custom_interval else input$interval
  })

  # ── Typical (single subject) simulation ──
  sim_data <- reactive({
    if (input$dose_mode == "custom") {
      ev_both <- parse_custom_doses(input$custom_doses, input$custom_interval, input$custom_end_day)
      if (is.null(ev_both)) return(data.frame(time = 0, CP = 0, time_d = 0))
      end_h <- input$custom_end_day * 24
    } else {
      n_doses <- floor(input$n_days * 24 / input$interval)
      ev_depot <- ev(amt = input$dose, cmt = 1, ii = input$interval, addl = max(n_doses - 1, 0))
      ev_cent  <- ev(amt = input$dose, cmt = 2, ii = input$interval, addl = max(n_doses - 1, 0))
      ev_both  <- ev_depot + ev_cent
      end_h <- input$n_days * 24
    }
    mod %>%
      param(
        WT = input$wt, POP = as.numeric(input$pop),
        IND = as.numeric(input$cyp_ind), INH = as.numeric(input$cyp_inh),
        ALBL = as.numeric(input$alb_low), ALPH = as.numeric(input$alp_high),
        ETA_CL = 0, ETA_V1 = 0, ETA_V2 = 0, ETA_V3 = 0,
        ETA_KA = 0, ETA_D1 = 0, ETA_F1 = 0, PROP_ERR = 0
      ) %>%
      ev(ev_both) %>%
      mrgsim(end = end_h, delta = 0.25) %>%
      as.data.frame() %>%
      mutate(time_d = time / 24)
  })

  # ── Population Monte Carlo simulation ──
  pop_data <- reactive({
    n_subj <- input$n_subj

    set.seed(12345)
    idata <- data.frame(
      ID = 1:n_subj,
      ETA_CL = rnorm(n_subj, 0, sqrt(omega2_CL)),
      ETA_V1 = rnorm(n_subj, 0, sqrt(omega2_V1)),
      ETA_V2 = rnorm(n_subj, 0, sqrt(omega2_V2)),
      ETA_V3 = rnorm(n_subj, 0, sqrt(omega2_V3)),
      ETA_KA = rnorm(n_subj, 0, sqrt(omega2_KA)),
      ETA_D1 = rnorm(n_subj, 0, sqrt(omega2_D1)),
      ETA_F1 = rnorm(n_subj, 0, sqrt(omega2_F1)),
      WT = input$wt, POP = as.numeric(input$pop),
      IND = as.numeric(input$cyp_ind), INH = as.numeric(input$cyp_inh),
      ALBL = as.numeric(input$alb_low), ALPH = as.numeric(input$alp_high),
      PROP_ERR = 0
    )

    if (input$dose_mode == "custom") {
      ev_both <- parse_custom_doses(input$custom_doses, input$custom_interval, input$custom_end_day)
      if (is.null(ev_both)) return(data.frame(time = 0, CP = 0, ID = 1))
      end_h <- input$custom_end_day * 24
    } else {
      n_doses <- floor(input$n_days * 24 / input$interval)
      ev_depot <- ev(amt = input$dose, cmt = 1, ii = input$interval, addl = max(n_doses - 1, 0))
      ev_cent  <- ev(amt = input$dose, cmt = 2, ii = input$interval, addl = max(n_doses - 1, 0))
      ev_both  <- ev_depot + ev_cent
      end_h <- input$n_days * 24
    }

    out <- mod %>%
      ev(ev_both) %>%
      idata_set(idata) %>%
      mrgsim(end = end_h, delta = 0.5) %>%
      as.data.frame() %>%
      filter(is.finite(CP))

    # Add RUV if requested
    sigma_prop <- if (input$pop_ruv == "patient") sigma_prop_patient else sigma_prop_clinpharm
    if (input$include_ruv) {
      set.seed(54321)
      out$CP <- out$CP * (1 + rnorm(nrow(out), 0, sigma_prop))
      out$CP <- pmax(out$CP, 0)
    }

    out
  })

  # Summarised PI bands
  pi_summary <- reactive({
    out <- pop_data()
    ci_pct <- input$pi_pct / 100
    lower_q <- (1 - ci_pct) / 2
    upper_q <- 1 - lower_q

    out %>%
      group_by(time) %>%
      summarise(
        median = median(CP, na.rm = TRUE),
        lower = quantile(CP, lower_q, na.rm = TRUE),
        upper = quantile(CP, upper_q, na.rm = TRUE),
        .groups = "drop"
      ) %>%
      filter(is.finite(median)) %>%
      mutate(time_d = time / 24)
  })

  # ── Individual tab metrics ──
  last_interval_data <- reactive({
    d <- sim_data()
    last_start <- max(0, sim_end_h() - last_interval_h())
    d %>% filter(time >= last_start)
  })

  output$cmax <- renderText(sprintf("%.1f", max(last_interval_data()$CP, na.rm = TRUE)))
  output$ctrough <- renderText({
    vals <- last_interval_data()$CP
    sprintf("%.1f", min(vals[vals > 0], na.rm = TRUE))
  })
  output$auc <- renderText({
    d <- sim_data()
    last_start <- max(0, sim_end_h() - last_interval_h())
    last_end <- last_start + last_interval_h()
    last <- d %>% filter(time >= last_start, time <= last_end)
    if (nrow(last) < 2) return("--")
    auc_val <- sum(diff(last$time) * (head(last$CP, -1) + tail(last$CP, -1)) / 2)
    sprintf("%.0f", auc_val)
  })
  output$thalf <- renderText({
    last <- last_interval_data()
    cmax_val <- max(last$CP, na.rm = TRUE)
    tmax_idx <- which.max(last$CP)
    desc <- last[tmax_idx:nrow(last), ]
    cross <- desc %>% filter(CP <= cmax_val / 2)
    if (nrow(cross) > 0 && tmax_idx > 0) {
      sprintf("%.1f", cross$time[1] - last$time[tmax_idx])
    } else { "~28" }
  })

  # ── Population tab metrics ──
  pop_ss_metrics <- reactive({
    out <- pop_data()
    last_start <- max(0, sim_end_h() - last_interval_h())
    last_end <- last_start + last_interval_h()
    ss <- out %>% filter(time >= last_start, time <= last_end)

    ci_pct <- input$pi_pct / 100
    lower_q <- (1 - ci_pct) / 2
    upper_q <- 1 - lower_q

    per_subj <- ss %>%
      group_by(ID) %>%
      summarise(cmax = max(CP, na.rm = TRUE), ctrough = min(CP[CP > 0], na.rm = TRUE), .groups = "drop")

    list(
      med_cmax = median(per_subj$cmax, na.rm = TRUE),
      med_ctrough = median(per_subj$ctrough, na.rm = TRUE),
      cmax_lo = quantile(per_subj$cmax, lower_q, na.rm = TRUE),
      cmax_hi = quantile(per_subj$cmax, upper_q, na.rm = TRUE),
      ctrough_lo = quantile(per_subj$ctrough, lower_q, na.rm = TRUE),
      ctrough_hi = quantile(per_subj$ctrough, upper_q, na.rm = TRUE)
    )
  })

  output$pop_cmax <- renderText(sprintf("%.1f", pop_ss_metrics()$med_cmax))
  output$pop_ctrough <- renderText(sprintf("%.1f", pop_ss_metrics()$med_ctrough))
  output$pop_cmax_range <- renderText({
    m <- pop_ss_metrics()
    sprintf("%.0f - %.0f", m$cmax_lo, m$cmax_hi)
  })
  output$pop_ctrough_range <- renderText({
    m <- pop_ss_metrics()
    sprintf("%.0f - %.0f", m$ctrough_lo, m$ctrough_hi)
  })

  # ── Individual PK Plot ──
  output$pkPlot <- renderPlot({
    d <- sim_data()
    if (input$dose_mode == "custom") {
      title_str <- "Lenvatinib Custom Dosing Schedule \u2014 Typical Patient"
    } else {
      dose_lab <- paste0(input$dose, " mg")
      int_lab <- paste0("Q", input$interval, "H")
      title_str <- paste0("Lenvatinib ", dose_lab, " ", int_lab, " \u2014 Typical Patient")
    }
    sub_parts <- c(paste0("Weight: ", input$wt, " kg"),
                   ifelse(as.numeric(input$pop) == 0, "Cancer Patient", "Healthy Volunteer"))
    if (input$cyp_ind) sub_parts <- c(sub_parts, "CYP3A4 Inducer")
    if (input$cyp_inh) sub_parts <- c(sub_parts, "CYP3A4 Inhibitor")
    if (input$alb_low) sub_parts <- c(sub_parts, "Low Albumin")
    if (input$alp_high) sub_parts <- c(sub_parts, "High ALP")

    p <- ggplot(d, aes(x = time_d, y = CP)) +
      geom_line(color = "#8b5cf6", linewidth = 1) +
      labs(x = "Time (days)", y = "Plasma Concentration (ng/mL)",
           title = title_str, subtitle = paste(sub_parts, collapse = " | ")) +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(face = "bold", size = 15),
            plot.subtitle = element_text(color = "#6c757d", size = 11),
            panel.grid.minor = element_blank())

    if (input$log_scale) p <- p + scale_y_log10()
    p
  })

  # ── Population PK Plot ──
  output$popPlot <- renderPlot({
    pi <- pi_summary()
    d <- sim_data()  # overlay typical patient
    pi_lab <- paste0(input$pi_pct, "% Prediction Interval")
    ruv_lab <- if (input$include_ruv) "IIV + RUV" else "IIV only"
    if (input$dose_mode == "custom") {
      title_str <- "Lenvatinib Custom Schedule \u2014 Population Simulation"
    } else {
      dose_lab <- paste0(input$dose, " mg")
      int_lab <- paste0("Q", input$interval, "H")
      title_str <- paste0("Lenvatinib ", dose_lab, " ", int_lab, " \u2014 Population Simulation")
    }
    sub_str <- paste0("N=", input$n_subj, " virtual subjects | ", pi_lab, " (", ruv_lab, ")")

    p <- ggplot() +
      geom_ribbon(data = pi, aes(x = time_d, ymin = lower, ymax = upper),
                  fill = "#8b5cf6", alpha = 0.2) +
      geom_line(data = pi, aes(x = time_d, y = median, linetype = "Population Median"),
                color = "#8b5cf6", linewidth = 0.7) +
      geom_line(data = d, aes(x = time_d, y = CP, linetype = "Typical Patient"),
                color = "#2c3e50", linewidth = 0.9) +
      scale_linetype_manual(values = c("Typical Patient" = "solid", "Population Median" = "dashed")) +
      labs(x = "Time (days)", y = "Plasma Concentration (ng/mL)",
           title = title_str, subtitle = sub_str, linetype = "") +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(face = "bold", size = 15),
            plot.subtitle = element_text(color = "#6c757d", size = 11),
            panel.grid.minor = element_blank(),
            legend.position = "bottom")

    if (input$log_scale) p <- p + scale_y_log10()
    p
  })

  # ── Dosing Schedule tab outputs ──
  output$standard_sched_table <- renderTable({
    data.frame(
      Day = 0,
      `Dose (mg)` = input$dose,
      Interval = paste0("Q", input$interval, "H"),
      Duration = paste0(input$n_days, " days"),
      check.names = FALSE
    )
  })

  output$custom_sched_table <- renderTable({
    lines <- strsplit(trimws(input$custom_doses), "\n")[[1]]
    lines <- lines[nchar(trimws(lines)) > 0]
    rows <- lapply(lines, function(ln) {
      parts <- as.numeric(strsplit(trimws(ln), "[,\\s]+")[[1]])
      if (length(parts) >= 2 && !any(is.na(parts))) {
        data.frame(Day = parts[1], `Dose (mg)` = parts[2], check.names = FALSE)
      }
    })
    rows <- rows[!sapply(rows, is.null)]
    if (length(rows) > 0) {
      df <- do.call(rbind, rows)
      df$Interval <- paste0("Q", input$custom_interval, "H")
      df
    } else {
      data.frame(Day = numeric(), `Dose (mg)` = numeric(), Interval = character(), check.names = FALSE)
    }
  })

  output$dose_timeline <- renderPlot({
    lines <- strsplit(trimws(input$custom_doses), "\n")[[1]]
    lines <- lines[nchar(trimws(lines)) > 0]
    rows <- lapply(lines, function(ln) {
      parts <- as.numeric(strsplit(trimws(ln), "[,\\s]+")[[1]])
      if (length(parts) >= 2 && !any(is.na(parts))) data.frame(day = parts[1], dose = parts[2])
    })
    rows <- rows[!sapply(rows, is.null)]
    if (length(rows) == 0) return(NULL)
    sched <- do.call(rbind, rows)
    # Add end point for step display
    sched$day_end <- c(sched$day[-1], input$custom_end_day)

    ggplot(sched, aes(xmin = day, xmax = day_end, ymin = 0, ymax = dose)) +
      geom_rect(fill = "#8b5cf6", alpha = 0.6, color = "#8b5cf6") +
      geom_text(aes(x = (day + day_end) / 2, y = dose + 1, label = paste0(dose, " mg")),
                size = 4, fontface = "bold") +
      scale_y_continuous(limits = c(0, max(sched$dose) * 1.2)) +
      labs(x = "Day", y = "Dose (mg)", title = "Dose Timeline") +
      theme_minimal(base_size = 13)
  })
}

shinyApp(ui = ui, server = server)