密度汎関数法に基づくタンパク質全電子計算プログラムProteinDFの現状と今後の展望

佐藤 文俊


Return

1 はじめに

タンパク質の働きは多岐に渡る。ミトコンドリア内膜の呼吸鎖では電子伝達やプロトンの移動を行っているし、光合成を行うバクテリアや樹木では光励起や光異性化が重要である。薬剤代謝に大きな役割を果たすP450(CYP)が行っているのは化学反応である。しかもこれらの反応は生理温度という穏やかなエネルギーで進む。自然は実に精巧な仕組みを私たちに授けたが、逆に言えばその解析にはkcal/molの精度が要求される。しかも、タンパク質の約3分の1はアミノ酸が連なってできたペプチド鎖以外に、活性中心や構造安定化モチーフとして金属イオンを持っている。タンパク質の機能解析にはアミノ酸、金属イオンを同程度の精度で扱える手段が必要である。
反応の主役である電子を取り扱い、かつ精度の高い解析を行うためには、タンパク質全体を取り扱う大規模な量子計算によるシミュレーションが必須である。1993年1月、九州工業大学情報工学部で、柏木浩九州工業大学名誉教授と筆者が中心となり、タンパク質の全電子を計算するプロジェクトを発足した。数千〜数万原子からなる巨大分子の量子化学計算と高い計算精度という相反する要求を両立させる方法として計算法には密度汎関数法を、長期にわたるプログラムの改良と拡張にはオブジェクト指向技術をそれぞれ採用した[1]。プロジェクトを初めてから7年半、2000年の夏に9,600次元のヘムタンパク質シトクロムcの全電子波動関数の計算がようやく達成された[2, 3]。計算機は15台のワークステーションを100Base-TX Ethernetで接続した理論ピーク性能が15 GFLOPSのクラスタである。実質計算には1年半を要したものの、試行錯誤の分を取り除いた正味の計算時間は約3ヶ月で、この計算により金属タンパク質でも一般分子と同じ計算方法で電子状態が解けることが実証された。本プログラムをProteinDFと名付けた[4]。
これとほぼ同時期にハートリー・フォック法によるChallacombeとSchweglerの先駆的なP53の計算(3,836軌道)[5]や、高田らのPKCの計算(8,672軌道)[6]などのタンパク質全電子計算が報告された。また、諸熊らのQM/MM法の外挿法であるONIOM法[7]、北浦らのフラグメント分割を利用したFMO法[8]といった実践的な計算手法も提案された。それぞれのプログラムや計算法は独自に発展を遂げており、これは本特集号で解説されているはずである。様々なソフトウェア、計算手法の充実に伴い、いよいよ生体分子の電子状態計算は実用的な基盤が整いつつある。
ProteinDFも、上記の計算で得られた経験が様々な形でフィードバックされ、2002-2004年の文部科学省ITプログラム「戦略的基盤ソフトウェアの開発」[9]、2005-2007年の次世代IT基盤構築のための研究開発「革新シミュレーションソフトウェアの研究開発」[10]を通して、計算規模・計算速度ともに大幅に向上している。そこで、本稿ではProteinDFシステムの現状を報告し、これから窺える近未来の展望について解説する。

2 ProteinDF

2. 1 手法

ProteinDFはKohn-Sham-Roothaan法に基づく密度汎関数法で独自に発展した高速計算法を採用している。Sambe、Felton [11]が提唱し、Andzelm、Wimmer [12]、Dunlapら[13]が改良した、RI(Resolution of Identity)法である。Hartree-Fock-Roothaan方程式は分子積分演算と一般行列演算からなる行列方程式である。一般行列演算は線形代数計算であり、コンピュータを用いて高速かつ精密に計算することが出来る標準的なライブラリが存在する。分子積分演算も、Boysの提唱[14]によって基底関数にガウス型関数を使用されるようになったため、完全に解析的に計算できる。しかも、小原と雑賀によって非常に高速な計算方法も開発されている[15]。ところが、Kohn-Sham-Roothaan方程式にはこれに加えて交換相関項の計算が必要で、全てを解析的には計算できない。そのため、解析的に計算できるはずの電子相互作用項の計算にも、本質を損なわないあるいは計算精度には悪影響を与えない程度の近似を導入し、これを高速に計算してしまおうという考え方が生まれた。その一つがRI法である。

2. 1. 1 クーロン項に関するRI法

この方法では、Roothaan法による分子軌道の基底関数展開と同様に電子密度にも補助の基底関数 を導入し、

と展開する。電子密度の補助基底関数の展開係数raのもっともシンプルな求め方は、電子密度の分布をフィッティングする方法である。しかし、Dunlapらは、クーロンエネルギーの誤差

が最小になるようフィッティングする方が、計算精度がよくなることを明らかにした[13]。これはクーロン項を

と展開することと等価である。ここで、分子積分は

であり、gp(r)は分子軌道の基底関数、-1乗は逆行列を意味する。実際には、電子密度は全空間で積分すると

全電子数となることから、ラグランジェの未定乗数Lを用いて

とフィッティングさせ、電子密度が持つ性質を保存させる。拘束条件を課さない方法も提唱されている[16]。(3)式から、基本的にはそれでも良いことがわかるが、その場合は一般により多くの補助基底関数の数を要する。

2. 1. 2 交換相関項に関するRI法

上記のクーロン項に関するRI法と同様のアイデアで、交換相関項を計算する方法も開発されている。すなわち、(1)式の電子密度と同様に、交換相関ポテンシャルにも補助の基底関数 を導入する。

交換相関ポテンシャル展開係数は、分布が最も合うようにフィットする。すなわち、

とする。なお、この積分計算にはガウスの求積法による数値計算法が使用される。(7), (8)式から交換相関項は

である。以上より、クーロン項と交換相関項にRI法を適用する場合、Kohn-Sham行列は

と計算される。

2. 1. 3 RI法の特徴

RI法における最大の恩恵は、4中心の2電子積分 が、2回の3中心のクーロン積分 で計算されることである。インデックスが減ることは計算量が減ることを意味する。実際の計算では、カットオフ法や高速多重極展開法を用いるため、4中心積分でも3中心積分も計算サイズ依存性には変わりがない。しかし、サイズ依存性に掛かる係数が格段に小さくなることにより、この方法は4中心を使用した方法よりも一般に数倍から数十倍は高速である。そのため、RI法はHartree-Fock-Roothaan法にも逆輸入されている。Fockの交換項もRI法で計算することができる[17]。
歴史的には、電子密度、交換相関ポテンシャルともにそれぞれの補助基底関数で展開する方法が提唱され、その後、密度汎関数法による分子軌道法では広くその手法が使用された。そのため、密度汎関数法固有の解法とする解説があるが、RI法はクーロン項、交換相関項の計算の一選択肢に過ぎない。また、交換相関項のRI法はハイブリッド法には使えないと考えられる向きがある。しかし、ハイブリッド法は局所交換相関、非局所補正、Fockの交換項の線形和で表現されるため、Fockの交換項を別に計算して足しこめばよい。
RI法の最大の欠点は、(6)式や(8)式のように電子密度や交換相関ポテンシャルをフィットするため、変分原理を満たさなくなることである。このため、RI法を用いた自己無撞着(SCF)計算の収束過程では、繰り返し計算される物理量は徐々に変化しなくなるものの、収束した全エネルギーが初期の値よりも必ず低くなるとは限らない。全エネルギーが上昇しながらSCF計算が収束する事態が頻繁に生じるため、全エネルギーが下降することを以って、SCF解が改善された証拠とすることができない。そのため、RI法では通常の方法以上に様々なSCF計算物理量で収束過程を詳細にチェックする必要がある。

2. 2 並列化

RI法に基づく密度汎関数法で律速となる計算ルーチンは、分子積分演算、交換相関ポテンシャル係数を求める数値積分、そして行列演算である。また、これらの計算ルーチンは行列要素の保持に大量のメモリ領域を要する。例えば、100,000軌道の計算では1つの行列のサイズが倍精度実数で80 GBにのぼる。
分子積分演算は、全ての行列要素を各CPUに分散保存させながら、積分カットオフによる計算プレスクリーニングを組み合わせたダイレクト法によって、計算に要するワークエリアを圧縮してSingle Program Multiple Data(SPMD)を達成するという独自のアルゴリズムで並列化した。例えば、Kohn-Sham行列の計算では、前回のSCF計算の変化量として

と計算される。ここで、Dは前回のSCF計算の変化量である。各分子積分はp, qのインデックスのみによる上限値の判定ができるので、これが計算精度よりも小さければ計算しなくて済む。この判定をあらかじめ行い、実際に(11)式の計算をおこなうタスクをスクリーニングしてしまう。計算結果は の形で保持すればワークエリアを大幅に圧縮できる。また、(11)式タイプのRI法による並列化では、分子積分に積和する量が補助基底関数の数のベクトル量(Dra,Dmg)である点も通信量削減に寄与する。本手法における並列効率向上の鍵は、積分カットオフ法で生き残ったp, qインデックスを軌道シェルタイプの組み合わせ毎に、積分計算を均等分散処理させることである。分子積分演算は(6)式や全エネルギーの計算でも使用されるが、並列化の方法は同じである。
交換相関ポテンシャル係数を決める積分は、交換相関汎関数が有理関数であるため、数値計算で求める必要がある。しかし、この多中心の数値積分は、一中心の数値積分の和に変換して計算する[18]。そのため並列化は容易で、原子毎にタスクを分割して計算することができる。
行列演算は対角化、行列積、逆行列である。これらの演算には、並列版行列演算ライブラリであるScaLAPACK [19]の密行列用の関数を用いた。本ライブラリも行列要素を分散保持するため、全ての演算は1CPU毎に超巨大メモリを必要としないプログラム構造となった[20]。

2. 3 主な計算結果

ProteinDFを用いて、306残基、原子数4,728、電子数18,552、基底関数の数26,790、補助基底関数の数48,684のインスリン6量体の全電子状態計算を行なった。使用した計算機は64個のItanium2 プロセッサがNUMAlinkTM(3.2 GB/s)で結合されたSGI Altix3700で、ピーク性能は333 GFLOPSである。全電子計算は擬カノニカル局在化軌道による初期値作成法[21, 22]とアンダーソンの収束加速法[23]により、26,790×26,790次元のSCF計算を17回の繰り返しで無事に収束させることに成功した[24]。Figure 1はSCF計算中のKohn-Sham軌道エネルギーの履歴である。本計算における初期値の良さと収束の速さを示すと同時に、この規模のタンパク質全電子計算が安全に実行できている様子が分かる。
全計算時間は65時間で(Table 1)、SCF一回転目では82%の並列化効率が得られた(Table 2)。これは99%以上の並列化率である。133残基のインターロイキン(原子数2,157、電子数8,210、基底関数の数11,909、補助基底関数の数21,751)を用いて、Altix3700とGbit Ethernetで接続したPCクラスタ(ピーク性能は154 GFLOPS)で同様の計算見積もりを行った結果、ネットワークスピードが遅い場合(PCクラスタ)、行列並列演算が大きく影響を受けて全体の効率が下がるが、本手法による分子積分、交換相関フィッティング並列計算に関してはネットワークスピードの影響が少ないことが明らかになった(データは省略する)[20]。本研究により、少なくとも100CPUクラスで有効に並列化されたRI密度汎関数法に基づくカノニカル分子軌道法プログラムが完成した。


Figure 1. History of the orbital energy on insulin hexamer per SCF iteration. The calculation was successfully converged at the 17-th iteration. The distribution of the orbital energy did not significantly change from the 11-th to the converged (17-th) iteration.

Table 1. Calculation time for insulin hexamer using 64 processors on 1.3-GHz Itanium 2 cluster (SGI Altix3700).
Calculation routinesElapse time(hour:minute)
Integrals before SCF loop0:27
X matrix0:56
Grid generation8:27
Inverse matrix0:34
    SCF loop
    Average3:12
    (First)(3:49)
    (Coverged: 17th)(2:36)
Etc0:21
Total64:58

インスリンは膵臓で作られるホルモンタンパク質で、糖尿病の治療には注射によるインスリン投与がよく用いられている。ところが、インスリンは強い自己集合能があるため、薬剤のような高濃度下では6量体がもっとも安定、続いて2量体が安定な構造となる。注射をしてから血液内に吸収されるまでの間に6量体から2量体、そして薬効のある単量体へと解離して、インスリン受容体に作用するため、効くまでに時間がかかってしまう。最近、遺伝子工学を駆使した新しいインスリン製剤が登場した。遺伝子置換によって自己集合能を弱める設計が施された超即効型のインスリン、逆にこれを強めることによって基礎代謝を維持させた持続型インスリンが開発されている。

Table 2. Timing-data in the first SCF calculation of insulin hexamer using 1.3 GHz Itanium 2 cluster (SGI Altix3700).

Figure 2はインスリン単量体QM/MMモデル計算と6量体全電子計算の電子分布の差である。タンパク質の電子はまるで海のように豊潤で、凝集によって波紋を作るように電子が再配置する様子がわかる。このような結果は、当グループが全電子計算にこだわる理由を端的に表している。電子分布が変わってしまうため、クーロンエネルギーの見積もりも大幅に変化する。古典論の計算で用いられている力場は大変良くできているが、最も大きな割合を占めるクーロンエネルギーが変わってしまっては結果の見積もりを誤ってしまうだろう。実際、6量体から3つの2量体に解離するのに必要なエネルギーは、古典論では解離エネルギーをだいぶ過小評価していた(投稿準備中)。さらにシミュレーションを積み重ねる必要があるが、量子化学計算によってそもそも単量体が安定なインスリンが設計できる可能性があり、もしこれが臨床に通るようなことがあれば、打てばすぐ効く超々即効型のインスリンが登場するかもしれない。


Figure 2. Three-dimensional distribution of electron density difference between insulin hexamer and its six monomers model by QM/MM calculation.

3 実用化の鍵:ProteinEditor

3. 1 ProteinDFシステム

当グループは、以下の5つの研究開発項目を統合し、ProteinDFシステムを構築している(Figure 3)。A. 密度汎関数法による超大規模タンパク質全電子計算エンジン: ProteinDF。B. タンパク質の構造最適化、分子動力学、自由エネルギー計算機能:ProteinMD、タンパク質の原子核座標の運動を古典論で記述するプログラムである。Eのインターフェースの機能と連動する以外に、ProteinDFとの連携部分を開発している。現行のab initio分子動力学(MD)計算のレベルは十残基程度である。C. 密度汎関数法−配置間相互作用法による電子励起・電子移動計算機能: ProteinCI、本パッケージは現在開発中である。D. タンパク質波動関数データベース:ProteinQR、現在60種類ほどのタンパク質波動関数データを計算し、その一部を公開している[25]。E. シミュレーションを支援し、物性を評価・解析できる高品位グラフィカル・ユーザー・インターフェース(GUI):ProteinEditor、100残基規模のタンパク質全電子計算を支援するバージョンを公開している[26]。
本システムは、生命現象の基礎過程解明といった基礎研究のみならず、古典論による手法が主であった創薬、生体触媒、生体分子素子、環境物質などへの応用研究に信頼性の高い基盤技術を提供することも目標としている。


Figure 3. Concept of ProteinDF System.

タンパク質のような巨大分子の量子化学計算ではシミュレーションが大規模で複雑になる。煩雑さやミスが増すとともに、これらの制御やプリ・ポスト処理にはより多くの専門的な知識と経験が要求される。これらを軽減するためのインターフェースの役割は実用化の観点から、非常に重要である。そこで、本章では統合環境GUIであるProteinEditorについて解説する。量子化学計算が専門の研究者のみならず、生物物理学や生化学の研究者も想定して開発を進めている。

3. 2 タンパク質全電子計算達成の基本要素

3. 2. 1 タンパク質全電子計算の流れ

Figure 4にタンパク質全電子計算達成までの流れを示す。タンパク質立体構造情報はProtein Data Bank(PDB)[27]に収集・登録されており、多くのタンパク質のシミュレーションで利用されている。しかし、その約8割がX線構造解析で得られた結果であり、水素原子の構造情報を適切に追加する必要ある。アミノ酸置換による比較実験も標準的なシミュレーション手段である(Figure 4には記述なし。Figure 7参照)。これ以外にも、PDBデータには化学的な見地から問題と思われるほど構造に歪みや異常な原子間距離を持っているものが多く見られる。一般にこのような構造は反応性の高い軌道を与えるため、そのままの構造を用いた量子化学計算は結果判断を誤らせる危険がある。構造最適化やモデリングを行うプリ処理機能が必須である。続く実際の計算実行には、タンパク質に適したSCF計算の初期値作成や、計算制御が鍵を握る。最後に出力される膨大な量の計算結果の表示・解析といったポスト処理も重要である。ProteinEditorはこれら全ての機能をサポートしている。


Figure 4. Flow of all-electron calculation on proteins.

3. 2. 2 QCLO法を用いたタンパク質全電子計算シナリオ

タンパク質全電子計算では求める行列要素数が膨大な上に、密度汎関数法では占有軌道と非占有軌道のギャップが狭く、非常に良い初期値からSCF計算を始めないと解が発散してしまうか、正しい解に収束しない。この良い初期値を作成する方法として、局在化軌道(Localized Orbital: LO)を応用した擬カノニカル局在化軌道(Quasi-Canonical LO: QCLO)法を開発した[21]。QCLOは、ある領域(フラグメント)に局在化しているが、その領域内ではカノニカル軌道に似た軌道を与える。QCLO法は小さなペプチド鎖のQCLOをつなぎ合わせて、より大きなペプチド鎖の高品位な初期値を作成するというもので、計算分子(フレーム分子と呼ぶ)を徐々に伸張して、その操作を繰り返し、タンパク質の全電子計算を達成する。


Figure 5. Convergence process for all-electron calculation on proteins.

Figure 5にQCLO法によるタンパク質全電子計算の流れを示す。伸長の諸段階をステップと呼び、まずStep 1では、タンパク質をアミノ酸残基ごとに分割し、各1残基のアミノ酸の計算を実行する。次にStep 2ではStep 1のアミノ酸残基の電子密度を合成した初期値から3残基の計算を行う。これでフレーム分子にペプチド結合が表現される。Step 3ではStep 2で計算された3残基の結果から、QCLOを計算し、それぞれ中央のアミノ酸残基のQCLOをつなぎ合わせて数残基のフレーム分子の初期値を作成し計算を行う。Step 4以降はStep 3と同様の過程でフレーム分子サイズを徐々に大きくして、フレーム分子がタンパク質全体となるまで繰り返す。QCLO法ではStep 3以降のフレーム分子の伸長の仕方やステップ数に任意性があり、これをシナリオと呼んでいる。計算者はタンパク質の構造を考慮に入れて、安全に全電子計算を達成できるシナリオを決める必要がある。

3. 3 ProteinEditorの機能

ProteinEditorはタンパク質全電子計算を達成するための様々な機能、およびそのGUIをサポートしている(Figure 6)。ProteinEditorは全てスクラッチから開発しており、産学を問わず提供することができる。OSはWindowsに対応、分子グラフィックスはOpenGLを用いて描画している。タンパク質の立体構造はPDB、グラフィックスデータの一部はAVS[28]のフィールドデータに準拠している。


Figure 6. Overview of ProteinEditor.


Figure 7. Example of mutation.

3. 3. 1 タンパク質計算構造の作成支援機能

PDBデータへの水素付加機能やMD計算と連動させて構造歪みを緩和する機能、および構造に問題がないかをRamachandranプロットなどで構造妥当性のチェックを支援する機能を実装した。なお、現在MDには古典法を採用している。水素付加は、PDB hetero group dictionaryの情報に基づいていおり、アミノ酸残基や核酸以外のヘテロ分子もPDBに登録されている構造には水素を付加することが可能である。
タンパク質内のアミノ酸残基置換やメチル化、脱アミノ化といった化学修飾を行うこともできる。Figure 7にアミノ酸残基置換と構造最適化を行った例を示す。(a)は1つのアミノ酸残基を選択している。(b)ではSequenceViewerによって(a)で選択されたアミノ酸残基を置換している。 (c)は置換後の構造を表し、原子間衝突が発生している様子が分かる。この結果に対して構造最適化を適用した結果が(d)である。原子間衝突が回避されている様子を確認できる。また、この構造が化学的に妥当な構造になっているかは、構造妥当性評価機能でチェックすることができる。

3. 3. 2 自動計算法支援機能

全電子計算のためのQCLOシナリオの自由な編集とそのジョブ実行を支援する機能としてScenarioEditorが実装されている(Figure 8)。自動シナリオ作成機能では共有結合が過不足なく考慮された基本シナリオが自動的に生成される。このシナリオは経験的に求めた値でピッチを固定した等差級数的なものであるが、数十残基程度のタンパク質の場合、十分安全に収束する。しかし、タンパク質は折り畳まれることで、一次構造上は離れていても高い相互作用を示す部位が多々存在する。このような場合はタンパク質が大きくなるほど顕著となり、安全に収束させるためには的確にシナリオを修正する必要がある。ScenarioEditorでは、マウスを用いて対話的、かつ容易にシナリオの編集が可能である。また、シナリオを修正する際、相互作用の強さに応じて、フレーム分子の定義を変更するべきかどうかを判断する情報が欲しい。そこで、重なり積分やクーロン相互作用に着目し、その表示機能も装備した(後述)。


Figure 8. Snapshot of ScenarioEditor. See also Figure 5.

3. 3. 3 大規模分子グラフィックス

各種表現に必要なタンパク質のための大規模分子グラフィックスを実装した。総原子数が120,792原子(重原子は58,884個)のシャペロンタンパク質はPentium4 3.06GHz + RADEON9700PROを使用して、全原子モデル構築に2秒で対応することができる。一度モデルが構築されれば、回転、拡大縮小などの操作はそれほど違和感のないスピードで行うことができ、十分な性能を達成している。その他にも、対話的な構造選択機能による選択部位複合表示、各種物理量表示、MD計算で得られた構造変化のアニメーション表示などをサポートしている。
Figure 9はインスリン(PDB ID: 1HLS)の等電子密度面に静電ポテンシャルの分布をマッピングした例である。一般にOpenGLを用いた表示では物体表面の材質情報や光源の設定により物体表面の任意の画素の光度を計算しているため、物体表面の光の反射により色分布を誤認する可能性がある。Figure 9では鏡面反射などの材質情報を変更し、静電ポテンシャルの分布の情報が正しく伝わるように表現している。このようにProteinEditorの分子グラフィックスは、高速にタンパク質の構造が表示できるだけでなく、結果の読み取りミスを軽減する表示の工夫を行っている。


Figure 9. Electrostatic potential mapped on isosurface of electron density.

3. 3. 4 機能の連携

ProteinEditorに実装されている各機能は、計算実行のミスの軽減や、より分かりやすい結果解析などを実現するために様々な箇所で連携している。Figure 10は自動計算法支援機能におけるScenarioEditor、MatrixViewer、分子グラフィックスの連動例である。ScenarioEditorは現在のシナリオを表示しているが、マウスを使ってこれらを対話的に編集できる。MatrixViewerはProteinDFが生成する行列のグラフィカル表示機能で、この例では重なり積分行列が表示されている。
(a)はScenarioEditorでマウスを用いて選択されたフレーム分子を短縮した結果である。点線で示した次のフレーム分子との重なっているフラグメント数を保存したままフレーム分子を伸縮できるようになっている。(b)はフレーム分子の編集に伴い表示が変わったMatrixViewerの表示結果である。一つのドットが一つのアミノ酸残基を示しており、他のアミノ酸残基との相互作用の強さを色合いの違いで表している。矩形の枠で示されているものが各フレーム分子である。編集前ではフレーム分子の切れ目が相互作用の高い部分に位置しているが、編集後では相互作用の低い部分が切れ目となるように調整されているのが分かる。(c)では選択されているフレーム分子がBall and Stick表示で示されており、選択されているフレーム分子がどこに位置しているかを確認することができる。
これら3つのウィンドウは連携して動作し、注目すべき残基の状態変化を様々な角度から観察することができる。これにより、より安全に計算シナリオを編集することが可能となっている。

3. 4 VR技術を取り入れたモデリング機能

ProteinEditorのタンパク質モデリング支援機能をベースに、さらに柔軟な処理に対応するようVirtual Reality(VR)技術を取り入れた高品位なモデリング支援機能の開発へ発展させている。


Figure 10. Cooperation of various functions in ProteinEditor. Top: Overview, (a) frame editing, (b) matrix viewing, and (c) molecular graphics.

3. 4. 1 タンパク質モデリングVRシステム

タンパク質立体構造の3次元空間における原子位置認識と原子間に生じる力の場を把握し、対話的な構造(原子、アミノ酸残基、基質、配位子)の編集を支援することを目的とした機能である。従来の分子モデリングでは、奥行き感があまり得られない2次元的な表示と平面的な動きしか表すことの出来ないマウスを用いるため、編集する構造位置決定が困難であることが多い。また数十原子以下の一般分子ではその必要性も低かった。数千原子からなるタンパク質の密構造を編集する場合には極めて有効であると考えている。
Figure 11に構成を示す。ProteinEditorのタンパク質構造のためのモデリング機能ModelingEditorに、反力デバイスPHANTOMによる触感インターフェース、および裸眼立体視ディスプレイによる視覚インターフェースを組み込んだ(SHARP製3D液晶を搭載したノートPCもしくは液晶ディスプレイ。SGI製立体視ライブラリISL対応ハードウェアでの立体視は開発中)。モデリングは反力デバイスを使用してタンパク質の構造を変異・修正させるが、その時に修正部位にかかる力をMDで計算し、デバイスにフィードバックさせる。立体視表示とこの触感の連動によって、直感的かつ精密に構造を修正することができる。なお、対話的なモデリングを実現するために、現在MDシミュレーション専用計算機MDGRAPE-3 [29]の組み込みも行っている。


Figure 11. Protein modeling VR system.

例えば、タンパク質のシミュレーションでは、アミノ酸置換がよく行われる。従来、このような操作で生じる構造歪みの除去にはMDによる構造最適化を行ってきたが、側鎖や官能基が自由に回転できない局所解に落ち込み、構造緩和が困難となる問題があった。結局これらは一つ一つ手作業で修正するケースが多い。本システムを使用すれば、より低いポテンシャルバリアを選んで強制的にこれを乗り越えさせる操作ができるため、このような作業の補助に適しているだろう。

3. 4. 2 タンパク質ポケットの認識

タンパク質(酵素)の表面や内部にはポケットと呼ばれるくぼみや空洞があり、作用を受ける基質分子はポケットに結合して適切な化学変化を受ける。タンパク質と基質の相互作用解析を行う上で、ポケットの認識は非常に重要である。しかし、ポケットの構造は複雑な上にタンパク質内部に隠れた状態になっていることがあり、粗視化による解析ではポケットと基質の結合状況、さらにはタンパク質立体構造との関係を正確に認識することが難しい。
本タンパク質モデリングVRシステムはこのような解析にも適している。さらに、没入型立体視システムを利用すれば、ポケット内部から全方向を見渡すことができるため、タンパク質活性部位の状況をより分かり易く認識できる。Figure 12はProteinEditorをCAVEシステムHoloStageTMに適用したタンパク質ポケット表示の実験風景である[30]。ポケット形状作成にはLIGSITEアルゴリズム[31]を用いた。紙面では分り難いが、試験者はポケット内部に視点中心を置き、立体メガネに付いた位置センサー(ヘッドトラッキングシステム)により、視点に同期して表示が移動・回転する。そのため、試験者はシステム内を自由に移動して構造の裏側もリアルタイムに確認できる。この実験によって、タンパク質のような巨大で複雑な系では、VR技術を用いた表示および状態の認識が非常に有効となることが明らかとなった。
ProteinEditorは量子化学計算を支援するGUIが出発点ではあるが、このように量子化学計算に直接御興味がなくとも、タンパク質構造を研究される方々に是非使っていただきたい。


Figure 12. Experiment of protein pocket display applying ProteinEditor to CAVE system (HoloStageTM).

4 おわりに

一般化学分子において量子化学計算は実用レベルで実験値を再現・予測する、なくてはならない手段となった。生体分子でもこの信頼の置ける方法が使えるようになりつつある。
Figure 13はProteinDFの各計算ルーチンが最初のSCF計算にかかった時間と軌道数をlog-logプロットしたものである[19]。分子積分、数値積分、対角化と行列積の傾きはそれぞれ2.1、1.8、2.9、2.9であり、SCF計算全体では2.3である。これらの値は計算がサイズの何乗に依存するかを表す。全電子カノニカル計算で倍精度によるシミュレーションの限界と目される100,000軌道計算では、行列演算がSCF計算の半分を占めるものの、このストレートな手法による計算達成はもはや時間の問題と言える。


Figure 13. Log-log plot between elapsed time of computational routines in the first SCF iteration and the number of orbitals. All calculations were carried out using 16 processors on a 1.3-GHz Itanium2 cluster (SGI Altix3700). The slopes of total calculation, molecular integrals, fitting routine of the XC terms, diagonalization, and multiplications are 2.3, 2.1, 1.8, 2.9, and 2.9, respectively.

当グループにおける、次のステップは化学反応を追跡するシミュレーションや、電子伝達や光励起反応を解析する励起状態計算の確立である。前者では薬剤やバイオ素材の機能と効率の解析が、後者ではスペクトルの計算から実験値との直接比較が行えるようになるため、ますます有用性が高まるだろう。もちろん、これらのためのGUIも適宜開発中である。実験と理論が両輪として働くようになったとき、その研究分野は飛躍的に発展する。計算機の弛まぬ発展も強力に本分野をサポートしている。

本稿は筆者が代表してまとめたものであるが、紹介した内容は開発に携わったチーム全員の研究成果である。本稿で紹介した成果は、文部科学省次世代IT基盤構築のための研究開発「革新的シミュレーションソフトウェアの研究開発」プロジェクト、および東京大学生産技術研究所の「選定研究」の支援を受けたものである。計算には統計数理研究所の統計科学スーパーコンピュータシステムを利用した(2006-ISM.CRP-1016)。

参考文献

[ 1] F. Sato, Y. Shigemitsu, I. Okazaki, S. Yahiro, M. Fukue, S. Kozuru and H. Kashiwagi, Int. J. Quant. Chem., 63, 245 (1997).
[ 2] F. Sato, T. Yoshihiro, M. Era and H. Kashiwagi, Chem. Phys. Lett., 341, 645 (2001).
Chem. Phys. Lett., 392, 565 (2004).
[ 3] T. Yoshihiro, F. Sato and H. Kashiwagi, Chem. Phys. Lett., 346, 313 (2001).
[ 4] 柏木浩, 佐藤文俊, 吉廣保, 稲葉亨, 西川宜孝, タンパク質量子化学計算 ProteinDFの夢と現実, 柏木浩, 佐藤文俊編, アドバンスソフト, 東京 (2004).
[ 5] M. Challacombe and E. Schwegler, J. Chem. Phys., 106, 5526 (1997).
[ 6] K. Tsuda, H. Kaneko, J. Shimada and T. Takada, Comp. Phys. Commun., 142, 140 (2001).
[ 7] F. Maseras and K. Morokuma, J. Comp. Chem., 16, 1170 (1995).
[ 8] K. Kitaura, T. Sawai, T. Asada, T. Nakano and M. Uebayasi, Chem. Phys. Lett., 312, 319 (1999).
[ 9] 文部科学省ITプログラム「戦略的基盤ソフトウェアの開発」,
URL:
http://www.fsis.iis.u-tokyo.ac.jp/
[10] 次世代IT基盤構築のための研究開発「革新シミュレーションソフトウェアの研究開発」,
URL:
http://www.rss21.iis.u-tokyo.ac.jp/
[11] H. Sambe and R. H. Felton, J. Chem. Phys., 62, 1122 (1975).
[12] J. Andzelm and E. Wimmer, J. Chem. Phys., 96, 1280 (1992).
[13] B. I. Dunlap, J. W. D. Connolly and J. R. Sabin, J. Chem. Phys., 71, 3396 (1979).
[14] S. F. Boys, Proc. Roy. Soc. [London], A, 200, 542 (1950).
[15] S. Obara and A. Saika, J. Chem. Phys., 84, 3963 (1986).
[16] K. Eichkorn, O. Treutler, H. Ohm, M. Haser and R. Ahlrichs, Chem. Phys. Lett., 242, 652 (1995).
[17] F. Weigend, Phys. Chem. Chem. Phys., 4, 4285 (2002).
[18] A. D. Becke, J. Chem. Phys., 88, 2547 (1988).
[19] L. S. Blackford, J. Choi, A. Cleary, E. D'Azevedo, J. Demmel, L. Dhillon, J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. Walker and R. C. Whaley, ScaLAPACK User's Guide, Philadelphia, PA (1997).
[20] T. Inaba and F. Sato, J. Comput. Chem., 28, 984 (2007).
[21] H. Kashiwagi, H. Iwai, K. Tokieda, M. Era, T. Sumita, T. Yoshihiro and F. Sato, Mol. Phys., 101, 81 (2003).
[22] T. Inaba, S. Tahara, N. Nishikawa, H. Kashiwagi and F. Sato, J. Comput. Chem., 26, 987 (2005).
[23] D. G. Anderson, J. ACM., 12, 547 (1965).
[24] T. Inaba, N. Tsunekawa, T. Hirano, T. Yoshihiro, H. Kashiwagi and F. Sato, Chem. Phys. Lett., 434, 331 (2007).
[25] 西村民男, 稲葉亨, 小池聡, 平野敏行, 田原才静, 吉廣保, 西村康幸, 西野典子, 佐藤文俊, 柏木浩, ProteinDFによるタンパク質量子化学計算事例集, アドバンスソフト, 東京 (2005).
[26] 西村康幸, 吉廣保, 佐藤文俊, 生産研究, 59, 93 (2007).
[27] The Protein Data Bank,
URL:
http://www.rcsb.org/pdb/
[28] Application Visualization System,
URL:
http://www.kgt.co.jp
[29] M. Taiji, Proceedings of Hot Chips 16, IEEE Computer Society (2004), in CD-ROM.
[30] 西村康幸, 吉廣保, 佐藤文俊, 第23回CG・可視化研究会(CAVE研究会) (2006).
[31] B. Huang and M. Schroeder, BMC Struct. Biol., 6, 19 (2006).


Return