赵东霞， 马 悦， 赵 健

(辽宁师范大学 化学化工学院, 辽宁 大连 116029)

应用ABEEMσπ方法估算含有羟基的分子的pKa值

赵东霞， 马 悦， 赵 健

(辽宁师范大学 化学化工学院, 辽宁 大连 116029)

使用B3LYP/6-311++G(d,p)方法，优化获得含有羟基的训练集分子的稳定几何结构，以HF/STO-3G方法所计算的Mulliken电荷为基准，采用线性回归和最小二乘法调试出ABEEMσπ方法计算电荷所需要的参数(价态电负性χ*和价态硬度η*).探讨拟合训练集分子中与羟基相连的C原子和与羟基的H原子的电荷差值与实验pKa值的线性方程.通过该线性方程和ABEEMσπ所计算的电荷，估算出一些含有羟基测试集分子的pKa值.这些分子包括了12个含有羟基的有机小分子；1个Tyr二肽、6个Ser二肽；质子化和中性的Trp-cage蛋白质.使用ABEEMσπ方法所估算的pKa值与实验值很接近.因此，ABEEMσπ方法能快速估算其他含有羟基分子的pKa值.

原子-键电负性均衡原理;电荷;pKa;二肽;Trp-cage蛋白质

1 计算方法

1.1 原子-键电负性均衡模型(ABEEMσπ)

(1)

根据电负性均衡原理，在形成分子时，由于各区域间的电子转移和重新分布，调整后的各区域的有效电负性相等，并且等于分子电负性χmol，即

χa=χb=…=χn=χmol.

(2)

各个区域电荷的和为分子的电荷，即

qa+qb+…+qn=qmol.

(3)

将式(1)和式(2)以及电荷限制条件式(3)联立，就能获得分子各个区域的电荷.

1.2 计算细节

选取的训练集分子为a01～a10分子，测试集分子包括b01～b11，其中b10有2种不同构象含有羟基的12个有机小分子、1个Tyr二肽和6个Ser二肽以及质子化和中性的Trp-cage1蛋白质.使用Gaussian09[10]程序，用B3LYP/6-311++G(d,p)方法，优化获得训练集和测试集分子的稳定几何结构.使用HF/STO-3G计算训练集分子的Mulliken电荷，拟合确定ABEEMσπ参数(包括价态电负性χ*和价态硬度η*).使用ABEEMσπ参数获得训练集分子的电荷分布.拟合得到训练集分子中与官能团相关的原子电荷与pKa值的线性方程，进而估算测试集分子的pKa值.

2 结果与讨论

2.1 ABEEMσπ方法的参数

以HF/STO-3G方法所计算的训练集分子的Mulliken电荷为基准，调试出ABEEMσπ方法中所涉及新增位点的价态电负性χ*和价态硬度2η*参数，列于表1中.

表1 训练集分子中新定义ABEEMσπ参数(价态电负性χ*和价态硬度2η*)

将表1中新定义的ABEEMσπ参数的标号标记在原子的右侧，示于图1中.在ABEEMσπ程序中，

图1 ABEEMσπ参数标号示意图

各个位点的标号用6位数来表示，以106 257为例，标号的第1位是1代表原子区域；第2和第3位分别代表原子序数，例如06代表碳原子；第4位是2，代表以双键的形式与其他原子相连，即碳原子的周围有1个双键；第5位和第6位一起表示新添加的标号，即57代表新添加的标号.

2.2 ABEEMσπ方法所计算的电荷

依据ABEEMσπ方法，使用已有的[5-7]和本文所拟合的ABEEMσπ参数，计算出训练集分子和测试集分子的电荷分布.表2列出了ABEEMσπ所计算的电荷与从头算电荷的线性方程的斜率k，截距b，线性相关系数R.

表2 训练集和测试集分子的ABEEMσπ与从头算电荷的线性相关方程以及相关系数

从表2可以看出，b01～b11、1个Tyr和6个Ser的ABEEMσπ电荷与从头算电荷线性方程的斜率在1.000 0左右，截距趋近于0，线性相关系数在0.96以上；质子化和中性Trp-cage的ABEEMσπ电荷与从头算电荷线性方程的斜率在1.100 0左右，截距趋近于0，线性相关系数在0.94以上.可以看出ABEEMσπ所计算的电荷与从头算方法有很好的线性关系，说明本文所拟合的ABEEMσπ参数是合理和可转移的.

2.3 训练集分子中原子的电荷与pKa线性关系

表3 ABEEMσπ模型下，训练集分子的C原子和H原子的电荷

2.4 ABEEMσπ方法所估算的测试集分子pKa值

使用上述的方程和ABEEMσπ所计算的电荷，检验了测试集分子的pKa值.测试集分子中与羟基相连的碳原子和羟基中氢原子的ABEEMσπ电荷以及pKa的预测值和实验值均列于表4中.

表4 ABEEMσπ模型下，测试集分子的C和H原子的电荷以及pKa值

从表4中可以看出，对于b01～b11中的含羟基的有机小分子，大部分分子的pKa估算值与实验值的误差都在1.1以内，除了b01和b1001.b01的误差为1.78，b1001的误差为2.88，其原因是分子的结构不同，所计算的电荷不同，从而所预测的pKa值不同.图2中的b1001分子和b1002分子是同一种分子的2个不同构象，其中，b1001分子中羟基上的氢原子偏向相邻甲氧基中的氧原子，形成分子内氢键，因此，分子中羟基上的H较难解离，所以pKa值较大.而b1002分子中羟基上的氢原子远离甲氧基中的氧原子及其孤对电子，不形成分子内氢键，所以相比b1001分子而言，b1002分子羟基中的H更容易解离.使用ABEEMσπ电荷得到的pKa的预测值也可推算出b1002分子比b1001分子更容易解离，与分析结果一致.因此pKa值的大小还与羟基上H原子与其他原子的相互作用有关，与氢原子所处的周围环境有关.

Tyr二肽的pKa估算值与实验值的误差是0.26，二者非常接近.选取的Ser二肽的6个构象中，Ser1二肽和Ser4二肽的pKa预测值远大于其他构象的pKa.虽然没有找到6个构象对应的pKa实验值，但是从结构分析，可以得出Ser1和Ser4的羟基的H原子与分子肽键的O原子有氢键作用，不易解离出来，pKa值较大.Ser4羟基的H原子与分子肽键的O原子的距离为0.197 91nm，其距离最近，氢键作用最大，它们间最难解离，pKa值应该最大，这与估算结果相吻合.

图2 部分测试集分子的结构

对于质子化的Trp-cage蛋白质(含有304个原子)和中性的Trp-cage1蛋白质，主要估算了其中3个含有羟基的氨基酸的pKa值，分别是暴露于溶剂部分的Ser13，被埋藏于蛋白质内部的Ser14，及暴露于溶剂部分的Tyr3，结构示于图2中.Matasui等人[19]通过热力学循环和动力学模拟等方法预测出Trp-cage蛋白质中Ser13、Ser14、Tyr3残基的pKa值分别是14.69、20.80和10.10，从而判断出这3种氨基酸解离能力的大小顺序是Tyr3>Ser13>Ser14.从表4可以看出，Trp-cage和Trp-cage1蛋白质中Tyr3、Ser13、Ser14的pKa预测值和Matasui等人[19]计算的pKa值很接近.对于质子化的Trp-cage蛋白质残基，二者的pKa预测值误差在1.00以内；对于中性的Trp-cage1蛋白质残基，二者的pKa预测值误差在0.59以内.并且2种情况下氨基酸的pKa预测值大小顺序均为Tyr3Ser13>Ser14，这与Matasui等人[19]预测结果一致.除此之外，计算Trp-cage蛋白质的电荷，ABEEMσπ方法的计算速度远快于从头算的速度，从头算需要18 min左右，而ABEEMσπ仅需要1 s.

3 结 论

应用ABEEMσπ方法对含羟基分子的pKa值进行了预测，得到的小分子和二肽分子的pKa预测值与实验值十分接近，蛋白质中氨基酸残基的pKa预测值与使用溶剂化能得到的结果相接近.因此，应用ABEEMσπ方法可以快速、准确的计算出分子中原子的电荷，进而可以很好地估算含有羟基分子的pKa值.

[1] TISHMACK P A,BASHFORD D,HARMS E,et al.Use of 1H NMR spectroscopy and computer simulations to analyze histidine pKachanges in a protein tyrosine phosphatase:experimental and theoretical determination of electrostatic properties in a small protein[J].Biochemistry,1997,36(39):11984-11994.

[2] ALI S T,WKARAMAT S,WKA J,et al.Theoretical prediction of pKavalues of seleninic,selenenic,sulfinic and carboxylic acids by quantum-chemical methods[J].J Phys Chem A,2010,114(47):12470-12478.

[3] SCHÜÜRMANN G,COSSI M,BARONE V,et al.Prediction of the pKaof carboxylic acids using the ab initio continuum-solvation model PCM-UAHF[J].J Phys Chem A,1998,17(33):6706-6712.

[4] UGUR I,MARION A,PARANT S,et al.Rationalization of the pKavalues of alcohols and thiols using atomic charge descriptors and its application to the prediction of amino acid pKa′s[J].J Chem Inf Model,2014,54(8):2200-2213.

[5] WANG C S,YANG Z Z.Atom-bond electronegativity equalization method.II.Lone-pair electron model[J].J Chem Phys,1999,110(13):6189-6197.

[6] YAO C,YANG Z Z.General atom-bond electronegativity equalization method and its application in prediction of charge distributions in polypeptide[J].Chem Phys Lett,2000,316(3):324-329.

[7] ZHAO D X,LIU C,WANG F F,et al.Development of a polarizable force field using multiple fluctuating charges per atom[J].J Chem Theory Comput,2010,6(3):795-804.

[8] 杨忠志,刘娇,宫利东,等.并行程序实现ABEEMσπ模型电荷分布计算[J].辽宁师范大学学报(自然科学版),2008,31(2):177-180.

[9] 宫利东,徐燕.应用ABEEMσπ模型研究含氮杂环化合物的电荷分布[J].辽宁师范大学学报(自然科学版),2008,31(1):64-66.

[10] FRISCH M J,TRUCKS G W,SCHLEGEL H B,et al.Gaussian 09[CP].Wallingford:Gaussian Inc,2009.

[11] ALBERT A,PHILLIPS J N.Ionization constants of heterocyclic substances.Part Ⅱ.Hydroxy-derivatives of nitrogenous six-membered ring-compounds[J].J Chem Soc,1956：1294-1304.

[12] BOLTON P D,HALL F M,REECE I H.Effects of substituents on the thermodynamic functions of ionisation of meta-substituted phenols[J].J Chem Soc B,1967：709-712.

[13] HAMANN S D,LINTON M.Influence of pressure on the ionization of substituted phenols[J].J Chem Soc Faraday T,1974,70(1):2239-2249.

[14] LIDE D R.CRC handbook of chemistry and physics[M].Boca Raton Florida:CRC Press,2006:87.

[15] KO H C,O′HARA W F,HU T,et al.Ionization of substituted phenols in aqueous solution[J].J Am Chem Soc,1964,86(5):1003-1004.

[16] TAKAHASHI S,COHEN L A,MILLER H K,et al.Calculation of the pKavalues of alcohols from .sigma.constants and from the carbonyl frequencies of their esters[J].J Org Chem,1971,36(9):1205-1209.

[17] VÖLGYI G,RUIZ R,BOX K,et al.Potentiometric and spectrophotometric pKadetermination of water-insoluble compounds:validation study in a new cosolvent system[J].Anal Chim Acta,2007,583(2):418-428.

[18] THURLKILL R L,GRIMSLEY G R,SCHOLTZ J M,et al.pKvalues of the ionizable groups of proteins[J].Protein Sci,2006,15(5):1214-1218.

[19] MATASUI T,BABA T,WKAMIYA K,et al.An accurate density functional theory based estimation of pKavalues of polar residues combined with experimental data:from amino acids to minimal proteins[J].Phys Chem Chem Phys,2012,14(12):4181-4187.

Estimation of pKavalues of hydroxyl molecules by using ABEEMσπmethod

ZHAODongxia,MaYue,ZHAOJian

(School of Chemistry and Chemical Engineering, Liaoning Normal University, Dalian 116029, China)

The B3LYP/6-311++G(d,p) method is used to optimize a training set of small molecules with an OH group and to obtain their geometry structures.ABEEMσπparameters including the valence-state electronegativityχ*and the valence-state hardnessη*have been determined through a linear regression and least-square optimization procedure to reproduce HF/STO-3G Mulliken charge distribution for the molecules mentioned above.We have obtained the linear equation of the difference between C connected with OH group atomic charges and H atom with experimental pKa.ThepKavaluesoftestsetmoleculeswithOHgroupareestimatedbythislinearequation.ThesemoleculesincludetwelveorganicsmallmoleculeswithanOHgroup,oneTyrdipeptide,sixSerdipeptides,protonatedandneutralTrp-cageproteins.ThepKavaluesestimatedbyABEEMσπmethod are close to the experimental pKavalues.Thus,ABEEMσπmethod can be used to fast estimate pKavaluesoftheothermoleculescontainingOHgroup.

ABEEMσπ;charge;pKa;dipeptide;Trp-cage protein

2016-05-08 基金项目：国家自然科学基金资助项目(21133005;21473083) 作者简介：赵东霞(1971-),女,河南长垣人,辽宁师范大学教授,博士,博士生导师.

1000-1735(2017)01-0052-06

10.11679/lsxblk2017010052

O641.121

A