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

# ── Vortioxetine PK/PD Model ──
# Source: Naik et al. 2016, Basic Clin Pharmacol Toxicol 118:344-355
# 2-compartment PK + Emax PK/Efficacy (MADRS) + Nausea risk

model_code <- '
$PARAM @annotated
TVCL_EU  : 39.0   : CL/F EU (L/hr)
TVCL_US  : 51.0   : CL/F USA (L/hr)
TVCL_ROW : 38.0   : CL/F RoW (L/hr)
V2F      : 2900   : Central volume V2/F (L)
QF       : 23.0   : Intercompartmental CL Q/F (L/hr)
V3F      : 670    : Peripheral volume V3/F (L)
KA       : 0.14   : Absorption rate (1/hr)
HT_EFF   : 0.40   : Height effect on CL (L/hr/cm)
CRCL_EFF : 0.18   : CrCL effect on CL (L/hr per mL/min)
REGION   : 0      : Region (0=EU, 1=USA, 2=RoW)
HT       : 167    : Height (cm)
CRCL     : 106    : Creatinine clearance (mL/min)

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

$MAIN
double TVCL = TVCL_EU;
if (REGION == 1) TVCL = TVCL_US;
if (REGION == 2) TVCL = TVCL_ROW;

double CLi = TVCL + (CRCL - 106.0) * CRCL_EFF + (HT - 167.0) * HT_EFF;
if (CLi < 1.0) CLi = 1.0;

double K20 = CLi / V2F;
double K23 = QF / V2F;
double K32 = QF / V3F;

$ODE
dxdt_GUT    = -KA * GUT;
dxdt_CENT   =  KA * GUT - (K20 + K23) * CENT + K32 * PERIPH;
dxdt_PERIPH =  K23 * CENT - K32 * PERIPH;

$TABLE
double CP = CENT / V2F;
double CP_ngml = CP * 1000.0;

$CAPTURE @annotated
CP_ngml : Plasma concentration (ng/mL)
CLi     : Individual CL/F (L/hr)
'

mod <- mcode("vortioxetine", model_code)

# ── PK/Efficacy Emax function ──
calc_dmadrs <- function(cav_ngml, baseline_madrs, region_usa = FALSE) {
  ec50 <- 24.9 + (-0.6) * baseline_madrs
  if (ec50 < 1) ec50 <- 1
  emax <- 7.0
  if (region_usa) emax <- emax - 4.07
  e0 <- -13.2
  dmadrs <- e0 + (cav_ngml * emax) / (ec50 + cav_ngml)
  return(dmadrs)
}

# ── Nausea risk (logistic) ──
calc_nausea_risk <- function(cav_ngml) {
  if (cav_ngml <= 0) return(0)
  logit <- -1.6 + 0.2 * log(cav_ngml)
  prob <- exp(logit) / (1 + exp(logit))
  return(prob * 100)
}

# ── 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; }
  .metric-danger .metric-value { color: #ef4444; }
")

# ── UI ──
ui <- page_sidebar(
  title = "Vortioxetine (Trintellix\u00AE) PK/PD Simulator",
  theme = app_theme,

  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,

    h6("Patient Characteristics"),
    selectInput("region", "Geographic Region",
      choices = c("European Union" = "0", "United States" = "1", "Rest of World" = "2"),
      selected = "0"
    ),
    sliderInput("height", "Height (cm)", min = 150, max = 200, value = 167, step = 1),
    sliderInput("crcl", "Creatinine Clearance (mL/min)", min = 30, max = 200, value = 106, step = 1),
    sliderInput("baseline_madrs", "Baseline MADRS Score", min = 18, max = 45, value = 31, step = 1),

    hr(),
    h6("Dosing"),
    selectInput("dose", "Vortioxetine Dose (mg/day)",
      choices = c(5, 10, 15, 20),
      selected = 10
    ),
    numericInput("n_days", "Treatment Duration (days)", value = 56, min = 7, max = 180),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE),

    hr(),
    card(
      card_header("Model Parameters", class = "bg-light"),
      p(strong("CL/F:"), "39-51 L/hr (region-dependent)"),
      p(strong("V2/F:"), "2,900 L"),
      p(strong("Q/F:"), "23 L/hr"),
      p(strong("V3/F:"), "670 L"),
      p(strong("Ka:"), "0.14 hr\u207B\u00B9"),
      p(strong("t\u00BD:"), "~57 hr"),
      p(strong("EC50:"), "24.9 ng/mL"),
      p(strong("Emax:"), "7.0 MADRS points"),
      p(em("Approved doses: 5-20 mg/day"))
    )
  ),

  # ── Metrics Row ──
  layout_column_wrap(
    width = 1/6,
    fill = FALSE,
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("cmax_val")),
      div(class = "metric-label", "Cmax,ss (ng/mL)")
    ),
    div(class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("ctrough_val")),
      div(class = "metric-label", "Ctrough,ss (ng/mL)")
    ),
    div(class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("cav_val")),
      div(class = "metric-label", "Cav,ss (ng/mL)")
    ),
    div(class = "metric-card metric-info",
      div(class = "metric-value", textOutput("thalf_val")),
      div(class = "metric-label", "t\u00BD (hr)")
    ),
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("dmadrs_val")),
      div(class = "metric-label", "\u0394MADRS")
    ),
    div(class = "metric-card metric-danger",
      div(class = "metric-value", textOutput("nausea_val")),
      div(class = "metric-label", "Nausea Risk (%)")
    )
  ),

  # ── PK Plot ──
  card(
    card_header("Plasma Concentration-Time Profile"),
    full_screen = TRUE,
    plotOutput("pk_plot", height = "420px")
  ),

  # ── Tabbed Content ──
  navset_card_tab(
    title = "Exposure-Response Analysis",
    nav_panel("Efficacy (MADRS)",
      plotOutput("er_efficacy_plot", height = "400px")
    ),
    nav_panel("Safety (Nausea)",
      plotOutput("er_safety_plot", height = "400px")
    ),
    nav_panel("Dose Comparison",
      plotOutput("dose_compare_plot", height = "400px")
    ),
    nav_panel("Model Information",
      markdown(paste0(
        "## Vortioxetine Population PK/PD Model\n\n",
        "**Source:** Naik H, Chan S, Vakilynejad M, et al. ",
        "A Population Pharmacokinetic-Pharmacodynamic Meta-Analysis of Vortioxetine ",
        "in Patients with Major Depressive Disorder. ",
        "*Basic Clin Pharmacol Toxicol.* 2016;118:344-355.\n\n",
        "### PK Model\n",
        "- **Structure:** 2-compartment with first-order absorption and linear elimination\n",
        "- **Data:** 10,498 concentrations from 3,160 MDD/GAD patients across 12 studies\n",
        "- **Covariates:** CrCL and height on CL/F; geographic region\n",
        "- **Software:** NONMEM 7 (FOCE-I)\n\n",
        "### PK/Efficacy Model\n",
        "- **Endpoint:** Change from baseline in MADRS total score (\u0394MADRS)\n",
        "- **Model:** Emax with baseline MADRS effect on EC50\n",
        "- **Data:** 2,537 patients from 7 short-term MDD studies\n",
        "- **EC50:** 24.9 ng/mL | **Emax:** 7.0 MADRS points\n",
        "- **Plateau:** Response plateaus at Cav > 19 ng/mL\n\n",
        "### PK/Safety Model\n",
        "- **Endpoint:** Treatment-emergent nausea (binary)\n",
        "- **Model:** Logistic regression with log(Cav)\n",
        "- **Predicted nausea risk:** 22-27% across 5-20 mg dose range\n\n",
        "### References\n",
        "- [FDA Label (Trintellix/Vortioxetine)](https://www.accessdata.fda.gov/drugsatfda_docs/label/2021/204447s021s022lbl.pdf)\n",
        "- [Naik et al. 2016 (PubMed)](https://doi.org/10.1111/bcpt.12513)\n",
        "- Brand name: Trintellix\u00AE (formerly Brintellix\u00AE)\n",
        "- Approved doses: 5, 10, 15, 20 mg once daily for MDD\n",
        "- Mechanism: Multimodal (SERT inhibitor + 5-HT receptor modulator)\n"
      ))
    )
  ),

  # ── 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"),
    " \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_days <- input$n_days
    region <- as.numeric(input$region)
    ht <- input$height
    crcl <- input$crcl

    out <- mod %>%
      param(REGION = region, HT = ht, CRCL = crcl) %>%
      ev(amt = dose_mg, cmt = 1, ii = 24, addl = n_days - 1) %>%
      mrgsim(end = n_days * 24, delta = 0.5) %>%
      as.data.frame()

    out
  })

  ss_metrics <- reactive({
    d <- sim_data()
    n_days <- input$n_days
    # Last dosing interval for steady-state
    ss_start <- (n_days - 1) * 24
    ss_end <- n_days * 24
    ss <- d %>% dplyr::filter(time >= ss_start & time <= ss_end)

    cmax <- max(ss$CP_ngml, na.rm = TRUE)
    ctrough <- min(ss$CP_ngml[ss$CP_ngml > 0], na.rm = TRUE)

    # Cav = Dose / (CL/F * 24)
    cli <- ss$CLi[1]
    dose_mg <- as.numeric(input$dose)
    cav <- (dose_mg / (cli * 24)) * 1e6  # mg -> ng conversion: mg/L * 1e6 = ng/mL... wait

    # Actually: Cav = Dose(mg) / (CL(L/hr) * tau(hr)) gives mg/L
    # mg/L * 1000 = ug/L = ng/mL? No. 1 mg/L = 1000 ug/L = 1000000 ng/L = 1000 ng/mL
    # So Cav(ng/mL) = Cav(mg/L) * 1000
    # But wait dose is in mg, CL in L/hr, tau=24hr
    # Cav(mg/L) = dose_mg / (CL_Lhr * 24)
    # Cav(ng/mL) = Cav(mg/L) * 1000
    # For 10mg, CL=42: 10/(42*24) = 0.00992 mg/L = 9.92 ng/mL... 
    # Paper says EC50=24.9 ng/mL which suggests Cav at 10mg ~10-20 ng/mL range
    # Actually let me recalc: paper says F=75%, these are apparent CLs
    # 10mg / (42 L/hr * 24hr) = 0.00992 mg/L = 9.92 ug/L... 
    # 1 mg/L = 1000 ug/L? No: 1 mg = 1000 ug, so 1 mg/L = 1000 ug/L
    # And 1 ug/L = 1 ng/mL. So 0.00992 mg/L = 9.92 ug/L = 9.92 ng/mL
    # That checks out with expected therapeutic range

    cav_ngml <- (dose_mg / (cli * 24)) * 1000.0

    # Half-life from 2-cpt
    k20 <- cli / 2900
    k23 <- 23 / 2900
    k32 <- 23 / 670
    beta_sum <- k20 + k23 + k32
    beta_prod <- k20 * k32
    disc <- beta_sum^2 - 4 * beta_prod
    disc <- max(disc, 0)
    lambda_z <- (beta_sum - sqrt(disc)) / 2
    thalf <- log(2) / lambda_z

    region_usa <- (as.numeric(input$region) == 1)
    dmadrs <- calc_dmadrs(cav_ngml, input$baseline_madrs, region_usa)
    nausea <- calc_nausea_risk(cav_ngml)

    list(cmax = cmax, ctrough = ctrough, cav = cav_ngml, thalf = thalf,
         dmadrs = dmadrs, nausea = nausea)
  })

  output$cmax_val <- renderText({
    sprintf("%.1f", ss_metrics()$cmax)
  })
  output$ctrough_val <- renderText({
    sprintf("%.1f", ss_metrics()$ctrough)
  })
  output$cav_val <- renderText({
    sprintf("%.1f", ss_metrics()$cav)
  })
  output$thalf_val <- renderText({
    sprintf("%.1f", ss_metrics()$thalf)
  })
  output$dmadrs_val <- renderText({
    sprintf("%.1f", ss_metrics()$dmadrs)
  })
  output$nausea_val <- renderText({
    sprintf("%.1f%%", ss_metrics()$nausea)
  })

  # ── PK Plot ──
  output$pk_plot <- renderPlot({
    d <- sim_data()
    if (input$log_scale) {
      d <- d %>% dplyr::filter(CP_ngml > 0.001)
    }
    p <- ggplot(d, aes(x = time / 24, y = CP_ngml)) +
      geom_line(color = "#8b5cf6", linewidth = 0.8) +
      labs(x = "Time (days)", y = "Vortioxetine Concentration (ng/mL)",
           title = paste0(input$dose, " mg QD — Plasma Concentration")) +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(face = "bold"))

    if (input$log_scale) p <- p + scale_y_log10()

    # Add EC50 reference line
    p <- p + geom_hline(yintercept = 24.9, linetype = "dashed", color = "#ef4444", alpha = 0.6) +
      annotate("text", x = 0, y = 26, label = "EC50 = 24.9 ng/mL",
               hjust = 0, color = "#ef4444", size = 3.5)

    p
  })

  # ── Efficacy E-R Plot ──
  output$er_efficacy_plot <- renderPlot({
    cav_range <- seq(0, 60, by = 0.5)
    baseline <- input$baseline_madrs
    region_usa <- (as.numeric(input$region) == 1)

    er_data <- data.frame(
      Cav = cav_range,
      dMADRS = sapply(cav_range, calc_dmadrs, baseline_madrs = baseline, region_usa = region_usa)
    )

    current_cav <- ss_metrics()$cav
    current_dm <- ss_metrics()$dmadrs

    ggplot(er_data, aes(x = Cav, y = dMADRS)) +
      geom_line(color = "#8b5cf6", linewidth = 1.2) +
      geom_point(data = data.frame(Cav = current_cav, dMADRS = current_dm),
                 aes(x = Cav, y = dMADRS), color = "#ef4444", size = 4) +
      geom_vline(xintercept = 24.9, linetype = "dashed", color = "#6c757d", alpha = 0.5) +
      annotate("text", x = 25.5, y = -10, label = "EC50", color = "#6c757d", size = 3.5) +
      labs(x = "Average Steady-State Concentration (ng/mL)",
           y = "\u0394MADRS (change from baseline)",
           title = "Exposure-Efficacy: Vortioxetine Cav vs \u0394MADRS",
           subtitle = paste0("Baseline MADRS = ", baseline,
                            " | Red dot = current dose (", input$dose, " mg)")) +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(face = "bold"))
  })

  # ── Safety E-R Plot ──
  output$er_safety_plot <- renderPlot({
    cav_range <- seq(1, 60, by = 0.5)
    safety_data <- data.frame(
      Cav = cav_range,
      Nausea = sapply(cav_range, calc_nausea_risk)
    )

    current_cav <- ss_metrics()$cav
    current_n <- ss_metrics()$nausea

    ggplot(safety_data, aes(x = Cav, y = Nausea)) +
      geom_line(color = "#f59e0b", linewidth = 1.2) +
      geom_point(data = data.frame(Cav = current_cav, Nausea = current_n),
                 aes(x = Cav, y = Nausea), color = "#ef4444", size = 4) +
      labs(x = "Average Steady-State Concentration (ng/mL)",
           y = "Predicted Nausea Risk (%)",
           title = "Exposure-Safety: Vortioxetine Cav vs Nausea Risk",
           subtitle = paste0("Red dot = current dose (", input$dose, " mg)")) +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(face = "bold")) +
      coord_cartesian(ylim = c(0, 50))
  })

  # ── Dose Comparison Plot ──
  output$dose_compare_plot <- renderPlot({
    region <- as.numeric(input$region)
    ht <- input$height
    crcl <- input$crcl
    doses <- c(5, 10, 15, 20)

    all_data <- lapply(doses, function(d_mg) {
      out <- mod %>%
        param(REGION = region, HT = ht, CRCL = crcl) %>%
        ev(amt = d_mg, cmt = 1, ii = 24, addl = 13) %>%
        mrgsim(end = 14 * 24, delta = 1) %>%
        as.data.frame() %>%
        mutate(Dose = paste0(d_mg, " mg"))
      out
    })
    combined <- bind_rows(all_data)
    combined$Dose <- factor(combined$Dose, levels = paste0(doses, " mg"))

    if (input$log_scale) {
      combined <- combined %>% dplyr::filter(CP_ngml > 0.001)
    }

    p <- ggplot(combined, aes(x = time / 24, y = CP_ngml, color = Dose)) +
      geom_line(linewidth = 0.8) +
      scale_color_manual(values = c("#10b981", "#8b5cf6", "#f59e0b", "#ef4444")) +
      geom_hline(yintercept = 24.9, linetype = "dashed", color = "#6c757d", alpha = 0.5) +
      annotate("text", x = 0, y = 26, label = "EC50", hjust = 0, color = "#6c757d", size = 3.5) +
      labs(x = "Time (days)", y = "Concentration (ng/mL)",
           title = "Dose Comparison: 5, 10, 15, 20 mg QD (2 weeks)",
           subtitle = "Dashed line = EC50 (24.9 ng/mL)") +
      theme_minimal(base_size = 14) +
      theme(plot.title = element_text(face = "bold"),
            legend.position = "top")

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

shinyApp(ui = ui, server = server)