登上Nature子刊！深势科技RiD方法解决生物大分子体系高维空间采样难题
近日，深势科技团队题为"Efficient Sampling of High-dimensional Free Energy Landscapes Using Adaptive Reinforced Dynamics"的工作被 Nature Computational Science 接收。
强化动力学（Reinforced Dynamics, RiD）方法借助神经网络的高维表示能力，可以处理几十甚至上百维的集合变量（Collective Variables, CVs），大大降低了CV选取的敏感性，使得增强采样方法能够更高效的应用在复杂的生物大分子体系中。本工作在强化动力学的基础上提出了自适应置信阈值的方法，有效地探索了上百个集合变量和具有高能垒体系的自由能面。
分子动力学模拟的
效率难题
分子动力学（Molecular Dynamics, MD）模拟是通过计算机模拟分子体系物理运动的数值计算方法。分子动力学理论在经典牛顿力学框架下（使用经典力场近似），通过哈密顿力学方程演化体系的经典轨线。这一理论假设体系满足各态遍历性，即通过无限长时间的模拟计算，我们可以通过MD算法遍历体系相空间中的任意状态。现实中我们不可能做到无限长时间的计算，往往以有限时间的模拟得到体系的近似性质。通过分子动力学模拟，我们可以得到体系的热力学性质的统计估计量；通过分析得到的轨迹，我们可以得到系统在有限时间下的动力学行为。
分子动力学模拟具有广泛的应用，如蛋白质折叠，蛋白质结构优化，自由能计算等；在材料体系，如相变理论，液体理论等也有大量研究。在CASP竞赛（Critical Assessment of Protein Structure Prediction ）中，基于分子动力学方法的蛋白结构优化策略取得了优异的结果[1]。
但分子动力学方法存在两大难题：力场精度和采样效率。目前所使用的经典力场是通过拟合量子化学计算和实验数据，得到的势能函数形式和参数，如CHARMM，AMBER等。力场参数对模拟的结果具有较大的影响，由于MD本身的经典近似和力场参数精度的影响，MD模拟的结果无法达到量子化学精度。在更微观的体系，或者涉及化学反应（化学键断裂）的体系，经典分子动力学模拟应用有限。另一方面，为了节省计算资源，我们往往只能进行有限时间的模拟，轨迹很难遍历体系的各个状态。对于含有较高能垒的复杂体系，经典分子动力学模拟效率较低，轨迹容易陷入势能面上的一个局部极小值势阱，难以观察到体系的明显变化。
对于两大问题，我们正在致力于发展具有量子化学精度的神经网络力场[2-4]和发展对复杂体系的增强采样方法[5]。
针对采样效率问题，很多增强采样方法（Enhanced Sampling）被提出，包括基于温度的方法和基于偏置势的方法。基于温度的方法，如温度副本交换动力学（Replica Exchange Molecular Dynamics, REMD），在模拟过程中交换高温下的轨迹与低温下的轨迹，通过升高温度加快体系越过自由能势垒的过程。基于偏置势的方法，如元动力学（Metadynamics），在探索过的反应坐标（也叫集合变量，Collective Variables，CVs）上施加高斯函数形式的偏置势，迫使体系向未探索的区域演化。这里的反应坐标是指一组关于体系直角坐标系的函数，它往往描述了体系演化过程中较慢的自由度。加快这些慢自由度能够大大加快体系的演化速度。在这两类方法中，能选定反应坐标的体系，基于偏置势的方法一般更高效。
基于偏置势的方法在复杂生物体系表现欠佳，这是因为无法在模拟前很好定义能够准确描述反应方向的一个集合变量。而传统方法同时处理多个（3个以上）的反应坐标时，模拟会遇到“维度灾难”，模拟效率大大降低。
强化动力学方法
应对复杂蛋白质体系
针对这一问题，作者开发了强化动力学(Reinforced Dynamics，RiD)方法。强化动力学借助神经网络能够处理几十甚至上百维的集合变量，大大降低了CV选取的敏感性，使得增强采样方法能够应用在复杂蛋白体系中。
RiD的工作流由三个部分组成：采样（Exploration）（图1 a）、平均力计算（Labeling）(图1 b)、网络训练（Training）(图1 c)。
图1 RiD的工作流
其中，RiD将使用模型偏差（即多个模型的输出的标准差）来指示当前自由能估计的优劣并指导偏置势的添加——模型偏差小，该区域学习到的自由能较为准确，那么在模拟时我们将在该区域添加偏置势；模型偏差大，表明数据集不充分，学习到的自由能不准确，则不会在该区域施加偏置势。这里便涉及到界定模型偏差最大可接受值的一个超参数，我们称之为置信阈值（Trust Level）。
这里，作者使用自适应调节的置信阈值——对采样结果进行层次聚类，通过聚类数自适应调整置信阈值，避免模拟中长时间陷入局部势阱。
RiD的“精雕细琢”
助力蛋白质折叠与优化
蛋白质折叠
蛋白质的折叠问题可以作为验证分子模拟方法效率的良好基准。我们使用强化动力学方法模拟了蛋白质chignolin的折叠过程。
使用蛋白质的18个二面角作为CV，进行31次迭代计算，结果如图2所示。图2 a,b 为模拟中两条轨迹相对于折叠态的RMSD，这表明RiD模拟过程既观察到了折叠过程（RMSD<0.15 nm），也观察到了去折叠过程。轨迹中最优构象相对于折叠态结构的RMSD为0.019 nm。图2 d-f展示出了RiD构建的高维自由能面的投影图，我们将这一结果与经典MD模拟的数据相对比（图2 g-i），结果表明RiD方法得到的自由能面与经典MD的结果吻合，同时可以采样到更广的自由能面。
图2 蛋白质chignolin的折叠与去折叠
同时作者也对比了元动力学框架下的其他方法:偏置交换元动力学（Bias-exchange Metadynamics, BEMetaD）和平行偏置元动力学（Parallel Bias Metadynamics, PBMetaD），表1为蛋白质的折叠与去折叠速率。总体而言，强化动力学模拟和元动力学框架下的不同方法都成功折叠了蛋白质chignolin，但只有RiD和PBMetaD可以较好的去折叠，这意味着RiD和PBMetaD面对较深的亚稳态时，可以很好的进行采样。
表1 蛋白质chignolin的折叠与去折叠速率（number/μs）
注：BE0.2, BE0.5, BE0.8为BEMetaD, PB0.5为PBMetaD，具体参数细节见原文
这一测试展示了RiD方法在采样效率和构建自由能面的能力。
蛋白质结构优化
过去数年间，通过引入多序列比对得到的共进化统计信息，和最新的深度学习方法，人们对无模板蛋白质结构预测取得了重大突破，标志性事件是2020年 AlphaFold2 在结构预测精度上的巨大提升。AlphaFold2 对 CASP14 比赛中有47% 蛋白的结构预测的结果 Global Distance Test–high Accuracy （GDT-HA）[6] 打分小于75，对这些蛋白质仍然需要进一步的结构优化。
目前主要的蛋白质优化（Refinement）手段中有两类表现优异：一类是通过图神经网络（Graph Neural Network，GNN）对蛋白编码后生成距离约束，进一步通过Rosetta等程序进行优化 [7]；另一类是从初始结构出发，借助分子动力学模拟（Molecular Dynamics，MD）对局域相空间采样，从所得构象系综中找出自由能面上的能量极小点作为优化结构。其中，基于MD方法的蛋白质优化结果表现最好[1]。基于经典MD的蛋白优化方法仍存在较大的局限性。其主要问题在于：难以对相空间充分采样找到自由能全局极小点。如果蛋白的两种构象之间存在较大的势垒，那么经典MD过程需要经过很长时间才能越过这一势垒，这导致了对蛋白质的结构优化有限。
在这里，作者将RiD方法应用于蛋白质结构优化情景中。作者使用了CASP13比赛中的三个例子（R0974s1、R0986s1、R1002-D2）对RiD进行了测试（图3）。对于三个测试用例，作者分别选取了136，182，116个二面角作为CV，采用了8条轨迹并行采样，使用了含有4层隐藏层、每层含有1200个神经元节点的深度神经网络（Deep Neural Network, DNN）进行学习，经过16个迭代后，得到精修后的结果。
图3 RiD在蛋白结构精修中的结果（以CASP13中的R0974s1, R0986s1 和 R1002-D2为例）
使用GDT-HA对结果进行打分，满分为100，分数越高表明结果越接近晶体结构。三个结构的初始GDT-HA分数为65.6，59.2和72.9，经过RiD精修后，它们的分数均达到78以上，其中最好分数分别为92.4，82.6和79.7（图3），结构优化提升显著，远好于目前分数最高的FEIGLAB（经典分子动力学方法）。通过计算三个例子的标准差可以看出（图3b），RiD的模拟结果标准差最小，方法具有较好的鲁棒性。尤其是对于R1002-D2一例（图3h），经典MD轨迹从初始结构（图3h 黑色X）明显陷入了局域势阱中（图3h 紫红色区域）而根本无法触及晶体结构（蓝色X），RiD方法则成功越过势垒并采样到了晶体结构附近。这充分体现了RiD在蛋白体系和高维空间中的采样能力。
图4 R1002-D2体系在RiD模拟中的构象变化，白色为晶体结构
具体而言，如图4d中的红色箭头所示，R1002-D2体系预测出来的结构（图4d，蓝色结构）N端β-strand中存在约0.57 nm的移位误差（图4d，红色箭头）。在RiD模拟过程中，β-strand之间的氢键首先断裂，并导致RMSD超过0.4 nm的构象变化（图4 a, e蓝绿色结构）。只有在这种大的构象变化之后，才能看到正确的β-strand，RMSD约为0.11 nm（图4f，红色结构）。因此，RiD的成功归因于沿着过渡路径跨越高自由能垒的能力。这一测试也表明RiD的高采样效率，具有潜力能够完成抗体优化，Loop区优化等任务。
RiD方法展望
为领域突破带来可能
蛋白质在生物体内是不断运动的，它们通过不断地构象变化，以及与其他分子的结合、解离行使生物学功能，而目前描述蛋白质构象系综的最佳方式之一是借助分子模拟手段。RiD作为生物大分子体系的高效采样问题的出色解决方案，为多个领域的突破带来了可能。
对于药物设计而言，RiD能够在更多的场景发挥重要作用：探寻构象变化中的隐藏口袋和别构口袋，困难药物靶点中固有无序蛋白的构象采样，抗体环区优化等。我们非常欢迎对别构和固有无序蛋白等场景感兴趣的伙伴一起进行更多的探索。
目前，我们已经将文中的RiD方法重构为开源软件供用户使用，可以在下方链接获取更详细的使用资料。
地址：
https://github.com/deepmodeling/rid-kit
我们鼓励读者直接看该工作原文，如有文章获取困难，或想与开发者进一步交流，请您关注深势科技公众号并回复关键词“DP-RiD”，我们期待您的交流反馈。
Reference
[1] Heo, L. & Feig, M. Experimental accuracy in protein structure refinement via molecular dynamics simulations. Proc. Natl Acad. Sci. USA 115, 13276–13281 (2018)
[2] Zhang, L., et al. "Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics." Physical review letters 120.14, 143001 (2018).
[3] Zhang, L., et al. "End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems." Advances in Neural Information Processing Systems, 4441-4451 (2018).
[4] Zhang, Y., et al. "DP-GEN: A concurrent learning platform for the generation of reliable deep learning based potential energy models." Computer Physics Communications 253 , 107206 (2020).
[5] Zhang, L., Wang, H. & E, W. Reinforced dynamics for enhanced sampling in large atomic and molecular systems. J. Chem. Phys. 148, 124113 (2018).
[6] Zemla, A. LGA: a method for finding 3D similarities in protein structures. Nucleic Acids Res. 31, 3370–3374 (2003)
[7] Jing, X., Xu, J., Fast and Effective Protein Model Refinement Using Deep Graph Neural Networks, Nature Computational Science. 1, 462–469 (2021).
上下滑动查看更多
关于深势科技
深势科技有限公司（“深势科技”）是一家成立于2018年的科技公司，于2019年开启规模化发展之路，致力于以新一代分子模拟技术解决微观尺度工业设计难题。 以打造切实服务于药企、材料商和科研机构的模拟研发平台为主要业务方向，以解放研发工作者的生产力为主要业务目标。
深势科技具有强大的科研与产业落地能力。其新一代分子模拟算法在保持量子力学精度的基础上，将分子动力学的计算速度提升了至少五个数量级，且对算力的需求与体系的原子数量呈线性依赖；结合高性能计算，能够对数十亿原子规模的体系进行量子力学精度的计算模拟。团队核心成员获得2020年全球计算机高性能计算领域的最高奖项“戈登·贝尔奖”，相关工作当选2020年中国十大科技进展，以及2020年全球人工智能十大科技进展。