Sampling Challenges of MM/PBSA Binding Energy Calculations
Xu, Zhou, Zheng, Wang, Peng*, Li* — J. Phys. Chem. B 2025 | DOI: 10.1021/acs.jpcb.5c04908
⚠️ MM/PBSA精度はサンプリング充足性が決め手。短時間MDの「擬似収束」に要注意 — 570μs大規模検証で実証
① 背景と課題

MM/PBSA法は計算コストとFEP精度の中間に位置する主力ツールだが、MDサンプリングの充足性が精度を根本的に左右する。従来研究は「長いMDほど良い」と信じていたが、実際には：

短時間MDが「擬似収束」を示し、誤った安心感を与えることがある
逆に長時間MDでも力場誤差が顕在化して実験値との乖離が広がる場合がある
遅い構造遷移（μsスケールのリガンドフリップ等）が収束を妨げる本質的原因

→ 4タンパク質×19複合体×3独立×10μs = 総570μs の超大規模検証でこの問題を解明

② 系統別収束結果
タンパク質収束状況最大差異
TNKS2✅ 全5系収束0.19〜1.13
PLpro⚠️ リガンド依存〜数kcal/mol
HIF-2α⚠️ リガンド依存〜数kcal/mol
c-Met❌ 全5系非収束最大12.9

単位: kcal/mol（10μs終点での3トラジェクトリ間最大差異）

570 μs
総シミュレーション時間（前例のない規模）
③ リガンド遅動変化の3分類
🔄 回転可能結合周りの内部回転
（hif2a-4 など）
↓
↔️ ポケット内全体的平行移動
（cmet-11 など）
↓
🔁 ポケット内リガンド完全フリップ
（plpro-7SQE: 4μs継続）

フリップ発生系でMM/PBSA値が急激に変動

PLpro-7SQEではリガンドが4μsにわたり反転状態を維持。この遅い動的変化が短時間MDでは捕捉不能。

④ PCA-based鍵相互作用同定

受容体-リガンド間原子距離（カットオフ8Å）をPCAにかけ、PC2投影とMM/PBSAエネルギーのPearson相関 = 0.73 を確認。

E166
H-bond形成時に結合エネルギー低下。存在/不在で収束困難の主因。
Y170
E166と水架橋を形成し間接的にH-bond形成を調節するアロステリック効果。
Y267
π-π蓋形成。開閉状態でエネルギーが変動。
自動同定フロー
PCA → PC2 → 上位ウェイト残基 → H-bond/π-π解析
⑤ 相互作用統計（19系合計）
相互作用種別同定数結合に寄与
H-bond108本77本（71%）
Salt bridge4本全て寄与
π-π interaction60本大半が寄与
11本のH-bond: H-bond形成時に結合エネルギーが逆に上昇（競合H-bondの置き換えによる）

水架橋（water bridge）: 短時間MDの高占有率が長時間MDで低下 → サンプリング不足の追加指標

⑥ lib/fep への実装提案
  • MMGBSAEngineにrun_replicates(n=3)メソッド追加 — 複数トラジェクトリ収束診断
  • 累積平均最大差異（kcal/mol）を収束指標として自動レポート
  • PCA-based鍵残基自動同定クラス実装（MDAnalysis + sklearn）
  • GaMD/REST2前処理ラッパー — 増強サンプリング統合
  • HBondAnalyzerへwater_bridge_occupancy()メソッド追加
実装優先度: HIGH
  • 収束診断なしのMM/PBSA結果は信頼性に疑問符
  • 短時間MDで「良い結果」でも擬似収束の可能性を自動警告
    • J. Phys. Chem. B 2025, 129, 11666–11678 | DOI: 10.1021/acs.jpcb.5c04908 | Xu, Zhou, Zheng, Wang, Peng, Li | 計算化学 | lib/fep, lib/md