First universal pharmacophore model for hERG1 K+ channel activators: acthER

Journal of Molecular Graphics and Modelling 74 (2017) 153–170

Contents lists available at ScienceDirect

Journal of Molecular Graphics and Modelling journal homepage: www.elsevier.com/locate/JMGM

Topical Perspectives

First universal pharmacophore model for hERG1 K+ channel activators: acthER Serdar Durdagi a,∗ , Ismail Erol a,b , Ramin Ekhteiari Salmas a , Matthew Patterson c , Sergei Y. Noskov c,∗ a

Computational Biology and Molecular Simulations Laboratory, Department of Biophysics, School of Medicine, Bahcesehir University, Istanbul, Turkey Department of Chemistry, Gebze Technical University, Kocaeli, Turkey c Centre for Molecular Simulation, Department of Biological Sciences, University of Calgary, Calgary, Alberta, Canada b

a r t i c l e

i n f o

Article history: Received 19 August 2016 Received in revised form 28 March 2017 Accepted 29 March 2017 Available online 5 April 2017 Keywords: hERG1 K channel Pharmacophore modeling 3D-QSAR Molecular docking Molecular Dynamics (MD) Simulations Structure-based Pharmacophore Modeling

a b s t r a c t The intra-cavitary drug blockade of hERG1 channel has been extensively studied, both experimentally and theoretically. Structurally diverse ligands inadvertently block the hERG1 K+ channel currents lead to drug induced Long QT Syndrome (LQTS). Accordingly, designing either hERG1 channel openers or current activators, with the potential to target other binding pockets of the channel, has been introduced as a viable approach in modern anti-arrhythmia drug development. However, reports and investigations on the molecular mechanisms underlying activators binding to the hERG1 channel remain sparse and the overall molecular design principles are largely unknown. Most of the hERG1 activators were discovered during mandatory screening for hERG1 blockade. To fill this apparent deficit, the first universal pharmacophore model for hERG1 K+ channel activators was developed using PHASE. 3D structures of 18 hERG1 K+ channel activators and their corresponding measured binding affinity values were used in the development of pharmacophore models. These compounds spanned a range of structurally different chemotypes with moderate variation in binding affinity. A five sites AAHRR (A, hydrogen-bond accepting, H, hydrophobic, R, aromatic) pharmacophore model has shown reasonable high statistical results compared to the other developed more than 1000 hypotheses. This model was used to construct steric and electrostatic contour maps. The predictive power of the model was tested with 3 external test set compounds as true unknowns. Finally, the pharmacophore model was combined with the previously developed receptor-based model of hERG1 K+ channel to develop and screen novel activators. The results are quite striking and it suggests a greater future role for pharmacophore modeling and virtual drug screening simulations in deciphering complex patterns of molecular mechanisms of hERG1 channel openers at the target sites. The developed model is available upon request and it may serve as basis for the synthesis of novel therapeutic hERG1 activators. © 2017 Published by Elsevier Inc.

1. Introduction Ion channels are membrane bound proteins that act as gated pathways for the movement of ions across the cell membrane and they play important roles in the physiology of all cells [1–4,49–51]. A large number of inherited or drug-induced human diseases are associated with malfunction in ion transport across ion channels and generally known as channelopathy [5]. Channelopathies mostly result from a loss of channel function, making channel acti-

∗ Corresponding authors. E-mail addresses: [email protected] (S. Durdagi), [email protected] (S.Y. Noskov). http://dx.doi.org/10.1016/j.jmgm.2017.03.020 1093-3263/© 2017 Published by Elsevier Inc.

vation an attractive approach to regain the function. Several classes of potassium channels are involved in regulating the heart rate by setting both the amplitude and duration of the action potential, and the resting membrane potential [5–7,52]. Abnormalities in function of these ion channels, due to inherited mutations or pharmacological blockage, can prolong the duration of the action potential, leading to the development of severe arrhythmias (i.e. Long QT syndromes – LQTS) [8,9]. Genetic analysis has revealed that mutations in potassium channels such as the human ether-a-go-go related gene (hERG) and KvLQT1 form the molecular basis of LQTS [10]. The hERG1 channels are among prime targets for the pharmacological management of arrhythmias [11,12]. One of the promising strategies in avoiding drug-induced QT prolongation due to hERG1 blockade is to develop set of small molecules that may promote

154

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

channel openings e.g. activate hERG1 current [5]. Potassium (K+ ) channels can be activated by a potent and selective class of compounds known as channel openers (aka activators) [49]. Most of the drugs that enhance the current of ion channels were discovered by mandatory screening of novel drug candidates for blocking potencies [13–16]. The ideal opener should be capable of augmenting currents without marked alterations or kinetics or, in the case of gating modifiers, as a result of faster activation, slower deactivation, altered voltage-dependence, or a combination of these mechanisms [17–19]. In order to test the potential clinical advantage of this approach, channel-specific activators would be needed. Unfortunately to date only limited number of potent hERG1 K+ channel activators are known [20]. Although activators represent a new road to the development of anti-arrhythmic drugs, the molecular mechanisms governing their action remain to be discovered [49]. A thorough examination of the effects of hERG1 activators in vivo and in intact cardiac tissue has not yet been performed [45]. The information about the structural mechanism of actions [53] of these drugs is very limited and highly controversial. It is known that the hERG1 potassium channel openers (otherwise known as hERG1 channel agonist) are not uniform in their biophysical mechanism of action. While some drugs shift the V0.5 of voltage-dependent activation to more hyperpolarized potentials, some drugs increase the hERG1 current by shifting the voltage dependence of C-type of inactivation to more depolarizing potentials [15,53]. Some familial mutations in the hERG1 result in channels with voltage–dependent inactivation thresholds shifted substantially to the right (depolarized potentials) [5,21]. This shift can result in an increase in the time-dependent current at the expense of a decrease in the anti-arrhythmic tail current. These mutations can lead to familial Short-QT Syndrome (SQTS) [22]. Indeed drugs that increase the magnitude of the tail current or slow its deactivation can be anti-arrhythmic [23]. Different mechanisms of action have been suggested for several hERG1 activators. For example, RPR260243 activator induces two distinct changes in hERG1 gating: (i) a slowing of deactivation and (ii) a positive shift in voltage-dependent inactivation [15]. In contrast, NS1643 and NS3623 openers enhance the hERG1 current by reducing the channel inactivation threshold [24]. The PD-307243 hERG1 activator enhances the current by slowing the rates of both inactivation and deactivation [25]. The mechanism of action for the hERG1 activator PD-118057 is still unknown. Zhou et al. reported that it might only reduce channel inactivation [20]. Perry et al. investigated the putative binding site and mechanism of action of the hERG1 activator RPR260243 [26]. They reported that the putative binding site of RPR260243 is distinct from the hERG1 binding site of the blockers. Perry and colleagues suggested that RPR260243 alters the deactivation and inactivation gating of the hERG1 by interaction with a cluster of residues located in the neighbourhood of the cytoplasmic ends of the S5 and S6 domains of a single channel subunits. In a previous study, our group also endorse the possibility of multiple binding sites for NS1643 activator [19,27]. Binding of NS1643 to one site- the central cavity- appears to mediate pharmacologic block of hERG1. We have also defined a second potential binding site in the neighbourhood of E544, which appears to be involved slowing of deactivation and shifting the voltage-dependence of activation. Therefore, determining the mechanisms of action and binding site information for known hERG1 channel modulators is important for the design of robust novel channel activators. However, X-ray structures of most of K channels have yet to be solved. Thus, either homology models of these proteins (through the use of a template from another solved structure) or ligand-based models have to be constructed. Hence, the aim of this study is to attempt to construct and validate a universal pharmacophore model for hERG1 channel openers. To the best of our knowledge, currently there is no any ligand-based modeling study for hERG1 channel openers using

their 3D minimized structures and their corresponding EC50 values. Thus, results are quite striking and derived information from the pharmacophore model can be used for (i) understanding the general pharmacophore sites of hERG1 channel openers; (ii) the design of new hERG1 channel openers with enhanced activity and other tailored properties. 2. Methods 2.1. Selection of hERG1 channel openers for the pharmacophore model generation The binding interactions between a protein and a ligand depend on structural and chemical complementarity between them. In this study, we used PHASE from Maestro molecular modeling package for generation of pharmacophore hypotheses. PHASE classifies these chemical features as hydrogen bond acceptors (labeled as A), hydrogen bond donors (labeled as D), hydrophobic groups (labeled as H), negatively charged groups (labeled as N), positively charged groups (labeled as P) and aromatic rings (labeled as R). In order to define a common pharmacophore or hypothesis, PHASE employs an analysis of k-point pharmacophores derived from the conformational sets of active compounds and then identifies all spatial arrangements with pharmacophore features that are shared by those compounds. In developing an universal pharmacophore model for hERG1 channel openers, our search was included compounds that activate hERG1 channel from available literature data. Table 1 includes 18 compounds, which are used in the development of pharmacophore model hypotheses. These compounds span a range of structurally different chemotypes, with quite varying binding affinities between 0.52–79.4 ␮M. For instance, compound 7 is over 2.2 log units more active than compound 6. The compounds in Table 1 were generated from fragments using the LigPrep toolset from the Schrodinger’s Maestro molecular modeling package [28]. After initial construction, energy minimization is performed for each compound. This was completed using the steepest descent (SD) and a conjugate gradient (CG) approaches until the gradient converged at 5.00 J/mol Å. In order to survey the potential energy surface, conformation searches were conducted on all minimized compounds. Conformers were created by the variation in the dihedral angles or atomic positions. To sufficiently search conformational space, the number of trials was set to 1000. If the trial conformation had an energy decrease from the previous iteration, it became an accepted conformation. In the end, upwards of 100 conformers were generated for each of the compound. The entire minimization process was performed by the Monte Carlo Multiple Minimum (MCMM) conformational search analysis, coupled with the MMFFs force field, implemented by the Maestro molecular modeling package [29]. 2.2. 3D QSAR analysis of the selected compounds Once the compounds 1–18 were energy minimized and had conformers generated, they were randomly split into around 75% training and 25% internal test sets for the 3D QSAR analysis, based on the PHASE pharmacophore modeling module [30]. More specifically, compounds 1 through 14 formed the selected training set, while 15 through 18 became the internal test set. Within the training set, a wide range of activities was present, with lowest pEC50 equal to 4.1 (from compound 6) and the highest being 6.28 (from compound 7). An activity threshold value was set to be between the pEC50 values of 4.8–5.3 (i.e., compounds that have pEC50 values less than 4.8 are considered as inactive and compounds that have pEC50 values higher than 5.3 are considered active). From this, the total number of actives was equal to 4, while the total number of inac-

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

155

Table 1 Activators used to develop the hERG1 K channel opener pharmacophore models. In the PHASE analysis, compounds (1–14) and (15–18) were the training and test sets, respectively. (Activity values were taken from literature [13,20,26,43–47]). Activity EC50 (␮M)

Observed Activity pEC50

Predicted Activity pEC50

Difference

LUF-7389

4.60

5.34

5.24

0.10

2

LUF-7346

12.00

4.92

5.22

−0.30

3

T531-332M

30.00

4.52

4.62

−0.10

4

RPR-260243

15.00

4.82

4.70

0.12

5

NS-1643

10.50

4.98

4.90

0.08

6

NS-3623

79.40

4.10

4.13

−0.03

7

Mallotoxin

0.52

6.28

6.31

−0.03

8

AZSMO-23

11.20

4.95

4.91

0.04

9

T531-297

13.60

4.87

4.93

−0.06

10

T531-400

15.40

4.81

4.80

0.01

Compound Number

Compound Name

1

2D Structure

156

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

Table 1 (Continued) Activity EC50 (␮M)

Observed Activity pEC50

Predicted Activity pEC50

Difference

T531-343

16.10

4.79

4.92

−0.13

12

ML-T531

3.10

5.51

5.29

0.22

13

T531-362

12.20

4.91

4.71

0.20

14

T531-367

30.00

4.52

4.67

−0.15

15

LUF-7244

4.60

5.34

5.29

0.05

16

T531-298A

10.70

4.97

5.01

−0.04

17

T531-364(R)

11.50

4.94

4.93

0.01

18

T531-346

7.30

5.14

5.07

0.07

Compound Number

Compound Name

11

2D Structure

tive compounds was equal to 4, and 10 molecules with moderate hERG1 activities. Furthermore, in developing the pharmacophore modeling, the following setting was considered: for a hypothesis to be accepted it must match at least 2 out of the 4 actives. To establish the optimal number of components to be used, R2 vs partial least square (PLS) number graphs were obtained (Supplementary material, Fig. S1). The number of components used was selected based on the convergence of the R2 values.

2.3. Validation of the pharmacophore models The pharmacophore models for hERG1 activators developed from the compounds in Table 1, was validated using an external test set compounds (19–21, Table 2). Compounds 19–21 were chosen because they are known hERG1 channel activators, with known EC50 values. These structures were constructed, minimized, and

had conformers generated based on the same procedure mentioned previously for training set and internal test set compounds. 2.4. Molecular docking Two different molecular docking programs and five different algorithms (Glide/Standard Precision (SP), Glide/Extra Precision (XP), Glide/Induced Fit Docking (IFD), Glide/Quantum-Polarized Ligand Docking (QPLD) and GOLD) [31–34] were used. (i) Maestro Docking Tools: In Glide/SP, Glide/XP, Glide/IFD, and Glide/QPLD docking calculations, default settings were used. Glide/IFD was used to take into account possible conformational sampling both for binding site of the receptor and ligand. Glide/QPLD was considered to take into account charge polarization around binding site by the calculation of quantum mechanics (QM) charges of ligands with the 6–31G*/LACVP* basis set, B3LYP density functional, and ‘Ultrafine’ SCF accuracy level. (ii) GOLD program was used to dock

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

157

Table 2 Observed and predicted pEC50 values of the external test set compounds 19–21. (Activity values were taken from literature [13–15,17,23,24,26,46,48]). Compound Number

Compound Name

19

2D Structure

Activity EC50 (␮M)

Observed Activity pEC50

Predicted Activity pEC50

Residual

ICA-105574

0.50

6.30

4.75

1.55

20

KB-130015

12.00

4.92

4.92

0.00

21

PD-118057

2.90

5.54

5.04

0.50

2.5. Molecular dynamics (MD) simulations

Fig. 1. AAHRR 5-sited pharmacophore model laid over one of the most active compound, LUF-7244 (15). (Grey = carbon, white = hydrogen, red = oxygen, blue = nitrogen).

same ligands to cross-validate performance of Glide/SP, Glide/XP, Glide/IFD, and Glide/QPLD. The GOLD program (version 5.1) was used with two default docking scores (the GOLD Fitness and ChemScore). The dockings produced using the GoldScore function were re-scored and re-ranked with the optimized ChemScore function. Partial flexibility of the receptor was gained with selected amino acids residues at the selected binding sites. The default genetic algorithm parameters were used.

Compound NS1643 (5) and its derivatives E1, E2, and E5 (see Fig. S2 and Table S1 at the Supplementary materials) were docked to the hERG1 open-state S1-S6 model and their topdocking poses at the E544 site were selected and submitted to 50 ns MD simulations using Desmond package [31]. The atomic interactions of the protein, ligands and water atoms were calculated by OPLS 2005 force field [32]. The systems were embedded into a dipalmitoylphosphatidylcholine (DPPC) membrane bilayer, and then placed in an orthorhombic box with layers of explicit TIP3P water molecules of 10 Å thickness. The particle-mesh Ewald method [33] was implemented to calculate the long-range electrostatic interactions. A cut-off radius of 9.0 Å was used for van der Waals and Coulombic interactions. Nose–Hoover thermostat [34] and Martyna–Tobias–Klein protocols [35] were used to restrain the temperature and pressure of the systems at 310 K and 1.01325 bar, respectively. The systems were minimized for a maximum of 5000 iterations, and a convergence threshold of 1 kcal mol−1 Å−1 was implemented. Then the systems were relaxed using a step-wise procedure and gradually equilibrated. Finally, production run for MD simulations were carried out for each system using 2.0 fs as time step. 2.6. E-Pharmacophore screening The advantages of both ligand and structure-based approaches can be achieved by combining these two different methods into

Fig. 2. The predicted pEC50 values versus the observed pEC50 values of the molecules used to generate the 5-site pharmacophore AAHRR model. (Training set R2 = 0.92, test set R2 = 0.95).

158

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

Fig. 3. One of the lowest active compound NS3623 (5 pEC50 observed = 4.10; pEC50 predicted = 4.09), aligned to AAHRR pharmacophore. The cubes represent favourable and unfavourable contributions to predicted activity. (a) Positive charge (blue is favourable, orange is negative). (b) Negative charge (red is favourable, yellow is unfavourable). (c) Electron withdrawing groups (Red is favourable, cyan is unfavourable). (d) Hydrogen bond donor (blue is favourable, orange is unfavourable). (e) Hydrophobic/non polar groups (green is favourable, purple is unfavourable).

E-Pharmacophore modeling (i.e., structure-based pharmacophore modeling). In this method, energetically minimized pharmacophore models are derived by generating docking poses at the binding sites of the hERG1 K+ channel. For this aim, two different schemes were applied. In the first part, a well-known hERG1 activator NS1643 was docked into three different activation sites of the open state model of hERG1 channel: K525, E544, and N629. Resulted docking poses were used to generate pharmacophore models. 9 different E-pharmacophore models were derived: 4 for the K525, 3 for the E544, and 2 for the N629 activator binding sites. Then, Otava Tangible Ligand Database (∼2.6 million compounds) was screened with generated pharmacophore models. Moreover, representative structure from collected trajectory frames of hERG1/NS1643 complex throughout MD simulations was used for the derivation of E-Pharmacophore models. Ligand (NS1643) position at the binding cavity was used to generate Epharmacophore model, and this model is then used to screen same Otava Tangible ligand library.

3. Results and discussion 3.1. Derivation of pharmacophore models The resulting hERG1 pharmacophore models, generated by the compounds in Table 1 with the highest Q2 and R2 , can be seen in Fig. 1 using compound 15 as a template. This 5-sited AAHRR model (acthER) is characterized by the presence of five different functional groups: two aromatic (R), two hydrogen-bond acceptors (A), and one hydrophobic moieties (H); with the geometries illustrated as shown. Distances between each site has been represented at Fig. S3 at the Supplementary materials. Fitness values of the training and test set compounds with respect to the AAHRR pharmacophore model has been given in Table S2 at the Supplementary Materials. While Fig. 2 illustrates the observed and predicted activities of training and test set compounds, Table 3 summarizes the obtained statistical results. The AAHRR pharmacophore model was the best predictor of pEC50 val-

Fig. 4. One of the highest active compound ML-T531 (12 pEC50 observed = 5.51; pEC50 predicted = 5.33), aligned to AAHRR pharmacophore. The cubes represent favourable and unfavourable contributions to predicted activity. (a) Positive charge (blue is favourable, orange is negative). (b) Negative charge (red is favourable, yellow is unfavourable). (c) Electron withdrawing groups (Red is favourable, cyan is unfavourable). (d) Hydrogen bond donor (blue is favourable, orange is unfavourable). (e) Hydrophobic/non polar groups (green is favourable, purple is unfavourable).

2 XP IFD QPLD

0.449 2.503 2.166

GOLD

0.483

-0.10 -0.356

Phase

-1.457

-0.77

2.03 0.10

-0.57 -0.213

0.12

9

-0.841 -0.256 -0.896

0.18 -0.066

10

1.865 1.627

SP

-0.005

8

0.433 0.275 0.073 0.40

1.255 0.03

7

Phase

0.18

-0.257 -0.567

6

GOLD

-0.058

0.358

-0.963

1.46

1.558 1.378

0.114 0.32

5

QPLD

0.607 0.174

Derivave Number

IFD

2.453 0.12

4

XP

-2.13 -1.998

-2.039

1.035

-0.085 -0.215

10

10

0.296 0.804 1.193 0.40

1.44

0.10

-0.093 -0.052

-1.718

1.106 0.36

9

9 -0.234

1.193 0.738 0.616 0.18

SP

-1.772 -1.82

3

-0.807

3.013 3.044

8

8

1.883 0.78 1.093 1.377 0.414 0.12

-0.105 -0.621 -1.094 -0.10

7

1.231

0.771 0.893 1.03 0.311 0.03

6

2.914

Phase

7 0.308 0.10

GOLD

2.19

0.32

5

QPLD

6 -0.10

1.608 1.494 1.971 2.052

-0.874

1.02

2.827 2.012

0.36

0.36

4

5 0.03

Derivave Number

IFD

2.783

0.127

-0.121

-1.438 -1.499 -1.598 -0.821

XP

0.377 0.16 0.584 0.715 0.32

-0.913 -1.029

3

1.605

SP

4 -1.204 -0.983 -0.578 -0.332

-1.189

0.26 0.744

1.547 1.526

0.029 0.17 0.36

3

3 -0.368

-1.81

0.351 0.26

1.858

-0.717

-1.262

-0.977 -1.82 -1.58 -0.453 -0.731

-2.029

2

2

1.805

0.26 -1.898

-2.784 -3.405 -3.12

0.796

0.18 0.11

2.06 1.473

2.202 1.781

0.314 0.801 0.40

1.685

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

-0.869

-0.547 -0.396 -0.58 -0.461

0.002 0.13 0.051 0.11

1

-0.776 -1.323

0.594 0.256 0.11

1

1

-1.708 -1.956

Predicted Binding Energy (kcal/mol)

Predicted Binding Energy (kcal/mol)

0.283

-1.106 -1.52

DERIVATIVE NUMBER

3.414

159

Fig. 5. Comparison of predicted binding energies of NS1643 derivatives from 10 enumaration site (please see Table S4 for structures) using developed 5 sites pharmacophore model AAHRR.5266 by PHASE, and docking scores of NS1643 and its derivatives at the hERG1 potential binding sites (top, K525; mid, E544; bottom, N629) using Glide/SP, Glide/XP, Glide/IFD, Glide/QPLD, and GOLD molecular docking algorithms.

PREDICTED BINDING ENERGY (KCAL/MOL)

160

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

Fig. 6. RMSDs of ligands (left) and RMSF of protein diagrams (right) for the studied MD systems (compounds NS1643, E1, E2, and E5 at the active sites of the protein). Ligand RMSD indicates how stable the ligand is with respect to the protein and its binding pocket. ‘ProFit’ shows the RMSD of a ligand when the protein-ligand complex is first aligned on the protein backbone of the reference and then the RMSD of the ligand heavy atoms is measured. If the values observed are significantly larger than the RMSD of the protein, then it is likely that the ligand has diffused away from its initial binding site. ‘LigFit’ shows the RMSD of a ligand that is aligned and measured just on its reference conformation. This RMSD value measures the internal fluctuations of the ligand atoms.

Table 3 Summary of the PHASE statistical parameters for the top 5-sited hypothesis, AAHRR. SD R2 F P RMSE Q2 Pearson-R

0.155 0.9218 64.8 8179E-07 0.0446 0.9206 0.9764

ues. The Q2 and R2 from this hypothesis was both found as 0.92. Secondly, this pharmacophore was able to predict the pEC50 value within 0.30 logarithm units of all (100%) compounds. Statistical parameters of other successful hypotheses generated throughout the pharmacophore modeling search analyses has been given at Table S3, Supplementary Materials. In our previous study, we have used calculated binding scores (from induced fit docking (IFD)/Glide XP scores at S4S5 linker site) of each compound as calculated binding affinity and these values were considered at the construction of 3D-QSAR for the synthesized 24 NS1643-like structures from our group [36,38]. In that study, the developed best hypothesis was 5-sited AADHR model. So, 4 out of 5 sites for the current developed model using experimental binding affinities of known hERG1 activators and previously developed hypothetical model using converted docking scores as calculated binding affinity are common. Figs. 3 and 4 show both the positive and negative contributions of the different moieties used by the pharmacophores in predicting the activity of each compound, illustrated in these cases, by the low active compound NS3623 (5, pEC50 = 4.10) and the high active compound ML-T531 (12, pEC50 = 5.51). These fig-

ures help to visually showcase how the pharmacophore predicts pEC50 values. These two compounds not only had the low and high observed pEC50 values, but also the low and high predicted pEC50s. In Fig. 3, the first thing to notice is that the pharmacophore model does not align very well with compound 5, the aromatic rings are not centred with pharmacophore R, the hydrophobic as well as hydrogen bond accepting pharmacophores do not fall on any hydrophobic functional groups and hydrogen bond acceptors, respectively. It is clear that compound 5 will have a low predicted pEC50 (4.09) (Fig. 3). In addition to poor alignment with the pharmacophores for compound 5, Fig. 3 also illustrates favourable and unfavourable contributions to hERG1 opener activity through the use of contour maps. Overall compound 5 shows either no or very small favourable contributions to activity in the form of electrostatic interactions, hydrogen bonding, electron withdrawing groups, or stabilizing hydrophobic moieties. In fact there are significant unfavourable contributions to activity, particularly in the case of both electron withdrawing and hydrophobic groups (Fig. 3c,e). On the other hand, the pharmacophores align very well with the functional groups (Fig. 4) for active compound 12. There are numerous favourable contributions to binding including stabilization due to negative charge and positions of electron withdrawing groups (Fig. 4a–e). The number of favourable contributions quite clearly outweigh the unfavourable interactions such as hydrophobic/steric penalties. Next, generated pharmacophore model was used to predict activities of 9 known hERG1 activators (R-L3, BMS-204352, Retigabine, PD-307243, A-935142, Diflunisal, Niflumic Acid, Flufenamic Acid, and Zoxazolamine). Although these molecules are wellknown hERG1 activators their EC50 values are not reported yet in

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

161

Fig. 7. Per-residue interactions of the compound E1, E5, NS1643 and E2 in complexes with the open state hERG1, in which the occupancy interaction values of the main amino acids around the ligands are calculated along the simulations.

162

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

Fig. 8. E-Pharmacophore models were generated by docking simulations of NS1643 at the K525 site. These E-pharmacophore models were used to screen ∼2.6 million compounds. Best 1000 compound ranked respect to their fitness score and these compounds were then docked to K525 binding site. Docking scores were plotted respect to the fitness scores.

Fig. 9. E-Pharmacophore models were generated by docking NS1643 at the E544 site. These E-pharmacophore models were used to screen ∼2.6 million compounds. Best 1000 compound ranked respect to their fitness score and these compounds were then docked to E544 binding site. Docking scores were plotted respect to the fitness scores.

literature. Thus, our aim was to test if the derived ligand-based model could identify these molecules as hERG1 activator as well as to predict their corresponding EC50 values. 2D structures of these compounds are collected in Table S4, Supplementary materials. Predicted pEC50 values were used to calculate binding free energies (Go ). These compounds were also docked to three previously

identified activator binding sites of hERG1 open channel centered at residues K525, E544 or N629, respectively. Predicted activity and binding free energies were given in Table S4. Predicted binding free energies were also compared with docking scores at three binding sites.

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

163

Table 4 Top fitted compounds and their docking scores respect to pharmacophores generated by docking poses at the K525 site. Site

Docking Pose

Hits

K525

1st

2nd

3rd

4th

2D Structure of Hits

Code

Fitness

Docking Score (kcal/mol)

1

8881797

2.594

−5.148

2

1102487

2.553

−7.174

3

4116631

2.551

−5.232

1

5439517

2.394

−4.986

2

5650867

2.348

−5.413

3

8850560

2.337

−5.143

1

5438338

2.564

−4.851

2

5436447

2.533

−4.737

3

5436588

2.524

−4.794

1

8882655

2.631

−5.417

2

8882656

2.620

−5.025

3

8882657

2.613

−5.028

3.2. Validation of pharmacophore models In order to validate the 5-sited AAHRR pharmacophore model, an external test set was also used. The results from these predictions are collected in Table 2. The derived 5-sited pharmacophore model performed reasonably well, predicting the pEC50 of KB-130015 (20)

exactly same as experimental value. Also the predicted pEC50 of the larger compound, 21 was well within 0.5 logarithmic units. The predicted pEC50 value for compound 19 by the model is overestimated, but still prediction was within around 1.5 logarithmic units. In order to test if our model (AAHRR) shows selectivity for hERG1 activators, non-blockers and blockers; 5 hERG blockers for each log

164

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

Table 5 Top fitted compounds and their docking scores respect to pharmacophores generated by docking poses at the E544 site. Site

Docking Pose

Hits

E544

1st

3rd

5th

2D Structure of Hits

Code

Fitness

Docking Score (kcal/mol)

1

6478349

2.294

−4.747

2

5621456

2.292

−6.078

3

5621430

2.290

−5.425

1

8882655

2.516

−5.304

2

4118938

2.504

−6.176

3

4121307

2.503

−6.080

1

8881797

2.577

−6.686

2

8881801

2.568

−6.092

3

8882655

2.536

−5.280

unit from 9 to 5 (e.g. ≈pIC50 9 (5 compounds); ≈pIC50 8 (5 compounds); ≈pIC50 7 (5 compounds); ≈pIC50 6 (5 compounds), ≈pIC50 5 (5 compounds)), 20 non-blockers with pIC50 < 4.5 were collected from ChEMBL database; as well as all activators that are used in generation and validation stages of our pharmacophore model AAHRR are considered. 5-sited pharmacophore model (AAHRR) failed as expected to predict range of 7, 8 and 9 pIC50 values, predicting their activities between 2 and 4 log-unit difference. Since generated model was build using known hERG1 activator structures, failure of model during prediction of the activities of the hERG1 blockers seems to be natural. Blockers that have pIC50 values between 6 and 5, predicted well. All blockers in these range, predicted within 1.2 log-unit. This finding shows that, weak binding hERG1 blockers may have similar pharmacophore sites with hERG1 activators (Tables S5–S7, Supplementary materials). 3.3. Combination of ligand-based models with receptor-based models Derivatives of one of the well-known hERG1 activators (NS1643) were generated using CombiGlide module and fragment library (∼4000 small-compound database) of Schrodinger Molecular mod-

eling package [35]. Together with different binding combinations of the fragments with 10 enumeration sites (Fig. S2, Supplementary materials) for the chosen template (NS1643), the total number of derivatives was around 40000. Predicted activities of these compounds were docked at the previously developed and reported S1-S6 hERG1 open state receptor model from our group [36,53]. Subsequently, selected top-scored molecules were also screened and compared at the developed ligand-based model AAHRR. The most active NS1643 derivatives from each of the experimentally identified binding sites implicated in channel activation are shown in Table S1 at the supplementary materials. These compounds were docked at the open state of hERG1 channel and docking scores were compared with wild-type NS1643. According to the previous binding site reports from the literature, three major binding pockets were considered: (i) E544 site (i.e., binding cavity within 15 Å from E544); (ii) N629 site (i.e., binding cavity within 15 Å from N629); (iii) K525 site (i.e., binding cavity within 15 Å from K525). Two different docking programs (Glide and GOLD) were used. 4 different algorithm of Glide (SP, XP, IFD, and QPLD) were considered. Fig. 5 shows residual docking scores of wild NS1643 and 10 derived NS1643 analogues that showed highest predicted pEC50 values from 10 enumeration sites. Predicted pEC50 values

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

165

Table 6 Top fitted compounds and their docking scores respect to pharmacophores generated by docking poses at the N629 site. Site

Docking Pose

Hits

N629

2nd

3rd

2D Structure of Hits

Code

Fitness

Docking Score (kcal/mol)

1

7263412

2.496

−6.204

2

7263263

2.494

−6.496

3

7263473

2.477

−7.227

1

5651383

1.991

−5.438

2

5438351

1.962

−5.906

3

5438332

1.941

−6.388

Fig. 10. E-Pharmacophore models were generated by docking NS1643 at the N629 site. These E-pharmacophore models were used to screen ∼2.6 million compounds. Best 1000 compound ranked respect to their fitness score and these compounds were then docked to N629 binding site. Docking scores were plotted respect to the fitness scores.

of these compounds by constructed 5-sited pharmacophore model AAHRR have been given in Table S1 at the supplementary materials. Comparison of predicted binding energies of NS1643 derivatives from 10 enumeration site using developed 5 sites pharmacophore model AAHRR by PHASE, and docking scores of NS1643 and its

derivatives at the hERG1 potential binding sites using five different molecular docking algorithms are summarized in Fig. 5. Results showed that derivative #2 (E2) shows a consistent common profile in both the five different docking algorithms and ligand-based model prediction at the binding sites of E544 and K525 (Fig. 5).

166

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

Fig. 11. E-Pharmacophore models were generated using representative NS1643/hERG1 complex structure from MD simulations at the E544 site. These E-pharmacophore models were then used to screen ∼2.6 million database compounds. Best 1000 compound ranked respect to their fitness score and these compounds were docked to E544 binding site. Docking scores were plotted respect to the fitness scores.

3.4. MD simulations

3.5. E-Pharmacophore screening

MD simulations of the NS1643 (5) and its derivatives E1, E2, and E5 in binding to the open state hERG1 model were performed to better understand the binding mechanisms as well as to identify crucial amino acids in binding. Fig. 6 represents RMSDs of ligands and RMSF of protein diagrams for the studied MD systems. Internal fluctuations of the ligand atoms were found small. As expected, while linker regions have higher RMSF values, trans-membrane domains show rigid conformations. The occupancy interactions for compounds NS1643, E1, E2 and E5 are given in Fig. 7, in which, not only the hydrogen and hydrophobic interactions are identified, but also water bridges, have formed between the ligand and protein are mentioned. In the case of the compound E1, it is found that Tyr542, and Ser654 are among the highly pronounced amino acids in the binding domain via forming of hydrophobic and hydrogenbond interactions, respectively. In addition, the other amino acids, which have participated in the stabilization of the ligand are determined in the diagram. In the case of the compound E5, it formed a strong hydrogen-bonding interaction with the Tyr542 at the voltage sensor domain. NS1643 has bound to Glu544 and Phe557 via forming polar and ␲-␲ stacking interactions, respectively. It seems that water-mediated interactions have given an enhanced role in binding of the Glu544 to the pair hydroxyl group of the ligand. The diagram expressed that Glu544 is highly pronounced and contributed in binding to the hydroxyl groups and as well as amide groups of the derivative E2. In addition, Tyr542 not only has bound as ␲-␲ stacking interactions to the ring of the compound but also established water-mediated bridges, as shown in the profile. Moreover, the other amino acids, which have participated in stabilization of the ligands, together with their occupancy interactions are determined in the diagrams.

Four different pharmacophore models (three 4-sited and one 3-sited models) were generated using docking poses at the K525 binding site: ADRR, ADDR, AARR and ADR, respectively. These models were used to screen Otava Tangible Ligand Library (∼2.6 million compounds). Best-fitted 1000 compounds were selected from screening and these ligands were docked into K525 binding site by using Glide/SP method. Docking poses were plotted respect to the fitness score, (Fig. 8). For each pharmacophore model, 3 ligands with highest fitness values and their corresponding docking scores are provided in Table 4. Same protocol was also applied for the E544 and N629 binding sites, successfully. Results are given in Tables 5 and 6, respectively. 3 high-scored E-pharmacophore models were obtained throughout the structurebased pharmacophore search analyses at the E544 site: DDD, ADDRR, and DDD (Fig. 9). The generated E-pharmacophores used in docking simulations at the corresponding binding sites. Representative NS1643/hERG1 complex pharmacophore models were ADDR and DDDHR for the N629 binding site (Fig. 10). All generated E-pharmacophore models have similar scaffolds, but they differ in their fitness values and docking scores. The representative structure (hERG1/NS1643 complex) from MD simulations was also used to generate pharmacophore model. 6-sited AADDDR model was found as the top-scored hypothesis and it is used in screening of ligand library. 2D structures of hit molecules with highest fitness values and their corresponding docking scores are collected in Table 7. Fig. 11 shows fitness values versus docking scores of selected top 1000 molecules with high fitness scores.

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

167

Table 7 Top fitted compounds and their docking scores respect to pharmacophore generated by docking poses (structures were taken from MD simulation) at the E544 site. Hits

2D Structure of Hits

Code

Fitness

Docking Score (kcal/mol)

1

7263314

2.556

−6.743

2

8882447

2.528

−7.472

3

7263473

2.528

−7.355

4

7263932

2.523

−7.391

5

5438046

2.515

−7.053

6

5438191

2.513

−7.261

7

5438184

2.498

−7.317

8

8882449

2.497

−7.234

9

5438303

2.495

−7.584

10

8882302

2.494

−7.722

3.6. Count of residue-ligand interactions A docking study is performed from E544 site for 21 hERG1 activators used in generation and validation studies of pharmacophore modeling and around 2000 docking poses were generated. Interaction Fingerprints script from Maestro is used to count residue-ligand interactions from backbone and from side chain atoms (Fig. 12). Results show that while Y542 is the most crucial residue from backbone interactions with these 21 hERG1 activators, E544 is the most essential amino acid for the side chain interactions for these ligands. Similarity calculation on docking poses were performed with Tanimato similarity metric using default settings. Clustering is then performed on the docking poses. Dendrogram representation of clustering and 3D and 2D ligand interaction diagrams for representative pose of NS1643 at the target binding site from clustering analyses of docking poses are shown in Fig. 13.

3.7. Limitations of the current models and avenues for the future development It is important to note that hERG activation mechanisms are very complex and often involve multiple binding sites and drug impact could not be reduced to one unique activator site. One may need to employ whole-cell kinetic modelling to de-convolute mechanisms of drug action and to isolate functional states impacted by binding of an activator [19,37] prior to structure-driven drug design. As we mention in the introduction many of the blockers are exhibiting multi-channel binding modes. Therefore, the viable strategy for developing safe hERG activators [38] has to include multi-channel screening and this would be a next step in the developing models of drug-induced cardiotoxicity [18,39]. An important limitation in the receptor-based model is quality of the hERG1 structures used for its development. Our study utilizes only open-pore model of

168

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

Fig. 12. Count of Residue and Ligand Interactions for backbone (top) and sidechain (bottom) atoms using derived 2000 docking poses of used 21 hERG1 activator. X and T represents chain name at the target protein.

hERG1 channel, which was developed using ROSETTA modelling and a multi-template approach with constraints derived from 2 trans-membrane (TM) bacterial K+ channels and 6 TM channels from Shakers family available at a time. Recently, the full channel structure of the highly homologous Eag1 channel has been solved by MacKinnon’s lab using Cryo-EM at ∼3.8 Å resolution [40]. While the voltage-sensor region is in the conformational corresponding to the open state of the channel, the pore domain appears to capture closed conformation with an intra-cellular gate closed. It is interesting to note that the pore domain models used in our study (S6 helix forming intracellular cavity, pore helix and selectivity filter regions) show an agreement to the published structure in positioning key residues for drug binding (T623, S624, Y652 and F656) [41]. The region that differs the most between models and solved structure is highly mobile S5-pore linker, unique for this family of proteins, but not involved in screening studies reported in this submission. Rapid developments of Cryo-EM technology is likely to drive fur-

ther refinement of hERG1 states and will have a major impact on the hERG1 channel activators design. 4. Conclusions Predicting the binding sites of hERG1 channel activators, analyses of their interactions with channel and biophysical mechanisms using computational techniques is a very challenging problem [9,17]. Even novel anti-arrhythmic drugs entering the market are often plagued by apparent cardio-toxicity [39]. The integration of data from various experimental and computational techniques is required to generate viable leads and to refine mechanisms of drug action [42]. In this study, ligand-based and structure-based modeling approaches have been combined with molecular docking and MD simulations to provide a new picture for the hERG1 channel activators field [5]. For this aim, a five-sited universal pharmacophore model (AAHRR) is successfully generated and validated

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170

169

Fig. 13. Dendogram representation of clustering of docking poses of 21 hERG1 openers and 3D and 2D ligand interaction diagrams for representative docking pose conformer of NS1643 at the target binding site.

with true unknowns for hERG1 channel openers. This model is combined with previously developed receptor-based hERG1 K channel model in the open state and novel hERG1 channel activators are proposed. The developed models may serve as basis for the synthesis of novel potential therapeutic hERG1 activators allowing for rapid isolation of most promising molecular leads.

Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jmgm.2017.03. 020.

References [1] B. Hille, Modulation of ion-channel function by G-protein-coupled receptors, Trends Neurosci. 17 (12) (2016) 531–536. [2] M. Sheng, D.T.S. Pak, Ligand-gated ion channel interactions with cytoskeletal and signaling proteins, Annu. Rev. Physiol. 62 (1) (2000) 755–778. [3] M.R. Blatt, G. Thiel, Hormonal control of ion channel gating, Annu. Rev. Plant Physiol. Plant Mol. Biol. 44 (1) (1993) 543–567. [4] B. Nilius, G. Droogmans, Ion channels and their functional role in vascular endothelium, Physiol. Rev. 81 (4) (2001) 1415–1459. [5] E. Grandi, M.C. Sanguinetti, D.C. Bartos, D.M. Bers, Y. Chen-Izu, N. Chiamvimonvat, H.M. Colecraft, B.P. Delisle, J. Heijman, M.F. Navedo, S. Noskov, C. Proenza, J.I. Vandenberg, V. Yarov-Yarovoy, Potassium channels in the heart: structure, function and regulation, J. Physiol. 595 (2017) 2209–2228. [6] J. Warmke, R. Drysdale, B. Ganetzky, A distinct potassium channel polypeptide encoded by the Drosophila eag locus, Science 252 (5012) (1991) 1560–1562. [7] D.J. Snyders, Structure and function of cardiac potassium channels, Cardiovasc. Res. 42 (2) (1999) 377–390.

[8] W.A. Volberg, B.J. Koci, W. Su, J. Lin, J. Zhou, Blockade of human cardiac potassium channel human ether-a-go-go-related gene (HERG) by macrolide antibiotics, J. Pharmacol. Exp. Ther. 302 (1) (2002) 320–327. [9] M.C. Sanguinetti, J.S. Mitcheson, Predicting drug–hERG channel interactions that cause acquired long QT syndrome, Trends Pharmacol. Sci. 26 (3) (2005) 119–124. [10] I. Splawski, J. Shen, K.W. Timothy, M.H. Lehmann, S. Priori, J.L. Robinson, A.J. Moss, P.J. Schwartz, J.A. Towbin, G.M. Vincent, M.T. Keating, Spectrum of mutations in long-QT syndrome genes. KVLQT1, HERG, SCN5A, KCNE1, and KCNE2, 102 (10) (2000) 1178–1185. [11] J.I. Vandenberg, B.D. Walker, T.J. Campbell, HERG K+ channels: friend and foe, Trends Pharmacol. Sci. 22 (5) (2001) 240–246. [12] D. Thomas, C.A. Karle, J. Kiehn, The cardiac hERG/IKr potassium channel as pharmacological target: structure, function, regulation, and clinical applications, Curr. Pharm. Des. 12 (18) (2006) 2271–2283. [13] H. Zeng, I.M. Lozinskaya, Z. Lin, R.N. Willette, D.P. Brooks, X. Xu, Mallotoxin is a novel human ether-a-go-go-related gene (hERG) potassium channel activator, J. Pharmacol. Exp. Ther. 319 (2) (2006) 957–962. [14] J. Kang, X.-L. Chen, H. Wang, J. Ji, H. Cheng, J. Incardona, W. Reynolds, F. Viviani, M. Tabart, D. Rampe, Discovery of a small molecule activator of the human ether-a-go-go-related gene (HERG) cardiac K+ Channel, Mol. Pharmacol. 67 (3) (2005) 827–836. [15] O. Casis, S.-P. Olesen, M.C. Sanguinetti, Mechanism of action of a novel human ether-a-go-go-related gene channel activator, Mol. Pharmacol. 69 (2) (2006) 658–665. [16] A. Anwar-Mohamed, K.H. Barakat, R. Bhat, S.Y. Noskov, D.L. Tyrrell, J.A. Tuszynski, M. Houghton, A human ether-a-go-go-related (hERG) ion channel atomistic model generated by long supercomputer molecular dynamics simulations and its use in predicting drug cardiotoxicity, Toxicol. Lett. 230 (3) (2014) 382–392. [17] S.D. Longman, T.C. Hamilton, Potassium channel activator drugs: mechanism of action, pharmacological properties, and therapeutic potential, Med. Res. Rev. 12 (2) (1992) 73–148. [18] P.C. Yang, J.D. Moreno, C.Y. Miyake, S.B. Vaughn-Behrens, M.T. Jeng, E. Grandi, X.H. Wehrens, S.Y. Noskov, C.E. Clancy, In silico prediction of drug therapy in catecholaminergic polymorphic ventricular tachycardia, J. Physiol. 594 (3) (2016) 567–593. [19] J. Guo, Y.M. Cheng, J.P. Lees-Miller, L.L. Perissinotti, T.W. Claydon, C.M. Hull, S. Thouta, D.E. Roach, S. Durdagi, S.Y. Noskov, H.J. Duff, NS1643 interacts around

170

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28] [29] [30]

[31]

[32]

[33] [34]

[35] [36]

[37]

S. Durdagi et al. / Journal of Molecular Graphics and Modelling 74 (2017) 153–170 L529 of hERG to alter voltage sensor movement on the path to activation, Biophys. J. 108 (6) (2015) 1400–1413. P.-Z. Zhou, J. Babcock, L.-Q. Liu, M. Li, Z.-B. Gao, Activation of human ether-a-go-go related gene (hERG) potassium channels by small molecules, Acta Pharmacol. Sin. 32 (6) (2011) 781–788. M.J. Perrin, R.N. Subbiah, J.I. Vandenberg, A.P. Hill, Human ether-a-go-go related gene (hERG) K+ channels: function and dysfunction, Prog. Biophys. Mol. Biol. 98 (2–3) (2008) 137–148. R. Brugada, K. Hong, R. Dumaine, J. Cordeiro, F. Gaita, M. Borggrefe, T.M. Menendez, J. Brugada, G.D. Pollevick, C. Wolpert, E. Burashnikov, K. Matsuo, Y. Sheng Wu, A. Guerchicoff, F. Bianchi, C. Giustetto, R. Schimpf, P. Brugada, C. Antzelevitch, Sudden death associated with short-QT syndrome linked to mutations in HERG, Circulation 109 (1) (2004) 30–35. Z. Su, J. Limberis, A. Souers, P. Kym, A. Mikhail, K. Houseman, G. Diaz, X. Liu, R.L. Martin, B.F. Cox, G.A. Gintant, Electrophysiologic characterization of a novel hERG channel activator, Biochem. Pharmacol. 77 (8) (2009) 1383–1390. X. Xu, M. Recanatini, M. Roberti, G.-N. Tseng, Probing the binding sites and mechanisms of action of two human ether-a-go-go-related gene channel activators, 1,3-bis-(2-hydroxy-5-trifluoromethyl-phenyl)-urea (NS1643) and 2-[2-(3,4-dichloro-phenyl)-2,3-dihydro-1H-isoindol-5-ylamino]-nicotinic acid (PD307243), Mol. Pharmacol. 73 (6) (2008) 1709–1721. E. Gordon, I.M. Lozinskaya, Z. Lin, S.F. Semus, F.E. Blaney, R.N. Willette, X. Xu, 2-[2-(3,4-Dichloro-phenyl)-2,3-dihydro-1H-isoindol-5-ylamino]-nicotinic Acid (PD-307243) causes instantaneous current through human ether-a-go-go-related gene (hERG) potassium channels, Mol. Pharmacol. 73 (3) (2007) 638–651. M. Perry, F.B. Sachse, M.C. Sanguinetti, Structural basis of action for a human ether-a-go-go-related gene 1 potassium channel activator, Proc. Natl. Acad. Sci. 104 (34) (2007) 13827–13832. S. Durdagi, J. Guo, J.P. Lees-Miller, S.Y. Noskov, H.J. Duff, Structure-guided topographic mapping and mutagenesis to elucidate binding sites for the human ether-a-go-go-related gene 1 potassium channel (KCNH2) activator NS1643, J. Pharmacol. Exp. Ther. 342 (2) (2012) 441–452. LigPrep, Schrödinger, LLC, New York, NY, 2009. Maestro, Schrödinger, LLC, New York, NY, 2009. S.L. Dixon, A.M. Smondyrev, E.H. Knoll, S.N. Rao, D.E. Shaw, R.A. Friesner, PHASE: a new engine for pharmacophore perception 3D QSAR model development, and 3D database screening: 1. Methodology and preliminary results, J. Comput. Aided Mol. Des. 20 (10) (2006) 647–671. R.A. Friesner, J.L. Banks, R.B. Murphy, T.A. Halgren, J.J. Klicic, D.T. Mainz, M.P. Repasky, E.H. Knoll, M. Shelley, J.K. Perry, D.E. Shaw, P. Francis, P.S. Shenkin, Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy, J. Med. Chem. 47 (7) (2004) 1739–1749. R.A. Friesner, R.B. Murphy, M.P. Repasky, L.L. Frye, J.R. Greenwood, T.A. Halgren, P.C. Sanschagrin, D.T. Mainz, Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-Ligand complexes, J. Med. Chem. 49 (21) (2006) 6177–6196. W. Sherman, H.S. Beard, R. Farid, Use of an induced fit receptor structure in virtual screening, Chem. Biol. Drug Des. 67 (1) (2006) 83–84. G. Jones, P. Willett, R.C. Glen, A.R. Leach, R. Taylor, Development and validation of a genetic algorithm for flexible docking1, J. Mol. Biol. 267 (3) (1997) 727–748. CombiGlide, Schrödinger, LLC, New York, NY, 2009. J.P. Lees-Miller, J.O. Subbotina, J. Guo, V. Yarov-Yarovoy, S.Y. Noskov, H.J. Duff, Interactions of H562 in the S5 helix with T618 and S621 in the pore helix are important determinants of hERG1 potassium channel structure and function, Biophys. J. 96 (9) (2009) 3600–3610. L.L. Perissinotti, J. Guo, P.M. De Biase, C.E. Clancy, H.J. Duff, S.Y. Noskov, Kinetic model for NS1643 drug activation of WT and L529I variants of Kv11. 1 (hERG1) potassium channel, Biophys. J. 108 (6) (2015) 1414–1424.

[38] J. Guo, S. Durdagi, M. Changalov, L.L. Perissinotti, J.M. Hargreaves, T.G. Back, S.Y. Noskov, H.J. Duff, Structure driven design of novel human ether-a-go-go-related-gene channel (hERG1) activators, PLoS One 9 (9) (2014) e105553. [39] J.P. Lees-Miller, J. Guo, Y. Wang, L.L. Perissinotti, S.Y. Noskov, H.J. Duff, Ivabradine prolongs phase 3 of cardiac repolarization and blocks the hERG1 (KCNH2) current over a concentration-range overlapping with that required to block HCN4, J. Mol. Cell. Cardiol. 85 (2015) 71–78. [40] J.R. Whicher, R. MacKinnon, Structure of the voltage-gated K(+) channel Eag1 reveals an alternative voltage sensing mechanism, Science 353 (6300) (2016) 664–669. [41] Y. Wang, J. Guo, L.L. Perissinotti, J. Lees-Miller, G. Teng, S. Durdagi, H.J. Duff, S.Y. Noskov, Role of the pH in state-dependent blockade of hERG currents, Sci. Rep. 6 (2016) 32536. [42] Y. Zhang, V.H. Barocas, S.A. Berceli, C.E. Clancy, D.M. Eckmann, M. Garbey, G.S. Kassab, D.R. Lochner, A.D. McCulloch, R. Tran-Son-Tay, N.A. Trayanova, Multi-scale modeling of the cardiovascular system: disease development, progression, and clinical intervention, Ann. Biomed. Eng. 44 (9) (2016) 2642–2660. [43] Z. Yu, J.P.D. van Veldhoven, I. ME’t Hart, A.H. Kopf, L.H. Heitman, A.P. Ijzerman, Synthesis and biological evaluation of negative allosteric modulators of the Kv11.1(hERG) channel, Eur. J. Med. Chem. 106 (2015) 50–59. [44] Z. Yu, J. Liu, J.P.D. van Veldhoven, A.P. IJzerman, M.J. Schalij, D.A. Pijnappels, L. Heitman, A.A.F. H.; de Vries, Allosteric modulation of Kv 11.1 (hERG) channels protects against drug-Induced ventricular arrhythmias, Circ. Arrhythm. Electrophysiol. 9 (4) (2016). [45] R.S. Hansen, T.G. Diness, T. Christ, J. Demnitz, U. Ravens, S.-P. Olesen, M. Grunnet, Activation of human ether-a-go-go-related gene potassium channels by the diphenylurea 1,3-bis-(2-hydroxy-5-trifluoromethyl-phenyl)-urea (NS1643), Mol. Pharmacol. 69 (1) (2006) 266–277. [46] R. Mannikko, M.H. Bridgland-Taylor, H. Pye, S. Swallow, N. Abi-Gerges, M.J. Morton, C.E. Pollard, Pharmacological and electrophysiological characterization of AZSMO-23, an activator of the hERG K+ channel, Br. J. Pharmacol. 172 (12) (2015) 3112–3125. [47] H. Zhang, B. Zou, H. Yu, A. Moretti, X. Wang, W. Yan, J.J. Babcock, M. Bellin, O.B. McManus, G. Tomaselli, F. Nan, K.-L. Laugwitz, M. Li, Modulation of hERG potassium channel gating normalizes action potential duration prolonged by dysfunctional KCNQ1 potassium channel, Proc. Natl. Acad. Sci. 109 (29) (2012) 11866–11871. [48] M. Perry, F.B. Sachse, J. Abbruzzese, M.C. Sanguinetti, PD-118057 contacts the pore helix of hERG1 channels to attenuate inactivation and enhance K+ conductance, Proc. Natl. Acad. Sci. 106 (47) (2009) 20075–20080. [49] S. Durdagi, J. Subbotina, J. Lees-Miller, J. Guo, H.J. Duff, S.Y. Noskov, Insights into the molecular mechanism of hERG1 channel activation and blockade by drugs, Curr. Med. Chem. 17 (30) (2010) 3514–3532. [50] S. Durdagi, C. Zhao, J.E. Cuervo, S.Y. Noskov, Atomistic models for free energy evaluation of drug binding to membrane proteins, Curr. Med. Chem. 18 (17) (2011) 2601–2611. [51] S. Durdagi, H.J. Duff, S.Y. Noskov, Combined Receptor and Ligand-Based Approach to the Universal Pharmacophore Model Development for Studies of Drug Blockade to the hERG1 Pore Domain, J. Chem. Inf. Model 51 (2) (2011) 463–474. [52] S. Durdagi, T. Randall, H.J. Duff, A. Chamberlin, S.Y. Noskov, Rehabilitating drug-induced long-QT promoters: In-silico design of hERG-neutral cisapride analogues with retained pharmacological activity, BMC Pharmacol. Toxicol. 15 (14) (2014) 1–15. [53] S. Durdagi, S. Deshpande, H.J. Duff, S.Y. Noskov, Modeling of Open, Closed, and Open-Inactivated States of the hERG1 Channel: Structural Mechanisms of the State-Dependent Drug Binding, J. Chem. Inf. Model 52 (10) (2012) 2760–2774.