Regular Article
Calculations of the binding free energies of the Comprehensive in vitro Proarrhythmia Assay (CiPA) reference drugs to cardiac ion channels
2023 年 20 巻 2 号 論文ID: e200016
詳細
2023 年 20 巻 2 号 論文ID: e200016
The evaluation of the inhibitory activities of drugs on multiple cardiac ion channels is required for the accurate assessment of proarrhythmic risks. Moreover, the in silico prediction of such inhibitory activities of drugs on cardiac channels can improve the efficiency of the drug-development process. Here, we performed molecular docking simulations to predict the complex structures of 25 reference drugs that were proposed by the Comprehensive in vitro Proarrhythmia Assay consortium using two cardiac ion channels, the human ether-a-go-go-related gene (hERG) potassium channel and human NaV1.5 (hNaV1.5) sodium channel, with experimentally available structures. The absolute binding free energy (ΔGbind) values of the predicted structures were calculated by a molecular dynamics-based method and compared with the experimental half-maximal inhibitory concentration (IC50) data. Furthermore, the regression analysis between the calculated values and negative of the common logarithm of the experimental IC50 values (pIC50) revealed that the calculated values of four and ten drugs deviated significantly from the regression lines of the hERG and hNaV1.5 channels, respectively. We reconsidered the docking poses and protonation states of the drugs based on the experimental data and recalculated their ΔGbind values. Finally, the calculated ΔGbind values of 24 and 19 drugs correlated with their experimental pIC50 values (coefficients of determination=0.791 and 0.613 for the hERG and hNaV1.5 channels, respectively). Thus, the regression analysis between the calculated ΔGbind and experimental IC50 data ensured the realization of an increased number of reliable complex structures.
The complex structures, as well as absolute binding free energies of 25 reference drugs proposed by the Comprehensive in vitro Proarrhythmia Assay consortium to two cardiac ion channels (the human ether-a-go-go-related gene (hERG) potassium channel and human NaV1.5 (hNaV1.5) sodium channel), were predicted. The regression analysis between the calculated binding-free-energy values and experimental half-maximal inhibitory concentration data revealed inaccurate complex structures. By reconsidering the docking poses and protonation states, we observed that the calculated values for 24 and 19 drugs correlated well with their experimental data for hERG and hNaV1.5, respectively.
The assessment of the cardiotoxicity of drugs is a crucial drug-development process because cardiotoxicity is among the most severe side effects of drugs. It is caused by the blockade of cardiac ion channels by drug molecules. The Comprehensive in vitro Proarrhythmia Assay (CiPA) consortium recently proposed a guideline that recommends multiple cardiac-ion-channel assays for the accurate assessment of proarrhythmic risks [1,2]. Particularly, the bindings of drugs to the three cardiac ion channels (the human ether-a-go-go related gene (hERG) potassium channel, human NaV1.5 (hNaV1.5) sodium channel, and human CaV1.2 calcium channel) are crucial for the assessment [3]. However, assessments via in vitro and in vivo experiments are very expensive and time-consuming. Therefore, the in silico predictions of the inhibitory activities of drugs on these channels can reduce the assessment costs and improve the efficiency of the drug-development process.
Channel-blocking molecules inhibit ion permeation by binding to the cavity inside the pore of an ion channel. Therefore, the inhibitory activities of drugs can be basically estimated by predicting their binding structures to the pore cavity and calculating their affinities. In this study, we focused on the hERG and the hNaV1.5 channels because their structures have been experimentally determined via cryo-electron microscopy (cryo-EM) [4–7]. The hERG channel exhibits a tetrameric structure in which each subunit comprises six transmembrane helices (S1–S6). In this structure, helices S5 and S6 form the ion-conducting pore domain, whereas helices S1–S4 form the voltage sensor domain (VSD). It is assumed that drugs penetrate the pore of the hERG channel from the cytoplasmic side. Previous mutagenesis studies revealed that the affinities of various hERG blockers are decreased by the mutations of T623, S624, and V625 near the cytoplasmic end of the selectivity filter, as well as Y652 and F656 in helix S6 [8–17]. Particularly, docking simulations with homology models of the hERG channel have revealed that S624, Y652, F656, and F557 interact directly with various drugs [12–19]. Additionally, a recent study employed cryo-EM to solve a complex structure comprising hERG and astemizole, indicating that the channel blockers directly and sterically inhibit K+ flux [5].
The structure of the hNaV1.5 channel has also been solved in the ligand-free form, as well as in a complex containing quinidine [6,7]. The α subunit of the hNaV1.5 channel is formed by four pseudosymmetric domains (domains I–IV), and each domain comprises six transmembrane helices (S1–S6) that are further divided into the VSDs (helices S1–S4) and pore region (helices S5–S6). The loops between S5 and S6 of the four domains form the selective filter. This complex structure indicates that the bound quinidine molecule interacts with F1760. Additionally, the extant mutagenesis studies revealed that L1462 and L1466 on helix S6 of domain III, as well as F1760 and Y1767 on helix S6 of domain IV, can participate in the binding of local anesthetic drugs [20–28]. Notably, F1760 is a common binding residue for most anesthetics. Therefore, small-molecule drugs can bind inside the pore of the hNaV1.5 channel, interacting with these residues.
In our previous study [29], we calculated the binding free energies (ΔGbind) of 12 structurally diverse drugs to the hERG channel by combining molecular docking simulation and a molecular dynamics (MD)-simulation-based free-energy-calculation method known as the massively parallel computation of absolute binding free energy with well-equilibrated states (MP-CAFEE) [30,31]. The calculated values correlated well with the experimental inhibition constants. In this study, we applied this method to CiPA reference drugs and calculated their ΔGbind values to the hERG and hNaV1.5 channels. The prediction accuracies were evaluated by comparing the calculated ΔGbind and experimental half-maximal inhibitory concentration (IC50) values.
As described in our previous study [29], the structural model of the hERG channel was constructed using its cryo-EM-determined structure (PDB ID: 5VA2). Briefly, the structures of the missing residues of the extracellular loops (residues 578–582 and 598–602) were modeled using Modeller 9.18 [32], and potassium ions were placed at the S1 and S3 sites of the selectivity filter. Further, the structural model of the hNaV1.5 channel was constructed using the cryo-EM-determined structure of the complex containing quinidine (PDB ID: 6LQA). Next, the quinidine was removed from the model, and disulfide bonds were formed between C280 and C326, C335 and C341, C906 and C915, C1363 and C1384, and C1728 and C1742. Furthermore, the protonation states of the ionizable residues at pH 7 were determined using the Protein Preparation Wizard of the Schrödinger Suite.
We utilized 25 of the 28 CiPA drugs for the docking simulations (Supplementary Table S1) [2,33]. We excluded three drugs, namely erythromycin (classified as macrolides), as well as nifedipine and nitrendipine (both classified as dihydropyridines), from the dataset, because they might modulate the channels via different mechanisms from blockade of the pore. It is believed that erythromycin allosterically modulates the blockade by the other pore inhibitors by binding to an extracellular site in the hERG channel [34–36]. The mutation of Y652 of the hERG channel does not affect the binding of erythromycin, indicating that it does not penetrate the pore cavity [34]. Dihydropyridines are known to allosterically modulate calcium channels [37]. The complex structures of voltage-gated calcium channels with several dihydropyridine drugs have demonstrated the binding of the drug molecules to the outer side of the pore domain [38–40].
The structures of the 25 drugs were obtained from the PubChem database [41]. The protonation state of each drug at pH 7 was predicted, and each structure was optimized by the OPLS3 force field [42]. We conducted these calculations with the LigPrep module of the Schrödinger Suite [43].
Docking SimulationsWe performed docking simulations using the Glide module of the Schrödinger Suite. The drug molecules were docked to the hERG channel by the induced fit docking (IFD) protocol [44]. We utilized the pore domain (residues 544–671) of the hERG-channel structures as the receptor, and S624, Y652, and F656 were set as the flexible residues of the IFD protocol. The center of the box that defines the search space was set to the center of masses of Y652 in the four subunits. The extra-precision (XP) mode [45,46] of the Glide module was used in the final redocking step of the IFD. The maximum number of outputs was set to 50 for each ligand, and the docking poses were ranked by the IFD score. The docking simulations for the hNaV1.5 channel were performed in the Glide XP mode. We utilized the structure of the pore domain (residues 233–429, 823–943, 1317–1481, and 1640–1781) as the receptor. The center of mass of F1760 was used as the center of the box that defines the search space.
MD SimulationsWe performed MD simulations for the ligand-bound structures of the hERG channel, as well as the ligand-free and ligand-bound structures of the hNaV1.5 channel, in a solvated lipid bilayer environment. Water molecules were added to the central cavity of the channel before embedding the channel structure in the lipid bilayer using the CHARMM-GUI server [47,48]. The hERG-channel system comprised a pore domain of the channel, about 290 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) molecules, 18,000 water molecules, K+, and Cl–. Regarding the hNaV1.5 channel, we constructed two systems with different protein lengths. One was the pore domain-only system comprising a pore domain of the channel, 290 POPC molecules, 27,000 water molecules, Na+, and Cl–, and the other was the VSD–pore system comprising VSDs and the pore domain of the channel structure (residues 119–429, 698–943, and 1188–1781), 610 POPC molecules, 58,000 water molecules, Na+, and Cl–. Further, we performed MD simulations for the ligand-only systems, each of which comprised one ligand molecule and approximately 6,700 water molecules. K+ and Cl– were added to each system to ensure that the net charge and concentration of the whole system were 0 and approximately 0.15 M, respectively. We employed the CHARMM36m-WYF [49,50] force field for the proteins [49], ions [51], and POPC [52]. The TIP3P model [53] was employed to model the water molecules. The force-field parameters of the drugs were extracted from the CHARMM General Force Field [54] (CGenFF); the topology and parameter files were generated by the CHARMM-GUI ligand modeler [55]. Each system was energy-minimized and equilibrated via the CHARMM-GUI protocol. We performed MD simulations under the constant-NPT condition to generate equilibrium ensembles (production runs). Additionally, the temperature was maintained at 303.15 K by a Nose–Hoover thermostat [56,57]. The pressure was maintained at 1×105 Pa by a Parrinello–Rahman barostat [58]. The electrostatic interactions were calculated via the particle-mesh Ewald (PME) method with a real-space cutoff of 1.4 nm [59]. Furthermore, the van der Waals interactions were calculated using a modified Lennard–Jones potential (where the force was smoothly switched to zero between 1.0 and 1.2 nm). The lengths of the bonds containing the hydrogen atoms were constrained by the linear constraint solver (LINCS) algorithm to allow the utilization of a time step of 2 fs [60]. The production runs for the ligand-free and ligand-only systems proceeded for 500 and 20 ns, respectively. All the MD simulations were performed using Gromacs 2019 [61].
Calculation of the Absolute Binding Free EnergiesThe ΔGbind values were calculated via the MP-CAFEE method, following our previous study [29]. We performed 50-ns MD simulations five times using different initial velocities for each channel with the top-ranked docking pose of each drug. Thereafter, the moving average of the intermolecular interaction energy (the sum of the van der Waals and electrostatic interaction energies) between the protein and ligand was calculated using a 2-ns window. Next, the snapshot with the lowest interaction energy was selected as the initial structure for the free-energy calculation. Regarding the ligand-only systems, the final structure of the 20-ns production run, which was performed once per ligand, was employed as the initial structure. The free-energy differences between the λ=0 (where the ligand interacts fully with its surroundings) and λ=1 (where the intermolecular interactions between the ligand and others are completely decoupled) states were calculated for the protein–ligand complexes (ΔGcomp) and ligand-only systems (ΔGlig). Here, we considered 32 intermediate states between the λ=0 and λ=1 states to ensure the sufficient proximity of the two adjacent states, between which the free-energy difference could accurately be calculated. The MD simulations of the free-energy calculation were performed six times at different initial velocities. The durations of the MD simulation of each state were 2 and 1.2 ns for the protein–ligand complex and ligand-only system, respectively. The temperature and pressure were maintained at 303.15 K and 1 bar by the Nose–Hoover thermostat and Berendsen barostat, respectively [62]. The free-energy difference was calculated via the Bennett acceptance-ratio method [63,64]. The ΔGbind value was calculated using the expression: ΔGbind=ΔGcomp–ΔGlig.
In our previous study, we demonstrated that the pore-domain structure of the hERG channel was maintained without the VSD and C-linker in a 500-ns MD simulation [29]. Here, we compared the stabilities of the pore-domain structures of the hNaV1.5 channel in an MD simulation of the two systems (the pore-domain-only and VSD–pore systems). Each system was embedded in a solvated lipid bilayer, after which a 500-ns MD simulation was performed. The Cα root-mean-square deviations (RMSDs) of the transmembrane regions of the pore domain (residues 251–273, 358–386, 390–423, 840–870, 884–910, 912–942, 1332–1356, 1405–1430, 1444–1481, 1655–1680, 1697–1721, and 1745–1778) from the initial structures were calculated (Supplementary Figure S1). The RMSD values of the VSD–pore system were around 2 Å throughout the simulation. Although the RMSD values of the pore domain-only system were slightly larger than those of the VSD–pore system, they remained at around 3 Å, indicating that the pore-domain structure was stably maintained without the VSD. Therefore, we utilized the pore-domain structures of the hERG and hNaV1.5 channels in the subsequent analyses.
Molecular Docking SimulationsWe performed the molecular docking simulations of the drugs to the pore cavities of the hERG and hNaV1.5 channels using the Glide module of the Schrödinger Suite. For the hERG channel, we applied the IFD protocol for the docking simulation (following our previous study), in which the center of the masses of Y652 in the four subunits was used as the center of the search space, and S624, Y652, and F656 were set as the flexible residues. We increased the maximum number of outputs to 50 because our previous study indicated that many conformations must be explored to ensure an accurate prediction [29]. Considering that six of the drugs analyzed in our previous study are among the CiPA reference drugs, they were analyzed again in this study using updated simulation settings. Figure 1a shows a plot of the Glide docking score (Gscore) of the top-ranked pose against the negative of the common logarithm of the experimental IC50 (pIC50) [33]. The scores correlated weakly with the experimental data (R2=0.314). We also analyzed the spatial distributions of the aromatic rings and positively charged amines since most of the drugs contained these groups (Figure 2), and the results revealed that the geometric centers of the aromatic rings were close to Y652 or F656, demonstrating the criticality of π–π interactions to ligand binding (Figure 2a). The distribution of the nitrogen atoms of the positively charged amines indicated that most of them were located near S624 or Y652, generating hydrogen bonds or cation–π interactions with the hydroxy group of S624 or aromatic ring of Y652, respectively (Figure 2b). Consistent with previous studies, these results indicated that hydrogen bonding with S624 and Y652, as well as π–π and cation–π interactions with Y652 and F656, respectively, contribute to drug binding.
Plots of the Gscore of the top-ranked poses against the experimental pIC50 values of the 25 compounds for the hERG (a) and hNaV1.5 (b) channels.
Distributions of the aromatic rings (a) and positively charged amines (b) in all the top-ranked poses for the hERG channel. The geometric centers of the aromatic rings, as well as the nitrogen atoms of the amines of the ligand molecules, are shown as magenta balls. The protein backbone structures are shown as ribbons. The sidechains of S624, Y652, and F656 are shown in a stick representation, and their carbon atoms are colored light blue, green, and orange, respectively. The hydrogen and oxygen atoms are colored white and red, respectively. The 90°-rotated views are shown on the right panels. All the molecular graphic images were produced with UCSF Chimera [65].
Regarding the hNaV1.5 channel, the rigid receptor docking was used with the Glide XP score because the preliminary calculations indicated that IFD performed worse than rigid receptor docking (Supplementary Figure S2). We set the center of the search space as the center of the masses of F1760 on helix S6 of domain IV, which is involved in the binding of local anesthetics [20–28]. The Gscores of the top-ranked poses did not correlate with the experimental pIC50 values (R2=8.00×10–4) (Figure 1b). Figure 3 shows the spatial distributions of the aromatic rings and positively charged amines, and the results revealed that the geometric centers of the aromatic rings were at the center of the pore or spaces between domains III and IV and between domains II and III, close to putative binding residues, L1462 and L1466 (on helix S6 of domain III) and F1760 and Y1767 (on helix S6 of domain IV). The distribution of the nitrogen atoms of the positively charged amines indicated that most of them were at the upper part of the pore. A previous study revealed that the binding of a positively charged molecule under the selectivity filter impairs that of Na+ to the selective filter [66]. Therefore, the predicted binding modes can inhibit ion permeation. Although the docking scores did not correlate with the experimental data, the predicted models were overall consistent with the previous findings.
Distributions of the aromatic rings (a) and positively charged amines (b) in all the top-ranked poses for the hNaV1.5 channel. The geometric centers of the aromatic rings, as well as nitrogen atoms of the amines of the ligand molecules, are shown as magenta balls. The protein backbone structures are shown as ribbons. The side chain of F1760 is shown in a green-stick representation. The 90°-rotated views are shown on the right panels. Roman numbers I–IV indicate the positions of domains I–IV, respectively.
Next, we calculated ΔGbind of the drugs via the MP-CAFEE method [30]. The complex structures exhibiting the top-ranked docking poses were used as the initial structures for the MD simulations. The ΔGbind values are listed in Supplementary Table S2 and plotted against the experimental pIC50 values in Figure 4. Notably, the calculated values for the six common drugs were not the same as those of the previous study because we used a newer force field and a different cutoff scheme for the non-bonded interactions [29]. The coefficient of determination for the hERG channel was 0.465, indicating that the correlation improved upon calculating the ΔGbind values. The root mean squared error (RMSE) from the regression line was 5.03 kcal mol–1. The calculated values for ibutilide, metoprolol, ondansetron, and quinidine deviated significantly from the regression line of the plot (Figure 4a).
Plots of the calculated ΔGbind values vs. the values of pIC50 of the 25 compounds for the hERG (a) and hNaV1.5 (b) channels. The black-colored data points represent the outliers. The solid line represents the regression line calculated from all the data. The dashed line represents a regression line after excluding 10 outliers.
The calculated affinity of ibutilide was weaker than what was anticipated from the regression line. To obtain a more energetically favorable structure, all the docking poses, which were generated for ibutilide, were classified into clusters via the average-linkage method using the RMSD values of the non-hydrogen atoms as the distance measure (Supplementary Figure S3). Thus, five clusters were obtained; the top-ranked and second-ranked poses (poses 1 and 2) belonged to the largest cluster, whereas the third-ranked pose (pose 3) belonged to the second-largest cluster. Therefore, we calculated ΔGbind of pose 3. The equilibrium structures of poses 1 and 3 are shown in Figures 5a and 6a, respectively. Pose 1 formed hydrogen bonds with S624 and T623, a π–π interaction with Y652, and cation–π interactions with Y652 and F656. Conversely, pose 3 formed cation–π interactions with Y652 and F656, as well as a hydrogen-bond interaction with Y652. The calculated ΔGbind of pose 3 was –22.50 kcal mol–1, which was lower than that of pose 1 (–15.38 kcal mol–1).
Equilibrium structures of ibutilide (a), ondansetron (b), quinidine (c), and metoprolol (d) for the hERG channel. The ligand molecules and sidechains of F557, S624, Y652, and F656 are shown in a stick representation. The carbon atoms of the ligand molecules, F557, S624, Y652, and F656, are colored magenta, dark gray, light blue, green, and orange, respectively. The hydrogen, nitrogen, oxygen, and sulfur atoms are colored white, blue, red, and yellow, respectively. The black dashed lines represent the hydrogen bonds, π–π interactions, or cation–π interactions, as identified by the Maestro module of the Schrödinger Suite.
Equilibrium structures of ibutilide (a), ondansetron (b), and quinidine (c) for the further calculation of the hERG channel. The ligand molecules and sidechains of F557, S624, Y652, and F656 are shown by a stick representation. The carbon atoms of the ligand molecules, F557, S624, Y652, and F656, are colored magenta, dark gray, light blue, green, and orange, respectively. The hydrogen, nitrogen, oxygen, and sulfur atoms are colored white, blue, red, and yellow, respectively. The black dashed lines represent the hydrogen bonds, π–π interactions, or cation–π interactions, as identified by the Maestro module of the Schrödinger Suite.
In the calculation for ondansetron, we first assumed that the nitrogen atom on the imidazole ring was protonated and that the drug exhibited a net charge of +1 (Figure 5b). However, ondansetron can exist in the neutral state at the physiological pH because its experimental pKa is 7.4 [67]. Therefore, we performed the docking simulation and the ΔGbind calculation for the neutral state of ondansetron. In the equilibrium structure, the imidazole and indole rings were located near S624 and F656, respectively (Figure 6b). The ΔGbind value calculated with the MP-CAFEE method was –17.84 kcal mol–1, which is closer to the regression line than that for the protonated state (–27.91 kcal mol–1).
In the equilibrium structure of quinidine, the ligand molecule was located at the upper part of the cavity away from F656, interacting with Y652 (Figure 5c). However, a previous experiment confirmed that the mutation of F656 more strongly affects the inhibitory activity of quinidine than the mutation of Y652 [9]. Since pose 5 obtained the best score in the poses interacting with F656, we calculated its ΔGbind. In the equilibrium structure, the positively charged amine was near S624, and the aromatic ring was near Y652 and F656 (Figure 6c). The calculated ΔGbind value was –16.42 kcal mol–1, which was also closer to the regression line than that of pose 1 (–27.36 kcal mol–1). Although pose 1 is more energetically favorable than pose 5, the binding site of pose 1 is located farther from the intracellular entrance of the pore than that of pose 5, which would make it difficult for quinidine to reach the binding site of pose 1.
Metoprolol was located at the center of the pore and formed π–π and π–cation interactions with Y652 in its equilibrium structure, corresponding with the general features of many other drugs (Figure 5d). However, a previous study revealed that the mutation of Y652 or F656 only exerts a small effect on the binding of metoprolol [68]. Therefore, metoprolol might bind outside the pore cavity, and further information on the binding site is necessary for an accurate calculation. By using the newly calculated ΔGbind values for ibutilide, ondansetron, and quinidine and excluding the data for metoprolol, the coefficient of determination of the correlation with the experimental pIC50 values increased to 0.791 and the RMSE was improved to 3.66 kcal mol–1 (Figure 7a).
Plot of ΔGbind vs. pIC50 after the additional calculations. The results for the hERG (a) and hNaV1.5 (b) channels. The data points colored black represent the outliers, and the data points colored orange represent the results of the additional calculations. The solid line represents the regression line calculated from the data points without the outliers.
Regarding the hNaV1.5 channel, the coefficient of determination of the correlation between the calculated ΔGbind and the experimental pIC50 values was 8.07×10–5; it was not improved by the ΔGbind calculation. The RMSE from the regression line was 7.85 kcal mol–1. In the plot shown in Figure 4b, ten drugs, namely astemizole, azimilide, cisapride, dofetilide, droperidol, mexiletine, quinidine, risperidone, tamoxifen, and vandetanib, exhibited significant deviations from the regression line.
The calculated value of astemizole was significantly lower than what was anticipated from the regression line. Astemizole was predicted to be in the doubly protonated state (net charge=+2), in which the nitrogen atoms in the piperidine and benzimidazole rings were protonated. However, the nitrogen atom in the benzimidazole ring can be deprotonated at the physiological pH because its experimental pKa was 6.2 [69]. Therefore, we performed the docking simulation and the ΔGbind calculation for the singly protonated state exhibiting a net charge of +1. Figures 8a and 9a show the equilibrium structures for the doubly and singly protonated states, respectively. In the equilibrium structure, astemizole was closer to F1760 in the singly protonated state than in the doubly protonated state. The calculated ΔGbind was –22.96 kcal mol–1, which is closer to the regression line than the original value (–51.39 kcal mol–1).
Equilibrium structures of astemizole (a), droperidol (b), quinidine (c), and mexiletine (d) for the hNaV1.5 channel. The ligand molecules, as well as the sidechain of F1760, are shown in a stick representation. The carbon atoms of the ligand molecules and F1760 are colored magenta and green, respectively. The hydrogen, fluorine, nitrogen, and oxygen atoms are colored white, yellowish green, blue, and red, respectively.
Equilibrium structures of astemizole (a), droperidol (b), and quinidine (c) for the additional calculation of the hNaV1.5 channel. The ligand molecules, as well as the sidechain of F1760, are shown in a stick representation. The carbon atoms of the ligand molecules and F1760 are colored magenta and green, respectively. The hydrogen, fluorine, nitrogen, and oxygen atoms are colored white, yellowish green, blue, and red, respectively.
The calculated value for droperidol was higher than what was anticipated from the regression line. This molecule was not protonated (net charge=0) in the initial calculation. However, the nitrogen atom of the piperidine ring of droperidol can be protonated at the physiological pH because its experimental pKa was 7.6 [70]. Therefore, we performed the docking simulation and the ΔGbind calculation for the protonated state exhibiting a net charge of +1. Figures 8b and 9b show the equilibrium structures of the unprotonated and protonated states, respectively. Droperidol became closer to F1760 in the protonated state than in the unprotonated state; its calculated ΔGbind was –20.83 kcal mol–1, which was energetically more favorable and closer to the regression line than that obtained for unprotonated droperidol (–14.51 kcal mol–1).
The calculated value for quinidine was also higher than what was anticipated from the regression line. Considering that the experimental structure of the hNaV1.5–quinidine complex is available, we performed its MD simulations. In the docking-structure-based equilibrium structure, the positively charged amine was near the rim of the pore (Figure 8c). However, in the experimental-structure-based equilibrium structure, the amine was located near the center of the pore (Figure 9c). The experimental-structure-based calculated ΔGbind was –20.77 kcal mol–1, indicating that the experimental structure was energetically more favorable than the docking structure (–9.51 kcal mol–1).
In the equilibrated structure of mexiletine (Figure 8d), the ligand molecule was located near F1760 and L1462, corresponding to the results of the previous study [20]. However, the calculated ΔGbind value was higher than what was experimentally anticipated. A previous study indicated that lidocaine, which shares a common structure with mexiletine, binds to the hNaV1.5 channel in a 2:1 molar ratio [66]. As mexiletine is smaller than lidocaine, the pore cavity of the hNaV1.5 channel could accommodate two mexiletine molecules. The interactions between the two mexiletine molecules might be necessary for achieving a strengthened affinity to the hNaV1.5 channel.
The calculated ΔGbind values of the other six drugs (azimilide, cisapride, dofetilide, risperidone, tamoxifen, and vandetanib) were lower than what was anticipated from the regression line. These drug molecules formed close interactions with protein atoms in the narrow space between domains II and III or III and IV (Supplementary Figure S4). One possible reason for the overestimation of their affinities was the extreme closeness of their interactions. To determine docking poses forming relatively weak interactions with the protein, cluster analysis was performed for each drug (Supplementary Figure S5). The results indicated that poses 2, 6, 9, and 5 of azimilide, cisapride, dofetilide, and vandetanib, respectively, were located near the center of the pore and formed relatively weak interactions with the protein. Additionally, we calculated the ΔGbind values of these complex structures (Supplementary Table S2, Supplementary Figure S6), although the calculated values were still much lower than the regression line, except for dofetilide. Therefore, the five drugs (azimilide, cisapride, risperidone, tamoxifen, and vandetanib) might bind to different sites from that considered in this study. For example, a complex structure comprising NaVMs and tamoxifen revealed that eight tamoxifen molecules bound outside the pore near its cytoplasmic exit [71]. Further information on binding sites is required for the accurate predictions of the ΔGbind values of these drugs. By using the newly calculated ΔGbind values of astemizole, dofetilide, droperidol, and quinidine, as well as excluding the data of the other six drugs, the coefficient of determination of the correlation with the experimental pIC50 values increased to 0.613 and the RMSE was improved to 2.08 kcal mol–1 (Figure 7b).
In this study, we performed docking simulations for 25 drugs, which were selected from the CiPA reference drugs, to generate complex structures with two cardiac ion channels, hERG and hNaV1.5. ΔGbind of the top-ranked pose of each drug–channel pair was calculated. The regression analysis between the calculated ΔGbind values and experimental pIC50 values revealed that the calculated values of four and ten drugs deviated significantly from the regression lines for the hERG and hNaV1.5 channels, respectively. Further, we reconsidered the docking poses, as well as protonation states of the drugs based on their experimental data and recalculated their ΔGbind values. Thus, we determined that the calculated ΔGbind values of 24 and 19 drugs correlated well with their experimental pIC50 values, with coefficients of determination of 0.791 and 0.613 for the hERG and hNaV1.5 channels, respectively.
The experimental IC50 values alone could not offer detailed information about the interactions between the drugs and ion channels. Moreover, it is still challenging to accurately predict the complex structures between them. As demonstrated in this study, the calculation of the ΔGbind values for the docking structures of many drugs, as well as analyses of the correlation between their calculated values and experimental data facilitated the isolation of reliable docking structures from incorrect ones as the ΔGbind values of the incorrect docking structures must exhibit significant deviations from the regression line. Notably, the calculated ΔGbind value of the experimental structure of the hNaV1.5–quinidine complex was close to the regression line, confirming the reliability of the complex structures that yielded ΔGbind values that correlated well with the experimental data.
Notably, the experimental data employed in this study were the IC50 values because experimental IC50 values are more available than the dissociation constants, which were used in the previous study [29]. As the dissociation constant can be converted into ΔGbind, it can be directly compared with the calculated ΔGbind. Conversely, pIC50 does not strictly correlate with ΔGbind because IC50 values depend on the concentration of the ion channel. This is among the causes of the slightly worse correlation obtained for the hERG channel in this study than that obtained in the previous study.
In this study, reliable complex structures could not be obtained for one (metoprolol) and six (azimilide, cisapride, mexiletine, risperidone, tamoxifen, and vandetanib) drugs for the hERG and hNaV1.5 channels, respectively. Additionally, three drugs (erythromycin, nifedipine, and nitrendipine) in the CiPA reference drugs were excluded from the calculations. When the experimental information about the binding sites of these drugs become available, reliable predictions of their complex structures will become feasible. As reliable data on the complex structures of various drugs are accumulated, it will become possible to predict the complex structures and inhibitory activities of newly developed drugs by applying the regression equation to the ΔGbind values calculated using the predicted complex structures.
The authors declare no conflict of interest.
T.T. designed the research. T.N. conducted the structural modeling and simulations, as well as analyzed the data. Both authors drafted the manuscript and gave final approval of the version to be published.
The evidence data that were generated and analyzed in the study are available from the corresponding author upon reasonable request.
The authors thank Drs. Toshiaki Hisada and Seiryo Sugiura at UT-Heart Inc. for providing information about the CiPA reference drugs and their experimental data. This work was supported by MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Overcoming heart failure pandemic with innovative integration of multi-scale heart simulator and large-scale clinical data, JPMXP1020200202) and used computational resources of supercomputer Fugaku provided by the RIKEN Center for Computational Science (Project IDs: hp200121, hp210180, and hp220178). This work was also supported in part by Research Support Project for Life Science and Drug Discovery (Basis for Supporting Innovative Drug Discovery and Life Science Research (BINDS)) from AMED under Grant Numbers JP21am0101107 and JP22ama121027.
