水分子替换导致BTK抑制剂活性下降100多倍的理论解释

摘要：本文介绍了OpenEye SZMAP计算水分子位置及其热力学性质的原理，还介绍了水-中性水探针自由能差(n_ddG)的概念及其用途。本文还以一对结构相似的BTK抑制剂为例，演示了如何用SZMAP在BTK与其抑制剂结合位点里水的热力学性质、以及如何在结合位点里绘制n_ddG等值图、然后再用水的热力学性质ΔG;判断水是否适合被替换，并用n_ddG等值图来理解结合位点SAR，解释为何一个甲基的差异会导致活性110多倍的差异。

肖高铿/2021-08-02

BTK抑制剂算例简介

图1所示的化合物1与2是一对BTK抑制剂^[1]，其中1对BTK的抑制活性是其甲基化衍生物2的112倍，IC₅₀分别为1.3与145nM。

。

图1. BTK抑制剂化合物1与2的化学结构及其活性^[1]

结构生物学数据表明，除了与HOH843的相互作用之外，抑制剂1、2在与BTK的相互作用模式几乎一模一样：1与结合位点的HOH843发生氢键相互作用（PDB 6AUA），而化合物2的甲基则置换了该水分子（PDB 6EP9），如图2所示。

图2. 化合物1、2与BTK相互作用比较示意图。左：化合物1与BTK的结合模式（PDB 6AUA）；右：化合物2与BTK的结合模式（PDB 6EP9）；左边白色圆圈的水为HOH843，在6AUA中与化合物1发生氢键相互作用，在6EP9中被化合物2的甲基置换而消失。

将PDB 6AUA与6EP9按Cα进行重叠，发现化合物1与2基本完全重合，所以这一对分子在构象能上也没有差异，其生物学活性活性差异可以归结为对水分子HOH843的置换带来的，这为水分子置换计算解释SAR提供了非常好的算例。

图3. PDB 6AUA的水分子HOH 843通过氢键桥连着配体与蛋白相互作用、介导着氢键网络

图3展示了PDB 6AUA的HOH 843在化合物1结合位点里氢键网络，它通过氢键桥连着配体与残基THR474以及GLU475，并与HOH892、871形成水的氢键网络。形成鲜明对比的是，在化合物2-BTK的共晶复合物（PDB 6EP9）里该水消失，其它的氢键网络水却还保留着，化合物2的活性却下降为1的100多倍。因此，理论预见HOH 843的重要性具有有非常重要的药化指导作用。

本文的一个主要目的是，用SZMAP分析PDB 6AUA的结合位点，从蛋白结合位点的“SAR”角度来说明化合物1、2对BTK的活性差异。

SZMAP的方法简介

SZMAP (solvent-Zap-mapping)^[2]是一种半显式溶剂映射方法，可用来研究蛋白-配体复合物中水的热力学性质，它兼具了显式与隐式溶剂模型方法的优势。如图4所示，SZMAP将一个显式水探针放入到蛋白-配体体系里，水探针是一个三原子组成、与许多全原子力场的水模型相似，具有偏电荷（partial charges）。而其余的溶剂用ZAP toolkit的Poisson-Boltzmann隐式溶剂模型表示。SZMAP可对放置于蛋白上网格的每个点采集水的能量，从而在结合位点里构建水的热力学结合参数的三维映射图；当然也可以仅对用户指定的坐标点进行计算。比起显式的溶剂方法，SZMAP显而易见的一个优势是计算速度，因为它避免了为了获得可靠结果而进行长时间的平衡与高强度采样。

图4. SZMAP的半连续溶剂化模型：在蛋白-配体研究体系里放置显式水分子探针，并用隐式水填充其它位置，SZMAP可以在格点上或用户指定的位置采集水的取向。

SZMAP既可在每个格点上、也可以在指定的位置上放置探针计算热力学性质。在给定的探针位置，SZMAP采集探针水的各种取向。对于每一个取向，计算体系的能量(E_j)，包含下面几个成分：水探针与蛋白以及配体的vdW相互作用(E_{vdw protein-water})，水探针与配体、蛋白的静电相互作用以及极化的PB连续介质（溶剂）的相互作用（E_{PB protein-water}），还有蛋白、水的静电去溶剂化项(E_psolv与E_wsolv)。

$$
E_{j} = E_{vdw\ protein-water} + E_{BP\ protein-water} + E_{psolv} + E_{wsolv} \cdots(1)
$$

水探针取向之系综及其能量可以用来构建配分函数Q：

$$
Q = \sum_{j=1}^{N_{orient}}e^{-\frac{E_{j}-E_{min}}{kT}} \cdots(2)
$$

从配分函数Q，SZMAP可以推导出探针处于采样的某一个方向时的概率Prob_j：

$$
prob_{j} = {e^-{E_j \over {kT}} \over Q} \cdots(3)
$$

概率与配分函数可用于SZMAP的热力学性质计算，比如与连续介质的水比较计算探针水取向的熵变（Eq 4），探针的自由能变化（∆G, Eq 5）以及焓变（Eq 6）:

$$
T\Delta{S} = (-kT\sum_{j=1}^{N_{orient}}prob_{j}ln{prob_{j}} – kTlnN_{orient}) \cdots(4)
$$

$$
\Delta{G} = -(kTlnQ-E_{min} – kTlnN_{orient}) \cdots(5)
$$

$$
\Delta{H} = \Delta{G} + T\Delta{S} \cdots(6)
$$

此外，SZMAP还可用来计算水探针相对于两个不同参比态的热力学性质（图5）。其中一个参比态为真空探针，是一个置换了连续介质溶剂的球形空腔，与水探针相似却不含任何原子，对vdW能量项没有任何贡献。这个参比态设计用来在同一个位置对水探针与真空空间进行比较。另一个参比态是中性水探针(neutral water probe)，它是对脂肪烃基的模拟。中性水探针与水探针一样由三原子组成，但是所有原子的偏电荷均为0。此外，对中性水探针进行与对水探针一样的取向采样。因此，当中性水探针的自由能减去水探针的自由能，就抵消了vdW项，那么能量差就反映了由于水分子偏电荷带来的结合能差异。

图5. SZMAP的参比态示意图。上：真空参比态；下：中性参比态。水探针的着色反映了氢氧原子的偏电荷，而中性参比态没有着色反映了这些原子是没有电荷的。

SZMAP的一个主要应用，如其名所暗示，是在蛋白结合位点里映射出水的偏好位置。SZMAP可以在蛋白-配体复合物的网格点上进行采样，计算的溶剂化自由能以及其它热力学性质项的等值面等经渲染后可用来可视化呈现结合位点内这些值的变化，SZMAP的水-中性探针自由能差（neutral probe free energy difference，以下简称n_ddG或n_ΔΔG）衡量了正常、带电的水探针比不带电荷而具有水形状探针的偏好程度。在结合位点里，如果n_ddG数值为负数则表示在亲水区有利，如果数值为正则表示在疏水区有利。如果为0或接近0，意味着对极性或非极性分子均无偏好。一般来说-0.5kcal/mol的等值图可以用来将水分子探针比非极性水探针显著有利的位置给识别出来。

用SZMAP分析BTK结合位点

鉴于我们感兴趣问题是在化合物1上引入一个甲基以置换结晶水HOH843的分子设计（化合物2）是否对活性有利，因此，用SZMAP对HOH843位置进行n_ΔΔG分析正是我们需要的：如果HOH843处的n_ΔΔG为负值，说明引入甲基对活性不利；如果HOH843处的n_ΔΔG为正值，说明引入甲基对活性有利。

首先从PDB上下载PDB 6AUA，然后用Spruce进行蛋白结构准备，取Chain A的Alternative A用于SZMAP计算stabilization格点：

szmap  -mpi_np 24 \
       -stabilization true\
       -prefix 6auaA_stbl \
       -p 6AUA_AaltA__DU__BXJ_A-704.oedu \
       -excludeDUsolvent true \
       -results_set max

其中，stabilization选项设置为true指示szmap进行稳定化格点计算；其中excludeDUsolvent选项设置为true指示szmap将体系里的水忽略。这一步需要15分钟左右，结果获得一个新文件6auaA_stbl.oeb.gz。该文件包括了我们需要的SZMAP稳定化格点以及各种热力学性质。

还可以将Apo结构里HOH843所在位置的热力学性质提取出来，比如n_ddG:

1
2
3

szmap_grid  -input_mol 6auaA_stbl.oeb.gz \
            -grid neut_diff_apo_free_energy_grid \
            -at_coords HOH843.pdb

结果如下：

note: Attaching values from 'neut_diff_apo_free_energy_grid' to atoms in szmap_grid_interp.oeb.gz tagged as 'szmap_neut_diff_free_energy'.
OpenEye SZMAP results from 'neut_diff_apo_free_energy_grid':
num     atom    x       y       z       mask    n_ddG
1.1      O      19.583  5.022   4.047   1       -2.370
1.2      H2     19.669  5.871   4.484   0       -0.276
1.3      H1     19.240  4.409   4.758   0       0.000
1.total 1       mask    >       0.800   1       -2.370
1.avg   1       mask    >       0.800   .       -2.370

我们可以发现，HOH843氧所在位置的n_ΔΔG = -2.370kcal/mol，因此引入甲基的设计是对活性不利的，这与实验结果一致，活性降低了大约110倍。

用VIDA进行可视化分析，重点是对蛋白Apo结构的结合位点进行“水-中性探针自由能差（n_ΔΔG）分析。这里，我绘制了两个等值图：若n_ΔΔG小于-1.5kcal/mol，则将等值面渲染为黄色，以示水结合有利区（或者是中性水代表的亲脂基团结合有不利区）；若n_ΔΔG大于0.8 kcal/mol，则将等值图渲染为紫红色，以示水结合不利区（或者是中性水代表的亲脂基团结合有利区）。

图6. SZMAP计算结果，红色球：结合位点里的结晶水；等值图用中性探针自由能差n_ΔΔG渲染：黄色为n_ΔΔG小于-1.5 kcal/mol等值面，紫红色为n_ΔΔG大于0.8 kcal/mol等值面。

对6AUA结合位点Apo结构进行的n_ΔΔG可视化分析结果如图6所示：黄色区域为n_ΔΔG小于-1.5kcal/mol，紫红色区域为n_ΔΔG大于0.8。与配体酰胺NH相互作用的HOH843位于黄色等值面里，为n_ΔΔG小于0的区域，HOH843的n_ΔΔG = -2.370kcal/mol。

中性探针自由能差(n_ΔΔG)提供了蛋白角度的SAR分析：黄色区域为引入亲水基团有利的区域，而紫色区域为引入亲脂基团有利区，因此可以从配体与蛋白的匹配性来看分子设计是否合理。化合物1的酰胺羰基氧落在n_ΔΔG小于0的黄色区域，说明化合物这部分结构特征对结合是有利的。我们会发现化合物1的很多亲脂片段与紫色网格重合，紫色区域意味着代表脂肪烃基的无电荷水有利区，这意味配体的亲脂性部分与蛋白结构特征匹配良好。同样的道理，用甲基取代酰胺氮上的氢，则亲脂性的甲基与黄色网格代表的极性基团有利区相冲突，因此化合物2的甲基与蛋白结构特征是不匹配、对活性不利的设计。

n__ΔΔG在HOH843位点处给出的SAR可被另一系列BTK抑制剂复合物结构所证明，比如PDB4ZLY与4Z3V。在这两个位点里，HOH843被极性的酰胺片段相当完美的替换。

如图7所示，作为分子量非常小的起始片段，也可通过该水的替换而具有非常高的配体效率，LE高达0.53。

图7. Christopher等人<sup>[4]</sup>的片段筛选苗头化合物（PDB 4ZLY）结构及其活性

图7. Christopher等人^[4]的片段筛选苗头化合物（PDB 4ZLY）结构及其活性

上述片段系列证明具有非常好的发展前景，在该片段苯环附近有一个子口袋，口袋里有可替换结晶水HOH954(PDB 4ZLZ)，如图8所示，该结晶水介导了蛋白-配体氢键相互作用。用SZMAP分析，可以发现结晶水HOH954的ΔG = 1.173lcal/mol,说明这是一个可被替换或置换的结晶水。进一步用中性水探针自由能差进行SAR分析，该结晶水的n_ddG = -1.611kcal/mol，这说明该结晶水适合用进行水分子替换设计。该结晶水附近也被黄色的n_ddG = -0.8 kcal/mol等值面包围，说明该区域亲水。

图8. PDB 4ZLZ结合位点n_ddG的SZMAP分析及HOH954热力学性质。n_ddG等值图：黄色为小于-0.8kcal/mol, 紫红色为大于0.8kcal/mol

在另一个博文中，我们进一步对HOH954进行水分子替换实验^[5]，可重现Christopher等人活性提高到IC₅₀=4nM的化合物，见图9。

图9. Christopher等人的BTK抑制剂复合物结构（PDB 4Z3V）结合模式示意图

分析其复合物晶体结构PDB 4Z3V的结合模式可发现对HOH843替换的结合模式保留不变（如图9所示），同时完美地对HOH954水分子进行替换，拥有完美的与PHE943以及GLY414的氢键网络相互作用（图10-右）：实现将4ZLZ的HOH954相互作用(图10-左)合并到新化合物里。这进一步证实了n_ΔΔG给出的SAR信息时可靠的，可用于药化工作的优化指导。

图10. PDB 4Z3V配体的一部分完美替换了PBD 4ZLZ水分子HOH954

总的来说，Christopher等人^[4]IC₅₀=3.5μM的高配体效率片段证明了n_ddG给出HOH843结合水可被亲水基团替换有利的SAR判断；对该系列片段进一步对子口袋进行水分子分析并对HOH954进行水分子替换的优化最终获得IC₅₀=4nM活性的先导化合物，这再次证明了n_ddG描述的SAR是可靠的。

水HOH843: 留着还是置换/替换?

虽然前一节回答了甲基取代对活性的影响，但是没有回答这个水是否可以被置换或替换的问题。这就涉及对Apo结构计算水结合自由能的问题。用szmap_grid将HOH843坐标处的热力学性质提取出来，比如提取结合自由能(ΔG)：

1
2
3

szmap_grid  -input_mol 6auaA_stbl.oeb.gz \
            -grid apo_free_energy_grid \
            -at_coords HOH843.pdb

会得到下面类似的输出：

note: Attaching values from 'apo_free_energy_grid' to atoms in szmap_grid_interp.oeb.gz tagged as 'szmap_free_energy'.
OpenEye SZMAP results from 'apo_free_energy_grid':
num     atom    x       y       z       mask    dG
1.1      O      19.583  5.022   4.047   1       2.337
1.2      H2     19.669  5.871   4.484   0       0.995
1.3      H1     19.240  4.409   4.758   0       0.000
1.total 1       mask    >       0.800   1       2.337
1.avg   1       mask    >       0.800   .       2.337

其中dG列的第一行给出了apo结构的水合自由能ΔG，同样也可以得到相应位置上的其它热力学性质ΔH与TΔS。

OpenEye SZMAP results from 'apo_enthalpy_grid':
num     atom    x       y       z       mask    dH
1.1      O      19.583  5.022   4.047   1       0.851
1.2      H2     19.669  5.871   4.484   0       0.438
1.3      H1     19.240  4.409   4.758   0       0.000
1.total 1       mask    >       0.800   1       0.851
1.avg   1       mask    >       0.800   .       0.851

OpenEye SZMAP results from 'apo_entropy_grid':
num     atom    x       y       z       mask    TdS
1.1      O      19.583  5.022   4.047   1       -1.486
1.2      H2     19.669  5.871   4.484   0       -0.557
1.3      H1     19.240  4.409   4.758   0       0.000
1.total 1       mask    >       0.800   1       -1.486
1.avg   1       mask    >       0.800   .       -1.486

我们关心水HOH843氧原子位置上的ΔG、ΔH、与TΔS，分别为2.337、0.851与-1.486kcal/mol。由于水合自由能ΔG大于0，这是个unhappy water，因此易于被配体替换或置换。可以进一步借助BROOD与SPARK等软件的水分子提换策略来设计新的化合物。

虽然从ΔG角度HOH843提供了被置换或替换的机会，但是不正确的置换或替换并不能带来预期的活性增强。如前所述，中性探针自由能差可以提供置换或替换的SAR指导。一方面HOH843位置n_ΔΔG小于0，化合物2的甲基在性质上与之不匹配，甲基置换是对活性不利的；另一方面，甲基不能模拟HOH843的相互作用。两个因素最终造成化合物2的活性下降为1的112倍。

上面是从Apo结构计算的自由能来考虑配体周围水的处理方式。如果当蛋白-配体复合物结构可获得时，还可以用稳定化自由能来觉得水的处理方式，请参见SZMAP | 水分子的稳定化自由能计算及其应用^[6]。

关于蛋白SAR的思考：与Watermap比较

在确定HOH843的水合自由能ΔG大于0 kcal/mol之后，接下来就要决定如何处理该水分子：用配体的一部分来模拟水的作用而替换它（也就是替换或replace），还是用配体的一个疏水基团把这个水顶替掉（也就是置换或displace）。这其实是从蛋白角度考察SAR的问题。有多种不同方法用来判断什么时候使用替换或置换策略，用ΔG与ΔH等热力学性质来判断、研究SAR是一种流行的方法。Cappel等人在“Watermap在药物发现的应用”一文中有详细描述：ΔG与ΔH均大于0的水合位点应该是疏水区，建议用甲基或卤素等原子对水分子进行置换^[3]。原文引用如下：

When a hydration site has both ΔG 》 0 and ΔH 》 0 it is generally possible to displace the water with a ligand to yield a gain in potency (if the ligand fits in the binding pocket without making unfavorable interactions). Hydration sites with this thermodynamic signature are typically in hydrophobic regions, and displacement can be readily obtained with a sterically complementary group such as a CH3 or a halogen.

比较有意思的是，在本算例中，如果将SZMAP计算的热力学性质（ΔG、ΔH与TΔS分别为2.337、0.851与-1.486kcal/mol）按照Watermap对水分子置换的标准来指导新化合物设计将得到甲基取代对活性有利的SAR建议。显然，在这个例子里，这个建议错了。当然，不同方法计算的热力学性质可能会完全不同，将watermap的判断标准迁移到SZMAP计算的热力学性质上可能并不恰当。

从蛋白角度研究SAR，我认为n_ddG（n_ΔΔG）是一个非常直观的量，如图5所示：

n_ΔΔG = ΔG^Water – ΔG^{uncharged water}

如前所述，水代表了亲水基团，而中性水代表了亲脂性基团，而它们差值n_ddG（n_ΔΔG）的正、负号直接反映了该位点是亲脂有利还是亲水有利。从这一点看，n_ddG比起Watermap的ΔG与ΔH符号判断法更加直观。此外，通过绘制n_ddG（n_ΔΔG）等值图（如图6所示），借助3D可视化也更加直观地理解SAR以指导设计。

总的来说，在这个算例中，n_ΔΔG给出正确的SAR建议，而Watermap的ΔG与ΔH符号判断法不适合直接借鉴过来用于SAR建议。此外，n_ΔΔG等值图直观地给出了分子设计的SAR建议。

致谢

感谢中国海洋大学徐锡明老师、OpenEye的Gunther与黄鹤博士提供的宝贵建议。

文献

Nittinger, E.; Gibbons, P.; Eigenbrot, C.; Davies, D. R.; Maurer, B.; Yu, C. L.; Kiefer, J. R.; Kuglstatter, A.; Murray, J.; Ortwine, D. F.; et al. Water Molecules in Protein–Ligand Interfaces. Evaluation of Software Tools and SAR Comparison. J. Comput. Aided. Mol. Des. 2019, 33 (3), 307–330. https://doi.org/10.1007/s10822-019-00187-y.
SZMAP (Version: 1.6.2.1). https://docs.eyesopen.com/applications/szmap
Cappel, D.; Sherman, W.; Beuming, T. Calculating Water Thermodynamics in the Binding Site of Proteins – Applications of WaterMap to Drug Discovery. Curr. Top. Med. Chem. 2017, 17 (23), 2586–2598. https://doi.org/10.2174/1568026617666170414141452.
Smith, C. R.; Dougan, D. R.; Komandla, M.; Kanouni, T.; Knight, B.; Lawson, J. D.; Sabat, M.; Taylor, E. R.; Vu, P.; Wyrick, C. Fragment-Based Discovery of a Small Molecule Inhibitor of Bruton’s Tyrosine Kinase. J. Med. Chem. 2015, 58 (14), 5437–5444. https://doi.org/10.1021/acs.jmedchem.5b00734.
肖高铿. SPARK替换结晶水分子. 墨灵格的博客. http://blog.molcalx.com.cn/2017/06/15/displacing-crystallographic-water-molecules-with-spark.html
肖高铿. 水分子的稳定化自由能计算及其应用. 墨灵格的博客. http://blog.molcalx.com.cn/2017/06/16/openeye-szmap-water.html