GROMACS+MM/PBSA计算蛋白-配体结合自由能——一个药物筛选项目的能量分解实录

这个项目的目标是筛选一批候选化合物对靶标激酶的亲和力。分子对接粗筛出了12个hits，但对接的打分函数在高通量下的假阳性率高得吓人——需要MD+MM/PBSA做一层精筛。

MM/PBSA的方法论看起来直截了当：跑一段MD，把轨迹按时间帧切片，每帧做ΔG = G_complex − G_protein − G_ligand，然后取平均。但真正上手跑的时候，几个决策点的代价远超预期。

决策点一：MD轨迹取多长？

对接pose的MD跑了50 ns，前20 ns拿出来做MM/PBSA计算。这个”前20 ns”的选择不是拍脑袋——头5 ns的RMSD曲线还在剧烈上升，蛋白和配体都在做结构调整。如果一个药物分子连20 ns的稳定结合都做不到，那后续的体外实验也大概率没有活性。

但项目里一个细节值得留意：前20 ns和后20 ns（30-50 ns）的ΔG值差了接近3 kcal/mol。后20 ns算出来的结合更强——因为配体在30 ns后有一次明显的构象调整，一个苯环从溶剂暴露区翻到了疏水口袋深处，ΔG的变化主要来自溶剂可及表面积（SASA）项的非极性贡献。

这不是说后20 ns”更准”——而是说你取哪一段，取决于你相信结合模式在哪个时段达到了稳定。项目最后的判断是取20-50 ns，因为前20 ns包含了太多非平衡态构象。但如果你做的是共价抑制剂的结合，可能恰恰要取前20 ns——成键之前的pose才是你关心的。

决策点二：介电常数取多少？

MM/PBSA公式里，极性溶剂化能是通过解Poisson-Boltzmann方程得到的，这个计算强烈依赖于溶质和溶剂的介电常数设置。默认设ε_in=1、ε_out=80——这是在模拟”蛋白内部=真空”和”水溶液=连续介质”。

问题是蛋白内部不是真空。蛋白骨架和侧链的极化效应在ε_in=1时被严重低估。这个项目对比了三组：

ε_in=2时ΔG=-14.7 kcal/mol，与等温滴定量热（ITC）实验的Ki换算值（ΔG_exp=-16.3 kcal/mol）偏差最小。ε_in=1把静电力夸大了近50%，导致结合被系统性高估——这在极性配体（含羧基、磷酸基团）上尤其严重。

这个ε_in=2的取值并非通用结论。它来自J. Chem. Inf. Model. 2019年一篇大量基准测试的系统验证——在对200+蛋白-配体体系的对比中，ε_in=2-2.5是让MM/PBSA结果最接近实验值的区间。但如果你的配体是高度疏水的，ε_in=1反而更准——因为非极性结合的主导贡献是范德华和SASA，极性项下调对整体影响很小。

决策点三：熵变要不要算？

MM/PBSA标准protocol包含熵变项（-TΔS），用简正模分析（normal mode analysis）计算。但这个计算的计算量极大——100帧轨迹做NMA在48核CPU上跑了将近四个小时——而且结果的统计噪音常常比信号还大。

这个项目的ΔS计算结果：-12.4 ± 8.7 kcal/mol。±8.7是什么概念？结合自由能ΔG_total才-14.7。也就是说熵变项的误差条是信号的两倍宽。

能不能不算熵变直接比较ΔH？在某些场景下可以——如果比对的是同一靶标的同一结合口袋，配体结构变化不大，熵变贡献趋于抵消。但在这个项目里，12个hits中有三个柔性链结构差异很大，熵变贡献不可忽略。

项目最后的处理方式：保留熵变项，但用ΔG（含熵）和ΔH（不含熵）做双重排序。只有在两者一致给出的top 5 hits中才推进后续实验。这个策略牺牲了灵敏度，换来了鲁棒性。

回过头看

MM/PBSA不是拿到ΔG数值就完事了。它的价值在于能量分解——你可以把结合自由能拆成范德华、静电、极性溶剂化、非极性溶剂化、熵五部分，哪个项是结合的主要驱动力，一拆就知道。

在这个项目里，top 1化合物结合的主导贡献是非极性溶剂化（占比约55%），说明结合模式以疏水匹配为主——这个信息比ΔG绝对值更有价值，直接指导了下一步的骨架优化：在保持疏水匹配的前提下增强氢键网络，提高结合选择性。

引用来源：Miller et al., JCTC 2012（MM/PBSA协议标准）；Sun et al., J. Chem. Inf. Model. 2019（ε_in基准测试）；gmx_MMPBSA官方文档 v1.5+。