# Top-level model assembly: wires PK + target engagement + IGF-1 biomarker +
# mechanism core + clinical translation into one deSolve system per drug.
# See model-spec.md for every equation and parameter source.

library(deSolve)

`%||%` <- function(a, b) if (is.null(a)) b else a

#' Run a TED QSP simulation.
#'
#' @param drug_id one of names(ted_drugs), e.g. "teprotumumab", "none".
#' @param regimen_override optional list to override default_regimen fields
#'   (dose_mg_per_kg, n_doses, interval_days).
#' @param body_weight_kg patient body weight, used for mg/kg dosing.
#' @param sex "male" or "female" -- affects teprotumumab Vc (published
#'   covariate effect).
#' @param trab_baseline patient's baseline disease-drive severity (1.0 =
#'   population-typical active TED at diagnosis).
#' @param disease_duration_days time already elapsed since TED onset before
#'   this simulation's t = 0 (treatment start) -- natural-history burnout is
#'   applied retroactively so long-standing patients start less active.
#' @param tau_burnout natural-history decay time constant, days.
#' @param w_igf1r,w_tshr transactivation weights (Section 5).
#' @param proptosis_baseline_mm,Pmax,wM,wA,M50,A50,G50 clinical-translation
#'   params (Section 7).
#' @param F_halflife,G_halflife,effector_halflife mechanism-core turnover
#'   half-lives in days (Section 6). effector_halflife applies to both the
#'   myofibroblast and adipocyte pools. Exposed as arguments (rather than
#'   hardcoded) so `scripts/calibrate_mechanism_core.R` can fit them against
#'   digitized trial data.
#' @param sim_days total simulation length.
#' @param dt output resolution, days.
#' @return data.frame with time + all state and derived-readout columns.
ted_simulate <- function(drug_id,
                          regimen_override = list(),
                          body_weight_kg = 70,
                          sex = c("female", "male"),
                          trab_baseline = 1.0,
                          disease_duration_days = 190,
                          tau_burnout = 730,
                          w_igf1r = 0.85,
                          w_tshr = 0.85,
                          proptosis_baseline_mm = 22,
                          Pmax = 9.696,
                          wM = 0.5,
                          wA = 0.5,
                          M50 = 1,
                          A50 = 1,
                          G50 = 1,
                          F_halflife = 3,
                          G_halflife = 21,
                          effector_halflife = 20.275,
                          sim_days = 250,
                          dt = 1) {

  sex <- match.arg(sex)
  drug <- ted_drugs[[drug_id]]
  if (is.null(drug)) stop("Unknown drug_id: ", drug_id)

  regimen <- utils::modifyList(if (is.null(drug$default_regimen)) list() else drug$default_regimen,
                                regimen_override)

  kF_off <- koff_from_halflife(F_halflife)
  kG_off <- koff_from_halflife(G_halflife)
  kM_off <- koff_from_halflife(effector_halflife)
  kA_off <- koff_from_halflife(effector_halflife)

  # --- dosing schedule -------------------------------------------------
  # IV doses are modeled as instantaneous bolus events (added to Ac) rather
  # than a finite-rate infusion forcing-function. The real infusions are
  # ~60-90 minutes -- negligible next to every other timescale here (weeks
  # between doses, days-to-weeks half-lives) -- and representing them as a
  # huge instantaneous rate (dose / ~0.06 days) over a very narrow window
  # was found to be numerically unreliable: the adaptive ODE solver
  # sometimes under-resolved that narrow high-rate pulse, silently losing
  # most of a dose depending on unrelated prior solver call history (see
  # git history / model-spec.md if this regresses again -- confirmed via
  # repeated identical calls giving different late-time trajectories).
  # Bolus events sidestep the issue entirely since there's no narrow
  # high-rate window for the solver to step over.
  has_pk <- !is.na(drug$pk_type)
  events_data <- NULL

  if (has_pk && drug$pk_type == "iv_2cmt_mm") {
    dose_per_kg <- rep_len(regimen$dose_mg_per_kg, regimen$n_doses)
    dose_mg <- dose_per_kg * body_weight_kg
    starts <- (seq_len(regimen$n_doses) - 1) * regimen$interval_days
    events_data <- data.frame(var = "Ac", time = starts, value = dose_mg, method = "add")
  } else if (has_pk && drug$pk_type == "oral_1cmt") {
    events_data <- build_oral_events(regimen$dose_mg, regimen$n_doses, regimen$interval_days)
  }
  breakpoints <- if (!is.null(events_data)) events_data$time else numeric(0)

  # --- initial conditions ------------------------------------------------
  drive0 <- disease_drive(disease_duration_days, trab_baseline, tau_burnout)
  y0 <- c(Fib = drive0, GAG = drive0, Myo = drive0, Fat = drive0, IGF1 = 1)
  if (has_pk && drug$pk_type == "iv_2cmt_mm") y0 <- c(y0, Ac = 0, Ap = 0)
  if (has_pk && drug$pk_type == "oral_1cmt")  y0 <- c(y0, Agut = 0, Ac = 0)

  vc <- if (!is.null(drug$pk_sex_vc)) drug$pk_sex_vc[[sex]] else if (has_pk) drug$pk$Vc else NA

  # --- ODE right-hand side -------------------------------------------------
  rhs <- function(t, y, parms) {
    with(as.list(y), {

      occ_igf1r <- 0; occ_tshr <- 0; trab_clear_eff <- 0
      dAc <- 0; dAp <- 0; dAgut <- 0

      if (has_pk && drug$pk_type == "iv_2cmt_mm") {
        pk <- pk_2cmt_mm_deriv(Ac, Ap, drug$pk$CL, drug$pk$Q, vc, drug$pk$Vp,
                                drug$pk$Vmax, drug$pk$Km)
        dAc <- pk$dAc; dAp <- pk$dAp
        occ <- occupancy(pk$C, drug$pk$Km)
        if (drug$mechanism == "igf1r") occ_igf1r <- occ
        if (drug$mechanism == "tshr")  occ_tshr  <- occ
        if (drug$mechanism == "trab_clearance") {
          trab_clear_eff <- (drug$trab_clearance_emax %||% 0.6) * occ
        }
      } else if (has_pk && drug$pk_type == "oral_1cmt") {
        pk <- pk_oral_1cmt_deriv(Agut, Ac, drug$pk$ka, drug$pk$CL, drug$pk$Vc)
        dAgut <- pk$dAgut; dAc <- pk$dAc
        occ_igf1r <- occupancy(pk$C, drug$pk$Km)
      }

      dIGF1 <- igf1_deriv(IGF1, occ_igf1r, kin = 1, kout = 1, Imax = 0.71)

      drive <- disease_drive(t + disease_duration_days, trab_baseline, tau_burnout)
      eff_drive <- effective_drive(drive, occ_igf1r, occ_tshr, trab_clear_eff, w_igf1r, w_tshr)

      core <- mechanism_core_deriv(Fib, GAG, Myo, Fat, eff_drive, kF_off, kG_off, kM_off, kA_off)

      derivs <- c(core$dFib, core$dGAG, core$dMyo, core$dFat, dIGF1)
      names(derivs) <- c("Fib", "GAG", "Myo", "Fat", "IGF1")
      if (has_pk && drug$pk_type == "iv_2cmt_mm") derivs <- c(derivs, Ac = dAc, Ap = dAp)
      if (has_pk && drug$pk_type == "oral_1cmt")  derivs <- c(derivs, Agut = dAgut, Ac = dAc)

      list(derivs)
    })
  }

  out_times <- sort(unique(c(seq(0, sim_days, by = dt), breakpoints[breakpoints <= sim_days])))

  events_arg <- if (!is.null(events_data)) list(data = events_data) else NULL

  sol <- as.data.frame(ode(y = y0, times = out_times, func = rhs, parms = NULL,
                            events = events_arg, method = "lsoda"))

  # --- derived readouts -----------------------------------------------------
  sol$CAS <- cas_score(sol$GAG, G50)
  sol$Proptosis_mm <- proptosis_mm(sol$Myo, sol$Fat, proptosis_baseline_mm, Pmax, wM, wA, M50, A50)
  sol$Diplopia <- diplopia_grade(sol$Myo)
  if ("Ac" %in% names(sol) && has_pk) {
    sol$Conc_ugmL <- sol$Ac / vc
  } else {
    sol$Conc_ugmL <- 0
  }

  attr(sol, "drug") <- drug
  attr(sol, "regimen") <- regimen
  sol
}