Quantifying the uncertainty of molecular dynamics simulations : Good-Turing statistics revisited

技術要約：グッド・チューリング統計を用いた分子動力学シミュレーションの不確実性の定量化

問題の定義
生体分子系の分子動力学（MD）シミュレーションは計算負荷が高く、マクロ分子が取り得るすべての構造を完全にサンプリングすることは事実上不可能です。すべてのアクセス可能な構造を忠実にサンプリングすることは必ずしも必要ではありませんが、研究者には、もしシミュレーションを延長した場合に、全く新しい（未観測の）構造が観察される確率を推定する手法が求められます。著者らは以前、グッド・チューリング統計をMDトラジェトリに適用してこの不確実性を推定できることを確立しました。しかし、初期の実装には重大なスケーラビリティの制限がありました。それは、全トラジェトリにわたる二次元のルート平均二乗偏差（RMSD）行列の計算と保存を必要とすることでした。このメモリ要件は構造数の二乗に比例して増加するため、非常に長いシミュレーション（例：数百万フレームを含むもの）への適用を阻害していました。

手法
著者らは、メモリ要件を二次スケーリングから線形スケーリングへと削減し、極めて長いシミュレーション（最大2,200万構造）への適用を可能にする、新しい変法としてのグッド・チューリング・アルゴリズムを提案しています。

• 核となる概念： この手法は、これまでに観察されたすべての構造から $x$ Å RMSD 以上離れた新しい構造が観察される確率（$P$）を推定します。グッド・チューリングの頻度推定の文脈では、$P = N_1 / N$ となります。ここで、$N$ は総既知のコンフォメーション数、$N_1$ は正確に一度だけ観察されたコンフォメーション数です。

• 新しいアルゴリズム： 全体の2D RMSD行列を構築してデンドログラムを生成するための階層的クラスタリングを行う代わりに、この新しいアプローチは行列自体を回避します。

  1. 行単位の処理： トラジェトリ内の各構造（サブサンプリング係数 $s$ を介して選択）について、その参照構造と他のすべての構造との間のRMSDを計算します。
  2. 最大化： その特定の参照構造に対する「最大」RMSD値のみを保存します。
  3. ソート： これら（各参照構造につき1つ）の最大RMSD値のリストを降順にソートします。
  4. 確率の割り当て： ソートされたリストを確率にマッピングします。最大のRMSDは $P = 1/N$ に対応し、2番目に大きいものは $P = 2/N$ に対応し、以下同様です（ここで $i$ はソートされたリストのインデックス）。これにより、$P$ 対 RMSD の曲線が生成されます。
  5. メモリ効率： このプロセスは、完全な $N \times N$ 行列を保持するのではなく、構造ごとに1つの数値（最大RMSD）のみを保持することを要求します。これは、物理的なメモリ効率のために実行速度をトレードオフとするものです。

• サブサンプリング係数（$s$）の決定： 連続する構造がメカニズム的に相関しないようにするために、サブサンプリング係数 $s$ を決定することは極めて重要なステップです。著者らは、時間間隔（$\delta t$）に対する最大RMSDの「成長曲線」に対して、区分線形回帰（Rパッケージ dpseg を使用）を用いるヘリスティックな手法を採用しています。アルゴリズムは、RMSDが安定したプラトーに収束する点を見出し、不確実性を過小評価することを避けるために、観察された構造変化と矛盾しない最も長いタイムスケールを選択します。

主な貢献

1. 線形メモリ・スケーリング： 主な貢献は、メモリ制約により以前は不可能であった、最大2,200万構造のトラジェトリに対するグッド・チューリング解析を可能にするアルゴリズムの再定式化です。
2. 等価性の検証： 著者らは、この新しいメモリ効率の高い手法が、従来の行列ベースの実装と本質的に同一の結果を与えることを示しています。
3. ソフトウェアの可用性： 新しいアルゴリズムを実装したコンピュータプログラムが、標準的なリポジトリ（GitHub）を通じてオープンソースソフトウェアとして提供されています。

結果
著者らは、安定した折り畳みタンパク質（ROP）、ミニタンパク質（FipWWドメイン）、およびペプチド（CLN025, 6NM2）を含む様々な系をカバーする、6.6から20 $\mu$sの範囲のトラジェトリにこの新手法を適用しました。

• 一貫性： ROPおよびCLN025のシミュレーションにおける、古典的な（2D行列）手法と新しい（行最大値）手法との直接比較において、優れた一致が見られました。
• 不確実性の定量化： 得られた確率曲線は、構造的不確実性を定量化することに成功しました。例えば、安定したROPタンパク質のシミュレーション時間を2倍にすると、約0.95 Å RMSD以内で異なる新しい構造が得られると推定されました。対照的に、34残基のFipWWタンパク質の折り畳みシミュレーションでは、時間を2倍にすると最大で約11 Å RMSD離れた構造が得られると予測され、折り畳みの広大なコンフィギュレーショナル空間を捉えています。
• ダイナミクスの識別： 曲線は、ネイティブ状態への再訪問が高い確率を示す安定したフォルダー（CLN025）と、より大きなRMSDにおいて有意に異なる構造を示す柔軟または無秩序なペプチド（6NM2）を区別しました。

意義と限界
本論文は、この手法がMDシミュレーションの不確実性を定量化するための「信頼でき、安定し、かつ検証可能な」方法を提供すると主張しています。これにより、研究者は「シミュレーション時間を2倍にした場合、期待される最も異なる新奇構造のRMSDはいくらか？」といった実用的な問いに答えることができます。この能力は、MDシミュレーションから導かれる結論の妥当性を強固にするための重要なステップとして提示されています。

しかし、著者らは根本的な限界についても明示的に認めています。グッド・チューリング統計は、未知のサイズのプールから個別のオブジェクトをサンプリングする場合にのみ厳密に有効です。MDトラジェトリは連続的であり、本手法はそれらを離散化するために単一のサブサンプリング係数（$s$）の選択に依存しています。著者らは、トラジェトリに混合されたタイムスケール（例：速い側鎖の揺らぎと遅い折り畳みイベント）が含まれている場合、アルゴリズムが単一のプラトーを特定することに苦慮する可能性があると指摘しています。デフォルトでは、アルゴリズムは不確実性を過小評価することを避けるために最も長い観察タイムスケールを選択しますが、これは、シミュレーション設計が異質なタイムスケールを混合している場合、すべての動的な詳細を捉えきれない可能性があることを意味します。ノイズや複雑なデータにおいて収束点をアルゴリズム的に定義することの難しさから、サンプリング係数の決定は、この手順の中で技術的に最も脆弱な部分となっています。