Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Early afterdepolarisation tendency as a simulated pro-arrhythmic risk indicator

Beth McMillan *a, David J. Gavaghan a and Gary R. Mirams b
aComputational Biology, Dept. of Computer Science, University of Oxford, Oxford, OX1 3QD, UK. E-mail: beth.mcmillan@cs.ox.ac.uk; Fax: +44 (0)1865 283532; Tel: +44 (0)1865 273838
bCentre for Mathematical Biology, School of Mathematical Sciences, University of Nottingham, Nottingham, NG7 2RD, UK

Received 18th May 2017 , Accepted 12th September 2017

First published on 14th September 2017


Drug-induced Torsades de Pointes (TdP) arrhythmia is of major interest in predictive toxicology. Drugs which cause TdP block the hERG cardiac potassium channel. However, not all drugs that block hERG cause TdP. As such, further understanding of the mechanistic route to TdP is needed. Early afterdepolarisations (EADs) are a cell-level phenomenon in which the membrane of a cardiac cell depolarises a second time before repolarisation, and EADs are seen in hearts during TdP. Therefore, we propose a method of predicting TdP using induced EADs combined with multiple ion channel block in simulations using biophysically-based mathematical models of human ventricular cell electrophysiology. EADs were induced in cardiac action potential models using interventions based on diseases that are known to cause EADs, including: increasing the conduction of the L-type calcium channel, decreasing the conduction of the hERG channel, and shifting the inactivation curve of the fast sodium channel. The threshold of intervention that was required to cause an EAD was used to classify drugs into clinical risk categories. The metric that used L-type calcium induced EADs was the most accurate of the EAD metrics at classifying drugs into the correct risk categories, and increased in accuracy when combined with action potential duration measurements. The EAD metrics were all more accurate than hERG block alone, but not as predictive as simpler measures such as simulated action potential duration. This may be because different routes to EADs represent risk well for different patient subgroups, something that is difficult to assess at present.


Introduction

Torsades de Pointes (TdP) is a particular type of polymorphic ventricular tachycardia, characterised by an unusual electrocardiogram, in which the QRS complex appears to be twisted around the baseline. TdP usually spontaneously resolves, sometimes causing syncope (sudden fainting due to a drop in blood pressure), but it can also cause cardiac arrest or sudden death.1,2

Anti-arrhythmic drugs, such as quinidine,3 are commonly linked to TdP. Terfenadine, a non-sedating anti-histamine, was one of the first non-cardiac drugs to be linked to increased TdP risk.4 Terfenadine was withdrawn from the market in 1997 after being linked to 41 cases of TdP, one of which was lethal.5 Cisapride is a drug that was used for treating gastroesophageal reflux disease.6 After causing 97 cases of TdP, of which six were fatal, cisapride was withdrawn from the market.7

Quinidine, terfenadine, and cisapride were all found to strongly block IKr, the rapid delayed rectifying potassium current in the heart, which is carried by the channel whose primary subunit is a product of the human ether-a-go-go related gene (hERG).4,6,8,9

Since these discoveries, testing for hERG block has become a mandatory requirement for new pharmaceuticals.10 hERG block as a measure of TdP risk is very sensitive (gives few false negatives) and has prevented torsadogenic drugs from entering the market. Certain marketed drugs, such as verapamil and ranolazine, block hERG but are not linked with TdP.11,12 There are therefore concerns that hERG block lacks specificity (gives false positives for TdP risk), preventing the development of potentially useful drugs.13 As such, elucidation of other factors that mediate TdP risk is needed.

Early afterdepolarisations (EADs), phenomena in which the membrane depolarises a second time during the action potential, are heavily implicated in the onset of TdP.15 Hearts suffering from TdP show EADs alongside transmural dispersion of repolarisation.16–19 EADs are seen in monophasic action potential recordings from dog hearts during TdP, as shown in Fig. 1.14,19


image file: c7tx00141j-f1.tif
Fig. 1 Induction of EADs and TdP in an anaesthetised dog by treatment with ibutilide.14 Reproduced with permission from P. G. A. Volders, et al., Progress in the understanding of cardiac early after depolarizations and torsades de pointes: time to revise current concepts, Cardiovasc. Res., 2000, 46(3), 376–392. Published by Oxford University Press on behalf of the European Society of Cardiology (ESC).

Our investigation is based on several mechanisms that are known to promote EADs. One mechanism is based on a form of Long QT syndrome (LQT8) that is caused by an increase in the L-type calcium current, caused by a gain of function mutation in the CACNA1C gene, that increases the density of L-type calcium channels at the cell surface up to three-fold.20,21 Patients with LQT8 frequently exhibit TdP.22 L-type calcium current agonists (agents which increase the conductance of the L-type calcium channel) have also been shown to cause: TdP in live mice;23 EADs and ventricular tachycardia in intact mouse hearts;24 EADs in sheep and dog Purkinje fibres;25 and EADs in ferret ventricular myocytes.26 L-type calcium current blockers such as verapamil are used as anti-arrhythmic agents, and have been shown to suppress EADs and TdP in rabbit models of LQT3 and chronic heart failure.11,27,28

Another EAD mechanism we investigate in this study is based on Brugada syndrome, which is a condition linked to ventricular fibrillation and elevation in the ST segment of the electrocardiogram.29 Estimates of the proportion of sudden deaths caused by Brugada vary from 4–12%, and it is estimated to be present in 0.05% of the world population. In 15–20% of Brugada cases, the cause is a mutation in the SCN5A gene, which codes for the alpha subunit of the fast sodium channel.29 A missense mutant version of SCN5A (T1620M), was shown to alter the inactivation curve of the fast sodium current, shifting it in the positive direction by 10 mV,30 and another missense mutant, L812Q, shifted the inactivation curve by 20 mV.31 Adding this inactivation curve shift to cardiac action potential models can cause EADs (Noble et al., personal communication).

In addition to being of interest in drug-induced TdP, the IKr potassium current is also involved in LQT2, and is also linked to EADs.9 Loss-of-function mutations in hERG cause a reduction in the IKr current of up to 97%,32 which causes QT prolongation, and usually increases risk of TdP. LQT mutations do not always present with significant clinical QT prolongation, but interaction with drugs that affect cardiac ion channels could cause an increased risk of TdP.33

Two other Long QT conditions that are linked to TdP are LQT1 and LQT3.34 LQT1 causes an increase in the slow delayed potassium current IKs, and LQT3 causes an increase in the persistent sodium current IpNa. We have not used these as EAD-provoking interventions in this work because we were unable to provoke EADs in the O'Hara cell model by blocking IKs or increasing IpNa.

The prediction of arrhythmias by in silico modelling of action potentials in response to ion channel block offers a new way to test novel compounds at the pre-clinical stage. A previous study by our group created an improved measure of a compound's propensity for causing TdP arrhythmias, using simulated action potential duration as a metric.35 The approach takes into account the contributions of multiple ion channels to the shape and length of the action potential, and classifies drugs into discrete risk categories, based on their effect on action potential duration. This method was more accurate than the commonly-used ‘hERG safety factor’, that is the ratio of hERG IC50 to effective free therapeutic plasma concentration (EFTPC), or log10(hERG IC50/EFTPCmax).36 Simulation studies have been extended to predict results of rabbit wedge studies and the Thorough QT study.37,38

A recent study used principal component analysis to assemble a large number of biomarkers from different models, the results suggested that a two-dimensional binary classification based on both the simulated diastolic calcium concentration and the APD50 was effective at separating drugs into positive or negative for torsadogenicity.39 We aim to extend these approaches to allow prediction of the TdP risk classes of drugs that increase action potential duration but are safe (particularly late/persistent sodium blockers), and to account for the interaction of drug block with underlying conditions such as ion channel mutations. The appearance of EADs at increased drug concentrations has been studied computationally as a risk indicator for TdP.40 Our study complements this work by exploring EADs as a risk indicator in the context of disease states at clinically-relevant concentrations.

Ion channel conductance modification as a cause of EADs has also been investigated recently in the context of atrial fibrillation: a global sensitivity analysis of atrial cell models was used to examine which ion channel changes lead to EADs and then logistic regression was used to estimate the probability of EADs as functions of conductances.41

The Comprehensive in vitro Pro-arrhythmia Assay (CiPA) is a proposal to use multi-ion channel screening in combination with human stem cell-derived cardiomyocytes and computational cardiac modelling to create new metrics for the prediction of drug-induced TdP, moving away from using QT interval prolongation as a surrogate marker and towards in vitro and in silico methods.42 This study investigates whether a marker linked mechanistically to EAD formation may be helpful in determining torsadogenic risk.

We applied ion channel block to cardiac cell models to simulate the effects of drugs of known torsadogenic risk, then combined these drug effects with simulated disease states. The level of disease state necessary to provoke an EAD for each drug was used as a marker for TdP risk. These markers were evaluated against clinical risk categories, separately, combined, and in combination with pre-existing TdP risk metrics. By combining these markers we hope to create a measure of risk for a population that contains people with each EAD-inducing condition. In this way, we bring together three key factors that influence arrhythmogenic risk: the ion channel blocking properties of the compounds, the concentrations of compound found in humans, and underlying risk factors from the patient.

Methods

We selected 41 drugs of known torsadogenic risk, and simulated their ion channel blocking effects in the O'Hara 2011 human ventricular cell model.43 Using a range of interventions, we determined the threshold of intervention at which an EAD could be provoked in the cell model, i.e. the lowest level of intervention that was necessary for an EAD to be produced. The differences in ‘threshold for EAD’ between different drugs were used to classify drugs by arrhythmic risk. These steps are detailed below.

Drug inclusion criteria

To select drugs to use as a training set, we used three criteria based on the amount of data available on both the pro-arrhythmic risk of the compound and its effect on ion currents in cardiac cells.

1. As a starting point, we included drugs that we previously studied in Mirams et al.35

2. Additional drugs were then included in our study if they had been included in five or more of the papers analysing TdP risk discussed in a recent summary paper,44 and if over 70% of studies agreed on high or low TdP risk.

3. Drugs that are on the CiPA list42 were automatically included if ion current block data were available for three or more channels of interest (even if they had been in fewer than five studies or had poor agreement in risk category between studies).

If there was disagreement in risk class using the above sources, the default category was the one used in Mirams et al.,35 and if not listed there then the one used in Redfern et al.36

To be included in our dataset, the drugs were also required to have IC50 values available in the literature for three or more of the ionic currents of interest, found by manual patch clamp. The ionic currents of interest were: the fast and late/persistent sodium currents, the L-type calcium current, the rapid and slow delayed rectifier potassium currents, the transient outward current, and the inward rectifier potassium current. The pIC50 values we used are given in Table 1, and the references for these can be found in the ESI. Where there were multiple IC50 values available, preference was given to those from human cells at body temperature (37 °C) and to the most recent study.

Table 1 pIC50 values for each compound in the dataset (shown as log[thin space (1/6-em)]M) for the fast sodium (INa), L-type calcium (ICaL), rapid delayed rectifier potassium (hERG or IKr), slow delayed rectifier potassium (IKs), persistent sodium (IpNa), transient outward (Ito) and inward rectifier potassium (IK1) currents, and effective free therapeutic plasma concentration (EFTPC) (nM). For references, please see the ESI. ‘n/a’ indicates that the channel has been screened and no effect was measured: either the IC50 value was above the maximum concentration being tested, or there is no drug-induced block of this channel
Compound pIC50 (log[thin space (1/6-em)]M) EFTPCmax (nM)
I Na I CaL hERG I Ks I pNa I to I K1
Ajmaline 5.0862 4.1487 5.9830 5.5751 17[thin space (1/6-em)]103
Amiodarone 4.3936 6.5686 7.5229 5.7595 5.1739 5.4250 0.5
Amitriptyline 4.6990 4.9355 5.4841 5.5627 5.3533 5.0000 41
Azimilide 4.7212 4.7496 7.0000 5.8539 70
Bepridil 5.4318 6.6757 6.2218 5.0000 5.7414 5.1805 33
Chlorpromazine 5.5528 5.5229 5.8297 5.3410 4.7959 5.2147 38
Cibenzoline 5.1079 4.5229 4.6459 4.9101 4.3318 976
Cisapride 4.8327 4.7959 8.1871 5.4698 5.0301 4.9
Clozapine 5.6021 1603
Desipramine 5.8182 5.7673 5.8570 108
Diltiazem 4.8477 7.2676 5.1314 4.5346 122
Disopyramide 5.7447 4.0910 4.6799 742
DL-Sotalol 3.5560 14[thin space (1/6-em)]733
Dofetilide 3.5229 4.5735 8.3010 2
Flecainide 5.1871 4.5918 5.7959 5.4685 5.1079 753
Fluvoxamine 4.4045 5.3098 5.4202 264
Halofantrine 7.6655 0.5
Imipramine 5.4437 5.0915 5.4685 4.3010 106
Loratadine 5.0809 0.45
Methadone 4.9508 4.5735 5.3188 507
Mexiletine 4.3665 4.0000 4.3010 4.7545 4129
Mibefradil 6.0088 6.8069 5.7447 5.4403 12
Nifedipine 4.4318 7.2218 3.5607 3.4437 4.8633 7.7
Nitrendipine 4.6655 7.6021 5.0000 5.1135 3.02
Ondansetron 4.0531 4.6468 6.0915 4.7171 899.88
Pentamidine 3.6882 6.7696 10
Phenytoin 4.3098 3.9872 4.0000 7000
Pimozide 7.2676 6.6198 7.8239 0.43
Prenylamine 5.5986 5.9066 7.1871 17
Propafenone 5.9245 5.7447 6.3565 5.3940 5.3188 5.1487 241
Propranolol 5.6778 4.7447 5.5485 26
Quetiapine 4.7721 4.9830 5.2392 33
Quinidine 4.7799 4.8069 6.5229 5.3099 3.0580 3.6990 3237
Ranolazine 3.5317 3.5287 4.9393 2.7212 5.1871 3200
Risperidone 3.9914 4.1367 6.2218 1.81
Sertindole 5.6383 5.0506 7.8539 6.0555 5.3979 1.59
Tedisamil 4.6990 5.6021 5.3565 85
Terfenadine 6.0128 6.4260 8.0506 5.6990 9
Terodiline 4.8182 8.3979 4.5229 5.1549 12
Thioridazine 5.7375 5.8861 7.4815 4.8539 979
Verapamil 4.3820 7.0000 6.8447 2.5796 81


Stratification into 4 risk categories was taken from the Mirams et al. (2011) paper where available, and the Redfern et al. (2003) paper otherwise. We pooled risk categories 1 and 2 together, as they represent the same risk level for different drug classes, leaving classes 2–5.35,36 The only compounds in the dataset not covered by this metric were: (i) ranolazine, which was assigned to category 4, i.e. the drug has isolated reports of TdP in humans; and (ii) cibenzoline, which was assigned to category 5 based on Lawrence et al.,45i.e. there are no reports of TdP in humans with this drug.46,47

Action potential and drug block models

We used two recent human ventricular myocyte electrophysiology models based on human datasets, the O'Hara (2011)43 endocardial and Grandi (2010)49 models. Drug block of ion channels was modelled as a reduction in channel conductance as a function of the concentration of the compound, [D], and the IC50 value.50 The change in maximum conductance for a channel j was described by a Hill equation:
 
image file: c7tx00141j-t1.tif(1)
where gj is the maximum conductance of the drug-blocked channel, and j is the conductance of the channel when there is no compound present. The Hill coefficient here is set to 1, as the variability in experimentally-inferred Hill coefficients from patch clamp can be so high that using 1 as the Hill coefficient may reduce error.51 Note that this formula was applied for all channel/drug combinations listed in Table 1, apart from late/persistent sodium in the Grandi model—as this model does not have a distinct late/persistent sodium current.

Drug concentrations were set to the maximum effective free therapeutic plasma concentrations (EFTPC) for each individual drug, to provide a realistic estimate of ion channel block in vivo. For a list of EFTPCs, see Table 1.

Detecting afterdepolarisations

The appearance of early afterdepolarisations (EADs) in a simulation was determined by the slope of the voltage trace between adjacent time points. Whenever the slope was greater than +1 mV ms−1, a depolarisation was reported. To remove depolarisations caused by the stimulus, depolarisations that occurred between 50 ms before and 100 ms after each stimulus were disregarded. The algorithm detects single EADs, multiple EADs, and EADs without repolarisation, as shown in Fig. 2. For the full algorithm, see the DetectAfterDepolarisations class in the code repository (see “Numerical methods and simulation procedure”).
image file: c7tx00141j-f2.tif
Fig. 2 Types of EADs that are detected using the EAD detection algorithm. Left: A single EAD with repolarisation. Centre: Multiple EADs with repolarisation. Right: Multiple EADs without repolarisation.

Provoking afterdepolarisations

The failure of hERG block alone to predict torsadogenic risk suggests that drug-induced TdP is mediated by more than one ionic mechanism. We hypothesise that the interaction of certain disease states with torsadogenic drugs could lead to greater susceptibility to TdP. We propose to look at the single-cell phenomenon of EADs, as studies suggest TdP and EADs are intrinsically linked.14,16–19 Spatial differences in ion channel expression throughout cardiac tissue would alter the susceptibility of each cell to EADs, meaning that some regions would have EADs while others did not, increasing the probability of re-entrant waves, but we constrain this study to looking at the ‘trigger’ event of EAD initiation.

To simulate disease states after drug block was applied, EADs were provoked using three interventions mimicking disease states discussed in the introduction, based on a manuscript from Noble et al.:52

• conductance of the L-type calcium current was increased by a scaling factor to simulate a gain-of-function mutation which causes LQT8.21

• the midpoint of the inactivation curve of the fast sodium current was shifted using an additive constant to simulate Brugada syndrome.30

• conductance of the rapid delayed rectifier current IKr was decreased by a scaling factor to simulate the effect of LQT2 mutations.9

The O'Hara et al.43 model was used for these EAD simulations, as it includes all of the currents of interest, and is being investigated for CiPA-related proarrythmic risk prediction.42

We used a slow pacing interval of 3 seconds because bradycardia is also known to facilitate onset of EADs and TdP.17 We performed the same simulations at faster pacing rates, which did not significantly change the classification results.

Previously suggested measures

In order to compare our EAD tendency measures with previously suggested risk indicators, we used our dataset of drug actions to calculate some risk metrics that have previously been proposed.

We calculated the hERG IC50/EFTPCmax safety margin proposed by Redfern et al.;36 the predicted APD90 calculated using the Grandi et al.49 human ventricular cell model, as in Mirams et al.;35 the predicted APD50 combined with diastolic calcium concentration calculated using both the Grandi et al.49 model, and the O'Hara et al.43 model as in Lancaster and Sobie;39 and the predicted sum of normalised total persistent sodium current and L-type calcium current over the course of the action potential at increasing drug concentrations (known as ‘cqInward’) using the O'Hara et al.43 model, as in Li et al.,48 but using the baseline model rather than their dynamic hERG block model (as we did not have hERG kinetic data for all compounds). The Grandi model was chosen for APD90 rather than the O'Hara model as it was used in the Mirams et al. (2011) metric.

Numerical methods and simulation procedure

Simulations were run using the Chaste C++ framework53,54 with the ApPredict bolt-on project55 and custom written code. Cell models and initial conditions were imported from CellML56 files using PyCML.57

The adaptive timestep solver CVODE58 was used to solve the model's differential equations, with a relative tolerance of 10−5 and absolute tolerance of 10−7. The output timestep was 0.1 ms and a stimulus current was applied every 3 s for 3 ms with a magnitude of −25.5 μA μF−1.

The procedure we used is outlined in Algorithm 1. First, drug block was applied to the model by modifying the conductance parameters for the appropriate channels, using eqn (1) to calculate updated conductances. The model was then run to steady state using the pacing protocol above. Steady state was reached when the norm of the change in the model's state variables was less than 10−6 between paces. This generally took 1000–10[thin space (1/6-em)]000 paces.

EAD-provoking interventions were applied by modifying either the ICaL conductance, the IKr conductance, or the INa shift parameter. After setting the intervention parameter, the simulation was run for 12 s, and the presence or absence of EADs was detected. Interval bisection was used to find the threshold of intervention that was necessary to cause any EADs. We started with interval ranges of: [0,20] mV for the INa shift intervention; [1,80] for the ICaL conductance scaling; and [0,1] for the IKr conductance scaling. Interval bisection was set to terminate when consecutive values were less than 10−4 units apart. EAD simulations for a single compound can be run in less than 10 minutes, and many compounds can be run in parallel. All the code is available to download from github at https://github.com/teraspawn/EadPredict.

Algorithm 1 Procedure to find the EAD threshold.
1: Set ion channel conductances in cell model to new values based on drug block.
2: Run cell model to steady state and then save this state to reset to later.
3: Set upper and lower limits α, β for intervention value.
4: Repeat
5: Reset cell model to earlier state.
6: Set intervention to (α + β) ÷ 2 (e.g. multiply ICaL conductance by this factor).
7: Run model for 12 s.
8: Check for early afterdepolarisation (upwards trajectory after initial depolarisation).
9: Adjust α or β using interval bisection.
10: Until |αβ| < 10−4

Linkage analysis

To visualize whether drugs with similar risk require similar EAD thresholds, linkage analysis was used to create dendrograms of drug similarity based on each of the metrics. The Euclidean distance between the metric values (e.g. APD90 for each of the drugs) was used to construct a dendrogram, grouping drugs which had similar values. To combine metrics, the results were shifted such that the control value was zero, and then scaled to be within [−1,1]. Combinations of metrics such as all the EAD metrics, the EAD metrics with APD90, and all the metrics together, were then used to create classification trees. The optimal leaf ordering was calculated using the R package “dendextend”, to best sort the leaves in descending order of risk category (without changing the branching structure).59

This method allows for the grouping of drugs by similarity rather than by rigid categories, allowing for new compounds to be visually ranked by closeness to torsadogenic and non-torsadogenic drugs. The output of the optimal leaf ordering algorithm was evaluated by the sum of the square difference between the risk category of each drug in the ranking and an optimal ordering (2,2,…,2,3,3,…,3,4,4,…etc.):

 
image file: c7tx00141j-t2.tif(2)
where o is the optimal ordering of risk category, a is the actual risk category, and N is the number of drugs in the dataset.

Pro-arrhythmic risk classification

As in Mirams et al.,35 risk categories were based on the Redfern et al.36 classes (with category 1 and category 2 pooled, as they represent equivalent levels of risk for different drug classes):

• Category 2: Either Class Ia and III antiarrhythmics, or drugs that have been withdrawn from market due to TdP.

• Category 3: Drugs with a measurable incidence or numerous reports of TdP in humans.

• Category 4: Drugs with isolated reports of TdP in humans.

• Category 5: Drugs with no reports of TdP in humans.

We also considered a binary classifiers, where categories 2 and 3 were grouped as torsadogenic and categories 4 and 5 were grouped as non-torsadogenic. Results from this were not materially different, and can be found in the ESI spreadsheet under “Binary”.

We used EAD thresholds, APDs, diastolic calcium concentration, and hERG IC50/EFTPCmax to classify compounds into one of the four TdP risk categories described earlier. We tested two classification methods: linear discriminant analysis (LDA), and support vector machines (SVM).

Leave-one-out cross-validation and five-group cross-validation were used to check the robustness and accuracy of both the classifiers. We evaluated performance by calculating the errors in classification (how many risk classes away from the correct class a drug was classified as), and comparing these means of absolute error, E:

 
image file: c7tx00141j-t3.tif(3)
where, a is the actual category (i.e. the real risk category for the drug), c is the category assigned using the classification method, and N is the number of drugs in the dataset.

In LDA, the metrics in categories are assumed to follow a Normal distribution, and the points where the inferred distributions overlap are used as the category boundaries.60 In SVM, a hyperplane is used to separate the data points, and the optimal hyperplane is found by maximising the distance between the hyperplane and the closest points to it.61

In leave-one-out cross-validation one drug was removed from the dataset and the boundaries were re-calculated. The left-out drug was then placed into a risk category based on these new boundaries. Five-group cross-validation was also used: drugs were randomly assigned to five groups, and then each group was removed from the dataset and the classification boundaries were re-calculated. The drugs in the removed group were then classified into a risk category based on these new boundaries, and accuracy was calculated as above.

Results

EAD thresholds

In general, more torsadogenic drugs caused a decrease in the threshold of intervention required to provoke an EAD, i.e. they made the electrophysiology models more vulnerable to EADs. Full tables of EAD thresholds can be found in the ESI.

Some examples of EADs induced by the ICaL increase protocol are shown in Fig. 3. The low risk drug, nitrendipine, required a greater increase in ICaL conductance (27.23×) than the control (24.13×) to induce an EAD. Conversely, the high risk drug cisapride required a much smaller increase (8.23×) to induce an EAD. Similarly, models whose ion channel conductances had been modified to simulate the effects of cisapride required an INa inactivation curve shift of only 15.45 mV to produce an EAD, whereas at control a larger shift of 17.41 mV was required, and the low risk drug diltiazem required a shift of 17.49 mV, greater than control.


image file: c7tx00141j-f3.tif
Fig. 3 Some examples of EADs provoked using the L-type calcium increase protocol, as described in the methods section. The first EAD caused by the increasing intervention for each drug is highlighted in red. The less torsadogenic drug (nitrendipine) requires more provocation than the control to cause an EAD (i.e. its EAD threshold is higher), and the more torsadogenic drug (cisapride) requires less provocation (i.e. its EAD threshold is lower).

In general, ICaL-provoked EAD thresholds were lower for drugs which strongly block the INa and hERG channels. ICaL EAD thresholds were not linear with ICaL drug block, showing that multi-ion channel effects affected the thresholds. Less IKr block was required to cause an EAD for strong hERG blockers, except for verapamil, which is a very strong hERG blocker, but also strongly blocks ICaL. Thioridazine needed no hERG block to cause an EAD, despite being a strong ICaL blocker. This is probably due to thioridazine's very strong effect on hERG, with an IC50 value of 0.034 × EFTPCmax.

I Na inactivation curve shift thresholds were lower for drugs which strongly block hERG. Unlike the other EAD metrics, block of ICaL and INa did not increase the EAD threshold for the INa inactivation curve shift protocol. For example, verapamil, which is a strong hERG and ICaL blocker, had a low INa shift EAD threshold.

All three of the EAD metrics consistently had the lowest thresholds for the following drugs: azimilide, cisapride, ondansetron, terfenadine, and terodiline. Three drugs (ajmaline, quinidine, and thioridazine) produced EADs without any intervention required. With the exception of thioridazine (category 3), all of these drugs are in the highest risk category 2.

A table of EAD thresholds and other metrics can be found in the ESI.

Linkage analysis

The dendrogram from the linkage analysis of the hERG IC50/EFTPCmax metric, APD90, and the combined EAD metrics are shown in Fig. 4. The optimal leaf ordering algorithm was moderately successful in sorting the drugs in ascending order of torsadogenic risk. Therefore, a classification scheme based on the similarity of new compounds to existing compounds of known torsadogenic risk could be a useful tool. The ranking error measure for all metrics can be found in the ESI.
image file: c7tx00141j-f4.tif
Fig. 4 Left: Linkage analysis based on the hERG IC50/EFTPCmax metric. Centre: Linkage analysis based on the APD90 metric. Right: Linkage analysis based on the ICaL EAD threshold metric. Risk categories are indicated by colour: category 2 (dangerous) drugs are shown in red, category 3 in yellow, category 4 in black, and category 5 (safe) in green.

Pro-arrhythmic risk classification

Table 2 shows the mean absolute error, E, in classification for each of the proarrythymic risk markers and classification methods. Overall, there was good agreement between the different classification schemes and validation procedures. Full tables of all results can be found in the ESI.
Table 2 Mean absolute errors in classification, calculated using eqn (3). The “5 LDA” and “5 SVM” rows are the results from the 5-group cross-validation of linear discriminant analysis classification and support vector machines classification, respectively. The “1 LDA” and “1 SVM” rows are the results from leave-one-out cross-validation. The columns are arranged by the lowest sum of errors from the four classification measures. Yellow cells show metrics with low errors and purple cells have high errors. ‘L&S’ refers to the diastolic calcium concentration-APD50 combination metric from a paper by Lancaster & Sobie,39 ‘cqInward’ refers to the metric from Li et al.48 and ‘Redfern hERG Safety Factor’ refers to the hERG IC50/EFTPCmax metric proposed in Redfern et al.,36 as discussed above
image file: c7tx00141j-u1.tif


The leave-one-out analysis showed that the Support Vector Machine method of classification was more accurate for the Grandi Lancaster–Sobie metric, the Grandi APD90 metric, the O'Hara Cqinward metric, and the hERG IC50/EFTPCmax metric. For all other metrics, Linear Discriminant Analysis was more accurate.

Using LDA, the most accurate metrics were the O'Hara ICaL increase EAD metric and the O'Hara ICaL increase EAD metric combined with the Grandi APD90 metric. Using SVM, the most accurate metric was the Grandi Lancaster–Sobie metric, followed by Grandi APD90.

The least accurate metric using LDA was the cqInward metric. Using SVM, the least accurate metric was the IKr EAD protocol, followed by the INa inactivation curve shift EAD metric.

The best metric over all classification methods was the Grandi Lancaster–Sobie metric using SVM, and the worst over all classification method was the IKr EAD protocol using SVM.

The ICaL increase EAD metric correctly classified cibenzoline, desipramine, and fluvoxamine into categories 5, 4, and 4, respectively, while the APD90 metric mis-classified them into the more dangerous categories 4, 3, and 3, respectively, and the hERG IC50/EFTPCmax metric mis-classified all three drugs into the most dangerous category 2.

Discussion

In this paper, we have presented a novel method for predicting pro-arrhythmic risk by provoking EADs in computational models of cardiac cells, in combination with simulated ion channel block. As expected, EAD thresholds were usually lower for drugs which strongly block hERG, except when hERG block was combined with other ion channel block. ICaL block removed the effect of hERG block on the IKr EAD threshold, for example, verapamil, a strong ICaL blocker, required a large decrease in IKr current to cause an EAD despite being a strong hERG blocker. This indicates that for patients with LQT2 or other hERG mutations, drugs which block ICaL could be beneficial in preventing TdP. ICaL EAD thresholds were lower for drugs which strongly block both the INa and hERG channels. Therefore, for patients with increased ICaL activity, due to ICaL agonists or genetic mutations, drugs which block both INa and hERG could increase the risk of TdP. INa inactivation curve shift thresholds were not increased by ICaL block—for example, verapamil had a low INa EAD threshold despite strongly blocking ICaL. This indicates that, for patients with Brugada syndrome, the torsadogenic effects of hERG block cannot be ameliorated by ICaL block. ICaL increase as a torsadogenic risk metric showed clear improvement over the hERG ‘safety margin’ marker of IC50/EFTPCmax.

In general, results were consistent over classification methods and evaluation scheme, as shown by leave-one-out cross-validation and five-group cross-validation.

The ICaL increase EAD metric sorted cibenzoline, desipramine, and fluvoxamine into their correct risk categories, while both hERG IC50/EFTPCmax and Grandi APD90 classified them into more dangerous categories. All three of these drugs are relatively strong ICaL blockers. Cibenzoline is also a strong INa blocker and a weak IpNa blocker. Block of the persistent sodium or L-type calcium currents have been shown to suppress EADs.62,63 These results show that EAD metrics can accurately predict the torsadogenicity of drugs which block hERG and increase APD without causing TdP.

Amiodarone was persistently misclassified by every metric except for hERG IC50/EFTPCmax using LDA, and the IKr EAD metric and cqInward metric using SVM. Amiodarone weakly blocks several ion channels, but also has active metabolites that were not considered here.64 This combination of ion channel block leads to only a small increase in APD, and slightly decreased EAD thresholds. Ion channel block by amiodarone has been shown to be use-dependent, and there is evidence that amiodarone binds only to the open state of the INa channel. Therefore, a more detailed model of amiodarone binding kinetics may be required. However, despite being a Class III antiarrythmic, amiodarone poses a much lower, if still measurable, TdP risk than most other compounds in risk category one65 and may therefore be a candidate for re-classification into a lower clinical risk category.

Following APD-based metrics, an EAD-based approach to pro-arrhythmic risk prediction offers a mechanistic link between ion channel block and TdP.35,37 Building on previous EAD-based metrics,40,41 our method looks at causes other than increased drug concentration as a cause of TdP, and incorporates the action of underlying disease states. Linkage analysis may allow for more fine-grained risk ranking and be helpful in showing compounds with similar properties.

The separation of drugs into risk categories based on TdP clinical incidence is a difficult problem. TdP can only be diagnosed when a patient is being monitored on an electrocardiogram, so episodes of TdP may be missed, giving underestimates of incidence. In addition, there is little publicly-available information on drug prescription numbers, meaning that it is difficult to calculate the number of TdP cases per prescription/dose. The lack of information about TdP incidence per dose makes risk classes uncertain, and a drug prescribed to people who are more likely to get electrocardiograms will have an increase in the number of cases of TdP diagnosed for that drug compared to a similarly torsadogenic drug. Wiśniowska and Polak44 showed how different risk classifications have been given to the same drugs in different studies. Our strategy of using compounds with uncontroversial risk categorisations could ameliorate this problem, but the differences in risk category are fundamentally due to a lack of data on incidence.

The approach we have presented suggests a strategy by which risk assessments might be made for different patient subgroups (analogous to the different disease-mimicking EAD provoking interventions we applied). But as discussed, at present incidence data for the population as a whole is lacking, and this problem is exacerbated for smaller patient subgroups, making evaluation of our patient-group-specific predictions impossible at present. To improve these and other TdP prediction efforts, more data on TdP incidence rates will be needed.

One weakness of our study is the variety of sources from which the IC50 values were obtained. The differences in experimental protocol, temperatures, cell types and equipment used in these experiments may add uncertainty to our simulations,66 which we have not included in the computations presented here. A dataset from a single set of manual patch-clamp experiments with low variability would allow standardisation across pro-arrhythmic risk classification methods for a range of drugs.

There are a large number of available cardiac cell models. For the afterdepolarization aspects of this study, we have used only the O'Hara et al.43 model, because it includes the persistent sodium current. The addition of an appropriate late sodium current to existing models would be a useful extension to this work, to allow us to look at predictions from a wider range of models. This is work we are pursuing. We did not use the new dynamic hERG block model from Li et al. (2017)48 which may explain why our implementation of the Cqinward metric was not as successful as in that paper.

Our study does not look at pro-arrhythmic markers at the tissue or whole heart level. EADs are cell-level phenomena that interact with several other factors in the onset of TdP in tissue.67 The effects of spatial heterogeneity in ion channel expression and in fibre organisation on the arrhythmogenic effects of EAD susceptibility induced by drug block are likely to be significant in the translation from single-cell EADs to tissue- and organ-level effects.

Instead of using a simple pore-block model for drug interactions with ion channels, it would be interesting to account for different hERG binding kinetics, which have been explored in a recent study.48 The differences in binding may alter the effect of the drugs on action potential duration and EAD susceptibility.

Conclusions

We have proposed and investigated novel metrics for predicting drug-induced torsadogenic risk based on early-stage pre-clinical data on ion channel block. Our EAD-based metrics combine ion channel block data with disease states in order to predict increased EAD susceptibility in cell models, as a marker for pro-arrhythmic risk. Some EAD metrics were an improvement on the hERG IC50/EFTPCmax safety margin as a predictor of clinical incidence of Torsades de Pointes. The ICaL increase metric performed well, and was more accurate than hERG block alone, but was not as predictive as simpler measures such as simulated action potential duration we have published previously. This may be because different routes to EADs mimic different diseases in patient subgroups and represent risk well for these patients only, but evaluating this is difficult without further data on clinical TdP incidence rates.

Conflicts of interest

The authors have no conflicts of interest to declare.

Acknowledgements

Thanks to Martin Fink, Alan Garny, Penny Noble & Denis Noble for sharing their ideas on the generation of EADs and a paper draft. BM was funded by the EPSRC (Grant Number EP/G03706X/1). GRM gratefully acknowledges support from a Sir Henry Dale Fellowship jointly funded by the Wellcome Trust and the Royal Society (Grant Number 101222/Z/13/Z).

References

  1. P. J. Schwartz, L. Crotti and R. Insolia, Circ.: Arrhythmia Electrophysiol., 2012, 5, 868–877 CrossRef PubMed .
  2. A. Fabiato and P. Coumel, Cardiovasc. Drugs Ther., 1991, 5, 167–169 CrossRef CAS PubMed .
  3. J. L. Bauman, R. A. Bauernfeind, J. V. Hoff, B. Strasberg, S. Swiryn and K. M. Rosen, Am. Heart J., 1984, 107, 425–430 CrossRef CAS PubMed .
  4. R. L. Woosley, Y. Chen, J. P. Freiman and R. A. Gillis, J. Am. Med. Assoc., 1993, 269, 1532–1536 CrossRef CAS PubMed .
  5. J. Bryan, Pharm. J., 2011, 287, 511–512 Search PubMed .
  6. D. Rampe, M. L. Roy, A. Dennis and A. M. Brown, FEBS Lett., 1997, 417, 28–32 CrossRef CAS PubMed .
  7. Committee on Safety of Medicines and Medicines Control Agency, Curr. Probl. Pharmacovigilance, 2000, 26, 9–14 Search PubMed .
  8. S. S. Po, D. W. Wang, I. C. Yang, J. P. Johnson, L. Nie and P. B. Bennett, J. Cardiovasc. Pharmacol., 1999, 33, 181–185 CrossRef CAS PubMed .
  9. M. E. Curran, I. Splawski, K. W. Timothy, G. M. Vincent, E. D. Green and M. T. Keating, Cell, 1995, 80, 795–803 CrossRef CAS PubMed .
  10. European Medicines Agency, ICH Topic S7B, 2005 Search PubMed .
  11. P. Milberg and W. Haverkamp, Basic Res. Cardiol., 2005, 371, 365–371 CrossRef PubMed .
  12. S. Cobbe, Eur. Heart J., 2004, 6, 9–11 Search PubMed .
  13. N. Stockbridge, J. Morganroth, R. R. Shah and C. Garnett, Drug Saf., 2013, 36, 167–182 CrossRef PubMed .
  14. P. G. A. Volders, M. A. Vos, B. Szabo, K. R. Sipido, S. H. M. De Groot, A. P. M. Gorgels, H. J. J. Wellens and R. Lazzara, Cardiovasc. Res., 2000, 46, 376–392 CrossRef CAS PubMed .
  15. G. Frommeyer and L. Eckardt, Nat. Rev. Cardiol., 2015, 13, 36–47 CrossRef PubMed .
  16. M. A. Habbab and N. El-Sherif, Am. J. Med., 1990, 89, 241–246 CrossRef CAS PubMed .
  17. N. El-Sherif, R. H. Zeiler, W. Craelius, W. B. Gough and R. Henkin, Circ. Res., 1988, 63, 286–305 CrossRef CAS PubMed .
  18. B.-R. Choi, F. Burton and G. Salama, J. Physiol., 2002, 543, 615–631 CrossRef CAS .
  19. D. M. Roden, Eur. Heart J., 1993, 14, 56–61 CrossRef PubMed .
  20. K. Wemhöner, C. Friedrich, B. Stallmeyer, A. J. Coffey, A. Grace, S. Zumhagen, G. Seebohm, B. Ortiz-Bonnin, S. Rinné, F. B. Sachse, E. Schulze-Bahr and N. Decher, J. Mol. Cell. Cardiol., 2015, 80, 186–195 CrossRef PubMed .
  21. P. Kang, M. Liao, M. R. Wester, J. S. Leeder and R. E. Pearce, Circ.: Cardiovasc. Genet., 2013, 6, 490–499 CrossRef PubMed .
  22. I. Splawski, K. W. Timothy, N. Decher, P. Kumar, F. B. Sachse, A. H. Beggs, M. C. Sanguinetti and M. T. Keating, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 8089–8096 CrossRef CAS PubMed  ; discussion 8086–8.
  23. A. Mazur, D. M. Roden and M. E. Anderson, Circulation, 1999, 100, 2436–2442 CrossRef .
  24. N. S. Ghais, Y. Zhang, A. A. Grace and C. L.-H. Huang, Exp. Physiol., 2009, 94, 240–254 CrossRef CAS PubMed .
  25. C. T. January, J. M. Riddle and J. J. Salata, Circ. Res., 1988, 62, 563–571 CrossRef CAS PubMed .
  26. E. Marban, S. W. Robinson and W. G. Wier, J. Clin. Invest., 1986, 78, 1185–1192 CrossRef CAS PubMed .
  27. L. Schamroth, D. Krikler and C. Garrett, Br. Med. J., 1972, 1, 660–662 CrossRef CAS PubMed .
  28. P. Milberg, M. Fink, C. Pott, G. Frommeyer, J. Biertz, N. Osada, J. Stypmann, G. Mönnig, M. Koopmann, G. Breithardt and L. Eckardt, Br. J. Pharmacol., 2012, 166, 557–568 CrossRef CAS PubMed .
  29. P. L. Hedley, P. Jørgensen, S. Schlamowitz, J. Moolman-Smook, J. K. Kanters, V. A. Corfield and M. Christiansen, Hum. Mutat., 2009, 30, 1256–1266 CrossRef CAS PubMed .
  30. Q. Chen, G. E. Kirsch, D. Zhang, R. Brugada, J. Brugadak, P. Brugada, D. Potenza, A. Moya, M. Borggrefe, R. Ortiz-lopez, Z. Wang, C. Antzelevitch, R. E. O. Brien, E. Schulze-bahr, M. T. Keating, A. Jeffrey and Q. Wang, Nature, 1998, 392, 293–296 CrossRef CAS PubMed .
  31. L. Wang, X. Meng, Z. Yuchi, Z. Zhao, D. Xu, D. Fedida, Z. Wang and C. Huang, Cell. Physiol. Biochem., 2015, 36, 2250–2262 CrossRef CAS PubMed .
  32. A. Kagan, Z. Yu, G. I. Fishman and T. V. McDonald, J. Biol. Chem., 2000, 275, 11241–11248 CrossRef CAS PubMed .
  33. J. E. Talmadge, Circulation, 1999, 99, 529–533 CrossRef .
  34. H. Morita, J. Wu and D. P. Zipes, Lancet, 2008, 372, 750–763 CrossRef CAS .
  35. G. R. Mirams, Y. Cui, A. Sher, M. Fink, J. Cooper, B. M. Heath, N. C. McMahon, D. J. Gavaghan and D. Noble, Cardiovasc. Res., 2011, 91, 53–61 CrossRef CAS PubMed .
  36. W. S. Redfern, L. Carlsson, A. S. Davis, W. G. Lynch, I. MacKenzie, S. Palethorpe, P. K. Siegl, I. Strang, A. T. Sullivan, R. Wallis, A. J. Camm and T. G. Hammond, Cardiovasc. Res., 2003, 58, 32–45 CrossRef CAS PubMed .
  37. K. A. Beattie, C. Luscombe, G. Williams, J. Munoz-Muriedas, D. J. Gavaghan, Y. Cui and G. R. Mirams, J. Pharmacol. Toxicol. Methods, 2013, 68, 88–96 CrossRef CAS PubMed .
  38. G. R. Mirams, M. R. Davies, S. J. Brough, M. H. Bridgland-Taylor, Y. Cui, D. J. Gavaghan and N. Abi-Gerges, J. Pharmacol. Toxicol. Methods, 2014, 70, 246–254 CrossRef CAS PubMed .
  39. M. C. Lancaster and E. A. Sobie, Clin. Pharmacol. Ther., 2016, 100, 371–379 CrossRef CAS PubMed .
  40. B. Christophe, Pharmacol. Rep., 2013, 65, 1281–1293 CrossRef CAS PubMed .
  41. S. Morotti and E. Grandi, MethodsX, 2016, 4, 25–34 CrossRef PubMed .
  42. P. T. Sager, G. Gintant, J. R. Turner, S. Pettit and N. Stockbridge, Am. Heart J., 2014, 167, 292–300 CrossRef PubMed .
  43. T. O'Hara, L. Virág, A. Varró and Y. Rudy, PLoS Comput. Biol., 2011, 7, e1002061 Search PubMed .
  44. B. Wiśniowska and S. Polak, Drug Discovery Today, 2017, 22, 10–16 CrossRef PubMed .
  45. C. L. Lawrence, M. H. Bridgland-Taylor, C. E. Pollard, T. G. Hammond and J. P. Valentin, Br. J. Pharmacol., 2006, 149, 845–860 CrossRef CAS PubMed .
  46. Z. Liu, R. B. Williams and B. D. Rosen, J. Cardiovasc. Dis. Res., 2013, 4, 187–190 CrossRef PubMed .
  47. S. R. Mittal, J. Assoc. Physicians India, 2014, 62, 426–427 CAS .
  48. Z. Li, S. Dutta, J. Sheng, P. N. Tran, W. Wu, K. Chang, T. Mdluli, D. G. Strauss and T. Colatsky, Circ.: Arrhythmia Electrophysiol., 2017, 10, e004628 CrossRef CAS PubMed .
  49. E. Grandi, F. S. Pasqualini and D. M. Bers, J. Mol. Cell. Cardiol., 2010, 48, 112–120 CrossRef CAS PubMed .
  50. D. Bottino, R. C. Penland, A. Stamps, M. Traebert, B. Dumotier, A. Georgieva, G. Helmlinger and G. S. Lett, Prog. Biophys. Mol. Biol., 2006, 90, 414–443 CrossRef CAS PubMed .
  51. R. C. Elkins, M. R. Davies, S. J. Brough, D. J. Gavaghan, Y. Cui, N. Abi-gerges and G. R. Mirams, J. Pharmacol. Toxicol. Methods, 2013, 68, 112–122 CrossRef CAS PubMed .
  52. P. J. Noble, D. Noble and A. Garny, 2017, in preparation.
  53. G. R. Mirams, C. J. Arthurs, M. O. Bernabeu, R. Bordas, J. Cooper, A. Corrias, Y. Davit, S.-J. Dunn, A. G. Fletcher and D. G. Harvey, et al. , PLoS Comput. Biol., 2013, 9, e1002970 CAS .
  54. J. Cooper, R. J. Spiteri and G. R. Mirams, Front. Physiol., 2015, 6, 1–16 Search PubMed .
  55. G. Williams and G. R. Mirams, J. Pharmacol. Toxicol. Methods, 2015, 75, 10–16 CrossRef CAS PubMed .
  56. A. A. Cuellar, C. M. Lloyd, P. F. Nielsen, D. P. Bullivant, D. P. Nickerson and P. J. Hunter, Simulation, 2003, 79, 740–747 CrossRef .
  57. J. Cooper, S. McKeever and A. Garny, Proceedings of the 2006 ACM SIGPLAN Symposium on Partial Evaluation and Semantics-based Program Manipulation, New York, NY, USA, 2006, pp. 12–20 Search PubMed .
  58. S. D. Cohen and A. C. Hindmarsh, Comput. Phys., 1996, 10, 138–143 CrossRef .
  59. T. Galili, Bioinformatics, 2015, 31(22), 3718–3720 CrossRef CAS PubMed .
  60. A. L. Tarca, V. J. Carey, X.-W. Chen, R. Romero and S. Drăghici, PLoS Comput. Biol., 2007, 3, e116 CrossRef PubMed .
  61. N. Cristianini and J. Shawe-Taylor, An Introduction to Support Vector Machines: And Other Kernel-based Learning Methods, Cambridge University Press, New York, NY, USA, 2000 Search PubMed .
  62. C. Antzelevitch, L. Belardinelli, A. C. Zygmunt, A. Burashnikov, J. M. Di Diego, J. M. Fish, J. M. Cordeiro and G. Thomas, Circulation, 2004, 110, 904–910 CrossRef CAS PubMed .
  63. J. N. Weiss, A. Garfinkel, H. S. Karagueuzian, P.-S. Chen and Z. Qu, Heart Rhythm, 2010, 7, 1891–1899 CrossRef PubMed .
  64. M. R. Ghovanloo, M. Abdelsayed and P. C. Ruben, Front. Pharmacol., 2016, 7, 1–11 Search PubMed .
  65. Y. G. Yap and A. J. Camm, Heart, 2003, 89, 1363–1372 CrossRef CAS PubMed .
  66. G. R. Mirams, P. Pathmanathan, R. A. Gray, P. Challenor and R. H. Clayton, J. Physiol., 2016, 594, 6833–6847 CrossRef CAS PubMed .
  67. N. Vandersickel, E. Van Nieuwenhuyse, G. Seemann and A. V. Panfilov, Front. Physiol., 2017, 8, 404 CrossRef PubMed .

Footnote

Electronic supplementary information (ESI) available. See DOI: 10.1039/C7TX00141J

This journal is © The Royal Society of Chemistry 2017