JP2006252485A - リガンド探索装置、リガンド探索方法、プログラム、および記録媒体 - Google Patents

Abstract

【課題】本発明は、誘導適合型受容体も含めた受容体／リガンド結合の解析を研究するためのリガンド探索装置、リガンド探索方法、プログラム、および記録媒体を提供することを目的とする。
【解決手段】本発明は、受容体の基準振動解析計算をまず行い、定常状態の受容体の主鎖二面角の揺らぎを計算する。そして、その揺らぎを基にした拘束を各原子にかけながら分子動力学計算を行うことでより精度の良い受容体の動的構造を予測する。また、本発明は、分子動力学計算より得た動的構造及び相互作用関数を用いることで、誘導適合型受容体にも応用できる受容体／リガンド結合を精度良く予測する。
【選択図】 図１

Description

本発明は、タンパク質の立体構造座標を用いたリガンド探索装置、リガンド探索方法、プログラム、および記録媒体に関し、特に、タンパク質の立体構造座標が既知の場合、相互作用すると考えられるリガンドを予測するリガンド探索装置、リガンド探索方法、プログラム、および記録媒体に関するものである。

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

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

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

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

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

Ｍｏｒｒｉｓ， Ｇ． Ｍ． Ｇｏｏｄｓｅｌｌ， Ｄ． Ｓ． Ｈａｌｌｉｄａｙ， Ｒ． Ｓ． Ｈｕｅｙ， Ｒ． Ｈａｒｔ， Ｗ． Ｅ． Ｂｅｌｅｗ， Ｒ． Ｋ． Ｏｌｓｏｎ， Ａ． Ｊ． （１９９８） Ａｕｔｏｍａｔｅｄ ｄｏｃｋｉｎｇ ｕｓｉｎｇ ａ Ｌａｍａｒｃｋｉａｎ ｇｅｎｅｔｉｃ ａｌｇｏｒｉｔｈｍ ａｎｄ ａｎ ｅｍｐｉｒｉｃａｌ ｂｉｎｄｉｎｇ ｆｒｅｅ ｅｎｅｒｇｙ ｆｕｎｃｔｉｏｎ． Ｊ． Ｃｏｍｐｕｔ． Ｃｈｅｍ． １９： １６３９−１６６２； Ｇｏｏｄｓｅｌｌ， Ｄ． Ｓ． Ｍｏｒｒｉｓ， Ｇ． Ｍ． Ｏｌｓｏｎ， Ａ． Ｊ． （１９９６） Ａｕｔｏｍａｔｅｄ ｄｏｃｋｉｎｇ ｏｆ ｆｌｅｘｉｂｌｅ ｌｉｇａｎｄｓ： ａｐｐｌｉｃａｔｉｏｎｓ ｏｆ ＡｕｔｏＤｏｃｋ． Ｊ． Ｍｏｌ． Ｒｅｃｏｇｎｉｔ． ９： １−５ Ｅｗｉｎｇ， Ｔ． Ｊ． Ｍａｋｉｎｏ， Ｓ． Ｓｋｉｌｌｍａｎ， Ａ． Ｇ． Ｋｕｎｔｚ Ｉ． Ｄ． （２００１） ＤＯＣＫ ４．０： ｓｅａｒｃｈ ｓｔｒａｔｅｇｉｅｓ ｆｏｒ ａｕｔｏｍａｔｅｄ ｍｏｌｅｃｕｌａｒ ｄｏｃｋｉｎｇ ｏｆ ｆｌｅｘｉｂｌｅ ｍｏｌｅｃｕｌｅ ｄａｔａｂａｓｅｓ． Ｊ． Ｃｏｍｐｕｔ． Ａｉｄｅｄ Ｍｏｌ． Ｄｅｓ． １５： ４１１−２８ Ｒａｒｅｙ， Ｍ， Ｋｒａｍｅｒ， Ｂ， Ｌｅｎｇａｕｅｒ， Ｔ， Ｋｌｅｂｅ， Ｇ． （１９９６） Ａ ｆａｓｔ ｆｌｅｘｉｂｌｅ ｄｏｃｋｉｎｇ ｍｅｔｈｏｄ ｕｓｉｎｇ ａｎ ｉｎｃｒｅｍｅｎｔａｌ ｃｏｎｓｔｒｕｃｔｉｏｎ ａｌｇｏｒｉｔｈｍ． Ｊ． Ｍｏｌ． Ｂｉｏｌ． ２６１： ４７０−８９； Ｒａｒｅｙ， Ｍ． Ｗｅｆｉｎｇ， Ｓ． Ｌｅｎｇａｕｅｒ， Ｔ． （１９９６） Ｐｌａｃｅｍｅｎｔ ｏｆ ｍｅｄｉｕｍ−ｓｉｚｅｄ ｍｏｌｅｃｕｌａｒ ｆｒａｇｍｅｎｔｓ ｉｎｔｏ ａｃｔｉｖｅ ｓｉｔｅｓ ｏｆ ｐｒｏｔｅｉｎｓ． Ｊ． Ｃｏｍｐｕｔ． Ａｉｄｅｄ Ｍｏｌ． Ｄｅｓ． １０： ４１−５４ Ｊｏｎｅｓ， Ｇ． Ｗｉｌｌｅｔｔ， Ｐ． Ｇｌｅｎ， Ｒ． Ｃ． （１９９５） Ｍｏｌｅｃｕｌａｒ ｒｅｃｏｇｎｉｔｉｏｎ ｏｆ ｒｅｃｅｐｔｏｒ ｓｉｔｅｓ ｕｓｉｎｇ ａ ｇｅｎｅｔｉｃ ａｌｇｏｒｉｔｈｍ ｗｉｔｈ ａ ｄｅｓｃｒｉｐｔｉｏｎ ｏｆ ｄｅｓｏｌｖａｔｉｏｎ． Ｊ． Ｍｏｌ． Ｂｉｏｌ． ２４５： ４３−５３； Ｊｏｎｅｓ， Ｇ． Ｗｉｌｌｅｔｔ， Ｐ． Ｇｌｅｎ， Ｒ． Ｃ． Ｌｅａｃｈ， Ａ． Ｒ． Ｔａｙｌｏｒ， Ｒ． （１９９７） Ｄｅｖｅｌｏｐｍｅｎｔ ａｎｄ ｖａｌｉｄａｔｉｏｎ ｏｆ ａ ｇｅｎｅｔｉｃ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｆｌｅｘｉｂｌｅ ｄｏｃｋｉｎｇ． Ｊ． Ｍｏｌ． Ｂｉｏｌ． ２６７： ７２７−４８ Ｍｉｚｕｔａｎｉ， Ｍ． Ｙ． Ｉｔａｉ， Ａ． （２００４） Ｅｆｆｉｃｉｅｎｔ ｍｅｔｈｏｄ ｆｏｒ ｈｉｇｈ−ｔｈｒｏｕｇｈｐｕｔ ｖｉｒｔｕａｌ ｓｃｒｅｅｎｉｎｇ ｂａｓｅｄ ｏｎ ｆｌｅｘｉｂｌｅ ｄｏｃｋｉｎｇ： ｄｉｓｃｏｖｅｒｙ ｏｆ ｎｏｖｅｌ ａｃｅｔｙｌｃｈｏｌｉｎｅｓｔｅｒａｓｅ ｉｎｈｉｂｉｔｏｒｓ． Ｊ． Ｍｅｄ． Ｃｈｅｍ． ４７： ４８１８−４８２８

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

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

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

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

このような目的を達成するために、本発明にかかる請求項１に記載のリガンド探索装置は、単数または複数鎖のタンパク質の座標データが与えられた場合に、当該タンパク質と結合するリガンドを探索するリガンド探索装置において、上記タンパク質の座標データに対して、誘導適合を反映した誘導適合パラメータを用いて動的挙動を考慮した構造変化を行い、構造変化後タンパク質座標データを選択する構造変化後タンパク質座標データ選択手段と、上記構造変化後タンパク質座標データ選択手段にて選択された上記構造変化後タンパク質座標データから、上記リガンドと重ね合わせを行う空間点を指定する空間点指定手段と、上記空間点指定手段にて指定された上記空間点と、上記リガンドのリガンド座標データとを用いて、上記タンパク質と上記リガンドとが結合した場合の相互作用関数を計算する相互作用関数計算手段と、上記相互作用関数計算手段により計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドを評価するリガンド評価手段と、を備えたことを特徴とする。

また、本発明にかかる請求項２に記載のリガンド探索装置は、請求項１に記載のリガンド探索装置において、上記相互作用関数計算手段は、以下の数式１に示すＳｓｃｏｒｅ（ｉ，ｊ）を用いて上記相互作用関数を計算することを特徴とする。

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

また、本発明にかかる請求項３に記載のリガンド探索装置は、請求項１または２に記載のリガンド探索装置において、上記相互作用関数計算手段は、上記相互作用関数のスコアが最大になるように最適化する相互作用関数最適化手段、をさらに備えたことを特徴とする。

また、本発明にかかる請求項４に記載のリガンド探索装置は、請求項３に記載のリガンド探索装置において、上記相互作用関数計算手段は、上記相互作用関数最適化手段により上記相互作用関数を最適化した後に、重ねあわせた上記リガンドに対して、上記タンパク質との相互作用エネルギーを計算し、当該相互作用エネルギーについてリガンド立体構造データのコンフォメーションを微調整しながら最適化する相互作用エネルギー最適化手段、をさらに備えたことを特徴とする。

また、本発明にかかる請求項５に記載のリガンド探索装置は、請求項４に記載のリガンド探索装置において、上記リガンド評価手段は、上記相互作用エネルギー最適化手段により最適化した後に、上記リガンド立体構造データのコンフォメーションを大きく変動させた後、再度、上記相互作用関数計算手段を実行し、上記相互作用関数計算手段により計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドの再評価を行う再評価手段、をさらに備えたことを特徴とする。

また、本発明にかかる請求項６に記載のリガンド探索装置は、請求項１から５のいずれか１つに記載のリガンド探索装置において、上記構造変化後タンパク質座標データ選択手段は、上記誘導適合パラメータおよび／または上記構造変化後タンパク質座標データを算出するときに、上記タンパク質座標データに対し基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、当該ゆらぎの大きさを拘束条件として、分子動力学計算を行うこと、を特徴とする。

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

Ｅｒｏｔ＝Ｋｒｏｔ（φ―φ０）² ・・・（数式２）

（Ｅｒｏｔはタンパク質の立体構造中において主鎖原子の２面角のエネルギーを示す。φは主鎖原子の２面角である。φ０は主鎖原子の２面角の標準値である。ここで、Ｋｒｏｔの値が大きい場合は、φはφ０に拘束される。）

Ｅｐｏｓ＝Ｋｐｏｓ（ｒ−ｒ０）² ・・・（数式３）

（Ｅｐｏｓはタンパク質の立体構造中において主鎖原子の位置のエネルギーを示す。ｒは主鎖原子の座標である。ｒ０は主鎖原子の座標の標準値である。ここで、Ｋｐｏｓの値が大きい場合は、ｒはｒ０に拘束される。）

また、本発明にかかる請求項８に記載のリガンド探索装置は、請求項１から７のいずれか１つに記載のリガンド探索装置において、上記相互作用関数計算手段は、上記相互作用関数として、相互作用エネルギー関数にタンパク質の動的性質を表現する動的性質関数を「弾性エネルギー」として加えて用いること、を特徴とする。

また、本発明にかかる請求項９に記載のリガンド探索装置は、請求項８に記載のリガンド探索装置において、上記相互作用関数計算手段は、「弾性エネルギー」として、タンパク質の局所的な柔らかさを考慮し、以下の数式４に示す関数「Ｕ衝突」を適応すること、を特徴とする。

ψ（ｉ，ｊ）＝Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²
・・・（数式４）
（Ｍは衝突を不許可とする活性部位の原子の数であり、Ｎはリガンドの原子の数である。活性部位における動的挙動の少ない側鎖原子及び主鎖原子のみのｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義する。）

また、本発明にかかる請求項１０に記載のリガンド探索装置は、請求項１から７のいずれか１つに記載のリガンド探索装置において、上記相互作用関数計算手段は、上記相互作用関数として、相互作用エネルギー関数に、当該タンパク質の基準振動解析結果または二次構造判定結果をタンパク質の動的性質を表現する動的性質関数として加えて用いること、を特徴とする。

また、本発明はリガンド探索方法に関するものであり、本発明にかかる請求項１１に記載のリガンド探索方法は、単数または複数鎖のタンパク質の座標データが与えられた場合に、当該タンパク質と結合するリガンドを探索するリガンド探索方法において、上記タンパク質の座標データに対して、誘導適合を反映した誘導適合パラメータを用いて動的挙動を考慮した構造変化を行い、構造変化後タンパク質座標データを選択する構造変化後タンパク質座標データ選択ステップと、上記構造変化後タンパク質座標データ選択ステップにて選択された上記構造変化後タンパク質座標データから、上記リガンドと重ね合わせを行う空間点を指定する空間点指定ステップと、上記空間点指定ステップにて指定された上記空間点と、上記リガンドのリガンド座標データとを用いて、上記タンパク質と上記リガンドとが結合した場合の相互作用関数を計算する相互作用関数計算ステップと、上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドを評価するリガンド評価ステップと、を含むことを特徴とする。

また、本発明にかかる請求項１２に記載のリガンド探索方法は、請求項１１に記載のリガンド探索方法において、上記相互作用関数計算ステップは、以下の数式１に示すＳｓｃｏｒｅ（ｉ，ｊ）を用いて上記相互作用関数を計算することを特徴とする。

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

また、本発明にかかる請求項１３に記載のリガンド探索方法は、請求項１１または１２に記載のリガンド探索方法において、上記相互作用関数計算ステップは、上記相互作用関数のスコアが最大になるように最適化する相互作用関数最適化ステップ、をさらに含むことを特徴とする。

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

また、本発明にかかる請求項１５に記載のリガンド探索方法は、請求項１４に記載のリガンド探索方法において、上記リガンド評価ステップは、上記相互作用エネルギー最適化ステップにより最適化した後に、上記リガンド立体構造データのコンフォメーションを大きく変動させた後、再度、上記相互作用関数計算ステップを実行し、上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドの再評価を行う再評価ステップ、をさらに含むことを特徴とする。

また、本発明にかかる請求項１６に記載のリガンド探索方法は、請求項１１から１５のいずれか１つに記載のリガンド探索方法において、上記構造変化後タンパク質座標データ選択ステップは、上記誘導適合パラメータおよび／または上記構造変化後タンパク質座標データを算出するときに、上記タンパク質座標データに対し基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、当該ゆらぎの大きさを拘束条件として、分子動力学計算を行うこと、を特徴とする。

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

Ｅｒｏｔ＝Ｋｒｏｔ（φ―φ０）² ・・・（数式２）

（Ｅｒｏｔはタンパク質の立体構造中において主鎖原子の２面角のエネルギーを示す。φは主鎖原子の２面角である。φ０は主鎖原子の２面角の標準値である。ここで、Ｋｒｏｔの値が大きい場合は、φはφ０に拘束される。）

Ｅｐｏｓ＝Ｋｐｏｓ（ｒ−ｒ０）² ・・・（数式３）

（Ｅｐｏｓはタンパク質の立体構造中において主鎖原子の位置のエネルギーを示す。ｒは主鎖原子の座標である。ｒ０は主鎖原子の座標の標準値である。ここで、Ｋｐｏｓの値が大きい場合は、ｒはｒ０に拘束される。）

また、本発明にかかる請求項１８に記載のリガンド探索方法は、請求項１１から１７のいずれか１つに記載のリガンド探索方法において、上記相互作用関数計算ステップは、上記相互作用関数として、相互作用エネルギー関数にタンパク質の動的性質を表現する動的性質関数を「弾性エネルギー」として加えて用いること、を特徴とする。

また、本発明にかかる請求項１９に記載のリガンド探索方法は、請求項１８に記載のリガンド探索方法において、上記相互作用関数計算ステップは、「弾性エネルギー」として、タンパク質の局所的な柔らかさを考慮し、以下の数式４に示す関数「Ｕ衝突」を適応すること、を特徴とする。

ψ（ｉ，ｊ）＝Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²
・・・（数式４）
（Ｍは衝突を不許可とする活性部位の原子の数であり、Ｎはリガンドの原子の数である。活性部位における動的挙動の少ない側鎖原子及び主鎖原子のみのｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義する。）

また、本発明にかかる請求項２０に記載のリガンド探索方法は、請求項１１から１７のいずれか１つに記載のリガンド探索方法において、上記相互作用関数計算ステップは、上記相互作用関数として、相互作用エネルギー関数に、当該タンパク質の基準振動解析結果または二次構造判定結果をタンパク質の動的性質を表現する動的性質関数として加えて用いること、を特徴とする。

また、本発明はプログラムに関するものであり、本発明にかかる請求項２１に記載のリガンド探索方法をコンピュータに実行させることを特徴とするプログラムは、単数または複数鎖のタンパク質の座標データが与えられた場合に、当該タンパク質と結合するリガンドを探索するリガンド探索方法をコンピュータに実行させるプログラムにおいて、上記タンパク質の座標データに対して、誘導適合を反映した誘導適合パラメータを用いて動的挙動を考慮した構造変化を行い、構造変化後タンパク質座標データを選択する構造変化後タンパク質座標データ選択ステップと、上記構造変化後タンパク質座標データ選択ステップにて選択された上記構造変化後タンパク質座標データから、上記リガンドと重ね合わせを行う空間点を指定する空間点指定ステップと、上記空間点指定ステップにて指定された上記空間点と、上記リガンドのリガンド座標データとを用いて、上記タンパク質と上記リガンドとが結合した場合の相互作用関数を計算する相互作用関数計算ステップと、上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドを評価するリガンド評価ステップと、を含むことを特徴とする。

また、本発明にかかる請求項２２に記載のプログラムは、請求項２１に記載のプログラムにおいて、上記相互作用関数計算ステップは、以下の数式１に示すＳｓｃｏｒｅ（ｉ，ｊ）を用いて上記相互作用関数を計算することを特徴とする。

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

また、本発明にかかる請求項２３に記載のプログラムは、請求項２１または２２に記載のプログラムにおいて、上記相互作用関数計算ステップは、上記相互作用関数のスコアが最大になるように最適化する相互作用関数最適化ステップ、をさらに含むことを特徴とする。

また、本発明にかかる請求項２４に記載のプログラムは、請求項２３に記載のプログラムにおいて、上記相互作用関数計算ステップは、上記相互作用関数最適化ステップにより上記相互作用関数を最適化した後に、重ねあわせた上記リガンドに対して、上記タンパク質との相互作用エネルギーを計算し、当該相互作用エネルギーについてリガンド立体構造データのコンフォメーションを微調整しながら最適化する相互作用エネルギー最適化ステップ、をさらに含むことを特徴とする。

また、本発明にかかる請求項２５に記載のプログラムは、請求項２４に記載のプログラムにおいて、上記リガンド評価ステップは、上記相互作用エネルギー最適化ステップにより最適化した後に、上記リガンド立体構造データのコンフォメーションを大きく変動させた後、再度、上記相互作用関数計算ステップを実行し、上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドの再評価を行う再評価ステップ、をさらに含むことを特徴とする。

また、本発明にかかる請求項２６に記載のプログラムは、請求項２１から２５のいずれか１つに記載のプログラムにおいて、上記構造変化後タンパク質座標データ選択ステップは、上記誘導適合パラメータおよび／または上記構造変化後タンパク質座標データを算出するときに、上記タンパク質座標データに対し基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、当該ゆらぎの大きさを拘束条件として、分子動力学計算を行うこと、を特徴とする。

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

Ｅｒｏｔ＝Ｋｒｏｔ（φ―φ０）² ・・・（数式２）

（Ｅｒｏｔはタンパク質の立体構造中において主鎖原子の２面角のエネルギーを示す。φは主鎖原子の２面角である。φ０は主鎖原子の２面角の標準値である。ここで、Ｋｒｏｔの値が大きい場合は、φはφ０に拘束される。）

Ｅｐｏｓ＝Ｋｐｏｓ（ｒ−ｒ０）² ・・・（数式３）

（Ｅｐｏｓはタンパク質の立体構造中において主鎖原子の位置のエネルギーを示す。ｒは主鎖原子の座標である。ｒ０は主鎖原子の座標の標準値である。ここで、Ｋｐｏｓの値が大きい場合は、ｒはｒ０に拘束される。）

また、本発明にかかる請求項２８に記載のプログラムは、請求項２１から２７のいずれか１つに記載のプログラムにおいて、上記相互作用関数計算ステップは、上記相互作用関数として、相互作用エネルギー関数にタンパク質の動的性質を表現する動的性質関数を「弾性エネルギー」として加えて用いること、を特徴とする。

また、本発明にかかる請求項２９に記載のプログラムは、請求項２８に記載のプログラムにおいて、上記相互作用関数計算ステップは、「弾性エネルギー」として、タンパク質の局所的な柔らかさを考慮し、以下の数式４に示す関数「Ｕ衝突」を適応すること、を特徴とする。

ψ（ｉ，ｊ）＝Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²
・・・（数式４）
（Ｍは衝突を不許可とする活性部位の原子の数であり、Ｎはリガンドの原子の数である。活性部位における動的挙動の少ない側鎖原子及び主鎖原子のみのｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義する。）

また、本発明にかかる請求項３０に記載のプログラムは、請求項２１から２７のいずれか１つに記載のプログラムにおいて、上記相互作用関数計算ステップは、上記相互作用関数として、相互作用エネルギー関数に、当該タンパク質の基準振動解析結果または二次構造判定結果をタンパク質の動的性質を表現する動的性質関数として加えて用いること、を特徴とする。

また、本発明は記録媒体に関するものであり、本発明にかかる請求項３１に記載の記録媒体は、請求項２１から３０のいずれか１つに記載されたプログラムを記録したことを特徴とする。

本発明によれば、（１）タンパク質の座標データに対して、誘導適合を反映した誘導適合パラメータを用いて動的挙動を考慮した構造変化を行い、構造変化後タンパク質座標データを選択し、（２）選択された構造変化後タンパク質座標データから、リガンドと重ね合わせを行う空間点を指定し、（３）指定された空間点と、リガンドのリガンド座標データとを用いて、タンパク質とリガンドとが結合した場合の相互作用関数を計算し、（４）計算された相互作用関数に基づいて当該タンパク質と結合するリガンドを評価する。これにより、リガンドのフレキシビリティだけでなく受容体側のフレキシビリティも考慮して、Ｉｎｄｕｃｅｄ−Ｆｉｔ型の受容体タンパク質に結合するリガンドを効率的に精度よく探索することができるという効果を奏する。

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

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

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

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

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

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

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

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

以下に、本発明にかかるリガンド探索装置、リガンド探索方法、プログラム、および記録媒体の実施の形態を図面に基づいて詳細に説明する。なお、この実施の形態により本発明が限定されるものではない。

本明細書において幾つかの用語を使用するが、特に明記しない限り、次の意味を有する。
「標的タンパク質」とは、立体構造の詳細がＸ線結晶解析やＮＭＲ解析、ホモロジーモデリング法により既に決定されており、リガンド探索の対象とするタンパク質を意味する。
「原子座標」とは、三次元空間上で立体構造を記述するものである。それは空間上のある点を原点とする互いに垂直な三方向の相対的な距離であり、タンパク質中に存在する水素原子を除く原子１つあたりに３個の数字からなるベクトル量である。

［本発明の基本原理］
ここでは、本発明の基本原理について、図１を参照して説明する。図１は、本発明の基本原理を示す概念図である。本発明は、単数または複数鎖のタンパク質立体構造が与えられた場合において、当該タンパク質の立体構造から誘導適合を反映したパラメータおよび構造変化した立体構造座標を例えば基準振動計算方法や分子動力学計算方法により予め算出し、当該パラメータおよび構造変化した立体構造座標を用いて該当タンパク質と別の物質（リガンド）が結合した場合の相互作用関数を定義し、当該相互作用関数によって該当タンパク質と結合する物質（リガンド）を評価し、選定するリガンド探索装置、リガンド探索方法、コンピュータプログラム、および記録媒体に関する。

まず、本発明は、化合物データベースからリガンドを１つ選択し、そのリガンドの立体構造データを取得する（ステップＳ−１）。また、対象となるタンパク質の立体構造データを取得する（ステップＳ−２）。

つぎに、本発明は、当該タンパクの立体構造データに基づいて、誘導適合を反映する誘導適合パラメータを用いて動的挙動を考慮した構造変化を行い、複数の構造変化後タンパク質座標データを用意し、ランダムに１つの構造変化後タンパク質座標データを選択する（ステップＳ−３）。

つぎに、本発明は、リガンドと重ね合わせを行う、構造変化後タンパク質座標データ中の空間点を指定する（ステップＳ−４）。ここで、当該空間点は、例えば、以下に示す（１）や（２）などの方法で指定してもよい。

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

そこで、本発明では、活性部位及びリガンドに水素原子が含まれていない場合、以下の規則によりダミー水素原子を発生させる。
１）ｓｐ²軌道原子を中心とする正三角形状にダミー原子を発生させる（図２）。すなわち、図２に示すように、ｓｐ²軌道原子の窒素原子（Ａ）を中心とする正三角形の空いている位置にダミー水素原子（Ｂ）を発生させる。
２）ｓｐ³軌道原子では、水素結合を形成する距離にある場合、水素原子を共有するように回転できると考え、水素結合相互作用を計算するときには距離のみを考慮する。このため、ｓｐ³軌道原子にはダミー原子を発生させない。

また、金属及び水では、活性部位・リガンド結合の仲介役となり得るので、相互作用する位置に以下のようにダミー原子を発生させる。
１）鉄のような金属には、正八面体状にダミー原子を発生させる（図３）。すなわち、図３に示すように、亜鉛（Ａ）を中心とする正八面体の空いている位置にダミー原子（Ｂ）を発生させる。
２）水には、正四面体状にダミー原子を発生させる。

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

（２）構造活性相関情報を利用した空間点の指定（ステップＳ−４−２）
また、本発明は、リガンドの構造活性相関（ＳＡＲ）情報に着目し、以下の（Ａ）〜（Ｄ）の項目を入力情報にすることにより、空間点を指定する。
（Ａ）ＳＡＲ情報から得られた活性部位の原子（以後、「Ａ原子」）。ＰＤＢ形式に従う。
（Ｂ）「Ａ原子」と相互作用するであろうリガンドの原子タイプ（以後、「Ｂタイプ」）。ＳＹＢＹＬのＭＯＬ２形式に従う。
（Ｃ）「Ａ原子」と「Ｂタイプ」との相互作用の強さ（以後、「Ｃ強さ」）。
（Ｄ）「Ａ原子」と「Ｂタイプ」との相互作用する距離（以後、「Ｄ距離」）（単位はÅ）。

なお、本発明では、上記（Ａ）〜（Ｄ）の入力情報に基づいて、タンパク質中の活性部位内におけるリガンドの初期座標を利用し、以下の１）〜４）に示す規則により空間点を作成してもよい。
１）「Ａ原子」がドナーまたは金属及び水の場合（ＳＡＲ情報の活性部位側の指定が水素結合ドナー、金属原子の場合）には、ステップＳ−４−１で発生させたダミー原子の方向に対して「Ａ原子」から「Ｄ距離」の位置及びその周囲を初期座標に選ぶ（図４、図５）。
２）「Ａ原子」がｓｐ³軌道原子の場合（ＳＡＲ情報の活性部位側の指定がｓｐ³軌道原子の場合）には、「Ａ原子」から「Ｄ距離」の周囲を初期座標に選ぶ（図６）。
３）「Ａ原子」が水素結合アクセプター（以後、アクセプター）の場合（ＳＡＲ情報の活性部位側の指定が水素結合アクセプターの場合）には、「Ａ原子」の結合延長上「Ｄ距離」の位置及びその周囲を初期座標に選ぶ（図７）。
４）その他の場合（ＳＡＲ情報の活性部位側の指定がその他の原子の場合）には、「Ａ原子」を中心とする半径が「Ｄ距離」の球表面上の点を初期座標に選ぶ（図８）。

ここで、上記１）〜４）とは異なり、リガンドの初期座標を直接指定してもよい。

再び図１に戻り、本発明は、ステップＳ−１において選択したリガンド座標データ中の原子と、ステップＳ−４で指定したタンパク質座標データ中の空間点とのペアを重複がないようにランダムに選択する（ステップＳ−５）。

つぎに、本発明は、以下の数式１に示す相互作用関数であるスコアＳｓｃｏｒｅ（ｉ，ｊ）を計算する（ステップＳ−６）。

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

つぎに、本発明は、ステップＳ−６により求めた相互作用関数のスコアが最大になるように調整（最適化）する（ステップＳ−７）。ここで、スコアを最大にする手法としては、例えば、シミュレーティッドアニーリング法が挙げられる。また、時間短縮にはステップＳ−５およびステップＳ−６を複数回（例えば１００００回）繰り返し、Ｓｓｃｏｒｅ（ｉ，ｊ）が最大になるペアを探し、そのペア情報をもとにリガンドを初期座標に重ね合わせる方法を適応することが好ましい。

つぎに、本発明は、ステップＳ−７により相互作用関数を最適化したときに、重ね合わせたリガンドに対して、タンパク質との相互作用エネルギーを計算し、当該相互作用エネルギーについてリガンド座標データのコンフォメーションを微調整しながら最適化する（ステップＳ−８）。リガンドのコンフォメーションの微調整は、ステップＳ−７で算出されたリガンド座標を中心に、例えば、並進、回転、または、シングルボンドまわりの角度をＲＳＭＤで０．３を越えない程度に座標変化などをさせることで、行ってもよい。

ここで、リガンド座標データのコンフォメーションの微調整は例えばランダムサーチで最適化することが好ましい。なお、ランダムサーチでは、以下の１）〜３）の項目に従って、タンパク質の活性部位とリガンドとの微小変化を例えば８０００回行い、最適エネルギー「Ｕ最適」が最小になるようにする。
１）回転可能な結合のうち最大で５つを乱数で選び、結合ごとにランダムに±１０．０°の範囲内で回転させ、リガンドのコンフォメーションを換える。この過程を例えば３回に一度行う。
２）ｘ、ｙ、ｚ軸方向それぞれにおいて、ランダムに±１．０Åの範囲内でリガンドの並進運動を行う。この過程を例えば２回に一度行う。
３）回転中心座標それぞれにおいて、ランダムに±１．０Åの範囲内で回転中心座標を移動させ、さらに３次元方向の角度それぞれに対して、ランダムに±５．０°の範囲内でリガンドの回転運動を行う。この過程を例えば５回に一度行う。

つぎに、本発明は、リガンド座標データのコンフォメーションを大きく変動して、ステップＳ−５から再スタートを行い、ステップＳ−８までを繰り返して最適化を繰り返し行う（ステップＳ−９）。なお、コンフォメーション改変は、ステップＳ−７で算出されたリガンド座標を中心に、例えば、並進、回転、または、シングルボンドまわりの角度をＲＳＭＤで０．３以上になるよう座標変化などをさせることで、行ってもよい。

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

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

つぎに、本発明は、ステップＳ−４〜ステップＳ−９までの過程を、ステップＳ−３で用意した複数の構造変化後タンパク質座標データに対して行い、最適なタンパク質とリガンドとの複合体座標、および最適エネルギー「Ｕ最適」を算出する（ステップＳ−１０）。

つぎに、本発明は、以上の過程を、ステップＳ−１で用意した化合物データベース中の全てのリガンドに対して行い、化合物データベース中から該当タンパク質と結合する可能性のあるリガンドを選択する（ステップＳ−１１）。

以上、本発明の基本原理について説明したが、本発明は、タンパク質の誘導適合を反映するパラメータおよび／または構造変化した立体構造座標を分子動力学計算方法を用いて算出する場合、該当タンパク質の立体構造に対し基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、そのゆらぎの大きさを拘束条件として、分子動力学計算を行うことで、タンパク質の立体構造をエネルギー最適構造より大きく離れないようにして分子動力学計算を行ってもよい。

また、本発明において、本分子動力学計算方法による分子動力学計算は、例えば、基準振動計算より主鎖原子の２面角のゆらぎの値を算出し、当該ゆらぎ値を以下の数式２または数式３のように分子動力学計算における力の定数Ｋの部分に入れてもよい。

Ｅｒｏｔ＝Ｋｒｏｔ（φ―φ０）² ・・・（数式２）

Ｅｒｏｔはタンパク質の立体構造中において主鎖原子の２面角のエネルギーを示す。φは主鎖原子の２面角である。φ０は主鎖原子の２面角の標準値である。ここで、Ｋｒｏｔの値が大きい場合は、φはφ０に拘束される。

Ｅｐｏｓ＝Ｋｐｏｓ（ｒ−ｒ０）² ・・・（数式３）

Ｅｐｏｓはタンパク質の立体構造中において主鎖原子の位置のエネルギーを示す。ｒは主鎖原子の座標である。ｒ０は主鎖原子の座標の標準値である。ここで、Ｋｐｏｓの値が大きい場合は、ｒはｒ０に拘束される。

また、本発明は、リガンドとタンパク質との相互作用を評価する際の目的関数（相互作用関数）として、従来の相互作用エネルギー関数に、タンパク質の動的性質を表現する動的性質関数を「弾性エネルギー」として加えてもよい。これにより、タンパク質の立体構造座標から相互作用エネルギーを高速に算出するとともに、タンパク質の動的挙動に関する物理化学的性質を明確に描写することができる。

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

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

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

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

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

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

また、本発明は、タンパク質の動的性質関数としてタンパク質の基準振動解析結果または二次構造判定結果を用いてもよい。なお、二次構造判定においては、タンパク質のへリックスやシート部分では揺らぎは小さいと考え、それ以外では揺らぎは大きいと考え、相互作用の評価関数、分子動力学計算の拘束条件に適応する。

また、本発明によれば、計算された該当タンパク質（タンパク質座標データ）が代表的な複数の立体構造座標である場合や、該当タンパク質の立体構造が例えば核磁気共鳴スペクトルの解析結果のような複数の立体構造座標である場合についても、該当タンパク質と結合するリガンドの探索について、複数座標すべてを全自動的にかつ短時間で同等に評価することができる。

［システム構成］
ここでは、本発明が適用される本システムの構成について、図１１６を参照して詳細に説明する。図１１６は、本発明が適用される本システムの構成の一例を示すブロック図であり、該構成のうち本発明に関係する部分のみを概念的に示している。

図１１６に示すように、本システムは、概略的に、タンパク質と結合する物質（リガンド）を評価し、選定するリガンド探索装置１００と、リガンド立体構造データやタンパク質立体構造データなどに関する外部データベースや各種の外部プログラムなどを提供する外部システム２００とを、ネットワーク３００を介して通信可能に接続して構成されている。

ネットワーク３００は、リガンド探索装置１００と外部システム２００とを相互に接続する機能を有し、例えばインターネットやＬＡＮなどである。

外部システム２００は、ネットワーク３００を介して、リガンド探索装置１００と相互に接続され、利用者に対してリガンド立体構造データやタンパク質立体構造データなどに関する外部データベースや各種の外部プログラムを実行するウェブサイトを提供する機能を有する。ここで、外部システム２００は、ＷＥＢサーバやＡＳＰサーバ等として構成してもよく、そのハードウェア構成は、一般に市販されるワークステーション、パーソナルコンピュータ等の情報処理装置およびその付属装置により構成してもよい。また、外部システム２００の各機能は、外部システム２００のハードウェア構成中のＣＰＵ、ディスク装置、メモリ装置、入力装置、出力装置、通信制御装置等、およびそれらを制御するプログラム等により実現される。

リガンド探索装置１００は、概略的に、リガンド探索装置１００の全体を統括的に制御するＣＰＵ等の制御部１０２と、通信回線等に接続されるルータ等の通信装置（図示せず）に接続される通信制御インターフェース部１０４と、各種のデータベースやファイルなどを格納する記憶部１０６と、入力装置１１２や出力装置１１４に接続される入出力制御インターフェース部１０８と、を備えて構成されており、これら各部は任意の通信路を介して通信可能に接続されている。さらに、リガンド探索装置１００は、ルータ等の通信装置および専用線等の有線または無線の通信回線を介して、ネットワーク３００に通信可能に接続されている。

記憶部１０６に格納される各種のデータベースやテーブルやファイル（リガンド座標データベース１０６ａ〜リガンド評価結果ファイル１０６ｆ）は、固定ディスク装置等のストレージ手段であり、各種処理に用いる各種のプログラムやテーブルやファイルやデータベースやウェブページ用ファイルなどを格納する。

これら記憶部１０６の各構成要素のうち、リガンド座標データベース１０６ａは、リガンド座標データを格納するリガンド座標データ格納手段である。タンパク質座標データベース１０６ｂは、タンパク質座標データを格納するタンパク質座標データ格納手段である。構造変化後タンパク質座標データファイル１０６ｃは、後述する構造変化後タンパク質座標データ選択部１０２ａにより選択された構造変化後タンパク質座標データを格納する構造変化後タンパク質座標データ格納手段である。指定空間点ファイル１０６ｄは、後述する空間点指定部１０２ｂにより指定された空間点に関する情報を格納する指定空間点格納手段である。相互作用関数計算結果ファイル１０６ｅは、後述する相互作用関数計算部１０２ｃにより計算された相互作用関数の計算結果に関する情報を格納する相互作用関数計算結果格納手段である。リガンド評価結果ファイル１０６ｆは、後述するリガンド評価部１０２ｄにより評価されたリガンドの評価結果に関する情報を格納するリガンド評価結果格納手段である。

通信制御インターフェース部１０４は、リガンド探索装置１００とネットワーク３００（またはルータ等の通信装置）との間における通信制御を行う。すなわち、通信制御インターフェース部１０４は、他の端末と通信回線を介してデータを通信する機能を有する。

入出力制御インターフェース部１０８は、入力装置１１２や出力装置１１４の制御を行う。ここで、出力装置１１４としては、モニタ（家庭用テレビを含む）の他、スピーカ等を用いることができる（なお、以下においては出力装置１１４をモニタとして記載する場合がある。）。また、入力装置１１２としては、キーボードやマウス、マイクなどを用いることができる。また、モニタも、マウスと協働してポインティングデバイス機能を実現する。

制御部１０２は、ＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）等の制御プログラム、および所要データを格納するための内部メモリを有し、これらのプログラム等により種種の処理を実行するための情報処理を行う。制御部１０２は、機能概念的に、構造変化後タンパク質座標データ選択部１０２ａと、空間点指定部１０２ｂと、相互作用関数計算部１０２ｃと、リガンド評価部１０２ｄと、を含んで構成されている。

これら制御部１０２の各構成要素のうち、構造変化後タンパク質座標データ選択部１０２ａは、タンパク質の座標データに対して、誘導適合を反映した誘導適合パラメータを用いて動的挙動を考慮した構造変化を行い、構造変化後タンパク質座標データを選択する構造変化後タンパク質座標データ選択手段である。空間点指定部１０２ｂは、構造変化後タンパク質座標データ選択部１０２ａにより選択された構造変化後タンパク質座標データから、リガンドと重ね合わせを行う空間点を指定する空間点指定手段である。

相互作用関数計算部１０２ｃは、空間点指定部１０２ｂにより指定された空間点とリガンドのリガンド座標データとを用いて、タンパク質とリガンドとが結合した場合の相互作用関数を計算する相互作用関数計算手段である。ここで、相互作用関数計算部１０２ｃは、図１１７に示すように、相互作用関数最適化部１０２ｃ１と、相互作用エネルギー最適化部１０２ｃ２と、をさらに含んで構成されている。相互作用関数最適化部１０２ｃ１は、相互作用関数のスコアが最大になるように最適化する相互作用関数最適化手段である。相互作用エネルギー最適化部１０２ｃ２は、相互作用関数最適化部１０２ｃ１により相互作用関数を最適化した後に、重ねあわせたリガンドに対してタンパク質との相互作用エネルギーを計算し、当該相互作用エネルギーについてリガンド立体構造データのコンフォメーションを微調整しながら最適化する相互作用エネルギー最適化手段である。

再び図１１６に戻り、リガンド評価部１０２ｄは、相互作用関数計算部１０２ｃにより計算された相互作用関数に基づいて当該タンパク質と結合するリガンドを評価するリガンド評価手段である。ここで、リガンド評価部１０２ｄは、図１１８に示すように、再評価部１０２ｄ１をさらに含んで構成されている。再評価部１０２ｄ１は、相互作用エネルギー最適化部１０２ｃ２により最適化した後に、リガンド立体構造データのコンフォメーションを大きく変動させた後、再度、相互作用関数計算部１０２ｃを実行し、相互作用関数計算部１０２ｃにより計算された相互作用関数に基づいて当該タンパク質と結合するリガンドの再評価を行う再評価手段である。

なお、これら各部によって行われる処理の詳細については、後述する。

［システムの処理］
ここでは、上述のように構成された本実施の形態における本システムのメイン処理の一例について、図１１５などを参照して詳細に説明する。図１１５は、本実施形態における本システムのメイン処理の一例を示すフローチャートである。ここでは、図１１５を参照して、タンパク質立体構造と誘導適合を利用したリガンドの探索に関して説明する。

まず、リガンド探索装置１００は、制御部１０２の処理により、３次元座標を含むリガンドのデータベースを用意し、用意したデータベースを記憶部１０６のリガンド座標データベース１０６ａに格納する（ステップＳ０）。ここで、リガンドのデータベースとしては、例えば、ＡＣＤ等のような市販化合物データベースや、化合物を描いて収集した仮想化合物データなどを用いてもよい。なお、リガンドのデータベースは、分子力学法等を用いて３次元化することが望ましい。

つぎに、リガンド探索装置１００は、制御部１０２の処理により、ステップＳ０で用意したリガンドデータベースから、特定リガンドを探索するための標的タンパク質の立体構造を選択し、当該選択したタンパク質の立体構造データ（３次元座標）を入手し、記憶部１０６のタンパク質座標データベース１０６ｂに格納する（ステップＳ１０）。なお、３次元座標には、公共データベースであるＰＤＢやホモロジーモデリング法等で作成した立体構造座標を用いることが望ましい。

つぎに、リガンド探索装置１００は、構造変化後タンパク質座標データ選択部１０２ａの処理により、ステップＳ１０により選択された標的タンパク質の基準振動計算を行い、主鎖原子の位置のゆらぎと２面角のゆらぎを求める（ステップＳ２０）。具体的には、まず、ステップＳ１０で定められた標的タンパク質の動的挙動を表すパラメータを基準振動解析法による計算結果のデータベースから取得する、もしくは当該パラメータを二次構造判定計算を行って取得する。

まず、ステップＳ２０において、タンパク質の動的挙動を表すパラメータを基準振動解析法により取得する方法について説明する。基準振動解析法とは、ポテンシャルエネルギーを変位の二次関数として近似し、運動方程式を厳密に解き、最適化構造の周りの微小な振動を解析する方法である。解くべき運動方程式は下記の式（１）または式（２）である。なお、基準振動解析法の詳細については、「Ｗｉｌｓｏｎ，Ｅ．Ｂ．，Ｄｅｃｉｕｓ，Ｊ．Ｃ．，ａｎｄ Ｃｒｏｓｓ，Ｐ． Ｃ． １９５５． Ｍｏｌｅｃｕｌａｒ Ｖｉｂｒａｔｉｏｎ． ＭｃＧｒａｗ−Ｈｉｌｌ．」に記載されている。

・・・（１）

・・・（２）

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

・・・（３）

・・・（４）
ここで、ｑ_iは振動の自由度に対応した座標である。また、ｑ_i’（式（
３）における「ｑｉドット」を意味する）はｑｉの時間による微分である。また、ｑ_jは、下記の式（５）の通りである。

・・・（５）
ここで、Ａ_jkは、集団運動Ｑ_kと個々の原子運動ｑ_jとを結ぶ係数である。また、ｑ_j ⁰は最適化座標である。ただし、Ｑ_kは、下記の式に示す基準振動である。

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

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

また、ステップＳ２０において、タンパク質の動的挙動を表すパラメータを二次構造判定計算を行って取得する方法について説明する。二次構造判定はタンパク質の立体構造座標から計算される。なお、ソフトウェアとしては、ＤＳＳＰ、ＳＴＲＩＤＥ等が好ましいが、基本的にはタンパク質の主鎖のねじれ角と水素結合パターンから判別される方法を用いればよい。ここで、「ＤＳＳＰ（Ｄｉｃｔｉｏｎａｒｙ ｏｆ ｐｒｏｔｅｉｎ ｓｅｃｏｎｄａｒｙ ｓｔｒｕｃｔｕｒｅ ｏｆ ｐｒｏｔｅｉｎ）」とは、ＰＤＢ書式のファイルを入力ファイルとして、主鎖の水素結合パターンと、内部回転角等を解析しαへリックスとβシートとを判定するソフトウェアである。ＤＳＳＰの詳細は、「Ｋａｂｓｃｈ，Ｗ． ＆ Ｓａｎｄｅｒ，Ｃ． （１９８３） Ｄｉｃｔｉｏｎａｒｙ ｏｆ ｐｒｏｔｅｉｎ ｓｅｃｏｎｄａｒｙ ｓｔｒｕｃｔｕｒｅ：ｐａｔｔｅｒｎ ｒｅｃｏｇｎｉｔｉｏｎ ｏｆ ｈｙｄｒｏｇｅｎ−ｂｏｎｄｅｄ ａｎｄ ｇｅｏｍｅｔｒｉｃａｌ ｆｅａｔｕｒｅｓ． Ｂｉｏｐｏｌｉｍｅｒｓ， ２２：２５７７−２６３７」に記載されている。また、「ＳＴＲＩＤＥ（Ｐｒｏｔｅｉｎ ｓｅｃｏｎｄａｒｙ ｓｔｒｕｃｔｕｒｅ ａｓｓｉｇｎｍｅｎｔ ｆｒｏｍ ａｔｏｍｉｃ ｃｏｏｒｄｉｎａｔｅ）」とは、ＰＤＢ書式のファイルを入力ファイルとして、主鎖の水素結合パターンと、内部回転角等を解析しαへリックスとβシートとを判定するソフトウェアである。ＳＴＲＩＤＥの詳細は、「Ｆｒｉｓｈｍａｎ， Ｄ ＆ Ａｒｇｏｓ，Ｐ．（１９９５） Ｋｎｏｗｌｅｄｇｅ−ｂａｓｅｄ ｓｅｃｏｎｄａｒｙ ｓｔｒｕｃｔｕｒｅ ａｓｓｉｇｎｍｅｎｔ． Ｐｒｏｔｅｉｎｓ： ｓｔｒｕｃｔｕｒｅ， ｆｕｎｃｔｉｏｎ ａｎｄ ｇｅｎｅｔｉｃｓ， ２３， ５６６−５７９」に記載されている。

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

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

再び図１１５のメイン処理の説明に戻り、リガンド探索装置１００は、構造変化後タンパク質座標データ選択部１０２ａの処理により、ステップＳ２０において求めた標的タンパク質のゆらぎを拘束条件として用いた分子動力学計算を行う（ステップＳ３０）。

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

Ｕ位置 ＝ Ｋ位置 ＊ Ｒ²
・・・（６）
ここで、「Ｋ位置」は例えば３００．０とし、Ｒは基準座標からのずれとする。

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

Ｕ二面角 ＝ Ｋ二面角 ＊ （θ−θ平衡）²
・・・（７）
θは二面角（単位ｒａｄ）である。「Ｋ二面角」には、最大値と最小値を指定することで、その範囲内で主鎖二面角揺らぎに対応するように各二面角に対して不均一な拘束がかかるようにする。以後、主鎖二面角を拘束しながら行うＭＤを二面角拘束ＭＤと呼ぶことにする。

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

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

再び図１１５のメイン処理の説明に戻り、リガンド探索装置１００は、空間点指定部１０２ｂの処理により、標的タンパク質のリガンド結合部位に、リガンドを配置するための空間点群を指定する（ステップＳ４０）。具体的には、ステップＳ３０で作成した複数のタンパク質立体構造座標のうち、１つをランダムに選択する。ここで、タンパク質座標中の空間点は、例えば以下の（１）または（２）などの方法で指定させる。

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

そこで、本実施形態では、活性部位及びリガンドに水素原子が含まれていない場合、以下の規則によりダミー水素原子を発生させる。
１）ｓｐ²軌道原子を中心とする正三角形状にダミー原子を発生させる（図２）。すなわち、図２に示すように、ｓｐ²軌道原子の窒素原子（Ａ）を中心とする正三角形の空いている位置にダミー水素原子（Ｂ）を発生させる。
２）ｓｐ³軌道原子では、水素結合を形成する距離にある場合、水素原子を共有するように回転できると考え、水素結合相互作用を計算するときには、距離のみを考慮する。このため、ｓｐ³軌道原子にはダミー原子を発生させない。

また、金属及び水では、活性部位・リガンド結合の仲介役となり得るので、相互作用する位置に以下の規則によりダミー原子を発生させる。
１）鉄のような金属には、正八面体状にダミー原子を発生させる（図３）。すなわち、図３に示すように、亜鉛（Ａ）を中心とする正八面体の空いている位置にダミー原子（Ｂ）を発生させる。
２）水には、正四面体状にダミー原子を発生させる。

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

（２）構造活性相関情報を利用した空間点の指定
また、本実施形態では、リガンドの構造活性相関（ＳＡＲ）情報に着目し、以下の（Ａ）〜（Ｄ）項目を入力情報にすることにより、空間点を指定する。
（Ａ）ＳＡＲ情報から得られた活性部位の原子（以後、「Ａ原子」）。ＰＤＢ形式に従う。
（Ｂ）「Ａ原子」と相互作用するであろうリガンドの原子タイプ（以後、「Ｂタイプ」）。ＳＹＢＹＬのＭＯＬ２形式に従う。
（Ｃ）「Ａ原子」と「Ｂタイプ」との相互作用の強さ（以後、「Ｃ強さ」）。
（Ｄ）「Ａ原子」と「Ｂタイプ」との相互作用する距離（以後、「Ｄ距離」）（単位はÅ）。

なお、本実施形態では、（Ａ）〜（Ｄ）の入力情報に基づいてタンパク質中の活性部位内におけるリガンドの初期座標を利用し、以下の１）〜４）に示す規則により空間点を作成してもよい。
１）「Ａ原子」がドナーまたは金属及び水の場合（ＳＡＲ情報の活性部位側の指定が水素結合ドナー、金属原子の場合）には、「（１）ダミー原子の発生による空間点の指定」で発生させたダミー原子の方向に対して「Ａ原子」から「Ｄ距離」の位置及びその周囲を初期座標に選ぶ（図４、図５）。
２）「Ａ原子」がｓｐ³軌道原子の場合（ＳＡＲ情報の活性部位側の指定がｓｐ³軌道原子の場合）には、「Ａ原子」から「Ｄ距離」の周囲を初期座標に選ぶ（図６）。
３）「Ａ原子」が水素結合アクセプター（以後、アクセプター）の場合（ＳＡＲ情報の活性部位側の指定が水素結合アクセプターの場合）には、「Ａ原子」の結合延長上「Ｄ距離」の位置及びその周囲を初期座標に選ぶ（図７）。
４）その他の場合（ＳＡＲ情報の活性部位側の指定がその他の原子の場合）には、「Ａ原子」を中心とする半径が「Ｄ距離」の球表面上の点を初期座標に選ぶ（図８）。

ここで、上記１）〜４）とは異なり、リガンドの初期座標を直接指定してもよい。

再び図１１５のメイン処理の説明に戻り、リガンド探索装置１００は、相互作用関数計算部１０２ｃの処理により、ステップＳ０で定められた１つのリガンドに対し、リガンドの各原子をステップＳ３０で定められた空間点群に重ね合わせる（ステップＳ５０）。具体的には、距離行列を用いたアライメント作成アルゴリズム（ＤＡＬＩ）［Ｈｏｌｍ， Ｌ．， ＆ Ｓａｎｄｅｒ， Ｃ． （１９９３） Ｐｒｏｔｅｉｎ Ｓｔｒｕｃｔｕｒｅ Ｃｏｍｐａｒｉｓｏｎ ｂｙ Ａｌｉｇｎｍｅｎｔ ｏｆ Ｄｉｓｔａｎｃｅ Ｍａｔｒｉｃｅｓ Ｊ． Ｍｏｌ． Ｂｉｏｌ． ２３３， １２３−１３８］を低分子用に改良した以下の（１）〜（４）に示す手法により初期座標とリガンドとを重ね合わせる。

（１）１つの「Ｂタイプ」には、リガンドの原子タイプが複数対応することが多い。そこで、乱数を用いて「Ｂタイプ」とリガンドの原子タイプで同一視できるペアを作成する。ただし、ペアにおいてリガンドの原子タイプが重複しないようにする。
（２）「Ｂタイプ」には、ステップＳ４０により複数の初期座標が含まれているので、初期座標も乱数を用いて選択する。
（３）選択された初期座標とリガンドそれぞれの距離行列を作成し、以下の相互作用関数であるＳｓｃｏｒｅ（ｉ，ｊ）を計算する。

ここで、ｄ_ij ^sは該当タンパク質中のｉ番目とｊ番目の空間点距離である。また、ｄ_ij ^Cは化合物中のｉ番目とｊ番目の原子間距離である。αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数である。βは重なりと定義できる限界値を与えるための定数である。
また、αは１．５、βは０．８とするのが好ましい。
（４）（１）〜（３）を複数回（例えば１００００回）繰り返し、Ｓｓｃｏｒｅ（ｉ，ｊ）が最大になるペアを探し、そのペア情報をもとにリガンドを初期座標に重ね合わせる。

つぎに、リガンド探索装置１００は、相互作用関数計算部１０２ｃの処理により、ステップＳ５０での重ね合わせに対して、ステップＳ２０およびステップＳ３０で定められた計算結果によりタンパク質の動的挙動を表すパラメータを取得し、該当パラメータを用いてリガンドとタンパク質との相互作用エネルギーを、リガンドのコンフォメーションを微調整しながら計算する（ステップＳ６０）。つまり、ステップＳ５０で重ね合わせたリガンドに対して、タンパク質との相互作用エネルギーを、コンフォメーションを微調整しながら最適化して計算する。リガンドのコンフォメーションの微調整は、ステップＳ５０で算出されたリガンド座標を中心に、例えば、並進、回転、または、シングルボンドまわりの角度をＲＳＭＤで０．３を越えない程度の座標変化などをさせることで、行ってもよい。

ここで、リガンド座標データのコンフォメーションの微調整は、例えばランダムサーチで最適化することが好ましい。ランダムサーチでは、以下の１）〜３）の項目に従って、タンパク質の活性部位とリガンドとの微小変化を例えば８０００回行い、最適エネルギー「Ｕ最適」が最小になるようにする。
１）回転可能な結合のうち最大５つ乱数で選び、結合ごとにランダムに±１０．０°の範囲内で回転させリガンドのコンフォメーションを換える。この過程を例えば３回に一度行う。
２）ｘ、ｙ、ｚ軸方向それぞれにおいて、ランダムに±１．０Åの範囲内でリガンドの並進運動を行う。この過程を例えば２回に一度行う。
３）回転中心座標それぞれにおいて、ランダムに±１．０Åの範囲内で回転中心座標を移動させ、さらに３次元方向の角度それぞれに対して、ランダムに±５．０°の範囲内でリガンドの回転運動を行う。この過程を例えば５回に一度行う。
ここで、最適エネルギー「Ｕ最適」は以下の式で定義する。なお、当該式の右辺に示した各エネルギー関数については以下で順に説明する。

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

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

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

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

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

（ｂ）水素結合に関するエネルギー関数「Ｕ水素」
リガンドの１つのドナー（アクセプター）に対して１つだけ水素結合を形成すると考え、最短にある活性部位側のアクセプター（ドナー）を選び、水素を介した結合角θ（図９参照。ただし、複数の水素原子が付加されているドナー原子の場合には、最小の水素結合角をθと定義する。）を算出して、次の条件により分岐して、ψ（ｉ）を計算するように定義する。なお、図９において、Ａはドナー、Ｂは水素、Ｃはアクセプター、θは水素結合角を示す。

（１）ドナー原子がｓｐ³軌道原子、または水素結合角θが±３０．０°以内のとき

（２）水素結合角θが±３０．０°以上のとき

Ｎはリガンドのドナー＋アクセプターの数であり、Ｒは水素結合を形成する二原子間距離であり、「Ｋ水素（ｉ）」及び「Ｒ水素（ｉ）」は原子タイプごとに決めた水素結合の相互作用の強さ及び距離である。

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

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

（ｄ）スタッキングエネルギー関数「Ｕスタッキング」
活性部位及びリガンドの芳香環を形成する原子に通し番号を付け、活性部位においては芳香環の中心座標を算出した。活性部位のｉ番目とリガンドのｊ番目との原子間距離Ｒがカットオフ以内にあるとき、ｉ番目の原子が形成する芳香環の中心座標をｉ’、ｊ番目の最短距離あり共に同じ芳香環を形成するリガンドの原子をｊ’としたとき、∠ｉｉ’ｊ＝θ_i'j、∠ｉ’ｉｊ＝θ_ij、∠ｉｉ’ｊ’＝θ_i'j'、∠ｉ’ｉｊ’＝θ_ij'を算出し（図１０）、θ_i'jとθ_ijとが９０．０°±１０．０°のとき、「Ｒ境界」と「θ境界」を求め次の条件により分岐してψ（ｉ，ｊ）を計算するように定義する。図１０において、ｉ’は活性部位の芳香環中心を示し、ｉは活性部位の芳香環原子を示し、ｊおよびｊ’はリガンドの芳香環原子を示す。

Ｉｆ Ｒ境界＜０．０，
ψ（ｉ，ｊ）＝−Ｋスタッキング（ｉ，ｊ） ＊ Ｒ境界
Ｅｌｓｅ，
ψ（ｉ，ｊ）＝−Ｋスタッキング（ｉ，ｊ） ＊ θ境界

Ｒ境界＝１．０−（Ｒスタッキング（ｉ，ｊ）−Ｒ）²
θ境界＝｜１．０−Θ｜

Ｍは活性部位の芳香環を形成する原子の数であり、Ｎはリガンドの芳香環を形成する原子の数であり、「Ｋスタッキング（ｉ，ｊ）」及び「Ｒスタッキング（ｉ，ｊ）」は原子タイプごとに決めたスタッキングの強さ及び距離である。πは円周率であり、θはθ_i'j'とθ_ij'においてΘが最小になる角度である。また、カットオフは例えば５．０Åとする。

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

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

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

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

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

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

（ｆ）リガンド内部エネルギー「Ｕ内部」
回転可能な結合を微小に変動させていくと、誤差で結合が切れる恐れがあるために「ψ結合長（ｉ）」を、また、リガンド内部で原子衝突が起こることも避けるために「ψ衝突（ｉ，ｊ）」を計算するように定義する。

ψ結合長（ｉ）＝Ｋ結合長＊｛１０００．０ ＊ （Ｒ結合長（ｉ）−Ｒ₁）｝²
ψ衝突（ｉ，ｊ）＝Ｋ衝突＊（Ｒ衝突−Ｒ₂）²

Ｌは回転可能な結合の数であり、Ｍはリガンドの原子数であり、Ｎはｉ番目の原子の非結合原子数である。また、「Ｋ結合長」は１００．０であることが好ましい。「Ｒ結合長（ｉ）」は初期構造の結合長である。「Ｋ衝突」は１５０．０であり、「Ｒ衝突」は２．２Åであることが好ましい。また、Ｒ₁とＲ₂は二原子間距離とする。

再び図１１５のメイン処理の説明に戻り、リガンド探索装置１００は、相互作用関数計算部１０２ｃの処理により、ステップＳ５０で定められたリガンドに対して、リガンド座標データのコンフォメーションを大きく変動して、ステップＳ５０から再スタートを行い、ステップＳ６０までを繰り返して最適化を繰り返し行う（ステップＳ７０）。なお、コンフォメーション改変は、ステップＳ５０で算出されたリガンド座標を中心に、例えば、並進、回転、または、シングルボンドまわりの角度をＲＳＭＤで０．３以上になるよう座標変化などさせることで、行ってもよい。

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

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

つぎに、リガンド探索装置１００は、相互作用関数計算部１０２ｃの処理により、ステップＳ７０まで得られた標的タンパク質とリガンドとの相互作用エネルギーを決定する（ステップＳ８０）。具体的には、ステップＳ４０からステップＳ７０まで「Ｕ最適」が最適値となる最適なタンパク質とリガンドとの複合体座標、および最適エネルギー「Ｕ最適」を算出する。

つぎに、リガンド探索装置１００は、制御部１０２の処理により、ステップＳ４０に戻り、ステップＳ０中の別のリガンドを選択し、各処理部の処理により、ステップＳ８０まで計算する（ステップＳ９０）。なお、ステップＳ４０からステップＳ９０まではステップＳ０中の化合物データベース中のリガンド全てについて行う。

つぎに、リガンド探索装置１００は、リガンド評価部１０２ｄの処理により、ステップＳ０中でのリガンドに対し、ステップＳ９０において定められた相互作用エネルギーを比較し、標的タンパク質に結合すると予想されるリガンドを選択する（ステップＳ１００）。具体的には、ステップＳ９０まで評価されたタンパク質とリガンドとの複合体座標および最適エネルギー「Ｕ最適」に基づいて、ステップＳ０中のデータベース中から該当タンパク質と結合する可能性のある化合物（リガンド）を選択する。

以上、本システムのメイン処理の説明を終了する。

以上、説明したように、リガンド探索装置１００によれば、相互作用関数によって、該当タンパク質と結合する物質（具体的にはリガンド）を評価し、選定することができる。具体的には、該当タンパク質の立体構造に際し、基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、そのゆらぎの大きさを拘束条件として分子動力学計算を行うことで、タンパク質の立体構造をエネルギー最適構造よりおおきく離れないようにして分子動力学計算を行うことができる。また、タンパク質の動的挙動に関する物理化学的性質を明確に描写することができる。また、タンパク質の動的性質を表現する関数として基準振動解析結果またはタンパク質の二次構造判定結果を用いることができる。また、核磁気共鳴スペクトルの解析結果のような複数の立体構造座標である場合についても、該当タンパク質と結合するリガンドの探索について、複数座標すべてを全自動的にかつ短時間で同等に評価することができる。

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

また、リガンド探索装置１００によれば、該当タンパク質に結合するリガンドを選択する際に（０）〜（８）に示した一連の処理を全自動または手動的に行う。
（０）化合物データベースからリガンドを１つ選択する。該当タンパク立体構造として、誘導適合を反映するパラメータを用いて動的挙動を考慮した複数の構造変化座標を用意し、ランダムに１つの構造を選択する。
（１）重ね合わせを行う該当タンパク中の空間点を指定する。
（２）（０）で選択したリガンド中の原子と（１）で指定した空間点とのペアを重複がないようにランダムに選択する。
（３）以下のスコアＳｓｃｏｒｅ（ｉ，ｊ）を計算する。

ｄ_ij ^Sは該当タンパク質中のｉ番目とｊ番目の空間点距離である。ｄ_ij ^Cは化合物中のｉ番目とｊ番目の原子間距離である。αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数である。βは重なりと定義できる限界値を与えるための定数である。
（４）（３）のスコアが最大になるように調整する。
（５）（４）で重ねあわせたリガンドに対してタンパク質との相互作用エネルギーをコンフォメーションを微調整しながら最適化計算する。
（６）リガンドのコンフォメーションを大きく動かして、（２）から再スタートを行い、（５）までを繰り返して最適化を行う。
（７）（１）〜（６）までの過程を（０）で用意した複数の構造変化座標に対して行い、最適なタンパク質とリガンドとの複合体座標、最適エネルギー「Ｕ最適」を算出する。
（８）（１）〜（７）までの過程を（０）で用意した化合物データベース中の全てのリガンドに対して行い、化合物データベース中から該当タンパク質と結合する可能性のあるリガンドを選択する。

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

また、リガンド探索装置１００によれば、リガンドとタンパク質との相互作用を評価する際の目的関数として、従来の相互作用エネルギー関数に、タンパク質の動的性質を表現する関数を弾性エネルギーとして加え、タンパク質の立体構造座標から相互作用エネルギーを高速に算出するとともに、タンパク質の動的挙動に関する物理化学的性質を明確に描写する。

また、リガンド探索装置１００によれば、タンパク質の動的性質を表現する関数として基準振動解析結果またはタンパク質の二次構造判定結果を用いる。

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

（二面角拘束ＭＤおよびクラスタリングにおけるパラメータ定数の決定）
上述した実施形態におけるリガンド探索装置１００を用いて、基準振動解析により二面角のゆらぎ値を計算した。なお、本実施例１においては、二面角のゆらぎ値を、分子動力学計算における拘束条件として、以下の式中の「Ｋ二面角」に適応した。

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

θは二面角（単位ｒａｄ）である。実際には、「Ｋ二面角」の最大値と最小値を指定することにより、その範囲内で主鎖二面角揺らぎに対応するように各二面角に対して不均一な拘束がかかるようにした。そのため、本実施例１では、「Ｋ二面角」の適切な最大値と最小値を決定することを目的とした。

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

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

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

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

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

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

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

図１１に１ＬＵＤのＭＯＤＥＬ１における基準振動解析結果を示し、図１２、図１３、図１５〜１８にスコアとパラメータとの比較の結果を示した。図１１では、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示した。紫色Ｃは活性部位である。図１２において、基準振動解析の結果、定数αは７０％が良いという結果になった。しかし、一般性を持たせときに７０％ではクラスタリングの精度が低下した場合もあったので、本実施例１では８０％を定数αの値にした。図１３、図１５〜１８において、クラスタリング定数をα＝８０．０％、β＝０．４Åに固定した。なお、黒色に近いほどスコアが小さいことを表わす。

これらの結果より、スコアが小さくなる拘束条件としては図１４の値が最適であると判断された。これらの値の妥当性は、例えば、主鎖原子のみではなく、Ｃα原子、側鎖原子、全原子について調査しても図１４のパラメータ値が最適であることが分かった。

（拘束パラメータ有無による分子動力学計算の相違）
上述したリガンド探索装置１００により算出された拘束パラメータを適応した分子動力学計算を２．０ｎｓｅｃまで行った。そして、活性部位の主鎖原子の動的挙動において、拘束パラメータを適応しない場合と比較して、構造がどの程度変化するかを調べた。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

図５３、図５５に１ＹＥＲを用いた拘束二面角分子動力学の結果を示した。正解構造１ＹＥＴの活性部位とのｒｍｓを表示した。Ａは主鎖原子、Ｂは側鎖原子、Ｃは全原子、である。図５５は、１ＹＥＲにおける活性部位・リガンド結合解析の結果である。ＭＤ計算を０．１ｎｓｅｃまで行ったときの結合解析及び１．０ｎｓｅｃまで行ったとき、また、受容体動的構造クラスタリングを行うときの母集団を１００ｆｓｅｃごと及び１０００ｆｓｅｃごとで行ったときの結合解析である。評価法は、正解構造との活性部位及びリガンドのｒｍｓで行った。図５４にリガンドドッキングにおける空間点指定のパラメータ値を示した。図５４は、１ＹＥＴより得られた構造活性相関情報である。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Ｃａｓｅ１）ＳＡＲ４ヵ所指定
図７９に１ＵＫ４の活性部以内での空間指定を示した。図７９は１ＵＫ４（Ｂ鎖）より得られた構造活性相関情報である。図８０に、１ＵＫ３（Ｂ鎖）におけるｉｎ ｓｉｌｉｃｏ スクリーニングの結果を示した。

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

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

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

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

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

Ｃａｓｅ３）ＳＡＲ５ヵ所指定
図８７には１ＵＫ３の活性部位内での空間指定を示した。図８７は１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報に関する図である。図８８には、ＳＡＲ５ヵ所指定で実行したｉｎ ｓｉｌｉｃｏスクリーニング（ハイスループットスクリーニング）結果を示した。

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

Ｃａｓｅ４）リガンド原子タイプ指定の変更
図９０に１ＵＫ３の活性部位内での空間指定を示した。図９０は１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報に関する図である。図９１には、リガンド原子タイプ指定変更で実行したｉｎ ｓｉｌｉｃｏスクリーニング（ハイスループットスクリーニング）の結果を示した。

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

Ｃａｓｅ５）受容体固定
図９３に活性部位内での空間指定を示した。図９３は１ＵＫ４（Ｂ鎖）より得られた構造活性相関情報に関する図である。図９４に、受容体を固定した状態で実行したｉｎ ｓｉｌｉｃｏ スクリーニング（ハイスループットスクリーニング）の結果を示した。

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

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

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

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

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

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

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

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

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

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

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

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

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

図１１２には、１ＢＺＦの活性部位内における空間指定を示した。図１１２は１ＢＺＦ用に変更した構造活性相関情報に関する図である。図１１３および図１１４には、１ＢＺＦにおける受容体／リガンド結合の結果を示した。図１１３および図１１４は１ＢＺＦ（ＭＯＤＥＬ１８）のリガンド結合解析の結果である。
（ｉ）最適化したときの受容体には初期構造を選択した。リガンドのｒｍｓは０．８８８４であった（図１１３）。
（ｉｉ）ｔｒｉｍｅｔｒｅｘａｔｅ、１ＢＺＦ（ＭＯＤＥＬ１８）のリガンドである（図１１４）。

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

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

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

符号の説明

１００ リガンド探索装置
１０２ 制御部
１０２ａ 構造変化後タンパク質座標データ選択部
１０２ｂ 空間点指定部
１０２ｃ 相互作用関数計算部
１０２ｄ リガンド評価部
１０４ 通信制御インターフェース部
１０６ 記憶部
１０６ａ リガンド座標データベース
１０６ｂ タンパク質座標データベース
１０６ｃ 構造変化後タンパク質座標データファイル
１０６ｄ 指定空間点ファイル
１０６ｅ 相互作用関数計算結果ファイル
１０６ｆ リガンド評価結果ファイル
１０８ 入出力制御インターフェース部
１１２ 入力装置
１１４ 出力装置
２００ 外部システム
３００ ネットワーク

Claims (31)

1. 単数または複数鎖のタンパク質の座標データが与えられた場合に、当該タンパク質と結合するリガンドを探索するリガンド探索装置において、
   上記タンパク質の座標データに対して、誘導適合を反映した誘導適合パラメータを用いて動的挙動を考慮した構造変化を行い、構造変化後タンパク質座標データを選択する構造変化後タンパク質座標データ選択手段と、
   上記構造変化後タンパク質座標データ選択手段にて選択された上記構造変化後タンパク質座標データから、上記リガンドと重ね合わせを行う空間点を指定する空間点指定手段と、
   上記空間点指定手段にて指定された上記空間点と、上記リガンドのリガンド座標データとを用いて、上記タンパク質と上記リガンドとが結合した場合の相互作用関数を計算する相互作用関数計算手段と、
   上記相互作用関数計算手段により計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドを評価するリガンド評価手段と、
   を備えたことを特徴とするリガンド探索装置。
2. 請求項１に記載のリガンド探索装置において、上記相互作用関数計算手段は、以下の数式１に示すＳｓｃｏｒｅ（ｉ，ｊ）を用いて上記相互作用関数を計算することを特徴とするリガンド探索装置。
   

   

   ・・・（数式１）
   （ここで、ｄ_ij ^sは該当タンパク質中のｉ番目とｊ番目の空間点距離である。また、ｄ_ij ^Cは化合物中のｉ番目とｊ番目の原子間距離である。また、αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数である。また、βは重なりと定義できる限界値を与えるための定数である。）
3. 請求項１または２に記載のリガンド探索装置において、上記相互作用関数計算手段は、
   上記相互作用関数のスコアが最大になるように最適化する相互作用関数最適化手段、
   をさらに備えたことを特徴とするリガンド探索装置。
4. 請求項３に記載のリガンド探索装置において、上記相互作用関数計算手段は、
   上記相互作用関数最適化手段により上記相互作用関数を最適化した後に、重ねあわせた上記リガンドに対して、上記タンパク質との相互作用エネルギーを計算し、当該相互作用エネルギーについてリガンド立体構造データのコンフォメーションを微調整しながら最適化する相互作用エネルギー最適化手段、
   をさらに備えたことを特徴とするリガンド探索装置。
5. 請求項４に記載のリガンド探索装置において、上記リガンド評価手段は、
   上記相互作用エネルギー最適化手段により最適化した後に、上記リガンド立体構造データのコンフォメーションを大きく変動させた後、再度、上記相互作用関数計算手段を実行し、上記相互作用関数計算手段により計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドの再評価を行う再評価手段、
   をさらに備えたことを特徴とするリガンド探索装置。
6. 請求項１から５のいずれか１つに記載のリガンド探索装置において、上記構造変化後タンパク質座標データ選択手段は、
   上記誘導適合パラメータおよび／または上記構造変化後タンパク質座標データを算出するときに、上記タンパク質座標データに対し基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、当該ゆらぎの大きさを拘束条件として、分子動力学計算を行うこと、
   を特徴とするリガンド探索装置。
7. 請求項６に記載のリガンド探索装置において、上記構造変化後タンパク質座標データ選択手段は、
   基準振動計算より主鎖原子の２面角のゆらぎの値を算出し、当該ゆらぎ値を以下の数式２または数式３に示す分子動力学計算における力の定数Ｋとすることにより、分子動力学計算を行うこと、
   を特徴とするリガンド探索装置。
   
   Ｅｒｏｔ＝Ｋｒｏｔ（φ―φ０）² ・・・（数式２）
   
   （Ｅｒｏｔはタンパク質の立体構造中において主鎖原子の２面角のエネルギーを示す。φは主鎖原子の２面角である。φ０は主鎖原子の２面角の標準値である。ここで、Ｋｒｏｔの値が大きい場合は、φはφ０に拘束される。）
   
   Ｅｐｏｓ＝Ｋｐｏｓ（ｒ−ｒ０）² ・・・（数式３）
   
   （Ｅｐｏｓはタンパク質の立体構造中において主鎖原子の位置のエネルギーを示す。ｒは主鎖原子の座標である。ｒ０は主鎖原子の座標の標準値である。ここで、Ｋｐｏｓの値が大きい場合は、ｒはｒ０に拘束される。）
8. 請求項１から７のいずれか１つに記載のリガンド探索装置において、上記相互作用関数計算手段は、
   上記相互作用関数として、相互作用エネルギー関数にタンパク質の動的性質を表現する動的性質関数を「弾性エネルギー」として加えて用いること、
   を特徴とするリガンド探索装置。
9. 請求項８に記載のリガンド探索装置において、上記相互作用関数計算手段は、
   「弾性エネルギー」として、タンパク質の局所的な柔らかさを考慮し、以下の数式４に示す関数「Ｕ衝突」を適応すること、
   を特徴とするリガンド探索装置。
   

   

   ψ（ｉ，ｊ）＝Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²
   ・・・（数式４）
   （Ｍは衝突を不許可とする活性部位の原子の数であり、Ｎはリガンドの原子の数である。活性部位における動的挙動の少ない側鎖原子及び主鎖原子のみのｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義する。）
10. 請求項１から７のいずれか１つに記載のリガンド探索装置において、上記相互作用関数計算手段は、
    上記相互作用関数として、相互作用エネルギー関数に、当該タンパク質の基準振動解析結果または二次構造判定結果をタンパク質の動的性質を表現する動的性質関数として加えて用いること、
    を特徴とするリガンド探索装置。
11. 単数または複数鎖のタンパク質の座標データが与えられた場合に、当該タンパク質と結合するリガンドを探索するリガンド探索方法において、
    上記タンパク質の座標データに対して、誘導適合を反映した誘導適合パラメータを用いて動的挙動を考慮した構造変化を行い、構造変化後タンパク質座標データを選択する構造変化後タンパク質座標データ選択ステップと、
    上記構造変化後タンパク質座標データ選択ステップにて選択された上記構造変化後タンパク質座標データから、上記リガンドと重ね合わせを行う空間点を指定する空間点指定ステップと、
    上記空間点指定ステップにて指定された上記空間点と、上記リガンドのリガンド座標データとを用いて、上記タンパク質と上記リガンドとが結合した場合の相互作用関数を計算する相互作用関数計算ステップと、
    上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドを評価するリガンド評価ステップと、
    を含むことを特徴とするリガンド探索方法。
12. 請求項１１に記載のリガンド探索方法において、上記相互作用関数計算ステップは、以下の数式１に示すＳｓｃｏｒｅ（ｉ，ｊ）を用いて上記相互作用関数を計算することを特徴とするリガンド探索方法。
    

    

    ・・・（数式１）
    （ここで、ｄ_ij ^sは該当タンパク質中のｉ番目とｊ番目の空間点距離である。また、ｄ_ij ^Cは化合物中のｉ番目とｊ番目の原子間距離である。また、αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数である。また、βは重なりと定義できる限界値を与えるための定数である。）
13. 請求項１１または１２に記載のリガンド探索方法において、上記相互作用関数計算ステップは、
    上記相互作用関数のスコアが最大になるように最適化する相互作用関数最適化ステップ、
    をさらに含むことを特徴とするリガンド探索方法。
14. 請求項１３に記載のリガンド探索方法において、上記相互作用関数計算ステップは、
    上記相互作用関数最適化ステップにより上記相互作用関数を最適化した後に、重ねあわせた上記リガンドに対して、上記タンパク質との相互作用エネルギーを計算し、当該相互作用エネルギーについてリガンド立体構造データのコンフォメーションを微調整しながら最適化する相互作用エネルギー最適化ステップ、
    をさらに含むことを特徴とするリガンド探索方法。
15. 請求項１４に記載のリガンド探索方法において、上記リガンド評価ステップは、
    上記相互作用エネルギー最適化ステップにより最適化した後に、上記リガンド立体構造データのコンフォメーションを大きく変動させた後、再度、上記相互作用関数計算ステップを実行し、上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドの再評価を行う再評価ステップ、
    をさらに含むことを特徴とするリガンド探索方法。
16. 請求項１１から１５のいずれか１つに記載のリガンド探索方法において、上記構造変化後タンパク質座標データ選択ステップは、
    上記誘導適合パラメータおよび／または上記構造変化後タンパク質座標データを算出するときに、上記タンパク質座標データに対し基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、当該ゆらぎの大きさを拘束条件として、分子動力学計算を行うこと、
    を特徴とするリガンド探索方法。
17. 請求項１６に記載のリガンド探索方法において、上記構造変化後タンパク質座標データ選択ステップは、
    基準振動計算より主鎖原子の２面角のゆらぎの値を算出し、当該ゆらぎ値を以下の数式２または数式３に示す分子動力学計算における力の定数Ｋとすることにより、分子動力学計算を行うこと、
    を特徴とするリガンド探索方法。
    
    Ｅｒｏｔ＝Ｋｒｏｔ（φ―φ０）² ・・・（数式２）
    
    （Ｅｒｏｔはタンパク質の立体構造中において主鎖原子の２面角のエネルギーを示す。φは主鎖原子の２面角である。φ０は主鎖原子の２面角の標準値である。ここで、Ｋｒｏｔの値が大きい場合は、φはφ０に拘束される。）
    
    Ｅｐｏｓ＝Ｋｐｏｓ（ｒ−ｒ０）² ・・・（数式３）
    
    （Ｅｐｏｓはタンパク質の立体構造中において主鎖原子の位置のエネルギーを示す。ｒは主鎖原子の座標である。ｒ０は主鎖原子の座標の標準値である。ここで、Ｋｐｏｓの値が大きい場合は、ｒはｒ０に拘束される。）
18. 請求項１１から１７のいずれか１つに記載のリガンド探索方法において、上記相互作用関数計算ステップは、
    上記相互作用関数として、相互作用エネルギー関数にタンパク質の動的性質を表現する動的性質関数を「弾性エネルギー」として加えて用いること、
    を特徴とするリガンド探索方法。
19. 請求項１８に記載のリガンド探索方法において、上記相互作用関数計算ステップは、
    「弾性エネルギー」として、タンパク質の局所的な柔らかさを考慮し、以下の数式４に示す関数「Ｕ衝突」を適応すること、
    を特徴とするリガンド探索方法。
    

    

    ψ（ｉ，ｊ）＝Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²
    ・・・（数式４）
    （Ｍは衝突を不許可とする活性部位の原子の数であり、Ｎはリガンドの原子の数である。活性部位における動的挙動の少ない側鎖原子及び主鎖原子のみのｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義する。）
20. 請求項１１から１７のいずれか１つに記載のリガンド探索方法において、上記相互作用関数計算ステップは、
    上記相互作用関数として、相互作用エネルギー関数に、当該タンパク質の基準振動解析結果または二次構造判定結果をタンパク質の動的性質を表現する動的性質関数として加えて用いること、
    を特徴とするリガンド探索方法。
21. 単数または複数鎖のタンパク質の座標データが与えられた場合に、当該タンパク質と結合するリガンドを探索するリガンド探索方法をコンピュータに実行させるプログラムにおいて、
    上記タンパク質の座標データに対して、誘導適合を反映した誘導適合パラメータを用いて動的挙動を考慮した構造変化を行い、構造変化後タンパク質座標データを選択する構造変化後タンパク質座標データ選択ステップと、
    上記構造変化後タンパク質座標データ選択ステップにて選択された上記構造変化後タンパク質座標データから、上記リガンドと重ね合わせを行う空間点を指定する空間点指定ステップと、
    上記空間点指定ステップにて指定された上記空間点と、上記リガンドのリガンド座標データとを用いて、上記タンパク質と上記リガンドとが結合した場合の相互作用関数を計算する相互作用関数計算ステップと、
    上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドを評価するリガンド評価ステップと、
    を含むリガンド探索方法をコンピュータに実行させることを特徴とするプログラム。
22. 請求項２１に記載のプログラムにおいて、上記相互作用関数計算ステップは、以下の数式１に示すＳｓｃｏｒｅ（ｉ，ｊ）を用いて上記相互作用関数を計算することを特徴とするプログラム。
    

    

    ・・・（数式１）
    （ここで、ｄ_ij ^sは該当タンパク質中のｉ番目とｊ番目の空間点距離である。また、ｄ_ij ^Cは化合物中のｉ番目とｊ番目の原子間距離である。また、αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数である。また、βは重なりと定義できる限界値を与えるための定数である。）
23. 請求項２１または２２に記載のプログラムにおいて、上記相互作用関数計算ステップは、
    上記相互作用関数のスコアが最大になるように最適化する相互作用関数最適化ステップ、
    をさらに含むことを特徴とするプログラム。
24. 請求項２３に記載のプログラムにおいて、上記相互作用関数計算ステップは、
    上記相互作用関数最適化ステップにより上記相互作用関数を最適化した後に、重ねあわせた上記リガンドに対して、上記タンパク質との相互作用エネルギーを計算し、当該相互作用エネルギーについてリガンド立体構造データのコンフォメーションを微調整しながら最適化する相互作用エネルギー最適化ステップ、
    をさらに含むことを特徴とするプログラム。
25. 請求項２４に記載のプログラムにおいて、上記リガンド評価ステップは、
    上記相互作用エネルギー最適化ステップにより最適化した後に、上記リガンド立体構造データのコンフォメーションを大きく変動させた後、再度、上記相互作用関数計算ステップを実行し、上記相互作用関数計算ステップにより計算された上記相互作用関数に基づいて当該タンパク質と結合する上記リガンドの再評価を行う再評価ステップ、
    をさらに含むことを特徴とするプログラム。
26. 請求項２１から２５のいずれか１つに記載のプログラムにおいて、上記構造変化後タンパク質座標データ選択ステップは、
    上記誘導適合パラメータおよび／または上記構造変化後タンパク質座標データを算出するときに、上記タンパク質座標データに対し基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、当該ゆらぎの大きさを拘束条件として、分子動力学計算を行うこと、
    を特徴とするプログラム。
27. 請求項２６に記載のプログラムにおいて、上記構造変化後タンパク質座標データ選択ステップは、
    基準振動計算より主鎖原子の２面角のゆらぎの値を算出し、当該ゆらぎ値を以下の数式２または数式３に示す分子動力学計算における力の定数Ｋとすることにより、分子動力学計算を行うこと、
    を特徴とするプログラム。
    
    Ｅｒｏｔ＝Ｋｒｏｔ（φ―φ０）² ・・・（数式２）
    
    （Ｅｒｏｔはタンパク質の立体構造中において主鎖原子の２面角のエネルギーを示す。φは主鎖原子の２面角である。φ０は主鎖原子の２面角の標準値である。ここで、Ｋｒｏｔの値が大きい場合は、φはφ０に拘束される。）
    
    Ｅｐｏｓ＝Ｋｐｏｓ（ｒ−ｒ０）² ・・・（数式３）
    
    （Ｅｐｏｓはタンパク質の立体構造中において主鎖原子の位置のエネルギーを示す。ｒは主鎖原子の座標である。ｒ０は主鎖原子の座標の標準値である。ここで、Ｋｐｏｓの値が大きい場合は、ｒはｒ０に拘束される。）
28. 請求項２１から２７のいずれか１つに記載のプログラムにおいて、上記相互作用関数計算ステップは、
    上記相互作用関数として、相互作用エネルギー関数にタンパク質の動的性質を表現する動的性質関数を「弾性エネルギー」として加えて用いること、
    を特徴とするプログラム。
29. 請求項２８に記載のプログラムにおいて、上記相互作用関数計算ステップは、
    「弾性エネルギー」として、タンパク質の局所的な柔らかさを考慮し、以下の数式４に示す関数「Ｕ衝突」を適応すること、
    を特徴とするプログラム。
    

    

    ψ（ｉ，ｊ）＝Ｋ衝突 ＊ （Ｒ衝突（ｉ，ｊ）−Ｒ）²
    ・・・（数式４）
    （Ｍは衝突を不許可とする活性部位の原子の数であり、Ｎはリガンドの原子の数である。活性部位における動的挙動の少ない側鎖原子及び主鎖原子のみのｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義する。）
30. 請求項２１から２７のいずれか１つに記載のプログラムにおいて、上記相互作用関数計算ステップは、
    上記相互作用関数として、相互作用エネルギー関数に、当該タンパク質の基準振動解析結果または二次構造判定結果をタンパク質の動的性質を表現する動的性質関数として加えて用いること、
    を特徴とするプログラム。
31. 請求項２１から３０のいずれか１つに記載のプログラムを記録したことを特徴とするコンピュータ読み取り可能な記録媒体。

Images