生体分子の分子動力学シミュレーションにおける効率的な時間刻みの範囲
− グリシン、アラニン、バリン、ロイシン、イソロイシン3量体 −

佐藤 麗, 寺前 裕之, 石元 孝佳, 長嶋 雲兵


Return

1 はじめに

分子動力学法(MD: Molecular Dynamics)は、分子集合系の運動方程式を構成原子のすべてに対して逐次的に解いて、分子集団の軌跡を求める方法である[1]。しかし、運動方程式を数値的に解く際、時刻tに関する数値積分を行う必要がある[2]。この時間刻みの値は小さい方が定性的に正確な計算結果を得ることができる。
時間刻みの値と計算精度の関係は、Berendesen and Gunsteren が、調和振動子の初期値問題に対してEuler法, Verlet法, Gear法に対して、おおまかに計算誤差がそれぞれ、時間刻み幅の1乗、2乗、4乗で増加することを報告している[3]。
全シミュレーション時間は時間刻みとステップ数の積で表されるから、全シミュレーション時間が一定ならば時間刻みの値が小さいと計算に要する時間はより多くかかってしまうことになる。つまり、大規模生体分子などに対するMD計算を実行する際、時間刻みの値を大きくできれば計算時間の大幅な短縮につながると期待される。
生体分子のMDシミュレーションでは、フレキシブルモデルの場合、1.0 fs前後の時間刻みに設定し、共有結合の長さを平衡長にするように拘束をかけた場合、2.0 fs前後に設定するというのが一般的である。しかしながら、この時間刻み幅の妥当性に関する報告を調べたが、系統的な解析の報告が見あたらず、多くの生体分子のMDシミュレーションで時間刻み幅の設定が誤差の妥当性を考慮することなく、過去の経験を元に設定されているようであった。
そのため本研究では、アミノ酸3量体を計算対象としたMDシミュレーションでの効率的な時間刻みの範囲に関して研究を行い、時間刻みをどの程度まで大きく取ることができるのか、その許容範囲について調べた。

2 計算方法

本研究では、生体分子のMDシミュレーションソフトであるPEACH ver.5.8プログラム[4]を用いた。PEACHの積分法はVerlet法である。 5種類のアミノ酸(グリシン、アラニン、バリン、ロイシン、イソロイシン)の3量体を計算対象とした(Figure 1)。各アミノ酸3量体の構造はN末端側には水素イオン、C末端側には水酸化物イオンをそれぞれ付加し末端処理している。すべての計算において初期温度を300 Kとし、全シミュレーション時間は3.0 psとした。1つの計算を行うための時間幅である時間刻みを0.02 fsから3.0 fsまで変化させ、各々に対応したステップ数分の計算を行った。この時0.02 fsを基準時間刻みと定め、基準時間刻みの計算結果とその他の時間刻みの計算結果について比較を行った。


Figure 1. Chemical formula of aminoacid trimers

3 結果と考察

Figure 2 にグリシンのエネルギートラジェクトリを示した。時間刻みを変化させてもグリシンのエネルギートラジェクトリの概形はほとんど変化が見られない。しかし、時間刻みが増加する毎に徐々に高エネルギー側へシフトしていることがわかった。また、時間刻みが3.0 fsの時に、エネルギートラジェクトリが発散した。
その他のアミノ酸もグリシンと同様に、時間刻みを変化させてもエネルギートラジェクトリの概形はほとんど変化が見られず、ただ時間刻みが増加する毎に高エネルギー側へシフトし、時間刻みが3.0 fsの時に誤差によりエネルギートラジェクトリが発散した。
MD計算の結果として、基準時間刻みとその他の時間刻みについてのエネルギー差の平均値、最大値、最小値、標準偏差を求めた。(Figure 3 (a-d)) 平均値は時間刻みを大きくしていく毎に急速に増加している。このことから、時間刻みを大きくしていくと徐々に誤差が蓄積していくことを表している。このため時間刻み3.0 fsでは誤差が蓄積することにより発散したことがわかる。最大値、最小値、標準偏差では、時間刻み0.5 fsまでは急激な増加を示し、時間刻み0.5 fsからは緩やかな増加を示している。
特に標準偏差は0.5 fsから1.5 fsまではほとんど差がないことがわかる。また、平均値、最大値、標準偏差では誤差の大きさがともに各アミノ酸の分子量の嵩高い順、すなわちイソロイシン,ロイシン,バリン,アラニン,グリシンの順に並んでいる。ただし、最小値では分子量の嵩高さについての一定の順序は見られなかった。
Table 1Figure 3(a)の誤差の平均値に対する累乗近似関数の係数と冪、そしてR2を示した。これより、ほぼ時間刻み値の2.1乗で誤差が拡大していくことが判る。本研究での積分法はVerlet法であるので、Berendesen and Gunsteren [3]の結果と一致する。また係数がGly < Ala < Val < Leu < Ileの順で大きく、この順で誤差が大きくなる。この5種類はR2が0.99を超えているので、この累乗近似は非常によい近似となっている。


Figure 2. Energy trajectory of Glycine trimer

最後に、嵩高さについての検討としてアミノ酸の2量体を最適化し最適化構造とペプチド結合の二面角を0.5°回転させた時のエネルギー差を求めた。(Figure 4)
Figure 4は、ペプチド結合部分の二面角を回転させた時の回転障壁の程度を表しているのでこれは各アミノ酸の二面角の曲がり易さを表していることになる。このグラフでも、各アミノ酸の分子量の嵩高い順に並んでいる。つまり、MDシミュレーションの計算誤差の大小と各アミノ酸の分子量には相関関係があり分子量を大きくしていくと計算誤差も大きくなる。


Figure 3. Dt dependency of energy difference

Table 1. Power regression of averaged error (Figure 3(a)) as Y=CXa
PowerCoefficientR2
Gly2.12670.30450.9958
Ala2.11220.48740.9971
Val2.13930.49720.9970
Leu2.10110.66700.9972
Ile2.16520.75860.9974


Figure 4. Potential surface along a dihedral angle

4 まとめ

本研究では、アミノ酸3量体を計算対象としたMDシミュレーションでの効率的な時間刻みの範囲に関して研究を行い、時間刻みをどの程度まで大きく取ることができるのか、その許容範囲について調べた。基準時間刻みとその他の時間刻みについてのエネルギー差の標準偏差には、時間刻み0.5 fsから1.5 fsまではほとんど変化がない。すなわち、時間刻み0.5 fsから1.5 fsまでは同程度の精度を有していることを表している。つまり、この範囲内ではどの時間刻みでシミュレーションを行ってもあまり誤差がない。ゆえに、時間刻み1.5 fs付近でシミュレーションを行うのが時間短縮の面から効果的であると考えられる。更に、生体分子の分子量の大きいものでは誤差の度合いが大きくなることもわかった。

参考文献

[ 1] B. J. Alder, T. E. Wainwright, J. Chem. Phys., 27, 1208-1209 (1957).
[ 2] L. Verlet, Phys. Rev., 98, 159 (1967).
[ 3] Berendesen and Gunsteren, Practical algorithms for dynamic simulation, Proceedings of the International School of Physics << Enrico Fermi>> Course XCVII, Molecular-Dynamics Simulation of Statistical- Mechanical Systems, ed. by edi Ciccotti and Hoover, North-Holland (1986), pp.43-65.
[ 4] Y. Komeiji, Program for Energetic Analysis of bioCHemical molecules.
http://staff.aist.go.jp/y-komeiji/peach/peach.html


Return