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

# ── Torsemide PopPK Model ──
# Reference: Jeong et al. Pharmaceutics 2022, 14, 771
# 2-compartment model with first-order absorption + lag-time
# Covariates: CYP2C9 (on CL/F, CL2/F), OATP1B1 (on V/F)
# Population: 112 healthy Korean males

model_code <- '
$PARAM @annotated
tvCLF   : 1.740  : Typical CL/F (L/h) - IM reference
tvVF    : 1.027  : Typical V/F central (L) - ET reference
tvV2F   : 3.914  : Typical V2/F peripheral (L)
tvCL2F  : 0.828  : Typical CL2/F (L/h) - IM reference
tvKA    : 0.861  : Absorption rate constant (1/h)
tvTLAG  : 0.274  : Absorption lag time (h)
dCLF_CYP2C9EM  : 0.510  : CL/F increase for CYP2C9 EM
dCL2F_CYP2C9EM : 0.365  : CL2/F increase for CYP2C9 EM
dVF_OATP1B1IT  : -0.410 : V/F change for OATP1B1 IT
dVF_OATP1B1PT  : -0.646 : V/F change for OATP1B1 PT
CYP2C9_EM  : 1   : CYP2C9 phenotype (1=EM, 0=IM)
OATP1B1    : 0   : OATP1B1 phenotype (0=ET, 1=IT, 2=PT)

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

$MAIN
double CLF = tvCLF * (1.0 + dCLF_CYP2C9EM * CYP2C9_EM);
double VF  = tvVF;
if(OATP1B1 == 1) VF = tvVF * (1.0 + dVF_OATP1B1IT);
if(OATP1B1 == 2) VF = tvVF * (1.0 + dVF_OATP1B1PT);
double V2F  = tvV2F;
double CL2F = tvCL2F * (1.0 + dCL2F_CYP2C9EM * CYP2C9_EM);
double KA   = tvKA;
ALAG_DEPOT  = tvTLAG;

$ODE
dxdt_DEPOT  = -KA * DEPOT;
dxdt_CENT   = KA * DEPOT - (CLF/VF) * CENT - (CL2F/VF) * CENT + (CL2F/V2F) * PERIPH;
dxdt_PERIPH = (CL2F/VF) * CENT - (CL2F/V2F) * PERIPH;

$TABLE
double CP = CENT / VF;
CP = CP > 0 ? CP : 0;

$CAPTURE @annotated
CP   : Plasma concentration (mg/L)
CLF  : Individual CL/F (L/h)
VF   : Individual V/F central (L)
CL2F : Individual CL2/F (L/h)
'

mod <- mcode("torsemide_popPK", model_code)

# ── Theme ──
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; }
")

# ── UI ──
ui <- page_sidebar(
  title = "Torsemide Population PK Simulator",
  theme = app_theme,
  
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    
    h6("Dosing"),
    selectInput("dose", "Dose (mg)", choices = c(5, 10, 20), selected = 10),
    numericInput("n_doses", "Number of Doses", value = 1, min = 1, max = 30),
    numericInput("interval", "Dosing Interval (h)", value = 24, min = 4, max = 48),
    
    hr(),
    h6("Pharmacogenomics"),
    selectInput("cyp2c9", "CYP2C9 Phenotype",
      choices = c("EM (Extensive Metabolizer)" = "EM",
                  "IM (Intermediate Metabolizer)" = "IM"),
      selected = "EM"),
    selectInput("oatp1b1", "OATP1B1 Phenotype",
      choices = c("ET (Extensive Transporter)" = "ET",
                  "IT (Intermediate Transporter)" = "IT",
                  "PT (Poor Transporter)" = "PT"),
      selected = "ET"),
    
    hr(),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE),
    checkboxInput("show_genotype_comparison", "Compare All Genotypes", value = FALSE),
    
    hr(),
    card(
      card_header("Population PK Parameters", class = "bg-light"),
      p(strong("CL/F:"), "1.74 L/h (IM) / 2.63 L/h (EM)"),
      p(strong("V/F:"), "1.03 L (ET) / 0.61 L (IT) / 0.36 L (PT)"),
      p(strong("V2/F:"), "3.91 L"),
      p(strong("CL2/F:"), "0.83 L/h (IM) / 1.13 L/h (EM)"),
      p(strong("Ka:"), "0.861 /h"),
      p(strong("Tlag:"), "0.274 h"),
      p(em("Source: Jeong et al. 2022"))
    ),
    
    card(
      card_header("Therapeutic Context", class = "bg-light"),
      p("Torsemide is a loop diuretic used for edema in CHF, liver disease, and kidney disease."),
      p(strong("Typical doses:"), "5-20 mg oral, once daily"),
      p(strong("Bioavailability:"), "~80-90%"),
      p(strong("Protein binding:"), ">99%"),
      p(strong("Half-life:"), "~2.5-3.5 h"),
      p(strong("Metabolism:"), "Primarily hepatic via CYP2C9")
    )
  ),
  
  # Metrics row
  layout_column_wrap(
        width = 1/4,
        fill = FALSE,
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("cmax_val")),
      div(class = "metric-label", "Cmax (mg/L)")
    ),
    div(class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("tmax_val")),
      div(class = "metric-label", "Tmax (h)")
    ),
    div(class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("auc_val")),
      div(class = "metric-label", "AUC (mg\u00B7h/L)")
    ),
    div(class = "metric-card metric-info",
      div(class = "metric-value", textOutput("thalf_val")),
      div(class = "metric-label", "t\u00BD (h)")
    )
  ),
  
  # PK plot
  card(
    card_header("Concentration-Time Profile"),
    full_screen = TRUE,
    plotOutput("pkPlot", height = "450px")
  ),
  
  # Tabbed content
  navset_card_tab(
    title = "Analysis",
    nav_panel("Genotype Comparison",
      plotOutput("genotypePlot", height = "400px")
    ),
    nav_panel("Parameter Table",
      tableOutput("paramTable")
    ),
    nav_panel("Model Information",
      markdown(paste0(
        "## Torsemide Population PK Model\n\n",
        "**Reference:** Jeong SH, Jang JH, Cho HY, Lee YB. Population Pharmacokinetic (Pop-PK) Analysis ",
        "of Torsemide in Healthy Korean Males Considering CYP2C9 and OATP1B1 Genetic Polymorphisms. ",
        "*Pharmaceutics* 2022; 14(4):771. [DOI: 10.3390/pharmaceutics14040771](https://doi.org/10.3390/pharmaceutics14040771)\n\n",
        "### Model Structure\n",
        "- **Type:** Two-compartment with first-order absorption and lag-time\n",
        "- **Software:** Phoenix NLME (confirmed with NONMEM)\n",
        "- **Population:** 112 healthy Korean males\n",
        "- **Doses studied:** 5, 10, and 20 mg oral single dose\n",
        "- **Estimation:** FOCE with \u03B7-\u03B5 interaction\n\n",
        "### Covariates\n",
        "- **CYP2C9** affects CL/F and CL2/F: EM has ~51% higher CL/F and ~37% higher CL2/F vs IM\n",
        "- **OATP1B1** affects V/F: IT has ~41% lower V/F, PT has ~65% lower V/F vs ET\n\n",
        "### Key Findings\n",
        "- Genetic polymorphisms in CYP2C9 and OATP1B1 significantly affect torsemide PK\n",
        "- CYP2C9 genotype impacts clearance (metabolism)\n",
        "- OATP1B1 genotype impacts volume of distribution (hepatic uptake)\n",
        "- The two processes appear to act independently\n\n",
        "### References\n",
        "- [FDA Label (Demadex/Torsemide)](https://www.accessdata.fda.gov/drugsatfda_docs/label/2017/020136s027lbl.pdf)\n",
        "- [Full Paper on PubMed Central](https://pmc.ncbi.nlm.nih.gov/articles/PMC9028991/)\n",
        "- [PBPK-PD Modeling of Torsemide (Pharmaceutics 2022)](https://www.mdpi.com/1999-4923/14/12/2720)\n\n",
        "### Disclaimer\n",
        "*For research and educational purposes only. Not for clinical decision-making.*"
      ))
    )
  ),
  
  # 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"),
    " \u2022 Built by Sunny \u2600\uFE0F (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) {
  
  sim_data <- reactive({
    dose_mg <- as.numeric(input$dose)
    n_doses <- input$n_doses
    ii <- input$interval
    
    cyp2c9_em <- ifelse(input$cyp2c9 == "EM", 1, 0)
    oatp1b1_val <- switch(input$oatp1b1, "ET" = 0, "IT" = 1, "PT" = 2)
    
    end_time <- ii * n_doses + 12
    
    ev_obj <- ev(amt = dose_mg, cmt = 1, ii = ii, addl = n_doses - 1)
    
    out <- mod %>%
      param(CYP2C9_EM = cyp2c9_em, OATP1B1 = oatp1b1_val) %>%
      ev(ev_obj) %>%
      mrgsim(end = end_time, delta = 0.05) %>%
      as.data.frame()
    
    out
  })
  
  # Metrics
  output$cmax_val <- renderText({
    d <- sim_data()
    sprintf("%.2f", max(d$CP, na.rm = TRUE))
  })
  
  output$tmax_val <- renderText({
    d <- sim_data()
    sprintf("%.2f", d$time[which.max(d$CP)])
  })
  
  output$auc_val <- renderText({
    d <- sim_data()
    auc <- sum(diff(d$time) * (head(d$CP, -1) + tail(d$CP, -1)) / 2)
    sprintf("%.1f", auc)
  })
  
  output$thalf_val <- renderText({
    d <- sim_data()
    # Calculate from model params
    cyp2c9_em <- ifelse(input$cyp2c9 == "EM", 1, 0)
    oatp1b1_val <- switch(input$oatp1b1, "ET" = 0, "IT" = 1, "PT" = 2)
    clf <- 1.740 * (1 + 0.510 * cyp2c9_em)
    vf <- 1.027
    if(oatp1b1_val == 1) vf <- 1.027 * (1 - 0.410)
    if(oatp1b1_val == 2) vf <- 1.027 * (1 - 0.646)
    v2f <- 3.914
    cl2f <- 0.828 * (1 + 0.365 * cyp2c9_em)
    
    # Terminal half-life from eigenvalues
    a <- clf/vf + cl2f/vf + cl2f/v2f
    b <- (clf * cl2f) / (vf * v2f)
    disc <- pmax(a^2 - 4*b, 0)
    beta <- (a - sqrt(disc)) / 2
    thalf <- 0.693 / beta
    sprintf("%.1f", thalf)
  })
  
  # PK Plot
  output$pkPlot <- renderPlot({
    d <- sim_data()
    
    if(input$show_genotype_comparison) {
      # Simulate all 6 genotype combos
      combos <- expand.grid(CYP2C9 = c("EM", "IM"), OATP1B1 = c("ET", "IT", "PT"))
      all_data <- lapply(1:nrow(combos), function(i) {
        cyp_em <- ifelse(combos$CYP2C9[i] == "EM", 1, 0)
        oatp_val <- switch(as.character(combos$OATP1B1[i]), "ET" = 0, "IT" = 1, "PT" = 2)
        dose_mg <- as.numeric(input$dose)
        ev_obj <- ev(amt = dose_mg, cmt = 1, ii = input$interval, addl = input$n_doses - 1)
        end_time <- input$interval * input$n_doses + 12
        out <- mod %>%
          param(CYP2C9_EM = cyp_em, OATP1B1 = oatp_val) %>%
          ev(ev_obj) %>%
          mrgsim(end = end_time, delta = 0.1) %>%
          as.data.frame()
        out$Genotype <- paste0("CYP2C9:", combos$CYP2C9[i], " / OATP1B1:", combos$OATP1B1[i])
        out
      })
      d_all <- do.call(rbind, all_data)
      
      d_plot <- if(input$log_scale) d_all %>% dplyr::filter(CP > 0.001) else d_all
      
      p <- ggplot(d_plot, aes(x = time, y = CP, color = Genotype)) +
        geom_line(linewidth = 0.8) +
        labs(x = "Time (h)", y = "Concentration (mg/L)",
             title = paste0("Torsemide ", input$dose, " mg - All Genotype Combinations")) +
        theme_minimal(base_size = 14) +
        theme(legend.position = "bottom", legend.text = element_text(size = 9)) +
        guides(color = guide_legend(ncol = 2))
      
    } else {
      d_plot <- if(input$log_scale) d %>% dplyr::filter(CP > 0.001) else d
      
      p <- ggplot(d_plot, aes(x = time, y = CP)) +
        geom_line(linewidth = 1.2, color = "#8b5cf6") +
        labs(x = "Time (h)", y = "Concentration (mg/L)",
             title = paste0("Torsemide ", input$dose, " mg PK Profile (",
                           input$cyp2c9, "/", input$oatp1b1, ")")) +
        theme_minimal(base_size = 14)
    }
    
    if(input$log_scale) p <- p + scale_y_log10()
    p
  })
  
  # Genotype comparison plot
  output$genotypePlot <- renderPlot({
    dose_mg <- as.numeric(input$dose)
    combos <- expand.grid(CYP2C9 = c("EM", "IM"), OATP1B1 = c("ET", "IT", "PT"))
    
    metrics <- lapply(1:nrow(combos), function(i) {
      cyp_em <- ifelse(combos$CYP2C9[i] == "EM", 1, 0)
      oatp_val <- switch(as.character(combos$OATP1B1[i]), "ET" = 0, "IT" = 1, "PT" = 2)
      ev_obj <- ev(amt = dose_mg, cmt = 1)
      out <- mod %>%
        param(CYP2C9_EM = cyp_em, OATP1B1 = oatp_val) %>%
        ev(ev_obj) %>%
        mrgsim(end = 24, delta = 0.05) %>%
        as.data.frame()
      data.frame(
        Genotype = paste0(combos$CYP2C9[i], "/", combos$OATP1B1[i]),
        Cmax = max(out$CP),
        AUC = sum(diff(out$time) * (head(out$CP, -1) + tail(out$CP, -1)) / 2)
      )
    })
    metrics_df <- do.call(rbind, metrics)
    
    p1 <- ggplot(metrics_df, aes(x = Genotype, y = Cmax, fill = Genotype)) +
      geom_col(show.legend = FALSE) +
      labs(x = "", y = "Cmax (mg/L)", title = paste0("Cmax by Genotype (", dose_mg, " mg)")) +
      theme_minimal(base_size = 14) +
      scale_fill_brewer(palette = "Set2") +
      theme(axis.text.x = element_text(angle = 30, hjust = 1))
    
    p1
  })
  
  # Parameter table
  output$paramTable <- renderTable({
    data.frame(
      Parameter = c("tvCL/F", "tvV/F", "tvV2/F", "tvCL2/F", "tvKa", "tvTlag",
                    "dCL/F (CYP2C9 EM)", "dCL2/F (CYP2C9 EM)", 
                    "dV/F (OATP1B1 IT)", "dV/F (OATP1B1 PT)"),
      Estimate = c(1.740, 1.027, 3.914, 0.828, 0.861, 0.274,
                   0.510, 0.365, -0.410, -0.646),
      `RSE (%)` = c(4.19, 7.92, 2.65, 6.71, 2.93, 6.64,
                    12.37, 15.23, 31.54, 31.56),
      `IIV (%)` = c(17.1, 74.9, 18.9, 14.9, NA, 58.4,
                    NA, NA, NA, NA),
      Unit = c("L/h", "L", "L", "L/h", "1/h", "h",
               "fraction", "fraction", "fraction", "fraction"),
      check.names = FALSE
    )
  }, digits = 3)
}

shinyApp(ui = ui, server = server)