NC 2026 | KIMMDY: 一种生物分子反应仿真平台
今天介绍的是一项发表在 Nature Communications 上的工作，展示了 KIMMDY（KInetic Monte Carlo Molecular DYnamics）平台的创新应用。KIMMDY 是一种能够模拟大规模生物分子反应的仿真工具，旨在解决传统分子动力学模拟在反应性生物分子系统中的局限性。传统的分子力学（MM）模型在处理生物分子反应时，通常无法模拟反应物之间的化学变化，尤其是在长时间尺度上的反应过程。KIMMDY 通过结合分子动力学（MD）和动力学蒙特卡罗（kMC）方法，采用图神经网络（GNN）来预测反应率，并在此基础上选择反应并进行模拟，从而有效克服了这一限制。
KIMMDY 可以应对多种反应类型，包括自由基反应、亲核取代反应以及光二聚化等。通过与实验数据的比较，KIMMDY 展示了其在蛋白质和核酸系统中的强大应用潜力，能够揭示复杂的反应路径，帮助我们重新解读实验数据，并为未来的实验提供新的思路。此外，KIMMDY 还能够模拟反应级联和竞争反应，为研究人员提供了前所未有的计算效率和灵活性。
获取详情及资源:
0 摘要
分子模拟在生物学研究中已成为不可或缺的工具。其准确性持续提高，但直接模拟所有生命过程中的生化反应仍然是计算上具有挑战性的。提出了一种生物分子反应仿真器，利用动力学蒙特卡洛方法对构象集合中的反应进行建模。方法KIMMDY能够处理动态、大规模系统中的连续竞争反应，即便是在秒级或更慢的时间尺度上。通过图神经网络预测反应速率，同时也能够使用更简单的基于物理的或启发式模型。方法已通过实验数据验证，并通过一系列应用展示了其强大功能和多样性，包括自由基反应、亲核取代反应和光二聚化反应。示例系统包括蛋白质和DNA。KIMMDY有助于理解复杂系统中的生化反应级联过程，帮助重新解释实验数据，并为未来的湿实验提供灵感。
1 引言
生物大分子本质上具有反应性。它们的反应性驱动了生命过程中必不可少的生物学过程。大自然利用大量巧妙的生化反应，从亲核取代和加成反应到重排和氧化还原反应，推动着代谢、信号传递和能量转移等生命过程中至关重要的环节。除这些严格控制且特定的反应外，生物大分子在其生命周期内会经历广泛的非特异性反应，这些反应受到变化条件（如pH、光照或氧化应激）的影响，在生理条件下会导致持续的分子损伤。
这些生化过程发生在并受到细胞内外复杂、拥挤和动态分子环境的驱动。分子模拟为深入理解细胞内部机制提供了前所未有的视角，并帮助克服了实验中时空分辨率的限制。生物分子模拟中的原子级或粗粒度模型的长度尺度和分子复杂性正接近甚至已达到细胞级别。然而，这些模拟仍然缺乏生命的关键特征——反应性，通常为了提高效率而牺牲了这一点。分子力学（MM）力场是处理大规模系统并在足够的时间尺度上进行模拟的首选工具，但这些模型无法模拟化学键的断裂或形成。然而，监测分子运动与化学反应之间的相互作用是理解和设计复杂分子系统的关键。
为了克服经典模拟的局限性，提出了反应力场、允许特定反应发生的方法，如经验价键方法和常规pH模拟。此外，近年来的机器学习原子间势能（MLIPs）正在重新定义原子级建模的角色，这些模型天生具备反应性。
所有这些模型——MLIPs、经典反应性模型，甚至高效的非反应性力场——都面临一个根本性问题，即其可达到的时间尺度远低于大多数生化反应的时间尺度。混合方法，如量子力学/分子力学（QM/MM）和机器学习/分子力学（ML/MM），也存在相同的限制，它们的性能受限于较慢的计算方法。即便这些混合方法在未来变得更高效，它们在处理生物大分子的复杂化学反应时仍然不可行，特别是在面对翻译后修饰或自由基链反应等反应时，反应区域可能遍布整个系统，且反应的时间尺度可达到秒级。克服这些限制将为对反应性生物分子系统的大规模模拟提供可能。
图1｜KIMMDY与QM/MM方案的对比 a. 在QM/MM（或类似的MLIP、MLIP/MM方案）中模拟反应（上图）与使用KIMMDY模拟反应（下图）进行比较。b. 与QM（或MLIP）方法进行过渡态搜索的高昂费用不同，KIMMDY通过从MD系统中提取反应物构型，并通过kMC步骤模拟反应，利用反应特定模型推断反应速率。
KIMMDY（动力学蒙特卡洛分子动态）是一种生物分子反应仿真器。KIMMDY在动态分子系统中搜索反应，并基于仿真器（即机器学习、启发式或物理替代模型）预测反应速率。随后，KIMMDY在动力学蒙特卡洛（kMC）步骤中选择一个反应，实施该反应并生成新的构象集合，以便进一步仿真反应。尽管早前已提出过kMC和MD的组合，但方法首次实现了化学动力学（即反应速率）与反应物构象动态的直接依赖。KIMMDY能够预测高反应性和大型生物分子系统中的反应动态，应对其中的反应级联和竞争反应，并且理论上能够处理几乎任何时间尺度上的反应。
图2｜KIMMDY作为生物分子反应仿真器 a. 中心工作流程的可视化。KIMMDY在动态分子系统中搜索反应，基于学习到的启发式（例如基于实验）或物理模型预测反应速率，通过kMC选择其中一种反应并执行该反应。b. KIMMDY的反应搜索步骤的详细方案，展示了可能的输入集生成器和反应类型。MD和SD分别指分子动力学和随机动力学。
尽管此前已使用过MD和kMC的组合，但方法直接考虑了构象集合，并利用图神经网络（GNN）的优势快速预测基于集合的反应势垒。同时，也利用了经典生物分子物理力场的已验证精度和速度，并通过机器学习简便地扩展至新的化学体系。KIMMDY经过实验验证，展示了其在开壳和闭壳反应以及蛋白质和核酸系统中的应用。KIMMDY被实现为一个模块化平台，可以扩展到任意反应，只要反应速率和力场可通过机器学习模型或实验数据获取。为充分展示这一扩展性，包括了从广泛反应和系统中获得的结果。
2 结果
2.1 模拟的反应动态能够正确预测小分子中的自由基反应
第一个应用领域选择了氢原子转移（HAT）。HAT在生物大分子和软物质中普遍存在，能够发生非特异性反应，并且通常会发生在多个后续自由基迁移步骤的级联反应中，且反应时间尺度通常为秒级，因此这是一个标准方法（如QM/MM或直接MD模拟结合机器学习势能）难以处理的反应类型。首先验证了KIMMDY在小型烃类自由基中的HAT反应，该反应的实验数据已有报道。所使用的系统包括从1-丙基自由基到1-辛基自由基的n-烷烃衍生自由基。用于推断kMC反应速率的机器学习势能障碍预测模型是通过在DFT水平上计算蛋白质HAT反应的势垒进行训练的，假设模型足够通用，能够推广到烃类。需要强调的是，在此处的速率是根据仅从反应物结构预测的反应势垒推断的，而不是通过进行昂贵的过渡态搜索（使用QM或MLIP）。
图3b展示了典型的KIMMDY工作流程，以一个辛基自由基为例。图3c展示了不同HAT反应在500K下的选择概率，并与实验结果进行了比较。选择了实验得出的速率的平均值作为给定氢转移反应的参考值，再除以所有考虑的氢转移反应的平均实验速率之和。同样，从KIMMDY模拟中提取了这些氢转移的选择概率。与实验结果一致，方法正确识别了1–5作为主要的HAT反应，这是由于过渡态的六元环配置较为有利，同时准确预测了1–n氢转移反应概率的整体趋势，特别是对于高环张力过渡态的转变的低概率。每次运行中，随着势垒预测数量的增加，氢转移速率很快收敛（图3d）。因此，KIMMDY稳健地再现了整体反应动态，且对于最主要的氢转移反应，其相对误差较小。
HAT反应的整体正确排序和稳健性在与绝对反应速率的直接比较中也得到了验证，并且选择概率在图3c中的20次独立KIMMDY模拟运行中保持了大致稳定。然而，基于GNN的仿真器低估了绝对反应速率。其原因可能在于GNN仿真器的训练数据，其中包括了受限的过渡态，且仅优化了直接参与的原子。此外，量子力学隧穿效应未被考虑，而隧穿效应已知能在300K下将反应速率提高两到三个数量级。将零点振动能量校正纳入仿真器的训练数据集，并对模型进行进一步调优，以适应烷基自由基，可能会进一步改善结果。尽管存在低估的情况，KIMMDY仍然准确再现了预期的反应概率，这是在反应仿真器应用场景中的关键量。同时，验证结果表明，KIMMDY与学习的反应势垒能够便捷地应用于未在训练数据集中明确包含的系统（如n-烷烃），证明了其广泛适用性。
图3｜验证KIMMDY：n-烷基自由基中的氢原子转移（HAT） a. 在1-辛基自由基中进行的1–4、1–5和1–6 HAT的示意图。未配对的电子用红点表示。b. KIMMDY HAT预测周期。KIMMDY在MD生成的构象集内搜索可能的HAT反应，通过仿真器推断其势垒，选择一个反应并通过kMC执行，调整拓扑结构，开始新的MD模拟。c. 使用KIMMDY推导的七烷基和辛烷基自由基的HAT反应概率，并与n-烷基的实验概率进行比较。对于给定的转移，实验速率取自所有可用的n-烷基自由基并计算平均值。圆圈表示单次100ns的KIMMDY运行（20次独立运行），彩色条表示平均值。在(c)中概率设置为0。原子间的氢转移无法在实验中检测到，因此被排除在外。d. 对每个反应进行不同数量的势垒预测的平均速率（4个辛烷基运行，10ns，500K）。
2.2 生物大分子反应仿真器KIMMDY揭示了蛋白质自由基的反应路径
在验证了反应仿真器在HAT反应中的应用后，接下来应用于生物系统中的第一个例子。自由基物种是由同裂解生成的，随后通过多个连续的HAT反应在系统中迁移（图4a）。胶原蛋白，由对齐并交联的三重螺旋构成（图4b），在受到机械应力时会产生同裂解产生的自由基。基于热力学推导和电子顺磁共振（EPR）数据，提出这些机械自由基可能迁移至DOPA，一种含有儿茶酚的翻译后修饰氨基酸。尽管同裂解位点和DOPA作为自由基清除剂的作用已被广泛研究，KIMMDY提供了能够仔细分析自由基从同裂解位点到自由基清除剂的转移路径的能力。
在一个含有260万个原子的胶原蛋白纤维系统中，观察到在600个模拟的HAT反应中自由基转移路径的种类繁多。通过HAT和构象变化，自由基从同裂解原子迁移了几十埃（图4c）。在模拟过程中，多种氨基酸承担了自由基物种，并且侧链自由基比骨架自由基更为频繁。HAT反应主要是分子间的，或者在氢供体与氢受体之间存在超过六个键，仅包含少数1–2、1–3或1–4的转移反应，这与烷烃的模拟结果一致。对于作为自由基清除剂的DOPA，与不涉及DOPA的反应相比，包含DOPA羟基的反应的反应势垒显著较低（图4d）。中位势垒差为4.5 kcal/mol，使用Eyring方程，这意味着反应速度提高了1800倍。还从同裂解位点到DOPA羟基的直接HAT反应进行了采样。因此，DOPA自由基物种在受力的胶原蛋白中是动力学上可及的。并未观察到DOPA羟基自由基作为氢受体的反应性低于其他蛋白质自由基，并将其归因于多个相邻自由基清除剂之间的快速HAT反应。
有趣的是，来自交联断裂（Cα和Cβ原子，三价吡啶交联PYD）的自由基大部分在反应性模拟过程中未发生迁移（图4c）。表明与在随后的反应中形成的更具反应性的氢受体相比，这些自由基相对稳定。在多个假设检验场景中，测试了具有显著降低的反应势垒的氨基酸氢供体。断裂的PYD的羟基被识别为一个常见的氢供体。此前未曾研究过这种基团作为潜在自由基清除剂，但根据KIMMDY，它的位置理想且动力学上可及。
通过比较PYD的氢抽取键解离能（BDE）与20种氨基酸及DOPA的BDE，探讨了PYD是否能在胶原蛋白中热力学地作为自由基汇。结果显示，PYD和DOPA的BDE最低，这意味着它们是非常优秀的自由基清除剂。先前对拉伸胶原蛋白进行的EPR吸收光谱测量也支持了PYD作为合理的自由基候选者。该光谱表明，DOPA是主要的自由基物种，但也展示了一个无法解释的区域，该区域与计算得到的PYD光谱相吻合）。包括PYD在内的模拟结果提高了与实验数据的整体拟合度，将R2值从0.89提升至0.98。因此，KIMMDY使得能够重新解读实验数据，并识别出胶原蛋白中一个被忽视的自由基稳定基团。
图4｜DOPA和PYD清除同裂解生成的自由基 a. KIMMDY中的自由基连续反应从同裂解开始，随后进行多次HAT。未配对的电子用红点表示。b. 胶原纤维的结构，背骨原子（银色）、潜在的DOPA位点（浅蓝色）和PYD交联（深红色）。放大图显示交联位点的HAT反应，黄色箭头表示。c. 36次模拟中，自由基原子到同裂解的PYD Cα和Cβ原子的距离，分别有八个自由基原子，最大进行20次连续的HAT反应。密度在第20次反应前显示。d. 反应中涉及DOPA羟基作为H供体（n=49，中位数=28.7kcal/mol）、H受体（n=41，中位数=31.4kcal/mol）或没有DOPA（n=467，中位数=33.2kcal/mol）的HAT势垒分布。盒子表示四分位数范围（IQR，Q1–Q3），中心线表示中位数，胡须延伸到1.5×IQR范围内的最极端数据点。使用单侧Welch独立样本t检验评估差异的显著性，未进行多重比较调整（*p<0.05,**p<0.01,***p<0.001）。H供体DOPA羟基与无DOPA之间的差异显著（使用Welch t检验，p=3.9e-7），而H受体DOPA羟基与无DOPA之间的差异不显著（p=1.0）。e. PYD（蓝线）和DOPA（橙线）在不同氨基酸的抽取过程中的键解离能（BDE）分布（来自参考文献27，灰条）。f. 大鼠尾胶原的实验EPR谱（黑色），经参考文献25许可复制，与DOPA的谱（橙色）进行比较，DOPA被提议为信号的关键贡献者。PYD交联的谱（蓝色），作为KIMMDY提出的候选者，可以解释实验数据中未解决的部分。
2.3 KIMMDY能够模拟生物相关系统中的竞争反应
在大多数生化系统中，不同的反应可以相互竞争，导致产生完全不同的产物。为了展示KIMMDY处理不同反应类型的能力，比较了胶原蛋白中的同裂解反应和异裂解反应，类似的模拟也可以用于其他生物分子系统。虽然同裂解会导致后续的链式反应并生成可检测的自由基，异裂解则是一个闭壳反应，在蛋白质中涉及水的攻击，即水解反应（图5a）。这两种反应都可以通过力的作用促进，且都已通过实验观察到，二者的竞争取决于环境因素，如溶剂可接近性或pH值。
通过结合量子力学和分子力学的模拟（QM/MM）可以对这两种反应进行直接比较，但根本上受到采样和计算成本的限制。这正是KIMMDY的模块化系统发挥作用的地方，它提供了一个仿真器，用于在现实的生物环境中模拟它们的竞争。对于同裂解反应，反应速率是根据基于物理的贝尔模型进行模拟的，而水解反应的速率则通过基于力夹具实验、pH值和表面可接近性的启发式模型进行模拟，这展示了即使没有依赖于构象的学习速率，KIMMDY也能轻松使用。
如预期的那样，在没有外部拉力的情况下，水解和同裂解的反应速率都较低，也就是说，在常温常压条件下，肽链的自发裂解不太可能发生。当蛋白质系统受到1 nN的拉力时，水解反应在单一肽链中优先于同裂解反应。然而，在密集且交联的胶原蛋白系统中，在其他条件相同的情况下，同裂解反应的速率显著加速，达到与水解反应相同量级的速率，并在比较最高速率时超越了水解反应。这是由于大胶原纤维系统中高应力集中区域的存在，将反应速率推向了力反应曲线的一个区域，在该区域水解反应速率已经下降，而同裂解反应仍在上升。这一结果对水解模型选择（基于物理的量子力学模型与基于实验的启发式模型）以及其他模型参数具有鲁棒性。
结果表明，在简单分子中反应速率差异较大的反应，在更复杂的分子系统中可能会变得具有竞争性，KIMMDY使得能够直接得出潜在反应结果的结论，例如在胶原蛋白等密集蛋白质网络中的自由基形成。
图5｜水解与同裂解 a. 示意图，展示肽链在碱催化水解或同裂解下的反应。b. 胶原纤维和单一肽链的分子渲染图。c. 两个竞争反应的反应速率比较，基于同裂解的物理模型31和水解的启发式模型28。每个设置的最高速率值显示在相应的分布上方。水解速率分布中的圆圈表示单一肽链的实验参考值。
2.4 KIMMDY揭示DNA折纸图案中意外较低的量子产率
KIMMDY不仅可以应用于蛋白质系统。在紫外线照射下，DNA中的嘧啶碱基可以形成环丁烷嘧啶二聚体（CPDs），这种二聚体在皮肤癌的发展中起着至关重要的作用。在DNA纳米技术中，CPDs被用作连接件之间的共价交联，增强DNA折纸结构的化学稳定性。分子动力学（MD）研究表明，DNA构象中有利于二聚化的频率与实验中观察到的量子产率ϕ之间存在直接关系，而KIMMDY现在提供了对这一情景的直接验证。
使用基于距离和角度的启发式模型，并调整参数以再现小型基准系统的实验量子产率，通过KIMMDY计算反应速率。当模拟框架中的速率超过定义的阈值时，将其计为能够二聚化的构象，从而贡献到量子产率。
在胸腺嘧啶二核苷酸（TdT）系统中，KIMMDY预测了两种产物，即预期的顺式同型（cis-syn）和反式同型（trans-syn）CPD异构体，尽管后者的量子产率要低得多（图6a）。对底层构象的检查显示，在反应发生之前，胸腺嘧啶碱基可以通过绕N-糖苷键的旋转进行syn-anti转变。在实验中，顺式同型异构体的形成频率是反式同型异构体的八倍。KIMMDY结果提出，较低的量子产率以及syn-anti前体状态较少被占据，可能是这一现象的解释。
KIMMDY预测了突出和交叉结构的量子产率（图6b），这些结构的量子产率显著低于双链（ds）或缺口图案（图6d），尽管突出和交叉常用于通过CPD交联来稳定DNA折纸。此差异可能有三个原因：用于折纸稳定化的长时间照射可能弥补了本身较低的量子产率；距离和二面角分布显示（图6c）交叉系统的灵活性增加，这在完整的折纸结构中可能更受限制；当前KIMMDY插件未考虑的其他光产物（如cis-anti、trans-anti和6-4光产物）可能会形成。因此，建议未来的实验研究应调查DNA折纸中量子产率和光产物的多样性。
在具有连续二聚化位点的系统中（图6b）观察到量子产率的预测没有变化（图6e）。表明CPD与第二个反应位点之间相隔五个碱基对的情况不会显著影响可用于额外二聚化的构象空间。表明，先前由于CPD形成而受损的DNA并不一定更容易发生额外的二聚化反应，可能对突变性损伤的形成产生影响。
图6｜使用KIMMDY研究的DNA系统 a. 从反式-反式构象或顺式-反式构象出发的胸腺嘧啶二聚化反应方案，生成不同预测量子产率ϕ的顺式-反式和反式-反式立体异构体。b. 测试的DNA基序的结构示意图。红色圆圈表示胸腺嘧啶。胸腺嘧啶之间的实线连接表示二聚体；胸腺嘧啶之间的虚线连接表示KIMMDY考虑的反应位点。c. 从三个MD模拟中随机采样的10,000个快照中确定的ds中心和交叉系统的距离和角度分布的核密度估计（KDE）。KDE绘制的是反应性胸腺嘧啶双键之间的距离和它们之间的二面角。d. DNA折纸中常见基序的预测量子产率（n=3，ds中心、缺口和悬垂链，n=6，交叉来自不同的100ns MD模拟）。条形图表示平均值±标准差；点表示单个模拟。e. 连续二聚化系统（及相应的控制系统）的预测量子产率，其中两个反应位点相距五个碱基对（n=3，左仅、左先、右仅、右先；n=4，右第二；n=5，左第二，来自不同的100ns MD模拟）。条形图表示平均值±标准差；点表示单个模拟。
讨论
KIMMDY能够模拟单一反应、反应级联和分子系统中的竞争反应，计算成本比其他方法低几个数量级，无论是在从头计算、混合QM/MM、反应力场还是MLIP水平上。与后者及其变体的直接模拟相比，对于大多数生物学相关过程，计算速度慢到几乎无法使用，因为它们目前的时间尺度仅限于微秒级（甚至更慢），且需要对选定反应进行增强采样。KIMMDY则通过仿真反应，即不显式建模反应，从而以较低的计算成本达到任意时间尺度。
创新之处在于利用最近开发的图神经网络（GNN）来完成两个新任务：（i）以构象感知的方式仿真反应及其势垒，采用混合MD/ML方法；（ii）学习分子力学力场，以在反应环境中进行物理基础的集合生成。通过与实验和DFT数据的对比验证了反应仿真器。
与使用马尔可夫状态模型（MSM）描述反应动力学的方法相似，KIMMDY不需要计算量大的从头计算或MLIP MD模拟，而是允许使用MM力场。这样大大减少了计算成本，并使得可以研究更大规模的生物分子系统，如胶原蛋白纤维。此外，KIMMDY可以方便地应用于具有相似化学特性的系统，而无需事先进行计算。然而，KIMMDY将可访问的反应空间限制为由反应插件预定义的反应类型，而从头计算方法可以发现所有可能的反应和过渡态。
在三个典型应用中，KIMMDY提出了意外的预测，重新诠释了以往的实验，并激发了新的实验。首先，在胶原蛋白中识别出一个新的高效自由基稳定基团，这为解释至今未解的光谱数据提供了可能，并对理解组织力学化学和衰老具有重要意义。其次，KIMMDY展示了在从简单到复杂的结构环境中，开放壳和闭壳生物化学之间的竞争可能发生显著变化，从而调和了以往观察中的表面矛盾。第三，KIMMDY预测DNA中光诱导二聚化的量子产率在常见的DNA折纸图案中出乎意料地低，并且与先前的邻近二聚化反应无关——这些结果是可验证的，且对生物医学应用具有高度相关性。
尽管KIMMDY能够模拟多种系统和反应，但仍然存在固有的局限性。KIMMDY的反应是为特定反应类型选择并实现的，它无法预测研究者未预见的反应类型。例如，二聚化反应可以扩展以考虑所有可能的光产物。为了在粗粒度模型中模拟反应网络的广泛范围，可以在模拟反应之前，广泛探索反应网络，当前的技术仍主要依赖启发式方法。通过使用机器学习引擎，将反应速率从全原子表示转移到粗粒度表示，KIMMDY可以在接近DFT精度的粗粒度模拟中实现生物分子反应性模拟，这将是未来研究的主题。另一个例子是用从头计算数据学习来替代启发式速率，在激发态能量表面上基于过渡态估计DNA中的二聚化速率，这应该能提高仿真器的准确性。水解模型使用与溶剂可接近表面面积的线性关系，仍然是一个假设，机器学习仿真器可以消除这一假设。然而，KIMMDY适用于可以识别高效且足够准确的启发式模型的问题。
总之，KIMMDY是一种多功能且高度可扩展的计算方法，用于探索动态分子系统中生物化学反应性。在KIMMDY的大规模模拟中，生物化学反应性几乎不增加额外成本。预计KIMMDY将在多个领域提供新的见解，从酶在细胞内如何在严格调控下工作，到生物分子如何被永久性地修饰、损伤、降解或修复。