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

profiles <- tibble::tribble(
  ~study,  ~age_group, ~CL,  ~F1,  ~KA,  ~KA_oral_fit, ~Q_fit, ~V1,   ~V2,   ~Vss,  ~thalf_ref, ~auc_iv_ref, ~cmax_iv_ref, ~auc_oral_ref, ~cmax_oral_ref, ~bsa_iv_ref, ~bsa_oral_ref,
  "renal", "0-5",      4.71, 0.44, 0.80, 0.80,         2.00,   7.20,  5.71,  13.19, 3.28,       22.18,       10.19,        22.22,         5.10,           0.52,        0.46,
  "renal", "6-11",     4.92, 0.54, 0.67, 0.80,         7.00,   15.03, 9.94,  24.99, 4.41,       37.86,       9.03,         43.78,         6.01,           0.93,        0.77,
  "renal", "12-16",    7.40, 0.53, 0.72, 0.75,         10.00,  22.07, 18.84, 40.81, 5.62,       38.58,       9.40,         39.88,         5.40,           1.43,        1.07,
  "liver", "0-5",      4.05, 0.52, 0.45, 0.80,         8.75,   1.66,  5.65,  7.62,  1.65,       24.30,       12.20,        23.40,         5.51,           0.49,        0.35,
  "liver", "6-11",     2.86, 0.71, 0.35, 0.55,         6.75,   5.74,  14.60, 20.40, 6.80,       35.20,       9.29,         46.80,         5.29,           0.50,        0.36,
  "liver", "12-16",    15.1, 0.64, 0.42, 0.85,         12.00,  12.80, 30.70, 43.50, 4.35,       23.40,       11.80,        25.80,         6.90,           1.77,        1.17
)

calc_bsa <- function(height_cm, weight_kg) {
  sqrt((height_cm * weight_kg) / 3600)
}

renal_fraction <- function(crcl) {
  if (crcl >= 70) {
    1
  } else if (crcl >= 50) {
    0.5
  } else if (crcl >= 40) {
    0.25
  } else {
    0
  }
}

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

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

model_code <- '
$PARAM @annotated
CLi   : 5.40 : Clearance (L/h)
V1i   : 5.00 : Central volume (L)
Qi    : 8.00 : Intercompartmental clearance (L/h)
V2i   : 15.0 : Peripheral volume (L)
KAi   : 0.42 : First-order absorption rate constant (1/h)
F1i   : 0.55 : Relative bioavailability of ganciclovir from valganciclovir
TLAGi : 0.22 : Absorption lag time (h)

$CMT @annotated
DEPOT  : Oral valganciclovir depot (mg)
CENT   : Central compartment (mg)
PERIPH : Peripheral compartment (mg)

$MAIN
F_DEPOT = F1i;
ALAG_DEPOT = TLAGi;

$ODE
dxdt_DEPOT  = -KAi * DEPOT;
dxdt_CENT   = KAi * DEPOT
              - (CLi / V1i) * CENT
              - (Qi / V1i) * CENT
              + (Qi / V2i) * PERIPH;
dxdt_PERIPH = (Qi / V1i) * CENT - (Qi / V2i) * PERIPH;

$TABLE
double CP = CENT / V1i;

$CAPTURE CP CLi V1i Qi V2i KAi F1i TLAGi
'

mod <- mcode("pescovitz_pediatric_ganciclovir", model_code, quiet = TRUE)

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; }
  .dose-box {
    background: #eef6ff;
    border-left: 4px solid #0d6efd;
    padding: 12px 16px;
    border-radius: 4px;
    font-size: 13px;
  }
  .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 <- page_sidebar(
  title = "Pediatric Ganciclovir / Valganciclovir PK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 360,

    h6("Clinical Scenario"),
    selectInput(
      "study",
      "Study Setting",
      choices = c(
        "Renal transplant study" = "renal",
        "Liver transplant study" = "liver",
        "Custom schedule" = "custom"
      ),
      selected = "renal"
    ),
    selectInput(
      "age_group",
      "Age Band (years)",
      choices = c("0-5", "6-11", "12-16"),
      selected = "6-11"
    ),

    hr(),
    h6("Route and Dosing"),
    radioButtons(
      "route",
      NULL,
      choices = c("IV Ganciclovir" = "iv", "Oral Valganciclovir" = "oral"),
      selected = "oral",
      inline = TRUE
    ),
    checkboxInput(
      "use_algorithm",
      "Use valganciclovir algorithm: 7 x BSA x CrCL (cap 900 mg)",
      value = TRUE
    ),
    uiOutput("dose_input"),
    checkboxInput(
      "apply_renal_adjustment",
      "Apply study renal dose reduction bands (CrCL >=70/50-69/40-49/<40)",
      value = TRUE
    ),
    selectInput(
      "interval",
      "Dosing Interval",
      choices = c("Every 12h (BID)" = 12, "Every 24h (QD)" = 24),
      selected = 24
    ),
    numericInput("n_days", "Duration (days)", value = 4, min = 1, max = 21),

    hr(),
    h6("Patient Characteristics"),
    sliderInput("height", "Height (cm)", min = 55, max = 190, value = 125, step = 1),
    sliderInput("weight", "Weight (kg)", min = 5, max = 100, value = 25, step = 0.5),
    sliderInput(
      "crcl",
      "Estimated CrCL (mL/min/1.73 m²)",
      min = 20,
      max = 220,
      value = 110,
      step = 5
    ),

    hr(),
    uiOutput("dose_note"),
    hr(),
    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_out")),
      div(class = "metric-label", "Cmax (μg/mL)")
    ),
    div(
      class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("ctrough_out")),
      div(class = "metric-label", "Ctrough (μg/mL)")
    ),
    div(
      class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("auc_out")),
      div(class = "metric-label", "AUC0-24 (μg·h/mL)")
    ),
    div(
      class = "metric-card metric-info",
      div(class = "metric-value", textOutput("thalf_out")),
      div(class = "metric-label", "t½ terminal (h)")
    )
  ),

  navset_card_underline(
    nav_panel(
      "PK Simulation",
      card(full_screen = TRUE, plotOutput("pkPlot", height = "500px"))
    ),
    nav_panel(
      "Model Information",
      markdown("
## Pescovitz et al. 2010 Pediatric Transplant PopPK Simulator

**Publication:** Pescovitz MD, Ettenger RB, Strife CF, et al. *Transpl Infect Dis.* 2010;12:195-203.  
**Population:** Pediatric renal and liver transplant recipients receiving IV ganciclovir or oral valganciclovir solution.  
**Structural model:** Two-compartment ganciclovir model with first-order formation after oral valganciclovir, lag time, and relative bioavailability.

### Final Combined Population Parameters (Table 3)

| Parameter | Estimate | Unit | Notes |
|-----------|----------|------|-------|
| Ka | 0.42 | h⁻¹ | Combined oral absorption / formation |
| F1 | 0.55 | fraction | Ganciclovir from oral valganciclovir |
| Tlag | 0.22 | h | Oral lag time |
| CL | 5.4 | L/h | Final combined model |
| Vss | 20 | L | Final combined model |
| Vperiph | 15 | L | Final combined model |
| Q | 8 | L/h | Final combined model |

### Paper-Anchored Pediatric Profiles (Table 5)

This simulator uses the published **age-band and study-specific post hoc medians** for CL, F1, Ka, Vcent, and Vperiph.  
That choice is deliberate: the paper states that **height** and **CrCL** were significant covariates, but it does **not publish the exact covariate coefficients** needed to recreate the original NONMEM covariate equations.  
So the app keeps the published 2-CMT structure and uses the paper's pediatric parameter profiles directly for each age band.

To better align the simulated peak concentrations with **Table 4**, the app uses lightly calibrated distribution and oral absorption terms for each age band while preserving the paper-anchored exposure values.

### Dosing Logic

- **Reference IV dose:** 200 mg/m²
- **Reference oral valganciclovir dose:** 520 mg/m²
- **Renal adjustment bands (Table 1):**
  - CrCL ≥ 70: 100% of full dose
  - CrCL 50-69: 50%
  - CrCL 40-49: 25%
  - CrCL < 40: withhold reference dose
- **Later pediatric oral algorithm discussed by the authors:** dose = **7 x BSA x CrCL**, capped at **900 mg**

### Key Findings from the Paper

- Oral valganciclovir 520 mg/m² produced ganciclovir exposure broadly similar to IV ganciclovir 200 mg/m².
- Younger children were underexposed when dosing used **BSA alone**.
- Adding **renal function (CrCL)** improved exposure matching across ages.
- Type of transplant was **not a significant structural PK covariate** in the final model.
- In the liver cohort, sampling was performed over 12 hours around the morning dose, but the manuscript reports **AUC0-24** using a modeled once-daily convention. The app defaults the liver comparison view to that published 24-hour convention.
      ")
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Primary Reference"),
        tags$p(
          "Pescovitz MD, Ettenger RB, Strife CF, Sherbotie JR, Thomas SE, McDiarmid S, ",
          "Bartosh S, Ives J, Bouw MR, Bucuvalas J. Pharmacokinetics of oral valganciclovir ",
          "solution and intravenous ganciclovir in pediatric renal and liver transplant recipients. ",
          tags$em("Transpl Infect Dis."),
          " 2010;12:195-203. ",
          tags$a(
            href = "https://doi.org/10.1111/j.1399-3062.2009.00478.x",
            target = "_blank",
            "DOI: 10.1111/j.1399-3062.2009.00478.x"
          )
        ),
        tags$h5("Drug Context"),
        tags$ul(
          tags$li(tags$strong("Drug:"), " Ganciclovir / oral prodrug valganciclovir"),
          tags$li(tags$strong("Class:"), " Antiviral nucleoside analogue"),
          tags$li(tags$strong("Clinical use:"), " CMV prophylaxis / treatment in transplant recipients"),
          tags$li(tags$strong("Elimination:"), " Predominantly renal excretion of unchanged ganciclovir"),
          tags$li(tags$strong("Study focus:"), " Pediatric BSA and renal function informed dosing")
        ),
        tags$h5("Important Limitation"),
        tags$p(
          "The original paper reports the significant covariates but does not publish the full covariate equations. ",
          "This app therefore anchors simulation to the paper's published age-band parameter medians and dosing logic, ",
          "rather than claiming an exact reconstruction of the authors' unpublished NONMEM covariate coefficients."
        )
      )
    )
  ),

  div(
    style = "text-align:center; padding:20px; margin-top:30px; border-top:1px solid #e9ecef; color:#6c757d; font-size:12px;",
    "Powered by ",
    tags$a(
      href = "https://www.pkpdbuilder.com",
      target = "_blank",
      style = "color:#8b5cf6; font-weight:500;",
      "PKPDBuilder.com"
    ),
    " • Pediatric transplant PopPK • Pescovitz et al. 2010",
    tags$br(),
    tags$span(
      style = "font-size:10px;",
      "For research and educational purposes only. Not for clinical decision-making."
    )
  )
)

server <- function(input, output, session) {
  observe({
    if (input$route == "iv") {
      updateCheckboxInput(session, "use_algorithm", value = FALSE)
    }
  })

  observeEvent(input$route, {
    if (input$route == "iv") {
      updateNumericInput(session, "dose_m2", value = 200)
    } else {
      updateNumericInput(session, "dose_m2", value = 520)
    }
  }, ignoreInit = TRUE)

  observeEvent(input$study, {
    if (input$study == "renal") {
      updateSelectInput(session, "interval", selected = "24")
      updateNumericInput(session, "n_days", value = 4)
    } else if (input$study == "liver") {
      updateSelectInput(session, "interval", selected = "24")
      updateNumericInput(session, "n_days", value = 14)
    }
  }, ignoreInit = FALSE)

  observeEvent(input$age_group, {
    if (input$age_group == "0-5") {
      updateSliderInput(session, "height", value = 90)
      updateSliderInput(session, "weight", value = 13)
      updateSliderInput(session, "crcl", value = 110)
    } else if (input$age_group == "6-11") {
      updateSliderInput(session, "height", value = 125)
      updateSliderInput(session, "weight", value = 25)
      updateSliderInput(session, "crcl", value = 110)
    } else {
      updateSliderInput(session, "height", value = 160)
      updateSliderInput(session, "weight", value = 50)
      updateSliderInput(session, "crcl", value = 110)
    }
  }, ignoreInit = FALSE)

  output$dose_input <- renderUI({
    shiny::req(input$route)
    if (input$route == "oral" && isTRUE(input$use_algorithm)) {
      return(NULL)
    }

    label <- if (input$route == "iv") "Dose (mg/m²)" else "Valganciclovir Dose (mg/m²)"
    numericInput("dose_m2", label, value = if (input$route == "iv") 200 else 520, min = 50, max = 1000, step = 10)
  })

  selected_profile <- reactive({
    study_key <- if (input$study %in% c("renal", "liver")) input$study else "renal"
    dplyr::filter(profiles, study == study_key, age_group == input$age_group)
  })

  dose_info <- reactive({
    shiny::req(input$height, input$weight, input$crcl, input$route)
    p <- selected_profile()
    bsa <- calc_bsa(input$height, input$weight)
    fraction <- if (isTRUE(input$apply_renal_adjustment)) renal_fraction(input$crcl) else 1

    if (input$route == "oral" && isTRUE(input$use_algorithm)) {
      dose_mg <- min(900, 7 * bsa * input$crcl)
      regimen <- "7 x BSA x CrCL oral algorithm"
      applied_bsa <- bsa
    } else {
      shiny::req(input$dose_m2)
      applied_bsa <- if (input$route == "iv") p$bsa_iv_ref else p$bsa_oral_ref
      dose_mg <- input$dose_m2 * applied_bsa * fraction
      regimen <- paste0(input$dose_m2, " mg/m² reference-style dose (paper-median BSA)")
    }

    list(
      current_bsa = bsa,
      applied_bsa = applied_bsa,
      fraction = fraction,
      dose_mg = dose_mg,
      regimen = regimen
    )
  })

  output$dose_note <- renderUI({
    d <- dose_info()
    p <- selected_profile()
    auc_ref <- if (input$route == "iv") p$auc_iv_ref else p$auc_oral_ref
    cmax_ref <- if (input$route == "iv") p$cmax_iv_ref else p$cmax_oral_ref

    div(
      class = "dose-box",
      tags$strong("Calculated regimen"),
      tags$br(),
      sprintf(
        "Current BSA %.2f m² | Applied BSA %.2f m² | Actual dose %.0f mg | Renal fraction %.2f",
        d$current_bsa, d$applied_bsa, d$dose_mg, d$fraction
      ),
      tags$br(),
      d$regimen,
      tags$hr(),
      tags$strong("Paper reference for selected profile"),
      tags$br(),
      sprintf(
        "AUC %.2f μg·h/mL | Cmax %.2f μg/mL | Reference t½ %.2f h",
        auc_ref, cmax_ref, p$thalf_ref
      )
    )
  })

  sim_data <- reactive({
    p <- selected_profile()
    d <- dose_info()
    tau <- as.numeric(input$interval)
    shiny::req(nrow(p) == 1, tau > 0, input$n_days >= 1)

    q_est <- p$Q_fit
    ka_est <- if (input$route == "oral") p$KA_oral_fit else p$KA
    tlag_est <- if (input$route == "oral") 0 else 0.22
    sim_mod <- mod %>%
      param(
        CLi = p$CL,
        V1i = p$V1,
        Qi = q_est,
        V2i = p$V2,
        KAi = ka_est,
        F1i = p$F1,
        TLAGi = tlag_est
      )

    n_doses <- max(1, ceiling((input$n_days * 24) / tau))
    if (input$route == "oral") {
      ev_obj <- ev(amt = d$dose_mg, cmt = 1, ii = tau, addl = n_doses - 1)
    } else {
      rate_val <- if (d$dose_mg > 0) d$dose_mg / 1 else 0
      ev_obj <- ev(amt = d$dose_mg, cmt = 2, ii = tau, addl = n_doses - 1, rate = rate_val)
    }

    sim_mod %>%
      ev(ev_obj) %>%
      mrgsim(end = input$n_days * 24, delta = 0.1) %>%
      as.data.frame()
  })

  last_interval <- reactive({
    d <- sim_data()
    tau <- as.numeric(input$interval)
    shiny::req(nrow(d) > 0)
    end_time <- max(d$time, na.rm = TRUE)
    start_time <- max(0, end_time - tau)
    dplyr::filter(d, time >= start_time, time <= end_time)
  })

  last_24h <- reactive({
    d <- sim_data()
    shiny::req(nrow(d) > 0)
    end_time <- max(d$time, na.rm = TRUE)
    start_time <- max(0, end_time - 24)
    dplyr::filter(d, time >= start_time, time <= end_time)
  })

  terminal_half_life <- reactive({
    p <- selected_profile()
    q_est <- p$Q_fit
    k10 <- p$CL / p$V1
    k12 <- q_est / p$V1
    k21 <- q_est / p$V2
    disc <- (k10 + k12 + k21)^2 - 4 * k10 * k21
    beta <- 0.5 * ((k10 + k12 + k21) - sqrt(max(disc, 0)))
    if (!is.finite(beta) || beta <= 0) {
      return(NA_real_)
    }
    log(2) / beta
  })

  output$cmax_out <- renderText({
    d <- last_24h()
    shiny::req(nrow(d) > 0)
    sprintf("%.2f", max(d$CP, na.rm = TRUE))
  })

  output$ctrough_out <- renderText({
    d <- last_interval()
    shiny::req(nrow(d) > 0)
    sprintf("%.2f", tail(d$CP, 1))
  })

  output$auc_out <- renderText({
    d <- last_24h()
    shiny::req(nrow(d) > 1)
    sprintf("%.2f", auc_trap(d$time, d$CP))
  })

  output$thalf_out <- renderText({
    t12 <- terminal_half_life()
    if (!is.finite(t12)) {
      return("NA")
    }
    sprintf("%.2f", t12)
  })

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

    plt <- d |>
      dplyr::mutate(time_days = time / 24) |>
      ggplot(aes(x = time_days, y = CP)) +
      geom_line(color = "#8b5cf6", linewidth = 1.1) +
      labs(
        x = "Time (days)",
        y = "Ganciclovir Concentration (μg/mL)",
        title = "Predicted Ganciclovir Plasma Concentration-Time Profile",
        subtitle = paste(
          if (input$route == "iv") "IV ganciclovir" else "Oral valganciclovir",
          "|",
          toupper(if (input$study == "liver") "liver study profile" else "renal study profile"),
          "| Age", input$age_group
        )
      ) +
      theme_minimal(base_size = 13) +
      theme(
        plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank()
      )

    if (isTRUE(input$log_scale) && any(d$CP > 0, na.rm = TRUE)) {
      plt <- d |>
        dplyr::filter(CP > 0) |>
        dplyr::mutate(time_days = time / 24) |>
        ggplot(aes(x = time_days, y = CP)) +
        geom_line(color = "#8b5cf6", linewidth = 1.1) +
        scale_y_log10() +
        labs(
          x = "Time (days)",
          y = "Ganciclovir Concentration (μg/mL)",
          title = "Predicted Ganciclovir Plasma Concentration-Time Profile",
          subtitle = paste(
            if (input$route == "iv") "IV ganciclovir" else "Oral valganciclovir",
            "|",
            toupper(if (input$study == "liver") "liver study profile" else "renal study profile"),
            "| Age", input$age_group
          )
        ) +
        theme_minimal(base_size = 13) +
        theme(
          plot.title = element_text(face = "bold"),
          panel.grid.minor = element_blank()
        )
    }

    plt
  })
}

shinyApp(ui, server)