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

# ── mrgsolve Model: Minimal PBPK-TMDD with HER2 Target Shedding ──
model_code <- '
$PARAM @annotated
CL_eff    : 0.0094  : Effective clearance (L/h)
V_P       : 3.60    : Plasma/central volume (L)
V_I       : 12.0    : Interstitial/peripheral volume (L)
Q         : 0.12    : Intercompartmental clearance (L/h)
kon_m     : 2520.0  : Association rate membrane HER2 ((umol/L)^-1 h^-1)
koff_m    : 1.26    : Dissociation rate membrane HER2 (h^-1)
RM_ss     : 0.0035  : Membrane-bound HER2 steady state (umol/L)
kdeg_m    : 0.10    : HER2 degradation rate (h^-1)
kint_m    : 0.10    : Drug-RM complex internalization (h^-1)
ksh_p     : 0.0023  : Shedding rate constant (h^-1)
lam       : 0.01    : Lymphatic ECD transport rate (h^-1)
kelm_s    : 0.01    : Drug-soluble complex elimination interstitial (h^-1)
kdeg_p    : 0.20    : ECD HER2 plasma degradation (h^-1)
kelm_p    : 0.20    : Drug-ECD complex elimination plasma (h^-1)
aff_ratio : 1.0     : koff_s/koff_m ratio (1=same affinity)

$CMT @annotated
DRUG_P : Trastuzumab plasma (umol)
DRUG_I : Trastuzumab interstitial (umol)
RM     : Membrane-bound HER2 (umol/L concentration)
RS_A   : Soluble HER2 interstitial (umol amount)
RP_A   : Soluble HER2 plasma (umol amount)
DRM_A  : Drug-RM complex interstitial (umol)
DRS_A  : Drug-RS complex interstitial (umol)
DRP_A  : Drug-RP complex plasma (umol)

$MAIN
double ksyn_m = (ksh_p + kdeg_m) * RM_ss;
if (NEWIND <= 1) {
  RM_0 = RM_ss;
  double RS0_conc = ksh_p * RM_ss / lam;
  RS_A_0 = RS0_conc * V_I;
  RP_A_0 = lam * RS0_conc * V_I / kdeg_p;
  DRM_A_0 = 0.0;
  DRS_A_0 = 0.0;
  DRP_A_0 = 0.0;
}

$ODE
double Cp = DRUG_P / V_P;
double Ci = DRUG_I / V_I;
double RS = RS_A / V_I;
double RP = RP_A / V_P;
double kon_s = kon_m;
double koff_s = koff_m * aff_ratio;
double ksyn_m_ode = (ksh_p + kdeg_m) * RM_ss;
double bind_RM = kon_m * Ci * RM - koff_m * (DRM_A / V_I);
double bind_RS = kon_s * Ci * RS - koff_s * (DRS_A / V_I);
double bind_RP = kon_s * Cp * RP  - koff_s * (DRP_A / V_P);
dxdt_DRUG_P = -CL_eff * Cp - Q * Cp + Q * Ci - bind_RP * V_P;
dxdt_DRUG_I = Q * Cp - Q * Ci - bind_RM * V_I - bind_RS * V_I;
dxdt_RM = ksyn_m_ode - kdeg_m * RM - ksh_p * RM - bind_RM;
dxdt_RS_A = ksh_p * RM * V_I - lam * RS_A - bind_RS * V_I;
dxdt_RP_A = lam * RS_A - kdeg_p * RP_A - bind_RP * V_P;
dxdt_DRM_A = bind_RM * V_I - kint_m * DRM_A;
dxdt_DRS_A = bind_RS * V_I - (kelm_s + lam) * DRS_A;
dxdt_DRP_A = bind_RP * V_P + lam * DRS_A - kelm_p * DRP_A;

$TABLE
double Cp_ugml = (DRUG_P / V_P) * 146.0;
double Ci_ugml = (DRUG_I / V_I) * 146.0;
double ECD_ngml = (RP_A / V_P) * 97000.0;
double RM_occ = 1.0 - RM / RM_ss;
double RM_free = RM / RM_ss;
double RM_conc = RM;

$CAPTURE Cp_ugml Ci_ugml ECD_ngml RM_occ RM_free RM_conc
'

mod <- mcode("trastuzumab_pbpk_tmdd", model_code)

# ── Helper Functions ──
ecd_to_ksh <- function(ecd_ngml, RM_ss = 0.0035, kdeg_p = 0.2,
                       V_P = 3.6, V_I = 12, MW_ECD = 97000) {
  ecd_ngml * kdeg_p * V_P / (RM_ss * V_I * MW_ECD)
}

make_events <- function(load_mgkg, maint_mgkg, wt_kg, interval_h, n_doses) {
  load_umol  <- load_mgkg  * wt_kg / 146
  maint_umol <- maint_mgkg * wt_kg / 146
  infusion_dur <- 1.5
  ev_list <- lapply(seq_len(n_doses), function(i) {
    dose_amt <- if (i == 1) load_umol else maint_umol
    t_start  <- (i - 1) * interval_h
    ev(amt = dose_amt, cmt = 1, time = t_start, rate = dose_amt / infusion_dur)
  })
  do.call(c, ev_list)
}

# ── PKPDBuilder Brand Theme ──
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-primary .metric-value { color: #8b5cf6; }
  .metric-slate .metric-value { color: #64748b; }
  .metric-success .metric-value { color: #10b981; }
  .metric-warning .metric-value { color: #f59e0b; }
  .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 ──
ui <- page_sidebar(
  title = "Trastuzumab PBPK-TMDD Simulator",
  theme = app_theme,

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

    h6("Dosing"),
    sliderInput("load_dose", "Loading Dose (mg/kg)",
                min = 2, max = 8, value = 4, step = 0.5),
    sliderInput("maint_dose", "Maintenance Dose (mg/kg)",
                min = 1, max = 6, value = 2, step = 0.5),
    selectInput("interval", "Dosing Interval",
                choices = c("Weekly (q1w)" = 168,
                            "Every 3 weeks (q3w)" = 504),
                selected = 168),
    sliderInput("n_doses", "Number of Doses",
                min = 1, max = 12, value = 8, step = 1),

    hr(),
    h6("Patient"),
    sliderInput("wt", "Body Weight (kg)",
                min = 40, max = 120, value = 73, step = 1),

    hr(),
    h6("HER2 Shedding"),
    sliderInput("ecd_her2", "ECD HER2 (ng/mL)",
                min = 0, max = 2200, value = 14.6, step = 10),
    tags$small(style = "color:#f59e0b;font-weight:600;",
               "500\u20131,000 ng/mL = clinical high-risk range"),
    selectInput("aff_ratio", "Soluble HER2 Affinity Ratio",
                choices = c("1x (same as membrane)" = 1,
                            "10x lower" = 10,
                            "100x lower" = 100,
                            "1000x lower" = 1000),
                selected = 1),

    hr(),
    h6("Display"),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE)
  ),

  # ── Metric Cards ──
  layout_column_wrap(
    width = 1 / 4,
    fill = FALSE,
    div(class = "metric-card metric-primary",
        div(class = "metric-value", textOutput("cmax_val")),
        div(class = "metric-label", "Plasma Cmax (\u00b5g/mL)")),
    div(class = "metric-card metric-slate",
        div(class = "metric-value", textOutput("trough_val")),
        div(class = "metric-label", "Trough (\u00b5g/mL)")),
    div(class = "metric-card metric-success",
        div(class = "metric-value", textOutput("ro_val")),
        div(class = "metric-label", "HER2 Receptor Occupancy (%)")),
    div(class = "metric-card metric-warning",
        div(class = "metric-value", textOutput("ecd_val")),
        div(class = "metric-label", "ECD HER2 (ng/mL)"))
  ),

  # ── PK Plot ──
  card(
    card_header("Trastuzumab Plasma Concentration-Time Profile"),
    full_screen = TRUE,
    min_height = "500px",
    plotOutput("pk_plot", height = "500px")
  ),

  # ── Tabbed Content ──
  navset_card_underline(
    title = "Details",
    nav_panel(
      "Receptor Dynamics",
      card(
        card_header(
          "Membrane HER2 Receptor Occupancy",
          tags$small(style = "color:#6c757d;margin-left:8px;",
                     "100% = fully occupied")
        ),
        plotOutput("ro_plot", height = "450px")
      )
    ),
    nav_panel(
      "HER2 Shedding",
      card(
        card_header("Plasma ECD HER2 Over Time"),
        plotOutput("ecd_plot", height = "450px")
      )
    ),
    nav_panel(
      "Model Information",
      div(
        style = "padding:16px;",
        tags$h4("Minimal PBPK-TMDD with HER2 Target Shedding"),
        tags$p(
          tags$strong("Structure:"),
          " 8 ODEs \u2014 2 drug compartments (plasma + interstitial),",
          " membrane-bound HER2 (RM), soluble HER2 in interstitial (RS)",
          " and plasma (RP), plus 3 drug-target complexes (DRM, DRS, DRP)."
        ),
        tags$h5("Key Finding"),
        tags$p(
          "ECD HER2 levels of 500\u20131,000 ng/mL cause a sharp decline in",
          " trastuzumab trough concentrations. Patients with high circulating",
          " soluble HER2 may fail to maintain therapeutic drug levels",
          " (target trough: 10\u201320 \u00b5g/mL)."
        ),
        tags$h5("Key Parameters"),
        tags$table(
          style = "width:100%;border-collapse:collapse;margin-top:8px;",
          tags$thead(
            tags$tr(
              style = "border-bottom:2px solid #dee2e6;",
              tags$th(style = "text-align:left;padding:6px;", "Parameter"),
              tags$th(style = "text-align:left;padding:6px;", "Value"),
              tags$th(style = "text-align:left;padding:6px;", "Description")
            )
          ),
          tags$tbody(
            tags$tr(
              tags$td(style = "padding:6px;", "kon"),
              tags$td(style = "padding:6px;", "2,520 (\u00b5mol/L)\u207b\u00b9 h\u207b\u00b9"),
              tags$td(style = "padding:6px;", "Association rate constant")
            ),
            tags$tr(
              style = "background:#f8f9fa;",
              tags$td(style = "padding:6px;", "KD"),
              tags$td(style = "padding:6px;", "0.5 nmol/L"),
              tags$td(style = "padding:6px;", "Dissociation constant (koff/kon)")
            ),
            tags$tr(
              tags$td(style = "padding:6px;", "ksh_p"),
              tags$td(style = "padding:6px;", "0\u20130.4 h\u207b\u00b9"),
              tags$td(style = "padding:6px;",
                      "Shedding rate range (0.0023 = 14.6 ng/mL ECD)")
            ),
            tags$tr(
              style = "background:#f8f9fa;",
              tags$td(style = "padding:6px;", "RM_ss"),
              tags$td(style = "padding:6px;", "0.0035 \u00b5mol/L"),
              tags$td(style = "padding:6px;",
                      "Steady-state membrane HER2 concentration")
            ),
            tags$tr(
              tags$td(style = "padding:6px;", "CL_eff"),
              tags$td(style = "padding:6px;", "0.0094 L/h"),
              tags$td(style = "padding:6px;", "Effective (linear) clearance")
            )
          )
        ),
        tags$h5(style = "margin-top:16px;", "Model Diagram"),
        tags$pre(
          style = paste0(
            "background:#f8f9fa;padding:16px;border-radius:6px;",
            "font-size:12px;line-height:1.6;overflow-x:auto;"
          ),
"Plasma                          Interstitial (tissue)
\u250c\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2510   Q    \u250c\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2510
\u2502  DRUG_P (Ab)    \u2502 \u2500\u2500\u2500\u2500\u2500 \u2502  DRUG_I (Ab)    \u2502
\u2502  CL_eff \u2193       \u2502        \u2502                 \u2502
\u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518        \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518
      \u2502  bind_RP                  \u2502  bind_RM + bind_RS
\u250c\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2510        \u250c\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2510
\u2502  RP (ECD HER2)  \u2502        \u2502  RM (memb HER2) \u2502
\u2502  kdeg_p \u2193       \u2502  lam   \u2502  kdeg_m \u2193       \u2502
\u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2500\u2500\u2500\u2500 \u2502  ksh_p \u2193 shed   \u2502
      \u2502  kelm_p              \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518
\u250c\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2510        \u250c\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2510
\u2502  DRP (complex)  \u2502  lam   \u2502  DRS (complex)  \u2502
\u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518 \u2500\u2500\u2500\u2500 \u2502  kelm_s \u2193       \u2502
                           \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518
                           \u250c\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2510
                           \u2502  DRM (complex)  \u2502
                           \u2502  kint_m \u2193       \u2502
                           \u2514\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2518"
        )
      )
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Key References"),
        tags$ol(
          tags$li(
            "Gibiansky E, et al. (2013). Incorporating Target Shedding Into",
            " a Minimal PBPK-TMDD Model for Monoclonal Antibodies.",
            tags$em(" CPT: Pharmacometrics & Systems Pharmacology"),
            ", 2(11):1-10. DOI: 10.1038/psp.2013.73"
          ),
          tags$li(
            "Mager DE, Jusko WJ (2001). General pharmacokinetic model for",
            " drugs exhibiting target-mediated drug disposition.",
            tags$em(" J Pharmacokinet Pharmacodyn"),
            ", 28(6):507-532."
          ),
          tags$li(
            "Baselga J, et al. (1996). Phase II study of weekly intravenous",
            " recombinant humanized anti-p185HER2 monoclonal antibody in",
            " HER2/neu-overexpressing metastatic breast cancer.",
            tags$em(" J Clin Oncol"),
            ", 14(3):737-744."
          )
        ),
        tags$h5("Therapeutic Context"),
        tags$ul(
          tags$li(tags$strong("Drug:"), " Trastuzumab (Herceptin\u00ae)"),
          tags$li(tags$strong("Class:"),
                  " Monoclonal Antibody (HER2-targeted), 146 kDa"),
          tags$li(tags$strong("Indication:"),
                  " HER2+ metastatic breast cancer"),
          tags$li(tags$strong("Target trough:"),
                  " 10\u201320 \u00b5g/mL (preclinical studies)"),
          tags$li(tags$strong("Route:"),
                  " IV infusion (90-minute)")
        )
      )
    )
  ),

  # ── Footer ──
  tags$footer(
    style = paste0(
      "text-align:center;padding:12px;color:#94a3b8;font-size:0.75rem;",
      "border-top:1px solid #e2e8f0;margin-top:16px;"
    ),
    tags$a(
      href = "https://pkpdbuilder.com", target = "_blank",
      style = "color:#8b5cf6;font-weight:600;text-decoration:none;",
      "PKPDBuilder.com"
    ),
    " \u00b7 Trastuzumab PBPK-TMDD Simulator \u00b7 Gibiansky et al. (CPT:PSP 2013)"
  )
)

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

  sim_data <- reactive({
    shiny::req(input$load_dose, input$maint_dose, input$wt,
               input$interval, input$n_doses)

    interval_h <- as.numeric(input$interval)
    aff        <- as.numeric(input$aff_ratio)
    ksh_val    <- ecd_to_ksh(input$ecd_her2)

    ev_doses <- make_events(
      load_mgkg  = input$load_dose,
      maint_mgkg = input$maint_dose,
      wt_kg      = input$wt,
      interval_h = interval_h,
      n_doses    = input$n_doses
    )

    end_time <- (input$n_doses - 1) * interval_h + interval_h + 24

    out <- mod %>%
      param(ksh_p = ksh_val, aff_ratio = aff) %>%
      ev(ev_doses) %>%
      mrgsim(end = end_time, delta = 0.5) %>%
      as.data.frame() |>
      mutate(time_days = time / 24)

    out
  })

  # ── Metric Cards ──
  output$cmax_val <- renderText({
    d <- sim_data()
    shiny::req(nrow(d) > 0)
    sprintf("%.1f", max(d$Cp_ugml, na.rm = TRUE))
  })

  output$trough_val <- renderText({
    d <- sim_data()
    shiny::req(nrow(d) > 0)
    interval_h <- as.numeric(input$interval)
    last_dose_time <- (input$n_doses - 1) * interval_h
    trough_time <- last_dose_time + interval_h
    trough_row <- d |> dplyr::filter(time >= trough_time - 1)
    if (nrow(trough_row) > 0) {
      sprintf("%.1f", min(trough_row$Cp_ugml, na.rm = TRUE))
    } else {
      sprintf("%.1f", tail(d$Cp_ugml, 1))
    }
  })

  output$ro_val <- renderText({
    d <- sim_data()
    shiny::req(nrow(d) > 0)
    last_occ <- tail(d$RM_occ, 1)
    sprintf("%.1f", last_occ * 100)
  })

  output$ecd_val <- renderText({
    d <- sim_data()
    shiny::req(nrow(d) > 0)
    sprintf("%.1f", tail(d$ECD_ngml, 1))
  })

  # ── PK Plot ──
  output$pk_plot <- renderPlot({
    d <- sim_data()
    shiny::req(nrow(d) > 0)

    d_plot <- if (input$log_scale) {
      d |> dplyr::filter(Cp_ugml > 0.01)
    } else {
      d
    }

    p <- ggplot(d_plot, aes(x = time_days, y = Cp_ugml)) +
      annotate("rect",
               xmin = -Inf, xmax = Inf,
               ymin = 10, ymax = 20,
               fill = "gray80", alpha = 0.3) +
      annotate("text",
               x = max(d_plot$time_days) * 0.98, y = 15,
               label = "Target trough\n10\u201320 \u00b5g/mL",
               hjust = 1, size = 3, color = "gray50") +
      geom_line(color = "#8b5cf6", linewidth = 1) +
      labs(
        x = "Time (days)",
        y = "Plasma Concentration (\u00b5g/mL)",
        title = paste0(
          "Trastuzumab ", input$load_dose, "/",
          input$maint_dose, " mg/kg ",
          ifelse(as.numeric(input$interval) == 168, "q1w", "q3w"),
          " | BW ", input$wt, " kg",
          " | ECD HER2 ", input$ecd_her2, " ng/mL"
        )
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(size = 13, face = "bold"),
        panel.grid.minor = element_blank()
      )

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

    p
  })

  # ── Receptor Occupancy Plot ──
  output$ro_plot <- renderPlot({
    d <- sim_data()
    shiny::req(nrow(d) > 0)

    ggplot(d, aes(x = time_days, y = RM_occ * 100)) +
      geom_line(color = "#10b981", linewidth = 1) +
      labs(
        x = "Time (days)",
        y = "Receptor Occupancy (%)",
        title = "Membrane HER2 Receptor Occupancy"
      ) +
      ylim(0, 100) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(size = 13, face = "bold"),
        panel.grid.minor = element_blank()
      )
  })

  # ── ECD HER2 Plot ──
  output$ecd_plot <- renderPlot({
    d <- sim_data()
    shiny::req(nrow(d) > 0)

    baseline_ecd <- input$ecd_her2

    ggplot(d, aes(x = time_days, y = ECD_ngml)) +
      geom_hline(yintercept = baseline_ecd,
                 linetype = "dashed", color = "red", linewidth = 0.6) +
      annotate("text",
               x = max(d$time_days) * 0.98, y = baseline_ecd,
               label = paste0("Baseline: ", baseline_ecd, " ng/mL"),
               hjust = 1, vjust = -0.5, size = 3, color = "red") +
      geom_line(color = "#f59e0b", linewidth = 1) +
      labs(
        x = "Time (days)",
        y = "ECD HER2 (ng/mL)",
        title = "Plasma Soluble ECD HER2 Concentration"
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(size = 13, face = "bold"),
        panel.grid.minor = element_blank()
      )
  })
}

shinyApp(ui = ui, server = server)