# Ganciclovir PK/PD Simulator - CMV Viral Load
# Based on: Koloskoff et al., Clin Pharmacokinet (2025)
# DOI: 10.1007/s40262-025-01526-z

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

# ── PK/PD Model Code ──
# Indirect viral turnover model with stimulated degradation
model_code <- '
$PARAM @annotated
kin    : 0.00087   : Virus production rate (log10 copies/mL/h)
kout   : 0.00023   : Virus elimination rate (1/h)
Emax   : 16.3      : Maximum drug effect (unitless)
EC50   : 23.5      : AUC for half-max effect (mg*h/L)
VL0    : 3.61      : Baseline viral load (log10 copies/mL)
AUC24  : 40        : Ganciclovir AUC0-24 (mg*h/L)

$CMT @annotated
VL : CMV Viral Load (log10 copies/mL)

$MAIN
// Initial viral load - set to baseline on first dose
VL_0 = VL0;

// AUC0-12 for the model (model uses 12h AUC internally)
double AUC12 = AUC24 / 2.0;

// Drug effect on viral elimination
double EFFECT = (Emax * AUC12) / (EC50 + AUC12);

$ODE
// Indirect response model: stimulation of viral elimination
// dR/dt = kin - kout * (1 + E) * R
dxdt_VL = kin - kout * (1.0 + EFFECT) * VL;

$TABLE
// Viral load output (constrained to LLOQ)
double VL_out = VL;
double VL_log = VL > 2.0 ? VL : 2.0;  // LLOQ = 2 log10 copies/mL (200 copies/mL)
double VL_detectable = VL > 2.0 ? 1.0 : 0.0;

$CAPTURE @annotated
VL_out       : Viral load (log10 copies/mL)
VL_log       : Viral load clamped at LLOQ (log10 copies/mL)
VL_detectable : Detectable flag (1 = yes)
'

# Compile model
mod <- mcode("ganciclovir_cmv", model_code, compile = TRUE)

# ── Theme with Custom Metric Cards ──
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; }
  .metric-danger .metric-value { color: #ef4444; }
  
  .target-zone {
    background: linear-gradient(135deg, #10b98110 0%, #8b5cf610 100%);
    border: 1px dashed #10b981;
    border-radius: 8px;
    padding: 12px;
    margin-top: 10px;
  }
  .target-title {
    font-weight: 600;
    color: #10b981;
    margin-bottom: 5px;
  }
")

# ── UI ──
ui <- page_sidebar(
  title = "Ganciclovir PK/PD Simulator: CMV Viral Load",
  theme = app_theme,
  
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    
    h6("Drug Exposure"),
    sliderInput("auc24", "Ganciclovir AUC₀₋₂₄ (mg·h/L)", 
                min = 10, max = 100, value = 50, step = 5),
    helpText("Target range: 40-60 mg·h/L"),
    
    hr(),
    h6("Baseline Infection"),
    sliderInput("vl0", "Initial Viral Load (log₁₀ copies/mL)", 
                min = 2.5, max = 6, value = 3.6, step = 0.1),
    helpText("Typical range: 2.5-5.0"),
    
    hr(),
    h6("Simulation Duration"),
    sliderInput("duration", "Treatment Duration (days)", 
                min = 7, max = 56, value = 28, step = 7),
    
    hr(),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE),
    checkboxInput("show_lloq", "Show LLOQ Threshold", value = TRUE),
    
    hr(),
    # Target reference card
    div(class = "target-zone",
      div(class = "target-title", "Therapeutic Targets"),
      p(style = "margin: 0; font-size: 12px;",
        strong("AUC₀₋₂₄:"), " 40-60 mg·h/L", br(),
        strong("Goal:"), " Undetectable VL (<200 copies/mL)", br(),
        strong("≥40 mg·h/L:"), " 85% clear by day 28", br(),
        strong(">60 mg·h/L:"), " Minimal added benefit"
      )
    ),
    
    hr(),
    card(
      card_header("Model Parameters", class = "bg-light"),
      p(style = "font-size: 11px; margin: 0;",
        strong("k", tags$sub("in"), ":"), " 0.00087 (virus production)", br(),
        strong("k", tags$sub("out"), ":"), " 0.00023 (virus elimination)", br(),
        strong("E", tags$sub("max"), ":"), " 16.3", br(),
        strong("EC", tags$sub("50"), ":"), " 23.5 mg·h/L", br(),
        strong("LLOQ:"), " 2.0 log₁₀ (200 copies/mL)"
      )
    )
  ),
  
  # ── Main Content ──
  
  # Metrics row
  layout_column_wrap(
        width = 1/4,
        fill = FALSE,
    div(class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("auc_display")),
      div(class = "metric-label", "AUC₀₋₂₄ (mg·h/L)")
    ),
    div(class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("vl_day14")),
      div(class = "metric-label", "VL at Day 14")
    ),
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("vl_day28")),
      div(class = "metric-label", "VL at Day 28")
    ),
    div(class = "metric-card metric-info",
      div(class = "metric-value", textOutput("time_lloq")),
      div(class = "metric-label", "Days to LLOQ")
    )
  ),
  
  # Main PD plot
  card(
    card_header("CMV Viral Load Decline"),
    full_screen = TRUE,
    plotOutput("pdPlot", height = "420px")
  ),
  
  # Tabbed content
  navset_card_tab(
    title = "Analysis",
    
    nav_panel("AUC Comparison",
      plotOutput("comparePlot", height = "420px")
    ),
    
    nav_panel("Probability of Response",
      layout_columns(
        col_widths = c(6, 6),
        card(
          card_header("Day 14: -1 log₁₀ Reduction"),
          tableOutput("probTable14")
        ),
        card(
          card_header("Day 28: Undetectable VL"),
          tableOutput("probTable28")
        )
      )
    ),
    
    nav_panel("Model Information",
      markdown("
## Pharmacodynamic Model

This simulator implements the population PK/PD model from Koloskoff et al. (2025) 
for ganciclovir/valganciclovir treatment of CMV infection in pediatric transplant recipients.

### Model Structure
**Indirect Response Model with Stimulated Viral Degradation:**

$$\\frac{dR}{dt} = k_{in} - k_{out} \\times \\left(1 + \\frac{E_{max} \\times AUC_{0-12}}{EC_{50} + AUC_{0-12}}\\right) \\times R$$

Where:
- **R** = CMV viral load (log₁₀ copies/mL)
- **k_in** = Virus production rate (0.00087 log₁₀ copies/mL/h)
- **k_out** = Virus elimination rate (0.00023 h⁻¹)
- **E_max** = Maximum drug effect (16.3)
- **EC₅₀** = AUC for half-maximal effect (23.5 mg·h/L)

### Clinical Targets
- **AUC₀₋₂₄ = 40-60 mg·h/L**: Optimal therapeutic window
- **AUC₀₋₂₄ ≥ 40 mg·h/L**: 85.4% probability of undetectable viral load at day 28
- **AUC₀₋₂₄ > 60 mg·h/L**: Minimal additional efficacy; may increase toxicity risk

### Reference
Koloskoff K, et al. *Pharmacokinetic/Pharmacodynamic Modelling and Monte Carlo Simulations 
to Predict Cytomegalovirus Viral Load in Pediatric Transplant Recipients Treated with 
(val)Ganciclovir.* Clin Pharmacokinet. 2025. DOI: [10.1007/s40262-025-01526-z](https://doi.org/10.1007/s40262-025-01526-z)

### Indication
- CMV infection/reactivation in SOT (solid organ transplant) or HSCT (hematopoietic stem cell transplant) recipients
- Preemptive therapy upon detection of CMV DNAemia

### Limitations
- Model based on pediatric data (n=36 occurrences in 29 children)
- Does not include toxicity endpoints (myelosuppression)
- Does not account for resistant CMV strains
      ")
    ),
    
    nav_panel("References",
      markdown("
## FDA Label
- [Ganciclovir Injection Label (FDA)](https://www.accessdata.fda.gov/drugsatfda_docs/label/2017/209347lbl.pdf)

## Key Publications
1. **Primary Source:** Koloskoff K, et al. Clin Pharmacokinet 2025. DOI: 10.1007/s40262-025-01526-z
2. Franck B, et al. Population pharmacokinetics of ganciclovir and valganciclovir in paediatric transplant recipients. Br J Clin Pharmacol. 2021;87:3105-14.
3. Stockmann C, et al. Clinical pharmacokinetics and pharmacodynamics of ganciclovir and valganciclovir in children with CMV infection. Expert Opin Drug Metab Toxicol. 2015;11:205-19.
4. Kotton CN, et al. Third International Consensus Guidelines on CMV Management in Solid-organ Transplantation. Transplantation. 2018;102:900-31.

## Therapeutic Context
- **Ganciclovir**: Nucleoside analogue; inhibits CMV DNA polymerase
- **Valganciclovir**: Oral prodrug with ~60% bioavailability
- **Primary indication**: CMV prophylaxis and treatment in transplant recipients
- **Monitoring**: TDM recommended due to high PK variability; target AUC₀₋₂₄ = 40-60 mg·h/L
      ")
    )
  ),

  # ── PKPDBuilder Branding Footer ──
  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"),
    " • 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 ──
server <- function(input, output, session) {
  
  # Run simulation
  sim_data <- reactive({
    shiny::req(input$auc24, input$vl0, input$duration)
    
    # Update model parameters
    sim <- mod %>%
      param(AUC24 = input$auc24, VL0 = input$vl0) %>%
      mrgsim(end = input$duration * 24, delta = 1) %>%  # hourly resolution
      as.data.frame() %>%
      mutate(
        day = time / 24,
        VL_display = pmax(VL_out, 2.0)  # Clamp at LLOQ for display
      )
    
    return(sim)
  })
  
  # Comparison data (multiple AUCs)
  compare_data <- reactive({
    shiny::req(input$vl0, input$duration)
    
    auc_values <- c(20, 30, 40, 50, 60, 70)
    
    all_sims <- lapply(auc_values, function(auc) {
      mod %>%
        param(AUC24 = auc, VL0 = input$vl0) %>%
        mrgsim(end = input$duration * 24, delta = 6) %>%
        as.data.frame() %>%
        mutate(
          day = time / 24,
          AUC_label = paste0("AUC = ", auc, " mg·h/L"),
          AUC = auc,
          VL_display = pmax(VL_out, 2.0)
        )
    })
    
    bind_rows(all_sims)
  })
  
  # ── Metrics ──
  output$auc_display <- renderText({
    sprintf("%.0f", input$auc24)
  })
  
  output$vl_day14 <- renderText({
    d <- sim_data()
    vl14 <- d %>% dplyr::filter(day >= 14) %>% dplyr::slice(1) %>% dplyr::pull(VL_out)
    if (length(vl14) == 0) vl14 <- NA
    if (is.na(vl14)) return("--")
    if (vl14 < 2.0) return("<LLOQ")
    sprintf("%.1f", vl14)
  })
  
  output$vl_day28 <- renderText({
    d <- sim_data()
    vl28 <- d %>% dplyr::filter(day >= 28) %>% dplyr::slice(1) %>% dplyr::pull(VL_out)
    if (length(vl28) == 0) vl28 <- NA
    if (is.na(vl28)) return("--")
    if (vl28 < 2.0) return("<LLOQ")
    sprintf("%.1f", vl28)
  })
  
  output$time_lloq <- renderText({
    d <- sim_data()
    below_lloq <- d %>% dplyr::filter(VL_out < 2.0)
    if (nrow(below_lloq) > 0) {
      first_day <- min(below_lloq$day)
      sprintf("%.0f", first_day)
    } else {
      ">56"
    }
  })
  
  # ── Main PD Plot ──
  output$pdPlot <- renderPlot({
    d <- sim_data()
    shiny::req(nrow(d) > 0)
    
    p <- ggplot(d, aes(x = day, y = VL_display)) +
      geom_line(color = "#8b5cf6", linewidth = 1.2) +
      labs(
        x = "Time (days)",
        y = expression("CMV Viral Load (log"[10]*" copies/mL)"),
        title = paste0("Ganciclovir AUC₀₋₂₄ = ", input$auc24, " mg·h/L")
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank()
      )
    
    # Add LLOQ threshold
    if (input$show_lloq) {
      p <- p + 
        geom_hline(yintercept = 2.0, linetype = "dashed", color = "#ef4444", linewidth = 0.8) +
        annotate("text", x = max(d$day) * 0.95, y = 2.15, 
                 label = "LLOQ (200 copies/mL)", hjust = 1, color = "#ef4444", size = 3.5)
    }
    
    # Add initial VL marker
    p <- p + 
      geom_point(data = d %>% dplyr::slice(1), aes(x = day, y = VL_display), 
                 color = "#f59e0b", size = 3) +
      annotate("text", x = 1, y = input$vl0 + 0.2, 
               label = paste0("Baseline: ", input$vl0), hjust = 0, color = "#f59e0b", size = 3.5)
    
    if (input$log_scale) {
      p <- p + scale_y_log10()
    } else {
      p <- p + ylim(1.5, max(input$vl0 + 0.5, 5))
    }
    
    p
  })
  
  # ── Comparison Plot ──
  output$comparePlot <- renderPlot({
    d <- compare_data()
    shiny::req(nrow(d) > 0)
    
    # Highlight therapeutic range
    d <- d %>%
      mutate(
        in_target = AUC >= 40 & AUC <= 60,
        AUC_label = factor(AUC_label, levels = paste0("AUC = ", c(20, 30, 40, 50, 60, 70), " mg·h/L"))
      )
    
    p <- ggplot(d, aes(x = day, y = VL_display, color = AUC_label, linetype = in_target)) +
      geom_line(linewidth = 1) +
      geom_hline(yintercept = 2.0, linetype = "dashed", color = "#ef4444", linewidth = 0.5) +
      scale_color_manual(
        values = c(
          "AUC = 20 mg·h/L" = "#94a3b8",
          "AUC = 30 mg·h/L" = "#64748b",
          "AUC = 40 mg·h/L" = "#10b981",
          "AUC = 50 mg·h/L" = "#8b5cf6",
          "AUC = 60 mg·h/L" = "#3b82f6",
          "AUC = 70 mg·h/L" = "#06b6d4"
        )
      ) +
      scale_linetype_manual(values = c("FALSE" = "dotted", "TRUE" = "solid"), guide = "none") +
      labs(
        x = "Time (days)",
        y = expression("CMV Viral Load (log"[10]*" copies/mL)"),
        title = "CMV Viral Load Decline by Ganciclovir Exposure",
        subtitle = "Solid lines = Target range (40-60 mg·h/L); Red dashed = LLOQ",
        color = "Exposure"
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(face = "bold"),
        legend.position = "right"
      )
    
    p
  })
  
  # ── Probability Tables ──
  output$probTable14 <- renderTable({
    # From paper Table 3: Probability of -1 log decrease after 2 weeks
    data.frame(
      `AUC₀₋₂₄ (mg·h/L)` = c(10, 20, 30, 40, 50, 60, 70, 80),
      `Probability (%)` = c(53.4, 73.1, 83.1, 88.8, 90.6, 92.6, 94.2, 95.5),
      check.names = FALSE
    )
  }, striped = TRUE, hover = TRUE)
  
  output$probTable28 <- renderTable({
    # From paper Table 4: Probability of undetectable VL at day 28
    data.frame(
      `AUC₀₋₂₄ (mg·h/L)` = c(10, 20, 30, 40, 50, 60, 70, 80),
      `Probability (%)` = c(52.2, 70.7, 80.2, 86.1, 89.2, 91.4, 92.8, 94.6),
      check.names = FALSE
    )
  }, striped = TRUE, hover = TRUE)
}

shinyApp(ui = ui, server = server)