amber pmemd free energy perturbation with charge correction
Amber 14,16
Amber14 force field for protein
GAFF for ligands
Particle Mesh Ewald MD
RESP charges, at HF/6-31G* level.
Torsion and angle parameters according to the Seminario approach
Windows per transformation: 24, lambda values 25
Assumed pH 7.4
Maestro Protein Preparation Wizard for protonation states of protein and ligand
open histidine conformation
Sampling time per lambda: 2 ns
Langevin dynamics, collision frequency of 2 ps-1
Cumulative sampling time per perturbation: 50 ns
Neutralizing plasma (No counter-ions)
SHAKE
TIP3P water
tishake=1
Avogadro
MBAR (pymbar)
Overlap convergence metrics
Cubic and Octahedral water boxes
Charge correction according to Rocklin et al. [1]
APBS
All molecular dynamics (MD) simulations and FEP simulations were performed with the Amber 14 (prerelease) and 16 softwares. Farnesoid X receptor and the ligands were solvated in a truncated octahedral box
of water molecules, extending at least 9 Å from the solute using the leap program in the Amber suite.
No counter-ions were used in the calculations and we have previously shown that they do not have any
influence on the calculated free-energy differences. Ligands were built from crystal structures from D3R GC2 Stage 1 as starting structures. The FEP simulations were run with the pmemd module of Amber, using the dual topology scheme with both ligands
in the topology files. We employed 25 states with λ = 0.000, 0.0025, 0.0050, 0.0075, …, 0.0850, 0.0900, 0.0900, and 1.0000, using a linear transformation of the potentials.
Electrostatic and van der Waals interactions were perturbed concomitantly, using soft-core potentials for both types of interactions.
The soft-core potentials were not only used for atoms differing between the ligands molecules, but for all atoms in the
ring systems neighboring the perturbed group.
To make the calculations comparable between the two versions of Amber, we used the keyword
tishake=1 for the Amber calculations.
For each λ value, we first performed 100 steps of minimisation, with the heavy atoms of the protein and ligands
restrained towards the starting structure with a force constant of 418.4 kJ/mol/Å2.
This was followed by 20 ps constant-volume equilibration with the same restraints and 2 ns constant-pressure equilibration
without any restraints.
Finally, a 2 ns production simulation was run, during which structures were sampled every 2 ps.
In the MD simulations, bonds involving hydrogen atoms were constrained with the SHAKE algorithm, allowing for a time-step of 2 fs.
In all simulations, the temperature was kept constant at 300 K using Langevin dynamics with a collision frequency
of 2 ps-1, and the pressure was kept constant at 1 atm using a weak-coupling isotropic algorithm with a relaxation time of 1 ps.
Long-range electrostatics were handled by particle-mesh Ewald (PME) summation with a fourth-order B spline
interpolation and a tolerance of 10–5. The cut-off for Lennard-Jones interactions was set to 8 Å.
Relative binding free energies between two ligands, L0 and L1 (∆∆Gbind),
were estimated using a thermodynamic cycle that relates ∆∆Gbind to the free energy of alchemically
transforming L0 into L1 when they are either bound to the host, ∆Gbound, or are free in solution, ∆Gfree
∆∆Gbind = ∆Gbind(L1) – ∆Gbind(L0) = ∆Gbound – ∆Gfree
where ∆Gbound and ∆Gfree were estimated by the Bennett acceptance-ratio method (BAR). Free
energies were also calculated by multi-state BAR, thermodynamic integration, and exponential
averaging, using the pymbar software.
Charge corrections for magnitude of free-energy and standard error. http://signe.teokem.lu.se/~ulf/Methods/ChargedFEPCorrections.html .
To estimate the convergence of the various perturbations, several different overlap measures were employed.
We calculated the Bhattacharyya coefficient for the energy distribution overlap (W), the Wu & Kofke overlap measures of the
energy probability distributions (KAB) and their bias metrics (PI),
the weight of the maximum term in the exponential average (wmax), the difference of the forward and backward exponential
average estimate (DDGEA), and the difference between the BAR and TI estimates (DDGTI ),
although this difference may also reflect the integration error in TI.
[1] Rocklin et al., Calculating the binding free energies of charged species based on explicit-solvent simulations employing
lattice-sum methods: An accurate correction scheme for electrostatic finite-size effects, 139, 2013.
Dummy
Maestro/Omega/MGLTools/AutoDock Vina
Assumed pH 7.4
Tautomers considered
Gasteiger charges
Maestro's prepwizard was used to add hydrogens using default parameters.
All the waters were retained while GAFF force field parameters were assigned to the protein and all
ligands. The binding site residues, which were considered to be residues 56, 65, 200 and 210, were energy minimized.
Ligand conformational libraries were generated using OMEGA. The most stable
conformation found for each ligand was then used as a starting point for docking.
Preparation for docking with Vina was done using MGLTools where the ligand PDB files from
OMEGA were converted to PDBQTs by assigning charges and atom types. The binding site residues
above were treated as flexible and the grid box determined.
Iterated Local Search global optimizer search method
Exhaustiveness=50 #exhaustiveness of global search (default=8)
Vina scoring funtion (empirical + knowledge-based function)
Num_modes=5 #max number of poses to generate
Energy_range=5 #energy difference (kcal/mol) between the best and worst binding mode
Docking runs were executed with the above specified parameters while default values
were applied for the rest of the variables. The top 5 poses from the Vina docking are submitted with this protocol.