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

model_code <- '
$PARAM @annotated
CLPOP     : 22.0   : Typical clearance (L/h)
VCPOP     : 87.8   : Typical central volume (L)
QPOP      : 1.16   : Typical intercompartmental clearance (L/h)
VPPOP     : 8.25   : Typical peripheral volume (L)
KAPOP     : 4.01   : First-order absorption rate constant (1/h)
WT        : 70     : Body weight (kg)
WT_CL_Q   : 0.453  : Weight exponent on CL and Q
WT_V      : 0.52   : Weight exponent on Vc and Vp
CL_RIF    : 1.0    : Rifampin effect on CL
CL_FLZ    : 1.0    : Fluconazole effect on CL
CL_FLU    : 1.0    : Fluvoxamine effect on CL
CL_TAFO_D : -0.186 : Max time effect on CL
CL_TAFO_H : 21.6   : Half-life of CL change (h)
EFFDOSR   : 200    : Effective daily dose (mg/day)
CL_DOSE_E : -0.169 : Effective daily dose exponent on CL

$CMT @annotated
DEPOT  : First-order absorption depot (mg)
CENT   : Central compartment (mg)
PERIPH : Peripheral compartment (mg)

$MAIN
double wt_scale_cl = pow(WT / 70.0, WT_CL_Q);
double wt_scale_v  = pow(WT / 70.0, WT_V);
double time_factor = 1.0 + CL_TAFO_D * (1.0 - exp(-1.0 * (log(2.0) / CL_TAFO_H) * TIME));
double dose_factor = pow(EFFDOSR / 200.0, CL_DOSE_E);

double CLi = CLPOP * wt_scale_cl * CL_RIF * CL_FLZ * CL_FLU * time_factor * dose_factor;
double VCi = VCPOP * wt_scale_v;
double Qi  = QPOP * wt_scale_cl;
double VPi = VPPOP * wt_scale_v;

double K10 = CLi / VCi;
double K12 = Qi / VCi;
double K21 = Qi / VPi;

$ODE
dxdt_DEPOT  = -KAPOP * DEPOT;
dxdt_CENT   = KAPOP * DEPOT - K10 * CENT - K12 * CENT + K21 * PERIPH;
dxdt_PERIPH = K12 * CENT - K21 * PERIPH;

$TABLE
double CP_NG_ML = (CENT / VCi) * 1000.0;

$CAPTURE CP_NG_ML CLi VCi Qi VPi
'

mod <- mcode("abrocitinib_pop_pk", model_code)

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

inv_logit <- function(x) {
  exp(x) / (1 + exp(x))
}

calc_abrocitinib_inputs <- function(
  dose_mg,
  interval_h,
  n_days,
  wt,
  sex,
  age_group,
  race_group,
  population,
  formulation,
  ddi,
  fed
) {
  popf <- 0.5977
  popft <- log(popf / (1 - popf))
  multiple_dosing <- as.integer((n_days * 24 / interval_h) > 1)

  cov_rif_f <- if (ddi == "rifampin") -2.08 else 0
  cov_cypinh_f <- if (ddi %in% c("fluconazole", "fluvoxamine")) 1.31 else 0
  cov_forms_f <- switch(
    formulation,
    "phase_ii_100" = 0,
    "phase_ii_b" = -1.02,
    "phase_iii" = -0.766,
    "suspension" = 0,
    0
  )
  cov_multi_f <- if (multiple_dosing == 1) 0.241 else 0
  cov_race_f <- if (race_group == "asian_other") 0.815 else 0
  cov_ptst_f <- switch(
    population,
    "healthy" = 0,
    "psoriasis" = 0.489,
    "ad" = 0.489,
    "hepatic_mild" = 1.3,
    "hepatic_moderate" = 1.3,
    0
  )
  cov_age_f <- if (age_group == "adolescent") -0.589 else 0
  cov_dose_f <- if (dose_mg >= 800) -0.778 else 0
  cov_sex_f <- if (sex == "female") 0.353 else 0

  fi <- inv_logit(
    popft + cov_rif_f + cov_cypinh_f + cov_forms_f + cov_multi_f +
      cov_race_f + cov_ptst_f + cov_age_f + cov_dose_f + cov_sex_f
  )

  cov_food_ak1 <- if (fed) 0 else 1
  cov_forms_ak1 <- switch(
    formulation,
    "suspension" = 1 + 1.17,
    "phase_ii_b" = 1 - 0.68,
    "phase_ii_100" = 1,
    "phase_iii" = 1,
    1
  )
  ak1_mg <- 121 * cov_food_ak1 * cov_forms_ak1
  first_order_fraction <- min(1, ak1_mg / dose_mg)
  alag1 <- if (formulation == "suspension") 0 else 0.183
  alag2 <- alag1

  eff_dose_daily <- fi * dose_mg * 24 / interval_h

  cl_rif <- if (ddi == "rifampin") 1 + 0.264 else 1
  cl_flz <- if (ddi == "fluconazole") 1 - 0.541 else 1
  cl_flu <- if (ddi == "fluvoxamine") 1 - 0.234 else 1

  list(
    fi = fi,
    first_order_fraction = first_order_fraction,
    zero_order_fraction = 1 - first_order_fraction,
    alag1 = alag1,
    alag2 = alag2,
    eff_dose_daily = eff_dose_daily,
    cl_rif = cl_rif,
    cl_flz = cl_flz,
    cl_flu = cl_flu,
    wt = wt
  )
}

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

  first_order_amt <- dose_mg * cov_inputs$fi * cov_inputs$first_order_fraction
  zero_order_amt <- dose_mg * cov_inputs$fi * cov_inputs$zero_order_fraction

  depot_events <- data.frame(
    ID = 1,
    time = dose_times + cov_inputs$alag1,
    amt = first_order_amt,
    cmt = 1,
    evid = 1,
    rate = 0
  )

  central_events <- data.frame(
    ID = 1,
    time = dose_times + cov_inputs$alag2,
    amt = zero_order_amt,
    cmt = 2,
    evid = 1,
    rate = 75.3
  )

  bind_rows(depot_events, central_events) |>
    arrange(time, cmt)
}

calc_terminal_half_life <- function(cl_value, vc_value, q_value, vp_value) {
  k10 <- cl_value / vc_value
  k12 <- q_value / vc_value
  k21 <- q_value / vp_value
  beta <- 0.5 * ((k10 + k12 + k21) - sqrt((k10 + k12 + k21)^2 - 4 * k21 * k10))
  log(2) / beta
}

simulate_abrocitinib <- function(
  dose_mg = 200,
  interval_h = 24,
  n_days = 7,
  wt = 80,
  sex = "male",
  age_group = "adult",
  race_group = "western",
  population = "ad",
  formulation = "phase_iii",
  ddi = "none",
  fed = FALSE,
  delta = 0.25
) {
  cov_inputs <- calc_abrocitinib_inputs(
    dose_mg = dose_mg,
    interval_h = interval_h,
    n_days = n_days,
    wt = wt,
    sex = sex,
    age_group = age_group,
    race_group = race_group,
    population = population,
    formulation = formulation,
    ddi = ddi,
    fed = fed
  )

  regimen <- build_regimen(
    dose_mg = dose_mg,
    interval_h = interval_h,
    n_days = n_days,
    cov_inputs = cov_inputs
  )

  sim <- mod %>%
    param(
      WT = wt,
      EFFDOSR = cov_inputs$eff_dose_daily,
      CL_RIF = cov_inputs$cl_rif,
      CL_FLZ = cov_inputs$cl_flz,
      CL_FLU = cov_inputs$cl_flu
    ) %>%
    data_set(regimen) %>%
    mrgsim(end = n_days * 24, delta = delta) %>%
    as.data.frame() |>
    mutate(
      time_h = time,
      cp_ng_ml = CP_NG_ML
    )

  attr(sim, "cov_inputs") <- cov_inputs
  sim
}

calc_auc <- function(time_vec, conc_vec) {
  sum(diff(time_vec) * (head(conc_vec, -1) + tail(conc_vec, -1)) / 2)
}

calc_metrics <- function(sim_data, interval_h) {
  end_time <- max(sim_data$time_h, na.rm = TRUE)
  window_start <- max(0, end_time - interval_h)
  last_window <- sim_data |>
    dplyr::filter(time_h >= window_start, time_h <= end_time)

  cl_end <- tail(sim_data$CLi, 1)
  vc_end <- tail(sim_data$VCi, 1)
  q_end <- tail(sim_data$Qi, 1)
  vp_end <- tail(sim_data$VPi, 1)

  list(
    cmax = max(last_window$cp_ng_ml, na.rm = TRUE),
    ctrough = min(last_window$cp_ng_ml, na.rm = TRUE),
    auc = calc_auc(last_window$time_h, last_window$cp_ng_ml),
    thalf = calc_terminal_half_life(cl_end, vc_end, q_end, vp_end)
  )
}

reference_scenarios <- list(
  custom = NULL,
  ref_healthy = list(dose = 200, interval = 24, days = 7, wt = 70, sex = "male", age = "adult", race = "western", population = "healthy", formulation = "phase_iii", ddi = "none", fed = FALSE),
  adult_ad = list(dose = 200, interval = 24, days = 7, wt = 80, sex = "male", age = "adult", race = "western", population = "ad", formulation = "phase_iii", ddi = "none", fed = FALSE),
  asian_ad = list(dose = 200, interval = 24, days = 7, wt = 66, sex = "male", age = "adult", race = "asian_other", population = "ad", formulation = "phase_iii", ddi = "none", fed = FALSE),
  adolescent_ad = list(dose = 200, interval = 24, days = 7, wt = 61, sex = "male", age = "adolescent", race = "western", population = "ad", formulation = "phase_iii", ddi = "none", fed = FALSE),
  low_weight_adolescent = list(dose = 200, interval = 24, days = 7, wt = 25, sex = "male", age = "adolescent", race = "western", population = "ad", formulation = "phase_iii", ddi = "none", fed = FALSE)
)

ui <- page_sidebar(
  title = "Abrocitinib PK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    selectInput(
      "preset",
      "Scenario Preset",
      choices = c(
        "Custom" = "custom",
        "Reference Healthy 70 kg" = "ref_healthy",
        "Adult AD 80 kg" = "adult_ad",
        "Asian AD 66 kg" = "asian_ad",
        "Adolescent AD 61 kg" = "adolescent_ad",
        "Adolescent 25 kg" = "low_weight_adolescent"
      ),
      selected = "adult_ad"
    ),
    accordion(
      accordion_panel(
        "Dosing",
        sliderInput("dose", "Dose (mg)", min = 25, max = 800, value = 200, step = 25),
        numericInput("interval", "Dosing Interval (h)", value = 24, min = 12, max = 48, step = 12),
        numericInput("n_days", "Duration (days)", value = 7, min = 1, max = 14, step = 1),
        selectInput(
          "formulation",
          "Formulation",
          choices = c(
            "Phase III tablet" = "phase_iii",
            "Phase II 100 mg tablet" = "phase_ii_100",
            "Phase IIb tablets" = "phase_ii_b",
            "Oral suspension" = "suspension"
          ),
          selected = "phase_iii"
        ),
        checkboxInput("fed", "High-fat meal", value = FALSE)
      ),
      accordion_panel(
        "Patient Factors",
        sliderInput("wt", "Weight (kg)", min = 25, max = 120, value = 80, step = 1),
        selectInput("sex", "Sex", choices = c("Male" = "male", "Female" = "female"), selected = "male"),
        selectInput("age_group", "Age Group", choices = c("Adult" = "adult", "Adolescent" = "adolescent"), selected = "adult"),
        selectInput("race_group", "Race", choices = c("Western (White/Black/Unknown)" = "western", "Asian / Other" = "asian_other"), selected = "western"),
        selectInput(
          "population",
          "Population",
          choices = c(
            "Healthy volunteer" = "healthy",
            "Atopic dermatitis" = "ad",
            "Psoriasis" = "psoriasis",
            "Mild hepatic impairment" = "hepatic_mild",
            "Moderate hepatic impairment" = "hepatic_moderate"
          ),
          selected = "ad"
        ),
        selectInput(
          "ddi",
          "Drug Interaction",
          choices = c(
            "None" = "none",
            "Rifampin" = "rifampin",
            "Fluconazole" = "fluconazole",
            "Fluvoxamine" = "fluvoxamine"
          ),
          selected = "none"
        )
      ),
      accordion_panel(
        "Display",
        checkboxInput("log_scale", "Log scale", value = FALSE)
      ),
      open = c("Dosing", "Patient Factors")
    ),
    div(
      class = "small-note",
      "Model based on Wojciechowski et al. Clin Pharmacokinet 2022;61:709-723.",
      br(),
      "Key publication anchor: White adult AD, 80 kg, 200 mg QD -> Cmax 1123 ng/mL; AUC24 5662 ng·h/mL."
    )
  ),
  navset_card_underline(
    title = "Abrocitinib Population PK",
    nav_panel(
      "Simulation",
      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 (ng/mL)")),
        div(class = "metric-card metric-warning",
            div(class = "metric-value", textOutput("ctrough")),
            div(class = "metric-label", "Ctrough (ng/mL)")),
        div(class = "metric-card metric-primary",
            div(class = "metric-value", textOutput("auc")),
            div(class = "metric-label", "AUC0-24 (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,
        card_header("Plasma Concentration-Time Profile"),
        plotOutput("pk_plot", height = "500px"),
        card_footer(textOutput("scenario_summary"))
      )
    ),
    nav_panel(
      "Model Information",
      div(
        class = "ref-box",
        tags$h5("Model Summary"),
        tags$p("Two-compartment population PK model with parallel zero- and first-order absorption, time-dependent clearance, and covariate effects on bioavailability and clearance."),
        tags$ul(
          tags$li("Typical CL: 22 L/h; Vc: 87.8 L; Q: 1.16 L/h; Vp: 8.25 L"),
          tags$li("Absorption: k0 75.3 mg/h, AK1 121 mg, ka 4.01 h^-1, tablet lag 0.183 h"),
          tags$li("Fixed absolute bioavailability anchor for oral suspension: 0.5977"),
          tags$li("Population and special-population effects are carried on the manuscript's final model coefficients")
        ),
        tags$h5("Publication Anchors"),
        tags$ul(
          tags$li("White adult AD, 80 kg, 200 mg QD: Cmax 1123 ng/mL; AUC24 5662 ng·h/mL"),
          tags$li("Asian adult AD, 66 kg vs White adult AD, 80 kg: Cmax ratio 1.43; AUC ratio 1.48"),
          tags$li("Adolescent AD vs adult AD after weight adjustment: exposure overall similar; typical adolescent lower by about 14% for Cmax and 19% for AUC"),
          tags$li("25 kg adolescent vs 70 kg adult: about 27% higher exposure")
        ),
        tags$h5("Implementation Notes"),
        tags$ul(
          tags$li("Hepatic impairment is implemented as the dedicated special-population scenario reported in the paper, not stacked additively on top of AD."),
          tags$li("Food effect is modeled through the manuscript's fixed suppression of the first-order absorption component under high-fat meal conditions."),
          tags$li("Dose-dependent clearance uses effective daily dose, so non-Q24 exploratory schedules are model extrapolations rather than direct paper scenarios.")
        )
      )
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Primary Reference"),
        tags$ol(
          tags$li("Wojciechowski J, Nicholas T, Malhotra BK, Wang X, Fostvedt L, Valdez H. Population Pharmacokinetics of Abrocitinib in Healthy Individuals and Patients with Psoriasis or Atopic Dermatitis. Clinical Pharmacokinetics. 2022;61:709-723.")
        ),
        tags$h5("Supporting Context from the Paper"),
        tags$ul(
          tags$li("Steady state was achieved by day 4 after once-daily dosing."),
          tags$li("The final model used NONMEM with a two-compartment disposition, parallel zero- and first-order absorption, time-dependent F, and time- and dose-dependent CL."),
          tags$li("The paper reports exposure shifts with race, age, hepatic impairment, and CYP-mediated drug-drug interactions.")
        )
      )
    )
  ),
  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",
    br(),
    tags$span(style = "font-size: 10px;", "For research and educational purposes only. Not for clinical decision-making.")
  )
)

server <- function(input, output, session) {
  observeEvent(input$preset, {
    preset <- reference_scenarios[[input$preset]]
    if (is.null(preset)) {
      return()
    }
    updateSliderInput(session, "dose", value = preset$dose)
    updateNumericInput(session, "interval", value = preset$interval)
    updateNumericInput(session, "n_days", value = preset$days)
    updateSliderInput(session, "wt", value = preset$wt)
    updateSelectInput(session, "sex", selected = preset$sex)
    updateSelectInput(session, "age_group", selected = preset$age)
    updateSelectInput(session, "race_group", selected = preset$race)
    updateSelectInput(session, "population", selected = preset$population)
    updateSelectInput(session, "formulation", selected = preset$formulation)
    updateSelectInput(session, "ddi", selected = preset$ddi)
    updateCheckboxInput(session, "fed", value = preset$fed)
  }, ignoreInit = TRUE)

  sim_data <- reactive({
    shiny::req(input$dose, input$interval, input$n_days, input$wt)
    simulate_abrocitinib(
      dose_mg = input$dose,
      interval_h = input$interval,
      n_days = input$n_days,
      wt = input$wt,
      sex = input$sex,
      age_group = input$age_group,
      race_group = input$race_group,
      population = input$population,
      formulation = input$formulation,
      ddi = input$ddi,
      fed = isTRUE(input$fed)
    )
  })

  metrics <- reactive({
    calc_metrics(sim_data(), input$interval)
  })

  output$cmax <- renderText(sprintf("%.0f", metrics()$cmax))
  output$ctrough <- renderText(sprintf("%.0f", metrics()$ctrough))
  output$auc <- renderText(sprintf("%.0f", metrics()$auc))
  output$thalf <- renderText(sprintf("%.1f", metrics()$thalf))

  output$scenario_summary <- renderText({
    cov_inputs <- attr(sim_data(), "cov_inputs")
    sprintf(
      "Estimated absolute F: %.3f | First-order fraction: %.2f | Effective daily dose: %.1f mg/day | Steady-state is typically reached by day 4 in the publication.",
      cov_inputs$fi,
      cov_inputs$first_order_fraction,
      cov_inputs$eff_dose_daily
    )
  })

  output$pk_plot <- renderPlot({
    plot_df <- sim_data()
    if (isTRUE(input$log_scale)) {
      plot_df <- plot_df |>
        dplyr::filter(cp_ng_ml > 0)
    }

    p <- ggplot(plot_df, aes(x = time_h, y = cp_ng_ml)) +
      geom_vline(xintercept = 96, linetype = "dotted", color = "#64748b") +
      geom_line(color = "#8b5cf6", linewidth = 0.9) +
      annotate("text", x = 96, y = max(plot_df$cp_ng_ml, na.rm = TRUE) * 0.92, label = "Day 4", hjust = -0.05, color = "#64748b", size = 3.5) +
      labs(
        x = "Time (hours)",
        y = "Plasma Concentration (ng/mL)",
        title = paste0("Abrocitinib ", input$dose, " mg every ", input$interval, " h"),
        subtitle = "Final population PK model with structural and covariate terms from Wojciechowski et al."
      ) +
      theme_minimal(base_size = 14)

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

    p
  })
}

app <- shinyApp(ui = ui, server = server)
app