WaterMap: Mapping Solvation Free Energies of Active Site Waters and Using This Data to Understand Binding Affinities

1. WaterMap: Mapping Solvation Free Energies of Active Site Waters and Using This Data to Understand Binding Affinities Richard A. Friesner Columbia University

2. Thermodynamic Decomposition of Ligand/Protein Binding ∆Hconf penalty ∆Sconfpenalty ∆Hconf penalty ∆Sconf penalty ∆G(2) ∆G(1) ∆H penalty ∆S reward ∆H ? ∆S ? ∆G(4) ∆G(3) ∆H reward ∆Srot/trans penalty ∆G(5) Solvated Ligand Solvated Apo Protein Solvated Ligand in Bioactive Conformation Solvated Protein in Ligand-Induced Conformation Desolvated Ligand Ligand-Induced Desolvated Protein Binding Site Protein/Ligand Complex For lead optimization and comparison between two ligands, ∆∆G(4) is typically largest component of ∆∆Gbind

3. Water displacement is key step that drives binding affinity • Hydrogen bonds are very important for specificity, but can be made by the ligand either in the protein or in solution • Some gain due to releasing waters into bulk, mostly entropic • A few exceptional cases based on unusual electrostatic configurations in the active site (neuramididase) • Most waters in the active site, particularly those in hydrophobic regions, do not have competitive entropy and/or hydrogen bonding as compared to bulk water • Hence if ligand group complementary to protein can displace these waters, significant free energy gains accrue • There are “special” regions of the active site where waters have particularly poor free energies • Hydrophobic enclosure • Hydrophobic enclosure plus correlated hydrogen bonds

4. p38 MAP Kinase Naphthyl ring is hydrophobically enclosed in allosteric pocket – huge drop in binding affinity if this group is removed 1kv2

5. HIV-1 Protease No enclosure, hydrophobic groups are only on one “side” of phenyl group of ligand – some benefit, captured in normal hydrophobic pair term, but not nearly as large as when enclosure is present – lack of enclosure explains why HIV protease ligands must be large to acquire substantial binding affinity 1hpx

6. CKD2 Staurosporine in 1aq1: special H bond NH CO pair 3 kcal/mol reward and central part of ring gets a 3 kcal/mol hydrophobic packing reward 1aq1

7. Streptavidin / biotin Classic problem in molecular recognition Enclosure + triple correlated H bond explains huge binding affinity of small ligand – result is within 1.5 kcal/mol of experiment Model has now been validated by all atom MD simulations 1stp

8. Mapping Solvent Chemical Potential – Overview • Based on inhomogenious solvation theory1 • <E> is taken directly from simulation • Se from a local expansion of orientational and spatial correlation functions • Combines MD simulation and trajectory analysis • ~10 ns simulation with explicit waters and restrained protein • Waters are spatially clustered and analyzed • Chemical potential (entropy and enthalpy) are computed for each hydration site • Energy terms are relative to bulk water • Preliminary applications have been published2,3 • Patent filed Sept 2007 r = position  = orientation  = integral over   = bulk water density g = correlation function sw = solute-solvent terms sww = solute-solvent-solvent terms 1Lazaridis T (1998) J Phys Chem B 102:3531-3541 2Young T et al. (2007) PNAS 104:808-813 3Abel R et al. (2008)J Am Chem Soc 130:2817-2831

9. WaterMap • A systematic approach using molecular dynamics (MD) methods to map out free energy of water clusters in the active site • Clustering techniques used to analyze MD trajectory and identify regions of high and low density • High density regions are assigned enthalpies and entropies using inhomogeous solvation theory • Protocol is simple and fast; remove the ligand, run the MD simulation (overnight on 8 processors), and WaterMap is created for that receptor conformation – can then be used to identify key areas for growing ligand, and in favorable circumstances to estimate change in binding free energy if a group displacing a new water is added • Can also run with ligand present – some technical issues which have been overcome in latest version (will show a few examples) • Numerous applications will be discussed in presentations that follow

10. Overview of WaterMap positioning as a theoretical method and as used in practice to date • Accurate “relative” structures are absolutely required and induced fit effects can be critical • Even small changes in structure can alter water network in nontrivial ways • For late stage lead optimization, building structures (vs. free docking) is often best as it promotes cancellation of error • Induced fit calculations are critical for novel structures, and then relative protein reorganization energy must be estimated • Much faster than FEP, can be applied to many ligands after one calculation, for favorable cases more accurate due to noise reduction (more discussion later) • Provides vivid, reasonably accurate pictures, useful for medicinal chemists in design process • Exciting successes already observed in explaining literature data and in drug discovery projects

11. Test of IHT/watermap on simple test cases; methane in hydrophobic model systems, comparison with FEP

13. Example solvent density distribution and clustering (Streptavidin)

14. Factor Xa – initial published watermap results for pharmaceutically interesting target (JACS 2008) • Use inhomogenous solvation theory to calculate binding free energy differences between pairs of factor Xa ligands • Pairs chosen so modification of ligand typically displaces waters near hydrophobic surface • Avoids changes in charge state, hydrogen bonding, etc. upon addition of new group (otherwise other effects must be added, work along these lines is in progress) • Map of water molecules built from single structure, ligands superimposed on this structure – so only one MD simulation is required (works in part because fxa is an unusually rigid protein) • Results explain key features of reverse binders currently in Phase III clinical trials (J&J+Bayer, BMS, etc.)

15. Hydration Site Thermodynamics

16. WaterMap Binding Energy Predictions – Factor Xa Inhibitors(no adjustable parameters) 3 2 1 0 -1 -2 -3 -4 -5 -6 -7 R2 = 0.63 Calculated ∆∆Gbind (kcal/mol) -7 -6 -5 -4 -3 -2 -1 0 1 Experimental ∆∆Gbind (kcal/mol) Abel et al. J. Am. Chem. Soc.2008, 130, 2817

17. FXa – SAR of neutral P1 substituents = stable waters & favorable enthalpy = most unstable water in binding site 10,000 nM 2500 nM 2100 nM 470 nM 89 nM 73 nM 3 nM

18. Case Study 1: thrombin R2 = 0.82 P.I. = 0.91 R2 = 0.26 P.I. = 0.64 R2 = 0.51 P.I. = 0.66 WM/MM ΔG bind (kcal/mol)‏ MM-GB/SA ΔG bind (kcal/mol)‏ WM ΔG bind (kcal/mol)‏ Experimental ΔG bind (297 K kcal/mol)‏ Experimental ΔG bind (297 K kcal/mol)‏ Experimental ΔG bind (297 K kcal/mol)‏ MM-GB/SA WM/MM WM R2MW=0.45

19. The ERBIN PDZ-domain Phage Display reveals that WETWV is the optimal Erbin binding motif (Skelton et al., JBC, 2003) Alanine scanning of the peptide suggest that P-0 and P-1 are most important for affinity Mutagenesis of the PDZ domain cannot explain the crucial role of Trp at P-1: “curiously, alanine substitutions of residues that contact Trp -1 did not decrease peptide binding” (Skelton et al., JBC, 2003)

20. WaterMap calculations explain origin of high affinity and selectivity of peptides Experimental ∆G (kcal/mol) y = 0.5647x + 19.144 R2 = 0.7314 WaterMap ∆G (kcal/mol) Unfavorable water molecules are found in the P-0 and P-1 pockets, consistent with experimental data Peptide affinities correlate very well with WaterMap energies (R2=0.73) Only water molecules with -TS and H larger than 1 kcal/mol (i.e., enthalpically and entropically unstable waters) are shown

21. Abl/c-Kit selectivity Gleevec does not hit this specific high energy water in Abl (red) or c-Kit (green) A methyl variant of Gleevec binds tighter to c-Kit and only displaces the c-Kit water

22. Bcl-xL WaterMap ABT-737 All highly unstable (unhappy) waters in binding site (∆G > 2.8 kcal/mol)

23. High potency ligands displace all highly unstable waters and no stable waters 9 waters in each bin

24. S-enantiomers bind much more weakly than R-enatiomers

25. 2.0 1.5 1.9 R R-enantiomer displaces 3 unstable waters

26. S-enantiomer does not displace 3 unstable waters S WaterMap predicts that R-enantiomer will bind with -pKi 3.5 greater than S-enantiomer Experiment: -pKi difference in binding is 2.5

27. CDK2 inhibitors 1ke6 1ke7 1ke5 1ke9 1ke8 ∆Gbind (kcal/mol) -11.2 -11.0 -8.5 -8.4 -8.2 ∆∆G between tight and weak binders is ~3 kcal/mol

28. Differences in receptor structures result in different location and energy of key water site ∆G = 2.1 kcal/mol Shifted 0.5 Å away from R-group 2iw9 – has cyclin A2 bound ∆G = 3.4 kcal/mol 1ke6 – does not have cyclin A2 bound

29. Using the “right” structure when predicting binding affinities Correlation between experimental ∆G of binding and WaterMap ∆G

30. Using WaterMap with Ligand Present • Should be more accurate since this takes into account how ligand modifies properties of water molecules surrounding it • However, technical problems arise due to creation of pockets of solvent not in contact with bulk – specialized algorithms are needed to populate and sample these regions (particle insertion methods based on Grand Canonical Monte Carlo) • New release has these algorithms and tests indicate substantial improvement; some examples follow

31. 2bhh: largest outlier of XP for CDK2 data set, no obvious explanation • Actual binding affinity ~ 6kcal/mole; calculations yield ~9 kcal/mole • Approach: run watermap with ligand present (currently subject of extensive studies,appears necessary in a number of interesting cases) • Analogue of 2bhh ligand is available with binding affinity data from literature • By comparing the results for the analogue and for the compound itself, a useful hypothesis is suggested: a modified group of 2bhh disrupts the water structure around the ligand; this effect can be quantified relative to the comparison compound, and the differential (which is NOT represented in the new XP/WM scoring function) is in reasonable agreement with the error in the XP/WM prediction for 2bhh

32. CDK2: DDG=~-2.5 kcal/mol (Exp.)‏ 2BHH Carbonyl variant

33. New version of Glide XP – not yet released • First goal is to rank order highly diverse compounds with different charge states, etc. – these data sets do not work well with methods like MM-PBSA, etc. – and at the same time achieve much higher enrichments than previously • Second goal is to eliminate high scoring decoys; current methods with “enrichment” of ~2-5x have huge numbers of incorrectly scored random compounds • Reasonable rank ordering is a necessary condition for making a large improvement • Systematically investigate errors and try to fix them • First develop scoring function that works with self-docking of PDB crystal structures, correcting for core reorganization energy • Then use induced fit methods to enable similar approach to work in cross docking • Key developments • New terms added to XP • Combination of XP and watermap; use of watermap to analyze remaining outliers • Analysis of reorganization energy

34. Integration of WaterMap with Glide XP • WaterMap run on a single structure • Waters inserted into all members of test set; process contains some noise due to induced fit effects; current algorithm is quite crude but works acceptably at targeted level of accuracy • Need to avoid double counting • Watermap results treated as a correction to XP hydrophobic packing score • Objective is to pick up large differences which indicate a qualitative error in empirical model –key improvements seen in rejection of false positives • Use WaterMap to identify rotatable bonds which are in bulk-like regions – here the implicit approximation in the pair hydrophobic term is not relevant and entropic penalty is applied

35. Results of integrating watermap with Glide – tests on PDB data sets (decoy data from docking of 1000 decoys, Schrodinger decoy set)

36. PDE4 decoys • Two fundamental problems identified; neither has anything to do with the scoring function • Problem 1: coordination of metal sites in PDE4 with water molecules • Modeling of metals in problematic in docking calculations • QM studies suggest Zn in PDE4 should have octahedral coordination, or energy is problematic • Water molecules supply this coordination • About 50% of decoys displace one or more key waters with inappropriate group; this never happens in actives • Problem 2: regulatory region blocks part of active site; only appears in recent PDB crystal structures • Affects many of remaining 50% of decoys • If new structure is used for docking, should eliminate this problem • Combination of two effects may explain virtually entire difficulty with PDE4 decoys (NOT a “failure” of scoring function in the traditional sense)

37. “Bumped” chelated solvent

38. “Bumped” chelated solvent

39. “Bumped” C-terminal regulatory region

40. “Bumped” C-terminal regulatory region

41. Example of a good native

42. Coworkers • Columbia • Bruce Berne, Tom Young, Robert Abel, Lingle Wang • Schrodinger • Ramy Farid, John Shelley, Woody Sherman, Thijs Beuman, Rob Murphy