I’m working on free energy calculations for a charged ligand (+1) bound to a protein. As I’m new to handling charge-changing perturbations, I would appreciate some guidance.
I attempted the co-alchemical ion approach by adding a chloride ion to the ligand’s topology file (TKT) and using the .mdp setting couple-moltype = TKT to ensure the ligand and the ion couple and decouple simultaneously. I performed 500 ps transitions across 95 snapshots extracted from a 10 ns trajectory.
However, the results show highly negative energy values for the state 1 to state 0 transition, leading to a very negative binding free energy (~ -16 kcal/mol). This result deviates significantly from the experimental value (~ -9 kcal/mol) and has a large standard deviation (~5 kcal/mol).
Did I implement the co-alchemical ion approach correctly? Alternatively, are there other recommended approaches for handling these calculations? I’ve read about post-calculation correction schemes (which I assume would not require introducing an ion during the calculation) and the dual-system single-box approach, but I haven’t tried either of these before.
Would it be advisable to explore other methods, or should I adjust the current approach? Any advice or suggestions would be greatly appreciated!
it appears that, in principle, the co-alchemical ion approach has worked, because if you would have had e.g. a (dis)appearing charge the obtained delta-deltaG would have likely been a lot larger in absolute value (an ion has a hydration free energy on the order of 100 kcal/mol).
there are 3 things that come to mind:
the (dis)appearing ion should not interact with the protein/ligand. Did you apply a position or distance restraint to make sure they stay out of each other’s way?
instead of a co-alchemical ion you could also try the double-system/single box approach, in which you do the protein leg and the (inverse) solvent leg in the same box. This saves the assumptions around the co-alchemical ion. Note that also here a weak restraint should be applied to make sure that the two subsystems remain at least 1 nm or so apart
the statistical error you get is substantial. You could try multiple repeats and/or longer transitions (say 1 or 2 ns) to see if reducing the statistical error would improve the accuracy.
Thank you for your prompt response and valuable insights.
No, I have simply placed the ligand and ion far apart from each other. The only restraints applied are Boresch restraints between the ligand and the protein. If I were to apply positional restraints, would it be necessary to introduce a correction term for these types of restraints as well?
Thank you for the suggestion; I will definitely try that.
When you refer to multiple repeats, are you suggesting repeating the equilibrium MD simulations used to extract the snapshots, or just repeating the transitions? Additionally, do you think it would be beneficial to increase the duration of the equilibrium MD, for instance, from 10 ns to 100 ns?
in principle yes, but it is usually a safe assumption that this restraint is not doing any work as a function of lambda, so this can typically be ignored.
yes, it would be most unbiased to repeat also the equilibrium runs (as they show chaotic behavior they will sample independent equilibrium trajectories). The major sampling issue is probably in the non-equilibrium transitions so increasing their length will likely help more than increasing the length of the equilibrium runs.