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

utils::globalVariables(c(
  "time", "compound", "concentration_nM", "TAM", "NDT", "OHTAM",
  "ENDO", "ZENDO", "NOXTAM", "OHTAM4P", "AAS", "dose", "endo",
  "aas", "phenotype", "target", "value", "metric", "ENDO_NM",
  "TAM_NM", "NDT_NM", "OHTAM_NM", "day", "dose_label",
  "endo_label", "aas_label", "endo_target", "aas_target"
))

MW <- c(
  TAM = 371.52,
  NDT = 357.49,
  OHTAM = 387.52,
  ENDO = 373.49,
  ZENDO = 373.49,
  NOXTAM = 387.52,
  OHTAM4P = 387.52
)

BASELINE <- list(
  KA = 0.90,
  VTAM = 1380,
  KTAM_NDT = 5.20e-3,
  KTAM_4OHTAM = 3.72e-5,
  KTAM_4POHTAM = 6.16e-8,
  KTAM_NOXTAM = 2.48e-7,
  KNDT_ENDO_NM = 5.42e-4,
  KNDT_ZENDO = 4.08e-7,
  K4OHTAM_ENDO = 1.81e-3,
  KE_NDT = 2.46e-3,
  KE_ENDO = 7.93e-3,
  KE_ZENDO = 7.93e-3,
  KE_4OHTAM = 7.93e-3,
  KE_4POHTAM = 7.93e-3,
  KE_NOXTAM = 7.93e-3
)

K_ENDO_BY_PHENO <- c(
  "Ultra-rapid metabolizer" = 6.87e-4,
  "Normal metabolizer" = 5.42e-4,
  "Intermediate metabolizer" = 2.86e-4,
  "Poor metabolizer" = 0.88e-4
)

K_ENDO_BY_DIPLOTYPE <- c(
  "NM/NM/IM or UM/UM" = 7.24e-4,
  "NM/NM" = 6.39e-4,
  "NM/IM" = 4.36e-4,
  "NM/PM" = 3.61e-4,
  "IM/IM" = 1.93e-4,
  "IM/PM" = 1.68e-4,
  "PM/PM" = 0.94e-4
)

model_code <- "
$PARAM @annotated
KA             : 0.90    : Absorption rate constant (1/h)
VTAM           : 1380    : Fixed apparent volume for tamoxifen and metabolites (L)
KTAM_NDT       : 5.20e-3 : TAM to N-desmethyltamoxifen formation (1/h)
KTAM_4OHTAM    : 3.72e-5 : TAM to 4-hydroxytamoxifen formation (1/h)
KTAM_4POHTAM   : 6.16e-8 : TAM to 4'-hydroxytamoxifen formation (1/h)
KTAM_NOXTAM    : 2.48e-7 : TAM to tamoxifen N-oxide formation (1/h)
KNDT_ENDO      : 5.42e-4 : NDT to endoxifen formation (1/h)
KNDT_ZENDO     : 4.08e-7 : NDT to Z'-endoxifen formation (1/h)
K4OHTAM_ENDO   : 1.81e-3 : 4-OHTAM to endoxifen formation (1/h)
KE_NDT         : 2.46e-3 : NDT elimination rate constant (1/h)
KE_ENDO        : 7.93e-3 : Endoxifen elimination rate constant (1/h)
KE_ZENDO       : 7.93e-3 : Z'-endoxifen elimination rate constant (1/h)
KE_4OHTAM      : 7.93e-3 : 4-OHTAM elimination rate constant (1/h)
KE_4POHTAM     : 7.93e-3 : 4'-OHTAM elimination rate constant (1/h)
KE_NOXTAM      : 7.93e-3 : NOX-TAM elimination rate constant (1/h)

$CMT @annotated
GUT     : Oral absorption depot (mg)
TAM     : Tamoxifen amount (mg)
NDT     : N-desmethyltamoxifen amount (mg equivalent)
OHTAM   : 4-hydroxytamoxifen amount (mg equivalent)
ENDO    : Endoxifen amount (mg equivalent)
ZENDO   : Z'-endoxifen amount (mg equivalent)
NOXTAM  : Tamoxifen N-oxide amount (mg equivalent)
OHTAM4P : 4'-hydroxytamoxifen amount (mg equivalent)

$ODE
double formation_total = KTAM_NDT + KTAM_4OHTAM + KTAM_4POHTAM + KTAM_NOXTAM;

dxdt_GUT = -KA * GUT;
dxdt_TAM = KA * GUT - formation_total * TAM;
dxdt_NDT = KTAM_NDT * TAM - (KNDT_ENDO + KNDT_ZENDO + KE_NDT) * NDT;
dxdt_OHTAM = KTAM_4OHTAM * TAM - (K4OHTAM_ENDO + KE_4OHTAM) * OHTAM;
dxdt_ENDO = KNDT_ENDO * NDT + K4OHTAM_ENDO * OHTAM - KE_ENDO * ENDO;
dxdt_ZENDO = KNDT_ZENDO * NDT - KE_ZENDO * ZENDO;
dxdt_NOXTAM = KTAM_NOXTAM * TAM - KE_NOXTAM * NOXTAM;
dxdt_OHTAM4P = KTAM_4POHTAM * TAM - KE_4POHTAM * OHTAM4P;

$TABLE
double TAM_NM = (TAM / VTAM) * 1000000 / 371.52;
double NDT_NM = (NDT / VTAM) * 1000000 / 357.49;
double OHTAM_NM = (OHTAM / VTAM) * 1000000 / 387.52;
double ENDO_NM = (ENDO / VTAM) * 1000000 / 373.49;
double ZENDO_NM = (ZENDO / VTAM) * 1000000 / 373.49;
double NOXTAM_NM = (NOXTAM / VTAM) * 1000000 / 387.52;
double OHTAM4P_NM = (OHTAM4P / VTAM) * 1000000 / 387.52;
double AAS = TAM_NM + NDT_NM + 100 * OHTAM_NM + 100 * ENDO_NM;

$CAPTURE TAM_NM NDT_NM OHTAM_NM ENDO_NM ZENDO_NM NOXTAM_NM OHTAM4P_NM AAS
"

mod <- mcode("tamoxifen_metabolite_covariate_model", model_code)

covariate_parameters <- function(
    phenotype = "Normal metabolizer",
    diplotype = "Use phenotype",
    inhibitor = "None",
    cyp3a4_22 = FALSE,
    cyp2c19_2 = FALSE,
    cyp2b6_6_homo = FALSE,
    age = 48,
    weight = 64) {
  pars <- BASELINE

  if (!identical(diplotype, "Use phenotype")) {
    pars$KNDT_ENDO <- K_ENDO_BY_DIPLOTYPE[[diplotype]]
  } else {
    pars$KNDT_ENDO <- K_ENDO_BY_PHENO[[phenotype]]
  }

  if (inhibitor == "Weak/moderate" && phenotype %in% c("Normal metabolizer", "Ultra-rapid metabolizer")) {
    pars$KNDT_ENDO <- pars$KNDT_ENDO * 0.680
  }
  if (inhibitor == "Potent" && phenotype %in% c("Normal metabolizer", "Ultra-rapid metabolizer")) {
    pars$KNDT_ENDO <- pars$KNDT_ENDO * 0.434
  }

  if (phenotype %in% c("Intermediate metabolizer", "Poor metabolizer")) {
    pars$KTAM_4OHTAM <- pars$KTAM_4OHTAM * 0.768
  }
  if (cyp3a4_22) {
    pars$KTAM_NDT <- pars$KTAM_NDT * 0.773
    pars$KE_NDT <- pars$KE_NDT * 0.812
  }
  if (cyp2c19_2) {
    pars$KTAM_4OHTAM <- pars$KTAM_4OHTAM * 0.866
  }
  if (cyp2b6_6_homo) {
    pars$KTAM_NOXTAM <- pars$KTAM_NOXTAM * 0.766
  }

  pars$KTAM_NDT <- pars$KTAM_NDT * (age / 48)^(-0.298)
  pars$KTAM_4OHTAM <- pars$KTAM_4OHTAM * (age / 48)^(-0.547)
  pars$KTAM_NOXTAM <- pars$KTAM_NOXTAM * (age / 48)^(-0.296)
  pars$KNDT_ENDO <- pars$KNDT_ENDO * (age / 48)^(-0.480)
  pars$KE_NDT <- pars$KE_NDT * (weight / 64)^(0.245)

  unlist(pars)
}

simulate_tamoxifen <- function(
    dose_mg = 20,
    tau_h = 24,
    days = 180,
    phenotype = "Normal metabolizer",
    diplotype = "Use phenotype",
    inhibitor = "None",
    cyp3a4_22 = FALSE,
    cyp2c19_2 = FALSE,
    cyp2b6_6_homo = FALSE,
    age = 48,
    weight = 64) {
  pars <- covariate_parameters(
    phenotype = phenotype,
    diplotype = diplotype,
    inhibitor = inhibitor,
    cyp3a4_22 = cyp3a4_22,
    cyp2c19_2 = cyp2c19_2,
    cyp2b6_6_homo = cyp2b6_6_homo,
    age = age,
    weight = weight
  )

  ev_dose <- ev(amt = dose_mg, cmt = 1, ii = tau_h, addl = ceiling(days * 24 / tau_h) - 1)

  mod |>
    param(as.list(pars)) |>
    ev(ev_dose) |>
    mrgsim(end = days * 24, delta = 12) |>
    as.data.frame() |>
    mutate(
      day = time / 24,
      ENDO_ng_ml = ENDO_NM * MW[["ENDO"]] / 1000,
      TAM_ng_ml = TAM_NM * MW[["TAM"]] / 1000
    )
}

steady_state_window <- function(sim, tau_h = 24) {
  max_time <- max(sim$time, na.rm = TRUE)
  dplyr::filter(sim, time >= max_time - tau_h)
}

summary_metrics <- function(sim, tau_h = 24) {
  ss <- steady_state_window(sim, tau_h)
  data.frame(
    endo_trough = min(ss$ENDO_NM, na.rm = TRUE),
    endo_peak = max(ss$ENDO_NM, na.rm = TRUE),
    tam_trough = min(ss$TAM_NM, na.rm = TRUE),
    aas_trough = min(ss$AAS, na.rm = TRUE),
    aas_peak = max(ss$AAS, na.rm = TRUE),
    time_to_steady_days = 5 * log(2) / BASELINE$KE_ENDO / 24
  )
}

build_dose_grid <- function(phenotype, inhibitor, age, weight, cyp3a4_22, cyp2c19_2, cyp2b6_6_homo) {
  doses <- c(20, 40, 60, 80)
  lapply(doses, function(dose) {
    sim <- simulate_tamoxifen(
      dose_mg = dose,
      phenotype = phenotype,
      inhibitor = inhibitor,
      age = age,
      weight = weight,
      cyp3a4_22 = cyp3a4_22,
      cyp2c19_2 = cyp2c19_2,
      cyp2b6_6_homo = cyp2b6_6_homo
    )
    metrics <- summary_metrics(sim)
    data.frame(
      dose = dose,
      endo = metrics$endo_trough,
      aas = metrics$aas_trough
    )
  }) |>
    bind_rows()
}

app_theme <- bs_theme(
  version = 5,
  bootswatch = "flatly",
  primary = "#2563eb"
) |> bs_add_rules("
  .metric-card {
    background: #f8fafc;
    border-radius: 8px;
    padding: 14px;
    border: 1px solid #dbe3ef;
    min-height: 96px;
  }
  .metric-value {
    font-size: 22px;
    font-weight: 700;
    color: #0f172a;
  }
  .metric-label {
    font-size: 12px;
    color: #64748b;
  }
  .metric-good .metric-value { color: #047857; }
  .metric-warn .metric-value { color: #b45309; }
  .ref-box {
    background: #eff6ff;
    border-left: 4px solid #2563eb;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
  .limit-box {
    background: #fff7ed;
    border-left: 4px solid #ea580c;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
  .pkpd-footer {
    margin-top: 24px;
    padding-top: 12px;
    border-top: 1px solid #e2e8f0;
    color: #64748b;
    font-size: 12px;
  }
")

ui <- page_sidebar(
  title = "Tamoxifen CYP2D6 Dose Adaptation Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Patient and regimen",
    width = 340,
    sliderInput("dose", "Tamoxifen dose (mg/day)", min = 20, max = 80, value = 20, step = 20),
    selectInput(
      "phenotype",
      "CYP2D6 phenotype",
      choices = names(K_ENDO_BY_PHENO),
      selected = "Normal metabolizer"
    ),
    selectInput(
      "diplotype",
      "Optional CYP2D6 diplotype",
      choices = c("Use phenotype", names(K_ENDO_BY_DIPLOTYPE)),
      selected = "Use phenotype"
    ),
    selectInput(
      "inhibitor",
      "CYP2D6 inhibitor",
      choices = c("None", "Weak/moderate", "Potent"),
      selected = "None"
    ),
    sliderInput("age", "Age (years)", min = 25, max = 84, value = 48, step = 1),
    sliderInput("weight", "Body weight (kg)", min = 40, max = 131, value = 64, step = 1),
    checkboxInput("cyp3a4_22", "CYP3A4*22 carrier", value = FALSE),
    checkboxInput("cyp2c19_2", "CYP2C19*2 carrier", value = FALSE),
    checkboxInput("cyp2b6_6_homo", "CYP2B6*6/*6", value = FALSE),
    hr(),
    checkboxInput("log_scale", "Log concentration axis", value = FALSE)
  ),
  layout_column_wrap(
    width = 1 / 4,
    uiOutput("endo_card"),
    uiOutput("aas_card"),
    uiOutput("tam_card"),
    uiOutput("ss_card")
  ),
  card(
    full_screen = TRUE,
    card_header("Steady-state concentration-time profiles"),
    plotOutput("pk_plot", height = "500px")
  ),
  navset_card_underline(
    nav_panel(
      "Dose Targets",
      layout_columns(
        col_widths = c(7, 5),
        plotOutput("dose_plot", height = "380px"),
        tableOutput("dose_table")
      )
    ),
    nav_panel(
      "Model Information",
      div(
        class = "ref-box",
        tags$b("Model basis: "),
        "Puszkiel et al. CPT 2021 Table 2. The app implements a seven-compartment parent-metabolite structure for TAM, NDT, 4-OHTAM, ENDO, Z'-ENDO, NOX-TAM, and 4'-OHTAM using the published population rate constants and covariate multipliers."
      ),
      div(
        class = "limit-box",
        tags$b("Important limitation: "),
        "This is a simplified educational simulator built from published aggregate parameters. It does not reproduce the full NONMEM model, IIV, residual variability, genotype imputation, or the original 1,000-subject simulation uncertainty."
      ),
      tags$ul(
        tags$li("ENDO target lines: 16 nM therapeutic threshold and 40 nM median reference in CYP2D6 normal metabolizers."),
        tags$li("AAS uses TAM + NDT + 100 x 4-OHTAM + 100 x ENDO as a practical reconstruction of the published antiestrogenic activity score concept."),
        tags$li("Covariates included: CYP2D6 phenotype/diplotype, CYP2D6 inhibitors, CYP3A4*22, CYP2C19*2, CYP2B6*6/*6, age, and body weight.")
      )
    ),
    nav_panel(
      "References",
      tags$p(
        "Puszkiel A, Arellano C, Vachoux C, et al. Model-Based Quantification of Impact of Genetic Polymorphisms and Co-Medications on Pharmacokinetics of Tamoxifen and Six Metabolites in Breast Cancer. Clinical Pharmacology & Therapeutics. 2021;109(5):1244-1255."
      ),
      tags$p(
        "de Vries Schultink AHM, et al. An Antiestrogenic Activity Score for tamoxifen and its metabolites is associated with breast cancer outcome. Breast Cancer Research and Treatment. 2017;161:567-574."
      )
    )
  ),
  div(class = "pkpd-footer", "Generated by PKPDBuilder")
)

metric_card <- function(value, label, class = "") {
  div(
    class = paste("metric-card", class),
    div(class = "metric-value", value),
    div(class = "metric-label", label)
  )
}

server <- function(input, output, session) {
  sim_data <- reactive({
    shiny::req(input$dose, input$phenotype, input$diplotype, input$inhibitor, input$age, input$weight)
    simulate_tamoxifen(
      dose_mg = input$dose,
      phenotype = input$phenotype,
      diplotype = input$diplotype,
      inhibitor = input$inhibitor,
      age = input$age,
      weight = input$weight,
      cyp3a4_22 = input$cyp3a4_22,
      cyp2c19_2 = input$cyp2c19_2,
      cyp2b6_6_homo = input$cyp2b6_6_homo
    )
  })

  metrics <- reactive({
    summary_metrics(sim_data())
  })

  output$endo_card <- renderUI({
    m <- metrics()
    cls <- if (m$endo_trough >= 16) "metric-good" else "metric-warn"
    metric_card(sprintf("%.1f nM", m$endo_trough), "Endoxifen trough", cls)
  })

  output$aas_card <- renderUI({
    m <- metrics()
    cls <- if (m$aas_trough >= 1798) "metric-good" else "metric-warn"
    metric_card(sprintf("%.0f", m$aas_trough), "AAS trough", cls)
  })

  output$tam_card <- renderUI({
    m <- metrics()
    metric_card(sprintf("%.0f nM", m$tam_trough), "Tamoxifen trough")
  })

  output$ss_card <- renderUI({
    m <- metrics()
    metric_card(sprintf("%.0f days", m$time_to_steady_days), "Approx. 5 half-lives")
  })

  output$pk_plot <- renderPlot({
    plot_data <- sim_data() |>
      dplyr::filter(time >= max(time) - 24 * 14) |>
      select(time, TAM_NM, NDT_NM, OHTAM_NM, ENDO_NM) |>
      pivot_longer(-time, names_to = "compound", values_to = "concentration_nM") |>
      mutate(
        day = time / 24,
        compound = recode(
          compound,
          TAM_NM = "Tamoxifen",
          NDT_NM = "N-desmethyltamoxifen",
          OHTAM_NM = "4-hydroxytamoxifen",
          ENDO_NM = "Endoxifen"
        )
      )

    p <- ggplot(plot_data, aes(x = day, y = concentration_nM, color = compound)) +
      geom_line(linewidth = 1) +
      geom_hline(yintercept = 16, linetype = "dashed", color = "#b45309") +
      geom_hline(yintercept = 40, linetype = "dotted", color = "#047857") +
      labs(
        x = "Time (days)",
        y = "Concentration (nM)",
        color = NULL
      ) +
      theme_minimal(base_size = 13) +
      theme(legend.position = "bottom")

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

    p
  })

  dose_grid <- reactive({
    build_dose_grid(
      phenotype = input$phenotype,
      inhibitor = input$inhibitor,
      age = input$age,
      weight = input$weight,
      cyp3a4_22 = input$cyp3a4_22,
      cyp2c19_2 = input$cyp2c19_2,
      cyp2b6_6_homo = input$cyp2b6_6_homo
    )
  })

  output$dose_plot <- renderPlot({
    plot_data <- dose_grid() |>
      pivot_longer(c(endo, aas), names_to = "metric", values_to = "value") |>
      mutate(metric = recode(metric, endo = "Endoxifen trough (nM)", aas = "AAS trough"))

    targets <- data.frame(
      metric = c("Endoxifen trough (nM)", "Endoxifen trough (nM)", "AAS trough", "AAS trough"),
      target = c("16 nM", "40 nM", "1798", "3753"),
      value = c(16, 40, 1798, 3753)
    )

    ggplot(plot_data, aes(x = dose, y = value)) +
      geom_line(color = "#2563eb", linewidth = 1) +
      geom_point(color = "#2563eb", size = 3) +
      geom_hline(data = targets, aes(yintercept = value, linetype = target), color = "#64748b") +
      facet_wrap(~metric, scales = "free_y") +
      scale_x_continuous(breaks = c(20, 40, 60, 80)) +
      labs(x = "Dose (mg/day)", y = NULL, linetype = "Reference") +
      theme_minimal(base_size = 13) +
      theme(legend.position = "bottom")
  })

  output$dose_table <- renderTable({
    table_data <- dose_grid() |>
      mutate(
        dose_label = dose,
        endo_label = sprintf("%.1f", endo),
        aas_label = sprintf("%.0f", aas),
        endo_target = ifelse(endo >= 16, "Yes", "No"),
        aas_target = ifelse(aas >= 1798, "Yes", "No")
      ) |>
      select(dose_label, endo_label, aas_label, endo_target, aas_target)

    names(table_data) <- c(
      "Dose (mg/day)",
      "Endoxifen trough (nM)",
      "AAS trough",
      "ENDO >=16 nM",
      "AAS >=1798"
    )
    table_data
  }, striped = TRUE, bordered = TRUE, spacing = "s")
}

shinyApp(ui, server)