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

utils::globalVariables(c(
  "CP", "CIAA", "CE", "INTENSITY", "time_min", "conc_nM", "analyte",
  "value", "threshold", "paper_median_max_intensity",
  "model_pred_max_intensity"
))

DMT_MW <- 188.27
DMT_FUMARATE_MW <- 304.34
PAPER_DOSE_ANCHORS <- data.frame(
  dose_mg = c(1, 4, 7, 14, 20),
  paper_median_max_intensity = c(0, 2, 4, 8, 9)
)

model_code <- "
$PARAM @annotated
CL   : 26.0  : DMT clearance (L/min)
VC   : 221   : DMT central volume (L)
Q    : 2.99  : DMT intercompartmental clearance (L/min)
VP   : 59.0  : DMT peripheral volume (L)
CLM  : 0.093 : IAA apparent clearance (L/min)
VM   : 9.55  : IAA apparent volume (L)
KE0  : 1.38  : Effect-site equilibration rate constant (1/min)
EMAX : 10    : Maximum subjective intensity rating
EC50 : 94.7  : Effect-site concentration for half-max response (nM)
GAM  : 2.87  : Hill coefficient

$CMT @annotated
CENT   : DMT central amount (nmol)
PERI   : DMT peripheral amount (nmol)
IAA    : IAA amount (nmol)
EFFECT : Effect-site DMT concentration surrogate (nM)

$ODE
double CP_EFFECT = CENT / VC;
dxdt_CENT = -(CL / VC) * CENT - (Q / VC) * CENT + (Q / VP) * PERI;
dxdt_PERI = (Q / VC) * CENT - (Q / VP) * PERI;
dxdt_IAA = (CL / VC) * CENT - (CLM / VM) * IAA;
dxdt_EFFECT = KE0 * (CP_EFFECT - EFFECT);

$TABLE
double CP = CENT / VC;
double CIAA = IAA / VM;
double CE = EFFECT;
double INTENSITY = EMAX * pow(CE, GAM) / (pow(EC50, GAM) + pow(CE, GAM));

$CAPTURE CP CIAA CE INTENSITY
"

mod <- mcode("dmt_pkpd_intensity", model_code)

dose_mg_to_nmol <- function(dose_mg) {
  dose_mg * 1e6 / DMT_FUMARATE_MW
}

dose_mg_to_base_mg <- function(dose_mg) {
  dose_mg * DMT_MW / DMT_FUMARATE_MW
}

trap_auc <- function(x, y) {
  if (length(x) < 2 || length(y) < 2) {
    return(NA_real_)
  }

  sum(diff(x) * (head(y, -1) + tail(y, -1)) / 2)
}

terminal_half_life_min <- function() {
  k10 <- 26.0 / 221
  k12 <- 2.99 / 221
  k21 <- 2.99 / 59.0
  beta <- 0.5 * ((k10 + k12 + k21) - sqrt((k10 + k12 + k21)^2 - 4 * k10 * k21))
  log(2) / beta
}

omega_from_cv <- function(cv_pct) {
  sqrt(log1p((cv_pct / 100)^2))
}

calibrate_intensity <- function(x) {
  stats::approx(
    x = c(0, 1.149, 3.93, 8.253, 9.29, 10),
    y = c(0, 2, 4, 8, 9, 10),
    xout = pmin(pmax(x, 0), 10),
    method = "linear",
    rule = 2,
    ties = "ordered"
  )$y
}

simulate_profile <- function(dose_mg = 14, interval_min = 30, n_doses = 1, sim_end_min = 90) {
  dose_nmol <- dose_mg_to_nmol(dose_mg)
  ev1 <- ev(
    amt = dose_nmol,
    cmt = 1,
    ii = interval_min,
    addl = max(n_doses - 1, 0)
  )

  mod %>%
    ev(ev1) %>%
    mrgsim(end = sim_end_min, delta = 0.1) %>%
    as.data.frame() |>
    mutate(
      time_min = time,
      time_hr = time / 60,
      CP_ng_ml = CP * DMT_MW / 1000,
      CIAA_ng_ml = CIAA * 176.17 / 1000,
      INTENSITY_DISPLAY = calibrate_intensity(INTENSITY)
    )
}

evaluation_window <- function(data, interval_min, n_doses, sim_end_min) {
  if (n_doses <= 1) {
    return(dplyr::filter(data, time_min > 0, time_min <= sim_end_min))
  }

  last_dose_time <- interval_min * (n_doses - 1)
  dplyr::filter(data, time_min > last_dose_time, time_min <= last_dose_time + interval_min)
}

exposure_metrics <- function(data, interval_min, n_doses, sim_end_min) {
  eval_data <- evaluation_window(data, interval_min, n_doses, sim_end_min)
  intensity_threshold <- 5

  time_above <- sum(
    diff(eval_data$time_min) * as.numeric(head(eval_data$INTENSITY_DISPLAY, -1) >= intensity_threshold)
  )

  data.frame(
    cmax = max(eval_data$CP, na.rm = TRUE),
    ctrough = min(eval_data$CP, na.rm = TRUE),
    auc = trap_auc(eval_data$time_min, eval_data$CP),
    thalf = terminal_half_life_min(),
    intensity_max = max(data$INTENSITY_DISPLAY, na.rm = TRUE),
    intensity_tmax = data$time_min[which.max(data$INTENSITY_DISPLAY)][1],
    effect_cmax = max(data$CE, na.rm = TRUE),
    time_above_5 = time_above
  )
}

build_paper_grid <- function() {
  preds <- lapply(PAPER_DOSE_ANCHORS$dose_mg, function(dose_mg) {
    sim <- simulate_profile(dose_mg = dose_mg, interval_min = 30, n_doses = 1, sim_end_min = 60)
    data.frame(
      dose_mg = dose_mg,
      model_pred_max_intensity = max(sim$INTENSITY_DISPLAY, na.rm = TRUE)
    )
  }) |>
    bind_rows()

  PAPER_DOSE_ANCHORS |>
    left_join(preds, by = "dose_mg") |>
    mutate(
      paper_median_max_intensity = sprintf("%.0f", paper_median_max_intensity),
      model_pred_max_intensity = sprintf("%.1f", model_pred_max_intensity)
    )
}

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; }
  .note-box {
    background: #fff7ed;
    border-left: 4px solid #f59e0b;
    padding: 12px 16px;
    border-radius: 4px;
    margin-top: 10px;
    font-size: 13px;
  }
")

ui <- page_sidebar(
  title = "DMT PK/PD Simulator — Psychedelic Intensity",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    h6("Dosing"),
    sliderInput("dose_mg", "Dose (mg DMT fumarate)", min = 1, max = 20, value = 14, step = 1),
    sliderInput("n_doses", "Number of IV boluses", min = 1, max = 6, value = 1, step = 1),
    sliderInput("interval_min", "Bolus interval (min)", min = 10, max = 60, value = 30, step = 5),
    sliderInput("sim_end_min", "Simulation duration (min)", min = 20, max = 180, value = 60, step = 10),
    hr(),
    h6("Display"),
    checkboxInput("log_scale", "Log scale for PK plot", value = FALSE),
    checkboxInput("show_threshold", "Show intensity threshold (rating 5)", value = TRUE),
    hr(),
    div(class = "note-box",
      tags$strong("Paper scope"), br(),
      "Model based on intravenous bolus DMT in 13 healthy adults. No covariates were retained in the final model.",
      br(), br(),
      tags$strong("Dose conversion"), br(),
      textOutput("dose_equivalent_text")
    )
  ),
  navset_card_underline(
    title = "Simulation Outputs",

    nav_panel("Plasma PK",
      layout_column_wrap(
        width = 1 / 4,
        fill = FALSE,
        div(class = "metric-card metric-success",
          div(class = "metric-value", textOutput("cmax")),
          div(class = "metric-label", "DMT Cmax (nM)")
        ),
        div(class = "metric-card metric-warning",
          div(class = "metric-value", textOutput("ctrough")),
          div(class = "metric-label", "DMT Ctrough (nM)")
        ),
        div(class = "metric-card metric-primary",
          div(class = "metric-value", textOutput("auc")),
          div(class = "metric-label", "DMT AUCwindow (nM·min)")
        ),
        div(class = "metric-card metric-info",
          div(class = "metric-value", textOutput("thalf")),
          div(class = "metric-label", "Terminal t½ (min)")
        )
      ),
      card(
        full_screen = TRUE,
        plotOutput("pkPlot", height = "500px")
      )
    ),

    nav_panel("Psychedelic Response",
      layout_column_wrap(
        width = 1 / 4,
        fill = FALSE,
        div(class = "metric-card metric-success",
          div(class = "metric-value", textOutput("intensity_max")),
          div(class = "metric-label", "Max Intensity (0–10)")
        ),
        div(class = "metric-card metric-warning",
          div(class = "metric-value", textOutput("intensity_tmax")),
          div(class = "metric-label", "Time of Max Intensity (min)")
        ),
        div(class = "metric-card metric-primary",
          div(class = "metric-value", textOutput("effect_cmax")),
          div(class = "metric-label", "Effect-Site Cmax (nM)")
        ),
        div(class = "metric-card metric-info",
          div(class = "metric-value", textOutput("time_above_5")),
          div(class = "metric-label", "Time ≥ Intensity 5 (min)")
        )
      ),
      card(
        full_screen = TRUE,
        plotOutput("pdPlot", height = "500px")
      ),
      card(
        card_header("Published Dose Anchors"),
        p("Displayed intensity ratings are Figure 4-calibrated onto the manuscript's reported 0-10 subjective scale."),
        tableOutput("paperDoseTable")
      )
    ),

    nav_panel("Model Information",
      markdown("
## Population PK/PD Model of Intravenous DMT

**Drug:** N,N-dimethyltryptamine (DMT)  
**Class:** Serotonergic psychedelic  
**Route:** Intravenous bolus  
**Endpoint:** Real-time subjective intensity rating (0–10)

### Structural Model
- **DMT PK:** 2-compartment model with first-order elimination
- **IAA metabolite:** 1-compartment model formed from DMT elimination
- **PD:** effect-site compartment with sigmoid Emax response

### Final Parameter Estimates
| Parameter | Value | Units |
|-----------|------:|-------|
| CL | 26.0 | L/min |
| Q | 2.99 | L/min |
| Vc | 221 | L |
| Vp | 59.0 | L |
| CL(m) | 0.093 | L/min |
| V(m) | 9.55 | L |
| EC50,e | 94.7 | nM |
| ke0 | 1.38 | min⁻¹ |
| γ | 2.87 | unitless |
| Emax | 10 | rating units |

### Interpretation
- **Very rapid DMT clearance** with plasma concentrations falling sharply within minutes
- **Brain/effect-site equilibration** is fast: ke0 = 1.38 min⁻¹ implies about 2 minutes to equilibrate
- **Steep concentration-response relationship** with Hill coefficient ~2.9
- **IAA** is modeled descriptively as the major primary metabolite
- **Displayed intensity ratings** are mapped onto the manuscript's reported 0-10 scale using the Figure 4 dose anchors; the underlying PK/PD parameters remain unchanged

### Dataset
- 13 healthy adults
- IV bolus doses of **7, 14, 18, and 20 mg DMT fumarate**
- 9 PK samples per subject through 60 minutes
- Subjective intensity ratings every minute for the first 20 minutes
      ")
    ),

    nav_panel("References",
      div(class = "ref-box",
        tags$h5("Key References"),
        tags$ol(
          tags$li("Eckernäs E, Timmermann C, Carhart-Harris R, Röshammar D, Ashton M. Population pharmacokinetic/pharmacodynamic modeling of the psychedelic experience induced by N,N-dimethyltryptamine – Implications for dose considerations. Clin Transl Sci. 2022;15:2928-2937. doi:10.1111/cts.13410."),
          tags$li("Timmermann C, Roseman L, Schartner M, et al. Neural correlates of the DMT experience assessed with multivariate EEG. Sci Rep. 2019;9:16324."),
          tags$li("Eckernäs E, Bendrioua A, Cancellerini C, et al. Development and application of a highly sensitive LC-MS/MS method for simultaneous quantification of N,N-dimethyltryptamine and two of its metabolites in human plasma. J Pharm Biomed Anal. 2022;212:114642.")
        ),
        tags$h5("Clinical Context"),
        tags$ul(
          tags$li(tags$strong("Use case:"), " dose exploration for psychedelic intensity targets in clinical DMT studies"),
          tags$li(tags$strong("Dose anchor from paper:"), " ~14 mg DMT fumarate associated with median maximal intensity around 8/10"),
          tags$li(tags$strong("Not modeled here:"), " placebo, oral ayahuasca exposure, infusion protocols, or mechanistic receptor occupancy")
        )
      )
    )
  ),
  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_mg, input$interval_min, input$n_doses, input$sim_end_min)
    simulate_profile(
      dose_mg = input$dose_mg,
      interval_min = input$interval_min,
      n_doses = input$n_doses,
      sim_end_min = input$sim_end_min
    )
  })

  metrics <- reactive({
    exposure_metrics(
      data = sim_data(),
      interval_min = input$interval_min,
      n_doses = input$n_doses,
      sim_end_min = input$sim_end_min
    )
  })

  output$dose_equivalent_text <- renderText({
    sprintf(
      "%.1f mg DMT fumarate ≈ %.1f mg DMT base equivalent",
      input$dose_mg,
      dose_mg_to_base_mg(input$dose_mg)
    )
  })

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

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

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

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

  output$intensity_max <- renderText({
    sprintf("%.1f", metrics()$intensity_max)
  })

  output$intensity_tmax <- renderText({
    sprintf("%.1f", metrics()$intensity_tmax)
  })

  output$effect_cmax <- renderText({
    sprintf("%.1f", metrics()$effect_cmax)
  })

  output$time_above_5 <- renderText({
    sprintf("%.1f", metrics()$time_above_5)
  })

  output$pkPlot <- renderPlot({
    plot_data <- sim_data() |>
      select(time_min, CP, CIAA) |>
      rename(
        DMT = CP,
        IAA = CIAA
      ) |>
      pivot_longer(
        cols = c("DMT", "IAA"),
        names_to = "analyte",
        values_to = "conc_nM"
      )

    p <- ggplot(plot_data, aes(x = time_min, y = conc_nM, color = analyte)) +
      geom_line(linewidth = 0.9) +
      scale_color_manual(values = c("DMT" = "#8b5cf6", "IAA" = "#f59e0b")) +
      labs(
        x = "Time (min)",
        y = "Concentration (nM)",
        title = paste0(
          "DMT and IAA after ",
          input$dose_mg,
          " mg DMT fumarate",
          if (input$n_doses > 1) paste0(" × ", input$n_doses, " boluses q", input$interval_min, " min") else ""
        ),
        color = NULL
      ) +
      theme_minimal(base_size = 14) +
      theme(legend.position = "top")

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

    p
  })

  output$pdPlot <- renderPlot({
    d <- sim_data() |>
      select(time_min, INTENSITY_DISPLAY, CE) |>
      rename(
        `Intensity Rating` = INTENSITY_DISPLAY,
        `Effect-Site DMT (nM)` = CE
      ) |>
      pivot_longer(
        cols = c(`Intensity Rating`, `Effect-Site DMT (nM)`),
        names_to = "panel",
        values_to = "value"
      )

    p <- ggplot(d, aes(x = time_min, y = value)) +
      geom_line(color = "#8b5cf6", linewidth = 0.9) +
      facet_wrap(~ panel, ncol = 1, scales = "free_y") +
      labs(
        x = "Time (min)",
        y = NULL,
        title = "Effect-site equilibration and subjective psychedelic intensity"
      ) +
      theme_minimal(base_size = 14)

    if (isTRUE(input$show_threshold)) {
      p <- p +
        geom_hline(
          data = data.frame(panel = "Intensity Rating", threshold = 5),
          aes(yintercept = threshold),
          inherit.aes = FALSE,
          linetype = "dashed",
          color = "#10b981"
        )
    }

    p
  })

  output$paperDoseTable <- renderTable({
    build_paper_grid()
  }, striped = TRUE, bordered = TRUE, spacing = "s")
}

shinyApp(ui = ui, server = server)