MM/PBSA：从轨迹估算结合自由能

MM/PBSA（Molecular Mechanics / Poisson-Boltzmann Surface Area）是介于对接打分与 FEP 之间的自由能估算方法：比对接分数准，比 FEP 便宜。理解它的公式构成与近似，才能知道结果能信到什么程度。

基本公式

# 结合自由能的分解：
#
#   ΔG_bind = ΔH - TΔS
#           ≈ ΔE_MM + ΔG_solv - TΔS
#
# 其中：
#   ΔE_MM  = ΔE_内部 + ΔE_静电 + ΔE_范德华
#            （分子力学能，用力场算）
#
#   ΔG_solv = ΔG_极性溶剂化 + ΔG_非极性溶剂化
#            极性部分：解泊松-玻尔兹曼方程（PB）
#                      或广义玻恩近似（GB，见 126）
#            非极性部分：通常与溶剂可及表面积成正比
#                      ΔG_nonpolar = γ · ΔSASA + β
#
#   -TΔS   = 构象熵损失
#            【最难算、最常被省略的一项】
#
# 单轨迹近似（single-trajectory approach）：
#   只跑复合物的 MD，
#   从轨迹中分别提取复合物、蛋白、配体的构象
#   → 省去两次额外的模拟
#   → 【隐含假设：结合前后各组分构象不变】
#   → 对刚性体系合理，对柔性体系有偏差
#
#   优点：误差部分抵消，结果通常更稳定
#   缺点：忽略了构象重组能

熵项：最大的近似

• 为什么难算：构象熵需要对整个构象空间采样并统计，MD 的采样往往远不充分；
• 常见做法：
  • 直接省略——这是最常见的做法。此时得到的不是自由能，而是「结合焓的估计」；
  • 简正模式分析：计算量极大（对每一帧做 Hessian 对角化），且基于谐振近似，误差大；
  • 准谐近似：从轨迹的协方差矩阵估计，比简正模式便宜但仍不准；
  • 相互作用熵：从静电与范德华能量的涨落估计，较新的方法。
• 省略熵项的后果：绝对值完全不可信（通常严重高估结合强度）。但如果比较的分子熵贡献相近（同系列类似物），相对排序仍可能有用——这正是 MM/PBSA 的正确定位；
• 实践建议：不要报告绝对的 ΔG_bind，只用相对值做排序，并明确说明是否包含熵项。

PB 与 GB 的选择

	PB（泊松-玻尔兹曼）	GB（广义玻恩）
原理	数值求解连续介质静电方程	解析近似
精度	更高	较低
速度	慢（每帧需解方程）	快 10～100 倍
适用	精细比较，帧数少	大量分子的筛选
对带电体系	更可靠	误差较大

实践中 GB 用得更多（见 126《MM/GBSA》），因为速度差异太大；PB 用于对关键分子的仔细比较。

用 gmx_MMPBSA 计算

# gmx_MMPBSA 支持 GROMACS/AMBER/NAMD 的轨迹
pip install gmx_MMPBSA

# 输入文件 mmpbsa.in:
#   &general
#     sys_name="complex",
#     startframe=100,      # 跳过平衡阶段
#     endframe=500,
#     interval=2,          # 每 2 帧取一次（减少相关性）
#     verbose=2,
#   /
#   &pb
#     istrng=0.15,         # 离子强度 0.15 M
#     inp=1,
#     radiopt=0,           # 用力场半径
#   /
#   &gb
#     igb=5,               # GB 模型（2、5、8 常用）
#     saltcon=0.15,
#   /
#   # 逐残基分解（很有用）
#   &decomp
#     idecomp=2,
#     dec_verbose=0,
#   /

gmx_MMPBSA -O -i mmpbsa.in \
  -cs complex.tpr -ci index.ndx -cg 1 13 \
  -ct trajectory.xtc \
  -o FINAL_RESULTS.dat \
  -eo FINAL_RESULTS.csv

# 关键参数：
#   startframe   【必须跳过平衡阶段】
#   interval     取样间隔，避免相邻帧高度相关
#   igb          GB 模型的选择会显著影响结果
#   istrng       离子强度应匹配实验条件

逐残基能量分解：最有实用价值的输出

# idecomp=2 会给出每个残基对结合能的贡献
#
# 输出类似：
#   Residue      Internal  vdW      Elec     Polar    Nonpolar  Total
#   LEU:145      0.00      -2.34    -0.12    0.45     -0.28     -2.29
#   ASP:86       0.00      -0.87    -8.92    6.71     -0.15     -3.23
#   PHE:80       0.00      -3.11    -0.34    0.52     -0.31     -3.24
#   ...
#
# 用途：
#   1) 识别贡献最大的残基（热点残基）
#      → 设计时应保住与这些残基的相互作用
#
#   2) 识别【不利】的残基
#      Total 为正 = 该残基对结合不利
#      → 这是明确的改造机会
#
#   3) 解释选择性
#      比较靶点与脱靶的残基贡献差异
#
#   4) 指导突变实验设计
#      预测哪些突变会显著影响结合
#      → 可用实验验证计算的合理性
#
# 【注意】：
#   残基分解的可靠性低于总能量
#   （分解方式本身有任意性）
#   → 用于产生假设，不作为定论

# 解析结果
import pandas as pd
df = pd.read_csv("FINAL_RESULTS.csv", skiprows=...)
hotspots = df[df["Total"] < -1.0].sort_values("Total")
print("热点残基:")
print(hotspots[["Residue", "Total"]].head(10))

准确性的现实

用途	可靠性
绝对结合自由能	不可靠（误差常达数 kcal/mol）
同系列类似物的相对排序	中等——主要用途
跨骨架比较	差
带电分子的比较	差（GB 对电荷敏感）
识别热点残基	中等——有实用价值
优于对接分数	通常是的（见 095《Docking Score》）
不如 FEP	明确不如（见 127《FEP、RBFE 与 ABFE》）

提高可靠性的做法

• 用多个独立副本：3 个模拟各算一次，看结果分散度——这是评估误差最实际的方法；
• 充分平衡后再取样：前 10～20 ns 通常应丢弃；
• 取样间隔要够大：相邻帧高度相关，等效样本数远小于帧数；
• 报告标准误而非标准差：并说明基于多少个独立副本；
• 在已知数据上校准：用该靶点已知活性的分子测试排序能力，得到该体系上的可信度；
• 避免比较电荷不同的分子：净电荷不同时，溶剂化能的误差会主导结果。

延伸资源

• MM/GBSA：126《MM/GBSA》；FEP：127《FEP、RBFE 与 ABFE》；MD 入门：123《分子动力学 MD 入门》；
• 打分函数：095《Docking Score》；力场参数化：202《OpenFF Toolkit》。