分子動力学計算システムPEACHのWindows環境への移植とその応用−タンパク質の熱変性シミュレーション

中田 吉郎, 滝沢 俊治, 荻野 峰正, 宇野 雅美


Return

1 はじめに

現在、分子動力学計算法がタンパク質や核酸の研究分野で広く用いられている。そのプログラムは非常に複雑で大きなものであり主としてワークステーション用に開発されている。しかし近年、パソコンの性能が急速に発展してきたので、このような巨大なシミュレーションもパソコンで実行可能となってきている。
そこで我々は、古明地らによって開発されたPEACH(Program for Energetic Analysis of biochemical molecules)[1 - 3]をWindows環境に移植して広く利用できるようにし、そのシステムを用いて数種の系について実際にシミュレーションを行った。また、我々がこれまで開発してきた3次元分子表示プログラム[4, 5]をこのシステム用に変更して、PEACH自身には含まれていない表示部分を補って、より使いやすいシステムにした。
パソコンで実際に計算を実行してみると、有意なデータを得るためには計算時間がまだ非常に長くかかることが判ったので、パソコンを数台用いた並列計算[6]を利用して時間を短縮することも可能にした。
移植したシステムを用いて2つのタンパク質モデルの水中でのシミュレーションをおこない、シミュレーション温度による構造の時間変化を検討した。

2 分子動力学プログラムの移植

2. 1 Windows環境への移植

オリジナルPEACH[7]はUnix環境のFORTRAN90で開発されている。Windows環境への移植はVisual Fortranを用いて行った。移植に際しては若干のコンパイラーの違いによる文法上の修正を行ったが、ほぼオリジナルどおりに移植することができた。
計算に際してオリジナルではシェルスクリプトファイルを用いて実行をしているが、Windows版[8]ではバッチファイルで計算を実行させるように変更した。利用者が用意すべきファイルは、入出力ファイル名を指定するためのファイル(*.csh)と計算に必要なパラメータを指定するためのファイル(*.in)と必要なデータファイルである。
PEACHシステムに付けられているDNA(312D.pdb)のデモ計算について、Windows版での実行結果(時間)をTable 1に示す。これからパソコンでもワークステーションと同様の計算が可能であることがわかった。また市販の最高速のパソコンを用いればさらに計算時間を短縮できるはずである。

Table 1. Comparison between Windows version and Unix version.
WindowsUnix
Initial structure312D.pdb
Periodic box39×42×48 (A)
Atoms(solute ion solvent)
439, 6(Mg2+), 6981(H2O)
ProgramRUNMD(md4)
Simulation time2 ps
ComputerSONY Vaio LX80Alphastation 600au
CPUPentiumIV 1.5 GHzAlpha21164A 600 MHz
OSWindows2000Unix
CompilarVisual FortranDEC Fortran
Execution time1.49 h2.51 h

2. 2 計算結果のグラフィック表示

分子動力学計算で得られた分子構造などの表示には、我々が開発した3次元分子表示プログラムPDBDISP[3]、MDDISP[4]を組み込んだ。このプログラムでは計算前後の構造を重ね合わせて表示し比較することも可能である。シミュレーション結果のグラフ表示については、表計算プログラムEXCELに簡単に計算結果のデータを移せるので、そのグラフ表示機能を用いることができる。

2. 3 並列計算の導入

対象とする系が巨大になると計算時間が非常にかかるようになる。そこで本システムでは計算能力を上げるために、何台かの計算機をつなぎ協調作業させ、一台のときよりも速く計算させる並列計算という方法も利用できるようにした。その方法の一つとしてMPI(Message-Passing Interface)というライブラリーがある。これはMPIフォーラムより1994年に発表されMPICHとして Argonne National Laboratories より配布されている[9]。
並列計算までの手順は、まず始めにプログラムをライブラリーを用いて並列化する。PEACH Ver.3.3 の runmd モジュールはMPIライブラリーを用いて並列化されているので、それをWindows用に若干の変更を行った。次に複数台のWindowsNTマシン(NT,2000,XP)を用意し、それらをLANでつなぎ、それぞれにMPICHをインストールする。このときすべてのマシンのユーザ名とパスワードは同一にしておく。そのうちの1台をホストマシンとして並列化されたプログラムを実行すると並列計算が行われる。MPIはマルチアーキテクチャには対応していないので、同一のOS/CPUアーキテクチャの計算機をクラスタリングして一つの計算機として使うのに向いている。
市販のWindowsパソコン(PentiumIV、2 GHz)を使用して、1AG2を用いた場合の計算速度は10 ps当たり16時間であった(Table 2)。つぎにパソコンを3台用いて並列計算を行ったところ、10 ps当たり6時間ぐらいに高速化することができた。並列計算はパソコンの台数を増やすとその分計算速度は速くなるが、一方でネットワークを介してのデータのやり取りの時間もかかるようになるので、無制限に増やせばよいというものではない。今後はネットワークの速度や機器の性能の進歩によってこの部分の時間は短縮されると思われる。

Table 2. Execution time and Speed-up ratio in number of processors.
Number of ProcessorsExecution time (hours)Speed-up ratio
116.21.000
29.331.736
36.402.531
45.303.057

3 分子動力学計算による計算例

具体的なタンパク質の系について、このプログラムを用いてシミュレーションを行った。水中のタンパク質の性質として熱変性と呼ばれる現象がある。一般に水溶性タンパク質は70〜90℃で変性し失活することが知られている。そこでシミュレーション温度を変えてその立体構造を見ることによって、熱による変性が起きるかどうかを調べることにした。
今回の計算には、ホルモンタンパク質であるグルカゴン(1GCN.pdb)[10]と、ねずみの脳よりとられたプリオンタンパク質の中心部分(1AG2.pdb)[11]を用いた。
普通、球状タンパク質の熱変性は秒の単位で、また部分的なヘリックス−コイル転位はナノ秒からマイクロ秒のあたりで起きることが知られている。ここで用いたグルカゴンはアミノ酸29個からなる小さなタンパク質で、その立体構造はαへリックス構造であることが知られている。そこで数ナノ秒のシミュレーションで構造変化が起きる可能性が考えられるのでモデルとして用いた。またプリオンタンパク質は現在狂牛病の原因物質として注目されているタンパク質で、その立体構造の変化が狂牛病の原因ではないかと考えられている。熱に対する性質に他のタンパク質との相違があるかを検討するためのモデルとした。
Table 3に分子動力学シミュレーションに用いたモデル系を示す。

Table 3. Simulation models.
Amino AcidsPeriodic box(A)H2O molecules
1GCN2963×36×312197
1AG210353×53×534004

3. 1 水中のグルカゴン

グルカゴンのX線解析構造(29アミノ酸残基、1GCN)を用いて、水中での分子動力学計算によるシミュレーションを行った。設定温度を300 Kと350 Kと400 Kとして、熱的な性質を比較検討した。力場はAMBER94[12]をもちい、溶媒は水(2197分子)とし、周期境界条件、NTVの箱型モデル、クーロン力の計算にはエワルド法を用いた。Figure 1 に PDBDISP によるグルカゴンの立体構造を示す。


Figure 1. The initial structure of the glucagon model.

Figure 2に、溶質分子全原子の初期構造からの位置の変位(根平均二乗変位、RMSD)の時間変化を示す。このデータから見ると、300 Kにおいては、水中に入れた結晶構造が始めの100 psぐらいまでの間で緩んで大きく変位し、さらに熱運動により若干の変動が生じている。1000 ps以降はほぼ一定値を保っている。350 K, 400 Kに温度を高くしてシミュレートした場合には、400 psあたりからさらに別の構造変位が起こり、400 Kでは3000 psでまた別の構造変位を起こしている。350 Kでは2000 psのあたりでもとの構造(300 Kでの平衡構造)に戻っているように見える。


Figure 2. The RMSD as a function of time for the glucagon model.

次に1GCNの時間に対する二次構造の変化をFigures 3 - 5に示す。グラフの縦軸は1GCNのアミノ酸番号、塗りつぶし部分(水色)がαへリックス構造を表している。


Figure 3. Secondary structure in glucagon model as a function of simulation time at 300 K. The blue zone represents a-helix structure.

Figure 3より、300Kでは立体構造は安定化していてほぼ変化していないことが確認できた。またFigure 4から、350 Kでは加熱を始めてしばらくたってから少しαへリックスが解け、2500 ps後には初めとほぼ同じ場所にαへリックスができていることが分かる。Figure 5から、400 Kでは加熱を始めるとすぐにαへリックスが解け始め、5 ns後にはαヘリックス部分が非常に短くなっていることが分かった。これは熱変性によって立体構造が変化し、変性したためだと思われる。


Figure 4. Secondary structure in glucagon model as a function of simulation time at 350 K.


Figure 5. Secondary structure in glucagon model as a function of simulation time at 400 K.

Figure 6に400 Kにおいて5000 psシミュレートした時点での構造を WebLab ViewerLite[13]を用いて示す。これを見るとほぼヘリックス構造が解けてコイル構造になっていることが分かる。このようにグルカゴンのような短いヘリックス鎖の場合、温度を350 Kから400 Kに上げることにより数ナノ秒の時間でヘリックスコイル転移が起こることが示された。


Figure 6. The structure of glucagon after 5 ns simulation at 400 K.

3. 2 水中のプリオン分子モデル

NMR解析により得られたプリオンタンパク質の中心部分の構造(103アミノ酸残基、1AG2)を用いて、水中での分子動力学計算によるシミュレーションを行った。 Figure 7 にその構造を示す。設定温度を300 Kと400 Kとして熱的な性質を比較検討した。力場はAMBER94をもちい、溶媒は水(4004分子)とし、周期境界条件、NTVの箱型モデル、クーロン力の計算にはエワルド法を用いた。シミュレーション時間は数ナノ秒を目標とした。


Figure 7. The structure of the prion model.

Figure 8には、溶質分子全原子の初期構造からの位置の変位(根平均二乗変位)の時間変化を示す。この図から見ると、始めの200 psぐらいまでの間で水中に入れた結晶構造が緩んだ構造変位が生じている。その後300 Kの場合には徐々に構造が緩んできていることが読み取れる。一方400 Kに温度を高くしてシミュレートした場合は、300 Kの場合に比べて構造の変位が早く現れている。しかし5 nsまでのシミュレーション結果では、双方の温度での構造の違いは現れていない。これは今回用いた構造がプリオン分子の中心部分の構造のみで、残りのペプチド鎖の中心部分への影響が入っていないこと、シミュレーション時間が系の大きさに対して不足していることなどが原因として考えられる。また300 Kの値が平衡値になっていないことなどから力場パラメータの影響も考えられるので、他のパラメータなどを用いてその影響を検討することも必要と思われる。


Figure 8. The RMSD as a function of time for the prion model.

4 まとめ

生体分子用分子動力学システムPEACHをWindows環境で動作するように移植することができた。PEACHシステムに付けられているDNAのデモ計算について、Windows版で実行した結果、UNIX版と同じ値を得ることができた。また計算時間に関しては、有用な結果を得るためには並列計算を用いることで対処できることがわかった。
さらに2つのモデルタンパク質について、このシステムを用いて分子動力学シミュレーションを行った。グルカゴンモデルは300 K, 350 K, 400 Kでシミュレーションを行い、プリオンモデルは300 K, 400 Kでシミュレーションを行った。またグルカゴンモデルの400 Kは350 Kで1000 psまで計算し、その後400 Kに温度を上げた。これらのグラフを比べると両方とも結晶構造からのずれ(RMSD)は、300 Kよりも400 Kの方が大きくなっていることが分かった。このことからタンパク質は300 Kでのシミュレーションでは安定であるが、350 K以上では2次構造を中心とした構造変化が起きていると考えられる。またグラフより、グルカゴンモデルはプリオンモデルよりも結晶構造からのずれが大きいことが分かった。つまりグルカゴンモデルは熱変化により構造が変化していると考えられる。このことからグルカゴンモデルはプリオンモデルよりも分子量が小さいため相互作用を受けやすいと考えられる。

参考文献

[ 1] Komeiji, Y., Uebayashi, M., Takata, R., Shimizu, A., Itsukashi, K., Taiji, M., J. Comp. Chem., 18, 1546-1563 (1997).
[ 2] 古明地勇人、上林正巳、長島雲兵, 生体分子の分子動力学シミュレーション(1)方法, J. Chem. Software, 6, 1-36 (2000).
[ 3] 古明地勇人、田島澄恵、原口誠、高橋伸幸、上林正巳、長島雲兵, 生体分子の分子動力学シミュレーション(2)応用, J. Chem. Software, 7, 1-28 (2001).
[ 4] 中田吉郎、滝沢俊治、上林正巳、中野達也, OpenGLを用いた生体分子の立体構造表示プログラム, J. Chem. Software, 6, 157-164 (2000).
[ 5] 中田吉郎、滝沢俊治, Windows版PEACHとそのシミュレーション表示, 化学ソフトウエア学会年会2001研究討論会, 埼玉 (2001).
[ 6] Komeiji, Y., Haraguchi, M., Nagashima, U., Parallel Computing, 27, 977-987 (2000).
[ 7] http://staff.aist.go.jp/y-komeiji/peach/peach.html
[ 8] http://www3.nibh.go.jp/~nakata/peach/peachw.html
[ 9] http://www-unix.mcs.anl.gov/~ashton/mpich.nt/
[10] Sasaki, K., Dockerill, S., Adamiak, D. A., Tickle, I. J., Blundell, T., Nature, 257, 751 (1975).
[11] Riek, R., Hornemann, S., Wider, G., Billeter, M., Glockshuber, R., Wuthrich, K., Nature, 382, 180 (1996).
[12] Cornell, W. D., Cieplak, P., Layly, C. I., Gould, I. R., Merz, K. M., Ferguson, D. M., Spellmeyer, D. C., Fox, T., Caldwell, J. W., Kollman, P. A., J. Am. Chem. Soc., 117, 5175-5197 (1995).
[13] WebLab ViewerLite Version 4.0, Molecular Simulation Inc. (2000).


Return