Multidimensional estimation of cardiac arrhythmia potential across space and time

All medications have adverse effects. Among the most serious of these are cardiac arrhythmias1. The current safety evaluation for cardiac toxicity involves interrogating effects of the drug on the delayed rectifier potassium current in single cells2 and the QT interval in healthy volunteers3. However, this paradigm for drug safety evaluation is costly4, lengthy5, conservative3, and impede efficient drug development6. Here we combine multiscale experiment and simulation, high-performance computing, and machine learning to create a multidimensional risk estimator to quickly and reliably stratify new and existing drugs according to their pro-arrhythmic potential. We capitalize on recent developments in machine learning and integrate information across ten 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 Gaus-

All pharmacological agents have the potential to impact cardiac repolarization and, with it, cause torsades de pointes, a ventricular arrhythmia characterized by rapid, irregular patterns in the electrocardiogram 10 . 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 11 . Predicting this potentially fatal heart rhythm is challenging given the complex interplay between genetic predisposition and medications.
While 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 12 . There is a general agreement between clinical researchers, pharmaceutical companies, and regulatory agencies that computational tools should play a more central role in the pro-arrhythmic risk assessment of new drugs 13-15 . However, 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 16 . We have recently proposed a novel exposure-response simulator that allows us visualize how different drugs modulate the cardiac electrophysiology across ten orders of magnitude in space and time 17 . Here, we combine this simulator with machine learning Figure 1: Hybrid computational-experimental approach to characterize the proarrhythmic potential of existing and new drugs. We characterize calcium transients in cardiomyocytes in response to drugs, both computationally (top) and experimentally (bottom) and identify the ion channels that most likely generate early afterdepolarizations (left). We screen the blockade space of the two most relevant channels and identify the pro-arrhythmic classification boundary using high performance computing and machine learning (center). We validate our approach using electrocardiograms in simulations and isolated perfused hearts (right). We demonstrate the potential of our classifier by risk stratifying 23 drugs and comparing the result against their reported risk categories.
4 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/545863 doi: bioRxiv preprint first posted online Feb. 10, 2019; Increasing evidence suggests that early afterdepolarizations are a precursor of torsades de pointes at the cellular level 20, 21 . To identify which ion channels have the most significant impact on the appearance of early afterdepolarizations, we model single midwall cells and block seven ion channels and monitoring for early afterdepolarizations. Figure 2, top left, illustrates these seven ion channels within the O'Hara Rudy model for ventricular cardiomyocytes 22 . We fit a logistic regression and extract the marginal effects to quantify the effect of each channel blockade on the probability of early afterdepolarizations. Our results show that of the seven channels, the rapid delayed rectifier potassium current I Kr and the L-type calcium current I CaL have the most pronounced effects on early afterdepolarizations. Yet, these two currents display opposite effects: The current I Kr significantly increases the risk of early afterdepolarizations, while I CaL decreases the risk. To validate these findings, we isolate rat ventricular cardiomyocytes and expose them to the drug dofetilide, which selectively blocks I Kr . We record calcium fluorescence and compare it to the calcium transients predicted by the computational model of human ventricular endocardial cells. Figure 2, middle row, 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. We select the two ion channels that most strongly enhance and prevent early afterdepolarizations, I Kr and I CaL to study torsades de pointes. We use our high resolution human heart model 23 to simulate the effect of combined I Kr and I CaL block at different concentrations 17 . We 5 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.  6 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/545863 doi: bioRxiv preprint first posted online Feb. 10, 2019; adopt a particle learning Gaussian process classifier with adaptive sampling to efficiently explore the parameter space, probing the points of maximum entropy 24 . Figure 2, bottom right, summarizes the results of our pro-arrhythmic risk classification. The vertical axis reveals the pro-arrhythmic risk for a selective block of the I Kr : At a critical I Kr 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 pro-arrhythmic risk for a combined block with I CaL : When combining I Kr and I CaL block, the critical I Kr block increases above 70%. Beyond an I CaL block of 60%, the heart will not develop arrhythmia, no matter how high the I Kr block.
To validate these predictions we took advantage of an isolated Langendorff perfused rat heart preparation exposed to two different drugs, dofetilide, which selectively blocks the I Kr and nifedipine, which selectively blocks I CaL . Figure 3 shows that for the baseline case without drugs, both the computational model and experimental system display normal sinus rhythm, first row. Blocking I Kr 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 Figure 2, middle row. Additionally blocking I CaL by co-administering a small concentration of nifedipine markedly alters the excitation pattern both computationally and experimentally, but still triggers irregular beats. Increasing the I CaL block by co-administering a large concentration of nifedipine removes the risk of arrhythmias both computationally and experimentally, the hearts excite at a regular pattern. 7 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/545863 doi: bioRxiv preprint first posted online Feb. 10, 2019; 8 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.  9 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/545863 doi: bioRxiv preprint first posted online Feb. 10, 2019; To validate our approach, we calculate the critical concentrations for 23 drugs using the risk assessment classifier in Figure 2

11
. CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

13
. CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. 14 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.

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 C m is the membrane capacitance and I ion (φ, q) is the ionic current, which we represented as a function of the transmembrane potential φ and a set of state variables q 1 .
The state variables obey ordinary differential equations,q = g(φ, q), as functions of the transmembrane potential φ and their current values q 2 . For our single cell simulations, we used ventricular cardiomyocytes with 15 ionic currents and 39 state variables 3 , I ion = I Kr + I Ks + I K1 + I CaL + I Na + I CaNa + I CaK + I Cab + I Nab + I Kb + I to + I NaK + I pCa + I NaCa,i + I NaCa,ss , with a minor modification 4 of the fast sodium current I NaP 5 . We parameterized the model for human midwall cells 3 , and modeled the effect of drugs by selectively blocking the relevant ionic currents I ion 6 . 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 16 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/545863 doi: bioRxiv preprint first posted online Feb. 10, 2019; electrophysiology 7, 8 , and scale the ionic current I i by this fractional block 9 , We studied the relative importance of seven ion channels, I CaL , I K1 , I Kr , I Ks , I NaL , I NaP , and I to on inducing early afterdepolarizations. To achieve a steady state, we paced the cells for 600 cycles at a frequency of 1Hz. We defined the presence of early afterdepolarizations as the occurrence of a change in potential greater than 0.1mV/ms between the 50 and 1000ms of the last two recorded cycles 10 . 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 fit 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. We normalized the results by the maximum 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 9 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 11 . The basis of this model is the classical monodomain model that characterizes the spatio-temporal evolution of the transmembrane potential φ through the following partial differential equation, In addition to the local source term I ion /C m from equation (1), the transmembrane potential depends on the global flux term div(D · ∇φ), where D is the conductivity tensor that accounts for a fast signal propagation of D parallel to the fiber direction f and a slow signal propagation of D ⊥ perpendicular to it 1 , We used the O'Hara Rudy model 3 from equation (2)  * 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 combine our computational models with machine learning techniques 15,16 . Briefly, to characterize ventricular fibrillation, we performed n = 40 human heart simulations and employed a particle learning method to systematically sample the classification boundary 18 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/545863 doi: bioRxiv preprint first posted online Feb. 10, 2019; within the parameter space. To identify the boundary that divides the arrhythmic and nonarrhythmic domains, we used a Gaussian process classifier and adaptively sampled the point of maximum entropy 17 . We generated the first n = 10 samples from a latin hypercube design, and sampled the remaining n = 30 samples adaptively. 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 pro-arrhythmic risk estimator in Figure 2, bottom right, and validated our approach against the known risk classification of these drugs. To select the compounds, we began with a list 31 drugs 7 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 18,19 , and did not consider the remaining eight controversial drugs. Table 1  (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/545863 doi: bioRxiv preprint first posted online Feb. 10, 2019; taneously beating hearts with ECG electrodes located at the apex and base. After ten minutes 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 five minutes.
For n≥6 hearts in each group, we recorded the ECG by Animal Bio Amp (AD Instruments, Colorado) and monitored it continuously throughout the experiment and the a washout period using a Power Lab system (AD Instruments, Colorado).
* 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 eight minutes, and as zero otherwise. We then quantified the relationship between the prevalence of arrhythmia and the concentration of dofetilide using a non-linear 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.

* Data Availability
The data that support the findings of this study are available from the corresponding 21 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/545863 doi: bioRxiv preprint first posted online Feb. 10, 2019; author upon reasonable request. 22 . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint . http://dx.doi.org/10.1101/545863 doi: bioRxiv preprint first posted online Feb. 10, 2019; . CC-BY 4.0 International license It is made available under a (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.