Burghardt, B.

Merlitz, H.

Wenzel, W.

English

KITopen

the active site of dihydrofolate reductase (4dfr).For this target, a highly optimized inhibitor isknown (methotrexate, MTX), including the crys-tallographic X-ray structure of its docked confor-mation. MTX emerges by far as the top scoringligand in the screened database, validating bothdocking approach and the choice of the scoringfunction. It is shown that MTX and other topscoring database leads coordinate with preferredamino acids for the creation of hydrogen bondsinside the cavity. We discuss how this observationcould be employed to identify hot-spots in theactive site of the protein and thus support thecreation of pharmacophore models.2. Computational methods2.1. Preparation of the three-dimensional databaseThe ligands were taken from the open part ofthe NCI database (nciopen.mol), which had beenprocessed to generate the three-dimensionalstructures with the help of Corina [19,20] (nci-open3d.mol). Hydrogen atoms and partial chargeswere added using the InsightII software packageand its esff forcefield. Rotatable bonds were iden-tified with a simple algorithm which searched forsingle bonds, excluding ring-structures, trivialsingle-atom end-groups and atoms with sp2 hybridorbitals. Among 125 000 compounds available, thefirst 10 000 were selected which satisfied the con-ditions of having not more than 100 atoms and notmore than 10 rotatable bonds.2.2. Preparation of the protein coordinates anddefinition of the active sitesThe protein coordinates were taken from theX-ray structure of Escherichia coli dihydrofolatereductase with MTX (pdb entry: 4dfr, monomerB) [21]. Hydrogen atoms and partial charges wereattached using InsightII with esff force field. Thebinding site was defined with the docked ligand inthe X-ray structure: A volume of 5 A radius wasdefined around the center of mass position of thedocked MTX. If the global minimum on the PESof any ligand was located such that either thecenter of mass of this ligand or more than 20 li-gand atoms were inside this volume, the ligand wasregarded as docked.2.3. Scoring functionFor the simulations discussed below we used thefollowing scoring functionS ¼XProteinXLigandRijr12ij Aijr6ijþ qiqjrij!þXH-bondscosHij~Rijr12ij~Aijr10ij!; ð1Þwhich contains the empirical Pauli repulsion (firstterm), the Van de Waals attraction (second term),the electrostatic Coulomb potential (third term)and the angular dependent hydrogen bond po-tential (term four and five). The Lennard Jonesparameters Rij and Aij were taken from OPLSAA[22], the partial charges qi were computed withInsightII and esff force field, and the hydrogenbond parameters ~Rij, ~Aij were taken from Auto-Dock [10]. The electrostatic potential was pre-calculated on a grid. The other contributions,however, turned out to be too rapidly varying tobe reliably stored on a grid and interpolated. In-stead we exploited the short-range character ofthese interactions. We stored a list of atoms next toeach grid point and computed the short-rangepotentials on the fly. This procedure obviates theneed for interpolation of rapidly changing termsand permits a straightforward implementation ofthe angle dependence of the hydrogen bond po-tentials. With these changes the accuracy of theevaluation of the scoring function could be sig-nificantly increased over an earlier entirely grid-based implementation [18].2.4. Optimization methodMany stochastic optimization methods getroutinely trapped in metastable minima of thepotential energy surface. This problem is particu-larly acute in receptor ligand docking, where thereis little specificity between ligand and receptoroutside the receptor pocket and little room forreorientation inside. The degree of affinity betweenreceptor and ligand is discernable only where closeto the ideal docking position, in particular becauseof the short-range nature of stabilizing hydrogenbonds. In the proximity of typical docked posi-tions, however, the ligand cannot be significantlyrotated or distorted without clashing with receptoratoms.The STUN method [16] avoids this problem byletting the particle in the minimization processtunnel forbidden regions of the PES. As in sim-ulated annealing [23] we retain the idea of a biasedrandom walk, but apply a non-linear transforma-tion to the potential energy surfaceESTUNðxÞ ¼ 1 exp½cðEðxÞ  E0Þ; ð2Þwhere E0 is the lowest minimum encountered bythe dynamical process so far. Alternately a suitableupper bound for the global minimum can be usedfor E0. This effective potential preserves the loca-tions of all minima, but maps the entire energyspace from E0 to the maximum of the potentialonto the interval [0,1]. At a given finite tempera-ture of O(1), the dynamical process can thereforepass through energy barriers of arbitrary height,while the low energy-region is resolved even betterthan in the original potential. The degree ofsteepness of the cutoff is controlled by the tun-neling parameter c.2.5. ComputationThe electrostatic potentials were stored on a30 30 30 A3grid with grid spacing of 0.5 A,which was positioned around the center of mass ofthe X-ray structure of MTX. Each of the 10 000ligands was shifted into the cavity and, if required,moved so that heavy clashes of the initial config-uration were avoided. We then performed 10STUN simulations, comprising 3 105 energyevaluations each. Only the repetition of the dock-ing procedure ensured that weakly bound ligandsare reliably docked. The (non-physical) simulationtemperature was fixed to 5 K, the STUN param-eter c was chosen as 0.01 for all ligands. A driftterm was applied in order to increase the samplingrate inside the cavity with respect to peripheralregions [18]. The size of the displacements in theconformational updates was automatically ad-justed during the simulation to keep the accep-tance rate around 0.5. In average, each ligandconsumed about 20 seconds of CPU time for 105updates on a PC (AMD Athlon 1.2 GHz).3. Results6100 ligands, out of 10 000 compounds, reachedconfigurations where the ligand was embedded inthe cavity and the external binding energy, scoredwith Eq. (1), was below )50 kJ/mol. It is wellknown that scoring functions as the one used hereare too inaccurate to yield quantitatively approx-imate the natural affinity of the ligand. Neverthe-less they often provide a useful relative ranking ofleads when docked under identical conditions.Fig. 1 shows the external binding energies of thedocked compounds. Among them, the natural li-gand MTX clearly scored best, with a minimumenergy conformation which differed only 1.4 Afrom the X-ray crystal structure (Fig. 2). Thisfigure shows that both conformations essentiallydiffer in the positions of only one of their car-boxylate groups. In the minimal conformation ofthe scoring function this group is turned around toFig. 1. The affinities of 6100 ligands docked to 4dfr. The natural ligand (MTX) was scoring best. The five top rankingdatabase leads are denoted as a; . . . ; e.create an extra hydrogen bond to Lys-32 which isabsent in the crystallographic conformation. In theexperimental configuration a conserved watermolecule (Wat-672) mediates competing hydrogenbonds in the natural environment, an effect whichcannot be accounted for with the present scoringfunction.The results of the screening process were en-couraging not only because the natural ligand wasselected with high specificity, but also becauseother high ranking ligands showed significantsimilarities in the binding mechanism. This ob-servation, more than the exact ranking in thebinding affinity, suggests that virtual screeningmethods can aid rational drug design. In Table 1we compare some of the relevant features of MTXwith the five ligands which scored next. It is ob-vious that no other ligand could create as manyhydrogen bonds with the protein as MTX. On theother hand, it appears that there exist certain res-idues inside the cavity which serve as preferredbinding sites, in particular Arg-52 and Asp-27.One may hypothesize that these residues, alongwith the general shape of the receptor pocket, areresponsible for the specificity of the receptor li-gand interaction. We therefore replaced bothamino acids with valine, which cannot act as ahydrogen bond donor/acceptor. The resulting af-finities of a repeated screening process are shownin Fig. 3. Not surprisingly, the binding energies arenow generally reduced and those ligands whichhad originally scored well have lost their topranking positions. This exercise demonstrates howvirtual screening can be employed to identify hot-spots inside the active site of the protein to supportthe creation of pharmacophore models.Finally, the reliability of docking MTX and thetop ranking leads was investigated. The dockingcycle was defined as successful if the lowest mini-mum led to a ligand conformation which wascloser than 2 A to the best known minimum whichFig. 3. The affinities of the ligands after docking to a modified4dfr, in which Arg 52 and Asp 27 are replaced by Valine.Ligands which originally scored best have now lost their topranking positions, including MTX.Fig. 2. MTX in its crystallographic X ray conformation (black)and the global minimum found with STUN (color). Bothconformations differ in the positions of their carboxylategroups.Table 1Properties of ligand MTX and the 5 top ranking leads of thedatabase, denoted as a eLead #Atoms #Rot. Affinity Hydrogen bonds to(kJ/mol)MTX 54 9 233 Arg 52, Asp 27, Arg 57,Lys 32, Ile 5a 46 10 173 Arg 52, Ile 14, Ala 7b 28 8 169 Asp 27, Trp 22, Leu 24c 38 7 164 Arg 52, Ile 50, Ala 7,Ile 94d 29 4 163 Arg 52, Asp 27, Trp 22e 37 9 156 Arg 52, Asp 27, Trp 22#Rot. means the number of rotatable bonds, i.e. the numberof internal degrees of freedom in the simulation.satisfied the docking criteria of Section 2.2. Fig. 4plots the average docking rates as a function ofenergy evaluations. We observe that after 106scorings none of the ligands were able to reach itsminimum conformation with more than 50%probability, but that significant docking rates werereached already after 105 energy evaluations. Gi-ven a finite amount of computational resources itis therefore advantageous to perform severalshorter docking runs instead of few long cycles tolocalize the global minimum. One exception is leada: Here, a docking position outside the cavity ex-ists which is 20 kJ/mol lower than the best con-formation inside the cavity. As a result the rate ofdocking inside the cavity remained very low, in-dependent on the length of the cycle. Interesting isthe fact that lead e was docking even faster andmore reliable than MTX. After only 104 scorings adocking rate of 25% was reached. This may indi-cate that the configuration space of this particularreceptor ligand system is especially favorable forreaching the global minimum.4. SummaryIn this investigation we have performed an au-tomated virtual screening of 10 000 compounds,which were taken from the openly available data-base nciopen3D, to the active site of E. coli dihy-drofolate reductase (4dfr). The global optimizationtechnique STUN, which was already successful inthe rigid ligand approximation, proved reliable andefficient also for the more demanding applicationto flexible ligands. It was able to match the naturalligand MTX (which contains 9 rotatable bonds) toan accuracy of 1.4 A with the X-ray structure. Thescoring function Eq. (1) was ranking this ligand asthe best among the 10 000 screened molecules.We have further demonstrated how a virtualdatabase screening may serve as a tool for phar-macophore modeling, because the top rankingleads are useful probes to identify hot spots insidethe cavity. Despite the fact that there appear toexist leads that docked faster and more reliablethan MTX, none of them reached a similarly highaffinity. We hope this study contributes to thesearch of new powerful techniques for virtualscreening especially under the perspective of futureprojects which ask for additional degrees of free-dom by means of rotatable groups in the receptorand the inclusion of solvent molecules. Our resultsmotivate further investigations which employSTUN as a global optimization technique in suchapplications.AcknowledgementsWe thank the Fond der Chemischen Industrie,the BMBF, the Deutsche Forschungsgemeinschaft(Grant WE 1863/11-1) and the Kurt EberhardBode Stiftung for financial support.References[1] W. Walters, M. Stahl, M. Murcko, Drug Discov. Today 3(1998) 160.[2] C. Bissantz, G. Folkerts, D. Rognan, J. Med. Chem. 43(2000) 4759.[3] I. Kuntz, J. Blaney, S. Oatley, R. Langridge, T. Ferrin,J. Mol. Biol. 161 (1982) 269.[4] E. Perola, K. Xu, T. Kollmeyer, S. Kaufmann, F.Prendergast, Y. Pang, J. Med. Chem. 43 (2000) 401.[5] M. Rarey, B. Kramer, T. Lengauer, G. Klebe, J. Mol. Biol.261 (1996) 470.[6] W. Welch, J. Ruppert, A. Jain, Chem. Biol. 3 (1996)449.Fig. 4. Docking rates of MTX and the five top ranking leads asa function of energy evaluations.[7] C. Baxter, C. Murray, D. Clark, D. Westhead, M.Eldridge, Proteins 33 (1998) 367.[8] T. Hou, L. Chen, X. Xu, Protein Eng. 12 (1999) 639.[9] G. Jones, P. Wilett, R. Glen, A. Leach, R. Taylor, J. Mol.Biol. 267 (1997) 721.[10] G. Morris, D. Goodsell, R. Halliday, R. Huey, W. Hart, R.Belew, A. Olson, J. Comput. Chem. 19 (1998) 1639.[11] P. Charifson, J. Corkery, M. Murcko, W. Walters, J. Med.Chem. 42 (1999) 5100.[12] D. Gehlhaar, G. Verkhivker, P. Rejto, C. Sherman, D.Fogel, L. Fogel, S. Freer, Chem. Biol. 2 (1995) 317.[13] D. Goodsell, A. Olson, Proteins 8 (1990) 195.[14] M. Liu, S. Wang, J. Comput. Aided Mol. Des. 13 (1999)435.[15] C. McMartin, R. Bohacek, J. Comput. Aided Mol. Des. 11(1997) 333.[16] W. Wenzel, K. Hamacher, Phys. Rev. Lett. 82 (1999) 3003.[17] U.H.E. Hansmann, Eur. Phys. J. B 12 (1999) 607.[18] H. Merlitz, W. Wenzel, Chem. Phys. Lett. 362 (2002) 271.[19] J. Sadowski, J. Gasteiger, Chem. Rev. 93 (1993) 2567.[20] J. Sadowski, J. Gasteiger, G. Klebe, J. Chem. Inf. Comput.Sci. 34 (1994) 1000.[21] J. Bolin, D. Filman, A. Matthews, R. Hamlin, J. Kraut,J. Biol. Chem. 257 (1982) 13650.[22] W. Jorgensen, N. McDonald, J. Mol. Struct. 424 (1997)145.[23] S. Kirkpatrick, C. Gelatt, M. Vecchi, Science 220 (1983)671.

Application of the stochastic tunneling method to high throughput database screening

https://publikationen.bibliothek.kit.edu/110054330/149887249

Application of the stochastic tunneling method to high throughput database screening

Abstract

Similar works

Full text

Available Versions

KITopen