蛋白質科学会アーカイブ, 14, e105 (2021)

GaREUS 法を用いた自由エネルギー地形の計算
理化学研究所・生命機能科学研究センター
尾嶋拓

Calculation of Free-Energy Landscape Using the GaREUS Method
RIKEN Center for Biosystems Dynamics ResearchHiraku Oshima（投稿日 2021/03/04、再投稿日 2021/06/16、受理日 2021/08/23）
キーワード：分子動力学法（Molecular Dynamics）、自由エネルギー計算法、拡張アンサンブル法、GaREUS 法、GENESIS
イントロダクション
分子動力学シミュレーション（Molecular Dynamics: MD）は、粒子間に働く相互作用および力を計算し、運動方程式を解くことで、粒子の運動を求める方法である。溶媒を含む全ての相互作用を考慮し、古典力学の範囲内で生体分子の動的構造を計算し、その計算結果を解析することで様々な物理量を求めることができる。タンパク質や核酸などの生体分子の動的構造や相互作用を原子解像度で調べる手法として生命科学において広く用いられている。特に、MD シミュレーションでは各状態での安定性または存在確率を求めることができ、これらは自由エネルギー地形と呼ばれている。自由エネルギー地形を詳細に分析することで、生体分子の構造変化や分子機構を理解することができる。
しかし、生体分子の自由エネルギー地形は多数の局所安定状態を持った複雑な形状をしており、限られた計算時間ではシミュレーションが局所安定状態に捕らわれてしまう。この問題を避けるために、拡張アンサンブル法などの効率的なサンプリング法が用いられている。拡張アンサンブル法では、エネルギーあるいは事前に定義された反応座標における一次元ランダムウォークを実現することで、状態空間を広く探索する。
様々な拡張アンサンブル法が提案されているが、計算コストを抑えつつ、タンパク質構造を効率的に探索する計算手法が望ましい。そのような計算手法として、我々はGaussian accelerated MD（GaMD）法（文献 1,2）と Replica-Exchange Umbrella
Sampling （ REUS ） 法 （ 文 献 3,4 ） の 2 つ の 拡 張 ア ン サ ン ブ ル 法 を 組 み 合 わ せ た
Gaussian accelerated Replica-Exchange Umbrella Sampling（GaREUS）法を提案した（文献 5）。GaREUS 法では、GaMD 法によってタンパク質の構造柔軟性を向上させ、同時に REUS 法によって反応座標空間を広く探索する。GaREUS 法で得られたトラジェクトリーから二段階の Reweighting を行うことで、正確な自由エネルギー地形を求めることが出来る。GaREUS 法は、理化学研究所で開発している MD シミュレーション
ソフトウェア「GENESIS」（文献 6,7）に実装されており、GPU ワークステーションからスーパーコンピュータまで様々な計算機を用いた計算が可能である。本稿ではGaREUS 法を用いた自由エネルギー地形の計算方法を説明する。

1

蛋白質科学会アーカイブ, 14, e105 (2021)

計算機・ソフトウェア
GaREUS 法の実行には MD シミュレーションソフトウェア GENESIS を用いる。
GENESIS はオープンソフトウェア（LGPLv3）として公開しており、無料で利用できる。
最新版は GENESIS のウェブサイト（https://www.r-ccs.riken.jp/labs/cbrt/）から入手可能である。
GENESIS をコンパイル・実行するには、OS として Linux、C コンパイラおよびFortran コンパイラとして gcc、gfortran または Intel Fortran/C++ compiler、開発環
境として GNU autoconf、並列計算ライブラリーとして OpenMPI、数学ライブラリーとして BLAS と LAPACK が必要である。GPU を利用する場合は CUDA ライブラリーも必要となる。詳細な必要条件やコンパイル方法は GENESIS のウェブサイトかマニュアルを参考にしてほしい。
GaREUS 法では複数のレプリカ（MD シミュレーション）を同時に実行する必要があるため、複数の計算ノードを持ったコンピュータクラスターかスーパーコンピュータを用いる。Intel 社、AMD 社、富士通（A64FX）の CPU を搭載しており、GPU を利用する
場合は、CUDA 対応の nVIDIA 社の GPU を搭載したものを用いる。例えば、10 レプリ
カの計算をする場合、8CPU コア以上を搭載した計算ノードが 10 個必要となる。複数ノードを用いた計算にはバッチジョブシステムの利用が必要だが、利用方法はソフトウェアや計算機システムごとに大きく異なるため、各システムのマニュアルを参照するか管理者に問い合わせてほしい。
計算手順
1）生体分子モデルのセットアップ
2）構造最適化と平衡化
3）GaMD パラメータの決定
4）本計算
5）解析

2蛋白質科学会アーカイブ, 14, e105 (2021)

実験の詳細
1）生体分子モデルのセットアップ
例として Alanine Tripeptide を用いる（図 1）。GaREUS の反応座標として、残基 1 の酸素と残基 3の水素の原子間距離を用いる（図 1の点線丸）。Alanine Tripeptideに水分子を加え、一辺 50 Å程度のボックスを作る。ここでは力場パ
1：Alanine Tripeptide
ラメータとして CHARMM 力場の 図
点線丸で囲った原子間の距離を反応座標として使用
C36 を用いる。タンパク質のセットアップについては、GENESIS のウェブサイト（https://www.r-ccs.riken.jp/labs/cbrt/）のチュートリアル（Level 1:Basic tutorials）および文献 8 で説明されているので、そちらを参照していただきたい。
2）構造最適化と平衡化
セットアップで得た構造は原子同士が衝突していたり、不安定な構造になっていたりすることがある。ポテンシャルエネルギーが極端に高すぎるため、そのまま MD 計算をしてしまうとシステムが崩壊してしまう。構造最適化によってポテンシャルエネルギーを下げ、原子の衝突を取り除く。このとき、対象としているタンパク質やリガンドの位置や距離に拘束をかけ、分子が動かないようにする。GENESIS では最急降下法で構造最適化を行う。
構造最適化が終わったシステムは分子に速度が与えられていない絶対零度の状態になっている。分子に初期速度を与え、システムの温度を標的温度までゆっくりと上昇させる。
GENESIS では Langevin 熱浴で等温等積（NVT）の MD 計算を行う。
次に、等温等圧（NPT）で MD 計算を行い、システムの体積を調整する。体積の時間変化を確認し、体積があまり変化しなくなるまでシミュレーションを行う。GENESIS では Bussi の方法で温度・圧力を調整する。
最後に、本計算と同じ熱力学条件で平衡化を行う。平衡化に必要なステップ数はシステ
ムによって異なるため、温度、体積、ポテンシャルエネルギーの時間変化がドリフトしていないか注意する。構造のトラジェクトリーを見て、構造が壊れていないかも確認する。
Root-Mean-Square Displacement （ RMSD ） や Root-Mean-Square Fluctuation（RMSF）等で時間変化を確認し、あまり変化しなくなったところで平衡化を終える。構造最適化および平衡化についても GENESIS のウェブサイトのチュートリアル（Level 1:Basic tutorials）および文献 8 で GENESIS のインプットファイルも含めて詳細な説明されているので、そちらを参照していただきたい。
3）GaMD パラメータの決定

3

蛋白質科学会アーカイブ, 14, e105 (2021)(8)

図 2：GaREUS で用いる数式

GaMD ではシステムの本来のポテンシャル U(x) にバイアス ΔUGaMD(U(x)) を加えることで、エネルギー障壁を下げる（図 2 の式(1)および図 3）。U(x) が閾値 E よりも小さい
ときのみ調和ポテンシャル型のバイアスを加え、それ以外は何も足さない（図 2 の式(2)）。このようにすることで、システムを局所安定な状態に捉えられることなく、状態空間を広くサンプリングすることができるようになる。U’(x) が連続関数かつ U(x) における大小関係（U(x1) < U(x2) ならば U’(x1) < U’(x2)）を維持するために、閾値 E はポテンシ
ャルの最大値 Umax と最小値 Umin によって条件づけられる（図 2 の式(3)）。また、正確なReweighting のために、U’(x) の分散は図 2 の式(4)を満たす必要がある。ここで、Uave と4

蛋白質科学会アーカイブ, 14, e105 (2021)

σU はそれぞれポテンシャルの U(x) の平均値と分散で、σ0 はユーザーが設定する分散の最大値である。k0 を k(Umax-Umin) と定義すると、閾値 E を下限 Umax に設定する場合、k0は図 2 の式(5)で決定する。ここでは最小限の説明しかしていないが、GaMD のパラメー
タの詳細は Miao の原著論文（文献 1,2）を参照していただきたい。
GaMD では二面角とそれ以外の全ポテンシャルエネルギーを別々にブーストする「デュアルブースト」が利用できる（文献 9）。二面角のブーストによってタンパク質の構造柔軟性を向上し、構造探索を加速する。同時に、二面角以外の全ポテンシャルエネルギーをブーストすることで、タンパク質近傍の水分子の運動等を加速し、ケージング効果等を回避する。デュアルブーストの場合は、上記
図 3：GaMD のイメージ図の GaMD パラメータが二面角とそれ以外ブーストポテンシャルによってエネルギー障壁が下がるの全ポテンシャルエネルギーに対してそれぞれ必要となる。
k と E を決めるには Umax, Umin, Uave, σU が必要となるが、系のポテンシャルエネルギー地形が複雑なため、一般的にはこれらは事前にわからない。まずブーストをかけずに MD計算を行い、ポテンシャルの形状を探索し、これらのパラメ ー タ を 決 め る 。 通 常 の MD の イ ン プ ッ ト に 加 え て 、[GAMD]セクションを追加する（図 4）。それぞれのパラメータの意味は、gamd は GaMD をするか否か、boost=noはブーストしない、boost̲type=DUAL はデュアルブース
トを用いる、thresh̲type=LOWER は式(3)の下限 Emax を用いる、sigma0̲dih と sigma0̲pot はそれぞれ二面角とそれ 以 外 の ポ テ ン シ ャ ル の σ0 と い う 意 味 で あ る 。
図 4：初期パラメータ設定のたupdate̲period は GaMD のパラメータ Umax, Umin, Uave, σUめの GaMD のインプットファを更新する頻度で、この段階では MD 計算のステップ数をイルの一部同じにしておく。
計算が終わると[OUTPUT]の gamdfile で指定したファイ
ルに図 5 のような結果が出力される。POT̲MAX, POT̲MIN, POT̲AVE, POT̲DEV が二面角以外のポテンシャルの Umax, Umin, Uave, σU であり、POT̲TH, POT̲KC, POT̲KC0 はそれらから計算した E と k となっている。同様に DIH が付いたものは二面角のパラメータである。

図 5：GaMD の初期パラメータの出力の例

5

蛋白質科学会アーカイブ, 14, e105 (2021)

初期パラメータが決まったら、実際にGaMD のブーストポテンシャルを加えながら
パ ラ メ ー タ を 更 新 し て い く 。 [GAMD] でboost=yes とし、pot̲max や dih̲max 等に先ほどの gamdfile に出力された初期パラメータを入力する（図 6）。ブーストポテンシャルを加えることで探索できるポテンシャル表面の領域が広がり、Umax, Umin, Uave, σU が変化する。update̲period=500 とすると 500 ステップごとに、pot̲max や dih̲max 等を新たなUmax, Umin, Uave, σU を基に更新する。パラメー
タの更新を繰り返すことで、探索できるポテンシャル表面が徐々に広くなっていき、より良い Umax, Umin, Uave, σU が見つかりやすくなる。
更新されたパラメータは update̲period ごとに gamdfile で指定されたファイルに出力される（図 7）。パラメータの時間変化を確認し、あまり変化しなくなったところで計算を終える。出力の最後の行の値で GaMD のパラメータを決定し、以降は変化させない。

図 6：パラメータを更新するための GaMD インプ
ットファイルの一部

図 7：GaMD のパラメータの更新の例

4）本計算
GaREUS では、GaMD のブーストポテンシャル ΔUGaMD(U(x)) を加えつつ、REUS シミュレーションを行う（図 2 の式(6)）。ΔUiREUS (ξ(x)) は REUS のバイアスポテンシャル、
ξ(x) は REUS で用いる反応座標、i はレプリカの番号である。ΔUiREUS (ξ(x)) は図 2 の式(7)のように調和ポテンシャル型になっており、レプリカごとに異なる反応座標 ξi に系を拘束するようになっている。ξi を定期的に交換することで、エネルギー障壁にとらわれず、反応座標空間を広く探索できる。

6

蛋白質科学会アーカイブ, 14, e105 (2021)

本計算用のインプットファイルでは、[GAMD]に 3）で決定した GaMD のパラメータを入力する（図 8）。
GaMD の パ ラ メ ー タ を 固 定 す る た め に 、update̲period は 0 にしておく。REUS も同時に行うので、[REMD]に REUS の設定を入力する。dimension
はレプリカ交換のパラメータの次元で、ここでは 1 とする 。 exchange̲period は 交 換 の 頻 度 で 、 こ こ で は1000 ステップに一度 REUS のパラメータを交換する。
REUS の 場 合 は type1 を RESTRAINT に す る 。
nreplica1 はレプリカの数である。rest̲function1 ではREUS で 用 い る バ イ ア ス ポ テ ン シ ャ ル と し て 、[RESTRAINTS]の拘束ポテンシャルの番号を指定する。
[OUTPUT]で GaREUS で出力されるファイルを指定する。remfile は各レプリカにおけるパラメータ番号が出力され、計算結果をソーティングする際に利用する。
REUS で用いるバイアスポテンシャルは図 9 のように指定する。まず、[SELECTION]でバイアスをかける原子を指定する。ここでは図 1 の点線赤丸が囲った酸素原子と水素原子をそれぞれ group1, group2 で指定している。[RESTRAINTS]ではバイアスの関数を指定する。 図 8：GaREUS のインプットファイルの一部
function が DIST とし、select̲index を 1 2 とすると、[SELECTION]で指定した group1 と group2 の原子間に距離の拘束ポテンシャルがかかる。調和ポテンシャルの係数と平衡値は、constant とreference で指定する。REUS で用いるレプリカ数に対応した数の値を指定する。図 9 のように reference の値をレプリカごとに変えることで、反応座標上のそれぞれの場所に拘束がかかることになる。

図 9：GaREUS で指定する反応座標パラメータ例

以上のインプットファイルを基に GENESIS で本計算を実行する。ログファイルに[REMD]の exchange̲period で指定したステップごとにレプリカの交換に関する情報が出力される（図 10）。この例では 4 と 5 が交換しており、AcceptanceRatio が 1 増えている。交換後のパラメータの番号とレプリカの番号（トラジェクトリーの番号）との対応は RepIDtoParmID と ParmIDtoRepID からわかる。[RESTRAINTS]で指定した反応座標のパラメータによっては交換率が低すぎたり、高すぎたりすることがあるが、その場

7蛋白質科学会アーカイブ, 14, e105 (2021)

合はパラメータを再調節する。全レプリカで AcceptanceRatio が均一に 20〜30％程度になるようにパラメータを選択するのが良いとされている。

図 10：GaREUS のログファイルの例
レプリカ交換の結果が出力されている

5）解析
GaREUS 法で得られたトラジェクトリーには GaMD のバイアス ΔUGaMD(U(x)) と REUSのバイアス ΔUiREUS (ξ(x)) がかかっている。トラジェクトリーから直接自由エネルギー地形を計算してしまうと、これらのバイアスがかかった不正確な地形が得られてしまう。正確な自由エネルギー地形を得るために、GaREUS のトラジェクトリーの確率分布をReweighting す る こ と で バ イ ア ス を 取 り 除 く 必 要 が あ る 。 こ こ で は 二 段 階 のReweighting を行う。まず REUS のバイアスを取り除き、続いて GaMD のバイアスを取り除く。
本計算で得られた各レプリカのトラジェクトリーでは、REUS のパラメータ交換によって 1 つのレプリカ内でバイアスのパラメータが時間によって変化している。つまり、トラジェクトリー内の各スナップショットには異なるバイアスがかかっている。REUS のバイアスの効果を取り除くにはバイアスパラメータごとにスナップショットがまとめられている方が扱いやすいので、まずトラジェクトリーをパラメータごとにソーティング。
GENESIS の remd̲convert に図 11 の 図 11：remd̲convert のインプット例
ようなインプットを入力し、ソーティン
グ す る 。 [INPUT] の dcdfile, remfile,logfile にそれぞれ GaREUS で得た各レプリカのトラジェクトリー（dcd ファイル）、

8

蛋白質科学会アーカイブ, 14, e105 (2021)

remfile、logfile を指定し、[OPTION]にステップ数や exchange̲period 等の対応する値を指定する。ソーティングされたファイルは[OUTPUT]で指定した trjfile や logfile に出力される。出力された dcd ファイル内では、1 つのバイアスパラメータのみが用いられている（例えば par1.dcd では 1 番目のパラメータのみ）。
ソーティングされたトラジェクトリーから反応座標の値を計算する。図 12 のようなインプットファイルを GENESIS の trj̲analysis に入力する。[TRAJECTORY]の trjfile に先ほどソーティングして得た dcd ファイルを指定し、md̲step, mdout̲period にGaREUS で用いた MD のステップ数と出力頻度を指定する。2 原子間の距離の場合は[OPTION]に distance を指定
すると、トラジェクトリーから距離を計算し、[OUTPUT]の disfile で指定したファイルに出力される。二面角を使う 場 合 は distance の 代 わ り にtorsion、角度の場合は angle とし、[OUTPUT]ではそれぞれ disfile の代わりに torfile, angfile でファイルを指定する。すべてのトラジェクトリー（つ 図 12：trj̲analysis のインプット例
まりすべてのレプリカの dcd ファイル）に対して同様の計算を行う。
MultiState Bennett Acceptance Ratio（MBAR）を用いて、REUS のバイアスを取り除く（文献 10）。図 13 のようなインプットファイルを GENESIS の mbar̲analysisに入力する。反応座標の結果を[INPUT]の cvfile に指定し、[RESTRAINTS]に REUS で用いた調和ポテンシャルバイアスのパラメータを指定する。[OUTPUT]で weightfile,fenefile, pmffile を指定すると、それぞ重み、レプリカ間の自由エネルギー変化、反応座

図 13：mbar̲analysis のインプット例

9蛋白質科学会アーカイブ, 14, e105 (2021)

標に沿った自由エネルギー地形（PMF）を出力できる。[MBAR]で MBAR の設定を指定する。temperature は計算で用いた温度にし、rest̲funciton で[RESTRAINTS]の関数を指定する。grids で PMF のグリッドを指定でき、0.0 15.0 301 とすると、反応座標が0.0 か ら 15.0 の 範 囲 で 等 間 隔 で 301 個 の グ リ ッ ド ご と に PMF を 出 力 す る 。
mbar̲analysis で特に重要なものは、weightfile で出力される重みで、これは GaREUSで得た各スナップショットごとの確率に対応している。この重みを用いることで、
GaREUS で指定していない反応座標空間にも射影することで、REUS のバイアスがかかてっていない自由エネルギー地形を求めることができる。ただ、mbar̲analysis で取り除けるバイアスは REUS のもののみで、GaMD のバイアスは取り除けていない。
最後に、GaMD のバイアスを取り除く。バイアスがかかっていない確率分布 p(ξ) は図2 の式(8)のように変形できる。ここでダッシュ記号がついているものは GaMD のバイアスがかかった分布や平均値である。MBAR によって REUS のバイアスは取り除かれているので、mbar̲analysis で得られた重みから p’(ξ) は計算できる。exp[βΔUGaMD] の平均は統計的に収束が非常に悪く、そのまま計算すると誤差が非常に大きくなる。そのため、この項はキュムラント展開をし、二次の項までで近似する（文献 11）。これは ΔUGaMD の分布を正規分布（Gaussian）で近似していることに対応する。ソーティングした logfileに出力されている GaMD のブーストポテンシャル ΔUGaMD と mbar̲analysis で得た重みから、p(ξ) を計算する。現在 GENESIS の Analysis Tools には GaMD バイアスのReweighting 用のプログラムは用意されていない。今後 GENESIS のウェブサイトで公開
予定の GaREUS のチュートリアルで Reweighting 用の Python スクリプトを配布する予定である。
図 14 は反応座標に沿った自由エネ
ルギー地形の結果である。GaMD のバイアスがかかったもの（つまりMBAR で REUS のバイアスだけ取り除いたもの）ではエネルギー障壁が低いが、全バイアスを取り除いたもの（Unbiased）ではエネルギー障壁が高くなっている。GaMD によって小さくなった障壁が Reweighting によって正しく戻されていることがわかる。
p(ξ) での反応座標 ξ は、GaREUS で用いた反応座標と同じである必要はな
く、任意の反応座標に対して計算する 図 14：GaREUS で用いた反応座標の自由エネルギー地形ことができ、複数の反応座標にも対応 GaMD のバイアスがかかったもの（GaMD biased）とバイア
し て い る 。 例 え ば 、 Alanine スを取り除いたもの（Unbiased）
Tripeptide の 2 つ の 二 面 角 φ（Ala1C-Ala2N-Ala2Cα-Ala2C）とψ（Ala2N-Ala2Cα-Ala2C-Ala3N）を用いて二次元の自由エネルギー地形を計算すると図 15 のようになる。GaMD のバイアスがかかっている状態（GaMD biased）だと全体的に浅い地形になっており、ベイスン間を移動しやすくなっている。バイアスを取り除いた地形（Unbiased）はエネルギー障壁が元通りに高くなっている。このように GaMD のバイアスをかけることで、REUS で用いた反応座10

蛋白質科学会アーカイブ, 14, e105 (2021)

標以外の自由度に対しても動きやすくなっており、広く構造探索ができるようになり、Reweighting によって正しい自由エネルギー地形に戻すことができる。

図 15：GaREUS で用いていない反応座標（二面角 φ と ψ）の自由エネルギー地形
GaMD のバイアスがかかったもの（GaMD biased）とバイアスを取り除いたもの（Unbiased）

以上が GaREUS による自由エネルギー計算のプロトコルであり、通常の REUS シミュレーションと比較して効率的に自由エネルギー地形を求めることができる。計算時間は系に含まれる原子数に大きく依存するため、一概に目安は述べることは難しい。対象とする系で通常の MD 計算を行い、計算時間を見積もってから GPU クラスターやスーパーコンピュータで計算することをお勧めする。また、GaREUS 法で用いた入力パラメータ等の詳細や一覧については、GENESIS のマニュアルにまとめられているので、そちらを参照してほしい。

11

蛋白質科学会アーカイブ, 14, e105 (2021)

文献[1] Miao, Y. et al., J. Chem. Theory Comput., 11, 3584‒3595 (2015)[2] Pang, Y. T. et al., J. Chem. Theory Comput., 13, 9‒19 (2017)[3] Sugita, Y. et al., J. Chem. Phys., 113, 6042­6051 (2000)[4] Fukunishi, H. et al., J. Chem. Phys., 116, 9058­9067 (2002)[5] Oshima, H. et al., J. Chem. Theory Comput., 15, 5199‒5208 (2019)[6] Jung, J. et al., Wiley Interdiscip. Rev. Comput.Mol. Sci., 5, 310­323 (2015)[7] Kobayashi, C. et al., J. Comput. Chem., 38, 2193­2206 (2017)[8] 小林千草, 蛋白質科学会アーカイブ, 11, e091 (2018)[9] Hamelberg, D. et al., J. Chem. Phys., 120, 11919‒11929 (2004)[10] Shirts, M. R. & Chodera, J. D., J. Chem. Phys., 129, 124105 (2008)[11] Miao, Y. et al., J. Chem. Theory Comput., 10, 2677‒2689 (2014)

このコンテンツはクリエイティブ・コモンズ 表示 - 非営利 - 改変禁止 4.0 国際 ライセンスの下に提供されています。
© 日本蛋白質科学会 2021 Licensed under クリエイティブ・コモンズ 表示 - 非営利 - 改変禁止 4.0 国際 ライセンス

12