有机小分子化合物解离常数pKa的预测研究进展

  对于弱酸和弱碱化合物, 解离常数(pKa)是最重要的理化性质参数之一, 其决定化合物的溶解度、亲脂性、生物富集性、毒性以及药物分子的吸收、分布、代谢和排泄(ADME) 性质. 通过实验方法测定化合物水溶液中的pKa受到物质稳定性、仪器测定范围以及人力物力消耗等多方面的限制, 因此过去几十年间发展了大量的pKa预测方法. 本文以有机小分子化合物为研究对象, 回顾了20年来pKa预测的研究成果, 包括pKa实验数据的来源、质量、 测定方法, 重点介绍3类预测方法(线性自由能关系模型、定量结构-性质关系模型和第一性 原理方法), 并简单总结了常用的商业软件, 最后提出未来pKa预测研究要关注的问题.

  解离常数(pKa)是有机物很重要的性质, 决定化合物在介质中的存在形态, 进而决定其溶解度、 亲脂性、生物富集性以及毒性[1~7]. 对于药物分子, pKa 还会影响其药代动力学和生物化学性质[8]. 精确预测有机物的pKa值在环境化学、生物化学、药物化学以及药物开发等领域都有重要意义. 本文系统调研了 近20年来针对有机小分子化合物在水溶液中pKa值的理论预测的研究成果, 包括实验数据的来源与质量分 析、理论预测方法的进展以及商业软件的开发情况等, 对该领域的研究做全面梳理, 以期对从事相关领域研究的科技工作人员起到一定的帮助. 本文所讨论的很多预测方法并不局限于水溶液, 只要有足够的训练集和预测 集数据, 也能够适用于别的环境条件, 比如气相或者其他溶液中化合物解离常数的预测.

  其中, A 代表酸分子解离后形成的共轭碱. 解离常数定义为该化学反应的平衡常数, 假定相关物质的活 度可以用浓度数据代替且水分子浓度为常数, 则:

  因此, pKa为50%化合物发生解离时溶液的pH. pKa越小, 化合物酸性越强. 通常情况下, 水中弱酸 pKa的范围为2~12. 相应地, 碱的解离常数使用其 共轭酸的pKa来表示 :

  对于每个物种, 溶液中的自由能都能分解成 气相自由能与溶解自由能之和, Gaq(AH)=G (AH)+ Gs (AH), 则整个质子反应过程:

  式(6)~(8)提供了一种pKa的理论计算方式, 即在不同的量子化学水平上计算气相质子转移能G 与溶解自由能Gs , 进而得到Gaq与pKa. 在气相和水相中的质子转移过程与反应自由能和溶解自由能的关系如图1所示.

  图1 气相和水相反应中质子从弱酸分子AH转移到水分子上的热动力学循环过程

  如果化合物含有n 个可解离的官能团, 溶液体系中总共可形成2n 个微物种(考虑到每个官能团解离与 否, 共有2n 种组合)和n2n−1个micro-pKas(2个微物种之 间的平衡常数). 所有带有相同数目的未解离质子的 微物种组成了体系的宏观状态, 共有n+1个(分别有0, 1,···, n个未解离质子), 即产生n个macro-pKas. 图2以 西替利嗪为例展示了带有3个可解离官能团的化合物 的解离过程, 表1列出了micro-pKas的实验测定值. 当 n3, 假如没有额外的信息或假设条件, micro-pKas无法通过滴定曲线西替利嗪的解离图解. 西替利嗪有3个解离官能团, 因此有 23 =8 个微物种和 3×22 =12 个 micro-pKas, 以及 3+1=4 个宏观状态和 3 个 macro-pKas. h 代表结合的氢原子个数[8]

  到目前为止, 已有文值, 例如, Beilstein数据库整合了很多小分子化合物pKa实验测定的文献, 共包含约15 万个pKa数据, 能够最终靠MDL Crossfire Commandar获得[11]; Langes化学手册也是pKa实验数据的一个重要来源[12]. 还能够正常的使用一些商业软件, 比如ACD和 SPARC, 查询化合物的pKa实验值, 其中, SPARC仅 能提供与所查询的化学结构完全一致的物质pKa数据[13], 而ACD允许用户进行相似性定义, 除了能给出精确匹配的物质, 还能提供一些结构相似物质的 pKa值[14]. 表S1总结了最新的免费和商业pKa数据库 以及包含pKa数据的书籍[15~33].

  除了数据来源(很多数据库是商业性质的), 数据 质量是另外一个重要问题. 即使是公开发表的数据, 也不一定都可靠, 有一些时候甚至相互矛盾. 例如, 在 Beilstein数据库中搜索苯酚, 会得到近30个数据, 并 呈现2个峰值, 有6个数据在1.0左右, 另外20几个数据分布在10.0左右. 很显然, 1.0是错误数据, 是再 次对pKa的负对数[34].

  可接受的实验测定误差大概为0.5个pKa [35]. 除 了实验本身的误差或者错误以外, 数据集的质量问 题还来自于以下几个方面[8]: (1) pKa数据被分配给错 误的化合物, 导致这个问题的原因可能是化合物名字有歧义、不标准, 或者是印刷错误等; (2) 错误的数值, 比如上面苯酚的例子, 或者pKa数据的印刷错误、 使用Ka或pKb(碱)代替pKa等; (3) 对含多个可解离中 心的化合物, pKa分配错误; (4) 同1个物质数据重复 输入, 即使数据是完全相同的, 在训练集中重复出现 会导致该化合物的重要性提高, 如果同时出现在训 练集和预测集中, 会提高验证结果的统计学性能; (5) 使用预测值代替实验值; (6) 不一致的实验条件, 比 如温度或者溶剂; (7) 不同实验室或实验人员的误差, 例如对二氯苯二磺胺, 有不同的文献数据, pKa1= 8.24, pKa2=9.50[36] vs. pKa1=7.40, pKa2=8.60[37]. 基于上述原因, 在使用这些数据时, 特别是用于 开发预测模型时, 要进行前处理, 例如筛选实验测定 条件、不同来源的数值进行统计分析、比照系列物的 pKa值、查阅原始文献等, 这些处理都会在一定程度上减少或消除数据质量上的问题.

  现代分析化学的发展提供了多种pKa测定方法, 有些已经实现自动化操作. 常用的实验手段包括电位滴定法、电导滴定法、分光光度法、毛细管电泳法 (CE)、核磁共振法(NMR)、高效液相色谱法(HPLC), 这些方法也可以联合使用进行pKa的高通量筛选, 例 如CE与质谱结合、HPLC与质谱结合等. 各种方法的 测定原理及优缺点, 详见文献[34, 38].

  Lee和Crippen[34]分析了Beilstein数据库中的pKa, 发现近一半的数据未提及测定方法, 另外一半数据, 绝大多数是使用电位滴定法、分光光度法和电导测定 法得到的, 占总量的98%以上. 此外, 尽管该数据库 包含了近15万pKa数据, 但对于单解离中心的化合物, 去除重复、可疑数据之后, 可靠的pKa数据不足2000个[5]. 由此可见, 作为化合物最基础的理化性质参 数, pKa的实验数据还相当缺乏, 并且存在着准确性 不高等问题, 特别是对于新兴的环境污染物以及待 开发的药物而言, 更是如此. 因此, 发展性能良好的 预测方法尤为重要.

  建立良好的pKa预测模型面临以下挑战: (1) 化合物构象灵活性, 有些化合物分子存在同分异构体或能够形成分子内氢键, 会显著影响其pKa值[39,40]; (2) 含有多个解离官能团, 因为要考虑micro-pKas, 使得pKa的预测变得更加复杂; (3) 可获得的用于建 立模型的pKa实验值. 尽管有些物质, 比如苯酚或羧 酸已有很多预测模型, 但是往往需要预测pKa的化合 物(新兴污染物或药物)不在现有模型的应用域内, 因此对于这些物质仍然需要精确的实验测定值. 在过去的20年间, 开发了大量的pKa预测模型, 涉及不同的计算水平、描述符类别以及统计分析方法 等. 当然, 回顾所有的研究方法不太现实, 因此本文试图总结大部分重要的研究成果, 列于表S2中[41~85].

  很多重要的pKa预测[41~45]都是基于线性自由能 关系, 使用Hammett方程, 认为取代基对pKa的贡献 是可以叠加的[42]:

  其中, pKa0 是母体化合物的解离常数;  是常数, 与化合物的种类有关; i表示第i个取代基对母体化合 物解离常数的影响; m是取代基的个数.

  该方法认为母体化合物携带了大部分化学信息, 取代基的贡献通过i来表征, 因此对于特定分子的 pKa预测受到如下限制: (1) 必须知道其母体化合物 的pKa值; (2) 必须已知所有取代基团的i 值.

  LFER模型是早期pKa预测研究中最常采用的方 法, 典型案例是Multi-CASE[44]. 共包含2464个有机 酸, 其中训练集1848个, 预测集616个, 也是最早使 用如此大规模数据的预测研究. 所有化合物按照中心解离官能团分成22个子集, 取代基(分子亚结构)对 中心官能团解离的影响通过碎片系数来表征, 例如, 对羧酸子集, 共定义了30个碎片系数. 此外, 还结合使用了分配系数、水溶解度、分子量、电负性、Hückel 分子轨道电荷密度以及硬度等分子结构参数. 基于训练集的模型对预测集的r2 =0.91, 均方根误差 RMSE=0.774; 基于全部数据的模型对214个药物分子的预测, r2 =0.70, RMSE=1.44[44].

  尽管是最早发展的方法, LFER模型目前仍用于商业软件或免费软件包, 例如EPIK[46,47]和SPARC[45]. SPARC结合LFER和扰动分子轨道理论, 将取代基对 pKa的影响分为共振效应、静电效应、溶剂化作用和氢键作用4个方面, 解决了键和电子离域对pKa的影 响. 这样处理使得SPARC可以计算micro-pKas, 并进而得到macro-pKas和其他相关性质, 例如滴定曲线个解离 中心的化合物), 总共4338个pKa值, 跨越30个数量级, 得到了很好的预测结果, RMSE为0.37. SPARC方 法最主要的应用局限在于其反应中心和取代基数据库有限, 有些性质还需要实验测定, 限制了对于新化合物的应用.

  由此可见, LFERs预测方法成熟、性能稳定、应用范围广、计算速度较快, SPARC和ACD等软件都可以实现micro-pKas和macro-pKas的同时计算, 并具有良好的预测结果, 因此LFERs成为目前商业或免费软件包最常采用的计算方法. 然而, LFERs最主要的问题在于对新物质的预测应用受到其基础数据库所涵盖范围的限制. 如果新物质中含有数据库所不包含的分子碎片或者取代基, 常见的处理方法是使用已有的分子量最接近的同族原子或基团替代, 例如 SPARC没有对所有原子类型参数化, 实际应用中使用碳原子代替硅原子, 硫原子代替硒原子. 很显然, 这种处理会影响预测的准确度. 因此, 扩大数据集的涵盖范围、提高预测精度是以LFERs理论为主的预测方法的主要发展方向.

  pKa QSPR模型的基本假设是pKa值的大小由化 合物的分子结构决定, 而分子结构信息可以通过计算或实验测定而获得, 称之为分子描述符. 由于化合物结构和性质多样, 目前开发使用的描述符种类非 常多, 专著“Handbook of Molecular Descriptors”[86]列 出了1600余种不同的描述符, 使用这一些描述符进而 通过统计学或者机器学习方法建立pKa预测模型. 预 测能力取决于模型所能够提取的描述符与pKa之间线 性或非线性相关关系的程度. 很多QSPR模型组合使 用了不同种类的描述符以及不同的统计分析方法, 具体详见表S2.

  在众多小分子化合物pKa的QSPR模型中, 量子化学描述符因其具有清晰的物理化学意义、有助于进行机理解释而得到广泛应用. 由于计算量子化学描 述符要进行分子结构优化, 早期的、大数据集的研究 多使用半经验(semi-empirical)方法[39,49,50,65,87,88]. 随着计算机软硬件水平的发展, 计算能力不断提高, 从头算(ab inito)和密度泛函理论(DFT)也得到越来越多 的应用[55,57,58].

  超离域度和极化率是常使用的量子化学描述符 [39,40,50]. Tehan等 人 [39,50]使用这 2类描述符以及Coulson净原子电荷、最低未占据分子轨道能建立了417个有机酸(包括苯酚和羧酸, 并按解离基团类型、 能够形成分子内氢键以及芳环上取代基的位置分为6 个子集)和282个有机碱(包括苯胺、脂肪胺和杂环化 合物共9个子集)的pKa预测模型. 分子结构优化和频 率分析都是在半经验AM1水平上开展的. 对有机酸 各子集, r2 0.80, RMSE0.76. 但该系列模型受到应 用域的限制, 对外部测试集的预测性能显著下降, 即 使将模型重新拟合(包括原始数据集和外部测试集), 统计学性能仍显著低于原始模型[40]. Yu等人[51]将上述数据集拓展为697个有机酸和502个有机碱, 使用 AM1方法计算了反映解离官能团原子的能量和电荷 信息的3类局域分子结构描述符, 建立了5类酸和7类碱的QSPR预测模型, 并进一步分析了影响化合物分 子解离的电子结构特征.作为化合物分子亲电性指数[89]的改进, 基团性质[90](解离基团各原子亲电性或亲核性之和)在pKa预 测中也得到了很好的应用[55]. 对63个羧酸、取代苯酚、苯胺、有机膦酸和脂肪醇分子, 基团亲核性的倒 数与pKa之间存在良好的线].

  分子静电势(MEP)也是常用的量子化学描述 符[56~58,91~93]. 3种MEP: 空间最小静电势(Vmin)、表面最小静电势(VS,min)和表面最大静电势(VS,max)与pKa有 较强的相关关系[57,58]. 近年来, 运用密度泛函活性理 论(DFRT)[91]发展的MEP和价自然原子轨道(NAO)描 述符也被创造性的应用于有机化合物pKa的预测中. DFRT方法认为分子的酸性来自于可解离原子(H)的 相邻原子, 即分子解离过程中的所有变化都能被相 邻原子的电子性质所反映, 从该意义上讲, DFRT方 法属于局域方法(不反映分子整体电子性质). Liu 和 Pedersen[56]使用MEP并减去表征类别的参考值后, 得到了关于硫醇、磺酸、脂肪醇、羧酸、苯胺和胺类物 质的单一的pKa预测模型. 在此基础上, 该方法也成 功应用预测了有不同数目取代基的羧酸以及苯酚的 pKa值[92,93]. 此外, Brinck等人[93,95,96]还使用可解离原子的最小表面解离能(IS,min)建立了QSPR模型. 尽管 MEP和IS,min都与pKa有较好的相关关系, 但他们是两 类不同的描述符, 而且通常来讲, VS,min和IS,min也不相 关[58,94]. 如果仅使用单一描述符, 则IS,min得到的模型 性能比使用MEP更好一些, 此外, 两类描述符结合 使用并不能提高模型性能[58].

  基于量子化学描述符所发展的QSPR模型最主要 的特点是模型透明、易于机理解释, 描述符个数一般 较少, 并且模型的建立也不会受到训练集化合物结 构差异性的影响. 很多量子化学描述符属于分子局 域性质参数, 可以计算micro-pKas. 模型的计算精确 度和时间取决于所使用的理论计算方法. 一般而言, 使用ab initio和DFT方法比半经验方法得到的结果准 确性和精确度更高, 但同时也意味着更多的计算资 源的消耗. 此外, 随着半经验方法理论的不断发展, 其计算性能越来越高, 因此预测较大分子或大数据 集化合物时, 要综合考虑计算时间和精确度, 选择合 适的理论水平和预测模型.

  除了量子化学描述符, 化学信息学方法也对pKa 预测研究做出了重要的贡献[59,60,63,97]. Xing等人[59,60] 认为官能团的解离依赖于其拓扑微环境, 使用树状 结构指纹谱图定量预测了有机化合物的pKa. 该研究 将含不同解离官能团的化合物分类, 使用50~300个 指纹图谱建立了预测模型 , 标准偏差 SE分别为 0.41(625个酸)和0.30(412个碱) [60]. 2007年, Jelfs等人[63]结合使用二维(2D)分子结构描述符和量子化学 参数, 并扩大了数据集范围, 得到了新的pKa预测模 型, 发现使用分子树(环形分子谱图)至少要考虑中心 原子的五级分子键.

  化学信息学方法计算速度快, 所包含的描述符 个数通常较多, 既有定量, 也有定性的描述符, 特别 适合处理化学结构和性质差距较大的大数据集化合 物. 其最主要的问题在于对含有多解离官能团的化 合物处理不好, 另外模型的稳定性和精确度也有待 提高.

  比较分子立场分析(CoMFA)方法也用于小数据集、结构相似化合物的pKa预测, 并使用偏最小二乘 方 法 (PLS) 建立了 4 个参数与 pKa 的线]. 需要注意的是, 尽管这些三维(3D)模型统 计学性能较好, 但是其基于2个前提: 训练集化合物 很少且化学结构高度相似以及现有结果依赖于所选 择的分子构象. 比较分子表面积分析(CoMSA)也使 用相似的3D描述符和PLS拟合了苯甲酸的pKa值[67]. 相比较而言, CoMFA方法对pKa的预测性能更好.

  与其他QSPRs模型不同的是, 3D-QSPRs仅能给出 影响解离过程的分子作用力场的信息, 而没有明确的 模型形式, 因此, 不利于深入的机理分析和预测应用.

  建立QSPR模型, 除了要选择合适的分子结构描 述符, 还要考虑统计分析方法. PLS[60,68,69]、多元线性 回归分析(MLR)[39,40,51]和主成分分析(PCA)[98]都是常 用的建模方法. 由于回归模型能够得到描述符和所 研究性质(pKa)之间明确的相关关系, 因此有利于机 理解释. 然而如果pKa与描述符非线性相关, 人工神 经网络(ANN)[99,100]、支持向量机(SVM)[101]等机器学习方法会更适用. 这一些方法, 更适合处理大量数据, 鉴别线性回归所不能揭示的非线性相关关系. HabiliYangjeh等人[70]使用PCA/遗传算法(GA)/ANN计算了 水溶液中282个含氮有机化合物的pKa值. 训练集、交 叉验证集和预测集分别包括170, 56和56个化合物, 使用406个描述符构建了179个主成分, 能够揭示 99.9%以上的pKa变化. 最终模型包含了15个主成分, 预测集的RMSE=0.0750 (ANN) vs 1.4863 (MLR). 该 模型最大的问题在于可能存在过度拟合以及无法鉴别 解离点, 即不能计算micro-pKas. ANN方法最成功的应用是Simulations Plus Inc.开发的ADMET Predictor[102], 训练集和预测集分别包括9075和2253个pKa值. 这是纯粹使用神经网络方法所建立的模型, 也是目前为数不多的能够计算micro-pKas的稳定性很高的应用程 序. 但该程序仅适用于含C, N, O, S, P, F, Cl和Br原子 的化合物. 机器学习的方法也越来越多地应用于 pKas预测[103~105]. Rupp等人[106]使用基于核的学习方 法预测了698个化合物的pKa值, 由于不需要分子结 构优化, 该方法运算速度较快.

  QSPRs模型得到广泛研究和应用的一个重要原 因在于大部分预测模型是简单、易于解释机理的线可以看到, 目前大多数QSPRs模型是基 于较少的训练集化合物或者按照解离官能团分类而 建立的, 化合物结构和性质差异较小. 这样的模型往 往具有较好的稳定性, 然而应用域也较小, 对域外化 合物的预测可能会产生很大误差, 这也是QSPRs模 型面临的主要问题. 由于pKa是非常复杂的分子性质, 针对化合物结构差异巨大的数据集建立单一的、性能优越的预测模型是不太现实的. 因此, 推荐的处理方法是先分再合, 即将整个数据集按照解离官能团先分类, 然后使用同一类描述符分别建立不同官能团 化合物的预测模型, 再将预测结果合并.

  此外, 在模型建立过程中要注意过度拟合的问 题, 特别是在含有非线性描述符时. 过拟合会导致模型的外部验证和预测应用能力显著下降. 描述符个 数过多会增加过度拟合的风险, 一般推荐的描述符 与化合物个数之比要小于1:5.

  使用热动力学循环(图1)或者直接在溶液中计算 能量是获得pKa值的另外一类常用方法. 该方法虽然 需要耗费大量的计算机资源, 但广泛适用于多种化 合物. 例如, Schüürmann等人[75,76,107]使用第一性原理计算了羧酸和苯酚在气相中的能量和溶解过程的贡献 , 进而得到 pKa 值 ; Shields 课题组 [77,108] 使用 CPCM[109]连续介质模型计算了20个苯酚和6个羧酸 的pKa. CPCM使用COMSO[110], 是在极化连续介质模 型下, 使用似导体屏蔽模型计算分子极化电荷. 在此 基础上, Klamt等人[111,112]开发了真实溶液下的似导体屏蔽模型(COMSO-RS), RS代表真实溶剂. COMSO-RS对COMSO方法进行了统计热力学上的后处理, 将 量子化学方法的应用拓展到了整个流体热动力学领 域, 包括混合溶剂和不同温度. Ho和Coote [113]比较了使用3种不同的溶剂模型 , CPCM-(UAKS/UAHF), COSMO-RS和SM6对13种中性、10种阳离子型和5种阴离子型碳氢酸的pKa的计算结果, 发现CPCM模型 对中性化合物pKa的计算效果最好, 且该方法对化合物结构不敏感, 可以广泛用于中性碳氢酸的pKa计算, 但对离子型化合物计算结果不理想, 而COMSO溶剂模型对全部化合物的整体计算效果最好.

  理论上来讲, 使用第一性原理的方法不会受到化合物分子结构的限制, 应用场景范围最广泛, 而且可以用来计算化合物在其他溶液中的解离常数. 然而该 方法需要大量的计算资源, 因此不适用于较大分子或者大数据集的化合物. 另外, 与其他方法相比, 其预测精度也没有明显提高, 主要原因是受溶剂模型的准确度的限制, 所以该方法的应用还需要依赖于理论化学的持续不断的发展.

  表2列出了目前主流的pKa预测软件. 已有一些工作比较了这些软件对不同化合物集的应用情 况[6,72,114~117]. 例如, Liao和Nicklaus[6]使用ACD/pKa等 9种软件计算了197个药物分子共261个pKa值(141个酸性基团和120个碱性基团), 发现ACD/pKa软件的精确度最高, 并且能获得所有解离官能团的pKa值. 其他的比较研究也发现ACD/pKa, ADME Boxes和 Marvin往往是预测性能最好的几种软件. 此外, 所有 基于统计学方法(LEER, QSPR)的软件, 计算速度远 远快于基于量子化学方法的应用程序, 并能用于大 数据集, 但SPARC软件由于包含了扰动分子轨道理 论, 速度相对较慢. 例如, 对197个化合物, ADMET Predictor计算时间不超过1 s, 而Jaguar对其中的1个叔胺分子(克冠二胺)就要花费2 d[6].

  然而现有的比较研究多侧重于预测性能, 事实上作为商业软件, 其文档编制、可用性、可靠性、自动化程度以及模型的改进等方面也是很重要的. 由于pKa是化合物最基本的理化性质参数, 对于开发新的、更好的pKa预测软件仍有巨大需求[111].

  在过去的几十年间, 发展了许多预测小分子化 合物pKa的方法, 总体上可大致分为LFER模型、基于不同描述符和统计分析方法的QSPR模型以及使用第一性原理计算的方法. 第一性原理的方法普遍适用于所有化合物, 但其预测精度会受到溶剂模型的限制, 而且需要耗费大量计算资源. LFER方法是发展历史最长、至今仍然在使用的非常成熟的方法, 最成功的应用即ACD/pKa. 近些年来, 量子化学和化学信息学 技术的发展促进了不一样的描述符的开发, 结合 统计分析方法的运用, 产生了很多性能优良的QSPR 模型, 也广泛应用于pKa预测软件的开发. 通常, LFER和QSPR方法运算速度较快, 但用于pKa预测会受到其应用域的限制.由于量子化学描述符具有清晰明确的物理化学 意义, 因此高精度的量化计算(ab initio/DFT)结合恰 当的统计分析方法来发展QSPR模型将是未来pKa计算的一个重要方向. 此外, 实验数据是发展预测方法 的基础, 然而最好的pKa数据集往往来自于发展商业软件的公司, 而方法创新最初总是源自于学术界, 因 此加强工业界与学术界的合作至关重要.

  平台声明：该文观点仅代表作者本人，搜狐号系信息发布平台，搜狐仅提供信息存储空间服务。