DRUGAI
今天为大家介绍的是来自John D. Chodera团队的一篇论文。开发可靠且可扩展的分子力学（MM）力场——这些快速的用于描述分子系统势能面的经验模型——对于生物分子模拟和计算辅助药物设计是不可或缺的。在此，作者介绍了一种通用且可扩展的机器学习MM力场——espaloma-0.3，以及一个使用图神经网络的端到端可微分框架以克服传统基于规则的方法的限制。espaloma-0.3在单个GPU上训练一天即可拟合一个包含超过110万次能量和力计算的多样化量子化学数据集，能再现与药物发现高度相关的化学领域（包括小分子、肽和核酸）的量子化学能量特性。此外，该力场保持了小分子的量子化学能量最小化几何结构，并保留了肽和折叠蛋白的凝聚相特性，自洽地参数化蛋白质和配体，生成稳定的模拟结果，从而高度准确地预测结合自由能。此方法展示了系统构建更准确且易于扩展到新化学领域的力场的显著潜力。
分子力学（MM）力场是一种快速的经验模型，通过将生物分子系统视为一组原子质量点来描述其势能面。这些质量点通过非键相互作用和价键相互作用（键、角和扭转），这些相互作用通常被参数化以再现量子化学构象能量学和物理特性。尽管其对基础物理模型的表示进行了简化，MM力场已被证明在生物分子模拟和计算辅助药物设计中是不可或缺的，适用于多种任务，如列举假定的生物活性构象、通过虚拟筛选进行命中识别、预测膜渗透性、生物分子动力学模拟以及通过炼金术自由能计算估算蛋白质-配体结合自由能。
I类MM力场是速度和精度之间广泛流行的折衷方法
I类MM力场因其简单的函数形式提供的计算效率，被广泛应用于蛋白质、脂质、核酸及其他相关生物分子。
分子系统的总势能UMM由坐标x定义，通过为系统中的每个原子或价键项（键、角、扭转）指定的力场参数ΦFF = {Kr, Kθ, r0, θ0, Kø,n, ø0, q, σ, ε}i来确定。可以引入平面外项（不正当扭转）与扭转项一起使用，以改善分子的平面性。范德华相互作用通常使用Lennard-Jones 12–6势和Lorentz–Berthelot结合规则来描述，不同原子类型之间的参数σ和ε由此确定，但也可能采用其他组合规则确定。实际上，这种相互作用通常需要结合独立开发的特定化学领域的不同力场参数，以补充生物分子系统的异质性。
传统MM力场参数化方法在难以应对复杂性，限制精度
传统上，构建MM力场需要物理有机化学方面的专家知识，以建立原子类型规则，将原子分类为表示不同化学环境的离散类别，从而能够从相关原子、键、角和扭转参数表中分配MM参数。这就创建了一个难以处理的混合离散-连续优化问题，构成了一个依赖人力的劳动密集型挑战。力场的准确性受化学感知分辨率的限制，而化学感知分辨率又受限于不同原子类型的数量。试图通过增加原子类型数量来提高准确性会导致键、角和扭转参数的组合爆炸，带来了很大的实际限制。因此，建模者经常求助于定制的参数生成工具，如Paramfit、FFBuilder或OpenFF BespokeFit，为需要高精度的特定分子分配参数，这需要进行昂贵的量子化学计算，从而降低了I类力场的速度优势。
传统MM力场参数化方法采用分而治之，而不是自洽性方法
为了控制原子类型复杂性的爆炸，生物分子力场通常采用为蛋白质、小分子和其他生物分子独立构建独立但据称兼容的模型。例如，最近发布的AmberTools 23建议结合独立开发的力场来模拟包含蛋白质、DNA、RNA、水、单价和双价反离子、脂质、碳水化合物、糖缀合物、小分子、翻译后修饰和核酸修饰的系统，这些力场的开发可能总共花费了超过100人年的努力。虽然这简化了对每类分子进行参数化的子问题，但将这些独立的力场一起用于处理复杂的异质系统既不简单也不理想。每个力场设计用于建模的化学空间通常存在重叠，无法保证这些区域的参数完全一致且兼容。这在多类生物分子相互作用时引入了显著的问题，存在准确性差的风险，并在不同类别的分子必须共价键合的情况下大大增加了难度。因此，扩展到新类别的生物分子或化学空间变得非常耗时，因为结合力场通常会导致可能的力场参数的巨大组合空间，其结果的力场质量很大程度上取决于用户的选择。
已经有许多努力试图系统化和自动化力场开发过程。例如，Open Force Field Initiative开发了一些现代的、开源的工具、数据集和力场，它们采用了一种直接的化学感知方法，使用标准的基于SMARTS的化学子结构查询以分层方式分配整个价键参数集（原子、键、角、扭转），试图缓解参数的组合爆炸。然而，大部分工作集中在小分子上，将力场扩展到新的化学领域仍然需要人为努力——共同优化离散的化学感知规则和连续的力场参数仍然是不可行的。
Espaloma 提供了一个灵活的、端到端可微分的框架，用于使用GNN分配分子力学（MM）参数
图1：Espaloma在三个连续的阶段中参数化分子系统
如图1所示，Espaloma类似于基于原子类型的力场操作，在其中化学感知被预定义以生成MM力场参数ΦFF。然而，Espaloma不是处理原子类型，而是使用图神经网络在化学图G上操作，该网络由神经网络模型参数ΦNN参数化：
生成的参数ΦFF随后可以在标准的分子力学软件包中使用，以计算任何构象的MM能量和力。
阶段1
图神经网络生成连续的向量化原子嵌入，取代离散的原子类型规则。首先，使用诸如RDKit之类的化学信息学工具包，将分子系统抽象为一个图，其中节点和边分别表示原子和共价键。Espaloma使用GNN替代基于规则的化学环境感知，在这个图上操作。这些神经架构通过消息传递迭代更新并池化嵌入向量到相邻原子，从简单的输入特征中学习有用的原子化学环境表示。因此，化学图中的对称性（化学等价性）在本质上得以保留，同时得到了丰富、连续且可微的原子环境表示。
阶段2
对称性保持池化生成连续的键、角和扭转嵌入。第一阶段中由GNN确定的表示被用来以保持对称性的方式生成键、角和扭转的表示，其中相关的等效原子排列通过Janossy池化列出并求和。
阶段3
原子、键、角和扭转的神经网络参数化取代了表格参数分配。在最后阶段，前馈神经网络学习从这些对称保持的不变原子、键、角和扭转嵌入到MM参数ΦFF的映射，这些参数反映了适用于这些项的特定化学环境。每个不同的参数类别（如原子、键、角和扭转参数）由单独的神经网络分配，使这个阶段完全模块化。这个阶段类似于传统力场构建中的最终表查找步骤，但由于捕捉特定于分配势能项的化学环境的连续嵌入，它提供了显著的灵活性。
最终输出的是由神经网络基于其相关权重ΦNN唯一确定的一组力场参数ΦFF。这意味着一旦优化了ΦNN，生物分子模拟可以像使用传统MM力场那样快速进行。在Espaloma框架内，还可以使用与几何无关的电荷平衡方法快速生成AM1-BCC质量的原子部分电荷。
详尽的开放量子化学数据集涵盖生物分子：小分子、蛋白质和核酸
为了开发一个广泛适用于生物分子建模的自洽MM力场，作者首先整理了一个高质量的气相量子化学数据集，并将其存储在QCArchive中。这个整理后的量子化学数据集由多个组件构建，为相关的生物化学提供了互补的覆盖范围。
为了捕捉生物分子的粗糙构象能量面，量子化学数据集从三种不同的QCArchive工作流程中提取：Dataset、OptimizationDataset和TorsionDriveDataset。Dataset包含结构的单点能量计算，这些结构不一定在其局部量子能量极小值处，由MD模拟或构象生成器生成。OptimizationDataset是给定结构的QM优化轨迹集合。TorsionDriveDataset涉及对一组可旋转的扭转进行扫描，然后进行QM优化。
Espaloma力场再现量子化学能量和力
作者增强了原始的Espaloma框架，以提高模型稳定性和数据效率：
- 在训练中加入了量子化学力，以提供更多关于量子化学势能面的信息；
- 对正扭转和不正扭转力常数应用L2正则化，以抑制扭转曲线中的虚假特征；
- 使用n = 1, 2周期性表达不正扭转项，以降低模型复杂性，并与通常采用n = 1, 2周期性的其他传统力场对齐；
- 消除对共振形式敏感的节点特征，以确保同一分子的化学等价表示接收相同的参数。
表1：再现量子化学能量与力方面性能比较
如表1所示，espaloma-0.3在再现量子化学能量和力方面显著优于所有基线力场（gaff-2.11, openff-2.0.0, openff-2.1.0, Amber ff14SB, Amber RNA. OL3），展示了espaloma-0.3比当前一代I类MM势能更准确地再现生物分子和有机化学的量子化学能量面。相反，广泛应用于生物分子模拟领域的基线力场相对于量子化学计算产生了相当大的能量误差和巨大的力误差（平均是espaloma-0.3的两到三倍）。这种性能优势在各种化学类别中都得到了验证，表明espaloma-0.3在广泛的化学和生物化学建模任务中的普遍适用性。
Espaloma力场保持量子化学能量极小值
作者接下来研究了espaloma-0.3定量再现量子化学平衡构象能量的能力是否能够延伸到定性保持量子化学局部能量极小值构象的能力——这对于准确表示诸如配体结合对接研究、模拟或自由能计算等现象的几何结构至关重要。为评估这一点，作者使用标准化的气相QM优化几何结构基准来比较使用espaloma-0.3和基线力场（openff-2.0.0、openff-2.1.0和gaff-2.11）优化的构象与其在B3LYP-D3BJ/DZVP理论水平下的QM优化几何结构之间的差异。
图2：Espaloma-0.3保留了量子化学能最小值的位置
如图2a和b所示，与量子化学参考值相比，使用espaloma-0.3优化的构象几何结构和相对能量比基线力场显示出更好的一致性。此外，使用espaloma-0.3获得的MM优化几何结构中的键、角和扭转与量子化学值显示出接近的一致性（图2c），总体性能与基线力场相当或略好。espaloma-0.3的键异常值（>0.1 Å）来自三个连接到脂肪族碳的磺酰胺，共包含30个构象，占整个基准数据集构象的0.04%，其S–N键距离在磺酰胺基团中比QM优化几何结构长约0.4 Å。
Espaloma力场再现肽和折叠蛋白的实验核磁共振观测
肽
为了定量评估espaloma-0.3对氨基酸内在主链偏好的建模能力，作者对十三种短的无序肽进行了MD模拟，这些肽的核磁共振（NMR）观测已通过实验测量。这些肽由3到5个残基组成，无保护基团，并由于NMR实验的低pH环境而具有质子化的C端。测量的邻近标量耦合揭示了这些肽的主链二面角偏好。使用Karplus模型从500 ns轨迹中计算标量耦合，并使用χ²值量化与实验可观测量的一致性。
图3：非结构化肽的实验核磁共振标量耦合性能比较
总体而言，espaloma-0.3与实验的吻合度比ff14SB更高，如低χ²值所示（图3a）。值得注意的是，ff14SB在具有短侧链的氨基酸（如甘氨酸和丙氨酸）上的实验一致性更高（图3b）。这并不令人惊讶，因为ff14SB的主链扭转参数被调整以再现数据集中包含的丙氨酸5聚肽的NMR标量耦合。然而，espaloma-0.3在更具挑战性的氨基酸（如带电的赖氨酸、体积大的蛋氨酸或β支链的缬氨酸）上与实验的一致性更高。
折叠蛋白质
图4：折叠球状蛋白的实验NMR标量耦合比较
通过对四种球状蛋白质（GB3、BPTI、溶菌酶和泛素）进行MD模拟，评估了折叠蛋白质中主链和侧链c1的内在动力学（图4）。这些蛋白质的标量耦合已通过NMR实验测量，并且在蛋白质力场开发中已被广泛研究。结果显示，espaloma-0.3准确再现了所有四种折叠蛋白质的实验NMR标量耦合，但其ANE值略高于ff14sb，这较大的偏差往往来自侧链标量耦合。
Espaloma力场准确描述了蛋白质-配体结合自由能
图5：蛋白质-配体结合自由能计算
在图5和表2中，作者展示了espaloma-0.3在蛋白质和配体自洽参数化的情况下，其蛋白质-配体结合自由能性能与ff14SB openff-2.1.0相当。espaloma-0.3在绝对（ΔG）和相对（ΔΔG）自由能上的RMSE分别为1.02 [95% CI: 0.74, 1.37] kcal mol−1和1.12 [95% CI: 0.88, 1.41] kcal mol−1。相应地，ff14SB openff-2.1.0的ΔG和ΔΔG RMSE分别为1.01 [95% CI: 0.73, 1.33] kcal mol−1和1.21 [95% CI: 0.93, 1.54] kcal mol−1。尽管报告的误差和相关统计数据的置信区间有重叠，但这些结果令人鼓舞，因为espaloma-0.3展示了其覆盖不同化学领域的能力，而传统力场在这方面已经奋斗了几十年却未能实现。
值得注意的是，如图5所示，所有三种情况下Mcl1系统均出现了一个较大的异常值。使用ff14SB espaloma-0.3计算的相对结合亲和力ΔΔG为4.05 kcal mol−1（图5b）。然而，当从一个反转的结合构象进行炼金结合自由能计算时，作者发现误差可以减少到2.60 kcal mol−1，这与实验差值（ΔΔG = −0.54 kcal mol−1）更为一致。
讨论
espaloma-0.3通过自洽参数化蛋白质和配体，展示了其在再现量子化学能量、实验NMR数据和蛋白质-配体结合自由能方面的出色表现。尽管存在一些异常情况，espaloma-0.3依然证明了其在生物分子建模和药物发现中的潜力和广泛适用性。本文还强调了espaloma-0.3在准确性、灵活性和扩展性方面的显著提升，为未来力场开发提供了重要参考。
编译 | 于洲
审稿 | 曾全晨
参考资料
Takaba K, Friedman A, Cavender C, et al. Machine-learned molecular mechanics force fields from large-scale quantum chemical data[J]. Chemical Science, 2024.
