Help with unexpectedly large free-energy contribution for a Cys-to-Arg mutation using pmx

pmx
1

Hi everyone!

I am a new user of pmx, and I am trying to calculate the free-energy difference associated with a Cys-to-Arg mutation. I have performed the calculations, but the results look unexpectedly large to me. Could someone help me interpret the results and identify potential issues in my protocol?

The Cys-to-Arg mutation changes the net charge of the protein. To compensate for this charge difference, I introduced a dummy Cl⁻ ion that is alchemically coupled to Arg and scales together with it.

Another complication is that the mutation site is close to a GDP molecule, which might also affect the results.

However, after performing the calculations, I found that the work associated with the mutation was much larger than I expected. I am wondering whether this could be due to an issue with my protocol or whether I am missing something in the interpretation of the results.

Here are the details of my setup:

  • GROMACS: 2023.5
  • Force field: CHARMM36m
  • Ligand parameterization: CGenFF 5.0
  • pmx: Installed from the develop branch (git checkout develop)
  • Equilibration: 20 ns, using the last 10 ns to extract initial structures
  • Initial structures: Extracted at 100 ps intervals
  • Nonequilibrium MD simulations: 100 ps per trajectory

Since dH/dλ is calculated for the entire system, the reported work may include contributions from both the mutation and the alchemical transformation of the compensating chloride ion, including changes in its solvation environment. However, the magnitude of the work still seems unusually large to me.

Could you please share your thoughts on what might be causing this behavior? In particular, I would appreciate any suggestions regarding the treatment of the charge-changing mutation, the use of the compensating ion, and the simulation protocol.

Thank you in advance for your help!

Here is the topology of the dummy Cl:
[ moleculetype ]
; Name nrexcl
DUM_CL 3

[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 383 Cl rtp Cl q -1.0
1 CLA 383 DCL DCL 1 -1 35.45 DUM_CL 0 35.45
And vdw part:
[ atomtypes ]
P3 15 30.974000 0.000 A 0.383086448800 2.4476400 ; pyrophosphate phosphorus
DUM_CL 0 35.450000 0.000 A 0.0 0.0

Below I am attaching plots with dH/dl over the time and final dG calculation:

image
image964×755 19.5 KB

wplot
wplot2400×1800 219 KB

And here is the .mdp of the non-equilibrium transformation
integrator = md
dt = 0.002
nsteps = 50000 ; 100 ps

nstxout-compressed = 500000
nstlog = 500000
nstenergy = 500000
nstcalcenergy = 1

;nstxtcout = 50000
;nstvout = 50000
;nstfout = 50000
;nstenergy = 50000
;nstlog = 50000
;

cutoff-scheme = Verlet
nstlist = 20
vdwtype = Cut-off
vdw-modifier = force-switch
DispCorr = no
rvdw = 1.2
rlist = 1.2
rvdw-switch = 1.0
rcoulomb = 1.2
coulombtype = PME
;
tcoupl = V-rescale
tc_grps = Protein Non-Protein
tau_t = 0.1 0.1
ref_t = 300 300
;
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 5.0
compressibility = 4.5e-5
ref_p = 1.0
;
constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes
;
nstcomm = 100
comm_mode = linear
comm_grps = Protein Non-Protein
;
free-energy = yes
init-lambda = 0 ; start from state A
delta-lambda = 2e-05 ; complete transition in this number of steps
sc-coul = yes ; use soft-core also for coulombic interactions
sc-alpha = 0.3 ; soft-core
sc-sigma = 0.25 ;
sc-power = 1 ;
nstdhdl = 1 ; write to dhdl at each step

2

Hey,

firstly, to know the actual ddG of this perturbation you need the second leg of the thermodynamic cycle.

To improve convergence: (i) do not annihilate LJ parameters for the co-alchemical CL ion; (ii) run longer transition simulations; (iii) try setting in the mdp sc-function=gapsys, the rest of the sc parameters can be kept at their default values

3

Hi vgapsys!

Thanks for your reply and suggestions!

My second leg involves a similar transformation (Arg → Cys), but in a different protein conformation. I want to investigate how the mutation affects the relative stability of the two protein states. I am trying to follow a similar approach to the one described in this publication.

Yes, I have already tried suggestion (III), and it helped partially by reducing the gap between the forward and backward work distributions. However, the gap is still quite large (see the attached figures). I hope that suggestions (I) and (II) will help further!

image
image1234×531 106 KB

I have two additional questions:

  1. Should I be concerned about the non-Gaussian shape of the forward work distribution? I understand that BAR does not require the work distributions to be Gaussian, but could a deviation from Gaussian behavior be a symptom of insufficient sampling or some other issue with my simulation setup?

  2. Am I correct in understanding that work values (W) in the range of 500–2000 kJ/mol can be expected for charge-changing mutations, or would values of this magnitude suggest a problem with my setup or protocol?

I would appreciate any further advice. Thanks again for your help!