结合自由能计算（MM/PB(GB)SA）

分子动力学
分子动力学方案
发布日期:   2023-06-27
更新日期:   2023-07-04
文章字数:   2.5k
阅读时长:   10 分
阅读次数:  

1. 用途

采用MM/PB(GB)SA方法对常规分子动力学轨迹进行结合自由能计算和能量分解，暂不支持计算熵、不支持膜蛋白。

2. 预备知识

MM/PB(GB)SA

MM/PB(GB)SA

计算结合自由能的方法有很多，例如，热力学积分（Thermodynamic Integration，TI）、自由能微扰（Free Energy Perturbation，FEP）、MM/PB(GB)SA、线性相互作用能（Linear Interaction Energy，LIE）。MM/PB(GB)SA方法是精度和速度折衷的方法，广泛应用于受体-配体结合自由能计算。该方法全称 分子力学/泊松-波尔兹曼（广义波恩）表面积（Molecular Mechanics / Poisson Boltzmann (Generalized Born) Surface Area）。顾名思义，该方法将结合自由能拆分成分子力学项和溶剂化能分别计算。

其基本原理是：计算两个溶剂化分子在结合（bound）和游离（unbound）状态的结合自由能之差或者比较同一个分子不同的溶剂化构象的自由能。这里以受体蛋白和配体小分子的结合为例，大致讲述其计算原理。

最简单的思路，莫过于按照下图所示的过程直接计算结合自由能，即可分别计算受体、配体以及它们的复合物各自在溶液中的能量，然后算差值：

ΔGbind0=ΔGcomplex0−(ΔGrecptor0+ΔGligand0)\Delta G^0_{\mathrm{bind}} = \Delta G^0_{\mathrm{complex}} - (\Delta G^0_{\mathrm{recptor}} + \Delta G^0_{\mathrm{ligand}})

但在实际计算中，会遇到一个很严重的问题——能量贡献主要来自溶液间相互作用，并且总能量的波动幅度远大于结合能，这样就需要非常长的时间才能收敛。

为此，通过下面的热力学循环，“兜个圈”来避免这种糟糕的情况：

ΔGbind,solv0=ΔGbind,vacuum0+ΔGsolv,complex0−(ΔGsolv,ligand0+ΔGsolv,receptor0)\Delta G^0_{\mathrm{bind, solv}} = \Delta G^0_{\mathrm{bind, vacuum}} + \Delta G^0_{\mathrm{solv, complex}} - (\Delta G^0_{\mathrm{solv, ligand}} + \Delta G^0_{\mathrm{solv, receptor}})

这个方案的含义是，将溶剂中的总结合自由能分拆成分子力学项（真空中的结合自由能）和溶剂化能两部分分别计算。

  • 分子力学项

该项仍然按照第一个公式计算，即：

ΔGbind,vacuum0=ΔGcomplex,vacuum0−(ΔGreceptor,vacuum0+ΔGligand,vacuum0)\Delta G^0_{\mathrm{bind, vacuum}} = \Delta G^0_{\mathrm{complex, vacuum}} - (\Delta G^0_{\mathrm{receptor, vacuum}} + \Delta G^0_{\mathrm{ligand, vacuum}})

根据结合自由能的定义，对于每一个分子，按照下面公式计算：

ΔGvacumm0=ΔEMM0−T⋅ΔS0=(ΔEint0+ΔEvdw0+ΔEele0)−T⋅ΔS0\Delta G^0_{\mathrm{vacumm}} = \Delta E^0_{\mathrm{MM}} - T · \Delta S^0 = (\Delta E^0_{\mathrm{int}} + \Delta E^0_{\mathrm{vdw}} + \Delta E^0_{\mathrm{ele}}) - T · \Delta S^0

其中，ΔEint0\Delta E^0_{\mathrm{int}}包含了键、键角和二面角能量，ΔEvdw0\Delta E^0_{\mathrm{vdw}}为非键范德华能量，ΔEele0\Delta E^0_{\mathrm{ele}}为非键静电能量，T⋅ΔS0T · \Delta S^0为熵贡献，可以采用简正模分析（normal mode analysis）得到。但实际计算中，通常忽略该部分贡献。因为采用MM/PB(GB)SA方案计算的体系通常为受体-配体结合前后构象变化不大的，这部分贡献可以在计算差值时抵消掉。此外，简正模分析是非常耗时且误差界限（margin of error）很大、会引入显著的不确定度。

  • 溶剂化能

溶剂化能可分成两部分：极性溶剂化能和非极性溶剂化能。

ΔGsolv0=ΔGsolv,polar0+ΔGsolv,nonpolar0\Delta G^0_{\mathrm{solv}} = \Delta G^0_{\mathrm{solv, polar}} + \Delta G^0_{\mathrm{solv, nonpolar}}

对于前者，主要有两种计算方法：Poisson-Boltzmann和Generalized Born方程，计算较为复杂，这里不展开叙述。对于后者，计算较为简单，它正比于溶剂可及表面积（solvent-accessible surface area，SASA），它体现了疏水效应，常用的表达式为：

ΔGsolv,nonpolar0=γ⋅SASA+β\Delta G^0_{\mathrm{solv, nonpolar}} = \gamma·\mathrm{SASA} + \beta

将上述全部公式稍加整理，总结如下：

ΔGbind0=ΔGcomplx0−(ΔGrecptor0+ΔGligand0)\Delta G^0_{\mathrm{bind}} = \Delta G^0_{\mathrm{complx}} - (\Delta G^0_{\mathrm{recptor}} + \Delta G^0_{\mathrm{ligand}})

对于等号右侧的每一项，按照下面公式计算：

ΔG0=(ΔGvdw0+ΔGele0)+(ΔGsolv,polar0+ΔGsolv,nonpolar0)\Delta G^0 = (\Delta G^0_{\mathrm{vdw}} + \Delta G^0_{\mathrm{ele}}) + (\Delta G^0_{\mathrm{solv, polar}} + \Delta G^0_{\mathrm{solv, nonpolar}})

对应于MM/PB(GB)SA计算结果中呈现的能量项。

3. 入口

平台地址：（https://cloud.yinfotek.com/)

平台左侧菜单栏【计算方案】->【大方案】->【分子动力学】->【结合自由能计算（MM/PB(GB)SA）】

4. 平台操作步骤

4.1 提交任务

  1. 从云盘选择拓扑文件和轨迹文件

    • 拓扑文件必须与轨迹文件原子对应；

    • 轨迹文件可多选，但须注意顺序（除非采样范围为全部）。

  2. 设置采样范围

    点击【计算总帧数】获得当前轨迹总帧数，再按需设置采样间隔。

    最终采样帧数 = (结束帧号 - 开始帧号 + 1) / 采样间隔。例如，(10000 - 1 + 1) / 50 = 200。

    建议对RMSD平稳的那段轨迹进行采样，跨度和间隔也应尽量大。

  3. 划分受体和配体，设置盐浓度

    填写amber mask（详见《【分子动力学模拟（Amber 20）】》中的预备知识），指定受体或配体。

    在本例中，受体为1-266号残基，配体为267号残基。因此，有下图所示两种写法。

    当体系包含更多成分时，对受体、配体的划分需要特别注意。例如，体系包含蛋白、辅酶、有机小分子、Mg离子（与辅酶配位），想计算有机小分子的结合力，可以这样划分：将蛋白、辅酶、Mg离子归为受体，有机小分子归为配体。

  4. 选择盐并设置盐浓度。

    虽然MM/GB(PB)SA采用隐式溶剂模型，并未真正添加离子，但不同的盐对应的离子强度不同。

4.2 修复PDB

  1. 查看结合自由能

    表格给出GB和PB模型的结合自由能数据（详情见预备知识），mmpbsa.xlsx文件包含了结合自由能数据和能量分解数据。

    结合自由能各项能量含义为：

    ΔGvdw\Delta G_{\mathrm{vdw}}：范德华能量项

    ΔGele\Delta G_{\mathrm{ele}}：静电能量项

    ΔGpolar\Delta G_{\mathrm{polar}}：极性溶剂化能

    ΔGnonpolar\Delta G_{\mathrm{nonpolar}}：非极性溶剂化能

    ΔGgas\Delta G_{\mathrm{gas}}：分子力学项（气相中的能量），ΔGgas=ΔGvdw+ΔGele\Delta G_{\mathrm{gas}} = \Delta G_{\mathrm{vdw}} + \Delta G_{\mathrm{ele}}

    ΔGsolv\Delta G_{\mathrm{solv}}：溶剂化能，ΔGsolv=ΔGpolar+ΔGnonpolar\Delta G_{\mathrm{solv}} = \Delta G_{\mathrm{polar}} + \Delta G_{\mathrm{nonpolar}}

    ΔGtotal\Delta G_{\mathrm{total}}：总结合自由能，ΔGtotal=ΔGgas+ΔGsolv\Delta G_{\mathrm{total}} = \Delta G_{\mathrm{gas}} + \Delta G_{\mathrm{solv}}

    GB和PB模型计算的是溶剂化能部分，ΔGgas\Delta G_{\mathrm{gas}}的计算是一样的。结合自由能数值越小，结合力越强。GB模型的ΔGtotal\Delta G_{\mathrm{total}}通常比PB模型的值小很多，这意味着两者在溶剂化能计算部分有较大差异。在应用这些数据解析实验现象或探究分子机制时，应注意与实验数据是否冲突，是否能解析得通，选择更符合实际的。

  2. 查看能量分解数据

    与结合自由能对应，该页也按GB和PB模型分别给出能量分解的数据。

    • 每个残基的总能量等于其骨架和侧链的加和。
    • 分子动力学模拟中残基编号是从1开始重新排序的，可能与实际编号不一致，请下载mmpbsa.xlsx文件，修改成实际编号，再进行分析会更为方便。

点击表头三角箭头按钮，勾选过滤条件，找到贡献突出的关键残基，再结合3D结构分析，就能精确地了解受体-配体间的相互作用了。

原则上，（正或负）贡献突出的关键残基是ΔGtotal\Delta G_{\mathrm{total}}绝对值较大的，但ΔGgas\Delta G_{\mathrm{gas}}更能反映出受体与配体之间相互作用大小，与3D结构分析中看到的相互作用力直接对应起来，而ΔGsolv\Delta G_{\mathrm{solv}}则反映了溶剂化有利于还是不利于结合。

例如，下图分别为GB和PB模型的计算结果。GB结果显示，GLU54的ΔGes\Delta G_{\mathrm{es}}高达-17.85 kcal/mol，导致其ΔGgas\Delta G_{\mathrm{gas}}达到-18.51 kcal/mol。这表明，它与配体分子有强烈的静电力作用，是个关键残基。但ΔGpolar\Delta G_{\mathrm{polar}}同样很大，致使ΔGsolv\Delta G_{\mathrm{solv}}达到19.66 kcal/mol，抵消了静电力作用带来的强结合力贡献（ΔGtotal\Delta G_{\mathrm{total}}为正数）。但这个残基并非是多余的，从3D结构中可以清楚看到，它与配体和邻近的赖氨酸形成强烈的盐桥作用。

再查看PB给出的ΔGsolv\Delta G_{\mathrm{solv}}，尽管仍然为正数，但比GB小很多，ΔGtotal\Delta G_{\mathrm{total}}已为负数（-5.77 kcal/mol）。由此可知，该残基确实是关键残基，PB计算值更符合结构分析结果，应采用PB的数据分析本体系。

实际上，侯廷军等人在2012年学术年会的一篇评估研究报告指出，MM/PBSA的预测结果更接近实验数值，而MM/GBSA则对相对结合自由能预测性能更好。对MM/PB(GB)SA方法计算的结合自由能，应结合实际进行取舍。

另外，需要注意的是，PB模型未实现对非极性溶剂化能（ΔGnonploar\Delta G_{\mathrm{nonploar}}）的分解，因此该项数值总是0，但通常绝对值很小，不影响分析。

参考文献

[1] Amber tutorial: http://ambermd.org/tutorials/advanced/tutorial3/index.php

[2] 侯廷军, 李有勇. MM/PBSA和MM/GBSA对蛋白-配体自由能计算精度的评估研究. 中国化学会学术年会, 2012.

文章作者: 殷赋量子氢
文章链接: http://blog.yinfotek.com/posts/10402/
版权声明: 本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源 殷赋量子氢 !
分子动力学
  目录