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

model_code <- '
$PARAM @annotated
CL     : 0.241 : Effective clearance (L/day)
V      : 8.0   : Effective distribution volume (L)
KA     : 1.04  : Absorption rate constant (1/day)
F1     : 0.70  : Absolute bioavailability (fraction)
WT     : 70    : Body weight (kg)
REFWT  : 70    : Reference body weight (kg)
EXPCL  : 0.75  : Allometric exponent on clearance
EXPV   : 1.00  : Allometric exponent on volume
INDCL  : 1.00  : Indication-specific clearance multiplier
INDV   : 1.00  : Indication-specific volume multiplier

$CMT @annotated
DEPOT : Subcutaneous depot (mg)
CENT  : Central compartment (mg)

$MAIN

double CLi = CL * pow(WT / REFWT, EXPCL) * INDCL;
double Vi = V * pow(WT / REFWT, EXPV) * INDV;
F_DEPOT = F1;

$ODE

dxdt_DEPOT = -KA * DEPOT;
dxdt_CENT = KA * DEPOT - (CLi / Vi) * CENT;

$TABLE

double CP = CENT / Vi;
double CL_OUT = CLi;
double V_OUT = Vi;
double F_SC = F1;

$CAPTURE @annotated
CP    : Plasma concentration (mg/L; numerically equal to ug/mL)
CL_OUT : Adjusted clearance (L/day)
V_OUT : Adjusted effective distribution volume (L)
F_SC  : Subcutaneous bioavailability
'

mod <- mcode("bimekizumab_316", model_code, quiet = TRUE)

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

get_indication_factors <- function(indication) {
  switch(
    indication,
    hs = list(cl = 1.31, v = 1.18, label = "HS exposure adjustment: CL +31%, V +18%"),
    list(cl = 1.00, v = 1.00, label = "Reference adult population PK")
  )
}

build_event <- function(dose_mg, interval_days, n_doses) {
  ev(
    amt = dose_mg,
    cmt = 1,
    ii = interval_days,
    addl = max(n_doses - 1, 0)
  )
}

auc_trap <- function(time, conc) {
  if (length(time) < 2 || length(conc) < 2) return(NA_real_)
  sum(diff(time) * (head(conc, -1) + tail(conc, -1)) / 2)
}

ui <- page_sidebar(
  title = "Bimekizumab (BIMZELX) PK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    h6("Clinical context"),
    selectInput(
      "indication",
      "Indication group",
      choices = c(
        "Psoriasis / PsA / nr-axSpA / AS" = "pso",
        "Hidradenitis suppurativa" = "hs"
      ),
      selected = "pso"
    ),
    sliderInput("dose", "Dose (mg)", min = 160, max = 320, value = 320, step = 160),
    sliderInput("interval_w", "Dosing interval (weeks)", min = 2, max = 8, value = 4, step = 2),
    sliderInput("n_doses", "Number of doses", min = 1, max = 16, value = 1, step = 1),
    hr(),
    h6("Patient"),
    sliderInput("wt", "Body weight (kg)", min = 40, max = 140, value = 70, step = 5),
    hr(),
    h6("Display"),
    checkboxInput("show_ref_band", "Show observed peak range (12–50 µg/mL)", value = TRUE),
    checkboxInput("log_scale", "Log scale (Y-axis)", value = FALSE)
  ),
  layout_column_wrap(
    width = 1 / 4,
    fill = FALSE,
    div(class = "metric-card metric-success",
      div(class = "metric-value", textOutput("cmax")),
      div(class = "metric-label", "Last-interval Cmax (µg/mL)")
    ),
    div(class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("ctrough")),
      div(class = "metric-label", "Last-interval Ctrough (µg/mL)")
    ),
    div(class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("auc")),
      div(class = "metric-label", "Last-interval AUC (µg·day/mL)")
    ),
    div(class = "metric-card metric-info",
      div(class = "metric-value", textOutput("thalf")),
      div(class = "metric-label", "Approx. t½ (days)")
    )
  ),
  card(
    card_header(textOutput("plot_title")),
    full_screen = TRUE,
    plotOutput("pkPlot", height = "500px")
  ),
  navset_card_underline(
    title = "Drug Information",
    nav_panel(
      "Model Information",
      markdown(
        "
## Bimekizumab population PK summary

**Drug:** Bimekizumab-bkzx (BIMZELX)  
**Class:** Humanized IgG1 monoclonal antibody against IL-17A / IL-17F  
**Route:** Subcutaneous  
**Model used here:** Pragmatic one-compartment SC biologic model with first-order absorption and body-weight scaling

### PK anchors from the prescribing information
- Median peak concentration: **25 µg/mL**
- Median Tmax: **3–4 days**
- Absolute bioavailability: **70%**
- Apparent clearance: **0.337 L/day**
- Terminal half-life: **23 days**
- Volume of distribution at steady state: **11.2 L**

### How this simulator is calibrated
- The model is calibrated to reproduce the label PK anchors for SC dosing: **Cmax ~25 µg/mL**, **Tmax ~3.5 days**, **half-life ~23 days**.
- Weight is implemented with standard biologic-style allometric scaling on clearance and distribution volume.
- Hidradenitis suppurativa applies the label-reported exposure shift: **CL +31%**, **V +18%**.
- The effective model volume is smaller than the label Vss because this is a practical exposure simulator rather than a full proprietary PopPK reconstruction.

### Practical dosing notes
- **Plaque psoriasis:** 320 mg Q4W induction, then Q8W maintenance
- **Psoriatic arthritis / axial spondyloarthritis:** 160 mg Q4W
- **Hidradenitis suppurativa:** 320 mg Q2W induction, then Q4W maintenance

This is a pragmatic label-based simulator, not a full proprietary PopPK reconstruction.
        "
      )
    ),
    nav_panel(
      "References",
      div(class = "ref-box",
        tags$h5("Key sources"),
        tags$ol(
          tags$li(
            "BIMZELX (bimekizumab-bkzx) US Prescribing Information. FDA label 761151s010."
          ),
          tags$li(
            "Bimzelx Summary of Product Characteristics (SmPC), medicines.org.uk — confirms class and population PK summary."
          )
        ),
        tags$h5("Model context"),
        tags$ul(
          tags$li(tags$strong("Drug class:"), " Interleukin inhibitor / immunosuppressant"),
          tags$li(tags$strong("Analyte:"), " Plasma bimekizumab concentration"),
          tags$li(tags$strong("Units:"), " mg/L, numerically equivalent to µg/mL"),
          tags$li(tags$strong("Use case:"), " Exposure visualization across approved SC regimens and body weights")
        )
      )
    )
  ),
  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) {
  sim_data <- reactive({
    factors <- get_indication_factors(input$indication)
    interval_days <- input$interval_w * 7
    end_days <- input$n_doses * interval_days
    ev1 <- build_event(input$dose, interval_days, input$n_doses)

    mod %>%
      param(
        WT = input$wt,
        INDCL = factors$cl,
        INDV = factors$v
      ) %>%
      ev(ev1) %>%
      mrgsim(end = end_days, delta = 0.25) %>%
      as.data.frame() |>
      mutate(
        time_days = time,
        conc_ugml = CP
      )
  })

  summary_interval <- reactive({
    d <- sim_data()
    shiny::req(nrow(d) > 1)

    interval_days <- input$interval_w * 7
    metric_start <- max(max(d$time_days) - interval_days, 0)
    d_last <- d |>
      dplyr::filter(time_days >= metric_start)
    shiny::req(nrow(d_last) > 1)

    cmax <- max(d_last$conc_ugml, na.rm = TRUE)
    ctrough <- d_last$conc_ugml[nrow(d_last)]
    auc_tau <- auc_trap(d_last$time_days, d_last$conc_ugml)

    factors <- get_indication_factors(input$indication)
    cl_now <- 0.241 * (input$wt / 70) ^ 0.75 * factors$cl
    v_now <- 8.0 * (input$wt / 70) * factors$v
    thalf <- log(2) * v_now / cl_now

    list(
      cmax = cmax,
      ctrough = ctrough,
      auc_tau = auc_tau,
      thalf = thalf,
      metric_start = metric_start,
      label = factors$label
    )
  })

  output$cmax <- renderText({
    s <- summary_interval()
    sprintf("%.1f", s$cmax)
  })

  output$ctrough <- renderText({
    s <- summary_interval()
    sprintf("%.1f", s$ctrough)
  })

  output$auc <- renderText({
    s <- summary_interval()
    sprintf("%.0f", s$auc_tau)
  })

  output$thalf <- renderText({
    s <- summary_interval()
    sprintf("%.1f", s$thalf)
  })

  output$plot_title <- renderText({
    indication_label <- if (input$indication == "hs") "Hidradenitis suppurativa" else "Psoriasis / PsA / axSpA"
    paste0(
      "Concentration-Time Profile — ",
      input$dose,
      " mg Q",
      input$interval_w,
      "W • ",
      indication_label
    )
  })

  output$pkPlot <- renderPlot({
    d <- sim_data()
    s <- summary_interval()
    d_plot <- if (input$log_scale) d |> dplyr::filter(conc_ugml > 0.001) else d
    shiny::req(nrow(d_plot) > 0)

    p <- ggplot(d_plot, aes(x = time_days, y = conc_ugml))

    if (isTRUE(input$show_ref_band)) {
      p <- p +
        annotate(
          "rect",
          xmin = -Inf,
          xmax = Inf,
          ymin = 12,
          ymax = 50,
          fill = "#8b5cf6",
          alpha = 0.08
        ) +
        geom_hline(
          yintercept = c(12, 25, 50),
          linetype = c("dashed", "solid", "dashed"),
          color = c("#8b5cf6", "#6d28d9", "#8b5cf6"),
          alpha = c(0.4, 0.6, 0.4)
        )
    }

    p <- p +
      geom_line(color = "#8b5cf6", linewidth = 0.9) +
      labs(
        x = "Time (days)",
        y = "Plasma concentration (µg/mL)",
        subtitle = paste0("Metrics reflect the most recent dosing interval • ", s$label)
      ) +
      theme_minimal(base_size = 14)

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

    p
  })
}

shinyApp(ui = ui, server = server)