2014 年 11 月 27 日
温度レプリカ交換分子動力学法についての Q&A
ここでは Sugita, Okamoto1)が開発した拡張アンサンブル法の一つである温度レプリカ交換分子動力学法（T-REMD）を実行する場合に考慮しなければならない点を記載します。
レプリカ交換法は今や、様々な MD プログラムパッケージ（GROMACS, NAMD, AMBER,CHARMM, GENESIS など）に組み込まれてきています。更に、他の MD プログラムを利用してレプリカ交換を専門に行うインターフェースプログラムも自由に使えるようになっています(MMTSB2), REIN3)など)。T-REMD の詳細については原論文 1)や総説 4,5)、解説 6)を参照して下さい。また、論文 7)にはレプリカ交換法の批評があります。
T-REMD は複数のシステムのコピー(レプリカ)を用意して独立に MD 計算を実行し、ある周期で隣り合う温度を持つレプリカ間で温度パラメータをメトロポリスクライテリアに従って交換する事で、それぞれのレプリカの温度がランダムに変化し、それがエネルギーのランダムウォークを誘起します。そのエネルギー空間上のランダムウォークがシステム（例：タンパク質、リガンド、水からなる系）の大きな構造変化を誘起し、構造空間上のサンプリング効率を上げるというものです。実行するにあたっては、どの様なシステムで、何を見るのかを決めた上で、サンプリング効率を上げる為に必要な温度範囲とレプリカ数を決めます。
１． どのような温度幅（最低温度と最高温度）を設定すべきか？
一般的には、最高温度は「エネルギー極小状態に留まらないよう、十分高くとる必要がある」ことが基準となります。
通常、Weighted Histogram Analysis Method (WHAM)を用いて、議論する温度での自由エネルギー解析を行います。簡易的に行う場合には論文などで議論する温度のトラジェクトリを集めて解析します。いずれにしても、議論する温度領域のサンプリングが十分取れる様に、最低温度と最高温度を選ぶ必要があります。最低温度は、最終的に解析したい系の温度に設定してもよいでしょう。最高温度は何を調べたいかによって変わります。シス
テムのフォールディングを調べたい場合にはできるだけ高い温度にする必要があります。
局所的な揺らぎ、構造変化、構造配置などを調べるときには、他の部分が壊れない様にある程度の温度にするか、ある程度高温にして束縛を入れます。
T-REMD を行った結果、最低温度から最高温度までが適切に交換されていることを確認する必要があります。それをトンネリングと言います。例えば、最低温度の系と最高温度の系とが T-REMD シミュレーションの間で全く交換されていない、すなわちトンネリングしていない場合は、温度とレプリカの数の再設定が必要になります。
参考までに、
Sugita らの論文 1)では、真空中での Met-enkephalin の構造サンプリングを扱っています1

が、その時の温度幅は 200K から 700K です。水分子を含んだシステムの場合には、「典型的な水溶液中のペプチドの例では、最低温度 300K ほど、最高温度 450K ほどです」との記述もあります 6)。
２．どれほどのレプリカを用意すべきか？
それぞれのレプリカは他のレプリカと交換されるために用意されているのですから、交換受容確率（exchange acceptance probability）が小さすぎると、T-REMD の意味がありません。最適な交換受容確率の理論的根拠はありますが、レプリカの交換方式によって異なります。ランダムに交換させる場合には 20%程度が適切であるとされています。ただし、実際にはシステムによる影響が大きく、計算毎にトライアンドエラーで決めます。参考の、0.149（239⇔286K）、為に、Sugita らの論文 1)での交換受容確率は、0.160（200⇔239K）
0.143（286⇔342K）、0.139（342⇔409K）、0.142（409⇔489K）、0.146（489⇔585K）、0.146（585⇔700K）となっています。この数値は T-REMD の実行後に得られますが、短いテスト実行から外挿します。全ての交換受容確率がだいたい同じ値であることが望ましいです。
最低温度と最高温度が設定されると、交換受容比率に応じて必要なレプリカ数はある程度推定できます。指数的に温度を決める簡単な方法もありますが、Patriksson ら 8)はモデ
ル粒子を用いたレプリカ数とその温度を設定するための方法のホームページを開いています 9)。いずれにしても、大ざっぱな目安を与えているだけなので、実際にはシステムによって調整する必要があります。
３．各レプリカに対して初期配座をどのように設定すべきか？
システムの大きさに依りますが、通常、普通に MD 計算を実施するときと同じように、最初にそれぞれの温度で minimization し、平衡状態シミュレーションをナノ秒程度行った後に実行します。レプリカ交換法を実行させる時、「適切なパラメータ（温度幅とレプリカ数など）を用い、REMD シミュレーションによってサンプリングが適切に行われているのなら、異なった初期配座から始めても、全て同じ初期配座から始めても、計算結果は同じはずである」と言われています。
Met-enkepharin の構造サンプリングを試みていますが、以下は例ですが、Sugita ら 1)は、その全てのレプリカの出発構造は同じで、Met-enkepharin の伸びきった構造です。ただ、
T-REMD 計算によるデータ収集を行う前に、各温度での通常の分子動力学計算を 100ps 行い、さらに平衡化のために 100ps の T-REMD 計算を行っています。
Nguyen ら 10)はアミロイド・ベータの断片 Aβ16-22 の２量体（や３量体）の構造サンプリングを行うにあたり、以下の方法で初期配座を決めています。「(1)結晶構造既知の Aβ10-35から Aβ16-22 を切り出し、ランダムに配置して２量体の構造を作る。(2)２量体を八面体の箱の中に入れ、3,089 個の水分子を加える。(3)最急降下法で系のエネルギーを極小化する。(4)2定圧（1 atm）下で 1ns 平衡化計算を行う。(5) 定温（300K）下で 1ns 平衡化計算を行う。
これらは Berendsen の方法で行った。(6)系をさらに定温（300K）、定積下で 1ns 平衡化計算する。
」この結果得られた配座が全てのレプリカの初期構造として T-REMD 計算に用いられています。Nguyen らが設定した温度幅は最低温度 290K、最高温度 400K でした。また、T-REMD 計算を 55ns 行っていますが、最初の 5ns のデータは解析から除外しています。
Park ら 11)は survivin２量体の構造サンプリングを行っています。用いた２量体構造は
PDB:1E31 で 142 残基からなる A 鎖と B 鎖の断片 survivin5-117 です。そして以下の方法で初期配座を決めています。
「(1)まずその断片２量体を切頂八面体の箱に入れ、水分子 TIP3Pを加え、エネルギー極小化を行う。 (2)温度幅を 300K から 340K とし、4K ごとに 10 個のレプリカを設定する。なお、T-REMD 後のターゲット温度は 308.9K とされている。(3)エネルギー極小化された構造を 10 個のレプリカの初期構造とし、50ps の温度平衡化分子動力学計算を行う。(4)ついで、200ps の圧力平衡化計算を行う。これらは Berendsen の方法で行った。(5)この結果得られた 10 個のレプリカごとに、２量体のみの配座と速度を取り出す。
」ここで取り出された 10 個の配座と速度が、水を陽には取り扱わない T-REMD 法の初期構造となっています。
Miyashita らは、Implicit solvent/Membrane model を用いて APP の膜貫通部位の膜中での構造予測 12) と、二量体構造予測 13) を行いました。構造予測 12) は、温度範囲を300-700K/32 replica にしており、この高温はシステムがアンフォールディングする温度にしています。初期構造はそれぞれの温度で、minimization と平衡化計算を 200ps 行っています。およそ 10ns の計算を行い初期構造依存性を避ける為に最後の 5ns で解析をしています。
二量体構造予測 13)は、温度範囲を 300-500K/32 replica としており、最高温度は二量体の配置を見るので、APP 自体のヘリックスが壊れすぎない温度にしています。32 個の異なるインターフェースを持つ初期構造を作り、minimization と平衡計算をそれぞれ 200ps 行いました。そして、32 レプリカそれぞれにその構造を１つずつ充てて、レプリカ交換シミュレーションを実施しています。およそ 20ns 程度の計算を行い、最後の 10ns を解析に用いています。
４．レプリカ交換の試みは何ステップごとに行うべきか？
理論的には１ステップに 1 回の交換でも全く問題ありません。しかし、実際の計算ではうまくいかない場合が多いです。結局、理論とは矛盾するのですが、緩和過程を考える必要がでてくるからです。また、レプリカ交換のプログラムにも依ります。杉田らのオリジナルのレプリカ交換分子動力学プログラムなど、温度スケーリングを用いているプログラ
ムでは短い交換周期でも問題ありません(REIN3)など)。一方、交換の際、熱浴を利用して温度変化させるプログラムの場合(MMTSB2)など)は、緩和過程も考慮する必要があります。
3

目安として、システムの大きさにもよりますが多くの場合、大体 2ps(1000 step)程度毎の交換の計算がなされています。経験上、大きなシステムでは緩和に時間がかかるので長めにした方が良いです。この周期でないといけないという決まりはありませんが、正しい統計を得るのに十分なステップ数が必要になるという事です。
Nguyen ら 10)は Berendsen の方法で温度制御をしていますが、その時の緩和時間τは0.1ps です。そして、分子動力学計算の時間刻みは 2fs、レプリカ交換を試みる時間間隔は1.5ps としています。つまり 750 ステップ毎にレプリカ交換を試みています。
Ostermeir ら 5)は、Periole ら 14)の研究を参照しつつ、レプリカ交換の後に十分な緩和の時間を得るためには、レプリカ交換は少なくとも~1ps の間隔を置いて行うべきだろうと述べています。
参考文献：
1) Sugita Y and Okamoto Y (1999) Replica-exchange molecular dynamics method forprotein folding, Chem Phys Lett, 314, 141-151
2) Feig M et al. (2004) MMTSB Tool Set: enhanced sampling and multiscale modelingmethods for applications in structural biology, J Mol Graph Model, 22, 377-95.
3) Miyashita N et al. (2014) REIN: Replica-exchange INterface for simulating protein
dynamics and function, Int J Quantum Chem, DOI: 10.1002/qua.24785
4) Mitsutake A et al. (2013) Enhanced sampling algorithms, Methods Mol Biol, 924,153-95
5) Ostermeir K and Zacharias M (2013) Advanced replica-exchange sampling to study
the flexibility and plasticity of peptides and proteins, Biochim Biophys Acta, 1834,847-536) http://www.gromacs.org/Documentation/Tutorials/にある 2013 年 GROMACS 利用者のための Workshop and Conference 資料や、
http://biomol.bme.utexas.edu/wiki/index.php/ABS:REMDなど（2014 年 11 月 27 日アクセス確認）
7) Zuckerman DM and Lyman E (2006) A Second Look at Canonical Sampling ofBiomolecules using Replica Exchange Simulation, J Chem Theory Comput,

2(4):12001202
8) Patriksson A and van der Spoel D (2008) A temperature predictor for paralleltempering simulations, Phys Chem Chem Phys, 10, 2073-2077
9) http://folding.bmc.uu.se/remd/（2014 年 11 月 27 日アクセス確認）
10) Nguyen PH et al. (2011) Effects of all-atom force fields on amyloid oligomerization:
replica exchange molecular dynamics simulations of the Aβ16-22 dimer and trimer, Phys
Chem Chem Phys, 13, 9778-9788
411) Park IH and Li C (2010) Dynamic ligand-induced-fit simulation via enhanced
conformational samplings and ensemble dockings: a survivin example, J Phys Chem B,114, 5144-5153
12) Miyashita N et al. (2009) Structures of β-amyloid peptide 1-40, 1-42, and 1-55-the672-726 fragment of APP-in a membrane environment with implications forinteractions with γ-secretase, J Am Chem Soc, 131, 17843-17852
13) Miyashita N et al. (2009) Transmembrane structures of amyloid precursor protein
dimer predicted by replica-exchange molecular dynamics simulations, J Am Chem Soc,131, 3438–3439
14) Periole X, Mark AE (2007) Convergence and sampling efficiency in replica exchange
simulations of peptide folding in explicit solvent, J Chem Phys, 126, 014903
以上
追記：
本文を作成するに当たり理化学研究所生命システム研究センターの宮下尚之博士のご協力を得た。ただし、文責は江口にある。

5