Optimization of an In silico Cardiac Cell Model for Proarrhythmia Risk Assessment

Drug-induced Torsade-de-Pointes (TdP) has been responsible for the withdrawal of many 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 and in 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 mathematical model 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 carried by 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 published human cardiomyocyte experimental data under control and drug block conditions. Using this optimized model and manual patch clamp data, we developed an updated version of the metric that quantifies the net electronic charge carried by major inward and outward 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 that is more consistent with experimental data, an improved metric that can classify drugs at concentrations both near and higher than clinical exposure, and a physiological framework 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.

The dynamic IKr model within the optimized IKr-dyn ORd model was the same as described in our previous paper and the details can be found there (Li et al. 2017). Briefly, it contains a physiological component representing channel gating and a pharmacodynamic component representing drug binding. The physiological component is based on our recently published temperature-dependent hERG model (Li et al., 2016) and further optimized to recapitulate native IKr behavior. It has two closed states (C), two inactivated-closed states ( ⁄ , where Kmax is the maximum drug effect at saturating concentrations, D is the concentration/dose of the drug, n is the Hill coefficient, and EC50 is the concentration when half maximum effect is achieved. The "trapping rate" (Kt) of channel closing with drug bound (IO*→C* and O*→C*) was manually fixed at 3.5 • 10 −5 ms −1 . The rate of channel opening when drug is bound (C*→O* and C*→IO*) is equal to Kt • X(V), where X is the steady-state gate activation function for IKr in ORd model (O'Hara et al., 2011) and V is the transmembrane potential in millivolts (mV). X(V) is defined as 1 �1 + exp (−(V − Vhalf ̵ trap) 6.789 ⁄ ))� ⁄ , where Vhalf-trap is the membrane voltage at which half of drug-bound channels are open. The IO* to IO transition rate is determined by microscopic reversibility. Some important considerations while developing this model are listed below.

IC50 vs. EC50
A Hill equation is often used to describe the concentration-dependent steady state fraction of block and the IC50 is the drug concentration where the steady-state inhibition reaches 50% of maximally possible block level. On the other hand, in the field of mechanism-based pharmacodynamics modeling a Hill equation-like Emax format is often used to describe concentration-dependent reaction rate (Felmlee et al., 2012) where EC50 is the drug concentration when the reaction rate reaches 50% of the maximum speed. We applied the Emax format to concentration-dependent binding rates governing

Unusually high EC50 values
When the EC50 is small (close to concentrations tested experimentally), the Emax equation used is the bona fide Emax model (i.e. the relationship between drug concentration and reaction rates O→O* or IO→IO* is sigmoidal with a saturating effect). When the EC50 is much larger than the highest concentration tested then the Emax model becomes a liner model because Kmax/ (D n + EC50 n ) ≈ Kmax EC50 n ⁄ . In this case the ratio Kmax EC50 n ⁄ describes the linear relationship between drug concentration and reaction rates for O→O* or IO→IO*. Under this "linear approximation" circumstances the two parameters Kmax and EC50 n are unidentifiable. However, it's their ratio that really matters, and the ratio is actually identifiable via a global identifiability test (Table S2 in (Li et al., 2017)): consistent ratio values were obtained from repeated fitting runs with different starting points (random seeds).

The choice of Emax model
Currently, the intended purpose of this model within a regulatory pipeline (Comprehensive in vitro Proarrhythmia Assay), constrains all new compounds to be evaluated/simulated by the same model. Therefore, the drug binding component of this model needs to be flexible enough to represent different mode of binding, for example linear vs Emax form of binding rate. Because most current drug-hERG dynamic models assumes only linear relationship between the drug concentration and binding rate, we first tried a linear form of binding rate where the reaction rates O→O* or IO→IO* are of the form Kf • D (Kf is diffusion rate while D is drug concentration). We found that a much worse fit was achieved for some compounds when assuming the linear drug block model (data not shown). We think it is partly because some drugs (like cisapride) have a saturating effect in its concentration-binding rate relationship, and only a model like Emax could capture that. Therefore, the current Emax formula was chosen to accommodate both a bona fide Emax model (with suturing effect) and (through linear approximation) a linear model (no saturating effect).

The choice of Kt value
When we developed the model a key "reality check" is to reproduce different drugs' ability to cause reverse use dependency (RUD) at APD prolongation (Fig 5 of our previous paper (Li et al., 2017)). We found that when allowing the Kt parameter to vary (through fitting) the RUD pattern is incorrect: dofetilide may have no RUD while verapamil may have it, which is contradictory to many published data. The value Kt was therefore fixed and was selected to best replicate these reported RUD patterns.

Number of bound states
An important consideration for CiPA is ease of implementation; therefore, we decided to find the simplest drug binding model to explain the data. The simplest drug binding format in the literature, which is widely used, is assuming two drug-bound states: open-bound and inactivated-bound, corresponding to O* and IO*. However, it didn't work well, probably because it didn't capture the "trapping" phenomenon (Li et al., 2017). After that we found the current drug bound model (three drug-bound states O*, IO*, and C*) was the simplest addition to the base IKr model which was able to provide the best fit to the fractional block data. While we could have explored the use of a more complex model including additional closed and inactivated bound states this would have increased the complexity of the model structure so we selected the simpler drug bound component with just one open, inactivated and closed state.

Supplemental Tables
The following Supplemental Tables are a reproduction of Tables 1, 2 and 3 from (Li et al., 2017). They are included here to give more context about the current work.  (Li et al., 2017)). For other drugs like cisapride and diltiazem, the estimated EC50 n s are not significantly higher than the tested concentrations (EC50 n < 10-fold of estimated IC50), suggesting there is indeed a sigmoidal doseresponse relationship. The free Cmax values are either directly taken from the referenced paper, or in the case of ondansetron, calculated using reported total Cmax and protein binding data from FDA drug labels.  (Redfern et al., 2003) Supplemental Table 3. IC50s used in AP simulation.

Supplemental
IC50s and Hill coefficients (h's) for IKr are estimated using Milnes protocol data in this study (Milnes et al., 2010;Li et al., 2017). Other channels' IC50 and h values are calculated using the blocking data from Crumb et al. (Crumb et al., 2016). Zeros in the table indicate no detectable blocking for that channel. Note that many drugs' IC50s are very high because only a narrow concentration range around Cmax was tested in Crumb study (Crumb et al., 2016). In those cases, the estimated IC50s may not be accurate in predicting block at high concentrations. Therefore, all simulations in this study were limited to ≤25x Cmax.  Table 4 (main text) for correlation coefficients.
Supplemental Figure 3: Relationship between the IKr reduction threshold and APD90 for 0.5x (dark blue) to 25x (light blue) Cmax for each drug. Results where EADs occur without added IKr reduction (quinidine ≥ 2.3x Cmax) and results where the maximum IKr reduction did not trigger an EAD (diltiazem all Cmax; verapamil ≥ 1.7x Cmax; mexiletine ≥ 3.8x Cmax) were excluded. See Table 4 (main text) for correlation coefficients ( Table 4 shows that mexiletine and verapamil have positive correlation coefficients as opposed to negative correlation coefficients for all other drugs).
Supplemental Figure 4: Relationship between the IKr reduction threshold and the normalized charge passed by INaL and ICaL (cqInward (Li et al., 2017)) for 0.5x (dark blue) to 25x (light blue) Cmax for each drug. Results where EADs occur without added IKr reduction (quinidine ≥ 2.3x Cmax) and results where the maximum IKr reduction did not trigger an EAD (diltiazem all Cmax; verapamil ≥ 1.7x Cmax; mexiletine ≥ 3.8x Cmax) were excluded. See Table 4 (main text) for correlation coefficients.
Supplemental Figure 5: Transmembrane voltage (Trans. voltage) through time for ranolazine (black solid line) and cisapride (dashed gray line) at 1x Cmax without (left panel) and with 91.6% IKr reduction (right panel) for a CL of 2000 ms. Corresponding APD90 (ms) and qNet (µC/µF) are reported in black for ranolazine and in gray for cisapride. 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).