[gmx-users] forward and backward states for the binding free energy energy

[Qasim Pars]

Dear users,

I was trying to calculate the binding free energy of the receptor-ligand
complex with the Crooks Gaussian Intersection (CGI) method
using GROMACS2016. The steps followed are here:

#FORWARD STATE:
I added the below free energy commands to the em.mdp, nvt.mdp, npt.mdp and
md.mdp for the state A:

; Free energy control  (state A)
free-energy              = yes
init-lambda              = 0
delta-lambda             = 0
sc-alpha                 = 0.5
sc-power                 = 1
sc-sigma                 = 0.3
sc-coul                  = yes
couple-moltype           = ligand
couple-intramol          = no
couple-lambda0           = vdw-q ;(state A)
couple-lambda1           = none  ;(state B)
nstdhdl                  = 100

After getting 100 ns (with 1000 frames) production run (md.mdp), I checked
whether the system has converged or not. The RMSD shows that the system has
converged after 50 ns. Then 50 starting snapshots were extracted from the
last 50 ns of the trajectory (gmx trjconv -s md.tpr -f md.xtc -o frame.gro
-pbc mol -ur compact -sep -b 50001 -skip 10). Then 50 forward simulations
of 100 ps were carried out for receptor-ligand complex with a 1-step
switching process. The free energy control part of the input file used for
100 ps forward state simulations is here:

; Free energy control  (state A-->state B: 0--->1)
free-energy              = yes
init-lambda              = 0
delta-lambda             = 2e-5 ;(1/nsteps)
sc-alpha                 = 0.5
sc-power                 = 1
sc-sigma                 = 0.3
sc-coul                  = yes
couple-moltype           = ligand
couple-intramol          = no
couple-lambda0           = vdw-q ;(state A)
couple-lambda1           = none  ;(state B)
nstdhdl                  = 100

#BACKWARD STATE:
I added the below lines to the em.mdp, nvt.mdp, npt.mdp and md.mdp for the
state B (only the init-lamba has been changed to 1):

; Free energy control  (state B)
free-energy              = yes
init-lambda              = 1
delta-lambda             = 0
sc-alpha                 = 0.5
sc-power                 = 1
sc-sigma                 = 0.3
sc-coul                  = yes
couple-moltype           = ligand
couple-intramol          = no
couple-lambda0           = vdw-q ;(state A)
couple-lambda1           = none  ;(state B)
nstdhdl                  = 100

Getting 100 ns production run, checking RMSD.... and so on. The steps are
the same as the forward state. The free energy part of input file used for
100 ps backward state simulations is as follows (the only difference is
delta-lambda):


; Free energy control  (state B-->state A: 1--->0)
free-energy              = yes
init-lambda              = 0
delta-lambda             = -2e-5 ;(-1/nsteps)
sc-alpha                 = 0.5
sc-power                 = 1
sc-sigma                 = 0.3
sc-coul                  = yes
couple-moltype           = ligand
couple-intramol          = no
couple-lambda0           = vdw-q ;(state A)
couple-lambda1           = none  ;(state B)
nstdhdl                  = 100

#ANALYSIS
The result of Crooks Gaussian Intersection (CGI) analysis (~20 kcal/mol) is
not compatible with literature (5,8 kcal/mol).


#QUESTIONS
1-) The mean of the forward state is ~15 kcal/mol. That is too big, right?
Maybe the ligand doesn't get decoupled in the forward state?

2-) Is the free energy control stuff used in the em.mdp, nvt.mdp, npt.mdp
and md.mdp correct for the forward and backward/reversible state
simulations? I must use the same free energy control parameters in both
states except init-lambda?

3-) Is the free energy control stuff used for 100 ps simulations of the
forward and backward/reversible state simulations correct? The only
difference between both states must be delta-lambda?

4-) Is the order of [ intermolecular_interactions ] (please see the below)
in the topology file belonging to complex structure right?

...
[ molecules ]
; Compound        #mols
ligand                 1
Protein             1
SOL             4184
NA                 15
CL                 12

[ intermolecular_interactions ]
[ bonds ]
....

[ angle_restraints ]
...
...

[ dihedral_restraints ]
...
...
...

5-) I added the state B parameters to the [ intermolecular_interactions ]
parts. That should be enough to define the state B of a molecule in
GROMACS2016?

Thanks in advance.

Qasim Pars