DoxToxOmics — Visual Explainer

The problem

We detect heart damage only after it is irreversible.

Doxorubicin is a backbone chemotherapy whose dose ceiling is set by cardiotoxicity. Troponin and ejection fraction are lagging indicators — the 2022 European Society of Cardiology (ESC) and 2024 Heart Failure Society of America (HFSA) cardio-oncology guidance call for earlier, mechanism-linked markers.

Clinical stakes

Doxorubicin-induced heart failure climbs from 3–5% at 400 mg/m² to 18–48% at 700 mg/m² (max lifetime dose 400–550 mg/m²)
Childhood survivors carry a 5–0× lifetime cardiovascular mortality risk
Troponin rises and left ventricular ejection fraction (LVEF) falls only after irreversible injury

Two gaps DoxToxOmics addresses

Biomarkers — no early, mechanism-linked plasma readout exists
Class comparison — no matched iPSC-CM anthracycline-vs-TKI transcriptomic comparison has been published
Why iPSC-CMs — human genetics, scale, and the matched DToxS (Drug Toxicity Signature) compendium (54 drugs, 6 donors)

The detection gap — injury accrues before any clinical signal

Schematic of the clinical problem: troponin and ejection fraction cross detection thresholds only once damage is irreversible — motivating an earlier, mechanism-linked readout (2022 ESC guidelines).

Approach

One hypothesis, three Specific Aims.

Hypothesis. Doxorubicin induces a dose-dependent transcriptomic cascade in hESC-CMs — from transcriptional suppression and iron/metabolic stress at low dose to extracellular-matrix (ECM) remodeling and ion-channel dysregulation at therapeutic dose — that is mechanistically distinct from TKI cardiotoxicity and predictive of electrophysiological dysfunction.

Aim 1 — Transcriptomics

Differential expression at 1 & 2.5 µM
Fold-change dose metrics
Singular value decomposition (SVD; Hansen benchmark)

Aim 2 — Mechanisms + electrophysiology (EP)

Pre-ranked Reactome gene-set enrichment analysis (GSEA)
Mito-gene → electrophysiology linkage
11-scenario metabolic / reactive oxygen species (ROS) simulation

Aim 3 — Validation + translation

Anthracycline-vs-TKI (DToxS Arm 2)
O'Hara-Rudy action-potential (AP) modeling
Plasma protein quantitative trait locus (pQTL) · signatureSearch reversal

Aim 1 · transcriptomics

A broad, dose-graded response.

DEG counts — False discovery rate (FDR) < 0.05, |log₂FC| > 1. ~7,810 differentially expressed genes (DEGs) at 1 µM and ~8,929 at 2.5 µM; 10,403 union, 6,336 shared (edgeR/DESeq2 NB-GLM, effective n = 2).

A broad response (~7,810 / 8,929 DEGs of 18,312 expressed) — the anthracycline footprint, on a method-correct count
Sub-proportional dose response (slope 0.823, R² 0.607): the core program is largely engaged by 1 µM
First principal component (PC1) = 94.1% variance is the shared dose-response (drug) signal — loadings r = 0.864 with 1 µM log₂FC — so GSEA runs on the uncorrected ranking (removing PC1 on a 2-dose matrix inverts signs)

0

dose-proportionality slope

0

expressed genes tested

Aim 2 · pathways

Transcription suppressed; translation & OXPHOS relatively enriched — never reversing with dose.

Top GSEA NES — Transcription (RNA Pol II, SUMOylation; NES ≤ −2.3) is the dominant suppressed axis; translation & OXPHOS transcripts are relatively enriched (NES up to +3.7).

Pathway direction — Across all 1,816 Reactome pathways, **zero** reverse direction between doses — the response intensifies, it does not reorganize.

0 / 0

significant pathways (1 µM / 2.5 µM)

0

pathways consistently up (75 down, 0 flippers)

0

direction-flippers

Ferroptosis / iron-handling module

A targeted 13-gene panel resolves an axis that pathway-level GSEA under-weights: 11 of 13 genes are differentially expressed at 2.5 µM
The NRF2 antioxidant arm is suppressed (NFE2L2 −1.0, SLC7A11, HMOX1 −2.3) while iron import/storage rises (TFRC +1.8, FTH1/FTL +1.8–2.0); iron-handling activity climbs dose-dependently (ρ = +0.89). ABCB8 is model-sensitive.
A coherent iron-overload / impaired-defense signature — reported as a hypothesis-generating readout, not demonstrated ferroptotic death

Every pathway intensifies with dose — none reverses sign

A line per representative Reactome program from 1 µM to 2.5 µM. Suppression and activation both deepen with dose; across all 1,816 pathways, none change sign — the response intensifies, it does not reorganize.

Aim 2 · electrophysiology

From mitochondrial loss to calcium overload.

IK1 ↓ alone → +58 ms (+17%); combined ATP + ROS nets +34 ms (+10%) on a 339.6 ms baseline
Calcium mishandling confirmed — non-recovering scenarios fail to restore the ~50 nM diastolic baseline (systolic peak ~320–380 nM is a normal ~7–8× of diastole; the pathology is sustained diastolic Ca²⁺, not a rise versus control)
Run on the GMU Hopper cluster (SLURM 7200345), supporting the May-1 calcium hypothesis

Calcium overload scenarios — Peak intracellular [Ca²⁺]ᵢ per scenario; open markers recover, filled markers do not.

APD prolongation vs literature — The predicted +10–17% action-potential duration (APD90) prolongation sits inside the published bracket (Wang & Korth +32%, Ducroq +18%, Knoll ~17%, Iachetta −3% negative control). Predictions are first-order — direction and relative magnitude.

What calcium overload looks like, beat to beat

diastolic ≈ 50 nM · systolic peak ~300–380 nM (normal ~7–8× ratio)doxorubicin: elevated, non-relaxing diastolic

Schematic transients — the systolic peak (~300–380 nM) sits a normal ~7–8× above the ~50 nM diastolic baseline; under doxorubicin the abnormality is the failure to relax back to that baseline before the next beat, leaving a sustained, elevated diastolic [Ca²⁺]ᵢ — the calcium-mishandling signature seen in every non-recovering simulation (not a 7–8× rise versus control).

Calcium overload and APD prolongation, coupled on one beat

baselinedoxorubicin Vmdoxorubicin [Ca²⁺]ᵢ

Time-aligned membrane voltage (top) and intracellular calcium (bottom) for one beat. Under doxorubicin the intracellular Ca²⁺ fails to relax to its ~50 nM diastolic baseline (the systolic peak itself is a normal ~7–8× of diastole); the sustained, elevated diastolic Ca²⁺ holds the plateau (+17% APD₉₀) and, via the Na⁺/Ca²⁺ exchanger (NCX), drives the secondary depolarization seen as an early afterdepolarization (EAD) — the pro-arrhythmic endpoint. Schematic; direction and relative magnitude only.

The 14 single-channel action-potential simulations

O'Hara–Rudy CiPA model · each conductance scaled by 2^log₂FC at 2.5 µM · 1000-beat steady state · cycling through all 14 →

CACNA1CL-type Ca²⁺ · I_CaL1 / 14

Phase-2 plateau · EC coupling

▼ log₂FC -1.40FDR < 0.05shortens plateau

KCNH2hERG · I_Kr2 / 14

Phase-3 repolarization

▼ log₂FC -1.11FDR < 0.05prolongs APD₉₀

SCN5ANa_v1.5 · I_Na3 / 14

Phase-0 upstroke

▼ log₂FC -2.11FDR < 0.05upstroke collapse

KCNQ1KvLQT1 · I_Ks4 / 14

Phase-3 repolarization

▲ log₂FC +2.09FDR < 0.05shortens APD₉₀

KCND3Kv4.3 · I_to5 / 14

Phase-1 early repolarization

▼ log₂FC -1.97FDR < 0.05raised early plateau

KCNE1minK · I_Ks β6 / 14

I_Ks auxiliary subunit

▼ log₂FC -0.89n.s.subtle in isolation

KCNJ2Kir2.1 · I_K17 / 14

Resting-potential maintenance

▼ log₂FC -0.93FDR < 0.05+17% APD₉₀ (I_K1 50%)

KCNIP2KChIP2 · I_to β8 / 14

I_to auxiliary subunit

▲ log₂FC +0.72n.s.subtle in isolation

ATP1A1Na⁺/K⁺ ATPase α19 / 14

Ion homeostasis

▲ log₂FC +2.01FDR < 0.05mild shortening

ATP2B1PMCA110 / 14

Ca²⁺ extrusion

▼ log₂FC -2.25FDR < 0.05Ca²⁺ handling

PLNPhospholamban11 / 14

SERCA2a regulation

▼ log₂FC -0.66n.s.subtle in isolation

RYR2Ryanodine receptor 212 / 14

SR Ca²⁺ release

▼ log₂FC -1.63FDR < 0.05Ca²⁺ handling

SLC8A1NCX1 · I_NaCa13 / 14

Na⁺/Ca²⁺ exchange

▼ log₂FC -2.16FDR < 0.05Ca²⁺ handling

NALCNNALCN14 / 14

Non-selective leak current

▲ log₂FC +0.16n.s.subtle in isolation

Single-channel perturbations are subtle in isolation — most shift APD₉₀ only slightly; SCN5A (I_Na ↓) collapses the upstroke, and KCNJ2 (I_K1 50%) gives the +17% APD₉₀. The pro-arrhythmic phenotype emerges only when all 14 conductances change together (the combined scenarios above). Schematic AP shapes — direction and relative magnitude.

Aim 3 · class specificity

Anthracyclines are mechanistically distinct from TKIs.

Anthracycline vs TKI — Per-drug correlation with the DoxToxOmics signature across 54 DToxS drugs on a matched iPSC-CM platform.

4 / 4 anthracyclines rank top-4 of 54 (label-shuffling null p < 10⁻⁴)
Kruskal–Wallis H = 11.18, p = 3.7×10⁻³; Dunn's post-hoc (Bonferroni/Holm); Mann–Whitney p = 9.5×10⁻⁴
Robust to epirubicin removal (p = 2.1×10⁻³) — a class property, not one outlier
The novel contribution — first matched-platform comparison (Paper 1)

Anthracyclines lead: EPI / DOX / DAU rank 1 / 2 / 3 of 54 drugs

◤ EPI / DOX / DAU top 3

rank 1 — strongest correlationrank 54

Each bar is one of 54 real DToxS drugs (GSE217421), sorted by Spearman correlation with the corrected doxorubicin signature. Anthracyclines (rust) lead — EPI/DOX/DAU rank 1/2/3 (IDA 6th); anthra > TKI Mann–Whitney p = 1.3×10⁻⁴; DToxS's own doxorubicin ranks #2 (positive control).

Aim 3 · biomarkers

What could be read out in plasma.

Under a Hansen-style per-source background (panel ∩ differentially-expressed universe), the HPA secretome is the strongly-enriched, plasma-accessible route (p = 6.3×10⁻⁴³)
Heart-failure genome-wide association study (HF-GWAS): 35/47 genes, p = 0.014 (modestly enriched); Human Phenotype Ontology heart-failure set (HPO-HF) 11/20, n.s. (p = 0.65); broad plasma panel not enriched (p ≈ 1.0)
Data-driven lead BAG3 (log₂FC −2.13); CAMK1D · MEF2C · TBX5 · TBX3

Biomarker nomination is anchored on the secretome — secreted, plasma-accessible, genetically-instrumented DEGs (227 with cis-pQTL instruments) — rather than the broad SomaScan panel.

pQTL enrichment — −log₁₀(hypergeometric p) per source; the broad plasma panel is not enriched, the HF-anchored panels are.

The secretome is the plasma-accessible route; HF-GWAS is modestly enriched

HPA secretome6.00

HF-GWAS (35/47)1.85

HPO-HF (n.s.)0.19

Broad plasma panel0.00

−log₁₀(hypergeometric p) per source under a Hansen-style per-source background (panel ∩ DE universe); dashed line = p 0.05. HPA secretome p = 6.3×10⁻⁴³; HF-GWAS p = 0.014 (35/47); HPO-HF n.s. (p = 0.65). Data-driven lead BAG3 (log₂FC −2.13). Robust to the gene-list cut: on the top-1000 DEGs (500↑/500↓) with BH-FDR across panels, secretome q = 3×10⁻¹⁷, HF-GWAS q = 0.020, HPO-HF n.s.

Aim 3 · compound reversal

Candidate signature reversers — reported with caveats.

signatureSearch LINCS reversers — signatureSearch/LINCS connectivity (corrected signature, 8,140 compounds). (A) Top reversers by normalized connectivity score; (B) mimicker positive control.

Kinase/signaling inhibitors (EGFR/erlotinib analog, QL-X-138/mTOR, QL-XII-47/BTK, GDC-0879/RAF, triciribine/AKT) and protein-synthesis inhibitors (emetine, cephaeline, homoharringtonine) as transcriptomic-signature reversers
Reported with caveats — transcriptomic-similarity hypotheses for functional testing, not cardioprotectants; protein-synthesis inhibitors score as reversers only because the corrected signature carries translation transcripts upward
Mimicker positive control passes cleanly: the anthracyclines/anthracenediones are recovered at the top of all 8,140 compounds (epirubicin #7, mitoxantrone #9, daunorubicin #12, doxorubicin #55, idarubicin #66)

Methods

How we did the math.

Each method shown with its defining equation and a small diagram of what the equation is computing — LaTeX rendered via MathJax.

Differential expression · OLS on VST counts · Wald · BH FDR

Per-gene linear model on variance-stabilized counts; per-dose treatment contrasts give log₂FC; Wald p with Benjamini–Hochberg FDR control.

In this work18,312 expressed genes across 12 hESC-CM libraries (effective n = 2/dose). ~7,810 DEGs @ 1 µM and ~8,929 @ 2.5 µM (6,336 shared, 10,403 union; edgeR/DESeq2 NB-GLM). Sub-proportional dose: slope 0.823 (R² 0.607).

$$y_{ij}\;=\;\beta_0 + \beta_1\,\mathbb{1}_{d=1\mu\mathrm{M}} + \beta_2\,\mathbb{1}_{d=2.5\mu\mathrm{M}} + \varepsilon_{ij}$$

GSEA (pre-ranked) · Reactome · running enrichment score

Walk the gene list ranked by signed significance; the running sum rewards hits inside the gene set and penalizes misses outside; maximum deviation = ES; normalize by gene-set size for NES; permutation-based p.

In this work1,786 Reactome pathways tested per dose; 315 / 385 significant (1 µM / 2.5 µM). Suppressed: RNA Pol II transcription NES −2.33, SUMOylation −2.45; translation/OXPHOS transcripts relatively enriched (NES +3.7). Zero pathways flip direction between doses.

$$\mathrm{NES}\;=\;\frac{\mathrm{ES}}{\overline{\bigl|\mathrm{ES}_{\text{null}}\bigr|}}\,,\qquad p_{\text{perm}}\;=\;\Pr\!\bigl[\mathrm{ES}_{\text{null}}\ge\mathrm{ES}\bigr]$$

SVD / PCA · Hansen pipeline · drug vs culture axis

Decompose the dose × gene matrix into orthogonal singular components; semantic interpretation by loadings (drug axis vs culture axis vs noise).

In this workOn the NB-GLM matrix, PC1 = 94.1% of variance, loadings r = 0.864 with 1 µM log₂FC — the drug signal. (Hansen's 252-condition LINCS P0 has PC1 = culture axis; here PC1 IS the drug, so removing it on a 2-dose matrix would invert signs — so pathway analysis uses the pre-SVD matrix.)

$$X \;=\; U\,\Sigma\,V^{\!\top}\,,\qquad \mathrm{PC}_k\,\text{score} = U_{:,k}\,\sigma_k$$

Pearson correlation · DToxS Arm 2 · matched-platform signature

Per-drug correlation of each DToxS drug's transcriptomic signature against the doxorubicin signature on a matched iPSC-CM platform; rank by r.

In this work54 drugs on the REAL public DToxS panel (GSE217421), corrected signature. EPI/DOX/DAU rank 1/2/3 (Spearman ρ 0.71/0.38/0.32); anthra > TKI Mann–Whitney p = 1.3×10⁻⁴; DToxS's own DOX ranks #2 (positive control).

$$r \;=\; \frac{\sum_i (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_i (x_i - \bar{x})^2\,\sum_i (y_i - \bar{y})^2}}$$

Kruskal–Wallis + Dunn's · Non-parametric ANOVA on ranks

Non-parametric ANOVA on ranks tests whether the three drug classes have different median signature correlations; Dunn's post-hoc with Bonferroni / Holm correction.

In this workK-W H = 11.18, p = 3.7 × 10⁻³ across the three classes; pairwise Mann–Whitney p = 9.5 × 10⁻⁴. Robust to epirubicin removal (p = 2.1 × 10⁻³) — a class property, not one outlier.

$$H \;=\; \frac{12}{N(N+1)} \sum_{i=1}^{k} \frac{R_i^{2}}{n_i} \,-\, 3(N+1)$$

Permutation / label-shuffle null · Distribution-free p-value

Shuffle the drug-class labels many times to build a null distribution for the test statistic; tail probability gives a distribution-free p-value.

In this workB = 10,000 shuffles. Under random labels, probability of all 4 anthracyclines landing in the top 4 of 54 is p < 10⁻⁴. Independent of distributional assumptions.

$$\hat{p} \;=\; \frac{1 \;+\; \#\bigl\{T^{(b)} \,\ge\, T_{\text{obs}}\bigr\}}{1 \;+\; B}$$

Hypergeometric enrichment · Panel ∩ DEGs vs background

Test overlap of the DEG set with a curated panel (HF-GWAS, HPO, plasma) under a Hansen-style per-source background — the panel intersected with the gene-detection universe of that source.

In this workAgainst the DEG union: HF-GWAS 35/47, p = 0.014; HPO-HF n.s. (p = 0.65). Broad plasma panel not enriched (p ≈ 1.0). HPA secretome strongly enriched (p = 6.3×10⁻⁴³) — the plasma-accessible route; lead BAG3 (log₂FC −2.13).

$$P(X = k)\;=\;\frac{\binom{K}{k}\,\binom{N-K}{\,n-k\,}}{\binom{N}{n}}$$

O'Hara–Rudy AP simulation · Myokit CiPA v1.0 · 14 conductances · 1000 beats

Hodgkin–Huxley-style ODEs over 14 ionic conductance parameters; each scaled by 2^log₂FC at 2.5 µM; integrated to steady state via Myokit (O'Hara–Rudy CiPA v1.0).

In this work1000-beat steady state; baseline APD₉₀ = 339.6 ms. I_K1 50% → +58 ms (+17%); combined ATP + ROS → +34 ms (+10%) — inside the literature bracket (Wang & Korth +32%, Ducroq +18%, Knoll ~17%; Iachetta negative control). Systolic peak ~320–380 nM over a ~49 nM diastolic baseline (a normal ~7–8× systole:diastole ratio, NOT a rise versus control); the abnormality is the failure to restore diastolic Ca²⁺ in every non-recovering scenario. (Hopper SLURM 7200345.)

$$C_m\,\frac{dV}{dt} \;=\; -\!\!\sum_{X}\! I_X\,,\qquad g_X^{\text{dox}} \;=\; g_X^{\text{base}}\cdot 2^{\log_2 \mathrm{FC}_X}$$

Synthesis

One coherent mechanistic cascade.

Every arrow is backed by a specific analysis above — and the chain is echoed by three independent 2025–26 models (Schaffert adverse-outcome-pathway [AOP] network, Dean p53-MDM2-p21 loop, Kolitsida Mfn2/NCLX).

Molecular trigger

DNA intercalation

topoisomerase-IIβ (TOP2B) poisoning · ROS generation

01

Organelle stress

Transcriptional suppression

RNA Pol II · SUMOylation ↓ (NES −2.3); translation/OXPHOS relatively ↑

02

Organelle stress

Iron dysregulation + TF shutdown

TFRC/FTH1/FTL ↑ · GATA4/E2F1/BRCA1 ↓

03

Cellular remodeling

ECM remodeling + UPR

unfolded protein response (UPR); stress programs ↑

04

Cellular remodeling

Ion-channel remodeling

16 channels → 14 model conductances

05

Electrical phenotype

APD prolongation

+58 ms (+17%) when I_K1 ↓

06

Electrical phenotype

Calcium overload

fails to relax to ~50 nM diastolic (sustained elevation)

07

Electrical phenotype

EAD emergence

early afterdepolarization → pro-arrhythmic

08

The phenotypic endpoint — a prolonged action potential with an early afterdepolarization

baseline APdoxorubicin: +APD₉₀ & EAD

Schematic — illustrative of direction and relative magnitude, not a quantitative O'Hara–Rudy trace. Baseline APD₉₀ = 339.6 ms; I_K1↓ alone → +58 ms (+17%); the EAD emerges under combined ATP-deficit + ROS when intracellular Ca²⁺ fails to relax to its ~50 nM diastolic baseline (sustained elevation).

Logistics

Path to Spring 2027.

Jun 16 2026Proposal defense · 9:00 AM ← scheduled

Fall 2026BINF 998 · candidacy processed · Paper 1 submit (Dec 2026)

Spring 2027BINF 999 (6 cr) · revise Ch 1–5 · figures

Mar 31 2027Format review + graduation application hard deadline

Apr 2027Dissertation defense

May 22 2027Degree conferral

1

Schedule the proposal defense

Target June 2026 — availability within the committee is the only dependency.

2

Raise methodological concerns now

Better addressed here than at the defense; the work is current and defensible.

3

Confirm the Spring 2027 path

BINF 998 (Fall) + BINF 999 (Spring); candidacy processed in fall; format review Mar 31; defense Apr–May 2027.