Absolute FEP result too high

User discussions
1

GROMACS version: 2024
GROMACS modification: No

I calculated protein-ligand absolute binding free energy for a library of 15 molecules. I got results as low as -1.5 kcal/mol and -40.24 kcal/mol. I am just wondering what might be the reason for such a high result. My protein is transmembrane protein and I prepared the system using CHARMM-GUI with DPPC lipid bilayer and 0.15 M NaCl aqueous layer. I equilibrate and simulate with the script provided there and get boresch restriction values from the last 2ns. Then I proceed with the FEP lambda simulations using GROMACS with mdp options noted below. Overall my results are in -8 to -10 kcal/mol range. In three cases I get abnormally large values as noted above while experimentally I thought binding free energy should plateau around -15 kcal/mol (picomolar Kd). I can observe in these cases that the ligands have carboxylic acid interacting with two arginines in the protein. All the ligands have a hydrophobic tail in common that docks into hydrophobic submerged in the lipid bilayer. Everything seems to be working fine so I don’t know if I should accept these results. Any possible explanation is welcome so that I can follow up on those.

; BONDS

;----------------------------------------------------

constraint_algorithm = lincs ; holonomic constraints

constraints = h-bonds ; hydrogens only are constrained

lincs_iter = 1 ; accuracy of LINCS (1 is default)

lincs_order = 6 ; also related to accuracy (4 is default)

lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)

continuation = no ; formerly known as ‘unconstrained-start’ - useful for exact continuations and reruns

; NEIGHBOR SEARCHING

;----------------------------------------------------

cutoff-scheme = Verlet

ns-type = grid ; search neighboring grid cells

nstlist = 20 ; 20 fs (default is 10) gromacs itself recommended increase in order to gain 10-15% performance

rlist = 1.9 ; short-range neighborlist cutoff (in nm) increased from 1.2

pbc = xyz ; 3D PBC

; ELECTROSTATICS & EWALD

;----------------------------------------------------

coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics

rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)

ewald_geometry = 3d ; Ewald sum is performed in all three dimensions

pme-order = 6 ; interpolation order for PME (default is 4)

fourierspacing = 0.10 ; grid spacing for FFT

ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb

; VAN DER WAALS

;----------------------------------------------------

vdw-type = Cut-off

vdw-modifier = Potential-shift-Verlet

verlet-buffer-tolerance = -1 ; To override the automated buffer setting , see below for explanation

rvdw = 1.2 ; short-range van der Waals cutoff (in nm)

2

I explored protein conformation change as a result of those two arginines interacting with carboxylic acid functional group and wondered if conformation change to some higher energy state can explain the jump in FEP. Visually the protein and the residue rotamers in the active site didn’t change much. I measured the average of three distances that makes up a triangle around the active site and observed that the values don’t change much for free protein or the 15 protein- ligand complexes. So I assume that this is not a viable explanation. I don’t know if I can make a better analysis in this regards using Gromacs functions

3

Also the above mdp is for free energy option couple-intramol=no. I increased the rlist to circumvent the “perturbed, excluded non-bonded pair interactions beyond the pair-list cut-off” error. I simulated the ligand solvation with these option. I also tried simulating all protein-ligand and ligand in water simulations with couple-intramol=yes without increased rlist and got -21 kcal/mol. It is an improvement over -40 but still an overestimation in my opinion. I am also not 100 % sure if the thermodynamic cycle is complete like this and if analytical restriction energy calculation stays the same.