Three-Dimensional Heart Model-Based Screening of Proarrhythmic Potential by in silico Simulation of Action Potential and Electrocardiograms

The proarrhythmic risk is a major concern in drug development. The Comprehensive in vitro Proarrhythmia Assay (CiPA) initiative has proposed the JTpeak interval on electrocardiograms (ECGs) and qNet, an in silico metric, as new biomarkers that may overcome the limitations of the hERG assay and QT interval. In this study, we simulated body-surface ECGs from patch-clamp data using realistic models of the ventricles and torso to explore their suitability as new in silico biomarkers for cardiac safety. We tested seven drugs in this study: dofetilide (high proarrhythmic risk), ranolazine, verapamil (QT increasing, but safe), bepridil, cisapride, mexiletine, and diltiazem. Human ventricular geometry was reconstructed from computed tomography (CT) images, and a Purkinje fiber network was mapped onto the endocardial surface. The electrical wave propagation in the ventricles was obtained by solving a reaction-diffusion equation using finite-element methods. The body-surface ECG data were calculated using a torso model that included the ventricles. The effects of the drugs were incorporated in the model by partly blocking the appropriate ion channels. The effects of the drugs on single-cell action potential (AP) were examined first, and three-dimensional (3D) body-surface ECG simulations were performed at free Cmax values of 1×, 5×, and 10×. In the single-cell and ECG simulations at 5× Cmax, dofetilide, but not verapamil or ranolazine, caused arrhythmia. However, the non-increasing JTpeak caused by verapamil and ranolazine that has been observed in humans was not reproduced in our simulation. Our results demonstrate the potential of 3D body-surface ECG simulation as a biomarker for evaluation of the proarrhythmic risk of candidate drugs.

The proarrhythmic risk is a major concern in drug development. The Comprehensive in vitro Proarrhythmia Assay (CiPA) initiative has proposed the JTpeak interval on electrocardiograms (ECGs) and qNet, an in silico metric, as new biomarkers that may overcome the limitations of the hERG assay and QT interval. In this study, we simulated body-surface ECGs from patch-clamp data using realistic models of the ventricles and torso to explore their suitability as new in silico biomarkers for cardiac safety. We tested seven drugs in this study: dofetilide (high proarrhythmic risk), ranolazine, verapamil (QT increasing, but safe), bepridil, cisapride, mexiletine, and diltiazem. Human ventricular geometry was reconstructed from computed tomography (CT) images, and a Purkinje fiber network was mapped onto the endocardial surface. The electrical wave propagation in the ventricles was obtained by solving a reaction-diffusion equation using finite-element methods. The body-surface ECG data were calculated using a torso model that included the ventricles. The effects of the drugs were incorporated in the model by partly blocking the appropriate ion channels. The effects of the drugs on single-cell action potential (AP) were examined first, and three-dimensional (3D) body-surface ECG simulations were performed at free Cmax values of 1×, 5×, and 10×. In the single-cell and ECG simulations at 5× Cmax, dofetilide, but not verapamil or ranolazine, caused arrhythmia. However, the non-increasing JTpeak caused by verapamil and ranolazine that has been observed in humans was not reproduced in our simulation. Our results demonstrate the potential of 3D body-surface ECG simulation as a biomarker for evaluation of the proarrhythmic risk of candidate drugs.

INTRODUCTION
The proarrhythmic effects of cardiac and non-cardiac drugs have comprised a major drug safety issue for the past 20 years (Zipes, 1987;De Ponti et al., 2002;Moro et al., 2010). The electrocardiogram (ECG) is an effective means of determining whether a drug is proarrhythmic. Under certain conditions, prolongation of the QT interval increases the risk of developing Torsades de pointes (TdP), which can lead to sudden cardiac death (Thomas and Behr, 2016). The Comprehensive in vitro Proarrhythmia Assay (CiPA) was recently proposed to improve the accuracy of drug safety prediction during preclinical and clinical development (Vicente et al., 2018;Wallis et al., 2018). The CiPA comprises in silico simulation of several ion-channel assays and ECG studies to identify biomarkers of false-positive results of single hERG channel assays and thorough QT (TQT) studies (Darpo, 2010) performed according to the International Council for Harmonisation (ICH) S7B and E14 guidelines. There have been a large number of studies that investigated the effect of drugs on ECG using in silicon three-dimensional (3D) heart model. Zemzemi et al. (2013) examined the effect of the block of ion channels on ECG parameters. Okada et al. (2018) generated an arrhythmic hazard map under multiple ion channel blocks. Sahli Costabal et al. (2019) investigated the critical drug concentration which induced torsade de pointes. Rivolta et al. (2017) performed sensitivity analysis of JT peak and T-wave morphology parameters.
In this study, we further examined the utility of 3D ECG simulation in evaluating drug safety by simulating ECG at relatively high concentrations of drugs using realistic models of the ventricles and torso. We tested dofetilide, bepridil, cisapride, ranolazine, verapamil, mexiletine, and diltiazem using the 3D model and examined the morphologies of the simulated ECG data according to drug concentration. Dofetilide, bepridil, and cisapride are high or intermediateproarrhythmic-risk drug that prolongs the QT by blocking hERG. Verapamil and ranolazine are "false positive" lowproarrhythmic-risk drugs; they induce prolonged QT by hERG blockade while simultaneously blocking inward Ca 2+ (verapamil) and Na + (ranolazine) ion channels (Vicente et al., 2018). Mexiletine and diltiazem are lowproarrhythmic-risk drugs that do not prolong the QT at all. Recently, CiPA researchers proposed a new ECG biomarker, JTpeak, which may enable the identification of drugs producing false-positive results (Vicente et al., 2018). Proarrhythmic drugs prolong the QT and the JTpeak due to hERG blockade, but not when the hERG blockade is offset by simultaneous blockade of other depolarizing ion channels (as by verapamil and ranolazine: only the QT prolonged but not the JTpeak). In this study, we explored the ability of the results of 3D ECG simulations Abbreviations: AP, action potential; CT, computed tomography; EAD, early afterdepolarization; ECG, electrocardiogram; ORd, O'Hara-Rudy dynamic; TdP, Torsades de pointes; TDR, transmural dispersion of repolarization. to identify false-positive results independently of clinically obtained ECG data.

ECG Simulation Using Models of the Ventricles and Torso
The model construction and ECG simulation are also described in our previous papers (Im et al., 2008;Lim et al., 2013;Ryu et al., 2019). Human ventricular geometry and torso were from our previous studies (Lim et al., 2013;Ryu et al., 2019) (Figures 1A-C). Human ventricular geometry was reconstructed from the computed tomography (CT) images obtained from the University of Ulsan Medical Center using a commercially available software Aquarius intuition (TeraRecon Inc., San Mateo, CA, United States). Tetrahedral mesh was generated inside the 3D ventricular model using an inhouse software ( Figure 1A). The number of grid element was 1,475,818. For the modeling of Purkinje fibers, the 2dimensional representation of the Purkinje network shown in the paper by Berenfeld and Jalife (1998) was digitized, scaled to the size of the 3D model, and mapped onto the endocardial surface of the 3D model of the ventricles ( Figure 1B). Pacing was applied at the location of His bundle. The model of Purkinje fibers simply transmits the electrical signal unidirectionally. The speeds of signal transmission at various regions were adjusted manually so that the simulated activation map matches that of clinical data (Durrer et al., 1970; Supplementary Figure S2). The end nodes of the Purkinje network stimulated myocardium by applying stimulation current of −80.0 A/F until the membrane potential exceeds −10 mV. The endocardial nodes connected to the node nearest to each end node of the Purkinje network were considered the Purkinje-muscle junction (PMJ). All the tetrahedral elements containing the PMJ nodes were stimulated. Signal propagation from the stimulation nodes throughout the tissue was obtained by solving a reaction-diffusion equation (Eq. 1) numerically: where V m is the transmembrane potential, t is time, D is the diffusion tensor, C m is the membrane capacitance, and I ion and I stim are ionic and stimulation currents, respectively. The equation was spatially discretized by using finite element and the time derivative was approximated by forward Euler method (Im et al., 2008).
To calculate I ion , the O'Hara-Rudy dynamic (ORd) human ventricular cell model was used (O'Hara et al., 2011). ORd model has three types of cells: endocardial, M, and epicardial cells. Each cell type was assigned to the ventricular wall with reference to the figure shown in the paper by Trudel et al. (2004). We also tested three recently published optimized cell models (Mann et al., 2016;Dutta et al., 2017;Krogh-Madsen et al., 2017). To calculate ECG values, the boundary element model of the human torso proposed by Potse et al. (2009) was used ( Figure 1C). The ECG was calculated by computing the potentials on the torso surface using the following equation (Potse et al., 2009): where φ ek (r) is potential at a point r on surface k. σ − k and σ + k are the conductivity inside and outside the surface k, respectively, J c is the source current density field, and r and r are variables.
The summation is over all surfaces l. d rr is the solid angle subtended at r by the infinitesimal surface element located at r . The key parameters of simulation are shown in Table 1.

Incorporation of the Effects of Drugs in the ECG Simulation
For the simulation of drug effect on ECG, we used the parameter values obtained by CiPA researchers to compare with their clinical data (Crumb et al., 2016;Li et al., 2017Li et al., , 2019. We tested seven drugs: dofetilide, bepridil, cisapride, verapamil, ranolazine, mexiletine, and diltiazem. The effects of drugs were incorporated in the ECG simulation by partly blocking the corresponding ion channels (I Na , I NaL , I CaL , and I Kr ) in the ionic-current model. The percentage of blockage of each ionic current was calculated using the Hill equation (Goutelle et al., 2008). The Cmax, IC 50 , and Hill coefficient values for each drug with respect to each ionic current were adopted from the literature ( Table 2) (Crumb et al., 2016;Li et al., 2017Li et al., , 2019. Cmax means free Cmax unless otherwise stated. The effects of each drug on single-cell action potentials  The percentages were calculated using the Hill equation applying the Cmax, IC50, and Hill coefficient values adopted from the literatgure (Crumb et al., 2016;Li et al., 2017Li et al., , 2019 [of endocardial (endo), epicardial (epi), and mid-myocardial (M) cells] were examined first, and 3D ECG simulations were performed at 1×, 5× and 10× Cmax.  Table S1). The Cmax values are listed in Table 2. The increases in the 90% AP duration (APD 90 ) for the three drugs in endocardial, M, and epicardial cells compared with the no-drug control are shown in Table 3. Among the three drugs, dofetilide induced the greatest increase in the APD 90 value, followed by ranolazine and verapamil. When Cmax was increased from 1× to 10×, APD 90 increased for all three drugs and all three cell types, with the exception of M cells, in the presence of dofetilide. Dofetilide induced  (Figure 2). Table 4 shows the transmural dispersion of repolarization (TDR) values, calculated as the difference between the largest and smallest APD 90 s among the endocardial, M, and epicardial cells. The TDR was largest in the case of dofetilide, and verapamil did not alter the TDR at a Cmax of 1× compared with the drug-free control ( Table 4). At a Cmax of 10×, verapamil increased the APD 90 of epicardial cells to a greater degree than that of M cells (Table 3), which resulted in a decreased TDR compared with the drug-free control ( Table 4).

Effects on 3D ECG Parameters
To examine the effects of drug concentration on ECG parameters, 3D ECG simulations were performed using the conditions shown Frontiers in Physiology | www.frontiersin.org  FIGURE 4 | Effects on body-surface ECG parameters. Body-surface ECG data (lead I) are shown for 7 drugs at Cmax values of 1×, 5×, and 10× compared with drug-free conditions. The QT interval was longest for dofetilide.
in Figure 2. Figure 4 shows the simulated ECGs (lead I) for the seven drugs according to concentration. Dofetilide resulted in the greatest increase in the QTc value at 1× Cmax (Table 4). At a Cmax value of 5×, dofetilide induced ventricular flutter; at 10× Cmax, dofetilide, cisapride, and ranolazine induced ventricular flutter. In contrast to findings reported by the CiPA researchers (Vicente et al., 2018), the JT peakC value increased with the QTc value for ranolazine and verapamil (  shows ECGs obtained from using different optimized cell models for the seven drugs at 1× Cmax. Dofetilide exhibited relatively long QT interval in all the cell models except for the model of Krogh-Madsen et al. (2017) in which the ECG morphology was irregular. The amplitude of the T wave was largest in the case of Dutta et al. (2017) while the model of Mann et al. (2016) exhibited the smallest T wave amplitude. Table 5 shows JT peak c prolongation of drugs for different optimized cell models. The optimized cell models resulted in JT peak c prolongations which are more consistent with clinical observations than the original ORd model. For ranolazine, metabolites seem to play a significant role in drug binding (Moreno et al., 2013). Because dofetilide induced ventricular tachycardia in the single-cell model at Cmax values of 5× and 10×, ventricular tachycardia was examined in the 3D model. Figure 6 shows the AP from the single-cell model, the AP at a point in the 3D model, and ECG data from the 3D model in the presence of dofetilide at a Cmax of 10×, which indicates the presence of ventricular tachycardia. Figure 6 also shows snapshots of ventricular AP propagation, which exhibits rotational activation. Table 6 lists the occurrences of ventricular tachycardia caused by the three drugs according to concentration.
The arrhythmia morphology was not polymorphic, which is a limitation of our model.
To test the suitability of the models to examine JT peak , we checked the rate dependence of JT peak using various models without any drug effect. All the models showed decreasing JT peak as heart rate increased with the model of Dutta et al. (2017) exhibiting the best agreement with clinical data (Johannesen et al., 2014 ; Figure 7). We also validated intercellular conduction  by comparing activation times obtained from our ventricular model with those in the literature (Supplementary Figure S1).

DISCUSSION
In evaluations of drug safety, the QT interval in ECGs has received much attention because drug-induced prolongation of the QT interval is, under certain conditions, associated with the risk of TdP, a fatal ventricular arrhythmia. However, prolongation of the QT interval does not always lead to TdP. The recently proposed CiPA initiative aims to enable more comprehensive evaluation of drug safety (Vicente et al., 2018;Wallis et al., 2018). In this study, we examined the effects of seven drugs on the QT interval using a realistic in silico 3D body-surface ECG model that included the ventricles and torso. Among the seven drugs, dofetilide resulted in the greatest increase in the QT interval, which is consistent with published data (Vicente et al., 2015). Dofetilide use entails a relatively high risk of TdP (Tisdale, 2016). However, the occurrence of TdP requires not only prolongation of the QT interval, but also early afterdepolarization (EAD) and TDR (Antzelevitch et al., 2004). In this study, dofetilide also exhibited the highest TDR values (Figure 2 and Table 4). In addition, dofetilide does not block the I NaL and I CaL channels ( Table 2), which increases the probability of EAD. The increases in APD 90 caused by ranolazine and verapamil were smaller than those induced by dofetilide, which resulted in smaller increases in the QT interval in the 3D ECG simulation. Ranolazine blocks I NaL channels almost as effectively as I Kr channels ( Table 2) and entails a low risk of TdP because blockade of I NaL channels decreases the risk of EAD (Hawwa and Menon, 2013). Interestingly, for ranolazine, the magnitude of the increase in the APD 90 was similar in endocardial and M cells, whereas for verapamil it was similar in M and epicardial cells, at 1× Cmax (Table 3). Dofetilide induced the greatest increase in the APD 90 of M cells; most QT interval-prolonging drugs increase the APD of M cells preferentially, thereby increasing the TDR value (Antzelevitch et al., 2004). Verapamil did not affect the TDR at a Cmax of 1× and decreased it at 10× Cmax, consistent with the low risk of TdP associated with its use (Milberg et al., 2005). In this study, the increases in JT peakC were similar to those in QTc, in contrast to the finding of Vicente et al. (2015) that the JT peakC is not increased by ranolazine or verapamil at a Cmax value of 1×. T peak corresponds to the time of epicardial repolarization, and most drugs that increase the epicardial APD also increase the JT peakC because of delayed epicardial layer repolarization. Ranolazine and verapamil increased the epicardial APD and the JT peakC in our simulation. Thus, the discrepancy in the JT peakC interval between our simulated results and human ECGs performed after administration of verapamil or ranolazine implies that the 3D model needs further improvement. In order to make the model accurately predict the prolongations of JT peak FIGURE 7 | Rate dependence of JTpeak. JTpeak was obtained for different heart rates with different cell models. Simulation results were compared with clinical data by Johannesen et al. (2014). and T peak -T end , the improvement of the cell models seems to be needed. If epicardial APD remains the same, and endocardial APD increases under the effects of a drug, JTpeak should remain the same, and T peak -T end should increase, which is expected in the cases of safe drugs. The current cell models do not exhibit these behaviors of APD changes under the effects of safe drugs, which disqualifies the current model as a biomarker. EAD was also observed in our simulated AP for ranolazine at a Cmax value of 10×, but not 5×. This result may reflect a limitation of the in vitro data-based simulation, in which the role of metabolites was not considered. However, simulation at a Cmax value of 10× was not recommended by the CiPA because of excessive variability in the values of the markers at higher concentrations (Li et al., 2019). Thus, the ventricular tachycardia induced by dofetilide (but not by verapamil or ranolazine) at 5× Cmax may have potential as an in silico biomarker for screening of the TdP risks posed by candidate molecules.
Although validation is needed to improve the predictive capability of the model, this study demonstrated the possibility of the model to become an effective biomarker to examine the effects of drugs on body-surface ECG parameters using realistic 3D models of the ventricles and torso. This step could lead to our ultimate goal of creating a comprehensive in silico drugsafety testing system.
This study also has several limitations. First, in the human ventricular model, we had difficulty in defining the spatial distribution of the sandwiched midcardial cell layer between the endocardial and epicardial cells. The distribution was adopted from the scheme presented in Trudel et al. (2004). Second, in the model of the ventricles and torso, we did not consider the lungs and other tissues between the body surface and the heart when solving for the body surface potentials. The bone located between the heart and body surface might influence ECG data due to its much higher electrical impedance than that of body fluids. We did not consider the effect of this bone in the ECG algorithm. In a future study, the electrical impedance of bone will be considered. Similarly, the effects of non-homogeneous properties of extracellular tissue should be incorporated into the heart model. However, we believe that these limitations did not greatly affect the major findings of this study.

DATA AVAILABILITY
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
D-SY and ES provided the main idea for this research. MH, SH, and MP obtained and analyzed the data. MH wrote the initial draft of the manuscript. D-SY and ES edited the manuscript. CL provided advice on physiological issues and contributed to the toxicity test protocol. All authors reviewed the manuscript.