質点系の挙動を予測するシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体
【課題】生体高分子を含む系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することを可能とする。
【解決手段】シミュレーション装置10において、質点系の構造を記述する質点座標xに基づいて、中継座標RSおよび非中継座標RFを設定する座標設定手段16と、非中継座標RFを中継座標RSに従属させることにより非中継座標の構成RF(RS)を求め、中継座標RSの変化に伴う非中継座標RFの変化が中継座標RSに与える影響を考慮しながら、集団座標q†の関数としての中継座標の構成RS(q†)を求める座標抽出手段18と、運動方程式の解として得られる集団座標q†(t)、中継座標の上記構成RS(q†)および非中継座標の上記構成RF(RS)に基づいて、質点座標xの時間発展を予測する逆変換手段20とを備える。
【解決手段】シミュレーション装置10において、質点系の構造を記述する質点座標xに基づいて、中継座標RSおよび非中継座標RFを設定する座標設定手段16と、非中継座標RFを中継座標RSに従属させることにより非中継座標の構成RF(RS)を求め、中継座標RSの変化に伴う非中継座標RFの変化が中継座標RSに与える影響を考慮しながら、集団座標q†の関数としての中継座標の構成RS(q†)を求める座標抽出手段18と、運動方程式の解として得られる集団座標q†(t)、中継座標の上記構成RS(q†)および非中継座標の上記構成RF(RS)に基づいて、質点座標xの時間発展を予測する逆変換手段20とを備える。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、コンピュータを用いて、モデル化された質点系の動力学的な挙動を予測するシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体に関するものである。
【背景技術】
【0002】
コンピュータ技術の発展に伴って、シミュレーションによる理論的な計算により、タンパク質、核酸、脂質および多糖等の生体高分子の動力学的な大変形の挙動を原子レベルで理論的に解明しようという研究が盛んである。具体的には、標的とするタンパク質と結合候補分子(当該タンパク質に対する結合の適合性を有するか否かについての解析の対象となる分子)との互いの親和性を理論的に予測する創薬スクリーニング、および、タンパク質の一次配列から立体構造を構築する原理を明らかにして一次構造から高次構造を理論的に構築するタンパク質のフォールディング(折り畳み機構)の解析が挙げられる。
【0003】
例えば、生体高分子の動力学的な挙動を予測するシミュレーション方法として、巨大分子でも原子レベルでシミュレーションを実行することが可能な分子動力学法またはモンテカルロ法等に基づくシミュレーション方法(非特許文献1から7)がある。分子動力学法によれば、運動方程式に従って多原子系の時間発展を微小時間単位で順次追跡することができる。しかし、生体高分子を含む多原子系のポテンシャル曲面には、多くのローカルミニマム(極小点)或いは多くのエネルギー障壁が存在しているため、上記のような方法では、多原子系の状態が初期構造に近いローカルミニマムに捕捉されて、計算に膨大な時間を要するという問題があった。また、この問題はモンテカルロ法でも同じである。
【0004】
そこで、上記のようなローカルミニマムによるトラップの問題を解消する方法として、分子動力学法に基づく計算において「強制的な力」を印加することによりローカルミニマムによるトラップを脱出させる計算方法が知られている。例えば、特許文献1から3並びに非特許文献8および9では、多原子系が一度経験した構造を二度と取らないように仮想的な相互作用を逐次導入することや、高温相の多原子系の加速された運動で見つかった構造を低温相の多原子系の運動に逐次反映させることが開示されている。このような相互作用の導入等により、多原子系が通過するポテンシャル曲面上の経路の探索を加速して、生体高分子の動力学的な挙動を計算することが可能となる。
【先行技術文献】
【特許文献】
【0005】
【特許文献1】特開2005−267592号公報
【特許文献2】特開2005−524129号公報
【特許文献3】再表2006/068271号公報
【非特許文献】
【0006】
【非特許文献1】Todd J. A. Ewing, Irwin D. Kuntz, Journal of Computational Chemistry, Vol. 18, Issue 9 (1997), p. 1175-1189
【非特許文献2】Garrett M. Morris, David S. Goodsell, et al., Journal of Computational Chemistry, Vol. 19, Issue 14 (1998), p. 1639-1662
【非特許文献3】Matthias Rarey, Bernd Kramer, et al., Journal of Molecular Biology, Vol. 261, Issue 3 (1996), p. 470-489
【非特許文献4】Ruben Abagyan, Maxim Totrov, et al., Journal of Computational Chemistry, Vol. 15, Issue 5 (1994), p. 488-506
【非特許文献5】Gareth Jones, Peter Willett, et al., Journal of Molecular Biology, Vol. 267, Issue 3 (1997), p. 727-748
【非特許文献6】Richard A. Friesner, Jay L. Banks, et al., Journal of Medicinal Chemistry, Vol. 47, Issue 7 (2004), p. 1739-1749
【非特許文献7】Thomas A. Halgren, Robert B. Murphy, et al., Journal of Medicinal Chemistry, Vol. 47, Issue 7 (2004), p. 1750-1759
【非特許文献8】Yoshifumi Fukunishi, Yoshiaki Mikami, et al., Journal of Physical Chemistry B, Vol. 107, Issue 47 (2003), p. 13201-13210
【非特許文献9】Yuji Sugita, Yuko Okamoto, Chemical Physics Letters, Vol. 314, Issue 1-2 (1999), p. 141-151
【発明の概要】
【発明が解決しようとする課題】
【0007】
しかしながら、生体高分子等がモデル化された質点系の運動を記述する運動方程式に基づいて、時間に関する積分で当該質点系の挙動を予測する分子動力学を用いた計算では、シミュレーションの計算時間および計算精度がトレードオフの関係になってしまうという問題がある。具体的には、特許文献1から3並びに非特許文献8および9のようなシミュレーション方法では、計算時間は短縮されるものの、印加した「強制的な力」によって質点系の本来有する運動秩序が破壊され、計算結果として非現実的な大変形の挙動を算出してしまう。そして、より効率よくローカルミニマムによるトラップを脱出させるために「強制的な力」の作用を大きくすればするほど、その傾向は顕著となる。
【0008】
本発明は上記問題に鑑みてなされたものであり、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することを可能とするシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体を提供することを目的とするものである。
【課題を解決するための手段】
【0009】
上記課題を解決するために、本発明に係るシミュレーション装置は、
モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション装置において、
質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定する座標設定手段と、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるK個の集団座標の関数としての中継座標の構成を求める座標抽出手段と、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測する逆変換手段とを備えることを特徴とするものである。
【0010】
なお本明細書において、K、MおよびNは、K<M<3Nの関係を満たし、それぞれ1以上の整数を表す。
【0011】
本明細書において、「質点系の構造」とは、当該質点系を構成するN個の質点が造り出す立体構造を意味する。
【0012】
「質点系の構造の変化を主として担う中継座標」とは、当該座標の変化が当該質点系の構造の立体構成に大きな影響を与えるような座標を意味する。
【0013】
「中継座標を設定し」とは、中継座標として、一部の質点座標そのもの、質点座標を組み合わせて定義できる座標又はこれらの組み合わせを設定することを意味する。非中継座標についても同様である。
【0014】
そして、本発明に係るシミュレーション装置において、座標抽出手段は、
中継座標および非中継座標の関数として表現された質点系のポテンシャルエネルギーを求める第1の工程を実施し、
ポテンシャルエネルギーを用いた断熱近似条件に従って非中継座標を中継座標に従属させる第2の工程を実施し、
現状の中継座標および非中継座標の下で、上記影響を考慮しながら、中継座標に関するポテンシャルエネルギーの微分係数を求める第3の工程を実施し、
ポテンシャルエネルギーの上記微分係数の下で、自己無撞着集団座標法におけるポテンシャルエネルギーの微分係数を用いた基礎方程式に従って、集団座標に関する中継座標の微分係数を求める第4の工程を実施し、
中継座標の上記微分係数の下で集団座標を微小量だけ更新することにより、更新された中継座標を求める第5の工程を実施し、
更新された中継座標の下で、中継座標に従属した非中継座標について構造緩和を行う第6の工程を実施し、
その後、非中継座標の構造緩和後の状態における中継座標および非中継座標に基づいて、第3の工程から第6の工程までを順次繰り返すことにより、中継座標の上記構成を求めるものであることが好ましい。
【0015】
「現状の中継座標および非中継座標の下で」とは、1回目の第3の工程においては、座標設定手段により設定された状態における中継座標および非中継座標の下で、2回目以降の第3の工程においては、直前の第6の工程で行われた非中継座標の構造緩和後の状態における中継座標および非中継座標の下で、その第3の工程を実施することを意味する。
【0016】
また、本発明に係るシミュレーション装置において、上記影響を考慮する方法は下記式1から式3の少なくとも1つを使用することであることが好ましい。
【0017】
【数1】
【0018】
本明細書において、i、jおよびkは、それぞれ1からMまでの整数を表し、
α、βおよびγは、それぞれ1から3N−Mまでの整数を表し、
RSiは、質点系におけるi番目の中継座標を表し、
RFαは、質点系におけるα番目の非中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
RFは(RF1,RF2,・・・,RF(3N−M))を表し、
RF(RS)は、中継座標に従属させられた非中継座標を表し、
V(RS,RF)は、中継座標と非中継座標で表された質点系のポテンシャルエネルギーを表し、
Veff(RS)は、V(RS,RF)にRF(RS)を代入して得られる有効ポテンシャルエネルギーを表す。
【0019】
式2において、第3項の(i←→j)は、第2項におけるiおよびjの文字を相互に入れ替えた項(つまり、第2項におけるiをjへ、jをiへ書き換えた項。以下同様である。)を表す。
【0020】
式3において、第3項の(i←→k)は、第2項におけるiおよびkの文字を相互に入れ替えた項を表し、第4項の(j←→k)は、第2項におけるjおよびkの文字を相互に入れ替えた項を表し、第6項の(i←→k)は、第5項におけるiおよびkの文字を相互に入れ替えた項を表し、第7項の(j←→k)は、第5項におけるjおよびkの文字を相互に入れ替えた項を表す。
【0021】
さらに、式1から式3において下記式4を用いた。
【0022】
【数2】
【0023】
Kαβ−1(Rs)はKαβ(Rs)の逆行列を表し、
Kαβ(Rs)およびJαi(Rs)はそれぞれ式5および式6で定義される。
【0024】
【数3】
【0025】
また、本発明に係るシミュレーション装置において、断熱近似条件は下記式7であることが好ましい。
【0026】
【数4】
【0027】
また、本発明に係るシミュレーション装置において、集団座標の個数KはK=1を満たすものであり、自己無撞着集団座標法における基礎方程式は下記式8および9であることが好ましい。
【0028】
【数5】
【0029】
本明細書において、q1は、集団座標を表し、
miは質点系におけるi番目の中継座標に関する質量を表し、
φi(RS)は、式9を満たす関数(固有ベクトル)のi番目の成分を表し、
Λ(RS)は、式9を満たす関数(固有値)を表す。
【0030】
或いは、本発明に係るシミュレーション装置において、集団座標の個数KはK=1を満たすものであり、自己無撞着集団座標法における基礎方程式は下記式10から12であることが好ましい。
【0031】
【数6】
【0032】
φi(RS,λ)およびκ(RS,λ)は式12を満たす関数であり、φi(RS,λ)はそのi番目の成分を表し、
λは、補助座標(中継座標から独立しかつ集団座標の関数として取り扱う変数)を表す。
【0033】
或いは、本発明に係るシミュレーション装置において、座標抽出手段は、第4の工程において、基礎方程式中の集団座標に関する中継座標または補助座標の微分係数が有する符号の任意性を解消するように、この基礎方程式を解く上で中継座標から独立しかつ集団座標の関数として取り扱う変数(つまり補助座標)の数を増やして、計算を行うものであることが好ましい。
【0034】
また、補助座標の数を増やす場合において、座標抽出手段は、その結果得られた下記式13で表される基礎方程式に従って計算を行うものであることが好ましい。さらにこの場合において集団座標の個数KはK=1を満たすことが好ましい。
【0035】
【数7】
【0036】
Yは、下記式14で定義されるMK+M+K次元のベクトルであり、vμは、下記式15の非斉次線形方程式の解ベクトルである。
【0037】
【数8】
【0038】
式15におけるCおよびsμはそれぞれ式16および式17で定義される。
【0039】
【数9】
【0040】
V,ij(RS)およびV,ijk(RS)はそれぞれ式18および式19で定義される。
【0041】
【数10】
【0042】
なお、μおよびνはそれぞれ1からKまでの整数を表し、
qμはμ番目の集団座標を表し、
φiμおよびΛμは補助座標を表す。
【0043】
或いは、補助座標の数を増やす場合において、座標抽出手段は、その結果得られた下記式20で表される基礎方程式に従って計算を行うものであることが好ましい。さらにこの場合において集団座標の個数KはK=1を満たすことが好ましい。
【0044】
【数11】
【0045】
Zは、下記式21で定義されるMK+M+2K次元のベクトルであり、cμνは、μごとに定義される下記式22で表される値をそれぞれ最小にするように一意に決められる定数であり、wμは、下記式23で定義された行列DのK次元の特異値空間を張るK個のMK+M+2K次元の単位ベクトルである。なお、λμおよびρμは、φiμと同様に補助座標を表す。
【0046】
【数12】
【0047】
また、本発明に係るシミュレーション装置において、座標抽出手段は、ポテンシャルエネルギーの3次微分の項の計算を下記式24に基づいて行うものであることが好ましい。φiは任意のベクトルのi番目の成分を表し、nは式25で定義される。
【0048】
【数13】
【0049】
また、本発明に係るシミュレーション装置において、座標設定手段は、質点系の構造のうち特徴的な部分構造ごとに抽出された当該部分構造を代表する代表座標を中継座標として設定するものであることが好ましい。
【0050】
本明細書において、「特徴的な部分構造」とは、質点系の構造のうち、形態的および/または機能的な特徴を有する部分構造を意味する。
【0051】
また、本発明に係るシミュレーション装置において、質点系は生体高分子を含む多原子系であり、
上記部分構造は、生体高分子の2次構造、生体高分子の構成単位または生体高分子の主鎖であり、
上記部分構造についての代表座標は、その部分構造を構成する原子の座標そのもの、その原子の座標を組み合わせて定義できる座標またはその部分構造の周期間隔とすることができる。
【0052】
また、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子はタンパク質であり、上記部分構造はタンパク質の2次構造であり、この2次構造についての代表座標は、その2次構造を構成する原子群の重心座標または2次構造の屈曲角であることが好ましい。この場合において、2次構造はヘリックス構造、βシート、ターン、ループおよびランダムコイルの少なくとも1つであることが好ましい。
【0053】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子がタンパク質であり、上記部分構造がタンパク質の残基であり、この残基についての代表座標が、その残基を構成する原子群の重心座標であることが好ましい。
【0054】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子がタンパク質であり、上記部分構造がタンパク質の主鎖であり、この主鎖についての代表座標が、その主鎖を構成する原子のそれぞれの座標であることが好ましい。
【0055】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の2次構造であり、この2次構造についての代表座標が、その2次構造を構成する原子群の重心座標または2次構造の屈曲角であることが好ましい。この場合において、2次構造は螺旋構造であることが好ましい。
【0056】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の残基であり、この残基についての代表座標が、その残基を構成する原子群の重心座標であることが好ましい。
【0057】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の主鎖であり、この主鎖についての代表座標が、その主鎖を構成する原子のそれぞれの座標であることが好ましい。
【0058】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の螺旋構造であり、この螺旋構造についての代表座標が、その螺旋構造の周期間隔であることが好ましい。
【0059】
また、本発明に係るシミュレーション装置において、多原子系はさらに、生体高分子に対する結合候補分子を含むものとすることができる。
【0060】
本発明に係るシミュレーション方法は、
上記に記載のシミュレーション装置によって実行される、モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション方法であって、
質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、
質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定し、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、
中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるK個の集団座標の関数としての中継座標の構成を求め、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測することを特徴とするものである。
【0061】
本発明に係るシミュレーションプログラムは、上記に記載のシミュレーション方法をコンピュータに実行させることを特徴とするものである。
【0062】
本発明に係るコンピュータ読み取り可能な記録媒体は、上記に記載のシミュレーションプログラムを記録したものである。
【発明の効果】
【0063】
本発明に係るシミュレーション装置は、上記のような座標設定手段と、座標抽出手段と、逆変換手段とを備え、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0064】
本発明に係るシミュレーション方法は、上記シミュレーション装置によって実行される方法であって、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0065】
本発明に係るプログラムおよび記録媒体は、上記シミュレーション方法を実行させることが可能であるから、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【図面の簡単な説明】
【0066】
【図1】本発明のシミュレーション方法により誘導適合を考慮した場合における、タンパク質と結合候補分子の結合過程を示す概略図である。
【図2】タンパク質および結合候補分子を含む多原子系の誘導適合を考慮した概念的なポテンシャル曲面における反応経路を示す概略図である。
【図3】集団座標の抽出を利用して多原子系の反応経路を求める手順を概念的に示す図である。
【図4】本発明における変数変換および逆変換の概念を示す図である。
【図5】実施形態のシミュレーション装置の構成を概略的に示すブロック図である。
【図6A】実施形態のシミュレーション方法における計算手順を概略的に示すフローチャート図である。
【図6B】実施形態のシミュレーション方法における計算手順を概略的に示すフローチャート図である。
【図7】実施形態における中継座標RS、非中継座標RFおよび原子の座標xの関係を示す概略図である。
【図8】実施形態のシミュレーション方法によって得られた所定の複合体の反応経路を示すグラフである。
【図9】図8のグラフにおける複合体の結合過程を示す概略図である。図9aは複合体の結合過程における開始状態を示し、図9bは複合体の結合過程における最終状態を示す。
【図10】符号の任意性が解消された計算の場合と符号の任意性が残ったままの計算の場合との比較結果を示すグラフである。
【発明を実施するための形態】
【0067】
以下、本発明の実施形態について図面を用いて説明するが、本発明はこれに限られるものではない。なお、視認しやすくするため、図面中の各構成要素の縮尺等は実際のものとは適宜異ならせてある。
【0068】
まず、詳細な実施形態の説明を行う前に、従来技術に対して本発明が有する技術上の意義を明らかにするべく、本発明の基本となる技術的思想およびそれに至った経緯を概説する。ここでは、説明を分かりやすくするために、具体的に生体高分子および結合候補分子からなる複合体を含む多原子系の動力学的な挙動を予測する場合を例にして説明する。
【0069】
上記のような多原子系の動力学的な挙動の予測は、例えば、開発の期間短縮およびコスト削減等の観点から新薬開発において重要な役割を果たす。薬物等の生理活性物質は特定のタンパク質に結合することによりその化学的性質を発揮する。したがって、新薬開発においては、標的とするタンパク質に結合しやすい薬物候補を数十万から数百万とも言われる膨大な化合物から絞り込む必要がある(創薬スクリーニング)。また究極的には、薬物候補を膨大な化合物から数個程度にまで絞り込むことが要求されている。そこで、薬物候補を膨大な化合物から効果的に絞り込むためには多原子系の挙動を原子レベルで取り扱い、分子間さらには原子間の相互作用を正確に評価する必要がある。
【0070】
図1は、タンパク質の動力学的な挙動を考慮した場合における、タンパク質および結合候補分子の結合過程を示す概略図である。タンパク質2は、結合候補分子6と結合するためのポケット4を有している。つまり、当該結合候補分子6は、当該タンパク質2の生理活性物質と言える。しかしながら、当該ポケット4は、必ずしもそのままで結合候補分子6がその中に収まることができるような形状を有しているわけではない(図1a)。例えば図1では、ポケット4の開口が結合候補分子6の大きさよりも小さいことが示されている。そこで、このような場合には、タンパク質2は、接近した結合候補分子6との相互作用に応じて、開口の形状および大きさを結合候補分子6の形状および大きさに適合させながらその構造を変化させる(図1bおよびc)。
【0071】
このように、タンパク質および核酸等の生体高分子が生理活性物質との相互作用によりその構造を変化させる現象は「誘導適合」と言われる。分子間の相互作用を正確に計算するためには、当然この誘導適合を考慮しなければならない。そして、これは、多原子系を原子レベルで取り扱うことによって実現できる。
【0072】
一方、図2は、多原子系の誘導適合を考慮した概念的なポテンシャル曲面上の反応経路を示す概略図である。図2の横軸X1はタンパク質の構造の変化を概念的に示す軸であり、縦軸X2はタンパク質および結合候補分子の距離を示す軸である。つまり、図2は、タンパク質の構造の変化および分子間の距離に応じた多原子系のポテンシャル曲面を表している。図2に示されるように、誘導適合が生じる系は、ポテンシャル曲面上の2つの安定点(つまり、複合体が解離体である状態を示す点Aおよび複合体の結合が完了した状態を示す点C)を結ぶ反応経路RP上にエネルギー障壁B(鞍点)を有することが解かる。
【0073】
したがって、誘導適合を考慮しながらタンパク質2と結合候補分子6との結合過程を解析するということは、タンパク質2および結合候補分子6に関する上記ポテンシャル曲面上において、エネルギー障壁Bを超えて2つの安定点AおよびCを結ぶ反応経路を探索することに帰着する。
【0074】
以上を考慮すると、多原子系の動力学的な挙動を予測するシミュレーションにおいて計算精度を向上させかつ計算時間を短縮するためには、多原子系の原子レベルでの取り扱いを可能にすること、および非経験的な反応経路の探索を容易にすることが必要となる。
【0075】
しかしながら、数千から数百万個の原子で構成される多原子系の挙動は、1秒〜1時間オーダーの時間規模で起こる緩慢な大変形であるため、従来のシミュレーション方法では、原子レベルでの取り扱いは可能であるものの、理論的な計算を行うために数十年もの膨大な時間が必要であるという問題があった。そこで、多原子系の挙動の理論的な計算を実用的な期間内で実施できるようにするために、様々なシミュレーション手法が検討・開発されている。前述した「強制的な力」を印加する分子動力学法もその1つであるが、その問題点は既に述べた通りである。また、標的とするタンパク質を剛体として近似して、取り扱う変数を減らすことで計算量を減らす手法も検討されている。しかし、このような粗い近似の下では、タンパク質の動力学的な挙動の影響を当然考慮できず、結果としてこれらの分子間に働く相互作用を正しく計算できない。このような場合、例えば100万個ある候補の化合物から薬物候補として1万個程度に絞り込むことしかできない。また、その他のシミュレーション手法でも計算精度および計算時間のトレードオフの問題は解決するに至っていない。
【0076】
そこで、本発明者は、生体高分子を含む多原子系をN個の質点から構成される質点系にモデル化した上で、この質点系の挙動を集団運動として取り扱い、質点座標(質点系の構造を一次的に記述する質点の座標)の自由度より小さい自由度の集団座標を、その質点座標に基づいた集団運動の中から抽出することに着目した。
【0077】
(集団座標)
集団座標について説明する。一般的に「集団座標」とは、集団変数の座標要素を意味する。そして、「集団変数」とは、質点系の正準変数(その質点系における質点座標および質点の運動量)と正準変換により対応づけられる一般変数(q,p)のうち、その正準変数よりも自由度が小さいものを意味する。なお、正準変換により対応づけられる一般変数(q,p)において、qは座標要素であり、pは運動量要素である。qおよびpは下記式26および式27によってそれぞれ定義される。
【0078】
【数14】
【0079】
式26および式27において、ηは1から3Nまでの整数を表す。
【0080】
したがって、集団変数は、下記式28および式29によってそれぞれ定義されるq†およびp†を用いて(q†,p†)と表すことができる。
【0081】
【数15】
【0082】
式28および式29は、集団変数(q†,p†)が式28および式29に示されたμ=1〜Kまでの2K個の変数により構成されることを表している。言い換えれば、集団変数(q†,p†)とは、時間に依存して変化する可変成分(q1,q2,・・・,qK,p1,p2,・・・,pK)および時間変化に対して定数となる不変成分(qK+1=0,qK+2=0,・・・,q3N=0,pK+1=0,pK+2=0,・・・,p3N=0)から構成される一般座標の当該可変成分であると言える。なお、式28および式29では、不変成分が定数としてゼロの値を取る場合を示しているが、定数はゼロに限られない。
【0083】
そして、集団座標は、集団変数(q†,p†)の座標要素であるから、q†と表すことができる。より具体的には、集団座標とは、(q1,q2,・・・,qK)から構成される変数の集合である。なお、集団座標自体は、質点座標のみから正準点変換により求めることもできる。
【0084】
上記のように、集団運動として取り扱う対象となる座標(例えば、質点座標および後述する中継座標)からより自由度の低い集団座標へと変数変換を行うことを「変数分離」ともいう。
【0085】
集団座標の自由度は、変数分離の対象となる座標の自由度よりも小さければ、特に限定されない。しかしながら、集団座標の自由度が小さい程、その後の計算処理が容易となるため、集団座標の自由度は1(つまりK=1)であることが好ましい。
【0086】
(集団運動としての取り扱い)
質点座標に対して変数分離を実施すると、集団座標q†によって、質点系の挙動における本質的な運動、つまり、より低次元的な集団運動の記述が可能となる。この結果、3N次元の質点座標による運動の記述がK次元の集団座標q†による運動の記述へと置き換えられ、質点系の運動の記述が単純化される。具体的には、質点座標で記述された運動方程式は、集団座標q†を引数として持つ質点座標xの構成x=x(q†)とq†についての運動方程式との組み合わせに単純化される。なお、xは(x1,x2,・・・,x3N)を表す。質点座標x(q†)の具体的な構成は質点系の集団運動としての挙動の仕方を表し、集団座標q†についての運動方程式はその挙動の基本法則を表す。したがって、運動方程式を解く上で、取り扱う変数の自由度がより減少することになるため、理論的な計算が容易になる。
【0087】
そして、質点座標xに対して変数分離を実施して集団座標q†に関するポテンシャル曲面上の反応経路を見つけ出し、座標を集団座標q†から質点座標xに逆変換して、質点座標xに関するポテンシャル曲面上で反応経路を表現することにより、質点系の反応経路を求めることができる。上記のような方法によれば、低次元の運動方程式を解いた後、集団座標q†から質点座標xへの変数変換を実行するだけであるため、本質的に質点系の状態がローカルミニマムに捕捉されることが起こり得ない。
【0088】
例えば、図2に示された複合体を含む多原子系の反応経路を求める場合の概念的な手順は以下のようになる。図3は、集団座標の抽出を利用して多原子系の反応経路を求める手順を概念的に示す図である。まず、質点座標X1およびX2に関するポテンシャル曲面(図3a)が、抽出された集団座標q1およびq2に関するポテンシャル曲面(図3b)に変換される。この際、集団座標として、2つの安定点AおよびC並びにエネルギー障壁Bが一直線上に並ぶような集団座標を採用することが重要である。このような集団座標q1およびq2に関するポテンシャル曲面上では、2つの安定点AおよびCを結ぶ反応経路RPを一気に構成することができる(図3b)。その後、座標の逆変換により、変数X1およびX2に関するポテンシャル曲面上での反応経路RPが求められる(図3c)。ここで得られた一次元的な軌跡こそ理論化学の分野で古くから用いられてきた「反応座標」と同定できる変数である。
【0089】
(生体高分子を含む多原子系の挙動を集団運動として取り扱う際の課題)
変数分離は集団運動理論に基づいて行われる。例えば、集団運動理論の1つとして、自己無撞着集団座標法(SCC法:Self-consistent Collective Coordinate method)がある。SCC法は、現在原子核の分野で研究が進められている物理学における集団運動理論の1つである。より具体的には、SCC法は、系の集団運動を記述する本質的な集団変数を求め、その系のハミルトニアンに含まれる正準変数が張る6N次元の位相空間から、動力学的な集団運動に関与する分離多様体を見出し、その分離多様体上またはその近傍で定義されたハミルトニアンによって集団運動を記述する運動理論である。「分離多様体」とは、その系のハミルトニアンに含まれる正準変数が張る6N次元の位相空間中の部分空間であって、その系の挙動を表す軌道が閉じ込められる部分空間を意味する。つまり、分離多様体中の任意の1点を初期値とした場合、上記軌道は常にこの分離多様体中に拘束されることになる。SCC法に関する詳細については、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らの文献を非特許文献10という。)等を参照されたい。
【0090】
そこで、例えば原子核反応において取り扱う核子(陽子および中性子)を化学反応において取り扱う原子に置き換えてSCC法を適用することにより、多原子系の動力学的な挙動を集団運動として取り扱ってこの多原子系の反応経路を求めることができるようにも思える。
【0091】
しかし、本発明者は、本発明の対象である生体高分子を含むような大規模な系に対してSCC法等の既存の集団運動理論をそのまま適用するだけでは、計算精度を向上させかつ計算時間を短縮するには不十分であるという問題を新たに見出した。これは、原子核反応の場合にはせいぜい数百個の核子を取り扱えばよいのに対して、水中における生体高分子を含む多原子系の化学反応の場合には数千から数百万個の原子を取り扱わなければならないことに起因する。つまり、集団運動理論とはいえ既存の手法では変数の制御、つまり変数分離の計算が複雑になるためである。
【0092】
そして、上記のような問題は、複合体を含む多原子系の反応経路を求める場合に限られるものではない。つまり、より一般的に、生体高分子を含む系の反応経路を求める場合に同様の問題が生ずることは容易に想像される。
【0093】
「発明の説明」
以上のような経緯を踏まえ、本発明者は、質点座標の自由度より小さい自由度の集団座標によって集団運動を記述する集団運動理論であって、巨大な原子群である生体高分子を含む多原子系の構造計算にも適用可能な新規な集団運動理論に基づいて、質点系の動力学的な挙動を予測するシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体を発明するに至った。
【0094】
より具体的には、本発明に係るシミュレーション装置は、
モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション装置において、
質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定する座標設定手段と、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるK個の集団座標の関数としての中継座標の構成を求める座標抽出手段と、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測する逆変換手段とを備えることを特徴とするものである。
【0095】
そして、本発明に係るシミュレーション方法は、
上記に記載のシミュレーション装置によって実行される、モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション方法であって、
質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、
質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定し、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、
中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるK個の集団座標の関数としての中継座標の構成を求め、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測することを特徴とするものである。
【0096】
また、本発明に係るシミュレーションプログラムは、上記に記載のシミュレーション方法をコンピュータに実行させることを特徴とするものである。
【0097】
また、本発明に係るコンピュータ読み取り可能な記録媒体は、上記に記載のシミュレーションプログラムを記録したものである。
【0098】
なお、「N個の質点から構成された質点系」とは、質点系の構成要素である質点の総数がN個であることを意味する。
【0099】
<座標設定手段>
座標設定手段は、質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定するものである。上記「3N個の質点座標」とは、より具体的には3次元空間におけるN個の質点の座標xである。
【0100】
まず、座標設定手段は、質点座標の構成に基づいて、質点系の巨視的な集団運動(大振幅集団運動)を記述する上で質点系の構造の変化を主として担うM個の中継座標を設定する。そして、後述する座標抽出手段はこの設定された中継座標に対して変数分離を実施することになる。言い換えれば、本発明では、質点系を記述する独立した座標のうち大振幅集団運動による質点系の構造変化に対して影響の少ない座標は既知であるので、このような座標は集団運動理論における変数分離の対象から予め除外される。つまり、3N次元の質点座標から集団座標を抽出するのではなく、M次元の中継座標から集団座標を抽出することにより、変数分離の計算が単純化される。Mは、K<M<3Nを満たす整数である。また、Kは、上記式28および式29での説明同様、集団座標q†の要素の個数である。
【0101】
なお、この中継座標からさらに変数の階層化を行ってもよい。つまり、上記のように設定された中継座標の構成に基づいて、さらに低次元の2次的な中継座標を設定し、この2次的な中継座標から集団座標を抽出することも可能である。
【0102】
中継座標は、大振幅集団運動による質点系の構造変化に対する影響の度合を基準に設定される。そして、この影響の度合が大きい座標ほど中継座標としてふさわしい。言い換えれば、どのような質点系の大振幅集団運動を解析の対象としているのかによって、中継座標の設定方法および具体的な座標の内容が変わる。
【0103】
具体的には、中継座標は、質点系に含まれる一部の質点の座標そのもの、質点の座標を組み合わせて定義できる座標又はこれらの組み合わせを用いて設定される。また中継座標は、例えば質点系の構造のうち特徴的な部分構造に関連して設定されることが好ましい。そしてこの場合において、中継座標は、特徴的な部分構造ごとに抽出された、当該部分構造を代表する代表座標であることがより好ましい。代表座標は、1つの特徴的な部分構造に複数あってもよい。
【0104】
「特徴的な部分構造」とは、質点系の構造のうち、形態的および/または機能的な特徴を有する部分構造を意味する。特徴的な部分構造は、例えば質点系が生体高分子を含む多原子系である場合には、生体高分子の2次構造(部分的な折り畳み構造)、生体高分子の構成単位および生体高分子の主鎖である。また、特徴的な部分構造の「代表座標」とは、その特徴的な部分構造を代表して表す座標を意味する。代表座標は、例えば、その特徴的な部分構造を構成する質点(原子)の座標そのもの、その質点(原子)の座標を組み合わせて定義できる座標またはその部分構造の周期間隔等である。
【0105】
より具体的には、上記のような特徴的な部分構造としては、生体高分子がタンパク質である場合には、タンパク質の2次構造が挙げられる。具体的な2次構造としては、ヘリックス構造(310ヘリックス、αヘリックス、πヘリックスおよびβヘリックス等)、βシート、ターン、ループおよびランダムコイル等が挙げられる。そして、タンパク質の2次構造についての代表座標としては、その2次構造を構成する原子群の重心座標または2次構造の屈曲角が挙げられる。例えば、一般的なタンパク質中に含まれる上記のような特徴的な部分構造(高次構造)の個数は、せいぜい数個から数十個である。これにより、変数分離の対象となる変数の個数を大幅に低減することができる。また、生体高分子と結合候補分子との結合反応における挙動を予想する場合には、結合候補分子の構造自体を多原子系の特徴的な部分構造の1つとして取り扱うことも可能である。
【0106】
或いは、特徴的な部分構造としては、生体高分子がタンパク質である場合には、タンパク質の残基(タンパク質の構成単位であるアミノ酸部分。N末端残基およびC末端残基を含む)とすることもできる。そして、タンパク質の残基についての代表座標としては、その残基を構成する原子群の重心座標が挙げられる。
【0107】
或いは、特徴的な部分構造としては、生体高分子がタンパク質である場合には、タンパク質の主鎖とすることもできる。そして、タンパク質の主鎖についての代表座標としては、主鎖を構成する原子のそれぞれの座標が挙げられる。
【0108】
或いは、特徴的な部分構造としては、生体高分子が核酸である場合には、核酸の2次構造とすることもできる。具体的な2次構造としては、螺旋構造等が挙げられる。螺旋構造の場合には、所定の周期ごとに特徴的な部分構造として取り扱う。そして、核酸の2次構造についての代表座標としては、その2次構造を構成する原子群の重心座標または2次構造の屈曲角が挙げられる。
【0109】
或いは、特徴的な部分構造としては、生体高分子が核酸である場合には、核酸の残基(核酸の構成単位であるヌクレオチド部分)とすることもできる。そして、核酸の残基についての代表座標としては、その残基を構成する原子群の重心座標が挙げられる。
【0110】
或いは、特徴的な部分構造としては、生体高分子が核酸である場合には、核酸の主鎖とすることもできる。そして、核酸の主鎖についての代表座標としては、主鎖を構成する原子のそれぞれの座標が挙げられる。
【0111】
或いは、生体高分子が核酸でありかつ核酸の螺旋構造である場合には、螺旋構造の周期間隔を代表座標とすることもできる。
【0112】
特徴的な部分構造を実際にどのように設定するかについては、シミュレーションを行う対象の質点系の構成および/または現象に応じて適宜決定される。例えば、タンパク質の折りたたみ構造の形成過程や核酸の螺旋構造の形成過程を解析の対象とする場合には、生体高分子の構成単位や生体高分子の基本的な骨格とも言える主鎖を特徴的な部分構造として採用することが好ましい。生体高分子の主鎖については、共有結合により連なった一本の主鎖全体を1つの特徴的な部分構造としてもよいし、上記一本の主鎖を分割して得られた部分的な構造のそれぞれを特徴的な部分構造としてもよい。また、タンパク質と結合候補分子との結合反応における挙動を解析の対象にする場合には、上記で挙げた主鎖の他、へリックス構造、βシート、ターン、ループおよびランダムコイル等のそれぞれ、又はこれらの組み合わせを特徴的な部分構造として採用することもできる。特徴的な部分構造は複数の2次構造を含んでもよい。つまり、2次構造のそれぞれを1つの特徴的な部分構造とすることもでき、主鎖に沿って隣接した複数の2次構造をまとめて1つの特徴的な部分構造とすることもできる。また、核酸の螺旋構造に阻害剤が挟まって結合する反応を解析の対象にする場合には、上記で挙げた主鎖の他、当該螺旋構造の周期構造、又はこれらの組み合わせを特徴的な部分構造として採用することもできる。螺旋構造については、全体を1つの特徴的な部分構造としてもよいし、螺旋構造を所定の周期で分割して得られた部分的な構造のそれぞれを特徴的な部分構造としてもよい。
【0113】
「非中継座標」は、質点系の構造を記述しかつ中継座標から独立した座標である。非中継座標も、中継座標と同様に、質点系に含まれる一部の質点の座標そのもの、質点の座標を組み合わせて定義できる座標又はこれらの組み合わせを用いて設定される。質点座標の組み合わせによっては、非中継座標の構成の一部を中継座標で表現することもできる。ただし、非中継座標は、中継座標から独立した座標であるため、中継座標のみによって表現されることはない。例えば非中継座標は、一部の質点の座標が中継座標として設定されている場合には、N個の質点のうちその質点を除く質点座標とすることができる。また、非中継座標は、質点座標を組み合わせて定義できる座標のみが中継座標として設定されている場合には、そのすべてまたは一部を、質点座標そのもの或いは質点座標を組み合わせて定義できる座標で、かつ中継座標とは独立な座標とすることができる。
【0114】
中継座標の具体的な設定方法の2つの例について説明する。
【0115】
1つは、多原子系に含まれるN個の原子を特徴的な部分構造ごとに集団分けして原子の集団ごとに重心を抽出し、これらの重心の座標のみを中継座標として設定する方法である。この場合には、非中継座標は、その原子が属する部分構造の重心座標に対する相対座標でかつ独立な座標とすることが好ましい。したがって、非中継座標の個数は3N−M個となる。これは、多原子系の3N個の質点座標をM個の重心座標と3N個の相対座標に分離したとき、相対座標間には定義からM個の関係式が発生するため、つまり相対座標のうちM個は独立な座標ではないためである。この方法は、特徴的な部分構造同士の位置関係が大振幅集団運動による多原子系の構造変化に大きく影響することを考慮したものである。特徴的な部分構造自体の変形の規模が小さく、かつその変形の時間スケールが化学反応の時間スケールに比べて短いため、その変形が多原子系の構造変化に与える影響は小さい。
【0116】
もう1つは、生体高分子の主鎖を構成する原子のそれぞれの座標を中継座標として設定する方法である。この場合には、非中継座標は、その生体高分子の側鎖および溶媒和等のその他の原子の座標とすることが好ましい。したがって、非中継座標の個数は3N−M個となる。この方法は、生体高分子の主鎖の形状が大振幅集団運動による多原子系の構造変化に大きく影響することを考慮したものである。側鎖は主鎖に従属して移動することから、側鎖の変形が多原子系の構造変化に与える影響は小さい。この場合において、結合候補分子等の生体高分子以外の分子が多原子系に含まれている場合には、生体高分子以外の分子については、その分子の重心を中継座標として設定してもよい。この場合には、例えばその分子を構成する原子の座標は非中継座標として設定される。
【0117】
なお、後述する第1の実施形態は前者の設定方法を用いた場合に相当する。
【0118】
中継座標は、必要に応じて、すべての質点を自由な状態(質点系に固定しない状態)にした上で質点系に対して全構造緩和が行われた後に、設定されることが好ましい。「全構造緩和」とは、適当な開始状態から質点系のポテンシャルエネルギーが常に減少していくように、そのポテンシャルエネルギーの勾配がゼロとなる状態まで、質点の座標を移動させることを意味する。全構造緩和はいわゆる一般的な構造緩和と同義であるが、後述する「部分構造緩和」と区別するためにこのように称することにする。この全構造緩和によって見出される質点系の状態は、上記開始状態近傍のローカルミニマムの一つである。つまり、この全構造緩和は、ポテンシャル曲面上の上下を繰り返しながら安定点間を結ぶ反応経路およびミニマムポイント(最小点)を見出す作業ではない。この全構造緩和は、特に、構造緩和する対象を構成する原子のすべてを移動させる点で後述する「部分構造緩和」と区別される。この全構造緩和は、例えば、共役勾配法、最急降下法、逆ヘッセ法等既存の方法を用いて実施される。
【0119】
<座標抽出手段>
座標抽出手段は、中継座標の関数としての非中継座標の構成を求め、集団座標の関数としての中継座標の構成を求める。「中継座標の関数としての非中継座標の構成」とは、中継座標RSで表された非中継座標RFについての関数RF(RS)の具体的な内容を意味し、「集団座標の関数としての中継座標の構成」とは、集団座標q†で表された中継座標RSについての関数RS(q†)の具体的な内容を意味する。
【0120】
(中継座標の関数としての非中継座標の構成)
非中継座標の構成RF(RS)は、非中継座標RFを中継座標RSに従属させること、つまり非中継座標RFおよび中継座標RSの関係を規定する条件式を与えることにより得られる。このような条件式は、特に限定されないが、例えば断熱近似条件から得られる下記式30の条件式を用いることができる。なお、「断熱近似」とは、非中継座標RFが中継座標RSの変化に対し即座に追随することができるとする近似を意味する(M. Born, and J. R. Oppenheimer, Ann. Phys. 84 (1927) 457およびH. Haken, Synergetics 2rd Edition, Springer, 1978参照)。式30は、与えられたRSに対し、RFが、ポテンシャルエネルギーVの局所安定点、つまりポテンシャルエネルギーVのRFに関する勾配がゼロの点にいるという条件を表している。
【0121】
【数16】
【0122】
式30において、V(RS,RF)は、中継座標RSおよび非中継座標RFを用いて表された質点系のポテンシャルエネルギーを表す。V(RS,RF)は、質点座標xで表されたポテンシャルエネルギーV(x)に、中継座標RSが設定された段階で定まる質点座標xの構成x(RS,RF)を代入することにより得ることができる。例えば、質点座標xの構成x(RS,RF)は、中継座標RSが設定された段階で定まる中継座標RSの構成RS(x)および非中継座標RFの構成RF(x)に基づいて、これらの逆関数を求めることにより得られる。
【0123】
なお、すべての非中継座標RF1〜RF(3N−M)を中継座標RSに従属させるため、すべての非中継座標RF1〜RF(3N−M)でV(RS,RF)を偏微分することが必要である。
【0124】
非中継座標RFが中継座標RSに従属した質点系においては、下記式31のように非中継座標RFが中継座標RSとして表されることになり、質点系のポテンシャルエネルギーVは下記式32のように表される。以下、上記のような非中継座標RFが中継座標RSに従属した条件の下、中継座標RSのみの関数として表された質点系のポテンシャルエネルギーVeffを有効ポテンシャルエネルギーともいう。
【0125】
【数17】
【0126】
また、式31のように非中継座標RFが中継座標RSに従属する場合には、質点座標xは下記式33のように中継座標の関数として表すことができる。
【0127】
【数18】
【0128】
(集団座標の関数としての中継座標の構成)
中継座標の構成RS(q†)は、非中継座標RFが中継座標RSに従属するという条件の下、中継座標RSの変化に伴う非中継座標RFの変化が中継座標RSに与える影響を考慮しながら、中継座標RSに対して集団運動理論における基礎方程式を解くことにより求められる。中継座標の構成RS(q†)を求めるに際し重要な点は、「中継座標RSに対して集団運動理論における基礎方程式を解くこと」およびその際に「中継座標RSの変化に伴う非中継座標RFの変化が中継座標RSに与える影響(フィードバック)を考慮すること」である。
【0129】
以下、この2点について説明する。
【0130】
中継座標RSに対して集団運動理論における基礎方程式を解く理由は、前述したように、3N次元の質点座標xから集団座標q†を抽出するのではなく、M次元の中継座標RSから集団座標q†を抽出することにより、変数分離の計算を単純化するためである。
【0131】
SCC法における最も基本的な基礎方程式は下記式34から式36である(非特許文献10)。
【0132】
【数19】
【0133】
miは質点系におけるi番目の中継座標に関する質量を表す。miは、i番目の中継座標が複数の質点の重心座標である場合には、それらの質点の質量の総量となる。中継座標が、質点座標を組み合わせた座標であって重心座標ではない場合、単純に質量の総量にはならない。この場合、通常の解析力学の処方で中継座標に対応する一般運動量を求め、ハミルトニアンを書き下し、そのハミルトニアンの一般運動量が関与する項の係数として、中継座標に関する質量を求める。また、φiμ(RS)、λμ(RS)およびρμ(RS)は式35および式36を満たすRSの関数であり、φiμ(RS)は(i,μ)番目、λμ(RS)およびρμ(RS)はμ番目の成分である。また、V,i(RS)およびV,ij(RS)はそれぞれ下記式37および式38で定義される。
【0134】
【数20】
【0135】
したがって、一般的には式34から式36に従い、M個の座標(RS1,RS2,・・・,RSM)で張られるM次元の超空間の中にK個の集団座標(q1,q2,・・・,qK)でパラメトライズ(媒介変数化)されるK次元の超平面(以下、単に平面という)を構成することで、関数RSi(qμ)が求められる。
【0136】
また、基礎方程式としては非特許文献10における近似計算のための基礎方程式を援用することもできる。このとき、例えば自由度が1(K=1)の集団座標q†の関数としての中継座標の構成RS(q1)を求める場合には、中継座標RSを用いた集団運動理論における基礎方程式は下記式39および式40のように表される。
【0137】
【数21】
【0138】
φi(RS)は、式40を満たす関数(固有ベクトル)のi番目の成分を表し、
Λ(RS)は、式40を満たす関数(固有値)を表す。
【0139】
一方、従来の非特許文献10におけるSCC法では、近似的に計算が行われており、変数分離を高精度に行うことができない場合がある。そこで、本発明者は、中継座標RSに対する変数分離をより正確に計算することができる新たな近似法に基づき基礎方程式を導出し、この基礎方程式に基づく新たなSCC法(SCC2法)の理論を構築した。
【0140】
例えばSCC2法に基づいて、自由度が1(K=1)の集団座標q†の関数としての中継座標の構成RS(q1)を求める場合には、中継座標RSを用いた集団運動理論における基礎方程式は下記式41から式43のように表される。
【0141】
【数22】
【0142】
φi(RS,λ)およびκ(RS,λ)は式43を満たす関数であり、φi(RS,λ)はそのi番目の成分を表し、
λは、補助座標を表す。
【0143】
最終的に、式39または式41の解として、式44で表されるような中継座標RSの構成が求められる。
【0144】
【数23】
【0145】
また、上記ではK=1の場合における基礎方程式について説明したが、本発明はK=1の場合に限られない。具体的には、K≧2の場合におけるSCC法の基礎方程式への式39および式40の拡張は、非特許文献10に記載されているように、フロベニウスの定理を用いて行うことができる。一方、K≧2の場合におけるSCC2法の基礎方程式への式41から式43の拡張も、同様の方法で自明に行うことができる。
【0146】
このように、中継座標という階層化された変数(つまり質点系の構造を二次的に記述するより低次元の変数)を導入することにより、3N次元からK次元への縮約の問題がM次元からK次元への縮約の問題に簡略化されるため、変数分離の計算が容易となる。
【0147】
しかしながら、本発明では、質点系の構造を正確に計算するために、フィードバックを考慮しながら変数分離を実施しなければならない。フィードバックを考慮しないということは、存在する質点(または非中継座標)を存在しないものとみなして計算することと同義であり、この影響を考慮せずに変数分離を行うと計算内に矛盾が生じて計算結果が破錠する問題が起こるためである。そこで、フィードバックを考慮するために、既存の集団運動理論に新たな理論に基づいた修正を加える。この新たな理論とは、非中継座標RFが中継座標RSに従属するという条件の下、集団運動理論の基礎方程式におけるポテンシャルエネルギーの座標に関する微分係数を正確に計算し、さらに、中継座標RS(q†)を微小更新(RS(q†)→RS(q†+Δq†))させるごとに非中継座標RFに対して部分構造緩和を行うというものである。
【0148】
従来の集団運動理論では、ポテンシャルエネルギーV(x)の質点座標xに関する微分係数に対する所定の固有ベクトルφの方向に、質点座標xを微小更新(x(q†)→x(q†+Δq†))することが繰り返されるのみである。
【0149】
一方本発明では、ポテンシャルエネルギーVの座標に関する微分係数について、その計算をする際、下記式45から式47のフィードバック方程式を適用することができる。
【0150】
【数24】
【0151】
式46において、第3項の(i←→j)は、第2項におけるiおよびjの文字を相互に入れ替えた項を表す。
【0152】
式47において、第3項の(i←→k)は、第2項におけるiおよびkの文字を相互に入れ替えた項を表し、第4項の(j←→k)は、第2項におけるjおよびkの文字を相互に入れ替えた項を表し、第6項の(i←→k)は、第5項におけるiおよびkの文字を相互に入れ替えた項を表し、第7項の(j←→k)は、第5項におけるjおよびkの文字を相互に入れ替えた項を表す。
【0153】
さらに、式45から式47において下記式48を用いた。
【0154】
【数25】
【0155】
Kαβ−1(Rs)はKαβ(Rs)の逆行列を表し、
Kαβ(Rs)およびJαi(Rs)はそれぞれ式49および式50で定義される。
【0156】
【数26】
【0157】
式45から式47それぞれにおいて、右辺の第1項以外の項が、中継座標RSの変化に伴う非中継座標RFの変化が中継座標RSに与える影響を表している。このようなフィードバック方程式の構成は従来知られておらず、本発明者が初めて見出したものである。フィードバック方程式は、1階微分、2階微分および3階微分の3つの式を含むが、集団運動理論における基礎方程式の内容に応じて必要なものを使用すればよい。例えば、集団運動理論における基礎方程式として、式39および式40を採用した場合には、フィードバック方程式は式46の2階微分の式のみで充分である。
【0158】
そして、非中継座標RFに対する部分構造緩和は、ポテンシャルエネルギーVの座標に関する微分係数であってフィードバック方程式で規定されるものに対する固有ベクトルφの方向に、中継座標RS(q†)を微小更新(RS(q†)→RS(q†+Δq†))させ、更新された中継座標RS(q†+Δq†)の下で、非中継座標RFが中継座標RSに従属するという条件に従って構造緩和を行うことにより行われる。「中継座標RS(q†+Δq†)の下で」とは、更新された中継座標RS(q†+Δq†)を質点系に固定した状態でという意味である。また、「非中継座標RFが中継座標RSに従属するという条件に従って構造緩和を行う」とは、質点系のポテンシャルエネルギーが常に減少していくように、そのポテンシャルエネルギーの勾配がゼロとなる状態まで、中継座標RSに従属した状態を維持しながら非中継座標RFを移動させることを意味する。上記のように中継座標RSを質点系に固定した状態で行う構造緩和を「部分構造緩和」という。前述したように、この部分構造緩和は、中継座標RSを質点系に固定した状態で構造緩和を行う点で全構造緩和と異なる。ただし、部分構造緩和においても、構造緩和を行う手法としては、全構造緩和と同様に、例えば、共役勾配法、最急降下法、逆ヘッセ法等既存の方法を用いることができる。なお、更新する方向を決める上記固有ベクトルφは、最小固有値に対応するものであることが好ましい。本発明は、質点系の大振幅集団運動を対象としているためである。
【0159】
K=1である場合、部分構造緩和における非中継座標RFの開始状態は、RFα(RS(q1))でもよいが、下記式51で与えられる微小量だけ変化させた状態、つまり下記式52で表される状態RFα†であることが好ましい。
【0160】
【数27】
【0161】
式52で表される座標RFα†は、微小更新が無限小ならば非中継座標RF(RS(q1+Δq1))と一致する。しかしながら、実際の計算では微小更新は有限である。部分構造緩和における非中継座標RFの開始状態をRFα†とすることにより、部分構造緩和を行う上での安全性を高めることができる。
【0162】
(座標抽出手段における処理手順)
座標抽出手段における一連の具体的な処理手順について説明する。
【0163】
まず、座標抽出手段は、中継座標Rsおよび非中継座標RFの関数として表現された質点系のポテンシャルエネルギーV(RS,RF)を求める第1の工程を実施する。
【0164】
次に、座標抽出手段は、非中継座標RFを中継座標RSに従属させる第2の工程を実施する。中継座標RSの関数としての非中継座標RFの構成が得られる。ここでは、例として式30の断熱近似条件に従って非中継座標RFを中継座標RSに従属させた場合を考える。
【0165】
次に、座標抽出手段は、座標設定手段によって設定された状態における中継座標RSおよび非中継座標RFの下で、例えば式46のフィードバック方程式に従って、中継座標RSに関するポテンシャルエネルギーVの微分係数を求める第3の工程を実施する。
【0166】
次に、座標抽出手段は、ポテンシャルエネルギーVの上記微分係数の下で、例えば式39および式40のSCC法におけるK=1の場合の基礎方程式に従って、集団座標qに関する中継座標RS(q)の微分係数を求める第4の工程を実施する。
【0167】
次に、座標抽出手段は、中継座標RS(q)の上記微分係数の下で集団座標qを微小量Δqだけ更新することにより、更新された中継座標RS(q+Δq)を求める第5の工程を実施する。
【0168】
次に、座標抽出手段は、更新された中継座標RS(q+Δq)の下で、上記断熱近似条件に従って中継座標RSに従属した非中継座標RFについて部分構造緩和を行う第6の工程を実施する。この際、式52で表される状態RFα†を部分構造緩和の開始状態とする場合には、第4の工程で得られた中継座標RS(q)の上記微分係数から式51に従って、集団座標qに関する非中継座標RF(q)の微分係数を求め、この微分係数の下で非中継座標を微小量Δqだけ更新させておく。
【0169】
そして、座標抽出手段は、中継座標RSおよび非中継座標RFを交互に更新しながら、第3の工程から第6の工程までを例えば次の極小点に至るまで順次繰り返すことにより、集団座標qの関数としての中継座標RSの構成を求める。つまり、2回目以降の第3の工程においては、直前の第6の工程で行われた非中継座標RFの部分構造緩和後の状態における中継座標RSおよび非中継座標RFの下で、上記第3の工程と同様に、中継座標RSに関するポテンシャルエネルギーVの微分係数を求めることになる。極小点を判定する方法は、特に限定されないが、例えばポテンシャルエネルギーVの上記微分係数がゼロ相当になったか否かを判定する方法を挙げることができる。
【0170】
<逆変換手段>
逆変換手段は、集団座標q†についての運動方程式の解として得られる時間tの関数としての集団座標q†(t)、式44で表される中継座標RSの上記構成RS(q)および式31で表される非中継座標RFの上記構成RF(RS)に基づいて、質点座標xの時間発展を予測する。
【0171】
具体的には以下の通りである。図4は、本発明における変数変換および逆変換の概念を示す図である。図4では、例として質点座標xの一部の座標そのものを中継座標RSとし、この中継座標RSから自由度が1の集団座標q†=q1を抽出した場合について示している。
【0172】
図4aは、質点座標xに基づいて中継座標RSおよび非中継座標RFを設定する様子を示している。図4aでは、M個の質点座標が中継座標RSとして設定され、それ以外の3N−M個の質点座標が非中継座標RFとして設定されている。中継座標RSおよび非中継座標RFの設定は、前述したとおり、座標設定手段において行われる。図4bは、非中継座標RFが中継座標RSの関数として与えられることを示している。図4cは、中継座標RSが集団座標q1の関数として与えられることを示している。非中継座標RFの構成および中継座標RSの構成は、前述したとおり、座標抽出手段において求められる。
【0173】
図4dは、集団座標q1が時間tの関数として与えられることを示している。集団座標q1についての運動方程式を解くことにより求められる。集団座標q1についての運動方程式は、質点座標xについての運動方程式に、中継座標RSが設定された段階で定まる質点座標xの構成x(RS,RF)を代入して、中継座標RSおよび非中継座標RFについての運動方程式を求め、この運動方程式に式31および式44を代入することにより求めることができる。
【0174】
以上の図4aからdまでの関係を用いることにより、すべての質点座標xを時間tの関数として求めることができる(図4e)。つまり、逆変換手段は、時間tの関数としての集団座標q1を求めた後、中継座標RSを介した質点座標xから集団座標q1への変数変換に対してこれを逆変換することにより、時間tの関数としての質点座標x(t)の構成を求めるものである。この関数x(t)は、時間変化に対する質点座標の挙動を表すものであるから、この関数x(t)こそが求めるべき反応経路に相当するものである。したがって、この関数x(t)の時間変化を見ることにより質点座標の時間発展を予測することが可能となる。
【0175】
「本発明の第1の実施形態」
図5は、本実施形態のシミュレーション装置10の構成を概略的に示すブロック図である。また、図6Aおよび図6Bは、本実施形態のシミュレーション方法における計算手順を示すフローチャート図である。なお、本実施形態では、生体高分子およびこれに対する結合候補分子を含むN個の原子から構成される多原子系において、多原子系の構造のうち特徴的な部分構造ごとに抽出された重心の座標をM個の中継座標として設定し、かつ、集団座標の自由度が1(K=1)である場合について説明する。
【0176】
本実施形態のシミュレーション装置10は、生体高分子および結合候補分子から構成される複合体の挙動を予測するシミュレーション装置であって、図5に示されるように、解析対象となる多原子系を規定するのに必要なデータ(入力データ)を入力するための入力手段12と、装置の各部を制御する制御手段14と、座標設定手段16と、座標抽出手段18と、逆変換手段20と、解析結果を表示する表示手段22とを備える。
【0177】
そして、本実施形態のシミュレーション方法は、シミュレーション装置10により実行されるシミュレーション方法であって、入力手段12によって入力データを入力し、座標設定手段16によってM個の中継座標RSおよび3N−M個の非中継座標RFを設定させ、座標抽出手段18によって、K=1の集団座標q1の関数としての中継座標RSの構成および中継座標RSの関数としての非中継座標RFの構成を求めさせ、逆変換手段20によって時間tの関数としての原子の座標x(t)を求めさせ、表示手段22によって解析結果を表示させるものである。
【0178】
また、本実施形態に係るシミュレーションプログラムは、上記に記載のシミュレーション方法をコンピュータに実行させることを特徴とするものである。
【0179】
また、本実施形態に係るコンピュータ読み取り可能な記録媒体は、上記に記載のシミュレーションプログラムを記録したものである。
【0180】
<入力手段>
入力手段12は、ユーザーにより入力データが入力される部分である。入力された入力データは制御手段14へ出力される。入力データの入力方法は、特に限定されず、例えばシミュレーション装置10をユーザーが操作しながら手動で行う方法、所定の記録媒体から読み込む方法等が挙げられる。
【0181】
<制御手段>
制御手段14は、各部とのデータのやり取り等の処理を制御する部分である。入力手段12から制御手段14へ入力された入力データは座標設定手段16に出力される。入力データが座標設定手段16へ出力されることにより、当該入力データに基づく多原子系の挙動の計算が開始される。制御手段14には、各部とのデータのやり取りの中で生成される途中のおよび最終的な計算結果を記録するため、メモリ等の記憶媒体が含まれる。
【0182】
また、制御手段14は、逆変換手段20で得られた原子の座標x(t)に関するデータに基づいて、その結果を図やグラフを用いて表示するよう表示手段22を制御する。
【0183】
<座標設定手段>
座標設定手段16は、制御手段14から座標設定手段16へ入力された入力データに基づいて、複合体の構造を最適化し、その後M個の中継座標RSおよび3N−M個の非中継座標RFを設定するものである。得られた中継座標RSおよび非中継座標RFに関するデータは、制御手段14を介して座標抽出手段18に出力される。
【0184】
<座標抽出手段>
座標抽出手段18は、座標設定手段16で得られた中継座標RSおよび非中継座標RFに関するデータに基づいて、多原子系全体のポテンシャルエネルギーV(RS,RF)を求め、このV(RS,RF)に基づいて集団座標q1の関数としての中継座標RSの構成および中継座標RSの関数としての非中継座標RFの構成を求めるものである。得られた中継座標RSの構成および非中継座標RFの構成に関するデータは、制御手段14を介して逆変換手段20に出力される。
【0185】
<逆変換手段>
逆変換手段20は、座標抽出手段18で得られた中継座標RSの構成および非中継座標RFの構成に関するデータに基づいて、時間tの関数としての原子の座標x(t)を求めるものである。得られた原子の座標x(t)に関するデータは制御手段14に出力される。
【0186】
<表示手段>
表示手段22は、制御手段14の指示に従い、例えば多原子系のポテンシャルエネルギーVをグラフで表示したり、複合体の結合過程における構造および/または最終状態の構造を画像として表示したり、複合体の結合過程を動画として表示したりするものである。表示方法としては特に限定されず、2D表示、3D表示等公知の方法を用いることが可能である。
【0187】
以下、図6Aおよび図6Bを参照しながら、本実施形態の装置10を用いたシミュレーション方法の手順を詳細に説明する。
【0188】
<STEP1>
まず、入力手段12から入力データがシミュレーション装置10に入力される(図6AのST1)。入力データの具体的な内容は、生体高分子およびこれに対する結合候補分子を含むN個の原子から構成される多原子系を構成する原子の種類、個数、初期配置の情報となる座標xおよび座標xからポテンシャルエネルギーV(x)を構成するルール等である。
【0189】
<STEP2>
次に、複合体を含む多原子系の構造が最適化されて、初期状態における構造が求められる(図6AのST2)。
【0190】
具体的には以下の通りである。まず、座標設定手段16は、生体高分子および結合候補分子に対してそれぞれ独立に全構造緩和を実施して、生体高分子および結合候補分子の構造をそれぞれ最適化する。そして、それぞれ構造が最適化された生体高分子および結合候補分子は、互いに相互作用を開始する程度(例えば、結合候補分子が生体高分子の表面にわずかに接触する程度)の距離に配置される。この状態が、生体高分子および結合候補分子からなる複合体の初期状態となる。
【0191】
<STEP3>
次に、STEP2で得られた初期状態における構造に基づいて、中継座標RSおよび非中継座標RFが設定される(図6AのST3)。
【0192】
座標設定手段16は、生体高分子の特徴的な部分構造ごとに(M−3)/3個の部分構造に分け、結合候補分子全体を1個の部分構造として取り扱い、M/3個の部分構造のそれぞれから重心を抽出する。そして、このM/3個の重心の座標がM個の中継座標RSとして設定され、M/3個の各部分構造に含まれる原子の各重心からの相対座標でかつ独立な座標が3N−M個の非中継座標RFとして設定される。本実施形態では、非中継座標RFとして、その原子の属する部分構造の重心座標からの相対座標が設定されるとする。
【0193】
<STEP4>
次に、中継座標RSおよび非中継座標RFで表された多原子系のポテンシャルエネルギーV(RS,RF)が求められる(図6AのST4)。
【0194】
図7は、本実施形態における中継座標RS、非中継座標RFおよび原子の座標xの関係を示す概略図である。なお、図7において、m(mは1≦m≦M/3を満たす整数)番目の特徴的な部分構造に含まれる原子数をNmと表した。つまり、N1+N2+・・・+Nm+・・・+NM/3=Nである。本発明では、Mは3次元空間における原子(質点)の座標の個数であるため、通常M/3は整数となる。また、本実施形態においては、中継座標RS、非中継座標RFおよび原子の座標xは、図7から具体的に下記式53および式54で表されるような関係を有していることがわかる。
【0195】
【数28】
【0196】
式53は、m番目の中継座標RS、つまりm番目の特徴的な部分構造の重心座標と原子の座標xとの関係を示した式である。式54は、m番目の特徴的な部分構造に属するn番目の原子の座標x、中継座標RSおよび非中継座標RFの関係を示した式である。式53および式54において、nは、特徴的な部分構造のそれぞれに含まれる原子の通し番号である。したがって、m番目の特徴的な部分構造についてのnは、1からNmまでの整数値をとる。また、3N1+3N2+・・・+3Nm−1はm=1のときゼロとする。
【0197】
したがって、ポテンシャルエネルギーV(RS,RF)は、ポテンシャルエネルギーV(x)に対して式54を適用することにより得ることができる。
【0198】
<STEP5>
次に、ポテンシャルエネルギーV(RS,RF)を用いた式30の断熱近似条件に従って、非中継座標RFが中継座標RSに従属させられる(図6AのST5)。本実施形態では、すべての非中継座標RF1〜RF3NでV(RS,RF)を偏微分することになる。
【0199】
<STEP6>
次に、現状の中継座標RSおよび非中継座標RFの下で、式46のフィードバック方程式に従い、中継座標RSに関するポテンシャルエネルギーVの2階微分係数が求められる(図6AのST6)。
【0200】
<STEP7>
次に、STEP6で得られた2階微分係数の下で、式39および式40のSCC法の基礎方程式に従い、集団座標q1に関する中継座標RSの微分係数dRS/dq1が求められる(図6AのST7)。
【0201】
<STEP8>
次に、STEP7で得られた微分係数dRS/dq1の下で、式51に従い集団座標q1に関する非中継座標RFの微分係数dRF/dq1が求められる(図6BのST8)。
【0202】
<STEP9>
次に、STEP7で得られた微分係数dRS/dq1およびSTEP8で得られた微分係数dRF/dq1の下で、集団座標q1が微小量Δq1だけ更新される(図6BのST9)。
【0203】
<STEP10>
次に、微小量Δq1だけ更新された中継座標RS(q1+Δq1)および更新された非中継座標RF(q1+Δq1)が求められる(図6BのST10)。
【0204】
例えば、STEP2における初期状態を集団座標q1の初期条件と設定した(つまり、q1=0のときに多原子系が上記初期状態をとるものとした)場合には、中継座標RS(q1+Δq1)および非中継座標RF(q1+Δq1)は下記式55のように表せられる。
【0205】
【数29】
【0206】
また、STEP13の判断結果に応じてSTEP6からSTEP12までを順次繰り返し実施する場合において、2回目以降のSTEP10で順次下記式56のように表せられるRSi(2Δq1)、RSi(3Δq1)、・・・が得られることになる。ただし、微分係数dRS/dq1は各回で同じ値になるとは限らない。上記の作業を繰り返すことにより、最終的に集団座標q1の関数としての中継座標RSの構成RS(q1)が求められる。
【0207】
【数30】
【0208】
<STEP11>
STEP11では、STEP10で更新された中継座標RS(q1+Δq1)の下で、STEP10で更新された非中継座標RF(q1+Δq1)を開始状態として非中継座標RFについて部分構造緩和が行われる(図6BのST11)。
【0209】
<STEP12>
次にSTEP12において、部分構造緩和後の状態における多原子系のポテンシャルエネルギーVが計算される(図6BのST12)。なお、STEP13の判断結果に応じて、繰り返しの各回でポテンシャルエネルギーVが計算される。ポテンシャルエネルギーの数値計算は、有効ポテンシャルエネルギーVeffに基づいた計算でもよい。
【0210】
<STEP13>
次にSTEP13において、STEP12で計算されたポテンシャルエネルギーVが極小点に至ったか否かが判断される(図6BのST13)。生体高分子および結合候補分子の化学反応を考えた場合、化学反応の終了点はポテンシャルエネルギーVが次に極小となるとき、つまりポテンシャル曲面上の次の安定点である。そこで、ポテンシャルエネルギーVの値を基準に化学反応の終了点が判断される。ポテンシャルエネルギーVが次の極小点に至ったと判断されればSTEP14へ進み、そうでなければ再度STEP6からSTEP12までが繰り返される。また、ポテンシャルエネルギーVが極小点に至ったと判断された場合であっても、必要に応じて、ポテンシャルエネルギーVがさらに次の極小点に至るまで解析を進めることも可能である。
【0211】
<STEP14>
STEP13においてポテンシャルエネルギーVが極小点に至ったと判断された場合、STEP14において、そのポテンシャルエネルギーVにおける多原子系の構造が、求めるべき最終状態の構造となる(図6BのST14)。
【0212】
そして、最終状態の構造が得られるとともに、式44で表される中継座標RSの構成RS(q1)および式31で表される非中継座標RFの構成RF(RS)が求められる。その後、時間tの関数としての座標x(t)が求められる。
【0213】
<作用効果>
図8は、本実施形態のシミュレーション方法によって得られた、所定のタンパク質2およびこれに対する生理活性物質6(薬物)の複合体の反応経路を示すグラフである。また、図9は、図8のグラフにおける複合体の結合過程を示す概略図である。図9aは複合体の結合過程における初期状態を示し、図9bは複合体の結合過程における最終状態を示す。
【0214】
図8のグラフにおいて、横軸は所定のタンパク質2の残基2aと残基2bの距離であり、縦軸は薬物6および上記タンパク質2の残基2cの距離である。残基2aおよび残基2bは、タンパク質2のポケット4の開口にある残基であり、残基2cは、タンパク質2および薬物6が結合した時に薬物6に対向する位置にある残基である(図9)。また、グラフ中のAおよびCが、例えば図2における安定点AおよびCにそれぞれ対応し、グラフ中のBが、例えば図2におけるエネルギー障壁Bに対応している。グラフ中のプロットは時間軸上で等間隔な点を示しており、このグラフからエネルギー障壁Bに至るまでは反応がゆっくり進み、エネルギー障壁Bを超えた後は反応が速やかに進む様子がわかる。本実施形態のシミュレーション方法により、2000個程度の原子を含む多原子系においてこの反応経路を求めるのに要した時間はおよそ30分である。
【0215】
また、図8のグラフ中の「X線解析」の位置は、上記タンパク質2および上記薬物6についてX線回折測定による結晶構造解析の結果実際に得られた実測値を示す。これにより、本発明は非経験的にシミュレーションを行っているにもかかわらず、非常に高い精度でこの複合体の結合状態を計算できていることがわかる。
【0216】
以上のように、本実施形態に係るシミュレーション装置は、上記のような変数設定手段と、変数抽出手段と、逆変換手段とを備え、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0217】
本実施形態に係るシミュレーション方法は、上記シミュレーション装置によって実行される方法であって、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0218】
本実施形態に係るプログラムおよび記録媒体は、上記シミュレーション方法を実行させることが可能であるから、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0219】
「本発明の第2の実施形態」
次に第2の実施形態について説明する。本実施形態は、非特許文献10におけるSCC法を援用した基礎方程式(式39および式40)またはSCC2法の基礎方程式(式41から式43)の解を求める上で、これらの基礎方程式とそれぞれ同値の方程式を解く点で、第1の実施形態と異なる。したがって、第1の実施形態と同様な構成要素の詳細な説明は必要のない限り省略する。
【0220】
本実施形態のシミュレーション装置10も、第1の実施形態と同様に、入力手段12と、制御手段14と、座標設定手段16と、座標抽出手段18と、逆変換手段20と、解析結果を表示する表示手段22とを備える。
【0221】
そして、本実施形態では座標抽出手段18は、式39および式40と同じ解を得るために、SCC法における基礎方程式として下記式57を使用する。
【0222】
【数31】
【0223】
Yは、下記式58で定義されるMK+M+K次元のベクトルである。つまり、Yは、iおよびμのそれぞれを走らせた各要素、m11/2RS1,・・・,mM1/2RSMと、φ11,・・・,φM1,・・・,φ1K,・・・,φMKと、Λ1,・・・,ΛKとを有する。なおφiμおよびΛμは補助座標を表す。
【0224】
【数32】
【0225】
vμは、下記式59の非斉次線形方程式の解ベクトルである。なお、vμはMK+M+K次元のベクトルである。
【0226】
【数33】
【0227】
式59におけるCおよびsμはそれぞれ式60および式61で定義される。なお、Cは(MK+M+K)×(MK+M+K)の正方行列である。Cは一見3×3の行列に見えるが、各成分に更に行列を有する。つまり、3×3の行列として表現されたCのうち、1行1列の成分はMK×M行列(つまり行数がMKで列数がM)であり、1行2列の成分はMK×MK行列であり、1行3列の成分はMK×K行列であり、2行2列の成分はK×MK行列であり、3行1列の成分はM×M行列であり、その他の成分はゼロ行列である。各成分の要素は、iあるいはμを行の足としてi=1〜Mとμ=1〜Kで走らせ、jあるいはνを列の足としてj=1〜Mとν=1〜Kで走らせて構成する。また、sμはMK+M+K次元のベクトルである。sμも、Cと同様に一見3次元のベクトルに見えるが、各成分に更にベクトルを有する。つまり、sμの3次元ベクトルとしての第3成分はM次元ベクトルであり、その他の成分はゼロベクトルである。第3成分の要素は、iをi=1〜Mで走らせて構成する。
【0228】
【数34】
【0229】
V,ijk(RS)は式62で定義される。
【0230】
【数35】
【0231】
上記式57はKが一般化された場合の基礎方程式である。そして、K=1の場合には、式57は下記式63と同値となる。また、上記非特許文献10に記載されているように、フロベニウスの定理を用いてK=1の場合の式63から、一般化されたKについての式57を導出することも可能である。
【0232】
【数36】
【0233】
式63において(χ1,・・・,χM,ψ1,・・・,ψM,κ)は、式64で表される非斉次線形方程式の2M+1次元の解ベクトルである。
【0234】
【数37】
【0235】
なお式64において、行列と固有ベクトルとの積を取り3つの方程式で表現したときにjおよびkについて1からMの総和をとる。
【0236】
式63および式64が式39および式40と同値であることは以下のようにして説明することができる。式38、式62および式63に注意して式64を(φi1,Λ1)に関してq1で積分すると式65が得られる。
【0237】
【数38】
【0238】
さらに、式64に含まれる恒等式χi=φi1を式63の第1式に代入すると式66が得られる。
【0239】
【数39】
【0240】
ここで、式65の第1式および式66におけるφi1およびΛ1をそれぞれφi(RS)およびΛ(RS)とおく(つまり、φi1およびΛ1をRSの関数として取り扱う)と、式39および式40が再現される。
【0241】
また、座標抽出手段18は、式41から式43と同じ解を得る場合には、SCC2法における基礎方程式として下記式67を使用する。
【0242】
【数40】
【0243】
Zは、下記式68で定義されるMK+M+2K次元のベクトルである。つまり、Zは、iおよびμのそれぞれを走らせた各要素、m11/2RS1,・・・,mM1/2RSMと、φ11,・・・,φM1,・・・,φ1K,・・・,φMKと、λ1,・・・,λKと、ρ1,・・・,ρKとを有する。なお、λμおよびρμは、φiμと同様に補助座標を表す。また、cμνは、μごとに定義される下記式69で表される値をそれぞれ最小にするように一意に決められる定数である。また、wμは、下記式70で定義された(MK+M+K)×(MK+M+2K)次元の退化行列DのK次元の特異値空間を張るK個のMK+M+2K次元の単位ベクトルである。Dは一見3×4の行列に見えるが、各成分に更に行列を有する。つまり、3×4の行列として表現されたDのうち、1行1列の成分はMK×M行列であり、1行2列の成分はMK×MK行列であり、1行3列の成分はMK×K行列であり、2行1列の成分はM×M行列であり、2行2列の成分はM×MK行列であり、2行4列の成分はM×K行列であり、3行2列の成分はK×MK行列であり、その他の成分はゼロ行列である。各成分の要素は、iあるいはμを行の足としてi=1〜Mとμ=1〜Kで走らせ、jあるいはνを列の足としてj=1〜Mとν=1〜Kで走らせて構成する。
【0244】
【数41】
【0245】
上記式67はKが一般化された場合の基礎方程式である。そして、K=1の場合には、式67は下記式71と同値となる。また、上記非特許文献10に記載されているように、フロベニウスの定理を用いてK=1の場合の式71から、一般化されたKについての式67を導出することも可能である。
【0246】
【数42】
【0247】
式71において(χ1,・・・,χM,ψ1,・・・,ψM,κ,σ)は、式72で表される(2M+1)×(2M+2)の退化行列による斉次線形方程式の2M+2次元の解ベクトル(不定解)である。
【0248】
【数43】
【0249】
なお式72において、行列と固有ベクトルとの積を取り3つの方程式で表現したときにjおよびkについて1からMの総和をとる。つまり、(χ1,・・・,χM,ψ1,・・・,ψM,κ,σ)は、式72の退化行列の1次元の特異値空間に属する単位ベクトルuを用いて、式73で表される。ここで、cは、式74で表せられる値を最小にするように一意に決められる定数である。
【0250】
【数44】
【0251】
式71および式72が式41から式43と同値であることは以下のようにして説明することができる。式38、式62、式71および式73に注意して式72を(φi1,λ1,ρ1)に関してq1で積分すると式75が得られる。
【0252】
【数45】
【0253】
式75の第1式および第2式からφi1およびρ1を消去すると式76が得られる。そして、RSiとλ1をq1の関数と見なして式76をq1で微分し、式71の第1式および第3式を用いると式77が得られる。
【0254】
【数46】
【0255】
ここで、式71の第1式および第3式並びに式77におけるλ1、χi、κおよびVをそれぞれλ、φi(RS,λ)、κ(RS,λ)およびVeff(RS)とおくと、式41から式43が再現される。
【0256】
本実施形態では、例えば式63または式71から分かるように、基礎方程式を解く上で独立した座標として取り扱うq1の関数としての変数(中継座標RSおよび補助座標の全体であり、以下独立した変数という。)の数が増えている。具体的には、式39および式40では独立した変数は、M個のRSi(i=1からM)のみであったのに対し、式39および式40と同値の式63では独立した変数は、M個のRSi(i=1からM)と、M個のφi1(i=1からM)と、Λ1とから構成され、その総数は2M+1個である。一方、式41から式43では独立した変数は、M個のRSi(i=1からM)およびλであったのに対し、式41から式43と同値の式71では独立した変数は、M個のRSi(i=1からM)と、M個のφi1(i=1からM)と、λ1と、ρ1とから構成され、その総数は2M+2個である。
【0257】
このように独立した変数(具体的には補助座標)の数を増やして方程式を同値変形することにより得られる利点は以下の通りである。
【0258】
例えば式40からφi(RS)を求める場合、式40がφi(RS)の符号反転に対して不変であることを反映し、φi(RS)の符号に任意性が残る。またこのことは、式43からφi(RS,λ)およびκ(RS,λ)を求める場合にも同じである。この符号の任意性は、系の状態が反応経路を逆走する経路もまた反応経路を与えることを表す。一度逆走しても再度逆走することにより順方向に戻れば、最終的に反応経路を見出すことができることを考慮すれば、この符号の任意性はむしろ自然な現象を表現していると言える。しかしながら、この符号の任意性は、上記逆走によって反応経路を効率よく構成することが難しくなるという問題を引き起こす場合がある。さらに、式40および式43は実用上差分化して解かれるが、その差分化の程度Δq1が有限値であるため、逆走が生じても系の状態が元の反応経路上に戻らない場合があり、このような逆走は誤った反応経路が構成されてしまうきっかけとなるという問題もある。
【0259】
そこで、本実施形態では、独立した変数の数を増やすことによりφiμ(RS)の符号の任意性を解消している。その結果、上記逆走の問題が解決され、より効率よくかつ正確に反応経路を構成することが可能となる。なお、独立した変数の数が2倍程度に増えたとしても、計算処理増加による負担はほとんど生じない。
【0260】
図10は、符号の任意性が解消された計算の場合と符号の任意性が残ったままの計算の場合との比較結果を示すグラフである。具体的には、符号24のグラフは、式71から式74を使用して、メリチンというタンパク質の構造変化ごとのポテンシャルエネルギーを計算した結果を示し、符号26のグラフは、式41から式43を使用して、同じタンパク質の構造変化ごとのポテンシャルエネルギーを計算した結果を示す。なお、横軸は、構造変化を見やすくするために導入されたメリチンの屈曲角θであり、縦軸は系のポテンシャルエネルギーである。θ=140degreeの状態が初期状態に相当する。
【0261】
符号の任意性が残ったままの計算では、θ=95degree程度でエネルギーが急上昇しており、この点で一度逆走してしまった。一方、符号の任意性が解消された計算ではそのような現象は見られず、θ=80degree程度でエネルギーの極大点に至った。グラフ24におけるθ=60degreeの離れた点は、上記極大点における状態から系に全体構造緩和を適用して得られたもう1つの局所安定点である。
【0262】
以上のように、本実施形態に係るシミュレーション装置およびシミュレーション方法も、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものであるから、第1の実施形態と同様の効果を奏する。
【0263】
<設計変更>
上記各実施形態の基礎方程式において、ポテンシャルエネルギーVeff(RS)の3次微分の項が登場する(例えば、式43、式60、式64、式70および式72)。この3次微分の計算処理に要する時間はオーダーでM3(Mは中継座標の個数)に比例するため、その2次微分の計算処理に要する時間がオーダーでM2に比例することを考慮すると、この3次微分の計算はシミュレーションを行う上で非常に大きな負担である。したがって、例えば式64におけるV,ijk(RS)・φk1について、V,ijk(RS)を先に計算しその計算値とφk1の積をとると多くの時間を要してしまう。そこで、本発明者は下記式78を開発し、2つの2次微分との差分をとることで直接V,ijk(RS)・φk1を計算できるようにした。式78を適用することにより、オーダーでM3程度となる計算時間をM2程度にまで短縮することが可能となる。さらに、式78の右辺の2次微分に対して式46のフィードバック方程式を適用することが好ましい。
【0264】
【数47】
【0265】
φiは任意のベクトルのi番目の成分を表し、nは式79で定義される。
【0266】
【数48】
【0267】
(SCC法の基礎方程式について)
最後に、SCC法における最も基本的な基礎方程式(式34から式36)と本明細書で例示した4組の基礎方程式(式39および式40の組、式41から式43の組、式57並びに式67)との関係について簡単に述べておく。
【0268】
前述したように一般的には、最も基本的な基礎方程式に従って平面を構成すれば、関数RSi(qμ)が求められる。しかしながら、上記平面を構成する際に実は2つの問題がある。1つは、式34から式36を同時に厳密に満たす平面が一般のポテンシャルエネルギーVに対して存在するとは限らないという原理的な問題である。そして、もう1つは、式34から式36に従い平面を構成しようとすると逆走が生じる場合があるという実用上の問題(前述した逆走の問題)である。特に、後者の問題の原因は明白で、式35と式36がφiμ(RS)とρμ(RS)の符号を決定することができない、つまり符号反転に対して不変であるからである。
【0269】
そこで、この原理的な問題に対しては、式34から式36を満たす平面を近似的に見出す方法が有効である。実際に非特許文献10では、式36をあきらめ、式34と式35のみで平面を構成しようとした。そしてこの策は、K=1の場合には式39および式40を与える。また、本発明者は、式34を最大限満たしつつ、式35と式36で平面を構成する方法を開発した。そしてこの策は、K=1の場合には式41から式43を与える。具体的には、式35と式36においてφi1(RS)とρ1(RS)を消去することで得られる下記式80を、RSの他、λ1もq1の関数とみなしてq1で微分すると式41から式43が得られる。
【0270】
【数49】
【0271】
一方、実用上の問題に対しては、式34から式36においてφiμ、λμおよびρμをRSから独立しかつq†の関数としての変数(すなわち補助座標)とみなす方法が有効である。これは、φiμ、λμおよびρμをRSから独立した変数とみなすと、φiμ、λμおよびρμは、RSと同様に、平面を構成する際に少しずつ連続的に変化する量となり、突然符号が変わるような事態が発生しないからである。φiμ、λμおよびρμを補助座標として取り扱うためには、φiμ、λμおよびρμについても式34のようなqνでの微分を含む式が必要である。そこで、式35および式36をそれぞれqνで全微分し、これをdqνで割ることにより、下記式81および式82が得られる。これにより、最も基本的な基礎方程式は、式34、式81および式82からなる実用上の問題のない基礎方程式へと変換されたことになる。
【0272】
【数50】
【0273】
そして、式34、式81および式82からなる基礎方程式における原理的な問題に対処することを考える。その対処方法は、例えば上記と同様に、式34、式81および式82を満たす平面を近似的に見出す方法である。具体的には、式34および式81のみで平面を構成する策は式57を与える。一方、式34を最大限満たしつつ、式81と式82で平面を構成する策は式67を与える。
【産業上の利用可能性】
【0274】
各実施形態では、本発明をタンパク質および薬物からなる複合体の結合過程に適用した場合について説明したが、本発明はこれに限られるものではない。つまり、本発明は、より一般的に生体高分子および結合候補分子からなる複合体の形成過程について適用することが可能である。また、本発明は、複合体の形成過程の他、複合体の解離過程、タンパク質のapo体の構造および生体高分子のフォールディング(折り畳み機構)等、生体高分子を含む多原子系の挙動の解析についても適用することが可能である。
【符号の説明】
【0275】
2 タンパク質
4 ポケット
6 結合候補分子
10 シミュレーション装置
12 入力手段
14 制御手段
16 座標設定手段
18 座標抽出手段
20 逆変換手段
22 表示手段
A 複合体が解離体である状態を示す点
B エネルギー障壁
C 結合が完了した複合体である状態を示す点
RP 複合体の挙動の反応経路
【技術分野】
【0001】
本発明は、コンピュータを用いて、モデル化された質点系の動力学的な挙動を予測するシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体に関するものである。
【背景技術】
【0002】
コンピュータ技術の発展に伴って、シミュレーションによる理論的な計算により、タンパク質、核酸、脂質および多糖等の生体高分子の動力学的な大変形の挙動を原子レベルで理論的に解明しようという研究が盛んである。具体的には、標的とするタンパク質と結合候補分子(当該タンパク質に対する結合の適合性を有するか否かについての解析の対象となる分子)との互いの親和性を理論的に予測する創薬スクリーニング、および、タンパク質の一次配列から立体構造を構築する原理を明らかにして一次構造から高次構造を理論的に構築するタンパク質のフォールディング(折り畳み機構)の解析が挙げられる。
【0003】
例えば、生体高分子の動力学的な挙動を予測するシミュレーション方法として、巨大分子でも原子レベルでシミュレーションを実行することが可能な分子動力学法またはモンテカルロ法等に基づくシミュレーション方法(非特許文献1から7)がある。分子動力学法によれば、運動方程式に従って多原子系の時間発展を微小時間単位で順次追跡することができる。しかし、生体高分子を含む多原子系のポテンシャル曲面には、多くのローカルミニマム(極小点)或いは多くのエネルギー障壁が存在しているため、上記のような方法では、多原子系の状態が初期構造に近いローカルミニマムに捕捉されて、計算に膨大な時間を要するという問題があった。また、この問題はモンテカルロ法でも同じである。
【0004】
そこで、上記のようなローカルミニマムによるトラップの問題を解消する方法として、分子動力学法に基づく計算において「強制的な力」を印加することによりローカルミニマムによるトラップを脱出させる計算方法が知られている。例えば、特許文献1から3並びに非特許文献8および9では、多原子系が一度経験した構造を二度と取らないように仮想的な相互作用を逐次導入することや、高温相の多原子系の加速された運動で見つかった構造を低温相の多原子系の運動に逐次反映させることが開示されている。このような相互作用の導入等により、多原子系が通過するポテンシャル曲面上の経路の探索を加速して、生体高分子の動力学的な挙動を計算することが可能となる。
【先行技術文献】
【特許文献】
【0005】
【特許文献1】特開2005−267592号公報
【特許文献2】特開2005−524129号公報
【特許文献3】再表2006/068271号公報
【非特許文献】
【0006】
【非特許文献1】Todd J. A. Ewing, Irwin D. Kuntz, Journal of Computational Chemistry, Vol. 18, Issue 9 (1997), p. 1175-1189
【非特許文献2】Garrett M. Morris, David S. Goodsell, et al., Journal of Computational Chemistry, Vol. 19, Issue 14 (1998), p. 1639-1662
【非特許文献3】Matthias Rarey, Bernd Kramer, et al., Journal of Molecular Biology, Vol. 261, Issue 3 (1996), p. 470-489
【非特許文献4】Ruben Abagyan, Maxim Totrov, et al., Journal of Computational Chemistry, Vol. 15, Issue 5 (1994), p. 488-506
【非特許文献5】Gareth Jones, Peter Willett, et al., Journal of Molecular Biology, Vol. 267, Issue 3 (1997), p. 727-748
【非特許文献6】Richard A. Friesner, Jay L. Banks, et al., Journal of Medicinal Chemistry, Vol. 47, Issue 7 (2004), p. 1739-1749
【非特許文献7】Thomas A. Halgren, Robert B. Murphy, et al., Journal of Medicinal Chemistry, Vol. 47, Issue 7 (2004), p. 1750-1759
【非特許文献8】Yoshifumi Fukunishi, Yoshiaki Mikami, et al., Journal of Physical Chemistry B, Vol. 107, Issue 47 (2003), p. 13201-13210
【非特許文献9】Yuji Sugita, Yuko Okamoto, Chemical Physics Letters, Vol. 314, Issue 1-2 (1999), p. 141-151
【発明の概要】
【発明が解決しようとする課題】
【0007】
しかしながら、生体高分子等がモデル化された質点系の運動を記述する運動方程式に基づいて、時間に関する積分で当該質点系の挙動を予測する分子動力学を用いた計算では、シミュレーションの計算時間および計算精度がトレードオフの関係になってしまうという問題がある。具体的には、特許文献1から3並びに非特許文献8および9のようなシミュレーション方法では、計算時間は短縮されるものの、印加した「強制的な力」によって質点系の本来有する運動秩序が破壊され、計算結果として非現実的な大変形の挙動を算出してしまう。そして、より効率よくローカルミニマムによるトラップを脱出させるために「強制的な力」の作用を大きくすればするほど、その傾向は顕著となる。
【0008】
本発明は上記問題に鑑みてなされたものであり、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することを可能とするシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体を提供することを目的とするものである。
【課題を解決するための手段】
【0009】
上記課題を解決するために、本発明に係るシミュレーション装置は、
モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション装置において、
質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定する座標設定手段と、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるK個の集団座標の関数としての中継座標の構成を求める座標抽出手段と、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測する逆変換手段とを備えることを特徴とするものである。
【0010】
なお本明細書において、K、MおよびNは、K<M<3Nの関係を満たし、それぞれ1以上の整数を表す。
【0011】
本明細書において、「質点系の構造」とは、当該質点系を構成するN個の質点が造り出す立体構造を意味する。
【0012】
「質点系の構造の変化を主として担う中継座標」とは、当該座標の変化が当該質点系の構造の立体構成に大きな影響を与えるような座標を意味する。
【0013】
「中継座標を設定し」とは、中継座標として、一部の質点座標そのもの、質点座標を組み合わせて定義できる座標又はこれらの組み合わせを設定することを意味する。非中継座標についても同様である。
【0014】
そして、本発明に係るシミュレーション装置において、座標抽出手段は、
中継座標および非中継座標の関数として表現された質点系のポテンシャルエネルギーを求める第1の工程を実施し、
ポテンシャルエネルギーを用いた断熱近似条件に従って非中継座標を中継座標に従属させる第2の工程を実施し、
現状の中継座標および非中継座標の下で、上記影響を考慮しながら、中継座標に関するポテンシャルエネルギーの微分係数を求める第3の工程を実施し、
ポテンシャルエネルギーの上記微分係数の下で、自己無撞着集団座標法におけるポテンシャルエネルギーの微分係数を用いた基礎方程式に従って、集団座標に関する中継座標の微分係数を求める第4の工程を実施し、
中継座標の上記微分係数の下で集団座標を微小量だけ更新することにより、更新された中継座標を求める第5の工程を実施し、
更新された中継座標の下で、中継座標に従属した非中継座標について構造緩和を行う第6の工程を実施し、
その後、非中継座標の構造緩和後の状態における中継座標および非中継座標に基づいて、第3の工程から第6の工程までを順次繰り返すことにより、中継座標の上記構成を求めるものであることが好ましい。
【0015】
「現状の中継座標および非中継座標の下で」とは、1回目の第3の工程においては、座標設定手段により設定された状態における中継座標および非中継座標の下で、2回目以降の第3の工程においては、直前の第6の工程で行われた非中継座標の構造緩和後の状態における中継座標および非中継座標の下で、その第3の工程を実施することを意味する。
【0016】
また、本発明に係るシミュレーション装置において、上記影響を考慮する方法は下記式1から式3の少なくとも1つを使用することであることが好ましい。
【0017】
【数1】
【0018】
本明細書において、i、jおよびkは、それぞれ1からMまでの整数を表し、
α、βおよびγは、それぞれ1から3N−Mまでの整数を表し、
RSiは、質点系におけるi番目の中継座標を表し、
RFαは、質点系におけるα番目の非中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
RFは(RF1,RF2,・・・,RF(3N−M))を表し、
RF(RS)は、中継座標に従属させられた非中継座標を表し、
V(RS,RF)は、中継座標と非中継座標で表された質点系のポテンシャルエネルギーを表し、
Veff(RS)は、V(RS,RF)にRF(RS)を代入して得られる有効ポテンシャルエネルギーを表す。
【0019】
式2において、第3項の(i←→j)は、第2項におけるiおよびjの文字を相互に入れ替えた項(つまり、第2項におけるiをjへ、jをiへ書き換えた項。以下同様である。)を表す。
【0020】
式3において、第3項の(i←→k)は、第2項におけるiおよびkの文字を相互に入れ替えた項を表し、第4項の(j←→k)は、第2項におけるjおよびkの文字を相互に入れ替えた項を表し、第6項の(i←→k)は、第5項におけるiおよびkの文字を相互に入れ替えた項を表し、第7項の(j←→k)は、第5項におけるjおよびkの文字を相互に入れ替えた項を表す。
【0021】
さらに、式1から式3において下記式4を用いた。
【0022】
【数2】
【0023】
Kαβ−1(Rs)はKαβ(Rs)の逆行列を表し、
Kαβ(Rs)およびJαi(Rs)はそれぞれ式5および式6で定義される。
【0024】
【数3】
【0025】
また、本発明に係るシミュレーション装置において、断熱近似条件は下記式7であることが好ましい。
【0026】
【数4】
【0027】
また、本発明に係るシミュレーション装置において、集団座標の個数KはK=1を満たすものであり、自己無撞着集団座標法における基礎方程式は下記式8および9であることが好ましい。
【0028】
【数5】
【0029】
本明細書において、q1は、集団座標を表し、
miは質点系におけるi番目の中継座標に関する質量を表し、
φi(RS)は、式9を満たす関数(固有ベクトル)のi番目の成分を表し、
Λ(RS)は、式9を満たす関数(固有値)を表す。
【0030】
或いは、本発明に係るシミュレーション装置において、集団座標の個数KはK=1を満たすものであり、自己無撞着集団座標法における基礎方程式は下記式10から12であることが好ましい。
【0031】
【数6】
【0032】
φi(RS,λ)およびκ(RS,λ)は式12を満たす関数であり、φi(RS,λ)はそのi番目の成分を表し、
λは、補助座標(中継座標から独立しかつ集団座標の関数として取り扱う変数)を表す。
【0033】
或いは、本発明に係るシミュレーション装置において、座標抽出手段は、第4の工程において、基礎方程式中の集団座標に関する中継座標または補助座標の微分係数が有する符号の任意性を解消するように、この基礎方程式を解く上で中継座標から独立しかつ集団座標の関数として取り扱う変数(つまり補助座標)の数を増やして、計算を行うものであることが好ましい。
【0034】
また、補助座標の数を増やす場合において、座標抽出手段は、その結果得られた下記式13で表される基礎方程式に従って計算を行うものであることが好ましい。さらにこの場合において集団座標の個数KはK=1を満たすことが好ましい。
【0035】
【数7】
【0036】
Yは、下記式14で定義されるMK+M+K次元のベクトルであり、vμは、下記式15の非斉次線形方程式の解ベクトルである。
【0037】
【数8】
【0038】
式15におけるCおよびsμはそれぞれ式16および式17で定義される。
【0039】
【数9】
【0040】
V,ij(RS)およびV,ijk(RS)はそれぞれ式18および式19で定義される。
【0041】
【数10】
【0042】
なお、μおよびνはそれぞれ1からKまでの整数を表し、
qμはμ番目の集団座標を表し、
φiμおよびΛμは補助座標を表す。
【0043】
或いは、補助座標の数を増やす場合において、座標抽出手段は、その結果得られた下記式20で表される基礎方程式に従って計算を行うものであることが好ましい。さらにこの場合において集団座標の個数KはK=1を満たすことが好ましい。
【0044】
【数11】
【0045】
Zは、下記式21で定義されるMK+M+2K次元のベクトルであり、cμνは、μごとに定義される下記式22で表される値をそれぞれ最小にするように一意に決められる定数であり、wμは、下記式23で定義された行列DのK次元の特異値空間を張るK個のMK+M+2K次元の単位ベクトルである。なお、λμおよびρμは、φiμと同様に補助座標を表す。
【0046】
【数12】
【0047】
また、本発明に係るシミュレーション装置において、座標抽出手段は、ポテンシャルエネルギーの3次微分の項の計算を下記式24に基づいて行うものであることが好ましい。φiは任意のベクトルのi番目の成分を表し、nは式25で定義される。
【0048】
【数13】
【0049】
また、本発明に係るシミュレーション装置において、座標設定手段は、質点系の構造のうち特徴的な部分構造ごとに抽出された当該部分構造を代表する代表座標を中継座標として設定するものであることが好ましい。
【0050】
本明細書において、「特徴的な部分構造」とは、質点系の構造のうち、形態的および/または機能的な特徴を有する部分構造を意味する。
【0051】
また、本発明に係るシミュレーション装置において、質点系は生体高分子を含む多原子系であり、
上記部分構造は、生体高分子の2次構造、生体高分子の構成単位または生体高分子の主鎖であり、
上記部分構造についての代表座標は、その部分構造を構成する原子の座標そのもの、その原子の座標を組み合わせて定義できる座標またはその部分構造の周期間隔とすることができる。
【0052】
また、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子はタンパク質であり、上記部分構造はタンパク質の2次構造であり、この2次構造についての代表座標は、その2次構造を構成する原子群の重心座標または2次構造の屈曲角であることが好ましい。この場合において、2次構造はヘリックス構造、βシート、ターン、ループおよびランダムコイルの少なくとも1つであることが好ましい。
【0053】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子がタンパク質であり、上記部分構造がタンパク質の残基であり、この残基についての代表座標が、その残基を構成する原子群の重心座標であることが好ましい。
【0054】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子がタンパク質であり、上記部分構造がタンパク質の主鎖であり、この主鎖についての代表座標が、その主鎖を構成する原子のそれぞれの座標であることが好ましい。
【0055】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の2次構造であり、この2次構造についての代表座標が、その2次構造を構成する原子群の重心座標または2次構造の屈曲角であることが好ましい。この場合において、2次構造は螺旋構造であることが好ましい。
【0056】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の残基であり、この残基についての代表座標が、その残基を構成する原子群の重心座標であることが好ましい。
【0057】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の主鎖であり、この主鎖についての代表座標が、その主鎖を構成する原子のそれぞれの座標であることが好ましい。
【0058】
或いは、本発明に係るシミュレーション装置において、質点系が生体高分子を含む多原子系である場合には、生体高分子が核酸であり、上記部分構造が核酸の螺旋構造であり、この螺旋構造についての代表座標が、その螺旋構造の周期間隔であることが好ましい。
【0059】
また、本発明に係るシミュレーション装置において、多原子系はさらに、生体高分子に対する結合候補分子を含むものとすることができる。
【0060】
本発明に係るシミュレーション方法は、
上記に記載のシミュレーション装置によって実行される、モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション方法であって、
質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、
質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定し、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、
中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるK個の集団座標の関数としての中継座標の構成を求め、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測することを特徴とするものである。
【0061】
本発明に係るシミュレーションプログラムは、上記に記載のシミュレーション方法をコンピュータに実行させることを特徴とするものである。
【0062】
本発明に係るコンピュータ読み取り可能な記録媒体は、上記に記載のシミュレーションプログラムを記録したものである。
【発明の効果】
【0063】
本発明に係るシミュレーション装置は、上記のような座標設定手段と、座標抽出手段と、逆変換手段とを備え、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0064】
本発明に係るシミュレーション方法は、上記シミュレーション装置によって実行される方法であって、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0065】
本発明に係るプログラムおよび記録媒体は、上記シミュレーション方法を実行させることが可能であるから、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【図面の簡単な説明】
【0066】
【図1】本発明のシミュレーション方法により誘導適合を考慮した場合における、タンパク質と結合候補分子の結合過程を示す概略図である。
【図2】タンパク質および結合候補分子を含む多原子系の誘導適合を考慮した概念的なポテンシャル曲面における反応経路を示す概略図である。
【図3】集団座標の抽出を利用して多原子系の反応経路を求める手順を概念的に示す図である。
【図4】本発明における変数変換および逆変換の概念を示す図である。
【図5】実施形態のシミュレーション装置の構成を概略的に示すブロック図である。
【図6A】実施形態のシミュレーション方法における計算手順を概略的に示すフローチャート図である。
【図6B】実施形態のシミュレーション方法における計算手順を概略的に示すフローチャート図である。
【図7】実施形態における中継座標RS、非中継座標RFおよび原子の座標xの関係を示す概略図である。
【図8】実施形態のシミュレーション方法によって得られた所定の複合体の反応経路を示すグラフである。
【図9】図8のグラフにおける複合体の結合過程を示す概略図である。図9aは複合体の結合過程における開始状態を示し、図9bは複合体の結合過程における最終状態を示す。
【図10】符号の任意性が解消された計算の場合と符号の任意性が残ったままの計算の場合との比較結果を示すグラフである。
【発明を実施するための形態】
【0067】
以下、本発明の実施形態について図面を用いて説明するが、本発明はこれに限られるものではない。なお、視認しやすくするため、図面中の各構成要素の縮尺等は実際のものとは適宜異ならせてある。
【0068】
まず、詳細な実施形態の説明を行う前に、従来技術に対して本発明が有する技術上の意義を明らかにするべく、本発明の基本となる技術的思想およびそれに至った経緯を概説する。ここでは、説明を分かりやすくするために、具体的に生体高分子および結合候補分子からなる複合体を含む多原子系の動力学的な挙動を予測する場合を例にして説明する。
【0069】
上記のような多原子系の動力学的な挙動の予測は、例えば、開発の期間短縮およびコスト削減等の観点から新薬開発において重要な役割を果たす。薬物等の生理活性物質は特定のタンパク質に結合することによりその化学的性質を発揮する。したがって、新薬開発においては、標的とするタンパク質に結合しやすい薬物候補を数十万から数百万とも言われる膨大な化合物から絞り込む必要がある(創薬スクリーニング)。また究極的には、薬物候補を膨大な化合物から数個程度にまで絞り込むことが要求されている。そこで、薬物候補を膨大な化合物から効果的に絞り込むためには多原子系の挙動を原子レベルで取り扱い、分子間さらには原子間の相互作用を正確に評価する必要がある。
【0070】
図1は、タンパク質の動力学的な挙動を考慮した場合における、タンパク質および結合候補分子の結合過程を示す概略図である。タンパク質2は、結合候補分子6と結合するためのポケット4を有している。つまり、当該結合候補分子6は、当該タンパク質2の生理活性物質と言える。しかしながら、当該ポケット4は、必ずしもそのままで結合候補分子6がその中に収まることができるような形状を有しているわけではない(図1a)。例えば図1では、ポケット4の開口が結合候補分子6の大きさよりも小さいことが示されている。そこで、このような場合には、タンパク質2は、接近した結合候補分子6との相互作用に応じて、開口の形状および大きさを結合候補分子6の形状および大きさに適合させながらその構造を変化させる(図1bおよびc)。
【0071】
このように、タンパク質および核酸等の生体高分子が生理活性物質との相互作用によりその構造を変化させる現象は「誘導適合」と言われる。分子間の相互作用を正確に計算するためには、当然この誘導適合を考慮しなければならない。そして、これは、多原子系を原子レベルで取り扱うことによって実現できる。
【0072】
一方、図2は、多原子系の誘導適合を考慮した概念的なポテンシャル曲面上の反応経路を示す概略図である。図2の横軸X1はタンパク質の構造の変化を概念的に示す軸であり、縦軸X2はタンパク質および結合候補分子の距離を示す軸である。つまり、図2は、タンパク質の構造の変化および分子間の距離に応じた多原子系のポテンシャル曲面を表している。図2に示されるように、誘導適合が生じる系は、ポテンシャル曲面上の2つの安定点(つまり、複合体が解離体である状態を示す点Aおよび複合体の結合が完了した状態を示す点C)を結ぶ反応経路RP上にエネルギー障壁B(鞍点)を有することが解かる。
【0073】
したがって、誘導適合を考慮しながらタンパク質2と結合候補分子6との結合過程を解析するということは、タンパク質2および結合候補分子6に関する上記ポテンシャル曲面上において、エネルギー障壁Bを超えて2つの安定点AおよびCを結ぶ反応経路を探索することに帰着する。
【0074】
以上を考慮すると、多原子系の動力学的な挙動を予測するシミュレーションにおいて計算精度を向上させかつ計算時間を短縮するためには、多原子系の原子レベルでの取り扱いを可能にすること、および非経験的な反応経路の探索を容易にすることが必要となる。
【0075】
しかしながら、数千から数百万個の原子で構成される多原子系の挙動は、1秒〜1時間オーダーの時間規模で起こる緩慢な大変形であるため、従来のシミュレーション方法では、原子レベルでの取り扱いは可能であるものの、理論的な計算を行うために数十年もの膨大な時間が必要であるという問題があった。そこで、多原子系の挙動の理論的な計算を実用的な期間内で実施できるようにするために、様々なシミュレーション手法が検討・開発されている。前述した「強制的な力」を印加する分子動力学法もその1つであるが、その問題点は既に述べた通りである。また、標的とするタンパク質を剛体として近似して、取り扱う変数を減らすことで計算量を減らす手法も検討されている。しかし、このような粗い近似の下では、タンパク質の動力学的な挙動の影響を当然考慮できず、結果としてこれらの分子間に働く相互作用を正しく計算できない。このような場合、例えば100万個ある候補の化合物から薬物候補として1万個程度に絞り込むことしかできない。また、その他のシミュレーション手法でも計算精度および計算時間のトレードオフの問題は解決するに至っていない。
【0076】
そこで、本発明者は、生体高分子を含む多原子系をN個の質点から構成される質点系にモデル化した上で、この質点系の挙動を集団運動として取り扱い、質点座標(質点系の構造を一次的に記述する質点の座標)の自由度より小さい自由度の集団座標を、その質点座標に基づいた集団運動の中から抽出することに着目した。
【0077】
(集団座標)
集団座標について説明する。一般的に「集団座標」とは、集団変数の座標要素を意味する。そして、「集団変数」とは、質点系の正準変数(その質点系における質点座標および質点の運動量)と正準変換により対応づけられる一般変数(q,p)のうち、その正準変数よりも自由度が小さいものを意味する。なお、正準変換により対応づけられる一般変数(q,p)において、qは座標要素であり、pは運動量要素である。qおよびpは下記式26および式27によってそれぞれ定義される。
【0078】
【数14】
【0079】
式26および式27において、ηは1から3Nまでの整数を表す。
【0080】
したがって、集団変数は、下記式28および式29によってそれぞれ定義されるq†およびp†を用いて(q†,p†)と表すことができる。
【0081】
【数15】
【0082】
式28および式29は、集団変数(q†,p†)が式28および式29に示されたμ=1〜Kまでの2K個の変数により構成されることを表している。言い換えれば、集団変数(q†,p†)とは、時間に依存して変化する可変成分(q1,q2,・・・,qK,p1,p2,・・・,pK)および時間変化に対して定数となる不変成分(qK+1=0,qK+2=0,・・・,q3N=0,pK+1=0,pK+2=0,・・・,p3N=0)から構成される一般座標の当該可変成分であると言える。なお、式28および式29では、不変成分が定数としてゼロの値を取る場合を示しているが、定数はゼロに限られない。
【0083】
そして、集団座標は、集団変数(q†,p†)の座標要素であるから、q†と表すことができる。より具体的には、集団座標とは、(q1,q2,・・・,qK)から構成される変数の集合である。なお、集団座標自体は、質点座標のみから正準点変換により求めることもできる。
【0084】
上記のように、集団運動として取り扱う対象となる座標(例えば、質点座標および後述する中継座標)からより自由度の低い集団座標へと変数変換を行うことを「変数分離」ともいう。
【0085】
集団座標の自由度は、変数分離の対象となる座標の自由度よりも小さければ、特に限定されない。しかしながら、集団座標の自由度が小さい程、その後の計算処理が容易となるため、集団座標の自由度は1(つまりK=1)であることが好ましい。
【0086】
(集団運動としての取り扱い)
質点座標に対して変数分離を実施すると、集団座標q†によって、質点系の挙動における本質的な運動、つまり、より低次元的な集団運動の記述が可能となる。この結果、3N次元の質点座標による運動の記述がK次元の集団座標q†による運動の記述へと置き換えられ、質点系の運動の記述が単純化される。具体的には、質点座標で記述された運動方程式は、集団座標q†を引数として持つ質点座標xの構成x=x(q†)とq†についての運動方程式との組み合わせに単純化される。なお、xは(x1,x2,・・・,x3N)を表す。質点座標x(q†)の具体的な構成は質点系の集団運動としての挙動の仕方を表し、集団座標q†についての運動方程式はその挙動の基本法則を表す。したがって、運動方程式を解く上で、取り扱う変数の自由度がより減少することになるため、理論的な計算が容易になる。
【0087】
そして、質点座標xに対して変数分離を実施して集団座標q†に関するポテンシャル曲面上の反応経路を見つけ出し、座標を集団座標q†から質点座標xに逆変換して、質点座標xに関するポテンシャル曲面上で反応経路を表現することにより、質点系の反応経路を求めることができる。上記のような方法によれば、低次元の運動方程式を解いた後、集団座標q†から質点座標xへの変数変換を実行するだけであるため、本質的に質点系の状態がローカルミニマムに捕捉されることが起こり得ない。
【0088】
例えば、図2に示された複合体を含む多原子系の反応経路を求める場合の概念的な手順は以下のようになる。図3は、集団座標の抽出を利用して多原子系の反応経路を求める手順を概念的に示す図である。まず、質点座標X1およびX2に関するポテンシャル曲面(図3a)が、抽出された集団座標q1およびq2に関するポテンシャル曲面(図3b)に変換される。この際、集団座標として、2つの安定点AおよびC並びにエネルギー障壁Bが一直線上に並ぶような集団座標を採用することが重要である。このような集団座標q1およびq2に関するポテンシャル曲面上では、2つの安定点AおよびCを結ぶ反応経路RPを一気に構成することができる(図3b)。その後、座標の逆変換により、変数X1およびX2に関するポテンシャル曲面上での反応経路RPが求められる(図3c)。ここで得られた一次元的な軌跡こそ理論化学の分野で古くから用いられてきた「反応座標」と同定できる変数である。
【0089】
(生体高分子を含む多原子系の挙動を集団運動として取り扱う際の課題)
変数分離は集団運動理論に基づいて行われる。例えば、集団運動理論の1つとして、自己無撞着集団座標法(SCC法:Self-consistent Collective Coordinate method)がある。SCC法は、現在原子核の分野で研究が進められている物理学における集団運動理論の1つである。より具体的には、SCC法は、系の集団運動を記述する本質的な集団変数を求め、その系のハミルトニアンに含まれる正準変数が張る6N次元の位相空間から、動力学的な集団運動に関与する分離多様体を見出し、その分離多様体上またはその近傍で定義されたハミルトニアンによって集団運動を記述する運動理論である。「分離多様体」とは、その系のハミルトニアンに含まれる正準変数が張る6N次元の位相空間中の部分空間であって、その系の挙動を表す軌道が閉じ込められる部分空間を意味する。つまり、分離多様体中の任意の1点を初期値とした場合、上記軌道は常にこの分離多様体中に拘束されることになる。SCC法に関する詳細については、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らの文献を非特許文献10という。)等を参照されたい。
【0090】
そこで、例えば原子核反応において取り扱う核子(陽子および中性子)を化学反応において取り扱う原子に置き換えてSCC法を適用することにより、多原子系の動力学的な挙動を集団運動として取り扱ってこの多原子系の反応経路を求めることができるようにも思える。
【0091】
しかし、本発明者は、本発明の対象である生体高分子を含むような大規模な系に対してSCC法等の既存の集団運動理論をそのまま適用するだけでは、計算精度を向上させかつ計算時間を短縮するには不十分であるという問題を新たに見出した。これは、原子核反応の場合にはせいぜい数百個の核子を取り扱えばよいのに対して、水中における生体高分子を含む多原子系の化学反応の場合には数千から数百万個の原子を取り扱わなければならないことに起因する。つまり、集団運動理論とはいえ既存の手法では変数の制御、つまり変数分離の計算が複雑になるためである。
【0092】
そして、上記のような問題は、複合体を含む多原子系の反応経路を求める場合に限られるものではない。つまり、より一般的に、生体高分子を含む系の反応経路を求める場合に同様の問題が生ずることは容易に想像される。
【0093】
「発明の説明」
以上のような経緯を踏まえ、本発明者は、質点座標の自由度より小さい自由度の集団座標によって集団運動を記述する集団運動理論であって、巨大な原子群である生体高分子を含む多原子系の構造計算にも適用可能な新規な集団運動理論に基づいて、質点系の動力学的な挙動を予測するシミュレーション装置およびシミュレーション方法並びにその方法を実行するためのプログラムおよび記録媒体を発明するに至った。
【0094】
より具体的には、本発明に係るシミュレーション装置は、
モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション装置において、
質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定する座標設定手段と、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるK個の集団座標の関数としての中継座標の構成を求める座標抽出手段と、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測する逆変換手段とを備えることを特徴とするものである。
【0095】
そして、本発明に係るシミュレーション方法は、
上記に記載のシミュレーション装置によって実行される、モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション方法であって、
質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、
質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定し、
非中継座標を中継座標に従属させることにより中継座標の関数としての非中継座標の構成を求め、
中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される一般座標の上記可変成分であるK個の集団座標の関数としての中継座標の構成を求め、
集団座標についての運動方程式の解として得られる時間の関数としての集団座標、中継座標の上記構成および非中継座標の上記構成に基づいて、質点座標の時間発展を予測することを特徴とするものである。
【0096】
また、本発明に係るシミュレーションプログラムは、上記に記載のシミュレーション方法をコンピュータに実行させることを特徴とするものである。
【0097】
また、本発明に係るコンピュータ読み取り可能な記録媒体は、上記に記載のシミュレーションプログラムを記録したものである。
【0098】
なお、「N個の質点から構成された質点系」とは、質点系の構成要素である質点の総数がN個であることを意味する。
【0099】
<座標設定手段>
座標設定手段は、質点系の構造を記述する3N個の質点座標に基づいて、質点系の構造の変化を主として担うM個の座標である中継座標を設定し、質点系の構造を記述しかつ中継座標から独立した座標である非中継座標を設定するものである。上記「3N個の質点座標」とは、より具体的には3次元空間におけるN個の質点の座標xである。
【0100】
まず、座標設定手段は、質点座標の構成に基づいて、質点系の巨視的な集団運動(大振幅集団運動)を記述する上で質点系の構造の変化を主として担うM個の中継座標を設定する。そして、後述する座標抽出手段はこの設定された中継座標に対して変数分離を実施することになる。言い換えれば、本発明では、質点系を記述する独立した座標のうち大振幅集団運動による質点系の構造変化に対して影響の少ない座標は既知であるので、このような座標は集団運動理論における変数分離の対象から予め除外される。つまり、3N次元の質点座標から集団座標を抽出するのではなく、M次元の中継座標から集団座標を抽出することにより、変数分離の計算が単純化される。Mは、K<M<3Nを満たす整数である。また、Kは、上記式28および式29での説明同様、集団座標q†の要素の個数である。
【0101】
なお、この中継座標からさらに変数の階層化を行ってもよい。つまり、上記のように設定された中継座標の構成に基づいて、さらに低次元の2次的な中継座標を設定し、この2次的な中継座標から集団座標を抽出することも可能である。
【0102】
中継座標は、大振幅集団運動による質点系の構造変化に対する影響の度合を基準に設定される。そして、この影響の度合が大きい座標ほど中継座標としてふさわしい。言い換えれば、どのような質点系の大振幅集団運動を解析の対象としているのかによって、中継座標の設定方法および具体的な座標の内容が変わる。
【0103】
具体的には、中継座標は、質点系に含まれる一部の質点の座標そのもの、質点の座標を組み合わせて定義できる座標又はこれらの組み合わせを用いて設定される。また中継座標は、例えば質点系の構造のうち特徴的な部分構造に関連して設定されることが好ましい。そしてこの場合において、中継座標は、特徴的な部分構造ごとに抽出された、当該部分構造を代表する代表座標であることがより好ましい。代表座標は、1つの特徴的な部分構造に複数あってもよい。
【0104】
「特徴的な部分構造」とは、質点系の構造のうち、形態的および/または機能的な特徴を有する部分構造を意味する。特徴的な部分構造は、例えば質点系が生体高分子を含む多原子系である場合には、生体高分子の2次構造(部分的な折り畳み構造)、生体高分子の構成単位および生体高分子の主鎖である。また、特徴的な部分構造の「代表座標」とは、その特徴的な部分構造を代表して表す座標を意味する。代表座標は、例えば、その特徴的な部分構造を構成する質点(原子)の座標そのもの、その質点(原子)の座標を組み合わせて定義できる座標またはその部分構造の周期間隔等である。
【0105】
より具体的には、上記のような特徴的な部分構造としては、生体高分子がタンパク質である場合には、タンパク質の2次構造が挙げられる。具体的な2次構造としては、ヘリックス構造(310ヘリックス、αヘリックス、πヘリックスおよびβヘリックス等)、βシート、ターン、ループおよびランダムコイル等が挙げられる。そして、タンパク質の2次構造についての代表座標としては、その2次構造を構成する原子群の重心座標または2次構造の屈曲角が挙げられる。例えば、一般的なタンパク質中に含まれる上記のような特徴的な部分構造(高次構造)の個数は、せいぜい数個から数十個である。これにより、変数分離の対象となる変数の個数を大幅に低減することができる。また、生体高分子と結合候補分子との結合反応における挙動を予想する場合には、結合候補分子の構造自体を多原子系の特徴的な部分構造の1つとして取り扱うことも可能である。
【0106】
或いは、特徴的な部分構造としては、生体高分子がタンパク質である場合には、タンパク質の残基(タンパク質の構成単位であるアミノ酸部分。N末端残基およびC末端残基を含む)とすることもできる。そして、タンパク質の残基についての代表座標としては、その残基を構成する原子群の重心座標が挙げられる。
【0107】
或いは、特徴的な部分構造としては、生体高分子がタンパク質である場合には、タンパク質の主鎖とすることもできる。そして、タンパク質の主鎖についての代表座標としては、主鎖を構成する原子のそれぞれの座標が挙げられる。
【0108】
或いは、特徴的な部分構造としては、生体高分子が核酸である場合には、核酸の2次構造とすることもできる。具体的な2次構造としては、螺旋構造等が挙げられる。螺旋構造の場合には、所定の周期ごとに特徴的な部分構造として取り扱う。そして、核酸の2次構造についての代表座標としては、その2次構造を構成する原子群の重心座標または2次構造の屈曲角が挙げられる。
【0109】
或いは、特徴的な部分構造としては、生体高分子が核酸である場合には、核酸の残基(核酸の構成単位であるヌクレオチド部分)とすることもできる。そして、核酸の残基についての代表座標としては、その残基を構成する原子群の重心座標が挙げられる。
【0110】
或いは、特徴的な部分構造としては、生体高分子が核酸である場合には、核酸の主鎖とすることもできる。そして、核酸の主鎖についての代表座標としては、主鎖を構成する原子のそれぞれの座標が挙げられる。
【0111】
或いは、生体高分子が核酸でありかつ核酸の螺旋構造である場合には、螺旋構造の周期間隔を代表座標とすることもできる。
【0112】
特徴的な部分構造を実際にどのように設定するかについては、シミュレーションを行う対象の質点系の構成および/または現象に応じて適宜決定される。例えば、タンパク質の折りたたみ構造の形成過程や核酸の螺旋構造の形成過程を解析の対象とする場合には、生体高分子の構成単位や生体高分子の基本的な骨格とも言える主鎖を特徴的な部分構造として採用することが好ましい。生体高分子の主鎖については、共有結合により連なった一本の主鎖全体を1つの特徴的な部分構造としてもよいし、上記一本の主鎖を分割して得られた部分的な構造のそれぞれを特徴的な部分構造としてもよい。また、タンパク質と結合候補分子との結合反応における挙動を解析の対象にする場合には、上記で挙げた主鎖の他、へリックス構造、βシート、ターン、ループおよびランダムコイル等のそれぞれ、又はこれらの組み合わせを特徴的な部分構造として採用することもできる。特徴的な部分構造は複数の2次構造を含んでもよい。つまり、2次構造のそれぞれを1つの特徴的な部分構造とすることもでき、主鎖に沿って隣接した複数の2次構造をまとめて1つの特徴的な部分構造とすることもできる。また、核酸の螺旋構造に阻害剤が挟まって結合する反応を解析の対象にする場合には、上記で挙げた主鎖の他、当該螺旋構造の周期構造、又はこれらの組み合わせを特徴的な部分構造として採用することもできる。螺旋構造については、全体を1つの特徴的な部分構造としてもよいし、螺旋構造を所定の周期で分割して得られた部分的な構造のそれぞれを特徴的な部分構造としてもよい。
【0113】
「非中継座標」は、質点系の構造を記述しかつ中継座標から独立した座標である。非中継座標も、中継座標と同様に、質点系に含まれる一部の質点の座標そのもの、質点の座標を組み合わせて定義できる座標又はこれらの組み合わせを用いて設定される。質点座標の組み合わせによっては、非中継座標の構成の一部を中継座標で表現することもできる。ただし、非中継座標は、中継座標から独立した座標であるため、中継座標のみによって表現されることはない。例えば非中継座標は、一部の質点の座標が中継座標として設定されている場合には、N個の質点のうちその質点を除く質点座標とすることができる。また、非中継座標は、質点座標を組み合わせて定義できる座標のみが中継座標として設定されている場合には、そのすべてまたは一部を、質点座標そのもの或いは質点座標を組み合わせて定義できる座標で、かつ中継座標とは独立な座標とすることができる。
【0114】
中継座標の具体的な設定方法の2つの例について説明する。
【0115】
1つは、多原子系に含まれるN個の原子を特徴的な部分構造ごとに集団分けして原子の集団ごとに重心を抽出し、これらの重心の座標のみを中継座標として設定する方法である。この場合には、非中継座標は、その原子が属する部分構造の重心座標に対する相対座標でかつ独立な座標とすることが好ましい。したがって、非中継座標の個数は3N−M個となる。これは、多原子系の3N個の質点座標をM個の重心座標と3N個の相対座標に分離したとき、相対座標間には定義からM個の関係式が発生するため、つまり相対座標のうちM個は独立な座標ではないためである。この方法は、特徴的な部分構造同士の位置関係が大振幅集団運動による多原子系の構造変化に大きく影響することを考慮したものである。特徴的な部分構造自体の変形の規模が小さく、かつその変形の時間スケールが化学反応の時間スケールに比べて短いため、その変形が多原子系の構造変化に与える影響は小さい。
【0116】
もう1つは、生体高分子の主鎖を構成する原子のそれぞれの座標を中継座標として設定する方法である。この場合には、非中継座標は、その生体高分子の側鎖および溶媒和等のその他の原子の座標とすることが好ましい。したがって、非中継座標の個数は3N−M個となる。この方法は、生体高分子の主鎖の形状が大振幅集団運動による多原子系の構造変化に大きく影響することを考慮したものである。側鎖は主鎖に従属して移動することから、側鎖の変形が多原子系の構造変化に与える影響は小さい。この場合において、結合候補分子等の生体高分子以外の分子が多原子系に含まれている場合には、生体高分子以外の分子については、その分子の重心を中継座標として設定してもよい。この場合には、例えばその分子を構成する原子の座標は非中継座標として設定される。
【0117】
なお、後述する第1の実施形態は前者の設定方法を用いた場合に相当する。
【0118】
中継座標は、必要に応じて、すべての質点を自由な状態(質点系に固定しない状態)にした上で質点系に対して全構造緩和が行われた後に、設定されることが好ましい。「全構造緩和」とは、適当な開始状態から質点系のポテンシャルエネルギーが常に減少していくように、そのポテンシャルエネルギーの勾配がゼロとなる状態まで、質点の座標を移動させることを意味する。全構造緩和はいわゆる一般的な構造緩和と同義であるが、後述する「部分構造緩和」と区別するためにこのように称することにする。この全構造緩和によって見出される質点系の状態は、上記開始状態近傍のローカルミニマムの一つである。つまり、この全構造緩和は、ポテンシャル曲面上の上下を繰り返しながら安定点間を結ぶ反応経路およびミニマムポイント(最小点)を見出す作業ではない。この全構造緩和は、特に、構造緩和する対象を構成する原子のすべてを移動させる点で後述する「部分構造緩和」と区別される。この全構造緩和は、例えば、共役勾配法、最急降下法、逆ヘッセ法等既存の方法を用いて実施される。
【0119】
<座標抽出手段>
座標抽出手段は、中継座標の関数としての非中継座標の構成を求め、集団座標の関数としての中継座標の構成を求める。「中継座標の関数としての非中継座標の構成」とは、中継座標RSで表された非中継座標RFについての関数RF(RS)の具体的な内容を意味し、「集団座標の関数としての中継座標の構成」とは、集団座標q†で表された中継座標RSについての関数RS(q†)の具体的な内容を意味する。
【0120】
(中継座標の関数としての非中継座標の構成)
非中継座標の構成RF(RS)は、非中継座標RFを中継座標RSに従属させること、つまり非中継座標RFおよび中継座標RSの関係を規定する条件式を与えることにより得られる。このような条件式は、特に限定されないが、例えば断熱近似条件から得られる下記式30の条件式を用いることができる。なお、「断熱近似」とは、非中継座標RFが中継座標RSの変化に対し即座に追随することができるとする近似を意味する(M. Born, and J. R. Oppenheimer, Ann. Phys. 84 (1927) 457およびH. Haken, Synergetics 2rd Edition, Springer, 1978参照)。式30は、与えられたRSに対し、RFが、ポテンシャルエネルギーVの局所安定点、つまりポテンシャルエネルギーVのRFに関する勾配がゼロの点にいるという条件を表している。
【0121】
【数16】
【0122】
式30において、V(RS,RF)は、中継座標RSおよび非中継座標RFを用いて表された質点系のポテンシャルエネルギーを表す。V(RS,RF)は、質点座標xで表されたポテンシャルエネルギーV(x)に、中継座標RSが設定された段階で定まる質点座標xの構成x(RS,RF)を代入することにより得ることができる。例えば、質点座標xの構成x(RS,RF)は、中継座標RSが設定された段階で定まる中継座標RSの構成RS(x)および非中継座標RFの構成RF(x)に基づいて、これらの逆関数を求めることにより得られる。
【0123】
なお、すべての非中継座標RF1〜RF(3N−M)を中継座標RSに従属させるため、すべての非中継座標RF1〜RF(3N−M)でV(RS,RF)を偏微分することが必要である。
【0124】
非中継座標RFが中継座標RSに従属した質点系においては、下記式31のように非中継座標RFが中継座標RSとして表されることになり、質点系のポテンシャルエネルギーVは下記式32のように表される。以下、上記のような非中継座標RFが中継座標RSに従属した条件の下、中継座標RSのみの関数として表された質点系のポテンシャルエネルギーVeffを有効ポテンシャルエネルギーともいう。
【0125】
【数17】
【0126】
また、式31のように非中継座標RFが中継座標RSに従属する場合には、質点座標xは下記式33のように中継座標の関数として表すことができる。
【0127】
【数18】
【0128】
(集団座標の関数としての中継座標の構成)
中継座標の構成RS(q†)は、非中継座標RFが中継座標RSに従属するという条件の下、中継座標RSの変化に伴う非中継座標RFの変化が中継座標RSに与える影響を考慮しながら、中継座標RSに対して集団運動理論における基礎方程式を解くことにより求められる。中継座標の構成RS(q†)を求めるに際し重要な点は、「中継座標RSに対して集団運動理論における基礎方程式を解くこと」およびその際に「中継座標RSの変化に伴う非中継座標RFの変化が中継座標RSに与える影響(フィードバック)を考慮すること」である。
【0129】
以下、この2点について説明する。
【0130】
中継座標RSに対して集団運動理論における基礎方程式を解く理由は、前述したように、3N次元の質点座標xから集団座標q†を抽出するのではなく、M次元の中継座標RSから集団座標q†を抽出することにより、変数分離の計算を単純化するためである。
【0131】
SCC法における最も基本的な基礎方程式は下記式34から式36である(非特許文献10)。
【0132】
【数19】
【0133】
miは質点系におけるi番目の中継座標に関する質量を表す。miは、i番目の中継座標が複数の質点の重心座標である場合には、それらの質点の質量の総量となる。中継座標が、質点座標を組み合わせた座標であって重心座標ではない場合、単純に質量の総量にはならない。この場合、通常の解析力学の処方で中継座標に対応する一般運動量を求め、ハミルトニアンを書き下し、そのハミルトニアンの一般運動量が関与する項の係数として、中継座標に関する質量を求める。また、φiμ(RS)、λμ(RS)およびρμ(RS)は式35および式36を満たすRSの関数であり、φiμ(RS)は(i,μ)番目、λμ(RS)およびρμ(RS)はμ番目の成分である。また、V,i(RS)およびV,ij(RS)はそれぞれ下記式37および式38で定義される。
【0134】
【数20】
【0135】
したがって、一般的には式34から式36に従い、M個の座標(RS1,RS2,・・・,RSM)で張られるM次元の超空間の中にK個の集団座標(q1,q2,・・・,qK)でパラメトライズ(媒介変数化)されるK次元の超平面(以下、単に平面という)を構成することで、関数RSi(qμ)が求められる。
【0136】
また、基礎方程式としては非特許文献10における近似計算のための基礎方程式を援用することもできる。このとき、例えば自由度が1(K=1)の集団座標q†の関数としての中継座標の構成RS(q1)を求める場合には、中継座標RSを用いた集団運動理論における基礎方程式は下記式39および式40のように表される。
【0137】
【数21】
【0138】
φi(RS)は、式40を満たす関数(固有ベクトル)のi番目の成分を表し、
Λ(RS)は、式40を満たす関数(固有値)を表す。
【0139】
一方、従来の非特許文献10におけるSCC法では、近似的に計算が行われており、変数分離を高精度に行うことができない場合がある。そこで、本発明者は、中継座標RSに対する変数分離をより正確に計算することができる新たな近似法に基づき基礎方程式を導出し、この基礎方程式に基づく新たなSCC法(SCC2法)の理論を構築した。
【0140】
例えばSCC2法に基づいて、自由度が1(K=1)の集団座標q†の関数としての中継座標の構成RS(q1)を求める場合には、中継座標RSを用いた集団運動理論における基礎方程式は下記式41から式43のように表される。
【0141】
【数22】
【0142】
φi(RS,λ)およびκ(RS,λ)は式43を満たす関数であり、φi(RS,λ)はそのi番目の成分を表し、
λは、補助座標を表す。
【0143】
最終的に、式39または式41の解として、式44で表されるような中継座標RSの構成が求められる。
【0144】
【数23】
【0145】
また、上記ではK=1の場合における基礎方程式について説明したが、本発明はK=1の場合に限られない。具体的には、K≧2の場合におけるSCC法の基礎方程式への式39および式40の拡張は、非特許文献10に記載されているように、フロベニウスの定理を用いて行うことができる。一方、K≧2の場合におけるSCC2法の基礎方程式への式41から式43の拡張も、同様の方法で自明に行うことができる。
【0146】
このように、中継座標という階層化された変数(つまり質点系の構造を二次的に記述するより低次元の変数)を導入することにより、3N次元からK次元への縮約の問題がM次元からK次元への縮約の問題に簡略化されるため、変数分離の計算が容易となる。
【0147】
しかしながら、本発明では、質点系の構造を正確に計算するために、フィードバックを考慮しながら変数分離を実施しなければならない。フィードバックを考慮しないということは、存在する質点(または非中継座標)を存在しないものとみなして計算することと同義であり、この影響を考慮せずに変数分離を行うと計算内に矛盾が生じて計算結果が破錠する問題が起こるためである。そこで、フィードバックを考慮するために、既存の集団運動理論に新たな理論に基づいた修正を加える。この新たな理論とは、非中継座標RFが中継座標RSに従属するという条件の下、集団運動理論の基礎方程式におけるポテンシャルエネルギーの座標に関する微分係数を正確に計算し、さらに、中継座標RS(q†)を微小更新(RS(q†)→RS(q†+Δq†))させるごとに非中継座標RFに対して部分構造緩和を行うというものである。
【0148】
従来の集団運動理論では、ポテンシャルエネルギーV(x)の質点座標xに関する微分係数に対する所定の固有ベクトルφの方向に、質点座標xを微小更新(x(q†)→x(q†+Δq†))することが繰り返されるのみである。
【0149】
一方本発明では、ポテンシャルエネルギーVの座標に関する微分係数について、その計算をする際、下記式45から式47のフィードバック方程式を適用することができる。
【0150】
【数24】
【0151】
式46において、第3項の(i←→j)は、第2項におけるiおよびjの文字を相互に入れ替えた項を表す。
【0152】
式47において、第3項の(i←→k)は、第2項におけるiおよびkの文字を相互に入れ替えた項を表し、第4項の(j←→k)は、第2項におけるjおよびkの文字を相互に入れ替えた項を表し、第6項の(i←→k)は、第5項におけるiおよびkの文字を相互に入れ替えた項を表し、第7項の(j←→k)は、第5項におけるjおよびkの文字を相互に入れ替えた項を表す。
【0153】
さらに、式45から式47において下記式48を用いた。
【0154】
【数25】
【0155】
Kαβ−1(Rs)はKαβ(Rs)の逆行列を表し、
Kαβ(Rs)およびJαi(Rs)はそれぞれ式49および式50で定義される。
【0156】
【数26】
【0157】
式45から式47それぞれにおいて、右辺の第1項以外の項が、中継座標RSの変化に伴う非中継座標RFの変化が中継座標RSに与える影響を表している。このようなフィードバック方程式の構成は従来知られておらず、本発明者が初めて見出したものである。フィードバック方程式は、1階微分、2階微分および3階微分の3つの式を含むが、集団運動理論における基礎方程式の内容に応じて必要なものを使用すればよい。例えば、集団運動理論における基礎方程式として、式39および式40を採用した場合には、フィードバック方程式は式46の2階微分の式のみで充分である。
【0158】
そして、非中継座標RFに対する部分構造緩和は、ポテンシャルエネルギーVの座標に関する微分係数であってフィードバック方程式で規定されるものに対する固有ベクトルφの方向に、中継座標RS(q†)を微小更新(RS(q†)→RS(q†+Δq†))させ、更新された中継座標RS(q†+Δq†)の下で、非中継座標RFが中継座標RSに従属するという条件に従って構造緩和を行うことにより行われる。「中継座標RS(q†+Δq†)の下で」とは、更新された中継座標RS(q†+Δq†)を質点系に固定した状態でという意味である。また、「非中継座標RFが中継座標RSに従属するという条件に従って構造緩和を行う」とは、質点系のポテンシャルエネルギーが常に減少していくように、そのポテンシャルエネルギーの勾配がゼロとなる状態まで、中継座標RSに従属した状態を維持しながら非中継座標RFを移動させることを意味する。上記のように中継座標RSを質点系に固定した状態で行う構造緩和を「部分構造緩和」という。前述したように、この部分構造緩和は、中継座標RSを質点系に固定した状態で構造緩和を行う点で全構造緩和と異なる。ただし、部分構造緩和においても、構造緩和を行う手法としては、全構造緩和と同様に、例えば、共役勾配法、最急降下法、逆ヘッセ法等既存の方法を用いることができる。なお、更新する方向を決める上記固有ベクトルφは、最小固有値に対応するものであることが好ましい。本発明は、質点系の大振幅集団運動を対象としているためである。
【0159】
K=1である場合、部分構造緩和における非中継座標RFの開始状態は、RFα(RS(q1))でもよいが、下記式51で与えられる微小量だけ変化させた状態、つまり下記式52で表される状態RFα†であることが好ましい。
【0160】
【数27】
【0161】
式52で表される座標RFα†は、微小更新が無限小ならば非中継座標RF(RS(q1+Δq1))と一致する。しかしながら、実際の計算では微小更新は有限である。部分構造緩和における非中継座標RFの開始状態をRFα†とすることにより、部分構造緩和を行う上での安全性を高めることができる。
【0162】
(座標抽出手段における処理手順)
座標抽出手段における一連の具体的な処理手順について説明する。
【0163】
まず、座標抽出手段は、中継座標Rsおよび非中継座標RFの関数として表現された質点系のポテンシャルエネルギーV(RS,RF)を求める第1の工程を実施する。
【0164】
次に、座標抽出手段は、非中継座標RFを中継座標RSに従属させる第2の工程を実施する。中継座標RSの関数としての非中継座標RFの構成が得られる。ここでは、例として式30の断熱近似条件に従って非中継座標RFを中継座標RSに従属させた場合を考える。
【0165】
次に、座標抽出手段は、座標設定手段によって設定された状態における中継座標RSおよび非中継座標RFの下で、例えば式46のフィードバック方程式に従って、中継座標RSに関するポテンシャルエネルギーVの微分係数を求める第3の工程を実施する。
【0166】
次に、座標抽出手段は、ポテンシャルエネルギーVの上記微分係数の下で、例えば式39および式40のSCC法におけるK=1の場合の基礎方程式に従って、集団座標qに関する中継座標RS(q)の微分係数を求める第4の工程を実施する。
【0167】
次に、座標抽出手段は、中継座標RS(q)の上記微分係数の下で集団座標qを微小量Δqだけ更新することにより、更新された中継座標RS(q+Δq)を求める第5の工程を実施する。
【0168】
次に、座標抽出手段は、更新された中継座標RS(q+Δq)の下で、上記断熱近似条件に従って中継座標RSに従属した非中継座標RFについて部分構造緩和を行う第6の工程を実施する。この際、式52で表される状態RFα†を部分構造緩和の開始状態とする場合には、第4の工程で得られた中継座標RS(q)の上記微分係数から式51に従って、集団座標qに関する非中継座標RF(q)の微分係数を求め、この微分係数の下で非中継座標を微小量Δqだけ更新させておく。
【0169】
そして、座標抽出手段は、中継座標RSおよび非中継座標RFを交互に更新しながら、第3の工程から第6の工程までを例えば次の極小点に至るまで順次繰り返すことにより、集団座標qの関数としての中継座標RSの構成を求める。つまり、2回目以降の第3の工程においては、直前の第6の工程で行われた非中継座標RFの部分構造緩和後の状態における中継座標RSおよび非中継座標RFの下で、上記第3の工程と同様に、中継座標RSに関するポテンシャルエネルギーVの微分係数を求めることになる。極小点を判定する方法は、特に限定されないが、例えばポテンシャルエネルギーVの上記微分係数がゼロ相当になったか否かを判定する方法を挙げることができる。
【0170】
<逆変換手段>
逆変換手段は、集団座標q†についての運動方程式の解として得られる時間tの関数としての集団座標q†(t)、式44で表される中継座標RSの上記構成RS(q)および式31で表される非中継座標RFの上記構成RF(RS)に基づいて、質点座標xの時間発展を予測する。
【0171】
具体的には以下の通りである。図4は、本発明における変数変換および逆変換の概念を示す図である。図4では、例として質点座標xの一部の座標そのものを中継座標RSとし、この中継座標RSから自由度が1の集団座標q†=q1を抽出した場合について示している。
【0172】
図4aは、質点座標xに基づいて中継座標RSおよび非中継座標RFを設定する様子を示している。図4aでは、M個の質点座標が中継座標RSとして設定され、それ以外の3N−M個の質点座標が非中継座標RFとして設定されている。中継座標RSおよび非中継座標RFの設定は、前述したとおり、座標設定手段において行われる。図4bは、非中継座標RFが中継座標RSの関数として与えられることを示している。図4cは、中継座標RSが集団座標q1の関数として与えられることを示している。非中継座標RFの構成および中継座標RSの構成は、前述したとおり、座標抽出手段において求められる。
【0173】
図4dは、集団座標q1が時間tの関数として与えられることを示している。集団座標q1についての運動方程式を解くことにより求められる。集団座標q1についての運動方程式は、質点座標xについての運動方程式に、中継座標RSが設定された段階で定まる質点座標xの構成x(RS,RF)を代入して、中継座標RSおよび非中継座標RFについての運動方程式を求め、この運動方程式に式31および式44を代入することにより求めることができる。
【0174】
以上の図4aからdまでの関係を用いることにより、すべての質点座標xを時間tの関数として求めることができる(図4e)。つまり、逆変換手段は、時間tの関数としての集団座標q1を求めた後、中継座標RSを介した質点座標xから集団座標q1への変数変換に対してこれを逆変換することにより、時間tの関数としての質点座標x(t)の構成を求めるものである。この関数x(t)は、時間変化に対する質点座標の挙動を表すものであるから、この関数x(t)こそが求めるべき反応経路に相当するものである。したがって、この関数x(t)の時間変化を見ることにより質点座標の時間発展を予測することが可能となる。
【0175】
「本発明の第1の実施形態」
図5は、本実施形態のシミュレーション装置10の構成を概略的に示すブロック図である。また、図6Aおよび図6Bは、本実施形態のシミュレーション方法における計算手順を示すフローチャート図である。なお、本実施形態では、生体高分子およびこれに対する結合候補分子を含むN個の原子から構成される多原子系において、多原子系の構造のうち特徴的な部分構造ごとに抽出された重心の座標をM個の中継座標として設定し、かつ、集団座標の自由度が1(K=1)である場合について説明する。
【0176】
本実施形態のシミュレーション装置10は、生体高分子および結合候補分子から構成される複合体の挙動を予測するシミュレーション装置であって、図5に示されるように、解析対象となる多原子系を規定するのに必要なデータ(入力データ)を入力するための入力手段12と、装置の各部を制御する制御手段14と、座標設定手段16と、座標抽出手段18と、逆変換手段20と、解析結果を表示する表示手段22とを備える。
【0177】
そして、本実施形態のシミュレーション方法は、シミュレーション装置10により実行されるシミュレーション方法であって、入力手段12によって入力データを入力し、座標設定手段16によってM個の中継座標RSおよび3N−M個の非中継座標RFを設定させ、座標抽出手段18によって、K=1の集団座標q1の関数としての中継座標RSの構成および中継座標RSの関数としての非中継座標RFの構成を求めさせ、逆変換手段20によって時間tの関数としての原子の座標x(t)を求めさせ、表示手段22によって解析結果を表示させるものである。
【0178】
また、本実施形態に係るシミュレーションプログラムは、上記に記載のシミュレーション方法をコンピュータに実行させることを特徴とするものである。
【0179】
また、本実施形態に係るコンピュータ読み取り可能な記録媒体は、上記に記載のシミュレーションプログラムを記録したものである。
【0180】
<入力手段>
入力手段12は、ユーザーにより入力データが入力される部分である。入力された入力データは制御手段14へ出力される。入力データの入力方法は、特に限定されず、例えばシミュレーション装置10をユーザーが操作しながら手動で行う方法、所定の記録媒体から読み込む方法等が挙げられる。
【0181】
<制御手段>
制御手段14は、各部とのデータのやり取り等の処理を制御する部分である。入力手段12から制御手段14へ入力された入力データは座標設定手段16に出力される。入力データが座標設定手段16へ出力されることにより、当該入力データに基づく多原子系の挙動の計算が開始される。制御手段14には、各部とのデータのやり取りの中で生成される途中のおよび最終的な計算結果を記録するため、メモリ等の記憶媒体が含まれる。
【0182】
また、制御手段14は、逆変換手段20で得られた原子の座標x(t)に関するデータに基づいて、その結果を図やグラフを用いて表示するよう表示手段22を制御する。
【0183】
<座標設定手段>
座標設定手段16は、制御手段14から座標設定手段16へ入力された入力データに基づいて、複合体の構造を最適化し、その後M個の中継座標RSおよび3N−M個の非中継座標RFを設定するものである。得られた中継座標RSおよび非中継座標RFに関するデータは、制御手段14を介して座標抽出手段18に出力される。
【0184】
<座標抽出手段>
座標抽出手段18は、座標設定手段16で得られた中継座標RSおよび非中継座標RFに関するデータに基づいて、多原子系全体のポテンシャルエネルギーV(RS,RF)を求め、このV(RS,RF)に基づいて集団座標q1の関数としての中継座標RSの構成および中継座標RSの関数としての非中継座標RFの構成を求めるものである。得られた中継座標RSの構成および非中継座標RFの構成に関するデータは、制御手段14を介して逆変換手段20に出力される。
【0185】
<逆変換手段>
逆変換手段20は、座標抽出手段18で得られた中継座標RSの構成および非中継座標RFの構成に関するデータに基づいて、時間tの関数としての原子の座標x(t)を求めるものである。得られた原子の座標x(t)に関するデータは制御手段14に出力される。
【0186】
<表示手段>
表示手段22は、制御手段14の指示に従い、例えば多原子系のポテンシャルエネルギーVをグラフで表示したり、複合体の結合過程における構造および/または最終状態の構造を画像として表示したり、複合体の結合過程を動画として表示したりするものである。表示方法としては特に限定されず、2D表示、3D表示等公知の方法を用いることが可能である。
【0187】
以下、図6Aおよび図6Bを参照しながら、本実施形態の装置10を用いたシミュレーション方法の手順を詳細に説明する。
【0188】
<STEP1>
まず、入力手段12から入力データがシミュレーション装置10に入力される(図6AのST1)。入力データの具体的な内容は、生体高分子およびこれに対する結合候補分子を含むN個の原子から構成される多原子系を構成する原子の種類、個数、初期配置の情報となる座標xおよび座標xからポテンシャルエネルギーV(x)を構成するルール等である。
【0189】
<STEP2>
次に、複合体を含む多原子系の構造が最適化されて、初期状態における構造が求められる(図6AのST2)。
【0190】
具体的には以下の通りである。まず、座標設定手段16は、生体高分子および結合候補分子に対してそれぞれ独立に全構造緩和を実施して、生体高分子および結合候補分子の構造をそれぞれ最適化する。そして、それぞれ構造が最適化された生体高分子および結合候補分子は、互いに相互作用を開始する程度(例えば、結合候補分子が生体高分子の表面にわずかに接触する程度)の距離に配置される。この状態が、生体高分子および結合候補分子からなる複合体の初期状態となる。
【0191】
<STEP3>
次に、STEP2で得られた初期状態における構造に基づいて、中継座標RSおよび非中継座標RFが設定される(図6AのST3)。
【0192】
座標設定手段16は、生体高分子の特徴的な部分構造ごとに(M−3)/3個の部分構造に分け、結合候補分子全体を1個の部分構造として取り扱い、M/3個の部分構造のそれぞれから重心を抽出する。そして、このM/3個の重心の座標がM個の中継座標RSとして設定され、M/3個の各部分構造に含まれる原子の各重心からの相対座標でかつ独立な座標が3N−M個の非中継座標RFとして設定される。本実施形態では、非中継座標RFとして、その原子の属する部分構造の重心座標からの相対座標が設定されるとする。
【0193】
<STEP4>
次に、中継座標RSおよび非中継座標RFで表された多原子系のポテンシャルエネルギーV(RS,RF)が求められる(図6AのST4)。
【0194】
図7は、本実施形態における中継座標RS、非中継座標RFおよび原子の座標xの関係を示す概略図である。なお、図7において、m(mは1≦m≦M/3を満たす整数)番目の特徴的な部分構造に含まれる原子数をNmと表した。つまり、N1+N2+・・・+Nm+・・・+NM/3=Nである。本発明では、Mは3次元空間における原子(質点)の座標の個数であるため、通常M/3は整数となる。また、本実施形態においては、中継座標RS、非中継座標RFおよび原子の座標xは、図7から具体的に下記式53および式54で表されるような関係を有していることがわかる。
【0195】
【数28】
【0196】
式53は、m番目の中継座標RS、つまりm番目の特徴的な部分構造の重心座標と原子の座標xとの関係を示した式である。式54は、m番目の特徴的な部分構造に属するn番目の原子の座標x、中継座標RSおよび非中継座標RFの関係を示した式である。式53および式54において、nは、特徴的な部分構造のそれぞれに含まれる原子の通し番号である。したがって、m番目の特徴的な部分構造についてのnは、1からNmまでの整数値をとる。また、3N1+3N2+・・・+3Nm−1はm=1のときゼロとする。
【0197】
したがって、ポテンシャルエネルギーV(RS,RF)は、ポテンシャルエネルギーV(x)に対して式54を適用することにより得ることができる。
【0198】
<STEP5>
次に、ポテンシャルエネルギーV(RS,RF)を用いた式30の断熱近似条件に従って、非中継座標RFが中継座標RSに従属させられる(図6AのST5)。本実施形態では、すべての非中継座標RF1〜RF3NでV(RS,RF)を偏微分することになる。
【0199】
<STEP6>
次に、現状の中継座標RSおよび非中継座標RFの下で、式46のフィードバック方程式に従い、中継座標RSに関するポテンシャルエネルギーVの2階微分係数が求められる(図6AのST6)。
【0200】
<STEP7>
次に、STEP6で得られた2階微分係数の下で、式39および式40のSCC法の基礎方程式に従い、集団座標q1に関する中継座標RSの微分係数dRS/dq1が求められる(図6AのST7)。
【0201】
<STEP8>
次に、STEP7で得られた微分係数dRS/dq1の下で、式51に従い集団座標q1に関する非中継座標RFの微分係数dRF/dq1が求められる(図6BのST8)。
【0202】
<STEP9>
次に、STEP7で得られた微分係数dRS/dq1およびSTEP8で得られた微分係数dRF/dq1の下で、集団座標q1が微小量Δq1だけ更新される(図6BのST9)。
【0203】
<STEP10>
次に、微小量Δq1だけ更新された中継座標RS(q1+Δq1)および更新された非中継座標RF(q1+Δq1)が求められる(図6BのST10)。
【0204】
例えば、STEP2における初期状態を集団座標q1の初期条件と設定した(つまり、q1=0のときに多原子系が上記初期状態をとるものとした)場合には、中継座標RS(q1+Δq1)および非中継座標RF(q1+Δq1)は下記式55のように表せられる。
【0205】
【数29】
【0206】
また、STEP13の判断結果に応じてSTEP6からSTEP12までを順次繰り返し実施する場合において、2回目以降のSTEP10で順次下記式56のように表せられるRSi(2Δq1)、RSi(3Δq1)、・・・が得られることになる。ただし、微分係数dRS/dq1は各回で同じ値になるとは限らない。上記の作業を繰り返すことにより、最終的に集団座標q1の関数としての中継座標RSの構成RS(q1)が求められる。
【0207】
【数30】
【0208】
<STEP11>
STEP11では、STEP10で更新された中継座標RS(q1+Δq1)の下で、STEP10で更新された非中継座標RF(q1+Δq1)を開始状態として非中継座標RFについて部分構造緩和が行われる(図6BのST11)。
【0209】
<STEP12>
次にSTEP12において、部分構造緩和後の状態における多原子系のポテンシャルエネルギーVが計算される(図6BのST12)。なお、STEP13の判断結果に応じて、繰り返しの各回でポテンシャルエネルギーVが計算される。ポテンシャルエネルギーの数値計算は、有効ポテンシャルエネルギーVeffに基づいた計算でもよい。
【0210】
<STEP13>
次にSTEP13において、STEP12で計算されたポテンシャルエネルギーVが極小点に至ったか否かが判断される(図6BのST13)。生体高分子および結合候補分子の化学反応を考えた場合、化学反応の終了点はポテンシャルエネルギーVが次に極小となるとき、つまりポテンシャル曲面上の次の安定点である。そこで、ポテンシャルエネルギーVの値を基準に化学反応の終了点が判断される。ポテンシャルエネルギーVが次の極小点に至ったと判断されればSTEP14へ進み、そうでなければ再度STEP6からSTEP12までが繰り返される。また、ポテンシャルエネルギーVが極小点に至ったと判断された場合であっても、必要に応じて、ポテンシャルエネルギーVがさらに次の極小点に至るまで解析を進めることも可能である。
【0211】
<STEP14>
STEP13においてポテンシャルエネルギーVが極小点に至ったと判断された場合、STEP14において、そのポテンシャルエネルギーVにおける多原子系の構造が、求めるべき最終状態の構造となる(図6BのST14)。
【0212】
そして、最終状態の構造が得られるとともに、式44で表される中継座標RSの構成RS(q1)および式31で表される非中継座標RFの構成RF(RS)が求められる。その後、時間tの関数としての座標x(t)が求められる。
【0213】
<作用効果>
図8は、本実施形態のシミュレーション方法によって得られた、所定のタンパク質2およびこれに対する生理活性物質6(薬物)の複合体の反応経路を示すグラフである。また、図9は、図8のグラフにおける複合体の結合過程を示す概略図である。図9aは複合体の結合過程における初期状態を示し、図9bは複合体の結合過程における最終状態を示す。
【0214】
図8のグラフにおいて、横軸は所定のタンパク質2の残基2aと残基2bの距離であり、縦軸は薬物6および上記タンパク質2の残基2cの距離である。残基2aおよび残基2bは、タンパク質2のポケット4の開口にある残基であり、残基2cは、タンパク質2および薬物6が結合した時に薬物6に対向する位置にある残基である(図9)。また、グラフ中のAおよびCが、例えば図2における安定点AおよびCにそれぞれ対応し、グラフ中のBが、例えば図2におけるエネルギー障壁Bに対応している。グラフ中のプロットは時間軸上で等間隔な点を示しており、このグラフからエネルギー障壁Bに至るまでは反応がゆっくり進み、エネルギー障壁Bを超えた後は反応が速やかに進む様子がわかる。本実施形態のシミュレーション方法により、2000個程度の原子を含む多原子系においてこの反応経路を求めるのに要した時間はおよそ30分である。
【0215】
また、図8のグラフ中の「X線解析」の位置は、上記タンパク質2および上記薬物6についてX線回折測定による結晶構造解析の結果実際に得られた実測値を示す。これにより、本発明は非経験的にシミュレーションを行っているにもかかわらず、非常に高い精度でこの複合体の結合状態を計算できていることがわかる。
【0216】
以上のように、本実施形態に係るシミュレーション装置は、上記のような変数設定手段と、変数抽出手段と、逆変換手段とを備え、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0217】
本実施形態に係るシミュレーション方法は、上記シミュレーション装置によって実行される方法であって、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものである。したがって、集団座標の抽出によりシミュレーションを原子レベルで実行することができ、かつ中継座標の導入により当該集団座標を抽出する際に取り扱う座標の個数を減少させることができる。この結果、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0218】
本実施形態に係るプログラムおよび記録媒体は、上記シミュレーション方法を実行させることが可能であるから、質点系の動力学的な挙動を予測するシミュレーションにおいて、計算精度を向上させかつ計算時間を短縮することが可能となる。
【0219】
「本発明の第2の実施形態」
次に第2の実施形態について説明する。本実施形態は、非特許文献10におけるSCC法を援用した基礎方程式(式39および式40)またはSCC2法の基礎方程式(式41から式43)の解を求める上で、これらの基礎方程式とそれぞれ同値の方程式を解く点で、第1の実施形態と異なる。したがって、第1の実施形態と同様な構成要素の詳細な説明は必要のない限り省略する。
【0220】
本実施形態のシミュレーション装置10も、第1の実施形態と同様に、入力手段12と、制御手段14と、座標設定手段16と、座標抽出手段18と、逆変換手段20と、解析結果を表示する表示手段22とを備える。
【0221】
そして、本実施形態では座標抽出手段18は、式39および式40と同じ解を得るために、SCC法における基礎方程式として下記式57を使用する。
【0222】
【数31】
【0223】
Yは、下記式58で定義されるMK+M+K次元のベクトルである。つまり、Yは、iおよびμのそれぞれを走らせた各要素、m11/2RS1,・・・,mM1/2RSMと、φ11,・・・,φM1,・・・,φ1K,・・・,φMKと、Λ1,・・・,ΛKとを有する。なおφiμおよびΛμは補助座標を表す。
【0224】
【数32】
【0225】
vμは、下記式59の非斉次線形方程式の解ベクトルである。なお、vμはMK+M+K次元のベクトルである。
【0226】
【数33】
【0227】
式59におけるCおよびsμはそれぞれ式60および式61で定義される。なお、Cは(MK+M+K)×(MK+M+K)の正方行列である。Cは一見3×3の行列に見えるが、各成分に更に行列を有する。つまり、3×3の行列として表現されたCのうち、1行1列の成分はMK×M行列(つまり行数がMKで列数がM)であり、1行2列の成分はMK×MK行列であり、1行3列の成分はMK×K行列であり、2行2列の成分はK×MK行列であり、3行1列の成分はM×M行列であり、その他の成分はゼロ行列である。各成分の要素は、iあるいはμを行の足としてi=1〜Mとμ=1〜Kで走らせ、jあるいはνを列の足としてj=1〜Mとν=1〜Kで走らせて構成する。また、sμはMK+M+K次元のベクトルである。sμも、Cと同様に一見3次元のベクトルに見えるが、各成分に更にベクトルを有する。つまり、sμの3次元ベクトルとしての第3成分はM次元ベクトルであり、その他の成分はゼロベクトルである。第3成分の要素は、iをi=1〜Mで走らせて構成する。
【0228】
【数34】
【0229】
V,ijk(RS)は式62で定義される。
【0230】
【数35】
【0231】
上記式57はKが一般化された場合の基礎方程式である。そして、K=1の場合には、式57は下記式63と同値となる。また、上記非特許文献10に記載されているように、フロベニウスの定理を用いてK=1の場合の式63から、一般化されたKについての式57を導出することも可能である。
【0232】
【数36】
【0233】
式63において(χ1,・・・,χM,ψ1,・・・,ψM,κ)は、式64で表される非斉次線形方程式の2M+1次元の解ベクトルである。
【0234】
【数37】
【0235】
なお式64において、行列と固有ベクトルとの積を取り3つの方程式で表現したときにjおよびkについて1からMの総和をとる。
【0236】
式63および式64が式39および式40と同値であることは以下のようにして説明することができる。式38、式62および式63に注意して式64を(φi1,Λ1)に関してq1で積分すると式65が得られる。
【0237】
【数38】
【0238】
さらに、式64に含まれる恒等式χi=φi1を式63の第1式に代入すると式66が得られる。
【0239】
【数39】
【0240】
ここで、式65の第1式および式66におけるφi1およびΛ1をそれぞれφi(RS)およびΛ(RS)とおく(つまり、φi1およびΛ1をRSの関数として取り扱う)と、式39および式40が再現される。
【0241】
また、座標抽出手段18は、式41から式43と同じ解を得る場合には、SCC2法における基礎方程式として下記式67を使用する。
【0242】
【数40】
【0243】
Zは、下記式68で定義されるMK+M+2K次元のベクトルである。つまり、Zは、iおよびμのそれぞれを走らせた各要素、m11/2RS1,・・・,mM1/2RSMと、φ11,・・・,φM1,・・・,φ1K,・・・,φMKと、λ1,・・・,λKと、ρ1,・・・,ρKとを有する。なお、λμおよびρμは、φiμと同様に補助座標を表す。また、cμνは、μごとに定義される下記式69で表される値をそれぞれ最小にするように一意に決められる定数である。また、wμは、下記式70で定義された(MK+M+K)×(MK+M+2K)次元の退化行列DのK次元の特異値空間を張るK個のMK+M+2K次元の単位ベクトルである。Dは一見3×4の行列に見えるが、各成分に更に行列を有する。つまり、3×4の行列として表現されたDのうち、1行1列の成分はMK×M行列であり、1行2列の成分はMK×MK行列であり、1行3列の成分はMK×K行列であり、2行1列の成分はM×M行列であり、2行2列の成分はM×MK行列であり、2行4列の成分はM×K行列であり、3行2列の成分はK×MK行列であり、その他の成分はゼロ行列である。各成分の要素は、iあるいはμを行の足としてi=1〜Mとμ=1〜Kで走らせ、jあるいはνを列の足としてj=1〜Mとν=1〜Kで走らせて構成する。
【0244】
【数41】
【0245】
上記式67はKが一般化された場合の基礎方程式である。そして、K=1の場合には、式67は下記式71と同値となる。また、上記非特許文献10に記載されているように、フロベニウスの定理を用いてK=1の場合の式71から、一般化されたKについての式67を導出することも可能である。
【0246】
【数42】
【0247】
式71において(χ1,・・・,χM,ψ1,・・・,ψM,κ,σ)は、式72で表される(2M+1)×(2M+2)の退化行列による斉次線形方程式の2M+2次元の解ベクトル(不定解)である。
【0248】
【数43】
【0249】
なお式72において、行列と固有ベクトルとの積を取り3つの方程式で表現したときにjおよびkについて1からMの総和をとる。つまり、(χ1,・・・,χM,ψ1,・・・,ψM,κ,σ)は、式72の退化行列の1次元の特異値空間に属する単位ベクトルuを用いて、式73で表される。ここで、cは、式74で表せられる値を最小にするように一意に決められる定数である。
【0250】
【数44】
【0251】
式71および式72が式41から式43と同値であることは以下のようにして説明することができる。式38、式62、式71および式73に注意して式72を(φi1,λ1,ρ1)に関してq1で積分すると式75が得られる。
【0252】
【数45】
【0253】
式75の第1式および第2式からφi1およびρ1を消去すると式76が得られる。そして、RSiとλ1をq1の関数と見なして式76をq1で微分し、式71の第1式および第3式を用いると式77が得られる。
【0254】
【数46】
【0255】
ここで、式71の第1式および第3式並びに式77におけるλ1、χi、κおよびVをそれぞれλ、φi(RS,λ)、κ(RS,λ)およびVeff(RS)とおくと、式41から式43が再現される。
【0256】
本実施形態では、例えば式63または式71から分かるように、基礎方程式を解く上で独立した座標として取り扱うq1の関数としての変数(中継座標RSおよび補助座標の全体であり、以下独立した変数という。)の数が増えている。具体的には、式39および式40では独立した変数は、M個のRSi(i=1からM)のみであったのに対し、式39および式40と同値の式63では独立した変数は、M個のRSi(i=1からM)と、M個のφi1(i=1からM)と、Λ1とから構成され、その総数は2M+1個である。一方、式41から式43では独立した変数は、M個のRSi(i=1からM)およびλであったのに対し、式41から式43と同値の式71では独立した変数は、M個のRSi(i=1からM)と、M個のφi1(i=1からM)と、λ1と、ρ1とから構成され、その総数は2M+2個である。
【0257】
このように独立した変数(具体的には補助座標)の数を増やして方程式を同値変形することにより得られる利点は以下の通りである。
【0258】
例えば式40からφi(RS)を求める場合、式40がφi(RS)の符号反転に対して不変であることを反映し、φi(RS)の符号に任意性が残る。またこのことは、式43からφi(RS,λ)およびκ(RS,λ)を求める場合にも同じである。この符号の任意性は、系の状態が反応経路を逆走する経路もまた反応経路を与えることを表す。一度逆走しても再度逆走することにより順方向に戻れば、最終的に反応経路を見出すことができることを考慮すれば、この符号の任意性はむしろ自然な現象を表現していると言える。しかしながら、この符号の任意性は、上記逆走によって反応経路を効率よく構成することが難しくなるという問題を引き起こす場合がある。さらに、式40および式43は実用上差分化して解かれるが、その差分化の程度Δq1が有限値であるため、逆走が生じても系の状態が元の反応経路上に戻らない場合があり、このような逆走は誤った反応経路が構成されてしまうきっかけとなるという問題もある。
【0259】
そこで、本実施形態では、独立した変数の数を増やすことによりφiμ(RS)の符号の任意性を解消している。その結果、上記逆走の問題が解決され、より効率よくかつ正確に反応経路を構成することが可能となる。なお、独立した変数の数が2倍程度に増えたとしても、計算処理増加による負担はほとんど生じない。
【0260】
図10は、符号の任意性が解消された計算の場合と符号の任意性が残ったままの計算の場合との比較結果を示すグラフである。具体的には、符号24のグラフは、式71から式74を使用して、メリチンというタンパク質の構造変化ごとのポテンシャルエネルギーを計算した結果を示し、符号26のグラフは、式41から式43を使用して、同じタンパク質の構造変化ごとのポテンシャルエネルギーを計算した結果を示す。なお、横軸は、構造変化を見やすくするために導入されたメリチンの屈曲角θであり、縦軸は系のポテンシャルエネルギーである。θ=140degreeの状態が初期状態に相当する。
【0261】
符号の任意性が残ったままの計算では、θ=95degree程度でエネルギーが急上昇しており、この点で一度逆走してしまった。一方、符号の任意性が解消された計算ではそのような現象は見られず、θ=80degree程度でエネルギーの極大点に至った。グラフ24におけるθ=60degreeの離れた点は、上記極大点における状態から系に全体構造緩和を適用して得られたもう1つの局所安定点である。
【0262】
以上のように、本実施形態に係るシミュレーション装置およびシミュレーション方法も、階層化された中継座標を導入し、中継座標の変化に伴う非中継座標の変化が中継座標に与える影響を考慮しながら、集団的かつ本質的な質点系の挙動を記述する集団運動理論における集団座標を抽出し、この集団座標についての運動方程式を解くことにより質点座標の時間発展を予測するものであるから、第1の実施形態と同様の効果を奏する。
【0263】
<設計変更>
上記各実施形態の基礎方程式において、ポテンシャルエネルギーVeff(RS)の3次微分の項が登場する(例えば、式43、式60、式64、式70および式72)。この3次微分の計算処理に要する時間はオーダーでM3(Mは中継座標の個数)に比例するため、その2次微分の計算処理に要する時間がオーダーでM2に比例することを考慮すると、この3次微分の計算はシミュレーションを行う上で非常に大きな負担である。したがって、例えば式64におけるV,ijk(RS)・φk1について、V,ijk(RS)を先に計算しその計算値とφk1の積をとると多くの時間を要してしまう。そこで、本発明者は下記式78を開発し、2つの2次微分との差分をとることで直接V,ijk(RS)・φk1を計算できるようにした。式78を適用することにより、オーダーでM3程度となる計算時間をM2程度にまで短縮することが可能となる。さらに、式78の右辺の2次微分に対して式46のフィードバック方程式を適用することが好ましい。
【0264】
【数47】
【0265】
φiは任意のベクトルのi番目の成分を表し、nは式79で定義される。
【0266】
【数48】
【0267】
(SCC法の基礎方程式について)
最後に、SCC法における最も基本的な基礎方程式(式34から式36)と本明細書で例示した4組の基礎方程式(式39および式40の組、式41から式43の組、式57並びに式67)との関係について簡単に述べておく。
【0268】
前述したように一般的には、最も基本的な基礎方程式に従って平面を構成すれば、関数RSi(qμ)が求められる。しかしながら、上記平面を構成する際に実は2つの問題がある。1つは、式34から式36を同時に厳密に満たす平面が一般のポテンシャルエネルギーVに対して存在するとは限らないという原理的な問題である。そして、もう1つは、式34から式36に従い平面を構成しようとすると逆走が生じる場合があるという実用上の問題(前述した逆走の問題)である。特に、後者の問題の原因は明白で、式35と式36がφiμ(RS)とρμ(RS)の符号を決定することができない、つまり符号反転に対して不変であるからである。
【0269】
そこで、この原理的な問題に対しては、式34から式36を満たす平面を近似的に見出す方法が有効である。実際に非特許文献10では、式36をあきらめ、式34と式35のみで平面を構成しようとした。そしてこの策は、K=1の場合には式39および式40を与える。また、本発明者は、式34を最大限満たしつつ、式35と式36で平面を構成する方法を開発した。そしてこの策は、K=1の場合には式41から式43を与える。具体的には、式35と式36においてφi1(RS)とρ1(RS)を消去することで得られる下記式80を、RSの他、λ1もq1の関数とみなしてq1で微分すると式41から式43が得られる。
【0270】
【数49】
【0271】
一方、実用上の問題に対しては、式34から式36においてφiμ、λμおよびρμをRSから独立しかつq†の関数としての変数(すなわち補助座標)とみなす方法が有効である。これは、φiμ、λμおよびρμをRSから独立した変数とみなすと、φiμ、λμおよびρμは、RSと同様に、平面を構成する際に少しずつ連続的に変化する量となり、突然符号が変わるような事態が発生しないからである。φiμ、λμおよびρμを補助座標として取り扱うためには、φiμ、λμおよびρμについても式34のようなqνでの微分を含む式が必要である。そこで、式35および式36をそれぞれqνで全微分し、これをdqνで割ることにより、下記式81および式82が得られる。これにより、最も基本的な基礎方程式は、式34、式81および式82からなる実用上の問題のない基礎方程式へと変換されたことになる。
【0272】
【数50】
【0273】
そして、式34、式81および式82からなる基礎方程式における原理的な問題に対処することを考える。その対処方法は、例えば上記と同様に、式34、式81および式82を満たす平面を近似的に見出す方法である。具体的には、式34および式81のみで平面を構成する策は式57を与える。一方、式34を最大限満たしつつ、式81と式82で平面を構成する策は式67を与える。
【産業上の利用可能性】
【0274】
各実施形態では、本発明をタンパク質および薬物からなる複合体の結合過程に適用した場合について説明したが、本発明はこれに限られるものではない。つまり、本発明は、より一般的に生体高分子および結合候補分子からなる複合体の形成過程について適用することが可能である。また、本発明は、複合体の形成過程の他、複合体の解離過程、タンパク質のapo体の構造および生体高分子のフォールディング(折り畳み機構)等、生体高分子を含む多原子系の挙動の解析についても適用することが可能である。
【符号の説明】
【0275】
2 タンパク質
4 ポケット
6 結合候補分子
10 シミュレーション装置
12 入力手段
14 制御手段
16 座標設定手段
18 座標抽出手段
20 逆変換手段
22 表示手段
A 複合体が解離体である状態を示す点
B エネルギー障壁
C 結合が完了した複合体である状態を示す点
RP 複合体の挙動の反応経路
【特許請求の範囲】
【請求項1】
モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション装置において、
前記質点系の構造を記述する3N個の質点座標に基づいて、前記質点系の構造の変化を主として担うM個の座標である中継座標を設定し、前記質点系の構造を記述しかつ前記中継座標から独立した座標である非中継座標を設定する座標設定手段と、
前記非中継座標を前記中継座標に従属させることにより前記中継座標の関数としての前記非中継座標の構成を求め、前記中継座標の変化に伴う前記非中継座標の変化が前記中継座標に与える影響を考慮しながら、前記中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される前記一般座標の前記可変成分であるK個の集団座標の関数としての前記中継座標の構成を求める座標抽出手段と、
前記集団座標についての運動方程式の解として得られる時間の関数としての前記集団座標、前記中継座標の前記構成および前記非中継座標の前記構成に基づいて、前記質点座標の時間発展を予測する逆変換手段とを備えることを特徴とするシミュレーション装置。
(K、MおよびNは、K<M<3Nの関係を満たし、それぞれ1以上の整数を表す。)
【請求項2】
前記座標抽出手段が、
前記中継座標および前記非中継座標の関数として表現された前記質点系のポテンシャルエネルギーを求める第1の工程を実施し、
前記ポテンシャルエネルギーを用いた断熱近似条件に従って前記非中継座標を前記中継座標に従属させる第2の工程を実施し、
現状の前記中継座標および前記非中継座標の下で、前記影響を考慮しながら、前記中継座標に関する前記ポテンシャルエネルギーの微分係数を求める第3の工程を実施し、
前記ポテンシャルエネルギーの前記微分係数の下で、自己無撞着集団座標法における前記ポテンシャルエネルギーの微分係数を用いた基礎方程式に従って、前記集団座標に関する前記中継座標の微分係数を求める第4の工程を実施し、
前記中継座標の前記微分係数の下で前記集団座標を微小量だけ更新することにより、更新された前記中継座標を求める第5の工程を実施し、
前記更新された前記中継座標の下で、前記中継座標に従属した前記非中継座標について構造緩和を行う第6の工程を実施し、
その後、前記非中継座標の構造緩和後の状態における前記中継座標および前記非中継座標に基づいて、前記第3の工程から前記第6の工程までを順次繰り返すことにより、前記中継座標の前記構成を求めるものであることを特徴とする請求項1に記載のシミュレーション装置。
【請求項3】
前記影響を考慮する方法が下記式1から式3の少なくとも1つを使用することであることを特徴とする請求項2に記載のシミュレーション装置。
【数1】
(i、jおよびkは、それぞれ1からMまでの整数を表し、
α、βおよびγは、それぞれ1から3N−Mまでの整数を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RFαは、前記質点系におけるα番目の前記非中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
RFは(RF1,RF2,・・・,RF(3N−M))を表し、
RF(RS)は、前記中継座標に従属させられた前記非中継座標を表し、
V(RS,RF)は、前記中継座標と前記非中継座標で表された前記ポテンシャルエネルギーを表し、
Veff(RS)は、V(RS,RF)にRF(RS)を代入して得られる有効ポテンシャルエネルギーを表す。
式2において、第3項の(i←→j)は、第2項におけるiおよびjの文字を相互に入れ替えた項を表す。
式3において、第3項の(i←→k)は、第2項におけるiおよびkの文字を相互に入れ替えた項を表し、第4項の(j←→k)は、第2項におけるjおよびkの文字を相互に入れ替えた項を表し、第6項の(i←→k)は、第5項におけるiおよびkの文字を相互に入れ替えた項を表し、第7項の(j←→k)は、第5項におけるjおよびkの文字を相互に入れ替えた項を表す。
さらに、式1から式3において下記式4を用いた。
【数2】
Kαβ−1(Rs)はKαβ(Rs)の逆行列を表し、
Kαβ(Rs)およびJαi(Rs)はそれぞれ式5および式6で定義される。
【数3】
)
【請求項4】
前記断熱近似条件が下記式7であることを特徴とする請求項2または3に記載のシミュレーション装置。
【数4】
(RFαは、前記質点系におけるα番目の前記非中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
RFは(RF1,RF2,・・・,RF(3N−M))を表し、
V(RS,RF)は、前記中継座標と前記非中継座標で表された前記ポテンシャルエネルギーを表す。)
【請求項5】
前記集団座標の個数KがK=1を満たすものであり、
前記基礎方程式が下記式8および9であることを特徴とする請求項2から4いずれかに記載のシミュレーション装置。
【数5】
(iおよびjは、それぞれ1からMまでの整数を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
q1は、前記集団座標を表し、
miは前記質点系におけるi番目の前記中継座標に関する質量を表し、
φi(RS)は、式9を満たす関数(固有ベクトル)のi番目の成分を表し、
Λ(RS)は、式9を満たす関数(固有値)を表し、
Veff(RS)は有効ポテンシャルエネルギーを表す。)
【請求項6】
前記集団座標の個数KがK=1を満たすものであり、
前記基礎方程式が下記式10から12であることを特徴とする請求項2から4いずれかに記載のシミュレーション装置。
【数6】
(i、jおよびkは、それぞれ1からMまでの整数を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
q1は、前記集団座標を表し、
miは前記質点系におけるi番目の前記中継座標に関する質量を表し、
φi(RS,λ)およびκ(RS,λ)は、式12を満たす関数のi番目の成分を表し、
λは、補助座標を表し、
Veff(RS)は有効ポテンシャルエネルギーを表す。)
【請求項7】
前記座標抽出手段が、前記第4の工程において、前記基礎方程式中の前記集団座標に関する前記中継座標または補助座標の微分係数が有する符号の任意性を解消するように、該基礎方程式を解く上で前記中継座標から独立しかつ前記集団座標の関数として取り扱う変数の数を増やして、計算を行うものであることを特徴とする請求項2から4いずれかに記載のシミュレーション装置。
【請求項8】
前記座標抽出手段が、前記中継座標から独立しかつ前記集団座標の関数として取り扱う変数の数を増やした結果得られた下記式13で表される基礎方程式に従って計算を行うものであることを特徴とする請求項7に記載のシミュレーション装置。
【数7】
(Yは、下記式14で定義されるMK+M+K次元のベクトルであり、
vμは、下記式15の非斉次線形方程式の解ベクトルである。
【数8】
式15におけるCおよびsμはそれぞれ式16および式17で定義される。
【数9】
V,ij(RS)およびV,ijk(RS)はそれぞれ式18および式19で定義される。
【数10】
また、i、jおよびkはそれぞれ1からMまでの整数を表し、
μおよびνはそれぞれ1からKまでの整数を表し、
qμはμ番目の前記集団座標を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
miは前記質点系におけるi番目の前記中継座標に関する質量を表し、
φiμおよびΛμは、RSから独立した補助座標を表し、
Veff(RS)は有効ポテンシャルエネルギーを表す。)
【請求項9】
前記集団座標の個数KがK=1を満たすことを特徴とする請求項8に記載のシミュレーション装置。
【請求項10】
前記座標抽出手段が、前記中継座標から独立しかつ前記集団座標の関数として取り扱う変数の数を増やした結果得られた下記式20で表される基礎方程式に従って計算を行うものであることを特徴とする請求項7に記載のシミュレーション装置。
【数11】
(Zは、下記式21で定義されるMK+M+2K次元のベクトルであり、
cμνは、μごとに定義される下記式22で表される値をそれぞれ最小にするように一意に決められる定数であり、
wμは、下記式23で定義された行列DのK次元の特異値空間を張るK個のMK+M+2K次元の単位ベクトルである。
【数12】
V,ij(RS)およびV,ijk(RS)はそれぞれ式24および式25で定義される。
【数13】
また、i、jおよびkはそれぞれ1からMまでの整数を表し、
μおよびνはそれぞれ1からKまでの整数を表し、
qμはμ番目の前記集団座標を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
miは前記質点系におけるi番目の前記中継座標に関する質量を表し、
φiμ、λμおよびρμは、RSから独立した補助座標を表し、
Veff(RS)は有効ポテンシャルエネルギーを表す。)
【請求項11】
前記集団座標の個数KがK=1を満たすことを特徴とする請求項10に記載のシミュレーション装置。
【請求項12】
前記座標抽出手段が、前記ポテンシャルエネルギーの3次微分の項の計算を下記式26に基づいて行うものであることを特徴とする請求項6から11いずれかに記載のシミュレーション装置。
【数14】
(φiは任意のベクトルのi番目の成分を表し、
V,ij(RS)、V,ijk(RS)およびnはそれぞれ式27から式29で定義される。
【数15】
)
【請求項13】
前記座標設定手段が、前記質点系の構造のうち特徴的な部分構造ごとに抽出された該部分構造を代表する代表座標を前記中継座標として設定するものであることを特徴とする請求項1から12いずれかに記載のシミュレーション装置。
【請求項14】
前記質点系が生体高分子を含む多原子系であり、
前記部分構造が、前記生体高分子の2次構造、前記生体高分子の構成単位または前記生体高分子の主鎖であり、
前記部分構造についての前記代表座標が、該部分構造を構成する原子の座標そのもの、該原子の前記座標を組み合わせて定義できる座標または前記部分構造の周期間隔であることを特徴とする請求項13に記載のシミュレーション装置。
【請求項15】
前記生体高分子がタンパク質であり、
前記部分構造が前記タンパク質の2次構造であり、
該2次構造についての前記代表座標が、該2次構造を構成する原子群の重心座標または前記2次構造の屈曲角であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項16】
前記2次構造がヘリックス構造、βシート、ターン、ループおよびランダムコイルの少なくとも1つであることを特徴とする請求項15に記載のシミュレーション装置。
【請求項17】
前記生体高分子がタンパク質であり、
前記部分構造が前記タンパク質の残基であり、
該残基についての前記代表座標が、該残基を構成する原子群の重心座標であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項18】
前記生体高分子がタンパク質であり、
前記部分構造が前記タンパク質の主鎖であり、
該主鎖についての前記代表座標が、該主鎖を構成する原子のそれぞれの座標であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項19】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の2次構造であり、
該2次構造についての前記代表座標が、該2次構造を構成する原子群の重心座標または前記2次構造の屈曲角であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項20】
前記2次構造が螺旋構造であることを特徴とする請求項19に記載のシミュレーション装置。
【請求項21】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の残基であり、
該残基についての前記代表座標が、該残基を構成する原子群の重心座標であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項22】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の主鎖であり、
該主鎖についての前記代表座標が、該主鎖を構成する原子のそれぞれの座標であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項23】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の螺旋構造であり、
該螺旋構造についての前記代表座標が、該螺旋構造の周期間隔であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項24】
前記多原子系が、前記生体高分子に対する結合候補分子を含むものであることを特徴とする請求項14から23いずれかに記載のシミュレーション装置。
【請求項25】
請求項1から24いずれかに記載のシミュレーション装置によって実行される、モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション方法であって、
前記質点系の構造を記述する3N個の質点座標に基づいて、前記質点系の構造の変化を主として担うM個の座標である中継座標を設定し、
前記質点系の構造を記述しかつ前記中継座標から独立した座標である非中継座標を設定し、
前記非中継座標を前記中継座標に従属させることにより前記中継座標の関数としての前記非中継座標の構成を求め、
前記中継座標の変化に伴う前記非中継座標の変化が前記中継座標に与える影響を考慮しながら、前記中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される前記一般座標の前記可変成分であるK個の集団座標の関数としての前記中継座標の構成を求め、
前記集団座標についての運動方程式の解として得られる時間の関数としての前記集団座標、前記中継座標の前記構成および前記非中継座標の前記構成に基づいて、前記質点座標の時間発展を予測することを特徴とするシミュレーション方法。
(K、MおよびNは、K<M<3Nの関係を満たし、それぞれ1以上の整数を表す。)
【請求項26】
請求項25に記載のシミュレーション方法をコンピュータに実行させためのシミュレーションプログラム。
【請求項27】
請求項26に記載のシミュレーションプログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項1】
モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション装置において、
前記質点系の構造を記述する3N個の質点座標に基づいて、前記質点系の構造の変化を主として担うM個の座標である中継座標を設定し、前記質点系の構造を記述しかつ前記中継座標から独立した座標である非中継座標を設定する座標設定手段と、
前記非中継座標を前記中継座標に従属させることにより前記中継座標の関数としての前記非中継座標の構成を求め、前記中継座標の変化に伴う前記非中継座標の変化が前記中継座標に与える影響を考慮しながら、前記中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される前記一般座標の前記可変成分であるK個の集団座標の関数としての前記中継座標の構成を求める座標抽出手段と、
前記集団座標についての運動方程式の解として得られる時間の関数としての前記集団座標、前記中継座標の前記構成および前記非中継座標の前記構成に基づいて、前記質点座標の時間発展を予測する逆変換手段とを備えることを特徴とするシミュレーション装置。
(K、MおよびNは、K<M<3Nの関係を満たし、それぞれ1以上の整数を表す。)
【請求項2】
前記座標抽出手段が、
前記中継座標および前記非中継座標の関数として表現された前記質点系のポテンシャルエネルギーを求める第1の工程を実施し、
前記ポテンシャルエネルギーを用いた断熱近似条件に従って前記非中継座標を前記中継座標に従属させる第2の工程を実施し、
現状の前記中継座標および前記非中継座標の下で、前記影響を考慮しながら、前記中継座標に関する前記ポテンシャルエネルギーの微分係数を求める第3の工程を実施し、
前記ポテンシャルエネルギーの前記微分係数の下で、自己無撞着集団座標法における前記ポテンシャルエネルギーの微分係数を用いた基礎方程式に従って、前記集団座標に関する前記中継座標の微分係数を求める第4の工程を実施し、
前記中継座標の前記微分係数の下で前記集団座標を微小量だけ更新することにより、更新された前記中継座標を求める第5の工程を実施し、
前記更新された前記中継座標の下で、前記中継座標に従属した前記非中継座標について構造緩和を行う第6の工程を実施し、
その後、前記非中継座標の構造緩和後の状態における前記中継座標および前記非中継座標に基づいて、前記第3の工程から前記第6の工程までを順次繰り返すことにより、前記中継座標の前記構成を求めるものであることを特徴とする請求項1に記載のシミュレーション装置。
【請求項3】
前記影響を考慮する方法が下記式1から式3の少なくとも1つを使用することであることを特徴とする請求項2に記載のシミュレーション装置。
【数1】
(i、jおよびkは、それぞれ1からMまでの整数を表し、
α、βおよびγは、それぞれ1から3N−Mまでの整数を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RFαは、前記質点系におけるα番目の前記非中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
RFは(RF1,RF2,・・・,RF(3N−M))を表し、
RF(RS)は、前記中継座標に従属させられた前記非中継座標を表し、
V(RS,RF)は、前記中継座標と前記非中継座標で表された前記ポテンシャルエネルギーを表し、
Veff(RS)は、V(RS,RF)にRF(RS)を代入して得られる有効ポテンシャルエネルギーを表す。
式2において、第3項の(i←→j)は、第2項におけるiおよびjの文字を相互に入れ替えた項を表す。
式3において、第3項の(i←→k)は、第2項におけるiおよびkの文字を相互に入れ替えた項を表し、第4項の(j←→k)は、第2項におけるjおよびkの文字を相互に入れ替えた項を表し、第6項の(i←→k)は、第5項におけるiおよびkの文字を相互に入れ替えた項を表し、第7項の(j←→k)は、第5項におけるjおよびkの文字を相互に入れ替えた項を表す。
さらに、式1から式3において下記式4を用いた。
【数2】
Kαβ−1(Rs)はKαβ(Rs)の逆行列を表し、
Kαβ(Rs)およびJαi(Rs)はそれぞれ式5および式6で定義される。
【数3】
)
【請求項4】
前記断熱近似条件が下記式7であることを特徴とする請求項2または3に記載のシミュレーション装置。
【数4】
(RFαは、前記質点系におけるα番目の前記非中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
RFは(RF1,RF2,・・・,RF(3N−M))を表し、
V(RS,RF)は、前記中継座標と前記非中継座標で表された前記ポテンシャルエネルギーを表す。)
【請求項5】
前記集団座標の個数KがK=1を満たすものであり、
前記基礎方程式が下記式8および9であることを特徴とする請求項2から4いずれかに記載のシミュレーション装置。
【数5】
(iおよびjは、それぞれ1からMまでの整数を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
q1は、前記集団座標を表し、
miは前記質点系におけるi番目の前記中継座標に関する質量を表し、
φi(RS)は、式9を満たす関数(固有ベクトル)のi番目の成分を表し、
Λ(RS)は、式9を満たす関数(固有値)を表し、
Veff(RS)は有効ポテンシャルエネルギーを表す。)
【請求項6】
前記集団座標の個数KがK=1を満たすものであり、
前記基礎方程式が下記式10から12であることを特徴とする請求項2から4いずれかに記載のシミュレーション装置。
【数6】
(i、jおよびkは、それぞれ1からMまでの整数を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
q1は、前記集団座標を表し、
miは前記質点系におけるi番目の前記中継座標に関する質量を表し、
φi(RS,λ)およびκ(RS,λ)は、式12を満たす関数のi番目の成分を表し、
λは、補助座標を表し、
Veff(RS)は有効ポテンシャルエネルギーを表す。)
【請求項7】
前記座標抽出手段が、前記第4の工程において、前記基礎方程式中の前記集団座標に関する前記中継座標または補助座標の微分係数が有する符号の任意性を解消するように、該基礎方程式を解く上で前記中継座標から独立しかつ前記集団座標の関数として取り扱う変数の数を増やして、計算を行うものであることを特徴とする請求項2から4いずれかに記載のシミュレーション装置。
【請求項8】
前記座標抽出手段が、前記中継座標から独立しかつ前記集団座標の関数として取り扱う変数の数を増やした結果得られた下記式13で表される基礎方程式に従って計算を行うものであることを特徴とする請求項7に記載のシミュレーション装置。
【数7】
(Yは、下記式14で定義されるMK+M+K次元のベクトルであり、
vμは、下記式15の非斉次線形方程式の解ベクトルである。
【数8】
式15におけるCおよびsμはそれぞれ式16および式17で定義される。
【数9】
V,ij(RS)およびV,ijk(RS)はそれぞれ式18および式19で定義される。
【数10】
また、i、jおよびkはそれぞれ1からMまでの整数を表し、
μおよびνはそれぞれ1からKまでの整数を表し、
qμはμ番目の前記集団座標を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
miは前記質点系におけるi番目の前記中継座標に関する質量を表し、
φiμおよびΛμは、RSから独立した補助座標を表し、
Veff(RS)は有効ポテンシャルエネルギーを表す。)
【請求項9】
前記集団座標の個数KがK=1を満たすことを特徴とする請求項8に記載のシミュレーション装置。
【請求項10】
前記座標抽出手段が、前記中継座標から独立しかつ前記集団座標の関数として取り扱う変数の数を増やした結果得られた下記式20で表される基礎方程式に従って計算を行うものであることを特徴とする請求項7に記載のシミュレーション装置。
【数11】
(Zは、下記式21で定義されるMK+M+2K次元のベクトルであり、
cμνは、μごとに定義される下記式22で表される値をそれぞれ最小にするように一意に決められる定数であり、
wμは、下記式23で定義された行列DのK次元の特異値空間を張るK個のMK+M+2K次元の単位ベクトルである。
【数12】
V,ij(RS)およびV,ijk(RS)はそれぞれ式24および式25で定義される。
【数13】
また、i、jおよびkはそれぞれ1からMまでの整数を表し、
μおよびνはそれぞれ1からKまでの整数を表し、
qμはμ番目の前記集団座標を表し、
RSiは、前記質点系におけるi番目の前記中継座標を表し、
RSは(RS1,RS2,・・・,RSM)を表し、
miは前記質点系におけるi番目の前記中継座標に関する質量を表し、
φiμ、λμおよびρμは、RSから独立した補助座標を表し、
Veff(RS)は有効ポテンシャルエネルギーを表す。)
【請求項11】
前記集団座標の個数KがK=1を満たすことを特徴とする請求項10に記載のシミュレーション装置。
【請求項12】
前記座標抽出手段が、前記ポテンシャルエネルギーの3次微分の項の計算を下記式26に基づいて行うものであることを特徴とする請求項6から11いずれかに記載のシミュレーション装置。
【数14】
(φiは任意のベクトルのi番目の成分を表し、
V,ij(RS)、V,ijk(RS)およびnはそれぞれ式27から式29で定義される。
【数15】
)
【請求項13】
前記座標設定手段が、前記質点系の構造のうち特徴的な部分構造ごとに抽出された該部分構造を代表する代表座標を前記中継座標として設定するものであることを特徴とする請求項1から12いずれかに記載のシミュレーション装置。
【請求項14】
前記質点系が生体高分子を含む多原子系であり、
前記部分構造が、前記生体高分子の2次構造、前記生体高分子の構成単位または前記生体高分子の主鎖であり、
前記部分構造についての前記代表座標が、該部分構造を構成する原子の座標そのもの、該原子の前記座標を組み合わせて定義できる座標または前記部分構造の周期間隔であることを特徴とする請求項13に記載のシミュレーション装置。
【請求項15】
前記生体高分子がタンパク質であり、
前記部分構造が前記タンパク質の2次構造であり、
該2次構造についての前記代表座標が、該2次構造を構成する原子群の重心座標または前記2次構造の屈曲角であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項16】
前記2次構造がヘリックス構造、βシート、ターン、ループおよびランダムコイルの少なくとも1つであることを特徴とする請求項15に記載のシミュレーション装置。
【請求項17】
前記生体高分子がタンパク質であり、
前記部分構造が前記タンパク質の残基であり、
該残基についての前記代表座標が、該残基を構成する原子群の重心座標であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項18】
前記生体高分子がタンパク質であり、
前記部分構造が前記タンパク質の主鎖であり、
該主鎖についての前記代表座標が、該主鎖を構成する原子のそれぞれの座標であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項19】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の2次構造であり、
該2次構造についての前記代表座標が、該2次構造を構成する原子群の重心座標または前記2次構造の屈曲角であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項20】
前記2次構造が螺旋構造であることを特徴とする請求項19に記載のシミュレーション装置。
【請求項21】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の残基であり、
該残基についての前記代表座標が、該残基を構成する原子群の重心座標であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項22】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の主鎖であり、
該主鎖についての前記代表座標が、該主鎖を構成する原子のそれぞれの座標であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項23】
前記生体高分子が核酸であり、
前記部分構造が前記核酸の螺旋構造であり、
該螺旋構造についての前記代表座標が、該螺旋構造の周期間隔であることを特徴とする請求項14に記載のシミュレーション装置。
【請求項24】
前記多原子系が、前記生体高分子に対する結合候補分子を含むものであることを特徴とする請求項14から23いずれかに記載のシミュレーション装置。
【請求項25】
請求項1から24いずれかに記載のシミュレーション装置によって実行される、モデル化されたN個の質点から構成された質点系の挙動を予測するシミュレーション方法であって、
前記質点系の構造を記述する3N個の質点座標に基づいて、前記質点系の構造の変化を主として担うM個の座標である中継座標を設定し、
前記質点系の構造を記述しかつ前記中継座標から独立した座標である非中継座標を設定し、
前記非中継座標を前記中継座標に従属させることにより前記中継座標の関数としての前記非中継座標の構成を求め、
前記中継座標の変化に伴う前記非中継座標の変化が前記中継座標に与える影響を考慮しながら、前記中継座標と正準変換により対応付けられる一般座標であって時間に依存して変化する可変成分および時間変化に対して定数となる不変成分から構成される前記一般座標の前記可変成分であるK個の集団座標の関数としての前記中継座標の構成を求め、
前記集団座標についての運動方程式の解として得られる時間の関数としての前記集団座標、前記中継座標の前記構成および前記非中継座標の前記構成に基づいて、前記質点座標の時間発展を予測することを特徴とするシミュレーション方法。
(K、MおよびNは、K<M<3Nの関係を満たし、それぞれ1以上の整数を表す。)
【請求項26】
請求項25に記載のシミュレーション方法をコンピュータに実行させためのシミュレーションプログラム。
【請求項27】
請求項26に記載のシミュレーションプログラムを記録したコンピュータ読み取り可能な記録媒体。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6A】
【図6B】
【図7】
【図8】
【図9】
【図10】
【図2】
【図3】
【図4】
【図5】
【図6A】
【図6B】
【図7】
【図8】
【図9】
【図10】
【公開番号】特開2013−84255(P2013−84255A)
【公開日】平成25年5月9日(2013.5.9)
【国際特許分類】
【出願番号】特願2012−208211(P2012−208211)
【出願日】平成24年9月21日(2012.9.21)
【出願人】(306037311)富士フイルム株式会社 (25,513)
【Fターム(参考)】
【公開日】平成25年5月9日(2013.5.9)
【国際特許分類】
【出願日】平成24年9月21日(2012.9.21)
【出願人】(306037311)富士フイルム株式会社 (25,513)
【Fターム(参考)】
[ Back to top ]