CN101419214B - 基于分层原子加和模型的分子酸碱解离常数的预测方法 - Google Patents

Abstract

本发明涉及一种全新的基于分层原子加和模型的有机小分子酸碱解离常数(pK_a)的预测方法。基于在给定温度下酸碱解离平衡的自由能变化的线性关系规律，该方法首先以Hammett-Taft方程和Cherkasov等人提出的“解离中心-其余部分”的处理取代基效应建立一种分层的原子加和模型，然后通过该模型由化合物具体结构计算相应化合物的酸碱解离常数pK_a值。该方法不存在常规方法涉及大量取代基电子效应常数及校正因子的困难，同时保证了预测的快速、准确，对多个样本集表现出良好的数据拟合和预测能力，打破传统的在新药研发的早期对药物pK_a进行研究，大规模降低了药物研发成本，提高了新药发现效率。

Description

基于分层原子加和模型的分子酸碱解离常数的预测方法

技术领域

本发明涉及一种新型的基于分层的原子加和模型的分子酸碱解离常数预测方法，适用于根据化合物分子结构信息对该化合物进行虚拟体外高通量吸收、分布、代谢、排泄等(ADME)性质的评价和筛选。

背景技术

传统的药物设计往往首先关注化合物的活性(一般为体内活性)，然后才是化合物的药代动力学特征(图1)。但随着活性化合物发现能力的不断增强以及ADME研究的不断深入，这种传统的、效率低下的串行药物开发流程已经不能适应现代的药物开发需要。取而代之的更具效率的并行化药物开发流程(图1)。在并行化流程中，药物的多种性质同步优化，这就可以大大加速药物开发的效率。但是这种并行化流程在现在多数情况下仍是一个理论上的概念，主要原因是建立高通量药代动力学筛选模型的难度较大。细胞层次的研究经过多年的发展，已经建立了一些高通量的体外ADME研究方法，比如测定肠吸收的细胞单层转运试验、基于肝细胞或提取肝微粒体的新陈代谢和药物-药物相互作用试验、基于肝细胞或其他组织细胞的生长抑制为指标的细胞毒性试验等。但是到目前为止，高通量筛选还仅仅局限于少数几种药代动力学性质。在试验实现高通量筛选存在困难的情况下，发展有效的药代动力学特征理论预测及高通量虚拟筛选模型具有非常重要的意义。

药物分子的酸碱解离常数(pK_a)是显著影响药物生物利用度的性质。任何药物在起作用之前一般至少会通过被动扩散或主动转运跨生物膜一次，根据给药途径和作用靶标的位置差异，药物介入的溶液环境会有明显不同。对酸性或碱性药物来说，分布和扩散很大程度依赖于环境pH值下该药物分子的解离程度。类似的，受解离程度影响的药物分子的静电性质和氢键形成性能，药物分子与结合位点或主动转运载体之间的亲和性(结合自由能)也与环境pH值下的分子解离程度有很重要的关系。

另外，辛醇/水分配系数(P，logP)是决定肠吸收和药物可利用度的重要属性之一。然而logP是指分子中性状态下的脂水分配比例，对于酸或碱因为其分配依赖于水相的pH值，必须考虑其解离状态。这种pH依赖的分布系数(D，logD)是P和pK_a的函数。对于最简单的一元酸来说：

log D＝logP-log(1+K_a/[H⁺])≈logP-(pH-pK_a)    (eq.1)

在实际应用中，logD是比logP更重要的评价指标，而其依赖于pK_a的计算。因此，pK_a是早期药物设计和筛选，尤其是用于指导化合物定购和合成的一项重要ADME性质相关指标。

到目前为止已经有一些模型和方法可以对pK_a提供不同精度的预测，这些模型基本可以分为以下三种：

1)量化模型；

2)半经验模型；

3)基于化学信息学的模型。

大部分pK_a经验模型是基于线性自由能关系的，即可用著名的Hammett和Taft方程描述如下：

${\mathrm{pK}}_a=p{K}_a^0-\mathrm{\ρ\Σ\σ}---(\mathrm{eq}.2)$

$p{K}_a=p{K}_a^0-{\mathrm{\ρ}}^{*}\mathrm{\Σ}{\mathrm{\σ}}^{*}---(\mathrm{eq}.3)$

在给定温度下，酸碱解离平衡的自由能变化遵循线性关系，插入一个取代基将会引起的自由能变化ΔG和相应的pK_a变化。因此，对于指定分子的pK_a可以通过已知分子的pK_a ⁰和量化与其相关的取代基效应(如识别特征性碎片)来进行预测。

基于碎片的方法在定量构效关系(Quantitative Structure-ActivityRelationship，QSAR)研究中的重要性正在得到越来越广泛得认识，在许多商业程序中也使用基于碎片的方法进行pK_a的预测，如PALLAS/PKALC(CompuDrug Chemistry Ltd.)和ACD/pK_a-DB(Advanced Chemistry Development Inc.)。但是如前文提到的，基于碎片方法的计算很大程度上取决于预先定义的碎片常数库的大小，另外通常只能处理有公共结构的同系物，扩展性有限。

根据某一基团的解离程度取决于相邻的原子和键所定义的子环境这一假设，Xing等人提出了一种2D树型分子指纹的方法进行pK_a预测。在该模型中，每个解离中心以一个由出现频率和在解离中心周围每个原子和基团类型的位置所组成的计数矢量所表征。这种方法只定义了少数几个简单的基团，更类似于原子加和方法：同种原子或基团类型的贡献根据其距离解离中心的不同而不同，从而避免了直接使用大量的预定义碎片值。但是，这种方法事实上并没有降低模型复杂度，最终的取代基效应也只考虑了距解离中心五个连接键的原子和碎片贡献，计算精度较差。

Cherkasov等人提出一种非碎片3D方法量化取代基效应来进行pK_a预测。该方法的理论基础主要来自于Taft方程(eq.3)和诱导效应“通过空间传递”的假设。预定义参数不完全的问题通过采用原子分类和“解离中心-其余部分”的思路进行解决，提高了模型预测的可扩展性。然而由于计算需要分子的3D结构信息，而在模建过程中没有考虑分子柔性和多重构象问题，该模型仍有待进一步改进。

新药研发过程费用昂贵、时间冗长、淘汰率高。平均一个新药的研发需要花费10亿美元以上，耗时约12年时间。大约有90％的候选药物在临床期间被淘汰，这是研发过程费用昂贵的主要因素。候选药物淘汰的原因中很大一部分是药物动力学方面引起的，例如口服吸收性差，生物利用度低。传统的新药研发流程中，药物动力学的研究处于研发的中后期。然而，如果在药物研发的早期对药物动力学进行研究，以尽早淘汰药代动力学参数不理想的候选药物，可以大规模降低药物研发成本，提高新药发现效率。计算机模拟药物动力学正在为全球各大制药公司应用，并将会在药物的研发过程中起到越来越重要的作用。

发明内容

本发明的目的主要是提供一种适合高通量筛选，同时预测准确、适用范围广泛的有机分子酸碱解离常数评价方法。

该目的主要是通过采用一种新型的基于分层原子加和模型和遗传算法的参数估计方法，对现有分子结构和相关pK_a数据进行训练而达到的。

在本发明中，我们首先基于Hammett-Taft方程和Cherkasov等人提出的“解离中心-其余部分”的处理取代基效应的思路，建立一种分层的原子加和模型。其中任一原子对给定解离中心的电子效应贡献可以认为是该原子的原子类型及其与解离中心的距离决定的。我们定义这种效应贡献为该原子的内在电子效应常数(σ)和其所处的位置的电子效应传递性能(T)的函数：

${\mathrm{pK}}_a=p{K}_a^0-\mathrm{\ρ\Σ}\left(\mathrm{\σT}\right)---(\mathrm{eq}.4)$

与Hammett-Taft方程相同，pK_a ⁰是参考化合物的解离常数值，ρ是解离中心依赖的参数，反映了不同解离中心对取代基效应的敏感性差异。

为了描述T，我们还需要建立一套分子结构框架来表征原子所处的位置。这里，我们对每个分子都构建了如图2所示的以酸碱解离中心为根的分子连接树。分支的层数和原子编号分别记做i和j，则位置(i，j)表示处于第i层的第j个原子。然后引入描述符γ表示每个原子类型的电子效应传递能力。根据电子效应“通过键传递”的假设(键的连续极化)，我们提出下面的电子效应传递性公式用于T的计算：

$T_{i,j}=\underset{k}{\overset{\mathrm{pathes}}{\mathrm{\Σ}}}\underset{l=1}{\overset{i-1}{\mathrm{\Π}}}{\mathrm{\γ}}_{i,j,k,l}---(\mathrm{eq}.5)$

其中k表示给定位置(i，j)到解离中心(根)的第k条通路，l表示通路上每个原子的层编号。这里我们认为位置(i，j)的电子效应贡献是所有到根最短路径的累积贡献，所以只考虑内部(i-1)层的原子。

使用公式eq.5可以定量计算每个原子的电子效应贡献。

因此，通过给传统原子累加方法引入一个位置项T，我们可以将其用于pK_a的计算。

为了简化模型，我们注意到参数ρ可以进一步合并到传递项中。因为传递因子(γ)根据原子类型而不同，而原子分类规则又考虑到了周围原子的类型，如果我们定义第一层的原子(直接与解离中心相连)具有不同的类型，反应参数可以被隐式的包含在原子分类中。因此，pK_a的计算公式可以重写为：

$p{K}_a=p{K}_a^0-\underset{i}{\overset{\mathrm{levels}}{\mathrm{\Σ}}}\underset{i,j}{\overset{\mathrm{atoms}}{\mathrm{\Σ}}}\left({\mathrm{\σ}}_{i,j}\underset{i,j,k}{\overset{\mathrm{pathes}}{\mathrm{\Σ}}}\underset{l=1}{\overset{i-1}{\mathrm{\Π}}}{\mathrm{\γ}}_{i,j,k,l}\right)---(\mathrm{eq}.6)$

进一步概括来说，本发明一种新型的基于分层的原子加和模型的分子酸碱解离常数预测方法包括如下步骤：

a)根据待测化合物的结构确定分子的解离中心，并根据其类型确定其基准解离常数pK_a ⁰，建立以解离中心为根的分子连接树；

b)根据分子中任一原子的类型和所处位置确定该原子的内在电子效应常数σ和电子效应传递性能γ；

c)根据下式计算该化合物的酸碱解离常数pK_a：

$p{K}_a=p{K}_a^0-\underset{i}{\overset{\mathrm{levels}}{\mathrm{\Σ}}}\underset{i,j}{\overset{\mathrm{atoms}}{\mathrm{\Σ}}}\left({\mathrm{\σ}}_{i,j}\underset{i,j,k}{\overset{\mathrm{pathes}}{\mathrm{\Σ}}}\underset{l=1}{\overset{i-1}{\mathrm{\Π}}}{\mathrm{\γ}}_{i,j,k,l}\right)$

其中，k表示给定原子位置(i，j)到解离中心的第k条通路，l表示第k条通路上原子的层编号，给定原子位置(i，j)原子处于第i层的第j个原子，pK_a为酸碱解离常数，pK_a ⁰为参考化合物解离常数。

作为取样标准，选择25℃，水溶液中测得1300个化合物的pK_a实验值，数据样本来自于Lange′s Handbook of Chemistry，其中选取最常见的41种解离中心，其pK_a ⁰值由其对应的最简单的有机分子确定，如醇羟基-OH的pK_a ⁰取值为甲醇的pK_a，并使用SMARTS格式描述该解离中心，由此构建得到解离中心参数表，再根据待测化合物结构与解离中心表对比即可确定基准解离常数pK_a ⁰。

步骤b中内在电子效应常数参数表和电子效应传递参数表可通过如下途径得到：对于每种原子类型定义引入两个原子类型依赖的变量σ和γ，分别表征该原子对解离中心的内在电子效应贡献和电子效应传递性能，然后用SMARTS对原子类型定义分类规则；使用可编程原子分类PATTY回溯算法完成原子类型指定，得到内在电子效应常数参数表和电子效应传递参数表。

本发明所述的分子酸碱解离常数预测方法中涉及的计算模型对应作为算法依据，可以使用C++开发模块化的分子酸碱解离常数预测程序。该程序可以自动识别常见的小分子文件格式，无需小分子的3D结构信息，具有快速(～500个化合物/秒，P4 2.4G)、稳定的特点。程序用C++完成，在多个平台编译测试通过，包括：Visual C++7.0(windows)、MIPS C++(SGI)、GCC(Linux)。可以快速、准确的预测有机分子的酸碱解离常数。

本发明方法中每个预定义参数值使用基于遗传算法(GA)的适应性参数优化进行拟合。该模型不考虑小分子的三维信息，同时也无需进行量化和分子机制等计算，可以快速准确的预测有机小分子的pK_a值。既可以保持基团加和方法数据拟合能力强，模型物理意义强的特点，同时可以较好的避免现有方法碎片类型预定义不足，扩展性不强的问题：该项技术适用于药物研发早期阶段预测尚未合成的化合物的ADME性质，辅助研究人员选择那些ADME性质较好的化合物进行合成、药理测试，并能指导临床试验的设计，从而大大减少人力、物力的浪费，降低药物研发成本，提高药物研发效率。

附图说明

图1、传统药物研发和现代药物研发流程的比较；

图2、以解离中心为根的氯苯胺分子连接树的构建；

图3、训练集分子pK_a预测值与实验值的相关性散点图；

图4、测试集分子pK_a预测值与实验值的相关性散点图；

图5、对包含22个上市药物的基准pK_a测试集的统计结果；

具体实施方式

实施例1.参数表构建

最小二乘法是传统基团/原子加和模型最常用的拟合组成性参数的方法，但对于非线性模型，并不适合使用最小二乘方法进行模建。已经证实自适应性优化是一种有效的处理对不了解搜索空间(解空间)问题的方法。在本研究中，发明人使用遗传算法进行自适应的参数估计，具体步骤如下：

数据样本来自于Lange′s Handbook of Chemistry。作为取样标准，选择25℃，水溶液中测得的pK_a实验值共1300个化合物，涵盖了广泛的化学多样性空间。其中选取最常见的41种解离中心，其pK_a ⁰值由其对应的最简单有机分子来确定，如醇羟基-OH的pK_a ⁰取值为甲醇的pK_a；并使用SMARTS格式描述该解离中心，由此构建得到表1。

表1、解离中心参数表。

ID解离中心SMARTS描述                                      pK_a ⁰   解离中心名

1  [#1][OH，Oh][SD4](＝O)(＝O)                            -1.86  磺酸(sulfonic_acid)

2  [#1][Oh，OH][n+]                                       0.688  氧化吡啶(pyridine_oxide)

3  [#1][OH，Oh][SD3]＝O                                   1.5    亚磺酸(sulfinic_acid)

4  [#1][OH，Oh][PX4]([OH0，Oh0])(＝O)[OH，Oh]             1.54   磷酸_1(phosphoric_acid_1)

5  [#1][OH，Oh][As]＝O                                    1.67   次胂酸(arsinic_acid)

6  [#1][NH3+，Nh3+，Nh2+，Nh2+，Nh+，Nh+][CX3](＝O)[OH0]  1.7    氨基甲酸酯_1(carbamate_1)

7  [#1][OH，Oh]C(＝O)[N+0]C(＝O)[OH，Oh]                  2.15   亚氨基二乙酸_1

                                                                 (iminodiacetic_acid_1)

8  [#1][OH，Oh][$([PX4][#6])](＝O)[OH，Oh]                2.38   膦酸_1(phosphonic__acid_1)

9  [#1][OX2h，OX2H][PX3]([！#8])＝O                       3.08   次膦酸(phosphinic_acid)

10 [#1][OH，Oh][AsX4](＝O)[OH，Oh]                        3.41   胂酸_1(arsonic_acid_1)

11 [#1][$([C^3]([NX3+](＝O)[O-])[NX3+](＝O)[O-])]         5.21   二硝基甲烷(dinitromethane)

12 [#1][OH，Oh][$([CX3][#6])]＝[Se]                       3.75   硒基醋酸(selenoacetic_acid)

13 [#1][OH，Oh][C]＝O                                     3.77   羧酸(carboxylic_acid)

14 [#1][OH，Oh]C(＝O)[AsH，Ash，AsH2，Ash2]               4.67   砷基醋酸_1

                                                                 (arsenoacetic_acid_1)

15 [#1][OH，Oh][SeX4](＝O)(＝O)[$(Oc)]                    4.79    硒酸(selenic_acid)

16  [#1][OH，Oh]B([OH，Oh])[$(O[#6])]                      5.5   硼酸_1(boric_acid_1)

17  [#1][#7^2+，n+，#7X3+]                                 5.17  硝基芳香烃_1

                                                                 (aromatic_nitro_1)

18  [#1][$(C1C(＝O)CCCC1(＝O))]                            5.26  烯巳二酮(cyclohexanedione)

19  [#1][$([#7]1C(＝O)N＝NC1(＝O))]                        5.8   三唑烷二酮

                                                                 (triazolidine_dione)

20  [#1][OH，Oh][PX4]([OH0，Oh0])(＝O)([O-])               6.31  磷酸_2phosphoric_acid_2

21  [#1][AsH2，Ash2，AsH，Ash]C(＝O)[O-]                   7.68  砷基醋酸_2

                                                                 (arsenoacetic_acid_2)

22  [#1][OH，Oh][$([PX4][#6])](＝O)[O-]                    7.74  膦酸_2(phosphonic_acid_2)

23  [#1][N+0][CX3](＝O)[OH0]                               7.88  氨基甲酸酯_2(carbamate_2)

24  [#1][OH，Oh][PX4]([N+0])(＝O)[O-]                      8.03  氨基磷酸_2

                                                                 (phosphoramidic_acid_2)

25  [#1][OH，Oh][AsX4](＝O)[OH0-，Oh0-]                    8.18  胂酸_2(arsonic_acid_2)

26  [#1][OX2h，OX2H]O[#6^2]＝O                             8.2   peroxic_acid

27  [#1][OH，Oh]C＝CC＝O                                   8.24  二酮_1(dione_1)

28  [#1]C(C(＝O))C(＝O)                                    8.3   二酮_2(dione_2)

29  [#1][NH0+，NH+，Nh+，NH2+，Nh2+，NH3+，Nh3+，NH4+，Nh  9.2   charged_amine

    4+]

30  [#1][OH，Oh][N][CX2]＝O                                9.4   羟基乙酰胺

                                                                 (hydroxyacetamide)

31  [#1][nH+0，nh+0]                                       9.972 硝基芳香烃_2

                                                                 (aromatic_nitro_2)

32  [#1][OH，Oh][B][OH，Oh]                                10    硼酸(boronic_acid)

33  [#1][N+0](C＝O)C＝O                10     二酰亚胺(imide)

34  [#1]C[NX3+1](＝O)[O-]              10.12  硝基甲烷(Nitromethane)

35  [#1][NH2，Nh2]C#N                  10.27  胺腈(cyanamide)

36  [#1][OH，Oh]B([O-])[$(O[#6])]      10.5   硼酸_2(boric_acid_2)

37  [#1][SX2h，SX2H]                   10.7   硫醇(thiol)

38  [#1][OX2h，OX2H][$(O[#6])]         11.5   过氧化氢物(hydroperoxide)

39  [#1]ON＝C                          13     肟(oxime)

40  [#1][NH2，Nh2，NH1，Nh1]C＝O       15.1   酰胺(amide)

41  [#1][OX2h，OX2H]                   15.54  醇(alcohol)

然后在这些数据中随机抽取1029个作为测试集，剩余271个用作测试集。另外发明人还对文献中常用的pK_a基准测试集进行了预测验证。

对于这种对解搜索空间不确定的问题，发明人使用遗传算法(GA)进行自适应的参数估计，具体步骤如下：

1)、染色体编码

在该问题中一个染色体组代表待定的一组原子参数集。为了确保参数的取值具有物理意义，使用实数数组编码的等位基因染色体，其中每个基因对应一个原子参数(有预先设定可取值范围的实数)。

2)、目标函数

GA的目标函数决定了种群遗传进化的方向，在该问题中进化的目标是提高模型的数据拟合能力，如最大化回归系数R²。实际应用中我们发现使用下式可以加快收敛：

$J=R^2/{\mathrm{SE}}^2=(n-m-1)(\frac{1}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(Y_i^{\mathrm{pred}}-Y_i^{\mathrm{obsv}})}^2}-\frac{1}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(Y_i^{\mathrm{pred}}-\stackrel{\‾}{Y})}^2})---(\mathrm{eq}.7)$

其中n代表训练集的样本数，m代表染色体长度。Y_ipred是使用公式(eq.7)和可能解对应的原子参数值计算得到的预测pK_a值。这样，种群中每个染色体的适应性使用公式eq.7进行评价，较高的J值代表了较好的数据拟合能力。

3)、遗传操作

为了得到最优解集，GA使用选择、交叉和突变等操作完成进化。选择操作我们使用较常见的轮盘赌方法；交叉操作使用混合交叉——先根据父值的距离产生均匀分布，然后按此分布产生子值；产生的子值又经过突变操作，这里我们使用高斯突变，即新值是基于当前值的高斯分布进行取值的。后两种特殊目的的操作尤其适合实数解空间的搜索。

经基于遗传算法的自适应参数估计之后得到的内在电子效应常数(σ)和电子效应传递参数(γ)的参数估计结果如表2和表3所示。

表2、内在电子效应常数(σ)参数表。

ID  原子类型的SMARTS描述	原子类型名	γ
			1    [*]2    [#6]3    [$([C^3+0]([#6^3，#1])([#6^3，#1])[#6^3，#1])]4    [$([C^3+0]([#6^3，#1])([#6^3，#1])([#6^3，#1])[#6^3，#1])]	DuC3XC3X1C3	-5.264850.745150.21621-0.34526

5     [$([C^3+0][#8，#16，N^3+0])]6     [#6^2]7     [$([C^2]([#6^3，#1])([#6^3，#1])＝[#6])]8     [C^1+0]9     [c]10    [$([c][#6^3，#1])]11    [$([c][#8，#16，N^3+0])]12    [$(c[N+])]13    [$([#8^3]([#6^3，#1])[#6^3，#1])]14    [#8^2]15    [#7^2]16    [$([#7^2]([#6])([#6^3，#1])[#6^3，#1])]17    [$([#7^2]C＝O)]18    [$([#7](＝[#6])([#6^3，#1])[#6^3，#1])]19    [N^1]20    [n]21    [$(n1naaa1)]22    [#7^2+]23    [$([#7](～[OX1])～[OX1])]24    [$([#16](～[OX1])～[OX1])]25    [#9]26    [#17]27    [Br]28    [I]

C3_OC2XC2C1CarXCarCar_OCar_N+Osp3Osp2XN2XNp1NamN2N1NarNar5_Nar5N2+NO2SO2FClBrI

-1.963331.62149-1.522790.583510.253620.53203-1.768871.299836.379466.058340.142890.96074-0.620634.946826.45076-0.1327221.933993.362526.818513.773794.974266.389436.955386.0952

表3、电子效应传递参数(γ)参数表

ID  原子类型的SMARTS描述	原子类型名	γ
			1    [*]2    [#6]3    [$([C^3+0][F，Cl，Br，I，$(N(～[OX1])～[OX1])])]4    [$([C^3+0][O，#16X2])]5    [$([C^3+0][N^3！H0+])]6    [$([C^3+0][$(C(＝O)O)，$(S(＝O)＝O)，$(P(＝O)＝O)])]7    [$([C^3+0]([N^3！H0+])[$(C(＝O)O)，$(S(＝O)＝O)，$(P(＝O)＝O)])]8    [$([C^3+0][#7^3+0])]9    [#6^2]10   [$([#6^2][#8])]11   [$([#6]＝[#8])]12   [$([#6](＝[#7+])[#7+0])]13   [$([#6](＝[#7+])[#7+])]14   [c]15   [$(c(c)(c)[#6^3，#1])]16   [$(cn)]17   [$(c1aaaa1)]18   [$(c1c[$(c[F，Cl，Br，I，C！^3，N+，$(N[C，#7，#8，#16]＝，#*)])，n]ccc1)]19   [$(c1cc[$(c[F,Cl，Br，I，C！^3，N+，$(N[C，#7，#8，#16]＝，#*)])，n]cc1)]20   [$(c1[$(c[F，Cl，Br，I，C！^3，N+，$(N[C，#7，#8，#16]＝，#*)])，n]cccc1)]21   [$(c(c)(c)[F，Cl，Br，I，#6！^3，N+，$(N[C，#7，#8，#16]＝，#*)])]	DuC3XC3_XC3)OC3_N4C3_COOCaC_NHC2XC2OCX＝OC(＝N+)NC(＝N+)N+CarXCarCar5_NarCar6_Narc6_7c6_8c6_6c6_5	0.08104510.4187540.3824630.4868010.5572180.5225170.6578921.591730.5033950.3201910.3710781.958280.5648050.6662770.4686490.6066460.7742630.5288460.6953650.5491210.564856

22  [$(c(c)(c)[#8])]23  [$(c(c)(c)[#16X2])]24  [$(c[NX4+])]25  [$(c1([NX4+])aaaaa1)]26  [$(c([n+]1)aaaa1)]27  [$(c([#1，#6^3])([n+]1)aaaal)]28  [$(c([F，Cl，Br，I，C！^3，N+，$(N[C，#7，#8，#16]＝，#*)])([n+]1)aaa[a]1)]29  [$(c([#8，#16，$([#7]([#1，#6^3])[#1，#6^3])])([n+]1)aaa[a]1)]30  [$(c([n+]1)aaa1)]31  [$(c([n+]1)([F，Cl，Br，I，#6！^3，N+，$(N[C，#7，#8，#16]＝，#*)])aaa1)]32  [$(c([n+]1)aa[a！c]1)]33  [$(cS(～O)(～O)O)]34  [$(cP(～O)(～O)O)]35  [#6^1]36  [#7]37  [N^2]38  [$(N[#6，#7，#8]＝，：，#*)]39  [$([N]([#6^3，#1])([#6^3，#1])[#6，#7，#8]＝，：，#*)]40  [n]41  [$(n1aaaaa1)]42  [#8]43  [#8^2]44  [#15]45  [#16]

Car_O3Car_S3Car5_N4Car6X_N4c6_1c6_2c6_3c6_4c5_1c5_2c5_3Car_SO3Car_PO3C1N3N2NplXNpln5_1n6_1Osp3XOsp2XPS

0.4062390.821540.2233890.7936490.8139850.7809960.680320.6078541.016710.9231720.5056041.094630.6179961.399030.396980.6845320.5690940.8850870.8293630.653550.09570510.1068750.2674610.468649

46    [$([#16](～O)～O)]

SO2

0.508219

对于训练集1029个化合物进行统计分析，结果显示本发明涉及的模型具有良好的统计学指标(R²＝0.9534，SE＝0.6504，图3)；同时对外部数据集(271个化合物)可以进行准确的预测(R²＝0.9331，SE＝0.7416，图4)。

为了进一步测试我们的模型对药物分子pK_a值的预测能力，发明人对经常作为pK_a模型评价的一个基准测试集进行验证(由22个药物分子组成，涵盖多种酸碱类型)，统计结果及化合物编号见表4，预测值与实验值的相关性如图5所示。

表4.对22个药物组成的测试集预测结果统计

ID药物名称	预测pKa	实验pKa	预测偏差
				1阿托品(Atropine)2氯霉素(Chloramphenicol)3氯噻嗪(Chlorothiazide)4氯丙嗪(Chlorpromazine)5西咪替丁(Cimetidine)6地西泮(Diazepam)7地尔硫卓(Diltiazem)8苯海拉明(Diphenhydramine)9丙吡胺(Disopyramide)10氟芬那酸(Flufenamic_Acid)11呋塞米(Furosemide)12氟哌啶醇(Haloperidol)	10.087511.65298.940049.917786.310284.03939.008689.4134510.46783.915873.453519.27588	9.911.039.59.36.83.38.91910.43.93.98.3	0.187490.62285-0.559960.61778-0.489720.73930.098680.413450.067790.01587-0.446490.97589

13丙米醇(Imipramine)14利多卡因(Lidocaine)15苯巴比妥(Phenobarbital)16苯妥英(Phenytoin)17普鲁卡因胺(Procainamide)18普罗帕酮(Propafenone)19普萘洛尔(Propranolol)20丁卡因(Tetracaine)21甲氧苄啶(Trimethoprim)22维拉帕米(Verapamil)

9.954477.795757.383419.215329.16569.047999.199378.990336.153459.62656

9.57.947.448.39.49.39.58.497.29.04

0.45447-0.14425-0.056590.91532-0.2344-0.25201-0.300630.50033-1.046550.58656

上表中药物1～22在图5中的数值点分布依次表示为O¹～O₂₂。

可以看出，对结构多样的药物分子，我们提出的分层原子累加的模型能保持很好的预测效果(R²＝0.9412，SE＝0.5052)，能够应用于pK_a的计算。

本发明中采用的分层原子累加的方法通过引入位置因子γ使得电子效应可以使用碎片加和贡献的方法进行计算，既保持了传统基团加和方法数据拟合能力强，模型物理意义强的特点，同时可以较好的避免现有方法碎片类型预定义不足，扩展性不强的问题。该方法在药物分子pKa预测方面取得了良好的应用。

实施例2

以对氯苯胺为例(图2)：

a、首先识别解离中心，使用基于子结构匹配的方法查询得到NH₃+，对应于表1中解离中心苯胺(ID：17)，pK_a ⁰值为5.17。构建如图2所示分子连接树，解离中心NH₃+为第0层(i＝0)。

b、根据表2、3中定义的电子效应常数(σ)表和电子效应传递参数(γ)表，使用基于子结构匹配的方法查询得到各层原子类型对应的σ和γ值。

c、以计算第5层(i＝5)原子对解离中心的电子效应贡献，首先确定处于第5层的原子(共1个，原子类型为氯，此处i＝5，j＝1)；查表2得到该原子的电子效应常数σ_5，1＝6.39(ID：26)；该原子涉及两条通路(k＝1，2)，通路上的各个原子查表3得到其传递参数，如通路1(k＝1，图2靠左的箭头标注)上的四个芳香碳(按所在层数标记，l＝1，2，3，4)对应的传递参数分别为：γ_{5，1，1，1}＝0.22(ID：24)，γ_{5，1，1，2}＝0.55(ID：20)，γ_{5，1，1，3}＝0.55(ID：20)，γ_{5，1，1，4}＝0.56(ID：21)。通路2上的电子效应传递性计算同通路1，则氯原子对解离中心的的电子效应贡献可由下式计算得到：

${\mathrm{\σ}}_{\mathrm{5,1}}{T}_{\mathrm{5,1}}={\mathrm{\σ}}_{\mathrm{5,1}}\underset{k=1}{\overset{2}{\mathrm{\Σ}}}\underset{l=1}{\overset{4}{\mathrm{\Π}}}{\mathrm{\γ}}_{\mathrm{5,1},k,l}={\mathrm{\σ}}_{\mathrm{5,1}}({\mathrm{\γ}}_{\mathrm{5,1,1,1}}{\mathrm{\γ}}_{\mathrm{5,1,1,2}}{\mathrm{\γ}}_{\mathrm{5,1,1,3}}{\mathrm{\γ}}_{\mathrm{5,1,1,4}}+{\mathrm{\γ}}_{5,1,2,1}{\mathrm{\γ}}_{\mathrm{5,1},2,2}{\mathrm{\γ}}_{\mathrm{5,1},2,3}{\mathrm{\γ}}_{\mathrm{5,1},2,4})$

$=6.39\×(0.56\×0.55\×0.55\×0.22+0.56\×0.55\×0.55\×0.22)$

$=0.49$

d、同理，计算其余各层原子的电子效应贡献，利用分层加和模型公式eq.7计算最后的pK_a值为4.14(实测值3.99)。

Claims (2)

1.一种基于分层的原子加和模型预测化合物分子酸碱解离常数的方法，其特征在于，包括如下步骤：

a)根据待测化合物的结构确定分子的解离中心，并根据待测化合物类型确定其参考化合物解离常数pK_a ⁰，建立以解离中心为根的分子连接树；其中参考化合物解离常数pK_a ⁰通过如下途径确定：

作为取样标准，选择25℃，水溶液中测得1300个化合物的pK_a实验值，数据样本来自于Lange′s Handbook of Chemistry，其中选取最常见的41种解离中心，其pK_a ⁰值由其对应的最简单的有机分子确定，并使用SMARTS格式描述该解离中心，由此构建得到解离中心参数表，再根据待测化合物结构与解离中心表对比即可确定参考化合物解离常数pK_a ⁰；

b)通过如下途径根据分子中任一原子的类型和所处位置确定该原子的内在电子效应常数σ和电子效应传递性能γ：对于每种原子类型定义引入两个原子类型依赖的变量σ和γ，分别表征该原子对解离中心的内在电子效应贡献和电子效应传递性能，然后用SMARTS对原子类型定义分类规则；使用可编程原子分类PATTY回溯算法完成原子类型指定，得到内在电子效应常数参数表和电子效应传递参数表；

c)根据下式计算该化合物的酸碱解离常数pK_a：

${\mathrm{pK}}_a={\mathrm{pK}}_a^0-\underset{i}{\overset{\mathrm{levels}}{\mathrm{\Σ}}}\underset{i,j}{\overset{\mathrm{atoms}}{\mathrm{\Σ}}}\left({\mathrm{\σ}}_{i,j}\underset{i,j,k}{\overset{\mathrm{pathes}}{\mathrm{\Σ}}}\underset{l=1}{\overset{i-1}{\mathrm{\Π}}}{\mathrm{\γ}}_{i,j,k,l}\right)---(\mathrm{eq}.6)$

其中，k表示给定原子位置(i，j)到解离中心的第k条通路，l表示第k条通路上原子的层编号，给定原子位置(i，j)表示原子处于第i层的第j个原子，pK_a为酸碱解离常数，pK_a ⁰为参考化合物解离常数。

2.根据权利要求1所述的分子酸碱解离常数预测方法，其特征在于步骤c所用公式eq.6由如下方法得到：

1)首先基于Hammett-Taft方程和Cherkasov等人提出的“解离中心-其余部分”的处理取代基效应的思路，建立一种分层的原子加和模型，得到：

$p{K}_a=p{K}_a^0-\mathrm{\ρ\Σ}\left(\mathrm{\σT}\right)---(\mathrm{eq}.4);$

其中ρ是解离中心依赖的参数；

2)为了描述T，建立一套分子结构框架来表征原子所处的位置；具体为：首先对每个分子都构建了以酸碱解离中心为根的分子连接树，分支的层数和原子编号分别记做i和j，则位置(i，j)表示处于第i层的第j个原子；然后引入描述符γ表示每个原子类型的电子效应传递能力；最后根据电子效应“通过键传递”的假设，得到如下电子效应传递性公式用于T的计算：

$T_{i,j}=\underset{k}{\overset{\mathrm{pathes}}{\mathrm{\Σ}}}\underset{l=1}{\overset{i-1}{\mathrm{\Π}}}{\mathrm{\γ}}_{i,j,k,l}---(\mathrm{eq}.5),$

其中(i，j)、k、l定义如权利要求1所述；

3)由eq.4和eq.5合并得到eq.6

${\mathrm{pK}}_a={\mathrm{pK}}_a^0-\underset{i}{\overset{\mathrm{levels}}{\mathrm{\Σ}}}\underset{i,j}{\overset{\mathrm{atoms}}{\mathrm{\Σ}}}\left({\mathrm{\σ}}_{i,j}\underset{i,j,k}{\overset{\mathrm{pathes}}{\mathrm{\Σ}}}\underset{l=1}{\overset{i-1}{\mathrm{\Π}}}{\mathrm{\γ}}_{i,j,k,l}\right).$

Images