電子相関を露わに導入した QM/MM レプリカ交換・分子動力学計算法の開発
萩原陽介，重田育照，舘野賢
筑波大学大学院数理物質科学研究科物理学専攻
mtateno@comas.frsc.tsukuba.ac.jp
概要： 近年，サイズの小さいタンパク質（ペプチド）に対しては，その変性状態から出発した分子動力学 (MD) 計算により，折畳まれた立体構造を理論的に予測する試みがなされている。ところが実際には，わずか 10 個のアミノ酸残基からなるペプチドの折畳み構造すら正確には再現できない場合があり，その主因として電子相関による van der Waals 相互作用が正確に記述できていない点があげられる。本研究では，そうした効果を MD 計算にはじめて露わに取入れる手法を開発し実装することにより，従来よりもはるかに高精度な計算法を開発した。

1

はじめに互作用を従来の MD 計算で正確に取り扱うことの近年におけるコンピュータできない理由でもある。さらには，こうした芳香の高速化に伴い，タンパク質の環同士の T 字型相互作用による立体構造の安定化

折 畳み 過程 ( フォー ルディ ン

は，chignolin に限らずタンパク質構造の至るとこグ)自体が，分子動力学(MD)ろに存在する。したがってその正確な取扱いは，計算のターゲットとなりつつ広く一般にタンパク質の機能解析においても本来ある。現在，タンパク質のサイ必須である。そこで本研究では，先の電子相関のズが～50 アミノ酸残基程度ま効果を露わに取入れた高精度な MD 計算を実現し，では，そうした折畳み過程の理タンパク質の折畳み過程を正確に追跡することを論的解析が行われつつある。ところが計算結果に目指した。

ついては，わずか 10 残基からなるタンパク質でさ2 理論・理論・手法

えも，実験的に得られた立体構造を MD 計算によ先の vdW 相互作用による芳香環同士の安定化って正確には再現できない場合もあるのが実情でエネルギは，単純な Hartree-Fock (HF) 法や密度ある。例えば，上図の chignolin の場合には，その汎関数法 (DFT) では評価することが不可能であ

分子内部に存在するふたつの芳香環 (Tyr2 およびり，高次の摂動法や Coupled Cluster (CCSD) 法なTrp9) が T 字型の結合を形成して，一定の立体構

ど，電子相関を高精度に考慮することの可能な手造を保持していることが NMR 分光法によって明法を用いなければならない。そこで本研究では，らかになっている (図 1)。しかし MD 計算によれ

van der Waals 相互作用をより効率的に求めるたば，これらの芳香環が非常に大きく揺らぎ，それめに，DFT（相関交換項に vdW 相互作用がまったらの位置関係は一定に保持されず，その結果ふたく含まれない）で得られた電子密度から vdW 相互つ以上の主たる立体構造が現れる。こうした理論作用を補正する，以下の ALL 汎関数を用いる [1]。と実験結果との深刻な矛盾は，強相関電子系から生じる van der Waals (vdW) 相互作用による芳香環Evdw = −6ρ 1 2 (r ) ρ 1 2 (r2 ) 1d 3r1 d 3r2 1 2 1f (r )3 2 ∫∫(4π )ρ (r1 ) + ρ 1 2 (r2 ) r12 6 damp 12

場におけるハミルトニアンには露わには記述され  α 6 f damp = exp −  AB    r12  ていない。他方で，こうした系を量子力学的に取ρは全電子密度，fdamp は発散する距離間にある電同士の安定性に起因するためであり，現在の古典り扱うには高度に高精度な理論を必要とし（後述），計算量の増大を招く結果となる。これが，先の相子からの寄与をカットオフするためのダンピング関数，αAB は原子種ごとに定められたダンピングの強さを制御するパラメタである。DFT 計算に

ALL 汎関数による補正を加えた計算は，ベンゼン

10.000等の分子に対して行われており，CCSD 法で得ら
stacking energy (kcal/mol)

8.000れた結果とよく一致することが示された[2]。
ただしこのままでは，MD 計算等において現れるコンフォメーションのすべてに対して，いちいちDFT を適用しなければ電子密度が求められない。

6.000
4.000

平行型(CCSD)

2.000
0.000
3.000
-2.000垂直型(CCSD)
平行型(ALL)
4.000

5.000

6.000

7.000

8.000

垂直型(ALL)-4.000-6.000

そこで先の汎関数をより効率的に用いるため，-8.000distance (angstrom)Tyr2 および Trp9 のそれぞれの電子密度を下式のガウス関数で近似することにより，MD 計算の前にあらかじめ電子密度を求めておき，これを用いて芳香環同士の vdW 相互作用を ALL functionalを用いて記述することを試みた。

ρ (r ) = qiπ−3 2(ai exp − ai (r − ri )
322)

まとめ
生体高分子内に頻繁に存在する芳香環同士の相互作用による立体構造の安定化は，強い電子相関から生じる vdW 相互作用に起因する。そうした相(π 2)互作用は古典場で露わに記述されておらず，それRiが従来おこなわれてきた MD 計算の誤りの主たる

12

ai =4

ここで，qi は原子 i の総電子数，Ri はその vdW 半原因になっている。そこで本研究では電子相関の径であり，これは Generalized Born 法における表効果を ALL functional を用いて露わに取入れた，式 と 一 致 し て い る [3] 。 つ ま り chignolin の従来よりもはるかに高精度な計算法の開発を行っTyr2-Trp9 の vdW 相互作用については，その電子た。これを基に今後レプリカ交換 MD シミュレー相関の効果を上記の形式によって露わに取入れるションを実行することで，極めて正確なフォールことにより高精度に記述しようというものである。

ディング過程の解明が期待される。

その他の相互作用は従来の古典場 (AMBER ff99)を用いて記述する。

3結果
図 2 のように Tyr2 とTrp9 の芳香環同士を互いに平行 (図 2 右) およびT 字型 (図 2 左) に結合させたモデル構造を構築し，ふたつの芳香環同士の距離を様々な値に設定し，CCSD 法により（重心間の距離に対する）ポテンシャル曲線を計算した。次に，このポテンシャル曲線を再現する電子密度分布を，上記のガウス関数のパラメタ ai を変数として Simulated Annealing(SA) 法により探索した。これによって得られたポテンシャル曲線を元に，さらにダンピング関数
fdamp のαAB および古典場におけるレナード・ジョーンズ (LJ) ポテンシャル反発項のふたつの係数を，SA 法の探索空間に含めて再度，最適化を行った。

以上の結果，両モデル(平行型および T 字型)において CCSD 法で得られたポテンシャル面をよく再現する電子密度分布およびパラメタ・セットを得ることに成功した (図 3)。これを基に ALL
functional 由来の力を求め，MD 計算を行うためにそれらを AMBER9 パッケジ内に実装した。

参考文献[1] Y. Andersson, D. C. Langreth and B. I.
Lundqvist, Phys. Rev. Lett. 76,76 102 (1996).
[2] T. Sato, T. Tsuneda, K, Hirao, J. Phys. Chem.
123,123 104307 (2005).
[3] M. Schaefer, M, Karplus, J. Phys. Chem. 100,001578 (1996).