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

model_code <- '
$PARAM @annotated
CL    : 0.0720 : Typical clearance (L/day) for 70 kg reference subject
VC    : 3.40   : Central volume (L) for 70 kg reference subject
Q     : 0.0900 : Intercompartmental clearance (L/day)
VP    : 1.60   : Peripheral volume (L) for 70 kg reference subject
WT    : 70     : Body weight (kg)
E0    : 1404   : Baseline LDH (U/L)
IMAX  : 0.825  : Maximum inhibition of LDH production
IC50  : 21     : Concentration for 50% of maximum LDH effect (mg/L ~= ug/mL)
HILL  : 1.52   : Hill coefficient
KOUTD : 0.2424 : LDH first-order degradation rate constant (1/day)

$CMT @annotated
CENT   : Central compartment (mg)
PERIPH : Peripheral compartment (mg)
LDH    : Serum LDH response (U/L)

$MAIN
double size_factor = WT / 70.0;
double CLi = CL * size_factor;
double VCi = VC * size_factor;
double Qi  = Q * size_factor;
double VPi = VP * size_factor;

LDH_0 = E0;

$ODE
double CP = CENT / VCi;
double CPpow = pow(CP > 0 ? CP : 0.0, HILL);
double IC50pow = pow(IC50, HILL);
double INH = IMAX * CPpow / (IC50pow + CPpow);
double KIN = E0 * KOUTD;

dxdt_CENT   = -(CLi / VCi) * CENT - (Qi / VCi) * CENT + (Qi / VPi) * PERIPH;
dxdt_PERIPH =  (Qi / VCi) * CENT - (Qi / VPi) * PERIPH;
dxdt_LDH    = KIN * (1.0 - INH) - KOUTD * LDH;

$TABLE
double CPtab = CENT / VCi;

$CAPTURE @annotated
CPtab : Serum ravulizumab concentration (mg/L ~= ug/mL)
CLi   : Individual clearance (L/day)
VCi   : Individual central volume (L)
VPi   : Individual peripheral volume (L)
'

mod <- mcode("ravulizumab_pkpd_322", model_code)

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; }
  .ref-box {
    background: #f0f4ff;
    border-left: 4px solid #8b5cf6;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
  .ref-box a { color: #8b5cf6; }
  .context-note {
    font-size: 12px;
    color: #6c757d;
    margin: 6px 0 0 6px;
  }
")

ui <- page_sidebar(
  title = "Ravulizumab (Ultomiris) PK/PD Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    accordion(
      open = c("Dosing", "Patient"),
      accordion_panel(
        "Dosing",
        icon = bsicons::bs_icon("capsule"),
        radioButtons(
          "dose_strategy",
          "Dosing strategy",
          choices = c(
            "FDA review / label weight-band regimen" = "label",
            "Custom maintenance regimen" = "custom"
          ),
          selected = "label"
        ),
        uiOutput("weight_band_text"),
        conditionalPanel(
          "input.dose_strategy === 'custom'",
          numericInput("loading_dose", "Loading dose (mg)", value = 2700, min = 1200, max = 5000, step = 100),
          numericInput("maintenance_dose", "Maintenance dose (mg)", value = 3300, min = 1200, max = 5000, step = 100)
        ),
        sliderInput("interval_days", "Maintenance interval (days)", min = 42, max = 70, value = 56, step = 1),
        numericInput("n_cycles", "Maintenance doses to simulate", value = 6, min = 1, max = 12),
        sliderInput("infusion_h", "Infusion duration (hours)", min = 1, max = 4, value = 2, step = 0.5)
      ),
      accordion_panel(
        "Patient",
        icon = bsicons::bs_icon("person"),
        sliderInput("wt", "Body weight (kg)", min = 40, max = 160, value = 70, step = 1),
        numericInput("baseline_ldh", "Baseline LDH (U/L)", value = 1404, min = 200, max = 4000, step = 25),
        tags$p(
          style = "font-size:12px; color:#6c757d; margin-top:8px;",
          "The FDA review identified body weight as a driver of both clearance and volume of distribution. This simulator uses linear size scaling around a 70 kg reference subject because the summary review did not publish the final covariate exponents."
        )
      ),
      accordion_panel(
        "Display",
        icon = bsicons::bs_icon("sliders"),
        checkboxInput("show_target", "Show therapeutic concentration target (146 ug/mL)", value = TRUE),
        checkboxInput("show_uln", "Show LDH normalization threshold (246 U/L)", value = TRUE),
        checkboxInput("log_scale", "Log scale concentration axis", value = FALSE)
      )
    )
  ),
  div(
    class = "context-note",
    "Default values reproduce the Phase 3 PNH dosing strategy summarized in FDA BLA 761108: loading dose on Day 1, maintenance starting Day 15, then every 8 weeks."
  ),
  layout_column_wrap(
    width = 1 / 4,
    fill = FALSE,
    div(
      class = "metric-card metric-success",
      div(class = "metric-value", textOutput("cmax")),
      div(class = "metric-label", "Cmax,ss (ug/mL)")
    ),
    div(
      class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("ctrough")),
      div(class = "metric-label", "Ctrough,ss (ug/mL)")
    ),
    div(
      class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("auc")),
      div(class = "metric-label", HTML("AUC<sub>tau,ss</sub> (ug.day/mL)"))
    ),
    div(
      class = "metric-card metric-info",
      div(class = "metric-value", textOutput("thalf")),
      div(class = "metric-label", "Terminal t1/2 (days)")
    )
  ),
  card(
    card_header("Serum Ravulizumab Concentration-Time Profile"),
    full_screen = TRUE,
    plotOutput("pkPlot", height = "500px")
  ),
  card(
    card_header("Predicted LDH Response"),
    full_screen = TRUE,
    plotOutput("ldhPlot", height = "420px")
  ),
  navset_card_underline(
    title = "Drug Information",
    nav_panel(
      "Model Information",
      markdown("
## Ravulizumab population PK/PD simulator

**Source document:** FDA BLA 761108 Multi-Disciplinary Review for Ultomiris (ravulizumab)  
**Population:** Adults with paroxysmal nocturnal hemoglobinuria (PNH)  
**Route:** Intravenous infusion  
**Target:** Complement component C5  

### Structural model used in this app
- Two-compartment linear PK model calibrated to the FDA review summary exposures
- Indirect response PD model for LDH lowering, using the final published PK-LDH parameters from the FDA review
- Weight-band dosing from the Phase 3 regimen: loading on Day 1, maintenance from Day 15, then every 8 weeks

### PK anchors from the FDA review
- Mean steady-state **Ctrough,ss**: **497 ug/mL**
- Mean steady-state **Cmax,ss**: **1460 ug/mL**
- Mean **Vss**: **5.35 L**
- Mean terminal **t1/2**: **49.7 days**

### PD parameters from Table 4 of the FDA review
| Parameter | Estimate | Units |
|-----------|----------|-------|
| E0 | 1404 | U/L |
| Imax | 0.825 | fraction |
| IC50 | 21.0 | ug/mL |
| IC95 | 146 | ug/mL |
| Hill coefficient | 1.52 | unitless |
| Kout | 0.0101 | h^-1 |

### Important modeling note
The FDA review provides the final PK-LDH model parameters for LDH, but it does **not** publish the full final population PK parameter table in the review text. To keep the simulator useful and honest, this app uses a calibrated 2-compartment PK approximation that prioritizes the reported Phase 3 concentration-time anchors and combines it with the published LDH indirect-response model.
")
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Primary references"),
        tags$ol(
          tags$li(
            tags$strong("FDA Center for Drug Evaluation and Research."),
            " BLA 761108 Multi-Disciplinary Review and Evaluation: Ultomiris (ravulizumab), 2018."
          ),
          tags$li(
            "ULTOMIRIS (ravulizumab-cwvz) prescribing information. Alexion Pharmaceuticals."
          )
        ),
        tags$h5("Therapeutic context"),
        tags$ul(
          tags$li(tags$strong("Class:"), " Complement C5 inhibitor monoclonal antibody"),
          tags$li(tags$strong("Indication in the review:"), " Adult PNH"),
          tags$li(tags$strong("Approved adult PNH regimen:"), " 2400/3000 mg, 2700/3300 mg, or 3000/3600 mg by weight band"),
          tags$li(tags$strong("Clinical PK message:"), " The Phase 3 regimen maintained trough concentrations above the therapeutic target through the full 8-week interval"),
          tags$li(tags$strong("Clinical PD message:"), " Ravulizumab produced immediate, complete, and sustained terminal complement inhibition with near-maximal LDH lowering")
        )
      )
    )
  ),
  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 <- function(input, output, session) {
  trap_auc <- function(time, conc) {
    if (length(time) < 2) {
      return(NA_real_)
    }
    sum(diff(time) * (head(conc, -1) + tail(conc, -1)) / 2)
  }

  approved_regimen <- reactive({
    shiny::req(input$wt)

    if (input$wt < 60) {
      list(group = "40 to <60 kg", loading = 2400, maintenance = 3000)
    } else if (input$wt < 100) {
      list(group = "60 to <100 kg", loading = 2700, maintenance = 3300)
    } else {
      list(group = ">=100 kg", loading = 3000, maintenance = 3600)
    }
  })

  regimen <- reactive({
    shiny::req(input$dose_strategy, input$interval_days, input$n_cycles, input$infusion_h)

    if (input$dose_strategy == "label") {
      band <- approved_regimen()
      list(
        label = paste0(band$loading, " mg load, then ", band$maintenance, " mg q", input$interval_days, "d"),
        loading = band$loading,
        maintenance = band$maintenance,
        group = band$group
      )
    } else {
      list(
        label = paste0(input$loading_dose, " mg load, then ", input$maintenance_dose, " mg q", input$interval_days, "d"),
        loading = input$loading_dose,
        maintenance = input$maintenance_dose,
        group = "Custom"
      )
    }
  })

  output$weight_band_text <- renderUI({
    band <- approved_regimen()
    tags$p(
      style = "font-size:12px; color:#6c757d; margin-top:6px;",
      paste0(
        "Current weight band: ", band$group,
        " (", band$loading, " mg loading / ", band$maintenance, " mg maintenance)"
      )
    )
  })

  sim_end_days <- reactive({
    shiny::req(input$interval_days, input$n_cycles)
    14 + input$interval_days * input$n_cycles
  })

  regimen_events <- reactive({
    shiny::req(input$interval_days, input$n_cycles, input$infusion_h)

    reg <- regimen()
    infusion_days <- input$infusion_h / 24
    maint_times <- 14 + seq(0, by = input$interval_days, length.out = input$n_cycles)
    event_df <- data.frame(
      time = c(0, maint_times),
      amt = c(reg$loading, rep(reg$maintenance, input$n_cycles)),
      cmt = 1,
      rate = c(reg$loading / infusion_days, rep(reg$maintenance / infusion_days, input$n_cycles))
    )

    mrgsolve::as.ev(event_df)
  })

  sim_data <- reactive({
    shiny::req(input$wt, input$baseline_ldh)

    mod %>%
      param(WT = input$wt, E0 = input$baseline_ldh) %>%
      ev(regimen_events()) %>%
      mrgsim(end = sim_end_days(), delta = 0.05) %>%
      as.data.frame() |>
      mutate(
        time_d = time,
        time_w = time / 7,
        CP = pmax(.data$CPtab, 0),
        LDH = pmax(.data$LDH, 0)
      )
  })

  last_interval <- reactive({
    shiny::req(input$interval_days, input$n_cycles)

    d <- sim_data()
    last_start <- 14 + (input$n_cycles - 1) * input$interval_days

    d |>
      dplyr::filter(.data$time_d >= last_start, .data$time_d <= last_start + input$interval_days)
  })

  individual_pk <- reactive({
    shiny::req(input$wt)

    size_factor <- input$wt / 70
    CLi <- 0.0720 * size_factor
    VCi <- 3.40 * size_factor
    Qi <- 0.0900 * size_factor
    VPi <- 1.60 * size_factor

    k10 <- CLi / VCi
    k12 <- Qi / VCi
    k21 <- Qi / VPi
    disc <- (k10 + k12 + k21)^2 - 4 * k10 * k21
    beta <- ((k10 + k12 + k21) - sqrt(disc)) / 2

    list(CLi = CLi, VCi = VCi, Qi = Qi, VPi = VPi, beta = beta)
  })

  output$cmax <- renderText({
    li <- last_interval()
    sprintf("%.0f", max(li$CP, na.rm = TRUE))
  })

  output$ctrough <- renderText({
    li <- last_interval()
    sprintf("%.0f", min(li$CP, na.rm = TRUE))
  })

  output$auc <- renderText({
    li <- last_interval()
    sprintf("%.0f", trap_auc(li$time_d, li$CP))
  })

  output$thalf <- renderText({
    pk <- individual_pk()
    sprintf("%.1f", log(2) / pk$beta)
  })

  output$pkPlot <- renderPlot({
    d <- sim_data()
    reg <- regimen()

    if (input$log_scale) {
      d <- d |>
        dplyr::filter(.data$CP > 0)
    }

    max_cp <- max(d$CP, na.rm = TRUE)
    target_top <- max(max_cp * 1.05, 146)

    p <- ggplot(d, aes(x = .data$time_d, y = .data$CP)) +
      geom_line(color = "#8b5cf6", linewidth = 0.9) +
      labs(
        x = "Time (days)",
        y = "Ravulizumab concentration (ug/mL)",
        title = paste0("Ravulizumab ", reg$label),
        subtitle = "Calibrated to FDA review Phase 3 exposure summary for adult PNH"
      ) +
      theme_minimal(base_size = 13)

    if (isTRUE(input$show_target)) {
      p <- p +
        annotate("rect", xmin = -Inf, xmax = Inf, ymin = 146, ymax = target_top, fill = "#10b981", alpha = 0.08) +
        geom_hline(yintercept = 146, linetype = "dashed", color = "#10b981", linewidth = 0.7) +
        annotate(
          "text",
          x = max(d$time_d) * 0.82,
          y = 146 * 1.08,
          label = "IC95 / therapeutic target ≈ 146 ug/mL",
          hjust = 0,
          size = 3.4,
          color = "#198754"
        )
    }

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

    p
  })

  output$ldhPlot <- renderPlot({
    d <- sim_data()
    reg <- regimen()

    max_ldh <- max(d$LDH, na.rm = TRUE)

    p <- ggplot(d, aes(x = .data$time_d, y = .data$LDH)) +
      geom_line(color = "#2563eb", linewidth = 0.9) +
      geom_hline(yintercept = input$baseline_ldh, color = "#6c757d", linetype = "dotted") +
      labs(
        x = "Time (days)",
        y = "Serum LDH (U/L)",
        title = paste0("Predicted LDH response: ", reg$label),
        subtitle = "Indirect response model using FDA Table 4 PK-LDH parameters"
      ) +
      theme_minimal(base_size = 13)

    if (isTRUE(input$show_uln)) {
      p <- p +
        annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = 246, fill = "#10b981", alpha = 0.10) +
        geom_hline(yintercept = 246, linetype = "dashed", color = "#10b981", linewidth = 0.7) +
        annotate(
          "text",
          x = max(d$time_d) * 0.80,
          y = min(260, max_ldh * 0.35),
          label = "LDH normalization threshold = 246 U/L",
          hjust = 0,
          size = 3.4,
          color = "#198754"
        )
    }

    p
  })
}

shinyApp(ui = ui, server = server)