Introduction

Torsades des Pointes is a potentially lethal ventricular arrhythmia that can be induced by many classes of drugs. These include antibiotics, antipsychotics, antihistamines, chemotherapeutics and antiarrhythmics (Yap and Camm, 2003). The majority of torsadogenic drugs block the hERG ion channel which carries the delayedrectifier potassium current (I_Kr) (Witchel, 2011). For this reason, international safety guidelines require hERG block to be assessed in living cells prior to conducting human trials (ICH Harmonised Tripartite Guideline, 2005). However, the standard hERG assay is overly sensitive. It does not accommodate multi-channel effects which render some drugs safe despite blocking hERG (Martin et al., 2004; Hoffmann and Warner, 2006). Consequently, many useful drugs are prematurely abandoned during pre-clinical trials. The safety pharmacology community are actively pursuing new in vitro and in silico assays that improve accuracy by targeting multiple ion channels (Pugsley et al., 2008; Colatsky et al., 2016).

In silico assays use computational models of cardiomyocyte electrophysiology in place of a living cell (Mirams et al., 2011; Lancaster and Sobie, 2016; Mann et al., 2016; Dutta et al., 2017; Passini et al., 2017; Ballouz et al., 2020; Llopis-Lorente et al., 2020). Drug blockade (Figure 1A) is simulated in the model (Figure 1B) by attenuating the conductivity of the relevant ion currents. The conductance parameters are labeled G by convention. The simulations can be repeated across a diverse range of cardiac phenotypes to ensure generality of the results (Britton et al., 2013; Ni et al., 2018; Gong and Sobie, 2018). The individual phenotypes are constructed by randomizing the conductance parameters to mimic natural variation in ion channel expression. The method produces a population of cardiac action potentials (Figure 1C). Selected biomarkers within the action potentials are then statistically analyzed for drug-induced changes (e.g. Hondeghem, 2005; Varshneya et al., 2021). Contemporary research is largely concerned with improving those biomarkers.

Figure 1.

Yet the main problem with the conventional approach is that it requires multitudes of computationally intensive simulations for every drug that is assessed. Pharmacology laboratories must invest heavily in specialist computing resources and expertise before they can apply the methods to their own drugs. We propose a new approach that allows drugs to be assessed without conducting drug-specific simulations. The method is initiated by simulating a diverse population of cardiomyocytes in the absence of drugs. That simulation need only be done once. The drug-free population is then used to identify the principal relationship between ionic conductances and ectopic phenotypes. We call that relationship the axis of arrhythmia because it describes the principal pathway for transforming benign cardiomyocytes into ecoptic cardiomyocytes. Thereafter, the axis serves as a yardstick for assessing the torsadogenic risk of drugs directly from their physiological signature of ion channel blockade alone.

Results

The action potentials of 100,000 randomized ventricular cardiomyocytes (Figure 2A) were simulated using a variant of the O’Hara-Rudy model (O’Hara et al., 2011) that was optimized for long-QT syndrome (Mann et al., 2016; Krogh-Madsen et al., 2017). Early afterdepolarizations were chosen as a biomarker for Torsades des Pointes. Cardiomyocytes that exhibited them were classified as ectopic (red) and those t hat d id not were classified as benign (gray).

Figure 2.

The four cardiac ion currents that we investigated (I_CaL, I_Kr, I_NaL, I_Ks) have previously been implicated in torsadogenic risk (Dutta et al., 2017; LlopisLorente et al., 2020). The myocytes were constructed by re-scaling the conductance parameters (G_CaL, G_Kr, G_NaL, G_Ks) with randomly selected multipliers that were drawn uniformly from a logarithmic scale (Figure 2B). The use of the logarithmic coordinate frame is crucial to our subsequent analysis.

The ectopic and benign phenotypes were found to be clearly segregated in parameter space (Figures 2C and 2D). Multivariate logistic regression was used to identify the linear boundary (dashed line) that best separated the two classes. We refer to it as the decision boundary following the conventions of classification statistics.

Decision boundary

The multivariate logistic regression equation,

describes the log-odds of a cardiomyocyte being ectopic where X_CaL = ln(G_CaL), X_Kr = ln(G_Kr), X_NaL = ln(G_NaL) and X_Ks = ln(G_Ks). The decision boundary is the hyperplane in four-dimensional parameter space (X_CaL, X_Kr, X_NaL, X_Ks) where p = 0.5. Figures 2C and 2D show the intersection of that hyperplane (dashed) with two-dimensional slices of parameter space. Although we illustrate the concept in two dimensions, the analysis itself is conducted in four dimensions.

Our estimates of the regression coefficients, β₀ = 6.416±.059, β_CaL = 2.509±.024, β_Kr =− 2.471± .024, β_NaL = 0.847±.012, β_Ks = −1.724±.018, were all statistically significant (p < .001). The reported confidence intervals represent 1±SE. The model itself was also significantly different from the null model, χ²(99995, n=100000) = 73900, p < .001.

The action of drugs in logarithmic coordinates

Drugs attenuate the conductivity of ion channels in a multiplicative fashion. The conductance of the drugged channel is defined as,

where G_ion is the baseline conductance of the ion channel and δ∈ [0, 1] is the fractional conductance imposed by the drug. Logarithmic coordinates transform the action of drug blockade into an additive process,

which can be expressed with vector geometry (Figure 3A). We denote the action of a drug by the vector,

Figure 3.

where α = ln(δ). The values of δ would ordinarily be obtained from patch-clamp experiments. In our case, we calculated them from the published potencies of 109 drugs collated by Llopis-Lorente et al. (2020).

Axis of arrhythmia

We define the axis of arrhythmia as a line that runs orthogonally to the decision boundary (Figure 3A). It represents the shortest path for shifting any cardiomyocyte into the ectopic regime by modifying the conductances of its ion channels. The basis vector of the axis,

is defined by the coefficients of the regression equation.

Drug risk metric

The arrhythmogenic component of a drug is obtained by projecting the action of the drug onto the axis of arrhythmia. The length of the projection is our metric of drug risk. Specifically,

where A is the action of the drug, B is the basis vector of the axis of arrhythmia, and A·B = Σi α_iβ_i is the dot product. The metric is normalized to the euclidean length of B, which is denoted ║B║. For our estimates of the regression coefficients,

Risk scores

According to our risk metric, drugs that shift the electrophysiology towards the ectopic region have positive scores. Conversely, drugs that shift the electrophysiology away from ectopy have negative scores. Drugs that do neither have scores near zero.

The a priori probability of ectopy

The probability of ectopy along the axis of arrhythmia is,

where x = β₀ + β_CaLX_CaL + β_KrX_Kr + β_NaLX_NaL + β_KsX_Ks. Eq. (3) is obtained by rearranging Eq. (1). It describes the cross-sectional profile of the decision boundary (Figure 3B). As such, it represents the a priori probability of a cardiomyocyte being ectopic, on average, based on its proximity to the decision boundary. The shallow slope of the profile is indicative of the statistical uncertainty in fitting a linear boundary to the data points (Figures 2C,D).

Susceptibility in the natural population

A drug alters the electrophysiology of all cardiomyocytes to the same extent but only some cardiomyocytes become ectopic. The most susceptible are those that are closest to the decision boundary. We calculated the proportion of the natural population that would be susceptible to a given drug by analyzing how the drug shifts the population density with respect to the decision boundary.

Following the methods of Sobie (2009), Sadrieh et al. (2013), Morotti and Grandi (2017), Gong and Sobie (2018) we assumed that ion channel conductances varied independently in nature and were log-normally distributed (μ = − 0.112, σ = 0.472). By definition, a lognormal distribution maps onto a normal distribution in logarithmic coordinates. The natural population in our parameter space is therefore a symmetric multivariate Gaussian density function (Figure 4A). In the absence of drugs, the natural population density in four dimensions is centered at the point O = {μ, μ, μ, μ}.

Figure 4.

The proportion of myocytes that become ectopic depends on how far the population is shifted along the axis of arrhythmia according to the risk metric (Figure 4B). The proportion is calculated as the product of the a priori probability of the ectopy (Eq. 3) and the density of the drugged population (Eq. 5; Methods). For the case of 10x therapeutic dose of Ibutilide — a potent blocker of I_Kr — the proportion of the natural population that is susceptible to drug-induced Torsades is 41.5% (Figure 4). Since the size of the susceptible population is monotonically related to the drug’s risk score, the torsadogenic risk can be described using either terminology.

Validation against known drugs

We used the Hill equation to reconstruct the drugresponse profiles of G_Kr, G_CaL, G_Ks and G_NaL from the half-maximal inhibitory concentrations (IC₅₀) of 109 drugs reported by Llopis-Lorente et al., 2020. The results for each drug are included in the Supplementary Information. Llopis-Lorente et al. (2020) labeled the clinical risks according to the Credible Meds list of QT drugs (Woosley et al., 2019). In that labeling scheme, Class 1 drugs carry a known torsadogenic risk; Class 2 drugs carry a possible risk; Class 3 drugs carry a risk but only in conjunction with other factors; Class 4 drugs have no evidence of risk at all.

Cases of Ajmaline and Linezolid

Ajmaline is an anti-arrhythmic drug that primarily blocks the fast sodium current (Kiesecker et al., 2004). Clinicians must administer Ajmaline with care because it is also a potent blocker of I_Kr (Figure 5A, left) with a known risk of Torsades (Class 1). In comparison, Linezolid is an antibacterial agent that has no clinical evidence of Torsades (Class 4) even though it too blocks I_Kr. Albeit less than it blocks I_CaL (Figure 5A, right).

Figure 5.

Indeed, the two drugs have nearly opposite effects on G_CaL and G_Kr, as can be seen in the case of 25x therapeutic dose (Figure 5B). Using vector notation, the blocking action of 25x dose of Ajmaline is,

and the blocking action of Linezolid is,

The two drugs shift the electrophysiology in opposite directions along the axis of arrhythmia (Figure 5C). The corresponding risk scores are +1.16 for Ajmaline and − 1.18 for Linezolid. The signs of the scores suggest that Ajmaline is pro-arrhythmic whereas Linezolid is not, which agrees with the clinical risk labels. Indeed, the sizable negative score for Linezolid suggests that it may actually have anti-arrhythmic properties.

The effect of both drugs on the natural population is shown in Figure 5D. Ajmaline shifts the population density by 1.16 units towards the ectopic regime, which suggests that 26% of the recipients would be susceptible to Torsades. Conversely, Linezolid shifts the population 1.18 units away from the ectopic regime. So only 0.0095% of those who received Linezolid would be susceptible. A substantial drop from the baseline rate of 0.93%.

Drug action varies with dose

The direction of a drug’s line of action is liable to vary with dosage because of differential binding to each ion channel. For example, Ajmaline disproportionately blocks more I_Kr than I_CaL with increasing dosage (Figure 6A). Hence the pathway of Ajmaline follows a gentle arc in parameter space rather than a straight line (Figure 6B). The effect is less pronounced for Propafenone which recruits the ion channels more uniformly (Figure 6C). Nonetheless, its path is still not strictly linear (Figure 6D). Linezolid follows a similar pattern to Ajmaline but in the opposite direction (Figures 6E and 6F). In general, the pathways of all 109 drugs (light gray traces) exhibited a variety of curvatures that influence their projection onto the axis of arrhythmia. Consequently the relative scores of some drugs change subtly with dosage.

Figure 6.

Testing the risk metric

The metric was tested by scoring all 109 drugs over a range of doses and comparing the results to the clinical risk labels from Credible Meds. The clinical labels were lumped into two categories for this purpose: UNSAFE (class 1 and 2) versus SAFE (class 3 and 4). Drugs that scored above a given threshold (risk > θ) were predicted to be UNSAFE and those that scored below the threshold (risk < θ) were predicted to be SAFE. The threshold was optimized for each dosage. For drugs assessed at 25x dose, a classification accuracy of 90.8% was achieved using a scoring threshold of θ = 0.195 (Figure 7).

Figure 7.

The procedure was repeated for doses ranging from 1x to 32x. The classification accuracy for all dosage levels was found to lay between 88.1% and 90.8% (Figure 8A). The differences were primarily due to classification errors in a small number of borderline cases. So we refrain from nominating any one dosage as being optimal. The conventional hERG assay, in comparison, has an accuracy of 78.9% with the same dataset (LlopisLorente et al., 2020).

Figure 8.

From a safety perspective, the trade-off between false negatives and false positives can be tuned by adjusting the scoring threshold, θ. This is illustrated by the Receiver Operating Characteristic (ROC) curves (Figures 8B and 8C). In our case, there is little difference between the ROC curves for drugs assessed at 5x dose versus 25x dose. The area under the ROC curves (AU-ROC) are nearly identical at 91.3% and 91.5% respectively. In comparison, the conventional hERG assay has an AUROC of 77 ± 7% (Kramer et al., 2013).

Discussion

In this study we have proposed a new metric of torsadogenic risk that is based on the axis of arrhythmia. The major benefit of the metric is that it can be applied to novel drugs without conducting any new computer simulations. The drug response profiles of four ion currents (δ_CaL, δ_Kr, δ_NaL, δ_Ks) are all that is needed to calculate the torsadogenic risk of the drug. The ion currents can be measured using standard patch clamp techniques and the risk metric can be calculated with pen and paper. This simplicity removes a technological hurdle to the adoption of computational assays in safety pharmacology.

Identifying the axis of arrhythmia

All of the simulations in the present study were conducted in the absence of drugs. They were only needed to identify the axis of arrhythmia and do not need to be repeated when applying the metric. The axis encapsulates the linear relationship between the four ion channels and the onset of early after-depolarizations. It corresponds to the most potent combination of proarrhythmic drug block that is theoretically possible. The assumption of linearity is crucial for averaging the risk across all individuals. So even though a nonlinear decision boundary might fit the data better, it would not be helpful because the results would be patient specific.

Susceptibility to drug-induced arrhythmia

Our method can also predict the proportion of the population that would be susceptible to a drug without explicitly simulating it. The analysis is possible because drug blockade and natural diversity both operate on the same properties of ion channels. The two biological processes are therefore mathematically interchangeable. As such, the distribution of the drugged population can be inferred from the drug-free population by shifting it according to the drug’s risk score. The size of the susceptible population is a function of the risk score. So the torsadogenic risk can be reported either as a raw score or as a percentage of the population at risk.

Comparison to conventional approaches

Conventional in silico safety assays are designed to apply the drug directly to simulated cardiomyocytes and then use biomarkers in the action potential to predict the torsadogenic risk. See Grandi et al. (2018) for a review. The predictions are optimized for established drugs using a training procedure that has been formalized by the Comprehensive In Vitro Proarrhythmia Assay (CiPA) initiative (Colatsky et al., 2016). The CiPA steering committee (Li et al., 2019) recommends the CiPAORdv1.0 variant of the O’Hara-Rudy model which incorporates the kinetics of drugs binding with the hERG channel (Li et al., 2017). The recommended biomarker is qNet (Dutta et al., 2017) which measures the net charge of six major ion currents (I_Kr, I_CaL, I_NaL, I_to, I_Ks, I_K1). The clinical risks of the CiPA drugs (n=28) are labeled low, intermediate or high. Low versus intermediate-or-high risk drugs were predicted with 84–95% accuracy using manual patch clamp techniques; or 93-100% accuracy using automated patch clamping (Li et al., 2019). Whereas low-or-intermediate versus high risk drugs were predicted to 92–100% accuracy using manual patch clamp; or 88-98% using automated patch clamp (Li et al., 2019). The measurement of the drug potencies is a source of considerable variability, which is exacerbated by the small number of test drugs (n=16). The accuracy ranges reported by Li et al. (2019) appear to outperform our results, although they cannot be compared directly because the risk labels are stratified very differently.

Passini et al. (2017) used a larger dataset (n=62) to obtain 89% accuracy at predicting class 1 versus class 2–4 torsadogenic risk labels from Credible Meds. The cardiomyocytes were simulated using the baseline O’HaraRudy model without dynamic drug-binding kinetics. The accuracy is within the range of our results, albeit using a slightly different risk stratification scheme. Passini et al. (2017) scored the torsadogenic risk independently of dosage by averaging number of early afterdepolarizations, weighted by the concentration of the drug. However we believe it is better to quantify the risk as a function of dosage, since even the most lethal drug is safe at zero dose.

The present study uses the same dataset (n=109) as Llopis-Lorente et al. (2020) who investigated a suite of biomarkers using a recalibrated O’Hara-Rudy model which also had revised gating kinetics for I_Na. Their best performing biomarker was qNet, which predicted class 1–2 versus class 3–4 torsadogenic risk with 92.7% accuracy. It exceeds our best result by approximately 2% on the same risk stratification scheme. LlopisLorente et al. (2020) increased the overall performance to 94.5% by combining the best biomarkers into a decision tree. Interestingly, early after-depolarizations proved to be their worst biomarker with only 78.9% accuracy. Based on our analysis, it is likely they undersampled the early after-depolarizations by constraining their cardiomyocytes to normal physiological limits.

Comparison to MICE models

MICE models (Kramer et al., 2013) use logistic regression to predict the torsadogenic risk of a drug directly from the half maximal inhibitory concentrations (IC₅₀) for hERG, Cav1.2 and Nav1.5. They are derived using statistical model fitting techniques. Kramer et al. (2013) trained six candidate models on a dataset of 55 drugs. Their best model was found to predict the torsadogenic labels of the drugs with 90.9% accuracy using only the difference between hERG and Cav1.2. Specifically,

where H= − ln(IC₅₀) for hERG, and C = − ln(IC₅₀) for Cav1.2. The form of the predictor variables are strikingly similar to our own except that we use the fractional conductance (δ) instead of IC₅₀. In our case, the predictors were derived from biophysical principals of drug block rather than through statistical model fitting. The accuracy of the two models are nearly identical, albeit on different datasets.

Comparison to Bnet

Bnet (Mistry, 2018) is a simple linear model that predicts the torsadogenic risk of a drug directly from the net blockade of inward and outward ion currents. Notably,

where R_i and D_j represent the percentage block of the repolarizing currents (I_Kr, I_Ks, I_to) and the depolarizing currents (I_CaL, I_Na, I_NaL, I_K1) respectively. Percentage block is akin to (1 −δ) in our model but without the logarithmic transform that we use to analyze attenuation as an additive process. Bnet predicts the clinical risk labels of the CiPA validation drugs just as well as the CiPAORdv1.0 model when adjusted for drug binding kinetics (Li et al., 2019; Mistry, 2019; Han et al., 2019). This has opened a debate between model complexity and biophysical realism in which proponents of biophysical models advocate their explanatory benefits (Lancaster and Sobie, 2017), whereas proponents of simple models advocate their predictive power without the computational expense (Mistry et al., 2015; Mistry, 2017).

Conclusion

Our approach resolves the debate between model complexity and biophysical realism by combining both approaches into the same enterprise. The method is initiated using biophysical models to explore the relationship between ion channels and torsadogenic risk — as it is best understood by theory. The findings are then condensed into a simple linear model that can be applied without recapitulating the biophysical details. That is not to say the linear model is mere statistical correlation. It still represents the biophysical action of multichannel drug blockade but only as far as it relates to early after-depolarizations. Our approach thus represents a convergence of biophysical and simple models. One that provides the benefits of biophysical modeling without the computational burden. We believe it will accelerate the adoption of computational assays in safety pharmacology and ultimately reduce the burden of animal testing.

Methods

The action potentials of human endocardial ventricular cardiomyocytes were simulated using a variant of the ORD11 model (O’Hara et al., 2011) in which the maximal conductances of G_Ks, G_Kr, G_CaL, G_NaL, P_NaCa and P_NaK were re-scaled (Table 1) to better reproduce the clinical phenotypes of long QT syndrome (Mann et al., 2016; Krogh-Madsen et al., 2017). The source code for the ORD11 model (O’Hara et al., 2011) was adapted to run in the Brain Dynamics Toolbox (Heitmann et al., 2018; Heitmann and Breakspear, 2022). The adapted source code is available from Heitmann (2023). The differential equations were integrated forward in time using the matlab ode15s solver with error tolerances AbsTol=1e-3 and RelTol=1e-6. The model was paced at 1 Hz with a stimulus of -70 mV and duration of 0.5 ms. All simulations were equilibrated for at least 1000 beats prior to analysis.

Table 1.

Classification of ectopic phenotypes

Blocks of four successive beats were analyzed to accommodate alternans. Cells were classified a s ectopic if any of those four beats contained an early afterdepolarization — as defined by any secondary peak that rose above -50 mV and was separated from other peaks by at least 100 ms. Cells that did not exhibit early afterdepolarizations were classified as benign.

Parameter domain

We chose the domain of parameter space (G e^±3) using trial and error. The domain must be large enough to cover the ectopic region but not so large as to be unduly influenced by biological e xtremes. Our chosen domain spans conductances between 0.05 and 20 times their baseline value.

Natural population density

Cell-to-cell variability was mimicked by scaling the conductances (G_CaL, G_Kr, G_NaL, G_Ks) by a random multiplier that was drawn from a log-normal distribution (Sobie, 2009). The parameters of the distribution (μ=0.112, σ = 0.472) were chosen to give the multipliers a mean of 1 and standard deviation of 0.5. The spread was based on our experience with previous simulations (Sadrieh et al., 2013, 2014; Ballouz et al., 2020; TeBay et al., 2022). The four ion channel species were assumed to vary independently.

By definition, the log-normal distribution maps onto the normal distribution under the logarithmic transform. The natural population was therefore represented in logarithmic parameter space by a symmetric multivariate normal distribution. Specifically,

Where

is the univariate normal distribution.

Joint probability

Rotational symmetry allowed the multivariate distribution to be projected onto the axis of arrhythmia as a univariate distribution. The probability of a cardiomyocyte in the natural population being ectopic is,

where p(x) is the a priori probability of ectopy for that phenotype (Eq. 3) and f (x) is the proportion of cells in the population with that phenotype (Eq. 4). Drugs serve to shift f (x) along the axis of arrhythmia. The size of that shift is defined by the risk metric (Eq. 2).

Drug dataset

The drug response curves were reconstructed using the Hill equation,

where δ∈ [0, 1] represents the fractional conductance of the ion channel. The concentration of the drug (C) was normalized to the effective free therapeutic plasma concentration (EFTPC). The Hill coefficient was fixed at h = 1 on the assumption that single drug molecules bound to the ion channels independently. The half-maximal inhibitory concentrations (IC₅₀) were taken from Supplementary Table S2 of Llopis-Lorente et al. (2020). That data was curated from multiple public databases and scientific publications to ameliorate the effects of crosssite and cross-platform variability (Kramer et al., 2020). The dataset uses the median IC₅₀ for cases where multiple values were reported (Llopis-Lorente et al., 2020). Similarly, it uses the highest therapeutic dosage for cases where multiple EFTPCs were reported (LlopisLorente et al., 2020). Ion channels with missing IC₅₀ values were assumed not to be blocked by the drug (δ = 1). The reconstructed response curves are given in Supplementary Dataset S1. The clinical risk labels for the drugs were transcribed from Table 1 of LlopisLorente et al. (2020).

Acknowledgements

This study was supported by Australian NHMRC grants app1182032, app1182623 and NSW Health cardiovascular capacity building grant. SH was also funded by The Medical Advances Without Animals Trust (MAWA). The aim of MAWA is to advance medical science and improve human health and therapeutic interventions without the use of animals or animal products. This research includes computations using the computational cluster Katana supported by Research Technology Services at UNSW Sydney.

Supplementary Videos

Supplementary Video SV1

Supplementary Information

Supplementary Dataset S1

Drugs with Class 1 Torsadogenic risk

Drugs with Class 2 Torsadogenic risk

Drugs with Class 3 Torsadogenic risk

Drugs with Class 4 Torsadogenic risk