DRUGAI
今天为大家介绍的是来自美国圣地亚哥加利福尼亚大学Rommie E. Amaro团队的一篇论文。分子动力学模拟在探索复杂的生物过程中已变得不可或缺，但其在捕捉罕见事件方面的局限性阻碍了作者对药物-靶点动力学的理解。在这篇观点文章中，作者研究了“里程碑模拟”（milestoning simulations）领域以理解这一挑战。里程碑（milestoning）方法将药物-靶点复合物的相空间划分为离散单元，从而提供了更长时间尺度的见解。本文追溯了里程碑模拟在药物-靶点动力学研究中的历史、应用和未来潜力。文章探讨了里程碑方法的基本原理，强调了概率转移和转移时间独立性的重要性。通过重新审视具有Voronoi镶嵌的马尔可夫里程碑方法来解决传统里程碑方法的挑战。在关注该领域进展的同时，本文也指出了通过里程碑模拟估计药物-靶点解离速率常数时面临的紧迫挑战，为更有效的药物设计策略铺平了道路。
原子级分子动力学（MD）模拟已经成为研究复杂生物过程（如蛋白质折叠、蛋白质-膜相互作用和药物-靶点相互作用）的一种强大工具。随着算法、软件和硬件能力的显著提升，模拟的时间尺度得以延长，从而能够观察这些过程。然而，严格表征这些过程的动力学和热力学仍然具有挑战性。由于MD模拟的时间步长受限于飞秒级，这限制了其捕捉生物系统中复杂构象转变或罕见事件（毫秒级或更长时间）的能力。为弥补这一时间尺度的差距，许多不同类型的增强采样方法被开发出来。增强采样方法大致可以分为两类。第一类通过引入偏置势来加速构象态之间的转变。在某些情况下，包括复制交换MD、选择性集成回火和高斯加速MD，偏置作用于整个系统。在其他情况下，包括元动力学和变分增强采样，偏置只选择性地应用于系统内的某些力（所谓的"集体变量"）。第二类增强采样方法称为路径采样方法，重点在于对转变区域进行采样。这些方法包括转变路径采样、加权系综（ensemble）模拟、转变界面采样和里程碑模拟。
药物与靶点的结合和解离动力学是多种生物现象的重要决定因素，包括酶催化、细胞信号传导和免疫系统激活。这些分子相互作用使细胞能够将外部刺激转化为生化信号，从而维持基本的生理功能。具有较长停留时间的药物会在靶点的活性位点内停留更久，从而延长其药理作用；相反，快速解离的药物通常药效较弱。因此，全面描述动力学特征对靶向治疗药物的合理设计和优化至关重要。尽管MD模拟在研究多种生物过程方面取得了重要进展，但其在药物-靶点动力学中的应用仍面临独特的挑战和机遇。药物与靶点的相互作用通常发生在不同的时间尺度上，可能持续数秒、数小时甚至更久。时间尺度的限制使得传统MD模拟难以捕捉这些罕见事件。此外，结合位点的较大构象变化、力场对结合态极化效应的准确性要求，以及实验验证的需求，进一步增加了研究药物-靶点动力学的难度。
最近，里程碑理论已成为研究药物-靶点动力学的一种稳健且高效的方法。里程碑方法通过将复杂过程划分为一系列称为里程碑的过渡状态来简化分析。这些分布在系统构型空间中的里程碑突破了常规MD模拟的局限性，可以在更长的时间范围内模拟生物分子过程的动力学。本文介绍了药物-靶点动力学里程碑模拟的进展和挑战，概述了里程碑模拟在药物-靶点动力学研究中的历史、应用、挑战和未来前景。
传统里程碑方法
里程碑方法是一种计算方法，它通过将复杂过程分解为系统构型空间中一系列称为“里程碑”的中间状态之间的转变来简化分析（图1a）。
图1
这种方法能够在更长的时间尺度上模拟系统的整体动力学，突破了传统MD模拟的限制。里程碑理论基于两个基本假设。首先，在特定模拟时间t'时，里程碑i(t')之间的连续转变遵循概率序列、、，这由在时间t内从一个里程碑i成功转变到下一个里程碑j的概率决定，确保这种进展的整体统计行为类似于离散时间马尔可夫过程（方程1）：
方程1
在这个过程中，步骤的顺序至关重要，每次转变的概率取决于前一次转变的结果，从而形成一个连贯且可预测的序列。其次，连续里程碑交叉之间的转变时间τ在统计上是相互独立的。基于这些假设，观察到给定轨迹（给定转变速率）的概率p取决于以下几个部分：模拟中观察到的转变次数、从跨越里程碑i之后到跨越其他里程碑之前所花费的总时间，以及具有分量的速率矩阵Q，这构成了方程2的里程碑理论基础。
方程2
虽然这个公式假设存在单一的长轨迹，但在实际应用中，里程碑法需要在每个里程碑上重新初始化短时分子动力学模拟。特定的里程碑组合（被定义为最优里程碑）能够精确满足第一个假设，这简化了对平均首达时间（MFPT）等动力学特性的分析。然而，最优里程碑的确定和近似非常困难。因此，一种替代方法是通过在里程碑之间反复传播轨迹来实现精确的里程碑方法，这样就不需要最优里程碑。尽管这种方法不依赖任何假设，但需要注意的是，这种计算的精确性需要付出更高的计算成本，通常要增加一个数量级。
在过阻尼系统中，里程碑法非常适合准确地捕捉状态之间的概率跃迁。这种方法考虑了系统动力学的物理细节，确保所选择的里程碑能够进行精确的动力学分析。如果近似放置最优里程碑，或使用精确里程碑法等方法来补偿非最优里程碑，就能确保里程碑之间具有特定的概率跃迁模式。最优里程碑被定义为提交（commitor）函数的等值面（isosurface），这是跃迁路径理论中的一个重要概念。使用等提交面（isocommittor surfaces）作为里程碑可以保证跃迁过程中所需要的概率行为。虽然第二个假设可能并不普遍适用，但使用这些最优里程碑仍然可以精确计算平均首达时间，使其在复杂系统的动力学分析中发挥重要作用。
通过里程碑重新初始化MD轨迹可以补偿里程碑位置设置的不理想情况。虽然理想的里程碑能够确保在里程碑之间的跃迁概率保持真正的恒定，无论系统是如何穿过里程碑表面的，但理想里程碑仅能被近似，且获取成本较高。因此，需要通过有效的方式从里程碑重新初始化MD模拟，以确保产生正确的跃迁概率和跃迁时间。
可以通过估算首次命中点分布（FHPD）来解决缺乏理想里程碑的问题。FHPD依赖于位置 xx处的密度，表示从里程碑i开始的轨迹首次命中里程碑j的概率。在理想的里程碑集合中，密度不依赖于起始里程碑i。为了获得FHPD，可以使用跃迁路径理论（TPT）作为框架，其中FHPD的密度通过公式3定义。
公式3
使用里程碑方法的药物-靶点动力学研究
动力学速率模拟估算（Simulation-Enabled Estimation of Kinetic Rates，SEEKR）是一个用于估计和的多尺度里程碑框架。SEEKR使用原子级的分子动力学（MD）模拟以及较简化但计算效率更高的布朗动力学（BD）模拟，来估算药物与靶点结合及解离的动力学速率常数。SEEKR工作流程从一个结合态（bound state）的受体-配体复合物的三维结构开始。定义一个集体变量（CV），通常是配体与结合态中受体的α碳原子之间的质心距离。基于该CV，将复合物的相空间划分为MD和BD区域。在MD区域中，进一步划分出里程碑之间的Voronoi单元。通过运行受控MD（SMD）或元动力学（metadynamics）模拟，从结合态逐步移动到非结合态，保存每个Voronoi单元中的受体-配体复合物结构。
在里程碑之间的空间内运行独立且并行的MD模拟，而在BD区域运行BD模拟。通过施加反射边界条件，将轨迹限制在各自的里程碑之间。根据每个Voronoi单元中的转变次数和停留时间，构建转变矩阵，进而计算药物-靶点的解离速率常数（）。BD模拟提供了一种高效的计算方法，用于估算药物-靶点复合物的结合速率常数（）。与MD模拟不同，BD模拟采用简化模型和近似，例如将分子视为刚体并使用隐式溶剂，这使得模拟能够采用更大的时间步长。BD模拟能够准确地模拟药物结合事件的初始阶段，并考虑显式的静电相互作用。
SEEKR框架集成了Browndye模拟包，该包采用Luty-McCammon-Zhou算法生成配体围绕受体的多条轨迹，用以估算结合速率。这些轨迹最终以分子接触或逃逸结束，通过接触概率计算结合速率常数。该方法在分子识别和结合中静电力显著的情况下特别有效。
SEEKR框架在多个受体-配体复合物中展示了与实验接近的koff速率，包括停留时间从微秒到数小时不等的情况。这些应用涵盖了从较为简单的系统（如β-环糊精-配体复合物和胰蛋白酶-苯胺基胍复合物）到更具挑战性的系统（如JAK-抑制剂复合物、一系列Hsp90-抑制剂复合物，以及一系列苏氨酸-酪氨酸激酶（TTK）抑制剂）。
对于药物-靶点动力学的一个有效问题是，Milestoning方法是否能够有效捕捉复杂生物系统中灵活性和构象变化的影响，而这些因素并未在Milestoning模型中明确表示。如果灵活运动的时间尺度相对于单次Milestoning模拟的持续时间较短，且这些运动能够在模拟时间内被充分采样，那么这些运动将自动包含在Milestoning模型中。然而，如果这些运动较慢且在单次Milestoning模拟中未被充分采样，而这些缓慢的灵活运动会直接影响所研究过程（如结合/解离）的时间尺度或能量学，那么这些运动应通过CV明确表示并加入Milestoning分析。如果研究某系统需要超过一个CV，则需要使用多维Milestoning，例如通过Voronoi镶嵌方法。这一功能在SEEKR软件中已被实现。
尽管SEEKR框架通过多尺度方法简化了和速率的计算方法，但值得注意的是，一些独立的里程碑研究专门针对速率常数的计算。这些研究通过探索不同的药物−靶点复合物和场景，进一步加深了对药物−靶点结合与解离动力学的理解，并确立了这些方法作为估算koff速率的高效手段。图2展示了在里程碑模拟中使用的各种基准系统中，受体及其对应抑制剂所遇到的复杂性水平。这些系统包括β-环糊精−配体、胰蛋白酶−苯胺胍、JAK−抑制剂、TTK−抑制剂以及Hsp90−抑制剂复合物。
图 2
多种因素导致受体−配体复合物的复杂性，并凸显了使用新颖模拟方法准确捕捉这种动态行为的必要性。配体的复杂性包括分子大小、化学多样性、立体化学、动态行为、结构灵活性以及其相互作用中的特异性和选择性等因素；而受体的复杂性则取决于其构象多样性、变构调节位点、蛋白−配体相互作用强度，以及结合位点及其相关环或铰链的灵活性。
一项最新研究使用里程碑方法研究了格列卫从Abl激酶解离的分子通路和绝对解离速率。Abl激酶是癌症治疗中的关键靶点。利用里程碑模拟方法，估算出平均MFPT（1/koff）为0.055秒，与实验测得的驻留时间0.04 ± 0.01秒接近，总模拟时间为1微秒。
里程碑的设置采用了一种逐步推进的方法，首先通过SMD模拟生成了初始的43个里程碑。随着从这些初始里程碑出发的非约束轨迹偏离了SMD路径，后续又发现了新的里程碑，并重复这一过程，直到获得一个连通的转移矩阵和有限的MFPT。MFPT的计算方法是汇总从反应物（结合态）到产物（解离态）的所有可能路径，同时考虑里程碑之间调整后的转移概率。
面临的挑战与未来方向
（1）力场的准确性。药物分子在结合位点的相互作用涉及电荷作用、范德华力、氢键和极化等多种分子力，是一种复杂的相互作用。结合态的极化现象（即电子密度响应外部电场的重新分布）通常会影响药物−靶点的相互作用。因此，研究并开发更优的固定电荷力场用于药物发现和设计已成为一个重要的研究方向，例如Open Force Field Initiative的相关工作。
尽管药物分子和结合位点残基会因为电子极化而改变其电荷分布，但在特定的里程碑计算中，非极化力场在某些系统中依然表现出色。精确估算极化效应非常关键，因为它会影响结合亲和力、动力学行为以及药物的治疗效果。极化力场（如AMOEBA和CHARMM Drude力场）专门用于捕捉受体−配体复合物中的电子极化现象，能够更真实地描述极化效应对结合与解离动力学的影响。
（2）确定配体解离自由能的最佳估算。在里程碑模拟中，初始条件的选择，尤其是起始里程碑的设置至关重要。传统方法通常采用SMD模拟，这可能导致当配体从结合口袋移出时，受体−配体的动态表现出较大的张力。因此，里程碑模拟的一个关键环节是准确估算配体从结合位点解离的自由能曲线。
耗散修正的目标分子动力学方法（Dissipation-corrected targeted MD，dcTMD）是一种新颖的计算慢过程自由能曲线的方法。它通过对能量耗散进行修正，解决了模拟复杂过程（如配体从结合位点解离）的准确性问题，从而为系统动态提供更精确的估算。dcTMD通过考虑模拟过程中摩擦的累积效应，修正能量耗散，确保计算出的自由能曲线能够准确反映因耗散引起的能量损失。
dcTMD对受体−配体复合物中的配体施加偏置力，作为一种引导机制，驱动配体沿预定义反应坐标进行解离。该方法利用了Jarzynski等式，通过对系统所做功的指数平均值与解离机制相关的自由能变化建立联系。dcTMD的一大特点是可以将总反应通量分解为多个路径，每条路径都有其独特的能量曲线和摩擦因子。这种分解能力使dcTMD能够处理复杂的自由能景观，尤其适用于能量景观存在多条路径的系统（如T4溶菌酶和HSP90等）。通过为模拟提供优化的初始条件，dcTMD确保里程碑模拟从更真实和具有代表性的状态开始，从而提高了模拟的准确性。
随机加速分子动力学（RAMD）是一种增强采样方法，通过向配体施加随机方向的力，加速蛋白−配体解离事件。在分子动力学模拟中，RAMD向配体施加一个大小恒定、方向随机的力。这个力的初始方向是随机选择的，并可能根据配体相对于受体蛋白的运动而变化。RAMD涉及几个关键参数，包括配体运动检查的时间间隔、配体−受体位移阈值、配体的最大位移以及所施加力的大小。利用RAMD模拟可以高效采样配体的解离路径，从而生成自由能曲线。配体解离路径的模拟快照可以保存为里程碑模拟的起始结构。最近，一种改进版的RAMD方法——τ-RAMD被用于估算药物−靶点相对解离动力学，并根据驻留时间对一系列Hsp90−抑制剂复合物进行排序。在里程碑模拟中引入RAMD或τ-RAMD，用于预测初始里程碑，可以帮助准确估算动力学速率常数。
（3）找到准确的集体变量估算方法并采用多维里程碑方案。复杂的生物过程，例如药物与靶点的结合或解离，通常伴随着多次构象变化。使用单一维度的集体变量（CV）可能会过度简化这些过程，从而可能忽略系统中更高维度的构象配置和转变。换句话说，如果选择的CV未能代表系统中最缓慢的过程，里程碑法可能会导致动力学结果产生偏差或不准确。多维里程碑法在某些情况下具有明显的优势。
在里程碑法中使用多维CV可以描述更复杂且具有曲度的反应路径。通过在多维网格中结合两个或更多的CV，作者能够展现一个更高维的自由能景观，从而对能量景观进行更准确和详细的探索。此外，某些能量障碍只有在使用多维反应路径时才会显现出来。多维CV还便于观察涉及多个维度协调运动的过程，从而为反应路径、转变速率和能量景观提供可靠的洞察。
图3展示了Muller势的能量景观。该图比较了一维CV方法（沿状态A和状态B之间的最小能量路径标记）与二维Voronoi剖分法。后者能够在复杂的能量景观中提供更可能且更详细的路径表示。
图 3
（4）量子/分子力学的动力学模拟（QM/MM MD）结合了量子力学计算的精确性和分子力学模拟的高效性，用于研究复杂环境中的生物过程。在这种方法中，一部分原子（通常是活性位点或感兴趣的区域）采用量子力学方法处理，而其余原子或残基则通过分子力学进行经典处理。
在药物与靶点结合或解离过程中，某些机制（如质子转移、电子重排、键的形成与断裂）可以通过靠近结合位点区域的QM/MM MD模拟进行定量估算。QM/MM MD模拟特别适合于确定过渡态，即沿集体变量（CV）高能量区域，从而提供有关解离过程限速步骤的重要信息，这对于理解过程的动力学至关重要。
CV沿线的最初内层里程碑通常位于药物与靶点结合位点附近，这为应用QM/MM MD模拟提供了机会。虽然在这些里程碑范围内的QM/MM模拟需要大量计算资源，但研究其在模拟时间尺度上提升精确性的潜力是一个有趣的方向。
（5）扩展里程碑生成轨迹的应用范围。现有的增强采样方法已被用于预测药物-靶点动力学，并适用于多种靶蛋白。然而，这些靶蛋白中许多在结构复杂性和化学组成上差异较大。胰蛋白酶-苯胍复合物作为一种基准系统，被广泛用于方法开发研究。然而，由于其作为药物-靶点复合物的简单性以及缺乏隐性集体变量（CV），该系统存在一定局限性。在这个系统中，解离过程相对简单，包括从结合口袋中的过渡，有时伴随表面扩散，最终使苯胍进入溶剂中。
然而，最近的一项研究揭示了胰蛋白酶-苯胍复合物中更多的复杂性，特别是水分子的作用。研究中使用的CV主要聚焦于结合腔内水分子的行为和排列。这一发现挑战了之前关于胰蛋白酶-苯胍作为简单药物-靶点复合物的理解，强调了在预测药物-靶点动力学时需要考虑这些附加变量，以获得更准确的结果。
为了全面评估新的增强采样方法，基准系统应包含具有不同化学骨架的多样化配体。此外，这些系统还需要高分辨率的实验结构和动力学数据作为支持，以便进行全面的评价。最近的研究在Hsp90-抑制剂复合物上取得了显著进展。这些配体具有多样化的骨架结构、不同的停留时间以及多条结合与解离路径。值得注意的是，这些复合物在配体解离过程中表现出显著的构象变化。这些特点使Hsp90-抑制剂复合物成为用于评估新型增强采样方法的最全面和稳健的基准系统之一。
此外，TTK系统的八种抑制剂也提供了一个医学相关性较高但规模较小的基准化合物集合，可用于动力学预测。图4展示了70种Hsp90抑制剂的散点图，根据其解离速率常数和骨架多样性进行了分类，展现了该基准系统的多样性和复杂性。
图 4
编译|黄海涛
审稿|王梓旭
参考资料
Ojha, A. A., Votapka, L. W., & Amaro, R. E. (2024). Advances and Challenges in Milestoning Simulations for Drug–Target Kinetics. Journal of Chemical Theory and Computation.
