分子对接和分子模拟：从AutoDock Vina到MD验证的完整工作流

分子对接和分子模拟，在药物设计工作流里是两个阶段——对接负责快速筛选结合模式，MD负责在溶剂环境下验证结合稳定性和计算结合自由能。问题在于，对接打分的Top 1构象，在后续MD里经常散架——蛋白侧链重排、溶剂分子挤进结合口袋、或者配体转了90度，让原本”完美”的氢键完全消失。

对接打分函数的陷阱：Affinity不等于结合自由能

AutoDock Vina的默认打分（affinity，单位kcal/mol）是一个经验性的打分函数，它把范德华、静电、氢键、去溶剂化等项用权重加起来，再拟合到实验结合常数上。这个打分在同源蛋白-配体对上表现不错，但对全新靶点或者别构口袋，假阳性率能到40%以上。

实际经验：对接做完之后，永远不要只看Vina打分排前10的构象。要用聚类分析（比如按RMSD聚类）把构象分成3-5个簇，每个簇选一个代表去跑MD——这样能覆盖到打分稍低但结合模式完全不同的可能性。

蛋白结构准备：别让PDB的坑毁掉整个对接

从PDB下载的晶体结构，直接拿来做对接的人我见过太多了。问题有几个：

1. 缺失侧链/_loop区：PDB里很多低分辨率结构的柔性区域是直接缺掉的。用Modeller或者Swiss-PDBViewer补完之后再对接，否则结合口袋的形状会被人为”封死”。
2. 水分子的处理：有的晶体结构里关键水分子（比如人Hsp90的ATP结合位点里的水桥）是结合必需的，直接删掉所有水会导致对接构象完全偏离。
3. 质子化状态：生理pH下，蛋白的His可以是HID/HIE/HIP三种质子化状态，配体的羧基可以是COO⁻或者COOH。用pdb2gmx（GROMACS）或者Reduce（PHENIX套件）来处理质子化，不要直接用PDB原始文件的质子化状态。

对接盒子（Grid Box）的设置原则

对接盒子的中心一般设在已知活性位点的几何中心，盒子尺寸要能覆盖配体的全部构象自由度。经验值：盒子边长 = 配体最长尺寸 + 8-10 Å。

但有一种情况需要特别处理：别构位点（allosteric site）。如果已经有实验证据或者同源建模表明别构口袋存在，盒子要另外设在这个位置，不能只盯着正构位点。更激进一点，可以用盲对接（blind docking）——把盒子设成覆盖整个蛋白表面（边长40-60 Å），让Vina自己去”找”结合位点。代价是计算时间暴增，而且对计算资源的要求高得多。

MD验证：从对接构象到RMSD稳定

对接得到的Top构象，导入GROMACS做溶剂化MD，标准流程是：

1. 能量最小化（steepest descent，500步起）→ 消除对接构象里原子的不良接触。
2. NVT预热（100 ps，从0 K升温到300 K）→ 让溶剂水逐步包围蛋白-配体复合物。
3. NPT平衡（1-2 ns）→ 盒子尺寸收敛，压强稳定在1 bar。
4. 生产运行（10-50 ns）→ 分析RMSD、RMSF、氢键占有率。

判断对接构象是否”真的稳定”，核心看配体重原子的RMSD是否已经收敛（通常<2 Å波动），以及初始对接时的关键氢键在MD里是否还能维持（占有率>30%才勉强算稳定）。

MM-PBSA/GBSA：从MD轨迹算结合自由能

MD跑完之后，可以用**MM-PBSA（Poisson-Boltzmann Surface Area）或者MM-GBSA（Generalized Born Surface Area）**从轨迹帧里算结合自由能ΔG_bind。

GROMACS用户可以用g_mmpbsa工具（第三方，需编译）或者pmx套件。核心公式拆解：

ΔGbind=ΔEMM+ΔGsolv−TΔSconfΔGbind​=ΔEMM​+ΔGsolv​−TΔSconf​

其中ΔE_MM可以在每帧轨迹里直接算（蛋白-配体相互作用能），ΔG_solv需要解Poisson-Boltzmann方程（PBSA）或者用广义玻恩模型（GBSA）近似，熵项TΔS最麻烦，一般用正态模分析（Normal Mode Analysis）近似算，或者用溶剂化自由的的温度导数来估算。

参考文献与数据来源

• AutoDock Vina官方文档 —— 打分函数原理、参数设置和输出解读
• GROMACS官方MM-PBSA教程 —— g_mmpbsa的安装和使用示例
• J. Med. Chem. 59, 6479 (2016) —— 分子对接与MD结合自由能计算的系统对比，假阳性分析

分子对接和分子模拟不是”先对接后MD验证”这种流水线关系。对接的打分函数局限性、蛋白结构的质量、MD验证的RMSD判据——这三关把住了，计算药物设计的结果才能和实验对上。