手机版

蛋白配体结合能计算 — MM/PBSA与FEP方法的选择与应用

发布时间：2026-07-16 来源：科研学术网

小中大

字号：

蛋白配体结合能计算是计算机辅助药物设计（CADD）里最核心也最让人头疼的环节。对接能给你一个大概的排序，但要做精细的lead optimization，你需要的是和实验IC50能有定量对标关系的结合自由能。目前最主流的两条路线——MM/PBSA（分子力学/泊松-玻尔兹曼表面积）和FEP（自由能微扰）——各有优劣，选哪个、怎么跑、结果怎么解读，每个环节都有让人翻车的陷阱。这篇文章是我做了三年多结合能计算项目后复盘出来的经验，重点讲MM/PBSA和FEP的实战决策。

一、蛋白配体结合能计算的两条路线：MM/PBSA还是FEP

蛋白配体结合能计算选择MM/PBSA还是FEP，取决于你的具体需求、资源约束和精度要求。没有绝对的优劣，只有适合不适合。

先说MM/PBSA。它的基本原理是从MD轨迹中抽取快照，对每个快照分别计算气相分子力学能量（ΔE_MM）、极性溶剂化自由能（ΔG_PB，用泊松-玻尔兹曼方程求解）和非极性溶剂化自由能（ΔG_SA，通常用溶剂可及表面积SASA估算），然后做统计平均。MM/PBSA的优势是计算效率高——一个配体体系跑完MD后，MM/PBSA分析通常在几小时内就能完成（取决于轨迹帧数）。而且MM/PBSA不依赖配体之间的结构相似性，可以对任意化学结构的配体打分。

但MM/PBSA有明显的精度上限。因为它不显式考虑结合过程中的构象熵变（除非做昂贵的normal mode分析），对配体结合口袋里的水分子处理也比较粗糙。在基准测试中，MM/PBSA的预测结合能和实验值的相关系数（R²）通常在0.3-0.6之间，RMSE在2-4 kcal/mol。这意味着MM/PBSA适合做”粗筛”——从几十个候选里筛出最强的几个——但不适合做精细的”这个甲基该放哪里”级别的优化。

FEP就不一样了。FEP通过一系列λ窗口在两个状态（比如配体A和配体B）之间做可逆的”变换”，直接计算自由能差。因为FEP考虑构象熵变、水分子重排等细节，精度远高于MM/PBSA。在理想条件下（充分采样、合理的λ窗口分布、力场参数准确），FEP的预测精度可以达到1 kcal/mol以内。文献里FEP和实验IC50的RMSE通常在0.8-1.5 kcal/mol，R²在0.6-0.8。

FEP的代价是计算量巨大。一个配体对（A→B）的FEP计算通常需要12-24个λ窗口，每个窗口跑5-20 ns的MD。对于蛋白-配体体系（通常2-5万原子），这意味着一对配体需要几万到十几万核时。而且FEP要求两个配体之间有明确的化学变换路径（perturbation map），不适合比较结构差异太大的配体。

实战中的策略通常是：用MM/PBSA做初筛（从100个候选筛到20个），然后用FEP对Top 20做精细排序。如果计算资源有限，可以用MM/PBSA+配体效率指标（LE, LLE等）做多维度评分，替代FEP。

二、分子对接后的MD平衡：不是跑够时间就行

蛋白配体结合能计算的第一步通常是用分子对接获得初始的蛋白-配体复合物结构，然后跑MD做平衡。很多人觉得”跑个10 ns MD平衡一下就行了”，这是最大的误区。平衡不只是跑够时间，而是跑对条件。

第一个关键问题是蛋白的准备。对接用的蛋白结构通常来自PDB晶体结构，里面可能有结晶水、辅因子、翻译后修饰、甚至配体本身。哪些水分子该保留？我的一般规则是：桥接水（同时和蛋白及配体形成氢键的水）必须保留，它们在结合中可能起关键作用。埋在蛋白内部的结构水也应该保留，因为它们维持蛋白的局部构象。表面水可以删掉，因为它们会在MD过程中快速交换。

辅因子和金属离子的处理更棘手。如果蛋白活性位点有Zn²⁺或Mg²⁺，你不能简单地用一个+2的点电荷表示——这样算出来的静电相互作用会严重高估金属离子的效应。对Zn²⁺，至少要用”锌指参数”（增加范德华半径来模拟配位饱和效应）或者用12-6-4 Lennard-Jones势（在LJ势中加一个C4/r⁴项来模拟电荷转移和诱导效应）。

第二个关键问题是平衡的判断标准。很多人只看RMSD——骨架RMSD平了就觉得平衡了。但RMSD平了不代表体系平衡了。配体在结合口袋里的构象可能在RMSD看起来稳定的时候还在缓慢调整——比如一个柔性侧链从一种rotamer翻转到另一种，RMSD可能只变0.3 Å，但结合自由能可能差1-2 kcal/mol。所以判断平衡要多维度看：骨架RMSD、配体RMSD、配体重原子RMSD、蛋白-配体氢键数目、配体的回转半径、以及结合口袋里的关键距离（比如催化残基和配体的距离）。

第三个关键问题：平衡阶段要不要加约束？很多人习惯在平衡的前几个ns给蛋白重原子加位置约束（position restraint），让溶剂和配体先适应。这个做法本身没问题，但约束解除的时机很关键。如果在体系还没充分溶剂化的时候就解除约束，蛋白可能发生大的构象漂移。我一般分三个阶段加约束：前1 ns强约束（力常数10 kcal/mol/Å²），中间1-2 ns弱约束（1-2 kcal/mol/Å²），最后完全放开。每个阶段都要确认温度和压力已经稳定。

三、MM/PBSA计算中极性溶剂化能的PB求解器选择

MM/PBSA计算中最关键的单项是极性溶剂化自由能ΔG_PB，它由泊松-玻尔兹曼方程求解得到。PB求解器的选择和相关参数设置，对最终结合能的影响可能大到几个kcal/mol。

Amber的MM/PBSA模块支持两种PB求解器：内置的pbsa和调用外部的Delphi。GROMACS用户通常用g_mmpbsa工具配合APBS求解器。不同求解器的核心算法差异不大（都是有限差分法求解PB方程），但网格设置、介电常数处理、离子强度处理等细节会导致结果差异。

最关键的两个参数是内部介电常数（ε_in）和外部介电常数（ε_out）。ε_in是溶质（蛋白和配体）内部的介电常数，ε_out是溶剂的介电常数。ε_out通常设为80（水的介电常数），争议不大。但ε_in的选择是一个持续了二十年的争论。

经典文献推荐ε_in=1——因为溶质内部是真空，没有电子极化。但ε_in=1算出来的极性溶剂化能通常过大，导致结合能远高于实验值。所以很多课题组用ε_in=2或4来”经验性”修正。实际上这个”经验修正”有物理依据：蛋白内部的局部环境不是真空，有肽键偶极、侧链极化等效应。ε_in=2-4相当于用一个等效均匀介质近似这些效应。

我做过的基准测试表明：对于带电配体（比如含羧酸或胺基的），ε_in=4给出的结合能和实验值对标最好。对于中性配体，ε_in=2就够了。但最严谨的做法是用ε_in做敏感性分析——分别用1、2、4算，看结合能的排序是否一致。如果排序一致，说明你的结论对ε_in不敏感；如果排序变了，需要更仔细地分析。

另一个重要参数是PB网格间距。默认值通常是0.5 Å，这对大多数体系够了。但如果配体很小（分子量<200）或者结合口袋非常狭窄，建议降到0.25 Å。网格太粗会在溶质-溶剂边界产生数值噪声，影响ΔG_PB的精度。代价是内存消耗增加（和网格间距的立方成反比）。

离子强度的影响相对较小。生理条件下的离子强度是0.15 M，PB方程里的Debye-Hückel屏蔽项在这个范围内对结合能的影响通常小于0.5 kcal/mol。但如果你的配体带多个电荷，离子屏蔽效应会更显著，建议做离子强度的敏感性分析（0到0.2 M）。

四、FEP计算中λ窗口分布和采样充分性验证

FEP计算的成败，60%取决于λ窗口分布和采样充分性。很多人用FEP的默认设置跑完就算，结果误差大到让FEP的优势荡然无存。

λ窗口分布决定了你把自由能变化ΔG”切成”多少段来计算。窗口太少，每个窗口的自由能变化太大，相邻窗口之间的相空间重叠不够，导致统计误差大。窗口太多，计算量线性增加。经验法则是：让相邻窗口之间的相空间重叠度（通过能量差分布的overlap来衡量）至少达到10-20%。

具体到λ窗口数目，取决于perturbation的类型。对于一个甲基→乙基的变化（增加一个CH₂），12个λ窗口通常够了。对于电荷变化（比如引入或消除一个电荷），需要更密的λ窗口（16-20个），因为静电相互作用的自由能变化更陡峭。对于环的变化（比如苯环→吡啶环），建议用18-24个窗口，因为需要同时处理范德华和静电两个维度的变化。

“软核势”（soft-core potential）的使用也很重要。当一个原子在λ=0或λ=1处”出现”或”消失”时，如果直接用标准的Lennard-Jones势，会产生奇异性（能量趋向无穷）。软核势用一个修正的LJ势函数，在原子出现/消失时避免奇异性。大多数FEP实现默认启用软核势，但软核参数（α和β）的选择会影响结果。默认值（α=0.5, β=12）对大多数体系工作良好，但如果你的perturbation涉及大体积变化（比如氢→叔丁基），可能需要调整。

采样充分性验证是FEP里最被低估的环节。常用的验证方法有几个：看正向（forward）和反向（backward）FEP的自由能是否闭合（hysteresis）；检查每个λ窗口的∂U/∂λ分布是否接近正态（如果不是，说明采样不够）；做不同初始结构的独立重复FEP计算看结果是否一致。

我遇到过的最离谱案例：一个FEP计算给出了-3.5 kcal/mol的ΔΔG，和实验值（-1.2 kcal/mol）差了2.3 kcal/mol。查了三天才发现某个λ窗口里蛋白的一个loop在两个亚稳态之间来回跳，导致该窗口的采样没有充分覆盖构象空间。延长该窗口的模拟时间到原来的3倍后，ΔΔG修正为-1.5 kcal/mol。

五、结合能结果与实验IC50的对标偏差分析

蛋白配体结合能计算做完后，最终要和实验IC50或Kd对标。对标偏差的来源很多，了解每一个来源能帮你判断你的计算结果到底可不可信。

第一个偏差来源：计算条件和实验条件的差异。计算是在”理想”条件下做的——固定pH、固定离子强度、不考虑蛋白聚集或降解。实验IC50的测定条件（pH、缓冲液、温度、孵育时间）和计算条件不完全一致，这个差异本身就能贡献0.5-1 kcal/mol的偏差。

第二个偏差来源：力场精度。Amber ff14SB/ff19SB蛋白力场和GAFF/GAFF2配体力场的结合能预测精度，在最好的情况下RMSE也在1-2 kcal/mol。如果你用的配体力场参数是GAFF自动生成的（比如用antechamber），有些原子类型的电荷和范德华参数可能不够准。对关键配体，建议用量子化学计算（HF/6-31G*级别的RESP电荷拟合）替代默认的AM1-BCC电荷。

第三个偏差来源：IC50到结合自由能的转换。关系是ΔG_bind = RT ln(IC50) + 常数，但这个转换假设了简单的竞争性结合机制。如果配体是非竞争性或变构抑制剂，IC50和结合自由能之间的关系就不是这么简单了。所以在对标之前，先确认你的配体是竞争性抑制剂。

第四个偏差来源：实验IC50本身的误差。同一个化合物在不同实验室测的IC50可以差2-5倍（对应0.4-1 kcal/mol的自由能差异）。所以当你看到计算和实验差1 kcal/mol的时候，先不要慌——这可能还在实验误差范围内。

第五个偏差来源：蛋白柔性。如果你的蛋白在结合配体后有大的构象变化（比如诱导契合），而你的MD没有采样到这个变化，算出来的结合能会偏高（不够负）。这是MM/PBSA和FEP共同的盲区——它们都假设起始结构能覆盖结合态的主要构象。

实用对标策略：不追求和单个配体的实验IC50精确对标，而是看一组配体（至少10-20个）的相对排序。如果计算能正确区分强结合（nM级）和弱结合（μM级）的配体，即使绝对值有系统偏差，对药物设计来说已经是可用的结果。排序正确率（rank correlation或enrichment factor）是比RMSE更有实际意义的评估指标。