Protein geometry distortion during MD simulation

GROMACS version:2020.4
GROMACS modification: No
Greetings to all!
I am performing a 100 ns MD trajectory calculation for the complex of two proteins using CHARMM36 FF (2020). As a result, I get a model of a complex with a rather strongly distorted protein geometry (violation of planarity for planar groups, etc.), so coot is quite unhappy with it :). I use distance restrains for some important hydrogen bonds (“bond type 6” inside one protein chain and “pull code” for intermolecular interactions). How can the geometry of the model, obtained during the MD calculation, be optimized? Are there some tools for this in gromacs? Alternatively, do I need to set the simulation initially by imposing some additional restrictions to achieve good geometry of the amino acid residues and protein main chain? Thank you very much!

Hi,

I get a model of a complex with a rather strongly distorted protein geometry (violation of planarity for planar groups, etc.)

Which molecules/functonal groups have planar distortion? If those are standard residues, then it could be that you did not set your mdp parameter in line with CHARMM36 force field. Or in alternative you did not properly relax the starting structure in order to release possible steric clashes.

Best regards
Alessandra

Dear Alessandra! Thank you very much for the answer.
Well these are the standard residues like His, Trp, Tyr etc. What do you mean to set mdp parameter in line with CHARMM36 force field? My md.mdp looks like this:

title = Protein MD simulation
define = -DDISREA ; distance restrain
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000000 ; 2 * 50000000 = 100000 ps (100 ns)
dt = 0.002 ; 2 fs
; Output control
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save coordinates every 10.0 ps
; Bond parameters
continuation = yes ; continuing from NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs-iter = 1 ; accuracy of LINCS
lincs-order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns-type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme-order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate
tau-t = 0.1 0.1 ; time constant, in ps
ref-t = 310 310 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau-p = 2.0 ; time constant, in ps
ref-p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen-vel = no ; continuing from NPT equilibration

; Pull code for intermolecular distance restrains
pull = yes
pull-ngroups = 5
pull-ncoords = 3
pull-group1-name = a_496
pull-group2-name = a_1234
pull-group3-name = a_1235
pull-group4-name = b_3440
pull-group5-name = b_3443
pull-coord1-type = umbrella
pull-coord2-type = umbrella
pull-coord3-type = umbrella
pull-coord1-geometry = distance ; simple distance increase
pull-coord2-geometry = distance ; simple distance increase
pull-coord3-geometry = distance ; simple distance increase
pull-coord1-groups = 1 5
pull-coord2-groups = 2 4
pull-coord3-groups = 3 5
pull-coord1-dim = Y Y Y
pull-coord2-dim = Y Y Y
pull-coord3-dim = Y Y Y
pull-coord1-rate = 0.0 ; not pull, just distance restraint at a reference distance bw 2 groups
pull-coord2-rate = 0.0 ; not pull, just distance restraint at a reference distance bw 2 groups
pull-coord3-rate = 0.0 ; not pull, just distance restraint at a reference distance bw 2 groups
pull-coord1-k = 1000 ; kJ mol^-1 nm^-2
pull-coord2-k = 1000 ; kJ mol^-1 nm^-2
pull-coord3-k = 1000 ; kJ mol^-1 nm^-2
pull-coord1-init = 0.30
pull-coord2-init = 0.28
pull-coord3-init = 0.28
pull-nstxout = 0
pull-nstfout = 0

Is there something wrong with it that doesn’t match CHARMM36 FF?

Hi,
I guess you use the pull code between atoms, the applied potential is in “competition” with the force field potential, thus the use of not properly calibrate values for pulling could results in underised distortion of other interactions (e.i aromatic ring).

Best regards
Alessandra

Hi!
OK, but the geometry is far from ideal even for a model calculated without any additional constraints (pull code, distance restrains or position restrains). So, may be there is a way to correct the geometry of the resulting model in gromacs? Or the way to make an averaged model with good geometry for the entire trajectory or part of it?
Thanks!

Yours sincerely
Oleg

Hi,

There are several steps to consider:

  1. generation of the starting structure (e.i from experimental structure or from homology modelling, or).
  2. once you have a reasonable starting structure, you can energy minimize it (to assure that you are in one of the minima of the energy landscape and to eliminate steric clash)
  3. you can equilibrate the structure to get to the desidered temperature and pressure (using using position restrain for the macromolecule)
  4. then you can start to sample configuration

Did you see distorsion in the geometry after all those steps? Then I suggest to check carefully the starting structure.

Concerning the average structure, it is normal that an time-averaged structure show distortion in the amino-acid side chain geometry. For example the time-average of a rotating aromatic ring will be points on a line (not a ring).

Best regards
Alessandra

Thank you so much for your answers. I will analyze the structures in detail at all stages.
With best wishes
Oleg.