蛋白质-配体结合自由能的MM/PBSA计算中采样不足如何影响结果

MM/PBSA在药物设计里用得非常多——把一个蛋白-配体复合物的MD轨迹拿出来，对每个帧算一次分子力学能+溶剂化自由能，取平均，复合物的值减去蛋白和配体各自的值就是结合自由能。流程成熟、工具成熟（gmx_MMPBSA脚本几乎一站式搞定），但真正让人头疼的不是”怎么算”，是”跑多久才算够”。

一个激酶（约320个残基）与小分子抑制剂（分子量420Da）的体系，实验Kd值在纳摩尔级（对应的ΔG约-12 kcal/mol）。我们用三种不同的模拟时长跑了一组对比——10ns、50ns、200ns——同样的力场（AMBER ff14SB蛋白+GAFF2配体）、同样的溶剂条件（TIP3P水、0.15M NaCl），唯一的变量就是采样时间。出来的结果差异远远超出了预期。

10ns：看起来收敛，实则陷阱

10ns的轨迹在RMSD上已经平了——蛋白Cα的RMSD在8ns之后稳定在1.8±0.2Å，配体重原子的RMSD也收敛到了0.5Å以内。如果你只检查轨迹的RMSD稳定性，10ns看起来足够了。

但这个”够了”是假象。配体在结合口袋里的构象采样远没达到平衡——10ns里配体的可旋转键（共6个）中有4个没有完成完整的构象翻转。这些翻转的自由能垒大概在3-5 kcal/mol之间，在300K下需要几十纳秒才能被充分采样。缺少翻转意味着配体的熵贡献被严重低估——MM/PBSA的准简谐近似对熵的估计高度依赖构象采样的完备性，缺了四个旋转自由度的贡献，计算出来的结合自由能比实验值浅了将近7 kcal/mol（模拟值-5.2，实验值约-12.0）。

不止是熵的问题。结合口袋里的几个关键水分子——尤其是与配体形成氢键网络的那两个桥接水——在10ns内根本没来得及交换。这些水分子的位置和取向直接影响PB（Poisson-Boltzmann）计算的介电边界和溶质-溶剂界面能。水分子没采样全，溶剂化自由能的统计误差不是一个常数，而是与溶质构象耦合的系统误差。

50ns：主要构象被采样了，但尾部不够

把模拟拉到50ns之后，RMSD没有明显变化（仍在1.8Å附近），但配体的构象空间明显展开了。6个可旋转键中5个完成了至少一次翻转，剩下那一个（与苯环相连的-C-O-键，翻转能垒约6 kcal/mol）在50ns内没有翻转。缺失这一个旋转自由度对熵的直接贡献只有约0.5 kcal/mol，但间接影响更大——这个键连接着配体的一个关键疏水基团，未翻转状态下该基团始终埋在结合口袋的疏水区；翻转之后这个基团会短暂暴露到溶剂里，引起附近3-4个水分子的重排。

结合能模拟值从10ns的-5.2提升到了-9.8 kcal/mol，离实验参考值（-12.0）还差约2.2 kcal/mol。差距来自三个来源：缺失的那个高能垒旋转自由度的熵贡献（~0.5 kcal/mol）、那些未交换的水分子的溶剂化自由能残余误差（~0.8 kcal/mol）、以及配体结合/解离路径上蛋白质构象变化的可逆功（~0.9 kcal/mol，MM/PBSA的单轨迹近似会忽略这部分）。

200ns：采样不足的最后一公里

200ns模拟里，所有6个可旋转键都至少翻转了2-3轮。结合口袋里的桥接水分子也有充分的交换——每个水分子在口袋里的平均驻留时间约2-5ns，200ns里有40-100轮的交换事件，统计上已经达到了收敛的标准（相对标准误差<5%）。

最终的MM/PBSA结合能是-11.4 kcal/mol，离实验值-12.0 kcal/mol的差距缩到了0.6 kcal/mol（约2.5 kJ/mol）。这个残余差距可以归因到力场精度本身和MM/PBSA方法的内在近似（比如PB连续介质模型的介电常数设定）。对于药物设计前期的筛选来说，0.6 kcal/mol的偏差完全可以接受——给定纳摩尔级Kd的实验测量误差通常在0.2-0.3 kcal/mol的级别，MD+MM/PBSA能做到0.6 kcal/mol已经是当前方法链的合理精度上限。

三组数据给出的经验法则

把这个经验推出去之前，需要申明一件事：200ns不是一个普适的”够了”的答案。这个体系只有320个残基和一个420Da的配体，计算量适中。如果换成500+残基的GPCR或1500Da的大环肽配体，200ns可能刚够蛋白质松弛，配体的构象平衡还远没开始。

但三组对比给出了一条判定准则而非固定数值：在MM/PBSA计算之前，检查配体的所有可旋转键是否至少完成了一轮翻转。如果有任何一个可旋转键从头到尾待在一个构象里没动，不要算MM/PBSA——算出来的ΔG偏差可以大到让一个纳摩尔级的高活性分子被误判为微摩尔级的弱活性分子。

决定精度的不是CPU小时数，是构象的遍历性。而遍历性，在RMSD收敛之后才刚刚开始。