We used the O’Hara Rudy ventricular cell model modified by Dutta et al. (2017) to determine the drug’s effect on myocardial ionic channels. The membrane potential V m ) of myocardial cells can be calculated using the formula represented by Eq. 1. dVm dt = − 1 C m I total + I stim (1) where I total is the sum of transmembrane ionic currents that consist of sodium current ( I Na ), transient outward potassium current ( I to ), late sodium current ( I NaL ), L-type calcium current ( I CaL ), sodium current through L-type calcium channel ( I CaNa ), potassium current through L-type calcium channel ( I CaK ), rapid delayed rectifier potassium current ( I Kr ), slow delayed rectifier potassium current ( I Ks ), inward rectifier potassium current ( I K 1 ), sodium-calcium exchange current ( I NaCa ), sodium-potassium ATPase current ( I NaK ), background currents ( I Nab , I Cab , I Kb ), and sarcolemma calcium pump current ( I pCa ). Meanwhile, I stim is the current induced by an external stimulus. C m is the cell membrane capacitance set at 1.0 μ F for the experiment in this study (O'Hara et al., 2011).

We utilized the model of drug effects based on the study from Mirams et al. (2011) that was inspired by the work of Hill. (1910). The inhibition effects of the drug on the ion channel could be modeled through a conduction-block formulation as expressed by Eq. 2. inhibition effect = 1 1 + IC 50 D h (2) where the IC 50 represents the concentration of 50% inhibition of ionic current, D represents the dosage of drugs, and h represents the Hill coefficient. The inhibition of the drug is assumed to affect multiple ion channels such as C aL , K 1 , K s , N a , N aL , t o , and K r or hERG. Finally, the ion channel’s conductance under the drug effect could be expressed as shown in Eq. 3. g i = g control , i 1 − inhibition effect (3) where the g i is the maximum conductance of ion channel i under drug effect and g control , i is the maximum conductance of ion channel i without drug.

In silico simulation of the drug’s effect was conducted to generate in silico features. The IC 50 and Hill coefficient were used in this study provided by Passini et al. (2017). The drug effects are simulated using various concentrations for each drug namely 1, 5, and 10. Initially, the cell’s voltage profile is simulated without adding drugs for 1,000 stimulations with a cycle length of 2,000 ms to reach the steady state condition. After that, drug effects were applied for 1,000 stimulations with the same cycle length. Following Chang et al. (2017a), the AP with the highest repolarization slope ( dVm dt repol ) within the last 250 stimulations is selected to generate in silico features. The illustration for AP and Ca profiles is shown in Figure 2. For AP beat that fully repolarized, the search of dVm dt repol is between 30%–90% repolarization; within 30% to the end of beat for AP beat that repolarises to 30% but not 90%, or between the peak of AP to the end of a cycle for AP beat that cannot repolarise by 30%. In silico features that we collected consist of the maximum rate of change of membrane potential during repolarization ( dVm dt repol ), maximum membrane potential rate ( dVm dt max ), peak of action potential ( Vm peak ), resting value of membrane potential ( Vm resting ), action potential duration from peak to 90% and 50% repolarization ( APD 90 and APD 50 ), APD tri ( APD 90 - APD 50 ), peak of calcium intracellular concentration ( Ca peak ), diastolic intracellular calcium concentration ( Ca diastole ), calcium duration from peak to 90% and 50% repolarization ( CaD 90 and CaD 50 ), Ca tri ( CaD 90 - CaD 50 ), qNet, and qInward. Following Dutta et al. (2017), the qNet was defined as the total ionic charge during AP from six ion channels as shown in Eq. 4: qNet = ∫ 0 BCL I Kr + I CaL + I to + I NaL + I Ks + I K 1 dt (4)

Furthermore, the qInward, as studied by Li et al. (2019), was defined as shown in Eq. 5: qInward = 1 2 ∫ 0 BCL I NaL , drug dt ∫ 0 BCL I NaL , control dt + ∫ 0 BCL I CaL , drug dt ∫ 0 BCL I CaL , control dt (5)

Grid search is a commonly used technique in machine learning for optimizing hyperparameters. It systematically explores all possible combinations of hyperparameter values by creating a grid configuration (Elgeldawi et al., 2021; Gressling, 2021; Belete and Huchaiah, 2022). Each combination is trained and evaluated using a validation set to assess its performance. The goal is to identify the hyperparameter values that yield the best performance. While grid search effectively finds the best hyperparameters, it becomes inefficient when dealing with high-dimensional hyperparameter spaces. As hyperparameters and their potential values increase, the number of evaluations required grows exponentially. Specifically, if there are k parameters with n distinct values, the complexity of the grid search is expected to increase at a rate of O nk . To address this issue, it is essential to carefully constrain the hyperparameter search space to improve the efficiency of grid search as an optimization approach (Elgeldawi et al., 2021). Limiting the range of possible hyperparameter values can make the grid search process more focused and computationally manageable.

This study used four classifier algorithms: KNN, XGBoost, RF, and ANN. Our study’s selection of these classifier models was driven by their specific strengths and suitability for our research objectives (Supplementary Table S4 Supplementary Material). KNN is a non-parametric algorithm characterized by its ability to operate without making assumptions about the underlying data distribution (Sha’abani et al., 2020). It is particularly suitable for situations where the data distribution is not explicitly known or may exhibit non-standard characteristics. XGBoost employs a boosting technique to improve model performance sequentially by correcting errors (Hendrawan et al., 2022; Arif Ali et al., 2023). It has robustness in handling linear and non-linear relationships, including missing data (Hendrawan et al., 2022; Arif Ali et al., 2023). RF combines multiple decision trees to improve overall prediction accuracy and reduce overfitting (Belgiu and Drăguţ, 2016). Moreover, RF is less sensitive to outliers and noise in the data (Parmar et al., 2019). ANN has the advantage of exploring complex, non-linear patterns and hierarchical features in the dataset (Hamzah and Mohamed, 2020). ANN can be adapted to various problem domains through adjustments in architecture and hyperparameters (Hamzah and Mohamed, 2020; Tuan Hoang et al., 2021). By using these individual algorithms separately, we aimed to contrast their performance and applicability, offering a holistic evaluation of their suitability for predicting drug-induced TdP risk.

Figure 3A presents the algorithm of the k-nearest neighbor (KNN) classifier algorithm employed in this study. In this approach, the training data underwent a projection into a multidimensional space, where each dimension denoted the in silico features obtained from the training data (Uddin et al., 2022; Fuadah et al., 2023). The training process encompassed the storage of feature vectors and associated labels. Meanwhile, during the prediction phase, the unlabeled testing data were labeled based on their proximity to the k nearest neighbors. Distances between feature vector positions in the training and testing data were computed using distance metrics within the multidimensional space, such as Euclidean, Chebyshev, and Minkowski. The prediction of the drug’s TdP risk is accomplished through majority voting based on the labels of the k-nearest neighbors. The optimization of the KNN model involved hyperparameter tuning utilizing the grid search method. The grid search method facilitated the selection of the best parameter values and the optimal k value from a range of options k = 1 , 3 , 5 , 7 , . . . , 31 , as well as the appropriate distance metric, including Euclidean, Minkowski, and Chebyshev, throughout the optimization process.

The Extreme Gradient Boosting (XGBoost) classification algorithm is an enhanced method based on gradient-boosting decision trees, which efficiently constructs boosted trees and operates in parallel (Chen and Guestrin, 2016; Ibrahem Ahmed Osman et al., 2021; Montomoli et al., 2021; Tarwidi et al., 2023). Figure 3B illustrates a topology representation of the XGBoost classification process, where K represents the number of decision trees and f k X i denotes the input function in the k th decision tree.

During training, the model continuously calculates node losses to identify leaf nodes with the most significant losses. XGBoost adds new decision trees by iteratively splitting input features. The objective of adding a new decision tree in XGBoost is to learn a new function, f k X , θ k , which complements the previous predictions. Once the training is completed and K decision trees are generated, each feature prediction sample corresponds to a leaf node in each decision tree, and each leaf node has an associated score. The scores from each tree are summed to obtain the final predicted value for that sample.

In this study, the XGBoost classification training model employs an ensemble of (50,100,150, and 200) decision trees. The complexity of the model increases with a higher number of decision trees. We set the options’ max depth parameter range (3, 5, 7, and 11). Additionally, the learning rate was also evaluated from 0.0001 to 0.1. The grid search method determined all the optimal tuning parameters used in this study.

The structure of the random forest (RF) classifier is shown in Figure 3C. RF comprises a group of decision tree classifier algorithms, which offer superior performance compared to using a single decision tree (Stavropoulos et al., 2020; Suhail et al., 2020; Xia, 2020). Random Forests combines two well-known classification tree approaches: boosting and bagging. It introduces an additional layer of randomness to the bagging technique. Both methods have distinct characteristics: boosting relies on the previous trees, assigning extra weight to misclassified points by earlier predictors, and making predictions based on weighted voting. On the other hand, bagging constructs each tree independently using a bootstrap sample of the dataset and makes predictions based on a simple majority vote.

RF incorporates two new strategies: Firstly, each tree is built using a different bootstrap sample of the data. Secondly, the splitting is performed at each node in the tree using the best predictor from a randomly selected subset of predictors rather than considering all variables as in standard trees. As a result, RF significantly modifies and improves upon the bagging approach by creating a diverse collection of uncorrelated trees and averaging their predictions.

In the classification process, all trees contribute by casting votes for their respective classes, and RF assigns the input to the class with the majority votes. The grid search technique identifies the optimal number of trees (50, 100, 150, and 200) and the best criterion (gini and entropy) that yields the highest performance outcome for the RF classifier.

An artificial neural network (ANN) is a fully connected architecture composed of three layers: input, hidden, and output layers (Shanbehzadeh et al., 2022; Pantic et al., 2023), as shown in Figure 3D. The input layer is responsible for receiving data from external sources. In this study, the input to the ANN architecture consisted of 14 in silico features. The hidden layers are responsible for processing the input from the preceding layer and transmitting the computed results to the output nodes. Specifically, the ANN utilized in this study incorporated one hidden layer comprising 14 nodes. The parametric rectified linear unit (PReLU) activation function was employed in the hidden layers, and a sigmoid activation function was utilized in the output layer to predict the TdP risk of drugs.

A grid search technique was employed to optimize the performance of the ANN model. The grid search aimed to identify the optimal choice of optimizer among Adam, Nadam, SGD, and RMSprop optimizers. Additionally, the grid search determined the optimal learning rate from 0.0001 to 0.1, yielding the highest performance for the ANN architecture.

Explainable AI, particularly in the context of machine learning, plays a crucial role in understanding the underlying factors driving predictions. In this study, we leveraged SHAP values to assess features’ importance in machine learning predictions. This approach is based on the concept of Shapley values from game theory, which was initially used to allocate rewards among players in a cooperative game (Lundberg et al., 2020). In the context of model interpretation, by calculating the SHAP values for each input feature, we gained insights into the contribution of individual features to the overall prediction.

In calculating SHAP values, the procedure initiates by establishing a baseline prediction, which is frequently determined by utilizing the model’s mean prediction over the entire dataset (Štrumbelj and Kononenko, 2014). The process involves systematically examining the impact of each feature by comparing the model’s prediction when including a particular feature and when excluding a particular feature. This difference reveals the extent to which a feature contributes to the prediction. Shapley values assign a credit to each feature based on its individual and collective impact on the prediction, ensuring that the contributions sum up correctly. Mathematically, the SHAP value ( ϕ ) for feature i on instance x is expressed as Eq. 6. ∅ i x = ∑ S ⊆ N \ i S ! N − S − 1 ! N ! f x s ∪ i − f x s (6)

Where N is the set of features, S is a subset of N excluding feature i , x s is the instance x with only the features included in S set. Meanwhile, f is the model’s prediction function.

In measuring evaluation metrics including accuracy, sensitivity, and specificity, we have to measure the true positive (TP), the true negative (TN), the false positive (FP), and the false negative (FN). TP represents a situation in which the model correctly predicts the high TdP risk as a high TdP risk. TN represents a situation in which no TdP risk is predicted as no TdP risk (Sharma et al., 2022). FP is when no TdP risk is wrongly predicted as high TdP risk, while the FN is when high TdP risk is incorrectly predicted as no TdP risk. In addition, we calculated accuracies, sensitivity, and specificity using Eqs 7–9, respectively. A c c u r a c y = TP + TN TP + TN + FP + FN (7) S e n s i t i v i t y = TP TP + FN (8) S p e c i f i c i t y = TN TN + FP (9)

In addition, this research also reported the area under the curve (AUC) score to assesses the classifier’s ability in distinguish between different classes. The AUC scored obtained by measuring the area under Receiver Operating Characteristic (ROC) curve that plotted two metrics including true positive rate ( s e n s i t i v i t y ) and false positive rate ( 1 − S p e c i f i c i t y .

In predicting the drug-induced TdP risk based on an electrophysiological model including inter-individual variability, we utilized 14 in silico features generated from in silico simulation of 67 drugs effect in 1,151 healthy control individuals. We provided the train and test set manually by adjusting 42 drugs as train set and 25 drugs as test set as shown in Table 1. The drugs were already categorized according to the TdP risk, which consists of 39 drugs of high TdP risk class and 28 drugs of no TdP risk class (Passini et al., 2017).

Furthermore, we performed a correlation analysis between 14 in silico features to know which features highly correlated with one another. According to the correlation heatmap between features, as shown in Supplementary Figure S1 (Supplementary Material), the highest correlation value showed between CaD 50 and CaD 90 with a correlation value of 0.77, followed by APD 50 and APD 90 with a correlation of 0.75, and Ca diastole and Vm resting with a correlation value of 0.74. However, there are no features that have correlation values between one another greater than 0.8 that are commonly used as a threshold for feature selection (Cunningham et al., 2021; Taylor, 1990; Zampieri et al., 2008). Therefore, we used 14 in silico features for predicting drug-induced TdP risk. We applied Z-score normalization to preprocess these features before using them as input to the machine learning models (Al-Faiz et al., 2018; Raju et al., 2020).

Figure 4 shows the 14 features from in silico simulations under various drug concentrations. Some features varied mainly in a narrow region, with only a few samples filling a more comprehensive range of data. For example, from the AP features, only dVm dt max and Vm peak showed a relatively wider distribution of data than other features. Most features from AP shape had some outlier samples with much larger or smaller values than most of the data samples, making the distribution plot look narrower. Furthermore, from the calcium dynamics features, Ca tri features showed a considerably narrow data distribution with the majority of data samples distributed primarily on the range of Ca tri at 400–500 ms while some outlier samples produced values more than 1750 ms. In addition, from the ionic charge features, the qInward yielded a narrow data distribution, with the majority of the data samples located around 0.5–1.5, while the qNet resulted in a significantly wider data distribution compared to qInward.

The distribution of training and testing drugs showed changes when various drug concentrations were deployed. Some features such as APD 90 , APD 50 , CaD 90 , CaD 50 , Ca tri , and qNet yielded shifted, wider data distribution when higher drug concentrations were applied mainly for high TdP-risk drugs. Among these features, only qNet produced data distribution that shifted towards negative values, i.e., more samples produced smaller or negative qNet values when higher drug concentrations were administered. However, in contrast, the features generated under no TdP-risk drugs showed minimal changes when higher drug concentrations were applied.

This study applied five-fold cross-validation with a grid search method to train 42 drugs using several machine learning models, including KNN, XGBoost, RF, and ANN. The grid search method provided the best hyperparameter setting to generate the best model of each machine learning model. Furthermore, the best model from each machine learning model was evaluated using the unseen dataset of 25 drugs, which consisted of 15 drugs of high TdP risk and 10 drugs of No TdP risk.

Our prediction models utilized simulations that observed the effect of drugs according to drug concentration variations at 1, 5, and 10 times EFTPCmax. Table 2 shows the performance in predicting drugs that induced TdP risk in 1, 5, and 10 × EFTPCmax concentrations using machine learning models. The highest performance for several machine learning models was provided at 10 × EFTPCmax. Among the several machine learning models in this study, the ANN model provided the best prediction of the drug’s TdP risk at 10 × EFTPCmax, followed by RF, XGBoost, and KNN.

We applied a grid search method for hyperparameter tuning automatically to select the best parameter of each classifier model. The grid search method selected the Adam optimizer with a learning rate of 0.001, 1,000 epochs, and a batch size 256 as the optimal parameter configuration for the ANN classifier model. Therefore, these parameters were employed to train the ANN model. Furthermore, for the RF model, the grid search selected entropy as the best criterion, with 100 trees as the optimal parameter of the RF model. For the XGBoost model, the grid search method determined 50 trees as the optimal number of estimators with three as maximum depth and learning rate 0.0001. Meanwhile, for the KNN algorithm, the grid search approach selected the Euclidean distance with a value of k = 1 as the best parameter for the KNN algorithm.

Furthermore, we have evaluated the model performance from each fold in five-fold cross-validation using the unseen dataset. The ANN model achieved the highest prediction performance with a 95% confidence interval when evaluated on the test data by obtaining an accuracy of 0.923 (0.908–0.937), sensitivity of 0.926 (0.909–0.942), specificity of 0.921 (0.906–0.935), and AUC score of 0.964 (0.954–0.975). The RF model prediction performance on test data obtained an accuracy of 0.918 (0.917–0.919), sensitivity of 0.940 (0.939–0.941), specificity of 0.888 (0.887–0.888), and AUC score of 0.972 (0.971–0.973).In addition, the XGboost model provided the highest performance on test data at 10 × EFTPCmax with an accuracy of 0.904 (0.901–0.907), sensitivity of 0.892 (0.886–0.899), specificity of 0.909 (0.903–0.914), and AUC score of 0.903 (0.902–0.905). Meanwhile, for the KNN algorithm, the highest classification performance at 5 × EFTPCmax. The best prediction obtained accuracy of 0.876 (0.872–0.882), sensitivity of 0.873 (0.864–0.881), specificity of 0.882 (0.881–0.884), and AUC score of 0.878 (0.874–0.881).

We investigated a dataset with 14 features, including dVm dt repol , dVm dt max , Vm peak , Vm resting , APD tri , APD 90 , APD 50 , Ca peak , Ca diastole , Ca tri , CaD 90 , CaD 50 , qNet, qInward. We used four classifier models: KNN, XGBoost, RF, and ANN, and conducted feature importance analysis using SHAP values to assess each feature’s contribution to the model’s predictions. As shown in Figure 5, the importance of features can vary depending on the classifier model. We observed some similarities and differences by comparing the feature importance across the different classifier models. For instance, qInward, APD 90 , APD 50 , and dVm dt repol consistently emerged as essential features in all classifier models. Interestingly, certain features displayed varying levels of importance across the models. This disparity suggests that different classifier models emphasize different feature importance.

Initially, we examined the feature importance rankings provided by the SHAP values for each classifier model using all 14 features. Subsequently, we conducted additional experiments by removing certain features that did not exhibit significant contributions in the SHAP value analysis. We analyzed the contribution of each feature based on SHAP value for all machine learning models. We considered five different groups of features: the top 3, the top 6, the top 9, the top 11, and all 14 features as input to the classifier models. The top 3 features include qInward, dVm dt repol , and APD 50 . The top 6 features include qInward, dVm dt repol , APD 50 , APD 90 , dVm dt max , and Ca peak . The top 9 features include qInward, dVm dt repol , APD 50 , APD 90 , dVm dt max , Ca peak , Ca tri , CaD 90 , and qNet. The 11 features including qInward, d V m d t r e p o l , A P D 50 , A P D 90 , d V m d t max , C a p e a k , C a t r i , C a D 90 , qNet, C a D 90 , and A P D t r i . However, our results indicated that removing certain features did not lead to significant differences in the performance of the classifier models. The models exhibited similar predictive performance and overall accuracy despite excluding features that did not demonstrate substantial contributions in the SHAP value analysis. This finding suggests that the excluded features may not have played a crucial role in the models’ decision-making processes. Table 3 compares the performance metrics achieved by the classifier models using the complete set of 14 features and the reduced sets of features. The metrics evaluated include accuracy, sensitivity, specificity, and AUC score with a 95% confidence interval. Notably, we observed that the performance of the ANN, RF, and XGBoost classifier models remained convergence for all feature groups with no statistically significant differences. Meanwhile, the KNN classifier model provided a smaller classification performance than other machine learning models.

As shown in Table 3, the ANN model provided the highest performance using 14 In silico features and still provided a good performance when only using 3 in silico features by providing an accuracy of 0.901 (0.898–0.903), sensitivity of 0.939 (0.937–0.940), specificity of 0.873 (0.870–0.875), and AUC score of 0.967 (0.966–0.968). However, removing several features in some cases leads to high sensitivity but low specificity or high specificity but low sensitivity. Based on the performance results using 3 in silico features, the higher sensitivity value compared to specificity value indicated that the model can predict high TdP risk, but is missclassified in predicting no TdP risk as high risk. Meanwhile, using 14 in silico features provided the highest sensitivity and specificity that show the model’s good capability in predicting both high TdP risk and no TdP risk.

Furthermore, we also evaluated the performance of machine learning models using the qNet feature proposed by Li et al. (2019) and APD 50 with Ca diastole proposed by Lancaster and Sobbie. (2016) The performance of machine learning models in predicting high and no TdP risk using qNet and APD 50 with Ca diastole is reported in Table 4. In predicting high and no TdP risk, all classifier models showed promissing results using APD 50 and Ca diastole as features with the highest performance provided by ANN model with an accuracy of 0.890 (0.882–0.897), sensitivity of 0.895 (0.887–0.902), specificity of 0.886 (0.879–0.893), and AUC score of 0.951 (0.944–0.957). Nevertheless, the performance decreased when only using qNet feature. The XGBoost and RF, model showed the similar classification performance as the ANN model that obtained an accuracy of 0.762 (0.761–0.763), sensitivity of 0.761 (0.760–0.762), specificity of 0.763 (0.762–0.764), and AUC score of 0.843 (0.842–0.844). In contrast, the KNN model obtained the lowest performance compared to others classifier models. This shows that using univariate feature could not capture the complete discriminatory factors contributing to the TdP risk of drugs compared to utilizing multiple in silico features.

According to the performance results, the proposed method using 14 in silico features with machine learning models obtained the highest performance in predicting high and no TdP risk of drugs. For further analysis, since CiPA categorized the TdP risk level of the drug into three categories, we also evaluated the performance of the proposed method in predicting high, intermediate, and low TdP risk as shown in Table 5. However, the proposed method could not optimally predict the TdP risk of drugs, especially due to the imbalanced class of the dataset between the high, intermediate, and low-risk categories of drugs (Supplementary Table S3. Supplementary Material).

The classification performance of three classes TdP risk of drugs provided the highest performance using 14 in silico features as input to the classifier models. The ANN model obtained an accuracy of 0.852 (0.846–0.857), sensitivity of 0.762 (0.758–0.766), specificity of 0.888 (0.886–0.889), and AUC score of 0.911 (0.907–0.914). The low sensitivity and high specificity indicated the model misclassified high and intermediate TdP risk as no TdP risk but most of no TdP risk classified correctly as no TdP risk. Therefore, for the overall performance, the model provided high accuracy and AUC score that showed the good capabilities of the model in differentiating between classes.

Furthermore, the RF model provided an accuracy of 0.835 (0.831–0.839), sensitivity of 0.705 (0.697–0.714), specificity of 0.869 (0.866–0.872), and AUC score of 0.905 (0.904–0.906). The XGBoost model provided the classification performance with an accuracy of 0.798 (0.785–0.810), sensitivity of 0.650 (0.620–0.679), specificity of 0.841 (0.829–0.854), and AUC score of 0.883 (0.867–0.899). Meanwhile, the KNN model obtained the lowest classification performance with an accuracy of 0.762 (0.761–0.763), sensitivity of 0.607 (0.606–0.608), specificity of 0.816 (0.815–0.817), and AUC score of 0.711 (0.710–0.712). In contrast with using 14 in silico features as input to the machine learning models for multiclass classification, utilizing qNet, or APD 50 with Ca diastole as features lead to decreasing classification performance of high, intermediate, and low TdP risk of drugs. It became evident that using only one or two features was insufficient in capturing the comprehensive set of factors contributing to the TdP risk of drugs when contrasted with the utilization of all 14 in silico features.