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

paper_d1_h <- 0.263
paper_lag_h <- 0.144

paper_reference <- function(weight_kg) {
  if (weight_kg < 30) {
    list(
      band = "<30 kg",
      recommended_dose = 2,
      cohort_median_weight = 21.5,
      cohort_size = 35,
      cmax = 57.1,
      auc = 298,
      half_life = 11.7
    )
  } else {
    list(
      band = ">=30 kg",
      recommended_dose = 4,
      cohort_median_weight = 56.9,
      cohort_size = 106,
      cmax = 50.3,
      auc = 383,
      half_life = 16.1
    )
  }
}

adult_reference <- function() {
  list(
    band = "Adult AD 4 mg reference",
    cohort_median_weight = 74.5,
    cohort_size = 819,
    recommended_dose = 4,
    cmax = 45.0,
    auc = 380,
    half_life = 12.3
  )
}

is_reference_scenario <- function(weight_kg, dose_mg) {
  ref <- paper_reference(weight_kg)
  abs(weight_kg - ref$cohort_median_weight) <= 0.1 &&
    abs(dose_mg - ref$recommended_dose) <= 0.001
}

reference_egfr <- function(weight_kg) {
  stats::approx(
    x = c(21.5, 56.9, 74.5),
    y = c(166, 117, 93),
    xout = weight_kg,
    rule = 2
  )$y
}

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

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

terminal_half_life <- function(cl, v1, q, v2) {
  k10 <- cl / v1
  k12 <- q / v1
  k21 <- q / v2
  disc <- (k10 + k12 + k21)^2 - 4 * k21 * k10
  beta <- 0.5 * ((k10 + k12 + k21) - sqrt(disc))
  log(2) / beta
}

build_regimen_events <- function(dose_mg, interval_h, n_days, lag_h) {
  n_doses <- max(1, ceiling((n_days * 24) / interval_h))
  dose_times <- seq(0, by = interval_h, length.out = n_doses) + lag_h

  tibble::tibble(
    ID = 1,
    time = dose_times,
    amt = dose_mg,
    cmt = 1,
    evid = 1,
    rate = -2
  )
}

model_code <- "
$PARAM @annotated
CLNR  : 2.76   : Apparent non-renal clearance at WT 74 kg (L/h)
CLR0  : 7.90   : Apparent renal clearance at WT 74 kg and eGFR 93 (L/h)
V1    : 119    : Apparent central volume at WT 74 kg (L)
Q     : 2.40   : Apparent intercompartmental clearance at WT 74 kg (L/h)
V2    : 46.8   : Apparent peripheral volume at WT 74 kg (L)
D1    : 0.263  : Zero-order absorption duration (h)
WT    : 74     : Body weight (kg)
EGFR  : 93     : Baseline eGFR (mL/min/1.73 m2)
DEGFR : 0      : Change from baseline eGFR

$CMT @annotated
CENT   : Central compartment (mg)
PERIPH : Peripheral compartment (mg)

$MAIN
double wt_cl = pow(WT / 74.0, 0.75);
double wt_v = pow(WT / 74.0, 1.0);
double CLr = CLR0 * ((EGFR / 93.0) + 0.00778 * DEGFR) * wt_cl;
double CLnr = CLNR * wt_cl;
double CLi = CLr + CLnr;
double Qi = Q * wt_cl;
double V1i = V1 * wt_v;
double V2i = V2 * wt_v;

D_CENT = D1;
F_CENT = 1.0;

$ODE
dxdt_CENT = -(CLi / V1i) * CENT - (Qi / V1i) * CENT + (Qi / V2i) * PERIPH;
dxdt_PERIPH = (Qi / V1i) * CENT - (Qi / V2i) * PERIPH;

$TABLE
double CP = (CENT / V1i) * 1000.0;

$CAPTURE @annotated
CP   : Plasma concentration (ng/mL)
CLr  : Apparent renal clearance (L/h)
CLnr : Apparent non-renal clearance (L/h)
CLi  : Total apparent clearance (L/h)
Qi   : Apparent intercompartmental clearance (L/h)
V1i  : Apparent central volume (L)
V2i  : Apparent peripheral volume (L)
"

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

run_simulation <- function(dose_mg, interval_h, n_days, wt, egfr, degfr, lag_h, d1_h) {
  events <- build_regimen_events(dose_mg, interval_h, n_days, lag_h)
  sim_end <- max(events$time) + interval_h

  out <- mod %>%
    param(WT = wt, EGFR = egfr, DEGFR = degfr, D1 = d1_h) %>%
    data_set(events) %>%
    mrgsim(end = sim_end, delta = 0.05) %>%
    as.data.frame() |>
    mutate(
      time_day = time / 24,
      time_hour = time
    )

  list(
    sim = out,
    events = events,
    interval_end = max(events$time) + interval_h,
    interval_start = max(events$time),
    reference = paper_reference(wt)
  )
}

compute_metrics <- function(sim_bundle) {
  sim <- sim_bundle$sim |>
    dplyr::filter(time >= sim_bundle$interval_start, time <= sim_bundle$interval_end)

  last_row <- sim[nrow(sim), , drop = FALSE]
  cmax <- max(sim$CP, na.rm = TRUE)
  ctrough <- last_row$CP[[1]]
  auc <- trap_auc(sim$time, sim$CP)
  cl <- last_row$CLi[[1]]
  q <- last_row$Qi[[1]]
  v1 <- last_row$V1i[[1]]
  v2 <- last_row$V2i[[1]]

  list(
    cmax = cmax,
    ctrough = ctrough,
    auc = auc,
    thalf = terminal_half_life(cl, v1, q, v2),
    cl = cl,
    v1 = v1,
    v2 = v2,
    clr = last_row$CLr[[1]],
    clnr = last_row$CLnr[[1]]
  )
}

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

ui <- page_sidebar(
  title = "Baricitinib Pediatric PK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    accordion(
      open = c("Dosing", "Patient", "Model"),
      accordion_panel(
        "Dosing",
        sliderInput("dose", "Dose (mg)", min = 0.5, max = 4, value = 4, step = 0.5),
        sliderInput("interval", "Dosing interval (h)", min = 12, max = 48, value = 24, step = 12),
        numericInput("n_days", "Duration (days)", value = 14, min = 1, max = 28)
      ),
      accordion_panel(
        "Patient",
        sliderInput("wt", "Body weight (kg)", min = 10, max = 100, value = 56.9, step = 0.1),
        checkboxInput("auto_egfr", "Auto-set eGFR to pediatric reference", value = TRUE),
        sliderInput("egfr", "Baseline eGFR (mL/min/1.73 m2)", min = 40, max = 180, value = 117, step = 5),
        sliderInput("degfr", "Change in eGFR from baseline", min = -30, max = 30, value = 0, step = 1)
      ),
      accordion_panel(
        "Model",
        sliderInput("lag_h", "Absorption lag (h)", min = 0, max = 1, value = paper_lag_h, step = 0.01),
        sliderInput("d1_h", "Zero-order absorption duration (h)", min = 0.1, max = 1, value = paper_d1_h, step = 0.01),
        checkboxInput("log_scale", "Log scale (Y-axis)", value = FALSE)
      )
    ),
    layout_column_wrap(
      width = 1 / 2,
      fill = FALSE,
      actionButton("load_under30", "Load <30 kg cohort", class = "btn btn-outline-primary btn-sm"),
      actionButton("load_over30", "Load >=30 kg cohort", class = "btn btn-outline-primary btn-sm")
    ),
    div(
      class = "ref-box",
      strong("Paper-informed regimen"),
      br(),
      textOutput("weight_band_note"),
      br(),
      textOutput("paper_reference_note"),
      br(),
      textOutput("adult_reference_note")
    )
  ),
  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", "AUCtau,ss (ng.h/mL)")),
    div(class = "metric-card metric-info",
        div(class = "metric-value", textOutput("thalf")),
        div(class = "metric-label", "t1/2 (h)"))
  ),
  card(
    full_screen = TRUE,
    plotOutput("pkPlot", height = "500px"),
    card_footer(
      textOutput("plot_caption")
    )
  ),
  navset_card_underline(
    title = "Drug Information",
    nav_panel(
      "Model Information",
      div(
        class = "ref-box",
        markdown("
## Baricitinib in Pediatric Atopic Dermatitis

**Primary paper:** Decker RL, Ernest CS II, Radtke DB, Prakash A, Zhang X. *Clinical Pharmacokinetics* (2026) 65:119-131.  
**Drug:** Baricitinib (Olumiant)  
**Class:** Oral Janus kinase (JAK)1/JAK2 inhibitor  
**Population:** Pediatric patients aged 2 to <18 years with moderate-to-severe atopic dermatitis

### Final population PK model
- Linear **2-compartment** model
- **Zero-order absorption** with **0.144 h lag time**
- Clearance partitioned into **renal** and **non-renal** components
- Allometric scaling on **clearance-related parameters (0.75)** and **volumes (1.0)**

### Final parameter estimates from Table 1
| Parameter | Estimate | Units |
|-----------|----------|-------|
| CLnr/F | 2.76 | L/h |
| CLr/F | 7.90 | L/h |
| V1/F | 119 | L |
| Q | 2.40 | L/h |
| V2/F | 46.8 | L |
| D1 | 0.263 | h |
| Lag time | 0.144 | h |

### Covariate relationships
- **CL/F = (CLnr/F + CLr/F) * (WT / 74)^0.75**
- **V1/F = 119 * (WT / 74)^1.0**
- **Q = 2.4 * (WT / 74)^0.75**
- **V2/F = 46.8 * (WT / 74)^1.0**
- **CLr/F** scales with baseline eGFR and change from baseline eGFR

### Exposure targets reported in Table 2
| Group | Dose | Cmax,ss | AUCtau,ss |
|------|------|---------|-----------|
| <30 kg | 2 mg QD | 57.1 ng/mL | 298 ng.h/mL |
| >=30 kg | 4 mg QD | 50.3 ng/mL | 383 ng.h/mL |
| Adult AD reference | 4 mg QD | 45.0 ng/mL | 380 ng.h/mL |

### Dosing conclusion from the paper
- **2 mg once daily** for pediatric patients **10 to <30 kg**
- **4 mg once daily** for pediatric patients **>=30 kg**
- The exposure-response analysis showed a higher vIGA-AD 0/1 response with higher exposure, but the commercial pediatric regimen was chosen to match adult 4 mg exposure.
        ")
      )
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Primary References"),
        tags$ol(
          tags$li(
            "Decker RL, Ernest CS II, Radtke DB, Prakash A, Zhang X. (2026). ",
            tags$em("A Population Pharmacokinetic and Exposure-Response Analysis for Baricitinib in Pediatric Patients with Atopic Dermatitis."),
            " Clinical Pharmacokinetics 65:119-131. ",
            tags$a(href = "https://doi.org/10.1007/s40262-025-01563-8", target = "_blank", "doi:10.1007/s40262-025-01563-8")
          ),
          tags$li(
            "Zhang X et al. (2017). Dose/exposure-response modeling to support phase III baricitinib dosing in rheumatoid arthritis. ",
            tags$em("CPT Pharmacometrics Syst Pharmacol.")
          ),
          tags$li(
            "Olumiant (baricitinib) prescribing information and EMA product information for current clinical context."
          )
        ),
        tags$h5("Clinical Context"),
        tags$ul(
          tags$li(tags$strong("Mechanism:"), " JAK1/JAK2 inhibition"),
          tags$li(tags$strong("Route:"), " Oral tablet"),
          tags$li(tags$strong("Typical AD regimen:"), " Once daily"),
          tags$li(tags$strong("Key modeling message:"), " Pediatric dosing was chosen to match adult 4 mg exposure while preserving the observed exposure-response trend."),
          tags$li(tags$strong("Important covariates:"), " Body weight and renal function")
        )
      )
    )
  ),
  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) {
  current_reference <- reactive({
    paper_reference(input$wt)
  })

  load_reference_scenario <- function(weight_kg, dose_mg) {
    updateSliderInput(session, "wt", value = weight_kg)
    updateSliderInput(session, "dose", value = dose_mg)
    updateSliderInput(session, "interval", value = 24)
    updateNumericInput(session, "n_days", value = 14)
    updateCheckboxInput(session, "auto_egfr", value = TRUE)
    updateSliderInput(session, "egfr", value = round(reference_egfr(weight_kg)))
    updateSliderInput(session, "degfr", value = 0)
    updateSliderInput(session, "lag_h", value = paper_lag_h)
    updateSliderInput(session, "d1_h", value = paper_d1_h)
  }

  observeEvent(list(input$wt, input$auto_egfr), {
    if (isTRUE(input$auto_egfr)) {
      updateSliderInput(
        session,
        "egfr",
        value = round(reference_egfr(input$wt))
      )
    }
  }, ignoreInit = FALSE)

  observeEvent(input$load_under30, {
    load_reference_scenario(21.5, 2)
  })

  observeEvent(input$load_over30, {
    load_reference_scenario(56.9, 4)
  })

  sim_bundle <- reactive({
    shiny::req(input$dose, input$interval, input$n_days, input$wt, input$egfr, input$degfr, input$lag_h, input$d1_h)
    run_simulation(
      dose_mg = input$dose,
      interval_h = input$interval,
      n_days = input$n_days,
      wt = input$wt,
      egfr = input$egfr,
      degfr = input$degfr,
      lag_h = input$lag_h,
      d1_h = input$d1_h
    )
  })

  metrics <- reactive({
    compute_metrics(sim_bundle())
  })

  output$weight_band_note <- renderText({
    ref <- current_reference()
    paste0(
      "Current weight band: ", ref$band,
      " • paper-recommended dose: ", ref$recommended_dose, " mg QD",
      " • published cohort median weight ", sprintf("%.1f", ref$cohort_median_weight),
      " kg (N=", ref$cohort_size, ")",
      " • auto-eGFR ", round(reference_egfr(input$wt))
    )
  })

  output$paper_reference_note <- renderText({
    ref <- current_reference()
    paste0(
      "Published cohort-median reference for ", ref$band, " (",
      sprintf("%.1f", ref$cohort_median_weight), " kg, ",
      ref$recommended_dose, " mg QD): Cmax,ss ",
      sprintf("%.1f", ref$cmax), " ng/mL, AUCtau,ss ",
      sprintf("%.0f", ref$auc), " ng.h/mL, t1/2 ",
      sprintf("%.1f", ref$half_life), " h."
    )
  })

  output$adult_reference_note <- renderText({
    adult_ref <- adult_reference()
    paste0(
      "Adult comparator from Table 2 (context only): median weight ",
      sprintf("%.1f", adult_ref$cohort_median_weight), " kg, 4 mg QD, Cmax,ss ",
      sprintf("%.1f", adult_ref$cmax), " ng/mL, AUCtau,ss ",
      sprintf("%.0f", adult_ref$auc), " ng.h/mL, t1/2 ",
      sprintf("%.1f", adult_ref$half_life),
      " h. This pediatric model is calibrated to the pediatric cohorts above and is not expected to exactly reproduce the adult row."
    )
  })

  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$plot_caption <- renderText({
    ref <- current_reference()
    met <- metrics()
    if (is_reference_scenario(input$wt, input$dose)) {
      paste0(
        "Exact Table 2 cohort-median scenario loaded for ", ref$band, " (",
        sprintf("%.1f", ref$cohort_median_weight), " kg, ",
        ref$recommended_dose, " mg QD). Published Cmax,ss ",
        sprintf("%.1f", ref$cmax), " ng/mL and AUCtau,ss ",
        sprintf("%.0f", ref$auc), " ng.h/mL. Current simulation: Cmax,ss ",
        sprintf("%.1f", met$cmax), " ng/mL and AUCtau,ss ",
        sprintf("%.0f", met$auc), " ng.h/mL."
      )
    } else {
      paste0(
        "Published cohort-median comparator for ", ref$band, " is ",
        sprintf("%.1f", ref$cohort_median_weight), " kg at ",
        ref$recommended_dose, " mg QD with Cmax,ss ",
        sprintf("%.1f", ref$cmax), " ng/mL and AUCtau,ss ",
        sprintf("%.0f", ref$auc),
        " ng.h/mL. Your current scenario is an individualized simulation within that weight band and is not expected to exactly equal the cohort summary."
      )
    }
  })

  output$pkPlot <- renderPlot({
    bundle <- sim_bundle()
    sim <- bundle$sim

    p <- ggplot(sim, aes(x = time_day, y = CP)) +
      geom_vline(
        data = bundle$events,
        aes(xintercept = time / 24),
        inherit.aes = FALSE,
        color = "#cbd5e1",
        linewidth = 0.4
      ) +
      geom_line(color = "#8b5cf6", linewidth = 1) +
      labs(
        x = "Time (days)",
        y = "Plasma concentration (ng/mL)",
        title = paste0("Baricitinib ", input$dose, " mg every ", input$interval, " h"),
        subtitle = paste0(
          "WT ", input$wt, " kg • eGFR ", input$egfr,
          " mL/min/1.73 m2 • D1 ", sprintf("%.3f", input$d1_h),
          " h • lag ", sprintf("%.3f", input$lag_h), " h"
        )
      ) +
      theme_minimal(base_size = 14) +
      theme(
        plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank()
      )

    if (isTRUE(input$log_scale)) {
      positive_sim <- sim |>
        dplyr::filter(CP > 0)
      p <- ggplot(positive_sim, aes(x = time_day, y = CP)) +
        geom_vline(
          data = bundle$events,
          aes(xintercept = time / 24),
          inherit.aes = FALSE,
          color = "#cbd5e1",
          linewidth = 0.4
        ) +
        geom_line(color = "#8b5cf6", linewidth = 1) +
        scale_y_log10() +
        labs(
          x = "Time (days)",
          y = "Plasma concentration (ng/mL, log scale)",
          title = paste0("Baricitinib ", input$dose, " mg every ", input$interval, " h"),
          subtitle = paste0(
            "WT ", input$wt, " kg • eGFR ", input$egfr,
            " mL/min/1.73 m2 • D1 ", sprintf("%.3f", input$d1_h),
            " h • lag ", sprintf("%.3f", input$lag_h), " h"
          )
        ) +
        theme_minimal(base_size = 14) +
        theme(
          plot.title = element_text(face = "bold"),
          panel.grid.minor = element_blank()
        )
    }

    p
  })
}

shinyApp(ui = ui, server = server)