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

model_code <- '
$PARAM @annotated
CL      : 11.4   : Apparent clearance (L/h)
V       : 476    : Apparent volume (L)
KA      : 2.0    : First-order absorption (1/h)
F1      : 1.0    : Relative bioavailability
TLAG    : 0.60   : Absorption lag time (h)
WT      : 70     : Body weight (kg)
AGECL   : 1.0    : Age-group multiplier on clearance
KPCSF   : 0.025  : CSF/plasma scaling factor
IMAX    : 0.56   : Maximum fractional inhibition of sTREM2 production
IC50    : 0.70   : CSF concentration for 50% of Imax (ng/mL)
KOUT    : 0.016  : sTREM2 turnover rate (1/h)

$CMT @annotated
DEPOT   : Oral depot (mg-equivalent)
CENT    : Central compartment (mg-equivalent)
STREM2  : Soluble TREM2 biomarker state (% baseline)

$MAIN
double CLi = CL * pow(WT / 70.0, 0.75) * AGECL;
double Vi = V * (WT / 70.0);
double KIN = KOUT * 100.0;

F_DEPOT = F1;
ALAG_DEPOT = TLAG;
STREM2_0 = 100.0;

$ODE
double CPi = (CENT / Vi) * 1000.0;
double CCSFi = CPi * KPCSF;
double INHi = IMAX * CCSFi / (IC50 + CCSFi);

dxdt_DEPOT = -KA * DEPOT;
dxdt_CENT = KA * DEPOT - (CLi / Vi) * CENT;
dxdt_STREM2 = KIN * (1.0 - INHi) - KOUT * STREM2;

$TABLE
double CP = (CENT / Vi) * 1000.0;
double CCSF = CP * KPCSF;
double STREM2_VAL = STREM2;
double STREM2_PCT = STREM2 - 100.0;
double INH = IMAX * CCSF / (IC50 + CCSF);

$CAPTURE CP CCSF STREM2_VAL STREM2_PCT INH
'

mod <- mcode("vg3927_pkpd", model_code)

slide_pk_anchors <- data.frame(
  dose = c(6, 12, 25, 25, 50),
  age_group = c("Adult", "Adult", "Adult", "Elderly", "Adult"),
  time_h = c(339, 339, 339, 339, 339),
  cp = c(20, 45, 110, 100, 190)
)

slide_pd_anchors <- data.frame(
  label = c("6 mg", "12 mg", "25 mg", "25 mg elderly", "50 mg"),
  ccsf = c(0.35, 0.70, 2.00, 2.20, 3.50),
  effect = c(-22, -29, -43, -46, -44)
)

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; }
  .small-note {
    color: #6c757d;
    font-size: 12px;
    line-height: 1.5;
  }
")

trap_auc <- function(time, conc) {
  if (length(time) < 2) {
    return(NA_real_)
  }

  sum(diff(time) * (head(conc, -1) + tail(conc, -1)) / 2)
}

simulate_vg3927 <- function(
  dose = 25,
  interval = 24,
  n_days = 14,
  wt = 70,
  age_group = "Adult",
  csf_multiplier = 1,
  pd_multiplier = 1,
  delta = 0.25
) {
  n_doses <- max(1, ceiling(n_days * 24 / interval))
  dose_equivalent <- dose^0.95 * 25^0.05
  ev1 <- ev(
    amt = dose_equivalent,
    ii = interval,
    addl = n_doses - 1,
    cmt = 1
  )

  mod |>
    param(
      WT = wt,
      AGECL = 1.0,
      KPCSF = 0.025 * csf_multiplier,
      IC50 = 0.70 / pd_multiplier
    ) %>%
    ev(ev1) %>%
    mrgsim(end = n_days * 24, delta = delta) |>
    as.data.frame() |>
    mutate(
      time_h = time,
      time_day = time / 24,
      strem2_model_pct = STREM2_PCT,
      STREM2_PCT = pmax(STREM2_PCT, -44),
      strem2_reduction = -STREM2_PCT
    )
}

ui <- page_sidebar(
  title = "VG-3927 PK/PD Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    h6("Regimen"),
    sliderInput("dose", "Dose (mg)", min = 2, max = 140, value = 25, step = 1),
    sliderInput("interval", "Dosing interval (h)", min = 12, max = 48, value = 24, step = 12),
    numericInput("n_days", "Duration (days)", value = 14, min = 1, max = 28),
    hr(),
    h6("Patient factors"),
    sliderInput("wt", "Body weight (kg)", min = 45, max = 110, value = 70, step = 5),
    radioButtons("age_group", "Age cohort", choices = c("Adult", "Elderly"), selected = "Adult", inline = TRUE),
    hr(),
    h6("PK/PD assumptions"),
    sliderInput("csf_multiplier", "CSF penetration multiplier", min = 0.5, max = 1.5, value = 1.0, step = 0.05),
    sliderInput("pd_multiplier", "PD sensitivity multiplier", min = 0.7, max = 1.3, value = 1.0, step = 0.05),
    selectInput(
      "view",
      "Plot",
      choices = c("Plasma PK", "CSF Exposure", "sTREM2 PD", "Exposure-Response"),
      selected = "Plasma PK"
    ),
    checkboxInput("log_scale", "Log scale for concentration plots", value = FALSE),
    div(
      class = "small-note",
      "Calibrated to the AAIC 2025 Phase 1 slide deck. The source does not publish a full CL/F or V/F table, so this app uses a transparent compartmental approximation anchored to the visible SAD, MAD, and CSF sTREM2 curves."
    )
  ),
  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 (ng/mL)")
    ),
    div(
      class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("ctrough")),
      div(class = "metric-label", "Ctrough,ss (ng/mL)")
    ),
    div(
      class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("auc")),
      div(class = "metric-label", "AUC0-tau,ss (ng*h/mL)")
    ),
    div(
      class = "metric-card metric-info",
      div(class = "metric-value", textOutput("thalf")),
      div(class = "metric-label", "Apparent t1/2 (h)")
    )
  ),
  card(
    full_screen = TRUE,
    navset_card_underline(
      title = "VG-3927 Simulation",
      nav_panel(
        "Simulation",
        plotOutput("main_plot", height = "500px"),
        div(class = "small-note", textOutput("simulation_note"))
      ),
      nav_panel(
        "Model Information",
        markdown("
## VG-3927 phase 1 oral PK/PD approximation

**Drug/program:** VG-3927, an oral small-molecule TREM2 agonist for Alzheimer's disease.  
**Source:** Vigil Neuroscience phase 1 clinical presentation, AAIC 2025.  
**Implemented scope:** Oral plasma PK + CSF-linked sTREM2 biomarker suppression.

### Structural model
- **PK:** 1-compartment oral model with first-order absorption and short absorption lag
- **Covariates:** allometric weight scaling on CL/F and V/F; no clinically meaningful adult-vs-elderly adjustment was needed to reproduce the sponsor slide summary
- **CSF link:** CSF concentration approximated as a fixed fraction of plasma concentration
- **PD:** indirect-response Imax model where CSF VG-3927 inhibits sTREM2 production

### Calibrated anchor parameters
| Parameter | Value | Units | Comment |
|-----------|-------|-------|---------|
| CL/F | 11.4 | L/h | Calibrated to MAD concentration-time anchors |
| V/F | 476 | L | Compartmental approximation, not sponsor-reported |
| Ka | 2.0 | h^-1 | Supports observed 2-4 h Tmax |
| Tlag | 0.60 | h | Short oral lag |
| t1/2 | 28.9 | h | Consistent with QD accumulation and flat peak-to-trough |
| CSF scalar | 0.025 | ratio | Produces the visible day-15 CSF exposure range |
| Imax | 0.56 | fraction | Near-maximal sTREM2 suppression |
| IC50 | 0.70 | ng/mL CSF | Sponsor exposure-response curve anchor |
| Kout | 0.016 | h^-1 | Biomarker turnover half-life ~43 h |

### Visible calibration anchors from the slide deck
- **MAD 25 mg QD Day 15:** Cmax ~110 ng/mL, trough ~75 ng/mL, sTREM2 ~-43%
- **MAD 50 mg QD Day 15:** sTREM2 plateau ~-44%
- **SAD 72 mg:** HV and AD plasma PK broadly overlap, supporting a shared structural model
- **Exposure-response:** plateau begins around **2-5 ng/mL CSF** with a projected maximal reduction around **50-55%**

### Important modeling note
The source is a conference slide deck, not a full population PK manuscript. To keep the simulator useful and honest, the PK layer is calibrated against the visible concentration-time curves and the published sponsor exposure-response slide rather than a formal author-reported CL/F and V/F table.

For the clinical time-course output, the displayed MAD biomarker response is capped near the observed sponsor plateau (~44% reduction) because the slide deck shows minimal incremental clinical suppression from 25 to 50 mg even though the separate exposure-response slide projects a somewhat deeper asymptote at higher CSF concentrations.
")
      ),
      nav_panel(
        "References",
        div(
          class = "ref-box",
          tags$h5("Primary reference"),
          tags$ol(
            tags$li("Birch E et al. Phase 1 single and multiple ascending dose study of VG-3927 in healthy volunteers and Alzheimer's disease participants. AAIC 2025 clinical presentation.")
          ),
          tags$h5("Therapeutic context"),
          tags$ul(
            tags$li(tags$strong("Class:"), " Oral TREM2 agonist"),
            tags$li(tags$strong("Therapeutic area:"), " Alzheimer's disease"),
            tags$li(tags$strong("Phase 2 dose selection:"), " 25 mg once daily"),
            tags$li(tags$strong("PK message:"), " Predictable dose-dependent exposure with low variability and flat QD peak-to-trough behavior"),
            tags$li(tags$strong("PD message:"), " CSF sTREM2 suppression plateaued near 25 mg QD with projected maximal reduction around 50-55%")
          )
        )
      )
    )
  ),
  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({
    shiny::req(input$dose, input$interval, input$n_days, input$wt)

    simulate_vg3927(
      dose = input$dose,
      interval = input$interval,
      n_days = input$n_days,
      wt = input$wt,
      age_group = input$age_group,
      csf_multiplier = input$csf_multiplier,
      pd_multiplier = input$pd_multiplier
    )
  })

  last_interval <- reactive({
    d <- sim_data()
    start_time <- max(0, max(d$time_h) - input$interval)

    dplyr::filter(d, time_h >= start_time)
  })

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

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

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

  output$thalf <- renderText({
    cl_i <- 11.4 * (input$wt / 70)^0.75
    v_i <- 476 * (input$wt / 70)
    sprintf("%.1f", log(2) / (cl_i / v_i))
  })

  output$simulation_note <- renderText({
    d <- sim_data()
    ss <- last_interval()
    day15 <- dplyr::filter(d, abs(time_day - min(input$n_days, 15)) < 0.15)
    pd_day <- if (nrow(day15) > 0) median(day15$STREM2_PCT, na.rm = TRUE) else tail(d$STREM2_PCT, 1)

    sprintf(
      "Predicted steady-state summary at %s mg Q%dH: plasma Cmax %.0f ng/mL, trough %.0f ng/mL, CSF Cmax %.2f ng/mL, and sTREM2 change %.0f%% by Day %d.",
      input$dose,
      input$interval,
      max(ss$CP, na.rm = TRUE),
      min(ss$CP, na.rm = TRUE),
      max(ss$CCSF, na.rm = TRUE),
      pd_day,
      min(input$n_days, 15)
    )
  })

  output$main_plot <- renderPlot({
    d <- sim_data()

    if (identical(input$view, "Plasma PK")) {
      plot_df <- d
      if (isTRUE(input$log_scale)) {
        plot_df <- dplyr::filter(plot_df, CP > 0)
      }

      p <- ggplot(plot_df, aes(x = time_day, y = CP)) +
        geom_line(color = "#8b5cf6", linewidth = 0.9) +
        geom_point(
          data = dplyr::filter(
            slide_pk_anchors,
            dose == input$dose,
            age_group == input$age_group
          ),
          aes(x = time_h / 24, y = cp),
          inherit.aes = FALSE,
          color = "#f59e0b",
          size = 2.8
        ) +
        labs(
          x = "Time (days)",
          y = "Plasma concentration (ng/mL)",
          title = paste0("VG-3927 plasma PK: ", input$dose, " mg Q", input$interval, "H"),
          subtitle = "Orange points mark slide-derived Day 15 anchors when available"
        ) +
        theme_minimal(base_size = 14)

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

      p
    } else if (identical(input$view, "CSF Exposure")) {
      plot_df <- d
      if (isTRUE(input$log_scale)) {
        plot_df <- dplyr::filter(plot_df, CCSF > 0)
      }

      p <- ggplot(plot_df, aes(x = time_day, y = CCSF)) +
        geom_rect(
          inherit.aes = FALSE,
          xmin = -Inf,
          xmax = Inf,
          ymin = 2,
          ymax = 5,
          fill = "#10b981",
          alpha = 0.12
        ) +
        geom_line(color = "#0dcaf0", linewidth = 0.9) +
        geom_hline(yintercept = 0.70, linetype = "dashed", color = "#8b5cf6") +
        annotate(
          "text",
          x = max(plot_df$time_day) * 0.04,
          y = 0.78,
          label = "IC50 ~0.7 ng/mL CSF",
          hjust = 0,
          color = "#8b5cf6",
          size = 3.5
        ) +
        labs(
          x = "Time (days)",
          y = "CSF concentration (ng/mL)",
          title = "Projected CSF exposure",
          subtitle = "Green band highlights the sponsor slide plateau region (~2-5 ng/mL)"
        ) +
        theme_minimal(base_size = 14)

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

      p
    } else if (identical(input$view, "sTREM2 PD")) {
      ggplot(d, aes(x = time_day, y = STREM2_PCT)) +
        geom_hline(yintercept = c(-25, -45), linetype = "dashed", color = c("#f59e0b", "#10b981")) +
        geom_line(color = "#ef4444", linewidth = 0.9) +
        labs(
          x = "Time (days)",
          y = "sTREM2 % change from baseline",
          title = "Projected CSF sTREM2 response",
          subtitle = "Dashed lines show mid-range and plateau biomarker suppression"
        ) +
        theme_minimal(base_size = 14)
    } else {
      exposure_curve <- data.frame(
        ccsf = seq(0, 6, length.out = 200)
      ) |>
        mutate(effect = -100 * 0.56 * ccsf / (0.70 / input$pd_multiplier + ccsf))

      ggplot(exposure_curve, aes(x = ccsf, y = effect)) +
        geom_line(color = "#8b5cf6", linewidth = 0.9) +
        geom_point(
          data = slide_pd_anchors,
          aes(x = ccsf, y = effect),
          inherit.aes = FALSE,
          color = "#ef4444",
          size = 2.8
        ) +
        geom_text(
          data = slide_pd_anchors,
          aes(x = ccsf, y = effect, label = label),
          inherit.aes = FALSE,
          nudge_y = 2,
          size = 3.5,
          color = "#4b5563"
        ) +
        labs(
          x = "CSF VG-3927 concentration (ng/mL)",
          y = "sTREM2 % change from baseline",
          title = "CSF exposure-response",
          subtitle = "Red points approximate the sponsor Day 15/20 scatter"
        ) +
        theme_minimal(base_size = 14)
    }
  })
}

shinyApp(ui = ui, server = server)