Bayesian Approach Enabled Objective Comparison of Multiple Human iPSC-derived Cardiomyocytes’ Proarrhythmia Sensitivities

The one-size-fits-all approach has been the mainstream in medicine, and the well-defined standards support the development of safe and effective therapies for many years. Advancing technologies, however, enabled precision medicine to treat a targeted patient population (e.g., HER2+ cancer). In safety pharmacology, computational population modeling has been successfully applied in virtual clinical trials to predict drug-induced proarrhythmia risks against a wide range of pseudo cohorts. In the meantime, population modeling in safety pharmacology experiments has been challenging. Here, we used five commercially available human iPSC-derived cardiomyocytes growing in 384-well plates and analyzed the effects of ten potential proarrhythmic compounds with four concentrations on their calcium transients (CaTs). All the cell lines exhibited an expected elongation or shortening of calcium transient duration with various degrees. Depending on compounds inhibiting several ion channels, such as hERG, peak and late sodium and L-type calcium or IKs channels, some of the cell lines exhibited irregular, discontinuous beating that was not predicted by computational simulations. To analyze the shapes of CaTs and irregularities of beat patterns comprehensively, we defined six parameters to characterize compound-induced CaT waveform changes, successfully visualizing the similarities and differences in compound-induced proarrhythmic sensitivities of different cell lines. We applied Bayesian statistics to predict sample populations based on experimental data to overcome the limited number of experimental replicates in high-throughput assays. This process facilitated the principal component analysis to classify compound-induced sensitivities of cell lines objectively. Finally, the association of sensitivities in compound-induced changes between phenotypic parameters and ion channel inhibitions measured using patch clamp recording was analyzed. Successful ranking of compound-induced sensitivity of cell lines was validated by visual inspection of raw data. Author Summary Cardiac safety is one of the most stringent regulatory risk monitoring during drug development. Regulatory agencies, including the Food and Drug Administration (FDA), require drug developers to conduct thorough preclinical and clinical studies to evaluate drug-induced proarrhythmia risks. The CiPA (Comprehensive in vitro Proarrhythmia Assay) initiative led by the FDA validated the applications of human cardiomyocytes derived from induced pluripotent stem cells in proarrhythmia risk assessments. The assay can be cost effective, and use of high-throughput approaches enables scale-up analysis. However, limited diversity in cell lines to recapitulate heterogeneity of human patients are still lacking in the cardiac safety analysis. We applied various computational tools to maximize capacity for analyzing diverse response of commercially available five cell lines against reference chemical compounds. Limited number of experimental data made it difficult to predict similarities and differences in drug responses among different cell lines. By applying Bayesian approach, our analysis could intuitively grasp the probabilities of expecting different responses from different cell lines.

The one-size-fits-all approach has been the mainstream in medicine, and the well-2 defined standards support the development of safe and effective therapies for many 3 years. Advancing technologies, however, enabled precision medicine to treat a 4 targeted patient population (e.g., HER2+ cancer). In safety pharmacology, 5 computational population modeling has been successfully applied in virtual clinical 6 trials to predict drug-induced proarrhythmia risks against a wide range of pseudo 7 cohorts. In the meantime, population modeling in safety pharmacology experiments 8 has been challenging. Here, we used five commercially available human iPSC-derived 9 cardiomyocytes growing in 384-well plates and analyzed the effects of ten potential 10 proarrhythmic compounds with four concentrations on their calcium transients 11 (CaTs). All the cell lines exhibited an expected elongation or shortening of calcium 12 transient duration with various degrees. Depending on compounds inhibiting several 13 ion channels, such as hERG, peak and late sodium and L-type calcium or IKs 14 channels, some of the cell lines exhibited irregular, discontinuous beating that  Kit (Molecular Devices LLC, CA, USA) by incubation for 2 hour in culture medium. 5 Simultaneous fluorescent signals from the entire plate were recorded on an FDSS-6 µCELL (Hamamatsu Photonics K.K., Shizuoka, Japan) fluorescent plate reader for 5 7 minutes. A baseline signal was recorded before compound treatment. The cell plate 8 was treated simultaneously with 4 doses of 10 compounds with each dose having 5 9 replicate samples. After 10 minutes in the above incubator, a compound-treated 10 signal was recorded. 11 Test compounds: 12 Ten compounds and their concentrations used for the analysis is listed in Table 1. 13

Waveform analysis: 14
After automatically finding all the CaT cycles using iVSurfer (InvivoSciences, Inc, 15 Madison, WI, USA). The six parameters including peak-width duration 80 (P80), concentrations of ten drugs with known potential proarrhythmic risks (Table 1) were 26 applied to the hiPSCMs growing in 384-well plates using protocols and media 27 provided by each vendor. Each well received a single dose of a drug, and CaT traces 28 before and after drug treatments were compared. As we and others [1, 13] utilized 29 durations of CaTs as a key metric to gauge drug-induced proarrhythmia sensitivity of 30 hiPSCMs, dose-dependent changes in CaT duration average of all cell lines were 31 plotted ( Figure S1). Most of the hERG inhibitors increased CaT duration 32 significantly, as expected. However, by inspecting individual CaT traces, we found 1 that some cell lines stopped beating or seldom but consistently exhibited irregular 2 beats at some dose levels ( Figure 1G). Therefore, simply analyzing only CaT duration 3 or beat rates could miss those proarrhythmic events. Consequently, the diligent 4 observations of raw data traces of many cell lines pointed to the importance of 5 multidimensional data analysis. 6 Six-parameter analysis 7 Multiple parameters were used to define differences in profiles of action potential 8 and CaT previously [14]. We applied a similar approach to parameterized profiles of 9 CaTs and their irregular beating to classify different trace patterns. For analyzing 10 >2,000 traces of each with ~10,000 time-points, we utilized a data analysis software 11 specialized in periodic data (see method), which automatically detected regular and 12 irregular cardiac beats ( Figure 1G Figure S2). Finally, the changes in the six-dimensional parameters 22 of before and after drug treatment (after drug/before drug-1) were plotted radially to 23 visualize compound-induced response patterns ( Figure 1J). 24

25
One of the major classes of proarrhythmia compounds is an hERG inhibitor. An 26 experimental class III anti-arrhythmia drug, E4031, inhibiting the hERG channel 27 with high specificity, was utilized to visualize drug-induced phenotypic response 28 patterns using the multiparameter approach. E4031 significantly prolonged the 29 duration of CaT and induced EADs in all cell lines with various ranges (2-5 fold 30 increase, Figure 2 A-F). Miracell was the most sensitive to E4031 and stopped 31 beating at its highest concentration ( Figure 2E). Asclestem and Axol also elongated 32 their CaT duration, but the PWD80/20 ratios did not increase as much as CDI, which 33 was confirmed by their ensemble average traces (Figure 2A, B, C). We also noticed 1 that E4031 not only elongated Myoride's CaT duration but also increased cycle length 2 SD ( Figure 2G). Among all, Myoridge was the least sensitive to E4031; PWD80 did 3 not change much until adding the highest dose of E4031 ( Figure 2F, K). Nonetheless, 4 the radial plots of multiple parameters visually captured the different response 5 patterns of five cell lines. Importantly, the response patterns of multiparameter 6 radial plots correctly represented the drug-induced changes in beat patterns, 7 providing an intuitive sense of how five cell lines responded to E4031. Ca channel inhibition and multi-channel inhibition 28 The radial plot visualization revealed cell line-specific sensitivity features unique to 29 the L-type calcium channel inhibitors with or without inhibiting other channels. We 30 used Verapamil & Terfenadine inhibiting L-type calcium (CaV1) and hERG channels 31 and Amiodarone inhibiting fast and late sodium channels and Ito current in addition 32 to hERG and CaV1 channels. All these drugs inhibited the hERG channel to a 1 different extent; however, almost no prolongation of CaT duration was recorded 2 ( Figure 4A-C). Asclestem and Myoridge lines first increased cycle lengths and their 3 SDs (i.e., irregular slow beats) and then stopped beating at higher concentrations. 4 Finally, when Chrom293B (Chrom) inhibiting IKs, Ito, and IKur (ultra-rapid 5 outward current) with marginal IKr but without affecting sodium and calcium 6 channels was applied, all cell lines showed slight prolongation of CaT duration and 7 stopped beating at the highest concentration in some cell lines ( Figure 4D). 8 and IKr, but Axol and Myoridge formed a cluster but were more responsive than the 30 other cell lines ( Figure S12, 3E). So, Axol and Myoridge could be a valuable tool for 31 analyzing the mechanism by which Moxifloxacin induced arrhythmia [28]. CDI and 32

Classification of cell lines
Miracell were the most sensitive to hERG inhibitors. A specific L-type calcium 33 channel inhibitor, Verapamil, lowered CaT amplitude of CDI and Axol and eventually 34 stopped their beating but increased cycle length without prolonging CaT duration in 1 Asclestem, Miracell, and Myoridge. A similar trend was seen in Terfenadine and 2 Amiodarone which inhibit IKr and sodium channels in addition to L-type calcium 3 channels. It is interesting to note that Verapamil and Amiodarone reduce heart rate 4 in patients with atrial fibrillation [29,30] [31], concurring slowing down beat-rate of 5 Asclestem, Miracell, and Myoridge cells 6 Terfenadine, was removed from the market due to its proarrhythmic risk for long 7 QT-related Torsades de Pointes (TdPs) [32]. The study data using isolated rabbit 8 heart and human atrial myocyte showed that Terfenadine induced non-TdP like 9 ventricular tachycardia/ventricular fibrillation due to a slowed conduction via        (A) Diagram showing CCA between phenotypic parameters and ion channel inhibition data. The CCA scatter plot shows the correlation between the first pair of canonical covariates: CCX1 and CCY1. CCX1 is the first canonical covariate for the phenotypic data and CCY1 is the first canonical covariate for the ion channel inhibition data. Each dose is plotted by a separate color: dose1 (blue), dose2 (orange), dose3 (green), dose4 (red). The correlation ecoefficiencies between the first pair of canonical covariates for each drug and each cell line were compared.  The radioplt showed no changes in P80 but increased cycle length (CL) and their standard deviation CSD). Figure S3. 2D plot of raw data after PCA.
Five sample data after PCA were plotted on the 1st and 2nd component plane. Data treated by the same compounds seem to form culsters but their distribution were not clearly visible. Figure S4. PCA analysis data in three-dimensional space. The 1st, 2nd, and 3rd Component represents each axis of the first three principal components. After Bayesian bootstrapping, DMSO control (1,000 black spheres) and chromanol 293B (Chrom, 1,000 purple spheres) data formed clusters in the 3D space. Standard deviation (SD) of each component axis of control and compound-treated data defined ellipsoids for control and chromanol 293B, enclosing data points. Sizes of ellipsoids represented data distribution with different SD sizes.