# ══════════════════════════════════════════════════════════════════════════════
# Simvastatin Population PK Simulator
# Based on: Jin et al. Pharm Res (2014) - PopPK of Simvastatin & Simvastatin Acid
# Built by: PKPDBuilder.com
# ══════════════════════════════════════════════════════════════════════════════

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

# ══════════════════════════════════════════════════════════════════════════════
# mrgsolve Model Code
# Two-compartment model for simvastatin with linked metabolite (simvastatin acid)
# Parameters from Jin et al. 2014 - Korean healthy adult males
# ══════════════════════════════════════════════════════════════════════════════

model_code <- '
$PARAM @annotated
// Simvastatin (parent) PK parameters
CL   : 571    : Clearance (L/h)
V1   : 199    : Central volume of distribution (L)
V2   : 2710   : Peripheral volume of distribution (L)
Q    : 199    : Intercompartmental clearance (L/h)
KA   : 0.964  : Absorption rate constant (1/h)
TLAG : 0.15   : Absorption lag time (h)
F1   : 1.0    : Bioavailability (fraction)

// Metabolite (simvastatin acid) parameters
// Note: VM was fixed at 1L in original study for identifiability
// Adjusted here for clinically realistic concentration ratios
FM   : 0.236  : Fraction metabolized to simvastatin acid
VM   : 100.0  : Metabolite central volume (L) - scaled for realistic conc
CLM  : 35.0   : Metabolite clearance (L/h) - scaled from 0.035 * 100
VM2  : 200.0  : Metabolite peripheral volume (L)
QM   : 20.0   : Metabolite intercompartmental clearance (L/h)

// Patient covariate (no significant covariates in this study)
WT   : 70     : Body weight (kg)

$CMT @annotated
DEPOT  : Oral absorption depot (mg)
CENT   : Parent central compartment (mg)
PERIPH : Parent peripheral compartment (mg)
CENTM  : Metabolite central compartment (mg)
PERIPHM: Metabolite peripheral compartment (mg)

$MAIN
// Set bioavailability and lag time
F_DEPOT = F1;
ALAG_DEPOT = TLAG;

$ODE
// Simvastatin disposition
double CP = CENT / V1;
double k10 = CL / V1;       // Elimination rate constant
double k12 = Q / V1;        // Transfer to peripheral
double k21 = Q / V2;        // Transfer from peripheral
double k_met = (CL * FM) / V1;  // Metabolic rate to simvastatin acid

// Simvastatin acid disposition
double CM = CENTM / VM;
double k10m = CLM / VM;
double k12m = QM / VM;
double k21m = QM / VM2;

// ODEs for parent drug
dxdt_DEPOT  = -KA * DEPOT;
dxdt_CENT   = KA * DEPOT - (k10 + k12) * CENT + k21 * PERIPH;
dxdt_PERIPH = k12 * CENT - k21 * PERIPH;

// ODEs for metabolite (simvastatin acid) - formed from parent metabolism
// Note: FM fraction of CL goes to metabolite, rest is other elimination pathways
double rate_formation = k_met * CENT;  // Rate of simvastatin acid formation
dxdt_CENTM   = rate_formation - (k10m + k12m) * CENTM + k21m * PERIPHM;
dxdt_PERIPHM = k12m * CENTM - k21m * PERIPHM;

$TABLE
// Output concentrations
double PARENT = CENT / V1;       // Simvastatin concentration (mg/L = μg/mL)
double PARENT_ng = PARENT * 1000; // Convert to ng/mL
double METAB = CENTM / VM;       // Simvastatin acid concentration
double METAB_ng = METAB * 1000;  // Convert to ng/mL

// Clamp to prevent negative values
PARENT_ng = PARENT_ng > 0 ? PARENT_ng : 0;
METAB_ng = METAB_ng > 0 ? METAB_ng : 0;

$CAPTURE @annotated
PARENT_ng : Simvastatin concentration (ng/mL)
METAB_ng  : Simvastatin acid concentration (ng/mL)
CP        : Parent conc (mg/L)
CM        : Metabolite conc (mg/L)
'

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

# ══════════════════════════════════════════════════════════════════════════════
# App 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; }
")

# ══════════════════════════════════════════════════════════════════════════════
# UI
# ══════════════════════════════════════════════════════════════════════════════

ui <- page_sidebar(
  title = "Simvastatin PopPK Simulator",
  theme = app_theme,
  
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 320,
    
    h6("Dosing"),
    selectInput("dose", "Dose (mg)", 
                choices = c(10, 20, 30, 40, 60, 80),
                selected = 40),
    
    numericInput("n_doses", "Number of Doses", value = 1, min = 1, max = 30),
    numericInput("interval", "Dosing Interval (hours)", value = 24, min = 12, max = 48),
    
    hr(),
    h6("Simulation"),
    numericInput("sim_hours", "Simulation Duration (hours)", value = 48, min = 24, max = 168),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE),
    checkboxInput("show_metabolite", "Show Simvastatin Acid", value = TRUE),
    
    hr(),
    card(
      card_header("Model Parameters", class = "bg-light"),
      p(strong("CL:"), "571 L/h"),
      p(strong("V1:"), "199 L (central)"),
      p(strong("V2:"), "2,710 L (peripheral)"),
      p(strong("Q:"), "199 L/h"),
      p(strong("Ka:"), "0.964 h⁻¹"),
      p(strong("F_M:"), "0.236 (→ SVA)"),
      hr(),
      p(em("Population: 133 healthy Korean males"))
    ),
    
    hr(),
    card(
      card_header("Clinical Context", class = "bg-light"),
      p(strong("Indication:"), "Hypercholesterolemia"),
      p(strong("Usual doses:"), "10-40 mg/day"),
      p(strong("Max dose:"), "80 mg/day (restricted)"),
      tags$a(href = "https://www.accessdata.fda.gov/drugsatfda_docs/label/2023/019766s104lbl.pdf", 
             target = "_blank", "FDA Prescribing Information ↗")
    )
  ),
  
  # ── Metrics Row ──
  layout_column_wrap(
        width = 1/4,
        fill = FALSE,
    
    div(class = "metric-card metric-success",
        div(class = "metric-value", textOutput("cmax_parent")),
        div(class = "metric-label", "Cmax Simvastatin (ng/mL)")
    ),
    div(class = "metric-card metric-info",
        div(class = "metric-value", textOutput("cmax_metab")),
        div(class = "metric-label", "Cmax SVA (ng/mL)")
    ),
    div(class = "metric-card metric-warning",
        div(class = "metric-value", textOutput("tmax_parent")),
        div(class = "metric-label", "Tmax (h)")
    ),
    div(class = "metric-card metric-primary",
        div(class = "metric-value", textOutput("auc_parent")),
        div(class = "metric-label", "AUC₀₋₂₄ (ng·h/mL)")
    )
  ),
  
  # ── Main Plot ──
  card(
    card_header("Concentration-Time Profile"),
    full_screen = TRUE,
    plotOutput("pk_plot", height = "450px")
  ),
  
  # ── Tabbed Content ──
  navset_card_tab(
    title = "Analysis",
    
    nav_panel("Dose Comparison",
      plotOutput("compare_plot", height = "400px")
    ),
    
    nav_panel("PK Metrics Table",
      h5("Single-Dose PK Metrics"),
      tableOutput("metrics_table"),
      hr(),
      h5("Multi-Dose Trough Levels (Cmin) at Steady State"),
      p(em("10 doses QD to approximate steady state, measured just before next dose")),
      tableOutput("trough_table")
    ),
    
    nav_panel("R Code",
      card(
        card_header("mrgsolve Model Code", class = "bg-light"),
        tags$pre(
          style = "background: #1e293b; color: #e2e8f0; padding: 20px; border-radius: 8px; font-size: 12px; overflow-x: auto; max-height: 600px;",
          model_code
        )
      ),
      card(
        card_header("Simulation Code Example", class = "bg-light"),
        tags$pre(
          style = "background: #1e293b; color: #e2e8f0; padding: 20px; border-radius: 8px; font-size: 12px; overflow-x: auto;",
'# Load libraries
library(mrgsolve)
library(dplyr)
library(ggplot2)

# Compile the model (paste model_code above)
mod <- mcode("simvastatin", model_code)

# Single dose simulation
ev_single <- ev(amt = 40, cmt = 1)  # 40 mg dose
sim_single <- mod %>%
  ev(ev_single) %>%
  mrgsim(end = 24, delta = 0.1) %>%
  as.data.frame()

# Get Cmin at 24h
cmin_24h <- sim_single$PARENT_ng[which.min(abs(sim_single$time - 24))]

# Multi-dose simulation (10 doses QD)
ev_multi <- ev(amt = 40, cmt = 1, ii = 24, addl = 9)
sim_multi <- mod %>%
  ev(ev_multi) %>%
  mrgsim(end = 240, delta = 0.1) %>%
  as.data.frame()

# Get trough levels at various timepoints
trough_times <- c(24, 48, 216, 240)  # Just before doses 2, 3, 10, and after dose 10
troughs <- sapply(trough_times, function(t) {
  sim_multi$PARENT_ng[which.min(abs(sim_multi$time - t))]
})

# Plot
ggplot(sim_multi, aes(x = time, y = PARENT_ng)) +
  geom_line(color = "#8b5cf6", linewidth = 1) +
  labs(title = "Simvastatin 40 mg QD - Multi-Dose Profile",
       x = "Time (hours)", y = "Concentration (ng/mL)") +
  theme_minimal()'
        )
      )
    ),
    
    nav_panel("Model Information",
      markdown("
## Model Structure

**Two-compartment model** for simvastatin (parent) with linked metabolite:

- **Simvastatin**: 2-compartment disposition with first-order absorption
- **Simvastatin Acid (SVA)**: Active metabolite, 2-compartment disposition

### Key Features
- Complex absorption kinetics simplified to first-order (Ka = 0.964 h⁻¹)
- Original study described 3 parallel absorption processes with multiple peaks
- Extensive first-pass metabolism (< 5% systemic bioavailability)
- Reversible interconversion between simvastatin and SVA

### Clinical Relevance
- Simvastatin is a **prodrug** — SVA is the active HMG-CoA reductase inhibitor
- Evening dosing recommended (cholesterol synthesis peaks at night)
- 80 mg dose restricted due to myopathy risk

### Reference
Jin SJ, Bae KS, Cho SH, et al. **Population Pharmacokinetic Analysis of Simvastatin and its Active Metabolite with the Characterization of Atypical Complex Absorption Kinetics.** *Pharm Res.* 2014. DOI: 10.1007/s11095-013-1284-0

### Population
- 133 healthy Korean male volunteers
- Age: 25.0 ± 2.6 years
- Weight: 69.8 ± 6.7 kg
- BMI: 22.8 ± 1.7 kg/m²
      ")
    )
  ),
  
  # ── 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$dose, input$n_doses, input$interval, input$sim_hours)
    
    dose_mg <- as.numeric(input$dose)
    n_doses <- input$n_doses
    interval <- input$interval
    sim_end <- input$sim_hours
    
    # Build dosing event
    ev <- ev(amt = dose_mg, cmt = 1, ii = interval, addl = n_doses - 1)
    
    # Run simulation
    out <- mod %>%
      ev(ev) %>%
      mrgsim(end = sim_end, delta = 0.1) %>%
      as.data.frame()
    
    out
  })
  
  # ── Calculate PK metrics ──
  metrics <- reactive({
    d <- sim_data()
    shiny::req(nrow(d) > 0)
    
    # Parent (simvastatin)
    cmax_parent <- max(d$PARENT_ng, na.rm = TRUE)
    tmax_parent <- d$time[which.max(d$PARENT_ng)]
    
    # AUC 0-24h (trapezoidal)
    d24 <- d %>% dplyr::filter(time <= 24)
    if (nrow(d24) > 1) {
      auc24 <- sum(diff(d24$time) * (head(d24$PARENT_ng, -1) + tail(d24$PARENT_ng, -1)) / 2)
    } else {
      auc24 <- 0
    }
    
    # Metabolite (simvastatin acid)
    cmax_metab <- max(d$METAB_ng, na.rm = TRUE)
    tmax_metab <- d$time[which.max(d$METAB_ng)]
    
    list(
      cmax_parent = cmax_parent,
      tmax_parent = tmax_parent,
      auc24_parent = auc24,
      cmax_metab = cmax_metab,
      tmax_metab = tmax_metab
    )
  })
  
  # ── Metric outputs ──
  output$cmax_parent <- renderText({
    m <- metrics()
    sprintf("%.2f", m$cmax_parent)
  })
  
  output$cmax_metab <- renderText({
    m <- metrics()
    sprintf("%.2f", m$cmax_metab)
  })
  
  output$tmax_parent <- renderText({
    m <- metrics()
    sprintf("%.1f", m$tmax_parent)
  })
  
  output$auc_parent <- renderText({
    m <- metrics()
    sprintf("%.1f", m$auc24_parent)
  })
  
  # ── Main PK Plot ──
  output$pk_plot <- renderPlot({
    d <- sim_data()
    shiny::req(nrow(d) > 0)
    
    dose_mg <- as.numeric(input$dose)
    
    # Prepare data for plotting
    if (input$show_metabolite) {
      d_long <- d %>%
        dplyr::select(time, PARENT_ng, METAB_ng) %>%
        pivot_longer(cols = c(PARENT_ng, METAB_ng),
                     names_to = "Analyte",
                     values_to = "Concentration") %>%
        mutate(Analyte = ifelse(Analyte == "PARENT_ng", "Simvastatin", "Simvastatin Acid"))
    } else {
      d_long <- d %>%
        dplyr::select(time, PARENT_ng) %>%
        dplyr::rename(Concentration = PARENT_ng) %>%
        mutate(Analyte = "Simvastatin")
    }
    
    # Filter for log scale
    if (input$log_scale) {
      d_long <- d_long %>% dplyr::filter(Concentration > 0.001)
    }
    
    p <- ggplot(d_long, aes(x = time, y = Concentration, color = Analyte)) +
      geom_line(linewidth = 1.2) +
      labs(
        title = paste0("Simvastatin ", dose_mg, " mg - Concentration vs Time"),
        subtitle = "Population PK simulation (Jin et al. 2014)",
        x = "Time (hours)",
        y = "Concentration (ng/mL)",
        color = "Analyte"
      ) +
      scale_color_manual(values = c("Simvastatin" = "#8b5cf6", "Simvastatin Acid" = "#10b981")) +
      theme_minimal(base_size = 14) +
      theme(
        legend.position = "bottom",
        plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank()
      )
    
    if (input$log_scale) {
      p <- p + scale_y_log10()
    }
    
    p
  })
  
  # ── Dose Comparison Plot ──
  output$compare_plot <- renderPlot({
    doses <- c(10, 20, 40, 60, 80)
    
    all_sims <- lapply(doses, function(dose_mg) {
      ev <- ev(amt = dose_mg, cmt = 1)
      out <- mod %>%
        ev(ev) %>%
        mrgsim(end = 24, delta = 0.1) %>%
        as.data.frame() %>%
        mutate(Dose = paste0(dose_mg, " mg"))
      out
    })
    
    d_compare <- bind_rows(all_sims)
    d_compare$Dose <- factor(d_compare$Dose, levels = paste0(doses, " mg"))
    
    ggplot(d_compare, aes(x = time, y = PARENT_ng, color = Dose)) +
      geom_line(linewidth = 1) +
      labs(
        title = "Simvastatin Exposure Across Doses",
        subtitle = "Single-dose simulation (10-80 mg)",
        x = "Time (hours)",
        y = "Simvastatin Concentration (ng/mL)",
        color = "Dose"
      ) +
      scale_color_viridis_d(option = "plasma", end = 0.85) +
      theme_minimal(base_size = 14) +
      theme(
        legend.position = "right",
        plot.title = element_text(face = "bold")
      )
  })
  
  # ── PK Metrics Table (Single Dose) ──
  output$metrics_table <- renderTable({
    doses <- c(10, 20, 40, 60, 80)
    
    metrics_list <- lapply(doses, function(dose_mg) {
      ev <- ev(amt = dose_mg, cmt = 1)
      out <- mod %>%
        ev(ev) %>%
        mrgsim(end = 24, delta = 0.1) %>%
        as.data.frame()
      
      cmax <- max(out$PARENT_ng)
      tmax <- out$time[which.max(out$PARENT_ng)]
      auc <- sum(diff(out$time) * (head(out$PARENT_ng, -1) + tail(out$PARENT_ng, -1)) / 2)
      
      # Cmin at 24h (single dose trough)
      cmin_24h <- out$PARENT_ng[which.min(abs(out$time - 24))]
      
      cmax_m <- max(out$METAB_ng)
      tmax_m <- out$time[which.max(out$METAB_ng)]
      
      data.frame(
        `Dose (mg)` = dose_mg,
        `Cmax SVS (ng/mL)` = round(cmax, 2),
        `Tmax SVS (h)` = round(tmax, 1),
        `AUC₀₋₂₄ SVS` = round(auc, 1),
        `Cmin 24h (ng/mL)` = round(cmin_24h, 4),
        `Cmax SVA (ng/mL)` = round(cmax_m, 2),
        `Tmax SVA (h)` = round(tmax_m, 1),
        check.names = FALSE
      )
    })
    
    bind_rows(metrics_list)
  }, striped = TRUE, bordered = TRUE, hover = TRUE)
  
  # ── Multi-Dose Trough Table ──
  output$trough_table <- renderTable({
    doses <- c(10, 20, 40, 60, 80)
    
    trough_list <- lapply(doses, function(dose_mg) {
      # 10 doses QD to reach approximate steady state
      ev <- ev(amt = dose_mg, cmt = 1, ii = 24, addl = 9)
      out <- mod %>%
        ev(ev) %>%
        mrgsim(end = 240, delta = 0.1) %>%
        as.data.frame()
      
      # Get trough (Cmin) at specified timepoints
      # These are measured just before the next dose
      cmin_24h <- out$PARENT_ng[which.min(abs(out$time - 24))]
      cmin_48h <- out$PARENT_ng[which.min(abs(out$time - 48))]
      cmin_216h <- out$PARENT_ng[which.min(abs(out$time - 216))]
      cmin_240h <- out$PARENT_ng[which.min(abs(out$time - 240))]
      
      data.frame(
        `Dose (mg)` = dose_mg,
        `Cmin 24h (ng/mL)` = round(cmin_24h, 4),
        `Cmin 48h (ng/mL)` = round(cmin_48h, 4),
        `Cmin 216h (ng/mL)` = round(cmin_216h, 4),
        `Cmin 240h (ng/mL)` = round(cmin_240h, 4),
        check.names = FALSE
      )
    })
    
    bind_rows(trough_list)
  }, striped = TRUE, bordered = TRUE, hover = TRUE)
  
}

# ══════════════════════════════════════════════════════════════════════════════
# Run App
# ══════════════════════════════════════════════════════════════════════════════

shinyApp(ui = ui, server = server)