Rational Drug Design - ACS Publications - American Chemical Society

Chapter 18

Further Development of a Genetic Algorithm for Ligand Docking and Its Application to Screening Combinatorial Libraries 1,5

1

2

3

Gareth Jones , Peter Willett Robert C . Glen , Andrew R. Leach , and Robin Taylor 4

Downloaded by EAST CAROLINA UNIV on December 25, 2017 | http://pubs.acs.org Publication Date: July 7, 1999 | doi: 10.1021/bk-1999-0719.ch018

1

Krebs Institute for Biomolecular Research and Department of Information Studies, University of Sheffield, Western Bank, Sheffield S10 2 T N , United Kingdom Tripos Inc., 1699 South Hanley Road, St. Louis, MO 63144 Glaxo Wellcome Medicines Research Centre, Gunnells Wood Road, Stevenage SG1 2NY, United Kingdom Cambridge Crystallographic Data Centre, 12 Union Road, Cambridge CB2 1EZ, United Kingdom 2

3

4

We have previously reported the development and validation of the program GOLD (Genetic Optimization for Ligand Docking), an automated ligand docking program that uses a genetic algorithm to explore the full range of ligand conformational flexibility with partial flexibility of the protein. The validation of the algorithm exposed certain weaknesses, particularly in its handling of hydrophobic ligands. We describe here a number of modifications designed to overcome these shortcomings. Using the improved algorithm we have investigated the relationship between the program's scoring function and ligand activity and applied the algorithm to screening combinatorial libraries. Using lipase as a target, GOLD docked a library of carboxylic acids and amines extracted from the Fine Chemicals Database. By predicting the binding modes of product from the predicted dockings of reactants the problem is considerably simplified. The success of this approach was verified by comparing the prediction obtained by docking the component reactants with the result obtained by docking the product. The procedure was repeated using a second library of sulphonyl chlorides and amines. Prediction o f small-molecule binding modes to macromolecules o f known threedimensional structure is a problem o f paramount importance i n rational drug design (the "docking" problem). There have been many different approaches to solving the docking problem (1,2), including: deterministic approaches (3,4); simulated annealing 5

Current address: Arena Pharmaceuticals, 6166 Nancy Ridge Drive, San Diego, CA 92121

© 1999 American Chemical Society

Parrill and Reddy; Rational Drug Design ACS Symposium Series; American Chemical Society: Washington, DC, 1999.

271

272


(5); genetic algorithms ( G A s , 6,7) and evolutionary prograrnming (8). G A s provide an evolutionary search paradigm that enables the rapid identification o f good, though not necessarily optimal, solutions to combinatorial optimization problems (9,10,11). This technique has proved remarkably successful i n almost all areas o f computational chemistry and molecular modeling, to which it has been applied (12). G O L D (6,13) is an automated ligand docking program that uses a G A to explore the full range o f ligand conformational flexibility. Conformational flexibility is encoded i n the G A i n binary bitstrings where the angle o f rotation about each rotatable bond is encoded i n a byte. In order that G O L D efficiently elucidates binding modes, hydrogen bonds have been directly encoded into the G A , using integer strings, such that each position on the string ecoded a particular ligand donor-hydrogen or acceptor while the value at that position encoded a complementary protein acceptor or donor-hydrogen. When decoding a G A chromosome least-squares fitting is employed to try and form some o f these encoded hydrogen bonds. A simple scoring function is used to rank generated binding modes. This comprised a term for hydrogen bonding (which took account o f the fundamental requirement that water must be displaced from both donor and acceptor before a bond is formed); a pairwise dispersion potential between the protein and the ligand and a molecular mechanics term for the internal energy for the ligand. The program has been extensively evaluated on a large number of protein-ligand complexes (13). The validation studies o f G O L D revealed weaknesses i n predicting the binding mode o f hydrophobic ligands (75). Here, we describe a number o f improvements, designed primarily to overcome this limitation. Using the improved algorithm a number o f activity studies have been carried out and results on sets o f influenza A neuraminidase and a-chymotrypsin ligands are presented. Finally, G O L D has been used to screen combinatorial libraries. Given the large size o f combinatorial libraries, the selection o f a small number o f reactants or products for further study against a known target is both a demanding and important problem i n rational drug discovery. Rather than enumerating a library and docking every library compound, we simplify the problem by docking only library monomers. The binding modes o f high affinity products can then be elucidated from the docked monomers. A knowledge o f how G O L D works would be helpful for the first section o f this chapter which deals with improvements to the program (see references 6 and 13). However, the sections on activity prediction and library screening can be read separately without any great understanding o f G O L D . E x t e n d i n g the Test Suite. In (73) we reported the validation o f G O L D on 100 test complexes extracted from the Protein Data Bank ( P D B ) (14). Shortly afterwards, one o f us ( R T ) expanded this test suite to 134 complexes. Further expansion without including complexes very similar to those selected, or biasing the dataset towards peptidic ligands did not seem possible. Using the same subjective classification o f Good (where all protein ligand interactions were correctly predicted), Close (the binding mode and all important interactions were elucidated), Errors (the prediction


273

was partially correct but had significant errors) and Wrong (where the predictioun was completely wrong) we obtained the results shown i n Table I. Overall we achieved prediction rate o f 72%, with 96 complexes being i n the Good and Close categories. A l l results can be viewed at http://www.ccdc.cam.ac.uk/prods/gold.htrnl.


Table I. Results o f docking predictions on a dataset o f 134 complexes. missing complex is 1 A C L , for which no binding mode could be elucidated. Subjective Result Good

Number

PDB identification

53

Close

43

Errors

15

Wrong

22

1ABE 1FKI 1P0C 2CTC 8GCH 1GLP 1BLH 1LDM 2PK4 1APT 1HFC 1BAF 2CMD 1AAQ IMUP 2MCP

1ACM 1HDY 1SRJ 2PHH 1AEC 1LAH 1DIE 1PHA 2YHX 1AZM 1IMB 1EAP 1CTR 1ACJ 2R07 ICDG

1AC0 1HEF 1STP 2SIM 1AHA 1LPM 1DR1 1PHG 3CPA 4EST 1LCP 1ETR 2LGS 1DID INIS ILMO

The

codes 1CBX 1HYT 1TPP 3AAH 1ASE 1MMQ 1DWD 1RNE 3GCH 1ATL 1NC0 1HDC 1LNA 1EED ITDB ITYL

1C0Y 1LST 1ULB 3PTB 1HSL 1MRG 1EPB 1SLT 3HVT 1BBP 1TNG 1LIC 1SNC 1ETA 2AK3

1CPS 1MDR 1XIE 3TPI 1BMA 1TRK 1GHB 1TKA 4CTS 1BYB 1TNI 1RDS 1UKZ 1HRI 2MTH

1DBB 1MRK 2 ADA 4DFR 1CIL 1TNL 1GLQ 1TMN 5P2P 1CBS 1TPH 1R0B

1DBJ 1PBD 2CGR 4PHV 1FRP 1WAP 1IDA 1XID 6ABP 1C0M

1FKG 1PHD 2CHT 7TIM 2GBP 1IVE 2DBL 6RNT 1FEN

6RSA 1ACK

1ICN 1IGJ 1MCR 2PLV 3CLA 4 FAB

In (75) several flaws were identified i n the algorithm. This extended verification confirmed these shortcomings, which were mainly related to the recognition o f hydrophobic interactions, and we now describe a number o f improvements that have been made to correct the algorithm. Improving the Algorithm Increased V D W Weighting. The fitness function used by G O L D comprised three main terms: a hydrogen bonding energy; a pairwise van der Waals steric energy between the protein and the ligand; and an internal energy term for the ligand. The total fitness score was the linear sum o f these terms, with each term being equally weighted. During testing it was observed that there was insufficient weight given to large areas o f hydrophobic contact and that there was insufficient probability o f charged ligand groups being solvent accessible. This was attributed to an imbalance i n the fitness function, with hydrophobic contact being undervalued relative to hydrogen bonding, since, while the van der Waals interaction between the protein and the ligand accounts for dispersive interactions, it fails to score for the favorable entropic effects o f water displacement. It was felt that this imbalace might be easily corrected by scaling the contribution the van der Waals energy between the protein and the ligand made to the fitness function by a weight. Using the results from the experiments i n (75) each test system was examined i n turn. The fitness scores o f the 20 solutions were recalculated using a new trail weighting for the the van der Waals energy


274


between the protein and the ligand. The binding mode with the highest fitness score was then selected as a prediction. This procedure was repeated using a number o f trial weights and a factor o f 1.375 was found to give the most correct predictions over the 100 complexes. This new weight was then applied within the algorithm and proved successful at correcting the imbalance between hydrophobic and polar interactions within the fitness function. Several complexes such as 4 F A B and 1 A C J were now correctly docked by G O L D . The ability o f G O L D to dock polar ligands was unaffected. A s an alternative to increasing the van der Waals weighting an approximate surface area calculation was used to determine desolvation energies (75). However, this proved no more discriminatory than increasing the van der Waals energy weighting and was very time consuming. Shape Fitting. Another problem observed during testing was a systematic problem i n placing hydrophobic groups i n deep cavities. The chromosome encoding used i n G O L D meant that the G A searched patterns o f hydrogen bonding motifs: the algorithm was thus directed to discover hydrogen bonding interactions, but found hydrophobic interactions by "accident". Moreover, given the tendency for bad bumps to occur as the ligand approaches deep cavities, cavity binding may actually have been discouraged. The chromosome encoding has now been extended to encompass hydrophobic interactions. A grid with inter-point spacing o f 0.25A was placed across the active site. A n augmented carbon probe atom was placed at each grid point and the van der Waals interaction between the probe and the active site measured. This probe atom had a van der Waals radius o f 2.5 A and the physical properties o f an sp carbon. If the interaction energy was less than -2.5kcal/mol then the grid point was considered to be a favorable position for ligand hydrophobic atoms and the point was labeled. Ligand carbons bonded to at least one hydrogen were considered to be hydrophobic and were also labeled. A new integer string was added to the chromosome which encoded possible hydrophobic interactions, such that the position on the string was a ligand hydrophobic atom label and the value on the string was a label o f a hydrophobic grid point. During least-squares fitting and the application o f genetic operators this string was handled i n an analogous fashion to the integer strings that encode for hydrogen bonding (as described i n detail by Jones et al. (6)). 3

This extension o f the chromosome encoding allowed G O L D to dock apolar ligands and ligands which did not form any hydrogen bonds to the protein, and removed the need for the small ligand model (75). The binding modes o f complexes such as 1HRI and 2R07, which were originally incorrectly predicted, were now successfully elucidated, without any adverse effects being observed when docking hydrophilic ligands. Solvent Exposed Charged Residues. Another problem observed i n the testing o f G O L D was the binding o f ligand groups to charged solvated residues. The exterior o f a protein generally contains many charged solvated residues, such as lysines, arginines



275

and glutamic and aspartic acids. These residues, being highly solvated and mobile, are not normally favorable for ligand binding. However, with the G O L D algorithm such groups, being charged, are very attractive. A n approximate surface area algorithm was used to identify such residues. Using this algorithm the ratio o f surface areas o f the residue i n the absence and presence o f the rest o f the protein was determined. If this ratio was greater than 0.5 then the residue was considered solvent exposed. Such residues had their hydrogen bonding energies recalculated using a water dielectric o f 78.5, which made such residues much less attractive to ligand binding. For example, the binding mode o f the ligand i n 6 R S A was originally incorrectly predicted, with the ligand binding to a solvent exposed residue, but is now predicted correctly. Torsional Distributions. A s currently documented G O L D searches all single acyclic torsions over 360°. However, searches o f small molecule crystallographic databases, such as the C S D (16), show that many common torsions are restricted. Figure 1 shows two common torsions and their distribution i n the C S D . G O L D can now read such torsional distributions and search only those ligand torsional angles which are seen in small molecule crystals. B y default a small torsional library o f 24 common torsions is available, but the user can also use the library o f over 200 torsions from the M I M U M B A program (17). Re-mapping the Chromosome. In the docking procedure described i n (75) the leastsquares fitting process is applied i n two passes. The chromosome is decoded to give a mapping i n 3 D space between points i n the ligand and points i n the protein active site. The first pass o f least-square fitting is applied to all points, while the second pass is applied to close contacts from the first pass. Thus a significant portion o f the chromosome w i l l not be used i n performing the final least-squares fitting. In the improved docking algorithm the unused chromosome mappings are reset to the dummy value o f ' - 1 ' (that is, they are un-mapped). The crossover operator has been altered so that if, at a particular position, one parent has a dummy value and the other parent does not, the child w i l l always inherit the non-dummy value. These changes have the effect that a child is much more likely to inherit meaningful chromosome mappings than previously and this was reflected i n improved fitness scores. Current Results on 100 Test Systems. The continued improvement o f the algorithm has resulted i n a superior rate o f prediction over the 100 test systems previously used to test our algorithm (75). Table II compares the root mean squared deviations ( R M S D ) o f G O L D predictions to the crystallographically observed binding mode obtained using the current algorithm with those previously reported (75). A s before, G O L D was run a number o f times and the solution with the highest score was retained as a prediction. U p to 20 dockings were performed per test system. Due to the increased reliability and reproducibility o f the algorithm it was felt that it was not necessary to generate 20 G A solutions, as, i n many cases all 20 predictions would be identical. So if, at any time, the best three solutions were all within 1.5A R M S D o f each other docking was terminated. However, i n at least two cases, i f the algorithm had run to 20 dockings a much better R M S D to the observed binding mode would have been obtained (previously the G A was always run 20 times). Examination o f the



276

-^c'

c

Acyclic

zooiTidWiaaa~raa~Ao

II

Acyclic X * H

olb 40.0 ado 120.018002005

Figure 1. Distributions for ester (on left) and alkyl torsions.


277

table shows that the new algorithm shows a clear improvement with 72% o f test systems having an R M S D o f 2.0A or less and 8 1 % having an R M S D o f 3.0A or less.


Table II. Current and previously reported results for G O L D on 100 P D B complexes listed i n (75). RMSD (A)

Total

0.5,1.0, 1.5, 2.0, 2.5, 3.0

15 12 7 2 19

Sum Total 11 45 60 72 79 81 100

Previous Sum Total 8 35 55 66 68 71 99

Comparison of Fitness Score With Binding Data A s we have shown i n validation experiments, G O L D is clearly capable o f reproducing the binding mode o f many ligands with considerable accuracy. However, it would also be desirable i f the program were able to predict the activity o f ligands. For this reason the relationship between G O L D fitness scores and the activity o f known ligands was explored for two test systems. Influenza A Neuraminidase. Activity and, where available, structural data for influenza A neuraminidase ( N A ) ligands was obtained (18-24). For many o f these compounds crystal structures were available (either publicly from the P D B , or as proprietary Glaxo Wellcome structures). The 34 compounds included 10 inactive molecules, which exhibited high structural similarity to some o f the active compounds (Neil Taylor, Personal Communication). Using G O L D , each ligand was docked 10 times (each docking comprised 100000 G A operations) into the active site o f N A and the solution with the highest G A fitness score retained. The results o f these experiments are shown i n Table III. From the table it can be seen that, with the exception o f the inactive B A N A ligands (18), G O L D always reproduced the observed binding mode. Figure 2 plots the relationship between observed IC50 values (logarithmic scale) and G O L D fitness scores. While the relationship is clearly not ideal there does appear to be a marked correlation and this is borne out by non-parametric statistical tests which show strong evidence o f rank correlation (Spearman test r =-0.65, / X 0 . 0 0 1 ; Kendall test, z=-0.48, p O . O O l ) . If we consider l O u M to be a cutoff for activity then there are 15 actives and 19 inactives i n the dataset. The question arises as to whether the G O L D score can be used to predict activity. Using a score o f 74 or greater to indicate activity we obtain Table I V . From the table it seems that the G O L D 5



278

Figure 2. Plot o f G O L D fitness score against

IC50

values for N A ligands.


279


Table III. Compound activities and Fitness scores for N A inhibitors. Where the crystal structure of the complex is available the PDB code (or Glaxo Wellcome structure id) and the RMSD between the GOLD prediction and observed binding mode are shown. Where no reference is given for a compound, the structure is proprietary. The compound id is that used in the reference. Compound Id

IC (pm)

3a(27) 41 (20) 2a(27) 2(27) 4j (20) 3b (21) 5e 2b (21) 4(27) 4g(20) la (27) DANA (19) DANA (18) EPANA (23) EPANA (23) 4d (20) la(27) \(24) 2(24) inactive4 inactive8 inactive 10 inactive3 inactive5 inactive7 inactive2 inactive 1 BANA105 (18) inactive6 inactive9 A P A N A (23) N A N A (22) BANA106 (18) BANA108 (18)

0.002 0.002 0.004 0.005 0.005 0.012 0.014 0.18 0.32 0.32 0.5 8.6 8.6 -10 -10 12 19 20 >100 >130 >130 >180 >210 >270 >390 >500 >640 750 >880 >900 -1000 -1000 10000 >20000

50

Fitness Score 82 .02 87 .39 76 .99 92 .96 81 .96 75 .52 85 .34 68 .93 91 .59 77 .78 74 .40 74,.73 75 .20 78 .98 77 .78 76,.24 67 .97 64 .14 83 .10 58 .86 71 .94 59 .31 66 .34 61 .80 72 .91 75 .18 64 .80 44 .85 62 .74 63 .13 78 .43 80 .66 46 .39 41 .03

RMSD Notes Crystal Structure Id 4 4 1.13 1 gs0023 1 gs0015 0.38 4 4 4 4 1 gs00154 0.32 4 4 lnnb 1 1.13 livf 1.18 2 liny 1 1.25 1.54 2 linx 4 4 4 5 4,6 4,6 4,6 4,6 4,6 4,6 4,6 4,6 2.55 3 livd 4,6 4,6 linw 1.28 2 2 2bat 1.21 2.30 2,3 live 2.84 2,3 live

1. Subtype N9 crystal structure. 2. Subtype N2 crystal structure. 3. GOLD failed to predict B A N A geometries correctly as these compounds contain non-planar amide bonds. However, the position of the benzene ring and acid group were correctly predicted. 4. Docked into the protein crystal structure of gs00023. 5. Docked into the protein crystal structure of gsOO 155. 6. Glaxo Wellcome inactives from assay data. These inactives have a high structural similarity to the actives.


280

score is a good indicator o f activity and it is most unlikely that this distribution could 2=

have arisen through chance (x 15.27,/?

Rational Drug Design - ACS Publications - American Chemical Society

Recommend Documents