# TED QSP Model — Development Documentation
This document is the project's narrative record: where every piece of data
came from, how it was obtained, how the code was designed and built
(including the bugs found along the way and how they were diagnosed), and a
complete bibliography. For the equations and parameter values themselves,
see `docs/model-spec.md` — this document explains *how those numbers were
obtained* and *why the code looks the way it does*; it does not repeat the
equations verbatim.
- **Data provenance** → Section 1
- **Code development history** → Section 2
- **Deployment** → Section 3
- **Full references** → Section 4
---
## 1. Data collection and provenance
### 1.1 What "data collection" meant for this project
No patient-level dataset was used or is available to this project. Every
number in the model comes from **published, publicly accessible sources** —
peer-reviewed journal articles and, in one case, the primary trial
publication's own text-layer figure labels. Three distinct collection
methods were used, and it matters which one produced which number:
1. **Direct quote from article text/tables** (highest confidence) — a
sentence or table cell in a paper states a parameter outright (e.g. "CL =
0.334 L/day"). Obtained via `WebSearch` and `WebFetch` against the
article's PMC/journal page.
2. **Text-layer extraction from a figure** (used once, for the calibration
data) — some PDFs embed the *actual numeric data labels* printed on a
chart as selectable text in the PDF's text layer, not just as pixels.
When present, this is exact, reported data — not an estimate — even
though it comes from a "figure." This is how the OPTIC trial's
visit-by-visit proptosis curve was obtained (Section 1.3).
3. **Structural/mechanistic assumption** (not data at all) — where no
quantitative source exists, a value was chosen for physical plausibility
and flagged `[ASSUMED]` in `model-spec.md`, never presented as if it were
measured.
Nothing in this project was obtained by visually reading pixel positions off
a rendered chart image (traditional plot-digitization, e.g. WebPlotDigitizer-
style). Method 2 above achieves the same *result* (a numeric time-series
extracted from a figure) but by a more precise mechanism, and is called out
explicitly in `data/digitized/README.md` and this document so that
precision claim can be checked rather than assumed.
### 1.2 Literature search process
Research was done via `WebSearch` across roughly a dozen targeted queries
covering: TED pathophysiology, teprotumumab mechanism and PK, the OPTIC
trial results, the CAS/proptosis/diplopia outcome instruments, the TED drug
development pipeline (veligrotug, VRDN-003, linsitinib, batoclimab, K1-70),
and R/Shiny/QSP engineering precedent (mrgsolve tutorials). Promising hits
were then opened with `WebFetch` (for HTML/PMC pages, which return cleanly)
or downloaded and read as PDF (for the NEJM OPTIC paper, whose HTML page
would not parse — the PDF was fetched via two independent mirror URLs,
cross-checked against each other, and read with the `Read` tool's native PDF
support).
An explicit search for prior QSP/mechanistic ODE models of TED was also run
(`"mathematical model Graves orbitopathy fibroblast adipogenesis cytokine
simulation ODE"`) and came back empty of any actual mechanistic model — only
in-vitro and mouse-model papers. **This model's mechanism core (Section 6 of
model-spec.md) is therefore an original structure proposed for this
project, not adapted from a published QSP model** — worth knowing when
judging how much to trust its untested parts (see Section 8 of
model-spec.md).
### 1.3 The OPTIC trial — the project's primary quantitative anchor
**Douglas RS, Kahaly GJ, Patel A, et al. "Teprotumumab for the Treatment of
Active Thyroid Eye Disease." N Engl J Med. 2020;382(4):341-352.**
DOI: 10.1056/NEJMoa1910434. NCT03298867.
This is a Phase 3, randomized, double-masked, placebo-controlled trial (N=41
teprotumumab / N=42 placebo) and is the single most important source in this
project — it supplied:
- **Baseline patient characteristics** (Table 1): mean age, sex ratio,
disease duration (6.2-6.4 months), baseline proptosis (22.6-23.2mm),
baseline CAS (5.1-5.3) → `data/digitized/optic_trial_baseline.csv`.
- **The visit-by-visit proptosis-change curve** (Figure 2B): least-squares
mean change from baseline in proptosis at weeks 0/6/12/18/24, both arms →
`data/digitized/optic_trial_proptosis_change.csv`. This is the **only**
continuous, repeated-measures endpoint in the paper that maps directly
onto a model output (`Proptosis_mm`), which is why it — and *only* it —
was used as the mechanism-core's curve-fitting target
(`scripts/calibrate_mechanism_core.R`).
- **Binary responder-rate curves** (Figure 2A/C/D/E): % proptosis responders,
% CAS-0-or-1, % overall responders, % diplopia responders, at each visit →
`data/digitized/optic_trial_responder_rates.csv`. These are population
threshold-crossing *rates*, not continuous mean scores, so they were
**not** used to fit `cas_score()` or `diplopia_grade()` — doing that
properly would require a virtual-population/threshold layer this model
doesn't have yet (flagged as future work in model-spec.md Section 8).
- **GO-QOL quality-of-life changes** (Figure 2F) →
`data/digitized/optic_trial_goqol.csv`, collected but **not used anywhere
in the current model** (no GO-QOL output exists) and flagged low-
confidence on visit-ordering because the PDF text extraction didn't
cleanly preserve which label belonged to which visit for that specific
panel.
**How Figure 2's exact numbers were obtained**: the PDF's text layer
includes the numeric data labels Douglas et al. printed directly on the
chart (e.g. the digits "−3.32", "78", "22" appear as selectable text in the
PDF, positioned where they visually sit on the figure). This was confirmed
by reading the full PDF with the `Read` tool (which returns PDFs as parsed
text+layout, not images) and cross-referencing every extracted number
against the same values stated in the paper's own body text and Table/
Figure 1B summary table (e.g. the abstract states "the mean change in
proptosis (−2.82 mm vs. −0.54 mm)" for the AUC-style through-week-24 summary,
which is consistent with, but numerically distinct from, Figure 2B's
week-24-specific values of −3.32mm/−0.53mm — both were kept, clearly
labeled which is which, rather than conflating them).
### 1.4 Teprotumumab population pharmacokinetics
**Xu et al. "Pharmacokinetics and Exposure-Response Relationship of
Teprotumumab, an Insulin-Like Growth Factor-1 Receptor-Blocking Antibody, in
Thyroid Eye Disease."** PMC8332554.
Population PK model built from 10 healthy subjects (single 1500mg IV dose)
plus 176 TED trial patients. Supplied every number in model-spec.md Section
3.1: two-compartment structure, CL, Q, Vc (with the male/female covariate
split), Vp, Vmax, Km, inter-individual variability (%CV) for each, and the
reported geometric-mean terminal half-life (19.9 days). Obtained via
`WebFetch` directly against the PMC page (parsed cleanly as HTML/text, no
PDF needed).
### 1.5 Veligrotug (VRDN-001) preclinical pharmacology
**"Preclinical pharmacology, pharmacokinetics, and pharmacodynamics of
veligrotug, a full antagonist antibody to the IGF-1 receptor in development
for thyroid eye disease."** PMC12629334.
Cynomolgus monkey non-compartmental PK data (single-dose and 13-week
repeat-dose) at 2/10/50 mg/kg — no human population PK model was found in
the literature search. Used only to establish the *order of magnitude*
(clearance, volume, half-life, and the ~2.5-fold IGF-1 rise) for veligrotug's
`confidence = "low"` registry entry, which reuses teprotumumab's human PK
structure as a scaled stand-in rather than a fitted human model. Explicitly
flagged in model-spec.md Section 3.2 as illustrative-only until real human
data is substituted.
### 1.6 TED pathophysiology (mechanism core's biological grounding)
Three sources anchor the disease-biology narrative in model-spec.md Section
1 and justify the model's structural choices (why TSHR/IGF-1R blockade share
a high transactivation weight, why three effector pools):
- **Krieger CC, et al. "Mechanisms in Thyroid Eye Disease: The TSH Receptor
Interacts Directly with the IGF-1 Receptor."** Endocrinology.
2025;166(2):bqaf009.
— direct ectodomain TSHR-IGF-1R interaction, the mechanistic basis for the
model's transactivation weights (`w_IGF1R`, `w_TSHR`).
- **"TSHR-IGF-IR complex drives orbital fibroblast misbehavior in thyroid
eye disease."** PubMed 39082947.
- **"A Comprehensive Review of Thyroid Eye Disease Pathogenesis: From Immune
Dysregulations to Novel Diagnostic and Therapeutic Approaches."**
Int J Mol Sci (MDPI). 2024;25(21):11628.
— Th1→Th2/profibrotic
cytokine shift over the disease course, the biphasic active→inactive
natural history, and the three-fate fibroblast differentiation (myo-
fibroblast / adipocyte / GAG-producer) that the model's F→{M,A,G}
structure is built around.
### 1.7 Clinical outcome instruments (CAS, proptosis, diplopia)
- **Mourits MP, et al.** "Clinical activity score as a guide in the
management of patients with Graves' ophthalmopathy." Clin Endocrinol
(Oxf). 1997;47:9-14. (The original CAS instrument — cited within the
OPTIC paper itself, not independently fetched, but its 7-item structure
is what `cas_score()`'s 0-7 range models.)
- **"Inter-observer Variability of Clinical Activity Score."** PubMed
37028695. — CAS's known
reliability limitations, noted in model-spec.md Section 1 as a reason the
model's CAS output shouldn't be over-interpreted at high precision even
once calibrated.
- **"Multidimensional predictive model for assessing clinical activity in
thyroid eye disease."** PMC12279714.
### 1.8 TED drug development pipeline (drug registry entries beyond teprotumumab)
Used to populate `R/drug_registry.R`'s mechanism/route/regimen metadata for
veligrotug, VRDN-003, linsitinib, batoclimab, and K1-70 — **contextual
sourcing only**, not quantitative parameterization (those entries are
`confidence = "low"` or `"stub"` and say so on-screen in the app):
- Viridian Therapeutics TED Programs page.
— veligrotug/
VRDN-003 mechanism and development stage.
- "An Efficacy, Safety, and Tolerability Study of Veligrotug (VRDN-001)…
(THRIVE-2)." ClinicalTrials.gov NCT06021054.
- "Emerging therapies in the medical management of thyroid eye disease."
PMC11182121. —
linsitinib, batoclimab, K1-70 mechanism summaries.
### 1.9 QSP/Shiny engineering precedent (informed the software-stack decision, not the disease model)
- Lu et al. "gPKPDviz: A flexible R shiny tool for pharmacokinetic/
pharmacodynamic simulations using mrgsolve." CPT Pharmacometrics Syst
Pharmacol. 2024. — confirmed Shiny+mrgsolve is the standard pattern for
this kind of tool; informed the decision to use `deSolve` instead (see
Section 2.2) once Rtools was found to be unavailable on this machine.
- PMX Solutions, "Creating a simple pharmacometric Shiny application with
mrgsolve in R" (Parts 1-2), 2018.
---
## 2. Code development history
This section documents *how the code came to be the way it is* — the
sequence of decisions, and, importantly, the bugs that were found during
development and exactly how they were diagnosed. This project hit two
genuine, non-obvious bugs during build-out; both are recorded here in full
because the debugging process itself is useful context for anyone extending
this code later.
### 2.1 Architecture
The model is built as five composable layers (full detail in model-spec.md
Section 2), so that adding a new drug never requires touching the disease-
mechanism code:
```
PK (per drug) → target engagement (occupancy) → disease drive (natural
history × drug effect) → mechanism core (fibroblast → effector pools,
shared across all drugs) → clinical translation (CAS/proptosis/diplopia)
```
`R/drug_registry.R` is the single place a new drug's PK parameters,
mechanism (`igf1r` / `tshr` / `trab_clearance`), and dosing regimen are
declared; `R/simulate.R`'s `ted_simulate()` dispatches on that declaration
and assembles one `deSolve::ode()` system per call.
### 2.2 Why `deSolve` instead of `mrgsolve`/`rxode2`
QSP models in R are conventionally built with `mrgsolve` (compiled C++ ODE
engine) — confirmed as the standard approach via the engineering-precedent
search in Section 1.9. This project uses plain `deSolve` instead, because
Rtools (the C++ toolchain `mrgsolve` needs to compile a model) was checked
and found not installed on the development machine
(`C:\Program Files\R\R-4.5.1` exists; no `rtools45` directory). `deSolve`
ships as a CRAN binary and needs no compiler, at the cost of slower
virtual-population simulation — acceptable for this app's current
single-patient interactive use case. The ODE right-hand-side in
`R/simulate.R` is written so a later port to `mrgsolve` (if population-scale
simulation becomes necessary) is mechanical, not a rewrite.
### 2.3 Calibration methodology and the identifiability problem it uncovered
`scripts/calibrate_mechanism_core.R` fits the mechanism core's turnover
rate(s) and `Pmax` against the OPTIC proptosis curve (Section 1.3) using
`optim()`. The full story, in order:
1. **First attempt** — fit all four of `F_halflife`, `effector_halflife`,
`Pmax`, `M50` jointly with unconstrained Nelder-Mead. It converged to a
numerically *better*-fitting but physically nonsensical point:
`F_halflife → 4.8e-7 days`, `Pmax → 8.8e5`, `M50 → 1.7e5`. Diagnosis: with
`Myo`/`Fat` ranging roughly 0-1 in this model (normalized to 1 at a
sustained baseline drive of 1), once `M50 >> 1` the saturating term
`M/(M+M50)` degenerates to the linear approximation `M/M50` — so `Pmax`
and `M50` are only identified through their *ratio*, never separately,
and the optimizer is free to send both to infinity along that ratio.
2. **Second attempt** — added box constraints (`L-BFGS-B`) to keep the
search in a plausible region. The fit then just sat at whichever bound
was given (`F_halflife` at its lower bound, `Pmax`/`M50` at their upper
bounds) — a clear sign the *bounds*, not the data, were now doing the
identifying, which is not a real calibration.
3. **Fix** — removed the degenerate parameters from the fit entirely rather
than fighting the optimizer: `M50 = A50 = 1` fixed (removes the ridge
structurally), `F_halflife = 3` days fixed (the 4 visit-level data points
can't distinguish "3 days" from "instantaneous" for an unobserved
upstream signal). Refitting just `effector_halflife` and `Pmax` (2
parameters, 8 data points) landed both cleanly inside their bounds —
`effector_halflife ≈ 20.3 days`, `Pmax ≈ 9.7` — a genuinely identified
result. Final RMSE ≈ 0.24mm across both arms and all 4 visits.
This sequence — and why the "better-fitting" first attempt was rejected —
is recorded in comments at the top of `scripts/calibrate_mechanism_core.R`
and in model-spec.md Sections 6-8, not just here, so the reasoning survives
independent of any one file.
### 2.4 Bug #1: an infusion-rate forcing function that silently lost doses
**Symptom**: a test asserting the calibrated model reproduces the OPTIC
proptosis curve within 0.3mm failed. Manual inspection showed the model,
run standalone, predicted the *disease relapsing* between weeks 18 and 24 —
physiologically wrong for a drug given continuously through week 21.
**Diagnosis process**: traced to the drug concentration trajectory itself,
which showed troughs collapsing to <0.05 µg/mL between later doses — far
below what the drug's own published Km (1.5 µg/mL) would suggest for a
saturating, accumulating regimen. The root cause turned out to be **how IV
dosing was represented**: each dose was modeled as a finite-rate infusion
(dose ÷ ~0.0625 days ≈ tens of thousands of mg/day, sustained for ~90
minutes) rather than an instantaneous bolus. Direct experiment confirmed
the bug was in the ODE solver's handling of that narrow, extreme-rate pulse:
**running the exact same simulation twice, back to back, in the same R
session, with byte-identical parameters, produced different late-time
trajectories** depending on what the solver had been asked to integrate
immediately beforehand in that session — i.e. the result depended on
solver call history, not on the model's actual inputs, which is a
correctness bug, not merely an accuracy limitation.
**Fix**: switched all IV dosing from a rate-based forcing function to
instantaneous bolus events (`deSolve`'s `events` mechanism, already used for
oral dosing) — physically justified since a 60-90 minute infusion is
negligible next to every other timescale in the model (weeks between doses,
days-to-weeks half-lives). After the fix, troughs build up realistically
across the 8-dose regimen (26 → 69 → 88 → 96 → 100 → 101.6 → 102 µg/mL
immediately before doses 2-8), consistent with veligrotug's preclinical
trough range (Section 1.5), and repeated identical simulations became
fully reproducible. See `R/simulate.R`'s dosing-schedule comment and
model-spec.md Section 3.1 ("Resolved numerical issue") for the full
technical account.
### 2.5 Bug #2: `rbind()` across drugs with different PK state columns
**Symptom**: after deployment, the live app showed "An error has occurred"
on every load. The Shiny Server log
(`/var/log/shiny-server/*.log` inside the container) showed: `Error in
rbind: numbers of columns of arguments do not match`, at `app.R`'s `sim()`
reactive.
**Diagnosis**: `ted_simulate()`'s output columns depend on the drug's PK
type — a drug with `pk_type = "iv_2cmt_mm"` (e.g. teprotumumab) has `Ac`/`Ap`
state columns in its output; the `"none"` (no-treatment) pseudo-drug has no
PK compartment at all and so has neither. The app's "overlay no-treatment"
feature did `rbind(treated, untreated)` on these two data frames directly,
which fails whenever the column *sets* differ (base R's `rbind.data.frame`
matches by name but still requires the same set of names on both sides).
**Fix**: subset both data frames to just the columns the app actually plots
(`time`, `scenario`, `CAS`, `Proptosis_mm`, `Diplopia`, `Conc_ugmL`, `IGF1`)
before combining, in `app/app.R`'s `sim()` reactive.
**A second, deployment-specific wrinkle surfaced while verifying this
fix**: re-requesting the app page right after editing `app.R` on disk still
showed the old error. Shiny Server keeps a persistent R worker process per
app across requests within its idle timeout, and that worker had already
`source()`'d the buggy code into memory before the fix landed — editing the
file on disk doesn't reach an already-running worker. The fix only took
effect after explicitly killing that worker process (`ps aux | grep R`
inside the container, `kill -9` the PID) so Shiny Server spawned a fresh one
on the next request. This is now documented as a standing operational note
(see the NAS deployment memory / `docs/DEVELOPMENT.md` Section 3) since it
will recur on every future code change deployed this way.
### 2.6 Testing strategy
`tests/testthat/` favors checking *properties that are actually true*
rather than asserting brittle exact numbers pulled from a single run:
- **PK** (`test-pk.R`): the pure-linear sub-model (Vmax=0) is checked
against the *closed-form* 2-compartment beta-phase half-life derived from
the published CL/Q/Vc/Vp (not against a remembered simulation output) —
this is what caught Bug #1 as a real discrepancy rather than plausible
noise. A second test explicitly asserts the *full* model (with TMDD)
decays faster than the linear-only one once concentration nears Km — the
expected, documented, real pharmacology (Section 3.1 of model-spec.md),
distinguished on purpose from the Bug #1 numerical artifact.
- **Calibration** (`test-calibration.R`): guards that the shipped default
parameters still reproduce the OPTIC proptosis curve within a stated RMSE
bound, so a future change to the mechanism core that silently breaks the
fit is caught by CI/test runs rather than discovered live.
---
## 3. Deployment
Full operational detail (container name, ports, SSH access, the "must kill
the stale R worker after any code edit" gotcha) lives in the auto-memory
file `project-ted-shiny-nas.md` (cross-referenced from `project-deploy-
policy.md`) rather than duplicated here, since that's operational state
that can change independently of the model code. In short: `Y:\apps` is
itself a network mount of the NAS path bind-mounted into a
`rocker/shiny-verse` Docker container (`rocker-shiny-verse-server`, port
3838), reachable at `https://shiny.kaseidon.com/ted/app/` via the NAS's Web
Station host routing — no separate deploy/copy step exists for this project,
unlike the `npm run deploy` static-site projects elsewhere in this repo.
---
## 4. Full references
### Primary quantitative sources (directly parameterized or calibrated the model)
1. Douglas RS, Kahaly GJ, Patel A, Sile S, Thompson EHZ, Perdok R, et al.
Teprotumumab for the Treatment of Active Thyroid Eye Disease. *N Engl J
Med*. 2020;382(4):341-352. doi:10.1056/NEJMoa1910434. NCT03298867.
2. Xu C, et al. Pharmacokinetics and Exposure-Response Relationship of
Teprotumumab, an Insulin-Like Growth Factor-1 Receptor-Blocking Antibody,
in Thyroid Eye Disease. PMC8332554.
3. Preclinical pharmacology, pharmacokinetics, and pharmacodynamics of
veligrotug, a full antagonist antibody to the IGF-1 receptor in
development for thyroid eye disease. PMC12629334.
### Pathophysiology / mechanistic grounding
4. Krieger CC, Boutin A, Jang D, et al. Mechanisms in Thyroid Eye Disease:
The TSH Receptor Interacts Directly with the IGF-1 Receptor.
*Endocrinology*. 2025;166(2):bqaf009.
5. TSHR-IGF-IR complex drives orbital fibroblast misbehavior in thyroid eye
disease. PubMed 39082947.
6. A Comprehensive Review of Thyroid Eye Disease Pathogenesis: From Immune
Dysregulations to Novel Diagnostic and Therapeutic Approaches. *Int J Mol
Sci*. 2024;25(21):11628.
7. Transcriptomic Profiling of Thyroid Eye Disease Orbital Fat Demonstrates
Differences in Adipogenicity and IGF-1R Pathway. PMC11665563.
8. PDGFRα+DPP4+ Fibroblasts-Macrophage Crosstalk Induces Orbital Fibrosis in
Treatment-Resistant Thyroid Eye Disease via the GAS6-AXL Pathway.
PMC12677676.
### Clinical outcome instruments
9. Mourits MP, Prummel MF, Wiersinga WM, Koornneef L. Clinical activity
score as a guide in the management of patients with Graves'
ophthalmopathy. *Clin Endocrinol (Oxf)*. 1997;47:9-14. (cited within
Douglas et al. 2020, ref. 1 above)
10. Inter-observer Variability of Clinical Activity Score: Assessments in
Patients With Thyroid Eye Disease. PubMed 37028695.
11. Multidimensional predictive model for assessing clinical activity in
thyroid eye disease. PMC12279714.
12. Teprotumumab for thyroid eye disease: early response is not required for
benefit. PMC9232498.
### TED drug pipeline (context for `R/drug_registry.R` stub entries)
13. Viridian Therapeutics. Thyroid Eye Disease Programs.
14. An Efficacy, Safety, and Tolerability Study of Veligrotug (VRDN-001), in
Participants With Chronic Thyroid Eye Disease (TED) (THRIVE-2).
ClinicalTrials.gov NCT06021054.
15. Emerging therapies in the medical management of thyroid eye disease.
PMC11182121.
### QSP/software engineering precedent
16. Lu J, et al. gPKPDviz: A flexible R shiny tool for pharmacokinetic/
pharmacodynamic simulations using mrgsolve. *CPT Pharmacometrics Syst
Pharmacol*. 2024.
17. PMX Solutions. Creating a simple pharmacometric Shiny application with
mrgsolve in R (Parts 1-2). 2018.
### Background reading (informed general understanding, not directly cited for a number)
18. Teprotumumab in the management of thyroid eye disease: mechanistic
insights and adverse reactions — a comprehensive review. *Front
Endocrinol*. 2025. doi:10.3389/fendo.2025.1480195.
19. Teprotumumab for Thyroid Eye Disease: Mechanism, Clinical Efficacy, and
Current Challenges. PMC12286191.
20. A Deep Learning Framework for Predicting Teprotumumab Treatment Response
in Thyroid Eye Disease. *ScienceDirect*. 2026.
doi:10.1016/j.xops.2026.100369 (S2666914526000369).