分子动力学跑完后,20+项常见指标到底说明什么?

阅读提示：本文篇幅较长，系统梳理了20多项常见MD分析指标，建议先收藏，再根据自己的研究问题选择性阅读——这些指标并不是要求全部计算一遍。

从轨迹可信度、配体保留、口袋机制到自由能证据链

100ns跑完，并不意味着结论自动出现。RMSD进入平台，可能只是体系暂时停在一个亚稳态；配体RMSD升高，可能是离开口袋，也可能只是柔性尾部翻转；氢键数量差不多，背后的供体-受体却可能已经全部替换；FEL出现一个深色能量盆，也可能只是轨迹没有走出初始状态。

很多报告的问题不是图少，而是把不同层级的问题混在一起：用全蛋白RMSD替配体位置下结论，用氢键数量替结合自由能下结论，再用一个MM/PB(GB)SA负值把前面的故事“封口”。这样做看起来指标齐全，实际上证据之间并没有闭环。

对于蛋白、蛋白-配体和蛋白-蛋白等生物类体系，RMSD、RMSF、Rg、SASA和氢键只是常见结果中的一部分；根据研究对象和科学问题，还可以继续统计距离与角度、残基接触频率、主链与侧链二面角、Ramachandran分布、二级结构、口袋开合、界面埋藏、水网络、离子分布、动态网络以及多种能量和自由能量。本文选取20多类具有代表性的分析逐项讲解，但MD后处理远不止这些，不可能在一篇文章中一一列完。没有列出不等于没有价值，列出来也不等于每个项目都必须计算。

下面仍按“它到底算什么、真正能回答什么、不能说明什么、必须与什么联用”逐项拆开。所有论文图只作为图形阅读示例，不代表本文接受原论文对该图的每一句解释。

读图之前，先把四个概念分开

1.运行正常：温度、压力、密度和能量没有明显异常，说明积分与系综控制基本正常。

2.达到平衡：所关注的观察量不再持续受初始条件驱动，生产期统计可以开始。

3.结构稳定：某个特定尺度的结构或几何在观察时间内保持在一定范围；必须说明是全蛋白、结构域、口袋还是配体。

4.采样收敛：不同时间块或独立重复对目标分布、状态人口或平均量给出一致结果。单条曲线走平通常不足以证明。

最容易忽视的统计问题：MD相邻帧高度相关。每10ps保存一帧，得到1万帧，并不代表有1万个独立样本。时间分块、有效样本数、重复轨迹和置信区间，往往比再增加一张漂亮的曲线更重要。

先有科学问题，再选分析指标

分子动力学分析不应预设一个固定的基础图组合，也不存在从基础指标一路做满到高级分析的统一等级。运行质控、结构描述、作用网络、构象状态和自由能方法回答的是不同问题，不能被压成一个固定模板。更合理的顺序是：先确定要证明什么，再选择最贴近问题的直接指标；必要时增加另一类独立证据，最后才考虑能量或更高成本的自由能方法。

本文主要面向蛋白、蛋白-配体和蛋白界面的常规全原子MD。膜体系还需膜厚、脂链序参数、蛋白倾角和脂质接触；核酸体系还需碱基配对、堆积和螺旋参数；聚合物、自组装与材料体系则常以RDF、配位数、密度剖面、MSD/扩散和团簇尺寸为主。即使只讨论生物体系，也还有构象熵、氢键动力学相关函数、动态网络、Markov状态模型等专项分析没有在本文逐一展开。因此，本文的20多项指标不能被机械套用到所有体系。

一、先判断轨迹能不能用于回答问题

第1-9项关注运行质量、整体结构和轨迹内部状态。这里最重要的不是找到一条“变平”的曲线，而是确定统计窗口、局部与整体的区别，以及结果是否在重复轨迹中可复现。

1. 温度、压力、密度与能量

基础质控｜常被一笔带过，但应先检查

图1 温度、压力与能量的基础质控时间序列。取Al-Anazi等补充图Figure S1a、S2a与S3a的完整坐标面板作版式拼合，曲线数据未重绘。

它到底算的是什么：这些量检查的是积分、温控、压控和体积调节是否按预期工作。温度反映动能分布，压力是瞬时维里量，密度反映NPT体系体积是否调整到合理区间，势能则反映体系在当前哈密顿量下是否仍有持续的松弛或漂移。它们是运行质量指标，不是蛋白功能指标。

它真正回答的问题：它们首先回答轨迹有没有明显的数值异常、平衡阶段是否结束，以及哪一段可以进入生产期统计。温度和密度在目标附近波动、势能不再持续单向漂移，通常说明基本热力学条件已建立；压力在原子尺度体系中瞬时波动很大并不罕见，应看长时间平均和分布。

读图时先看：先把能量最小化、NVT/NPT平衡和生产模拟分开；标出舍弃的平衡区间；报告均值、标准差与趋势。NVE中应重点检查总能漂移，NVT/NPT中则不能要求总能严格守恒，因为热浴或压浴会与体系交换能量。周期性水相体系还应核对密度是否接近所用水模型和温压条件下的合理范围。

不能这样下结论：温度和压力正常，只能说明模拟器在正常工作，不能证明蛋白已经构象收敛，更不能证明配体结合稳定。反过来，压力曲线不平滑也不等于模拟失败。

最合适的联用方式：与蛋白RMSD、关键距离和分块统计联用。若体系出现能量突跳，还应回查约束、时间步长、邻居表、周期性边界和原子重叠，而不是继续解释后续生物学指标。

更稳妥的结果表述：生产阶段温度与密度围绕目标值稳定波动，势能未见持续漂移，表明轨迹具备进一步结构分析的基本条件；这一步不单独用于判断构象或结合收敛。

2. 蛋白骨架RMSD

核心指标｜看整体偏离，不是看结合强弱

图2 蛋白骨架RMSD时间序列示例，原论文Figure 6，来源见参考文献[1]。

它到底算的是什么：RMSD是在选定原子完成最小二乘叠合后，当前帧相对参考结构的均方根位移。结果高度依赖三个选择：用什么原子对齐、用什么原子计算，以及参考是初始结构、平均结构还是某个实验状态。所谓“蛋白RMSD”最好明确写成Cα RMSD、主链RMSD或结构域RMSD。

它真正回答的问题：它回答蛋白整体相对参考构象偏离了多少，以及偏离是否进入相对稳定区间。持续爬升提示体系仍在松弛或发生较大转换；台阶式变化提示状态切换；平台只表示在这个原子选择和参考系下，整体偏离暂时不再系统增加。

读图时先看：先看是否去除了整体平移和转动，再看平台出现的时间、平台宽度、是否发生多次跃迁。多结构域蛋白应同时计算全局RMSD和各结构域RMSD：两个结构域本身都稳定，但彼此转动时，全蛋白RMSD仍可能升高。若比较apo、配体、突变体，原子映射、参考结构和统计窗口必须一致。

不能这样下结论：RMSD低不等于结构更稳定，RMSD高也不等于解折叠；它更不等于亲和力。一个轨迹在50 ns后走平，仍可能只是困在一个亚稳态。

最合适的联用方式：与RMSD矩阵、聚类和重复模拟判断状态是否可重复；与Rg、SASA、DSSP区分刚体运动和真实展开；与配体RMSD及口袋距离区分蛋白整体和结合位点行为。

更稳妥的结果表述：主链RMSD在约X ns后进入稳定波动区间，说明该时间窗内整体构象相对参考结构不再持续漂移；是否达到充分采样仍需由重复轨迹和构象状态分析确认。

3. 配体RMSD：位置变化与自身变形必须分开

核心指标｜最容易因对齐方式而被误判

图3 蛋白与配体RMSD的同步展示。图例同时区分“(Lig) fit on Prot”和“(Lig) fit on Lig”，原论文Figure 6A，来源见参考文献[2]。

它到底算的是什么：至少存在两种不同的配体RMSD。其一是先用蛋白或口袋残基对齐，再计算配体重原子RMSD，它保留了配体相对口袋的平移、转动和内部变形；其二是再对配体自身拟合后的内部RMSD，只反映配体构象变化。柔性配体还可分解为核心骨架RMSD和可旋转基团二面角。

它真正回答的问题：蛋白对齐后的配体RMSD用于判断初始结合姿势是否保留、是否翻转、换位或离开口袋；配体自身拟合后的RMSD用于判断它有没有改变内部构象。两者一高一低时，常能区分“整体换了位置”和“原地发生变形”。

读图时先看：检查对齐对象、配体原子映射和对称性校正。具有等价原子的芳香环或对称基团可能在化学等价翻转后产生虚高RMSD。台阶式升高应回看对应帧：如果质心距离不变、接触网络重排，可能是口袋内换姿势；若距离、SASA同时增大且接触消失，才更支持脱离。

不能这样下结论：不能对配体自身最佳拟合后，再用这个低RMSD证明它一直留在口袋；也不能把柔性尾部摆动造成的高RMSD直接写成解离。

最合适的联用方式：与配体-口袋质心距离、关键锚点距离、接触占有率、配体SASA和代表构象叠合联用。判断解离时至少需要位置、接触和溶剂暴露三类证据同向。

更稳妥的结果表述：蛋白对齐后的配体RMSD出现一次姿势转换，但质心距离和核心接触保持稳定，提示配体在口袋内重排，而非完全解离。

4. RMSF：残基柔性

核心指标｜高峰需要映射回结构

图4 按残基编号绘制的RMSF示例，原论文Figure 9，来源见参考文献[3]。

它到底算的是什么：RMSF是每个原子或残基围绕其时间平均位置的均方根涨落。最常见的是Cα RMSF。它不是简单地把RMSD按残基拆开：RMSF需要先选择对齐区域，再在选定时间窗内定义平均结构。若用全蛋白对齐，多结构域相对运动可能被混入局部柔性。

它真正回答的问题：它回答哪些残基、loop、末端、口袋边缘或界面区在当前轨迹中更容易移动。比较apo与复合物、WT与突变体时，可定位某个处理是否改变了局部波动轮廓。

读图时先看：先检查高峰位于无序末端还是功能区域，再比较同一残基编号、同一原子集和同一生产期的分布。应关注可重复的区域性变化，而不是单个尖峰。若研究界面，最好把界面残基、催化残基或门控loop直接标在横轴上，并报告重复轨迹的均值和误差。

不能这样下结论：高RMSF不等于不稳定，也不自动代表配体诱导变化；低RMSF也可能是体系被困住或约束过强。没有对照时，最多只能描述柔性分布。

最合适的联用方式：与DSSP判断高波动是否伴随二级结构改变，与关键距离和侧链二面角判断功能几何，与B-factor比较实验和模拟的柔性轮廓。

更稳妥的结果表述：配体结合后，口袋门控loop的RMSF在重复轨迹中一致下降，而其他结构域变化有限，提示结合主要限制了局部运动，而非使整个蛋白普遍刚化。

5. B-factor：把MD柔性与实验结构联系起来

易被忽视｜适合做实验-模拟对照，不是RMSF换单位这么简单

图5 实验晶体、结构恢复与MD换算B-factor的比较，原论文Figure 3，来源见参考文献[4]。

它到底算的是什么：在各向同性近似下，MD涨落可换算为B= 8π²RMSF²/3。公式看似只是单位变换，但实验晶体B-factor还混合了热运动、静态无序、晶格接触、分辨率、精修模型和占有率等因素；溶液MD的B-factor则依赖力场、温度、对齐和采样。

它真正回答的问题：它最适合回答模拟是否重现实验结构中“哪些区域更活、哪些区域更稳”的相对轮廓，而不是要求绝对数值逐点一致。若某些高B-factor区同时在MD中表现为高RMSF，可支持这些区域具有真实柔性；明显不一致则可能提示晶体接触或模型条件差异。

读图时先看：优先比较峰位、排序和相关系数，并明确Cα、主链还是侧链。必要时对模拟B-factor做线性缩放或标准化，但必须说明处理。还要检查实验结构中是否存在缺失残基、低占有率和晶体接触。

不能这样下结论：不能把实验B-factor当作纯热运动，也不能因为MD与晶体绝对值不等就判定力场错误。B-factor本身不告诉你柔性变化的方向和功能后果。

最合适的联用方式：与RMSF、DSSP、晶体接触分析及NMR/HDX等溶液实验联用。若研究突变，优先比较柔性轮廓是否系统改变，而不是比较某一个残基的绝对值。

更稳妥的结果表述：模拟与晶体B-factor在主要柔性峰位上具有一致趋势，但绝对幅度不同；这支持区域性柔性分布的一致性，不代表两种条件下原子位移完全等价。

6. 回转半径Rg

条件性指标｜只有研究整体紧致或展开时才重要

图6 回转半径Rg时间序列，完整保留原论文Figure 3C的坐标轴与四组轨迹，来源见参考文献[5]。

它到底算的是什么：Rg是所有选定原子相对质心的质量加权均方距离平方根，描述质量分布的空间尺度。它对整体旋转和平移不敏感，但对所选原子、缺失片段、寡聚状态和结构大小敏感。

它真正回答的问题：它回答蛋白或聚合物是否发生整体膨胀、收缩、塌缩或展开。对于球状单体，持续上升并同时伴随SASA增加和二级结构丢失，才较有力地提示展开；对于多结构域蛋白，Rg变化也可能来自结构域相对排布，而非折叠核心破坏。

读图时先看：看长期趋势、分布和重复轨迹，而不是只比两条曲线的均值。不同大小的蛋白绝对Rg不可直接横比；同一蛋白不同体系的小差异，也应报告误差和效应量。若研究局部口袋，Rg通常过于粗糙。

不能这样下结论：Rg低不等于更稳定，0.01-0.03 nm的差别更不能脱离波动和重复性写成“结构明显更紧致”。紧致也可能来自非天然塌缩。

最合适的联用方式：与全蛋白SASA、DSSP、接触图和结构域距离联用。对聚合物或无序蛋白，还应结合端到端距离、形状张量和标度关系。

更稳妥的结果表述：各体系Rg分布高度重叠，未见持续膨胀或塌缩，因此不能据此声称某一体系显著更紧致；局部差异需由口袋或结构域指标判断。

7. 全蛋白SASA

条件性指标｜看总体溶剂暴露，不看单个口袋

图7 全蛋白SASA时间序列，完整保留原论文Figure 3B的坐标轴与四组轨迹，来源见参考文献[5]。

它到底算的是什么：SASA是以一定探针半径滚过范德华表面得到的可溶剂接触面积。结果依赖探针半径、原子半径和所选原子。全蛋白SASA把疏水、极性、口袋、表面和界面全部加在一起，是一个全局量。

它真正回答的问题：它回答蛋白总体暴露面积是否系统改变。SASA持续上升可能对应展开或界面打开，下降可能对应压紧、埋藏或聚集；但方向本身没有好坏，需要结合结构背景。按极性/非极性或结构域分组后，往往比单一总值更有解释力。

读图时先看：先看趋势和分布，再定位变化来自哪些残基或区域。蛋白-蛋白复合物可进一步计算埋藏表面积；膜蛋白则应区分水相和脂相暴露。不同大小体系不能用绝对SASA直接排序。

不能这样下结论：SASA有波动不等于不稳定，总SASA下降也不自动等于结合更强。全蛋白SASA不能替代配体SASA、界面SASA或口袋水化分析。

最合适的联用方式：与Rg和DSSP判断整体折叠，与残基SASA定位暴露变化，与界面接触和口袋体积解释局部埋藏。

更稳妥的结果表述：总SASA无持续漂移，说明整体溶剂暴露保持在相近范围；口袋是否封闭以及配体是否埋藏仍需使用局部SASA和几何指标判断。

8. 配体SASA与埋藏程度

易被忽视｜判断配体是否向溶剂迁移很有价值

图8 配体性质时间序列与边缘分布；其中SASA位于倒数第二行。原论文Figure 12B，来源见参考文献[6]。

它到底算的是什么：配体SASA描述配体有多少表面可被溶剂探针接触。更严格的界面埋藏可用复合前后面积差计算，如SASAprotein + SASAligand - SASAcomplex；不同文献对是否除以2、是否使用同一构象有不同定义，必须写清。

它真正回答的问题：它回答配体在模拟中是持续深埋、部分暴露，还是逐渐进入溶剂。对于RMSD升高但是否解离不清楚的体系，配体SASA是很实用的第二证据。它还能显示同一口袋内由深结合转为浅结合的状态。

读图时先看：观察SASA变化是否与质心距离、接触数和口袋开合在同一时间发生。比较不同配体时，可使用暴露比例或相对埋藏面积，减少分子大小偏差。极性很强或长尾配体即使稳定结合，也可能保持较高SASA。

不能这样下结论：SASA低只能说明更埋藏，不等于自由能一定更有利；埋得深还可能付出更大的去溶剂化代价。不同大小、不同柔性的配体不能只靠绝对SASA排名亲和力。

最合适的联用方式：与配体RMSD、口袋距离、疏水接触、水桥和能量分析联用。若研究界面形成，应优先报告埋藏面积而非仅报告复合物总SASA。

更稳妥的结果表述：配体SASA在姿势转换后升高，同时口袋接触减少、质心距离增大，支持其由深埋状态转向更暴露的浅结合状态。

9. 帧间RMSD矩阵

易被忽视｜比单条RMSD更能看见状态与回访

图9 上三角为帧间RMSD矩阵，下三角为GROMOS聚类归属，原论文Additional file 6，来源见参考文献[7]。

它到底算的是什么：帧间RMSD矩阵计算任意两帧i与j之间的RMSD，得到Dij。它不再只把所有帧与初始结构比较，而是直接比较轨迹内部每个时间点之间的相似性，因此可以显示状态块、转换和回访。

它真正回答的问题：沿对角线出现连续低RMSD方块，表示一段时间内构象彼此相似；多个方块表示多个状态；相隔很远的时间段仍彼此相似，说明体系曾回访同一状态。它能揭示单条RMSD平台内部是否其实包含多个不同构象。

读图时先看：先读色标，再看块状结构、块与块之间的距离以及状态出现的时间顺序。若后半程只停留在一个新方块，而未回到此前状态，说明发生了转换，但不能据此证明平衡。矩阵结果同样依赖对齐、原子集和RMSD阈值。

不能这样下结论：高RMSD颜色只表示两帧差异大，不表示“坏”或高能；一个整齐方块也不自动等于充分采样，它可能只是体系长期困在一个状态。

最合适的联用方式：与聚类、PCA、FEL和状态时间线联用；用独立重复检查相同状态块是否再次出现。对于配体，也可单独构建配体或口袋残基的帧间矩阵。

更稳妥的结果表述：轨迹可分为两个内部相似、彼此差异明显的构象块，且后期未见充分往返，因此更适合描述为一次状态转换，而不是已经达到全局收敛。

二、再判断配体或界面到底发生了什么

第10-19项把问题从“蛋白有没有大变化”推进到“配体是否留在口袋、功能几何是否维持、相互作用如何替换”。对于蛋白-配体研究，这一组通常比多画几条全局稳定性曲线更关键。

10. 配体-口袋质心距离与关键几何

核心指标｜直接把分析连接到科学问题

图10 八组关键原子距离随时间变化，完整保留原论文Figure 7的全部距离面板，来源见参考文献[5]。

它到底算的是什么：质心距离描述两个原子组的整体相对位置；关键几何则包括催化原子距离、供受体距离、配位距离、进攻角、门控残基间距或蛋白界面锚点。相比通用RMSD，它们直接对应具体结构假设。

它真正回答的问题：整体距离回答配体是否仍在目标区域，局部距离和角度回答功能几何是否成立。对于酶、金属配位、共价反应前构象或门控通道，关键几何通常比全局RMSD更接近研究问题。

读图时先看：至少选择一个整体位置指标和2-3个机制锚点，并在轨迹中同步检查。距离应在处理周期性边界后计算，原子定义不能中途改变。对于催化构象，距离合适但角度不合适仍可能是非反应构象，因此应看二维分布或联合占有率。

不能这样下结论：一对原子距离稳定不能代表整个姿势稳定：配体可能围绕该锚点旋转。距离阈值也不是普适常数，应根据作用类型、实验结构和力场定义。

最合适的联用方式：与配体RMSD、接触指纹、关键侧链二面角和代表构象联用。若要证明某个催化几何更常出现，应报告满足全部几何条件的帧比例，而不是分别报告几个平均距离。

更稳妥的结果表述：配体整体位置保持不变，但催化距离和进攻角仅在部分时间窗同时满足要求，说明体系保留结合，却只间歇访问反应就绪构象。

11. 口袋体积、入口宽度与门控

易被忽视｜研究口袋重塑时比全蛋白Rg更直接

图11 MMP-2、MMP-3与MMP-9结合口袋体积随时间的变化，原论文Figure 3，来源见参考文献[8]。

它到底算的是什么：口袋体积由探针、网格或几何算法识别空腔并计算大小；入口宽度通常由若干门控残基或结构元素间距表示。结果依赖口袋种子、探针半径、网格间距、是否包含配体以及对瞬时水分子的处理。

它真正回答的问题：它回答结合位点是否发生开合、塌缩、扩张或通道门控，以及这些变化是否与配体姿势、突变或别构状态同步。对于口袋重塑、耐药突变和进入/退出路径，它往往比全蛋白RMSD和Rg更敏感。

读图时先看：固定算法和口袋定义，比较体积分布而非只比较平均值；寻找多峰分布或平台转换，并回到对应代表构象观察由哪些侧链或loop驱动。口袋体积增加可能来自入口打开，也可能只是某个侧链旋转后算法识别了相邻空腔。

不能这样下结论：口袋更小不等于结合更强，更大也不等于更容易结合。不同软件、探针和初始种子的绝对体积不能直接横向比较。

最合适的联用方式：与入口距离、关键侧链χ角、配体SASA、水占有、聚类和状态时间线联用。若体积变化是主结论，应展示至少两个状态的三维口袋而非只给曲线。

更稳妥的结果表述：突变体口袋体积分布转向较大区间，并伴随门控残基χ角和入口距离同步改变，支持局部口袋重塑；该结果本身不等同于亲和力升高。

12. 最小距离与接触数量

核心指标｜快速识别界面保持与完全分离

图12 每帧蛋白-配体总接触数时间序列，原论文Figure 6D，来源见参考文献[2]。

它到底算的是什么：最小距离是两个原子组之间最接近的一对原子的距离；接触数则统计低于指定阈值的原子对或残基对。两者都依赖原子选择、是否排除氢原子、距离阈值以及按原子还是按残基计数。

它真正回答的问题：最小距离适合快速发现蛋白与配体是否完全分开，接触数适合描述界面的总体接触丰富度和随时间的重排。稳定的最小距离加上持续接触，支持两者仍处于同一结合区域。

读图时先看：先报告阈值，再看接触数是否发生持续下降、突变式变化或在不同平台间切换。大配体天然产生更多接触，因此跨分子比较应归一化或比较共同锚点。总接触数不变时，还要检查接触身份是否已完全替换。

不能这样下结论：接触越多不等于结合越强，短暂碰撞也会被计数。最小距离很小只说明某一对原子接近，不能保证整体姿势正确。

最合适的联用方式：与按残基接触占有率、IFP、配体RMSD和SASA联用。蛋白-蛋白体系还应增加界面埋藏面积、盐桥和界面水分析。

更稳妥的结果表述：复合物在生产阶段始终保持近距离接触，但接触身份发生两次系统重排，提示界面未解离，却在不同结合亚态间切换。

13. 氢键数量

基础概览｜看网络规模，不识别具体锚点

图13 每帧蛋白-配体氢键数量，原论文Figure 11，来源见参考文献[3]。

它到底算的是什么：氢键数量是在每一帧按给定供体-受体距离和角度条件统计满足几何要求的氢键条数。不同软件可能使用不同距离、角度和氢原子定义，因此数值只在同一协议下可比。

它真正回答的问题：它回答界面在每一时刻存在多少极性接触，以及网络是相对连续还是频繁断裂重建。它适合做概览，尤其能识别配体离开后氢键数整体归零的时间段。

读图时先看：看分布、平均值和时间模式，不看最大值。平均始终为2条，可能是同两条氢键一直存在，也可能是十几条不同氢键轮流替换；这两种机制完全不同。还应区分蛋白-配体氢键、配体内部氢键和蛋白内部氢键。

不能这样下结论：氢键多不等于亲和力强。极性基团从水中进入口袋会付出去溶剂化代价，而计数不包含这一点。不能用某一帧的最高氢键数作为整条轨迹的证据。

最合适的联用方式：必须进一步统计具体供体-受体对的占有率和寿命，并与水桥、质子化状态、配体SASA和静电/溶剂化能综合解释。

更稳妥的结果表述：界面氢键总数保持在相近范围，但这一结果只能说明极性接触持续存在；具体锚定残基需由逐对占有率和寿命确定。

14. 氢键占有率与寿命

核心指标｜比“平均几条氢键”更接近稳定锚点

图14 按残基与作用类型堆叠的interaction fraction；绿色代表氢键，原论文Figure 6B，来源见参考文献[2]。寿命仍需由连续事件或相关函数另行统计。

它到底算的是什么：占有率是某一供体-受体对满足氢键几何条件的帧比例；寿命描述一次连续存在或允许短暂断裂后再次形成的持续时间。连续寿命与间歇寿命定义不同，必须说明是否允许短暂断键。

它真正回答的问题：占有率用于区分长期锚定作用和短暂接触，寿命则进一步区分“经常形成但每次很短”和“形成后持续很久”。某残基多个原子或多个作用可同时出现，因此按残基汇总的interaction fraction有时会超过100%，这不代表统计错误。

读图时先看：列出供体、受体、阈值、统计窗口和重复轨迹。对关键氢键同时画距离和角度分布，检查高占有率是否来自合理几何。多个中等占有率氢键可能构成可替换网络，其功能不一定弱于单条高占有率氢键。

不能这样下结论：高占有率不等于该氢键对结合自由能贡献最大，也不能忽略质子化、互变异构和水竞争。不同软件的占有率不能不加说明直接比较。

最合适的联用方式：与氢键总数、水桥、残基接触、MM/PB(GB)SA分解和突变验证联用。真正的“关键氢键”最好同时满足持续、几何合理、重复可见和化学机制相关。

更稳妥的结果表述：AspX与配体形成的直接氢键在多条重复轨迹中均保持较高占有率，并表现出较长连续寿命，支持其作为持续极性锚点；其能量重要性仍需其他方法验证。

15. 盐桥距离与占有率

机制指标｜先确定质子化状态，再谈离子作用

图15 盐桥距离及0/1占有时间线，完整保留原论文Figure 5A，来源见参考文献[9]。

它到底算的是什么：盐桥通常指带相反形式电荷的基团在一定距离内形成的离子接触。可用带电原子最小距离、羧酸氧到碱性氮的距离或两组电荷中心距离定义；不同定义会改变占有率。

它真正回答的问题：它回答两个带电基团是否长期靠近、是否发生断裂重建，以及突变、pH或离子强度是否改变了离子网络。盐桥在蛋白界面、核酸结合和门控网络中常有明确结构意义。

读图时先看：在计算前确认Asp/Glu、His、Lys/Arg、末端和配体的质子化状态。看距离分布是否有单峰或双峰，再以阈值计算占有率。羧酸两个氧之间频繁切换时，按单一原子统计会低估持续性，宜按基团定义。

不能这样下结论：若质子化状态错了，盐桥统计再漂亮也没有物理意义。盐桥占有率高也不等于净自由能一定有利，因为电荷去溶剂化和离子屏蔽可能抵消直接静电作用。

最合适的联用方式：与pKa/质子化评估、氢键和水桥、离子分布及静电-极性溶剂化分量联用。研究pH效应时，固定质子化的常规MD只代表预设状态，必要时需常pH MD。

更稳妥的结果表述：在所设质子化状态下，盐桥距离呈稳定近距离分布并具有较高占有率；该结论限定于当前pH模型和离子条件。

16. 残基接触频率与占有率

核心机制指标｜总接触数看规模，接触频率看身份

图16 四条独立重复中配体-残基疏水接触频率及配体占据体积，原论文Figure 2的完整相关面板，来源见参考文献[10]。该图以疏水接触为例展示“某残基在多少比例帧中形成接触”。

它到底算的是什么：接触频率通常定义为某一残基对或某一残基-配体关系满足预设接触条件的帧数，占统计总帧数的比例，即f = Ncontact/Nframes。接触可以按任意重原子距离、Cα/Cβ距离、范德华接触或特定疏水规则定义。它与“每帧总接触数”不是同一个量：总接触数告诉你界面当时有多少接触，接触频率告诉你究竟是哪一对残基在多长比例的轨迹中持续出现。

它真正回答的问题：它回答哪些残基是持续锚点、哪些只是偶尔碰撞、接触网络是否在不同亚态之间替换，以及WT/突变体、apo/复合物或不同重复轨迹是否保留同一组核心接触。蛋白-蛋白界面可统计残基对频率；蛋白-配体体系可统计每个口袋残基与配体的频率；疏水接触频率只是其中一种常见子类型。

读图时先看：先明确接触对象、原子级还是残基级、距离/几何阈值和生产期。逐残基曲线要映射回三维结构，并比较独立重复中的峰位是否重现；残基对矩阵则要寻找持续块和状态特异接触。单个残基若可同时接触配体多个原子，按相互作用条目累加后的interaction fraction可能超过1；而按“该残基本帧是否至少接触一次”定义的占有率应在0到1之间。

疏水接触应怎样包含在这里：疏水接触通常按非极性原子或疏水残基与配体碳原子的距离统计，适合描述口袋包围和形状互补。跨配体比较时应考虑分子大小，可比较共同核心或归一化接触比例；若要讨论疏水驱动，还需结合配体SASA、埋藏非极性面积和口袋水，而不能只数碳原子接触。

不能这样下结论：高接触频率只说明几何关系经常出现，不等于该残基对自由能贡献最大，也不等于突变后一定显著降低亲和力。MD相邻帧高度相关，90%的帧占有并不等于90个独立证据；阈值改变也会改变绝对频率。

最合适的联用方式：与最小距离/总接触数、IFP、配体RMSD、SASA和代表构象联用。筛选候选热点时，再与MM/PB(GB)SA残基分解或实验突变交叉验证；对疏水接触，还应结合范德华打包和口袋水分析。

更稳妥的结果表述：在统一接触定义下，残基X、Y和Z在多条独立轨迹中均保持较高接触频率，并映射到同一口袋核心，支持它们构成持续接触网络；这一统计描述接触保留，不单独等同于结合能热点。

17. π-π与阳离子-π作用

机制指标｜距离、夹角与偏移必须同时满足

图17 芳香环相互作用的双参数几何模型：环心距离d与法向夹角。原论文Figure 1，来源见参考文献[11]；MD中应逐帧计算并统计占有率。

它到底算的是什么：π-π作用需描述两个芳香环质心距离、法向夹角和横向偏移，可分近平行堆积和T形作用；阳离子-π作用则涉及带正电基团与芳香环中心及法向的几何关系。仅靠最近原子距离不足以定义。

它真正回答的问题：它回答芳香识别是否以合理几何持续存在，以及不同配体或突变是否改变了环的取向和作用类型。对芳香口袋、核酸碱基堆积和某些多酚体系，这类几何常是关键机制变量。

读图时先看：绘制质心距离与夹角的联合分布，或按严格几何计算占有率；对不同作用模式分别统计。还应检查环是否因对称翻转造成RMSD升高，但π几何仍保留。阳离子-π分析需确认阳离子的质子化和力场电荷。

不能这样下结论：静态对接图中两个芳香环靠近，不等于100 ns内存在稳定π-π作用。单一距离阈值会把不合理取向误算为堆积，经典力场对极化效应的描述也有限。

最合适的联用方式：与芳香侧链二面角、IFP、配体占据体积和能量/量化计算联用。若π作用是论文核心结论，建议给出二维几何分布而不只给最终帧截图。

更稳妥的结果表述：芳香环质心距离多数时间处于合理范围，但夹角分布显示平行与T形模式交替，因此更适合描述为动态芳香识别，而非单一固定堆积。

18. 水桥占有率、口袋水密度与驻留

易被忽视｜直接氢键消失后，水可能仍在维持识别

图18 口袋目标水合位点及分块水占有率，完整保留原论文Figure 8A-B，来源见参考文献[12]。

它到底算的是什么：水桥指一个或多个水分子同时满足蛋白-水和水-配体的氢键几何。桥接占有率关注某个位置是否经常由水介导；水密度图关注哪些空间位置常被水占据；驻留时间则关注某个水分子或水位点的交换动力学。三者不是同一概念。

它真正回答的问题：它回答直接氢键断开后是否存在水介导补偿、哪些结构水位点被保留或置换，以及口袋开合是否伴随水进入或排出。对极性口袋、金属周围、核酸结合和配体选择性，水网络经常是决定性却被忽略的部分。

读图时先看：区分“同一个水分子长期不动”和“同一空间位点不断由不同水交换占据”。高水桥占有率未必伴随长驻留时间。应报告水模型、几何阈值、统计窗口，并比较晶体水位点或多个重复轨迹。

不能这样下结论：存在水桥不自动等于结合增强；有些水提供稳定桥接，有些水的释放反而有利于结合。仅统计水桥数量无法判断热力学方向。

最合适的联用方式：与直接氢键占有率、配体SASA、三维水密度、驻留相关函数和口袋体积联用；若水是核心机制，可进一步采用GIST、WaterMap类分析或自由能方法。

更稳妥的结果表述：直接氢键占有率下降时，同一极性位点保持较高水桥占有，说明识别网络由直接作用转为水介导补偿；其自由能贡献仍需专项水分析。

19. 相互作用指纹IFP

易被忽视｜把“作用很多”变成可比较的动态模式

图19 按残基与作用类型绘制的相互作用指纹时间线，完整保留原论文Figure 5a及其图例，来源见参考文献[13]。

它到底算的是什么：IFP把每一帧中某残基与配体是否形成氢键、疏水、离子、π作用或水桥编码成离散特征，随后统计时间线、占有率或指纹相似度。它不是新的相互作用力，而是一种组织多类接触的表示方法。

它真正回答的问题：它回答哪些共同锚点在多个体系中保留、哪些作用只在某个配体或某个状态出现，以及姿势转换时相互作用网络怎样整体替换。相比单独列十几条作用，IFP更适合比较候选配体和重复轨迹。

读图时先看：先确认每类相互作用的几何定义，再看持续条带、协同出现和状态转换。低占有率作用不一定无意义，它可能是进入路径或替换网络的一部分；高占有率也要回到三维结构检查是否合理。比较多体系时应使用相同残基编号和同一规则。

不能这样下结论：总作用数或指纹相似度不能直接排序亲和力。IFP是几何和分类统计，不包含去溶剂化、熵和长程能量，也会受到阈值选择影响。

最合适的联用方式：与配体RMSD、聚类和状态时间线联用：先用IFP识别相互作用状态，再提取各状态代表构象；对多个重复轨迹报告均值和变异。

更稳妥的结果表述：两种配体共享同一组核心锚点，但第二种配体额外形成稳定水桥并减少姿势切换，说明其动态识别模式不同；是否亲和力更强仍需能量或实验支持。

三、只有机制问题明确时，再解析构象空间

第20-26项从主链φ/ψ、DSSP和关键侧链rotamer推进到PCA、DCCM、聚类与FEL，用于识别局部构象变化、集体运动、协同关系和状态人口。它们并不是所有项目的标准套餐；没有明确残基、状态变量和对照体系时，PCA、DCCM与FEL很容易只剩一张好看的图。

20. Ramachandran图：主链φ/ψ构象分布

经典但常被漏掉｜看主链扭转角采样，不是只做结构质检

图20 MD轨迹中蛋白主链φ/ψ角的二维Ramachandran密度图，原论文Figure 9，来源见参考文献[24]。图中颜色和等高线表示采样密度。

它到底算的是什么：Ramachandran图把每个非末端残基的主链二面角φ与ψ投影到二维平面。经典结构验证图常统计单个结构中残基是否落在favored、allowed或outlier区域；用于MD时，则可把轨迹中大量帧的φ/ψ角汇总成散点、密度或自由能式分布，并可针对全部残基、某类残基或某个关键残基分别作图。Gly和Pro的允许区域与普通残基不同，通常应单独处理。

它真正回答的问题：它回答主链实际采样了哪些构象区域、某个关键残基是否在α、β、polyproline II或其他构象盆之间转换，以及突变、配体或环境是否使局部主链人口发生偏移。它能补足DSSP的离散分类：DSSP告诉你某帧被分为哪类二级结构，Ramachandran图则直接展示构成这种分类的连续φ/ψ分布。

读图时先看：先区分静态结构质检图和轨迹分布图。若研究MD，不能只报告“最终结构有多少残基位于favored region”；应比较同一残基或同一残基集合在生产期的二维分布、主要峰位置、状态人口和重复轨迹。全蛋白汇总图可能被大量稳定螺旋残基淹没，真正的loop转换应单独画关键残基的φ/ψ时间序列或二维分布。跨体系比较必须采用相同残基集合、分箱和统计窗口。

不能这样下结论：轨迹中出现少量outlier不自动等于模拟失败，loop、Gly、末端和过渡态本来可能访问低人口区域；反过来，大多数点位于favored region也不能证明全蛋白稳定或采样收敛。若力场把主链过度限制在某一区域，图看起来“很集中”反而可能暴露采样偏差。

最合适的联用方式：与DSSP、局部RMSF、关键主链二面角时间序列和代表构象联用。若某个φ/ψ跃迁与功能有关，应进一步检查它是否同步改变氢键、门控距离、接触网络或催化几何，并在独立重复中重现。

更稳妥的结果表述：关键loop残基X的φ/ψ分布由单一α样区域转为α样与β样双峰，并与DSSP状态切换和门控距离变化同步，提示该残基参与局部主链构象转换；该结论不等同于全蛋白发生变性。

21. DSSP二级结构时间图

核心结构指标｜定位变化发生在哪些残基

图21 DSSP残基-时间图，不同颜色编码二级结构类别，原论文Figure 4，来源见参考文献[14]。

它到底算的是什么：DSSP依据主链氢键和几何为每个残基逐帧分配二级结构类别，如α螺旋H、3₁₀螺旋G、β链E、β桥B、turn T、bend S和coil/none。颜色只是类别标签，不代表能量或稳定性高低。

它真正回答的问题：它回答某段螺旋或β结构是否保留、在哪个时间点发生转换、变化是短暂呼吸还是持续重排。相较只给α/β总含量，残基-时间图能直接定位变化位置。

读图时先看：先读图例和类别合并规则，再寻找沿时间方向连续成片的变化。零星一两帧切换常是分类边界抖动；连续残基在较长时间内从H变为coil，才更支持局部解螺旋。对膜蛋白、酶活性loop或突变附近，应单独标出功能区域。

不能这样下结论：任何颜色变化都不等于变性，coil也不是“坏结构”。DSSP描述局部主链几何，不说明变化是否功能有利，也不能单独证明全蛋白稳定。

最合适的联用方式：与RMSF、氢键、Rg/SASA和代表构象联用。若某段二级结构变化是主结论，应在独立重复中重现，并检查是否受末端截断、质子化或膜环境影响。

更稳妥的结果表述：突变附近连续残基由螺旋转为coil并在后半程持续存在，同时局部RMSF升高，支持局部二级结构重排，而非由单帧分类波动造成。

22. 关键侧链二面角与rotamer状态

易被忽视｜全蛋白RMSD不变，口袋侧链仍可能换挡

图22 Glu150的χ1、χ2侧链二面角、对应盐桥距离以及Phe153主链ψ角随时间变化，原论文Figure 5，来源见参考文献[25]。

它到底算的是什么：侧链二面角χ1、χ2等决定残基rotamer。多数侧链并不是连续均匀旋转，而是在gauche+、trans和gauche-等构象盆之间停留与跳变。分析可用角度时间序列、圆周分布、二维χ1/χ2分布或rotamer占有率。由于角度具有-180°/180°周期，普通算术平均可能产生错误，例如179°与-179°的平均不应被解释为0°。

它真正回答的问题：它回答口袋侧链是否发生翻转、某个盐桥或氢键为何突然形成/断开、入口门控残基是否在开闭rotamer间切换，以及突变或配体是否改变关键侧链状态人口。许多机制事件只涉及少数侧链，蛋白骨架RMSD和Rg几乎不变，因此这类指标往往比全局曲线更直接。

读图时先看：先从结构假设中选择真正相关的残基，不宜把所有侧链都画一遍。标出rotamer区间和状态切换时间，并同步查看与其相连的距离、接触或口袋体积。对于Phe、Tyr、Asp、Glu等具有对称或近等价原子的侧链，应注意角度定义和对称性；统计分布时使用圆周方法，比较重复轨迹中的状态人口和转换是否一致。

不能这样下结论：单次χ角翻转不等于功能转换，它可能只是正常热涨落；停在一个rotamer也不自动说明稳定，可能是采样时间不足。不能把不同残基的角度绝对值直接比较，也不能用线性均值处理跨越周期边界的角度。

最合适的联用方式：与关键距离、氢键/盐桥占有率、口袋体积、IFP和聚类联用。若某一rotamer被认为控制功能，应展示它与功能几何的联合分布或条件占有率，而不是只给一条角度曲线。

更稳妥的结果表述：门控残基X的χ1由gauche+转入trans状态后，入口距离与配体SASA同步升高，且该状态在主要构象簇中重复出现，支持侧链rotamer切换参与口袋开放；单次翻转本身不作为独立机制证据。

23. 主成分分析PCA

机制专项｜看最大方差方向，不是直接看自由能

图23 MD位移向量投影到前两个主成分后的构象空间，原论文Figure 8a，来源见参考文献[15]。

它到底算的是什么：PCA对去除整体平移和转动后的坐标协方差矩阵进行特征分解。特征向量给出集体运动方向，特征值给出沿该方向的方差。PC1、PC2只是数据驱动坐标，不是预先定义的反应坐标，也不是能量轴。

它真正回答的问题：它回答当前轨迹中哪些集体运动贡献了最大方差，以及不同状态在这些方向上如何分离。将特征向量映射为箭头或极端构象，可解释PC对应结构域开合、loop摆动还是整体扭转。

读图时先看：报告对齐方式、原子集和前几个PC解释的方差比例。比较apo与复合物时，最好把合并轨迹投影到同一协方差基底，或定量比较子空间重叠；若各自单独做PCA，两个体系的PC1方向通常并不相同，散点范围不能直接比较。

不能这样下结论：PC1数值大不等于不稳定，散点更集中也不自动等于能量更低。PCA偏向大振幅线性运动，慢过程不一定就是最大方差方向。

最合适的联用方式：与FEL和聚类提取状态，与结构动画解释运动方向，与重复轨迹检查主子空间是否稳定。若关注最慢动力学过程，可考虑tICA或Markov状态模型，而不是只做PCA。

更稳妥的结果表述：在共同PCA基底上，两体系沿与结构域开合相关的PC1占据区域不同，提示构象集合发生偏移；该差异仍需由状态人口和重复轨迹验证。

24. DCCM动态交叉相关矩阵

易被误读｜相关不等于因果，也不等于促进或抑制

图24 DCCM残基-残基相关矩阵；正负相关必须按原图色标判断，原论文Figure 13，来源见参考文献[3]。

它到底算的是什么：DCCM通常计算两个原子位移向量的归一化点积Cij，范围-1到+1。接近+1表示沿相同方向协同位移，接近-1表示相反方向位移，接近0表示线性同向/反向耦合较弱。矩阵对称，主对角线为自身相关。

它真正回答的问题：它回答哪些残基或结构域在当前采样中表现为协同或反向运动，是否出现远距离耦联和可能的别构通信区域。远离主对角线的连续色块通常比零散像素更值得关注。

读图时先看：必须先读色标，不能默认红正蓝负。不同体系应使用相同对齐、原子集、时间窗和色标范围。全蛋白对齐方式会影响相关模式，短轨迹也可能产生偶然相关；应使用分块或重复轨迹检查色块是否稳健。线性DCCM还可能漏掉非线性或相位错开的耦合。

不能这样下结论：正相关不是“促进结合”，负相关不是“破坏稳定”，DCCM也不提供信息传播方向和因果关系。仅凭一张相关矩阵不能建立别构通路。

最合适的联用方式：与PCA、距离/二面角状态、广义相关或互信息以及动态网络分析联用，并把显著区域映射回三维结构。若声称别构通路，应有突变、实验或更严格因果分析支持。

更稳妥的结果表述：配体结合后两个远端结构域间的相关模式发生可重复改变，提示动态耦联重组；这一结果提出了候选别构联系，但不单独证明信号传递方向。

25. 构象聚类与代表构象

机制专项｜把数万帧归纳成少数可解释状态

图25 SOM构象聚类、簇占比与代表构象，原论文Figure 7，来源见参考文献[7]。

它到底算的是什么：聚类按照选定特征和距离度量把相似帧分组。特征可以是主链RMSD、口袋原子、配体姿势、接触指纹、二面角或降维坐标；代表构象宜选簇中心附近的真实帧medoid，而不是坐标平均后可能失真的“平均结构”。

它真正回答的问题：它回答轨迹中有哪些主要状态、各状态人口和出现时间，以及哪个真实构象可以代表后续口袋、相互作用或能量分析。它能把“RMSD有波动”转成具体的状态转换。

读图时先看：必须报告算法、特征、距离和阈值。检查簇数对阈值是否敏感，主要簇是否在重复轨迹中出现，以及状态人口是否被某一条轨迹支配。人口最高只表示在当前采样和定义下最常出现，不自动等于全局稳定态。

不能这样下结论：聚类结果不是唯一答案；换原子集或阈值，簇数可能明显改变。不能把一个100 ns轨迹中占比最大的簇称为“唯一稳定构象”。

最合适的联用方式：与帧间RMSD矩阵、PCA/FEL和IFP联用。先用与问题相关的特征聚类，再对每个主要状态分别统计口袋、作用网络和关键几何。

更稳妥的结果表述：轨迹可稳健归纳为三个主要结合亚态，其差别主要来自配体取向和门控侧链，而非蛋白整体折叠；各亚态人口需结合重复轨迹报告。

26. 自由能景观FEL

机制专项｜它是概率投影，不是独立于轨迹的新证据

图26 基于前两个主成分构建的自由能景观及各能谷代表构象，原论文Figure 5a，来源见参考文献[16]。

它到底算的是什么：常规MD中的FEL通常按

把所选坐标上的采样概率换算为相对自由能。低能盆本质上就是高概率区域。它依赖坐标选择、分箱或核密度参数、温度和轨迹是否接近平衡。

它真正回答的问题：它回答在所选一到两个坐标投影下，哪些区域最常访问、是否出现多个概率盆，以及候选状态之间是否存在低概率区。它适合组织构象状态，不应被当成与轨迹人口无关的额外能量计算。

读图时先看：优先选择与研究问题有关且能区分状态的坐标，而不是机械使用PC1-PC2。检查换坐标、换分箱和分块后盆的位置是否稳健，并从各盆提取代表构象。若轨迹没有发生充分往返，盆深和盆间差值可能严重依赖初始条件。

不能这样下结论：一个深蓝低能盆不等于充分收敛或全局最低能构象；二维投影可能把隐藏在其他自由度中的状态叠在一起。常规短MD上的FEL也不宜直接当作定量能垒。

最合适的联用方式：与聚类、状态时间线、重复模拟和增强采样联用。若要定量讨论能垒，应使用合适反应坐标和自由能方法，并提供收敛与往返证据。

更稳妥的结果表述：在所选两个构象坐标上，轨迹主要分布于两个概率盆；由于状态往返有限，该图用于描述候选亚稳态，而不将盆深解释为已收敛的定量能垒。

四、能量放在证据链最后

第27项把MM/PBSA或MM/GBSA的总能、统计误差、能量分量和残基分解合并为同一个主分析；第28项则讨论沿指定反应坐标的PMF。两者回答的问题不同，且都不能替代位置、接触和构象证据。

27. MM/PBSA与MM/GBSA：一个主分析，四层输出

能量专项｜总能、收敛、分量和残基分解不应拆成四个平行指标

图27A 端点结合自由能的running average与cumulative average随模拟时间变化，原论文Figure 8A，来源见参考文献[17]。

图27B GGAS、GSOLV与TOTAL等能量分量的比较，原论文Figure 8A-B，来源见参考文献[18]。

图27C 逐残基MM-PBSA贡献与候选热点在结构表面的映射，原论文Figure 9，来源见参考文献[19]。

它到底算的是什么：MM/PBSA和MM/GBSA都属于端点能量估算：从复合物、受体和配体的分子力学能与连续介质溶剂化能估算结合自由能，区别主要在极性溶剂化模型使用Poisson-Boltzmann还是Generalized Born。常见单轨迹协议可降低受体和配体内部能的噪声，但也隐含结合前后构象变化较小的假设。所谓总结合能、时间/分块收敛、能量分量和残基分解，并不是四种彼此独立的方法，而是同一套MM/PB(GB)SA分析的四个输出层次。

它真正回答的问题：第一层是总能与体系间相对趋势，回答在同一协议下哪个体系给出更有利的端点能量估计；第二层是时间依赖、分块、重复轨迹和误差，回答这个排序是否被某一小段轨迹或少数极端帧主导；第三层是范德华、气相静电、极性与非极性溶剂化等分量，回答差异在哪些计算项中体现以及它们如何抵消；第四层是逐残基分解，用于提出候选接触热点。后面三层是解释和质控总结果，不应被包装成三个额外的“结合能结论”。

读图时先看：先证明比较对象在生产期仍处于可比的结合状态，再选定一致的轨迹区间、帧抽样、力场、PB/GB模型、介电常数、离子强度和单/多轨迹协议。统计上应估计自相关或采用足够大的时间块，报告分块结果、独立重复差异和置信区间。累计均值变平并不能单独证明收敛，因为后期新数据对累计值的影响天然越来越小；必须同时看独立区块与replica。

能量分量应该怎样解释：常见分量包括范德华、气相静电、极性溶剂化和非极性溶剂化，总估计还可能包含构象熵。先核对公式、符号、单位和分量求和，再看抵消关系。静电项很负时，往往伴随不利的极性去溶剂化；范德华项更负可以提示打包更紧，但不能单独写成“疏水作用决定亲和力”。若省略熵，结果更接近焓样端点评分，而不是完整标准结合自由能。

残基分解应该怎样解释：残基分解把部分能量项按残基归属，用于从持续接触中筛选候选热点。应同时检查均值、误差、分块/重复稳定性、接触频率和三维位置。极性溶剂化与多体效应不能总被唯一分配给单个残基，因此分解值不是把该残基突变为Ala后的ΔΔG，也不能把最负残基直接称为已经验证的关键位点。

不能这样下结论：只取最后10 ns得到一个很负的数值，不足以证明收敛；同一轨迹抽取数千帧也不会自动产生数千个独立样本。MM/PB(GB)SA通常更适合同一协议下的相对比较，不宜直接当作精确实验绝对亲和力。带电配体、金属、显著构象重排和高度柔性体系对参数与采样尤其敏感，小于误差或重复间波动的差异不能过度解释。

最合适的联用方式：与配体RMSD/关键距离、接触频率、IFP、水化、主要构象状态和实验结果联用。若要验证残基突变效应或很小的相对亲和力差异，应考虑突变模拟、FEP/TI或实验，而不是继续拆更多MM/PBSA图或增加小数位。

更稳妥的结果表述：在相同协议下，体系B的MM/GBSA总能在独立分块和重复轨迹中均更有利，差异大于统计不确定性；分量分析显示范德华优势被部分极性溶剂化代价抵消，残基X与Y同时具有持续接触和较有利分解值。因此该结果作为相对能量和候选热点支持，不解释为精确绝对结合自由能或已验证突变ΔΔG。

28. 伞形采样PMF

高级自由能分析｜反应坐标、窗口重叠和隐藏自由度决定可信度

图28 同一体系沿不同解离路径得到的伞形采样PMF，原论文Figure 12，来源见参考文献[20]。

它到底算的是什么：伞形采样在一系列反应坐标窗口中施加偏置势，再用WHAM/MBAR等方法重构沿该坐标的PMF。曲线给出的是该坐标上的相对自由能，不是脱离路径所有自由度的完整描述。

它真正回答的问题：它回答结合、解离、膜穿越或构象转换沿指定路径的稳定区和相对能垒。与端点方法相比，它能描述过程，但只有当反应坐标捕捉关键慢自由度、各窗口充分采样并互相重叠时才可靠。

读图时先看：除PMF曲线外，必须查看窗口直方图重叠、不同时间段的收敛、误差以及正反方向或不同初始路径的一致性。若某窗口存在正交慢运动或姿势翻转，单纯延长所有窗口未必解决问题。对配体解离，还要注意方向、取向和标准态修正。

不能这样下结论：曲线最高点不自动等于真实动力学过渡态，PMF也不直接给出解离速率。把结合态最低点与远端平台简单相减，若缺少体积、取向和约束修正，通常不能直接称为标准结合自由能。

最合适的联用方式：与窗口重叠图、误差、二维反应坐标、不同路径和代表构象联用。若关注速率，需要milestoning、MSM、加速动力学或其他动力学方法，而非只看PMF高度。

更稳妥的结果表述：PMF显示沿所定义路径存在相对自由能垒，且窗口重叠与分块结果支持曲线稳定；该能垒限定于当前反应坐标，不直接解释为实验解离速率。

最后：怎样把指标组织成证据链？

不是先把20多项全部算完，再从中挑几张“好看的图”；而是先写出要证明的句子，再为这句话配置相互独立的证据。下面几条可以直接作为项目分析的最小逻辑。

配体是否稳定保留：蛋白对齐后的配体RMSD + 质心/锚点距离 + 接触/IFP + 配体SASA。四类结果同向，才适合讨论“仍处于口袋”。

突变是否重塑口袋：局部RMSF + 口袋体积/入口距离 + 关键侧链二面角 + 接触网络 + 代表构象。全蛋白RMSD和Rg只能作背景。

是否发生功能构象转换：与功能直接相关的距离或角度 + DSSP/局部结构 + 聚类或共同基底PCA + 状态时间线；FEL只负责组织状态人口。

是否存在别构耦联：局部扰动证据 + 远端构象变化 + DCCM/广义相关或动态网络 + 多条重复轨迹；最好再由突变或实验验证。

是否支持亲和力差异：先证明结合姿势和采样可比，再给MM/PB(GB)SA均值、误差、分块与重复；小差异不要过度解读，必要时升级到FEP/TI。

是否支持解离能垒：合理反应坐标 + 窗口重叠 + PMF误差和收敛 + 不同路径/方向验证；PMF能垒不直接等于实验速率。

全文最重要的一句话：指标不是越多越好。每一项分析都应该回答一个明确问题，并且它的结论边界要写清楚。真正有说服力的MD报告，往往图不一定最多，但位置、接触、构象和能量之间能够互相验证。

一个更实际的分析顺序

1.先做轨迹预处理：处理周期性边界、去除整体平移转动、确定生产期和对齐原子。

2.先看与研究问题最接近的几何量，而不是先堆RMSD、RMSF、Rg和SASA。

3.把时间序列转成分布、占有率、状态人口和误差，避免只描述“有波动”。

4.用代表构象把数值变化映射回三维结构，确认指标对应真实可解释的结构事件。

5.最后做能量，并检查分块、相关性和重复轨迹；能量只支持前面已经成立的结构机制。

参考文献

[1] Al-Anazi M, Al-Najjar BO, Khairuddean M. Structure-Based Drug Design Studies Toward the Discovery of Novel Chalcone Derivatives as Potential EGFR Inhibitors. Molecules. 2018;23:3203. 

[2] Harini K, et al. Ligand Docking Methods to Recognize Allosteric Inhibitors for G-Protein-Coupled Receptors. Bioinformatics and Biology Insights. 2021;15:11779322211037769. 

[3] Hussain M, et al. Computational modeling of cyclotides as antimicrobial agents against Neisseria gonorrhoeae PorB porin protein. Frontiers in Chemistry. 2024;12:1493165. 

[4] Liu N, Mikhailovskii O, Skrynnikov NR, Xue Y. Simulating diffraction photographs based on molecular dynamics trajectories of a protein crystal. IUCrJ. 2023;10:16-26. 

[5] Pirolli D, et al. Insights from Molecular Dynamics Simulations: Structural Basis for the V567D Mutation-Induced Instability of Zebrafish Alpha-Dystroglycan. PLOS ONE. 2014;9:e103866. 

[6] Ali A, et al. Network Pharmacology Integrated Molecular Docking and Dynamics to Elucidate Saffron Compounds Targeting Human COX-2 Protein. Medicina. 2023;59:2058. 

[7] Fraccalvieri D, et al. Conformational and functional analysis of molecular dynamics trajectories by self-organising maps. BMC Bioinformatics. 2011;12:158. 

[8] Durrant JD, de Oliveira CAF, McCammon JA. POVME: an algorithm for measuring binding-pocket volumes. Journal of Molecular Graphics and Modelling. 2011;29:773-776. 

[9] Mohanty P, Chatterjee KS, Das R. NEDD8 Deamidation Inhibits Cullin RING Ligase Dynamics. Frontiers in Immunology. 2021;12:695331. 

[10] Bhati AP, et al. Long Time Scale Ensemble Methods in Molecular Dynamics: Ligand-Protein Interactions and Allostery in SARS-CoV-2 Targets. Journal of Chemical Theory and Computation. 2023;19:3359-3378. 

[11] Brylinski M. Aromatic interactions at the ligand-protein interface: implications for the development of docking scoring functions. Chemical Biology & Drug Design. 2018;91:380-390. 

[12] Ge Y, et al. Enhancing Sampling of Water Rehydration on Ligand Binding: A Comparison of Techniques. Journal of Chemical Theory and Computation. 2022;18:1359-1381. 

[13] Ivanova A, Mokshyna O, Polishchuk P. StreaMD: the toolkit for high-throughput molecular dynamics simulations. Journal of Cheminformatics. 2024;16:123. 

[14] Haji-Allahverdipoor K, et al. Investigation of the new substitution glycine to alanine within the Kringle-2 domain of reteplase: a molecular dynamics study. BioTechnologia. 2024;105:201-213. 

[15] David CC, Jacobs DJ. Principal component analysis: a method for determining the essential dynamics of proteins. Methods in Molecular Biology. 2014;1084:193-226. 

[16] Maisuradze GG, Liwo A, Scheraga HA. Principal component analysis for protein folding dynamics. Journal of Molecular Biology. 2009;385:312-329. 

[17] Chung MKJ, Miller RJ, Novak B, Wang Z, Ponder JW. Accurate Host-Guest Binding Free Energies Using the AMOEBA Polarizable Force Field. Journal of Chemical Information and Modeling. 2023;63:2769-2782. 

[18] Budha B, et al. Structure-guided discovery and characterization of novel FLT3 inhibitors for acute myeloid leukemia treatment. PLOS ONE. 2025;20:e0334415. 

[19] Byun J, Lee J. Identifying the Hot Spot Residues of the SARS-CoV-2 Main Protease Using MM-PBSA and Multiple Force Fields. Life. 2022;12:54. 

[20] You W, Tang Z, Chang CEA. Potential Mean Force from Umbrella Sampling Simulations: What Can We Learn and What Is Missed? Journal of Chemical Theory and Computation. 2019;15:2433-2443. 

[21] Knapp B, Ospina L, Deane CM. Avoiding False Positive Conclusions in Molecular Simulation: The Importance of Replicas. Journal of Chemical Theory and Computation. 2018;14:6127-6138. 

[22] Grossfield A, et al. Best Practices for Quantification of Uncertainty and Sampling Quality in Molecular Simulations. Living Journal of Computational Molecular Science. 2019;1:5067. 

[23] Genheden S, Ryde U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opinion on Drug Discovery. 2015;10:449-461. 

[24] Margreitter C, Oostenbrink C. MDplot: Visualise Molecular Dynamics. The R Journal. 2017;9:164-186. 

[25] Salsbury FR Jr, Yuan Y, Knaggs MH, Poole LB, Fetrow JS. Structural and Electrostatic Asymmetry at the Active Site in Typical and Atypical Peroxiredoxin Dimers. Journal of Physical Chemistry B. 2012;116:6832-6843. 

（如有疏漏欢迎交流指正。）