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

# ── Bedaquiline + M2 PopPK Model ──
# Svensson et al. CPT:PSP 2016; 5:682-691
# 3-compartment BDQ + 1-compartment M2
# Covariates: weight (allometric), albumin (protein binding), age, race

model_code <- '
$PARAM @annotated
CL    : 2.62   : BDQ apparent clearance (L/h)
VC    : 198    : BDQ central volume (L) at 70 kg
Q1    : 3.66   : BDQ intercompartmental CL 1 (L/h)
VP1   : 8550   : BDQ peripheral volume 1 (L)
Q2    : 7.34   : BDQ intercompartmental CL 2 (L/h)
VP2   : 2690   : BDQ peripheral volume 2 (L)
KA    : 1.1    : Absorption rate constant (1/h)
CLM2  : 10.0   : M2 apparent clearance (L/h)
VM2   : 2200   : M2 volume (L)
FM    : 0.5    : Fraction metabolized to M2 (assumed)
WT    : 70     : Body weight (kg)
ALB   : 4.04   : Albumin (g/dL)
AGE   : 32     : Age (years)
BLACK : 0      : Black race (0/1)
ALB_SS : 4.04  : Albumin at steady state (g/dL)

$CMT @annotated
DEPOT  : Oral depot (mg)
CENT   : BDQ central (mg)
P1     : BDQ peripheral 1 (mg)
P2     : BDQ peripheral 2 (mg)
M2CPT  : M2 compartment (mg)

$MAIN
// Allometric scaling
double CLi  = CL * pow(WT / 70.0, 0.18);
double VCi  = VC * (WT / 70.0);
double Q1i  = Q1 * pow(WT / 70.0, 0.18);
double VP1i = VP1 * (WT / 70.0);
double Q2i  = Q2 * pow(WT / 70.0, 0.18);
double VP2i = VP2 * (WT / 70.0);
double CLM2i = CLM2 * pow(WT / 70.0, 0.18);
double VM2i  = VM2 * (WT / 70.0);

// Albumin effect on protein binding (power model coeff = 1)
double alb_ratio = ALB / ALB_SS;
double fu_ratio = 1.0 / alb_ratio;  // lower albumin -> higher fu -> higher apparent params
CLi  = CLi * fu_ratio;
VCi  = VCi * fu_ratio;
Q1i  = Q1i * fu_ratio;
VP1i = VP1i * fu_ratio;
Q2i  = Q2i * fu_ratio;
VP2i = VP2i * fu_ratio;
CLM2i = CLM2i * fu_ratio;
VM2i  = VM2i * fu_ratio;

// Age effect: 0.88% decrease per year from median 32
CLi  = CLi * (1.0 - 0.0088 * (AGE - 32.0));
CLM2i = CLM2i * (1.0 - 0.0088 * (AGE - 32.0));

// Black race: 84% higher CL
if (BLACK == 1) {
  CLi  = CLi * 1.84;
  CLM2i = CLM2i * 1.84;
}

// Floor values
CLi   = CLi > 0.01 ? CLi : 0.01;
CLM2i = CLM2i > 0.01 ? CLM2i : 0.01;
VCi   = VCi > 1.0 ? VCi : 1.0;
VM2i  = VM2i > 1.0 ? VM2i : 1.0;

$ODE
double CP_BDQ = CENT / VCi;
dxdt_DEPOT = -KA * DEPOT;
dxdt_CENT  = KA * DEPOT - (CLi/VCi + Q1i/VCi + Q2i/VCi) * CENT + (Q1i/VP1i) * P1 + (Q2i/VP2i) * P2;
dxdt_P1    = (Q1i/VCi) * CENT - (Q1i/VP1i) * P1;
dxdt_P2    = (Q2i/VCi) * CENT - (Q2i/VP2i) * P2;
dxdt_M2CPT = FM * CLi * CP_BDQ - (CLM2i / VM2i) * M2CPT;

$TABLE
double BDQ_conc = CENT / VCi;
BDQ_conc = BDQ_conc > 0 ? BDQ_conc : 0;
double M2_conc = M2CPT / VM2i;
M2_conc = M2_conc > 0 ? M2_conc : 0;
// Convert mg/L to ug/mL (same units, 1:1)

$CAPTURE @annotated
BDQ_conc : Bedaquiline concentration (mg/L)
M2_conc  : M2 metabolite concentration (mg/L)
'

mod <- mcode("bedaquiline_m2", 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; }
")

# ── Helper: build dosing regimen ──
build_regimen <- function(dose_load, dose_maint, duration_weeks) {
  events <- list()
  # Loading phase: 400 mg QD x 14 days
  events[[1]] <- ev(amt = dose_load, cmt = 1, ii = 24, addl = 13, time = 0)
  # Maintenance: 200 mg TIW (Mon/Wed/Fri = every 56h avg, use ii=56)
  # More accurately: doses at 48, 48, 72h intervals repeating
  # Simplified: 3x/week for remaining weeks
  maint_weeks <- duration_weeks - 2
  if (maint_weeks > 0) {
    maint_doses <- maint_weeks * 3 - 1
    # TIW: Mon(0h), Wed(48h), Fri(96h), Mon(168h)...
    # Pattern: 48, 48, 72 repeating
    t_start <- 14 * 24  # after loading
    times <- c()
    t <- t_start
    for (i in seq_len(maint_doses + 1)) {
      times <- c(times, t)
      cycle_pos <- (i - 1) %% 3
      if (cycle_pos == 0) t <- t + 48
      else if (cycle_pos == 1) t <- t + 48
      else t <- t + 72
    }
    for (tt in times) {
      events[[length(events) + 1]] <- ev(amt = dose_maint, cmt = 1, time = tt)
    }
  }
  do.call(c, events)
}

# ── UI ──
ui <- page_sidebar(
  title = "Bedaquiline + M2 PopPK Simulator",
  theme = app_theme,

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

    h6("Patient Characteristics"),
    sliderInput("weight", "Body Weight (kg)", min = 30, max = 120, value = 57, step = 1),
    sliderInput("albumin", "Albumin (g/dL)", min = 1.5, max = 5.0, value = 3.65, step = 0.05),
    sliderInput("age", "Age (years)", min = 18, max = 70, value = 32, step = 1),
    checkboxInput("black_race", "Black Race", value = FALSE),

    hr(),
    h6("Dosing Regimen"),
    numericInput("dose_load", "Loading Dose (mg)", value = 400, min = 100, max = 800, step = 50),
    numericInput("dose_maint", "Maintenance Dose (mg)", value = 200, min = 50, max = 600, step = 50),
    sliderInput("duration", "Treatment Duration (weeks)", min = 4, max = 36, value = 24, step = 1),
    sliderInput("sim_duration", "Simulation Duration (weeks)", min = 8, max = 72, value = 48, step = 4),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE),

    hr(),
    card(
      card_header("Model Parameters", class = "bg-light"),
      p(strong("CL/F:"), "2.62 L/h"),
      p(strong("V/F:"), "198 L (70 kg)"),
      p(strong("VP1/F:"), "8,550 L"),
      p(strong("VP2/F:"), "2,690 L"),
      p(strong("t½ terminal:"), "> 5 months"),
      p(em("Ref: Svensson et al. CPT:PSP 2016"))
    )
  ),

  # ── Metrics Row ──
  layout_column_wrap(
        width = 1/4,
        fill = FALSE,
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("cmax_bdq")),
      div(class = "metric-label", "BDQ Cmax (mg/L)")
    ),
    div(class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("ctrough_bdq")),
      div(class = "metric-label", "BDQ Ctrough Wk24 (mg/L)")
    ),
    div(class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("cmax_m2")),
      div(class = "metric-label", "M2 Cmax (mg/L)")
    ),
    div(class = "metric-card metric-info",
      div(class = "metric-value", textOutput("m2_bdq_ratio")),
      div(class = "metric-label", "M2:BDQ Ratio at Wk24")
    )
  ),

  # ── Main PK Plot ──
  card(
    card_header("Bedaquiline & M2 Concentration-Time Profiles"),
    full_screen = TRUE,
    plotOutput("pkPlot", height = "500px")
  ),

  # ── Tabbed Content ──
  navset_card_tab(
    title = "Analysis",
    nav_panel("Covariate Effects",
      plotOutput("covPlot", height = "400px")
    ),
    nav_panel("Dosing Schedule",
      plotOutput("dosePlot", height = "400px")
    ),
    nav_panel("Model Information", 
      markdown(paste0(
        "## Bedaquiline + M2 Population PK Model\n\n",
        "**Source:** Svensson EM, Dosne A-G, Karlsson MO. *Population Pharmacokinetics of ",
        "Bedaquiline and Metabolite M2 in Patients With Drug-Resistant Tuberculosis.* ",
        "CPT Pharmacometrics Syst Pharmacol. 2016;5(12):682-691.\n\n",
        "### Model Structure\n",
        "- **Bedaquiline:** 3-compartment model with transit absorption\n",
        "- **M2 Metabolite:** 1-compartment model\n",
        "- **Absorption:** First-order with ka ≈ 1.1 h⁻¹\n\n",
        "### Key Features\n",
        "- Extremely long terminal half-life (>5 months)\n",
        "- Loading phase (400 mg QD × 2 weeks) + maintenance (200 mg TIW × 22 weeks)\n",
        "- Highly protein bound (>99.9% BDQ, >99.7% M2)\n",
        "- Albumin affects disposition via protein binding changes\n",
        "- Black race: 84% higher apparent CL\n",
        "- Age: 0.88% CL decrease per year\n\n",
        "### Clinical Context\n",
        "- Approved for MDR-TB treatment\n",
        "- M2 metabolite associated with QT prolongation\n",
        "- Concentrations highest at end of loading phase (week 2)\n",
        "- CYP3A4-mediated metabolism\n\n",
        "### References\n",
        "- [FDA Label (SIRTURO)](https://www.accessdata.fda.gov/drugsatfda_docs/label/2024/204384s019lbl.pdf)\n",
        "- [McLeay et al. AAC 2014 - Prior PopPK](https://pmc.ncbi.nlm.nih.gov/articles/PMC4135833/)\n",
        "- [Frontiers 2023 - Dosing Evaluation](https://www.frontiersin.org/journals/pharmacology/articles/10.3389/fphar.2023.1022090/full)\n\n",
        "### Disclaimer\n",
        "*For research and educational purposes only. Not for clinical decision-making.*"
      ))
    )
  ),

  # ── 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({
    shiny::req(input$weight, input$albumin, input$age)

    ev_data <- build_regimen(input$dose_load, input$dose_maint, input$duration)
    end_h <- input$sim_duration * 7 * 24

    out <- mod %>%
      param(
        WT = input$weight,
        ALB = input$albumin,
        AGE = input$age,
        BLACK = as.numeric(input$black_race)
      ) %>%
      ev(ev_data) %>%
      mrgsim(end = end_h, delta = 6) %>%
      as.data.frame()

    out$time_weeks <- out$time / (24 * 7)
    out
  })

  # ── Metrics ──
  output$cmax_bdq <- renderText({
    d <- sim_data()
    sprintf("%.2f", max(d$BDQ_conc, na.rm = TRUE))
  })

  output$ctrough_bdq <- renderText({
    d <- sim_data()
    wk24_h <- 24 * 7 * 24
    idx <- which.min(abs(d$time - wk24_h))
    sprintf("%.3f", d$BDQ_conc[idx])
  })

  output$cmax_m2 <- renderText({
    d <- sim_data()
    sprintf("%.3f", max(d$M2_conc, na.rm = TRUE))
  })

  output$m2_bdq_ratio <- renderText({
    d <- sim_data()
    wk24_h <- 24 * 7 * 24
    idx <- which.min(abs(d$time - wk24_h))
    bdq <- d$BDQ_conc[idx]
    m2 <- d$M2_conc[idx]
    if (bdq > 0.001) sprintf("%.2f", m2 / bdq) else "N/A"
  })

  # ── Main PK Plot ──
  output$pkPlot <- renderPlot({
    d <- sim_data()

    d_long <- rbind(
      data.frame(time_weeks = d$time_weeks, Concentration = d$BDQ_conc, Analyte = "Bedaquiline"),
      data.frame(time_weeks = d$time_weeks, Concentration = d$M2_conc, Analyte = "M2 Metabolite")
    )

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

    p <- ggplot(d_long, aes(x = time_weeks, y = Concentration, color = Analyte)) +
      geom_line(linewidth = 1) +
      scale_color_manual(values = c("Bedaquiline" = "#8b5cf6", "M2 Metabolite" = "#f59e0b")) +
      geom_vline(xintercept = 2, linetype = "dashed", color = "grey50", alpha = 0.7) +
      geom_vline(xintercept = input$duration, linetype = "dashed", color = "red", alpha = 0.5) +
      annotate("text", x = 2, y = Inf, label = "End Loading", vjust = 1.5, hjust = -0.1, size = 3, color = "grey40") +
      annotate("text", x = input$duration, y = Inf, label = "End Treatment", vjust = 1.5, hjust = -0.1, size = 3, color = "red") +
      labs(x = "Time (weeks)", y = "Concentration (mg/L)", color = NULL) +
      theme_minimal(base_size = 14) +
      theme(legend.position = "top")

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

  # ── Covariate Effects Plot ──
  output$covPlot <- renderPlot({
    weights <- c(40, 55, 70, 85, 100)
    ev_data <- build_regimen(400, 200, 24)
    end_h <- 24 * 7 * 24

    results <- lapply(weights, function(w) {
      out <- mod %>%
        param(WT = w, ALB = 4.04, AGE = 32, BLACK = 0) %>%
        ev(ev_data) %>%
        mrgsim(end = end_h, delta = 12) %>%
        as.data.frame()
      out$Weight <- paste0(w, " kg")
      out$time_weeks <- out$time / (24 * 7)
      out
    })
    d <- do.call(rbind, results)

    ggplot(d, aes(x = time_weeks, y = BDQ_conc, color = Weight)) +
      geom_line(linewidth = 0.8) +
      scale_color_brewer(palette = "Set1") +
      labs(x = "Time (weeks)", y = "BDQ Concentration (mg/L)",
           title = "Effect of Body Weight on Bedaquiline Exposure",
           color = "Body Weight") +
      theme_minimal(base_size = 14) +
      theme(legend.position = "top")
  })

  # ── Dosing Schedule Plot ──
  output$dosePlot <- renderPlot({
    d <- sim_data()
    # Show dose amounts over time
    dose_times <- d %>% dplyr::filter(BDQ_conc > 0)

    ggplot(d, aes(x = time_weeks, y = BDQ_conc)) +
      geom_area(fill = "#8b5cf6", alpha = 0.3) +
      geom_line(color = "#8b5cf6", linewidth = 0.8) +
      annotate("rect", xmin = 0, xmax = 2, ymin = -Inf, ymax = Inf,
               fill = "#ef4444", alpha = 0.08) +
      annotate("rect", xmin = 2, xmax = input$duration, ymin = -Inf, ymax = Inf,
               fill = "#3b82f6", alpha = 0.05) +
      annotate("text", x = 1, y = max(d$BDQ_conc) * 0.95,
               label = paste0("Loading\n", input$dose_load, " mg QD"), size = 3.5, color = "#ef4444") +
      annotate("text", x = (2 + input$duration) / 2, y = max(d$BDQ_conc) * 0.95,
               label = paste0("Maintenance\n", input$dose_maint, " mg TIW"), size = 3.5, color = "#3b82f6") +
      labs(x = "Time (weeks)", y = "BDQ Concentration (mg/L)",
           title = "Bedaquiline Dosing Phases") +
      theme_minimal(base_size = 14)
  })
}

shinyApp(ui = ui, server = server)