多尺度神经网络势函数结合各向异性信息传递用于蛋白质动力学与酶促反应的快速精确模拟

《Journal of the American Chemical Society》：Multiscale Neural Network Potential with Anisotropic Message Passing for the Fast and Accurate Simulation of Protein Dynamics and Enzymatic Reactions

　　

研究人员提出了一种新一代AMP神经网络势函数（NNP），该势函数采用各向异性信息传递（AMP）架构，旨在通过类似量子力学/分子力学（QM/MM）结合静电嵌入的多尺度方法，以密度泛函理论（DFT）精度研究凝聚相中的大型生物分子系统。研究人员在新近发布的生物分子多尺度模拟数据集BMS25上训练了AMPv3模型，证明了该模型的高效率，使得研究人员能够模拟涉及数千个原子的蛋白质系统，并在显式分子力学（MM）溶剂存在下达到DFT精度，模拟时长可达100纳秒，这对于当代NNP而言是一项重大飞跃。研究人员观察到该模型在单GPU上对大系统具有出色的扩展性。AMPv3-BMS25（简称AMP-BMS）在基准测试上表现 promising，研究人员证明该模型可用于准确估算实验性质，包括小分子的溶剂化自由能和蛋白质的结构特征。最后，研究人员采用AMP-BMS/MM预测了分支酸变位酶和氟乙酸脱卤酶催化反应的自由能面。总体而言，AMP-BMS/MM在凝聚相中对蛋白质进行了累计23微秒模拟时间或480亿积分步的模拟。该工作确立了AMP-BMS作为一种高效、准确的生物分子多尺度模拟模型。

计算化学领域长期将蛋白质动力学与酶促反应的精确建模视为核心目标之一。分子动力学（MD）模拟通过原子尺度相互作用与宏观性质的联系为生物过程提供洞见，在药物计算机辅助设计或定制化催化功能与选择性酶的开发中有广泛应用。量子力学（QM）提供了化学的理论基础，催生了M?ller-Plesset微扰理论和Kohn-Sham密度泛函理论（DFT）等电子结构方法。然而，尽管硬件能力与电子结构代码取得重大进展，由于这些方法随系统规模的不利标度，对涉及数十万个原子的系统进行从头算MD模拟目前仍不可行。因此，自20世纪70年代以来，经典力场被开发用于以经典力学描述原子间相互作用来研究此类大型系统。这种原子相互作用的简化描述相比QM方法具有有利的计算标度，但牺牲了精度与可迁移性。

在典型的蛋白质水溶液MD模拟中，水分子间相互作用占系统总相互作用的大部分，导致大量计算资源被分配至通常不太关注的水-水相互作用。因此，早期便提出了量子力学/分子力学（QM/MM）等多尺度方法，将系统划分为QM区和MM区。该形式化方案可在感兴趣区域保留QM方法的高精度，同时利用经典力场对环境区域（MM区）的更好标度。此类方案显著降低了从头算MD的计算成本，但采用DFT处理QM区的QM/MM方案对于溶液中蛋白质的长时程模拟而言仍然过于昂贵。

神经网络势函数（NNP）作为解决精度与计算效率之间权衡的潜在方案应运而生，经过相应参考数据训练后可作为昂贵QM哈密顿量的廉价替代。近年来，在底层机器学习（ML）模型的理论与架构发展方面付出了巨大努力。随着领域成熟，工作重心逐渐转向通过将NNP与现有MD及多尺度模拟软件包集成、为更广泛社区开发用户友好界面来使能实际应用。

涉及NNP的MD模拟分为两类。第一类中，整个系统由NNP描述，类似从头算MD。虽然比DFT便宜得多，此类别中NNP模拟仍受计算成本与内存需求限制，采样时间和大型系统应用受限。需要多GPU设置才能访问更大系统，因此常应用过短的图截断，这在凝聚相系统中可能因长程相互作用的重要性而限制精度。解决这些局限的尝试包括增加长程相互作用项。

第二类则将系统类似QM/MM方案划分为由NNP评估的高精度区域和由经典力场描述的降低精度区域。由于与QM/MM的相似性，此类方案常被称为ML/MM，允许以更低计算成本模拟更大系统。这些设计要取得成功，不同子系统间的耦合至关重要。目前该类别大多数方法依赖机械嵌入方案或近似静电嵌入。尽管机械嵌入计算高效且易于NNP实现，但其因缺乏MM区对QM区的极化而为人诟病。为避免这些局限，研究人员和其他人已开发了直接使用静电嵌入的ML/MM方案。

2023年，研究人员发表了专为ML/MM模拟设计的各向异性信息传递（AMP）架构。AMP架构的第二版（AMPv2）用于传播小肽的MD模拟，并估算包括带电物种和过渡金属在内的挑战性分子体系的解离自由能。该系统仍为每个系统训练AMP模型，并探索了模型对类似分子的可迁移性。受AMPv2 promising性能的鼓舞，研究人员着手开发可应用于模拟蛋白质乃至涉及数万个原子的酶促反应的基础性AMP模型。

本工作中，研究人员介绍了AMP架构的第三版（AMPv3），该版本在约60000个分子系统、150万个构象的短肽和小分子（包括显式溶剂化）的QM/MM参考数据（BMS25数据集）上进行训练，为生物分子凝聚相系统的静电嵌入ML/MM MD模拟提供通用NNP。AMP模型的命名约定采用[AMP架构]-[嵌入方案]-[训练集]格式，因此该模型为AMPv3-EE-BMS25，简称为AMP-BMS。除了蛋白质的100 ns模拟外，本工作还展示了AMP-BMS用于描述水中酶催化反应的能力，这是ML/MM方案的重要应用。

AMPv3模型继承并扩展了先前工作中的SO(3)等变信息传递架构，用于静电嵌入形式下的ML/MM系统描述。本工作中的修改旨在实现更大系统的高效扩展和更高精度。核心创新在于ML与MM区域间的耦合机制：所有蛋白质和配体原子分配至ML子系统，溶剂分子分配至MM子系统。总势能表示为ML/ML、ML/MM和MM/MM相互作用之和。

ML区内相互作用分解为四个组分：AMP原子贡献总和、D4模型的双体色散校正、采用移位反应场方法的长程静电相互作用，以及Ziegler-Biersack-Littmark（ZBL）短程排斥势。AMP模型在每个原子上构建k阶原子多极矩M_ik，用于在后续信息传递步骤中纳入方向信息。通过张量积形成局部基矢，多极矩由边特征和隐藏节点特征预测。每个信息传递步中，多极矩用于构建各向异性特征g_ij，进而更新隐藏特征。相比初版AMP，多项简化降低了计算成本并提高了精度，最重要变化涉及多极张量的维度：不再使用单一多极矩组，而是在每个原子上展开多组多极矩，以提高各向异性信息传递方案的分辨率。

长程色散相互作用采用SpookyNet提出的处理方式，即D4色散模型的双体组分。由于准确描述分子偶极子对ML与MM粒子间相互作用至关重要，C₆系数不从用于静电相互作用的相同部分电荷获取，而是由神经网络参数化函数预测的标度因子与元素的真空C₆系数参考值相乘得到。长程静电相互作用通过移位反应场方法纳入，使用单调态假设，认为短程高阶贡献由ML势捕获。

ML与MM区域间的相互作用是AMP架构的决定性特征之一，能够高效描述周围电荷分布引入的极化效应。ML/MM相互作用分为静电相互作用和Lennard-Jones（LJ）项。为考虑MM部分电荷对ML区的极化，在信息传递最后一步展开额外多极矩组：极化贡献描述为极化率、MM粒子引起的电场、距离相关贡献和阻尼因子的乘积。所得多极矩与信息传递期间展开的其他多极矩拼接以引入ML与MM区的耦合。静电相互作用描述为ML区各向异性多极矩与MM区部分电荷之间的相互作用，包括单极-单极、单极-偶极和单极-四极相互作用。LJ相互作用参数取自经典力场，用于计入ML与MM区之间的交换排斥和色散相互作用。

研究人员在BMS25数据集上训练了AMPv3架构，该数据集包含约60000个独特分子系统的150万个静电嵌入QM/MM参考计算，涵盖短肽（最多三个残基）、单肽二聚体、小分子和反应路径快照，理论水平为ωB97M-D4/ma-def2-TZVPP，溶剂为TIP4P-FB水。所得AMP-BMS模型随后用于执行复杂度递增的多个ML/MM MD模拟系统。

静电嵌入校准方面，QM/MM计算中经典力场参数与QM哈密顿量的兼容性仍是长期挑战，这也延伸至ML/MM设置。此问题源于经典力场的受限函数形式及其独立于与之组合的QM方法进行的参数化。非键力场项通常优化以复现液体热力学性质如密度和气化焓，而非匹配QM计算的能量、几何和其他微观性质。因此，力场参数是有效参数，可能补偿受限函数形式中忽略的和互作用项。例如，TIP3P电荷过度极化以补偿显式极化建模的缺失。这种过度极化的MM部分电荷反过来导致QM区域电子密度的过度极化，引入非物理的"溢出"效应，即因缺乏Pauli排斥导致QM电子密度错误定域到MM原子上；或者高估QM与MM区域间的静电相互作用能量。

对于AMP-BMS/MM方法，研究人员研究了两种解决潜在过度极化及由此导致的ML/ML与ML/MM相互作用失衡的策略。第一种策略称为极化标度，通过参数λ_pol∈[0,1]阻尼ML区中由MM区诱导的多极矩，在极限情况λ_pol = 0时恢复机械嵌入。第二种策略为电荷标度，遵循QM/MM模拟中的成熟实践，MM电荷通过λ_ch∈[0,1]标度，阻尼诱导和永久多极矩的幅度。两种方案仅在推理时应用于ML与MM粒子间的相互作用，涉及单一参数，与距离、电荷大小或极化率无关。

两种策略参数均通过25个有机分子的水合自由能（HFE）进行校准，这些分子因与氨基酸侧链的化学相似性而被选用。HFE量化溶质从气相转移到水溶液的热力学成本，广泛用于经典力场的验证和参数化。电荷标度的最优性能在λ_ch = 0.9时发现，相对于实验值的平均绝对误差（MAE）为3.0 ± 0.7 kJ·mol^-1。极化标度的最优值在λ_pol = 0.35处，MAE为2.6 ± 0.5 kJ·mol^-1。值得注意的是，使用各自最优系数时，两种标度方案均优于FreeSolv数据库中报告的GAFF计算值（MAE为5.3 kJ·mol^-1）。这些发现进一步由ML水分子氧原子与周围MM水分子氧原子间的径向分布函数（RDF）、对称适应微扰理论（SAPT）二聚体相互作用能，以及色氨酸笼蛋白的NOE衍生残基间距离证实。

两种标度方案总体结果非常相似，因此主文本仅报告极化标度λ_pol = 0.35的结果。对于未来应用，研究人员推荐采用极化标度方案，因其在需要时允许更简单的系统特定校准。总之，这些发现表明无阻尼情况下溶质-溶剂ML/MM相互作用被显著高估，凸显了多尺度方案中ML/MM相互作用校准的必要性。

基准数据集验证方面，高精度势能面的准确表示是开发预测性NNP的必要条件。研究人员首先在BMS25数据集的寡肽和迷你蛋白外部测试集上验证AMP-BMS，确立其准确复现参考方法（ωB97M-D4/ma-def2-TZVPP）能量和梯度的能力；其次在真空中的若干文献基准上评估其相对于半经验、DFT和CCSD(T)水平QM计算的精度。

凝聚相基准中，NNP训练数据集很少包含扩展化学结构，因参考QM计算随系统规模增加的标度不利。NNP必须从训练时呈现的较小系统泛化到推理时的较大系统。为研究这种泛化能力，在短肽、小分子和反应路径上训练的AMP-BMS模型在BMS25数据集的寡肽和迷你蛋白子集上进行评估。这些子集旨在揭示NNP中可能存在的尺寸或电荷偏差，包括不同大小和电荷（4-12个残基，-2至+2电荷）的寡肽随机序列和迷你蛋白（166-487个原子，-2至+8电荷），这些分子显著大于训练集中最大分子（84个原子）。

寡肽结果表明，ML区梯度的MAE在不同分子尺寸间保持恒定，但随肽总电荷略有增加。梯度分量中位偏差在低个位数kJ·mol^-1 ?^-1范围内，但观察到少数大异常值，推测可能源于硬自由度的小偏差。每原子相对势能的MAE随系统规模增加而略有下降，表明误差增速慢于总能量增加，这是可取的性质。迷你蛋白子集结果与寡肽相当，未观察到系统规模与误差大小之间的明显相关性。总之，应用AMP-BMS/MM于日益增大的系统时未检测到误差显著增加，这对电荷和系统规模均成立。

对AMPv3-EE-BMS25模型静电嵌入质量的评估表明，模型以合理精度复现了短程相互作用的诱导偶极矩方向和大小（中位角偏差14°），而长程性能因涉及的诱导矩较小而如预期较弱。观察到诱导偶极矩大小的系统性低估（MAE = 0.035 e?），与分量解析嵌入能分析一致——该分析揭示λ_pol = 1.0时准确的总嵌入能源于高估的静电相互作用项E₁₂（平均符号误差MSE = +24.6 kJ mol^-1）与低估的极化代价E₃（MSE = -24.0 kJ mol^-1）之间的近乎完美抵消。将λ_pol降至0.35抑制了这种抵消，产生有效近似总嵌入能为E₁₂单独的模型，相对于DFT参考存在系统性低估。这一系统性低估导致与实验观测的更好一致性，构成直接证据表明DFT嵌入能本身因非极化MM水模型的过大固定电荷而过高。

真空基准中，AMP-BMS模型在凝聚相数据（即静电嵌入QM/MM）上训练。为评估其在无溶剂存在时的性能，模型在不同文献气相基准上进行评估。假设AMP架构的模块化特性——仅在信息传递最后一步考虑外场——使向真空环境的泛化变得直接。

首先考虑BMS25数据集中所有氨基酸及其质子化变体，基准测试AMP-BMS对抗参考DFT方法。为获得当前基准的气相氨基酸构象，从原始BMS25研究中的氨基酸水溶液SD轨迹采样50帧。尽管未在气相数据上训练，AMP-BMS对氨基酸达到一致低误差，平均相对构象能量MAE为1.93 kJ·mol^-1，平均各力分量MAE为2.07 kJ·mol^-1 ?^-1。

接着在文献气相数据集上基准测试AMP-BMS的相对构象能量，并与三种额外方法比较：BMS25 DFT方法、基于广泛OMol25数据集训练的NNP（MACE-OMol-0），以及基线半经验方法GFN2-xTB。在GMTKN55子集上，AMP-BMS对大部分数据集（ACONF、Amino20x4、BUT14DIOL、MCONF、PCONF21、SCONF）的MAE值低于化学精度阈值（4.184 kJ·mol^-1），仅UPU23略超。37CONF8的MAE恰在阈值。对于葡萄糖和α-麦芽糖构象数据集，以及Wiggle150数据集，AMP-BMS产生约7 kJ·mol^-1的较高MAE值，其中葡萄糖和Wiggle150的BMS25 DFT参考方法也显示出与耦合簇参考的较大偏差。总体而言，AMP-BMS性能仅略逊于MACE-OMol-0，尽管后者NNP特征130倍更多参数（总计5200万可训练参数）并在OMol25数据集的1亿个气相计算上训练。这些结果令人鼓舞，表明未来通过在AMP训练集中纳入气相数据来改进AMP的直接途径。

计算复杂度评估方面，对不同生物分子系统评估了AMP-BMS模型的计算成本。在单NVIDIA H200 GPU（141 GB VRAM）上应用于涉及多达十万个ML原子和超过一千万个MM原子的系统时，观察到看似亚线性标度（k = 0.92）。所用大系统包括SARS-CoV-2刺突蛋白δ变体（PDB: 7VHH）、炭疽毒素受体（PDB: 3HVD）和完整卫星烟草花叶病毒衣壳（PDB: 4NIA）。

有利的标度表明现代加速器提供的大规模并行性和专用硬件组件随系统规模增加仅缓慢饱和。CPU（AMD Ryzen 9 7950X，196 GB系统RAM）上观察到线性对数??(N log N)标度（k = 1.15），与理论考虑一致。基于GPU获得的标度，除以ML原子数后的模拟速度为每ML原子每步5.1 × 10^-6秒（仅ML原子）和2.7 × 10^-8秒（含MM原子）。这对应于1000原子系统、0.5 fs时间步长约8.5 ns/天的模拟速度。模拟约需每1000 ML原子1.4 GB VRAM（MM原子对VRAM影响可忽略）。

溶液中蛋白质动力学应用方面，先前涉及NNP的蛋白质MD模拟范围相当有限，通常仅持续数纳秒。研究人员在三种随机初速重复下对四种具有多样二级结构基序组成的蛋白质进行了AMP-BMS/MM MD模拟：鸡卵白溶菌酶（HEWL，PDB: 1AKI，1960 ML原子）、主要冷休克蛋白（CspA，PDB: 1MJC，1009 ML原子）、蛋白G（GP，PDB: 1PGB，855 ML原子）和来自结核分枝杆菌的分支酸变位酶（CM，PDB: 2FP2，2564 ML原子）。

三个系统（HEWL、GP和CM）的骨架RMSD在约20纳秒后稳定至1-1.5 ?范围，与经典力场观察到的RMSD值相似。CspA的三个重复显示相对较大变化，发生在连接β-折叠片的柔性环残基中。DSSP分析验证结构特征。HEWL的NOE质子-质子距离和¹H-¹⁵N序参数（S²））与实验数据比较，AMP-BMS/MM产生的结果与实验测量紧密匹配，总NOE违规率6%，S²的MAE为0.045。这些结果证实QM参考数据验证中观察到的低误差转化为数十纳秒的稳定MD模拟，能够复现蛋白质结构特征。

肽和迷你蛋白熔解曲线应用方面，蛋白质熔解处于展开和折叠结构之间，对分子模拟方法构成特殊挑战。对于这一现象，相互作用平衡的微小扰动可导致二级结构稳定性的显著变化，仅关注稳定折叠时不易检测。先前对蛋白质熔解的研究揭示了经典力场在复现熔解曲线方面的局限，表现出过度稳定某些二级结构基序的趋势。

为研究AMP-BMS/MM描述稳定折叠之外蛋白质的能力，研究人员在278至348 K温度范围内对三组重复的四条肽和迷你蛋白进行了模拟：(AAQAA)₃肽、 chignolin（PDB: 1UAO）、GB1发夹（PDB: 1PGB）和villin头结构域（PDB: 1VII）。这些分子先前用于评估力场CHARMM36IDPSFF。由于每温度有限模拟时间（3 × 100 ns，前50 ns作为平衡），模拟从折叠结构开始（展开通常比折叠更快观察到）。与经典力场的可比研究相比，此模拟时间更短，因此结果可能未完全收敛，部分数据点的三重重复间存在较大变化。尽管如此，清晰趋势已显现：所有系统均观察到实验与计算熔解行为间的良好一致性。相比之下，经典力场难以复现这些熔解曲线。有趣的是，熔解曲线是观察到极化标度和电荷标度方案差异的一个应用场景，电荷标度的实验一致性较低。

这些结果凸显了在QM数据上训练的NNP最显著的优势之一：无需针对蛋白质和系综性质的特定训练，AMP-BMS/MM即能复现折叠蛋白质动力学和蛋白质熔解特征，表明该方法提供了比经典力场对这些系统更好的物理表征，以及ML/MM方案中溶质-溶质和溶质-溶剂相互作用的准确平衡。

酶促反应应用方面，化学反应研究是NNP重要且有吸引力的应用领域，因其承诺以大幅降低的计算成本达到与QM参考方法相同的精度。研究人员和其他人已证明NNP可用于研究溶液中的化学反应，但这些研究聚焦于小系统尺寸，常需针对所研究反应的特定训练数据。研究人员探究AMP-BMS/MM是否可直接用于酶促反应自由能面的定性和定量建模。

在已建立的酶QM/MM MD设置中，QM区仅包涵反应位点及选定距离截止内的周围蛋白质残基。这种系统划分需要使用连接原子或其他近似处理跨越QM/MM边界的化学键截断，引入额外截断伪影。为减少截断伪影，已在QM区纳入数千原子并采用静电嵌入方案被证明对获得收敛能量是必需的。鉴于AMP-BMS/MM的架构和出色标度，研究人员的方法为两种局限提供了有吸引力的解决方案，因为酶的完整原子可纳入ML区，仅溶剂原子保留在MM区，以静电嵌入方式纳入。

研究人员研究了两种酶催化反应。第一个是分支酸转化为预苯酸，由分支酸变位酶催化，是莽草酸途径中芳香族氨基酸合成的关键步骤之一。此反应是QM/MM研究最广泛的酶促反应之一，因为底物在反应期间不与酶形成共价键，意味着完整酶可指派给MM区，无需处理断裂化学键。其次，研究人员研究了氟乙酸脱卤酶催化降解氟乙酸的决速步。

研究人员采用伞形采样计算了水中（24 ML原子，28000 MM粒子）和酶催化（5687 ML原子，85000 MM粒子）的分支酸转化为预苯酸的自由能面。与实验数据和先前研究一致，观察到放能反应（水中ΔG = -55.2 kJ·mol^-1，酶环境-95.7 kJ·mol^-1）。水中计算自由能差与实验数据（ΔG = -56.0 kJ·mol^-1）一致。虽然AMP-BMS未在该反应上训练，底物沿反应坐标的结构变化如预期所示。酶促反应屏障高度显著降低（50.1 kJ·mol^-1），与实验数据和先前报告匹配，解释了约10⁶的速率加速。然而，水至酶的ΔG^?相对降低虽匹配实验，绝对值与实验测定值仍有相当偏差。研究人员将此偏差追溯到训练数据对过渡态结构的有限覆盖。通过微调测试此假设，验证集表明默认AMP-BMS模型在过渡态区域相比反应物和产物区存在显著更大误差。微调至参考理论水平后，势垒降低47.2 kJ·mol^-1，但仍高估实验值。进一步研究表明，势垒高度随密度泛函选择变化约65 kJ·mol^-1，最关键因素为精确交换分数。更高精确交换分数导致更预测更高势垒，这与DFT近似中固有的自相互作用误差紧密相关。精确交换分数较低的泛函可能因幸运误差抵消而与实验更好一致。微调至B3LYP-D4/6-31G(d)水平后，活化势垒进一步降低，与实验值偏差在1 kJ·mol^-1以内，确认了部分偏差源自参考泛函本身而非ML模型的缺陷。

氟乙酸脱卤酶中，提出的酶促去氟化机制涉及至少三步和催化三联体（Asp110亲核试剂、His280碱、Asp134羧酸酯）。研究人员聚焦此酶催化水解的决速步——水分子对碳酰基的亲核攻击。通过Asp134Ala突变评估模型预测突变结构势垒高度的能力。采用AMP-BMS/MM对野生型酶（9257 ML原子）及其Asp134Ala突变体（9253 ML原子）进行了伞形采样，每组>230000 MM粒子。野生型酶预测势垒ΔG^?为113 kJ·mol^-1，Asp134Ala突变体增加至143 kJ·mol^-1，相对增加对应约10⁵倍，与类似机制酶催化皂化反应中观察到的相对催化活性降低一致。由于Asp134距离反应水和碳酰基较远（约8 ?），结果表明AMP-BMS/MM能够描述此类较长程效应。与分支酸异构化发现相似，绝对势垒高度相比实验值78.4 kJ·mol^-1仍过高，预期可通过扩展BMS25数据集和/或微调AMP-BMS模型来校正。

研究结论表明，默认AMP-BMS模型已可直接用于ML/MM模拟，无需微调即可研究酶催化反应中的相对趋势，这是NNP开发和应用中的重要里程碑。

AMP-BMS训练细节包括：两步行信息传递，图/短程截断4 ?，ML/ML长程截断14 ?，ML/MM极化截断8 ?，ML/MM LJ和静电截断9 ?，MM/MM相互作用截断9 ?，节点尺寸128，边尺寸32，ML/ML贝塞尔函数8个，ML/MM贝塞尔函数2个，多极通道8个，多极阶数2，总计412074个可训练参数。模型用PyTorch/2.2.1实现，采用He初始化、全连接神经网络和Swish非线性，Adam优化器训练24个epoch，指数衰减学习率。损失函数包括能量、ML和MM区梯度、分子偶极矩和四极矩、以及C₆系数的加权和。