I am planning on doing some ligand relative binding free energy calculations using pmx that involve charged mutations. The recommended protocol seems ‘one box two systems’. If I apply this to my problem, I would build two boxes in total:
If my understanding is correct, for one fast transition in box1 you basically sum the work for Protein-Ligand 1 → Protein- Ligand 2 and Ligand 2 → Ligand 1. This allows you to calculate the free energy difference by only two simulations, essentially calculating ΔGPL1,PL2 + ΔGL2,L1 instead of ΔGPL1,PL2 - ΔGL1,L2.
To keep the ligand in water system seperate from the protein-ligand system, positional restraints are applied on one central atom. This way translational entropy is impacted, but not the rotational entropy. In literature it is then explained:
The contribution of such position restraints to the translational partition function will be the same in both legs of the thermodynamic cycle and will cancel from the ΔΔG estimate.
Based on this, I have the following questions:
1 ) Am I correct in understanding that the contibution of the positional restraints cancel out because the free energy difference of restraining Ligand 1 is almost the same as the free energy difference of restraining Ligand 2? I have provided a pdf as attachment with the thermodynamic cycle and some equations to clarify this question. Is the analysis in the pdf correct?
2) If my understanding is correct, then I assume that “one box two systems” is only valid for small transformations. Is this correct?
3) If my understanding is not correct, how exactly does the effect of restraining translation cancel out?
Regarding point 2, suppose I restraint a central atom of the ligand for the “ligand in water” calculation and the same central atom of the ligand for the “ligand in protein” part of the box. Intuitively, I would now think that the translational movements of the ligand in the protein are slower, since only the protein can move while the ligand translation is hindered by the restraint. For small perturbations, I can understand that the ΔG for restraining the ligand position in the protein in state 1 (let us call this ΔG_PL1,PL1R) is very similar to -ΔG of removing the positional restraint in state 2 (let us call this ΔG_PL2R,PL2). So in mathemathical terms: ΔG_PL1,PL1R + ΔG_PL2R,PL2 = 0.
But if the first and second state are wildly different, then can the movements of the ligands in the protein site not also wildly differ in reality? As the translational relative velocities of the ligand in the pocket are impacted by the restraints in the calculations I would then think that the assumption of ΔG_PL1,PL1R + ΔG_PL2R,PL2 = 0 will not hold anymore… What is your view on this situation?
I think there are several concepts convoluted here:
regarding “slower translational movements”, while the restraint might have an effect on kinetics, in alchemical calculations, we are interested in thermodynamics only. The restraint may influence the dynamics, but the difference in kinetic partition function cancels out exactly;
the translational partition function function depends on volume, so if i) the box changes a lot between states A and B, and ii) the volume change is different between the legs of the cycle, there will be some contribution from the restraints. This contribution, however, is much smaller than the uncertainty in the alchemical calculations