Predicting the Binding of Fatty Acid Amide Hydrolase Inhibitors by

Oct 5, 2018 - Janssen Research and Development , 3210 Merryfield Row, San Diego , California 92121 , United States. J. Chem. Theory Comput. , Article ...
0 downloads 0 Views 1MB Size
Subscriber access provided by University of Sunderland

Molecular Mechanics

Predicting the Binding of Fatty Acid Amide Hydrolase Inhibitors by Free Energy Perturbation (FEP) Arjun Saha, Amy Y. Shih, Taraneh Mirzadegan, and Mark Seierstad J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00672 • Publication Date (Web): 05 Oct 2018 Downloaded from http://pubs.acs.org on October 6, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Predicting the Binding of Fatty Acid Amide Hydrolase Inhibitors by Free Energy Perturbation (FEP) Arjun Saha, Amy Y. Shih, Tara Mirzadegan, and Mark Seierstad Janssen Research and Development, 3210 Merryfield Row, San Diego, CA, 92121 Corresponding author phone: (858) 320-3361; email: [email protected]

ABSTRACT Since a goal of most drug discovery projects in either academia or industry is to design molecules that selectively bind to the desired protein, determination of protein-ligand binding free energies is of utmost importance in computer aided drug design. With the help of significant improvements in computer power, enhanced sampling techniques and accuracy of force fields, FEP (free energy perturbation) is becoming an important tool to estimate binding free energies in many drug discovery projects both retrospectively and prospectively. We have evaluated the ability of Schrödinger’s FEP+ to predict relative binding free energies of a congeneric series of noncovalent fatty acid amide hydrolase (FAAH) inhibitors using an in-house crystal structure. This study shows that although an impressively accurate correlation can be obtained with experimental IC50s considering small perturbations on the deeper side of the pocket, the same was not observed with small perturbations on the relatively more open-ended and solvent-accessible side of the pocket. To understand these observations, we thoroughly investigated several key factors including the sampling of asymmetrically substituted rings, different perturbation maps, impact of simultaneous perturbations at two different ends of the ligand, and selecting the perturbations in a “chemically sensible” way.

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 23

Introduction Potent binding of molecules to a target protein is one of the main goals of a typical drug discovery campaign.1 Several rigorous approaches can predict the difference in binding free energies of two structurally similar ligand molecules utilizing alchemical structural modifications,2 including free-energy perturbation (FEP) or thermodynamic integration (TI).2 The fundamentals of FEP date back several decades,3,4,5,6 although the impact of these approaches to industrial drug discovery projects has been very limited until recently.7,8,9,10 Requirements of long simulation times and accurate force fields have limited the use of this method in a fast-paced drug discovery setting.11 With the help of significant improvements in computer power, enhanced sampling techniques, low-cost parallel computing using GPUs and accuracy of force fields, FEP is becoming an important tool to predict relative binding free energies in many drug discovery projects.12 In this report, we have evaluated the ability of Schrodinger’s FEP+ to predict relative binding free energies of a noncovalent congeneric series of fatty acid amide hydrolase (FAAH)13,14 inhibitors using an in-house crystal structure. There are several literature reports on modeling FAAH inhibitors,15-26, 27 including the use of FEP.19,20 The serine hydrolase FAAH is a principal enzyme catalyzing the hydrolysis of certain bioactive lipids, including anandamide (AEA),14 an endogenous neurotransmitter that binds to cannabinoid receptors CB1 and CB2 which are expressed in the central nervous system.28 Inhibition of FAAH results in elevated amounts of AEA and other fatty acid amides. Several FAAH inhibitors have been reported in the literature; 14 many of these inhibit FAAH via covalent attachment to the catalytic serine. Janssen R&D carried out a high-throughput screening campaign on FAAH and observed an interesting hit compound (Figure 1a).29 Encouraged by the fact that all FAAH activity stems from the R-enantiomer, a detailed SAR exploration of substituted phenyl groups produced the inhibitors shown in Tables 1 and 2.29,

30

All the compounds in Table 1 have an identical

(unsubstituted) phenyl group depicted on the right, with variations only on the aryl ring on the left. These will be collectively called left-hand side (LHS) compounds. Likewise, Table 2 has righthand side (RHS) compounds with an identical chloro-trifluoromethyl phenyl attached to the pyrimidine on the left-hand side, and variants on the right. Note that ligand 31 is included in both Table 1 and Table 2; ligand 31 is the link between these two sets. All of these ligands are believed

2 ACS Paragon Plus Environment

Page 3 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

to bind in the FAAH active site non-covalently.29 IC50 values used in this calculation are from a functional assay using human FAAH.29

Figure 1: (a) Janssen HTS hit. (b) Compound used for crystallography study. An X-ray crystal structure of humanized rat FAAH was determined (PDB ID: ZZZZ) with the ligand in Figure 1b; the binding site is shown in Figure 2. The ligand is making hydrophobic interactions with Met-495, Leu-380 and Ile-407, π-π interactions with Phe-381, Phe-432 and Trp531, and an H-bond interaction with Thr-488. The LHS of the ligand is buried deep inside the pocket and the RHS is relatively open and solvent accessible.

3 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 23

Table 1: Left-hand side compounds considered in this report.

4 ACS Paragon Plus Environment

Page 5 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 2: Right-hand side compounds considered in this report.

Cpd 33

Ar =

IC50 (nM) 73 ± 23

Cpd 43*

Ar =

IC50 (nM) 500

34*

61 ± 7

44*

165

35*

110 ± 44

45*

28

36*

228 ± 79

46*

20

37*

271 ± 83

47*

>10,000

38*

57 ± 5

31

5.3 ± 3.2

39*

50 ± 19

40

14 ± 3

41*

59 ± 14

42*

33 ± 3

* Denotes racemate.

5 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 23

Figure 2: FAAH binding pocket with the crystal structure compound in Figure 1(b). Ligand is shown in green and the protein in grey. Key residues are labeled.

6 ACS Paragon Plus Environment

Page 7 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Method All FEP calculations were performed with version 2017-4 of the Schrödinger molecular modeling suite31 and its FEP+ framework. The FEP+ package includes GPU-enabled parallel molecular dynamics engine Desmond 3.9,32 with OPLS3 force field ,12 a cycle-closure correction algorithm33 to include redundant information into free energy estimates, and a REST-enhanced sampling technique.34,35 The protein structure was prepared in Maestro36 using the protein preparation wizard.37 All the maps were set up using FEP+ mapper38 to make perturbations between ligands. Missing force field parameters were calculated using the ffbuilder module.12 The following simulation parameters were used for FEP calculations: 12 lambda windows are simulated with the NPT ensemble, at 300 K and 1.01325 bar for 5ns and 10 ns. The IC50 of the racemic ligands was considered as half of the experimental IC50. MM-GBSA calculations were performed with Schrödinger's module;31 VSGB solvation model was used along with force field (OPLS3) minimization of the ligand. FEP input poses were prepared by alignment with the crystal structure ligand. For ligands with asymmetric substitution on either the left- or right-hand side phenyl, the selection of ring orientation was guided by the crystal structure; these considerations are illustrated in Figures 3 and 4. In Figure 3, ligand 16 (a LHS ligand) and ligand 37 (a RHS ligand) are shown aligned with the crystal structure ligand, which provides a precedent for alignment of the left phenyl ring; the 3- or 3,4-substituents are oriented analogous to the orientation of the dioxane. In addition to this crystallographic precedent, modeling of ligand binding into the FAAH structure suggests the same orientation for the left phenyl ring and also predicts a preferred orientation for the right phenyl ring. Figure 4a shows ligand 6 with the left phenyl ring in the opposite orientation, resulting in a clash with Leu-380. Similarly, in Figure 4b, ligand 37 has the right phenyl ring oriented such that the trifluoromethoxy group clashes with Leu-433. The preferred orientation of the right phenyl ring is indicated by the arrow, and shown in Figure 3.

7 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 23

Figure 3: Example of alignment of ligands with the crystal structure.

Figure 4: Protocol for orienting asymmetric rings for FEP+ input structures. (a) Position of the methyl group in ligand 16 is possible where one orientation (facing downward, at the position indicated by the arrow) is favored and the other (facing upward, shown explicitly) clashes with Leu 380 (red dashed lines indicate clash). (b) Position of the -OCF3 group in ligand 37 is possible where one orientation (facing downward, at the position indicated by the arrow) is favored and the other (facing upward, shown explicitly) clashes with Leu 433.

8 ACS Paragon Plus Environment

Page 9 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Results and discussion Initial FEP simulations used 47 ligands (Tables 1 and 2) with Schrödinger’s default mapper option (which uses Tanimoto similarity to evaluate the differences in a pair of molecules) to create perturbations between them. The correlation between experimental and calculated ΔG values is shown in Figure 5 (the map is shown in Figure S1); the correlation is very poor, with an R2 value of 0.30. The blue points in the plot indicate all the LHS ligands, while the green annotated points indicate all the RHS ligands. The correlation plot indicates some degree of separation between RHS and LHS ligands, with the RHS ligands consistently overpredicted in binding free energy compared to the LHS ones. We have classified the 84 edges in Figure 5 into two categories: a “pure edge” is defined as a perturbation between any two RHS or any two LHS ligands, whereas “cross edge” is defined as a perturbation between one RHS and one LHS ligand (any edge containing ligand 31 is a pure edge). The average ΔΔG error for the 66 pure edges is 0.89 kcal/mol, whereas the average ΔΔG error for the 18 cross edges is 3.23 kcal/mol (Table 3). Because cross edges involve simultaneous perturbations at two different ends of the ligand, they might be responsible for the poor results shown in Figure 5. We also calculated the ΔΔG errors associated with all the pure LHS edges and all the pure RHS edges; 66 pure edges consist of 52 pure LHS edges with an average ΔΔG error of 0.77 kcal/mol and 14 pure RHS edges having an average ΔΔG error of 1.33 kcal/mol (Table 3). With this analysis in mind, we next examined FEP maps that minimize or eliminate cross edges. Figure 6 shows the correlation between experimental and calculated ΔG values for a map (Figure S2) created biasing ligand 31, which differs from other ligands on only one phenyl ring. Even with this biased map (where a majority of the perturbations include ligand 31), RHS ligands are again overpredicted by FEP. Errors between FEP predicted and experimentally determined binding free energies are shown in Table 4, separated according to: (a) perturbation between ligand 31 and any LHS ligand and (b) perturbation between ligand 31 and any RHS ligand. FEP predictions are consistently worse for perturbations between ligand 31 and any RHS ligand; the average ΔΔG error for the 28 edges connecting ligand 31 and any LHS ligand is 0.65 kcal/mol, whereas the average ΔΔG error for the 13 edges connecting ligand 31 and any RHS ligand is 3.11 kcal/mol. Moreover, unlike with edges connecting ligand 31 and any LHS ligand (Table 4a), FEP

9 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 23

inaccurately predicts all the RHS ligands to be more potent than ligand 31 (FEP predicted values are negative in all the cases in Table 4b).

Figure 5: Correlation between experimental and calculated ΔG values for all 47 ligands using default FEP+ map (blue annotated points indicate the LHS ligands and green annotated points indicate the RHS ligands). The corresponding map (Figure S1) contains 84 edges. R2 value is 0.30. MUE and RMSE of 84 edges (ΔΔG) are 1.35 and 1.87 kcal/mol, respectively. The dark shaded area denotes 1 kcal/mol error range while the light shaded area denotes 2 kcal/mol error range.

10 ACS Paragon Plus Environment

Page 11 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 3: Analysis of 84 edges from Figure 5.

Types Total Pure Cross Pure left Pure right

Number of edges 84 66 18 52 14

MUE 1.35 0.89 3.23 0.77 1.33

Figure 6: Correlation between experimental and calculated ΔG values for a map biasing ligand 31 (blue annotated points indicate the LHS ligands and green annotated points indicate the RHS ligands). R2 value is 0.26. The MUE and RMSE are 1.37 and 1.86 kcal/mol, respectively. The corresponding map is shown in Figure S2. The dark shaded area denotes 1 kcal/mol error range while the light shaded area denotes 2 kcal/mol error range.

11 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 23

Table 4: Analysis of the 41 edges containing ligand 31 in the 'biased' map. (a) “Ligand1” refers to ligand 31 and “Ligand2” refers to LHS ligands. (b) “Ligand1” refers to ligand 31 and Ligand2 refers to RHS ligands. Experimental and predicted ΔΔG, and errors, are reported in kcal/mol. (a)

(b)

Ligand1 Ligand2 Exp. FEP Error 31 3 0.65 2.2 1.55 31 4 -0.03 1.33 1.36 31 21 2.25 3.6 1.35 31 2 2.62 3.78 1.16 31 25 -0.75 0.36 1.11 31 14 3.44 4.47 1.03 31 15 2.7 1.73 0.97 31 8 0.96 1.81 0.85 31 12 0.24 -0.58 0.82 31 23 0.62 -0.17 0.79 31 17 2.39 3.14 0.75 31 7 1.96 2.65 0.69 31 5 1.13 0.46 0.67 31 18 2.01 2.66 0.65 31 6 0.87 1.45 0.58 31 32 0.62 1.19 0.57 31 11 -0.83 -0.27 0.56 31 1 4.05 4.45 0.40 31 26 0.62 1 0.38 31 9 0.72 0.36 0.36 31 19 3.31 3.63 0.32 31 22 1.32 1.61 0.29 31 20 4.46 4.74 0.28 31 27 0.1 0.33 0.23 31 29 0.48 0.69 0.21 31 28 -0.3 -0.1 0.20 31 13 2.54 2.69 0.15 31 10 0.97 0.89 0.08 MUE 0.65

Ligand1 Ligand2 31 47 31 43 31 37 31 35 31 44 31 36 31 45 31 46 31 42 31 38 31 39 31 40 31 33 MUE

Exp. 4.05 2.28 1.92 1.38 1.62 1.81 0.57 0.37 0.67 0.99 0.92 0.57 1.14

FEP Error -1.54 5.59 -2.13 4.41 -2.26 4.18 -2.13 3.51 -1.83 3.45 -1.62 3.43 -2.71 3.28 -2.66 3.03 -2.11 2.78 -1.75 2.74 -1.41 2.33 -1.26 1.83 -0.45 1.59 3.11

12 ACS Paragon Plus Environment

Page 13 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

We next predicted the binding of RHS and LHS ligands separately by setting up their own maps within each set of ligands. The correlation plot for all the 32 LHS ligands (including ligand 31) is shown in Figure 7a. FEP predicts the relative binding energies of all the LHS ligands with high accuracy, having an R2 value of 0.79 (mean unsigned error (MUE) and root mean square error (RMSE) of all the edges (ΔΔG) are 0.67 and 0.91 kcal/mol, respectively). We observed similar results in extended simulations (10 ns, Figure S3). In a similar way, we also calculated the relative binding energies for all the RHS ligands (including ligand 31). The correlation plot in Figure 7b indicates that, unlike with LHS ligands, FEP cannot predict the relative binding energies of RHS ligands with high accuracy; the R2 value is 0.02 in a 5ns simulation. MUE and RMSE of these simulations are also worse than LHS ligands, and no improvement was observed with a 10ns simulation (Figure S4).

Figure 7: (a) Correlation between experimental and calculated ΔG values for all LHS ligands, 5ns simulation, R2 value is 0.79. MUE and RMSE are 0.67 and 0.91 kcal/mol, respectively. (b) Correlation between experimental and calculated ΔG values for all RHS ligands, 5ns simulation, R2 value is 0.02. MUE and RMSE are 1.31 and 1.64 kcal/mol, respectively. The dark shaded area denotes 1 kcal/mol error range while the light shaded area denotes 2 kcal/mol error range.

13 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 23

To understand the inability of FEP to predict binding free energies accurately for RHS ligands, we looked at the most representative structures for edges in Table 4. Figure 8a shows the most representative structures over the trajectory for an edge connecting ligand 31 (magenta) to ligand 1 (green). It suggests that there is no translational shifting of the binding mode of the most representative pose for the LHS ligands. On the other hand, Figure 8b shows a significant translational movement in the most representative structures for an edge connecting RHS ligand 40 to ligand 31. We wanted to investigate this translational movement to determine any possible correlation with the predicted binding free energies. Thus, we looked at the rest of the edges in Figure S2 and noticed that the presence or absence of translational movement was not consistent for RHS or LHS edges. To quantitatively evaluate the change in binding modes, we have calculated the average RMSD between binding poses of the ligand over the MD trajectory for all the edges in Figure S2. The average RMSD values for LHS ligands (2.40 ± 0.54Å, Table S1) are slightly more than RHS ligands (1.74A ± 0.37Å, Table S2), suggesting any translational movement is not the cause of the over prediction of RHS binding energies. We also calculated RMSD values between binding poses for the core region of the ligand (the area not included in the red boxes in Figure S5), only to reach at the same conclusion. Since the binding energy values (ΔG, which are calculated using “cycle closure”7,8) are highly dependent on the perturbation maps where different perturbation maps can lead to different ΔG values for the same ligands (Table S3), it is important to investigate the performance of RHS ligands from the ΔΔG values (Table 4b). The pattern we observed in Table 4b is that RHS ligands are all overpredicted, and the larger the perturbation, the higher the error in predicted ΔΔG values. This is not the case for LHS ligands, where ligands with large perturbations are predicted with high accuracy (Table 4a). These observations demonstrate that perturbation on the open-ended and solvent accessible side of the ligand can result in over prediction of the binding energies (ΔG).

14 ACS Paragon Plus Environment

Page 15 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 8: (a) Most representative structures (along the trajectory) of an edge connecting ligand 31 (magenta) and ligand 1 (green) in the 'biased' map (Figure S2). (b) Most representative structures (along the trajectory) of an edge connecting ligand 31 (magenta) and ligand 40 (green) in the 'biased' map (Figure S2). Sampling of asymmetrically substituted phenyl rings might be an important factor impacting the accuracy of our calculations for the RHS ligands. Thus, to investigate ligand conformation, we set up a smaller map to calculate the change in free energy between both the “upward” and “downward” facing Cl substitutions in ligand 31 and the crystal structure ligand. Regardless of the starting orientation of the Cl substitution in ligand 31, the FEP-predicted ΔΔG values between these pairs are fairly close; the ΔΔG value between “upward facing Cl substituted ligand 31” and the crystal structure ligand is 2.54 ± 0.11 kcal/mol and the ΔΔG value between “downward facing Cl substituted ligand 31” and the crystal structure ligand is 2.15 ± 0.11 kcal/mol. This suggests that FEP-predicted ΔΔG values are independent of the ligand starting conformations. 15 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 23

We next investigated the ligand conformation for all edges in the biased map. Throughout the 10 ns MD trajectory, the right-hand side aryl rings rarely deviated from the input conformation; this is consistent with modeling into the static crystal structure, which suggests a clash with Leu-433 in the flipped orientation. There was appropriate sampling on the left-hand side of the ligand in both 5 ns and 10 ns calculations. For example, in Figure 8a (the most representative structure of ligand 31, colored magenta), the orientation of the Cl substituted phenyl ring on the LHS is flipped relative to its input structure. Unlike in Figure 4, where Leu-380 clashes with the ligand, the same residue has now moved sufficiently in the most representative structure to facilitate the flip of this phenyl ring. This occurrence of flipping appears to be consistently true in most of the edges in the biased map, where ligand 31 has the Cl group facing downward in the FEP input structure but facing upward in the most representative structure. We also explored inclusion of asymmetrically substituted phenyl rings on both sides of the ligand in the “hot atoms” region36 (illustrated in Figure S5). In simulations with multiple durations (5, 10 and 20 ns), sampling was considerably improved on both sides of the ligand, with the RHS ligands now adequately sampling ring flipping. However, this strategy did not improve the prediction of binding energies for RHS ligands. We have also tried taking the entire ligand in the “hot atom” region, but again it did not improve FEPpredicted binding energies. To compare the performance of FEP, we also performed MM-GBSA calculations on these FAAH inhibitors. The correlation between experimental ΔG and MM-GBSA predicted ΔG values for all 47 ligands is shown in Figure 9 (Figure S6 shows the correlation for only LHS compounds). We observed that FEP-predicted binding free energies were significantly more accurate. In addition to the prediction of the experimental results, we were also interested to compare the rank order of the compounds in terms of binding free energies both in FEP and MM-GBSA. Table S4 shows the rank order for all 47 ligands and Table S5 shows the same for only LHS ligands. Although in a qualitative sense MM-GBSA has some predictive ability comparable to FEP (e.g. choosing the top ten compounds for synthesis), FEP is much more accurate in the quantitative prediction of ligand binding for these compounds.

16 ACS Paragon Plus Environment

Page 17 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 9: Correlation between MM-GBSA predicted and experimental ΔG values for all 47 ligands. R2 value is 0.10. Red and yellow denote LHS and RHS compounds, respectively. Conclusion The FEP method was not easily available at the time of this discovery project, and we designed this retrospective analysis to mimic the real-world use of FEP+ on these noncovalent FAAH inhibitors. The ability of FEP+ to predict relative binding free energies for compounds that differ in the deeper part of the FAAH binding pocket was excellent, and certainly accurate enough to provide a significant impact on a drug discovery project. This was true even though our choice of input orientation of the left-hand side aryl ring, guided by the only crystal structure in this series, may have been incorrect. Results with analogs that varied on the other end of the molecule were strikingly different. This part of the FAAH binding pocket is larger and solvent accessible, and all 17 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 23

inhibitors in this study that contained a substitution on that aryl ring were less potent than the unsubstituted parent. Yet FEP+ consistently overpredicted the potency of those ligands. This was true even with a variety of perturbation maps, simulation times, and sampling techniques. Supporting Information Supporting information available: Other FEP+ maps and correlation plots, RMSD of the ligand between binding poses over the MD trajectory, hot atom region illustration, additional MMGBSA results, and rank order comparison of MM-GBSA and FEP+ results. This information is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgement We thank Proteros Biostructures GmbH for the crystal structure. We thank Dr. Aaron Thompson for help with PDB submission and Dr. Nadeem Vellore for a critical review of the manuscript.

Notes The authors declare no competing financial interest.

18 ACS Paragon Plus Environment

Page 19 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

For Table of Contents Only

19 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 23

References

1

Jorgensen, W. L. Efficient Drug Lead Discovery and Optimization. Acc. Chem. Res. 2009, 42,

724−733. 2

Free Energy Calculations: Theory and Applications in Chemistry and Biology; Chipot, C.,

Pohorille, A., Eds.; Springer Series in Chemical Physics, Vol. 86; Springer: Berlin and Heidelberg, Germany, 2007. 3

McCammon, J. A.; Gelin, B. R.; Karplus, M. Dynamics of Folded Proteins. Nature 1977, 267,

585−590. 4

Jorgensen, W. L.; Ravimohan, C. Monte Carlo Simulation of Differences in Free Energies of

Hydration. J. Chem. Phys. 1985, 83, 3050−3054. 5

Wong, C. F.; McCammon, J. A. Dynamics and Design of Enzymes and Inhibitors. J. Am. Chem.

Soc. 1986, 108, 3830−3832. 6

Merz, K. M.; Kollman, P. A. Free Energy Perturbation Simulations of the Inhibition of

Thermolysin: Prediction of the Free Energy of Binding of a New Inhibitor. J. Am. Chem. Soc. 1989, 111, 5649−5658. 7

Abel, R.; Wang, L.; Harder, E. D.; Berne, B. J.; Friesner, R. A. Advancing Drug Discovery

through Enhanced Free Energy Calculations. Acc. Chem. Res. 2017, 50, 1625−1632. 8

Keränen, H.; Perez-Benito, L.; Ciordia, M.; Delgado, F.; Steinbrecher, T. B.; Oehlrich, D.; van

Vlijmen, H. W. T.; Trabanco, A. A.; Tresadern, G. Acylguanidine Beta Secretase 1 Inhibitors: A Combined Experimental and Free Energy Perturbation Study. J. Chem. Theory Comput. 2017, 13, 1439–1453. 9

Kuhn, B.; Tichý, M.; Wang, L.; Robinson, S.; Martin, R. E.; Kuglstatter, A.; Benz, J.; Giroud,

M.; Schirmeister, T.; Abel, R.; Diederich, F.; Hert, J. Prospective Evaluation of Free Energy Calculations for the Prioritization of Cathepsin L Inhibitors. J. Med. Chem. 2017, 60, 2485–2497. 10

Perez-Benito, L.; Keränen, H.; van Vlijmen, H. W. T.; Tresadern, G. Predicting Binding Free

Energies of PDE2 Inhibitors. The Difficulties of Protein Conformation. Scientific Reports, 2018, 8, 1−10. 11

Deng, Y.; Roux, B. Computations of Standard Binding Free Energies with Molecular Dynamics

Simulations. J. Phys. Chem. B 2009, 113, 2234−2246.

20 ACS Paragon Plus Environment

Page 21 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

12

Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J. Y.; Wang, L.; Lupyan, D.;

Dahlgren, M. K.; Knight, J. L.; Kaus, J. W.; Cerutti, D. S.; Krilov, G.; Jorgensen, W. L.; Abel, R.; Friesner, R. A. OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins. J. Chem. Theory Comput. 2016, 12, 281−296. 13

Cravatt, B. F.; Demarest, K.; Patricelli, M. P.; Bracey, M. H.; Giang, D. K.; Martin, B. R.;

Lichtman, A. H. Supersensitivity to Anandamide and Enhanced Endogenous Cannabinoid Signaling in Mice Lacking Fatty Acid Amide Hydrolase. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 9371–9376. 14

Breitenbucher, J. G.; Seierstad, M. Discovery and Development of Fatty Acid Amide Hydrolase

(FAAH) Inhibitors. J. Med. Chem. 2008, 51, 7327−7343. 15

Zhao, D. S.; Wang, H. Y.; Lian, Z. H.; Han, D. X.; Jin, X.; Pharmacophore Modeling and Virtual

Screening for the Discovery of New Fatty Acid Amide Hydrolase Inhibitors. Acta Pharmaceutica Sinica B 2011, 1, 27−35. 16

Lodola, A.; Capoferri, L.; Rivara, S.; Tarzia, G.; Piomelli, D.; Mulholland, A. J.; Mor, M.

Quantum Mechanics/Molecular Mechanics Modeling of Fatty Acid Amide Hydrolase Reactivation Distinguishes Substrate from Irreversible Covalent Inhibitors. J. Med. Chem. 2013, 56, 2500−2512. 17

Bowman, A. L.; Makriyannis, A. Approximating Protein Flexibility through Dynamic

Pharmacophore Models: Application to Fatty Acid Amide Hydrolase (FAAH). J. Chem. Inf. Model. 2011, 51, 3247–3253. 18

Palermo, G.; Favia, A. D.; Convertino, M., De Vivo, M. The Molecular Basis for Dual Fatty

Acid Amide Hydrolase (FAAH)/Cyclooxygenase (COX) Inhibition. ChemMedChem 2016, 11, 1252–1258. 19

Guimaraes, C. R. W.; Boger, D. L.; Jorgenson, W. L. Elucidation of Fatty Acid Amide Hydrolase

Inhibition by Potent ɑ-Ketoheterocycle Derivatives from Monte Carlo Simulations. J. Am. Chem. Soc. 2005, 127, 17377−17384. 20

Brohman, I. T.; Acevedo, O.; Jorgenson, W. L. Elucidation of Hydrolysis Mechanisms for Fatty

Acid Amide Hydrolase and Its Lys142Ala Variant via QM/MM Simulations. J. Am. Chem. Soc. 2006, 128, 16904−16913.

21 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

21

Page 22 of 23

Saario, S. M.; Poso, A.; Juvonen, R. O.; Jarvinen, T.; Salo-Ahen, O. M. H. Fatty Acid Amide

Hydrolase Inhibitors from Virtual Screening of the Endocannabinoid System. J. Med. Chem. 2006, 49, 4650−4656. 22

Acevedo, O.; Jorgenson, W. L. Advances in Quantum and Molecular Mechanical (QM/MM)

Simulations for Organic and Enzymatic Reactions. Acc. Chem. Res. 2010, 43, 142−151. 23

Lodola, A.; Mor, M.; Rivara, S.; Christov, C.; Tarzia, G.; Piomelli, D.; Mulholland, A. J.

Identification of Productive Inhibitor Binding Orientation in Fatty Acid Amide Hydrolase (FAAH) by QM/MM Mechanistic Modelling. Chem. Commun. 2008, 2, 214–216. 24

Palermo, G.; Branduardi, D.; Masetti, M.; Lodola, A.; Mor, M.; Piomelli, D.; Cavalli, A.; De

Vivo, M. Covalent Inhibitors of Fatty Acid Amide Hydrolase: A Rationale for the Activity of Piperidine and Piperazine Aryl Ureas. J. Med. Chem. 2011, 54, 6612–6623. 25

Lodola, A.; Mor, M.; Zurek, J.; Tarzia, G.; Piomelli, D.; Harvey, J. N.; Mulholland, A. J.

Conformational Effects in Enzyme Catalysis: Reaction via a High Energy Conformation in Fatty Acid Amide Hydrolase. Biophys. J. 2007, 92, L20– L22. 26

Lodola, A.; Mor, M.; Hermann, J. C.; Tarzia, G.; Piomelli, D.; Mulholland, A. J. QM/MM

Modelling of Oleamide Hydrolysis in Fatty Acid Amide Hydrolase (FAAH) Reveals a New Mechanism of Nucleophile Activation. Chem. Commun. 2005, 0, 4399–4401. 27

Palermo, G.; Rothlisberger, U.; Cavalli, A.; De Vivo, M. Computational insights into function

and inhibition of fatty acid amide hydrolase. Eur. J. Med. Chem. 2015, 91, 15−26. 28

Willoughby, K. A.; Moore, S. F.; Martin, B. R.; Ellis, E. F. The Biodisposition and Metabolism

of Anandamide in Mice. J. Pharmacol. Exp. Ther. 1997, 282, 243−247. 29

Keith, J. M.; Hawryluk, N.; Apodaca, R. L.; Chambers, A.; Pierce, J. M.; Seierstad, M.;

Palmer, J. A.; Webb, M.; Karbarz, M. J.; Scott, B. P.; Wilson, S. J.; Luo, L.; Wennerholm, M. L.; Chang, L.; Rizzolio, M.; Chaplan, S. R.; Breitenbucher, J. G. 1-Aryl-2-((6-aryl)pyrimidin-4yl)amino)ethanols as Competitive Inhibitors of Fatty Acid Amide Hydrolase. Bioorg. Med. Chem. Lett. 2014, 24, 1280–1284. 30

Apodaca, R.; Breitenbucher, J. G.; Chambers, A. L.; Xiaohu, D.; Hawryluk, N. A.; Keith, J. M.;

Neelakhandha, M. S.; Merit, J. E.; Pierce, J. M.; Seierstad, M.; Wei, X. Aryl HydroxyethylaminoPyrimidines and Triazines as Modulators of Fatty Acid Amide Hydrolase. U. S. Patent WO 2009/105220 A1, August 27, 2009. 31

Small-Molecule Drug Discovery Suite 2017-4, Schrödinger, LLC, New York, NY, 2018. 22 ACS Paragon Plus Environment

Page 23 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

32

Bowers, K. J.; Chow, E.; Xu, H.; Dror, R. O.; Eastwood, M. P.; Gregersen, B. A.; Klepeis, J. L.;

Kolossvary, I.; Moraes, M. A.; Sacerdoti, F. D.; Salmon, J. K.; Shan, Y.; Shaw, D. E. Scalable Algorithms for Molecular Dynamics Simulations on Commodity Clusters. Proceedings of the ACM/IEEE Conference on Supercomputing (SC06) 2006, 43. 33

Wang, L.; Deng, Y.; Knight, J. L.; Wu, Y.; Kim, B.; Sherman, W.; Shelley, J. C.; Lin, T.; Abel,

R. Modeling Local Structural Rearrangements Using FEP/REST: Application to Relative Binding Affinity Predictions of CDK2 Inhibitors. J. Chem. Theory Comput. 2013, 9, 1282−1293. 34

Liu, P.; Kim, B.; Friesner, R. A.; Berne, B. J. Replica Exchange with Solute Tempering: A

Method for Sampling Biological Systems in Explicit Water. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 13749−13754. 35

Wang, L.; Friesner, R. A.; Berne, B. J. Replica Exchange with Solute Scaling: A More Efficient

Version of Replica Exchange with Solute Tempering (REST2). J. Phys. Chem. B 2011, 115, 9431−9438. 36 37

Schrödinger LLC. Maestro, 2017-4. http://www.schrodinger.com/Maestro/. Madhavi Sastry, G. M.; Adzhigirey, M.; Day, T.; Annabhimoju, R.; Sherman, W. Protein and

Ligand Preparation: Parameters, Protocols, and Influence on Virtual Screening Enrichments. J. Comput.-Aided Mol. Des. 2013, 27, 221−234. 38

Liu, S.; Wu, Y.; Lin, T.; Abel, R.; Redmann, J. P.; Summa, M. C.; Jaber, V. R.; Lim, N. M.;

Mobley, D. L. Lead Optimization Mapper: Automating Free Energy Calculations for Lead Optimization. J. Comput.-Aided Mol. Des. 2013, 27, 755−770.

23 ACS Paragon Plus Environment