The HREMD-based simulations utilized a modified version of the open-source Python alchemical free energy code YANK, which is built on the OpenMM GPU-accelerated molecular simulation library [20 , 39 (link)]. We performed our simulations using a generalized Born (GB) implicit representation of water [25 ]. A Langevin dynamics integrator with a 2 fs time step and a 0.1 ps−1 collision frequency was used, with a bath temperature of 298 K, and bonds to hydrogen were constrained by the CCMA method [55 (link)]. A flat-bottom restraint was implemented to keep the ligand in the vicinity of the protein while allowing it to sample in an unbiased way all spatially available and physically reasonable conformational space consistent with binding. The specific choices made for this potential are described below. Hamiltonian replica exchange [32 ] was used to improve sampling, along with a number of improvements described below. Simulations were run on GPU computing resources provided by XSEDE, including the NCSA Forge and Lincoln clusters.
All preliminary tests of simulation parameters and the 10-fold replicate test of simulation consistency were performed with 1-methylpyrrole, a known binder. The ability of our approach to differentiate binders from non-binders was then examined by introducing another three ligands: benzene, a small binder, p-xylene, a larger binder which requires conformational change of Val111, and phenol, a nonbinder, as a control [15 (link)]. By using p-xylene, the ability of the method to sample all relevant biomolecular motions of the protein can be examined. The system used in our simulations is shown in Fig. 1.
With sufficient sampling of all relevant binding states, the simulations here can also be used to calculate the protein-ligand free energy of binding. For this purpose, we additionally performed HREMD simulations of the ligand alone, in implicit solvent, with the same parameters as described above.