The definition of drug groups according to their torsadogenic risk is critical for the development of the TdP risk classifiers. Different torsadogenic definitions from previous classification studies are listed below.

Redfern et al. assigned drugs to five categories based on the number of reports of TdP in humans, QT prolongation and TdP associated withdrawal from the market (Redfern et al., 2003). The five categories are:

R1: Class Ia and III antiarrhythmics with QT prolongation as intended effect.

Arizona Center for Education and Research on Therapeutics (AZCERT) maintains a list of drugs associated with QT prolongation/TdP risk (Woosley et al., 2017) (https://crediblemeds.org/) and has categorized the drugs into three groups:

CM1: Drugs with known risk of TdP. These drugs prolong the QT interval and are clearly associated with a known risk of TdP, even when taken as recommended.

Champeroux et al. assigned drugs into three categories based on the number of reports of TdP cases associated with the drug (Champeroux et al., 2005):

CH1: Drugs with numerous reports of TdP.

Based on a general consensus, a working group formed under the Comprehensive in vitro Proarrhythmia Assay (CiPA) initiative picked 28 compounds and categorized them into three groups (Colatsky et al., 2016; Fermini et al., 2016) for testing/training of the classifiers under the new CiPA paradigm:

CP1: Drugs with high risk.

To consistently compare with other methods we attempt to use the binary TdP definitions [i.e., a drug is either torsadogenic (TdP+) or non-torsadogenic (TdP−)] as in the original publications (see Table 2 and the Supplemental Material for further details regarding the exact risk definition used for the particular datasets). For dataset 7 where tertiary definition is reported, binary TdP definition was defined by placing CP1 and CP2 drugs into torsadogenic (TdP+) and CP3 drugs to non-torsadogenic (TdP−) categories. In the case of the merged dataset or the datasets that lacked binarized TdP definitions (Crumb et al., 2016), we assigned drugs as TdP+ or TdP− using a similar approach as in Lancaster and Sobie (2016). The drugs which fell in the known risk category (CM1) in the CredibleMeds database, R1, R2, and R3 category in Redfern et al. (2003) or drugs with several reports for TdP (CH1) (Champeroux et al., 2005) were assigned to TdP+. For the remaining drugs we categorized the drug as TdP+ if a warning for TdP associated with QT prolongation appeared on its package label (http://dailymed.nlm.nih.gov/). The risk categorization for the drugs are provided in the Supplemental Material. Paroxetine has a warning for TdP in the label and Imipramine was assigned to CH1 category in Champeroux et al. (2005). These two drugs were also assigned as TdP+ in Lancaster and Sobie (2016). Here, we defined them as TdP− as these two compounds are not directly associated with QT prolongation or TdP. These drugs inhibit CYP2D6 and increase plasma concentrations of TdP positive drugs such as Thioridazine (http://dailymed.nlm.nih.gov/). Sometimes alternative definitions were also considered and are explicitly defined in the manuscript.

The in-vitro ion channel assay data is converted to drug-induced block of ion channel (direct features) using

where IC_50,channel is the drug concentrations at which the whole-cell current through particular channels is reduced by half, C_Drug is the concentration of the drug and h is the Hill-coefficient. The Hill coefficient values were taken as reported in the original datasets. The IC₅₀, Hill coefficients, and EFTPC values for each of the datasets are given in the Supplemental Material. Note that Hill coefficient values had little impact, and even fixing Hill coefficient to 1 for all the drugs did not significantly alter the observed classification accuracies (results not shown). The drug-induced blocks of ion channels are used as input features for the machine-learning based classifiers or utilized to scale the maximum conductance (g_channel) of the particular ion channels in the in-silico models, assuming that drug-induced effect on multiple ion-channels are well represented by a simple conductance-block model (Mirams et al., 2011),

The effects of drug-induced modulation of multi-channel conductance (Equation 2) on the action potentials (AP) and calcium transients are simulated for all of the compounds using two versions of human ventricular myocyte models (O'Hara et al., 2011; Dutta et al., 2016). Ohara et al. model (OHR) (O'Hara et al., 2011) was picked as it has been chosen as the consensus base model for proarrhythmic risk assessment (Colatsky et al., 2016). OHR model was also shown to have the best predictive capability for TdP risk classification among the few tested models in Lancaster and Sobie (2016). We also utilized the modified version of the OHR model (OHRmv) which has been shown to better fit APD-rate dependence experimental data under drug block conditions, particularly improving the effect of I_NaL block on action potential prolongation (Dutta et al., 2016). Several derived features are extracted from AP and calcium transients, and TdP risk classifiers are constructed using these derived features. The details on the simulation protocols and the computations of the derived features are reported in section 2.6 that describes in-silico simulations.

We propose a two-step approach for TdP risk prediction. In the first step, we classified the drugs into non-torsadogenic and potentially torsadogenic categories. We performed the classification in the first step considering solely the block of hERG channels. Using a Redfern-like criteria (Redfern et al., 2003), we obtained the ratio between the drug concentration that produces 60% block of the hERG channel current (IC_60,hERG) and the drug EFTPC (i.e., hERGratio=IC60,hERGEFTPC). The motivation for using IC_60,hERG was that the 60% block of hERG channel currents triggers pause-induced EAD in the mid cell type of OHR and OHRmv models (2 Hz pacing rate) at a quiescent interval greater than 700 ms. For example, a drug with hERG IC₆₀ of 500 nM and EFTPC of 1 nM would yield a threshold of 500. Since this drug's EFTPC would be far from the critical hERG block concentration, we would classify this compound as non-torsadogenic. Previous classifiers based on EAD appearance for different datasets have shown a wide range of thresholds (approximately 30 × −200 × EFTPC) for achieving best TdP risk predictions (Christophe, 2013; Okada et al., 2015). Hence, we tested four different thresholds for hERG ratio of 50, 100, 150, and 200 for all the datasets and chose the one that gives the best classification accuracy. For the remaining drugs with EFTPC above the critical hERG block concentrations (hERG ratio less than the threshold), the role of multi-channel block was examined in the second step using logistic regression classifiers ignoring the EFTPC values of the drugs. The regression classifiers employed as inputs either the Block_channel of additional non-hERG ion channels (direct features) or the features derived from the simulated calcium transient and the AP in the ventricular myocyte models, at drug concentrations equal to IC_60,hERG. Such two-step classifier partially solves the problem of EFTPC variability, restricting EFTPC usage only to the first step to primarily discard the drugs that produce insufficient hERG block at extremely high concentrations by classifying them as non-torsadogenic. Therefore, the moderate variations in EFTPC values of the drugs would only matter for a very small population of the drugs with hERG ratio close to the threshold in the classification. A summary of the two-step approach is given in Figure 1.

In order to compare the performance of the classifiers based on the two-step approach to the performance of the classifiers based on features obtained at actual drug EFTPC concentrations, we also constructed TdP risk classifiers using the direct and derived features at reported maximum EFTPC of the drugs. They are referred in the current paper as one-step classifiers (hERG ratio is not utilized for these classifiers).

We utilized statistical/machine-learning models for binary classification of the drugs into TdP+ or TdP− categories. The binarized torsadogenic definitions for each drug were used to train/test the classifier models. Here, we used Logistic regression model. SVM and neural network models were also tested and resulted in comparable classification accuracies (results are shown in the Supplemental Material). Python's scikit-learn package (Pedregosa et al., 2011) (http://scikit-learn.org/stable/) was used to train/test different models. Here, we present results for logistic regression models only as other methods produced similar results. The generalized model equation is described as

where Feature represents the input metrics to the model (either direct feature or the derived feature), n is the number of input features used to train/test the model, and β₀ and βi (i = 1, 2, .., n) are the parameters to be determined. The predictive power of the model was evaluated by the leave-one-out (LOO) cross validation technique.

A two-dimensional TdP risk map with hERG ratio (IC60,hERGEFTPC) on the x-axis and summation of one or more features (block of I_CaV and I_NaL) on the y-axis were constructed for visualization of the two-step (MCB@EAD) classifier. The hERG ratio threshold and regression coefficients from the two-step classifier are used to generate the two-dimensional risk maps. The hERG ratio (step 1 in the two-step classifier) that provides the best classification in the two-step classifier is used to set the threshold along the x-axis. Drugs that fall in the region above the hERG ratio threshold are considered to be non-torsadogenic. For the drugs with hERG ratio less than the threshold, the coefficients of the logistic regression model in step 2 of the two-step classifier are used to determine the threshold along the y-axis of the risk map. An example of separating hyperplane that would be obtained from the second step of the two-step classifier is given by

where block_ICaV and block_INaL are the blocks of I_CaV and I_NaL, respectively. feature_i are additional input features of the model, such as drug trapping parameters. β_ICaV, β_INaL, and β_{feature, i} represent the regression coefficients. The regression coefficients obtained from the step 2 are normalized to the coefficient β_ICaV for the I_CaV block. This gives

Thus, the left hand side of Equation (5) is plotted on the y-axis, and the ratio βinterceptβICaV determines the threshold along this axis. For example, if only block of I_CaV is taken into consideration, the risk map has the values of I_CaV along the y-axis and hERG ratio along the x-axis.

Ternary classification for Dataset 7 requires multiple hyperplanes with different regression coefficients to separate the high, low and intermediate risk drugs. To represent the ternary classification in a two-dimensional risk map similar to the binary classification, we summed the different features (Feature_sum) assuming identical weights for each of the features (i.e., β₁ = β₂ = .. = β_n equal to β_f) reducing the classification model to 11+exp-(β0+βfFeaturesum). For Dataset 7, assuming identical weights for different features resulted in similar accuracy to a multinomial logistic regression classifier while providing a simpler visualization in one two-dimensional plot as in the binary classification. The ratio -β0βf obtained after training the model is used to set the two thresholds along the y-axis. Along the x-axis arbitrary hERG ratio of 25, the value slightly greater than the maximum hERG ratio observed for high risk drugs in Li et al. (2017) and Fermini et al. (2016), was utilized to separate low and intermediate risk drugs from the high risk drugs. The hERG ratio of 150 obtained from the merged dataset was utilized to separate the high and intermediate risk drugs from the low risk drugs.

The alteration in the action potential and calcium transients at the cellular level arising from drug-induced multi-channel block were simulated for all the compounds (Dataset 8) using the OHR model. Using a similar approach as in Lancaster and Sobie (2016), simulations were carried out at three pacing rates (0.5, 1, and 2 Hz) for each of the endo, mid and epi cell types resulting in 9 simulations per drug, and 13 metrics were obtained from the AP and Ca²⁺ transients. Simulations were carried out for 1,000 beats to allow the models to reach the steady state. The 13 metrics obtained from the in-silico models are listed below:

Upstroke velocity

We systematically construct the classifiers on each of the 13 derived features extracted from the action potentials and calcium transients at two different drug concentrations (i.e., at C_Drug = EFTPC and C_Drug = hERG IC₆₀).

The onset of TdP is usually preceded by a sudden reduction of the heart rate (i.e., by pauses or long cycle lengths) (Neal Kay et al., 1983; Viskin et al., 2000). Here, we test the generation of pause-induced EADs, that are implicated as triggers of TdP (Viswanathan and Rudy, 1999; Liu and Laurita, 2005), in simulations of drug-induced multi-channel blockage in the ventricular myocytes models (O'Hara et al., 2011; Dutta et al., 2016). The basic protocol was similar to that in Viswanathan and Rudy (1999), where stimulation of the cell is carried out 200 times at a constant cycle length of 500 ms. After 200 stimuli, an additional stimulus was applied following a pause of 1,000 ms. Drug-induced EAD development was tested at drug concentrations = IC_60,hERG in the mid cell type. EAD analysis was also performed at C_drug = IC_60,hERG for combinations of I_CaV, I_NaL, and I_Ks blocks ranging from 0–100% with step of 10% resulting in a set of 1,000 simulations. TdP risk prediction was carried out using ability of drugs to induce EADs as a classification criteria (EAD+: drugs that induce EADs at 60% hERG block concentrations, EAD−: drugs that do not induce EADs at 60% hERG block concentrations). Figure 1 illustrates the classification based on the drug EAD risk. The drugs with hERG ratio greater than the hERG ratio threshold determined for the particular dataset (using the two-step approach on the direct features) were considered EAD−. For the remaining drugs the block of ion-channels was calculated at drug concentration equal to IC_60,hERG and overlaid on the parametric space obtained from EAD analysis at varying combinations of I_CaV, I_NaL and I_Ks blocks (Figure 2A), for both the OHR and OHRmv models, to determine whether a drug will induce EAD or not at 60% hERG block concentrations.

The system of ordinary differential equations were solved using the rapid integration scheme (a combination of forward Euler, Rush-Larsen method, Rush and Larsen, 1978 and adaptive time-step) proposed in the original model (O'Hara et al., 2011). For the EAD simulations rapid integration scheme proposed in O'Hara et al. (2011) yielded different results to gold standard simulations with fixed time step of 0.001 ms. Hence, for the EAD simulations we utilized forward Euler method with a time step of 0.001 ms. Execution scripts were written in C++.

A short-long cycle length (i.e., pause) often precedes the onset of TdP (Neal Kay et al., 1983; Viskin et al., 2000). The pause is known to facilitate the formations of EADs (Viswanathan and Rudy, 1999; Liu and Laurita, 2005). Here, we test the effects of drug-induced block of different channels on triggering of pause-induced EADs. Block of hERG channel causes prolongation of action potential and can result in the generation of EADs. In the OHR and the OHRmv models, the amount of hERG block required to induce pause-induced EADs is reduced with increase in the duration of the pause. hERG block by 57 and 55% induced EADs in the mid-cell paced at 2 Hz (500 ms pacing cycle length) following a 700 ms pause in the OHR and OHRmv models, respectively. The amount of hERG block required to induce EADs in the OHR and OHRmv model following 1,000 ms pause is reduced to a 47 and 46% block, respectively. Under these critical blocks of hERG channel, modifications of other channels may promote or inhibit the EADs. First, we tested individually the effects of block of six ion-channel currents (I_CaV, I_NaL, I_Ks, I_K1, I_Na,fast, and I_to) in promoting or inhibiting pause-induced EADs at fixed blocks of hERG channels that are marginally above (60%) and below (40%) the critical value (48%) of hERG block that is required to induced EADs in the models. At 60% hERG block, block of I_CaV (> 30% for both the OHR and OHRmv models) and I_NaL (> 60% for the OHRmv model) resulted in suppression of EADs. The remaining five channels had no inhibitory effects. At 40% hERG block, blocks of I_Ks and I_K1 led to induction of EADs with lower block of I_Ks currents promoting triggering of EADs compared to I_K1. The block of remaining five channels did not result in EADs at 40% hERG block.

For visualization of combined effects of the channels on EAD induction, we performed EAD analysis for varying combinations of block of the three most sensitive non-hERG channels (I_CaV, I_NaL, and I_Ks) regulating EAD generation, each ranging from 0 to 100% at C_drug = IC_60,hERG. Figure 2A represents the EAD test for the models. The red region represents the set of combinations of blocks for these currents that resulted in the EADs in the OHR model (EAD+ region). The region in green covers the parameter subspace where no EADs were observed in the OHR model (EAD− region). Separation of the EAD+ region and EAD− region is outlined by the blue and yellow surfaces for the OHR and OHRmv model, respectively (Figure 2A). Among the non-hERG channels, block of I_CaV had the highest modulatory effects on EAD generation under normal conditions. Under critical hERG block (60%), block of I_CaV by more than 30% resulted in suppression of the EADs in both the OHR and OHRmv models (Figure 2A). Block of I_NaL in the absence of the block of I_CaV did not result in suppression of the EADs in the OHR model. Even in the OHRmv model with improved I_NaL formulation, block of more than 60% of the late sodium current was required to suppress pause-induced EADs. Block of I_Ks currents increased the amount of block of I_NaL and I_CaV currents required to prevent EADs. A similar parametric space for EAD generation was also analyzed under enhanced late sodium currents (Figure 2B), associated with LQT3. In the OHRmv model the conductance of the late sodium current was doubled. Block of the hERG channel by 30% was enough to induce EADs in the OHRmv model under enhanced late sodium currents. EAD generation was examined at this 30% hERG block and is shown in Figure 2B. For the simulations with enhanced late sodium currents, the difference in the effects of I_NaL and I_CaV was significantly reduced. At 30% hERG block in the OHRmv model, block of greater than 20 and 30% was required for EAD suppression by I_CaV and I_NaL, respectively (Figure 2B).

Table 3 lists the accuracy of TdP risk prediction using EADs as the classification criteria. As an illustration, the data points for the drugs in one of the datasets (Dataset 5) are overlaid on the the parametric space in Figure 2A at hERG blocks of 60% to visualize the agreement between drugs EAD and TdP risks. The drugs with positive TdP risk are shown in red, while the drugs with negative TdP risk are shown in green. Actual AP profiles were also simulated for drugs in Datasets 5 and 7 at 60% hERG block concentrations taking into account the block of all seven channels. Figure 2C gives representative examples of the AP profiles for the multi-channel block of six of the drugs from Dataset 7 simulated in the OHR model. Our results show a good concordance between drugs torsadogenic risk and its propensity to induce EADs (i.e., most of the torsadogenic drugs resulted in pause-induced EADs in the models while no EADs are observed for majority of the non-torsadogenic drugs at 60% hERG block drug concentrations) across all the datasets. Majority of the datasets (Datasets 1, 2, 4, 6) contains values of blocks for two non-hERG channels (I_CaV, I_Na,fast) of which only I_CaV had an impact on EAD occurrence in the OHR and OHRmv models. Hence, for these datasets, OHR and OHRmv models give identical accuracies as the drugs with greater than 30% I_CaV block would result in EAD suppression in both of these models. For the Datasets 3, 5, and 7, drug-induced blocks of multiple ion channels are reported. However, very few drugs among these datasets are located in the region between the EAD risk decision surfaces for the OHR (blue surface) and OHRmv (yellow surface) models (i.e., the region in the parameter space that would result in a different prediction between the OHR and OHRmv models, Figure 2A). Hence, similar accuracies are observed for both of the models across all datasets. For the Dataset 7, Ranolazine (TdP−) was the only drug that was located on the negative side of the EAD risk decision surface for the OHRmv and to the positive side of the decision surface for the OHR, and hence predicted correctly by the OHRmv model but not by the OHR model. On the contrary for Dataset 5, Quinine (TdP+) ended up on the negative side of the EAD risk separating surface for the OHRmv and to the positive side of the EAD risk decision surface of the OHR model resulting in its incorrect prediction using the OHRmv model (Figure 2A).

Although the biophysical models can provide mechanistic insights underlying TdP genesis, the benefits of using biophysical models in terms of classification is unclear. Here, we wanted to examine the performance of the classifiers built on direct features using the proposed method. The predictive power of the TdP risk classifiers built on the direct features using the proposed method (MCB@EAD) is shown in Table 4 (two-step classifier column). Classification scores on the direct features at EFTPC are reported for comparison (one-step classifier column). The predictive ability of the two-step classifier was comparable or better than those for the classifiers built on the various features in the original datasets and also better than those for the one-step classifier. Most datasets comprise in-vitro assay data for drug-induced block of I_Kr, I_Na,fast, and I_CaV (Datasets 1, 2, 4, and 6). Drug-induced blocks of additional channels were reported in three datasets (Datasets 3, 5, and 7). The classifiers built using the block of I_Kr and I_CaV as inputs provided high TdP risk prediction scores (Table 4). Utilizing the block of I_NaL as an additional input feature improved the prediction for the two datasets (Dataset 3 and Dataset 7) by classifying correctly one more drug, Ranolazine, as TdP−. On the contrary, in Dataset 5, addition of the block of I_NaL to the features reduced the number of correctly classified drugs by one (classifying Ritonavir incorrectly as TdP−). Taking into account the block of additional ion channels (I_Na,fast for Datasets 1, 2, 4, 6 and I_Na,fast, I_Ks, I_to, I_K1 for Datasets 3, 5, 7) did not provide any further improvement in the performance of the classifier for all datasets.

For visualization of our two-step approach we presented two-dimensional risk maps. Drug-induced block of the I_CaV or the sum of I_CaV and I_NaL blocks is plotted against the hERG ratio (IC60,hERGEFPTCdrug) for each drug in a two-dimensional risk map (Figure 3). The values of the hERG ratio (x-axis) and the block of I_CaV (y-axis) that provide the best discrimination vary widely across the datasets (Figure 3). Using the threshold of 22 and 100 for the block of I_CaV and hERG ratio, respectively, resulted in the perfect classification for Dataset 1 (Mirams et al., 2011). However, for Dataset 2 the hERG ratio threshold of 200 and threshold for I_CaV of 57 provided the best classification. Looking at the risk maps in Figure 3 we can see that the thresholds that provide best classification accuracy vary across datasets resulting in variations in the obtained high risk zones across datasets. For datasets 3, 5, 7 where we also considered block of I_NaL as one of the input feature in addition to the block of I_CaV the value on y-axis of the risk map were reported using Equation (5). A two-dimensional risk map for Dataset 5 for two alternate TdP definitions is shown in Figure 4. On the two-dimensional risk map with only I_CaV block on the y-axis, we also highlight, using a blue rectangular outline, the separation between the region with EAD presence and absence observed in the in-silico simulation at varying I_CaV blocks. For example, at critical hERG block simulations in the model at block of I_CaV less than 30% would result in EADs. Block of I_CaV by more than 30% results in suppression of the observed EADs. There is approximate correspondence between the EAD observance region in the in-silico model and the high torsadogenic region obtained from classifying the direct-features using machine learning.

We examine the performance of classifiers from various derived features extracted from the simulated action potential and calcium transient in predicting TdP risk. Drug-induced multi-channel block was simulated in the OHR model for all the compounds in the merged dataset (Dataset 8) considering the block I_Kr and I_CaV channels. Simulations were also carried out taking into account drug-induced block of all the channels with available IC₅₀ values. Taking into account the block of additional channels resulted in similar accuracies and are reported in the Supplemental Material. Figure 5 shows the performance of the logistic-regression classifier to discriminate TdP+ and TdP− drugs (Figures 5A,B) and to predict drug-induced EADs (Figures 5C,D). The classifiers were built on 13 features extracted from the steady-state APs and Ca²⁺ transients (for 3 cell types at 3 pacing rates). The derived features were extracted at either EFTPC of the drugs (Figures 5A,C) or at concentrations at which each drug would produce 60% hERG block (C_drug = IC_60,hERG) (Figures 5B,D). Among the various derived features obtained at C_drugs = EFTPC, diastolic Ca²⁺ levels provided the best discrimination score between the TdP+ and TdP− drugs (~79% accuracy at 1 Hz in epi cell type, Figure 5A) in agreement with a previous report (Lancaster and Sobie, 2016). Classifier was also built using APD₅₀ and Diastolic Ca²⁺ together as inputs (the combination that provided the best prediction in Lancaster and Sobie, 2016) but did not give an improved classification for the merged dataset. Several derived features performed well (>90% accuracy) for EAD risk prediction (Figure 5). However, at drug concentrations equal to IC_60,hERG (C_drug = IC_60,hERG), each of the 13 features from the Ca²⁺ transient and AP provided high classification scores (~85% maximum) for TdP risk assessment (Figure 5C). The highest accuracy for each of the features was obtained at different pacing rates (0.5, 1, and 2 Hz) and for different cell types (endo, mid and epi). The maximum classification scores to discriminate the drugs that induced EAD in the model from the drugs that did not induce EADs was 100% (Figure 5D) for each of the features. Our results suggest that at fixed hERG block concentrations where trigger events such as EADs arise, several derived features obtained from the model including features from the Ca²⁺ transient can provide good TdP risk prediction. These derived features also highly correlate with EADs (Figure 5). However, the derived features extracted from simulations of the drug-induced effects in ventricular myocyte model did not result in much improvement in TdP risk assessment over the classifiers built on the direct-features using the proposed method (Table 4). Moreover, combining the direct and derived features to build the TdP risk classifiers also did not improve the classification performance (results are not shown).

Our results show that classifiers built on the direct features serve as excellent predictors of TdP risk of the drugs categorized into binary risk groups. However, a working committee under the CiPA initiative led by FDA has recently categorized 28 drugs into tertiary risk categories (low, medium, and high risk compounds) (Colatsky et al., 2016; Fermini et al., 2016). Hence, we test the predictive capabilities of the classifier based on direct features to classify the drugs categorized into tertiary risk categories. Figures 6A,B show a two-dimensional risk map for the 12 drugs in Li et al. (2017). These drugs have been analyzed previously using modified OHR model that incorporate dynamic hERG channel interactions. The 12 drugs are a subset of the 28 drugs categorized into three risk categories (CP1, CP2, and CP3) under the CiPA initiative. As a first step, we developed a TdP risk map using only the block of I_CaV channel. The hERG ratio threshold of 150 and the threshold of I_CaV block of 45%, the values that provided best classification for the merged dataset, were utilized for classifying the low high and intermediate risk drugs from low risk drugs Figure 6A. Arbitrary value of hERG ratio and I_CaV block of 25 and 15, respectively (a value greater than the maximum hERG ratio and maximum block of I_CaV among the high risk drugs in Datasets 7 and 9) was utilized to separate high risk drugs from the low and intermediate risk drugs. Three of the four drugs in low risk (CP3) category were classified correctly. Ranolazine was the only misclassified drug. The boundaries of red zone were defined to include all high risk drugs and hence all the drugs from CP1 category were correctly classified. However, several drugs in intermediate risk were incorrectly classified. Next, we built regression classifier using as input metric the sum of block of I_CaV, I_NaL channels and the degree of drug-trapping parameter which was shown to be essential to improve risk prediction of intermediate risk drugs in the original dataset. Figure 6B shows the risk map built using this metric. The threshold along the y-axis was obtained from the regression coefficients. The hERG ratio threshold of 25 and 150 were utilized as before. Including the degree of drug trapping characterized by open-bound/closed-bound ratio for the drugs at steady-state (Li et al., 2017) as one of the features in addition to the blocks of CaV and NaL channels, resulted in the perfect separation of the 12 drugs in 3 categories (Figure 6B).

We employed the same approach to test all of the 28 drugs (Fermini et al., 2016) categorized in CP1, CP2, and CP3 categories under the CiPA initiative. In-vitro assays for 12 of these 28 drugs were reported in Crumb et al. (2016) and analyzed in Li et al. (2017). To augment this dataset, we extract the IC₅₀ values for hERG and CaV blocks from Datasets 1, 2, 3, 5, 6, 7 resulting in characterization of 26 of the 28 drugs categorized under the CiPA initiative. Two drugs (Azimilide and Loratidine) were absent in all of the datasets analyzed here and hence were not taken into consideration. We used the mean value of the block if the drug was present in more than one dataset. The final dataset of the 26 drugs is reported in the Supplemental Material. Figure 6C show the two-dimensional risk maps for the 26 drugs with hERG ratio on the x-axis and the block of CaV on the y-axis. The 45% I_CaV block threshold and hERG ratio threshold of 150 yielded almost perfect binary classification (high and intermediate vs. low risk drugs) with only 2 (Ranolazine and Tamoxifen) drugs of the 8 from the CP3 category and one drug (Clozapine) of the 16 drugs from CP1 and CP2 cateogry (high and intermediate risk drugs) classifying incorrectly. However, no clear separation was observed amoung the drugs in CP1 and CP2 categories, with several drugs from CP2 category ending up in the high risk region. Considering additional features such as I_NaL and the degree of drug trapping (if either of the features were not available their value was set to zero) and utilizing the threshold obtained from training the Dataset 7 (Figure 6B) resulted in only 3 intermediate risk drugs (Clarithromycin, Domperidone, Droperidol) and 1 low risk drug (Tamoxifen) of the 26 drugs to be misclassified (Figure 6D). The classifier already performs well in spite of testing data from heterogeneous sources with some missing values. Further refinement of the method may be possible when a dataset is available with all 28 drugs characterized with a consistent methodology.

Different binary definitions for the drug's torsadogenic risk have been used across the literature (Lancaster and Sobie, 2016; Ando et al., 2017; Wiśniowska and Polak, 2017). Tables 5, 6 list the different classification accuracy scores obtained for four different binary TdP definitions (Datasets 5 and 8, respectively). Classifiers were constructed on the block of one or multiple-ion channels as inputs at critical hERG block concentrations (C_drug = IC_60,hERG). The various definitions not only resulted in variability of classification scores (Tables 5, 6), but also changed the role of different ion channels in accurate TdP risk classification (Table 5). Using the block of I_CaV currents provided the best accuracy scores under two of the four definitions (Target 1 and Target 2) for Dataset 5. Including the effects on additional ion-channels did not improve classification scores for these two definitions. Taking into account the block of late sodium currents in addition to CaV channels, provided the best classification accuracy for remaining two of the four TdP definitions (Target 3 and 4) (Table 5).

Although the role of multiple ion-channels have been suggested for improved TdP risk prediction, classifiers have been primarily built on the blocks of I_Kr, I_CaV, and I_Na,peak currents (Mirams et al., 2011; Christophe, 2013, 2015; Kramer et al., 2013; Lancaster and Sobie, 2016). A recent assay reports the drug-induced effects on seven ion channels (Crumb et al., 2016) providing an opportunity to identify the ion channels that are important for pro-arrhythmic risk assessment. The results of our parametric simulations of EAD indicate a potential role of block of I_CaV, I_NaL, and I_Ks currents, in addition to I_Kr, for determination of torsadogenic risk of the drugs (Figure 2). It should be noted that the EAD simulations results are highly dependent on the ventricular myocyte model. For example, block of late sodium current in the OHRmv plays a more prominent role in regulation of AP sensitivity to EADs as compared to the OHR model (Figure 2A). The I_Ks plays a much bigger role in Ten Tusscher and Panfilov model (Ten Tusscher and Panfilov, 2006) than OHR model in regulation of APD as shown in Mirams et al. (2014). The present datasets have limited examples of block of I_CaV, I_NaL, and I_Ks currents in the same compounds. Moreover, among the non-hERG channels, the regulatory effect of I_CaV block was the highest with the block of I_CaV by only 30% resulting in EAD suppression at critical hERG current block (Figure 2). The classifiers constructed on the block of I_Kr and I_CaL provided the best discrimination between torsadogenic and non-torsadogenic drugs for the majority of the datasets tested here, including the dataset where drug-induced effects on seven ion channels were reported (Figure 2, Tables 4, 5). Our results suggest that among different channels, examination of block of I_CaV and I_Kr might be the most critical for TdP risk prediction. Relative block of I_CaV and I_Kr (among the three currents measured in the in-vitro assay) was shown to provide the best risk prediction, with no role of peak/fast sodium currents in improving the classification in Kramer et al. (2013).

Examination of late sodium block can be important for the drugs with low to moderate I_CaV block as these drugs would be predicted TdP+ if only I_Kr and I_CaV block are considered for risk prediction. Moderate to high block of I_NaL by these drugs can result in suppression of EADs (Figure 2) indicating lower TdP risk. Earlier datasets did not report values for I_NaL block. Dataset 5 reports the value of drug-induced block of seven ion-channels, including the block of I_NaL. Among the drugs with low to moderate I_CaV in Dataset 5, only three drugs [Ranolazine(CM3), Toremifene(CM2) and Quinine(CM3)] have greater than 30% I_NaL block at critical hERG block concentrations. For the limited data with inconsistent risk categorization, taking into account the I_NaL block did not improve predictive power of the classifiers (Tables 4, 5). A small improvement in TdP prediction was observed for Datasets 3 and 7 (Okada et al., 2015; Li et al., 2017) when considering drug-induced block of I_NaL as one of the input features by correctly classifying Ranolazine (the only drug with high late sodium block in absence of I_CaV block) in both the datasets (Table 4 and Figure 3). The limited data and inconclusive/minor improvement in torsadogenic risk classification make it difficult to ascertain the role of ion channels such as I_NaL and I_Ks in predicting TdP risk.

In-silico biophysical models can be thought of as a complex non-linear transfer function, which translates the drug-induced multi-channel block effects at channel level (input) to alterations in APs and calcium transients at cellular/tissue levels (output). Several in-silico electrical biophysical models of human ventricular cell models have been published over the last decade (e.g., Ten Tusscher and Panfilov, 2006; Grandi et al., 2010; O'Hara et al., 2011; Himeno et al., 2015). TdP risk classification on features extracted from the drug-induced responses in isolated cell (Mirams et al., 2011; Christophe, 2013, 2015; Lancaster and Sobie, 2016), tissue (Trenor et al., 2013; Kubo et al., 2017) or organ level (Okada et al., 2015) computational models can provide physiological/mechanistic insights. Moreover, in-silico models serve as an excellent tool for evaluation of drug-safety in diseased conditions (Trenor et al., 2013; Kubo et al., 2017). Our simulations in the OHRmv model under pathological conditions (enhanced late sodium currents) reveal that EAD can appear at significantly lower drug concentrations as only 30% hERG block was required to induce pause-induced EADs under pathological conditions compared to 60% hERG block under normal conditions. Moreover, the modulatory effects of non-hERG channels on EAD induction was also significantly different for simulations under pathological conditions (Figure 2). On the other hand, biophysical models show considerable differences in their formulations and can lead to different predictions based on the chosen model. Simulations of EAD generation show a significant difference between the surfaces separating EAD+ region from EAD− region, that are obtained from the OHR and OHRmv models (Figure 2). Recently, several efforts have been carried out for optimization of in-silico cardiac cell models for pro-arrhythmia risk assessment (Dutta et al., 2016; Mann et al., 2016; Li et al., 2017).

Statistical/machine learning classifiers that use measured in-vitro block of multiple cardiac channels (direct features) as their input (Kramer et al., 2013; Mistry et al., 2015) demonstrated comparable accuracy as compared to TdP risk classifiers built on derived features, questioning the need of additional complexity provided by the in-silico models. On the contrary, the study by Okada et al. suggested the need of highly detailed three-dimensional cardiac models for pro-arrhythmic risk assessment and showed relatively low predictive ability using the direct features and also certain derived features from the in-silico lumped parameter (zero-dimensional) cellular models (Okada et al., 2015). Derived features from in-silico simulations that incorporate dynamic drug-hERG channel interactions were shown to improve prediction of TdP risk (Li et al., 2017). For all the datasets tested here, including the datasets in Okada et al. (2015) and Li et al. (2017), we showed that the classifiers built on the direct features performed equally or better than the previously developed classifiers on the derived features (Table 4). Our results show that for currently available in-vitro assay datasets simple models based on the direct features can provide similar accuarcy to more complex models based on derived features. It should be noted that our two-dimensional risk classifiers on the direct features also utilized insights gained from the computational models (the direct features are examined at critical hERG block concentrations where EAD can arise in the in-silico models). Our parametric simulation for EAD induction highlights one of the possible reasons for the insignificant improvement in predictive power of classifiers built on the derived features from the in-silico models. Although a non-linear surface is obtained from the in-silico models separating the EAD+ and EAD− regions (Figure 2), a hyperplane

constructed using direct features can result in nearly identical separation, where a, b, c, and d are the parameters of the hyperplane and block_ICaV, block_INaL, and block_IKs are the values of block of I_CaV I_NaL, and I_Ks, respectively. Moreover, with most of the datasets comprising values for block of few channels (Mirams et al., 2011; Kramer et al., 2013) and the much higher incidence of drug-induced block of particular ion channels (I_NaL, I_CaV, and I_Kr) even when drug-induced modulation of several channels are examined (Crumb et al., 2016), the result is a congregation of majority of the data in a small region of the plausible high-dimensional risk space (e.g., see Figure 2). For example, the data in Mirams et al. (2011) and Kramer et al. (2013) would fall on a single edge of the 3D EAD space in Figure 2 in the absence of values for drug-induced block on I_NaL and I_Ks in these datasets. This allows risk classification to be performed by a hyperplane with a single parameter, such as block_ICaV. Here, we utilized an additional metric, i.e., the hERG ratio, to further improve the classification performance of the direct-feature based classifiers (Figures 3, 4, 6). For the limited data currently available, risk classification using simple statistical models built on the direct features as the one presented here may suffice.

The classifiers built on derived features obtained from the in-silico models are based on certain underlying physiological phenomenon (APD, increase in calcium levels, etc.). Hence, derived features are thought to allow better extrapolation to examine drug targets other than those in the training set. However, diverse derived features from the in-silico models have been suggested as possible candidate metrics. Several features from the biophysical models, such as APD₅₀, APD₉₀, calcium level peak, and CaD₉₀ provided the best classification depending on the selected in-silico model (Mirams et al., 2011). Other derived features (EADs, TDR, change in I_CaV & I_NaL) extracted from the AP and calcium transient (Christophe, 2013, 2015; Li et al., 2017) have also been suggested as possible candidate metrics for TdP risk prediction. Rather than examining the individual features separately, a recent study performed a comprehensive feature selection among 331 metrics and determined that two metrics, APD₅₀ and diastolic Ca²⁺ in the OHR model at 1 Hz pacing, provided the best discrimination between torsadogenic and non-torsadogenic drugs (Lancaster and Sobie, 2016). The overall diversity in reported plausible candidate metrics for TdP risk classification can be attributed to different simulation protocols, drug concentrations and biophysical models. We showed that several derived features obtained from the in-silico models may track together and provide equal predictive power for risk classification when examined independent of drug EFTPC (Figure 5C). Equal predictive ability of several features makes it difficult to determine the underlying causal mechanism. In addition, the identical performance of several derived features limits the extensibility of the classifier to untrained targets, as classification results depend on the specific set of features chosen to perform the classification. For example, examination of untrained ion-channel targets using a classifier with diastolic Ca²⁺ level as the primary risk discriminating feature predicts a decrease in torsadogenic risk for increased Na⁺ − Ca²⁺ currents (Lancaster and Sobie, 2016). On the contrary, TdP risk prediction under Na⁺ − Ca²⁺ modulation using a classifier with APD₅₀ or APD₉₀ as the primary discriminating feature would predict opposite effects, with decreased Na⁺ − Ca²⁺ exchanger current being associated with decreased TdP risk. Moreover, the derived features obtained from the highly complex biophysical models did not result in improved prediction over the classfiers built on the direct features using the proposed method.

One of the primary limitations is the quality of the datasets itself. The variability in the IC₅₀ values among the several datasets can be one of the reasons for the observation of different thresholds for the hERG ratio and I_CaV block that resulted in the best discrimination between TdP+ and TdP− drugs (Figure 3). Quantification of the uncertainties in the in-vitro channel screening data and their effects on risk prediction are presented in Johnstone et al. (2016). Inconsistencies in risk definition presents another important challenge for torsadogenic risk assessment. Wiśniowska and Polak (2017) reports a comprehensive list of compounds that have been inconsistently defined as TdP+ or TdP− in different studies to develop torsadogenic risk classifiers. The different categorizations can lead to different interpretations and accuracy scores for TdP risk determination (Table 5). Standardization of torsadogenicity definition, which would allow comparison of the performance of different classifiers/features, is required. Certain steps in this direction have been started. Based on a general consensus, a working group formed under CiPA initiative picked 28 compounds and categorized each into three groups (Colatsky et al., 2016; Fermini et al., 2016) for testing/training of the classifiers. In-silico simulations of dynamic drug-channel interactions might be essential to further improve the TdP risk assessment (Li et al., 2017). Inclusion of the drug-binding parameter, in addition to the amount of block of ion-channels, resulted in 100% prediction using our approach. Sufficient I_Kr block was assumed to be necessary for TdP generation in our method. The effects of non-hERG channels are thought to enhance or mitigate the torsadogenic effects of I_Kr block. The method resulted in excellent predictive performance across several datasets that report drug-induced block of various ion channels only. However, drug-induced enhancement of ion-channel currents such as I_NaL can result in increased TdP risk in the absence of hERG block (Lacerda et al., 2008; Yang et al., 2014). The method could be further extended to examine such effects when more data are available. The present work not only provides a new method for in-vitro ion-channel screening based TdP risk classification but also highlights several important issues in regards to the use of drug-induced multi-channel blockage for torsadogenic risk prediction.

All authors are employees of IBM Research. The 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. The reviewer EP and handling Editor declared their shared affiliation.