GROMACSで結合自由エネルギーを計算する

#はじめに

下記のチュートリアルを実行してGROMACSで結合自由エネルギーを求めていきます。

http://www.alchemistry.org/wiki/Absolute_Binding_Free_Energy_-_Gromacs_2016


詳しいことはチュートリアルを読んでいただきたいですが、このようなスキームで結合エネルギーを求めていきます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/27003605-1610-f9b0-1efc-de0005d5fbaf.png)


#準備

インプットファイルは下記リンクから落としてきます。

http://www.alchemistry.org/wiki/images/1/17/InputFiles_ABFE_GMX2016.zip

zipを解凍します。

```bash
unzip InputFiles_ABFE_GMX2016.zip
```

リガンドだけの系と、リガンド＋タンパクの系があります。

![Screenshot from 2021-07-19 14-56-54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/b02fe2dd-0e6c-b24b-cd4d-00ecb456d4fd.png)

![Screenshot from 2021-07-20 10-39-29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/9a4bc399-b72b-8020-b433-e91f4b4bec01.png)



#複合体からリガンドをデカップリングする(状態F→E→D)

状態F→E→Dに対応する計算、つまりリガンドに位置の拘束をかけて、リガンドと周囲との静電相互作用およびvdw相互作用を無くしていきます。

拘束(bonded),静電相互作用(coulomb),vdwの3つのラムダを使います。
以下のグラフのようにラムダを変化させます。

![Screenshot from 2021-08-05 16-07-51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/cdde362e-e0b4-66d3-7167-5dc4ec5d730f.png)

ちなみに位置の拘束はcomplex.top中で

```python
[ 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のディレクトリに移動します。

```bash
cd complex
```

mdpファイルも既に用意されています。
あとはrun.shを実行するだけですが、そのままでエラーが出る場合や計算が遅すぎる場合は修正を加えましょう。

私の場合は、run.shの中のgmx_avxをgmxに置換、また並列数を指定するため、mdrunのオプションに-nt 4を指定しました。
また、そのまま実行するとPMEの計算にGPUが使われなかったので、下記のように変更しました。

```bash
#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
```
以下のコマンドを使って、まとめて置換しました。

```bash
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）



```bash
cd ligand
```

あとはcomplexのときと同じようにrun.shやMDPに修正をしてrun.shを実行します。
これで必要な計算は終了です。


#結果解析

ligand, complexそれぞれで実行する

##xvgファイルを集める

まず、リンクからCp_xvg.shを落として実行し、各ラムダで出力されたxvgファイルを集めてきます。

http://www.alchemistry.org/wiki/images/7/71/Cp_xvg.sh.zip

```
./Cp_xvg.sh
```

dHdl_filesディレクトリができていて、dhdlファイルがコピーされている。


##自由エネルギーの計算

下記からalchemical-analysis-master.zipをダウンロードする。

https://github.com/MobleyLab/alchemical-analysis

残念ながらpython2系で書かれていました。
私はAnacondaで仮想環境を構築しましたが、お好きな方法で使ってください。

```bash
conda create -n py27 python=2.7
conda activate py27
conda install numpy
conda install matplotlib
conda install -c omnia pymbar
```

alchemical_analysis.pyを実行していきます。

```bash
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
でした


他の出力ファイルもなかなか面白いです。

![Screenshot from 2021-07-24 17-03-53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/ee702e1d-85d2-a779-7407-2e5befe93edc.png)


![Screenshot from 2021-07-24 17-04-43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/ba8946ac-c90f-81e5-298b-274cd60399d0.png)


#リガンドの拘束に関わるエネルギーを求める（状態B→C）

状態B→Cに移るエネルギーは先に行った計算には入っていませんでした。
これはBoreschらが提唱した式により解析的に求めることができます。

状態C→Bの自由エネルギーが
-6.906 kJ/mol
となるそうです。

restraint.pyというもので出せるそうなのですが、どこにあるのかわかりませんでした。
検索した感じですと、これ↓でしょうか。もしわかる方いたら教えてください。

https://github.com/choderalab/yank



#最終結果

さて、これまでに行った計算結果をまとめて、結合自由エネルギーを出していきましょう。
分散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)を加えています。

#おわりに

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

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

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

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