Applying distance restraints to maintain bifurcated H-bond observed in experimental data

GROMACS version: 2020
GROMACS modification: No

Hi GROMACS experts,

I am simulating a protein with covalently attached inhibitor described using CHARM36m for proteins. This protein features a characteristic oxyanion hole created in part by two main chain amides. The crystal structure shows the inhibitor possess a alkoxide occupying the oxyanion hole. In an effort to simulate the protein-inhibitor complex in this experimentally derived configuration, I’d like to ensure it is maintained throughout the simulation. What would be the most appropriate approach to achieve this goal? I’m still fairly new to this but my first consideration is to apply simple harmonic distance restraints (bond type 6) for the two amide hydrogen with the alkoxide.

Do you know that a biasing potential is necessary? Is the bifurcation not preserved in an unbiased simulation?

Hi Justin,

In an unbiased test simulation the bifurcation was not preserved. There is a protonated Histidine about ~5A from the alkoxide that seems to encourage the group to rotate toward it. This protonation state is informed from neutron diffraction.

For these simulations, I will be primarily be looking at secondary and tertiary structure information.

Then the distance restraints as proposed would be appropriate.

Thanks Justin,

So then to implement this, it looks like I’ll need to add a few lines near the end of my .itp topology.
Formatted something like this?

[ distance_restraints ]
; ai aj type index type’ low high bo kb
2310 2215 1 0 6 0.0 0.4 2.4 ?
2310 2233 1 0 6 0.0 0.4 2.4 ?

ai and aj are atoms will be the appropriate index
type is 1 because its always 1?
index is both 0 because these should be controlled together
bond type is 6 for simple harmonic restraints
low is the min distance
high is the max distance
bo is the distance is the distance (in nm) is taken from the minimized structure
kb is the force constant to be applied. Is there a resource where I might find an appropriate magnitude for the force constant?


I think it’s mostly an empirical effort to make the geometry work. In principle, one could probably do a QM vibrational analysis of the interaction to see at what frequency that mode of interaction vibrates. I’ve never seen anyone do that, though.

Thank you for your help Justin.

For anyone trying to do this in the future, please remember to enable the use of distance restraints in your .mdp file.