In-depth mechanistic analysis including high-throughput RNA sequencing in the prediction of functional and structural cardiotoxicants using hiPSC cardiomyocytes

ABSTRACT

Background

Cardiotoxicity remains one of the most reported adverse drug reactions that lead to drug attrition during pre-clinical and clinical drug development. Drug-induced cardiotoxicity may develop as a functional change in cardiac electrophysiology (acute alteration of the mechanical function of the myocardium) and/or as a structural change, resulting in loss of viability and morphological damage to cardiac tissue.

Research design and methods

Non-clinical models with better predictive value need to be established to improve cardiac safety pharmacology. To this end, high-throughput RNA sequencing (ScreenSeq) was combined with high-content imaging (HCI) and Ca²⁺ transience (CaT) to analyze compound-treated human-induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs).

Results

Analysis of hiPSC-CMs treated with 33 cardiotoxicants and 9 non-cardiotoxicants of mixed therapeutic indications facilitated compound clustering by mechanism of action, scoring of pathway activities related to cardiomyocyte contractility, mitochondrial integrity, metabolic state, diverse stress responses and the prediction of cardiotoxicity risk. The combination of ScreenSeq, HCI and CaT provided a high cardiotoxicity prediction performance with 89% specificity, 91% sensitivity and 90% accuracy.

Conclusions

Overall, this study introduces mechanism-driven risk assessment approach combining structural, functional and molecular high-throughput methods for pre-clinical risk assessment of novel compounds.

KEYWORDS:

• Cardiotoxicity
• HT-transcriptomics
• high-content imaging
• CaT
• hiPSC-CMs
• in vitro
• predictive

1. Introduction

Cardiotoxicity remains a predominant cause of pre-clinical and clinical drug failure [1,2]. Consequently, pre-clinical strategies with better predictive power need to be developed to guarantee the efficacy and safety of novel pharmaceutical drugs on cardiovascular functions, thus reducing later-stage attrition. Cutting-edge techniques and non-clinical models that closely represent the in vivo situation are needed to de-risk cardiotoxicity in early drug discovery and development.

The high incidence of drug-induced cardiotoxicity led to the adoption of the International Conference of Harmonization (ICH) S7B guideline in 2005, which called for the preclinical evaluation of new drug entities on cardiac electrophysiology using in vitro electrophysiology studies (typically human Ether-a-go-go Related-Gene (hERG) channel test) and in vivo QT assays in animal models [3]. This highly sensitive approach has been successful in reducing the percentage of proarrhythmic drug submissions to the U.S. Food and Drug Administration (FDA) [4]; however, it comes with low specificity and presents some limitations. Not all hERG channel blockers, for example, cause QT prolongation or Torsades de Pointes (TdP) (e.g. verapamil, due to its concomitant blockade of the depolarizing inward calcium current) [5] and, similarly, QT prolongation (often via blocking of the hERG channel) does not necessarily elicit proarrhythmia (e.g. ranolazine, due to the concomitant blockade of the depolarizing late inward sodium current) [6,7]. As a result, some valuable drug candidates may have been discarded due to false-positive risks or, if approved, their clinical use is limited by inappropriate warnings. Consequently, the Comprehensive in vitro Proarrhythmia Assay (CiPA) initiative was established to develop a new non-clinical safety paradigm, intended to move safety pharmacology toward an in silico and in vitro approach building on the emergence of new technologies such as stem-cell-derived cardiomyocytes [8].

Drug-induced cardiotoxicity can be functional in nature (defined as an acute alteration in the mechanical function of the myocardium), or it can also develop as a structural change due to morphological damage to cardiomyocytes and/or loss of viability [9]. In vitro cardiotoxicity screening strategies have predominantly focused on the detection of ECG abnormalities and QT interval prolongation, however assays that detect structural cardiotoxicity have emerged in the recent years, shedding light on the mechanisms leading to structural damage [10–12]. Likewise, high-throughput assays that collect transcriptional profiles and evaluate pathological pathway activities for drug-induced cardiotoxicity are currently limited in the scientific literature. Recent advances in the generation and large-scale production of human-induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) and human embryonic stem cell-derived cardiomyocytes (hESC-CMs) have allowed the development of novel high-throughput in vitro assays to detect cardiotoxicity. The utility of hESC-CMs to predict drug-induced structural cardiotoxicity was first shown in a high-content screening assay by live-cell fluorescent imaging of mitochondrial membrane potential (MMP), endoplasmic reticulum (ER) integrity, Ca²⁺ mobilization and membrane permeability combined with cellular adenosine triphosphate (ATP), which offered a sensitivity and specificity of 74% [13]. Subsequently, other studies have demonstrated the amenability of hiPSC-CMs to multi-parametric automated assays using both functional (beating activity) and structural (e.g. cell viability, reactive oxygen species generation, glutathione (GSH) depletion, lipid formation and troponin secretion) readouts for cardiotoxicity prediction [10,11,14–18]. hiPSC-CMs have been shown to hold great potential as an in vitro model for cardiotoxicity assays since they show fundamental electrophysiological and pharmacological characteristics, a cardiac-specific transcription profile, gene expression of key ion channels, excitation wave propagation and excitation-contraction coupling [19–23].

The aim of this publication was threefold: 1) to develop and validate high-throughput ribonucleic acid (RNA) sequencing (ScreenSeq) for the detection of alterations in molecular pathways in hiPSC-CMs as a new in vitro strategy to predict drug-induced cardiotoxicity 2) to demonstrate the feasibility of using hiPSC-CMs in structural assays such as high-content imaging (HCI) as well as functional assays, such as Ca²⁺ transient (CaT) analysis 3) to combine ScreenSeq, HCI and CaT analysis to provide a highly accurate cardiotoxicity prediction platform, while integrating complementary levels of compound response information and thereby providing a better mechanistic understanding of compounds. To this end, hiPSC-CMs were dosed with a set of 42 reference compounds typically cited in the literature for assay validation including 33 well-known functional and structural cardiotoxicants and 9 non-cardiotoxicants. Cardiotoxicity was then assessed using a high-content imaging (HCI), cellular ATP and Ca²⁺ transient assays, in addition to high-throughput RNA-sequencing (ScreenSeq).

2. Materials and methods

2.1. Materials

Human-induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) and cell culture media (CD-4) were provided by Evotec International GmbH. CD-4 media was prepared using Iscove’s modified Dulbecco’s Medium (IMDM) supplemented with 10% fetal bovine serum (FBS), 0.1% BioXtra human insulin (Sigma Cat. No. I9278) and 0.04% 1-thioglycerol (Sigma Cat. No. M6145). Experiments were conducted using two genetically distinct iPSC lines: EVOiPS0254 and EVOiPS0274. Two batches of the iPSC line EVOiPS0274 were also used for reproducibility studies. Both iPS lines were genetically modified to harbor a puromycin selection cassette driven by a cardiac-specific promotor using a proprietary protocol (puromycin concentration for selection was optimized as puromycin resistance varies between different batches). iPSCs were dissociated and aggregated into embryoid bodies for 1 day followed by differentiation into cardiomyocytes for a subsequent 17 days. After cardiomyocytes were differentiated, EBs were dissociated and the differentiated cells were seeded into T175 flasks for an additional 10 days for maturation and final selection of cardiomyocytes by puromycin treatment. The resulting cardiomyocyte preparations were then frozen and stored in liquid nitrogen until further use. Compounds were purchased from Sigma-Aldrich (Dorset, UK) or Cambridge Biosciences and were of the highest purity available.

2.2. Cell culture and compound treatment

hiPSC-CMs were seeded in 384-well plates (Corning® Cat. No.: 3764) pre-coated with a 0.1 µg/mL fibronectin solution (Sigma Cat. No. F1141) in phosphate buffered saline (PBS) containing Mg²⁺ and Ca²⁺. Fibronectin pre-coated plates were incubated at 37°C for 3 h (5% CO₂, 95% humidity), after which the coating was carefully removed using an automated liquid handling robot (Bravo Agilent). Subsequently, hiPSC-CMs were seeded at 10,000 cells/well and cultured for 10 days in CD-4 media. Media was initially refreshed 24 h post-seeding and then, every 48 h thereafter, while the cells were in culture. Compounds were selected across several therapeutic indications to cover a broad chemical space and, for assay validation, included a variety of mechanisms: 12 structural cardiotoxicants, 14 functional cardiotoxicants, 7 structural/functional cardiotoxicants and 9 non-cardiotoxicants (see Table 1). Top concentrations of compounds were based on 100 × maximum total human plasma concentrations (C_max) or solubility limit, while the lowest concentrations were lower than C_max or 50 × free concentration, calculated from the mean unbound fraction (Table 1). All compounds were dissolved in 100% dimethylsulfoxide (DMSO) (stock solutions). The stock solutions were further diluted in CD-4 media, and cells were dosed in triplicate at an 8-point dose response range using a half-log dilution series to achieve a final DMSO concentration of 0.5% (v/v). Vehicle controls were dosed with 0.5% (v/v) DMSO, while positive controls were dosed with sunitinib, a functional and structural cardiotoxicant.

Table 1. Summary of the 42-reference drug set. Table shows C_max values (µM), estimated free concentrations (µM), tested dose range (µM), approved/withdrawn year, clinical use and mechanism of action with primary targets.

Compound		Cmax (µM)	Estimated free concentrations (µM)	Tested dose range (µM)	First approved-withdrawn	Clinical use	Mechanism of action with primary targets
Non-cardiotoxicants	Acetaminophen	150	121	1.6–4000	1950-NA	Central Nervous System agent	Cyclooxygenase inhibitor
	Acetylsalicylic acid	6.68	3.77	0.4–1000	1899-NA	Anti-inflammatory agent	Cyclooxygenase inhibitor
	Acyclovir	2.87	2.43	0.4–1000	1982-NA	Antiviral agent	Thymidine kinase inhibitor
	Amoxicillin	10.6	8.89	0.12–300	1973-NA	Antibiotic	Penicillin-binding protein inhibitor
	Buspirone	0.01	0.001	0.004–10	1986-NA	Anxiolytic	Serotonin (5-HT1A) agonist
	Enalapril	0.83	0.37	0.04–100	1985-NA	Cardiovascular agent	Angiotensin-converting-enzyme (ACE) inhibitor
	Furosemide	4.45	0.074	0.2–500	1966-NA	Cardiovascular agent	Sodium-potassium-chloride cotransporter (NKCC2) inhibitor
	Sildenafil	0.24	0.009	0.016–40	1998-NA	Sexual dysfunction	Phosphodiesterase 5 (PDE5) inhibitor
	Tolbutamide	217	9.95	2.4–6000	1957-NA	Antidiabetic	Sulfonylurea receptor 1 inhibitor
Cardiotoxicants	Amitriptyline	14.5	0.83	0.024–60	1961-NA	Antidepressant	Serotonin-norepinephrine reuptake inhibitor
	Atenolol	2.9	2.78	0.1–250	1975-NA	Cardiovascular agent	β1-adrenergic receptor antagonist
	Bepridil	3.18	N/A	0.012–30	1990-NA	Cardiovascular agent	Ca2+ and Na+ channel blocker
	Cisapride	0.11	0.028	0.004–10	1980–2000	Gastrointestinal agent	Serotonin-4 (5-HT4) receptor agonist
	Digoxin	0.0035	0.003	0.00004–0.1	1954-NA	Cardiovascular agent	Na+/K+ ATPase inhibitor
	Diltiazem	0.3	0.058	0.02–50	1982-NA	Cardiovascular agent	L-type Ca2+ channel blocker
	Dobutamine	1.4	0.32	0.04–100	1978-NA	Cardiovascular agent	β1-adrenergic agonist
	Dopamine	0.12	0.12	0.008–20	1974-NA	Central Nervous System agent	Dopamine receptor agonist
	Epinephrine	0.002	0.002	0.0004–1	1901-NA	Anaphylaxis	Adrenergic and dopamine D4 receptor agonist
	Levosimendan	0.132	0.003	0.02–50	2000-NA	Cardiovascular agent	Cardiac troponin C binding agent, ATP-sensitive potassium channels opener
	Lidocaine	12.5	3.90	0.5–1250	1948-NA	Central Nervous System agent	Na+ channel blocker and epidermal growth factor receptor (EGFR) inhibitor
	Nifedipine	0.28	0.011	0.02–50	1981-NA	Cardiovascular agent	L-type Ca2+ channel blocker
	Propranolol	0.2	0.027	0.08–200	1967-NA	Cardiovascular agent	β-adrenergic receptor blocker
	Sotalol	12.3	9.71	0.5–1250	1992-NA	Cardiovascular agent	β-adrenergic receptor blocker, hERG blocker
	3’-azido-3’-deoxythymidine (AZT)	44.9	N/A	08–2000	1987-NA	Antiviral agent	Nucleoside reverse transcriptase inhibitor (NRTI)
	Amphotericin b	62.9	5.39	0.06–165	1958-NA	Antifungal	Binds with ergosterol
	Bortezomib	0.30	0.050	0.02–50	2003-NA	Antineoplastic agent	26S proteasome inhibitor
	Clozapine	0.25	0.013	0.02–50	1989-NA	Antipsychotic	Inhibitor of dopamine D2 and serotonin 2A receptors
	Cyclophosphamide	153	132.7	1.6–4000	1959-NA	Antineoplastic agent, immunosuppressant	DNA cross-linking
	Dasatinib	0.59	0.030	0.016–40	2006-NA	Antineoplastic agent	Inhibitor of kinases BCR-ABL, SCR, EPHA2, c-KIT and PDGFRβ
	5-Fluorouracil	4.6	2.95	0.2–500	1962-NA	Antineoplastic agent	DNA cross-linking, thymidylate synthase inhibitor
	Imatinib	4	0.23	0.03–75	2001-NA	Antineoplastic agent	Inhibitor of BCR-ABL tyrosine kinase, PDGF and SCF kinases
	Mitomycin c	7.1	N/A	0.4–1000	1974-NA	Antineoplastic agent	DNA cross-linking
	Mitoxantrone	3.3	0.83	0.004–10	1987-NA	Anti-neoplastic agent	DNA cross-linking, DNA topoisomerase II inhibitor
	Rofecoxib	0.03	0.004	0.01–25	1999–2004	Anti-inflammatory agent	Cyclooxigenase-2 inhibitor
	Rosiglitazone	1.44	0.003	0.08–200	1999–2010 (EMA)	Anti-diabetic	PPARG agonist
	Amiodarone	1.27	0.008	0.016–40	1985-NA	Cardiovascular agent	hERG, Na^+, and Ca^2+ channel blocker; noncompetitive α- and β-adrenergic inhibition
	Doxorubicin	11.7	2.97	0.02–50	1974–NA	Antineoplastic agent	DNA cross-linking, DNA topoisomerase II inhibitor
	Idarubicin	0.12	0.009	0.008–20	1990-NA	Antineoplastic agent	DNA cross-linking, DNA topoisomerase II inhibitor
	Isoproterenol	0.008	0.006	0.004–10	1948-NA	Cardiovascular agent	Nonselective β adrenoreceptor agonist and TAAR1 agonist
	Lapatinib	1.70	0.017	0.008–20	2007-NA	Antineoplastic agent	EGFR and HER2 receptor tyrosine kinases inhibitor
	Sunitinib	35.12	N/A	0.03–75	2006-NA	Antineoplastic agent	Multi-targeted receptor tyrosine kinase inhibitor (PDGFR, VEGFR, FLT, KIT, CSF-1 R)
	Verapamil	0.5	0.054	0.02–50	1978-NA	Cardiovascular agent	L-type Ca^2+, channel blocker

2.3. High-content imaging (HCI) and Ca²⁺ transient (CaT) measurements

Upon 10 days culture, hiPSC-CMs were incubated with the EarlyTox^TM Cardiotoxicity fluorescent dye (Molecular Devices). Following a 2-h incubation, hiPSC-CMs were dosed with compounds for 0 h (acute injection) or 24 h. Upon compound treatment, fast kinetic fluorescent reading was performed to detect individual Ca²⁺ transient peaks, which provided a multi-parametric transient profile (Ex 485 nm/Em 528 nm). The following parameters were measured and recorded over a 12-s period: amplitude, frequency, full peak width, full width at half maximum (FWHM), full rise time, rise time from 10%, full decay time, decay time to 10%, peak count, peak width at 10%, peak width at 50%, peak spacing. Readings were taken on a Cytation 3 Cell Imaging Multi-Mode Reader (BioTek), and raw data was analyzed using a multiscale-based peak detection (AMPD) method as described by Scholkmann et al. [24]. High-content imaging (HCI) was performed after the Ca²⁺ flux measurements using an ArrayScan® HCI reader (ThermoScientific) and a series of fluorescent dyes to quantify cell count, nuclear size, deoxyribonucleic acid (DNA) structure (Hoechst), calcium homeostasis (EarlyTox^TM), mitochondrial mass and mitochondrial membrane potential (Tetramethylrhodamine, ethyl ester, TMRE). Finally, cellular ATP was measured to assess gross cytotoxicity using the CellTiter-Glo Cell Viability Assay (Promega) as per manufacturer’s instructions. Minimum effective concentrations (MEC) and AC₅₀ values were calculated for each parameter and were used for the calculation of prediction metrics (sensitivity, specificity, and accuracy).

2.4. Whole genome high-throughput transcriptomics (ScreenSeq)

High-throughput transcriptomics was performed using Evotec’s (semi)-automated ScreenSeq platform. Cell lysis was performed in a 384-well plate format followed by well-specific labeling of samples. Purified mRNAs were labeled with a unique molecular identifier (UMI). After pooled cDNA synthesis, adapters containing Illumina’s unique dual indexes were attached for final library preparation and sequenced on a NovaSeq 6000 system. Transcriptomic profiling was performed in batches of 382 samples. Each batch contains 22 negative-control (i.e. DMSO-treated) samples and 15 blocks of compound-treated samples. Each block included eight compound concentrations, each in triplicate. Reads were mapped to the human genome (genome-build GRCh38.p13 – accession GCA_000001405.28) using Spliced Transcripts Alignment to a Reference (STAR) [25] (version 2.7.9a). Multi-mapped reads, reads with mapping to intergenic regions (with respect to Ensembl gene annotation information version 104) and reads with ambiguous mapping location (i.e. location associated with more than one gene) were discarded, and UMI-based deduplication was performed. Log-scaling of the resulting counts used CP10K (counts per 10 thousand) as scaling and 1 as pseudo-count unless indicated otherwise.

Differential expression analysis was performed separately for each sequencing batch, and – within batch – separately for each individual comparison (i.e. compound in each concentration vs DMSO) using DESeq2 [26] [DESeq2] (version 1.34) R package. To increase robustness by excluding individual outliers while keeping the procedure identical for all sequencing batches, 20 DMSO-treated samples were used as a control in these analyses for each batch, and the two with the lowest correlation with batch-specific averaged DMSO profiles were excluded. Here, Pearson correlation coefficients and log-scaled expression profiles were used, and a robust version of averaging was applied: for each gene, three lowest and three highest values were excluded, and the arithmetic mean was calculated for the remaining 16 values.

For triplicates of compound-treated samples, the following exclusion criteria were applied. Firstly, for each sequencing batch, the lower quartile of total counts for DMSO-treated samples was calculated, and a sample was excluded if it had total counts lower than one-fourth of this lower quartile. Secondly, for each sample in a triplicate, a Pearson correlation coefficient was calculated between the log-scaled expression profiles of this sample and each of the other two samples; additionally, a correlation coefficient between the expression profiles of this sample and an averaged profile of the other two samples was calculated. The maximum of these three numbers was used to quantify the consistency of a sample with its replicates. A sample was excluded if this quantification fell below 0.975.

2.5. Sample clustering

Principal Component Analysis (PCA) was performed using R package irlba (version 2.3.3), with centering and no scaling; the analysis was limited to 500 top features. For comparison of DMSO samples, log-scaled expression levels were used as feature values, and feature ranking was based on variance (i.e. 500 genes with the highest variability of expression level were selected). For studying compound effects, replicates were treated as one data point, and log-fold-changes identified within differential expression analysis were used as feature values. Feature selection in this case was based on the 10^th percentile of adjusted p-values associated with a gene in all performed comparisons (one comparison for each compound-concentration pair): genes with minimum (i.e. closest to zero) value of this percentile were selected.

2.6. Pathway enrichment

GO (gene ontology release 2020–09-10; GO term annotation for genes: UniProt release 2020_05) and WikiPathway (release 20,210,410) enrichment analyses used genes with adjusted p-values below 0.05 (and no fold-change threshold) as sets of differentially expressed genes; Fisher’s exact test (or, equivalently, hypergeometric test) was applied to identify enriched pathways and GO terms. For GO, the analysis used the package evoGO (version 0.1.0) to deprioritize terms that are redundant or unspecific.

2.7. Minimum effective concentration analysis of pathway responses

For each WikiPathway and compound, the lowest concentration associated with significant enrichment (p < 0.05) was determined. For each compound concentration with a significant pathway enrichment, enrichment consistency at higher compound concentrations for the same pathway was tested; a pathway was considered to be consistently significant across concentrations, if the median of log₁₀-transformed p values at the given concentration and all higher concentrations was below log₁₀ (0.05). The MEC was determined to be the lowest compound concentration with consistent pathway enrichment and used for the calculation of prediction metrics (sensitivity, specificity and accuracy). Pathway directionality was determined from the fraction of up- and down-regulated pathway genes at the MEC and was visualized with the ComplexHeatmap R package (version 2.10.0).

2.8. hiPSC-CMs quality control with human protein atlas expression data

Single-cell data and cell-type gene expression specificity data were obtained from the Human Protein Atlas (HPA) (version 22) [27]. Signature genes were selected from cardiomyocyte cell type enriched and group enriched genes with a least baseline expression of 15 transcripts per million (TPM), and a least expression fold-change of 10 in comparison with the median of all other cell types. A transcriptome-wide quantile regression model was generated from log₂ transcripts per million from ScreenSeq data against HPA data. Gene expression differences between hiPSC-CMs and HPA cardiomyocytes were calculated as regression residuals. A cardiomyocyte expression score was then calculated as the difference between the mean residuals of signature and non-signature genes with a least baseline expression of 15 TPM. Genes were grouped by selected evoGO-simplified GO terms with recursive inclusion of child terms and visualized with the ComplexHeatmap R package.

2.9. ScreenSeq cluster signature analysis

The Seurat R package (version 4.0.5) was used for comparison clustering, signature gene signature identification and dimension reductions. Genes significant in at least three comparisons in the complete dataset were used as background. Log₂ fold-changes vs. intra-plate DMSO samples were used as input. PCA was applied to all features, and comparison neighbors were identified considering the first 20 principal components (k = 9). Clusters were identified by shared nearest neighbor modularity optimization-based clustering (resolution = 1) and visualized by Uniform Manifold Approximation and Projection (UMAP). Signature gene identification was based on the Wilcoxon Rank Sum test (p < 0.05), by comparing each cluster with the pool of all other clusters, and with the very low differentially expressed gene (DEG) clusters 1 and 2. Signature features were filtered for an absolute fold-change of at least 1.25, and a log₂ fold-change with at least one standard deviation distance from 0. Enrichment of WikiPathways per cluster signature was performed with the clusterProfiler R package (version 4.2.2). For significant pathways (p_{Benjamini-Hochberg} < 0.01) represented by at least two features in the signature, the regulation direction was calculated as a fraction of up-regulated vs. total regulated genes. Pathway similarity was calculated as pairwise Jaccard index, and the rrvgo R package (version 1.6.0) was used to derive two-dimensional space coordinates.

2.10. Software

RStudio Server (2022.02.0 Build 443) running R (version 4.1.2, 2021–11-01) on an x86_64-conda-linux-gnu (64-bit) platform, Inkscape (version 1.1).

3. Results

3.1. Cardiotoxicity assessment with Ca²⁺ transient and high-content imaging assays

High-content imaging (HCI) and Ca²⁺ transient (CaT) analysis of 33 cardiotoxicants and 9 non-cardiotoxicants (Table 1) was used to validate the suitability of human-induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) for high-throughput cardiotoxicity prediction. Cardiotoxicants tested included those that caused either functional or structural cellular changes and those that elicited both outcomes (Tables 2,3 and 4), and hiPSC-CMs were dosed in triplicate at an 8-point dose response range covering total human plasma C_max concentrations (Table 1). CaT analysis (peak amplitude, frequency, width and decay time) was performed after acute and 24 h treatments to monitor functional alterations, while HCI analysis (cell count, cellular ATP, mitochondrial mass, MMP, cellular Ca²⁺ levels, DNA structure and nuclear size) was performed after 24 h treatments to assess structural changes [28].

Table 2. Summary of functional cardiotoxicants, reported functional effects (-ve: negative; +ve: positive) and in vivo cardiovascular effects. CHF: congestive heart failure, LVEF: left ventricular ejection fraction, RVEF: right ventricular ejection fraction, TdP: torsades de Pointes.

Compound	Functional effects	Functional effects reference	In vivo cardiovascular effects	Cardiac effects reference
Amitriptyline	-ve inotrope, -ve chronotrope, -ve dromotrope	[29]	Acute myocarditis, dilated cardiomyopathy, QT prolongation, CHF	[30–34]
Atenolol	-ve inotrope, -ve chronotrope, -ve dromotrope	[35,36]	Cardiac failure, bradycardia, hypotension, QT shortening	[17,144,149]
Bepridil	-ve inotrope, -ve chronotrope, -ve dromotrope	[37]	Ventricular arrhythmia, TdP, QT prolongation	[34,38,39]
Cisapride	+ve inotrope	[40]	Ventricular tachycardia, ventricular fibrillation, TdP, QT prolongation, cardiac arrest and sudden death	[34,41–43]
Digoxin	+ve inotrope, -ve chronotrope, -ve dromotrope	[44,45]	Cardiac rhythm disturbances, atrioventricular block, QT prolongation and shortening	[46–49,149]
Diltiazem	-ve inotrope, -ve chronotrope, -ve dromotrope	[37,50]	Bradycardia, hypotension, QT shortening	[17,18,37,50,51,144]
Dobutamine	+ve inotrope, +ve chronotrope	[52]	Increase in heart rate, QT prolongation and shortening	[34,53,54,144,149]
Dopamine	+ve inotrope, +ve chronotrope, +ve dromotrope	[52,55,56]	Cardiac arrhythmias, hypotension, angina pain	[57]
Epinephrine	+ve inotrope, +ve chronotrope, +ve dromotrope	[56,58]	Increased heart rate, increased RVEF and LVEF	[59,60]
Levosimendan	+ve inotrope	[61]	Increased heart rate, systolic blood pressure, fractional shortening and ejection fraction	[62]
Lidocaine	-ve inotrope, -ve chronotrope, -ve dromotrope	[63]	Bradycardia, cardiac arrest, QT shortening	[64,65,149]
Nifedipine	-ve inotrope, -ve chronotrope, -ve dromotrope	[37,50]	Hypotension, increased angina, myocardial infarction, CHF (rare)	[66,67]
Propranolol	-ve inotrope, -ve chronotrope, -ve dromotrope	[36]	Bradycardia, cardiac arrest, LV dysfunction, hypotension, QT prolongation	[15,16,49,68]
Sotalol	-ve inotrope, -ve chronotrope, -ve dromotrope	[36]	Life threatening ventricular tachycardia associated with QT prolongation, TdP	[11,16,49,68,69,142]

CaT analysis was validated with Ca²⁺ channel blockers, negative chronotropic adrenergic antagonists, positive chronotropic adrenergic agonists and QT-prolonging tyrosine kinase inhibitors (TKIs), which caused the expected reduction of Ca²⁺ wave amplitude, decreased frequency, increased frequency and increased peak width, respectively (Supplementary Figure S1A-D). Anthracyclines depleted Ca²⁺ waves in a time-dependent manner (Supplementary Figure S1E), which is consistent with their structural mechanism of action (MoA) and leads to a cumulative loss of functionality.

Dose response curves were applied to all individual assay readouts to determine the minimum effective concentration (MEC, mean value exceeding the vehicle control limits) for each readout in relation to C_max (Figure 1(a), Supplementary Data 1). Cardiotoxicity predictions using dynamic concentration thresholds as multiples of C_max were performed based on HCI readouts, CaT readouts or the combination of both assay types (Figure 1(b-c), Supplementary Data 2). The highest cardiotoxicity predictions were obtained using concentration thresholds between 10x and 25x C_max, by either the combination of CaT and HCI, or CaT alone. The combination of CaT and HCI with a fixed concentration threshold of 10x C_max resulted in 100% specificity (true negative: 9/9), 79% sensitivity (true positive: 26/33), and 83% accuracy. Elevated concentration thresholds allowed the true positive classification of rosiglitazone (CaT: 12x C_max), epinephrine (CaT: 13x C_max), and rofecoxib (HCI: 23x C_max), but resulted in false-positive classification of sildenafil (HCI: 14x C_max, CaT: 18x C_max). Although CaT analysis was predominantly superior over HCI, integration of the HCI method in the prediction improved the positive classification of idarubicin (HCI: below 1x C_max, CaT: 14x C_max) and facilitated the positive classification of rofecoxib (Figure 1(a)). Acute CaT analysis was more efficient at the detection of positive chronotropic compounds (epinephrine and isoproterenol) than 24 h CaT analysis, whereas 24 h CaT analysis correctly classified mitomycin C and bortezomib, which were classified as false negative in the acute treatment. Overall, the highest sensitivity (88%) and accuracy (88%) were obtained by combining HCI and CaT (acute and 24 h) at a 25x C_max concentration threshold with a single false-positive prediction (Figure 1(c)). Together, our data show that the combination of CaT at two time points (acute and 24 h) and HCI provides the best sensitivity for cardiotoxicity prediction.

Figure 1. Characterization of hiPSC-CM compound responses with HCI and CaT assays. a) Compound response map. Compounds are shown in columns grouped by cardiotoxicity. Assays readouts are shown in rows grouped by assays system and treatment duration. The MEC and the response direction (blue: up, red: down) were determined for each readout. MECs below 10x or 25x C_max are indicated by intense or light shading, respectively. b) Cardiotoxicity classification metrics. Compounds were classified as cardiotoxicants if the MEC of at least one readout in a) was below a dynamically selected C_max threshold. Compound cardiotoxicity predictions were then classified as false/true positive/negative based on Table 1. Prediction metrics are shown for the use of HCI, CaT or the combination of both assays. c) Cardiotoxicity prediction metrics at fixed 10x and 25x C_max thresholds indicated as vertical gray lines in b).

3.2. Establishment of high-throughput transcriptomics for cardiotoxicity risk assessment

To explore high-throughput transcriptomics for cardiotoxicity risk assessment, the transcriptional responses of hiPSC-CMs to 24 h treatments with a structural cardiotoxicant (imatinib), a functional cardiotoxicant (lidocaine), a combined structural and functional cardiotoxicant (sunitinib), and a non-cardiotoxicant (enalapril) were characterized by high-throughput RNA plate sequencing (ScreenSeq) (Figure 2(a), Supplementary Data 3). hiPSC-CMs derived from three independent hiPSCs were cultured on separate 384-well plates and dosed in triplicate at an 8-point dose response range or treated with DMSO (control). For the assessment of positional and plate effects, three intra-plate treatment repetitions were performed on each plate, and three replicate plates were processed for the first hiPSC-CM batch ScreenSeq identified transcripts for 6000 to 10,000 genes per sample with at least five unique molecular identifier (UMI) counts, without significant batch or plate effects (Figure 2(b)). Principal component analysis (PCA) of control samples (DMSO) revealed that the hiPSC-CM origin (hiPSC batch) was the main contributor to variance in baseline gene expression, whereas there was no observed technical plate effect (see plates 1–3) (Figure 2(c)). Cardiomyocyte (CM) signature genes from the Human Protein Atlas [137] were expressed at similar levels in hiPSC-CMs from different origins, including CM-specific components of contractile fibers, mitochondrial and membrane proteins such as gap junction subunits and transporters (Supplementary Figure S2A-E). Hence, all tested hiPSC-CMs were similarly CM-like despite the observed batch effects on sample clustering. The three cardiotoxicants induced reproducible dose response curves of DEGs (Figure 2(d), Supplementary Data 4) and dose response trajectories in PCA plots (Figure 2(e)) across hiPSC-CM batches and assay plates, whereas enalapril elicited nearly no transcriptional response. Sunitinib, lidocaine and imatinib altered pathways related to hiPSC-CM identity and functionality, mitochondrial, cholesterol and glycolytic metabolism, protein folding and turnover stress responses and signaling pathways relevant for hiPSC-CM viability (Figure 2(f), Supplementary Data 5, 6). This is in agreement with the well-characterized effects of sunitinib and imatinib on mitochondrial integrity [138,139] these structural toxicants, but not lidocaine, also altered mitochondrial pathways. Technical variability was observed only at the highest imatinib concentration due to severe cytotoxicity, and samples with compromised quality due to viability loss were excluded from pathway enrichment. In summary, ScreenSeq facilitates integrated model system validation and characterization of cardiotoxic dose responses in hiPSC-CMs.

Figure 2. Robustness of compound treatment responses in hiPSC-derived CMs. a) Scheme of CM origin and assay design. b) Number of identified genes by ScreenSeq with at least five UMI vs. transcript counts, including compound-treated and DMSO samples. c) PCA of normalized gene-level counts of the top 500 variable genes in DMSO control samples. d) DEGs over C_max-normalized compound concentrations. Lines represent plate means, and points represent DEGs per treatment sector. e) PCA of expression fold-changes of the top 500 variable genes vs. intra-plate DMSO controls. Principle components were calculated with all treatments, and individual treatments are highlighted in the individual plots. f) WikiPathway enrichment heatmap of compound treatments. Pathways with significant enrichment in at least one treatment below 25x C_max are displayed. The most significant false discovery rate (FDR) value in the tested concentration range is shown for each treatment sector.

3.3. High-throughput transcriptomics cardiotoxicity screening with ScreenSeq

ScreenSeq technology was applied for cardiotoxicity screening of the compounds tested by CaT and HCI (Figure 1), focusing on the 24 h treatment (Supplementary Data 3). Compound-induced differential gene expression was calculated in comparison with DMSO-treated controls. Cardiotoxicants had a stronger effect on differential gene expression than non-cardiotoxicants at comparable concentrations in relation to C_max values (Figure 3(a), Supplementary Data 7). Shared nearest neighbor clustering by gene fold-changes vs. DMSO-treated controls and UMAP representation segregated compound treatments into 13 clusters (Figure 3(b-d), Supplementary Figure S3, Supplementary Data 8). Most non-cardiotoxicants and low concentrations of cardiotoxicants, both associated with very low DEGs, clustered separately from the majority of cardiotoxicant treatments (clusters 1, 2). TKIs (cluster 3), α/β-adrenergic agonists (cluster 5) and DNA-damaging agents (clusters 11, 12) were grouped according to the distinct modes of action of their respective compound classes. A mixed group of compounds targeting various channels and receptors established a cluster group (clusters 7–9), in which each individual cluster was not strictly associated with the primary compound target. Three compound-specific clusters were formed by sunitinib (cluster 10), amphotericin b (cluster 6b) and the proteasome inhibitor bortezomib (cluster 13). Mixed cardiotoxicants and non-cardiotoxicants with weak effects on gene expression were dispersed across clusters 4 and 6. Very high concentrations of several compounds were grouped together in a high DEG cluster (10b), potentially reflecting a general cell stress response.

Figure 3. Differential gene expression and comparison clustering in the ScreenSeq transcriptomics cardiotoxicity screen. a) DEGs (fill color) detected per compound (y axis) by tested concentrations normalized to C_max values (x axis). Compounds were grouped by cardiotoxicity classification. b) Shared nearest neighbor clustering and UMAP organization by expression fold-changes of all genes regulated in at least three comparisons vs. intra-plate DMSO controls. Treatment conditions associated with each cluster are indicated. c-d) Mapping of DEGs (c), representative compounds (d), HCI data (e) and CaT data (f) on the UMAP from B. Fill colors indicate the number of DEGs (c), discrete compounds concentrations (yellow = low, red = high), or assay readouts normalized to intra-plate DMSO-treated controls (e, f).

HCI and CaT readouts were overlaid on ScreenSeq analysis clusters to assess the relation between compound-induced phenotypes and transcriptome state (Figure 3(e-f), Supplementary Data 9). For most HCI and CaT readouts, treatments with particularly high or low values were distributed across several clusters, showing that the biological states interrogated by HCI and CaT were not the major drivers of cluster formation. Such readouts were calcium, mitochondrial mass, MMP, nuclear size (HCI) and CaT readouts (amplitude, frequency, decay time and peak width). Increased contraction frequency, for example, mapped to the α/β-adrenergic agonist cluster (5), and isoproterenol, an α/β-adrenergic agonist that showed no clear frequency increase at 24 h treatment, did not localize to this cluster, arguing that cluster 5 is characterized by positive chronotropy; however, increased frequency was also apparent in the clusters 8 (diltiazem), 9 (verapamil) and 13 (bortezomib); hence, the impact of positive chronotropy in the global transcriptome state is not sufficient to cause condition aggregation into a single cluster. There were several exceptions of cluster-specific readouts. First, cell count and cellular ATP, correlated readouts of viability, were moderately reduced in various clusters, and cell count was severely reduced in cluster 10b. Hence, it is likely that a general toxicity response associated with loss of viability underlies cluster 10b. Second, reduction of ATP in comparison with cell count (HCI), indicative of energy depletion, was specific to cluster 10, suggesting a contribution of energy depletion to the distinct clustering of sunitinib, which is known to suppress mitochondrial energy metabolism [14,140,141]. Third, alterations of DNA structure (HCI) were highly pronounced in cluster 12, implying that structural DNA alterations induced by high concentrations of anthracyclines are a strong contributor of transcriptional identity.

ScreenSeq data clustering grouped compounds targeting ion channels and receptors (clusters 7–9), but the exact compound grouping suggests that the primary compound target is not the major clustering determinant (Figure 3(b)). Specifically, cluster 7 contained the two active β-adrenergic agonists (sotalol, propranolol), a Ca²⁺ channel blocker (nifedipine), a Na⁺ channel blocker (lidocaine), and an inhibitor of dopamine D2 and serotonin 2A receptors (clozapine). Cluster 8 contained two Ca²⁺ channel blockers (bepridil and diltiazem), one serotonin-norepinephrine reuptake inhibitor (amitriptyline) and a serotonin-4 (5-HT4) receptor agonist (cisapride). Cluster 9 contained a Ca²⁺ channel blocker (verapamil) and a mixed Na⁺/Ca²⁺ channel blocker and noncompetitive α/β-adrenergic inhibitor (amiodarone). Hence, the four Ca²⁺ channel blockers contained in our study are segregated into three groups. Despite the primary target diversity, most of the compounds cause CaT amplitude depletion and/or Ca²⁺ peak widening, indicating similar functional impairment. Signature gene extraction (Supplementary Data 10) and enrichment analysis of signature genes (Figure 4, Supplementary Figure S4, Supplementary Data 11) revealed that clusters 7–9 were characterized by largely overlapping biological responses: repression of glycolytic metabolism and mitochondrial respiration machinery, and elevated expression of DNA damage response pathways, and they differed in particular response pathways (7: glycogen metabolism, 8: TGF-β signaling and focal adhesion components, 9: cholesterol and heat shock proteins). Cluster 10 (sunitinib) was distinct by the absence of DNA damage response regulation, but shows a pronounced repression of additional mitochondrial components and induction of sterol regulatory element-binding protein (SREBP)-mediated cholesterol biosynthesis. The low viability/ATP-associated sub-cluster 10b was enriched in terms related to cell death (p53 pathway, ferroptosis) and processes indicative of low energy state (autophagy, AMPK signaling). The α/β-adrenergic agonist cluster (5) showed an opposite pathway regulation of the adrenergic antagonist-containing cluster (7), with a pronounced induction of glycolytic metabolism. The anthracycline/DNA damage clusters (11, 12) were different in terms of pathway regulation: while the low anthracycline concentration cluster 11 showed elevated activity of genotoxicity pathways, the high anthracycline concentration cluster 12 also exhibited an induction of mitochondrial respiratory chain components, mitochondrial metabolic processes (TCA cycle, fatty acid β-oxidation) and an induction of glycolytic metabolism, which is consistent with an adaptive response to impaired mitochondrial functionality. The proteasome inhibition cluster (13) showed a compensatory proteasome induction. Amphotericin B cluster (6b) regulated unfolded protein response (UPR), ER stress and fatty acid metabolism pathways consistent with membrane stress. The remaining clusters showed no (3, 4) or few disperse (6) enriched terms. In summary, ScreenSeq complements the HCI and CaT approaches by providing mechanistic information on compound activities and cellular responses.

Figure 4. Pathway enrichment landscape in ScreenSeq clusters. Gene signatures were identified for clusters from Figure 3(b) by comparison of cluster gene expression against the pool of all other clusters, and against the very low DEG clusters (1, 2). Significant WikiPathways were identified (p < 0.01) and organized into a similarity space based on the pairwise Jaccard index describing pathway gene intersection. Enrichment results for different clusters are split into facets, and the cluster number is shown in the top left corner of each facet. The direction of significantly regulated pathway genes is indicated by the fill color. The pathway enrichment significance is indicated by the circle size. The position of all pathways in the two-dimensional space is indicated as gray background points. Pathways important for biological interpretation are labeled with simplified names. All pathway enrichment results are shown in Supplementary Figure S4.

GNG: gluconeogenesis, ETC: electron transport chain, OXPHOS: oxidative phosphorylation, TCA: tricarboxylic acid, UPR: unfolded protein response, ER: endoplasmic reticulum, DDR: DNA damage response, HSP: heat shock protein, CI/III/IV: respiratory complex I/III/IV, Ub: Ubiquitin.

3.4. Cardiotoxicity prediction with ScreenSeq analysis

The DEGs in compound-treated cells were analyzed for pathway enrichment, and MECs for significant pathway enrichment, normalized by compound-specific C_max values, were determined (Supplementary Figure S5, Supplementary Data 12, 13). Enrichment was simplified by eliminating pathways with redundant enrichment pattern, focusing on representative physiological and pathological processes (Figure 5(a)). The enriched pathways were grouped into nine functional clusters (Figure 5(a)). A general cardiomyocyte functional cluster covering contractility and representative of viability-associated pathways (IL-18 signaling, VEGFR signaling) and cardiac disease states was enriched by most structural and functional cardiotoxicants at low concentrations. Glycolysis, gluconeogenesis and mitochondrial pathways (electron transport chain; ETC, oxidative phosphorylation; OXPHOS) were highly responsive to a similar spectrum of compounds. Specific aspects of mitochondrial processes (ETC complex assembly, tricarboxylic acid (TCA) cycle, fatty acid beta-oxidation) responded to a smaller subset of compounds. Finally, compound-specific pathways with high toxicology relevance included genotoxic stress, unfolded protein response, NRF2-mediated oxidative stress response and inflammatory response. Selected compounds affected cell cycle and differentiation-related pathways and ribosome homeostasis. Overall, both functional and structural cardiotoxicants caused more frequent pathway enrichment than non-cardiotoxicants.

Figure 5. Cardiotoxicity prediction with ScreenSeq analysis. a) The heatmap indicates if the indicated pathways (rows) are significantly regulated (Fisher exact test FDR < 0.05 with concentration consistency filtering) in hiPSC-CMs by the indicated compound treatments (columns). Enrichment by concentrations below 10x or 25x C_max are indicated by color with or without asterisk, respectively. The direction of significantly regulated pathway genes is indicated by the fill color at the MEC. The top annotation bar indicates the classification of compounds by cardiotoxicity. The side annotation bar indicates pathway grouping by biological function. A complete heatmap of all enriched WikiPathways is provided in Supplementary Figure S5. A table for the conversion of simplified pathway names to official WikiPathway names is provided in Supplementary Data 6. b) Sensitivity, specificity and accuracy (y axis) of cardiotoxicity prediction with varying concentration thresholds (x axis). Significance of any of the pathways in A with (Complete) or without (-Mitochondria) mitochondrial terms were used for prediction. c) Cardiotoxicity prediction metrics at fixed 10x and 25x C_max thresholds indicated as vertical gray lines in b).

GNG: gluconeogenesis, ETC: electron transport chain, OXPHOS: oxidative phosphorylation, TCA: tricarboxylic acid, UPR: unfolded protein response.

Cardiotoxicity prediction based on pathway enrichment MEC values was investigated, with or without inclusion of mitochondrial terms (Figure 5(b,c), Supplementary Data 14). The most accurate cardiotoxicity prediction was obtained without mitochondrial terms and a concentration threshold between 20x and 25x C_max (true negative: 9/9, true positive: 27/33, specificity: 100%, sensitivity: 82%, accuracy: 86%). Reduction of the concentration threshold to 10x C_max reduced the accuracy by false-negative classification of epinephrine. Inclusion of mitochondrial terms reduced the specificity by false-positive classification of acetaminophen. Mitochondrial terms were not required for the correct classification of any of the cardiotoxicants. A small group of cardiotoxicants (atenolol, levosimendan, rofecoxib, rosiglitazone, 3’-azido-3’-deoxythymidine (AZT), cyclophosphamide) was consistently classified as false negative (Figure 5(a)). These cardiotoxicants had either a low amount of DEGs in their respective toxicity class (atenolol, levosimendan, rofecoxib and AZT) or showed DEGs exclusively at the highest tested concentration (rosiglitazone, cyclophosphamide) (Figure 3(a)). Comparison of ScreenSeq cardiotoxicity prediction with HCI and CaT revealed that ScreenSeq facilitated the true positive classification of fluorouracil as cardiotoxicants and reduced the MEC of several other cardiotoxicants (Figure 6(a)).

Figure 6. Summary of cardiotoxicity prediction. a) Cardiotoxicity prediction results for HCI/CaT and ScreenSeq are derived from Figures 1 and 5 respectively. The numbers indicate the concentrations in relation to C_max at which the respective compounds are classified as positive. b) Sensitivity, specificity and accuracy (y axis) of cardiotoxicity prediction with varying concentration thresholds (x axis). Significance of at least one pathway from Figure 5(a) without mitochondrial terms or at least one assay readout from Figure 1(a). c) Cardiotoxicity prediction metrics at fixed 10x and 25x C_max thresholds indicated as vertical gray lines in b).

4. Discussion

In this study, three high-throughput technologies for cardiotoxicity prediction were utilized to analyze the mechanistic response of 42 pharmaceutical drugs, including 33 cardiotoxicants on hiPSC-CMs as the in vitro model system for human cardiomyocytes. These were HCI with six structural readouts and cellular viability, CaT analysis with four functional readouts, and the RNA-Seq method ScreenSeq with 6000–10,000 gene expression readouts per sample. All methods were applied after 24 h treatment except for CaT analysis, which was also applied after acute treatment to assess immediate functional responses. hiPSC-CMs recapitulated expected structural and functional toxicity mechanisms, validating their use as in vitro model system to detect and characterize modes of toxicity. Significant cardiotoxicity prediction metrics were obtained with a combination of HCI and CaT analysis combining both assay time points with fixed minimal effective concentration thresholds between 10x and 25x C_max (10x C_max: 100% specificity, 79% sensitivity, 83% accuracy; 25x C_max: 88% specificity, 88% sensitivity, 88% accuracy). Independently, CaT performed better at cardiotoxicity prediction than HCI, even when focusing on structural cardiotoxicants alone, which suggests that, at sufficiently long incubation times, various structural cardiotoxicants might affect cardiomyocyte functionality and be reliably classified as cardiotoxicants by CaT. ScreenSeq was established as a robust method to characterize cardiotoxicity-related pathway activation and compound mechanism of action, with pathway-based cardiotoxicity prediction performance comparable to the combination of HCI and CaT analysis (10x C_max: 100% specificity, 79% sensitivity, 83% accuracy; 25x C_max: 100% specificity, 82% sensitivity, 86% accuracy). Together, HCI, CaT and ScreenSeq provided the best cardiotoxicity prediction metrics (10x C_max: 100% specificity, 82% sensitivity, 86% accuracy; 25x C_max: 89% specificity, 91% sensitivity, 90% accuracy). Future compound set expansions will allow the refinement of classification thresholds and the establishment of optimal method and readout selections.

Previous studies have also used hiPSC-CMs as an in vitro model for cardiotoxicity risk assessment. Pointon et al. showed that hiPSC-CMs were suitable to detect drug-induced changes in cardiomyocyte contraction by testing 31 inotropic and 20 non-inotropic compounds using a video-based system to detect changes in cell contraction and a calcium transient assay. An IC₅₀ value of 50 µM as threshold for the peak count parameter reported the optimal sensitivity and specificity of 87% and 70%, respectively, while the same threshold for the transient profile parameters reported a sensitivity of 90% and specificity of 50% [17]. Ando et al. used the multi-electrode array system to investigate the effects of 57 drugs with various clinical torsadogenic risks on extracellular field potentials. Therein, drugs were categorized as high, intermediate or low risk according to the degree of field potential duration prolongation, which was compared to the torsadogenic risk categorization in CredibleMeds®. A concordance analysis of 36 drugs that showed either low or high risk gave a sensitivity of 81%, a specificity of 87% and an accuracy of 83% [142]. Dyballa et al. tested 92 compounds with known human cardiac liabilities in a high-throughput Ca²⁺ transient assay at three different time points: 5 min, 30 min and 90 min. In their study, a compound showing a significant effect in any of the features at any of the time points was classified as positive. This resulted in a sensitivity of 58%, specificity of 55% and accuracy of 57%. In the same study, the zebrafish model showed an improved sensitivity of 87%, a difference that, according to the authors, could be explained by the 13 cardiotoxicants not detected by hiPSC-CMs, due to, at least in part, metabolic reasons [15]. Among these, three compounds were investigated in the present study: doxorubicin, rofecoxib and sotalol. Doxorubicin, a structural/functional cardiotoxicant, was correctly identified here in all assay systems using both 10x C_max and 25x C_max thresholds. Rofecoxib and sotalol, a structural and a functional cardiotoxicant, respectively, were correctly identified in the HCI assay (25x C_max threshold) and the CaT assay and RNA-seq (10x C_max and 25x C_max thresholds), respectively. These results indicate that the toxic effects of the named drugs can effectively be detected using hiPSC-CMs as an in vitro model. In a more recent study, Palmer et al. used hiPSC-CMs in a metabolomics-based assay to identify both structural and functional cardiotoxicity. First, they used 66 drugs evaluated at a single, non-cytotoxic concentration to identify biomarkers that could identify functional and structural cardiotoxicants. Then, 81 drugs (including the original 66 drugs) were used to develop a model to predict cardiotoxicity based on changes in hiPSC-CM metabolism, which looked at the perturbation of four metabolites (arachidonic acid, lactic acid, 2’-deoxycytidine and thymidine) and cell viability. Depending on the metabolite-specific thresholds, the authors developed a model trained to maximize balance accuracy (BAC-trained model; 83% sensitivity, 90% specificity and 86% accuracy) or sensitivity (SEN-trained model; 90% sensitivity, 79% specificity and 85% accuracy) [143]. Finally, other publications have looked into the suitability of using hiPSC-CMs in cardiac microtissues for the prediction of cardiotoxicity. For example, Archer et al. used 3D cardiac microtissues with hiPSC-CMs, cardiac endothelial cells and cardiac fibroblasts in a high-throughput system to evaluate changes in MMP, endoplasmic reticulum integrity and cell viability. In their study, the authors tested a total of 29 drugs (15 structural cardiotoxicants and 14 non-structural cardiotoxicants) and reported that an IC₅₀ value of 10 μM as a cutoff provided optimal sensitivity (73%) and specificity (86%) for the detection of changes in cardiac structure. Pointon et al. used a similar model to detect drug-induced changes in cardiomyocyte contraction, by testing 29 inotropic and 13 non-inotropic compounds. An IC₅₀ value of 50 µM as threshold for the peak count parameter reported a sensitivity and specificity of 80% and 91%, respectively, for the detection of functional changes in cardiomyocytes [144].

In this study, high-throughput RNA-Seq (ScreenSeq) was established for the prediction and mechanistic characterization of compound-induced cardiotoxicity, and the synergism of ScreenSeq, HCI and CaT in detecting diverse cardiotoxicity mechanisms was demonstrated to predict overall cardiotoxicity risk. Cardiotoxicity predictions using our combined method compare favorably with the predictive value provided by other approaches that use hiPSC-CMs (see Table 5). In summary, hiPSC-CMs represent an important in vitro model system to identify drug-induced structural and functional toxicity [11] and this study further validates their use as an in vitro tool for cardiotoxicity prediction using HCI, CaT and high-throughput RNA-Seq to interrogate pathway-level cardiotoxicity responses and to derive mechanistic insight.

Table 5. Cross-comparison of prediction metrics and number of compounds tested (those overlapping with the present study are shown in brackets) with other publications that have used hiPSC-CMs as a tool for cardiotoxicity risk assessment using various methods and cutoffs.

Methods and cutoffs used	Sensitivity (%)	Specificity (%)	Accuracy (%)	Total compounds (overlapping)	Reference
HCI, CaT and ScreenSeq (25 x Cmax)	91	89	90	42	*
HCI and CaT analysis (25 x Cmax)	88	88	88	42	*
ScreenSeq (25 x Cmax)	82	100	86	42	*
HCI, CaT and ScreenSeq (10 x Cmax)	82	100	86	42	*
HCI and CaT analysis (10 x Cmax)	79	100	83	42	*
ScreenSeq (10 x Cmax)	79	100	83	42	*
Peak count (IC50 ≤ 50 µM)	87	70	NR	50 (21)	[17]
Transient profile (IC50 ≤ 50 µM)	90	50	NR	50 (21)	[17]
MEA system (CredibleMeds®)	81	87	83	36 (6)	[142]
CaT assay (any significant effect)	58	55	57	92 (16)	[15]
Metabolomics (BAC-trained model)	83	90	86	81 (26)	[143]
Metabolomics (SEN-trained model)	90	79	85	81 (26)	[143]
ER, MMP & ATP (IC50 ≤ 10 μM)	73	86	80	29 (18)	[140]
Peak count (IC50 ≤ 50 μM)	80	91	NR	42 (20)	[144]

*The present study, NR = not reported

The use of high-throughput RNA-Seq for the prediction and characterization of compound-induced cardiotoxicity in an in vitro model system was utilized in this study to provide an understanding of the underlying mode of action in addition to cardiac liability. Deriving insight into cardiotoxicity-related pathway activities thereby extends the selected information gained from classic HCI and CaT approaches, while providing comparable cardiotoxicity prediction performance based on pathway activities. The robustness of ScreenSeq was validated at several levels. First, read depth was confirmed to be stable across sequencing plates with 6000 to 10,000 genes per sample with at least five unique molecular identifier counts (Figure 2(b)). Second, the robustness of baseline expression was validated with inter-plate and intra-plate replicates, showing that plate and positional effects are negligible in comparison with the baseline transcriptome differences between the different iPSC lines or batches (Figure 2(c)). Third, the concentration-dependent gene regulation vs. intra-plate DMSO-treated controls was comparable in hiPSC-CMs from different iPSC lines or batches (Figure 2(d-f)), allowing the robust analysis of differential gene expression, sample identity and cardiotoxicity-related pathways. Fourth, compounds with similar mechanism of action were shown to elicit similar transcriptional states (Figure 3(b)) and similar pathway enrichment patterns (Figure 4, Figure 5(a)). In the same line, concentration-dependent mechanisms of the same compound were also distinguishable as exemplified by distinct clustering of low and high concentrations of DNA-damaging agents. Fifth, the directionality of pathway regulation was consistent with compound activities, such as the opposite regulation of glycolysis by β-adrenergic agonists and antagonists and the induction of genotoxic stress responses by genotoxic agents (Figure 5(a)). Importantly, pathways robustly regulated by the reference cardiotoxicants reflected expected mechanisms of toxicity. Sunitinib, a widely reported structural/functional cardiotoxicant, elicited transcriptional responses related to mitochondrial alterations, adaptation of energy metabolism, VEGF signaling, lipid homeostasis and protein folding (Figure 2(f)), which are consistent with its pleiotropic mechanisms of structural toxicity covered in the literature [14,140]. The structural toxicant imatinib, but not the functional toxicant lidocaine, elicited similar alterations in mitochondrial pathways [145]. Both types of cardiotoxicants consistently altered pathways related to general cardiomyocyte functionality, such as the striated muscle contraction pathway (Figures 2(f) and 5(a)). These general pathways were not affected by non-cardiotoxicants, indicating that their measurement by ScreenSeq could serve as robust readout of a physiological cardiomyocyte condition.

Additionally, heterogeneity and differentiation state are major sources of experimental variance that need to be considered in large-scale screens with hiPSC-derived model systems. Hence, the use of high-throughput RNA-Seq in both model system selection and model system quality control during screening was demonstrated in this study. By comparing a large panel of signature genes to in vivo expression levels, a similar cardiomyocyte identity of hiPSC-CMs from different iPSC lines and batches was observed (Supplementary Figure S2B-D). Furthermore, ScreenSeq allowed the comparison of reference compound responses in all hiPSC-CMs despite global gene expression baseline differences, showing only modest line-specific response differences restricted to lipid metabolism (statin, cholesterol, omega-9 fatty acids) (Figure 2(f)). Hence, applying ScreenSeq in model system selection supports the assessment of in vivo resemblance and extrapolation of crucial toxicity-related pathways. While this study shows largely homogeneous cardiotoxicant responses in cells derived from two donors, future efforts will include the screening of genetically diverse donor panels to cover the population heterogeneity of cardiotoxicant responses.

Another important aspect of this study was the detection of the key characteristics of cardiotoxicity through the combination of HCI, CaT and ScreenSeq [146]. These assays cover a specific spectrum of phenotypes relevant for cardiotoxicity prediction (Figure 6(a)), which explains the optimal cardiotoxicity prediction performance when combining all methods (Figure 6(c)). The CaT analysis detects changes in contraction behavior and Ca²⁺ wave dynamics commonly altered by functional cardiotoxicants. The HCI assay monitors changes in cell loss, global alterations in Ca²⁺ levels, mitochondrial state and DNA structure and, when paired with the measurement of gross ATP levels, it provides direct evidence for structural alterations. ScreenSeq quantifies gene-wise transcript levels and transcriptional responses, from which cell identity (differentiation state), pathway activities and signatures of compound targets and toxicity mechanisms can be derived, providing the highest level of mechanistic insights into molecular events. A combined application of HCI, CaT and ScreenSeq revealed that hiPSC-CMs were suitable to investigate all cell-autonomous, cardiac key characteristics (KCs) of the 12 previously proposed KCs of cardiovascular toxicity [146] (Figures 1(a), 4 and 5(a)). For instance, impaired regulation of cardiac excitability (KC 1), contractility and relaxation (KC 2) were detected by CaT readouts, which integrate contraction cycles and Ca²⁺ fluctuations (Figure 1(a), Supplementary Figure S1A-E). ScreenSeq detected altered regulation of the striated muscle cell contraction pathway and regulation of energy metabolism as frequently associated transcriptional responses (Figure 5(a)), which are likely of an adaptive nature. Induction of cardiomyocyte injury and death (KC 3) was measured as cell count and ATP content (Figure 1(a)). ScreenSeq distinguishes different types of injury through regulation of distinct response pathways (proteotoxic stress, nuclear genotoxic stress, NRF2-mediated oxidative stress response, ferroptosis) (Figures 4 and Figure 5(a)). Acute toxicity at high compound concentrations was further associated with autophagy and AMPK signaling (Figure 4, cluster 10b). Refinement of stress signatures with mechanistic compound sets will extend the resolution of this method in the future. Impaired mitochondrial functionality (KC 8) was monitored by HCI and ScreenSeq (Figures 1(a), Figure 5(a)). Induction of oxidative stress (KC 10) was observed as activation of the NRF2 pathway (Figure 5(a)). Compound-specific alterations in multiple pathways of paracrine factor and hormone signaling (KC 12) were further detected by ScreenSeq (PPAR-α, PPAR-γ, TGF-β, TNF-α, EGF/EGFR and VEGF/VEGFR) (Figure 4). Dyslipidemia (KC 7) is mainly associated with vascular toxicity; however, ScreenSeq analysis identified cholesterol regulation as a common lipid metabolism response to structural and functional cardiotoxicants in hiPSC-CMs (Figure 5(a)). The implication of this regulation on lipid homeostasis remains to be further investigated. General conclusions on the predictive power of individual key characteristics and pathways of cardiotoxicity require the expansion of the chemical space of the compound set. Notably, cell death-related pathways and autophagy were not common among the response mechanisms (Supplementary Figure S5), but future expansion of the compound panel will address their potential use as cardiotoxicity predictors.

In summary, HCI, CaT analysis and ScreenSeq together facilitate optimal cardiotoxicity prediction in hiPSC-CMs by integrating complementary levels of compound response information to maximize the coverage of cardiotoxicity mechanisms. Herein, a range of reference compounds with well-characterized structural (S) and functional (F) alterations were tested. Negative cardiotoxicants that were consistently classified as true negatives across all assays and thresholds (10x C_max and 25x C_max) included acyclovir, amoxicillin, acetylsalicylic acid, buspirone, enalapril, furosemide and tolbutamide. The two exceptions were acetaminophen (false positive with ScreenSeq at 10x C_max) and sildenafil (false positive with CaT and HCI at 25x C_max). Acetaminophen was detected as a false positive in transcriptomics due to the upregulation of mitochondrial pathways (Figure 5). Acetaminophen is known to be metabolized by the hepatic cytochrome CYP2E1 into a highly reactive toxic metabolite, N-acetyl-para-benzo-quinone imine (NAPQI), which binds to a number of mitochondrial proteins inducing mitochondrial oxidative stress [147]. The responses here observed cannot be attributed to NAPQI since the cardiac tissue does not typically express CYP2E1, therefore this points out to other mechanisms leading to mitochondrial dysfunction potentially caused by the parent drug. Sildenafil is a phosphodiesterase 5 inhibitor that relaxes vascular smooth muscle and does not affect the force of cardiac contraction or increase the risk of ventricular arrhythmia in vivo [148]. However, here sildenafil showed significant in vitro effects by reducing the peak amplitude and increasing the beating frequency, in addition to decreasing the cellular ATP while increasing the free cytosolic Ca²⁺. These data are consistent with another Ca²⁺ transient study published by Pointon et al., where sildenafil was also categorized as a false positive in hiPSC-CMs due to causing significant changes in a number of peak parameters at similar concentrations [17]. Positive cardiotoxicants that were consistently classified as true positives across all assays and thresholds included amiodarone (S/F), amitriptyline (F), amphotericin B (S), bepridil (F), bortezomib (S), cisapride (F), digoxin (F), diltiazem (F), doxorubicin (S/F), imatinib (S), lapatinib (S/F), mitoxantrone (S), nifedipine (F), sunitinib (S/F) and verapamil (S/F). Positive cardiotoxicants that were not detected in any of our assay systems included atenolol (F), AZT (S), and levosimendan (F).

Within the CaT assay, frequency, peak width and peak decay time were primarily affected by functional or mixed structural/functional cardiotoxicants (Figure 1). The frequency readout correctly identified expected positive (epinephrine, dopamine, dobutamine and isoproterenol) and negative (amitriptyline, sotalol, lidocaine, lapatinib, sunitinib and propranolol) chronotropic effects. Other compounds with reported negative chronotropic effects did not show a significant effect in the frequency parameter within the concentration range tested at any of the time points (atenolol, bepridil, nifedipine and amiodarone). Diltiazem and verapamil showed a significant increase in frequency at both time points at low concentrations only, while at higher concentrations they caused complete amplitude depletion (Figure 1). Finally, increased peak width and/or decay time was detected for various cardiotoxicants with inhibitory activity on the hERG channel (amitriptyline, cisapride, sotalol, propranolol, lidocaine, dobutamine, dasatinib, lapatinib, sunitinib) (Figure 1(a)). These compounds have been reported to cause QT prolongation in the literature (Tables 2, 3 and 4), except for lidocaine, which has been reported to induce QT shortening [149]. Ca²⁺ antagonists such as bepridil, diltiazem, nifedipine, verapamil and amiodarone (Table 1) showed an expected decrease in Ca²⁺ wave amplitude, as can be seen in Figure 1 [150,151]. A range of other cardiotoxicants with different primary targets also decreased amplitude in the absence of loss of cell viability (Na⁺ channel blocker lidocaine, β-adrenergic receptor blockers dobutamine, propranolol and sotalol), while, in some cases, a decrease in amplitude was paired with the loss of cell viability after 24 h (digoxin, bortezomib), implying that indirect suppression of Ca²⁺ dynamics is a common secondary effect of cardiotoxicants. Previous publications have provided evidence for beating suppression by these drugs using hiPSC-CMs [17,18].

Table 3. Summary of structural cardiotoxicants, proposed mechanisms of cardiac toxicity and in vivo cardiovascular effects. CHF: congestive heart failure, ER: endoplasmic reticulum, LVEF: left ventricular ejection fraction, TdP: torsades de Pointes, LVD: left ventricular dysfunction.

Compound	Proposed mechanisms of cardiac toxicity	Structural effects reference	In vivo cardiovascular effects	Cardiac effects reference
3’-azido-3’-deoxythymidine (AZT)	mtDNA and L-carnitine depletion, oxidative stress, inhibition of mitochondrial bioenergetic machinery	[70,71,153]	Cardiomyopathy, QT prolongation	[72,73]
Amphotericin b	Loss of ER integrity	[140]	Arrhythmia, atrial fibrillation, bradycardia, cardiac arrest, cardiomegaly, QT prolongation	[74–77]
Bortezomib	Induction of caspase-3/7 activity, oxidative stress	[78,79]	Arrhythmia, CHF, decreased LVEF, QT prolongation, hypotension	[80,81]
Clozapine	IgE-mediated pathways, cytokine-driven responses and oxidative stress	[82]	Myocardial ischemia, myocarditis, arrhythmia, QT prolongation, TdP, cardiomyopathy	[83–85]
Cyclophosphamide	Toxic metabolites cause depletion of antioxidants/ATP levels, altered contractility, damaged endothelium and enhanced pro-inflammatory/pro-apoptotic activities	[86,155]	Acute cardiac toxicity, CHF, myocarditis, myocardial necrosis, QT prolongation	[87,88]
Dasatinib	Inhibition of VEGF signaling and RAF/MEK/ERK pro-survival pathway, ER stress, mitochondrial ROS production	[89,90]	QT prolongation, CHF, LVD, myocardial infarction, cardiomyopathy, arrhythmia, cardiomegaly, pulmonary artery hypertension, pleural effusion	[91–94]
5-Fluorouracil	Citrate accumulation, depletion of high-energy phosphates, protein kinase C-mediated vasoconstriction, myocardial and endothelial cells apoptosis	[95]	Myocardial ischemia, CHF, coronary vasospasm, QT prolongation, (supra)ventricular tachycardia, LVD, cardiac fibrillation, arrhythmia, angina, myocardial infarction, cardiomyopathy	[88,96–99]
Imatinib	Mitochondrial toxicity, activation of ER stress, reduction of cellular ATP	[100,145]	Arrhythmias, CHF, decreased LVEF, edema, severe fluid retention, QT prolongation	[94,101,138]
Mitomycin c	Enhancement of the cardiac damage incurred by prior doxorubicin therapy (synergistic effect), oxidative stress	[102,103]	CHF	[102,104]
Mitoxantrone	Mitochondrial toxicity, ROS production	[105]	CHF, decreased LVEF, tachycardia, arrhythmia	[106,107]
Rofecoxib	Lipid peroxidation, oxidative damage	[108,109]	Cardiovascular thrombotic events, sudden death, unstable angina, ischemic attack/stroke, QT shortening	[110,149]
Rosiglitazone	Oxidative stress-induced mitochondrial dysfunction	[111,112]	CHF, myocardial infarction	[113–115]

Table 4. Summary of structural/functional cardiotoxicants, proposed mechanisms of structural cardiotoxicity, reported functional effects (-ve: negative; +ve: positive) and in vivo cardiovascular effects. CHF: congestive heart failure, ER: endoplasmic reticulum, MMP: mitochondrial membrane potential, LVEF: left ventricular ejection fraction, TdP: torsades de Pointes.

Compound	Proposed mechanisms of structural toxicity	Structural effects reference	Functional effects	Functional effects reference	Cardiovascular effects	Cardiac effects reference
Amiodarone	Loss of ER integrity and cell viability, MMP disruption	[140]	-ve inotrope, -ve chronotrope, -ve dromotrope	[150]	Arrhythmia, heart block, sinus bradycardia, CHF, ventricular fibrillation	[11,13]
Doxorubicin	Mitochondrial toxicity, lipid peroxidation, oxidative stress, loss of ER integrity and cell viability, MMP disruption	[116,117,140]	Biphasic inotropic effects	[118]	CHF, decreased LVEF, sinus tachycardia, myocarditis, QT prolongation, cardiomyopathy	[116,117,119–121]
Idarubicin	Mitochondrial toxicity, oxidative stress, loss of ER integrity and cell viability, MMP disruption	[122,140]	-ve inotrope	[123]	CHF, arrhythmia, cardiomyopathy, decreased LVEF, QT prolongation	[88,120,122]
Isoproterenol	Oxidative stress	[124–126]	+ve inotrope, +ve chronotrope, +ve dromotrope	[52,127,128]	Tachycardia, palpitations, ventricular arrhythmias, myocarditis	[129]
Lapatinib	Disruption of cardiomyocyte survival via EGF signaling inhibition, loss of ER integrity and cell viability, MMP disruption	[130,140]	-ve inotrope	[131]	Decreased LVEF, QT prolongation	[90,130,132,133]
Sunitinib	Disruption of cardiomyocyte survival via VEGF signaling inhibition, disruption of energy homeostasis, lipid accu-mulation, mitochondrial damage, loss of ER integrity	[14,140,141]	-ve inotrope	[14,131,134]	Decreased LVEF, CHF, QT prolongation, TdP, cardiomyopathy	[11,14,135,136,141]
Verapamil	Lipid accumulation	[11]	-ve inotrope, -ve chronotrope,-ve dromotrope	[37,50,144,151]	CHF, pulmonary edema, hypotension, ventricular fibrillation	[11,17,68,144]

Some structural cardiotoxicants such as 5-fluorouracil, AZT, cyclophosphamide, clozapine, dasatinib, rosiglitazone and isoproterenol were not detected by HCI alone. 5-fluorouracil interrupts DNA replication by inhibiting thymidylate synthase, and it is consistently categorized as a cardiotoxicant in vivo in the literature. However, the in vitro results have not always been successful at predicting such toxicities [16,18,140,143]. For instance, Guo et al. showed the challenges for obtaining correct predictions for 5-fluorouracil using impedance measurements in hiPSC-CMs [16,152]. Archer et al. were able to determine an IC₅₀ value for three parameters related to structural changes in cardiac microtissues (cellular viability (IC₅₀ = 27.5 µM), ER integrity (IC₅₀ = 60 µM) and MMP (IC₅₀ = 11.2 µM)), however, their 10 µM cutoff for optimal prediction metrics categorized 5-fluorouracil as a false negative [140]. Sirenko et al. characterized the beating profiles of hiPSC-CMs using intracellular Ca²⁺ flux readouts and were able to identify significant changes in the beating count, peak spacing and cell viability caused by 5-fluorouracil at concentrations < 30 µM, which allowed identifying 5-fluorouracil as a true positive [18]. In this study, both the CaT and the HCI assay failed to detect any significant changes at 10x C_max and 25x C_max thresholds, while ScreenSeq facilitated the identification of 5-fluorouracil as true cardiotoxicant at 10x C_max due to the genotoxicity and cell cycle entry pathways being significantly regulated, which agrees with 5-fluorouracil’s main mechanism of action (Figure 5).

In this study, no significant effects were observed for AZT at any end-point. It has been proposed that the cardiotoxicity shown by AZT may be mediated by disrupting the substrate supply of 5′-triphosphate for mt-DNA replication. Therefore, for mitochondrial biogenesis to be affected, long-term treatment with AZT may be required before mitochondria are depleted enough to significantly reduce ATP synthesis [153]. Clements et al., however, investigated the effects of AZT on the hESC-CMs electrophysiology (multielectrode array assay), beating (impedance assay) and subcellular structure (high content imaging) and could not identify any effects at the concentrations tested at 72 h incubation treatment [12].

Another interesting compound is cyclophosphamide, whose in vivo cardiotoxicity classification in the literature is not always unanimous. Archer et al. categorized cyclophosphamide as a false negative due to the IC₅₀ values of the structural parameters being ≥ 10 µM in cardiac microtissues, while Sirenko et al. found it to be a true positive through changes in some beating parameters but as a false negative using the cell viability assay at 24 h in hiPSC-CMs (cutoff of IC₅₀ < 30 µM) [18,140]. Conversely, Guo et al. reported cyclophosphamide as a true negative through impedance measurements in hiPSC-CMs due to negative or equivocal observations in the clinic for hERG inhibition, clinical QT prolongation and TdP/arrhythmia [16]. Here, cyclophosphamide was considered a positive structural cardiotoxicant due to the reported cardiac events in the literature [154]. Cyclophosphamide acts as a DNA cross-linking agent, mainly acting through its metabolite, phosphoramide mustard. Additionally, cyclophosphamide undergoes hepatic metabolism that results in the production of acrolein, a toxic metabolite believed to act on the myocardium and endothelial cells causing depletion of antioxidants/ATP levels, altered contractility and damaged endothelium and enhanced pro-inflammatory/pro-apoptotic activities resulting in cardiomyopathy [155]. The lack of hepatic metabolism in our in vitro system could protect against both metabolites and could explain the absence of detrimental structural effects observed in the HCI assay, as well as the low transcriptional response elicited by cyclophosphamide (Figure 6(a)). In future studies, this could be addressed by co-culturing cardiomyocytes and hepatocytes in a multi-organ human-on-a-chip system [156].

ScreenSeq improved the overall sensitivity by identifying various molecular mechanisms of structural toxicity, such as alterations in cardiac pathways, genotoxicity, ER stress and mitochondrial toxicity (Figure 4, Figure 5(a)). Genotoxicity pathway enrichment was associated with genotoxic compounds (mitoxantrone, mitomycin C, fluorouracil and idarubicin) and hence more specific than HCI chromatin readouts (Figure 1(a)). Several ER stress inducers (sunitinib, lapatinib, amphotericin B) elicited protein misfolding responses (UPR, proteasome degradation, Parkin-Ubiquitin system). Cardiotoxicants with known mitochondrial toxicity component targeted mitochondrial ETC/OXPHOS and respiratory complex assembly pathways (mitoxantrone, amiodarone, dasatinib, imatinib, doxorubicin, idarubicin and sunitinib). While transcriptional responses are a promising readout for mitochondrial damage, it should be noted that other non-cardiotoxicants with other organ toxicities could present mitochondrial liabilities, therefore signature refinement would be beneficial to specifically identify cardiotoxicity-associated mitochondrial damage. A clear example of this is acetaminophen, a well-known drug for liver injury but not cardiotoxicity which showed mitochondrial pathway regulation (Figure 5(a)) [147].

In summary, this study demonstrates the potential of HCI, CaT analysis and ScreenSeq in predicting compound-induced cardiotoxicity in hiPSC-CMs, and thereby serves as a proof of concept for an integrated approach aiming at maximum coverage of cardiotoxicity mechanisms. Expansion of the compound collection will further enhance the diversity of toxicity mechanisms, optimize the definition of relevant concentration ranges for cardiotoxicity classification, and facilitate the definition of refined transcriptomic signatures to predict compound targets and diverse toxicity mechanisms. The use of additional endpoints could further support pathway activities monitored with ScreenSeq and distinguish between primary and adaptive effects. The current in vitro model system is designed to cover cell-autonomous aspects of cardiotoxicity. Cardiac microtissues containing the three major cardiac cell types, cardiomyocytes, fibroblasts and endothelial cells, recapitulate better the cardiac organ and could therefore add predictive value to our assays [140,157]. Non-cardiomyocytes play an important role in cardiac physiology and their interactions are important not only in normal heart functionality but also in the development of disease phenotypes and cardiotoxicity [158,159]. For instance, it has been suggested that sunitinib’s cardiotoxicity is mediated through fibrotic processes [160]. In summary, compound set expansion and model system optimization will extend the application of combined HCI, CaT analysis and ScreenSeq in the integrated prediction of compound cardiotoxicity and mechanism of action.

5. Conclusion

This study introduces a mechanism-driven risk assessment approach combining structural (high-content imaging; HCI), functional (Ca²⁺ transience; CaT) and high-throughput RNA-sequencing (ScreenSeq) for the pre-clinical risk assessment of novel compounds in hiPSC-CMs. Together, HCI, CaT and ScreenSeq covered a broad spectrum of phenotypes and mechanisms relevant for cardiotoxicity prediction, providing excellent cardiotoxicity prediction metrics (10x C_max: 100% specificity, 82% sensitivity, 86% accuracy; 25x C_max: 89% specificity, 91% sensitivity, 90% accuracy).

Declaration of interests

A Rosell-Hidalgo, C Bruhn, P Walker, E Shardlow, R Barton, S Ryder, T Samatov, A Hackmann, G Aquino, M Fernandes dos Reis, V Galatenko, R Fritsch and C Dohrmann are employed at Cyprotex Discovery Ltd., where these assays are performed for our partners as a service. The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript. This includes employment, consultancies, honoraria, stock ownership or options, expert testimony, grants or patents received or pending, or royalties.

Reviewer Disclosures

Peer reviewers on this manuscript have no relevant financial or other relationships to disclose.

Author Contributions

A Rosell-Hidalgo: Paper writing, experimental work and data analysis; C Bruhn: Paper writing and data analysis; P Walker: Paper writing, experimental planning and data analysis; E Shardlow: Experimental work; R Barton: Experimental work; S Ryder: Experimental work and data analysis; T Samatov: Experimental work; A Hackmann: Experimental work; G Aquino: Experimental work; M Fernandes dos Reis: Data analysis; V Galatenko: Data analysis; R Fritsch: Data analysis and experimental planning; C Dohrmann: Experimental planning.

Author Note

^*These authors contributed equally to this work.

Acknowledgments

We thank V Kari and A Filipchyk for discussions and help in data processing.

Data availability statement

ScreenSeq data were deposited into the Gene Expression Omnibus database under accession number GSE244740 (GEO Accession viewer (nih.gov)).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/17425255.2023.2273378

Additional information

Funding

This paper was funded by Evotec International GmbH.