分子对接的结合能：评分函数精度、构象采样与MM/PBSA后处理的完整评估

分子对接的结合能是整个计算机辅助药物设计(CADD)流程中最容易被误用的数值。对接软件输出一个结合能数值——AutoDock Vina给出一个”binding affinity”(通常以kcal/mol为单位)，很多非计算背景的用户把它当作实验结合自由能的直接估计，拿来做候选化合物的排序和取舍。这个操作如果不对接出一个理解——数值与实际的差距——比不用对接更危险：它会放大一个系统性偏差，让排序在几百个候选物中指向错误的化合物。

评分函数的先天局限：精度天花板在哪

市面上主流的对接评分函数可以分为三类：

**基于力场的评分函数**(AutoDock 4. DOCK 6)：把结合能拆分为范德华作用(Lennard-Jones项)、静电作用(库仑项——通常用距离依赖介电常数或GB模型屏蔽)和去溶剂化能(SASA比例项)。这种拆分的物理基础坚实——分子间相互作用力的确是这些项的加和——但精度受限于参数的粗粒化：同一组LJ参数要描述蛋白主链酰胺氮和配体甲基氢之间的作用，也要描述蛋白侧链羧基氧和配体芳香环之间的作用。力场的”普适性”是以精度为代价的。

**基于经验的评分函数**(AutoDock Vina, ChemScore, X-Score)：用大量已知蛋白-配体复合物结构和对应的实验结合亲和力数据(如PDBbind数据库)，训练一个多元线性回归模型。X-Score的训练集包含约200个蛋白-配体复合物，通过拟合范德华项、氢键项、疏水项、熵项的加权组合来重现实验结合能。优势是不含不明物理意义的项——每一项对应一个可解释的物理作用;劣势是训练集之外的配体化学空间——特别是含有非典型官能团(硼酸、磺酰胺等)的配体——预测误差可能从±2 kcal/mol跳升到±5 kcal/mol以上。

**基于知识的评分函数**(PMF, DrugScore)：从蛋白-配体复合物晶体结构数据库中统计原子对距离的分布概率，用反转Boltzmann关系推导出”平均力势”(PMF, Potential of Mean Force)。DrugScore的构建使用了约600个蛋白-配体复合物的结构数据。优势是训练过程不依赖实验结合能数据，不受实验测量误差的影响;劣势是统计精度高度依赖数据库对特定原子对类型的覆盖密度——F·O(氟-氧)这类组合在PDB中样本量偏小，相应的PMF分布不够收敛(noise-to-signal ratio偏高)。

这三种评分函数的共同天花板：它们在单点结构上计算结合能，没有考虑蛋白和配体在结合过程中的构象变化和熵损失——而熵损失通常贡献结合自由能的20~50%。一个只看焓、不看熵的评分函数，就像只记收入不算支出——结果的方向可能对，但数值绝不精确。

对接构象 ≠ 结合构象：采样充分性检验

分子对接的结合能依赖于一个前提：对接算法找到的低能构象接近真实的结合构象(即晶体结构中的pose)。这个前提在”自对接”(self-docking，把晶体结构中的配体拿出来再对回同一个蛋白)中成立的概率约70-80%——这个数字来自多个对接软件的benchmark研究(如CSAR 2014 benchmark)。但在”交叉对接”(cross-docking，把配体A对到配体B的蛋白构象上)中，成功率骤降到30-50%。

原因不复杂：蛋白在结合不同配体时，结合口袋的构象会有适应性变化——侧链旋转、loop位移、甚至结构域的微小重排。如果对接时把蛋白当作刚体(大多数对接软件的默认设置)，配体只能在一个”不完美”的口袋里寻找最佳位置，找到的构象当然不是全局最优——它是”在这个错误口袋形状下的局部最优”。

对于对接结果的构象可靠性评估，有几个实操标准：

1. 配体的内部几何是否合理——键长、键角、二面角不能偏离力场标准值太多(>3σ) 2. 配体与蛋白的互补性——疏水基团对疏水口袋，极性基团对极性口袋或暴露溶剂 3. 氢键的几何参数——供体-受体距离 2.5~3.5Å，角度偏差<30°(偏离太多，这个氢键在物理上是弱的) 4. 如果存在同一靶点的多个已知配体的晶体结构，比较对接pose与已知配体结合模式的一致性——结合模式(binding mode)是否重现，是比结合能数值更可靠的验证指标

一个项目里的实际做法：对接时跑三套参数——默认、exhaustiveness加倍、增加额外扭转自由度扫描——三套参数产出的top-1构象如果一致，这个构象的可信度更高;如果不一致，说明构象采样不充分，需要加大采样力度(增加exhaustiveness值，或者换用能处理蛋白柔性的对接软件如AutoDockFR、RosettaLigand)。

结合能重打分：不要用对接评分函数做最终排序

从对接软件直接输出的结合能数值来排序候选化合物，是一个普遍性错误。对接评分函数的设计目标是”快速区分结合物和非结合物”——它在筛选早期阶段的虚筛中有效(以AUC-ROC约0.7~0.8的水平)，但在”精细排序”这个精度层级上不成立。

分子对接的结合能评估需要一个两步走：对接阶段用评分函数做粗筛(从几千~几万候选物筛到50~100个)，精细排序阶段用更高精度的评分方法重打分。

**重打分的方法选择：**

**基于机器学习的评分函数(RF-Score、NNScore)**：用随机森林或神经网络在PDBbind的约5000~15000个蛋白-配体复合物上训练，输入特征是蛋白-配体原子对计数的向量化表示。RF-Score v3在CSAR benchmark上Pearson相关系数约0.8(vs实验结合能)，显著优于Vina评分函数(约0.55~0.6)。代价：模型的可解释性低——你能看到一个预测值和置信区间，但不知道为什么这个配体得分高——这在高通量筛选中可以接受，在需要结构-活性关系(SAR)优化的场景下不够有用。

**MM/GBSA和MM/PBSA**：用分子动力学模拟生成的构象系综来计算结合自由能，基于分子力学能量(MM)加上隐式溶剂计算的溶剂化自由能(GBSA或PBSA)。平均每套MM/PBSA对一个配体的计算成本约几十至几百CPU核时(取决于MD时长和帧数)，比对接评分函数重了约100~1000倍的计算量——但这恰好是它存在的理由：用计算量换取精度。

MM/GBSA在几个benchmark上的Pearson相关系数约0.6~0.7.MM/PBSA略高(约0.65~0.75)，两者都优于对接评分函数，但依然不是实验精度(实验结合能测量的重复性误差在±0.5 kcal/mol)。关键瓶颈在熵贡献的计算——大多数MM/PBSA计算使用简正模分析(normal mode analysis)来估计振动熵变化，而振动熵在蛋白-配体系统中的贡献被普遍低估。实验上测得的结合熵贡献(-TΔS)可能在-5~+5 kcal/mol的区间内波动，简正模分析给出的通常只有-2~0 kcal/mol。

从对接到MD：三级递进的结合能评估框架

分子对接的结合能评估在项目中逐步积累了一套三级递进框架：

**第一级：高通量虚筛**——用对接评分函数(AutoDock Vina, Glide SP)做初筛，从数万到数十万个候选物中筛选出前200~500个。这一级只看”能不能结合”，不关心排序——AUC 0.7~0.8意味着前500个里大约有350~400个是”可能结合”的，其中混杂了100~150个假阳性。这个代价在算力上可以接受——因为下一步会把这些假阳性过滤掉。

**第二级：一致性对接+重打分**——对第一级筛出的200~500个候选物，用至少两套不同的对接软件(如Vina + Glide或Vina + DOCK)做对接，只保留两套软件都给出类似结合构象(RMSD<2Å)的候选物。这一轮过后数量大约压缩到50~100个。然后对这50~100个用ML评分函数(如RF-Score)+ MM/GBSA做重打分，取两个评分的加权排名前列(各占50%权重，相关性不足的评分降权)，选出前10~20个。

**第三级：MD验证**——对前10~20个候选配体-蛋白复合物分别跑50~200ns的显式溶剂MD模拟，用最后50ns的轨迹做MM/PBSA结合自由能计算(取100~200帧，每帧间隔~500ps，保证帧之间的统计独立性)。在这一级，如果前三级给出的排序和前两级差异很大——例如对接评分最高、重打分第二高的配体在MD精算后掉到第15名——一以MD精算为准。MD产生的构象系综包含了蛋白口袋在结合后的弛豫效应，这是对接和单点评分函数都无法捕捉的物理过程。

结合能数值的解读底线

分子对接的结合能在最终报告中不应该只是一个数值。它需要一个完整的上下文：

这个结合能来自哪个计算方法(对接评分/MM-GBSA/MM-PBSA)、使用了哪个力场参数集、MD采样时长、以及与已知活性化合物(阳性对照)的对比。没有阳性对照的结合能排序，缺乏验证基准——你无从判断排序的前几名究竟是”真有活性”还是”计算误差的受益者”。

此外，对接结合能的排序在1~2 kcal/mol范围内的精细差异(如-8.5 vs -8.3 kcal/mol)不具统计意义——评分函数和MM/PBSA的误差范围在±2~3 kcal/mol量级。只有当差距>3 kcal/mol(约100倍结合常数差异)时，排序才是可靠的。

这也许是对分子对接的结合能最诚实的交代：它能告诉你谁是排名前10%的候选物、谁在后50%，但10%内部的精细排序，需要实验去完成。把对接当作一个筛子而不是一台天平——这个定位明确了，用起来才不会犯错。