JP4314128B2 - タンパク質立体構造と誘導適合を利用したリガンド探索方法 - Google Patents

Description

本発明はタンパク質の立体構造座標を用いたリガンド探索方法、詳しくはタンパク質立体構造座標が既知の場合、相互作用すると考えられるリガンドを予測する方法に関する。
更に、本発明は、この方法で得られるリガンド、前期方法に使用可能なデータベース、データベース構造、及びコンピューターソフトプログラム、これを搭載したコンピュータやインターフェース等にも関する。

酵素や受容体等の生体機能を維持するために必要なタンパク質には、基質特異性と呼ばれる性質があり、活性部位が基質分子構造の細部にわたり常に一致しているＬｏｃｋ＆Ｋｅｙ型と、基質が無いときには活性部位が不活性なランダムな状態にあり、基質が来るとこれを取り込むために活性部位が活性な状態に変化するＩｎｄｕｃｅｄ−Ｆｉｔ（誘導結合）型がある。誘導適合型とは、リガンドと結合する際にリガンド結合部位の立体構造が変化しリガンドを取り込むことが可能になる受容体をいう。

タンパク質の立体構造を用いたリガンド分子探索のための計算化学的手法としてはまずＤＯＣＫ、ＦｌｅｘＸ、Ｌｕｄｉ、ＧＯＬＤ、といった３次元化合物データベースサーチ（Ｖｉｒｔｕａｌ Ｓｃｒｅｅｎｉｎｇ）が知られている。これらは高速ドッキングスタディーとも呼ばれ大規模な化合物ライブラリサーチが可能である。しかし、本手法では評価に粗い近似を用いるため、結合配座や結合エネルギーの予測能は低い。さらにタンパク質とリガンドとの結合に大変重要な「誘導結合」に対応する計算式パラメータを充分に取り込んでいないので、たとえあったとしても、乱数を発生させ受容体の側鎖を動かす程度であり、計算結果の精度に充分なものとはいえない。

タンパク質とリガンドとの結合に重要な「誘導結合」をシミュレーションする方法としてはＭＤ（分子動力学計算）、ＭＭ（分子力学計算）ＭＣ（モンテカルロ法）が知られている。これらの方法は比較的精度良く、結合配座や結合エネルギーの予測が可能である。ここで、分子動力学法（ＭＤ）と呼ばれる手法に関しては、ある分子を構成する各原子において、古典力学に基づく運動方程式を逐次的に解くことにより、その分子の動的構造を計算する方法であり、タンパク質の動的挙動を高精度でシュミレーションすることが可能である。しかし、計算に時間を要するため、多数の分子を扱うことは困難であり必ずしも有用な手法とはなっていない。さらに、従来法では該当タンパク質に対して分子動力学計算を行うとタンパク質立体構造はＸ線、ＮＭＲ等で解析された座標から大きくズレる。こうしたズレはタンパク質の動的挙動の物理化学的描写を含んでいるがＮＭＲ等で示される動的挙動の実験的な結果と矛盾する挙動となる場合があり、必ずしも精度の高いシミュレーションとならないことが多い。

このように従来のｉｎ ｓｉｌｉｃｏ ｓｃｒｅｅｎｉｎｇ関連では、タンパク質とリガンドとの結合に大変重要な「誘導結合」に対応する計算式パラメータを充分に取り込んでいないので、計算結果の精度に充分なものとはいえない。

一方、分子シミュレーションでは上記の誘導結合を表現し、解析することは可能であるが、高精度の結果を得るためには相当の時間を必要とする。多くの結果は、初期構造座標に依存してしまう。

本発明者等は、任意のタンパク質の立体構造が与えられたとき、該当タンパク質に結合するリガンドを探索する方法について検討をおこなった。現在流通している受容体・リガンド結合解析ソフトには、リガンドのフレキシビリティを考慮しているものは多くあるが、受容体側のフレキシビリティを考慮しているものはほとんどない。たとえあったとしても、乱数を発生させ受容体の側鎖を動かす程度であり、Ｌｏｃｋ＆Ｋｅｙ型の受容体に対応しているものばかりであった。そこで、Ｉｎｄｕｃｅｄ−Ｆｉｔ型の受容体を対象にした受容体・リガンド結合解析ソフト開発することにした。

本発明が解決しようとする課題は、農薬、医薬品等の開発に特に重要な鍵となる、該当タンパク質に結合するリガンドを探索する方法を精度よく、かつ従来法に比べてはるかに効率的に精度よく探索する方法を提供することである。また、リガンド分子の多様な改変や受容体等のタンパク質の改変を迅速かつ効率的に行う方法を提供することにもある。更に本発明により、リガンド−タンパク質間の相互作用様式を解明し、それ等相互作用の認識機構を明確化することで、疾病の原因を特定したり、それ等に関連する薬物の開発を促進すること等を目的とする。

本発明者等は、任意のタンパク質立体構造が与えられたとき該当タンパク質に結合するリガンドを探索する方法について検討を重ねた結果、下記［１］〜［６］の方法およびそのためのコンピュータプログラムを見出し、或いは開発した。

ここで、分子動力学法（ＭＤ）と呼ばれる手法があり、これはある分子を構成する各原子において、古典力学に基づく運動方程式を逐次的に解くことにより、その分子の動的構造を計算する方法である。つまり、これはある分子を構成する各原子における古典力学を土台とした動的挙動を計算する方法である。従って、この手法をうまく取り込むことができれば、リガンドを取り込んでいない状態のＩｎｄｕｃｅｄ−Ｆｉｔ型受容体を初期状態に選んでも受容体・リガンド結合を再現できると考えた。ＭＤ計算は古典力学を土台にしているため、各原子にある程度の拘束をかける必要がある。そこで、まず初めに受容体の基準振動解析を行い受容体の主鎖二面角揺らぎを計算し、この主鎖二面角揺らぎに基づいて各原子に拘束をかけてＭＤを計算する手法を開発した。具体的には、基準振動解析計算をまず行い、定常状態の主鎖二面角の揺らぎを計算する。そして、その揺らぎを基にした拘束を各原子にかけながら分子動力学計算を行うことでより精度の良い受容体の動的構造を予測する。また、分子動力学計算より得た動的構造及び相互作用関数を用いることで、誘導適合型受容体にも応用できる受容体／リガンド結合を精度良く予測できる。本発明の方法は、より真に近い受容体／リガンド結合を予測する方法である。従って、本発明の方法は、医農薬分子の設計に極めて有用である。

［１］ 任意の単数を含む複数鎖のタンパク質立体構造が与えられた場合において、該当タンパク質の立体構造から誘導適合を反映したパラメータおよび構造変化した立体構造座標を例えば基準振動計算方法や分子動力学計算方法よりあらかじめ算出し、当該パラメータおよび構造変化した立体構造座標を用いて該当タンパク質と別の物質が結合した場合の相互作用関数を定義し、当該相互作用関数によって該当タンパク質と結合する物質をコンピュータプログラムにより評価し、選定する方法。

［２］ ［１］によって記載された方法において、該当タンパク質に結合するリガンドを選択する際に（０）〜（８）に示した一連の処理を全自動または手動的に行うことを特徴とするコンピュータプログラム。

（０） 化合物データベースからリガンドを１つ選択する。該当タンパク立体構造として、誘導適合を反映するパラメータを用いて動的挙動を考慮した複数の構造変化座標を用意し、ランダムに１つの構造を選択する。

（１） 重ね合わせを行う該当タンパク中の空間点を指定する。空間点は例えば以下のような方法で発生させる。

１． ダミー原子の発生による空間点の発生
リガンドとタンパク質との相互作用における水素結合に着目し、タンパク質中の水素結合サイトを空間点として指定する。水素結合における重要事項は距離と角度である。つまり、角度を計算するためには、水素結合ドナー（以後、ドナー）に水素原子が必要になる。

そこで、活性部位及びリガンドに水素原子が含まれていない場合、以下の規則によりダミー水素原子を発生させた。

１） ｓｐ²軌道原子を中心とする正三角形状にダミー原子を発生（図２）。すなわち、図２に示すように、ｓｐ²軌道原子の窒素原子（Ａ）を中心とする正三角形の空いている位置にダミー水素原子（Ｂ）を発生させた。

２） ｓｐ³軌道原子では、水素結合を形成する距離にある場合、水素原子を共有するように回転できると考え水素結合相互作用を計算するときには、距離のみを考慮することにした。このため、ｓｐ³軌道原子にはダミー原子を発生させない。

金属及び水では、活性部位・リガンド結合の仲介役となり得るので、相互作用する位置に以下のようにダミー原子を発生させた。

１） 鉄のような金属には、正八面体状にダミー原子を発生（図３）。すなわち、図３に示すように、亜鉛（Ａ）を中心とする正八面体の空いている位置にダミー原子（Ｂ）を発生させた。

２） 水は、正四面体状にダミー原子を発生。

ただし、活性部位と相互作用している方向にはダミー原子を発生させないことにした。

２． 構造活性相関情報を利用した空間点の発生
リガンドの構造活性相関（ＳＡＲ）情報に着目し、以下の項目を入力情報にすることにした。

（Ａ） ＳＡＲから得られた活性部位の原子（以後、「Ａ原子」）。ＰＤＢ形式に従う。

（Ｂ） 「Ａ原子」と相互作用するであろうリガンドの原子タイプ（以後、「Ｂタイプ」）。ＳＹＢＹＬのＭＯＬ２形式に従う。

（Ｃ） 「Ａ原子」と「Ｂタイプ」との相互作用の強さ（以後、「Ｃ強さ」）。

（Ｄ） 「Ａ原子」と「Ｂタイプ」との相互作用する距離（以後、「Ｄ距離」）（単位はÅ）。

Ａ）〜Ｄ）をもとにタンパク質中の活性部位内におけるリガンドの初期座標を利用し、以下の規則により空間点を作成することにした。

１） 「Ａ原子」がドナーまたは金属及び水の場合（ＳＡＲ情報の活性部位側の指定が水素結合ドナー、金属原子の場合）、１．で発生させたダミー原子の方向に対して「Ａ原子」から「Ｄ距離」の位置及びその周囲を初期座標に選んだ（図４、図５）。

２） 「Ａ原子」がｓｐ³軌道原子の場合（ＳＡＲ情報の活性部位側の指定がｓｐ³軌道原子の場合）には、「Ａ原子」から「Ｄ距離」の周囲を初期座標に選んだ（図６）。

３） 「Ａ原子」が水素結合アクセプター（以後、アクセプター）の場合（ＳＡＲ情報の活性部位側の指定が水素結合アクセプターの場合）、「Ａ原子」の結合延長上「Ｄ距離」の位置及びその周囲を初期座標に選んだ（図７）。

４） その他の場合（ＳＡＲ情報の活性部位側の指定がその他の原子の場合）には、「Ａ原子」を中心とする半径が「Ｄ距離」の球表面上の点を初期座標に選んだ（図８）。

５） （１）〜（４）とは、異なりリガンドの初期座標を直接指定することもできるようにした。

（２） （０）で選択したリガンド中の原子と（１）で指定した空間点とのペアを重複がないようにランダムに選択する。

（３） 以下のスコアＳｓｃｏｒｅ（ｉ，ｊ）を計算する。

ｄ_ij ^sは該当タンパク質中のｉ番目とｊ番目の空間点距離。ｄ_ij ^Cは化合物中のｉ番目とｊ番目の原子間距離。αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数。βは重なりと定義できる限界値を与えるための定数。

αは１．５、βは０．８とするのが好ましい。

（４） （３）のスコアが最大になるように調整する。スコアを最大にする手法としては、例えば、シミュレーティッドアニーリング法が挙げられる。または、時間短縮には（２）、（３）を１００００回繰り返し、Ｓｓｃｏｒｅ（ｉ，ｊ）が最大になるペアを探し、そのペア情報をもとにリガンドを初期座標に重ね合わせる方法を適応することが好ましい。

（５） （４）で重ねあわせたリガンドに対してタンパク質との相互作用エネルギーをコンフォメーションを微調整しながら最適化計算する。リガンドのコンフォメーションの微調整は、（４）で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３を越えない程度に座標変化させる。

微調整は例えばランダムサーチで最適化することが好ましい。ランダムサーチは以下の項目に従ってタンパク質の活性部位とリガンドとの微小変化を８０００回行い、最適エネルギー「Ｕ最適」が最小になるようにする。

１〕 回転可能な結合のうち最大５つ乱数で選び、結合ごとにランダムに±１０．０°の範囲内で回転させリガンドのコンフォメーションを換える。この過程を３回に一度行う。

２〕 ｘ、ｙ、ｚ軸方向それぞれにおいて、ランダムに±１．０Åの範囲内でリガンドの並進運動を行う。この過程を２回に一度行う。

３〕 回転中心座標それぞれにおいて、ランダムに±１．０Åの範囲内で回転中心座標を移動させ、さらに３次元方向の角度それぞれに対して、ランダムに±５．０°の範囲内でリガンドの回転運動を行う。この過程を５回に一度行う。

（６） リガンドのコンフォメーションを大きく動かして、（２）から再スタートを行い、（５）までを繰り返して最適化を行う。コンフォメーション改変は、（３）で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３以上になるよう座標変化させる。

リガンドのコンフォメーションを大きく動かした最適化は、例えば（５）で最適化したエネルギー「Ｕ最適」でのコンフォメーションに対して、回転可能な結合をランダムに５つ選び、原子タイプごとに決められた回転角度間隔に従ってランダムに回転させる。その後（２）、（３）の過程を、５０００回繰り返し行う。

ただし、リガンドのコンフォメーションを変化させた後、リガンドの内部エネルギー「Ｕ内部」を計算しその値が５００．０以上のときはその後の計算をスキップし、次のリガンドコンフォメーションを発生させるようにする。

（７） （１）〜（６）までの過程を（０）で用意した複数の構造変化座標に対して行い、最適なタンパク質とリガンドとの複合体座標、最適エネルギー「Ｕ最適」を算出する。

（８） （１）〜（７）までの過程を（０）で用意した化合物データベース中の全てのリガンドに対して行い、化合物データベース中から該当タンパク質と結合する可能性のあるリガンドを選択する。

［３］ ［１］〜［２］によって記載された方法において、タンパク質の誘導適合を反映するパラメータおよび構造変化した立体構造座標を分子動力学計算方法を用いて算出する場合、該当タンパク質の立体構造に際し、基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、そのゆらぎの大きさを拘束条件として、分子動力学計算を行うことで、タンパク質の立体構造をエネルギー最適構造よりおおきく離れないようにして分子動力学計算を行う方法。

本手法による分子動力学計算は、例えば、基準振動計算より主鎖原子の２面角のゆらぎの値を算出し、該当のゆらぎ値を以下のように分子動力学計算における力の定数Ｋの部分に入れる。
Ｅｒｏｔ＝Ｋｒｏｔ（φ―φ０）²
Ｅｒｏｔはタンパク質の立体構造中において主鎖原子の２面角のエネルギーのことを示す。
φは主鎖原子の２面角。
φ０は主鎖原子の２面角の標準値。
Ｋｒｏｔの値が大きい場合はφはφ０に拘束される。
Ｅｐｏｓ＝Ｋｐｏｓ（ｒ−ｒ₀）²
Ｅｐｏｓはタンパク質の立体構造中において主鎖原子の位置のエネルギーのことを示す。
ｒは主鎖原子の座標。
ｒ０は主鎖原子の座標の標準値。
Ｋｐｏｓの値が大きい場合はｒはｒ０に拘束される。

[４] ［１］〜［２］によって記載された方法において、リガンドとタンパク質との相互作用を評価する際の目的関数として、従来の相互作用エネルギー関数に、タンパク質の動的性質を表現する関数を弾性エネルギーとして加え、タンパク質の立体構造座標から相互作用エネルギーを高速に算出するとともに、タンパク質の動的挙動に関する物理化学的性質を明確に描写することを特徴とする手法。

弾性エネルギーとして、タンパク質の局所的な柔らかさを考慮し、以下の関数「Ｕ衝突」として適応する。活性部位における動的挙動の少ない側鎖原子及び主鎖原子のみのｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義した。

ψ（ｉ，ｊ）＝Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²

Ｍは衝突を不許可とする活性部位の原子の数、Ｎはリガンドの原子の数、「Ｋ衝突」は１０００．０であることが好ましい。「Ｒ衝突（ｉ，ｊ）」は活性部位のｉ番目の原子とリガンドのｊ番目の原子それぞれのＶａｎ ｄｅｒ Ｗａａｌｓ半径の和とした。

ここで、活性部位の各原子に対し衝突を許す重み付けｗ（ｉ）が定義された場合以下の式を用いる。ただし、ｗ（ｉ）は、０〜１の範囲の実数とする。

ψ（ｉ，ｊ）＝ｗ（ｉ） ＊ Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²

Ｍは活性部位の原子の数、Ｎはリガンドの原子の数、「Ｋ衝突」は１０００．０であることが好ましい。「Ｒ衝突（ｉ，ｊ）」は活性部位のｉ番目の原子とリガンドのｊ番目の原子それぞれのＶａｎ ｄｅｒ Ｗａａｌｓ半径の和とした。

また、弾性エネルギーとしては、以下の関数を用いて定義することも可能である。

Ｅｖ＝ｗ（ｈａｒｄ ｓｈａｐｅ ｒｅｇｉｏｎ）、Ｅ＝０（ｓｏｆｔ ｓｈａｐｅ ｒｅｇｉｏｎ）

ｈａｒｄ ｓｈａｐｅ ｒｅｇｉｏｎとは、タンパク質の立体構造中、動的挙動の小さい部分であり、ｓｏｆｔ ｓｈａｐｅ ｒｅｇｉｏｎとは、動的挙動の大きい部分のことを指す。Ｗは定数で１００であることが好ましい。

[５] [１]〜[４]によって記載された方法において、タンパク質の動的性質を表現する関数として基準振動解析結果またはタンパク質の二次構造判定結果を用いることを特徴とした方法。

二次構造判定においては、タンパク質のへリックス、シート部分は揺らぎは小さいと考え、それ以外は揺らぎは大きいと考え相互作用の評価関数、分子動力学計算の拘束条件に適応する。

［６］ ［１］〜［５］によって記載された方法において、計算された該当タンパク質が代表的な複数の立体構造座標である場合や、該当タンパク質の立体構造が例えば核磁気共鳴スペクトルの解析結果のような複数の立体構造座標である場合についても、該当タンパク質と結合するリガンドの探索について、複数座標すべてを全自動的にかつ短時間で同等に評価することを可能とした方法。

本発明によれば、本発明はタンパク質の立体構造座標を用いたリガンド探索方法、詳しくはタンパク質立体構造座標が既知の場合、相互作用すると考えられるリガンドを予測する方法に関して、タンパク質とリガンドとの結合に大変重要であるタンパク質の動的挙動を反映したパラメータを取得し、かつタンパク質の動的挙動を反映したリガンドとの新規な相互作用評価関数を用いて、該当タンパク質の立体構造と結合する新規リガンドを予測を行うことができる。これにより従来法と比較して、より信頼性の高い、かつ医薬品設計等に適したタンパク質の立体構造を世界中で解析されている大量のゲノム配列に関しても対応するスピードで構築することができる。従来はｉｎ ｓｉｌｉｃｏ スクリーニングにおいては、タンパク質とリガンドとの相互作用に重要な誘導結合を充分に取り扱うことのできるアルゴリズムが見出されていなかった時点に対し、タンパク質とリガンドとの相互作用エネルギー関数に、基準振動計算結果、もしくは二次構造予測から得られるタンパク質の「ゆらぎ」を表すパラメータを簡易にとりこむ計算式を導入した。

さらに、分子動力学シミュレーションにおいては、この方法により、該当タンパク質の動的挙動を反映したパラメータとリガンドとの相互作用評価関数に関して、該当タンパク質についての基準振動計算を行い、その結果を分子動力学計算に反映させることを特徴とする。従来はタンパク質の動的挙動のシミュレーションを行うためには、分子動力学計算を用いていたが、従来法で該当タンパク質に対して分子動力学計算を行うとタンパク質立体構造はＸ線、ＮＭＲ等で解析された座標から大きくズレる。こうしたズレはタンパク質の動的挙動の物理化学的描写を含んでいるがＮＭＲ等で示される動的挙動の実験的な結果と矛盾する挙動となる場合があり、必ずしも精度の高いシミュレーションとならないことが多い。そこで、分子動力学計算を行う際には、タンパク質の立体構造をある程度固定しシミュレーションを行う必要があり、本手法では分子動力学計算におけるエネルギー関数中で主鎖原子の２面角に拘束をかける手法を開発した。さらに２面角の拘束条件としては、そのパラメータとして予め該当タンパク質の基準振動計算を行い、主鎖原子の２面角のゆらぎを算出し、そのゆらぎの大きさにより例えばゆらぎの大きい部分は拘束条件を緩め、ゆらぎの小さい部分は拘束条件を強めるパラメータとして用いることとした。こうした条件でタンパク質の分子動力学シミュレーションを行うことで、精度よく動的挙動を描写することができる。加えてこうして算出された分子シミュレーションからタンパク質の動的挙動を描写した座標を取得することができ、これを利用することでさまざまなリガンド結合部位の形状を用いたリガンド探索を行うことができる。

これらの結果、今までのｉｎ ｓｉｌｉｃｏ スクリーニングでは見出すことができなかった新規なリガンドを発見することを可能するとともに、今までは長時間を必要とする分子シミュレーションでしか解析できなかった「誘導結合」を含めたタンパク質−リガンドとの相互作用解析を短時間で行うことを可能にした。

こうしたアルゴリズムは、既存ソフトウェアよりも誘導結合現象をより深く考慮した“ｉｎ ｓｉｌｉｃｏ ｓｃｒｅｅｎｉｎｇ”に対応可能とし、誘導結合現象と疎水相互作用の正しい理解のもと単純化している。アルゴリズムは単純化されているので、自動化により多くのターゲットタンパク質を処理可能とする。その結果、例えば１００万以上の化合物データベースから、新規で、もっともらしい化合物を探索することができるので実験では対応できない規模のデータベースからもっともらしい化合物を現実的な時間内に探索することが出来る。

また、タンパク質−リガンドとの相互作用解析を短時間で行うことが可能になるので、例えば代謝、毒性の原因となる数多くのタンパク質と薬物との相互作用解析が可能となり、ｉｎ ｓｉｌｉｃｏでの薬物の代謝、毒性予測を行うことができる。

本発明において、リガンドとして取り扱うことのできる分子は使用するリガンドの種類や数を限定しないため、蛋白質、ペプチド、ＤＮＡ、薬剤成分、金属、イオン、糖類、核酸成分、ホルモンを含む全ての物質を当該リガンドと見なすことができる。この方法によって、具体的に農薬、医薬品等の分子設計を行うことができる。

リガンドとタンパク質との相互作用エネルギー評価関数には、従来ドッキング法では静電エネルギー項、ｖａｎ ｄｅｌ ｗａａｌｓ項、さらにはソフトドッキング法等に見られてる動的挙動を表現するための調整項が主に用いられているが、本手法においてはタンパク質とリガンドとの相互作用中にはソフトドッキング法等に見られてる動的挙動を表現するための調整項を用いる代わりに古典力学で用いられている弾性衝突の理論を適応し、タンパク質とリガンドとの相互作用に関して、その物理化学的性質をより明確にした。このことによりタンパク質の構造変化と相互作用との関係を得ることができ、リガンドの機能の理解を迅速かつ正確に行うための手助けとなる。

尚、本発明で利用するタンパク質の立体構造は、その際Ｘ線結晶構造解析等により、タンパク質の立体構造として３次元座標が決定されたもの以外に、タンパク質の経験的なモデリング法、特にホモロジーモデリング法或いはスレッディング法を利用して作成した立体構造座標をも適応することができる。

以下、本発明の実施の形態について説明する。好適かつ代表的な例について説明するものであり、本発明はこれに限定されるものではない。

本明細書において幾つかの用語を使用するが、特に明記しない限り、次の意味を有する。

「標的タンパク質」とは、立体構造の詳細がＸ線結晶解析やＮＭＲ解析、ホモロジーモデリング法により既に決定されており、リガンド探索の対象とするタンパク質を意味する。

「原子座標」とは、三次元空間上で立体構造を記述するものである。それは空間上のある点を原点とする互いに垂直な三方向の相対的な距離であり、タンパク質中に存在する水素原子を除く原子一つあたりに３個の数字からなるベクトル量である。

図１は、本発明によるタンパク質立体構造と誘導適合を利用したリガンド探索方法の一例を示すフローチャートである。図１は、本発明によるタンパク質の立体構造予測を示すフローチャート。

図１に示す通り、この方法においては、先ず、ステップ０において、リガンドのデータベースを用意する。リガンドのデータベースは分子力学法等に用いて３次元化することが望ましい。

ステップ１０において標的タンパク質の立体構造を選定し、座標を入手する。

ステップ２０において、標的タンパク質の基準振動計算を行い、主鎖原子の位置のゆらぎと２面角のゆらぎを求める。

ステップ３０において、ステップ２０において求めた標的タンパク質のゆらぎを拘束条件として用いた分子動力学計算を行う。

ステップ４０において、標的タンパク質のリガンド結合部位に、リガンドを配置するための点郡を指定する。

ステップ５０において、ステップ０で定められた一つのリガンドに対し、リガンドの各原子をステップ３０で定められた点群に重ね合わせる。

ステップ６０において、ステップ５０で定められた重ねあわせに対し、ステップ２０およびステップ３０で定められた計算結果によりタンパク質の動的挙動を表すパラメータを取得し、該当パラメータを用いてリガンドとタンパク質との相互作用エネルギーをリガンドコンフォメーションを微調整しながら計算する。コンフォメーションの微調整は、ステップ５０で定められたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３を越えない程度に座標変化させる。

ステップ７０において、ステップ５０で定められたリガンドに対して、コンフォメーションを大きく動かして、ステップ５０から再スタートを行い、ステップ７０までを繰り返して最適化を行う。コンフォメーション改変は、ステップ５０で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３以上になるよう座標変化させる。

ステップ８０において、ステップ７０まで得られた標的タンパク質とリガンドとの相互作用エネルギーを決定する。

ステップ９０において、ステップ４０に戻り、ステップ０中の別のリガンドを選択し、ステップ８０まで計算する。

ステップ１００において、ステップ０中でのリガンドに対し、ステップ９０において定められた相互作用エネルギーを比較し、標的タンパク質に結合すると予想されるリガンドを選択する。

本発明では、従来は標的タンパク質の立体構造が与えられた場合におけるリガンド探索が、タンパク質とリガンドとの相互作用にとって重要なタンパク質の動的な性質を反映させることが困難であった時点に対し、ステップ２０から９０までを行うことで、タンパク質とリガンドとの相互作用において、タンパク質の動的性質に関する運動エネルギー部分を弾性率として扱い、リガンドの近接に伴うタンパク質の弾性的性質を取り込んだエネルギー評価関数（安定化に寄与する疑似弾性エネルギー）を用いることで、相互作用エネルギーをｖａｎ ｄｅｌ ｗａａｌｔｈ相互作用エネルギーと疑似弾性エネルギーとの和の形式とし、物理化学的現象をはっきりさせた評価関数を用いることが従来法と異なる点である。以下、各ステップについて更に詳細に説明する。

（ステップ０：リガンドデータベースの準備）
３次元座標を含むリガンドデータベースを用意する。リガンドデータベースとしては、例えば、ＡＣＤ等のような市販化合物データベース、化合物を描いて収集した仮想化合物データを用いることができる。ガンドのデータベースは分子力学法等に用いて３次元化することが望ましい。

（ステップ１０：標的タンパク質立体構造の選択と取得）
ステップ０で定められたリガンドデータベースから特定リガンドを探索するための標的タンパク質を選択し、３次元座標を入手する。３次元座標は、公共データベースであるＰＤＢやホモロジーモデリング法等で作成した立体構造座標を用いることが望ましい。

（ステップ２０：標的タンパク質の基準振動計算）
ステップ１０で定められた参照タンパク質の動的挙動を表すパラメータを基準振動解析法による計算結果のデータベースもしくは二次構造判定計算をおこない取得する。まず、基準振動解析法によるタンパク質の動的挙動を表すパラメータ取得方法について下記に示す。

基準振動解析法とは、ポテンシャルエネルギーを変位の二次関数として近似し、運動方程式を厳密に解き、最適化構造の周りの微小な振動を解析する方法を意味する。解くべき運動方程式は下記式（１）または（２）である。

・・・（１）

・・・（２）
ただし、

ここでω_kは固有値、Ｕ_ikは固有ベクトルであり、δ_ijはクロネッカーのデルタである。Ｔ_ijとＶ_ijはそれぞれ運動エネルギーＥ_kとポテンシャルエネルギーＶに関係し、下記式（３）および（４）の通りである。

・・・（３）

・・・（４）

ここで、ｑ_iは振動の自由度に対応した座標、ｑ_i ⁰は最適化座標、ｑ_i’（式（３）における「ｑ_iドット」を意味する）はｑ_iの時間による微分である。Ａ_jkは集団運動Ｑ_kと個々の原子運動ｑ_jを結ぶ係数であり、下記式（５）の通りである。

・・・（５）
ただし、基準振動

である。ここで、α_kとδ_kは初期条件で定められる。

上記した基準振動解析法の詳細は「Ｗｉｌｓｏｎ，Ｅ．Ｂ．，Ｄｅｃｉｕｓ，Ｊ．Ｃ．，ａｎｄ Ｃｒｏｓｓ，Ｐ． Ｃ． １９５５． Ｍｏｌｅｃｕｌａｒ Ｖｉｂｒａｔｉｏｎ． ＭｃＧｒａｗ−Ｈｉｌｌ．」に記載されている。

参照タンパク質に対して、上記で得られた固有値、固有ベクトルを用いて、ある温度・ある固有値での各Ｃα原子の位置ゆらぎを計算し、このゆらぎの値をＣαが含まれるアミノ酸のゆらぎの値とする。目的タンパク質の各アミノ酸のゆらぎの値は、ステップ５０におけるアライメントを利用して、目的配列と参照配列の比較から対応するアミノ酸残基ペアにおいて、目的タンパク質のゆらぎの値として参照タンパク質と同一のものを当てはめておく。ゆらぎの値を求められなかったものについては、予め設定しておいた値をあてはめる。こうして得た目的タンパク質の各アミノ酸のゆらぎの値を目的タンパク質の動的な挙動を表すパラメータとする。

次に、二次構造判定計算によるタンパク質の動的挙動を表すパラメータ取得方法について下記に示す。

二次構造判定はタンパク質の立体構造座標から計算される。ソフトウェアとしては、ＤＳＳＰ、ＳＴＲＩＤＥ等が好ましいが、基本的にはタンパク質の主鎖のねじれ角と水素結合パターンから判別される方法を用いる。

ここで、「ＤＳＳＰ（Ｄｉｃｔｉｏｎａｒｙ ｏｆ ｐｒｏｔｅｉｎ ｓｅｃｏｎｄａｒｙ ｓｔｒｕｃｔｕｒｅ ｏｆ ｐｒｏｔｅｉｎ）」とは、ＰＤＢ書式のファイルを入力ファイルとして、主鎖の水素結合パターンと、内部回転角等を解析しαへリックスとβシートとを判定するソフトウェアである。ＤＳＳＰの詳細は、「Ｋａｂｓｃｈ，Ｗ． ＆ Ｓａｎｄｅｒ，Ｃ． （１９８３） Ｄｉｃｔｉｏｎａｒｙ ｏｆ ｐｒｏｔｅｉｎ ｓｅｃｏｎｄａｒｙ ｓｔｒｕｃｔｕｒｅ：ｐａｔｔｅｒｎ ｒｅｃｏｇｎｉｔｉｏｎ ｏｆ ｈｙｄｒｏｇｅｎ−ｂｏｎｄｅｄ ａｎｄ ｇｅｏｍｅｔｒｉｃａｌ ｆｅａｔｕｒｅｓ． Ｂｉｏｐｏｌｉｍｅｒｓ， ２２：２５７７−２６３７」に記載されている。

「ＳＴＲＩＤＥ（Ｐｒｏｔｅｉｎ ｓｅｃｏｎｄａｒｙ ｓｔｒｕｃｔｕｒｅ ａｓｓｉｇｎｍｅｎｔ ｆｒｏｍ ａｔｏｍｉｃ ｃｏｏｒｄｉｎａｔｅ）」とは、ＰＤＢ書式のファイルを入力ファイルとして、主鎖の水素結合パターンと、内部回転角等を解析しαへリックスとβシートとを判定するソフトウェアである。ＳＴＲＩＤＥの詳細は、「Ｆｒｉｓｈｍａｎ， Ｄ ＆ Ａｒｇｏｓ，Ｐ．（１９９５） Ｋｎｏｗｌｅｄｇｅ−ｂａｓｅｄ ｓｅｃｏｎｄａｒｙ ｓｔｒｕｃｔｕｒｅ ａｓｓｉｇｎｍｅｎｔ． Ｐｒｏｔｅｉｎｓ： ｓｔｒｕｃｔｕｒｅ， ｆｕｎｃｔｉｏｎ ａｎｄ ｇｅｎｅｔｉｃｓ， ２３， ５６６−５７９」に記載されている。

参照タンパク質に対して、上記ソフトウェア等を用いて、二次構造計算を行い、各アミノ酸がとるαへリックス構造、βシート構造、ループ構造を判定する。目的タンパク質の各アミノ酸の二次構造は、ステップ５０におけるアライメントを利用して、目的配列と参照配列の比較から対応するアミノ酸残基ペアにおいて、目的タンパク質の二次構造判定として参照タンパク質と同一のものを当てはめておく。二次構造判定を求められなかったものについては、予め設定しておいた結果をあてはめる。こうして得た目的タンパク質の各アミノ酸の二次構造判定結果を目的タンパク質の動的な挙動を表すパラメータとする。

上記の目的タンパク質の動的挙動を表すパラメータとしては、参照タンパク質の基準振動解析法より取得した結果を用いることが好ましく、該当計算結果は別途データベースとして保存されているものを使用する。二次構造判定計算結果は、好ましくは、基準振動解析が行われていない参照タンパク質を用いる際に基準振動解析計算の代用として使用する。

（ステップ３０：標的タンパク質の分子動力学計算（拘束つき））
主鎖の位置拘束エネルギー「Ｕ位置」を導入し、初期の受容体骨格の変動を抑えながらＡＰＲＩＣＯＴ［Ｙｏｎｅｄａ Ｓ． ＆ Ｕｍｅｙａｍａ Ｈ． （１９９２） Ｆｒｅｅ ｅｎｅｒｇｙ ｐｅｒｔｕｒｂａｔｉｏｎ ｃａｌｃｕｌａｔｉｏｎｓ ｏｎ ｍｕｌｔｉｐｌｅ ｍｕｔａｔｉｏｎ ｂａｓｅｓ Ｊ． Ｃｈｅｍ． Ｐｈｙｓ． ９７， ６７３０−６７３６］を用いて最小化（条件：温度３００Ｋ、受容体の表面から水分子が最低２分子配置できる箱状水槽、力場：ＡＭＢＥＲ［Ｓ． Ｊ． Ｗｅｉｎｅｒ， Ｐ．Ａ． Ｋｏｌｌｍａｎ， Ｄ．Ａ． Ｃａｓｅ， Ｕ．Ｃ． Ｓｉｎｇｈ， Ｃ． Ｇｈｉｏ， Ｇ． Ａｌａｇｏｎａ， Ｓ． Ｐｒｏｆｅｔａ， ＆， Ｐ． Ｗｅｉｎｅｒ （１９８４） Ａ ｎｅｗ ｆｏｒｃｅ ｆｉｅｌｄ ｆｏｒ ｍｏｌｅｃｕｌａｒ ｍｅｃｈａｎｉｃａｌ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｎｕｃｌｅｉｃ ａｃｉｄｓ ａｎｄ ｐｒｏｔｅｉｎｓ Ｊ． Ａｍ． Ｃｈｅｍ． Ｓｏｃ． １０６， ７６５−７８４］）を行った。

Ｕ位置 ＝ Ｋ位置 ＊ Ｒ²

ここで、「Ｋ位置」は３００．０、Ｒは基準座標からのずれとした。

続いて、ＡＰＲＩＣＯＴに二面角拘束エネルギー「Ｕ二面角」を導入して、最小化した受容体のＭＤ計算（条件：温度３００Ｋ、受容体の表面から水分子が最低２分子配置できる箱状水槽、力場：ＡＭＢＥＲ）を行った。

Ｕ二面角 ＝ Ｋ二面角 ＊ （θ−θ平衡）²

θは二面角（単位ｒａｄ）。

「Ｋ二面角」には、最大値と最小値を指定することで、その範囲内で主鎖二面角揺らぎに対応するように各二面角に対して不均一な拘束がかかるようにした。以後、主鎖二面角を拘束しながら行うＭＤを二面角拘束ＭＤと呼ぶことにする。

二面角拘束ＭＤ計算によりタンパク質構造座標を入手するには、受容体動的構造のクラスタリングを行う。

あらかじめ指定した活性部位に対して、ＭＤの途中経過１００ｆｓｅｃごとの受容体を重ね合わせた構造及び初期構造の活性部位を母集団とした。まず初めに、クラスタリングすることにより側鎖の動的情報が失われる可能性が高いことから、側鎖の二面角χにおいて母集団のα％が平均角度±２０．０°の範囲で保存されている全側鎖二面角を収集した。ただし、主鎖の根元に近い方からχが保存されていないと判定された場合はそれ以降のχは保存されていないものとした。

次に、収集した保存側鎖二面角をすべて網羅している構造を母集団から抽出した。そして、抽出した構造の類似性を比較するために全原子ｒｍｓ（ｒｏｏｔ ｍｅａｎ ｓｑｕａｒｅ）がβÅ以下の場合、同一構造と判断して一方を削除、最終的に選ばれた構造をもとに受容体動的構造クラスターを作成した。また、保存されていなかった二面角χを構成する原子では、変動する可能性が高いことから活性部位・リガンド結合計算において衝突しても良いことにした。ただし、α、βは定数。

（ステップ４０：基準となる標的タンパク質上の空間点を選択する。）
ステップ３０で作成した複数のタンパク質立体構造座標のうち１つをランダムに選択する。タンパク質座標中の空間点は例えば以下のような方法で発生させる。

１． ダミー原子の発生による空間点の発生
リガンドとタンパク質との相互作用における水素結合に着目し、タンパク質中の水素結合サイトを空間点として指定する。水素結合における重要事項は距離と角度である。つまり、角度を計算するためには、水素結合ドナー（以後、ドナー）に水素原子が必要になる。そこで、活性部位及びリガンドに水素原子が含まれていない場合、以下の規則によりダミー水素原子を発生させた。

１） ｓｐ²軌道原子を中心とする正三角形状にダミー原子を発生（図２）。

２） ｓｐ³軌道原子では、水素結合を形成する距離にある場合、水素原子を共有するように回転できると考え水素結合相互作用を計算するときには、距離のみを考慮することにした。このため、ｓｐ³軌道原子にはダミー原子を発生させない。

金属及び水では、活性部位・リガンド結合の仲介役となり得るので、相互作用する位置に以下のようにダミー原子を発生させた。

１） 鉄のような金属には、正八面体状にダミー原子を発生（図３）。

２） 水は、正四面体状にダミー原子を発生。

ただし、活性部位と相互作用している方向にはダミー原子を発生させないことにした。

２． 構造活性相関情報を利用した空間点の発生
リガンドの構造活性相関（ＳＡＲ）情報に着目し、以下の項目を入力情報にすることにした。

Ａ） ＳＡＲから得られた活性部位の原子（以後、「Ａ原子」）。ＰＤＢ形式に従う。

Ｂ） 「Ａ原子」と相互作用するであろうリガンドの原子タイプ（以後、「Ｂタイプ」）。ＳＹＢＹＬのＭＯＬ２形式に従う。

Ｃ） 「Ａ原子」と「Ｂタイプ」との相互作用の強さ（以後、「Ｃ強さ」）。

Ｄ） 「Ａ原子」と「Ｂタイプ」との相互作用する距離（以後、「Ｄ距離」）（単位はÅ）。

Ａ）〜Ｄ）をもとにタンパク質中の活性部位内におけるリガンドの初期座標を利用し、以下の規則により空間点を作成することにした。

１） 「Ａ原子」がドナーまたは金属及び水の場合、１．で発生させたダミー原子の方向に対して「Ａ原子」から「Ｄ距離」の位置及びその周囲を初期座標に選んだ（図４、図５）。

２） 「Ａ原子」がｓｐ³軌道原子の場合には、「Ａ原子」から「Ｄ距離」の周囲を初期座標に選んだ（図６）。

３） 「Ａ原子」が水素結合アクセプター（以後、アクセプター）の場合、「Ａ原子」の結合延長上「Ｄ距離」の位置及びその周囲を初期座標に選んだ（図７）。

４） その他の場合には、「Ａ原子」を中心とする半径が「Ｄ距離」の球表面上の点を初期座標に選んだ（図８）。

５） （１）〜（４）とは、異なりリガンドの初期座標を直接指定することもできるようにした。

（ステップ５０：ステップ４０の点群にステップ０中でリガンド座標の一つを重ね合わせる。）
距離行列を用いたアライメント作成アルゴリズム（ＤＡＬＩ）［Ｈｏｌｍ， Ｌ．， ＆ Ｓａｎｄｅｒ， Ｃ． （１９９３） Ｐｒｏｔｅｉｎ Ｓｔｒｕｃｔｕｒｅ Ｃｏｍｐａｒｉｓｏｎ ｂｙ Ａｌｉｇｎｍｅｎｔ ｏｆ Ｄｉｓｔａｎｃｅ Ｍａｔｒｉｃｅｓ Ｊ． Ｍｏｌ． Ｂｉｏｌ． ２３３， １２３−１３８］を低分子用に改良した手法で初期座標とリガンドを重ね合わせた。

（１） １つの「Ｂタイプ」には、リガンドの原子タイプが複数対応することが多い。そこで、乱数を用いて「Ｂタイプ」とリガンドの原子タイプで同一視できるペアを作成した。ただし、ペアにおいてリガンドの原子タイプが重複しないようにした。

（２） 「Ｂタイプ」には、ステップ４０−２．により複数の初期座標が含まれているので、初期座標も乱数を用いて選択した。

（３） 選択された初期座標とリガンドそれぞれの距離行列を作成し、Ｓｓｃｏｒｅ（ｉ，ｊ）を計算する。

ｄ_ij ^sは該当タンパク質中のｉ番目とｊ番目の空間点距離。ｄ_ij ^Cは化合物中のｉ番目とｊ番目の原子間距離。αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数。βは重なりと定義できる限界値を与えるための定数。

αは１．５、βは０．８とするのが好ましい。

（４） （１）〜（３）を１００００回繰り返し、Ｓｓｃｏｒｅ（ｉ，ｊ）が最大になるペアを探し、そのペア情報をもとにリガンドを初期座標に重ね合わせた。

（ステップ６０：標的タンパク質との相互作用関数をリガンドのコンフォメーションを微調整しながら最適化する。）
ステップ５０で重ねあわせたリガンドに対してタンパク質との相互作用エネルギーをコンフォメーションを微調整しながら最適化計算する。リガンドのコンフォメーションの微調整は、ステップ５０で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３を越えない程度に座標変化させる。

微調整は例えばランダムサーチで最適化することが好ましい。ランダムサーチは以下の項目に従ってタンパク質の活性部位とリガンドとの微小変化を８０００回行い、最適エネルギー「Ｕ最適」が最小になるようにする。

１〕 回転可能な結合のうち最大５つ乱数で選び、結合ごとにランダムに±１０．０°の範囲内で回転させリガンドのコンフォメーションを換える。この過程を３回に一度行う。
２〕 ｘ、ｙ、ｚ軸方向それぞれにおいて、ランダムに±１．０Åの範囲内でリガンドの並進運動を行う。この過程を２回に一度行う。

回転中心座標それぞれにおいて、ランダムに±１．０Åの範囲内で回転中心座標を移動させ、さらに３次元方向の角度それぞれに対して、ランダムに±５．０°の範囲内でリガンドの回転運動を行う。この過程を５回に一度行う。

最適エネルギーは以下のように定義する。

Ｕ最適＝Ｕ_SAR＋Ｕ水素＋Ｕ疎水＋Ｕスタッキング＋Ｕ衝突＋Ｕ内部

ここで、原子のＶａｎ ｄｅｒ Ｗａａｌｓ半径及び原子間相互作用距離はＡＭＢＥＲ９９［Ｊ． Ｗａｎｇ， Ｐ． Ｃｉｅｐｌａｋ ＆ Ｐ．Ａ． Ｋｏｌｌａｍ （２０００） Ｈｏｗ ｗｅｌｌ ｄｏｅｓ ａ ｒｅｓｔｒａｉｎｅｄ ｅｌｅｃｔｒｏｓｔａｔｉｃ ｐｏｔｅｎｔｉａｌ （ＲＥＳＰ） ｍｏｄｅｌ ｐｅｒｆｏｒｍ ｉｎ ｃａｌｃｕｌａｔｉｎｇ ｃｏｎｆｏｒｍａｔｉｏｎａｌ ｅｎｅｒｇｉｅｓ ｏｆ ｏｒｇａｎｉｃ ａｎｄ ｂｉｏｌｏｇｉｃａｌ ｍｏｌｅｃｕｌｅｓ？ Ｊ． Ｃｏｍｐｕｔ． Ｃｈｅｍ． ２１， １０４９−１０７４］およびＭＭ３パラメータ［Ｍａ Ｂ．， Ｌｉｉ Ｊ．−Ｈ．， Ａｌｌｉｎｇｅｒ Ｎ．Ｌ． （２０００） Ｍｏｌｅｃｕｌａｒ ｐｏｌａｒｉｚａｂｉｌｉｔｉｅｓ ａｎｄ ｉｎｄｕｃｅｄ ｄｉｐｏｌｅ ｍｏｍｅｎｔｓ ｉｎ ｍｏｌｅｃｕｌａｒ ｍｅｃｈａｎｉｃｓ Ｊ． Ｃｏｍｐｕｔ． Ｃｈｅｍ． ２１， ８１３−８２５］を参考にした。

（ａ） ＳＡＲに関するエネルギー関数
ＳＡＲ情報に従う指標としてエネルギーＵ_SARを定義した。

ψ（ｉ）＝Ｋ_SAR（ｉ） ＊ （Ｒ_SAR（ｉ）−Ｒ）²−δ

ＮはＳＡＲ情報の数、Ｒは「Ａ原子」からリガンド側の相互作用原子までの距離、Ｋ_SAR（ｉ）はｉ番目の「Ｃ強さ」、Ｒ_SAR（ｉ）はｉ番目の「Ｄ距離」、δは２０．０とした。

（ｃ） 疎水相互作用エネルギー
活性部位（ＡＬＡ、ＣＹＳ、ＰＨＥ、ＩＬＥ、ＬＥＵ、ＭＥＴ、ＰＲＯ、ＶＡＬ、ＴＲＰ、ＴＹＲの側鎖。ただし、ＴＹＲの水酸基は除く）及びリガンド（炭素原子）の疎水相互作用し得る原子に通し番号を付け、活性部位のｉ番目とリガンドのｊ番目との原子間距離Ｒがカットオフ以内にあるときψ（ｉ，ｊ）を計算するように定義した。

Ｍは活性部位の疎水相互作用し得る原子の数、Ｎはリガンドの疎水相互作用し得る原子の数、「Ｋ疎水（ｉ，ｊ）」及び「Ｒ疎水（ｉ，ｊ）」は原子タイプごとに決めた疎水相互作用の強さ及び距離とした。また、カットオフは８．０Åとした。

（ｅ） 分子間衝突（弾性衝突エネルギー）
活性部位（クラスタリングの際に定義した、保存されている側鎖原子及び主鎖原子のみ）のｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義した。

ψ（ｉ，ｊ）＝Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²

Ｍは衝突を不許可とする活性部位の原子の数、Ｎはリガンドの原子の数、「Ｋ衝突」は１０００．０とした。「Ｒ衝突（ｉ，ｊ）」は活性部位のｉ番目の原子とリガンドのｊ番目の原子それぞれのＶａｎ ｄｅｒ Ｗａａｌｓ半径の和とした。

ここで、活性部位の各原子に対し衝突を許す重み付けｗ（ｉ）が定義された場合以下の式を用いる。ただし、ｗ（ｉ）は、０〜１の範囲の実数とする。

ψ（ｉ，ｊ）＝ｗ（ｉ） ＊ Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²

Ｍは活性部位の原子の数、Ｎはリガンドの原子の数、「Ｋ衝突」は１０００．０であることが好ましい。「Ｒ衝突（ｉ，ｊ）」は活性部位のｉ番目の原子とリガンドのｊ番目の原子それぞれのＶａｎ ｄｅｒ Ｗａａｌｓ半径の和とした。

（ステップ７０：リガンドのコンフォメーションを大きく改変する。）
リガンドのコンフォメーションを大きく動かして、ステップ５０から再スタートを行い、ステップ６０までを繰り返して最適化を行う。コンフォメーション改変は、ステップ６０で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３以上になるよう座標変化させる。

リガンドのコンフォメーションを大きく動かした最適化は、例えばステップ５０で最適化したエネルギー「Ｕ最適」でのコンフォメーションに対して、回転可能な結合をランダムに５つ選び、原子タイプごとに決められた回転角度間隔に従ってランダムに回転させる。その後ステップ５０、ステップ６０の過程を、５０００回繰り返し行う。ただし、リガンドのコンフォメーションを変化させた後、リガンドの内部エネルギー「Ｕ内部」を計算しその値が５００．０以上のときはその後の計算をスキップし、次のリガンドコンフォメーションを発生させるようにする。

（ステップ８０：標的タンパク質とリガンドとの相互作用エネルギーの決定）
ステップ４０から７０まで「Ｕ最適」が最適値となる最適なタンパク質とリガンドとの複合体座標、最適エネルギー「Ｕ最適」を算出する。

（ステップ９０：ステップ４０に戻り、ステップ０中の別のリガンドを選択し、ステップ８０まで計算する。）
ステップ４０から９０まではステップ０中の化合物データベース中のリガンド全てについて行われる。

（ステップ１００：リガンドの選択）
ステップ９０まで評価されたタンパク質とリガンドとの複合体座標、最適エネルギー「Ｕ最適」を基にステップ０中のデータベース中のリガンドから該当タンパク質と結合する可能性のある化合物を選ぶ。

以下に、本発明にかかる相互作用関数によって該当タンパク質と結合する物質をコンピュータプログラムにより評価し、選定する方法及びその装置、該当タンパク質に結合するリガンドを選択する際に一連の処理を全自動または手動的に行うことを特徴とするコンピュータプログラム及びその装置、該当タンパク質の立体構造に際し、基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、そのゆらぎの大きさを拘束条件として、分子動力学計算を行うことで、タンパク質の立体構造をエネルギー最適構造よりおおきく離れないようにして分子動力学計算を行う方法及びその装置、タンパク質の動的挙動に関する物理化学的性質を明確に描写することを特徴とする手法及びその装置、タンパク質の動的性質を表現する関数として基準振動解析結果またはタンパク質の二次構造判定結果を用いることを特徴とした方法及びその装置、核磁気共鳴スペクトルの解析結果のような複数の立体構造座標である場合についても、該当タンパク質と結合するリガンドの探索について、複数座標すべてを全自動的にかつ短時間で同等に評価することを可能とした方法及びその装置を具体的な例を用いて詳細に説明する。
下記の実施例は、本発明の具体的な認識を得る一助と見るべきであり、本発明の範囲を何ら制限するものではない。

（二面角拘束ＭＤおよびクラスタリングにおけるパラメータ定数の決定）
基準振動解析により二面角のゆらぎ値が計算される。本発明においては、二面角のゆらぎ値を分子動力学計算における拘束条件として、

Ｕ二面角 ＝ Ｋ二面角 ＊ （θ−θ平衡）²
θは二面角（単位ｒａｄ）。

中の「Ｋ二面角」に適応する。実際には「Ｋ二面角」の最大値と最小値を指定することで、その範囲内で主鎖二面角揺らぎに対応するように各二面角に対して不均一な拘束がかかるようにしている。本実施例では、「Ｋ二面角」の適切な最大値と最小値を決定することを目的とする。

さらに分子動力学計算後、構造変化した座標をクラスター解析し代表構造を選択する。その際、あらかじめ指定した活性部位に対して、ＭＤの途中経過１００ｆｓｅｃごとの受容体を重ね合わせた構造及び初期構造の活性部位を母集団とし、まず初めに、クラスタリングすることにより側鎖の動的情報が失われる可能性が高いことから、側鎖の二面角χにおいて母集団のα％が平均角度±２０．０°の範囲で保存されている全側鎖二面角を収集した。ただし、主鎖の根元に近い方からχが保存されていないと判定された場合はそれ以降のχは保存されていないものとした。次に、収集した保存側鎖二面角をすべて網羅している構造を母集団から抽出した。そして、抽出した構造の類似性を比較するために全原子ｒｍｓ（ｒｏｏｔ ｍｅａｎ ｓｑｕａｒｅ）がβÅ以下の場合、同一構造と判断して一方を削除、最終的に選ばれた構造をもとに受容体動的構造クラスターを作成した。また、保存されていなかった二面角χを構成する原子では、変動する可能性が高いことから活性部位・リガンド結合計算において衝突しても良いことにした。ここでα、βは定数となるが、本実施例では、適切なα、βを決定することを目的とする。

ここでは、リガンドと接触している活性部位において最も良い主鎖の動的構造を得ることが目的であるので、ｒｍｓ（ｒｏｏｔ ｍｅａｎ ｓｑｕａｒｅ）を計算するときは活性部位における主鎖原子（Ｎ、Ｃα、Ｃ、Ｏ）の４原子のみが対象とした。

「Ｋ二面角」の最大値と最小値、クラスタリング定数α、βはＮＭＲで解析される構造を再現できる値が適切であると考えられる。ジヒドロ葉酸還元酵素（ＤＨＦＲ、ＰＤＢ ｃｏｄｅ：１ＬＵＤ）はＮＭＲで解析された構造である。そこでまず、ＮＭＲ構造のジヒドロ葉酸還元酵素（ＤＨＦＲ、ＰＤＢ ｃｏｄｅ：１ＬＵＤ）のＭＯＤＥＬ１を初期構造に基準振動解析を行い、ゆらぎ値を求めその後分子動力学計算をおこなった。また、１ＬＵＤ（ＭＯＤＥＬ１）に含まれていたリガンドの各原子から半径６Å以内に含まれる受容体残基を活性部位と定義した。分子動力学計算後は、受容体動的構造クラスタリングまで行った。ただし、ＭＤは０〜０．１ｎｓｅｃまでの結果を使用した。ここで、ＭＤでの拘束の最小値と最大値は０から１０００まで（１００ごと）、クラスタリングにおいては、定数αは０％から９０％まで（１０％ごと）、定数βはＮＭＲ構造平均値ｒｍｓを参考に、０．１Åから０．６Åまで（０．１Åごと）の数値を網羅的に行い、１ＬＵＤにおけるＮＭＲ構造すべてと比較することにより定数を決定することにした。

ＮＭＲ構造平均値は、受容体動的構造クラスタリングにおける定数βを決定するための参考として求めた。ｔｈｅ Ｐｒｏｔｅｉｎ Ｄａｔａ Ｂａｎｋ（ＰＤＢ）のＮＭＲ構造のうち、受容体が単純タンパク質で、１つのＰＤＢファイル内に記載されていたＮＭＲ構造が１０パターン以上あり、リガンドを含む１１７種類を対象に、活性部位のＮＭＲ構造平均値ｒｍｓを求めることにした。

まず、ＭＯＤＥＬ１においてリガンドの各原子から半径６Å以内に含まれる受容体残基を活性部位と定義した。ＭＯＤＥＬ１以外の構造において、ＭＯＤＥＬ１の活性部位とのｒｍｓをそれぞれ求め、さらにその平均ｒｍｓを求めた。ここで、平均ｒｍｓが１．０Å以上の場合は明らかな動的構造と見なせるので、そのようなＰＤＢファイルを対象からはずした。これにより、対象となるＰＤＢファイルは７１種類となった。７１種類の平均ｒｍｓをさらに平均化した値をＮＭＲ構造平均値ｒｍｓとした。このようにして得られたＮＭＲ構造平均値ｒｍｓは０．６２となった。

「Ｋ二面角」の適切な最大値と最小値、クラスタリングにおける定数α、βの決定に関しては、各パラメータ値とＮＭＲ構造との比較を行った。

１ＬＵＤには２４種類のＭＯＤＥＬが含まれており、ＭＯＤＥＬ１を対象にしたので、これを除く２３種類のＭＯＤＥＬの活性部位を正解構造とした。計算の結果出力された各受容体動的構造クラスターにおいて各正解構造とｒｍｓを計算しその中で最小のｒｍｓを「ＲＭＳ最小」として、各受容体動的構造クラスターから得られた「ＲＭＳ最小」の平均値をスコアとし、このスコアが最小となるパラメータを採用することにした。

図１１に１ＬＵＤのＭＯＤＥＬ１における基準振動解析結果を、図１２、図１３、図１５〜１８にスコアとパラメータとの比較の結果を示す。

図１１では、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示した。紫色Ｃは活性部位である。

図１２において、この結果、定数αは７０％が良いという結果になったが、一般性を持たせときに７０％ではクラスタリングの精度が低下した場合もあったので、これよりマイルドにした８０％を定数αの値にした。

図１３および図１５〜１８において、クラスタリング定数をα＝８０．０％、β＝０．４Åに固定した。なお、黒色に近いほどスコアが小さい。

これらの結果より、スコアが小さくなる拘束条件としては図１４の値が最適であると判断された。

これらの値の妥当性は、例えば、主鎖原子のみではなく、Ｃα原子、側鎖原子、全原子について調査しても図１４のパラメータ値が最適であることが分かる。

（拘束パラメータ有無による分子動力学計算の相違）
本発明によって開発された拘束パラメータを適応した分子動力学計算を２．０ｎｓｅｃまで行い、活性部位の主鎖原子の動的挙動が拘束パラメータを適応しない場合と比較して構造がどの程度変化するかを調べた。

Ｃａｓｅ１）
ジヒドロ葉酸還元酵素（１ＬＵＤのＭＯＤＥＬ１）を対象に検証した。結果を図１９〜図２１に示す。基準振動計算結果は。実施例１で求めた値を適応した。

図１９では、１ＬＵＤのＰＤＢファイル内に記載されている２４種類の各モデル構造をＭＯＤＥＬ１と活性部位の主鎖原子においてｒｍｓを計算し、その平均ｒｍｓを点線で表示。二面角拘束があるとき（Ａ）とないとき（Ｂ）において、活性部位の主鎖原子の初期構造からのずれをｒｍｓで表示。

図２０に二面角拘束なしでＭＤを計算させた時のＮＭＲ構造との比較を示す。図２０において、白色はＮＭＲ構造（１ｌｕｄ）であり、黒色はＭＤ構造（１ｌｕｄ）である。

表１に二面角拘束なしでＭＤを計算させた時のＮＭＲ構造との比較を示す。

図２１に二面角拘束ありでＭＤを計算させた時のＮＭＲ構造との比較を示す。

表２に二面角拘束ありでＭＤを計算させた時のＮＭＲ構造との比較を示す。

Ｃａｓｅ２）
ここでは、ＦＭＡＳ［Ｏｇａｔａ Ｋ．， Ｕｍｅｙａｍａ Ｈ． （２０００） Ａｎ ａｕｔｏｍａｔｉｃ ｈｏｍｏｌｏｇｙ ｍｏｄｅｌｉｎｇ ｍｅｔｈｏｄ ｃｏｎｓｉｓｔｉｎｇ ｏｆ ｄａｔａｂａｓｅ ｓｅａｒｃｈｅｓ ａｎｄ ｓｉｍｕｌａｔｅｄ ａｎｎｅａｌｉｎｇ Ｊ． Ｍｏｌ． Ｇｒａｐｈｉｃｓ Ｍｏｄ． １８， ２５８−２７２］によりモデリングした構造（モデル構造）とＸ線構造を初期構造に選び、初期構造及び拘束の有無に依存することを検証した。また、リガンドの各原子から半径１０Å以内に含まれる受容体残基を活性部位と定義した。

ｃｅｌｌｕｌａｒ ｒｅｔｉｎｏｉｃ ａｃｉｄ ｂｉｎｄｉｎｇ ｐｒｏｔｅｉｎ ｔｙｐｅ ＩＩ（ＣＲＡＢＰ−ＩＩ）（ＰＤＢ ｃｏｄｅ：１ＣＢＱ）のＸ線構造（立体構造）を利用した。また、参照タンパク質にホモロジー３２．１％のｉｎｔｅｓｔｉｎａｌ ｆａｔｔｙ ａｃｉｄ ｂｉｎｄｉｎｇ ｐｒｏｔｅｉｎ（ＰＤＢ ｃｏｄｅ：１ＩＣＭ ）を選び、図２２のアライメントでモデル構造を作成した。図２３、図２４、図２５にＸ線とモデルの構造比較を示す。

図２３には、１ＣＢＱの立体構造（Ｘ線構造（赤色Ａ）およびモデル構造（青色Ｂ））を示す。図２４には、図２３の緑色Ｃで示される物質である６−（２，３，４，５，６，７−ｈｅｘａｈｙｄｒｏ−２，４，４−ｔｒｉｍｅｔｈｙｌ−１−ｍｅｔｈｙｌｅｎｅｉｎｄｅｎ−２−ｙｌ）−３−ｍｅｔｈｙｌｈｅｘａ−２，４−ｄｉｅｎｏｉｃ ａｃｉｄの構造を示す。

図２５には、１ＣＢＱのＸ線構造とモデル構造の相違をｒｍｓで表示した。

図２６は１ＣＢＱのＸ線構造の基準振動解析の結果を示す図であり、図２７は１ＣＢＱのモデル構造の基準振動解析の結果を示す図である。図２６および図２７では、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。

図２８には、１ＣＢＱのＸ線構造とモデル構造の分子動力学法（ＭＤ）計算の結果を示す。Ｘ線構造の活性部位の主鎖原子とのｒｍｓを求めた。図２８において、Ａは初期構造がＸ線構造で二面角拘束なし、Ｂは初期構造がＸ線構造で二面角拘束あり、Ｃは初期構造がモデル構造で二面角拘束なし、Ｄは初期構造がモデル構造で二面角拘束あり、である。

Ｃａｓｅ３）
ＦｌａｖｏｄｏｘｉｎのＸ線構造（ＰＤＢ ｃｏｄｅ：１Ｊ９Ｇ）を利用した。また、参照タンパク質にホモロジー２９．２％のｆｌａｖｏｄｏｘｉｎ（ＰＤＢ ｃｏｄｅ：１ＡＨＮ）を選び、図２９のアライメントでモデル構造を作成した。図２９には、１Ｊ９Ｇおよび１ＡＨＮのアライメントを示す。

図３０には、１Ｊ９Ｇの立体構造（Ｘ線構造（赤色Ａ）およびモデル構造（青色Ｂ））を示す。図３１には、図３０における緑色Ｃで示される物質であるｆｌａｖｉｎ ｍｏｎｏｎｕｃｌｅｏｔｉｄｅの構造を示す。

図３２には、１Ｊ９ＧのＸ線構造とモデル構造の相違をｒｍｓで表示する。

図３３には１Ｊ９ＧのＸ線構造の基準振動解析の結果を、図３４には１Ｊ９Ｇのモデル構造の基準振動解析の結果を示す。図３３および図３４において、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥ［８］による二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。

図３５には、１Ｊ９ＧのＸ線構造とモデル構造の分子動力学法（ＭＤ）計算の結果を示す。Ｘ線構造の活性部位の主鎖原子とのｒｍｓを求めた。Ａは初期構造がＸ線構造で二面角拘束なし。Ｂは初期構造がＸ線構造で二面角拘束あり。Ｃは初期構造がモデル構造で二面角拘束なし。Ｄは初期構造がモデル構造で二面角拘束あり。

Ｃａｓｅ４）
Ｍａｔｒｉｘ ｍｅｔａｌｌｏｐｒｏｔｅｉｎａｓｅ−８（ＭＭＰ−８）のＸ線構造（ＰＤＢ ｃｏｄｅ：１ＭＭＢ）を利用した。また、参照タンパク質にホモロジー５５．０％のＭＭＰ−３（ＰＤＢ ｃｏｄｅ：１Ｂ３Ｄ）を選び、図３６のアライメントでモデル構造を作成した。図３６には、１ＭＭＢおよび１Ｂ３Ｄ＿Ａのアライメントを示す。

図３７には、１ＭＭＢの立体構造（Ｘ線構造（赤色Ａ）およびモデル構造（青色Ｂ））を示す。図３８には、図３７における緑色Ｃで示される物質であるｂａｔｉｍａｓｔａｔの構造を示す。

図３９には、１ＭＭＢのＸ線構造とモデル構造の相違をｒｍｓで表示する。

図４０には１ＭＭＢのＸ線構造の基準振動解析の結果を、図４１には１ＭＭＢのモデル構造の基準振動解析の結果を示す。図４０および図４１において、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥ［８］による二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。

図４２には、１ＭＭＢのＸ線構造とモデル構造の分子動力学法（ＭＤ）計算の結果を示す。Ｘ線構造の活性部位の主鎖原子とのｒｍｓを求めた。Ａは初期構造がＸ線構造で二面角拘束なし。Ｂは初期構造がＸ線構造で二面角拘束あり。Ｃは初期構造がモデル構造で二面角拘束なし。Ｄは初期構造がモデル構造で二面角拘束あり。

Ｃａｓｅ１）〜Ｃａｓｅ４）に示したとおり、拘束パラメータを適応した分子動力学計算結果は、拘束パラメータを適応しない場合と比較して、大きな構造変化は少ない。このことは古典力学を適応しているため大きな構造変化をしてしまう分子動力学法において、拘束パラメータを適応することで大きな構造変化を合理的に拘束することができ，理想的な構造座標を得ることが可能であるということを示している。また、ホモロジーが高ければ、ＦＭＡＳの構造構築精度もあがる。すなわち、Ｘ線に近い構造を得られるので、アミノ酸の数個異なるミューテーションタンパク質にもこの手法は利用できる。

（タンパク質／リガンド複合体モデルの検証）
本発明により該当タンパク質に結合するリガンドの複合体立体構造が予測される。本実施例では、こうして予測された複合体立体構造座標の予測精度を検証する。検証には、複合体の立体構造が既知で、リガンドの有無もしくはリガンドの種類により活性部位の形が異なるＩｎｄｕｃｅｄ−Ｆｉｔ型のタンパク質を用いた。ここで、リガンドの各原子から半径１０Å以内の残基をタンパク質の活性部位と定義した。また、Ｘ線構造またはＮＭＲ構造を初期構造に選んだＭＤでは、ほぼ一定の構造を保ち続けることが分かったのでＭＤを１．０ｎｓｅｃまで行うことにした。ただし、水素原子を除いて計算した。複合体モデル構築は、発明実施の形態に従って行った。

Ｃａｓｅ１）
ジヒドロ葉酸還元酵素（ＤＨＦＲ）である１ＢＺＦと１ＬＵＤとはホモロジーが１００．０％でかつ結合しているリガンドが異なることにより活性部位の形が異なる。そこで、１ＢＺＦ（ＭＯＤＥＬ１８）を初期構造として選択し、リガンドとして２，４−ｄｉａｍｉｎｏ−５−（３，４，５−ｔｒｉｍｅｔｈｏｘｙ−ｂｅｎｚｙｌ）−ｐｙｒｉｍｉｄｉｎ−１−ｉｕｍ（図４９）を用い、本発明によるプログラムによってタンパク質／リガンド複合体モデルを作成し、正解構造である１ＬＵＤ（ＭＯＤＥＬ４）と比較することで検証した（図４３）。

図４３に、ジヒドロ葉酸還元酵素の立体構造を示す。図４３では、１ＬＵＤ（ＭＯＤＥＬ４）受容体（緑色Ａ）とリガンド（赤色Ｂ）、１ＢＺＦ（ＭＯＤＥＬ１８）の受容体（青色Ｃ）とリガンド（水色Ｄ）を示す。

図４４に１ＢＺＦの基準振動計算解析を示す。図４４では、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。

図４５、図４７に１ＢＺＦを用いた拘束二面角分子動力学の結果を示す。図４５には、正解構造１ＬＵＤ（ＭＯＤＥＬ４）の活性部位とのｒｍｓを表示。Ａは主鎖原子。Ｂは側鎖原子。Ｃは全原子。図４７は、１ＢＺＦ（ＭＯＤＥＬ１８）における活性部位・リガンド結合解析であり、ＭＤ計算を０．１ｎｓｅｃまで行ったときの結合解析及び１．０ｎｓｅｃまで行ったとき、また、受容体動的構造クラスタリングを行うときの母集団を１００ｆｓｅｃごと及び１０００ｆｓｅｃごとで行ったときの結合解析。評価法は、正解構造との活性部位及びリガンドのｒｍｓで行った。

図４６にリガンドドッキングにおける空間点指定のパラメータ値を示す。図４６は、１ＬＵＤ（ＭＯＤＥＬ４）より得られた構造活性相関情報。

図４８〜図５０にタンパク／リガンド複合体と正解構造との比較を示す。図４８は、０〜１．０ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団により作られたクラスターを用いて行った活性部位・リガンド結合。緑色Ａは正解構造の１ＬＵＤ（ＭＯＤＥＬ４）。青色Ｂは初期構造の１ＢＺＦ（ＭＯＤＥＬ１８）。そして赤色Ｃはリガンド結合における最適構造。要素色Ｄはリガンドの正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは０．９６１４。図５０のリガンドを結合させることにより活性部位の主鎖原子のｒｍｓで０．２７９１の誘導が生じた。図４９は黒色：正解（１ｌｕｄのｍｏｄｅｌ ４）、灰色：初期構造（１ｂｚｆのｍｏｄｅｌ １８）、白色：最適構造。図５０は、２，４−ｄｉａｍｉｎｏ−５−（３，４，５−ｔｒｉｍｅｔｈｏｘｙ−ｂｅｎｚｙｌ）−ｐｙｒｉｍｉｄｉｎ−１−ｉｕｍ。１ＬＵＤのリガンド。

Ｃａｓｅ２）
ｈｅａｔ ｓｈｏｃｋ ｐｒｏｔｅｉｎ ９０（ＨＳＰ９０）である１ＹＥＲと１ＹＥＴはホモロジーが１００．０％でリガンド結合の有無により活性部位の形が異なる。そこで、リガンド結合していない１ＹＥＲを初期構造に選び、リガンドとしてｇｅｌｄａｎａｍｙｃｉｎを用い、正解構造である１ＹＥＴと比較することで検証した（図５１）。図５１は、ｈｅａｔ ｓｈｏｃｋ ｐｒｏｔｅｉｎ ９０の立体構造である。１ＹＥＴの受容体（緑色Ａ）とリガンド（赤色Ｂ）。１ＹＥＲの受容体（青色Ｃ）。

図５２に１ＹＥＲの基準振動解析を示す。二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。

図５３、図５５に１ＹＥＲを用いた拘束二面角分子動力学の結果を示す。正解構造１ＹＥＴの活性部位とのｒｍｓを表示。Ａは主鎖原子。Ｂは側鎖原子。Ｃは全原子。図５５は、１ＹＥＲにおける活性部位・リガンド結合解析。ＭＤ計算を０．１ｎｓｅｃまで行ったときの結合解析及び１．０ｎｓｅｃまで行ったとき、また、受容体動的構造クラスタリングを行うときの母集団を１００ｆｓｅｃごと及び１０００ｆｓｅｃごとで行ったときの結合解析。評価法は、正解構造との活性部位及びリガンドのｒｍｓで行った。

図５４にリガンドドッキングにおける空間点指定のパラメータ値を示す。図５４は、１ＹＥＴより得られた構造活性相関情報。

図５６および図５７にタンパク／リガンド複合体と正解構造との比較を示す。図５６および図５７は１ＹＥＲ・リガンド結合。図５６は、０〜０．１ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団により作られたクラスターを用いて行った活性部位・リガンド結合。緑色Ａは正解構造の１ＹＥＴ。青色Ｂは初期構造の１ＹＥＲ、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄはリガンドの正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは１．２０８１。図５７のリガンドを結合させることにより活性部位の主鎖原子のｒｍｓで０．１６１９の誘導が生じた。図５７は、ｇｅｌｄａｎａｍｙｃｉｎ。１ＹＥＴのリガンド。

Ｃａｓｅ３）
ｍｉｔｏｇｅｎ−ａｃｔｉｖａｔｅｄ ｐｒｏｔｅｉｎ ｋｉｎａｓｅ（ＭＡＰ ｋｉｎａｓｅ）である１Ａ９Ｕと１ＯＵＫはホモロジー１００．０％でかつ結合しているリガンドが異なることにより活性部位の形が異なる。そこで、１Ａ９Ｕを初期構造に選び、リガンドとして１ＯＵＫ中に含まれるリガンドを用い、正解構造である１ＯＵＫと比較することでと検証した（図５８）。図５８はｍｉｔｏｇｅｎ−ａｃｔｉｖａｔｅｄ ｐｒｏｔｅｉｎ ｋｉｎａｓｅの立体構造。図５８において、１ＯＵＫの受容体（緑色Ａ）とリガンド（赤色Ｂ）。１Ａ９Ｕの受容体（青色Ｃ）とリガンド（水色Ｄ）。

図５９に基準振動解析結果を示す。図５９は１Ａ９Ｕの基準振動解析。二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。

図６０、図６２に１ＹＥＲを用いた拘束二面角分子動力学の結果を示す。図６０は１Ａ９ＵのＭＤ。図６０に正解構造１ＯＵＫの活性部位とのｒｍｓを表示。Ａは主鎖原子。Ｂは側鎖原子。Ｃは全原子。図６２は１Ａ９Ｕにおける活性部位・リガンド結合解析。ＭＤ計算を０．１ｎｓｅｃまで行ったときの結合解析及び１．０ｎｓｅｃまで行ったとき、また、受容体動的構造クラスタリングを行うときの母集団を１００ｆｓｅｃごと及び１０００ｆｓｅｃごとで行ったときの結合解析。評価法は、正解構造との活性部位及びリガンドのｒｍｓで行った。

図６１にリガンドドッキングにおける空間点指定のパラメータ値を示す。図６１は１ＯＵＫより得られた構造活性相関情報。

図６３〜図６５にタンパク／リガンド複合体と正解構造との比較を示す。図６３〜図６５は１Ａ９Ｕ・リガンド結合。図６３は、０〜０．１ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団により作られたクラスターを用いて行った活性部位・リガンド結合。緑色Ａは正解構造の１ＯＵＫ。青色Ｂは初期構造の１Ａ９Ｕ、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄはリガンドの正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは１．６１１２。図６５のリガンドを結合させることにより活性部位の主鎖原子のｒｍｓで０．１８７１の誘導が生じた。図６４において、黒色：正解（１ｏｕｋ）、灰色：初期構造（１ａ９ｕ）、白色：最適構造。図６５は１ＯＵＫのリガンド。４−［５−［２−（１−ｐｈｅｎｙｌ−ｅｔｈｙｌａｍｉｎｏ）−ｐｙｒｉｍｉｄｉｎ−４−ｙｌ］−１−ｍｅｔｈｙｌ−４−（３−ｔｒｉｆｌｕｏｒｏｍｅｔｈｙｌｐｈｅｎｙｌ）−１Ｈ−ｉｍｉｄａｚｏｌ−２−ｙｌ］−ｐｉｐｅｒｉｄｉｎｅ。

Ｃａｓｅ１）〜Ｃａｓｅ３）に示す通り、本発明により作成されるタンパク質／リガンド複合体モデルは、誘導結合型のタンパク質／リガンド複合体の立体構造を精度よく予測可能であることが分かる。

（Ｆｘａを用いたｉｎ ｓｉｌｉｃｏ Ｓｃｒｅｅｎｉｎｇへの応用例）
本発明により、セリンプロテアーゼの１種であるＦｘａの立体構造（図６６）を用い、化合物データベースからＦｘａに結合する可能性のあるリガンドを探索した。立体構造には、１ＡＩＸを用い、リガンドデータベースとして、ＰＤＢデータベースより収集した３６３３種類のリガンドを用いた。発明実施の形態に従い、ｉｎ ｓｉｌｉｃｏ ｓｃｒｅｅｎｉｎｇを行った。その結果を図６７に示す。

図６７は、化合物データベース中のリガンドのうち、１ＡＩＸとの相互作用エネルギーの上位１００個を示している。図６７で、太字は１ＡＩＸ中に含まれているリガンドで、斜線はセリンプロテアーゼ。ＰＤＢ ｃｏｄｅとは、リガンドが含まれているもとのＰＤＢｃｏｄｅを示す。図６７には、１ＡＩＸにもともと含まれているリガンドがランキング１９位に入っている。

ランキング１９位におけるタンパク質／リガンド複合体構造とを図６８および図６９に示す。図６８において、白色は受容体、黒色は１ＡＩＸのリガンド。図６９は１ＡＩＸ中のリガンド。

図６７中のランキング３５位、３８位、８０位はすべてセリンプロテアーゼに結合するリガンドである。

これらの構造とおけるタンパク質／リガンド複合体構造を図７０および図７１、図７２および図７３、図７４および図７５に示す。図７０および図７１はランキング３５位におけるタンパク質／リガンド複合体構造。図７０において、白色は受容体、黒色は１ＡＵＪのリガンド。図７１は１ＡＵＪ中のリガンド。図７２および図７３はランキング３８位におけるタンパク質／リガンド複合体構造。図７２において、白色は受容体、黒色は正解（１ＦＯＲ）のリガンド、ＲＭＳは１．５００。図７３は１ＦＯＲ中のリガンド。図７４および図７５はランキング８０位におけるタンパク質／リガンド複合体構造。図７４において、白色は受容体、黒色は１Ｋ１Ｍのリガンド。図７５は１Ｋ１Ｍ中のリガンド。

これらの結果から、本発明により、化合物データベースからもっともらしい化合物を選択することが可能であることが分かる。

（異なる条件でのｉｎ ｓｉｌｉｃｏスクリーニング）
構造活性相関（ＳＡＲ）の情報により順位が変動することを検証する。また、受容体を固定した場合の順位の変動も検証する。

ここでは、ｓｅｖｅｒｅ ａｃｕｔｅ ｒｅｓｐｉｒａｔｏｒｙ ｓｙｎｄｒｏｍｅ（ＳＡＲＳ）のプロテアーゼを用いたｉｎ ｓｉｌｉｃｏ ｓｃｒｅｅｎｉｎｇを行った。初期構造にはリガンドを含まない１ＵＫ３（Ｂ鎖）を、またリガンドを含む１ＵＫ４（Ｂ鎖）のリガンド結合様式を構造活性相関情報として利用した。活性部位は１ＵＫ４（Ｂ鎖）のリガンドの各原子から半径１０Å以内に含まれる受容体残基部位。リガンドデータベースとして、ＰＤＢより収集した３６３３種類のリガンドを用いた。ただし、結合解析で利用する受容体動的構造クラスターには０〜０．１ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団で作られたものを使用した。また、水素原子を除いて計算した。

図７６はＳＡＲＳプロテアーゼの立体構造。１ＵＫ４（Ｂ鎖）の受容体（緑色Ａ）とリガンド（赤色Ｂ）。１ＵＫ３（Ｂ鎖）の受容体（青色Ｃ）。

図７７に１ＵＫ３の基準振動解析の結果を示す。図７７は１ＵＫ３（Ｂ鎖）の基準振動解析。二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。

図７８に１ＵＫ３の分子動力学計算の結果を示す。図７８は１ＵＫ３（Ｂ鎖）のＭＤ。１ＵＫ４（Ｂ鎖）の活性部位とのｒｍｓを表示。Ａは主鎖原子。Ｂは側鎖原子。Ｃは全原子。

Ｃａｓｅ１） ＳＡＲ４ヵ所指定
図７９に活性部位内での空間指定を示す。図７９は１ＵＫ４（Ｂ鎖）より得られた構造活性相関情報。

図８０にｉｎ ｓｉｌｉｃｏ スクリーニングの結果を示す。図８０は１ＵＫ３（Ｂ鎖）におけるｉｎ ｓｉｌｉｃｏスクリーニングの結果。

図８１に正解構造との比較を示す。図８１は１ＵＫ３と１ＵＫ４との比較。順位２５。緑色Ａは１ＵＫ４（Ｂ鎖）。青色Ｂは初期構造の１ＵＫ３（Ｂ鎖）、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄは１ＵＫ４のペプチド性リガンド（ＡＳＮ−ＳＥＲ−ＴＨＲ−ＬＥＵ−ＧＬＮ）の正解構造。Ｅは計算結果によるリガンド。リガンドのｒｍｓは２．５７２１。正解構造との活性部位の主鎖原子のｒｍｓは、初期構造では１．０２４８、最適構造では１．０７９２。

図８２および図８３に、ｉｎ ｓｉｌｉｃｏ スクリーニングの順位１を示す。図８３は１ＱＦ４のリガンド（Ｃ８−Ｒ）−ｈｙｄａｎｔｏｃｉｄｉｎ ５’−ｐｈｏｓｐｈａｔｅ。

Ｃａｓｅ２） ＳＡＲ３ヵ所指定
図８４に１ＵＫ３の活性部位内での空間指定を示す。図８４は１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。

図８５には、１ＵＫ３での最適構造と正解構造との比較を示す。図８５は１ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。順位４９。緑色Ａは１ＵＫ４（Ｂ鎖）。青色Ｂは初期構造の１ＵＫ３（Ｂ鎖）、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄは１ＵＫ４のペプチド性リガンド（ＡＳＮ−ＳＥＲ−ＴＨＲ−ＬＥＵ−ＧＬＮ）の正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは２．００５７。正解構造との活性部位の主鎖原子のｒｍｓは、初期構造では１．０２４８、最適構造では１．０４６９。

図８６には、ＳＡＲ３ヵ所指定でのｉｎ ｓｉｌｉｃｏスクリーニング結果を示す。図８６はＳＡＲ３ヵ所指定で実行したｉｎ ｓｉｌｉｃｏスクリーニングの結果。

Ｃａｓｅ３） ＳＡＲ５ヵ所指定
図８７に１ＵＫ３の活性部位内での空間指定を示す。図８７は１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。

図８８にはＳＡＲ５ヵ所指定でのｉｎ ｓｉｌｉｃｏスクリーニング結果を示す。図８８はＳＡＲ５ヵ所指定で実行したハイスループットスクリーニングの結果。

図８９には、１ＵＫ３での最適構造と正解構造との比較を示す。図８９は１ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。順位２。緑色Ａは１ＵＫ４（Ｂ鎖）。青色Ｂは初期構造の１ＵＫ３（Ｂ鎖）、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄは１ＵＫ４のペプチド性リガンド（ＡＳＮ−ＳＥＲ−ＴＨＲ−ＬＥＵ−ＧＬＮ）の正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは１．２５７８。正解構造との活性部位の主鎖原子のｒｍｓは、初期構造では１．０２４８、最適構造では１．１６２０。

Ｃａｓｅ４） リガンド原子タイプ指定の変更
図９０に１ＵＫ３の活性部位内での空間指定を示す。図９０は１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。

図９１には、リガンド原子タイプ指定変更でのｉｎ ｓｉｌｉｃｏスクリーニング結果を示す。図９１はリガンド原子タイプ指定変更で実行したハイスループットスクリーニングの結果。

図９２には、１ＵＫ３での最適構造と正解構造との比較を示す。図９２は１ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。順位７７４。緑色Ａは１ＵＫ４（Ｂ鎖）。青色Ｂは初期構造の１ＵＫ３（Ｂ鎖）、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄは１ＵＫ４のペプチド性リガンド（ＡＳＮ−ＳＥＲ−ＴＨＲ−ＬＥＵ−ＧＬＮ）の正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは２．５２１６。正解構造との活性部位の主鎖原子のｒｍｓは、初期構造では１．０２４８、最適構造では１．０７９２。

Ｃａｓｅ５） 受容体固定
図９３に活性部位内での空間指定を示す。図９３は１ＵＫ４（Ｂ鎖）より得られた構造活性相関情報。

図９４に、受容体を固定したｉｎ ｓｉｌｉｃｏ スクリーニングの結果を示す。図９４は受容体を固定した状態で実行したハイスループットスクリーニングの結果。

図９５には、１ＵＫ３と１ＵＫ４を重ね合わせたリガンドと計算結果のリガンドとの比較を示す。順位３９。灰色は１ＵＫ３の活性部位構造。黒色は１ＵＫ３と１ＵＫ４を重ね合わせたリガンド。白色は計算結果のリガンド。

Ｃａｓｅ１）〜Ｃａｓｅ４）を見ると、ＳＡＲの指定が多いほど参考にしたリガンドの順位が良くなる。つまり、参考にできるリガンドの結合情報が信頼できる場合には、ＳＡＲの情報を多くしたｉｎ ｓｉｌｉｃｏスクリーニングを行い、信頼性に欠ける場合には、ＳＡＲの情報数を減らし、さらに、リガンド原子タイプ指定の幅を広げることで、様々なリガンドがランキング上位に分布する。そして、その分布情報をもとにＳＡＲ情報を作り変えてｉｎ ｓｉｌｉｃｏスクリーニングを実行するとより信頼性の持てる結果が出力されるはずである。

Ｃａｓｅ１）とＣａｓｅ５）を見ると、受容体の動的構造の有無による順位変動を示している。これは、リガンドの動きのみの最適化に比べ、リガンド及び受容体それぞれが動く最適化の方が原子のぶつかりをさけることに優れている。従って、同じ位置に配置するための最適化エネルギーに差が生じる。

（二面角拘束分子動力学計算のパラメータに関するＭＤパラメータの分布）
ＦＭＮ−ｂｉｎｄｉｎｇ ｐｒｏｔｅｉｎにおける二面角拘束ＭＤパラメータの分布。
ここでは、１ＬＵＤ以外のＮＭＲ構造でも二面角拘束分子動力学計算及びクラスタリングのパラメータが同様の結果を生じるのかを検証する。そこで、ＦＭＮ−ｂｉｎｄｉｎｇ ｐｒｏｔｅｉｎのＮＭＲ構造（ＰＤＢ ｃｏｄｅ：１ＡＸＪ）のＭＯＤＥＬ１を初期構造に選んだ。評価法は、受容体動的構造クラスタリングのパラメーター（α＝８０．０％、β＝０．４Å）を固定したこと以外は実施例１に従った。

図９６に１ＡＸＪにおける二面角拘束分子動力学計算のパラメータ決定のスコアの分布状況を示す。図９６は１ＡＸＪにおける二面角拘束ＭＤパラメータの分布。Ａに近い部分ほどスコアで小さい。１ＬＵＤの時と同様に二面角拘束の最大値８００、最小値０では良い結果を示す。

（二面角拘束ＭＤ）
ここでは、主鎖二面角拘束ＭＤで各原子の動的構造を検証する。また、時として基準振動解析が収束せず二面角揺らぎ情報が得られないことがある。そこで、図１３により主鎖二面角に対して均一な拘束（５００）でＭＤを行ったときでも良い結果になっていることから、この場合における動的構造も検証する。実施例１に従い、拘束なし、二面角揺らぎを用いた拘束及び均一な拘束（５００）の条件のもとＭＤを行った。

図９７〜図１０８に１ＬＵＤに対して行った分子動力学計算の各原子における動的挙動の結果を示す。図９７〜図１０８は１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。図９７および図９８は活性部位の主鎖原子、図９９および図１００は受容体の主鎖原子、図１０１および図１０２は活性部位の側鎖原子、図１０３および図１０４は受容体の側鎖原子、図１０５および図１０６は活性部位の全原子、図１０７および図１０８は受容体の全原子、における動的挙動の結果を示す図である。１ＬＵＤのＰＤＢファイル内に記載されている２４種類の各モデル構造をＭＯＤＥＬ１と活性部位の主鎖原子においてｒｍｓを計算し、その平均ｒｍｓを点線で表示。二面角拘束があるとき（Ａ）とないとき（Ｂ）及び二面角拘束が５００で一定（Ｃ）において、活性部位の主鎖原子の初期構造からのずれをｒｍｓで表示。

ここには、記載しないが１ＣＢＱ、１Ｊ９Ｇ、１ＭＭＢ、１ＢＺＦ（ＭＯＤＥＬ１８）、１ＹＥＲ、１Ａ９Ｕ及び１ＵＫ３（Ｂ鎖）に関しても主鎖二面角揺らぎに基づく拘束ＭＤの結果を見ると、図１０９〜図１１１と同様に主鎖原子の抑制があると、拘束のない側鎖原子にも一定動きを示す。受容体の動きにおいて主鎖原子の動きの比重が大きいことが理解できる。

（異なる条件での結合解析）
ここでは、二面角拘束ＭＤ及びクラスタリングのパラメータが異なっても誘導が生じることを検証する。拘束の最大値１００、最小値０及び受容体動的構造クラスタリング定数α＝８０．０％、β＝１．０Åに設定して、その他は実施例２に従った。ただし、受容体動的構造クラスターには０〜０．１ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団で作られたものを使用した。また、活性部位の定義は、リガンドの各原子から半径６Å以内にある受容体残基とした。

図１０９〜図１１１には、異なる条件で受容体／リガンド結合の結果を示す。
（ｉ） １ＢＺＦ（ＭＯＤＥＬ１８）で、リガンド結合により活性部位の主鎖原子ｒｍｓで０．２６８６の誘導が生じた（図１０９）。活性部位全体のｒｍｓでは０．１２２４の誘導。リガンドのｒｍｓは０．８５２６。
（ｉｉ） １ＹＥＲで、リガンド結合により活性部位の主鎖原子ｒｍｓで０．２３７６の誘導が生じた（図１１０）。活性部位全体のｒｍｓでは０．０８１６の誘導。リガンドのｒｍｓは０．７２４６。
（ｉｉｉ） １Ａ９Ｕで、リガンド結合により活性部位の主鎖原子ｒｍｓで０．２１５０の誘導が生じた（図１１１）。活性部位全体のｒｍｓでは０．０４６４の誘導。リガンドのｒｍｓは０．９４６４。

ただし、緑色は正解構造、青色は初期構造、赤色は最適構造。要素色は正解リガンド、水色は最適リガンド。

図１０９〜図１１１で示すように、各条件が異なっていても、与えられた条件の中で最適な結果を生じることができる。

（正解構造を初期構造に選んだときの結合解析）
ここでは、ＤＨＦＲの１ＢＺＦ及び１ＬＵＤはリガンドの結合様式が似ているので、構造活性相関情報を一部変更し１ＢＺＦのリガンドの結合解析を行う。条件としては、初期構造に１ＢＺＦ（ＭＯＤＥＬ１８）、０〜０．１ｎｓｅｃまでの母集団より作成したクラスターを使用した。

図１１２には、１ＢＺＦの活性部位内における空間指定を示す。図１１２は１ＢＺＦ用に変更した構造活性相関情報。

図１１３および図１１４には、１ＢＺＦにおける受容体／リガンド結合の結果を示す。図１１３および図１１４は１ＢＺＦ（ＭＯＤＥＬ１８）のリガンド結合解析。
（ｉ） 最適化したときの受容体には初期構造を選択し、リガンドのｒｍｓは０．８８８４。（図１１３）
（ｉｉ） ｔｒｉｍｅｔｒｅｘａｔｅ。１ＢＺＦ（ＭＯＤＥＬ１８）のリガンド。（図１１４）

初期構造が元々ＰＤＢに登録されていた構造、つまり最適構造であったため、計算結果でもそれが図１１３および図１１４のように再現できる。

以上、本発明の方法は、医農薬の分子設計等を中心に、受容体/リガンド結合の解析を行う分野（医薬品設計）において、極めて有用であると考えられる。本発明は、産業上多くの分野、特に医薬品、食品、化粧品、医療、構造解析、機能解析等の分野で広く実施することができ、故に極めて有用である。

本発明によるタンパク質立体構造と誘導適合を利用したリガンド探索方法の一例を示すフローチャートである。 ｓｐ²軌道原子におけるダミー水素原子発生。 金属原子におけるダミー原子発生。 構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。 構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。 構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。 構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。 構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。 水素結合角の定義。 スタッキングにおける角度の定義。 １ＬＵＤ（ＭＯＤＥＬ１）の基準振動解析の結果。 ＭＤ及びクラスタリングのパラメータとスコア。 クラスタリング定数を固定した時のＭＤでの二面角拘束の最大値、最小値の分布。 拘束パラメータ。 クラスタリングのパラメータを固定した時の二面角拘束パラメータの分布。 クラスタリングのパラメータを固定した時の二面角拘束パラメータの分布。 クラスタリングのパラメータを固定した時の二面角拘束パラメータの分布。 クラスタリングのパラメータを固定した時の二面角拘束パラメータの分布。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤの結果。 二面角拘束なしでＭＤを計算させた時のＮＭＲ構造との比較。 二面角拘束ありでＭＤを計算させた時のＮＭＲ構造との比較。 １ＣＢＱのアライメント。 １ＣＢＱの立体構造。 １ＣＢＱの立体構造。 １ＣＢＱのＸ線構造とモデル構造の相違をｒｍｓで表示。 １ＣＢＱのＸ線構造とモデル構造の基準振動解析。 １ＣＢＱのＸ線構造とモデル構造の基準振動解析。 １ＣＢＱのＸ線構造とモデル構造のＭＤ。 １Ｊ９Ｇのアライメント。 １Ｊ９Ｇの立体構造。 １Ｊ９Ｇの立体構造。 １Ｊ９ＧのＸ線構造とモデル構造の相違をｒｍｓで表示。 １Ｊ９ＧのＸ線構造とモデル構造の基準振動解析。 １Ｊ９ＧのＸ線構造とモデル構造の基準振動解析。 １Ｊ９ＧのＸ線構造とモデル構造のＭＤ。 １ＭＭＢのアライメント。 １ＭＭＢの立体構造。 １ＭＭＢの立体構造。 １ＭＭＢのＸ線構造とモデル構造の相違をｒｍｓで表示。 １Ｊ９ＧのＸ線構造とモデル構造の基準振動解析。 １Ｊ９ＧのＸ線構造とモデル構造の基準振動解析。 １ＭＭＢのＸ線構造とモデル構造のＭＤ。 ジヒドロ葉酸還元酵素の立体構造。 １ＢＺＦ（ＭＯＤＥＬ１８）の基準振動解析。 １ＢＺＦ（ＭＯＤＥＬ１８）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ４）より得られた構造活性相関情報。 １ＢＺＦ（ＭＯＤＥＬ１８）における活性部位・リガンド結合解析。 １ＢＺＦ（ＭＯＤＥＬ４）・リガンド結合。 １ＢＺＦ（ＭＯＤＥＬ４）・リガンド結合。 １ＢＺＦ（ＭＯＤＥＬ４）・リガンド結合。 ｈｅａｔ ｓｈｏｃｋ ｐｒｏｔｅｉｎ ９０の立体構造。 １ＹＥＲの基準振動解析。 １ＹＥＲのＭＤ。 １ＹＥＴより得られた構造活性相関情報。 １ＹＥＲにおける活性部位・リガンド結合解析。 １ＹＥＲ・リガンド結合。 １ＹＥＲ・リガンド結合。 ｍｉｔｏｇｅｎ−ａｃｔｉｖａｔｅｄ ｐｒｏｔｅｉｎ ｋｉｎａｓｅの立体構造。 １Ａ９Ｕの基準振動解析。 １Ａ９ＵのＭＤ。 １ＯＵＫより得られた構造活性相関情報。 １Ａ９Ｕにおける活性部位・リガンド結合解析。 １Ａ９Ｕ・リガンド結合。 １Ａ９Ｕ・リガンド結合。 １Ａ９Ｕ・リガンド結合。 １ＡＩＸの立体構造。 ｉｎ ｓｉｌｉｃｏ ｓｃｒｅｅｎｉｎｇの結果。 タンパク質／リガンド複合体構造。 タンパク質／リガンド複合体構造。 タンパク質／リガンド複合体構造。 タンパク質／リガンド複合体構造。 タンパク質／リガンド複合体構造。 タンパク質／リガンド複合体構造。 タンパク質／リガンド複合体構造。 タンパク質／リガンド複合体構造。 ＳＡＲＳプロテアーゼの立体構造。 １ＵＫ３（Ｂ鎖）の基準振動解析。 １ＵＫ３（Ｂ鎖）のＭＤ。 １ＵＫ４（Ｂ鎖）より得られた構造活性相関情報。 １ＵＫ３（Ｂ鎖）におけるｉｎ ｓｉｌｉｃｏスクリーニングの結果。 １ＵＫ３と１ＵＫ４との比較。 ｉｎ ｓｉｌｉｃｏスクリーニングの順位１。 ｉｎ ｓｉｌｉｃｏスクリーニングの順位１。 １ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。 １ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。 ＳＡＲ３ヵ所指定で実行したｉｎ ｓｉｌｉｃｏスクリーニングの結果。 １ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。 ＳＡＲ５ヵ所指定で実行したハイスループットスクリーニングの結果。 １ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。 １ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。 リガンド原子タイプ指定変更で実行したハイスループットスクリーニングの結果。 １ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。 １ＵＫ４（Ｂ鎖）より得られた構造活性相関情報。 受容体を固定した状態で実行したハイスループットスクリーニングの結果。 １ＵＫ３と１ＵＫ４を重ね合わせたリガンドと計算結果のリガンドとの比較。 １ＡＸＪにおける二面角拘束ＭＤパラメータの分布。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 １ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。 異なる条件で受容体／リガンドの結合結果。 異なる条件で受容体／リガンドの結合結果。 異なる条件で受容体／リガンドの結合結果。 １ＢＺＦ用に変更した構造活性相関情報。 １ＢＺＦ（ＭＯＤＥＬ１８）のリガンド結合解析。 １ＢＺＦ（ＭＯＤＥＬ１８）のリガンド結合解析。

Claims (12)

1. 単数または複数鎖のタンパク質の座標データが与えられた場合に、当該タンパク質と結合するリガンドをコンピュータで下記のステップを実行することにより探索するリガンド
   探索方法であって、
   上記コンピュータは、処理手段と記憶手段とを少なくとも備え、
   上記記憶手段は、
   上記タンパク質の座標データと、
   上記リガンドのリガンド座標データと、
   上記タンパク質の活性部位の原子、当該原子と相互作用する上記リガンドの原子タイプ、上記原子と上記原子タイプとの相互作用の強さ、および、上記原子と上記原子タイプとの相互作用する距離、を少なくとも含む構造活性相関情報と、
   を少なくとも記憶し、
   上記処理手段において実行される、
   上記記憶手段に記憶された上記タンパク質の上記座標データに対して、基準振動計算を行うことにより、各アミノ酸のゆらぎ値を、誘導適合を反映した誘導適合パラメータとして求め、当該誘導適合パラメータを拘束条件として用いた分子動力学計算を行うことにより、上記タンパク質の動的挙動を考慮した構造変化後タンパク質座標データを選択する構造変化後タンパク質座標データ選択ステップと、
   上記構造変化後タンパク質座標データ選択ステップにて選択された上記構造変化後タンパク質座標データ、および、上記記憶手段に記憶された上記構造活性相関情報に基づいて、上記構造変換後タンパク質座標データ上に上記リガンドと重ね合わせを行う空間点を指定する空間点指定ステップと、
   上記空間点指定ステップにて指定された上記空間点と、上記記憶手段に記憶された上記リガンドの上記リガンド座標データとを用いて、上記タンパク質と上記リガンドとが結合した場合の相互作用関数を計算する際、上記タンパク質の上記動的挙動を表現する動的性質関数を「弾性エネルギー」として加えて計算する相互作用関数計算ステップと、
   を含み、
   上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドを評価し、選定することを特徴とするリガンド探索方法。
2. 請求項１に記載のリガンド探索方法において、上記相互作用関数計算ステップは、
   上記「弾性エネルギー」として、上記タンパク質の局所的な柔らかさを考慮した以下の数式１に示す上記動的性質関数「Ｕ衝突」を適応すること、
   を特徴とするリガンド探索方法。
   

   
3. 請求項１または２に記載のリガンド探索方法において、上記構造変化後タンパク質座標データ選択ステップは、
   上記基準振動計算より主鎖原子の２面角の上記ゆらぎ値を算出し、当該ゆらぎ値を以下の数式２または数式３に示す上記分子動力学計算における力の定数Ｋとすることにより、上記分子動力学計算を行うこと、
   を特徴とするリガンド探索方法。
   

   
4. 請求項１から３のいずれか１つに記載のリガンド探索方法において、上記相互作用関数計算ステップは、
   上記相互作用関数のスコアが最大になるように最適化する相互作用関数最適化ステップ、
   をさらに含むことを特徴とするリガンド探索方法。
5. 請求項４に記載のリガンド探索方法において、上記相互作用関数計算ステップは、
   上記相互作用関数最適化ステップにより上記相互作用関数を最適化した後に、重ねあわせた上記リガンドに対して、上記タンパク質との相互作用エネルギーを、以下の数式４に示す最適エネルギー関数「Ｕ最適」を用いて計算し、当該相互作用エネルギーについてリガンド立体構造データのコンフォメーションを微調整しながら最適化する相互作用エネルギー最適化ステップ、
   をさらに含むことを特徴とするリガンド探索方法。
   

   
6. 請求項５に記載のリガンド探索方法において、上記相互作用エネルギー最適化ステップは、上記最適エネルギー関数「Ｕ最適」に、以下の数式５に示す上記エネルギー関数「Ｕ _ＳＡＲ 」を適応すること、
   を特徴とするリガンド探索方法。
   

   
7. 請求項５または６に記載のリガンド探索方法において、上記相互作用エネルギー最適化ステップは、上記最適エネルギー関数「Ｕ最適」に、以下の数式６、数式７または数式８に示す上記エネルギー関数「Ｕ水素」を適応すること、
   を特徴とするリガンド探索方法。
   

   
8. 請求項５から７のいずれか１つに記載のリガンド探索方法において、上記相互作用エネルギー最適化ステップは、上記最適エネルギー関数「Ｕ最適」に、以下の数式９、および、数式１０に示す上記エネルギー関数「Ｕ疎水」を適応すること、
   を特徴とするリガンド探索方法。
   

   
9. 請求項５から８のいずれか１つに記載のリガンド探索方法において、上記相互作用エネルギー最適化ステップは、上記最適エネルギー関数「Ｕ最適」に、以下の数式１１、および、数式１２に示す上記エネルギー関数「Ｕスタッキング」を適応すること、
   を特徴とするリガンド探索方法。
   

   
10. 請求項５から９のいずれか１つに記載のリガンド探索方法において、上記相互作用エネルギー最適化ステップは、
    上記最適エネルギー関数「Ｕ最適」に、以下の数式１３に示す上記エネルギー関数「Ｕ内部」を適応すること、
    を特徴とするリガンド探索方法。
    

    
11. 請求項５から１０のいずれか１つに記載のリガンド探索方法において、
    上記処理手段において実行される、
    上記相互作用エネルギー最適化ステップにより最適化した後に、上記リガンド立体構造データのコンフォメーションを大きく変動させた後、再度、上記相互作用関数計算ステップを実行し、上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドの再評価を行う再評価ステップ、
    をさらに含むことを特徴とするリガンド探索方法。
12. 請求項１から１１のいずれか１つに記載のリガンド探索方法において、上記構造変化後タンパク質座標データ選択ステップは、
    上記誘導適合パラメータに、上記タンパク質の基準振動解析結果または二次構造判定結果を用いること、
    を特徴とするリガンド探索方法。

Images