Rotamers & Packer¶

@Author: 吴炜坤

@email：weikun.wu@xtalpi.com/weikunwu@163.com

一、旋转异构体(Rotamer)¶

最早在有机化学的课本上就着重介绍过旋转异构体的概念，经典的例子是一个乙烷分子的不同二面角状态下能量的分布曲线，由于氢原子(红色标记)存在位阻效应，在重叠式状态下能量最高，而交叉式的能量最低。当能垒足够高的时候，构象异构体之间就不能轻易互相转换，呈现出几个极小值。具有极小值势能的几个构象互为旋转异构体(Rotamer)。

蛋白质中的氨基酸侧链构像也呈现出离散型的分布, 对于每一个sp3-sp3杂化间为中心的$\chi$二面角，都有g-(-60°),t(180°),g+(+60°)三种离散低能状态存在（类似乙烷分子的交叉式和重叠式）。

例如，谷氨酸的$\chi_{1}$二面角分布中，CG-CB原子形成的二面角是sp3-sp3杂化键为中心的，因此也有3种低能状态，但与上述的乙烷分子不同，CA原子上相连是C=O、NH、H三种不同的基团，这些基团有着不同的体积（位阻效应不同）。对于$\chi_{1}$角来讲g-,t,g+三种二面角的出现概率是不相等的。换言之，谷氨酸侧链$\chi_{1}$角的分布是依赖于蛋白骨架二面角的（backbone-dependent rotamer），因为骨架上NH、C=O基团的朝向是由$\phi$和$\psi$角，不同的骨架二面角具有不同的位阻效应。

对于sp3-sp3杂化键为中心的$\chi$二面角构象通常是对称的，也称为具有rotameric性质的。既然氨基酸侧链的$\chi$角状态分布是可数的，那么将每个$\chi$角的状态组合起来，就可以得到若干个离散的能量构象。如对于谷氨酸而言， 𝜒1 和 𝜒2 是rotameric的，分别有3种状态g-,t,g+。那组合起来便有9种状态，分别是<g-,g->, <g-,t>, <g-,g+>, <t,g->, <t,g+>, <t,t>, <g+,g->, <g+,t>, <g+,g+>。每一种Rotamer的出现概率是依赖于蛋白骨架二面角，如果将蛋白骨架二面角$\phi$和$\psi$以每10°x10°作为统计的单位区间，就可以得到Rotamer出现的概率分布图(共计9个)。

细心的读者可能会发现，上述例子中的谷氨酸最后一个$\chi_{3}$二面角的中心是sp3-sp2杂化的，CG原子上连接的是羧基团平面。这种二面角是Non-rotameric的。Non-rotameric二面角往往是非对称性的，其分布状态受到先前原子二面角构象的影响。如对于谷氨酸而言，$\chi_{1}$和$\chi_{2}$是rotameric的，构象异构体共有9种状态，$\chi_{3}$二面角的分布式概率依赖于9种状态进行统计的。如下图中直观展示，谷氨酸$\chi_{1}$$\chi_{2}$处于<g+,g->状态下，$\chi_{3}$二面角的概率分布:

二、Rotamer Library¶

在之前的章节中，我们提及过Rosetta使用MCMC的方法对骨架构象进行采样。但是在全原子模型下，氨基酸残基的侧链构象的建模是尤为重要的。不同的侧链构象具有不同的能量，侧链建模本质上就是根据Rotamer能量的高低进行构象优化和搜索。Rotamer Library在这个过程中即负责记录每种Rotamer构象在特定$\phi$和$\psi$的bin格中的出现概率，随后将概率取-log对数得到Rotamer自身的单体能量(当然在Rosetta中处理稍微更复杂一些)。

以下是Rosetta中内置的Dunbrack Backbone-dependent Rotamer Library中谷氨酸部分数据的展示:

字段的含义:

• Phi, Psi字段记录了骨架二面角$\phi$和$\psi$的bin格的区间;
• r1/r2/r3/4代表$\chi$角的id，如果是rotameric的，r1=1、2或3时，对应的就是g-,t,g+三个区间。以第一行数据为例，谷氨酸的$\chi_{3}$是Non-rotameric，r3=1时，对应的角度是-1.2°。但是r3 id对应的角度不是唯一的，它是根据$\chi_{1}$和$\chi_{2}$的id组合不同而不同的，如第11行中，r1=3，r2=3时，r3同样等于1，但是其对应的角度却是-37.2°;
• Probabil代表了某个Rotamer出现概率，在一个$\phi$和$\psi$的bin格的区间中，所有Rotamer出现概率之和为1;
• chi1Val等代表$\chi_{1}$的二面角平均值;
• chi1Sig等代表$\chi_{1}$的二面角在这个bin格中的分布方差。

在Rosetta中Rotamer的能量计算公式如下，主要由三项组成:

$E_{\text {fa_dun }}=\sum_{r}-\ln \left(P\left(\operatorname{rot}_{r} \mid \phi_{r}, \psi_{r}, a a_{r}\right)\right)+\sum_{k<T_{r}} \frac{1}{2}\left(\frac{\chi_{k, r}-\mu_{\chi_{k}}\left(\phi_{r}, \psi_{r} \mid \operatorname{rot}_{r}, \mathrm{aa}_{r}\right)}{\sigma_{\chi_{k}}\left(\phi_{r}, \psi_{r} \mid \operatorname{rot}_{r}, \mathrm{aa}_{r}\right)}\right)^{2}+$ $-\ln \left(P\left(\chi_{T_{r}, r} \mid \phi_{r}, \psi_{r}, \operatorname{rot}_{r}, \mathrm{aa}_{r}\right)\right)$

第一项是$\sum_{r}-\ln \left(P\left(\operatorname{rot}_{r} \mid \phi_{r}, \psi_{r}, \mathrm{aa}_{r}\right)\right)$，这一项其实就是对一个Rotamer的Probabil取-ln对数，可直接从数据库中获取，直接评估当前rotamer类型的标准化平均能量。

第二项是$\sum_{k<T_{r}} \frac{1}{2}\left(\frac{\chi_{k, r}-\mu_{\chi_{k}}\left(\phi_{r}, \psi_{r} \mid \operatorname{rot}_{r}, \mathrm{aa}_{r}\right)}{\sigma_{\chi_{k}}\left(\phi_{r}, \psi_{r} \mid \operatorname{rot}_{r}, \mathrm{aa}_{r}\right)}\right)^{2}$，该项的作用是对当$\chi$角偏离平均值时做出的能量惩罚/校正(能量升高)，因为一个rotamer的$\chi$角不可能正好等于数据库中的平均值。该项惩罚递增的速率由方差决定。这一惩罚项只对rotameric$\chi$角有效。该函数形状类似高斯分布，$\mu_{\chi_{k}}$指的是chi1Val、chi2Val等为$\chi$角的平均值，$\sigma_{\chi_{k}}$即为chi1Sig等值，代表方差。

第三项是$-\ln \left(P\left(\chi_{T_{r}, r} \mid \phi_{r}, \psi_{r}, \operatorname{rot}_{r}, \mathrm{aa}_{r}\right)\right)$，该项主要对最后一个$\chi$角起作用，主要用于处理Non-rotameric二面角。当Rotamer最后一个$\chi$角是rotameric时，该项恒等于0。反之对最后一个$\chi$角做出相应的能量校正。

综上使用三项既可以分别描述每一个$\chi$角的能量分布曲线，也可以同时用于描述每种Rotamer构型r在整个$\phi$和$\psi$bin中的概率分布，如下效果:

三、Rosetta Packer¶

截止上述，我们已经完整地阐述了在Rosetta中侧链Rotamer的相关概念以及Rotamer Library和能量的计算方式。在这一节，我们将会介绍Rosetta Packer的基本算法以及具体应用方式。

试想下如下场景: 在一个连续的两股螺旋组成的体系中，第i号位点上为精氨酸残基(R)，第j号位点上是谷氨酸残基(E)。现在需要求解i和j位点上构象能量最低的状态？ 一个可能的笨方法是遍历所有i和j位点上的rotamer构象，分别计算每一组rotamer状态下的能量值，最终选取能量最低的两个rotamer即可。

在Rosetta中，侧链优化的基本原理的确如此，每两个氨基酸位点之间，不同Rotamer组合的能量就是如此遍历计算，并储存在一个矩阵(InteractionGraph)当中。不过当氨基酸位点数量的增加时，侧链构象优化的问题将变得极为复杂(NP-hard)。Rosetta Packer算法采用的是mcmc+模拟退火的方法，配合InteractionGraph的预计算的数据，通过查询表的方式，非常快速地得到目标构象的能量值，而不需要每次都反复重新计算，该机制大幅度提升了侧链优化和搜索的速度。

Packer模拟退火的C++代码¶

# Run simulated annealing:
for ( int i = 1; i <= num_outer_iterations; ++i ) {
  for ( int j = 1; j <= num_inner_iterations; ++j ) {
     int newrot = pick_random_rotamer();
     compute deltaE = E( newrot ) - E( oldrot )
     if ( metropolis_accept( temp, deltaE )) {
          accept random_rotamer;
      }
  }
  decrease_temperature( temp );
}

Rotamer构象模拟退火的包括两个部分:

第一个部分是外部的模拟退火控制，根据outer循环的迭代次数不断地降低温度，使得整个蛋白质的构象被“困”在能量的极小值区域。

第二个部分是mcmc采样，pick_random_rotamer函数即随机选择一个位点，并随机生成该氨基酸残基的一个rotamer。接着是熟悉的Metropolis准则，如果当前构象能量下降接受新的构象，反之以一定的概率接受新构象。

完成整个模拟退火过程后，将能量最低的构象作为输出。

Packer的运行模式¶

Packer有三种运行模式，分别是“Repacking”、“Rotamer Trial”和“Design”，这三种模式有着应用上的一定区别。

“Repacking”直译的含义是“重新打包”，顾名思义，需要将以前侧链的状态先“忘记”，重头优化构象。具体的做法是:

1. 将所有的侧链原子删除
2. 逐步添加一些侧链结构，并运行模拟退火
3. 重复这个过程，得到一系列低能量的侧链构象集合

“Rotamer Trial”直译是“尝试Rotamer”，顾名思义，需要在当前的构象上继续优化，优点是快速。具体的做法是:

1. 随机选择一个位点
2. 运行模拟退火搜索能量更低的构象
3. 如此往复，直到所有侧链能量比较高的位点都得到了优化

“Design”直译就是“序列设计”，因此这个模式下可以运行氨基酸类型发生变化。其具体的做法与“Repacking”类似，但最大的不同在于其考虑了所有被允许出现的氨基酸所有的Rotamer构象。这种模式往往需要配合骨架的Relax来避免过快收敛。