Protein Flexibility in Ligand Docking and Virtual Screening to Protein Kinases

doi:10.1016/j.jmb.2004.01.003

J. Mol. Biol. (2004) 337, 209–225

Protein Flexibility in Ligand Docking and Virtual Screening to Protein Kinases Claudio N. Cavasotto* and Ruben A. Abagyan Molsoft LLC, 3366 N Torrey Pines Ct. Suite 300, La Jolla CA 92037, USA

The main complicating factor in structure-based drug design is receptor rearrangement upon ligand binding (induced fit). It is the induced fit that complicates cross-docking of ligands from different ligand –receptor complexes. Previous studies have shown the necessity to include protein flexibility in ligand docking and virtual screening. Very few docking methods have been developed to predict the induced fit reliably and, at the same time, to improve on discriminating between binders and nonbinders in the virtual screening process. We present an algorithm called the ICM-flexible receptor docking algorithm (IFREDA) to account for protein flexibility in virtual screening. By docking flexible ligands to a flexible receptor, IFREDA generates a discrete set of receptor conformations, which are then used to perform flexible ligand –rigid receptor docking and scoring. This is followed by a merging and shrinking step, where the results of the multiple virtual screenings are condensed to improve the enrichment factor. In the IFREDA approach, both side-chain rearrangements and essential backbone movements are taken into consideration, thus sampling adequately the conformational space of the receptor, even in cases of large loop movements. As a preliminary step, to show the importance of incorporating protein flexibility in ligand docking and virtual screening, and to validate the merging and shrinking procedure, we compiled an extensive small-scale virtual screening benchmark of 33 crystal structures of four different protein kinases sub-families (cAPK, CDK-2, P38 and LCK), where we obtained an enrichment factor fold-increase of 1.85 ^ 0.65 using two or three multiple experimental conformations. IFREDA was used in eight protein kinase complexes and was able to find the correct ligand conformation and discriminate the correct conformations from the “misdocked” conformations solely on the basis of energy calculation. Five of the generated structures were used in the small-scale virtual screening stage and, by merging and shrinking the results with those of the original structure, we show an enrichment factor fold increase of 1.89 ^ 0.60, comparable to that obtained using multiple experimental conformations. Our cross-docking tests on the protein kinase benchmark underscore the necessity of incorporating protein flexibility in both ligand docking and virtual screening. The methodology presented here will be extremely useful in cases where few or no experimental structures of complexes are available, while some binders are known. q 2004 Elsevier Ltd. All rights reserved.

*Corresponding author

Keywords: protein flexibility and induced-fit; ligand docking; structurebased drug design; virtual screening; protein kinases

Supplementary data associated with this article can be found at doi: 10.1016/j.jmb.2004.01.003 Abbreviations used: BPMC, biased probability Monte Carlo; DA, docking accuracy; EF, enrichment factor; GB/SA, generalized Born/surface area; ICM, internal coordinate mechanics; IFREDA, ICM-flexible receptor docking algorithm; LBP, ligand-binding pocket; MD, molecular dynamics; PK, protein kinase; RMSD, root-mean-square deviation; VS, virtual screening. E-mail address of the corresponding author: [email protected] 0022-2836/$ - see front matter q 2004 Elsevier Ltd. All rights reserved.

210

Introduction The rapid progress of genomics explosion will result in a dramatic increase of novel yet biologically validated targets for drug discovery. Structure-based drug design is now established as a key first step in the lengthy process of developing new drugs.1 Thus, the role of computer-aided drug design through virtual screening (VS) of available or virtual chemical libraries will thus continue to grow.2 – 5 However, advances in this technology are badly needed to improve the accuracy of the predicted geometries and scores. Induced molecular flexibility is fundamental to understanding the principles of molecular recognition between ligand and receptor. Upon ligand binding, many systems undergo rearrangements, which range from local motions of side-chains to large domain movements. In any case, receptor flexibility might have a dramatic impact in the ligand docking problem and VS. It has been shown that even small changes in the receptor conformation can be important in computing binding affinities.6 The importance of receptor flexibility and its implication in drug discovery has been highlighted in two excellent reviews in this field,7,8 and everything points in the direction that protein mobility will have an increasing role in computeraided drug design in the future. Dealing with protein flexibility is essential to predict the orientation and interactions of a ligand within a binding pocket in the absence of experimental structural information. Prediction of mutation resistance to drugs can benefit from reliable docking algorithms that include conformational sampling of the receptor. There have been a few reports concerning the impact of protein flexibility in ligand docking to different protein families. Considering two inhibitors of HIVp, Bouzida et al.9 demonstrated convincingly the limitation introduced by considering a single and rigid receptor structure. An analysis of the sensitivity of the docking results to protein flexibility in thrombin, thermolysis and neuraminidase6 showed that only 49% of the ligands are cross-docked correctly to a receptor structure bound to a different ligand, while small movements in the receptor structure can lead to errors up to 14 kJ/mol in the binding energy prediction. The authors pointed out that sidechain flexibility is not sufficient for accounting for the mis-docking of inhibitors. While evaluating different docking methods on three different receptors,10 the authors showed that lead docking to a single receptor conformation significantly reduces the chances of finding the correct pose. There have been some attempts in the past to include protein flexibility in the ligand docking procedure. These include the early attempts using soft docking,11 partial side-chain flexibility,12,13 continuous side-chain sampling14 and rotameric libraries.15,16 The hinge-bending concept was also used to model receptor flexibility.17,18 Although

Protein Flexibility in Docking and Virtual Screening

incorporating side-chain flexibility was a big step forward, current methods should go beyond this point to include backbone rearrangements. As has been pointed out recently,19 the use of several receptor structures seems to be the best choice to date to incorporate flexibility in the docking problem, but many questions arise. What should be the source of these structures? How many are needed? How should they be used or results be combined? And, more importantly, how should this flexibility be incorporated in the VS procedure? To date, an explicit and direct consideration of the receptor plasticity in the VS procedure is still computationally unattainable. Different approaches on how to use an ensemble of receptor structures in ligand docking have been presented. Knegtel et al. used NMR and crystal structures to generate combined interaction grids by averaging with respect to energy and geometry.20 These composite grids were used with a rigid ligand approach to re-dock native ligands and to identify known ligands from a small compound database. FlexE incorporates flexibility through discrete alternative conformations of varying parts of the protein taken from structures that have a very similar backbone trace, which are merged in a combinatorial way and considered directly during the rigid ligand docking.21 The method is evaluated for root-mean-square deviations (RMSD) on ten different proteins containing 105 crystal structures and 60 different ligands. When the top ten solutions for each ligand are con˚ of sidered, FlexE finds the ligands within 2.0 A their native pose in 67% of the cases. An analysis of four different choices for combining many structures into a single representative energy grid has been performed recently.22 The authors used 21 crystal structures of the HIV-1 protease with diverse inhibitors. These inhibitors were docked into the generated combined grids using AutoDock23 – 25 and the RMSD from native pose and the binding energy were evaluated and it was found that the weight-averaged grids perform best. Frimurer et al. used a rotamer library of four key residues to improve predictions of binding geometry and affinities in protein tyrosine phosphatase 1B (PTP1B).26 The library was built on the basis of observations on three crystal structures of PTP1B bound to different inhibitors and resulted in 96 models. Docking of the three inhibitors to these models improved their geometry and binding energy predictions. Side-chain flexibility has also been included in SLIDE, which uses a set of template points constructed for hydrogen bond donors and acceptors, and hydrophobic regions to represent the binding pocket.27 Molecular dynamics (MD) simulations have been used to generate an ensemble of different receptor conformations as input for the generation of a composite interaction weight-averaged grid.28 This method was applied to VS against

211

Protein Flexibility in Docking and Virtual Screening

dihydrofolate reductase and found improvements in the top-ranked 10% of a database of drug-like molecules. Multiple structures generated through MD simulations were used to build a receptorbased pharmacophore model for the HIV-1 integrase.29 In the search for the correct ligand – receptor conformation, the relaxed complex method30,31 uses an ensemble of structures generated through MD simulations of the unliganded receptor to dock a mini-library of binders using a fast rigid receptor docking method. This rapid docking is used as a filter and selected complex conformations are rescored with a more accurate energy function using the molecular mechanics/Poisson – Boltzmann surface area approach. The observed experimental complexes are found within the lowest free energy complexes. Recently, we have used a continuous stochastic global optimization method using the internal coordinate mechanics (ICM)14 methodology to incorporate side-chain and backbone flexibility to analyze the structural binding determinants of an RXR antagonist to the receptor (C.N.C. et al., unpublished results). With a similar methodology, the binding mode of the ligand and pocket sidechain conformations could be predicted from a random starting conformation, with and without the ligand present, in the seven-transmembrane proteins rhodopsin and bacteriorhodopsin.32 The ˚ RMSD and the ligand was predicted within 0.2 A ˚ . In a RMSD for the pocket side-chains was 0.3 A recent article, a procedure was presented to include side-chain receptor flexibility and continuum solvation in flexible ligand docking.33 Using a Monte Carlo simulation and the generalized Born/surface area (GB/SA) continuum solvent model, the authors used a test set of 14 complexes to evaluate ligand docking with and without sidechain flexibility. Although the RMSD values of the ligand were comparable in both approaches, energy discrimination of the binding mode actually deteriorated in the flexible receptor docking procedure. The authors concluded that including protein flexibility may result in a rugged energy landscape with less distinguishable multiple minima. Despite previous studies, to our knowledge no systematic attempt has been undertaken to examine and incorporate the influence of protein flexibility on docking geometries and on VS and enrichment factors (EFs). In a real-life VS experiment, RMSD values cannot be computed and we are left with a collection of binding scores of the docked compound library. We present a novel algorithm called the ICMflexible receptor docking algorithm (IFREDA), which incorporates protein flexibility in ligand docking and VS, especially in the cases where multiple experimental structures representative of the conformational space of the target protein are not available. IFREDA generates an ensemble of receptor conformations by performing flexible ligand docking of selected known binders to a

flexible receptor. The conformational ensemble thus generated is then used to perform flexible ligand – rigid receptor docking and scoring, and results from the multiple VS are then condensed using a merging and shrinking procedure. IFREDA accounts for both side-chains and key backbone movements, thus sampling the conformational space of the target receptor. IFREDA was used in eight protein kinase (PK) complexes and was able to identify the correct ligand pose as the bestenergy ranking conformation. By merging and shrinking the results from VS against the generated structures, the average EF fold increase was comparable to that obtained using multiple experimental structures (, 1.9). To underscore the necessity of incorporating protein flexibility in ligand docking and VS, and to validate our merging and shrinking procedure to improve EFs using multiple experimental structures, we report extensive flexible ligand – grid receptor docking and small-scale VS tests against 33 crystal structures of four different PK subfamilies. PKs have been implicated in proliferation, invasion and metastasis of many types of cancer (the importance of PKs as relevant drug targets and the progress in developing new PK inhibitors have been reviewed recently34). Induced-fit effects in PKs may be one of the factors that explain why the common ATP-binding sites are good drug targets.7 The accuracy of native ligand docking and cross-docking, and the impact of receptor flexibility in ligand docking and scoring are reported, together with the validation of the merging and shrinking procedure to condense the VS results of multiple receptor conformations to improve the EFs. In Results, we report the analysis on the PK benchmark of 33 crystal structures and the use of IFREDA on eight PK complexes. This is followed by Discussion, where we address the best strategy for the choice of multiple receptor conformations, and the limitations and further improvements of our algorithm. A summary of our results and achievements followed by a detailed description of the Methods used are then presented.

Results The 100% accurate docking of native PK ligands The PK family is a difficult test set for any docking-scoring method. Upon ligand binding, the side-chains in the binding pocket may adopt different conformational states and, in some cases, loop rearrangements are observed (induced-fit). A part of the ligand is usually solvent-exposed and held in position by many hydrophobic contacts. It has been pointed out that induced-fit effects in PKs together with hydrophobicity of the binding pocket may explain why the common ATP-binding sites are good drug targets.7 Also, some ligand – protein interactions are mediated by water

212

Protein Flexibility in Docking and Virtual Screening

Table 1. Protein kinase complexes used in ligand docking and virtual screening

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

PDB entry

Ligand code

Ligand name

1BKX 1BX6 1FMO 1STC 1YDR 1YDS 1YDT 1JLU 1AQ1 1DI8 1DM2 1E1X 1E9H 1FVT 1FVV 1G5S 1H1P

adn ba1 adn sto iqp iqs iqb

Adenosine Balanol Adenosine Staurosporine H7 H8 H89

stu dtq hmd nw1 inr 106 107 i17 cmg

Staurosporine Anilinoquinazoline 2 Hymenialdisine NU6027 E226 Oxindole 16 Oxindole 91 H717 NU2058

Kinase family cAPK cAPK cAPK cAPK cAPK cAPK cAPK cAPK CDK2 CDK2 CDK2 CDK2 CDK2 CDK2 CDK2 CDK2 CDK2

molecules. The four protein kinase sub-families used in this study were: three serine/threonine PKs, the cAMP-dependent protein kinase (cAPK), the cyclin-dependent kinase 2 (CDK-2), the mitogen-activated protein P38, and one tyrosine PK, the lymphocyte-specific kinase (LCK). They were chosen because there are more than four crystal structures of complexes with diverse compounds (see Table 1 for the crystal structures used in this study). The ligand of 1QPC (AMP-PNP) was not used for ligand docking, since its disordered g-phosphate group complicates the determination of comparable RMSD values. However, the receptor structure of 1QPC was used for cross-docking and VS. The first step in addressing the problem of structural flexibility of PKs was to perform a small-scale ligand docking and VS against the PK receptors. A 1000 compound library of random molecules seeded with the corresponding PK co-crystallized ligands (see Materials and Methods for details) was screened against each of the receptors (see Table 1) using the ICM flexible ligand –grid receptor docking algorithm (see Materials and Methods).35 – 38 During the energy optimization of the ligand in the field of the receptor, a conformational set of low-energy states is generated and the best-energy conformation is scored. In this way, the reported RMSD values for ligands in flexible ligand– rigid receptor docking always refer to the best energy-ranked solution. The screening is repeated four times and the best score of each ligand is kept. The native ligands were docked to the PK structures with remarkable accuracy. As it is shown in ˚ and 100% of Table 2, the average RMSD is 0.74 A the native ligands in the 29 holo complexes are ˚ RMSD (individual RMSD docked within 1.5 A values for each receptor are detailed in Figures 1 and 2; and see the Supplementary Material). RMSD values are always calculated between heavy atoms of the docked ligand with the com-

18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

PDB entry

Ligand code

Ligand name

Kinase Family

1H1Q 1H1S 1JSV 1HCL 1A9U 1BL6 1BL7 1DI9 1BMK 1M7Q 1P38 1QPC 1QPD 1QPE 1QPJ 3LCK

2a6 4sp u55

NU6094 NU6102 PNU112455A

sb2 sb6 sb4 msq sb5 dqo

SB203580 SB216695 SB220025 Anilinoquinazoline 3 SB218655 14e

anp stu pp2 stu

ANP-PNP Staurosporine PP2 Staurosporine

CDK2 CDK2 CDK2 CDK2 P38 P38 P38 P38 P38 P38 P38 LCK LCK LCK LCK LCK

plex in the crystal structure, after superposition of the backbone atoms within the ligand-binding pocket. When symmetric moieties in the molecules are present, all possible atom numberings within the symmetric portion are generated and that with the lowest RMSD is kept. Evaluating the strong impact of receptor flexibility in ligand docking and scoring The influence of the induced conformational changes on docking results becomes apparent when a ligand is docked to a receptor complexed with another compound. Since, sometimes, small variations in the structure of the binding pocket might have a large impact on docking geometries, cross-docking experiments are useful to assess the magnitude of this influence. As it is shown in Table 2, on average, only 70% of the ligands are cross-docked correctly (apo structures were excluded from this calculation), in agreement with a recent study on docking strategies, where ICM was ranked first for pose prediction among the Table 2. RMS deviation for the docking of native ligands and docking accuracy in cross-docking experiments Complete set

cAPK

CDK2

P38

LCK

0.71

0.84

0.77

0.63

Docking of native ligands 76.0 fRMSD,1:0 A a fRMSD,1:5 A a

˚) Average RMSD (A Cross-docking Average DAb

100.0 0.74 70.0

65.0

65.0

72.0

88.0

RMSD values refer to the best energy docking solution. Fraction of ligands with RMSD values below the indicated ˚ ). threshold (A b Docking accuracy (DA) calculated as: DA ¼ fRMSD,2:0 A þ 0:5ð fRMSD,3:0 A 2 fRMSD,2:0 A Þ: a

213

Protein Flexibility in Docking and Virtual Screening

˚ ) for heavy atoms is calculated after superposition of the backbone Figure 1. Cross-docking on CDK-2. RMSD (A atoms within the ligand-binding pocket. RMSD values refer to the best energy solution. Thick border cells indicate native complex. DA, docking accuracy. pThe disordered and solvent-exposed benzyl moiety of ligand 107 was excluded from the RMSDs calculations.

other three docking methods.10 The RMSD values of ligand cross-docking for each PK sub-family are shown in Table 2 (detailed cross-docking matrices are detailed in Figures 1 and 2; and see the Supplementary Material). In a real-life case, the goal of VS is to select a small number (, 2%) of potential binders to the receptor of interest from a large source library. For a particular focused library built with the top s compounds of the ranked database, the enrichment

˚ ) for heavy Figure 2. Cross-docking on cAPK. RMSD (A atoms is calculated after superposition of the backbone atoms within the ligand-binding pocket. RMSD values refer to the best energy solution. Thick border cells indicate native complex. DA, docking accuracy.

factor can be calculated as: EFðsÞ ¼

Hitss Hitstotal = NCs NCtotal

where NC is the number of compounds. It is evident, however, that pose prediction, while being a necessary component of the VS procedure, is not sufficient for accurate scoring and thus high EFs. In the small-scale VS of the 1000 compound library, about 80% of the native binders are scored in the top 1.5% of the screened database when docked to their co-crystallized structure, showing the accuracy of our scoring function (throughout this work we did not try to optimize or improve the scoring function specifically for PKs). In order to examine the relation between RMSD values and ranking, and to guide us to assess the impact of structural diversity on ligand docking and scoring, the correlation between compound scoring and RMSD deviation is plotted in Figure 3. Only 4% of the points are in the upper-left corner (bad-pose, good-score), and , 50% of these are close to the ˚ . The common threshold for a borderline of 2.5 A compound to be considered docked correctly is ˚ . However, due to variations in the position ,2 A of the backbone when overlaying the structures for RMSD calculation, we preferred to use a ˚ . In the lower-right corner threshold value of 2.5 A we have the good-pose, bad-score compounds (17%). The fraction of ligands that are docked correctly

214

Figure 3. RMS deviations versus ranking of protein kinase ligands. Colored marks represent outliers: yellow circle, cAPK; blue triangles, CDK-2; red squares, P38. No outlier was found for LCK. The black lines delimit ˚ RMSD deviations. the areas of 10% ranking and 2.5 A RMSD values refer to the best energy docking solution.

and ranked within the top 10% is , 49% (lower-left corner). If we compare this number with the , 70% of compounds that are cross-docked correctly, we see that the PK scoring and ranking is more sensitive to the induced-fit effects than ligand docking. Still, both 49% of good ranking compounds and , 70% of correct geometries is indeed a good result. In the next section we improve the performance in both the geometry prediction and the EFs by using many experimental receptor conformations. How to merge the screening results of known multiple receptor conformations to improve the enrichment factor? Inspection of the cross-docking RMSD values and EFs can help to determine which crystal structure will be used or if it is meaningful to do the full VS against more than one receptor structure to better represent conformational flexibility. For example, there are two diverse cAPK structures with adenosine bound (PDB codes 1BKX and 1FMO). No distinction as to which to choose can be made from the cross-docking RMSD, but inspection of the EFs shows that 1FMO might be a better target for VS (16.7 versus 0.0 for the top 2% screened database), partially linked probably to its ˚ versus 2.6 A ˚ ). CDK-2 better resolution (2.2 A structure 1H1S links good EF (58.3 for the top 1% selection) with high docking accuracy (DA). When two or more crystal structures are available, small-scale VS including non-native known binders could be very useful. In this case, results of the screening against each conformation are merged and the best rank for each compound is kept, thus shrinking the scoring list to the size corresponding to a VS against a single receptor

Protein Flexibility in Docking and Virtual Screening

(referred to here as the merging-shrinking procedure). We show in Table 3 how the combination of screening results against two and three different PK structures can yield better RMSD and EF values than each one alone. Note that the better performance is not self evident because the mergingshrinking procedure can dilute the correct results. Within each of the PK sub-families, a few representatives with chemically different ligands were selected (Table 3). The LCK sub-family members were not included, since they had only two different ligands. In addition, the 1D19 receptor was excluded from the P38 sub-family, since the initial EF was zero for 1% or 2% and by calculating the EF increase would lead to a singularity. To characterize the improvement in EF through merging and shrinking the screening results, we introduce the EF fold-increase, calculated as: EF fold-increase ¼

EFmerged EFinvididual

The average EF fold increase was 1.85 ^ 0.65 calculated for the 21 groups of Table 3 using 141 EF foldincrease values. The EF increase considering the top 1%, 2% and 10% of the screened database was 1.89 ^ 0.60, 1.83 ^ 0.60 and 1.83 ^ 0.74, respectively, which shows that the EF fold-increase was roughly independent of the size of the focused selection. The EF fold-increase was less than one in only four out of the 141 cases (, 3%), which highlights the advantage of the merging-shrinking technique. In those cases, the EF drop occurred because of the high EF values for structures 1H1S and 1YDS, which were diluted with inferior EF values of the other receptor conformations. The distribution of the EF fold-increase values is shown in Figure 4. It is evident that combining even two structures is enough to have most of the ligands docked correctly (in 19 out of the 21 groups more than 75% of the ligands are docked correctly when the merged results are considered). The mergingshrinking procedure has a more significant impact on the RMSD values than on EFs, showing again that scoring is more sensitive to protein flexibility than ligand docking. However, an average EF fold increase of , 2 should be always considered as significant. Since the computing time for the flexible ligand – rigid receptor docking and scoring is , 1 – 2 minutes using a 700 MHz processor (1 Mb RAM for the dual-processor node), screening against a few structures is affordable and has the advantage of using actual receptor conformations, which may differ in side-chain conformation, and in the rearrangement of loops. The ICM-flexible receptor docking algorithm (IFREDA) When several different crystal structures of the same receptor are available, a small-scale VS

215

Protein Flexibility in Docking and Virtual Screening

Table 3. Improvements in enrichment factors and RMSD values for focused libraries by using the merging and shrinking procedure

PDB entry

EF (top 1%)

EF (top 2%)

EF (top 10%)

%Binders with ˚ RMSD , 2.5 A

1DM2 1G5S 1DM2 þ 1G5S

25.0 41.7 58.3

16.7 25.0 29.2

3.3 5.8 6.6

41.7 91.7 91.7

1DM2 1JSV 1DM2 þ 1JSV

25.0 33.3 50.0

16.7 20.8 29.2

3.3 5.0 7.5

41.7 41.7 75.0

1DM2 1H1P 1DM2 þ 1H1P

25.0 16.7 41.7

16.7 8.3 20.8

3.3 4.2 5.8

41.7 58.3 83.3

1DM2 1G5S 1JSV 1DM2 þ 1G5S þ 1JSV

25.0 41.7 33.3 66.7

16.7 25.0 20.8 37.5

3.3 5.8 5.0 9.2

41.7 91.7 41.7 100.0

1DM2 1H1P 1JSV 1DM2 þ 1H1P þ 1JSV

25.0 16.7 33.3 50.0

16.7 8.3 20.8 25.0

3.3 4.2 5.0 7.5

41.7 58.3 41.7 100.0

1AQ1 1E1X 1AQ1 þ 1E1X

41.7 33.3 50.0

29.2 20.8 33.3

6.7 4.2 7.5

75.0 25.0 91.7

1E1X 1G5S 1E1X þ 1G5S

33.3 41.7 58.3

20.8 25.0 33.3

4.2 5.8 6.7

25.0 91.7 91.7

1E1X 1G5S 1DM2 1E1X þ 1G5S þ 1DM2

33.3 41.7 25.0 66.7

20.8 25.0 16.7 33.3

4.2 5.8 3.3 8.3

25.0 91.7 41.7 91.7

1FVT 1E1X 1FVT þ 1E1X

25.0 33.3 50.0

16.7 20.8 29.2

5.0 4.2 6.7

58.3 25.0 66.7

1H1S 1E1X 1H1S þ 1E1X

58.3 33.3 66.7

29.2 20.8 37.5

9.2 4.2 7.5

100.0 25.0 100.0

1DM2 1H1P 1DM2 þ 1H1P

25.0 16.7 41.7

16.7 8.3 20.8

3.3 4.2 5.8

41.7 58.3 83.3

1FMO 1YDS 1FMO þ 1YDS

16.7 50.0 50.0

16.7 25.0 33.3

5.0 5.0 8.3

50.0 83.3 83.3

1FMO 1BX6 1FMO þ 1BX6

16.7 16.7 33.3

16.7 8.33 16.7

5.0 1.67 6.67

50.0 50.0 66.7

1BX6 1YDS 1BX6 þ 1YDS

16.7 50.0 33.3

8.33 25.0 16.7

1.67 5.0 5.0

50.0 83.3 83.3

1FMO 1STC 1FMO þ 1STC

16.7 16.7 33.3

16.7 8.33 16.7

5.0 5.0 8.33

50.0 83.3 83.3

1BX6 1STC 1BX6 þ 1STC

16.7 16.7 33.3

8.33 8.33 16.7

1.67 5.0 3.33

50.0 83.3 83.3

1FMO 1BX6 1STC 1FMO þ 1BX6 þ 1STC

16.7 16.7 16.7 50.0

16.7 8.33 8.33 25.0

5.0 1.67 5.0 6.67

50.0 50.0 83.3 83.3

1A9U 1BMK

16.7 16.7

8.33 8.33

5.0 3.33

83.3 66.7 (continued)

Table 3 Continued

PDB entry

EF (top 1%)

EF (top 2%)

EF (top 10%)

%Binders with ˚ RMSD , 2.5 A

1A9U þ 1BMK

16.7

16.7

8.33

83.3

1A9U 1M7Q 1A9U þ 1M7Q

16.7 16.7 33.3

8.33 8.33 16.7

5.0 1.67 5.0

83.3 83.3 83.3

1BMK 1M7Q 1BMK þ 1M7Q

16.7 16.7 16.7

8.33 8.33 16.7

3.33 1.67 5.0

66.7 83.3 83.3

1A9U 1M7Q 1BMK 1A9U þ 1M7Q þ 1BMK

16.7 16.7 16.7 33.3

8.33 8.33 8.33 25.0

5.0 1.67 3.33 5.0

83.3 83.3 66.7 83.3

Virtual screening is performed against two and three multiple experimental crystal structures of the same protein kinase subfamily.

followed by analysis of cross-docking RMSDs and EFs can help to determine how to include potential structural diversity of the receptor-binding pocket in the VS using large compound libraries as discussed above. Merging and shrinking the results of two or three screenings could be a solution to improve the EF. The IFREDA is especially useful in cases where only one holo (or apo) crystal structure and known binders are available. It has three main steps: (i) de novo receptor structure generation by performing flexible ligand docking of known binders to a flexible receptor, where side-chain and essential backbone movements are taken into consideration; (ii) VS against the generated conformations using flexible ligand –grid receptor docking; (iii) merging and shrinking of the VS results in the same way as described above. The computational procedure to de novo generation of alternative multiple receptor conformations This procedure (see Materials and Methods for details) consists of: Seeding. The ligand was placed within the ligand-binding pocket and four conformations are generated by flipping the ligand 1808 with respect to its principal axes of inertia. Each of these starting poses is then followed by ten random displacements and rigid rotations of the ligand, thus generating 40 starting complex conformations. Soft-van der Waals structure relaxation. Each complex conformation was allowed to relax through in vacuo minimization using a variable weight soft van der Waals algorithm. Selected backbone loops, pocket side-chains and the ligand were regarded as fully flexible in this stage. Stochastic global energy optimization. The 40 complexes generated through seeding and in vacuo minimization were subjected to a stochastic global energy optimization of the side-chain and ligand

216

Protein Flexibility in Docking and Virtual Screening

Figure 4. Distribution of the enrichment factor fold-increase obtained by the merging-shrinking procedure on 21 groups of two or three protein kinases (Table 3).

torsion angles using the double-energy minimization scheme.14 During simulations, a set of geometrically diverse low-energy states are stored, which are then clustered by comparing the RMSD of the heavy-atom coordinates of the ligand to eliminate redundant conformations. The top-ranking conformations were subjected to full minimization, keeping flexible the selected backbone loops, pocket side-chains and the ligand. The energy was then re-evaluated with a more accurate solvation energy term by solving the Poisson equation using the boundary element algorithm.39

because of ligand –receptor clashes (induced-fit upon binding) or because non-optimal ligand – receptor contacts, which lead to poor binding scores. The complexes chosen were A, 1FMO þ balanol; B, 1JLU þ staurosporine; C, 1JLU þ H7; D, 1DM2 þ NU6027; E, 1E1X þ oxindole 16; F, 1JSV þ oxindole 16; G, 1BMK þ 14e; H, 3LCK þ PP2. It has been observed that binding pockets have very poorly and very highly mobile regions.41 Indeed, this could be observed within each PK

Optimization of the simulation temperature for the de novo receptor generation procedure

Table 4. Lowest free energy conformations of the eight protein kinase structures generated by using a de novo generation of multiple receptor conformations

Four complexes with ligands belonging to different chemical spaces (1YDT, 1H1Q, 1FVT and 1DM2) were chosen to run a Monte Carlo energy optimization procedure at different temperatures. ˚ from the ligand and torSide-chains within 6.5 A sion angles of the ligand were randomized at 458 amplitude and the position of the ligand center of ˚ . For each commass was displaced randomly 2 A plex, on average ten simulations were performed at 300 K, 600 K, 1000 K and 2000 K. For a given complex and temperature we recorded the number of low-energy conformations found and if the global minimum was found or not. The low-energy conformations were defined as those within 10 kcal/mol of the global minimum. We found that 600 K was the best choice for Monte Carlo simulations, in agreement with what has been already found for conformational searches in proteins.40 De novo receptor generation of accurate PK structures using IFREDA We selected eight PK complexes for modeling through docking to a flexible receptor either

Complex

Ranking

Ligand RMSD ˚) (A

DDG (kcal/mol)

A

1FMO þ balanol

1 4

1.4 2.7

0.0 7.7

B

1JLU þ staurosporine

1 2

0.8 3.0

0.0 8.2

C

1JLU þ H7

1 3

0.7 4.0

0.0 11.9

D

1DM2 þ NU6027

2 1

1.8 8.1

þ0.5 0.0

E

1E1X þ oxindole 16

1 2

2.4 4.4

0.0 2.2

F

1JSV þ oxindole 16

1 2

2.1 4.3

0.0 1.3

G

1BMK þ 14e

1 2

1.9 7.7

0.0 5.6

H

3LCK þ PP2

1 3

0.9 6.1

0.0 2.2

˚ and kcal/mol, respectRMSD and free energy values are in A ively. DDG is the difference in free energy with respect to the bestenergy conformation.

Protein Flexibility in Docking and Virtual Screening

sub-family. In the eight PK complexes, the parts regarded as flexible during the procedure (loops) include the mobile glycine-rich flap plus other residues that depend on the particular PK subfamily. They were chosen on the basis of the observed flexible parts, comparing crystal structures of the same PK sub-family. These selections are somewhat arbitrary and could be expanded. The following residues were selected: (i) cAPK 49– 58, 70– 74, 120 –127, 170 –173, 181– 187 and 322– 330; (ii) CDK-2 8 –18, 80 – 86, 131– 134, 143 – 147; (iii) P38 30– 38, 104 –112, 50– 54, 167 –169; (iv) LCK 250 –260, 316– 323. The eight complexes were optimized using the de novo structure generation tool of the IFREDA procedure as described in Materials and Methods. In Table 4 we show the RMSD of the best-ranking solution together with the rank, RMSD and relative free energy of the first significantly structurally different conformation. We selected the best conformations solely on the basis of energy calculation, and it is seen ˚ of their that most of the ligands are within 2 A native structure. 1FMO þ balanol complex. In the cAPK sub-family, balanol cannot dock correctly in the adenosinebound pocket (PDB entries 1BKX and 1FMO), since it clashes with the glycine flap (Gly50-ThrGly-Ser-Phe-Gly-Arg-Val). This loop moves about ˚ upwards upon balanol binding, as it is seen in 2A structure 1BX6. In Figure 5 we show how model A reproduces this effect. Balanol superimposes very well with its native structure (1BX6) making a large number of contacts with the generated receptor structure. 1JLU þ staurosporine complex. Staurosporine, a

217

flat and rigid compound that could be docked correctly only to its native structure, induces conformational changes in its neighboring residues, especially in the glycine flap and the C-terminal loop containing Phe327, which is pushed away by ˚ . We show in Figure 6 staurosporine more than 2 A how this induced fit was indeed reproduced in model B, while staurosporine superimposes excellently with its native structure. 1JLU þ H7 complex. Ligand H7 could dock correctly to apo structure 1JLU but scored poorly. The flexible receptor docking of H7 to 1JLU has a very good RMSD (see Table 4). It is shown below that during the small-scale VS, while the scoring of H7 in model C improved (H7 is now placed within the top-ranking 2%), H7 was not placed at the top of the hit list, probably due to the omission of water molecules in the VS. 1DM2 þ NU6027 complex. Superposition of structures 1DM2 and 1E1X shows that ligand NU6027 cannot dock correctly into 1DM2 because the amino group in position 2 clashes with Leu83. This clash pushes the ligand away and the optimal contacts for binding cannot be achieved. In the optimized model D, the backbone of the hinge region including Leu83 rearranges, restoring the triplet hydrogen bond pattern of NU6027 with Glu81 and Leu83. The pyrimidine ring superimposes very well with that of 1E1X, and the ˚ comes mainly from the deviation RMSD of 1.8 A of the high B-factor cyclohexyl moiety. In fact, the crystal structure of 1E1X does not show an optimal conformation for the hydrogen bond between the 2-amino group of NU6027 and the O in Leu83, since the angle O –H –N is , 1148, while in the

Figure 5. Flexible docking of balanol into the 1FMO binding pocket (model A). The modeled structrure is in grey, while 1FMO is displayed in yellow. Balanol carbon atoms are displayed in yellow (corresponding to native structure ˚ 1BX6) and white (modeled complex). Notice the displacement in the modeled structure of the glycine-rich flap , 2 A upwards to allow balanol to bind. (The picture was constructed using ICM version 3.0.)

218

Protein Flexibility in Docking and Virtual Screening

Figure 6. Flexible docking of staurosporine into the 1JLU binding pocket (model B). The modeled structrure is in grey, while 1JLU is displayed in yellow. The glycine-rich flap has been cut for the sake of clarity. Staurosporine carbon atoms are displayed in yellow (corresponding to native structure 1STC) and white (modeled complex). Notice the displacement in the modeled structure of the C-terminal loop, which contains Phe327. (The picture was constructed using ICM version 3.0.)

optimized model this angle is , 1508, closer to the optimal 1808. However, in model D, the distance between NH of Leu83 and N-1 of NU6027 is some˚ ), which may account for the what large (, 2.4 A fact that NU6027 does not rank at the top in the small-scale VS (see below). In model D, the energy of the best pose was slightly worse than that of a misdocked one, by 0.5 kcal/mol (see Table 4), which we consider within the accuracy of our energy function (1 cal ¼ 4.184 J). Of course, the misdocked conformation could be discarded based on structural considerations of the binding mode to PKs. 1E1X þ oxindole 16 complex. Compound oxindole 16 is unable to dock to structure 1E1X, mainly due to the positions of Lys33 and Lys89. Superposition with the native structure of oxindole 16 (1FVT) shows that Lys33 clashes with the Br-benzyl part of oxindole 16, while Lys89 clashes with the SO2 part. Modeling of 1E1X with oxindole 16 by flexible receptor docking (model E) shows a rearrangement of both lysine residues. However, Lys33 was anchored by a salt-bridge with Asp45 and could not adopt the completely buried and folded position that is exhibited in 1FVT, whereby it has no close contact. In this way, oxindole 16 could not fit completely into the binding pocket, with a displacement towards the exposed solvent region ˚ . The hydrogen bond between the indoliof , 1.5 A none O and HN of Leu83 is conserved but that with the carbonyl group of Glu81 is lost. 1JSV þ oxindole 16 complex. The clash between Lys33 and oxindole 16 is the reason why this com-

pound could not be docked to 1JSV. Flexible receptor docking of oxindole 16 into 1JSV (model F) shows the same effects as in model E. In this ˚ , and the average case, the ligand is shifted , 1.2 A RMSD is marginally better (see Table 4). In this case, the hydrogen bond with HN of Leu83 is conserved and that with Glu81 is lost, but a new one is formed between the O of Leu83 and the H at N17 in oxindole 16. 1BMK þ 14e complex. Compound 14e is another example of correct docking to 1BMK and poor scoring. The triplet hydrogen-bonding pattern between the ligand and the backbone (O of Leu107, NH of Met109 and NH of Gly110) could not be reproduced through flexible ligand – rigid receptor docking; only the hydrogen bond with Met109 could be achieved. Using IFREDA to generate the complex structure 1BMK þ compound 14e (model G) restores the hydrogen bond with Leu107 and keeps that with Met109. The explanation of why the hydrogen bond with Gly110 was not restored lies in the fact that in the native crystal structure of compound 14e (1M7Q) a peptide flip between Met109 and Gly110 is observed, which could not be reproduced in our model. With that flip, the HN of Gly110 faces the pocket, thus enabling the formation of the hydrogen bond (this peptide flip effect has been studied recently by Fitzgerald et al.42). Despite this fact, the ligand RMSD of model G is good, especially taking into account that the shift comes mainly from the high B-factor piperazine moiety. 3LCK þ PP2 complex. The main reason why the

Protein Flexibility in Docking and Virtual Screening

Table 5. Improvements in enrichment factors and RMSD values for focused libraries by using the merging and shrinking procedure when virtual screening is performed against the native structure and another one generated using IFREDA

PDB entry

EF (top 1%)

EF (top 2%)

EF (top 10%)

%Binders with ˚ RMSD , 2.5 A

1FMO 1FMO þ balanol Combined EF

16.7 50.0 50.0

16.7 33.3 33.3

5.0 8.3 8.3

50.0 83.3 83.3

1JLU 1JLU þ staurosporine Combined EF

16.7 16.7 33.3

8.3 8.3 16.7

3.3 1.7 5.0

66.7 50.0 83.3

1JLU 1JLU þ H7 Combined EF

16.7 16.7 16.7

8.3 25.0 25.0

3.3 6.7 8.3

66.7 66.7 83.3

1DM2 1DM2 þ NU6027 Combined EF

25.0 16.7 41.7

16.7 12.5 25.0

3.3 5.0 6.7

41.7 50.0 75.0

1BMK 1BMK þ 14e Combined EF

16.7 16.7 16.7

8.3 8.3 16.7

3.3 5.0 5.0

66.7 66.7 66.7

RMSD values refer to the best-energy solution.

grid docking of PP2 to structure 3LCK fails is not due to clashes but to a weak hydrogen bond interaction pattern between PP2 and the receptor, as it is seen when superimposing the native structure of PP2 (1QPE) with 3LCK. In model H, this pattern is improved slightly and in the ligand docking using a rigid receptor approach of PP2 to model H ˚. the ligand is placed within 1.3 A Improvements in RMSD values and enrichment factors using IFREDA The receptor structures thus obtained were used to perform small-scale VS in the same way as was done for the crystal structures (see above). The RMSD and EF values for the corresponding native receptors and those generated through IFREDA are shown in Table 5, together with the values obtained by merging and shrinking of the screening results. The EF increase for five analyzed groups with respect to their single original structure was as high as 1.89 ^ 0.60 (averaged over 15 EF values), which is even slightly higher than the value obtained through using multiple experimental structures (1.85 ^ 0.65), while no decrease in the enrichment factor was observed. Although tested on a limited set, this result shows that the IFREDA is an efficient tool to incorporate receptor conformational diversity in the VS procedure. For the model A, balanol got the best score. However, two other ligands among six (H8 and H89) were placed in the top 1% ranking. It should be remarked that H89 was not even able to dock correctly to 1FMO. All of the ligands except for ˚ . This shows staurosporine are docked within 0.9 A that modeling by flexible receptor docking does

219

not make a tight custom-pocket that accommodates only the ligand used in the optimization, as will be seen in other cases. Combination of EFs of both VS does not make any difference in this case: the generated receptor structure alone is enough to improve the EF significantly. The structure of 1JLU complexed with staurosporine (model B) is expected to behave in a fashion similar to that of native structure 1STC. And this is exactly the case. Staurosporine gets the top ranking score in the database, but no other ligand is present in the top 10% ranking. However, beside staurosporine, adenosine and H7 were ˚ but scored poorly. The comdocked within 0.8 A bined EF is the arithmetic sum of each VS against the native structure and the model and the increase is significant. As stated previously, compound H7 docked correctly in 1JLU, but has a rather poor score (it ranked within the top 30%). The VS on model C shows H7 now in the top 2% ranking. Compound H89, which could not dock to 1JLU, is now able to dock in the generated binding pocket and both H8 and H89 also improve their score significantly compared to 1JLU. Adenosine, H7, H8 and H89 ˚ , and balanol within 2.7 A ˚, are docked within 0.9 A somewhat worse than to 1JLU, but it still ranks in the top 1%. One possible explanation for the fact that H7 does not score better in model C could be found in the fact that water molecules were not considered in the VS. Compound NU6027, which could not dock to ˚ in model 1DM2, was placed within 1.8 A D. Although its score was not optimal and NU6027 ranked in the top 7%, other compounds that were either “misdocked” in 1DM2 (NU2058, NU6094 and NU6102) or docked correctly with poor score (H717) were ranked in model D within the top 1.5%, while compound NU2058 was within the top 5% ranking. As a consequence, the EF of the combined set is higher than the individual EF values for each of the three cut-offs, and the RMSD values are improved. Ligand docking and small-scale VS against models E and F (1E1X þ oxindole 16 and 1JSV þ oxindole 16, respectively) using a rigid receptor model failed to place oxindole 16 in its native pose but, surprisingly the similar compound oxindole 91 was docked in place with a top ranking score. Although the RMSD of oxindole 16 in models E and F is within an acceptable range (see Table 4), we believe that further refinement might be necessary to get better structures to be used in VS. In the model of 1BMK with compound 14e, this ˚, compound ranks first with an RMSD of 2 A slightly worse than in 1BMK but with a much better score. In this case, the improvement in the EF values obtained by merging and shrinking is modest. As expected, there is no improvement in RMSD values, since compound 14e docked already correctly into 1BMK and the side-chain rearrangement was supposed to be small, just to improve ligand – receptor contacts.

220

No EF improvement is reported for LCK structures, since only two different ligands were tested. Compound PP2, which could not dock in apo struc˚ and ture 3LCK, could dock in model H within 1.3 A ranked in the top 5%. Staurosporine could dock cor˚ ) and ranked in the top 5%. rectly (RMSD 0.9 A

Discussion Receptor flexibility: which is the best choice to deal with multiple receptor conformations? Potential binders can be lost (ranked poorly) in the VS for the following reasons. (i) They are misdocked because of clashes with the receptor or because the ligand –protein contacts are not strong enough to hold the ligand in a correct position (for example, if the quality of the crystal structure is bad or residues within the pocket have high B-factors). (ii) They are docked correctly, but they do not score properly because the ligand – receptor contacts are not optimal. (iii) They are misdocked or they do not score properly because water molecules or ions are not included in the receptor model. (iv) Uncertainty in the ionization state of the ligand or the receptor, due to receptor-induced (ligand-induced) pKa changes in the ligand (receptor). (v) They are misdocked because of insufficient sampling or they are docked correctly, but they do not score properly because of failures in the scoring function. The first two reasons are related to rearrangements of the binding pocket upon ligand binding. Direct incorporation of protein flexibility in screening of large compound libraries is computationally too expensive. The use of many receptor structures has been characterized as the best choice to incorporate receptor flexibility in ligand docking.19 However, there are many unanswered questions regarding how these structures should be used, generated and, in particular, how the flexibility consideration impacts the success of the VS procedure. The merging-and-shrinking of screening results from multiple but fixed receptor conformations (either experimental or computergenerated) can deal with both side-chain and backbone flexibility, as is the case of balanol binding to 1FMO (Figure 5) or staurosporine binding to 1JLU (Figure 6). Using the merge-and-shrink procedure with two and three receptor structures (either experimental or generated) we obtained an EF fold-increase of , 1.90, showing that a small ensemble of PK structures can satisfactorily span the conformational space of the binding pocket. The computational effort for flexible ligand– rigid receptor docking is about one to two minutes (700 MHz processor, 1 Mb RAM for the dual-processor node), so screening against many different structures (say five) is computationally still affordable. In cases where a large number of experimental structures are available, a small-scale ligand docking and screening

Protein Flexibility in Docking and Virtual Screening

may be used as a tool to help to eliminate redundant structures from the ensemble, while keeping the most structurally diverse. Comparison with other methods is difficult, since most of them analyze the impact of receptor flexibility on docking accuracy only, and all of them, including that presented here, have been tested on a limited number of protein families, sometimes even in just one protein with many available crystal structures, thus leaving an open question regarding transferability. However, some analogies and differences with other methods can be pointed out. By using many receptor conformations, we avoid directly the use of averaged composite grid maps, in which the structural diversity of the pocket is represented by only one structure.20,22,28 While in those studied cases some of the weighted-average grid maps performed satisfactorily in evaluating RMSD values and binding energies of known binders, there might be other cases in which a single grid map will not be able to represent receptor flexibility, specially if loop rearrangements are involved. Other alternatives, like generating an ensemble of structures through combination of rotameric states of key residues,26 will lead to a combinatorial multiplication of structures, which makes it unfeasible to be used in VS of hundreds of thousand of compounds. The incorporation of flexibility through discrete alternative conformations of varying parts of the protein that are merged in a combinatorial way and considered directly during the ligand docking (FlexE21) is definitely an interesting approach. The method has been evaluated using a rigid ligand docking algorithm and is currently limited to the use of source structures that have a similar backbone trace. The time for protein preparation and ligand docking (per ligand) is roughly equivalent to dock and score one ligand to five structures using IFREDA. Limitations of IFREDA: a guide for further improvements It is important to remark here that the de novo generation of receptor structures through IFREDA does not make a tight custom-pocket that accommodates only the ligand used in the optimization, as shown in Results. Two of the eight generated structures were not suitable for VS (models E and F), since it was not possible to re-dock the ligand used in the optimization using a rigid receptor approach. Those structures were far from useless, since other known binders could be docked and scored correctly. However, failing to re-dock the ligand used in the complex structure generation should be taken as evidence that the model should be improved. This could be linked to the comparatively high RMSD ˚ ) in both models values of the ligand (, 2 A (see Table 4). Two main improvements could be incorporated in our de novo generation procedure. (i)

221

Protein Flexibility in Docking and Virtual Screening

Performing the stochastic energy optimization on a wider selection of side-chain torsional angles (for example, randomly perturbating and minimizing with respect to L2; see Stochastic global energy optimization in Materials and Methods). (ii) Using a better solvation energy term during the stochastic optimization (for example, GB/SA). This could help to improve the energy discrimination of mis-docked conformations (for example, in model D, where two conformations were energetically indistinguishable). Although backbone flexibility is incorporated in the IFREDA, we were not able to reproduce the peptide flip between Met109 and Gly110 of 1M7Q compared to 1BMK. To overcome the large potential barrier of this flip, a Monte Carlo sampling of the backbone atoms would be necessary and might be incorporated in future developments.

Conclusions We present an algorithm, called IFREDA, to incorporate protein flexibility in ligand docking and scoring, specially useful in a drug discovery scenario where few or no holo experimental structures are available, while some binders are known. Initially, IFREDA generates a set of receptor conformations by docking flexible ligands to a flexible receptor through a global energy optimization of the complexes. Both side-chains and essential backbone distortions are included in this optimization, thus sampling the conformational space of the receptor, even in cases of loop rearrangements. The receptor structural set is then used to perform flexible ligand– grid receptor docking and scoring, followed by merging the VS scores and keeping the best rank for each compound. The scoring list is thus shrunk to the size corresponding to a VS against a single receptor. To underscore the impact of protein flexibility on RMSD and EF values, and the necessity of somehow including induced-fit effects in VS, in the first stage of our study we report ligand docking and VS tests against 33 crystal structures of four different sub-families of protein kinases; cAPK, CDK-2, P38 and LCK. We used a 1000 compound library built using a random set of molecules seeded with co-crystallized PK ligands. Native ligands were ˚ to their docked with an average RMSD of 0.74 A co-crystallized structure (maximum RMSD was ˚ , see Table 2), while , 70% of the ligands 1.5 A ˚ , and 80% of the were cross-docked within 2.5 A native ligands scored in the top-ranking 1.5%, which altogether can be taken as a validation of our docking and scoring algorithm, even though we did not make any attempt to improve or tune the scoring function for PKs. The fraction of ˚ and ranking within docked ligands within 2.5 A the top 10% is , 49% (see Figure 3). Comparing this with the , 70% of compounds that are crossdocked correctly, it is observed that induced-fit effects have stronger effect on scoring and ranking than on ligand docking accuracies.

We then showed that merging the VS results against two and three diverse experimental receptor conformations by retaining the best ranking position for each ligand, led to an average EF foldincrease of 1.85 ^ 0.65, with less than 3% cases in which the EF actually deteriorates due to an increase of false positives (see Figure 3). In the merged set, most of the ligands are predicted ˚ RMSD from their native pose (more within 2.5 A than 75% in 19 out of 21 of the groups). Since the computing time for the flexible ligand –rigid receptor docking is about one to two minutes using a 700 MHz processor (1 Mb RAM for the dualprocessor node), screening against a few number of structures is affordable and has the advantage of using actual receptor conformations that may differ in side-chain conformation, and in the rearrangement of loops. The IFREDA procedure has been validated in eight PK complexes and was able to generate the correct receptor and ligand bound conformation and energetically discriminate it from mis-docked conformations within the accuracy of the energy function, even in cases where loop rearrangement ˚ ) was necessary for the ligand binding (see (, 2 A Table 4). Since it has been shown that binding sites have regions of very high and very low stability,41 some portions of the backbone were considered as rigid, while others and the side-chains within the binding pocket were flexible. The choice of flexible backbone parts has been undertaken by inspecting the mobile parts in the available crystal structures. In some cases this is coincident with high B-factor regions. Results from the small-scale VS against the generated set of structures were condensed using the merging and shrinking procedure. We observed an EF fold-increase of 1.89 ^ 0.60 (see Table 5), slightly better than using multiple experimental structures, which clearly shows that multiple receptor conformations generated through IFREDA represent the structural space of the binding pocket and that these conformations can be used when multiple experimental structures are not available, or when the conformational space of the receptor is not represented by the experimental structures available. Interestingly, our de novo structure generation does not make a custom-pocket that is suitable only for the ligand used in pocket generation. Our methodology was successful in incorporating protein flexibility even in cases of loop rearrangements and improving the results of ligand docking and VS. A better solvation energy term for the Monte Carlo sampling and a methodology to allow backbone flipping will be incorporated in future developments.

Materials and Methods Receptor preparation The coordinates of protein kinase complexes were

222 taken from the RCSB Protein Data Bank (PDB).43 Kinase ˚ or incomplete complexes with resolution higher than 2.8 A ligand structures were discarded. Structures containing charged ligands with counter ions were not included. The 33 selected complexes are listed in Table 1. In each subfamily, one apo structure was included. Two adenosinebound structures of cAPK (1BKX and 1FMO) and two staurosporine-bound structures of LCK (1QPD and 1QPJ) were included, since they exhibit different conformations for the glycine-rich flap. P38 apo structure 1P38 was preferred over 1WFC due to a higher resolution. Hydrogen and missing heavy atoms were added to the receptor structure followed by local minimization using the conjugate gradient algorithm and analytical derivatives in the internal coordinates space. The energetically most favorable tautomeric state of His was chosen. Positions of Asn and Gln were optimized to maximize hydrogen bonding. In cases where the occupancy was not equal to 1 we chose the conformation with the lowest energy. The position of polar hydrogen atoms (including those of water molecules) in the ˚ of the ligand ligand-binding pocket (LBP) within 5 A were optimized but water molecules and peptide PKI(5 – 24) were then removed whenever present for subsequent calculations. No structures with alternative hydroxyl rotamers or histidine tautomers were generated for VS, since receptor flexibility was addressed directly in this work. Small-molecule library preparation

Protein Flexibility in Docking and Virtual Screening

library was built by merging the corresponding ligands with the compounds from the random database. Since the number of ligands depends on the PK sub-family, some compounds from the random database were deleted in a random manner to ensure a total of 1000 compounds. Energy calculation and optimization The molecular system is described in terms of internal coordinates variables, using a modified ECEPP/351 forcefield with distance-dependent dielectric constant for the energy calculations as implemented in ICM.14,37 The biased probability Monte Carlo (BPMC) minimization procedure was used for global energy optimization.50,52 The BMPC global energy optimization method consists of the following steps: (1) a random conformation change of the free variables according to a predefined continuous probability distribution;50,52 (2) local energy minimization of analytical differentiable terms; (3) calculation of the complete energy including nondifferentiable terms such us entropy and solvation energy; (4) acceptance or rejection of the total energy based on the Metropolis criterion53 and return to step 1. Flexible ligand/grid receptor docking The VS module as implemented in ICM35 – 38 was used. The energy function used during the flexible ligand –grid receptor docking simulations consisted of:36 E ¼ EFFint þ Evw þ Eel þ Ehb þ Ehp

The 3D coordinates of the ligands were taken from the crystal structure and their correct stereochemistry and formal charges were assigned. The ligands are listed in Table 1 together with their PDB three-letter code. The protonation state was determined according to an environment at pH 7.4. The following groups were regarded as charged: staurosporine amino group,44 nitrogen atoms of terminal amino substituents in isoquinolinesulfonamide derivative inhibitors H7, H8 and H89,45 the guanidine ring of hymenialdisine,46 the adenine-derived compound H717 amino group at its diaminocyclohexane substituent at C2,47 nitrogen in the cyclohexyl amine group of compound SB220025, piperidine nitrogen in dyhydroquinazoline derivative (PDB code dqo).48 Each ligand was then assigned the MMFF atom types and charges,49 and subjected to a global energy optimization using the ICM stochastic optimization algorithm.50 A random database of compounds was generated using every 272nd molecule of the Diverse Set database of ChemBridge (ChemBridge, Inc., San Diego, CA) consisting of 272,938 compounds. Carboxylic acids and primary, secondary and tertiary amine groups were regarded as being charged. These compounds were prepared and energy optimized in a fashion similar to that used for the ligand compounds. The overlap of basic chemical descriptors between the ligand and the random libraries is essential for meaningful results using the VS protocol. Four properties were calculated: molecular mass ðMÞ, number of rotatable bonds ðNrot Þ, number of hydrogen bond acceptors ðNHBacc Þ and donors ðNHBdon Þ: The average and standard deviation values for the ligand and random libraries were respectively as follows: M, 351.5 ^ 83.0 and 364.0 ^ 73.8; Nrot , 2.3 ^ 1.9 and 2.4 ^ 2.0; NHBacc , 3.9 ^ 1.9 and 3.6 ^ 1.5; NHBdon , 3.2 ^ 1.6 and 1.1 ^ 0.9. For each PK sub-family, a 1000 compound sample

where EFFint is the internal force-field energy of the ligand, and Evw ,Eel ,Ehb and Ehp are the van der Waals, electrostatic, hydrogen bond and hydrophobic potential terms, respectively. The last four terms are pre-calculated ˚ to accelerate energy on a grid spacing of 0.5 A evaluation. The van der Waals potential in its standard 12-6 form is too sensitive and may introduce noise into the energy function. For intermolecular interactions, a softer truncated van der Waals potential is used according to: 8 o if Eovw # 0 > < Evw Evw ¼ Eovw Emax > if Eovw . 0 : o Evw þ Emax where the Emax value was chosen as 1.5 kcal/mol. The total van der Waals map is composed by three maps representing hydrogen atoms, first-row atoms (C, N, O and F) and heavy atoms (P, S, Cl, Br and I). The electrostatic map is calculated using a distancedependent dielectric constant 1 ¼ 4r: The hydrogen bonding potential is represented by spherical Gaussians centered at the donor or acceptor sites, according to: Ehb ðrÞ ¼ Eohb e2kr2rep k

2

=d2hb

where the points r denote the locations of the interaction ˚ from the atom center. In the centers, which are at 1.7 A case of hydrogen atoms, the center is placed along the axis of the covalent bond attaching the hydrogen atom to the rest of the molecule. In the case of heavy sp2 atoms, one (for nitrogen) or two (for oxygen) centers are placed at an angle of 1208 to the existing covalent bond. For sp3 oxygen and sulfur atoms, two centers are placed in tetrahedral geometry at 1098 to the existing covalent bonds and to each other. The radius of the interaction

223

Protein Flexibility in Docking and Virtual Screening

˚ . The Eo for donor and sphere is set to dhb ¼ 1:4 A hb acceptor atoms is assumed to be 2.5 kcal/mol. The hydrophobic potential on point r on the grid is calculated as: 2

2

Ehp ðrÞ ¼ Eohp e2dsurf ðrÞ=dw where dsurf ðrÞ is the distance from r to the closest point on the hydrophobic surface, and dw is the effective radius of the hydrophobic interaction, which is taken as the ˚ . The value of Eo diameter of the water molecule, 2.8 A hp 3 kcal/mol was chosen to approximate the surface 2 ˚ tension of 30 cal/(mol A ) for extended hydrophobic surfaces in test cases. The energy of the free ligand is minimized in the absence of the receptor prior to docking. The ligand is then docked into the grid representation of the receptor using the global energy optimization method described above. The best energy conformation was scored using an empirical scoring function based on its fit into the binding pocket. In this way, the reported RMSD values of ligands always refer to the best-energy conformation. The ICM scoring function consists of the following terms:35 – 38 Escore ¼ DEIntFF þ TDSTor þ a1 DEHBond þ a2 DEHBDesol þ a3 DESolEl þ a4 DHPhob þ a5 QSize where DEIntFF includes the van der Waals interaction of the ligand with the receptor as well as the internal force-field energy of the ligand; TDSTor is the ligand conformational entropy loss contribution upon binding, which is assumed to be proportional to the number of free torsions; DEHBond is the hydrogen bonding term; DEHBDesol is the term that accounts for the disruption of hydrogen bonds with solvent upon ligand binding (desolvation of hydrogen bond donors and acceptors). DESolEl is the solvation electrostatic energy change upon binding, calculated solving the Poisson equation using the boundary element method as implemented in the REBEL module of ICM.54 DEHPhob is the hydrophobic free energy gain assumed to be proportional to the accessible surface area of the hydrophobic atoms of the receptor and ligand, buried upon binding. The surface ˚ 2. QSize is a tension parameter was set at 0.012 kcal/A size correction term to avoid bias towards larger ligands. While this term has no direct physical meaning, it may account for otherwise omitted interactions of the ligand with the solvent, primarily the van der Waals dispersion interaction. The weights a1 – a5 give each interaction the appropriate strength and were optimized on a diverse benchmark of complexes.35 Since many complexes present non-satisfied hydrogen bonds for the oxygen atoms of the sulfonamide moiety, no penalty for hydrogen bond desolvation was considered in these type of oxygen atoms. The docking and scoring process was repeated four times, and the best score for each ligand was retained. De novo receptor structure generation using IFREDA Seeding The ligand was placed in the LBP by overlaying its co-crystallized receptor structure with the one to be modeled. The ligand was then flipped 1808 with respect to its principal axes of inertia, thus generating four starting conformations. This was undertaken to avoid any bias towards the preferred ligand conformation, though in

real-life calculations the knowledge of preferred ligand– receptor interactions in PK might help. Conformations of the ligand transversal to the PK pockets were avoided based on physical evidence of the binding mode to PKs. Starting with each of these four poses, the position of the center of mass of the ligand was displaced randomly ˚ and the orientation of the molecule within a sphere of 3 A in space was randomized around a random axis with an ˚ /the radius of gyration of the molecule. amplitude of 3 A This procedure was repeated ten times for each starting pose, so that 40 conformations were obtained. The ampli˚ for random displacement was chosen for two tude of 3 A reasons: to avoid generating useless conformations with the ligand outside the pocket and because the main interest is to generate alternative interaction modes of the ligand with the neighboring residues. It should be made clear that the objective of this seeding and the subsequent minimization is merely to obtain an ensemble of different starting receptor þ ligand conformations. The full global optimization that predicts the pocket conformation and the ligand pose is described below. Soft-van der Waals structure relaxation The definition of the binding pocket is flexible and, of course, can include the whole protein. Based on observations that binding sites have regions of very high and very low stability,41 here, we took the following approach: Some loops (which definition depends on the specific PK) were regarded as flexible, which means that both backbone and side-chains torsion angles were free (referred to here as Loops). The rest of the backbone atoms of the receptor were kept fixed. We then defined three layers of atoms (L1, L2 and L3) of increasing size as follows: L1, side-chains that have at least one atom in ˚ ; L2, contact with any atom of the ligand within 4.5 A L1 þ side-chain atoms of Loops þ residue side-chains that have at least one heavy-atom in contact (within ˚ ) with L1 or with Loops’ side-chains; L3, L2 þ 4.5 A backbone atoms of Loops þ residue side-chains that ˚ ) with have at least one heavy-atom contact (within 4.5 A Loops’ backbone. Regarding L3 as flexible by unfreezing its torsion angles, each of the 40 complexes was subjected to five cycles of in vacuo minimization in the internal coordinates space using a soft-van der Waals potential,36 with increasing weight in each cycle. In this way, we avoid high-energy gradients at the beginning of the minimization, and thus the binding site does not fall apart. Bi-quadratic distance restraints14 between the ligand heavy atoms and neighboring Ca are imposed so that a penalty is added to the energy when the ligand displaces ˚ from the initial position. Again, the more than 4.5 A objective of these restraints is to keep the ligand within the physical limits of the pocket, avoiding the generation of unrealistic complex conformations for the next stages. Stochastic global energy optimization The 40 complexes generated through seeding and in vacuo minimization were subjected to global energy optimization using the BPMC procedure, during which torsional variables associated with L1 atoms were perturbated randomly and local minimization of differentiable energy terms was performed with respect to L2. The number of energy evaluations during the BPMC and local minimizations was limited to 2,000,000 and 2000, respectively. Entropy50 and solvation energy terms were

224

then calculated and added to the in vacuo energy to be used in the acceptance/rejection stage. The solvation energy term was based on atomic solvation parameters.55 Torsional angles were sampled between 2 1808 and 1808, while the amplitude for ligand random displacement ˚ and 3 A ˚ /ligand radius of and rigid rotation was 3 A gyration, respectively. The optimal temperature for the simulation was found to be 600 K, as described in Results. During each of the 40 simulations, 40 geometrically diverse low-energy states are stored in a conformational set.40 These 40 conformational sets (each with 40 low-energy conformations) are then merged and their electrostatic and solvation energy contributions reevaluated. The electrostatic part was calculated using MMFF charges,49 and the interaction with the solvent by solving the Poisson equation using the boundary element algorithm.39 The non-polar contribution to the solvation energy was assumed to be proportional to the solvent-accessible surface. Energy values for isolated ligand and receptor were not calculated, since we were interested only in relative values of the free energy. The conformational set is then clustered to eliminate redundant conformations by comparing the RMSD of the atomic coordinates of the ligand heavy atoms. The ˚ . The top ranking conformations within cutoff was 0.4 A 30 kcal/mol (with a minimum of 15 and maximum 40) were subjected to full minimization with regard to L3 and then the energy re-evaluated as described above. The best-energy conformation was used for subsequent VS calculations.

Protein Flexibility in Docking and Virtual Screening

10. 11. 12. 13.

14.

15.

16.

17.

Acknowledgements 18.

We thank Andrew Orry & Maxim Totrov for many useful discussions. 19.

References 1. Amzel, L. M. (1998). Structure-based drug design. Curr. Opin. Biotechnol. 9, 366– 369. 2. Rosenfeld, R., Vajda, S. & DeLisi, C. (1995). Flexible docking and design. Annu. Rev. Biophys. Biomol. Struct. 24, 677– 700. 3. Taylor, R. D., Jewsbury, P. J. & Essex, J. W. (2002). A review of protein-small molecule docking methods. J. Comput. Aided Mol. Des. 16, 151–166. 4. Abagyan, R. & Totrov, M. (2001). High-throughput docking for lead generation. Curr. Opin. Chem. Biol. 5, 375– 382. 5. Shoichet, B. K., McGovern, S. L., Wei, B. & Irwin, J. J. (2002). Lead discovery using molecular docking. Curr. Opin. Chem. Biol. 6, 439– 446. 6. Murray, C. W., Baxter, C. A. & Frenkel, A. D. (1999). The sensitivity of the results of molecular docking to induced fit effects: application to thrombin, thermolysin and neuraminidase. J. Comput. Aided Mol. Des. 13, 547– 562. 7. Davis, A. M. & Teague, S. J. (1999). Hydrogen bonding, hydrophobic interactions, and failure of the rigid receptor hypothesis. Angew. Chem. Int. Ed. Engl. 38, 736– 749. 8. Teague, S. J. (2003). Implications of protein flexibility for drug discovery. Nature Rev. Drug Discov. 2, 527– 541. 9. Bouzida, D., Rejto, P. A., Arthurs, S., Colson, A. B.,

20. 21.

22.

23.

24. 25.

26.

Freer, S. T., Gehlhaar, D. K. et al. (1999). Computer simulations of ligand-protein binding with ensembles of protein conformations: a Monte Carlo study of HIV-1 protease binding energy landscapes. Int. J. Quantum Chem. 72, 73 – 84. Cheney, D. & Mueller, L. (2003). Evaluation of strategies for molecular docking. Abstr. Pap. Am. Chem. Soc. 226, 144. Jiang, F. & Kim, S. H. (1991). Soft docking: matching of molecular surface cubes. J. Mol. Biol. 219, 79 – 102. Leach, A. R. (1994). Ligand docking to proteins with discrete side-chain flexibility. J. Mol. Biol. 235, 345– 356. Jones, G., Willett, P. & Glen, R. C. (1995). Molecular recognition of receptor sites using a genetic algorithm with a description of desolvation. J. Mol. Biol. 245, 43 – 53. Abagyan, R., Totrov, M. & Kuznetsov, D. (1994). ICM—a new method for protein modeling and design—applications to docking and structure prediction from the distorted native conformation. J. Comput. Chem. 15, 488– 506. Desmet, J., Wilson, I. A., Joniau, M., De Maeyer, M. & Lasters, I. (1997). Computation of the binding of fully flexible peptides to proteins with flexible side-chains. FASEB J. 11, 164– 172. Schaffer, L. & Verkhivker, G. M. (1998). Predicting structural effects in HIV-1 protease mutant complexes with flexible ligand docking and protein side-chain optimization. Proteins: Struct. Funct. Genet. 33, 295– 310. Sandak, B., Nussinov, R. & Wolfson, H. J. (1998). A method for biomolecular structural recognition and docking allowing conformational flexibility. J. Comput. Biol. 5, 631– 654. Sandak, B., Wolfson, H. J. & Nussinov, R. (1998). Flexible docking allowing induced fit in proteins: insights from an open to closed conformational isomers. Proteins: Struct. Funct. Genet. 32, 159– 174. Carlson, H. A. (2002). Protein flexibility is an important component of structure-based drug discovery. Curr. Pharm. Des. 8, 1571– 1578. Knegtel, R. M., Kuntz, I. D. & Oshiro, C. M. (1997). Molecular docking to ensembles of protein structures. J. Mol. Biol. 266, 424–440. Claussen, H., Buning, C., Rarey, M. & Lengauer, T. (2001). FlexE: efficient molecular docking considering protein structure variations. J. Mol. Biol. 308, 377– 395. ¨ sterberg, F., Morris, G. M., Sanner, M. F., Olson, A. J. O & Goodsell, D. S. (2002). Automated docking to multiple target structures: incorporation of protein mobility and structural water heterogeneity in AutoDock. Proteins: Struct. Funct. Genet. 46, 34 –40. Morris, G. M., Goodsell, D. S., Huey, R. & Olson, A. J. (1996). Distributed automated docking of flexible ligands to proteins: parallel applications of AutoDock 2.4. J. Comput. Aided Mol. Des. 10, 293– 304. Goodsell, D. S., Morris, G. M. & Olson, A. J. (1996). Automated docking of flexible ligands: applications of AutoDock. J. Mol. Recogn. 9, 1– 5. Morris, G. M., Goodsell, D. S., Halliday, R. S., Huey, R., Hart, W. E., Belew, R. K. & Olson, A. J. (1998). Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. J. Comput. Chem. 19, 1639– 1662. Frimurer, T. M., Peters, G. H., Iversen, L. F., Andersen, H. S., Moller, N. P. & Olsen, O. H. (2003). Ligand-induced conformational changes: improved

Protein Flexibility in Docking and Virtual Screening

27. 28.

29.

30.

31.

32.

33.

34. 35.

36.

37. 38. 39. 40. 41.

42.

43.

predictions of ligand binding conformations and affinities. Biophys. J. 84, 2273– 2281. Schnecke, V. & Kuhn, L. A. (2000). Virtual screening with solvation and ligand-induced complementarity. Persp. Drug Discov. Des., 20, 171– 190. Broughton, H. B. (2000). A method for including protein flexibility in protein-ligand docking: improving tools for database mining and virtual screening. J. Mol. Graph. Model. 18, 247– 257. Carlson, H. A., Masukawa, K. M., Rubins, K., Bushman, F. D., Jorgensen, W. L., Lins, R. D. et al. (2000). Developing a dynamic pharmacophore model for HIV-1 integrase. J. Med. Chem. 43, 2100–2114. Lin, J. H., Perryman, A. L., Schames, J. R. & McCammon, J. A. (2002). Computational drug design accommodating receptor flexibility: the relaxed complex scheme. J. Am. Chem. Soc. 124, 5632–5633. Lin, J. H., Perryman, A. L., Schames, J. R. & McCammon, J. A. (2003). The relaxed complex method: accommodating receptor flexibility for drug design with an improved scoring scheme. Biopolymers, 68, 47 – 62. Cavasotto, C. N., Orry, A. J. W. & Abagyan, R. (2003). Structure-based identification of binding sites, native ligands and potential inhibitors for G-protein coupled receptors. Proteins: Struct. Funct. Genet. 51, 423–433. Taylor, R. D., Jewsbury, P. J. & Essex, J. W. (2003). FDS: flexible ligand and receptor docking with a continuum solvent model and soft-core energy function. J. Comput. Chem. 24, 1637– 1656. Dancey, J. & Sausville, E. A. (2003). Issues and progress with protein kinase inhibitors for cancer treatment. Nature Rev. Drug Discov. 2, 296– 313. Totrov, M. & Abagyan, R. (1999). Derivation of sensitive discrimination potential for virtual screening. RECOMB ‘99. Proceedings of the Third Annual International Conference on Computational Molecular Biology, Lyon, France, pp. 37 –38, ACM Press, New York. Totrov, M. & Abagyan, R. (2001). Protein– ligand docking as an energy optimization problem. In Drug-Receptor Thermodynamics: Introduction and Experimental Applications (Raffa, R. B., ed.), pp. 603–624, Wiley, New York. Molsoft LLC (2003). ICM Manual 3.0, Molsoft LLC, La Jolla, CA. Schapira, M., Abagyan, R. & Totrov, M. (2003). Nuclear hormone receptor targeted virtual screening. J. Med. Chem. 46, 3045– 3059. Totrov, M. & Abagyan, R. (1996). The contourbuildup algorithm to calculate the analytical molecular surface. J. Struct. Biol. 116, 138– 143. Abagyan, R. & Argos, P. (1992). Optimal protocol and trajectory visualization for conformational searches of peptides and proteins. J. Mol. Biol. 225, 519–532. Luque, I. & Freire, E. (2000). Structural stability of binding sites: consequences for binding affinity and allosteric effects. Proteins: Struct. Funct. Genet., Suppl. 4, 63 – 71. Fitzgerald, C. E., Patel, S. B., Becker, J. W., Cameron, P. M., Zaller, D., Pikounis, V. B. et al. (2003). Structural basis for p38a MAP kinase quinazolinone and pyridol-pyrimidine inhibitor specificity. Nature Struct. Biol. 10, 764– 769. Bernstein, F. C., Koetzle, T. F., Williams, G. J., Meyer, E. F., Jr, Brice, M. D., Rodgers, J. R. et al. (1977). The Protein Data Bank: a computer-based archival file for macromolecular structures. J. Mol. Biol. 112, 535–542.

225

44. Toledo, L. M. & Lydon, N. B. (1997). Structures of staurosporine bound to CDK2 and cAPK—new tools for structure-based design of protein kinase inhibitors. Structure, 5, 1551– 1556. 45. Engh, R. A., Girod, A., Kinzel, V., Huber, R. & Bossemeyer, D. (1996). Crystal structures of catalytic subunit of cAMP-dependent protein kinase in complex with isoquinolinesulfonyl protein kinase inhibitors H7, H8, and H89. Structural implications for selectivity. J. Biol. Chem. 271, 26157– 26164. 46. Meijer, L., Thunnissen, A. M., White, A. W., Garnier, M., Nikolic, M., Tsai, L. H. et al. (2000). Inhibition of cyclin-dependent kinases, GSK-3beta and CK1 by hymenialdisine, a marine sponge constituent. Chem. Biol. 7, 51 – 63. 47. Dreyer, M. K., Borcherding, D. R., Dumont, J. A., Peet, N. P., Tsay, J. T., Wright, P. S. et al. (2001). Crystal structure of human cyclin-dependent kinase 2 in complex with the adenine-derived inhibitor H717. J. Med. Chem. 44, 524– 530. 48. Stelmach, J. E., Liu, L., Patel, S. B., Pivnichny, J. V., Scapin, G., Singh, S. et al. (2003). Design and synthesis of potent, orally bioavailable dihydroquinazolinone inhibitors of p38 MAP kinase. Bioorg. Med. Chem. Letters, 13, 277– 280. 49. Halgren, T. (1995). Merck molecular force field I– V. J. Comput. Chem. 17, 490– 641. 50. Abagyan, R. & Totrov, M. (1994). Biased probability Monte-Carlo conformational searches and electrostatic calculations for peptides and proteins. J. Mol. Biol. 235, 983–1002. 51. Nemethy, G., Gibson, K. D., Palmer, K. A., Yoon, C. N., Paterlini, M. G., Zagari, A. et al. (1992). Energy parameters in polypeptides. 10. Improved geometrical parameters and nonbonded interactions for use in the ECEPP/3 algorithm, with application to proline-containing peptides. J. Phys. Chem. 96, 6472–6484. 52. Abagyan, R. A. & Totrov, M. (1999). Ab initio folding of peptides by the optimal-bias Monte Carlo minimization procedure. J. Comput. Phys. 151, 402– 421. 53. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. (1953). Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087– 1092. 54. Totrov, M. & Abagyan, R. (2001). Rapid boundary element solvation electrostatics calculations in folding simulations: successful folding of a 23-residue peptide. Biopolymers, 60, 124– 133. 55. Wesson, L. & Eisenberg, D. (1992). Atomic solvation parameters applied to molecular dynamics of proteins in solution. Protein Sci. 1, 227– 235.

Edited by J. Thornton (Received 9 October 2003; received in revised form 26 December 2003; accepted 6 January 2004)

Supplementary Material comprising two Figures is available on Science Direct