摘要：2012年辉瑞公司的Trujillo等人在 “双水记（a tale of two waters）”一文描述了对hH-PGDS结合口袋的两个氢键网络水分子进行替换与置换实验，结果表明对主水分子（W1）的替换设计导致化合物的活性降低几百倍，而对辅助水分子（W2）的置换仅导致活性稍有降低。本文用Flare的GIST分析结合口袋的水能量特征，尝试从水的能量学角度来解释“双水记”的实验现象。结果发现主水分子W1是happy水，而辅助水分子W2是unhappy水，这可解释“双水记”观察到的实验现象。Flare的GIST水分析理论计算对实验设计具有正确的指导作用。

肖高铿/2021-07-25

双水记的由来

之前，我们分享了用FLARE 3D-RISM探测BTK激酶结合位点（PDB 4ZLZ）识别“unhappy”的桥连抑制剂与蛋白相互作用的水^[1]，然后我们用SPARK的水分子替换（replace）策略将"unhappy"水分子与蛋白间的相互作用引入到配体中，从达到提高新化合物活性实现结构优化的目的^[2]。

然而，并非对桥接配体与蛋白相互作用的水分子的替换（将水分子用配体的一部分替换，并将水分子与蛋白的相互作用合并到新配体里）或置换（只是用配体原子替换了水分子的位置，水分子被排除到溶剂中，水分子原有的相互作用没被整合到配体里）皆可提高化合物活性。比如，当该水分子为“happy”的时候，意味着该水分从热力学性质上讲更愿意呆在结合位点里，对这样的水分子进行替换或置换可能对化合物结合亲和力的提高没有帮助。2012年，辉瑞公司的Trujillo等人在 “双水记（a tale of two waters）”一文就分享了对hH-PGDS抑制剂先导化合物尝试进行水分子置换与替换但不能改善化合物结合亲和力的例子^[3]。

表1. 双水记涉及的4个主要化合物

图1展示了化合物9（见表1）与hH-PGDS共晶结构（PDB 4EE0）揭示的分子间相互作用。在结合位点里有两个发生氢键网络的水分子W1（即双水记之primary water, A HOH 320）和W2（即双水记之auxiliary water, A HOH 305）与化合物9的R1基团有紧密接触，这两个水正是“双水记”的主角。其中W1通过氢键网络桥连着配体吡啶环的氮以及蛋白LEU199的羧酸氧以及水W2发生氢键网络相互作用；W2水分子介导了GLY13与其它水分子的氢键网络。Trujillo等人^[3]以化合物9为起点分别设计了对W2水分子置换的化合物11（表1）、对W1水分子替换的化合物13与14（表1）以探讨这种置换或替换对活性的影响。

图1. 化合物9与hH-PGDS的共晶（PDB 4EE0）的结合模式以及关键的水分子示意图。其中被虚线红圈标记两个水是我们关注的，为了方便起见，将在下面的与吡啶环氮相互作用的那个称为W1(A HOH 320)，在上面的那个水称为W2(A HOH 305)。W1与W2水分子在双水记中分别被称为primary water与auxiliary water。

“双水记”的一个实验设计是对W2水分子进行置换，通过在化合物9吡啶氮的邻位引入一个甲基（化合物11）来替换W2水。如图2所示，晶体学研究（PDB 4EDZ）表明，化合物11的甲基确实占据了W2水分子的位置，虽然不是完全重合，但是原子间距离为2Å（小于两者的半径之和），因此可视为甲基置换了水。但该甲基并不拥有W2介导氢键网络的能力，因此这是个水分子置换（display）而非水分子替换（replace）。该水分子置换使得W2水分子从结合位点里消失了，但是W1介导的氢键网络还存在。比起化合物9（IC₅₀=2.34nM），水分子置换策略设计的化合物11对hH-PGDS的活性没有提高，而是下降了3.5倍（IC₅₀=8.26nM）。此外，通过叠合PDB 4EE0与4EDZ会发现化合物11比之化合物9发生了轻微的旋转以适应新引入的甲基。

图2. 水分子置换策略设计的化合物11与hH-PGDS共晶（PDB 4EDZ）的结合模式

“双水记”的另一个实验设计是对W1水分子的替换，该策略将化合物9的吡啶环氮替换为碳，并引入胺甲基与羟甲基得到化合物13、14，以期望用阳离子N或羟基占据水并继承其氢键网络实现将W1水的相互作用引入到新化合物。如图3所示，化合物13、14与hH-PGDS的共晶结构表明，羟基与胺基确实占据了W1水分子位置，并介导了W1水分子一样的氢键网络；W2介导的氢键网络还存在，该设计完美地实现了水分子替换策略。

图3. 水分子替换策略设计的化合物13、14与hH-PGDS共晶结合模式。

然而酶学实验结果令人出乎意料，与化合物9（IC₅₀=2.34nM）相比，对水分子W1替换策略设计的化合物13与14对hH-PGDS的抑制剂活性不但没有改善，反而分别大幅下降了630倍（IC₅₀=1480nM）与360倍（IC₅₀=845nM）。这与对水分W2置换策略设计的化合物11，活性仅下降3.5倍形成了鲜明对比。

作者指出，在晶体结构中，化合物13、14连接替换W1水分子的杂原子（O与N）与平面萘环间连接臂的两面角近乎平面（分别为21与27°），而优势两面角为90°。此外，化合物13、14的萘环与苯环连接臂两面角偏离低能构象两面角约20°。以化合物14为例，其两面角为117°，而低能构象的两面角为97°，偏离了约20°。作者认为，这种过高的构象内能似乎大大超过了溶剂置换带来的熵增收益，从而导致化合物13、14的活性降低。

作者总结到，其分析和结构数据强烈地表明，若想获得对hH-PGDS最大的结合亲和力需要抑制剂与埋在活性部位的主要水分子（W1）形成氢键相互作用。此外，看来对W1溶剂分子的环境、几何以及化学限制妨碍了用抑制剂原子对其进行置换。 相反，配体-蛋白质复合物对辅助水分子（W2）的置换则耐受良好，抑制剂的结合亲和力只是温和下降一点。

黄传龙^[5]也对该水分子替换与置换也做了详细的分析，对替换与置换导致活性降低的解释与Trujillo等人^[3]基本一致。认为对水分子（W1）本身形成很好的氢键网络，即使成功被替换可能也没什么收益；相反为了替换这个水分子而引入的基团，导致分子的低能构象改变。一方面替换无收益，另一方面引入构象惩罚，可能是分子活性变差的原因(之一)。

显然，如果能在合成实验前用计算方法预见此类不可替换的水是非常有用的，因为这可以为我们避免不必要的尝试。因此，本文的目的主要是尝试从能量学角度，讨论水分子W1与W2的可替换/置换性，寻找可靠的理论指导方法。

Flare GIST水分子分析

在最近发布的Flare V5中引入了一个新的模块GIST^[4]。GIST是一种水分析方法，用于评估结合口袋的水合作用，并通过在分子动力学运行结束时对显式溶剂分布采样来计算相关的水热力学。 GIST分析的结果显示为活性位点水合作用的三维等值面映射“happy”（绿色，与负∆G相关）和“unhappy”（红色，与正∆G相关）区域。unhappy区域映射的是蛋白质活性位点内更可能的可成药区域。happy区域映射的是蛋白结合位点内较低可能成药的区域，在那里水分子更稳定，因而更难置换。在本算例中，我用Flare GIST分析化合物9结合位点，重点是确认水分子W1与W2是否可被替换。

图4. 化合物9与hH-PGDS共晶结构（PDB 4EE0）结合位点的GIST水分子分析结果。其中红色等值图ΔG大于0.5kcal/mol,绿色等值图ΔG小于-0.5kcal/mol

用Flare GISP对化合物9与hH-PGDS共晶结构（PDB 4EE0）结合位点进行的水分子分析结果如图4所示，结合位点里大部分区域被红色unhappy水覆盖，这说明该结合位点具有良好的成药性。其中，W1水分子与绿色的happy区重合，因此W1为happy水；W2水分子与红色的happy区大致重合，因此W2为unhappy水。从计算结果可以预计，替换happy的W1水是能量上不利的，对W1水的处理策略通常采用作为蛋白的一部分保留不动、设计化合物与之相互作用，或者用配体的一部分完美模拟该水（但要付出代价以踢走该水），化合物13、14对W1替换后的获益可能不如踢走水所付出的代价，因而活性大幅降低。置换W2水分子的化合物11从能量上看是有利的，但是这个置换没有“继承”到W2的氢键网络，因此活性没有增加，而是温和地下降了。

总的来说，用Flare GIST对化合物9结合位点的分析表明，W1是happy水，应当在设计中予以保留或者用配体的一部分进行替换以模拟水的作用；W2是unhappy水，可以被置换或替换。这与观察到实验结果基本一致，对happy的W1替换设计化合物13与14活性大幅降低了几百倍，而对unhappy的W2水分子进行置换设计的化合物11仅在活性上稍微降低（降低了3.5倍）；这也与Trujillo等人^[3]对这这两个水的总结一致。在本算例中，Flare的GIST水分析理论计算对实验设计具有正确的指导作用。

基于QM的两面角扫描分析

为了验证作者提到的两面角问题，用Flare的Custom FF模块在B3LYP/6-31G*理论水平下对每个两面角进行24个构象的能量扫描分析，为了加快计算速度，还将图5、图6所示的阴影部分结构裁剪掉。

对水分子替换杂原子-萘环连接臂的两面角扫描结果如图5所示，该两面角的生物活性构象比低能构象的构象能高了不仅4kcal/mol，这是极其不利的。

图5. 水分子替换杂原子-萘环连接臂的两面角扫描

对萘环-苯环的两面角扫描结果如图6所示，该两面角的生物活性构象比低能构象的构象能仅高约0.45kcal/mol左右，相当于高、低能构象Boltzmann权重比大约为30:70，从能量上讲不是特别的不利。

图6. 萘环-苯环连接臂的两面角扫描

总的来说，基于QM的两面角扫描分析表明， 化合物13、14的杂原子-萘环连接臂两面角处于高能状态，是化合物13、14活性大幅降低的另一个因素。萘环-苯环的两面角在能量上虽然也不利，但不是特别的严重。对两面角进行QM扫描分析结果与作者Trujillo^[3]及黄传龙^[5]解释基本一致。

蛋白相互作用势(Protein Interaction Potential, PIP)分析

进行水分子替换设计，不仅要与水的能量学特征匹配，化合物还要有合理的构象，而且还要形状与静电与蛋白互补。因此，对4EE0的结合位点进行了蛋白相互作用势分析。

蛋白相互作用势(Protein Interaction Potential, PIP)是Cresset分子相互作用势对蛋白质的延伸，两者都是使用XED力场计算的。该方法在原理上类似于配体场的计算：蛋白质的活性位点被探针原子充满，并且计算每个格点上的相互作用势。该方法利用了Mehler等人^[9]的距离依赖的介电函数来更好地处理蛋白质结构中的大量带电基团。鉴于我们研究的对象为同系物，因此仅计算和显示活性位点的蛋白质相互作用电势。

为了理解化合物1与2的活性差异，在Flare里计算4EE0结合位点的蛋白相互作用势，包括干口袋(不包含结合结晶水)与湿口袋（包含结晶水）的两种情况的蛋白相互作用势。其中湿口袋仅保留GIST被预测为稳定的那个结晶水。

图7. 蛋白相互作用势分析。左：湿的4EE0结合位点；右：干的4EE0结合位点。红色：正蛋白相互作用势等值面（ev=+8）；蓝色：正蛋白相互作用势等值面（ev=-8）

结果如图7所示，"干"（不包括结晶水分子）和"湿"（包含了活性位点稳定的结晶水分子）4EE0活性位点蛋白相互作用电势以令人满意的方式与化合物9的配体场匹配。尤其是苯并吡啶片段的含孤对电子氮的吡啶环位于4EE0活性中心的正相互作用势区域，显然，湿结合位点比干结合位点的红色正相互作用区更大，并覆盖了吡啶环。

图8展示了化合物9的配体场，吡啶氮方向上的负静电势与蛋白的正相互作用势互补，这是对化合物的活性是有利的。

图8. 化合物9的配体场。红色：正静电势；蓝色：负静电势

hH-PGDS抑制剂的SAR

图9. 化合物9、11、13、14的配体场比较

化合物9苯并吡啶的吡啶环氮方向上配体场的符号直接与化合物的活性相关。如图9所示的配体场，化合物9、11在吡啶环方向上有强的负配体场，与蛋白的正相互作用势匹配（见图10），9与11的IC₅₀分别为2.34与8.26nM。化合物13由于具有正电荷的氨基，在同样方向配体场全部为正，与蛋白的正相互作用势完全冲突（图图10）；此外13酰胺片段处还与蛋白在静电上发生冲突（见图10），因此活性最低，其IC₅₀为1480nM；化合物14的羟基具有部分负配体场与部分正配体场，在与蛋白的正相互作用势区有冲突的同时也有一定的匹配性（见图10），因此其活性强于化合物13而弱于化合物9与11，其IC₅₀为845nM。

图10. 化合物9、11、13、14的配体场（实心等值图）与4EE0“干”结合位点的蛋白相互作用势（网格线等值图）比较

上面的静电比较也存在一些问题。首先，四个化合物不仅只是在静电上有差异，而静电又不是影响活性的唯一因素。其次，四个化合物与同一个蛋白结构4EE0的“干”结合口袋进行比较虽然使静电比较变得简单、易行，但是没有考虑蛋白构象变化的因素。同一个结合口袋与不同配体结合时，结合口袋的残基侧链可能会因此而变化，哪怕仅仅是末端羟基的取向变化都可能导致蛋白的相互作用势变化。鉴于此，进一步进行了静电互补性打分，单独对每个化合物与各自结合口袋进行定量的静电比较。

总的来说，基于配体场与4EE0干结合位点的蛋白相互作用势似乎可以解释化合物9、11、13、14的SAR，这提醒进行水分子替换的分子设计需要同时考虑化合物与结合位点之间的静电特性。

hH-PGDS抑制剂的静电互补性（Electrostatic Complementarity）分析

Flare的静电互补性打分（Electrostatic Complementarity score，EC score）的分析方法^[10]将配体静电的分析与蛋白静电的分析结合起来，以产生配体-蛋白质复合物的静电互补性(EC)的视觉评估和数值评估。基本思想非常简单：当配体和受体的静电势匹配(即具有相同的量值和相反的符号)时，实现配体和受体之间的最大静电亲和力。通过将该方法应用于一系列文献数据集，我们分析了视觉和数值分量，表明它与报道的生物活性差异相关，并能够预测报道的生物活性差异^[10]。我们用Flare的EC技术对化合物9,11,13,14等进行了进行静电互补性比较。

图11. 化合物9、11、13、14的静电互补性打分（EC score）与pIC₅₀具有极其显著的相关性，Pearsion相关性系数R²=0.897。

XED力场的静电互补性分析表明，互补性EC打分值与IC₅₀的负对数(pIC₅₀)具有极其显著的相关性，Pearsion相关性系数R²=0.897。结合上面的蛋白相互作用势与配体场分析，静电相互作用可能是影响该系列化合物活性的主要因素。

当然，这不是真正严格意义上的QSAR分析。首先，数据点太少（才4个数据点），不具备严格的统计学意义。其次，数据点的活性分布并不均匀，也不满足QSAR的前提。如图11所示，数据主要分布在两端而缺少中间的点，此类数据“很容易”仅仅巧合而具有高的R²。

此外，配体与蛋白的结合取决于许多不同的物理效应。静电是其中之一，也是非常重要的一个，但它不是唯一的一个，所以EC分数只会预测其他效应不变时的活性差异。在本算例中，四个化合物之间的差异不仅仅只是在空间上某一点的静电差异，因此，EC与蛋白相互作用势在本算例中用于定量解释尚需更多的数据来建模并进行严格的验证，但并不影响水能量学特征以及构象能对活性产生影响的结论。

总结

本文从水的热力性质出发，用Flare GIST方法分析了Trujillo等人双水记的两个水分子，结果发现W1为ΔG小于0的happy water；而W2为ΔG大于0的unhappy water。根据一般的处理原则，应该将W1水做为蛋白的一部分设计新化合物与之相互作用或不要动它，或者用配体的一部分完美模拟W1同时付出踢走W1的代价；对于unhappy的W2，可以用替换或置换的策略。水分子的能量学特征与化合物9,11,13,14活性差异是一致的。

本文还在B3LYP/6-31G*理论水平下对关键两面角进行了分析，验证了作者关于过高内能的假设。水分子分析与基于QM的两面角特征分析解释了为什么对W1的替换设计化合物13与14活性大幅降低，而对W2置换设计的化合物11仅稍有降低。化合物13、14的关键两面角能量特征与化合物9,11,13,14活性差异也是一致的。

本文还进行了蛋白相互作用势与配体场的分析，结果表明化合物9、11、13、14的SAR与场互补性有直接的关系；进一步用XED力场对4个化合物进行了静电互补性分析，发现静电互补性打分值(EC score)与IC₅₀的负对数(pIC₅₀)具有极其显著的相关性，Pearsion相关性系数R²=0.897。需要注意的是，这种相关性存在数据点太少、活性数据分布不均匀等统计学的缺陷，相关性也许仅仅只是一种巧合。

总的来说，本文从水的能量学特征与构象稳定性解释了双水记分子的活性差异；并发现蛋白相互作用势与配体场以及静电互补性对活性差异也有极强的相关性，但这种相关性不具备严格的统计学意义，需要更多数据进行建模与验证。这也提醒了采用水分子替换策略进行分子设计的时候，需要同时考虑分子的构象能与静电特性。

算例下载

对PDB 4EE0进行GIST计算的Flare算例：PGDS-4EE0-GIST.flr

文献

1. 肖高铿.3D-RISM预测水分子的位置与稳定性.墨灵格的博客. http://blog.molcalx.com.cn/2017/06/18/flare-3d-rism-tutorial.html
2. 肖高铿.SPARK替换结晶水分子.墨灵格的博客. http://blog.molcalx.com.cn/2017/06/15/displacing-crystallographic-water-molecules-with-spark.html
3. Trujillo, J. I.; Kiefer, J. R.; Huang, W.; Day, J. E.; Moon, J.; Jerome, G. M.; Bono, C. P.; Kornmeier, C. M.; Williams, M. L.; Kuhn, C.; et al. Investigation of the Binding Pocket of Human Hematopoietic Prostaglandin (PG) D2 Synthase (HH-PGDS): A Tale of Two Waters. Bioorg. Med. Chem. Lett. 2012, 22 (11), 3795–3799. https://doi.org/10.1016/j.bmcl.2012.04.004.
4. Nguyen, C.; Gilson, M. K.; Young, T. Structure and Thermodynamics of Molecular Hydration via Grid Inhomogeneous Solvation Theory. 2011. arXiv:1108.4876v1. https://arxiv.org/abs/1108.4876
5. 黄传龙. replace水分子而活性降低的例子. 微信公众号：CADD谈药，2021年3月7日
6. Tanger, J. C.; Pitzer, K. S. Calculation of the Thermodynamic Properties of Aqueous Electrolytes to 1000.Degree.C and 5000 Bar from a Semicontinuum Model for Ion Hydration. J. Phys. Chem. 1989, 93 (12), 4941–4951. https://doi.org/10.1021/j100349a053.
7. Rashin, A. A.; Bukatin, M. A. Continuum-Based Calculations of Hydration Entropies and the Hydrophobic Effect. J. Phys. Chem. 1991, 95 (8), 2942–2944. https://doi.org/10.1021/j100161a002.
8. Grant, J. A.; Pickup, B. T.; Nicholls, A. A Smooth Permittivity Function for Poisson-Boltzmann Solvation Methods. J. Comput. Chem. 2001, 22 (6), 608–640. https://doi.org/10.1002/jcc.1032.
9. E. L. Mehler, The Lorentz-Debye-Sack theory and dielectric screening of electrostatic effects in proteins and nucleic acids, in Molecular Electrostatic Potentials: Concepts and Applications, Theoretical and Computational Chemistry Vol. 3, 1996
10. Bauer, M. R.; Mackey, M. D. Electrostatic Complementarity as a Fast and Effective Tool to Optimize Binding and Selectivity of Protein-Ligand Complexes. J. Med. Chem. 2019, acs.jmedchem.8b01925. https://doi.org/10.1021/acs.jmedchem.8b01925.

需要Flare、计算服务或亲自试用，请联系我们