GROMACS Calculation of MM-PBSA Binding Free Energy Tutorial
There are many methods for calculating binding free energy, such as Thermodynamic Integration (TI), Free Energy Perturbation (FEP), MM/PB(GB)SA, Linear Interaction Energy. The MM/PB(GB)SA method is a compromise between accuracy and speed, and is widely used in the calculation of receptor-ligand binding free energy. The full name of this method is Molecular Mechanics/Poisson Boltzmann (Generalized Born) Surface Area. As the name suggests, this method splits the binding free energy into molecular mechanics terms and solvation energies to be calculated separately. The basic principle is to calculate the difference between the binding free energy of two solvated molecules in the bound and unbound states or to compare the free energy of different solvation conformations of the same molecule. Taking the binding of receptor protein and small ligand molecule as an example, the calculation principle will be roughly described.
For GROMACS users, there has not been a convenient tool for calculating MM-PBSA free energy for a long time. Since there is a script mmpbsa.py in AMBER to calculate MM-PBSA, the previous GROMACS users often first convert the obtained trajectory to AMBER format, and then use AMBER tools for free energy calculations. With the gradual increase of GROAMCS users, at present, there are already two tools that can directly use GROMACS tools for MM-PBSA free energy calculations.
GMXPBSAtool
Version 2.0: GMXPBSA 2.0: A GROMACS tool to perform MM/PBSA and computational alanine scanning
Version 2.1: GMXPBSA 2.1: A GROMACS tool to perform MM/PBSA and computational alanine scanning
Related papers: Exploring PHD Fingers and H3K4me0 Interactions with Molecular Dynamics Simulations and Binding Free Energy Calculations: AIRE-PHD1, a Comparative Study
Download link: http://gdriv.es/gmxpbsa
g_mmpbsa
The binding free energy calculated by g_mmpbsa consists of the following parts:
Part MM: energy in molecular mechanics, energy in vacuum, van der Waals and Coulomb electrostatic interaction energy
Part PB: Polarization solvation energy
SA part: non-polar solvation energy, can be calculated using three different methods SASA, SAV and Weeks-Chandler-Andersen (WCA)
When calculating, first use g_mmpbsa to calculate each part separately. After all calculations are completed, run the 3 scripts attached to the program in order for post-processing: MmPbSaStat.py calculates the average binding energy and its standard deviation/error, and MmPbSaDecomp.py performs energy decomposition, energy2fbc converts energy into a temperature factor. The specific process is as follows:
Calculate the SA part:
g_mmpbsa -f traj.xtc -s topol.tpr -n input.ndx -i apolar_sasa.mdp -nomme -pbsa -decomp -apol sasa.xvg -apcon sasa_contrib.dat
Use SAV model:
g_mmpbsa -f traj.xtc -s topol.tpr -n input.ndx -i apolar_sav.mdp -nomme -pbsa -decomp -apol sav.xvg -apcon sav_contrib.dat
Calculate the average binding energy:
Use bootstrap method: python MmPbSaStat.py -bs -nbs 2000 -m energy_MM.xvg -p polar.xvg -a apolar.xvg (or sasa.xvg)
Energy breakdown:
python MmPbSaDecomp.py -m contrib_MM.dat -p contrib_pol.dat -a contrib_apol.dat -bs -nbs 2000 -ct 999 -o final-contrib-energy.dat -om energymapin.dat
You must write -m, -p, -a in order. Metafile.dat file can also be used, where the content is the path of 3 files (must follow the order of MM, polar, apolar). The energy we need is written in Here. The final-contrib-energy.dat file is the specific content of energy decomposition, and another output file energymapin.dat is prepared for the next command.
Convert to temperature factor
python energy2fbc -s md.tpr -i energymapin.dat -n input.ndx -c complex.pdb -s1 s1.pdb -s2 s2.pdb
-s1, -s2 are for chain display, the sum of the number of atoms of s1 and s2 must be equal to the total number of atoms in the system in the previous calculation. The script will color the residues according to the temperature factor (b-factor), in PyMol Choose to color by b-factor to see the effect.
* For Research Use Only.
Related Services
