(57)【特許請求の範囲】
【請求項１】
モデル化されたＮ個の質点から構成された質点系の挙動を予測するシミュレーション装置において、
前記質点系の構造を記述する３Ｎ個の質点座標に基づいて、前記質点系の構造の変化を主として担うＭ個の座標である中継座標を設定し、前記質点系の構造を記述しかつ前記中継座標から独立した座標である非中継座標を設定する座標設定手段と、
前記非中継座標を前記中継座標に従属させることにより前記中継座標の関数としての前記非中継座標の構成を求め、前記中継座標の変化に伴う前記非中継座標の変化が前記中継座標に与える影響を考慮しながら、前記中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される前記一般座標の前記可変成分であるＫ個の集団座標の関数としての前記中継座標の構成を求める座標抽出手段と、
前記集団座標についての運動方程式の解として得られる時間の関数としての前記集団座標、前記中継座標の前記構成および前記非中継座標の前記構成に基づいて、前記質点座標の時間発展を予測する逆変換手段とを備えることを特徴とするシミュレーション装置。
（Ｋ、ＭおよびＮは、Ｋ＜Ｍ＜３Ｎの関係を満たし、それぞれ１以上の整数を表す。）
【請求項２】
前記座標抽出手段が、
前記中継座標および前記非中継座標の関数として表現された前記質点系のポテンシャルエネルギーを求める第１の工程を実施し、
前記ポテンシャルエネルギーを用いた断熱近似条件に従って前記非中継座標を前記中継座標に従属させる第２の工程を実施し、
現状の前記中継座標および前記非中継座標の下で、前記影響を考慮しながら、前記中継座標に関する前記ポテンシャルエネルギーの微分係数を求める第３の工程を実施し、
前記ポテンシャルエネルギーの前記微分係数の下で、自己無撞着集団座標法における前記ポテンシャルエネルギーの微分係数を用いた基礎方程式に従って、前記集団座標に関する前記中継座標の微分係数を求める第４の工程を実施し、
前記中継座標の前記微分係数の下で前記集団座標を微小量だけ更新することにより、更新された前記中継座標を求める第５の工程を実施し、
前記更新された前記中継座標の下で、前記中継座標に従属した前記非中継座標について構造緩和を行う第６の工程を実施し、
その後、前記非中継座標の構造緩和後の状態における前記中継座標および前記非中継座標に基づいて、前記第３の工程から前記第６の工程までを順次繰り返すことにより、前記中継座標の前記構成を求めるものであることを特徴とする請求項１に記載のシミュレーション装置。
【請求項３】
前記影響を考慮する方法が下記式１から式３の少なくとも１つを使用することであることを特徴とする請求項２に記載のシミュレーション装置。
【数１】
0005697638.tif 000052
（ｉ、ｊおよびｋは、それぞれ１からＭまでの整数を表し、
α、βおよびγは、それぞれ１から３Ｎ−Ｍまでの整数を表し、
ＲＳｉは、前記質点系におけるｉ番目の前記中継座標を表し、
ＲＦαは、前記質点系におけるα番目の前記非中継座標を表し、
ＲＳは（ＲＳ１，ＲＳ２，・・・，ＲＳＭ）を表し、
ＲＦは（ＲＦ１，ＲＦ２，・・・，ＲＦ（３Ｎ−Ｍ））を表し、
ＲＦ（ＲＳ）は、前記中継座標に従属させられた前記非中継座標を表し、
Ｖ（ＲＳ，ＲＦ）は、前記中継座標と前記非中継座標で表された前記ポテンシャルエネルギーを表し、
Ｖｅｆｆ（ＲＳ）は、Ｖ（ＲＳ，ＲＦ）にＲＦ（ＲＳ）を代入して得られる有効ポテンシャルエネルギーを表す。
式２において、第３項の（ｉ←→ｊ）は、第２項におけるｉおよびｊの文字を相互に入れ替えた項を表す。
式３において、第３項の（ｉ←→ｋ）は、第２項におけるｉおよびｋの文字を相互に入れ替えた項を表し、第４項の（ｊ←→ｋ）は、第２項におけるｊおよびｋの文字を相互に入れ替えた項を表し、第６項の（ｉ←→ｋ）は、第５項におけるｉおよびｋの文字を相互に入れ替えた項を表し、第７項の（ｊ←→ｋ）は、第５項におけるｊおよびｋの文字を相互に入れ替えた項を表す。
さらに、式１から式３において下記式４を用いた。
【数２】
0005697638.tif 000053
Ｋαβ−１（Ｒｓ）はＫαβ（Ｒｓ）の逆行列を表し、
Ｋαβ（Ｒｓ）およびＪαｉ（Ｒｓ）はそれぞれ式５および式６で定義される。
【数３】
0005697638.tif 000054
）
【請求項４】
前記断熱近似条件が下記式７であることを特徴とする請求項２または３に記載のシミュレーション装置。
【数４】
0005697638.tif 000055
（ＲＦαは、前記質点系におけるα番目の前記非中継座標を表し、
ＲＳは（ＲＳ１，ＲＳ２，・・・，ＲＳＭ）を表し、
ＲＦは（ＲＦ１，ＲＦ２，・・・，ＲＦ（３Ｎ−Ｍ））を表し、
Ｖ（ＲＳ，ＲＦ）は、前記中継座標と前記非中継座標で表された前記ポテンシャルエネルギーを表す。）
【請求項５】
前記集団座標の個数ＫがＫ＝１を満たすものであり、
前記基礎方程式が下記式８および９であることを特徴とする請求項２から４いずれかに記載のシミュレーション装置。
【数５】
0005697638.tif 000056
（ｉおよびｊは、それぞれ１からＭまでの整数を表し、
ＲＳｉは、前記質点系におけるｉ番目の前記中継座標を表し、
ＲＳは（ＲＳ１，ＲＳ２，・・・，ＲＳＭ）を表し、
ｑ１は、前記集団座標を表し、
ｍｉは前記質点系におけるｉ番目の前記中継座標に関する質量を表し、
φｉ（ＲＳ）は、式９を満たす関数（固有ベクトル）のｉ番目の成分を表し、
Λ（ＲＳ）は、式９を満たす関数（固有値）を表し、
Ｖｅｆｆ（ＲＳ）は有効ポテンシャルエネルギーを表す。）
【請求項６】
前記集団座標の個数ＫがＫ＝１を満たすものであり、
前記基礎方程式が下記式１０から１２であることを特徴とする請求項２から４いずれかに記載のシミュレーション装置。
【数６】
0005697638.tif 000057
（ｉ、ｊおよびｋは、それぞれ１からＭまでの整数を表し、
ＲＳｉは、前記質点系におけるｉ番目の前記中継座標を表し、
ＲＳは（ＲＳ１，ＲＳ２，・・・，ＲＳＭ）を表し、
ｑ１は、前記集団座標を表し、
ｍｉは前記質点系におけるｉ番目の前記中継座標に関する質量を表し、
φｉ（ＲＳ，λ）およびκ（ＲＳ，λ）は、式１２を満たす関数のｉ番目の成分を表し、
λは、補助座標を表し、
Ｖｅｆｆ（ＲＳ）は有効ポテンシャルエネルギーを表す。）
【請求項７】
前記座標抽出手段が、前記第４の工程において、前記基礎方程式中の前記集団座標に関する前記中継座標または補助座標の微分係数が有する符号の任意性を解消するように、該基礎方程式を解く上で前記中継座標から独立しかつ前記集団座標の関数として取り扱う変数の数を増やして、計算を行うものであることを特徴とする請求項２から４いずれかに記載のシミュレーション装置。
【請求項８】
前記座標抽出手段が、前記中継座標から独立しかつ前記集団座標の関数として取り扱う変数の数を増やした結果得られた下記式１３で表される基礎方程式に従って計算を行うものであることを特徴とする請求項７に記載のシミュレーション装置。
【数７】
0005697638.tif 000058
（Yは、下記式１４で定義されるＭＫ＋Ｍ＋Ｋ次元のベクトルであり、
ｖμは、下記式１５の非斉次線形方程式の解ベクトルである。
【数８】
0005697638.tif 000059
式１５におけるＣおよびｓμはそれぞれ式１６および式１７で定義される。
【数９】
0005697638.tif 000060
Ｖ,ｉｊ（ＲＳ）およびＶ,ｉｊｋ（ＲＳ）はそれぞれ式１８および式１９で定義される。
【数１０】
0005697638.tif 000061
また、ｉ、ｊおよびｋはそれぞれ１からＭまでの整数を表し、
μおよびνはそれぞれ１からＫまでの整数を表し、
ｑμはμ番目の前記集団座標を表し、
ＲＳｉは、前記質点系におけるｉ番目の前記中継座標を表し、
ＲＳは（ＲＳ１，ＲＳ２，・・・，ＲＳＭ）を表し、
ｍｉは前記質点系におけるｉ番目の前記中継座標に関する質量を表し、
φｉμおよびΛμは、ＲＳから独立した補助座標を表し、
Ｖｅｆｆ（ＲＳ）は有効ポテンシャルエネルギーを表す。）
【請求項９】
前記集団座標の個数ＫがＫ＝１を満たすことを特徴とする請求項８に記載のシミュレーション装置。
【請求項１０】
前記座標抽出手段が、前記中継座標から独立しかつ前記集団座標の関数として取り扱う変数の数を増やした結果得られた下記式２０で表される基礎方程式に従って計算を行うものであることを特徴とする請求項７に記載のシミュレーション装置。
【数１１】
0005697638.tif 000062
（Ｚは、下記式２１で定義されるＭＫ＋Ｍ＋２Ｋ次元のベクトルであり、
ｃμνは、μごとに定義される下記式２２で表される値をそれぞれ最小にするように一意に決められる定数であり、
ｗμは、下記式２３で定義された行列ＤのＫ次元の特異値空間を張るＫ個のＭＫ＋Ｍ＋２Ｋ次元の単位ベクトルである。
【数１２】
0005697638.tif 000063
Ｖ,ｉｊ（ＲＳ）およびＶ,ｉｊｋ（ＲＳ）はそれぞれ式２４および式２５で定義される。
【数１３】
0005697638.tif 000064
また、ｉ、ｊおよびｋはそれぞれ１からＭまでの整数を表し、
μおよびνはそれぞれ１からＫまでの整数を表し、
ｑμはμ番目の前記集団座標を表し、
ＲＳｉは、前記質点系におけるｉ番目の前記中継座標を表し、
ＲＳは（ＲＳ１，ＲＳ２，・・・，ＲＳＭ）を表し、
ｍｉは前記質点系におけるｉ番目の前記中継座標に関する質量を表し、
φｉμ、λμおよびρμは、ＲＳから独立した補助座標を表し、
Ｖｅｆｆ（ＲＳ）は有効ポテンシャルエネルギーを表す。）
【請求項１１】
前記集団座標の個数ＫがＫ＝１を満たすことを特徴とする請求項１０に記載のシミュレーション装置。
【請求項１２】
前記座標抽出手段が、前記ポテンシャルエネルギーの３次微分の項の計算を下記式２６に基づいて行うものであることを特徴とする請求項６から１１いずれかに記載のシミュレーション装置。
【数１４】
0005697638.tif 000065
（φｉは任意のベクトルのｉ番目の成分を表し、
Ｖ,ｉｊ（ＲＳ）、Ｖ,ｉｊｋ（ＲＳ）およびｎはそれぞれ式２７から式２９で定義される。
【数１５】
0005697638.tif 000066
）
【請求項１３】
前記座標設定手段が、前記質点系の構造のうち特徴的な部分構造ごとに抽出された該部分構造を代表する代表座標を前記中継座標として設定するものであることを特徴とする請求項１から１２いずれかに記載のシミュレーション装置。
【請求項１４】
前記質点系が生体高分子を含む多原子系であり、
前記部分構造が、前記生体高分子の２次構造、前記生体高分子の構成単位または前記生体高分子の主鎖であり、
前記部分構造についての前記代表座標が、該部分構造を構成する原子の座標そのもの、該原子の前記座標を組み合わせて定義できる座標または前記部分構造の周期間隔であることを特徴とする請求項１３に記載のシミュレーション装置。
【請求項１５】
前記生体高分子がタンパク質であり、
前記部分構造が前記タンパク質の２次構造であり、
該２次構造についての前記代表座標が、該２次構造を構成する原子群の重心座標または前記２次構造の屈曲角であることを特徴とする請求項１４に記載のシミュレーション装置。
【請求項１６】
前記２次構造がヘリックス構造、βシート、ターン、ループおよびランダムコイルの少なくとも１つであることを特徴とする請求項１５に記載のシミュレーション装置。
【請求項１７】
前記生体高分子がタンパク質であり、
前記部分構造が前記タンパク質の残基であり、
該残基についての前記代表座標が、該残基を構成する原子群の重心座標であることを特徴とする請求項１４に記載のシミュレーション装置。
【請求項１８】
前記生体高分子がタンパク質であり、
前記部分構造が前記タンパク質の主鎖であり、
該主鎖についての前記代表座標が、該主鎖を構成する原子のそれぞれの座標であることを特徴とする請求項１４に記載のシミュレーション装置。
【請求項１９】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の２次構造であり、
該２次構造についての前記代表座標が、該２次構造を構成する原子群の重心座標または前記２次構造の屈曲角であることを特徴とする請求項１４に記載のシミュレーション装置。
【請求項２０】
前記２次構造が螺旋構造であることを特徴とする請求項１９に記載のシミュレーション装置。
【請求項２１】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の残基であり、
該残基についての前記代表座標が、該残基を構成する原子群の重心座標であることを特徴とする請求項１４に記載のシミュレーション装置。
【請求項２２】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の主鎖であり、
該主鎖についての前記代表座標が、該主鎖を構成する原子のそれぞれの座標であることを特徴とする請求項１４に記載のシミュレーション装置。
【請求項２３】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の螺旋構造であり、
該螺旋構造についての前記代表座標が、該螺旋構造の周期間隔であることを特徴とする請求項１４に記載のシミュレーション装置。
【請求項２４】
前記多原子系が、前記生体高分子に対する結合候補分子を含むものであることを特徴とする請求項１４から２３いずれかに記載のシミュレーション装置。
【請求項２５】
請求項１から２４いずれかに記載のシミュレーション装置によって実行される、モデル化されたＮ個の質点から構成された質点系の挙動を予測するシミュレーション方法であって、
前記質点系の構造を記述する３Ｎ個の質点座標に基づいて、前記質点系の構造の変化を主として担うＭ個の座標である中継座標を設定し、
前記質点系の構造を記述しかつ前記中継座標から独立した座標である非中継座標を設定し、
前記非中継座標を前記中継座標に従属させることにより前記中継座標の関数としての前記非中継座標の構成を求め、
前記中継座標の変化に伴う前記非中継座標の変化が前記中継座標に与える影響を考慮しながら、前記中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される前記一般座標の前記可変成分であるＫ個の集団座標の関数としての前記中継座標の構成を求め、
前記集団座標についての運動方程式の解として得られる時間の関数としての前記集団座標、前記中継座標の前記構成および前記非中継座標の前記構成に基づいて、前記質点座標の時間発展を予測することを特徴とするシミュレーション方法。
（Ｋ、ＭおよびＮは、Ｋ＜Ｍ＜３Ｎの関係を満たし、それぞれ１以上の整数を表す。）
【請求項２６】
請求項２５に記載のシミュレーション方法をコンピュータに実行させためのシミュレーションプログラム。
【請求項２７】
請求項２６に記載のシミュレーションプログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、コンピュータを用いて、モデル化された質点系の動力学的な挙動を予測するシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体に関するものである。
【背景技術】
【０００２】
コンピュータ技術の発展に伴って、シミュレーションによる理論的な計算により、タンパク質、核酸、脂質および多糖等の生体高分子の動力学的な大変形の挙動を原子レベルで理論的に解明しようという研究が盛んである。具体的には、標的とするタンパク質と結合候補分子（当該タンパク質に対する結合の適合性を有するか否かについての解析の対象となる分子）との互いの親和性を理論的に予測する創薬スクリーニング、および、タンパク質の一次配列から立体構造を構築する原理を明らかにして一次構造から高次構造を理論的に構築するタンパク質のフォールディング（折り畳み機構）の解析が挙げられる。
【０００３】
例えば、生体高分子の動力学的な挙動を予測するシミュレーション方法として、巨大分子でも原子レベルでシミュレーションを実行することが可能な分子動力学法またはモンテカルロ法等に基づくシミュレーション方法（非特許文献１から７）がある。分子動力学法によれば、運動方程式に従って多原子系の時間発展を微小時間単位で順次追跡することができる。しかし、生体高分子を含む多原子系のポテンシャル曲面には、多くのローカルミニマム（極小点）或いは多くのエネルギー障壁が存在しているため、上記のような方法では、多原子系の状態が初期構造に近いローカルミニマムに捕捉されて、計算に膨大な時間を要するという問題があった。また、この問題はモンテカルロ法でも同じである。
【０００４】
そこで、上記のようなローカルミニマムによるトラップの問題を解消する方法として、分子動力学法に基づく計算において「強制的な力」を印加することによりローカルミニマムによるトラップを脱出させる計算方法が知られている。例えば、特許文献１から３並びに非特許文献８および９では、多原子系が一度経験した構造を二度と取らないように仮想的な相互作用を逐次導入することや、高温相の多原子系の加速された運動で見つかった構造を低温相の多原子系の運動に逐次反映させることが開示されている。このような相互作用の導入等により、多原子系が通過するポテンシャル曲面上の経路の探索を加速して、生体高分子の動力学的な挙動を計算することが可能となる。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００５−２６７５９２号公報
【特許文献２】特開２００５−５２４１２９号公報
【特許文献３】再表２００６／０６８２７１号公報
【非特許文献】
【０００６】
【非特許文献１】Todd J. A. Ewing, Irwin D. Kuntz, Journal of Computational Chemistry, Vol. 18, Issue 9 (1997), p. 1175-1189
【非特許文献２】Garrett M. Morris, David S. Goodsell, et al., Journal of Computational Chemistry, Vol. 19, Issue 14 (1998), p. 1639-1662
【非特許文献３】Matthias Rarey, Bernd Kramer, et al., Journal of Molecular Biology, Vol. 261, Issue 3 (1996), p. 470-489
【非特許文献４】Ruben Abagyan, Maxim Totrov, et al., Journal of Computational Chemistry, Vol. 15, Issue 5 (1994), p. 488-506
【非特許文献５】Gareth Jones, Peter Willett, et al., Journal of Molecular Biology, Vol. 267, Issue 3 (1997), p. 727-748
【非特許文献６】Richard A. Friesner, Jay L. Banks, et al., Journal of Medicinal Chemistry, Vol. 47, Issue 7 (2004), p. 1739-1749
【非特許文献７】Thomas A. Halgren, Robert B. Murphy, et al., Journal of Medicinal Chemistry, Vol. 47, Issue 7 (2004), p. 1750-1759
【非特許文献８】Yoshifumi Fukunishi, Yoshiaki Mikami, et al., Journal of Physical Chemistry B, Vol. 107, Issue 47 (2003), p. 13201-13210
【非特許文献９】Yuji Sugita, Yuko Okamoto, Chemical Physics Letters, Vol. 314, Issue 1-2 (1999), p. 141-151
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかしながら、生体高分子等がモデル化された質点系の運動を記述する運動方程式に基づいて、時間に関する積分で当該質点系の挙動を予測する分子動力学を用いた計算では、シミュレーションの計算時間および計算精度がトレードオフの関係になってしまうという問題がある。具体的には、特許文献１から３並びに非特許文献８および９のようなシミュレーション方法では、計算時間は短縮されるものの、印加した「強制的な力」によって質点系の本来有する運動秩序が破壊され、計算結果として非現実的な大変形の挙動を算出してしまう。そして、より効率よくローカルミニマムによるトラップを脱出させるために「強制的な力」の作用を大きくすればするほど、その傾向は顕著となる。
【０００８】
本発明は上記問題に鑑みてなされたものであり、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することを可能とするシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体を提供することを目的とするものである。
【課題を解決するための手段】
【０００９】
上記課題を解決するために、本発明に係るシミュレーション装置は、
モデル化されたＮ個の質点から構成された質点系の挙動を予測するシミュレーション装置において、
質点系の構造を記述する３Ｎ個の質点座標に基づいて、質点系の構造の変化を主として担うＭ個の座標である中継座標を設定し、質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定する座標設定手段と、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるＫ個の集団座標の関数としての中継座標の構成を求める座標抽出手段と、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測する逆変換手段とを備えることを特徴とするものである。
【００１０】
なお本明細書において、Ｋ、ＭおよびＮは、Ｋ＜Ｍ＜３Ｎの関係を満たし、それぞれ１以上の整数を表す。
【００１１】
本明細書において、「質点系の構造」とは、当該質点系を構成するＮ個の質点が造り出す立体構造を意味する。
【００１２】
「質点系の構造の変化を主として担う中継座標」とは、当該座標の変化が当該質点系の構造の立体構成に大きな影響を与えるような座標を意味する。
【００１３】
「中継座標を設定し」とは、中継座標として、一部の質点座標そのもの、質点座標を組み合わせて定義できる座標又はこれらの組み合わせを設定することを意味する。非中継座標についても同様である。
【００１４】
そして、本発明に係るシミュレーション装置において、座標抽出手段は、
中継座標および非中継座標の関数として表現された質点系のポテンシャルエネルギーを求める第１の工程を実施し、
ポテンシャルエネルギーを用いた断熱近似条件に従って非中継座標を中継座標に従属させる第２の工程を実施し、
現状の中継座標および非中継座標の下で、上記影響を考慮しながら、中継座標に関するポテンシャルエネルギーの微分係数を求める第３の工程を実施し、
ポテンシャルエネルギーの上記微分係数の下で、自己無撞着集団座標法におけるポテンシャルエネルギーの微分係数を用いた基礎方程式に従って、集団座標に関する中継座標の微分係数を求める第４の工程を実施し、
中継座標の上記微分係数の下で集団座標を微小量だけ更新することにより、更新された中継座標を求める第５の工程を実施し、
更新された中継座標の下で、中継座標に従属した非中継座標について構造緩和を行う第６の工程を実施し、
その後、非中継座標の構造緩和後の状態における中継座標および非中継座標に基づいて、第３の工程から第６の工程までを順次繰り返すことにより、中継座標の上記構成を求めるものであることが好ましい。
【００１５】
「現状の中継座標および非中継座標の下で」とは、１回目の第３の工程においては、座標設定手段により設定された状態における中継座標および非中継座標の下で、２回目以降の第３の工程においては、直前の第６の工程で行われた非中継座標の構造緩和後の状態における中継座標および非中継座標の下で、その第３の工程を実施することを意味する。
【００１６】
また、本発明に係るシミュレーション装置において、上記影響を考慮する方法は下記式１から式３の少なくとも１つを使用することであることが好ましい。
【００１７】
【数１】
0005697638.tif 000002
【００１８】
本明細書において、ｉ、ｊおよびｋは、それぞれ１からＭまでの整数を表し、
α、βおよびγは、それぞれ１から３Ｎ−Ｍまでの整数を表し、
ＲＳｉは、質点系におけるｉ番目の中継座標を表し、
ＲＦαは、質点系におけるα番目の非中継座標を表し、
ＲＳは（ＲＳ１，ＲＳ２，・・・，ＲＳＭ）を表し、
ＲＦは（ＲＦ１，ＲＦ２，・・・，ＲＦ（３Ｎ−Ｍ））を表し、
ＲＦ（ＲＳ）は、中継座標に従属させられた非中継座標を表し、
Ｖ（ＲＳ，ＲＦ）は、中継座標と非中継座標で表された質点系のポテンシャルエネルギーを表し、
Ｖｅｆｆ（ＲＳ）は、Ｖ（ＲＳ，ＲＦ）にＲＦ（ＲＳ）を代入して得られる有効ポテンシャルエネルギーを表す。
【００１９】
式２において、第３項の（ｉ←→ｊ）は、第２項におけるｉおよびｊの文字を相互に入れ替えた項（つまり、第２項におけるｉをｊへ、ｊをｉへ書き換えた項。以下同様である。）を表す。
【００２０】
式３において、第３項の（ｉ←→ｋ）は、第２項におけるｉおよびｋの文字を相互に入れ替えた項を表し、第４項の（ｊ←→ｋ）は、第２項におけるｊおよびｋの文字を相互に入れ替えた項を表し、第６項の（ｉ←→ｋ）は、第５項におけるｉおよびｋの文字を相互に入れ替えた項を表し、第７項の（ｊ←→ｋ）は、第５項におけるｊおよびｋの文字を相互に入れ替えた項を表す。
【００２１】
さらに、式１から式３において下記式４を用いた。
【００２２】
【数２】
0005697638.tif 000003
【００２３】
Ｋαβ−１（Ｒｓ）はＫαβ（Ｒｓ）の逆行列を表し、
Ｋαβ（Ｒｓ）およびＪαｉ（Ｒｓ）はそれぞれ式５および式６で定義される。
【００２４】
【数３】
0005697638.tif 000004
【００２５】
また、本発明に係るシミュレーション装置において、断熱近似条件は下記式７であることが好ましい。
【００２６】
【数４】
0005697638.tif 000005
【００２７】
また、本発明に係るシミュレーション装置において、集団座標の個数ＫはＫ＝１を満たすものであり、自己無撞着集団座標法における基礎方程式は下記式８および９であることが好ましい。
【００２８】
【数５】
0005697638.tif 000006
【００２９】
本明細書において、ｑ１は、集団座標を表し、
ｍｉは質点系におけるｉ番目の中継座標に関する質量を表し、
φｉ（ＲＳ）は、式９を満たす関数（固有ベクトル）のｉ番目の成分を表し、
Λ（ＲＳ）は、式９を満たす関数（固有値）を表す。
【００３０】
或いは、本発明に係るシミュレーション装置において、集団座標の個数ＫはＫ＝１を満たすものであり、自己無撞着集団座標法における基礎方程式は下記式１０から１２であることが好ましい。
【００３１】
【数６】
0005697638.tif 000007
【００３２】
φｉ（ＲＳ，λ）およびκ（ＲＳ，λ）は式１２を満たす関数であり、φｉ（ＲＳ，λ）はそのｉ番目の成分を表し、
λは、補助座標（中継座標から独立しかつ集団座標の関数として取り扱う変数）を表す。
【００３３】
或いは、本発明に係るシミュレーション装置において、座標抽出手段は、第４の工程において、基礎方程式中の集団座標に関する中継座標または補助座標の微分係数が有する符号の任意性を解消するように、この基礎方程式を解く上で中継座標から独立しかつ集団座標の関数として取り扱う変数（つまり補助座標）の数を増やして、計算を行うものであることが好ましい。
【００３４】
また、補助座標の数を増やす場合において、座標抽出手段は、その結果得られた下記式１３で表される基礎方程式に従って計算を行うものであることが好ましい。さらにこの場合において集団座標の個数ＫはＫ＝１を満たすことが好ましい。
【００３５】
【数７】
0005697638.tif 000008
【００３６】
Yは、下記式１４で定義されるＭＫ＋Ｍ＋Ｋ次元のベクトルであり、ｖμは、下記式１５の非斉次線形方程式の解ベクトルである。
【００３７】
【数８】
0005697638.tif 000009
【００３８】
式１５におけるＣおよびｓμはそれぞれ式１６および式１７で定義される。
【００３９】
【数９】
0005697638.tif 000010
【００４０】
Ｖ,ｉｊ（ＲＳ）およびＶ,ｉｊｋ（ＲＳ）はそれぞれ式１８および式１９で定義される。
【００４１】
【数１０】
0005697638.tif 000011
【００４２】
なお、μおよびνはそれぞれ１からＫまでの整数を表し、
ｑμはμ番目の集団座標を表し、
φｉμおよびΛμは補助座標を表す。
【００４３】
或いは、補助座標の数を増やす場合において、座標抽出手段は、その結果得られた下記式２０で表される基礎方程式に従って計算を行うものであることが好ましい。さらにこの場合において集団座標の個数ＫはＫ＝１を満たすことが好ましい。
【００４４】
【数１１】
0005697638.tif 000012
【００４５】
Ｚは、下記式２１で定義されるＭＫ＋Ｍ＋２Ｋ次元のベクトルであり、ｃμνは、μごとに定義される下記式２２で表される値をそれぞれ最小にするように一意に決められる定数であり、ｗμは、下記式２３で定義された行列ＤのＫ次元の特異値空間を張るＫ個のＭＫ＋Ｍ＋２Ｋ次元の単位ベクトルである。なお、λμおよびρμは、φｉμと同様に補助座標を表す。
【００４６】
【数１２】
0005697638.tif 000013
【００４７】
また、本発明に係るシミュレーション装置において、座標抽出手段は、ポテンシャルエネルギーの３次微分の項の計算を下記式２４に基づいて行うものであることが好ましい。φｉは任意のベクトルのｉ番目の成分を表し、ｎは式２５で定義される。
【００４８】
【数１３】
0005697638.tif 000014
【００４９】
また、本発明に係るシミュレーション装置において、座標設定手段は、質点系の構造のうち特徴的な部分構造ごとに抽出された当該部分構造を代表する代表座標を中継座標として設定するものであることが好ましい。
【００５０】
本明細書において、「特徴的な部分構造」とは、質点系の構造のうち、形態的および／または機能的な特徴を有する部分構造を意味する。
【００５１】
また、本発明に係るシミュレーション装置において、質点系は生体高分子を含む多原子系であり、
上記部分構造は、生体高分子の２次構造、生体高分子の構成単位または生体高分子の主鎖であり、
上記部分構造についての代表座標は、その部分構造を構成する原子の座標そのもの、その原子の座標を組み合わせて定義できる座標またはその部分構造の周期間隔とすることができる。
【００５２】
また、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子はタンパク質であり、上記部分構造はタンパク質の２次構造であり、この２次構造についての代表座標は、その２次構造を構成する原子群の重心座標または２次構造の屈曲角であることが好ましい。この場合において、２次構造はヘリックス構造、βシート、ターン、ループおよびランダムコイルの少なくとも１つであることが好ましい。
【００５３】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子がタンパク質であり、上記部分構造がタンパク質の残基であり、この残基についての代表座標が、その残基を構成する原子群の重心座標であることが好ましい。
【００５４】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子がタンパク質であり、上記部分構造がタンパク質の主鎖であり、この主鎖についての代表座標が、その主鎖を構成する原子のそれぞれの座標であることが好ましい。
【００５５】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の２次構造であり、この２次構造についての代表座標が、その２次構造を構成する原子群の重心座標または２次構造の屈曲角であることが好ましい。この場合において、２次構造は螺旋構造であることが好ましい。
【００５６】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の残基であり、この残基についての代表座標が、その残基を構成する原子群の重心座標であることが好ましい。
【００５７】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の主鎖であり、この主鎖についての代表座標が、その主鎖を構成する原子のそれぞれの座標であることが好ましい。
【００５８】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の螺旋構造であり、この螺旋構造についての代表座標が、その螺旋構造の周期間隔であることが好ましい。
【００５９】
また、本発明に係るシミュレーション装置において、多原子系はさらに、生体高分子に対する結合候補分子を含むものとすることができる。
【００６０】
本発明に係るシミュレーション方法は、
上記に記載のシミュレーション装置によって実行される、モデル化されたＮ個の質点から構成された質点系の挙動を予測するシミュレーション方法であって、
質点系の構造を記述する３Ｎ個の質点座標に基づいて、質点系の構造の変化を主として担うＭ個の座標である中継座標を設定し、
質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定し、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、
中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるＫ個の集団座標の関数としての中継座標の構成を求め、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測することを特徴とするものである。
【００６１】
本発明に係るシミュレーションプログラムは、上記に記載のシミュレーション方法をコンピュータに実行させることを特徴とするものである。
【００６２】
本発明に係るコンピュータ読み取り可能な記録媒体は、上記に記載のシミュレーションプログラムを記録したものである。
【発明の効果】
【００６３】
本発明に係るシミュレーション装置は、上記のような座標設定手段と、座標抽出手段と、逆変換手段とを備え、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【００６４】
本発明に係るシミュレーション方法は、上記シミュレーション装置によって実行される方法であって、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【００６５】
本発明に係るプログラムおよび記録媒体は、上記シミュレーション方法を実行させることが可能であるから、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【図面の簡単な説明】
【００６６】
【図１】本発明のシミュレーション方法により誘導適合を考慮した場合における、タンパク質と結合候補分子の結合過程を示す概略図である。
【図２】タンパク質および結合候補分子を含む多原子系の誘導適合を考慮した概念的なポテンシャル曲面における反応経路を示す概略図である。
【図３】集団座標の抽出を利用して多原子系の反応経路を求める手順を概念的に示す図である。
【図４】本発明における変数変換および逆変換の概念を示す図である。
【図５】実施形態のシミュレーション装置の構成を概略的に示すブロック図である。
【図６Ａ】実施形態のシミュレーション方法における計算手順を概略的に示すフローチャート図である。
【図６Ｂ】実施形態のシミュレーション方法における計算手順を概略的に示すフローチャート図である。
【図７】実施形態における中継座標ＲＳ、非中継座標ＲＦおよび原子の座標ｘの関係を示す概略図である。
【図８】実施形態のシミュレーション方法によって得られた所定の複合体の反応経路を示すグラフである。
【図９】図８のグラフにおける複合体の結合過程を示す概略図である。図９ａは複合体の結合過程における開始状態を示し、図９ｂは複合体の結合過程における最終状態を示す。
【図１０】符号の任意性が解消された計算の場合と符号の任意性が残ったままの計算の場合との比較結果を示すグラフである。
【発明を実施するための形態】
【００６７】
以下、本発明の実施形態について図面を用いて説明するが、本発明はこれに限られるものではない。なお、視認しやすくするため、図面中の各構成要素の縮尺等は実際のものとは適宜異ならせてある。
【００６８】
まず、詳細な実施形態の説明を行う前に、従来技術に対して本発明が有する技術上の意義を明らかにするべく、本発明の基本となる技術的思想およびそれに至った経緯を概説する。ここでは、説明を分かりやすくするために、具体的に生体高分子および結合候補分子からなる複合体を含む多原子系の動力学的な挙動を予測する場合を例にして説明する。
【００６９】
上記のような多原子系の動力学的な挙動の予測は、例えば、開発の期間短縮およびコスト削減等の観点から新薬開発において重要な役割を果たす。薬物等の生理活性物質は特定のタンパク質に結合することによりその化学的性質を発揮する。したがって、新薬開発においては、標的とするタンパク質に結合しやすい薬物候補を数十万から数百万とも言われる膨大な化合物から絞り込む必要がある（創薬スクリーニング）。また究極的には、薬物候補を膨大な化合物から数個程度にまで絞り込むことが要求されている。そこで、薬物候補を膨大な化合物から効果的に絞り込むためには多原子系の挙動を原子レベルで取り扱い、分子間さらには原子間の相互作用を正確に評価する必要がある。
【００７０】
図１は、タンパク質の動力学的な挙動を考慮した場合における、タンパク質および結合候補分子の結合過程を示す概略図である。タンパク質２は、結合候補分子６と結合するためのポケット４を有している。つまり、当該結合候補分子６は、当該タンパク質２の生理活性物質と言える。しかしながら、当該ポケット４は、必ずしもそのままで結合候補分子６がその中に収まることができるような形状を有しているわけではない（図１ａ）。例えば図１では、ポケット４の開口が結合候補分子６の大きさよりも小さいことが示されている。そこで、このような場合には、タンパク質２は、接近した結合候補分子６との相互作用に応じて、開口の形状および大きさを結合候補分子６の形状および大きさに適合させながらその構造を変化させる（図１ｂおよびｃ）。
【００７１】
このように、タンパク質および核酸等の生体高分子が生理活性物質との相互作用によりその構造を変化させる現象は「誘導適合」と言われる。分子間の相互作用を正確に計算するためには、当然この誘導適合を考慮しなければならない。そして、これは、多原子系を原子レベルで取り扱うことによって実現できる。
【００７２】
一方、図２は、多原子系の誘導適合を考慮した概念的なポテンシャル曲面上の反応経路を示す概略図である。図２の横軸Ｘ１はタンパク質の構造の変化を概念的に示す軸であり、縦軸Ｘ２はタンパク質および結合候補分子の距離を示す軸である。つまり、図２は、タンパク質の構造の変化および分子間の距離に応じた多原子系のポテンシャル曲面を表している。図２に示されるように、誘導適合が生じる系は、ポテンシャル曲面上の２つの安定点（つまり、複合体が解離体である状態を示す点Ａおよび複合体の結合が完了した状態を示す点Ｃ）を結ぶ反応経路ＲＰ上にエネルギー障壁Ｂ（鞍点）を有することが解かる。
【００７３】
したがって、誘導適合を考慮しながらタンパク質２と結合候補分子６との結合過程を解析するということは、タンパク質２および結合候補分子６に関する上記ポテンシャル曲面上において、エネルギー障壁Ｂを超えて２つの安定点ＡおよびＣを結ぶ反応経路を探索することに帰着する。
【００７４】
以上を考慮すると、多原子系の動力学的な挙動を予測するシミュレーションにおいて計算精度を向上させかつ計算時間を短縮するためには、多原子系の原子レベルでの取り扱いを可能にすること、および非経験的な反応経路の探索を容易にすることが必要となる。
【００７５】
しかしながら、数千から数百万個の原子で構成される多原子系の挙動は、１秒〜１時間オーダーの時間規模で起こる緩慢な大変形であるため、従来のシミュレーション方法では、原子レベルでの取り扱いは可能であるものの、理論的な計算を行うために数十年もの膨大な時間が必要であるという問題があった。そこで、多原子系の挙動の理論的な計算を実用的な期間内で実施できるようにするために、様々なシミュレーション手法が検討・開発されている。前述した「強制的な力」を印加する分子動力学法もその１つであるが、その問題点は既に述べた通りである。また、標的とするタンパク質を剛体として近似して、取り扱う変数を減らすことで計算量を減らす手法も検討されている。しかし、このような粗い近似の下では、タンパク質の動力学的な挙動の影響を当然考慮できず、結果としてこれらの分子間に働く相互作用を正しく計算できない。このような場合、例えば１００万個ある候補の化合物から薬物候補として１万個程度に絞り込むことしかできない。また、その他のシミュレーション手法でも計算精度および計算時間のトレードオフの問題は解決するに至っていない。
【００７６】
そこで、本発明者は、生体高分子を含む多原子系をＮ個の質点から構成される質点系にモデル化した上で、この質点系の挙動を集団運動として取り扱い、質点座標（質点系の構造を一次的に記述する質点の座標）の自由度より小さい自由度の集団座標を、その質点座標に基づいた集団運動の中から抽出することに着目した。
【００７７】
（集団座標）
集団座標について説明する。一般的に「集団座標」とは、集団変数の座標要素を意味する。そして、「集団変数」とは、質点系の正準変数（その質点系における質点座標および質点の運動量）と正準変換により対応づけられる一般変数（ｑ，ｐ）のうち、その正準変数よりも自由度が小さいものを意味する。なお、正準変換により対応づけられる一般変数（ｑ，ｐ）において、ｑは座標要素であり、ｐは運動量要素である。ｑおよびｐは下記式２６および式２７によってそれぞれ定義される。
【００７８】
【数１４】
0005697638.tif 000015
【００７９】
式２６および式２７において、ηは１から３Ｎまでの整数を表す。
【００８０】
したがって、集団変数は、下記式２８および式２９によってそれぞれ定義されるｑ†およびｐ†を用いて（ｑ†，ｐ†）と表すことができる。
【００８１】
【数１５】
0005697638.tif 000016
【００８２】
式２８および式２９は、集団変数（ｑ†，ｐ†）が式２８および式２９に示されたμ＝１〜Ｋまでの２Ｋ個の変数により構成されることを表している。言い換えれば、集団変数（ｑ†，ｐ†）とは、時間に依存して変化する可変成分（ｑ１，ｑ２，・・・，ｑＫ，ｐ１，ｐ２，・・・，ｐＫ）および時間変化に対して定数となる不変成分（ｑＫ＋１＝０，ｑＫ＋２＝０，・・・，ｑ３Ｎ＝０，ｐＫ＋１＝０，ｐＫ＋２＝０，・・・，ｐ３Ｎ＝０）から構成される一般座標の当該可変成分であると言える。なお、式２８および式２９では、不変成分が定数としてゼロの値を取る場合を示しているが、定数はゼロに限られない。
【００８３】
そして、集団座標は、集団変数（ｑ†，ｐ†）の座標要素であるから、ｑ†と表すことができる。より具体的には、集団座標とは、（ｑ１，ｑ２，・・・，ｑＫ）から構成される変数の集合である。なお、集団座標自体は、質点座標のみから正準点変換により求めることもできる。
【００８４】
上記のように、集団運動として取り扱う対象となる座標（例えば、質点座標および後述する中継座標）からより自由度の低い集団座標へと変数変換を行うことを「変数分離」ともいう。
【００８５】
集団座標の自由度は、変数分離の対象となる座標の自由度よりも小さければ、特に限定されない。しかしながら、集団座標の自由度が小さい程、その後の計算処理が容易となるため、集団座標の自由度は１（つまりＫ＝１）であることが好ましい。
【００８６】
（集団運動としての取り扱い）
質点座標に対して変数分離を実施すると、集団座標ｑ†によって、質点系の挙動における本質的な運動、つまり、より低次元的な集団運動の記述が可能となる。この結果、３Ｎ次元の質点座標による運動の記述がＫ次元の集団座標ｑ†による運動の記述へと置き換えられ、質点系の運動の記述が単純化される。具体的には、質点座標で記述された運動方程式は、集団座標ｑ†を引数として持つ質点座標ｘの構成ｘ＝ｘ（ｑ†）とｑ†についての運動方程式との組み合わせに単純化される。なお、ｘは（ｘ１，ｘ２，・・・，ｘ３Ｎ）を表す。質点座標ｘ（ｑ†）の具体的な構成は質点系の集団運動としての挙動の仕方を表し、集団座標ｑ†についての運動方程式はその挙動の基本法則を表す。したがって、運動方程式を解く上で、取り扱う変数の自由度がより減少することになるため、理論的な計算が容易になる。
【００８７】
そして、質点座標ｘに対して変数分離を実施して集団座標ｑ†に関するポテンシャル曲面上の反応経路を見つけ出し、座標を集団座標ｑ†から質点座標ｘに逆変換して、質点座標ｘに関するポテンシャル曲面上で反応経路を表現することにより、質点系の反応経路を求めることができる。上記のような方法によれば、低次元の運動方程式を解いた後、集団座標ｑ†から質点座標ｘへの変数変換を実行するだけであるため、本質的に質点系の状態がローカルミニマムに捕捉されることが起こり得ない。
【００８８】
例えば、図２に示された複合体を含む多原子系の反応経路を求める場合の概念的な手順は以下のようになる。図３は、集団座標の抽出を利用して多原子系の反応経路を求める手順を概念的に示す図である。まず、質点座標Ｘ１およびＸ２に関するポテンシャル曲面（図３ａ）が、抽出された集団座標ｑ１およびｑ２に関するポテンシャル曲面（図３ｂ）に変換される。この際、集団座標として、２つの安定点ＡおよびＣ並びにエネルギー障壁Ｂが一直線上に並ぶような集団座標を採用することが重要である。このような集団座標ｑ１およびｑ２に関するポテンシャル曲面上では、２つの安定点ＡおよびＣを結ぶ反応経路ＲＰを一気に構成することができる（図３ｂ）。その後、座標の逆変換により、変数Ｘ１およびＸ２に関するポテンシャル曲面上での反応経路ＲＰが求められる（図３ｃ）。ここで得られた一次元的な軌跡こそ理論化学の分野で古くから用いられてきた「反応座標」と同定できる変数である。
【００８９】
（生体高分子を含む多原子系の挙動を集団運動として取り扱う際の課題）
変数分離は集団運動理論に基づいて行われる。例えば、集団運動理論の１つとして、自己無撞着集団座標法（ＳＣＣ法：Self-consistent Collective Coordinate method）がある。ＳＣＣ法は、現在原子核の分野で研究が進められている物理学における集団運動理論の１つである。より具体的には、ＳＣＣ法は、系の集団運動を記述する本質的な集団変数を求め、その系のハミルトニアンに含まれる正準変数が張る６Ｎ次元の位相空間から、動力学的な集団運動に関与する分離多様体を見出し、その分離多様体上またはその近傍で定義されたハミルトニアンによって集団運動を記述する運動理論である。「分離多様体」とは、その系のハミルトニアンに含まれる正準変数が張る６Ｎ次元の位相空間中の部分空間であって、その系の挙動を表す軌道が閉じ込められる部分空間を意味する。つまり、分離多様体中の任意の１点を初期値とした場合、上記軌道は常にこの分離多様体中に拘束されることになる。ＳＣＣ法に関する詳細については、Sin-itiro Tomonaga, Progress of Theoretical Physics, Vol. 13, No. 5 (1955), p. 482-496、Toshio Marumori, Toshihide Maskawa, et al., Progress of Theoretical Physics, Vol. 64, No. 4 (1980), p. 1294-1314、およびGiu Do Dang, Abraham Klein, et al., Physics Reports, Vol. 335, Issues 3-5 (2000), p. 93-274（以下、DangおよびKleinらの文献を非特許文献１０という。）等を参照されたい。
【００９０】
そこで、例えば原子核反応において取り扱う核子（陽子および中性子）を化学反応において取り扱う原子に置き換えてＳＣＣ法を適用することにより、多原子系の動力学的な挙動を集団運動として取り扱ってこの多原子系の反応経路を求めることができるようにも思える。
【００９１】
しかし、本発明者は、本発明の対象である生体高分子を含むような大規模な系に対してＳＣＣ法等の既存の集団運動理論をそのまま適用するだけでは、計算精度を向上させかつ計算時間を短縮するには不十分であるという問題を新たに見出した。これは、原子核反応の場合にはせいぜい数百個の核子を取り扱えばよいのに対して、水中における生体高分子を含む多原子系の化学反応の場合には数千から数百万個の原子を取り扱わなければならないことに起因する。つまり、集団運動理論とはいえ既存の手法では変数の制御、つまり変数分離の計算が複雑になるためである。
【００９２】
そして、上記のような問題は、複合体を含む多原子系の反応経路を求める場合に限られるものではない。つまり、より一般的に、生体高分子を含む系の反応経路を求める場合に同様の問題が生ずることは容易に想像される。
【００９３】
「発明の説明」
以上のような経緯を踏まえ、本発明者は、質点座標の自由度より小さい自由度の集団座標によって集団運動を記述する集団運動理論であって、巨大な原子群である生体高分子を含む多原子系の構造計算にも適用可能な新規な集団運動理論に基づいて、質点系の動力学的な挙動を予測するシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体を発明するに至った。
【００９４】
より具体的には、本発明に係るシミュレーション装置は、
モデル化されたＮ個の質点から構成された質点系の挙動を予測するシミュレーション装置において、
質点系の構造を記述する３Ｎ個の質点座標に基づいて、質点系の構造の変化を主として担うＭ個の座標である中継座標を設定し、質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定する座標設定手段と、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるＫ個の集団座標の関数としての中継座標の構成を求める座標抽出手段と、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測する逆変換手段とを備えることを特徴とするものである。
【００９５】
そして、本発明に係るシミュレーション方法は、
上記に記載のシミュレーション装置によって実行される、モデル化されたＮ個の質点から構成された質点系の挙動を予測するシミュレーション方法であって、
質点系の構造を記述する３Ｎ個の質点座標に基づいて、質点系の構造の変化を主として担うＭ個の座標である中継座標を設定し、
質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定し、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、
中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるＫ個の集団座標の関数としての中継座標の構成を求め、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測することを特徴とするものである。
【００９６】
また、本発明に係るシミュレーションプログラムは、上記に記載のシミュレーション方法をコンピュータに実行させることを特徴とするものである。
【００９７】
また、本発明に係るコンピュータ読み取り可能な記録媒体は、上記に記載のシミュレーションプログラムを記録したものである。
【００９８】
なお、「Ｎ個の質点から構成された質点系」とは、質点系の構成要素である質点の総数がＮ個であることを意味する。
【００９９】
＜座標設定手段＞
座標設定手段は、質点系の構造を記述する３Ｎ個の質点座標に基づいて、質点系の構造の変化を主として担うＭ個の座標である中継座標を設定し、質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定するものである。上記「３Ｎ個の質点座標」とは、より具体的には３次元空間におけるＮ個の質点の座標ｘである。
【０１００】
まず、座標設定手段は、質点座標の構成に基づいて、質点系の巨視的な集団運動（大振幅集団運動）を記述する上で質点系の構造の変化を主として担うＭ個の中継座標を設定する。そして、後述する座標抽出手段はこの設定された中継座標に対して変数分離を実施することになる。言い換えれば、本発明では、質点系を記述する独立した座標のうち大振幅集団運動による質点系の構造変化に対して影響の少ない座標は既知であるので、このような座標は集団運動理論における変数分離の対象から予め除外される。つまり、３Ｎ次元の質点座標から集団座標を抽出するのではなく、Ｍ次元の中継座標から集団座標を抽出することにより、変数分離の計算が単純化される。Ｍは、Ｋ＜Ｍ＜３Ｎを満たす整数である。また、Ｋは、上記式２８および式２９での説明同様、集団座標ｑ†の要素の個数である。
【０１０１】
なお、この中継座標からさらに変数の階層化を行ってもよい。つまり、上記のように設定された中継座標の構成に基づいて、さらに低次元の２次的な中継座標を設定し、この２次的な中継座標から集団座標を抽出することも可能である。
【０１０２】
中継座標は、大振幅集団運動による質点系の構造変化に対する影響の度合を基準に設定される。そして、この影響の度合が大きい座標ほど中継座標としてふさわしい。言い換えれば、どのような質点系の大振幅集団運動を解析の対象としているのかによって、中継座標の設定方法および具体的な座標の内容が変わる。
【０１０３】
具体的には、中継座標は、質点系に含まれる一部の質点の座標そのもの、質点の座標を組み合わせて定義できる座標又はこれらの組み合わせを用いて設定される。また中継座標は、例えば質点系の構造のうち特徴的な部分構造に関連して設定されることが好ましい。そしてこの場合において、中継座標は、特徴的な部分構造ごとに抽出された、当該部分構造を代表する代表座標であることがより好ましい。代表座標は、１つの特徴的な部分構造に複数あってもよい。
【０１０４】
「特徴的な部分構造」とは、質点系の構造のうち、形態的および／または機能的な特徴を有する部分構造を意味する。特徴的な部分構造は、例えば質点系が生体高分子を含む多原子系である場合には、生体高分子の２次構造（部分的な折り畳み構造）、生体高分子の構成単位および生体高分子の主鎖である。また、特徴的な部分構造の「代表座標」とは、その特徴的な部分構造を代表して表す座標を意味する。代表座標は、例えば、その特徴的な部分構造を構成する質点（原子）の座標そのもの、その質点（原子）の座標を組み合わせて定義できる座標またはその部分構造の周期間隔等である。
【０１０５】
より具体的には、上記のような特徴的な部分構造としては、生体高分子がタンパク質である場合には、タンパク質の２次構造が挙げられる。具体的な２次構造としては、ヘリックス構造（３１０ヘリックス、αヘリックス、πヘリックスおよびβヘリックス等）、βシート、ターン、ループおよびランダムコイル等が挙げられる。そして、タンパク質の２次構造についての代表座標としては、その２次構造を構成する原子群の重心座標または２次構造の屈曲角が挙げられる。例えば、一般的なタンパク質中に含まれる上記のような特徴的な部分構造（高次構造）の個数は、せいぜい数個から数十個である。これにより、変数分離の対象となる変数の個数を大幅に低減することができる。また、生体高分子と結合候補分子との結合反応における挙動を予想する場合には、結合候補分子の構造自体を多原子系の特徴的な部分構造の１つとして取り扱うことも可能である。
【０１０６】
或いは、特徴的な部分構造としては、生体高分子がタンパク質である場合には、タンパク質の残基（タンパク質の構成単位であるアミノ酸部分。Ｎ末端残基およびＣ末端残基を含む）とすることもできる。そして、タンパク質の残基についての代表座標としては、その残基を構成する原子群の重心座標が挙げられる。
【０１０７】
或いは、特徴的な部分構造としては、生体高分子がタンパク質である場合には、タンパク質の主鎖とすることもできる。そして、タンパク質の主鎖についての代表座標としては、主鎖を構成する原子のそれぞれの座標が挙げられる。
【０１０８】
或いは、特徴的な部分構造としては、生体高分子が核酸である場合には、核酸の２次構造とすることもできる。具体的な２次構造としては、螺旋構造等が挙げられる。螺旋構造の場合には、所定の周期ごとに特徴的な部分構造として取り扱う。そして、核酸の２次構造についての代表座標としては、その２次構造を構成する原子群の重心座標または２次構造の屈曲角が挙げられる。
【０１０９】
或いは、特徴的な部分構造としては、生体高分子が核酸である場合には、核酸の残基（核酸の構成単位であるヌクレオチド部分）とすることもできる。そして、核酸の残基についての代表座標としては、その残基を構成する原子群の重心座標が挙げられる。
【０１１０】
或いは、特徴的な部分構造としては、生体高分子が核酸である場合には、核酸の主鎖とすることもできる。そして、核酸の主鎖についての代表座標としては、主鎖を構成する原子のそれぞれの座標が挙げられる。
【０１１１】
或いは、生体高分子が核酸でありかつ核酸の螺旋構造である場合には、螺旋構造の周期間隔を代表座標とすることもできる。
【０１１２】
特徴的な部分構造を実際にどのように設定するかについては、シミュレーションを行う対象の質点系の構成および／または現象に応じて適宜決定される。例えば、タンパク質の折りたたみ構造の形成過程や核酸の螺旋構造の形成過程を解析の対象とする場合には、生体高分子の構成単位や生体高分子の基本的な骨格とも言える主鎖を特徴的な部分構造として採用することが好ましい。生体高分子の主鎖については、共有結合により連なった一本の主鎖全体を１つの特徴的な部分構造としてもよいし、上記一本の主鎖を分割して得られた部分的な構造のそれぞれを特徴的な部分構造としてもよい。また、タンパク質と結合候補分子との結合反応における挙動を解析の対象にする場合には、上記で挙げた主鎖の他、へリックス構造、βシート、ターン、ループおよびランダムコイル等のそれぞれ、又はこれらの組み合わせを特徴的な部分構造として採用することもできる。特徴的な部分構造は複数の２次構造を含んでもよい。つまり、２次構造のそれぞれを１つの特徴的な部分構造とすることもでき、主鎖に沿って隣接した複数の２次構造をまとめて１つの特徴的な部分構造とすることもできる。また、核酸の螺旋構造に阻害剤が挟まって結合する反応を解析の対象にする場合には、上記で挙げた主鎖の他、当該螺旋構造の周期構造、又はこれらの組み合わせを特徴的な部分構造として採用することもできる。螺旋構造については、全体を１つの特徴的な部分構造としてもよいし、螺旋構造を所定の周期で分割して得られた部分的な構造のそれぞれを特徴的な部分構造としてもよい。
【０１１３】
「非中継座標」は、質点系の構造を記述しかつ中継座標から独立した座標である。非中継座標も、中継座標と同様に、質点系に含まれる一部の質点の座標そのもの、質点の座標を組み合わせて定義できる座標又はこれらの組み合わせを用いて設定される。質点座標の組み合わせによっては、非中継座標の構成の一部を中継座標で表現することもできる。ただし、非中継座標は、中継座標から独立した座標であるため、中継座標のみによって表現されることはない。例えば非中継座標は、一部の質点の座標が中継座標として設定されている場合には、Ｎ個の質点のうちその質点を除く質点座標とすることができる。また、非中継座標は、質点座標を組み合わせて定義できる座標のみが中継座標として設定されている場合には、そのすべてまたは一部を、質点座標そのもの或いは質点座標を組み合わせて定義できる座標で、かつ中継座標とは独立な座標とすることができる。
【０１１４】
中継座標の具体的な設定方法の２つの例について説明する。
【０１１５】
１つは、多原子系に含まれるＮ個の原子を特徴的な部分構造ごとに集団分けして原子の集団ごとに重心を抽出し、これらの重心の座標のみを中継座標として設定する方法である。この場合には、非中継座標は、その原子が属する部分構造の重心座標に対する相対座標でかつ独立な座標とすることが好ましい。したがって、非中継座標の個数は３Ｎ−Ｍ個となる。これは、多原子系の３Ｎ個の質点座標をＭ個の重心座標と３Ｎ個の相対座標に分離したとき、相対座標間には定義からＭ個の関係式が発生するため、つまり相対座標のうちＭ個は独立な座標ではないためである。この方法は、特徴的な部分構造同士の位置関係が大振幅集団運動による多原子系の構造変化に大きく影響することを考慮したものである。特徴的な部分構造自体の変形の規模が小さく、かつその変形の時間スケールが化学反応の時間スケールに比べて短いため、その変形が多原子系の構造変化に与える影響は小さい。
【０１１６】
もう１つは、生体高分子の主鎖を構成する原子のそれぞれの座標を中継座標として設定する方法である。この場合には、非中継座標は、その生体高分子の側鎖および溶媒和等のその他の原子の座標とすることが好ましい。したがって、非中継座標の個数は３Ｎ−Ｍ個となる。この方法は、生体高分子の主鎖の形状が大振幅集団運動による多原子系の構造変化に大きく影響することを考慮したものである。側鎖は主鎖に従属して移動することから、側鎖の変形が多原子系の構造変化に与える影響は小さい。この場合において、結合候補分子等の生体高分子以外の分子が多原子系に含まれている場合には、生体高分子以外の分子については、その分子の重心を中継座標として設定してもよい。この場合には、例えばその分子を構成する原子の座標は非中継座標として設定される。
【０１１７】
なお、後述する第１の実施形態は前者の設定方法を用いた場合に相当する。
【０１１８】
中継座標は、必要に応じて、すべての質点を自由な状態（質点系に固定しない状態）にした上で質点系に対して全構造緩和が行われた後に、設定されることが好ましい。「全構造緩和」とは、適当な開始状態から質点系のポテンシャルエネルギーが常に減少していくように、そのポテンシャルエネルギーの勾配がゼロとなる状態まで、質点の座標を移動させることを意味する。全構造緩和はいわゆる一般的な構造緩和と同義であるが、後述する「部分構造緩和」と区別するためにこのように称することにする。この全構造緩和によって見出される質点系の状態は、上記開始状態近傍のローカルミニマムの一つである。つまり、この全構造緩和は、ポテンシャル曲面上の上下を繰り返しながら安定点間を結ぶ反応経路およびミニマムポイント（最小点）を見出す作業ではない。この全構造緩和は、特に、構造緩和する対象を構成する原子のすべてを移動させる点で後述する「部分構造緩和」と区別される。この全構造緩和は、例えば、共役勾配法、最急降下法、逆ヘッセ法等既存の方法を用いて実施される。
【０１１９】
＜座標抽出手段＞
座標抽出手段は、中継座標の関数としての非中継座標の構成を求め、集団座標の関数としての中継座標の構成を求める。「中継座標の関数としての非中継座標の構成」とは、中継座標ＲＳで表された非中継座標ＲＦについての関数ＲＦ（ＲＳ）の具体的な内容を意味し、「集団座標の関数としての中継座標の構成」とは、集団座標ｑ†で表された中継座標ＲＳについての関数ＲＳ（ｑ†）の具体的な内容を意味する。
【０１２０】
（中継座標の関数としての非中継座標の構成）
非中継座標の構成ＲＦ（ＲＳ）は、非中継座標ＲＦを中継座標ＲＳに従属させること、つまり非中継座標ＲＦおよび中継座標ＲＳの関係を規定する条件式を与えることにより得られる。このような条件式は、特に限定されないが、例えば断熱近似条件から得られる下記式３０の条件式を用いることができる。なお、「断熱近似」とは、非中継座標ＲＦが中継座標ＲＳの変化に対し即座に追随することができるとする近似を意味する（M. Born, and J. R. Oppenheimer, Ann. Phys. 84 (1927) 457およびH. Haken, Synergetics 2rd Edition, Springer, 1978参照）。式３０は、与えられたＲＳに対し、ＲＦが、ポテンシャルエネルギーＶの局所安定点、つまりポテンシャルエネルギーＶのＲＦに関する勾配がゼロの点にいるという条件を表している。
【０１２１】
【数１６】
0005697638.tif 000017
【０１２２】
式３０において、Ｖ（ＲＳ，ＲＦ）は、中継座標ＲＳおよび非中継座標ＲＦを用いて表された質点系のポテンシャルエネルギーを表す。Ｖ（ＲＳ，ＲＦ）は、質点座標ｘで表されたポテンシャルエネルギーＶ（ｘ）に、中継座標ＲＳが設定された段階で定まる質点座標ｘの構成ｘ（ＲＳ，ＲＦ）を代入することにより得ることができる。例えば、質点座標ｘの構成ｘ（ＲＳ，ＲＦ）は、中継座標ＲＳが設定された段階で定まる中継座標ＲＳの構成ＲＳ（ｘ）および非中継座標ＲＦの構成ＲＦ（ｘ）に基づいて、これらの逆関数を求めることにより得られる。
【０１２３】
なお、すべての非中継座標ＲＦ１〜ＲＦ（３Ｎ−Ｍ）を中継座標ＲＳに従属させるため、すべての非中継座標ＲＦ１〜ＲＦ（３Ｎ−Ｍ）でＶ（ＲＳ，ＲＦ）を偏微分することが必要である。
【０１２４】
非中継座標ＲＦが中継座標ＲＳに従属した質点系においては、下記式３１のように非中継座標ＲＦが中継座標ＲＳとして表されることになり、質点系のポテンシャルエネルギーＶは下記式３２のように表される。以下、上記のような非中継座標ＲＦが中継座標ＲＳに従属した条件の下、中継座標ＲＳのみの関数として表された質点系のポテンシャルエネルギーＶｅｆｆを有効ポテンシャルエネルギーともいう。
【０１２５】
【数１７】
0005697638.tif 000018
【０１２６】
また、式３１のように非中継座標ＲＦが中継座標ＲＳに従属する場合には、質点座標ｘは下記式３３のように中継座標の関数として表すことができる。
【０１２７】
【数１８】
0005697638.tif 000019
【０１２８】
（集団座標の関数としての中継座標の構成）
中継座標の構成ＲＳ（ｑ†）は、非中継座標ＲＦが中継座標ＲＳに従属するという条件の下、中継座標ＲＳの変化に伴う非中継座標ＲＦの変化が中継座標ＲＳに与える影響を考慮しながら、中継座標ＲＳに対して集団運動理論における基礎方程式を解くことにより求められる。中継座標の構成ＲＳ（ｑ†）を求めるに際し重要な点は、「中継座標ＲＳに対して集団運動理論における基礎方程式を解くこと」およびその際に「中継座標ＲＳの変化に伴う非中継座標ＲＦの変化が中継座標ＲＳに与える影響（フィードバック）を考慮すること」である。
【０１２９】
以下、この２点について説明する。
【０１３０】
中継座標ＲＳに対して集団運動理論における基礎方程式を解く理由は、前述したように、３Ｎ次元の質点座標ｘから集団座標ｑ†を抽出するのではなく、Ｍ次元の中継座標ＲＳから集団座標ｑ†を抽出することにより、変数分離の計算を単純化するためである。
【０１３１】
ＳＣＣ法における最も基本的な基礎方程式は下記式３４から式３６である（非特許文献１０）。
【０１３２】
【数１９】
0005697638.tif 000020
【０１３３】
ｍｉは質点系におけるｉ番目の中継座標に関する質量を表す。ｍｉは、ｉ番目の中継座標が複数の質点の重心座標である場合には、それらの質点の質量の総量となる。中継座標が、質点座標を組み合わせた座標であって重心座標ではない場合、単純に質量の総量にはならない。この場合、通常の解析力学の処方で中継座標に対応する一般運動量を求め、ハミルトニアンを書き下し、そのハミルトニアンの一般運動量が関与する項の係数として、中継座標に関する質量を求める。また、φｉμ（ＲＳ）、λμ（ＲＳ）およびρμ（ＲＳ）は式３５および式３６を満たすＲＳの関数であり、φｉμ（ＲＳ）は（ｉ，μ）番目、λμ（ＲＳ）およびρμ（ＲＳ）はμ番目の成分である。また、Ｖ，ｉ（ＲＳ）およびＶ,ｉｊ（ＲＳ）はそれぞれ下記式３７および式３８で定義される。
【０１３４】
【数２０】
0005697638.tif 000021
【０１３５】
したがって、一般的には式３４から式３６に従い、Ｍ個の座標（ＲＳ１，ＲＳ２，・・・，ＲＳＭ）で張られるＭ次元の超空間の中にＫ個の集団座標（ｑ１，ｑ２，・・・，ｑＫ）でパラメトライズ（媒介変数化）されるＫ次元の超平面（以下、単に平面という）を構成することで、関数ＲＳｉ（ｑμ）が求められる。
【０１３６】
また、基礎方程式としては非特許文献１０における近似計算のための基礎方程式を援用することもできる。このとき、例えば自由度が１（Ｋ＝１）の集団座標ｑ†の関数としての中継座標の構成ＲＳ（ｑ１）を求める場合には、中継座標ＲＳを用いた集団運動理論における基礎方程式は下記式３９および式４０のように表される。
【０１３７】
【数２１】
0005697638.tif 000022
【０１３８】
φｉ（ＲＳ）は、式４０を満たす関数（固有ベクトル）のｉ番目の成分を表し、
Λ（ＲＳ）は、式４０を満たす関数（固有値）を表す。
【０１３９】
一方、従来の非特許文献１０におけるＳＣＣ法では、近似的に計算が行われており、変数分離を高精度に行うことができない場合がある。そこで、本発明者は、中継座標ＲＳに対する変数分離をより正確に計算することができる新たな近似法に基づき基礎方程式を導出し、この基礎方程式に基づく新たなＳＣＣ法（ＳＣＣ２法）の理論を構築した。
【０１４０】
例えばＳＣＣ２法に基づいて、自由度が１（Ｋ＝１）の集団座標ｑ†の関数としての中継座標の構成ＲＳ（ｑ１）を求める場合には、中継座標ＲＳを用いた集団運動理論における基礎方程式は下記式４１から式４３のように表される。
【０１４１】
【数２２】
0005697638.tif 000023
【０１４２】
φｉ（ＲＳ，λ）およびκ（ＲＳ，λ）は式４３を満たす関数であり、φｉ（ＲＳ，λ）はそのｉ番目の成分を表し、
λは、補助座標を表す。
【０１４３】
最終的に、式３９または式４１の解として、式４４で表されるような中継座標ＲＳの構成が求められる。
【０１４４】
【数２３】
0005697638.tif 000024
【０１４５】
また、上記ではＫ＝１の場合における基礎方程式について説明したが、本発明はＫ＝１の場合に限られない。具体的には、Ｋ≧２の場合におけるＳＣＣ法の基礎方程式への式３９および式４０の拡張は、非特許文献１０に記載されているように、フロベニウスの定理を用いて行うことができる。一方、Ｋ≧２の場合におけるＳＣＣ２法の基礎方程式への式４１から式４３の拡張も、同様の方法で自明に行うことができる。
【０１４６】
このように、中継座標という階層化された変数（つまり質点系の構造を二次的に記述するより低次元の変数）を導入することにより、３Ｎ次元からＫ次元への縮約の問題がＭ次元からＫ次元への縮約の問題に簡略化されるため、変数分離の計算が容易となる。
【０１４７】
しかしながら、本発明では、質点系の構造を正確に計算するために、フィードバックを考慮しながら変数分離を実施しなければならない。フィードバックを考慮しないということは、存在する質点（または非中継座標）を存在しないものとみなして計算することと同義であり、この影響を考慮せずに変数分離を行うと計算内に矛盾が生じて計算結果が破錠する問題が起こるためである。そこで、フィードバックを考慮するために、既存の集団運動理論に新たな理論に基づいた修正を加える。この新たな理論とは、非中継座標ＲＦが中継座標ＲＳに従属するという条件の下、集団運動理論の基礎方程式におけるポテンシャルエネルギーの座標に関する微分係数を正確に計算し、さらに、中継座標ＲＳ（ｑ†）を微小更新（ＲＳ（ｑ†）→ＲＳ（ｑ†＋Δｑ†））させるごとに非中継座標ＲＦに対して部分構造緩和を行うというものである。
【０１４８】
従来の集団運動理論では、ポテンシャルエネルギーＶ（ｘ）の質点座標ｘに関する微分係数に対する所定の固有ベクトルφの方向に、質点座標ｘを微小更新（ｘ（ｑ†）→ｘ（ｑ†＋Δｑ†））することが繰り返されるのみである。
【０１４９】
一方本発明では、ポテンシャルエネルギーＶの座標に関する微分係数について、その計算をする際、下記式４５から式４７のフィードバック方程式を適用することができる。
【０１５０】
【数２４】
0005697638.tif 000025
【０１５１】
式４６において、第３項の（ｉ←→ｊ）は、第２項におけるｉおよびｊの文字を相互に入れ替えた項を表す。
【０１５２】
式４７において、第３項の（ｉ←→ｋ）は、第２項におけるｉおよびｋの文字を相互に入れ替えた項を表し、第４項の（ｊ←→ｋ）は、第２項におけるｊおよびｋの文字を相互に入れ替えた項を表し、第６項の（ｉ←→ｋ）は、第５項におけるｉおよびｋの文字を相互に入れ替えた項を表し、第７項の（ｊ←→ｋ）は、第５項におけるｊおよびｋの文字を相互に入れ替えた項を表す。
【０１５３】
さらに、式４５から式４７において下記式４８を用いた。
【０１５４】
【数２５】
0005697638.tif 000026
【０１５５】
Ｋαβ−１（Ｒｓ）はＫαβ（Ｒｓ）の逆行列を表し、
Ｋαβ（Ｒｓ）およびＪαｉ（Ｒｓ）はそれぞれ式４９および式５０で定義される。
【０１５６】
【数２６】
0005697638.tif 000027
【０１５７】
式４５から式４７それぞれにおいて、右辺の第１項以外の項が、中継座標ＲＳの変化に伴う非中継座標ＲＦの変化が中継座標ＲＳに与える影響を表している。このようなフィードバック方程式の構成は従来知られておらず、本発明者が初めて見出したものである。フィードバック方程式は、１階微分、２階微分および３階微分の３つの式を含むが、集団運動理論における基礎方程式の内容に応じて必要なものを使用すればよい。例えば、集団運動理論における基礎方程式として、式３９および式４０を採用した場合には、フィードバック方程式は式４６の２階微分の式のみで充分である。
【０１５８】
そして、非中継座標ＲＦに対する部分構造緩和は、ポテンシャルエネルギーＶの座標に関する微分係数であってフィードバック方程式で規定されるものに対する固有ベクトルφの方向に、中継座標ＲＳ（ｑ†）を微小更新（ＲＳ（ｑ†）→ＲＳ（ｑ†＋Δｑ†））させ、更新された中継座標ＲＳ（ｑ†＋Δｑ†）の下で、非中継座標ＲＦが中継座標ＲＳに従属するという条件に従って構造緩和を行うことにより行われる。「中継座標ＲＳ（ｑ†＋Δｑ†）の下で」とは、更新された中継座標ＲＳ（ｑ†＋Δｑ†）を質点系に固定した状態でという意味である。また、「非中継座標ＲＦが中継座標ＲＳに従属するという条件に従って構造緩和を行う」とは、質点系のポテンシャルエネルギーが常に減少していくように、そのポテンシャルエネルギーの勾配がゼロとなる状態まで、中継座標ＲＳに従属した状態を維持しながら非中継座標ＲＦを移動させることを意味する。上記のように中継座標ＲＳを質点系に固定した状態で行う構造緩和を「部分構造緩和」という。前述したように、この部分構造緩和は、中継座標ＲＳを質点系に固定した状態で構造緩和を行う点で全構造緩和と異なる。ただし、部分構造緩和においても、構造緩和を行う手法としては、全構造緩和と同様に、例えば、共役勾配法、最急降下法、逆ヘッセ法等既存の方法を用いることができる。なお、更新する方向を決める上記固有ベクトルφは、最小固有値に対応するものであることが好ましい。本発明は、質点系の大振幅集団運動を対象としているためである。
【０１５９】
Ｋ＝１である場合、部分構造緩和における非中継座標ＲＦの開始状態は、ＲＦα（ＲＳ（ｑ１））でもよいが、下記式５１で与えられる微小量だけ変化させた状態、つまり下記式５２で表される状態ＲＦα†であることが好ましい。
【０１６０】
【数２７】
0005697638.tif 000028
【０１６１】
式５２で表される座標ＲＦα†は、微小更新が無限小ならば非中継座標ＲＦ（ＲＳ（ｑ１＋Δｑ１））と一致する。しかしながら、実際の計算では微小更新は有限である。部分構造緩和における非中継座標ＲＦの開始状態をＲＦα†とすることにより、部分構造緩和を行う上での安全性を高めることができる。
【０１６２】
（座標抽出手段における処理手順）
座標抽出手段における一連の具体的な処理手順について説明する。
【０１６３】
まず、座標抽出手段は、中継座標Ｒｓおよび非中継座標ＲＦの関数として表現された質点系のポテンシャルエネルギーＶ（ＲＳ，ＲＦ）を求める第１の工程を実施する。
【０１６４】
次に、座標抽出手段は、非中継座標ＲＦを中継座標ＲＳに従属させる第２の工程を実施する。中継座標ＲＳの関数としての非中継座標ＲＦの構成が得られる。ここでは、例として式３０の断熱近似条件に従って非中継座標ＲＦを中継座標ＲＳに従属させた場合を考える。
【０１６５】
次に、座標抽出手段は、座標設定手段によって設定された状態における中継座標ＲＳおよび非中継座標ＲＦの下で、例えば式４６のフィードバック方程式に従って、中継座標ＲＳに関するポテンシャルエネルギーＶの微分係数を求める第３の工程を実施する。
【０１６６】
次に、座標抽出手段は、ポテンシャルエネルギーＶの上記微分係数の下で、例えば式３９および式４０のＳＣＣ法におけるＫ＝１の場合の基礎方程式に従って、集団座標ｑに関する中継座標ＲＳ（ｑ）の微分係数を求める第４の工程を実施する。
【０１６７】
次に、座標抽出手段は、中継座標ＲＳ（ｑ）の上記微分係数の下で集団座標ｑを微小量Δｑだけ更新することにより、更新された中継座標ＲＳ（ｑ＋Δｑ）を求める第５の工程を実施する。
【０１６８】
次に、座標抽出手段は、更新された中継座標ＲＳ（ｑ＋Δｑ）の下で、上記断熱近似条件に従って中継座標ＲＳに従属した非中継座標ＲＦについて部分構造緩和を行う第６の工程を実施する。この際、式５２で表される状態ＲＦα†を部分構造緩和の開始状態とする場合には、第４の工程で得られた中継座標ＲＳ（ｑ）の上記微分係数から式５１に従って、集団座標ｑに関する非中継座標ＲＦ（ｑ）の微分係数を求め、この微分係数の下で非中継座標を微小量Δｑだけ更新させておく。
【０１６９】
そして、座標抽出手段は、中継座標ＲＳおよび非中継座標ＲＦを交互に更新しながら、第３の工程から第６の工程までを例えば次の極小点に至るまで順次繰り返すことにより、集団座標ｑの関数としての中継座標ＲＳの構成を求める。つまり、２回目以降の第３の工程においては、直前の第６の工程で行われた非中継座標ＲＦの部分構造緩和後の状態における中継座標ＲＳおよび非中継座標ＲＦの下で、上記第３の工程と同様に、中継座標ＲＳに関するポテンシャルエネルギーＶの微分係数を求めることになる。極小点を判定する方法は、特に限定されないが、例えばポテンシャルエネルギーＶの上記微分係数がゼロ相当になったか否かを判定する方法を挙げることができる。
【０１７０】
＜逆変換手段＞
逆変換手段は、集団座標ｑ†についての運動方程式の解として得られる時間ｔの関数としての集団座標ｑ†（ｔ）、式４４で表される中継座標ＲＳの上記構成ＲＳ（ｑ）および式３１で表される非中継座標ＲＦの上記構成ＲＦ（ＲＳ）に基づいて、質点座標ｘの時間発展を予測する。
【０１７１】
具体的には以下の通りである。図４は、本発明における変数変換および逆変換の概念を示す図である。図４では、例として質点座標ｘの一部の座標そのものを中継座標ＲＳとし、この中継座標ＲＳから自由度が１の集団座標ｑ†＝ｑ１を抽出した場合について示している。
【０１７２】
図４ａは、質点座標ｘに基づいて中継座標ＲＳおよび非中継座標ＲＦを設定する様子を示している。図４ａでは、Ｍ個の質点座標が中継座標ＲＳとして設定され、それ以外の３Ｎ−Ｍ個の質点座標が非中継座標ＲＦとして設定されている。中継座標ＲＳおよび非中継座標ＲＦの設定は、前述したとおり、座標設定手段において行われる。図４ｂは、非中継座標ＲＦが中継座標ＲＳの関数として与えられることを示している。図４ｃは、中継座標ＲＳが集団座標ｑ１の関数として与えられることを示している。非中継座標ＲＦの構成および中継座標ＲＳの構成は、前述したとおり、座標抽出手段において求められる。
【０１７３】
図４ｄは、集団座標ｑ１が時間ｔの関数として与えられることを示している。集団座標ｑ１についての運動方程式を解くことにより求められる。集団座標ｑ１についての運動方程式は、質点座標ｘについての運動方程式に、中継座標ＲＳが設定された段階で定まる質点座標ｘの構成ｘ（ＲＳ，ＲＦ）を代入して、中継座標ＲＳおよび非中継座標ＲＦについての運動方程式を求め、この運動方程式に式３１および式４４を代入することにより求めることができる。
【０１７４】
以上の図４ａからｄまでの関係を用いることにより、すべての質点座標ｘを時間ｔの関数として求めることができる（図４ｅ）。つまり、逆変換手段は、時間ｔの関数としての集団座標ｑ１を求めた後、中継座標ＲＳを介した質点座標ｘから集団座標ｑ１への変数変換に対してこれを逆変換することにより、時間ｔの関数としての質点座標ｘ（ｔ）の構成を求めるものである。この関数ｘ（ｔ）は、時間変化に対する質点座標の挙動を表すものであるから、この関数ｘ（ｔ）こそが求めるべき反応経路に相当するものである。したがって、この関数ｘ（ｔ）の時間変化を見ることにより質点座標の時間発展を予測することが可能となる。
【０１７５】
「本発明の第１の実施形態」
図５は、本実施形態のシミュレーション装置１０の構成を概略的に示すブロック図である。また、図６Ａおよび図６Ｂは、本実施形態のシミュレーション方法における計算手順を示すフローチャート図である。なお、本実施形態では、生体高分子およびこれに対する結合候補分子を含むＮ個の原子から構成される多原子系において、多原子系の構造のうち特徴的な部分構造ごとに抽出された重心の座標をＭ個の中継座標として設定し、かつ、集団座標の自由度が１（Ｋ＝１）である場合について説明する。
【０１７６】
本実施形態のシミュレーション装置１０は、生体高分子および結合候補分子から構成される複合体の挙動を予測するシミュレーション装置であって、図５に示されるように、解析対象となる多原子系を規定するのに必要なデータ（入力データ）を入力するための入力手段１２と、装置の各部を制御する制御手段１４と、座標設定手段１６と、座標抽出手段１８と、逆変換手段２０と、解析結果を表示する表示手段２２とを備える。
【０１７７】
そして、本実施形態のシミュレーション方法は、シミュレーション装置１０により実行されるシミュレーション方法であって、入力手段１２によって入力データを入力し、座標設定手段１６によってＭ個の中継座標ＲＳおよび３Ｎ−Ｍ個の非中継座標ＲＦを設定させ、座標抽出手段１８によって、Ｋ＝１の集団座標ｑ１の関数としての中継座標ＲＳの構成および中継座標ＲＳの関数としての非中継座標ＲＦの構成を求めさせ、逆変換手段２０によって時間ｔの関数としての原子の座標ｘ（ｔ）を求めさせ、表示手段２２によって解析結果を表示させるものである。
【０１７８】
また、本実施形態に係るシミュレーションプログラムは、上記に記載のシミュレーション方法をコンピュータに実行させることを特徴とするものである。
【０１７９】
また、本実施形態に係るコンピュータ読み取り可能な記録媒体は、上記に記載のシミュレーションプログラムを記録したものである。
【０１８０】
＜入力手段＞
入力手段１２は、ユーザーにより入力データが入力される部分である。入力された入力データは制御手段１４へ出力される。入力データの入力方法は、特に限定されず、例えばシミュレーション装置１０をユーザーが操作しながら手動で行う方法、所定の記録媒体から読み込む方法等が挙げられる。
【０１８１】
＜制御手段＞
制御手段１４は、各部とのデータのやり取り等の処理を制御する部分である。入力手段１２から制御手段１４へ入力された入力データは座標設定手段１６に出力される。入力データが座標設定手段１６へ出力されることにより、当該入力データに基づく多原子系の挙動の計算が開始される。制御手段１４には、各部とのデータのやり取りの中で生成される途中のおよび最終的な計算結果を記録するため、メモリ等の記憶媒体が含まれる。
【０１８２】
また、制御手段１４は、逆変換手段２０で得られた原子の座標ｘ（ｔ）に関するデータに基づいて、その結果を図やグラフを用いて表示するよう表示手段２２を制御する。
【０１８３】
＜座標設定手段＞
座標設定手段１６は、制御手段１４から座標設定手段１６へ入力された入力データに基づいて、複合体の構造を最適化し、その後Ｍ個の中継座標ＲＳおよび３Ｎ−Ｍ個の非中継座標ＲＦを設定するものである。得られた中継座標ＲＳおよび非中継座標ＲＦに関するデータは、制御手段１４を介して座標抽出手段１８に出力される。
【０１８４】
＜座標抽出手段＞
座標抽出手段１８は、座標設定手段１６で得られた中継座標ＲＳおよび非中継座標ＲＦに関するデータに基づいて、多原子系全体のポテンシャルエネルギーＶ（ＲＳ，ＲＦ）を求め、このＶ（ＲＳ，ＲＦ）に基づいて集団座標ｑ１の関数としての中継座標ＲＳの構成および中継座標ＲＳの関数としての非中継座標ＲＦの構成を求めるものである。得られた中継座標ＲＳの構成および非中継座標ＲＦの構成に関するデータは、制御手段１４を介して逆変換手段２０に出力される。
【０１８５】
＜逆変換手段＞
逆変換手段２０は、座標抽出手段１８で得られた中継座標ＲＳの構成および非中継座標ＲＦの構成に関するデータに基づいて、時間ｔの関数としての原子の座標ｘ（ｔ）を求めるものである。得られた原子の座標ｘ（ｔ）に関するデータは制御手段１４に出力される。
【０１８６】
＜表示手段＞
表示手段２２は、制御手段１４の指示に従い、例えば多原子系のポテンシャルエネルギーＶをグラフで表示したり、複合体の結合過程における構造および／または最終状態の構造を画像として表示したり、複合体の結合過程を動画として表示したりするものである。表示方法としては特に限定されず、２Ｄ表示、３Ｄ表示等公知の方法を用いることが可能である。
【０１８７】
以下、図６Ａおよび図６Ｂを参照しながら、本実施形態の装置１０を用いたシミュレーション方法の手順を詳細に説明する。
【０１８８】
＜ＳＴＥＰ１＞
まず、入力手段１２から入力データがシミュレーション装置１０に入力される（図６ＡのＳＴ１）。入力データの具体的な内容は、生体高分子およびこれに対する結合候補分子を含むＮ個の原子から構成される多原子系を構成する原子の種類、個数、初期配置の情報となる座標ｘおよび座標ｘからポテンシャルエネルギーＶ（ｘ）を構成するルール等である。
【０１８９】
＜ＳＴＥＰ２＞
次に、複合体を含む多原子系の構造が最適化されて、初期状態における構造が求められる（図６ＡのＳＴ２）。
【０１９０】
具体的には以下の通りである。まず、座標設定手段１６は、生体高分子および結合候補分子に対してそれぞれ独立に全構造緩和を実施して、生体高分子および結合候補分子の構造をそれぞれ最適化する。そして、それぞれ構造が最適化された生体高分子および結合候補分子は、互いに相互作用を開始する程度（例えば、結合候補分子が生体高分子の表面にわずかに接触する程度）の距離に配置される。この状態が、生体高分子および結合候補分子からなる複合体の初期状態となる。
【０１９１】
＜ＳＴＥＰ３＞
次に、ＳＴＥＰ２で得られた初期状態における構造に基づいて、中継座標ＲＳおよび非中継座標ＲＦが設定される（図６ＡのＳＴ３）。
【０１９２】
座標設定手段１６は、生体高分子の特徴的な部分構造ごとに（Ｍ−３）／３個の部分構造に分け、結合候補分子全体を１個の部分構造として取り扱い、Ｍ／３個の部分構造のそれぞれから重心を抽出する。そして、このＭ／３個の重心の座標がＭ個の中継座標ＲＳとして設定され、Ｍ／３個の各部分構造に含まれる原子の各重心からの相対座標でかつ独立な座標が３Ｎ−Ｍ個の非中継座標ＲＦとして設定される。本実施形態では、非中継座標ＲＦとして、その原子の属する部分構造の重心座標からの相対座標が設定されるとする。
【０１９３】
＜ＳＴＥＰ４＞
次に、中継座標ＲＳおよび非中継座標ＲＦで表された多原子系のポテンシャルエネルギーＶ（ＲＳ，ＲＦ）が求められる（図６ＡのＳＴ４）。
【０１９４】
図７は、本実施形態における中継座標ＲＳ、非中継座標ＲＦおよび原子の座標ｘの関係を示す概略図である。なお、図７において、ｍ（ｍは１≦ｍ≦Ｍ／３を満たす整数）番目の特徴的な部分構造に含まれる原子数をＮｍと表した。つまり、Ｎ１＋Ｎ２＋・・・＋Ｎｍ＋・・・＋ＮＭ／３＝Ｎである。本発明では、Ｍは３次元空間における原子（質点）の座標の個数であるため、通常Ｍ／３は整数となる。また、本実施形態においては、中継座標ＲＳ、非中継座標ＲＦおよび原子の座標ｘは、図７から具体的に下記式５３および式５４で表されるような関係を有していることがわかる。
【０１９５】
【数２８】
0005697638.tif 000029
【０１９６】
式５３は、ｍ番目の中継座標ＲＳ、つまりｍ番目の特徴的な部分構造の重心座標と原子の座標ｘとの関係を示した式である。式５４は、ｍ番目の特徴的な部分構造に属するｎ番目の原子の座標ｘ、中継座標ＲＳおよび非中継座標ＲＦの関係を示した式である。式５３および式５４において、ｎは、特徴的な部分構造のそれぞれに含まれる原子の通し番号である。したがって、ｍ番目の特徴的な部分構造についてのｎは、１からＮｍまでの整数値をとる。また、３Ｎ１＋３Ｎ２＋・・・＋３Ｎｍ−１はｍ＝１のときゼロとする。
【０１９７】
したがって、ポテンシャルエネルギーＶ（ＲＳ，ＲＦ）は、ポテンシャルエネルギーＶ（ｘ）に対して式５４を適用することにより得ることができる。
【０１９８】
＜ＳＴＥＰ５＞
次に、ポテンシャルエネルギーＶ（ＲＳ，ＲＦ）を用いた式３０の断熱近似条件に従って、非中継座標ＲＦが中継座標ＲＳに従属させられる（図６ＡのＳＴ５）。本実施形態では、すべての非中継座標ＲＦ１〜ＲＦ３ＮでＶ（ＲＳ，ＲＦ）を偏微分することになる。
【０１９９】
＜ＳＴＥＰ６＞
次に、現状の中継座標ＲＳおよび非中継座標ＲＦの下で、式４６のフィードバック方程式に従い、中継座標ＲＳに関するポテンシャルエネルギーＶの２階微分係数が求められる（図６ＡのＳＴ６）。
【０２００】
＜ＳＴＥＰ７＞
次に、ＳＴＥＰ６で得られた２階微分係数の下で、式３９および式４０のＳＣＣ法の基礎方程式に従い、集団座標ｑ１に関する中継座標ＲＳの微分係数ｄＲＳ／ｄｑ１が求められる（図６ＡのＳＴ７）。
【０２０１】
＜ＳＴＥＰ８＞
次に、ＳＴＥＰ７で得られた微分係数ｄＲＳ／ｄｑ１の下で、式５１に従い集団座標ｑ１に関する非中継座標ＲＦの微分係数ｄＲＦ／ｄｑ１が求められる（図６ＢのＳＴ８）。
【０２０２】
＜ＳＴＥＰ９＞
次に、ＳＴＥＰ７で得られた微分係数ｄＲＳ／ｄｑ１およびＳＴＥＰ８で得られた微分係数ｄＲＦ／ｄｑ１の下で、集団座標ｑ１が微小量Δｑ１だけ更新される（図６ＢのＳＴ９）。
【０２０３】
＜ＳＴＥＰ１０＞
次に、微小量Δｑ１だけ更新された中継座標ＲＳ（ｑ１＋Δｑ１）および更新された非中継座標ＲＦ（ｑ１＋Δｑ１）が求められる（図６ＢのＳＴ１０）。
【０２０４】
例えば、ＳＴＥＰ２における初期状態を集団座標ｑ１の初期条件と設定した（つまり、ｑ１＝０のときに多原子系が上記初期状態をとるものとした）場合には、中継座標ＲＳ（ｑ１＋Δｑ１）および非中継座標ＲＦ（ｑ１＋Δｑ１）は下記式５５のように表せられる。
【０２０５】
【数２９】
0005697638.tif 000030
【０２０６】
また、ＳＴＥＰ１３の判断結果に応じてＳＴＥＰ６からＳＴＥＰ１２までを順次繰り返し実施する場合において、２回目以降のＳＴＥＰ１０で順次下記式５６のように表せられるＲＳｉ（２Δｑ１）、ＲＳｉ（３Δｑ１）、・・・が得られることになる。ただし、微分係数ｄＲＳ／ｄｑ１は各回で同じ値になるとは限らない。上記の作業を繰り返すことにより、最終的に集団座標ｑ１の関数としての中継座標ＲＳの構成ＲＳ（ｑ１）が求められる。
【０２０７】
【数３０】
0005697638.tif 000031
【０２０８】
＜ＳＴＥＰ１１＞
ＳＴＥＰ１１では、ＳＴＥＰ１０で更新された中継座標ＲＳ（ｑ１＋Δｑ１）の下で、ＳＴＥＰ１０で更新された非中継座標ＲＦ（ｑ１＋Δｑ１）を開始状態として非中継座標ＲＦについて部分構造緩和が行われる（図６ＢのＳＴ１１）。
【０２０９】
＜ＳＴＥＰ１２＞
次にＳＴＥＰ１２において、部分構造緩和後の状態における多原子系のポテンシャルエネルギーＶが計算される（図６ＢのＳＴ１２）。なお、ＳＴＥＰ１３の判断結果に応じて、繰り返しの各回でポテンシャルエネルギーＶが計算される。ポテンシャルエネルギーの数値計算は、有効ポテンシャルエネルギーＶｅｆｆに基づいた計算でもよい。
【０２１０】
＜ＳＴＥＰ１３＞
次にＳＴＥＰ１３において、ＳＴＥＰ１２で計算されたポテンシャルエネルギーＶが極小点に至ったか否かが判断される（図６ＢのＳＴ１３）。生体高分子および結合候補分子の化学反応を考えた場合、化学反応の終了点はポテンシャルエネルギーＶが次に極小となるとき、つまりポテンシャル曲面上の次の安定点である。そこで、ポテンシャルエネルギーＶの値を基準に化学反応の終了点が判断される。ポテンシャルエネルギーＶが次の極小点に至ったと判断されればＳＴＥＰ１４へ進み、そうでなければ再度ＳＴＥＰ６からＳＴＥＰ１２までが繰り返される。また、ポテンシャルエネルギーＶが極小点に至ったと判断された場合であっても、必要に応じて、ポテンシャルエネルギーＶがさらに次の極小点に至るまで解析を進めることも可能である。
【０２１１】
＜ＳＴＥＰ１４＞
ＳＴＥＰ１３においてポテンシャルエネルギーＶが極小点に至ったと判断された場合、ＳＴＥＰ１４において、そのポテンシャルエネルギーＶにおける多原子系の構造が、求めるべき最終状態の構造となる（図６ＢのＳＴ１４）。
【０２１２】
そして、最終状態の構造が得られるとともに、式４４で表される中継座標ＲＳの構成ＲＳ（ｑ１）および式３１で表される非中継座標ＲＦの構成ＲＦ（ＲＳ）が求められる。その後、時間ｔの関数としての座標ｘ（ｔ）が求められる。
【０２１３】
＜作用効果＞
図８は、本実施形態のシミュレーション方法によって得られた、所定のタンパク質２およびこれに対する生理活性物質６（薬物）の複合体の反応経路を示すグラフである。また、図９は、図８のグラフにおける複合体の結合過程を示す概略図である。図９ａは複合体の結合過程における初期状態を示し、図９ｂは複合体の結合過程における最終状態を示す。
【０２１４】
図８のグラフにおいて、横軸は所定のタンパク質２の残基２ａと残基２ｂの距離であり、縦軸は薬物６および上記タンパク質２の残基２ｃの距離である。残基２ａおよび残基２ｂは、タンパク質２のポケット４の開口にある残基であり、残基２ｃは、タンパク質２および薬物６が結合した時に薬物６に対向する位置にある残基である（図９）。また、グラフ中のＡおよびＣが、例えば図２における安定点ＡおよびＣにそれぞれ対応し、グラフ中のＢが、例えば図２におけるエネルギー障壁Ｂに対応している。グラフ中のプロットは時間軸上で等間隔な点を示しており、このグラフからエネルギー障壁Ｂに至るまでは反応がゆっくり進み、エネルギー障壁Ｂを超えた後は反応が速やかに進む様子がわかる。本実施形態のシミュレーション方法により、２０００個程度の原子を含む多原子系においてこの反応経路を求めるのに要した時間はおよそ３０分である。
【０２１５】
また、図８のグラフ中の「Ｘ線解析」の位置は、上記タンパク質２および上記薬物６についてＸ線回折測定による結晶構造解析の結果実際に得られた実測値を示す。これにより、本発明は非経験的にシミュレーションを行っているにもかかわらず、非常に高い精度でこの複合体の結合状態を計算できていることがわかる。
【０２１６】
以上のように、本実施形態に係るシミュレーション装置は、上記のような変数設定手段と、変数抽出手段と、逆変換手段とを備え、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【０２１７】
本実施形態に係るシミュレーション方法は、上記シミュレーション装置によって実行される方法であって、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【０２１８】
本実施形態に係るプログラムおよび記録媒体は、上記シミュレーション方法を実行させることが可能であるから、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【０２１９】
「本発明の第２の実施形態」
次に第２の実施形態について説明する。本実施形態は、非特許文献１０におけるＳＣＣ法を援用した基礎方程式（式３９および式４０）またはＳＣＣ２法の基礎方程式（式４１から式４３）の解を求める上で、これらの基礎方程式とそれぞれ同値の方程式を解く点で、第１の実施形態と異なる。したがって、第１の実施形態と同様な構成要素の詳細な説明は必要のない限り省略する。
【０２２０】
本実施形態のシミュレーション装置１０も、第１の実施形態と同様に、入力手段１２と、制御手段１４と、座標設定手段１６と、座標抽出手段１８と、逆変換手段２０と、解析結果を表示する表示手段２２とを備える。
【０２２１】
そして、本実施形態では座標抽出手段１８は、式３９および式４０と同じ解を得るために、ＳＣＣ法における基礎方程式として下記式５７を使用する。
【０２２２】
【数３１】
0005697638.tif 000032
【０２２３】
Yは、下記式５８で定義されるＭＫ＋Ｍ＋Ｋ次元のベクトルである。つまり、Ｙは、ｉおよびμのそれぞれを走らせた各要素、ｍ１１／２ＲＳ１，・・・，ｍＭ１／２ＲＳＭと、φ１１，・・・，φＭ１，・・・，φ１Ｋ，・・・，φＭＫと、Λ１，・・・，ΛＫとを有する。なおφｉμおよびΛμは補助座標を表す。
【０２２４】
【数３２】
0005697638.tif 000033
【０２２５】
ｖμは、下記式５９の非斉次線形方程式の解ベクトルである。なお、ｖμはＭＫ＋Ｍ＋Ｋ次元のベクトルである。
【０２２６】
【数３３】
0005697638.tif 000034
【０２２７】
式５９におけるＣおよびｓμはそれぞれ式６０および式６１で定義される。なお、Ｃは（ＭＫ＋Ｍ＋Ｋ）×（ＭＫ＋Ｍ＋Ｋ）の正方行列である。Ｃは一見３×３の行列に見えるが、各成分に更に行列を有する。つまり、３×３の行列として表現されたＣのうち、１行１列の成分はＭＫ×Ｍ行列（つまり行数がＭＫで列数がＭ）であり、１行２列の成分はＭＫ×ＭＫ行列であり、１行３列の成分はＭＫ×Ｋ行列であり、２行２列の成分はＫ×ＭＫ行列であり、３行１列の成分はＭ×Ｍ行列であり、その他の成分はゼロ行列である。各成分の要素は、ｉあるいはμを行の足としてｉ＝１〜Ｍとμ＝１〜Ｋで走らせ、ｊあるいはνを列の足としてｊ＝１〜Ｍとν＝１〜Ｋで走らせて構成する。また、ｓμはＭＫ＋Ｍ＋Ｋ次元のベクトルである。ｓμも、Ｃと同様に一見３次元のベクトルに見えるが、各成分に更にベクトルを有する。つまり、ｓμの３次元ベクトルとしての第３成分はＭ次元ベクトルであり、その他の成分はゼロベクトルである。第３成分の要素は、ｉをｉ＝１〜Ｍで走らせて構成する。
【０２２８】
【数３４】
0005697638.tif 000035
【０２２９】
Ｖ,ｉｊｋ（ＲＳ）は式６２で定義される。
【０２３０】
【数３５】
0005697638.tif 000036
【０２３１】
上記式５７はＫが一般化された場合の基礎方程式である。そして、Ｋ＝１の場合には、式５７は下記式６３と同値となる。また、上記非特許文献１０に記載されているように、フロベニウスの定理を用いてＫ＝１の場合の式６３から、一般化されたＫについての式５７を導出することも可能である。
【０２３２】
【数３６】
0005697638.tif 000037
【０２３３】
式６３において（χ１，・・・，χＭ，ψ１，・・・，ψＭ，κ）は、式６４で表される非斉次線形方程式の２Ｍ＋１次元の解ベクトルである。
【０２３４】
【数３７】
0005697638.tif 000038
【０２３５】
なお式６４において、行列と固有ベクトルとの積を取り３つの方程式で表現したときにｊおよびｋについて１からＭの総和をとる。
【０２３６】
式６３および式６４が式３９および式４０と同値であることは以下のようにして説明することができる。式３８、式６２および式６３に注意して式６４を（φｉ１，Λ１）に関してｑ１で積分すると式６５が得られる。
【０２３７】
【数３８】
0005697638.tif 000039
【０２３８】
さらに、式６４に含まれる恒等式χｉ＝φｉ１を式６３の第１式に代入すると式６６が得られる。
【０２３９】
【数３９】
0005697638.tif 000040
【０２４０】
ここで、式６５の第１式および式６６におけるφｉ１およびΛ１をそれぞれφｉ（ＲＳ）およびΛ（ＲＳ）とおく（つまり、φｉ１およびΛ１をＲＳの関数として取り扱う）と、式３９および式４０が再現される。
【０２４１】
また、座標抽出手段１８は、式４１から式４３と同じ解を得る場合には、ＳＣＣ２法における基礎方程式として下記式６７を使用する。
【０２４２】
【数４０】
0005697638.tif 000041
【０２４３】
Ｚは、下記式６８で定義されるＭＫ＋Ｍ＋２Ｋ次元のベクトルである。つまり、Ｚは、ｉおよびμのそれぞれを走らせた各要素、ｍ１１／２ＲＳ１，・・・，ｍＭ１／２ＲＳＭと、φ１１，・・・，φＭ１，・・・，φ１Ｋ，・・・，φＭＫと、λ１，・・・，λＫと、ρ１，・・・，ρＫとを有する。なお、λμおよびρμは、φｉμと同様に補助座標を表す。また、ｃμνは、μごとに定義される下記式６９で表される値をそれぞれ最小にするように一意に決められる定数である。また、ｗμは、下記式７０で定義された（ＭＫ＋Ｍ＋Ｋ）×（ＭＫ＋Ｍ＋２Ｋ）次元の退化行列ＤのＫ次元の特異値空間を張るＫ個のＭＫ＋Ｍ＋２Ｋ次元の単位ベクトルである。Ｄは一見３×４の行列に見えるが、各成分に更に行列を有する。つまり、３×４の行列として表現されたＤのうち、１行１列の成分はＭＫ×Ｍ行列であり、１行２列の成分はＭＫ×ＭＫ行列であり、１行３列の成分はＭＫ×Ｋ行列であり、２行１列の成分はＭ×Ｍ行列であり、２行２列の成分はＭ×ＭＫ行列であり、２行４列の成分はＭ×Ｋ行列であり、３行２列の成分はＫ×ＭＫ行列であり、その他の成分はゼロ行列である。各成分の要素は、ｉあるいはμを行の足としてｉ＝１〜Ｍとμ＝１〜Ｋで走らせ、ｊあるいはνを列の足としてｊ＝１〜Ｍとν＝１〜Ｋで走らせて構成する。
【０２４４】
【数４１】
0005697638.tif 000042
【０２４５】
上記式６７はＫが一般化された場合の基礎方程式である。そして、Ｋ＝１の場合には、式６７は下記式７１と同値となる。また、上記非特許文献１０に記載されているように、フロベニウスの定理を用いてＫ＝１の場合の式７１から、一般化されたＫについての式６７を導出することも可能である。
【０２４６】
【数４２】
0005697638.tif 000043
【０２４７】
式７１において（χ１，・・・，χＭ，ψ１，・・・，ψＭ，κ，σ）は、式７２で表される（２Ｍ＋１）×（２Ｍ＋２）の退化行列による斉次線形方程式の２Ｍ＋２次元の解ベクトル（不定解）である。
【０２４８】
【数４３】
0005697638.tif 000044
【０２４９】
なお式７２において、行列と固有ベクトルとの積を取り３つの方程式で表現したときにｊおよびｋについて１からＭの総和をとる。つまり、（χ１，・・・，χＭ，ψ１，・・・，ψＭ，κ，σ）は、式７２の退化行列の１次元の特異値空間に属する単位ベクトルｕを用いて、式７３で表される。ここで、ｃは、式７４で表せられる値を最小にするように一意に決められる定数である。
【０２５０】
【数４４】
0005697638.tif 000045
【０２５１】
式７１および式７２が式４１から式４３と同値であることは以下のようにして説明することができる。式３８、式６２、式７１および式７３に注意して式７２を（φｉ１，λ１，ρ１）に関してｑ１で積分すると式７５が得られる。
【０２５２】
【数４５】
0005697638.tif 000046
【０２５３】
式７５の第１式および第２式からφｉ１およびρ１を消去すると式７６が得られる。そして、ＲＳｉとλ１をｑ１の関数と見なして式７６をｑ１で微分し、式７１の第１式および第３式を用いると式７７が得られる。
【０２５４】
【数４６】
0005697638.tif 000047
【０２５５】
ここで、式７１の第１式および第３式並びに式７７におけるλ１、χｉ、κおよびＶをそれぞれλ、φｉ（ＲＳ，λ）、κ（ＲＳ，λ）およびＶｅｆｆ（ＲＳ）とおくと、式４１から式４３が再現される。
【０２５６】
本実施形態では、例えば式６３または式７１から分かるように、基礎方程式を解く上で独立した座標として取り扱うｑ１の関数としての変数（中継座標ＲＳおよび補助座標の全体であり、以下独立した変数という。）の数が増えている。具体的には、式３９および式４０では独立した変数は、Ｍ個のＲＳｉ（ｉ＝１からＭ）のみであったのに対し、式３９および式４０と同値の式６３では独立した変数は、Ｍ個のＲＳｉ（ｉ＝１からＭ）と、Ｍ個のφｉ１（ｉ＝１からＭ）と、Λ１とから構成され、その総数は２Ｍ＋１個である。一方、式４１から式４３では独立した変数は、Ｍ個のＲＳｉ（ｉ＝１からＭ）およびλであったのに対し、式４１から式４３と同値の式７１では独立した変数は、Ｍ個のＲＳｉ（ｉ＝１からＭ）と、Ｍ個のφｉ１（ｉ＝１からＭ）と、λ１と、ρ１とから構成され、その総数は２Ｍ＋２個である。
【０２５７】
このように独立した変数（具体的には補助座標）の数を増やして方程式を同値変形することにより得られる利点は以下の通りである。
【０２５８】
例えば式４０からφｉ（ＲＳ）を求める場合、式４０がφｉ（ＲＳ）の符号反転に対して不変であることを反映し、φｉ（ＲＳ）の符号に任意性が残る。またこのことは、式４３からφｉ（ＲＳ，λ）およびκ（ＲＳ，λ）を求める場合にも同じである。この符号の任意性は、系の状態が反応経路を逆走する経路もまた反応経路を与えることを表す。一度逆走しても再度逆走することにより順方向に戻れば、最終的に反応経路を見出すことができることを考慮すれば、この符号の任意性はむしろ自然な現象を表現していると言える。しかしながら、この符号の任意性は、上記逆走によって反応経路を効率よく構成することが難しくなるという問題を引き起こす場合がある。さらに、式４０および式４３は実用上差分化して解かれるが、その差分化の程度Δｑ１が有限値であるため、逆走が生じても系の状態が元の反応経路上に戻らない場合があり、このような逆走は誤った反応経路が構成されてしまうきっかけとなるという問題もある。
【０２５９】
そこで、本実施形態では、独立した変数の数を増やすことによりφｉμ（ＲＳ）の符号の任意性を解消している。その結果、上記逆走の問題が解決され、より効率よくかつ正確に反応経路を構成することが可能となる。なお、独立した変数の数が２倍程度に増えたとしても、計算処理増加による負担はほとんど生じない。
【０２６０】
図１０は、符号の任意性が解消された計算の場合と符号の任意性が残ったままの計算の場合との比較結果を示すグラフである。具体的には、符号２４のグラフは、式７１から式７４を使用して、メリチンというタンパク質の構造変化ごとのポテンシャルエネルギーを計算した結果を示し、符号２６のグラフは、式４１から式４３を使用して、同じタンパク質の構造変化ごとのポテンシャルエネルギーを計算した結果を示す。なお、横軸は、構造変化を見やすくするために導入されたメリチンの屈曲角θであり、縦軸は系のポテンシャルエネルギーである。θ＝１４０ｄｅｇｒｅｅの状態が初期状態に相当する。
【０２６１】
符号の任意性が残ったままの計算では、θ＝９５ｄｅｇｒｅｅ程度でエネルギーが急上昇しており、この点で一度逆走してしまった。一方、符号の任意性が解消された計算ではそのような現象は見られず、θ＝８０ｄｅｇｒｅｅ程度でエネルギーの極大点に至った。グラフ２４におけるθ＝６０ｄｅｇｒｅｅの離れた点は、上記極大点における状態から系に全体構造緩和を適用して得られたもう１つの局所安定点である。
【０２６２】
以上のように、本実施形態に係るシミュレーション装置およびシミュレーション方法も、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものであるから、第１の実施形態と同様の効果を奏する。
【０２６３】
＜設計変更＞
上記各実施形態の基礎方程式において、ポテンシャルエネルギーＶｅｆｆ（ＲＳ）の３次微分の項が登場する（例えば、式４３、式６０、式６４、式７０および式７２）。この３次微分の計算処理に要する時間はオーダーでＭ３（Ｍは中継座標の個数）に比例するため、その２次微分の計算処理に要する時間がオーダーでＭ２に比例することを考慮すると、この３次微分の計算はシミュレーションを行う上で非常に大きな負担である。したがって、例えば式６４におけるＶ，ｉｊｋ（ＲＳ）・φｋ１について、Ｖ，ｉｊｋ（ＲＳ）を先に計算しその計算値とφｋ１の積をとると多くの時間を要してしまう。そこで、本発明者は下記式７８を開発し、２つの２次微分との差分をとることで直接Ｖ，ｉｊｋ（ＲＳ）・φｋ１を計算できるようにした。式７８を適用することにより、オーダーでＭ３程度となる計算時間をＭ２程度にまで短縮することが可能となる。さらに、式７８の右辺の２次微分に対して式４６のフィードバック方程式を適用することが好ましい。
【０２６４】
【数４７】
0005697638.tif 000048
【０２６５】
φｉは任意のベクトルのｉ番目の成分を表し、ｎは式７９で定義される。
【０２６６】
【数４８】
0005697638.tif 000049
【０２６７】
（ＳＣＣ法の基礎方程式について）
最後に、ＳＣＣ法における最も基本的な基礎方程式（式３４から式３６）と本明細書で例示した４組の基礎方程式（式３９および式４０の組、式４１から式４３の組、式５７並びに式６７）との関係について簡単に述べておく。
【０２６８】
前述したように一般的には、最も基本的な基礎方程式に従って平面を構成すれば、関数ＲＳｉ（ｑμ）が求められる。しかしながら、上記平面を構成する際に実は２つの問題がある。１つは、式３４から式３６を同時に厳密に満たす平面が一般のポテンシャルエネルギーＶに対して存在するとは限らないという原理的な問題である。そして、もう１つは、式３４から式３６に従い平面を構成しようとすると逆走が生じる場合があるという実用上の問題（前述した逆走の問題）である。特に、後者の問題の原因は明白で、式３５と式３６がφｉμ（ＲＳ）とρμ（ＲＳ）の符号を決定することができない、つまり符号反転に対して不変であるからである。
【０２６９】
そこで、この原理的な問題に対しては、式３４から式３６を満たす平面を近似的に見出す方法が有効である。実際に非特許文献１０では、式３６をあきらめ、式３４と式３５のみで平面を構成しようとした。そしてこの策は、Ｋ＝１の場合には式３９および式４０を与える。また、本発明者は、式３４を最大限満たしつつ、式３５と式３６で平面を構成する方法を開発した。そしてこの策は、Ｋ＝１の場合には式４１から式４３を与える。具体的には、式３５と式３６においてφｉ１（ＲＳ）とρ１（ＲＳ）を消去することで得られる下記式８０を、ＲＳの他、λ１もｑ１の関数とみなしてｑ１で微分すると式４１から式４３が得られる。
【０２７０】
【数４９】
0005697638.tif 000050
【０２７１】
一方、実用上の問題に対しては、式３４から式３６においてφｉμ、λμおよびρμをＲＳから独立しかつｑ†の関数としての変数（すなわち補助座標）とみなす方法が有効である。これは、φｉμ、λμおよびρμをＲＳから独立した変数とみなすと、φｉμ、λμおよびρμは、ＲＳと同様に、平面を構成する際に少しずつ連続的に変化する量となり、突然符号が変わるような事態が発生しないからである。φｉμ、λμおよびρμを補助座標として取り扱うためには、φｉμ、λμおよびρμについても式３４のようなｑνでの微分を含む式が必要である。そこで、式３５および式３６をそれぞれｑνで全微分し、これをｄｑνで割ることにより、下記式８１および式８２が得られる。これにより、最も基本的な基礎方程式は、式３４、式８１および式８２からなる実用上の問題のない基礎方程式へと変換されたことになる。
【０２７２】
【数５０】
0005697638.tif 000051
【０２７３】
そして、式３４、式８１および式８２からなる基礎方程式における原理的な問題に対処することを考える。その対処方法は、例えば上記と同様に、式３４、式８１および式８２を満たす平面を近似的に見出す方法である。具体的には、式３４および式８１のみで平面を構成する策は式５７を与える。一方、式３４を最大限満たしつつ、式８１と式８２で平面を構成する策は式６７を与える。
【産業上の利用可能性】
【０２７４】
各実施形態では、本発明をタンパク質および薬物からなる複合体の結合過程に適用した場合について説明したが、本発明はこれに限られるものではない。つまり、本発明は、より一般的に生体高分子および結合候補分子からなる複合体の形成過程について適用することが可能である。また、本発明は、複合体の形成過程の他、複合体の解離過程、タンパク質のａｐｏ体の構造および生体高分子のフォールディング（折り畳み機構）等、生体高分子を含む多原子系の挙動の解析についても適用することが可能である。
【符号の説明】
【０２７５】
２     タンパク質
４     ポケット
６     結合候補分子
１０   シミュレーション装置
１２   入力手段
１４   制御手段
１６   座標設定手段
１８   座標抽出手段
２０   逆変換手段
２２   表示手段
Ａ     複合体が解離体である状態を示す点
Ｂ     エネルギー障壁
Ｃ     結合が完了した複合体である状態を示す点
ＲＰ   複合体の挙動の反応経路
【図１】
0005697638.tif 000067
【図２】
0005697638.tif 000068
【図３】
0005697638.tif 000069
【図４】
0005697638.tif 000070
【図５】
0005697638.tif 000071
【図６Ａ】
0005697638.tif 000072
【図６Ｂ】
0005697638.tif 000073
【図７】
0005697638.tif 000074
【図８】
0005697638.tif 000075
【図９】
0005697638.tif 000076
【図１０】
0005697638.tif 000077