通过string-method(弦理论)搜索蛋白质构象之间最小能量路径(MEP)

1.背景相关

弦方法(string-method)现在已经在几种常见的自由能计算插件中都已经实现了基本的功能，目前为止常见的自由能插件有colvars和plumed，plumed相比于colvars而言可以和各大主流的分子动力学模拟程序并用进行重要性采样计算，而colvars目前主要集成于NAMD，GROMACS和LAMMPS之中，且个人认为colvars的计算效率相对于plumed来说更加快，在自己的笔记本上使用gromacs搭配colvars插件对于一个10万左右的原子体系进行3反应坐标的增强采样的时候，GPU的利用率维持在百分之70左右(在进行cMD的过程中)，每隔500000步使用colvars进行计算平均耗时不超过1秒，所以效率还是可以的。在这里还会涉及到使用自由能计算的时候的续跑问题，这里我就简单通过gromacs(colvars)做一个简单的示例，也欢迎各位和我讨论和NAMD搭配colvars插件的续跑问题，由于我对于plumed用的不多，所以在此就不做主要介绍了，使用gromacs进行自由能计算时的续跑命令如下：

gmx mdrun -v -deffnm eabf -cpi eabf.cpt -colvars colvars_Normal_dihedral.dat -colvars_restart eabf.restart.colvars.state.dat -ntmpi 1 -ntomp 24 -pin on

在这个命令中我们可以看到的是只需要同时加上colvars和colvars_restart两个参数，colvars后面本身加的配置文件就是为colvars编写的输入配置文件，如果要进行续跑，则两个参数必须要同时出现，colvars_restart后面则应该加上模拟所产生的restart文件，由于在gromacs中colvars_restart后面只能加上后缀为.dat的文件，所以必须要将eabf.restart.colvars.state文件拷贝一份为.dat的格式，当你按照如上的，命令进行续跑的时候一般会出现restart中所展示的截止反应坐标和当前帧的反应坐标数值相差过大的问题，此时需要你同时修改eabf.restart.colvars.state.dat和eabf.restart.colvars.state两个文件中对应反应坐标的数值。具体情况可以参考这一个colvars的官方issue：

issue#535:Restarting simulation with colvars

在我们接下来要讲的内容中，我将结合一个经典的示例分子和一个蛋白质构象变化体系具体讲一下如何计算由state1转变成state2的有意义的转换路径，寻找一条最有可能的最小作用路径，一般来讲，都是取两个状态之间沿着自由能面的一条直线猜测路径，在估计计算沿着这一条路径上的自由能变化，在这一个过程中我们将会使用到gspath和gzpath两个反应坐标并将其应用到自适应偏置力重要性采样算法之中，也就是ABF，之前在知乎上看到有人回答说的是string-method并不支持ABF而是支持metadynamics，可能是因为软件版本不同的原因，目前使用NAMD或者gromacs都可以实现string-method方法。

2.String-Method的理论背景

如上图所示为一张仅供测试的自由能景观图，在此图中可以看到的是有一条沿着state1至state2的一张自定义构象转换路径，在这张图上可以划分几个构象点，每一个点则可以为，在这里我们需要描述一组集体变量（反应坐标)或者是粗变量：,上图中的公式则是将表示为了一组反应坐标或者是粗略变量，这里就需要简要介绍一下这个对象的具体含义了，通常在物理化学或者统计力学之中，一般来说是系统在微观状态下在某一个时间节点上的函数，可以用于描述系统的宏观性质变化，这里的n的数值大小一般远小于N，N定义为系统的原子总数，所以这个很方便理解，其实也就是在自由能景观图上自定义的一个初猜路径，那么对于自定义的变量所对应的平均力势即PMF则可以通过下述公式进行表示：。在这个公式中为系统中一个给定微观状态的玻尔兹曼因子，则是对应状态的势能，在在这里可以将理解为上述粗变量集合公式中所对应的, 和两者等价, 是热力学因子，是玻尔兹曼常数，是系统的绝对温度。有意思的是在分子 中会观察到一个函数，此函数在数学物理中具有特殊的意义，也很有意思, 定义在0点的积分为一且在0点处趋向于无穷大但是在0以外的所有地方都等于0，笔者记得自己第一次接触到函数的时候也很疑惑，可以理解为一个广义函数或者分布，所以在分子中起到了一定的粗变量筛选的作用，接下来我们在回过去细讲一下分子的具体含义，分子表示了对于所有可能的微观状态的 或者 的积分，这个积分用于计算在特定下的所有可能配置的总权重。然后分母则是所有可能微观状态下的玻尔兹曼因子的积分，代表了整个配分函数。既然一般用于描述路径上的某一个状态。所以相应的函数在分子处则可以理解为一个过滤函数, 函数本身具有正交性和筛选性，其筛选性则主要体现在下述的公式之中：

其筛选性可以通过函数的极限性质证明可得，其实反过来想如果将理解为反应坐标或者是定义的反应坐标的集合表示会更加合理一点。还有就是当第一次见到上面的第二个公式的时候可能会有很多人有疑问这个是不是一个导数，将其理解为状态的变化量，但是在我个人的理解看来应该是一个不同于的函数，也就是说的目的适用于筛选特定的​值，可以将其抽象的理解为某一个位置或者原子坐标。回到正题，当我们假设定义的反应坐标可以在自由能表面随机游走，基于布朗运动的形式，就可以通过时间变化量来离散化运动，下面的这个公式则是描述了反应坐标随时间变化的动力学模型，一般可以用于描述分子动力学模拟过程中反应坐标是如何随时间演化的:

在上述公式中, 是在时间之后反应坐标的值，则是在初始时间(t=0)的时候反应坐标的值, 这一部分内容中是由温度因子加权的扩散张量与保守力的乘积之和。这部分内容主要是代表了体系中的粒子或者原子由于保守力而产生的确定性移动，这里在接下来我会按照我的理解来主要解释一下扩散张量和保守力的相关概念，这两个词汇一般主要出现在理论物理与分子动力学之中，其中扩散张量，这个比较常见的地方会出现在DTI中(不是我们长认为的Drug Target Interaction，一般为字面意思直译为Diffusivity Tensor)，即扩散张量成像技术，其主要是一个度量分子在各个方向上扩散能力的张量。在各向同性介质中，分子在所有方向上的扩散能力是相同的，扩散可以用一个标量（即扩散系数）来描述。然而，在各向异性介质中，分子在不同方向上的扩散能力不同，因此需要用一个张量来完整地描述这种扩散行为。扩散张量的每一个元素描述了粒子在方向上的扩散如何受到方向上扩散的影响。例如，在各向异性材料或者复杂流体中，分子在某一方向上的运动可能会影响或者被另一方向上的运动影响。关于保守力，主要是一个物理力学上的概念，指的是做功只与物体的初始位置和最终位置有关，而与物体运动的路径无关的力。保守力可以由一个标量势能函数导出，即力是势能的负梯度。在分子动力学模拟中，保守力通常指的是由分子间相互作用势能（如Lennard-Jones势、库仑势等）导出的力。例如，如果一个分子在势能场中移动，那么在任意位置上，保守力指向的是势能减小最快的方向。在动力学方程中，保守力会对原子的运动产生决定性的贡献，而扩散张量主要描述了原子由于热运动和介质异质性所导致的随机扩散过程，那么在式中会发现出现了扩散张量和保守力的乘积，那么乘积一般也具有特定的物理意义，例如在上述的等式中，则是代表了原子由于梯度势能所经历的定向运动的强度，同时也代表了原子是如何响应体系中的势能变化的，也影响着原子在自由能表面上的移动方向。因此公式中的第二部分发现主要描述了不同粗变量，回过头我们再来看一下公式中的第三部分是扩散张量关于反应坐标的梯度，表示扩散率是如何随着反应坐标而发生变化的，最后一项为随机噪声，特性是高斯分布，这部分也主要代表了由于热涨落而导致的随机运动。

因此基于上述的内容和自由能景观图，如果我们想要探究A点(假设)到B点之间的演化过程，那么我们首先需要生成一条初猜路径，这一条初猜路径上会分布很多的点，在这里我将其称为过渡点，当然有别的叫法只要合理应该都没有什么问题，然后根据每一个点生成一批轨迹群，从这一批轨迹值中我们计算得到了一个对应过渡点的平均偏移量，所以说如果你的初猜路径中包含有个点，然后每个点所对应的轨迹群的数量为，假设每一条轨迹的模拟时间尺度为，那么截至到目前为止的内容我们针对于一个体系所需要跑的模拟时间尺度就为，后续还会涉及到迭代收敛问题，为了将初猜路径收敛到最小能量路径，从一个任意路径开始，将会对每一个初猜路径中的过渡点进行演化，直到其在垂直于初猜路径上的垂直分量变化是0为止，因为在这样子的迭代过程中我们需要保持对应的过渡点之间的距离不变，同样在这样子的一个过程中会有一些点会趋向于自由能低点运动，这也是目前string-method的一个主要问题之一。刚才在上述的内容中我们提到了平均偏移量这个概念，同时这个变量也是导致初猜路径中过渡点运动转化的主要影响因素，下面的这一个公式表示了平均偏移量的主要构成:

这个公式其实和反应坐标集合随着时间变化的公式十分相像，其中在时间内反应坐标集合的变化量，则是反应坐标在初始时刻所对应的值。这个平均漂移量可用于更新和优化弦方法中的路径，使其逼近最小作用路径。在实践中，这意味着通过在弦的每个影像点上基于无偏轨迹集合计算出的平均偏移，我们可以逐步迭代地改进初始路径，直至它收敛到最小作用路径。这个过程通常涉及重复的模拟和分析，以确保路径最终符合所求的最优路径。所以在这里我们设置每一个过渡影像点需要迭代的次数为，结合上述的内容，那么string-method的需要跑的总时间尺度为，因此一般来讲没有足够算力的情况下不要尝试使用string-method，并且相比于Dijkstra算法，string-method对于涉及到周期性变量的MEP结算结果并不是特别好。

接下来涉及到pathcv，即路径集合变量，pathcv其实主要分为两个部分, , 这两部分所对应的公式如下：

在这之中为均方根偏差RMSD的平方，所以在定义pathcv的时候还会涉及到rmsd部分的编写，详细的编写细节可以参考plumed或者colvars插件手册，其实和定义别的反应坐标的过程差不多相差不大，是比较简单的。

上述的大部分内容基于自己的理解，如果有理解错误的地方欢迎大家指正;)