Pharmaceutical case study Page 1 of 11 Smart Virtual Screening Protocol, Part III: Using Structure Based Drug Design to Enhance an In Silico Workflow. Industry Sector Pharmaceutical Luke Fisher, Shikha Varma, Teresa Lyons, Deqi Chen Organization Introduction McMaster University’s High Throughput Screening Lab organized a competitive research project for all who wished to participate . The idea was to design 1 Accelrys Key Products 2 Cerius ®·LigandFit 2 Cerius ·LibDock 2 Cerius ·LigandScore Discovery Studio® (DS) Accord™ for Excel TSAR® (CORINA) Catalyst® in silico experiments that rationalized the results from a replicate primary screen for E. coli dihydrofolate reductase (DHFR) inhibitors. In all, there were 49,995 molecules in the primary screen that yielded a set of 32 actives for DHFR . 2 In a recent case study using recursive partitioning (RP) , the full 3 dataset was used to build a predictive, decision tree model that yielded a 10fold enrichment above the primary screen results (32 actives out of 49, 995 = 0.06% yield). In the model, 30 of the 32 actives (94%) were correctly predicted as actives out of a set of only 3829 molecules (7% of full database = 0.7% yield). While the result of this RP case study is very useful as a standalone highthroughput vHTS model, we were also interested in seeing if an additional step using another in silico technique would further refine the results. The RP model in this study was built using a binary criterionInactive (HIGH) or Active (LOW)so all molecules in the study were predicted as one class or the other. As such, our goal in this case study is to take the results of the RP study and perform a structure-based drug design (SBDD) approach using the DHFR protein structure as a target. All molecules were used as input for docking into the protein active site and prioritized (ranked) using a series of scoring functions for protein-ligand interactions. Methods In a recent case study, an exhaustive protocol was performed to determine the best protein target for this SBDD work . 1RA2 was determined to be the preferred protein target for this project. Since we were interested in both competitive and non-competitive inhibitors for DHFR, we opted to remove the cofactor present in the complex, which led to a large binding pocket in the protein. We chose to use the Site Features Docking (SFD) tool in Cerius ·LigandFit v 4.9 (also known as LibDock) for our docking and scoring calculations. 6 2 5 4 Accelrys Corporate Headquarters 10188 Telesis Court, Suite 100 San Diego, CA 92121 United States Tel: +1 858 799 5000 Accelrys European Headquarters 334 Cambridge Science Park Cambridge, CB4 0WN, UK Tel: +44 1223 228500 Accelrys Asia Headquarters Nishi-shimbashi TS Bldg 11F Nishi-shimbashi 3-3-1, Minato-ku, Tokyo, 105-0003, Japan Tel: +81 3 3578 3860 Pharmaceuticals case study continued Preparation of the Input Molecules To ensure a clean dataset for our SBDD work, we used a combination of Accelrys' ‘Accord Enterprise Informatics (AEI) and the TSAR packages’ CORINA (v1.82) program to prepare the input file. AEI was used to check for duplicate structures in the original file given to us by McMaster University. CORINA was then used to convert the 2D format of each molecule in the SD file into a 3D format that was used as input for the next step of the project. This procedure ensured that all the structural redundancies were eliminated and also guaranteed that the structures were geometrically optimized for further study. Catalyst Database Creation To use SFD, the library of compounds to be docked in the active site must be in a Catalyst format either a Catalyst 3D, multi-conformer database or a series of CPD files (one for each molecule in the library). Here, we converted the 49,995 compound collection into a Catalyst database format. The database was built using the fast method, with 250 as the maximum number of conformers of each molecule. Statistics on this database reveal that the total number of conformers for each molecule was about 48. The average molecular weight for the dataset was 332 g/mol. Site Features Docking (SFD) The full database collection was docked into the 1RA2 active site. The first step in using SFD is to select a binding pocket and define a 3D box of protein atoms surrounding the site. With the Cerius GUI, this is straightforward using a number of techniques. In this work, we opted to generate a site model using the flood-filling algorithm in Cerius ·LigandFit to visualize the cavity. Next, we selected the center position of all the site points in order to define the center of the box for SFD. The box was then enlarged to ensure we captured all atoms and residues of interest for docking (see Figure 1). 2 2 9 7 8 Page 2 of 11 Figure 1: Site Features Docking Box Pharmaceuticals case study continued The next step for setting up SFD calculations is to prepare a series of ‘hotspots’ in the binding pocket. The hotspots are a series of interaction points, defined exclusively by the protein atoms. The hotspots are used to dock ligands that contain the equivalent binding features. For example, a ligand must have an apolar feature aligned with the apolar hotspot. Conversely, a polar feature in the ligand will need to overlay with the polar hotspots. In this work, we defined the hotspots using the default settings: 50 hotspots for both the apolar and polar sites. Figure 2 illustrates the hotspots, which have been overlaid with the x-ray bound folate molecule for reference. In this figure, the blue dots correspond to the polar hotspots and the yellow dots represent the apolar hotspots. Page 3 of 11 Figure 2: Hotspots created with Site Features Docking in the 1RA2 Active Site. Blue Dots: Polar Hotspots. Yellow Dots: Apolar Hotspots. Running SFD from the Command Line SFD can be run in either the GUI or by the command line. Due to the time required for this docking calculation, we decided to partition the job across multiple processors. All calculations were run on five SGI Octanes with dual R10000 processors. We used the following command line operations: Example: $C2DIR/exec/DockCatDB -p dockcatdb.par -o trs1.ldb -l trs1.log -b 1 -e 10001 -g 1 -f 1 -2 1 $C2DIR/exec/ldb2sd -i INPUT -o OUTPUT.SD -m 1 -s DJD_SCORE The first executable (DockCatDB) was used to partition the input using the −b and −e flags. In the example just cited, the job was set up to dock all molecules between 1 and 10,000. The −g flag indicates that we opted to calculate a grid for the docking calculations. The −f option is used to flag the desired level of quality and speed of the calculations. All of our calculations were then run using the medium speed setting (value=1). Pharmaceuticals case study continued The -2 option flag tells SFD how to handle dihedral angles between two sp2 atoms during the minimization step. The second executable (ldb2sd) was used to convert the binary .ldb file into a standard 3D SD file. Additional parameters were set via a hidden .LibDock file in the run directory. The only alteration to this set of parameters was changing the tolerance value to 0.7. This higher value for tolerance allowed for a more thorough search based on the SFD matching procedure. During the course of the project, we initially decided to dock the entire collection of compounds through SFD. We used the medium speed mode in SFD, which yielded 25 poses per molecule, or 1.25 million poses in the full collection. However, since we were interested in only the output from RP, we used a Perl script to extract all poses from the output SD files for the 3829 molecules that had been predicted as active for the RP model . Also important to note is that 167 molecules were not docked with SFD since the minimum overlap between hotspots and the ligand was not achieved. This lead to a final total of 3662 compounds that were included in the analysis. One caveat of the SFD algorithm is that all hydrogens are cleaved from all poses. Since the primary scoring function for SFD (DJD) is calculated without hydrogen dependency, the score is not affected. However, for this work, we were interested in utilizing all of the scoring functions in Cerius in order to prioritize and rank our results. As such, we needed to reattach the hydrogens to the docked poses. This was accomplished with the Open Babel v1.100.1 program because it is able to handle large SD files using a command line executable. After the hydrogens were reattached, there was a distinct possibility that the new hydrogen atoms would be too close to protein atoms (van der Waals clashes). Because of this fact, we minimized each pose rigidly in the active site to ensure these clashes would not occur. Working with Accelrys development scientists, we were provided with a patch release of Cerius 2 2 10 Page 4 of 11 LigandFit that allowed for Rigid Body Minimization (RBM) of all of the poses that contained the re-attached hydrogens. This procedure was extremely fast and produced more accurate docking poses that are needed for scoring. Following the RBM of all poses, we calculated the following protein-ligand scoring functions in Cerius LigandFit: LigScore1, LigScore2, PLP1, PMF , and Jain. 14 15 2 12 PLP2 , 13 Results and Discussion The set of molecules in the training set contained 32 active molecules defined by a 75% residual activity measurement. Figure 3 lists the 32 active compounds, which were generated with DS Accord for Excel (v 5.2). 16 Pharmaceuticals case study continued Page 5 of 11 Figure 3: Known 32 Actives in the Full Dataset To analyze the effectiveness of all of the Cerius scoring functions, we first extracted the best-scored pose for a molecule using the Cerius SD File Filter option. The scores were output to Microsoft Excel to generate hit rate plots for all scoring functions. The question we needed to address was as follows: Which scoring functions are performing well for this DHFR set of actives? 2 2 Pharmaceuticals case study continued Page 6 of 11 We approached this question by examining the rank order of all actives versus that of all the molecules in the full dataset. In other words, a scoring function is performing well if a large percentage of actives is ranked high (retrieved) in a small percentage of the full set. The enrichment plot for all scoring functions is shown in Figure 4. Figure 4: Enrichment Plot for All Scoring Functions Figure 4 represents the enrichment of all scoring functions for the subset of 3662 molecules returned from RP that was used in the analysis. The x-axis represents the rank order of all compounds in the RP dataset. The y-axis is the number of actives in the hit list. The plot in Figure 4 indicates that two scoring functionsLS1 and PMF performed well above random (percent actives returned >> percent of the database screened). The remaining scores were all close to the random line. With this observation, we calculated a consensus score using only LS1 and PMF. Consensus scoring in Cerius allows for a user-defined cutoff based on the rank order of all poses in the dataset. In other words, for a single scoring function, a user-defined cutoff of 10% indicates that a pose will be assigned a value of 1 if the overall rank of that pose is in the top 10% of all poses. When multiple scoring functions are selected, the consensus score is a simple sum of each score, when the percent cutoff is satisfied. In our work here, we calculated the consensus score based on three separate cutoff values10%, 20%, and 30%using both the individual poses and the best pose only (BPO) option for a molecule. The maximum consensus score for all options is a value of 2: top X% for both LS1 and PMF. 2 Pharmaceuticals case study continued All scores calculated in this work are shown in Figure 5. Page 7 of 11 Figure 5: Statistical Results from All Scoring Functions In Figure 5: • The Method column is the score used in the analysis. Consensus scores are labeled as ‘CON’ followed by the percent cutoff (e.g. CON10 is a Consensus score with a 10% cutoff). BPO stands for best pose only calculation of the consensus score. • % Cov Act is the percent coverage of actives in the dataset. Since a single scoring function can provide a complete rank order of all molecules, we chose to use the 3 layers here50% coverage (16 out of 32 actives), 80% coverage (25/32) and 100% (32/32). However, for consensus scoring, we are not able to get a full rank order list; the consensus scores gave us categorical results (consensus score of 2, 1, or 0). In this figure, all statistics are based on consensus score of 2. • Rank is the highest rank value (largest number) for the score. • % DB Sc is the percentage of the database screenedthe RP output set of 3662. • Enrich is the enrichment value for the RP output set. • GH is the Goodness of Hit (Güner-Henry) score for the RP output set. • % Full DB Sc is the percentage of the full McMaster database screened 49,995 molecules. • Full DB E is the enrichment value for the full McMaster database. • Full DB GH is the Goodness of Hit score for the full McMaster database. • The table entries are ranked by the Full DB Enrichment column. Review of the statistics in Figure 5 reveals some interesting observations. First, all of the consensus scores performed well when ranked by enrichment. When considering the full dataset, a minimum 20-fold enrichment was observed, with the maximum enrichment value at over 60-fold. This represents a 2- to6-fold increase in enrichment relative to RP alone (10-fold). Pharmaceuticals case study continued Second, LS1 and PMF also ranked high in the list of scores, especially at the 50% coverage level. These statistics were, of course, expected after visualizing the enrichment plot in Figure 4. One can see that LS1 and PMF have a steep increase in the number of actives ranked high in the output list before reaching a threshold. Third,none of the scores were able to extract all of the known actives (100% Cov Act). Nearly the entire set of 3662 molecules was needed to rank the full set of 32 actives high in the results. This is not surprising as this is a very restrictive cutoff value for in Page 8 of 11 silico screening, and is generally considered beyond the scope of such models. Retrieving all 32 actives out in an experiment would be quite challenging for any technique. Fourth,while enrichment is a valuable metric for evaluating the quality of results,we feel it is very important to also consider the coverage of actives in the dataset. The enrichment is a measure of the selectivity (high percentage of actives covered in a low percentage of the full dataset covered). We can observe this on our results by viewing the CON10BPO statistics in Figure 5. As illustrated in the figure, this score resulted in a very small percentage of molecules (0.1%) with a consensus score of 2 (45 out of 49,995); however, the coverage of actives was also very small (3% or 1 out of 32). As such, for this work, we felt strongly that it was important to find a good intersection between the coverage of actives and the enrichment. By inspection of the results,we determined that the CON10 set was most desirable. In this result,we see the following: 1) Nearly half of the known actives are retrieved (15 of the 32 actives —47%). 2) Only 1% of the full database was screened (484 out of 49,995). These results far exceed a random screening approach, where one would expect percentages for the actives retrieved and the database screened to be roughly equal. While other results may have had a higher enrichment score (or less molecules in the output), the determining factor was that the coverage of actives was not high enough, in our opinion. It certainly could be debated that a more selective result, such as CON20BPO, would yield a better hit rate. However, a loss of 81% of the actives (19% coverage) was too great of a penalty and this must be considered in the final analysis,when determining the best protocol to use for this type of vHTS workflow. Of the 484 molecules with a consensus score of 2 from the CON10 calculation, there were 15 known actives in the output: MAC-0002330 MAC-0002600 MAC-0002654 MAC-0002655 MAC-0010097 Pharmaceuticals case study continued MAC-0012674 MAC-0012675 MAC-0022248 MAC-0030856 MAC-0030857 MAC-0030861 MAC-0030862 MAC-0030863 MAC-0037144 MAC-0038969 The range of scores was 3.8 to 4.9 for LSI and 10.6 to 30.3 for PMF. The remaining 469 inactives in the CON10 results had a range of scores of 3.7 to 6.0 for LS1 and 10.6 to 60.4 for PMF. Of course, each of the scores was in the 10th percentile of the overall LS1 and PMF scores for the RP dataset. Figure 6 shows the CON10 pose for two of the known actives. Figure 6: Poses of Two Known Actives, MAC-0002600 (top) and MAC-0030857 (bottom,), with CON10 Scores of 2 Page 9 of 11 Figure 6 illustrates the differences in poses of two more structurally dissimilar molecules (by visual inspection only) from the 15 actives that had a CON10 score of 2. The red molecule is the X-ray bound ligand (folate) for 1RA2 placed here only as a reference. As one may expect, both of the molecules extend into the region occupied by both the folate and the open cavity left from removal of the cofactor (near TYR100). To reiterate, with the knowledge that the collection of actives were both competitive and non-competitive inhibitors, we opted to remove the cofactor to allow for additional interactions to occur with the protein. This decision turned out to be quite fruitful for this project due to this result: Of the 15 known actives in the CON10 set, eight were competitive inhibitors, which represents a near 50/50 split between competitive and noncompetitive inhibitors. Therefore, the sequential screening protocol developed here not only produced a substantial 48-fold enrichment, but it also retrieved compounds with multiple modes of interaction to inhibit DHFR. Pharmaceuticals case study continued Summary The McMaster University virtual screening competition put forth an interesting challenge: developing in silico techniques that could be used to identify a database of 32 actives out of a full collection of 50,000 molecules. Using several tools from Accelrys Software Inc, we employed a sequential screening workflow using recursive partitioning and structure based drug design that greatly reduced the complexity of the problem. If we were to employ a standard screening approach, we would expect a hit rate of only 0.06% (32/50,000). However, using the in silico approach presented here, the hit rate was over 3% (15/484). Another way to illustrate the utility of this procedure is by highlighting this fact: Using standard or random screening, one would expect to screen 50% of the dataset in order to retrieve 50% of the actives. However, using the procedure outlined here, we would expect to retrieve nearly half of the actives by only screening 1% of the full collection. Clearly, this rational approach based on computational techniques and the series of steps presented here can quickly focus on a small percentage of the dataset predicted to be active, thereby greatly reducing the time required to biologically screen this compound collection. Page 10 of 11 References 1. McMaster University High Throughput Screening Lab, Hamilton, Ontario, Canada, http://hts.mcmaster.ca/HTSDataMiningCompetition.htm. 2. Zolli-Junan, M., Cechetto, J. D., Hartlen, R, Diagle, D. M., Brown, E. D., ‘High throughput screening identifies novel inhibitors of Escherichia coli dihydrofolate reductase that are competitive with dihydrofolate,’ Bioorg. Med. Chem. Lett., 2003, 13, 2493-2496. 3. Varma, S. et. al. ‘Creating a Smart Virtual Screening Protocol II: Recursive Partitioning for Sequential Screening,’ http://www.accelrys.com/cases/ smart_virtual_screening_.pdf. 4. Lyons, T. et. al. ‘Creating a Smart Virtual Screening Protocol I: Preparing the Target Protein,’ http://www.accelrys.com/cases. 5. Venkatachalam, C.M., Jiang, X., Oldfield, T., Waldman, M., ‘LigandFit: A Novel Method for the Shape-Directed Rapid Docking of Ligands to Protein Active 2 Sites,’ J Mol Graph Model, 2003, 21, 289-307. Cerius ·LigandFit Software; Accelrys Software Inc, http://www.accelrys.com/cerius2/c2ligandfit.html. 6. Diller, D. J., Merz, K. M., ‘High throughput docking for library design and library prioritization,’ Proteins, 2001, 43, 113-124. 7. Accord Entreprise Informatics Software, Accelrys Software Inc, http://www.accelrys.com/aei/. Pharmaceuticals case study continued 8. TSAR Software, Accelrys Software Inc, http://www.accelrys.com /products/tsar/corina.html. 9. Catalyst Software, Accelrys Software Inc, http://www.accelrys.com/catalyst/. 10. Perl script written By Amit Kulkarni, Accelrys Software Inc. http://www.perl.org. 11. OpenBabel Program, SourceForge, http://openbabel.sourceforge.net/. 12. Gehlhaar, D.K., Moeerder, K.E., Zichi, D., Sherman, C.J., Ogden, R.C., Freer, S.T., ‘De Novo Design of Enzyme Inhibitors by Monte Carlo Ligand Generation,’ J. Page 11 of 11 Med. Chem., 1995, 38, 466-472. 13. Gehlhaar, D.K.; Bouzida, D.; Rejto, P.A. 1999. Rational Drug Design: Novel Methodology and Practical Applications. L. Parrill, M. Rami Reddy, [editors]. Washington, DC: American Chemical Society, 1999, Series title: ACS symposium series 719, 292-311. 14. Muegge, I., Martin, Y.C., ‘A General and Fast Scoring Function for ProteinLigand Interactions: A Simplified Potential Approach,’ J. Med. Chem. 1999, 42, 791-804. 15. Jain, A.N., ‘Scoring noncovalent protein-ligand interactions: A continuous differentiable function tuned to compute binding affinities,’ J. Comp. Aided Molec. Design, 1996, 10, 427-440. 16. DS Accord for Excel program, Accelrys Software Inc, http://www.accelrys.com/dstudio/ds_accord/excel.html. 17. Güner, O. F.; Henry, D. R., In Pharmacophore Perception, Development, and Use in Drug Design; Güner, O. F. ed., International University Line: La Jolla, CA, 2000, p. 191. Smart Virtual Screening Protocol, Part III: Using Structure Based Drug Design to Enhance an In Silico Workflow.