(58)【調査した分野】(Int.Cl.,DB名)
∫Δν(φ)P(φ)dφ計算手段が、第2の相互作用エネルギーφ出現確率計算手段によって計算されたP(φ)と、第1の相互作用エネルギーε出現確率計算手段によって計算されたP0'(ε;φ)と、第2の相互作用エネルギーε出現確率計算手段によって計算されたP'(ε)とに基づいて、相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφをエネルギー表示法により計算することを特徴とする、請求項1又は2に記載の計算装置。
フラグメントBが仮想原子である点電荷を含む原子から構成され、第2の座標取得手段が、フラグメントBの点電荷を、原子集合体Aの構造aを構成する原子の電荷パラメータに付加することを特徴とする、請求項1〜3のいずれか1項に記載の計算装置。
コンピュータの制御部が、∫Δν(φ)P(φ)dφ計算ステップにおいて、第2の相互作用エネルギーφ出現確率計算ステップで計算されたP(φ)と、第1の相互作用エネルギーε出現確率計算ステップで計算されたP0'(ε;φ)と、第2の相互作用エネルギーε出現確率計算ステップで計算されたP'(ε)とに基づいて、相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφをエネルギー表示法により計算することを特徴とする、請求項5又は6に記載の計算方法。
フラグメントBが仮想原子である点電荷を含む原子から構成され、コンピュータの制御部が、第2の座標取得ステップにおいて、フラグメントBの点電荷を、原子集合体Aの構造aを構成する原子の電荷パラメータに付加することを特徴とする、請求項5〜7のいずれか1項に記載の計算方法。
∫Δν(φ)P(φ)dφ計算手段が、第2の相互作用エネルギーφ出現確率計算手段によって計算されたP(φ)と、第1の相互作用エネルギーε出現確率計算手段によって計算されたP0'(ε;φ)と、第2の相互作用エネルギーε出現確率計算手段によって計算されたP'(ε)とに基づいて、相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφをエネルギー表示法により計算することを特徴とする、請求項9又は10に記載のプログラム。
フラグメントBが仮想原子である点電荷を含む原子から構成され、第2の座標取得手段が、フラグメントBの点電荷を、原子集合体Aの構造aを構成する原子の電荷パラメータに付加することを特徴とする、請求項9〜11のいずれか1項に記載のプログラム。
【発明を実施するための形態】
【0028】
一般に、古典統計力学理論において、反応式(1)で表される変化に関し、変化前の原子集合体Aの自由エネルギー及びフラグメントBの自由エネルギーの和と、変化後の原子集合体ABの自由エネルギーとの差は、数式(2)で表される。
【0029】
ここで、変化前の原子集合体A及びフラグメントBは熱平衡状態にあり、変化後の原子集合体ABも熱平衡状態にある。また、以下の数式(2)〜(5)において、積分記号のインテグラル(∫)は、各積分変数に関する積分記号を集合的に表すものとする。
【0030】
一方、変化後の原子集合体ABにおける構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφの出現確率(つまりエネルギー分布)であるP(φ)は、数式(3)で表されるとともに、変化前の原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABにおける構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφの出現確率(つまりエネルギー分布)であるP
0(φ)は、数式(4)で表され、これらの比は数式(5)と整理することができる。数式(1)は、数式(5)を変形することにより導くことができる。ここで、数式(3)及び(4)における、
【数4】
は、ディラックのデルタ関数であって、φが
【数5】
と等しい場合に1であり、その他の場合は零となる関数である。
【0033】
数式(1)の第1項は、変化後の原子集合体ABにおける構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφの平均値を表す。
【0034】
数式(1)の第1項及び第2項は、反応式(1)で表される変化における、フラグメントBと原子集合体Aの構造aとの間に結合が生じることに起因する自由エネルギー変化量を表す。
【0035】
数式(1)の第3項のΔν(φ)は、構造aと該構造aと結合したフラグメントBとの相互作用エネルギーがφである場合の、反応式(1)で表される変化における自由エネルギー変化量であって、フラグメントBと、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除いた原子集合体の一部又は全部との相互作用に起因した自由エネルギー変化量を表す。
【0036】
数式(1)の第2項について、変化後の原子集合体ABにおける構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφの出現確率であるP(φ)及び変化前の原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABにおける構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφの出現確率であるP
0(φ)は、P(φ)とP
0(φ)の比P(φ)/P
0(φ)が発散しないことが必要となる。つまり、P(φ)が、0<P(φ)≦1の範囲で出現確率を有するφでは、P
0(φ)も0<P
0(φ)≦1の出現確率を持つ必要がある。
【0037】
数式(1)の第3項の∫Δν(φ)P(φ)dφは、反応式(1)で表される変化における、相互作用エネルギーεに起因する自由エネルギー変化量を表し、以下の数式(6)を用いて、従来のエネルギー表示法を適用して計算することができる。また、数式(6)左辺の表記は、従来のエネルギー表示法における表記と同様、分布関数を表す。ここで、本発明において、フラグメントBが後述する点電荷等の仮想原子を含む場合、原子集合体ABからフラグメントBを除くとは、例えば、原子集合体ABがアニソール(C
6H
5OCH
3)であり、フラグメントBがメトキシ基(−OCH
3)及び点電荷(該メトキシ基が結合するベンゼン環の炭素原子が有する電荷パラメータを持つ仮想原子)である場合、アニソールからメトキシ基及び該点電荷が持つ電荷パラメータを除去することである。つまり、原子集合体Aは、C
6H
5で表される原子集合体であって、原子集合体ABにおいてメトキシ基が結合しているベンゼン環上の炭素原子が有する電荷パラメータが零である原子集合体である。また、原子集合体AにフラグメントBが結合するとは、例えば、フラグメントBがメトキシ基(−OCH
3)及び点電荷(該メトキシ基が結合するベンゼン環の炭素原子が有する電荷パラメータを持つ仮想原子)であり、原子集合体Aが、フラグメントBのメトキシ基が結合するベンゼン環上の炭素原子が有する部分電荷が零であって、C
6H
5で表される原子集合体である場合、フラグメントBのメトキシ基の酸素原子とベンゼン環上の炭素原子が共有結合を形成し、さらに、フラグメントBの点電荷が持つ部分電荷(電荷パラメータ)をメトキシ基の酸素原子が結合するベンゼン環の炭素原子の電荷パラメータに付加することである。なお、本発明では、原子が有する電荷パラメータを、部分電荷と称することがある。また、例えば、原子集合体ABを構成する原子の電荷パラメータからフラグメントBの点電荷が持つ電荷パラメータを除去することを、フラグメントBの点電荷を原子集合体ABから除去するとも称することがある。同様に、フラグメントBの点電荷が持つ電荷パラメータを、原子集合体ABを構成する原子の電荷パラメータに付加することをフラグメントBの点電荷を原子集合体ABに付加する、又は原子集合体ABを構成する原子の電荷パラメータに付加するとも称することがある。
【0038】
数式(1)の第3項の∫Δν(φ)P(φ)dφは、変化後の原子集合体ABにおける相互作用エネルギーφの各階級の出現確率であるP(φ)と、変化前の原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABにおける相互作用エネルギーφの各階級における相互作用エネルギーεの各階級の出現確率であるP
0'(ε;φ)と、変化後の原子集合体ABにおける相互作用エネルギーεの各階級の出現確率であるP'(ε)とに基づいて、Δν(φ)を計算することなく、計算することができる。その他に、∫Δν(φ)P(φ)dφは、例えば、相互作用エネルギーφの各階級における、相互作用エネルギーεに起因する自由エネルギー変化量Δν(φ)を計算し、計算されたΔν(φ)と、変化後の原子集合体ABにおける相互作用エネルギーφの各階級の出現確率であるP(φ)とを用いて計算することもできる。
【0040】
本発明においては、詳細を後述する通り、フラグメントBを構成する原子を、原子集合体Aを構成する原子に付加することにより、サンプリング効率を飛躍的に向上させることができる。よって、Δν(φ)を計算する方法は、特に限定されず、例えば従来の粒子挿入法やエネルギー表示法を用いることができる。なお、本発明で計算される数式(1)の自由エネルギー差ΔGは、溶媒和自由エネルギー差や結合自由エネルギー差等、任意の自由エネルギー差の計算に適用することができる。
【0041】
本発明において計算されるΔGは、反応式(1):
A+B→AB (1)
[式中、Aは、構造aからなる又は構造aを含んでなる原子集合体を表し、Bは、1個以上の原子からなるフラグメントを表し、ABは、原子集合体Aと該原子集合体Aの構造aに結合したフラグメントBとからなる原子集合体を表す。]
で表される変化に関する自由エネルギー変化量、すなわち、変化前の原子集合体Aの自由エネルギー及びフラグメントBの自由エネルギーの和G
1と、変化後の原子集合体ABの自由エネルギーG
2との差(G
2−G
1)である。なお、変化前では、原子集合体AとフラグメントBとは相互作用していない。
【0042】
原子集合体A及び原子集合体ABは、1個又は2個以上の原子で構成される原子の集合体である。原子集合体A及び原子集合体ABには、他の原子と結合していない1個又は2個以上の原子が含まれてもよい。ここで、原子には、イオンも含まれる。また、仮想原子(実在しない原子)も含まれる。原子集合体A及び原子集合体ABを構成する原子間結合としては、例えば、共有結合、配位結合、水素結合、静電相互作用、疎水性相互作用等が挙げられる。原子集合体A及び原子集合体ABは、1種の原子間結合(例えば、共有結合)で構成されてもよいし、2種以上の原子間結合の組合せ(例えば、共有結合と、1種又は2種以上の他の原子間結合との組合せ)で構成されてもよい。
【0043】
原子集合体A及び原子集合体ABは、1種又は2種以上の分子から構成されてもよい。
分子としては、例えば、リガンド、該リガンドが結合する蛋白質、水分子等の溶媒分子等が挙げられる。原子集合体Aが1種の分子で構成される場合、原子集合体ABも1種の分子で構成され、原子集合体Aが2種以上の分子で構成される場合、原子集合体ABも2種以上の分子で構成される。原子集合体A及び原子集合体ABが2種以上の分子で構成される場合としては、例えば、原子集合体A及び原子集合体ABが、リガンドの他に、該リガンドが結合する蛋白質及び/又は水分子等の溶媒分子を含む場合等が挙げられる。原子集合体A及び原子集合体ABが、リガンドの他に、該リガンドが結合する蛋白質を含む場合、リガンドと蛋白質とは複合体を形成していてもよい。該複合体におけるリガンドと蛋白質との間の結合様式としては、例えば、配位結合、水素結合、静電相互作用、共有結合、疎水性相互作用等が挙げられる。原子集合体A及び原子集合体ABが、リガンド又はリガンドと蛋白質との複合体を含む場合、リガンド又は複合体と水分子等の溶媒分子との間の結合様式としては、例えば、配位結合、水素結合、静電相互作用、共有結合、疎水性相互作用等が挙げられる。
【0044】
原子集合体Aは、構造aからなるか、又は、構造aを含んでなる。原子集合体Aが1個の分子で構成される場合、構造aは、当該分子の構造の全体であってもよいし(この場合、原子集合体Aは構造aからなる)、当該分子の構造の一部であってもよい(この場合、原子集合体Aは構造aを含んでなる)。原子集合体Aが2種以上の分子で構成される場合、構造aは、当該2種以上の分子の構造の全体であってもよいし(この場合、原子集合体Aは構造aからなる)、当該2種以上の分子の構造の一部であってもよい(この場合、原子集合体Aは構造aを含んでなる)。原子集合体Aが、蛋白質とリガンドとの複合体を含む場合、当該複合体の構造の一部としては、例えば、リガンドの構造の全体又は一部、蛋白質の構造の全体又は一部、リガンドの全体又は一部と蛋白質の一部とからなる部分、受容体の全体又は一部とリガンドの一部とからなる部分等が挙げられる。
【0045】
原子集合体Aは、原子集合体ABからフラグメントBが除去された原子集合体に相当する。例えば、原子集合体ABがアニソール(C
6H
5OCH
3)であり、フラグメントBがメトキシ基(−OCH
3)である場合、原子集合体AはC
6H
5となる。また、例えば、原子集合体ABがアニソール(C
6H
5OCH
3)であり、フラグメントBがメトキシ基(−OCH
3)及び点電荷(該メトキシ基が結合するベンゼン環の炭素原子が有する電荷パラメータを持つ仮想原子)である場合、原子集合体Aは、アニソールからメトキシ基及び該点電荷を除去することにより得られる。つまり、原子集合体AはC
6H
5で表される原子集合体であって、原子集合体ABにおいてメトキシ基が結合していたベンゼン環上の炭素原子が有していた電荷パラメータが零である原子集合体である。
【0046】
原子集合体ABは、原子集合体Aと該原子集合体Aの構造aに結合したフラグメントBとからなる。原子集合体ABは、原子集合体Aに対してフラグメントBが結合した原子集合体に相当する。例えば、原子集合体Aがアニソールからメトキシ基及び点電荷(該メトキシ基が結合するベンゼン環の炭素原子が有する部分電荷)を除去した構造であって、フラグメントBがメトキシ基及び該点電荷より構成される場合、原子集合体ABはアニソールとなる。
【0047】
構造aを構成する原子の数は、1個以上である限り特に限定されないが、好ましくは3個以上、さらに好ましくは5個以上である。構造aを構成する原子の数が3個以上である場合、原子集合体Aの構造aへのフラグメントBの結合を効率よく実施することができるため、サンプリング効率を向上させることができる。また、構造aを構成する原子の数の上限値は特に限定されない。
【0048】
構造aを構成する原子の種類は特に限定されない。構造aを構成する原子としては、例えば、炭素原子、水素原子、窒素原子、酸素原子等が挙げられる。構造aは、1種の原子で構成されてもよいし、2種以上の原子で構成されてもよい。
【0049】
構造aを構成する原子間の結合様式は特に限定されない。構造aを構成する原子間の結合様式としては、例えば、共有結合、配位結合、水素結合、静電相互作用、疎水性相互作用等が挙げられる。構造aは、1種の原子間結合(例えば、共有結合)で構成されてもよいし、2種以上の結合の組合せ(例えば、共有結合と、1種又は2種以上の他の結合との組合せ)で構成されてもよい。
【0050】
構造aを構成する原子には、1個又は2個以上の仮想原子(実在しない原子)が含まれていてもよい。構造aは、実在原子のみで構成されていてもよいし、仮想原子のみで構成されてもよいし、実在原子と仮想原子との組合せで構成されてもよい。仮想原子としては、例えば、原子価殻が形式的に8電子有していない原子が挙げられる。なお、周期表の1〜2属及び13〜18属に属する典型元素に該当する原子は、通常、その原子価殻に形式的に8電子を有した電子配置をとってオクテット則を満足して、化学的に安定に存在することが知られている。その他の仮想原子としては、例えば、他の原子との間に、結合長、結合角、捻れ角等が関係する結合に関する相互作用やファンデールワールス相互作用を有しない(つまり、ファンデールワールスポテンシャルがゼロ)が、他の原子間との静電相互作用を有する(つまり、電荷パラメータを持ち、クーロンポテンシャルを有する)点電荷、他の原子間との静電相互作用はなく(つまり、クーロンポテンシャルがゼロ)、他の原子との間に結合長、結合角、捻れ角等に起因する結合に関する相互作用及びファンデールワールス相互作用を有する原子、イオン状態で存在する原子、他の原子間に結合長、結合角、捻れ角等に起因する結合に関する相互作用、ファンデールワールス相互作用及び静電相互作用のいずれも有さないダミー原子等が挙げられる。なお、以後、ファンデールワールス相互作用に起因するエネルギーをファンデールワールスポテンシャルと称し、静電相互作用に起因するエネルギーをクーロンポテンシャルと称する。
【0051】
フラグメントBを構成する原子の数は、1個以上である限り特に限定されない。フラグメントBを構成する原子の数の上限は特に限定されないが、通常25個、好ましくは10個である。
【0052】
フラグメントBを構成する原子の種類は特に限定されない。フラグメントBを構成する原子としては、例えば、炭素原子、水素原子、窒素原子、酸素原子等が挙げられる。フラグメントBは、1種の原子で構成されてもよいし、2種以上の原子で構成されてもよい。
【0053】
フラグメントBを構成する原子間の結合様式は特に限定されない。フラグメントBを構成する原子間の結合様式としては、例えば、共有結合、配位結合、水素結合、静電相互作用、疎水性相互作用等が挙げられる。フラグメントBは、1種の原子間結合(例えば、共有結合)で構成されてもよいし、2種以上の結合の組合せ(例えば、共有結合と、1種又は2種以上の他の結合との組合せ)で構成されてもよい。
【0054】
フラグメントBを構成する原子には、1個又は2個以上の仮想原子(実在しない原子)が含まれていてもよい。フラグメントBは、実在原子のみで構成されてもよいし、仮想原子のみで構成されてもよいし、実在原子と仮想原子との組合せで構成されてもよい。仮想原子としては、構造aと同様の具体例が挙げられる。フラグメントBが、点電荷を仮想原子として含み、実在原子と仮想原子との組合せで構成される場合、該原子間には共有結合はなく、結合長、結合角、捻れ角等に起因する結合に関する相互作用はない。また、コンピュータシミュレーション上、フラグメントBの仮想原子(点電荷)と実在原子との間の静電相互作用は計算する方が好ましいが、自由エネルギーを計算する過程において統一されていれば、計算しても計算しなくてもどちらでも構わない。フラグメントBの仮想原子(点電荷)と実在原子との間に共有結合があると仮定し、該原子間に3個以上の共有結合を有する場合に、分子動力学シミュレーション等のコンピュータシミュレーションにおける、所謂「1−4相互作用」を考慮し、仮想原子(点電荷)と実在原子との間の静電相互作用を計算するのが好ましい。例えば、構造aBがアニソール(C
6H
5OCH
3)であり、フラグメントBがメトキシ基(−OCH
3)及び点電荷(該メトキシ基が結合するベンゼン環の炭素原子が有する電荷パラメータを有する仮想原子)である場合、該点電荷とメトキシ基の水素原子それぞれとの間の静電相互作用を計算するのが好ましい。
【0055】
以下、本発明の一実施形態に係る計算装置について、図面に基づいて説明する。
図1に示すように、計算装置1は、入力部11と、制御部12と、記憶部13と、出力部14とを備えており、これらはシステムバスで相互に接続されている。計算装置1は、汎用のコンピュータを基本ハードウェアとして用いることにより実現可能である。
【0056】
入力部11は、例えば、オペレータが操作するキーボード、マウス等のポインティングデバイスから構成され、オペレータからの指示(例えば、プログラムの実行の指示、処理結果の表示指示)、各処理に必要なデータの入力等の各種操作信号を入力する。入力されたデータは、制御部12により記憶部13に記憶される。
【0057】
制御部12は、例えば、CPU、RAM、ROM等から構成され、記憶部13に記憶されている各種データ、各種プログラム等に基づいて、ΔGを計算する。この際、制御部12は、ΔGを計算するための計算プログラムを実行することにより、第1の原子集合体モデル作成手段12A、第1の座標取得手段12B、第2の座標取得手段12C、第1の相互作用エネルギーφ度数分布作成手段12D、第1の相互作用エネルギーφ出現確率計算手段12E、第1の相互作用エネルギーε度数分布作成手段12F、第1の相互作用エネルギーε出現確率計算手段12G、第2の原子集合体モデル作成手段12H、第3の座標取得手段12I、第2の相互作用エネルギーφ度数分布作成手段12J、第2の相互作用エネルギーφ出現確率計算手段12K、第2の相互作用エネルギーε度数分布作成手段12L、第2の相互作用エネルギーε出現確率計算手段12M、∫Δν(φ)P(φ)dφ計算手段12N、及びΔG計算手段12Oとして機能する。なお、制御部12のCPUは、複数の演算コアを有し、各演算コアが、ΔGを計算するための計算プログラムの実行によって、各種手段として機能する構成であってもよい。こうした構成によれば、複数の計算を並列して処理することができ、これにより、ΔGの計算に要する時間を短縮することができる。
【0058】
記憶部13は、例えば、RAM、ハードディスク等のストレージから構成され、各種データ、各種プログラム等を記憶する。記憶部13に記憶されるデータ及びプログラムとしては、例えば、変化前の原子集合体Aに関するデータ13A、第1のシミュレーション条件に関するデータ13B、第1のシミュレーションプログラム13C、変化後の原子集合体ABに関するデータ13D、第2のシミュレーション条件に関するデータ13E、第2のシミュレーションプログラム13F等が挙げられる。なお、本実施形態では、第1のシミュレーションプログラム13C及び第2のシミュレーションプログラム13Fを別個のプログラムとして記載しているが、第1のシミュレーションプログラム13C及び第2のシミュレーションプログラム13Fは1個のプログラムであってもよい。
【0059】
出力部14は、例えば、ディスプレイ等から構成され、各種手段によって取得又は計算された結果(例えば、ΔG計算手段12Oによって計算されたΔG等)を出力する。
【0060】
第1の原子集合体モデル作成手段12Aは、変化前の原子集合体Aをモデル化した第1の原子集合体モデルを作成する。第1の原子集合体モデル作成手段12Aによって作成された第1の原子集合体モデルに関するデータ(例えば、原子集合体Aを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)は、制御部12により記憶部13に記憶される。
【0061】
第1の原子集合体モデル作成手段12Aは、第1の原子集合体モデルを作成する際、記憶部13に記憶されている変化前の原子集合体Aに関するデータ13Aを用いる。変化前の原子集合体Aに関するデータ13Aは、例えば、制御部12が読み取り可能なファイルの形態で記憶部13に記憶されている。変化前の原子集合体Aに関するデータ13Aは、変化前の原子集合体Aをモデル化した第1の原子集合体モデルを作成可能である限り特に限定されない。変化前の原子集合体Aに関するデータ13Aとしては、例えば、原子集合体Aを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等が挙げられる。
【0062】
第1の原子集合体モデルは、コンピュータシミュレーションが可能である限り特に限定されない。原子集合体モデルとしては、全原子モデル、ビーズスプリングモデル、ユナイテッドアトムモデル等が挙げられる。なお、ビーズスプリングモデルは、原子集合体Aを構成するモノマーユニットを1つのビーズ(セグメント)とみなし、仮想的なバネで結合させたモデルであり、ユナイテッドアトムモデルは、水素原子を重原子(例えば、炭素原子)に含めて1つの原子(質点)として取り扱うモデルである。
【0063】
第1の座標取得手段12Bは、第1の原子集合体モデル作成手段12Aによって作成された第1の原子集合体モデルに対するコンピュータシミュレーションの結果として出力されるスナップショットより、第1〜第iの状態F
1〜F
i(iは2以上の整数である。)の各状態における原子集合体Aの座標を取得する。ここで、スナップショットには、原子集合体Aを構成する全ての原子の座標が含まれる。つまり、状態F
1〜F
iの各状態におけるスナップショットには、それぞれ、状態F
1〜F
iの各状態における原子集合体Aを構成する全ての原子の座標が含まれる。以下、状態F
1、F
2、・・・、F
iにおける原子集合体Aの座標を、それぞれ、座標R
A(F
1)、R
A(F
2)、・・・、R
A(F
i)という場合がある。第1の座標取得手段12Bによって取得された原子集合体Aを構成する全ての原子の座標は、制御部12により記憶部13に記憶される。
【0064】
第1の座標取得手段12Bが実行するコンピュータシミュレーションは、統計力学理論に基づく限り特に限定されない。代表的なコンピュータシミュレーションとしては、例えば、分子動力学法、モンテカルロ法、分子動力学法又はモンテカルロ法と第一原理計算とを組み合わせた方法等が挙げられる。
【0065】
第1の座標取得手段12Bは、記憶部13に記憶されている第1の原子集合体モデルに関するデータ、第1のシミュレーション条件に関するデータ13B、第1のシミュレーションプログラム13C等に基づいてコンピュータシミュレーションを実行する。第1のシミュレーション条件に関するデータ13Bは、制御部12が読み取り可能なファイルの形態で記憶部13に記憶されている。
【0066】
シミュレーション条件は、コンピュータシミュレーションを実行可能である限り特に限定されない。シミュレーション条件としては、例えば、温度条件、圧力条件、ポテンシャルパラメータの種類及び値、発生させるアンサンブルの種類、境界条件、スナップショットの数等の出力条件等が挙げられる。
【0067】
ポテンシャルパラメータとしては、例えば、結合長、結合角、捻れ角等が関係する原子間の結合に関するパラメータ、及び原子間に働くファンデールワールス相互作用及び静電相互作用等に関するパラメータが挙げられる。
【0068】
ポテンシャルパラメータとしては、例えば、Amber、GAFF、CHARMm(登録商標)、CGenFF(CHARMm36)、DISCOVER、GROMOS、DREIDING、OPLS等の公知のポテンシャルパラメータを用いることができる。
【0069】
例えば、ポテンシャルパラメータとしてAmberを用いる場合、原子集合体Aを構成する各原子に対して割与えられる原子タイプ及び各原子間の原子タイプのペアに基づいて、結合長、結合角、ファンデールワールス等に関するパラメータの値が決定される。
【0070】
ポテンシャルパラメータは、原子間の結合に関するポテンシャルパラメータと結合の有無に関わらない非結合ポテンシャルパラメータに分けられる。計算精度を考慮した場合、以下に述べる方法によりポテンシャルパラメータを決定することが好ましい。
【0071】
原子間の結合に関するポテンシャルパラメータは、第一原理計算に基づく量子化学計算によって決定することが好ましく、量子化学計算としては、ハートレーフォック法(以下「HF法」という。)を用いることが好ましく、汎関数としてB3LYPを用いた密度汎関数法(以下「B3LYP法」という。)を用いることがさらに好ましい。HF法又はB3LYP法を用いる場合、1電子軌道を展開する基底関数を指定する必要がある。計算時間と計算精度を勘案すると、6−31G(d)基底関数、6−31G(d、p)基底関数、6−31+G(d、p)基底関数を好適に用いることができる。
【0072】
非結合ポテンシャルパラメータは、更に、ファンデールワールスポテンシャル及びクーロンポテンシャルに関するパラメータに分けられる。前者は、密度と蒸発熱の実験値を再現するように決定することが好ましい。後者のパラメータは、原子が有する部分電荷であり、量子化学計算から求めた静電ポテンシャル(ESP)を再現するように決定することが好ましい。ESPを再現する方法としては、公知のMK法やCHelpG法を好適に用いることができる。
【0073】
以上の量子化学計算は、有償又は無償で入手可能な種々の量子化学計算プログラムパッケージを利用することができる。例えば、「Gaussian98(登録商標)」、「Gaussian03(登録商標)」、「Gaussian09(登録商標)」、「GAMESS」といった製品名で市販又は公開されている汎用プログラムを好適に用いることができる。
【0074】
コンピュータシミュレーションとして分子動力学法と第一原理計算とを組み合わせた方法を用いる場合、原子集合体Aを構成する原子の全てに対して第一原理計算を実施してもよいし、原子集合体Aを構成する原子の一部に対して第一原理計算を実施し、かつ、原子集合体Aを構成する残りの原子に対しては分子力場法を実施してもよい。前者では、各原子に関する結合情報、ポテンシャルパラメータ、原子タイプは不要であり、後者では、第一原理計算を適用する原子範囲に従って、各原子に関する結合情報、ポテンシャルパラメータ、原子タイプ等を適宜調節すればよい。
【0075】
コンピュータシミュレーションによって発生させるアンサンブルとしては、例えば、NVEアンサンブル(粒子数N、体積V、エネルギーEが一定)、NVTアンサンブル(粒子数N、体積V、温度Tが一定)、NPTアンサンブル(粒子数N、圧力P、体積Vが一定)等が挙げられる。なお、アンサンブルは、系がとり得る微視的な状態の集団(統計集団)である。
【0076】
境界条件としては、例えば、周期境界条件等が挙げられる。分子動力学シミュレーションにおいて周期境界条件が採用される場合、基本セルの周囲にイメージセルが配置される。
【0077】
温度条件は、反応式(1)で表される変化が生じる際の温度に設定される。なお、温度は絶対温度であり、単位はケルビン(K)である。
【0078】
第1の座標取得手段12Bは、コンピュータシミュレーションによって発生するアンサンブルに含まれる状態F
1〜F
iの各状態における原子集合体Aの座標を取得する。状態F
1〜F
iの各状態は、コンピュータシミュレーションの実行により発生する微視的な状態である。状態F
1〜F
iは、コンピュータシミュレーションによって発生せしめた微視的状態のうちの一部であってもよいし、全部であってもよい。
【0079】
iの値は、2以上である限り特に限定されず、計算対象である自由エネルギーの種類(例えば、溶媒和自由エネルギー、結合自由エネルギー等)に応じて適宜選択することができる。計算対象の自由エネルギーが溶媒和自由エネルギーである場合、iの値は、好ましくは10000以上、より好ましくは100000以上、さらに好ましくは1000000以上である。計算対象の自由エネルギーが結合自由エネルギーである場合、iの値は、好ましくは100000以上、さらに好ましくは1000000以上である。なお、iの上限値は特に限定されないが、通常10000000、好ましくは5000000である。
【0080】
コンピュータシミュレーションとして分子動力学シミュレーションを用いる場合、原子集合体Aの状態は、初期状態(時間T
0)から、状態F
1(時間T
1)、状態F
2(時間T
2)、・・・、状態F
i(時間T
i)と経時的に変化する。そして、時系列に沿って、状態F
1〜F
iの各状態における原子集合体Aの座標が取得される。取得された状態F
1〜F
i(時間T
1〜T
i)の各状態(各時間)における原子集合体Aの座標は、それが取得された各状態(各時間)と対応付けられて、制御部12により記憶部13に記憶される。
【0081】
コンピュータシミュレーションとしてモンテカルロ法シミュレーションを用いる場合、乱数を発生させることにより状態F
1〜F
iを作り出し、状態F
1〜F
iの各状態における原子集合体Aの座標が取得される。取得された状態F
1〜F
iの各状態における原子集合体Aの座標は、それが取得された各状態と対応付けられて、制御部12により記憶部13に記憶される。
【0082】
コンピュータシミュレーションとして分子動力学シミュレーションを用いる場合、第1の座標取得手段12Bは、第1の原子集合体モデルに対して、分子動力学シミュレーションの実行前にエネルギー極小化計算を実行することが好ましい。エネルギー極小化計算は、例えば、記憶部13に記憶されている第1のシミュレーションプログラム13Cの機能に基づいて実行することができる。エネルギー極小化計算によって、原子集合体モデルの初期構造が有する不自然な構造の歪みが除去され、コンピュータシミュレーションの初期段階における時間積分の発散を回避することができる。
【0083】
コンピュータシミュレーションとして分子動力学シミュレーションを用いる場合、第1の座標取得手段12Bは、第1の原子集合体モデル(好ましくはエネルギー極小化計算後の原子集合体モデル)を平衡化し、平衡化後の状態F
1〜F
iにおける原子集合体Aの座標を取得することが好ましい。例えば、第1の座標取得手段12Bは、第1の原子集合体モデル(好ましくはエネルギー極小化計算後の原子集合体モデル)に対してコンピュータシミュレーションを実行し、例えば、ある物理量の値の変動幅が一定に達した段階で、又は、ある一定時間が経過した段階で、原子集合体モデルが平衡化したと判断し、平衡化後の状態F
1〜F
iにおける原子集合体Aの座標を取得する。
【0084】
第2の座標取得手段12Cは、第1の座標取得手段12Bによって取得された原子集合体Aの座標に基づいて、状態F
1〜F
iの各状態における原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABの座標を取得する。
【0085】
第2の座標取得手段12Cによって取得される原子集合体ABの座標は、原子集合体ABを構成する全ての原子の座標を含む。以下、状態F
x(x=1,2,・・・,i)における原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABの座標を「R
AB(F
x)」と表記する場合がある。取得された原子集合体ABの座標は、その基礎となった原子集合体Aの状態と対応付けられて、制御部12により記憶部13に記憶される。すなわち、原子集合体ABの座標R
AB(F
1)、R
AB(F
2)、・・・、R
AB(F
i)は、それぞれ、状態F
1、F
2、・・・、F
iと対応付けられて、制御部12により記憶部13に記憶される。
【0086】
状態F
x(x=1,2,・・・,i)における原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABは、1種類の原子集合体であってもよいし、2種類以上の原子集合体であってもよい。すなわち、座標R
AB(F
x)は、状態F
xにおける原子集合体Aに対して、構造aに対する相対位置が異なる第1、第2、・・・、第pのフラグメントBが結合することにより生じる原子集合体ABの座標を意味する場合がある。pは、xの値から独立して選択された、2以上の整数である。したがって、pはx=1,2,・・・,iにおいて同一の値であってもよいし、x=1,2,・・・,iにおいて異なる値であってもよい。以下、第1、第2、・・・、第pのフラグメントBを、それぞれ、「フラグメントB
1」、「フラグメントB
2」、・・・、「フラグメントB
p」と表記し、状態F
xにおける原子集合体Aに対して、フラグメントB
1、B
2、・・・、B
pが結合することにより生じる原子集合体ABの座標を、それぞれ、「R
AB(F
x,B
1)」、「R
AB(F
x,B
2)」、・・・、「R
AB(F
x,B
p)」と表記する場合がある。
【0087】
状態F
xにおける原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABの座標は、実際に原子集合体ABを発生させて取得してもよいし、実際に原子集合体ABを発生させることなく取得してもよい。後者の場合、例えば、状態F
xにおける原子集合体Aに対して結合させるべきフラグメントBの構造aに対する相対位置を決定し、状態F
xにおける原子集合体Aの座標と、決定されたフラグメントBの構造aに対する相対位置とに基づいて、状態F
xにおける原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABの座標R
AB(F
x)を取得することができる。
【0088】
状態F
xにおける原子集合体Aに対して結合させるべきフラグメントBの構造aに対する相対位置は特に限定されるものではないが、原子集合体Aを構成する原子と、フラグメントBを構成する原子とが衝突する構造を有する原子集合体ABが発生する頻度を低減する観点から、構造aを構成する原子の座標に対してフラグメントBを構成する原子の座標を限定する拘束条件を課すことが好ましい。
【0089】
拘束条件は、構造aを構成する原子の座標に対してフラグメントBを構成する原子の座標を限定し得る限り特に限定されない。例えば、フラグメントBを構成する一部の原子について、その座標を、構造aを構成する原子の座標に対して限定する条件であってもよいし、フラグメントBを構成する全ての原子について、その座標を、構造aを構成する原子の座標に対して限定する条件であってもよい。
【0090】
一実施形態において、第2の座標取得手段12Cは、以下の処理:
状態F
xにおける原子集合体Aに結合させるべきフラグメントBに関し、フラグメントBを構成する原子から選択された原子bの、構造aに対する相対位置条件(例えば、構造aを構成する1個又は2個以上の原子に対する距離、角度等)を規定する処理W101;
規定された相対位置条件に基づいて、原子bの座標を取得する処理W102;
取得された原子bの座標と、状態F
xにおける原子集合体Aの座標R
A(F
x)とに基づいて、状態F
xにおける原子集合体Aに対して原子bが結合することにより生じる原子集合体Abの座標R
Ab(F
x)を取得する処理W103;並びに
原子bの座標を基準として、コンピュータシミュレーション、公知の配座発生方法等により、フラグメントBを構成する原子b以外の原子の座標を発生させ、状態F
xにおける原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABの座標R
AB(F
x)を取得する処理W104
を、x=1〜iで繰り返して実行する。これにより、状態F
1〜F
iの各状態における原子集合体Aに対して、構造aに対する相対位置が限定されたフラグメントBが結合することによって生じる原子集合体ABの座標R
AB(F
1)〜R
AB(F
i)を取得することができる。
【0091】
第2の座標取得手段12Cは、処理W101において、状態F
xにおける原子集合体Aに結合させるべきフラグメントBに関し、互いに異なる第1、第2、・・・、第pの相対位置条件を規定してもよい。pは、xの値から独立して選択された、2以上の整数である。したがって、pはx=1,2,・・・,iにおいて同一の値であってもよいし、x=1,2,・・・,iにおいて異なる値であってもよい。
【0092】
処理W101において、状態F
xにおける原子集合体Aに結合させるべきフラグメントBに関し、互いに異なる第1、第2、・・・、第pの相対位置条件が規定される場合、第2の座標取得手段12Cは、第1〜第pの相対位置条件のそれぞれに関し、処理W102〜W104を繰り返す。これにより、状態F
xにおける原子集合体Aに対して、構造aに対する相対位置が異なるフラグメントB
1、B
2、・・・、B
pが結合することにより生じる原子集合体ABの座標R
AB(F
x,B
1)、R
AB(F
x,B
2)、・・・、R
AB(F
x,B
p)を取得することができる。
【0093】
別の実施形態において、第2の座標取得手段12Cは、以下の処理:
構造aと該構造aに結合したフラグメントBとからなる又は構造aと該構造aに結合したフラグメントBとを含んでなる原子集合体Cをモデル化した第3の原子集合体モデルを作成する処理W201;並びに
作成された第3の原子集合体モデルに対するコンピュータシミュレーションにより、第1〜第kの状態H
1〜H
k(kは2以上の整数である。)の各状態における原子集合体Cの座標(以下、状態H
1、H
2、・・・、H
kにおける原子集合体Cの座標を、それぞれ、「R
C(H
1)」、「R
C(H
2)」、・・・、「R
C(H
k)」と表記する場合がある。)を取得する処理W202
を実行した後、
構造aを構成する原子から選択された1個又は2個以上の原子からなる選択原子群(なお、2個以上の原子が選択される場合だけでなく、1個の原子が選択される場合も含めて便宜上「選択原子群」という。)を選択して、状態F
xにおける原子集合体Aの選択原子群の座標に対し、状態H
1〜H
kから選択された1又は2以上の状態における原子集合体Cの選択原子群の座標を回転及び/又は並進させることにより、原子集合体Aの選択原子群と原子集合体Cの選択原子群との間で、対応する原子間の距離の二乗和が最小となる原子集合体Cの座標を作成(原子集合体Aの選択原子群の座標に原子集合体Cの選択原子群の座標をフィッティング)し、作成された原子集合体Cの座標に基づいて、原子集合体Aに状態H
1〜H
kから選択された1又は2以上の状態における原子集合体Cを重ね合わせる処理W203;並びに
状態F
xにおける原子集合体Aの座標R
A(F
x)と、状態F
xにおける原子集合体Aに重ね合わせた原子集合体CのフラグメントBの1又は2以上の座標とに基づいて、状態F
xにおける原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABの座標を取得する処理W204
を、x=1〜iで繰り返して実行する。これにより、状態F
1〜F
iの各状態における原子集合体Aに対して、構造aに対する相対位置が限定されたフラグメントBが結合することによって生じる原子集合体ABの座標R
AB(F
1)〜R
AB(F
i)を取得することができる。
【0094】
原子集合体Cは、構造aと該構造aに結合したフラグメントBとからなる分子C’のみから構成されてもよいし(この場合、分子C’の周囲は真空である)、構造aと該構造aに結合したフラグメントBとからなる分子C’と1種又は2種以上の他の分子(例えば、水分子等の溶媒分子)とから構成されてもよい。
【0095】
処理W201において、第2の座標取得手段12Cは、原子集合体Cをモデル化した第3の原子集合体モデルを作成する。第2の座標取得手段12Cによって作成された第3の原子集合体モデルに関するデータ(例えば、原子集合体Cを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)は、制御部12により記憶部13に記憶される。
【0096】
第2の座標取得手段12Cは、第3の原子集合体モデルを作成する際、記憶部13に記憶されている原子集合体Cに関するデータ(不図示)を用いる。原子集合体Cに関するデータは、例えば、制御部12が読み取り可能なファイルの形態で記憶部13に記憶されている。原子集合体Cに関するデータは、原子集合体Cをモデル化した第3の原子集合体モデルを作成可能である限り特に限定されない。原子集合体Cに関するデータとしては、例えば、分子C’に関するデータ(例えば、分子C’を構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)、分子C’の周辺環境に関するデータ(例えば、分子C’の周囲に位置する溶媒分子の有無、種類、個数、溶媒分子を構成する原子の座標等)等が挙げられる。
【0097】
第3の原子集合体モデルは、コンピュータシミュレーションが可能である限り特に限定されない。原子集合体モデルとしては、全原子モデル、ビーズスプリングモデル、ユナイテッドアトムモデル等が挙げられる。
【0098】
第3の原子集合体モデルに対するコンピュータシミュレーションは、統計力学理論に基づく限り特に限定されない。代表的なコンピュータシミュレーションとしては、例えば、分子動力学法、モンテカルロ法、分子動力学法又はモンテカルロ法と第一原理計算とを組み合わせた方法等が挙げられる。第3の原子集合体モデルに対するコンピュータシミュレーションは、第1の原子集合体モデルに対するコンピュータシミュレーションと同様に実行することができる。
【0099】
コンピュータシミュレーションとして分子動力学シミュレーションを用いる場合、第2の座標取得手段12Cは、第3の原子集合体モデルに対して、分子動力学シミュレーションの実行前にエネルギー極小化計算を実行することが好ましい。エネルギー極小化計算は、例えば、記憶部13に記憶されている第1のシミュレーションプログラム13Cの機能に基づいて実行することができる。エネルギー極小化計算によって、原子集合体モデルの初期構造が有する不自然な構造の歪みが除去され、コンピュータシミュレーションの初期段階における時間積分の発散を回避することができる。
【0100】
コンピュータシミュレーションとして分子動力学法を用いる場合、第2の座標取得手段12Cは、第3の原子集合体モデル(好ましくはエネルギー極小化計算後の原子集合体モデル)を平衡化し、平衡化後の状態H
1〜H
kにおける原子集合体Cの座標を取得することが好ましい。例えば、第2の座標取得手段12Cは、第3の原子集合体モデル(好ましくはエネルギー極小化計算後の原子集合体モデル)に対してコンピュータシミュレーションを実行し、例えば、ある物理量の値の変動幅が一定に達した段階で、又は、ある一定時間が経過した段階で、原子集合体モデルが平衡化したと判断し、平衡化後の状態H
1〜H
kにおける原子集合体Cの座標を取得する。
【0101】
第2の座標取得手段12Cは、処理W203において、状態F
xにおける原子集合体Aの座標R
A(F
x)に対して、状態H
1〜H
kにおける原子集合体Cの座標から、構造aに対するフラグメントBの相対位置が互いに異なるp組の座標を選択してもよい。pは、xの値から独立して選択された、2以上の整数である。したがって、pはx=1,2,・・・,iにおいて同一の値であってもよいし、x=1,2,・・・,iにおいて異なる値であってもよい。
【0102】
処理W203において、状態H
1〜H
kにおける原子集合体Cの座標から、構造aに対するフラグメントBの相対位置が互いに異なるp組の座標が選択される場合、第2の座標取得手段12Cは、p組の座標のそれぞれに関し、処理W203〜W204を繰り返す。これにより、状態F
xにおける原子集合体Aに対して、構造aに対する相対位置が異なるフラグメントB
1、B
2、・・・、B
pが結合することにより生じる原子集合体ABの座標R
AB(F
x,B
1)、R
AB(F
x,B
2)、・・・、R
AB(F
x,B
p)を取得することができる。
【0103】
処理W203において、選択原子群の座標を基準として、状態F
xにおける原子集合体Aの座標R
A(F
x)に対して、状態H
1〜H
kから選択された1又は2以上の状態における原子集合体Cの座標を回転及び/又は並進させることにより、原子集合体Aに状態H
1〜H
kから選択された1又は2以上の状態における原子集合体Cを重ね合わせる方法は、最小二乗法が好ましい。該最小二乗法を用いることによって、状態H
1〜H
kから選択された1又は2以上の状態における原子集合体Cの座標のうち、選択原子群を構成する原子と、フラグメントBを構成する原子との相対的な位置関係を維持しつつ、状態H
1〜H
kから選択された1又は2以上の状態における原子集合体Cの座標のうち、フラグメントBを構成する原子の座標を、状態F
xにおける原子集合体Aの座標R
A(F
x)に対して付加することができる。これにより、フラグメントBを構成する原子と構造aを構成する原子とが衝突する構造の発生頻度を下げて、原子集合体Aに対してフラグメントBを結合させることができる。つまり、サンプリング効率を向上させることができ、統計的に高い精度で計算結果を得ることができる。
【0104】
最小二乗法の具体的な実施方法は特に限定されない。例えば、最小二乗法は、以下のように実施することができる。状態F
1〜F
iの各状態における原子集合体Aの座標のうち、選択原子群の座標を「(x
s_A:n、y
s_A:n、z
s_A:n)」と表記し、状態H
1〜H
kから選択された1つの状態における原子集合体Cの座標のうち、選択原子群の座標を「(x
C:n、y
C:n、z
C:n)」と表記する場合、原子集合体Cの座標のうち、選択原子群の原子間の相対座標(内部座標)を保持する条件の下、原子集合体Cの選択原子群の座標を回転及び/又は並進させることにより、数式(7)が最小となるように、(x
T:n、y
T:n、z
T:n)を求める。ここで、座標の添え字で与えられているnは、選択原子群を構成する個々の原子に対して与えられるものであり、座標(x
s_A:n、y
s_A:n、z
s_A:n)及び(x
T:n、y
T:n、z
T:n)は選択原子群を構成する個々の原子の座標を表す。さらに、数式(7)のCnは、原子集合体A及び原子集合体Cにおける選択原子群を構成する各原子に対して設定可能であって、零以上の実数である。選択原子群を構成する各原子に対して設定するCnの相対比によって、数式(7)に基づいた最小二乗法による選択原子群を構成する原子の座標の一致の重要度を設定することができる。例えば、選択原子群を構成する原子が重原子(水素原子以外の原子)と水素原子の両方を含む場合、水素原子に比して該重原子のCnを大きく設定することにより、原子集合体A及び原子集合体Cにおける選択原子群の該重原子の座標を重点的に一致させることができる。
【数9】
【0105】
次に、最小二乗法により得られた座標(x
T:n、y
T:n、z
T:n)に、原子集合体Cにおける選択原子群を構成する原子の座標が一致するように、原子集合体Cを構成する原子を回転及び/又は並進させることによって、原子集合体Cを構成する原子間の相対座標(内部座標)を保持しつつ、原子集合体Cを構成する原子の座標を移動せしめ、該原子集合体CにおけるフラグメントBを構成する原子の座標(ただし、点電荷等の仮想原子は除く)を原子集合体ABにおけるフラグメントBを構成する原子の座標とすることにより、原子集合体ABの座標R
AB(F
1)〜R
AB(F
i)を取得することができる。
【0106】
フラグメントBが点電荷(仮想原子)を含む場合、第2の座標取得手段12Cは、フラグメントBの点電荷を、原子集合体Aの構造aを構成する原子の電荷パラメータに付加することが好ましい。コンピュータシミュレーション上、仮想原子を除いたフラグメントBを構成する原子を、原子集合体Aの構造aを構成する原子へ結合させる際、原子集合体Aの構造aを構成する原子の部分電荷(電荷パラメータ)を変える必要が生じる場合がある。例えば、仮想原子を除いたフラグメントBを構成する原子が、電子吸引基又は電子供与基を構成し、フラグメントBを構成する原子を、原子集合体Aの構造aを構成する原子へ結合させる際、原子集合体Aの構造aを構成する原子と、仮想原子を除いたフラグメントBを構成する原子との間に結合が生じる場合、通常、原子集合体Aの構造aを構成する原子のうち、フラグメントBと結合する原子及び/又はフラグメントBと結合する原子の周囲の原子(例えば、原子集合体Aを構成する原子のうち、仮想原子を除いたフラグメントBを構成する原子と結合する原子と結合している原子)の部分電荷(電荷パラメータ)は変化する。フラグメントBとの結合の有無により部分電荷(電荷パラメータ)が変化する構造aを構成する原子を、フラグメントBを構成する原子(実在原子)として取り扱い、構造aを構成する原子として取り扱わない方法も考えられるが、該原子の部分電荷(電荷パラメータ)の変化量と同等の電荷を有した点電荷を、フラグメントBを構成する原子(仮想原子)として含めることにより、原子集合体AB(終状態)と原子集合体A(始状態)との変化を小さくすることができるため、計算精度を向上させることが可能となる。
【0107】
第1の相互作用エネルギーφ度数分布作成手段12Dは、第2の座標取得手段12Cによって取得された原子集合体ABの座標に基づいて、構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφを計算し、相互作用エネルギーφの各階級の度数を表す度数分布を作成する。計算された相互作用エネルギーφ及び作成された度数分布は、制御部12により記憶部13に記憶される。相互作用エネルギーφの各階級は、その計算の基礎とされた原子集合体ABの座標と対応付けられて、制御部12により記憶部13に記憶される。
【0108】
相互作用エネルギーφは、原子集合体ABの座標R
AB(F
1)〜R
AB(F
i)のそれぞれに関して計算される。相互作用エネルギーφは、例えば、記憶部13に記憶されている第1のシミュレーションプログラム13Cの機能に基づいて計算することができる。相互作用エネルギーφを計算する際、用いられるシミュレーションプログラム及びポテンシャルパラメータは特に限定されるものではないが、原子集合体ABの座標は、第1の座標取得手段12Bが実行するコンピュータシミュレーションにより求められるものであるので、第1の座標取得手段12Bが実行するコンピュータシミュレーションで用いられるものと同一のシミュレーションプログラム及びポテンシャルパラメータを用いることが好ましい。
【0109】
相互作用エネルギーφの各階級の度数を表す度数分布は、例えば、横軸が相互作用エネルギーφの各階級を表し、縦軸が相互作用エネルギーφの各階級の度数を表すヒストグラムとして作成することができる。
【0110】
相互作用エネルギーφの各階級は、例えば、相互作用エネルギーφを中心とし、相互作用エネルギー間隔Δφを持った相互作用エネルギー区間[φ−Δφ/2〜φ+Δφ/2]を階級間隔として作成することができ、各相互作用エネルギー区間におけるΔφは、すべての区間について一定のΔφであってもよいし、相互作用エネルギー区間によって適宜変更したΔφを用いても構わない。数式(1)の第2項の計算を行う上で、P
0(φ)の階級間隔ΔφとP(φ)の階級間隔Δφとを揃えることが好ましく、P
0(φ)の階級間隔Δφに対してP(φ)の階級間隔Δφを揃えることが好ましい。相互作用エネルギーφの分割数、つまり階級間隔Δφの個数は、特に限定されないが、好ましくは50〜500、さらに好ましくは250〜500である。
【0111】
数式(1)の第2項には、相互作用エネルギーφの出現確率P(φ)を相互作用エネルギーφの出現確率P
0(φ)によって除する計算過程を含むため、P(φ)が零とは異なる値であって、P
0(φ)が零となる相互作用エネルギー区間が存在すると、数式(1)の第2項は発散してしまい、数式(1)に基づいた自由エネルギー計算は出来なくなってしまう。そこで、すべての相互作用エネルギー区間[φ−Δφ/2〜φ+Δφ/2]において、P
0(φ)が一定となるように相互作用エネルギー区間毎にΔφを適宜設定することが好ましい。相互作用エネルギー区間毎にΔφを適宜設定し、P
0(φ)の階級間隔Δφに対してP(φ)の階級間隔Δφを揃えることにより、統計的に精度の高い計算結果を得ることができる。
【0112】
第1の相互作用エネルギーφ出現確率計算手段12Eは、第1の相互作用エネルギーφ度数分布作成手段12Dによって作成された度数分布に基づいて、相互作用エネルギーφの各階級の出現確率P
0(φ)を計算する。計算されたP
0(φ)は、制御部12により記憶部13に記憶される。
【0113】
第1の相互作用エネルギーφ出現確率計算手段12Eは、相互作用エネルギーφの各階級の度数を表す度数分布を規格化することにより、相互作用エネルギーφの各階級の出現確率P
0(φ)を計算することができる。ここで、規格化とは、度数分布における相互作用エネルギーφの各階級の度数を、相互作用エネルギーφの各階級の度数の総和によって除することである。
【0114】
第1の相互作用エネルギーε度数分布作成手段12Fは、第2の座標取得手段12Cによって取得された原子集合体ABの座標に基づいて、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、第1の相互作用エネルギーφ度数分布作成手段12Dによって作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。計算された相互作用エネルギーε及び作成された度数分布は、制御部12により記憶部13に記憶される。
【0115】
一実施形態において、第1の相互作用エネルギーε度数分布作成手段12Fは、第2の座標取得手段12Cによって取得された原子集合体ABの座標に基づいて、原子集合体ABの座標R
AB(F
1)〜R
AB(F
i)のそれぞれに関して、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算した後、第1の相互作用エネルギーφ度数分布作成手段12Dによって作成された度数分布に基づいて、相互作用エネルギーφの各階級における、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを抽出し、第1の相互作用エネルギーφ度数分布作成手段12Dによって作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0116】
別の実施形態において、第1の相互作用エネルギーε度数分布作成手段12Fは、第1の相互作用エネルギーφ度数分布作成手段12Dによって作成された度数分布に基づいて、第2の座標取得手段12Cによって取得された原子集合体ABの座標から、相互作用エネルギーφの各階級に属する原子集合体ABの座標を抽出した後、抽出された原子集合体ABの座標に基づいて、相互作用エネルギーφの各階級における、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、第1の相互作用エネルギーφ度数分布作成手段12Dによって作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0117】
原子集合体ABがリガンドと溶媒分子(例えば、水分子)とから構成され、構造a及び該構造aに結合したフラグメントBからなる構造aBが該リガンドである場合には、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεとは、フラグメントBと該フラグメントB周辺に存在する各水分子との相互作用エネルギーを指す。また、原子集合体ABがリガンド、蛋白質及び水分子より構成され、構造a及び該構造aに結合したフラグメントBからなる構造aBが該リガンドである場合には、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεとは、フラグメントBと各水分子との間の相互作用エネルギー及びフラグメントBと蛋白質との間の相互作用エネルギーを指す。ここで挙げられた相互作用エネルギーεの具体例は、第1の相互作用エネルギーε度数分布作成手段12Fによって計算される相互作用エネルギーεだけでなく、第2の相互作用エネルギーε度数分布作成手段12Lによって計算される相互作用エネルギーεについても当てはまる。
【0118】
相互作用エネルギーεは、第1の相互作用エネルギーφ度数分布作成手段12Dによって作成された度数分布(例えばヒストグラム)に含まれる相互作用エネルギーφの全ての階級に関して計算される。相互作用エネルギーεは、例えば、記憶部13に記憶されている第1のシミュレーションプログラム13Cの機能に基づいて計算することができる。相互作用エネルギーεを計算する際、用いられるシミュレーションプログラム及びポテンシャルパラメータは特に限定されるものではないが、原子集合体ABの座標は、第1の座標取得手段12Bが実行するコンピュータシミュレーションにより求められるものであるので、第1の座標取得手段12Bが実行するコンピュータシミュレーションで用いられるものと同一のシミュレーションプログラム及びポテンシャルパラメータを用いることが好ましい。
【0119】
相互作用エネルギーεの各階級の度数を表す度数分布は、相互作用エネルギーφの各階級における、相互作用エネルギーεに関して作成される。相互作用エネルギーεの各階級の度数を表す度数分布は、例えば、横軸が相互作用エネルギーεの各階級を表し、縦軸が相互作用エネルギーεの各階級の度数を表すヒストグラムとして作成することができる。
【0120】
相互作用エネルギーεの各階級は、例えば、相互作用エネルギーεを中心とし、相互作用エネルギー間隔Δεを持った相互作用エネルギー区間[ε−Δε/2〜ε+Δε/2]を階級間隔として作成することができ、各相互作用エネルギー区間におけるΔεは、すべての区間について一定のΔεであってもよいし、相互作用エネルギー区間によって適宜変更したΔεを用いても構わない。相互作用エネルギーεの分割数、つまり階級間隔Δεの個数は、シミュレーションにより得られるスナップショットの数、フラグメントB周囲の水分子等の数、相互作用エネルギーεのエネルギー区間幅等に依るため、特に限定されないが、好ましくは500〜5000、さらに好ましくは2000〜5000である。
【0121】
第1の相互作用エネルギーε出現確率計算手段12Gは、第1の相互作用エネルギーε度数分布作成手段12Fによって作成された度数分布に基づいて、相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の出現確率P
0'(ε;φ)を計算する。計算されたP
0'(ε;φ)は、制御部12により記憶部13に記憶される。なお、P
0'(ε;φ)中の「φ」は、相互作用エネルギーεの各階級の出現確率P
0'(ε;φ)が、相互作用エネルギーφの階級ごとに計算されることを明確にするための表記である。
【0122】
第1の相互作用エネルギーε出現確率計算手段12Gは、相互作用エネルギーεの各階級の度数を表す度数分布を規格化することにより、相互作用エネルギーεの各階級の出現確率P
0'(ε;φ)を計算することができる。ここで、規格化とは、度数分布における相互作用エネルギーεの各階級の度数を、相互作用エネルギーεの各階級の度数の総和によって除することである。
【0123】
第2の原子集合体モデル作成手段12Hは、変化後の原子集合体ABをモデル化した第2の原子集合体モデルを作成する。第2の原子集合体モデル作成手段12Hによって作成された第2の原子集合体モデルに関するデータ(例えば、原子集合体ABを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)は、制御部12により記憶部13に記憶される。
【0124】
第2の原子集合体モデル作成手段12Hは、第2の原子集合体モデルを作成する際、記憶部13に記憶されている変化後の原子集合体ABに関するデータ13Dを用いる。変化後の原子集合体ABに関するデータ13Dは、例えば、制御部12が読み取り可能なファイルの形態で記憶部13に記憶されている。変化後の原子集合体ABに関するデータ13Dは、変化後の原子集合体ABをモデル化した第2の原子集合体モデルを作成可能である限り特に限定されない。変化後の原子集合体ABに関するデータ13Dとしては、例えば、原子集合体ABを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等が挙げられる。
【0125】
第2の原子集合体モデルは、コンピュータシミュレーションが可能である限り特に限定されない。原子集合体モデルとしては、全原子モデル、ビーズスプリングモデル、ユナイテッドアトムモデル等が挙げられる。なお、ビーズスプリングモデルは、原子集合体ABを構成するモノマーユニットを1つのビーズ(セグメント)とみなし、仮想的なバネで結合させたモデルであり、ユナイテッドアトムモデルは、水素原子を重原子(例えば、炭素原子)に含めて1つの原子(質点)として取り扱うモデルである。
【0126】
第3の座標取得手段12Iは、第2の原子集合体モデル作成手段12Hによって作成された第2の原子集合体モデルに対するコンピュータシミュレーションの結果として出力されるスナップショットより、第1〜第jの状態G
1〜G
j(jは2以上の整数である。)の各状態における原子集合体ABの座標を取得する。ここで、スナップショットには、原子集合体ABを構成する全ての原子の座標が含まれる。つまり、状態G
1〜G
jの各状態におけるスナップショットには、それぞれ、状態G
1〜G
jの各状態における原子集合体ABを構成する全ての原子の座標が含まれる。以下、状態G
1、G
2、・・・、G
iにおける原子集合体ABの座標を、それぞれ、座標R
AB(G
1)、R
AB(G
2)、・・・、R
AB(G
j)という場合がある。第3の座標取得手段によって取得された原子集合体ABの座標は、制御部12により記憶部13に記憶される。
【0127】
第3の座標取得手段12Iが実行するコンピュータシミュレーションは、統計力学理論に基づく限り特に限定されない。代表的なコンピュータシミュレーションとしては、例えば、分子動力学法、モンテカルロ法、分子動力学法又はモンテカルロ法と第一原理計算とを組み合わせた方法等が挙げられる。
【0128】
第3の座標取得手段12Iは、記憶部13に記憶されている第2の原子集合体モデルに関するデータ、第2のシミュレーション条件に関するデータ13E、第2のシミュレーションプログラム13F等に基づいてコンピュータシミュレーションを実行する。シミュレーション条件に関するデータ13Eは、制御部12が読み取り可能なファイルの形態で記憶部13に記憶されている。
【0129】
シミュレーション条件は、コンピュータシミュレーションを実行可能である限り特に限定されない。コンピュータシミュレーションに関する具体的な説明は、第1の座標取得手段12Bが実行するコンピュータシミュレーションと同様であるので、省略する。
【0130】
第3の座標取得手段12Iは、コンピュータシミュレーションによって発生するアンサンブルに含まれる状態G
1〜G
jの各状態における原子集合体ABの座標を取得する。状態G
1〜G
jの各状態は、コンピュータシミュレーションの実行により発生する微視的な状態である。状態G
1〜G
jは、コンピュータシミュレーションによって発生する微視的状態のうちの一部であってもよいし、全部であってもよい。
【0131】
jの値は、2以上である限り特に限定されず、計算対象である自由エネルギーの種類(例えば、溶媒和自由エネルギー、結合自由エネルギー等)に応じて適宜選択することができる。計算対象の自由エネルギーが溶媒和自由エネルギーである場合、jの値は、好ましくは10000以上、より好ましくは100000以上である。計算対象の自由エネルギーが結合自由エネルギーである場合、jの値は、好ましくは100000以上、さらに好ましくは1000000以上である。なお、jの上限値は特に限定されないが、通常は10000000、好ましくは5000000である。
【0132】
コンピュータシミュレーションとして分子動力学法を用いる場合、原子集合体ABの状態は、初期状態(時間T
0)から、状態G
1(時間T
1)、状態G
2(時間T
2)、・・・、状態G
j(時間T
j)と経時的に変化する。そして、時系列に沿って、状態G
1〜G
jの各状態における原子集合体ABの座標が取得される。取得された状態G
1〜G
j(時間T
1〜T
j)の各状態(各時間)における原子集合体ABの座標は、それが取得された各状態(各時間)と対応付けられて、制御部12により記憶部13に記憶される。
【0133】
コンピュータシミュレーションとしてモンテカルロ法を用いる場合、乱数を発生させることにより状態G
1〜G
jを作り出し、状態G
1〜G
jの各状態における原子集合体ABの座標が取得される。取得された状態G
1〜G
jの各状態における原子集合体ABの座標は、それが取得された各状態と対応付けられて、制御部12により記憶部13に記憶される。
【0134】
コンピュータシミュレーションとして分子動力学法を用いる場合、第3の座標取得手段12Iは、第2の原子集合体モデルに対して、分子動力学シミュレーションの実行前にエネルギー極小化計算を実行することが好ましい。エネルギー極小化計算は、例えば、記憶部13に記憶されている第2のシミュレーションプログラム13Fの機能に基づいて実行することができる。エネルギー極小化計算によって、原子集合体モデルの初期構造が有する不自然な構造の歪みが除去され、コンピュータシミュレーションの初期段階における時間積分の発散を回避することができる。
【0135】
コンピュータシミュレーションとして分子動力学法を用いる場合、第3の座標取得手段12Iは、第2の原子集合体モデル(好ましくはエネルギー極小化計算後の原子集合体モデル)を平衡化し、平衡化後の状態G
1〜G
jにおける原子集合体ABの座標を取得することが好ましい。例えば、第3の座標取得手段12Iは、第2の原子集合体モデル(好ましくはエネルギー極小化計算後の原子集合体モデル)に対してコンピュータシミュレーションを実行し、例えば、ある物理量が閾値に達した段階で、又は、ある一定時間が経過した段階で、原子集合体モデルが平衡化したと判断し、平衡化後の状態G
1〜G
jにおける原子集合体ABの座標を取得する。
【0136】
第2の相互作用エネルギーφ度数分布作成手段12Jは、第3の座標取得手段12Iによって取得された原子集合体ABの座標に基づいて、構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφを計算し、相互作用エネルギーφの各階級の度数を表す度数分布を作成する。相互作用エネルギーφの各階級は、その計算の基礎とされた原子集合体ABの座標と対応付けられて、制御部12により記憶部13に記憶される。
【0137】
相互作用エネルギーφは、原子集合体ABの座標R
AB(G
1)〜R
AB(G
j)のそれぞれに関して計算される。相互作用エネルギーφは、例えば、第2のシミュレーションプログラム13Fの機能に基づいて計算することができる。相互作用エネルギーφを計算する際、用いられるシミュレーションプログラム及びポテンシャルパラメータは特に限定されるものではないが、原子集合体ABの座標は、第3の座標取得手段12Iが実行するコンピュータシミュレーションにより求められるものであるので、第3の座標取得手段12Iが実行するコンピュータシミュレーションで用いられるものと同一のシミュレーションプログラム及びポテンシャルパラメータを用いることが好ましい。
【0138】
相互作用エネルギーφの各階級の度数を表す度数分布は、例えば、横軸が相互作用エネルギーφの各階級を表し、縦軸が相互作用エネルギーφの各階級の度数を表すヒストグラムとして作成することができる。
【0139】
相互作用エネルギーφの各階級は、例えば、相互作用エネルギーφを中心とし、相互作用エネルギー間隔Δφを持った相互作用エネルギー区間[φ−Δφ/2〜φ+Δφ/2]を階級間隔として作成することができ、各相互作用エネルギー区間におけるΔφは、すべての区間について一定のΔφであってもよいし、相互作用エネルギー区間によって適宜変更したΔφを用いても構わない。数式(1)の第2項の計算を行う上で、P
0(φ)の階級間隔ΔφとP(φ)の階級間隔Δφとを揃えることが好ましく、P
0(φ)の階級間隔Δφに対してP(φ)の階級間隔Δφを揃えることがさらに好ましい。相互作用エネルギーφの分割数、つまり階級間隔Δφの個数は、特に限定されないが、好ましくは50〜500、さらに好ましくは250〜500である。
【0140】
数式(1)の第2項には、相互作用エネルギー出現確率P(φ)を相互作用エネルギー出現確率P
0(φ)によって除する計算過程を含むため、P(φ)が零とは異なる値であって、相互作用エネルギー出現確率P
0(φ)が零となる相互作用エネルギー区間が存在すると、数式(1)の第2項は発散してしまい、数式(1)に基づいた自由エネルギー計算は出来なくなってしまう。そこで、すべての相互作用エネルギー区間[φ−Δφ/2〜φ+Δφ/2]において、P
0(φ)が一定となるように相互作用エネルギー区間毎にΔφを適宜設定することが好ましい。相互作用エネルギー区間毎にΔφを適宜設定し、P
0(φ)の階級間隔Δφに対してP(φ)の階級間隔Δφを揃えることにより、統計的に精度の高い計算結果を得ることができる。
【0141】
第2の相互作用エネルギーφ出現確率計算手段12Kは、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布に基づいて、相互作用エネルギーφの各階級の出現確率P(φ)を計算する。計算されたP(φ)は、制御部12により記憶部13に記憶される。
【0142】
第2の相互作用エネルギーφ出現確率計算手段12Kは、相互作用エネルギーφの各階級の度数を表す度数分布を規格化することにより、相互作用エネルギーφの各階級の出現確率P(φ)を計算することができる。ここで、規格化とは、度数分布における相互作用エネルギーφの各階級の度数を、相互作用エネルギーφの各階級の度数の総和によって除することである。
【0143】
第2の相互作用エネルギーε度数分布作成手段12Lは、第3の座標取得手段12Iによって取得された原子集合体ABの座標に基づいて、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。計算された相互作用エネルギーε及び作成された度数分布は、制御部12により記憶部13に記憶される。
【0144】
第2の相互作用エネルギーε度数分布作成手段12Lは、相互作用エネルギーεを、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布における相互作用エネルギーφの各階級と関連付けることなく計算してもよいし、関連付けて計算してもよい。相互作用エネルギーεを相互作用エネルギーφの各階級と関連付けて計算する場合、相互作用エネルギーεは、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布(例えばヒストグラム)に含まれる相互作用エネルギーφの全ての階級に関して計算され、相互作用エネルギーεの各階級の度数を表す度数分布は、相互作用エネルギーφの各階級における、相互作用エネルギーεに関して作成される。
【0145】
相互作用エネルギーεを相互作用エネルギーφの各階級と関連付けることなく計算する場合、第2の相互作用エネルギーε度数分布作成手段12Lは、第3の座標取得手段12Iによって取得された原子集合体ABの座標に基づいて、原子集合体ABの座標R
AB(G
1)〜R
AB(G
j)のそれぞれに関して、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0146】
相互作用エネルギーεを相互作用エネルギーφの各階級と関連付けて計算する場合、第2の相互作用エネルギーε度数分布作成手段12Lは、第3の座標取得手段12Iによって取得された原子集合体ABの座標に基づいて、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0147】
一実施形態において、第2の相互作用エネルギーε度数分布作成手段12Lは、第3の座標取得手段12Iによって取得された原子集合体ABの座標に基づいて、原子集合体ABの座標R
AB(G
1)〜R
AB(G
j)のそれぞれに関して、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算した後、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布に基づいて、相互作用エネルギーφの各階級における、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを抽出し、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0148】
別の実施形態において、第2の相互作用エネルギーε度数分布作成手段12Lは、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布に基づいて、第3の座標取得手段12Iによって取得された原子集合体ABの座標から、相互作用エネルギーφの各階級に属する原子集合体ABの座標を抽出した後、抽出された原子集合体ABの座標に基づいて、相互作用エネルギーφの各階級における、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0149】
相互作用エネルギーεは、例えば、第2のシミュレーションプログラム13Fの機能に基づいて計算することができる。相互作用エネルギーεを計算する際、用いられるシミュレーションプログラム及びポテンシャルパラメータは特に限定されるものではないが、原子集合体ABの座標は、第3の座標取得手段12Iが実行するコンピュータシミュレーションにより求められるものであるので、第3の座標取得手段12Iが実行するコンピュータシミュレーションで用いられるものと同一のシミュレーションプログラム及びポテンシャルパラメータを用いることが好ましい。
【0150】
相互作用エネルギーεの各階級の度数を表す度数分布は、例えば、横軸が相互作用エネルギーεの各階級を表し、縦軸が相互作用エネルギーεの各階級の度数を表すヒストグラムとして作成することができる。
【0151】
相互作用エネルギーεの各階級は、例えば、相互作用エネルギーεを中心とし、相互作用エネルギー間隔Δεを持った相互作用エネルギー区間[ε−Δε/2〜ε+Δε/2]を階級間隔として作成することができ、各相互作用エネルギー区間におけるΔεは、すべての区間について一定のΔεであってもよいし、相互作用エネルギー区間によって適宜変更したΔεを用いても構わない。また、P
0'(ε;φ)における階級間隔ΔεとP'(ε)における階級間隔Δεは揃えることが好ましい。相互作用エネルギーεの分割数、つまり階級間隔Δεの個数は、シミュレーションにより得られるスナップショットの数、フラグメントB周囲の水分子の数、相互作用エネルギーεのエネルギー区間幅等に依るため、特に限定されないが、好ましくは500〜5000、さらに好ましくは2000〜5000である。
【0152】
第2の相互作用エネルギーε出現確率計算手段12Mは、第2の相互作用エネルギーε度数分布作成手段12Lによって作成された度数分布に基づいて、相互作用エネルギーεの各階級の出現確率P'(ε)を計算する。計算されたP'(ε)は、制御部12により記憶部13に記憶される。
【0153】
第2の相互作用エネルギーε出現確率計算手段12Mは、相互作用エネルギーεの各階級の度数を表す度数分布を規格化することにより、相互作用エネルギーεの各階級の出現確率P'(ε)を計算することができる。ここで、規格化とは、度数分布における相互作用エネルギーεの各階級の度数を、相互作用エネルギーεの各階級の度数の総和によって除することである。
【0154】
第2の相互作用エネルギーε度数分布作成手段12Lが、相互作用エネルギーεを、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布における相互作用エネルギーφの各階級と関連付けることなく計算する場合、第2の相互作用エネルギーε出現確率計算手段12Mは、第2の相互作用エネルギーε度数分布作成手段12Lによって作成された度数分布に基づいて、相互作用エネルギーφの各階級と関連付けることなく、相互作用エネルギーεの各階級の出現確率P'(ε)を計算する。
【0155】
第2の相互作用エネルギーε度数分布作成手段12Lが、相互作用エネルギーεを、第2の相互作用エネルギーφ度数分布作成手段12Jによって作成された度数分布における相互作用エネルギーφの各階級と関連付けて計算する場合、第2の相互作用エネルギーε出現確率計算手段12Mは、第2の相互作用エネルギーε度数分布作成手段12Lによって作成された度数分布に基づいて、相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の出現確率P'(ε;φ)を計算する。なお、P'(ε;φ)中の「φ」は、相互作用エネルギーεの各階級の出現確率P'(ε;φ)が、相互作用エネルギーφの階級ごとに計算されることを明確にするための表記である。
【0156】
∫Δν(φ)P(φ)dφ計算手段12Nは、第2の相互作用エネルギーφ出現確率計算手段12Kによって計算されたP(φ)と、第1の相互作用エネルギーε出現確率計算手段12Gによって計算されたP
0'(ε;φ)と、第2の相互作用エネルギーε出現確率計算手段12Mによって計算されたP'(ε)とに基づいて、相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφを計算する。
【0157】
「∫Δν(φ)P(φ)dφ」中の「Δν(φ)」は、相互作用エネルギーφの各階級における、相互作用エネルギーεに起因する自由エネルギー変化量を表す。原子集合体ABがリガンドと溶媒分子(例えば、水分子)とから構成され、構造a及び該構造aに結合したフラグメントBからなる構造aBが該リガンドである場合、相互作用エネルギーεに起因する自由エネルギー変化量とは、フラグメントBと該フラグメントB周辺に存在する水分子との相互作用エネルギーに起因する自由エネルギー変化量を指す。また、原子集合体ABがリガンド、蛋白質及び水分子より構成され、構造a及び該構造aに結合したフラグメントBからなる構造aBが該リガンドである場合、相互作用エネルギーεに起因する自由エネルギー変化量とは、フラグメントBと水分子との間の相互作用エネルギーに起因する自由エネルギー変化量と、フラグメントBと蛋白質との間の相互作用エネルギーに起因する自由エネルギー変化量との和を指す。
【0158】
∫Δν(φ)P(φ)dφ計算手段12Nは、好ましくは、エネルギー表示法により、第2の相互作用エネルギーφ出現確率計算手段12Kによって計算されたP(φ)と、第1の相互作用エネルギーε出現確率計算手段12Gによって計算されたP
0'(ε;φ)と、第2の相互作用エネルギーε出現確率計算手段12Mによって計算されたP'(ε)とに基づいて、相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφを計算する。
【0159】
一実施形態において、∫Δν(φ)P(φ)dφ計算手段12Nは、Δν(φ)を計算することなく、第2の相互作用エネルギーφ出現確率計算手段12Kによって計算されたP(φ)と、第1の相互作用エネルギーε出現確率計算手段12Gによって計算されたP
0'(ε;φ)と、第2の相互作用エネルギーε出現確率計算手段12Mによって計算されたP'(ε)とに基づいて、エネルギー表示法により、∫Δν(φ)P(φ)dφを計算する。ここで用いられるP'(ε)は、相互作用エネルギーφの各階級と関連付けられることなく計算された相互作用エネルギーεの各階級の出現確率である。
【0160】
別の実施形態において、∫Δν(φ)P(φ)dφ計算手段12Nは、第1の相互作用エネルギーε出現確率計算手段12Gによって計算されたP
0'(ε;φ)と、第2の相互作用エネルギーε出現確率計算手段12Mによって計算されたP'(ε;φ)とに基づいて、相互作用エネルギーφの各階級における、相互作用エネルギーεに起因する自由エネルギー変化量Δν(φ)を計算し、計算されたΔν(φ)と、第2の相互作用エネルギーφ出現確率計算手段12Kによって計算されたP(φ)とに基づいて、∫Δν(φ)P(φ)dφを計算する。この際、∫Δν(φ)P(φ)dφ計算手段12Nは、好ましくは、エネルギー表示法により、第1の相互作用エネルギーε出現確率計算手段12Gによって計算されたP
0'(ε;φ)と、第2の相互作用エネルギーε出現確率計算手段12Mによって計算されたP'(ε;φ)とに基づいて、相互作用エネルギーφの各階級における、相互作用エネルギーεに起因する自由エネルギー変化量Δν(φ)を計算する。
【0161】
ΔG計算手段12Oは、第1の相互作用エネルギーφ出現確率計算手段12Eによって計算されたP
0(φ)と、第2の相互作用エネルギーφ出現確率計算手段12Kによって計算されたP(φ)と、∫Δν(φ)P(φ)dφ計算手段12Nによって計算された∫Δν(φ)P(φ)dφと、数式(1):
【数10】
[式中、Rは気体定数を表し、Tは反応式(1)で表される変化が生じる際の絶対温度を表す。]
とに基づいて、ΔGを計算する。
【0162】
以下、本発明の一実施形態に係る計算処理手順について、図面に基づいて説明する。
図2は、本実施形態に係る計算装置1の処理手順を示すフローチャートである。
制御部12は、変化前の原子集合体Aに関する計算処理(ステップS110〜190)と、変化後の原子集合体ABに関する計算処理(ステップS210〜280)とを実行する際、一方を実行した後、他方を実行してもよいし、両方を並行して実行してもよい。
【0163】
以下、ステップS110〜S150について説明する。
制御部12は、ステップS110〜S150を順次実施する。
ステップS110において、制御部12は、記憶部13に記憶されている変化前の原子集合体Aに関するデータ13Aを用いて、変化前の原子集合体Aをモデル化した第1の原子集合体モデルを作成する。作成された第1の原子集合体モデルに関するデータは、制御部12により記憶部13に記憶される。ステップS110において、制御部12は、第1の原子集合体モデル作成手段12Aとして機能する。
【0164】
ステップS120において、制御部12は、記憶部13に記憶されている第1の原子集合体モデルに関するデータ、第1のシミュレーション条件に関するデータ13B、第1のシミュレーションプログラム13C等を用いて、第1の原子集合体モデルに対するコンピュータシミュレーションにより、第1〜第iの状態F
1〜F
i(iは2以上の整数である。
)の各状態における原子集合体Aの座標を取得する。ステップS120において、制御部12は、第1の座標取得手段12Bとして機能する。
【0165】
コンピュータシミュレーションとして分子動力学法を用いる場合に制御部12が実行するステップS120の一実施形態を
図3に示す。この実施形態において、制御部12は、ステップS110の後、第1の原子集合体モデルに対してエネルギー極小化計算を実行し(ステップS121)、次いで、第1の原子集合体モデルを平衡化し(ステップS122)、平衡化後の第1の原子集合体モデルに対するコンピュータシミュレーションにより、平衡化後の状態F
1〜F
iの各状態における原子集合体Aの座標を取得する。ステップS122において、制御部12は、例えば、第1の原子集合体モデルにおけるある物理量が閾値に達した段階で、又は、コンピュータシミュレーションの開始から一定時間が経過した段階で、第1の原子集合体モデルが平衡化したと判断する。
【0166】
ステップS130において、制御部12は、ステップS120で取得された原子集合体Aの座標に基づいて、状態F
1〜F
iの各状態における原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABの座標を取得する。ステップS130において、制御部12は、第2の座標取得手段12Cとして機能する。
【0167】
制御部12が実行するステップS130の一実施形態を
図4に示す。この実施形態において、制御部12は、ステップS120の後、以下のステップ:
状態F
xにおける原子集合体Aに結合させるべきフラグメントBに関し、フラグメントBを構成する原子から選択された原子bの、構造aに対する相対位置条件(例えば、構造aを構成する1個又は2個以上の原子に対する距離、角度等)を規定するステップS131a;
規定された相対位置条件に基づいて、原子bの座標を取得するステップS132a;
取得された原子bの座標と、状態F
xにおける原子集合体Aの座標R
A(F
x)とに基づいて、状態F
xにおける原子集合体Aに対して原子bが結合することにより生じる原子集合体Abの座標R
Ab(F
x)を取得するステップS133a;並びに
原子bの座標を基準として、コンピュータシミュレーション、公知の配座発生方法等により、フラグメントBを構成する原子b以外の原子の座標を発生させ、状態F
xにおける原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABの座標R
AB(F
x)を取得するステップS134aを実行する。そして、ステップS135aにおいて、制御部12は、ステップS131a〜S134aがx=1,2,・・・,iで実行されたか否かを判定し、「NO」の場合にはステップS131a〜S134aを繰り返す。すなわち、制御部12は、ステップS135aで「YES」と判定されるまで、x=1,2,・・・,iで繰り返しステップS131a〜S134aを実行する。これにより、状態F
1〜F
iの各状態における原子集合体Aに対して、構造aに対する相対位置が限定されたフラグメントBが結合することによって生じる原子集合体ABの座標R
AB(F
1)〜R
AB(F
i)を取得することができる。
【0168】
ステップS131aにおいて、制御部12は、状態F
xにおける原子集合体Aに結合させるべきフラグメントBに関し、互いに異なる第1、第2、・・・、第pの相対位置条件を規定してもよい。pは、xの値から独立して選択された、2以上の整数である。したがって、pはx=1,2,・・・,iにおいて同一の値であってもよいし、x=1,2,・・・,iにおいて異なる値であってもよい。
【0169】
ステップS131aにおいて、状態F
xにおける原子集合体Aに結合させるべきフラグメントBに関し、互いに異なる第1、第2、・・・、第pの相対位置条件が規定される場合、制御部12は、第1〜第pの相対位置条件のそれぞれに関し、ステップS132a〜S134aを繰り返す。これにより、状態F
xにおける原子集合体Aに対して、構造aに対する相対位置が異なるフラグメントB
1、B
2、・・・、B
pが結合することにより生じる原子集合体ABの座標R
AB(F
x,B
1)、R
AB(F
x,B
2)、・・・、R
AB(F
x,B
p)を取得することができる。
【0170】
制御部12が実行するステップS130の好ましい実施形態を
図5に示す。この実施形態において、制御部12は、ステップS120の後、以下のステップ:
構造aと該構造aに結合したフラグメントBとからなる又は構造aと該構造aに結合したフラグメントBとを含んでなる原子集合体Cをモデル化した第3の原子集合体モデルを作成するステップS131b;並びに、
作成された第3の原子集合体モデルに対するコンピュータシミュレーションにより、第1〜第kの状態H
1〜H
k(kは2以上の整数である。)の各状態における原子集合体Cの座標を取得するステップS132b
を実行した後、
構造aを構成する原子から選択された1個又は2個以上の原子からなる選択原子群(なお、2個以上の原子が選択される場合だけでなく、1個の原子が選択される場合も含めて便宜上「選択原子群」という。)を選択して、状態F
xにおける原子集合体Aの選択原子群の座標に対し、状態H
1〜H
kから選択された1又は2以上の状態における原子集合体Cの選択原子群の座標を回転及び/又は並進させることにより、原子集合体Aの選択原子群と原子集合体Cの選択原子群との間で、対応する原子間の距離の二乗和が最小となる原子集合体Cの座標を作成(原子集合体Aの選択原子群の座標に原子集合体Cの選択原子群の座標をフィッティング)し、作成された原子集合体Cの座標に基づいて、原子集合体Aに状態H
1〜H
kから選択された1又は2以上の状態における原子集合体Cを重ね合わせるステップS133b;並びに
状態F
xにおける原子集合体Aの座標と、該原子集合体Aに対して重ね合わせた原子集合体CのフラグメントBの1又は2以上の座標とに基づいて、状態F
xにおける原子集合体Aに対してフラグメントBが結合することにより生じる原子集合体ABの座標を取得するステップS134b
を、x=1〜iで繰り返して実行する。そして、ステップS135bにおいて、制御部12は、ステップS133b〜S134bがx=1,2,・・・,iで実行されたか否かを判定し、「NO」の場合にはステップS133b〜S134bを繰り返す。すなわち、制御部12は、ステップS135bで「YES」と判定されるまで、x=1,2,・・・,iで繰り返しステップS133b〜S134bを実行する。これにより、状態F
1〜F
iの各状態における原子集合体Aに対して、構造aに対する相対位置が限定されたフラグメントBが結合することによって生じる原子集合体ABの座標R
AB(F
1)〜R
AB(F
i)を取得することができる。
【0171】
ステップS133bにおいて、制御部12は、状態H
1〜H
kにおける原子集合体Cの座標から、構造aに対するフラグメントBの相対位置が互いに異なるp組の座標を選択してもよい。pは、xの値から独立して選択された、2以上の整数である。したがって、pはx=1,2,・・・,iにおいて同一の値であってもよいし、x=1,2,・・・,iにおいて異なる値であってもよい。
【0172】
ステップS133bにおいて、状態H
1〜H
kにおける原子集合体Cの座標から、構造aに対するフラグメントBの相対位置が互いに異なるp組の座標が選択される場合、制御部12は、p組の座標のそれぞれに関し、ステップS133b〜S134bを繰り返す。これにより、状態F
xにおける原子集合体Aに対して、構造aに対する相対位置が異なるフラグメントB
1、B
2、・・・、B
pが結合することにより生じる原子集合体ABの座標R
AB(F
x,B
1)、R
AB(F
x,B
2)、・・・、R
AB(F
x,B
p)を取得することができる。
【0173】
ステップS140において、制御部12は、ステップS130で取得された原子集合体ABの座標に基づいて、構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφを計算し、ステップS150において、相互作用エネルギーφの各階級の度数を表す度数分布を作成する。ステップS140及びS150において、制御部12は、第1の相互作用エネルギーφ度数分布作成手段12Dとして機能する。
【0174】
以下、ステップS160〜S190について説明する。
ステップS150の後、制御部12は、相互作用エネルギーφに関する処理(ステップS160)と、相互作用エネルギーεに関する処理(ステップS170〜S190)とを実行する。この際、制御部12は、ステップS160を実行した後、ステップS170〜S190を実行してもよいし、ステップS170〜S190を実行した後、ステップS160を実行してもよいし、ステップS160とステップS170〜S190とを並行して実行してもよい。いずれの場合も、制御部12は、ステップS170〜S190を順次実施する。
【0175】
ステップS160において、制御部12は、ステップS150で作成された度数分布に基づいて、相互作用エネルギーφの各階級の出現確率P
0(φ)を計算する。ステップS160において、制御部12は、第1の相互作用エネルギーφ出現確率計算手段12Eとして機能する。
【0176】
ステップS170において、制御部12は、ステップS130で取得された原子集合体ABの座標に基づいて、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、ステップS180において、ステップS150で作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。ステップS170及びS180において、制御部12は、第1の相互作用エネルギーε度数分布作成手段12Fとして機能する。
【0177】
ステップS170において、制御部12は、ステップS130で取得された原子集合体ABの座標に基づいて、原子集合体ABの座標R
AB(F
1)〜R
AB(F
i)のそれぞれに関して、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算した後、ステップS130で取得された原子集合体ABの座標とステップS150で作成された度数分布とに基づいて、相互作用エネルギーφの各階級における、原子集合体ABから該構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを抽出し、ステップS150で作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成してもよいし、ステップS150で作成された度数分布に基づいて、ステップS130で取得された原子集合体ABの座標から、相互作用エネルギーφの各階級に属する原子集合体ABの座標を抽出した後、抽出された原子集合体ABの座標に基づいて、相互作用エネルギーφの各階級における、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、ステップS150で作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成してもよい。
【0178】
ステップS190において、制御部12は、ステップS180で作成された度数分布に基づいて、相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の出現確率P
0'(ε;φ)を計算する。ステップS190において、制御部12は、第1の相互作用エネルギーε出現確率計算手段12Gとして機能する。
【0179】
以下、ステップS210〜S240について説明する。
制御部12は、ステップS210〜S240を順次実施する。
ステップS210において、制御部12は、記憶部13に記憶されている変化後の原子集合体ABに関するデータ13Dを用いて、変化後の原子集合体ABをモデル化した第2の原子集合体モデルを作成する。作成された第2の原子集合体モデルに関するデータは、制御部12により記憶部13に記憶される。ステップS210において、制御部12は、第2の原子集合体モデル作成手段12Hとして機能する。
【0180】
ステップS220において、制御部12は、記憶部13に記憶されている第2の原子集合体モデルに関するデータ、第2のシミュレーション条件に関するデータ13E、第2のシミュレーションプログラム13F等を用いて、第2の原子集合体モデルに対するコンピュータシミュレーションにより、第1〜第jの状態G
1〜G
j(jは2以上の整数である。
)の各状態における原子集合体ABの座標を取得する。ステップS220において、制御部12は、第3の座標取得手段12Iとして機能する。
【0181】
コンピュータシミュレーションとして分子動力学法を用いる場合に制御部12が実行するステップS220の一実施形態を
図6に示す。この実施形態において、制御部12は、ステップS210の後、第2の原子集合体モデルに対してエネルギー極小化計算を実行し(ステップS221)、次いで、第2の原子集合体モデルを平衡化し(ステップS222)、平衡化後の第2の原子集合体モデルに対するコンピュータシミュレーションにより、平衡化後の状態G
1〜G
jの各状態における原子集合体ABの座標を取得する。ステップS222において、制御部12は、例えば、第2の原子集合体モデルにおけるある物理量が閾値に達した段階で、又は、コンピュータシミュレーションの開始から一定時間が経過した段階で、第2の原子集合体モデルが平衡化したと判断する。
【0182】
ステップS230において、制御部12は、ステップS220で取得された原子集合体ABの座標に基づいて、構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφを計算し、ステップS240において、相互作用エネルギーφの各階級の度数を表す度数分布を作成する。ステップS230及びS240において、制御部12は、第2の相互作用エネルギーφ度数分布作成手段12Jとして機能する。
【0183】
以下、ステップS250〜S280について説明する。
ステップS240の後、制御部12は、相互作用エネルギーφに関する処理(ステップS250)と、相互作用エネルギーεに関する処理(ステップS260〜S280)とを実行する。この際、制御部12は、ステップS250を実行した後、ステップS260〜S280を実行してもよいし、ステップS260〜S280を実行した後、ステップS250を実行してもよいし、ステップS250とステップS260〜S280とを並行して実行してもよい。いずれの場合も、制御部12は、ステップS260〜S280を順次実施する。
【0184】
ステップS250において、制御部12は、ステップS240で作成された度数分布に基づいて、相互作用エネルギーφの各階級の出現確率P(φ)を計算する。ステップS250において、制御部12は、第2の相互作用エネルギーφ出現確率計算手段12Kとして機能する。
【0185】
ステップS260において、制御部12は、ステップS220で取得された原子集合体ABの座標に基づいて、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、ステップS270において、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。ステップS260及びS270において、制御部12は、第2の相互作用エネルギーε度数分布作成手段12Lとして機能する。
【0186】
ステップS260において、制御部12は、相互作用エネルギーεを、ステップS240で作成された度数分布における相互作用エネルギーφの各階級と関連付けることなく計算してもよいし、関連付けて計算してもよい。相互作用エネルギーεが相互作用エネルギーφの各階級と関連付けられて計算される場合、相互作用エネルギーεは、ステップS240で作成された度数分布(例えばヒストグラム)に含まれる相互作用エネルギーφの全ての階級に関して計算され、相互作用エネルギーεの各階級の度数を表す度数分布は、相互作用エネルギーφの各階級における、相互作用エネルギーεに関して作成される。
【0187】
制御部12が、相互作用エネルギーεを相互作用エネルギーφの各階級と関連付けることなく計算する場合、制御部12は、ステップS220で取得された原子集合体ABの座標に基づいて、原子集合体ABの座標R
AB(G
1)〜R
AB(G
j)のそれぞれに関して、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0188】
制御部12が、相互作用エネルギーεを相互作用エネルギーφの各階級と関連付けて計算する場合、制御部12は、ステップS220で取得された原子集合体ABの座標に基づいて、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、ステップS240で作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0189】
一実施形態において、制御部12は、ステップS220で取得された原子集合体ABの座標に基づいて、原子集合体ABの座標R
AB(G
1)〜R
AB(G
j)のそれぞれに関して、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算した後、ステップS240で作成された度数分布に基づいて、相互作用エネルギーφの各階級における、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを抽出し、ステップS240で作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0190】
別の実施形態において、制御部12は、ステップS240で作成された度数分布に基づいて、ステップS220で取得された原子集合体ABの座標から、相互作用エネルギーφの各階級に属する原子集合体ABの座標を抽出した後、抽出された原子集合体ABの座標に基づいて、相互作用エネルギーφの各階級における、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、ステップS240で作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0191】
制御部12が、相互作用エネルギーεを相互作用エネルギーφの各階級と関連付けることなく計算する場合、ステップS260及びS270において、制御部12は、ステップS220で取得された原子集合体ABの座標に基づいて、原子集合体ABの座標R
AB(G
1)〜R
AB(G
j)のそれぞれに関して、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、相互作用エネルギーεの各階級の度数を表す度数分布を作成する。
【0192】
制御部12が、相互作用エネルギーεを相互作用エネルギーφの各階級と関連付けて計算する場合、ステップS260及びS270において、制御部12は、ステップS220で取得された原子集合体ABの座標に基づいて、原子集合体ABの座標R
AB(G
1)〜R
AB(G
j)のそれぞれに関して、原子集合体ABから構造a及び該構造aに結合したフラグメントBからなる構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算した後、ステップS240で作成された度数分布に基づいて、相互作用エネルギーφの各階級における、原子集合体ABから構造aBを除くことにより生じる原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを抽出し、ステップS240で作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成してもよいし、ステップS240で作成された度数分布に基づいて、ステップS220で取得された原子集合体ABの座標から、相互作用エネルギーφの各階級に属する原子集合体ABの座標を抽出した後、抽出された原子集合体ABの座標に基づいて、相互作用エネルギーφの各階級における、原子集合体ABから構造aBを除いた原子集合体の一部又は全部と、フラグメントBとの相互作用エネルギーεを計算し、ステップS240で作成された度数分布の相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数を表す度数分布を作成してもよい。
【0193】
ステップS280において、制御部12は、ステップS270で作成された度数分布に基づいて、相互作用エネルギーεの各階級の出現確率P'(ε)を計算する。ステップS280において、制御部12は、第2の相互作用エネルギーε出現確率計算手段12Mとして機能する。
【0194】
ステップS260及びS270において、制御部12が、相互作用エネルギーεを、ステップS240で作成された度数分布における相互作用エネルギーφの各階級と関連付けることなく計算し、相互作用エネルギーεの各階級の度数分布を作成する場合、ステップS280において、制御部12は、ステップS270で作成された度数分布に基づいて、相互作用エネルギーφの各階級と関連付けることなく、相互作用エネルギーεの各階級の出現確率P'(ε)を計算する。
【0195】
ステップS260及びS270において、制御部12が、相互作用エネルギーεを、ステップS240で作成された度数分布における相互作用エネルギーφの各階級と関連付けて計算し、相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数分布を作成する場合、ステップS280において、制御部12は、ステップS270で作成された度数分布に基づいて、相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の出現確率P'(ε;φ)を計算する。
【0196】
以下、ステップS300について説明する。
制御部12は、変化前の原子集合体Aに関する計算処理(ステップS110〜S190)と、変化後の原子集合体ABに関する計算処理(ステップS210〜S280)とを実行した後、自由エネルギー変化量∫Δν(φ)P(φ)dφに関する計算処理(ステップS300)を実行する。
【0197】
ステップS300において、制御部12は、ステップS250で計算されたP(φ)と、ステップS190で計算されたP
0'(ε;φ)と、ステップS280で計算されたP'(ε)とに基づいて、相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφを計算する。ステップS300において、制御部12は、∫Δν(φ)P(φ)dφ計算手段12Nとして機能する。
【0198】
一実施形態において、制御部12は、Δν(φ)を計算することなく、ステップS250で計算されたP(φ)と、ステップS190で計算されたP
0'(ε;φ)と、ステップS280で計算されたP'(ε)とに基づいて、エネルギー表示法により、∫Δν(φ)P(φ)dφを計算する。ここで用いられるP'(ε)は、相互作用エネルギーφの各階級と関連付けられることなく計算された相互作用エネルギーεの各階級の出現確率である。
【0199】
別の実施形態において、制御部12は、ステップS190で計算されたP
0'(ε;φ)と、ステップS280で計算されたP'(ε;φ)とに基づいて、相互作用エネルギーφの各階級における、相互作用エネルギーεに起因する自由エネルギー変化量Δν(φ)を計算し、計算されたΔν(φ)と、ステップS250で計算されたP(φ)とに基づいて、∫Δν(φ)P(φ)dφを計算する。この際、∫Δν(φ)P(φ)dφ計算手段12Nは、好ましくは、エネルギー表示法により、ステップS190で計算されたP
0'(ε;φ)と、ステップS280で計算されたP'(ε;φ)とに基づいて、相互作用エネルギーφの各階級における、相互作用エネルギーεに起因する自由エネルギー変化量Δν(φ)を計算する。
【0200】
以下、ステップS400について説明する。
制御部12は、自由エネルギー変化量∫Δν(φ)P(φ)dφに関する計算処理(ステップS300)を実行した後、ΔGに関する計算処理(ステップS400)を実行する。
【0201】
ステップS400において、制御部12は、ステップS160で計算されたP
0(φ)と、ステップS250で計算されたP(φ)と、ステップS300で計算された∫Δν(φ)P(φ)dφと、数式(1):
【数11】
[式中、Rは気体定数を表し、Tは反応式(1)で表される変化が生じる際の絶対温度を表す。]
とに基づいて、ΔGを計算する。
【0202】
本実施形態に係る計算装置1の機能は、制御部12を、第1の原子集合体モデル作成手段12A、第1の座標取得手段12B、第2の座標取得手段12C、第1の相互作用エネルギーφ度数分布作成手段12D、第1の相互作用エネルギーφ出現確率計算手段12E、第1の相互作用エネルギーε度数分布作成手段12F、第1の相互作用エネルギーε出現確率計算手段12G、第2の原子集合体モデル作成手段12H、第3の座標取得手段12I、第2の相互作用エネルギーφ度数分布作成手段12J、第2の相互作用エネルギーφ出現確率計算手段12K、第2の相互作用エネルギーε度数分布作成手段12L、第2の相互作用エネルギーε出現確率計算手段12M、∫Δν(φ)P(φ)dφ計算手段12N、及びΔG計算手段12Oとして機能させるプログラム(すなわち、制御部12に、ステップS110〜S190、ステップS210〜S280、ステップS300及びステップS400を実行させるプログラム)が記録されたコンピュータ読み取り可能な記録媒体を用いて、コンピュータに該プログラムをインストールすることにより実現することができる。プログラムが記録されたコンピュータ読み取り可能な記録媒体としては、例えば、ROM、フロッピーディスク(登録商標)、ハードディスク、光ディスク、光磁気ディスク、CD−ROM、磁気テープ、不揮発性メモリカード等が挙げられる。
【0203】
なお、本発明は、コンピュータがプログラムを実行することよって本実施形態に係る計算装置1の機能を実現する実施形態だけでなく、プログラムがコンピュータにおいて稼働しているOS(オペレーティングシステム)、他のアプリケーションソフト等と共同することによって本実施形態に係る計算装置1の機能を実現する実施形態も包含する。また、本発明は、プログラムがコンピュータの機能拡張ボード又はコンピュータに接続された機能拡張ユニットに備わるメモリに格納された後、該プログラムの指示に基づいて、該機能拡張ボード又は該機能拡張ユニットに備わるCPU等が実際の処理の一部又は全部を実行することによって本実施形態に係る計算装置1の機能を実現する実施形態も包含する。
【実施例】
【0204】
本発明の具体的実施態様を以下に実施例をもって述べるが、本発明はこれに限定されるものものではない。
本実施例では、反応式(2)及び(3):
【化1】
において、ΔG2−ΔG1で定義されるΔΔGを計算した。
【0205】
反応式(2)及び(3)中、Rは、蛋白質を表し、Rを囲む□は、蛋白質Rの周囲に水分子が存在することを表し、L
1は、蛋白質Rと結合するリガンドを表し、L
1を囲む□は、リガンドL
1の周囲に水分子が存在することを表し、L
1Rは、蛋白質Rと該蛋白質Rに結合したリガンドL
1との複合体を表し、L
1Rを囲む□は、複合体L
1Rの周囲に水分子が存在することを表し、L
2は、蛋白質Rと結合するリガンドであって、リガンドL
1とは異なるリガンドを表し、L
2を囲む□は、リガンドL
2の周囲に水分子が存在することを表し、L
2Rは、蛋白質Rと該蛋白質Rに結合したリガンドL
2との複合体を表し、L
2Rを囲む□は、複合体L
2Rの周囲に水分子が存在することを表す。なお、反応式(2)の左辺において、L
1とRとは相互作用しておらず、反応式(3)の左辺において、L
2とRとは相互作用していない。
【0206】
ΔΔGは、反応式(4):
【化2】
より、ΔΔG=ΔG4−ΔG3と表すことができる。反応式(4)の左辺において、蛋白質RとリガンドL
1及びL
2とは相互作用していないため、ΔG3は、周囲に水分子が存在した状態にあるリガンドL
1と、周囲に水分子が存在した状態にあるリガンドL
2との自由エネルギー差である。
【0207】
ΔG4は、反応式(5):
【化3】
より、ΔG4=ΔG4a+ΔG4bと表すことができる。なお、反応式(5)中、L
1R→L
0Rの変化において除去されるフラグメントと、L
0R→L
2Rの変化において付加されるフラグメントの両方を便宜上「B」で表すが、実際には、L
1R→L
0Rの変化において除去されるフラグメントと、L
0R→L
2Rの変化において付加されるフラグメントとは異なる。L
1R→L
0Rの変化において除去されるフラグメント及びL
0R→L
2Rの変化において付加されるフラグメントの詳細については、後述する。
【0208】
反応式(5)中、L
0Rは、蛋白質Rと該蛋白質Rに結合したリガンドL
0との複合体を表し、L
0Rを囲む□は、複合体L
0Rの周囲に水分子が存在することを表す。
【0209】
ΔG3は、反応式(6):
【化4】
より、ΔG3=ΔG3a+ΔG3bと表すことができる。なお、反応式(6)中、L
1→L
0の変化において除去されるフラグメントと、L
0→L
2の変化において付加されるフラグメントの両方を便宜上「B」で表すが、実際には、L
1→L
0の変化において除去されるフラグメントと、L
0→L
2の変化において付加されるフラグメントとは異なる。L
1→L
0の変化において除去されるフラグメント及びL
0→L
2の変化において付加されるフラグメントの詳細については、後述する。
【0210】
反応式(6)中、L
0は、蛋白質Rと結合するリガンドであって、リガンドL
1及びリガンドL
2とは異なるリガンドを表し、L
0を囲む□は、リガンドL
0の周囲に水分子が存在することを表す。なお、反応式(6)中、Rは、L
1、L
2、L
0のいずれのリガンドとも相互作用していない。よって、反応式(6)中のΔG3は、周囲に水分子が存在した状態にあるリガンドL
1と、周囲に水分子が存在した状態にあるリガンドL
2との自由エネルギー差を意味する。
【0211】
〔実施例1〕
(1)ΔG4aの計算
反応式(7):
【化5】
で表される変化に関し、複合体L
0R及び該複合体L
0Rの周囲に存在する水分子を構成する変化前の原子集合体を「原子集合体A」とし、複合体L
1R及び該複合体L
1Rの周囲に存在する水分子を構成する変化後の原子集合体を「原子集合体AB」とし、変化前の原子集合体Aの自由エネルギー及びフラグメントBの自由エネルギーの和と、変化後の原子集合体ABの自由エネルギーとの差ΔG4a’を数式(1)に基づいて計算した。
【0212】
なお、ΔG4aは、反応式(8):
【化6】
で表される変化に関する自由エネルギー変化量であるので、数式(1)に基づいて計算されるΔG4a’にマイナスを乗じたものがΔG4aに相当する。
【0213】
蛋白質Rとして、セリンプロテアーゼの1種であるトリプシン(Trypsin)を選択した。
【0214】
リガンドL
1として、トリプシン阻害剤であるベンズアミジンを選択した。リガンドL
1の化学構造は、以下の通りである。
【化7】
【0215】
リガンドL
1を構成する各原子の部分電荷は、以下の通りである。部分電荷は、小数点以下第4位を四捨五入して、小数点以下第3位までを記載している。以後に記載の部分電荷も同様である。なお、これら部分電荷は、後述するMATCHプログラムにより作成されたものである。
【化8】
【0216】
リガンドL
0として、リガンドL
1から、アミジン基のパラ位に位置するベンゼン環の炭素原子に結合する水素原子と、該水素原子が結合するベンゼン環の炭素原子の部分電荷(−0.115)を持つ点電荷とからなるフラグメントBが除去された分子を選択した。リガンドL
0の化学構造は、以下の通りである。
【化9】
【0217】
リガンドL
0を構成する各原子の部分電荷は、以下の通りである。
【化10】
【0218】
フラグメントBとして、水素原子及び点電荷(部分電荷;−0.115)を選択した。水素原子は、リガンドL
1において、ベンゼン環に結合していた水素原子であり、点電荷は、リガンドL
1において、水素原子が結合していた炭素原子が有していた部分電荷(−0.115)を持つ。リガンドL
1から、フラグメントBが除去されることにより、ベンゼン環を構成する炭素原子のうち、フラグメントBの水素原子が結合していた炭素原子の部分電荷が零である分子がリガンドL
0として得られる。フラグメントBの構造は以下の通りである。ここで、部分電荷(−0.115)を持つ点電荷と、水素原子とは、結合していない、つまり、結合長、結合角、捻れ角等に起因する相互作用はない。よって、以下では、仮想原子と水素原子との間は、破線にて図示している。また、フラグメントBの点電荷と、点電荷以外のフラグメントBの原子との静電相互作用は、該原子と点電荷との間に共有結合を有すると仮定した際に、該原子間に2個以上の他の原子を含む場合に生じることとした。以下の他実施例においても、同様に記載することとする。
【化11】
【0219】
リガンドL
0及びリガンドL
1の構造から、リガンドL
0を構成する全ての原子を選択し、「構造a」とした。したがって、リガンドL
0は、構造aからなり、リガンドL
1は、構造aと該構造aに結合したフラグメントBとからなる。
【0220】
本実施例では、
図1に示される計算装置を用いて、
図2に示されるステップS110〜S190、S210〜S280、S300及びS400に対応するステップS110A〜S190A、S210A〜S280A、S300A及びS400Aを実行することにより、ΔG4a’を計算した。ステップS120Aでは、
図3に示されるステップS121〜S123に対応するステップS121A〜S123Aを実行し、ステップS130Aでは、
図5に示されるステップS131b〜S134bに対応するステップS131A〜S134Aを実行し、ステップS220Aでは、
図6に示されるステップS221〜S223に対応するステップS221A〜S223Aを実行した。
【0221】
[変化前の原子集合体Aに含まれる複合体L
0R及び変化後の原子集合体ABに含まれる複合体L
1Rに関するデータの入力]
ステップS110A及びS210Aを実行する前の準備ステップとして、入力部11は、変化前の原子集合体Aに含まれる複合体L
0R及び変化後の原子集合体ABに含まれる複合体L
1Rに関するデータ(複合体L
1Rを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等、複合対L
0Rを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)を制御部12に入力した。ここで、原子の種類とは、周期表における元素の種類及び原子タイプであって、各原子の持つ原子タイプ又は各原子間の原子タイプのペアは、ポテンシャルパラメータの値と1対1に対応して与えられるものである。入力されたデータは、制御部12により記憶部13に記憶された。
【0222】
変化後の原子集合体ABに含まれる複合体L
1Rに関するデータは、CCG(Chemical Computing Group)社より提供されている統合計算化学システムMOE(Molecular Operating Environment)を用いて、立体構造データベースであるProtein Data Bank(PDB;http://www.rcsb.org/pdb/home/home.do)よりダウンロードした立体構造データファイル(PDBコード:3PTB)に対して、該立体構造データファイルに含まれている、水分子に帰属される酸素原子を削除した後、MOEに実装されている「Protonate 3D」機能を用いて水素原子及び該水素原子の座標を追加することにより作成した立体構造データファイル(3PTB_AB)及び可視化ソフトVMD(Visual Molecular Dynamics:http://www.ks.uiuc.edu/Research/vmd/)に実装されている「Automatic PSF Builder」機能を用いて得た。なお、立体構造データ(PDBコード:3PTB)は、カルシウムに関する座標データを含んでおり、該カルシウムはイオンとして、また蛋白質の一部として取り扱った。また、変化前の原子集合体Aに含まれる複合体L
0Rに関するデータは、MOEの「Builder」機能を用いて、立体構造データファイル(3PTB_AB)よりフラグメントBを構成する仮想原子である点電荷以外の原子を除くことによって作成した立体構造データファイル(3PTB_A)に基づいて、VMDの「Automatic PSF Builder」機能を用いて作成した。
【0223】
なお、立体構造データファイル(3PTB_AB)には、複合体L
1Rに関する座標データが含まれており、立体構造データファイル(3PTB_AB)より、MOEを用いて、リガンドL
1に関するものを選択して、リガンドL
1に関する立体構造データファイル(3PTB_L
1)を作成し、VMDの「Automatic PSF Builder」機能により、リガンドL
1に関するデータを得ることができる。同様に、立体構造データファイル(3PTB_A)よりガンドL
0に関する立体構造データファイル(3PTB_L
0)を作成し、さらに、部分電荷に関しては、仮想原子を除いたフラグメントBを構成する原子(つまり、水素原子)と結合するベンゼン環上の炭素原子の部分電荷を零と設定することにより、リガンドL
0に関するデータを得た。
【0224】
[シミュレーション条件に関するデータの入力]
後述するステップS110A及びS210Aを実行する前の準備ステップにおいて、入力部11は、エネルギー極小化計算に関するシミュレーション条件(各原子に掛かる力の計算回数、ポテンシャルパラメータの種類及び値、境界条件、ファンデールワールスポテンシャル及びクーロンポテンシャル計算に関するスイッチング関数及びカットオフ半径、クーロンポテンシャル計算における長距離相互作用、1−4相互作用に関する条件等)、及び分子動力学シミュレーションに関するシミュレーション条件(シミュレーション時間、温度条件、圧力条件、ポテンシャルパラメータの種類及び値、発生させるアンサンブルの種類、境界条件、運動方程式の数値解法、数値積分の時間刻み、ファンデールワールスポテンシャル及びクーロンポテンシャル)計算に関するスイッチング関数及びカットオフ半径、クーロンポテンシャル計算における長距離相互作用に関する条件、1−4相互作用、スナップショットの数等の出力条件等)をそれぞれ制御部12に入力した。入力されたデータは、制御部12により記憶部13に記憶された。
【0225】
蛋白質Rを構成する各原子のポテンシャルパラメータの種類としては、CHARMm22を、リガンドL
0を構成する各原子のポテンシャルパラメータの種類及びリガンドL
1を構成する各原子のポテンシャルパラメータの種類としては、CGenFF(Charmm General ForceField)を、水分子を構成する各原子のポテンシャルパラメータの種類としては、TIP3Pを用いた。なお、リガンドL
0の各原子のポテンシャルパラメータ及びリガンドL
1の各原子のポテンシャルパラメータは、MATCHプログラム(http://brooks.chem.lF.umich.edu/index.php?page=match&subdir=articles/resources/soFtware)を用いて決定した。
【0226】
[ステップS110A:第1の原子集合体モデルの作成]
ステップS110Aにおいて、制御部12は、可視化ソフトVMDを用いて、記憶部13に記憶されている立体構造データファイル(3PTB_A)を読み込み、次いで、該VMDに実装されている「solvate」機能を用いて、タンパク質Rと該タンパク質Rに結合したリガンドL
0との複合体L
0Rの周囲に水分子を9492分子発生させた直方体型基本セルを作成し、変化前の原子集合体Aをモデル化した第1の原子集合体モデルを作成した。作成された第1の原子集合体モデルに関するデータ(原子集合体Aを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)は、制御部12により記憶部13に記憶された。ここで、原子の種類とは、周期表における元素の種類及びコンピュータシミュレーションに用いるポテンシャルパラメータを指定する原子タイプである。
【0227】
[ステップS120A:コンピュータシミュレーションによる原子集合体Aの座標の取得]
ステップS120Aにおいて、制御部12は、以下のステップS121A〜S123Aを実行した。
【0228】
[ステップS121A:第1の原子集合体モデルに対するエネルギー極小化計算]
ステップS121Aにおいて、制御部12は、第1の原子集合体モデルに対して、エネルギー極小化計算に関するシミュレーション条件(各原子に掛かる力の計算回数、ポテンシャルパラメータの種類及び値、境界条件、ファンデールワールスポテンシャル及びクーロンポテンシャル計算に関するスイッチング関数及びカットオフ半径、クーロンポテンシャル計算における長距離相互作用、1−4相互作用に関する条件等)を分子動力学シミュレーションプログラムであるNAMDプログラム(http://www.ks.uiuc.edu/Research/namd/)に読み取らせ、NAMDの「minimize」機能を用いたエネルギー極小化計算を実施した。エネルギー極小化計算後の第1の原子集合体モデルに関するデータ(原子集合体Aを構成する各原子の座標、種類、質量、部分電荷、原子タイプ、原子間の結合情報等)は、制御部12により記憶部13に記憶された。ここで、原子の種類とは、周期表における元素の種類及びコンピュータシミュレーションに用いるポテンシャルパラメータを指定する原子タイプである。
【0229】
エネルギー極小化計算の条件は、以下の通りである。
境界条件として周期境界条件を採用し、計算回数に関しては、NAMDプログラムの入力ファイルにおいて、「minimize」を10000とした。
スイッチング関数及びカットオフ半径に関しては、NAMDプログラムの入力ファイルにおいて、ファンデールワールスポテンシャル及びクーロンポテンシャル計算に用いるスイッチング関数に関するキーワード「switching」を10オングストローム、カットオフ半径に関するキーワード「cutoff」を12オングストロームと設定した。
クーロンポテンシャルの計算の長距離相互作用に関しては、「PMEGridSpacing」を1.0オングストローム及び「PMEInterpOrder」を4に設定したパーティクルメッシュエワルド(PME:Particle Mesh Ewald)法を用いた。
1−4相互作用に関しては、「exclude」をscaled1−4、「1−4scaling」を1とした。
【0230】
[ステップS122A:第1の原子集合体モデルの平衡化]
ステップS122Aにおいて、制御部12は、記憶部13に記憶されているエネルギー極小化計算後の第1の原子集合体モデルに関するデータ(原子集合体Aを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)、シミュレーション条件に関するデータ(シミュレーション時間、温度条件、圧力条件、ポテンシャルパラメータの種類及び値、発生させるアンサンブルの種類、境界条件、運動方程式の数値解法、数値積分の時間刻み、ファンデールワールスポテンシャル及びクーロンポテンシャル計算に関するスイッチング関数及びカットオフ半径、クーロンポテンシャル計算における長距離相互作用、1−4相互作用に関する条件、スナップショットの数等の出力条件等)を、NAMDプログラムに読み取らせ、エネルギー極小化計算後の第1の原子集合体モデルに対して分子動力学シミュレーションを実行し、第1の原子集合体モデルを平衡化させた。平衡化後の第1の原子集合体モデルに関するデータ(原子集合体Aを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)は、制御部12により記憶部13に記憶された。ここで、原子の種類とは、周期表における元素の種類及びコンピュータシミュレーションに用いるポテンシャルパラメータを指定する原子タイプである。
【0231】
なお、分子動力学シミュレーションの条件は、以下の通りに設定した。
境界条件として、周期境界条件を採用した。圧力条件に関しては、NAMDプログラムの入力ファイルにおいて、圧力制御のためのキーワード「langevinPiston」をonに設定し、「langevinPistonTarget」を1.01325と設定することで1atmに制御した。温度条件に関しては、NAMDプログラムの入力ファイルにおいて、温度制御のためのキーワード「langevin」をonに設定し、「langevinTemp」を50ケルビン〜300ケルビンまで50ケルビン毎に設定することにより段階的に昇温させた。シミュレーション時間は、50ケルビン〜250ケルビンではそれぞれ100ピコ秒とし、300ケルビンでは1000ピコ秒とした。運動方程式の数値解法としてshake法を採用して、時間刻みは2フェムト秒とした。出力条件として、200フェムト秒毎に第1の原子集合体モデルの座標を出力させた。
【0232】
スイッチング関数及びカットオフ半径に関しては、NAMDプログラムの入力ファイルにおいて、ファンデールワールスポテンシャル及びクーロンポテンシャル)計算に用いるスイッチング関数に関するキーワード「switching」を10オングストローム、カットオフ半径に関するキーワード「cutoff」を12オングストロームと設定した。
【0233】
クーロンポテンシャルの計算の長距離相互作用に関しては、「PMEGridSpacing」を1.0オングストローム及び「PMEInterpOrder」を4に設定したパーティクルメッシュエワルド(PME:Particle Mesh Ewald)法を用いた。
1−4相互作用に関しては、「exclude」をscaled1−4、「1−4scaling」を1とした。
【0234】
アンサンブルに関しては、NAMDプログラムに実装されている上記温度制御及び圧力制御を実行することによりNPTアンサンブルを構成した。
【0235】
[ステップS123A:コンピュータシミュレーションによる原子集合体Aの座標の取得]
ステップS123Aにおいて、制御部12は、記憶部13に記憶されている平衡化後の第1の原子集合体モデルに関するデータ(原子集合体Aを構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)、シミュレーション条件に関するデータ(シミュレーション時間、温度条件、圧力条件、ポテンシャルパラメータの種類及び値、発生させるアンサンブルの種類、境界条件、運動方程式の数値解法、数値積分の時間刻み、ファンデールワールスポテンシャル及びクーロンポテンシャル計算に関するスイッチング関数及びカットオフ半径、クーロンポテンシャル計算における長距離相互作用に関する条件、1−4相互作用スナップショットの数等の出力条件等)を、NAMDプログラムに読み取らせ、平衡化後の第1の原子集合体モデルに対して分子動力学シミュレーションを実行した。ここで、原子の種類とは、周期表における元素の種類及びコンピュータシミュレーションに用いるポテンシャルパラメータを指定する原子タイプである。
【0236】
なお、分子動力学シミュレーションの条件は、以下の通りに設定した。
境界条件として、周期境界条件を採用した。圧力条件に関しては、NAMDプログラムの入力ファイルにおいて、圧力制御のためのキーワード「langevinPiston」をonに設定し、「langevinPistonTarget」を1.01325と設定することで1atmに制御した。温度条件に関しては、NAMDプログラムの入力ファイルにおいて、温度制御のためのキーワード「langevin」をonに設定し、「langevinTemp」を300ケルビンに設定した。シミュレーション時間は、5000ピコ秒とし、運動方程式の数値解法としてshake法を採用して、時間刻みは2フェムト秒とした。出力条件に関しては、50フェムト秒毎に第1の原子集合体モデルの座標を出力させた。
【0237】
スイッチング関数及びカットオフ半径に関しては、NAMDプログラムの入力ファイルにおいて、ファンデールワールスポテンシャル及びクーロンポテンシャル計算に用いるスイッチング関数に関するキーワード「switching」を10オングストローム、カットオフ半径に関するキーワード「cutoff」を12オングストロームと設定した。
【0238】
クーロンポテンシャルの計算の長距離相互作用に関しては、「PMEGridSpacing」を1.0オングストローム及び「PMEInterpOrder」を4に設定したパーティクルメッシュエワルド(PME:Particle Mesh Ewald)法を用いた。
1−4相互作用に関しては、「exclude」をscaled1−4、「1−4scaling」を1とした。
【0239】
アンサンブルに関しては、NAMDプログラムに実装されている上記温度制御及び圧力制御を実行することによりNPTアンサンブルを構成した。
【0240】
制御部12は、分子動力学シミュレーションの開始から50フェムト秒ごと、すなわち、50フェムト秒後(時間T
1)、100フェムト秒後(時間T
2)、・・・・、5000ピコ秒後(時間T
100000)における原子集合体Aの座標を取得した。これにより、時間T
i(i=1,2,・・・,100000)における状態F
iにおける原子集合体Aの座標R
A(F
i)が取得された。取得された原子集合体Aの座標は、R
A(F
1)、R
A(F
2)、・・・、R
A(F
100000)の合計100000組の座標である。
【0241】
制御部12は、原子集合体Aの座標R
A(F
1)〜R
A(F
100000)を時間軸に沿って並べて記憶部13に記憶させた。すなわち、制御部12は、時間T
i(i=1,2,・・・,100000)における状態F
iにおける原子集合体Aの座標R
A(F
i)を、時間T
iと対応付けて記憶部13に記憶させた。
【0242】
[ステップS130A:原子集合体ABの座標の取得]
ステップS130Aにおいて、制御部12は、ステップS131A〜S134Aを実行した。
【0243】
[ステップS131A:第3の原子集合体モデルの作成]
原子集合体Cとして、構造aと該構造aに結合したフラグメントBとからなるリガンドL
1を選択した。原子集合体Cにおいて、リガンドL
1の周囲は真空とした。
ステップS131Aにおいて、制御部12は、記憶部13に記憶されている第3の原子集合体モデルに関するデータ(リガンドL
1を構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)、エネルギー極小化計算に関するシミュレーション条件(各原子に掛かる力の計算回数、ポテンシャルパラメータの種類及び値、境界条件、ファンデールワールスポテンシャル及びクーロンポテンシャル計算に関するスイッチング関数及びカットオフ半径、クーロンポテンシャル計算における長距離相互作用、1−4相互作用に関する条件等)に関するデータを、NAMDプログラムに読み取らせ、第3の原子集合体モデルに対してエネルギー極小化計算を実施した。エネルギー極小化計算後の第3の原子集合体モデルに関するデータ(リガンドL
1を構成する各原子の座標、種類、質量、部分電荷、原子間の結合情報等)は、制御部12により記憶部13に記憶された。ここで、原子の種類とは、周期表における元素の種類及びコンピュータシミュレーションに用いるポテンシャルパラメータを指定する原子タイプである。
【0244】
なお、エネルギー極小化計算の条件は、以下の通りである。
境界条件に関しては、ステップS131Aでは、真空中におけるリガンドL
1のエネルギー極小化計算であるため、境界条件及びクーロンポテンシャル計算の長距離相互作用に関する条件は特に指定しなかった。各原子に掛かる力の計算回数に関しては、NAMDプログラムの入力ファイルにおいて、「minimize」を1000とした。
スイッチング関数及びカットオフ半径に関しては、NAMDプログラムの入力ファイルにおいて、ファンデールワールスポテンシャル及びクーロンポテンシャル計算に用いるスイッチング関数に関するキーワード「switching」を10オングストローム、カットオフ半径に関するキーワード「cutoff」を12オングストロームと設定した。
1−4相互作用に関しては、「exclude」をscaled1−4、「1−4scaling」を1とした。
【0245】
[ステップS132A:コンピュータシミュレーションによる原子集合体Cの座標の取得]
ステップS132Aにおいて、制御部12は、エネルギー極小化計算後の第3の原子集合体モデルに対して分子動力学シミュレーションを実行した。
【0246】
分子動力学シミュレーション条件は、以下の通りに設定した。
温度条件に関しては、NAMDプログラムの入力ファイルにおいて、温度制御のためのキーワード「langevin」をonに設定し、「langevinTemp」を300ケルビンに設定した。シミュレーション時間は、10ナノ秒とし、運動方程式の数値解法としてshake法を採用し、時間刻みは2フェムト秒とした。出力条件に関して、100フェムト秒毎に第3の原子集合体モデルの座標を出力させ、原子集合体Cの座標を100000組生成せしめた。
【0247】
スイッチング関数及びカットオフ半径に関しては、NAMDプログラムの入力ファイルにおいて、ファンデールワールスポテンシャル及びクーロンポテンシャル計算に用いるスイッチング関数に関するキーワード「switching」を10オングストローム、カットオフ半径に関するキーワード「cutoff」を12オングストロームと設定した。
1−4相互作用に関しては、「exclude」をscaled1−4、「1−4scaling」を1とした。
【0248】
制御部12は、分子動力学シミュレーションの開始から100フェムト秒ごと、すなわち、100フェムト秒後(時間T
1)、200フェムト秒後(時間T
2)、・・・・、10ナノ秒後(時間T
100000)における原子集合体Cの座標を取得した。これにより、時間T
k(k=1,2,・・・,100000)に対応する状態H
kにおける原子集合体Cの座標R
C(H
k)が取得された。取得された原子集合体Cの座標は、R
C(H
1)、R
C(H
2)、・・・、R
C(H
100000)の合計100000組の座標である。
【0249】
制御部12は、原子集合体Cの座標R
C(H
1)〜R
C(H
100000)を時間軸に沿って並べて記憶部13に記憶させた。すなわち、制御部12は、時間T
k(k=1,2,・・・,100000)に対応する状態H
kにおける原子集合体Cの座標R
C(H
k)を、時間T
kと対応付けて記憶部13に記憶させた。
【0250】
[ステップS133A:状態F
xにおける原子集合体Aの座標に対する原子集合体Cの座標のフィッティング]
ステップS133Aにおいて、制御部12は、原子集合体Aの座標R
A(F
1)と原子集合体Cの座標R
C(H
1)とのペアQ(F
1,H
1)、原子集合体Aの座標R
A(F
2)と原子集合体Cの座標R
C(H
2)とのペアQ(F
2,H
2)、・・・・、原子集合体Aの座標R
A(F
100000)と原子集合体Cの座標R
C(H
100000)とのペアQ(F
100000,H
100000)の合計100000組のペアを構成した。
【0251】
各ペアにおいて、構造a(リガンドL
0)を構成する全ての原子を選択し、「選択原子群」とした。
【0252】
制御部12は、ペアQ(F
1,H
1)〜Q(F
100000,H
100000)の各ペアについて、原子集合体Cの座標の各原子間の相対座標(内部座標)を保持しつつ、数式(7)を用いて、原子集合体Aのうち選択原子群を構成する原子に対して、原子集合体Cのうち選択原子群を構成する原子をフィッティングさせ、原子集合体Aに原子集合体Cを重ね合わせた。この際、数式(7)において各原子に対して設定されたCnは、以下の通りである。
【0253】
【化12】
【0254】
[ステップS134A:原子集合体ABの座標の取得]
ステップS134Aにおいて、制御部12は、ペアQ(F
1,H
1)〜Q(F
100000,H
100000)の各ペアについて、原子集合体Aに重ね合わせた原子集合体Cの座標から、原子集合体Cが有するフラグメントBの座標(ただし、仮想原子である点電荷の座標は除く)を取得し、これを原子集合体Aの座標に追加して、原子集合体ABの座標R
AB(F
i)を取得した。これにより、原子集合体ABの座標R
AB(F
1)〜R
AB(F
100000)が取得された。
【0255】
[ステップS140A:相互作用エネルギーφの計算]
ステップS140Aにおいて、制御部12は、原子集合体ABの座標R
AB(F
1)〜R
AB(F
100000)のそれぞれに関して、構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφを計算した。この際、制御部12は、以下のステップを実行した。
【0256】
制御部12は、原子集合体ABの座標R
AB(F
i)から、構造a及び該構造aに結合したフラグメントBからなる構造aB(リガンドL
1)を構成する各原子の座標R
L1(F
i)を抽出する処理をi=1,2,・・・,100000で繰り返して実行し、原子集合体ABの座標R
AB(F
1)〜R
AB(F
100000)から抽出された、構造aB(リガンドL
1)を構成する各原子の座標R
L1(F
1)〜R
L1(F
100000)を取得した。
【0257】
制御部12は、構造aB(リガンドL
1)を構成する各原子の座標、種類、質量、部分電荷、ポテンシャルパラメータの種類及び値、原子間の結合情報等をNAMDに読み取らせ、座標R
L1(F
i)における構造aB(リガンドL
1)のエネルギーE
L1(F
i)を計算する処理をi=1,2,・・・,100000で繰り返して実行し、座標R
L1(F
1)〜R
L1(F
100000)のそれぞれに関して、構造aB(リガンドL
1)のエネルギーE
L1(F
1)〜E
L1(F
100000)を計算した。ここで、E
L1(F
i)は、構造aB(リガンドL
1)の原子間の相互作用に起因するエネルギーであり、リガンドL
1の内部エネルギーである。
【0258】
制御部12は、原子集合体ABの座標R
AB(F
i)から、構造a(リガンドL
0)を構成する各原子の座標R
a(F
i)を抽出する処理をi=1,2,・・・,100000で繰り返して実行し、原子集合体ABの座標R
AB(F
1)〜R
AB(F
100000)から抽出された、構造aを構成する各原子の座標R
a(F
1)〜R
a(F
100000)を取得した。
【0259】
制御部12は、構造aを構成する各原子に関する座標、種類、質量、部分電荷、ポテンシャルパラメータの種類及び値、原子間の結合情報等をNAMDに読み取らせ、座標R
a(F
i)における構造aのエネルギーE
a(F
i)を計算する処理をi=1,2,・・・,100000で繰り返して実行し、座標R
a(F
1)〜R
a(F
100000)のそれぞれに関して、構造aのエネルギーE
a(F
1)〜E
a(F
100000)を計算した。ここで、E
a(F
i)は、構造aの原子間の相互作用に起因するエネルギーであり、構造aの内部エネルギーである。
【0260】
制御部12は、原子集合体ABの座標R
AB(F
i)から、仮想原子を除いたフラグメントBを構成する原子(構造aBがリガンドL
1である本実施例の場合、水素原子)の座標R
d(F
i)を取得する処理をi=1,2,・・・,100000で繰り返して実行し、原子集合体ABの座標R
AB(F
1)〜R
AB(F
100000)から抽出された、仮想原子を除いたフラグメントBを構成する原子の座標R
d(F
1)〜R
d(F
100000)を取得した。
【0261】
制御部12は、仮想原子を除いたフラグメントBを構成する原子に関する座標、種類、質量、部分電荷、ポテンシャルパラメータの種類及び値、原子間の結合情報等をNAMDに読み取らせ、座標R
d(F
i)に関し、仮想原子を除いたフラグメントBを構成する原子のエネルギーE
d(F
i)及びクーロンポテンシャルEC
d(F
i)を計算する処理を、i=1,2,・・・,100000で繰り返して実行し、座標R
d(F
1)〜R
d(F
100000)のそれぞれに関して、仮想原子を除いたフラグメントBを構成する原子のエネルギーE
d(F
1)〜E
d(F
100000)及びクーロンポテンシャルEC
d(F
1)〜EC
d(F
100000)を取得した。ここで、E
d(F
i)は、仮想原子を除いたフラグメントBを構成する原子間の相互作用に起因するエネルギーであり、EC
d(F
i)は仮想原子を除いたフラグメントBを構成する原子間の静電相互作用に起因するエネルギーである。なお、構造aBがリガンドL
1の場合、仮想原子を除いたフラグメントBを構成する原子は水素原子のみの1原子であるので、いずれの値も零であった。
【0262】
制御部12は、座標R
d(F
i)に関し、仮想原子を除いたフラグメントBを構成する原子のエネルギーE
d(F
i)から、クーロンポテンシャルEC
d(F
i)を引いて、エネルギーE
e(F
i)を計算する処理を、i=1,2,・・・,100000で繰り返して実行し、エネルギーE
e(F
1)〜E
e(F
100000)を計算した。
【0263】
制御部12は、原子集合体ABの座標R
AB(F
i)から、仮想原子を除いたフラグメントBを構成する原子と、構造aB(リガンドL
1)を構成する原子のうち、フラグメントBを構成する点電荷の部分電荷を有する原子(リガンドL
1のベンゼン環を構成する炭素原子のうち、フラグメントBを構成する水素原子と結合を有する炭素原子)とからなる構造fを構成する各原子(つまり、仮想原子を含めたフラグメントBを構成する原子)の座標R
f(F
i)を取得する処理をi=1,2,・・・,100000で繰り返して実行し、原子集合体ABの座標R
AB(F
1)〜R
AB(F
100000)から抽出された、構造fを構成する各原子の座標R
f(F
1)〜R
f(F
100000)を取得した。
【0264】
制御部12は、構造fを構成する各原子の座標、種類、質量、部分電荷、ポテンシャルパラメータの種類及び値、原子間の結合情報等をNAMDに読み取らせ、座標R
f(F
i)に関し、クーロンポテンシャルEC
f(F
i)を計算する処理を、i=1,2,・・・,100000で繰り返して実行し、クーロンポテンシャルEC
f(F
1)〜EC
f(F
100000)を取得した。ここで、EC
f(F
i)は、構造fを構成する原子間の静電相互作用に起因するエネルギーである。なお、構造aBがリガンドL
1である本実施例では、フラグメントBを構成する水素原子と点電荷の間に共有結合を有すると仮定した際に、該原子間に他原子を含まないため、EC
f(F
i)は零であった。
【0265】
制御部12は、エネルギーE
e(F
i)と、クーロンポテンシャルEC
f(F
i)とを足して、エネルギーE
B(F
i)を計算する処理を、i=1,2,・・・,100000で繰り返して実行し、エネルギーE
B(F
1)〜E
B(F
100000)を計算した。
【0266】
制御部12は、次式:
φ(F
i)=E
L1(F
i)−E
a(F
i)−E
B(F
i) 数式(8)
に基づいて、構造aとフラグメントBとの相互作用エネルギーφを計算する処理を、i=1,2,・・・,100000で繰り返して実行した。
【0267】
[ステップS150A:相互作用エネルギーφの各階級の度数分布の作成]
ステップS150Aにおいて、制御部12は、横軸が相互作用エネルギーφの各階級を表し、縦軸が相互作用エネルギーφの各階級の度数を表すヒストグラムを作成した。なお、相互作用エネルギーφに関するヒストグラムにおいて、階級間隔Δφは、各階級間隔Δφにおける該Δφとその度数H
0(φ)の積が一定となるように、各階級間隔Δφを設定し、相互作用エネルギーφの分割数、つまり階級間隔Δφの個数は、100とした。
【0268】
[ステップS160A:相互作用エネルギーφの各階級の出現確率P
0(φ)の計算]
ステップS160Aにおいて、制御部12は、ステップS150Aで作成されたヒストグラムを規格化し、相互作用エネルギーφの各階級の出現確率P
0(φ)を計算した。
【0269】
[ステップS170A:相互作用エネルギーεの計算]
入力部11は、エネルギー表示法プログラム群(http://sourceForge.net/p/ermod/wiki/Home/)の一部であるgen_structure、gen_inputプログラムを用いて作成された、パラメータファイル等の、計算に必要な各種ファイルを、制御部12に入力した。入力されたファイルは、制御部12により記憶部13に記憶させた。なお、本発明において用いたエネルギー表示法プログラム群は、バージョン0.2.3のものを用いた。
【0270】
なお、gen_inputプログラムにより作成された、構造aB(リガンドL
1)を構成する各原子の部分電荷及びファンデールワールスに関するパラメータが入力されているSltInfoファイルに関して、構造aを構成する原子のファンデールワールスに関するパラメータは零に変更した。また、仮想原子を含めたフラグメントBを構成する原子の部分電荷以外は零に変更した。
【0271】
制御部12は、エネルギー表示法プログラム群の一部であるermodプログラムを用いて、原子集合体ABの座標R
AB(F
i)に関して、フラグメントBと蛋白質R(本実施例では1個)との相互作用エネルギーεp(i)、及び、フラグメントBと各水分子(本実施例では9492個の水分子のそれぞれ)との相互作用エネルギーεw(i)を計算する処理を、i=1,2,・・・,100000で繰り返し実行した。ここで、各εw(i)は、フラグメントBと各水分子との合計9492個の相互作用エネルギーの値を持つ集合である。
【0272】
[ステップS180A:相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数分布の作成]
制御部12は、ステップS170Aで計算されたεp(1)〜εp(100000)及びεw(1)〜εw(100000)から、ステップS150Aで作成されたヒストグラムにおける相互作用エネルギーφの各階級に属する座標を用いて計算されたεp及びεwを選択し、選択されたεp及びεwに関するヒストグラムH
0'(εp;φ)、H
0'(εw;φ)を作成した。
【0273】
[ステップS190A:相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の出現確率P
0'(ε;φ)の計算]
制御部12は、ヒストグラムH
0'(εp;φ)及びH
0'(εw;φ)における各階級の出現確率P
0'(εp;φ)及びP
0'(εw;φ)を計算した。
ここで、ステップS180A及びS190Aにおける、ヒストグラムH
0'(εp;φ)、H
0'(εw;φ)及び出現確率P
0'(εp;φ)、P
0'(εw;φ)は、エネルギー表示法プログラム群の一部であるermodプログラムにより作成した。なお、ermodプログラムのパラメータに関して、ヒストグラムH
0'(εp;φ)及び出現確率P
0'(εp;φ)の作成では、ermodプログラムのパラメータを入力するparameter_erファイル中において、デフォルト条件よりecfbin及びec0binの設定を0.05、engdivの設定を10、maxinsの設定を1に変更した。また、H
0'(εw;φ)及びP
0'(εw;φ)の作成は、デフォルト条件よりengdivの設定を10、maxinsの設定を1に変更して実施した。
【0274】
制御部12は、ステップS180A〜S190Aを、ステップS150Aにて作成されたヒストグラムにおける相互作用エネルギーφの各階級に対して実施し、相互作用エネルギーφの各階級に関して、相互作用エネルギーεの出現確率P
0'(εp;φ)及びP
0'(εw;φ)を計算した。
【0275】
[ステップS210A:第2の原子集合体モデルの作成]
[ステップS220A:コンピュータシミュレーションによる原子集合体ABの座標の取得]
ステップS210A及びS220Aにおいて、制御部12は、原子集合体Aを原子集合体ABに変更した点を除き、ステップS110A及びS120Aと同様にしてステップS210A及びS220Aを実行し、ステップS223Aの分子動力学シミュレーションの開始から50フェムト秒ごと、すなわち、50フェムト秒後(時間T
1)、100フェムト秒後(時間T
2)、・・・・、5000ピコ秒後(時間T
100000)における原子集合体ABの座標を取得した。これにより、時間T
j(j=1,2,・・・,100000)における状態G
jにおける原子集合体ABの座標R
AB(G
j)が取得された。取得された原子集合体ABの座標は、R
AB(G
1)、R
AB(G
2)、・・・、R
AB(G
100000)の合計100000組の座標である。
【0276】
制御部12は、原子集合体ABの座標R
AB(G
1)〜R
AB(G
100000)を時間軸に沿って並べて記憶部13に記憶させた。すなわち、制御部12は、時間T
j(j=1,2,・・・,100000)における状態G
jにおける原子集合体ABの座標R
AB(G
j)を、時間T
jと対応付けて記憶部13に記憶させた。
【0277】
[ステップS230A:相互作用エネルギーφの計算]
ステップS230Aにおいて、制御部12は、ステップS140Aと同様にしてステップS230Aを実行し、原子集合体ABの座標R
AB(G
1)〜R
AB(G
100000)のそれぞれに関して、構造aと該構造aに結合したフラグメントBとの相互作用エネルギーφを計算した。
【0278】
[ステップS240A:相互作用エネルギーφの各階級の度数分布の作成]
ステップS240Aにおいて、制御部12は、横軸が相互作用エネルギーφの各階級を表し、縦軸が相互作用エネルギーφの各階級の度数を表すヒストグラムを作成した。各階級間隔Δφは、ステップS150Aで設定した各Δφと同じΔφとし、相互作用エネルギーφの分割数、つまり階級間隔Δφの個数は、100とした。
【0279】
[ステップS250A:相互作用エネルギーφの各階級の出現確率P(φ)の計算]
ステップS250Aにおいて、制御部12は、ステップS240Aで作成されたヒストグラムを規格化し、相互作用エネルギーφの各階級の出現確率P(φ)を計算した。
【0280】
[ステップS260A:相互作用エネルギーεの計算]
ステップS260Aにおいて、制御部12は、エネルギー表示法プログラム群の一部であるermodプログラムを用いて、原子集合体ABの座標R
AB(G
j)に関して、フラグメントBと蛋白質R(本実施例では1個)との相互作用エネルギーεp(i)、及び、フラグメントBと各水分子(本実施例では9492個の水分子のそれぞれ)との相互作用エネルギーεw(i)を計算する処理を、i=1,2,・・・,100000で繰り返し実行した。ここで、各εw(i)は、フラグメントBと各水分子との合計9492個の相互作用エネルギーの値を持つ集合である。この際、ステップS170Aと同様にして、計算に必要な各種ファイルを作成するとともに、SltInfoファイルを変更した。
【0281】
[ステップS270A:相互作用エネルギーεの各階級の度数分布の作成]
ステップS270Aにおいて、制御部12は、ステップS260Aで計算されたεp(1)〜εp(100000)及びεw(1)〜εw(100000)に関するヒストグラムH'(εp)、H'(εw)を作成した。
【0282】
[ステップS280A:相互作用エネルギーεの各階級の出現確率P'(ε)の計算]
ステップS280Aにおいて、制御部12は、ヒストグラムH'(εp)及びH'(εw)における各階級の出現確率P'(εp)及びP'(εw)を計算した。
ここで、ステップS270A及びS280Aにおける、ヒストグラムH'(εp)、H'(εw)及び出現確率P'(εp)、P'(εw)は、エネルギー表示法プログラム群の一部であるermodプログラムにより作成した。ermodプログラムのパラメータに関して、ヒストグラムH'(εp)及びP'(εp)の作成では、ermodプログラムのパラメータを入力するparameter_erファイル中において、デフォルト条件よりecfbin及びec0binの設定を0.05と変更した。また、H'(εw)及びP'(εw)の作成は、デフォルト条件により実施した。
【0283】
制御部12は、ステップS260A〜S280Aにより、相互作用エネルギーεの各階級の出現確率P'(εp)及びP'(εw)を計算した。
【0284】
[ステップS300A:相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφの計算]
ステップS300Aにおいて、制御部12は、以下のステップにより、相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφを計算した。
制御部12は、ステップS250Aで計算された、相互作用エネルギーφの各階級の出現確率P(φ)と、ステップS190Aで計算された、相互作用エネルギーφの各階級における相互作用エネルギーεpの各階級の出現確率P
0'(εp:φ)と、ステップS280Aで計算された、相互作用エネルギーεpの各階級の出現確率P'(εp)と、エネルギー表示法プログラム群の一部であるslvfeプログラムとを用いて、∫Δνp(φ)P(φ)dφを計算した。
【0285】
制御部12は、ステップS250Aで計算された、相互作用エネルギーφの各階級の出現確率P(φ)と、ステップS190Aで計算された、相互作用エネルギーφの各階級における相互作用エネルギーεwの各階級の出現確率P
0'(εw:φ)と、ステップS280Aで計算された、相互作用エネルギーεwの各階級の出現確率P'(εw)と、エネルギー表示法プログラム群の一部であるslvfeプログラムとを用いて、∫Δνw(φ)P(φ)dφを計算した。
【0286】
制御部12は、slvfeプログラムの出力として得られるフラグメントBの自己相互作用エネルギーEself(Self−energy)と、上記で計算された∫Δνp(φ)P(φ)dφ及び∫Δνw(φ)P(φ)dφとを用いて、数式(9)により、相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφを計算した。なお、Eselfは、エネルギー表示法において、周期境界条件を課した分子動力学シミュレーションを実施した際に、補正項として算出される値である。
∫Δν(φ)P(φ)dφ=∫Δνp(φ)P(φ)dφ+∫Δνw(φ)P(φ)dφ+Eself 数式(9)
【0287】
[ステップS400A:ΔG4aの計算]
ステップS400Aにおいて、制御部12は、ステップS160Aで計算された、相互作用エネルギーφの各階級の出現確率P
0(φ)と、ステップS250Aで計算された、相互作用エネルギーφの各階級の出現確率P(φ)と、ステップS300Aで計算された∫Δν(φ)P(φ)dφと、数式(1):
【数12】
[式中、Rは気体定数を表し、Tは分子動力学シミュレーションで設定した絶対温度(300ケルビン)を表す。]
とに基づいて、ΔG4a’を計算し、計算されたΔG4a’にマイナスを乗じてΔG4aを計算した。
【0288】
なお、数式(1)の第1項では、相互作用エネルギーφを単純平均することによって算出した。
【0289】
計算されたΔG4aは、4.18(kcal/mol)であった。
【0290】
(2)ΔG3aの計算
反応式(9):
【化13】
で表される変化に関し、リガンドL
0及び該リガンドL
0の周囲に位置する水分子を構成する変化前の原子集合体を「原子集合体A」とし、リガンドL
1及び該リガンドL
1の周囲に位置する水分子を構成する変化後の原子集合体を「原子集合体AB」とし、変化前の原子集合体Aの自由エネルギー及びフラグメントBの自由エネルギーの和と、変化後の原子集合体ABの自由エネルギーとの差ΔG3a’を数式(1)に基づいて計算した。
【0291】
なお、ΔG3aは、反応式(10):
【化14】
で表される変化に関する自由エネルギー変化量であるので、数式(1)に基づいて計算されるΔG3a’にマイナスを乗じたものがΔG3aに相当する。
【0292】
本実施例では、
図1に示される計算装置を用いて、
図2に示されるステップS110〜S190、S210〜S280、S300及びS400に対応するステップS110B〜S190B、S210B〜S280B、S300B及びS400Bを実行することにより、ΔG3a’を計算した。ステップS120Bでは、
図3に示されるステップS121〜S123に対応するステップS121B〜S123Bを実行し、ステップS130Bでは、
図5に示されるステップS131b〜S134bに対応するステップS131B〜S134Bを実行し、ステップS220Bでは、
図6に示されるステップS221〜S223に対応するステップS221B〜S223Bを実行した。
【0293】
[変化前の原子集合体Aに含まれるリガンドL
0及び変化後の原子集合体ABに含まれるリガンドL
1に関するデータの入力]
ステップS110B及びS210Bを実行する前の準備ステップにおいて、入力部11は、変化前の原子集合体Aに含まれるリガンドL
0及び変化後の原子集合体ABに含まれるリガンドL
1に関するデータを制御部12に入力した。入力されたデータは、制御部12により記憶部13に記憶された。入力されたデータは、ΔG4aの計算で用いた複合体L
0R及びL
1Rに関するデータより、リガンドL
0及びリガンドL
1に関するデータを抽出したものであって、座標データは、それぞれ立体構造ファイル3PTB_L
0、3PTB_L
1に存在する。
【0294】
[シミュレーション条件に関するデータの入力]
ステップS110B及びS210Bを実行する前の準備ステップにおいて、入力部11は、シミュレーション条件に関するデータを制御部12に入力した。入力されたデータは、制御部12により記憶部13に記憶された。入力されたデータは、蛋白質Rに関するものを除いた以外は、ΔG4aの計算と同様である。
【0295】
[ステップS110B:第1の原子集合体モデルの作成]
ステップS110Bにおいて、制御部12は、可視化ソフトVMDを用いて、記憶部13に記憶されているPDBファイル(3PTB_L
0)を読み込み、次いで、該VMDに実装されている「solvate」機能を用いて、リガンドL
0を構成する原子の重心座標を原点として、815分子の水分子を発生させて、直方体型の基本セルを作成し、変化前の原子集合体Aをモデル化した第1の原子集合体モデルを作成した。なお、各原子のポテンシャルパラメータに関して、リガンドL
0及びリガンドL
1を構成する各原子にはCGenFF、水分子を構成する各原子にはTIP3Pを用いた。また、リガンドL
0及びリガンドL
1の各原子のポテンシャルパラメータは、MATCHプログラムを用いて決定した。
【0296】
[ステップS120B:コンピュータシミュレーションによる原子集合体Aの原子の取得]
[ステップS130B:原子集合体ABの座標の取得]
[ステップS140B:相互作用エネルギーφの計算]
[ステップS150B:相互作用エネルギーφの各階級の度数分布の作成]
[ステップS160B:相互作用エネルギーφの各階級の出現確率P
0(φ)の計算]
[ステップS170B:相互作用エネルギーεの計算]
[ステップS180B:相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の度数分布の作成]
[ステップS190B:相互作用エネルギーφの各階級における、相互作用エネルギーεの各階級の出現確率P
0'(ε;φ)の計算]
制御部12は、ステップS120A〜S190Aと同様にしてステップS120B〜S190Bを実行した。
なお、分子動力学シミュレーションの条件は、段階昇温による平衡化を実施せず、エネルギー極小化後に温度300K/シミュレーション時間100ピコ秒の平衡化シミュレーションを実施したこと、及び平衡化後シミュレーション時間を1000ピコ秒とし、10フェムト秒刻みに原子集合体Aの座標を出力した以外はステップS120Aと同様にした。
【0297】
[ステップS210B:第2の原子集合体モデルの作成]
[ステップS220B:コンピュータシミュレーションによる原子集合体ABの座標の取得]
[ステップS230B:相互作用エネルギーφの計算]
[ステップS240B:相互作用エネルギーφの各階級の度数分布の作成]
[ステップS250B:相互作用エネルギーφの各階級の出現確率P(φ)の計算]
[ステップS260B:相互作用エネルギーεの計算]
[ステップS270B:相互作用エネルギーεの各階級の度数分布の作成]
[ステップS280B:相互作用エネルギーεの各階級の出現確率P'(ε)の計算]
制御部12は、ステップS210A〜S280Aと同様にしてステップS210B〜S280Bを実行した。
なお、分子動力学シミュレーションの条件は、段階昇温による平衡化を実施せず、エネルギー極小化後に温度300K/シミュレーション時間100ピコ秒の平衡化シミュレーションを実施したこと、及び平衡化後シミュレーション時間を1000ピコ秒とし、10フェムト秒刻みに原子集合体ABの座標を出力した以外はステップS220Aと同様にした。
【0298】
[ステップS300B:相互作用エネルギーεに起因する自由エネルギー変化量∫Δν(φ)P(φ)dφの計算]
[ステップS400B:ΔG3aの計算]
制御部12は、ステップS300A〜S400Aと同様にしてステップS300B〜S400Bを実行した。
【0299】
計算されたΔG3aは、1.59(kcal/mol)であった。
【0300】
(3)ΔG4bの計算
反応式(11):
【化15】
で表される変化に関し、複合体L
0R及び該複合体L
0Rの周囲に位置する水分子を構成する変化前の原子集合体を「原子集合体A」とし、複合体L
2R及び該複合体L
2Rの周囲に位置する水分子を構成する変化後の原子集合体を「原子集合体AB」とし、変化前の原子集合体Aの自由エネルギー及びフラグメントBの自由エネルギーの和と、変化後の原子集合体ABの自由エネルギーとの差ΔG4bを計算した。
【0301】
本実施例では、
図1に示される計算装置を用いて、
図2に示されるステップS110〜S190、S210〜S280、S300及びS400に対応するステップS110C〜S190C、S210C〜S280C、S300C及びS400Cを実行することにより、ΔG4bを計算した。ステップS120Cでは、
図3に示されるステップS121〜S123に対応するステップS121C〜S123Cを実行し、ステップS130Cでは、
図5に示されるステップS131b〜S134bに対応するステップS131C〜S134Cを実行し、ステップS220Cでは、
図6に示されるステップS221〜S223に対応するステップS221C〜S223Cを実行した。
【0302】
リガンドL
2の化学構造は、以下の通りである。
【化16】
【0303】
リガンドL
2を構成する各原子の部分電荷は、以下の通りである。
【化17】
【0304】
フラグメントBの構造は以下の通りである。
【化18】
【0305】
ステップS110C〜S190C、S210C〜S280C、S300C及びS400Cは、リガンドL
1をリガンドL
2に変更した点及びフラグメントBを変更した点を除き、ΔG4aの計算におけるステップS110A〜S190A、S210A〜S280A、S300A及びS400Aと同様に実行した。計算されたΔG4bは、−17.84(kcal/mol)であった。
【0306】
(4)ΔG3bの計算
反応式(12):
【化19】
で表される変化に関し、リガンドL
0及び該リガンドL
0の周囲に位置する水分子を構成する変化前の原子集合体を「原子集合体A」とし、リガンドL
2及び該リガンドL
2の周囲に位置する水分子を構成する変化後の原子集合体を「原子集合体AB」とし、変化前の原子集合体Aの自由エネルギー及びフラグメントBの自由エネルギーの和と、変化後の原子集合体ABの自由エネルギーとの差ΔG3bを計算した。
【0307】
本実施例では、
図1に示される計算装置を用いて、
図2に示されるステップS110〜S190、S210〜S280、S300及びS400に対応するステップS110D〜S190D、S210D〜S280D、S300D及びS400Dを実行することにより、ΔG3bを計算した。ステップS120Dでは、
図3に示されるステップS121〜S123に対応するステップS121D〜S123Dを実行し、ステップS130Dでは、
図5に示されるステップS131b〜S134bに対応するステップS131D〜S134Dを実行し、ステップS220Dでは、
図6に示されるステップS221〜S223に対応するステップS221D〜S223Dを実行した。
【0308】
ステップS110D〜S190D、S210D〜S280D、S300D及びS400Dは、リガンドL
1をリガンドL
2に変更した点及びフラグメントBを変更した点を除き、ΔG3aの計算におけるステップS110B〜S190B、S210B〜S280B、S300B及びS400Bと同様に実行した。計算されたΔG3bは、−14.91(kcal/mol)であった。
【0309】
(5)ΔΔGの計算
ΔΔG=ΔG4−ΔG3
ΔG4=ΔG4a+ΔG4b
ΔG3=ΔG3a+ΔG3b
より計算されたΔΔGは、−0.34(kcal/mol)であった。
【0310】
公知文献(Journal of Chemical Theory and Computation,8,3686−3695(2012)及びJouranl of the American Chemical Society,99,2331−2336(1977))に記載されている、実験により求められたΔΔGは、−0.10(kcal/mol)であり、本発明の計算方法によって計算されたΔΔGは、実験により求められたΔΔGとの差が1kcal/mol以内であって良く一致した。
【0311】
〔実施例2〕
リガンドL
2の代わりに、リガンドL
3を用いて実施例1と同様にしてΔΔGを計算した。
リガンドL
3の化学構造は、以下の通りである。
【化20】
【0312】
リガンドL
3を構成する各原子の部分電荷は、以下の通りである。
【化21】
【0313】
フラグメントBの構造は、以下の通りである。
【化22】
【0314】
ここで、以下の通り、フラグメントBの原子として点電荷(部分電荷;0.000)を便宜上取り扱うことによって、実施例1と同様にΔΔGを計算した。
【化23】
【0315】
ΔG4a及びΔG3aは、実施例1と同一の値を用いた。
計算されたΔG4bは、−8.12(kcal/mol)であった。
計算されたΔG3bは、−5.45(kcal/mol)であった。
計算されたΔΔGは、−0.08(kcal/mol)であった。
公知文献(Journal of Chemical Theory and Computation,8,3686−3695(2012)及びJouranl of the American Chemical Society,99,2331−2336(1977))に記載されている、実験により求められたΔΔGは、0.27(kcal/mol)であり、本発明の計算方法によって算出したΔΔGは、実験により求められたΔΔGとの差が1kcal/mol以内であって良く一致した。
【0316】
〔実施例3〕
リガンドL
2の代わりに、リガンドL
4を用いて実施例1と同様にしてΔΔGを計算した。
リガンドL
4の化学構造は、以下の通りである。
【化24】
【0317】
リガンドL
4を構成する各原子の部分電荷は、以下の通りである。
【化25】
【0318】
フラグメントBの構造は、以下の通りである。
【化26】
【0319】
ΔG4a及びΔG3aは、実施例1と同一の値を用いた。
計算されたΔG4bは、−23.00(kcal/mol)であった。
計算されたΔG3bは、−19.88(kcal/mol)であった。
計算されたΔΔGは、−0.53(kcal/mol)であった。
【0320】
公知文献(Journal of Chemical Theory and Computation,8,3686−3695(2012)及びJouranl of the American Chemical Society,99,2331−2336(1977))に記載されている、実験により求められたΔΔGは、−0.40(kcal/mol)であり、本発明の計算方法によって算出したΔΔGは、実験により求められたΔΔGとの差が1kcal/mol以内であって良く一致した。
【0321】
〔実施例4〕
リガンドL
2の代わりに、リガンドL
5を用いて実施例1と同様にしてΔΔGを計算した。
リガンドL
5の化学構造は、以下の通りである。
【化27】
【0322】
リガンドL
5を構成する各原子の部分電荷は、以下の通りである。
【化28】
【0323】
フラグメントBの構造は、以下の通りである。
【化29】
【0324】
ここで、以下の通り、フラグメントBの原子とて点電荷(部分電荷;0.000)を便宜上取り扱うことによって、実施例1と同様にΔΔGを計算した。
【化30】
【0325】
ΔG4a及びΔG3aは、実施例1と同一の値を用いた。
計算されたΔG4bは、1.78(kcal/mol)であった。
計算されたΔG3bは、3.55(kcal/mol)であった。
計算されたΔΔGは、0.82(kcal/mol)であった。
【0326】
公知文献(Jouranl of the American Chemical Society, 125,10570−10579(2003))に記載されている、実験により求められたΔΔGは、0.93(kcal/mol)であり、本発明の計算方法によって算出したΔΔGは、実験により求められたΔΔGとの差が1kcal/mol以内であって良く一致した。
【0327】
〔実施例5〕
リガンドL
2の代わりに、リガンドL
6を用いて実施例1と同様にしてΔΔGを計算した。
リガンドL
6の化学構造は、以下の通りである。
【化31】
【0328】
リガンドL
6を構成する各原子の部分電荷は、以下の通りである。
【化32】
【0329】
フラグメントBの構造は、以下の通りである。
【化33】
【0330】
ここで、以下の通り、フラグメントBの原子とて点電荷(部分電荷;0.000)を便宜上取り扱うことによって、実施例1と同様にΔΔGを計算した。
【化34】
【0331】
ΔG4a及びΔG3aは、実施例1と同一の値を用いた。
計算されたΔG4bは、−3.79(kcal/mol)であった。
計算されたΔG3bは、−0.95(kcal/mol)であった。
計算されたΔΔGは、−0.25(kcal/mol)であった。
【0332】
公知文献(Jouranl of the American Chemical Society,125,10570−10579(2003))に記載されている、実験により求められたΔΔGは、0.25(kcal/mol)であり、本発明の計算方法によって算出したΔΔGは、実験により求められたΔΔGとの差が1kcal/mol以内であって良く一致した。