JP6441301B2 - 相対結合親和性および誤差のサイクルクロージャ推定 - Google Patents

Description

本発明は、バイオインフォマティクスの一般的な分野におけるものであり、より詳細には、受容体／リガンドの相対結合親和性およびそれらの誤差の推定または計算におけるものである。

生物学的プロセスは多くの場合、タンパク質−リガンド結合事象（例えば、リガンドとこのリガンドの受容体との結合）に依拠し、このため、関連するエネルギー論の正確な計算が計算構造に基づく薬剤設計の中心となる目標である^１〜３。

速度と精度との様々なトレードオフを有するタンパク質−リガンド結合自由エネルギーを計算する多岐にわたる様々な手法が開発されている。これらの手法は、仮想スクリーニングおよび分子力学法／一般化ボルン表面積（ＭＭ／ＧＢＳＡ）または分子力学法／ポアソン−ボルツマン（ＭＭ／ＰＢＳＡ）モデル等の高速エンドポイント法、ＧＢ法よりもタンパク質の脱溶媒和からの寄与を明示的に考慮するＷＭ／ＭＭ法、ならびに基礎をなす力場の精度限界内での結合の完全なまたは本質的に完全なエネルギー記述および完全なサンプリングを考慮するように設計される自由エネルギー摂動（ＦＥＰ）法および熱力学的積分（ＴＩ）法を含む^{１，２，４〜１０}。

タンパク質−リガンド結合親和性を計算する様々な方法の中でも、ＴＩ、ＦＥＰ、ラムダ動力学、アルケミカルＯＳＲＷ（ａｌｃｈｅｍｉｃａｌ ＯＳＲＷ）等の自由エネルギー計算は一般的に、結合事象の熱力学的に完全な記述を提供し、少なくとも理論上、用いられる力場の限界内での正確な予測および位相空間の完全なサンプリングをもたらす。しかしながら、産業との関連における自由エネルギー計算の実際の用途では、特に複数のコンフォメーションをサンプリングするときに必要とされる大きな計算リソース、生体プロジェクト（ｌｉｖｅ ｐｒｏｊｅｃｔ）に影響を与えるのに必要なタイムスケール、およびこのような自由エネルギー計算に関連するサンプリング誤差を決定するためのロバストで有意味な手法の必要性によって制限される場合がある^{４，１１，１２}。

本発明者らは、特に、リガンドの同種組の相対結合親和性を推定し、これらの推定に関連するサンプリング誤差を決定するために、受容体−リガンドの相対結合親和性の自由エネルギー計算を改善するための一般化可能なコンピュータ実装方法を発見した。

本発明の１つの態様は、包括的には、個々のリガンド組メンバーと受容体との間で複合体を形成するための、受容体と１組のリガンドの個々のメンバーとの間の結合の相対強度を決定するコンピュータ実装方法を特徴とする。本方法は、少なくとも１つの閉じたサイクルを形成するリガンド対の組についての複数の相対結合自由エネルギー差のコンピュータ実装された決定と、各閉じた熱力学的サイクルに関するヒステリシスの大きさのコンピュータ実装された決定とを含む。概して、本方法は、
ａ．閉じた熱力学的サイクルについて観察されるヒステリシスの値を確率的に導く、閉じた熱力学的サイクルの区画の各々に沿った自由エネルギー差およびこれらの自由エネルギー差に関する誤差分布のコンピュータ実装された確率的決定と、
ｂ．ａにおいて決定された確率モデルに含まれる閉じた熱力学的サイクルにおける区画ごとの最も確からしい自由エネルギー差のコンピュータ実装された決定と、
ｃ．ｂにおける確率的決定からの、閉じた熱力学的サイクルにおける各区画に沿ったリガンド対ごとの最も確からしい自由エネルギー差に関連する最も確からしい誤差のコンピュータ実装された決定と、
ｄ．コンピュータユーザに対する、ｃにおける相対結合自由エネルギーおよび誤差を表す値のコンピュータ実装された出力と、を備える。

好ましい実施形態では、誤差を推定する工程は、２つ以上の閉じた熱力学的サイクルの区画に沿ったリガンド間での結合自由エネルギー差のコンピュータ実装された解析と、閉じた熱力学的サイクルの各々に関するヒステリシスの大きさのコンピュータ実装された確率的決定とを備える。また、好ましい実施形態では、工程ｃおよびｄは、関数

を最小にする区画ごとの１組の自由エネルギー値を、制約

を用いて決定することを含み、
ここで、Ｅｉは、所与の区画ｉについての計算される自由エネルギー差であり、
Ｆ_ｉは、所与の区画ｉについての理論上の自由エネルギー差であり、
σ_ｉは、区画ｉについての計算される自由エネルギー差の標準偏差であり、全ての閉じたサイクルについての理論上の自由エネルギー差の和は０である。

他の好ましい実施形態では、受容体はタンパク質であり、リガンドは同種であり、工程ａにおける観察されるヒステリシスの確率モデルの構築においてガウス分布が仮定され、工程ａにおける自由エネルギーシミュレーションに関連する誤差分布は一様であると仮定され、誤差分布は、工程ａにおいてベネット誤差に加算されると仮定される。閉じた熱力学的サイクルの連結性はグラフまたは行列として表される。確率的決定は、制限するものではないが、グラフ理論法、行列代数法、ベイズ法、最大尤度法を含む様々な方法によって行うことができる。

本発明の他の態様は、上記で記載した方法を実行するための有形の非一時的命令を備えるコンピュータ可読媒体、およびこれらの方法を実行するための非一時的コンピュータ可読命令を有する汎用グラフィック処理ユニットを特徴とする。

以下で、ＦＥＰ／ＲＥＳＴ（自由エネルギー摂動／溶質テンパリングを用いるレプリカ交換）^{１２〜１５}、ＦＥＰ／ＭＤ、およびＴＩ相対結合親和性自由エネルギー計算に対するサイクルクロージャ親和性推定および誤差解析機構の発明を例証するが、本発明は、ラムダ動力学、アルケミカルメタダイナミクス（ａｌｃｈｅｍｉｃａｌ ｍｅｔａｄｙｎａｍｉｃｓ）、アルケミカルＯＳＲＷ（ａｌｃｈｅｍｉｃａｌ ＯＳＲＷ）等を含む複数の
他の相対結合自由エネルギー計算プラットフォームと共に用いられ得る。また、本発明は、特定の相対自由エネルギー計算プラットフォームが、分子動力学、モンテカルロ、グランドカノニカルモンテカルロ、レプリカ交換分子動力学、加速分子動力学または任意の他のサンプリングプロトコルに基づくサンプリングアルゴリズムを利用して実施されるか否かに関わらず利用され得る。

特に、サイクリン依存性キナーゼＣＤＫ２−サイクリンＡ受容体^１６に結合する一連の同種リガンドからのサブセットを調査することによって本発明を例証する。ＣＤＫ２は、細胞増殖およびＲＮＡポリメラーゼＩＩ転写サイクルの調節において様々な機能を実行するＣＤＫファミリーのメンバーである。ＣＤＫ２は、腫瘍選択性治療戦略のための重要な薬剤標的としても同定されている^{１７，１８}。

本発明を例証するために提示する特定の系に関わらず、当業者であれば、本発明をリガンド結合調査のための一般的なツールとして用いることができることを理解するであろう。この点を説明するために、本発明を、ＪＮＫキナーゼおよびＢＡＣＥに結合するリガンドの相対結合親和性の算出にも適用した。これらもまた、医薬的に関心の高いタンパク質標的である。

本発明の１つまたは複数の実施形態の詳細が添付の図面および以下の記載において述べられている。本発明の他の特徴、目的および利点は、記載および図面、ならびに特許請求の範囲から明らかとなろう。

この特許または出願ファイルは、色彩を付して作成された少なくとも１つの図面を含んでいる。彩色図面（複数の場合もある）が付された、この特許または特許出願公開の写しは、請求および所要の手数料の納付により、特許商標庁によって提供される。

ＦＥＰ計算における一貫性および信頼性の検査の問題を示す単純なモデルの図。１組の３つのリガンドＬ１、Ｌ２およびＬ３について、それぞれ推定自由エネルギー差Ｅ２１、Ｅ３２およびＥ１３を用いた３つのＦＥＰ計算が実行された。３つのＦＥＰ計算からどのように相対自由エネルギー（Ｆ２１、Ｆ３２およびＦ１３）を入手するか、そして計算に信頼性があるか否かをどのように判定するか？ ＣＤＫ２阻害剤結合親和性計算において用いられる、同じエッジを共有する２つのサイクルを有するより複雑なパスの図。 一般的な用途における相対自由エネルギー計算のための例示的なグラフ。グラフ上の各頂点はリガンドを表す。グラフ上の各エッジは相対自由エネルギー計算を表し、エッジの矢印は相対結合自由エネルギー計算の方向を示し、各エッジ上の値は、その特定のリガンド対についての自由エネルギー計算から計算された自由エネルギー差である。 ３つの独立した部分グラフへの図３の全体グラフの分解を示す図。部分グラフは同じ頂点を共有する場合があるが、同じエッジを共有しない。 結晶構造を有する６つのＣＤＫ２リガンドの構造、ならびにＦＥＰ／ＲＥＳＴ計算およびＦＥＰ／ＭＤ計算のために用いられるパスの図。３つの赤い矢印は、計算の一貫性および信頼性を評価するための３つの別個のサイクルを完成させる相対結合自由エネルギー計算パスを表す。５つの青で示す相対結合自由エネルギー計算のためにリガンドｌｈｌｑが用いられた。なぜなら、他の化合物は全て、リガンドｌｈｌｑのベンゼン環に異なる置換基を単に付加することによって生成され得るためである。 結晶構造を有するＣＤＫ２阻害剤の組についての、ＦＥＰ／ＭＤおよびＦＥＰ／ＲＥＳＴの予測結果対実験による結合親和性データの図。ＦＥＰ／ＭＤの予測結合親和性は大幅に変動せず、これにより、これらの値を用いてリガンドｌｈｌｒ、ｌｈｌｓ、ｌｏｉｕ、ｌｏｉ９およびｌｏｉｙを確実にランク付けすることが不可能になっている。対照的に、ＦＥＰ／ＲＥＳＴの予測結果は、実験データと非常に高い相関を有し（Ｒ２＝０．９１）、全ての対について偏差は３．３４７ｋＪ／ｍｏｌ（０．８ｋｃａｌ／ｍｏｌ）未満である。この図において報告されるＦＥＰ／ＲＥＳＴ予測のためのエラーバーは、サイクルクロージャ解析を用いて本発明の実施形態によって計算される。 リガンドｌｈｌｒのためのＦＥＰ／ＲＥＳＴシミュレーションのための結晶構造に対応する開始コンフォメーション（左）、ＲＥＳＴ／ＲＥＳＴ軌道において観察される、クロロ置換フェニル環がフリップした代替結合モード（中央）、およびＸ線結晶学におけるリガンドｌｈｌｒの電子密度（右）の図。タンパク質は、主鎖について彩色されたリボンで表され、側鎖についてスティック表現で表される一方、リガンドはチューブ表現で描かれ、電子密度は無地の黄色の表面として表示される。密度の視覚化のために用いられるシグマカットオフは１．０である。リガンドｌｈｌｒの場合、電子密度は、クロロ置換ベンゼン環について２つのピークを有し、リガンドがＦＥＰ／ＲＥＳＴシミュレーションに一致して２つの結合モードを有することを示す。 リガンドｌｈｌｓのためのＦＥＰ／ＲＥＳＴシミュレーションのための結晶構造に対応する開始コンフォメーション（左）、およびＲＥＳＴ／ＲＥＳＴ軌道において観察される、代替結合モード（右）の図。タンパク質は、主鎖について彩色されたリボンで表され、側鎖についてスティック表現で表される一方、リガンドはチューブ表現で描かれる。開始コンフォメーションと比較して、スルホンアミド基はＣ−Ｓ結合の回りを約１２０°回転し、それに応じて、結合ポケットを取り囲む残基、特にＧｌｎ８５およびＡｓｐ８６が再配列する。スルホンアミド基とタンパク質残基との間に２つの水素結合が形成される。 ＣＤＫ２受容体への結晶構造結合を有しないリガンドのためのＦＥＰ／ＲＥＳＴ計算に用いられる相対結合自由エネルギー計算パスの図。 結晶構造を有する６つのリガンドおよび結晶構造を有しない１０個のリガンドの双方を含むＣＤＫ２阻害剤の全体組についての、ＦＥＰ／ＲＥＳＴの予測結果対実験による結合親和性データの図。結晶構造を有する６つのリガンド（図６に表示される）は、図において赤色で彩色され、結晶構造を有しない１０個のリガンドは黒色で彩色される。この図において報告されるＦＥＰ／ＲＥＳＴ予測のためのエラーバーは、付録Ａに詳述される方法に基づいて、サイクルクロージャを用いて計算される。 ２つの独立したサイクルクロージャを有するＪｎｋ１受容体に結合する７つのリガンドのための相対結合自由エネルギー計算パスの図。 ２つの独立したサイクルクロージャを有するＢＡＣＥに結合する６つのリガンドのための相対結合自由エネルギー計算パスの図。 活性部位における水結合（ｗａｔｅｒ ｂｏｕｎｄ）を有するシタロンデヒドラターゼおよび有しないシタロンデヒドラターゼに結合する３つのリガンドのための相対結合自由エネルギー計算パスの図。活性部位における水結合を有するシミュレーションは「（Ｈ２Ｏ）」によって示される。相対結合自由エネルギー計算パスにおいて４つの独立したサイクルクロージャが構築される。 メタン−プレート結合親和性および水中でプレートに帯電する自由エネルギーを連結する熱力学的サイクルの図。灰色の粒子は、モデルの疎水性囲い（ｈｙｄｒｏｐｈｏｂｉｃ ｅｎｃｌｏｓｕｒｅ）を形成するレナード・ジョーンズ原子を表し、赤色の粒子は負に帯電したイオンを表し、青色の粒子は正に帯電したイオンを表し、緑色の粒子は、囲いに結合する融合原子メタンを表す。 本方法を実行するコンピュータシステムのブロック図。 本方法における工程のフローチャート。

本発明によりサイクルクロージャを用いて計算される相対結合親和性は、実験データと高度に相関し、計算の収束を示す有意味なエラーバーが生成される。
これは、本発明を利用してＦＥＰ／ＲＥＳＴ自由エネルギー計算の値を求めるときに特
に明らかである。このＦＥＰ／ＲＥＳＴ自由エネルギー計算は、受容体の活性部位におけるリガンドの様々なコンフォメーションの完全なサンプリングを提供する。特定の理論に縛られることなく、本発明者らは、リガンドのうちの２つについて、複数の結合モードが存在すること、およびこれらのモードをサンプリングすることが、結合親和性の正しい予測に重要であることを割り出した。ＦＥＰ／ＲＥＳＴは比較的短いシミュレーション時間内で重要なコンフォメーションをサンプリングすることができるが、ＦＥＰ／ＭＤを用いると、リガンドは初期コンフォメーションにトラップされた。本方法は、ＦＥＰ／ＭＤ自由エネルギー計算に関連する高い誤差、およびＦＥＰ／ＲＥＳＴ計算に関連する低い誤差を示すこの特徴を説明する。このため、本発明は、計算における関連誤差の妥当な推定、および予測の信頼性を提供することによって相対自由エネルギー計算技法を改善することができる^１９。

最初に、本方法は、いくつかの他のグループによって行われたような相対自由エネルギー計算においてサイクルクロージャを構築し^{１０，２０，２１}、次に、この入力から、本方法は、閉じた熱力学サイクルごとに計算された自由エネルギーの和が理論値０からどれだけ変化しているかを評価することによって、相対自由エネルギー計算の一貫性および信頼性を決定する。次のセクションにおいて、この方法を詳細に説明する。これによって、当業者は、本方法が、どのようにサイクルクロージャ相対自由エネルギー計算によって取得された複数の自由エネルギー推定値から相対自由エネルギーの信頼性のある予測をもたらし、また、これらの計算に関連する予測誤差限界をもたらし、また、系統誤差を有する相対自由エネルギー計算にフラグを付けるメカニズムをもたらすかを直ちに理解することができる。

サイクルクロージャ相対自由エネルギー一貫性の発明の詳細な説明
Ｉ．相対結合自由エネルギー計算パスにおけるサイクルクロージャを用いた一貫性検証および最良自由エネルギー推定量
以下の単純なモデルは、本方法および本方法の概念を説明する。１組の３つのリガンドＬ_１、Ｌ_２およびＬ_３を検討する。３つのリガンド間の実験により測定された相対結合自由エネルギーの差が、

であると仮定する。ここで、

は、それぞれリガンドＬ_１、Ｌ_２およびＬ_３について実験により測定された結合自由エネルギーである。自由エネルギーは熱力学的特性であり、したがって、

である（図１を参照）。仮に、３つの相対結合自由エネルギー計算が、Ｌ_１からＬ_２まで、Ｌ_２からＬ_３まで、およびＬ_３からＬ_１まで実行されるとすると、３つの相対結合自由エネルギー計算パスについて計算される自由エネルギーの差はそれぞれ、Ｅ_２ｌ、Ｅ_３２およびＥ_１３である。シミュレーションが完全に収束し、力場が完全である場合、理想的には、

である。しかしながら、実際には、計算される相対自由エネルギーに関連する誤差が存在し、通常、

である。Δを、このサイクルクロージャに関連する自由エネルギーのヒステリシスと呼ぶ。これらの誤差は、不規則変動に起因する不偏統計誤差、および位相空間の不完全なサンプリングに起因するバイアス誤差（タンパク質および／またはリガンドはコンフォメーション空間の極小にトラップされる）、および力場における誤差を含む。本方法は、これらの計算の一貫性および信頼性を評価する。

実験値と比較した相対自由エネルギー計算における誤差は、２つのカテゴリに分けられ得る。２つのカテゴリとは、系の真の潜在的なエネルギー面と比較した、シミュレーションにおいて用いられる力場の違いから生じる系統誤差と、位相区間の不規則に不完全なサンプリングもしくは系統的に不完全なサンプリングのいずれかに起因する非収束、または自由エネルギー推定量自体、例えば、ＴＩ対ＦＥＰ等から生じる誤差である。Ｆ_２１が、無限長の不偏シミュレーションおよび不偏相対自由エネルギー推定量からの基礎をなす力場に関する２つの熱力学的状態間（すなわち、リガンドＬ_１からリガンドＬ_１まで）の理論的自由エネルギー差を表すものとする。力場に系統誤差が存在しない場合、

である。例えば、Ｌ_１からＬ_２までの計算の場合の実際の相対自由エネルギー計算では、シミュレーションは有限のサンプリング時間量を用いて実行され、サンプリングは何らか
のバイアスを有する場合があり、このため、計算される自由エネルギーＥ_２１は、自由エネルギーＥ_２１の理論値Ｆ_２１から離れる場合があり、シミュレーションの初期構成に依存することになる。

様々な初期構成、および速度の様々な乱数種から開始して無限回同じ相対自由エネルギー計算を繰り返すと、計算される自由エネルギーは分布を有する。一般性を損なうことなく、リガンドＬ_ｉおよびリガンドＬ_ｊ間の理論的自由エネルギー差がＦ_ｊｉであるリガンドＬ_ｉからリガンドＬ_ｊまでの相対結合自由エネルギーの計算に関して、Ｐ（Ｅ_ｊｉ｜Ｆ_ｊｉ）は、計算される自由エネルギーＥ_ｊｉの分布を表す。Ｐ（Ｅ_ｊｉ｜Ｆ_ｊｉ）は、原理的に、ガウス分布、ローレンツ分布、一様分布、デルタ分布等の任意の種類の分布であり得る。

この方法は、相対結合自由エネルギー計算パスにおいて構築されたサイクルクロージャを有する分布Ｐ（Ｅ_ｊｉ｜Ｆ_ｊｉ）に基づいて、リガンドＬ_ｉからリガンドＬ_ｊまでの自由エネルギー差を予測し、この予測に関する関連誤差を予測する。すなわち、

である。ここで、

はリガンドＬ_ｉとリガンドＬ_ｊとの間の予測自由エネルギー差を表す、ｅ_ｊｉは予測に関連付けられた誤差を表す。
予測は、多くの様々な方法、例えば、観察の確率を最大にする最大尤度法、計算される自由エネルギーに基づいてパラメーターを最適化するベイズ統計法等を用いて行われ得る。

ここで、最大尤度法を用いる１つの特定の例を通じて、かつ計算される自由エネルギー
についてガウス分布を仮定して、反転を例示する。当業者であれば、本発明が、他のタイプの統計的解析法および他のタイプの分布のための一般的なツールとして用いられ得ることを理解するであろう。ベイズ統計を用いた導出例も、このセクションの末尾において与えられる。

計算される自由エネルギーが、平均Ｆ_２１（系統的バイアスなし）および標準偏差σ_２１を有するガウス分布であると仮定する。このとき、単一の相対自由エネルギー計算がこのパスについてＥ_２１の値を与える確率密度は以下となる。

同様に、Ｌ_２からＬ_３へのパス、およびＬ_３からＬ_１へのパスについて、相対自由エネルギー計算がＥ_３２およびＥ_１３の値を与える確率密度はそれぞれ以下となる。

理論的自由エネルギー差Ｆ_２１、Ｆ_３２およびＦ_１３の所与の組について、３つの相対自由エネルギー計算がＥ_２１、Ｅ_３２およびＥ_１３の値を与える全体尤度Ｌは以下となる。

最大尤度法によれば、Ｆ_２１、Ｆ_３２およびＦ_１３の最も可能性の高い値は、上記の尤度を最大にする値の組である。上記の式の対数をとると、尤度を最大にする値の組は、以下の関数を最小にする値の組である。

ここで、制約：
Ｆ_２１＋Ｆ_３２＋Ｆ_１３＝０
を有する。

ラグランジュ乗数を用いると、尤度を最大にする値の組は以下となる。

上記の推定量は系統的バイアスを有さず、様々なパスからの自由エネルギー予測間には実質的に相違がない。本方法は、上記の推定量を２つのパスからの加重平均として解釈する。例えば、リガンドＬ_１とＬ_２との間の自由エネルギー差Ｆ_２１は、Ｅ_２１または−（Ｅ_３２＋Ｅ_１３）から推定することができ、最良の推定量は、２つの予測の加重平均である。パスに沿って計算された自由エネルギー差の標準偏差が小さいほど、最良推定量に対する重みが大きくなり、逆もまた同様である。

更に、上記のモデルによれば、サイクルクロージャのヒステリシスＥ_２１＋Ｅ_３２＋Ｅ_１３＝Δも、平均０および標準偏差

を有するガウス分布である。３つの計算される自由エネルギーの和が、２ｓよりも大きく０から離れている場合、計算が収束しておらず、結果に信頼性がないことがほぼ確実（Ｐ＝０．９５）であり、３つの計算される自由エネルギーの和が、ｓよりも大きく０から離れている場合、計算が収束しておらず、結果に信頼性がない場合がある可能性が高い（Ｐ＝０．６８）。このため、本方法は、予測に信頼性があるか否かを決定する。

実際の相対自由エネルギー計算では、高い計算コストに起因して、所与の計算は通常一回のみ実行され、予測ごとに標準偏差が推定され得ない。この場合、本方法は、計算ごと
の標準偏差が同じである、例えばパスごとに３．３４７ｋＪ／ｍｏｌ（０．８ｋｃａｌ／ｍｏｌ）であると仮定する。そして、２つのパスからの自由エネルギー結果の相違が５．７５７ｋＪ／ｍｏｌ（１．４ｋｃａｌ／ｍｏｌ）よりも大きい場合、これは計算が収束しておらず、おそらく信頼性がないことを示す。この仮定の下で、最良自由エネルギー推定量は下式である。

更に、計算ごとの標準偏差が同じであるという仮定の下で、サイクルクロージャのヒステリシス自体が、各自由エネルギー予測に関連する誤差の推定値を提供する。上記で検討したように、サイクルクロージャのヒステリシスΔは、平均０および標準偏差

を有するガウス分布である。このため、所与の計算がサイクルクロージャについて値Δを有するヒステリシスを生成する確率は以下となる。

上記の確率を最大にするシグマの値は、各自由エネルギー計算に関連する標準偏差の最大尤度推定値を与える。

上記の導出において、本方法が推定する自由エネルギーおよび関連誤差は、相対自由エネルギー計算から得た全ての情報に基づく最適推定値である。計算における基礎をなす力場が系統的バイアスを有する場合、完全に収束した相対自由エネルギー計算結果であっても、実験により測定された値から離れる場合がある。これは上記の解析に基づいて補正され得ないものである。

上記で記載した方法を、より多くのメンバーをサイクルクロージャ内に有し、いくつかのサイクルが同じエッジを共有する、より複雑な事例に適用する。ここで、結晶構造を有する６つのＣＤＫ２リガンド（図２）のためのＦＥＰ／ＭＤ計算およびＦＥＰ／ＲＥＳＴ計算において用いられる、パスのための最良自由エネルギー推定量を取得する方法を説明する。任意の種類の相対結合自由エネルギー計算パスのためのより一般的なアルゴリズムが以下のセクションで与えられる。

ベイズ統計を用いた最良自由エネルギー推定量
リガンドＬ_ｉからリガンドＬ_ｊに、理論的自由エネルギー差Ｆ_ｊｉを有して突然変異するＦＥＰ計算について、仮に、計算される自由エネルギーＥ_ｊｉが、平均Ｆ_ｊｉおよび標準偏差σ_ｊｉを有するガウス分布であるものとする（最大尤度法の場合に用いられるのと同じモデル）。このとき、１つのＦＥＰ計算から計算される自由エネルギーがＥ_ｊｉの値を与える確率密度は以下となる。

Ｆ_ｊｉの多くの様々な値は、様々な確率を有する同じ計算された自由エネルギーＥ_ｊｉを導くことができる。ベイズ統計によれば、本方法は、ＦＥＰ計算を実行する前に、２つのリガンド間の自由エネルギー差の分布（事前分布）の推定値

を有し、ＦＥＰ計算を実行した後、ＦＥＰからの計算結果に基づいて分布（事後分布

を調整する。
ＦＥＰ計算を実行する前に、相対結合自由エネルギー計算に関する情報が知られていない場合、Ｆ_ｊｉの事前分布に関する妥当な推定は、Ｆ_ｊｉが−∞と＋∞との間で一様分布するというものであり、すなわち、

である。事前分布は正規化されないことに留意されたい。
ＦＥＰ計算がＥ_ｊｉの値を与える場合、ベイズ統計によれば、Ｆ_ｊｉの事後確率は以下となる。

上記の事例の場合、３つのリガンドＬ_１、Ｌ_２およびＬ_３を用いて、本方法は、計算された自由エネルギー差Ｅ_２１、Ｅ_３２およびＥ_１３を有する閉じたサイクルを形成するＦＥＰ計算を実行する。このとき、ベイズ統計によれば、（Ｆ_２１，Ｆ_３２，Ｆ_１３）の事後分布は、
Ｆ_２１＋Ｆ_３２＋Ｆ_１３＝０
という制約の下で以下となる。

事後分布のピックは、

に位置する。
これにより、上記で提示した最大尤度法から取得した結果を確認する。
ＩＩ．サイクルクロージャにより導出された自由エネルギー推定値および誤差推定値の発明の一般的な実施のプログラミング
１組のＬ個のリガンドについて、本方法は全てのリガンドをサイクルクロージャ（複数の場合もある）と連結するＮ個の相対自由エネルギー計算を実行する。このセクションに記載される本発明の特定の実施形態は、リガンド間の相対自由エネルギーの最良予測およびこれらの予測に関連するエラーバーをもたらす。入力は、相対自由エネルギーシミュレーションが実行されるリガンドのＮ個の対間のＮ個の計算される自由エネルギー差であり、出力は、相対自由エネルギーシミュレーションが実行されたリガンドのＮ個の対間のＮ個の予測相対自由エネルギー差および関連するエラーバー、ならびにまた、組内のリガンドの任意の対間の相対自由エネルギー差およびこれらの関連するエラーバーである。予測相対自由エネルギーはヒステリシスを有さず、これは、任意の２つのリガンド間の予測相対自由エネルギー差が、これらを連結するパスから独立しており、このため内部的に矛盾しないことを意味する。

図３に示す例示的なパスにおいて、８つ全てのリガンドをいくつかのサイクルクロージャに連結する、１組の８つのリガンドのための１１個の相対自由エネルギー計算が存在する。プログラムの入力は、リガンドの１１個の対についての１１個の相対自由エネルギー計算結果であり、出力は、リガンドの１１個の対についての予測相対自由エネルギー差、これらの関連するエラーバー、ならびにまた、組内の８つのリガンドの任意の対間の相対自由エネルギー差およびこれらの関連するエラーバーである。

入力：
｛Ｅ_ｉ；Ｌ_ｊ→Ｌ_ｋ｝、ここでｉ∈｛１，２，３，．．．Ｎ｝
出力：
｛Ｆ_ｉ，ｅ_ｉ；Ｌ_ｊ→Ｌ_ｋ｝、ここでｉ∈｛１，２，３，．．．Ｎ｝
および｛Ｆ_ｉｊ，ｅ_ｉｊ｝、ここでｉ，ｊ∈｛Ｌ_２．．．Ｌ｝
アルゴリズム
工程１：実行される全ての相対結合自由エネルギー計算について計算される自由エネルギー、ならびにこれらの対応する初期リガンドおよび最終リガンドを列挙する。

この工程から入手するものは以下のリストである。
｛Ｅ_ｉ；Ｌ_ｊ→Ｌ_ｋ｝、ここでｉ∈｛１，２，３，．．．Ｎ｝
ここで、Ｅ_ｉは第ｉの相対結合自由エネルギー計算についての、相対自由エネルギーシミュレーションにより計算される自由エネルギーである。この相対自由エネルギー計算についての初期リガンドおよび最終リガンドはＬ_ｊおよびＬ_ｋである。これはプログラムの入力である。

工程２：相対結合自由エネルギー計算パスにおける全ての閉じたサイクルを列挙し、全体グラフについてサイクルクロージャ連結性行列（Ｃ３Ｍ）を組み立てる。
この工程において、本方法は、全てのリガンド対にわたって探索し、閉じたサイクルを形成する全ての対を発見し、この情報を連結性行列に変換する。

例えば、上記に示すリガンド対パスにおいて、対１、２および５は閉じたサイクルを形成する。自由エネルギーは状態関数であるので、そのような閉じたサイクルに沿った自由エネルギーの和は０となるべきである。換言すれば、
Ｆ_１＋Ｆ_２−Ｆ_５＝０
である。

一般的な事例では、閉じたサイクルに沿った自由エネルギーに対する拘束は、以下の一般的な拘束関数に書かれ得る。

ここで、Ｆ_ｉは、閉じたサイクルに沿ったエッジの２つの頂点上の２つのリガンド間の自由エネルギー差であり、Ｍ_ｉは、エッジが閉じたサイクルと同じ方向にある場合、１であり、エッジが閉じたサイクルの反対方向にある場合、−１であり、エッジが閉じたサイクル内にない場合、０である。

この工程において、本方法は、リガンド対パス内の全ての閉じたサイクルおよび全ての拘束を列挙する。上記で示す例において、閉じたサイクルについて以下の拘束を有する。
Ｆ_１＋Ｆ_２−Ｆ_５＝０
Ｆ_３＋Ｆ_４−Ｆ_５＝０
Ｆ_４−Ｆ_７−Ｆ_６＝０
Ｆ_１＋Ｆ_２−Ｆ_４−Ｆ_３＝０
Ｆ_５−Ｆ_７−Ｆ_６−Ｆ_３＝０
Ｆ_１＋Ｆ_２−Ｆ_７−Ｆ_６−Ｆ_３＝０
Ｆ_９＋Ｆ_１１−Ｆ_１０＝０
一般的な事例において、仮に、リガンド対パス内に合計Ｍ個の閉じたサイクルがあるとすると、Ｍ個の拘束方程式を有する。

Ｍ_ｊｉ＝Ｍ_ＧはＭ×Ｎの行列であり、行列の各行は閉じたサイクルのための拘束方程式を表す。Ｍ_Ｇを用いて、全体グラフのサイクルクロージャ連結性行列（Ｃ３Ｍ）を表す。
上記で示す例において、

である。
工程３：全体グラフを独立した部分グラフに分解する。
相対自由エネルギー計算グラフにおいて、２つの部分グラフを連結する閉じたサイクルがない場合、リガンドのサブセット間の予測相対自由エネルギーは、リガンドの他のサブセット間の予測相対自由エネルギーから独立している。上記に示す例において、リガンドの組［１，２，３，４，５］間の予測相対自由エネルギーは、リガンドの組［６，７，８］間の相対自由エネルギーから独立している。なぜなら、２つの部分グラフを連結する閉じたサイクルがないためである。この工程において、全体グラフは各独立した部分グラフに分解され、次に、各独立した部分グラフが連続して解かれる。

工程３ａ：Ｍ_Ｇの各行を検査し、各行における非ゼロ要素の列数を記録し、列数を以下のリストに入れる。

この工程において、本方法はＭ個のリストを生成する。各リストは、Ｃ３Ｍ行列の各行の非ゼロ要素の列数を含む。これを以下の記載について明確にするために、この工程において生成される、未加工リストであるＭ個のリストを含む組は、組Ｓ_ｒと呼ばれる。

上記で示す例において、本方法のこの工程は、７つのリストを含む未加工リスト組を生成する。
Ｌ_１＝｛１，２，５｝
Ｌ_２＝｛３，４，５｝
Ｌ_３＝｛４，６，７｝
Ｌ_４＝｛１，２，３，４｝
Ｌ_５＝｛３，５，６，７｝
Ｌ_６＝｛１，２，３，６，７｝
Ｌ_７＝｛９，１０，１１｝
工程３ｂ：未加工リスト組Ｓ_ｒにおいて第２のリストから最後のリストまで反復し、Ｌ_ｊに含まれる要素を第１のリストＬ_１内の要素と比較する。ここで、ｊ∈｛２，３，．．．ｍ｝である。これらが１つの共通要素（または数個の共通要素）を含む場合、Ｌ_ｊ内の全ての要素をマージして第１のリストＬ_１にし、未加工リスト組からＬ_ｊを削除する。これらがコメント要素を一切有しない場合、リストＬ_ｊを未加工リスト組内に維持する。

上記に示す例において、Ｌ_２、Ｌ_４、Ｌ_５、Ｌ_６は、Ｌ_１との共通要素（複数の場合もある）を含むので、これらのリストはこの工程において削除され、これらのリストの要素はＬ_１にマージされる。この工程から未加工リスト組内に残るものは、以下である。

Ｌ_１＝｛１，２，３，４，５，６，７｝
Ｌ_２＝｛４，６，７｝
Ｌ_７＝｛９，１０，１１｝
工程３ｃ：工程３ｂにおいて削除されていないリストについて工程３ｂを繰り返し、要素がＬ_１との共通要素を含む場合、これらの要素をＬ_１に統合する。これ以上リストが削除され得なくなるまでこの工程を繰り返す。

上記で示す例において、この工程から未加工リスト組Ｓ_ｒ内に残るものは、以下である。
Ｌ_１＝｛１，２，３，４，５，６，７｝
Ｌ_７＝｛９，１０，１１｝
工程３ｄ：未加工リスト組Ｓ_ｒ内の第１のリストＬ_１を最終リスト組Ｓ_ｆ＝｛Ｌ_１．．．｝内に移し、これを未加工リスト組から削除する。

上記で示す例において、この工程の後、最終リスト組Ｓにおいて、
Ｓ_ｆ＝｛Ｌ_１＝｛１，２，３，４，５，６，７｝｝
であり、未加工リスト組内に１つのリストが残される。

Ｓ_ｒ＝｛Ｌ_７＝｛９，１０，１１｝｝
工程３ｅ：工程３ｄの後、未加工リスト組内に２つ以上のリストが残されている場合、残されているリストが１つのみになるまで、工程３ｂ〜工程３ｄを繰り返す。

上記に示す例では、この工程において動作は必要ない。なぜなら、工程３ｄの後、未加工リスト組内には１つのリストしか残されていないからである。
工程３ｆ：未加工リスト組Ｓ_ｒ内の最後のリストを、最終リスト組Ｓ_ｆに移す。この工程の後、最終リスト組内にいくつかのリストが存在し、未加工リスト組内にはリストが残されていない。最終リスト組内のリストは、共通要素を一切含まない。

上記に示す例において、最終リスト組内に以下を有する。
Ｓ_ｆ＝｛Ｌ_ｉ＝｛１，２，３，４，５，６，７｝，Ｌ_２＝｛９，１０，１１｝｝
工程３ｇ：列番号１〜Ｎにわたって反復する。列番号が最終リスト組内の任意のリストにある場合、これをスキップし、そうでない場合、その列番号しか含まないリストを最終リスト組に加える。

これはこのプロセスにおける最終工程である。この工程の後、最終リスト組はいくつかのリストを含み、リストの内のいずれの間においても共通要素は存在しない。列番号１〜Ｎは、最終リスト組における１つのリスト内にのみある。

上記に示す例において、この工程の後、
Ｓ_ｆ＝｛Ｌ_１＝｛１，２，３，４，５，６，７｝，Ｌ_２＝｛９，１０，１１｝，Ｌ_３＝｛８｝｝
である。

このプロセスの後、全体グラフは、各個々の独立した部分グラフに分解される。各独立した部分グラフはエッジを含み、エッジ番号は最終リスト組Ｓ_ｆの各リスト内に列挙されている。これらの独立した部分グラフは共有エッジを含まないが、共有頂点を有する場合がある。

上記に示す例において、全体グラフは３つの独立した部分グラフに分解される。第１の部分グラフは７つの結合自由エネルギー計算パス［Ｅｌ，Ｅ２，Ｅ３，Ｅ４，Ｅ５，Ｅ６，Ｅ７］を含み、第２の部分グラフは３つの相対結合自由エネルギー計算パス［Ｅ９，Ｅ１０，Ｅ１１］を含み、第３の部分グラフは１つの相対結合自由エネルギー計算パス［Ｅ８］を含む。これらの部分グラフは図４に示される。

工程４：独立した部分グラフごとのＣ３ＭであるＭ_Ｓを入手する。
工程３から、全体連結性グラフは既に独立した部分グラフに分解されており、各部分グラフのエッジ番号は最終リスト組内のリストに記憶されている。この工程において、独立した部分グラフごとのＣ３Ｍを入手することに着手する。

仮に、リスト内にＮ_ｉ個の要素が存在するとして、最終リスト組Ｓ_ｆから各リストＬ_ｉを取る。全体グラフのＣ３ＭであるＭ_Ｇの全ての行にわたって探索する。Ｍ_Ｇの行ごとに、列番号がリストＬ_ｉ内にある全てのｎ個の要素を探索し、これらのうちのいずれかがゼロでない場合、この行のｎ個の列をＭ^ｉ _Ｓに入れる。この部分グラフのためのＣ３ＭはリストＬ_ｉに対応している。これを、Ｍ_Ｇの全ての行、および最終リスト組Ｓ_ｆ内の全てのリストに対して行う。

この工程から、全体グループのＣ３Ｍである、Ｍ×Ｎの行列であるＭ_Ｇが、それぞれが独立した各部分グラフに対応する各独立したＣ３Ｍである、Ｍ_ｉ＝Ｎ_ｉ行列であるＭ^ｉ _Ｓに分解される。これらの行列の次元は以下の特性を満たす。

この工程の後、全体グラフは独立した部分グラフに分解され、全体グラフのＣ３ＭであるＭ_Ｇも、独立した部分グラフごとのＣ３ＭであるＭ^ｉ _Ｓに分解される。
Ｍ^ｉ _Ｓ内の各行が、パス内の閉じたサイクルを表し、ある独立した部分グラフのための連結性行列Ｍ^ｉ _Ｓは、ゼロ要素を含む場合があり、これは、その部分グラフ内に閉じたサイクルがないことを意味することに留意されたい。

上記の例において、３つの独立した部分グラフのためのＣ３Ｍは、以下である。

第３の部分グラフのためのＭ_Ｓは空であることに留意されたい。
工程５：部分グラフＭ_ＳごとにＣ３Ｍを独立した行列Ｍ_ＳＲに縮約する。
部分グラフごとに工程４から生成されたＣ３Ｍにおいて、いくつかの行は互いに独立しておらず、これは、Ｍ_Ｓ内に列挙された閉じたサイクルが独立していないことを意味する。例えば、上記で示す例において、部分グラフ１について、３つの独立したサイクルが存在するが、６つ全ての閉じたサイクルが連結性行列Ｍ_Ｓ ^１に列挙されている。それらのうちの３つは必ずしも必要でない。

この工程において、各部分グラフのＣ３ＭであるＭ^ｉ _Ｓは、独立した行ベクトルのみを含む行列に縮約される。独立した行ベクトルのみを含む変換された行列をＲＣ３Ｍ（縮約されたサイクルクロージャ連結性行列（Ｒｅｄｕｃｅｄ Ｃｙｃｌｅ Ｃｌｏｓｕｒｅ Ｃｏｎｎｅｃｔｉｖｉｔｙ Ｍａｔｒｉｘ））と呼び、Ｍ^ｉ _ＳＲと表す。

行列Ｍ^ｉ _Ｓを行列Ｍ^ｉ _ＳＲに縮約する多くの方法が存在する。ここでは単純な方法を列挙する。
工程５ａ：Ｍ^ｉ _Ｓの最初の２つの行に新たな２×Ｎ_ｉ行列Ｍを形成させる。

工程５ｂ：ｄｅｔ｜ＭＭ^Ｔ｜を計算する。抑制がゼロでない場合、２つの行は独立している。抑制がゼロの場合、２つの行は相互依存しており、第２の行を削除する。
工程５ｃ：Ｍ^ｉ _Ｓの第３の行をＭに入れ、工程５ｂを繰り返す。

工程５ｄ：Ｍ^ｉ _Ｓ内の全ての残りの行について工程５ｃを繰り返す。
この工程の後、部分グラフごとにＭ^ｉ _ＳＲが取得される。Ｍ^ｉ _ＳＲは独立した行ベクトルのみを含み、これは、独立した閉じたサイクルのみが行列内に列挙されていることを意味する。

上記の例において、第１の部分グラフのためのＭ_ＳＲは以下である。

Ｃ３Ｍ内の行の順序および行列を縮約するアルゴリズムに依拠して、部分グラフごとのＭ_ＳＲは一意ではないが、Ｍ_ＳＲ内の行数は、行の順序、および縮約を行うアルゴリズムに対し独立している。

この工程の後、縮約されたサイクルクロージャ連結性行列Ｍ_ＳＲが独立した部分グラフごとに取得される。Ｍ_ＳＲは各部分グラフ内に独立したサイクルのみを含む。以下において、独立した部分グラフごとにサイクルクロージャ方程式を解く。

工程６：独立した部分グラフごとにサイクルクロージャ方程式を解く。
この工程は、部分グラフごとに、ヒステリシスを有しないリガンドのＮ_ｉ個の対について、自由エネルギーの１組の最適推定値を生成し、また、確率を最大にする。

２つの異なる事例が存在する。
事例ａ：部分グラフ内に閉じたサイクルが存在しない。このとき、この事例に対する解は以下となる。

Ｆ_ｉ＝Ｅ_ｉ
ｅ_ｉ＝０
上記に示す例において、第３の部分グラフは閉じたサイクルを一切有しない。このため、第３の部分グラフのための解は以下となる。

Ｆ_８＝Ｅ_８
ｅ_Ｒ＝０
事例ｂ：部分グラフ内にｍ個の独立した閉じたサイクルが存在する。

この事例では、ｎ個の要素および対応するＲＣ３Ｍである、Ｍ＝Ｍ_ＳＲ（ｍ×ｎ行列）を有する最終リスト組Ｓ_ｆ内のリストＬ_ｓに対応する部分グラフごとに、以下の関数、

ここで、ｉ∈Ｌ_ｓ
が、以下の制約、

ここで、ｊ∈｛１，２．．．ｍ｝、ｉ∈Ｌ_ｓ
の下で最小化される。
ラグランジュ乗数を用いて、以下の関数、

を、以下の制約、

ここで、ｊ∈｛１，２．．．ｍ｝、およびｉ∈Ｌ_ｓ
の下で最小化する必要がある場合がある。ここで、Ｃ_ｊはラグランジュ乗数の係数である。

式

を解いて、本方法は、

を確立する。
ここで、Ｍ^Ｔ _ｊｉはＭ_ｊｉの転置である。
行列式において書かれると、

であり、以下の２組の式がもたらされる。

式（１）を式（２）に代入して、以下の式を得る。

ＲＣ３Ｍ行列ＭおよびベクトルＥの双方が知られているので、式が解かれ、係数ベクトルＣを得ることができる。係数が解かれた後、式（１）は、シミュレーションのｎ個の対について最良の自由エネルギー推定量をもたらす。

工程７：可能な相対結合自由エネルギー計算パスごとに誤差推定値を得る。
上記で導入された仮説によれば、閉じたサイクルに沿った相対自由エネルギー計算の和は、標準偏差

を有するガウス分布である。ここで、ｎは閉じたサイクル上のメンバー数である。この閉じたサイクル上のリガンドの対間の計算される差ごとの誤差は、

であり、ここで、Δは閉じたサイクルのヒステリシスである。
閉じたサイクル上のリガンドの対ごとに、閉じたサイクル上のリガンドの対間の計算される自由エネルギー差の誤差が、そのサイクルのヒステリシスによって推定され得る。

より誘導形の形式で書かれると、Ｃ３ＭであるＭ_Ｓ（ここでは縮約されていないサイクルクロージャ連結性行列であることに留意されたい）の行ごとに、以下の誤差推定値がサイクル上の各メンバーに割り当てられ得る。

全ての閉じたサイクルについて、すなわちＭ_Ｓ内の全ての行について式（４）に従って誤差計算を行った後、相対結合自由エネルギー推定値ｉのための最終誤差推定値が全ての行からの推定値の最大値として選択される。換言すれば、

である。
工程６および工程７の後、ｎ個の相対自由エネルギーのための最良自由エネルギー推定値、およびまた独立した部分グラフごとのこれらの関連するエラーバーを得る。

工程８：組内のリガンドの任意の対間の自由エネルギー差の推定値およびこれらの関連するエラーバーを報告する。
前の工程において、本方法は既に、各独立した部分グラフのリガンドのｎ個の対について、最良自由エネルギー推定値およびこれらの関連するエラーバーを提供している。全ての独立した部分グラフからこれらのデータを共にマージすることによって、以下の出力が生じる。

上記で列挙されていないリガンドの任意の対について、本方法は、２つのリガンドを連結することができるグラフ上の全ての可能なパスを列挙し、次に、パスの各エッジ上の自由エネルギーおよびまた誤差を総和する。自由エネルギーはパスから独立しているが、誤差はパスに依拠する場合がある。報告する全ての誤差の最大値は、具体的な式で書かれ得る。

ここで、

であり、最大値は、リガンドＬ_ｉおよびリガンドＬ_ｊを連結する全ての可能なパスについて取られる。この工程において、自由エネルギー推定値および誤差は、自由エネルギー計算が実行されない（グラフにおいて直接連結されていない）リガンドの対についてのみ計算されることに留意されたい。自由エネルギー計算が実行されるリガンドの対についての自由エネルギー推定値および誤差は、前の工程において計算される。

当業者にとって、上記において、閉じたサイクル（複数の場合もある）について観察されるヒステリシス（複数の場合もある）を引き起こす、確からしい自由エネルギー差およびこれらの誤差の確率モデルのコンピュータ実装を開発し、また、そのような確率モデルを所与として、自由エネルギー差およびこれらの誤差の最適な予測を行うための、多くの様々な手法が取られ得ることが明らかであるはずである。例えば制限するものではないが、行列またはグラフによる方法を用いて、必要とされるデータ構造を作成することができ、行列代数またはグラフ理論による方法を用いて確率モデルを構築することができ，最大尤度法またはベイズ法を用いて、ヒステリシス（複数の場合もある）と矛盾しないようにモデルをパラメーター化することができる。先行するセクションにおける代数行列式は、本発明の特定の実施形態にすぎず、経路の区画（ｌｅｇ）の最も確からしい誤差がベネット誤差、ブートストラップ誤差または双方に加算されると仮定することを含む、多くの方法で容易に一般化され得る。同様に、上記代数行列式は、ガウス分布だけでなく、制限するものではないが、ローレンツ分布、レヴィ分布、レイリー分布、デルタ分布および一様分布等の他の分布を用いるように一般化され得る。

したがって、閉じたサイクル（複数の場合もある）について観察されるヒステリシス（複数の場合もある）を引き起こす、確からしい自由エネルギー差およびこれらの誤差の特定の確率モデルを構築し、また、そのような確率モデルを所与として、自由エネルギー差およびこれらの誤差の最適な予測を行うために検討され得る多岐にわたる手法は、本明細書に記載される本発明の興味深い実施形態であろう。

様々な特定の受容体およびリガンドに対する様々な特定の相対結合自由エネルギー計算への本発明の適用の例示的な実施例
ＦＥＰ／ＲＥＳＴ相対自由エネルギー計算アルゴリズム
本発明は、相対自由エネルギー計算（例えば、明確な溶媒モデル計算）の精度および信頼性を改善する。そのような精度は、異なるリガンドが結合するときに結合ポケットを取り囲むタンパク質残基が大きく動く場合、または２つのリガンドが異なる結合モードを採用するか、もしくは１つもしくは複数の変性結合モードが現れる場合に、達成するのが特に困難である^{４，１１，１９}。このような事例では、全ての関係するコンフォメーションをサンプリングすることによって、自由エネルギー差の推定値の精度が改善するが、異なるコンフォメーション間で交換するための移行時間が、ほとんどのＦＥＰ／ＭＤシミュレーション長のタイムスケールでアクセス可能となるには長すぎる場合がある。

ここで、本発明者らは、近年開発された相対自由エネルギー計算法であるＦＥＰ／ＲＥ
ＳＴを用いる方法の有用性を説明する、相対結合自由エネルギー計算の特定の組を提示する。ワンら（Ｗａｎｇ ｅｔ ａｌ．）^１２によって導入されたＦＥＰ／ＲＥＳＴは、強化されたサンプリング法であるＲＥＳＴを、効率的なラムダホッピングプロトコルを通じてＦＥＰ計算に組み込み、容易にアクセス可能なシミュレーション時間内で相対タンパク質−リガンド結合親和性計算における関係するローカル構造再配列をサンプリングする。本発明者らはまた、特定の選択された相対結合自由エネルギー計算方法と独立して、本発明の有用性を同様に説明する、より定型的なＦＥＰ／ＭＤ相対自由エネルギー計算も提示する。

ＦＥＰは、λ値によって表される一連の別個の工程において系間でアルケミカル的に（ａｌｃｈｅｍｉｃａｌｌｙ）変換することによって、２つの系間の自由エネルギー差を計算することを伴う。ここで、λは、初期状態の場合の０から最終状態の場合の１まで変動する。ＦＥＰ／ＲＥＳＴは、タンパク質残基およびリガンドを含むことができる結合ポケット（「ホット」領域と呼ぶ）を取り囲む局所領域のための潜在エネルギーを変更する。中間λウィンドウについて、局所領域のための潜在エネルギーは１未満の係数でスケーリングされる。このようにして、エネルギー障壁が下げられ、これらの中間λウィンドウにおいて様々なコンフォメーションの効率的なサンプリングが可能になり、様々なコンフォメーションがハミルトニアンレプリカ交換法を通じて最終状態に伝搬される^{１５，２２，２３}。潜在的なエネルギーがスケーリングされる領域を「ホット」領域と呼ぶ。なぜなら、潜在的なエネルギーへのローカル寄与が、その領域の有効温度が高くなるように、より小さな値にスケーリングされているためである。しかし、この類推は、シミュレーション内においてより高い運動エネルギーを実際に有する粒子には引き継がない。ＦＥＰ／ＲＥＳＴアルゴリズムに関する詳細は、本願明細書に援用される、ワンら（Ｗａｎｇ ｅｔ ａｌ．）^１２によるＦＥＰ／ＲＥＳＴの研究において記載されている。

ＲＥＳＴ強化型サンプリングのための「ホット」領域は、ＳＡＲの研究において変異したリガンド官能基を含むように構成される。リガンド官能基が芳香環に付着されると、芳香環も「ホット」領域内に含まれ、環がＲＥＳＴ強化型サンプリングを通じてフリップすることが可能になる。ＲＥＳＴ「ホット」領域のための効率的な温度プロファイルは、参考文献２６の式９に従って計算され^２４、レプリカ交換の予測受容比を０．３になるように設定し、最適値は、必要なレプリカ数と交換効率との間のトレードオフを反映する。通常のＦＥＰ計算およびＦＥＰ／ＲＥＳＴ計算の双方について合計１２個のλウィンドウが用いられる。初期状態に一意な静電相互作用は、レナード・ジョーンズ（ＬＪ）相互作用の前にオフにされ、最終状態に一意なＬＪ相互作用は静電相互作用の前にオンにされた。結合された相互作用は、初期状態から平滑に最終状態にされる。ＬＪ相互作用のコアは、シミュレーションにおける特異性および不安定性を回避するように軟化された^２５。本発明による現行のシミュレーションにおいて、ＳＡＲの研究^１６において置換基が変更されたベンゼン環が「ホット」領域に含められ、合計で１２個のλウィンドウが達成される最高有効温度は概ね９００Ｋである。

ＦＥＰ／ＲＥＳＴおよびＦＥＰ／ＭＤの様々な相対自由エネルギー計算の詳細
デズモンド（Ｄｅｓｍｏｎｄ）^{２６，２７}において実施されるリガンド官能基変異モジュールが計算をセットアップするのに用いられた^２８。ＯＰＬＳ２００５力場^{２９，３０}がＳＰＣ^３１水モデルと共にタンパク質およびリガンドのために用いられた。シミュレーションのための開始構造は、ＩＤ ｌｈｌｑ、ｌｈｌｒ、ｌｈｌｓ、ｌｏｉｕ、ｌｏｉ９およびｌｏｉｙを有するＰＤＢ構造から取られた^１６。ＰＤＢ内に複数のコンフォメーションが存在した場合はいつでも、シミュレーションにおいて用いるために第１のコンフォメーションを選択した。タンパク質がタンパク質調剤ウィザードを用いて調剤され^{３２，３３}、この間、７．０のｐＨを仮定してプロトン化状態が割り当てられた^３４。サイクリンＡはＣＤＫ２と共に保持され、全ての計算について、受容体の一部として扱われた。Ｎ
ａ^＋イオンまたはＣＩ⁻イオンが、電気的中性を維持するために加えられた。系は、デフォルトのデズモンド緩和プロトコルを用いて緩和および平衡化された。デズモンド緩和プロトコルは、一連の最小化および短い分子動力学シミュレーションを伴う。通常のＦＥＰ計算およびＦＥＰ／ＲＥＳＴ計算の双方について合計１２個のλウィンドウが用いられた。製造段階は、複雑なシミュレーションの場合、２ｎｓ続き、溶媒シミュレーションの場合、５ｎｓ続いた。近傍のラムダウィンドウ間のレプリカ交換が１．２ｐｓおきに試行された。製造段階は、ＦＥＰ／ＭＤおよびＦＥＰ／ＲＥＳＴを用いてサンプリングされた。ベネット受容比法（ＢＡＲ）^３５を用いて自由エネルギーが計算された。ブートストラッピング^{１９，３６}およびＢＡＲ解析誤差予測^{１９，３６，３７}の双方を用いて自由エネルギー計算ごとに誤差が推定され、２つの誤差のうちの大きい方が報告された。

本発明の特定の適用による予測の一貫性および信頼性の評価
最も単純には、共通参照リガンドから全てのＦＥＰシミュレーションを実行することによって、相対タンパク質−リガンド結合親和性ＦＥＰ計算におけるリガンドの全体組について相対結合親和性を順位付けすることができる。更に、計算ごとに完全に収束した結果を取得する場合、最終的な予測は理論的に、相対結合自由エネルギー計算を行うために選択されるパスと独立している。例えば、３つのリガンドＡ、ＢおよびＣの場合、リガンドＡとリガンドＣとの間で相対結合親和性を取得するための２つの戦略、すなわち（１）ＡからＣまでパスを直接サンプリングする、または（２）ＡからＢおよびＢからＣの２工程でパスをサンプリングし、次に、２つの取得した自由エネルギーを合算する、が存在する。理論的に、上記の２つの方法からの最終的な収束した自由エネルギー推定値は同じであるはずである。一方、実際には、ベネット誤差限界によって特徴付けることができる位相空間のサンプリングにおける不規則変動に起因する不偏誤差と、ベネット誤差限界によって特徴付けされない位相空間の系統的に不完全なサンプリングに起因したバイアス誤差とを含む、各計算における誤差に起因して、上記の２つのパスからの分散および平均値の双方が、通例、幾分異なっている^３８。これらの制限を所与として、取得した結果の信頼性を評価することが重要である。

予測の一貫性および信頼性に関する追加情報を得るために、全てのパスについてサイクルクロージャを含む。上記で与える例において、２つのパスからの結果が指定された閾値以内まで一致する場合、これは予測に一貫性があり、利用される力場の制限内で信頼性がある可能性が高いことを示す。一方、２つのパスからの結果が、ある閾値を超えて異なっている場合、これは、結果が収束していない場合があり、予測に信頼性がない場合があることを示す場合がある。閾値は、各サイクル内のメンバーの数、およびサイクル内の区画について予期する誤差に依拠する。本発明のこの特定の実施形態は、予測の信頼性、およびサイクルごとの適切な閾値をどのように適応的に割り当てるかを厳密に評価する。更に、パス内にサイクルクロージャを含めると、２つの状態間の自由エネルギー差を計算する複数の独立した方法が存在する。本発明のこの特定の実施形態では、閉じたサイクルに沿って取得された複数の自由エネルギー推定値からの最適予測も報告される。この実施形態において、本発明はまた、各サイクルのヒステリシスからの自由エネルギー計算に関連する収束誤差に関する情報も直接提供する。

本明細書において検討される方法は、これらの工程を実行するようにプログラムされたコンピュータシステムによって実行される。図１５において、システムは、コンピュータ（複数の場合もある）１１０と、データストア１２０と、プログラムモジュール１３０とを備える。コンピュータ命令は、１つまたは複数の非一時的コンピュータ可読媒体に記憶される。図１６は、全ての自由エネルギー計算パスについての計算される自由エネルギーの列挙（３０２）と、全ての閉じたサイクルの列挙およびサイクルクロージャ連結性行列の組み立て（３０４）と、独立した部分グラフへの全体グラフの分解（３０６）とを含む本方法における工程を示す。

以下で、ある特定の例を検討し、以下のように関係する表を提示する。
表１：参照リガンドｌｈｌｑと比較した、結晶構造を有する５つのリガンドについてＦＥＰ／ＭＤおよびＦＥＰ／ＲＥＳＴにより予測される相対結合自由エネルギー。
表２：本発明の実施形態によって導出される式を用いるＦＥＰ／ＭＤ計算およびＦＥＰ／ＲＥＳＴ計算からの、結晶構造を有する６つのリガンドについての最終自由エネルギー予測。
表３：ＦＥＰ／ＲＥＳＴを用いて研究されるＣＤＫ２阻害剤およびこれらの化合物＃の構造。
表４：結晶構造を有しない追加の１０個の化合物についての、ＦＥＰ／ＲＥＳＴにより予測される相対結合自由エネルギー。
表５：本発明の実施形態によって導出される式を用いる、リガンドの全体組についての最終自由エネルギー予測。
表６：ｊｎｋ１リガンドの８つの相対結合自由エネルギー計算パスについてのＦＥＰ／ＲＥＳＴ結果。
表７：ｊｎｋ１に結合している７つのリガンドについての予測結合親和性およびそれらの関連誤差。
表８：ＢＡＣＥリガンドのための７つの７相対結合自由エネルギー計算パスについてのＦＥＰ／ＲＥＳＴ結果。
表９：ＢＡＣＥに結合している６つのリガンドについての予測結合親和性およびそれらの関連誤差。
表１０：シタロンデヒドラターゼリガンドの９個の相対結合自由エネルギー計算パスについてのＭＣ／ＦＥＰ結果。
表１１：シタロンデヒドラターゼに結合している３つのリガンドについての予測結合親和性およびそれらの関連誤差。
表１２：図１４において描かれる４つの状態間のサイクルクロージャにより予測される自由エネルギー差およびそれらの関連誤差。

結果および検討
結晶構造を有するリガンドにおける性能
図５に表示される、ＣＤＫ２−サイクリンＡ受容体に結合している６つのリガンドについて、相対結合親和性を予測するためのＦＥＰ／ＭＤシミュレーションおよびＦＥＰ／ＲＥＳＴシミュレーションが行われた。これらは、ＰＤＢ ＩＤ ｌｈｌｑ、ｌｈｌｒ、ｌｈｌｓ、ｌｏｉ９、ｌｏｉｕおよびｌｏｉｙを有するホロ複合体ごとに利用可能な、結晶構造を有する６つのみのリガンドであり^１６、これらのリガンドのための実験による結合親和性データは、同じ方法を用いる同じ公表文献を基にした。このため、リガンドの結合モードは結晶学から既知であるので、ここで報告される計算は、本質的に純粋にスコアリングの演習である。これらの結晶構造において、いずれのリガンドが結合しているかに関わらず、ＣＤＫ２受容体・コンフォメーションは本質的に同じである。このため、シミュレーション収束を大きく制限するタンパク質コンフォメーション変化のサンプリングは予期しなかった。計算において、５つ全ての摂動のための開始構造として、ｌｈｌｑからのタンパク質コンフォメーションおよびリガンドを用いることを選択した（図５において青い矢印で描画される）。なぜなら、組内の全ての他のリガンドは、リガンドｌｈｌｑのベンゼン環上に様々な置換基を付加することによって生成されるためでる。

ＦＥＰ／ＭＤ結果とＦＥＰ／ＲＥＳＴ結果との比較
結晶構造を有するリガンド間の相対結合自由エネルギー計算パスについてのＦＥＰ／ＭＤ結果およびＦＥＰ／ＲＥＳＴ結果が表１に与えられる。予測と、実験による結合親和性データとの相関が図６に表示される。実験による結合親和性の〜１４．６４ｋＪ／ｍｏｌ（〜３．５ｋｃａｌ／ｍｏｌ）の範囲に関わらず、ＦＥＰ／ＭＤは、５つの複合体（参照
複合体ｌｈｌｑを除外する）がほぼ同じ結合親和性を有し、結果として０．３２のＲ^２値が得られることを予測する。比較により、ＦＥＰ／ＲＥＳＴ結果は、実験値と高く相関し（Ｒ^２＝０．９１）、実験値からのＦＥＰ／ＲＥＳＴ結果の偏差は、全てのリガンド対について３．３４７ｋＪ／ｍｏｌ（０．８ｋｃａｌ／ｍｏｌ）未満である。

ＦＥＰ／ＭＤ結果およびＦＥＰ／ＲＥＳＴ結果の比較から、２つの方法からの複合体ｌｏｉ９、ｌｏｉｕおよびｌｉｏｙのための予測は、本質的に同じであり、実験結果との改善された相関の主な理由は、ＦＥＰ／ＲＥＳＴを用いたリガンドｌｈｌｓおよびｌｈｌｒについてのより良好な結合親和性推定値を基にしているようである。以下で、これらの２つの事例についてＦＥＰ／ＲＥＳＴがＦＥＰ／ＭＤよりも大幅に性能が優れている理由を検討する。

リガンドｌｈｌｒについてのＦＥＰ／ＭＤ軌道およびＦＥＰ／ＲＥＳＴ軌道の検査は、２−クロロ置換ベンゼン環がＦＥＰ／ＲＥＳＴ軌道において数回フリップする一方で、ＦＥＰ／ＭＤ軌道では環は初期コンフォメーションにトラップされたままであることを示した。結晶構造のコンフォメーションのうちの１つに対応する初期コンフォメーション、およびＦＥＰ／ＲＥＳＴ計算において観察される代替フリップコンフォメーションが図７に表示される。興味深いことに、リガンドｌｈｌｒについて、堆積されたＸ線結晶構造における電子密度は、ＦＥＰ／ＲＥＳＴ計算において見られる交互の結合モードに対応する代替リガンド結合モード（図７の右パネル）も示した。

ｌｈｌｓについて、結晶構造においてモデル化される結合モードに加えて、ＦＥＰ／ＲＥＳＴ軌道では別のコンフォメーションがサンプリングされるが、ＦＥＰ／ＭＤシミュレーションではサンプリングされない。図８は、結晶構造において見られる結合モードおよびＦＥＰ／ＲＥＳＴ計算においてサンプリングされる代替結合モードを表示する。結晶構造では、スルホンアミド基はＡｓｐ８６との２つの水素結合（ＨＢ）を行い、一方は主鎖との水素結合であり、もう一方は側鎖との水素結合である。ＦＥＰ／ＲＥＳＴシミュレーションにおいてサンプリングされる代替結合モードでは、スルホンアミド基は結晶構造に対しＣＳ結合の回りを約１２０°回転する。それに応じて、結合ポケットを取り囲む残基、特にＧｌｎ８５およびＡｓｐ８６が再配列し、スルホンアミド基とタンパク質残基との間に２つの水素結合が形成され、一方はＡｓｐ側鎖との水素結合であり、他方はＧｌｎ８５側鎖との水素結合である。Ｘ線結晶学における結合ポケットおよびリガンドを取り囲むこれらの数個の残基のＢ係数は、これらの基が非常に不安定であることを示す^１６。したがって、リガンドが溶媒内で２つの結合モード間で切り替わり得る可能性が非常に高い。ＦＥＰ／ＲＥＳＴシミュレーション軌道では、結晶学的に観察される結合モードと代替結合モードとの占有比はほぼ６０：４０であり、結晶学的に決定されたモードに有利である。一方、この代替結合モードのサンプリングは、結合親和性の正確な予測のために非常に重要なようであり、ＦＥＰ／ＭＤと比較して、この種の相対結合自由エネルギーを正確に算出する際に、ＦＥＰ／ＲＥＳＴのより高い有効性に対し責任を担う。興味深いことに、本発明は、実験データを参照することなくこれを決定することを可能にする。なぜなら、ＦＥＰ／ＭＤに関連する算出誤差は、本発明のこの特定の実施形態によって、高くなるように算出され（ｌｏｉｕおよびｌｈｌｓの場合、４．２２６ｋＪ／ｍｏｌ（１．０１ｋｃａｌ／ｍｏｌ）、ならびにｌｏｉｙ、ｌｏｉ９およびｌｈｌｓの場合、〜２．９２９ｋＪ／ｍｏｌ（〜０．７ｋｃａｌ／ｍｏｌ））、それに対し、ＦＥＰ／ＲＥＳＴに関連する算出誤差は、本発明によって低くなる（全てのリガンドについて０．８３６８ｋＪ／ｍｏｌ（０．２ｋｃａｌ／ｍｏｌ）〜１．６７４ｋＪ／ｍｏｌ（０．４ｋｃａｌ／ｍｏｌ））ように算出されるためである。

本発明の特定の実施形態を用いた相対自由エネルギー計算の一貫性および信頼性の評価
ＦＥＰ／ＲＥＳＴの一貫性および信頼性を評価するために、相対結合自由エネルギー計
算パス（図５）において３つのサイクルクロージャを形成する３つの更なるシミュレーションを実行した。これらのサイクルを形成する３つのパスのＦＥＰ／ＲＥＳＴ結果が表１の最後の３つの行に与えられる。本発明の特定の実施形態によって決定されたこれらのリガンドの結合親和性のための最終予測が表２に与えられる。３つ全てのサイクルについて、サイクルは、異なるパスから３．３４７ｋＪ／ｍｏｌ（０．８ｋｃａｌ／ｍｏｌ）未満の差異で自己矛盾なく閉じており、これは、自由エネルギーが良好に収束し、基礎をなす力場の制限内で予測の信頼性が高いことを示す。

一方、サイクルが十分に閉じており、様々なパスから計算される自由エネルギーがややタイトな閾値内で一致している場合であっても、これらは正確に同じではない。このため、同じリガンド対について、リガンド間のパスをどのようにトラバースするかに依拠して、自由エネルギー差のための複数の様々な推定値を有する。次に、自由エネルギーのための統計的に最良の推定値を取得するために、計算において生成される全ての情報をどのように最良に組み込むか？という問が残される。本発明の特定の実施形態を用いて、自由エネルギーの最適推定値、およびこれらの予測に関連する誤差を決定する。ここで、全てのサイクルは良好に閉じており、サイクルクロージャを有する最終予測は、サイクルクロージャを有しない上記の予測と良好に一致し、これらの計算の高い信頼性を更に立証する。最終予測と実験結果との間の偏差は、全て２．９２９ｋＪ／ｍｏｌ（０．７ｋｃａｌ／ｍｏｌ）未満であり、これらは実験データと非常に高い相関を有する（Ｒ^２＝０．９２）。

本発明のこの特定の実施形態が、ＦＥＰ／ＭＤ相対自由エネルギー計算を用いて算出された同様の閉じたサイクルに適用されたときと対象的に、リガンドｌｈｌｒおよびｌｈｌｓを伴う閉じたサイクルのヒステリシスは非常に大きく（それぞれ７．３２２ｋＪ／ｍｏｌ（１．７５ｋｃａｌ／ｍｏｌ）および５．６４８ｋＪ／ｍｏｌ（１．３５ｋｃａｌ／ｍｏｌ）の大きさ）、ＦＥＰ／ＭＤ計算が収束していないことを示す。本発明のこの特定の実施形態によるサイクルクロージャ自由エネルギー推定値を用いると、これらのＦＥＰ／ＭＤ自由エネルギー計算の誤差は、ＦＥＰ／ＲＥＳＴ結果（それぞれ１．４６４ｋＪ／ｍｏｌ（０．３５ｋｃａｌ／ｍｏｌ）および１．０４６ｋＪ／ｍｏｌ（０．２５ｋｃａｌ／ｍｏｌ）の誤差）よりもはるかに大きい（それぞれ４．２２６ｋＪ／ｍｏｌ（１．０１ｋａｌ／ｍｏｌ）および２．８４５ｋＪ／ｍｏｌ（０．６８ｋａｌ／ｍｏｌ））。

結晶構造を有しないリガンドにおける性能
本発明によるモデルを更に実証するために、ホロ複合体に対し利用可能な結晶構造を有しない同じ受容体に結合する別の１０個のリガンドについて更なるＦＥＰ／ＲＥＳＴ計算が実行された。これらの１０個のリガンドは、ハンドキャッスル（Ｈａｎｄｃａｓｔｌｅ）^１６によって同じ公表文献からランダムに選択された。上記で言及した結晶構造を有するリガンドと同じ方法を用いて結合親和性が測定され、これらは同じコアを保つ。元の公表文献^１６のように、リガンドの構造およびこれらの対応する複合体番号が表３に与えられる。これらのリガンドについてのＦＥＰ／ＲＥＳＴ予測相対結合親和性が、実験データと共に表４において報告される。サイクルクロージャは、これらのリガンドのためのパスにおいて構築され（図９）、サイクルクロージャのヒステリシスは全て４．１８４ｋＪ／ｍｏｌ（１．０ｋｃａｌ／ｍｏｌ）未満であり、これは、計算が高度に収束していることを示す。本発明のこの特定の実施形態による最終的な予測相対結合親和性が表５において報告される。相対結合自由エネルギー計算パスにおいて構築されるサイクルクロージャを有するリガンドについて、自由エネルギー予測に関連する誤差も表５において報告される。予測相対結合親和性と、リガンドの全体組のための実験データとの間の相関（結晶構造を有する６つのリガンドおよび結晶構造を有しない１０個のリガンド）が図１０に表示される。結晶構造を有しないこれらのリガンドについて、ＦＥＰ／ＲＥＳＴ自由エネルギー計算に適用されたときに本発明のこの特定の実施形態の性能が良好であることは明らかであり、手法の移行可能性を更に実証している。実験データと比較した予測結合親和性の平
均偏差は１．８４１ｋＪ／ｍｏｌ（０．４４ｋｃａｌ／ｍｏｌ）であり、Ｒ^２値は０．８２である。

結論
上記において、複数の実施形態において、本方法が、近年開発された自由エネルギー計算ＦＥＰ／ＲＥＳＴ法およびより定型的に用いられるＦＥＰ／ＭＤ自由エネルギー計算法に適用され、ＣＤＫ２−サイクリンＡ受容体に結合する一連の同種リガンド間の相対結合親和性を計算することができることを示した。本発明は、ＦＥＰ／ＲＥＳＴ自由エネルギー計算予測とＦＥＰ／ＭＤ自由エネルギー計算予測との比較を容易にした。結晶構造を有する６つのリガンドについて相対結合親和性を順位付けする際に、より定型的なＦＥＰ／ＭＤ法は性能が不良であったが、本発明は、ＦＥＰ／ＭＤ予測に関連する誤差が高いことを示す値を提供した。対照的に、ＦＥＰ／ＲＥＳＴにより計算された相対結合親和性は、実験データと極めて一致し、本発明の適用により収束することが示された。結晶構造を有しない１０個のリガンドを含むより大きなデータセットに対しＦＥＰ／ＲＥＳＴ法によって実行される相対自由計算への本発明の適用は、同種リガンドの相対結合親和性を正しく順位付けする本発明の能力を更に実証した。

本発明は、任意のタイプの相対自由エネルギー計算手順に関連するサイクルクロージャ計算から導出された情報を用いて自由エネルギーの一貫性および信頼性を評価する手段を提供する。本発明は、予測の一貫性および信頼性を判定する明確な判断基準、および良好に形成された誤差限界を有する自由エネルギーの最適推定値を提供する。ＣＤＫ２−サイクリンＡ系列への本発明の適用は、ＦＥＰ／ＲＥＳＴ相対自由エネルギー計算を用いる本発明による予測の信頼性が高いことを示し、これは、予測が実験データと極めて一致していることによっても証明される。同様に、本発明はＦＥＰ／ＭＤ予測の信頼性がより低いことを妥当に決定した。

本発明の一般性および有用性において、他の系および他の相対結合自由エネルギー計算への本発明の様々な実施形態の更なる適用を含む。
多岐にわたる相対結合自由エネルギー計算方式を用いた多岐にわたる系への本発明の様々な実施形態の更なる適用
Ｉ．ＪＮＫ１キナーゼへの本発明の一実施形態の適用
相対結合自由エネルギー計算パスにおける２つのサイクルクロージャを有するｊｎｋｌ受容体と結合するリガンドのサブセットについてＦＥＰ／ＲＥＳＴ相対結合自由エネルギーシミュレーションが実行された（図１１）。リガンドのサブセットのためのリガンド番号は［２，６，７，９，１０，１２，１３］であり、８つの相対結合自由エネルギー計算パスのための８つのＦＥＰ／ＲＥＳＴ結果が表６に列挙される。

サイクルクロージャの最良エネルギー推定量および誤差推定量を用いて、７つのリガンドについての結合親和性およびこれらの関連するエラーバーが表７において報告される。
サイクルクロージャ推定値は、シミュレーションの力場において系統誤差を補正することができないことに留意されたい。このため、リガンドのうちの２つ（リガンド１２および１３）のための予測結合親和性は、実験値と比較して６．２７６ｋＪ／ｍｏｌ（１．５ｋｃａｌ／ｍｏｌ）よりも大きく離れる。この偏差は、力場における誤差または実験測定値における不確実性のいずれかに起因する場合がある。しかし、予測結果と実験値との全体相関は非常に良好である。

ＩＩ．ＢＡＣＥへの本発明の一実施形態の適用
相対結合自由エネルギー計算パスにおける２つの独立したサイクルクロージャを有するＢＡＣＥに結合するリガンドのサブセットについてＦＥＰ／ＲＥＳＴシミュレーションが実行された（図１２に示す）。リガンドのサブセットのためのリガンド番号は［４ｊ，４
ｏ，４ｐ，１７ｄ，１７ｇ，１７ｈ］であり、７つの相対結合自由エネルギー計算パスのための７つのＦＥＰ／ＲＥＳＴ結果が表８に列挙される。

サイクルクロージャの最良自由エネルギー推定量および誤差推定量を用いて、６つのリガンドについての結合親和性およびこれらの関連するエラーバーが表９において報告される。

ここでもまた、リガンド１７ｈについての結合親和性は、約８．３６８ｋＪ／ｍｏｌ（２ｋｃａｌ／ｍｏｌ）だけ過剰予測されるが、予測結果と実験値と全体相関は非常に良好であることに留意されたい。

ＩＩＩ．シタロンデヒドラターゼに結合するリガンドのＦＥＰ／ＭＣ相対結合自由エネルギー計算への本発明の一実施形態の適用
活性部位における水結合を有するシタロンデヒドラターゼおよび有しないシタロンデヒドラターゼに結合する３つのリガンド間の相対結合親和性がＦＥＰ／ＭＣを用いて計算される^３９。４つの独立サイクルクロージャを有するこれらのシミュレーションのための相対結合自由エネルギー計算パスが図１３に表示される。相対結合自由エネルギー計算のための計算される結合親和性が表１０において報告される。

サイクル［Ｌ１（Ｈ２Ｏ），Ｌ２（Ｈ２Ｏ），Ｌ３（Ｈ２Ｏ）］のヒステリシスは、５．８５８ｋＪ／ｍｏｌ（１．４ｋｃａｌ／ｍｏｌ）よりもはるかに大きい１７．９９ｋＪ／ｍｏｌ（４．３ｋｃａｌ／ｍｏｌ）であり、これは、これらの計算が収束していないことを示し、このサイクルの内部の計算される結合親和性も実験データからの大きな偏差を有する。同様に、サイクル［Ｌ３（Ｈ２Ｏ），Ｌ３，Ｌ２，Ｌ２（Ｈ２Ｏ）］について、ヒステリシスは、５．８５８ｋＪ／ｍｏｌ（１．４ｋｃａｌ／ｍｏｌ）よりも大きい６．２７６ｋＪ／ｍｏｌ（１．５ｋｃａｌ／ｍｏｌ）であり、これは、これらの計算にも信頼性がないことを示す。

他方で、他のサイクル［Ｌ１（Ｈ２Ｏ），Ｌ１，Ｌ３，Ｌ３（Ｈ２Ｏ）］および［Ｌ１，Ｌ２，Ｌ３］のヒステリシスは非常に小さく（それぞれ、０．８３６８ｋＪ／ｍｏｌ（０．２ｋｃａｌ／ｍｏｌ）および２．０９２ｋＪ／ｍｏｌ（０．５ｋｃａｌ／ｍｏｌ））、計算に信頼性があることを示す。上記に記載のサイクルクロージャ推定値を有する３つのリガンド間の最終的な予測相対結合親和性が表１１において報告される。

ＩＶ．モデルエンクロージャーのための試験粒子挿入による相対結合自由エネルギーサンプリングへの本発明の一実施形態の適用
融合原子メタンと電荷を有するモデルエンクロージャーおよび電荷を有しないモデルエンクロージャーとの間の相対結合自由エネルギーが試験粒子挿入法を用いて計算される^４０。結合メタンを有するモデルエンクロージャーおよび結合メタンを有しないモデルエンクロージャーに帯電する自由エネルギーも、ＦＥＰ／ＭＤサンプリングを用いて計算される。相対結合自由エネルギー計算パスは、閉じたサイクルを形成し（図１４）、サイクル上のエッジごとの関連自由エネルギーがここで報告される。

上記で記載したサイクルクロージャ解析を用いて、図１６において描く４つの状態間の予測自由エネルギー差およびこれらの関連誤差が表１２において報告される。
複数の実施形態が記載された。それにもかかわらず、本発明の趣旨および範囲から逸脱することなく種々の変更を行うことができることが理解されよう。

表
表１：参照リガンドｌｈｌｑと比較した、結晶構造を有する５つのリガンドについてＦ
ＥＰ／ＭＤおよびＦＥＰ／ＲＥＳＴにより予測される相対結合自由エネルギー。結晶構造を有する６つのリガンドについて３つのサイクルクロージャを形成する３つの相対結合自由エネルギー計算パスの結果も表に含まれる。ＢＡＲ解析誤差推定を用いて計算される自由エネルギーの誤差も表に含まれる。ＦＥＰ／ＭＤによる予測結果は、リガンドｌｈｌｒおよびｌｈｌｓについて実験データからの大きな偏差を有し（表における太字）、一方、ＦＥＰ／ＲＥＳＴによる予測結果は、全て実験データと良好に一致する。自由エネルギーはｋＪ／ｍｏｌ（ｋｃａｌ／ｍｏｌ）単位で報告される。

表２：本発明の実施形態によって導出される式を用いるＦＥＰ／ＭＤ計算およびＦＥＰ／ＲＥＳＴ計算からの、結晶構造を有する６つのリガンドについての最終自由エネルギー予測。サイクルクロージャを用いる最終自由エネルギー予測についての、本発明の実施形態によって導出される式を用いる誤差推定値も表に含まれる。リガンドＩｈｌｑが参照として用いられ、このため、リガンドＩｈｌｑの偏差は設計により０である。リガンドＩｈｌｑがこの研究のために用いられる参照リガンドであるので、リガンドＩｈｌｑについてサイクルクロージャエラーバーは報告されず、リガンドＩｈｌｑの結合自由エネルギーは、プロット目的でリガンドＩｈｌｑの実験値に全く同一に設定されることに留意されたい。自由エネルギーはｋＪ／ｍｏｌ（ｋｃａｌ／ｍｏｌ）単位で報告される。

表３：ＦＥＰ／ＲＥＳＴを用いて研究されるＣＤＫ２阻害剤およびこれらの化合物＃の構造。

表４：結晶構造を有しない追加の１０個の化合物についての、ＦＥＰ／ＲＥＳＴにより予測される相対結合自由エネルギー。ＢＡＲ解析誤差推定を用いて計算される自由エネルギーについての誤差も表において報告される。自由エネルギーはｋＪ／ｍｏｌ（ｋｃａｌ
／ｍｏｌ）単位で報告される。

表５：本発明の実施形態によって導出される式を用いる、リガンドの全体組についての最終自由エネルギー予測。サイクルクロージャを用いる最終自由エネルギー予測についての、本発明の実施形態を用いる誤差推定値も表に含まれる。リガンドＩｈｌｑが参照として用いられ、このため、リガンドＩｈｌｑの偏差は設計により０である。リガンドＩｈｌｑがこの研究のために用いられる参照リガンドであるので、リガンドＩｈｌｑについてサイクルクロージャエラーバーは報告されず、リガンドＩｈｌｑの結合自由エネルギーは、プロット目的でリガンドＩｈｌｑの実験値に全く同一に設定されることに留意されたい。リガンド２０、３１および３２についてもサイクルクロージャエラーバーは報告されない。なぜならこれらの種は、誤差について閉じたサイクルを形成することが要求されるように、参照リガンド以外の他のリガンドに連結することが技術的に困難であったためである。自由エネルギーはｋＪ／ｍｏｌ（ｋｃａｌ／ｍｏｌ）単位で報告される。

表６：ｊｎｋ１リガンドの８つの相対結合自由エネルギー計算パスについてのＦＥＰ／ＲＥＳＴ結果。

表７：ｊｎｋ１に結合している７つのリガンドについての予測結合親和性およびそれらの関連誤差。

表８：ＢＡＣＥリガンドのための７つの相対結合自由エネルギー計算パスについてのＦＥＰ／ＲＥＳＴ結果。

表９：ＢＡＣＥに結合している６つのリガンドについての予測結合親和性およびそれらの関連誤差。

表１０：シタロンデヒドラターゼリガンドの９個の相対結合自由エネルギー計算パスについてのＭＣ／ＦＥＰ結果。

表１１：シタロンデヒドラターゼに結合している３つのリガンドについての予測結合親和性およびそれらの関連誤差。３つのリガンドについての予測結合親和性は比較的大きな誤差を有するにもかかわらず、非常に大きな誤差を有する相対結合自由エネルギー計算パスはフィルタリング除去され、正しい順位が予測される。

＊注：表において報告される結合親和性はリガンド１と比較した相対結合親和性である。
表１２：図１４において描かれる４つの状態間のサイクルクロージャにより予測される自由エネルギー差およびそれらの関連誤差。Ｍは融合原子メタンを表す、Ｐは疎水性プレートを表す、ＣＰは帯電プレートを表す、ＭＰは融合原子メタンと結合した疎水性プレートを表す、ＭＣＰは融合原子メタンと結合した帯電プレートを表す。

Claims (17)

1. 個々のリガンド組メンバーと受容体との間で複合体を形成するための、前記受容体と１組のリガンドの個々のメンバーとの間の結合の相対強度を決定するコンピュータ実装された方法であって、
   ａ．コンピュータが、少なくとも１つの閉じたサイクルを形成するリガンド対の組について、前記閉じた熱力学的サイクルについて観察される前記ヒステリシスの大きさを確率的に導く、前記閉じた熱力学的サイクルの区画の各々に沿った結合自由エネルギー差と前記結合自由エネルギー差に関する誤差分布とを決定する工程と、
   ｂ．前記コンピュータが、前記ａの工程において決定された確率モデルに含まれる前記閉じた熱力学的サイクルにおける区画ごとの最も確からしい結合自由エネルギー差を決定する工程と、
   ｃ．前記コンピュータが、前記ｂの工程における前記閉じた熱力学的サイクルにおける各区画に沿ったリガンド対ごとの前記最も確からしい結合自由エネルギー差に関連する最も確からしい誤差を決定する工程と、
   ｄ．前記コンピュータが、コンピュータユーザに対する、前記ｃの工程における前記結合自由エネルギー差および前記誤差を表す値を出力する工程と、を備える、方法であって、同方法はさらに、
   前記コンピュータが、前記結合自由エネルギー差の計算のための自由エネルギー、ならびにこれらの対応する初期リガンドおよび最終リガンドのリガンド対を列挙する工程と、
   前記コンピュータが、閉じたサイクルを形成するリガンド対を探索して発見し、同発見したリガンド対を表す連結性行列を生成する工程と、
   前記コンピュータが、全体グラフを独立した部分グラフに分解する工程であって、前記独立した部分グラフの各々は、他の独立した部分グラフと連結する閉じたサイクルを有していない、工程と、
   前記コンピュータが、前記独立した部分グラフの各々に対応する連結性行列を入手する工程と、
   前記コンピュータが、前記部分グラフの各々に対応するそれぞれの連結性行列を、独立した行ベクトルのみを含む独立した行列に縮約する工程と、
   前記コンピュータが、前記部分グラフごとに、前記リガンド対について最も確からしい結合自由エネルギー推定値を生成する工程と、
   前記コンピュータが、前記最も確からしい結合自由エネルギー推定値に関連する戦記最も確からしい誤差を決定する工程と
   を備える方法。
2. 前記誤差を決定する工程は、前記コンピュータが、２つ以上の閉じた熱力学的サイクルの区画に沿ったリガンド間での結合自由エネルギー差を解析する工程と、前記コンピュータが、閉じた熱力学的サイクルの各々に関するヒステリシスの大きさを決定する工程とを含む、請求項１に記載の方法。
3. 前記ｃおよびｄの工程は、前記コンピュータが、関数
   を最小にする区画ごとの１組の自由エネルギー値を、制約
   を用いて決定することを含み、
   Ｅｉは、所与の区画ｉについての計算される自由エネルギー差であり、
   Ｆ_ｉは、所与の区画ｉについての理論上の自由エネルギー差であり、
   σ_ｉは、区画ｉについての前記計算される自由エネルギー差の標準偏差であり、全ての閉じたサイクルについての前記理論上の自由エネルギー差の和は０である、請求項１に記載の方法。
4. 前記受容体はタンパク質である、請求項１に記載の方法。
5. 前記リガンドは同種である、請求項１に記載の方法。
6. 工程１ａ．前記ａの工程における前記観察されるヒステリシスの前記確率モデルの構築においてガウス分布が仮定される、請求項１に記載の方法。
7. 自由エネルギーシミュレーションに関連する前記誤差分布は、工程１ａにおいて一様であると仮定される、請求項１に記載の方法。
8. 自由エネルギーシミュレーションに関連する前記誤差分布は、工程１ａにおいてベネット誤差に加算されると仮定される、請求項１に記載の方法。
9. 前記閉じた熱力学的サイクルの連結性はグラフとして表される、請求項１に記載の方法。
10. 前記閉じた熱力学的サイクルの連結性は行列として表される、請求項１に記載の方法。
11. 前記確率的決定はグラフ理論法を実行することを含む、請求項１に記載の方法。
12. 前記ａの工程における前記決定は行列代数法を実行することを含む、請求項１に記載の方法。
13. 前記ａの工程における前記決定はベイズ法を実行することを含む、請求項１に記載の方法。
14. 前記ａの工程における前記決定は最大尤度法を実行することを含む、請求項１に記載の方法。
15. 請求項１に記載の方法を実行するための有形の非一時的な命令を含むコンピュータ可読記憶媒体。
16. コンピュータと、請求項１に記載の方法を実行するための非一時的なコンピュータ可読命令を記憶するコンピュータ可読記憶媒体とを備えるコンピュータシステムであって、前記コンピュータは、前記媒体から前記命令を読み出し前記命令を実行するように構成されている、コンピュータシステム。
17. プロセッサと、請求項１に記載の方法を実行するための非一時的なコンピュータ可読命令を記憶するメモリとを備える汎用グラフィック処理ユニットであって、前記プロセッサは、前記メモリから前記命令を読み出し前記命令を実行するように構成されている、汎用グラフィック処理ユニット。