Skip to content

Latest commit

 

History

History
108 lines (87 loc) · 5.54 KB

README.md

File metadata and controls

108 lines (87 loc) · 5.54 KB

FEP-ABFE: free-energy perturbation with absolute binding free energy calculation

What's this?

FEP-ABFE calculates the ligand - the receptor binding free energy. Suppose the value obtained by the pipeline being ΔG. The dissociation constant K_d is then obtained by:

K_d = C_0 exp(ΔG/RT)

where C_0 = 1 mol / L and RT is the gasconstant * temperature (K). Typically K_d is around fmol/L to μmol/L order.

Prepration

Before running the script, prepare mdtraj and pyedr package for python3, in the system you run fepsuite.

% pip3 install cython mdtraj pyedr --user   # cython is necessary for mdtraj
% python3 -c "import mdtraj"   # check whether mdtraj works

Also, you will need GROMACS 2022.5 or later (2022.4 or -2021 does not work correctly).

Finally, if your ligand is not neutral (charged), you will need to install APBS. If you are using Ubuntu and you are an administrator you can install by:

% sudo apt install apbs

For other cases, manual installation of APBS is difficult (due to bunch of dependencies). I recommend using binary releases of APBS if you can use.

Calculating the free energy

Then, follow the next order to run the FEP-ABFE.

  1. Prepare topology with protein + ligand + water + ions of you want to calculate. If you are starting from the structure generated by docking, please see the latter section before proceeding the calculation. We recommend to use rhombic dodecahedron box for the solvation. The system must be neutralized.
  2. Set topology file name as topol_ionized.top and coorinate file as conf_ionized.pdb or conf_ionized.gro. (You can use symlink for that).
  3. Prepare a directory with two depths, e.g. abfecalc/mol1. Put topology and structure files under abfecalc/mol1.
  4. Clone this git repository to somewhere, e.g. $HOME/repos/fep-suite.
  5. Copy an entire contents under rundir-template to abfecalc, by e.g.:
cd $HOME/abfecalc
cp -r $HOME/repos/fep-suite/abfe/rundir_template/* .
  1. Modify run.zsh's FEPSUITE_ROOT to git repository root such as $HOME/repos/fep-suite. Change JOBSYSTEM to your system.
  2. Modify para_conf.zsh according to your calculation environment.
  3. Chdir to abfecalc directory, e.g. cd $HOME/abfecalc.
  4. Run the submission script with secondary directory name (in this case mol1) and all, e.g. ./run.zsh mol1 all. The script will submit all necessary jobs.
  5. Wait until all calculations finish.
  6. If everything works fine, output should be generated at abfecalc/mol1/result.txt.

How to interpret the result

An example result look like:

annihilation-complex -5.449 0.837
annihilation-lig -4.264 0.253
charging-complex -26.856 0.268
charging-lig 23.439 0.079
lr-annihilation-complex 6.986 0.054
lr-annihilation-lig 0.095 0.002
lr-complex -7.244 0.062
lr-lig -0.248 0.002
restrain -0.872 0.010
restrain-analytical 2.246 0.000
----
annihilation -9.713 0.874
charging -3.418 0.279
long-range-correction -0.411 0.083
restrain 1.374 0.010
----
total -12.168 0.921

Each line comprises of name, mean value and standard deviation of each value in kcal/mol. First section represents contributions from individual terms of the thermodynamic cycle. These signs of these values are set so that the total value, in the last section, is just the summation of these values. Second section represents subtotals of the first section, separated by major contributin factors. The total value represents the final free-energy difference. In the above example, it's -12.168 kcal/mol +/- 0.921 kcal/mol.

Starting from docking results

Typically, poses obtained from docking simulation are not stable. We recommend the following way to sample configurations. (TODO)

Bugs and caveats

Currently charged correction for charged lingands are in "testing" phase. If you find any strange baheviour please let me know.

Ligand with a large diameter requires long rlist length and therefore consumes much calculation time.

We assume the ratio of symmetry number of the molecule to be 1. This will not be a big problem unless you are running simulations for molecules with extremely high symmetry, e.g. benzene. In case you are predicting delta-G of such molecules, you need to add log-symmetry terms (see e.g. Boresch paper below for details). If anyone know a good (and practical) algorithm for getting symmetry number please let me know.

Current script does not work with MacOS (gsed instead of sed is needed)

References

Overall algorithm implemented in this pipeline is similar to:

Accurate calculation of the absolute free energy of binding for drug molecules
Matteo Aldeghi, Alexander Heifetz, Michael J. Bodkin, Stefan Knapp and Philip C. Biggin
Chem. Sci. 7, 207-218 (2016). doi:10.1039/c5sc02678d

We wrote the initial version based on SI of the following paper (but later switched to Aldegi et al. paper)

Identify potent SARS-CoV-2 main protease inhibitors via accelerated free energy perturbation-based virtual screening of existing drugs
Zhe Li, Xin Li, Yi-You Huang, Yaoxing Wu, Runduo Liu, Lingli Zhou, Yuxi Lin, Deyan Wu, Lei Zhang, Hao Liu, Ximing Xu, Kunqian Yu, Yuxia Zhang, Jun Cui, Chang-Guo Zhan, Xin Wang and Hai-Bin Luo
Proc. Natl. Acad. Sci. 117 (44), 27381-27387 (2020). doi:10.1073/pnas.2010470117

Core part of this program is absolute free energy correction, and is implemented based on

Absolute Binding Free Energies: A Quantitative Approach for Their Calculation
Stefan Boresch, Franz Tettinger, Martin Leitgeb and Martin Karplus
J. Phys. Chem. B 107, 9535-9551 (2003). doi:10.1021/jp0217839