All the in silico drug trials presented in this study were performed in a population of 1,213 human ventricular control models, built using the O'Hara-Rudy dynamic (ORd) model (O'Hara et al., 2011) as baseline and the methodology described by Britton et al. (2013) and further discussed by Muszkiewicz et al. (2016). The ORd human ventricular AP model was chosen for this study because of: (i) the large number of human ventricular experimental data obtained from more than 140 hearts used in its construction and evaluation; (ii) its ability to reproduce and probe pro-arrhythmic mechanisms, including repolarization abnormalities and APD alternans, as shown in multiple studies and reviewed by Britton et al. (2017); (iii) its choice within the CiPA initiative (Sager et al., 2014; Colatsky et al., 2016).

Ionic conductances were sampled in the [0–200]% range of the baseline model values, to include both healthy and potentially abnormal ionic current profiles (with low/high ion channel densities corresponding to loss/gain-of-function of specific ionic channels due to e.g., genetic mutations), but still with a healthy-looking AP. These are important as they have been implicated in increased pro-arrhythmic risk (Sanguinetti and Tristani-Firouzi, 2006; Itoh et al., 2016; Wang et al., 2016). Nine ionic conductances were considered: fast and late Na⁺ current (G_Na and G_NaL respectively), transient outward K⁺ current (G_to), rapid and slow delayed rectifier K⁺ current (G_Kr and G_Ks), inward rectified K⁺ current (G_K1), Na⁺-Ca²⁺ exchanger (G_NCX), Na⁺-K⁺ pump (G_NaK), and the L-type Ca²⁺ current (G_CaL).

Only the AP models exhibiting a phenotype in agreement with human experimental data from undiseased hearts (Britton et al., 2017) were selected for the control population, which consists of 1,213 models. A more detailed description of the control population used in this study is included in the Supplementary Material, together with the experimental AP biomarker ranges used for the calibration process (Table S1) and the scaling factors of the ionic conductances for the 1,213 models (Table S2).

A total of 62 reference compounds were considered in this study. The list includes antiarrhythmic drugs in Classes I, III, and IV, as well as other drugs used for different purposes but with known effects on cardiac ion channels. Most of these drugs have a multichannel action, which makes prediction and interpretation of cardiotoxicity challenging. Drugs were selected to include all the ones in Kramer et al. (2013), as well as 15 compounds widely used as reference compounds, which were characterized in more depth both in simulations and experiments and are listed in Table 1.

Each drug was assigned to a TdP risk category, based on the classification by CredibleMeds® (Woosley and Romer, 1999), available on www.crediblemeds.org (as of July 2017): 1 (high risk), the drug prolongs the QT interval and is clearly associated with a known TdP risk, even when taken as recommended; 2 (possible risk), the drug prolongs the QT interval, but there is a lack of evidence of TdP risk when taken as recommended; 3 (conditional risk), the drug is associated with TdP but only under certain circumstances, e.g., excessive dose or interaction with other drugs; NC (not classified), the drug was reviewed by CredibleMeds® but the evidence available was not enough to assign it to any of the previous categories, and therefore no action was taken. Of the 62 compounds, 24 are classified as high risk and 13 as potential/conditional risk, for a total of 37 drugs associated with TdP risk. Verapamil (classified as NC) and the remaining 24 compounds (not listed) are considered as TdP category 0 (no TdP risk) for the purpose of this study.

Drug effects were simulated using a simple pore-block model consistent with data available for drug/ion channel interactions, consisting of IC₅₀ and Hill coefficient (h) for each drug/ion channel. Up to 7 ion channels were considered for this study: fast Na⁺ current (I_Na), rapid/slow delayed rectified K⁺ current (I_Kr/I_Ks), transient outward K⁺ current (I_to), L-type Ca²⁺ current (I_CaL), inward rectifier K⁺ current (I_K1), and late Na⁺ current (I_NaL). The experimental IC₅₀ and h used for the drug assays were collected mainly from three different sources: (i) our internal database, measured with either manual or automated patch-clamp techniques (when the IC₅₀ concentration was not reached in the experiments, an estimate was computed from the percentage of block at the maximum tested concentration, with h equal to 1); (ii) (Kramer et al., 2013), data acquired with automated patch-clamp; (iii) (Crumb et al., 2016), data acquired with manual patch-clamp.

For the compounds that were included in more than one of these datasets, multiple inhibitory profiles were considered to investigate the impact of variability in drug characterization. Each IC₅₀ and h set was simulated separately, resulting in 87 different drug trials: each trial is referred to with the name of the compound together with a roman numeral, to differentiate multiple entries (e.g., Bepridil I, Bepridil II, Bepridil III).

Multiple concentrations were investigated for each compound, chosen to match those used in the experimental drug assays, as well as to explore different multiples of the maximal effective free therapeutic concentration (EFTPC_max), up to 100-fold. The EFTPC_max values were taken from literature, mainly from Kramer et al. (2013) or Crumb et al. (2016). When multiple values were found for the same compound, the higher one was considered for simulations.

The full list of compounds, together with the IC₅₀, Hill coefficient and the EFTPC_max used for in silico drug trials is provided in Table S3.

All the simulations presented in this study were conducted using Virtual Assay (v.1.3.640 2014 Oxford University Innovation Ltd. Oxford, UK), a user-friendly C⁺⁺ based software package with a graphical user interface for in silico drug assays, to facilitate its use by non-experts in computational modeling, and available upon request. Virtual Assay uses the ordinary differential equation solver CVODE, part of the open-source Sundials suite (Hindmarsh et al., 2005; Serban and Hindmarsh, 2005), implementing time adaptive Backward Differentiation Formulas with relative and absolute tolerance equal to 1e-5 and 1e-7, respectively. Our results could therefore be replicated using other software products, as Matlab (Mathworks Inc. Natwick, MA, USA) or Chaste (Pitt-Francis et al., 2008). As an example, a comparison of simulations obtained with Virtual Assay and with the Matlab solver ode15s (Shampine and Reichelt, 1997) is shown in the Supplemental Material, Figure S2. Simulation results were analyzed both in Virtual Assay and in Matlab.

Following drug application, all models were stimulated at 1 Hz for 150 beats, and the last AP trace in each simulation was analyzed. AP biomarkers were extracted, including: AP duration at 40, 50, and 90% of repolarization (APD₄₀, APD₅₀, APD₉₀); APD₉₀ dispersion, defined as the difference between the maximum and minimum value of APD₉₀ in the population (ΔAPD₉₀), AP triangulation, defined as the difference between APD₉₀ and APD₄₀ (Tri₉₀₋₄₀); maximum upstroke velocity (dV/dt_MAX); peak voltage (V_peak); resting membrane potential (RMP), computed as in Britton et al. (2017). Drug-induced changes in those biomarkers are presented as percentage change in median with drug, compared to control (no drug).

All AP traces were automatically checked for repolarization and depolarization abnormalities (RA and DA, respectively). RA were defined as the presence of a positive derivative of the membrane potential (V_m) 150 ms after the AP peak (representative of early after-depolarizations, EADs), or when the membrane potential did not reach the resting condition following an AP upstroke (V_m > −40 mV) by the end of the beat. DA were defined as AP traces in which the upstroke phase was compromised, i.e., when the max upstroke V_m was lower than 0 mV, or when the time needed to reach 0 mV was longer than 100 ms.

Drugs were classified as risky when RA occurred in the population of models at different concentrations, based on: true positives (drug with reports of TdP risk classified as risky); true negatives (drug with no reports of TdP risk classified as safe), false positives (drug with no reports of TdP risk, classified as risky); false negatives (drugs with reports of TdP risk, classified as safe). The performances of the classification were evaluated based on: sensitivity, defined as the number of true positives divided by the sum of true positives and false negatives; specificity, defined the number of true negatives divided by the sum of the true negatives and false positives; accuracy, defined as the sum of true positives and true negatives divided by the total number of drugs. Classification results based on RA were compared against the ones obtained for APD prolongation at 10x EFTPCmax. APD₉₀ prolongation threshold to define risk was fixed to 6%, considering the correspondence between QTc and APD₉₀, and the current guidelines suggesting QTc prolongation >20 ms (which correspond to 5.7% for a normal QT of 350 ms) as a definite risk factor for TdP (Salvi et al., 2010). Results for the population of models were also compared against the same results obtained with the single ORd model.

A scoring system was developed by integrating RA occurrence at multiple concentrations. The fraction of models developing RA was multiplied by a factor inversely related to the drug concentration at which those RA occur (e.g., 1/100 for RA occurring at 100x EFTPC_max). Contributions from all the different concentrations were added together, and the total score was normalized, according to the following formula (where nRA_i is the number of models showing RA at the tested concentration i, w_i = EFTPC_max / i is the weight inversely related to the tested concentration i, and n_mod is the total number of models in the population).

The TdP score thus obtained varies between 0 and 1, where 0 corresponds to a drug with no RA, and 1 to a drug which develops 100% of RA at every concentration. By using the proposed score, RA are considered more severe when occurring at low concentrations and/or affecting a high fraction of the population of models.

In silico results were compared with properties computed from ECGs in rabbit wedge preparations and CT recordings from hiPS-CMs. Experimental data were acquired for the 15 compounds listed in Table 1 at multiple concentrations, as described below.

Recordings of ECGs from left ventricular rabbit wedge preparations have been previously described and partly published in Lu et al. (2016). The biomarkers extracted include QRS complex and QT interval duration, defined as the time from the onset of the QRS complex to the point at which the final downslope of the T wave crossed the isoelectric line.

CT recordings from hiPS-CMs (Cor 4U) were acquired as part of this study on pre-plated preparations from Axiogenesis (Cologne, Germany). Full method details are included in the Supplementary Material. Quantified biomarkers included CT beat rate (CTBR) and duration at 90% of the initial base value (CTD₉₀), known to be correlated with APD (Gauthier et al., 2012; Spencer et al., 2014), similarly to other studies (Lu et al., 2015; Zeng et al., 2016).

All experimental results are presented as median percentage changes with respect to the baseline. Drug-induced changes in experimental values need to be compared against the effect measured without drugs, i.e., with vehicle, defining cut-off values. In the rabbit wedge, the changes measured with the vehicle were always quite small: <5% for QT and <3% for QRS. On the other hand, in CT assays using hiPS-CMs, the lower and upper limit were 19% and 24% of the baseline, with 95% confidence interval (n = 222 vehicle controls). Therefore, only CTD₉₀ prolongations >25% were considered relevant for this assay. In silico, no biomarker changes are observed when the drug effect is not included, since the models are paced until steady state: therefore, the cut-off value is equal to 0%. Statistical analysis was performed with the Wilcoxon-Mann-Whitney Test by using R Project for Statistical Computing, and p < 0.05 was considered as significant. Very small p-values (e.g., p < 1e-6) were obtained in all simulation results due to the high number of models considered, and therefore we focus on the magnitude of the differences observed (White et al., 2014).

In silico drug trials for a total of 62 reference compounds were performed in the control population of 1,213 human ventricular AP control models, based on the ORd model (O'Hara et al., 2011) and constructed as described in Methods. Before drug application, all the models exhibit a healthy-looking AP phenotype, in agreement with human experimental recordings from un-diseased hearts (Figure S1A). When including drug effect for concentration up to 100-fold the EFTPC_max, each model responds in a different way, depending on its underlying ionic properties. We first evaluated drug-induced changes in the AP biomarkers.

Figure 1 shows APD₉₀ and dV/dt_MAX distributions from the in silico population for 5 compounds (Dofetilide I, Flecainide I, Nimodipine, Ranolazine I, and Verapamil II) at multiple concentrations. Extended results for additional compounds, including all AP biomarkers, are shown in Figures S3–S31.

Most drugs (all except Nimodipine and Nisoldipine) result in APD prolongation at 40, 50, and 90% of repolarization, as well as increased AP triangulation (Tri₉₀₋₄₀), mainly as a result of hERG channel block (e.g., Figure 1, left column and A–C in Figures S3–S31). For 30x EFTPC_max dose, Flecainide III, Bepridil I, and Dofetilide III showed the largest APD₉₀ prolongations (+180%, +175%, and +157%, respectively).

APD prolongation caused by BaCl₂ (mainly due to I_K1 block) is stronger in APD₉₀ compared to APD₄₀ and APD₅₀ (Figures S3A–C). As a secondary effect of I_K1 block, BaCl₂ leads to a decrease in RMP (Figure S3H), which also exhibits larger variability, while for all the other compounds it remains almost constant (H in Figures S3–S31).

Consistent with their expected mode of action, class I anti-arrhythmic drugs (Procainamide, Lidocaine, Mexiletine, Phenytoin, Flecainide) as well as other drugs affecting Na⁺ channels as secondary effect (e.g., Bepridil) show a strong decrease of upstroke velocity (e.g., Figures 1D,H), together with a decrease of V_peak (F, G in Figures S3–S31).

Verapamil II represents an interesting example of the combined block of I_Kr and I_CaL (Figure 1I). For low concentrations (0.01–0.1 μM) the Ca²⁺ block is predominant, and APD₉₀ is slightly decreased (−1 and −4% respectively), whereas I_Kr block compensates its effects for higher concentrations (>0.5 μM), resulting in a clear APD prolongation (e.g., +38% for 1 μM). Verapamil II also leads to slower AP upstroke for high concentrations (Figure 1J).

Results obtained using the baseline ORd model (Figure 1 and Figures S3–S31, black diamonds) are in overall agreement with the range of AP biomarkers in the population of human models. This is with the exception of cases in which the baseline ORd yields abnormal APs for high doses of certain drugs (e.g., Dofetilide I 0.1 and 0.2 μM, Figures 1A,B). In those cases, the human population of models still allows exploration of the full concentration range for each compound.

Drug action resulted in an increase in the phenotypic variability yielded by the human ventricular population, as illustrated in Figure 2 for Moxifloxacin III (Figure 2A) and Dofetilide I (Figure 2B), mainly blocking I_Kr, and for Flecainide I (Figure 2C), inhibiting both I_Kr and I_Na. Following drug application, some models in the in silico population display normal but prolonged APs (gray traces) while a fraction develop repolarization abnormalities (RA, pink traces). Due to its strong effect on I_Na, Flecainide I (Figure 2C) also cause an overall reduction of dV/dt_MAX, visible in the upstroke phase of the AP, and depolarization abnormalities (DA, green traces) in specific models.

Quantitative analysis of the underlying ionic mechanisms reveals consistency on the mechanisms underlying RA and DA across different drugs and concentrations. Models displaying RA are characterized by low G_Kr, and G_NaK, and high G_CaL and G_NCX, i.e., a reduced repolarization reserve (Figures 2D–F, pink vs. gray boxplots). Low G_Ks also plays a role when a larger fraction of the population displays RA (Figures 2D,E). Models displaying DA are characterized mainly by low G_Na/G_CaL/G_NaL, i.e. the net inward current in the initial phase of the AP is reduced (Figure 2F, green vs. gray boxplots).

We hypothesized that occurrence of RA following drug application in the human population would be predictive of in vivo TdP, given the potential mechanistic link between them (El-sherif et al., 1990; Dutta et al., 2016). In silico drug trial predictions were evaluated against clinical reports of TdP using the TdP risk categories provided by CredibleMeds® (Woosley and Romer, 1999) and further described in Methods. When multiple inhibitory profiles were simulated for the same compound, the worse scenario was considered, i.e., the higher occurrence of RA, and the larger APD₉₀ prolongation.

Figure 3 shows the classification results for the 49 compounds with either high or no TdP risk (Figure 3A), and for the full set of 62 compounds (Figure 3B) based on RA occurrence up to 100x EFTPC_max and APD₉₀ prolongation >6% at 10x EFTPC_max for both the population of models (top row), and the single ORd model (bottom row). Accuracy reached 96% for the classification of high vs. no TdP risk compounds using the RA-based classification with the in silico population of models, compared to 80% based on APD prolongation (Figure 3A). Using the single ORd model, the higher accuracy was 76% and was obtained using APD prolongation as biomarker.

When including also compounds with possible/conditional risk (Figure 3B), accuracy with the RA-based classification for the in silico population was 89%, compared to 81% based on APD prolongation (48 and 77% based on RA and APD prolongation, respectively, for the single ORd model).

In summary, all drugs with high TdP risk (category 1) were correctly identified as risky with the RA-based classification in the in silico population, and only 5 drugs with possible/conditional TdP risk (category 2 and 3) resulted as false negatives: Clozapine, Dasatinib, Paroxetine, Saquinavir, Voriconazole. It is worth noting that drugs with possible/conditional TdP risk lack evidence of pro-arrhythmic risk when taken as recommended: TdP reports are usually related to excessive dose or interaction with other drugs.

Overall, classification based on APD prolongation exhibited high sensitivity, but low specificity. Indeed, many compounds prolong APD without being associated with TdP risk, e.g., Verapamil, thus resulting in false positives. On the contrary, false positives are rare in the RA-based classification (Lidocaine and Mexiletine), and they only develop RA at the maximum tested concentration (100x EFTPC_max). Indeed, both Lidocaine and Mexiletine are Class Ib anti-arrhythmic drugs, which have been associated with cardiotoxicity in case of overdose (Denaro and Benowitz, 1989).

RA-based classification is dependent on the maximum tested concentration: a higher concentration is more likely to provoke RA in the in silico population, thus increasing sensitivity but possibly decreasing specificity, since even safe drugs might lead to RA at very high doses. We reported here the classification results obtained for concentrations up to 100x EFTPC_max, and a comparison between results for 30x and 100x EFTPC_max is included in the Supplementary Material (Figure S32).

Figure 4 shows all the tested compounds classified using the TdP score computed from the in silico drug trials using the fraction of models displaying RA at each tested concentration, as described in Methods. The TdP score varies from 0 to 1, and is higher when RA occur at low concentrations and/or affecting a high fraction of the population of models.

The distribution of compounds in the safe zone (TdP equal to 0, left side) and risky zone (TdP > 0, right side) reflects the classification results summarized in the confusion matrix with the higher (89%) accuracy (Figure 3). All safe compounds (no reported TdP risk, green dots) have a TdP score equal to 0, except Lidocaine and Mexiletine. All compounds with known risk of TdP (TdP risk category 1, red dots) have a positive score, and tend to be distributed toward the right end of the plot. Most of the compounds with possible or conditional TdP risk (TdP risk category 2 or 3, orange and yellow dots, respectively) have a positive score, except Clozapine, Dasatinib, Paroxetine, Saquinavir, Voriconazole. The TdP score is also dependent on the maximum considered concentrations. Again, higher concentrations lead to an increase in sensitivity while decreasing specificity. A comparison between the TdP scores computed up to 30x and 100x EFTPC_max is included in the Supplementary Material (Figure S33).

In silico drug trials are likely to be used as an additional tool for drug safety assessment in combination with experimental methods. It is therefore important to evaluate the consistency between in silico results and experimental data. Thus, simulation results for the 15 reference compounds with varied actions on ion channels (Table 1) were compared against recordings obtained using rabbit wedge and hiPS-CM preparations, as two techniques considered in safety pharmacology. In Figure 5, changes in QT interval duration in rabbit wedge, CTD₉₀ in hiPS-CMs and in silico APD₉₀ (from red to green, left side) were compared to evaluate drug-induced changes in repolarization. Changes in QRS complex duration in rabbit wedge and in silico dV/dt_MAX were quantified to evaluate drug effects on depolarization (Figure 5, from purple to blue, right side). Negative variations in dV/dt_MAX are considered as positive changes in depolarization time (opposite sign), to facilitate the comparison.

Drug-induced effects on biomarkers are in overall agreement for all three methodologies. Figure 5 presents consistent increase/decrease of QT, CTD₉₀ and APD₉₀, as well as consistency between positive changes in QRS and reduction of dV/dt_MAX. Variations are generally larger in the in silico APD₉₀ than in QT interval in rabbit wedge, indicating higher sensitivity of the in silico assay, and the wider range of ionic scenarios evaluated in the virtual population than in the limited number of experiments.

For Verapamil at 0.1 μM, small changes were observed in both QT interval in rabbit wedge and simulated APD, whereas CTD₉₀ in hiPS-CMs was reduced. Such a reduction in CTD₉₀ in the hiPS-CMs is also accompanied by a significant increase in beating rate (CTBR +61%), which does not occur in silico and in the rabbit wedge experiments as these are paced externally.

For Lidocaine and Mexiletine, their main effect is fast I_Na block, which results in a decreased dV/dt_MAX in the in silico models, and a wider QRS complex in the rabbit wedge ECG. Furthermore, both in vitro CTD₉₀ and in silico APD₉₀ are prolonged, whereas QT interval decreases slightly (Lidocaine) or remains unchanged (Mexiletine) in rabbit wedge, suggesting that in vitro CT and in silico AP are more prone to display prolongation for non-selective Class I anti-arrhythmic drugs. It is worth noticing that the AP prolongation in silico is reduced when inhibition of the I_NaL current is taken into account in the simulations, corresponding to Lidocaine II and Mexiletine II (APD₉₀ +84% vs. +68%, Mexiletine I vs. Mexiletine II, 100 μM). The tendency for prolongation under Na⁺ block of in silico AP and hiPS-CMs CT is confirmed also for another class I anti-arrhythmic drug (Phenytoin), which causes negligible changes in both CTD₉₀ and in silico APD₉₀, and a decrease in QT interval in rabbit wedge.

Bepridil constitutes a good example of multichannel block, intended to block I_CaL, but also affecting I_Kr and I_Na. In the rabbit wedge, the QT interval is relatively prolonged at low concentrations (0.3–1 μM), compared to more selective Ca²⁺ blockers (e.g., Nimodipine and Nisoldipine), and it goes back to normal (+4%) at 10 μM. Both hiPS-CMs CT and in silico AP are prolonged in a dose-dependent manner, confirming once again the higher sensitive to I_Kr block of these techniques.

Interestingly, BaCl₂ effects are of smaller magnitude in hiPS-CMs compared to rabbit wedge preparation and in silico models: relevant CTD₉₀ prolongation was detected only at 100 μM (+47%), while up to 10 μM (more than 2-fold BaCl₂ IC₅₀ for I_K1) the drug-induced effects on CT were negligible. This may be due to differences in I_K1 expression between the cell types considered (Liang et al., 2013; Kim et al., 2015).

In silico results for compounds with multiple sets of IC₅₀ and h are in overall agreement with each other, even if the magnitude of drug-induced changes may vary. As an example, the decrease in dV/dt_MAX for the three variations of Flecainide is almost the same at 10 μM (−64, −61, and −68%, respectively), while for lower concentrations the differences between the three inhibition profiles are more noticeable (e.g., −25, −7, and 0% at 1 μM, respectively).

HL, JR, AH, and DG are employees of Janssen Pharmaceutica NV. RG is an employee of Oxford Computer Consultants, Oxford, UK. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.