メタダイナミクス

メタダイナミクス（MTD、略してMETADまたはMetaDとも呼ばれる）は、計算物理学、化学、生物学におけるコンピュータシミュレーション手法である。システムの自由エネルギーやその他の状態関数を推定するために使用される。システムのエネルギーランドスケープの形状によってエルゴード性が妨げられる。 2002年にAlessandro LaioとMichele Parrinelloによって初めて提案され^[1]、通常は分子動力学シミュレーション内に適用されます。MTDは、適応バイアス分子動力学^[2]、適応反応座標力^[3]、局所仰角アンブレラサンプリング[ ⁴ ]など、いくつかの新しい手法とよく似ています。最近では、元のメタダイナミクスとよく調整されたメタダイナミクス^[5]の両方が重要度サンプリングのコンテキストで導出され、適応バイアスポテンシャル設定の特殊なケースであることが示されました。^{[6] MTDは、} Wang–Landauサンプリングに関連しています。^[7]

この技術は、（時系列順に）デフレーション、^[8] トンネリング、^[9] タブーサーチ、^[10] 局所標高、^[11] コンフォメーションフラッディング、^[12] エンクヴィスト・カールストローム^[13]および 適応バイアス力法を含む多数の関連手法に基づいています。^[14]

メタダイナミクスは非公式に「自由エネルギーの井戸を計算砂で埋める」と説明されている。^[15]このアルゴリズムは、システムが少数の集団変数(CV) で記述できると想定している。シミュレーション中、集団変数によって決定される空間におけるシステムの位置が計算され、システムの実際のエネルギーランドスケープに正のガウスポテンシャルが追加される。このようにして、システムは前のポイントに戻ることを阻止される。シミュレーションの展開中、より多くのガウスが合計され、システムが前のステップに戻ることを阻止する度合いがますます高くなり、システムが完全なエネルギーランドスケープを探索するまで続く。この時点で、修正された自由エネルギーは集団変数の関数として定数になり、これが集団変数が大きく変動し始める理由である。この時点で、エネルギーランドスケープは、すべてのガウスの合計の逆として回復できる。

2 つのガウス関数の加算間隔、およびガウス関数の高さと幅は、精度と計算コストの比率を最適化するために調整されます。ガウス関数のサイズを変更するだけで、メタダイナミクスをフィッティングして、大きなガウス関数を使用してエネルギーランドスケープの大まかなマップを非常に迅速に生成したり、小さなガウス関数を使用してより細かい粒度の記述に使用したりできます。^[1]通常、ガウス関数のサイズを適応的に変更するには、よく調整されたメタダイナミクス^[5]が使用されます。また、ガウス関数の幅は、適応型ガウス関数メタダイナミクスを使用して適応できます。^[16]

メタダイナミクスは、適応アンブレラサンプリングなどの方法に比べて、探索するエネルギーランドスケープの初期推定を必要としないという利点がある。 ^[1]しかし、複雑なシミュレーションに適切な集団変数を選択することは簡単ではない。通常、適切な集団変数のセットを見つけるには数回の試行が必要であるが、いくつかの自動手順が提案されている：必須座標、^[17]スケッチマップ、^[18]および非線形データ駆動型集団変数。^[19]

独立したメタダイナミクスシミュレーション（レプリカ）を結合して、使いやすさと並列パフォーマンスを向上させることができます。提案されている方法はいくつかあります。マルチウォーカーMTD、^[20] 、並列テンパリングMTD、^[21]、バイアス交換MTD、^[22]、集団変数テンパリングMTD ^[23]です。最後の3つは並列テンパリング法に似ており、レプリカ交換を使用してサンプリングを改善します。通常、メトロポリス-ヘイスティングスアルゴリズムがレプリカ交換に使用されますが、無限スワッピング^[24]とSuwa-Todo ^[25]アルゴリズムを使用すると、レプリカ交換率が向上します。^[26]

典型的な（単一レプリカ）MTDシミュレーションには最大3つのCVを含めることができますが、マルチレプリカアプローチを使用しても、実際には8つのCVを超えることは困難です。この制限は、ガウス関数（カーネル）を追加することによって構築されるバイアスポテンシャルに起因します。これは、カーネル密度推定器（KDE）の特殊なケースです。一定のKDE精度に必要なカーネルの数は、次元数とともに指数関数的に増加します。したがって、バイアスポテンシャルの同じ精度を維持するには、MTDシミュレーションの長さをCVの数とともに指数関数的に増やす必要があります。また、バイアスポテンシャルは、高速評価のために、通常は規則的なグリッドで近似されます。^[27]グリッドを格納するために必要なメモリも、次元（CV）の数とともに指数関数的に増加します。

メタダイナミクスの高次元一般化はNN2Bである。^[28]これは、2つの機械学習アルゴリズム、すなわち最近傍密度推定器（NNDE）と人工ニューラルネットワーク（ANN）に基づいている。NNDEはKDEに代わるもので、短いバイアスシミュレーションからバイアスポテンシャルの更新を推定し、ANNは結果として得られるバイアスポテンシャルを近似するために使用される。ANNは、高次元関数のメモリ効率の高い表現であり、微分（バイアス力）はバックプロパゲーションアルゴリズムで効果的に計算される。^{[28] [29]}

ANNを適応バイアスポテンシャルに利用する代替法では、推定に平均ポテンシャル力を使用します。 ^{[30]この方法も、}適応バイアス力（ABF）法の高次元一般化です。 ^[31]さらに、ANNのトレーニングはベイズ正則化を使用して改善され、^[32] ANNのアンサンブルをトレーニングすることで近似誤差を推測できます。^[30]

2015年にホワイト、ダマ、ヴォスは、分子動力学シミュレーションを望ましい自由エネルギー面に合わせて形作る方法である実験指向メタダイナミクスを導入しました。この技術は、シミュレーションを実験データと一致するコンフォメーションに導き、複雑な分子システムとその挙動に対する理解を深めます。^[33]

2020年には、メタダイナミクスの進化形であるオンザフライ確率強化サンプリング法（OPES）が提案され、^{[34] [35] [36]、現在では}Michele Parrinelloの研究グループが採用している。 ^[37] OPES法は、堅牢なパラメータがわずかしかなく、メタダイナミクスよりも収束が速く、再重み付けスキームが単純である。^[38] 2024年には、OPESのレプリカ交換バリアントであるOneOPESが開発され、^[39]熱勾配と複数のCVを利用して、複数の自由度を持つ大規模な生化学システムをサンプリングするように設計されています。このバリアントは、多数の自由度を少数のCVで捉えることが困難な場合が多いこのようなシステムを記述するという課題に対処することを目的としています。OPESは、バージョン2.7以降、PLUMEDライブラリに実装されています。 ^[40]

直交座標の位置にある古典的な - 粒子系があると仮定します。粒子の相互作用は、ポテンシャル関数で記述されます。ポテンシャル関数形式 (たとえば、高エネルギー障壁によって分離された 2 つの極小値) は、分子動力学法やモンテカルロ法によるエルゴードサンプリングを防止します。 ${\textstyle N}$${\textstyle \{{\vec {r}}_{i}\}}$ ${\textstyle (i\in 1...N)}$ ${\textstyle ({\vec {r}}_{i}\in \mathbb {R} ^{3})}$${\textstyle V\equiv V(\{{\vec {r}}_{i}\})}$

MTD の一般的な考え方は、サンプリングされた状態の再訪問を阻止することでシステム サンプリングを強化することです。これは、システムハミルトニアン にバイアス ポテンシャルを追加することで実現されます。 ${\textstyle H}$${\displaystyle V_{\text{バイアス}}}$

${\displaystyle H=T+V+V_{\text{バイアス}}}$。

バイアス ポテンシャルは集団変数 の関数です。集団変数は粒子の位置の関数です。バイアス ポテンシャルは、バイアスを速度 で追加することによって継続的に更新されます。ここで、 は 時間 における瞬間的な集団変数値です。 ${\textstyle (V_{\text{バイアス}}\equiv V_{\text{バイアス}}({\vec {s}}\,))}$${\displaystyle ({\vec {s}}\equiv {\vec {s}}(\{{\vec {r}}_{i}\}))}$${\displaystyle \omega}$${\displaystyle {\vec {s}}_{t}}$${\displaystyle t}$

${\displaystyle {\frac {\partial V_{\text{bias}}({\vec {s}}\,)}{\partial t}}=\omega \,\delta (|{\vec {s}}-{\vec {s}}_{t}|)}$。

無限に長いシミュレーション時間では、蓄積されたバイアスポテンシャルは反対の符号（および無関係な定数）を持つ自由エネルギーに収束します。 ${\displaystyle t_{\text{sim}}}$${\displaystyle C}$

${\displaystyle V_{\text{bias}}({\vec {s}}\,)=\!\!\int _{0}^{t_{\text{sim}}}\!\!\!\omega \,\delta (|{\vec {s}}-{\vec {s}}_{t}|)\;dt\quad \Rightarrow \quad F({\vec {s}}\,)=-\!\!\!\!\lim _{t_{\text{sim}}\to \infty }\!\!V_{\text{bias}}({\vec {s}}\,)+C}$

計算効率の高い実装では、更新プロセスは時間間隔（床関数を表す）に離散化され、関数は局所的な正のカーネル関数に置き換えられます。バイアスポテンシャルは、時間における瞬間的な集団変数値を中心とするカーネル関数の合計になります。 ${\displaystyle \tau }$${\displaystyle \lfloor \;\rfloor }$${\displaystyle \delta }$ ${\displaystyle K}$${\displaystyle {\vec {s}}_{j}}$${\displaystyle \tau j}$

${\displaystyle V_{\text{bias}}({\vec {s}}\,)\approx \tau \!\!\!\sum _{j=0}^{\left\lfloor {\frac {t_{\text{sim}}}{\tau }}\right\rfloor }\!\!\omega \,K(|{\vec {s}}-{\vec {s}}_{j}|)}$。

通常、カーネルは多次元ガウス関数であり、その共分散行列には対角の非ゼロ要素のみが含まれます。

${\displaystyle V_{\text{bias}}({\vec {s}}\,)\approx \tau \!\!\!\sum _{j=0}^{\left\lfloor {\frac {t_{\text{sim}}}{\tau }}\right\rfloor }\!\!\omega \exp \!\!\left(\!-{\frac {1}{2}}\left|{\frac {{\vec {s}}-{\vec {s}}_{j}}{\vec {\sigma }}}\right|^{2}\right)}$。

パラメータ、、は事前に決定され、シミュレーション中は一定に保たれます。 ${\displaystyle \tau }$${\displaystyle \omega }$${\displaystyle {\vec {\sigma }}}$

以下は分子動力学（MD）に基づく MTD の擬似コードです。ここで、 とはそれぞれ -粒子系の位置と速度です。バイアスはMD ステップごとに更新され、システムの力への寄与はです。 ${\displaystyle \{{\vec {r}}\}}$${\displaystyle \{{\vec {v}}\}}$${\displaystyle N}$${\displaystyle V_{\text{bias}}}$${\displaystyle n=\tau /\Delta t}$${\displaystyle \{{\vec {F}}\,\}}$${\displaystyle \{{\vec {F}}_{\text{bias}}\}}$

初期値を設定し、設定する${\displaystyle \{{\vec {r}}\}}$${\displaystyle \{{\vec {v}}\}}$ 
 ${\displaystyle V_{\text{bias}}({\vec {s}}\,):=0}$

各MDステップ:
     CV値
        を計算:${\displaystyle {\vec {s}}_{t}:={\vec {s}}(\{{\vec {r}}\})}$
    
    各 MDステップ:
        バイアスポテンシャル
            を更新:${\displaystyle n}$${\displaystyle V_{\text{bias}}({\vec {s}}\,):=V_{\text{bias}}({\vec {s}}\,)+\tau \omega \exp \!\!\left(\!-{\frac {1}{2}}\left|{\frac {{\vec {s}}-{\vec {s}}_{t}}{\vec {\sigma }}}\right|^{2}\right)}$
    
    原子力
        を計算する:${\displaystyle {\vec {F}}_{i}:=-{\frac {\partial V(\{{\vec {r}}\,\})}{\partial {\vec {r}}_{i}}}\overbrace {\left.-{\frac {\partial V_{\text{bias}}({\vec {s}}\,)}{\partial {\vec {s}}}}\right|_{{\vec {s}}_{t}}\!\!\!{\frac {\partial {\vec {s}}(\{{\vec {r}}\,\})}{\partial {\vec {r}}_{i}}}} ^{{\vec {F}}_{{\text{bias}},i}}}$
    
    伝播 し、${\displaystyle \{{\vec {r}}\}}$${\displaystyle \{{\vec {v}}\}}$${\displaystyle \Delta t}$

カーネルの有限サイズにより、バイアスポテンシャルは平均値の周りで変動します。収束した自由エネルギーは、バイアスポテンシャルを平均化することで得られます。平均化は、集団変数に沿った動きが拡散的になる から始まります。 ${\displaystyle t_{\text{diff}}}$

${\displaystyle {\bar {F}}({\vec {s}})=-{\frac {1}{t_{\text{sim}}-t_{\text{diff}}}}\int _{t_{\text{diff}}}^{t_{\text{sim}}}\!\!\!\!\!V_{\text{bias}}({\vec {s}},t)\,dt+C}$

メタダイナミクスは以下の研究に使用されています:

• タンパク質の折り畳み^[22]
• 化学反応^[41]
• 分子ドッキング^{[42] [43]}
• 相転移^[44 ]
• ^{疎水性[45]}および親水性^[46]単層カーボンナノチューブへのDNAのカプセル化。

PLUMED ^[47]は、多くのMTDアルゴリズムと集合変数を実装したオープンソース ライブラリです。柔軟なオブジェクト指向設計^{[48] [49]を採用しており、いくつかのMDプログラム（} AMBER、GROMACS、LAMMPS、NAMD、Quantum ESPRESSO、DL_POLY_4、CP2K、OpenMM）とインターフェースすることができます。 ^{[50] [51]}

他のMTD実装としては、Collective Variables Module ^[52] ( LAMMPS、NAMD、GROMACS用)、ORAC、CP2K、^[53] EDM、^[54] Desmondなどがあります。

• メタダイナミクス入門
• 羽毛
• Colvars モジュール Web サイト (NAMD、LAMMPS、Gromacs)
• メタダイナミクスのビジュアルムービー
• オンザフライ確率強化サンプリング (OPES)

• ローカル標高
• 平行焼戻し
• アンブレラサンプリング