DECEMBER 2017
ISTANBUL TECHNICAL UNIVERSITY GRADUATE SCHOOL OF SCIENCE ENGINEERING AND TECHNOLOGY
Ph.D. THESIS
IN SILICO DESIGN OF HERG NON-BLOCKER COMPOUNDS WITH RETAINED PHARMACOLOGICAL ACTIVITY USING MULTI-SCALE
MOLECULAR MODELING APPLICATIONS
Gülru KAYIK
Chemistry Department Chemistry Programme
DECEMBER 2017 Chemistry Department Chemistry Programme
ISTANBUL TECHNICAL UNIVERSITY GRADUATE SCHOOL OF SCIENCE ENGINEERING AND TECHNOLOGY
IN SILICO DESIGN OF HERG NON-BLOCKER COMPOUNDS WITH RETAINED PHARMACOLOGICAL ACTIVITY USING MULTI-SCALE
MOLECULAR MODELING APPLICATIONS
Ph.D. THESIS Gülru KAYIK (509112008)
Thesis Advisor: Prof. Dr. Nurcan TÜZÜN Thesis Co-Advisor: Assoc. Prof. Dr. Serdar DURDAĞI
ARALIK 2017 Kimya Anabilim Dalı
Kimya Programı
İSTANBUL TEKNİK ÜNİVERSİTESİ FEN BİLİMLERİ ENSTİTÜSÜ
HERG BLOKER OLMAYAN FARMAKOLOJİK AKTİVİTESİ KORUNMUŞ BİLEŞİKLERİN ÇOK BOYUTLU MOLEKÜLER MODELLEME
UYGULAMALARI İLE İN SİLİKO TASARIMI
DOKTORA TEZİ Gülru KAYIK
(509112008)
Tez Danışmanı: Prof. Dr. Nurcan TÜZÜN Eş Danışman: Doç. Dr. Serdar DURDAĞI
Gülru KAYIK, a Ph.D. student of İTU Graduate School of Science Engineering and Technology student ID 509112008, successfully defended the thesis entitled “IN SILICO DESIGN OF HERG NON-BLOCKER COMPOUNDS WITH RETAINED PHARMACOLOGICAL ACTIVITY USING MULTI-SCALE MOLECULAR MODELING APPLICATIONS”, which she prepared after fulfilling the requirements specified in the associated legislations, before the jury whose signatures are below.
Thesis Advisor : Prof. Dr. Nurcan TÜZÜN ... Istanbul Technical University
Co-advisor : Assoc. Prof. Dr. Serdar DURDAĞI ... Bahçeşehir University
Jury Members : Prof. Dr. Mine YURTSEVER ... Istanbul Technical University
Prof. Dr. Kemal YELEKÇİ ... Kadir Has University
Assis. Prof. Dr. Bülent BALTA ... Istanbul Technical University
Assoc.Prof. Dr. Fethiye Aylin SUNGUR ... Istanbul Technical University
Prof. Dr. Safiye ERDEM ... Marmara University
FOREWORD
First of all, I would like to thank my Ph.D. thesis advisor Prof. Dr. Nurcan Tüzün and co-advisor Assoc. Prof. Dr. Serdar Durdağı for their kind concern, recommendations and supports during the course of my Ph.D. studies.
I would like to present my acknowledgements to Istanbul Technical University Research Fund BAP (Project numbers: 38208 and 30492) and the National Center for High Performance Computing of Turkey (UHEM) under Grant 10982010 for supporting this thesis and providing the related computer resources.
The numerical calculations reported in this thesis were partially performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources).
I would also like to thank The Scientific and Technological Research Council of Turkey (TUBITAK) for granting me the 2214-A Research Grant and providing financial support during the course of my Ph.D. thesis.
December 2017 Gülru KAYIK
TABLE OF CONTENTS Page FOREWORD ... vii TABLE OF CONTENTS ... ix ABBREVIATIONS ... xiii SYMBOLS ... xv
LIST OF TABLES ... xvii
LIST OF FIGURES ... xix
SUMMARY ... xxv
ÖZET... ... xxvii
1. INTRODUCTION ... 1
2. IN SILICO DESIGN OF NOVEL HERG-NEUTRAL SILDENAFIL LIKE PDE5 INHIBITORS ... 3
2.1 Introduction ... 3
2.2 Methods ... 5
2.2.1 Molecular docking simulations ... 6
2.2.2 Fragment-based de novo drug design & virtual screening ... 7
2.2.3 MD simulations and post-processing MD analyses ... 8
2.2.4 MD simulations of the target receptor: PDE5 in its apo state and bound with its inhibitors ... 8
2.2.5 MD simulations of hERG K+ ion channel: Apo state and bound with PDE5 inhibitors ... 9
2.2.6 Molecular Mechanics/Generalized Born surface area (MM/GBSA) calculations ... 11
2.3 Results and Discussion ... 11
2.3.1 Analysis of the key interactions of sildenafil with the target receptor (PDE5) ... 13
2.3.2 Comparison of used docking tools in terms of predicting the binding positions of “Sildenafil” in the central cavities of hERG1 channel ... 15
2.3.2.1 GOLD ... 16
2.3.2.2 AutoDock ... 21
2.3.2.3 MOE ... 24
2.3.3 In silico Alanine mutagenesis study ... 25
2.3.4 General statements on the binding energy predictions derived from GOLD, AutoDock and MOE ... 27
2.3.5 Binding interactions of Vardenafil and Tadalafil with the hERG K+ ion channel ... 27
2.3.6 Virtual screening results ... 29
2.3.7 MM/GBSA analyses ... 30
2.4 Conclusions ... 30 3. INVESTIGATION OF PDE5/PDE6 AND PDE5/PDE11 SELECTIVE
POTENT TADALAFIL-LIKE PDE5 INHIBITORS USING
COMBINATION OF MOLECULAR MODELING APPROACHES, MOLECULAR FINGERPRINT-BASED VIRTUAL SCREENING
PROTOCOLS AND STRUCTURE-BASED PHARMACOPHORE
DEVELOPMENT ... 37
3.1 Introduction ... 37
3.2 Methods ... 41
3.2.1 Ligand and protein preparations ... 41
3.2.2 Virtual library screening... 41
3.2.3 Flexible molecular docking simulations ... 41
3.2.4 Molecular Dynamics simulations ... 42
3.2.5 Molecular Mechanics Generalized Born Solvation (MM/GBSA) Calculations ... 42
3.3 Results and discussion ... 43
3.3.1 Validation of the docking methodology ... 43
3.3.2 Constructing the homology models of the catalytic domains of PDE6 (amino acid residues: 482-816) and PDE11 (amino acid residues: 587-910) .... 43
3.3.3 Binding affinity and binding pattern analysis of the hit compounds and tadalafil with PDE5, PDE6 and PDE11 ... 47
3.3.4 MD simulations of apo and holo states of PDE5, PDE6 and PDE11 bound with the selected hit compounds (ZINC02120502 and ZINC16031243) and tadalafil ... 59
3.3.5 MM-GBSA calculations ... 62
3.3.6 hERG K+ ion channel activity of the compounds ... 63
3.3.7 E-Pharmacophore studies ... 65
3.4 Conclusions ... 68
4. STRUCTURAL INVESTIGATION OF VESNARINONE AT THE PORE DOMAINS OF OPEN AND OPEN-INACTIVATED STATES OF HERG1 K+ CHANNEL ... 71
4.1 Introduction ... 71
4.2 Computational Methods ... 73
4.2.1 Protein-Ligand docking calculations ... 73
4.2.2 MD simulations ... 74
4.2.3 Principal component analysis (Covariance analysis) of the MD trajectories ... 75
4.3 Results and Discussion ... 76
4.3.1 Protein-ligand docking calculations ... 76
4.3.2 MD simulations ... 79
4.3.2.1 General statements on the backbone and ligand RMSD evaluations and ‘Short Range’ energetics ... 83
4.3.2.2 MD simulations initiated with GOLD/GoldScore fitness functions: analysis of the trajectories of vesnarinone-hERG1 K+ channel model in its open-inactivated state ... 86
4.3.2.3 MD simulations initiated with GOLD/GoldScore fitness functions: analysis of the trajectories of vesnarinone-hERG1 K+ channel model in its open-state ... 89
4.3.2.4 Principal Component Analysis (PCA) & comparison of time-dependent behaviors of apo states and vesnarinone-bound hERG1 systems . 93 4.3.2.5 MM/PBSA (Molecular Mechanics/Poisson Boltzmann surface area) calculations ... 95
4.3.2.6 Extension of the MD simulations of hERG1-Vesnarinone complex systems ... 102
5. CONCLUSIONS AND RECOMMENDATIONS ... 111 REFERENCES ... 113 CURRICULUM VITAE ... 127
ABBREVIATIONS
B3LYP : Becke-Three Parameter Lee, Yang, Parr Density Functional cAMP : 3′, 5′-cyclic adenosine monophosphate
CG : Conjugate Gradient
cGMP : 3′, 5′-cyclic guanosine monophosphate CS : Conformational Search
DFT : Density Functional Theory ED : Erectyle dysfunction
EMEA : European Medicines Agency
FDA : U.S Food And Drug Administration
GBVI/WSA dG : Generalized Born Volume Integral/Weighted Surface Area hERG : human ether-à-go-go-related gene
LGA : Lamarckian Genetic Algorithm LINCS : Linear Constraint Method LQTS : Long QT Syndrome MD : Molecular Dynamics
MM-GBSA : Molecular Mechanics Generalized Born Solvation Area MM-PBSA : Molecular Mechanics Poisson-Boltzmann Solvation Area NPT : Isothermal-Isobaric Ensemble
NVT : Isochoric-Isothermal Ensemble
PDE11 : PDE11-phosphodiesterase 5-type enzyme PDE5 : PDE5-phosphodiesterase 5-type enzyme PDE6 : PDE6-phosphodiesterase 5-type enzyme PME : Particle Mesh Ewald
QSAR :Quantitative Structure-Activity Relationships RMSD : Root mean square deviation
RMSF : Root mean square fluctuation SD : Steepest Descent
SAR : Structure-Activity Relationships TPSA : Topological Surface Area
SYMBOLS
Å : Angstrom
Cmax : Maximum Concentration fs : femto second
IC50 : Inhibitor concentration-half maximum MW : Molecular Weight nm : nanometer ns : nano second μM : micro molar μs : micro second ps : pico second K : Kelvin T : Temperature
LIST OF TABLES
Page Table 2.1 : Comparison of different docking programs: GOLD, MOE and AutoDock (“os” and “ois” stands for open-state and open-inactivated-state of hERG K+ ion channel, respectively. Binding Energies are given in kcal/mol). 18 Table 2.2 : Summary of interactions of sildenafil with the hERG channel. (Top 10 docking poses of GOLD program are considered in generating the specific interactions, Gray Area: Open-State, White Area: hERG-Open-Inactivated-State) ... 22 Table 2.3 : In silico Alanine Mutagenesis Results. The binding free energies (ΔGbinding) are expressed in kcal/mol. The results are calculated according
to the top-ranking pose generated with GOLD program. ... 26 Table 2.4 : 2D structures and predicted binding free energies of the novel sildenafil analogs. ... 31 Table 2.5 : Physicochemical properties of the molecules, calculated by MOE. ... 34 Table 2.6 : Comparison of protein-ligand free energy results of Sildenafil and Frag_21 Molecules Using MM/GBSA Analyses. ... 34 Table 3.1 : Predicted binding free energies (Chemscore.dG) and 2D structures of ZINC compounds against the principal target, PDE5 and off-target enzymes, PDE6 and PDE11 (binding scores are expressed in kJ/mol and calculated by Chemscore fitness function implemented in GOLD Docking Program). Chemscore.dG values were converted to calculated IC50 values -for the purpose of selectivity comparison-according to the
formula; ΔGbinding =RTlnIC50, where T is taken as 300 K. ... 48 Table 3.2 : Comparison of protein-ligand free energy results of tadalafil with selected hit compounds using MM/GBSA calculations. ... 63 Table 3.3 : Predicted binding affinities of the selected compounds within the hERG K+ channel. Each compound was docked into the central cavities of the channel by GOLD docking software with ChemScore fitness function. Docking studies were realized by considering the two known conformational states of the channel. OS and OIS states stand for the open and open-inactivated states. dG.ChemScore values are expressed in kJ/mol. Tadalafil and the two selected potent and selective PDE5 inhibitor compounds ZINC02120502 and ZINC16031243 are shown in bold in the Table. ... 66 Table 4.1 : Docking scores (Fitness) of vesnarinone binding inside the central cavities hERG channels, generated with GOLD program. ... 83 Table 4.2 : Average Lennard-Jones (LJ), Average Coulomb (CL) and the total value of the averages of these energies between the hERG1 channel and vesnarinone, extracted from the whole simulation time (50 ns). ASP, ChemScore, CHEMPLP and GoldScore denote the fitness functions that served as the docking algorithms to deliver the starting input geometries for the MD simulations. The energies (expressed in kJ/mol unit), given in bold, represent the lowest values of their corresponding columns. .... 87
Table 4.3 : MM-PBSA Results of Vesnarinone Binding to hERG Channels. Interaction Binding Energy (ΔGbind) components with their standart
deviations are shown in kJ/mol. ... 99 Table 4.4 : MM-PBSA Results of Vesnarinone Binding to hERG1 Channels (for GoldScore trajectories) during the 0.5 μs simulation time. Interaction binding energies with different components with their standart deviations are shown in kJ/mol. ... 106 Table 4.5 : Comparison of docking poses and scores using model and cryo-EM structures. ... 108
LIST OF FIGURES
Page Figure 2.1 : 2D chemical structures of sildenafil, vardenafil and tadalafil and their corresponding IC50 values for the hERG K+ ion channel. ... 6 Figure 2.2 : Flowchart of the current study in designing the novel sildenafil-like molecules based on fragment replacement strategy in order to obtain promising molecules whose binding affinity to hERG K+ ion channel is decreased as the pharmacological activity against their therapeutic receptors, PDE5, is retained... 9 Figure 2.3 : Virtual screening analysis. Binding energy values are plotted against the frequency (number of molecules in the corresponding energy range). Values are obtained at the end of the VS protocol via docking the sildenafil analogs against hERG channel in its open state by means of MOE software. ... 10 Figure 2.4 : 3D structures of homology models of hERG K+ ion channels used in the current study. S1-S6 and S5-S6 domains are shown for the open state (a and b, left side) and open-inactivated state of the channel (a and b, right side), respectively. ... 12 Figure 2.5 : (A) 3D ligand interactions diagram of the sildenafil generated with GOLD program using rigid target treatment. Bidentate hydrogen bondings interactions with the invariant Gln817 residue are shown with red-dashed lines. (B) Superposition of sildenafil orientations obtained with different docking programs is shown: White color represents the rigid receptor treatment; turquaz color represents the flexible receptor treatment; brown color represents the flexible residue treatment on the H-loop domain, obtained with GOLD); green color represents the MOE program binding pose prediction and yellow color represents the AUTODOCK program binding pose prediction. The coloring scheme is based on carbon atoms. (C) The group definition of sildenafil, used in the current work. ... 14 Figure 2.6 : 2D protein-ligand interaction diagrams of top docking poses produced by three different programs are represented. Polar and hydrophobic resiues are circled in purple and green colors, respectively. Basic residues are shown with purple color with an exterior blue circle. Green-dashed arrow represents hydrogen bonding interaction between ligand and the side chain of a residue. Solvent-accessible surface area for ligand and receptor is shown with blue smudge and turquoise halo, respectively. Arene-arene and arene-hydrogen bonding interactions are shown with green colors. 3D protein–ligand interaction diagrams are also shown. .. 18 Figure 2.7 : Superpositions of the top ten docking solutions for sildenafil orientations at the central cavity of hERG1 (OS), generated by AutoDock and their corresponding binding energy predictions. ... 23 Figure 2.8 : A schematic view of the top docking poses of vardenafil and tadalafil at the central cavity of hERG1 channel. Ribbon representations and side
chains of the channel are covered within the 2.5 Å space around the ligand. Hydrogen bonds are shown with red dotted line with their distances. Left and right columns represent the interactions with the open state and open-inactivated state of the hERG1 channel, respectively. ... 29 Figure 2.9 : Comparison of RMSD graphs of sildenafil and sildenafil_frag21 molecules at the PDE5 and hERG1 open and open inactivated state targets. ... 34 Figure 3.1 : Flowchart of the current study in the effort for identifying novel and selective PDE5 inhibitors. ... 40 Figure 3.2 : Electrostatic maps of the active sites of the enzymes (left panel). Blue, red and white colours represent positively charged, negatively charged and hydrophobic preferences built at the drug-binding cavity site (near 5 Å distance from the ligand). Tadalafil fulfills the positively charged electrostatic requirement created by the acceptor atom (Oɛ) of the invariant Glutamine side chain in each active sites via carrying a hydrogen-bond donating moiety (-NH) at the amide fragment (namely, Glutamine Switch). On the other hand, hydrophobic residues, Phe820, Phe776 and Trp 820, sandwich the ligand (namely, Hydrophobic Clamp). 2D ligand-protein interaction diagrams (right panel). Green arrows indicate hydrogen bonding interactions. Green and purple discs show hydrophobic and polar residues, respectively. Representations are created with MOE molecular modeling package. ... 45 Figure 3.3 : Top docking poses of tadalafil with PDE enzymes. Only polar hydrogens are shown for clarity. Protein residues within 2.5 Å distance around tadalafil are depicted in the figures. Hydrogen bonds between amide hydrogen of tadalafil and Oɛ atom of invariant Glutamine amino acid residue are represented with red dashed lines. ... 46 Figure 3.4 : Aminoacid sequence alignments of PDE6 andPDE11 over PDE5 catalytic site residues. Alignment procedure was achieved with BLOSUM62 matrix via MOE software. ... 46 Figure 3.5 : Ramachandran plots of the homology models of PDE6 and PDE11. ... 47 Figure 3.6 : Superposition of PDE5, PDE6 and PDE11, illustrated with blue, cyan and red colors, respectively. The counterions, Zn2+ and Mg2+, are shown with blue circles at the metal binding side. ... 53 Figure 3.7 : Contact energy profiles of the catalytic side residues that correspond to PDE5, PDE6 and PDE11. The x axis and y axis represent the aminoacid residues numbering and atom-atom contact pair energies in kcal/mol unit, respectively. ... 54 Figure 3.8 : 2D protein-ligand interaction diagrams of the selected compounds (Table 3.1) with the catalytic site residues of PDE5. ... 55 Figure 3.9 : Superposition of 27 selected compounds at the end of the docking simulations. ... 56 Figure 3.10 : (A) Docked pose of ZINC02120502 at the active site of PDE6. (B) Docked pose of ZINC02120502 at the active site of PDE11. (C) Overlay of ZINC02120502 (tan) onto the crystal orientation [118] of tadalafil (blue). (D) Overlay of ZINC16031243 docked poses at the active site of PDE5 (blue), PDE6 (pink) and PDE11 (tan). Protein–ligand interaction diagrams of ZINC02120502 and ZINC16031243 with PDEs (shown in the bottom of the figure). ... 57
Figure 3.11 : Traces of protein backbone RMSD (root-mean-squared-deviation) evaluation during the whole production stages of the MD Simulations. 60 Figure 3.12 : The RMSD evaluation of the ligands during the simulation time. ... 61 Figure 3.13 : Traces of hydrogen bonding interactions throughout simulation time (x and y axis represents the simulation time and distances between Oε atoms of Gln817 (PDE5), Gln773 (PDE6), Gln869 (PDE11) and indole fragments hydrogen in the ligands, respectively). Color codes: orange: PDE5+ZINC16031243; purple PDE11+tadalafil; blue: PDE5+ZINC02120502; red: PDE5+tadalafil; green: PDE6+tadalafil. .. 62 Figure 3.14 : Simulated structures of ZINC0210502 and ZINC16031243 at the substrate pockets of PDE5, PDE6 and PDE11. ... 64 Figure 3.15 : Overlay of docking pose (blue) and representative structure of ZINC02120502 (white) in the catalytic pocket of PDE11. ... 65 Figure 3.16 : Root-mean-square fluctuation (RMSF) values per residue during the MD simulations. ... 67 Figure 3.17 : (Top) Derived top-scored six-sited (RRRHHH) E-pharmacophore model; (bottom) 176 000 compounds from Otava small-molecules database are screened against derived pharmacophore model and top-1000 compounds that have high Fitness scores with these sites are then docked at the PDE5 binding pocket using Glide/SP (standard precision). Compounds that show high docking scores as well as high fitness scores are shown in the figure. 2D ligand interaction diagram of selected Otava compound (1094821) is also represented in the figure. ... 69 Figure 4.1 : Topology of the hERG channel in its open-state conformation. The channel is a coassembled of four identical chains; each chain is shown with different colours (top). S5 and S6 helices line the pore domain of the channel. A close look at the some selected S6 and SF (selectivity filter) residues that are important for drug binding (bottom). ... 73 Figure 4.2 : Open state hERG Channel (S5 and S6 chains, only, represented with green ribbon) inserted in DPPC lipid-bilayer. Water and lipid molecules are illustrated with red and olive-orange lines. Potassium cations (K+) are located at the energetically favorable S0-S2-S4 sites at the selectivity filter (SF) and also pore domain (PD) region, depicted in pink spheres. 75 Figure 4.3 : The superposition of the top docking poses of vesnarinone, generated by four different scoring functions in the GOLD program... 79 Figure 4.4 : 2D and 3D depiction of the interaction of vesnarinone with the aminoacid residues that line the pore domain and central inner cavity of the hERG1 channels. HOIS and HOS stands for the open-inactivated and open states of the hERG1 channel, respectively. The red, purple and green lines represent the hydrogen bonding interactions and π-π stacking interactions within the residues and vesnarinone at around 5 Å distance, respectively. Green, turquoise and gray colours indicate hydrophobic, polar and Glycine residues, respectively. ... 80 Figure 4.5 : Backbone RMSD values of hERG1 open (top) and open-inactivated (bottom) states during the MD simulations. ... 83 Figure 4.6 : Ligand RMSD traces during the production stages of the MD simulations (x and y axes represent the simulation time and ligand RMSD values in ns and Å units, respectively). (The least-square fit was done on protein Cα atoms for RMSD measurements.). ... 84
Figure 4.7 : Ligand RMSD traces during the production stages of the MD simulations. X and y axes represent the simulation time and ligand RMSD values in ns and Å units, respectively. The minimum, maximum and average RMSD values are also shown on the graphs, in Å unit. (The least-square fit was done based on the starting position of vesnarinone at the begining of simulation in order to measure the conformational flexibility of the ligand.) ... 85 Figure 4.8 : Short range energetics (van der Waals and Electrostatics) interactions between vesnarinone and hERG1 channel throughout the simulation time. Left and right panels represent the open and open-inactivated states of the channel, respectively; starting from the docking outputs of the ASP, ChemScore, CHEMPLP and GoldScore fitness functions (from top to down). CL and LJ are abbreviations for Coulomb (Electrostatics) and Lennard Jones (van der Waals), respectively. ... 86 Figure 4.9 : The crucial residue-vesnarinone interactions of the initial geometry and MD-simulated structure of the GoldScore trajectory for the open-inactivated state of hERG1. Hydrogen bonding and π-π stacking interactions are depicted in the figure.The distances of stacking interactions are measured based on the centroids of the aromatic rings. 87 Figure 4.10 : Distances of remarkable intermolecular interactions between vesnarinone and hERG1 channels residues in its open-inactivated state, during the MD simulations. “C=O group’s oxygen atom” and “quinolin moiety’s oxygen atom” phrases refer to the atom groups in vesnarinone. ... 88 Figure 4.11 : Distances of remarkable intermolecular interactions between vesnarinone and hERG channels residues in its open-state during the MD simulations. “Methoxy oxygens and quinolin moiety –NH group’s H atom” phrases refer to the atom groups in vesnarinone. ... 90 Figure 4.12 : The crucial residue-vesnarinone interactions of the starting geometry and MD-simulated structure of the GoldScore trajectory for the open state of hERG1.Hydrogen bonding and π-π stacking interactions are depicted in the figure. The distances of stacking interactions are measured based on the centroids of the aromatic rings. ... 91 Figure 4.13 : Snapshots from MD simulations (namely, ChemScore trajectory) of hERG1open state-vesnarinone complexed system. Vesnarinone and hERG1 are represented with sticks and ribbon presentations, respectively. Each chain is coloured with different colours and one of the chains of hERG is deleted, for clarity (top). Overlay of the snapshots from the simulations were shown. The representation was created with a color code, from red to blue according to the simulation time and only 100 frames (taken each 50 frames in the evenly distributed time course) are shown for clarity (bottom). ... 92 Figure 4.14 : Overlay of the snapshots from the MD simulations (namely, ChemScore trajectory) of hERG1 Open-inactivated state-vesnarinone complexed systems.Vesnarinone and hERG1 are represented with sticks and ribbon presentations, respectively. Each chain is coloured with different colours and two of the chains of hERG are deleted, for clarity (top). Snapshots from the simulations were shown. The representation was created with a color code, from red to blue according to the
simulation time and only 100 frames (taken each 50 frames in the evenly distributed time course) are shown for clarity. (bottom). ... 93 Figure 4.15 : Distance between hydrogen atom of Ser649 and dihydroquinolin moiety’s oxygen atom of vesnarinone, in the ChemScore trajectory for the open state of hERG1 channel during the 50 ns simulation production period. Just at the beginning periods of the simulation (within the 0-5 ns time period), H-bond between these atoms was broken and never formed later for the rest of the simulation. ... 94 Figure 4.16 : First ten eigenvalues of the protein backbone movement along the simulation time. Note that, slopes of the ChemPLP and GoldScore curves overlap for the hERG1 open inactivated state. ... 95 Figure 4.17 : Per-residue energy decompositions of MM-PBSA binding energy calculations. (Contributions of each residue were presented based on the cumulative contributions from each chain.). ... 100 Figure 4.18 : The course of distances established between methoxy oxygens of vesnarinone and Ser649 H atom. The positions of Tyr652 and Phe656 residues at the beginning and end of the MD simulation. Vesnarinone is shown in magenta color. ... 101 Figure 4.19 : Vesnarinone superposition for representative frames of GoldScore (green) and ChemPLP (red) trajectories in the hERG1 open inactivated state. Only polar hydrogen is shown for clarity. ... 102 Figure 4.20 : Protein backbone RMSD traces during the 550-ns production stages of the MD simulations, started from the GoldScore docking output geometries. The overlay of vesnarinone at 50-ns and 550-ns are also depicted (below). ... 103 Figure 4.21 : Ligand RMSD traces during the 0.55 μs production stages of the MD simulations, started from the GoldScore docking output geometries. .. 104 Figure 4.22 : Short range energetics (van der Waals and Electrostatics) interactions between vesnarinone and hERG1 channel throughout the 0.55 μs simulation time (for GoldScore trajectories). CL and LJ represent the Coloumb (Electrostatics ) and Lennard Jones interactions (van der Waals), respectively. Average Coloumb and Lennard Jones interaction energies are -44.37 kJ/mol and -195.44 kJ/mol for the open state of hERG1 whereas average Coloumb and Lennard Jones interaction energies are -91.50 kJ/mol and -233.76 kJ/mol for the open inactivated state of hERG1. ... 105 Figure 4.23 : Per-residue interaction energies of MM-PBSA binding energy calculations of vesnarinone binding at the hERG1 channels during the 0.5 μs simulation time. (Contributions of each residue were presented based on the cumulative contributions from each chain.). ... 106 Figure 4.24 : Alignment of hERG channel cryo-EM (cyan) and 3D model
(open-state) (orange) structures... 107 Figure 4.25 : Per-residue energy decompositions of MM-PBSA binding energy calculations of vesnarinone binding to open-state hERG1 channels (using Cryo-EM structure) during the 50 ns simulation time. (Energetic contributions of each residue for binding were presented based on the cumulative contributions from each chain.). (Residues between 635-668 are shown only, for clarity.) Note that, the calculated per-residue contribution energy profile demonstrates merely the preliminary outcome for 50 ns simulation time. The relatively high positive energetic
contribution to binding enthalpy–that came from polar Ser624 residues-need to be better refined via longer simulation time which may eventually further rectify the polar solvation energy component of ΔGbinding that also militate against binding process. For all that, with respect to the aminoacid residues that are situated at the pore domains of
the channel, the relevant outcome still presents a qualitative data whether which residues tend to mostly contribute to binding; such as the bulky and aromatic Tyr652 residues on the S6 helix of the channel. ... 109
IN SILICO DESIGN OF HERG NON-BLOCKER COMPOUNDS WITH RETAINED PHARMACOLOGICAL ACTIVITY USING MULTI-SCALE
MOLECULAR MODELING APPLICATIONS SUMMARY
hERG K+ ion channels are responsible in the regulation process of the action potential of human ventricular myocyte by contributing the rapid component of delayed rectifier K+ current (IKr) component of the cardiac action potential. The direct
inhibition of hERG channels/current raises cardiac diseases, i.e., serious life-threatening arrhythmias further leading to Torsades de pointes (TdP) and long QT syndrome (LQTS). Various drugs have been withdrawn from the drug markets such as terfenadine, cisapride, astemizol or restricted in their use icluding thioridazine, haloperidol, sertindole, and pimozide due to hERG-related life-threatening cardiac arrhythmias.
In this thesis, molecular interactions between human ether-à-go-go-related gene (hERG) blocker compounds and pore domains of hERG1 K+ ion channels (in its
voltage driven open and open-inactivated conformational states) together with their
target proteins were investigated using computational modeling methods. In silico drug rehabilitation studies were conducted for the selected compounds (i.e., some phosphodiesterase 5-type enzyme (PDE5) inhibitors) as well with an aim of reducing their hERG1 blocking affinity while keeping their principal enzymes’ activity. Using computational approaches such as homology modeling, protein-ligand docking, Alanine mutagenesis studies, E-pharmacophore modeling, structure- and ligand-based virtual screening strategies (such as rapid screening of drug-like molecular database with molecular similarity approaches and fragment based in silico design of novel compounds using fragment libraries), molecular dynamics (MD) simulations, post-MD analyses via considering the detailed analysis of the trajectories-especially conformational behaviours of both the apo and holo systems-and post-MD computations ((e.g., Molecular Mechanics Poisson-Boltzmann Surface Area and Molecular Mechanics Generalized Born Surface Area (MM-PBSA and MM-GBSA, respectively)) for evaluating the binding energies of the studied ligands towards the hERG channels and also towards PDE enzymes have enabled to get insight into the important blocking elements of hERG channels in terms of drug association with its central inner cavities. Additionally, some drug-like hits (with reduced hERG blocking affinity and kept principal target activity) were presented in the consequence of the multi-step virtual screening protocols.
Hence, the consequences of the derived outcomes from this thesis may be helpful in the effort for designing novel and safe drugs.
The studied drugs, herein, (e.g., sildenafil-ViagraTM,vardenafil-LevitraTM, tadalafil-CialisTM and vesnarinone) are all PDE5 inhibitors used in the treatment of erectyle dysfunction except for vesnarinone which is used as a PDE3 inhibitor). Phosphodiesterase (PDE) type enzymes are important biological elements of human organismsand found in many tissues and organs. Their biological function in human
organism is to regulate the cytoplasmic concentration levels of intracellular second messengers, 3’,5’-cyclic guanosine monophosphate (cGMP) and/or 3’,5’-cyclic adenosine monophosphate (cAMP) via catalytic degradation reactions occuring in their substrate pockets. At the end of the applied multi-scale modeling applications, tadalafil and sildenafil-like molecules having drug-like properties (tested using pharmacokinetic modeling applications) have been identified and proposed as novel and safe potential PDE5 inhibitors. In addition to the hERG direct blocking investigations, selectivity of tadalafil-like hits has also been tested against PDE6 and PDE11 isoenzymes which have been identified as the causes of the side-effects of the used PDE5 inhibitor drugs in the market.
In summary, through the duration of this thesis, various computational modeling tools have been utilized and applied in the investigation of novel and safe drugs (with reduced side-effects towards hERG channels and also towards the related PDE-isoenzymes) whilst keeping their principal target activities. Besides, the critical protein-ligand interactions have been illuminated at the molecular level along with the dynamics of the studied drugs and their target proteins.
This thesis comprises of three chapters that present three scientific articles that have been published at peer-reviewed journals.
In the first article (Chapter 2), various computer aided drug-designing and computational molecular modeling techniques were used for the investigation of the action mechanisms of the FDA approved drugs, sildenafil-ViagraTM, vardenafil-LevitraTM and tadalafil-CialisTM with their target protein, PDE5 and also hERG1 channels. In addition, fragment-based virtual screening strategy was employed in order to obtain potent and safe sildenafil-like molecules.
Second article (Chapter 3) focuses on the development of potent tadalafil-like molecules via combination of ligand-based screening and structure-based modeling protocols. hERG1 binding affinities of the selected compounds together with tadalafil were evaluated via flexible molecular docking computations. In addition, PDE5/PDE6 & PDE5/PDE11 selectivities of the compounds were studied and important structural binding patterns such as the critical residue-ligand interactions were highlighted.
Third article (Chapter 4) deals with vesnarinone (used a PDE3 inhibitor agent)-hERG complex systems and time-dependent dynamical behaviour of vesnarinone at the pore domains of hERG channels in its open and open-inactivated states. By the use of molecular docking, MD simulations and detailed post-MD analysis computations, possible binding modes of vesnarinone within the central cavities of the channels were proposed. Also, crucial hERG residues in terms of vesnarinone binding were further highlighted which may help to design safe and novel drug-like molecules. In Chapter 5, overall interpretation of the results and potential further works are briefly presented.
HERG BLOKER OLMAYAN FARMAKOLOJİK AKTİVİTESİ KORUNMUŞ BİLEŞİKLERİN ÇOK BOYUTLU MOLEKÜLER MODELLEME
UYGULAMALARI İLE İN SİLİKO TASARIMI ÖZET
Bu tezde, bilgisayar destekli moleküler modelleme metodları kullanılarak, hERG blokörü olan bileşiklerin hERG1 K+
iyon kanalları (farklı membran potansiyellerinde açık ve açık-inaktif konformasyonel halleri göz önünde bulundurularak) ve ilgili hedef proteinleri ile olan etkileşimleri moleküler seviyede detaylı olarak incelenmiştir. İn siliko yöntemlerden yararlanılarak, literatürden seçilmiş bazı fosfodiesteraz tip 5 (PDE5) enzim inhibitörlerine yönelik olarak hERG1 kanalları için bağlanma afinitelerinin azaltılması amacı ile rehabilitasyon çalışmaları yapılmıştır. Aynı zamanda, bu bileşiklerin hedef proteinleri ile olan etkileşimlerinin aydınlatılması ve bağlanma enerjilerinin korunmasına yönelik olarak da hesaplamalar yapılmış ve incelemelerin sonuçları tezde sunulmuştur. Protein-ligant moleküler kenetlenme, homoloji modelleme, Alanin mutasyon çalışması, E-Farmakofor modelleme, yapısal ve ligant bazlı olarak sanal molekül kütüphanelerinin çeşitli teknikler ile tarama stratejileri, moleküler dinamik (MD) simülasyonlar, MD sonrası analizler (simülasyon trajektörilerinin detaylı olarak incelenmesi-özellikle apo ve holo sistemlerinin konformasyonal davranışları olmak üzere) ve ligantların kompleks yaptığı proteinlere (hERG1 iyon kanalı ve aynı zamanda PDE enzimlerine karşı olan) bağlanma enerjilerinin tayin edilmesi amacı ile MD sonrası yapılan hesaplamalar ((Moleküler Mekanik Poisson Boltzmann Yüzey Alanı PBSA) ve Moleküler Mekanik Genelleştirilmiş Born Yüzey Alanı (MM-GBSA)) neticesinde, hERG kanallarının merkez iç bağlanma bölgesinde gerçekleşen ilaç-benzeri küçük moleküller ile etkileşimlerinde rol oynayan önemli bloke edici elementler ve aminoasit rezidüleri ortaya konulmuştur. Aynı zamanda, uygulanan çok kademeli sanal tarama protokolleri sonucunda, bazı ilaç-benzeri moleküller (ana reseptör bağlanma afinitesi korunarak hERG bağlanma afinitesi düşürülmüş) sunulmuştur.
hERG K+ iyon kanalları, insan kalp karıncığı kas hücrelerinin görev aldığı aksiyon potansiyelinin düzenlenmesi işleminde rol üstlenirler. Bu aşamada, hERG kanalları, hücre içi-dışı iyon konsantrasyon gradientlerinden kaynaklanan hücresel elektriksel potansiyeline bağlı olarak genel olarak üç aşamalı bir konformasyonel aşamadan (açık hal, açık-inaktif hal ve kapalı hal olmak üzere) geçerek potasyum akışını sağlarlar. hERG kanallarının içerisinden gerçekleşen potasyum iyonları akışı, hERG kanallarının sahip olduğu seçici filtrenin (SF) sahip olduğu aminoasit rezidülarının (S624, V625,G626, F627, G628) K+ iyonları ve su molekülleri ile oluşturabildiği koordinasyon motifi sayesinde gerçekleşir. Bu potasyum akışı, kalp aksiyon potansiyelinin, IKr komponentine (geciktirilmiş düzeltici K+ akışı hızlı bileşeni)
önemli bir katkı yapar. hERG kanallarının ve görev aldıkları akışın doğrudan inhibe edilmesi, ciddi ve ölümcül kalp hastalıklarına (örneğin Torsades de pointes (TdP) ve
uzun QT sendromu (LQTS)) neden olabilmektedir. Yıllar içerisinde, yol açtıkları kardiyak aritmiyaları sebebi ile birçok ilaç (bazı antihistaminler, antibakteriyeller, antipsikotikler ve antidepresanlar dahil olmak üzere) piyasadan çekilmiştir ve/veya kullanımları kısıtlandırılmıştır. Bu ilaçlara örnek olarak terfenadine, cisapride, astemizol, thioridazine, haloperidol, sertindole ve pimozide verilebilir. Bu yüzden, FDA (Amerikan Gıda ve İlaç Dairesi) yeni ilaç geliştirme aşamalarında yapılan testlerde ilaçların hERG taramasından geçirilmesini 2007 yılında zorunlu kılmıştır. Tezde, seçilen ilaçlara uygulanan in siliko rehabilitasyon çalışmalarının sonuçları, bu sebeple yeni ve güvenli ilaç tasarım çalışmalarına bir katkı olarak değerlendirilebilir.
Bu tezde çalışılan ilaçlar (sildenafil-ViagraTM
, vardenafil-LevitraTM, tadalafil-CialisTM ve vesnarinone)-vesnarinone haricinde (PDE3 inhibitörü)-PDE5 inhibitörüdürler ve erektil fonksiyon bozukluğunu tedavi amacı ile kullanılırlar. PDE enzimleri, insan organizmasında çok geniş doku ve organlara yayılmış ve birçok hastalığın da tedavisinde kullanılan ilaçların etki ettiği ana reseptör olma özelliğine sahiptirler. Bu enzimlerin biyolojik fonksiyonları, substrat (3’,5’-siklik guanozin monofosfat (cGMP) ve/veya 3’,5’-siklik adenozin monofosfat (cAMP)) bağlanma bölgelerinde gerçekleştirdikleri katalitik degradasyon reaksiyonları üzerinden bu substratların sitoplazmik konsantrasyon seviyelerini düzenlemektir. Bu tezde uygulanan çok kademeli moleküler modelleme teknikleri neticesinde, ilaç olma özelliğine sahip (farmakokinetik özellikleri moleküler modelleme uygulamaları ile tahmin edildikten sonra) yeni ve güvenilir (hERG1 ve PDE5 bağlanma enerjileri moleküler kenetlenme ve MM-PBSA ve MM-GBSA hesaplamaları ile kontrol edilmiştir) PDE5 inhibitörleri olarak tadalafil ve sildenafil benzeri moleküller ortaya konulmuştur. İlaçların, ana reseptörlerine karşı olan bağlanma afinitelerinin korunarak hERG bloke etme özelliklerinin düşürülmesi çalışmalarının yanısıra, PDE6 ve PDE11 enzimlerine karşı olan aktiviteleri de çalışılmıştır. Çünkü, PDE6 ve PDE11 enzimlerinin katalitik kısımlarının PDE5 ile aminoasit sekansı ve 3-boyutlu sekonder yapılarındaki benzerlikleri ve bu benzerliklerin neticesinde ise PDE5 inhibitörü olarak kullanılan ilaçların bu izoenzimlere karşı aktivite gösterdiği bilinmektedir. Bu izoenzimlere karşı olan bağlanma yatkınlığı neticesinde bazı yan etkiler (görme bozukluğu, kas ağrıları, vb.) görülebilmektedir.
Özetle, tez süresince, birçok bilgisayar destekli modelleme araçları ile birlikte moleküler modelleme teknik ve yaklaşımlarından yararlanılmış, yeni ve yan etkileri azaltılmış (ana reseptöre olan-PDE5-etki mekanizması korunarak) güvenilir ilaç-benzeri moleküller ortaya konulmuştur. Bununla birlikte, önemli ve kritik protein-ligant etkileşimleri detaylı bir şekilde moleküler düzeyde çalışılmış ve tartışılmıştır. Moleküllerin ve etki ettikleri proteinlerin dinamik özellikleri (zaman içerisindeki konformasyonel değişimleri vb.) uygulanan MD simülasyon teknikleri ve MD sonrası trajektörü analizleri ile ortaya çıkarılmıştır.
Bu tez, SCI ve SCI-Expanded indeksinde yeralan dergilerde yayını gerçekleştirilmiş üç adet makalenin derlenmesinden oluşturulmuştur. Yayınlanma tarihlerine göre makalelerin içerikleri ve amaçları şu şekilde detaylandırılabilir:
Birinci makalenin konusu; sildenafil, vardenafil ve tadalafil ilaçlarının (FDA onaylı birinci jenerasyon PDE5 inhibitörleri), hedef reseptörleri olan PDE5 enzimi ve hERG1 kanal modelleri (açık ve açık-inaktif halleri) ile etkileşimlerinin birden fazla moleküler kenetlenme programı (esnek, rijit ve indüklenmiş protein-ligand kenetlenme yaklaşımları uygulanarak) kullanılarak incelenmesi, bağlanma özellikleri
ve bağlanma enerjilerinin ortaya konulması olarak özetlenebilir. Çalışmanın amacı ise, sildenafil molekülünün fragment-bazlı sanal tarama tekniği kullanılarak belirli bir fragmentin çıkarılıp bu fragment yerine eklenen fragmentler ile oluşturulan bir moleküler sanal kütüphanenin moleküler kenetlenme hesapları ile taranması ve MD simülasyon ile MM-GBSA analizleri neticesinde hERG afinitesi düşürülmüş ve PDE5 aktivitesi korunmuş sildenafil-benzeri moleküllerin açığa çıkarılması şeklinde verilebilir. Ayrıca, hERG kanalının bağlanma bölgesinden seçilmiş bazı aminoasit rezidüleri için sildenafilin bağlanmasına yönelik olarak Alanin mutasyon çalışması da yapılmıştır. Bu aşamadaki amaç hangi amoasitlerin sildenafilin bağlanmasına daha çok katkı yaptığının aydınlatılmasıdır. Seçilen rezidüler teker teker Alanin aminoasit rezidüsüne dönüştürülmüş ve moleküler kenetlenme çalışmaları elde edilen bu mutant hERG kanalları ile tekrarlanmıştır. Ortaya çıkan bulgular (moleküler kenetlenme skorları) ile doğal-tipteki (wild) sonuçlar karşılaştırılmış ve sildenafilin bağlanması bağlamındaki kritik rezidüler açığa çıkarılmıştır.
İkinci makalede ise, diğer bir PDE5 inhibitörü olarak kullanılan tadalafil molekülü ele alınmıştır. Geniş bir ilaç-benzeri moleküler veri tabanına ligand-bazlı bir tarama tekniği uygulanarak (Tanimoto benzerlik katsayısı ile MACCS yapısal moleküler parmakizi araçlarının kombine edilmesi vasıtasıyla) tadalafile %80 oranında yapısal benzerlik gösteren moleküller filtrelenmiştir. Filtrenen bu moleküllerin, PDE5, PDE6, PDE11 ve hERG kanallarına esnek moleküler doking çalışmaları yürütülmüştür. Bundan önce, PDE6 ve PDE11 enzimlerinin deneysel üç boyutlu yapılarının henüz keşfedilememesi sebebiyle 3D homoloji modelleri yapılmıştır. Seçilen bazı moleküllerin MD simülasyonları gerçekleştirilmiştir. Buradaki temel amaç, filtrelenen tadalafil-benzeri moleküllerin PDE5/PDE6 ve PDE5/PDE11 seçiciliklerinin korunabilmesi ve ayrıca bazı hERG afinitesi düşürülmüş moleküllerin önerilmesidir. Bunun için, seçilen moleküller ile hERG kanallarının hem açık hem de açık-inaktif hallerine moleküler kenetlenme çalışmaları yapılmış ve sonuçlar yayınlanmıştır.
Üçüncü makalede ise, vesnarinon molekülünün (PDE3 inhibitörü) hERG kanallarının iç boşluk kısmına (açık ve açık-inaktif halleri gözönünde bulundurularak) bağlanması ile ilgili olarak, değişik skorlama fonksiyonları kullanılarak yapılan moleküler kenetlenme çalışmaları sonuçlarının incelenmesi, birbirleri ile karşılaştırılması ve peşinden gerçekleştirilen MD simülasyonların ve MM-PBSA analizlerin yardımı ile vesnarinon molekülünün hERG kanalları ile olan etkileşimlerinin analiz edilmesi ve en olası bağlanma oriyantasyonunun açığa çıkartılması hedeflenmiştir. Elde edilen moleküler kenetlenme skorları, MM-PBSA hesaplamalarından gelen bağlanma enerjileri, simülasyonlardan elde edilen protein-ligant arasındaki Lennard-Jones (van der Waals) ve Elektrostatik (Coloumb) enerjileri ile karşılaştırılmıştır. Ayrıca her bir simülasyon yörüngesi için yapılan simülasyon süresindeki detaylı ligand/protein RMSD (karekök ortalama sapma değeri) ve PCA (Temel Bileşenler Analizi) analizleri ile elde edilen bulgular detaylandırılmıştır. Buradan elde edilen sonuçlar, hERG-ilaç molekülü etkileşimlerinin daha iyi anlaşılabilmesi, bu bağlamda kritik önem taşıyan aminoasit rezidülerinin aydınlatılması (hem protein-ligant kenetlenme hem de MM-PBSA sonrası yapılan rezidü başına dekompozisyon hesaplamaları vasıtasıyla) ve bunların neticesinde elde edilen bulgular ile beraber uygulanan in siliko modelleme tekniklerinin, hERG bağlanması rehabilite edilmek istenen moleküllere genel anlamda bir örnek teşkil edebilmesi bakımından önem taşıması olarak özetlenebilir.
1. INTRODUCTION
The subject of this Ph.D. thesis is to investigate the molecular interactions between hERG blocker compounds and the central inner cavities of hERG K+ ion channels (in its open and open-inactivated conformational states which are voltage-dependent processes driven by the ion concentration/flux across the cell membrane bilayers) along with their corresponding main target (e.g., PDE5) and off-target proteins (e.g., PDE6 and PDE11) via various in silico techniques. Also, computer-aided drug design strategies have been employed with aim of identifying safe and novel drug-like molecules with reduced hERG blocking affinity and cross-reactivity with off-target PDE enzymes while keeping their principal off-target activities.
The main concept and purpose of the studies within this Ph.D thesis is about understanding and illuminating the blocking mechanisms and blocking elements of some selected PDE inhibitors (sildenafil,vardenafil, tadalafil and vesnarinone) towards hERG1 channel models and also within their main targets and further develop novel and safe drug-like hits by utilizing multi-scale molecular modeling applications.
Bringing into novel and safe drugs to drug market is a long time requiring process with its own rigour sides and expensive multi-step pre-clinical and clinical stages. Besides, thanks to the rapidly developing computational molecular modeling techniques, (such as virtual screening strategies of large molecular databases, protein-ligand docking, QSAR (Quantitaive-Structure Activity Relationships) developments and many other structure- and ligand-based methods in order to predict the accurate bioactive conformations and correct binding free energy trends of small molecules with their target proteins), cost-effective alternative routes could be supplemented to the early phases of the conventional drug-discovery (developing new molecules) and/or rehabilitation efforts of the undesired side-effects of existing drugs. Since 1990s, mandatory hERG toxicity screening tests have been incorporated into the new drug development steps by FDA (U.S Food and Drug Administration) and EMEA (European Medicines Agency). Computational screening techniques and
further in silico investigation of protein-ligand holo systems, as in this study, present a branch of useful scientific data in the area of a drug design field. Hence, the outcomes of this thesis which have been published in peer-reviewed journals may contribute to the design/rehabilitation efforts of hERG non-blocker compounds and also could be beneficial in the design of more potent and safe PDE enzymes.
This thesis comprises of three chapters that present three scientific articles that have been published at peer-reviewed journals.
In the first article (Chapter 2), various computer aided drug-designing and computational molecular modeling techniques were used for the investigation of the action mechanisms of the FDA approved drugs, sildenafil-ViagraTM, vardenafil-LevitraTM and tadalafil-CialisTM with their target protein PDE5, and also hERG1 channels. In addition, fragment-based virtual screening strategy was employed in order to obtain potent and safe sildenafil-like molecules.
Second article (Chapter 3) focuses on the development of potent tadalafil-like molecules via combination of ligand-based screening and structure-based modeling protocols. hERG1 binding affinities of the selected compounds together with tadalafil were evaluated via flexible molecular docking computations. In addition, PDE5/PDE6 & PDE5/PDE11 selectivities of the compounds were studied and important structural binding patterns such as the critical residue-ligand interactions were highlighted.
Third article (Chapter 4) deals with vesnarinone (used a PDE3 inhibitor agent)-hERG complex systems and time-dependent dynamical behaviour of vesnarinone at the pore domains of hERG channels in its open and open-inactivated states. By the use of molecular docking, MD simulations and detailed post-MD analysis computations, possible binding modes of vesnarinone within the central cavities of the channels were proposed. Also, crucial hERG residues in terms of vesnarinone binding were further highlighted which may help to design safe and novel drug-like molecules. In Chapter 5, overall interpretation of the results and potential further works are briefly presented.
2. IN SILICO DESIGN OF NOVEL HERG-NEUTRAL SILDENAFIL LIKE PDE5 INHIBITORS1
2.1 Introduction
Cyclic nucleotide phosphodiesterase enzymes (PDEs) have functions in regulating the levels of intracellular second messengers, 3′,5′-cyclic adenosine monophosphate (cAMP) and 3′,5′-cyclic guanosine monophosphate (cGMP), via hydrolysis and decomposing mechanisms in cells [1]. These enzymes are abundant in diverse tissues and systems such as immune and cardiovascular systems in human body. They take essential roles in modulating various cellular activities such as memory and smooth muscle functions [2-4]. The diversity and abundancy of these enzymes in human body render them potential drug targets for many diseases (i.e. asthma, depression, erectile dysfunction, pulmonary disease) [5-7]. PDEs are classified into 11 superfamilies according to their biochemical functions, structural and kinetic profiles, or their target substrates. These superfamilies have more than 60 subtypes and isoforms that are generated by different promoters in distinct tissues. Specifically, while cGMP is decomposed by PDE5, PDE6, and PDE9; and cAMP is hydrolyzed by PDE4, PDE7, and PDE8; PDE1, PDE2, PDE3, and PDE10 hydrolyze both cAMP and cGMP. PDE5A is the only subtype of PDE5 which has been identified so far. It has four isoforms (PDE5A1–4) that are mainly expressed in the corpus cavernosum [8]. The catalytic domain of PDE5A (residues: 535–860) complexed with sildenafil, vardenafil, and tadalafil – the drugs that are used in the treatment of erectile dysfunction (ED) was first isolated in 2002 by Sung et al. [8].
In silico studies have been used in better understanding of the molecular
determinants of drug-PDE5 binding interactions so far with a profound attention to the side effects of sildenafil (i.e. visual problems, headache, etc.) which stemmed from the cross-reactivity of this drug to the PDE6 and PDE11 enzymes, particularly.
1 This chapter is based on the article “Kayık, G., Tüzün, N. Ş., Durdagi, S. (2016). In silico design of hERG-neutral sildenafil-like PDE5 inhibitors. Journal of Bimolecular Structure and Dynamics, doi: 10.1080/07391102.2016.1231634.
Hence, enormous effort for developing “second generation PDE5 inhibitors” has been mainly focused on discovery of novel drugs with better selectivity against other enzymes within the PDE family, so far. For example, Huang et al. [9] performed homology modeling, molecular docking, and molecular dynamics (MD) simulations studies for understanding the better selectivity of tadalafil targeting to the PDE6 enzyme compared to sildenafil. Cichero et al. [10] worked on tadalafil analogs in order to gain a detailed perspective for understanding the key interactions between the catalytic side residues of PDE5 and PDE11 enzymes that are responsible for the cross-reactivity within these referred receptors. In addition to the structure-based design studies, ligand-based approaches were also applied for designing novel PDE5 inhibitors [11]. Many in vitro studies have also been carried out for the same aim. [12-20].
In addition to the aforementioned sanitary concerns, another unpleasant side effect of sildenafil that has been related to cardiovascular problems has also been considered. Patch clamp experiments showed that the IC50 values of the human
ether-à-go-go-related gene (hERG1) potassium (K) ion channel blocking affinity of sildenafil, vardenafil, and tadalafil as 33, 12, and 100 μM, respectively [21]. In vitro studies highlight the concentration-dependent blockage of hERG channel by these drugs [21,22]. hERG channel is responsible for the regulation of the action potential of human ventricular myocyte by contributing the rapid component of delayed rectifier K+ current (IKr) component of the cardiac action potential [23]. Blockage of the
hERG channel by several drugs is known to cause a loss of function leading to serious life-threatening disorders such as long QT prolongation (LQTS). Several drugs including antihistamines, antibacterials, antipschotics, and antidepressants have been identified as LQTS inducer and some drugs have been either withdrawn from the market (terfenadine, cisapride, astemizol, etc.) or restricted in use (thioridazine, haloperidol, sertindole, and pimozide, etc.) due to life-threatening cardiac arrhythmias [24,25]. Although experimental information (IC50, Cmax, etc.) of
these drugs against hERG channel [21,22] are available, the absence of the crystal structure of hERG ion channel in apo form as well as X-ray structures of these compounds with the channel limits the perception of molecular determinants on the interactions of these PDE5 inhibitors with the hERG channel. The mandatory incorporation of hERG assays on the developing stage of new drugs prior to their
marketing leads researchers to test their molecules in terms of hERG blockage. For instance, the newly developed two PDE5 molecules were tested for their ability to cross-react with PDE1-PDE4, PDE6, and PDE11; and also their inhibition potency of the hERG current was checked by in vivo analysis [26]. The investigation of the electrophysiological features into an another PDE5 inhibitor (ER-118,585) has underlined the risk of cardiac action potential duration according to the applied in
vivo experiments [27].
The aim of the present study is to identify the interaction patterns and binding affinity predictions of the selected PDE5 inhibitors (sildenafil, vardenafil and tadalafil) against the hERG1 channel in its open and open-inactivated conformational states as well as attempting to identify PDE5 inhibitor analogs with lower binding affinity to hERG1 ion channel while keeping the pharmacological activity against its principal target PDE5 using in silico methods such as molecular docking, virtual screening, MD simulations, and Alanine mutagenesis studies. The IC50 values of
these selected PDE5 inhibitors at the hERG1 channel as well as their 2D chemical structures are shown in Figure 2.1. In this sense, current study will particularly contribute to the understanding of the underlying molecular mechanism and molecular determinants (key residue-ligand interactions) that are responsible for hERG1 binding affinity in its two distinct conformational states (i.e. open and open-inactivated states).
2.2 Methods
2D structures of the ligands were drawn with the Builder tool in MOE molecular modeling package (MOE, 2015) [28]. Energy minimization and conformational search of the ligands were conducted with MMFF94x force field. Ligand and protein preparation steps were performed in the wash and protein preparation modules of the MOE program, respectively. The unmutated and full-length crystal structure of PDE5 (PDB ID: 2H42) [29] was used as target structure in this study. While co-crystallized ligand was removed from the active site of the enzyme, crystal water molecules within the binding pocket were retained for the docking simulations. Three different docking algorithms, namely GOLD [30], AUTODOCK [31], and MOE docking were used in the current study.
Figure 2.1 : 2D chemical structures of sildenafil, vardenafil and tadalafil and their corresponding IC50 values for the hERG K+ ion channel.
2.2.1 Molecular docking simulations
GOLD (v.5.3.0): Consensus docking protocol was used to generate protein-ligand complexes with GOLD 5.3.0 software. In this respect, two docking scoring functions were combined: GoldScore (for docking) and ChemScore (for re-scoring). Ten amino acids within the binding cavity of PDE5 (Asn661, Asn662, Ser663, Tyr664, Ile665, Ile778, Phe786, Met816, Gln817 and Phe820) were handled with full flexibility during the docking simulations. The first five aminoacids belong to the H-loop (660–683) whereas the others are part of the Glutamine switch pocket (Q pocket) and Hydrophobic clamp pocket (P pocket) where the critical drug-protein interactions take place. The flexible dihedral angles of Phe656 and Tyr652 of hERG ion channel were allowed full rotation with 10o increment. The number of genetic algorithm (GA) run was set to 100. Search efficiency was set to its maximum value (200%)-exploring the search space as wide as possible-in order to increase the reliability of the docking results.
MOE (v.2014.09): Therapeutic receptor (PDE5) is treated rigid during the docking simulations. However, side chains of the central inner cavity of the off-target receptor (hERG) were tethered with a weight factor of 10 via induced fit docking protocol. Two different docking scores were utilized in order to evaluate the predicted binding free energy of the ligands to these receptors: 1-London dG and 2-GBVI/WSA dG (Generalized-Born Volume Integral/Weighted Surface Area dG). Triangle Matcher was chosen as the ligand placement methodology. MMFF94x force field is used to refine the free energy of binding in the second refinement step. 100 poses were generated in each re-scoring steps.
AUTODOCK (v.4.2): The active site of PDE5 enzyme was defined based on the sildenafil coordinates in the crystal structure allowing enough space for ligand movement during the docking pose prediction. The grid spacing of 0.3 Å was used in grid generation calculations for hERG K+ ion channel and PDE5. The Gasteiger charges were added to the receptors and ligands. The docking search parameters were used with their default values during the docking simulations except that 100 docking poses were generated with GA for each ligand and Lamarckian Genetic Algorithm (LGA) was used for the conformational search step. The aromatic Tyr652 and Phe656 aminoacid residues on S6 helix regions of the hERG potassium ion channel were set as flexible during docking simulations. However, the coordinates of the crystal structure of PDE5 enzyme were treated as rigid.
2.2.2 Fragment-based de novo drug design & virtual screening
The sulfone and methylpiperazine groups are truncated from the molecule and replaced with the fragment library of MOE. Scaffold replacement in MOE resulted more than 100.000 sildenafil analogs. Some physiochemical restrictions were applied as filter during the virtual library generation in order to give sufficient drug-like properties to the derived molecules. After the attachment of fragments to the remaining core of the sildenafil, the generated molecules were filtered according to their lipophilicity (SlogP), molecular weight (MW), and topological surface area (TPSA) values. These descriptors are calculated from the atomic connection of the molecules and can be rapidly used for the indication of the absorption and bioavailability features of the drug-like candidate molecules. The following criteria were used in designing the virtual library: MW <500; 40 <TPSA <140; −4 < SlogP <8
where SlogP measures logarithm of the octanol/water partition coefficient. Molecules that fail to pass these physiochemical requirements were rejected and eliminated from the library. Finally, compounds that are passed from the Lipinski-drug-like test (around 100.000 molecules) were docked into the inner cavity of the hERG1 channel in its open-conformational state using virtual screening (VS) mode of the docking protocol implemented in MOE. Additionally, molecules having 70% similarity with sildenafil (i.e.,118 molecules) were downloaded from ZINC database (http://zinc.docking.org) in order to test their activity to PDE5 and hERG1 targets. Energy minimizations and partial charge calculations of the molecules were
performed by MMF94X force field with 0.001 kcal.mol−1.Å−2 RMSD gradient sensitivity. 10 conformers for each molecule were generated with 0.005 RMSD gradient by the Low-Mod MD method. The lowest energy conformer of each molecule was chosen to be docked into central cavity of the hERG1 ion channel in its open state with GOLD docking program. The flow chart for the whole VS protocol is summarized in Figure 2.2. For structure-based virtual screening, triangle matcher placement was used to generate 30 poses for each ligand combined with London dG rescoring function using MOE software. hERG open state conformation for the virtual screening protocol was used. Around 900 molecules were then selected (according to their relative low London dG energy of binding) for further docking simulations in the GOLD program at the end of the virtual screening step. The statistical analysis of the VS results in terms of binding energy prediction is shown in Figure 2.3.
2.2.3 MD simulations and post-processing MD analyses
MD simulations were performed for selected compounds by Gromacs 4.6.5 package [32].
2.2.4 MD simulations of the target receptor: PDE5 in its apo state and bound with its inhibitors
Energy minimizations of the systems were conducted by steepest descent (SD) integrator using an initial step size of 0.01 nm with 10.000 iterations and 1000 kJ/mol.nm minimum force tolerance. 5 ns position restrained dynamics in NPT ensemble with constraints on all bonds were applied prior to production MD stage. 50 ns production MD simulations were realized. 2 fs time-step was used for the simulations. Electrostatic interactions were calculated with Particle Mesh Ewald technique with a cut-off distance of 9.0 Å and van der Waals interactions were considered with a cut-off distance of 14.0 Å. Berendsen thermostat and Berendsen barostat were used as temperature and pressure coupling algorithms. The Linear Constraint method (LINCS) was used in order to fix all the bond lengths. Gromos43A1 force field was used for the simulations with leap-frog algorithm. Ligand topologies were prepared with the PRODRG Server [33]. Simulations were
conducted in cubic simulation boxes solvated with SPC waters and neutralized with Cl− ions.
Figure 2.2 : Flowchart of the current study in designing the novel sildenafil-like molecules based on fragment replacement strategy in order to obtain promising
molecules whose binding affinity to hERG K+ ion channel is decreased as the pharmacological activity against their therapeutic receptors, PDE5, is retained. 2.2.5 MD simulations of hERG K+ ion channel: Apo state and bound with PDE5
inhibitors
S5–S6 domains of the hERG K+ ion channel were considered in MD simulations. CHARMM-GUI Service [34] was used for the preparation of the protein-membrane systems. hERG1 channels in its open state and open-inactivated state were inserted into lipid bilayer membrane that contains DPPC-type lipid molecules. After solvation of the systems with TIP3P waters and adding neutralizing ions, systems energy were minimized using steepest descent algorithm with Verlet cut-off scheme. 5000 iteration steps and 1000 kJ/mol.nm minimum force tolerance were used for this step. LINCS algorithm was used in order to put constraints on hydrogen bonds. The Particle Mesh Ewald method was used in order to treat the long-range electrostatics
with a cutoff distance of 12 Å. van der Waals interactions were also considered within 12 Å distance. Systems were equilibrated with six steps of default equilibration scheme, provided by CHARMM-GUI, using leap-frog integrator with Verlet cut-off scheme. Berendsen thermostat and Berendsen barostat were used during the equilibration period for the temperature and pressure control. Nose-Hoover and Parrinello-Rahman coupling algorithms were utilized in order to maintain 310 K and 1 atm simulation conditions, respectively, during the 50 ns production stages. CHARMM36 force field was used for the simulations.
Representative frames of the simulations were generated from the last 10 ns trajectories for the hERG/sildenafil, PDE5/sildenafil, and hERG1/sildenafil_frag21, PDE5/sildenafil_frag21 bounded systems and further subjected to MM-GBSA free energy calculations with Prime program [35,36].
Figure 2.3 : Virtual screening analysis. Binding energy values are plotted against the frequency (number of molecules in the corresponding energy range). Values are obtained at the end of the VS protocol via docking the sildenafil analogs against
