以下に、本発明の実施形態に係るエネル� ��ー関数最適化システム及び方法、並びにプ� ��グラムとエネルギー関数最適化方法を利用� ��た分子シミュレーションシステム及び方法� ��並びにプログラムについて、図1乃至図18を� ��いて説明する。

 <第1の実施形態>
 第1の実施形態におけるエネルギー関数最適 化システム及び方法について説明する。

 本実施形態におけるエネルギー関数最適� ��システムを図1に示す。本システムは、キー ボード等の入力装置1と、プログラム制御に� �り動作する処理装置2と、ディスプレイ装置� ��印刷装置などの出力装置3と、を含む。

 処理装置2は、多変数最適化指標最適解探 索部201を含む。多変数最適化指標最適解探索 部201は、エネルギー関数の最適解を探索する 。その際、多変数最適化指標を最適化指標と して用いる。多変数最適化指標最適解探索部 201は、分子の安定構造とその分子の不安定構 造とにおけるエネルギーを用いる第1の指標� �、複数種類の分子の安定構造におけるエネ� �ギーを用いる第2の指標との多変数関数であ� ��多変数最適化指標を生成する。そして、生� ��した多変数最適化指標を評価し、それが最� ��解であるか否かを判断する。最適解を探索� ��終えたら,多変数最適化指標最適解探索部201 は最適解を与えるパラメータを含むエネルギ ー関数を出力装置3に出力する。最適解を与� �るパラメータを含むエネルギー関数とは、� �適化されたエネルギー関数のことである。

 本システムは、第1の指標と第2の指標の� �変数関数である多変数最適化指標を最適化� �標として前記エネルギー関数の最適解を探� �し、エネルギー関数の最適化を行う。した� �って、本システムは第1と第2の指標の両方を 考慮してエネルギー関数の最適化を行ってい る。

 このため、本システムによって最適化さ� ��たエネルギー関数は、第1及び第2の指標の� �ずれか一方のみで最適化されたエネルギー� �数よりも汎用性が高く、優れている。すな� �ち、分子シミュレーションを実行するにあ� �って、第1及び第2の指標のいずれか一方のみ で最適化されたエネルギー関数を使用するよ りも、本システムによって最適化されたエネ ルギー関数を使用する方が実際の分子の挙動 を正確に表現でき、予測精度の高いシミュレ ーションを行うことができる。

 本システムによって最適化されたエネル� ��ー関数を用いて行う分子シミュレーション� ��高い予測精度で実際の分子の挙動を表現で� ��る理由について説明する。第1の指標の関数 をf(p)、第2の指標の関数をg(p)、多変数関数最 適化指標であるf(p)とg(p)の多変数関数をh(f,g)� ��する。pはエネルギー関数に含まれるパラメ ータを表す。ここで、f(p)の最適解を与える� �ラメータ{p1}と、g(p)の最適解を与えるパラメ ータ{p2}は、それぞれh(h,g)空間のある断面に� �ける局所的最適解である。対して、本シス� �ムによって求められたh(f,g)の最適解を与え� �パラメータ{p3}は、h(h,g)空間における大域的� ��適解である。それゆえ、本システムによっ� ��得られるエネルギー関数は、第1と第2の指� �両方を最大限に考慮して最適化されている� �よって、分子シミュレーションにおいて、� �システムによって最適化されたエネルギー� �数を使用する方が、第1及び第2の指標のいず れか一方のみで最適化されたエネルギー関数 を使用するよりも、シミュレーションの予測 精度が高くなる。

 次に、各構成について説明する。

 多変数最適化指標最適解探索部201で用い� ��第1と第2の指標は、第1の実施形態に係るシ� ��テムに含まれておらず、本システムの外か� ��入力される。例えば、第1及び第2の指標を� �成する外部システムを本システムとは別に� �意して、その外部システムに第1及び第2の指 標を生成させる。多変数最適化指標最適解探 索部201は外部システムから第1と第2の指標を� ��け取った後、計算等を行ったりハッシュテ� ��ブルを用いたりして、多変数最適化指標を� ��成する。或いは、ユーザーが入力装置1を利 用して第1と第2の指標を本システムに入力し� ��も良い。

 分子の安定構造とその分子の不安定構造� ��におけるエネルギーを用いる第1の指標とは 、分子の安定構造とその分子の不安定構造と におけるエネルギーを、エネルギー関数の最 適化の際に考慮するという意味である。同様 に、複数種類の分子の安定構造におけるエネ ルギーを用いる第2の指標とは、複数種類の� �子の安定構造におけるエネルギーを、エネ� �ギー関数の最適化の際に考慮するという意� �である。

 第1と第2の指標は、既知の分子情報を利� �して求められる。第1の実施形態に係るシス� ��ムには既知の分子情報は含まれておらず、� ��知の分子情報を本システムの外から取得す� ��。既知の分子情報は、第1と第2の指標の生� �に用いることができれば良く、その情報の� �得方法は特に限定されない。例えば、第3の� ��施形態にて説明するように、システム内に� ��知の分子情報を記憶する記憶装置を設けて� ��良い。或いは、ユーザーが入力装置1を利用 して分子構造に関する情報を本システムに入 力しても良い。

 既知の分子情報とは、分子構造及びエネ� ��ギー値またはエネルギー由来の物理量等を� ��す。

 本システムで使用するエネルギー関数は� ��分子シミュレーション用であって経験的手� ��によって最適化されるものであれば良い。� ��えば、ドッキングシミュレーション用のソ� ��トウェアであるFlexXに含まれるエネルギー� �数を用いることができる。

 第1の指標は、分子の安定構造とその分子 の不安定構造とにおけるエネルギーを考慮し て求められ、最適化指標として利用できるも のであればその求め方は特に限定しない。例 えば、式1によってzスコアを計算し、求めら� ��たzスコアを第1の指標として用いることが� �きる。

式1中に含まれるEexpは分子の安定構造におけ� ��エネルギーを表す。<Ecalc>は分子の不安 定構造におけるエネルギーの平均値を表す。 σcalcは分子の不安定構造におけるエネルギー の標準偏差を表す。

 zスコアを計算する場合、分子の不安定構 造におけるエネルギーの平均値及び標準偏差 を計算するために、分子の不安定構造は複数 あるのが望ましい。

 分子の安定構造とその分子の不安定構造� ��含まれるデータセットが複数ある場合は、� ��々の分子について上記zスコアを計算し、そ れらzスコアの平均を第1の指標として用いる� ��が望ましい。

 第1の指標として、分子の安定構造と不安 定構造におけるエネルギーの差、RMSD等を用� �ることも可能である。

 第2の指標は、複数種類の分子の安定構造 におけるエネルギーを考慮して求められ、最 適化指標として利用できるものであればその 求め方は特に限定しない。例えば、結合状態 である分子の結合エネルギーの実験値と、エ ネルギー関数によって求められる分子の結合 構造の計算値との相関係数を第2の指標とし� �用いることができる。別の例として、第1の� ��標の例として挙げたzスコアがある。

 多変数最適化指標最適解探索部201は、第1 と第2の指標の多変数関数である多変数最適� �指標を求めて評価するものであれば良く、� �の関数の形は特に限定しない。ユーザーが� �意に一方の最適化指標に重みをつけたりす� �こともできる。多変数最適化指標の例とし� �、第1と第2の指標の2変数関数である、第1の� ��標と第2の指標の二乗との積が挙げられる。

 多変数最適化指標最適解探索部201は、多� ��数最適化指標を用いて、本実施形態で用い� ��エネルギー関数の最適解を探索するもので� ��れば良く、その探索方法は特に限定しない� ��例えば、エネルギー関数の最適解を与えう� ��全パラメータに対して計算を行うことによ� ��て、エネルギー関数の最適解を求めること� ��できる。

 また、最適解候補を与えるパラメータの� ��が膨大で最適解探索のための時間を抑えた� ��場合は、最適化アルゴリズムを用いて近似� ��を求めて、その近似解を最適解としても良� ��。最適化アルゴリズムは、シミュレーティ� ��・アニーリング手法、遺伝的アルゴリズム� ��タブ・サーチ手法、シミュレーティド・エ� ��リューション、確率的進化手法など、一般� ��なものであれば良く、用いる最適化アルゴ� ��ズムは特に限定しない。

 前述した方法を用いて最適解を探索する� ��合、必ずしも全工程を終えてから最適解を� ��める必要はない。最適解を探索する途中で� ��索を中断し、それまでに求めた解の中から� ��も良い解を最適解として用いても良い。

 エネルギー関数に含まれるパラメータは� ��数の場合もある。その場合は、全てのパラ� ��ータを考慮してエネルギー関数の最適化を� ��う。

 次に、図1及び図2を参照して本システム� �動作について説明する。図2は本システムの� ��作を示す。

 入力装置1から実行の命令が与えられると 、処理装置2は、多変数最適化指標最適解探� �部201によって、第1の指標と第2の指標を用い て多変数最適化指標の生成し評価する(ステ� �プ501)。

 そして、ステップ501の動作は、多変数最� ��化指標最適解探索部201によって、最適化が� ��了するまで繰り返される(ステップ502)。最� �化の終了はユーザーが選んだ最適解探索方� �によって異なる。エネルギー関数の最適解� �与えうる全パラメータに対して計算を行い� �ネルギー関数の最適解を探索する場合は、� �てのパラメータに対して計算し終えたとき� �最適化が終了する。最適化アルゴリズムを� �いてエネルギー関数の近似解を求める場合� �、ユーザーが指定した最適化終了条件を満� �したときに最適化が終了する。

 最適化が終了した後、多変数最適化指標� ��適解探索部201が最適化されたエネルギー関� ��を出力装置3に出力する(ステップ503)。

 以上、説明したとおり、本システムは、� ��1の指標と第2の指標の多変数関数である多� �数最適化指標を最適化指標として、前記エ� �ルギー関数の最適解を探索し、エネルギー� �数の最適化を行う。本システムは、タンパ� �質等の高分子の分子シミュレーションにお� �て特に有効である。タンパク質等の高分子� �は、その立体構造が非常に複雑である。そ� �ため、第1と第2の指標のいずれかのみを用い て最適化されたエネルギー関数では、実際の 分子の挙動を正確に再現するのは難しい。予 測精度の高い分子シミュレーションを行うた めには、第1及び第2の指標両方を考慮してエ� ��ルギー関数を最適化するのが望ましい。本� ��ステムは第1及び第2の指標の両方を用いて� �ネルギー関数を最適化することによって、� �頼性の高い分子シミュレーションを実現す� �。

 <第2の実施形態>
 第2の実施形態におけるエネルギー関数最適 化システム及び方法であって、処理装置内に 第1の指標生成部と第2の指標生成部を含むも� ��について説明する。本実施形態は、第1の実 施形態の応用であるため、第1の実施形態と� �様な点については説明を省略する。

 本実施形態と第1の実施形態の異なる点に ついて図3、図4を使って説明する。図3は本実 施形態におけるエネルギー関数最適化システ ムを示す。図4は本システムの動作を示す。

 処理装置2は、第1の指標生成部202と、第2� ��指標生成部203と、多変数最適化指標最適解� ��索部201と、を含む。第1の指標生成部202はエ ネルギー関数の第1の指標を生成し、求めら� �た第1の指標を多変数最適化指標最適解探索� ��201に出力する。第2の指標生成部203はエネル ギー関数の第2の指標を生成し、求められた� �2の指標を多変数最適化指標最適解探索部201� ��出力する。多変数最適化指標最適解探索部2 01は第1の指標と第2の指標を受け取り、多変� �最適化指標を求め、評価する。なお、ここ� �言う第1の指標及び第2の指標とは、第1の実� �形態で説明したものと同じである。

 本システムの処理装置2は、第1の指標生� �部202と第2の指標生成部203とを含む。多変数� ��適化指標最適解探索部201は、第1の指標の値 と第2の指標の値とを用いて多変数最適化指� �を求めて評価する。第1の指標生成部202と第2 の指標生成部203を本システムに含むことによ って、第1の指標の値と第2の指標の値を本シ� ��テム内で準備することができる。そのため� ��多変数最適化指標の評価及び多変数最適化� ��標最適解の探索を円滑に行うことができる� ��

 ただし、第2の実施形態のシステムは第1� �実施形態のシステムと同様に既知の分子情� �は含んでいないため、第1の指標生成部202と� ��2の指標生成部203は、第1及び第2の指標を生� ��するにあたって、既知の分子情報を本シス� ��ムの外から取得する。或いは、ユーザーが� ��力装置1を利用して本システムに既知の分子 情報を入力しても良い。

 次に、図3及び図4を参照して本システム� �動作について説明する。

 入力装置1から実行の命令が与えられると 、処理装置2は、第1の指標生成部202によって� ��一つのパラメータに対して第1の指標を生成 する(ステップ500A)。そして、生成された値を 多変数最適化指標最適解探索部201に出力する 。

 第2の指標生成部203は、ステップ500Aで用� �たパラメータと同じものに対して第2の指標� ��生成する(ステップ500B)。そして、生成され� ��値を多変数最適化指標最適解探索部201に出� ��する。

 次に、多変数最適化指標最適解探索部201� ��、ステップ500A、ステップ500Bで得られた第1� ��び第2の指標の値を用いて多変数最適化指標 を求め、評価する(ステップ501)。

 図4のステップ500Aからステップ501までの� �作は、多変数最適化指標最適解探索部201に� �って、最適化が終了するまで繰り返される(� ��テップ502)。

 そして、最適化が終了した後、多変数最� ��化指標最適解探索部201が最適化されたエネ� ��ギー関数を出力装置3に出力する(ステップ50 3)。

 なお、上述の動作説明においては第1の指 標を生成(ステップ500A)した後に第2の指標を� �成(ステップ500B)したが、生成する順序はど� �らが先でも構わない。第2の指標を生成した� ��に第1の指標を生成しても良い。或いは、第 1と第2の指標の生成を並列に実行しても良い� ��

 以上、説明したとおり、本システムの多� ��数最適化指標最適解探索部201は第1の指標生 成部202と第2の指標生成部203とを含む。これ� �よって、第1の指標の値と第2の指標の値を本 システム内で準備することができ、多変数最 適化指標の評価及び多変数最適化指標最適解 の探索を円滑に行うことができる。

 <第3の実施形態>
 第3の実施形態におけるエネルギー関数最適 化システム及び方法であって、記憶装置を含 むものについて説明する。本実施形態は、第 1及び第2の実施形態の応用であるため、第1及 び第2の実施形態と同様な点については説明� �省略する。

 本実施形態と第1及び第2の実施形態の異� �る点について図5、図6を使って説明する。図 5は本実施形態におけるエネルギー関数最適� �システムを示す。図6は本システムの動作を� ��す。

 本システムは記憶装置4を含む。記憶装置 4は、トレーニングデータ記憶部401と探索過� �発生データ記憶部402を含む。トレーニング� �ータ記憶部401は、第1の指標生成部202及び第2 の指標生成部203でトレーニングデータとして 用いる既知の分子情報を記憶している。第1� �び第2の指標生成部は各指標を生成する際に� ��レーニングデータ記憶部401から既知の分子� ��報を読み込む。探索過程発生データ記憶部4 02は、多変数最適化指標最適解探索部201が多� ��数最適化指標の最適解を探索する際に生じ� ��一時的な情報を記憶する。多変数最適化指� ��最適解探索部201は必要に応じて探索過程発� ��データ記憶部402に記憶されている情報を読� ��取る。

 本システムは、トレーニングデータ記憶� ��401を含む。トレーニングデータ記憶部401に� ��憶されるトレーニングデータとは、経験的� ��法によりエネルギー関数を最適化する際に� ��用する情報である。システムにトレーニン� ��データ記憶部401を含むことによって、トレ� ��ニングデータをシステム外に準備したり、� ��ーザーが入力したりする必要がなくなる。� ��たがって、本システムでは、第1の指標生成 部202及び第2の指標生成部203による第1及び第2 の指標の生成を円滑に行うことができる。

 本システムは、探索過程発生データ記憶� ��402を含む。探索過程発生データ記憶部402は� ��多変数最適化指標最適解探索部201で行う多� ��数最適化指標の最適解の探索時における一� ��的な情報を保存する。例えば最適化アルゴ� ��ズムを用いて近似解を求める場合、探索中� ��おいて、あるパラメータを用いて生成され� ��解が最適解であるか否かが判断できない。� ��こで、そのようなパラメータは最適解候補� ��与えるものとして記憶しておく必要がある� ��そのようなパラメータの情報量が少量であ� ��ば処理装置2内に記憶しても良いが、莫大に なる場合には、探索過程発生データ記憶部402 に保存するのが望ましい。

 上述の方法で近似解を求める場合、最適� ��が終了した後に探索過程発生データ記憶部4 02に記憶されているデータの中から最も良い� ��を与えるパラメータを抽出すれば良い。

 トレーニングデータ記憶部401にはトレー� ��ングデータとして既知の分子情報が記憶さ� ��る。これは、第1の指標及び第2の指標を求� �るために利用される情報である。すなわち� �分子の安定構造とその不安定構造とにおけ� �エネルギーや、複数種類の分子の安定構造� �おけるエネルギーの情報である。例えば、� �数種類の分子の結合エネルギーの実験値や� �ネルギー関数によって求められる複数種類� �分子の結合構造の計算値、或いはそれらを� �めるための情報や、結合後の分子の安定構� �におけるエネルギーや結合後の分子の不安� �構造におけるエネルギー、或いはそれらを� �めるための情報が挙げられる。

 探索過程発生データ記憶部402は本システ� ��に含まれなくても良い。例えば、エネルギ� ��関数の最適解を与えうる全パラメータに対� ��て計算を行ってエネルギー関数の最適解を� ��索する際に、あるパラメータに対して求め� ��れた最適解が、その直前のパラメータに対� ��て求められた最適解よりも良かった場合に� ��み、記録してある最適解を与えるパラメー� ��の情報を更新するようにする。全てのパラ� ��ータに対して前述のように計算することに� ��って、探索過程において発生する一時的な� ��報を最小限に抑えることが可能となる。こ� ��ような場合においては、探索過程発生デー� ��記憶部402を本システムから除くことができ� ��。

 次に、図5及び図6を参照して本実施形態� �動作について説明する。

 入力装置1から実行の命令が与えられると 、処理装置2は第1の指標生成部202により、ト� ��ーニングデータ記憶部401からトレーニング� ��ータを読み込み(ステップ500a)、エネルギー� ��数を用いて一つのパラメータに対して第1の 指標を生成する(ステップ500A)。そして、第1� �指標生成部202は、生成した第1の指標を多変� ��最適化指標最適解探索部201に出力する。

 第2の指標生成部203は、トレーニングデー タ記憶部401からトレーニングデータを読み込 み(ステップ500b)、エネルギー関数を用いてス テップ500Aと同じパラメータに対して第2の指� ��を生成する(ステップ500B)。そして、第2の指 標生成部203は、生成した第2の指標を多変数� �適化指標最適解探索部201に出力する。

 ステップ501aからステップ501までの動作は 、多変数最適化指標最適解探索部201によって 、最適化が終了するまで繰り返される(ステ� �プ502)。その際、多変数最適化指標最適解探� ��部201は必要に応じて探索中に生じる一時的� ��情報を探索過程発生データ記憶部402に記憶� ��る。

 以上、説明したとおり、本システムは記� ��装置4を含む。記憶装置4は、トレーニング� �ータ記憶部401と探索過程発生データ記憶部40 2とを含む。本システムはトレーニングデー� �記憶部401をシステム内に含むことによって� �トレーニングデータをシステム外に準備し� �り、ユーザーが入力したりする必要がなく� �る。したがって、本システムでは、第1の指� ��生成部202及び第2の指標生成部203による第1� �び第2の指標の生成を円滑に行うことができ� ��。また、本システムは探索過程発生データ� ��憶部402をシステム内に含むことによって、� ��適化アルゴリズムを用いて近似解を求める� ��に発生する大量の一時的な情報を保存する� ��とが可能となる。

 <第4の実施形態>
 第4の実施形態では、上述したエネルギー関 数最適化方法によって最適化されたエネルギ ー関数を利用する分子シミュレーションシス テム及び方法について説明する。本実施形態 は、第1、第2、及び第3の実施形態の応用であ るため、第1、第2、及び第3の実施形態と同様 な点については説明を省略する。

 本実施形態と第1、第2、及び第3の実施形� ��の異なる点について図7、図8を使って説明� �る。図7は本実施形態に係るエネルギー関数� ��適化方法によって最適化されたエネルギー� ��数を利用する分子シミュレーションシステ� ��を示す。図8は本システムの動作を示す。

 本システムの処理装置2は分子シミュレー ション部204を含む。分子シミュレーション部 204は、エネルギー関数最適化方法によって最 適化されたエネルギー関数403を読み込む。そ して、分子シミュレーション部204は、その最 適化されたエネルギー関数を用いて分子シミ ュレーションを実行する。分子シミュレーシ ョン部204は分子シミュレーションを実行後、 その結果を出力装置3に出力する。

 本システムは、エネルギー関数最適化方� ��によって最適化されたエネルギー関数を利� ��する分子シミュレーション部204を含む。こ� ��によって、本システムは、第1の指標のみで 最適化されたエネルギー関数または第2の指� �のみで最適化されたエネルギー関数のいず� �か一方のみを用いる分子シミュレーション� �ステムよりも、実際の分子の挙動を正確に� �現でき、予測精度の高いシミュレーション� �行うことができる。

 第4の実施形態に係るシステムにはエネル ギー関数を最適化する機能が含まれていない ため、最適化されたエネルギー関数403は、本 システムの外で用意する。例えば、本システ ムとは別のシステムを用意し、その別のシス テムにエネルギー関数の最適化を行わせても 良い。

 最適化されたエネルギー関数403は、本シ� ��テムで用いる分子シミュレーションのため� ��エネルギー関数であって、上述した方法に� ��って最適化されたエネルギー関数であれば� ��でも良い。

 処理装置2には、第1の指標生成部、第2の� ��標生成部、多変数最適化指標最適解探索部� ��トレーニングデータ記憶部、探索過程発生� ��ータ記憶部が含まれていても良い。これら� ��全て本システム内に含まれていれば、エネ� ��ギー関数の最適化から分子シミュレーショ� ��までの一連の動作を本システムのみで実行� ��きる。このような場合においては、最適化� ��れたエネルギー関数403を本システムの外で� ��意する必要はない。

 分子シミュレーション部204で用いる分子� ��ミュレーション技術は、経験的手法を用い� ��最適化されるエネルギー関数を利用するも� ��であれば特にその方法は限定しない。例と� ��て、ドッキングシミュレーション、分子力� ��法などがある。

 次に、図7及び図8を参照して本実施形態� �動作について説明する。

 入力装置1から実行の命令が与えられると 処理装置2は分子シミュレーション部204によ� �、最適化されたエネルギー関数403を読み込� �(ステップ506)。そして、その最適化されたエ ネルギー関数を用いて分子シミュレーション を実行する(ステップ507)。分子シミュレーシ� ��ン部204は分子シミュレーションの実行後、� ��の結果を出力装置3に出力する(ステップ508)� ��

 以上、説明したとおり、本システムは、� ��述したエネルギー関数最適化方法によって� ��適化されたエネルギー関数を利用する分子� ��ミュレーション部を含んでいる。これによ� ��て、第1の指標のみで最適化されたエネルギ ー関数及び第2の指標のみで最適化されたエ� �ルギー関数のいずれか一方のみを用いる分� �シミュレーションシステムよりも、予測精� �が高いシミュレーションを行うことができ� �。

 <実施例1>
 次に、本発明の実施例を、図9,10を参照して 説明する。

 本実施例では、エネルギー関数の最適化� ��、本発明の実施形態に係る最適化方法、第1 の指標のみを用いる方法、第2の指標のみを� �いる方法の3通りで行い、その結果を比較し� ��。

 本実施例で用いる分子シミュレーション� ��ステムは、第3の実施形態における分子シミ ュレーションシステムと同様の構成を成して いる。

 記憶装置4に含まれるトレーニングセット 記憶部41には、エネルギー関数の最適化に用� ��る38種類のタンパク質の情報が記憶されて� �る。それぞれのタンパク質には、そのタン� �ク質と結合することが既知の化合物(リガン� ��)が存在する。タンパク質とリガンドが結合 することによって形成される複合体は、タン パク質‐リガンド複合体と呼ばれる。それぞ れのタンパク質の情報には、そのタンパク質 ‐リガンド複合体のX線結晶構造(結合構造)、 結合構造生成シミュレーションによって人工 的に生成された非結晶構造、及び結合エネル ギーを求めるための実験データが含まれる。 タンパク質‐リガンド複合体の非結晶構造の 数は各タンパク質によって異なり、その数は 6から92である。

 分子シミュレーション用エネルギー関数� ��、FlexXに含まれるエネルギー関数の一部を� �いた。FlexXのエネルギー関数は以下のように 表される(Rarey, M.ら J. Mol. Biol. 1996, 261, 47 0-89.)。

 δGbindはタンパク質‐リガンド複合体の結合 状態を表す結合構造スコア、Fiは位置に依存� ��る関数を、δGiはエネルギーパラメータを、 σは相互作用に関わる全ての原子ペアの和を� ��す。matchは水素結合と、金属コンタクトと� �芳香族間の相互作用から成る項である。lipo� ��疎水性相互作用を表す項、ambigは極性原子� �非極性原子の相互作用を表す項、clashは原子 の衝突に対するペナルティ項、rotは化合物が タンパク質と結合することによって失うエン トロピー項を表す。nrotは化合物の回転可能� �結合数である。δG0はこのエネルギー関数特� ��の定数項である。

 本実施例で用いたエネルギー関数を以下� ��示す。

 式3は、式2からclashとrotの項を除いたもので ある。これらの項は、実際の分子間の相互作 用エネルギーに直接関わる項ではなく、分子 シミュレーションの予測精度を上げるために 付加された項である。本実施例では、式2か� �clashとrotの項を除いた、相互作用エネルギー に関わる項のみで形成されるエネルギー関数 を用いた。

 本実施例では、エネルギー関数の最適な� ��ラメータを生成するために以下のスコアリ� ��グ関数を定義した。すなわち、エネルギー� ��数の最適化では、このスコアリング関数の� ��適化を行った。

 式4中の(a、b、c)は(δGmatch、δGlipo、δGambig)の スケーリングパラメータである。

 最適化の際には、各スケーリングパラメ� ��タを-3から3まで0.1毎に区切り、全ての組み� ��わせについて計算した。この場合、(a、b、c )の全ての組み合わせは613通りある。

 第1の指標と第2の指標の生成では、タン� �ク質‐リガンド複合体の結合エネルギー、� �験結合構造スコア、非結合構造スコアを用� �た。

 タンパク質‐リガンド複合体の結合エネ� ��ギーとして、結合エネルギーを求めるため� ��実験データの平衡定数から求められる-logKd� ��用いた。

 実験結合構造スコアとして、タンパク質� ��リガンド複合体のX線結晶構造を用いてエネ ルギー関数を解くことによって得られる値を 用いた。

 非結合構造スコアとして、結合構造生成� ��ミュレーションによって人工的に生成され� ��非結晶構造を用いてエネルギー関数を解く� ��とによって得られる値を用いた。

 エネルギー関数の最適化に用いる38個の� �ンパク質を図9に示す。タンパク質の情報はP rotein Data Bank (http://www.wwpdb.org/)に登録され� �いるものを利用した。図9中に記載されたPDB� ��はProtein Data BankのIDを表す。図9中に記載さ れた‐logKdは結合エネルギーを表す。

 第1の指標は、タンパク質‐リガンド複合 体の平均zスコア(zave)とした。zaveを求めるた� ��に、まずは式1に従ってタンパク質‐リガン ド複合体におけるzスコアを計算した。その� �、Eexpとして実験結合構造スコアを用いた。& lt;Ecalc>として非結合構造スコアの平均値を 用いた。σcalcとして非結合構造スコアの標準 偏差を用いた。そして、それぞれのタンパク 質‐リガンド複合体において得られたzスコ� �の平均を求め、zaveを求めた。

 この場合、実験によって得られるタンパ� ��質‐リガンド複合体の構造を正確に再現す� ��ためには、zaveは負に大きい方が良い。した がって、エネルギー関数の最適化はzaveが負� �大きくなるように実行された。

 第2の指標として、複数種類の分子の結合 エネルギーの実験値と、エネルギー関数によ って求められる複数種類の分子の結合構造の 計算値との相関係数を用いた。複数種類の分 子の結合エネルギーの実験値として、タンパ ク質‐リガンド複合体の結合エネルギーを用 いた。エネルギー関数によって求められる複 数種類の分子の結合構造の計算値として、実 験結合構造スコアを用いた。結合エネルギー と実験結合構造スコアのデータ列をそれぞれ x={xi}、y={yi}、(i=1,2,…,38)とすると、相関係数r は次の式で表される。

 xaveとyaveは各データ列の相加平均を表す。

 この場合、実験によって得られるタンパ� ��質‐リガンド複合体の構造を正確に再現す� ��ためには、相関係数は高い方が良い。した� ��って、エネルギー関数の最適化は相関係数� ��高くなるように実行された。

 多変数最適化指標(opt)は、次の式で表さ� �るものを用いた。

 上述したrとzaveの特徴から、実験によって� �られるタンパク質‐リガンド複合体の構造� �正確に再現するためには、optは負に大きい� �が良い。したがって、エネルギー関数の最� �化は、optを負に最大化させるパラメータを� �索するように実行された。

 本実施例では、本発明の実施形態に係る� ��適化指標としてoptを用いてエネルギー関数� ��最適化を行った。そして、比較例1として第 1の指標であるzaveのみを用いてエネルギー関� ��の最適化を行った。加えて、比較例2として 第2の指標であるrのみを用いてエネルギー関� ��の最適化を行った。

 説明の便宜上、optを用いて最適化された� ��ネルギー関数を、本発明の実施形態に係る� ��ネルギー関数と呼ぶ。第1の指標のみを用い て最適化されたエネルギー関数を第1のエネ� �ギー関数と呼ぶ。第2の指標のみを用いて最� ��化されたエネルギー関数を第2のエネルギー 関数と呼ぶ。

 それぞれの最適化の結果について、図10� �用いて説明する。図10は、それぞれの最適化 指標を用いて最適化されたエネルギー関数を 用いて、zave、r、optを求めた結果を示す。図� ��の第1の指標の値はzaveを表す。図中の第2の� ��標の値はrを表す。図中の多変数最適化指標 の値はoptを表す。図中の「本実施形態」は、 本発明の実施形態に係るエネルギー関数を表 す。比較例1は、第1のエネルギー関数を表す� ��比較例2は第2のエネルギー関数を表す。

 比較例1では、他のものと比べて、第1の� �標の値が最も良かった。しかし、この場合� �最適化の際に第2の指標については一切考慮� ��れていないため、第2の指標の値は悪かった 。一方、比較例2では第2の指標の値が最も良� ��った。しかし、この場合、最適化の際に第1 の指標については一切考慮されていないため 、第1の指標の値は悪かった。

 本発明の実施形態に係るエネルギー関数� ��は、比較例1と比較例2に比べて、第1の指標� ��値及び第2の指標の値共に良い結果を得られ た。本発明の実施形態に係るエネルギー関数 の第2の指標の値においては、比較例2におけ� ��第2の指標の値と同等であった。多変数最適 化指標の値は本発明の実施形態に係るエネル ギー関数によって求めたものが最も良かった 。

 第1または第2の指標のみを用いて最適化� �た場合、他の指標は全く考慮されていない� �対して、本発明の実施形態では第1の指標と� ��2の指標の多変数関数を最適化指標として用 いたため、両指標を考慮して最適化が行われ ている。そのため、本発明の実施形態に係る エネルギー関数では両指標において良い結果 を得ることができた。

 したがって、分子シミュレーションにお� ��て本発明に係るエネルギー関数を用いれば� ��第1の指標と第2の指標との両方を考慮して� �ミュレーションを実行することができる。

 以上の結果から、本発明の実施形態によ� ��て、エネルギー関数の第1の指標である、分 子の安定構造とその分子の不安定構造とにお けるエネルギーと、エネルギー関数の第2の� �標である、複数種類の分子の安定構造にお� �るエネルギーと、の両方を考慮してエネル� �ー関数の最適化を行うことが可能となる。� �発明の実施形態によって最適化されたエネ� �ギー関数は前述した二つの要素を考慮して� �適化されているため、それぞれを単独で用� �て最適化されたエネルギー関数よりも、分� �シミュレーションに適している。すなわち� �本システムによって最適化されたエネルギ� �関数を用いることによって、特に、両方の� �ネルギーを考慮する必要のある分子シミュ� �ーションにおいて、高い予測精度が実現で� �る。

 <実施例2>
 本実施例では、実施例1にて最適化された本 発明の実施形態に係るエネルギー関数の精度 を評価するために以下のドッキングシミュレ ーションのテストを行った。

 テストに用いる標的タンパク質としてEstr ogen Receptor(ER)を使用した。また、ERにドッキ� ��グさせる化合物として、結合することが分� ��っている化合物(既知リガンド)32個、結合し ないことが分かっている化合物1340個の計1373� ��を用いた。

 化合物の多数の配座を発生させるプログ� ��ムとして、市販のドッキングシミュレーシ� ��ンソフトウェアFlexSISを用いた。

 評価方法を以下に説明する。

 はじめにFlexSISを用いて全ての化合物に対 して配座を発生させた。

 次に、実施例1において最適化された本発 明の実施形態に係るエネルギー関数を用いて 、ある化合物の全配座に対してERとの結合エ� ��ルギーを計算した。そして、結合エネルギ� ��が最も低い配座を選抜し、その結合エネル� ��ーをERと化合物の結合エネルギーとした。� �意した1373個の全化合物に対して同様の計算� ��し、全化合物に対して結合エネルギーを決� ��した。その後、結合エネルギーが低い順に� ��合物を並べ、化合物に対して昇順に順位付� ��を行った。

 結合時の分子の構造においては結合エネ� ��ギーが低い方が安定であると考えられる。� ��たがって、上述のように化合物に対して順� ��付けを行った場合、上位の化合物ほどERと� �結合状態が安定であるといえる。すなわち� �上位の化合物ほど、シミュレーションによ� �予測上ではERと結合しやすいと考えられる。

 さらに、第1のエネルギー関数、第2のエ� �ルギー関数を用いて、同様の計算及び動作� �行った。

 結果を図11に示す。図11はエンリッチメン トプロットである。x軸は、順位付けされた� �化合物の上位からの割合(%)を表す。y軸は、� ��既知リガンドに対する、全化合物の上位x%� �含まれる既知リガンドの割合(%)を表す。図� �の「本実施形態」は本発明の実施形態に係� �エネルギー関数の結果を表す。図中の比較� �1は第1のエネルギー関数の結果を表す。図中 の比較例2は第2のエネルギー関数の結果を表� ��。

 本実施例の場合、xの値が低いところでy� �値が高くなるような結果を示すエネルギー� �数が最も優れていると評価できる。

 本発明の実施形態に係るエネルギー関数� ��比較例1、比較例2のそれぞれの結果を比較� �たところ、図11に示すとおり、本発明の実施 形態に係るエネルギー関数の結果が最も良か った。すなわち、ERと化合物の結合状態にお� ��るエネルギーを予測して、ERと結合しやす� �化合物の順位付けをそれぞれのエネルギー� �数を用いて行ったところ、最も多くの既知� �ガンドが本発明の実施形態に係るエネルギ� �関数による順位付けの上位に含まれていた� �したがって、本発明の実施形態に係るエネ� �ギー関数を用いて行った分子シミュレーシ� �ンの予測精度が最も高かった。

 以上、説明したとおり、第1の指標のみを 用いて最適化されたエネルギー関数、第2の� �標のみを用いて最適化されたエネルギー関� �、本発明の実施形態に係る最適化方法を用� �て最適化されたエネルギー関数のそれぞれ� �用いて分子シミュレーションを行ったとこ� �、本発明の実施形態に係る最適化方法を用� �て最適化されたエネルギー関数を用いた分� �シミュレーションの予測精度が最も高かっ� �。本発明の実施形態に係る最適化方法を用� �て最適化されたエネルギー関数を分子シミ� �レーションに用いることによって、実際の� �子の挙動をより正確に表現することが可能� �なる。

 <実施例3>
 実施例3では、下記式によって求められるopt を用いて最適化されたエネルギー関数の精度 を評価した。下記式に含まれるrの指数nは任� ��の数である。本実施例では指数nを変化させ ることによるエネルギー関数の性能の変化を 調べた。なお、n=2の場合における下記式は実 施例1で用いたoptを求める式(式6)と同じもの� �なる。

 初めにエネルギー関数の最適化を行った。� ��ネルギー関数は実施例1で用いたものと同じ ものを用いた。エネルギー関数の最適化は実 施例1と同様に行った。ただし、optとして、n= 1、2,3,4,5の場合における上記式を用いた。最� ��化のために用いるタンパク質データとして� ��実施例1で用いた38個のタンパク質データか� ��labe、3fx2、4claを除いた計35個のデータを用� �た。

 次に、最適化したエネルギー関数の精度� ��価を実施例2と同様の方法にて行い、比較例 1,2と本発明の実施形態によって最適化された エネルギー関数の結果とを比較した。

 精度評価の結果を図12~図18に示す。図12~� �16はn=1,2、3,4,5の場合の結果を表すエンリッ� �プロットメントである。図17はn=1,2,3のそれ� �れの場合における順位付けされた全化合物� �上位25%の結果を表すエンリッチプロットメ� �トである。従って、図17は図12、13,14に示す� �位付けされた全化合物の上位25%を抽出した� �である。図18はn=2,4、5のそれぞれの場合にお ける順位付けされた全化合物の上位25%の結果 を表すエンリッチプロットメントである。従 って、図18は図13、15、16に示す順位付けされ� ��全化合物の上位25%を抽出した図である。x軸 は順位付けされた全化合物の上位からの割合 (%)を表す。y軸は、全既知リガンドに対する� �全化合物の上位x%含まれる既知リガンドの割 合(%)を表す。図中のzr1は、本実施形態のエネ ルギー関数であってn=1の場合におけるoptを用 いて最適化されたものの結果である。zr2~zr5� �zr1と同様に、本実施形態のエネルギー関数� �あって、n=2~5の場合におけるoptを用いて最適 化されたものの結果である。比較例1は第1の� ��ネルギー関数の結果、比較例2は第2のエネ� �ギー関数の結果を表す。

 図17、18では、x軸は25%、y軸は50%までしか� ��していない。これは、薬剤開発の初期工程� ��ある化合物スクリーニング等における実用� ��を考慮すると、上位20%程度までの結果が性� ��評価の対象として相応しいからである。

 本発明の実施形態に係るエネルギー関数� ��あってn=1~5のもの、比較例1、比較例2のそれ ぞれの結果を比較したところ、図12~図16に示� ��とおり、本発明の実施形態に係るエネルギ� ��関数であってn=1~5のものの結果が最も良か� �た。すなわち、ERと化合物の結合状態におけ るエネルギー関数を予測して、ERと結合しや� ��い化合物の順位付けをそれぞれのエネルギ� ��関数を用いて行ったところ、最も多くの既� ��リガンドが本発明の実施形態に係るエネル� ��ー関数による順位付けの上位に含まれてい� ��。

 本実施例では、optの計算式(式7)に含まれ� ��rnの指数nを1~5として、エネルギー関数の最� ��化を行った。そして、n=1~5の全ての場合に� �けるエネルギー関数にて、比較例よりも優� �た結果を得ることができた。したがって、� �実施形態のエネルギー関数の第2の指標とし� ��、rnのような単調関数を用いることは有用� �あると考えられる。

 なお、本実施例におけるrnはnが大きくな� ��とrnは小さくなることから、rnは単調減少関 数であると言える。

 加えて、n=1~5の場合におけるそれぞれの� �果を比較したところ、図17、18に示すとおり� ��n=2の場合における結果が最も良かった。

 なお、本実施例では最適化指標としてn=1~ 5におけるopt(式7)を用いたが、nの値はこれら� ��限定されるものではない。例えば、n=6やn=7� ��おけるoptを用いても良い。

 また、最適なoptは、エネルギー関数や用� ��るデータ等によって変わる。従って、ユー� ��は適宜、nの値や第1の指標、第2の指標等を� ��更して、最適なoptを設定するのが望ましい� ��

 実施例3では、本発明の実施形態に係るエ ネルギー関数の最適化方法にて、単調減少関 数rnを第2の指標として用いた。そして、それ によって最適化されたエネルギー関数の性能 と比較例の性能とを比べて、本発明の実施形 態に係る方法に従って最適化されたエネルギ ー関数の方が優れていることを確認した。従 って、本実施形態のエネルギー関数の最適化 方法で第2の指標として単調減少関数を用い� �ことによって、高い性能を有するエネルギ� �関数が得られるようになる。

 以上、実施形態および実施例を参照して� ��発明を説明したが、本発明は実施形態およ� ��実施例に限定されるものではない。クレー� ��に定義された本発明の構成や詳細には、発� ��のスコープ内で当業者が理解し得る様々な� ��更をすることができる。

 この出願は、2008年1月16日に出願された日 本出願特願2008-006426と、2008年4月22日に出願さ れた日本出願特願2008-110871とを基礎として優� ��権の利益を主張するものであり、その開示� ��全てを引用によってここに取り込む。