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

趙東霞， 馬 悅， 趙 健

(遼寧師范大學 化學化工學院, 遼寧 大連 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