gmx_MMPBSA的计算-part6
gmx_MMPBSA的计算
前言：这一part将在闲时写完。下一part将介绍线性相互作用能(LIE)方法，我感觉这个方法一直没有受到重视，但是基于已有分子进行训练来说，我觉得它可能也很靠谱。
背景
gmx_MMPBSA是基于Amber中的MMPBSA.py进行的结合自由能计算，因此需要事先安装AmberTool20或者更高版本。AmberTool20可由官方网站免费下载，大家可根据官方手册说明自行安装。它是基于AMBER MMPBSA.py用于计算终态自由能(gmx格式文件)的工具。包括两个程序：gmx_MMPBSA：主程序，执行计算功能。gmx_MMPBSA_ana：图形界面分析数据并保存高质量图片。
目前已知的支持力场的文件是：AMBER和CHARMM，所以在使用gromacs运算时，大家一定要选好力场。主推AMBER，GROMOS力场大家可以放弃了。
分子力场泊松玻尔兹曼表面积（Poisson−Boltzmann surface area，PBSA）和广义波恩表面积（Generalized Born surface area，GBSA）是非常受欢迎的计算结合自由能的两个模型。它们计算精确度高于对接函数，计算时间少于炼金术结合自由能（Alchemical binding free energy ABFE）计算。MM/PBSA和MM/GBSA被广泛运用生物大分子的折叠、蛋白配体相互作用、蛋白-蛋白相互作用。
热力学上的自由能指的是体系内部所具有的能量，决定了热动力学过程的走向以及体系所处状态的概率。研究自由能是研究生物分子运动和相互作用的手段之一，实验测试自由能非常耗时且经费投入较多，因此利用动力学模拟获取分子体系的自由能变化是目前主流的计算方案，尤其是在药物设计、蛋白质折叠和蛋白-蛋白相互作用研究领域。在药物与蛋白结合的热力学过程中，在等温等压条件下体系的自由能通过如下公式定义： F= -kBT lnZ。 其中kB为玻尔兹曼常数，T为体系温度，Z为体系分布函数，如果假定体系处于热动力学稳定状态，则Z根据如下公式定义：
其中V0是体积常数，N是离子数，h是普朗克常数，P和V是压力和体积。H(p,r)为指定动能（p）与位置（r）的粒子的汉密尔顿函数。有上述公式可知当原子数较少时，可以直接计算其绝对自由能，然而当原子数增高大时，计算绝对自由能变得非常困难。因此研究多原子体系的时候更多的是通过计算相对自由能变化来研究热力学过程的走向。将研究状态A与B之间的能量变化通过如下公式定义：
其中PA和PB为状态A与B的概率。
基于以上理论开发了很多计算药物与蛋白结合自由能的计算机工具，对接函数速度快精度低，而更为复杂的炼金术结合自由能计算，例如自由能微扰（FEP）和热力学积分算法（TI）精度高但计算非常耗时，而且由于生物分子体系的复杂度高，计算往往无法收敛，实际操作过程中另一个缺陷是准备步骤繁琐，容易出错。虽然近来有利用GPU加速来提升其计算速度，但是仍然不适合大规模筛选，只在具有较高计算资源项目组进行个别分子的优化研究中得以运用。MM/PBSA和MM/GBSA方法则介于上述两者之间，兼具精度速度，是广泛被使用的计算结合自由能的工具。而且这些工具还能进行能量分解，每个氨基酸残基对能量的贡献值对于指导药物设计也是极具帮助意义。
MM/PB(GB)SA在药物设计中的应用
1. 虚拟筛选
为了快速计算，通常采用简单的打分函数来估计结合亲和力，但是在大多数对接打分函数中简化或者忽略了溶剂化自由能。在许多情况下，将分子对接和MM/PB(GB)SA结合是一种合理的策略，可正确识别结合构象以及正确排序。有科研同仁同时对98个靶标的“对接能力”用MM/PB(GB)SA进行评估，发现在大多数情况下可以显著提高找到配体的正确结合构象的比例（成功率为69.4%）。但是预测的结果受到原子电荷、溶质介电常数、MD模拟长度以及熵计算等等的影响。
2. 蛋白-蛋白相互作用
蛋白-蛋白相互作用（PPIs）在活细胞的大多数生化过程中发挥关键作用。通过计算方法探索两个蛋白间的相互作用，通常在初始阶段使用精度低的打分函数快速筛选大分子数据库，随后使用更严格但更耗时的打分函数进行重新打分排序。研究表明，MM/GBSA比ZDOCK能更好地预测两个蛋白结合的构象。因此，考虑到低计算成本和相对较高的预测精度MM/GBSA可能是预测蛋白-蛋白系统正确结合构象和结合亲和力的强大工具。
3. 蛋白-核酸体系
蛋白-核酸（DNA/RNA）相互作用在许多生物过程起到重要作用，如基因表达的调控、RNA剪切等。但是现有的蛋白-RNA相互作用（PRIs）打分函数并不可靠，因此准确预测PRIs的3D结构和结合亲和力仍然相当困难。MM/PB(GB)SA可被用于PRIs的研究。Orr等人使用MM/GBSA预测RNA修饰是否有利于与目标蛋白结合，预测结果与实验数据有很高的相关性(r2> 0.9)。目前感觉不是很纯熟，因此这方面的合作，一概不接。
问题集绵来啦，这都是重点：
[1] gmx_MMPBSA计算时帧数多少合适？
在使用 gmx_MMPBSA 计算时，帧数的多少会影响计算结果的准确性和稳定性。一般来说，帧数越多，结果越准确，但计算时间也会变长。因此，需要在时间和准确性之间取得一个平衡。通常建议使用至少1000个均匀分布的快照进行计算。这意味着从分子动力学模拟中选取至少1000个时间点，这些时间点应该覆盖整个模拟过程。如果模拟时间非常短，则可以使用更多的帧来覆盖整个模拟过程。另外，如果您希望更加准确的结果，可以使用更多的帧来进行计算。如果计算时间不是主要问题，建议使用至少10,000个帧进行计算。如果需要更高的准确性，则可以使用更多的帧，但要记住，使用过多的帧可能会导致计算过于耗时。
需要注意的是，使用过少的帧可能会导致结果不准确，而使用过多的帧可能会导致计算时间过长。因此，需要根据具体情况来确定合适的帧数。
回答来自GPTchat，回答很中肯，采纳。
后记：如何确定自己选了1000帧？
答：一般我们步长是0.002ps，输出精度我常选1000。这样就是2ps保留一帧，1000ps=500帧，即1ns=500帧。我通常选后期10ns进行计算，即5000帧。在gmx_MMPBSA计算时选择间隔5帧，即计算了1000帧。
[2] 关于多线程自己电脑实测：
普通MMGBSA和MMPBSA计算可以多线程。
IE无法多线程。
24核超线程，最多用12，选多报错。
[3]GB和PB选择哪个？
Ⅰ：在应用这些数据解析实验现象或探究分子机制时，应注意与实验数据是否冲突，是否能解析得通，选择更符合实际的。
Ⅱ： 侯廷军等人在2012年学术年会的一篇评估研究报告指出，MM/PBSA的预测结果更接近实验数值，而MM/GBSA则对相对结合自由能预测性能更好。对MM/PB(GB)SA方法计算的结合自由能，应结合实际进行取舍。侯廷军, 李有勇. MM/PBSA和MM/GBSA对蛋白-配体自由能计算精度的评估研究. 中国化学会学术年会, 2012.
Ⅲ：MM/PB(GB)SA性能比较总结，选哪个自己决定：
（1）MM/PB(GB)SA预测精度严重依赖于力场、电荷模型、连续溶剂模型、内部介电常数、采样方法以及构象熵。
（2）力场在分子模拟中起到关键作用，决定了系统的所有相互作用。由于目前隐式溶剂模型都是非极化力场发展起来的，而MM/PBSA更适合极化力场。测试结果表明在MD模拟和MM/PB(GB)SA计算中使用相同的力场，不同力场的混合导致预测不准确。
（3）结果表明PB和GB方法预测能力受到半径的影响。
（4）GB模型具有最好的预测结合亲和力的能力。
MM/PBSA方法比MM/GBSA方法表现出更好的排序性能。
RESP电荷对MM/PB(GB)SA方法具有最佳性能，而基于AM1-BCC和ESP电荷预测的结果也相当令人满意。
（5）在计算绝对结合自由能的性能上，MM/GBSA预测能力比MM/PBSA差，但是MM/GBSA在不含金属的体系中具有更好的结合亲和力的排序能力。
（6）个别体系中MMGBSA的精度已经达到了FEP的计算精度，表明MMGBSA在计算结合自由能上的优势，相比之下MMPBSA的计算性能对体系敏感度更高，而且实际操作中MMPBSA的计算更加耗时，由于溶剂模型的影响MMPBSA经常出现结合自由能为正数的异常情况，也是导致MMPBSA使用不如MMGBSA频繁的原因之一（网上常问问题，原因就是这个，你算正无外乎选择错分组或者换MMGBSA）。
（7）MD时长不是越长越好，对于MMGBSA计算1ns的动力学模拟很多时候也可以满足计算需求。所以说还是有时候。
（8）PB的原理是最严格的，所以大量GB模型的文章都用PB的结果当参照判断GB的设置、求解有效Born半径的方法的合理性。-来自卢天。
参考文献：[1] Wang E, Sun H, Wang J, et al. End-point binding free energy calculation with MM/PBSA and MM/GBSA: strategies and applications in drug design. Chemical Reviews, 2019, 119(16): 9478-9508.
