(57)【要約】【課題】本発明は、誘導適合型受容体も含めた受容体／リガンド結合の解析を研究する方法の提供を目的としてなされたものである。【解決手段】本発明は、受容体の基準振動解析計算をまず行い、定常状態の受容体の主鎖二面角の揺らぎを計算する。そして、その揺らぎを基にした拘束を各原子にかけながら分子動力学計算を行うことでより精度の良い受容体の動的構造を予測する。また、本発明は、分子動力学計算より得た動的構造及び相互作用関数を用いることで、誘導適合型受容体にも応用できる受容体／リガンド結合を精度良く予測する。【選択図】        図１
2005242493.tif 000002
【特許請求の範囲】
【請求項１】
任意の単数を含む複数鎖のタンパク質立体構造が与えられた場合において、該当タンパク質の立体構造から誘導適合を反映したパラメータおよび構造変化した立体構造座標を例えば基準振動計算方法や分子動力学計算方法よりあらかじめ算出し、当該パラメータおよび構造変化した立体構造座標を用いて該当タンパク質と別の物質が結合した場合の相互作用関数を定義し、当該相互作用関数によって該当タンパク質と結合する物質をコンピュータプログラムにより評価し、選定する方法及びその装置。
【請求項２】
請求項１によって記載された方法において、該当タンパク質に結合するリガンドを選択する際に（０）〜（８）に示した一連の処理を全自動または手動的に行うことを特徴とするコンピュータプログラム及びその装置。
（０）  化合物データベースからリガンドを１つ選択する。該当タンパク立体構造として、誘導適合を反映するパラメータを用いて動的挙動を考慮した複数の構造変化座標を用意し、ランダムに１つの構造を選択する。
（１）  重ね合わせを行う該当タンパク中の空間点を指定する。
（２）  （０）で選択したリガンド中の原子と（１）で指定した空間点とのペアを重複がないようにランダムに選択する。
（３）  以下のスコアＳｓｃｏｒｅ（ｉ，ｊ）を計算する。
【数１】
2005242493.tif 000027
ｄijSは該当タンパク質中のｉ番目とｊ番目の空間点距離。ｄijCは化合物中のｉ番目とｊ番目の原子間距離。αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数。βは重なりと定義できる限界値を与えるための定数。
（４）  （３）のスコアが最大になるように調整する。
（５）  （４）で重ねあわせたリガンドに対してタンパク質との相互作用エネルギーをコンフォメーションを微調整しながら最適化計算する。
（６）  リガンドのコンフォメーションを大きく動かして、（２）から再スタートを行い、（５）までを繰り返して最適化を行う。
（７）  （１）〜（６）までの過程を（０）で用意した複数の構造変化座標に対して行い、最適なタンパク質とリガンドとの複合体座標、最適エネルギー「Ｕ最適」を算出する。
（８）  （１）〜（７）までの過程を（０）で用意した化合物データベース中の全てのリガンドに対して行い、化合物データベース中から該当タンパク質と結合する可能性のあるリガンドを選択する。
【請求項３】
請求項１〜請求項２によって記載された方法において、タンパク質の誘導適合を反映するパラメータおよび構造変化した立体構造座標を分子動力学計算方法を用いて算出する場合、該当タンパク質の立体構造に際し、基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、そのゆらぎの大きさを拘束条件として、分子動力学計算を行うことで、タンパク質の立体構造をエネルギー最適構造よりおおきく離れないようにして分子動力学計算を行う方法及びその装置。
【請求項４】
請求項１〜請求項２によって記載された方法において、リガンドとタンパク質との相互作用を評価する際の目的関数として、従来の相互作用エネルギー関数に、タンパク質の動的性質を表現する関数を弾性エネルギーとして加え、タンパク質の立体構造座標から相互作用エネルギーを高速に算出するとともに、タンパク質の動的挙動に関する物理化学的性質を明確に描写することを特徴とする手法及びその装置。
【請求項５】
請求項１〜請求項４によって記載された方法において、タンパク質の動的性質を表現する関数として基準振動解析結果またはタンパク質の二次構造判定結果を用いることを特徴とした方法及びその装置。
【請求項６】
請求項１〜請求項５によって記載された方法において、計算された該当タンパク質が代表的な複数の立体構造座標である場合や、該当タンパク質の立体構造が例えば核磁気共鳴スペクトルの解析結果のような複数の立体構造座標である場合についても、該当タンパク質と結合するリガンドの探索について、複数座標すべてを全自動的にかつ短時間で同等に評価することを可能とした方法及びその装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明はタンパク質の立体構造座標を用いたリガンド探索方法、詳しくはタンパク質立体構造座標が既知の場合、相互作用すると考えられるリガンドを予測する方法に関する。
更に、本発明は、この方法で得られるリガンド、前期方法に使用可能なデータベース、データベース構造、及びコンピューターソフトプログラム、これを搭載したコンピュータやインターフェース等にも関する。
【背景技術】
【０００２】
酵素や受容体等の生体機能を維持するために必要なタンパク質には、基質特異性と呼ばれる性質があり、活性部位が基質分子構造の細部にわたり常に一致しているＬｏｃｋ＆Ｋｅｙ型と、基質が無いときには活性部位が不活性なランダムな状態にあり、基質が来るとこれを取り込むために活性部位が活性な状態に変化するＩｎｄｕｃｅｄ−Ｆｉｔ（誘導結合）型がある。誘導適合型とは、リガンドと結合する際にリガンド結合部位の立体構造が変化しリガンドを取り込むことが可能になる受容体をいう。
【０００３】
タンパク質の立体構造を用いたリガンド分子探索のための計算化学的手法としてはまずＤＯＣＫ、ＦｌｅｘＸ、Ｌｕｄｉ、ＧＯＬＤ、といった３次元化合物データベースサーチ（Ｖｉｒｔｕａｌ  Ｓｃｒｅｅｎｉｎｇ）が知られている。これらは高速ドッキングスタディーとも呼ばれ大規模な化合物ライブラリサーチが可能である。しかし、本手法では評価に粗い近似を用いるため、結合配座や結合エネルギーの予測能は低い。さらにタンパク質とリガンドとの結合に大変重要な「誘導結合」に対応する計算式パラメータを充分に取り込んでいないので、たとえあったとしても、乱数を発生させ受容体の側鎖を動かす程度であり、計算結果の精度に充分なものとはいえない。
【０００４】
タンパク質とリガンドとの結合に重要な「誘導結合」をシミュレーションする方法としてはＭＤ（分子動力学計算）、ＭＭ（分子力学計算）ＭＣ（モンテカルロ法）が知られている。これらの方法は比較的精度良く、結合配座や結合エネルギーの予測が可能である。ここで、分子動力学法（ＭＤ）と呼ばれる手法に関しては、ある分子を構成する各原子において、古典力学に基づく運動方程式を逐次的に解くことにより、その分子の動的構造を計算する方法であり、タンパク質の動的挙動を高精度でシュミレーションすることが可能である。しかし、計算に時間を要するため、多数の分子を扱うことは困難であり必ずしも有用な手法とはなっていない。さらに、従来法では該当タンパク質に対して分子動力学計算を行うとタンパク質立体構造はＸ線、ＮＭＲ等で解析された座標から大きくズレる。こうしたズレはタンパク質の動的挙動の物理化学的描写を含んでいるがＮＭＲ等で示される動的挙動の実験的な結果と矛盾する挙動となる場合があり、必ずしも精度の高いシミュレーションとならないことが多い。
【０００５】
このように従来のｉｎ  ｓｉｌｉｃｏ  ｓｃｒｅｅｎｉｎｇ関連では、タンパク質とリガンドとの結合に大変重要な「誘導結合」に対応する計算式パラメータを充分に取り込んでいないので、計算結果の精度に充分なものとはいえない。
【０００６】
一方、分子シミュレーションでは上記の誘導結合を表現し、解析することは可能であるが、高精度の結果を得るためには相当の時間を必要とする。多くの結果は、初期構造座標に依存してしまう。
【発明の開示】
【発明が解決しようとする課題】
【０００７】
本発明者等は、任意のタンパク質の立体構造が与えられたとき、該当タンパク質に結合するリガンドを探索する方法について検討をおこなった。現在流通している受容体・リガンド結合解析ソフトには、リガンドのフレキシビリティを考慮しているものは多くあるが、受容体側のフレキシビリティを考慮しているものはほとんどない。たとえあったとしても、乱数を発生させ受容体の側鎖を動かす程度であり、Ｌｏｃｋ＆Ｋｅｙ型の受容体に対応しているものばかりであった。そこで、Ｉｎｄｕｃｅｄ−Ｆｉｔ型の受容体を対象にした受容体・リガンド結合解析ソフト開発することにした。
【０００８】
本発明が解決しようとする課題は、農薬、医薬品等の開発に特に重要な鍵となる、該当タンパク質に結合するリガンドを探索する方法を精度よく、かつ従来法に比べてはるかに効率的に精度よく探索する方法を提供することである。また、リガンド分子の多様な改変や受容体等のタンパク質の改変を迅速かつ効率的に行う方法を提供することにもある。更に本発明により、リガンド−タンパク質間の相互作用様式を解明し、それ等相互作用の認識機構を明確化することで、疾病の原因を特定したり、それ等に関連する薬物の開発を促進すること等を目的とする。
【課題を解決するための手段】
【０００９】
本発明者等は、任意のタンパク質立体構造が与えられたとき該当タンパク質に結合するリガンドを探索する方法について検討を重ねた結果、下記［１］〜［６］の方法およびそのためのコンピュータプログラムを見出し、或いは開発した。
【００１０】
ここで、分子動力学法（ＭＤ）と呼ばれる手法があり、これはある分子を構成する各原子において、古典力学に基づく運動方程式を逐次的に解くことにより、その分子の動的構造を計算する方法である。つまり、これはある分子を構成する各原子における古典力学を土台とした動的挙動を計算する方法である。従って、この手法をうまく取り込むことができれば、リガンドを取り込んでいない状態のＩｎｄｕｃｅｄ−Ｆｉｔ型受容体を初期状態に選んでも受容体・リガンド結合を再現できると考えた。ＭＤ計算は古典力学を土台にしているため、各原子にある程度の拘束をかける必要がある。そこで、まず初めに受容体の基準振動解析を行い受容体の主鎖二面角揺らぎを計算し、この主鎖二面角揺らぎに基づいて各原子に拘束をかけてＭＤを計算する手法を開発した。具体的には、基準振動解析計算をまず行い、定常状態の主鎖二面角の揺らぎを計算する。そして、その揺らぎを基にした拘束を各原子にかけながら分子動力学計算を行うことでより精度の良い受容体の動的構造を予測する。また、分子動力学計算より得た動的構造及び相互作用関数を用いることで、誘導適合型受容体にも応用できる受容体／リガンド結合を精度良く予測できる。本発明の方法は、より真に近い受容体／リガンド結合を予測する方法である。従って、本発明の方法は、医農薬分子の設計に極めて有用である。
【００１１】
［１］  任意の単数を含む複数鎖のタンパク質立体構造が与えられた場合において、該当タンパク質の立体構造から誘導適合を反映したパラメータおよび構造変化した立体構造座標を例えば基準振動計算方法や分子動力学計算方法よりあらかじめ算出し、当該パラメータおよび構造変化した立体構造座標を用いて該当タンパク質と別の物質が結合した場合の相互作用関数を定義し、当該相互作用関数によって該当タンパク質と結合する物質をコンピュータプログラムにより評価し、選定する方法。
【００１２】
［２］  ［１］によって記載された方法において、該当タンパク質に結合するリガンドを選択する際に（０）〜（８）に示した一連の処理を全自動または手動的に行うことを特徴とするコンピュータプログラム。
【００１３】
（０）  化合物データベースからリガンドを１つ選択する。該当タンパク立体構造として、誘導適合を反映するパラメータを用いて動的挙動を考慮した複数の構造変化座標を用意し、ランダムに１つの構造を選択する。
【００１４】
（１）  重ね合わせを行う該当タンパク中の空間点を指定する。空間点は例えば以下のような方法で発生させる。
【００１５】
１．  ダミー原子の発生による空間点の発生
リガンドとタンパク質との相互作用における水素結合に着目し、タンパク質中の水素結合サイトを空間点として指定する。水素結合における重要事項は距離と角度である。つまり、角度を計算するためには、水素結合ドナー（以後、ドナー）に水素原子が必要になる。
【００１６】
そこで、活性部位及びリガンドに水素原子が含まれていない場合、以下の規則によりダミー水素原子を発生させた。
【００１７】
１）  ｓｐ2軌道原子を中心とする正三角形状にダミー原子を発生（図２）。すなわち、図２に示すように、ｓｐ2軌道原子の窒素原子（Ａ）を中心とする正三角形の空いている位置にダミー水素原子（Ｂ）を発生させた。
【００１８】
２）  ｓｐ3軌道原子では、水素結合を形成する距離にある場合、水素原子を共有するように回転できると考え水素結合相互作用を計算するときには、距離のみを考慮することにした。このため、ｓｐ3軌道原子にはダミー原子を発生させない。
【００１９】
金属及び水では、活性部位・リガンド結合の仲介役となり得るので、相互作用する位置に以下のようにダミー原子を発生させた。
【００２０】
１）  鉄のような金属には、正八面体状にダミー原子を発生（図３）。すなわち、図３に示すように、亜鉛（Ａ）を中心とする正八面体の空いている位置にダミー原子（Ｂ）を発生させた。
【００２１】
２）  水は、正四面体状にダミー原子を発生。
【００２２】
ただし、活性部位と相互作用している方向にはダミー原子を発生させないことにした。
【００２３】
２．  構造活性相関情報を利用した空間点の発生
リガンドの構造活性相関（ＳＡＲ）情報に着目し、以下の項目を入力情報にすることにした。
【００２４】
（Ａ）  ＳＡＲから得られた活性部位の原子（以後、「Ａ原子」）。ＰＤＢ形式に従う。
【００２５】
（Ｂ）  「Ａ原子」と相互作用するであろうリガンドの原子タイプ（以後、「Ｂタイプ」）。ＳＹＢＹＬのＭＯＬ２形式に従う。
【００２６】
（Ｃ）  「Ａ原子」と「Ｂタイプ」との相互作用の強さ（以後、「Ｃ強さ」）。
【００２７】
（Ｄ）  「Ａ原子」と「Ｂタイプ」との相互作用する距離（以後、「Ｄ距離」）（単位はÅ）。
【００２８】
Ａ）〜Ｄ）をもとにタンパク質中の活性部位内におけるリガンドの初期座標を利用し、以下の規則により空間点を作成することにした。
【００２９】
１）  「Ａ原子」がドナーまたは金属及び水の場合（ＳＡＲ情報の活性部位側の指定が水素結合ドナー、金属原子の場合）、１．で発生させたダミー原子の方向に対して「Ａ原子」から「Ｄ距離」の位置及びその周囲を初期座標に選んだ（図４、図５）。
【００３０】
２）  「Ａ原子」がｓｐ3軌道原子の場合（ＳＡＲ情報の活性部位側の指定がｓｐ3軌道原子の場合）には、「Ａ原子」から「Ｄ距離」の周囲を初期座標に選んだ（図６）。
【００３１】
３）  「Ａ原子」が水素結合アクセプター（以後、アクセプター）の場合（ＳＡＲ情報の活性部位側の指定が水素結合アクセプターの場合）、「Ａ原子」の結合延長上「Ｄ距離」の位置及びその周囲を初期座標に選んだ（図７）。
【００３２】
４）  その他の場合（ＳＡＲ情報の活性部位側の指定がその他の原子の場合）には、「Ａ原子」を中心とする半径が「Ｄ距離」の球表面上の点を初期座標に選んだ（図８）。
【００３３】
５）  （１）〜（４）とは、異なりリガンドの初期座標を直接指定することもできるようにした。
【００３４】
（２）  （０）で選択したリガンド中の原子と（１）で指定した空間点とのペアを重複がないようにランダムに選択する。
【００３５】
（３）  以下のスコアＳｓｃｏｒｅ（ｉ，ｊ）を計算する。
【数１】
2005242493.tif 000003
ｄijsは該当タンパク質中のｉ番目とｊ番目の空間点距離。ｄijCは化合物中のｉ番目とｊ番目の原子間距離。αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数。βは重なりと定義できる限界値を与えるための定数。
【００３６】
αは１．５、βは０．８とするのが好ましい。
【００３７】
（４）  （３）のスコアが最大になるように調整する。スコアを最大にする手法としては、例えば、シミュレーティッドアニーリング法が挙げられる。または、時間短縮には（２）、（３）を１００００回繰り返し、Ｓｓｃｏｒｅ（ｉ，ｊ）が最大になるペアを探し、そのペア情報をもとにリガンドを初期座標に重ね合わせる方法を適応することが好ましい。
【００３８】
（５）  （４）で重ねあわせたリガンドに対してタンパク質との相互作用エネルギーをコンフォメーションを微調整しながら最適化計算する。リガンドのコンフォメーションの微調整は、（４）で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３を越えない程度に座標変化させる。
【００３９】
微調整は例えばランダムサーチで最適化することが好ましい。ランダムサーチは以下の項目に従ってタンパク質の活性部位とリガンドとの微小変化を８０００回行い、最適エネルギー「Ｕ最適」が最小になるようにする。
【００４０】
１〕  回転可能な結合のうち最大５つ乱数で選び、結合ごとにランダムに±１０．０°の範囲内で回転させリガンドのコンフォメーションを換える。この過程を３回に一度行う。
【００４１】
２〕  ｘ、ｙ、ｚ軸方向それぞれにおいて、ランダムに±１．０Åの範囲内でリガンドの並進運動を行う。この過程を２回に一度行う。
【００４２】
３〕  回転中心座標それぞれにおいて、ランダムに±１．０Åの範囲内で回転中心座標を移動させ、さらに３次元方向の角度それぞれに対して、ランダムに±５．０°の範囲内でリガンドの回転運動を行う。この過程を５回に一度行う。
【００４３】
（６）  リガンドのコンフォメーションを大きく動かして、（２）から再スタートを行い、（５）までを繰り返して最適化を行う。コンフォメーション改変は、（３）で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３以上になるよう座標変化させる。
【００４４】
リガンドのコンフォメーションを大きく動かした最適化は、例えば（５）で最適化したエネルギー「Ｕ最適」でのコンフォメーションに対して、回転可能な結合をランダムに５つ選び、原子タイプごとに決められた回転角度間隔に従ってランダムに回転させる。その後（２）、（３）の過程を、５０００回繰り返し行う。
【００４５】
ただし、リガンドのコンフォメーションを変化させた後、リガンドの内部エネルギー「Ｕ内部」を計算しその値が５００．０以上のときはその後の計算をスキップし、次のリガンドコンフォメーションを発生させるようにする。
【００４６】
（７）  （１）〜（６）までの過程を（０）で用意した複数の構造変化座標に対して行い、最適なタンパク質とリガンドとの複合体座標、最適エネルギー「Ｕ最適」を算出する。
【００４７】
（８）  （１）〜（７）までの過程を（０）で用意した化合物データベース中の全てのリガンドに対して行い、化合物データベース中から該当タンパク質と結合する可能性のあるリガンドを選択する。
【００４８】
［３］  ［１］〜［２］によって記載された方法において、タンパク質の誘導適合を反映するパラメータおよび構造変化した立体構造座標を分子動力学計算方法を用いて算出する場合、該当タンパク質の立体構造に際し、基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、そのゆらぎの大きさを拘束条件として、分子動力学計算を行うことで、タンパク質の立体構造をエネルギー最適構造よりおおきく離れないようにして分子動力学計算を行う方法。
【００４９】
本手法による分子動力学計算は、例えば、基準振動計算より主鎖原子の２面角のゆらぎの値を算出し、該当のゆらぎ値を以下のように分子動力学計算における力の定数Ｋの部分に入れる。
Ｅｒｏｔ＝Ｋｒｏｔ（φ―φ０）2
Ｅｒｏｔはタンパク質の立体構造中において主鎖原子の２面角のエネルギーのことを示す。
φは主鎖原子の２面角。
φ０は主鎖原子の２面角の標準値。
Ｋｒｏｔの値が大きい場合はφはφ０に拘束される。
Ｅｐｏｓ＝Ｋｐｏｓ（ｒ−ｒ0）2
Ｅｐｏｓはタンパク質の立体構造中において主鎖原子の位置のエネルギーのことを示す。
ｒは主鎖原子の座標。
ｒ０は主鎖原子の座標の標準値。
Ｋｐｏｓの値が大きい場合はｒはｒ０に拘束される。
【００５０】
[４]  ［１］〜［２］によって記載された方法において、リガンドとタンパク質との相互作用を評価する際の目的関数として、従来の相互作用エネルギー関数に、タンパク質の動的性質を表現する関数を弾性エネルギーとして加え、タンパク質の立体構造座標から相互作用エネルギーを高速に算出するとともに、タンパク質の動的挙動に関する物理化学的性質を明確に描写することを特徴とする手法。
【００５１】
弾性エネルギーとして、タンパク質の局所的な柔らかさを考慮し、以下の関数「Ｕ衝突」として適応する。活性部位における動的挙動の少ない側鎖原子及び主鎖原子のみのｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義した。
【数２】
2005242493.tif 000004
ψ（ｉ，ｊ）＝Ｋ衝突  ＊  （Ｒ衝突（ｉ，ｊ）−Ｒ）2
Ｍは衝突を不許可とする活性部位の原子の数、Ｎはリガンドの原子の数、「Ｋ衝突」は１０００．０であることが好ましい。「Ｒ衝突（ｉ，ｊ）」は活性部位のｉ番目の原子とリガンドのｊ番目の原子それぞれのＶａｎ  ｄｅｒ  Ｗａａｌｓ半径の和とした。
【００５２】
ここで、活性部位の各原子に対し衝突を許す重み付けｗ（ｉ）が定義された場合以下の式を用いる。ただし、ｗ（ｉ）は、０〜１の範囲の実数とする。
【数３】
2005242493.tif 000005
ψ（ｉ，ｊ）＝ｗ（ｉ）  ＊  Ｋ衝突  ＊  （Ｒ衝突（ｉ，ｊ）−Ｒ）2
Ｍは活性部位の原子の数、Ｎはリガンドの原子の数、「Ｋ衝突」は１０００．０であることが好ましい。「Ｒ衝突（ｉ，ｊ）」は活性部位のｉ番目の原子とリガンドのｊ番目の原子それぞれのＶａｎ  ｄｅｒ  Ｗａａｌｓ半径の和とした。
【００５３】
また、弾性エネルギーとしては、以下の関数を用いて定義することも可能である。
Ｅｖ＝ｗ（ｈａｒｄ  ｓｈａｐｅ  ｒｅｇｉｏｎ）、Ｅ＝０（ｓｏｆｔ  ｓｈａｐｅ  ｒｅｇｉｏｎ）
ｈａｒｄ  ｓｈａｐｅ  ｒｅｇｉｏｎとは、タンパク質の立体構造中、動的挙動の小さい部分であり、ｓｏｆｔ  ｓｈａｐｅ  ｒｅｇｉｏｎとは、動的挙動の大きい部分のことを指す。Ｗは定数で１００であることが好ましい。
【００５４】
[５]  [１]〜[４]によって記載された方法において、タンパク質の動的性質を表現する関数として基準振動解析結果またはタンパク質の二次構造判定結果を用いることを特徴とした方法。
【００５５】
二次構造判定においては、タンパク質のへリックス、シート部分は揺らぎは小さいと考え、それ以外は揺らぎは大きいと考え相互作用の評価関数、分子動力学計算の拘束条件に適応する。
【００５６】
［６］  ［１］〜［５］によって記載された方法において、計算された該当タンパク質が代表的な複数の立体構造座標である場合や、該当タンパク質の立体構造が例えば核磁気共鳴スペクトルの解析結果のような複数の立体構造座標である場合についても、該当タンパク質と結合するリガンドの探索について、複数座標すべてを全自動的にかつ短時間で同等に評価することを可能とした方法。
【発明の効果】
【００５７】
本発明によれば、本発明はタンパク質の立体構造座標を用いたリガンド探索方法、詳しくはタンパク質立体構造座標が既知の場合、相互作用すると考えられるリガンドを予測する方法に関して、タンパク質とリガンドとの結合に大変重要であるタンパク質の動的挙動を反映したパラメータを取得し、かつタンパク質の動的挙動を反映したリガンドとの新規な相互作用評価関数を用いて、該当タンパク質の立体構造と結合する新規リガンドを予測を行うことができる。これにより従来法と比較して、より信頼性の高い、かつ医薬品設計等に適したタンパク質の立体構造を世界中で解析されている大量のゲノム配列に関しても対応するスピードで構築することができる。従来はｉｎ  ｓｉｌｉｃｏ  スクリーニングにおいては、タンパク質とリガンドとの相互作用に重要な誘導結合を充分に取り扱うことのできるアルゴリズムが見出されていなかった時点に対し、タンパク質とリガンドとの相互作用エネルギー関数に、基準振動計算結果、もしくは二次構造予測から得られるタンパク質の「ゆらぎ」を表すパラメータを簡易にとりこむ計算式を導入した。
【００５８】
さらに、分子動力学シミュレーションにおいては、この方法により、該当タンパク質の動的挙動を反映したパラメータとリガンドとの相互作用評価関数に関して、該当タンパク質についての基準振動計算を行い、その結果を分子動力学計算に反映させることを特徴とする。従来はタンパク質の動的挙動のシミュレーションを行うためには、分子動力学計算を用いていたが、従来法で該当タンパク質に対して分子動力学計算を行うとタンパク質立体構造はＸ線、ＮＭＲ等で解析された座標から大きくズレる。こうしたズレはタンパク質の動的挙動の物理化学的描写を含んでいるがＮＭＲ等で示される動的挙動の実験的な結果と矛盾する挙動となる場合があり、必ずしも精度の高いシミュレーションとならないことが多い。そこで、分子動力学計算を行う際には、タンパク質の立体構造をある程度固定しシミュレーションを行う必要があり、本手法では分子動力学計算におけるエネルギー関数中で主鎖原子の２面角に拘束をかける手法を開発した。さらに２面角の拘束条件としては、そのパラメータとして予め該当タンパク質の基準振動計算を行い、主鎖原子の２面角のゆらぎを算出し、そのゆらぎの大きさにより例えばゆらぎの大きい部分は拘束条件を緩め、ゆらぎの小さい部分は拘束条件を強めるパラメータとして用いることとした。こうした条件でタンパク質の分子動力学シミュレーションを行うことで、精度よく動的挙動を描写することができる。加えてこうして算出された分子シミュレーションからタンパク質の動的挙動を描写した座標を取得することができ、これを利用することでさまざまなリガンド結合部位の形状を用いたリガンド探索を行うことができる。
【００５９】
これらの結果、今までのｉｎ  ｓｉｌｉｃｏ  スクリーニングでは見出すことができなかった新規なリガンドを発見することを可能するとともに、今までは長時間を必要とする分子シミュレーションでしか解析できなかった「誘導結合」を含めたタンパク質−リガンドとの相互作用解析を短時間で行うことを可能にした。
【００６０】
こうしたアルゴリズムは、既存ソフトウェアよりも誘導結合現象をより深く考慮した“ｉｎ  ｓｉｌｉｃｏ  ｓｃｒｅｅｎｉｎｇ”に対応可能とし、誘導結合現象と疎水相互作用の正しい理解のもと単純化している。アルゴリズムは単純化されているので、自動化により多くのターゲットタンパク質を処理可能とする。その結果、例えば１００万以上の化合物データベースから、新規で、もっともらしい化合物を探索することができるので実験では対応できない規模のデータベースからもっともらしい化合物を現実的な時間内に探索することが出来る。
【００６１】
また、タンパク質−リガンドとの相互作用解析を短時間で行うことが可能になるので、例えば代謝、毒性の原因となる数多くのタンパク質と薬物との相互作用解析が可能となり、ｉｎ  ｓｉｌｉｃｏでの薬物の代謝、毒性予測を行うことができる。
【００６２】
本発明において、リガンドとして取り扱うことのできる分子は使用するリガンドの種類や数を限定しないため、蛋白質、ペプチド、ＤＮＡ、薬剤成分、金属、イオン、糖類、核酸成分、ホルモンを含む全ての物質を当該リガンドと見なすことができる。この方法によって、具体的に農薬、医薬品等の分子設計を行うことができる。
【００６３】
リガンドとタンパク質との相互作用エネルギー評価関数には、従来ドッキング法では静電エネルギー項、ｖａｎ  ｄｅｌ  ｗａａｌｓ項、さらにはソフトドッキング法等に見られてる動的挙動を表現するための調整項が主に用いられているが、本手法においてはタンパク質とリガンドとの相互作用中にはソフトドッキング法等に見られてる動的挙動を表現するための調整項を用いる代わりに古典力学で用いられている弾性衝突の理論を適応し、タンパク質とリガンドとの相互作用に関して、その物理化学的性質をより明確にした。このことによりタンパク質の構造変化と相互作用との関係を得ることができ、リガンドの機能の理解を迅速かつ正確に行うための手助けとなる。
【００６４】
尚、本発明で利用するタンパク質の立体構造は、その際Ｘ線結晶構造解析等により、タンパク質の立体構造として３次元座標が決定されたもの以外に、タンパク質の経験的なモデリング法、特にホモロジーモデリング法或いはスレッディング法を利用して作成した立体構造座標をも適応することができる。
【発明を実施するための最良の形態】
【００６５】
以下、本発明の実施の形態について説明する。好適かつ代表的な例について説明するものであり、本発明はこれに限定されるものではない。
【００６６】
本明細書において幾つかの用語を使用するが、特に明記しない限り、次の意味を有する。
【００６７】
「標的タンパク質」とは、立体構造の詳細がＸ線結晶解析やＮＭＲ解析、ホモロジーモデリング法により既に決定されており、リガンド探索の対象とするタンパク質を意味する。
【００６８】
「原子座標」とは、三次元空間上で立体構造を記述するものである。それは空間上のある点を原点とする互いに垂直な三方向の相対的な距離であり、タンパク質中に存在する水素原子を除く原子一つあたりに３個の数字からなるベクトル量である。
【００６９】
図１は、本発明によるタンパク質立体構造と誘導適合を利用したリガンド探索方法の一例を示すフローチャートである。図１は、本発明によるタンパク質の立体構造予測を示すフローチャート。
【００７０】
図１に示す通り、この方法においては、先ず、ステップ０において、リガンドのデータベースを用意する。リガンドのデータベースは分子力学法等に用いて３次元化することが望ましい。
【００７１】
ステップ１０において標的タンパク質の立体構造を選定し、座標を入手する。
【００７２】
ステップ２０において、標的タンパク質の基準振動計算を行い、主鎖原子の位置のゆらぎと２面角のゆらぎを求める。
【００７３】
ステップ３０において、ステップ２０において求めた標的タンパク質のゆらぎを拘束条件として用いた分子動力学計算を行う。
【００７４】
ステップ４０において、標的タンパク質のリガンド結合部位に、リガンドを配置するための点郡を指定する。
【００７５】
ステップ５０において、ステップ０で定められた一つのリガンドに対し、リガンドの各原子をステップ３０で定められた点群に重ね合わせる。
【００７６】
ステップ６０において、ステップ５０で定められた重ねあわせに対し、ステップ２０およびステップ３０で定められた計算結果によりタンパク質の動的挙動を表すパラメータを取得し、該当パラメータを用いてリガンドとタンパク質との相互作用エネルギーをリガンドコンフォメーションを微調整しながら計算する。コンフォメーションの微調整は、ステップ５０で定められたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３を越えない程度に座標変化させる。
【００７７】
ステップ７０において、ステップ５０で定められたリガンドに対して、コンフォメーションを大きく動かして、ステップ５０から再スタートを行い、ステップ７０までを繰り返して最適化を行う。コンフォメーション改変は、ステップ５０で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３以上になるよう座標変化させる。
【００７８】
ステップ８０において、ステップ７０まで得られた標的タンパク質とリガンドとの相互作用エネルギーを決定する。
【００７９】
ステップ９０において、ステップ４０に戻り、ステップ０中の別のリガンドを選択し、ステップ８０まで計算する。
【００８０】
ステップ１００において、ステップ０中でのリガンドに対し、ステップ９０において定められた相互作用エネルギーを比較し、標的タンパク質に結合すると予想されるリガンドを選択する。
【００８１】
本発明では、従来は標的タンパク質の立体構造が与えられた場合におけるリガンド探索が、タンパク質とリガンドとの相互作用にとって重要なタンパク質の動的な性質を反映させることが困難であった時点に対し、ステップ２０から９０までを行うことで、タンパク質とリガンドとの相互作用において、タンパク質の動的性質に関する運動エネルギー部分を弾性率として扱い、リガンドの近接に伴うタンパク質の弾性的性質を取り込んだエネルギー評価関数（安定化に寄与する疑似弾性エネルギー）を用いることで、相互作用エネルギーをｖａｎ  ｄｅｌ  ｗａａｌｔｈ相互作用エネルギーと疑似弾性エネルギーとの和の形式とし、物理化学的現象をはっきりさせた評価関数を用いることが従来法と異なる点である。以下、各ステップについて更に詳細に説明する。
【００８２】
（ステップ０：リガンドデータベースの準備）
３次元座標を含むリガンドデータベースを用意する。リガンドデータベースとしては、例えば、ＡＣＤ等のような市販化合物データベース、化合物を描いて収集した仮想化合物データを用いることができる。ガンドのデータベースは分子力学法等に用いて３次元化することが望ましい。
【００８３】
（ステップ１０：標的タンパク質立体構造の選択と取得）
ステップ０で定められたリガンドデータベースから特定リガンドを探索するための標的タンパク質を選択し、３次元座標を入手する。３次元座標は、公共データベースであるＰＤＢやホモロジーモデリング法等で作成した立体構造座標を用いることが望ましい。
【００８４】
（ステップ２０：標的タンパク質の基準振動計算）
ステップ１０で定められた参照タンパク質の動的挙動を表すパラメータを基準振動解析法による計算結果のデータベースもしくは二次構造判定計算をおこない取得する。まず、基準振動解析法によるタンパク質の動的挙動を表すパラメータ取得方法について下記に示す。
【００８５】
基準振動解析法とは、ポテンシャルエネルギーを変位の二次関数として近似し、運動方程式を厳密に解き、最適化構造の周りの微小な振動を解析する方法を意味する。解くべき運動方程式は下記式（１）または（２）である。
【数４】
2005242493.tif 000006
・・・（１）
【数５】
2005242493.tif 000007
・・・（２）
ただし、
【数６】
2005242493.tif 000008
【００８６】
ここでωkは固有値、Ｕikは固有ベクトルであり、δijはクロネッカーのデルタである。ＴijとＶijはそれぞれ運動エネルギーＥkとポテンシャルエネルギーＶに関係し、下記式（３）および（４）の通りである。
【数７】
2005242493.tif 000009
・・・（３）
【数８】
2005242493.tif 000010
・・・（４）
【００８７】
ここで、ｑiは振動の自由度に対応した座標、ｑi0は最適化座標、ｑi’（式（３）における「ｑiドット」を意味する）はｑiの時間による微分である。Ａjkは集団運動Ｑkと個々の原子運動ｑjを結ぶ係数であり、下記式（５）の通りである。
【数９】
2005242493.tif 000011
・・・（５）
ただし、基準振動
【数１０】
2005242493.tif 000012
である。ここで、αkとδkは初期条件で定められる。
【００８８】
上記した基準振動解析法の詳細は「Ｗｉｌｓｏｎ，Ｅ．Ｂ．，Ｄｅｃｉｕｓ，Ｊ．Ｃ．，ａｎｄ  Ｃｒｏｓｓ，Ｐ．  Ｃ．  １９５５．  Ｍｏｌｅｃｕｌａｒ  Ｖｉｂｒａｔｉｏｎ．  ＭｃＧｒａｗ−Ｈｉｌｌ．」に記載されている。
【００８９】
参照タンパク質に対して、上記で得られた固有値、固有ベクトルを用いて、ある温度・ある固有値での各Ｃα原子の位置ゆらぎを計算し、このゆらぎの値をＣαが含まれるアミノ酸のゆらぎの値とする。目的タンパク質の各アミノ酸のゆらぎの値は、ステップ５０におけるアライメントを利用して、目的配列と参照配列の比較から対応するアミノ酸残基ペアにおいて、目的タンパク質のゆらぎの値として参照タンパク質と同一のものを当てはめておく。ゆらぎの値を求められなかったものについては、予め設定しておいた値をあてはめる。こうして得た目的タンパク質の各アミノ酸のゆらぎの値を目的タンパク質の動的な挙動を表すパラメータとする。
【００９０】
次に、二次構造判定計算によるタンパク質の動的挙動を表すパラメータ取得方法について下記に示す。
【００９１】
二次構造判定はタンパク質の立体構造座標から計算される。ソフトウェアとしては、ＤＳＳＰ、ＳＴＲＩＤＥ等が好ましいが、基本的にはタンパク質の主鎖のねじれ角と水素結合パターンから判別される方法を用いる。
【００９２】
ここで、「ＤＳＳＰ（Ｄｉｃｔｉｏｎａｒｙ  ｏｆ  ｐｒｏｔｅｉｎ  ｓｅｃｏｎｄａｒｙ  ｓｔｒｕｃｔｕｒｅ  ｏｆ  ｐｒｏｔｅｉｎ）」とは、ＰＤＢ書式のファイルを入力ファイルとして、主鎖の水素結合パターンと、内部回転角等を解析しαへリックスとβシートとを判定するソフトウェアである。ＤＳＳＰの詳細は、「Ｋａｂｓｃｈ，Ｗ．  ＆  Ｓａｎｄｅｒ，Ｃ．  （１９８３）  Ｄｉｃｔｉｏｎａｒｙ  ｏｆ  ｐｒｏｔｅｉｎ  ｓｅｃｏｎｄａｒｙ  ｓｔｒｕｃｔｕｒｅ：ｐａｔｔｅｒｎ  ｒｅｃｏｇｎｉｔｉｏｎ  ｏｆ  ｈｙｄｒｏｇｅｎ−ｂｏｎｄｅｄ  ａｎｄ  ｇｅｏｍｅｔｒｉｃａｌ  ｆｅａｔｕｒｅｓ．  Ｂｉｏｐｏｌｉｍｅｒｓ，  ２２：２５７７−２６３７」に記載されている。
【００９３】
「ＳＴＲＩＤＥ（Ｐｒｏｔｅｉｎ  ｓｅｃｏｎｄａｒｙ  ｓｔｒｕｃｔｕｒｅ  ａｓｓｉｇｎｍｅｎｔ  ｆｒｏｍ  ａｔｏｍｉｃ  ｃｏｏｒｄｉｎａｔｅ）」とは、ＰＤＢ書式のファイルを入力ファイルとして、主鎖の水素結合パターンと、内部回転角等を解析しαへリックスとβシートとを判定するソフトウェアである。ＳＴＲＩＤＥの詳細は、「Ｆｒｉｓｈｍａｎ，  Ｄ  ＆  Ａｒｇｏｓ，Ｐ．（１９９５）  Ｋｎｏｗｌｅｄｇｅ−ｂａｓｅｄ  ｓｅｃｏｎｄａｒｙ  ｓｔｒｕｃｔｕｒｅ  ａｓｓｉｇｎｍｅｎｔ．  Ｐｒｏｔｅｉｎｓ：  ｓｔｒｕｃｔｕｒｅ，  ｆｕｎｃｔｉｏｎ  ａｎｄ  ｇｅｎｅｔｉｃｓ，  ２３，  ５６６−５７９」に記載されている。
【００９４】
参照タンパク質に対して、上記ソフトウェア等を用いて、二次構造計算を行い、各アミノ酸がとるαへリックス構造、βシート構造、ループ構造を判定する。目的タンパク質の各アミノ酸の二次構造は、ステップ５０におけるアライメントを利用して、目的配列と参照配列の比較から対応するアミノ酸残基ペアにおいて、目的タンパク質の二次構造判定として参照タンパク質と同一のものを当てはめておく。二次構造判定を求められなかったものについては、予め設定しておいた結果をあてはめる。こうして得た目的タンパク質の各アミノ酸の二次構造判定結果を目的タンパク質の動的な挙動を表すパラメータとする。
【００９５】
上記の目的タンパク質の動的挙動を表すパラメータとしては、参照タンパク質の基準振動解析法より取得した結果を用いることが好ましく、該当計算結果は別途データベースとして保存されているものを使用する。二次構造判定計算結果は、好ましくは、基準振動解析が行われていない参照タンパク質を用いる際に基準振動解析計算の代用として使用する。
【００９６】
（ステップ３０：標的タンパク質の分子動力学計算（拘束つき））
主鎖の位置拘束エネルギー「Ｕ位置」を導入し、初期の受容体骨格の変動を抑えながらＡＰＲＩＣＯＴ［Ｙｏｎｅｄａ  Ｓ．  ＆  Ｕｍｅｙａｍａ  Ｈ．  （１９９２）  Ｆｒｅｅ  ｅｎｅｒｇｙ  ｐｅｒｔｕｒｂａｔｉｏｎ  ｃａｌｃｕｌａｔｉｏｎｓ  ｏｎ  ｍｕｌｔｉｐｌｅ  ｍｕｔａｔｉｏｎ  ｂａｓｅｓ  Ｊ．  Ｃｈｅｍ．  Ｐｈｙｓ．  ９７，  ６７３０−６７３６］を用いて最小化（条件：温度３００Ｋ、受容体の表面から水分子が最低２分子配置できる箱状水槽、力場：ＡＭＢＥＲ［Ｓ．  Ｊ．  Ｗｅｉｎｅｒ，  Ｐ．Ａ．  Ｋｏｌｌｍａｎ，  Ｄ．Ａ．  Ｃａｓｅ，  Ｕ．Ｃ．  Ｓｉｎｇｈ，  Ｃ．  Ｇｈｉｏ，  Ｇ．  Ａｌａｇｏｎａ，  Ｓ．  Ｐｒｏｆｅｔａ，  ＆，  Ｐ．  Ｗｅｉｎｅｒ  （１９８４）  Ａ  ｎｅｗ  ｆｏｒｃｅ  ｆｉｅｌｄ  ｆｏｒ  ｍｏｌｅｃｕｌａｒ  ｍｅｃｈａｎｉｃａｌ  ｓｉｍｕｌａｔｉｏｎ  ｏｆ  ｎｕｃｌｅｉｃ  ａｃｉｄｓ  ａｎｄ  ｐｒｏｔｅｉｎｓ  Ｊ．  Ａｍ．  Ｃｈｅｍ．  Ｓｏｃ．  １０６，  ７６５−７８４］）を行った。
Ｕ位置  ＝  Ｋ位置 ＊  Ｒ2
ここで、「Ｋ位置」は３００．０、Ｒは基準座標からのずれとした。
【００９７】
続いて、ＡＰＲＩＣＯＴに二面角拘束エネルギー「Ｕ二面角」を導入して、最小化した受容体のＭＤ計算（条件：温度３００Ｋ、受容体の表面から水分子が最低２分子配置できる箱状水槽、力場：ＡＭＢＥＲ）を行った。
Ｕ二面角  ＝  Ｋ二面角  ＊  （θ−θ平衡）2
θは二面角（単位ｒａｄ）。
【００９８】
「Ｋ二面角」には、最大値と最小値を指定することで、その範囲内で主鎖二面角揺らぎに対応するように各二面角に対して不均一な拘束がかかるようにした。以後、主鎖二面角を拘束しながら行うＭＤを二面角拘束ＭＤと呼ぶことにする。
【００９９】
二面角拘束ＭＤ計算によりタンパク質構造座標を入手するには、受容体動的構造のクラスタリングを行う。
【０１００】
あらかじめ指定した活性部位に対して、ＭＤの途中経過１００ｆｓｅｃごとの受容体を重ね合わせた構造及び初期構造の活性部位を母集団とした。まず初めに、クラスタリングすることにより側鎖の動的情報が失われる可能性が高いことから、側鎖の二面角χにおいて母集団のα％が平均角度±２０．０°の範囲で保存されている全側鎖二面角を収集した。ただし、主鎖の根元に近い方からχが保存されていないと判定された場合はそれ以降のχは保存されていないものとした。
【０１０１】
次に、収集した保存側鎖二面角をすべて網羅している構造を母集団から抽出した。そして、抽出した構造の類似性を比較するために全原子ｒｍｓ（ｒｏｏｔ  ｍｅａｎ  ｓｑｕａｒｅ）がβÅ以下の場合、同一構造と判断して一方を削除、最終的に選ばれた構造をもとに受容体動的構造クラスターを作成した。また、保存されていなかった二面角χを構成する原子では、変動する可能性が高いことから活性部位・リガンド結合計算において衝突しても良いことにした。ただし、α、βは定数。
【０１０２】
（ステップ４０：基準となる標的タンパク質上の空間点を選択する。）
ステップ３０で作成した複数のタンパク質立体構造座標のうち１つをランダムに選択する。タンパク質座標中の空間点は例えば以下のような方法で発生させる。
【０１０３】
１．  ダミー原子の発生による空間点の発生
リガンドとタンパク質との相互作用における水素結合に着目し、タンパク質中の水素結合サイトを空間点として指定する。水素結合における重要事項は距離と角度である。つまり、角度を計算するためには、水素結合ドナー（以後、ドナー）に水素原子が必要になる。そこで、活性部位及びリガンドに水素原子が含まれていない場合、以下の規則によりダミー水素原子を発生させた。
【０１０４】
１）  ｓｐ2軌道原子を中心とする正三角形状にダミー原子を発生（図２）。
【０１０５】
２）  ｓｐ3軌道原子では、水素結合を形成する距離にある場合、水素原子を共有するように回転できると考え水素結合相互作用を計算するときには、距離のみを考慮することにした。このため、ｓｐ3軌道原子にはダミー原子を発生させない。
【０１０６】
金属及び水では、活性部位・リガンド結合の仲介役となり得るので、相互作用する位置に以下のようにダミー原子を発生させた。
【０１０７】
１）  鉄のような金属には、正八面体状にダミー原子を発生（図３）。
【０１０８】
２）  水は、正四面体状にダミー原子を発生。
【０１０９】
ただし、活性部位と相互作用している方向にはダミー原子を発生させないことにした。
【０１１０】
２．  構造活性相関情報を利用した空間点の発生
リガンドの構造活性相関（ＳＡＲ）情報に着目し、以下の項目を入力情報にすることにした。
【０１１１】
Ａ）  ＳＡＲから得られた活性部位の原子（以後、「Ａ原子」）。ＰＤＢ形式に従う。
【０１１２】
Ｂ）  「Ａ原子」と相互作用するであろうリガンドの原子タイプ（以後、「Ｂタイプ」）。ＳＹＢＹＬのＭＯＬ２形式に従う。
【０１１３】
Ｃ）  「Ａ原子」と「Ｂタイプ」との相互作用の強さ（以後、「Ｃ強さ」）。
【０１１４】
Ｄ）  「Ａ原子」と「Ｂタイプ」との相互作用する距離（以後、「Ｄ距離」）（単位はÅ）。
【０１１５】
Ａ）〜Ｄ）をもとにタンパク質中の活性部位内におけるリガンドの初期座標を利用し、以下の規則により空間点を作成することにした。
【０１１６】
１）  「Ａ原子」がドナーまたは金属及び水の場合、１．で発生させたダミー原子の方向に対して「Ａ原子」から「Ｄ距離」の位置及びその周囲を初期座標に選んだ（図４、図５）。
【０１１７】
２）  「Ａ原子」がｓｐ3軌道原子の場合には、「Ａ原子」から「Ｄ距離」の周囲を初期座標に選んだ（図６）。
【０１１８】
３）  「Ａ原子」が水素結合アクセプター（以後、アクセプター）の場合、「Ａ原子」の結合延長上「Ｄ距離」の位置及びその周囲を初期座標に選んだ（図７）。
【０１１９】
４）  その他の場合には、「Ａ原子」を中心とする半径が「Ｄ距離」の球表面上の点を初期座標に選んだ（図８）。
【０１２０】
５）  （１）〜（４）とは、異なりリガンドの初期座標を直接指定することもできるようにした。
【０１２１】
（ステップ５０：ステップ４０の点群にステップ０中でリガンド座標の一つを重ね合わせる。）
距離行列を用いたアライメント作成アルゴリズム（ＤＡＬＩ）［Ｈｏｌｍ，  Ｌ．，  ＆  Ｓａｎｄｅｒ，  Ｃ．  （１９９３）  Ｐｒｏｔｅｉｎ  Ｓｔｒｕｃｔｕｒｅ  Ｃｏｍｐａｒｉｓｏｎ  ｂｙ  Ａｌｉｇｎｍｅｎｔ  ｏｆ  Ｄｉｓｔａｎｃｅ  Ｍａｔｒｉｃｅｓ  Ｊ．  Ｍｏｌ．  Ｂｉｏｌ．  ２３３，  １２３−１３８］を低分子用に改良した手法で初期座標とリガンドを重ね合わせた。
【０１２２】
（１）  １つの「Ｂタイプ」には、リガンドの原子タイプが複数対応することが多い。そこで、乱数を用いて「Ｂタイプ」とリガンドの原子タイプで同一視できるペアを作成した。ただし、ペアにおいてリガンドの原子タイプが重複しないようにした。
【０１２３】
（２）  「Ｂタイプ」には、ステップ４０−２．により複数の初期座標が含まれているので、初期座標も乱数を用いて選択した。
【０１２４】
（３）  選択された初期座標とリガンドそれぞれの距離行列を作成し、Ｓｓｃｏｒｅ（ｉ，ｊ）を計算する。
【数１１】
2005242493.tif 000013
ｄijsは該当タンパク質中のｉ番目とｊ番目の空間点距離。ｄijCは化合物中のｉ番目とｊ番目の原子間距離。αは、該当タンパク質中の空間点群と化合物が完全に重なりあった場合にＳｓｃｏｒｅ（ｉ，ｊ）を最大値とするための定数。βは重なりと定義できる限界値を与えるための定数。
【０１２５】
αは１．５、βは０．８とするのが好ましい。
【０１２６】
（４）  （１）〜（３）を１００００回繰り返し、Ｓｓｃｏｒｅ（ｉ，ｊ）が最大になるペアを探し、そのペア情報をもとにリガンドを初期座標に重ね合わせた。
【０１２７】
（ステップ６０：標的タンパク質との相互作用関数をリガンドのコンフォメーションを微調整しながら最適化する。）
ステップ５０で重ねあわせたリガンドに対してタンパク質との相互作用エネルギーをコンフォメーションを微調整しながら最適化計算する。リガンドのコンフォメーションの微調整は、ステップ５０で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３を越えない程度に座標変化させる。
【０１２８】
微調整は例えばランダムサーチで最適化することが好ましい。ランダムサーチは以下の項目に従ってタンパク質の活性部位とリガンドとの微小変化を８０００回行い、最適エネルギー「Ｕ最適」が最小になるようにする。
【０１２９】
１〕  回転可能な結合のうち最大５つ乱数で選び、結合ごとにランダムに±１０．０°の範囲内で回転させリガンドのコンフォメーションを換える。この過程を３回に一度行う。
２〕  ｘ、ｙ、ｚ軸方向それぞれにおいて、ランダムに±１．０Åの範囲内でリガンドの並進運動を行う。この過程を２回に一度行う。
【０１３０】
回転中心座標それぞれにおいて、ランダムに±１．０Åの範囲内で回転中心座標を移動させ、さらに３次元方向の角度それぞれに対して、ランダムに±５．０°の範囲内でリガンドの回転運動を行う。この過程を５回に一度行う。
【０１３１】
最適エネルギーは以下のように定義する。
Ｕ最適＝ＵSAR＋Ｕ水素＋Ｕ疎水＋Ｕスタッキング＋Ｕ衝突＋Ｕ内部
【０１３２】
ここで、原子のＶａｎ  ｄｅｒ  Ｗａａｌｓ半径及び原子間相互作用距離はＡＭＢＥＲ９９［Ｊ．  Ｗａｎｇ，  Ｐ．  Ｃｉｅｐｌａｋ  ＆  Ｐ．Ａ．  Ｋｏｌｌａｍ  （２０００）  Ｈｏｗ  ｗｅｌｌ  ｄｏｅｓ  ａ  ｒｅｓｔｒａｉｎｅｄ  ｅｌｅｃｔｒｏｓｔａｔｉｃ  ｐｏｔｅｎｔｉａｌ  （ＲＥＳＰ）  ｍｏｄｅｌ  ｐｅｒｆｏｒｍ  ｉｎ  ｃａｌｃｕｌａｔｉｎｇ  ｃｏｎｆｏｒｍａｔｉｏｎａｌ  ｅｎｅｒｇｉｅｓ  ｏｆ  ｏｒｇａｎｉｃ  ａｎｄ  ｂｉｏｌｏｇｉｃａｌ  ｍｏｌｅｃｕｌｅｓ？  Ｊ．  Ｃｏｍｐｕｔ．  Ｃｈｅｍ．  ２１，  １０４９−１０７４］およびＭＭ３パラメータ［Ｍａ  Ｂ．，  Ｌｉｉ  Ｊ．−Ｈ．，  Ａｌｌｉｎｇｅｒ  Ｎ．Ｌ．  （２０００）  Ｍｏｌｅｃｕｌａｒ  ｐｏｌａｒｉｚａｂｉｌｉｔｉｅｓ  ａｎｄ  ｉｎｄｕｃｅｄ  ｄｉｐｏｌｅ  ｍｏｍｅｎｔｓ  ｉｎ  ｍｏｌｅｃｕｌａｒ  ｍｅｃｈａｎｉｃｓ  Ｊ．  Ｃｏｍｐｕｔ．  Ｃｈｅｍ．  ２１，  ８１３−８２５］を参考にした。
【０１３３】
（ａ）  ＳＡＲに関するエネルギー関数
ＳＡＲ情報に従う指標としてエネルギーＵSARを定義した。
【数１２】
2005242493.tif 000014
ψ（ｉ）＝ＫSAR（ｉ）  ＊  （ＲSAR（ｉ）−Ｒ）2−δ
ＮはＳＡＲ情報の数、Ｒは「Ａ原子」からリガンド側の相互作用原子までの距離、ＫSAR（ｉ）はｉ番目の「Ｃ強さ」、ＲSAR（ｉ）はｉ番目の「Ｄ距離」、δは２０．０とした。
【０１３４】
（ｂ）  水素結合に関するエネルギー関数
リガンドの１つのドナー（アクセプター）に対して１つだけ水素結合を形成すると考え、最短にある活性部位側のアクセプター（ドナー）を選び、水素を介した結合角θ（図９。ただし、複数の水素原子が付加されているドナー原子の場合には、最小の水素結合角をθと定義した。）を算出して、次の条件により分岐してψ（ｉ）を計算するように定義した。図９において、Ａはドナー、Ｂは水素、Ｃはアクセプター、θは水素結合角を示す。
【数１３】
2005242493.tif 000015
（１）  ドナー原子がｓｐ3軌道原子、または水素結合角θが±３０．０°以内のとき
【数１４】
2005242493.tif 000016
（２）  水素結合角θが±３０．０°以上のとき
【数１５】
2005242493.tif 000017
Ｎはリガンドのドナー＋アクセプターの数、Ｒは水素結合を形成する二原子間距離、「Ｋ水素（ｉ）」及び「Ｒ水素（ｉ）」は原子タイプごとに決めた水素結合の相互作用の強さ及び距離とした。
【０１３５】
（ｃ）  疎水相互作用エネルギー
活性部位（ＡＬＡ、ＣＹＳ、ＰＨＥ、ＩＬＥ、ＬＥＵ、ＭＥＴ、ＰＲＯ、ＶＡＬ、ＴＲＰ、ＴＹＲの側鎖。ただし、ＴＹＲの水酸基は除く）及びリガンド（炭素原子）の疎水相互作用し得る原子に通し番号を付け、活性部位のｉ番目とリガンドのｊ番目との原子間距離Ｒがカットオフ以内にあるときψ（ｉ，ｊ）を計算するように定義した。
【数１６】
2005242493.tif 000018
【数１７】
2005242493.tif 000019
Ｍは活性部位の疎水相互作用し得る原子の数、Ｎはリガンドの疎水相互作用し得る原子の数、「Ｋ疎水（ｉ，ｊ）」及び「Ｒ疎水（ｉ，ｊ）」は原子タイプごとに決めた疎水相互作用の強さ及び距離とした。また、カットオフは８．０Åとした。
【０１３６】
（ｄ）  スタッキングエネルギー
活性部位及びリガンドの芳香環を形成する原子に通し番号を付け、活性部位においては芳香環の中心座標を算出した。活性部位のｉ番目とリガンドのｊ番目との原子間距離Ｒがカットオフ以内にあるとき、ｉ番目の原子が形成する芳香環の中心座標をｉ’、ｊ番目の最短距離あり共に同じ芳香環を形成するリガンドの原子をｊ’としたとき、∠ｉｉ’ｊ＝θi'j、∠ｉ’ｉｊ＝θij、∠ｉｉ’ｊ’＝θi'j'、∠ｉ’ｉｊ’＝θij'を算出し（図１０）、θi'jとθijとが９０．０°±１０．０°のとき、「Ｒ境界」と「θ境界」を求め次の条件により分岐してψ（ｉ，ｊ）を計算するように定義した。図１０において、ｉ’は活性部位の芳香環中心、ｉは活性部位の芳香環原子、ｊ，ｊ’はリガンドの芳香環原子を示す。
【数１８】
2005242493.tif 000020
Ｒ境界＝１．０−（Ｒスタッキング（ｉ，ｊ）−Ｒ）2
【数１９】
2005242493.tif 000021
θ境界＝｜１．０−Θ｜
Ｉｆ  Ｒ境界＜０．０，
ψ（ｉ，ｊ）＝−Ｋスタッキング（ｉ，ｊ）  ＊  Ｒ境界
Ｅｌｓｅ，
ψ（ｉ，ｊ）＝−Ｋスタッキング（ｉ，ｊ）  ＊  θ境界
Ｍは活性部位の芳香環を形成する原子の数、Ｎはリガンドの芳香環を形成する原子の数、「Ｋスタッキング（ｉ，ｊ）」及び「Ｒスタッキング（ｉ，ｊ）」は原子タイプごとに決めたスタッキングの強さ及び距離とした。πは円周率、θはθi'j'とθij'においてΘが最小になる角度とした。また、カットオフは５．０Åとした。
【０１３７】
（ｅ）  分子間衝突（弾性衝突エネルギー）
活性部位（クラスタリングの際に定義した、保存されている側鎖原子及び主鎖原子のみ）のｉ番目の原子とリガンドのｊ番目の原子との原子間距離Ｒが衝突距離「Ｒ衝突（ｉ，ｊ）」以内のとき、ψ（ｉ，ｊ）を計算するように定義した。
【数２０】
2005242493.tif 000022
ψ（ｉ，ｊ）＝Ｋ衝突  ＊  （Ｒ衝突（ｉ，ｊ）−Ｒ）2
Ｍは衝突を不許可とする活性部位の原子の数、Ｎはリガンドの原子の数、「Ｋ衝突」は１０００．０とした。「Ｒ衝突（ｉ，ｊ）」は活性部位のｉ番目の原子とリガンドのｊ番目の原子それぞれのＶａｎ  ｄｅｒ  Ｗａａｌｓ半径の和とした。
【０１３８】
ここで、活性部位の各原子に対し衝突を許す重み付けｗ（ｉ）が定義された場合以下の式を用いる。ただし、ｗ（ｉ）は、０〜１の範囲の実数とする。
【数２１】
2005242493.tif 000023
ψ（ｉ，ｊ）＝ｗ（ｉ）  ＊  Ｋ衝突  ＊  （Ｒ衝突（ｉ，ｊ）−Ｒ）2
Ｍは活性部位の原子の数、Ｎはリガンドの原子の数、「Ｋ衝突」は１０００．０であることが好ましい。「Ｒ衝突（ｉ，ｊ）」は活性部位のｉ番目の原子とリガンドのｊ番目の原子それぞれのＶａｎ  ｄｅｒ  Ｗａａｌｓ半径の和とした。
【０１３９】
（ｆ）  リガンド内部エネルギー
回転可能な結合を微小に変動させていくと、誤差で結合が切れる恐れがあるために「ψ結合長（ｉ）」を、また、リガンド内部で原子衝突が起こることも避けるために「ψ衝突（ｉ，ｊ）」を計算するように定義した。
【数２２】
2005242493.tif 000024
ψ結合長（ｉ）＝Ｋ結合長＊｛１０００．０  ＊  （Ｒ結合長（ｉ）−Ｒ1）｝2
ψ衝突（ｉ，ｊ）＝Ｋ衝突＊（Ｒ衝突−Ｒ2）2
Ｌは回転可能な結合の数、Ｍはリガンドの原子数、Ｎはｉ番目の原子の非結合原子数、「Ｋ結合長」は１００．０、「Ｒ結合長（ｉ）」は初期構造の結合長、「Ｋ衝突」は１５０．０、「Ｒ衝突」は２．２Å、またＲ1とＲ2は二原子間距離とした。
【０１４０】
（ステップ７０：リガンドのコンフォメーションを大きく改変する。）
リガンドのコンフォメーションを大きく動かして、ステップ５０から再スタートを行い、ステップ６０までを繰り返して最適化を行う。コンフォメーション改変は、ステップ６０で算出されたリガンド座標を中心に並進、回転、シングルボンドまわりの角度をＲＳＭＤで０．３以上になるよう座標変化させる。
【０１４１】
リガンドのコンフォメーションを大きく動かした最適化は、例えばステップ５０で最適化したエネルギー「Ｕ最適」でのコンフォメーションに対して、回転可能な結合をランダムに５つ選び、原子タイプごとに決められた回転角度間隔に従ってランダムに回転させる。その後ステップ５０、ステップ６０の過程を、５０００回繰り返し行う。ただし、リガンドのコンフォメーションを変化させた後、リガンドの内部エネルギー「Ｕ内部」を計算しその値が５００．０以上のときはその後の計算をスキップし、次のリガンドコンフォメーションを発生させるようにする。
【０１４２】
（ステップ８０：標的タンパク質とリガンドとの相互作用エネルギーの決定）
ステップ４０から７０まで「Ｕ最適」が最適値となる最適なタンパク質とリガンドとの複合体座標、最適エネルギー「Ｕ最適」を算出する。
【０１４３】
（ステップ９０：ステップ４０に戻り、ステップ０中の別のリガンドを選択し、ステップ８０まで計算する。）
ステップ４０から９０まではステップ０中の化合物データベース中のリガンド全てについて行われる。
【０１４４】
（ステップ１００：リガンドの選択）
ステップ９０まで評価されたタンパク質とリガンドとの複合体座標、最適エネルギー「Ｕ最適」を基にステップ０中のデータベース中のリガンドから該当タンパク質と結合する可能性のある化合物を選ぶ。
【０１４５】
以下に、本発明にかかる相互作用関数によって該当タンパク質と結合する物質をコンピュータプログラムにより評価し、選定する方法及びその装置、該当タンパク質に結合するリガンドを選択する際に一連の処理を全自動または手動的に行うことを特徴とするコンピュータプログラム及びその装置、該当タンパク質の立体構造に際し、基準振動計算を行い、各アミノ酸のゆらぎの大きさを求め、そのゆらぎの大きさを拘束条件として、分子動力学計算を行うことで、タンパク質の立体構造をエネルギー最適構造よりおおきく離れないようにして分子動力学計算を行う方法及びその装置、タンパク質の動的挙動に関する物理化学的性質を明確に描写することを特徴とする手法及びその装置、タンパク質の動的性質を表現する関数として基準振動解析結果またはタンパク質の二次構造判定結果を用いることを特徴とした方法及びその装置、核磁気共鳴スペクトルの解析結果のような複数の立体構造座標である場合についても、該当タンパク質と結合するリガンドの探索について、複数座標すべてを全自動的にかつ短時間で同等に評価することを可能とした方法及びその装置を具体的な例を用いて詳細に説明する。
下記の実施例は、本発明の具体的な認識を得る一助と見るべきであり、本発明の範囲を何ら制限するものではない。
【実施例１】
【０１４６】
（二面角拘束ＭＤおよびクラスタリングにおけるパラメータ定数の決定）
基準振動解析により二面角のゆらぎ値が計算される。本発明においては、二面角のゆらぎ値を分子動力学計算における拘束条件として、
Ｕ二面角  ＝  Ｋ二面角  ＊  （θ−θ平衡）2
θは二面角（単位ｒａｄ）。
中の「Ｋ二面角」に適応する。実際には「Ｋ二面角」の最大値と最小値を指定することで、その範囲内で主鎖二面角揺らぎに対応するように各二面角に対して不均一な拘束がかかるようにしている。本実施例では、「Ｋ二面角」の適切な最大値と最小値を決定することを目的とする。
【０１４７】
さらに分子動力学計算後、構造変化した座標をクラスター解析し代表構造を選択する。その際、あらかじめ指定した活性部位に対して、ＭＤの途中経過１００ｆｓｅｃごとの受容体を重ね合わせた構造及び初期構造の活性部位を母集団とし、まず初めに、クラスタリングすることにより側鎖の動的情報が失われる可能性が高いことから、側鎖の二面角χにおいて母集団のα％が平均角度±２０．０°の範囲で保存されている全側鎖二面角を収集した。ただし、主鎖の根元に近い方からχが保存されていないと判定された場合はそれ以降のχは保存されていないものとした。次に、収集した保存側鎖二面角をすべて網羅している構造を母集団から抽出した。そして、抽出した構造の類似性を比較するために全原子ｒｍｓ（ｒｏｏｔ  ｍｅａｎ  ｓｑｕａｒｅ）がβÅ以下の場合、同一構造と判断して一方を削除、最終的に選ばれた構造をもとに受容体動的構造クラスターを作成した。また、保存されていなかった二面角χを構成する原子では、変動する可能性が高いことから活性部位・リガンド結合計算において衝突しても良いことにした。ここでα、βは定数となるが、本実施例では、適切なα、βを決定することを目的とする。
【０１４８】
ここでは、リガンドと接触している活性部位において最も良い主鎖の動的構造を得ることが目的であるので、ｒｍｓ（ｒｏｏｔ  ｍｅａｎ  ｓｑｕａｒｅ）を計算するときは活性部位における主鎖原子（Ｎ、Ｃα、Ｃ、Ｏ）の４原子のみが対象とした。
【０１４９】
「Ｋ二面角」の最大値と最小値、クラスタリング定数α、βはＮＭＲで解析される構造を再現できる値が適切であると考えられる。ジヒドロ葉酸還元酵素（ＤＨＦＲ、ＰＤＢ  ｃｏｄｅ：１ＬＵＤ）はＮＭＲで解析された構造である。そこでまず、ＮＭＲ構造のジヒドロ葉酸還元酵素（ＤＨＦＲ、ＰＤＢ  ｃｏｄｅ：１ＬＵＤ）のＭＯＤＥＬ１を初期構造に基準振動解析を行い、ゆらぎ値を求めその後分子動力学計算をおこなった。また、１ＬＵＤ（ＭＯＤＥＬ１）に含まれていたリガンドの各原子から半径６Å以内に含まれる受容体残基を活性部位と定義した。分子動力学計算後は、受容体動的構造クラスタリングまで行った。ただし、ＭＤは０〜０．１ｎｓｅｃまでの結果を使用した。ここで、ＭＤでの拘束の最小値と最大値は０から１０００まで（１００ごと）、クラスタリングにおいては、定数αは０％から９０％まで（１０％ごと）、定数βはＮＭＲ構造平均値ｒｍｓを参考に、０．１Åから０．６Åまで（０．１Åごと）の数値を網羅的に行い、１ＬＵＤにおけるＮＭＲ構造すべてと比較することにより定数を決定することにした。
【０１５０】
ＮＭＲ構造平均値は、受容体動的構造クラスタリングにおける定数βを決定するための参考として求めた。ｔｈｅ  Ｐｒｏｔｅｉｎ  Ｄａｔａ  Ｂａｎｋ（ＰＤＢ）のＮＭＲ構造のうち、受容体が単純タンパク質で、１つのＰＤＢファイル内に記載されていたＮＭＲ構造が１０パターン以上あり、リガンドを含む１１７種類を対象に、活性部位のＮＭＲ構造平均値ｒｍｓを求めることにした。
【０１５１】
まず、ＭＯＤＥＬ１においてリガンドの各原子から半径６Å以内に含まれる受容体残基を活性部位と定義した。ＭＯＤＥＬ１以外の構造において、ＭＯＤＥＬ１の活性部位とのｒｍｓをそれぞれ求め、さらにその平均ｒｍｓを求めた。ここで、平均ｒｍｓが１．０Å以上の場合は明らかな動的構造と見なせるので、そのようなＰＤＢファイルを対象からはずした。これにより、対象となるＰＤＢファイルは７１種類となった。７１種類の平均ｒｍｓをさらに平均化した値をＮＭＲ構造平均値ｒｍｓとした。このようにして得られたＮＭＲ構造平均値ｒｍｓは０．６２となった。
【０１５２】
「Ｋ二面角」の適切な最大値と最小値、クラスタリングにおける定数α、βの決定に関しては、各パラメータ値とＮＭＲ構造との比較を行った。
【０１５３】
１ＬＵＤには２４種類のＭＯＤＥＬが含まれており、ＭＯＤＥＬ１を対象にしたので、これを除く２３種類のＭＯＤＥＬの活性部位を正解構造とした。計算の結果出力された各受容体動的構造クラスターにおいて各正解構造とｒｍｓを計算しその中で最小のｒｍｓを「ＲＭＳ最小」として、各受容体動的構造クラスターから得られた「ＲＭＳ最小」の平均値をスコアとし、このスコアが最小となるパラメータを採用することにした。
【０１５４】
図１１に１ＬＵＤのＭＯＤＥＬ１における基準振動解析結果を、図１２、図１３、図１５〜１８にスコアとパラメータとの比較の結果を示す。
【０１５５】
図１１では、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示した。紫色Ｃは活性部位である。
【０１５６】
図１２において、この結果、定数αは７０％が良いという結果になったが、一般性を持たせときに７０％ではクラスタリングの精度が低下した場合もあったので、これよりマイルドにした８０％を定数αの値にした。
【０１５７】
図１３および図１５〜１８において、クラスタリング定数をα＝８０．０％、β＝０．４Åに固定した。なお、黒色に近いほどスコアが小さい。
【０１５８】
これらの結果より、スコアが小さくなる拘束条件としては図１４の値が最適であると判断された。
【０１５９】
これらの値の妥当性は、例えば、主鎖原子のみではなく、Ｃα原子、側鎖原子、全原子について調査しても図１４のパラメータ値が最適であることが分かる。
【実施例２】
【０１６０】
（拘束パラメータ有無による分子動力学計算の相違）
本発明によって開発された拘束パラメータを適応した分子動力学計算を２．０ｎｓｅｃまで行い、活性部位の主鎖原子の動的挙動が拘束パラメータを適応しない場合と比較して構造がどの程度変化するかを調べた。
【０１６１】
Ｃａｓｅ１）
ジヒドロ葉酸還元酵素（１ＬＵＤのＭＯＤＥＬ１）を対象に検証した。結果を図１９〜図２１に示す。基準振動計算結果は。実施例１で求めた値を適応した。
【０１６２】
図１９では、１ＬＵＤのＰＤＢファイル内に記載されている２４種類の各モデル構造をＭＯＤＥＬ１と活性部位の主鎖原子においてｒｍｓを計算し、その平均ｒｍｓを点線で表示。二面角拘束があるとき（Ａ）とないとき（Ｂ）において、活性部位の主鎖原子の初期構造からのずれをｒｍｓで表示。
【０１６３】
図２０に二面角拘束なしでＭＤを計算させた時のＮＭＲ構造との比較を示す。図２０において、白色はＮＭＲ構造（１ｌｕｄ）であり、黒色はＭＤ構造（１ｌｕｄ）である。
【０１６４】
表１に二面角拘束なしでＭＤを計算させた時のＮＭＲ構造との比較を示す。
【表１】
2005242493.tif 000025
【０１６５】
図２１に二面角拘束ありでＭＤを計算させた時のＮＭＲ構造との比較を示す。
【０１６６】
表２に二面角拘束ありでＭＤを計算させた時のＮＭＲ構造との比較を示す。
【表２】
2005242493.tif 000026
【０１６７】
Ｃａｓｅ２）
ここでは、ＦＭＡＳ［Ｏｇａｔａ  Ｋ．，  Ｕｍｅｙａｍａ  Ｈ．  （２０００）  Ａｎ  ａｕｔｏｍａｔｉｃ  ｈｏｍｏｌｏｇｙ  ｍｏｄｅｌｉｎｇ  ｍｅｔｈｏｄ  ｃｏｎｓｉｓｔｉｎｇ  ｏｆ  ｄａｔａｂａｓｅ  ｓｅａｒｃｈｅｓ  ａｎｄ  ｓｉｍｕｌａｔｅｄ  ａｎｎｅａｌｉｎｇ  Ｊ．  Ｍｏｌ．  Ｇｒａｐｈｉｃｓ  Ｍｏｄ．  １８，  ２５８−２７２］によりモデリングした構造（モデル構造）とＸ線構造を初期構造に選び、初期構造及び拘束の有無に依存することを検証した。また、リガンドの各原子から半径１０Å以内に含まれる受容体残基を活性部位と定義した。
【０１６８】
ｃｅｌｌｕｌａｒ  ｒｅｔｉｎｏｉｃ  ａｃｉｄ  ｂｉｎｄｉｎｇ  ｐｒｏｔｅｉｎ  ｔｙｐｅ  ＩＩ（ＣＲＡＢＰ−ＩＩ）（ＰＤＢ  ｃｏｄｅ：１ＣＢＱ）のＸ線構造（立体構造）を利用した。また、参照タンパク質にホモロジー３２．１％のｉｎｔｅｓｔｉｎａｌ  ｆａｔｔｙ  ａｃｉｄ  ｂｉｎｄｉｎｇ  ｐｒｏｔｅｉｎ（ＰＤＢ  ｃｏｄｅ：１ＩＣＭ  ）を選び、図２２のアライメントでモデル構造を作成した。図２３、図２４、図２５にＸ線とモデルの構造比較を示す。
【０１６９】
図２３には、１ＣＢＱの立体構造（Ｘ線構造（赤色Ａ）およびモデル構造（青色Ｂ））を示す。図２４には、図２３の緑色Ｃで示される物質である６−（２，３，４，５，６，７−ｈｅｘａｈｙｄｒｏ−２，４，４−ｔｒｉｍｅｔｈｙｌ−１−ｍｅｔｈｙｌｅｎｅｉｎｄｅｎ−２−ｙｌ）−３−ｍｅｔｈｙｌｈｅｘａ−２，４−ｄｉｅｎｏｉｃ  ａｃｉｄの構造を示す。
【０１７０】
図２５には、１ＣＢＱのＸ線構造とモデル構造の相違をｒｍｓで表示した。
【０１７１】
図２６は１ＣＢＱのＸ線構造の基準振動解析の結果を示す図であり、図２７は１ＣＢＱのモデル構造の基準振動解析の結果を示す図である。図２６および図２７では、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。
【０１７２】
図２８には、１ＣＢＱのＸ線構造とモデル構造の分子動力学法（ＭＤ）計算の結果を示す。Ｘ線構造の活性部位の主鎖原子とのｒｍｓを求めた。図２８において、Ａは初期構造がＸ線構造で二面角拘束なし、Ｂは初期構造がＸ線構造で二面角拘束あり、Ｃは初期構造がモデル構造で二面角拘束なし、Ｄは初期構造がモデル構造で二面角拘束あり、である。
【０１７３】
Ｃａｓｅ３）
ＦｌａｖｏｄｏｘｉｎのＸ線構造（ＰＤＢ  ｃｏｄｅ：１Ｊ９Ｇ）を利用した。また、参照タンパク質にホモロジー２９．２％のｆｌａｖｏｄｏｘｉｎ（ＰＤＢ  ｃｏｄｅ：１ＡＨＮ）を選び、図２９のアライメントでモデル構造を作成した。図２９には、１Ｊ９Ｇおよび１ＡＨＮのアライメントを示す。
【０１７４】
図３０には、１Ｊ９Ｇの立体構造（Ｘ線構造（赤色Ａ）およびモデル構造（青色Ｂ））を示す。図３１には、図３０における緑色Ｃで示される物質であるｆｌａｖｉｎ  ｍｏｎｏｎｕｃｌｅｏｔｉｄｅの構造を示す。
【０１７５】
図３２には、１Ｊ９ＧのＸ線構造とモデル構造の相違をｒｍｓで表示する。
【０１７６】
図３３には１Ｊ９ＧのＸ線構造の基準振動解析の結果を、図３４には１Ｊ９Ｇのモデル構造の基準振動解析の結果を示す。図３３および図３４において、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥ［８］による二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。
【０１７７】
図３５には、１Ｊ９ＧのＸ線構造とモデル構造の分子動力学法（ＭＤ）計算の結果を示す。Ｘ線構造の活性部位の主鎖原子とのｒｍｓを求めた。Ａは初期構造がＸ線構造で二面角拘束なし。Ｂは初期構造がＸ線構造で二面角拘束あり。Ｃは初期構造がモデル構造で二面角拘束なし。Ｄは初期構造がモデル構造で二面角拘束あり。
【０１７８】
Ｃａｓｅ４）
Ｍａｔｒｉｘ  ｍｅｔａｌｌｏｐｒｏｔｅｉｎａｓｅ−８（ＭＭＰ−８）のＸ線構造（ＰＤＢ  ｃｏｄｅ：１ＭＭＢ）を利用した。また、参照タンパク質にホモロジー５５．０％のＭＭＰ−３（ＰＤＢ  ｃｏｄｅ：１Ｂ３Ｄ）を選び、図３６のアライメントでモデル構造を作成した。図３６には、１ＭＭＢおよび１Ｂ３Ｄ＿Ａのアライメントを示す。
【０１７９】
図３７には、１ＭＭＢの立体構造（Ｘ線構造（赤色Ａ）およびモデル構造（青色Ｂ））を示す。図３８には、図３７における緑色Ｃで示される物質であるｂａｔｉｍａｓｔａｔの構造を示す。
【０１８０】
図３９には、１ＭＭＢのＸ線構造とモデル構造の相違をｒｍｓで表示する。
【０１８１】
図４０には１ＭＭＢのＸ線構造の基準振動解析の結果を、図４１には１ＭＭＢのモデル構造の基準振動解析の結果を示す。図４０および図４１において、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥ［８］による二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。
【０１８２】
図４２には、１ＭＭＢのＸ線構造とモデル構造の分子動力学法（ＭＤ）計算の結果を示す。Ｘ線構造の活性部位の主鎖原子とのｒｍｓを求めた。Ａは初期構造がＸ線構造で二面角拘束なし。Ｂは初期構造がＸ線構造で二面角拘束あり。Ｃは初期構造がモデル構造で二面角拘束なし。Ｄは初期構造がモデル構造で二面角拘束あり。
【０１８３】
Ｃａｓｅ１）〜Ｃａｓｅ４）に示したとおり、拘束パラメータを適応した分子動力学計算結果は、拘束パラメータを適応しない場合と比較して、大きな構造変化は少ない。このことは古典力学を適応しているため大きな構造変化をしてしまう分子動力学法において、拘束パラメータを適応することで大きな構造変化を合理的に拘束することができ，理想的な構造座標を得ることが可能であるということを示している。また、ホモロジーが高ければ、ＦＭＡＳの構造構築精度もあがる。すなわち、Ｘ線に近い構造を得られるので、アミノ酸の数個異なるミューテーションタンパク質にもこの手法は利用できる。
【実施例３】
【０１８４】
（タンパク質／リガンド複合体モデルの検証）
本発明により該当タンパク質に結合するリガンドの複合体立体構造が予測される。本実施例では、こうして予測された複合体立体構造座標の予測精度を検証する。検証には、複合体の立体構造が既知で、リガンドの有無もしくはリガンドの種類により活性部位の形が異なるＩｎｄｕｃｅｄ−Ｆｉｔ型のタンパク質を用いた。ここで、リガンドの各原子から半径１０Å以内の残基をタンパク質の活性部位と定義した。また、Ｘ線構造またはＮＭＲ構造を初期構造に選んだＭＤでは、ほぼ一定の構造を保ち続けることが分かったのでＭＤを１．０ｎｓｅｃまで行うことにした。ただし、水素原子を除いて計算した。複合体モデル構築は、発明実施の形態に従って行った。
【０１８５】
Ｃａｓｅ１）
ジヒドロ葉酸還元酵素（ＤＨＦＲ）である１ＢＺＦと１ＬＵＤとはホモロジーが１００．０％でかつ結合しているリガンドが異なることにより活性部位の形が異なる。そこで、１ＢＺＦ（ＭＯＤＥＬ１８）を初期構造として選択し、リガンドとして２，４−ｄｉａｍｉｎｏ−５−（３，４，５−ｔｒｉｍｅｔｈｏｘｙ−ｂｅｎｚｙｌ）−ｐｙｒｉｍｉｄｉｎ−１−ｉｕｍ（図４９）を用い、本発明によるプログラムによってタンパク質／リガンド複合体モデルを作成し、正解構造である１ＬＵＤ（ＭＯＤＥＬ４）と比較することで検証した（図４３）。
【０１８６】
図４３に、ジヒドロ葉酸還元酵素の立体構造を示す。図４３では、１ＬＵＤ（ＭＯＤＥＬ４）受容体（緑色Ａ）とリガンド（赤色Ｂ）、１ＢＺＦ（ＭＯＤＥＬ１８）の受容体（青色Ｃ）とリガンド（水色Ｄ）を示す。
【０１８７】
図４４に１ＢＺＦの基準振動計算解析を示す。図４４では、二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。
【０１８８】
図４５、図４７に１ＢＺＦを用いた拘束二面角分子動力学の結果を示す。図４５には、正解構造１ＬＵＤ（ＭＯＤＥＬ４）の活性部位とのｒｍｓを表示。Ａは主鎖原子。Ｂは側鎖原子。Ｃは全原子。図４７は、１ＢＺＦ（ＭＯＤＥＬ１８）における活性部位・リガンド結合解析であり、ＭＤ計算を０．１ｎｓｅｃまで行ったときの結合解析及び１．０ｎｓｅｃまで行ったとき、また、受容体動的構造クラスタリングを行うときの母集団を１００ｆｓｅｃごと及び１０００ｆｓｅｃごとで行ったときの結合解析。評価法は、正解構造との活性部位及びリガンドのｒｍｓで行った。
【０１８９】
図４６にリガンドドッキングにおける空間点指定のパラメータ値を示す。図４６は、１ＬＵＤ（ＭＯＤＥＬ４）より得られた構造活性相関情報。
【０１９０】
図４８〜図５０にタンパク／リガンド複合体と正解構造との比較を示す。図４８は、０〜１．０ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団により作られたクラスターを用いて行った活性部位・リガンド結合。緑色Ａは正解構造の１ＬＵＤ（ＭＯＤＥＬ４）。青色Ｂは初期構造の１ＢＺＦ（ＭＯＤＥＬ１８）。そして赤色Ｃはリガンド結合における最適構造。要素色Ｄはリガンドの正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは０．９６１４。図５０のリガンドを結合させることにより活性部位の主鎖原子のｒｍｓで０．２７９１の誘導が生じた。図４９は黒色：正解（１ｌｕｄのｍｏｄｅｌ  ４）、灰色：初期構造（１ｂｚｆのｍｏｄｅｌ  １８）、白色：最適構造。図５０は、２，４−ｄｉａｍｉｎｏ−５−（３，４，５−ｔｒｉｍｅｔｈｏｘｙ−ｂｅｎｚｙｌ）−ｐｙｒｉｍｉｄｉｎ−１−ｉｕｍ。１ＬＵＤのリガンド。
【０１９１】
Ｃａｓｅ２）
ｈｅａｔ  ｓｈｏｃｋ  ｐｒｏｔｅｉｎ  ９０（ＨＳＰ９０）である１ＹＥＲと１ＹＥＴはホモロジーが１００．０％でリガンド結合の有無により活性部位の形が異なる。そこで、リガンド結合していない１ＹＥＲを初期構造に選び、リガンドとしてｇｅｌｄａｎａｍｙｃｉｎを用い、正解構造である１ＹＥＴと比較することで検証した（図５１）。図５１は、ｈｅａｔ  ｓｈｏｃｋ  ｐｒｏｔｅｉｎ  ９０の立体構造である。１ＹＥＴの受容体（緑色Ａ）とリガンド（赤色Ｂ）。１ＹＥＲの受容体（青色Ｃ）。
【０１９２】
図５２に１ＹＥＲの基準振動解析を示す。二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。
【０１９３】
図５３、図５５に１ＹＥＲを用いた拘束二面角分子動力学の結果を示す。正解構造１ＹＥＴの活性部位とのｒｍｓを表示。Ａは主鎖原子。Ｂは側鎖原子。Ｃは全原子。図５５は、１ＹＥＲにおける活性部位・リガンド結合解析。ＭＤ計算を０．１ｎｓｅｃまで行ったときの結合解析及び１．０ｎｓｅｃまで行ったとき、また、受容体動的構造クラスタリングを行うときの母集団を１００ｆｓｅｃごと及び１０００ｆｓｅｃごとで行ったときの結合解析。評価法は、正解構造との活性部位及びリガンドのｒｍｓで行った。
【０１９４】
図５４にリガンドドッキングにおける空間点指定のパラメータ値を示す。図５４は、１ＹＥＴより得られた構造活性相関情報。
【０１９５】
図５６および図５７にタンパク／リガンド複合体と正解構造との比較を示す。図５６および図５７は１ＹＥＲ・リガンド結合。図５６は、０〜０．１ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団により作られたクラスターを用いて行った活性部位・リガンド結合。緑色Ａは正解構造の１ＹＥＴ。青色Ｂは初期構造の１ＹＥＲ、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄはリガンドの正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは１．２０８１。図５７のリガンドを結合させることにより活性部位の主鎖原子のｒｍｓで０．１６１９の誘導が生じた。図５７は、ｇｅｌｄａｎａｍｙｃｉｎ。１ＹＥＴのリガンド。
【０１９６】
Ｃａｓｅ３）
ｍｉｔｏｇｅｎ−ａｃｔｉｖａｔｅｄ  ｐｒｏｔｅｉｎ  ｋｉｎａｓｅ（ＭＡＰ  ｋｉｎａｓｅ）である１Ａ９Ｕと１ＯＵＫはホモロジー１００．０％でかつ結合しているリガンドが異なることにより活性部位の形が異なる。そこで、１Ａ９Ｕを初期構造に選び、リガンドとして１ＯＵＫ中に含まれるリガンドを用い、正解構造である１ＯＵＫと比較することでと検証した（図５８）。図５８はｍｉｔｏｇｅｎ−ａｃｔｉｖａｔｅｄ  ｐｒｏｔｅｉｎ  ｋｉｎａｓｅの立体構造。図５８において、１ＯＵＫの受容体（緑色Ａ）とリガンド（赤色Ｂ）。１Ａ９Ｕの受容体（青色Ｃ）とリガンド（水色Ｄ）。
【０１９７】
図５９に基準振動解析結果を示す。図５９は１Ａ９Ｕの基準振動解析。二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。
【０１９８】
図６０、図６２に１ＹＥＲを用いた拘束二面角分子動力学の結果を示す。図６０は１Ａ９ＵのＭＤ。図６０に正解構造１ＯＵＫの活性部位とのｒｍｓを表示。Ａは主鎖原子。Ｂは側鎖原子。Ｃは全原子。図６２は１Ａ９Ｕにおける活性部位・リガンド結合解析。ＭＤ計算を０．１ｎｓｅｃまで行ったときの結合解析及び１．０ｎｓｅｃまで行ったとき、また、受容体動的構造クラスタリングを行うときの母集団を１００ｆｓｅｃごと及び１０００ｆｓｅｃごとで行ったときの結合解析。評価法は、正解構造との活性部位及びリガンドのｒｍｓで行った。
【０１９９】
図６１にリガンドドッキングにおける空間点指定のパラメータ値を示す。図６１は１ＯＵＫより得られた構造活性相関情報。
【０２００】
図６３〜図６５にタンパク／リガンド複合体と正解構造との比較を示す。図６３〜図６５は１Ａ９Ｕ・リガンド結合。図６３は、０〜０．１ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団により作られたクラスターを用いて行った活性部位・リガンド結合。緑色Ａは正解構造の１ＯＵＫ。青色Ｂは初期構造の１Ａ９Ｕ、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄはリガンドの正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは１．６１１２。図６５のリガンドを結合させることにより活性部位の主鎖原子のｒｍｓで０．１８７１の誘導が生じた。図６４において、黒色：正解（１ｏｕｋ）、灰色：初期構造（１ａ９ｕ）、白色：最適構造。図６５は１ＯＵＫのリガンド。４−［５−［２−（１−ｐｈｅｎｙｌ−ｅｔｈｙｌａｍｉｎｏ）−ｐｙｒｉｍｉｄｉｎ−４−ｙｌ］−１−ｍｅｔｈｙｌ−４−（３−ｔｒｉｆｌｕｏｒｏｍｅｔｈｙｌｐｈｅｎｙｌ）−１Ｈ−ｉｍｉｄａｚｏｌ−２−ｙｌ］−ｐｉｐｅｒｉｄｉｎｅ。
【０２０１】
Ｃａｓｅ１）〜Ｃａｓｅ３）に示す通り、本発明により作成されるタンパク質／リガンド複合体モデルは、誘導結合型のタンパク質／リガンド複合体の立体構造を精度よく予測可能であることが分かる。
【実施例４】
【０２０２】
（Ｆｘａを用いたｉｎ  ｓｉｌｉｃｏ  Ｓｃｒｅｅｎｉｎｇへの応用例）
本発明により、セリンプロテアーゼの１種であるＦｘａの立体構造（図６６）を用い、化合物データベースからＦｘａに結合する可能性のあるリガンドを探索した。立体構造には、１ＡＩＸを用い、リガンドデータベースとして、ＰＤＢデータベースより収集した３６３３種類のリガンドを用いた。発明実施の形態に従い、ｉｎ  ｓｉｌｉｃｏ  ｓｃｒｅｅｎｉｎｇを行った。その結果を図６７に示す。
【０２０３】
図６７は、化合物データベース中のリガンドのうち、１ＡＩＸとの相互作用エネルギーの上位１００個を示している。図６７で、太字は１ＡＩＸ中に含まれているリガンドで、斜線はセリンプロテアーゼ。ＰＤＢ  ｃｏｄｅとは、リガンドが含まれているもとのＰＤＢｃｏｄｅを示す。図６７には、１ＡＩＸにもともと含まれているリガンドがランキング１９位に入っている。
【０２０４】
ランキング１９位におけるタンパク質／リガンド複合体構造とを図６８および図６９に示す。図６８において、白色は受容体、黒色は１ＡＩＸのリガンド。図６９は１ＡＩＸ中のリガンド。
【０２０５】
図６７中のランキング３５位、３８位、８０位はすべてセリンプロテアーゼに結合するリガンドである。
【０２０６】
これらの構造とおけるタンパク質／リガンド複合体構造を図７０および図７１、図７２および図７３、図７４および図７５に示す。図７０および図７１はランキング３５位におけるタンパク質／リガンド複合体構造。図７０において、白色は受容体、黒色は１ＡＵＪのリガンド。図７１は１ＡＵＪ中のリガンド。図７２および図７３はランキング３８位におけるタンパク質／リガンド複合体構造。図７２において、白色は受容体、黒色は正解（１ＦＯＲ）のリガンド、ＲＭＳは１．５００。図７３は１ＦＯＲ中のリガンド。図７４および図７５はランキング８０位におけるタンパク質／リガンド複合体構造。図７４において、白色は受容体、黒色は１Ｋ１Ｍのリガンド。図７５は１Ｋ１Ｍ中のリガンド。
【０２０７】
これらの結果から、本発明により、化合物データベースからもっともらしい化合物を選択することが可能であることが分かる。
【実施例５】
【０２０８】
（異なる条件でのｉｎ  ｓｉｌｉｃｏスクリーニング）
構造活性相関（ＳＡＲ）の情報により順位が変動することを検証する。また、受容体を固定した場合の順位の変動も検証する。
【０２０９】
ここでは、ｓｅｖｅｒｅ  ａｃｕｔｅ  ｒｅｓｐｉｒａｔｏｒｙ  ｓｙｎｄｒｏｍｅ（ＳＡＲＳ）のプロテアーゼを用いたｉｎ  ｓｉｌｉｃｏ  ｓｃｒｅｅｎｉｎｇを行った。初期構造にはリガンドを含まない１ＵＫ３（Ｂ鎖）を、またリガンドを含む１ＵＫ４（Ｂ鎖）のリガンド結合様式を構造活性相関情報として利用した。活性部位は１ＵＫ４（Ｂ鎖）のリガンドの各原子から半径１０Å以内に含まれる受容体残基部位。リガンドデータベースとして、ＰＤＢより収集した３６３３種類のリガンドを用いた。ただし、結合解析で利用する受容体動的構造クラスターには０〜０．１ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団で作られたものを使用した。また、水素原子を除いて計算した。
【０２１０】
図７６はＳＡＲＳプロテアーゼの立体構造。１ＵＫ４（Ｂ鎖）の受容体（緑色Ａ）とリガンド（赤色Ｂ）。１ＵＫ３（Ｂ鎖）の受容体（青色Ｃ）。
【０２１１】
図７７に１ＵＫ３の基準振動解析の結果を示す。図７７は１ＵＫ３（Ｂ鎖）の基準振動解析。二面角φ（橙色Ａ）、ψ（緑色Ｂ）の揺らぎの大きさを示した。揺らぎが０．０に近いほど分子動力学法（ＭＤ）計算において二面角拘束が強くなる。また、ＳＴＲＩＤＥによる二次構造判定でα−へリックス（赤色Ｄ）、β−シート（青色Ｄ）を表示。紫色Ｃは活性部位。
【０２１２】
図７８に１ＵＫ３の分子動力学計算の結果を示す。図７８は１ＵＫ３（Ｂ鎖）のＭＤ。１ＵＫ４（Ｂ鎖）の活性部位とのｒｍｓを表示。Ａは主鎖原子。Ｂは側鎖原子。Ｃは全原子。
【０２１３】
Ｃａｓｅ１）  ＳＡＲ４ヵ所指定
図７９に活性部位内での空間指定を示す。図７９は１ＵＫ４（Ｂ鎖）より得られた構造活性相関情報。
【０２１４】
図８０にｉｎ  ｓｉｌｉｃｏ  スクリーニングの結果を示す。図８０は１ＵＫ３（Ｂ鎖）におけるｉｎ  ｓｉｌｉｃｏスクリーニングの結果。
【０２１５】
図８１に正解構造との比較を示す。図８１は１ＵＫ３と１ＵＫ４との比較。順位２５。緑色Ａは１ＵＫ４（Ｂ鎖）。青色Ｂは初期構造の１ＵＫ３（Ｂ鎖）、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄは１ＵＫ４のペプチド性リガンド（ＡＳＮ−ＳＥＲ−ＴＨＲ−ＬＥＵ−ＧＬＮ）の正解構造。Ｅは計算結果によるリガンド。リガンドのｒｍｓは２．５７２１。正解構造との活性部位の主鎖原子のｒｍｓは、初期構造では１．０２４８、最適構造では１．０７９２。
【０２１６】
図８２および図８３に、ｉｎ  ｓｉｌｉｃｏ  スクリーニングの順位１を示す。図８３は１ＱＦ４のリガンド（Ｃ８−Ｒ）−ｈｙｄａｎｔｏｃｉｄｉｎ  ５’−ｐｈｏｓｐｈａｔｅ。
【０２１７】
Ｃａｓｅ２）  ＳＡＲ３ヵ所指定
図８４に１ＵＫ３の活性部位内での空間指定を示す。図８４は１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。
【０２１８】
図８５には、１ＵＫ３での最適構造と正解構造との比較を示す。図８５は１ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。順位４９。緑色Ａは１ＵＫ４（Ｂ鎖）。青色Ｂは初期構造の１ＵＫ３（Ｂ鎖）、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄは１ＵＫ４のペプチド性リガンド（ＡＳＮ−ＳＥＲ−ＴＨＲ−ＬＥＵ−ＧＬＮ）の正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは２．００５７。正解構造との活性部位の主鎖原子のｒｍｓは、初期構造では１．０２４８、最適構造では１．０４６９。
【０２１９】
図８６には、ＳＡＲ３ヵ所指定でのｉｎ  ｓｉｌｉｃｏスクリーニング結果を示す。図８６はＳＡＲ３ヵ所指定で実行したｉｎ  ｓｉｌｉｃｏスクリーニングの結果。
【０２２０】
Ｃａｓｅ３）  ＳＡＲ５ヵ所指定
図８７に１ＵＫ３の活性部位内での空間指定を示す。図８７は１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。
【０２２１】
図８８にはＳＡＲ５ヵ所指定でのｉｎ  ｓｉｌｉｃｏスクリーニング結果を示す。図８８はＳＡＲ５ヵ所指定で実行したハイスループットスクリーニングの結果。
【０２２２】
図８９には、１ＵＫ３での最適構造と正解構造との比較を示す。図８９は１ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。順位２。緑色Ａは１ＵＫ４（Ｂ鎖）。青色Ｂは初期構造の１ＵＫ３（Ｂ鎖）、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄは１ＵＫ４のペプチド性リガンド（ＡＳＮ−ＳＥＲ−ＴＨＲ−ＬＥＵ−ＧＬＮ）の正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは１．２５７８。正解構造との活性部位の主鎖原子のｒｍｓは、初期構造では１．０２４８、最適構造では１．１６２０。
【０２２３】
Ｃａｓｅ４）  リガンド原子タイプ指定の変更
図９０に１ＵＫ３の活性部位内での空間指定を示す。図９０は１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。
【０２２４】
図９１には、リガンド原子タイプ指定変更でのｉｎ  ｓｉｌｉｃｏスクリーニング結果を示す。図９１はリガンド原子タイプ指定変更で実行したハイスループットスクリーニングの結果。
【０２２５】
図９２には、１ＵＫ３での最適構造と正解構造との比較を示す。図９２は１ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。順位７７４。緑色Ａは１ＵＫ４（Ｂ鎖）。青色Ｂは初期構造の１ＵＫ３（Ｂ鎖）、そして赤色Ｃはリガンド結合における最適構造。要素色Ｄは１ＵＫ４のペプチド性リガンド（ＡＳＮ−ＳＥＲ−ＴＨＲ−ＬＥＵ−ＧＬＮ）の正解構造。水色Ｅは計算結果によるリガンド。リガンドのｒｍｓは２．５２１６。正解構造との活性部位の主鎖原子のｒｍｓは、初期構造では１．０２４８、最適構造では１．０７９２。
【０２２６】
Ｃａｓｅ５）  受容体固定
図９３に活性部位内での空間指定を示す。図９３は１ＵＫ４（Ｂ鎖）より得られた構造活性相関情報。
【０２２７】
図９４に、受容体を固定したｉｎ  ｓｉｌｉｃｏ  スクリーニングの結果を示す。図９４は受容体を固定した状態で実行したハイスループットスクリーニングの結果。
【０２２８】
図９５には、１ＵＫ３と１ＵＫ４を重ね合わせたリガンドと計算結果のリガンドとの比較を示す。順位３９。灰色は１ＵＫ３の活性部位構造。黒色は１ＵＫ３と１ＵＫ４を重ね合わせたリガンド。白色は計算結果のリガンド。
【０２２９】
Ｃａｓｅ１）〜Ｃａｓｅ４）を見ると、ＳＡＲの指定が多いほど参考にしたリガンドの順位が良くなる。つまり、参考にできるリガンドの結合情報が信頼できる場合には、ＳＡＲの情報を多くしたｉｎ  ｓｉｌｉｃｏスクリーニングを行い、信頼性に欠ける場合には、ＳＡＲの情報数を減らし、さらに、リガンド原子タイプ指定の幅を広げることで、様々なリガンドがランキング上位に分布する。そして、その分布情報をもとにＳＡＲ情報を作り変えてｉｎ  ｓｉｌｉｃｏスクリーニングを実行するとより信頼性の持てる結果が出力されるはずである。
【０２３０】
Ｃａｓｅ１）とＣａｓｅ５）を見ると、受容体の動的構造の有無による順位変動を示している。これは、リガンドの動きのみの最適化に比べ、リガンド及び受容体それぞれが動く最適化の方が原子のぶつかりをさけることに優れている。従って、同じ位置に配置するための最適化エネルギーに差が生じる。
【実施例６】
【０２３１】
（二面角拘束分子動力学計算のパラメータに関するＭＤパラメータの分布）
ＦＭＮ−ｂｉｎｄｉｎｇ  ｐｒｏｔｅｉｎにおける二面角拘束ＭＤパラメータの分布。
ここでは、１ＬＵＤ以外のＮＭＲ構造でも二面角拘束分子動力学計算及びクラスタリングのパラメータが同様の結果を生じるのかを検証する。そこで、ＦＭＮ−ｂｉｎｄｉｎｇ  ｐｒｏｔｅｉｎのＮＭＲ構造（ＰＤＢ  ｃｏｄｅ：１ＡＸＪ）のＭＯＤＥＬ１を初期構造に選んだ。評価法は、受容体動的構造クラスタリングのパラメーター（α＝８０．０％、β＝０．４Å）を固定したこと以外は実施例１に従った。
【０２３２】
図９６に１ＡＸＪにおける二面角拘束分子動力学計算のパラメータ決定のスコアの分布状況を示す。図９６は１ＡＸＪにおける二面角拘束ＭＤパラメータの分布。Ａに近い部分ほどスコアで小さい。１ＬＵＤの時と同様に二面角拘束の最大値８００、最小値０では良い結果を示す。
【実施例７】
【０２３３】
（二面角拘束ＭＤ）
ここでは、主鎖二面角拘束ＭＤで各原子の動的構造を検証する。また、時として基準振動解析が収束せず二面角揺らぎ情報が得られないことがある。そこで、図１３により主鎖二面角に対して均一な拘束（５００）でＭＤを行ったときでも良い結果になっていることから、この場合における動的構造も検証する。実施例１に従い、拘束なし、二面角揺らぎを用いた拘束及び均一な拘束（５００）の条件のもとＭＤを行った。
【０２３４】
図９７〜図１０８に１ＬＵＤに対して行った分子動力学計算の各原子における動的挙動の結果を示す。図９７〜図１０８は１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。図９７および図９８は活性部位の主鎖原子、図９９および図１００は受容体の主鎖原子、図１０１および図１０２は活性部位の側鎖原子、図１０３および図１０４は受容体の側鎖原子、図１０５および図１０６は活性部位の全原子、図１０７および図１０８は受容体の全原子、における動的挙動の結果を示す図である。１ＬＵＤのＰＤＢファイル内に記載されている２４種類の各モデル構造をＭＯＤＥＬ１と活性部位の主鎖原子においてｒｍｓを計算し、その平均ｒｍｓを点線で表示。二面角拘束があるとき（Ａ）とないとき（Ｂ）及び二面角拘束が５００で一定（Ｃ）において、活性部位の主鎖原子の初期構造からのずれをｒｍｓで表示。
【０２３５】
ここには、記載しないが１ＣＢＱ、１Ｊ９Ｇ、１ＭＭＢ、１ＢＺＦ（ＭＯＤＥＬ１８）、１ＹＥＲ、１Ａ９Ｕ及び１ＵＫ３（Ｂ鎖）に関しても主鎖二面角揺らぎに基づく拘束ＭＤの結果を見ると、図１０９〜図１１１と同様に主鎖原子の抑制があると、拘束のない側鎖原子にも一定動きを示す。受容体の動きにおいて主鎖原子の動きの比重が大きいことが理解できる。
【実施例８】
【０２３６】
（異なる条件での結合解析）
ここでは、二面角拘束ＭＤ及びクラスタリングのパラメータが異なっても誘導が生じることを検証する。拘束の最大値１００、最小値０及び受容体動的構造クラスタリング定数α＝８０．０％、β＝１．０Åに設定して、その他は実施例２に従った。ただし、受容体動的構造クラスターには０〜０．１ｎｓｅｃの範囲内で１００ｆｓｅｃごとの母集団で作られたものを使用した。また、活性部位の定義は、リガンドの各原子から半径６Å以内にある受容体残基とした。
【０２３７】
図１０９〜図１１１には、異なる条件で受容体／リガンド結合の結果を示す。
（ｉ）  １ＢＺＦ（ＭＯＤＥＬ１８）で、リガンド結合により活性部位の主鎖原子ｒｍｓで０．２６８６の誘導が生じた（図１０９）。活性部位全体のｒｍｓでは０．１２２４の誘導。リガンドのｒｍｓは０．８５２６。
（ｉｉ）  １ＹＥＲで、リガンド結合により活性部位の主鎖原子ｒｍｓで０．２３７６の誘導が生じた（図１１０）。活性部位全体のｒｍｓでは０．０８１６の誘導。リガンドのｒｍｓは０．７２４６。
（ｉｉｉ）  １Ａ９Ｕで、リガンド結合により活性部位の主鎖原子ｒｍｓで０．２１５０の誘導が生じた（図１１１）。活性部位全体のｒｍｓでは０．０４６４の誘導。リガンドのｒｍｓは０．９４６４。
【０２３８】
ただし、緑色は正解構造、青色は初期構造、赤色は最適構造。要素色は正解リガンド、水色は最適リガンド。
【０２３９】
図１０９〜図１１１で示すように、各条件が異なっていても、与えられた条件の中で最適な結果を生じることができる。
【実施例９】
【０２４０】
（正解構造を初期構造に選んだときの結合解析）
ここでは、ＤＨＦＲの１ＢＺＦ及び１ＬＵＤはリガンドの結合様式が似ているので、構造活性相関情報を一部変更し１ＢＺＦのリガンドの結合解析を行う。条件としては、初期構造に１ＢＺＦ（ＭＯＤＥＬ１８）、０〜０．１ｎｓｅｃまでの母集団より作成したクラスターを使用した。
【０２４１】
図１１２には、１ＢＺＦの活性部位内における空間指定を示す。図１１２は１ＢＺＦ用に変更した構造活性相関情報。
【０２４２】
図１１３および図１１４には、１ＢＺＦにおける受容体／リガンド結合の結果を示す。図１１３および図１１４は１ＢＺＦ（ＭＯＤＥＬ１８）のリガンド結合解析。
（ｉ）  最適化したときの受容体には初期構造を選択し、リガンドのｒｍｓは０．８８８４。（図１１３）
（ｉｉ）  ｔｒｉｍｅｔｒｅｘａｔｅ。１ＢＺＦ（ＭＯＤＥＬ１８）のリガンド。（図１１４）
【０２４３】
初期構造が元々ＰＤＢに登録されていた構造、つまり最適構造であったため、計算結果でもそれが図１１３および図１１４のように再現できる。
【産業上の利用可能性】
【０２４４】
以上、本発明の方法は、医農薬の分子設計等を中心に、受容体/リガンド結合の解析を行う分野（医薬品設計）において、極めて有用であると考えられる。本発明は、産業上多くの分野、特に医薬品、食品、化粧品、医療、構造解析、機能解析等の分野で広く実施することができ、故に極めて有用である。
【図面の簡単な説明】
【０２４５】
【図１】本発明によるタンパク質立体構造と誘導適合を利用したリガンド探索方法の一例を示すフローチャートである。
【図２】ｓｐ2軌道原子におけるダミー水素原子発生。
【図３】金属原子におけるダミー原子発生。
【図４】構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。
【図５】構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。
【図６】構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。
【図７】構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。
【図８】構造活性相関（ＳＡＲ）情報をもとに活性部位内にリガンドを入れるための初期座標（Ｂ）発生。
【図９】水素結合角の定義。
【図１０】スタッキングにおける角度の定義。
【図１１】１ＬＵＤ（ＭＯＤＥＬ１）の基準振動解析の結果。
【図１２】ＭＤ及びクラスタリングのパラメータとスコア。
【図１３】クラスタリング定数を固定した時のＭＤでの二面角拘束の最大値、最小値の分布。
【図１４】拘束パラメータ。
【図１５】クラスタリングのパラメータを固定した時の二面角拘束パラメータの分布。
【図１６】クラスタリングのパラメータを固定した時の二面角拘束パラメータの分布。
【図１７】クラスタリングのパラメータを固定した時の二面角拘束パラメータの分布。
【図１８】クラスタリングのパラメータを固定した時の二面角拘束パラメータの分布。
【図１９】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤの結果。
【図２０】二面角拘束なしでＭＤを計算させた時のＮＭＲ構造との比較。
【図２１】二面角拘束ありでＭＤを計算させた時のＮＭＲ構造との比較。
【図２２】１ＣＢＱのアライメント。
【図２３】１ＣＢＱの立体構造。
【図２４】１ＣＢＱの立体構造。
【図２５】１ＣＢＱのＸ線構造とモデル構造の相違をｒｍｓで表示。
【図２６】１ＣＢＱのＸ線構造とモデル構造の基準振動解析。
【図２７】１ＣＢＱのＸ線構造とモデル構造の基準振動解析。
【図２８】１ＣＢＱのＸ線構造とモデル構造のＭＤ。
【図２９】１Ｊ９Ｇのアライメント。
【図３０】１Ｊ９Ｇの立体構造。
【図３１】１Ｊ９Ｇの立体構造。
【図３２】１Ｊ９ＧのＸ線構造とモデル構造の相違をｒｍｓで表示。
【図３３】１Ｊ９ＧのＸ線構造とモデル構造の基準振動解析。
【図３４】１Ｊ９ＧのＸ線構造とモデル構造の基準振動解析。
【図３５】１Ｊ９ＧのＸ線構造とモデル構造のＭＤ。
【図３６】１ＭＭＢのアライメント。
【図３７】１ＭＭＢの立体構造。
【図３８】１ＭＭＢの立体構造。
【図３９】１ＭＭＢのＸ線構造とモデル構造の相違をｒｍｓで表示。
【図４０】１Ｊ９ＧのＸ線構造とモデル構造の基準振動解析。
【図４１】１Ｊ９ＧのＸ線構造とモデル構造の基準振動解析。
【図４２】１ＭＭＢのＸ線構造とモデル構造のＭＤ。
【図４３】ジヒドロ葉酸還元酵素の立体構造。
【図４４】１ＢＺＦ（ＭＯＤＥＬ１８）の基準振動解析。
【図４５】１ＢＺＦ（ＭＯＤＥＬ１８）のＭＤ。
【図４６】１ＬＵＤ（ＭＯＤＥＬ４）より得られた構造活性相関情報。
【図４７】１ＢＺＦ（ＭＯＤＥＬ１８）における活性部位・リガンド結合解析。
【図４８】１ＢＺＦ（ＭＯＤＥＬ４）・リガンド結合。
【図４９】１ＢＺＦ（ＭＯＤＥＬ４）・リガンド結合。
【図５０】１ＢＺＦ（ＭＯＤＥＬ４）・リガンド結合。
【図５１】ｈｅａｔ  ｓｈｏｃｋ  ｐｒｏｔｅｉｎ  ９０の立体構造。
【図５２】１ＹＥＲの基準振動解析。
【図５３】１ＹＥＲのＭＤ。
【図５４】１ＹＥＴより得られた構造活性相関情報。
【図５５】１ＹＥＲにおける活性部位・リガンド結合解析。
【図５６】１ＹＥＲ・リガンド結合。
【図５７】１ＹＥＲ・リガンド結合。
【図５８】ｍｉｔｏｇｅｎ−ａｃｔｉｖａｔｅｄ  ｐｒｏｔｅｉｎ  ｋｉｎａｓｅの立体構造。
【図５９】１Ａ９Ｕの基準振動解析。
【図６０】１Ａ９ＵのＭＤ。
【図６１】１ＯＵＫより得られた構造活性相関情報。
【図６２】１Ａ９Ｕにおける活性部位・リガンド結合解析。
【図６３】１Ａ９Ｕ・リガンド結合。
【図６４】１Ａ９Ｕ・リガンド結合。
【図６５】１Ａ９Ｕ・リガンド結合。
【図６６】１ＡＩＸの立体構造。
【図６７】ｉｎ  ｓｉｌｉｃｏ  ｓｃｒｅｅｎｉｎｇの結果。
【図６８】タンパク質／リガンド複合体構造。
【図６９】タンパク質／リガンド複合体構造。
【図７０】タンパク質／リガンド複合体構造。
【図７１】タンパク質／リガンド複合体構造。
【図７２】タンパク質／リガンド複合体構造。
【図７３】タンパク質／リガンド複合体構造。
【図７４】タンパク質／リガンド複合体構造。
【図７５】タンパク質／リガンド複合体構造。
【図７６】ＳＡＲＳプロテアーゼの立体構造。
【図７７】１ＵＫ３（Ｂ鎖）の基準振動解析。
【図７８】１ＵＫ３（Ｂ鎖）のＭＤ。
【図７９】１ＵＫ４（Ｂ鎖）より得られた構造活性相関情報。
【図８０】１ＵＫ３（Ｂ鎖）におけるｉｎ  ｓｉｌｉｃｏスクリーニングの結果。
【図８１】１ＵＫ３と１ＵＫ４との比較。
【図８２】ｉｎ  ｓｉｌｉｃｏスクリーニングの順位１。
【図８３】ｉｎ  ｓｉｌｉｃｏスクリーニングの順位１。
【図８４】１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。
【図８５】１ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。
【図８６】ＳＡＲ３ヵ所指定で実行したｉｎ  ｓｉｌｉｃｏスクリーニングの結果。
【図８７】１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。
【図８８】ＳＡＲ５ヵ所指定で実行したハイスループットスクリーニングの結果。
【図８９】１ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。
【図９０】１ＵＫ３（Ｂ鎖）より得られた構造活性相関情報。
【図９１】リガンド原子タイプ指定変更で実行したハイスループットスクリーニングの結果。
【図９２】１ＵＫ３（Ｂ鎖）と１ＵＫ４（Ｂ鎖）との比較。
【図９３】１ＵＫ４（Ｂ鎖）より得られた構造活性相関情報。
【図９４】受容体を固定した状態で実行したハイスループットスクリーニングの結果。
【図９５】１ＵＫ３と１ＵＫ４を重ね合わせたリガンドと計算結果のリガンドとの比較。
【図９６】１ＡＸＪにおける二面角拘束ＭＤパラメータの分布。
【図９７】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図９８】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図９９】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１００】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１０１】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１０２】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１０３】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１０４】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１０５】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１０６】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１０７】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１０８】１ＬＵＤ（ＭＯＤＥＬ１）のＭＤ。
【図１０９】異なる条件で受容体／リガンドの結合結果。
【図１１０】異なる条件で受容体／リガンドの結合結果。
【図１１１】異なる条件で受容体／リガンドの結合結果。
【図１１２】１ＢＺＦ用に変更した構造活性相関情報。
【図１１３】１ＢＺＦ（ＭＯＤＥＬ１８）のリガンド結合解析。
【図１１４】１ＢＺＦ（ＭＯＤＥＬ１８）のリガンド結合解析。
【図１】
2005242493.tif 000028
【図２】
2005242493.tif 000029
【図３】
2005242493.tif 000030
【図４】
2005242493.tif 000031
【図５】
2005242493.tif 000032
【図６】
2005242493.tif 000033
【図７】
2005242493.tif 000034
【図８】
2005242493.tif 000035
【図９】
2005242493.tif 000036
【図１０】
2005242493.tif 000037
【図１１】
2005242493.tif 000038
【図１２】
2005242493.tif 000039
【図１３】
2005242493.tif 000040
【図１４】
2005242493.tif 000041
【図１５】
2005242493.tif 000042
【図１６】
2005242493.tif 000043
【図１７】
2005242493.tif 000044
【図１８】
2005242493.tif 000045
【図１９】
2005242493.tif 000046
【図２０】
2005242493.tif 000047
【図２１】
2005242493.tif 000048
【図２２】
2005242493.tif 000049
【図２３】
2005242493.tif 000050
【図２４】
2005242493.tif 000051
【図２５】
2005242493.tif 000052
【図２６】
2005242493.tif 000053
【図２７】
2005242493.tif 000054
【図２８】
2005242493.tif 000055
【図２９】
2005242493.tif 000056
【図３０】
2005242493.tif 000057
【図３１】
2005242493.tif 000058
【図３２】
2005242493.tif 000059
【図３３】
2005242493.tif 000060
【図３４】
2005242493.tif 000061
【図３５】
2005242493.tif 000062
【図３６】
2005242493.tif 000063
【図３７】
2005242493.tif 000064
【図３８】
2005242493.tif 000065
【図３９】
2005242493.tif 000066
【図４０】
2005242493.tif 000067
【図４１】
2005242493.tif 000068
【図４２】
2005242493.tif 000069
【図４３】
2005242493.tif 000070
【図４４】
2005242493.tif 000071
【図４５】
2005242493.tif 000072
【図４６】
2005242493.tif 000073
【図４７】
2005242493.tif 000074
【図４８】
2005242493.tif 000075
【図４９】
2005242493.tif 000076
【図５０】
2005242493.tif 000077
【図５１】
2005242493.tif 000078
【図５２】
2005242493.tif 000079
【図５３】
2005242493.tif 000080
【図５４】
2005242493.tif 000081
【図５５】
2005242493.tif 000082
【図５６】
2005242493.tif 000083
【図５７】
2005242493.tif 000084
【図５８】
2005242493.tif 000085
【図５９】
2005242493.tif 000086
【図６０】
2005242493.tif 000087
【図６１】
2005242493.tif 000088
【図６２】
2005242493.tif 000089
【図６３】
2005242493.tif 000090
【図６４】
2005242493.tif 000091
【図６５】
2005242493.tif 000092
【図６６】
2005242493.tif 000093
【図６７】
2005242493.tif 000094
【図６８】
2005242493.tif 000095
【図６９】
2005242493.tif 000096
【図７０】
2005242493.tif 000097
【図７１】
2005242493.tif 000098
【図７２】
2005242493.tif 000099
【図７３】
2005242493.tif 000100
【図７４】
2005242493.tif 000101
【図７５】
2005242493.tif 000102
【図７６】
2005242493.tif 000103
【図７７】
2005242493.tif 000104
【図７８】
2005242493.tif 000105
【図７９】
2005242493.tif 000106
【図８０】
2005242493.tif 000107
【図８１】
2005242493.tif 000108
【図８２】
2005242493.tif 000109
【図８３】
2005242493.tif 000110
【図８４】
2005242493.tif 000111
【図８５】
2005242493.tif 000112
【図８６】
2005242493.tif 000113
【図８７】
2005242493.tif 000114
【図８８】
2005242493.tif 000115
【図８９】
2005242493.tif 000116
【図９０】
2005242493.tif 000117
【図９１】
2005242493.tif 000118
【図９２】
2005242493.tif 000119
【図９３】
2005242493.tif 000120
【図９４】
2005242493.tif 000121
【図９５】
2005242493.tif 000122
【図９６】
2005242493.tif 000123
【図９７】
2005242493.tif 000124
【図９８】
2005242493.tif 000125
【図９９】
2005242493.tif 000126
【図１００】
2005242493.tif 000127
【図１０１】
2005242493.tif 000128
【図１０２】
2005242493.tif 000129
【図１０３】
2005242493.tif 000130
【図１０４】
2005242493.tif 000131
【図１０５】
2005242493.tif 000132
【図１０６】
2005242493.tif 000133
【図１０７】
2005242493.tif 000134
【図１０８】
2005242493.tif 000135
【図１０９】
2005242493.tif 000136
【図１１０】
2005242493.tif 000137
【図１１１】
2005242493.tif 000138
【図１１２】
2005242493.tif 000139
【図１１３】
2005242493.tif 000140
【図１１４】
2005242493.tif 000141