Journal of Computer Chemistry, Japan, Vol. 1, No. 1 (2002)

生体高分子のモンテカルロシミュレーションにおける分子変形方法の工夫

目黒俊幸, 山登一郎

1 はじめに

現在、ヒトを始めとして様々な生物のゲノムの解読が精力的に進められており[1]、その結果として膨大な量のタンパク質の一次配列情報が容易に入手できるようになってきている。このような状況において、タンパク質の一次配列から立体構造を構築する原理を明らかにし、一次構造から三次構造を予測することを実現できたならば、ゲノム情報を３次元の空間内で起きる生命現象とダイレクトに関連付けながら利用することが可能になる[2, 3]。このような理由からタンパク質の構築原理の解明と構造予測法を確立することは現代生命科学における重要なテーマとなっている。
タンパク質の構造予測は、そのエネルギー的側面に注目すると、配列の指定されたペプチド鎖の自由エネルギー最小の構造を求めることと言い換えることができる。しかしペプチド鎖の取りうる構造の数は莫大であり、全ての構造について自由エネルギーを求め、その中から最小のものを選び出すような操作は現在の計算機の能力では不可能である[4]。更に、自由エネルギー曲面には部分的な極小値が多数存在しているため、最適化をスムーズに行うことができない。これらの事情によって、この最小値探索問題を一般的に解くことは未解決で困難な問題となっている。この解決を目指して、様々な構造予測法が試みられている。それらを大別すると、既知の構造情報を利用して精度の高い予測を目指す経験的な方法と、原子間相互作用といった物理化学的な原理に基づき演繹的に折り畳み過程を再現しようとする非経験的な方法とに分けられる。
経験的な方法としてはホモロジーモデリング[5, 6]や3D-1D法[7]などがあげられる。これらの方法は、予測すべき構造に類似の構造がすでに実験的に解かれている場合に有効であるが、そういった構造情報がない場合は使用できず、また折り畳み機構の原理の解明を直接目指すものではないという欠点がある。
非経験的な方法としては、シンプルなモデルを用いて明快かつ厳密な理論的解析を目指す格子モデル法[8]や、本来はタンパク質の動力学的性質を調べるためのモレキュラーダイナミクス法（MD法）[9, 10]や統計力学的性質を調べるためのモンテカルロ法（MC法）[11, 12]を折り畳みシミュレーションに応用したもの等があげられる。
格子モデル法はタンパク質の構造をより単純な格子モデルに置き換え、様々なエネルギー計算やMC法と組み合わせてシミュレーションを行う方法であって、理論的扱いが他の方法と比較して容易である。しかし、現実のタンパク質の構造形成には側鎖同士のタイトパッキング[13] などのような構造形成に関与する多数の要因があると考えられるが、格子モデルではそれらを直接扱うことができないという弱点がある。
単純化されたモデルではない実際の分子構造に忠実なモデルを用いるシミュレーション方法としては、全原子モデルを用いたMD法やMC法が行われている。実際のタンパク質の構造エネルギー曲面は、天然構造のエネルギーが他の構造より顕著に低い折り畳みに都合のよいファネルモデル[14]と呼ばれる形状になっていると考えられている。そのようなエネルギー局面を計算機上で再現することができれば、シラミつぶし的な構造探索を行わなくてもMD法やMC法を使うことでタンパク質の天然構造にたどり着けると期待される。
MD法による構造予測では、ニュートン力学に基づいたシミュレーションにより、タンパク質を含む系がとりうるエネルギーの低い構造を求める。この方法は現実の物理化学的性質をよく反映したものであるが、その反面、長大な計算時間が必要であるという問題を抱えている。実際のタンパク質の折り畳みはミリ秒から秒、分のオーダーの現象であるのにもかかわらず、現代の計算機能力ではサブマイクロ秒相当のMDシミュレーションがようやく行えるというのが現状である[15]。
MC法による構造予測では[16, 17]計算上で分子を人為的に変形させそのエネルギー変化に応じて構造を変えていく。原子の直交座標を変数とした場合は、MD法のほうがMC法よりもサンプリング効率が高い。だが、分子変形の自由度を減らして計算時間を短縮したい場合、たとえば二面角を変数にする場合などは、MC法のほうが扱いやすい。よって、今回の方法では、MC法を採用した。MC法は原子レベルで表現されたペプチド鎖モデルのシミュレーションを比較的短い計算時間で実行できる、優れた方法といえる。このMC法のアルゴリズムや使用する力場関数の改良を続けていけば、ペプチド鎖の構造予測法を確立することができると期待される。本論文ではその一環としてMC法の分子変形方法の改良を試みた。
一般的に、MC法の分子変形操作ではタンパク質の二面角をランダムに変更させる二面角法を採用する場合が多い。この二面角法は、ペプチド鎖が短い場合は、扱いやすく効率的に安定構造を探索しうる有効な方法であるが、長いペプチド鎖の天然構造を再現するための最良の分子変形方法とは限らない。
二面角法以外のより適切な分子変形方法を用いてMCシミュレーションを行えば、実現しやすい分子構造もそれに応じて変化し、シミュレーション全体の天然構造再現効率が向上する可能性もある。我々は平行移動法と名づけた変形方法を新たに導入し、二面角法と平行移動法を用いて実際にMCシミュレーションを行い、その特徴を比較検討した。本論文ではその結果について報告する。

2 材料と方法

2. 1 MCシミュレーションに用いたポリペプチド鎖の由来とシミュレーションの条件

MCシミュレーションの対象として2種類のポリペプチドを用いた。
一つはRibonucleaseA のN末から13残基分にあたるCぺプチドであり、実験的にαヘリックス形成能があることが報告されている[18]。Cペプチドの配列は
LYS GLU THR ALA ALA ALA LYS PHE GLU ARG GLN HIS MET　である。
もう一つは28残基のポリペプチド（以後28merと呼ぶ）で、NMRによる構造解析で２本のストランドと一本のヘリックスをもつββα構造をしていることが報告されている [19] 。この配列は、LYS PRO TYR THR ALA ARG ILE LYS GLY ARG THR PHE SER ASN GLU LYS GLU LEU ARG ASP PHE LEU GLU THR PHE THR GLY ARG である。

2. 2 MCシミュレーションに用いたソフトウェアとハードウエア

計算に用いたMC法プログラムORTは当研究室でC言語を用いて開発した。ORTは二面角法と平行移動法の両方法を自由に組み合わせて分子の変形を行うことができる。このプログラムについての問い合わせは東京理科大学基礎工学部山登研究室で受け付けている。
計算はLinux(RedHat6.1)をインストールしたDOS/Vパソコン（cpuはpentiumⅢ800MHz相当）上で実行した。

2. 3 MCシミュレーションに用いた分子モデル

計算に用いる分子モデルは全原子モデルを用い、N末はNH₃⁺、C末はCOO^-とした。分子モデルを変形する際は、共有結合の長さは固定しておき、結合角度および二面角を変更するようにした。
全てのシミュレーションは直鎖状の初期構造から開始した。

2. 4 ポテンシャルエネルギーの計算方法

分子のポテンシャルエネルギー計算用の力場パラメータはAMBERのparm94.datファイルのものを用いた[20]。ただし、二面角の力場パラメータによりαヘリックス構造やβストランド構造への偏向が生じること[21]を回避するため、主鎖の二面角に関するエネルギー項は計算から除外した。
各エネルギー項目を計算する際のカットオフ半径は8.0 Åとした。
また水素結合に関しては　N－H…O　のなす角度が90゜以下では引力が弱くなり（0.0倍）、180゜に近くなると引力が強くなる（2.5倍）ように補正した。
今回の計算では水分子は明示的には取り扱っていない。そこで、溶媒としての水を連続体近似し、誘電率を工夫した。具体的には、微視的に見てペプチド鎖が水に囲まれているような系では均一な誘電率を定義できないことを考慮し、原子間距離に応じて変化する誘電率関数を用いた。一般的にタンパク質-水系のシミュレーションでは、ε＝ｒ　（ε は誘電率、ｒは原子間距離、この関数は原子間距離が3～5Åの範囲で、誘電率が3～5とタンパク質内部の誘電率と近い値をとる）という誘電率関数が多用される。しかし、本論文では誘電率を距離の2乗とする ε＝ｒ²　を採用した。この関数は比較的計算が容易であり、Mehlerらが水溶液中のイオンに関する実験値から求めた誘電率関数[22]に近い値をとる。30残基以下の小さいペプチド鎖のアンフォールディング状態では、電荷を持った原子の多くは水に露出していると考えられるため、誘電率関数は ε＝ｒとするよりもε＝ｒ²　として原子間距離が 3～8Å程度の範囲で水溶液中のイオンに関する誘電率を近似したほうが現実的であると期待される。
疎水性相互作用をエネルギー計算式に取り入れるために、疎水性側鎖間のポテンシャルエネルギーを見積もる際には、通常のvan der Waals相互作用や静電相互作用の他に、疎水性相互作用の効果に対応するエネルギー項を新たに定義し、計算に加えた。以下に、2つの側鎖間の疎水性相互作用エネルギーE_np（kcal/mol）を、各側鎖のCβ原子同士間の距離ｒ（Å）から求める式を示す。

ｒ ≦ 5.0 E_np= k_w×k_np×｛100.0*(r-5.0)²-1.0｝

5.0 ＜ｒ ≦ 6.5 E_np= k_w×k_np×（-1.0）

6.5 ＜ｒ E_np= k_w×k_np×MIN（5.0*(r-6.5)²-1.0 , 0.0）

（MIN(a,b)は、aとbの値を比較して小さい方の値を返す関数）

ここで、k_wは疎水性相互作用エネルギーの計算上の重み係数であり、今回はk_w＝18.0とした。また、k_npは各側鎖の疎水性の程度に応じて決められる重み係数である。今回の計算では計算対象の2つの側鎖が、強疎水性と強疎水性の組み合わせになっているときはk_np＝1.0、強疎水性と弱疎水性の組み合わせになっているときはk_np＝0.7、弱疎水性と弱疎水性の組み合わせになっているときはk_np＝0.5とした。このようにk_npを決めるには、20種類のアミノ酸側鎖を強疎水性か弱疎水性か非疎水性かに、あらかじめ分類しておく必要があるが、今回はEngelmanらの疎水性尺度[23]を参考にして、PHE、MET、ILE、LEU、VALを強疎水性グループ、CYS、TRP、ALA、ARG、LYSを弱疎水性グループ、残りのアミノ酸を非疎水性グループに分類した（ARGとLYSは側差の付け根部分がアリファティックな炭素鎖となっているため、ここでは弱疎水性グループに分類した）。
温度に対応するパラメータTは多くの場合は2000 Kで定温とした。この値は室温である300Kに較べて大きくなっているが、これは本論文では水を陽に含まない系を用いているので、水素結合や静電相互作用の形成によるエネルギーの利得が大きく見積もられてしまうことに対してバランスを取った結果である。

2. 5 平行移動法と二面角法の変形操作

本論文では分子の変形操作として平行移動法と二面角法という２つの異なった変形操作法を必要に応じて採用した。
平行移動法は、ペプチド鎖を適当な共有結合の前後で二つの部分に分け、その一方を任意の方向に平行移動させることによって分子を変形させる（Figure 1, left）。移動方向は全方位の中から乱数を用いて偏りなく決定する。移動距離も乱数を用いて決定する。一回の変形操作あたりの動かしうる最大距離をシミュレーション試行ごとに定め、これを最大移動距離と呼ぶことにすると、実際の移動距離は０より大きく最大移動距離以下の範囲でランダムに決定される。
一方、二面角法では適当な共有結合部分を軸として二面角を変化させることで変形操作を行う（Figure 1, right）。一回の変形操作で動かしうる最大の回転角度を以後、最大回転角度と呼ぶことにすると、実際の回転角度は０より大きく最大回転角度以下の範囲で決定される。
二面角法では回転軸から離れるにつれ、変形操作による原子の移動距離が大きくなり、原子同士の衝突が起き易くなるが、平行移動法では変形部位からの距離が近くても遠くても各位置での原子の移動距離は等しく、原子同士の衝突が起きにくいと考えられる。なお、二面角法では構造パラメータのうち二面角だけが変化するが、平行移動法ではそれに加えて結合角も変化する。この結合角に起因するポテンシャルエネルギーがあるため平行移動法を採用した場合の全エネルギーは二面角法を採用した場合よりやや高くなる傾向がある。

Figure 1. Schematic illustration of conformational deformation procedures. Left: Parallel translation method. Atoms B and C move the same distance in the same direction. Right: Dihedral angle rotation method. The turn axis is d. Dislocation of C is larger than B.

2. 6 MCシミュレーションの実行手順の概要

2. 6. 1 主鎖部分の構造変形方法

プログラムORTによるMCシミュレーションでは、主鎖部分の構造変形をするために下記の手順１から３までの操作が繰り返し行われる。なお手順１では変形方法が二面角法と平行移動法で異なるため、両方の変形方法を並記しておく。
（手順１）分子の任意の変形:
各ペプチド鎖の主鎖上にはペプチド結合による平面部が（残基数－１）個存在するが、その中から一平面を選択し変形操作を行う。ペプチド平面の主鎖部には、Figure 2で示すように、Ｎ末から順にCA1－C、C－N、N－CA2という３本の共有結合がある。
平行移動法では、その中のCA1－CとC－Nの２本の結合に対して独立かつ同時に変形操作を行う(Figure 2, left)。この変形操作は、具体的にはCA1からCに向かう方向ベクトルおよびCからNに向かう方向ベクトルの向きをそれぞれ定められた範囲内でランダムに変化させることによって行う。なお、NからCA2に向かう方向ベクトルは独立に動かさずに、常にCA1からCに向かう方向ベクトルと一致させる。
一方、二面角法ではCA1－C結合を軸とした回転操作とN－CA2結合を軸とした回転操作を同時に実行する(Figure 2, right)。２つの回転操作の回転角度はそれぞれ独立に乱数で決める。
(手順２）変形前後の構造エネルギー差の算出:
上記の<ポテンシャルエネルギーの計算方法>に記したやり方で、変形前後の２つの構造について、それぞれのポテンシャルエネルギーを求める。
(手順３)エネルギー差に基づいて分子の変形を採用するか否かの判定:
メトロポリス法[24] に基づき変形後の構造を採用するか棄却するかを判定する。

Figure 2. Illustration of a peptide bond subjected to deformation. Left, parallel translation method: We applied our deformation procedure to the directions of vectors of three covalent bonds, CA1-C (shown in black), C-N (grey), and N-CA2 (black), simultaneously at one step. Varying only a vector CA1-C, we can change the relative positions of atoms CA1 and CA2. Therefore, the direction of vector N-CA2 is maintained to be the same as that of vector CA1-C. However, we move the direction of vector C-N independently of vector CA1-C. The purpose of changing the vector C-N is to rotate the peptide plane. Right, dihedral rotation method: We applied our deformation procedure to the two covalent bonds (shown in black) on both sides of a peptide bond (shown in grey) simultaneously at one step. The deformation degrees of the dihedral angles y and f are independently decided. Therefore, the value of y is different from f.

この（手順１）から（手順３）までの一巡の作業をシミュレーション進行過程の１つのまとまりと考え以後では１stepと呼ぶことにする。
今回の全てのシミュレーションでは（手順１）における主鎖上のペプチド平面部の選択は、stepの進行に伴ってＮ末から順々に選んでいきC末に達したら再びN末に戻って繰り返すという順序で行った。ここで、ペプチド平面の選択個所がN末からC末まで一巡する間に行われる一連の操作をまとめて１sweepと呼ぶことにする。例えば28merのシミュレーションでは１sweepには27step分の操作が含まれることになる。

2. 6. 2 側鎖部分の構造変形方法

今回行ったシミュレーションでは主鎖部を２sweeps連続で変形しそれから側鎖部をN末からC末まで１残基づつ順番に変形するという組み合わせでシミュレーションを行った。
ORTでは基本的に側鎖部の変形操作は平行移動法の場合も二面角法の場合も同じ方法で行う。具体的には、ある残基を変形する場合は一度に側鎖部の全共有結合部に関して平行移動法（最大移動距離0.5Å）でランダムに変形し、メトロポリスの方法でその変形を採用するかどうかを判定する。ただし側鎖に存在するメチル基や芳香環のように現実には回転運動を起こしていると見なせる部分については平行移動法の代わりに二面角法（最大回転角度5゜）で変形操作を行う。

2. 6. 3 全シミュレーション過程を通しての設定事項

１試行のシミュレーションあたり100000 sweeps（主鎖部の変形操作だけカウントした場合）の計算を行った。主鎖部を2 sweeps計算する毎に側鎖部をN末からC末まで１巡分変形操作し、さらにペアリストの更新を行った。
本論文では様々な条件下でのシミュレーションを行っているが、ある条件に関する計算結果について議論する場合は、できるだけ偶然に左右されないように、各条件で乱数系列だけが異なるシミュレーションを10試行繰り返して行い、それらの結果を総合して解釈するようにした。

3 結果および考察

3. 1 CペプチドのMCシミュレーション

二面角法と平行移動法を用いてCペプチドのMCシミュレーションを行い、その結果を比較検討した。構造形成の有無を判断する指標としてヘリックス構造を形成している残基数を用いた。各残基がヘリックス構造になっているかどうかの判定はDSSP法[25] で行った。MC法では分子変形の規模の大きさが構造探索の効率に大きな影響を与えるため、どちらの分子変形方法を用いる場合でも`の程度を様々に変えてシミュレーションを行った。
Figure 3は横軸に変形の大きさ、縦軸にMCを100000 sweeps行った時点でのヘリックス形成残基数をとってグラフにしたもので、(a)が平行移動法、(b)が二面角法を用いた場合の結果である。それぞれのグラフには４通りの温度パラメータ（1000 K、2000 K、3000 K、4000 K）の下でのシミュレーションの結果が重ね書きされている。平行移動法を用いた場合はヘリックスが効率よく形成されると思われる適切な条件は、温度パラメータが1000 Kから2000 K程度、変形の大きさを表す最大移動距離は、温度パラメータが1000 Kの時は0.3Åから1.0Å、2000 Kの時は0.3Åから0.7Å程度であった。温度パラメータが3000 K以上ではほとんどヘリックスは形成されなかった。一方、二面角法を用いた場合は、ヘリックスを形成するのに適切な条件は、温度パラメータが2000 K、最大回転角度20゜から120゜であった。1000 Kや3000 Kでもヘリックスは形成されるが、明らかに2000 Kの時よりヘリックス形成残基数は減少していた。4000 Kではほとんどヘリックスは形成されなかった。

Figure 3. Dependence of the number of residues that form alpha helix on the magnitude of deformation (at 100000 sweep). a: Parallel translation method. b: Dihedral rotation method.

両方法での最もヘリックス形成効率がよいと思われる条件同士を比較すると、平行移動法（温度2000 K－最大移動距離0.25Å）では6残基程度しかヘリックスを形成していなかったのに対して、二面角法（温度パラメータ2000 K－最大回転角度35゜）では全13残基中10残基以上がヘリックスを形成していた。また二面角法（温度パラメータ2000 K）のヘリックス形成残基数は最大回転角度が 15゜から120゜までのほとんどの範囲で8残基以上となっており、平行移動法の最適条件での値を上回っていた。このことから、αヘリックスの形成に関しては二面角法の方が平行移動法より有利に働くことがわかる。これは、一本鎖ヘリックスのような単純な構造を対象にした場合には原子移動に伴う衝突が問題にならないため、二面角法の一回の変形操作による原子移動が大きいという利点がそのままシミュレーション全体での探索効率の良さに反映されるためと考えられる。
なお平行移動法では100000 sweep目のヘリックス形成残基数が少ないが、これは平行移動法を用いた場合にヘリックス構造が発達しないということを意味している訳ではない。今回は１条件あたり乱数系列のみを変えた10試行分のMCシミュレーションを行っているが、各試行において最初から100000 sweep時点までの全過程の間でヘリックス形成残基数が最大でどこまで大きくなるかを調べてみると、平行移動法を用いた場合でも、10回の試行のうち、どの試行でもヘリックス形成残基数が8残基以上になっている時期が存在しており、さらに、そのうちの5回の試行ではヘリックス形成残基数が10残基以上になる時期があった。すなわち、平行移動法でも瞬間的には長いヘリックスが形成されるが、そのヘリックスがすぐに崩れてしまうため、適当なsweep時点での10試行分の平均をとった場合には低いヘリックス形成残基数しか得られないのである。この、ヘリックスが別の形に変わり易いという問題は、エネルギー計算に用いる力場関数とも密接に関係しており、将来より正しい力場パラメータが開発されていくにつれて改善されていくことが期待される

3. 2 各変形方法を用いた場合のCペプチドのヘリックス形成速度

次いで、各変形方法を用いた場合のヘリックス形成速度を知るために、“平行移動法－2000 K－最大移動距離0.3Å”という条件と“二面角法－2000 K－最大回転角度30゜”という条件でのMCシミュレーションについて、開始時から終了時までのヘリックス形成残基数の経時変化をグラフ化した（Figure 4）。

Figure 4. Evolution of the number of residues that form alpha helix.Black circle: Dihedral rotation method(2000 K. The maximum rotation angle is 30°).Black square: Parallel translation method(2000 K. The maximum movement distance is 0.3A).

二面角法では500 sweeps、平行移動法では5000 sweeps程度のMCシミュレーションを行えば、グラフはほぼ一定の値になり以後は目立った増減は見られなかった。以上のことから、ヘリックス形成速度に関しても、二面角法は平行移動法より効率的で、10倍以上速くヘリックス構造を形成することができることがわかった。
平衡状態に到達するまでの実際の計算時間は、両方法とも200 sweepsあたり12秒ほど計算時間を要することから推定すると、二面角法で30秒、平行移動法でも300秒ほどであり、極めて短い計算時間でシミュレーションは平衡状態に到達するといえる。

3. 3 28merポリペプチドのMCシミュレーション

次いでやや複雑な天然構造をもつ28merポリペプチドを対象にして、二面角法と平行移動法を用いてMCシミュレーションを行い、その結果を比較検討した。このポリペプチドはββα構造をとることがNMRによる実験から知られている。
一般的に、天然構造に対してより類似度の高い構造を頻繁に発生させることができる分子変形方法が、構造予測を目指す上で効率的と考えられる。28merポリペプチドのシミュレーションでは、立体構造の類似度を表す指標としてRMSD(root mean square displacement,根平均二乗変位)を用いた。すなわち、天然構造との類似度が高い構造は、天然構造に対するRMSD値が低いものとし、各シミュレーション手法の有効性を、それから得られる構造のRMSD値がどの程度低くなるかに注目して判断した。同時に、RasMol[26]等のグラフィックソフトで得られた構造を表示し、数値に表れにくい立体的特徴にも注意を払った。
最初に平行移動法と二面角法を用いた場合に、100000 sweep時点でRMSDがどの程度まで低下しているかを調べた。Figure 5(a)は平行移動法かつ2000 Kの定温条件下で最大移動距離を様々に変えてMCを行った結果である。同様に、Figure 5(b)は二面角法かつ2000 Kの定温条件下で最大回転角度を様々に変えてMCを行った結果である。平行移動法を用いた場合は、最大移動距離が0.4ÅのときにRMSD値が5.6Åとなり最低値になった。また、最大移動距離が0.2Åから0.9Åの範囲でRMSD値は6.5Å付近の値をとっていた。一方、二面角法を用いた場合は、最大回転角度が40゜の場合にRMSD値が6.7Åとなり最低値になった。また、最大回転角度が10゜から110゜の範囲でRMSD値は7Å前後の値をとっていた。このことから、28merのMCシミュレーションでは、平行移動法の方が二面角法よりもRMSD値の低い構造を効率よく発生させうると考えられる。

Figure 5. Dependence of RMSD (at 100000 sweep) between calculated and NMR structures of 28-mer on the magnitude of deformation.Figure (a) shows RMSD calculated using the parallel translation method and (b) shows RMSD calculated using the dihedral rotation method. RMSDs of 10 MC trials were averaged and shown for each simulation condition.

なおここでのRMSD値は、それぞれの条件ごとに10回の試行を平均した値であり、また現状ではモンテカルロ法に用いられる力場関数の精度が十分なものではないため、5Å以上という大きな値となっているが、個別のシミュレーションでは3Å台のRMSD値も現れている。Table 1は100000 sweep時点の10試行分のRMSD値を平均しないで小さい順に並べたもので、上段は平行移動法で最大移動距離0.4Å、下段は二面角法で最大回転角度が40゜での結果である。このTable 1中で最も低いRMSD値は、平行移動法では3.4Å、二面角法では4.6Åであった。

Table 1. Distribution of the RMSD values (A) between calculated and NMR structures of 28-mer at 100000 sweep.

maximum minimum average

Parallel translation 3.4 4.3 4.3 4.5 5.0 6.5 6.8 7.0 7.0 7.5 5.6

Dihedral rotation 4.6 4.8 5.9 6.4 6.8 7.1 7.5 7.8 7.9 8.4 6.7

	maximum		minimum	average
Parallel translation	3.4	4.3	4.3	4.5	5.0	6.5	6.8	7.0	7.0	7.5	5.6
Dihedral rotation	4.6	4.8	5.9	6.4	6.8	7.1	7.5	7.8	7.9	8.4	6.7

ところで、ここでの議論は、「天然構造に対するRMSD値が低い構造は、天然構造と似ている」という前提のもとになされているが、厳密にはRMSDの数字だけで類似度を判断することはできない。たとえば、ある2つの構造のRMSD値が同じであったとしても、一方はより天然構造に近い構造に容易に変形できるが、他方は立体障害等によってそれができないような場合は、両者の天然構造に対する類似度は異なると考えることもできる。そのため、シミュレーションで計算された構造が実際にどのような形態的特徴を持つかを確認することは重要である。
例としてFigure 6に、Table 1に示したRMSD値に対応する20個の構造の中から6個の構造を選び表示した。上段の（A）～（C）は平行移動法によって得られた構造、下段の（D）～（F）は二面角法によって得られた構造で、右にいくほど天然構造（NMR構造）に対するRMSD値が低くなる。これら全ての構造はFigure 6の左端の直鎖状初期構造から計算された。C末側のαへリックスは多くの構造で再現されるが、N末側のβストランド部は再現される場合とされない場合とがあり、その構造上の差異がRMSD値に大きな影響を与えていると考えられた。このような事情はFigure 6に示した以外の構造に関しても同様であった。このような状況では、モデル分子の構造変化に対する特別な立体障害等の要因が存在しないため、「RMSD値が低いにもかかわらず天然構造に近づきにくく、その意味で天然構造との類似度が低い」というような構造は発生しにくいと考えられる。以上のことから、28merポリペプチドのシミュレーションでは天然構造に対する類似度は、RMSDで大まかに近似できると考えられる。

Figure 6. The calculated 28-mer polypeptide conformations. Only CA atoms in the backbone are shown. (A)-(C) are obtained using the parallel translation method. (D)-(F) are obtained using the dihedral rotation method. The RMSD value is shown at the left side of each structure. The figures were created with RasMol [26].

3. 4 各変形方法を用いた場合の28merポリペプチドのヘリックス形成速度

次いで、構造が形成される速度を調べるため、“平行移動法－最大移動距離0.4Å”の条件と、“二面角法－最大回転角度40゜”の条件でのRMSDの平均値と最小値の経時変化を調べた。Figure 7はRMSDの平均値の経時変化である。グラフから、平行移動法は初期には二面角法よりもRMSD値が高いが、最終的には二面角法よりRMSD値が低くなることがわかる。平行移動法では変形操作によって移動する全ての原子が同じ距離だけ少しずつ移動するので、シミュレーションが進んで構造が複雑になってくると二面角法と比較して原子の衝突が少なく、効率よく分子変形が行えるものと考えられる。

Figure 7. Evolution of the RMSDs between calculated and NMR structures of 28-mer. The black line shows RMSDs obtained using the parallel translation method (2000 K. The maximum movement distance is 0.4A) and the grey line shows RMSDs obtained using the dihedral rotation method (2000 K. The maximum rotation angle is 40°). RMSDs of 10 MC trials were averaged and shown.

Figure 8はそれぞれの条件での10試行分のRMSDの最小値の経時変化を示している。平行移動法での最小RMSD値は、シミュレーションの初期には二面角法よりも高いが、10000 sweep以降はほとんどの場合で二面角法より低くなり、20000 sweep以降からはしばしば 2Å台の値をとるようになった（Figure 9）。

Figure 8. Evolution of the RMSDs between calculated and NMR structures of 28-mer. In this figure, we show minimum RMSD values among ten MC trials at each sweep, giving an estimation of the formation speed of better structure without the influence of failing MC trials. The black line refers to the parallel translation method (2000 K. The maximum movement distance is 0.4A); the grey line relates to the dihedral rotation method (2000 K. The maximum rotation angle is 40°).

Figure 9. An example of low RMSD (2.85A) conformation of 28-mer polypeptide obtained by MC simulation using the parallel translation method (grey sticks) with the corresponding NMR structure (white sticks) superimposed. Only CA atoms in the backbone are shown. The figure was created with RasMol [26].

なお、28merの平行移動法では100 sweepsあたりの計算時間は30秒ほどであるので、RMSD値が2Å台の構造が初めて出現するまでに要する計算時間は約100分である。

3. 5 両変形方法を組み合わせたMCシミュレーション

これまでは、各変形方法を単独で用いていたが、今度は両変形方法を組み合わせた時の効果を調べるために平行移動法と二面角法を1 sweepごとに交互に行うシミュレーションを実行した。最大移動距離と最大回転角度が様々に異なる45通りの条件を決め、１条件あたり乱数系列だけが異なる10回のMCシミュレーションを行った。温度パラメータは全ての試行において2000 Kで定温にした。その結果得られた各条件での100000 sweep時のRMSDの平均値をFigure 10に示す。最大移動距離0.65Å－最大回転角度50゜の組み合わせで10試行分の平均RMSD値が最も低下し5.8Åとなった。この値は平行移動法（2000 K、最大移動距離0.4Å）を単独で行った場合のRMSD値である5.6Åよりはやや大きい値であった。次いで、Figure 10で低いRMSD値を達成する最大移動距離と最大回転角度の組み合わせ条件を３通り選び、シミュレーション開始時からのRMSDの平均値と最小値の経時変化をグラフ化し構造形成速度を調べた（Figure 11）。　3通りの組み合わせ条件のグラフはそれぞれ、“平行移動法のみ－最大移動距離0.4Å”という条件のグラフと、“二面角法のみ－最大回転角度40゜”という条件のグラフとの中間の形状をしていた（Figures 7, 8参照）。この結果から判断する限りは、両方法を混合して用いることが、平行移動法を単独で用いることより特別に有利であるとは思えない。ただし、今後、平行移動法と二面角法の混合方式を工夫していけば、両者の長所を備えたより良い分子変形方法が開発できる可能性は残っている。

Figure 10. The RMSDs between 28-mer structures calculated using the combined deformation method (at 100000 sweep, 2000K) and NMR structure. The combined method means simulation using parallel translation and dihedral rotation methods alternately each sweep. The X-axis shows the maximum movement distance for the parallel translation method. The Y-axis shows the maximum rotation angle for the dihedral rotation method. RMSDs of 10 MC trials at 100000 sweep were averaged and shown. The thick curves in the figure indicate contour lines (for RMSD= 6.0, 7.0, 8.0 and 9.0). The low RMSD region in this figure shows the optimal combinations of deformation parameters.

Figure 11. Evolution of the averaged RMSD values between calculated 28-mer structures using combined deformation methods (Figure 10) and NMR structure. We selected three arbitrary pairs of parameters at low RMSD region in Figure 10. Maximum movement distance and maximum rotation angle are 0.65A and 50°(black line), 0.80A and 50°(dark grey line), and 0.50A and 20° (light grey line), respectively. In figure (a), average RMSD values of 10 MC calculations are shown. In figure (b), minimum RMSD values among 10 MC calculations are shown.

以上のことより、28merのMCシミュレーションにおいて天然構造付近の構造を探索するためには、平行移動法を単独で用いることが最も効率が良いと考えられる。この結果はCペプチドの構造形成には二面角法の方が効率がよいという結果と対照的である。これは目指す構造が複雑になるにつれて、変形操作による原子の衝突が起きにくいという平行移動法の特徴が有利に働いてくるものと思われる。
Cペプチドや28mer以外のペプチド鎖の折り畳みシミュレーションではどちらの分子変形方法が有利になるかは、実際にやってみなければ不明であるが、今回の例で見たように、平行移動法は条件によっては二面角法よりも有利な方法になると期待される。

4 結論

MCミュレーションに用いる分子変形方法をより効率の良いものに改良することを目指し、平行移動法と二面角法とを比較検討した。
実際にCペプチドのMCシミュレーションを行ったところ、二面角法の方が平行移動法よりも迅速に天然構造に近いと思われるヘリックス構造を形成した。これは、構造の比較的単純なCペプチドでは二面角法の変形操作にともなう原子移動距離が大きいという利点が強く現れたためと考えられる。
この反対に、28merのMCシミュレーションでは、平行移動法を用いた方が二面角法を用いた場合よりも天然構造に近いRMSD値の低い構造が効率よく出現した。これは、平行移動法を採用した場合は変形操作によって動かされる各原子の移動距離が皆等しく小さいことから、複雑なペプチド鎖に二面角法を適用した場合に問題となる原子衝突による新規構造採用率の低下が回避されたためと考えられる。この確認のため、実際にあるフォールディング中間体のペプチド鎖中央部の共有結合（13番のSERのψ、14番のASNのφ、ψ、15番のGluのφ）に関して、平行移動法と二面角法とを適用した場合の新規構造の採択率（400steps中の採択回数の％）を比較した。平行移動法では最大移動距離が0.3Åの場合は採択率14.3％、0.4Åの場合は8.5％、0.5Åの場合は5.5％であったのに対して、二面角法では最大回転角度が30№.5％、40№.0％、50№.0％であった。このことは上記の考えを支持するものである。
平行移動法を用いたMCシミュレーションでは、protein Aに由来する33mer[27]のヘリックス-ループ-ヘリックス構造等も再現できることを確認しており、現在、各種構造に対する有効性を検討している。
今回の結果から判断すると、Cペプチドのような単純なαヘリックス構造を形成するものに対しては二面角法、28merのようなやや複雑な構造を持つものに対しては平行移動法が適していると考えられる。なお28merに関しては相乗効果を期待して平行移動法と二面角法を組み合わせたMCシミュレーションも行ったが、平行移動法を単独で行った場合以上の効果は得られなかった。
将来、平行移動法と正確な力場関数とを組み合わせて用いることが可能になれば、多様な分子変形経路が効率よく提供され、その中から天然構造に向かう適切な経路を選び出すことのできる理想的なMCシミュレーションが実行できるようになると期待される。

参考文献

[ 1] International Human Genome Sequencing Consortium, Initial sequencing and analysis of the human genome, Nature, 409, 860-921 (2001).
[ 2] Fetrow, J. S. and Skolnick, J, Method for prediction of protein function from sequence using the sequence-to-structure-to-function paradigm with application to Glutaredoxins/Thioredoxins and T1 Ribonucleases, J. Mol. Biol., 281, 949-968 (1998).
[ 3] Simons, K. T., Strauss, C. and Baker, D., Prospects for ab initio protein structural genomics, J. Mol. Biol., 307, 429-445 (2001).
[ 4] Levinthal, C., Are there pathways for protein folding?, J. Chim. Phys., 65, 44-45 (1968).
[ 5] Blundell, T. L., Sibanda, B. L., Sternberg, M. J. E. and Thornton, J. M., Knowledge-based prediction of protein structures and the design of novel molecules, Nature, 326, 347-352 (1987).
[ 6] Greer, J., Comparative modeling of homologous proteins, Methods Enzymol., 202, 239-252 (1991).
[ 7] Bowie, J. U., Luthy, R. and Eisenberg, D. A, method to identify protein sequences that fold into a known three-dimensional structure, Science, 253, 164-170 (1991).
[ 8] Hinds, D. A. and Levitt, M. A, lattice model for protein structure prediction at low resolution, Proc. Natl. Acad. Sci. USA, 89, 2536-2540 (1992).
[ 9] Momany, F. A., McGuire, R. F., Burgess, A. W. and Scheraga, H. A, Energy parameters in polypeptides. VII. geometric parameters, partial charges, non-bonded interactions, hydrogen bond interactions and intrinsic torsional potentials for naturally ocurring amino acids, J. Phys. Chem., 79, 2361-2381 (1975).
[10] McCammon, J. A., Gelin, B. R. and Karplus, M, Dynamics of folded proteins, Nature, 267, 585-590 (1977).
[11] Northrup, S. H. and McCammon, J. A, Simulation methods for protein structure fluctuations, Biopolymers, 19, 1001-1016 (1980).
[12] Noguti, T. and Go, N., Efficient Monte Carlo method for simulation of fluctuating conformations of native proteins, Biopolymers, 24, 527-546 (1985).
[13] Shakhnovich, E. I. and Finkelstein, A. V., Theory of cooperative transitions in protein molecules. I. Why denaturation of globular protein is a first-order phase transition, Biopolymers, 28, 1667-1680 (1989).
[14] Onuchic, J. N., Nicholas, D. S., Luthey-Schulten, Z. and Wolynes, P. G., Protein folding funnels: the nature of the transition state ensemble, Folding & Design, 1, 441-450 (1996).
[15] Duan, Y., Wang, L. and Kollman, P. A., The early stage of folding of villin headpiece subdomain observed in a 200-nanosecond fully solvated molecular dynamics simulation, Proc. Natl. Acad. Sci. USA, 95, 9897-9902 (1998).
[16] Hansmann, U. H. E. and Okamoto, Y., New Monte Carlo algorithms for protein folding, Current Opinion in Structural Biology, 9, 177-183 (1999).
[17] Shimada, J., Kussell, E. L. and Shakhnovich, E. I., The folding thermodynamics and kinetics of crambin using an all-atom Monte Carlo simulation, J. Mol. Biol., 308, 79-95 (2001).
[18] Brown, J. E. and Klee, W. A., Helix-coil transition of the isolated amino terminals of ribonuclease, Biochemistry, 10, 470-476 (1971).
[19] Dahiyat, B. I., Sarisky, C. A. and Mayo, S. L., De Novo protein design: Towards fully automated sequence selection, J. M. Biol., 273, 789-796 (1997).
[20] Cornell, W. D., Cieplak, P., Bayly, C. I., Gould, I. R., Merz, K. M. Jr., Ferguson, D. M., Spellmeyer, D. C., Fox, T., Caldwell, J. W. and Kollman, P. A., A second generation force field for the simulation of proteins, nucleic acids and organic molecules, J. Am. Chem. Soc., 117, 5179-5197 (1995).
[21] Ono, S., Nakajima, N., Higo, J. and Nakamura, H., Peptide free energy profile is strongly dependent on the force field: Comparison of C96 and AMBER95., J. Comput. Chem., 9, 748-762 (2000).
[22] Mehler, E. L. and Eichele, G., Electrostatic effects in water-accessible regions of proteins, Biochemistry, 23, 3887-3891 (1984).
[23] Engelman, D. M., Steitz, T. A. and Goldman, A., Identifying nonpolar transbilayer helices in amino acid sequences of membrane proteins, Annu. Rev. Biophys. Biophys. Chem., 15, 321-353 (1986).
[24] Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. and Teller, E., Equation of state calculations by fast computer machines, J. Chem. Phys., 21, 1087-1092 (1953).
[25] Kabsch, W. and Sander, C., Dictionary of protein secondary structure: Pattern recognition of hydrogen-bonded and geometrical features, Biopolymers, 22, 2577-2637 (1983).
[26] Sayle, R.A. and Milner-White, E. J, RasMol: Biomolecular graphics for all, Trends Biochem. Sci., 20, 373-375 (1995).
[27] Braisted, A.C. and Wells, J. A., Minimizing a binding domain from protein A, Proc. Natl. Acad. Sci. USA., 93, 5688-5692 (1996).

Return

ｒ ≦ 5.0	E_np= k_w×k_np×｛100.0*(r-5.0)²-1.0｝
5.0 ＜ｒ ≦ 6.5	E_np= k_w×k_np×（-1.0）
6.5 ＜ｒ	E_np= k_w×k_np×MIN（5.0*(r-6.5)²-1.0 , 0.0）