Visual Overview
Abstract
Abnormal cardiac electrical activity is a common side effect caused by unintended block of the promiscuous drug target human ether-à-go-go-related gene (hERG1), the pore-forming domain of the delayed rectifier K+ channel in the heart. hERG1 block leads to a prolongation of the QT interval, a phase of the cardiac cycle that underlies myocyte repolarization detectable on the electrocardiogram. Even newly released drugs such as heart-rate lowering agent ivabradine block the rapid delayed rectifier current IKr, prolong action potential duration, and induce potentially lethal arrhythmia known as torsades de pointes. In this study, we describe a critical drug-binding pocket located at the lateral pore surface facing the cellular membrane. Mutations of the conserved M651 residue alter ivabradine-induced block but not by the common hERG1 blocker dofetilide. As revealed by molecular dynamics simulations, binding of ivabradine to a lipophilic pore access site is coupled to a state-dependent reorientation of aromatic residues F557 and F656 in the S5 and S6 helices. We show that the M651 mutation impedes state-dependent dynamics of F557 and F656 aromatic cassettes at the protein-lipid interface, which has a potential to disrupt drug-induced block of the channel. This fundamentally new mechanism coupling the channel dynamics and small-molecule access from the membrane into the hERG1 intracavitary site provides a simple rationale for the well established state-dependence of drug blockade.
SIGNIFICANCE STATEMENT The drug interference with the function of the cardiac hERG channels represents one of the major sources of drug-induced heart disturbances. We found a novel and a critical drug-binding pocket adjacent to a lipid-facing surface of the hERG1 channel, which furthers our molecular understanding of drug-induced QT syndrome.
Introduction
The cardiac action potential is primarily generated by sodium and calcium channels, which depolarize the membrane potential, and by potassium channels, which repolarize the membrane potential and terminate the action potential (Nerbonne and Kass, 2005). Gene mutations, age-related factors, and drug-induced toxicity are all linked to various perturbation of action potentials, leading to potential lethal disorders of heart rhythm (arrhythmias) (Chiamvimonvat et al., 2017). Several K+-selective channels were identified as major determinants of proarrhythmic activity and main targets in antiarrhythmic drug development. Perhaps the most impactful of these is the human ether-a-go-go-related gene 1 (hERG1 or Kv11.1) channel, the K+-selective channel carrying the rapid delayed rectifier current (IKr) in myocytes (Trudeau et al., 1995). The physiologic role of IKr is to repolarize the late phase of cardiac action potential; hence, currents carried by hERG1 that contain human mutations are linked to arrhythmias (Gustina and Trudeau, 2009; Vandenberg et al., 2012). Likewise, pharmacological blockage of the IKr can lead to a prolongation of the QT interval (phase 3 repolarization), causing a drug-induced long-QT syndrome (Compton et al., 1996; Roden et al., 1996; Splawski et al., 1997; Huang et al., 2001; Vandenberg et al., 2001). Drug-induced block of hERG1 and its associated prolongation of the QT interval and proarrhythmia has resulted in hERG1 being one of the most studied ion channels (Numaguchi et al., 2000; Witchel, 2011; Sanguinetti, 2014). Recent high-throughput screening studies have provided additional evidence for the central role of hERG1 blockade in drug-safety assessments (Di Veroli et al., 2013a).
Despite the plethora of experimental data available for drug interactions with the hERG1 channel, a-priori prediction of the cardiotoxic potential of a novel compound during a preclinical developmental stage is a complex and a challenging task. The torsadogenicity is an emergent and a complicated property that depends on several factors: the conformational state of the channel being targeted by the compound (Chen et al., 2002; Stork et al., 2007; Lees-Miller et al., 2015; Wu et al., 2015), general kinetics of the drug access mechanisms (Guo et al., 2005; Di Veroli et al., 2013b; Hill et al., 2014), the solubility of the drug in the membrane (Wacker and Noskov, 2018), patient-to-patient variance in isoform composition of hERG1 (hERG1a/hERG1b) (Sale et al., 2008), and isoform-specific interactions with blockers (Abi-Gerges et al., 2011). Given these factors, understanding molecular determinants of drug-induced QT-prolongation remains one of the fundamental and elusive problems in the field of molecular pharmacology. Even newly released “cardio-safe” drugs entering the market, such as ivabradine, are found to block the IKr, prolong action potential duration, and can contribute to a potentially lethal drug-induced arrhythmia known as torsades de pointes (Duff et al., 1995; Lees-Miller et al., 2000, 2015; Chen et al., 2002; Melgari et al., 2015).
Ivabradine is commonly prescribed as a blocker of If (funny current) and acts as a heart-rate lowering agent for symptomatic management of the chronic heart failure. However, over a similar range of concentrations, the drug also blocks hERG; thus, ivabradine prolongs phase 3 of the action potential and has been reported to induce torsades de pointes when applied in a poly-pharmaceutical context (Hancox et al., 2015; Melgari et al., 2015; Frommeyer et al., 2017). In a previous study, by combining electrophysiology and molecular modeling techniques, we emphasized the importance of the lipophilic interface and high-affinity state-dependent blockade of hERG1 by ivabradine (Lees-Miller et al., 2015). The molecular simulations performed on the homology model of hERG1 pore domain in open and closed states indicated that ivabradine may bind to a lipid-facing binding pocket centered at the M651 residue (Lees-Miller et al., 2015). However at the time, findings were limited owing to the lack of structural information on the organization of the pore domain of the hERG1 channel. In 2017, a high-resolution cryoelectron microscopy (cryo-EM) structure (3.4 Å) of a highly homologous neuronal human ether à go-go 1 gene (hEAG1) channel was resolved with a voltage-sensor in a depolarized (open) state, but with the pore domain closed. Later in the same year, the high-resolution structure of the hERG1 channel (3.2 Å) was reported in the open state (Whicher and MacKinnon, 2016; Wang and MacKinnon, 2017) (Fig. 1). These new structural atomic details provided the opportunity to unravel the potential lipophilic access mechanisms for the ivabradine-induced block of hERG1, to assess the impact of mutations at the M651 site in the distal S6, and to understand the dynamics of allosteric coupling between residues involved in the high-affinity channel block.
Herein, we report an in-depth study that provides direct evidence of ivabradine’s interaction in a state-dependent manner with lipid-facing residues as part of the hERG1 blockade process. We assessed whether mutations of the lipid-facing residues, when coupled with C-type inactivation impacting mutations, affect the concentration of ivabradine required to block the hERG1 current. To reveal mechanisms involved in coupling dynamics of the lipid-facing residues with respect to the aromatic cassette involved in high-affinity drug block by various compounds (Fig. 1), we performed several multi-microsecond molecular dynamics (MD) simulations of the wild-type (WT) and mutated forms of hERG1 in conjunction with molecular biology and electrophysiological studies. The MD simulations in combination with the ensemble docking simulations allowed us to map molecular details of the probable lipophilic access pathway of ivabradine and potentially explain the drug’s dependence on the inactivation process of hERG1 channel. Our results show that ivabradine’s binding at the vicinity of the lipid-facing residue M651 is directly coupled to the conformational dynamics (reorientation) of the aromatic cassettes (F656 and F557) in the S6 (Fig. 1, A and B). Disruption of this allosteric coupling between drug binding on the lipid-facing surface and conformational dynamics of F656/F557 was shown to greatly attenuate ivabradine blockade.
Materials and Methods
Molecular Biology.
Methods for site-directed mutagenesis have been previously reported (Lees-Miller et al., 2015; Wang et al., 2016). The hERG1 constructs were transfected into mammalian human embryonic kidney cells. Single- and double-mutant constructs of hERG1 were produced using conventional overlap polymerase chain reaction with primers synthesized by Sigma Genosys (Oakville, ON, Canada) and sequenced using Eurofins MWG Operon (Huntsville, AL). Constructs were linearized with XbaI restriction endonuclease, and cRNA was transcribed in vitro using the mMessage mMachine T7 Ultra cRNA transcription kit (Ambion, Austin, TX).
General Setup for Electrophysiological Recordings.
The extracellular solution contained (in millimolars) NaCl 140, KCl 5.4, CaCl2 1, MgCl2 1, HEPES 5, and glucose 5.5; the pH of the solution was adjusted and kept at 7.4 with NaOH. Micropipettes were pulled from borosilicate glass capillary tubes on a programmable horizontal puller (Sutter Instruments, Novato, CA). The pipette solution contained the following: 10 mM KCl, 110 mM K-aspartate, 5 mM MgCl2, 5 mM Na2ATP, 10 mM EGTA (ethylene glycol-bis(-aminoethyl ether)-N,N,N,N tetraacetic acid), 5 mM HEPES, and 1 mM CaCl2; the solution was adjusted to pH 7.2 with KOH. Standard patch-clamp methods were used to measure the whole-cell currents of hERG1 mutants expressed in human embryonic kidney 293 cells using the AXOPATCH 200B amplifier (Axon Instruments). Unless otherwise indicated, the tail currents were recorded when the voltage was returned to −100 mV from +50 mV. Transfected human embryonic kidney cells were patched to record the hERG1 currents. Ivabradine was directly dissolved into the Tyrode solution right before the experiments. The solutions were used for the next 2 hours during the experiments. The stock solution of 100 μM ivabradine was prepared in the extracellular solution. Fresh stock solutions of ivabradine were prepared weekly.
Voltage-Dependence of Activation.
From a holding potential of −80 mV cells were depolarized for 1 second to a range of voltages from –100 to +40 mV followed by a step to −100 mV (1 second) to record the tail currents. The isochronal tail-current–voltage plots were fit to a single Boltzmann function (1):(1)where I/Imax is the normalized current, V1/2 is voltage of the half-maximal activation, k is the slope factor, and Vm is the membrane potential.
Analysis of Deactivation.
Deactivation of hERG1 tail currents was measured by activating channels at +40 mV, followed with a short (5-millisecond) repolarization step to −120 mV and deactivating steps at −120, −100, −60, and −40 mV. Currents at different voltages were normalized and fitted. The fitted data were averaged (n = 10).
Statistical Analysis of Electrophysiological Experiments.
Statsview (Abacus Concepts, Berkeley, CA) or QTI plot (Vasilef, 2013), Grace (http://plasma-gate.weizmann.ac.il/Grace/) were used to analyze the data. The null hypothesis of this study predicted no difference between the IC50 values of the single mutation to the double mutations assessed. The null hypothesis was rejected when the P value was <0.05 as evaluated by a one-way analysis of variance with Tukey test. The exact P values for data in Fig. 4 are presented in the Supplemental Table 1. All variance measures (bars) for electrophysiological data are shown as S.D. The study was exploratory, thus there was no a priori reason to consider whether there was an additive or subtractive interaction. In addition, we acknowledge that there may be other mutations, unknown to us at this time, which could be relevant. All of the comparisons were prespecified and all of the comparisons are reported. A priori, we generally required a minimum of n = 2 independent experiments for each point on the IC50 curve. However, the n values near the IC50 point and near the maximum blockade have a major impact on the reliability of the IC50 estimate. Frequently we increased the n values at these putative points to be more certain about the reliability of the measurement. The n values for each point of the concentration-response relationship are presented in the figure legends. J.G. executed the experiments and analyzed the experimental data, hence he was not blinded. However, each experiment was reviewed by a small committee of individuals to assess quality of the records at weekly laboratory meetings (J.G. and H.J.D.) and on a monthly basis (H.J.D. and S.Y.N.) to review raw experimental data.
Dofetilide and Ivabradine Water/Hexane Partitioning Experiments.
The distribution coefficient for dofetilide and ivabradine were determined by the use of a classic shake-flask method detailed in our previous publications (Perlovichl and Bauer-Brandl, 2003; Perlovich et al., 2006; Blokhina et al., 2016). Prior to conducting the experiments, both solvents were mutually saturated to reach equilibrium by slow stirring into a biphasic system for 2 days. Ivabradine and dofetilide were dissolved in a buffer at pH=7.4 and then were added to a hexane solution. To ensure complete equilibration of the system as indicated by the absence of turbidity on each phase, the flasks were shaken for about 48 hours in a thermostatic water bath at 293.15, 298.15, 303.15, 310.15, and 313.15 K. After reaching equilibration in the system, the samples from the lower phase were carefully removed with syringes for analysis. The molar concentrations of the dofetilide and ivabradine in the buffer phase were measured by a Cary-50 spectrophotometer (Varian) with an accuracy of 2%–4%. The experimental results are reported as an average value of at least three replicated experiments. The accuracy of the distribution coefficients were verified by comparing the starting mass of a compound and the total mass of the compound distributed in the two phases. The reproducibility of the measured concentrations was under 0.1%, and the maximum deviations from the average value were <0.15%. The ratio of the compound equilibrium concentration in mole fraction in the hexane phase (xH) to those in the aqueous phase (xB) was determined using the hexane/buffer distribution coefficient in the following form:(2)DH/B is the phase equilibrium constant for a drug distributed in the hexane phase and the saturated buffer phase.
The standard Gibbs energy of transfer ∆trGo from the buffer into an organic solvent was calculated by using:
(3)The temperature dependence of distribution (van’t Hoff method) was employed to obtain the enthalpy of transfer ∆trHo:
(4)The entropy of transfer ∆trSo was calculated from:
(5)In-depth details on the partitioning experiments and quality controls are provided in the Supplemental Materials.
Molecular Dynamics Simulations.
The three-dimensional structure for the open state of the channel in this study is the recently published high-resolution cryo-EM structure of the transmembrane domain of hERG1channel (PDB ID 5VA2). The structure was truncated before the Per-Arnt-Sim domain (PAS) and after cyclic nucleotide-binding domain domains as described in our previous publication (Wacker et al., 2017; Perissinotti et al., 2018). The basis for the three-dimensional structure for the closed-state hERG1 channel is the homology modeling of hEAG1 with template derived from the cryo-EM structure (PDB ID 5K7L) solved at 3.78 Å resolution. The SWISS-MODEL program (Kopp and Schwede, 2004) was used to develop the hERG1 closed homology model from the available hEAG1 channel structure. Sequence alignment was performed using the CLUSTALW algorithm (Thompson et al., 1994). The sequence similarity between hERG1 and hEAG1 channels for the pore domain (S5–S6) is over 75% (Wacker et al., 2017). The detailed analysis of two structures was published recently (Vandenberg et al., 2017; Wacker et al., 2017). The main differences between pore domains of the two channels are located in the extended turret region connecting S5 to the pore helix. The following three-step protocol was adapted to model missing residues and flexible elements: 1) threading for generation of initial models on the basis of template structure by copying coordinates over the aligned regions (for closed states), 2) low-resolution ROSETTA loop modeling using the cyclic coordinate descent method, 3) high-resolution all-atom refinement and selection of models on the basis of ROSETTA clustering (Bender et al., 2016).
Models of the protein were generated from the alignment in a stepwise manner. CHARMM-GUI (Jo et al., 2008) was used to prepare protein—dipalmitoylphosphatidyl choline lipid bilayer complexes solvated in 150 mM KCl aqueous solution using CHARMM-36 force-field and TIP3P water model (Jorgensen et al., 1983; MacKerell et al., 1998; Noskov et al., 2004; Noskov and Roux, 2008; Klauda et al., 2010; Best et al., 2012). The fully assembled systems were equilibrated for 10 nanoseconds using NAMD2.10 (Phillips et al., 2005) and then subjected to production runs with the Anton 2 supercomputer. The production runs were performed for 1.0–2.5 microseconds, each with CHARMM36M (Huang et al., 2017) force-field to assess structural dynamics of residues involved in lipophilic access in hERG1 WT and selected mutants. The production runs were executing in a semi-isotropic (NPaT) ensemble at a temperature of 315 K maintained by the Nosé-Hoover thermostat (Martyna et al., 1994). The time-step for production runs was set to 2 femtoseconds, and trajectories were saved every 240 picoseconds. Nonbonded and long-range electrostatic interactions were evaluated every 2 and 6 femtoseconds, respectively. Long-range electrostatics was calculated using the k-Gaussian–Ewald method implemented to enhance performance on Anton 2 platform (Shan et al., 2005; Shaw et al., 2014) with a 64 Å × 64 Å × 64 Å grid. SHAKE was used to constrain all bonds involving hydrogen atoms. All of the subsequent trajectory analysis was performed using the CHARMM c40b2 program package (Brooks et al., 2009).
Molecular Docking Protocols.
Ivabradine in neutral and charged states was docked in-silico to the hERG1 open- and closed-pore structures to establish binding modes. This was performed via the Induced-Fit Docking (IFD) approach available in the Schrödinger Suite (Schrodinger LLC, 2016). Schrödinger’s IFD protocol uses Glide and Prime to exhaustively consider possible binding modes and the associated conformational changes within the receptor’s active sites (Schrodinger LLC, 2016). In addition to the well established high-affinity binding pocket in the intracellular cavity of hERG1 channels, previous blinded docking studies revealed several alternative binding regions present in the hERG1 homology models and the cryo-EM-derived structures (Lees-Miller et al., 2015; Saxena et al., 2016; Wacker et al., 2017). To map these alternative binding pockets as accurately as possible, we adapted the following two-step protocol. First, a blind docking was performed covering an entire receptor (e.g., pore and voltage-sensing domains regions of hERG1), and then elucidated common binding sites were selected for further studies with high-precision grid mapping. The SiteMap module of the Schrödinger molecular modeling package was used to generate fine-grids for the subsequent precision docking. Following previously tested protocols, the grid was defined to 15 “site-points” for each pocket found during the blinded docking run (Schrodinger LLC, 2016). Then, site maps are cropped 10 Å from the nearest site-point (Halgren, 2009). The scoring was carried out using Schrödinger’s discretized version of the ChemScore empirical scoring function, resulting in a small number of best-refined poses, following which the best-docked protein-ligand complex was determined on the basis of a model energy score (Emodel) that combined the energy of the grid-score, the binding affinity predicted by GlideScore, and (for flexible docking) the internal strain energy for the model potential used to direct the conformational-search algorithm (Schrodinger LLC, 2016). The basis for the partial charges for the neutral form of the drug was the nonbonded parameters from the Optimized Potentials for Liquid Simulation Force Field (OPLS3) with parameters specifically optimized for drug-like molecules (Harder et al., 2016). For all docking simulations the funnel width was increased by adjusting the energy window to 5.0, the CvdW cutoff was set to 10.0 kcal/mol, and the clustering criteria was set to 0.75 and extra-precision (XP) (Friesner et al., 2006). In addition to the single-structure docking procedure, MD-generated ensembles for the open hERG1 transmembrane and mutant forms were used in an ensemble docking with ivabradine to account for the conformational dynamics of the pore domain. A similar protocol was used in a comparative study of a cationic ivabradine binding to the main intracavitary site present in the open state of hERG1.
Ensemble Docking Procedure.
From the last 1.0 microseconds of the production MD trajectories, we randomly selected 25 frames spaced every 40 nanoseconds. The basis of the alignment of each frame was the position of the backbone atoms from the pore domain (residues 545–572 and 635–669). Glide was used with the XP ensemble docking (Friesner et al., 2006) with Schrodinger Small-Molecule Drug Discovery Suite 2018-2 (2016) as described above for the single-structure docking protocol. The ligand binding site defined in a single-structure receptor docking was the basis for the generation of each receptor grid. Each generated grid is made of two boxes: the inner box for searching docking space that defines acceptable volume for the ligand center to explore and the outer box for searching the space of valid poses that must contain all ligand atoms for which grid potentials are computed. The inner cubic box was fixed at the dimension of 15 Å × 15 Å × 15 Å and the outer box was set to 30 Å × 30 Å × 30 Å. The number of selected docking poses per docking simulation was set to approx. 50, and only poses with binding affinities of < −3 kcal/mol were clustered. It is important to mention that grids generated by the SiteMap protocol were overlapping. Therefore, poses obtained from separate docking simulations were clustered into three sites on the basis of where the ligand was bound.
The first site corresponds to all poses found inside the intracellular cavity of the receptor (within approx. 3 Å of the Cα atoms of Y652 or F656); poses on the lipid-facing transmembrane domain (lipophilic, within approx. 3 Å of the Cα atom of M651) maps lipophilic site 2 and poses found between transmembrane segments of the pore domain (S5–S6) and the voltage sensing domain (S1–S4) map a putative “access path” site 3.
Results
Rationale for the Selection of Mutants.
Our previous in-silico screening of ivabradine binding to the hERG1 model representing open and closed states of the channel identified several potential binding modes (Lees-Miller et al., 2015). The best-scored binding poses for neutral and cationic ivabradine were clustered around a well known intracavitary site in the open state of hERG1; however, almost no binding in the internal cavity (Y652-F656) was observed for the closed-state model. Two aromatic residues in the S6 helix lining the intracavitary site that are most commonly associated with high-affinity blockade of hERG1 are Y652 and F656. Both residues are established as a critical determinant of hERG1-induced block associated with proarrhythmia (Ficker et al., 1998; Lees-Miller et al., 2000; Mitcheson et al., 2000; Perry et al., 2010). Recent work of Saxena et al. (2016) emphasized an important role for another aromatic residue (F557) from the S5 helix. It was shown that F557 may be involved in the drug stabilization of hERG1 binding pocket along with the residues in the S6 helix. Therefore, we chose to include F557X along with classic aromatic cassette Y652X and F656X to screen for ivabradine binding to the well established intracellular site.
As for the M651 mutation, previous blinded docking studies showed the presence of a potential binding domain for ivabradine in close proximity to this residue, only in the open state of the channel (Lees-Miller et al., 2015). However, no electrophysiological experimental data on M651X mutants was reported in the previous report. The MD simulations also showed a favorable energetic partitioning of the neutral state ivabradine into the lipid bilayers. The comparison of in-silico models used by Lees-Miller et al. (2015) to the recently solved cryo-EM structures showed that the key structural elements in the pore domain (S5–S6) were accurately captured by ROSETTA-generated models of hERG1 (RMSD < 3.5 Å) (Wang et al., 2016) providing additional support for docking studies on hERG1 models. In this study we created the M651T mutation to examine its effects on ivabradine binding (Fig. 2). We reasoned that the substitution of methionine by the small polar amino acid (threonine) at the M651 site could disrupt the lipophilic binding or access route of ivabradine. We discovered that M651T expressed well in cells but the mutation slowed deactivation kinetics of the channel and suppressed the block induced by ivabradine (Fig. 2). This raised a possibility that M651T modifies ivabradine-induced block by altering the deactivation time-course of the channel. To account for this covariable and to directly address whether the slowing of deactivation contributed to the shift in the concentration-response relationship to ivabradine, we created a double mutation, M651T/T618I. The rationale for its creation relates to a previous study wherein we reported that T618I accelerated deactivation kinetics. We empirically created the double mutation M651T/T618I to restore deactivation to WT values while retaining the key M651T substitution.
Mutations of Lipid-Facing Residue M651 Significantly Alters Thermodynamics and Kinetics of Ivabradine.
Figure 2, A and B, show raw data of the concentration-dependent ivabradine-induced block of WT hERG1 (Fig. 2A) in comparison with the M651T mutant (Fig. 2B). Figure 2C shows a representative time-course of block of IhERG1 by ivabradine at various concentrations. Iherg1 designates the Ikr-like current elicited by transfecting the HEK cells with the hERG1 cDNA construct. Figure 2D shows the mean concentration-induced block comparing ivabradine (open circles for the WT, closed circles for the M651T) to dofetilide (open triangles, WT, and closed triangles, M651T; n values are shown in the legend). At 10 μM ivabradine suppressed 75% of the WT current, whereas at the same concentration, the drug blocked only 10% of the M651T current (Fig. 2D). The mean IC50 was 6.5 ± 10 μM for ivabradine in WT (n = 10), for M651T the IC50 was >120 ±10 μM (n = 7), and beyond the drug’s solubility limit. Note that to accurately measure IC50, a concentration higher than its solubility limit is essential. The impact of the lipid-facing residue M651 to ivabradine was also compared with dofetilide, a prototypical class III antiarrhythmic drug, to assess whether the mutant had a greater impact on response to ivabradine.
Interestingly, the M651T mutation had no impact on the concentration-response relationships to dofetilide. The mean IC50 concentration for dofetilide in WT was 41 nM (n = 4), whereas it was 27 nM for the M651T mutation (n = 5). Thus, the M651T alters the concentration-response relation for ivabradine but not for dofetilide.
The stark contrast in the drug sensitivity for a lipid-facing mutation suggests that the two drugs are accessing the main intracavitary site of hERG1 via different pathways. The pKa values reported for the strongest basic chemical moiety in ivabradine (9.37) is similar but somewhat lower than of dofetilide (9.6). Nonetheless, the cationic form of the drug is a dominant form in the bulk aqueous solution at the physiologic pH values. However, the pKa values of weak-cationic compounds such as ivabradine are not the most accurate predictor of partitioning thermodynamics, since the equilibria between a neutral and a cationic form is a dynamic and environment-dependent process (DeMarco et al., 2018; Dickson et al., 2019).
To investigate thermodynamics of partitioning between aqueous phase and a model hydrophobic environment that mimics the membrane bilayer interior, we studied partitioning of ivabradine and dofetilide in an aqueous buffer and n-hexane. The thermodynamics of transfer for the two drugs showed significant differences. The partitioning free energy data summarized in Table 1 directly demonstrate that ivabradine has probabilities of partitioning nearly equal between aqueous buffer and a bulk hexane (Table 1). The free energy of transfer () between aqueous phase and a bulk hexane for T = 313.15 K are 0.5 ± 0.1 and 8.4 ± 0.2 kJ mol−1 for ivabradine and dofetilide, respectively. Therefore, partitioning data provide direct evidence that ivabradine is a lipophilic molecule, which partitions into the lipid compartment of the plasma membrane compared with dofetilide (Table 1). These findings are in excellent agreement with our previously published modeling studies (Lees-Miller et al., 2015).
Although M651T Slows Deactivation, Its Impact on Ivabradine-Induced Block Appears Independent of Deactivation.
In comparison with WT, M651T mutation slows deactivation; the taus of WT hERG1 and the mutant were 160 ± 40 and 400 ± 40 milliseconds (P < 0.01), respectively. Figure 2, A, B, and E, show the raw deactivation time-course in response to ivabradine in WT (Fig. 2A) and M651T (Fig. 2B) and in a double mutation which rescues deactivation to the WT values (Fig. 2E). Figure 2B shows the drug free time course of deactivation of the M651T mutation compared to the WT (Fig. 2A). The double mutant channel, M651T/T618I restored the deactivation time-course of the channel to values similar to WT hERG1 (compare Fig. 2, A, B, C, and E), both drug-free and with ivabradine. Even the double mutant of the hERG1 channel that contains M651T still shifted the IC50 concentration-response of ivabradine by more than two orders of magnitude (Fig. 2F). The mean IC50 values of M651T were nearly identical to that of the double mutation M651T/T618I (Fig. 2F). These data indicate that decreased pharmacological sensitivity to ivabradine was not the result of slowing deactivation kinetics induced by the M651T mutation.
Impact of Other M651X Substitutions on Ivabradine-Induced Block.
The relationship between baseline drug-free electrophysiological characteristics and IC50 values were evaluated by creating various substitutions at the M651 site (Supplemental Figs. 1–4). We sought to address whether there was a relationship between electrophysiological characteristics and IC50 response to ivabradine. We observed no significant correlation between drug-free voltage-dependence of activation (V1/2), deactivation kinetics, or voltage-dependence of C-type inactivation and the concentration-responsiveness to ivabradine. However, there was a modest correlation between the drug-free time-constant of recovery from inactivation and the mean IC50 values of ivabradine-induced block (R2 = 0.4). We next examined the effect of an IC50 concentration of ivabradine on ion currents elicited in various substitutions at the M651 site. The voltage-dependence of deactivation (during ivabradine treatment) and the kinetics of recovery from inactivation and their associated IC50 values were significant, but these correlations were quite modest. Although many substitutions at this lipid-facing residue (Supplemental Table 2) have significant impact on responsiveness to ivabradine, correlation between “the lipophilicity of the residues” and the mean IC50 values for ivabradine block was limited or even absent (R2 < 0.1–0.4). There were no apparent correlations between the IC50 for ivabradine binding or volume or solvent-accessible area of the residues in the 651 position (Supplemental Fig. 3). It is important to mention that in the absence of ivabradine, several of the M651X mutants exhibited altered gating behavior (Supplemental Figs. 3 and 4).
Therefore, ivabradine-induced block of hERG1 current could not be reduced to a simple drug binding in the vicinity of M651X. The drug action may be coupled to the gating process of the channel or altered accessibility to the intracellular cavity of the channel via some form of allosteric regulation. Essentially, drug-induced shifts in the time-course of recovery from inactivation appear to be a probable determinant of ivabradine potency. This coupling between presence of binding pockets, drug-channel interactions at the lipid-facing surface, and gating dynamics is often described in terms of a complex allosteric mechanism in which change in the topology of the binding pocket or inter-residue interactions modulates access or affinity for the substrate binding pocket (Gordon and Zagotta, 1995).
Temperature Dependence of Ivabradine Blockade.
The measured IC50 values for ivabradine in WT hERG1 differ somewhat between various laboratories. We note that our originally published IC50 value of 6.8 μM differs from that of Melgari et al. (2015). Their measurements were reported at 37°C, whereas our measurements were made at 22°C. Accordingly, to address whether experimental conditions accounted for some of these differences in measured IC50 values we directly compared IC50 values of ivabradine at 37°C versus 22°C. Our measured wild-type IC50 value at 37°C of 3.3 μM is relatively close to the values of 2.07 μM for hERG1 and 3.31 μM for WT-hERG1 1a/1b isoform reported by Melgari et al. (2015). Most importantly, we further re-established that the impact of M651T was not significantly altered by increasing the temperature to 37°C (Fig 3B). In fact, higher temperature may induce more potent block by ivabradine. The measured IC50 of ivabradine decreased from 6.8 to 3.3 μM with an increase in temperature, paralleling more favorable partitioning of ivabradine into the hexane at a higher temperature (Fig. 3; Table 1). Given that ivabradine is a highly lipophilic compound, the drug precipitates at higher concentrations as in these experiments, thereby making it impossible to calculate an exceedingly accurate IC50 value for ivabradine block of M651T, because the solubility of ivabradine prevents evaluations at concentrations >100 μM.
Interplay Between Aromatic Cassette (Y652, F656, and F557) Dynamics and Ivabradine-Induced Block.
A vast number of studies describing molecular determinants of high- to mid-affinity blockade have emphasized an important role played by aromatic residues in the pore-domain cavity of the channel (Y652 and F656) (Duff et al., 1995; Ficker et al., 1998; Lees-Miller et al., 2000, 2015; Mitcheson et al., 2000; Perry et al., 2010; Saxena et al., 2016). These residues are pivotal determinants of class III drug-induced block of IhERG; thus a range of double mutations were created to enhance our understanding of the interplay between those key residues in the pore domain and the lipid-facing residue M651 (Ficker et al., 1998; Lees-Miller et al., 2000; Mitcheson et al., 2000; Perry et al., 2010; Saxena et al., 2016). Figure 4 shows the pharmacological responses of ivabradine to the single and double mutations. Figure 4, A–C, show the effects of adding the M651T mutation to F656C, Y652, or S620T. S620T is a mutant with impaired C-type inactivation phenotype (Herzberg et al., 1998; Perry et al., 2007) that also has a major impact on channel-induced activation or block by small molecules (Ficker et al., 1998; Herzberg et al., 1998; Wu et al., 2014, 2015). We previously reported that the single mutation S620T decreases the pharmacologic responsiveness of IhERG1 to ivabradine (Lees-Miller et al. 2015).
A large concentration of ivabradine (100 μM) only partially blocks the F656C (72% block), Y652A (30%), M651T (30%), and S620T currents (32%), but almost no block was observed at the same concentration of ivabradine with either the M651T/F656C (<1% block), M651T/S620T (1% block), or M651/Y652A (12%) compared with their respective individual mutations as seen in Fig. 4 (P and n values are shown in the legend and statistical analysis of the data are provided in the Supplemental Table 1). The analysis of electrophysiological recordings in Fig. 4 show no additive or synergistic interplay between F656 and Y652 (Fig. 4E). These data indicate that the M651 site substantially modifies the impact of known aromatic mutations in the distal S6. Mutations of F557L and M651T, both, impaired ivabradine-induced block, but surprisingly the block of the double mutant F557L/M651T was similar to the WT values (Fig. 4D and the raw data are shown in Fig. 5A,B). Panel A shows the raw data before and after ivabradine addition. Panel B shows the mean IC50 concentration-reponse plots. Panel C shows the washout. These data indicate F557 in the S5 and M651 in the distal S6 negatively interplay to rescue pharmacologic response of ivabradine. However, no significant negative or positive interplay is observed with single or double mutations of F557L and S620T (Fig. 4F). In review, the M651 residue interplays positively with many other key residues that are structural determinants of ivabradine-induced block. In contrast, M651T interplays negatively with F557L (Fig. 5B). These data suggest a presence of a complex allosteric interaction of ivabradine with residues in hERG1—a novel mechanism established in this work.
Discussion
State-Dependence in the Ivabradine Blockade of hERG1 Currents.
One of the most important determinants of hERG1 block–associated proarrhythmia is the state-dependent kinetics of drug interactions with the channel (Di Veroli et al., 2013a, 2013b; Hill et al., 2014). The landmark feature of hERG1 channel kinetics is the rapid C-type inactivation (Sanguinetti et al., 1995; Schönherr and Heinemann, 1996; Spector et al., 1996). A number of mutations significantly alter or shift the voltage-dependence of C-type inactivation, including the S620T mutant. Therefore, the contribution of inactivation of M651T mutation on ivabradine-induced block was assessed by comparing the extent of block in M651T, S620T, and the double mutant M651T/S620T. By combining a noninactivating pore (S620T) with the M651T mutation, drug block was virtually eliminated. These data indicate that both the process of C-type inactivation and the M651 site in the distal S6 interact to produce at least additive impairment of ivabradine-induced block. Collectively, the results described for mutations in the positions F557, M651, Y652, and F656 indicate that drug blockade depends on the coupling between conformational dynamics of the residues in the pore domain and on a mechanism that involves repacking of the lipid- and/or pore-facing residues. How are state-dependent conformational dynamics of F557, M651, Y652, and F656 coupled to drug access/binding to a pocket? To understand dynamics of these residues we performed 1.2- to 2.5-microsecond all-atom molecular dynamics simulations in an explicit water/membrane system for selected mutants.
Orientations of F557 and F656 are State-Dependent Properties: From Pore-Lining to Lipid-Facing.
The availability of new cryo-EM structures for hERG1 and hEAG1 channels in combination with recent developments in MD simulations allowed us to investigate structural dynamics of the open and closed states along with mutants of interest, hence allowing us to test directly the hypothesis presented above. In the previous and widely accepted mechanistic models, the aromatic residues (Y652 and F656) were postulated to face the permeation pathway and coordinate drugs bound to the intracellular cavity (Chen et al., 2002; Perry et al., 2010). Various structural models with bacterial K+ channels or Shaker-family of K+ channels as their bases emphasized the importance of direct interactions between drugs bound in the water-filled cavity with F656 and/or Y652 (Perry et al., 2010; Wacker et al., 2017). The recent progress in cryo-EM technologies allowed us to model open and closed states of hERG1 channel. Although the model based on hEAG1 represents a closed pore, a cryo-EM structure for hERG1 corresponds to the open state. Importantly, spatial orientation of F656 and F557 side-chains display remarkable state dependence. As predicted by a majority of structural modeling studies, the F656 is pointing toward the inside of the cavity in the closed-state model (Wang et al., 2016; Wacker et al., 2017)). All-atom MD simulations indicate that F656 can rotate away from the permeation pathway and in the open-state model flickers between the intracavitary and the lipid-facing orientations (Fig. 6, A–C). In the open state, the position of F656 aromatic ring is in proximity to F557 and is close enough to form stable π-π stacking interactions at the lipophilic entry pathway.
To assess properties of the putative lipid-entry pathways and the impact of various mutations on the dynamics of the F557-M651X-F656 motif, we performed topological analysis of the production portion of MD simulations (for traces, see Supplemental Figs. 5 and 6) using the MOLEonline pathway analyzer (Berka et al., 2012). Two potential entry pathways were mapped from analysis of MD trajectories. Both pathways are defined by conformational states of F656 and F557 residues and illustrated in Fig. 6A. We analyzed the conformational dynamics of the F557-F656 pair to gain additional insight regarding the flexibility of a tentative lipid-facing binding site and its impact on the accessible volume required for drug diffusion into the primary intracavitary site (Supplemental Fig. 5). The conformational space of F557-F656 interacting pair and the impact of M651T mutation in different states of hERG1 is illustrated in Fig. 6, B–D. The direct interactions between F656 and F557 appear to be only quasi-stable (flickering state) in the WT hERG1 (Fig. 6, B, C, D, E). These interactions establish a large accessible volume for drug binding on the lipophilic site of the intracellular cavity. The interaction of the pair is state-dependent, which may directly support the previously postulated role of F656 in the state-dependent drug blockade of hERG1channel (Chen et al., 2002). Therefore, it is tempting to assign the F656-F557 pair a pivotal role in state-dependent gate control of ivabradine diffusion from the lipid bilayer to the main binding site in the intracellular cavity. The additional analysis of M651T mutant shows that the conformational flexibility of the F656 and F557 pair is significantly impeded with the mutant (Fig. 6C, D, E). Evidently, the M651T mutation significantly stabilizes the pair (F656 and F557) by decreasing the F656 “flickering” frequency (Fig. 6, B and C). Other mutations affecting conformational flexibility of F656 or F557 (Supplemental Figs. 7–14) are also associated with inhibition of hERG1 blockade by ivabradine. The structural states of F656 appear to be important modulators of high-affinity binding for major hERG1 blockers (Supplemental Fig. 14; Supplemental Table 4).
In review, the state-dependent orientation of F656 appears to be an important determinant of a putative lipophilic entry pathway explored by some of the hERG1 blockers, including ivabradine. The all-atom MD-refined models of hERG1 channel in its open and closed states demonstrate that F656 can rotate away from the intracellular cavity toward the lipid bilayer and form hydrophobic interactions with F557 and M651 residues (Fig. 1). Interestingly, the conformational flexibility of F656 (it can be directly modulated by mutations in the position M651 or F557) is in excellent correlation with findings observed experimentally in ivabradine blockade to WT and mutant forms of hERG1.
Induced-Fit Docking and MD-Ensemble Docking Supports Coexistence of Binding Sites in the Pore Cavity and Pore-Lipid Interface.
The extensive sampling of various conformational states of WT and mutant forms of hERG1 channel in different states achieved with Anton 2 platform allowed us to map tentative binding sites for ivabradine. Potential binding sites around the mutated residues were explored through docking of ivabradine to different sections of the protein using hERG1 in different states, also incorporating ensemble docking, where cluster-representatives from all-atom MD simulations were used. Two main binding sites, “lipophilic site” and “internal cavity site,” were further analyzed and compared for WT hERG1 and mutants. The stable and populated binding sites from Induced-Fit Docking are shown in Fig. 7, and a summary of docking energetics for all studied systems are shown in Tables 2 and 3. There were no poses found for ivabradine binding (either neutral of cationic form) to the main internal cavity of the hERG1 for the closed state of the receptor (Table 2). For the closed state of hERG1, ivabradine shows stable binding to the pocket centered at M651, in agreement with the previously proposed lipophilic access site (Lees-Miller et al., 2015) (Table 2). In contrast to the closed state of the channel, the docking to the open state revealed three stable binding modes. Ivabradine was favorably bound to the main pocket in the intracellular cavity, to an area around a fenestration window (in-between α-helices of S6), or to a lipid-facing binding pocket (Fig. 6A). The IFD docking performed to an open state of the hERG1 channel or the MD simulations to an ensemble of open structures show similar binding affinities for ivabradine (Table 3). The results of IFD docking qualitatively agrees with the experimental data (Table 2). At the same time, results from the MD-ensemble docking inherently have large uncertainties in computed binding affinities, rendering comparisons between mutants rather difficult (Table 3). However, the ensemble-based docking simulations allow for better quantification and comparison of relative populations found in each of the three binding sites. The population analysis is essential for understanding the binding processes that involve large and flexible ligands such as ivabradine (Zhao et al., 2010; Shoichet and Kobilka, 2012). The relative populations of binding poses found for identified sites are shown in Fig. 7 and Table 3. Comparing number of poses from the docking simulations, ivabradine preferentially binds to the lipophilic site or the internal cavity of the WT hERG1. However, a substantial number of poses are also found near a potential fenestration window or the “access pathway” (Fig. 7A; Supplemental Fig. 6).
For additional insights, we performed docking simulations on selected mutant systems, specifically M651T, F557L, and F557L/M651T hERG1 mutants. It was found that in both M651T and F557L mutants, ivabradine showed a preferential binding to the intra-cellular cavity of the channel and a very small occupancy (number of poses) in the lipophilic binding pocket. In conclusion, ensemble docking for both mutants M651T and F557L indicate disruption of the lipophilic access site owing to the repacking of this binding pocket caused by the mutations. The mutations in 651 and 557 reduced relative population of binding poses in the binding pocket near a potential fenestration window (access path to the internal cavity path) by 2-fold compared with the WT system. These findings are in excellent agreement with the electrophysiological data presented herein. The data combined for the double mutant M651T/F557L indicate almost WT-like binding affinities in all three mapped binding sites with modest increase in binding affinity for intracellular cavity- and the access binding sites (Fig. 7C). Therefore, the double mutant is expected to have WT-like binding properties in the open state of the channel. We conclude that the flexible nature of F656 side-chain, which can rotate toward the lipid-facing side and the internal cavity of the channel seems to be a determining factor in ivabradine-induced block. MD simulations have also shown that the double mutant F557L/M651T restores the WT-like dynamic of F656 flickering, therefore re-establishing drug occupancy in the lipophilic access site for ivabradine. Subsequent ivabradine docking to WT hERG1 and selected mutants (M651T, F557L, and F557L/M651T) indicated that the lipophilic site is accessible to both open and closed states of the channel. However, the closed state of the channel contains no binding pockets for ivabradine in the intracellular cavity and the “access” binding site near the fenestration window. The simulation data suggest that the structural rearrangements of this binding pocket during the channel’s activation may play a role in ivabradine binding to the lipophilic site and subsequent access to the entry site and ultimately into the intracellular cavity of the channel. Molecular docking simulations performed for the cationic form of ivabradine targeting primary site in the intracellular cavity show comparable binding free energies (to the lipophilic site) and also underline the crucial role of F656 dynamics (Supplemental Table 4) coupled to the complex interactions with F557 and M651 residues.
Conclusions
In summary, we show that the pore-lipid interface mutation (M651T) significantly impairs ivabradine-induced block of the hERG1 current but does not alter dofetilide-induced block. Thus the impact of M651T appears to be specific to ivabradine and emphasizes the underlying important role of a new lipophilic access pathway. The structural mechanisms of the observed lipophilic binding of ivabradine were discerned from a combination of microseconds-long MD and traditional IFD/MD-ensemble docking simulations. The modeling data emphasize the role of M651 as an allosteric modulator of state-dependent hydrophobic interaction between F557 and F656, both of which are well known determinants of ivabradine-induced block. The F557 residue was shown to interact with F656 by forming π-π stacking interactions, an interaction that is disrupted by M651. The overall process of M651 tugging on F656 results in a mobile and flickering F656, which flickers to face the lipid-facing side and the internal cavity. Whereas M651T mutant is unable to interact with F656, it enhances F656 and F557 π-π stacking interactions, resulting in a rigid F656. Thus, M651 controls the quasi-flickering state of F656 and shapes the topology of the binding pocket by controlling orientation of the F656/F557 hydrophobic cassette. MD simulations provided direct evidence that the mutations at these two positions (F557 and M651) cause changes in the topology in the vicinity of the proposed access and/or lipophilic binding sites by altering F656 flexibility and rotations, thereby impacting the number of binding poses and hence limiting accessibility to the main binding pocket in the intracellular cavity of the open-hERG1. We conclude that in WT-hERG1 channel F656, M651, F557 residues act as a dynamic gate that controls the pathway of drug access at the lipid-facing domain into the intracellular cavity of the channel. The experimental data provided firm evidence that the extent of ivabradine blockade recorded for F557L/M651T mutant with WT-like conformational dynamics was essentially similar to the WT hERG1.
Acknowledgments
We thank Dr. Marcela Madrid and Dr. Philip Blood of Pittsburgh Supercomputing Center for their assistance with running molecular dynamics simulations on Anton 2 platform. We express our gratitude to Martin Karplus and CHARMM project for waiving license fees for the CHARMM program package used for analysis of simulation data.
Authorship Contributions
Participated in research design: Perissinotti, Guo, Lees-Miller, Noskov, Duff.
Conducted experiments: Perissinotti, Guo, Kudaibergenova, Ol’khovich, Sharapova, Perlovich.
Contributed new reagents or analytic tools: Guo, Muruve, Gerull, Noskov, Duff.
Performed data analysis: Perissinotti, Guo, Kudaibergenova, Noskov, Duff.
Wrote or contributed to the writing of the manuscript: Perissinotti, Guo, Kudaibergenova, Muruve, Gerull, Noskov, Duff.
Footnotes
- Received December 21, 2018.
- Accepted May 28, 2019.
↵1 L.P., J.G., and M.K. contributed equally to this work.
This work was supported by the Canadian Institutes of Health Research Project Grant FRN-CIHR: 156236 (to J.G., J.L.-M., and H.J.D.); the National Institutes of Health [Grant R01HL128537-01 (S.Y.N., J.G., L.P., and M.K.)]; and the Discovery grant from the Natural Scientific and Engineering Research Council of Canada (to H.J.D.). Computational support for this work was partially provided by West-Grid Canada through a resource allocation award to S.Y.N. M.K. was supported by a Queen Elizabeth II graduate scholarship. Anton 2 computer time was provided by the Pittsburgh Supercomputing Center (PSC) through Grant R01GM116961 from the National Institutes of Health. The Anton 2 machine at PSC was generously made available by D.E. Shaw Research.
↵This article has supplemental material available at molpharm.aspetjournals.org.
Abbreviations
- cryo-EM
- cryoelectron microscopy
- hEAG1
- human ether à go-go 1 gene in reference to α-subunit of Kv10.1 K+ channel
- hERG1
- the human ether-à-go-go-related gene in reference to α-subunit of Kv11.1 K+ channel
- IFD
- Induced-Fit Docking
- Iherg1
- current carried by hERG1
- IKr
- rapid delayed rectifier current
- MD
- molecular dynamics
- WT
- wild-type
- Copyright © 2019 by The American Society for Pharmacology and Experimental Therapeutics