source("../../R/load_all.R", chdir = TRUE)

test_that("linear-only sub-model (Vmax=0) matches the analytic 2-cpt beta-phase half-life (~18.8 days)", {
  # Isolates the published CL/Q/Vc/Vp from the MM/TMDD term so the linear
  # disposition parameters alone can be checked against closed-form 2-cpt
  # theory -- see model-spec.md Section 3.1 for why the *full* model (with
  # TMDD) decays faster than this once concentration falls near Km.
  CL <- 0.334; Q <- 0.859; Vc <- 3.94; Vp <- 4.21
  k10 <- CL / Vc; k12 <- Q / Vc; k21 <- Q / Vp
  b <- k10 + k12 + k21
  beta <- (b - sqrt(b^2 - 4 * k10 * k21)) / 2
  analytic_halflife <- log(2) / beta
  expect_true(analytic_halflife > 15 & analytic_halflife < 22)

  # Single 1400 mg bolus at t=0 (real infusions are ~60-90 min, negligible
  # next to the ~19-day half-life -- see R/simulate.R's dosing note).
  rhs <- function(t, y, p) {
    pk <- pk_2cmt_mm_deriv(y["Ac"], y["Ap"], CL, Q, Vc, Vp, 0, 1.5)
    list(c(dAc = pk$dAc, dAp = pk$dAp))
  }
  times <- seq(0, 150, by = 1)
  sim <- as.data.frame(deSolve::ode(c(Ac = 1400, Ap = 0), times, rhs, NULL, method = "lsoda"))
  sim$C <- sim$Ac / Vc
  fit <- lm(log(C) ~ time, data = sim[sim$time >= 100 & sim$C > 0, ])
  simulated_halflife <- log(2) / -coef(fit)[["time"]]
  expect_equal(simulated_halflife, analytic_halflife, tolerance = 0.05)
})

test_that("full model (with TMDD) decays faster than the linear-only sub-model once C approaches Km", {
  # Known, documented behavior of the parallel linear+MM approximation
  # (model-spec.md Section 8): NOT yet validated against real observed
  # trough concentrations -- flagged as an open calibration item.
  sol <- ted_simulate("teprotumumab", sim_days = 300, dt = 1)
  late <- sol[sol$time >= 200 & sol$Conc_ugmL > 0, ]
  fit <- lm(log(Conc_ugmL) ~ time, data = late)
  full_halflife <- log(2) / -coef(fit)[["time"]]
  expect_true(full_halflife < 15)
})

test_that("teprotumumab reaches near-full IGF-1R occupancy at clinical doses", {
  sol <- ted_simulate("teprotumumab", sim_days = 150, dt = 1)
  peak_conc <- max(sol$Conc_ugmL)
  occ <- occupancy(peak_conc, ted_drugs$teprotumumab$pk$Km)
  expect_true(occ > 0.9)
})

test_that("IGF-1 rises 2-4x over baseline under sustained IGF-1R blockade", {
  sol <- ted_simulate("teprotumumab", sim_days = 150, dt = 1)
  expect_true(max(sol$IGF1) > 2 & max(sol$IGF1) < 4.5)
})

test_that("no-treatment natural history shows disease burnout over time", {
  sol <- ted_simulate("none", trab_baseline = 1, disease_duration_days = 0,
                       tau_burnout = 730, sim_days = 1000, dt = 10)
  expect_true(sol$CAS[sol$time == 1000] < sol$CAS[sol$time == 0])
})

test_that("treatment reduces CAS and proptosis vs no treatment at day 150", {
  treated <- ted_simulate("teprotumumab", sim_days = 150, dt = 1)
  untreated <- ted_simulate("none", sim_days = 150, dt = 1)
  d <- 150
  expect_true(treated$CAS[treated$time == d] < untreated$CAS[untreated$time == d])
  expect_true(treated$Proptosis_mm[treated$time == d] < untreated$Proptosis_mm[untreated$time == d])
})