Gromacs——教程学习(6)

12. 谈谈怎么判断分子动力学模拟是否达到了平衡

在计算RMSD之前必须先通过最小二乘法将各帧结构相对于参考结构进行最大程度叠合，从而消除体系的整体运动而令RMSD只体现生物分子内部结构的变化，这称为align或者least squares fit。

需要注意的是，蛋白质骨架的RMSD曲线并不总能反映出蛋白结构特征变化的全貌。一方面，RMSD的计算涉及到一大堆原子的平均，若一个很大的蛋白某个局部出现一定程度的显著的、不可忽视的结构变化，在RMSD计算过程中它会被很大程度淹没掉。

另一方面，骨架的RMSD不能反映侧链上出现的一些关键变化。例如上面这篇文章发现，在423 ps的时候磷酸化酪氨酸574残基侧链出现了显著的翻转，这使得被研究的酪氨酸激酶活性位点打开，这在分子生物学角度是一个很重要的现象，而从alpha碳的RMSD曲线上则完全体现不出来这一点。文中图4(A)绘制了这个残基侧链二面角随时间的变化，这才让这关键信息被展现出来。

如果蛋白质的RMSD曲线在已经跑的模拟时间内还没有稳定在一个值附近且持续了较长时间（有零点零几nm的上下波动是可以接受的），而你又需要蛋白质已达到平衡的结构，通常就需要延长模拟时间续跑，直到RMSD基本平稳。

有些人试图用衡量小分子液体的方式衡量生物分子模拟体系是否达到稳定，即还用温度、密度、势能等参数来衡量，这明显不行。因为这些量的收敛速度远快于蛋白质的RMSD的收敛，在这些指标已平稳的时候蛋白质部分通常还明显没达到平衡。虽然蛋白质整体结构的进一步变化必然会对势能等因素有多多少少的影响，但在势能曲线上基本体现不出来。通常模拟蛋白质都是在显式水模型下模拟的，势能绝大部分都是水所贡献的，显然蛋白质结构的改变对总势能的影响通常会被淹没；

下面两张图是计算化学公社论坛上我答疑时看到的，可见中途RMSD都出现过剧烈的变化，完全超越了平衡状态下的正常波动的程度。这种突变必须给予明确的解释，要不然审稿人看见肯定会问你是怎么回事。始终要记住，RMSD曲线是轨迹变化的定量反映，想解释抽象的RMSD曲线的突然变化是什么导致的，当然要去用VMD程序看轨迹动画才能充分弄清楚原因。绝对不可能RMSD曲线出现厉害的突变时在轨迹动画上却看不出任何端倪。

上图这种情况只要在VMD里观看轨迹动画自然能明白原因。就算不看轨迹也100%能断定肯定是因为周期性方式记录轨迹，而被计算RMSD的对象在模拟过程中越过了盒子所致。当某一帧突然体系被弄到了盒子另一头，RMSD自然会出现瞬间的大幅变化。解决这种RMSD曲线不连续性的方法很简单，消除被计算RMSD的对象的整体运动，并保持此对象的完整性就完了。

比如在GROMACS里，可以用gmx trjconv结合-pbc mol使得转换出的新轨迹中的分子保持完整，计算RMSD前再经过叠合来消除蛋白质的整体运动，其RMSD曲线就不可能突变。还有一种情况是被计算RMSD的对象是由多个分子组成的，比如有多条链的蛋白、蛋白+配体复合物、含有两条链的核酸，这种情况可能模拟中途只有其中一部分跑出盒子而另一部分还在盒子内，光是用-pbc mol保持单分子完整的话此对象的不同部分还是分家的，对这种情况需要用gmx trjconv结合-pbc cluster让指定的组在轨迹中一直保持完整，组的定义就对应于你要计算的RMSD的对象。顺带一提，对生物分子的模拟，我一般建议让程序在模拟过程中就一直对生物分子消除平动转动，对于GROMACS跑蛋白质来说就是设comm-grps = protein和comm-mode = angular，这样就避免了生物分子中途跑出盒子的问题。

不同的生物分子体系的RMSD达到平稳时RMSD的数值会有一定不同。柔性越大的体系，平衡时RMSD越大。比如DNA体系骨架柔性比一般较刚性的蛋白质要大，这会体现在达到平衡时RMSD相对较大上。

通过RMSD考察是否平衡不是必须得用初始帧做参考结构。J. Chem. Phys., 109, 10115 (1998)认为以初始结构当参考帧来观察RMSD会高估达到平衡的时间，因此会把一部分其实已经达到平衡的轨迹误认为还没平衡，导致浪费了一些轨迹。作者的观点是相对于初始结构，模拟开始后RMSD的增加由两部分组成：(1)非平衡态到达平衡态 (2)平衡态时在相空间采样。实际上(1)满足了的时候本质上生物分子就已经平衡了，这比起以初始结构为参考帧计算的RMSD曲线达到平稳时更早。为了判断什么时候就已经满足了(1)，作者建议以已经平衡了的结构作为参考帧来计算RMSD。下图是文中的图，作者取500、600、700、800、900 ps分别作为参考帧绘制了RMSD曲线，并假定这些时刻体系都已经达到了平衡。由图可见，每条曲线在参考帧附近都是迅速上升，这体现出已平衡的结构在相空间采样期间由于热运动自然而然地相对于参考结构发生的结构改变；而在图的前200 ps，这几个曲线都明显上升，因此体现出前200 ps是确确实实结构没有达到平衡、有整体性变化的阶段。

所以，像温度、势能、密度那些很容易达到平稳的量，只是体系达到平衡的必要而非充分条件，不能只拿这些量说事，除非是本文第1节说的均匀的小分子液体之类的简单体系。
势能只是能量角度上的一个标准，但凡有可能，最好还是再结合一个结构方面的标准来进一步判断平衡，明显会更有说服力。