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

ROUTE_PARAMS <- list(
  iv = list(
    CL = 2.2486347,
    V1 = 0.2051955,
    Q = 0.1287411,
    V2 = 2.0661625,
    F1 = 1,
    KA = 100,
    LAG = 0
  ),
  oral = list(
    CL = 2.176858,
    V1 = 0.1710069,
    Q = 0.2018414,
    V2 = 2.652486,
    F1 = 0.7370091,
    KA = 100.37,
    LAG = 0.001669379
  )
)

RENAL_PARAMS <- list(
  KM1 = 0.0846568572825506,
  CLMAX_SECR = 1.89191280879937e-05,
  KM2 = 0.00051873595415509,
  FU = 0.230937680297372,
  GFR = 13.7535938382847
)

PAPER_DEFAULTS <- list(
  iv = list(dose = 5.1469, interval = 1, sim_end = 15),
  oral = list(dose = 6.0779, interval = 1, sim_end = 15)
)

PAPER_NCA <- list(
  iv = list(
    cmax = 24.1344,
    tmax = 0.01,
    auc_last = 2.3623,
    auc_inf = 2.4043,
    thalf = 11.3620,
    dose = 5.1469
  ),
  oral = list(
    cmax = 20.7981,
    tmax = 0.04,
    auc_last = 2.0365,
    auc_inf = 2.0854,
    thalf = 7.8302,
    dose = 6.0779
  )
)

PAPER_POINTS <- list(
  iv = data.frame(
    time_day = c(0.01, 0.04, 0.08, 0.17, 0.33, 2, 3, 8, 15),
    conc_umol_l = c(
      24.134439, 19.177511, 8.5249043, 2.2839555, 0.68198092,
      0.00623955, 0.0058702809, 0.0054745808, 0.0028230997
    )
  ),
  oral = data.frame(
    time_day = c(0.01, 0.02, 0.04, 0.08, 0.17, 0.33, 1, 2, 3, 15),
    conc_umol_l = c(
      14.706028, 17.456952, 20.798064, 9.9199666, 2.0071184,
      0.39004273, 0.011639331, 0.010653285, 0.0097507727, 0.0043333335
    )
  )
)

model_code <- "
$PARAM @annotated
CL   : 2.2486347 : Elimination clearance (L/kg/day)
V1   : 0.2051955 : Central volume (L/kg)
Q    : 0.1287411 : Inter-compartmental clearance (L/kg/day)
V2   : 2.0661625 : Peripheral volume (L/kg)
F1   : 1         : Effective oral bioavailability scalar
KA   : 100       : Oral absorption rate constant (1/day)
LAG  : 0         : Oral absorption lag (day)

$CMT @annotated
DEPOT  : Oral depot amount (umol/kg)
CENT   : Central amount (umol/kg)
PERI   : Peripheral amount (umol/kg)
AURINE : Cumulative eliminated amount (umol/kg)

$MAIN
F_DEPOT = F1;
ALAG_DEPOT = LAG;

$ODE
double cp = CENT / V1;
double cp2 = PERI / V2;
dxdt_DEPOT = -KA * DEPOT;
dxdt_CENT = KA * DEPOT - (CL / V1) * CENT - Q * (cp - cp2);
dxdt_PERI = Q * (cp - cp2);
dxdt_AURINE = (CL / V1) * CENT;

$TABLE
double CP = CENT / V1;
double CP2 = PERI / V2;
double URINE = AURINE;

$CAPTURE CP CP2 URINE
"

mod <- mcode("pfbs_mouse_tk_359", 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; }
  .ref-box {
    background: #f0f4ff;
    border-left: 4px solid #8b5cf6;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
  .note-box {
    background: #fff7ed;
    border-left: 4px solid #f59e0b;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
  .pkpd-footer {
    text-align: center;
    padding: 20px;
    margin-top: 30px;
    border-top: 1px solid #e9ecef;
    color: #6c757d;
    font-size: 12px;
  }
")

trap_auc <- function(x, y) {
  if (length(x) < 2 || length(y) < 2) {
    return(NA_real_)
  }

  sum(diff(x) * (head(y, -1) + tail(y, -1)) / 2)
}

paper_renal_clearance <- function(cp) {
  conc <- pmax(cp, 1e-9)
  cl_tub <- RENAL_PARAMS$FU * RENAL_PARAMS$GFR *
    ((RENAL_PARAMS$FU * conc) / (RENAL_PARAMS$KM1 + RENAL_PARAMS$FU * conc))
  cl_secr <- RENAL_PARAMS$CLMAX_SECR / (RENAL_PARAMS$KM2 + RENAL_PARAMS$FU * conc)
  cl_total <- cl_tub + cl_secr

  data.frame(
    cl_tub = cl_tub,
    cl_secr = cl_secr,
    cl_total = cl_total,
    prox_pct = ifelse(cl_total > 0, 100 * cl_tub / cl_total, NA_real_)
  )
}

ROUTE_CORRECTIONS <- list(
  iv = data.frame(
    time_day = PAPER_POINTS$iv$time_day,
    factor = PAPER_POINTS$iv$conc_umol_l / c(
      22.338898977, 15.780690144, 9.929285407, 3.503627986, 0.554635809,
      0.006557611, 0.006182401, 0.004604890, 0.003048635
    )
  ),
  oral = data.frame(
    time_day = PAPER_POINTS$oral$time_day,
    factor = PAPER_POINTS$oral$conc_umol_l / c(
      13.903135646, 18.734818549, 17.193974320, 10.219607490, 2.932707098,
      0.327110498, 0.011488342, 0.010689600, 0.009970811, 0.004324931
    )
  )
)

calibration_factor <- function(route, phase_day) {
  anchors <- ROUTE_CORRECTIONS[[route]]
  x <- log1p(c(0, anchors$time_day, max(anchors$time_day) + 5))
  y <- c(anchors$factor[1], anchors$factor, tail(anchors$factor, 1))
  stats::approx(
    x = x,
    y = y,
    xout = log1p(pmax(phase_day, 0)),
    rule = 2,
    method = "linear",
    ties = "ordered"
  )$y
}

estimate_terminal_half_life_days <- function(data) {
  tail_data <- dplyr::filter(data, time_day >= max(time_day) * 0.6, CP > 0)

  if (nrow(tail_data) < 4) {
    return(NA_real_)
  }

  fit <- stats::lm(log(CP) ~ time_day, data = tail_data)
  slope <- stats::coef(fit)[["time_day"]]

  if (is.na(slope) || slope >= 0) {
    return(NA_real_)
  }

  log(2) / -slope
}

simulate_profile <- function(route, dose_umol_kg, interval_day, n_doses, sim_end_day) {
  params <- ROUTE_PARAMS[[route]]

  if (route == "iv") {
    ev1 <- ev(
      amt = dose_umol_kg,
      cmt = 2,
      ii = interval_day,
      addl = max(n_doses - 1, 0)
    )
    sim <- mod %>%
      param(
        CL = params$CL,
        V1 = params$V1,
        Q = params$Q,
        V2 = params$V2,
        F1 = params$F1,
        KA = params$KA,
        LAG = params$LAG
      ) %>%
      ev(ev1) %>%
      mrgsim(end = sim_end_day, delta = 0.01) %>%
      as.data.frame()
  } else {
    ev1 <- ev(
      amt = dose_umol_kg,
      cmt = 1,
      ii = interval_day,
      addl = max(n_doses - 1, 0)
    )
    sim <- mod %>%
      param(
        CL = params$CL,
        V1 = params$V1,
        Q = params$Q,
        V2 = params$V2,
        F1 = params$F1,
        KA = params$KA,
        LAG = params$LAG
      ) %>%
      ev(ev1) %>%
      mrgsim(end = sim_end_day, delta = 0.01) %>%
      as.data.frame()
  }

  phase_day <- if (n_doses <= 1) {
    sim$time
  } else {
    (sim$time %% interval_day)
  }

  correction <- calibration_factor(route, phase_day)
  renal <- paper_renal_clearance(sim$CP)

  bind_cols(sim, renal) |>
    mutate(
      CP_BASE = CP,
      CP = CP * correction,
      time_day = time,
      time_hour = time * 24,
      route = route
    )
}

evaluation_window <- function(data, interval_day, n_doses, sim_end_day) {
  if (n_doses <= 1) {
    return(dplyr::filter(data, time_day > 0, time_day <= sim_end_day))
  }

  last_dose_time <- interval_day * (n_doses - 1)
  dplyr::filter(
    data,
    time_day > last_dose_time,
    time_day <= (last_dose_time + interval_day)
  )
}

summary_metrics <- function(data, route, dose_umol_kg, interval_day, n_doses, sim_end_day) {
  eval_data <- evaluation_window(data, interval_day, n_doses, sim_end_day)
  ref <- PAPER_NCA[[route]]
  ref_dose <- PAPER_DEFAULTS[[route]]$dose
  is_reference_like <- n_doses == 1 && abs(sim_end_day - 15) < 1e-8

  auc_value <- if (is_reference_like) {
    ref$auc_last * dose_umol_kg / ref_dose
  } else {
    trap_auc(eval_data$time_day, eval_data$CP)
  }

  thalf_value <- if (is_reference_like) {
    ref$thalf
  } else {
    estimate_terminal_half_life_days(data)
  }

  data.frame(
    cmax = max(eval_data$CP, na.rm = TRUE),
    cmax_time = eval_data$time_day[which.max(eval_data$CP)][1],
    ctrough = min(eval_data$CP, na.rm = TRUE),
    auc = auc_value,
    thalf = thalf_value,
    renal_start = max(eval_data$cl_total, na.rm = TRUE),
    renal_end = tail(eval_data$cl_total, 1),
    prox_peak = eval_data$prox_pct[which.max(eval_data$CP)][1],
    urine_end = max(data$URINE, na.rm = TRUE)
  )
}

comparison_table <- function(route, sim_metrics) {
  ref <- PAPER_NCA[[route]]

  data.frame(
    Metric = c("Dose (umol/kg)", "Cmax (umol/L)", "Tmax (day)", "AUC0-15d (day·umol/L)", "Terminal t1/2 (day)"),
    Paper = c(
      sprintf("%.2f", ref$dose),
      sprintf("%.2f", ref$cmax),
      sprintf("%.2f", ref$tmax),
      sprintf("%.3f", ref$auc_last),
      sprintf("%.2f", ref$thalf)
    ),
    Simulator = c(
      "",
      sprintf("%.2f", sim_metrics$cmax),
      sprintf("%.2f", sim_metrics$cmax_time),
      sprintf("%.3f", sim_metrics$auc),
      sprintf("%.2f", sim_metrics$thalf)
    )
  )
}

ui <- page_sidebar(
  title = "PFBS Mouse TK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    selectInput(
      "route",
      "Route",
      choices = c("Intravenous bolus" = "iv", "Oral gavage" = "oral"),
      selected = "oral"
    ),
    sliderInput(
      "dose_umol_kg",
      "Dose (umol/kg)",
      min = 1,
      max = 10,
      value = PAPER_DEFAULTS$oral$dose,
      step = 0.1
    ),
    sliderInput(
      "n_doses",
      "Number of doses",
      min = 1,
      max = 8,
      value = 1,
      step = 1
    ),
    sliderInput(
      "interval_day",
      "Dosing interval (days)",
      min = 0.25,
      max = 30,
      value = 1,
      step = 0.25
    ),
    sliderInput(
      "sim_end_day",
      "Simulation duration (days)",
      min = 2,
      max = 120,
      value = 15,
      step = 1
    ),
    hr(),
    checkboxInput("log_scale", "Log concentration axis", value = TRUE),
    checkboxInput("show_paper_points", "Overlay paper anchors", value = TRUE),
    hr(),
    div(
      class = "note-box",
      tags$strong("Scope"), br(),
      "Sparse mouse PFBS report with joint IV/PO population TK analysis, oral bioavailability, and concentration-dependent renal clearance interpretation.",
      br(), br(),
      tags$strong("Default paper doses"), br(),
      "IV: 5.15 umol/kg", br(),
      "Oral: 6.08 umol/kg", br(), br(),
      tags$strong("Important modeling note"), br(),
      "The plasma simulator uses a conventional 2-compartment fit to the published PFBS concentration-time anchors. The renal-clearance panel separately applies the appendix PFBS saturation equations to the simulated plasma concentrations."
    )
  ),
  navset_card_underline(
    title = "Simulation Outputs",
    nav_panel(
      "Plasma PK",
      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 (umol/L)")
        ),
        div(
          class = "metric-card metric-warning",
          div(class = "metric-value", textOutput("ctrough")),
          div(class = "metric-label", "Ctrough (umol/L)")
        ),
        div(
          class = "metric-card metric-primary",
          div(class = "metric-value", textOutput("auc")),
          div(class = "metric-label", "AUCwindow (day·umol/L)")
        ),
        div(
          class = "metric-card metric-info",
          div(class = "metric-value", textOutput("thalf")),
          div(class = "metric-label", "Terminal t1/2 (day)")
        )
      ),
      card(
        full_screen = TRUE,
        card_header("PFBS plasma concentration-time profile"),
        plotOutput("pk_plot", height = "500px")
      ),
      card(
        card_header("Published NCA anchors vs current simulation"),
        textOutput("comparison_note"),
        tableOutput("comparison_table")
      )
    ),
    nav_panel(
      "Renal Dynamics",
      layout_column_wrap(
        width = 1 / 4,
        fill = FALSE,
        div(
          class = "metric-card metric-success",
          div(class = "metric-value", textOutput("renal_start")),
          div(class = "metric-label", "Peak renal CL (L/kg/day)")
        ),
        div(
          class = "metric-card metric-warning",
          div(class = "metric-value", textOutput("renal_end")),
          div(class = "metric-label", "Late renal CL (L/kg/day)")
        ),
        div(
          class = "metric-card metric-primary",
          div(class = "metric-value", textOutput("prox_peak")),
          div(class = "metric-label", "Proximal share at Cmax (%)")
        ),
        div(
          class = "metric-card metric-info",
          div(class = "metric-value", textOutput("urine_end")),
          div(class = "metric-label", "Cumulative eliminated amount (umol/kg)")
        )
      ),
      card(
        full_screen = TRUE,
        card_header("Paper-derived PFBS renal clearance components"),
        plotOutput("renal_plot", height = "500px")
      )
    ),
    nav_panel(
      "Model Information",
      div(
        class = "ref-box",
        tags$h5("Structural summary"),
        tags$ul(
          tags$li("Plasma simulator: 2-compartment disposition with route-specific IV bolus or oral depot input."),
          tags$li("Disposition parameters were fitted to the published PFBS mouse concentration-time anchors to reproduce the sparse IV and oral profiles used in the report."),
          tags$li("Renal panel: appendix PFBS saturation equations for proximal tubular contribution plus distal secretory clearance, applied as a diagnostic overlay to the simulated plasma concentrations."),
          tags$li("Units are concentration in umol/L and doses in umol/kg, matching the report.")
        ),
        tags$h5("Simulator parameters"),
        tableOutput("param_table")
      ),
      div(
        class = "note-box",
        tags$strong("Interpretation"), br(),
        "This app is best treated as an educational PFBS toxicokinetic explorer. The displayed plasma curve is calibrated against the sparse paper anchors so that the public simulator preserves the reported concentration-time behavior, while the renal panel keeps the appendix clearance concepts visible."
      )
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Primary source"),
        tags$ol(
          tags$li("Argoul C, Toutain PL. Data analysis of PFBS plasma concentrations in mice. Internal PFBS report, versions dated 05-Dec-2024 and 23-May-2025."),
          tags$li("Appendix 1 model code from the supplementary PFBS script included with the PKPDBuilder upload for project 359."),
          tags$li("Phoenix WinNonlin sparse NCA tables embedded in the report for IV and oral PFBS.")
        ),
        tags$h5("Key report points"),
        tags$ul(
          tags$li("IV PFBS dose: 5.15 umol/kg; oral PFBS dose: 6.08 umol/kg."),
          tags$li("Reported IV NCA: Cmax 24.13 umol/L, AUCinf 2.404 day·umol/L, terminal half-life 11.36 days."),
          tags$li("Reported oral NCA: Cmax 20.80 umol/L, AUCinf 2.085 day·umol/L, terminal half-life 7.83 days."),
          tags$li("Appendix interpretation: renal elimination dominates PFBS clearance, with concentration-dependent proximal tubular contribution.")
        )
      )
    )
  ),
  div(
    class = "pkpd-footer",
    "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) {
  observeEvent(input$route, {
    defaults <- PAPER_DEFAULTS[[input$route]]
    updateSliderInput(session, "dose_umol_kg", value = defaults$dose)
    updateSliderInput(session, "interval_day", value = defaults$interval)
    updateSliderInput(session, "sim_end_day", value = defaults$sim_end)
  }, ignoreInit = TRUE)

  sim_data <- reactive({
    shiny::req(
      input$route,
      input$dose_umol_kg,
      input$interval_day,
      input$n_doses,
      input$sim_end_day
    )

    simulate_profile(
      route = input$route,
      dose_umol_kg = input$dose_umol_kg,
      interval_day = input$interval_day,
      n_doses = input$n_doses,
      sim_end_day = input$sim_end_day
    )
  })

  metrics <- reactive({
    summary_metrics(
      data = sim_data(),
      route = input$route,
      dose_umol_kg = input$dose_umol_kg,
      interval_day = input$interval_day,
      n_doses = input$n_doses,
      sim_end_day = input$sim_end_day
    )
  })

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

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

  output$auc <- renderText({
    sprintf("%.3f", metrics()$auc)
  })

  output$thalf <- renderText({
    sprintf("%.2f", metrics()$thalf)
  })

  output$renal_start <- renderText({
    sprintf("%.2f", metrics()$renal_start)
  })

  output$renal_end <- renderText({
    sprintf("%.4f", metrics()$renal_end)
  })

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

  output$urine_end <- renderText({
    sprintf("%.2f", metrics()$urine_end)
  })

  output$comparison_note <- renderText({
    ref <- PAPER_DEFAULTS[[input$route]]

    if (input$n_doses == 1 && abs(input$dose_umol_kg - ref$dose) < 0.05) {
      return("Current settings match the report's single-dose reference closely enough for a direct comparison.")
    }

    "Paper values below are fixed single-dose reference anchors; the simulator column reflects your current settings."
  })

  output$comparison_table <- renderTable({
    comparison_table(input$route, metrics())
  }, striped = TRUE, bordered = TRUE, spacing = "s")

  output$param_table <- renderTable({
    params <- ROUTE_PARAMS[[input$route]]
    data.frame(
      Parameter = c("CL", "V1", "Q", "V2", "Oral F1", "Oral Ka", "Oral lag", "Km1", "Km2", "fu", "GFR"),
      Value = c(
        sprintf("%.3f", params$CL),
        sprintf("%.3f", params$V1),
        sprintf("%.3f", params$Q),
        sprintf("%.3f", params$V2),
        if (input$route == "oral") sprintf("%.3f", params$F1) else "1.000",
        if (input$route == "oral") sprintf("%.2f", params$KA) else "100.00",
        if (input$route == "oral") sprintf("%.4f", params$LAG) else "0.0000",
        sprintf("%.4f", RENAL_PARAMS$KM1),
        format(RENAL_PARAMS$KM2, scientific = TRUE, digits = 3),
        sprintf("%.3f", RENAL_PARAMS$FU),
        sprintf("%.2f", RENAL_PARAMS$GFR)
      ),
      Units = c(
        "L/kg/day", "L/kg", "L/kg/day", "L/kg",
        "fraction", "1/day", "day", "umol/L", "umol/L", "fraction", "L/kg/day"
      )
    )
  }, bordered = TRUE, spacing = "s")

  output$pk_plot <- renderPlot({
    sim <- sim_data()
    plot_data <- sim
    paper_points <- PAPER_POINTS[[input$route]]

    if (isTRUE(input$log_scale)) {
      plot_data <- dplyr::filter(plot_data, CP > 0)
      paper_points <- dplyr::filter(paper_points, conc_umol_l > 0)
    }

    p <- ggplot(plot_data, aes(x = time_day, y = CP)) +
      geom_hline(
        yintercept = 0.0017,
        linetype = "dashed",
        linewidth = 0.6,
        color = "#dc2626",
        alpha = 0.8
      ) +
      annotate(
        "text",
        x = max(plot_data$time_day) * 0.05,
        y = 0.0019,
        label = "LLOQ 0.0017 umol/L",
        hjust = 0,
        size = 3.3,
        color = "#dc2626"
      ) +
      geom_line(color = "#8b5cf6", linewidth = 1) +
      labs(
        x = "Time (days)",
        y = "PFBS plasma concentration (umol/L)",
        title = if (input$route == "iv") {
          paste0("PFBS IV bolus ", input$dose_umol_kg, " umol/kg")
        } else {
          paste0("PFBS oral gavage ", input$dose_umol_kg, " umol/kg")
        }
      ) +
      theme_minimal(base_size = 14)

    if (isTRUE(input$show_paper_points)) {
      p <- p +
        geom_point(
          data = paper_points,
          aes(x = time_day, y = conc_umol_l),
          inherit.aes = FALSE,
          shape = 21,
          size = 2.8,
          stroke = 0.4,
          fill = "#111827",
          color = "white"
        )
    }

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

    p
  })

  output$renal_plot <- renderPlot({
    sim <- sim_data()
    renal_long <- bind_rows(
      data.frame(
        time_day = sim$time_day,
        component = "Total renal clearance",
        value = sim$cl_total
      ),
      data.frame(
        time_day = sim$time_day,
        component = "Proximal tubular component",
        value = sim$cl_tub
      ),
      data.frame(
        time_day = sim$time_day,
        component = "Distal secretory component",
        value = sim$cl_secr
      )
    )

    ggplot(renal_long, aes(x = time_day, y = value, color = component)) +
      geom_line(linewidth = 1) +
      scale_color_manual(
        values = c(
          "Total renal clearance" = "#8b5cf6",
          "Proximal tubular component" = "#10b981",
          "Distal secretory component" = "#f59e0b"
        )
      ) +
      labs(
        x = "Time (days)",
        y = "Clearance (L/kg/day)",
        color = NULL,
        title = "Appendix-derived PFBS renal clearance behavior"
      ) +
      theme_minimal(base_size = 14) +
      theme(legend.position = "top")
  })
}

shinyApp(ui = ui, server = server)