分子动力学轨迹分析：RMSD 与 MM/PBSA
分子动力学（MD）模拟按力场对体系中每个原子的运动做数值积分，得到由成千上万帧组成的轨迹。轨迹分析就是从这些坐标中提取信息：RMSD 看整体是否进入平衡，RMSF 看哪些残基柔性大，氢键和接触分析看配体与哪些残基持续作用，MM/PBSA 或 MM/GBSA 对结合自由能做近似估算。输入通常是 GROMACS 的 tpr 与 xtc，或 AMBER 的 prmtop 与 nc 等拓扑加轨迹文件。药物研究中最常见的用途是检验对接姿势在溶液中是否稳定，并比较候选化合物或突变体的结合差异。
什么时候用
- 检验分子对接得到的复合物姿势在模拟中是否稳定
- 比较野生型与突变体的整体稳定性和局部柔性差异
- 找出配体与蛋白之间占有率高的氢键、盐桥和关键残基
- 对结合同一口袋的一组化合物做结合自由能排序
- 用主成分分析描述蛋白的主要运动模式和构象状态
不适合的情况
- 用 MM/PBSA 绝对值直接对应实验 ΔG：方法有系统偏差，适合相对比较
- 只有一条很短的轨迹却要下构象变化的结论：采样不足
- 配体电荷或质子化态设置有误的体系：分析无法弥补参数问题
需要准备的数据
- 拓扑文件
- GROMACS 的 tpr 或 top，AMBER 的 prmtop，CHARMM 的 psf；原子顺序必须与轨迹一致
- 轨迹文件
- xtc、trr、dcd 或 nc 格式；注明保存间隔和总时长，最好有多条独立重复轨迹
- 原子分组
- 蛋白、配体、口袋残基等索引组，用于叠合和计算；需明确配体的残基名
- 配体力场参数
- GAFF2、CGenFF 等生成的参数文件；计算结合自由能时必须与模拟所用参数一致
- 参考结构（可选）
- 起始结构或晶体结构，用作 RMSD 参考，并与模拟中的代表构象对比
分析步骤
- 1
处理周期性边界
模拟中分子会跨越盒子边界而显得断开。先用 gmx trjconv 的 -pbc mol 和 -center 把复合物居中成整，再做叠合；可视化确认蛋白和配体没有被拆散。
- 2
叠合去除整体运动
以蛋白主链或 Cα 为参照叠合每一帧，去掉整体平移和转动。算配体 RMSD 时同样按蛋白叠合，不再对配体自身叠合，这样才能反映配体在口袋中的位移。
- 3
判断平衡
画出蛋白主链 RMSD 随时间的曲线，找到进入平台的时间点，之前的部分作为平衡期，不纳入后续统计。多条重复轨迹的平台值应彼此接近。
- 4
柔性与紧密度
计算 Cα RMSF、回转半径和溶剂可及表面积。RMSF 峰通常在环区和末端，可与晶体 B-factor 对照；回转半径持续增大提示去折叠或结构域分离。
- 5
相互作用统计
统计配体与蛋白之间氢键、盐桥、疏水接触在所有帧中的占有率，并画出配体到关键残基的距离随时间变化。持续存在的作用比单帧截图更有说服力。
- 6
构象分析
对主链坐标做 PCA，把轨迹投影到前两个主成分上画自由能形貌图；对轨迹聚类得到代表构象，用于展示结合模式或作为后续对接的受体构象。
- 7
结合自由能计算
从平衡后的轨迹等间隔抽帧，用 gmx_MMPBSA 或 AmberTools 的 MMPBSA.py 计算 MM/PBSA 或 MM/GBSA，并做逐残基能量分解，找出贡献最大的残基。
关键参数与决策
分析时间窗
只用平衡后的轨迹做统计和自由能计算。起点按 RMSD 平台判断，不建议一律取后半段；不同体系的平衡时间可能差别很大。
叠合原子组
蛋白通常用主链或 Cα 叠合。末端或长环柔性很大时，只用二级结构核心区域叠合，避免柔性区抬高整体 RMSD、掩盖真实的结构变化。
抽帧间隔
MM/PBSA 不需要用全部帧，相邻帧高度相关。通常等间隔抽取数十到数百帧，增加帧数后结果不再明显变化，说明已经收敛。
溶剂模型与介电常数
PB 比 GB 更耗时，GB 有多个变体（如 igb=2、igb=5）。溶质内介电常数常用 1，带电或极性口袋有研究改用 2–4；选择会影响排序，同一研究内要保持一致。
熵项处理
正则模式或准简谐方法估算熵计算量大、误差也大，许多研究省略熵项只做相对比较；配体结构差异较大时忽略熵可能带来偏差，可考虑相互作用熵等近似。
独立重复次数
用不同初始速度跑多条重复轨迹，比单条更长的轨迹更能反映结果的稳定性。常见做法是至少 3 次重复，报告均值和标准差。
常见错误
- 未处理周期性边界就算 RMSD：配体或某条链跨越盒子边界，会让 RMSD 突然跳到几十埃，被误读为配体解离。
- 计算配体 RMSD 时对配体自身叠合：这只反映配体内部的构象变化，看不出它是否已经离开口袋，应以蛋白为参照叠合。
- 只用一条轨迹比较两个化合物：单次模拟的偶然性很大，曲线差异可能完全来自随机初始速度，需要多次重复。
- 把 MM/PBSA 得到的几十 kcal/mol 与实验 ΔG 直接对比：MM/PBSA 绝对值通常明显偏大，只适合同类体系内部排序。
- 不同化合物使用不同的电荷方法或力场参数：结合能差异会混入参数化差异，比较失去意义。
- 只展示 RMSD 曲线就宣称结合稳定：RMSD 平稳不等于关键相互作用保持，还应报告氢键占有率和配体到口袋的距离。
结果怎么看
RMSD 曲线
蛋白主链和配体（按蛋白叠合）的 RMSD 随时间变化。理想情况是早期上升后进入平台；配体 RMSD 持续上升或突然跳变时，结合可视化确认是否离开口袋；多条重复走势应一致。
RMSF 曲线
每个残基 Cα 在平衡期的波动幅度。对比野生型与突变体、结合前后，看口袋和突变位点附近柔性是否改变；整条曲线都很高时，先检查结构是否去折叠。
氢键与接触占有率表
每对相互作用在全部帧中出现的比例。占有率高的氢键可视为稳定作用；起始结构中存在、模拟开始后很快消失的作用，说明对接姿势在这里不稳定。
自由能形貌图
以前两个主成分或 RMSD 与回转半径为坐标的二维能量图。单一深谷说明构象集中；多个分开的低能区提示存在亚稳态，可分别提取代表结构比较。
MM/PBSA 能量分解表
把结合能拆成范德华、静电、极性溶剂化和非极性溶剂化各项，并给出逐残基贡献。看主导项是否合理，贡献最大的残基是否与相互作用统计和突变数据一致。
和相近方法怎么选
MM/PBSA 与 MM/GBSA
两者区别在极性溶剂化能的计算：PB 数值求解泊松-玻尔兹曼方程，GB 用广义玻恩近似，速度更快。大量化合物排序常用 GB，PB 理论上更细但对参数敏感；两者排序一致时结论更稳。
MM/PBSA 与自由能微扰（FEP）
FEP 和热力学积分通过炼金术路径计算相对结合自由能，精度通常更高，但需要大量采样和精细设置，适合先导优化中比较结构相近的类似物；MM/PBSA 成本低，更适合初步排序和姿势验证。
GROMACS 自带工具与 MDAnalysis、MDTraj
gmx 命令开箱即用，适合 RMSD、RMSF、氢键等标准分析，但定制能力有限；MDAnalysis 和 MDTraj 是 Python 库，能读多种轨迹格式，便于写自定义指标和批量处理。标准指标用 gmx，复杂分析用 Python 库。
RMSD 与 RMSF
RMSD 描述每一帧整体偏离参考结构的程度，随时间变化，用来判断平衡和整体稳定；RMSF 是每个原子或残基在时间上的平均波动，用来定位柔性区域。两者回答不同问题，通常一起报告。
常见问题
分子动力学 RMSD 多少算稳定？
没有统一阈值，关键是平衡后是否围绕某个值小幅波动，而不是绝对数值。常见的球状蛋白主链 RMSD 平台多在 1–3 Å，多结构域或柔性区多的蛋白会更高。配体 RMSD 重点看是否持续上升或出现跳变。
分子动力学模拟需要跑多长时间？
取决于要回答的问题，判断依据是关注的量是否收敛：RMSD 是否进入平台，氢键占有率在前后两半轨迹中是否一致，MM/PBSA 在增加帧数后是否稳定。研究大尺度构象变化需要更长时间或增强采样，多次独立重复比单纯延长一条轨迹更有说服力。
MM/PBSA 结合自由能是正值说明什么？
常见原因有：配体在模拟中离开了口袋，所用帧包含未平衡阶段，配体电荷或质子化态设置有误，或极性溶剂化惩罚过大。先看轨迹确认配体是否留在口袋，再检查参数和计算设置，不要直接解释为不结合。
GROMACS 轨迹分析常用哪些命令？
trjconv 处理周期性边界和叠合，rms、rmsf 算 RMSD 与 RMSF，gyrate 算回转半径，sasa 算溶剂可及面积，hbond 统计氢键，covar 加 anaeig 做 PCA，sham 画自由能形貌图，cluster 做构象聚类，均以 gmx 开头调用。
MM/PBSA 结果能和实验 IC50 对应吗？
不能直接换算。MM/PBSA 忽略或粗略处理熵，溶剂模型也有近似，绝对值通常偏大。它在同一靶点、同类化合物中的排序有时与实验趋势一致，因此适合用来排序和筛选，具体活性仍以实验测定为准。
在 Codar 里怎么做
在 Codar 里，你提供拓扑和轨迹文件，说明体系组成和比较对象。Codar 先确认原子分组和平衡时间窗，再做轨迹分析（周期性处理、RMSD、RMSF、回转半径、氢键占有率、PCA），并计算 MM/PBSA 或 MM/GBSA，验证对接姿势稳定性。交付曲线图、能量分解表、代表构象和可复现脚本。
我有 3 条 100 ns 的激酶-抑制剂复合物 GROMACS 轨迹，帮我算 RMSD、RMSF、氢键占有率和 MM/GBSA，判断对接姿势是否稳定。