use("shiny")
use("bslib")
use("mrgsolve")
use("dplyr")
use("ggplot2")

model_code <- '
$PARAM @annotated
CLF_IR   : 37.665 : Apparent clearance for immediate-release levodopa (L/h)
VF_IR    : 36.746 : Apparent volume for immediate-release levodopa (L)
KA_IR    : 3.456  : Absorption rate constant for immediate-release levodopa (1/h)
LAG_IR   : 0.00   : Absorption lag for immediate-release levodopa (h)
CLF_CR   : 49.875 : Apparent clearance for controlled-release levodopa (L/h)
VF_CR    : 76.379 : Apparent volume for controlled-release levodopa (L)
KA_CR    : 0.953  : Absorption rate constant for controlled-release levodopa (1/h)
LAG_CR   : 0.84   : Absorption lag for controlled-release levodopa (h)
ENTA_CLF : 0.85   : Relative apparent clearance with entacapone
ENTA_KA  : 0.85   : Relative absorption rate with entacapone
ENTA_LAG : 0.10   : Additional lag with entacapone (h)
FORM     : 1      : Formulation flag (1 = IR, 2 = CR)
ENTA     : 0      : Entacapone flag (0 = no, 1 = yes)
STAGE    : 2      : Hoehn and Yahr stage (1-4)

$CMT @annotated
DEPOT  : Oral depot (mg)
CENT   : Central compartment (mg)
EFFECT : Effect compartment (ug/mL)

$MAIN
double LAG_main = FORM < 1.5 ? LAG_IR : LAG_CR;
if (ENTA > 0.5) LAG_main = LAG_main + ENTA_LAG;
ALAG_DEPOT = LAG_main;

$ODE
double CLF_ode = FORM < 1.5 ? CLF_IR : CLF_CR;
double VF_ode  = FORM < 1.5 ? VF_IR  : VF_CR;
double KA_ode  = FORM < 1.5 ? KA_IR  : KA_CR;
double t12eq_ode = STAGE < 1.5 ? (133.0 / 60.0) : (STAGE < 2.5 ? (78.0 / 60.0) : (STAGE < 3.5 ? (28.0 / 60.0) : (20.0 / 60.0)));

if (ENTA > 0.5) {
  CLF_ode = CLF_ode * ENTA_CLF;
  KA_ode  = KA_ode  * ENTA_KA;
}

double KEL_ode = CLF_ode / VF_ode;
double KEO_ode = log(2.0) / t12eq_ode;
double CP_ode = CENT / VF_ode;

dxdt_DEPOT  = -KA_ode * DEPOT;
dxdt_CENT   = KA_ode * DEPOT - KEL_ode * CENT;
dxdt_EFFECT = KEO_ode * (CP_ode - EFFECT);

$TABLE
double CLFi = FORM < 1.5 ? CLF_IR : CLF_CR;
double VFi  = FORM < 1.5 ? VF_IR  : VF_CR;
double KAi  = FORM < 1.5 ? KA_IR  : KA_CR;
double LAGi = FORM < 1.5 ? LAG_IR : LAG_CR;
if (ENTA > 0.5) {
  CLFi = CLFi * ENTA_CLF;
  KAi  = KAi  * ENTA_KA;
  LAGi = LAGi + ENTA_LAG;
}

double EC50;
double EMAX;
double E0;
double HILL;
double DYSK_EC50;
double DYSK_EMAX;
double DYSK_HILL;

if (STAGE < 1.5) {
  EC50 = 0.20;
  EMAX = 45.0;
  E0 = 144.0;
  HILL = 2.0;
  DYSK_EC50 = 999.0;
  DYSK_EMAX = 0.0;
  DYSK_HILL = 1.0;
} else if (STAGE < 2.5) {
  EC50 = 0.29;
  EMAX = 35.0;
  E0 = 141.0;
  HILL = 5.0;
  DYSK_EC50 = 0.70;
  DYSK_EMAX = 1.0;
  DYSK_HILL = 14.0;
} else if (STAGE < 3.5) {
  EC50 = 0.60;
  EMAX = 41.0;
  E0 = 114.0;
  HILL = 7.0;
  DYSK_EC50 = 0.65;
  DYSK_EMAX = 3.0;
  DYSK_HILL = 27.0;
} else {
  EC50 = 0.94;
  EMAX = 55.0;
  E0 = 117.0;
  HILL = 18.0;
  DYSK_EC50 = 0.98;
  DYSK_EMAX = 3.0;
  DYSK_HILL = 30.0;
}

double KEL = CLFi / VFi;
double KEO = log(2.0) / (STAGE < 1.5 ? (133.0 / 60.0) : (STAGE < 2.5 ? (78.0 / 60.0) : (STAGE < 3.5 ? (28.0 / 60.0) : (20.0 / 60.0))));
double CP = CENT / VFi;
double CE = EFFECT > 0 ? EFFECT : 0.0;
double MOTOR = E0 + EMAX * pow(CE, HILL) / (pow(EC50, HILL) + pow(CE, HILL));
double DYSK = DYSK_EMAX > 0 ?
  DYSK_EMAX * pow(CE, DYSK_HILL) / (pow(DYSK_EC50, DYSK_HILL) + pow(CE, DYSK_HILL)) :
  0.0;
double ABOVE_BRIEF = CP >= 1.2 ? 1.0 : 0.0;
double ABOVE_ON = CP >= 1.6 ? 1.0 : 0.0;
double ABOVE_SUSTAINED = CP >= 2.0 ? 1.0 : 0.0;

$CAPTURE @annotated
CP              : Plasma levodopa concentration (ug/mL)
CE              : Effect-site levodopa concentration (ug/mL)
MOTOR           : Tapping motor effect score (arbitrary units)
DYSK            : Dyskinesia effect index (arbitrary units)
CLFi            : Apparent clearance (L/h)
VFi             : Apparent volume (L)
KAi             : Absorption rate constant (1/h)
LAGi            : Absorption lag time (h)
KEL             : Elimination rate constant (1/h)
KEO             : Effect-site equilibration rate constant (1/h)
ABOVE_BRIEF     : Above brief-response threshold
ABOVE_ON        : Above on-threshold
ABOVE_SUSTAINED : Above sustained-response threshold
'

mod <- mrgsolve::mcode("levodopa_pkpd_333", model_code)

metric_card <- function(output_id, label, class = "metric-primary") {
  div(
    class = paste("metric-card", class),
    div(class = "metric-value", textOutput(output_id)),
    div(class = "metric-label", label)
  )
}

reference_pk <- data.frame(
  formulation = c(
    "Immediate release (Sinemet 25/100 x2)",
    "Controlled release (Sinemet CR 50/200)"
  ),
  tmax_h = c(0.5, 2.1),
  cmax_ug_ml = c(3.26, 1.15),
  auc_ug_h_ml = c(5.31, 4.01),
  context = c("Healthy elderly, fasting", "Healthy elderly, fasting")
)

reference_pd <- data.frame(
  stage = c("I", "II", "III", "IV"),
  t12eq_min = c(133, 78, 28, 20),
  ec50_ug_ml = c(0.20, 0.29, 0.60, 0.94),
  emax = c(45, 35, 41, 55),
  baseline = c(144, 141, 114, 117),
  hill = c(2, 5, 7, 18)
)

dysk_reference <- data.frame(
  stage = c("I", "II", "III", "IV"),
  t12eq_min = c(NA, 22, 19, 17),
  ec50_ug_ml = c(NA, 0.70, 0.65, 0.98),
  emax = c(NA, 1, 3, 3),
  hill = c(NA, 14, 27, 30)
)

simulate_levodopa <- function(
  dose_mg = 200,
  interval_h = 4,
  n_days = 1,
  formulation = "IR",
  entacapone = FALSE,
  stage = 2,
  delta_h = 0.05
) {
  form_code <- ifelse(formulation == "CR", 2, 1)
  ent_flag <- ifelse(entacapone, 1, 0)
  n_doses <- max(1, ceiling((n_days * 24) / interval_h))
  last_dose_time <- (n_doses - 1) * interval_h
  sim_end_h <- last_dose_time + interval_h
  ev_obj <- mrgsolve::ev(amt = dose_mg, cmt = 1, ii = interval_h, addl = n_doses - 1)

  mod %>%
    mrgsolve::param(
      FORM = form_code,
      ENTA = ent_flag,
      STAGE = stage
    ) %>%
    mrgsolve::ev(ev_obj) %>%
    mrgsolve::mrgsim(end = sim_end_h, delta = delta_h) %>%
    as.data.frame() |>
    dplyr::mutate(
      time_h = time,
      last_dose_time = last_dose_time,
      interval_h = interval_h,
      dose_mg = dose_mg,
      formulation = formulation,
      entacapone = entacapone,
      stage = stage
    )
}

last_interval_metrics <- function(data) {
  last_start <- data$last_dose_time[[1]]
  interval_h <- data$interval_h[[1]]
  last_window <- data |>
    dplyr::filter(time_h >= last_start, time_h <= last_start + interval_h)

  cmax <- max(last_window$CP, na.rm = TRUE)
  ctrough <- tail(last_window$CP, 1)
  tmax_h <- last_window$time_h[[which.max(last_window$CP)]] - last_start
  auc_tau <- sum(
    diff(last_window$time_h) *
      (head(last_window$CP, -1) + tail(last_window$CP, -1)) / 2
  )
  motor_peak <- max(last_window$MOTOR, na.rm = TRUE)
  dysk_peak <- max(last_window$DYSK, na.rm = TRUE)
  above_16_h <- sum(diff(last_window$time_h) * head(last_window$ABOVE_ON, -1))
  above_20_h <- sum(diff(last_window$time_h) * head(last_window$ABOVE_SUSTAINED, -1))
  thalf_h <- log(2) / last_window$KEL[[1]]

  list(
    cmax = cmax,
    ctrough = ctrough,
    tmax_h = tmax_h,
    auc_tau = auc_tau,
    thalf_h = thalf_h,
    motor_peak = motor_peak,
    dysk_peak = dysk_peak,
    above_16_h = above_16_h,
    above_20_h = above_20_h,
    lag_h = last_window$LAGi[[1]],
    ka_h = last_window$KAi[[1]],
    clf_h = last_window$CLFi[[1]],
    vf_l = last_window$VFi[[1]]
  )
}

stage_label <- function(stage) {
  c("I", "II", "III", "IV")[stage]
}

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;
      min-height: 100px;
    }
    .metric-value {
      font-size: 24px;
      font-weight: bold;
      color: #2c3e50;
    }
    .metric-label {
      font-size: 12px;
      color: #7f8c8d;
      text-transform: uppercase;
      letter-spacing: 0.04em;
    }
    .metric-success .metric-value { color: #10b981; }
    .metric-warning .metric-value { color: #f59e0b; }
    .metric-primary .metric-value { color: #8b5cf6; }
    .metric-info .metric-value { color: #0dcaf0; }
    .summary-box {
      background: #f8fafc;
      border: 1px solid #e2e8f0;
      border-radius: 8px;
      padding: 14px 16px;
      margin-top: 12px;
      font-size: 13px;
      line-height: 1.6;
    }
    .ref-box {
      background: #f0f4ff;
      border-left: 4px solid #8b5cf6;
      padding: 12px 16px;
      border-radius: 4px;
      margin-top: 10px;
      font-size: 13px;
    }
    .note-box {
      background: #fff7ed;
      border-left: 4px solid #f59e0b;
      padding: 12px 16px;
      border-radius: 4px;
      margin-top: 14px;
      font-size: 13px;
    }
  ")

ui <- page_sidebar(
  title = "Levodopa PK/PD Simulator",
  theme = app_theme,
  sidebar = sidebar(
    title = "Simulation Settings",
    width = 340,
    accordion(
      accordion_panel(
        "Dosing",
        selectInput(
          "formulation",
          "Formulation",
          choices = c(
            "Immediate release (IR)" = "IR",
            "Controlled release (CR)" = "CR"
          ),
          selected = "IR"
        ),
        sliderInput("dose", "Levodopa dose (mg)", min = 50, max = 400, value = 200, step = 25),
        sliderInput("interval", "Dosing interval (h)", min = 3, max = 8, value = 4, step = 1),
        numericInput("n_days", "Simulation duration (days)", value = 1, min = 1, max = 7)
      ),
      accordion_panel(
        "Disease State",
        selectInput(
          "stage",
          "Hoehn & Yahr stage",
          choices = c(
            "Stage I" = 1,
            "Stage II" = 2,
            "Stage III" = 3,
            "Stage IV" = 4
          ),
          selected = 2
        ),
        checkboxInput("entacapone", "Add entacapone effect", value = FALSE),
        checkboxInput("log_scale", "Log scale PK plot", value = FALSE)
      ),
      accordion_panel(
        "Interpretation",
        div(
          class = "summary-box",
          HTML(
            paste(
              "<strong>Published thresholds from the review:</strong><br>",
              "Below 1.2 ug/mL: often no or only brief response.<br>",
              "Around 1.6 ug/mL: typical switch to 'on' in fluctuating patients.<br>",
              "Around 2.0 ug/mL: sustained response during infusion studies."
            )
          )
        )
      ),
      open = c("Dosing", "Disease State")
    )
  ),
  navset_card_underline(
    id = "main_tabs",
    title = "Levodopa Simulator",
    nav_panel(
      "Simulation",
      layout_column_wrap(
        width = 1 / 4,
        fill = FALSE,
        metric_card("cmax", "Cmax (ug/mL)", "metric-success"),
        metric_card("ctrough", "Ctrough (ug/mL)", "metric-warning"),
        metric_card("auc", "AUC0-tau (ug*h/mL)", "metric-primary"),
        metric_card("thalf", "Apparent t1/2 (h)", "metric-info")
      ),
      card(
        full_screen = TRUE,
        card_header("Plasma Levodopa Concentration-Time Profile"),
        plotOutput("pkPlot", height = "500px")
      ),
      card(
        card_header("Motor Effect and Dyskinesia Index"),
        plotOutput("pdPlot", height = "320px")
      ),
      div(class = "summary-box", htmlOutput("simSummary"))
    ),
    nav_panel(
      "Model Information",
      div(
        class = "summary-box",
        tags$h5("Model scope"),
        tags$p(
          "This simulator reproduces the review by Khor and Hsu (2007) using a",
          " formulation-specific oral one-compartment PK model linked to a stage-specific",
          " effect-compartment PD model."
        ),
        tags$ul(
          tags$li("IR and CR oral PK were calibrated to published Tmax, Cmax, and AUC values from healthy elderly subjects under fasting conditions."),
          tags$li("PD uses the reported effect-site equilibration half-life, EC50, Emax, baseline effect, and Hill coefficient by Hoehn & Yahr stage."),
          tags$li("The entacapone toggle is qualitative: it increases late exposure and slows early absorption, consistent with the review narrative."),
          tags$li("The displayed half-life is the model-derived apparent oral disposition half-life for the selected formulation.")
        )
      ),
      card(
        card_header("Formulation PK Anchors from the Review"),
        tableOutput("pkRefTable")
      ),
      card(
        card_header("Stage-Specific Motor Response Parameters"),
        tableOutput("pdRefTable")
      ),
      card(
        card_header("Stage-Specific Dyskinesia Parameters"),
        tableOutput("dyskRefTable")
      ),
      div(
        class = "note-box",
        tags$strong("Why this is a simplified build:"),
        tags$p(
          "The manuscript is a pharmacology review, not a full population PK report with reusable raw code. ",
          "This app therefore prioritizes published exposure and response landmarks rather than pretending to reconstruct unreported compartments."
        )
      )
    ),
    nav_panel(
      "References",
      div(
        class = "ref-box",
        tags$h5("Primary source"),
        tags$ol(
          tags$li("Khor SP, Hsu A. The Pharmacokinetics and Pharmacodynamics of Levodopa in the Treatment of Parkinson's Disease. Current Clinical Pharmacology. 2007;2:234-243.")
        ),
        tags$h5("Key figures and findings encoded here"),
        tags$ul(
          tags$li("Table 1: IR vs CR Tmax, Cmax, and AUC in healthy elderly subjects under fasting conditions."),
          tags$li("Table 2: Stage-specific t1/2eq, EC50, Emax, E0, and Hill coefficient for tapping motor effect and dyskinesia."),
          tags$li("Section 3.3: Optimal clinical response near 2.6-2.8 ug/mL with or without carbidopa."),
          tags$li("Section 3.5: Approximate threshold concentration near 1.6 ug/mL in fluctuating patients.")
        ),
        tags$h5("Clinical context"),
        tags$ul(
          tags$li(tags$strong("Drug class:"), " Dopamine precursor / antiparkinsonian agent"),
          tags$li(tags$strong("Indication:"), " Motor symptom control in Parkinson's disease"),
          tags$li(tags$strong("Route modeled:"), " Oral carbidopa/levodopa immediate- and controlled-release formulations"),
          tags$li(tags$strong("Use case:"), " Comparing how formulation choice and disease stage shift time above response thresholds")
        )
      )
    )
  ),
  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$formulation, input$stage)
    simulate_levodopa(
      dose_mg = input$dose,
      interval_h = input$interval,
      n_days = input$n_days,
      formulation = input$formulation,
      entacapone = input$entacapone,
      stage = as.numeric(input$stage)
    )
  })

  sim_metrics <- reactive({
    last_interval_metrics(sim_data())
  })

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

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

  output$auc <- renderText({
    sprintf("%.2f", sim_metrics()$auc_tau)
  })

  output$thalf <- renderText({
    sprintf("%.2f", sim_metrics()$thalf_h)
  })

  output$pkPlot <- renderPlot({
    d <- sim_data()
    if (input$log_scale) {
      d <- d |> dplyr::filter(CP > 0)
    }

    peak_x <- max(d$time_h, na.rm = TRUE) * 0.02

    p <- ggplot(d, aes(x = time_h, y = CP)) +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = 1.2, fill = "#fee2e2", alpha = 0.6) +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 1.2, ymax = 1.6, fill = "#fef3c7", alpha = 0.5) +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = 1.6, ymax = 2.0, fill = "#dcfce7", alpha = 0.45) +
      geom_hline(yintercept = c(1.2, 1.6, 2.0), linetype = "dashed", color = c("#ef4444", "#f59e0b", "#10b981")) +
      annotate("text", x = peak_x, y = 1.23, label = "Brief / no response", hjust = 0, vjust = 0, size = 3.2, color = "#b91c1c") +
      annotate("text", x = peak_x, y = 1.63, label = "Switch-on threshold", hjust = 0, vjust = 0, size = 3.2, color = "#b45309") +
      annotate("text", x = peak_x, y = 2.03, label = "Sustained response", hjust = 0, vjust = 0, size = 3.2, color = "#047857") +
      geom_line(color = "#8b5cf6", linewidth = 0.9) +
      labs(
        x = "Time (hours)",
        y = "Levodopa concentration (ug/mL)",
        title = paste0(
          ifelse(input$formulation == "IR", "Immediate-release", "Controlled-release"),
          " levodopa ", input$dose, " mg Q", input$interval, "H"
        ),
        subtitle = ifelse(
          input$entacapone,
          "Entacapone enabled: higher late exposure, slower early absorption",
          paste0("Hoehn & Yahr stage ", stage_label(as.numeric(input$stage)))
        )
      ) +
      theme_minimal(base_size = 14)

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

    p
  })

  output$pdPlot <- renderPlot({
    d <- sim_data()
    pd_df <- rbind(
      data.frame(time_h = d$time_h, value = d$MOTOR, endpoint = "Motor effect"),
      data.frame(time_h = d$time_h, value = d$DYSK, endpoint = "Dyskinesia index")
    )

    ggplot(pd_df, aes(x = time_h, y = value, color = endpoint)) +
      geom_line(linewidth = 0.9) +
      facet_wrap(~endpoint, ncol = 1, scales = "free_y") +
      scale_color_manual(values = c("Motor effect" = "#8b5cf6", "Dyskinesia index" = "#ef4444")) +
      labs(
        x = "Time (hours)",
        y = "Predicted effect",
        title = paste0("Stage ", stage_label(as.numeric(input$stage)), " pharmacodynamic response")
      ) +
      theme_minimal(base_size = 13) +
      theme(legend.position = "none")
  })

  output$simSummary <- renderUI({
    m <- sim_metrics()
    formulation_name <- ifelse(input$formulation == "IR", "Immediate-release", "Controlled-release")
    reference_row <- reference_pk[reference_pk$formulation == ifelse(
      input$formulation == "IR",
      "Immediate release (Sinemet 25/100 x2)",
      "Controlled release (Sinemet CR 50/200)"
    ), ]
    anchor_data <- simulate_levodopa(
      dose_mg = 200,
      interval_h = 24,
      n_days = 1,
      formulation = input$formulation,
      entacapone = FALSE,
      stage = 2
    )
    anchor_metrics <- last_interval_metrics(anchor_data)

    HTML(
      paste0(
        "<strong>", formulation_name, " summary:</strong> ",
        "Simulated Tmax is ", sprintf("%.2f", m$tmax_h), " h with Cmax ",
        sprintf("%.2f", m$cmax), " ug/mL and AUC0-tau ", sprintf("%.2f", m$auc_tau),
        " ug*h/mL. ",
        "Time above 1.6 ug/mL is ", sprintf("%.2f", m$above_16_h), " h; time above 2.0 ug/mL is ",
        sprintf("%.2f", m$above_20_h), " h. ",
        "Peak motor-effect score is ", sprintf("%.1f", m$motor_peak),
        " and peak dyskinesia index is ", sprintf("%.2f", m$dysk_peak), ".<br><br>",
        "<strong>Single-dose calibration check (200 mg, fasting, no entacapone):</strong> model Tmax ",
        sprintf("%.2f", anchor_metrics$tmax_h), " h, model Cmax ",
        sprintf("%.2f", anchor_metrics$cmax), " ug/mL, model AUC ",
        sprintf("%.2f", anchor_metrics$auc_tau), " ug*h/mL. ",
        "Published anchor: Tmax ",
        sprintf("%.2f", reference_row$tmax_h), " h, Cmax ",
        sprintf("%.2f", reference_row$cmax_ug_ml), " ug/mL, AUC ",
        sprintf("%.2f", reference_row$auc_ug_h_ml), " ug*h/mL.<br>",
        "<strong>Model parameters in use:</strong> CL/F ", sprintf("%.2f", m$clf_h),
        " L/h, V/F ", sprintf("%.2f", m$vf_l), " L, Ka ", sprintf("%.2f", m$ka_h),
        " 1/h, lag ", sprintf("%.2f", m$lag_h), " h."
      )
    )
  })

  output$pkRefTable <- renderTable({
    reference_pk
  }, striped = TRUE, bordered = TRUE, digits = 2)

  output$pdRefTable <- renderTable({
    reference_pd
  }, striped = TRUE, bordered = TRUE, digits = 2, na = "")

  output$dyskRefTable <- renderTable({
    dysk_reference
  }, striped = TRUE, bordered = TRUE, digits = 2, na = "")
}

shinyApp(ui = ui, server = server)