Case-specific performance of MM-PBSA, MM-GBSA, and SIE in virtual screening

Journal of Molecular Graphics and Modelling 62 (2015) 303–318

Contents lists available at ScienceDirect

Journal of Molecular Graphics and Modelling journal homepage: www.elsevier.com/locate/JMGM

Case-specific performance of MM-PBSA, MM-GBSA, and SIE in virtual screening Salla I. Virtanen 1,2 , Sanna P. Niinivehmas 1 , Olli T. Pentikäinen ∗ University of Jyvaskyla, Computational Bioscience Laboratory, Department of Biological and Environmental Science & Nanoscience Center, P.O. Box 35, FI-40014 University of Jyvaskyla, Finland

a r t i c l e

i n f o

Article history: Received 5 August 2015 Received in revised form 23 October 2015 Accepted 26 October 2015 Available online 30 October 2015 Keywords: Molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) Molecular mechanics generalized Born surface area (MM-GBSA) Solvated interaction energy (SIE) Virtual screening Molecular dynamics simulation Progesterone receptor Aldose reductase 2 Ampc beta-lactamase Heat shock protein 90 Phosphodiesterase Docking Shape comparison

a b s t r a c t In drug discovery the reliable prediction of binding free energies is of crucial importance. Methods that combine molecular mechanics force fields with continuum solvent models have become popular because of their high accuracy and relatively good computational efficiency. In this research we studied the performance of molecular mechanics generalized Born surface area (MM-GBSA), molecular mechanics Poisson–Boltzmann surface area (MM-PBSA), and solvated interaction energy (SIE) both in their virtual screening efficiency and their ability to predict experimentally determined binding affinities for five different protein targets. The protein-ligand complexes were derived with two different approaches important in virtual screening: molecular docking and ligand-based similarity search methods. The results show significant differences between the different binding energy calculation methods. However, the length of the molecular dynamics simulation was not of crucial importance for accuracy of results. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Virtual screening (VS) of drug molecules is becoming a standard approach in drug discovery. The major benefit of VS is that large compound databases can be screened time- and cost-effectively, allowing the laborious and expensive experimental tests to focus only on the top fraction of the screened molecules. Of particular interest in the drug discovery process is the reliable prediction of binding affinities of small molecules to the target protein, which is important both in the early VS phase and later in the optimization of drug candidates. In VS the binding energy prediction has traditionally relied upon empirical scoring functions implemented in molecular docking software, which can be used to screen relatively

∗ Corresponding author. Fax: +358 14 617 239. E-mail address: olli.t.pentikainen@jyu.fi (O.T. Pentikäinen). 1 Authors contributed equally to this work. 2 Present address: School of Chemical Technology, Aalto University, P.O. Box 16100, FI-00076 Aalto, Finland. http://dx.doi.org/10.1016/j.jmgm.2015.10.012 1093-3263/© 2015 Elsevier Inc. All rights reserved.

large numbers of molecules in a reasonable computation time. Many studies, however, have shown that the approximated scoring functions often produce results that poorly correlate with experimental binding data [1–4]. More rigorous methods for binding affinity prediction have been developed, such as thermodynamic integration (TI) [5] and free energy perturbation (FEP) [6], but computationally more efficient methods are desperately needed for VS of large compound databases. A reasonable compromise between computational efficiency and accuracy has been achieved in methodologies that combine molecular mechanics force fields and implicit solvation models. These methods include molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) [7], molecular mechanics generalized Born surface area (MM-GBSA) [8–10], and a more recently developed approach using a similar methodology called solvated interaction energy (SIE) [11]. In these methods, the solvation energy of molecules is taken into account by replacing the solvent with a dielectric continuum to account for the electrostatic contribution to solvation, along with terms that describe the nonpolar contribution. The major difference between SIE and MM-GB(PB)SA is

304

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

that the former has been calibrated on binding affinities in solution using a diverse dataset of 99 protein-ligand complexes [11]. The MM-GB(PB)SA and SIE methods have shown promising results in their ability to correctly rank the molecules with known affinity to their target proteins [12–20], and in the case of MMGB(PB)SA, results in the ability to distinguish between active and inactive molecules [14,21–24]. Furthermore, MM-GBSA has been successfully used in the prediction of the site of metabolism in ligands [25]. However, it has also been shown that the results are often case-specific, and depend heavily on the protocol used to create the protein–ligand complexes, parameters used in the molecular dynamics (MD) simulations, and parameters in the energy calculations [13,14]. In this research, the following issues in MMPBSA, MM-GBSA, and SIE were studied: (1) different approaches in preparation of the protein–ligand complex, (2) comparison of VS efficiencies of the methods, (3) prediction of the binding affinity, and (4) the effect of the length of the MD simulation to the results. We chose five different types of proteins as targets: aldose reductase-2 (ALR-2), AmpC ␤-lactamase (AmpC), human heat shock protein 90 (HSP90), phosphodiesterase-5 (PDE-5), and progesterone receptor (PR). The target proteins are known to undergo conformational changes upon binding, which are sometimes beyond short MD simulations [26–31]. For docking the conformational change poses a problem, because often only one static protein conformation is used, which might not be optimal for all the ligands. Accordingly, the joint usage of MD simulations with MM-GB(PB)SA could identify a more favorable conformation for the protein–ligand complex [32]. Thus, we also studied whether MMGB(PB)SA can improve the results for flexible protein targets. To be able to study these issues, appropriate datasets of known active molecules and decoy molecules are required. The Directory of Useful Decoys (DUD) [33] was used in the VS studies, and molecules with known IC50 values for the protein targets were chosen from the CHEMBL database [34].

The active ligands with known IC50 values for the target proteins were downloaded from the CHEMBL database [34]. The structures for each protein target were chosen from a single research article to ensure similar testing environments. The molecules and their experimental affinities are shown in Table S1 and the structures in Fig. S1 in the Supporting information. The molecules were converted from SMILES to 3D with LigPrep in Schrödinger molecular modeling package. All tautomers were generated in addition to ionization states within pH 7.4 ± 0.0. Defined stereochemistry was preserved in the structures. Low-energy conformers for the ligandbased superimposition were generated with ConfGen as described above. 2.3. Preparation of the protein–ligand complexes MM-GB(PB)SA and SIE methods are well suited for postprocessing the initial VS results. Therefore we wanted to test the performance of these methods in rescoring complexes generated with VS methods. Molecular docking is traditionally used as a tool to generate the protein–ligand complexes for MM-GB(PB)SA. However, we wanted to test also another method for creating the starting conformations, namely ligand-based superimpositioning of molecules. These superimposed or “docked” structures could then be used as an alternative to structures derived with docking. 2.4. Docking For molecular docking, GOLD 4.1 was employed. Default settings for the genetic algorithm were used allowing 10 docking runs for each molecule, out of which the best conformation according to GOLD score was chosen for post-processing. The binding site was defined by choosing an atom in a central location within the binding pocket. The water-fill radius of 10 Å for the binding site was used. 2.5. SHAEP

2. Materials and methods 2.1. Protein structures The protein structures for ALR-2 (pdb: 1ah3 [35], resolution 2.30 Å), AmpC (pdb: 1xgj [36], resolution 1.97 Å), HSP90 (pdb: 1uy6 [37], resolution 1.90 Å), PDE-5 (pdb: 1xp0 [38], resolution 1.79 Å), and PR (pdb: 1sr7 [27], resolution 1.46 Å) were downloaded from the PDB. All water molecules were removed with BODIL [39], except for the seven structural waters preserved also in the DUD PDE-5 protein structure. The metal ions in PDE-5, Zn2+ and Mg2+ , were left in their places for docking, as was the cofactor nicotinamide adenine dinucleotide phosphate (NADP+ ) in the ALR-2 structure. Hydrogen atoms were added to the structures with TLEAP available in ANTECHAMBER [40].

2.2. Ligand and decoy sets Ligands and decoy molecules for VS studies were downloaded from the DUD database [33]. The database contains 946 molecules for ALR2 (26 ligands, 920 decoys), 755 molecules for AmpC (21 ligands, 734 decoys), 885 molecules for HSP90 (24 ligands, 861 decoys), 1861 molecules for PDE5 (51 ligands, 1810 decoys), and 994 molecules for PR (27 ligands, 967 decoys). For docking the molecules were used as provided. For ligand-based superimposition low-energy conformers were generated for each molecule with ConfGen routine [41] available in Schrödinger molecular modeling package. The Fast option was used with Merck Molecular Force Field (MMFF) [42].

SHAEP [43] compares the shape and electrostatic potential similarity of molecules, and superimposes them according to the maximum overlap. The ligand from the protein crystal structure was extracted to be used as the target molecule in superimposition. MMFF94 charges were assigned to the crystal structure ligands. Default settings in SHAEP, where shape and electrostatics are equally contributing to the results were used. 2.6. MD simulations The SANDER module distributed in AMBER 10 [44] package was used for MD simulations. The protein–ligand complex generated with either GOLD or SHAEP was used as a starting conformation for MD simulations. Charges for the ligands and for NADP+ were derived with AM1-BCC [45] available in ANTECHAMBER. TLEAP was used to create force field parameters for the protein (ff03) [46] and the ligand/cofactor (gaff) [47], add hydrogens, and solvate the ligand–protein complex with a rectangular box of transferable intermolecular potential three-point water molecules (TIP3P) [48] 14 Å in all directions from the protein. Histidine residues were assigned the HIE state (proton at the -nitrogen), except the PDE-5 residues His617 and His653, which were treated as HID (proton at the ␦-nitrogen). The MD simulations were run as previously described [23]. Briefly, the system was first minimized with conjugate-gradient method for 1000 steps without restraints. This was followed by an equilibration step at constant volume by allowing the system to heat from 100 K to 300 K for 1000 steps with NMR restraints. The production simulation without restraints was run for 256,000 steps (simulation time 512 ps) at constant pressure. Temperature and pressure were maintained with the

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

Berendsen weak coupling method [49]. Electrostatics were treated with Particle-Mesh Ewald (PME) method [50,51], and cutoff value of 12 Å for nonbonded interactions was employed. The equilibration step and the production simulation were run under periodic boundary conditions. The SHAKE algorithm [52] was used for bonds involving hydrogens, which allowed the use of 2 fs time step.

2.8. SIE SIE calculations were performed with the software SIETRAJ. The calculations are done with the Eq. (3). Gbind (␳, Din , ␣, ␥, C) = ␣ × [EC (Din ) + Gbind R (␳, Din ) + E vdw + ␥ × MSA(␳)] + C

2.7. MM-GBSA and MM-PBSA In the MM-GBSA and MM-PBSA calculations the free energy of the ligand binding is estimated by taking into account the solvation energies of the interacting molecules in addition to molecular mechanics (MM) energies. The contribution of polar solvation energy is calculated with the implicit solvent model (GB or PB); whereas the nonpolar part of the solvation energy is dependent on the solvent accessible surface area (SA). The free energies of binding (Gbind ) can be estimated from the free energies of the three reactants with the Eq. (1). Gbind = < Gcomp > − < Grec > − < Glig >

(1)

where < > denotes an average over a set of snapshots along the MD trajectory, and comp, rec, and lig stand for complex, receptor, and ligand, respectively. The free energy for each of the reactants is estimated with the Eq. (2). G = EMM + Gsolv − TSsolute

(2)

where EMM is the molecular mechanics contribution in vacuo consisting of the sum of internal, electrostatics, and van der Waals energies; Gsolv is the contribution of solvation free energies expressed as the sum of polar (GGB ) and nonpolar (Gnonpolar ) solvation free energies; T is the temperature; and Ssolute is the solute entropy. Nonpolar solvation free energy can be estimated with Gnonpolar = ␥ SAS, where ␥ is the surface tension constant and SAS is the solvent accessible surface. For the binding energy analyses snapshots at 4 ps intervals were extracted from the MD trajectory, and the free energies were averaged over the ensemble of conformers produced (128 snapshots for each reactant). The polar solvation energies were calculated with the PB and GB approaches available in AMBER 10. Dielectric constants of 1 and 80 were used for the solute and the solvent, respectively. Although it has been reported previously that a higher solute dielectric constant 2 or 4 may produce more accurate results for some targets [18,20], here we wanted to use the default settings for each of the tested methods, since this is the common first approach, particularly if there is no experimental data available. The three different GB models available in AMBER 10 were tested: GB model with parameters developed by Tsui and Case [9] (IGB = 1, later referred to IGB1), and two GB models developed by Onufriev et al. [8] (IGB = 2, later referred to IGB2, and IGB = 5, later referred to IGB5). The hydrophobic contribution to the solvation free energy was estimated by calculating the solvent accessible surface area with Molsurf [53] and using a probe radius of 1.4 Å. The surface tension constant ␥ was set to 0.0072 kcal/mol/Å2 . The entropy term in the binding free energy is usually estimated with normal mode analysis (NMA) [54] of the vibration frequencies. Because of the high computational cost of NMA, the effect of entropy was neglected in the calculations. Earlier studies show that the estimation of entropy can be highly misleading, as the motion of protein cannot be properly studied, especially in the nanosecond timescale.

305

(3)

where EC and Evdw are the intermolecular Coulomb and van der Waals interaction energies in the bound state, respectively. Gbind R is the change in the reaction free energy between the bound and free states respectively, calculated by the Poisson equation. The MSA term is the change in molecular surface area upon binding. The coefficients ␳ (AMBER van der Waals radii linear scaling coefficient), ␣ (global proportionality coefficient relating to the loss of configurational entropy upon binding), ␥ (molecular surface area coefficient), Din (solute internal dielectric constant) and constant C were optimized with the set of 99 protein-ligand complexes. The optimized values are:  = 1.1, ˛ = 0.1048,  = 0.0129 kcal/mol/Å2 , Din = 2.25, and C = −2.89 kcal/mol [9]. The same 128 MD snapshot structures extracted in 4 ps intervals were used also in SIE calculations. Default settings in SIETRAJ were used. 2.9. Metrics To evaluate the VS efficiency of the methods, Receiver Operating Characteristics (ROC) curves and Area Under the Curve (AUC) values were used. The ROC curves and AUC values with standard errors were calculated with SPSS, version 21 (IBM Corp.). To evaluate the early enrichment of the VS results, enrichment factors were calculated for the top 1%, 5%, and 10% of the ranked molecular databases. The enrichment factors (EFn% ) were calculated according to the following Eq. (4) EFn% =

Ligsn% × Molsall Molsn% × Ligsall

(4)

where Ligsn% and Molsn% are the number of ligands and molecules in the top n% of the ranked compounds. Ligsall and Molsall are the number of ligands and molecules, respectively, in the entire database. In order to determine the effect of the length of the MD simulation on the VS success, the AUC values and enrichment factors were calculated for structures obtained from the beginning of the simulation to the following time steps: at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps. 3. Results 3.1. Virtual screening efficiency for ALR-2 In ALR-2 the binding site is near the surface of the protein and exposed to solvent (Fig. 1A and B). The ALR-2 binding pocket is positively charged because of the oxidized nicotinamide ring of the cofactor NADP+ . Thus, the inhibitors contain either negatively charged or polar substructures. Important hydrogen bonding occurs with residues Tyr48, His110, and Trp111 (Fig. 1A). Besides the positive charge, the binding pocket also contains hydrophobic amino acid residues. Fig. 1B shows the PB potential surface of the binding site. The docking results with GOLD and GoldScore show very low VS efficiency overall when the AUC value is considered (AUC = 0.57 ± 0.08; Fig. 2A). However, the early enrichments at 1% and 5% are fairly good (EF1% = 20.2, EF5% = 7.0, EF10% = 4.2). The enrichment obtained with SHAEP was worse compared to docking, giving an AUC value of 0.36 ± 0.05 and very low early enrichment (EF1% = 4.0, EF5% = 1.6, EF10% = 1.2; Fig. 2B).

306

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

Fig. 1. The binding sites of the targets ALR-2, AmpC, HSP90, PDE-5, and PR (A, C, E, G, and I), and the PB potential surfaces calculated in silico (B, D, F, H, and J). The inhibitors are shown with black carbon atoms, and the key amino acids important for binding are highlighted in purple. Hydrogen bonds are depicted with orange dashes. (A) The binding site of ALR-2, with the inhibitor tolrestat. The cofactor NADP+ is shown with green carbon atoms. (B) The surface of the binding site of ALR-2 shown with the PB potential. (C) The binding site of AmpC with the inhibitor 3-[[4-carboxy-2-hydroxyaniline]sulfonyl]thiophene-2-carboxylic acid. (D) The binding site of AmpC shown with PB potential surface. (E) The binding site of HSP90 with the inhibitor 9-butyl-8-(3,4,5-trimethoxybenzyl)-9h-purin-6-amine. (F) The surface of the binding site of HSP90 shown with the PB potential. (G) The binding site of PDE-5 with the inhibitor vardenafil. Also visible are the metal ions Mg2+ (yellow CPK) and Zn2+ (red CPK), along with the structural water molecules. (H) The binding site of PDE-5 shown with PB potential surface. (I) The binding site of PR with the inhibitor mometasone furoate. (J) The binding cavity of PR shown with the PB potential surface. The figures were prepared with PYMOL [55], and the PB potential surfaces were calculated with Adaptive Poisson–Boltzmann Solver (APBS) [56] implemented in PYMOL. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

307

Fig. 2. ROC curves for (A and B) ALR-2, (C and D) AMPC, (E and F) HSP90, (G and H) PDE-5, and (I and J) PR for the free energy calculation methods for structures obtained at 4 ps MD simulation. Se = sensitivity, 1-Sp = specificity. Figure was prepared with Rocker [57].

308

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

Table 1 AUC values and enrichment factors for ALR-2 with MD simulations of different lengths. The enrichment factors were calculated at 1%, 5%, and 10% of the ranked database. The number of found hits is in parentheses after the enrichment factor. The results were calculated for structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps. The best result(s) in each row is/are highlighted in bold. ALR-2

4 ps 1% 5% 10% 32 ps 1% 5% 10% 64 ps 1% 5% 10% 128 ps 1% 5% 10% 256 ps 1% 5% 10% 384 ps 1% 5% 10% 512 ps 1% 5% 10%

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.43 ± 0.06 4.0 (1) 2.3 (3) 1.2 (3) 0.45 ± 0.07 4.0 (1) 2.3 (3) 1.5 (4) 0.47 ± 0.07 0 1.6 (2) 1.5 (4) 0.46 ± 0.06 4.0 (1) 1.6 (2) 1.5 (4) 0.48 ± 0.06 0 1.6 (2) 1.2 (3) 0.48 ± 0.06 4.0 (1) 2.3 (3) 1.2 (3) 0.50 ± 0.06 4.0 (1) 2.3 (3) 1.5 (4)

0.48 ± 0.06 8.1 (2) 1.6 (2) 1.9 (5) 0.49 ± 0.06 4.0 (1) 2.3 (3) 1.2 (3) 0.47 ± 0.05 4.0 (1) 2.3 (3) 1.2 (3) 0.43 ± 0.05 4.0 (1) 1.6 (2) 1.2 (3) 0.43 ± 0.05 8.1 (2) 1.6 (2) 1.2 (3) 0.44 ± 0.06 8.1 (2) 1.6 (2) 1.9 (5) 0.45 ± 0.06 8.1 (2) 2.3 (3) 1.5 (4)

0.39 ± 0.06 0 1.6 (2) 1.2 (3) 0.41 ± 0.06 0 1.6 (2) 1.2 (3) 0.43 ± 0.06 0 1.6 (2) 0.8 (2) 0.43 ± 0.06 0 1.6 (2) 0.8 (2) 0.44 ± 0.06 0 0.8 (1) 0.8 (2) 0.45 ± 0.06 0 1.6 (2) 1.2 (3) 0.46 ± 0.06 0 1.6 (2) 1.2 (3)

0.46 ± 0.06 4.0 (1) 1.6 (2) 0.8 (2) 0.45 ± 0.06 4.0 (1) 1.6 (2) 1.2 (3) 0.44 ± 0.05 4.0 (1) 1.6 (2) 1.2 (3) 0.40 ± 0.05 0 0.8 (1) 0.8 (2) 0.41±0.05 4.0 (1) 1.6 (2) 0.8 (2) 0.42 ± 0.05 8.1 (2) 1.6 (2) 1.2 (3) 0.43 ± 0.05 4.0 (1) 1.6 (2) 1.2 (3)

0.64 ± 0.07 12.1 (3) 7.7 (10) 3.8 (10) 0.65 ± 0.07 12.1 (3) 7.7 (10) 3.8 (10) 0.65 ± 0.07 12.1 (3) 7.0 (9) 3.8 (10) 0.64 ± 0.07 12.1 (3) 7.0 (9) 3.8 (10) 0.64 ± 0.07 8.1 (2) 7.0 (9) 4.2 (11) 0.64 ± 0.07 8.1 (2) 7.0 (9) 4.2 (11) 0.63 ± 0.07 4.0 (1) 6.2 (8) 4.2 (11)

0.59 ± 0.07 4.0 (1) 4.6 (6) 3.5 (9) 0.54 ± 0.07 4.0 (1) 4.6 (6) 3.5 (9) 0.55 ± 0.07 4.0 (1) 4.6 (6) 3.1 (8) 0.54 ± 0.07 0 4.6 (6) 3.1 (8) 0.54 ± 0.07 0 4.6 (6) 3.1 (8) 0.54 ± 0.07 0 4.6 (6) 3.1 (8) 0.54 ± 0.07 0 3.9 (5) 3.1 (8)

0.60 ± 0.07 20.2 (5) 5.4 (7) 3.1 (8) 0.61 ± 0.07 16.2 (4) 6.2 (8) 3.8 (10) 0.61 ± 0.07 16.2 (4) 5.4 (7) 3.5 (9) 0.60 ± 0.07 16.2 (4) 6.2 (8) 3.8 (10) 0.60 ± 0.07 20.2 (5) 6.2 (8) 3.1 (8) 0.60 ± 0.07 20.2 (5) 6.2 (8) 3.1 (8) 0.60 ± 0.07 20.2 (5) 6.2 (8) 3.1 (8)

0.58 ± 0.06 8.1 (2) 4.6 (6) 2.3 (6) 0.59 ± 0.06 8.1 (2) 2.3 (3) 2.3 (6) 0.60 ± 0.06 4.0 (1) 2.3 (3) 2.3 (6) 0.57 ± 0.06 4.0 (1) 2.3 (3) 2.3 (6) 0.56 ± 0.06 4.0 (1) 2.3 (3) 2.7 (7) 0.55 ± 0.06 4.0 (1) 3.1 (4) 2.7 (7) 0.55 ± 0.06 4.0 (1) 3.1 (4) 3.1 (8)

0.60 ± 0.07 8.1 (2) 5.4 (7) 3.5 (9) 0.60 ± 0.08 4.1 (1) 6.2 (8) 4.2 (11) 0.58 ± 0.07 12.1 (3) 7.0 (9) 3.8 (10) 0.59 ± 0.07 16.2 (4) 7.0 (9) 3.8 (10) 0.59 ± 0.07 16.2 (4) 6.2 (8) 3.8 (10) 0.59 ± 0.07 16.2 (4) 7.0 (9) 3.8 (10) 0.59 ± 0.07 16.2 (4) 6.2 (8) 3.8 (10)

0.56 ± 0.06 12.1 (3) 5.4 (7) 3.1 (8) 0.58 ± 0.06 8.1 (2) 5.4 (7) 3.1 (8) 0.57 ± 0.06 4.0 (1) 3.9 (5) 2.3 (6) 0.55 ± 0.06 4.0 (1) 3.9 (5) 1.9 (5) 0.54 ± 0.06 4.0 (1) 3.1 (4) 1.9 (5) 0.53 ± 0.06 4.0 (1) 3.1 (4) 1.9 (5) 0.53 ± 0.06 4.0 (1) 3.1 (4) 1.9 (5)

For the complexes generated with GOLD, MM-GBSA with IGB5 parameters produced the highest AUC value of 0.65 ± 0.07 (Table 1, Fig. 2A), which is, in contrast to the used computational efforts, a relatively poor result. The very early enrichment at 1% is worse compared to docking, but the enrichments at 5% and 10% are comparable (Table 1). The other two MM-GBSA methods, IGB1 and IGB2, respectively, did not succeed in VS producing only AUC values of 0.50 ± 0.06 and 0.46 ± 0.06 at best, which is comparable to random selection of the molecules. The early enrichments are also very low for IGB1 and IGB2 (Table 1). MM-PBSA and SIE produced similar results compared to each other, giving AUC values of 0.61 ± 0.07 and 0.60 ± 0.07, respectively, with rather good early enrichment (Table 1). When SHAEP was used for creation of the protein–ligand complexes the results were similar for IGB1 and IGB2; however, the results obtained with IGB5, MM-PBSA, and SIE were worse than those obtained from docked ligand conformations (Table 1, Fig. 2B). Accordingly, from used methods the most suitable ones for ALR-2 are docking alone, and docking with post-processing with IGB5. 3.2. Virtual screening efficiency for AmpC The binding site for AmpC inhibitors is close to the surface and exposed to the solvent (Fig. 1C and D). The binding site contains both hydrophobic and hydrophilic residues. The VS efficiency with docking produced a very low AUC value, 0.42 ± 0.06, with also very poor early enrichment (EF1% = 0, EF5% = 0.9, EF10% = 0.5; Fig. 2C). Shape-based VS with SHAEP performed better producing an AUC value of 0.76 ± 0.06 and the enrichment factors were clearly higher compared to docking (EF1% = 17.9, EF5% = 6.6, EF10% = 3.8; Fig. 2D). All post-processing methods clearly improved the VS efficiency in terms of AUC values for the GOLD-generated structures (Table 2). Of these, MM-GBSA with the IGB5 parameters produced the highest AUC value 0.82 ± 0.04 as well as the highest enrichment factors.

Similarly to the GOLD-generated structures, the IGB5 produced the best VS results also for the structures generated with SHAEP (AUC = 0.89 ± 0.03; Table 2). Also the early enrichment was the best with IGB5, except at EF1% which was the best for the shape-based VS without any post-processing (17.9 vs 13.5). Taken together, the best results were obtained with SHAEP-generated structures and post-processing with IGB5. 3.3. Virtual screening efficiency for HSP90 The HSP90 inhibitors bind to the nucleotide binding pocket, which is located close to the surface of the protein (Fig. 1E and F). The binding site is composed of hydrophobic as well as hydrophilic residues, of which Asp93 forms key hydrogen bonding interactions with the inhibitors. The VS efficiency with docking produced an AUC value of 0.50 ± 0.05, which equals random picking of molecules. Additionally, the enrichment factors show that none of the active ligands were found among the top 10% of the ranked molecules (EF1% = 0, EF5% = 0, EF10% = 0; Fig. 2E). SHAEP, on the other hand, performed better in terms of the AUC value (0.64 ± 0.07) and produced reasonably good early enrichment (EF1% = 24.6, EF5% = 8.4, EF10% = 4.2; Fig. 2F). The post-processing of the structures generated with GOLD improved the AUC results only in the case of IGB5 (AUC = 0.64 ± 0.05; Table 3). In fact, the VS efficiency decreases with the longer MD simulation for all the other methods, going as low as below AUC = 0.40 (Table 3). Nevertheless, according to the early enrichment, SIE and IGB1 performed better. For the SHAEPgenerated structures the results show similar effects: IGB5 was the most successful in VS in terms of AUC values (AUC = 0.70 ± 0.05), whereas the early enrichment was the best with IGB1. However, in terms of the EF values the shape-based VS was clearly the most efficient without any post-processing.

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

309

Table 2 AUC values and enrichment factors for AmpC with MD simulations of different lengths. The enrichment factors were calculated at 1%, 5%, and 10% of the ranked database. The number of found hits is in parenthesis after the enrichment factor. The results were calculated for structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps. The best result(s) in each row is/are highlighted in bold. AmpC

4 ps 1% 5% 10% 32 ps 1% 5% 10% 64 ps 1% 5% 10% 128 ps 1% 5% 10% 256 ps 1% 5% 10% 384 ps 1% 5% 10% 512 ps 1% 5% 10%

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.67 ± 0.06 0 1.9 (2) 2.4 (5) 0.65 ± 0.05 0 0.9 (1) 1.4 (3) 0.65 ± 0.06 0 1.9 (2) 1.4 (3) 0.67 ± 0.06 0 1.9 (2) 1.9 (4) 0.60 ± 0.06 0 1.9 (2) 1.4 (3) 0.57 ± 0.06 0 0.9 (1) 1.4 (3) 0.56 ± 0.05 0 1.9 (2) 1.4 (3)

0.79 ± 0.03 0 2.8 (3) 2.9 (6) 0.83 ± 0.03 4.5 (1) 1.9 (2) 3.4 (7) 0.77 ± 0.04 4.5 (1) 2.8 (3) 3.8 (8) 0.77 ± 0.05 4.5 (1) 1.9 (2) 2.9 (6) 0.70 ± 0.05 0 3.8 (4) 2.4 (5) 0.67 ± 0.05 4.5 (1) 2.8 (3) 1.9 (4) 0.66 ± 0.05 4.5 (1) 2.8 (3) 1.9 (4)

0.70 ± 0.06 4.5 (1) 3.8 (4) 2.9 (6) 0.69 ± 0.05 0 1.9 (2) 1.9 (4) 0.69 ± 0.05 0 2.8 (3) 1.9 (4) 0.71 ± 0.05 0 3.8 (4) 2.4 (5) 0.64 ± 0.05 0 2.8 (3) 1.4 (3) 0.60 ± 0.05 4.5 (1) 0.9 (1) 1.0 (2) 0.59 ± 0.05 0 0.9 (1) 1.4 (3)

0.86 ± 0.03 0 4.7 (5) 4.3 (9) 0.88 ± 0.02 0 4.7 (5) 4.8 (10) 0.83 ± 0.04 0 3.8 (4) 4.3 (9) 0.82 ± 0.04 0 3.8 (4) 4.8 (10) 0.75 ± 0.05 0 3.8 (4) 2.4 (5) 0.72 ± 0.05 0 2.8 (3) 2.4 (5) 0.71 ± 0.05 0 3.8 (4) 1.9 (4)

0.81 ± 0.04 9.0 (2) 4.7 (5) 5.3 (11) 0.82 ± 0.04 4.5 (1) 3.8 (4) 3.8 (8) 0.81 ± 0.04 4.5 (1) 3.8 (4) 3.8 (8) 0.82 ± 0.04 4.5 (1) 4.7 (5) 4.8 (10) 0.80 ± 0.04 0 3.8 (4) 4.3 (9) 0.76±0.05 4.5 (1) 3.8 (4) 2.9 (6) 0.74 ± 0.05 0 2.8 (3) 2.4 (5)

0.88 ± 0.03 9.0 (2) 9.4 (10) 5.7 (12) 0.89 ± 0.03 13.5 (3) 9.4 (10) 5.7 (12) 0.87 ± 0.03 13.5 (3) 8.5 (9) 5.3 (11) 0.86 ± 0.03 13.5 (3) 8.5 (9) 5.7 (12) 0.83±0.04 4.5 (1) 5.7 (6) 4.8 (10) 0.82 ± 0.04 4.5 (1) 4.7 (5) 3.4 (7) 0.81 ± 0.04 4.5 (1) 2.8 (3) 3.8 (8)

0.63 ± 0.07 0 2.8 (3) 3.4 (7) 0.62 ± 0.06 0 3.8 (4) 2.4 (5) 0.62 ± 0.06 4.5 (1) 3.8 (4) 2.4 (5) 0.64 ± 0.06 9.0 (2) 2.8 (3) 2.4 (5) 0.59 ± 0.06 9.0 (2) 1.9 (2) 1.4 (3) 0.57 ± 0.06 4.5 (1) 2.8 (3) 1.4 (3) 0.57 ± 0.06 0 1.9 (2) 1.9 (4)

0.73 ± 0.05 4.5 (1) 4.7 (5) 3.4 (7) 0.78 ± 0.04 4.5 (1) 3.8 (4) 3.8 (8) 0.76 ± 0.04 4.5 (1) 4.7 (5) 3.4 (7) 0.76 ± 0.05 4.5 (1) 3.8 (4) 3.4 (7) 0.72 ± 0.05 4.5 (1) 2.8 (3) 2.9 (6) 0.70 ± 0.05 9.0 (2) 2.8 (3) 1.4 (3) 0.69 ± 0.05 4.5 (1) 2.8 (3) 1.4 (3)

0.62 ± 0.07 0 0.9 (1) 2.9 (6) 0.60 ± 0.06 4.5 (1) 2.8 (3) 1.9 (4) 0.60 ± 0.06 4.5 (1) 1.9 (2) 1.9 (4) 0.61 ± 0.06 4.5 (1) 1.9 (2) 1.9 (4) 0.55 ± 0.05 4.5 (1) 1.9 (2) 1.0 (2) 0.54 ± 0.05 0 0.9 (1) 0.5 (1) 0.53 ± 0.05 0 0 0.5 (1)

0.74 ± 0.05 0 3.8 (4) 2.9 (6) 0.75 ± 0.04 9.0 (2) 5.7 (6) 2.9 (6) 0.72 ± 0.05 13.4 (3) 5.7 (6) 2.9 (6) 0.71 ± 0.05 4.5 (1) 4.7 (5) 2.4 (6) 0.66 ± 0.05 4.5 (1) 1.9 (2) 2.4 (5) 0.64 ± 0.05 9.0 (2) 2.8 (3) 1.4 (3) 0.63 ± 0.05 4.5 (1) 2.8 (3) 1.4 (3)

Table 3 AUC values and enrichment factors for HSP90 with MD simulations of different lengths. The enrichment factors were calculated at 1%, 5%, and 10% of the ranked database. The number of found hits is in parenthesis after the enrichment factor. The results were calculated for structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps. The best result(s) in each row is/are highlighted in bold. HSP90

4 ps 1% 5% 10% 32 ps 1% 5% 10% 64 ps 1% 5% 10% 128 ps 1% 5% 10% 256 ps 1% 5% 10% 384 ps 1% 5% 10% 512 ps 1% 5% 10%

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.51 ± 0.07 12.3 (3) 3.4 (4) 2.5 (6) 0.44 ± 0.07 12.3 (3) 3.4 (4) 1.7 (4) 0.42 ± 0.07 12.3 (3) 3.4 (4) 2.1 (5) 0.40 ± 0.07 12.3 (3) 4.2 (5) 2.1 (5) 0.40 ± 0.07 12.3 (3) 3.4 (4) 2.1 (5) 0.39 ± 0.07 16.4 (4) 3.4 (4) 2.1 (5) 0.39 ± 0.07 16.4 (4) 3.4 (4) 2.1 (5)

0.50 ± 0.07 4.1 (1) 3.4 (4) 2.1 (5) 0.46 ± 0.07 8.2 (2) 1.7 (2) 1.3 (3) 0.45 ± 0.07 8.2 (2) 1.7 (2) 1.3 (3) 0.43 ± 0.07 8.2 (2) 1.7 (2) 1.3 (3) 0.41 ± 0.07 8.2 (2) 1.7 (2) 1.3 (3) 0.41 ± 0.07 8.2 (2) 1.7 (2) 1.3 (3) 0.41 ± 0.07 4.1 (1) 1.7 (2) 1.3 (3)

0.54 ± 0.07 8.2 (2) 3.4 (4) 2.1 (5) 0.45 ± 0.06 12.3 (3) 3.4 (4) 1.7 (4) 0.43 ± 0.06 12.3 (3) 3.4 (4) 2.1 (5) 0.41 ± 0.07 8.2 (2) 3.4 (4) 2.1 (5) 0.40 ± 0.07 8.2 (2) 3.4 (4) 2.1 (5) 0.40 ± 0.07 8.2 (2) 2.5 (3) 2.1 (5) 0.39 ± 0.07 8.2 (2) 2.5 (3) 1.7 (4)

0.49 ± 0.07 4.1 (1) 1.7 (2) 1.3 (3) 0.44 ± 0.07 4.1 (1) 0.8 (1) 1.3 (3) 0.44 ± 0.07 4.1 (1) 0.8 (1) 1.3 (3) 0.41 ± 0.07 4.1 (1) 0.8 (1) 0.8 (2) 0.40 ± 0.07 4.1 (1) 0.8 (1) 0.8 (2) 0.39 ± 0.07 4.1 (1) 0.8 (1) 0.8 (2) 0.38 ± 0.07 4.1 (1) 0.8 (1) 0.4 (1)

0.64 ± 0.06 0 1.7 (2) 1.7 (4) 0.64 ± 0.05 4.1 (1) 0.8 (1) 1.3 (3) 0.64 ± 0.05 4.1 (1) 0.8 (1) 0.8 (2) 0.64 ± 0.05 4.1 (1) 0.8 (1) 2.1 (5) 0.64 ± 0.05 4.1 (1) 0.8 (1) 2.1 (5) 0.64 ± 0.05 4.1 (1) 1.7 (2) 1.7 (4) 0.64 ± 0.05 4.1 (1) 1.7 (2) 2.1 (5)

0.70 ± 0.05 0 0.8 (1) 0.8 (2) 0.69 ± 0.05 0 0.8 (1) 0.4 (1) 0.69 ± 0.05 0 0 0.8 (2) 0.68 ± 0.05 0 0 1.3 (3) 0.68 ± 0.05 0 0 0.8 (2) 0.67 ± 0.05 0 0 0.8 (2) 0.66 ± 0.05 0 0 0.8 (2)

0.54 ± 0.07 4.1 (1) 1.7 (2) 1.7 (4) 0.49 ± 0.06 8.2 (2) 3.4 (4) 1.7 (4) 0.46 ± 0.06 12.3 (3) 2.5 (3) 1.7 (4) 0.45 ± 0.07 8.2 (2) 2.5 (3) 1.7 (4) 0.45 ± 0.06 12.3 (3) 2.5 (3) 1.7 (4) 0.44 ± 0.06 12.3 (3) 2.5 (3) 1.7 (4) 0.44 ± 0.06 8.2 (2) 2.5 (3) 1.7 (4)

0.50 ± 0.06 0 0.8 (1) 1.3 (3) 0.47 ± 0.07 4.1 (1) 0.8 (1) 1.7 (4) 0.46 ± 0.07 4.1 (1) 0.8 (1) 1.3 (3) 0.43 ± 0.06 4.1 (1) 0.8 (1) 0.8 (2) 0.41 ± 0.06 4.1 (1) 0.8 (1) 0.4 (1) 0.40 ± 0.06 4.1 (1) 0.8 (1) 0.4 (1) 0.39 ± 0.06 4.1 (1) 0.8 (1) 0.4 (1)

0.53 ± 0.06 12.3 (3) 4.2 (5) 2.5 (6) 0.46 ± 0.07 16.4 (4) 3.4 (4) 2.1 (5) 0.44 ± 0.07 16.4 (4) 4.2 (5) 2.1 (5) 0.42 ± 0.07 12.3 (3) 3.4 (4) 2.1 (5) 0.41 ± 0.07 12.3 (3) 2.5 (3) 1.7 (4) 0.40 ± 0.07 12.3 (3) 2.5 (3) 1.7 (4) 0.40 ± 0.06 12.3 (3) 2.5 (3) 1.7 (4)

0.53 ± 0.06 4.1 (1) 1.7 (2) 1.7 (4) 0.55 ± 0.06 4.1 (1) 1.7 (2) 1.7 (4) 0.53 ± 0.06 4.1 (1) 0.8 (1) 1.3 (3) 0.50 ± 0.06 4.1 (1) 0.8 (1) 1.3 (3) 0.48 ± 0.06 4.1 (1) 0.8 (1) 0.8 (2) 0.47 ± 0.06 4.1 (1) 0.8 (1) 0.4 (1) 0.46 ± 0.06 4.1 (1) 0.8 (1) 0.8 (2)

3.4. Virtual screening efficiency for PDE-5 In PDE-5, the binding site is very close to the surface of the protein (Fig. 1G and H). The binding site contains both polar and

hydrophobic amino acids that are important for binding and for interactions with explicit water molecules. Within the binding pocket Gln817 (conserved among all PDEs) forms key binding interactions with the inhibitors. PDE-5 is also a metalloenzyme con-

310

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

taining two metal ions, Mg2+ and Zn2+ , within the binding site. In VS studies the docking produced an AUC value of 0.69±0.04 with very low early enrichment (EF1% = 3.8, EF5% = 3.9, EF10% = 2.0; Fig. 2G). SHAEP was not as successful as GOLD in VS, producing only an AUC value of 0.49 ± 0.04. However, the early enrichment at 1% was better compared to docking (EF1% = 9.6, EF5% = 2.0, EF10% = 1.6; Fig. 2H). For the GOLD-generated structures the free energy calculations made with IGB1, IGB2, MM-PBSA, and SIE are roughly comparable in terms of the overall enrichment with AUC values between 0.70 and 0.75 (Table 4, Fig. 2G). The early enrichment, however, was the best for IGB1 and IGB2. In contrast to ALR-2, the MM-GBSA with the parameters from IGB5 produced very poor results compared to the other tested methods, with an AUC value of 0.53 ± 0.04 at best and the enrichment values mostly below 1.0 (Tables 2 and 5). With the SHAEP-created complexes IGB1 and IGB2 methods performed comparably to the GOLD-created structures when the AUC values were compared (Table 4, Fig. 2H). Notably, the early enrichment benefitted from the usage of SHAEP to create the starting geometries. For IGB5, MM-PBSA, and SIE, the results generally worsened somewhat. Taken together and considering the early enrichment, it is concluded that SHAEP, especially with post-processing with IGB1 or IGB2, yielded the best result. 3.5. Virtual screening efficiency for PR Progesterone receptor is a nuclear hormone receptor in which the binding site for the natural agonist ligand is buried inside the protein (Fig. 1I and J). The ligand-binding site contains mainly hydrophobic residues; however, important hydrogen bonding interactions occur with residues Gln725 and Arg766. For PR, the docking did not succeed in VS, with an AUC value of 0.37 ± 0.05, and the method could not identify any of the ligands in the first 10% of the ranked molecules (EF1% = 0, EF5% = 0, EF10% = 0; Fig. 2I). Although the enrichments for SHAEP were low (AUC = 0.59 ± 0.05; EF1% = 3.7, EF5% = 4.4, EF10% = 2.2; Fig. 2J), SHAEP performed better compared to GOLD. For the complexes generated with GOLD, all used postprocessing methods improved the results significantly, of which the best method, according to AUC values, was MM-GBSA with IGB5 parameters (0.72 ± 0.04) (Table 5, Fig. 2I). However, the early enrichment was rather poor for IGB5 (Table 5). The AUC values produced by MM-PBSA were nearly as good as with the IGB5 (0.69 ± 0.05), but the early enrichment was clearly better for MMPBSA (Table 5). For the other tested methods, IGB1, IGB2, and SIE, the AUC values showed weak overall performance; however, the early enrichments were better compared to IGB5. In contrast to GOLD, the structures generated with SHAEP produced a higher AUC value for IGB5 improving it from 0.72 ± 0.04 to 0.78 ± 0.05 (Table 5, Fig. 2J). Unfortunately, this improvement is not seen in the early enrichments, but is visible after the 10% in the form of a rapid rise of the AUC curve (Fig. 2J), which shows that more than half of the ligands are identified at that point. The results with IGB1, MM-PBSA, and SIE show relatively good early enrichment at 1%, and comparable enrichment up to 10%, but after that the performance is clearly worse than with IGB5. Overall, only the docking with GOLD clearly failed. This is most likely because the ligand set for PR contains molecules with highly versatile structures, and thus the one protein conformation is not optimal for all the ligands. 3.6. Effect of the MD simulation length The effect of the MD simulation length on the VS results was studied by calculating ROC AUC values and enrichment factors for structures extracted at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps (see Tables 2–6). Fig. 3 shows how the AUC values develop

along the 512 ps simulation time. The VS efficiency of the methods in terms of the AUC values stays at a fairly stable level throughout the 512 ps MD simulation, except for AmpC. The results thus imply that a short MD simulation can also be used in VS with the MMGB(PB)SA and SIE. 4. Estimation of binding free energies The methods were further evaluated by predicting binding affinities of molecules with known IC50 values. Both docking and SHAEP were used in the generation of the protein-ligand complexes and the correlations of the results to the experimental values were compared at the same MD simulation lengths. 4.1. ALR-2 The ALR-2 inhibitor set contains six structurally diverse molecules with wide range of activities (Table S1, Fig. S1). [26] In general, the protein–ligand complexes generated with SHAEP produced better results compared to GOLD-generated structures (Table 6). However, the best correlation of 0.94 was obtained with GOLD-generated complex with IGB2 parameters for the single conformation at 4 ps. Similar correlations were obtained for the complexes generated with SHAEP with slightly longer MD simulations (Table 6). Generally the correlations of the binding affinities for the inhibitors were very good with IGB1, IGB2, MM-PBSA, and SIE with complexes generated with SHAEP. Interestingly, in each studied case the results obtained from the structures generated with GOLD significantly worsened during the 512 ps simulation. With the complexes made with SHAEP, with the exception of MMPBSA, which stays really stable throughout the simulation, the correlations increase at first until 64–128 ps and then drop back close to the starting values. It is also notable that while IGB5 performed well in the identification of active ligands, its performance to predict the relative ligand-binding affinities for the selected ligands was not successful. Also, it is surprising that the starting conformations produced with SHAEP generated very good correlation, while the identification of active ligands from inactive molecules was less effective to that of GOLD. 4.2. AmpC The molecular set used as AmpC inhibitors contains 26 cephalosporin-derived structures with wide range of activities (Table S1, Fig. S1) [58]. The results show that the correlation coefficients are negative in almost every studied case (Table 7). Only the starting geometries generated with GOLD and the free energies of binding calculated at the first time point at 4 ps could produce a positive correlation between the experimental and predicted affinities. Still, the correlations were very poor, 0.31 at best (IGB1 in Table 7). The result could imply that during the MD simulation the protein conformation changes to enable favorable interactions with the weaker binding compounds. 4.3. 6HSP90 The HSP90 inhibitor set contains 14 molecules that have IC50 values within a relatively narrow range (Table S1, Fig. S1) [59]. The results show that the best correlations are obtained with MM-PBSA and with the starting geometries generated with GOLD, however, the correlation was still only 0.21 at best (Table 8). In contrast to the VS studies for HSP90, where the IGB5 performed the best with the structures generated with SHAEP, here the same approach produced clearly the poorest results.

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

311

Table 4 AUC values and enrichment factors for PDE-5 with MD simulations of different lengths. The enrichment factors were calculated at 1%, 5%, and 10% of the ranked database. The number of found hits is in parenthesis after the enrichment factor. The results were calculated for structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps. The best result(s) in each row is/are highlighted in bold. PDE-5

4 ps 1% 5% 10% 32 ps 1% 5% 10% 64 ps 1% 5% 10% 128 ps 1% 5% 10% 256 ps 1% 5% 10% 384 ps 1% 5% 10% 512 ps 1% 5% 10%

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.70 ± 0.04 5.8 (3) 4.7 (12) 3.5 (18) 0.72 ± 0.04 3.8 (2) 5.1 (13) 3.7 (19) 0.73 ± 0.04 3.8 (2) 5.1 (13) 3.5 (18) 0.75 ± 0.03 5.8 (3) 4.3 (11) 3.9 (20) 0.73 ± 0.04 5.8 (3) 3.9 (10) 3.5 (18) 0.72 ± 0.04 5.8 (3) 3.5 (9) 3.3 (17) 0.73 ± 0.04 7.7 (4) 3.9 (10) 3.5 (18)

0.72 ± 0.04 19.2 (10) 5.5 (14) 3.3 (17) 0.72 ± 0.04 15.4 (8) 5.5 (14) 4.1 (21) 0.73 ± 0.04 17.3 (9) 5.9 (15) 3.9 (20) 0.73 ± 0.04 17.3 (9) 5.5 (14) 3.7 (19) 0.73 ± 0.04 17.3 (9) 5.9 (15) 3.5 (18) 0.72 ± 0.04 17.3 (9) 6.3 (16) 3.7 (19) 0.71 ± 0.04 15.4 (8) 5.9 (15) 3.5 (18)

0.68 ± 0.04 5.8 (3) 3.9 (10) 3.1 (16) 0.70 ± 0.04 7.7 (4) 3.1 (8) 3.5 (18) 0.70 ± 0.04 7.7 (4) 2.8 (7) 3.1 (16) 0.70 ± 0.04 7.7 (4) 3.5 (9) 2.9 (15) 0.70 ± 0.04 7.7 (4) 3.5 (9) 3.1 (16) 0.70 ± 0.04 7.7 (4) 3.5 (9) 2.9 (15) 0.70 ± 0.04 7.7 (4) 3.1 (8) 2.6 (13)

0.71 ± 0.04 7.7 (4) 4.3 (11) 2.9 (15) 0.72 ± 0.04 11.5 (6) 3.1 (8) 3.1 (16) 0.73 ± 0.04 11.5 (6) 4.7 (12) 3.1 (16) 0.73 ± 0.04 11.5 (6) 4.3 (11) 3.3 (17) 0.72 ± 0.04 11.5 (6) 3.9 (10) 3.1 (16) 0.71 ± 0.04 11.5 (6) 4.3 (11) 3.3 (17) 0.70 ± 0.04 13.4 (7) 4.3 (11) 3.3 (17)

0.53 ± 0.04 0 0.4 (1) 0.8 (4) 0.52 ± 0.04 0 0.8 (2) 1.0 (5) 0.51 ± 0.04 0 1.2 (3) 1.0 (5) 0.52 ± 0.04 0 0.8 (2) 1.0 (5) 0.52 ± 0.04 0 0.8 (2) 0.8 (4) 0.52 ± 0.04 0 0.8 (2) 0.6 (3) 0.52 ± 0.04 0 0.8 (2) 0.8 (4)

0.48 ± 0.04 1.9 (1) 0.4 (1) 0.4 (2) 0.48 ± 0.04 1.9 (1) 0.8 (2) 0.8 (4) 0.48 ± 0.04 1.9 (1) 0.8 (2) 0.6 (3) 0.47 ± 0.04 0 0.8 (2) 0.4 (2) 0.47 ± 0.04 0 0.8 (2) 0.4 (2) 0.47 ± 0.04 1.9 (1) 0.8 (2) 0.4 (2) 0.47 ± 0.04 0 0.8 (2) 0.4 (2)

0.70 ± 0.03 3.8 (2) 2.0 (5) 2.2 (11) 0.73 ± 0.03 1.9 (1) 2.0 (5) 1.6 (8) 0.72 ± 0.03 1.9 (1) 1.6 (4) 1.6 (8) 0.73 ± 0.04 3.8 (2) 2.0 (5) 1.8 (9) 0.73 ± 0.03 3.8 (2) 2.0 (5) 1.6 (8) 0.73 ± 0.03 3.8 (2) 2.0 (5) 1.6 (8) 0.73 ± 0.03 3.8 (2) 2.0 (5) 1.6 (8)

0.60 ± 0.04 1.9 (1) 2.0 (5) 1.8 (9) 0.60 ± 0.04 0 2.0 (5) 1.6 (8) 0.60 ± 0.04 5.8 (3) 2.0 (5) 1.4 (7) 0.61 ± 0.04 5.8 (3) 1.6 (4) 1.6 (8) 0.61 ± 0.04 5.8 (3) 1.6 (4) 1.2 (6) 0.61 ± 0.04 5.8 (3) 1.6 (4) 1.4 (7) 0.60 ± 0.04 5.8 (3) 1.6 (4) 1.4 (7)

0.67 ± 0.04 3.8 (2) 2.8 (7) 2.4 (12) 0.70 ± 0.04 3.8 (2) 2.0 (5) 2.2 (11) 0.69 ± 0.04 1.9 (1) 2.4 (6) 2.2 (11) 0.69 ± 0.04 3.8 (2) 2.4 (6) 2.6 (13) 0.70 ± 0.04 3.8 (2) 2.4 (6) 2.4 (12) 0.70 ± 0.04 3.8 (2) 1.6 (4) 2.2 (11) 0.70 ± 0.04 3.8 (2) 1.6 (4) 2.4 (12)

0.64 ± 0.04 9.6 (5) 3.1 (8) 2.4 (12) 0.61 ± 0.04 7.7 (4) 3.5 (9) 2.2 (12) 0.62 ± 0.04 9.6 (5) 3.5 (9) 2.6 (13) 0.62 ± 0.04 9.6 (5) 3.5 (9) 2.2 (11) 0.63 ± 0.04 9.6 (5) 3.1 (8) 2.0 (10) 0.63 ± 0.04 5.8 (3) 1.6 (4) 1.4 (7) 0.60 ± 0.04 5.8 (3) 1.6 (4) 1.4 (7)

Table 5 AUC values and enrichment factors for PR with MD simulations of different lengths. The enrichment factors were calculated at 1%, 5%, and 10% of the ranked database. The number of found hits is in parentheses after the enrichment factor. The results were calculated for structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps. The best result(s) in each row is/are highlighted in bold. PR

4 ps 1% 5% 10% 32 ps 1% 5% 10% 64 ps 1% 5% 10% 128 ps 1% 5% 10% 256 ps 1% 5% 10% 384 ps 1% 5% 10% 512 ps 1% 5% 10%

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.51 ± 0.08 11.0 (3) 5.2 (7) 3.0 (8) 0.53 ± 0.08 11.0 (3) 4.4 (6) 3.0 (8) 0.53 ± 0.08 11.0 (3) 5.2 (7) 3.0 (8) 0.53 ± 0.08 11.0 (3) 5.2 (7) 3.0 (8) 0.53 ± 0.08 11.0 (3) 4.4 (6) 3.0 (8) 0.52 ± 0.08 11.0 (3) 4.4 (6) 3.0 (8) 0.52 ± 0.08 11.0 (3) 4.4 (6) 3.0 (8)

0.50 ± 0.06 7.4 (2) 3.7 (5) 1.9 (5) 0.47 ± 0.07 7.4 (2) 3.7 (5) 2.2 (6) 0.48 ± 0.07 7.4 (2) 3.7 (5) 2.2 (6) 0.49 ± 0.07 7.4 (2) 3.7 (5) 2.2 (6) 0.49 ± 0.07 7.4 (2) 3.7 (5) 2.2 (6) 0.50 ± 0.07 7.4 (2) 3.7 (5) 2.2 (6) 0.50 ± 0.07 7.4 (2) 3.7 (5) 2.2 (6)

0.53 ± 0.08 14.7 (4) 5.9 (8) 3.4 (9) 0.55 ± 0.08 14.7 (4) 5.2 (7) 3.4 (9) 0.56 ± 0.08 14.7 (4) 5.9 (8) 3.4 (9) 0.55 ± 0.08 18.4 (5) 5.9 (8) 3.0 (8) 0.55 ± 0.08 14.7 (4) 5.9 (8) 3.0 (8) 0.55 ± 0.08 18.4 (5) 5.9 (8) 3.0 (8) 0.54 ± 0.08 18.4 (5) 5.9 (8) 3.0 (8)

0.55 ± 0.07 14.7 (4) 3.7 (5) 2.6 (7) 0.52 ± 0.07 18.4 (5) 4.4 (6) 3.0 (8) 0.54 ± 0.07 14.7 (4) 5.2 (7) 3.0 (8) 0.55 ± 0.07 14.7 (4) 4.4 (6) 3.4 (9) 0.54 ± 0.07 14.7 (4) 4.4 (6) 3.4 (9) 0.55 ± 0.07 14.7 (4) 5.2 (7) 3.0 (8) 0.55 ± 0.07 11.0 (3) 5.2 (7) 2.6 (7)

0.72 ± 0.04 7.4 (2) 3.0 (4) 1.9 (5) 0.71 ± 0.04 3.7 (1) 3.7 (5) 2.6 (7) 0.72 ± 0.04 3.7 (1) 3.7 (5) 2.6 (7) 0.71 ± 0.05 3.7 (1) 3.7 (5) 2.6 (7) 0.71 ± 0.05 7.4 (2) 3.7 (5) 2.6 (7) 0.72 ± 0.04 7.4 (2) 3.0 (4) 2.6 (7) 0.72 ± 0.04 7.4 (2) 3.7 (5) 2.6 (7)

0.76 ± 0.05 0 3.7 (5) 3.4 (9) 0.78 ± 0.05 7.4 (2) 4.4 (6) 3.7 (10) 0.78 ± 0.05 7.4 (2) 3.0 (4) 4.5 (12) 0.77±0.05 7.4 (2) 3.0 (4) 4.5 (12) 0.77 ± 0.05 7.4 (2) 3.7 (5) 4.1 (11) 0.77 ± 0.05 7.4 (2) 4.4 (6) 3.7 (10) 0.77 ± 0.05 7.4 (2) 3.0 (4) 3.4 (9)

0.69 ± 0.05 11.0 (3) 4.4 (6) 3.4 (9) 0.69 ± 0.06 18.4 (5) 7.4 (10) 4.1 (11) 0.69 ± 0.05 18.4 (5) 7.4 (10) 4.5 (12) 0.68 ± 0.06 18.4 (5) 7.4 (10) 4.5 (12) 0.67 ± 0.07 18.4 (5) 6.6 (9) 4.5 (12) 0.66 ± 0.07 18.4 (5) 6.6 (9) 4.5 (12) 0.66 ± 0.07 18.4 (5) 5.9 (8) 4.5 (12)

0.64 ± 0.05 11.0 (3) 4.4 (6) 2.6 (7) 0.64 ± 0.06 14.7 (4) 4.4 (6) 3.4 (9) 0.65 ± 0.06 14.7 (4) 5.2 (7) 3.7 (10) 0.64 ± 0.06 14.7 (4) 6.6 (9) 3.7 (10) 0.63 ± 0.06 11.0 (3) 5.2 (7) 3.7 (10) 0.64 ± 0.06 14.7 (4) 5.2 (7) 3.7 (10) 0.65 ± 0.06 14.7 (4) 5.2 (7) 3.7 (10)

0.53 ± 0.08 14.7 (4) 5.9 (8) 3.4 (9) 0.53 ± 0.08 22.1 (6) 5.9 (8) 3.0 (8) 0.54 ± 0.08 22.1 (6) 5.9 (8) 3.4 (9) 0.54 ± 0.08 18.4 (5) 5.9 (8) 3.4 (9) 0.54 ± 0.08 18.4 (5) 5.9 (8) 3.4 (9) 0.54 ± 0.08 18.4 (5) 5.9 (8) 3.4 (9) 0.53 ± 0.08 18.4 (5) 5.9 (8) 3.0 (8)

0.51 ± 0.07 7.4 (2) 3.7 (5) 2.6 (7) 0.49 ± 0.07 14.7 (4) 3.7 (5) 2.2 (6) 0.48 ± 0.07 11.0 (3) 3.7 (5) 2.6 (7) 0.49 ± 0.07 14.7 (4) 4.4 (6) 3.0 (8) 0.48 ± 0.07 14.7 (4) 4.4 (6) 3.0 (8) 0.49 ± 0.07 14.7 (4) 4.4 (6) 3.0 (8) 0.50 ± 0.07 14.7 (4) 4.4 (6) 3.0 (8)

4.4. PDE-5 The molecular dataset for PDE-5 contains tadalafil and eight analogs (Table S1, Fig. S1). [60] It is an interesting set because the

IC50 values are available for several individual stereoisomers of the analogs. The total number of different structures in the set is 22. Considering that the molecules are analogous and contain different stereoisomers, it is probably not surprising that the correlations

312

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

Table 6 Correlation coefficients (r) for the pIC50 values versus the predicted binding affinities for ALR-2 with structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps of the MD simulation. The best result in each row is highlighted in bold. ALR-2

4 ps 32 ps 64 ps 128 ps 256 ps 384 ps 512 ps

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.93 0.85 0.71 0.60 0.53 0.51 0.47

0.80 0.88 0.89 0.89 0.87 0.85 0.81

0.94 0.84 0.71 0.60 0.53 0.48 0.44

0.82 0.92 0.93 0.92 0.89 0.87 0.84

−0.14 −0.22 −0.17 −0.16 −0.17 −0.24 −0.28

0.49 0.56 0.58 0.55 0.53 0.53 0.53

0.73 0.76 0.71 0.61 0.37 0.31 0.26

0.73 0.71 0.72 0.72 0.72 0.72 0.72

0.41 0.12 0.09 −0.01 −0.05 −0.01 −0.03

0.74 0.76 0.77 0.77 0.75 0.71 0.64

Fig. 3. The effect of the MD simulation length on the VS results in terms of the ROC AUC values. Solid lines depict the results obtained for GOLD-generated structures and the dashed lines depict the results for SHAEP-generated structures.

Table 7 Correlation coefficients (r) for the pIC50 values versus the predicted binding affinities for AmpC with structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps of the MD simulation. The best result in each row is highlighted in bold. AmpC

4 ps 32 ps 64 ps 128 ps 256 ps 384 ps 512 ps

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.31 −0.09 −0.08 −0.06 −0.08 −0.04 −0.02

−0.22 −0.26 −0.22 −0.28 −0.37 −0.39 −0.39

0.24 −0.17 −0.16 −0.13 −0.14 −0.10 −0.08

−0.23 −0.26 −0.23 −0.28 −0.36 −0.38 −0.39

0.08 −0.26 −0.30 −0.25 −0.20 −0.16 −0.14

−0.23 −0.25 −0.25 −0.27 −0.31 −0.32 −0.34

0.13 −0.25 −0.29 −0.28 −0.28 −0.24 −0.22

−0.36 −0.35 −0.34 −0.37 −0.48 −0.49 −0.49

0.27 −0.17 −0.19 −0.16 −0.15 −0.12 −0.10

−0.26 −0.27 −0.25 −0.31 −0.40 −0.41 −0.42

Table 8 Correlation coefficients (r) for the pIC50 values versus the predicted binding affinities for HSP90 with structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps of the MD simulation. The best result in each row is highlighted in bold. HSP90

4 ps 32 ps 64 ps 128 ps 256 ps 384 ps 512 ps

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

−0.09 −0.02 0.01 0.02 −0.03 −0.06 −0.06

0.13 0.12 0.09 0.08 0.01 -0.02 -0.05

−0.10 −0.01 0.04 0.05 0.00 −0.04 −0.05

−0.08 −0.09 −0.07 −0.07 −0.01 0.02 0.05

0.09 −0.12 −0.11 −0.12 −0.15 −0.17 −0.17

−0.42 −0.35 −0.38 −0.41 −0.41 −0.42 −0.41

0.18 0.16 0.15 0.21 0.04 −0.08 −0.13

0.02 −0.06 −0.05 −0.11 −0.04 −0.01 0.00

−0.16 −0.10 −0.09 −0.12 −0.13 −0.16 −0.18

−0.07 −0.09 −0.09 −0.07 −0.02 0.01 0.03

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

313

Table 9 Correlation coefficients (r) for the pIC50 values versus the predicted binding affinities for PDE-5 with structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps of the MD simulation. The best result in each row is highlighted in bold. PDE-5

4 ps 32 ps 64 ps 128 ps 256 ps 384 ps 512 ps

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.26 0.12 0.15 0.24 0.28 0.31 0.34

−0.22 0 0.10 0.17 0.18 0.17 0.18

0.23 0.07 0.09 0.2 0.23 0.26 0.29

−0.21 0.01 0.13 0.19 0.21 0.21 0.22

0.22 0.27 0.26 0.26 0.29 0.28 0.27

0.35 0.31 0.32 0.32 0.35 0.33 0.33

−0.01 0.01 0.11 0.15 0.14 0.17 0.20

−0.17 0.03 0.08 0.15 0.01 −0.05 −0.05

0.16 −0.18 −0.12 0.03 0.12 0.17 0.25

−0.17 0.07 0.11 0.18 0.16 0.19 0.22

Table 10 The correlation coefficients (r) for the PDE-5 inhibitor stereoisomers 4a–d, 7a–d, and 8a–d calculated for the structures obtained from the 512 ps MD simulation. The best result in each row is highlighted in bold. IGB1

4a–d 7a–d 8a–d

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.75 −0.22 0.55

0.59 0.08 0.67

0.76 −0.62 0.61

0.61 0.12 0.76

0.61 −0.64 0.94

0.73 −0.07 0.80

0.28 −0.06 0.21

0.70 −0.03 0.43

0.55 −0.10 0.33

0.67 0.19 0.59

ity. For isomers 8a–d the correlation was very high (0.94), and again the rank-order was correct for 3 of the 4 isomers. However, the two least active compounds, for which the order of binding was incorrectly predicted, have the measured as well as predicted activities very close to each other. As is obvious from the plot, the binding free energies could not be predicted correctly for isomers 7a–d, and the correlation is negative. 4.5. PR

Fig. 4. The predicted G values by IGB5 versus the experimental pIC50 values for the GOLD generated structures for the PDE-5 inhibitor stereoisomers of 4a–d, 7a–d, and 8a–d. The structures were obtained from the 512 ps MD simulation.

produced with the different methods are weak (Table 9). The most efficient in predicting the binding affinities was IGB5 with complexes generated with SHAEP. However, the best correlation was still only 0.35. Otherwise, in general, the complexes produced with GOLD ended up with slightly better correlations, and the MM-GBSA methods outperformed MM-PBSA and SIE. Also, the correlations seem to improve with longer MD simulations. Even though the overall correlations were very low for the whole set of PDE-5 inhibitors, a closer examination of the results shows that in fact the correlations for the different stereoisomers are rather good in some cases, when the molecules with four different stereoisomers were considered (Fig. 4, Table 10). Especially for the stereoisomers 4a–d and 8a–d (Table S1, Fig. S1), the results show relatively good correlations (0.76 and 0.94 at best, for 4a–d and 8a–d, respectively). For the stereoisomers 7a–d (Table S1, Fig. S1), however, the correlations were very low, or even negative. This is somewhat surprising since the only difference between 7a–d and 8a–d is that a methyl group in 7a–d is replaced by an ethyl group in 8a–d (Fig. S1). An example of the binding free energy predictions is shown in Fig. 4, which shows a plot of predicted G values from IGB5 with the GOLD-generated structures and the measured pIC50 values. The IGB5 could correctly rank-order 3 out 4 of the isomers 4a–d, with one outlier, which was predicted to be less active compared to the actual measured activ-

The molecular set for PR contains 12 molecules; mifepristone used as a reference, and 11 nonsteroidal chromene-based antagonists (Table S1, Fig. S1). [61] For this dataset SIE method produced the highest correlation of 0.83 with complexes generated with GOLD and the long 512 ps simulation (Table 11). In addition, the IGB2 parameters produced a reasonably good correlation when the GOLD-generated complexes were studied (0.71). In contrast to the VS studies with the DUD database, where IGB5 and MM-PBSA were among the best methods according to AUC values and enrichment factors, here these methods produced only weak correlations for the binding affinities (0.28 and 0.44, respectively). In fact, for the IGB5 the correlations were mostly negative. An interesting result is that for the complexes generated with GOLD the correlations generally increase during the MD simulation, whereas for complexes generated with SHAEP the correlations remain similar or decrease. However, for IGB5 the situation is the opposite. We further investigated whether the highly active molecules used in binding affinity predictions could be identified among the DUD molecules for the studied targets. For this, the results of the binding free energy calculations for the molecules with pIC50 ≥ 6 (Table S1) were combined to the results of DUD VS, and the number of these active molecules found in the top 1%, 5%, and 10% of the databases was determined (Table 12). For ALR-2 there is only one molecule with pIC50 ≥ 6 (Table S1). This molecule, tolrestat, was identified within the top 1% with only IGB5 and MM-PBSA with the SHAEP generated structures (Table 12). With SIE the active molecule was found within the top 5% with the GOLD created structures, and generally for the other methods it was found in the top 10%. The relatively poor performance overall for the free energy calculation methods is somewhat surprising, since the protein crystal structure, which was used in the studies, was originally bound to the same inhibitor, and the starting conformations with both GOLD and SHAEP were

314

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

Table 11 Correlation coefficients (r) for the pIC50 values versus the predicted binding affinities for PR with structures obtained at 4 ps, 32 ps, 64 ps, 128 ps, 256 ps, 384 ps, and 512 ps of the MD simulation. The best result in each row is highlighted in bold. PR

4 ps 32 ps 64 ps 128 ps 256 ps 384 ps 512 ps

IGB1

IGB2

IGB5

PB

SIE

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

0.48 0.57 0.57 0.57 0.58 0.62 0.63

0.54 0.44 0.45 0.47 0.46 0.47 0.48

0.53 0.63 0.64 0.66 0.68 0.70 0.71

0.65 0.56 0.57 0.58 0.56 0.57 0.58

−0.36 −0.35 −0.40 −0.47 −0.61 −0.60 −0.57

0.04 −0.15 0.01 −0.04 0.06 0.19 0.28

0.27 0.36 0.33 0.37 0.37 0.38 0.37

0.44 0.31 0.31 0.34 0.35 0.34 0.33

0.61 0.67 0.77 0.79 0.81 0.82 0.83

0.61 0.53 0.58 0.60 0.61 0.62 0.63

Table 12 The number of molecules with pIC50 ≥ 6 found in the top 1%, 5%, and 10% of the ranked list when the results were combined with the DUD VS results. The results were taken from the 512 ps simulations. The number in the parentheses implies the number of molecules with pIC50 ≥ 6 in the given dataset. IGB1

IGB2

IGB5

PB

SIE

ALR-2 (1)

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

GOLD

SHAEP

1% 5% 10 %

0 0 1

0 0 1

0 0 0

0 0 1

0 0 1

1 1 1

0 1 1

1 1 1

0 1 1

0 1 1

AmpC (13) 1% 5% 10 %

2 5 5

1 4 6

0 5 5

0 2 4

0 3 5

0 2 4

1 2 2

0 2 3

1 4 6

2 5 6

HSP90 (11) 1% 5% 10 %

0 0 0

0 0 0

0 0 1

0 0 0

1 1 1

1 1 1

0 0 0

0 0 0

0 0 0

0 0 0

PDE-5 (11) 1% 5% 10 %

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

0 0 3

0 1 1

0 0 2

0 0 0

PR (8) 1% 5% 10 %

1 5 7

1 3 5

1 4 5

1 2 4

2 8 8

5 8 8

1 5 5

1 3 5

4 8 8

2 6 7

very close to the correct binding conformation with RMSD values 1.90 and 0.48, respectively. It is also notable that the IGB1 and IGB2 methods, which were the most successful in rank-ordering the active molecules according to their binding affinities, were the least successful of the methods in identifying tolrestat from the DUD molecules. For AmpC there are 13 molecules with pIC50 ≥ 6 (Table S1). The results show relatively low efficiency in identifying the active molecules among the DUD molecules (Table 12). IGB1 and SIE with both GOLD and SHAEP generated structures, and MM-PBSA with GOLD generated structure find one to two molecules within the top 1%, whereas IGB2, IGB5, and MM-PBSA with SHAEP generated structure do not find any molecules within the top 1%. Within the top 5% IGB1, SIE, and IGB2 with GOLD generated structure identify four to five active molecules among the DUD molecules. No significant change is observed when the examination portion is increased from 5% to 10%. Within the top 10% the best result is achieved with IGB1 with SHAEP generated structure and with SIE both GOLD and SHAEP generated structures, six molecules found. For HSP90 there are 11 molecules with pIC50 ≥ 6 (Table S1). However, only one molecule is found within the top 1% with both GOLD and SHAEP generated structures with IGB5 (Table 12). Additionally, one molecule was found within the top 10% of GOLDgenerated structures with IGB2. Overall, recognition of highly active HSP90 molecules from the DUD molecules was really poor, which is in line with binding affinity predictions (Table 8). For PDE-5 there are a total of 11 molecules with pIC50 ≥ 6 (Table S1). The results show very low efficiency in identifying these active

molecules among the DUD molecules (Table 12). The MM-GBSA methods could not identify any of these molecules within the top 10%. Further, MM-PBSA could correctly identify only one tadalafil analog in the top 5% with the structures generated with SHAEP, and MM-PBSA and SIE could only find 3 and 2 molecules, respectively, within the top 10% of GOLD-generated structures. Tadalafil, which is the most potent inhibitor among the CHEMBL dataset used here, was found within the top 10% only with MM-PBSA and with the GOLD-generated structures (data not shown). The binding modes predicted by GOLD and SHAEP were not entirely optimal with RMSD values 3.02 and 2.43 for the GOLD- and SHAEP-generated structures, respectively, when the binding mode is compared to the PDE-5 structure with bound tadalafil (pdb: 1xoz) (Fig. 5). To test whether the poor performance of the MM-GB(PB)SA and SIE methods in ranking tadalafil is caused by the poor binding mode prediction, the free energy calculations were also performed for the tadalafil extracted from the structure 1xoz. The results for tadalafil did not improve, which implies that the poor binding affinity prediction is presumably caused by, e.g., the incorrect parameters in free energy calculations. Unfortunately, to our knowledge, there are no experimentally determined 3D structures available for the tadalafil analogs bound to PDE-5, so we were not able to test whether the results could have improved with a better binding conformation to the target protein. From the molecular dataset for PR, eight structures were included in this study (Table S1). The results show that IGB5 with both tested complex generation methods and SIE with the GOLDgenerated protein-ligand complexes could identify all of the active

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

Fig. 5. The binding mode of tadalafil predicted by GOLD (orange carbon atoms) and SHAEP (purple carbon atoms). The tadalafil from the crystal structure 1xoz is shown for comparison (green carbon atoms). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

molecules within the top 5% (Table 12). Also the number of identified hits at 1% was higher for these methods compared to IGB1, IGB2, and MM-PBSA. The most potent inhibitor, mifepristone, was found within the top 1% in every studied case (data not shown). It is noteworthy that even though IGB5 did not succeed in determining the correct order of molecules according to binding affinities (Table 11), the highly active molecules were nevertheless found among the top 5% of the ranked list of DUD database. 5. Discussion Earlier studies have indicated that MM-GB(PB)SA methods can be useful and accurate tools for binding affinity predictions. MMPBSA has been used to predict the binding affinities and binding modes of HIV-RT inhibitors with high accuracy [62]. MM-GBSA has been tested for various pharmaceutically important protein targets to study the ability of the method to calculate binding free energies, and sometimes the results have been in remarkable correlation with the experimental results [63]. The SIE method was studied in predicting binding free energies for phosphodiesterase 2 (PDE-2) inhibitors, and the results showed excellent correlation with experimental binding data [15]. In spite of these promising results, it has also been shown that the results with MM-GB(PB)SA methods can be highly case-specific and depend heavily on the parameters used to calculate the free energies and the methods by which the geometries of the protein–ligand complexes have been derived [13,64]. In this study we compared the performance of MM-PBSA, MMGBSA, and SIE in their ability to separate active molecules from inactive and the ability to predict experimentally determined binding affinities. Our results showed significant differences for the results obtained with these free energy calculation methods. The VS enrichment studies showed that for ALR-2 the MM-GBSA with IGB5 parameters produced the highest AUC values and enrichment factors, followed by MM-PBSA and SIE. Also for AmpC MM-GBSA with IGB5 achieved the best results according to AUC values and enrichment factors. According to the AUC values, MM-GBSA with

315

IGB5 is clearly the best option for HSP-90, however, the enrichment factors are poor. The best enrichments are achieved with IGB1 and SIE. For PDE-5, however, the IGB1 was the best considering both the AUC values and the early enrichment. The overall enrichment for PR according to the AUC values was clearly the best with IGB5, but when considering the early enrichment the best method was MMPBSA. When evaluating the methods by the accuracy of the binding affinity predictions the IGB2 was the most efficient for ALR-2, IGB1 for AmpC, MM-PBSA for HSP-90, IGB5 for PDE-5, and SIE for PR. So, for every test case a different method was the best in the accuracy of the binding affinity predictions. In theory, MM-PBSA is more rigorous than MM-GBSA. Thus, in principle, it is expected to perform better than the MM-GBSA methods studied here. In the VS studies MM-PBSA performed relatively well producing similar or better AUC values and enrichment factors for all the targets when compared to the other methods. However, considering that MM-PBSA requires more computational effort than the other tested methods, it might not be suitable for screening of large compound databases. Even though the performance of MM-PBSA was relatively good in VS compared to the other methods, in the binding free energy estimation studies MM-PBSA did not succeed as well. Although in some studies MM-PBSA has given more accurate results compared to MM-GBSA [12], always the increased computational effort of MM-PBSA does not yield better results than MM-GBSA [13]. The GB models make approximations to the PB equation, and thus they are also more efficient in computational terms than PB. The two GB models created by Onufriev et al. (IGB2 and IGB5) were developed to better predict the energies of buried residues for both small molecules and macromolecules as compared to IGB1 [8]. Even though the theory behind the IGB2 and IGB5 is the same, the results show significant differences in the VS performance and in their ability to rank-order molecules with known affinities. Further, IGB1 and IGB2 gave more similar results compared to each other. IGB5 performed well in producing the highest AUC values in VS for ALR-2, AmpC, HSP-90, and PR. This is in contrast with its ability to predict the binding affinities for the same targets. The most striking example is ALR-2, for which the other methods could produce reasonably good or excellent correlations compared to experimental data, but the results with IGB5 produced either weak or negative correlation. Earlier studies have also shown large variations between the MM-GB(PB)SA methods when different protein systems are used for evaluation. According to Hou and colleagues, IGB2 was the most accurate compared to other GB models and MM-PBSA [13]. SIE is similar to MM-PBSA in that it uses the Poisson–Boltzmann equation for the calculation of the polar solvation energy. However, the SIE method employs a boundary element method (BEM) for solving the Poisson–Boltzmann equation, which is computationally more efficient than the finite-difference method used in MM-PBSA [11]. Accordingly, SIE is computationally more efficient compared to the MM-PBSA method in free energy calculations, but at the same time less efficient than the GB methods. A further difference between MM-PBSA and SIE is that the latter contains terms that have been fitted using a set of protein–ligand complexes with known affinities and crystal structures. Perhaps owing to the similar method in calculating the polar solvation in both SIE and MM-PBSA, the VS and binding affinity prediction results show similar enrichments and correlations for ALR-2, AmpC, HSP-90, and PDE-5. However, for PR more differences appeared in the results. In VS for PR, the MM-PBSA method performed better compared to SIE in terms of AUC values. However, when the early enrichments were considered, SIE produced similar results for the most part as IGB5. In the binding affinity predictions for PR, SIE was clearly better compared to MM-PBSA with both methods used in complex creation. A major difference between SIE and MM-GB(PB)SA methods is that the predicted G values obtained with SIE are closer in magnitude

316

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

to the experimental G values, whereas for MM-GB(PB)SA the G values are often overestimated [15,65]. The MM-GB(PB)SA and SIE methods are comparatively costly when their computational efficiency is compared to traditional VS methods, such as docking and ligand-based methods. The majority of the computational load comes from the MD simulation, which is traditionally used to generate different geometries for the protein–ligand complexes. Therefore, if a shorter simulation could produce accurate results, it would enable screening of larger molecular databases with these methods. Earlier studies have shown that short MD simulations can be adequate in determining binding affinities, and even a single minimized snapshot can be used in free energy calculations and in VS [12,18,20]. Hou et al. have studied the effect of solute dielectric constant and the length of MD simulation on the binding affinity estimation by MM-PBSA for six different protein targets, and they found that a longer MD simulation is not always necessary for the accuracy of results [20]. This is supported by a study by Sun et al., where they tested the effect of solute dielectric constant and the length of MD simulation for the VS of three different tyrosine kinases, and they concluded that an MD simulation is not necessarily needed for efficient virtual screening, but a minimized docking structure can be sufficient [18]. Our results show that the length of the MD simulation did not affect much the VS efficiency in terms of AUC values and enrichment factors. In fact, often the highest values were obtained with relatively short MD simulations of less than 128 ps. Additionally, in many cases even a single snapshot from the beginning of the MD simulation could produce similar and sometimes better results compared to longer simulations. This is an encouraging result for VS applications of MM-GB(PB)SA and SIE methods. The length of the MD simulation had a more pronounced effect on the prediction of the experimentally determined binding affinities, and the results were case-specific. For ALR-2, shorter simulations of fewer than 128 ps produced better results, whereas the prediction of PDE-5 inhibitor binding free energies benefitted from long 512 ps simulations. For PR the results were a bit more complicated, because they depended on the method by which the protein-ligand complexes were produced; in general, the structures created with GOLD produced better correlation with longer simulations, and structures generated with SHAEP produced the highest correlation with the single conformation at 4 ps. The reason for this is probably the initial conformation of the ligand in relation to the protein. The results suggest that with SHAEP, the starting conformation for MD simulation was better compared to GOLD. However, in contrast to the SHAEP-generated complexes, the structures obtained by docking benefitted from the conformational search accomplished by MD and eventually produced better results. The method by which the protein-ligand complex is generated and how the geometries for free energy calculations have been derived can also influence the outcome. Here we studied two different methods for the complex generation, which are both important and widely used in VS, namely molecular docking and ligand-based similarity comparison. The results show that in many cases the results obtained with the different approaches are quite similar to each other. However, there are also differences. For AmpC all SHAEP generated starting complexes are clearly better compared to GOLD generated starting complexes according to AUC values and enrichment factors. Also, for PDE-5 in the VS studies with IGB1 and IGB2 the structures generated with SHAEP clearly produced better early enrichment compared to structures generated with GOLD. Overall the results show, however, that the ligand-based superimposition can be used in place of docking as a fast and robust way to generate protein–ligand complexes for post-processing with MM-GB(PB)SA or SIE.

The advantage of MM-GB(PB)SA and SIE over docking is that the protein–ligand complex is relaxed with energy-minimization and MD simulation, whereas docking is usually done to only one static protein structure. The protein targets studied here are known to undergo ligand-dependent conformational changes upon binding. For AmpC and PR the MM-GB(PB)SA and SIE methods were clearly better compared to GOLD in VS. According to the AUC values also PDE-5 benefitted slightly from the post-processing, however, when enrichment factors are scrutinized the improvement is clear. For HSP-90 the enrichment factors improved with all MM-GB(PB)SA and SIE compared to docking, and furthermore, the AUC values considerably improved with IGB5. However, for ALR-2 the postprocessing results did not show any significant improvement over docking. Thus, it appears that also the advantage of the conformational search to the binding energy calculations is case-specific. When considering the affinity predictions with the used methods, there can be many issues that affect the reliability of the results. For example, our main aim was to understand the usability of the implicit solvent models in the VS campaign. Accordingly, we chose to use single starting conformations obtained separately with docking and ligand-based superimpositioning. Obviously, the used conformation can be different from biologically active conformation. Furthermore, to minimize the computational cost in VS campaign, it is quite often necessary to rely on a single simulation, which should be as short as possible. As shown here and elsewhere [18,20], the length of the MD simulation is often irrelevant for the reliability of binding affinity predictions. In addition, tautomerization of small molecules may prevent the identification of biologically active conformation, which again leads into incorrect binding affinity predictions. For example, the method that we have employed for generation of ligand structures does not identify the correct protonation for PDE-5 inhibitor sildenafil [57]. Taken together, the choice between the usage of single binding conformation and MD simulation and usage of multiple starting conformations with repeated MD simulations is balancing between efficient VS campaign and statistical significance. 6. Conclusions Our results show that both docking and ligand-based molecular superimposition can be used for the generation of the starting structure for free energy calculations. The length of the MD simulation is not an important factor in VS applications of the methods and often even a single snapshot can be used for the binding energy calculations. However, in the prediction of the actual binding affinities, a longer MD simulation may be required. When predicting binding affinities, more accurate results can be obtained for ligands that have a wide range of activities and that are structurally relatively different. Despite the fact that the ranking of stereoisomers (PDE-5) was not highly successful, the methods could in some cases distinguish between the high-affinity and low-affinity stereoisomers. Overall, we found that with single MD simulations and with the protein targets studied here the results for both the VS and the binding affinity estimations were highly case-specific. The generation of the starting conformation creates another hurdle, because there are several ways to accomplish this, and without structural information, as is the case in VS, it is difficult to determine which method will be the best for a particular target. Thus, our study should remind the users of these techniques to use caution when applying these methods to VS or to binding affinity estimations. Acknowledgements This work was financially supported by FinPharma Doctoral Programme—Drug Discovery (SIV) and National Doctoral Pro-

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

gramme in Nanoscience (SPN). CSC, The Finnish IT Center for Science is acknowledged for computational resources (OTP; project Nos. jyy2516 and jyy2585).

[24] [25]

Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jmgm.2015.10. 012.

[26]

[27]

References [1] G.B. McGaughey, R.P. Sheridan, C.I. Bayly, J.C. Culberson, C. Kreatsoulas, S. Lindsley, et al., Comparison of topological, shape, and docking methods in virtual screening, J. Chem. Inf. Model. 47 (2007) 1504–1519. [2] G.L. Warren, C.W. Andrews, A.M. Capelli, B. Clarke, J. LaLonde, M.H. Lambert, et al., A critical assessment of docking programs and scoring functions, J. Med. Chem. 49 (2006) 5912–5931. [3] M. von Korff, J. Freyss, T. Sander, Comparison of ligand- and structure-based virtual screening on the DUD data set, J. Chem. Inf. Model. 49 (2009) 209–231. [4] S.I. Virtanen, O.T. Pentikainen, Efficient virtual screening using multiple protein conformations described as negative images of the ligand-binding site, J. Chem. Inf. Model. 50 (2010) 1005–1011. [5] T.P. Lybrand, J.A. McCammon, G. Wipff, Theoretical calculation of relative binding affinity in host-guest systems, Proc. Natl. Acad. Sci. U. S. A. 83 (1986) 833–835. [6] P. Kollman, Free-energy calculations—applications to chemical and biochemical phenomena, Chem. Rev. 93 (1993) 2395–2417. [7] P.A. Kollman, I. Massova, C. Reyes, B. Kuhn, S. Huo, L. Chong, et al., Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models, Acc. Chem. Res. 33 (2000) 889–897. [8] A. Onufriev, D. Bashford, D.A. Case, Exploring protein native states and large-scale conformational changes with a modified generalized born model, Proteins 55 (2004) 383–394. [9] V. Tsui, D.A. Case, Theory and applications of the generalized Born solvation model in macromolecular simulations, Biopolymers 56 (2000) 275–291. [10] J.M. Wang, T.J. Hou, X.J. Xu, Recent advances in free energy calculations with a combination of molecular mechanics and continuum models, Curr. Comput.-Aid Drug 2 (2006) 287–306. [11] M. Naim, S. Bhat, K.N. Rankin, S. Dennis, S.F. Chowdhury, I. Siddiqi, et al., Solvated interaction energy (SIE) for scoring protein-ligand binding affinities. 1. Exploring the parameter space, J. Chem. Inf. Model. 47 (2007) 122–133. [12] A.M. Ferrari, G. Degliesposti, M. Sgobba, G. Rastelli, Validation of an automated procedure for the prediction of relative free energies of binding on a set of aldose reductase inhibitors, Bioorg. Med. Chem. 15 (2007) 7865–7877. [13] T. Hou, J. Wang, Y. Li, W. Wang, Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations, J. Chem. Inf. Model. 51 (2011) 69–82. [14] G. Rastelli, A. Del Rio, G. Degliesposti, M. Sgobba, Fast and accurate predictions of binding free energies using MM-PBSA and MM-GBSA, J. Comput. Chem. 31 (2010) 797–810. [15] B. Yang, A. Hamza, G. Chen, Y. Wang, C.G. Zhan, Computational determination of binding structures and free energies of phosphodiesterase-2 with benzo[1,4]diazepin-2-one derivatives, J. Phys. Chem. B 114 (2010) 16020–16028. [16] M. Ylilauri, O.T. Pentikainen, MMGBSA as a tool to understand the binding affinities of filamin-peptide interactions, J. Chem. Inf. Model. 53 (2013) 2626–2633. [17] H.Y. Sun, Y.Y. Li, S. Tian, L. Xu, T.J. Hou, Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set, Phys. Chem. Chem. Phys. 16 (2014) 16719–16729. [18] H.Y. Sun, Y.Y. Li, M.Y. Shen, S. Tian, L. Xu, P.C. Pan, et al., Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring, Phys. Chem. Chem. Phys. 16 (2014) 22035–22045. [19] L. Xu, H.Y. Sun, Y.Y. Li, J.M. Wang, T.J. Hou, Assessing the performance of MM/PBSA and MM/GBSA methods. 3. The Impact of force fields and ligand charge models, J. Phys. Chem. B 117 (2013) 8408–8421. [20] T.J. Hou, J.M. Wang, Y.Y. Li, W. Wang, Assessing the performance of the molecular mechanics/Poisson Boltzmann surface area and molecular mechanics/generalized born surface area methods. II. The accuracy of ranking poses generated from docking, J. Comput. Chem. 32 (2011) 866–877. [21] X. Zhang, S.E. Wong, F.C. Lightstone, Toward fully automated high performance computing drug discovery: a massively parallel virtual screening pipeline for docking and molecular mechanics/generalized Born surface area rescoring to improve enrichment, J. Chem. Inf. Model. 54 (2014) 324–337. [22] D.C. Thompson, C. Humblet, D. Joseph-McCarthy, Investigation of MM-PBSA rescoring of docking poses, J. Chem. Inf. Model. 48 (2008) 1081–1091. [23] S.P. Niinivehmas, S.I. Virtanen, J.V. Lehtonen, P.A. Postila, O.T. Pentikainen, Comparison of virtual high-throughput screening methods for the

[28]

[29]

[30]

[31]

[32] [33] [34]

[35]

[36]

[37]

[38]

[39]

[40]

[41]

[42]

[43]

[44] [45]

[46]

[47]

[48] [49]

[50] [51]

317

identification of phosphodiesterase-5 inhibitors, J. Chem. Inf. Model. 51 (2011) 1353–1363. A.R. Anighoro, G, BEAR, a Molecular docking refinement and rescoring method, Comput. Mol. Biosci. 3 (2013) 27–31. R.O. Juvonen, M. Kuusisto, C. Fohrgrup, M.H. Pitkanen, T.J. Nevalainen, S. Auriola, et al., Inhibitory effects and oxidation of 6-methylcoumarin, 7-methylcoumarin and 7-formylcoumarin via human CYP2A6 and its mouse and pig orthologous enzymes, Xenobiotica; Fate Foreign Compd. Biol. Syst. (2015), http://dx.doi.org/10.3109/00498254.2015.1048327, in press. A.M. Ferrari, B.Q. Wei, L. Costantino, B.K. Shoichet, Soft docking and multiple receptor conformations in virtual screening, J. Med. Chem. 47 (2004) 5076–5084. K.P. Madauss, S.J. Deng, R.J. Austin, M.H. Lambert, I. McLay, J. Pritchard, et al., Progesterone receptor ligand binding pocket flexibility: crystal structures of the norethindrone and mometasone furoate complexes, J. Med. Chem. 47 (2004) 3381–3387. H. Wang, M. Ye, H. Robinson, S.H. Francis, H. Ke, Conformational variations of both phosphodiesterase-5 and inhibitors provide the structural basis for the physiological effects of vardenafil and sildenafil, Mol. Pharmacol. 73 (2008) 104–110. U. Pentikainen, O.T. Pentikainen, A.J. Mulholland, Cooperative symmetric to asymmetric conformational transition of the apo-form of scavenger decapping enzyme revealed by simulations, Proteins 70 (2008) 498–508. U. Pentikainen, L. Settimo, M.S. Johnson, O.T. Pentikainen, Subtype selectivity and flexibility of ionotropic glutamate receptors upon antagonist ligand binding, Org. Biomol. Chem. 4 (2006) 1058–1070. P.A. Postila, G.T. Swanson, O.T. Pentikainen, Exploring kainate receptor pharmacology using molecular dynamics simulations, Neuropharmacology 58 (2010) 515–527. M. Ylilauri, O.T. Pentikainen, Structural mechanism of N-methyl-d-aspartate receptor type 1 partial agonism, PLoS One 7 (2012) e47604. N. Huang, B.K. Shoichet, J.J. Irwin, Benchmarking sets for molecular docking, J. Med. Chem. 49 (2006) 6789–6801. A. Gaulton, L.J. Bellis, A.P. Bento, J. Chambers, M. Davies, A. Hersey, et al., ChEMBL: a large-scale bioactivity database for drug discovery, Nucleic Acids Res. 40 (2012) D1100–D1107. A. Urzhumtsev, F. Tete-Favier, A. Mitschler, J. Barbanton, P. Barth, L. Urzhumtseva, et al., A ‘specificity’ pocket inferred from the crystal structures of the complexes of aldose reductase with the pharmaceutically important inhibitors tolrestat and sorbinil, Structure 5 (1997) 601–612. D. Tondi, F. Morandi, R. Bonnet, M.P. Costi, B.K. Shoichet, Structure-based optimization of a non-beta-lactam lead results in inhibitors that do not up-regulate beta-lactamase expression in cell culture, J. Am. Chem. Soc. 127 (2005) 4632–4639. L. Wright, X. Barril, B. Dymock, L. Sheridan, A. Surgenor, M. Beswick, et al., Structure-activity relationships in purine-based inhibitor binding to HSP90 isoforms, Chem. Biol. 11 (2004) 775–785. G.L. Card, B.P. England, Y. Suzuki, D. Fong, B. Powell, B. Lee, et al., Structural basis for the activity of drugs that inhibit phosphodiesterases, Structure 12 (2004) 2233–2247. J.V. Lehtonen, D.J. Still, V.V. Rantanen, J. Ekholm, D. Bjorklund, Z. Iftikhar, et al., BODIL: a molecular modeling environment for structure-function analysis and drug design, J. Comput. Aided Mol. Des. 18 (2004) 401–419. J.M. Wang, W. Wang, P.A. Kollman, Antechamber: an accessory software package for molecular mechanical calculations, Abstr. Pap. Am. Chem Soc. 222 (2001) U403. K.S. Watts, P. Dalal, R.B. Murphy, W. Sherman, R.A. Friesner, J.C. Shelley, ConfGen a conformational search method for efficient generation of bioactive conformers, J. Chem. Inf. Model. 50 (2010) 534–546. T.A. Halgren, Merck molecular force field.1. Basis, form, scope, parameterization, and performance of MMFF94, J. Comput. Chem. 17 (1996) 490–519. M.J. Vainio, J.S. Puranen, M.S. Johnson, ShaEP molecular overlay based on shape and electrostatic potential, J. Chem. Inf. Model. 49 (2009) 492–502. D.A. Case, T.A. Darden, T.E. Cheatham III, C.L. Simmerling, J. Wang, R.E. Duke, et al., Amber10, University of California, San Francisco, 2008. A. Jakalian, B.L. Bush, D.B. Jack, C.I. Bayly, Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method, J. Comput. Chem. 21 (2000) 132–146. Y. Duan, C. Wu, S. Chowdhury, M.C. Lee, G. Xiong, W. Zhang, et al., A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations, J. Comput. Chem. 24 (2003) 1999–2012. J.M. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kollman, D.A. Case, Development and testing of a general amber force field, J. Comput. Chem. 25 (2004) 1157–1174. J. Aqvist, Ion water interaction potentials derived from free-energy perturbation simulations, J. Phys. Chem.-Us. 94 (1990) 8021–8024. H.J.C. Berendsen, J.P.M. Postma, W.F. Vangunsteren, A. Dinola, J.R. Haak, Molecular-dynamics with coupling to an external bath, J. Chem. Phys. 81 (1984) 3684–3690. T. Darden, D. York, L. Pedersen, Particle mesh Ewald—an N. Log(N) method for Ewald sums in large systems, J. Chem. Phys. 98 (1993) 10089–10092. H.G. Petersen, Accuracy and efficiency of the particle mesh Ewald method, J. Chem. Phys. 103 (1995) 3668–3679.

318

S.I. Virtanen et al. / Journal of Molecular Graphics and Modelling 62 (2015) 303–318

[52] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, Numerical-integration of cartesian equations of motion of a system with constraints—molecular-dynamics of N-alkanes, J. Comput. Phys. 23 (1977) 327–341. [53] M.L. Connolly, Analytical molecular-surface calculation, J. Appl. Crystallogr. 16 (1983) 548–558. [54] B.R. Brooks, D. Janezic, M. Karplus, Harmonic-analysis of large systems: 1. Methodology, J. Comput. Chem. 16 (1995) 1522–1542. [55] W.L. DeLano, Use of PYMOL as a communications tool for molecular science, Abstr. Pap. Am. Chem. Soc. 228 (2004) U313–U314. [56] N.A. Baker, D. Sept, S. Joseph, M.J. Holst, J.A. McCammon, Electrostatics of nanosystems: application to microtubules and the ribosome, Proc. Natl. Acad. Sci. U. S. A. 98 (2001) 10037–10041. [57] S.P. Niinivehmas, K. Salokas, S. Lätti, H. Raunio, O.T. Pentikäinen, Ultrafast protein structure-based virtual screening with Panther, J. Comput. Aid. Mol. Des. 29 (2015) 989–1006. [58] J.D. Buynak, L. Vogeti, V.R. Doppalapudi, G.M. Solomon, H. Chen, Cephalosporin-derived inhibitors of beta-lactamase. Part 4: the C3 substituent, Bioorg. Med. Chem. Lett. 12 (2002) 1663–1666. [59] B.W. Dymock, X. Barril, P.A. Brough, J.E. Cansfield, A. Massey, E. McDonald, et al., Novel, potent small-molecule inhibitors of the molecular chaperone Hsp90 discovered through structure-based design, J. Med. Chem. 48 (2005) 4212–4215.

[60] A.H. Abadi, B.D. Gary, H.N. Tinsley, G.A. Piazza, M. Abdel-Halim, Synthesis, molecular modeling and biological evaluation of novel tadalafil analogues as phosphodiesterase 5 and colon tumor cell growth inhibitors, new stereochemical perspective, Eur. J. Med. Chem. 45 (2010) 1278–1286. [61] I. Akritopoulou-Zanze, J.R. Patel, K. Hartandi, J. Brenneman, M. Winn, J.K. Pratt, et al., Synthesis and biological evaluation of novel, selective, nonsteroidal glucocorticoid receptor antagonists, Bioorg. Med. Chem. Lett. 14 (2004) 2079–2082. [62] J.M. Wang, P. Morin, W. Wang, P.A. Kollman, Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA, J. Am. Chem. Soc. 123 (2001) 5221–5230. [63] C.R.W. Guimaraes, M. Cardozo, MM-GB/SA rescoring of docking poses in structure-based lead optimization, J. Chem. Inf. Model. 48 (2008) 958–970. [64] A. Weis, K. Katebzadeh, P. Soderhjelm, I. Nilsson, U. Ryde, Ligand affinities predicted with the MM/PBSA method: dependence on the simulation method and the force field, J. Med. Chem. 49 (2006) 6596–6606. [65] T. Sulea, Q.Z. Cui, E.O. Purisima, Solvated Interaction Energy (SIE) for Scoring Protein-Ligand Binding Affinities: 2. Benchmark in the CSAR-2010 scoring exercise, J. Chem. Inf. Model. 51 (2011) 2066–2081.