(57)【要約】分子のコンピュータモデリングのために、換算座標を用いたモデルが、このモデルの運動方程式を積分する十分に安定な陰的積分法と共に使用される。この積分法における時間刻みは、１００を超える範囲で変動し得、コンピュータの効率を著しく増大させ、かつ、計算結果を促進する。静的解析シミュレーションおよび分動力学シミュレーションの両方が、手近な適用である。好適な換算座標分子モデルは、デカルト全原子座標ではなく、ねじれ角座標を組み込んだ、剛体分割法である。好適な十分に安定な積分方法は、誤差制御された動力学的な計算のためのＬ安定性一段階Ｒａｄａｕ５法、およびエネルギー極小化（静的）計算のためのＬ安定性陰的Ｅｕｌｅｒ法を含む。厳密性が少ない安定性必要条件を有する用途に対して、高安定性および効率的な他段階ＤＡＳＳＬ法が好ましい。【特許請求の範囲】
【請求項１】
分子の挙動をモデリングする方法であって、該方法は、以下：
該分子についてモデルを選択する工程であって、該モデルは、該分子についての運動方程式を有する、工程；および
大きな時間刻みでＬ安定性陰的積分器を用いて該モデルの方程式を積分して、該分子の該挙動の計算を得る工程
を包含する、方法。
【請求項２】
前記大きな時間刻みが少なくとも２００フェムト秒の区間を含む、請求項１に記載の方法。
【請求項３】
前記積分する工程が変動する時間刻みで実行される、請求項２に記載の方法。
【請求項４】
前記積分する工程における誤差を補正し、経時的な前記分子の状態履歴を得る工程をさらに包含する、請求項１に記載の方法。
【請求項５】
前記選択する工程が、スティッフ系モデルを選択して、経時的に前記分子の状態履歴を得る工程を含む、請求項１に記載の方法。
【請求項６】
前記積分する工程がエネルギー保存を回避し、前記分子の極小エネルギー状態を得る、請求項１に記載の方法。
【請求項７】
請求項１に記載の方法であって、前記Ｌ安定性積分器が、陰的Ｅｕｌｅｒ法、Ｒａｄａｕ５法、ＳＤＩＲＫ３法、ＳＤＩＲＫ４法および他の陰的Ｒｕｎｇｅ−Ｋｕｔｔａ法からなる群より選択される、積分器を含む、方法。
【請求項８】
請求項３に記載の方法であって、該方法は、
前記積分する工程における誤差を補正し、経時的に前記分子の状態履歴を得る工程を、さらに包含する、方法。
【請求項９】
前記選択する工程は、スティッフ系モデルを選択して、経時的に前記分子の状態履歴を得る工程を包含する、請求項３に記載の方法。
【請求項１０】
前記積分する工程がエネルギー保存を回避し、前記分子について極小エネルギー状態を得る、請求項３に記載の方法。
【請求項１１】
前記モデルが、前記分子の前記挙動の計算を速めるように選択された内部座標で記述される、請求項１に記載の方法。
【請求項１２】
前記モデルが、前記分子の、ねじれ角、剛体モデルを含む、請求項１１に記載の方法。
【請求項１３】
分子の挙動をモデリングする方法であって、該方法は、以下：
該分子についてのモデルを選択する工程であって、該モデルが運動方程式を有する工程；および
Ｌ安定性積分器を選択する工程；
少なくとも１００の範囲にわたって変動する間隔の時間刻みにおいて、前記Ｌ安定性積分器で該モデル方程式を積分して、該分子の挙動の計算を得る工程
を包含する、方法。
【請求項１４】
前記時間刻みが少なくとも２００フェムト秒の間隔を含む、請求項１３に記載の方法。
【請求項１５】
前記Ｌ安定性積分器を、該分子からエネルギーを除去するように選択し；そして、前記モデル方程式を、エネルギー保存なしで積分して、前記分子の極小エネルギー状態を得る、請求項１４に記載の方法。
【請求項１６】
前記Ｌ安定性積分器が陰的Ｅｕｌｅｒ積分器を含む、請求項１５に記載の方法。
【請求項１７】
前記モデル方程式を、誤差補正を行って積分し、経時的に前記分子の状態履歴を得る、請求項１４に記載の方法。
【請求項１８】
前記モデルをスティッフ運動方程式について選択し、経時的に前記分子の状態履歴を得る、請求項１４に記載の方法。
【請求項１９】
前記モデルをスティッフ運動方程式について選択し、そして前記モデル方程式を、誤差補正を行って積分し、経時的に前記分子の状態履歴を得る、請求項１４に記載の方法。
【請求項２０】
前記Ｌ安定性積分器がＲａｄａｕ５積分器を含む、請求項１９に記載の方法。
【請求項２１】
前記Ｌ安定性積分器が、陰的Ｅｕｌｅｒ法、Ｒａｄａｕ５法、ＳＤＩＲＫ３法、ＳＤＩＲＫ４法、および陰的Ｒｕｎｇｅ−Ｋｕｔｔａ法からなる群より選択される、請求項１４に記載の方法。
【請求項２２】
前記モデルが、前記分子の前記挙動の計算を速めるように選択された内部座標で記述される、請求項１４に記載の方法。
【請求項２３】
前記モデルが、前記分子の、ねじれ角、剛体モデルを含む、請求項２２に記載の方法。
【請求項２４】
複数の第２の分子を伴う第１の分子の挙動をモデリングする方法であって、該方法は、以下：
該第１の分子についての第１のモデルを選択する工程であって、該モデルが該第１の分子についての運動方程式を有する工程；
該第２の分子の各々についての第２のモデルを選択する工程であって、該モデルが該第２の分子についての運動方程式を有する工程；
Ｌ安定性積分器を選択する工程；
少なくとも１００の範囲にわたって変化する間隔の時間刻みにおいて、前記Ｌ安定性積分器で該モデル方程式を積分して、該複数の第２の分子を伴う第１の分子の挙動の計算を得る工程
を包含する、方法。
【請求項２５】
前記モデル方程式が、前記挙動の計算を速めるように選択された内部座標で記述される、請求項２４に記載の方法。
【請求項２６】
前記第２の分子が、塩、溶媒ならびに他の有機化合物および無機化合物からなる群より選択される、請求項２４に記載の方法。
【請求項２７】
前記第２の分子が水を含む、請求項２６に記載の方法。
【請求項２８】
前記第１の分子がタンパク質である、請求項２５に記載の方法。
【請求項２９】
前記大きな時間刻みが、少なくとも２００フェムト秒の区間である、請求項２５に記載の方法。
【請求項３０】
前記Ｌ安定性積分器を、該分子からエネルギーを除去するように選択し；そして前記モデル方程式を、エネルギー保存なしで積分して、前記分子の極小エネルギー状態を得る、請求項２９に記載の方法。
【請求項３１】
前記Ｌ安定性積分器が陰的Ｅｕｌｅｒ積分器を含む、請求項３０に記載の方法。
【請求項３２】
前記モデル方程式を、誤差補正を行って積分して、経時的に前記分子の状態履歴を得る、請求項２５に記載の方法。
【請求項３３】
前記モデルをスティッフ運動方程式について選択し、経時的な前記分子の状態履歴を得る、請求項２５に記載の方法。
【請求項３４】
前記モデルをスティッフ運動方程式について選択し、そして前記モデル方程式を、誤差補正を行って積分して、経時的な前記分子の状態履歴を得る、請求項２５に記載の方法。
【請求項３５】
コンピュータ上にて分子の挙動をモデリングするためのコンピュータコードであって、該コードは、以下：
該分子ついてモデルを規定するための第１のモジュールであって、該モデルは、該分子についての運動方程式を含む、モジュール、および
Ｌ安定性陰的積分器を用いて該モデルの方程式を積分して、該分子の該挙動の計算を得る第２のモジュール、
を備える、コード。
【請求項３６】
前記第２のモジュールが変動する時間刻みで前記運動方程式を積分する、請求項３５に記載のコンピュータコード。
【請求項３７】
前記時間刻みが少なくとも１００の範囲にわたる大きさで変化する、請求項３６に記載のコンピュータコード。
【請求項３８】
前記第１のモジュールが、内部座標を用いて前記モデルを定義する、請求項３５に記載のコンピュータコード。
【請求項３９】
前記内部座標が一般化座標および一般化速度を含む、請求項３８に記載のコンピュータコード。
【請求項４０】
前記第１のモジュールが、前記分子についての、剛体多体系、ねじれ角モデルを規定する、請求項３９に記載のコンピュータコード。
【請求項４１】
標的との相互作用について化合物ライブラリーをスクリーニングする方法であって、該方法は、以下：
（ａ）標的と化合物の相互作用についてモデルを選択する工程であって、該モデルは、該化合物および該標的についての運動方程式を有する、工程；
（ｂ）第１の化合物ライブラリーについてのデータを該運動方程式に入力する工程；
（ｃ）Ｌ安定性積分器を用いて、大きな時間刻みで該モデル方程式を積分し、該標的および該化合物の運動の計算を得て、それによって該標的と該化合物の相互作用を得る工程；
（ｄ）ライブラリーにおける、各化合物について、工程（ｂ）および（ｃ）を反復する工程；
（ｅ）該標的と該化合物の相互作用を比較する工程；
（ｆ）該標的とのその相互作用に基づいて選択された化合物を合成する工程；
を包含する、方法。
【請求項４２】
前記化合物ライブラリーが、前記標的と相互作用することが既知であるリード化合物および該標的との相互作用について試験される試験化合物を含む、請求項４１に記載の方法。
【請求項４３】
前記リード化合物がポリペプチドであり、そして前記試験化合物は低分子である、請求項４２に記載の方法。
【請求項４４】
前記リード化合物が抗体である、請求項４３に記載の方法。
【請求項４５】
請求項４２に記載の方法であって、前記化合物のうちの１つが、該標的と相互作用することが既知であるリード化合物であり、かつ前記比較する工程によって、前記試験化合物と前記標的との間の相互作用を、該標的と該リード化合物の相互作用と比較して、該リード化合物の相互作用と類似の該標的との相互作用を有する試験化合物を選択する、方法。
【請求項４６】
請求項４２に記載の方法であって、該方法は、前記標的と前記リード化合物を接触させ、そして該リード化合物と該標的との間の相互作用を検出することによって、１次ライブラリーから該リード化合物を同定する工程をさらに包含する、方法。
【請求項４７】
請求項４１に記載の方法であって、ここで、前記工程（ｂ）および（ｃ）の異なった反復が、第１の化合物および第２の化合物において実施され、該第２の化合物が、前記標的と第１の化合物の相互作用に基づいて選択される、方法。
【請求項４８】
前記標的との相互作用について、前記合成された化合物を試験する工程をさらに含む、請求項４１に記載の方法。
【請求項４９】
前記試験する工程が、非ヒト動物およびヒトにおいて、インビトロで実施される、前記４８に記載の方法。
【請求項５０】
薬学的組成物として前記合成された化合物を処方する工程をさらに包含する、請求項４１に記載の方法。
【請求項５１】
化合物および／または前記標的のライブラリーの少なくとも１つの構造に関連するデータを決定する工程をさらに包含する、請求項４１の方法。
【請求項５２】
前記データが、Ｘ線結晶解析によって決定される、請求項５１に記載の方法。
【請求項５３】
前記データが、赤外線分光法もしくは紫外線分光法、またはＮＭＲによって決定される、請求項５１に記載の方法。
【請求項５４】
請求項４１に記載の方法であって、前記化合物がタンパク質、核酸、多糖類、リン脂質、ホルモン、プロスタグランジン、ステロイド、および低分子からなる群より選択される、方法。
【請求項５５】
請求項５４に記載の方法であって、ここで、前記化合物が、βターン模倣物、芳香族性化合物、ヘテロ環式化合物、ベンゾジアゼピン、オリゴマー性Ｎ置換グリシン、およびオリゴカルバメートからなる群より選択される低分子である、方法。
【請求項５６】
請求項４１に記載の方法であって、ここで、前記標的がタンパク質、核酸、炭水化物および脂質からなる群より選択される、方法。
【請求項５７】
前記標的がレセプターである、請求項５６に記載の方法。
【請求項５８】
前記標的が膜結合レセプターである、請求項５７に記載の方法。
【請求項５９】
請求項４１に記載の方法であって、前記標的および／または該標的と相互作用する化合物を含む溶媒または基質についてのデータを、前記運動方程式に入力する工程をさらに包含する、方法。
【請求項６０】
前記基質がリン脂質膜である、請求項５９に記載の方法。
【請求項６１】
前記溶媒が水性溶媒である、請求項４１に記載の方法。
【請求項６２】
前記溶媒が有機溶媒である、請求項４１に記載の方法。
【請求項６３】
前記データが、前記化合物の成分の同一性を含む、請求項４１に記載の方法。
【請求項６４】
前記データが、前記化合物の原子の同一性を含む、請求項６３に記載の方法。
【請求項６５】
前記データが、Ｘ線結晶解析データを含む、請求項４１に記載の方法。
【請求項６６】
前記データが、前記運動方程式に、環境要因を入力する工程をさらに含む、請求項４１に記載の方法。
【請求項６７】
前記環境要因が、前記化合物と標的との間の相互作用が決定される温度または圧力である、請求項４１に記載の方法。
【請求項６８】
前記化合物ライブラリーが少なくとも１０１０のメンバーを含む、請求項４１に記載の方法。
【請求項６９】
前記化合物ライブラリーが少なくとも１０５０のメンバーを含む、請求項４１に記載の方法。
【請求項７０】
請求項４１に記載の方法であって、前記積分する工程が前記化合物と前記標的との間の結合親和性を決定し、そして、前記比較する工程が該標的との異なる化合物の結合親和性を比較し、そして前記合成する工程は、該標的について最も高い親和性を有する化合物を合成する、方法。
【請求項７１】
請求項４１に記載の方法であって、前記積分する工程が前記化合物と前記標的との間の相互作用を決定し、ここで、該相互作用は、少なくとも１０９Ｍ−１の親和性で該化合物が該標的に結合する、方法。
【請求項７２】
請求項４１に記載の方法であって、前記積分する工程が前記化合物と前記標的との間の相互作用を決定し、ここで、該相互作用は、該化合物が該標的を介してシグナルを変換する、方法。
【請求項７３】
請求項４１に記載の方法であって、ここで、前記化合物が潜在的な界面活性剤であって、そして、前記積分する工程が、該化合物と前記標的との間の相互作用を決定し、該相互作用が、該化合物が該標的を変性することを示す、工程。
【請求項７４】
所望の機能特性を有するようにタンパク質を進化させる方法であって、該方法は、以下：
（ａ）タンパク質の参照形態についてモデルを選択する工程であって、該モデルは、該タンパク質についての運動方程式を有する、工程；
（ｂ）タンパク質のアミノ酸置換についてのデータを、該運動方程式に入力する工程；
（ｃ）Ｌ安定性積分器を用いて、大きな時間刻みで該モデル方程式を積分し、該アミノ酸置換を有するタンパク質の運動の計算を得る工程；
（ｄ）さらなるアミノ酸置換について、工程（ｂ）および（ｃ）を反復する工程；
（ｅ）異なるアミノ酸置換を有するタンパク質の運動を比較する工程；
（ｆ）該比較に基づいて選択されたアミノ酸を有するタンパク質を合成する工程；
を包含する、方法。
【請求項７５】
請求項７４に記載の方法であって、ここで、所望する機能特性について、選択された該合成されたタンパク質を試験する工程を、さらに包含する、方法。
【請求項７６】
請求項７４に記載の方法であって、ここで、前記所望の機能特性が標的に結合する能力である、方法。
【請求項７７】
前記所望の特性は酵素活性である、請求項７４に記載の方法。
【請求項７８】
免疫グロブリン鎖をヒト化する方法であって、該方法は、以下：
（ａ）マウス抗体由来のＣＤＲ領域およびヒト抗体由来の可変領域フレームワークを含む免疫グロブリン鎖についてのアミノ酸配列を提供する工程：
（ｂ）該免疫グロブリン鎖についてモデルを選択する工程であって、該モデルは該免疫グロブリン鎖についての運動方程式を有する、工程；
（ｃ）Ｌ安定性積分器を用いて、大きな時間刻みで該モデル方程式を積分し、該免疫グロブリン鎖の運動の計算を得る、工程；
（ｄ）該可変フレームワーク領域にある、どのアミノ酸残基が該ＣＤＲ領域と相互作用するかを、該モデルから決定する工程；
（ｅ）該マウス抗体由来の対応するアミノ酸で、該ＣＤＲ領域と相互作用する該可変フレームワーク領域におけるアミノ酸残基のうちの１以上を、置換する工程；
（ｆ）１つ以上のアミノ酸残基を含む免疫グロブリンを合成する工程
を包含する、方法。
【請求項７９】
標的との結合について、前記合成された免疫グロブリン鎖を試験する工程をさらに含む、請求項７８に記載の方法。
【請求項８０】
特定の環境において、１以上の分子の挙動および特性を計算するための方法であって、該方法は、以下：
（ａ）該分子をおよびそれらの環境を数学的にモデリングする工程であって、該モデルは、換算された座標のセットによって表現される該分子についての運動方程式を有する、工程；および
（ｂ）大きな時間刻みを用いる陰的積分器を使用して、該モデル方程式を積分する工程であって、該積分器は、安定性特性を有して、そして刻み幅選択法が、該環境について十分な精度で該挙動または該特性を計算するとき、該大きな時間刻みの使用を可能にする、工程
を包含する、方法。
【請求項８１】
前記大きな時間刻みが、少なくとも２００フェムト秒の区間を含む、請求項８０に記載の方法。
【請求項８２】
前記積分する工程が、変動する時間刻みで実行される、請求項８０に記載の方法。
【請求項８３】
前記変動する時間刻みが、少なくとも２００フェムト秒の区間を含む、請求項８２に記載の方法。
【請求項８４】
前記刻み幅選択法が、精度推定を包含する請求項８０に記載の方法。
【請求項８５】
前記刻み幅選択法が、収束の必要条件を包含する、請求項８０に記載の方法。
【請求項８６】
前記刻み幅選択法が、エネルギー散逸の必要条件を包含する、請求項８０に記載の方法。
【請求項８７】
前記積分器がＬ安定性特性を有する、請求項８０に記載の方法。
【請求項８８】
請求項８０に記載の方法であって、前記積分器は、Ｒａｄａｕ、ＳＤＩＲＫ、ＳＩＲＫまたはＲｏｓｅｎｂｒｏｋのファミリーの積分法の２次以上のＬ安定性メンバーからなる群より選択される積分器を含む、方法。
【請求項８９】
Ｌ安定性積分器がＲａｄａｕ５積分法である、請求項８７に記載の方法。
【請求項９０】
前記積分器が、ＤＡＳＳＬおよびスティッフ系または微分代数系について設計された他の陰的多段法からなる群より選択される積分器を含む、請求項８０に記載の方法。
【請求項９１】
前記座標が、２以上の原子の各々を含む１以上の剛体および内部座標の使用によって換算される、請求項８０に記載の方法。
【請求項９２】
前記内部座標が、ねじれ角である、請求項９１に記載の方法。
【請求項９３】
前記座標が、剛体サブコンポーネントまたは伸縮サブコンポーネントへと下部構造分解することによって換算される、請求項８０に記載の方法。
【請求項９４】
前記環境が真空状態である、請求項８０に記載の方法。
【請求項９５】
前記環境が溶媒を含む、請求項８０に記載の方法。
【請求項９６】
前記溶媒が非明示的な表現物である、請求項９５に記載の方法。
【請求項９７】
前記非明示的な溶媒が、膜領域のような不均一な溶媒特性を含む、請求項９６に記載の方法。
【請求項９８】
前記環境が、動力学シミュレーションを含む、請求項８０に記載の方法。
【請求項９９】
前記環境がニュートン力学を含む、請求項９８に記載の方法。
【請求項１００】
前記環境が、Ｌａｎｇｅｖｉｎ動力学を含む、請求項９８に記載の方法。
【請求項１０１】
前記環境が、前記分子の換算エネルギー状態についての検索を含む、請求項８０に記載の方法。
【請求項１０２】
前記検索が、初期配置の局所エネルギーベイスンのみを含む、請求項１０１に記載の方法。
【請求項１０３】
前記検索が、初期配置の局所エネルギー以外のエネルギーベイスンを含む、請求項１０１に記載の方法。
【請求項１０４】
前記分子が、非ネイティブ環境にある単一の生体高分子を含み、そして前記特性が、該生体高分子の折り畳まれたネイティブ構造を含む、請求項８０に記載の方法。
【請求項１０５】
前記生体高分子が、ポリペプチドまたはタンパク質である、請求項１０４に記載の方法。
【請求項１０６】
前記生体高分子が核酸である、請求項１０４に記載の方法。
【請求項１０７】
請求項８０に記載の方法であって、前記分子が標的分子およびリガンド分子を含み、ここで、前記挙動は、標的に対するリガンドの結合を含むか、または、前記特性は、結合親和性、結合選択性、結合速度または他の結合特性を含む、方法。
【発明の詳細な説明】
【０００１】
（関連出願への相互参照）
本出願は、仮特許出願第６０／２４５，６８８号（２０００年１１月２日に出願）、そしてさらに第６０／２４５，７３０号（２０００年１１月２日に出願）、第６０／２４５，７３１号（２０００年１１月２日に出願された）、および第６０／２４５，７３４号（２０００年１１月２日に出願された）の優先出願日の利益が与えられる。これらの出願の全てを本明細書中で参考として援用する。
【０００２】
（発明の背景）
本発明は分子モデリングの分野に関し、より詳細には、分子の挙動および性質あるいは溶液中にある分子の相互作用系の予測のためのコンピュータにインプリメントされた方法に関する。本発明は、これらの所望された予測を行うための、分子力学モデルおよび時間積分を利用する計算法に関する。
【０００３】
分子力学における物体（ｂｏｄｙ）の運動は、ニュートンの運動法則によって決定される。力Ｆを受ける、質量ｍの物体について、ニュートンの第２法則は、以下：
Ｆ＝ｍａ
のように記述される。すなわち、物体の加速度ａは、物体に及ぼされる力全体に比例する。この単純な方程式には、高分子の動的モデリングについて膨大な複雑さが隠されている。この物体の加速度は、物体の速度の時間微分であり、物体の速度を決定するためには、その加速度が時間に関して積分されなければならない。同様にして、物体の速度は、物体の位置の時間微分であり、物体の位置を決定するためには、その速度が時間に関して積分されなければならない。従って、物体に及ぼす力の知識によると、積分演算が、所与の時間における物体の速度および位置を決定するために行われなければならない。
【０００４】
１つの分子内に、運動が考慮されなければならない複数の物体が存在する。代表的な分子力学モデルにおいて、分子の各原子が１つの物体として考えられ、各原子は、複数かつ複雑な力を受け、この力は、この系内のいずれの分子における他の原子のそれぞれの、現在の位置、ならびに、環境または溶媒の影響に、潜在的に関係する。従って、この分子の運動および形状の計算は、この系内のそれぞれの原子の位置および運動の決定を必要とする。従って、数千の原子を有する複雑な分子を含む分子の構造、動力学、および熱力学の計算は、コンピュータに十分に適合した課題のように考えられている。
【０００５】
実際には、分子モデリングの領域は、首尾よく運動をシミュレート（分子動力学、すなわち（ＭＤ））し、コンピュータによって多くの複雑な分子系のエネルギー極小状態または静止状態（静的解析）を決定してきた。代表的な分子モデリング用途は、酵素−リガンドドッキング、分子拡散、反応経路、相転移、およびタンパク質のフォールディングの研究を含んでいる。生命科学、ならびに、製薬、ポリマー、および化学産業の研究者が、複雑な分子の化学的過程の性質を理解し、それに従って新規の薬剤および材料を設計するために、それらの技術を使用し始めている。本来、これらのツールの受容は、現実のものを表現するときの結果の精度、モデルリングされるべき分子系のサイズおよび複雑さ、および解が得られる速度を含む、いくつかの要因に基づく。多くの計算の精度は実験と比較され、一般的に特定の境界内では適切であると理解される。しかし、従来技術におけるこれらのツールの使用は、適度の大きさの分子または分子系でさえも、それらをモデリングして有用である十分な長さの分子の時間履歴を得るために膨大な計算パワーを、必要としてきた。
【０００６】
時間積分を含む分子モデリング課題についての計算の複雑さの２つの原因が存在する。
【０００７】
１．特定の分子モデルであって、構成原子の位置、速度、および質量特性、構成原子間の原子間力、ならびにそれらの原子と構成原子の周囲環境との間の相互作用を記述するために使用される、分子モデル。
【０００８】
２．経時的にモデルを前進させるために使用される、特定の数値的方法。時間は、最終時刻に達するまで、非常に短い間隔（これは、時間刻みと呼称される）分だけ、繰り返し進められる。
一般的な実施において、分子モデルは、デカルト座標（ｘ，ｙ，ｚ）および溶質分子の各個々の原子の速度からなり、明示的な分子モデルは、個々の溶媒分子（明示的溶媒）または溶媒のバルク性質の解析的近似（非明示的溶媒）のいずれかから構成された溶媒環境のモデルと結合する。数値的方法は、かえる飛び法（ｌｅａｐｆｒｏｇ）積分器、Ｖｅｒｌｅｔ積分器または同様の簡単な積分法からなる（この方法は、Ｖｅｒｌｅｔによる「Ｃｏｍｐｕｔｅｒ’Ｅｘｐｅｒｉｍｅｎｔｓ’ｏｎ  Ｃｌａｓｓｉｃａｌ  Ｆｌｕｉｄｓ：Ｉ．Ｔｈｅｒｍｏｄｙｎａｍｉｃａｌ  Ｐｒｏｐｅｒｔｉｅｓ  ｏｆ  Ｌｅｎｎａｒｄ−Ｊｏｎｅｓ  Ｍｏｌｅｃｕｌｅｓ」，Ｐｈｙｓ．Ｒｅｖ．，１５９（１）：９８−１０３，１９６７年７月によって初めて議論された）。
【０００９】
実質的な研究は、例えば、以下のように分子モデルについての計算負荷を低減することにおいて完結した：剛体仮定条件で高次のモードを拘束することによって、モデルの複雑性の低減すること、剛体性または伸縮性の下部構造分解（ｓｕｂｓｔｒｕｃｔｕｒｉｎｇ）によってモデルを単純化すること、Ｎ次動力学、有効性のある非明示的な溶媒モデル、および力場モデルについての多重極法（ｍｕｌｔｉｐｏｌｅ  ｍｅｔｈｏｄ）（例えば、市販のＭＢＯ（Ｎ）Ｄソフトウエアパッケージに関する米国特許第５，４２４，９６３号を参照のこと））。
【００１０】
これまでは、現在の計算方法が非常に小さな時間刻み（代表的に、１〜１０フェムト秒（１０−１５〜１０−１４秒））を必要とするため、分子シミュレーションが非常に遅かった。採用された各時間刻みは、特定の分子モデルの新しい状態（各原子の位置および運動）の計算を必要とし、次いで、力の新しいセットの計算は、この新しい状態から生じる。例えば、高分子の複雑な挙動（例えば、タンパク質のフォールディング）の分子動力学シミュレーションは、代表的には、少なくとも１マイクロ秒〜数秒（数分でさえも）の時間範囲に及ぶ必要がある。現在、通常用いられている技術によって、これは、コンピュータシミュレーションにおいて１０９〜１０１６の時間刻みの必要性を生じる。この状態、および特に力の、刻み当たりの計算は、問題の規模が大きくなるにつれて非常に高くつく。今日利用可能な最高速のコンピュータを用いたとしても、数ヶ月、数年、または数百年のコンピュータ処理時間が、適度な大きさの系であってもこのような問題を解決するために必要とされる。
【００１１】
化学的過程および物理学的過程の正確なモデルを維持しつつ、時間刻みが顕著に増加する場合に解決され得る、分子モデリング問題の、速度およびサイズにおける多大な改善を達成し得る。これらの小さな時間刻みが、分子結合の振動において見出される非常に高い振動数の存在下で精度を維持する必要性という不可避な必要条件であることが分子動力学者によって広く信じられている。例えば、Ｌｅａｃｈによる、Ｍｏｌｅｃｕｌａｒ  Ｍｏｄｅｌｌｉｎｇ  Ｐｒｉｎｃｉｐｌｅｓ  ａｎｄ  Ａｐｐｌｉｃａｔｉｏｎｓ，１９９６，ｐ．３２８と、Ｂｅｒｅｎｄｓｅｎによる、ｉｎ  Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｍｏｌｅｃｕｌａｒ  Ｄｙｎａｍｉｃｓ：Ｃｈａｌｌｅｎｇｅｓ，Ｍｅｔｈｏｄｓ，Ｉｄｅａｓと、Ｄｅｕｆｌｈａｒｄら（編）による、Ｓｐｒｉｎｇｅｒ，１９９９，ｐｇ．１８と、Ｒａｐａｐｏｒｔによる、Ｔｈｅ  Ａｒｔ  ｏｆ  Ｍｏｌｅｃｕｌａｒ  Ｄｙｎａｍｉｃｓ  Ｓｉｍｕｌａｔｉｏｎ，Ｃａｍｂｒｉｄｇｅ，１９９５（訂正と共に再版された１９９８年、ｐ．５７）と米国特許第５，４２４，９６３号を参照のこと。
【００１２】
しかし、この常識的な確信は誤りである。数値解析のうちのコンピュータサイエンス下位区分によって、高振動数が存在するが、それがほとんど関心事とならない問題に対する、数値積分法の広範な理論が生まれてきた。これらの問題は、「スティッフ（ｓｔｉｆｆ）」問題と呼称される（例えば、ＨａｉｒｅｒおよびＷａｎｎｅｒ、Ｓｏｌｖｉｎｇ  Ｏｒｄｉｎａｒｙ  Ｄｉｆｆｅｒｅｎｔｉａｌ  Ｅｑｕａｔｉｏｎｓ  ＩＩ：Ｓｔｉｆｆ  ａｎｄ  Ｄｉｆｆｅｒｅｎｔｉａｌ−Ａｌｇｅｂｒａｉｃ  Ｐｒｏｂｌｅｍｓ，第２版、Ｓｐｒｉｎｇｅｒ，１９９６を参照のこと）。これらの場合においては、時間刻みを制限するものは、これらの積分法の安定性であって、その必要とされる解の精度ではない。積分器が安定性の性質において広く変動し、この性質は、安定性領域または安定性区間によって厳密に特徴付けられ得る。陽的積分法は、インプリメントされることが容易で、Ｖｅｒｌｅｔ法がその例であり、常に、非常に限定された安定性領域を有している。
【００１３】
一方で、陰的積分方法は、陽的方法よりもずっと複雑であるが、はるかに大きい安定性領域を有し得る。実際には、無条件の安定性を有する陰的積分方法が存在する。これは、理論上、この方法が、任意の大きな時間刻みを採用し得ることを意味する。このような方法は、「Ｌ安定性」と呼ばれる数学的性質を有する。それゆえに、所与のモデルおよび所望の計算について、「十分に安定な」積分方法の選択は、刻み幅が固有の精度の必要条件によってのみ制限されることを可能にする。実際には、陰的方法のみが十分に安定である。Ｌ安定性の方法は、常に十分安定である。さらに、陰的積分方法のみがＬ安定であり得るが、実際には、非常に少数の陰的積分方法がＬ安定である。言い換えると、Ｌ安定性積分方法は、十分安定な陰的積分方法のサブセットであり、その方法自体が全ての陰的積分方法のサブセットである。
【００１４】
この議論において、「大きい時間刻み」は、その幅が固有の精度の必要条件または本質的な収束必要条件によって制限されるが、その積分法の安定性の限界によっては制限されない、時間刻みである。実際には、分子動力学において遭遇する２００フェムト秒（ｆｓ）以上の任意の時間刻みは、この定義によってほとんど「大きい」ことが確実であるが、ほとんどの適用では、多くの非常に小さい時間刻みが大きいと考えられるべきである。共有結合伸縮項を組み込むシステムに対して、Ｖｅｒｌｅｔ安定性の関係によって２ｆｓに制限される。剛体モデルの使用によって除外された結合伸縮を有するシステムに対して、Ｖｅｒｌｅｔ安定性は、典型的には、４０ｆｓ未満に刻み幅を制限する。
【００１５】
数人の分子動力学者が陰的方法を用いて実験し、この方法を非実用的として拒否した。例えば、Ｓｃｈｌｉｃｋ、Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｍｏｌｅｃｕｌａｒ  Ｄｙｎａｍｉｃｓ：Ｃｈａｌｌｅｎｇｅｓ，Ｍｅｔｈｏｄｓ，Ｉｄｅａｓ，Ｄｅｕｆｌｈａｒｄら（編），Ｓｐｒｉｎｇｅｒ，１９９９，ｐ２３８を参照のこと。特に、誘導された減衰を介して、シミュレーションからエネルギーを除去する安定な方法の傾向は、致命的な欠陥であると考えられていた。なぜならば、この方法は、各時間刻みで、非線形系のために必要とされる膨大な計算時間を有しているからである。Ｓｃｈｌｉｃｋによる、ｏｐ．ｃｉｔ．，ｐｐ２３８−９および２４４を参照のこと。この減衰効果は、重要な欠陥であると考えられた。なぜなら、ほとんどの分子動力学シミュレーションは、エネルギーを保存することが必要であるためである。上記で引用されたＳｃｈｌｉｃｋ総説では、分子モデルは、明示的な減衰および安定な積分方法に起因するエネルギー損失を人為的な力によって取り戻すように導入したＬａｎｇｅｖｉｎ項を含んでいた。これらの力によって、事実上、安定な方法が所望されるような大きい時間刻みを取ることが妨げられる。陰的方法はこのような計算に実際に使用され得るが、エネルギー保存を必要としない、分子モデリング計算がまた、多く存在し、そして、本発明者らの方法は特に、そのような問題に対して有効である。本発明者らは、賢明なモデリング選択および慎重なインプリメンテーションを介して、実際の計算において有効な方法を教示する。
【００１６】
従来技術の陰的方法による成功がない結果、現在の分子モデリングシミュレーションツールは、主に、Ｖｅｒｌｅｔによって１９６７年に最初に議論されたエネルギー保存、シンプレクティック（ｓｙｍｐｌｅｃｔｉｃ）陽的積分方法に依存する。これらの積分方法の改変（例えば、ｌｅａｐｆｒｏｇ法、または速度Ｖｅｒｌｅｔ法および改変されたＢｅｅｍａｎ法など）は、Ｔｉｎｋｅｒ（Ｊａｙ  Ｐｏｎｄｅｒ，ＴＩＮＫＥＲ  Ｕｓｅｒ’ｓ  Ｇｕｉｄｅ，Ｖｅｒｓｉｏｎ３．８，Ｏｃｔｏｂｅｒ  ２０００，Ｗａｓｈｉｎｇｔｏｎ  Ｕｎｉｖｅｒｓｉｔｙ，Ｓｔ．Ｌｏｕｉｓ，ＭＯ）などの現在の分子動力学コードにおいて利用可能である。
【００１７】
低振動数成分および高振動数成分を分離することによってか、または、特別なＶｅｒｌｅｔ誘導法（例えば、ＳＨＡＫＥおよびＲＡＴＴＬＥ）と組み合わせて高振動数の結合振動を拘束することによって、時間刻みを大きくする他の最近の試みは、時間刻み幅を増加することにおいて限定された成功が有していた。２倍〜５倍だけの速度増加が達成された（Ｅｒｉｃ  Ｂａｒｔｈらによる、「Ａ  ｓｅｐａｒａｔｉｎｇ  ｆｒａｍｅｗｏｒｋ  ｆｏｒ  ｉｎｃｒｅａｓｉｎｇ  ｔｈｅ  ｔｉｍｅｓｔｅｐ  ｉｎ  ｍｏｌｅｃｕｌａｒ  ｄｙｎａｍｉｃｓ」，Ｃｏｍｐｕｔｅｒ  Ｓｉｍｕｌａｔｉｏｎ  ｏｆ  Ｂｉｏｍｏｌｅｃｕｌａｒ  Ｓｙｓｔｅｍｓ，Ｖｏｌ３．，ｐｐ９７−１２１，１９９７）。
【００１８】
要約すると、分子モデリング（特に分子動力学シミュレーション）において、努力の成果が微小な刻み幅によって妨げられてきた。積分が、依然として非常に小さい時間刻みで実施され、その結果、極端に手間のかかる計算を生じ、結果を得るのに時間がかかる。分子の研究における有用用途に対する障害であることが明らかである。結果を得るのに１年を費やす分子動力学シミュレーションなどは、実際の研究のために使用され得ない。対照的に、本発明は、有用かつ正確な計算結果を迅速に生じるように大きな時間刻みでの積分を可能にする方法を教示する。
【００１９】
これらの問題を回避するために、本発明は、特定の挙動または関心のある性質を計算する場合に計算時間を低減するための方法を教示する。
【００２０】
（発明の要旨）
本発明は、環境における分子のシステムの挙動または性質を計算する方法を教示する。本発明は、環境効果および換算座標（ｒｅｄｕｃｅｄ  ｃｏｏｒｄｉｎａｔｅ）で表現された分子について運動方程式を用いて分子系を数学的にモデリングする工程；および所望の挙動および性質の正確な計算を得るように大きい時間刻みで十分安定な積分器を用いてモデル方程式を積分する工程を含む。本方法は、計算時間の最適な使用のために、精度かつ収束の必要条件に従って、時間刻みのサイズを変動させる工程を含む。この時間刻みのサイズ（幅）は、少なくとも１００の範囲内で変動し得る。
【００２１】
好適な換算座標分子モデルは、デカルト全原子座標ではなく、ねじれ角座標を組み込んだ、剛体分割法である。好適な十分に安定な積分方法は、誤差制御された動力学的な計算のためのＬ安定一段階Ｒａｄａｕ５法、およびエネルギー極小化（静的）計算のためのＬ安定性陰的Ｅｕｌｅｒ法を含む。厳密性が少ない安定性必要条件を有する用途に対して、高安定性および効率的な多段階ＤＡＳＳＬ法が好ましい。
【００２２】
（特定の実施形態の説明）
本発明に従うモデリングモジュールのための、ソフトウェアの一般的なシステムアーキテクチャ４８およびそのいくつかのプロセスを、図１に示す。大きい矩形ブロックの各々は、ソフトウエアモジュールを表し、矢印はソフトウエアモジュール間を通過する情報を表す。このソフトウエアシステムアーキテクチャは、モデラー（ｍｏｄｅｌｅｒ）モジュール５０、生化学的コンポーネントモジュール５２、物理学的モデルモジュール５４、解析モジュール５６、および可視化モジュール５８を有する。これらのモジュールのいくつかの詳細は、以下に説明される。他のモジュールが、公で利用可能である。
【００２３】
モデラーモジュール５０は、ユーザが特定の分子系を規定する物理学的パラメータを入力するインターフェースを提供する。インターフェースはグラフィックまたはデータファイル入力（またはその両方）を有し得る。生物学的コンポーネントモジュール５２は、分子系の特定の数学モデルのためのモデラー入力を変換し、モデリングされる系の、分子、力場、および溶媒のそれぞれが、数学的モデリングされる、変換サブモジュール６０、６２、および６４に、分割される。例えば、Ｔｉｎｋｅｒ（Ｊａｙ  Ｐｏｎｄｅｒ，「ＴＩＮＫＥＲ  Ｕｓｅｒ’ｓ  Ｇｕｉｄｅ」，Ｖｅｒｓｉｏｎ３．８，２０００年１０月，Ｗａｓｈｉｎｇｔｏｎ  Ｕｎｉｖｅｒｓｉｔｙ，Ｓｔ．Ｌｏｕｉｓ，ＭＯ）を含む利用可能ないくつかのモデラーおよび生物学的コンポーネントモジュールが存在する。
【００２４】
生物学的コンポーネントモジュール５２から変換された物理パラメータを用いて、物理学的モデルモジュール５４は、分子系を数学的に定義する。多体系サブモジュール６６が、モジュール５４の中心に存在する。物理学的モデルモジュール５４および多体系サブモジュール６６が、以下に詳細に説明される。同時係属出願「ＭＥＴＨＯＤ  ＦＯＲ  ＡＮＡＬＹＴＩＣＡＬ  ＪＡＣＯＢＩＡＮ  ＣＯＭＰＵＴＡＴＩＯＮ  ＩＮ  ＭＯＬＥＣＵＬＡＲ  ＭＯＤＥＬＩＮＧ」と題された米国特許出願第＿＿号、および「ＭＥＴＨＯＤ  ＦＯＲ  ＲＥＳＩＤＵＡＬ  ＦＯＲＭ  ＩＮ  ＭＯＬＥＣＵＬＡＲ  ＭＯＤＥＬＩＮＧ」と題された米国特許出願第＿＿号の両出願は同日に出願され、以前に引用した仮特許出願からの優先権を主張し、両出願は、本譲受人に譲渡され、本明細書中で参考として援用され、本発明の観点からこれらの特許出願に開示された物理学的モデルモジュール５４および多体系モデルモジュール６６のさらなる説明を有する。
【００２５】
物理モデルモジュール５４および可視化モジュール５８と通信する解析モジュール５６は、物理学的モデルモジュール５４によって規定された分子系の計算モデルに解を提供する。解析モジュール５６は、物理学的モデルモジュール５４の異なる方程式を積分する、積分器サブモジュール６８のセットからなる。積分器サブモジュール６８は、経時的に分子系を前進させ、さらに分子系の極小エネルギー配置を決定する際に使用される静的解析を提供する。本発明の内容のほとんどを含むのは、解析モジュール５６および積分器サブモジュールジュール６８であり、以下に詳細に説明される。
【００２６】
可視化モジュール５８は、生物学的コンポーネントモジュール５２および解析モジュール５６からの入力情報を受信し、ユーザに分子系の三次元画像表示および分子系に対して得られた解を提供する。多くの可視化モジュールが現在利用可能であり、例えば、ＶＭＤ（Ａ．Ｄａｌｋｅら、「ＶＭＤ  Ｕｓｅｒ’ｓ  Ｇｕｉｄｅ」，Ｖｅｒｓｉｏｎ１．５，２０００年６月、Ｔｈｅｏｒｅｔｉｃａｌ  Ｂｉｏｐｈｙｓｉｃｓ  Ｇｒｏｕｐ，Ｕｎｉｖｅｒｓｉｔｙ  ｏｆ  Ｉｌｌｉｎｏｉｓ，Ｕｒｂａｎａ，Ｉｌｌｉｎｏｉｓ）である。
【００２７】
（分子モデルおよび多体系の説明）
以下で説明される積分器は、多体系（ｍｕｌｔｉｂｏｄｙ  ｓｙｓｔｅｍ；ＭＢＳ）に関する分子モデルの運動を記述する方程式のセットについて演算する。以下に詳細に説明された積分法の計算を支援するために、本発明に従って、ねじれ角、剛体モデルが、対象とする分子系を記述するために使用される。内部座標（選択された、一般化座標および一般化速度）が分子の状態を記述するために使用される。
【００２８】
このＭＢＳは、モデリングされる分子系を構成する、原子および有効な剛体性の結合の抽象化であり、かつ、シミュレーションによって扱われる問題についての重要な特徴を失うことなく、実際の物理学系、その環境下にある分子を単純化するように選択される。図１に示された一般的なシステムアーキテクチャに関して、このＭＢＳは、原子間の静電荷または他のエネルギー的相互作用を含まないだけでなく、この分子を囲む溶媒のモデルも含まない。力場は、生化学的コンポーネントのサブモジュール６２でモデリングされ、そして、溶媒は、生化学的コンポーネントのサブモジュール６４でモデリングされる。
【００２９】
図２は、対象分子のＭＢＳのツリー構造（ｔｒｅｅ  ｓｔｒｕｃｔｕｒｅ）を示す。このＭＢＳの基本的な抽象化は、ヒンジ接続された剛体１７０の１つ以上の集合の抽象化である。剛体は、その物体を構成する全ての粒子が互いに対して固定された位置を有する物理的な物体の数学的抽象化である。伸縮または他の相対運動は許容されない。ヒンジ接続は、２つの剛体間の許容される相対運動を規定する数学的抽象化である。これらの剛体およびヒンジ接続の例は、以下に説明される。
【００３０】
基本物体（ｂａｓｅ  ｂｏｄｙ）１７２と呼ばれる１以上の物体は、運動力学的状態が、基準（ｇｒｏｕｎｄ）１７４上の基準点に対して直接参照される点において固有の状態を有する。系のグラフは、１つ以上の「ツリー」である。ツリーの重要な特性は、任意の物体から任意の他の物体への経路が固有であり、すなわち、このグラフはループを含まないということである。ツリー内の物体は、ｎ個存在する（この出発点（ｂａｓｅ）は、標識１を有する）。ツリー内の物体は、規則的な標識が割り当てられ、このことによって、これらの標識が、基本物体から任意のリーフ物体（ｌｅａｆ  ｂｏｄｙ）１７６までの任意の経路上で決して減少しないことが意味される。リーフ物体は、単一の他の物体のみに接続される物体である。規則的な標識は、標識ｎをリーフ物体１７８（少なくとも１つでなければならない）に割り当てることによって達成され得る。この物体がグラフから除去される場合、その場合のツリーは、ｎ−１物体を有する。次いで、標識ｎ−１は、そのリーフ物体１８０のうちの１つに割り当てられ、全ての物体が標識化されるまで、このプロセスが繰り返される。これはまた、系内の任意の残っているツリーに対して行われる。
【００３１】
物体間の関係を維持することを支援するために、整数型の関数が使用され、系の各物体について内側物体を記録する。各出発点（ｂａｓｅ）対しては、内側物体（ｉｎｂｏａｒｄ  ｂｏｄｙ）が基準（ｇｒｏｕｎｄ）であり、物体ｋ（１８４）についてのペアレント物体すなわち内側物体１８２（ｉ）は、ｉ＝ｉｎｂ（ｋ）として参照される。さらに、記号Ｎは、慣性系、すなわちグラウンド系１７４をいう。上付き文字Ｏは、座標原点（０，０，０）をいう。
【００３２】
ある点から別の点までのベクトルについての記号は２点の名前を含む。従って、ｒＰＱは、点Ｐから点Ｑまでのベクトルである。基準系内の点の速度を表すベクトルは、以下のように点および基準系の名称を含む：ＮｖＰ。以後導入される特定の記号は、２つの基準系を関連付ける。この場合では、記号は２つの系の名前を含む。従って、ｉＣｋは、系ｉにおける系ｋの配向のための方向余弦行列である。この記号はそのペアレント系における代表的物体のための方向余弦行列を指す。従って、ｉＣｋ（ｊ）は、問題にしている現実の物体ｊを示す。左および右の上付き文字は、物体のインデックスを変更しない。これはまた、他の記号についても真である。アスタリスクは、転置を示す（例えばＨ＊（ｋ））。ベクトルの上のチルダ記号は、３×３歪対称クロス乗積、すなわち
【００３３】
【数１】
2004519026.tif 000002
を表す。
【００３４】
【数２】
2004519026.tif 000003
は、ｉ×ｉの恒等行列であり、そして
【００３５】
【数３】
2004519026.tif 000004
は、長さｉの零ベクトルであり、
【００３６】
【数４】
2004519026.tif 000005
は、ｉ×ｉの零行列である。
【００３７】
（モデルの剛体）
図３は、ＭＢＳのサンプル「ツリー」の基準配置１９０を示す。１つを超えるツリーが許容される。各物体の点は、そのヒンジポイントＱとして表される。例えば、点Ｑｋ  １８６は、物体ｋ  １８４のヒンジポイントである。座標軸の固定された集合は、慣性系１９８に確立される。ＭＢＳの任意の配置は、その基準配置１９０として選択される。この配置における慣性座標軸の写像は、各物体内で物体に固定された軸の集合を確立するために使用される。基準配置では、各ヒンジポイントＱはＰ（そのペアレント物体（または拡張された物体）の点）と一致する。各物体に対して、点Ｐは、物体の内側ヒンジポイントと呼ばれる。そのため、物体ｋ  １８４についての内側ヒンジポイント  Ｐｋ  １８８は、ペアレント物体ｉ  １８２に固定された点である。各基本物体の内側ヒンジポイントは、グラウンドに固定された点Ｏ  １９２である。図２に示された拡大図は、点Ｑｋ  １８６が、本体ｋ  １８４に固定され、点Ｐｋ１８８はペアレント物体ｉ  １８２に固定されることをより明確に示す。
【００３８】
ヒンジポイントの位置は、各物体についての定ベクトルである
【００３９】
【数４Ａ】
2004519026.tif 000006
１９４を規定し、これはまた
【００４０】
【数５】
2004519026.tif 000007
と記載され得る。物体ｋについてのベクトルは、そのペアレント物体ｉに固定される。このベクトルは、物体ｉのヒンジポイントから物体ｋの内側ヒンジポイントの範囲にある。ベクトル
【００４１】
【数５Ａ】
2004519026.tif 000008
１９６は、慣性系原点から第１の基本体の内側ヒンジポイント（またグラウンド内に固定された点）の範囲にあり、
【００４２】
【数６】
2004519026.tif 000009
と記載され得る。
【００４３】
物体に対して、ｍ（ｋ）、
【００４４】
【数７】
2004519026.tif 000010
は、そのヒンジ点Ｑｋに対する物体ｋの質量の特性を規定する。これらはそれぞれ、質量、第１質量モーメント、その物体の座標系にある、そのヒンジポイントについての物体の慣性行列である。粒子の分布を構成する剛体について、その質量特性は、予備処理モジュールによって計算される定数である。それらの計算の詳細は、ＫａｎｅによるＴ．Ｒ．Ｄｙｎａｍｉｃｓ，第３版、１９７８年１月、Ｓｔａｎｆｏｒｄ  Ｕｎｉｖｅｒｓｉｔｙ，Ｓｔａｎｆｏｒｄ，ＣＡ等の標準的な参考文献において見出され得る。
【００４５】
ヒンジポイントＱｋに対する物体ｋの空間的な慣性Ｍ（ｋ）は、対称６×６行列
【００４６】
【数８】
2004519026.tif 000011
によって与えられる。
【００４７】
この系における各ジョイントは幾何学的データによって説明される。例えば、ピンジョイントは、そのジョイントによって接続された２つの物体に固定された軸によって特徴付けられる。ジョイントついての具体的なデータは、その型に依存する。番号ｎ、ｉｎｂ関数、系の質量特性、ベクトル
【００４８】
【数８Ａ】
2004519026.tif 000012
、およびジョイント幾何学データ（ジョイント型を含む）は系パラメータを構成する。
【００４９】
（モデルのジョイントおよび一般化座標）
図４は、ＭＢＳの好適な実施形態のジョイントの定義（スライダージョイント１００、ピンジョイント１０２、およびボールジョイント１０４）を示す。各ジョイントは、内側ヒンジポイントＰｋ  １０８に関してヒンジポイントＱｋ  １０６の並進変位または回転変位を可能にする。これらの変位は、物体ｋについての一般化座標ｑ（ｋ）１１０によりパラメータ表示される。ちなみに、一般化座標は一般化量の一例であり、回転特性および並進特性の両方を有する量を指すことに留意されるべきである。例えば、点において作用する一般化力は、力ベクトルおよびトルクベクトルの両方からなる。スライダージョイント１００についての一般化座標ｑ（ｋ）は滑り変位ｘ１１２である。ピンジョイント１０２についての一般化座標ｑ（ｋ）は、角変位θ１１４である。ボールジョイント１０４についての一般化座標ｑ（ｋ）は、Ｅｕｌｅｒパラメータ（ε１，ε２，ε３，ε４）１１６である。
【００５０】
ジョイントの各々は、ピン、スライダー、またはボールジョイント；あるいはこれらのジョイントの組み合わせであり得る。多くの他のジョイント型は、フリージョイント、Ｕジョイント、シリンダージョイント、およびベアリングジョイントを含むがこれらに限定されない、これらのジョイントのタイプの組み合わせによって可能となり得る。例えば、基本物体内側ヒンジポイントから基本物体ヒンジポイントまでのベクトルの慣性測定数、ｑ（ｋ）＝（ｘ，ｙ，ｚ）は、３つの直交スライダージョイントとして、グラウンドにおける基本物体の変位を表現する。フリージョイントは、ボールジョイントと結合された３つの直交スライダージョイントからなり、全部で６つの自由度を有する。
【００５１】
全ての物体に対する一般化座標の集合は、システムの一般化座標であるベクトル
【００５２】
【数８Ｂ】
2004519026.tif 000013
を含む。
【００５３】
特定のジョイントについての一般化座標が与えられると、以下の２つの量：
ジョイント並進ベクトル
【００５４】
【数９】
2004519026.tif 000014
および
ペアレントにおける本体ｋに対する方向余弦行列
ｉＣｋ（ｋ）
が形成される。
この並進ベクトル
【００５５】
【数１０】
2004519026.tif 000015
は、ペアレント物体の座標系における物体ｋの内側ヒンジポイントＰから物体ｋのヒンジポイントＱへのベクトルを表す。これらの計算の詳細は、ジョイント型に依存し、容易に導出される。この記述の目的のために、系の一般化座標で与えられた
【００５６】
【数１１】
2004519026.tif 000016
を生成する関数の呼び出しが、前提となる。
【００５７】
導入されたように、各物体についてのヒンジポイントの選択は任意である。しかし、賢明な選択によって、問題が大幅に簡略化される。例えば、ピンジョイントに対して、このヒンジポイントは、ジョイントの軸上にある点として選択されるべきである。この選択に対して、点Ｐおよび点Ｑは、ジョイント角度の全ての値に対して一致したままであるため、このジョイントの並進は０である。点Ｑがこの軸から離れて選択される場合、点ＰおよびＱは、以下：
【００５８】
【数１２】
2004519026.tif 000017
のように互いに対して移動し、ここで、
【００５９】
【数１２Ａ】
2004519026.tif 000018
はジョイント軸の単位ベクトルであり、θはそのジョイント角度であり、
【００６０】
【数１３】
2004519026.tif 000019
は、軸上の任意の点から点Ｑへのベクトルである。
【００６１】
ピンジョイントおよびボールジョイントについて、本発明者は、軸上の点をヒンジポイントとして常に選択する。これらのジョイントに対して、並進ベクトル
【００６２】
【数１４】
2004519026.tif 000020
は０である。
【００６３】
スライダージョイントについて、並進ベクトル
【００６４】
【数１５】
2004519026.tif 000021
は
【００６５】
【数１５Ａ】
2004519026.tif 000022
である。
【００６６】
ピンに対する方向余弦行列は、
【００６７】
【数１６】
2004519026.tif 000023
であり、スライダーに対する方向余弦行列は、
【００６８】
【数１７】
2004519026.tif 000024
である。
【００６９】
（モデルの一般化座標）
ペアレントｉにおいて測定された物体ｋのヒンジポイントの一般化速度を、一般化速度の集合ｕ（ｋ）によってパラメータ表示されたｉＶｋ（ｋ）とすると
【００７０】
【数１８】
2004519026.tif 000025
であり、
ここで、行列Ｈ（ｋ）は、このジョイントに対してジョイントマップと呼ばれる。これは６×ｎｕ（ｋ）の行列であり、ここで、ｎｕ（ｋ）は、そのジョイントに対する自由度の数（ピンまたはスライダーに対して１、ボールに対して３、フリージョイントに対して６）である。Ｈ（ｋ）は、一般的には座標ｑに対する依存性を有する。ジョイントについての一般化速度が与えられると、ジョイントマップは、チャイルド（子）物体の系において表現されたジョイントの線形速度および角速度を生成する。これらのジョイントとして、本発明者らは、
【００７１】
【数１９】
2004519026.tif 000026
を使用する。
この全ての物体についての一般化速度の集合は、ベクトル
【００７２】
【数２０】
2004519026.tif 000027
、すなわち、この系についての一般化座標を含む。先述のように、（ｑ，ｕ）および特定のジョイント型が与えられた場合、ベクトルｉＶｋ（ｋ）を生成し得る関数の呼び出しが前提となる。以下の微分：
【００７３】
【数２１】
2004519026.tif 000028
を計算し得る関数の呼び出しがまた、前提となる。このルーチンは、以下の一般化位置座標の時間微分：
【００７４】
【数２２】
2004519026.tif 000029
を生成し、ここで、Ｗ（ｑ）は、
【００７５】
【数２３】
2004519026.tif 000030
およびｕに関連するブロック対角行列であり、以下のジョイント型：
【００７６】
【数２４】
2004519026.tif 000031
およびフリーのジョイントに依存するブロックの各々は、３つのシリンダージョイントおよび１つのボールジョントの組み合わせである。ボールジョイントについて、３つのｕと関連する４つの
【００７７】
【数２５】
2004519026.tif 000032
（Ｅｕｌｅｒパラメータの微分）が存在することに留意のこと。
【００７８】
同様に、ｉＡｋ（ｋ）、すなわち、そのペアレント（ｐａｒｅｎｔ）の物体ｋのヒンジポイントの一般化加速度は、以下：
【００７９】
【数２６】
2004519026.tif 000033
によって与えられる。
【００８０】
ここで計算されるのは、この記述目的のための、分子系の、これらの一般化座標
【００８１】
【数２７】
2004519026.tif 000034
、一般化速度
【００８２】
【数２８】
2004519026.tif 000035
、内部座標である。代表的な慣性座標（ｘ，ｙ，ｚ）およびこれらの慣性座標系における速度を用いて取り組むよりも、対象分子系についての計算が軽減される。
【００８３】
（第１運動学的計算）
分子系の内部座標、すなわち、
【００８４】
【数２９】
2004519026.tif 000036
およびこの系のパラメータが与えられると、各物体ｋについての、以下の位置、速度および加速度の運動学的状態が計算される。
【００８５】
物体ｋの各々について、
【００８６】
【数３０】
2004519026.tif 000037
を計算する。これらの計算は、反復的に実行され、基本物体の各々から開始され、そのリーフへと進行する。
【００８７】
ＮＣｋ（ｋ）、すなわち、基準（ｇｒｏｕｎｄ）におけるその物体ｋについての方向余弦行例が、以下のように定義される：
ＮＣｋ（１）＝ｉＣｋ（１）
ＮＣｋ（ｋ）＝ＮＣｋ（ｋ）ｉＣｋ（ｋ），ｋ＝２，．．．ｎ，ｉ＝ｉｎｂ（ｋ）
。
ｉＣｋ（ｋ）は、上記のジョイントルーチンから得られる。
【００８８】
そのペアレント系において表現される場合、物体ｋのペアレントのヒンジポイントＱｉからの、物体ｋのヒンジポイントＱｋへの位置ベクトルである、
【００８９】
【数３１】
2004519026.tif 000038
は、以下：
【００９０】
【数３２】
2004519026.tif 000039
のように定義される。
【００９１】
【数３３】
2004519026.tif 000040
は、ジョイントルーチンから得られる。
【００９２】
グローバル系として表現される場合、慣性系の原点Ｏから物体ｋのヒンジポイントＱｋへの位置ベクトルである
【００９３】
【数３４】
2004519026.tif 000041
は、以下：
【００９４】
【数３５】
2004519026.tif 000042
のように定義される。
【００９５】
物体ｋについての剛体の変換演算子ｉφｋは、以下：
【００９６】
【数３６】
2004519026.tif 000043
のように定義される。
【００９７】
物体ｋの系で表現される場合、そのヒンジポイントでの物体ｋについて空間ベクトルＶ（ｋ）は、以下：
【００９８】
【数３７】
2004519026.tif 000044
のように定義される。
【００９９】
物体ｋの系で表現される場合の、そのヒンジポイントにおける、物体ｋについての空間的加速度Ａ（ｋ）は、以下：
【０１００】
【数３８】
2004519026.tif 000045
で定義され、ここで、
【０１０１】
【数３９】
2004519026.tif 000046
である。
当然、これらの計算は、所望される場合、全て、単一の経路で計算され得る。
【０１０２】
時間刻み１増分についてのこれらの工程を完結した後に、ＭＢＳは、任意の物体の任意のポイントについて、（一般化）位置、（一般化）速度または（一般化）加速度の情報を計算するための運動論的要求に助力を提供する。これは、標準的な剛体式を使用して、その物体についてのヒンジ量に関する任意のポイントについて必要とされる情報を計算することによってなされる。
【０１０３】
（動力学的残差工程）
分子モデルの所定状態、すなわち、所定の
【０１０４】
【数４０】
2004519026.tif 000047
および系のパラメータを用いて開始して、プログラムルーチンは、ＭＢＳの「環境」をモデリングする。このようなルーチンは、コンピュータモデリング分野の当業者にとって容易に利用可能であり、作製可能であり得る。このルーチンは、積分サブモジュール６８の形態によって決定され、そして通過する、値（ｑ，ｕ）をとり、（状態依存的な）そのヒンジポイントＱｋにおける物体ｋについて適応された空間的力
【０１０５】
【数４１】
2004519026.tif 000048
、物体ｋについてのヒンジトルクσ（ｋ）を返す。図１に示される生化学的コンポーネントモジュール５２における力場モジュール６２および溶媒モジュール６４に基づいて、Ｔ（ｋ）およびσ（ｋ）は物理学的モデルモジュール５４において計算される。その後、物体ｋについての一般化速度ｕ（ｋ）に関するこの動力学的残差ρｕ（ｋ）は、以下の工程１および２：
１．各物体についての空間的負荷収支
【０１０６】
【数４２】
2004519026.tif 000049
の生成：
【０１０７】
【数４３】
2004519026.tif 000050
２．ρｕ（ｋ）の計算
【０１０８】
【数４４】
2004519026.tif 000051
によって、計算される。
【０１０９】
この動力学的残差ρｕ（ｋ）は、モデルについての運動方程式の（直接形式（ｄｉｒｅｃｔ  ｆｏｒｍ）とは対照的な）残差形式のために、現れる。差分方程式の残渣形式および直接形式ならびにそれらの積分の詳細な説明は、「ＭＥＴＨＯＤ  ＦＯＲ  ＲＥＳＩＤＵＡＬ  ＦＯＲＭ  ＩＮ  ＭＯＬＥＣＵＬＡＲ  ＭＯＤＥＬＩＮＧ」と題する、同一日付けの上記の同時係属中の米国特許出願番号＿＿＿＿に見出される。
【０１１０】
（第２運動学的計算）
Ｐ（ｋ）、Ｄ（ｋ）、ｉΨｋ（ｋ）、ｉＫｋ（ｋ）を計算する：
１．各物体の有節（ａｒｔｉｃｕｌａｔｅｄ）物体の慣性Ｐ（ｋ）の初期化。
【０１１１】
Ｐ（ｋ）＝Ｍ（ｋ）、ｋ＝１，．．．，ｎ
２．対象の生成
【０１１２】
【数４５】
2004519026.tif 000052
これらの量の関数依存性は、ｑについてのみである。
【０１１３】
（前進的動力学計算）
以下を計算する：
【０１１４】
【数４６】
2004519026.tif 000053
。
【０１１５】
（直接形式法）
直接形式法は、現在の状態（ｑ，ｕ）を取り、そしてその微分
【０１１６】
【数４７】
2004519026.tif 000054
を上記のアルゴリズムを使用して計算し、その後、このアルゴリズムは、時間を前進する積分法によって使用される。状態（ｑ，ｕ）を用いて開始した場合、以下の
【０１１７】
【数４８】
2004519026.tif 000055
を計算する：
１．上記のジョイント特異的なルーチンを使用して
【０１１８】
【数４９】
2004519026.tif 000056
を計算する
２．上記の第１運動学的計算を、
【０１１９】
【数５０】
2004519026.tif 000057
を用いて実施する
３．動力学残差計算を使用して残差ρｕを生成し、以下：
ρｕ＝−ρｕ
を否定する
４．第２運動学的計算を実施する
５．前進的動力学計算を実施し、
【０１２０】
【数５１】
2004519026.tif 000058
を計算する。
直接形式法は、この系について作用する適用された力に応答してヒンジ加速度
【０１２１】
【数５２】
2004519026.tif 000059
を生成する。ここで、
【０１２２】
【数５３】
2004519026.tif 000060
は、数値的方法を経過して、この分子モデルの運動方程式を積分する。
【０１２３】
（分子モデルの運動方程式を積分するための数値的方法）
上記で説明したように、分子系をモデルリングする努力は、これまでは、尋常でない量のコンピュータパワーおよび時間を要求していた。慎重に選択された分子モデルを用い、そして内部座標を使用したとしても、上述のように、運動方程式が積分されなければならない。これまで、これらの努力は、この分子系を規定するのに使用される差分方程式の、微小な時間刻みにおける積分に集中してきた。しかし、この差分方程式を大きな時間刻みにおいて積分することの単純な必要条件によって、分子モデリングの複雑な問題を解けない。より筋の通ったアプローチが必要とされる。
【０１２４】
（スティッフＭＤシミュレーションを解く）
初期値問題として提示される、常微分方程式（ＯＤＥ）または微分代数方程式（ＤＡＥ）の系を数値積分することを試みる場合、この最も大きな時間刻みが、所望される解の精度によってか、または使用される積分モデルの安定性によって制限され得る。陽的積分法（ｅｘｐｌｉｃｉｔ  ｉｎｔｅｇｒａｔｉｏｎ）を使用する場合の時間刻みは、所望される解の精度のみによって制限され、検討中の系が、「非スティッフ」とであると考えられる。しかし、この積分法は、「ブローアップ」する傾向が有るか、または検討中の系についての予測され得た時間刻みよりずっと微小である時間刻みにおいて不安定になる場合、用語「スティッフ（硬い；Ｓｔｉｆｆ）」は、この状態（すなわち、最も大きな時間刻みが具体的な積分法の安定性によって制限される）を記述するために使用される。
【０１２５】
本発明は、減衰していない高振動数（および、それによる、極めて小さい時間スケールにおける正確な解）が、重要なものではなく、かつ、この分子系のモデリングの長い時間スケールでの解に影響を及ぼさない系のモデリングに関する。所謂「スティッフ」系の問題の例は、１秒間で前後に振動する単振子のモデリングであり得る。ここで、極めて小さい質量が、非常に硬いバネを用いてこの振子の端に結合している。この質量およびバネの系の固有振動数は、例えば、１０００サイクル／秒である。すなわちこの振子の振動の各々について、この質量は、１０００回振動する。さらに、この小さな質量の高振動数の振動は、その小さな振幅のためにほとんど感知されず、そして本発明者らが研究している挙動に対する有意な方法のいずれにおいても大きなスケールの振動運動に影響を及ぼさない。時間刻みおよび誤差制御を伴った陽的積分法は、この振動する振子のモデルを解くために適用される。この高振動数の振動が誤差許容よりもずっと小さいときでさえもこの積分器が極めて微小な時間刻みをとる場合、この系は、「スティッフ」である。
【０１２６】
実施される単純な実験は、既知量、例えば、１０倍に誤差許容を緩めて同じ研究を再び実行することである。採用した時間刻みサイズが、積分器のオーダーが与えられて期待されたおおよその量だけ、増加しない場合、この問題はスティッフとなる。大きな時間刻みをとる試みによって、「ブローアップ」する積分法が生まれた。この性質は、この積分法の単なる人為的な結果である。本発明は、多くのこれまでの分子モデリングシミュレーションで本質的であった時間刻みのサイズに対するこのスティッフ限界を回避する。分子モデリングの問題におけるこの部類に取り組むために、本発明は、図１の積分器サブモジュール６８のための「十分に安定性である」陰的積分法を使用する。本発明者らは、「十分に安定性である」ことの厳密な定義を以下のように提供するが、上記の誤差許容の調整実験は、実際問題として機能し、この時間刻みサイズが、誤差許容の設定に対して期待されるように応答する場合、この方法は、当面の問題について十分に安定である。あるいは、本発明者らは、Ｌ安定性の方法を選択し得る。なぜならば、これらは、常に十分に安定性であるからである。
【０１２７】
陰的方法に対する導入として、単純なＥｕｌｅｒ積分法を考察する。以下：
【０１２８】
【数５４】
2004519026.tif 000061
のＯＤＥを積分するために、このＥｕｌｅｒ法の陽的バージョンは、その前の（過去の）解の周りで打ち切られたＴａｙｌｏｒ級数展開式、つまり、
ｙｎ＝ｙｎ−１＋ｈｎｆ（ｙｎ−１）
を使用する。すなわち、次の時間刻み幅ｈｎに対するｙｎについての解は、その前の解であるｙｎ−１にのみ依存する。従って、ｙｎは、この方程式の左辺のみにあり、直接的かつ陽的に解かれ得る。対照的に、Ｅｕｌｅｒ積分法の陰的バージョンは、その次の解の周りで打ち切られたＴａｙｌｏｒ級数展開式、つまり、
ｙｎ＝ｙｎ−１＋ｈｎｆ（ｙｎ）
を使用して、この展開式によって、所望の解ｙｎがその方程式の両辺にある方程式を生じ（故に、ｙｎについて陰的である）、従って、以下の方程式：
ｇ（ｙｎ）＝ｙｎ−ｙｎ−１−ｈｎｆ（ｙｎ）＝０
を解くために非線形の反復法（通常は、Ｎｅｗｔｏｎ法のいずれかのバージョン）が要求される。この積分法における一見したところ単純な変化は、非線形性の反復工程の実施を必要とするという多大な損失をして、この方法の安定性における劇的な変化をもたらす。
【０１２９】
安定性関数Ｒ（ｚ）の評価をすることによって、積分法の安定性を決定することが可能であり、この関数は、任意の積分法について書き示めされる。これらの安定性関数の導出は、直接的であるが、極めて複雑でもある。詳細は、ＨａｉｒｅｒおよびＷａｎｎｅｒ、Ｓｏｌｖｉｎｇ  Ｏｒｄｉｎａｒｙ  Ｄｉｆｆｅｒｅｎｔｉａｌ  Ｅｑｕａｔｉｏｎ  ＩＩ：Ｓｔｉｆｆ  ａｎｄ  Ｄｉｆｆｅｒｅｎｔｉａｎｌ−Ａｌｇｅｂｒａｉｃ  Ｐｒｏｂｌｅｍｓ、第２版、Ｓｐｒｉｎｇｅｒ、１９９６において見出される。本発明に従って、Ｌ安定性として既知の安定性の強力な形式によって、任意の分子モデリング問題についての十分な安定性が保証される。Ｌ安定性積分法は、Ａ安定性積分法として知られている、より低い安定性積分法から、強力な部分集合を形成する。多くの場合、Ａ安定性の方法またはＡ（α）安定性のようなより弱い方法はまた、十分に安定性である。
【０１３０】
数学的に、安定性関数Ｒ（ｚ）を伴う安定性領域（ｄｏｍａｉｎ）は、以下：
【０１３１】
【数５５】
2004519026.tif 000062
であり、ここで、
【０１３２】
【数５６】
2004519026.tif 000063
は、複素平面を表し、ここで、ｚは、ｚ＝ｘ＋ｉｙの形式の複素数である。具体的な問題の安定性は、ｚ＝ｈλと割り当てることによって大まかに試験され得、ここで、ｈは時間刻みであり、かつ、
【０１３３】
【数５７】
2004519026.tif 000064
は、積分される系の線形モデルの固有値であり、ここで、ωは、減衰固有振動数であり、ζは、減衰定数である。通常は、この方法の安定性を制限する固有値λは、この系における最高の振動数固有値である。一般に、この振動数が高くなればなるほど、この安定性の限界に達する前に使用される得る時間刻みｈは、小さくなる。大きなコンフォメーション変化を受ける特定の非線形性モデルについての十分な安定性を正確に決定するために、そのコンフォメーションの各々の周りで線形化されたこの系の全ての固有値が、その安定性領域に存在することを決定しなければならない。
【０１３４】
安定性関数の安定性変域Ｓから、この陰的積分法が、Ａ安定性であること、つまり、以下：
【０１３５】
【数５８】
2004519026.tif 000065
であるか否か（すなわち、少なくとも複素平面
【０１３６】
【数５９】
2004519026.tif 000066
の左半分全体が包含される場合、この方法はＡ安定性であるということ）を決定することが可能である。複素平面
【０１３７】
【数６０】
2004519026.tif 000067
における安定性領域Ｓの範囲を使用して、この積分法がＡ安定性か否かを決定する。
【０１３８】
この方法がＡ安定性である場合、この方法は、以下のようなＬ安定性の基準を満たし得る：
もし、
【０１３９】
【数６１】
2004519026.tif 000068
である場合、この方法はＬ安定性であり、かつ任意の問題について十分に安定性である。
【０１４０】
図５Ａ〜図５Ｃは、種々の既知の積分法についての安定性を例示する。これらの図面において、特定の積分法は、その安定性関数Ｒ（ｚ）と共に左側に与えられ、複素平面
【０１４１】
【数６２】
2004519026.tif 000069
におけるその安定性領域Ｓが、Ａ安定性である（または、Ａ安定性でない）ことの決定とともに中央に示され、そしてＬ安定性であることの決定が右側に示される。
【０１４２】
この陰的Ｅｕｌｅｒ積分法（これの安定性は図５Ａに例示される）が、最も強力なＬ安定性積分法であると認識され、これは、その大きな安定性変域およびシミュレーションにおける高振動数の迅速な減衰に起因する。この陰的な中点（ｍｉｄ−ｐｏｉｎｔ）法は、図５Ｂが示すように、明らかにＡ安定性であるが、Ｌ安定性ではない。このＲａｄａｕ５積分法は、図５Ｃが示すようにＬ安定性であり、その解における誤差の特に良好な制御を有するというさらなる特性を有している。スティッフ性（硬さ；ｓｔｉｆｆｎｅｓｓ）のさらなる説明、陰的積分解の技術ならびにＡ安定性およびＬ安定性は、上記で引用したＨａｉｒｅｒおよびＵ．Ａｓｃｈｅｒ、Ｃｏｍｐｕｔｅｒ  Ｍｅｔｈｏｄｓ  ｆｏｒ  Ｏｒｄｉｎａｒｙ  Ｄｆｆｅｒｅｎｔｉｎａｌ  Ｅｑｕａｔｉｏｎｓ  ａｎｄ  Ｄｉｆｆｅｒｅｎｔｉａｌ−Ａｌｇｅｂｒａｉｃ  Ｅｑｕａｔｉｏｎｓ、ＳＩＡＭ、Ｐｈｉｌａｄｅｌｐｈｉａ、ＰＡ、１９９８において見出され得る。
【０１４３】
興味深いことに、分子動力学シミュレーションにおいて使用される一般的な積分器である、Ｖｅｒｌｅｔ法は、陽的方法であるが、Ａ安定性もＬ安定性も持ち合わせてはいない。この方法についての安定性「区間（ｉｎｔｅｒｖａｌ）」は、大まかに、以下のように与えられる（Ｌｏｐｅｚ−Ｍａｒｃｏｓ，Ａｎ  ｅｘｐｌｉｃｉｔ  ｓｙｍｐｌｅｃｔｉｃ  ｉｎｔｅｇｒａｔｏｒ  ｗｉｔｈ  ｍａｘｉｍａｌ  ｓｔａｂｉｌｉｔｙ  ｉｎｔｅｒｖａｌ，Ｒｅｐｏｒｔ  ｏｆ  ｔｈｅ  Ｄｅｐａｒｔｍｅｎｔ  ｏｆ  Ａｐｐｌｉｅｄ  Ｍａｔｈｅｍａｔｉｃｓ，Ｕｎｉｖｅｒｓｉｄａｄ  ｄｅ  Ｖａｌｌａｄｏｌｉｄ，Ｓｐａｉｎ，１９９５）：
ｈ＜Ｌ
ここで、
【０１４４】
【数６３】
2004519026.tif 000070
の形式で計算されるＭＤ方程式についてＬ＝２／ωであり、かつ、ωは線形化モデルの最も高い振動数固有値である。ほとんどのＭＤシミュレーションについて、分子結合振動の高振動数によって、ｈが、約１〜２フェムト秒未満に制限される。ＳＨＡＫＥまたはＲＡＴＴＬＥを使用して、この最も高い振動数の結合振動を固定することによって、僅かならがら状況を改善し、約１０フェムト秒までの時間刻みを可能にする。しかし、この安定性問題は依然として存在する。
【０１４５】
本発明は、進歩が遅かった分子モデリングの少なくとも２つの分野で有意な前進を与える。第１の分野は、「静的解析（ｓｔａｔｉｃ  ａｎａｌｙｓｉｓ）」の分野であり、この分野は、所定の配置から開始して局所的なエネルギー極小を決定する問題を扱う分野である。これは、大域的な極小値を検索する間に遭遇する副次的問題を解くために使用され得る。すなわち、複雑な分子の化学成分が与えられた場合、例えば、この分子の安定な極小エネルギー配置はいかなるものであろうか。このような解が極めて有用な分子系の例は、タンパク質の、最終または中間の折り畳み配置である。本発明が直接的に有用な第２の分野は、ＭＤとしばしば呼称される、分子動力学の分野であり、この分野においては、分子系の時間履歴が所望される。分子系の初期座標が与えられた場合、分子動力学によって、経時間的にこの系の変化が調べられる。例えば、タンパク質の結合ポケットとの薬物リガンドの動力学的相互作用が決定され得る。
【０１４６】
（静的解析）
静的解析は、検討中の分子系の極小エネルギー配置を決定するために使用される。重要な極小エネルギー配置は、局所的な極小値または大域的な極小値であり得、そして、しばしば、この系についての機能的な配置（例えば、酵素または他の折り畳まれたタンパク質についての作動可能な配置）を表す。
【０１４７】
静的解析について好ましい実施形態は、換算座標（ｒｅｄｕｃｅｄ  ｃｏｏｒｄｉｎａｔｅ）分子モデルに、Ｌ安定性積分器を適用することであり、これは、この系から最大エネルギーを受け、そして、安定な配置に到達することが可能な最も大きい時間刻みを取ることが可能である。剛体系およびねじれ角換算モデルに適用される、陰的Ｅｕｌｅｒ（ＩＥ）積分法は、本発明に従う静的解析にとって好ましい実施形態である。単純に１次法である場合、この陰的Ｅｕｌｅｒ法は、大きな誤差が生じ、この誤差によって各時間刻みにおいて大きなエネルギー吸収を生じる。この安定性領域は、最も知られたものの１つであり、従って、大変大きな時間刻みを可能にする。この時間刻みは、一般的に、この非線形系の解が収束できる能力のみによって限定される。これが探索された極小エネルギー配置であり、経時的なこの分子系の特定の挙動ではないので、この方法によって生じた大きな誤差は、この結果の精度を妨げるほどではない。第２の可能な実施形態は、その誤差制御が無力であるＲａｄａｕ５である。
【０１４８】
陰的Ｅｕｌｅｒ積分法は、ベクトル関数
【０１４９】
【数６４】
2004519026.tif 000071
（ここで、ｙ＝（ｑ，ｕ）であり、ｑはこの分子系の位置状態を表し、そしてｕはこの分子系の速度状態を表している）についての図６のフローチャートで例示される。この関数ｆは、例えば、多体系動力学ならびに力（例えば、静電的な引力および斥力、ファンデルワールス力、および溶媒和力）の両方を含む。エントリー工程７９の後に、第１の操作工程８０は、反復行列Ｇを更新する。全ての陰的積分法について、この反復行列Ｇは、Ｇ＝Ｉ−αＪの形式を有し、ここで、Ｉは、単位行列であり、αは時間刻みｈｎのスカラー関数であり、ここで、この時間刻みは、ｔｎとｔｎ−１との間であり、そして、
Ｊは、以下：
【０１５０】
【数６５】
2004519026.tif 000072
によって、与えられるヤコビアンである。陰的Ｅｕｌｅｒ法について、α＝ｈｎである。ちなみに、コンピュータ使用時間をさらに節約するために、残差形式（ｒｅｓｉｄｕａｌ  ｆｏｒｍ）の方程式からヤコビの行列を計算する極めて効率的方法は、「ＭＥＴＨＯＤ  ＦＯＲ  ＡＮＡＬＹＴＩＣＡＬ  ＪＡＣＯＢＩＡＮ  ＣＯＭＰＵＴＡＴＩＯＮ  ＩＮ  ＭＯＬＥＣＵＬＡＲ  ＭＯＤＥＬＩＮＧ」と題する先述の同時係属中の米国特許出願第＿＿＿＿＿号において扱われ、そして現在の指定代理人に帰属されることに留意すべきである。本発明の場合と同様に、同一の参照された特許出願はまた、この分子系の状態を記述するために内部座標の使用を記載している。例えば、この分子の一部分の回転は、外部の参照座標に対してではなく、（この分子の）別の部分に対して記述される。このことによって、計算効率がさらに向上する。
【０１５１】
Ｎｅｗｔｏｎ反復法（ニュートン法の説明については、Ａｓｃｈｅｒ，ｏｐ．ｃｉｔ．を参照のこと）に従う工程の８２の系列によって、時刻ｔｎにおけるこの分子系の位置状態および速度状態を反復的に見出す。この状態ｙは、全ての位置状態および速度状態を代表する。ｙｎを見出すための反復は、ｙの変化が許容差Ｔｏｌ１の内にあるか、または可能な最大反復数であるｉｍａｘに到達するかのいずれかの場合に、終了する。許容差Ｔｏｌ１および最大反復数ｉｍａｘは、全体の効率を最大化するために経験的に調整される。代表的な値は、Ｔｏｌ１＝１０−４およびｉｍａｘ＝１０である。
【０１５２】
記号
【０１５３】
【数６６】
2004519026.tif 000073
は、ベクトルの２−ノルムを取ることを表現する。
【０１５４】
【数６７】
2004519026.tif 000074
を解くために反復行列Ｇの逆行列を計算するよりも、数値解析において周知の技術である、ＬＵ因数分解のようなより安定な線形解法技術を使用することが慣例的であることに留意すべきである。工程８４は、収束について試験する。この収束が満たされた場合に、状態ｙおよび時刻ｔは更新され、この時間刻みｈｎは、工程８８によって示されるように増加する。別な方法において、この時間刻みｔｎは、工程８６によって低減され、そして改変されたＮｅｗｔｏｎ反復法である系列８２が再び試みられる。この静的解析は、試験工程８７において、時間刻みが微小でありすぎる場合、作動しなくなる。工程８８において時間区間を２倍にすることまたは工程８６において時間区間を半分にすることは、この時間積分区間が変化する様子の単純な例であることに留意のこと。しばしば、より洗練されたアルゴリズムが、公に利用可能である積分法として使用される。
【０１５５】
状態ｙおよび時間区間が更新された後に、決定工程９０は、工程の最大許容数が実施されたか否かを試験する。工程の最大数ｎｍａｘがとられた場合、静的解析は作動しなくなる。他の場合、状態ｙｎにおけるその速度ｕｎは、工程９１においてゼロ（０）に設定され、そしてこの加速度
【０１５６】
【数６８】
2004519026.tif 000075
が、これらが許容差Ｔｏｌ２よりも小さいか否かを確認するために工程９２において試験される。もしそうであるならば、静的解析は成功である。他の場合、この工程は、工程９４においてインクリメントされ、そしてこのプロセスは、工程８０に返り、その反復行列などを更新する。代表的な値は、Ｔｏｌ２＝１０−５およびｎｍａｘ＝５００である。
【０１５７】
（分子動力学法）
分子モデリングの別の目的は、分子動力学、つまり、タンパク質のフォールディングまたはタンパク質における活性部位とリガンドとのドッキングのような、分子系における物理学的過程の時間履歴を正確に決定するシミュレーションである。
【０１５８】
本発明に従って、問題としている分子系をモデリングするＯＤＥは、誤差制御を伴う十分に安定な積分法によって時間について積分される。誤差制御を伴う高次（少なくとも２次）の十分に安定な積分器によって、必要とされる精度が提供されるが、これは、モデルにおける関連性のない高振動数を迅速に減衰させる。最も大きい可能な時間刻みが取られ、所望の制度が達成され；積分は、安定性問題によって制限されない。積分におけるける安定性問題に起因する時間刻みのサイズに対する制限なしに、精度と計算時間との間の諸条件の考慮がなされ得る。
【０１５９】
好ましい実施形態は、陰的Ｒａｄａｕ５積分法（特に、Ｔｙｐｅ  Ｒａｄａｕ  ＩＩＡ，５次元の陰的Ｒｕｎｇｅ−Ｋｕｔｔａ積分器）である。上記で参照のことされたＨａｉｒｅｒ、１１８〜１２７頁を参照のこと。Ｒａｄａｕ５はＬ安定性であり、故に、全てのモデルおよび環境について十分に安定性である。この積分法のインプリメントのフローチャート概観は、図７に示される。Ｒａｄａｕ５法は、３ステージを有する単一工程の陰的積分器である。従って、これは、図３に示される陰的Ｅｕｌｅｒ法と類似の構造を有するが、１ステージであはなく、３ステージを有し、複素代数学および行列変換を含む種々の方法を取り込み、操作総数および丸め誤差を低減する。このＲａｄａｕ５法はまた、ユーザ特異的な精度の要件に従って時間刻み幅を制御するための誤差推定器を有する。
【０１６０】
開始工程１１０の後に、ヤコビの行列Ｊが工程１１２で更新される。陰的Ｅｕｌｅｒ法と同様に、改変されたＮｅｗｔｏｎ反復が工程１１４で実施され、ここで、反復行例Ｇは、以下：
【０１６１】
【数６９】
2004519026.tif 000076
であり、そして残差関数
【０１６２】
【数７０】
2004519026.tif 000077
は、行列Ａおよび行列関数Ｆを含み、ここで、Ｒａｄａｕ５法の３ステージが展開される。記号
【０１６３】
【数７１】
2004519026.tif 000078
は、テンソル積を意味する。示された用語の詳細な説明、ならびに以下で説明される誤差推定器の用語については、Ｈａｉｒｅｒ，ｏｐ．ｃｉｔ．を参照のこと。
【０１６４】
この反復行列の収束は、工程１１６において試験される。この反復が最大反復数ｉｍａｘにおける許容差Ｔｏｌ１を満たさない場合、刻み幅（ｓｔｅｐｓｉｚｅ）ｈｎは、工程１１８において低減されて、そして、この反復は再度試され、試験工程１２０においてこの最小の刻み幅ｈｍｉｎに到達しなければ、解析的は作動しなくなる。代表的な値は、Ｈａｉｒｅｒに提供される。
【０１６５】
一旦、この反復が許容されると、この状態は工程１２２において更新され、そして、新しい刻み幅ｈｎが、誤差推定ｅｒｒに基づいて計算され、このｅｒｒは、Ｈａｉｒｅｒにおいて説明されるような絶対許容差および相対許容差の関数である。最終時刻ｔｆｉｎａｌに試験１２４で到達した場合、この動力学解析は首尾よく成功したことになる。他の条件がまた、ｔｆｉｎａｌに到達することの代りに、またはｔｆｉｎａｌに到達することに加えて、停止について試験され得る。そうでなければ、この刻み（ｓｔｅｐ）ｎが、工程１２６によってインクリメントされ、そして、このループを継続する。実際問題として、ｔｆｉｎａｌに到達すること以外の条件（例えば、運動エネルギーまたはポテンシャルエネルギーの指示されたレベル）が使用され、完了が示され得る。
【０１６６】
（本発明の適用例）
本発明の利点を例示するために、陰的Ｅｕｌｅｒ積分法、Ｒａｄａｕ５積分法、および従来技術であるＶｅｒｌｅｔ積分法が、本発明者らによって分子シミュレーション問題に適用された。図８は、２残基のアラニンジペプチド１５０を用いたタンパク質フラグメントの構造を例示し、このアラニンジペプチドについて安定、すなわち「静的」な、極小エネルギー配置が存在することが公知である。アラニンジペプチドは、アミノ酸の式Ａｌａ−Ａｌａおよび化学式ＮＨ３＋−ＣＨ−ＣαＨ−ＣＨ３−ＣＯＮＨ−ＣαＨ−ＣＨ３−ＣＯＯ−を有しており、ここで、Ｃαは、各残基中のα炭素であり、ここで、各残基の間のＣＯＮＨは硬いペプチド結合１５４である。この多体系の記述は、１つの物体において７つの原子１５２を有する７つの物体を含む。各物体は、共に硬く結合していると考えられる１以上の原子から構成される。この７つの物体は、総数２３の原子を表現している。この剛体の間の結合は、この物体が互いに対して回転することを可能にするピンジョイントで表現される共有結合である。このペプチド結合１５４のいずれかのサイドにおけるピンジョイントのうちの２つは、立体配置角（ｃｏｎｆｉｇｕｒａｔｉｏｎ  ａｎｇｌｅ）、φ１５６およびΨ１５８によって表現される。アラニンジペプチドのこのモデルは、可能な極小エネルギー配置（φが約−１４７°であり、Ψが約１６２°である）を有する。
【０１６７】
図９Ａ〜９Ｆにおけるグラフは、この３つの積分法の結果を例示する。図９Ａ〜図９Ｃは、それぞれ、Ｖｅｒｌｅｔ積分法、Ｒａｄａｕ５積分法、および陰的Ｅｕｌｅｒ積分法における立体配置角Ψについて結果を示し、そしてこれらの全ては、比較目的で同一の軸を有する。この縦軸は、角度である。同様に、図９Ｄ〜図９Ｆは、これらの３方法における立体配置角φについての結果を示し、そして、これらの全てはまた、比較目的のため同一軸を有する。これらの縦軸は、角度である。この横軸は、ＣＰＵ時間（８００ＭＨｚ  Ｐｅｎｔｉｕｍ（登録商標）ＩＩＩマイクロプロセッサを有するパーソナルコンピュータおける、秒）の対数目盛であり、各々のシミュレーションを完結させるために必要な時間を比較した。３つの全てのシミュレーションは、立体配置角、つまりΨ＝１３５°かつφ＝−１３５°について同一の初期座標から開始し、大まかに同じ結果を伴って終了した。
【０１６８】
標準的なＶｅｒｌｅｔ積分法は、この問題を解くのに約２９００秒を必要とするが、この陰的Ｅｕｌｅｒは、約２．５秒のみを必要とし、これは、同一のコンピュータ上において１０００倍以上速い。この陰的Ｅｕｌｅｒ解法はよりなめらかであり、Ｖｅｒｌｅｔ積分法が示したようなアラニンジペプチド分子系の不必要な高振動数成分を追跡することはないことに留意のこと。期待されたように、最終的な正しい解はこの高振動数成分に依存しない。
【０１６９】
Ｒａｄａｕ５積分法は、Ｖｅｒｌｅｔ法より７０倍速く、４０秒間を必要とした。この陰的Ｒａｄａｕ５解法は、「雑音性（ｎｏｉｓｙ）」であって、そして重要な挙動を追跡したが、Ｖｅｒｌｅｔ法が示したタンパク質のフラグメントの不必要な高振動数成分の追跡をおこなうことはなかった。期待されたように、最終的な解は、この不必要な高振動数成分に依存しなかった。
【０１７０】
図１０Ａ〜図１０Ｃによって、図９Ａ〜図９Ｆにおいて議論した３つのシミュレーションの各々について刻み幅（フェムト秒）に対するＣＰＵ時間（秒）が図示される。図１０Ａ〜１０Ｃのグラフにおいて、両方の軸は対数目盛であることに留意のこと。図１０Ａは、陽的Ｖｅｒｌｅｔ積分器によって達成され得る一定の１０フェムト秒の刻み幅を示す。図１０Ｂは、Ｒａｄａｕ５の刻み幅がシミュレーションの開始時の約１００フェムト秒から１０８フェムト秒（すなわち１００ナノ秒）へと増加することを示す。図１０Ｃは、この陰的Ｅｕｌｅｒの刻み幅がシミュレーション開始時の約１フェムト秒から１０４フェムト秒へと増加することを示す。これらの長期の刻み幅は、従来技術のＭＤシミュレーションにおいて前例のないものである。
【０１７１】
十分に安定な積分法（例えば、Ｌ安定性の方法）は、換算座標分子モデルの任意の形態に対して適用され得、本発明に従って分子モデリングにおける問題を解決するために使用され得る。このようなモデルとしては、以下が挙げられるが、これらに限定されない：
１）緊密なループおよび他の代数的な拘束ならびに開放性の木構造を用いた、分子の拘束系モデル；
２）上記ねじれ角動力学モデル以外の分子モデルの他の換算形式（例えば、下位構造化モデル（ｓｕｂｓｔｒｕｃｔｕｒｅｄ  ｍｏｄｅｌ））；
３）常微分方程式または微分代数方程式の残差形式ならびにその直接形式；
４）完全Ｎｅｗｔｏｎ法および他の反復技術、ならびに非線形方程式を解くために使用されるこの反復技術のための改変Ｎｅｗｔｏｎ法の使用；
５）数値的に誘導したヤコビアンならびに解析的に誘導したヤコビアンの使用；
６）部分的な全原子モデル、剛体モデル、伸縮モデル（ｆｒｅｘｉｂｌｅ  ｍｏｄｅｌ）、それらの組み合わせ、または分子の原子構造の他の任意の表現の使用；
７）全原子モデル（例えば、水または他の明示的な溶媒、薬物ならびに他の低分子）との換算座標モデルとの組み合わせの使用；
８）好ましい実施形態において示したような方法を含むが、これらに限定されない時間刻みサイズを調整するための種々の方法の使用；
９）Ｒａｄａｕ５および陰的ＥｕｌｅｒのＬ安定性積分器に加えて、積分器のＳＤＩＲＫ、ＳＩＲＫおよびＲｏｓｅｎｂｒｏｃｋのファミリーを含むがこれらに限定されない、誤差制御を伴うかまたは伴わない他のＬ安定性陰的積分法；
１０）ＤＡＳＳＬおよびＯＤＥまたは微分代数方程式（ＤＡＥ）のための他の多段法を含むが、これらに限定されない、他の十分に安定な方法。
【０１７２】
本発明に従う適切に換算分子モデルを使用する十分に安定な積分器を用いて、正確な分子モデリングがコンピュータ上で実行され得る速度は、劇的に改善され、そして本発明の有利な点が明白である。特に、本発明は、タンパク質のフォールディングに適用される場合に非常に有用であり、なぜならば、これらは、完結するのに極めて長い時間（通常は、事実上、ミリ秒から秒のオーダーである）を要する、長期間のスケールの反応であるためである。分子動力学に対する現行のアプローチは、大変遅いために、数ナノ秒より長い、低分子タンパク質以外の全てのタンパク質フォールディング操作をシミュレートすることはできない。本発明は、タンパク質の構造を決定するためのタンパク質のフォールディング問題を解くために重要なツールを提供する。現存の計算技術および実験技術用いて決定し得ない構造のタンパク質（例えば、膜結合タンパク質）が、本発明を用いて取り組まれ得る。何百万もの既知タンパク質の構造を経験的決定するための膨大な時間およびコストが回避される。本発明は、タンパク質のネイティブ構造が迅速に決定され、かつ薬物および他のタンパク質とのこれらの相互作用がシミュレートされるために、合理的な薬物およびタンパク質の設計を強化する。タンパク質のフォールディング経路、構造、および機能に対する研究が有意に増強される。
【０１７３】
本発明は、多くの他の生体分子（例えば、ＲＮＡ、ＤＮＡ、多糖類、および脂質）をシミュレートするために使用され得る。また、タンパク質−ＲＮＡ複合体（例えば、リボゾーム）およびタンパク質−ＤＮＡ複合体（例えば、ヒストンおよびクロマチン中のＤＮＡ）のようなこれらの生体分子の組み合わせの分子構造がシミュレートされ得る。タンパク質の構造を改変する過程（例えば、シャペロンタンパク質によるタンパク質の翻訳後修飾）は、シミュレートされ得る。
【０１７４】
（さらなる適用）
本発明は、計算分子モデリングに関係する多くのアルゴリズムにおける中核となる計算として使用され得る。例えば、アルゴリズムは、いくつかの所望の基準（例えば、統計的分布）に従う初期条件のセットを選択し得、そして多くの別個の分子動力学の実行にて各々の初期配置としてこのセットの１つのメンバーをとる。各実行が大規模な並列計算の一部分として別個のコンピュータ上にてなされ得、そのいくつかまたは全てが単一のコンピュータにおいて行われ得る。本発明は、分子動力学法を実行するために使用され得；次いで、さらなる処理のための高次のアルゴリズムによって決定が得られる。別のアルゴリズムは、リボソームの配置またはタンパク質の押し出しのシミュレーション（ここで、この分子モデルは、アミノ酸として増加し、そして物理的に現実的な速度または他に選択された速度でこのタンパク質に加えられる）であり、本発明が発生しているタンパク質の各長さの挙動および特性をシミュレートするのに使用される。別のクラスのアルゴリズムは、本発明に使用して実行される、時々発生するエネルギー増加事象とエネルギー保存シミュレーションまたはエネルギー散逸シミュレーションとを混合するシミュレーションである。このようなアルゴリズムは、代表的には溶媒によって生成される温度−浴効果（例えば、Ｌａｎｇｅｖｉｎ項）または、モデル温度効果を機能的または統計的にモデリングするように設計される他のエネルギー増加効果を獲得するように設計された入力を含む。
【０１７５】
本発明はまた、分子系の設計または改善を実施することを試みるアルゴリズムにおいて中核となる計算して有用である。これらのアルゴリズムにおいて、本発明を使用して特定の系の特性を計算し得る。これらの特性は、「設計パラメータ」と呼称される、特定された変化または変化の型（これらは、設計過程または改善過程の一部分として系に対して生成される）によって、変化され得る。本発明を使用して分析される場合、設計パラメータに対する変化の結果として生じる特性に対する変化について得られる情報を使用して、この設計パラメータに対するさらなる変化を指向し、所望の特性についての改善を導く。例えば、特定のリガンドと緊密に結合するタンパク質が所望される。はじめに、タンパク質−リガンド系が、本発明によって解析され、結合親和性特性がその結果として計算される。このタンパク質の個々のアミノ酸は、設計パラメータであると考えら得る。１つ以上のアミノ酸の対する変化が、いくつかのアルゴリズムに従って実施され、そのアルゴリズムは、無作為であるか、またはより洗練されたものであり得る。次いで、この結合親和性が、本発明を使用して再度計算され得る。結合親和性に対する得られた変化を使用して、特定のリガンドに対する所望の結合親和性に対する改善を生み出す配列が発見されるまで、アミノ酸に対するさらなる改変が誘導される。新規タンパク質が、実験室において、合成され、そしてリガンドに対して試験され、結果の妥当性を立証し、かつこの新規のタンパク質が医薬的用途または商業的な用途を有し得る可能性を決定し得る。
【０１７６】
他の設計アルゴリズムは、この分子モデルの任意のパラメータおける改善（特に、経験的に導出された力場、および溶媒特性を含む）を含む。これらのアルゴリズムは、異なった種類の換算座標モデル（例えば、アミノ酸を、目的の特性（例えば、電荷または疎水性）によって特徴付けられる単純なエレメントに抽象化したモデル）において実行される。
【０１７７】
分子構造が既知である場合、本発明の方法は、代替物、または従来の生物化学的スクリーニング法に対する補助薬としての標的との相互作用について、化合物ライブラリーをスクリーニングするために特に有用である。次いで、本発明のモデリング法により同定される、所望の様式において標的と相互作用と考えられる化合物または化合物のサブセットは、合成され、従来の生化学的アッセイによって試験される。従って、本発明は、合成される必要のある化合物数および他の場合において、所望の機能特性を有する化合物を同定するために必要とされる生化学的アッセイの数を低減する。本発明は、この用途について他のコンピュータ技術によりも優れている。なぜならば、これは、スクリーニングの間における標的およびリガンドの両方のコンフォメーション変化（フレキシビリティ）を可能にし、従って、予測の精度を著しく向上するからである。
【０１７８】
上述した一般的なアプローチに従って、この方法は、標的と化合物の相互作用についてのモデル（この化合物およびこの標的についての運動方程式を含む）を提供する。陰的積分の効果的な使用のために、このモデルは換算座標を使用すべきである。
【０１７９】
スクリーニングされるこの化合物およびこの標的に関するデータは、入力のために、この運動方程式に提供される。これらのデータは、ユーザによって提供され得るか、または記憶ファイル、リモートデータベースから、もしくは測定機器から獲得され得る。いくつかの例において、この化合物および／または標的は、化学的名称によって記述される。他の例において、この化合物または標的は、構成分子（例えば、アミノ酸またはヌクレオチドの配列）によって記述される。他の例において、これの化合物または標的は、構成原子およびこれらの原子を共に保持する結合の性質によって記述される。さらに、また、あるいは、化合物および／またはその標的が、実験データ（Ｘ線パターン、赤外線スペクトル、紫外線スペクトルまたは核磁気共鳴スペクトル）または、これらのデータに基づいて計算された情報（例えば、原子間距離、回転自由度、および励起状態）によって記述され得る。いくつかの方法において、化合物がその標的と相互作用する溶媒または他の環境（例えば、リン脂質基質）の同一性および／または組成のような、さらなるデータが提供され得る。いくつかの方法において、化合物および標的が相互作用する他の環境要因（温度または圧力）が提供される。
【０１８０】
運動方程式が解かれ、その標的との化合物の相互作用モデルを生成する。このモデルは、スクリーン上に表示される。相互作用に関する種々のパラメータはまた、（その標的との化合物の結合親和性、その標的とその化合物の会合についての速度定数、およびその化合物の特定原子とその標的の特定原子間の距離）出力であり得る。いくつかの例において、スクリーニングされる化合物とその標的との相互作用は、その標的と所望の様式で相互作用することがすでに知られている化合物の相互作用と比較される。その標的との好ましい相互作用は、結合親和性の強度、結合運動の速度、化合物と標的との適合の緊密度、シグナル伝達を示す標的におけるコンフォメーション変換の誘導、標的中の特定の原子に対する化合物中の特定原子の近接性によってか、または、その標的と所望の様式で相互作用することがすでに知られているコントロール化合物に対する化合物の適合の類似性によって、アッセイされ得る。いくつかの方法において、界面活性について化合物をスクリーニングするように、好ましい相互作用が、スクリーニングされる化合物によってそれが変性されることを示す標的の特異的な構造の消失によって示される。いくつかの方法において、モデルまたはモデルに基づくデータは、各化合物がスクリーニングされた後に、表示される。他の方法において、複数または全ての化合物は、スクリーニングされ、そしてあるサブセットのみに対してモデルまたはデータが表示される。
【０１８１】
本方法が使用されて、従来の方法においてスクリーニングされた化合物に同一または類似の型の化合物をスリーニングするのに使用され得る。このような化合物としては、ペプチド、抗体を含むタンパク質、低分子（５００ｋＤａ以下）、βターン模倣物、多糖類、リン脂質、ホルモン、プロスタグランジン、ステロイド、芳香性化合物、複素環式化合物、ベンゾジアゼピン、オリゴマー性Ｎ置換グリシンおよびオリゴカルバメートが挙げられる。大規模の化合物のコンビナトリアルライブラリーは、Ａｆｆｙｍａｘ、ＷＯ９５／１２６０８、Ａｆｆｙｍａｘ  ＷＯ９３／０６１２１、Ｃｏｌｕｍｂｉａ  Ｕｎｉｖｅｒｓｉｔｙ、ＷＯ９４／０８０５１、Ｐｈａｒｍａｃｏｐｅｉａ、ＷＯ９５／３５５０３およびＳｃｒｉｐｐｓ、ＷＯ９５／３０６４２（これらの各々は、全ての目的について、参考として援用される）に記載されるコード化合成ライブラリー（ｅｎｃｏｄｅｄ  ｓｙｎｔｈｅｔｉｃ  ｌｏｂｒａｒｙ；ＥＳＬ）法である。ペプチドライブラリーはまた、ファージディスプレイ法によって生成され得る。例えば、Ｄｅｖｌｉｎ、ＷＯ９１／１８９８０を参照のこと。構造データが、海洋、微生物、藻類、植物および真菌のような供給源から利用可能な天然化合物がまた、スクリーニングされ得る。いくつかの例において、スクリーニングされる化合物は、生化学的アッセイによって達成されるか、または他の場合、標的と所望の相互作用を有する、１以上の化合物を含む。このような化合物が、類似の相互作用を有する他の化合物を同定するためのコントロールとしての役割を担う。例えば、ファージディスプレイ技術を使用し、多くの抗体または他のポリペプチドを得て、そして、それらを標的との相互作用についてスクリーニングすることは、比較的容易である。しかし、抗体またはポリペプチドは、それらの大きなサイズおよび腸管で分解される傾向に起因して、時折、それら自体、治療用途（特に、経口投与）について適切ではない。本方法は、抗体または他のポリペプチドに対する類似の相互作用治療用途についてさらに改善された特徴（例えば、経口バイオアベイラビリティー）を有する、低分子の等価物を同定することを可能にする。
【０１８２】
いくつかの方法において、スクリーニングされる化合物の同一性は、任意のモデリングが実行される前に、予め決定され得る。他の方法において、１つの化合物と標的との間の相互作用が決定され得、次いで、スクリーニングされる次の化合物が、第２の化合物が標的との相互作用を改善することが期待されるような様式で、設計される。いくつかの方法において、スクリーニングされる化合物は、カーネル化合物（ｋｅｒｎｅｌ  ｃｏｍｐｏｕｎｄ）またはリード化合物の改変体を表す。他の方法において、化合物は、本質的に無作為にスクリーニングされる（例えば、ランダムペプチドのコレクションである）。スクリーニングされ得る化合物数は、従来の方法によりも有意に大きい。合成を必要とする従来のスクリーニングまたは化合物の個別化されたスクリーニングにおいて、数千もの化合物をスクリーニングすることでさえも極めて困難であり得る。対照的に、標的と化合物の相互作用をモデリングする本発明は、ずっと短い時間しか要さない本発明において、より多くのオーダー（例えば、１０４、１０６、１０８、１０１０、または１０１５）の化合物がスクリーニングされ得る。
【０１８３】
それに対して化合物がスクリーニングされる標的は、とりわけ、タンパク質、核酸、多糖類、脂質、または有機化学構造体であり得る。しばしば、その標的は、生物学的高分子であり、そしてその標的とその化合物の相互作用によって、この標的を作動するか、またはその標的に拮抗することを介して薬理学的効果を誘導することが所望される。この方法は、そのネイティブ環境から単離されたときにネイティブコンフォメーションを失う標的（例えば、膜結合タンパク質）の相互作用についてスクリーニングすることにおいて特に有用である。目的の標的としては、抗体（抗イディオタイプ抗体および自己免疫疾患（例えば、糖尿病、多発性硬化症、および慢性関節リウマチ）にて存在する自己抗体を含む）が挙げられる。他の目的の標識は、増殖因子レセプター（例えば、ＦＧＦＲ、ＰＤＧＦＲ、ＥＦＧ、ＮＧＦＲ、およびＶＥＧＦ）ならびに、それらのリガンドである。他の標的は、Ｇタンパク質レセプターであり、そして、サブスタンスＫレセプター、アンギオテンシンレセプター、αアドレナリン作動性レセプターおよびβアドレナリン作動性レセプター、セロトニン作動性レセプターならびにＰＡＦレセプターを含む。Ｇｉｌｍａｎ，Ａｎｎ．Ｒｅｖ．Ｂｉｏｃｈｅｍ．５６：６２５〜６４９（１９８７）を参照のこと。他の標識としては、イオンチャネル（例えば、カルシウムチャネル、ナトリウムチャネル、カリウムチャネル）、ムスカリンレセプター、アセチルコリンレセプター、ＧＡＢＡレセプター、グルタミン酸レセプター、およびドーパミンレセプターが挙げられる（Ｈａｒｐｏｌｄの５，４０１，６２９およびＵＳ５，４３６，１２８を参照のこと）。他の標識は、接着タンパク質（例えば、インテグリン、セレクチン）および免疫グロブリンスーパーファミリーのメンバーを含む（Ｓｐｒｉｎｇｅｒ，Ｎａｔｕｒｅ  ３４６：４２５〜４３３（１９９０），Ｏｓｂｏｒｎ、Ｃｅｌｌ  ６２：３（１９９０）；Ｈｙｎｅｓ，Ｃｅｌｌ  ６９：１１（１９９２）を参照のこと）。他の標識は、サイトカイン（例えば、インターロイキンのＩＬ−１〜ＩＬ−１３、腫瘍壊死因子αおよび腫瘍壊死因子β、インターロイキンα、インターロイキンβ、およびインターロイキンγ、腫瘍増殖因子β（ＴＧＦ−β）、コロニー刺激因子（ＣＳＦ）および顆粒球単球コロニー刺激因子（ＧＭ−ＣＳＦ））である。Ｈｕｍａｎ  Ｃｙｔｏｋｉｎｅｓ：Ｈａｎｄｂｏｏｋ  ｆｏｒ  Ｂａｓｉｃ  ＆  Ｃｌｉｎｉｃａｌ  Ｒｅｓｅａｒｃｈ（Ａｇｇｒａｗａｌら編、Ｂｌａｃｋｗｅｌｌ  Ｓｃｉｅｎｔｉｆｉｃ，Ｂｏｓｔｏｎ，ＭＡ  １９９１）を参照のこと。他の標的は、ホルモン、酵素、ならびに分子内メッセンジャーおよび分子間メッセンジャー（例えば、アデニルシクラーゼ、グアニルシクラーゼおよびホスホリパーゼＣ）である。薬物はまた、目的の標的である。標的分子は、ヒト由来、哺乳動物由来、細菌由来であり得る。他の標的は、微生物病原体（ウイルス性および細菌性の両方）ならびに腫瘍由来の、タンパク質、糖タンパク質、および炭水化物）である。なお他の標的は、ＵＳ４，３６６，２４１に記載される。標的によってスクリーニングされるいくつかの薬剤は、標的に結合するにすぎない。他の薬剤は、この標��を作動するかまたはそれに拮抗する。
【０１８４】
本発明の簡単な例として、タンパク質が、標的について改善された結合親和性を有するように、進化され得る。この方法は、その１次アミノ酸配列が既知であり、また、そのＸ線結晶解析に基づいた３次元構造が既知である、タンパク質の野生型形態または参照形態から開始する。このタンパク質は、その１次配列および３次構造が同様に既知であるタンパク質標的に結合することが知られている。このタンパク質および標的の相互作用は、上述のように運動方程式を解くことによって決定される。従って、この相互作用が評価され、野生型タンパク質と比較して、１以上のアミノ酸置換有するタンパク質の主要な接触残基およびその標的が決定される。次いで、運動方程式が、野生型と比較して１以上のアミノ酸置換を有する改変体にために再度解かれる。重要な接触が、野生型タンパク質の接触と比較される。さらなる接触および同じ接触に対するより短い結合距離の存在によって、より強い結合親和性が示唆される。換言すると、より少ない接触または長い結合距離の存在によって、より弱い結合親和性が示唆される。この過程が、さらなる改変体に対して反復される。次いで、その標的について最も強い親和性を有すると思われる、この改変体またはこの改変体のサブセットが合成され、実験的に試験される。
【０１８５】
別の例において、本発明の方法が抗体をヒト化するのに使用され得る。抗体は相補性決定領域（ＣＤＲ）を有し、このＣＤＲは、可変領域フレームワーク配列によって隔てられる結合の主要な原因であり得る。従来のヒト化過程において、ヒトアクセプター抗体および非ヒト（代表的にはマウス）ドナー抗体から開始される。目的は、ヒト抗体由来のフレームワーク領域と非ヒト抗体由来のＣＤＲを結合することである（Ｑｕｅｅｎら、Ｐｒｏｃ．Ｎａｔｌ．Ａｃａｄ．Ｓｃｉ．ＵＳＡ  ８６：１００２９〜１００３３（１９８９）およびＷＯ９０／０７８６１、ＵＳ５，６９３，７６２、ＵＳ５，６９３，７６１、ＵＳ５，５８５，０８９、ＵＳ５，５３０，１０１およびＷｉｎｔｅｒ、ＵＳ５，２２５，５３９（あらゆる目的のために、これらの全体が参考として援用される））。マウスＣＤＲ領域とヒト可変領域残基との非天然の近接した並置によって、非天然の配置的な拘束が生じ得、この拘束は、特定のアミノ酸残基の置換によって修正されなければ、結合親和性の消失を導く。置換のためのアミノ酸残基の選択は、コンピュータモデリングによって決定され得る。モデリングは、抗体単独の１次アミノ酸配列に基づいて実施され得るか、また、開始点として関連する抗体の鎖またはドメインについて解かれた構造を含み得る。この運動方程式は、抗体鎖について解かれ、３次元構造を決定する。このモデルは、どのフレームワークアミノ酸が最もＣＤＲ領域と緊密に相互作用するかということを示唆する。一般に、このモデルにおけるＣＤＲ領域の６Å以内にあるフレームワークのアミノ酸はＣＤＲ領域と相互作用すると考えられている。次いで、ヒトアクセプター抗体における対応するアミノ酸は、マウスドナー抗体由来の対応するアミノ酸で置換される。
【０１８６】
モデリングおよび、その標的との異なる化合物の相互作用の評価および比較の後に、スクリーニングされた化合物のうちの１つまたはそのサブセットが、合成および生化学的アッセイについて選択される。合成の性質は、化合物の性質に依存する。例えば、従来の有機化学、組換えＤＮＡ発現、固相ペプチド合成または固相合成は、この化合物に依存して使用され得る。次いで、それらの化合物は、標的との相互作用についてスクリーニングされ得る。いくつかの化合物が同時に試験される場合、そのアッセイはマイクロウェルプレートにおいて実施され得る。このアッセイは、その標的との、化合物の結合親和性または速度論を測定し得る。いくつかの方法において、このアッセイは、その標的と所望の様式で相互作用することが既知のコントロール化合物と競合する、その標的に対する化合物の結合特異性を測定する。いくつか方法において、このアッセイはその標的に対する化合物の触媒活性およびその逆を測定する。いくつかの方法において、この標的は、細胞性レセプターであり、そして、そのアッセイが、この化合物がレセプターを介してシグナルを変換する能力を測定する。いくつかの方法において、そのアッセイは、疾患動物モデル（例えば、ヒト疾患の症状を示すように設計されたトランスジェニック齧歯類）において実施される。この化合物の活性は、この齧歯類における疾患の症状の、阻止、軽減、または除去から決定される。次いで、インビトロ実験または動物実験において首尾良い結果を示す化合物は、ヒト臨床試験において試験されるか、さらなる誘導化合物設計の基礎としての役割を担い得る。臨床試験を通過した化合物は、臨床用途のために薬学的キャリアと共に処方される。この薬学的キャリアは、米国ＦＤＡまたは他国の類似の部局の十分な製造の慣例に従って製造される。非経口投与について、このキャリアは滅菌性であり、そして実質的に等張性である。
【０１８７】
従って、前述したことは、本発明の実施形態の完全な記述であるが、種々の実施形態、代替物または等価物が製造され得、そして使用されることは自明である。従って、上の記述は、本発明の範囲を限定するものとして扱うべきではなく、本発明の範囲は、添付の特許請求の範囲の境界および範囲によって規定される。
【図面の簡単な説明】
【図１】
図１は、本発明に従うソフトウエアシステムアーキテクチャのブロックモジュール図である。
【図２】
図２は、本発明に従う分子モデルの多体系のツリー構造を示す。
【図３】
図３は、図２の多体系の基準配置を示す。
【図４Ａ】
図４Ａは、図２の多体系における２つの物体間のスライダージョイントを示す。
【図４Ｂ】
図４Ｂは、図２の多体系における２つの物体間のピンジョイントを示す。
【図４Ｃ】
図４Ｃは、図２の多体系における２つの物体間のボールジョイントを示す。
【図５Ａ】
図５Ａは、陰的Ｅｕｌｅｒ積分法のＡ安定性試験およびＬ安定性試験である、安定性関数を示す。
【図５Ｂ】
図５Ｂは、陰的中点積分法のＡ安定性テストおよびＬ安定性テストである、安定性関数を示す。
【図５Ｃ】
図５Ｃは、Ｒａｄａｕ５積分法のＡ安定性テストおよびＬ安定性テストである安定性関数を示す。
【図６】
図６は、本発明の１つの実施形態に従う陰的Ｅｕｌｅｒ積分法の工程を示すフローチャートである。
【図７】
図７は、本発明の１つの実施形態に従うＲａｄａｕ５積分法の工程を示すフローチャートである。
【図８】
図８は、タンパク質フラグメントのアラニンジペプチドの分子構造の図である。
【図９Ａ】
図９Ａは、Ｖｅｒｌｅｔ積分法によって計算された場合の図８のアラニンジペプチドモデルに対する座標角度ψ  対  時間のプロットである。
【図９Ｂ】
図９Ｂは、Ｒａｄａｕ５積分法によって計算された場合の図８のアラニンジペプチドモデルに対する、座標角度ψ  対  時間のプロットである。
【図９Ｃ】
図９Ｃは、陰的Ｅｕｌｅｒ積分法によって計算された場合の図８のアラニンジペプチドモデルに対する、座標角度ψ  対  時間のプロットである。
【図９Ｄ】
図９Ｄは、Ｖｅｒｌｅｔ積分法によって計算された場合の図８のアラニンジペプチドモデルに対する、座標角度φ  対  時間のプロットである。
【図９Ｅ】
図９Ｅは、Ｒａｄａｕ５積分法によって計算された場合の図８のアラニンジペプチドモデルに対する、座標角度φ  対  時間のプロットである。
【図９Ｆ】
図９Ｆは、陰的Ｅｕｌｅｒ積分法によって計算された場合の図８のアラニンジペプチドモデルに対する、座標角度φ  対  時間のプロットである。
【図１０Ａ】
図１０Ａは、Ｖｅｒｌｅｔ積分法による、図９Ａおよび９Ｄのアラニンジペプチド座標シミュレーションについての、時間刻み幅  対  時間のプロットである。
【図１０Ｂ】
図１０Ｂは、Ｒａｄａｕ５積分法による図９Ｂおよび９Ｅのアラニンジペプチド座標シミュレーションについての、時間刻みの幅  対  時間のプロットである。
【図１０Ｃ】
図１０Ｃは、陰的Ｅｕｌｅｒ積分法による図９Ｃおよび９Ｆのアラニンジペプチド座標シミュレーションについての、時間刻みの幅  対  時間のプロットである。
【図１】
2004519026.tif 000079
【図５Ａ】
2004519026.tif 000080
【図５Ｂ】
2004519026.tif 000081
【図５Ｃ】
2004519026.tif 000082
【図６】
2004519026.tif 000083
【図７】
2004519026.tif 000084
【図９Ａ】
2004519026.tif 000085
【図９Ｂ】
2004519026.tif 000086
【図９Ｃ】
2004519026.tif 000087
【図９Ｄ】
2004519026.tif 000088
【図９Ｅ】
2004519026.tif 000089
【図９Ｆ】
2004519026.tif 000090
【図１０Ａ】
2004519026.tif 000091
【図１０Ｂ】
2004519026.tif 000092
【図１０Ｃ】
2004519026.tif 000093
【国際公開パンフレット】
2004519026.tif 000094
2004519026.tif 000095
2004519026.tif 000096
2004519026.tif 000097
2004519026.tif 000098
2004519026.tif 000099
2004519026.tif 000100
2004519026.tif 000101
2004519026.tif 000102
2004519026.tif 000103
2004519026.tif 000104
2004519026.tif 000105
2004519026.tif 000106
2004519026.tif 000107
2004519026.tif 000108
2004519026.tif 000109
2004519026.tif 000110
2004519026.tif 000111
2004519026.tif 000112
2004519026.tif 000113
2004519026.tif 000114
2004519026.tif 000115
2004519026.tif 000116
2004519026.tif 000117
2004519026.tif 000118
2004519026.tif 000119
2004519026.tif 000120
2004519026.tif 000121
2004519026.tif 000122
2004519026.tif 000123
2004519026.tif 000124
2004519026.tif 000125
2004519026.tif 000126
2004519026.tif 000127
2004519026.tif 000128
2004519026.tif 000129
2004519026.tif 000130
2004519026.tif 000131
2004519026.tif 000132
2004519026.tif 000133
2004519026.tif 000134
2004519026.tif 000135
2004519026.tif 000136
2004519026.tif 000137
2004519026.tif 000138
2004519026.tif 000139
2004519026.tif 000140
2004519026.tif 000141
2004519026.tif 000142
2004519026.tif 000143
2004519026.tif 000144
2004519026.tif 000145
2004519026.tif 000146
2004519026.tif 000147
2004519026.tif 000148
2004519026.tif 000149
【国際公開パンフレット（コレクトバージョン）】
2004519026.tif 000150
【国際公開パンフレット（コレクトバージョン）】
2004519026.tif 000151
2004519026.tif 000152
2004519026.tif 000153
2004519026.tif 000154
2004519026.tif 000155
2004519026.tif 000156
2004519026.tif 000157
2004519026.tif 000158
2004519026.tif 000159
2004519026.tif 000160
2004519026.tif 000161
2004519026.tif 000162
2004519026.tif 000163
2004519026.tif 000164
2004519026.tif 000165
2004519026.tif 000166
2004519026.tif 000167
2004519026.tif 000168
2004519026.tif 000169
2004519026.tif 000170
2004519026.tif 000171
2004519026.tif 000172
2004519026.tif 000173
2004519026.tif 000174
2004519026.tif 000175
2004519026.tif 000176
2004519026.tif 000177
2004519026.tif 000178
2004519026.tif 000179
2004519026.tif 000180
2004519026.tif 000181
2004519026.tif 000182
2004519026.tif 000183
2004519026.tif 000184
2004519026.tif 000185
2004519026.tif 000186
2004519026.tif 000187
2004519026.tif 000188
2004519026.tif 000189
2004519026.tif 000190
2004519026.tif 000191
2004519026.tif 000192
2004519026.tif 000193
2004519026.tif 000194
2004519026.tif 000195
2004519026.tif 000196
2004519026.tif 000197
2004519026.tif 000198
2004519026.tif 000199
2004519026.tif 000200
2004519026.tif 000201
2004519026.tif 000202
2004519026.tif 000203
【国際調査報告】
2004519026.tif 000204
2004519026.tif 000205
2004519026.tif 000206