GNM和ANM的渊源

        弹性网络模型 (elastic network model, ENM)将蛋白质的每个残基简化为一个Node, 用其Cα来代替, 残基之间的相互作用简化为谐振势, 用弹簧来代替, 并且所有弹簧的强度系数都相同。ENM完全忽略了残基之间精细的相互作用, 将传统原子层次的NMA方法简化到残基层次, 大大降低了计算量和内存消耗。Bahar、Hinsen等人的研究表明采用该残基水平上的简化模型计算得到的蛋白质大幅度运动模式与传统原子水平上的NMA方法的计算结果能够很好地吻合, 说明蛋白质的整体运动模式由其拓扑结构所决定, 而不是由残基之间精细的相互作用所决定。在Bahar等人所提出的模型中, 残基在平衡位置附近的振动认为是各向同性的高斯型的运动, 该模型称为高斯网络模型 (Gaussian network model, GNM)。GNM能够计算得到不同运动模式中各个残基的运动幅度, 但是无法获得其运动的方向。随后, Atilgan等人将GNM进行了扩展, 考虑了残基运动的方向性信息, 将各向同性模型发展为各向异性模型, 建立了各向异性网络模型 (anisotropic network model, ANM) 。

GNM的基本原理

        GNM将蛋白质的三维空间结构简化为一个弹性网络。蛋白质中的每个残基简化为一个点, 用其Cα来代替, 残基之间的相互作用简化为谐振势, 用弹簧来代替。为了判断残基之间是否存在相互作用, 选取一个距离截断半径 (一般取0.7~0.8nm) , 如果2个残基的Cα之间的距离小于截断半径, 则认为它们离得足够近, 存在相互作用, 用一弹簧相连, 反之, 则不存在相互作用。并且, 在GNM中所有弹簧的强度系数都相同。基于上述简化方法, 体系的势能等于所有弹簧的势能之和, 可以写成

        式中：γ为弹簧的强度系数；ΔRij为残基i和j之间距离的涨落；ΔRi和ΔRj分别表示残基i和残基j偏离各自平衡位置的涨落；N为蛋白质体系内残基的个数；ΔRT={ΔR1, ΔR2, …, ΔRN}表示各个残基偏离平衡位置的涨落, 上标T表示转置；Γ为Kirchhoff矩阵, 其矩阵元为

        式中：Rij为残基i和j之间的平衡距离；rc为距离截断半径。

 在GNM中, 由于体系的势能为谐振势, 体系的配分函数可以解析求解, 即

        式中: kB为玻尔兹曼常数；T为绝对温度。

       残基涨落的概率分布可以写成

根据统计物理学理论, 可以计算物理量的统计平均值。残基i和残基j之间涨落的交叉关联可通过以下公式进行计算：

残基i的均方涨落可以写为：

        通过上面的讨论可以发现, 上述物理量的计算都与Γ-1有关。Γ-1可以采用矩阵分解的方法进行求解, 即

        式中：U为正交矩阵, U的每一列uk (1≤k≤N) 为一个特征向量, 表示体系的一个运动模式；Λ-1为对角矩阵, 每一个对角元λk (1≤k≤N) 为体系的一个特征值, 表示相应运动模式的频率。特征值可以按照从小到大进行排列, 即0 = λ1<λ2<…<λN, 其中, λ1 = 0对应于体系的整体平动。

        一般情况下, 关心的是蛋白质体系的内部运动, 因此, 常常将λ1 = 0所对应的运动舍去。根据能量均分定理, 频率越小的运动模式对应的运动幅度越大, 因此, 特征值比较小的低频慢运动模式对应于大幅度的整体运动, 常常表示体系的功能性运动。而特征值比较大的高频快运动模式一般表示体系局部结构的不规则性, 常常对应于体系的关键位点。根据矩阵分解, Γ-1可以表示为所有特征向量的叠加, 即

        式中k = 1, 表示体系的整体平动, 舍去不考虑.通过Γ-1进而可以计算配分函数和物理量的统计平均值。

ANM的基本原理

        在GNM中, 残基在平衡位置附近的运动是各向同性的, 不区分其运动的方向性, 因此, GNM可以计算得到体系构象运动的幅度, 但是无法得到其运动的方向。为了解决这一问题, Atilgan等人对GNM进行了扩展, 把残基运动在x、y、z三个方向进行了区分, 建立了ANM。相比于GNM, 在ANM中, 3N×3N维的Hessian矩阵取代了原来的N×N维的Kirchhoff矩阵, Hessian矩阵的矩阵元为

        式中：1≤i, j≤N；V为体系的势能函数。采用与GNM类似的计算方法, 可以通过矩阵分解计算Hessian矩阵的逆矩阵, 获得体系的各个运动模式以及相应的频率, 进而, 基于统计物理学理论计算体系的配分函数和所关心的物理量的系综平均值。

个人理解的几个基本概念：

1. Mode：

意为模式，计算GNM或ANM中间过程得到的Kirchhoff矩阵或Hessian矩阵在进行特征值和特征向量分解时，以N×N维的Kirchhoff矩阵为例，会得到N个特征值（在弹性网络模型中定义为频率或强度）和N个对应的特征向量（在弹性网络模型中定义为Fluctuation Vector，即波动向量）。一个Mode即为一个特征值和一个特征向量的对应组合，此时的特征向量简单而言是一个N×1维的矩阵（N行1列），第 i 行的元素代表在某一个Mode下，第 i 个残基的波动情况。波动情况的说明如下图所示。

2. Motion:

意为运动，是一个或者多个Mode的集合。一般而言，前几个Mode对于蛋白质的动力学性质贡献最大。在实际应用中，需要具体问题具体分析，ProDy进行分析时，一般默认取前20个Mode，可根据自身需要进行调整。

3. ENM:

弹性网络模型，是用生物物理的方法来处理分析蛋白质的有效手段，可以理解为类似胡克弹性定律的一种变相拓展，从网络的高维度层面去表征蛋白质的内在动力学性质。具体的原理和GNM的推导在本文中已有表述，ANM 的推导鉴于涉及偏导数的内容，此处不再赘述。ENM 所涉及的具体后续分析在 ProDy 的官网均可找到对应教程，建议阅读 Ivet Bahar 实验室已发表的文献深入了解弹性网络的分析流程，从而理解作图结果背后隐藏的生物学意义。