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

# Salinger et al. 2013 BJOG — Magnesium Sulphate Population PK
# Magnesium sulphate for prevention of eclampsia: IV vs IM regimens
# Final PopPK Model: 1-CMT with IM absorption, weight & creatinine covariates

# Conversion utility: g MgSO4·7H2O → mmol Mg2+ (MW = 246.47 g/mol)
g_to_mmol <- function(g) (g * 1000) / 246.47

model_code <- r"(
$PARAM @annotated
CL   : 4.81  : Clearance (L/h) - typical at Cr = 0.8 mg/dL, WT = 55 kg
V    : 15.6  : Volume of distribution (L) - typical at WT = 55 kg
KA   : 0.317 : IM absorption rate constant (1/h)
FIM  : 0.862 : IM bioavailability (fraction)
BL   : 0.85  : Baseline endogenous Mg concentration (mmol/L)
WT   : 55    : Body weight (kg)
CR   : 0.8   : Serum creatinine (mg/dL)
h1   : 0.692 : Weight exponent on V (V ~ WT^0.692)
h2   : 1.480 : Creatinine exponent on CL (CL ~ (0.8/Cr)^1.48)

$CMT @annotated
DEPOT : IM absorption depot (mmol)
CENT  : Central compartment - Mg above baseline (mmol)

$MAIN
double CLi = CL * pow(0.8 / CR, h2);
double Vi  = V  * pow(WT / 55.0, h1);
F_DEPOT    = FIM;

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

$TABLE
double CP_mmol = CENT / Vi + BL;
double CP_mgdl = CP_mmol * 2.431;

$CAPTURE CP_mmol CP_mgdl CLi Vi
)"

mod <- mcode("magnesium_salinger2013", 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; }
  .metric-danger .metric-value { color: #ef4444; }
  .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; }
  .param-table td { padding: 4px 10px; }
  .param-table tr:nth-child(even) { background: #f8f9fa; }
")

ui <- page_sidebar(
  title = "Magnesium Sulphate (MgSO\u2084) PK Simulator — Eclampsia Prevention",
  theme = app_theme,
  sidebar = sidebar(
    title = "Dosing & Patient Settings", width = 340,
    h6("IV Loading Dose (Adjustable)"),
    sliderInput("iv_load_g",
      "IV Loading Dose (g MgSO\u2084\u00b77H\u2082O)",
      min = 4, max = 8, value = 4, step = 1
    ),
    tags$small(
      class = "text-muted",
      "Standard = 4g. Increase to 6-8g to match IM initial concentrations (Fig 4)."
    ),
    hr(),
    h6("Simulation Settings"),
    selectInput("dur_h", "Duration",
      choices = c("12 hours" = 12, "24 hours" = 24),
      selected = 24
    ),
    hr(),
    h6("Patient Characteristics"),
    sliderInput("wt", "Body Weight (kg)", min = 40, max = 90, value = 55, step = 1),
    tags$small(class = "text-muted", "Median study weight = 55 kg (5th-95th%ile: 42-78 kg)"),
    br(), br(),
    sliderInput("cr", "Serum Creatinine (mg/dL)", min = 0.5, max = 1.5, value = 0.8, step = 0.05),
    tags$small(class = "text-muted", "Median = 0.8 mg/dL (5th-95th%ile: 0.6-1.2 mg/dL)"),
    hr(),
    h6("Display"),
    radioButtons("units", "Concentration Units",
      choices = c("mg/dL (clinical)" = "mgdl", "mmol/L" = "mmol"),
      selected = "mgdl", inline = FALSE
    ),
    checkboxInput("log_scale", "Log Scale (Y-axis)", value = FALSE)
  ),
  navset_card_underline(
    title = "Magnesium Sulphate PK Simulator",
    full_screen = TRUE,
    nav_panel(
      "Simulation",
      layout_column_wrap(
        width = 1 / 4, fill = FALSE,
        div(
          class = "metric-card metric-primary",
          div(class = "metric-value", textOutput("cmax_iv")),
          div(class = "metric-label", "IV Cmax (mg/dL)")
        ),
        div(
          class = "metric-card metric-info",
          div(class = "metric-value", textOutput("cmax_im")),
          div(class = "metric-label", "IM+IV Cmax (mg/dL)")
        ),
        div(
          class = "metric-card metric-warning",
          div(class = "metric-value", textOutput("c12h")),
          div(class = "metric-label", "IV Conc at 12h (mg/dL)")
        ),
        div(
          class = "metric-card metric-success",
          div(class = "metric-value", textOutput("thalf")),
          div(class = "metric-label", "t\u00bd (h)")
        )
      ),
      card(
        full_screen = TRUE,
        height = "500px",
        plotOutput("pkPlot", height = "460px")
      )
    ),
    nav_panel(
      "Model Information",
      markdown("
## Magnesium Sulphate (MgSO\u2084) — Population PK Model

**Source:** Salinger DH et al. *BJOG* 2013;120:894-900. DOI: 10.1111/1471-0528.12222

**Study:** Randomised trial (n = 258 pregnant women with pre-eclampsia, India) comparing
IM and IV MgSO\u2084 regimens for eclampsia prevention.

---

### Model Structure
- **Type:** 1-compartment with first-order IM absorption
- **Software:** NONMEM 7
- **Residual error:** Proportional (25% CV Base; 22.9% CV Final)

### Final Model Parameters (Table 2)

| Parameter | Value | SE | Units |
|---|---|---|---|
| CL | 48.1 (4.81) | 10% | dL/h (L/h) |
| V | 156 (15.6) | 8% | dL (L) |
| KA (IM) | 0.317 | 32% | h\u207b\u00b9 |
| F (IM) | 86.2% | 3.2% | — |
| Baseline Mg | 0.85 (2.08) | 3.2% | mmol/L (mg/dL) |
| Weight exponent (h\u2081) | 0.692 | 18% | — |
| Creatinine exponent (h\u2082) | 1.48 | 11% | — |
| Residual CV | 22.9% | — | — |

### Covariate Effects
- **V:** V\u1d62 = V \u00d7 (WT\u1d62 / 55)\u2070\u00b7\u2076\u2079\u00b2
- **CL:** CL\u1d62 = CL \u00d7 (0.8 / Cr\u1d62)\u00b9\u00b7\u2074\u2078

### Dosing Regimens (MAGPIE Trial protocol)
**IV Regimen:**  
\u2022 Loading: 4g IV over 20 min  
\u2022 Maintenance: 1 g/hr continuous IV infusion

**IM Regimen:**  
\u2022 Loading: 4g IV over 20 min + 10g IM (5g each buttock)  
\u2022 Maintenance: 5g IM q4h

### Therapeutic Drug Monitoring
| Endpoint | Mg Concentration | Equivalence |
|---|---|---|
| Therapeutic minimum | \u22654.8 mg/dL | \u22652 mmol/L; \u22654 mEq/L |
| Therapeutic target | 4.8\u20138.4 mg/dL | 2\u20133.5 mmol/L |
| Patellar reflex loss | 8.5\u201312.2 mg/dL | 3.5\u20135 mmol/L |
| Respiratory paralysis | 12.2\u201315.8 mg/dL | 5\u20136.5 mmol/L |
| Cardiac arrest risk | >30.4 mg/dL | >12.5 mmol/L |

### Key Findings
- IM regimen produces **higher early concentrations** (larger loading dose)
- At steady state (12h+), both regimens produce **comparable concentrations**  
- **~25% of IV women** had concentrations <3.5 mg/dL (subtherapeutic)
- Increasing IV loading dose to **6g** closely matches IM early PK profile
- Bioavailability of IM MgSO\u2084 = **86.2%** (higher than assumed in protocol design)

### Unit Conversions
- 1 mmol/L Mg = 2.431 mg/dL = 2 mEq/L
- Therapeutic minimum 4 mEq/L = 2 mmol/L = **4.86 mg/dL**
      ")
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("\U0001F4DA Key References"),
        tags$ol(
          tags$li(
            tags$strong("Salinger DH et al. (2013)"),
            " Magnesium sulphate for prevention of eclampsia: are IM and IV regimens equivalent? A population pharmacokinetic study. ",
            tags$em("BJOG"), " 120:894\u2013900. ",
            tags$a(href = "https://doi.org/10.1111/1471-0528.12222", target = "_blank", "DOI: 10.1111/1471-0528.12222")
          ),
          tags$li(
            "Chuan FS et al. (2001) Population pharmacokinetics of magnesium in preeclampsia. ",
            tags$em("Am J Obstet Gynecol"), " 185:593\u20139."
          ),
          tags$li(
            "Lu J et al. (2002) PK-PD modelling of magnesium plasma concentration and blood pressure in pre-eclamptic women. ",
            tags$em("Clin Pharmacokinet"), " 41:1105\u201313."
          ),
          tags$li(
            "Magpie Trial Collaboration Group (2002) Do women with pre-eclampsia benefit from MgSO\u2084? ",
            tags$em("Lancet"), " 359:1877\u201390."
          )
        ),
        tags$h5("\U0001F48A Therapeutic Context"),
        tags$ul(
          tags$li(tags$strong("Drug class:"), " Electrolyte / anticonvulsant"),
          tags$li(tags$strong("Indication:"), " Eclampsia prevention and treatment in pre-eclampsia"),
          tags$li(tags$strong("Population:"), " Pregnant women with hypertensive disease (pre-eclampsia)"),
          tags$li(tags$strong("Route:"), " IV infusion or deep IM injection"),
          tags$li(tags$strong("Monitoring:"), " Deep tendon reflexes, respiratory rate, urine output, serum Mg")
        )
      )
    )
  ),
  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 Model: Salinger et al. BJOG 2013 \u2022 Built with mrgsolve", 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$iv_load_g, input$wt, input$cr, input$dur_h)

    dur   <- as.numeric(input$dur_h)
    use_mgdl <- input$units == "mgdl"

    # Build covariate-adjusted model
    mod_p <- mod %>% param(WT = input$wt, CR = input$cr)

    # --- IV Regimen ---
    iv_load_mmol  <- g_to_mmol(input$iv_load_g)
    iv_maint_rate <- g_to_mmol(1) # 1 g/hr = X mmol/hr

    iv_load_ev  <- ev(amt = iv_load_mmol,  cmt = 2,
                      rate = iv_load_mmol / (20 / 60), time = 0)
    iv_maint_ev <- ev(amt = iv_maint_rate * dur, cmt = 2,
                      rate = iv_maint_rate, time = 0)
    iv_events   <- c(iv_load_ev, iv_maint_ev)

    sim_iv <- mod_p %>%
      ev(iv_events) %>%
      mrgsim(end = dur, delta = 0.1) %>%
      as.data.frame() %>%
      mutate(
        regimen   = paste0("IV (", input$iv_load_g, "g load + 1g/hr)"),
        CP_display = if (use_mgdl) CP_mgdl else CP_mmol
      )

    # --- IM Regimen (fixed 4g IV + 10g IM loading, 5g IM q4h) ---
    im_iv_load_mmol <- g_to_mmol(4) # IM protocol always uses 4g IV load
    n_im_maint      <- max(1, floor((dur - 4) / 4))

    im_iv_load_ev <- ev(amt = im_iv_load_mmol, cmt = 2,
                        rate = im_iv_load_mmol / (20 / 60), time = 0)
    im_load_ev    <- ev(amt = g_to_mmol(10), cmt = 1, time = 0)
    im_maint_ev   <- ev(amt = g_to_mmol(5), cmt = 1,
                        ii = 4, addl = n_im_maint - 1, time = 4)
    im_events     <- c(im_iv_load_ev, im_load_ev, im_maint_ev)

    sim_im <- mod_p %>%
      ev(im_events) %>%
      mrgsim(end = dur, delta = 0.1) %>%
      as.data.frame() %>%
      mutate(
        regimen    = "IM (4g IV + 10g IM load; 5g IM q4h)",
        CP_display = if (use_mgdl) CP_mmol * 2.431 else CP_mmol
      )

    bind_rows(sim_iv, sim_im)
  })

  # Helper: extract data for one regimen
  regimen_data <- function(d, label_prefix) {
    dplyr::filter(d, grepl(label_prefix, regimen, fixed = FALSE))
  }

  output$cmax_iv <- renderText({
    d <- sim_data()
    iv_d <- regimen_data(d, "^IV")
    if (input$units == "mgdl") {
      sprintf("%.1f", max(iv_d$CP_display, na.rm = TRUE))
    } else {
      sprintf("%.2f", max(iv_d$CP_display, na.rm = TRUE))
    }
  })

  output$cmax_im <- renderText({
    d <- sim_data()
    im_d <- regimen_data(d, "^IM")
    if (input$units == "mgdl") {
      sprintf("%.1f", max(im_d$CP_display, na.rm = TRUE))
    } else {
      sprintf("%.2f", max(im_d$CP_display, na.rm = TRUE))
    }
  })

  output$c12h <- renderText({
    d <- sim_data()
    iv_d <- regimen_data(d, "^IV")
    val <- iv_d |> dplyr::filter(abs(time - 12) < 0.15) |>
      dplyr::pull(CP_display)
    if (length(val) == 0) val <- tail(iv_d$CP_display, 1)
    if (input$units == "mgdl") sprintf("%.1f", mean(val, na.rm = TRUE)) else sprintf("%.2f", mean(val, na.rm = TRUE))
  })

  output$thalf <- renderText({
    shiny::req(input$wt, input$cr)
    CLi <- 4.81 * (0.8 / input$cr)^1.48
    Vi  <- 15.6 * (input$wt / 55)^0.692
    sprintf("%.1f", log(2) / (CLi / Vi))
  })

  output$pkPlot <- renderPlot({
    d     <- sim_data()
    dur   <- as.numeric(input$dur_h)
    units <- input$units

    y_lab      <- if (units == "mgdl") "Serum Mg Concentration (mg/dL)" else "Serum Mg Concentration (mmol/L)"
    thr_low    <- if (units == "mgdl") 4.86  else 2.0    # therapeutic minimum
    thr_high   <- if (units == "mgdl") 8.44  else 3.47   # upper therapeutic
    tox_warn   <- if (units == "mgdl") 8.44  else 3.47   # patellar reflex risk (same as thr_high)
    tox_resp   <- if (units == "mgdl") 12.15 else 5.0    # respiratory paralysis
    bl_line    <- if (units == "mgdl") 2.07  else 0.85   # baseline

    colors <- c(
      "IV"  = "#8b5cf6",
      "IM"  = "#0ea5e9"
    )

    # Simplify regimen labels for legend
    d_plot <- d |> mutate(
      regimen_short = dplyr::case_when(
        grepl("^IV", regimen) ~ "IV",
        grepl("^IM", regimen) ~ "IM",
        TRUE                  ~ regimen
      )
    )

    p <- ggplot(d_plot, aes(x = time, y = CP_display, color = regimen_short)) +
      # Therapeutic window (green band)
      annotate("rect",
        xmin = -Inf, xmax = Inf,
        ymin = thr_low, ymax = thr_high,
        fill = "#10b981", alpha = 0.10
      ) +
      geom_hline(yintercept = thr_low,  linetype = "dashed",
                 color = "#10b981", linewidth = 0.5, alpha = 0.8) +
      geom_hline(yintercept = thr_high, linetype = "dashed",
                 color = "#f59e0b", linewidth = 0.5, alpha = 0.8) +
      geom_hline(yintercept = tox_resp, linetype = "dotted",
                 color = "#ef4444", linewidth = 0.5, alpha = 0.7) +
      geom_hline(yintercept = bl_line,  linetype = "solid",
                 color = "#9ca3af", linewidth = 0.4, alpha = 0.6) +
      # Threshold labels
      annotate("text", x = dur * 0.02, y = thr_low + 0.15,
               label = "Therapeutic min (~4.9 mg/dL)", hjust = 0,
               size = 3.2, color = "#10b981") +
      annotate("text", x = dur * 0.02, y = thr_high + 0.15,
               label = "Patellar reflex risk (8.4 mg/dL)", hjust = 0,
               size = 3.2, color = "#f59e0b") +
      annotate("text", x = dur * 0.02, y = tox_resp + 0.15,
               label = "Respiratory paralysis risk (12.2 mg/dL)", hjust = 0,
               size = 3.2, color = "#ef4444") +
      # PK curves
      geom_line(linewidth = 1.1) +
      scale_color_manual(
        values = c("IV" = "#8b5cf6", "IM" = "#0ea5e9"),
        labels = c(
          "IV" = paste0("IV (", input$iv_load_g, "g load + 1g/hr maint.)"),
          "IM" = "IM (4g IV + 10g IM load; 5g IM q4h)"
        ),
        name = "Regimen"
      ) +
      labs(
        x     = "Time (hours)",
        y     = y_lab,
        title = paste0(
          "Magnesium Sulphate PK — WT = ", input$wt, " kg | ",
          "Cr = ", input$cr, " mg/dL"
        ),
        subtitle = paste0(
          "CL = ", round(4.81 * (0.8 / input$cr)^1.48, 2), " L/h | ",
          "V = ", round(15.6 * (input$wt / 55)^0.692, 1), " L | ",
          "t\u00bd = ", round(log(2) / (4.81 * (0.8 / input$cr)^1.48 /
                                       (15.6 * (input$wt / 55)^0.692)), 1), " h"
        )
      ) +
      theme_minimal(base_size = 13) +
      theme(
        legend.position = "bottom",
        plot.title    = element_text(face = "bold", size = 13),
        plot.subtitle = element_text(color = "#6c757d", size = 11)
      )

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

shinyApp(ui = ui, server = server)