生体分子の分子動力学シミュレーション(1)方法

古明地 勇人, 上林 正巳, 長嶋 雲兵

目次

1. 　古典分子動力学シミュレーションの概要
   
2. 　初期構造と初期速度
   
3. 　力とポテンシャルの計算
   
4. 　時間積分
   
5. 　アンサンブル（温度制御など）
   
6. 　トラジェクトリーの解析方法
   
7. 　おわりに
   

1 古典分子動力学シミュレーションの概要

　本総説では、生体分子分子動力学用ソフトウェアPEACH（Program for Energetic Analysis of bioCHemical molecules,古明地, 1995）に導入されたアルゴリズムを中心に、タンパク質やDNAなど生体高分子の分子動力学シミュレーション（MD）の方法を解説する。
　古典分子動力学法（古典MD）とは、一般に「N体多体系の状態の時間発展を、ニュートンの運動方程式の数値解を差分法により求めてシミュレーションすること」と定義できる。つまり、連立微分方程式

　を、時間tに関する差分法に直して解くことである。ここで、m_i、r_i、は粒子iの質量と位置、U は系全体のポテンシャルエネルギーである。この場合、粒子は生体高分子の構成原子である。
　最近は、量子力学と組み合わせたMDも盛んであるが、本総説では省略し、古典力学に基づいたMDについてのみ解説する。古典MDの教科書を巻末に何冊か挙げておいた[

Figure 1. 分子動力学シミュレーション(Molecular Dynamics, MD)

2 初期構造と初期速度

2. 1 初期構造

2. 2 境界条件と溶媒の扱い方

Figure 2.2. 生体分子のMDにおける溶媒の扱い

　(2)はShell waterと呼ばれ、生体分子の周りに何層かの水を配置したもの。溶媒の分子数が少ない割には、比較的よい結果が得られる（Guenot & Kollman [13]など）。ただし、水分子に束縛を掛けないと、エントロピー効果のため水分子が蒸発してしまうことがある。
　(3)はCap waterと呼ばれ、ある原子（あるいは溶質の重心など）を中心にして球状に水を発生させる。球の中の水は自由だが、外に出て行こうとすると、軽い束縛力（バネの力）により中に押し戻されるようにする。
　従来はFigure 2.2(3)Aのように、タンパク質の一部分だけを覆って活性部位のダイナミクスを追ったり自由エネルギー計算を行ったりする場合に使うことが多かった。水が野球帽のような形になることからCap waterと呼ばれる。だが、PEACHでは、Figure 2.2(3)Bのように生体分子をすっぽりと覆った、言わば、ゴム風船の中に水と分子を閉じ込めた形にして使うことが多く、CapよりもむしろDropletと呼ぶほうがふさわしい。このDropletでも(2)のShellでも、水の表面が真空に露出しているという問題点がある。それが、水や溶質の挙動にどの程度影響を及ぼすかは、1999年の時点では系統的な報告は見あたらない。
　周期境界条件（Periodic boundary condition、Figure 2.2(4)）とは、溶質と溶媒を基本セル（実像）の中に入れ、そのレプリカ（虚像）からの寄与も計算に考慮する方法である。この方法では、(2)や(3)と違って表面が真空に露出しない。溶液や固体のMDでは、周期境界条件を使うのが普通であり、また、計算機の進歩や高速アルゴリズムの普及により（Darden et al. [16]など）生体分子のMDにも多く使われるようになった。
　周期的でない生体分子を周期的に扱うことによる問題点については、まだ、詳細な調査はない。しかし、Smith et al. [17, 18]、De Souza & Ornstein [19]などは、問題なさそうだとしている。溶質部分に注目するか、溶媒部分に注目するかでも、結論は異なるだろう。
　なお、Figure 2.2のどの溶媒モデルがよいかは一概に言えない。対象とする分子の性質と大きさ、シミュレーションしたい条件、求めたい物理量、利用できる計算機資源等により、選択は変わるからである。現在、生体高分子のMDでよく使われるのは、droplet((3)B)と周期境界(4)である。

2. 2 .1 周期境界条件での最近接粒子の選択

　周期境界条件の場合(Figure 2.2(4))、Figure 2.2.1のように数多い虚像の中から、最近接粒子のみを選択して相互作用を計算することが多い。これは、粒子iとjの相互作用を計算するときに、粒子jの実像と26個の虚像（図は２次元だが、実際は３次元）の計27個の中から、最もiに近いものを選んで計算する方法である。Allen & Tildesley [6]の教科書にいくつか方法が載っているが、PEACHでは、以下のアルゴリズムを利用している。

　xboxは、箱の一辺の長さである。また、anintはfortranの組み込み関数で、anint(a)はaを四捨五入した整数である。y、zについても同様に計算する。

2. 2. 2 系の電荷の中和

1. 　中和できる数の電荷を溶質全体に加える
   　たとえば、タンパク質全体が＋２の電荷を持っている場合、合計−２の電荷を全原子に均等にばらまけば、系は中和する。筆者の経験では、タンパク質分子量１万につき、電荷が１〜２までなら、この方法が有効である。しかし、電荷がこれ以上大きくなると、しばしば構造が崩れる場合があるので、あまり勧められない。
   
2. 　アミノ酸の荷電状態を変えて調節する
   　一部のアミノ酸の電荷を調節して、電荷を０にする。不自然な電荷分布になるからと嫌がる人もいるが、分子全体のダイナミクスに限れば、普通は大きな影響はないようである。系の一部の挙動の解析のためにMDを行うならば、無関係な部分に余計な電荷を押し付けるのも一つの方法であろう。
   
3. 　イオンを発生させる
   　系の電荷が大きい場合は、中和するだけの数のイオンを発生させる。一般的には、これが最も現実的で汎用性がある方法である。PEACHでは、発生させたイオンと溶質の間のクーロン相互作用が最小になるような位置の溶媒分子を選んで、それとイオンと置換することで発生させている。結果的に電荷をもった残基の近くに発生されることが多い。なお、電荷を中和させるためのイオンを発生させる方法には、それらを全くランダムに配置する、あるいは電荷を持った残基の近くに配置する、などの方法もある。
   

2. 3 エネルギー極小化（Energy Minimization）

2. 3. 1 最急降下法（Steepest descent method）

        DR = DRini
        do i = 1, nstep        
                compute F and ENEnew    !　力とエネルギーを計算
                DENE = ENEnew - ENEold
                dRMS = sqrt(SF2)/sqrt(Ndfree)       !   Ndfreeは系の自由度
                EXIT if either DENE or dRMS is converged
                if DENE > 0 then
                        DR = 1/2 DR
                else
                        DR = min(DRmax, 1.2DR)
                end if
                R = R + DR /fRMS *F
                ENEold = ENEnew 
        end do

・nstep：

　ステップ数。

・DR_ini：

　DRの初期値。

・DR_max：

　DRの最大値。

・DE_conv：

　|DE|<DE_convになったら収束とみなす。

・dRMS_conv：

　|dRMS|<dRMS_convになったら収束とみなす。

2. 4 初期速度

2. 5 その他

2. 5. 1 系の一部だけを動かすBelly法

1. 　計算時間を短縮する。
   
2. 　系の一部を止めて他の部分だけ動かして最適化等をする。
   

2. 5. 2 重心運動の凍結

Figure 3. 原子の結合トポロジーの定義

3. 1 共有結合力（bond）

Figure 3.1. 共有結合

3. 2 結合角（angle）

Figure 3.2. 結合角

　ただし、

　である。

3. 3 捻れ角（torsion angle）

Figure 3.3. ねじれ角

　上式の捻れ角によるポテンシャルは

　とあらわされる。ここで、与えられたパラメータはVがbarrier height、hが周期、gが位相（０かp）である。このとき、原子i,j,k,lにかかる力は以下の４つの式で計算される。

　ただし、f_a、f_b、 f_cは

　である。この式では、faはsinF≒０のとき計算できない。だが、sing＝０（∵g= 0 or p）を用いて、h＝１〜４の場合にのみ求めると、

　として計算できる。
　なお、i j k lが、この順序に並んでいない場合は、improper torsionと呼ばれているが、同様の式で計算できる。

3. 4 ファン・デル・ワールス（VDW）力および水素結合

Figure 3.4. VDWポテンシャル

　VDW力はFigure 3に示された結合トポロジーでの距離が1-5の位置以上に遠い原子間に適用される。共有結合していない原子ペアは、近い距離だと強く反発し合うが、遠くになると小さな力で引き合う（Figure 3.4)。これをファン・デル・ワールス（VDW）力と呼ぶ。この力はLennard-Jonesが提唱した以下のような簡単なポテンシャルで近似されることが多い。

　この式で、引力項が６乗なのは量子力学計算に基づいているが、反発項が12乗であることに理論的必然性はない（本来は、引力項のみをVDW力と呼ぶのだが、生体分子の力場では反発項も含めてVDWと呼ぶことが多い）。力は、

　となる。
　なお、古いAMBER力場では(Weiner et al. [22, 23]）、１対の原子ペアについて、VDW力か水素結合力かのどちらか一方を計算することになっている。水素結合のポテンシャルは下の3.4(3)式を使っている。

　VDW力や水素結合力は、一般的には原子間距離8-15Å程度で切断してしまって構わない。切断方法は3.5.2節参照。
　PEACHでは、計算機としてGRAPEを使う場合はVDW力とクーロン力は別々に計算するが、汎用機を使う場合、VDW力とクーロン力は同じペアリストを使って（3.5.2節）ひとつのサブルーチンで同時に計算している。

3. 5 クーロン力（Coulomb、Electrostatic、静電相互作用）

3. 5. 1 非周期境界条件における直接和

3. 5. 2 周期境界条件および非周期境界条件におけるカットオフ計算

Figure 3.5.2. Residue-based cutoff

　何ステップかに一回、どの原子とどの原子が相互作用するかを残基間の距離に基づいて調べてペアリストを作っておき、そのリストを使って力を計算する。PEACH では、以下のように作成する。MDを開始する前に、それぞれのアミノ酸残基や核酸塩基の原子の中で、残基の重心にもっとも近い原子を調べておく。その原子の座標を用いて、50 fsに一回ほどの頻度で残基―残基間ペアリストを作る（最終的に作りたいリストの切断半径よりも10Å以上長い切断半径を使う）。その残基―残基ペアリストを利用して、最終的に利用する原子―原子ペアリストを10 fs に一回ぐらいの頻度で作る。この方法により、ペアリスト作成時間を短縮できる。
　なお、VDW力や水素結合力（3.4）も同じペアリストを使って、クーロン力と同じサブルーチンで計算する。
　周期境界の場合、距離を最近接粒子間距離(2.2.1）によって測ること以外は非周期境界と同様である。ただし、切断半径は、箱の一辺の長さの半分以下にしなければならない。

3. 5. 3 周期境界条件におけるエワルド和 -表式 -

3. 5. 4 周期境界におけるエワルド和 -条件決め-

3. 5. 5 非共有結合力の計算法について

Figure 3.5.5. クーロン力の計算

　クーロン力を切断せずに計算する場合、非周期境界の場合は式3.5.1(1)、周期境界の場合は式3.5.3(1)を計算することになる。粒子数をNとすれば、前者はO(N²)、後者はO(M N²)（Mは虚像の数で、理想的には∞）の計算量が必要となる。このような計算をまともにしていたら、計算機がどんなに速くても足りない。そこで、Figure 3.5.5のように、クーロン力を近距離の力と、遠距離からの力にわけて計算することになる。
　カットオフ法は、遠距離からの寄与を無視した場合である。また、エワルド法の場合、（数学的に正しい表現ではないが）実空間の計算が近距離力、波数空間が遠距離力に対応する。
　近距離の寄与の計算は、決められた距離内に存在する全ての原子間について、クーロン力の直接和を計算する。普通はVDW力と一緒に計算するため、最低8Å、余裕を持って15Å以内の距離にある原子ペアについて力を計算する必要がある。時間刻みは、1-2 fs程度であり、これ以上長くできない。
　一方、遠距離の粒子からの寄与の計算は、様々な工夫がされている。フーリエ変換（エワルド法）や多重極子展開（PPPC法やFMM法）を利用して空間的にまとめることで、計算量を短縮できる。また、遠距離力は近距離力に比較して長い時間刻み幅をとることができる（多時間刻み幅法、4.2節参照）。空間的にまとめる方法と時間刻みを長くする方法は併用できるので、遠距離からの力の計算時間は精度を下げずにかなり短縮することができる。もちろん近似としての限界はある（Bishop et al. [30]）。
　近距離粒子間の力は直接和で計算せざるを得ず、高速化は難しい（もちろん、並列化などでの高速計算は可能である）。よって、ある系を8ÅカットオフでMDをした場合と、高速化手法で遠距離力を正確に計算した場合では、原理的に、後者は前者の計算速度を越えることができない。また、高速化アルゴリズムにより、計算時間はO(N^3/2)、O(NlogN)、O(N)等になるが、そのオーダーの係数は一般に大きく、粒子数が少ないと効果がない。高速化アルゴリズムを使うのならば、自分の扱いたい系の大きさと必要な精度を考えて、計算速度と精度の上昇が見込めるようなアルゴリズムを選択すべきである。

3. 6 1-4非共有結合力（1-4VDW、1-4Coulomb）

Table 3.6 1-4非共有結合のスケールファクター

力場	AMBER84 [22]	AMBER86 [23]	AMBER95 [24]	OPLS [25]
SC14VDW	2.0	2.0	2.0	8.0
SC14ELC	2.0	2.0	1.2	2.0

3. 7 排除原子の扱い方（Excluded atoms）

3. 7. 1 VDW力およびクーロン力の直接和の場合

3. 7. 2 クーロン力のエワルド和の場合

4. 1 速度Verlet法（Velocity-Verlet method）

        do i = 1, nstep
                v = v + Dt/2 * f/m
                r = r + Dt * v 
                call calforce(f)        ! 力を計算する 
                v = v + Dt /2 * f/m 
        end do

4. 2 多時間刻み幅法 RESPA（REference System Propagator Algorithm）

        do i = 1, nstep
                v = v + Dts/2 * fs/m
                do j1 = 1, L
                        v = v + Dtm/2 * fm/m
                        do j2 = 1, K
                                v = v + Dth/2 * fh/m
                                r = r + Dth * v
                                call calforce(fh)      ! 固い力の計算
                                v = v + Dth/2 * fh/m
                        end do
                        call calforce(fm)      ! 中間の力の計算
                        v = v + Dtm/2 * fm/m
                end do
                call calforce(fs)      ! 柔らかい力の計算
                v = v + Dts/2 * fs/m
        end do

5. 1 NVEアンサンブル（ミクロカノニカル ）

5. 2 NVTアンサンブル（Woodcock）

5. 3 NVTアンサンブル（Nose-Hoover）

        do i = 1, nstep
                ph = ph  + Dts/2 * Fh
                h = h  + Dts * ph /Q
                v = v * exp(-Dts/2 * ph /Q)

                ************ begin V-Verlet *************
                do i = 1, nstep
                        v  = v + Dt/2 * f/m
                        r  = r + Dt * v
                        call calforce(f)      ! calculate force
                        v = v + Dt /2 * f/m
                end do
                ************* end V-Verlet ***************

                v = v * exp(-Dts/2 * ph /Q)
                call calfehta(Fh)        ! calculate Fh
                                                 ! (=S v2/m - 3NkT) 
                ph = ph  + Dts/2 * Fh
        end do

Figure 6. MDによる物理量の求め方

　MDのトラジェクトリーから得られる情報は様々であり、全てを記すことはできない。巻末の教科書や文献を参照して欲しい。PEACHに導入した解析方法についてはKomeiji et al. [14, 26, 28, 36]に記述しておいたが、重要なものについて説明する。
　以下の説明では、ある物理量Aの時間平均を <A>と書く。すなわち、

　である。ただし、添字jは時間ステップを示し、N_stepは全時間ステップ数である。また、系の粒子数をN、個々の粒子iの座標、速度、質量は、r_i、 v_i、 m_iであらわす。

6. 1 エネルギーに関する物理量

6. 1. 1 エネルギーの平均、ゆらぎ

6. 1. 2 温度

6. 2 分子構造とダイナミクス

6. 2. 1 根平均二乗変位(Root Mean Square Deviation、RMSD)

6. 2. 2 溶媒露出面積（Accessible Surface Area, ASA)　

Figure 6.2.2. ASA

6. 2. 3 回転半径

6. 2. 4 根平均二乗揺らぎ(Root Mean Square Fluctuation, RMSF)　

6. 2. 5 Distance Map (DM)

Figure 6.2.5. DM

　DM（Distance Map）は、残基間（あるいは原子間）の距離に基づいて、タンパク質などのドメイン構造を解析する方法である。二つの残基間の距離R_ijを以下のように定義する。

　この値に基づいて描いたFigure 6.2.7のような二次元の地図をDMと呼ぶ。
　この図は、結晶構造やNMR構造に基づいて描いたものと通常は大差なく、その意味ではわざわざMDをして描く程のものではない。だが、以下のDFMを解析する場合のレファレンスとなる。

6. 2. 6 Distance Fluctuation Map (DFM)

6. 2. 7 動的相関図(Dynamical Cross Correlation Map, DCCM)

Figure 6.2.7. 動的相関

6. 2. 8 Essential Dynamics および Principal Component Analysis