In this tutorial we will use pmx and Gromacs to set up alchemical free energy calculations for a relative free energy calculation study. We will further perform the simulations and estimate the actual ddG value.
The protein and ligand topologies for this case have already been generated and described in the study:
Gapsys, Perez-Benito et al. Large scale relative protein ligand binding affinities using non-equilibrium alchemy. 2020. Chemical Science
https://pubs.rsc.org/en/content/articlelanding/2020/SC/C9SC03754C
All the input/output files from this work can be found on the pmx github page:
https://github.com/deGrootLab/pmx/tree/master/protLig_benchmark
In the tutorial, we will use a pair of ligands (18625-1: Cl in ortho position and 18626-1: Cl in meta position) binding to the JNK1 kinase and estimate the relative free energy difference $\Delta\Delta G$ between these molecules. The experimentally measured double free energy difference for this edge is -3.22 kJ/mol.
For the calculations, we will construct the following thermodynamic cycle:
conda create -n PMX python=3.7
conda activate PMX
conda install numpy matplotlib scipy pip jupyterlab pandas rdkit
git clone https://github.com/deGrootLab/pmx
cd pmx
git checkout develop
pip install .We will set up the calculations in a non-equilibrium framework.
Firstly, two 6 ns long equilibrium simulations will be performed in the physical states A and B. This will be done for the ligand solvated in water and protein-ligand complex.
Afterwards, we will extract snapshots from the equilibrium simulations and run 80 rapid (50 ps each) transitions in the directions A->B and A<-B.
During the transitions, we will record the work required for the ligand morphing. Finally, the work will be related to $\Delta G$ by applying Crooks fluctuation theorem.
The overall workflow consists of the following steps:
The workflow in sketched in the following figure. The protocol is followed for estimation of both $\Delta G_{L}$ and $\Delta G_{PL}$.