ORIGINAL RESEARCH
published: 23 August 2017
doi: 10.3389/fphys.2017.00616

Optimization of an In silico Cardiac
Cell Model for Proarrhythmia Risk
Assessment
Sara Dutta, Kelly C. Chang, Kylie A. Beattie, Jiansong Sheng, Phu N. Tran, Wendy W. Wu,
Min Wu, David G. Strauss, Thomas Colatsky † and Zhihua Li *Division of Applied Regulatory Science, Office of Clinical Pharmacology, Office of Translational Sciences, Center for Drug
Evaluation and Research, U.S. Food and Drug Administration, Silver Spring, MD, United States

Edited by:
Eleonora Grandi,
University of California, Davis,
United States
Reviewed by:
Eric A. Sobie,
Icahn School of Medicine at Mount
Sinai, United States
Colleen E. Clancy,
University of California, Davis,
United States*Correspondence:
Zhihua Li
zhihua.li@fda.hhs.gov†

Present Address:
Thomas Colatsky,
Marshview Life Science Advisors,
Seabrook Island, SC, United StatesSpecialty section:
This article was submitted toComputational Physiology andMedicine,a section of the journal
Frontiers in Physiology
Received: 31 March 2017
Accepted: 09 August 2017
Published: 23 August 2017
Citation:
Dutta S, Chang KC, Beattie KA,
Sheng J, Tran PN, Wu WW, Wu M,
Strauss DG, Colatsky T and Li Z
(2017) Optimization of an In silico
Cardiac Cell Model for Proarrhythmia
Risk Assessment.
Front. Physiol. 8:616.
doi: 10.3389/fphys.2017.00616

Frontiers in Physiology | www.frontiersin.org

Drug-induced Torsade-de-Pointes (TdP) has been responsible for the withdrawal ofmany drugs from the market and is therefore of major concern to global regulatory
agencies and the pharmaceutical industry. The Comprehensive in vitro Proarrhythmia
Assay (CiPA) was proposed to improve prediction of TdP risk, using in silico models andin vitro multi-channel pharmacology data as integral parts of this initiative. Previously,
we reported that combining dynamic interactions between drugs and the rapid delayed
rectifier potassium current (IKr) with multi-channel pharmacology is important for TdP risk
classification, and we modified the original O’Hara Rudy ventricular cell mathematicalmodel to include a Markov model of IKr to represent dynamic drug-IKr interactions(IKr-dynamic ORd model). We also developed a novel metric that could separate drugs
with different TdP liabilities at high concentrations based on total electronic charge carriedby the major inward ionic currents during the action potential. In this study, we further
optimized the IKr-dynamic ORd model by refining model parameters using publishedhuman cardiomyocyte experimental data under control and drug block conditions. Using
this optimized model and manual patch clamp data, we developed an updated versionof the metric that quantifies the net electronic charge carried by major inward andoutward ionic currents during the steady state action potential, which could classify
the level of drug-induced TdP risk across a wide range of concentrations and pacing
rates. We also established a framework to quantitatively evaluate a system’s robustness
against the induction of early afterdepolarizations (EADs), and demonstrated that the
new metric is correlated with the cell’s robustness to the pro-EAD perturbation of IKr
conductance reduction. In summary, in this work we present an optimized model thatis more consistent with experimental data, an improved metric that can classify drugs
at concentrations both near and higher than clinical exposure, and a physiologicalframework to check the relationship between a metric and EAD. These findings provide
a solid foundation for using in silico models for the regulatory assessment of TdP risk
under the CiPA paradigm.
Keywords: Torsade-de-Pointes (TdP), Comprehensive in vitro Proarrhythmia Assay (CiPA), rapid delayed rectifierpotassium current (IKr), in silico cardiac cell model, drug block, proarrythmia risk, model optimization

1

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

INTRODUCTION

with multi-channel pharmacology data. This IKr-dynamic ORd
model (hereinafter referred to as the original IKr-dyn ORd
model) was calibrated based on the original ORd model sothat it can reproduce experimentally recorded adult human left
ventricular cardiomyocyte action potential (AP) morphology andrate dependency under control (drug-free) conditions. However,this model calibration process in our previous work did notinclude experimental AP changes under the influence of different
channel blocking drugs. This may negatively affect the model’s
predictive power as this model is intended for simulating drug
effects under channel blocking conditions.
In this study we further optimized the original IKr-dyn
ORd model by adjusting channel conductance values of major
ionic currents according to human ventricular cardiomyocyte
experimental data in the presence and absence of various drugs
with different channel blocking activities. We show that thisoptimization procedure allowed the model (hereinafter referred
to as optimized IKr-dyn ORd model) to quantify more accurately
the impact of each individual current on the AP. We thenscreened a series of published and novel metrics computed bythis model based on their capability of stratifying CiPA training
compounds into their corresponding TdP risk categories using
drug-IKr binding kinetics and multi-channel pharmacology datacollected earlier through manual patch clamp systems (Li et al.,2017). The best metric identified to date is based on drug-inducedchanges to the net charge carried by ionic currents (qNet) during
the AP, which can stratify the 12 CiPA training drugs into three
TdP risk levels across various conditions. We also show thatthe increased predictive power of this metric is mechanisticallylinked to the incorporation of IKr-drug binding dynamics andthe improved representation of the block effects of individual
currents, two important features of the optimized IKr-dyn ORd
model. Finally, we developed a framework to evaluate a cell’srobustness against EAD generation, and demonstrated that the
new qNet metric is correlated with the system’s repolarization
robustness to external pro-EAD perturbations that could reducethe membrane density of the hERG channel (IKr conductance).

Drug-induced Torsade-de-Pointes (TdP) is a lethal arrhythmiathat has caused removal of several drugs from the market(Gintant, 2008). The current cardiac safety paradigm (describedby the ICH E14 and S7B guidelines) focuses on two markers
to assess TdP risk: in vitro block of the hERG (human
Ether-à-go-go-Related Gene) channel (representing the rapidly
activating delayed rectifier potassium current, or IKr), andprolongation of the QTc interval in clinical studies (Sager et al.,2014). However, while eliminating the incidence of TdP inmarketed drugs, this testing regime primarily aims at detecting
delayed ventricular repolarization rather than the clinical end
point TdP, and may be assigning proarrhythmia liability to
drugs that could in fact be safe (Sager et al., 2014). Therefore,
the Comprehensive in vitro Proarrhythmia Assay (CiPA) was
proposed as a new regulatory paradigm that assesses drug TdP
risk by combining measurements of drug effects on multiple
cardiac ionic currents in vitro with in silico modeling of drug
effects on the ventricular myocyte (Sager et al., 2014). The O’Hara
Rudy cardiac cell model (ORd) (O’Hara et al., 2011) was chosen
as the consensus base in silico model and a set of 28 drugs withknown levels of clinical TdP risk (high, intermediate, low/none)
were identified for the development and evaluation of the CiPA
paradigm (Colatsky et al., 2016; Fermini et al., 2016). The three
TdP risk categories were assigned by a Clinical Translation
Working Group comprised of clinical cardiologists, safety
pharmacologists, and clinical electrophysiologists according to
published and publically available data and expert opinion. The
28 CiPA drugs were separated into a training set of 12 compounds
to be used for calibration of the in silico model and the remaining
16 compounds are to be used later for validating the model. Boththe training and validation compound sets comprise drugs thatcover the full range of TdP risk categories and demonstrate varied
electrophysiological profiles.
Previous studies have presented computational frameworks
to assess TdP risk (Mirams et al., 2011; Kramer et al., 2013;
Lancaster and Sobie, 2016), but their use within the CiPA
framework is limited due to their differing TdP risk categories
from those defined in CiPA. In addition, prior studies simulated
drug effects using the half-maximal blocking concentration(IC50) for different drugs, which assumes simple pore blockof the ion channels and neglects any intricacies of drugion channel interactions that may be important factors in
predicting relative TdP risk. The importance of incorporatinga kinetic representation of drug-ion channel interactions has
been demonstrated in previous publications (Di Veroli et al.,2013, 2014; Li et al., 2017). In the Li et al. (2017) study werecently reported the development of a novel IKr dynamic
model that can capture drug-channel dynamic interactions, andthe integration of this IKr model into the ORd cardiac model

METHODS
Optimization of the IKr-Dynamic ORd
Model
The original IKr-dyn ORd model (described in Li et al., 2017and Expanded Methods in the Supplemental Material) was
further modified (optimized IKr-dyn ORd) by scaling five ioniccurrent conductances [IKr, the slow rectifier potassium current(IKs), inwardly rectifying potassium current (IK1), the L-type
calcium current (ICaL) and the late sodium current (INaL)]
so that the model provides a good fit to published APD ratedependence experimental data for control and five channel
blockers (O’Hara et al., 2011). The optimization was performed
using the model parameterization algorithm described in Li et al.
(2014). Briefly, an initial set of scaling factors was defined withina certain range (between 0.001 and 9) and their goodness of
fit was assessed using an objective cost function defined as theweighted sum of the squared errors between model simulations

Abbreviations: CiPA, Comprehensive in vitro proarrhythmia assay; TdP, Torsadede-Pointes; ORd, O’Hara Rudy dynamic cell model (O’Hara et al., 2011); IKr-dyn
ORd, ORd model with dynamic IKr; Inet, net current (sum of currents ICaL, INaL,
IKr, IKs, IK1, Ito); qNet, charge passed by Inet from the beginning to the end of the
AP beat (same for qCaL and ICaL, qNaL and INaL...); cqInward, change in charge
passed by ICaL and INaL; Cmax, free maximum plasma clinical drug exposures.

Frontiers in Physiology | www.frontiersin.org

2

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

separation between risk categories in our previous report (Liet al., 2017). In addition, we considered a new metric (qNet)
calculated as the net charge constituting (the integral or area
under the curve of) the net current (Inet) from the beginning to
the end of the simulated beat (defined as Inet = ICaL + INaL +IKr + IKs + IK1 + Ito). The currents making up Inet within ourstudy play an important role in modulating arrhythmic risk andhave been chosen based on input from pharmaceutical company
scientists and safety pharmacology experts as the main currentsof interest within the CiPA paradigm, as outlined in Fermini et al.
(2016).
To assess the robustness of a cell against EAD generation,
we simulated an added perturbation by reducing the maximum
conductance of IKr and reporting the minimum IKr reductionneeded to trigger an EAD. Simulations were run for varying
degrees of IKr conductance reduction (using a binary searchalgorithm) at a CL of 2,000 ms with a precision of 0.01%.
For each IKr reduction tested, EADs were defined as having a
positive differential (dV/dt) during the plateau phase of the AP(between APD30 and APD90) after 100 beats. The cell model
was pulsed for a 100 beats before checking for EADs to allow thesystem to reach quasi-steady state, as in Kurata et al. (2017). The
minimum IKr conductance reduction needed to trigger an EAD
was named IKr reduction threshold, which reflects the system’s
repolarization robustness against, or distance from, EADs. To
assess the relationship between the metric and the repolarization
robustness, we calculated the correlation coefficients (using thepearson method) between the metric at steady state after 1,000
beats (without added IKr reduction) and the IKr reduction
threshold for each drug across a series of concentrations (0.5x–
25x Cmax). Situations where no IKr reduction threshold could
be calculated (no EADs could be induced for the highest IKr
reduction tested) or IKr reduction threshold is 0 (an EAD
occurred without any added perturbation) were excluded fromthe correlation calculation.

and experimental measurements. The set of scaling factors thenunderwent iterative changes (i.e., mutation and recombination)
to create new generations of parameters and this process wascontinued until the convergence criterion was met (when thechange in the minimum error of the new parameters is less than
5% over the last 30 generations). The experimental data used forfitting were taken from Figure 8 of the ORd model paper (O’Hara
et al., 2011) and comprise APD rate dependence data for control
and 5 drug blocking conditions: 1 µM E-4031 (70% IKr block),
1 µM HMR-1556 (90% IKs block), 1 µM nisoldipine (90% ICaL
block), 100 µM BaCl2 (90% IK1 block), 10 µM mexiletine (54%INaL, 9% IKr, and 20% ICaL block). The simulated percentageof block for all drugs was kept the same as in the ORd model
paper (O’Hara et al., 2011), apart from mexiletine, which used
new pharmacology data from manual patch clamp systems at
physiological temperatures (Crumb et al., 2016). The algorithm
was run using in-house developed R scripts (R Core Team, 2014)
and C programs using the Snow, Rmpi and deSolve packages(lsoda solver with a 10−6 relative and absolute tolerance) (Yu,
2002; Soetaert et al., 2010; Tierney et al., 2015) on the FDA High
Performance Computer (HPC) with 160 cores.

Simulation Protocol for Metric Evaluation
All simulations were run from control steady state conditions(after 1,000 beats) at varying cycle lengths (CLs) 1,000, 2,000,
and 4,000 ms and stimulus of −80 µA/µF for 0.5 ms (as in theoriginal model). Block of ion channels at various concentrations
were simulated and run for another 1,000 beats to reach a
new steady state with drug. The last two beats were analyzed
to check for alternans, which was observed in the presence ofearly afterdepolarizations (EADs), defined as having a positive
derivative during the repolarization phase of the AP. The
pharmacology data for the 12 CiPA training compounds (thefull list and their corresponding risk categories can be found in
Supplemental Table 1) were the same as in our previous report(Li et al., 2017), where drug-IKr binding kinetic parameters were
estimated using an in vitro IKr dynamic protocol and IC50/Hillcoefficients based on Crumb et al. (2016) were used for theremaining channels [the peak sodium current (INa), INaL, ICaL,IK1, IKs and transient outward potassium current (Ito); all theparameters can be found in the Supplemental Tables 2 and 3].
Simulations were run for a range of drug concentrations: from
0.5x up to 25x free maximum plasma clinical drug exposures(Cmax). Simulations were run in R and C using the deSolvepackage (Soetaert et al., 2010).
We assessed a range of standard metrics as also consideredin Mirams et al. (2011), Lancaster and Sobie (2016): resting
membrane potential (resting Vm), maximum upstroke velocity(dV/dtmax), peak membrane potential (peak Vm), APD at
50% of the amplitude (APD50), APD at 90% of the amplitude(APD90), APD triangulation (APDtri) defined as APD90APD50, diastolic intracellular calcium concentration ([Ca2+ ]i )(diastolic Ca), peak [Ca2+ ]i (peak Ca), calcium transient duration
at 50% (CaD50) and 90% (CaD90) of the amplitude, calcium
transient triangulation (Catri) defined as CaD90-CaD50, as wellas the cqInward metric that quantifies the change in the amountof charge carried by INaL and ICaL, which demonstrated good

Frontiers in Physiology | www.frontiersin.org

Classification Methods
To assess the ability of the metrics to identify each drug’s
TdP risk level, we performed a proportional odds logistic
regression classification and a leave-one-out validation, as in
Mirams et al. (2011). If EADs were observed, the metric value
at the concentration prior to EAD generation was used for theclassification. We used the R lrm function from the rms package(https://CRAN.R-project.org/package=rms) and calculated theclassification training error for each metric as follows: the mean(across 12 drugs) of the absolute error (difference between
predicted and known risk categories), with the risk categories
defined as 1 = low risk, 2 = intermediate risk and 3 = high risk.
The proportional odds logistic regression model is a regression
model for ordinal dependent variables, and accounts for theordering by using cumulative probabilities defined as the oddsof (Y ≤ i) = P(Y ≤ i)/(1 − P(Y ≤ i)) for each risk categoryi, where Y is the variable that represents a drug’s risk category
and P(X) is the probability of X. The function uses maximum
likelihood estimates to calculate the probability of each drug
being a member of each risk category, and the drug is assigned
to the risk category corresponding to its highest probability. We

3August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

then performed a leave-one-out validation by removing one drug
from the data set and then predicting its risk category based on
the classification of the remaining drugs. This was performed
in turn for each drug within the data set and its leave-one-out
prediction error was calculated the same way as the training
error.

model compared to the original IKr-dyn ORd model (fitting
error 30.81 vs. 72.33). The main improvements in the qualityof fit to the experimental data are observed for drug blocking
conditions, especially with mexiletine (INaL blocker) and E4031 (IKr blocker) (Figures 1B,E respectively). In the case of
mexiletine, a reduction in the fitting error from 91.09 (original
IKr-dyn ORd model) to 18.36 (optimized IKr-dyn ORd model)
was achieved (Table 2). Figure 1B shows that, with the original
IKr-dyn ORd model, the simulated APD prolongation with
mexiletine is significantly longer than experimental data. A
similar pattern can be seen for the IKr blocker E-4031 (Table 2
and Figure 1E). Due to the opposite roles of INaL and IKr in
prolonging AP (Johannesen et al., 2016), this suggests that blockof INaL is underestimated and that of IKr is overestimated in
the original IKr-dyn ORd model. The optimized model corrected
these inaccuracies with better fitting to the experimental data,
which is important for TdP risk assessment as it is known
that INaL block plays an important role in counteracting proarrhythmic APD prolongation of IKr block (Orth et al., 2006;
Johannesen et al., 2016).
To further understand the contribution of various ioniccurrents to AP profile after the optimization process, we
compared the simulated AP traces and the ionic currents
over the time course of the steady state AP between theoriginal and our optimized IKr-dyn ORd model at different
cycle lengths. As described earlier in this section all of thecurrent conductances are increased in the optimized IKr-dyn
ORd model (Table 1). However, the AP shapes from bothmodels under control conditions are very similar, as shown in
Figure 2A. This is consistent with the fact that both models fit
the control AP morphology parameters (Figure 1A) reasonably
well. On the other hand, while only a small change in current
amplitude is observed for ICaL (Figure 2C), which only has
a 0.7% change in conductance (Table 1), clear differences areobserved for all other currents (IKr, INaL, IKs and IK1) withthe biggest changes occurring for INaL (conductance is increased
by 166.1% between the optimized and original models as shownin Table 1). This further demonstrates that INaL plays a bigger
role in the optimized model than the original IKr-dyn ORd
model.

RESULTS
Optimized IKr-dyn ORd Model
The optimized IKr-dyn ORd model was built by scaling the
conductance of the main ion currents (IKr, IK1, IKs, INaL, ICaL)of the original IKr-dyn ORd model (presented in Li et al., 2017)
to fit the APD rate dependence experimental data in control
and drug block conditions from O’Hara et al. (2011). The setof scaling factors that gives the best fit for the optimized modelis as follows: scaling IKr by 1.013, IKs by 1.870, IK1 by 1.698,
ICaL by 1.007, and INaL by 2.661, as summarized in Table 1.
A comparison of the simulation results from both the original
and optimized models to the experimental data is shown in
Figure 1 and the corresponding sum of squares errors (between
simulation and experimental data) are shown in Table 2. Sum
of squares error for the original ORd model as presented intheir paper (O’Hara et al., 2011) are also shown in Table 2 forcomparison purposes. We see that although for control and somecurrent blocking conditions the original IKr-dyn ORd model has
errors similar to the original ORd model, for other conditions
the errors were worsened (IK1 and ICaL blocking experiments),
resulting in an average error bigger than the original ORd model(72.33 vs. 57.77). However, the discrepancy between simulations
and experiments was significantly reduced in the new optimized
IKr-dyn ORd model.
As can be seen in Figure 1A, under control conditions boththe original IKr-dyn ORd and optimized IKr-dyn ORd models
display similar behavior. Although for control data points, theoptimized IKr-dyn ORd model fitting is slightly worse thanthe original IKr-dyn ORd model (fitting error 22.63 vs. 18.82
in Table 2), the average fit across both control and all drug
block conditions is much better for the optimized IKr-dyn ORd

Candidate Metrics
We then investigated whether the optimized IKr-dyn ORd model
could be used to stratify proarrhythmia risk levels. As a firststep we explored the changes in AP and individual currents
induced by three representative drugs (one taken from eachone of the CiPA TdP risk categories), using pharmacology data
as used in Li et al. (2017). Since it is known that the subtle
balance between inward (such as INaL and ICaL) and outward(such as IKr, IKs, IK1, and Ito) currents underlies the generationof EADs, a mechanistic precursor to TdP (Vos et al., 1995;
Volders et al., 2000; Wu et al., 2002; Weiss et al., 2010), we
also examined the net current between inward and outward
currents (Inet) in addition to individual currents. Figure 3 shows
simulations of AP, Inet, ICaL, INaL, IKr, IKs, IK1, and Ito forranolazine (low risk), cisapride (intermediate risk) and dofetilide(high risk), for a CL of 2,000 ms and a dose of 25x Cmax

TABLE 1 | Conductance scaling factors for the original and optimized IKr-dynamic
O’Hara-Rudy models (original and optimized IKr-dyn ORd): current conductancesof the rapid (IKr) and slow (IKs) rectifier potassium current, inwardly rectifying
potassium current (IK1), L-type calcium current (ICaL) and late sodium current(INaL) in the model are multiplied by the corresponding scaling factor.
Scaled currents

Original IKr-dyn
ORd model

Optimized IKr-dyn
ORd model

IKr

0.9

1.013

IKs

1.0

1.870

IK1

1.0

1.698

ICaL

1.0

1.007

INaL

1.0

2.661

Note that the IKr conductance in the original IKr-dyn ORd model was scaled as describedin Li et al. (2017).

Frontiers in Physiology | www.frontiersin.org

4August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

FIGURE 1 | Steady state action potential duration (APD) rate dependency under the conditions of control (A), 10 µM mexiletine [late sodium current INaL block (B)],
1 µM nisoldipine [L-type calcium current ICaL block (C)], 1 µM HMR-1556 [slow rectifier potassium current IKs block (D)], 1 µM E-4031 [rapid rectifier potassium
current IKr block (E)] and 100 µM BACl2 [inwardly rectifying potassium current IK1 block (F)] at varying cycle lengths (CLs) for the original dynamic IKr O’Hara Rudy
model (original IKr-dyn ORd; dashed lines) and the optimized dynamic IKr O’Hara Rudy model (optimized IKr-dyn ORd model; solid line). Experimental data mean
(symbol) and standard deviation (error bars) are from O’Hara et al. (2011). Control (A) shows APD at 90% (APD90; filled circles), 70% (APD70; filled triangles), 50%(APD50; filled squares) and 30% (APD30; plus sign) repolarization. All other panels show APD90.

concentration (25x Cmax) is used to highlight the potential
differences between various risk levels. The amount of electronic
charge carried by each current is calculated as the area under the
curve (AUC) of the individual current trace and is plotted for Inetin Figure 3C.
We see in Figure 3A that all three drugs cause prolongationof APD and the low risk drug, ranolazine, shows a greater
prolongation of APD compared to the intermediate risk drug,
cisapride (266.78 vs. 176 ms). The performance of APD90 asa metric for all the drugs from 0.5 to 25x Cmax, can be seenin Supplemental Figure 1. In fact, verapamil and ranolazine(both low risk) display APDs greater than most intermediaterisk drugs over a wide range of doses. Therefore, the amountof APD prolongation is not a good indicator of the TdP riskof a drug, demonstrating the unsuitability of APD alone as amarker for TdP risk. However, we notice that Inet (Figure 3B),
calculated as the sum of the five main currents that modulate
the plateau phase of the action potential (ICaL, INaL, IK1, IKr,
IKs, and Ito, shown in Figures 3D–I), does correlate with the
TdP risk category. As shown in Figure 3C, the order of qNet
(charge carried by Inet integrated from the beginning to the endof the AP beat) is consistent with the rank order of TdP risk
levels for the three drugs. At the end of the CL, ranolazine has
a qNet of 0.061 µC/µF while cisapride and dofetilide have aqNet of 0.037 µC/µF and 0.013 µC/µF, respectively. A detailed

TABLE 2 | Sum of squares error (divided by 100) between experimental action
potential duration (APD) rate dependence mean data (from Figure 8 in O’Hara
et al., 2011) and the original O’Hara Rudy model (original ORd) (O’Hara et al.,2011), the original IKr-dyn ORd (Li et al., 2017) as well as the optimized IKr-dyn
ORd under different conditions: control, mexiletine (blocks mainly INaL), HMR
1556 (blocks IKs), E-4031 (blocks IKr), BaCl2 (blocks IK1) and nisoldipine (blocks
ICaL).
Experiment

Sum of squares error
Original ORd

Original IKr-dyn
ORd

Optimized IKr-dyn
ORd

Control

17.20

18.82

22.63

Mexiletine (INaL)

92.92*

91.09

18.36

HMR 1556 (IKs)

56.35

57.08

55.34

E4031 (IKr)

145.03144.8772.33

Bacl2 (IK1)

29.83

67.47

11.41

Nisoldipine (ICaL)

5.29

54.62

4.76

Average

57.77

72.33

30.81

*Error was calculated using the updated mexiletine IC50 data (Crumb et al., 2016); using
the block suggested in the ORd paper of 90% INaL block (O’Hara et al., 2011), the sumof squares error is of 38.48, changing the average error to 48.70.

using our optimized model. A slow pacing rate (CL 2,000 ms)
is used here because bradycardia is a known risk factor for
TdP (Kurita et al., 1992; Kallergis et al., 2012), and a highFrontiers in Physiology | www.frontiersin.org

5August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

FIGURE 2 | Action potential (AP) (A), INaL (B), ICaL (C), IKr (D), IKs (E), IK1 (F) traces under control conditions for the original IKr-dyn ORd (dashed line) and the
optimized IKr-dyn ORd (solid line) for CLs of 500 (red), 1,000 (green) and 2,000 (blue) ms.

metric over a range of doses (0.5–25x Cmax) and a range of CLs(1,000, 2,000, and 4,000 ms) for the 12 CiPA training compounds.
This error quantifies the mean (across the 12 CiPA drugs)
difference between known and predicted risk levels for each
metric. We can see that across the full range of concentrations
and all CLs the qNet metric shows the smallest classification
training error. Notably, the qNet metric shows a classification
training error of 0 for concentrations greater than or equal to1x Cmax, meaning it consistently classifies each of the 12 CiPA
training compounds into the correct TdP risk category. The
cqInward metric performance is comparable to that of qNet at
low pacing rates (4,000 ms) and high drug concentrations. All of
the other standard metrics we considered show training errorsthat never come down to 0, which fluctuate across the range ofdoses.
The results presented in Figure 4 are consistent with the leaveone-out validation described in Table 3 performed on a subsetof the doses tested (1, 10, and 20x Cmax) for a CL of 2,000
ms; the cqInward and qNet show the smallest prediction errors
with values of 0.33 and 0.08 respectively at 20x Cmax. The othernext best performing metrics are peak Vm with an error of
0.42 and APD50, APD90 with errors of 0.5 at Cmax 20x. Of
note, at 1x Cmax, qNet and APD90 all have the same prediction
error of 0.17. This is because at lower concentrations (1x Cmax)
the effects of each drug are harder to differentiate due to there
often being only subtle effects on the AP morphology. However,
the CiPA paradigm assumes that the assessment of TdP risk
may occur at any time during drug discovery and development,even prior to the time the clinical effective drug concentrations
are known with any certainty. In addition, the incidence ofclinical TdP is limited and not necessarily related strictly tonormal (1x) clinical exposure (i.e., concomitant factors may playa role in expressing clinical TdP events). Therefore, we proposethat a metric should be evaluated under multiple physiologicaland pharmacological conditions. The overall evidence suggest

examination of the individual current profiles reveals whyranolazine caused the least amount of qNet decrease. As shown in
Figures 3D, G, ranolazine (green lines) caused a marked decreaseof the absolute amount of charge carried by IKr (qKr decrease of
0.119 µC/µF) and INaL (qNaL decrease of 0.07 µC/µF) at theend of the AP beat compared to control (black lines). Because
the outward current IKr and inward current INaL have opposite
directions, ranolazine-induced reduction (in absolute values) of
the two currents balanced each other and resulted in only asmall change of the net charge at the end of the AP (qNet,
Figure 3C). In contrast, dofetilide (Figure 3D, red lines) andcisapride (Figure 3D, blue lines) caused a significant reduction ofqKr (0.135 and 0.063 µC/µF respectively) through direct channel
blocking, and a slight increase of qNaL through prolonged APD.
These two effects changed Inet in the same direction and worked
together to decrease qNet significantly, with dofetilide causing
the biggest decrease due to more significant blocking of IKr. Notethat these drugs have some effects on other currents (Ito, IKs, andIK1) as well, but those changes are relatively small and will notchange the rank order of qNet values significantly for the threedrugs tested here. However, these other currents may become
important for drugs that directly block them. For example, theeffects on ICaL may be critical in determining the qNet change
and risk level for a calcium blocker.
These initial promising results prompted us to calculate thisnew Inet-based metric, qNet, for all 12 CiPA training compounds
and systematically compare its capability of separating the three
TdP risk levels to a range of commonly tested metrics (describedin the Methods section). The risk categories, IC50 and IKr
dynamic parameters for each drug are listed in Supplemental
Tables 1–3. Included in the comparison is also the cqInward
metric, described in our previous study and defined as thenormalized drug-induced change of the charge carried by the
inward currents INaL and ICaL (Li et al., 2017). As shown in
Figure 4, we calculated the classification training error for each

Frontiers in Physiology | www.frontiersin.org

6

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

FIGURE 3 | Transmembrane voltage [Trans. voltage (A)], net current [Inet (B)], charge carried by Inet (C), INaL (D), Ito (E), ICaL (F), IKr (G), IKs (H), IK1 (I) traces forcontrol (black solid line), ranolazine, a low TdP risk drug (green dashed line narrow spacing), cisapride, an intermediate TdP risk drug (blue dashed line normal
spacing), and dofetilide, a high TdP risk drug (red dashed line wide spacing), at 25x Cmax for 2,000 ms CL using the optimized IKr-dyn ORd. Charge carried by Inet,
INaL and IKr integrated from the beginning to the end of the AP beat (qX) are displayed on the graph.

to shed light on possible mechanistic differences among the drugs
tested, we used the best candidate metric qNet as a benchmark,and compared the performance of the optimized IKr-dyn ORd
model with model variations where each of the changes wasremoved in turn. Figure 5 shows computed qNet values for the
12 CiPA training drugs calculated over a range of drug doses
from 0.5x to 25x Cmax when using the optimized IKr-dyn ORd
model (Figure 5A), a model variation without incorporating
the IKr dynamic model (Figure 5B) and a model variation
incorporating the IKr dynamic model but without optimizing
channel conductances (Figure 5C). In line with results from
Figure 4 and Table 3, the metric qNet shows clear separation
between the 3 TdP risk categories across the range of doses
tested with the optimized IKr-dyn ORd model (Figure 5A);
however, this is not the case for the other two model variations(Figures 5B,C).

that qNet is the best among all the metrics tested, because it
has a training error of 0 across a wide range concentrations(1–25x Cmax) at various pacing frequencies (2,000 and 4,000ms), and the lowest leave-one-out error at all concentrations
tested.

The Impact of IKr-Drug Binding Kinetics
and Channel Conductance Optimization on
Risk Level Stratification
Compared to the original ORd (i.e., the consensus base model
for CiPA), the optimized IKr-dyn ORd model presented in thiswork has two important changes: the incorporation of a dynamic
IKr model to capture drug binding kinetics (Li et al., 2017), andan improved set of channel conductances to better represent thecontribution of individual currents to AP (Figures 1, 2). In order

Frontiers in Physiology | www.frontiersin.org

7

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

FIGURE 4 | Classification training error for a range of metrics [resting membrane potential (resting Vm), maximum upstroke velocity (dV/dtmax), peak membranepotential (peak Vm), APD50, APD90, action potential (AP) triangulation (APDtri), diastolic intracellular calcium concentration ([Ca2+ ]i ) (diastolic Ca), peak [Ca2+ ]i (peakCa), calcium transient duration at 50 and 90% of the amplitude (CaD50 and CaD90), calcium transient triangulation (Catri), change in amount of charge carried by
INaL and ICaL (cqInward) and the charge carried by Inet at the end of the AP beat normalized to control (qNet)] for varying drug doses (0.5–25x Cmax) and varying
CLs (1,000, 2,000, and 4,000 ms). Each box represents the mean (across 12 drugs) error (between predicted and known risk levels) for each metric at eachconcentration (0.5–25X Cmax). A training error of 0 represents perfect separation between the risk categories.

2.6 nM see Supplementary Material), with a safety margin (IKr
IC50/Cmax) of 3.8 (Redfern et al., 2003), which is close to thatof the high risk drug dofetilide (IC50 4.87 nM and Cmax 2
nM, safety margin 2.4) for example. So if IC50 data are usedwith an assumption of simple pore drug block, cisapride isgrouped with the high risk drugs. However, when we consider
the IKr-drug binding dynamic data (Li et al., 2017), cisapride,
but not high risk drugs like dofetilide, can rapidly dissociate fromthe hERG channel during diastolic intervals because it is nottrapped in the closed channel. Consequently, cisapride has an
actual block potency lower than high risk drugs despite similar
IKr IC50/Cmax ratio, which may explain why it belongs to theintermediate rather than high risk level. On the other hand,

The first model variation we tested does not have the IKr
dynamic model incorporated but instead uses simple IC50s to
represent channel block (Figure 5B). Note that this model variant
has gone through a channel conductance optimization processsimilar to that presented in this article, as described in Dutta
et al. (2016), so the difference observed between this model
variant (Figure 5B) and the full optimized IKr-dyn ORd model(Figure 5A) is mainly due to the different representation of IKr
block (dynamic vs. IC50s). From Figure 5B we can see that
there are two intermediate risk drugs that are not correctly
categorized: cisapride that is mixed with the high risk drugs, andchlorpromazine that is mixed with the low risk drugs. Cisapride
is a potent and selective IKr blocker (IC50 10.1 nM and Cmax

Frontiers in Physiology | www.frontiersin.org

8

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

than low risk drug. This demonstrates that including a dynamic
representation of IKr-ion channel interactions is important for
categorizing TdP risk of drugs and IC50 data alone are notsufficient.
The second model variation we tested has the IKr dynamic
model included, but without optimized channel conductances
to reproduce AP changes under channel blocking conditions(Figure 5C). Note that this model variant is the same as
the original IKr-dyn ORd model (Li et al., 2017) and, as
demonstrated in Figure 1, has an inaccurate quantification of theblock effects of individual currents compared to experimental
data. In this scenario the low risk drug ranolazine is misclassified
as a high risk compound (Figure 5C). Ranolazine is a potent IKr
and INaL current blocker and these two effects can balance eachother to reduce ranolazine’s TdP risk (Antzelevitch et al., 2004;
Johannesen et al., 2016; Saad et al., 2016). Because the INaL effectis underestimated and the IKr effect is overestimated without
channel conductance optimization (Figure 1), ranolazine has
a dominant IKr block effects when simulated by this model
variant and thus will be mistakenly put in the high risk category(Figure 5C). Taken together, this suggests that the two added
features are both important for TdP risk stratification and may
mechanistically explain why a certain drug belongs to a specificTdP risk level.

chlorpromazine is not a potent IKr blocker (safety margin 24.4,
similar to other low risk drugs) so when we look at IC50 only it
is classified closer to the low risk drugs. But when IKr dynamic
data are considered, chlorpromazine is highly trapped in theclosed hERG channel and very slow in unbinding during diastolic
intervals (Li et al., 2017). This makes it more dangerous than its
IKr IC50 suggests and thus classified as an intermediate rather

TABLE 3 | Leave-one-out prediction error for a range of metrics at a CL of 2,000
ms and 3 doses (1, 10, and 20x Cmax): resting membrane potential (resting Vm),
maximum upstroke velocity (dV/dtmax), peak membrane potential (peak Vm),
APD at 50 and 90% of the amplitude (APD50 and APD90), action potential (AP)triangulation (APDtri), diastolic intracellular calcium concentration ([Ca2+ ]i )(diastolic Ca), peak [Ca2+ ]i (peak Ca), calcium transient duration at 50 and 90% of
the amplitude (CaD50 and CaD90), calcium transient triangulation (Catri), change
in amount of charge carried by INaL and ICaL (cqInward) (Li et al., 2017) andcharge carried by the net current (qNet).
MetricLeave-one out prediction error
1x Cmax

10x Cmax

20x Cmax

qNet

0.17

0.08

0.08

cqInward

0.25

0.33

0.33

Catri

1.25

1.08

1.08

CaD90

1.42

1.42

0.83

CaD50

1.42

1.0

0.83

Physiological Significance of qNet

peak Ca

0.92

0.75

0.83

diastolic Ca

0.92

0.83

0.83

APDtri

0.5

0.67

0.58

APD90

0.17

0.5

0.5

APD50

0.33

0.33

0.5

peak Vm

0.42

0.42

0.42

In order to assess the physiological significance of the metric,
we borrowed some concepts from non-linear dynamic theory,
where EADs appear as membrane voltage oscillations whenthe equilibrium state at the plateau phase (membrane voltage
between 0 and −40 mV) changes its stability via bifurcation (Quet al., 2013; Kurata et al., 2017). The robustness of the system
could be evaluated by applying a specific perturbation with aseries of strengths and measuring the range of the perturbation

dV/dtmax

0.42

0.83

1.17

resting Vm

1.00

0.75

0.75

FIGURE 5 | qNet for the 12 CiPA training compounds for a range of doses (0.5–25x Cmax) at a pacing rate of 2,000 ms. (A) Optimized IKr-dyn ORd; (B) A model
variation without the incorporation of the IKr dynamic model (note this is the same model as in Dutta et al., 2016) and; (C) A model variation without the optimized
channel conductances to accurately quantify block effects of individual currents (note this is the same model as in Li et al., 2017). Different TdP risk levels are color
coded (high risk in red, intermediate risk in blue and low/no risk in green). Results are not shown once drug concentrations are high enough to induce early after
depolarizations (EADs) (i.e., quinidine).

Frontiers in Physiology | www.frontiersin.org

9August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

the system can tolerate without changing stability (i.e., emergenceor annihilation of oscillations) (Kurata et al., 2008). We applied
this concept to our model using IKr maximum conductance
reduction as a perturbation. In this case the minimum IKr
reduction required to induce an EAD (IKr reduction threshold)reflects the system’s robustness against, or distance from, EADs.
Therefore, for each drug over a range of concentrations from
0.5 to 25x Cmax we calculated the IKr reduction thresholds, andchecked their correlation with the metrics qNet, APD90, andcqInward respectively. Detailed correlation plots for each metriccan be found in the Supplemental Figures 2–4. Table 4 shows thecorrelation coefficients for each drug across all concentrations
for IKr reduction threshold vs. qNet, APD90 and cqInward
respectively. We see that qNet shows a strong correlation acrossall drugs (close to 1). As qNet increases the IKr reduction
threshold (and the system’s robustness against EAD) increases
and vice versa as qNet decreases. The bigger the qNet value thesafer the system is and the harder it is to induce EAD.
For APD90, in most cases there is a strong negative correlation
with IKr reduction threshold (close to −1) as expected, indicating
the longer the APD the lower the repolarization robustness(i.e., the closer to EAD) and vice versa (Table 4). However, thistrend reverses completely for some drugs like verapamil andmexiletine, where the correlation is positive (Table 4), suggesting
the longer the APD90 the higher the repolarization robustness(the further away from EAD). This is contradictory to thegeneral perception that longer APD90 (and QTc) signals a higher
EAD/TdP liability. These unexpected relationships between APD
and EAD can be seen more clearly in Figure 6, where the APtraces before and after the perturbation are shown. As can beseen from Figure 6A (left panel), using APD90 as a metric a cellunder mexiletine at 1x Cmax seems safer (APD less prolonged)
than at 10x Cmax, while qNet suggests otherwise (1x Cmax
is more dangerous due to a smaller qNet value). When thesame perturbation was applied (95% IKr reduction), the cellwith 1x Cmax of mexiletine but not 10x, generated an EAD(Figure 6A right panel), indicating the cell with lower mexiletine
concentration is actually closer to EAD generation, consistent
with the prediction of qNet but not APD90. The same pattern
can be seen in Figure 6B, where verapamil at 1x Cmax is shown
to be closer to EAD than at 3x Cmax through perturbation assays(right panel), contradictory to the prediction using APD90 but
not qNet (left panel). This pattern holds true when comparing
ranolazine and cisapride as compared in Figure 3. As describedearlier, a cell under ranolazine has a longer APD90 (indicating
higher risk) and also a higher qNet value (indicating lower risk)
than cisapride at 25x Cmax (Figure 6C left panel). An added
perturbation of 75% IKr reduction will trigger an EAD withcisapride but not ranolazine (Figure 6C right panel), supporting
the prediction of qNet but not APD90. Note that here we used
25x Cmax to match the concentrations used in Figure 3. When 1x
Cmax was used, the same pattern was seen for the two drugs (seeSupplemental Figure 5). This suggests under most circumstances
qNet is a better metric than APD90 in marking the repolarization
robustness to added perturbation of IKr reduction.
Finally, cqInward does not correlate well with robustness
against EAD generation, measured as IKr reduction threshold

Frontiers in Physiology | www.frontiersin.org

TABLE 4 | Correlation (using pearson method) between qNet, APD90 andcqInward and IKr reduction threshold for the 11 drugs (diltiazem is not included
because EADs could not be induced for the highest IKr reduction tested 99.99%)
for a CL of 2,000 ms across all doses from 0.5 to 25x Cmax. Simulations where
the IKr reduction threshold is 0 (EADs occur without added IKr reduction, as forquinidine ≥ 2.3x Cmax) and results where the IKr reduction threshold could not becalculated (the maximum IKr reduction tested, 99.99%, did not trigger an EAD, as
for diltiazem at all Cmax, verapamil ≥ 1.7x Cmax, and mexiletine ≥ 3.8x Cmax)
were excluded.
Drug

Correlation with IKr reduction threshold
qNet

APD90

cqInward

Quinidine

0.996−0.994−0.197

Bepridil

0.948−0.992

0.432

Sotalol

0.979−0.992−0.971

Dofetilide

0.96−0.993−1

Cisapride

0.988−0.996−0.994

Ondansetron

0.997−0.999−0.595

Terfenadine

0.968−0.9440.804

Chlorpromazine

0.995

−1

0.895

Ranolazine

0.87−0.992

0.961

Verapamil

0.977

0.998

−0.991

Mexiletine

0.974

0.989

−0.983(Table 4), despite a good performance (next to only qNet)on separating the risk categories for the training compounds(Figure 4). This suggests cqInward does not indicate therepolarization robustness to a perturbation of hERG channel
density decrease. Whether cqInward is correlated with therobustness to another perturbation, or its separating power on
the 12 training drugs is a non-physiological artifact, remains to beinvestigated. If the latter this highlights the importance to assessa metric in not only a pre-defined drug classification system,
but also a physiological framework to quantitatively evaluate thecorrelation between the metric and EAD.

DISCUSSION
In this study we present an optimized version of the ORdmodel (O’Hara et al., 2011), which incorporates a dynamic
representation of IKr to allow modeling of drug-IKr channelinteractions (Li et al., 2017) as well as providing a better fit toexperimental data in both control and drug blocking conditions
by rescaling ionic current conductances. Most notably, INaL
current is increased compared to the original model. We also
demonstrate that our optimized model, used in combination witha mechanistic net charge metric (qNet), enables good separationof 12 CiPA training compounds into their respective risk
categories over a range of drug concentrations and pacing rates.
Furthermore, we show that this is because qNet is correlated witha system’s repolarization robustness to external perturbation ofhERG channel density decrease, or IKr maximum conductance
reduction.

Optimization of the O’Hara Rudy Model
To optimize the model we rescaled ionic current conductancesin the model presented by Li et al. (2017). We demonstrate in

10

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment(Figure 5). This highlights the need for more precise model
representation to simulate drug effects and stratify TdP risk
levels. Additional human cardiomyocyte data may help to further
refine this model.

Figure 2 that recalibration of the model conductances leads tovery little changes in the AP in control conditions across a rangeof pacing frequencies. It did however shift the effect of someof the different currents in the cardiac cell model, most notably
increasing the role of INaL in contributing to AP. A recent study
also optimized the ORd model to fit various LQTS profiles (Mannet al., 2016). Mann et al. present an optimized version of theORd model by scaling the conductances of IKs (by 5.75), IKr (by1.00), ICaL (by 2.01), INaL (by 1.00), INaCa (by 2.95), and INaK(by 9.12). The scaling factors are different to the ones observedin our model: IKs is increased in both models although ours
is increased by a smaller amount 1.87, in their model IKr isunchanged while it is slightly increased in our model, ICaL is
increased significantly in their model but only very slightly inours and INaL is unchanged in their model while it is scaled by
2.661 in our model. The differences in INaL can be explainedby the differences in context of use of the model: Mann et al.
investigate the effects of increased INaL (LQTS3) as opposed todrug block of INaL, as in this study. Furthermore, a key difference
between our model optimization process and Mann et al. is that
we used human cardiomyocyte experimental data with various
channel blockers, while they used clinical LQTS data. However,
one of their findings was that the ORd model over predicts theeffect of IKr block (50% IKr block produced a 42% increase in
APD90 as opposed to the 16.5% observed clinically), which is
concurrent with our findings. An awareness of this property of
the ORd model is important as the model is often considereda consensus gold standard model for simulating drug effects
on cardiac cells, and properties such as the over prediction ofblock of IKr may lead to inaccurate predictions of drug effects
on cardiac electrophysiology. Our manuscript further highlights
this point and provides an alternative model with improved
balance of the effect of the different ionic currents in drug block
conditions.

qNet Correlates with the System’s
Robustness against EADS
Based on ideas from non-linear dynamic theory and studies
demonstrating mechanisms of EAD generation (Guo et al., 2007;
Weiss et al., 2010; Xie et al., 2010; Chang et al., 2013; Kurata et al.,2017), we established a theoretical framework to quantitatively
evaluate the physiological consequences of the change of theqNet (and in principle any) metric. A key concept here is thesystem robustness (Kurata et al., 2008), which is defined as thelevel of a specific perturbation the system can tolerate without aqualitative change of stability (e.g., emergence or annihilation ofoscillations). We applied that concept here using IKr maximum
conductance reduction as a perturbation. Note that in ourmodel all drugs’ hERG/IKr block is modeled as binding todifferent channel states without changing the IKr conductance.
Thus the IKr conductance decrease applied here reflects extra
pro-EAD perturbations independent of each drug’s direct ion
channel block activities, for example inter-subject variability(hERG channel density variation due to genetic background),intra-subject variability (regional difference in hERG channel
density), chronic drug effects (to block hERG maturation), or
drug-drug interaction. We found that qNet is correlated withthe cell’s repolarizing robustness to the perturbation of IKr
conductance reduction. When qNet increases, the cell’s IKr
reduction threshold also increases, meaning the cell is moving
away from EAD and needs a more severe perturbation of IKr
conductance reduction to trigger an EAD. The opposite happens
when qNet decreases. This positive correlation is consistentacross all the compounds tested in this study. In contrast, APD90
does not show a consistent correlation with the repolarization
robustness across all the drugs, suggesting for some drugs(mainly compounds with balanced inward and outward current
blocking activities) APD90 may not be a good indicator of
distance from EAD.
The concept of robustness to pro-EAD perturbations is highly
related to that of repolarization reserve, developed by Roden(1998) to describe the redundant cellular mechanisms to effectorderly and rapid repolarization, which can be disrupted by anadded stressor (perturbation), resulting in APD prolongation
and/or EAD. We chose to use the term robustness insteadof repolarization reserve because the latter has been widelyused to describe a cell’s repolarization mechanism against
both delayed repolarization (APD prolongation) and voltage
oscillation (EAD), which we show in Figure 6 are not necessarily
correlated with each other. In contrast, robustness of a system,a concept borrowed from non-linear dynamic theory (Kurata
et al., 2008), is directly related to emergence or annihilationof oscillations (EADs) in the presence of perturbations. There
are different types of perturbations that could be used to testa system’s robustness, for instance an applied bias current(Gray and Huelsing, 2008; Kurata et al., 2008), or an increased
conductance for ICaL and/or INaL. We chose IKr conductance

Performance of the qNet Metric Using the
Optimized Model
Using pharmacology data for the 12 CiPA training compounds(Li et al., 2017), we assessed the suitability of a range of standard
metrics based on AP morphology properties, as well as therecently published cqInward metric (Li et al., 2017) and our
new qNet metric. We demonstrated that the commonly usedAP-based metrics are poor indicators of TdP risk and foundthat our qNet metric allowed best separation of the CiPA
training compounds into their risk categories. Our new metricoutperformed the cqInward metric presented in Li et al. which
may be a consequence of optimized channel conductances tobetter quantify the block effects of individual currents.
Our optimized IKr-dyn ORd model has two importantfeatures compared to the original ORd model: incorporationof modeling drug-IKr interaction kinetics based on dynamichERG binding data (Li et al., 2017) and better characterizationof individual currents’ role in AP based on channel blocking
data. We demonstrate the importance of simulating drug-IKr
dynamics and accurate drug block conditions by providing
rationale for misclassification of compounds when either oneof the features were removed during TdP risk classification

Frontiers in Physiology | www.frontiersin.org

11

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

FIGURE 6 | AP traces for mexiletine (A) at 1x Cmax (black solid line) and 10x Cmax (gray dashed line) without (left panel) and with 95% IKr reduction (right panel);verapamil (B) at 1x Cmax (black solid line) and 3x Cmax (gray dashed line) without (left panel) and with 98% IKr reduction (right panel); and ranolazine (black solid line)
and cisapride (dashed gray line) (C) at 25x Cmax without (left panel) and with 75% IKr reduction (right panel) for a CL of 2,000 ms. Corresponding APD90 (ms) and
qNet (µC/µF) values are reported in black for mexiletine 1x Cmax, verapamil 1x Cmax and ranolazine 25x Cmax and in gray for mexiletine 10x Cmax, verapamil 3x
Cmax and cisapride 25x Cmax. Note the IKr reduction (simulated by scaling the IKr maximum conductance) is applied in addition to the drug block effect and is used
to assess the system’s robustness against EADs (see Results section).

reduction as a perturbation because it is independent of thedirect drug effects (the dynamic IKr model allows us to model
IKr blockers without changing IKr conductance), and also it
naturally reflects many physiological and pharmacological factors(hERG channel density variability, hERG channel trafficking
block, etc.). It is possible that using different perturbations thesame system can show different robustness against EADs. For
example, the second best metric cqInward in terms of risk

Frontiers in Physiology | www.frontiersin.org

category separation does not correlate with the robustness to IKr
conductance reduction, but could potentially correlate with therobustness to other perturbations. We also note that even qNet
is not perfectly correlated with the robustness to IKr reduction.
The correlation between qNet and IKr reduction threshold waschecked only for 12 drugs at selected concentrations (0.5–25xCmax), and it is not known if the strong correlation holds true
beyond the drugs and concentrations tested. Even within the

12

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessmenton our current knowledge, they provide enough information
to correctly separate drugs into their TdP risk categories.
Additionally, calcium transient properties in the ORd model
differ from other models, such as the Grandi et al. model(Grandi et al., 2010); therefore, changes to the calcium transient
could improve prediction of TdP risk. In fact, Cummins et al.
incorporated the Grandi et al. model [along with the ORd andthe ten Tusscher et al. model (ten Tusscher and Panfilov, 2006)]
in their TdP risk classification and found diastolic intracellular
calcium and APD to be good markers of TdP risk (Cummins
et al., 2014). However, as mentioned earlier in this study
Cummins et al. define a binary TdP risk stratification that does
not follow the same categorization as defined by CiPA.
A number of different avenues for further improvementof the model and TdP risk prediction approach presented
here are currently being explored. We are examining theuse of thresholds for TdP risk level classification, as well as
incorporating both variability and uncertainty within the model
predictions. In conclusion, in this manuscript we present an
optimized version of the IKr-dyn ORd model presented in
Li et al. (2017) that is able to accurately separate the CiPA
training compounds into their respective risk categories andcorrelates well with the system’s robustness against EADs. An
independent validation of this approach is limited, but moreongoing work will see further refinement of this model andincreasing its suitability to be used routinely within the CiPA
paradigm.

concentrations tested, some drugs (for instance ranolazine) do
not have a consistent correlation across all the concentrations(Supplemental Figure 2). This suggests it may be beneficial touse the repolarization robustness (for instance IKr reduction
threshold) directly as a metric so that it has clear and direct
physiological meaning. However, this method is much morecomputationally intensive: for each drug at each concentration,
the computing time for the IKr reduction threshold is more than
10 times that for qNet, as multiple levels of perturbations areneeded to find the threshold. In addition, it is hard to define
the metric if different perturbations to the same system lead todifferent thresholds (robustness). Thus, using a highly correlated
surrogate metric qNet is a practical choice currently.

Limitations and Ongoing Work
While the model and metric combination presented here havebeen able to separate all the CiPA training compounds into theirrespective TdP risk categories, we have yet to test this approachon the CiPA validation compounds or any compounds that were
not used in the training of the model, which would provide an
independent validation of this framework. A key limitation of thisapproach that prevents an independent validation study is that
we have not provided thresholds for the qNet metric, which could
be used to place an unknown compound within a specific TdP
risk category. Instead we would only be able to group together
compounds which would be expected to pose similar TdP risk.
As suggested in previous studies the sodium potassium pump
(Lancaster and Sobie, 2016; Britton et al., 2017) and the sodium
calcium exchanger (Armoundas et al., 2003; Nagy et al., 2004)
play an important role in EAD generation. Simulations ofhypothetical drugs by Lancaster and Sobie (2016) show that boththe sodium potassium pump and sodium calcium exchanger
were ranked as having the greatest influence on TdP risk, above
IKs, IK1, and Ito (but excluding IKr, ICaL, and INa). Further
experiments and simulations are needed to assess how CiPA
drugs affect these currents and whether they should be directly
taken into account in our net current calculation to improve TdP
risk prediction.
Another key factor to consider is that while we havedemonstrated the success of our approach using gold standard
manual patch clamp data. At least in a pre-regulatory setting, theCiPA framework will likely rely on the use of high-throughput
ion channel screening data acquired from different platforms
routinely used within the pharmaceutical industry. We would
therefore need to further refine this model to fit to highthroughput system generated data and demonstrate that themodel and metric combination identified perform equally wellin this case. Furthermore, dynamic modeling of other channels(such as ICaL) may be needed as the project moves forward;
however, at this stage detailed kinetic drug block data forother channels is not available, nor are the protocols to extract
the necessary parameters. A priority of CiPA is to keep theframework simple and constrain the cost of data generation;
therefore, we use only IC50 data for other channels as, based

Frontiers in Physiology | www.frontiersin.org

AUTHOR CONTRIBUTIONS
JS, PT, MW, and WW contributed to the acquisition of data andprovided experimental data guidance. SD and ZL contributed
to designing the work and carried out simulations. TC and DS
contributed to revising the manuscript. SD, ZL, KC, and KB
contributed to the analysis and interpretation of the data andwriting of the manuscript.

FUNDING
This project was funded by the FDA Critical Path Initiative.

ACKNOWLEDGMENTS
We would like to thank Drs. Thomas O’Hara, Norman
Stockbridge, Richard Gray, Leonid Livshitz, Jules Hancox, Najah
Abi-Gerges, Jamie Vandenberg, Gary Mirams, Bernard Fermini,
Adam Hill, Meisam Hosseini and Jose Vicente for useful
discussions and input.

SUPPLEMENTARY MATERIALThe Supplementary Material for this article can be foundonline at: http://journal.frontiersin.org/article/10.3389/fphys.
2017.00616/full#supplementary-material

13

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

REFERENCES

in predicting torsade de pointes. Sci. Rep. 3:2100. doi: 10.1038/srep
02100
Kurata, Y., Matsuda, H., Hisatome, I., and Shibamoto, T. (2008). Regional
difference in dynamical property of sinoatrial node pacemaking: role of na+channel current. Biophys. J. 95, 951–977. doi: 10.1529/biophysj.107.112854
Kurata, Y., Tsumoto, K., Hayashi, K., Hisatome, I., Tanida, M., Kuda, Y.,
et al. (2017). Dynamical mechanisms of phase-2 early afterdepolarizations
in human ventricular myocytes: insights from bifurcation analyses of twomathematical models. Am. J. Physiol. Heart Circ. Physiol. 312, H106–H127.
doi: 10.1152/ajpheart.00115.2016
Kurita, T., Ohe, T., Marui, N., Aihara, N., Takaki, H., Kamakura, S., et al. (1992).
Bradycardia-induced abnormal QT prolongation in patients with complete
atrioventricular block with torsades de pointes. Am. J. Cardiol. 69, 628–633.
doi: 10.1016/0002-9149(92)90154-Q
Lancaster, C., and Sobie, E. A. (2016). Improved prediction of drug-induced
torsades de pointes through simulations of dynamics and machine learningalgorithms. Clin. Pharmacol. Ther. 100, 371–379. doi: 10.1002/cpt.367
Li, Z., Dutta, S., Sheng, J., Tran, P. N., Wu, W., Chang, K., et al. (2017).
Improving the in silico assessment of proarrhythmia risk by combining
herg (human ether-à-go-go-related gene) channel–drug binding kinetics
and multichannel pharmacology. Circ. Arrhythm. Electrophysiol. 10:e004628.
doi: 10.1161/CIRCEP.116.004628
Li, Z., Zhou, H., Lu, Y., and Colatsky, T. (2014). A critical role for
immune system response in mediating anti-influenza drug synergies assessedby mechanistic modeling. CPT Pharmacometrics Syst. Pharmacol. 3, 1–9.
doi: 10.1038/psp.2014.32
Mann, S., Imtiaz, M., Winbo, A., Rydberg, A., Perry, M., Couderc, J.P., et al. (2016). Convergence of models of human ventricular myocyte
electrophysiology after global optimization to recapitulate clinical long QT
phenotypes. J. Mol. Cell. Cardiol. 100, 25–34. doi: 10.1016/j.yjmcc.2016.09.011
Mirams, G., Cui, Y., Sher, A., Fink, M., Cooper, J., Heath, B., et al.
(2011). Simulation of multiple ion channel block provides improved early
prediction of compounds’ clinical torsadogenic risk. Cardiovasc. Res. 91, 53–61.
doi: 10.1093/cvr/cvr044
Nagy, Z., Virág, L., Tóth, A., Biliczki, P., Acsai, K., Bányász, T., et al. (2004).
Selective inhibition of sodium–calcium exchanger by SEA-0400 decreases earlyand delayed afterdepolarization in canine heart. Br. J. Pharmacol. 143, 827–831.
doi: 10.1038/sj.bjp.0706026
O’Hara, T., Virág, L., Varró, A., and Rudy, Y. (2011). Simulation of
the undiseased human cardiac ventricular action potential: model
formulation and experimental validation. PLoS Comput. Biol. 7:e1002061.
doi: 10.1371/journal.pcbi.1002061
Orth, P., Hesketh, C., Mak, C., Yang, Y., Lin, S., Beatch, G., et al. (2006).
RSD1235 blocks late INa and suppresses early afterdepolarizations andtorsades de pointes induced by class III agents. Cardiovasc. Res. 70, 486–496.
doi: 10.1016/j.cardiores.2006.01.026
Qu, Z., Xie, L. H., Olcese, R., Karagueuzian, H. S., Chen, P. S., Garfinkel, A.,et al. (2013). Early afterdepolarizations in cardiac myocytes: beyond reduced
repolarization reserve. Cardiovasc. Res. 99, 6–15. doi: 10.1093/cvr/cvt,104
R Core Team (2014). R: A Language and Environment for Statistical Computing.
Vienna: R Foundation for Statistical Computing. Available online at: http://
www.R-project.org/.
Redfern, W. S., Carlsson, L., Davis, A. S., Lynch, W. G., MacKenzie, I., Palethorpe,S., et al. (2003). Relationships between preclinical cardiac electrophysiology,
clinical QT interval prolongation and torsade de pointes for a broad rangeof drugs: evidence for a provisional safety margin in drug development.
Cardiovasc. Res. 58, 32–45. doi: 10.1016/S0008-6363(02)00846-5
Roden, D. M. (1998). Taking the “idio” out of “idiosyncratic”: predicting
torsades de pointes. Pacing Clin. Electrophysiol. 21, 1029–1034.
doi: 10.1111/j.1540-8159.1998.tb00148.x
Saad, M., Mahmoud, A., Elgendy, I., and Richard Conti, C. (2016). Ranolazine in
cardiac arrhythmia. Clin. Cardiol. 39, 170–178. doi: 10.1002/clc.22476
Sager, P., Gintant, G., Turner, R., Pettit, S., and Stockbridge, N. (2014).
Rechanneling the cardiac proarrhythmia safety paradigm: a meeting reportfrom the Cardiac Safety Research Consortium. Am. Heart J. 167, 292–300.
doi: 10.1016/j.ahj.2013.11.004
Soetaert, K., Petzoldt, T., and Setzer, R. W. (2010). Solving differential equations in
R: package deSolve. J. Stat. Softw. 33, 1–25. doi: 10.18637/jss.v033.i09

Antzelevitch, C., Belardinelli, L., Zygmunt, A., Burashnikov, A., Di Diego, J.,
Fish, J., et al. (2004). Electrophysiological effects of ranolazine, a novelantianginal agent with antiarrhythmic properties. Circulation 110, 904–910.
doi: 10.1161/01.CIR.0000139333.83620.5D
Armoundas, A., Hobai, I., Tomaselli, G., Winslow, R., and O’Rourke, B. (2003).
Role of sodium-calcium exchanger in modulating the action potential ofventricular myocytes from normal and failing hearts. Circ. Res. 93, 46–53.
doi: 10.1161/01.RES.0000080932.98903.D8
Britton, O. J., Bueno-Orovio, A., Virág, L., Varró, A., and Rodriguez, B. (2017).
The Electrogenic Na+ /K+ pump is a key determinant of repolarization
abnormality susceptibility in human ventricular cardiomyocytes: a populationbased simulation study. Front. Physiol. 8:278. doi: 10.3389/fphys.2017.00278
Chang, M., de Lange, E., Calmettes, G., Garfinkel, A., Qu, Z., and Weiss, J.
(2013). Pro- and antiarrhythmic effects of ATP-sensitive potassium current
activation on reentry during early afterdepolarization-mediated arrhythmias.
Heart Rhythm 10, 575–582. doi: 10.1016/j.hrthm.2012.12.017
Colatsky, T., Fermini, B., Gintant, G., Pierson, J., Sager, P., Sekino, Y.,
et al. (2016). The Comprehensive in Vitro Proarrhythmia Assay (CiPA)
initiative—Update on progress. J. Pharmacol. Toxicol. Methods 81, 15–20.
doi: 10.1016/j.vascn.2016.06.002
Crumb, W., Vicente, J., Johannesen, L., and Strauss, D. (2016). An evaluation of 30
clinical drugs against the comprehensive in vitro proarrhythmia assay (CiPA)
proposed ion channel panel. J. Pharmacol. Toxicol. Methods. 81, 251–262.
doi: 10.1016/j.vascn.2016.03.009
Cummins, M., Dalal, P., Bugana, M., Severi, S., and Sobie, E. (2014).
Comprehensive analyses of ventricular myocyte models identify targets
exhibiting favorable rate dependence. PLoS Comput. Biol. 10:e1003543.
doi: 10.1371/journal.pcbi.1003543
Di Veroli, G., Davies, M., Zhang, H., Abi-Gerges, N., and Boyett, M. (2013).
High-throughput screening of drug-binding dynamics to HERG improves earlydrug safety assessment. Am. J. Physiol. Heart Circ. Physiol. 304, H104–H117.
doi: 10.1152/ajpheart.00511.2012
Di Veroli, G., Davies, M., Zhang, H., Abi-Gerges, N., and Boyett, M. (2014). hERG
inhibitors with similar potency but different binding kinetics do not pose thesame proarrhythmic risk: implications for drug safety assessment. J. Cardiovasc.
Electrophysiol. 25, 197–207. doi: 10.1111/jce.12289
Dutta, S., Strauss, D., Colatsky, T., and Li, Z. (2016). “Optimization of an in silico
cardiac cell model for proarrhythmia risk assessment,” in 2016 Computing in
Cardiology Conference (CinC) (Vancouver, BC), 869–872.
Fermini, B., Hancox, J., Abi-Gerges, N., Bridgland-Taylor, M., Chaudhary, K.,
Colatsky, T., et al. (2016). A new perspective in the field of cardiac safety testing
through the comprehensive in vitro proarrhythmia assay paradigm. J. Biomol.
Screen. 21, 1–11. doi: 10.1177/1087057115594589
Gintant, G. (2008). Preclinical torsades-de-pointes screens: advantages andlimitations of surrogate and direct approaches in evaluating proarrhythmic
risk. Pharmacol. Ther. 119, 199–209. doi: 10.1016/j.pharmthera.2008.04.010
Grandi, E., Pasqualini, F., and Bers, D. (2010). A novel computational model of thehuman ventricular action potential and Ca transient. J. Mol. Cell. Cardiol. 48,
112–121. doi: 10.1016/j.yjmcc.2009.09.019
Gray, R., and Huelsing, D. (2008). Excito-oscillatory dynamics as amechanism of ventricular fibrillation. Heart Rhythm 5, 575–584.
doi: 10.1016/j.hrthm.2008.01.011
Guo, D., Zhao, X., Wu, Y., Liu, T., Kowey, P., and Yan, G.-X. (2007). Ltype calcium current reactivation contributes to arrhythmogenesis associatedwith action potential triangulation. J. Cardiovasc. Electrophysiol. 18, 196–203.
doi: 10.1111/j.1540-8167.2006.00698.x
Johannesen, L., Vicente, J., Mason, J. W., Erato, C., Sanabria, C., Waite-Labott, K.,
et al. (2016). Late sodium current block for drug-induced long QT syndrome:
results from a prospective clinical trial. Clin. Pharmacol. Ther. 99, 214–223.
doi: 10.1002/cpt.205
Kallergis, E., Goudis, C., Simantirakis, E., Kochiadakis, G., and Vardas, P.
(2012). Mechanisms, risk factors, and management of acquired long QT
syndrome: a comprehensive review. ScientificWorldJournal. 2012:212178.
doi: 10.1100/2012/212178
Kramer, J., Obejero-Paz, C., Myatt, G., Kuryshev, Y., Bruening-Wright, A.,
Verducci, J., et al. (2013). MICE models: superior to the HERG model

Frontiers in Physiology | www.frontiersin.org

14

August 2017 | Volume 8 | Article 616

Dutta et al.

In silico Proarrhythmia Risk Assessment

Xie, Y., Sato, D., Garfinkel, A., Qu, Z., and Weiss, J. (2010). So little source, so muchsink: requirements for afterdepolarizations to propagate in tissue. Biophys. J. 99,
1408–1415. doi: 10.1016/j.bpj.2010.06.042
Yu, H. (2002). Rmpi: parallel statistical computing in R. R News 2, 10–14. Available
online at: https://cran.r-project.org/web/packages/Rmpi/citation.html

ten Tusscher, K. H., and Panfilov, A. V. (2006). Cell model for efficientsimulation of wave propagation in human ventricular tissue undernormal and pathological conditions. Phys. Med. Biol. 51, 6141–6156.
doi: 10.1088/0031-9155/51/23/014
Tierney, L., Rossini, A. J., Li, N., and Sevcikova, H. (2015). Snow: Simple Network ofWorkstations. Available online at: https://CRAN.R-project.org/package=snow;
https://cran.r-project.org/web/packages/snow/snow.pdf
Volders, P., Vos, M., Szabo, B., Sipido, K., Marieke de, G., Gorgels, A., et al.
(2000). Progress in the understanding of cardiac early afterdepolarizations
and torsades de pointes: time to revise current concepts. Cardiovasc. Res. 46,
376–392. doi: 10.1016/S0008-6363(00)00022-5
Vos, M., Verduyn, C., Gorgels, A., Lipcsei, G., and Wellens, H. (1995).
Reproducible induction of early afterdepolarizations and torsade depointes arrhythmias by D-sotalol and pacing in dogs with chronic
atrioventricular block. Circulation 91, 864–872. doi: 10.1161/01.CIR.
91.3.864
Weiss, J., Garfinkel, A., Karagueuzian, H., Chen, P. S., and Qu, Z. (2010). Early
afterdepolarizations and cardiac arrhythmias. Heart Rhythm 7, 1891–1899.
doi: 10.1016/j.hrthm.2010.09.017
Wu, J., Wu, J., and Zipes, D. (2002). Early Afterdepolarizations, U Waves, andTorsades de Pointes. Circulation 105, 675–676. doi: 10.1161/hc0602.104458

Frontiers in Physiology | www.frontiersin.org

Conflict of Interest Statement: The authors declare that the research was
conducted in the absence of any commercial or financial relationships that could
be construed as a potential conflict of interest.
The reviewer CC and handling Editor declared their shared affiliation, andthe handling Editor states that the process met the standards of a fair and objective
review.
Copyright © 2017 Dutta, Chang, Beattie, Sheng, Tran, Wu, Wu, Strauss, Colatsky
and Li. This is an open-access article distributed under the terms of the Creative
Commons Attribution License (CC BY). The use, distribution or reproduction inother forums is permitted, provided the original author(s) or licensor are credited
and that the original publication in this journal is cited, in accordance with accepted
academic practice. No use, distribution or reproduction is permitted which does notcomply with these terms.

15

August 2017 | Volume 8 | Article 616