A predictive in vitro risk assessment platform for pro-arrhythmic toxicity using human 3D cardiac microtissues

Cardiotoxicity of pharmaceutical drugs, industrial chemicals, and environmental toxicants can be severe, even life threatening, which necessitates a thorough evaluation of the human response to chemical compounds. Predicting risks for arrhythmia and sudden cardiac death accurately is critical for defining safety profiles. Currently available approaches have limitations including a focus on single select ion channels, the use of non-human species in vitro and in vivo, and limited direct physiological translation. We have advanced the robustness and reproducibility of in vitro platforms for assessing pro-arrhythmic cardiotoxicity using human induced pluripotent stem cell-derived cardiomyocytes and human cardiac fibroblasts in 3-dimensional microtissues. Using automated algorithms and statistical analyses of eight comprehensive evaluation metrics of cardiac action potentials, we demonstrate that tissue-engineered human cardiac microtissues respond appropriately to physiological stimuli and effectively differentiate between high-risk and low-risk compounds exhibiting blockade of the hERG channel (E4031 and ranolazine, respectively). Further, we show that the environmental endocrine disrupting chemical bisphenol-A (BPA) causes acute and sensitive disruption of human action potentials in the nanomolar range. Thus, this novel human 3D in vitro pro-arrhythmic risk assessment platform addresses critical needs in cardiotoxicity testing for both environmental and pharmaceutical compounds and can be leveraged to establish safe human exposure levels.

www.nature.com/scientificreports/ 5d. Use of CellTracker dyes to label hiPSC-CM (green) and hCF (red) prior to assembly as well as immunostaining with antibodies against cardiac troponin I (cTnI) and vimentin (vim) after fixation and sectioning confirmed that hiPSC-CMs and hCFs self-assembled to be highly interspersed (Fig. 1D,E), consistent with our prior work in rat cardiac microtissues 35,36 . Timeline shows cardiomyocyte differentiation from human-induced pluripotent stem cells (hiPSCs) in high density 2D culture. Cardiac directed differentiation was achieved with Wnt activation at day 1 and inhibition at day 3 (see "Methods"). Cardiac phenotype, visually confirmed by beating cells, appeared between days 8 and 12. Cardiomyocytes differentiated from hiPSCs were used for the production of microtissues or were further purified with a lactate-based metabolic selection protocol. Microtissues self-assembled in microwells after 14-28 days of differentiation of hiPSC-CMs (i.e., without or with lactate purification) and addition of human cardiac fibroblasts (hCFs), and they were cultured in the presence of 1 Hz electrical field stimulation (estim).
(B) Schematic of three-dimensional (3D) cardiac microtissue generation shows non-adhesive agarose gels with cylindrical recesses with hemispherical bottoms that guide self-assembly. Cardiac microtissues were cultured for 6-8 days with 1 Hz pacing. (C) Phase contrast image shows consistent spherical microtissue formation after 5 days of 3D culture in all 35 microwells. Scale bar, 800 μm. (D) Confocal image shows a representative cardiac tissue with hiPSC-CM (green) and hCF (red) stained with CellTracker dyes. Scale bar, 200 μm. (E) Confocal image shows a representative cardiac troponin I (red), vimentin (green), and DAPI stained cryosection (10 μm thick) of a microtissue fixed after 7 days in 3D culture. Scale bar, 50 μm. (F) Fluorescence image of microtissues at 3.2 × magnification was obtained during optical mapping. Typically, the action potentials (APs) from 4-9 microtissues were recorded simultaneously. (G-I) Schematics of the AP metrics of that were defined (with units) as: (G) "excitability" (%) measured from the percentage of captured APs during 10 s duration of recording with 2 s pacing cycle length, (H) "stimulation time delay" (ms; stim delay) between stimulation pulse and evoked AP upstroke (dF/dt max ), "rise time" (ms) of AP, "AP duration" (ms) to 30%, 50%, and 80% repolarization (APD 30 , APD 50 , APD 80 ), "APD to the maximum repolarization rate" (ms; APD MxR ) defined as time between AP upstroke and the end of rapid repolarization marked by d 2 F/dt 2 max , "APD triangulation" (ms; APD tri ) defined as APD MxR -APD 50 , and (I) occurrence of "early afterdepolarization" (EAD) reported as (%) of microtissues showing EADs. www.nature.com/scientificreports/ Electrophysiological characterization of 3D cardiac microtissues. In order to evaluate excitability and sensitivity to pro-arrhythmic toxicants, APs were recorded from the human cardiac 3D microtissues with voltage sensitive dyes (Fig. 1F). Cardiac microtissues were paced at 2 s basic cycle length (0.5 Hz) and APs were recorded for 8-10 s. Human cardiac 3D microtissues showed robust fluorescence changes tracing V m signal of APs during spontaneous and paced beats (Movie S2). The spontaneous rate could not be reliably measured because spontaneous action potentials were very rare and often only one spontaneous beat per scan interval of 8 ~ 10 s could be acquired. The data analysis of APs was automated to reduce bias and achieve increased throughput (Supplemental Fig. S1). Pro-arrhythmic risk metrics were defined as depicted in Fig. 1G-I: excitability, stimulation time delay between stimulation pulse and peak AP upstroke (stimulation delay time, stimΔ), rise time of AP upstroke, AP duration (APD) to 30%, 50%, and 80% repolarization (APD 30 , APD 50 , APD 80 ), AP duration to end of maximum repolarization rate (APD MxR , determined with d 2 F/dt 2 max ), APD triangulation (APD tri = APD MxR − APD 50 ), and presence of EADs (% of microtissues showing EAD).
Biological variation was low across biological replicates and experimental time (Fig. 2). Example traces show reproducibility ( Fig. 2A), and an example of APD distribution across different microtissues from a single mold was normally distributed (Fig. 2B). The Kolmogorov-Smirnov test for 3 molds with a total of 105 cardiac microtissues indicated normal APD distribution (p = 0.92). In order to determine the statistical power of this 3D in vitro methodology, we examined the standard deviation as the variability in APD. In microtissues formed with hiPSC-CMs that did not undergo lactate purification (70.3 ± 3.1% cardiac troponin T positive for n = 4 HiPSC-CMs (without lactate purification) were seeded in microwells without or with 5% hCFs. After 6-8 days in 3D culture, microtissues were loaded with di-4-ANEPPS for optical mapping of membrane potentials (V m ) under paced conditions. (A) Representative V m traces from individual microtissues paced at 0.5 Hz in different molds (biological replicates) and different batches (experimental replicates). Values were averaged from all pixels of single microtissues. Only microtissues greater than 60 pixels were analyzed (typically ~ 170 pixels/microtissue). (B) APD 80 distribution from 35 microtissues in a single mold is normally distributed. (C) APD variation was assessed for technical, biological, and experimental replicates for APDs measured at 0.5 Hz. The standard deviation from beat-to-beat in the same microtissue (technical replicates) was 8.7 ± 4.8 ms (n > 70 microtissues per batch, n = 7 batches), from microtissue-to-microtissue in the same mold (biological replicates) was 20.9 ± 8.4 (n = 2-3 molds per batch, n = 7 batches), from mold-to-mold in the same batch (also biological variability) was 11.7 ± 8.9 ms (notably less than microtissue-to-microtissue variability; n = 3 molds per batch, n = 7 batches), and from batch-to-batch of differentiated hiPSC-CMs (experimental variability and another level of biological variability) was 33.1 ms (n = 7 batches). Note that the standard deviation reported for batch-to-batch variability does not itself have variation (and therefore no error bar is appropriate). (D) The APD restitution curve displays typical rate-dependent increases in APD 80  www.nature.com/scientificreports/ batches), the mean APD 80 was 160.0 ± 33.1 ms (Fig. 2C, bottom). The standard deviation from batch to batch ("between batches") was the largest of our biological variability metrics at > 30 ms. Microtissue-to-microtissue and mold-to-mold variability were smaller, and variability in technical replicates beat-to-beat (same microtissue) was < 10 ms (Fig. 2C). Action potential duration showed the expected physiological rate dependence at physiological heart rates (60-75 bpm; Fig. 2D). Power analysis of sample size suggested that n = 17 cardiac microtissues are required to detect a 10% change in APD with 95% confidence, suggesting that testing within a single mold is sufficient for evaluating small but meaningful changes in APD measured in vitro 58 . The sample size can be further reduced by using paired testing (n = 8 based on beat-to-beat variation) on the same mold before and after toxicant exposure. In order to screen cardiotoxicity from the same microtissues before and after toxicant exposure, APDs should be stable between successive recordings. We chose a 20 min incubation period to ensure that the test compound had ample time to diffuse into the 3D cardiac microtissue 59 , which is the limiting factor since typical ion channel blockers bind and act relatively quickly, in less than a minute (for example, E4031 has a binding time constant of 0.8 sec 60 .). Figure 2E,F shows comparable traces and group data for APD 80 in our model over 20 min, which enabled testing of several doses of toxicant with this acute exposure protocol on the same mold to perform paired t-tests. The addition of 5% hCF to microtissues upon formation improved not only the consistency of formation and clean edges of microtissues (data not shown), but also significantly increased the number of excitable microtissues by more than 10% (Fig. 2G), while slightly shortening APD by 12% (Fig. 2H).
Our model was predictive of expected arrhythmia risk in pilot experiments with known drugs such as I Kr blocker E4031 (data not shown), yet the APD in these microtissues was under 200 ms, and while within the range reported for hPSC-CMs 61 . this was shorter than APD typically recorded in human 62,63 . To improve the reliability and robustness of our arrhythmia model, we selected for and matured our hiPSC-CMs with metabolic-based lactate purification of cardiomyocytes (hiPSC-CM LP ). Lactate purification lengthened APD (Fig. 3A) Figure 3. Lactate purification of hiPSC-CMs prior to microtissue formation with 5% hCFs improves excitability and lengthens APD. (A) Representative V m traces from microtissues formed from hiPSC-CMs with (red) and without lactate purification (black) plus 5% hCFs. HiPSC-CMs from the same differentiation batch but grown in 2D culture without the lactate selection protocol for 28 days (to match the age of the cells) were used as the most appropriate control for comparison. (B) Percent of hiPSC-CM microtissues showing APs during 0.5 Hz pacing increases with lactate purification of hiPSC-CMs (n = 4 batches, minimum 2 molds per batch). (C) APD 80 distribution from 35 microtissues in a single mold with hiPSC-CM LP and 5% hCF. (D) The average APD 80 was 259.0 ± 42.0 ms (n = 201). APD 80 variation between batches is small with hiPSC-CM LP : APDs were measured at 0.5 Hz cycle length pacing. The standard deviation from beat-to-beat in the same microtissue (technical replicates) was 6.4 ± 6.5 ms (n > 35 microtissues per batch, n = 4 batches), from microtissue-tomicrotissue in the same mold (biological replicates) was 36.1 ± 21.4 ms (n = 2-3 molds per batch, n = 4 batches), from mold-to-mold in the same batch (also biological variability) was 31.8 ± 3.8 ms (n = 2-3 molds per batch, n = 4 batches), and from batch-to-batch, the standard deviation (another level of biological variability) was 19.3 ms (n = 4 batches). (E) Cumulative probability plot for APD 80 (n = 3-4 batches, at least 2 molds per batch) shows consistently longer APD 80 with hiPSC-CM LP (red) vs. hiPSC-CM (black). (F) Comparison of APDs from microtissues from the same experimental batch (each color) with and without lactate purification (n = 4 batches, ≥ 2 molds per batch). Values are means ± SD. * P < 0.05). Note that in batch 4, the microtissues generated from hiPSC-CM without lactate purification were not excitable. www.nature.com/scientificreports/ and improved the excitability of microtissues when directly compared to experimental controls that were agematched but not lactate purified (84.4 ± 6.5% for hiPSC-CM LP vs. 25.2 ± 25.8% for hiPSC-CM) (Fig. 3B). Similar to hiPSC-CM microtissues with hCFs, the hiPSC-CM LP microtissues with hCFs have a normal distribution (Kolmogorov-Smirnov test p = 0.580; Fig. 3C) and low variability (Fig. 3D). Batch-to-batch standard deviation was reduced by 42% compared to hiPSC-CM microtissues and beat-to-beat APD variation from the same microtissue (technical replicates) was also smaller, whereas microtissue-to-microtissue variability was similar, and mold-to-mold variability was slightly higher than variation in hiPSC-CM tissues ( Fig. 3D vs. Figure 2C). The cumulative probability plot for APD 80 shifted to the right with lactate purification (Fig. 3E), and APD was doubled from 130.6 ± 54 ms in hiPSC-CM microtissues to 263.1 ± 54.9 ms with lactate purification (Fig. 3F). All subsequent experiments were performed with microtissues generated from hiPSC-CMLP and 5% hCFs.

Model qualification: APD and arrhythmia test with known high-risk and low-risk potassium channel (I Kr ) blockers.
To validate the capability of our model to predict arrhythmia risk, we exposed microtissues to the pro-arrhythmic hERG channel blocker E4031. APD 80 was prolonged under 2 µM E4031 ( Fig. 4A,B, Movies S3, S4). Cumulative distribution plots show a rightward shift with E4031 treatment (Fig. 4C). When analyzing single batches of microtissues (n = 35), we confirmed that APD 80 had a normal distribution of values with Q-Q plots (Supplemental Fig. 2). EADs were present in 54% of microtissues with 2 μM E4031 exposure in the absence of isoproterenol (ISO; Fig. 4A Fig. 4C). E4031 prolonged APD in a dose-dependent manner, and in some batches, additional beta-adrenergic stimulation using 100 nM ISO was needed to trigger EADs (Supplemental Fig. 3). Importantly, ISO alone shortened APDs and did not evoke EADs without E4031 (Supplemental Fig. 4A). We further tested blockade of the outward potassium current I to (responsible for phase 1 repolarization) with 4-Aminopyridine (4-AP) and stimulation of L-type calcium channels (responsible for maintaining phase 2 plateau) with the agonist BayK8644, and both prolonged APDs (Supplemental Fig. 4B,C). These results support that cardiac microtissues using hiPSC-CMs and hCFs are an effective platform for screening pro-arrhythmic toxicants. We further investigated how well this human cardiac microtissue model identifies compounds used medically and already determined to be low risk for arrhythmia in clinically-validated doses in patients even when there is a known capacity to block the hERG channel. One such compound is ranolazine, which is a known hERG channel blocker that has therapeutic value as an anti-arrhythmic drug at mean steady-state blood serum concentrations (C max ) of 6 μM (2600 ng/mL) with 95% confidence limits of 1 and 14 μM (400 and 6100 ng/mL, www.nature.com/scientificreports/ respectively) 64, 65 . In this concentration range, ranolazine had a minimal effect on APD ( Fig. 4D-F), which is in stark contrast to the specific hERG (I Kr ) blocker E4031 ( Fig. 4A-C). Figure 5 shows a more complete assessment of the dose-dependent effects of ranolazine at 0, 2, 10, and 100 μM concentrations on eight AP metrics, and traditional dose response curves are shown in Supplemental Fig. 5A. Figure 5A shows representative examples of AP traces. Figure 5B shows a color map presentation of seven of the AP metrics from 35 microtissues, and group data are shown for each of them plus excitability in Fig. 5C. Ranolazine at therapeutic concentrations (2-10 μM) had statistically significant effects but very small on all AP metrics ( Fig. 5C; see Supplemental Table 1 (top) for complete statistical analysis). For example, average APD MxR changes were negligible (− 9 ms and + 15 ms at 2 and 10 μM, respectively) versus control 0 μM (Fig. 5C,F). However, at high concentration (100 μM ranolazine), APD MXR more than doubled (ΔAPD = 357 ms, Fig. 5C,E). Despite prolongation, no EADs were observed with ranolazine even in the presence of ISO (see Fig. 5A bottom trace; 1-2 molds per batch, 2 batches) in contrast to the highly selective hERG channel blocker, E4031 (Fig. 4A,B and Supplemental Fig. 3). 100 μM ranolazine reduced excitability from 97.1% at 0-10 μM to 85.1% at 100 μM, as some of the pacing beats were missed ( Fig. 5A-C) and was associated with increased stimulation delay and AP upstroke rise time (Fig. 5B-D) consistent with its therapeutic effects on I Na . Statistically significant changes are detected across all metrics and often at all concentrations of ranolazine (Supplemental Table 1), suggesting that this platform is highly sensitive.
Screening pro-arrhythmic toxicity of environmental toxicant exposure. We investigated whether acute exposure to the environmental pollutant bisphenol A (BPA) alters AP parameters using our human cardiac microtissue platform. BPA has previously been shown to affect the cardiovascular system 66 , but its specific effects on arrhythmia are not well characterized. The response of cardiac microtissues to BPA was complex, altering multiple AP parameters including excitability, stimulation time delay, APDs at different levels of repolarization, APD MxR , and APD tri in a dose-dependent manner (Fig. 6). A small but significant increase in APD metrics appeared at the lowest concentration of BPA (1 nM), but BPA gradually shortened APD metrics at higher con- Figure 5. Dose-dependent effects of ranolazine on AP manifested at high concentrations above the therapeutic window in hiPSC-CM LP with 5% hCF microtissues. (A) Representative V m traces from 0, 10, 100 μM ranolazine, and 100 μM ranolazine plus 50 nM ISO show AP stability at 10 μM and instability in microtissues exposed to 100 μM Ran with or without ISO . (B) AP metric color map shows each microtissue (rows) and seven AP metrics (major columns) with dose-dependent response at indicated concentrations (minor columns). The color displayed shows shifts as measured by the standard deviation (σ) from the average value of the AP metric under the control condition (0 μM). Increases appear in the warm color spectrum (progressing yellow-orangered) and decreases appear in the cool color spectrum (progressing cyan-blue-purple). Black indicates nonexcitable microtissues. When excitability is lost in a microtissue, other metrics could not be measured. www.nature.com/scientificreports/ centrations (Fig. 6B,C). This biphasic response was most pronounced for APD MxR , which increased by 10.5 ms with 1 nM BPA, was similar to no BPA with 10 nM BPA, and decreased by 14.2 ms and 57.5 ms with 100 and 1000 nM BPA (Fig. 6C,F). Similar to APD, a biphasic effect of BPA was observed in a fraction of the microtissues (~ 6 of 35) on stimulation time delay (Fig. 6B, red/orange boxes at 1 nM in StimΔ column), suggesting that complex interactions of BPA on multiple cardiac ion channels involved in both excitation and repolarization was revealed during acute exposure across biological replicates in our model. Interestingly, a divergence between APD 80 and APD MxR was observed at the highest concentration (1000 nM; Fig. 6E,F) due to the appearance of slowed repolarization near the end of AP repolarization (see 1000 nM trace in Fig. 6A), which indicates that very high BPA exposure may also alter inward rectifier current.
The changes in AP metrics were markedly different between ranolazine and BPA, as reflected in the color maps of AP metric changes in response to each compound (Figs. 5B, 6B, respectively). These drug response signatures are visualized in Fig. 7, where significant changes (red, increased; blue, decreased) and their magnitude (color intensity) of the average AP metrics are shown for each compound at indicated concentrations (see also Supplemental Table 2). Ranolazine has a monophasic increase of rise time and APD 30/50/80 , whereas APD MxR and APD tri are biphasic, suggesting nuanced dose-dependent effects (Fig. 7A,C). These results align well with the known clinical effects at low doses and the drug label warnings of fatal arrhythmias (which appear at high doses). In contrast, BPA shows a biphasic response in APD 30/50/80 over the concentration range tested (Fig. 7B). Toward understanding the role of hERG channel blockade in these contrasting results, we applied computer simulations to our system. HERG channel blockade has been shown to increase triangulation of AP (APD tri ) since I Kr peaks during phase 3 repolarization. Our experimental data and computer simulations in Supplemental Fig. 6 show that E4031 (or I Kr blockade in computer modeling) increases the slope of APD tri vs. APD MxR while I to block by 4-AP or augmentation of I Ca by BayK8644 does not alter the slope of APD tri vs. APD MxR , supporting that APD tri is an important indicator of hERG channel blockade. Ranolazine increases the slope of APD tri vs. APD MxR associated with significant APD prolongation at 100 µM (Fig. 7C), indicating that hERG blockade becomes significant at this concentration. However, BPA shows biphasic effects on APD tri vs. APD MxR , increasing the slope in the range of 1-100 nM but not at 1000 nM. This suggests that BPA potentially blocks hERG starting at very low concentration (1-10 nM) but later alters other ion channels to mask hERG channel blockade and overall these integrated ion channel changes shorten APD. The results presented here support the notion a contention that our 3D human cardiac microtissue model and comprehensive data analysis provide important information regarding cardiotoxicity of unknown compounds. www.nature.com/scientificreports/

Discussion
In this study, we address a critical need in risk assessment of compounds for pro-arrhythmic cardiotoxicity with a 3D human cardiac microtissue model containing purified hiPSC-derived cardiomyocytes and human cardiac fibroblasts. We show robust quantification of eight distinct parameters derived from the action potential waveform and qualify our NAM by confirming physiological responses and differentiation between the high-risk compound E4031 and low-risk drug ranolazine. Importantly, this differentiation enables reduction of false-positive results, as this 3D microtissue model integrates all ion channels and signals in the human cardiomyocytes. Cardiotoxicity can arise from a variety of aberrant signaling so predictive platforms must be suited for situations in which the mechanism of action is not well defined 53 . Testing of a compound with poorly described or unknown cardiac effects such as BPA reveals nuanced responses (Fig. 6), creating a drug response signature, which enables arrhythmic risk and safety evaluation that yields insights into mechanisms of action and is relevant to physiological metrics of health. Successful pro-arrhythmic drug screening should be able to detect changes in action potential shapes associated with increased triggered activity (EADs and DADs). A major challenge in developing a predictive risk assessment platform is that multiple factors can contribute to initiation of cardiac arrhythmias [67][68][69][70][71] , necessitating that the platform be able to capture multiple pro-arrhythmic changes in AP shape. Here, we report the ability to detect altered excitability (% of captured beats), AP initiation and upstroke, beat-to-beat APD stability, and triangulation of AP shape as well as more standard metrics of prolongation of APD and EAD activity that indicate very high risk for arrhythmia. Our platform goes beyond what is reported by other platforms (such as spontaneous beat rate) that focus on cellular toxicity using live/dead, mitochondrial, and ER imaging 38 , RNAseq for gene expression 72 , or contractile amplitude and kinetics from microtissue edge detection 73 . rather than arrhythmia. Due to our high temporal resolution and signal-to-noise ratio, just a single voltage recording without signal averaging can be used to analyze APD metrics and even detect beat-to-beat instability in APD. These parameters can be compared on the same microtissue before and after compound exposure (20 min) using a paired t-test, which increases statistical power and reduces the number of experiments. Other platforms report the calcium transient, using it as a surrogate metric for the voltage signal, but unlike our approach, this surrogate does not provide mechanistic clues underlying arrhythmia induction. The AP metrics defined in this study are more numerous and less variable than those reported from 2D monolayer assays (using voltage/calcium dyes or microelectrode arrays) 47,49,50 , 3D spheroid assays 45 , or engineered heart tissues (EHTs) 43,44 . These metrics provide integrated data on how the compound impacts electrical activation through changing the AP or Ca handling (Supplemental Table 3). The metrics in our data analysis were selected to account for potential toxicity effects on major ion www.nature.com/scientificreports/ channels that form cardiac action potentials including sodium channels responsible for AP upstroke (stim delay and rise time), transient outward potassium channels responsible for phase 1 repolarization (APD 30 ), L-type Ca channels responsible for maintaining AP plateau and overall duration (APD 30 , APD 50 , APD 80 ), delayed-rectifier potassium currents responsible for phase 3 repolarization (APD MxR , APD tri ), and inward rectifier underlying final repolarization (APD MxR , APD 80 ) and maintaining resting membrane potential that determines excitability together with sodium channels (excitability). The 3D microenvironment of our microtissues has several advantages to provide quantitative data for predicting pro-arrhythmic toxicity of chemical compounds. The AP recordings here are the averaged behavior from 12,000-24,000 cells that can mask individual cell-to-cell variation to increase robustness of the response. Triggered activity is thought to be difficult to initiate in the 3D environment due to electrotonic dissipation of depolarizing currents during EADs, known as 'source-sink mismatch' 74 . Therefore, in order to initiate triggered activity, a large number of cells should independently elicit EADs/DADs simultaneously. Our 3D microtissues have a sufficient number of CMs and do generate EADs (Fig. 4 and Supplemental Fig. 3), suggesting that this model replicates the true risks of EAD formation with its multi-cellular 3D environment as is present at the organ level in the human heart. In a direct comparison of the impact of microenvironment on the predictive capacity of an in vitro model using hESC-CMs, Archer et al. showed with a ROC analysis of drug responses to structural cardiotoxins that a 3D microtissue platform has increased specificity versus the same measurements made in 2D monolayers 38 .
Our 3D human cardiotoxicity model has several advantages that enable customization and robust proarrhythmic risk evaluation. These advantages are: (1) experimentally determined cell ratios and cell phenotypes (Fig. 1B); (2) self-assembly with organotypic heterocellular interspersion (Fig. 1D,E) while eliminating unnatural substrate or matrix stiffness; (3) robust 3D electrical coupling; (4) rapid optical measurements to track AP shape for highly accurate quantitative metric extraction (Fig. 1F-H); (5) in addition to spontaneous beats, paced action potentials can be analyzed to increase predictability of APD changes under drug; and (6) generation and functional analysis of a large number of individual but consistent microtissues to provide greater throughput and high statistical power in analyses (Figs. 2C, 3D). We have carefully considered the biology of the cardiac cells and employed culturing techniques like moderately longer culture periods and the broadly adopted methods for metabolic selection to purify and mature the hiPSC-CMs 75 , which demonstrate a more mature ventricular physiology in their electrophysiology compared to unpurified populations. While more mature functional responses are observed in our study, future work will require analysis of sarcomeric protein and ion channel expression to characterize this maturation. Further, there is evidence that the presence of non-myocytes from hiPSCs or other sources promote more mature electromechanical function 76 . Incorporation of human cardiac fibroblasts (hCFs) enables heterotypic cell-cell interactions characteristic of the intact myocardium 36,77,78 . While 5% primary adult normal hCF content is used for stable, healthy cardiac electromechanical function based on our previous 79 . and current work, the platform allows for alterations in cell composition (such as other cell types, including endothelial cells) and ratios to mimic physiological and pathophysiological conditions that may be important for toxicity evaluation. While our platform can produce hundreds of microtissues in standard cell culture plates, the low variability beat-to-beat, microtissue-to-microtissue, mold-to-mold, and notably hiPSC-CM differentiation batch-to-batch (Figs. 2C, 3D), that has not been reported in other established models, reduces the required sample size to tens of microtissues per condition, effectively increasing experimental throughput.
Screening of compounds with diverse and known mechanisms of action is widely accepted as critical practice for assessing the validity of in vitro models, but developing metrics to understand biological mechanisms of action is necessary to move towards screening compounds with unknown mechanisms of action 53 . Many established cardiac spheroids and elongated microtissue models have not tried to distinguish responses to high-risk and low-risk drugs 34,45 . Irregular beating patterns in response E-4031 has been shown, but from tissue contraction, not voltage changes to show EADs [42][43][44] . Some established MEA models have been able to differentiate between the high-risk compound E-4031 and the low-risk drug ranolazine, and our results align well with previous MEA studies that have shown prolonged APD and prevalent EADs with E-4031 80,81 . Ranolazine is consistently shown to modestly increase APD, especially with increasing concentrations, but the presence of EADs is not consistently seen and may be dependent on the hiPSC cell line utilized and the commercial source of MEA platform 32,82 . Computational modeling also predicts increasing APD with ranolazine and does not suggest EADs 83 , and preclinical and clinical data suggest that ranolazine suppresses EADs and DADs 65 . While prolonged APs and EADs are important metrics in arrhythmia risk assessment, previous work in a 2D MEA platform of hiPSC-CMs alone has shown an inability to distinguish compounds that elicit TdP from those that are benign 84 . Because EADs at the cellular scale can be dampened at the tissue scale due to source-sink mismatch with coupled myocytes that are not prone to EADs 85 , 3D platforms may be more predictive of arrhythmia risk. In the clinic, ranolazine rarely causes TdP and is frequently used as anti-arrhythmic treatment despite label warnings of QT prolongation 57,64,86 , suggesting that cardiotoxicity screening must differentiate between high risk compounds and low risk compounds (that could pass safety standards and provide clinical benefit) as we show for E4031 and ranolazine, respectively (Fig. 4). Our platform effectively captures expected responses to standard metrics like APD, yet is less prone to EAD (for example, with ranolazine), suggesting that our 3D microtissue platform and expanded metrics may capture more nuanced responses to more accurately assess risks, predict acceptable exposure levels, and identify multiple modes of toxicity. With ranolazine, we show nuanced dose-dependent effects on the AP quantified by our eight metrics (Fig. 5). This altered excitation reflects ranolazine's blockage of multiple ion currents including late I Na , I Ca , I Na-Ca , I Ks , and I Kr 56,87 . To clearly and quantitatively evaluate arrhythmic risk, we further analyzed dose-dependent changes in triangulation of the AP, APD tri (defined by APD MxR -APD 50 ), which has been implicated as a pro-arrhythmic predictor [88][89][90] . Previous studies indicate that an increase in APD tri may originate from blockade of either I Kr or I Ks , resulting in delayed phase 3 repolarization 88,91 . Our data also showed that specific I Kr blockade by E4031 increases APD tri but that I to blockade by 4-AP or I Ca activation by BayK8644 had no effect on www.nature.com/scientificreports/ APD tri , which we confirm with computer modeling (Supplemental Fig. 6). The trend of increasing APD tri with a similarly increasing APD would suggest risky I Kr block and can be visualized by plotting APD tri versus APD MxR . Previous computer modeling studies predicted that I Kr blockade increases APD tri 92,93 . I Kr has unique activation and inactivation kinetics that cause a large rush of outward current, which in turn causes rapid repolarization during phase 3 repolarization. Blocking I Kr , therefore, delays phase 3 repolarization and increases APD tri . In our experiments, we used APD tri versus APD MxR , which is a normalized APD tri , and found that APD tri is not always associated with APD prolongation experimentally (e.g., with 4-AP an I to blocker or BayK an I Ca activator). I Kr block (by E4031) specifically increases APD tri in line with predictions by computer simulation studies. The slope of the regression line thus increases when I Kr is blocked and is a clear indicator for dose-dependent arrhythmic toxicity. At lower concentration, ranolazine did not affect slope but at 100 μM it significantly increased (Fig. 7C, bar plot), suggesting that ranolazine blocks I Kr , which is in agreement with previous voltage clamp studies that found IC 50 is 11.5 μM for I Kr 87 . However, the less well-known compound BPA showed a biphasic effect on slope of the APD tri vs. APD MxR relationship, peaking at 10 nM and decreasing at higher concentrations, while APD MxR decreases (Fig. 7D). These data suggest that I Kr may be affected at BPA concentrations starting at very low concentration (10 nM) and yet other ion channels and currents may compensate to counteract this effect at higher concentrations. Ultimately, this integrated AP response in human 3D cardiac microtissues with hiPSC-CMs is a defining feature of our model and a significant advantage to be able to predict pro-arrhythmic cardiotoxicity with high sensitivity that will be useful for targeted pharmaceutical development earlier in the drug development pipeline or predictive IVIVE to define safe chemical exposure levels.
It was our aim with this study and in our human cardiotoxicity tissue-engineered model to provide a fit-forpurpose NAM for advancing toxicity risk assessment for pharmaceutical development and safety evaluation of chemicals and environmental toxicants while reducing use of animals, improving mechanistic understanding of cardiotoxicity, and facilitating broad adoption of this simple human in vitro platform. Our platform showed expected changes in arrhythmia metrics to known toxicants and revealed alterations in response to the endocrine disrupting chemical BPA with high sensitivity that suggest cardiac risk and point toward disrupted ion channels as mechanisms of action that may not have been detected in other platforms with more surrogate metrics. Further, automated image analysis and computational approaches including machine learning that have mainly been applied to 2D systems and focus on different outcome metrics such as Ca transient and contractile motion tracking [94][95][96] . can be applied for AP analyses in our platform to increase its predictive power. Cardiotoxicity risk likely depends on multiple factors including sex, race, and disease status 97 , and toxic effects may be observed to a lesser degree in healthy human subjects compared to patients with comorbidities, defined as 'hidden toxicity' 18 . These diverse risk factors can be incorporated into our 3D microtissue risk assessment platform through the use of different hiPSC lines; thus exciting opportunities exist to better understand many fundamental mechanisms of toxicity. The predictive value of this platform has the potential to accelerate broad adoption and advances in the regulatory landscape for cardiotoxicity evaluation.
Fabrication of hydrogels and 3D culture. Scaffold-free 3D spherical microtissues were generated using non-adhesive agarose gels with cylindrical microwells with hemispherical bottoms to guide self-assembly (Fig. 1B). Sterilized 2% (wt/vol) agarose was pipetted into molds designed for 24-well plates with 800-μm-diameter rounded pegs (3D Petri Dish, MicroTissues). After being cooled to room temperature (~ 5 min), the agarose gels were separated from the molds and transferred to single wells of 24-well plates. For equilibration, 1 mL RPMI + B27 + 1% P/S medium was added to each well. Hydrogels were equilibrated for at least 1 h at 37 °C in a humidified incubator with 5% CO 2 . Molds were transferred to 6-well plates for electrical stimulation, and hiPSC-CMs with or without additional 5% hCFs in suspension were added to the center of the hydrogel seeding chamber (420-840 K cells/mold in 35 recesses) and allowed to settle into the recesses for 30 min. Medium (RPMI + B27 with 1% P/S and 10% FBS with 5 μM rock inhibitor (Y27632)) was then added to each well (5 mL). Medium was changed to RPMI + B27 with 1% P/S and 10% FBS at 24-48 h, and cells were cultured for 6-8 days www.nature.com/scientificreports/ with media changes every other day. During the 3D culture period, the self-assembled microtissues were field stimulated with a 1 Hz, 10.0 V, 4.0 ms duration bipolar pulse train with an Ionoptix C-Pace EP. To study hCF interspersion within tissues, hCF were incubated for 30 min in CellTracker Orange working solution and hiPSC-CM were incubated for 30 min in CellTracker Green working solution before seeding in 3D hydrogels. Spheroids were recovered from hydromolds and transferred to glass bottom dishes. Images were taken with an Olympus FV3000 confocal microscope and processed using ImageJ.
Microtissue size analysis. Stitched 4 × phase-contrast images of whole 24-well microtissue hydrogels were captured with a Nikon TE2000-U and a black and white/color digital camera (MicroVideo Instruments, Avon, MA). NIS Elements software was used for image acquisition and analysis. Image thresholding and particle size analysis was performed to determine the top view cross-sectional area of individual microtissues across each mold.
3D tissue sections and immunohistochemistry. Microtissues in 24-well hydrogels were fixed with 4% (vol/vol) paraformaldehyde (Electron Microscopy Sciences, Hatfield, PA) and 8% (wt/vol) sucrose in PBS overnight at room temperature. Molds were then rinsed twice with PBS and fully equilibrated (as indicated by their sinking, usually over 12 h) with 15% and then 30% (wt/vol) sucrose in PBS. Whole agarose gels containing microtissues were removed from sucrose solution, blotted dry, and embedded in Tissue-Tek CRYO-OCT Compound (Sakura). Blocks were stored at − 80 °C, sectioned on a Leica CM3050 cryostat microtome (Leica Biosystems, Buffalo Grove, IL) into 10-μm-thick sections, and placed on Superfrost Plus slides. After being air dried for 15 min, sections were post-fixed in 4% paraformaldehyde in PBS. For immunofluorescent staining at room temperature, frozen sections were rinsed 3 times for 5 min with 1 × PBS wash buffer. Non-specific binding was blocked with 1.5% goat serum for 1 h, followed by 24 h incubation in primary and followed by a 1 h incubation in secondary antibodies diluted in 1.5% goat serum. Primary antibodies were directed against cardiac troponin I (cTnI, 1:100, Abcam ab47003) and vimentin (1:100, Sigma V6630), and secondary antibodies were conjugated to Alexa Fluor 488 or Alexa Fluor 594 (1:200, Invitrogen). Coverslips were mounted with Vectashield mounting medium with DAPI. Images were taken with an Olympus FV3000 Confocal Microscope and processed using ImageJ.
Optical mapping and automated action potential analysis. Microtissues were loaded with voltagesensitive di-4-ANEPPS (5 μM for 10 min at 35 °C) for measurements of membrane potential (V m ). Fluorescence images were acquired at 979 frames/s using a Photometrics Evolve + 128 EMCCD camera (2 × 2 binning to 64 × 64 pixels, 18.7 × 18.7-μm 2 resolution, 1.2 × 1.2-mm 2 field of view) and an Olympus MXV10 macroview optical system 101 . Typically, four microtissues were recorded simultaneously per scan at this magnification. A step-by-step illustration of automated data analysis is available in Supplemental Fig. 1. Briefly, the pixels with APs were identified from Fast Fourier transformation (FFT) of fluorescence signals. After appropriate thresholding and image segmentation, the region of each microtissue was grouped and the fluorescence signals from the pixels in the same microtissue were average and used for AP analysis (Supplemental Fig. 1).
Validation and screening of toxicants for arrhythmogenic risk. A single mold of microtissues was mounted on a temperature-controlled chamber (Dual Automatic Temperature Controller TC-344B, Warner Instrument) to maintain 35 ± 1 °C and bathed with Tyrode's solution containing (in mM) 140 NaCl, 5.1 KCl, 1 MgCl 2 , 1 CaCl 2 , 0.33 NaH 2 PO 4 , 5 HEPES, and 7.5 glucose. Microtissues were stimulated with a platinum field stimulation electrode (Supplemental Fig. 7, Myopacer EP field stimulator, IonOptix, Milton, MA). The test compounds including E4031, 4-AP, BayK8644, ISO were purchased from Sigma Aldrich and dissolved in 100% DMSO to prepare 0.01-0.5 mM stock solutions. BPA was dissolved in 10% ethanol stock solution and diluted in Tyrode solution to the final concentration. Microtissues were exposed to vehicle (DMSO or ethanol) and the indicated concentrations of test compounds for 20 min, and action potentials were measured as described above.
Statistical analyses. The data from optical mapping are expressed as mean ± SD for n microtissues unless otherwise indicated. Statistical analyses were performed using Student's two tailed paired and unpaired t test. P values of 0.05 were considered statistically significant. Normality test was done using Kolmogorov-Smirnov test. The test for the equality of regression coefficients was done using Z statistics of two slopes and SE as described 102 www.nature.com/scientificreports/