フラグメント分子軌道計算ソフトウェア ADBS(Advance/BioStation) の解析事例

フラグメント分子軌道(FMO) 計算を活用した受託解析サービス（相互作用解析、タンパク質の安定性解析、ドッキング解析など）やソフトウェア導入のご相談は、下記の問合せフォームより、お気軽にお問い合わせください。

HIV-1プロテアーゼとリガントとの相互作用解析

HIV-1プロテアーゼは、HIVウィルスの増殖過程で重要な働きをする酵素タンパク質です。ロピナビルは、このHIV-1プロテアーゼに結合し、その働きを阻害することで、HIVウィルス増殖を抑える医薬品です。ここでは、HIV-1プロテアーゼダイマー-ロピナビル複合体（PDB code 1MUI;198アミノ酸残基+ロピナビル分子、全3224原子から構成）のフラグメント間相互作用解析によって、それらの結合様式を解析した事例についてご紹介します。

計算方法には、MLFMO法を利用し、ロピナビルとその近傍アミノ酸残基との相互作用解析は、MP2法に基づいて行いました。その結果、ロピナビルはプロテアーゼダイマー中のアミノ酸残基Asp25(A-chain、B-chain) およびAsp29(B-chain) とは静電相互作用、Ile50(A-chain、B-chain) とは主にファンデルワールス相互作用により、強く結合的に相互作用することが分かりました（図参照）。

このような医薬品-アミノ酸残基間の結合に関する詳細な情報は、ファーマコフォアの同定や活性の向上を目指した医薬品設計に活用することができます。

Trp-Cageの構造安定性解析

フラグメント間相互作用エネルギー情報の活用例として、ミニタンパク質Trp-Cage (PDB code 1L2Y) を題材に、Trp残基の側鎖配向構造の安定性解析をご紹介します。

このタンパク質は、Trp6を中心に周りを他のアミノ酸残基が取り囲んだ構造をとっています。タンパク質の核となるTrpの側鎖の向き（側鎖ロータマー）はタンパク質内でほぼ一意に固定されています。では、Trp側鎖はどのような理由で安定な側鎖ロータマーを保つでしょうか？ Trp側鎖がある決まった配向構造をとるのは、周辺から特定の分子相互作用を受けているためです。このようなタンパク質の側鎖配向の安定化要因の情報は、タンパク質の立体構造の特性や機能を理解するのに役に立ちます。 では、その答えをADBS(Advance/BioStation) を使って調べてみることにしましょう。

上図に示したのは、FMO-MP2/6-31Gレベルで計算したTrp-CageのTrp6残基と周辺アミノ酸残基との相互作用エネルギーです（括弧内のエネルギー単位はkcal/mol）。Trp6側鎖は電荷的に中性なインドール環で，Trp側鎖の構造を決定する主要因は図に示した7個の側鎖近傍アミノ酸残基との相互作用となります。Trp6と近傍アミノ酸残基とのフラグメント間相互作用エネルギーは全て負値となっており、Trp6側鎖は周辺アミノ酸残基より引き合っていると解釈できます。これはTrp-CageがTrp6を核とした立体構造をとっていることと理にかなった結果になっています。この計算から、Pro17はTrp6側鎖と最も強く相互作用しており、Trp6の側鎖の構造安定に最も重要なアミノ酸残基であることが分かりました。

このように、ADBS(Advance/BioStation) を活用することで、タンパク質のアミノ酸残基の側鎖が周辺アミノ酸残基との分子相互作用のどのようなバランスによって固有の配向構造を保つのか定量的に理解できます。

脂質輸送タンパク質のドッキング解析

タンパク質－リガンド間のドッキング解析におけるADBS(Advance/BioStation) の活用事例として、イネの脂質輸送蛋白質nsLTP1 (PDB code 1UVB)と脂質の一種であるパルミチン酸(PAL) とのドッキング解析をご紹介します。

nsLTP1はリン脂質、糖脂質、脂肪酸などの脂質分子と非特異的に結合し、細胞中の膜の間で脂質を輸送する機能をもつ蛋白質です。脂質は疎水性の長鎖炭化水素をもっており、炭化水素部位と蛋白質との結合はCH/π相互作用やファンデルワース力と呼ばれる弱い相互作用が支配的となります。この相互作用は分散力という量子力学的な挙動に由来し、相互作用の強さの正確な評価には高精度な量子化学計算が必要となります。このような弱い相互作用の結合形態を経験的な関数としてすべて網羅するのは困難なため、炭化水素鎖をもつ脂質分子と蛋白質とのドッキング構造予測は、経験的評価関数に基づいたドッキングソフトが苦手とするもののひとつです。

ADBS(Advance/BioStation) は、量子化学計算によって分散力を高精度にエネルギー評価できますので、ドッキングソフトが結合評価を苦手とする脂質のようなリガンドでもドッキング構造の結合の強さの大小を精密に求めることが可能です。

中性pH条件下でのnsLTP1とPALのドッキング構造をドッキングソフトを使って計算すると、ドッキングソフトが評価したスコア順位の第1位～第5位までのPALの各ポーズ（モデル候補）は図1のようになりました。

ベストポーズ（第1位候補）は、X線結晶構造解析で得られた実験結果のPALの配座（図1赤色）と異なり、さらに結合する位置も違っていました。この計算結果の別のポーズを調べると、スコア順位が第2位のポーズ（図1青色）は実験で得られたPALの結合状態と良く似ていました。また、第3～5位のポーズは、実験の配座と異なっていました。従って、このドッキング計算では、ベストポーズは実験結果の正しいPALの配座を再現しておらず、むしろ第2位のポーズが最も正解に近い配座を予測していた言えます。

これらのドッキング構造に対して、ADBS(Advance/BioStation) を使ってFMO-MP2/6-31Gレベルで計算した結合エネルギー結果を図2に示します。結合エネルギーが最も強い（最も大きい負値）のは第2位のポーズとなりました。この蛋白質の場合、ADBS(Advance/BioStation) のFMO計算で最も結合エネルギーの強いドッキング構造を探すと、実験結果を近い構造が得られることになります。結合エネルギーが強いほどタンパク質とリガンドとの結合はより安定であることを意味しています。現実のタンパク質とリガンドの結合はエネルギー的に最も安定な構造をとっていると思われますので、量子化学計算で求められた強い結合エネルギーの構造が実際の構造を反映しているのは、理にかなった結果であると言えます。

FMO計算による結合エネルギーで評価した結果が実際のドッキング構造と一致していたのは、量子化学計算により分散力を高精度に計算し、正確なPALとnsLTP1との結合エネルギーが得られたためであると考えられます。

このように、ドッキングソフトで計算された構造の再評価にADBS(Advance/BioStation) を活用すると、ドッキングソフトで予測されたベストポーズが本当に正しいかどうか判断する基準が得られ、より正しいドッキング構造を量子論の視点から探すことができます。

Natural Population Analysis (自然ポピュレーション解析)

分子の各原子の電荷を正しく評価することは、分子の電子状態、分子同士の電荷移動量の見積りや分子相互作用など、分子物性の正しい理解のために必要です。しかしながら、原子の電荷は量子力学的に観測可能な量ではないため、量子化学計算で得られた電子密度から一意に原子電荷を決めることはできません。このため、原子電荷の求め方には、様々な方法が提案されています。Mulliken（マリケン）電荷解析は、幅広く使われている電荷密度分割法です。この方法は最も簡便ですが、場合によって、定性的な解釈と相反する値を生じることが知られています。これに対して、最近利用されているのが、Frank Weinhold教授が提案したNatural Population Analysis (NPA) です [1]。NPAは、電荷密度から求められた各原子の電子占有数を重みにして原子軌道(Atomic Orbital, AO) を直交規格化させた軌道(自然原子軌道 Natural Atomic Orbital,NAOと呼びます) をベースに、原子の電荷を見積もる方法です。NPAによる原子の電荷(NPA電荷) は定性的な解釈と合うことが多いため、よく活用されています。

NPAによる電荷解析の事例をご紹介します。下図に示すのは、水溶液中のカルシウムイオンのモデル、カルシウム＋水分子クラスター [Ca(H₂O)₂₉]²⁺ です。定性的には、水溶液中のカルシウムはイオン化され、2価のカチオン(Ca²⁺) として存在していると考えられます。

図.　カルシウムイオン＋水分子クラスター [Ca(H₂O)₂₉]²⁺

表は、ADBS(Advance/BioStation) でFMO-HF/6-31G*レベルで計算した[Ca(H₂O)₂₉]²⁺のカルシウムイオンとその極近傍に存在する5個の水分子のMulliken電荷とNPA電荷です。カルシウムイオンのMulliken電荷は+0.675で、形式電荷+2と大きな差があります。Mulliken電荷で結果を解釈すると、この分子系は+1.325に相当する大きな正電荷がカルシウムイオンから水分子へ流れて、近傍の水分子のMulliken電荷はどれも+0.1e程度の正電荷を帯びていることになります。このカルシウムイオンのMulliken電荷の値は定性的に不自然で、水分子への電荷移動量も大きすぎますので、計算結果の理解に悩むところです。一方、NPA解析ではカルシウムイオンの電荷は+1.357となり、形式電荷である+2に近い値となっています。また周辺水分子のNPA電荷も、Mulliken電荷に比べると中性分子に近い値になっており、定性的に納得できる値です。このことから、このモデルの原子分子の電荷の値や、カルシウムイオンから水分子への電荷移動量は、Mulliken電荷よりNPA電荷で評価するほうが、理にかなっていると言えます。

表. [Ca(H₂O)₂₉]⁺²の電荷解析。FMO-HF/6-31G*レベルの計算。

ADBS(Advance/BioStation) のNatural Population Analysis機能は、FMO計算ならびに通常の量子化学計算に対応しています。また、MP2電子密度にも対応しています。蛋白質の各アミノ酸残基やリガンドなどのNPA電荷を求め、FMO計算で得られる相互作用エネルギーの解釈にご活用できます。

電荷解析は電子密度を各原子へ分配する方法ですが、それに類似して電子エネルギーを各原子へ分配することも可能です。これはエネルギー密度解析(Energy Density Analysis, EDA) [2,3]と呼ばれています。エネルギー密度解析から、分子の電子状態の特徴をより深く理解することができます。ADBS(Advance/BioStation) では、原子軌道(AO) と自然原子軌道(NAO) に基づいたエネルギー密度解析が可能です。

[1] A.E. Reed, R. B. Weinstock, and F. Weinhold (1985) J. Chem. Phys. 83, 735.

[2] H. Nakai (2002) Chem. Phys. Lett. 363, 73.

[3] T. Baba, M. Takeuchi, and H. Nakai (2006) Chem. Phys. Lett. 424, 193.

RI (Resolution of the Indentity) 計算機能

量子化学計算の中で最も時間がかかるのは、電子反発積分(4中心積分) と呼ばれる積分の計算です。RI (Resolution of the Identity) 法は、電子密度関数ρを補助基底関数χで展開し、4中心積分の評価を2つの3中心積分と１つの2中心積分の計算へ分解する方法です。通常、電子反発積分に占める計算コストは、分子系の基底関数の数Nの4乗 - O(N⁴) - に比例しますが、RI法を用いることによってO(N³) となり、大幅な計算速度の向上が期待されます。

蛋白質など生体分子内では水素結合やCH/π結合など弱い相互作用である分散力が数多く見られ、この分散力が生体機能の微妙な制御に大きな役割を持っていることが知られています。このため、分散力を記述できるMP2法による計算の需要が高まっています。MP2計算におけるMP2補正エネルギーを求めるには、莫大な数の電子反発積分の計算が必要となり、計算時間が非常にかかります。一般に、MP2補正エネルギーの計算時間スケールは基底関数Nの5乗 - O(N⁵) - ですが、RI法を適用すると(RI-MP2計算と呼ばれます）、O(N⁴) になり、計算時間が劇的に短縮されます。また、RI-MP2法は、通常のMP2法に比べて、計算機に要求されるメモリーの量が少ないため、計算機資源の制約でこれまで計算できなかった大きな分子系でのMP2計算が可能となります。

新バージョンADBS(Advance/BioStation) Version 3.2では高速なRI-MP2計算エンジンを組み込み、また、FMO計算(FMO-RIMP2)にも対応しておりますので、MP2計算による高精度な相互作用エネルギー計算をより短時間で行うことが可能となりました。さらに、Dimer-ES近似計算と呼ばれる遠距離フラグメントペアに対するクーロン積分の評価にRI法を適用することによって、FMO計算のさらなる高速化も可能となりました。なお、ADBS(Advance/BioStation) では、RI法の補助基底関数に、LCAD (Linear Combination of Atomic electron Distribution)法[1]を利用しています。

[1] S. Ten-no and S. Iwata, (1996) J. Chem. Phys. 105, 3604.

RI-MP2計算の計算速度パフォーマンス

RI-MP2計算による計算事例をご紹介します。インフルエンザウイルスの細胞膜上に存在するノイラミニダーゼ(Neuraminidase) は、ウイルスが感染細胞膜を通過する際に機能する酵素です。商品名タミフルとして有名なオセルタミビル（Oseltamivir）は、この酵素の機能を阻害し、ウイルスを感染細胞膜内に閉じこめ、ウイルスの感染拡大を防ぎます。ノイラミニダーゼ(PDB code 2HT7) はβストランド構造を多く含む球状蛋白質（アミノ酸残基数368）で、リガンドとしてタミフルが結合しています。この蛋白質に対して、FMO法による通常のMP2(FMO-MP2) 計算と、RI-MP2(FMO-RIMP2) 計算を実行し、両者の計算時間の比較を行いました。蛋白質の原子数は約5600（水素原子を含む）で、基底関数系は6-31G を用い全系の基底関数の数は約33000です。計算はAMD Opteronクラスター(16コア) で行いました。

FMO-MP2では、フラグメントモノマー計算で約2時間、フラグメントダイマー計算で約98時間かかり、合計で計算時間は約100時間かかりました。用いた計算機資源では、蛋白質のMP2計算に4日以上かかることになります。一方、FMO-RIMP2計算では、FMO-MP2計算に比べて、10倍程度高速化され、計算時間はわずか10時間と、半日足らずで計算が終了しました。これまで、FMO法によるMP2計算では、たくさんの計算リソースを必要とすることがボトルネックになっていましたが、RI-MP2法を活用することで非常に短時間で計算結果が得られるようになり、研究開発への活用が期待されます。

密度汎関数理論(DFT) とは？

最近、量子化学計算を用いた分子研究では、密度汎関数理論(Density Functional Theory, DFT) に基づいた計算が数多く行われています。量子化学計算でDFTがメジャーになったのは、Hartree-Fock(HF) 法では考慮できない電子相関効果をDFTで含めることができ、なおかつ、DFTの計算時間スケールはHF計算とほぼ同等という点にあります。実際、B3LYPのようなDFT計算で得られた分子物性・化学反応は実験データとよく合うことが数多く報告されており、ローコストでハイリターンな結果が得られるDFT計算は、特に実験研究者に人気があります。また、金属を含む分子系はHOMOとLUMOのエネルギーギャップが小さく、MP2計算では不安定な解に陥る場合がありますが、DFT計算ではこのような場合でも安定な解が得られるメリットもあります。

HF理論は電子間のクーロン相互作用を平均場近似に置き換えますが、この近似では電子間相互作用を厳密に評価することができません。この平均場近似を越えた電子間相互作用の影響は電子相関効果と呼ばれています。電子相関効果は分子物性に重要な役割をもつ場合が多く、たとえば、分子相互作用の分散力はまさに電子相関効果によるものです。したがって、より精密に実験値と比較したい場合、電子相関効果を含めた計算がどうしても必要になってきます。電子相関効果を扱うには、MP2などの多体摂動論・配置間相互作用理論・結合クラスター理論などが必要ですが、HF計算に比べると膨大な計算資源と計算時間を必要とし、大きな分子系では実用的ではありません。一方、DFTは計算スケールがHF計算とほぼ同等で、しかも、HF計算より電子相関効果を含めた結果が得られます。このため、MP2計算ができないほど大きな分子系の場合でも、DFT計算により信頼性の高い分子物性の情報を得ることができます。

DFTは、基底状態の電子密度からあらゆる物性が求められることを立証したHohenberg-Kohn定理を基礎におく多電子系理論です。これは、HF法やMP2法などの波動関数理論とは異なる理論体系をなしています。Kohn-Sham法はDFTで最もよく使われる方法で、電子エネルギーを運動エネルギー・核引力エネルギー・クーロン相互作用エネルギー・交換相関エネルギー項に分けて計算します。このうち、運動・核引力・クーロン相互作用エネルギー項は、HF理論と類似して、電子密度による汎関数型が与えられています。しかしながら、交換相関エネルギー項の厳密な汎関数型は現在わかっておりません。もし、交換相関エネルギー項が厳密に求まれば、Hohenberg-Kohn定理が保証するように、電子相関効果を含めた厳密な電子エネルギーが求められるはずです。このため、実験データを利用して実際の電子相関効果が正しく含まれるよう様々な分子で検証された近似的な交換相関汎関数が、数多くの研究者によって考案されており、今もなお工夫改良がなされています。さらに、B3LYPのような、HF理論の交換エネルギーを交換相関エネルギーに混ぜるハイブリッド型も提案されています。このように、DFT計算は、ある意味、実験値とあうように経験的に求めた交換相関汎関数を用いる点がありますので、"非経験的(ab initio)"にシュレディンガー方程式を解く、純粋な波動関数理論とは趣が違う面があります。このような事情から、DFT計算の分子物性値は、MP2計算と比べると、より実験測定値に近いものが得られるケースがあることが知られています。

新バージョンADBS(Advance/BioStation) Version 3.2ではDFT計算機能を取り入れました。もちろん、FMO計算にも対応しています。ラジカルなどの開殻系での計算も可能です。これにより、DFTレベルでの相互作用エネルギー解析が可能となりました。B3LYPやSVWNなど、Gaussian・GAMESSなどで用いられる主要な汎関数にほぼ対応しています。なお、通常、交換相関エネルギー項は解析的に解けないため、数値積分を用いて求めますが、ADBS(Advance/BioStation) では、この数値積分計算にMura-Knolesの動径グリッド[1]ならびにLebedevの立体角グリッド[2]を採用しています。

[1] M E. Mura and P. J. Knowles (1996) J. Chem. Phys. 104, 9848.

[2] V. I. Lebedev (1992) Russian Acad. Sci. Dokl. Math. 45, 587.

FMO-DFTの計算事例をご紹介します。図はグリシンが5つ結合したグリシンペプチド鎖（グリシン5量体）です。DFTの汎関数にB3LYP、基底関数に6-31Gを用いて、FMO-DFT(FMO-B3LYP) 計算を行いました。CPU数を変えて計算を行い、並列計算による計算速度の高速化（高速化率）を調べました。

表は、グリシンペプチド鎖のFMO-DFT計算の高速化率です。シングルCPUで計算したときに比べて、デュアルCPU(2コア)では計算時間が約2.5倍速くなっており、さらに16CPUの計算でも、シングルCPUの計算に比べて20倍近く速くなっています。このように、ADBS(Advance/BioStation) には非常に高速なFMO-DFT計算エンジンが組み込まれており、PCクラスターのような並列環境で、たいへん効率よい計算が実行できます。

ADBS(Advance/BioStation) で対応している密度汎関数

ADBS(Advance/BioStation) で対応しているDFTの主な密度汎関数を、Gaussianのキーワードと合わせて、ご紹介します。

Advance/ BioStation の 汎関数キーワード	説明	Gaussianの 汎関数キーワード
SLATER	Slater交換汎関数	HFS
BECKE88	Becke 88交換汎関数	HFB
PW91	Perdew-Wang 91交換汎関数	PW91
SVWN	Slater交換汎関数＋Vosko-Wilk-Nusair相関汎関数	SVWN5
SVWNrpa	Slater交換汎関数＋Vosko-Wilk-Nusair (RPA)相関汎関数	SVWN
BVWN	Becke88交換汎関数＋Vosko-Wilk-Nusair相関汎関数	BVWN5
BVWNrpa	Becke 88交換汎関数＋Vosko-Wilk-Nusair (RPA)相関汎関数	BVWN
SLYP	Slater交換汎関数＋Lee-Yang-Parr相関汎関数	SLYP
BLYP	Becke 88交換汎関数＋Lee-Yang-Parr相関汎関数	BLYP
B3LYP	(0.2×Hartree-Fock+0.8×Slater＋0.72×Becke88) 交換汎関数+(0.19×Vosko-Wilk-Nusair (RPA) + 0.81×Lee-Yang-Parr)相関汎関数	B3LYP
BHandHLYP	(0.5×Hartree-Fock+0.5×Slater＋0.5×Becke88) 交換汎関数+Lee-Yang-Parr相関汎関数	BHandHLYP
BHandH	(0.5×Hartree-Fock+0.5×Slater) 交換汎関数+Lee-Yang-Parr相関汎関数	BHandH

その他、ADBS(Advance/BioStation) で対応しているDFTの汎関数キーワードは次の通りです－SPW91・SP86・BPW91・B86・PBE・B3PW91・B3P86・X3LYP・B97・B971・B972・HCTH・HCTH120・HCTH147・HCTH407