BioServerによる高精度結合自由エネルギー計算

全文

(1)2005−BIO−1（6） 2005／7／25. 社団法人情報処理学会研究報告 IPSJ SIG Technical Report. BioServer による高精度結合自由エネルギー計算藤谷秀章富士通研究所概要富士通が開発した 8 並列 VLIW アーキテクチャの高性能 FR-V プロセッサーを 1 ラックに 1920 台塔載した BioServer 超並列コンピュータを用いて、FKBP 蛋白質とリガンドの結合自由エネルギーの絶対値を計算した。AMBER 力場とベネット受容比法を用いる事により、創薬応用に適用可能な計算精度で、結合自由エネルギーの実験値と良い一致をみた。これにより BioServer の様な超並列計算機が、薬化合物を合理的に設計する為に有用な事を示した。. The absolute binding free energy for FKBP protein and ligands were calculated using the Fujitsu BioServer massively parallel computer, which has 1920 FR-V processors (high performance and low power processor with 8 way VLIW architecture) in a rack. Using the AMBER force field and the Bennett acceptance ratio we obtained pharmaceutically useful predictability between the calculated and experimental absolute binding energies. The ability to perform these calculations demonstrates the utility of massively parallel computation resources such as the BioServer in rational drug design.. 1. 生体高分子のシミュレーション. 生体機能はそれを構成する種々の生体高分子が有機的に機能することで実現されている。個々の生体高分子の機能を理解する事は生体機能を理解する出発点になるが、例えば最も小さい蛋白質の一つである 36 個のアミノ酸基からなる villin でも、その総原子数は 389 個にもなり、villin の性質や挙動を計算化学の手法を用いて解析する事は極めて難しい問題になる。 1997 年に Kollman が超大䍫䍶䍍䍖䍛 ಽሶౝ ᳓ಽሶ䈱䍞䍻䍨䍽䍖⾰䈱䍪䍉䍎䍷䍡䍼䍆䍻䍖䍼型スーパーコンピュータを ↢ᚑ ᝄേ ㆇേ 数か月間独占して、villin の折畳み現象を分子動力学の 10-15 10-12 10-9 10-6 10-3 100 方法で 1 マイクロ秒までシ femto⑽. pico⑽. nano⑽. micro⑽. milli⑽. ⑽ ミュレートしたが、原子レベルの時間スケールがフェ図 1 時間スケール −15 ムト (10 ) 秒に対して villin が実際に折畳まれる時間が 10 マイクロ (10−5) 秒ほどであるために、実際に畳み込まれる時間までのシミュレーションには到達出来なかった [1]。この事もあり、蛋白質分子の畳み込み計算は現代計算科学の超難問の一つとして位置付けられている。 2002 年に Pande は世界中のパーソナルコンピュータにシミュレーションプログラムをダウンロードしてもらう Distributed Computing の方法を用いて、villin の畳み込み計算を実行して実験値と一致する畳込みレートを求める事に世界で初めて成功した。Pande はスーパーコ. −35− 1.

(2) ンピュータなどで用いられている従来の並列アルゴリズムではなく、世界中のインターネットで繋がっている数万台のコンピュータにシミュレーションを並列分散させる独自の方法を考案して超難問の蛋白質の畳み込み計算を初めて実現した [2]。スーパーコンピュータなどで用いられている Message Passing Interface (MPI) などの並列化アルゴリズムでは CPU 台数が多くなるに従い CPU 間のデータ通信が占める割合が飛躍的に増加して効率的な計算を妨げる事になる。例えば数万台の CPU で数万倍の効率的なシミュレーションを実行するためには、シミュレーションの根本アルゴリズムから変更する必要がある。蛋白質がリボゾームで生成され、畳み込みにより立体構造を獲得した後に、他の生体高分子やリガンドとどの様に相互作用をするかを調べる事は、基礎的研究の大きな課題であると同時に応用の観点からも極めて重要である。特に蛋白質と薬候補化合物であるリガンドとの結合自由エネルギーを正確に求める事は、創薬分野における永年の課題である。我々はこの課題に対して独自のアルゴリズムと超並列専用計算機を開発して、FKBP 蛋白質とリガンドの結合エネルギーを正確に求める事に成功した [3]。. 2. 結合自由エネルギー計算方法. あるリガンドがターゲットとする蛋白質に強く結合する場合は、その蛋白質本来の機能をリガンドが阻害して薬効果が現れる可能性が高い。この為、種々のリガンドが蛋白質に結合するエネルギーを精確に求める事が、計算機を使った薬設計 (Computational Drug Design) を実現する為の必須条件になる。一般に薬候補化合物と蛋白質の結合エネルギーは -10kcal/mol 程度であり、薬候補化合物間の優劣を判断するには 1kcal/mol ほどの精度で結合自由エネルギーを求める必要がある。室温での熱揺らぎエネルギー (kB T kB :ボルツマン定数) が約 0.6kcal/mol なので、1kcal/mol の精度で結合エネルギーを計算するのは至難の業であり、従来の方法では実現不可能であった。今迄に蛋白質とリガンドの結合エネルギーを求める様々の計算方法が提案されているが、その多くが実験データに基いた経験的な方法であり、実験データのない未知の系についても適用可能になる為に必要な “強固な物理理論に基づいた方法” は行われていなかった。そもそも結合エネルギーとは蛋白質とリガンド分子が水中 (細胞液) で結合した状態と、蛋白質とリガンドがそれぞれ別に水中に溶けている二つの異った平衡状態の間の自由エネルギー差である。この異った二つの熱力学状態の間を遷移する時に、蛋白質やリガンドは熱揺らぎによって様々に構造が変化したり、弱く結合する状態から自由エネルギーがより低くて強く結合する状態へと結合状態を変化させる。一般に薬化合物は数十から二百くらいまでの原子からなる分子であり、蛋白質は数百から数十万原子からなる巨大分子なので、結合エネルギーを精確に求めるには、これらの構造変化によるエントロピー変化を精密に計算する必要がある。蛋白質とリガンドが離れた状態から安定な結合状態に遷移するのにマイクロ秒のオーダを要するが、この間のエントロピー変化を計算するにはこの間の分子動力学シミュレーションを行って経巡る相空間のサンプリングを行う必要がある。分子動力学の 1 ステップを 1 秒で計算する高速のコンピュータを用いても、1 マイクロ秒までの 109 ステップを計算するには 32 年を要する。この不可能な計算を実行するには、まったく新しい超並列化アルゴリズムが必要である。 1997 年に Jarzynski は熱平衡状態間の非平衡仕事量の指数関数平均が二つの状態間の自由エネルギー差になるという基礎的な等式 ΔG = −kB T lnexp(−W/kB T )av を証明した [4]。このアイデアに基いて蛋白質とリガンドが別々に水に溶けて蛋白質とリガンド間に相互作用が −36− 2.

(3) 無い状態 H[λ = 0] と、蛋白質とリガンド間の相互作用が働いて結合している状態 H[λ = 1] を考える。相互作用の強さ λ をパラメータとして種々の中間状態を考えて、その隣合う状態間の遷移に必要な仕事量を求める。この仕事量の単純な指数関数平均を取るだけでは、仕事分布の裾の誤差が大きく現れて、精確な自由エネルギー差は得られない。これを克服する為に、図 2 の様に λ が小くなる方向と大きくなる二方向の仕事量を独立に求めて、その二つの分布からベネット受容比法 (Bennet Acceptance Ratio: BAR) を用いてそれぞれの間の自由エネルギー差を求める [5]。. 㱗䋽1. 㱗䋽0.75. P F (W ). P R( W). G1. P F (W ). +. 㱗䋽0.5. P R( W). G2. P F (W ). +. 㱗䋽0.25. P R( W). G3. P F (W ). +. 㱗䋽0. P R( W). G4. =G. 図 2 非平衡仕事量の測定図 2 では五つの λ 状態を図示しているが、計算精度を向上させる為に λ の刻み間隔をより小くして最適化する必要がある。またそれぞれの λ で独立に分子動力学を実行するが、その初期運動量分布によってサンプリングされる相空間が異なり、異なる自由エネルギーが得られるので、異なる初期運動量分布を持った λ のセットを多数用意する必要がある。λ の間隔と異なる初期運動量セットの数を最適化した処、33 点の λ 状態を 60 セット用いることで最もよい自由エネルギー分布が得られ、これ以上初期運動量セットを増やしても同じ分布が得られる事が分った。これで 33 × 60 = 1980 の分子動力学を同時並列に実行して、個々の分子動力学計算を 1 ナノ秒実行すると、全体ではマイクロ秒オーダの探索を行ったのと同等の計算量になる。. 3. FKBP リガンド. FK506 は臓器移植やアトピーなどの病気治療に用いられる免疫抑制剤で、FK506 に結合する蛋白質 FKBP(FK506 Biding Protein) は殆どの生物種に存在して蛋白質の主鎖の向きを大きく変換出来る酵素活性を持っている。FKBP は 107 個のアミノ酸基からなる蛋白質で、 1666 個の原子からなる。図 3 は今回結合エネルギー計算を行った FK506 とそのアナログである。図の右側の数値は実験で得られた抑制定数 (inhibition constatant: Ki ) であり、kB T ln Ki が結合エネルギーの実験値になる [6]。 FKBP 蛋白質とこれらのリガンドの結合自由エネルギーの絶対値を計算で求めるには、まずその結合している時の原子座標を全て知る必要がある。L8 、L9 、FK506 の三種類のリガンドについては FKBP と結合している状態の X 線解析データがあり、原子座標値が Protein Data Bank に登録されている。他のリガンドについては結合状態の X 線データが無く、構造が分っている三つのデータから結合状態を推測して初期座標を决めた。この FKBP とリガンドが結合している回りに水分子をランダムに配して分子動力学計算に用いる初期構造を作 −37− 3.

(4) 成した。. L2. L8. 10. L9. 7. 2000. L3. 660. L5. 110. L6. 12. L12. 30. FK506. 0.4. 図 3 FKBP リガンド全ての系に対して初期構造を緩和する為に室温での分子動力学計算を 10 から 20 ナノ秒ほど行った。水分子が初期のランダム配置から、室温で蛋白質の回りに構造をもった状態に緩. (a). (b). 図 4 FKBP と L2 リガンドの結合状態和するのに通常 1 から 3 ナノ秒掛るが、これに比較すると遥かに長い時間の緩和計算が必要 −38− 4.

(5) である。理由の一つには、初期構造が安定ではなく緩和シミュレーションを行っている間に構造が変化する事による。図 3 の全てのリガンドはピペコリトと α-ケトアミノを持っており、これが L8 、L9 、FK506 では FKBP の疎水ポケットの底に入る様に結合している。その為、他のリガンドも同様の結合形状を持っていると推測した。図 4(a) は L2 の初期構造でピペコリトが FKBP の疎水ポケットの底にあり、黄色で示されている FKBP の 59 番目アミノ酸であるトリプトファンのインドールの直ぐ上にある。これが分子動力学シミュレーションを行っている間に、疎水ポケットの中で回転して (b) の様にピペコリトが上になって安定した。この間約 14 ナノ秒である。L3 でも同様にピペコリトが疎水ポケットの上になったが、他のリガンドではこの様な回転は起らずに初期構造と同じ疎水ポケットの底にピペコリトは滞った。. 図 5 BioServer による結合エネルギー計算この緩和された構造を超並列結合エネルギー計算の初期構造として BioServer を使用して実行する。BioServer は一ラックに 1920 個の FR-V プロセッサを塔載した計算エンジンで [7] 、磁気ディスクを持ったサーバー機とペアになって機能する。実際に計算を進めるには、サーバ機で相互作用強度 λ や初期運動量分布を変えて初期の入力データを作成して、各 FR-V プロセッサーで 100 ピコ秒の分子動力学シミュレーションを実行させる。その結果をサーバ計算機に戻して、次の 100 ピコ秒のシミュレーションを行う為の入力データを作って FR ー V で次の計算を行わせる。この様にしてサーバ機には全ての独立な分子動力学シミュレーションの 100 ピコ秒毎の結果が蓄えられる。この 100 ピコ秒毎の二方向の仕事量から BAR 法で結合エネルギーを計算し、数千の FR-V プロセッサーで同時進行中の計算の進み具合をモニターする。 BioServer の中では相互作用の強さ λ が異なるシミュレーションが独立に行われているので、ある FR-V プロセッサを見ると蛋白質とリガンドが密に結合した状態の相空間のサンプリングが起こなわれており、他の FR-V プロセッサでは蛋白質とリガンドが離れた状態が、また他では中間の強さの相空間のサンプリングが行われている (図 5)。この方法による結合エネルギーの収束速度はリガンドの種類にも依存するが、大体 1000 個オーダの FR-V プロセッサを使用して 0.5 から 1.5 ナノ秒の分子動力学によるサンプリングによって良好な収束 −39− 5.

(6) が得られる。この様にして蛋白質とリガンドが水中で離れた状態から安定な結合状態を形成するまでのマイクロ秒オーダの相空間サンプリングを、1000 個オーダの並列化を行って同時処理する事が可能になった。. 4. 計算精度. BioServer の中の個々の FR-V プロセッサで行っている分子動力学シミュレーションには単精度版の改造した GROMACS(http://www.gromacs.org) を使用した。有限要素法を使った. Calculation (kcal/mol). -3 構造解析や電気回路シミュ L2 レーションなど殆どの数値 -4 L3 シミュレーションには倍精 -5 度演算が必須であるが、分 L5 -6 子動力学シミュレーション L12 は例外的に単精度で殆どの -7 L8 演算を行っても計算精度に -8 L6 は支障がない。これは運動 L9 -9 方程式を解く時間ステップが短い為に、原子間に働く -10 FK506 力の計算精度を上げる必要 -11 がない為である。 0.1 1 10 100 1000 10000 図 6 は BioServer によって Experimental Inhibition constant (nM) 得られた FKBP 蛋白質と 8 図 6 FKBP の結合自由エネルギー種類のリガンドの結合自由エネルギーを縦軸に、抑制定数の実験値を横軸に示した。計算値の直線からのずれの RMS エラーは 0.4 kcal/mol であり、従来の精度と比較すると一桁以上の改善が見られる。倍精度の GROMACS を用いて同じ計算を行ったが、単精度の計算値からのずれは 0.2 kcal/mol 以内で、結合自由エネルギー計算には単精度版 GROMACS による計算で十分である事が実証された。. 5. 計算機の新たな方向. 蛋白質とリガンドの結合自由エネルギーの絶対値を 1 kcal/mol 以内の精度で計算する事はコンピュータによる創薬を行う上で永年の悲願であった。BioServer を使用した超並列計算による結合エネルギー計算方法は、AMBER 力場を用いておりその基礎は分子の電子構造に基づく。また解離状態から結合状態に致るエントロピーを正確に計算する方法は物理の基礎理論に基いている。この様に “強固な物理理論に基づいた方法” を用いたシミュレーション技術は実験結果に新たな光明を当てて、より深く事象を理解する大きな助けになる。一個の蛋白質とリガンドの結合エネルギーを求めるのに膨大な計算機資源を必要とするが、得られる知見は何者にも代えがたい深い意味を含んでいる。この方法が創薬の現場にどの程度のインパクトを与えるかは今後の進展を見守る必要があるが、まずこの様な従来は不可能であったシミュレーションを実行可能にする計算機とソフトウェアを提供する事が重要である。図 7 は BioServer の製品版プロセッサボックスである。 −40− 6.

(7) 一枚のボードに 8 並列 VLIW アーキテクチャの高性能 FR-V プロセッサーが 4 個載り、それが 2U のボックスの中に 32 枚実装されている。ある蛋白質とリガンドの一つのセットの自由エネルギーを求めるのに少なく見積ってもプロセッサボックス 5 個くらいを使うのが効率的である。一つのターゲット蛋白質に対して複数のリガンドの結合エネルギーを平行して求めるのが一般的なので、実際に自由エネルギー計算をスタートすると計算機を休み無く動かす様になり、いとも簡単に 1920 個のプロセッサが入ったラックを自由エネルギー計算のジョブで一杯にして仕舞う。動くのは単精度 GROMACS だけである。. Processor Box (PB) Number of PE. 128PE. Peak performance. 184GFLOPS. Memory capacity. 32GB. Power consumption. 720VA. LAN. 1000Base-T㬍2port 100Base-TX㬍1port. Size (W㬍D㬍H). 482mm㬍785mm㬍87mm (2U). 図 7 AC55: 製品版 BioServer プロセッサボックス. SPARC や X86 などの高速汎用プロセッサでは、スーパースカラーなどのハードウェアによる並列化処理がサポートされており、多くのプログラムを効率的に動作させるのが比較的容易であるが、純粋な演算以外の処理が多くて電力の消費が大きくなる。これに対して VLIW アーキの FR-V ではプログラム構造とコンパイラ技術だけで並列化効率が決るので、任意のプログラムを効率的に動作させるのは難しくなる。しかし FR-V では純粋な並列演算処理だけが動いて、前処理などによる余分な電力消費は無く、同一演算処理量にたいする消費電力は FR-V の方が格段に少くなる。高性能 X86 プロセッサでは CPU 一個だけで 100 ワット以上を消費するのに対して、FR-V では 1.5 ワット以下である。プロセッサのクロック周波数比は 10 倍程であるのに対して、1 個のプロセッサ当りの消費電力は 60 倍近く異なる。この為に FR-V の VLIW 並列度が活用出来るアプリケーションを動かす場合には、消費電力や設置面積で優位になる。今回の場合の様に一つのプログラムが沢山のプロセッサで四六時中動いている状況では、VLIW に対す −41− 7.

(8) る並列最適化に多少の労力を割いても最終的パフォーマンスで多大なメリットが得られる。地球温暖化が大きな社会問題になっている現代では、低消費電力並列計算機が活躍する場面が今後多く現れると考えられる。. 謝辞共同研究者である Stanford 大学の Guha Jayachandran 、Christopher D. Snow 、Michael R. Shirts 、Eric J. Sorin 、Vijay S. Pande 、富士通研究所の谷田義明、伊藤正勝の各氏と、富士通 BioServer プロジェクトのメンバーに謝意を表します。GROMACS の改造では Stockholm 大学の Erik Lindahl 氏にお世話になりました。この研究は NEDO の『バイオ・IT 融合機器開発プロジェクト』の支援を受けています。. 参考文献 [1] Y. Duan, L. WAng, and P. A. Kollman, Proc. Natl. Acad. Sci. USA 95, 9897 (1998). [2] C. D. Snow, H. Nguyen, V. S. Pande, and M. Gruebele, Nature 420, 102 (2002). [3] H. Fujitani, Y. Tanida, M. Ito, G. Jayachandran, C. D. Snow, M. R. Shirts, E. J. Sorin, and V. S. Pande, J. Chem. Phys. in press. [4] C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997). [5] M. R. Shirts, E. Bair, G. Hooker, and V. S. Pande, Phys. Rev. Lett. 91, 140601 (2003). [6] D. A. Holt, J. I. Luengo, D. S. Yamashita, H. J. Oh, A. L. Konialian, H. K. Yen, L. W. Rozamus, M. Brandt, M. J. Bossard, M. A. Levy, D. S. Eggleston, J. Liang, L. W. Schultz, T. J. Stout, and J. Clardyt, J. Am. Chem. Soc. 115, 9925 (1993). [7] H. Okano, A. Suga, T. Shiota, Y. Takebe, Y. Nakamura, N. Higaki, H. Kimura, H. Miyake, T. Satoh, K. Kawasaki, R. Sasagawa, W. Shibamoto, M. Sasaki, N. Ando, T. Yamana, I. Fukushi, S. Tago, F. Hayakawa, T. Kamigata, S. Imai, A. Satoh, Y. Hatta, N. Nishimura, Y. Asada, T. Satoh, T. Sukemura, S. Ando, and H. Takahashi, ISSCC Dig. Tech. Papers, 374 (2002).. −42− 8.

(9)