#はじめに
下記のチュートリアルを実行してGROMACSで結合自由エネルギーを求めていきます。
詳しいことはチュートリアルを読んでいただきたいですが、このようなスキームで結合エネルギーを求めていきます。
#準備
インプットファイルは下記リンクから落としてきます。
zipを解凍します。
unzip InputFiles_ABFE_GMX2016.zip
リガンドだけの系と、リガンド＋タンパクの系があります。
#複合体からリガンドをデカップリングする(状態F→E→D)
状態F→E→Dに対応する計算、つまりリガンドに位置の拘束をかけて、リガンドと周囲との静電相互作用およびvdw相互作用を無くしていきます。
拘束(bonded),静電相互作用(coulomb),vdwの3つのラムダを使います。
以下のグラフのようにラムダを変化させます。
ちなみに位置の拘束はcomplex.top中で
[ intermolecular_interactions]
[ bonds ]
; ai aj type bA kA bB kB
1391 2615 6 0.654 0.0 0.654 4184.0
[ angles ]
; ai aj ak type thA fcA thB fcB
1393 1391 2615 1 88.8 0.0 88.8 41.84
1391 2615 2614 1 32.9 0.0 32.9 41.84
[ dihedrals ]
; ai aj ak al type thA fcA thB fcB
1410 1393 1391 2615 2 -159.7 0.0 -159.7 41.84
1393 1391 2615 2614 2 122.6 0.0 122.6 41.84
1391 2615 2614 2610 2 12.8 0.0 12.8 41.84
として定義されています。
complexのディレクトリに移動します。
cd complex
mdpファイルも既に用意されています。
あとはrun.shを実行するだけですが、そのままでエラーが出る場合や計算が遅すぎる場合は修正を加えましょう。
私の場合は、run.shの中のgmx_avxをgmxに置換、また並列数を指定するため、mdrunのオプションに-nt 4を指定しました。
また、そのまま実行するとPMEの計算にGPUが使われなかったので、下記のように変更しました。
#before
vdw-type = PME
integrator = sd
pme-order = 6
nstlist = 10
fourierspacing = 0.10
#after
vdw-type = Cut-off
integrator = md
pme-order = 4
nstlist = 40
fourierspacing = 0.16
以下のコマンドを使って、まとめて置換しました。
cd ./MDP/PROD
#ディレクトリ中の全mdpファイルについて、置換を行う。必要に応じてNPTやNVTも置換してください
ls ./*.mdp | xargs sed -i.bak -e 's/vdw-type = PME/vdw-type = cut-off/g'
ls ./*.mdp | xargs sed -i.bak -e 's/integrator = sd/integrator = md/g'
ls ./*.mdp | xargs sed -i.bak -e 's/pme-order = 6/pme-order = 4/g'
ls ./*.mdp | xargs sed -i.bak -e 's/nstlist = 10/nstlist = 40/g'
ls ./*.mdp | xargs sed -i.bak -e 's/fourierspacing = 0.10/fourierspacing = 0.16/g'
ls ./*.mdp | xargs sed -i.bak -e 's/constraints = all-bonds/constraints = h-bonds/g'
cd ../../
あとはrun.shを実行しましょう
chmod u+x run.sh
./run.sh
内容としては、エネルギー最小化、NVT、NPT、Production runの計算を各ラムダについて順次実行してくれます。
私の自宅PCで8時間程度はかかったと思います。
完了したらリガンドのほうの計算を行いましょう。
#溶媒からリガンドをデカップリングする（状態A→B）
cd ligand
あとはcomplexのときと同じようにrun.shやMDPに修正をしてrun.shを実行します。
これで必要な計算は終了です。
#結果解析
ligand, complexそれぞれで実行する
##xvgファイルを集める
まず、リンクからCp_xvg.shを落として実行し、各ラムダで出力されたxvgファイルを集めてきます。
./Cp_xvg.sh
dHdl_filesディレクトリができていて、dhdlファイルがコピーされている。
##自由エネルギーの計算
下記からalchemical-analysis-master.zipをダウンロードする。
残念ながらpython2系で書かれていました。
私はAnacondaで仮想環境を構築しましたが、お好きな方法で使ってください。
conda create -n py27 python=2.7
conda activate py27
conda install numpy
conda install matplotlib
conda install -c omnia pymbar
alchemical_analysis.pyを実行していきます。
unzip alchemical-analysis-master.zip
cd alchemical-analysis-master
python setup.py install
cd ../complex/dHdl_files
#xvgファイル中の温度設定がおかしかったので修正する
ls *.xvg | xargs sed -i.bak -e 's/T = 0 (K)/T = 300 (K)/g'
#自由エネルギー計算の実行。-dでxvgファイルがある場所の指定、-pでxvgファイル名を指定
python ../../alchemical-analysis-master/alchemical_analysis/alchemical_analysis.py -d . -p dhdl -t 300 -s 100 -u kcal -w -g
results.txtに結果が表示されます。
Complexのresult.txtは以下でした。
------------ --------------------- --------------------- --------------------- --------------------- --------------------- ---------------------
States TI (kcal/mol) TI-CUBIC (kcal/mol) DEXP (kcal/mol) IEXP (kcal/mol) BAR (kcal/mol) MBAR (kcal/mol)
------------ --------------------- --------------------- --------------------- --------------------- --------------------- ---------------------
0 -- 1 3.988 +- 0.019 3.976 +- 0.020 3.942 +- 0.046 3.952 +- 0.037 3.962 +- 0.019 3.962 +- 0.019
1 -- 2 2.529 +- 0.015 2.503 +- 0.017 2.512 +- 0.031 2.483 +- 0.031 2.504 +- 0.016 2.498 +- 0.015
2 -- 3 1.345 +- 0.012 1.323 +- 0.014 1.336 +- 0.022 1.297 +- 0.023 1.321 +- 0.012 1.324 +- 0.010
3 -- 4 0.403 +- 0.009 0.396 +- 0.011 0.411 +- 0.018 0.441 +- 0.020 0.403 +- 0.009 0.398 +- 0.006
4 -- 5 0.512 +- 0.003 0.512 +- 0.003 0.512 +- 0.006 0.514 +- 0.003 0.512 +- 0.003 0.515 +- 0.002
5 -- 6 0.499 +- 0.003 0.500 +- 0.003 0.498 +- 0.005 0.501 +- 0.004 0.500 +- 0.003 0.502 +- 0.002
6 -- 7 0.950 +- 0.006 0.951 +- 0.008 0.956 +- 0.014 0.962 +- 0.009 0.953 +- 0.006 0.956 +- 0.004
7 -- 8 0.872 +- 0.007 0.876 +- 0.008 0.861 +- 0.019 0.883 +- 0.008 0.876 +- 0.006 0.870 +- 0.005
8 -- 9 0.748 +- 0.008 0.751 +- 0.009 0.766 +- 0.031 0.751 +- 0.012 0.751 +- 0.008 0.747 +- 0.007
9 -- 10 0.560 +- 0.010 0.576 +- 0.012 0.582 +- 0.030 0.576 +- 0.014 0.571 +- 0.010 0.570 +- 0.009
10 -- 11 0.212 +- 0.015 0.226 +- 0.021 0.341 +- 0.030 0.206 +- 0.022 0.240 +- 0.015 0.245 +- 0.013
11 -- 12 -0.096 +- 0.011 -0.087 +- 0.012 -0.106 +- 0.023 -0.076 +- 0.015 -0.090 +- 0.010 -0.087 +- 0.009
12 -- 13 -0.348 +- 0.017 -0.333 +- 0.019 -0.309 +- 0.027 -0.334 +- 0.025 -0.330 +- 0.016 -0.330 +- 0.014
13 -- 14 -0.728 +- 0.025 -0.734 +- 0.028 -0.769 +- 0.054 -0.728 +- 0.044 -0.732 +- 0.025 -0.734 +- 0.021
14 -- 15 -1.032 +- 0.025 -1.058 +- 0.028 -1.074 +- 0.053 -1.040 +- 0.047 -1.069 +- 0.025 -1.093 +- 0.020
15 -- 16 -1.064 +- 0.017 -1.087 +- 0.019 -1.061 +- 0.026 -1.113 +- 0.026 -1.092 +- 0.016 -1.097 +- 0.011
16 -- 17 -0.837 +- 0.010 -0.847 +- 0.011 -0.840 +- 0.015 -0.848 +- 0.011 -0.844 +- 0.009 -0.831 +- 0.006
17 -- 18 -0.492 +- 0.005 -0.490 +- 0.006 -0.495 +- 0.007 -0.484 +- 0.007 -0.489 +- 0.005 -0.486 +- 0.004
18 -- 19 -0.154 +- 0.004 -0.151 +- 0.004 -0.149 +- 0.005 -0.152 +- 0.005 -0.151 +- 0.003 -0.152 +- 0.003
------------ --------------------- --------------------- --------------------- --------------------- --------------------- ---------------------
Coulomb: 8.266 +- 0.037 8.198 +- 0.038 8.202 +- 0.062 8.174 +- 0.057 8.190 +- 0.029 8.181 +- 0.035
vdWaals: -0.398 +- 0.070 -0.396 +- 0.071 -0.289 +- 0.106 -0.382 +- 0.083 -0.393 +- 0.050 -0.405 +- 0.066
TOTAL: 7.868 +- 0.079 7.802 +- 0.080 7.912 +- 0.123 7.791 +- 0.100 7.797 +- 0.058 7.775 +- 0.074
つまり、MBAR法で
23.827 +- 0.065 kJ/mol
でした
また、リガンドはMBAR法で
7.775 +- 0.074 kJ/mol
でした
他の出力ファイルもなかなか面白いです。
#リガンドの拘束に関わるエネルギーを求める（状態B→C）
状態B→Cに移るエネルギーは先に行った計算には入っていませんでした。
これはBoreschらが提唱した式により解析的に求めることができます。
状態C→Bの自由エネルギーが
-6.906 kJ/mol
となるそうです。
restraint.pyというもので出せるそうなのですが、どこにあるのかわかりませんでした。
検索した感じですと、これ↓でしょうか。もしわかる方いたら教えてください。
#最終結果
さて、これまでに行った計算結果をまとめて、結合自由エネルギーを出していきましょう。
分散sigma1とsigma2を足すと、(simga1^2 + sigma2^2)^1/2になることに注意すると、
(状態A→F)
=（状態A→B） - (状態C→B) ＋ （状態C→D） - （状態F→D）
= (7.775 +- 0.074) - (-6.906) + (0) - (23.827 +- 0.065) kJ/mol
= -9.146 +- 0.098 kJ/mol
チュートリアルの結果とずれてしまったのは、MDPファイルを変更したのがあるのかもしれません。
チュートリアルでは更にリガンドの対称性を考慮し、フェニル基の180度回転ポーズを想定してさらに-ktln(2)を加えています。
#おわりに
以上、チュートリアルを実行してみましたが、
- タンパク質-リガンド間の拘束の設定
- リガンド単体の拘束の自由エネルギー計算
を自分で設定・計算する際にどうしたらいいかが疑問点として残りました。
もしわかる方いらっしゃったら教えていただけると嬉しいです。
あと、レプリカ交換と組み合わせて精度を上げるのも（マシンパワーがあれば・・・）やってみたいですね。