Reconnaissance metadynamics and MD refinement of docking poses
Gromacs 4.6.2, Amber 14, Acpype
Assumed pH 7.4
GAFF force field for ligands
FF14SB force field for proteins
TIP3P water model
Octahedral solvation box extending 8 Angstrom outside complex
Particle mesh Ewald with a cutoff of 10 Angstrom
The antechamber module of Amber 14 was used to parametrise the ligands with GAFF force field. The tLeap module integrated with Amber 14 was used to generate complex topologies where the protein was treated with ff14sb
force field. The complexes were immersed into a truncated orthorhombic box with TIP3P water molecules such that no protein atom was within 8 Angstrom of the box edge. Acpype was used to convert Amber generated topology and coordinates to Gromacs compatible files.
Standard MD simulation
Length 50 ns
Clustering using 2 Angstrom cutoff
Starting poses were generated by Autodock Vina (for details see "Multi-target Docking using Autodock Vina" submission), with docking attempted for all FXR structures in the PDB. For each ligand, the complex predicted to have the lowest binding free energy was selected.
Energy minimization was performed with steepest-descent for 200 steps. NPT equlibration (1 ns) was performed by using a position restraint on the protein C-alpha atoms and ligand heavy atoms. The temperature was kept at 310 K with a velocity rescaling thermostat and the pressure was kept at 1 bar with a berendsen barostat.
The restraints were then removed and the simulation was continued as long as the ligand stayed close to its starting conformation, to verify the stability of the pose. More specifically, the simulation was terminated when the ligand heavy-atom RMSD exceeded 0.25 nm from a reference configuration taken as the average over the first nanosecond of unrestrained simulation (to allow for relaxation with the new force field compared to docking), or when the simulation length reached 51 ns. All simulations were performed using Gromacs 4.6.2.
For three ligands (22, 27, and 32), which were stable in the MD simulations, the last MD configuration was used as a starting point for enhanced sampling simulations using the reconnaissance metadynamics (Recon) method, which explores the space spanned by a number of collective variables (CVs) much faster than regular MD simulations. The Recon simulations were performed using Gromacs 4.6.2 together with the PLUMED 1.3 plugin for free-energy calculations.
As CVs, one dihedral representing each rotable bond was used. The number of CVs were 8, 7, and 8, respectively, for the three ligands (22, 27, and 32), and they were found by visual inspection. Each simulation was run for ~32 ns and the following Recon parameters were used:
onion height: 1.0 kJ
onion width: 1.5
onion deposition stride: 1 ps
basin tolerance: 0.2
basin expand parameter: 0.3
basin initial size: 1.5
clustering frequency: 100 ps (1000 samples per clustering)
Each Recon trajectory was clustered together with the corresponding MD trajectory using the GROMACS g_cluster utility, using the "Gromos" clustering method and a cutoff of 2 A. Clusters with less than 10 structures were discarded. For each of the ~5 largest clusters, an MD simulation was started to investigate the stability, followed by cluster analysis of the MD trajectory. A manual sorting of the final poses was performed, based on the apparent stability in the Recon simulation and in the MD simulation.
For the remaining 33 ligands, the same poses were submitted as for "MD refinement of docking poses".