Classifying Drugs by their Arrhythmogenic Risk Using Machine Learning.

All medications have adverse effects. Among the most serious of these are cardiac arrhythmias. Current paradigms for drug safety evaluation are costly, lengthy, conservative, and impede efficient drug development. Here, we combine multiscale experiment and simulation, high-performance computing, and machine learning to create a risk estimator to stratify new and existing drugs according to their proarrhythmic potential. We capitalize on recent developments in machine learning and integrate information across 10 orders of magnitude in space and time to provide a holistic picture of the effects of drugs, either individually or in combination with other drugs. We show, both experimentally and computationally, that drug-induced arrhythmias are dominated by the interplay between two currents with opposing effects: the rapid delayed rectifier potassium current and the L-type calcium current. Using Gaussian process classification, we create a classifier that stratifies drugs into safe and arrhythmic domains for any combinations of these two currents. We demonstrate that our classifier correctly identifies the risk categories of 22 common drugs exclusively on the basis of their concentrations at 50% current block. Our new risk assessment tool explains under which conditions blocking the L-type calcium current can delay or even entirely suppress arrhythmogenic events. Using machine learning in drug safety evaluation can provide a more accurate and comprehensive mechanistic assessment of the proarrhythmic potential of new drugs. Our study paves the way toward establishing science-based criteria to accelerate drug development, design safer drugs, and reduce heart rhythm disorders.

Significance

Drugs can have serious side effects and cause cardiac arrhythmias. Drug safety evaluation is expensive and lengthy. Here, we establish an easy-to-use diagram to stratify the risk of new and existing drugs. We use machine learning to integrate knowledge across 10 orders of magnitude in space and time and provide a holistic picture of the effect of drugs. Our approach identifies a pair of agonist-antagonist ionic currents that dominate arrhythmogenic events. For any combinations of these two currents, we create a single classifier that stratifies safe and arrhythmic regimes. We demonstrate that our classifier correctly identifies the risk categories of 22 common drugs. Our study could help accelerate drug development, design safer drugs, and reduce rhythm disorders.

Introduction

Developing a new drug is an expensive and lengthy process. The estimated average cost to design and approve a new drug is $2.5 billion (1), and the time to market from the initial discovery into the pharmacy is at least 10 years (2). Many drugs, not just cardiac drugs, interact with specific ion channels in the heart and can induce serious rhythm disorders (3). Indeed, the major focus of FDA toxicity testing is the proarrhythmic potential of a drug, as determined by its effect on repolarization. Specifically, the approval of a new drug requires assessing its impact on the rapid component of the delayed rectifier potassium current in single-cell experiments (4) and on the duration of ventricular activity in animal models and in healthy human volunteers (5). Unfortunately, the high cost and long time to test new compounds acts as an impediment to the discovery of new drugs (6). Further, the limited window provided by these criteria onto proarrhythmic potential generates false positives while at the same time preventing many potentially useful drugs from ever reaching the market (7). A combined approach of machine learning and multiscale modeling could significantly accelerate the early stages of drug development, guide the design of safe drugs, and help reduce drug-induced rhythm disorders (8).

Torsades de pointes is a serious side effect of many drugs

All pharmacological agents have the potential to impact cardiac repolarization and, with it, the QT interval. The most serious manifestation of both genetic and drug-induced long QT intervals is torsades de pointes, a ventricular arrhythmia characterized by rapid, irregular patterns in the electrocardiogram (9). Most episodes of torsades de pointes begin spontaneously and revert to normal sinus rhythm within a few seconds, but some persist, degenerate into ventricular fibrillation, and lead to sudden cardiac death, even in patients with structurally normal hearts (10). In the United States, more than 350,000 sudden cardiac deaths occur each year, but the true incidence of torsades de pointes is largely unknown (11). Predicting this potentially fatal heart rhythm is challenging given the complex interplay between genetic predisposition and medications, both prescription and over the counter. Increasing evidence suggests that early afterdepolarizations play a critical role in generating of torsades de pointes (12). Early afterdepolarizations are oscillations during the repolarization phase of the cellular action potential that result from a reduced outward current, an increased inward current, or both (13). The theory of nonlinear dynamics can help explain the ionic basis of early afterdepolarizations (14); yet, it remains unclear how early afterdepolarizations translate into three-dimensional arrhythmias. A better quantitative understanding of the relevant ionic currents would significantly reduce the design space and accelerate drug screening in the early stages of drug development.

Machine learning could help accelerate drug development

Leading pharmaceutical companies have long recognized the potential of machine learning, especially during the early stages of drug development: on the protein and cellular levels, machine learning can help identify efficient drug targets, confirm hits, optimize leads, and explain the molecular basis of therapeutic activity (15). On the tissue and organ levels, machine learning can guide pharmacological profiling and predict how a drug that was designed in the lab will affect an entire organ (16). Although using machine learning in the early stages of drug design, target selection, and high throughput screening is almost standard today, the potential of machine learning in the later stages of drug development, toxicity screening, and risk stratification has not been recognized to its full extent (17). A promising application of machine learning in the context of cardiotoxicity is to combine several experimentally measured and computationally simulated features into a unifying classifier for torsadogenic risk assessment (18). A recent study demonstrated that a machine-learning classifier that combines cellular action potentials and intracellular calcium waveforms provides better torsadogenic risk prediction than one focused on potassium channel block alone (19). Although there is a general agreement between clinical researchers, pharmaceutical companies, and regulatory agencies that computational tools should play a more central role in the proarrhythmic risk assessment of new drugs (20), current efforts focus exclusively on classifiers at the single-cell level and ignore ventricular heterogeneity and the interaction of different cell types across the entire heart (21). Previous computational studies have successfully focused on analyzing the QT interval, which, however, is known to have a low specificity to discriminate proarrhythmic risk (22). We have recently proposed a novel, exposure-response simulator that allows us to quickly and reliably visualize how different drugs—either individually or in combination—modulate ion channel dynamics, cellular electrophysiology, and electrocardiogram recordings across 10 orders of magnitude in space and time (23). Combining this simulator with machine-learning techniques (24) would allow us to seamlessly integrate experimental and computational data from the protein, cellular, tissue, and organ scales to assess cardiac toxicity during pharmacological profiling (25).

Fig. 1 illustrates how we use machine learning to combine computational (top) and experimental (bottom) tools and technologies at the single-cell (left) and whole-heart (right) levels. First, we probe how different ion channels modulate early afterdepolarizations on the single-cell level. Using a hybrid computational and experimental approach, we identify the two most relevant channels and systematically screen the two-channel parameter space to quantify the critical blockade that initiates torsades de pointes. Then, we use high-performance computing and machine learning to identify the classification boundary between the arrhythmic and nonarrhythmic domains in this space. We validate our approach using computational and experimental electrocardiograms from whole-heart simulations and isolated Langendorff perfused hearts. Finally, we demonstrate the potential of our classifier by risk stratifying 23 common drugs and comparing the result against the reported risk categories from the literature.

Figure 1

Hybrid computational-experimental approach to quickly and reliably characterize the proarrhythmic potential of existing and new drugs. We characterize calcium transients in ventricular cardiomyocytes in response to drugs both computationally (top) and experimentally (bottom) and identify the ion channels that most likely generate early afterdepolarizations (left). We then screen the concentration space of the two most relevant channels and identify the classification boundary between the arrhythmic and nonarrhythmic domains using high-performance computing and machine learning (center). We validate our approach using electrocardiograms, both computationally and experimentally, in whole-heart simulations and isolated Langendorff perfused hearts (right). We demonstrate the potential of our new classifier by risk stratifying 23 common drugs and comparing the result against the reported risk categories of these compounds. To see this figure in color, go online.

Materials and Methods

All studies were approved by the Stanford Administrative Panel on Laboratory Animal Care and conform to the Guide for the Care and Use of Laboratory Animals published by the National Institutes of Health.

Simulating action potentials in ventricular cardiomyocytes

We modeled the temporal evolution of the transmembrane potential ϕ using an ordinary differential equation,where Cm is the membrane capacitance and Iion(ϕ, ) is the ionic current, which we represented as a function of the transmembrane potential ϕ and a set of state variables (26). The state variables obey ordinary differential equations,  = (ϕ, ), as functions of the transmembrane potential ϕ and their current values (27). For our single-cell simulations, we used ventricular cardiomyocytes with 15 ionic currents and 39 state variables (28),with a minor modification (29) of the fast sodium current INaP (30). We parameterized the model for human midwall cells (28) and modeled the effect of drugs by selectively blocking the relevant ionic currents Iion (22). For a desired concentration C, for each current i, we calculate the fractional block βi using a Hill-type model parameterized with data from patch-clamp electrophysiology (20,31) and scale the ionic current Ii by this fractional block (23),

We studied the relative importance of seven ion channels—ICaL, IK1, IKr, KKs, INaL, INaP, and Ito—on inducing early afterdepolarizations. To achieve a steady state, we paced the cells for 600 cycles at a frequency of 1 Hz. We defined the presence of early afterdepolarizations as the occurrence of a change in potential greater than 0.1 mV/ms between the 50 and 1000 ms of the last two recorded cycles (21). We used a Latin hypercube design to perform 500 simulations and systematically varied the block of the seven ion channels between 0 and 95%. Then, we labeled the results depending on the presence or absence of early afterdepolarizations. We fitted a logistic regression and computed the marginal effects, which correspond to the derivative of the output of the regression with respect to the ion channel block (32,33). We normalized the results by the maximal value.

Simulating electrocardiograms in human hearts

To pass information across the scales, we created an ultra-high-resolution finite element model of the human heart (23) that represents individual ion channel dynamics through local ordinary differential equations at the integration point level and action potential propagation through global partial differential equations at the node point level (34). The basis of this model is the classical monodomain model that characterizes the spatiotemporal evolution of the transmembrane potential ϕ through the following partial differential equation,

In addition to the local source term Iion/Cm from Eq. 1, the transmembrane potential depends on the global flux term , where is the conductivity tensor that accounts for a fast signal propagation of D|| parallel to the fiber direction and a slow signal propagation of perpendicular to it (26),

We used the O’Hara-Rudy model (28) from Eq. 2 for all ventricular cells and the Stewart model (35) for all Purkinje cells. We discretized the monodomain Eq. 4 in time using finite differences and in space using finite elements (26) and introduced the transmembrane potential as a degree of freedom at the node point level and all state variables as local degrees of freedom at the integration point level (27). We solved the resulting system of equations using the finite element software package Abaqus (36) with an explicit time integration scheme. We discretized our simulation window of five healthy heart beats in time using 1.0 M equidistant time steps of Δt = 0.005 ms. We created our human heart model from magnetic resonance images of a healthy, 21-year-old, 50th percentile U.S. male (37), which we discretized in space using 6.9 M regular trilinear hexagonal elements with a constant edge length of h = 0.3 mm. This resulted in 7.5 M global degrees of freedom and 0.3 G local internal variables (34).

Using machine-learning tools to sample the parameter space

To quickly and efficiently sample the parameter space for a wide range of conditions and a wide variety of drugs we combined our computational models with machine-learning techniques (18,25). To characterize ventricular fibrillation, we performed n = 40 human heart simulations, which we classified into arrhythmic and nonarrhythmic depending on whether the model showed irregular or regular activations. To identify the classification boundary that divides the arrhythmic and nonarrhythmic domains, we used a Gaussian process classifier (38), which predicts an output distribution. From this distribution, we determined the point of maximal entropy in the parameter space, where we placed the next sample (39). We generated the first n = 10 samples from a Latin hypercube design and adaptively identified the locations of the remaining n = 30 samples. Our results suggest that n = 40 simulations are sufficient to reliably identify the classification boundary.

Classifying drugs into risk categories

We classified 23 drugs into high and low risk, based on our proarrhythmic risk estimator in Fig. 4, and validated our approach against the known risk classification of these drugs. To select the compounds, we began with a list of 31 drugs (20) for which the concentration block is thoroughly characterized. From these 31 drugs, we only considered those for which 70% or more of the published studies agreed on their risk classification (40,41) and did not consider the remaining eight controversial drugs. Table 1 summarizes the IC50 values used to compute the degree of blockade of the L-type calcium current ICaL and the rapid delayed rectifier potassium current IKr (20).

Figure 4

Proarrhythmic risk classification. Screening the parameter space of rapid delayed rectifier potassium current IKr and the L-type calcium current ICaL block reveals the classification boundary beyond which arrhythmias spontaneously develop. Blue electrocardiograms associated with the blue region displayed normal sinus rhythm; red electrocardiograms associated with the red regions spontaneously developed an episode of torsades de pointes. To see this figure in color, go online.

Table 1

Effect of Drugs on Ion Channels

Drug	ICaLIC50 [nM]	IKrIC50 [nM]	Cmax [nM]
Ajmaline	71,000	1040	900
Amiodarone	270	30	0.3
Bepridil	211	33	21.5
Chlorpromazine	–	1470	20.5
Cibenzoline	30,000	22,600	739
Cisapride	–	6.5	3.8
Diltiazem	450	17,300	87.5
Dofetilide	60,000	5	1.2
Fluvoxamine	4900	3100	196
Haloperidol	1700	27	2.4
Mexiletine	100,000	50,000	2787
Nifedipine	60	275,000	5.4
Nitrendipine	0.3	10,000	1.6
Phenytoin	103,000	100,000	4250
Pimozide	162	20	0.6
Prenylamine	1240	65	13
Propranolol	18,000	2828	19
Quinidine	15,600	300	2080.5
Sertindole	8900	14	0.8
Tedisamil	–	2500	80
Terfenadine	375	8.9	4.5
Thioridazine	1300	33	593.5
Verapamil	100	143	53

IC50-values and effective free therapeutic concentration Cmax for the 23 drugs used in this study (20).

Measuring calcium transients in isolated cardiomyocytes

To characterize calcium transients, we isolated ventricular cardiomyocytes from the hearts of male Sprague Dawley rats with a weight of 250–300 g (Charles River Laboratories, Wilmington, MA). We anesthetized the rats with inhaled isoflurane and quickly removed the hearts from the chest after euthanasia. We retrograde perfused the hearts with Ca2+-free Tyrode buffer (140 mM NaCl, 5.4 mM KCl, 0.33 mM NaH2PO4, 0.5 mM MgCl2, 11 mM glucose, and 5 mM HEPES (pH7.4)) at 1.0 mL/min for 3 min, followed by an enzyme solution containing collagenase (1.0 mg/mL collagenase type II; Worthington Biochemical, Lakewood, NJ), protease (0.05 mg/mL, type XIV; Sigma-Aldrich, St. Louis, MO), and 0.1 mM Ca2+ for 7 min. To harvest the cardiomyocytes, we cut the ventricular tissue into small pieces and filtered it with a 250 μm nylon mesh. We gradually increased the calcium concentration of the Tyrode solution to 1.0 mM for the physiologic analysis and incubated the cardiomyocytes for 15 min with 1 μM Fura-2-AM (Invitrogen, Carlsbad, CA) in Tyrode (1.0 mM, Ca2+). We mounted the cardiomyocytes into a recording chamber on the stage of an Olympus IX-71 inverted microscope (Olympus, Center Valley, PA), where we stimulated them electrically at a frequency of 0.5 Hz. Using a galvanometer-driven mirror (HyperSwitch; IonOptix, Westwood, MA), we excited Fura-2 at a wavelength of 340 or 380 nm and recorded the emission at 510 nm using a photomultiplier (IonOptix ). After 5 min of incubation with the drug dofetilide at concentrations of 4, 8, 16, 38, and 130 nM, we recorded cardiomyocyte calcium fluorescence at 250 Hz for 8 min for n = 6 cells each and analyzed the recordings in real time.

Recording electrocardiograms in perfused Langendorff hearts

To record electrocardiograms, we harvested the hearts of male Sprague Dawley rats with a weight of 250–300 g (Charles River Laboratories). We excised the hearts from anesthetized rats (2.5% isoflurane in 95% oxygen and 5% carbon dioxide), immediately cannulated the aorta, connected it to a constant pressure perfusion Langendorff system (Harvard Apparatus, Cambridge, MA) with Krebs solution (118 mM NaCl, 4.75 mM KCl, 25 mM NaHCO3, 1.2 mM KH2PO4, 1.2 mM MgSO4, 1.5 mM CaCl2, 11 mM glucose, and 2 mM pyruvate), warmed it to 37°C, and bubbled it with 95% oxygen and 5% carbon dioxide. We instrumented the spontaneously beating hearts with electrocardiogram electrodes located at the apex and base. After 10 min of equilibration, we switched the perfusion system to a reservoir to expose the hearts to selected concentrations of dofetilide and nifedipine for a period of 5 min. For n ≥ 6 hearts in each group, we recorded the electrocardiogram by Animal Bio Amp (ADInstruments, Colorado Springs, CO) and monitored it continuously throughout the experiment and the a washout period using a Power Lab system (ADInstruments).

Experimentally characterizing the effect of drugs

We characterized the occurrence of arrhythmias in both the isolated cardiomyocytes and the perfused hearts. For the isolated cardiomyocytes, we counted an arrhythmia episode as one if at least one early afterdepolarization occurred within the recording period of 8 min and as zero otherwise. We then quantified the relationship between the prevalence of arrhythmia and the concentration of dofetilide using a nonlinear regression curve with a two-parameter equation. For the perfused hearts, we calculated the percentage of premature ventricular contractions of all heart beats during the last minute of drug administration. We defined ventricular tachycardia as three or more consecutive premature ventricular contractions. We analyzed the data using Student’s t-test for normally distributed data with equal variance between groups and the Mann-Whitney U test for all other data. For all analyses, we used Prism 7.

Results

IKr and ICaL enhance and prevent early afterdepolarizations

It has been shown that early afterdepolarizations are a precursor of torsades de pointes at the cellular level (14). To identify which ion channels have the most significant impact on the appearance of early afterdepolarizations, we perform 500 simulations of single midwall cells and systematically blocked seven ion channels: the L-type calcium current ICaL, the inward rectifier potassium current IK1, the rapid and slow delayed rectifier potassium currents IKr and IKs, the fast and late sodium currents INaP and INaL, and the transient outward potassium current Ito. Fig. 2 illustrates these seven ion channels within the O’Hara-Rudy model for ventricular cardiomyocytes (28). After determining the presence or absence of early afterdepolarizations for all simulations, we fit a logistic regression and extracted the marginal effects, a measure that quantifies the effect of each channel blockade on the probability of early afterdepolarizations. Our results in Fig. 2 show that of the seven channels, the rapid delayed rectifier potassium current IKr and the L-type calcium current ICaL have the most pronounced effects on early afterdepolarizations. Yet, these two currents display opposite effects: the rapid delayed rectifier potassium current IKr significantly increased the risk of early afterdepolarizations, whereas the L-type calcium current ICaL decreases the risk.

Figure 2

Effect of different ion channels on the probability of early afterdepolarizations. Positive values imply that blocking this ion channel enhances early afterdepolarizations; negative values imply that blocking prevents early afterdepolarizations. Blocking the rapid delayed rectifier potassium current IKr and the L-type calcium current ICaL has the strongest effect on enhancing and preventing early afterdepolarizations. To see this figure in color, go online.

IKr blockade triggers early afterdepolarizations in simulation and experiment

To validate our findings of the computational model, we use isolated rat ventricular cardiomyocytes and expose them to the drug dofetilide, which selectively blocks the rapid delayed rectifier potassium current IKr. We record calcium fluorescence and compare it to the calcium transients predicted by the computational model of human ventricular endocardial cells. Fig. 3 shows the development of early afterdepolarizations in the presence of the drug dofetilide, both in isolated rat cardiomyocytes and in the single-cell model. In both cases, the relationship between the probability of early afterdepolarizations and the concentration of the drug is dose dependent: increasing the dose of dofetilide increases the probability of early afterdepolarizations.

Figure 3

Early afterdepolarizations in single-cell simulation and experiment. Isolated rat cardiomyocyte (top left) and the probability of developing early afterdepolarizations in response to the drug dofetilide at concentrations of 4, 8, 16, 38, and 130 nM (n = 6 cells each; top right) are shown. Calcium transients in response to the drug dofetilide at 0, 16, and 130 nM in the computational simulation (bottom left) and experiment (bottom right) are shown. To see this figure in color, go online.

Machine learning classifies the boundary beyond which arrhythmias develop

According to our simulated probability of early afterdepolarizations at the single-cell level in Fig. 2, we select the two ion channels that most strongly enhance and prevent early afterdepolarizations, the rapid delayed rectifier potassium current IKr and the L-type calcium current ICaL. We use our high-fidelity human heart model (37) to simulate the effect of combined IKr and ICaL block at different concentrations (23). Our human heart model typically takes 40 h to run using 160 CPUs for 5 s of simulation. To reduce the computational cost, we adopt a particle-learning Gaussian process classifier with adaptive sampling to efficiently explore the parameter space. This method automatically places new samples near the classification boundary that divides the arrhythmic and nonarrhythmic domains to increase its resolution.

Fig. 4 summarizes the results of our proarrhythmic risk classification. The blue electrocardiograms were sampled at points in the blue region and display normal sinus rhythm. The red electrocardiograms were sampled at points in the red region and spontaneously develop torsades de pointes. The white contour indicates the classification boundary. The vertical axis reveals the proarrhythmic risk for a selective block of the rapid delayed rectifier potassium current IKr: at a critical IKr block of 70%, the risk classification changes from low (shown in blue) to high (shown in red), and the heart will spontaneously develop torsades de pointes. Moving horizontally to the right modulates the proarrhythmic risk for a combined block with the L-type calcium current ICaL: when combining IKr and ICaL block, the critical IKr block decreases below 70%. Strikingly, beyond an ICaL block of 60%, the heart will not develop fibrillation, no matter how high the IKr block. In agreement with our observations on the cellular level in Fig. 2, Fig. 4 supports the notion that certain channels can have a positive effect and mitigate torsadogenic risk upon rapid delayed rectifier potassium current block.

IKr and ICaL enhance and reduce the risk of ventricular arrhythmias

To explore the interaction between the rapid delayed rectifier potassium current IKr and the L-type calcium current ICaL at the organ level, we combine computational modeling and isolated Langendorff perfused rat heart preparations using two different drugs: dofetilide, which selectively blocks the rapid delayed rectifier potassium current IKr, and nifedipine, which selectively blocks the L-type calcium current ICaL. We probe different concentrations of these two drugs and determine the presence of arrhythmias from the computational and experimental electrocardiograms. Fig. 5, top illustrates our Langendorff perfused heart, our four drug concentrations visualized in the proarrhythmic risk estimator, and the risk of premature ventricular contractions and arrhythmias for these four cases. Fig. 5, bottom shows the electrocardiograms in response to dofetilide at 0 and 20 nM combined with nifedipine at 0, 60, and 480 nM both for the computational simulation (left) and the experiment (right). For the baseline case without drugs, both the computational model and experimental system display normal sinus rhythm (first row). Blocking the rapid delayed rectifier potassium current IKr by administering dofetilide beyond a critical concentration induces arrhythmias both computationally and experimentally, second row, an observation that agrees well with the single-cell simulation and experiment in Fig. 3. Additionally, blocking the L-type calcium current ICaL by coadministering a small concentration of nifedipine markedly alters the excitation pattern both computationally and experimentally but still triggers irregular beats. Increasing the L-type calcium current ICaL block by coadministering a large concentration of nifedipine removes the risk of arrhythmias both computationally and experimentally; the hearts excite at a regular pattern, but at a slightly different rate than for the baseline case without drugs.

Figure 5

Ventricular arrhythmias in whole-heart simulation and Langendorff perfused hearts. Preparation of isolated rat heart (top left), four drug concentrations visualized in the proarrhythmic risk classification estimator (top middle), and risk of premature ventricular contractions and arrhythmias in response to varying concentrations of drugs dofetilide and nifedipine (n ≥ 6, ∗p < 0.05 compared to (1), #p < 0.05 compared to (2); top right) are shown. Dofetilide selectively blocks the rapid delayed rectifier potassium current IKr; nifedipine selectively blocks the L-type calcium current ICaL. Electrocardiograms in response to dofetilide at 0 and 20 nM combined with nifedipine at 0, 60, and 480 nM in the computational simulation (bottom left) and experiment (bottom right) are shown. To see this figure in color, go online.

Critical drug concentrations are a predictor of drug toxicity

To validate our approach, we calculate the critical concentrations for 23 common drugs using the risk assessment tool in Fig. 4. In essence, the individual block-concentration characteristics for each drug (3,31) map onto a trajectory in the IKr/ICaL plane of the risk assessment diagram. The intersection of this trajectory with the classification boundary defines the critical drug concentration. Curves that never cross the classification boundary indicate a safe drug.

Fig. 6 demonstrates that our classification boundary in Fig. 4 can reliably stratify the risk of 23 common drugs using their computationally predicted critical concentration: 14 drugs are classified as high-risk drugs. Of those, thioridazine and quinidine cross the classification boundary at the lowest concentrations of 0.1× and 0.3×, and chlorpromazine and amiodarone at the highest concentrations of 154.9× and 282.6×. Nine drugs are classified as low-risk drugs. Of those, propranolol crosses the classification boundary at 474.6×, and all other drugs never cross the classification boundary. To validate our simulation-based classifier, we used the predicted critical concentrations from Fig. 6 to train data-based Gaussian process classifiers using leave-one-out cross validation (20). Specifically, we trained 23 classifiers by excluding one drug from the training data and predicting the risk category of the excluded drug. With this method, we were able to correctly classify 22 of the 23 compounds into the risk categories predicted by our simulation-based classifier from Fig. 4 using the critical drug concentrations. For validation with the literature, each drug trajectory graph in Fig. 6 also displays the reported torsadogenic risk categories from 1 to 5 (20) and, in red and blue colors, the reported classification into torsadogenic and nontorsadogenic compounds (19).

Figure 6

Risk stratification of 23 drugs using our proarrythmic risk classification. Black and white regions indicate fibrillating and nonfibrillating regimes; red and blue curves represent the IKr/ICaL profiles of high- and low-risk drugs at varying concentrations; gray dots and numbers indicate the critical concentration at which the drug curves cross the classification boundary as predicted by our proarrhythmic risk classification in Fig. 4. For comparison, numbers from 1 to 5 indicate the reported torsadogenic risk (20); red and blue colors of the numbers indicate torsadogenic and nontorsadogenic compounds (19). To see this figure in color, go online.

Fig. 7 illustrates an experimental validation of our risk stratification for 12 drugs in isolated rabbit hearts. The dots represent our predicted critical concentrations according to Fig. 6; the black lines represent the experimentally measured critical concentrations at which those drugs induced early afterdepolarizations, torsades de pointes, or ventricular tachycardia. All experimental values were taken from (42), except for haloperidol (43) and quinidine (44). The seven drugs that our method classified as proarrhythmic—quinidine, cisapride, terfenadine, bepridil, dofetilide, haloperidol, and sertindole—also induced arrhythmias in the rabbit heart. Our computationally predicted critical concentrations for three drugs—quinidine, dofetilide, and haloperidol—coincide in the experimentally observed critical concentration range; the other four lie within one order of magnitude of this range. Four out of five drugs that our method classified as safe—nitrendipine, nifedipine, diltiazem, and verapamil—were also safe in the rabbit heart. The one drug for which we see a discrepancy, cibenzoline, is the only drug that induced ventricular tachycardia before early afterdepolarizations and torsades de pointes. This arrhythmic event could have been initiated by another mechanism because cibenzoline has been classified as safe for torsades de pointes (4).

Figure 7

Experimental validation of risk stratification for 12 drugs. Black lines represent the experimentally measured critical concentration in isolated rabbit hearts (42, 43, 44). Stars indicate classification by early afterdepolarization, torsades de pointes, and ventricular tachycardia. Dots represent our predicted critical concentration, with red indicating proarrhythmic and blue safe drugs (20). To see this figure in color, go online.

Fig. 8 illustrates a computational validation of our risk stratification for three drugs: terfenadine, bepidril, and verapamil. Our stratification classifies terfenadine and bepidril as high risk and verapamil as safe. To validate this classification, we apply all three drugs at 10× their effective free therapeutic concentration. Terfenadine, with a critical concentration of 4.4×, triggers an arrhythmia immediately after the first beat; bepidril, with a critical concentration of 4.9×, triggers an arrhythmia after the second beat; and verapamil, which never crosses the classification boundary, is nonarrhythmogenic. Although all three drugs initiate a similar degree of blockade of the rapid delayed rectifier potassium current IKr of 84, 86, and 79%, their blockade of the L-type calcium current ICaL of 11, 50, and 84% varies significantly. These three examples, now with a complete simulation, highlight the interaction of different channels and confirm the predictive power of our proarrhythmic risk estimator in Fig. 4 and its resulting risk stratification in Fig. 6.

Figure 8

Computational validation of risk stratification for three drugs applied at the same concentration. At 10× the effective free therapeutic concentration, terfenadine blocks 84% of IKr and 11% of ICaL, bepidril blocks 86% of IKr and 50% of ICaL, and verapamil blocks 79% of IKr and 84% of ICaL. The different degrees of blockade result in arrhythmic patterns for terfenadine and bepidril, but not for verapamil, for which the high degree of ICaL block prevents the development of arrhythmia and slows the beating rate. To see this figure in color, go online.

Discussion

Current drug screening paradigms are expensive, time consuming, and conservative. Here, we propose a new approach that integrates knowledge from the ion channel, single-cell, and whole-heart levels via computational modeling and machine learning to reliably predict the cardiac toxicity of new and existing drugs. Our results are based on a sensitivity analysis that identifies a pair of counteracting ion channels, IKr and ICaL, that play the most significant role in enhancing and reducing arrhythmogenic risk. We combine multiscale experiments, multiscale simulation, high-performance computing, and machine learning to create a risk classifier that allows us to identify the proarrhythmic potential of existing and new drugs, either in isolation or combined with other drugs. Naturally, our multiscale model is more complex than the single-cell models that are currently used to study drug-induced arrhythmias. However, the underlying idea of our approach is not to advocate for a battery of complex simulations for a single-purpose use but to create a general-purpose proarrhythmic risk classification diagram to directly identify the critical concentration of any drug and classify the drug into proarrhythmic or safe simply by means of its trajectory within this diagram. Collectively, our study provides new insights that are significant in the development of new compounds. Our efforts markedly extend current initiatives by pharmaceutical industries, clinical researchers, and regulatory agencies with the common goal to develop a new testing paradigm for a more accurate and comprehensive mechanistic assessment of new drugs.

Early afterdepolarizations are a multiple-channel phenomenon

At the single-cell level, we have shown that early afterdepolarizations are triggered when the rapid delayed rectifier potassium current IKr is blocked above a certain level. This is in line with the current regulatory framework, which identifies this channel as the most relevant for QT interval prolongation and torsades de pointes initiation (4). However, through computational modeling, we have demonstrated that early afterdepolarizations are better conceptualized as a multichannel phenomenon. Our sensitivity analysis in Fig. 2 identifies the rapid delayed rectifier potassium current IKr and the L-type calcium current ICaL as the most relevant currents for the formation of early afterdepolarizations. These two channels have opposing effects: blocking IKr can initiate and blocking ICaL can prevent early afterdepolarizations. In a recent study, we have found a similar trend at the QT interval level (25), which is also considered in current regulations (5). These results are in line with other studies that have highlighted the importance of altered calcium dynamics during early afterdepolarizations (14,45, 46, 47) and, more recently, also during delayed afterdepolarizations (48). These multichannel effects between the rapid delayed rectifier potassium current IKr and the L-type calcium current ICaL observed in Fig. 5 open the door toward a systematic search for blockade combinations that can offset the torsadogenic effects of IKr block alone (49).

IKr and ICaL modulate the onset of torsades de pointes

Our study shows that the rapid delayed rectifier potassium current IKr and the L-type calcium current ICaL not only determine the onset of early afterdepolarizations but also the development of torsades de pointes. Our results in Fig. 5 suggest that blocking the L-type calcium current ICaL can prevent the development of arrhythmias, even at high levels of rapid delayed rectifier potassium current IKr blockade, both in our high-resolution model and in isolated rat hearts. Recent studies have pointed out this preventive role of ICaL. An analysis of 55 compounds showed that adding the effects of ICaL blockade to IKr block improved the predictive potential, whereas adding the effects of INaL did not (50). However, this study only demonstrated correlation without a mechanistic explanation. A recent machine-learning-based approach suggested that risk prediction of torsades de pointes could be improved by including intracellular calcium currents (19). This trend was confirmed by a recent study that classified drugs in terms of IKr and ICaL blockade metrics (18). At the cellular level, these findings reflect the importance of these currents in the development of early afterdepolarizations (14). At the whole-heart level, the presence of these action potential abnormalities is a necessary but not sufficient condition to initiate torsades de pointes; here, heterogeneities (51,52) and electrotonic effects (12,21) play a major role in the propagation of this type of arrhythmia.

The degree of toxicity correlates with the critical drug concentration

We have classified drugs based on their critical concentration, the concentration at which they cross the classification boundary of our risk estimator in Fig. 4. Critical-concentration-based methods have been used both in rabbit models (42) and in computational models (53). Here, we successfully employed this concept by inducing arrhythmias at elevated drug concentrations both computationally and experimentally. Critical concentrations can be interpreted as the distance from an event of torsades de pointes: the higher the normalized concentration, the further away is the baseline concentration, and thus, the safer the compound. When using the critical drug concentration to stratify the risk of drugs in Fig. 6, we correctly identify quinidine, bepridil, dofetilide, chlorpromazine, cisapride, and terfenadine as high-risk and diltiazem, mexiletine, and verapamil as low-risk drugs, similar to a classifier based on net current (54). Fig. 8 confirms the high-risk action of terfenadine and bepridil and the low-risk action of verapamil, which is widely known as a calcium channel blocker with antifibrillatory effects (55). Moreover, we correctly identified 22 compounds as high and low risk in Fig. 6, compared to the reported high-risk categories 1–3 and low-risk categories 4–5 (20). For these 22 compounds, our classifier also agrees exactly with a recent machine-learning classifier based on action potential duration and diastolic calcium (19). To eliminate sources of noise in the evaluation of our model, we have only considered those drugs for which 70% or more of the published studies agreed on their risk classification (41). The only drug that our approach classifies incorrectly is propanolol, which has a critical concentration of 474.6× of the effective free therapeutic concentration. Although Fig. 6 suggests that this concentration is significantly higher than for all other high-risk drugs, the classifier is trained without any other compound similar to propanolol when performing leave-one-out cross validation. If more data were available, the predictive power of our classifier could be improved. Nonetheless, the potential of our approach lies in supporting the successful progression of compounds that have a poor selectivity to the rapid delayed rectifier potassium current alone and would, under current paradigms, be falsely discontinued through the drug discovery and development process. Our study suggest that our approach correctly identifies those drugs. Our risk estimator in Fig. 4 allows us to quickly and reliably screen the proarrhythmic potential of any drug, either in isolation or in combination with other drugs. Additionally, our method provides a direct and mechanistic biomarker to predict drug-induced arrhythmias to complement ad hoc indicators derived from single-cell models (20,56).

Limitations

Although our proposed method holds promise to rapidly assess the risk of a new drug, it has a few limitations: first, our major focus was on combining computational modeling and machine learning to create efficient risk estimators (8,57); long term, more experiments will be needed to better validate the method and broaden its scope and use. Second, our model is only as good as its input, the concentration-block curves; we have addressed this limitation in a separate study (25), similar to other groups (31,58), and found that there is a mismatch between the drugs that have been well characterized experimentally (3)—the input of the classifier—and the drugs that we agree on in their risk classification—the output of the classifier; to mitigate this limitation, we used a deterministic approach to classify the set of compounds. Third, we have used a simple pore-block model to include the effect of drugs on ion channel currents. In the future, we will investigate the effects of more sophisticated kinetics models (56). However, incorporating these models will limit the applicability of our method to a few compounds for which kinetic data are available. Fourth, our work has mainly followed recommendations by the CiPA initiative (6); it will be important to validate our model against other cell and heart models and, probably most importantly, against other compounds. Here, we selected a reduced set of compounds to test our classifier. There is a broad agreement on the risk classification of these drugs, which makes them an ideal data set to validate our model. Nonetheless, a larger data set of both measured drug effects and risks would be desirable to fully validate our approach (59). Fifth, we have based our initial studies on reported experiments and clinical observations, supplemented with our own cell-level and isolated heart studies with rodent hearts and reported rodent heart studies from the literature. It is an ongoing discussion to what extent studies of rat hearts can provide insight into human QT interval prolongation (60). There is reported evidence of the presence of rapid delayed rectifier potassium current IKr (61, 62, 63) and electrophysiological changes in response to dofetilide (64) in rats. Moreover, here we focus exclusively on the regime of early afterdepolarizations and arrhythmias, which we directly observe in our experiments. A critical and logical next step would be to validate our method using our own independent experiments with human adult cardiomyocytes, in larger animals, and, ideally, in healthy human volunteers. Ultimately, with a view toward precision cardiology, our approach has the potential to combine the personalized block-concentration characteristics and personalized cardiac geometries toward identifying the optimal course of care for each individual patient (65,66).

Conclusions

We propose a novel strategy toward drug screening. This was only possible by combining cutting edge technologies of multiscale exposure-response simulation, machine learning, and high-performance computing. Using systematic sensitivity analyses, we identified the L-type calcium channel as a critical antagonist to the rapid delayed rectifier potassium current in modulating arrhythmogenic risk. Our simulations highlight the mechanisms by which drug-induced arrhythmias propagate across scales, from modifications at the ion channel level via early afterdepolarizations at the cellular level to rapid oscillations in the electrocardiogram at the whole-heart level. Using machine learning, we integrate information from different scales and sources, experimental and computational, into a single risk estimator. Our results suggest that this proarrhythmic risk estimator can rapidly and reliably stratify any drug based on block-concentration characteristics from single-cell recordings. Our study provides a more holistic insight into the generation of drug-induced arrhythmias than current single-cell studies alone. We envision that our findings will help accelerate drug development and reduce the cost to deliver safe and effective drugs to patients.

Author Contributions

F.S.-C., K.S., E.A., and E.K. designed the research. F.S.-C. and performed the simulations. K.S. performed the experiments. F.S.-C., K.S., E.A., and E.K. analyzed the data. F.S.-C., K.S., E.A., and E.K. wrote the manuscript.