use("shiny")
use("bslib")
use("mrgsolve")
use("dplyr")
use("ggplot2")

model_code <- '
$PARAM @annotated
KA      : 0.0373 : Absorption rate constant (1/h)
CL70    : 0.0329 : Clearance for a 70 kg individual (L/h)
VC70    : 2.47   : Central volume for a 70 kg individual (L)
Q70     : 0.126  : Intercompartmental clearance for a 70 kg individual (L/h)
VP70    : 3.98   : Peripheral volume for a 70 kg individual (L)
FREF    : 0.8    : Absolute bioavailability in reference study
FREL    : -0.181 : Relative bioavailability effect in pooled studies
FRACFMV : 0.482  : Fraction of fat mass contributing to volume scaling
WT      : 90     : Body weight (kg)
BMI     : 32     : Body mass index (kg/m^2)
SEXF    : 1      : Sex flag (1 = female, 0 = male)
FSCEN   : 0      : Bioavailability scenario (0 = pooled studies, 1 = reference study)

$CMT @annotated
DEPOT  : SC depot amount (mg)
CENT   : Central amount (mg)
PERIPH : Peripheral amount (mg)

$MAIN
double FFM;
if (SEXF > 0.5) {
  FFM = (9270.0 * WT) / (8780.0 + 244.0 * BMI);
} else {
  FFM = (9270.0 * WT) / (6680.0 + 216.0 * BMI);
}
double FAT = WT - FFM;
double NFMC = FFM + FAT;
double NFMV = FFM + FRACFMV * FAT;

double CLi = CL70 * pow(NFMC / 70.0, 0.8);
double Qi  = Q70  * pow(NFMC / 70.0, 0.8);
double VCi = VC70 * (NFMV / 70.0);
double VPi = VP70 * (NFMV / 70.0);
double F1i = FREF * (FSCEN > 0.5 ? 1.0 : (1.0 + FREL));

F_DEPOT = F1i;

$ODE
double K10 = CLi / VCi;
double K12 = Qi / VCi;
double K21 = Qi / VPi;

dxdt_DEPOT  = -KA * DEPOT;
dxdt_CENT   = KA * DEPOT - K10 * CENT - K12 * CENT + K21 * PERIPH;
dxdt_PERIPH = K12 * CENT - K21 * PERIPH;

$TABLE
double CP = CENT / VCi;
double CP_ngml = CP > 0 ? CP * 1000.0 : 0.0;
double VSS = VCi + VPi;
double CLF = CLi / F1i;
double VSSF = VSS / F1i;

$CAPTURE @annotated
CP_ngml : Plasma concentration (ng/mL)
CLi     : Clearance (L/h)
Qi      : Intercompartmental clearance (L/h)
VCi     : Central volume (L)
VPi     : Peripheral volume (L)
VSS     : Steady-state volume (L)
F1i     : Bioavailability
CLF     : Apparent clearance (L/h)
VSSF    : Apparent steady-state volume (L)
'

mod <- mrgsolve::mcode("tirzepatide_pop_pk_330", model_code)

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

simulate_tirzepatide <- function(
  dose_mg = 5,
  interval_h = 168,
  n_weeks = 16,
  wt = 90,
  bmi = 32,
  sex = "Female",
  bioavailability = "pooled",
  delta_h = 2
) {
  n_doses <- max(4, ceiling((n_weeks * 168) / interval_h))
  last_dose_time <- (n_doses - 1) * interval_h
  sim_end_h <- last_dose_time + interval_h
  ev_obj <- mrgsolve::ev(amt = dose_mg, cmt = 1, ii = interval_h, addl = n_doses - 1)

  mod %>%
    mrgsolve::param(
      WT = wt,
      BMI = bmi,
      SEXF = ifelse(sex == "Female", 1, 0),
      FSCEN = ifelse(bioavailability == "reference", 1, 0)
    ) %>%
    mrgsolve::ev(ev_obj) %>%
    mrgsolve::mrgsim(end = sim_end_h, delta = delta_h) %>%
    as.data.frame() |>
    dplyr::mutate(
      time_h = time,
      time_days = time / 24,
      last_dose_time = last_dose_time,
      interval_h = interval_h,
      dose_mg = dose_mg,
      sex = sex,
      wt = wt,
      bmi = bmi,
      bioavailability = bioavailability
    )
}

build_weekly_schedule <- function(mode, n_weeks, fixed_dose, maintenance_dose,
                                  weeks_per_step, custom_doses) {
  if (mode == "fixed") {
    return(rep(fixed_dose, n_weeks))
  }

  if (mode == "titration") {
    step_doses <- seq(2.5, maintenance_dose, by = 2.5)
    schedule <- rep(step_doses, each = weeks_per_step)
    return(c(schedule, rep(maintenance_dose, n_weeks)) |> head(n_weeks))
  }

  parsed <- strsplit(custom_doses, "[,[:space:]]+")[[1]]
  parsed <- parsed[nzchar(parsed)] |> as.numeric()
  shiny::validate(
    shiny::need(length(parsed) > 0 && all(is.finite(parsed)),
                "Enter weekly doses separated by commas."),
    shiny::need(all(parsed >= 0 & parsed <= 15),
                "Each weekly dose must be between 0 and 15 mg.")
  )
  parsed
}

simulate_weekly_schedule <- function(doses, wt = 90, bmi = 32, sex = "Female",
                                     bioavailability = "pooled", delta_h = 2) {
  dose_times <- (seq_along(doses) - 1) * 168
  active <- doses > 0
  sim_end_h <- max(length(doses) * 168, 168)
  events <- mrgsolve::ev(
    time = dose_times[active],
    amt = doses[active],
    cmt = 1
  )

  mod %>%
    mrgsolve::param(
      WT = wt,
      BMI = bmi,
      SEXF = ifelse(sex == "Female", 1, 0),
      FSCEN = ifelse(bioavailability == "reference", 1, 0)
    ) %>%
    mrgsolve::ev(events) %>%
    mrgsolve::mrgsim(end = sim_end_h, delta = delta_h) %>%
    as.data.frame() |>
    dplyr::mutate(
      time_h = time,
      time_days = time / 24,
      last_dose_time = sim_end_h - 168,
      interval_h = 168,
      dose_mg = tail(doses, 1),
      sex = sex,
      wt = wt,
      bmi = bmi,
      bioavailability = bioavailability
    )
}

last_interval_metrics <- function(data) {
  last_start <- data$last_dose_time[[1]]
  interval_h <- data$interval_h[[1]]
  last_window <- data |>
    dplyr::filter(time_h >= last_start, time_h <= last_start + interval_h)

  cmax <- max(last_window$CP_ngml, na.rm = TRUE)
  tmax_h <- last_window$time_h[[which.max(last_window$CP_ngml)]] - last_start
  ctrough <- tail(last_window$CP_ngml, 1)
  auc_tau <- sum(
    diff(last_window$time_h) *
      (head(last_window$CP_ngml, -1) + tail(last_window$CP_ngml, -1)) / 2
  )

  vci <- last_window$VCi[[1]]
  vpi <- last_window$VPi[[1]]
  cli <- last_window$CLi[[1]]
  qi <- last_window$Qi[[1]]
  k10 <- cli / vci
  k12 <- qi / vci
  k21 <- qi / vpi
  beta <- 0.5 * ((k10 + k12 + k21) - sqrt((k10 + k12 + k21)^2 - 4 * k21 * k10))
  thalf_days <- log(2) / beta / 24

  list(
    cmax = cmax,
    tmax_h = tmax_h,
    ctrough = ctrough,
    auc_tau = auc_tau,
    thalf_days = thalf_days,
    clf = last_window$CLF[[1]],
    vssf = last_window$VSSF[[1]],
    f1 = last_window$F1i[[1]]
  )
}

paper_parameters <- data.frame(
  parameter = c(
    "Ka",
    "CL/70 kg",
    "Q/70 kg",
    "Vc/70 kg",
    "Vp/70 kg",
    "Bioavailability",
    "Fat-mass fraction on Vd",
    "Half-life",
    "Apparent CL (T2DM mean)",
    "Apparent Vd (T2DM mean)"
  ),
  estimate = c(
    "0.0373",
    "0.0329",
    "0.126",
    "2.47",
    "3.98",
    "0.8 fixed",
    "0.482",
    "5.4",
    "0.061",
    "10.3"
  ),
  units = c(
    "1/h",
    "L/h",
    "L/h",
    "L",
    "L",
    "fraction",
    "fraction",
    "days",
    "L/h",
    "L"
  )
)

model_information_content <- function() {
  tagList(
    tags$p(
      "This simulator implements the final tirzepatide population PK model from Schneck and Urva. The manuscript pooled 39,644 concentration records from 5,802 participants across 19 studies and concluded that tirzepatide PK is well described by a two-compartment model with first-order SC absorption and elimination."
    ),
    tags$h5("Model structure"),
    tags$ul(
      tags$li(tags$strong("Structure:"), " SC depot -> central -> peripheral, with first-order absorption and linear elimination"),
      tags$li(tags$strong("Body-size scaling:"), " CL and Q scale with total body weight to the 0.8 power; Vc and Vp scale with fat-free mass plus 48.2% of fat mass"),
      tags$li(tags$strong("Bioavailability:"), " absolute F fixed at 0.8, with a pooled-study relative factor available from the paper"),
      tags$li(tags$strong("Clinical conclusion:"), " no dose adjustment required for age, renal function, hepatic markers, sex, race, or ethnicity after accounting for body weight")
    ),
    tags$h5("Published parameters from Table 3"),
    tags$table(
      class = "table table-sm",
      tags$thead(tags$tr(tags$th("Parameter"), tags$th("Estimate"), tags$th("Units"))),
      tags$tbody(
        lapply(seq_len(nrow(paper_parameters)), function(i) {
          tags$tr(
            tags$td(paper_parameters$parameter[[i]]),
            tags$td(paper_parameters$estimate[[i]]),
            tags$td(paper_parameters$units[[i]])
          )
        })
      )
    ),
    tags$p("The paper reports dose-proportional exposure from 0.25 to 15 mg, Tmax between 8 and 72 hours after SC dosing, mean accumulation of 1.7-fold, and a mean half-life of 5.4 days."),
    tags$p("The published +22%/-33% exposure summary across 70-120 kg came from virtual-patient resampling with correlated demographics. Because this simulator lets weight, BMI, and sex be adjusted independently, isolated slider moves will not necessarily reproduce that summary exactly.")
  )
}

references_content <- function() {
  tagList(
    tags$h5("Primary reference"),
    tags$ol(
      tags$li(
        "Schneck KB, Urva S. Population pharmacokinetics of the GIP/GLP receptor agonist tirzepatide. ",
        tags$em("CPT: Pharmacometrics & Systems Pharmacology"),
        ". 2024;13(4):645-657. ",
        tags$a(href = "https://pubmed.ncbi.nlm.nih.gov/38356317/", target = "_blank", "PMID: 38356317")
      )
    ),
    tags$h5("Drug summary"),
    tags$ul(
      tags$li(tags$strong("Generic:"), " Tirzepatide"),
      tags$li(tags$strong("Brands:"), " Mounjaro, Zepbound"),
      tags$li(tags$strong("Class:"), " Dual GIP/GLP-1 receptor agonist"),
      tags$li(tags$strong("Route:"), " Subcutaneous"),
      tags$li(tags$strong("Approved maintenance doses:"), " 5, 10, or 15 mg once weekly after 2.5 mg initiation and 2.5 mg step-up titration"),
      tags$li(tags$strong("PK takeaway:"), " Exposure varies mainly with body size, but the observed differences were not considered clinically meaningful enough to require dose adjustment")
    ),
    tags$hr(),
    tags$p(
      style = "font-size: 12px; color: #6c757d; margin-bottom: 0;",
      "For research and educational purposes only. Not for clinical decision-making. Powered by ",
      tags$a(href = "https://www.pkpdbuilder.com", target = "_blank", "PKPDBuilder.com"),
      "."
    )
  )
}

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; }
    .summary-box {
      background: #f8fafc;
      border: 1px solid #e2e8f0;
      border-radius: 8px;
      padding: 14px 16px;
      margin-top: 12px;
      font-size: 13px;
      line-height: 1.5;
    }
    .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; }
    .regimen-note {
      background: #eef6ff;
      border-left: 3px solid #2563eb;
      border-radius: 4px;
      color: #334155;
      font-size: 12px;
      line-height: 1.45;
      padding: 9px 11px;
      margin-top: 8px;
    }
  ")

ui <- page_sidebar(
  title = "Zepbound (Tirzepatide) PK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    accordion(
      open = c("Dosing", "Patient"),
      accordion_panel(
        "Dosing",
        radioButtons(
          "regimen_mode",
          "Regimen",
          choices = c(
            "Zepbound titration" = "titration",
            "Fixed weekly dose" = "fixed",
            "Custom week-by-week" = "custom"
          ),
          selected = "titration"
        ),
        conditionalPanel(
          "input.regimen_mode == 'titration'",
          selectInput(
            "maintenance_dose",
            "Target maintenance dose",
            choices = c(5, 10, 15),
            selected = 10
          ),
          sliderInput("weeks_per_step", "Weeks at each dose", min = 4, max = 12, value = 4, step = 1),
          div(class = "regimen-note", "Starts at 2.5 mg weekly and increases by 2.5 mg after the selected number of weeks until the target dose is reached.")
        ),
        conditionalPanel(
          "input.regimen_mode == 'fixed'",
          sliderInput("dose_mg", "Weekly dose (mg)", min = 2.5, max = 15, value = 5, step = 2.5)
        ),
        conditionalPanel(
          "input.regimen_mode == 'custom'",
          textAreaInput(
            "custom_doses",
            "Dose for each week (mg)",
            value = "2.5, 2.5, 2.5, 2.5, 5, 5, 5, 5, 7.5, 7.5, 7.5, 7.5, 10, 10, 10, 10",
            rows = 5,
            resize = "vertical"
          ),
          div(class = "regimen-note", "Separate doses with commas. Use 0 for a skipped or held dose. The simulation duration follows the number of entries.")
        ),
        conditionalPanel(
          "input.regimen_mode != 'custom'",
          sliderInput("n_weeks", "Simulation duration (weeks)", min = 4, max = 52, value = 24, step = 1)
        )
      ),
      accordion_panel(
        "Patient",
        sliderInput("wt", "Body weight (kg)", min = 50, max = 160, value = 90, step = 1),
        sliderInput("bmi", "Body mass index (kg/m²)", min = 20, max = 50, value = 32, step = 0.5),
        radioButtons("sex", "Sex used for FFM equation", choices = c("Female", "Male"), selected = "Female")
      ),
      accordion_panel(
        "Model",
        selectInput(
          "bioavailability",
          "Bioavailability scenario",
          choices = c(
            "Pooled phase II/III studies (default)" = "pooled",
            "Reference absolute-bioavailability study" = "reference"
          ),
          selected = "pooled"
        ),
        tags$p(
          style = "font-size: 12px; color: #6c757d; margin-top: 8px;",
          "The pooled-study setting applies the manuscript's relative bioavailability term and reproduces the reported mean apparent clearance and apparent volume more closely."
        )
      ),
      accordion_panel(
        "Display",
        checkboxInput("log_scale", "Log scale (Y-axis)", value = FALSE),
        actionButton("show_model_info", "Model information", class = "btn-outline-primary w-100 mb-2"),
        actionButton("show_references", "References & disclaimer", class = "btn-outline-secondary w-100")
      )
    )
  ),
  tags$p(
    style = "font-size: 11px; color: #6c757d; margin: 4px 0 0 6px;",
    "Two-compartment SC population PK model from Schneck and Urva (CPT:PSP 2024) with semimechanistic body-composition scaling."
  ),
  layout_column_wrap(
    width = 1 / 4,
    fill = FALSE,
    metric_card("cmax", "Final-week Cmax (ng/mL)", "metric-success"),
    metric_card("ctrough", "Final-week Ctrough (ng/mL)", "metric-warning"),
    metric_card("auc", "Final-week AUC (ng·h/mL)", "metric-primary"),
    metric_card("thalf", "t½ (days)", "metric-info")
  ),
  card(
    card_header("Concentration-Time Profile"),
    full_screen = TRUE,
    div(
      style = "height: clamp(420px, 58vh, 680px); min-height: 420px;",
      plotOutput("pkPlot", height = "100%")
    ),
    div(
      class = "summary-box",
      strong("Current simulation summary"),
      br(),
      textOutput("simSummary")
    )
  )
)

server <- function(input, output, session) {
  observeEvent(input$show_model_info, {
    showModal(modalDialog(
      title = "Model information",
      model_information_content(),
      size = "l",
      easyClose = TRUE,
      footer = modalButton("Close")
    ))
  })

  observeEvent(input$show_references, {
    showModal(modalDialog(
      title = "References & disclaimer",
      references_content(),
      size = "l",
      easyClose = TRUE,
      footer = modalButton("Close")
    ))
  })

  dose_schedule <- reactive({
    shiny::req(input$regimen_mode)
    build_weekly_schedule(
      mode = input$regimen_mode,
      n_weeks = input$n_weeks %||% 24,
      fixed_dose = input$dose_mg %||% 5,
      maintenance_dose = as.numeric(input$maintenance_dose %||% 10),
      weeks_per_step = input$weeks_per_step %||% 4,
      custom_doses = input$custom_doses %||% ""
    )
  })

  sim_data <- reactive({
    shiny::req(input$wt, input$bmi)
    simulate_weekly_schedule(
      doses = dose_schedule(),
      wt = input$wt,
      bmi = input$bmi,
      sex = input$sex,
      bioavailability = input$bioavailability
    )
  })

  sim_metrics <- reactive({
    last_interval_metrics(sim_data())
  })

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

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

  output$auc <- renderText({
    sprintf("%.0f", sim_metrics()$auc_tau)
  })

  output$thalf <- renderText({
    sprintf("%.1f", sim_metrics()$thalf_days)
  })

  output$simSummary <- renderText({
    metrics <- sim_metrics()
    doses <- dose_schedule()
    sprintf(
      "%d-week regimen; final dose %.1f mg; Tmax in final week %.1f h; apparent CL/F %.3f L/h; apparent Vss/F %.1f L; F %.3f.",
      length(doses),
      tail(doses, 1),
      metrics$tmax_h,
      metrics$clf,
      metrics$vssf,
      metrics$f1
    )
  })

  output$pkPlot <- renderPlot({
    plot_data <- sim_data()
    doses <- dose_schedule()
    dose_data <- data.frame(
      time_days = (seq_along(doses) - 1) * 7,
      dose = doses
    ) |>
      dplyr::filter(dose > 0)
    if (isTRUE(input$log_scale)) {
      plot_data <- plot_data |>
        dplyr::filter(CP_ngml > 0)
    }

    p <- ggplot(plot_data, aes(x = time_days, y = CP_ngml)) +
      geom_line(color = "#2563eb", linewidth = 0.9) +
      geom_rug(data = dose_data, aes(x = time_days), inherit.aes = FALSE,
               sides = "b", color = "#dc2626", linewidth = 0.6) +
      labs(
        x = "Time (days)",
        y = "Plasma concentration (ng/mL)",
        title = "Zepbound concentration across the selected regimen",
        subtitle = "Red ticks mark administered doses"
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank()
      )

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

    p
  })
}

shinyApp(ui = ui, server = server)