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

# Molecular weights for unit conversion (mg/L → nmol/L)
MW_OSI   <- 500.59   # Osimertinib (g/mol)
MW_AZ5104 <- 486.57  # AZ5104 metabolite (g/mol)

# mrgsolve model: 1-CMT osimertinib + sequential 1-CMT AZ5104 metabolite
# Ref: Planchard et al., Br J Clin Pharmacol 2016 (doi:10.1111/bcp.13147)
model_code <- '
$PARAM @annotated
CL    : 14.2  : Apparent clearance osimertinib (L/h)
V     : 986   : Apparent volume of distribution osimertinib (L)
Ka    : 0.24  : Absorption rate constant (h-1)
CLm   : 31.5  : Apparent clearance AZ5104 metabolite (L/h)
Vm    : 207   : Apparent volume of distribution AZ5104 (L)
fm    : 0.25  : Fraction of osimertinib converted to AZ5104 (fixed)
WT    : 62    : Body weight (kg)
ALB   : 38.4  : Serum albumin (g/L)
ETHNF : 1.0   : Ethnicity factor on AZ5104 clearance (Caucasian=1.0)

$CMT @annotated
DEPOT  : Oral depot (mg)
PARENT : Osimertinib central (mg)
METAB  : AZ5104 metabolite central (mg)

$MAIN
// Covariate-adjusted parameters
// Body weight on CL (exp 0.56), V (exp 0.65), CLm (exp 0.99); albumin on V (exp 1.33)
double CLadj  = CL  * pow(WT / 62.0, 0.56);
double Vadj   = V   * pow(WT / 62.0, 0.65) * pow(ALB / 38.4, 1.33);
double CLmadj = CLm * pow(WT / 62.0, 0.99) * ETHNF;

$ODE
dxdt_DEPOT  = -Ka * DEPOT;
dxdt_PARENT =  Ka * DEPOT - (CLadj / Vadj) * PARENT;
dxdt_METAB  =  fm * (CLadj / Vadj) * PARENT - (CLmadj / Vm) * METAB;

$TABLE
// Convert mg/L → nmol/L using MW
double CP_nM = (PARENT / Vadj)  * (1e6 / 500.59);  // osimertinib nmol/L
double CM_nM = (METAB  / Vm)    * (1e6 / 486.57);  // AZ5104 nmol/L

// c-QTc: ΔQTcF (ms) = 0.271 ms per 10 nmol/L increase (linear mixed-effects model)
double dQTcF = (CP_nM / 10.0) * 0.271;

$CAPTURE CP_nM CM_nM dQTcF
'

mod <- mcode("osimertinib_popPK", 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; 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; }
  .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; }
  .footer-bar { text-align: center; font-size: 11px; color: #adb5bd; margin-top: 12px; padding-top: 8px; border-top: 1px solid #dee2e6; }
")

ui <- page_sidebar(
  title = "Osimertinib (Tagrisso) PopPK Simulator",
  theme = app_theme,

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

    h6("Dosing"),
    sliderInput("dose",     "Dose (mg)",               min = 20,  max = 240, value = 80,   step = 20),
    sliderInput("interval", "Dosing Interval (h)",     min = 12,  max = 48,  value = 24,   step = 12),
    numericInput("n_days",  "Simulation Duration (days)", value = 21, min = 7, max = 30),

    hr(),
    h6("Patient Characteristics"),
    sliderInput("wt",  "Body Weight (kg)",    min = 30,  max = 130, value = 62,   step = 1),
    sliderInput("alb", "Serum Albumin (g/L)", min = 15,  max = 55,  value = 38.4, step = 0.5),
    selectInput("ethn", "Ethnicity",
      choices = c(
        "Caucasian"                = "1.00",
        "Chinese"                  = "1.17",
        "Japanese"                 = "1.20",
        "Asian (other)"            = "1.21",
        "Non-Asian / Non-Caucasian"= "1.10"
      ),
      selected = "1.00"
    ),

    hr(),
    checkboxInput("show_metab", "Show AZ5104 metabolite", value = TRUE),
    checkboxInput("ss_only",    "Show steady-state window only (last 48 h)", value = FALSE),
    checkboxInput("log_scale",  "Log scale (Y-axis)", value = FALSE)
  ),

  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 (nmol/L)")),
    div(class = "metric-card metric-warning",
        div(class = "metric-value", textOutput("ctrough")),
        div(class = "metric-label", "Ctrough / Cmin (nmol/L)")),
    div(class = "metric-card metric-primary",
        div(class = "metric-value", textOutput("auc_ss")),
        div(class = "metric-label", "AUCss (nmol/L·h)")),
    div(class = "metric-card metric-info",
        div(class = "metric-value", textOutput("thalf")),
        div(class = "metric-label", "t½ osimertinib (h)"))
  ),

  navset_card_underline(
    title = "Osimertinib Population PK",
    full_screen = TRUE,

    nav_panel("PK Profile",
      card(
        full_screen = TRUE, height = "500px",
        plotOutput("pkPlot", height = "460px")
      ),
      div(class = "footer-bar",
          "PKPDBuilder · Powered by mrgsolve · Model: Planchard et al., Br J Clin Pharmacol 2016")
    ),

    nav_panel("Exposure-Response / QTcF",
      layout_columns(
        col_widths = c(6, 6),
        card(
          card_header("Predicted ΔQTcF at Steady State"),
          card_body(
            layout_column_wrap(
              width = 1/2, fill = FALSE,
              div(class = "metric-card metric-danger",
                  div(class = "metric-value", textOutput("qtcf_mean")),
                  div(class = "metric-label", "Mean ΔQTcF (ms)")),
              div(class = "metric-card metric-warning",
                  div(class = "metric-value", textOutput("qtcf_upper")),
                  div(class = "metric-label", "Upper 90% CI ΔQTcF (ms)"))
            ),
            div(class = "ref-box",
                tags$b("c-QTc Linear Model:"), br(),
                "ΔQTcF (ms) = 0.271 ms × (C_osimertinib / 10 nmol/L)",
                br(), "90% CI of slope: [0.241 – 0.301 ms per 10 nmol/L]",
                br(), "At 80 mg QD: predicted mean (upper 90% CI) = 14.2 (15.8) ms"
            )
          )
        ),
        card(
          card_header("Rash & Diarrhoea Risk vs. Exposure"),
          card_body(
            tableOutput("er_table"),
            div(class = "ref-box",
                "Probability modelled as logistic regression on log(AUCss).",
                br(), "Reference values are for typical NSCLC patient (WT 62 kg, Caucasian)."
            )
          )
        )
      )
    ),

    nav_panel("Model Information",
      card(
        card_body(
          markdown("
## Osimertinib (Tagrisso, AZD9291) — Population PK Model

### Drug
Osimertinib is a 3rd-generation, irreversible EGFR tyrosine kinase inhibitor (TKI) approved for:
- **1st-line NSCLC** with EGFR exon 19 deletions or L858R mutations (80 mg QD)
- **T790M-positive NSCLC** after prior EGFR TKI therapy (80 mg QD)

### Model Structure
- **Parent (osimertinib):** First-order oral absorption → 1-compartment disposition
- **Metabolite (AZ5104):** Sequential 1-compartment; formed from osimertinib via CYP3A4/5  
  - Fraction metabolised to AZ5104: **fm = 0.25** (fixed for mass balance)
  - AZ5104 circulates at ~10% of osimertinib exposure at steady state
- **Residual error:** Combined proportional + additive

### Population PK Parameters (typical NSCLC patient, 62 kg)

| Parameter | Estimate | %RSE | IIV (%) |
|-----------|----------|------|---------|
| CL/F (osimertinib) | 14.2 L/h | 1.8 | 45.6 |
| V/F (osimertinib) | 986 L | 2.8 | 51.8 |
| Ka | 0.24 h⁻¹ | 4.9 | 89.4 |
| CL/F (AZ5104) | 31.5 L/h | 2.7 | 52.3 |
| V/F (AZ5104) | 207 L | 4.4 | 51.8 |

### Covariate Effects

| Covariate | Parameter | Effect (exponent/factor) |
|-----------|-----------|--------------------------|
| Body weight | CL/F (parent) | (WT/62)^0.56 |
| Body weight | V/F (parent) | (WT/62)^0.65 |
| Serum albumin | V/F (parent) | (ALB/38.4)^1.33 |
| Body weight | CL/F (AZ5104) | (WT/62)^0.99 |
| Ethnicity (Chinese) | CL/F (AZ5104) | ×1.17 → ≈15% ↓ AUCss |
| Ethnicity (Japanese) | CL/F (AZ5104) | ×1.20 → ≈17% ↓ AUCss |
| Ethnicity (Asian other) | CL/F (AZ5104) | ×1.21 → ≈17% ↓ AUCss |
| Ethnicity (non-Asian/non-Caucasian) | CL/F (AZ5104) | ×1.10 → ≈9% ↓ AUCss |

*None of the covariates identified required dose adjustment (AUC change <2-fold across 5th–95th percentile ranges).*

### Typical Steady-State Exposures at 80 mg QD (NSCLC population)

| Metric | Osimertinib | AZ5104 |
|--------|-------------|--------|
| AUCss | 11,258 nmol/L·h | 1,271 nmol/L·h |
| Cssmax | 501 nmol/L | 56 nmol/L |
| Cssmin | 417 nmol/L | 52 nmol/L |

### Exposure–Response

- **Efficacy:** No relationship between osimertinib/AZ5104 AUCss and objective response rate, duration of response, or tumour shrinkage over the 20–240 mg dose range.
- **Safety (Rash/Diarrhoea):** Probability increases linearly with log(AUCss).
- **QTcF:** Linear relationship; mean ΔQTcF = 14.2 ms (upper 90% CI: 15.8 ms) at Cssmax for 80 mg QD. Slope = 0.271 ms per 10 nmol/L.

### Population
n = 780 subjects (748 NSCLC patients + 32 healthy volunteers); 21,930 plasma concentrations.  
Dose range: 20–240 mg. Median weight: 62 kg. Studies: AURA (phase I/II), AURA2 (phase II), Study 5 (HV).
          ")
        )
      )
    ),

    nav_panel("References",
      card(
        card_body(
          div(class = "ref-box",
            tags$b("Primary Reference:"), br(),
            "Planchard D, et al. Population pharmacokinetics and exposure-response of osimertinib in patients with non-small cell lung cancer. ",
            tags$i("Br J Clin Pharmacol."), " 2016;83(6):1216–1226. ",
            tags$a("doi:10.1111/bcp.13147", href="https://doi.org/10.1111/bcp.13147", target="_blank")
          ),
          br(),
          div(class = "ref-box",
            tags$b("Drug Information:"), br(),
            "Tagrisso (osimertinib) Prescribing Information. AstraZeneca Pharmaceuticals LP. ",
            tags$a("FDA Label", href="https://www.accessdata.fda.gov/drugsatfda_docs/label/2017/208065s006lbl.pdf", target="_blank")
          ),
          br(),
          div(class = "ref-box",
            tags$b("Modelling Software:"), br(),
            "mrgsolve: Elmokadem A, Riggs MM, Baron KT. Pharmacometric Modelling with mrgsolve. ", tags$i("CPT Pharmacometrics Syst Pharmacol."), " 2019.",
            br(), br(),
            tags$b("Built with PKPDBuilder — "), tags$a("pkpdbuilder.com", href="https://pkpdbuilder.com", target="_blank")
          )
        )
      )
    )
  )
)

server <- function(input, output, session) {

  sim_data <- reactive({
    shiny::req(input$dose, input$wt, input$alb, input$ethn, input$interval, input$n_days)

    dose_val     <- input$dose
    tau          <- input$interval
    n_doses      <- ceiling(input$n_days * 24 / tau)
    ethnf_val    <- as.numeric(input$ethn)

    # Build event table: n_doses administrations
    ev <- ev(amt = dose_val, ii = tau, addl = n_doses - 1, cmt = 1)

    # Simulation time grid: every 0.5 h
    total_h <- n_doses * tau
    tgrid_vec <- seq(0, total_h, by = 0.5)

    out <- mod |>
      param(WT = input$wt, ALB = input$alb, ETHNF = ethnf_val) |>
      mrgsim(events = ev, tgrid = tgrid_vec, output = "df") |>
      dplyr::filter(time > 0)

    out
  })

  # Derived steady-state metrics (use last full dosing interval)
  ss_metrics <- reactive({
    df <- sim_data()
    shiny::req(nrow(df) > 0)

    tau     <- input$interval
    n_doses <- ceiling(input$n_days * 24 / tau)
    total_h <- n_doses * tau

    # Last complete dosing interval
    t_start <- total_h - tau
    t_end   <- total_h

    ss_df <- dplyr::filter(df, time >= t_start, time <= t_end)

    if (nrow(ss_df) < 2) return(list(cmax=NA, ctrough=NA, auc=NA, thalf=48))

    cmax   <- max(ss_df$CP_nM, na.rm = TRUE)
    ctrough <- min(ss_df$CP_nM, na.rm = TRUE)

    # Trapezoidal AUC over last dosing interval
    auc_val <- sum(diff(ss_df$time) * (head(ss_df$CP_nM, -1) + tail(ss_df$CP_nM, -1))) / 2

    # Terminal t½ from population parameters (CL/V, adjusted for WT)
    wt_val   <- input$wt
    alb_val  <- input$alb
    cl_adj   <- 14.2 * (wt_val / 62.0)^0.56
    v_adj    <- 986  * (wt_val / 62.0)^0.65 * (alb_val / 38.4)^1.33
    ke       <- cl_adj / v_adj
    t12      <- log(2) / ke

    list(cmax = cmax, ctrough = ctrough, auc = auc_val, thalf = t12)
  })

  output$cmax    <- renderText({ m <- ss_metrics(); if (is.na(m$cmax))    "—" else format(round(m$cmax,    0), big.mark=",") })
  output$ctrough <- renderText({ m <- ss_metrics(); if (is.na(m$ctrough)) "—" else format(round(m$ctrough, 0), big.mark=",") })
  output$auc_ss  <- renderText({ m <- ss_metrics(); if (is.na(m$auc))     "—" else format(round(m$auc,     0), big.mark=",") })
  output$thalf   <- renderText({ m <- ss_metrics(); format(round(m$thalf, 1)) })

  # QTcF outputs
  output$qtcf_mean <- renderText({
    m    <- ss_metrics()
    cmax <- if (is.na(m$cmax)) 501 else m$cmax
    dqtcf <- (cmax / 10) * 0.271
    paste0(round(dqtcf, 1), " ms")
  })
  output$qtcf_upper <- renderText({
    m    <- ss_metrics()
    cmax <- if (is.na(m$cmax)) 501 else m$cmax
    dqtcf_upper <- (cmax / 10) * 0.301
    paste0(round(dqtcf_upper, 1), " ms")
  })

  # Exposure-response table
  output$er_table <- renderTable({
    m       <- ss_metrics()
    auc_sim <- if (is.na(m$auc)) 11258 else m$auc

    # Logistic ER: log-linear on AUCss (from paper, approximate from reported probabilities)
    # Intercept and slope fitted to match reported probabilities at each dose level
    # Rash:     P = plogis(-0.938 + 0.436 * log(AUCss/1000))
    # Diarrhea: P = plogis(-0.920 + 0.420 * log(AUCss/1000))
    p_rash  <- round(plogis(-0.938 + 0.436 * log(auc_sim / 1000)), 2)
    p_diarr <- round(plogis(-0.920 + 0.420 * log(auc_sim / 1000)), 2)

    data.frame(
      Endpoint     = c("Any Rash", "Any Diarrhoea"),
      `AUCss (nmol/L·h)` = c(round(auc_sim), round(auc_sim)),
      `Predicted P` = c(p_rash, p_diarr),
      `Reference (80 mg)` = c("0.46 [0.43–0.50]", "0.45 [0.42–0.49]"),
      check.names = FALSE
    )
  }, striped = TRUE, hover = TRUE, bordered = TRUE)

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

    tau     <- input$interval
    n_doses <- ceiling(input$n_days * 24 / tau)
    total_h <- n_doses * tau

    # Optionally trim to steady-state window (last 48 h)
    if (isTRUE(input$ss_only)) {
      df <- dplyr::filter(df, time >= (total_h - 48))
    }

    # Reshape for plotting
    plot_df <- data.frame(
      time  = c(df$time, df$time),
      conc  = c(df$CP_nM, df$CM_nM),
      Analyte = rep(c("Osimertinib", "AZ5104 (metabolite)"), each = nrow(df))
    )

    if (!isTRUE(input$show_metab)) {
      plot_df <- dplyr::filter(plot_df, Analyte == "Osimertinib")
    }

    # Therapeutic reference lines at 80 mg QD steady-state (published values)
    ref_cmax   <- 501   # nmol/L (paper: Cssmax osimertinib at 80 mg)
    ref_ctrough <- 417  # nmol/L (paper: Cssmin osimertinib at 80 mg)

    p <- ggplot(plot_df, aes(x = time / 24, y = conc, color = Analyte)) +
      geom_line(linewidth = 1.0) +
      scale_color_manual(
        values = c("Osimertinib" = "#8b5cf6", "AZ5104 (metabolite)" = "#f59e0b"),
        name   = NULL
      ) +
      labs(
        x     = "Time (days)",
        y     = "Plasma Concentration (nmol/L)",
        title = paste0("Osimertinib ", input$dose, " mg Q", input$interval, "H — Population PK Simulation"),
        subtitle = paste0("WT: ", input$wt, " kg | Albumin: ", input$alb, " g/L | Ethnicity: ",
                          c("1.00"="Caucasian","1.17"="Chinese","1.20"="Japanese",
                            "1.21"="Asian (other)","1.10"="Non-Asian/Non-Caucasian")[
                            as.character(as.numeric(input$ethn))])
      ) +
      theme_bw(base_size = 13) +
      theme(
        legend.position   = "top",
        plot.title        = element_text(size = 13, face = "bold"),
        plot.subtitle     = element_text(size = 11, color = "grey40"),
        panel.grid.minor  = element_blank()
      )

    # Reference lines (80 mg QD published steady-state, only show when dose=80 and osimertinib shown)
    if (input$dose == 80 && (!isTRUE(input$ss_only))) {
      p <- p +
        geom_hline(yintercept = ref_cmax,    linetype = "dashed", color = "#8b5cf6", alpha = 0.4, linewidth = 0.7) +
        geom_hline(yintercept = ref_ctrough, linetype = "dotted", color = "#8b5cf6", alpha = 0.4, linewidth = 0.7) +
        annotate("text", x = max(plot_df$time / 24) * 0.02, y = ref_cmax    + 15, label = "Published Cssmax (501)", color = "#8b5cf6", size = 3, alpha = 0.7, hjust = 0) +
        annotate("text", x = max(plot_df$time / 24) * 0.02, y = ref_ctrough - 20, label = "Published Cssmin (417)", color = "#8b5cf6", size = 3, alpha = 0.7, hjust = 0)
    }

    if (isTRUE(input$log_scale)) {
      p <- p + scale_y_log10(labels = scales::comma)
    } else {
      p <- p + scale_y_continuous(labels = scales::comma, limits = c(0, NA))
    }

    p
  })
}

shinyApp(ui = ui, server = server)