I am performing binding free energy calculations for a protein-ligand system using fep_lambdas for decoupling. Among them, I turned on the couple-intramol = yes, which led me to need to perform decoupled calculations on the ligand under the four systems:
After the calculation, I found that there was almost no difference in the energy difference from the BAR calculation for system A and B (for each system I did three repeated calculations). Meanwhile, I checked the BAR calculation results of system C and D, and the ligand dissolution free energy subtracted from the two systems was similar to the result obtained after I turned off the couple-intramol.
I think the small energy difference between A and B stems from the fact that the ligand-binding site I studied is a binding pocket that is deeply buried in the interior of the protein (PDB: 6B6G). This leads to the fact that in complex system, the environment around ligand is similar both in solution and in vacuum (I counted that in system A, the number of accessible water molecules around ligand is less than 10).
So I think that for this complex system, even though the couple-intramol is turned on, only calculating the energy difference between the A and C systems can still return valuable binding free energy information.
I would like to know whether my analysis and conclusion are reasonable, welcome to give comments, thank you!
Simulation B is not necessary to perform (as you almost conclude).
Think of it like this. When calculating the absolute binding free energy in the binding site you want the process:
Protein + Ligand → Ligand in vacuum.
and when you calculate the absolute solvation free energy you want the process:
Solvent + Ligand → Ligand in vacuum.
However, when you use intramolecular decoupling, the end state in the simulations is not Ligand in vacuum, but rather vacuum. So, for absolute binding/solvation free energies you need the extra step:
Ligand in vacuum → vacuum, i.e. simulation D in your list.
So, (dG from C) - (dG from D) gives you the absolute solvation free energy, or the transfer free energy from vacuum into solution. (dG from A) - (dG from D) gives you the absolute binding free energy, or the transfer free energy from vacuum into the binding site.
What you are usually interested in, when calculating the binding free energy, is the relative free energy difference when going from solution to the binding site. For that, you only need simulations A and C. The elimination in vacuum is no longer relevant (as you imply above).
However, if you use restraints in simulation A, which is almost certainly necessary to keep the ligand in place, things get more complicated … I still don’t think you need simulation B, but I’m not quite sure at the moment.
I treat the C-D simulation and the A-D simulation as calculating the dissolution free energy of ligand in the solution phase and the protein phase, respectively.
But I want to make sure that we can omit B only if this binding site is a buried protein site. If the binding site is locating on protein surface, then the simulation of B is still necessary, right?
In addition, you mentioned about restraints, which I saw discussed at the end of the GROMACS solvation free energy tutorial (Calculating free energy — GROMACS tutorials https://tutorials.gromacs.org documentation), and the tutorial suggested pulling code to keep the ligand in place. The restraint you mentioned should be position restraint or pulling codes? I used POSRES for this system. In fact, I am also confused about “the contribution of quadratic potential to free energy” mentioned in the tutorial. In my initial reference a absolute freedom can calculate the tutorial is completely without using any of the constraints of small molecule location: http://www.alchemistry.org/wiki/Absolute_Binding_Free_Energy_-_Gromacs_2016
I may be wrong, but I don’t see how B ever comes into play. Decoupling the ligand from the protein while the complex is in vacuum does not seem relevant to me. Or am I misunderstanding the setup of B?
Indeed, http://www.alchemistry.org/wiki/Absolute_Binding_Free_Energy_-_Gromacs_2016 is a good source of information for the whole calculation process and also how to handle restraints. You can use pulling to keep the ligand in place, but I think the orientational restraints recommended on that page are good, i.e., one distance, two angles and three torsions. Also, remember to decouple the restraints in the Protein - Ligand leg. The compensation for the “contribution of the quadratic potential to the free energy” corresponds to the analytical compensation described in the Restraints section on the Alchemistry page, with the reference to Boresch et al. But in the tutorial only a single pulling umbrella is used to restrain the ligand, instead of the set of restraints.
I asked Professor Hess under another question, when I also asked whether it was necessary to decoupling the ligand of the complex system under vacuum condition. TI - Free Energy - couple-intramol=no Prof.Hess replied yes.
I’m doing calculations on different types of systems. I think my understanding of the parameter couple-intramol is still not enough.
Why can’t we just calculate dGA - dGC to calculate the energy difference of process: Protein-Ligand in Solvent → Ligand in Solvent?
Does the process corresponding to the binding free energy that people usually think of must refer to the process: Protein-Ligand in Solvent → Ligand in Vacuum?