Strange PMF profile after umbrella sampling of the peptide pulling

[Solly]

GROMACS version: 2022.3
GROMACS modification: No

Dear GROMACS users!

I’ve gor a problem with calculating my PMF profile after umbrella sampling of protein-peptide pulling process. I’ve made an umbrella simulation for the complex of the beta-amyloid peptide with a receptor protein and calculated the binding free energy from PMF profile. It was very close to experimental values. Then i’ve thied to do umbrella sampling with the same receptor and modified peptide. I’ve made two isoforms which were shorter than initial. The pulling was ok and the umbrella sampling too as I can see looking at the trajectories. The distribution histograms were overlaping pretty good with no free spacing between them except a couple of small intervals that cannot reflect dramatically the whole profile and can be easily removed with additional md in new windows. But when I’ve tried to build the pmf profile for each of these new systems, I’ve got something weird in both systems. Here i’ve attached the data for one of the peptides. For the second the situation is very similar.
I’ve attached the figure of the PMF profile and histo.xvg.

Does anybody know how can such strange profile emerge?

histo.xvg (84.6 KB)

My mdp file was :
title = Umbrella pulling simulation
define = -DPOSRES_A
; Run parameters
integrator = md
dt = 0.002
tinit = 0
nsteps = 5000000 ; 10 ns
nstcomm = 10
; Output parameters
nstxout-compressed = 5000 ; every 10 ps
nstenergy = 5000
; Bond parameters
constraint_algorithm = lincs
constraints = h-bonds
continuation = yes
; Single-range cutoff scheme
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
rlist = 1.4
rcoulomb = 1.4
rvdw = 1.4
; PME electrostatics parameters
coulombtype = PME
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in two groups
Tcoupl = V-rescale
tc_grps = Protein Non-Protein
tau_t = 1.0 1.0
ref_t = 310 310
; Pressure coupling is on
Pcoupl = Berendsen
pcoupltype = isotropic
tau_p = 2.0
compressibility = 4.5e-5
ref_p = 1.0
refcoord_scaling = com
; Generate velocities is off
gen_vel = no
; Periodic boundary conditions are on in all directions
pbc = xyz
; Long-range dispersion correction
DispCorr = EnerPres
; Pull code
pull = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_group1_name = abeta
pull_group2_name = rage
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_dim = N N Y
pull_coord1_groups = 1 2
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.0 ; restrain in place
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull-group1-pbcatom = 4801
pull-pbc-ref-prev-step-com = yes

Thank you in advance for the answers!
Best regards,
Anna

[scinikhil]

This is probably poor sampling, you dont have overlaping smooth bell-shaped curve.

[Solly]

Hello again everybody!
I did an additional sampling and now my histogramm looks that way :
histo.xvg (181.1 KB)

But again my profile is bad :

So I suppose the problem is not in the sampling

[scinikhil]

What you are trying to do? what you are pulling?

[Solly]

I am attempting to calculate the binding coefficient for the RAGE protein with small peptides. I have successfully achieved this with the beta-amyloid peptide, and the resulting coefficient was close to the experimental value (IJMS | Free Full-Text | Docking and Molecular Dynamics-Based Identification of Interaction between Various Beta-Amyloid Isoforms and RAGE Receptor - my work, including all data on previous umbrella sampling, can be found here). Now, I need to perform umbrella sampling for shorter peptides, specifically halves of the full-length peptide. I am following the same procedure as in the previous successful case. I begin by taking the complex of RAGE and a peptide, namely Aβ17-42, then perform a pulling simulation to extract initial conformations for windows from this trajectory. I then perform umbrella MD for 10 ns, which was sufficient in the previous case. I have checked five random trajectories and they appear fine, with no artifacts. The distributions of the center of mass (COM) are good, as evidenced by my previous message. However, the resulting potential of mean force (PMF) is quite strange, especially considering that it was good for the whole peptide and the same receptor.

Thank you for your help,
Anna

[Solly]

Hello Gromacs users!
I still have problems with my PMF graph.
I’ve rerun all the md from the very start. I’ve removed pull-group1-pbcatom value from my mdp file.
[md_umbrella.mdp|attachment]
(upload://hXnORMuVSpDmQQAJ5DJ12JVGl3P.mdp) (1.9 KB)

The plots look better, but the energy range is still to low.
Could anybody provide any guidance on how to solve the problem?

Best regards,
Anna

[oSalad]

What is the experimental dG? You may have a very weak binder.

[Solly]

The dissociation constant (Kd) of Aβ42 with the soluble extracellular part of RAGE (sRAGE), determined using microscale thermophoresis (MST), was 1.0 ± 0.2 µM. That corresponds to — 8,171 kcal/mol. The dissociation constant (Kd) of Aβ17-42 was almost the same, while Kd of Aβ1-16 was too small to measure. So I’ve got half the dG value as in the experiment.