Article pubs.acs.org/JCTC
Acylguanidine Beta Secretase 1 Inhibitors: A Combined Experimental and Free Energy Perturbation Study Henrik Keran̈ en,† Laura Pérez-Benito,‡,§ Myriam Ciordia,⊥ Francisca Delgado,⊥ Thomas B. Steinbrecher,∥ Daniel Oehlrich,# Herman W. T. van Vlijmen,† Andrés A. Trabanco,⊥ and Gary Tresadern*,§ †
Computational Chemistry, Janssen Research & Development, Janssen Pharmaceutica N. V., Turnhoutseweg 30, B-2340 Beerse, Belgium ‡ Laboratori de Medicina Computacional Unitat de Bioestadistica, Facultat de Medicina, Universitat Autonoma de Barcelona, 08193, Bellaterra, Spain § Computational Chemistry, Janssen Research and Development, Janssen-Cilag, c/ Jarama 75A, 45007, Toledo, Spain ⊥ Neuroscience Medicinal Chemistry, Janssen Research and Development, Janssen-Cilag, c/ Jarama 75A, 45007, Toledo, Spain ∥ Schrödinger GmbH, Dynamostrasse 13, 68165 Mannheim, Baden-Württemberg, Germany, # Neuroscience Medicinal Chemistry, Janssen Research & Development, Janssen Pharmaceutica N. V., Turnhoutseweg 30, B-2340 Beerse, Belgium S Supporting Information *
ABSTRACT: A series of acylguanidine beta secretase 1 (BACE1) inhibitors with modified scaffold and P3 pocket substituent was synthesized and studied with free energy perturbation (FEP) calculations. The resulting molecules showed potencies in enzymatic BACE1 inhibition assays up to 1 nM. The correlation between the predicted activity from the FEP calculations and the experimental activity was good for the P3 pocket substituents. The average mean unsigned error (MUE) between prediction and experiment was 0.68 ± 0.17 kcal/mol for the default 5 ns lambda window simulation time improving to 0.35 ± 0.13 kcal/mol for 40 ns. FEP calculations for the P2′ pocket substituents on the same acylguanidine scaffold also showed good agreement with experiment and the results remained stable with repeated simulations and increased simulation time. It proved more difficult to use FEP calculations to study the scaffold modification from increasing 5 to 6 and 7 membered-rings. Although prediction and experiment were in agreement for short 2 ns simulations, as the simulation time increased the results diverged. This was improved by the use of a newly developed “Core Hopping FEP+” approach, which also showed improved stability in repeat calculations. The origins of these differences along with the value of repeat and longer simulation times are discussed. This work provides a further example of the use of FEP as a computational tool for molecular design.
■
applications are emerging,12−18 including work from our laboratories investigating FEP applied in lead optimization.19,20 Here, we explore FEP predictions of the binding energies of β-secretase 1 (BACE1) inhibitors in the different subpockets of the active site. Cleavage of amyloid precursor protein by this aspartyl protease leads to increases of β-amyloid (Aβ) peptides that aggregate forming senile plaques, one of the neuropathological features in Alzheimer’s disease (AD).21 Contemporary BACE1 inhibitors contain an amidine/guanidine moiety within a heterocycle of varying size.22 This protonated group forms an optimal hydrogen-bonding network to the catalytic aspartate dyad, Figure 1. Also, a quaternary alpha sp3 carbon provides an ideal vector to fill the P1−P3 and P2′ pockets of the binding site (Figures 1 and 2).23,24 These molecules have improved drug like properties and multiple examples are in clinical trials.25 Given the huge pharmaceutical interest in finding new treatments for AD, BACE1 is a very well-studied
INTRODUCTION
The accurate prediction of protein−ligand binding affinities is of major interest.1 Rigorous approaches can calculate the binding free energy difference between two structurally similar ligands by making use of alchemical structural modifications. Free-energy perturbation (FEP) or thermodynamic integration (TI), using molecular dynamics (MD) or Monte Carlo simulations, are among the widely used approaches.2 These methods are ideal for drug discovery lead optimization where close structural analogues are synthesized and their properties compared to previous best leads. Computation of accurate relative binding affinities can make a big impact in this costly phase. Also, it avoids the computationally expensive prediction of absolute binding free energies. Calculating relative protein− ligand binding affinities in this way dates back at least thirty years.3−8 Lately, new sampling algorithms, improved force fields and low-cost parallel computing (often graphics processing units GPU), have improved accuracy and turnaround time.9−11 Reports of large scale and industrial © 2017 American Chemical Society
Received: November 22, 2016 Published: January 19, 2017 1439
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation
molecular design, where an error below 1 kcal/mol is preferred.26 FEP binding energy predictions are performed for three different regions of the BACE1 active site: the P3 and P2′ subpockets as well as computationally more challenging scaffold modifications. Amidine/guanidine heterocyclic scaffolds of different sizes provide leads with varying potency and physicochemical properties.27 We focus on acylguanidines (i− iii) that display moderate pKa values (calc. basic pKa 7.5) important to achieve drug like BACE1 inhibitors.28−30 The trends in activity are studied via the systematic synthesis, screening and computational evaluation of the 5, 6 and 7 membered-ring molecules along with further examination with FEP. We look in detail at the FEP performance comparing with other methods and extending, repeating and averaging the FEP results. The simulation times were selected ranging from 2 up to 40 ns per lambda window of each perturbation. We also compare the standard FEP+ implementation to the recently developed “Core Hopping FEP+” on scaffold perturbations, showing beneficial application of FEP beyond its traditional application for close analogues in lead optimization. The extensive computational study and analysis provides further evidence of the attractiveness of FEP and MD techniques for drug discovery but also highlights some of the challenges. In particular, we show that additional sampling is important when the flexibility of the ligand-protein complex needs to be further explored to achieve accurate calculated binding free energies.
Figure 1. Scaffold structures of acylguanidine BACE1 inhibitors (i− iii). Schematic binding mode of amidine BACE1 inhibitors and general structure for molecules in this study.
■
METHOD Computational Docking Procedure. To select appropriate initial protein coordinates all BACE1 protein data bank (PDB) structures were searched to identify those containing ligands most similar to this study. Thus PDB ID 3ZOV31 that had an acylguanidine (ii) and PDB 3IN432 containing an aminohydantoin (i) were identified. The former was used for P3 and scaffold pocket perturbations and the latter for the P2′ pocket perturbations. The protein structure was prepared as follows. First, 3ZOV (P3 and scaffold pocket perturbations) and 3IN4 (P2′ pocket) were imported into Maestro33 and structure preparation was performed using the Protein Preparation Wizard34 with default settings to fix missing sidechains/atoms, missing loops were modeled, protein protonation states were assigned with PROPKA, the hydrogen bonding network was optimized, and ligand charges assigned with the OPLSv3 force field. To relieve local clashes a restrained minimization was performed with a 0.5 Å heavy-atom RMSD displacement cutoff, below which the minimization is terminated. The catalytic aspartates were treated in their ionized states. FEP starting geometries for the P3 and scaffold perturbations were chosen via docking. For docking Glide software (release 2015-4) from Schrödinger was used. The ligand crystalized in the 3ZOV was used to place the grid box. All active site waters were retained and treated as part of the receptor grid. Two Hbonds constraints on the Asp from the catalytic center (Asp 289 and Asp 93, often referred to as Asp228 and Asp32, respectively) were chosen to perform the docking. All ligand molecules were prepared for docking using the LigPrep tool. All default settings were used except ionization was explicitly set to ensure all ligands were protonated on the guanidine ring. ConfGen35 with fast settings was used to derive multiple 3D conformers for each molecule involving P3 and scaffold region perturbations (molecules 8a−g, 17a−g, and 27a−g). All
Figure 2. Docked binding mode of 17d showing important features of the BACE1 binding site such as the active site flap (green), 10s loop (blue), S1, S3, and S2′ subpockets, as well as key amino acids.
target and to some extent a test case for medicinal chemistry and molecular design aspects. We recently reported20 that FEP was capable of predicting the binding energies for substituents in the P3 subpocket to an accuracy of 0.58 and 0.91 kcal/mol with 5 ns FEP simulations for retrospective and prospective studies respectively. This improved to 0.57 and 0.59 kcal/mol with calculations performed with 20 ns simulation times. These data are well within the required accuracy to have an impact on 1440
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation
Figure 3. FEP results for the scaffold perturbations (legend in panel h): Each panel shows the same P3 pocket R-group (a−g, in Table 1). Simulations were performed for 2 (orange), 5 (blue), 10 (purple), and 40 (green) ns. The binding energies are normalized with respect to the least active of the three molecules (always the 5 membered-ring example) with experimental values shown in dark and light gray. The calculated binding energies are the average of three repeats and error bars show the standard deviations.
(Figures 3, S4, and S5) and define the perturbation between molecules. The scaffold perturbations for instance, perturbed between each 5, 6, and 7 membered-ring. Default simulation parameters were used: 12 lambda windows are simulated with the NPT ensemble, at 300 K and 1.01325 bar (for more details on FEP+ methodology see Supporting Information). Simulation time was varied with 2, 5, 10, 20, and 40 ns lambda window simulation lengths being performed for ligands both in the complex and in solution. Different time lengths were run independently, that is, they were not restarted from shorter simulations. All the simulations were repeated three times to monitor system stability. Each repeat used a different random seed to assign alternative random velocities at the start of each MD simulation. Starting conformations for the FEP calculations were chosen as described above based on docking for the 3ZOV BACE1 crystal structure case (P3 and scaffold pocket) and in a manually placed orientation for the P2′ perturbations in the 3IN4 crystal structure. For the P2′ pocket FEP two possible dihedral orientations were considered (Tables S4 and S5). Molecules were treated in an ionized form and missing force field parameters were calculated using the ffbuilder module. Proteins were used as prepared above, but all resolved crystal water molecules were retained for the free energy simulations. We report theoretical error estimates based on standard deviation between repeated simulations, and also the coefficient of determination (R2) and the mean unsigned error compared to experiment. Experimental binding free energies (dGexp) were computed using IC50 values by
conformers were then passed as input to the Glide XP docking thereby producing multiple docking solutions for each conformer of each molecule. The Glide XP scoring function was used, but sampling was increased through modifying a number of parameters within Glide: expanded sampling was turned on, and 15 initial poses were passed to postdocking minimization. All other docking parameters were set to the defaults. Results were then aggregated at the level of each molecule and the best poses inspected and chosen, inputs are provided in Supporting Information. Ligand molecules involving P2′ perturbations (molecules 28−39) were prepared as above but instead of docking their starting poses were determined by superimposing to the ligand in the 3IN4 crystal structure. Two possible orientations of P2′ substituents were considered (Table S4 and S5). Inputs are provided in Supporting Information. For comparison of performance of FEP with other methods, we used the Glide XP scoring and MM-GBSA calculations. The Glide XP scores were used from the protocol described above. The MM-GBSA for the P3 and scaffold pockets were run using the same Glide XP docking poses that were also used as input for the FEP calculations, whereas for the P2′ the same superimposed structures were used. The VSGB solvation model was used along with force field minimization of the ligand and a default approach with no protein minimization and a second approach with minimization of an 8 Å radius of the surrounding binding site (using the same active region for all ligands). Computational Free Energy Perturbation Procedure. All calculations were conducted using version 2015-4 of the Schrodinger molecular modeling suite and their FEP+ implementation, except “Core Hopping FEP+”, which used version 2016-3. The FEP+ product combines an accurate modern force field, OPLSv3,36 the efficient GPU-enabled parallel molecular dynamics engine Desmond, version 3.9,37 the REST enhanced sampling technique38,39 and the cycle-closure correction algorithm40 to incorporate redundant information into free energy estimates. Calculations were conducted using the FEP+ mapper technology41 to automate setup and analysis
dGexp ≈ −RT ln(IC50)
(1)
where R is the gas constant and T is temperature in K. This introduces a constant offset for all dGexp values computed from IC50 data compared to hypothetical values obtained from Ki data. As assay conditions are the same for all compounds, the following procedure ensures this constant has no effect on our results. The ddGFEP values for each compound are output from the calculations with respect to a reference compound within 1441
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation Scheme 1a
a
(a) TMSCN, NH4Cl, NH3 in MeOH, rt, 4 days, 82%. (b) HCl, reflux, 18 h, 63%. (c) SOCl2, MeOH, reflux, 18 h, 63%. (d) Chiral SFC purification. (e) CSCl2, DIPEA, DCM, rt, 18 h. (f) CH3NH2, DIPEA, ACN, 90 C, 3 h, 47% (two steps). (g) TBHP, NH3 aq., MeOH, rt, 7 h, 94%. (h) NaN3, CuI, Na2CO3, N,N′-dimethylethylenediamine, 110 C, 3 h, 64%. (i) ArCOOH, DMTMM, MeOH, rt, 2-26 h, 36−95%.
Scheme 2a
a (a) (R)-tert-butanesulfinamide, Ti(iPro)4, THF, 70 C, 24 h, 74%. (b) CH3COOCH3, LDA, TiCl(iPro)3, THF, −78 C, 4 h, 65%. (c) HCl (4M in 1,4-dioxane), MeOH, rt, 2 h, 83%. (d) Tert-butyl N-(methylcarbamothioyl)carbamate, EDCI, DIPEA, DMF, rt, 5 h, 78%. (e) H2, Pd/C, THF, rt, 5 h, 70%. (f) HNO3/H2SO4, −10 C, 15 min. Then, NaHCO3, THF, Boc2O, rt, 24 h, 69%. (g) H2, Pd/C, THF, 40 C, H-Cube reactor, 49%. (h) ArCOOH, DMTMM, MeOH, rt, 3-24 h, 15-73%. (i) HCl (4 M in 1,4-dioxane) or TFA in DCM, rt, 6-48 h, 41−99%.
the corresponding (2-fluoro-5-bromo-phenyl)-ethanone (1) in 82% yield by the Strecker reaction. Hydrolysis followed by esterification gave α-aminoester 3, which was purified by chiral separation to yield the desired enantiomer. Reaction with thiophosgene and subsequent cyclization with methylamine gave thioacylguanidine (R)-5. Guanidine (R)-6 was synthesized via reaction with tert-butyl hydroperoxide in aqueous ammonia. Finally, an Ullman reaction allowed access to the aniline (R)-7, that was coupled with different carboxylic acids using 4-(4,6dimethoxy-2,3,5-triazin-2-yl)-4-methylmorpholinium chloride
each mapper, which is assigned as ddGFEP = 0. The corresponding dGFEP were then generated by normalizing the ddGFEP to the mean of the dGexp according to eq 2. ⎛ ∑ ddG ∑ dGexp ⎞ FEP ⎟⎟ dGFEP = ddGFEP − ⎜⎜ − n n ⎝ ⎠
(2)
As such, mean unsigned errors and correlation were compared between dGexp and dGFEP. Synthesis. The synthesis of the 5 membered cycles is depicted in Scheme 1. The α-aminonitrile 2 was prepared from 1442
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation Scheme 3a
a
(a) NaBH4, MeOH, rt, 2 h, 96%. (b) Boc2O, NaHCO3, THF, rt, 16 h, 95%. (c) Dess-Martin, DCM, rt, 1 h, 93%. (d) Triethylphosphonoacetate, NaH, THF, rt, 3 h, 89%. (e) H2, Pt/C, EtOH, 35 C, H-Cube, 95%. (f) TFA, DCM, rt, 1 h, 40%. (g) p-anisaldehyde, NaBH4, EtOH, rt, 18 h, 41%. (h) Tert-butyl N-(methylcarbamothioyl)carbamate, EDCI, DIPEA, DMF, rt, 22 h, 96%. (i) NaOH, water, EtOH, 50 C, 2 h. (j) HATU, DIPEA, DMF, rt, 18 h, 54% (two steps). (k) NaN3, CuI, Na2CO3, N,N′-dimethylethylenediamine, DMSO, 110 C, 2 h, 74%. (l) H2, Pd/C, EtOH, rt, 16 h, 86%. (m) ArCOOH, DMTMM, MeOH, rt, 3−5 h, 73−92%. (n) HCl 4 M in dioxane or TFA in DCM, rt, 4−6 h, 31−99%.
(DMTMM) as the coupling agent, obtaining final compounds (R)-8a-g enantiomerically pure. Six-membered derivatives were prepared as shown in Scheme 2. Condensation of 1 with chiral (R)-tert-butanesulfinamide in the presence of Ti(OiPr)4 gave sulfinyl imine (R)-9.42 The Nsulfinyl activates the CN bond for the diastereoselective nucleophilic addition of methyl acetate and was deprotected to obtain (S)-11. Reaction with tert-butyl N(methylcarbamothioyl)carbamate gave the acylguanidine (S)12.43 Substitution of the bromine by an amino group failed as such the aniline was introduced by removing the bromine atom by hydrogenation. The F-substituent in the aryl group allowed the regioselective p-nitration. After neutralization, the scaffold was protected in situ with a Boc group yielding (S)-14. Reduction to the corresponding amine gave (S)-15, which was coupled with the carboxylic acids. Final compounds (S)-17a-g were obtained as pure enantiomers and hydrochloride or TFA salts after deprotection of the Boc group in acidic conditions. The α-aminoester (R)-3 was the starting point for the synthesis of the 7 membered cycles, Scheme 3). Reduction, followed by reaction with Boc anhydride and subsequent oxidation furnished aldehyde (R)-19. A Wittig−Horner− Wadsworth-Emmons reaction provided the E-double bond species which was hydrogenated, providing amino ester (S)-20. Deprotection followed by reprotection with a p-methoxybenzyl group, then condensation with tert-butyl N-(methylcarbamothioyl)carbamate yielded (S)-22. Cyclization failed so the ester was hydrolyzed to the carboxylic acid and the intramolecular peptidic coupling was performed using HATU.
Then, synthesis of the aniline using sodium azide, followed by removal of the p-methoxybenzyl gave acylguanidine (S)-25. Coupling of this intermediate with the carboxylic acids and deprotection of the Boc group gave final compounds (S)-27a-g as hydrochloride or TFA salts and enantiomerically pure. An example 8 membered scaffold was synthesized, it had never been reported in literature. Since the bromo in the para position to the fluoro proved problematic, we developed synthesis of the 8 membered scaffold with the aniline in place (Scheme 4). Acetophenone 28 was condensed with the Ellman’s reagent as described in scheme 2, followed by nucleophilic addition to yield 29. During this reaction, the trifluorocarbamate was deprotected. To continue, the aniline was protected with a carboxybenzyl (Cbz) group. The sulfinil moiety was deprotected in acidic conditions and the amino group was protected with a p-methoxybenzyl to yield 30. Condensation with tert-butyl N-(methylcarbamothioyl)carbamate, followed by addition of acryloyl chloride gave 32, which was cyclized by ring closure metathesis using Grubbs catalyst second generation to yield 33. Removal of both protecting groups (p-methoxybenzyl and Cbz) and the double bond was achieved by hydrogenation using Pd/C as catalyst. Coupling of 34 with a carboxylic acid and deprotection of the Boc group gave final compound 36 as a racemate. Enzymatic BACE1 assay: BACE1 enzymatic activity was assessed by a FRET assay using an amyloid precursor protein (APP) derived 13 amino acids peptide that contains the ‘Swedish’ Lys-Met/Asn-Leu mutation of the APP beta-secretase cleavage site as a substrate (Bachem cat No. M-2465) and 1443
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation Scheme 4a
(a) Ellman’s reagent, Ti(iPro)4, THF, 70 C, 36 h, 49%. (b) Allylmagnesium bromide, THF, −78 C, 1 h, 51%. (c) Benzyl chloroformate, NaHCO3 aq., THF, rt, 1 h, 82%. (d) HCl 4N dioxane, MeOH, rt, 1 h, 58%. (e) p-anisaldehyde, NaBH4, EtOH, 40 C, 18 h, 46%. (f) Tert-butyl N(methylcarbamothioyl)carbamate, EDCI, DIPEA, DMF, rt, 24 h, 91%. (g) Acryloyl chloride, DIPEA, DCM, rt, 1 h, 81%. (h) Grubbs 2nd generation, DCM, 50 C, 16 h, 71%. (i) H2, Pd/C, AcOEt, rt, 18 h, 50%. (j) ArCOOH, DMTMM, MeOH, rt, 3 h, 98%. (k) HCl 4 M in dioxane, rt, 2 h, 61%. a
perturbations between 3 molecules). 2) P3 pocket perturbations between 7 different R-groups (a to g in Table 1) with the same scaffold (i.e., 3 separate FEP mapper simulations of perturbations between 7 molecules). 3) P2′ pocket perturbations between 12 different substituents on the same 5 membered-ring scaffold (i.e., 1 FEP mapper simulation of perturbation between 12 molecules). The biological activity of the molecules and then the computational results are presented. The 21 molecules 8a−8g, 17a−17g, and 27a−27g from the three different scaffold sizes were successfully synthesized and screened for BACE1 inhibitory activity in both enzymatic and cellular assays, Table 1.44 The 6 membered cycles were most active, IC50 below 10 nM in both assays, followed by the 7 membered and finally the 5 membered-rings. The range of enzymatic activity for the P3 group substituents was largest for the 5 and 7 membered cycles, between 26 nM and 3.4 μM for the 5 membered-ring and 7 to 200 nM in the 7 membered cycle but only 1 to 10 nM for the 6 membered-ring examples. The activity of P3 pocket R-groups were not ranked the same for the three scaffolds, however, there were trends. The oxazole a and the 2-(methoxyethoxy)pyrazine f were among the least active while the alkyne derivative g was among the most potent reaching subnanomolar cellular IC50 for the 6 membered cycle. Most of the 7 membered cyclic derivatives showed lower cellular than enzymatic activity, which may be due to difficulties with cell permeation for the larger scaffold. An 8 membered-
soluble BACE1(1−454) (Aurigene, Custom made). This substrate contains two fluorophores, (7-methoxycoumarin-4yl) acetic acid (Mca) is a fluorescent donor with excitation wavelength at 320 nm and emission at 405 nm and 2,4dinitrophenyl (Dnp) is a proprietary quencher acceptor. The increase in fluorescence is linearly related to the rate of proteolysis. In a 384-well format, BACE1 is incubated with the substrate and the inhibitor. The amount of proteolysis is directly measured by fluorescence measurement in the Fluoroskan microplate fluorometer (Thermo scientific). For the low control no enzyme was added to the reaction mixture.
■
RESULTS Planning FEP Calculations, Synthesis, and Screening. We designed acylguanidines with 5, 6, and 7 membered-ring scaffolds and P3 R-group substituents similar to those reported in BACE inhibitors, Table 1.22 The P2′ was explored using the 5 membered-ring acylaguanidine (aminohydantoin, i) and different substituents taken from Malamas et al. 24 A combination of docking and manual placement were used to define the input for the FEP calculations. An example of 17d docked into BACE1 and describing the binding site features is shown in Figure 2. The FEP calculations were performed in three ways. (1) Scaffold perturbations between the 5, 6 and 7 membered-ring within groups of three molecules containing the same P3 R-group (i.e., 7 separate FEP mapper simulations of 1444
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation Table 1. Enzymatic and Cellular Activity of Acylguanidine BACE1 Inhibitorsa
n
comp. n = 0; 5 membered-rings 8a 8b 8c 8d 8e 8f 8g n = 1; 6 membered-rings 17a 17b 17c 17d 17e 17f 17g n = 2; 7 membered-rings 27a 27b 27c 27d 27e 27f 27g a
enzymatic BACE1 IC50 (nM)
cellular BACE1 IC50 (nM)
0 0 0 0 0 0 0
1950 575 813 759 676 3390 25.7
± ± ± ± ± ± ±
826 83.1 113 305 177 167 4.19
1700 417 708 589 588 2190 25.1
± ± ± ± ± ± ±
308 54.4 46.1 76.8 19.2 470 4.10
1 1 1 1 1 1 1
9.55 3.24 2.45 4.47 1.74 2.95 1.05
± ± ± ± ± ± ±
1.56 0.27 1.90 1.56 0.09 0.44 0.38
6.61 0.45 0.91 0.93 3.24 8.71 0.29
± ± ± ± ± ± ±
2.84 0.31 0.08 0.18 1.95 2.57 0.26
2 2 2 2 2 2 2
195 49.0 45.7 158 138 214 7.08
± ± ± ± ± ± ±
38.2 12.8 6.79 57.4 20.5 8.19 1.99
437 55.0 77.6 117 309 372 9.55
± ± ± ± ± ± ±
18.0 2.92 6.20 9.68 70.7 67.5 3.01
All data are mean of at least three repeats, with standard deviations provided.
Table 2. MUE and Standard Deviations for the Different FEP, Docking, and MM-GBSA Runs versus Experiment, Performed on the Scaffold Modifications
MUEa (kcal/mol) standard deviation (kcal/mol) correct ranking
XP Glide Docking
MM-GBSA 1b
MM-GBSA 2c
2 ns (3 repeats)
1.9 ± 0.64
5.3 ± 1.4
4.4 ± 1.4
0.81 ± 0.31 1.2 7/7
5 ns (3 repeats)
10 ns (3 repeats)
40 ns (3 repeats)
5 ns (3 repeats)d
1.3 ± 0.53
0.82 ± 0.30 0.97
0.91 ± 0.41 1.1
0.92 ± 0.39 1.1
0.57 ± 0.11 0.29
3/7
5/7
6/7
3/7
6/7
5 ns 1 sim
MUE is reported with ±90% confidence interval. MM-GBSA performed with only ligand minimization in a static receptor. MM-GBSA performed with ligand and an 8 Å radius of protein in the minimization. dPerformed with “Core Hopping FEP+”. a
b
c
ring prototype containing the methoxypyrazine tail d (molecule 99) had BACE1 enzymatic inhibition activity of 813 nM (not shown in Table 1) which compared to 0.93 nM for the equivalent 6 membered-ring analogue (17d). Given the relatively low activity we did not further consider this larger scaffold. FEP Scaffold Perturbations. FEP calculations were performed with sets of molecules containing different 5, 6, and 7 membered-ring scaffolds but the same R-group, Table 2 and Figure 3. The FEP Mapper comprised three perturbations replacing the entire cyclic amidine scaffold, a substantial structural change for calculations of this type. The calculations were run in triplicate for 2, 5, 10, and 40 ns lambda window simulation time, and the results are presented as averages with standard deviations. The calculated relative binding energies
were compared to experiment with all values normalized with respect to the least active, 5 membered-ring scaffold. Overall, the calculations suggested the 6 membered-ring scaffold would be most potent with the 5 membered-rings the least. The 2 ns simulations correctly predicted the 6 membered scaffolds to be the most active followed by the 7 and then 5 membered-rings (Figure 3, Table S1). This was the case for all the seven different P3 substituents. However, the standard deviations for 27a−27c cross zero making the direction of relative binding energy predictions less reliable. The 5 ns simulations successfully predicted the rank order of scaffold potency for 5 of the 7 different P3 R-groups. The order of activity for the a and g R-groups were incorrectly ranked and the standard deviations crossed zero for relative binding energy predictions of 27a, 27c−27e, and 27g. For the 10 ns 1445
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation
FEP P3 Pocket Perturbations. Calculations were performed perturbing the 7 different R-group substituents on the same scaffold. The simulations were run for 2, 5, 10, 20, and 40 ns, for all three scaffolds, and were repeated three times. Results for each scaffold are shown together in Figure 5 and separately in Figures S1−S3. The calculated data for each compound with the 7 membered-ring scaffold are provided in Table 3, along with the analogous data for the corresponding 5 and 6 membered-rings in Tables S2 and S3, respectively. Statistics for the FEP performance on all three scaffolds are summarized in Table 4. The performance of FEP on the 7 membered-ring scaffold shown in Figure 5 and Table 3 shows an improvement for both repeated simulations and increasing simulation time. The standard 5 ns single run showed an MUE with experiment of 0.68 ± 0.38 kcal/mol, Table 4. Averaging over three repeats delivered an MUE of 0.48 ± 0.37 kcal/mol and an average standard deviation of 0.40 kcal/mol while the standard deviation worsened for some molecules such as 27c, Table 3. The performance at 10 ns simulation was MUE of 0.39 ± 0.16 kcal/mol, but 27c still showed very large standard deviation between repeat simulations. Interestingly, this is resolved with the longer 40 ns simulations that showed significantly better MUE with experiment (0.09 ± 0.05 kcal/mol), suggesting further sampling is required for the complex to find an optimal binding pose. In general, the standard deviation also improved with longer simulation times, suggesting conformational sampling may be a source of error. One molecule, 27a, still showed a large standard deviation (0.87 kcal/mol, Table 3) in the 40 ns simulations. The same trends tend to hold for the 5 and 6 membered-ring scaffolds. For the 5 membered-ring the MUE reduced slightly from 0.48 ± 0.33, 0.43 ± 0.18, 0.41 ± 0.18 to 0.36 ± 0.15 kcal/ mol (Figure 5 and Table S2). Notably, the standard deviation was much smaller, 0.19, 0.14, 0.06 and 0.13 for 2, 10, 20, and 40 ns respectively. Hence, comparing with the 7 membered-ring, the scaffold affects the calculations despite being distal from the perturbation site. Performance for the 6 membered-ring was good and comparable for all the time length simulations (Figure 5 and Table S3). The R2 has little significance for such a small spread in experimental data (just 1.24 kcal/mol) but the MUE was stable at 0.61 ± 0.35, 0.53 ± 0.24, 0.60 ± 0.25 and 0.60 ± 0.28 kcal/mol for 2, 10, 20, and 40 ns FEP simulations. Also, the standard deviation was small for all simulation lengths, 0.14, 0.18, 0.15, and 0.09 for 2, 10, 20, and 40 ns. Hence, increasing simulation time did not have a negative impact on quality and the calculations were consistent. As was the case for the scaffold perturbations, a single 5 ns FEP simulation produced worse results, MUE of 0.69 ± 0.29, 0.67 ± 026, and 0.68 ± 0.38 kcal/mol for 5, 6, and 7 membered-rings respectively. Comparison was again performed with Glide XP docking, and the binding energy prediction from two different MM-GBSA protocols. In the former case the MUE for 5, 6 and 7 membered-ring sets of molecules was 1.6 ± 1.0, 1.4 ± 0.69 and 1.7 ± 1.0 kcal/mol respectively. The second of the two MM-GBSA protocols performed worse delivering MUE values of 3.9 ± 2.9, 4.3 ± 2.6, and 3.1 ± 2.6 for the increasing scaffolds from 5 to 7 membered-rings. Hence FEP performed better than docking and MM-GBSA. In summary, the FEP results for the P3 R-group perturbations showed good performance with low MUE. The quality of results for the 7 membered-ring and to a lesser extent the 5 membered-ring benefited from increasing simulation time, whereas the 6 membered-ring did not.
simulations, the relative binding energies of the R-groups for 6 of the 7 scaffolds were successfully predicted. The g R-group was again incorrectly predicted, with 8g incorrectly predicted to be more active than 27g. Predicted binding energies of the a Rgroup agreed with experiment but the standard deviation of 27a crossed zero. For the 40 ns simulations only 3 of the 7 different R-group rankings were correctly predicted, the a, c, e, and g Rgroups were incorrect and the standard deviations crossed zero for all R-groups except d. The mean unsigned errors (MUE) between predicted and experimental binding energies improved substantially for the average of three repeat 5 ns simulations compared to one single 5 ns simulation (0.82 ± 0.30 compared to 1.30 ± 0.53 kcal/ mol), Table 2. The longer simulation times of 10 and 40 ns did not lead to any improvement in MUE compared to 2 or 5 ns simulations (0.91 ± 0.41 and 0.92 ± 0.39 compared to 0.81 ± 0.31 and 0.82 ± 0.30 kcal/mol respectively). Hence, three repeat shorter simulations can provide better average results than a single longer simulation. In contrast to what we may expect, the standard deviation also does not improve for longer simulation time. Indeed, the predicted binding energy from the 10 and 40 ns simulations for specific molecules such as 27f and 27g show an especially large standard deviation, Figure 3 and Table S1. This suggests there are differences between repeat simulations which are not resolved via sampling to this time limit. Comparing with other methods, Glide XP docking and MM-GBSA protocols delivered MUE’s of 1.9 ± 0.64, 5.3 ± 1.4, and 4.4 ± 1.4 respectively, suggesting worse performance than FEP. Overall, the FEP calculations showed a trend for preference of the 6 membered-ring above the 7 and, subsequently, 5 membered-ring scaffolds. The recently developed “Core Hopping FEP+” was used to understand if the calculations on scaffold perturbations could be further improved. The same Mapper setup of 7 independent simulations of three molecules (Figure 3 panel h) was used along with a 5 ns lambda window simulation time and three repeats. The results of the calculations were improved and correlated well with experimental data, Figure 4. Particularly the
Figure 4. “Core Hopping FEP+” results for the scaffold perturbations (legend in Figure 3, panel h). Simulations were performed for 5 ns and repeated three times. The binding energies are normalized with respect to the least active of the three molecules (always the 5 memberedring) with experimental data shown in dark and light gray, and calculated binding data shown in dark and light blue. The error bars show the standard deviation.
stability of the calculations improved greatly, with an average standard deviation of 0.29 kcal/mol compared to 0.97 kcal/mol for the same simulation time with standard FEP. The average MUE also improved to 0.57 ± 0.11 kcal/mol, Table 2. With the “Core Hopping FEP+” 6 out of 7 different P3 R-groups were ranked correctly, where only 27e was ranked less potent than 8e. 1446
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation
Figure 5. Correlation of FEP calculated binding energies and experiment for seven different R-group substituents all on the same 5 (green spheres), 6 (blue diamonds) and 7 (red squares) membered-ring scaffolds. Calculations were performed for 2, 5, 10, 20, and 40 ns and were repeated three times, with standard deviations in both calculations and experiment shown with error bars. Black solid line corresponds to x = y, and black dashed lines corresponds to x = y ± 1. Correlation between scaffolds should not be considered in this figure as different offset terms are used to be able to plot the data in the same graphs. For separate plots, see Figures S1−S3.
Table 3. Details for FEP Calculations on 7 Membered-Ring Scaffold with Different P3 R-Group Perturbationsa
a
Calculated data are averages of three repeat simulations with standard deviations shown in parenthesis.
1447
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation
Table 4. Performance of FEP, Docking, and MM-GBSA Predictions versus Experiment for P3 R-Group Perturbations on Each Scaffolda XP Glide Docking
MM-GBSA 1c
MM-GBSA 2d
2 ns (3 repeats)
5 ns 1 simulation
5 ns (3 repeats)
10 ns (3 repeats)
20 ns (3 repeats)
40 ns (3 repeats)
R2 MUEb (kcal/mol) standard deviation (kcal/mol)
0.38 1.6 ± 1.0
0.08 4.8 ± 2.3
5 membered-ring scaffold, molecules 8a−8g 0.03 0.38 0.16 0.49 3.9 ± 2.9 0.48 ± 0.33 0.69 ± 0.29 0.57 ± 0.23 0.19 0.17
0.66 0.43 ± 0.18 0.14
0.68 0.41 ± 0.18 0.05
0.79 0.36 ± 0.15 0.13
R2 MUE (kcal/mol) standard deviation (kcal/mol)
0.21 1.4 ± 0.69
0.25 3.9 ± 2.5
6 membered-ring scaffold, molecules 17a−17g 0.08 0.45 0.14 0.27 4.3 ± 2.6 0.61 ± 0.35 0.67 ± 0.26 0.63 ± 0.27 0.14 0.19
0.27 0.53 ± 0.24 0.18
0.19 0.60 ± 0.25 0.15
0.18 0.60 ± 0.28 0.09
R2 MUE (kcal/mol) standard deviation (kcal/mol)
0.24 1.7 ± 1.0
0.43 6.2 ± 2.8
7 membered-ring scaffold, molecules 27a−27g 0.07 0.58 0.49 0.57 3.1 ± 2.6 0.58 ± 0.26 0.68 ± 0.38 0.48 ± 0.37 0.37 0.40
0.82 0.39 ± 0.16 0.53
0.74 0.37 ± 0.22 0.39
0.98 0.09 ± 0.05 0.36
MUE (kcal/mol)
1.6 ± 0.51
5.0 ± 1.4
0.45 ± 0.11
0.46 ± 0.13
0.35 ± 0.13
R2 MUE (kcal/mol) standard deviation (kcal/mol)
0.00 0.57 ± 0.20
0.06 2.4 ± 0.77
0.04 0.36 ± 0.15 0.06
0.00 0.37 ± 0.15 0.03
0.00 0.38 ± 0.16 0.07
averages over three scaffolds 3.8 ± 1.5 0.56 ± 0.18 0.68 ± 0.17 0.56 ± 0.16 5 membered-ring scaffold, molecules 28−39 0.02 0.08 0.05 0.05 7.8 ± 2.0 0.38 ± 0.14 0.38 ± 0.15 0.37 ± 0.15 0.06 0.07
a 2 R correlation coefficient, MUE, and standard deviations are shown for the different simulation times and averaged over three repeats. bMUE is reported with ± 90% confidence interval. cMM-GBSA performed with only ligand minimization in a static receptor. dMM-GBSA performed with ligand and an 8 Å radius of protein in the minimization.
Figure 6. Correlation of FEP calculated binding energies and experiment for twelve different P2′ R-group substituents all on the same 5 memberedring scaffold. Calculations were performed for 2, 5, 10, 20, and 40 ns and were repeated three times, with standard deviations in both calculations and experiment shown with error bars.
FEP P2′ Pocket Perturbations. Calculations were performed on the 12 molecules shown in Figure 6 (and Table S4) which shared the same 5 membered-ring scaffold but different P2′ pocket substituents. Modeling the initial binding
mode to launch the FEP calculations revealed that some compounds can adopt two alternative binding modes. They were distinguished by performing FEP calculations to find the lowest energy conformation (Table S5). Simulations were again 1448
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation
Figure 7. Understanding large standard deviations in the P3 R-group simulations on the 7 membered-ring scaffold. (A) The starting structure of the active site flap (gray) can be in an open (cyan) or closed (magenta) conformation . The measured distance between Cα of Asp93 (also referred to as Asp32) and Cα of Gln134 is indicated with a dashed line. (B) The CαAsp93 CαGln134 distance as a function of time during 10 ns simulations of 27c. During the simulation the flap either opens or closes depending on the repeat. (C) In the 40 ns simulations the opening and closure of the flap is sampled several times.
(methoxyethoxy)pyrazine f were typically least active and predicted as such in most FEP protocols. The former does not fully occupy the P3 pocket whereas the latter fills the pocket but likely with a higher entropic cost. The g R-group was usually the most potent as it can fill the pocket without a conformational penalty because of the rigid nature of the alkyne. For the scaffold perturbations, the sets of just three molecules did not allow comparison of correlations. However, the MUE were in the range of 0.81−0.92 kcal/mol when averaged over three repeats. No improvement in MUE or in the standard deviation for the three repeat simulations was seen for the simulation times attempted. Such very large perturbations (obliteration and reconstruction of a flexible ring) may require even longer simulations or alternative coupling protocols. To understand if these effects can be overcome we have studied the scaffold perturbation with a new “Core Hopping FEP+” implementation designed to specifically tackle and improve results for such cases.45 Medicinal chemistry explorations frequently consider alternative heterocyclic scaffolds with similar shape and electrostatic properties, or as here, homologation with additional CH2 groups to increase ring size. Such approaches represent a challenge for traditional FEP, because despite being small structural changes, they often require removal and reconstruction of a crucial motif, sometimes in the center of a lead molecule. This is a longknown problem in traditional single-topology FEP, because the energy contribution by decoupled “dummy” atoms does not automatically cancel if they are connected to the rest of the system by more than one bond. This makes transformations like ring-breaking and linker insertion hard to accomplish. Hence, the recent development of “Core Hopping FEP+”, a new approach designed to avoid the sampling problems caused
run for 2, 5, 10, 20, and 40 ns and in each case repeated three times (Figure 6 and Table S4). The performance of the calculations was not affected by simulation time or repeats, consistently giving very low MUE of 0.38 ± 0.14, 0.37 ± 0.15, 0.36 ± 0.15, 0.37 ± 0.15, 0.38 ± 0.16 kcal/mol and standard deviation of 0.1, 0.1, 0.1, 0.0, 0.1 for 2, 5, 10, 20, and 40 ns. This is in contrast to perturbations in the P3 pocket where molecules from the same 5 membered-ring scaffold showed improvements in both repeats and simulation time. Indeed, perturbations in the P2′ pocket show very good results with the default 5 ns simulation. This highlights that different parts of the same protein benefit differently from increased sampling. Also for this set of compounds comparison was made to Glide XP docking and MM-GBSA. For Glide XP docking MUE was 0.57 ± 0.20 kcal/mol, in a similar range as FEP. The two MMGBSA protocols performed worse, the MUE was 2.44 ± 0.77 and 7.8 ± 2.0 kcal/mol for protocol 1 and 2. In summary, the FEP calculations for substituents of R-groups in P2′ pocket performed very well with low MUE’s. Good performance was seen already with default settings, and no benefit of increased simulation time was seen.
■
DISCUSSION AND CONCLUSIONS Among the synthesized molecules the 6 membered-ring examples were the most potent, 17g being 1 nM in the enzymatic assay, the 7 membered-ring showed intermediate activity and the 5 membered-ring was least potent. The standard FEP calculations on the scaffold perturbations generally predicted the 6 membered-ring to be the most potent, however, the standard deviations in repeats were large. The synthetic exploration introduced the same P3 R-groups into each of the three different scaffolds, enabling the systematic computational study. The oxazole a or the 21449
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation
Figure 8. Movement of the 10s loop during FEP simulations: (A) Starting structure of 27a and 10s loop in gray, and open conformation of the 10s loop is shown in orange and closed conformation in blue. The measured distance between Cα of Ser71 and Cα of Gly291 is shown as a dashed line. (B) The distance of the 10s loop to the bottom of the pocket as a function of time. The movement of the 10s loop in simulations of 40 ns on outlier 27a (dark markers) are compared with 27g (light blue), which was not an outlier. Only one simulation of 27a goes to the closed conformation of the 10s loop, while the closed conformation is sampled substantially for 27g.
by such modifications. This is achieved by introducing a modified covalent bond potential energy term which can be rigorously introduced and removed in free energy transformations; see ref 45. For details on this soft bond stretch potential. Note that the properties of the unperturbed systems (lambda values 0 and 1) remain unchanged by this force field modification. A validation study45 of six pharmaceutically relevant test sets showed ‘Core Hopping FEP+’ could predict relative binding free energies with a small root-mean-squared error of 0.5 kcal/mol for a variety of scaffold modifications. In this case we used the same Mapper setup of 7 independent simulations of three molecules and the MUE and standard deviation for a 5 ns simulation repeated three times was found to be 0.57 ± 0.11 and 0.29 kcal/mol. These results show the beneficial application of FEP for scaffold hopping which can be considered outside the traditional lead optimization focus. For the P3 pocket R-groups FEP performance was consistent with another recent BACE1 study from our group,20 where we observed an MUE of 0.57 ± 0.11 and 0.59 ± 0.29 kcal/mol for 20 ns FEP simulations on two alternative BACE1 inhibitor chemical series. MUE performance improved with increasing simulation time for 5 and 7 membered-rings. Also, the standard deviation was larger for the 7 membered-rings, which improved with longer simulation times but was low for the 6 memberedring scaffold. This suggests the 7 membered-ring simulations start in either a suboptimal ligand or protein conformation that benefit from increased sampling. Hence, extra simulation time can be needed not only to overcome insufficient sampling in the vicinity of the perturbation but also distal to the region of modification. It is reassuring that by proper sampling of the protein and ligand movements MD/FEP can give excellent
correlation to experiment, even when performing the calculations on flexible systems and sets of ligands beyond the cocrystalized ligand space. To understand the change in performance with increased simulation time for the P3 pocket perturbations on the 7 membered-ring scaffold the FEP/MD trajectories were analyzed in more detail. Molecule 27c had the largest standard deviation from the repeat 10 ns simulations. For one repeat the active site flap (comprised of amino acids 66 to 77) went from the starting structure to a closed conformation and finished open, whereas for another the loop made the opposite movement (Figure 7). As such these two repeats gave very different binding energies, while the average of three was correct. With longer simulations the flap movement could be sampled correctly several times, suggesting why the standard deviation for these 7 membered-ring molecules improved at 20 and 40 ns. Second, we considered 27a which was the only substituent with a large standard deviation in the longer, 40 ns simulations. In this case we noticed that in only a few of the simulations the 10s loop, at the base of the P3 pocket (residues 9−14),46 would briefly close, while in most the loop remained open. Comparing to 27g, which did not have a large standard deviation and had good agreement between predicted and experimental binding energy, revealed that the 10s loop would open and close throughout the simulations, Figure 8. Hence, it seems that the 10s loop conformations may become trapped open in the case of 27a and not be fully sampled, whereas for other molecules such as 27g this was not the case and the loop showed a much wider freedom of movement. Compounds with the same 5 membered-ring scaffold but varying substituents going into the P2′ pocket were also 1450
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation
transformation between two molecules with 40 ns lambda window takes roughly 40 h using 4 processors. Combining this with a set of transformations between many compounds, repeated simulations, and limited number of GPU processors these calculations still takes many weeks to finish. In this study, FEP calculations were applied to 21 synthesized acylguanidine BACE1 inhibitors with modified scaffold size and P3 pocket substituents. The molecules had good potency in enzymatic BACE1 inhibition assays, achieving nanomolar activity. We tested the performance of the FEP approach to predict relative binding energies for structural modifications in different regions of the binding site. In all cases three repeat calculations were performed. Naturally this greatly increases the simulation time. However, it leads to improved quality of results for difficult cases, such as the standard approach for scaffold perturbations and the P3 pocket R-groups where protein and hence ligand flexibility are relevant. In contrast, for the P2′ pocket repeats and longer simulation time offered little benefit, but did not lead to worsening of results. Meanwhile, the scaffold perturbations did not improve with longer simulation time, suggesting even up to 40 ns energies were still not converged due to the difficultly of the perturbation. In some cases three repeats of 5 ns simulations can perform better than longer simulations of 40 ns, a useful conclusion given the increased opportunities for parallelization of repeats versus single trajectories. Overall, calculated binding energy changes associated with modifying the scaffold size gave large differences between repeated calculations, and the correlation to experiment was not satisfying, although this improved with a new “Core Hopping FEP+” approach. However, when modifying the P3 and P2′ pocket substituents, good agreement with experiment was seen. The correlation to experiment was further improved by repeated simulations with different initial velocities, and longer simulation times and is well within the error range useful for quantitative molecular design.
studied. The FEP calculations on these compounds showed very low MUEs with the default settings (single simulation, 5 ns), and no benefit on increased simulation time was seen. Additional sampling was not beneficial because the ligands 28 to 39 started the FEP calculations in accurate initial binding poses based on the crystallized ligand (PDB 3IN4) which is structurally very similar containing a 2,6 diethyl pyridine in place of the R-substituted phenyl (Figure 6). Overall, predicting binding energies in this pocket was a more achievable task because even the docking approach had relatively good performance and low MUE. Hence, although showing some better performance the benefit of the MD/FEP method was less. Many of the FEP strategies, particularly longer simulation times, outperformed Glide and MM-GBSA by showing better MUE. Despite the scaffold perturbations being larger than the R-group modifications, they also had better correlation to experiment than the alternative methods. Worse performance of Glide and MM-GBSA compared to FEP has now been seen in multiple reports including fragment sized compounds in eight different systems,47 and multiple drug discovery projects including BACE116 which had an MUE of 0.84 kcal/mol for 36 molecules, subsequently improved further with the application of the OPLSv3 force field.36 Details of FEP applications are emerging but there is a need to examine difficult cases and develop strategies to improve results (such as what we have shown for using repeat short simulations). A detailed study of 21 c-Jun N-terminal kinase-1 inhibitors showed that accounting for alternative dihedral orientations of input ligands and improving force field charges could improve the MUE to 0.9 kcal/mol.48 When protein reorganization is required, simulations need to be longer than 5 ns simulation and the replica exchange with solute tempering (REST) region should be extended to include protein residues.49 These applications show the potential of FEP but also highlight the difficulties which can be encountered. A basic but often neglected requirement in statistical thermodynamics is proper sampling.50 By the simple approach of repeating the calculations, meaningful statistics can be obtained, and the reproducibility of the calculations can be quantified by standard deviations. Our results showed that in some difficult cases variation between repeats can be very large, for instance applying standard FEP to scaffold perturbations. In addition, as we stressed in our previous report,20 relatively small data sets mean that confidence intervals often overlap and caution is needed with over interpretation of differences in results as robust statistical significance is often lacking. In general, averaging over three repeats delivered improvements with respect to experiment. Also longer simulations than today’s standard, of 2 or 5 ns, improved the statistics. By combining these two approaches proper error and trajectory analysis can be made. It also allows suboptimal starting structures of the ligands or protein to have a fair chance of finding the optimal pose. Even though the core of the molecule within a series is consistently placed, small substituents might have a large impact on the overall binding pose. Furthermore, flexibility and movement of ligand-protein complexes are a natural part of their behavior, such the 10s loop in the P3 pocket, or the flap in the scaffold region. Sampling these movements is vital to correctly calculate high quality binding energies useful for molecular design in lead optimization projects. Admittedly, these calculations do not come without a cost. On a cluster with Tesla K80 graphics cards, each full
■
ASSOCIATED CONTENT
* Supporting Information S
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b01141. P3 pocket R-group substituent, P3 pocket FEP results for 5 membered-ring scaffold, 6 membered-ring scaffold, and 7 membered-ring scaffold, detailed FEP results for the scaffold perturbations and for 5 and 6 membered-ring scaffolds with different P3 R-group perturbations, connection map (mapper) used for the P3 pocket substituent FEP calculations, detailed FEP results on 5 membered-ring scaffold with different P2′ R-group perturbations, initial FEP calculations on P2′ substituents to decide starting pose, connection map (mapper) used for the P2′ pocket substituent FEP calculations, experimental and computational protocols, general methods, and synthetic protocols (PDF) NMR analytical data (PDF) Ligand coordinates in separate sd files for each of the datasets (ZIP)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Tel.: +34 925 24 5777. ORCID
Andrés A. Trabanco: 0000-0002-4225-758X 1451
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Article
Journal of Chemical Theory and Computation
A.; Papaioannou, N.; Richard, D.; Ryan, M. S.; Wan, Z.-K.; Thorarensen, A. Imidazotriazines: Spleen Tyrosine kinase (Syk) Inhibitors Identified by Free Energy Perturbation (FEP). ChemMedChem 2016, 11, 217−233. (18) De Vivo, M.; Masetti, M.; Bottegoni, G.; Cavalli, A. Role of Molecular Dynamics and Related Methods in Drug Discovery. J. Med. Chem. 2016, 59, 4035−4061. (19) Rombouts, F. J. R.; Tresadern, G.; Buijnsters, P.; Langlois, X.; Tovar, F.; Steinbrecher, T.; Vanhoof, G.; Somers, M.; Andrés, J.-I.; Trabanco, A. A. Pyrido[4,3-e][1,2,4]triazolo[4,3-a]pyrazines as Selective, Brain Penetrant Phosphodiesterase 2 (PDE2) Inhibitors. ACS Med. Chem. Lett. 2015, 6, 282−286. (20) Ciordia, M.; Pérez-Benito, L.; Delgado, F.; Trabanco, A. A.; Tresadern, G. Application of Free Energy Perturbation for the Design of BACE1 Inhibitors. J. Chem. Inf. Model. 2016, 56, 1856−1871. (21) Hardy, J. A.; Higgins, G. A. Alzheimer’s Disease: The Amyloid Cascade Hypothesis. Science 1992, 256, 184−185. (22) Oehlrich, D.; Prokopcova, H.; Gijsen, H. J. M. The Evolution of Amidine-Based Brain Penetrant BACE1 Inhibitors. Bioorg. Med. Chem. Lett. 2014, 24, 2033−2045. (23) Malamas, M. S.; Erdei, J.; Gunawan, I.; Barnes, K.; Johnson, M.; Hui, Y.; Turner, J.; Hu, Y.; Wagner, E.; Fan, K.; Olland, A.; Bard, J.; Robichaud, A. J. Aminoimidazoles as Potent and Selective Human βSecretase (BACE1) Inhibitors. J. Med. Chem. 2009, 52, 6314−6323. (24) Malamas, M. S.; Erdei, J.; Gunawan, I.; Turner, J.; Hu, Y.; Wagner, E.; Fan, K.; Chopra, R.; Olland, A.; Bard, J.; Jacobsen, S.; Magolda, R. L.; Pangalos, M.; Robichaud, A. J. Design and Synthesis of 5,5′-Disubstituted Aminohydantoins as Potent and Selective Human β-Secretase (BACE1) Inhibitors. J. Med. Chem. 2010, 53, 1146−1158. (25) Godyn, J.; Jonczyk, J.; Panek, D.; Malawska, B. Therapeutic Strategies for Alzheimer’s Disease in Clinical Trials. Pharmacol. Rep. 2016, 68, 127−138. (26) Shoichet, B. K.; Walters, W. P.; Jiang, H.; Bajorath, J. Advances in Computational Medicinal Chemistry: A Reflection on the Evolution of the Field and Perspective Going Forward. J. Med. Chem. 2016, 59, 4033−4034. (27) Woltering, T. J.; Wostl, W.; Hilpert, H.; Roger-Evans, M.; Pinard, E.; Mayweg, A.; Göbel, M.; Banner, D. W.; Benz, J.; Travagli, M.; Pollastrini, M.; Marconi, G.; Gabellieri, E.; Guba, W.; Mauser, H.; Andreini, M.; Jacobsen, H.; Power, E.; Narquizian, R. BACE1 Inhibitors: A Head Group Scan on a Series of Amides. Bioorg. Med. Chem. Lett. 2013, 23, 4239−4243. (28) Tresadern, G.; Delgado, F.; Delgado, O.; Gijsen, H.; Macdonald, G. J.; Moechars, D.; Rombouts, F.; Alexander, R.; Spurlino, J.; Van Gool, M.; Vega, J. A.; Trabanco, A. A. Rational Design and Synthesis of Aminopiperazinones as Beta-Secretase (BACE) Inhibitors. Bioorg. Med. Chem. Lett. 2011, 21, 7255−7260. (29) Ginman, T.; Viklund, J.; Malmström, J.; Blid, J.; Emond, R.; Forsblom, R.; Johansson, A.; Kers, A.; Lake, F.; Sehgelmeble, F.; Sterky, K. J.; Bergh, M.; Lindgren, A.; Johansson, P.; Jeppsson, F.; Fälting, J.; Gravenfors, Y.; Rahm, F. Core Refinement Towards Permeable β-Secretase (BACE-1) Inhibitors with Low hERG Activity. J. Med. Chem. 2013, 56, 4181−4205. (30) Rombouts, F. J. R.; Tresadern, G.; Delgado, O.; MartinezLamenca, C.; Van Gool, M.; Garcia-Molina, A.; Alonso de Diego, S. A.; Oehlrich, D.; Prokopcova, H.; Alonso, J. M.; Austin, N.; Borghys, H.; Van Brandt, S.; Surkyn, M.; De Cleyn, M.; Vos, A.; Alexander, R.; Macdonald, G.; Moechars, D.; Gijsen, H.; Trabanco, A. A. 1,4-Oxazine β-Secretase 1 (BACE1) Inhibitors: From Hit Generation to Orally Bioavailable Brain Penetrant Leads. J. Med. Chem. 2015, 58, 8216− 8235. (31) Banner, D. W.; Gsell, B.; Benz, J.; Bertschinger, J.; Burger, D.; Brack, S.; Cuppuleri, S.; Debulpaep, M.; Gast, A.; Grabulovski, D.; Hennig, M.; Hilpert, H.; Huber, W.; Kuglstatter, A.; Kusznir, E.; Laeremans, T.; Matile, H.; Miscenic, C.; Rufer, A.; Schlatter, D.; Steyaert, J.; Stihle, M.; Thoma, R.; Weber, M.; Ruf, A. Mapping the Conformational Space Accessible to Bace2 Using Surface Mutants and Co-Crystals with Fab-Fragments, Fynomers, and Xaperones. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2013, 69, 1124−1137.
Gary Tresadern: 0000-0002-4801-1644 Author Contributions
H.K. and L.P.-B. contributed equally. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We thank the wider Janssen BACE1 research team, as well as members of screening and analysis and purification departments.
■
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, 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) Bash, P.; Singh, U.; Langridge, R.; Kollman, P. Free Energy Calculations by Computer Simulation. Science 1987, 236, 564−568. (6) Kollman, P. Free Energy Calculations: Applications to Chemical and Biochemical Phenomena. Chem. Rev. 1993, 93, 2395−2417. (7) Wong, C. F.; McCammon, J. A. Dynamics and Design of Enzymes and Inhibitors. J. Am. Chem. Soc. 1986, 108, 3830−3832. (8) 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. (9) Deng, Y.; Roux, B. Computations of Standard Binding Free Energies with Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113, 2234−2246. (10) Durrant, J.; McCammon, J. A. Molecular Dynamics Simulations and Drug Discovery. BMC Biol. 2011, 9, 71−79. (11) Michel, J.; Essex, J. W. Hit Identification and Binding Mode Predictions by Rigorous Free Energy Simulations. J. Med. Chem. 2008, 51, 6654−6664. (12) Bollini, M.; Domaoal, R. A.; Thakur, V. V.; Gallardo-Macias, R.; Spasov, K. A.; Anderson, K. S.; Jorgensen, W. L. Computationally Guided Optimization of a Docking Hit to Yield Catechol Diethers as Potent Anti-HIV Agents. J. Med. Chem. 2011, 54, 8582−8591. (13) Christ, C. D.; Fox, T. Accuracy Assessment and Automation of Free Energy Calculations for Drug Design. J. Chem. Inf. Model. 2014, 54, 108−120. (14) Homeyer, N.; Stoll, F.; Hillisch, A.; Gohlke, H. Binding Free Energy Calculations for Lead Optimization: Assessment of Their Accuracy in an Industrial Drug Design Context. J. Chem. Theory Comput. 2014, 10, 3331−3344. (15) Mikulskis, P.; Genheden, S.; Ryde, U. A Large-Scale Test of Free-Energy Simulation Estimates of Protein-Ligand Binding Affinities. J. Chem. Inf. Model. 2014, 54, 2794−2806. (16) Wang, L.; Wu, Y. J.; Deng, Y. Q.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; Romero, D. L.; Masse, C.; Knight, J. L.; Steinbrecher, T.; Beuming, T.; Damm, W.; Harder, E.; Sherman, W.; Brewer, M.; Wester, R.; Murcko, M.; Frye, L.; Farid, R.; Lin, T.; Mobley, D. L.; Jorgensen, W. L.; Berne, B. J.; Friesner, R. A.; Abel, R. Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug Discovery by Way of a Modern Free-Energy Calculation Protocol and Force Field. J. Am. Chem. Soc. 2015, 137, 2695−2703. (17) Lovering, F.; Aevazelis, C.; Chang, J.; Dehnhardt, C.; Fitz, L.; Han, S.; Janz, K.; Lee, J.; Kaila, N.; McDonald, J.; Moore, W.; Moretto, 1452
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453
Journal of Chemical Theory and Computation
■
(32) Malamas, M. S.; Barnes, K.; Johnson, M.; Hui, Y.; Zhou, P.; Turner, J.; Hu, Y.; Wagner, E.; Fan, K.; Chopra, R.; Olland, A.; Bard, J.; Pangalos, M.; Reinhart, P.; Robichaud, A. J. Di-substituted Pyridinyl Aminohydantoins as Potent and Highly Selective Human BetaSecretase (BACE1) Inhibitors. Bioorg. Med. Chem. 2010, 18, 630−639. (33) Schrodinger LLC. Maestro, 2014. http://www.schrodinger. com/Maestro/. (34) Madhavi Sastry, G.; 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. (35) Watts, K. S.; Dalal, P.; Murphy, R. B.; Sherman, W.; Friesner, R. A.; Shelley, J. C. ConfGen: A Conformational Search Method for Efficient Generation of Bioactive Conformers. J. Chem. Inf. Model. 2010, 50, 534−546. (36) 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. (37) 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., In Proceedings of the 2006 ACM/IEEE Conference on Supercomputing; ACM: Tampa, FL, 2006; p 84. (38) 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. (39) 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. (40) 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. (41) Liu, S.; Wu, Y.; Lin, T.; Abel, R.; Redmann, J.; Summa, C.; Jaber, V.; Lim, N.; Mobley, D. Lead Optimization Mapper: Automating Free Energy Calculations for Lead Optimization. J. Comput.-Aided Mol. Des. 2013, 27, 755−770. (42) Tang, T. P.; Ellman, J. A. The tert-Butanesulfinyl Group: An Ideal Chiral Directing Group and Boc-Surrogate for the Asymmetric Synthesis and Applications of β-Amino Acids. J. Org. Chem. 1999, 64, 12−13. (43) Banner, D., Hilpert, H.; Humm, R.; Mauser, H.; Mayweg, A. V.; Ricklin, F., Rogers-Evans, M. Dihydropyrimidinones for use as BACE2 Inhibitors. WO-2010/128058, May 5, 2010. (44) For biological assays details, see Supporting Information. (45) Wang, L.; Deng, Y.; Wu, Y.; Kim, B.; LeBard, D. N.; Wandschneider, D.; Beachy, M.; Friesner, R. A.; Abel, R. Accurate Modeling of Scaffold Hopping Transformations in Drug Discovery. J. Chem. Theory Comput. 2017, 13, 42−54. (46) Patel, S.; Vuillard, L.; Cleasby, A.; Murray, C. W.; Yon, J. Apo and Inhibitor Complex Structures of BACE (β-Secretase). J. Mol. Biol. 2004, 343, 407−416. (47) Steinbrecher, T. B.; Dahlgren, M.; Cappel, D.; Lin, T.; Wang, L.; Krilov, G.; Abel, R.; Friesner, R.; Sherman, W. Accurate Binding Free Energy Predictions in Fragment Optimization. J. Chem. Inf. Model. 2015, 55, 2411−2420. (48) Kaus, J. W.; Harder, E.; Lin, T.; Abel, R.; McCammon, J. A.; Wang, L. How to Deal with Multiple Binding Poses in Alchemical Relative Protein-Ligand Binding Free Energy Calculations. J. Chem. Theory Comput. 2015, 11, 2670−2679. (49) Lim, N. M.; Wang, L.; Abel, R.; Mobley, D. L. Sensitivity in Binding Free Energies Due to Protein Reorganization. J. Chem. Theory Comput. 2016, 12, 4620−4631. (50) Coveney, P. D.; Wan, S. On the Calculation of Equilibrium Thermodynamic Properties from Molecular Dynamics. Phys. Chem. Chem. Phys. 2016, 18, 30236−30240.
Article
NOTE ADDED AFTER ASAP PUBLICATION This paper was published ASAP on February 3, 2017 with mistakes in Schemes 1, 2, and 3. The corrected paper was reposted on February 7, 2017.
1453
DOI: 10.1021/acs.jctc.6b01141 J. Chem. Theory Comput. 2017, 13, 1439−1453