基于片段的药物发现（FBDD）在药物开发早期阶段变得越来越流行。计算工具在这些活动中起着关键作用，为库设计、虚拟筛选、小分子潜在结合位点的识别、结合几何结构的解析以及精确结合亲和力的预测提供了途径。在此背景下，基于分子动力学的模拟日益受到欢迎，但通常受到采样问题的限制。本文提出了一种新的巨正则非平衡候选蒙特卡洛（GCNCMC）方法以克服这些限制。GCNCMC尝试将片段插入或从感兴趣区域删除；每个提议的移动都需通过基于系统热力学性质的严格接受测试。我们证明了基于片段的GCNCMC能够高效找到遮蔽的片段结合位点，并准确采样多种结合模式。最终成功计算了片段的结合亲和力，而无需使用限制条件、处理多种结合模式或进行对称性校正。
图1给出了GCNCMC协议的高层次概述。简言之，整个协议包括运行常规分子动力学（MD）来推进系统，其间穿插GCNCMC移动。以相等的概率选择插入或删除移动，对于插入移动，非交互分子被随机放置到用户定义的GCMC区域并“打开”。对于删除移动，完全交互的分子在GCMC区域内随机选择然后“关闭”。这些开关由一系列炼金术扰动和每次扰动之间的少量MD控制，使得整个开关过程处于非平衡状态，即模拟在移动过程中没有时间达到平衡。在GCNCMC移动结束时，执行蒙特卡洛接受测试，根据系统的热力学性质，包括过剩化学势和所需浓度，决定移动被接受或拒绝。概念上，较高浓度会促进插入移动的接受，反之亦然。如果移动被接受，模拟继续从该状态开始；如果移动被拒绝，模拟则从尝试移动前的状态重新开始。
GCNCMC模拟准确反映水中的片段浓度
在之前的工作中，我们通过再现TIP3P水盒的质量密度验证了GCMC和GCNCMC方法。对于片段-水混合物而言，测量每种片段在水中的总体浓度更为合适。我们选择了0.5 M丙酮和0.1 M嘧啶溶液进行测试，因为它们在这些浓度下不会在水中聚集。这些测试的初始浓度为纯水盒不含其他物质，以及含有1 M丙酮和0.5 M嘧啶的溶液盒。为了进行这些模拟，我们设置了适合目标浓度的GCNCMC参数（见方法部分）。
图2显示了两种片段随模拟时间变化的浓度。在每种情况下，经过适当的平衡期后，浓度围绕目标值波动，表明GCNCMC模拟不仅能够维持定义的浓度，还能快速使系统达到平衡。丙酮的最终平均浓度为0.55 ± 0.02 M和0.56 ± 0.02 M，分别从0 M和1 M开始。虽然略高于预期的0.5 M，但两个系统结果的一致性令人欣慰，表明配体或水的过剩化学势值可能不足以精确采样0.5 M。对于嘧啶，选择了较低的0.1 M浓度并得到了良好重现。有关这些结果的进一步讨论见补充信息。
GCNCMC快速准确识别遮蔽结合位点
一种流行的寻找潜在可药用位点的计算机方法是通过MSMD模拟，其中蛋白质在含有高浓度有机探针分子的水中溶解。探针停留时间较长的区域被认为具有潜在可药用性。这些模拟通常适用于溶剂暴露或位于蛋白质表面的蛋白质结合位点，但对于遮蔽结合位点或更深层埋藏的结合位点，所需的模拟时间超过了方便模拟的时间尺度。使用GCNCMC，可以通过在蛋白质附近直接进行插入和删除操作来消除扩散时间的限制。将GCMC区域设置为覆盖整个或部分蛋白质意味着可以仅集中在这些区域进行采样。此协议本质上是传统MSMD模拟与GCNCMC增强采样的结合。在此，我们应用这种方法找到了T4-溶菌酶（T4L99A）和主要尿蛋白1（MUP1）的遮蔽结合口袋。为了做到这一点，我们将GCMC区域锚定在一个中心蛋白质残基上，以覆盖整个蛋白质。
T4L99A是一个广泛研究的测试系统，拥有丰富的实验数据，常用于自由能方法的发展。点突变（L99A）导致了一个蛋白质腔，该腔结合了一系列小芳香配体。尽管T4L99A系统相对简单，但有一些复杂性使其成为测试增强采样方法的特别有趣的对象。例如，由于结合位点完全遮蔽于溶剂之外，MSMD模拟难以通过经典MD映射该位点，因为需要的时间尺度超出了实际模拟的范围。此外，一些配体如甲苯在T4L99A上的结合存在显著的动力学障碍，阻碍其容易互换。在此，我们首先测试GCNCMC作为位点查找工具的能力，并将其映射遮蔽结合口袋的能力与经典MSMD模拟进行比较。
在我们的六次重复实验中，平均在34次GCNCMC移动后轻松找到了苯的结合位点。作为参考，我们在GTX1080 GPU上运行24小时墙钟时间，每次重复进行700次移动。一旦找到位点，配体在整个模拟过程中保持结合状态，因为在这个浓度下每次删除提议都被认为是不利的。通过对采样的苯重原子坐标进行网格化分析，我们可以统计每个网格点上苯原子出现的帧数，然后基于总帧数进行平均。通过将网格轮廓设定为表示至少90%的帧占有率（图3），我们看到晶体学结合姿态周围的清晰信号，表明在至少90%的模拟时间内苯分子存在于这个位点。此外，由于没有其他网格点占据如此高的百分比，结合位点清晰且无假阳性。
这些结果与基本MSMD模拟的T4L99A在0.5 M苯溶液中的结果进行了比较，如图3所示。MD网格在30%处轮廓化，表示苯原子停留时间超过30%的仿真帧的网格点，显示苯结合位点未被采样。即使将轮廓水平降低到1%，依然没有任何结合。这种基于网格的分析可以用来估算自由能，通过将其占有率与大块溶剂进行比较，类似其他MSMD方法。
GCNCMC自动捕获多种结合模式，无需先验知识
一旦已知结合位点，了解片段在该位点的所有可能结合方式是有帮助的。实验数据，包括晶体结构和结合亲和力，是对所有可能构型的系综平均。然而，在MD模拟中，改变结合模式相关的动能屏障非常高，使得这一过程难以采样。如果药物化学家不完全理解片段如何结合，可能会丢失如何优化或利用某些结合相互作用的关键信息。在此，我们探索GCNCMC在无需任何先验知识的情况下采样多种结合模式的能力。
主客体系统，如β-环糊精（βCD），是方便且易于处理的复合物，常用于测试新的模拟方法，并表现出许多与蛋白质结合片段相同的特性。由于其体积小，主客体系统的模拟通常收敛迅速。据报道，具有单一极性基团的客体分子以极性基团指向宿主两端的不同取向结合到βCD上，配体在次级醇取向上更有利地结合（图5）。为了检查这些结合模式是否出现在我们的GCNCMC模拟中，我们从滴定研究中提取框架（参见“GCNCMC滴定可以准确排名结合亲和力”部分），针对两个代表性片段——苯甲腈（图5）和对甲酚（补充图10）。具体来说，我们查看了大约50%占有率的B值（B_50）下的模拟，因为这是对应于解离常数_K D的B值（参见“方法”部分），并在最大数量的结合和解离事件下产生。我们叠加这些框架，看到正如预期那样，这两个配体的极性基团更倾向于指向βCD较宽的次级醇开口（图5）。
GCNCMC滴定可以准确排名结合亲和力
最后，既然我们已经验证了GCNCMC无需先验知识即可预测结合位点和特定位点内的相对结合模式群体，我们现在可以应用我们的滴定协议来计算不同片段系列与β-环糊精、T4L99A和MUP1的结合亲和力。这些片段的结构见补充信息。
对于预计算的每个配体过剩化学势，如同任何自由能计算所必需的，我们通过Adams值扫描（见方法部分）模拟一系列配体浓度。然后测量每个_B_ eq值下的平均占有率，并将逻辑函数拟合到这些数据（方程[18]）。对应的平均占有率为50%的配体浓度是解离常数_K_ D，这很容易与结合自由能相关联（方程[13]）。
图8显示了从GCNCMC滴定中提取的主客体结合自由能，与实验数据和使用平面底部约束的传统ABFE方法进行了比较。一般来说，相对于实验观察到轻微的结合亲和力高估，提示使用的力场参数可能不是最优的，这一趋势之前已有报道。尽管如此，计算得到的平均绝对误差（MAE）和均方根误差（RMSE）相对于实验分别为0.7和0.6 kcal mol−1，几乎所有数据点都在实验值的1 kcal mol−1以内。此外，与实验的相关性（R 2 = 0.94）和排名（τ = 0.84）表明该方法可以可靠且准确地按结合亲和力排列片段。与更成熟的方法FEP几乎完美的相关性提供了有希望的验证。
总结
在这项通讯中，我们进一步发展并应用了我们的巨正则非平衡蒙特卡洛方法，以采样小分子与蛋白质系统的结合，目的是提供一种可以帮助发现新片段先导化合物的工具。该方法已成功应用于检测结合位点、阐明结合几何结构和计算结合亲和力。
我们通过重现简单的系综属性（即水中的片段浓度）验证了该方法。我们展示了我们的GCNCMC实现能够准确重现用户定义的片段浓度，尽管在高浓度下，结果对所用过剩化学势值变得敏感。
我们应用该方法于两个蛋白质-片段系统T4L99A和MUP1，以增强传统的MSMD模拟。在这两种情况下，我们假设不了解结合区域，并随后运行GCNCMC模拟以生成一组可能的结合位点。我们进行了GCNCMC与基本MSMD协议的比较分析，揭示了GCNCMC在两个系统中都能轻易识别实验结合位点，而MSMD则不能。
然后我们展示了该方法如预期般重现片段结合模式。我们发现GCNCMC在主客体系统中采样简单片段的两种结合模式，我们的模拟显示客体分子以与实验和其他报告数据一致的方向结合。对甲苯与T4L99A结合的更详细检查显示，该方法成功复制了四种结合模式，其群体与先前发表的研究一致。
方法
理论
巨正则非平衡候选蒙特卡洛
基于配体的巨正则非平衡候选蒙特卡洛（GCNCMC）建立在Samways等人和Melling等人的工作基础上。最终结果在这里展示，详细推导见补充信息。
立方系统中GCNCMC插入和删除移动的接受概率如下（使用Adams公式）：
[ P_{\text{insert}} = \min \left[1, \frac{1}{N+1} e^{B_{\text{eq}}} e^{-\beta w(X| \Lambda_p)} \right] ]
(1)
[ P_{\text{delete}} = \min \left[1, N e^{-B_{\text{eq}}} e^{-\beta w(X| \Lambda_p)} \right] ]
(2)
其中 ( w(X| \Lambda_p) ) 是由协议 ( \Lambda_p ) 生成轨迹 ( X ) 所做的功，或换句话说，非平衡开关所做的功，( N ) 是移动前初始状态中分子的数量，( B_{\text{eq}} ) 是Adams值。对于特定的片段浓度，( B_{\text{eq}}(c) ) 定义为：
[ B_{\text{eq}}(c) = \beta \mu_{\text{sol}}' + \ln \left(\frac{V_{\text{GCMC}}}{V(c)}\right) ]
(3)
其中 ( \mu_{\text{sol}}' ) 是配体分子的过剩化学势，通常近似为其无限稀释的标准溶剂化自由能。( c ) 是参考溶液中分子的浓度，( V(c) ) 是指定浓度下溶液中每个分子的平均体积，( V_{\text{GCMC}} ) 是GCMC区域的体积。
对于给定的参考浓度，每个分子的平均体积可以简单计算为：
[ V(c) = \frac{1}{c N_A} ]
(4)
其中 ( N_A ) 是阿伏伽德罗常数。
最后，非平衡移动的长度由两个最终状态（λ = 0 和 λ = 1）之间的总扰动步数 ( n_{\text{pert}} )，以及每次扰动之间的MD传播步数 ( n_{\text{prop}} ) 控制。总的切换时间为：
[ \tau = (n_{\text{pert}} + 1) n_{\text{prop}} \delta t ]
(5)
其中 ( \delta t ) 是积分器的时间步长。为了保持详细的平衡，每次GCNCMC移动开始和结束时都有一个传播步骤。
GCMC球体
为了提高特定感兴趣区域的采样，用户可以通过定义一个球形的“GCMC区域”来定位GCNCMC移动。例如，这可以是一个已知的结合位点或整个蛋白质。然而，使用球体进行GCNCMC移动需要特别小心，因为在整个切换过程中，分子可能扩散进出球体，这意味着必须调整之前定义的接受率：
[ P_{\text{insert}} = \min \left[1, \frac{1}{N_T} e^{B_{\text{eq}}} e^{-\beta w(X| \Lambda_p)} \right] ]
(6)
[ P_{\text{delete}} = \min \left[1, N_0 e^{-B_{\text{eq}}} e^{-\beta w(X| \Lambda_p)} \right] ]
(7)
其中 ( N_0 ) 是移动开始时球体内GCMC分子的数量，( N_T ) 是移动结束时的分子数量。应该注意的是，如果被切换的分子在移动结束时位于球体外部，则必须自动拒绝该移动，因为无法提出反向过程，破坏了详细平衡条件。根据结合位点的特征，有时会导致大量移动被拒绝。
Adams值
在前面的部分中，我们介绍了Adams值 ( B_{\text{eq}} )。这是GCNCMC中的控制参数，最终决定了是否接受GCNCMC移动，因此其定义需要特别关注。
在我们之前的水GCMC和GCNCMC研究中，Adams值定义为：
(8)
其中 ( \mu_{\text{sol}}' ) 是分子的过剩化学势，通常近似为分子的无限稀释标准溶剂化自由能。水和小分子的标准状态通常定义为55 M和1 M，但在许多情况下，模拟一个分子（如片段）在非标准状态下的情况更具实验相关性。例如，类片段分子通常在微摩尔至毫摩尔范围内与其靶标结合。在这种情况下，当参考溶液中的分子偏离标准状态时，我们定义了具有特定浓度依赖性的Adams值：
[ B_{\text{eq}}(c) = \beta \mu_{\text{sol}}' + \ln \left(\frac{V_{\text{GCMC}}}{V(c)}\right) ]
(3)
其中 ( V(c) ) 是现在在浓度 ( c ) 下每个分子的平均体积。该方程也可以直接用浓度表示为：
[ B_{\text{eq}} = \beta \mu_{\text{sol}}' + \ln (N_A c_L V_{\text{GCMC}}) ]
(9)
其中
[ c_L = \frac{1}{N_A V(c)} ]
(10)
配体的过剩化学势定义为将分子添加到给定溶液中的自由能，在确定特定片段是否会偏好溶液或结合位点方面至关重要。对于本研究中的结合研究，我们近似特定片段的过剩化学势等于该片段的“无限稀释”水合自由能，或者换句话说，将片段分子添加到水盒子中的自由能。不过，应该强调的是，( \mu_{\text{sol}}' ) 的值也可能受到与其他同类型分子相互作用的影响，并且理论上也有浓度依赖性。然而，对于足够稀的浓度，如这些研究所使用的浓度，假设这些相互作用的影响可以忽略不计，并且 ( \mu_{\text{sol}}' ) 与浓度无关。因此，对于给定分子，只需要计算一个 ( \mu_{\text{sol}}' ) 值，并且该值可以应用于该分子的任何稀释浓度。关于此点的完整讨论和数据可以在补充信息中找到。
在 ( \mu_{\text{sol}}' ) 固定的情况下，Adams值仅依赖于参考溶液中的片段浓度。现在直观地说，较高的参考溶液浓度会导致蛋白质GCMC区域中更多的结合，如图10所示。随着浓度增加，结合事件的最大数量将保持在50%的占用率，从而导致最大的结合和解离事件。
图10：GCNCMC滴定示意图。
在多个B值下进行模拟，并绘制每个模拟的最终平均结合位点占用率，以便可以使用方程(15)计算 _B_50。简单的转换允许相同的数据在浓度尺度上绘制（右侧），以便可以使用方程(18)计算 _K_D。
(全文结束)
