CBI学会2017年大会

GENESISチュートリアル２
自由エネルギー計算 (実習)
理化学研究所 計算科学研究機構
粒子系生物物理研究チーム
講演者：神谷基司
motoshi.kamiya@riken.jp
2017年10月5日 タワーホール船堀（東京）1今日の内容 15:45 – 17:00 GENESISでの自由エネルギー計算••••

自由エネルギー計算とは?
アンブレラサンプリング法について
レプリカ交換アンブレラサンプリング法(実習)
WHAM法による自由エネルギー計算(実習)2アンブレラサンプリング法

3

自由エネルギー計算

F  E  TS……とかいう式は特に意識する必要ありません(してもいいです)

自由エネルギーは(意外かもしれませんが)直観的に理解できる量でもある― 結合した状態と結合していない状態の比率

p (確率分布)― 何らかの反応座標軸の分布FkT log p F exp  
 kT 

反応座標(原子間距離等)

反応座標(原子間距離等)

自由エネルギーはむしろ現象論的な観測量(っぽいもの)と直接的に関係する(エントロピーはもちろんエネルギーよりも直観的に理解しやすい)4

色々な自由エネルギー計算方法(より正確には相対的自由エネルギーの計算方法)
 自由エネルギー摂動法(free energy perturbation)― 仮想的な状態を間に挿入し、状態間の自由エネルギー差を計算する― 粒子を生成消滅するような計算も可能― 結合自由エネルギー計算に良く用いられる(REST/FEP等含む)

 アンブレラサンプリング(umbrella sampling)― 拘束ポテンシャルを導入し、反応座標に沿った自由エネルギー変化を計算する(反応座標は粒子間距離や角度等を用いる)― 粒子を生成消滅させるような計算は通常はできない― 結合ポーズの高精度な推定、改善― リガンド結合の自由エネルギー障壁の高さ(~反応速度)推定

他にも様々な自由エネルギー計算法が存在する。
さらには、REST/FEP や REST/REUS 法のように複数組み合わせる場合もある。
5

p (確率分布)

アンブレラサンプリング法(1)
通常のMD計算でこのような分布が得られれば問題ない……が、普通はそうはうまくいかない

反応座標(原子間距離等)

F

MDトラジェクトリ

不安定領域の構造サンプルはとりにくい F 
安定性にFの差がある場合の存在比 = exp   kT 
F = 0 kcal/mol の場合 = 1 : 1

3 kcal/mol の場合 = 1 : exp( -3 / 0.6 ))= 1 : 0.0067
5 kcal/mol の場合 = 1 : 0.00024
10 kcal/mol の場合 = 1 : 0.000000058(※ 300 K)

どうにかしたい……

アンブレラサンプリング法(2)①(通常直接計算できない)

②

③

④
アンブレラ
ポテンシャル+真の値
反応座標(原子間距離等)①
真の値+アンブレラ
ポテンシャル

真の値

②

③

④

それぞれのアンブレラ
ポテンシャル存在下で分布を計算する(ウィンドウと呼ばれる)

7

アンブレラサンプリング法(3)各ウィンドウで計算を行って分布(ヒストグラム)を取得するF

アンブレラ
ポテンシャルが存在しない状態の自由エネルギー(真の値)

WHAM(最尤推定法の一つ)
で結果をまとめる(アンブレラポテンシャルも除去する)
Weighted Histogram
Analysis Method8アンブレラサンプリング法(4)分布1分布2 分布4
分布3
OK

分布5

NG

分布にはある程度重なりがなければならない- アンブレラポテンシャルの場所や形の調整- 計算を長くする
反応座標以外の方向「も」重なりが必須

反応座標ではない軸

反応座標(原子間距離等)

反応座標が状態量である必要性= 反応座標がわかれば分布が一意に決まる(ものすごく厳密な表現では)この条件が満たされなくてもWHAMを使えば自由エネルギーは計算できてしまう(正しい理由は無い)

NGOK
元々ダメ
反応座標(原子間距離等)条件をきちんと満たした正しい計算を行うことは専門家であっても極めて難しい(用いる手法によって色々と例外有り)(他にも気にするべきことはたくさん有り)レプリカ交換アンブレラサンプリング法10

レプリカ交換アンブレラサンプリング法(1)
アンブレラ法の精度を向上させる方法の一つ(計算コストは多少犠牲になる)
計算の途中でアンブレラ
ポテンシャルを交換する

交換！

アンブレラ
ポテンシャル

真の値
真の値+アンブレラポテンシャル
アンブレラポテンシャルの数だけ同時に計算を行う
交換を伴うために全ての計算が実質的に同期されていなければならない(通常のアンブレラサンプリングの場合には計算を個別に行える)(トータルの計算量自体は大きく増えているわけではない)
Sugita, Kitao, Okamoto, JCP 113, 6042 (2000).

11アンブレラサンプリング
NGOK
元々ダメ
反応座標(原子間距離等)

レプリカ交換アンブレラ
サンプリング(REUS)

反応座標ではない軸

反応座標ではない軸

レプリカ交換アンブレラサンプリング法(2)OK

やっぱりダメ反応座標(原子間距離等)交換によってサンプリング空間が拡大する
前述のめんどくさい問題の改善が期待できる
様々な安定構造が得られるようになるかも
交換は詳細つり合いを満たす方式で行う
無条件に交換できるわけではない
反応座標方向の分布に重なりがあり、交換が行われていれば良いサンプリングができると期待できるOK

やっぱりNG

反応座標(原子間距離等)

レプリカ交換アンブレラサンプリング法(3)
交換

レプリカ1

レプリカ2

交換

一つのレプリカ(トラジェクトリ)が複数のアンブレラポテンシャルを経験する

レプリカ3(色はアンブレラポテンシャルを意味する)
時間
通常のアンブレラ法を使うよりも広い範囲の構造探索が可能
リガンド結合の場合では、一旦脱離、再び結合するような構造変化が探索できる(通常のアンブレラ法や自由エネルギー摂動法では難しい)実習(レプリカ交換アンブレラサンプリング)

14

計算の流れ サンプル系の紹介
 反応座標の設定
 (システムのセットアップ(省略)) レプリカ交換アンブレラサンプリング(平衡化) レプリカ交換アンブレラサンプリング(サンプリング) トラジェクトリの変換
 トラジェクトリから距離の計算
 自由エネルギー計算
 (時間があれば)より長いトラジェクトリを使った解析(赤字部分が実習)
GENESIS ウェブサイト(ただし英語)にもレプリカ交換アンブレラサンプリングのチュートリアルがあります
http://www.aics.riken.jp/labs/cbrt/tutorial/remd-tutorials/tutorial-2-2/
今回扱う系よりもシンプルな系でのサンプルです(計算ステップ数は多いので適宜調節してください)15

簡単なサンプル系
PDB ID: 2DYFFBP11/HYPAに含まれるWWドメインの一つ
Proリッチなペプチドと結合する

Pro

GSTAPPLPRという配列を持つペプチドが結合している
全体で40残基程度で非常に小さい(今回の設定では溶媒の水をいれても10,000原子程度)

Kato, Miyakawa, Tanokura, JBC 281, 40321 (2006).

16

結合の反応座標の設定(1)
Tyr14
His6

Trp25
２個のProが溝にはまり込むような形になっている
右: Tyr14, Ser23, Trp25 の間 => 結合サイト１
左: His6, Tyr14, Tyr16 の間

Pro
Tyr16

Ser23

=> 結合サイト２両方の結合サイトを独立に扱うこともできるが、今回は二つをまとめて扱う
Pro(右)-結合サイト１の距離 +Pro(左)-結合サイト２の距離を反応座標として用いる

反応座標(原子間距離等)
コレ独立に扱う場合に比べて計算量的に大きく有利
17
ただし、注意するべき点が多少増える

結合の反応座標設定(2)
Tyr14
His6Trp25

Pro
Tyr16Ser23

Tyr14 Cg, Trp25 Ce2, Ser23 Cb の重心とPro(右)の５員環炭素の重心の距離+->group2->group1

His6 Cg, Tyr14 Cz, Tyr16 Cg の重心とPro(左)の５員環炭素の重心の距離->group4->group3[RESTRAINTS]
nfunctions
function1
select_index1
constant1
reference1

アンブレラポテンシャル:= 1= DISTMASS= 1 2 3 4= …
= …

k  r  r0 constant[SELECTION]
group1 = rno:38 and (an:CG or an:CB or an:CD or an:CA)group2 = (rno:25 and an:CE2) or (rno:14 and an:CG) or ¥(rno:23 and an:CB)
group3 = rno:35 and (an:CG or an:CB or an:CD or an:CA)
group4 = (rno:6 and an:CG) or (rno:14 and an:CZ) or ¥(rno:16 and an:CG)

2reference(今回関わる原子は全て炭素なので重心でも中心でも同じ。DISTMASSのかわりにDISTも可)
18

システムセットアップ(実習は省略)(AMBER形式; 実習では時間の関係で省略)
NMR構造の一つを取り出す
タンパク質周囲に水を加えた構造データ(トポロジー)の作成

Pro

構造最適化
(contact_check は gpu 未対応なので注意
外すか cpu だけを使うようにする)MDによる平衡化(粒子数体積温度一定条件(NVT),粒子数圧力温度一定条件(NPT))(ここまでの平衡化は通常のMDと同様、インプットファイル等は置いてあるので参考にしてください)19

実習(レプリカ交換アンブレラサンプリング；平衡化)計算の実行(spdyn)

ログの確認

時間等の制限により
今回は 4 レプリカだけ% cd ~/Tutorial_2/3_reus% sh run_reus_eq.sh &メインのログファイル(インプットの処理)% less reus0.log
各レプリカのログファイル(一番目のレプリカ)% less reus0_1.log

初期構造は４つのレプリカ全てが同じ構造を使う

reus0.log (抜粋):

reus0_1.log (抜粋):[STEP4] Compute Single Point Energy for Molecules

INFO:
STEP
TIME
TOTAL_ENE …
--------------- --------------- --------------- …
INFO:
100
0.1000-25471.7400 …

STEP
BOND
ANGLE …
--------------- --------------- --------------- …
0122.4199
314.7267 …
[STEP5] Perform Replica-Exchange MD Simulation[STEP6] Deallocate Arrays
Output_Time> Averaged timer profile (Min, Max)
total time=
98.154
setup
=
5.117

INFO:

200

0.2000-25481.7864 …

INFO:

300

0.3000-25497.9852 …

INFO:

400

0.4000-25498.1773 …

INFO:

500

0.5000-25473.3131 …

INFO:

600

0.6000-25424.1691 …

INFO:

700

0.7000-25450.2135 …

INFO:

800

0.8000-25471.3453 …

20

実習(レプリカ交換アンブレラサンプリング；平衡化)
reus0.inp: (spdyn用インプットファイルの抜粋と説明)[OUTPUT]
logfilerstfiledcdfile= reus0_{}.log= reus0_{}.rst= reus0_{}.dcd# -> {} には自動的にレプリカIDが入る# ファイルパスの ~ (チルダ)はintel fortran# ならホームディレクトリに置換されるが、# gfortran では置換されない(バージョン依存あるかも)[REMD]
dimension= 1exchange_period = 0
iseed= 4398237# 交換頻度; ０の場合は交換しない# 交換につかう乱数の種

type1
nreplica1
rest_function1= RESTRAINT= 4= 1# レプリカ交換アンブレラサンプリング# 4 レプリカを用いる# 下の 1 番の restraint を用いる

[RESTRAINTS]
nfunctions
function1
select_index1
constant1
reference1= 1= DISTMASS# 今回は DIST でも問題無い(はず)= 1 2 3 4
# [SELECTION]のgroupに対応= 1.0 1.0 1.0 1.0 # 力の定数 k; レプリカ数だけ並べる= 6.0 7.0 8.0 9.0 # ポテンシャルの中心 r0; レプリカ数並べる

レプリカあたりのプロセス数は自動で設定される(合計/レプリカ数; 今回 16 / 4)
プロセス数は必ず均等。各レプリカには最低４プロセス必要
21(OpenMPは適宜) (GPUはあってもなくても可)実習(レプリカ交換アンブレラサンプリング)計算の実行(spdyn)
reus1_1.rem (抜粋):% cd ~/Tutorial_2/3_reus% sh run_reus_prod.sh &ログの確認
メインのログファイル(インプット処理、レプリカ交換確率表示等)% less reus1.log
各レプリカのログファイル(一番目のレプリカ)% less reus1_1.logレプリカの状態(一番目のレプリカ)% less reus1_1.rem[RESTRAINTS]
nfunctions
function1
select_index1
constant1
reference1= 1= DISTMASS= 1 2 3 4
= 1.0 1.0 1.0 1.0= 6.0 7.0 8.0 9.0

12

340
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
11000
12000
13000
14000
15000
16000
17000
18000
19000

111
1
1
1
1
11
1
1
1
1
1
1
122
2
3

reus1_2.rem (抜粋):
01000
2000
3000
4000
5000
6000
7000
8000
9000
10000
11000
12000
13000
14000
15000
16000
17000
18000
19000

22222223333333333332ステップ数
22

RESTRAINTS の順番に対応実習(レプリカ交換アンブレラサンプリング)
reus1.inp: (抜粋)[OUTPUT]
logfileremfilerstfiledcdfile= reus1_{}.log= reus1_{}.rem= reus1_{}.rst= reus1_{}.dcd[REMD]
dimension= 1exchange_period = 1000
iseed
= 4398237# -> {} には自動的にレプリカIDが入る# 現在のパラメータを表すトラジェクトリ# 交換頻度; 1000 ステップに一回# 交換につかう乱数の種# 継続計算時は無視されます# レプリカ交換アンブレラサンプリング# 4 レプリカを用いる# 下の 1 番の restraint を用いる

type1
nreplica1
rest_function1= RESTRAINT= 4= 1[RESTRAINTS]
nfunctions
function1
select_index1
constant1
reference1= 1= DISTMASS= 1 2 3 4= 1.0 1.0 1.0 1.0 # 力の定数 k; レプリカ数だけ並べる
23= 6.0 7.0 8.0 9.0 # ポテンシャルの中心 r0; レプリカ数並べる

実習(レプリカ交換アンブレラサンプリング)
reus1.log (抜粋):
REMD> Step:
Replica
123449000
Dimension:
ExchangeTrial
1 >
0N
3 >
2A
2 >
3A
4 >
0N

Parameter:
RepIDtoParmID:
ParmIDtoRepID:
REMD> Step:
Replica
1234

1
1
1

12
2
250000Dimension:
ExchangeTrial
1 >
2A
2 >
1A
3 >
4R
4 >
3R

Parameter:
RepIDtoParmID:
ParmIDtoRepID:2223331111

ExchangePattern:
AcceptanceRatio
0 /
011 /
2511 /
250 /0After12341Before
1234After
2134444

ExchangePattern:
AcceptanceRatio
10 /2510 /
250 /
250 /
25
3332Before
1324444平衡化の時とは異なり、レプリカの交換確率が表示される

実習(レプリカ交換アンブレラサンプリング)
パターン1 パターン2

交換試行の結果
A: 交換成功(Accepted)
R: 交換失敗(Rejected)
N: 何もしない
REMD> Step:
Replica1234

交換のパターン

49000
Dimension:
ExchangeTrial
1 >
0N
3 >
2A2 >
3A4 >0N

Parameter:
RepIDtoParmID:
ParmIDtoRepID:

現在のアンブレラ
ポテンシャルのID

1
1
1

1222

ExchangePattern:
AcceptanceRatio
0 /
011 /
2511 /
250 /0333

交換のターゲット(対称)0の場合は交換しない

443322

1

12
Before1324(交互に行う)After
1234444

交換確率を表す(成功回数 / 試行回数)
計算をリスタートした際には、この回数も引き継がれる実習(トラジェクトリの変換; remd_convert)
レプリカ交換アンブレラサンプリング計算が終わった段階では、各トラジェクトリは色々なアンブレラポテンシャルでのデータが混じっている構造の再配列を行う（同じアンブレラポテンシャルのデータを集める）
水の部分を取り除いて取り扱いを楽にしておく
計算の実行 (remd_convert)% cd ~/Tutorial_2/4_conv% sh conv.sh &ログの確認
ログ確認% less conv1.log
トラジェクトリのファイルサイズ% ls –shconv1.log (抜粋):MD step:
16000 Replica:
1
Open trajectory file: "../reus/reus1_3.dcd"
MD step:
17000 Replica:3
MD step:
18000 Replica:3
MD step:
19000 Replica:3
MD step:
20000 Replica:3
MD step:
21000 Replica:3
MD step:
22000 Replica:3
MD step:
23000 Replica:3
MD step:
24000 Replica:3
MD step:
25000 Replica:
1
MD step:
26000 Replica:
1
Open trajectory file: "../reus/reus1_2.dcd"
MD step:
27000 Replica:
2MD step:
28000 Replica:2実習(トラジェクトリの変換; remd_convert)
conv1.inp:[INPUT]
prmtopfile
ambcrdfile
dcdfileremfile= ../0_str/ww.prmtop= ../0_str/ww.crd= ../3_reus/reus1_{}.dcd= ../3_reus/reus1_{}.rem[OUTPUT]
trjfile= reus1_{}.dcd# 出力するファイル名。{}を指定するのを忘れずに[SELECTION]
group1= ai:1-608# ここで指定した原子だけを書き出す# 下のtrjout_atomで実際に指定# 先頭から順番に取り出すのではなく、途中の原子だけを取り出すことも# できるが、以後の解析で面倒なことが起こる可能性があるので# できれば避けた方が良いかもしれない[OPTION]
check_only
convert_type= NO= PARAMETER

dcd_md_period = 100
trjout_format = DCD
trjout_type= COOR
trjout_atom
pbc_correct= 1= MOLECULE# トラジェクトリーファイル(入力)# レプリカの状態# パラメータでそろえる。REUSの場合は# 同じアンブレラポテンシャルのデータまとめる# crdout_periodに相当# remfileと合わせるために必須# 出力時に周期境界の設定を除去する# 入力は周期境界有りで読まれる# 原子指定はここで行う# 周期境界で分子がちぎれないようにする処理27

実習(距離の計算; trj_analysis)
自由エネルギー計算のために反応座標の値(重心間距離の和)を各アンブレラポテンシャルのデータについて計算する(今回の場合はエネルギーの出力から直接計算できなくもない)計算の実行(trj_analysis)(各レプリカについて計算)% cd ~/Tutorial_2/5_trjana% sh trjana.sh &ログの確認
ログ確認(1番目)% less trjana1.log
データ確認(1番目)% less dist1-1_1.dis% less distsum1-1_1.dis

dist1-1_1.dis:
12
3
4
56
7
8
9
10
11
12
13
14
15
16
17

distsum1-1_1.dis:
2.235
2.115
1.951
2.169
2.016
2.373
2.348
2.410
2.257
2.270
2.119
2.007
1.964
2.077
1.869
1.876
1.843

3.510
3.624
3.744
3.696
3.479
3.329
3.287
3.265
3.383
3.471
3.709
3.696
3.630
3.633
3.737
3.985
3.899

1 5.745
2 5.739
3 5.695
4 5.865
5 5.495
6 5.702
7 5.635
8 5.675
9 5.6410 5.741
11 5.828
12 5.703
13 5.594
14 5.71
15 5.606
16 5.861
17 5.742

実習(距離の計算; trj_analysis)trjana1.inp: (通常のアンブレラサンプリングの場合でも同じインプットを使えます)[INPUT]
prmtopfile= ../0_str/ww.prmtop

ambcrdfile= ../0_str/ww.crd[OUTPUT]
comdisfile= dist1-1_1.dis[TRAJECTORY]
trjfile1= ../conv/reus1_1.dcd#trjfile2= ../conv/reus2_1.dcd
md_step1= 50000
mdout_period1 = 100
ana_period1= 100
repeat1#repeat1
trj_format
trj_type
trj_natom= 1= 2= DCD= COOR= 608# 質量を使うタイプの計算では必須# ambcrdやpdbファイルだけでは代用できない# comdisfile は重心間距離出力用ファイル# 他の値を計算する場合は違う名前を使う# 二つ同時に扱う場合の例# REUS計算時のステップ数# REUS計算時のcrdout_periodに合わせる# ステップ数の計算はmd_stepに準拠する# mdout_period=ana_periodで全て計算するようになる# 複数同時に処理する場合には1以外の数字# 二つ目のファイルがある場合# 変換前のトラジェクトリを使うならCOOR+BOX# トラジェクトリ中の原子数指定# [INPUT]での入力から変えるために必要[SELECTION](REUS計算の時と同じ。スペースが足りなかったので省略)# 今回は先頭から順番に取り出しているのでこれでOK[OPTION]
check_only= NO
com_distance1 = 1 2com_distance2 = 3 4# analysis_tool での共通パラメータ# SELECTION で指定したグループの# 重心距離を計算する(距離の和は計算できない)実習(自由エネルギー計算; wham_analysis)

WHAM(Weighted Histogram Analysis Method)法で自由エネルギー計算% cd ~/Tutorial_2/6_wham% sh wham.sh &ログの確認
ログ確認% less wham.log
自由エネルギーデータ% less output.pmf

自由エネルギー (kcal/mol)

計算の実行(wham_analysis)

NMR構造での値はこの辺結合サイト距離の和 (Å)(今回は計算が短すぎるので精度的には非常にアレです)実習(自由エネルギー計算; wham_analysis)wham.inp:[INPUT]
cvfile= ../5_trjana/distsum1-1_{}.dis# {} は REUS(MD)と同じ意味# ここで与えた軸に沿った# 自由エネルギープロファイルが# 計算される[OUTPUT]
pmffile= output.pmf# 自由エネルギーの出力ファイル[WHAM]
dimension
nblocks
temperature
tolerance
rest_function1
grids1= 1= 1= 300.0= 10E-08= 1= 4.0 22.0 91[RESTRAINTS]
constant1
reference1= 1.0 1.0 1.0 1.0= 6.0 7.0 8.0 9.0# どの拘束関数を使うかの指定# 自由エネルギープロファイルの# グリッドの定義# (最小値) (最大値) (グリッド数)# REUSの時と同じ値をいれる(MBAR(Multistate Bennett Acceptance Ratio)法を用いた場合でも同様の計算は可能)

まとめ

GENESIS でレプリカ交換アンブレラサンプリング
自由エネルギーを計算するまでの流れを簡単に実習

GENESIS ウェブサイト(※英語)
日本語での質問は可能になっています(要ユーザ登録)
http://www.aics.riken.jp/labs/cbrt/
もう少し簡単な系でのチュートリアル(※英語)http://www.aics.riken.jp/labs/cbrt/tutorial/remd-tutorials/tutorial-2-2/

おまけ（もうすこし大きなデータを使った場合）
Tutorial_2/x_appendix 以下にもう少し大きな計算のデータがあります。
これを使って自由エネルギー計算を試してみましょう。
(使用した中間ファイル等もある程度置いてあるので参考にしてください)
距離範囲 6 Å ~ 21 Å まで(16レプリカ), レプリカ当たり 50 ns の計算

自由エネルギー (kcal/mol)

距離の和のファイル(Tutorial_2/x_appendix/trjana)
2.5 ns までのデータ: distsum1-5_1.dis … distsum1-5_16.dis
10 ns までのデータ: distsum1-20_1.dis … distsum1-20_16.dis
25 ns までのデータ: distsum1-50_1.dis … distsum1-50_16.dis
50 ns 全体のデータ: distsum1-100_1.dis … distsum1-100_16.dis

結合サイト距離の和 (Å)

距離長いところの構造