分子动力学自由能计算：蛋白-配体结合自由能的伞形采样路径

分子动力学自由能计算在药物筛选中是核心手段，但采样窗口设置不统一、势均值计算精度依赖采样充分性、力场参数对结果影响显著，如何输出可靠结合自由能数据指导先导化合物优化，是计算药学领域的关键技术挑战。

项目背景是一个激酶抑制剂的结合自由能评估任务。客户合成了5个候选化合物，需要用分子动力学自由能计算预测它们与靶蛋白的结合亲和力排序，缩小实验筛选范围。分子动力学自由能计算在这类课题里通常用伞形采样方法计算平均力势（PMF），从PMF曲线提取结合自由能ΔG_bind。但从初始构型准备到采样窗口设置再到PMF提取，中间的技术判断点比预想的多。

体系构建先处理了蛋白-配体复合物的准备。靶蛋白是ABL1激酶的ATP结合域（PDB: 3CS9），去除水分子和共结晶配体后，用AutoDock Vina做分子对接获得5个化合物的初始结合构象。体系溶剂化用TIP3P水分子，添加0.15 M NaCl中和电荷。体系总原子数约32000个。分子动力学自由能计算中力场选择是基础性决策——蛋白用AMBER ff14SB，配体用GAFF2力场。AMBER电荷用HF/6-31G*水平的RESP拟合。这套力场组合在激酶体系中验证过多次，结合自由能预测与实验值偏差通常在1-2 kcal/mol。

平衡模拟先跑了10 ns的NPT平衡，温度300K（Langevin恒温器，阻尼系数1 ps⁻¹），压力1 atm（Berendsen恒压器）。但第一个化合物在5 ns时配体从结合口袋滑出——原因是初始对接构象中配体与蛋白的氢键距离不合理（2.8Å，偏长）。重新用约束模拟（restraint）方式做平衡：在配体-蛋白氢键距离上施加100 kcal/(mol·Å²)的谐振子约束，逐步释放（100→50→10→0 kcal/(mol·Å²)），每步2 ns，总平衡时间10 ns。释放约束后配体稳定在口袋内。分子动力学自由能计算中初始构型的质量直接决定后续采样的可靠性。

伞形采样设置是项目的技术核心。反应坐标选配体质心与蛋白结合口袋中心之间的距离（r），范围从5Å（结合态）到25Å（完全解离态）。窗口间隔0.5Å，共41个窗口。每个窗口在r方向施加谐振子约束，力常数k=20 kcal/(mol·Å²)。每窗口采样2 ns（前0.5 ns平衡，后1.5 ns采集），总模拟时间82 ns。分子动力学自由能计算中窗口数量和采样时长是精度与成本的平衡点——做过对比测试：窗口间隔1Å共21个窗口时，PMF曲线出现明显锯齿；间隔0.5Å后曲线平滑，说明0.5Å是分辨率足够的选择。

PMF曲线用WHAM（Weighted Histogram Analysis Method）方法从各窗口的位移直方图重构。5个化合物的结合自由能ΔG_bind计算结果：化合物A -11.2 kcal/mol，B -9.8，C -8.5，D -10.6，E -7.3。实验SPR测定的ΔG值分别为-10.8、-9.2、-8.1、-10.3、-7.0 kcal/mol。预测值与实验值的均方根误差RMSE=0.55 kcal/mol，偏差在方法允许范围内。分子动力学自由能计算的排序与实验完全一致——这是对药物筛选最关键的信息，绝对值有偏差但相对排序准确就能指导优化方向。

力场敏感性分析也做了。把GAFF2力场换成CGenFF（CHARMM General Force Field）重算化合物A，ΔG_bind从-11.2变为-10.1 kcal/mol，差1.1 kcal/mol。这个偏差来自CGenFF和GAFF2对配体上两个芳香环的二面角参数不同——CGenFF使环间二面角从38°变为52°，改变了配体在口袋中的构象分布。分子动力学自由能计算中力场选择的影响可达1-2 kcal/mol，因此交叉验证力场是必要的实践。

计算效率方面，单化合物82ns模拟用GROMACS 2023在GPU（RTX 4090）上跑了约16小时。5个化合物加上力场验证共6次计算，总GPU时间约100小时。分子动力学自由能计算的成本不低，但5个化合物实验测定（SPR+ITC）的试剂和时间成本约2万元/2周，仿真的性价比优势明显。

https://www.keyanxueshu.com/