# TED QSP Model — Specification Quantitative systems pharmacology model of Thyroid Eye Disease (TED / Graves' orbitopathy). This document is the source of truth for every equation, parameter, and its provenance. Code should not contain numbers that aren't traceable to a line in this file. For the narrative version — how each source was found and collected, why the code is architected this way, and the bugs found while building it, plus a full bibliography — see `docs/DEVELOPMENT.md`. Confidence tags used throughout: - **[LIT]** — taken directly from a published human PK/PD estimate. - **[CALIB]** — chosen so the model reproduces a published summary observation (e.g. fold-change, response rate), not itself independently measured. - **[ASSUMED]** — a structural or numerical modeling choice with no direct literature support; a placeholder to refine as better data is found or as the app is calibrated against digitized trial curves. ## 1. Disease biology encoded TSH-receptor (TSHR) and IGF-1 receptor (IGF-1R) autoantibodies co-engage on orbital fibroblasts; the two receptors transactivate each other (direct ectodomain interaction), so blocking either receptor substantially blunts the combined signal, not just its own fraction (Krieger 2024/2025, PMC11665563; academic.oup.com/endo/166/2/bqaf009). Activated fibroblasts differentiate into myofibroblasts (muscle enlargement), adipocytes (orbital fat expansion), and GAG/hyaluronan producers (edema/chemosis), amplified by a Th1 cytokine loop (TNF-α, IFN-γ, IL-6) early and Th2/profibrotic signaling later (MDPI 25:11628). Disease is classically biphasic: an active inflammatory phase followed by spontaneous "burnout" into an inactive/fibrotic phase over roughly 1-3 years **[ASSUMED timescale, qualitative from clinical literature]**. Clinical readouts: CAS (7-item EUGOGO activity score, 0-7, known to have poor inter-rater reliability — PMC12279714), proptosis in mm (Hertel exophthalmometer, ≥2mm change considered clinically meaningful), diplopia grade (0-3, Gorman scale convention). ## 2. Architecture ``` [Dosing regimen] -> [Drug PK] -> [Target engagement] -+ v [Baseline TRAb, natural-history decay] -> [Effective disease drive] | v [Fibroblast activation, F] (cytokine-amplified) / | \ v v v [Myofibroblast,M] [Adipocyte,A] [GAG/hyaluronan,G] \ | / v v v [Proptosis] [Diplopia] [CAS] ``` Every drug plugs into this core via one of three mechanism front-ends (`R/target_engagement.R`): - `igf1r_occupancy` — mAb blocks IGF-1R (teprotumumab, veligrotug, VRDN-003). - `tshr_occupancy` — mAb blocks TSHR directly (K1-70-like). - `trab_clearance` — FcRn antagonism accelerates autoantibody (IgG) clearance, reducing the upstream drive rather than blocking a receptor (batoclimab-like). The mechanism core (`R/mechanism_core.R`) and clinical translation (`R/clinical_translation.R`) are shared and never change per-drug — only the front-end swaps. This is the "swappable multi-mechanism" property the app is built around. ## 3. PK module ### 3.1 Teprotumumab — **[LIT]**, source: Xu et al., population PK, PMC8332554 (10 healthy subjects + 176 TED patients). Two-compartment, linear distribution + parallel Michaelis-Menten (saturable, target-mediated) elimination from the central compartment. Amounts in mg, concentration in µg/mL (= mg/L, consistent with Vc/Vp in L). ``` dAc/dt = - CL*(Ac/Vc) - Q*(Ac/Vc) + Q*(Ap/Vp) - Vmax*(Ac/Vc) / (Km + Ac/Vc) dAp/dt = Q*(Ac/Vc) - Q*(Ap/Vp) C(t) = Ac/Vc ``` Dosing enters as an instantaneous bolus addition to `Ac` (via deSolve `events`), not as a finite-rate infusion term in the ODE — see the note in `R/simulate.R`. The real infusions are ~60-90 minutes, negligible next to every other timescale here (weeks between doses, days-to-weeks half-life); modeling them as a rate forcing-function over that narrow a window turned out to be numerically unreliable (see "Resolved numerical issue" below). | Parameter | Value | IIV (%CV, lognormal) | |---|---|---| | CL | 0.334 L/day | 43.0% | | Q | 0.859 L/day | 47.8% | | Vc | 3.94 L (male) / 3.25 L (female) | 27.7% | | Vp | 4.21 L | 25.4% | | Vmax | 3.02 mg/day | — | | Km | 1.5 µg/mL | — | Resulting geometric-mean terminal half-life ≈ 19.9 days (26.2% CV) at clinical doses. TMDD saturates ≥3 mg/kg; clinical regimen (10 mg/kg x1, then 20 mg/kg x7, q3w IV infusion, label regimen) is well past saturation at Cmax, so PK is near-linear immediately post-infusion. **Resolved numerical issue** — an earlier version of this model represented each IV dose as a finite-rate infusion (dose / ~0.06 days, i.e. a rate of several thousand mg/day sustained for ~90 minutes) rather than a bolus. That narrow high-rate pulse turned out to be unreliably resolved by deSolve's adaptive-step `lsoda` integrator: identical parameters sometimes produced wildly different late-time trajectories depending on unrelated prior solver call history within the same R session (confirmed by direct reproduction — repeated calls with byte-identical inputs gave different results only depending on what the ODE solver had been asked to integrate immediately beforehand). Simulated troughs before later q3w doses were falling below 1 µg/mL (<20% occupancy) as an artifact of this — silently losing most of a dose's effect. Switching every IV dose to an instantaneous bolus event (adding directly to `Ac`, no rate term in the derivative) resolved it: troughs now build up smoothly across the 8-dose regimen (26 → 69 → 88 → 96 → 100 → 101.6 → 102 µg/mL immediately before doses 2-8), consistent with the veligrotug preclinical trough range (50-350 µg/mL, PMC12629334) for a similar mAb/dose, and the apparent post-last-dose half-life (~14 days measured on a linear-model fit, days 150-170) is much closer to the published ~19.9 days than the ~5-13 days seen under the buggy rate-based approach. Residual gap between 14 and 19.9 days is expected and fine — it's the real, published TMDD behavior described in Section 3.1 above (elimination genuinely speeds up once concentration approaches Km), not a numerical artifact. No significant exposure-response relationship was found in the source paper between teprotumumab trough concentration and CAS/proptosis/diplopia response — noted here so the app does not overclaim a dose-response for this drug's clinical effect; response was high even in the lowest exposure quartile (66.7% vs 14.9% placebo). ### 3.2 Veligrotug (VRDN-001) — **[ASSUMED]**, low confidence. Only preclinical (cynomolgus monkey) NCA data is published (PMC12629334): CL ~7.6-13.6 mL/day/kg, Vss ~59-107 mL/kg, t½ 3.2-7.2 days (dose-dependent, consistent with TMDD at low dose). No human population-PK model has been located. The registry uses the *teprotumumab* 2-cpt+MM structure with allometrically scaled CL/V (scaled by a placeholder 70 kg reference weight) as a stand-in, flagged `confidence = "low"` in `R/drug_registry.R`. **Treat veligrotug outputs as illustrative only until a real human PK model is substituted.** ### 3.3 VRDN-003, linsitinib, batoclimab, K1-70 — **[ASSUMED]**, stubs. Module structures are wired (1-cpt oral absorption for linsitinib; 2-cpt IV for the mAbs) but parameters are placeholders pending literature review. Flagged `confidence = "stub"` — the app should visibly warn when one of these is selected. ## 4. Target engagement Receptor occupancy modeled as a simple saturable binding function reusing the PK model's own Km as the effective affinity constant (consistent with the quasi-steady-state TMDD approximation, where Km approximates target binding affinity) — **[CALIB]**, avoids introducing an independent Kd with no data to fit it to: ``` Occupancy(t) = C(t) / (C(t) + Km) ``` ### IGF-1 biomarker (validation target for the IGF-1R front-end) Indirect-response (turnover) model: receptor blockade reduces IGF-1 receptor-mediated consumption, raising free IGF-1: ``` d(IGF1)/dt = kin - kout * (1 - Imax * Occupancy(t)) * IGF1 IGF1(0) = kin / kout (baseline) ``` `Imax = 0.71` **[CALIB]** chosen so full occupancy yields a ~3.4-fold rise in steady-state IGF-1, bracketing the published 100-350% (2-4.5-fold) increase range (PMC8332554; veligrotug preclinical ~2.5-fold). ## 5. Effective disease drive ``` Drive(t) = TRAb_baseline * exp(-t / tau_burnout) [natural history] EffDrive(t) = Drive(t) * (1 - w_IGF1R * Occ_IGF1R(t)) * (1 - w_TSHR * Occ_TSHR(t)) * (1 - TRAbClearanceEffect(t)) ``` `tau_burnout` **[ASSUMED]** default 730 days (~2 years), reflecting the classic biphasic active→inactive natural history. `w_IGF1R = w_TSHR = 0.85` **[ASSUMED]** — set high (rather than 0.5/0.5) to reflect the published TSHR/IGF-1R transactivation biology: blocking either receptor removes most of the synergistic signal, not just an independent additive half. This is a structural modeling choice, not a fitted value — flagged for sensitivity analysis. ## 6. Mechanism core All four states are turnover (production/loss) models driven by `EffDrive(t)` or by F. | State | Meaning | Turnover t½ | Feeds | |---|---|---|---| | F | Fibroblast activation signal | 3 days **[ASSUMED]** | M, A, G production rate | | G | GAG / hyaluronan pool | 21 days **[ASSUMED]** | CAS | | M | Myofibroblast pool | 20.3 days **[CALIB]** | Proptosis, diplopia | | A | Adipocyte pool | 20.3 days **[CALIB]** (shared with M) | Proptosis | ``` dF/dt = kF_on * EffDrive(t) - kF_off * F dG/dt = kG_on * F - kG_off * G dM/dt = kM_on * F - kM_off * M dA/dt = kA_on * F - kA_off * A ``` `k*_off = ln(2)/t½`; `k*_on` set so baseline (untreated, steady disease drive = 1) states equilibrate at a normalized value of 1.0, i.e. `k*_on = k*_off` at drive=1 steady state. `effector_halflife` (M/A) was fit (Section 8) against the OPTIC trial's digitized proptosis curve and landed well inside its search bounds (5-150 days), away from either edge — a genuinely identified parameter, not an artifact of wherever the bounds happened to be. `F_halflife` is **fixed** at 3 days rather than fit: when left free, it repeatedly pushed to whatever the lower search bound was (down to fractions of a day) without improving the fit meaningfully — the 4 visits (weeks 6/12/18/24) are too coarsely spaced to distinguish "3 days" from "instantaneous" for an upstream signal that's itself unobserved, so it isn't its own identifiable timescale from this data. `G_halflife` (CAS) remains **[ASSUMED]** — no continuous mean-CAS-by-visit data was available to fit it at all (see Section 8). `A` shares `effector_halflife` with `M`: with `M50 = A50`, this makes `Myo(t) ≡ Fat(t)` identically, which is why `wM`/`wA` are not separately identifiable from proptosis data alone (see Section 7). ## 7. Clinical translation ``` CAS(t) = 7 * G / (G + G50), G50 = 1 [ASSUMED, bounds CAS to 0-7] Proptosis(t) = Proptosis_baseline + Pmax * (wM * M/(M+M50) + wA * A/(A+A50)), M50=A50=1 [ASSUMED], Pmax=9.696 [CALIB] wM=wA=0.5 [ASSUMED, unidentifiable from proptosis alone -- see Section 6] Diplopia(t) = round(3 * M / (M + D50)), D50 = 1.5 [ASSUMED, 0-3 ordinal] ``` `M50`/`A50` are **fixed at 1**, not fit: an earlier attempt to fit them jointly with `Pmax` hit a hard identifiability wall. `Myo`/`Fat` range roughly 0-1 in this model (normalized to 1 at a sustained baseline drive of 1), so once `M50 >> 1`, `M/(M+M50) ≈ M/M50` for the entire trajectory — `Pmax` and `M50` are then only identified jointly via their ratio, never separately (the fit drove both toward diverging, physically meaningless values — `F_halflife → 5e-7 days`, `Pmax → 8.8e5`, `M50 → 1.7e5` — before bounds were added, and pushed to whatever bound was given afterward, which is itself a sign the *bound* was doing the identifying, not the data). Fixing `M50 = A50 = 1` removes that degenerate direction; `Pmax` is then cleanly identified (Section 8). ## 8. Calibration status (honest accounting) Done / literature-grounded: - Teprotumumab PK (Section 3.1) — validated against published half-life; dosing switched from a rate-based infusion to bolus events after finding the former was numerically unreliable (see "Resolved numerical issue," Section 3.1) — troughs now build up realistically across the 8-dose regimen instead of collapsing between doses. - IGF-1 biomarker fold-change (Section 4) — calibrated to published range. - **Mechanism core effector turnover + Pmax — fit to the OPTIC trial's digitized proptosis-change-from-baseline curve** (`data/digitized/optic_trial_proptosis_change.csv`, Douglas et al. NEJM 2020, Figure 2B) via `scripts/calibrate_mechanism_core.R` (bounded `optim`/L-BFGS-B, joint fit across both treatment arms and all 4 visits). Only 2 parameters were fit (`effector_halflife`, `Pmax`) against 8 points — `F_halflife` and `M50`/`A50` are fixed, not fit, for the identifiability reasons in Sections 6-7. RMSE = 0.242 mm across the 8 points, with both fitted parameters landing well inside their bounds (not pinned to an edge). Predicted vs. observed (mm): placebo wk6 −0.054 vs −0.38, wk12 −0.169 vs −0.64, wk18 −0.299 vs −0.59, wk24 −0.431 vs −0.53; teprotumumab wk6 −1.95 vs −2.00, wk12 −2.90 vs −2.70, wk18 −3.20 vs −3.26, wk24 −3.31 vs −3.32. The teprotumumab arm fits well throughout; the placebo/ natural-history arm captures the right order of magnitude and general upward trend but not the real data's small non-monotonic dip (wk18 −0.59 → wk24 −0.53, a slight *worsening* then improvement) — a smooth ODE can't reproduce that wiggle, plausibly just sampling noise at N=42. `disease_duration_days` default updated to 190 (trial mean ≈ 6.3 months) to match the calibration's patient assumption. Not yet done (flagged for follow-up, not hidden): - **CAS is still not curve-fit** — OPTIC only reports a binary "% CAS 0 or 1" responder rate by visit (`data/digitized/optic_trial_responder_rates.csv`), not a continuous mean-CAS trajectory, so there's nothing for `cas_score()` to be fit against directly without adding a virtual-population/threshold layer (Phase 6 stretch goal). `G_halflife` (21 days) is still a qualitative guess — worth revisiting once a continuous CAS source is found (candidates: the phase 2 Smith et al. 2017 NEJM trial, or the OPTIC supplementary appendix Fig. S2, neither pulled in this pass). - Diplopia (Gorman grade) is not fit to anything — OPTIC's diplopia responder curve (22/68/68/68% teprotumumab vs 4/18/21/29% placebo, `optic_trial_responder_rates.csv`) exists but is binary/ordinal, same issue as CAS. - Veligrotug, VRDN-003, linsitinib, batoclimab, K1-70 PK/PD parameters are still placeholders (Section 3.2-3.3) — none of this calibration touched them, since OPTIC is teprotumumab-only. - `w_IGF1R`, `w_TSHR`, `tau_burnout` are still structural assumptions, not fitted — good candidates for a sensitivity-analysis tab (Phase 6). The placebo arm's near-flat 24-week trajectory is at least consistent with `tau_burnout = 730` days not being wildly wrong (a much shorter burnout time constant would show more spontaneous improvement in the placebo arm than the ~0.5 mm actually observed). ## 9. Software `deSolve` (pure R, no C++ compiler required — Rtools is not installed on this machine) rather than `mrgsolve`/`rxode2`. The ODE right-hand-side is written as one function combining PK + target engagement + IGF-1 + mechanism core + translation, in `R/simulate.R`, so migrating to `mrgsolve` later is a mechanical port of the same equations, not a rewrite.