#はじめに

下記のチュートリアルを実行して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)を加えています。

#おわりに

以上、チュートリアルを実行してみましたが、

タンパク質-リガンド間の拘束の設定
リガンド単体の拘束の自由エネルギー計算

を自分で設定・計算する際にどうしたらいいかが疑問点として残りました。
もしわかる方いらっしゃったら教えていただけると嬉しいです。

あと、レプリカ交換と組み合わせて精度を上げるのも（マシンパワーがあれば・・・）やってみたいですね。