鄭子廷，閆賽紅,#，查金苗,*

1. 中國科學(xué)院生態(tài)環(huán)境研究中心，環(huán)境水質(zhì)學(xué)國家重點實驗室，北京 100085 2. 中國科學(xué)院大學(xué)，北京 100049

有機磷類化合物(organophosphorus compounds, OPs)包括有機磷阻燃劑(organophosphorus flame retardants, OPFRs)、有機磷殺蟲劑(organophosphorus pesticides, OPPs)和其他含磷的有機化合物。工業(yè)生產(chǎn)上，OFPRs作為溴化阻燃劑的低毒替代品被應(yīng)用于家電、建筑和電子領(lǐng)域；農(nóng)業(yè)上，OPPs則被廣泛用于殺蟲和害蟲控制[1-2]。此外在醫(yī)學(xué)領(lǐng)域，一些抗病毒的藥物中也含有OPs[3]。因此，OPs容易通過磨損、揮發(fā)和溶解等方式進入各種環(huán)境介質(zhì)中。目前科研人員已在諸如河流、地下水、海水、空氣和室內(nèi)等環(huán)境介質(zhì)中檢測到了OPs[4-6]。甚至在東亞地區(qū)的人類母乳樣本中檢測到10余種OPs，最高濃度為70 ng·g-1[7]。越來越廣泛存在的OPs可能造成的危害和健康風險引起了研究者們的關(guān)注。

OPs具有神經(jīng)毒性[8-9]、內(nèi)分泌毒性[10]和致癌性[11]。研究表明，接觸OPPs(如毒死蜱)會引起人類的神經(jīng)系統(tǒng)發(fā)育缺陷[12]。OPFRs如磷酸三甲苯(TMPP)、磷酸三苯酯(TPHP)和磷酸三(2,3-二溴丙基)酯(TDBPP)則顯示出對核受體(nuclear receptors, NRs)的親和力和內(nèi)分泌的干擾作用[13-14]。上述研究及更多OPs的毒性毒理研究往往都是通過傳統(tǒng)的體內(nèi)(invivo)或者體外(invitro)實驗進行的，這些方法需要耗費大量的時間與金錢。此外，實驗手段無法完整覆蓋一整類化合物的毒性數(shù)據(jù)，并且測試在實驗室間存在誤差[15]。因此，迫切需要一種有效的替代方法，對化合物毒性數(shù)據(jù)進行評估。

隨著計算科學(xué)的發(fā)展，計算毒理學(xué)已經(jīng)成為一種具有成本效益且快速的實驗替代方法[16]。計算毒理學(xué)整合各種來源的數(shù)據(jù)，開發(fā)基于數(shù)學(xué)、統(tǒng)計科學(xué)和機器學(xué)習(xí)的模型，對化學(xué)品進行高效和高通量的定性或定量預(yù)測[17]。定量結(jié)構(gòu)活性關(guān)系(quantitative structure-activity relationship, QSAR)模型便是計算毒理學(xué)領(lǐng)域被廣泛利用的模型之一，其原理是建立化合物的化學(xué)結(jié)構(gòu)與各種毒理學(xué)終點效應(yīng)和生物活性之間的函數(shù)關(guān)系[18]。QSAR模型被研究者們廣泛應(yīng)用于藥物發(fā)現(xiàn)、毒理學(xué)和生物學(xué)等領(lǐng)域[19]。除了建立化學(xué)結(jié)構(gòu)與毒性數(shù)據(jù)之間的關(guān)系模型外，QSAR模型還能通過分子描述符的解釋從一定程度進行毒理機制的闡述。Alves等[20]使用QSAR方法為COVID-19病毒建立了抑制性藥物預(yù)測模型，并應(yīng)用該模型從小分子數(shù)據(jù)庫中高通量地篩選了11個候選藥物，其中3個被美國轉(zhuǎn)化科學(xué)推動中心(National Center for Advancing Translational Sciences, NCATS)實驗證實為有效的病毒抑制劑。Jeon等[21]則收集大量實驗數(shù)據(jù)，為魚胚胎急性毒性建立了一個通用的基線毒性QSAR模型，并使用該模型篩選毒性化合物。大量的研究表明QSAR模型用于化合物的毒性預(yù)測行之有效。

分子對接是研究化合物與目標蛋白結(jié)合部位行為的工具，能夠通過計算手段探究化合物與特定靶點蛋白的結(jié)合特異性，從而從一定程度上反映化合物的毒性。Zou等[22]使用分子對接得到的配體-受體相互作用能量(binding energy)來反映化學(xué)品的毒性作用。Li等[23]通過對接羥基多溴聯(lián)苯醚(HO-PBDEs)與人類甲狀腺激素受體，從一定程度上揭示了HO-PBDEs的毒性機制。

現(xiàn)有的一些針對大鼠急性口服毒性的QSAR模型存在沒有使用測試集[24]、選用的分子描述符數(shù)量不規(guī)范[25]等問題。在本研究中，完善地建立了基于逐步算法的線性回歸(SW-MLR)和基于遺傳算法的線性回歸(GA-MLR)QSAR模型用于預(yù)測大鼠的急性口服毒性。研究收集的數(shù)據(jù)集中含有9種沒有毒性實驗數(shù)據(jù)的OPs，最終開發(fā)完成的QSAR模型將用于預(yù)測這9種OPs的急性毒性。分子對接也將作為一個輔助工具來評估這9種OPs潛在的神經(jīng)毒性。

1 材料與方法(Materials and methods)

1.1 數(shù)據(jù)集的收集和劃分

本研究首先收集了OPs對大鼠急性口服毒性的數(shù)據(jù)。檢索并從文獻中提取了62種環(huán)境中常見的OPs清單[26-27]。通過查詢EPA CompTox數(shù)據(jù)庫(https://comptox.epa.gov/dashboard)，收集了同一實驗條件下獲得的53種OPs急性毒性數(shù)據(jù)LD50(導(dǎo)致一組實驗動物中剛好死亡一半的化合物濃度)。獲取的OPs的急性毒性數(shù)值使用log函數(shù)進行規(guī)范化處理。整個數(shù)據(jù)集按照4∶1的比例隨機劃分為訓(xùn)練集(n=43)和測試集(n=10)，沒有毒性實驗數(shù)據(jù)的9種OPs作為最后的預(yù)測集(n=9)。在之后的計算實驗中，訓(xùn)練集僅用于建立模型，測試集則用于評估模型的預(yù)測能力，這種劃分方式有利于驗證QSAR模型的擬合和泛化能力。

1.2 分子描述符的計算

數(shù)據(jù)集中所有OPs的結(jié)構(gòu)文件(.sdf)來自PubChem(https://pubchem.ncbi.nlm.nih.gov/)數(shù)據(jù)庫。并使用PaDEL開源軟件[28]計算所有化合物的分子結(jié)構(gòu)描述符，共得到1 444個分子結(jié)構(gòu)描述符，包括(1)一維結(jié)構(gòu)描述符：如原子和化學(xué)鍵的數(shù)量、分子量等；(2)二維拓撲描述符：如拓撲指數(shù)、拓撲電荷和自相關(guān)系數(shù)等；(3)電荷描述符、原子電狀態(tài)描述符等。本研究中使用了PaDEL計算的一維和二維描述符，因為這些類型的描述符計算速度快，可重復(fù)性好，并且不需要分子結(jié)構(gòu)能量最小化過程，該過程中不同實驗室間可能會在計算上產(chǎn)生誤差。計算結(jié)束后，我們對得到的分子描述符數(shù)據(jù)進行了預(yù)處理，剔除了方差為0的單列分子描述符數(shù)據(jù)，相關(guān)系數(shù)>0.99的多組描述符最終僅保留一組。經(jīng)過預(yù)處理后，最終保留了833組分子描述符用于后續(xù)分析。

1.3 模型的構(gòu)建

分別使用逐步算法(SW)和遺傳算法(GA)進行分子描述符的選擇。2種算法選擇的分子描述符變量通過基于普通最小二乘法(OLS)的多元線性回歸(MLR)方法關(guān)聯(lián)毒性終點并建立QSAR模型，之后進行模型的內(nèi)部和外部評估。研究中使用SPSS 21.0軟件建立了SW-MLR模型，使用DTC-QSAR開源軟件[29]建立了GA-MLR模型。所有參數(shù)的統(tǒng)計分析過程在EXCEL軟件中完成。本研究詳細的QSAR模型構(gòu)建及OPs的LD50預(yù)測工作流程如圖1所示。

圖1 定量結(jié)構(gòu)活性關(guān)系(QSAR)模型構(gòu)建及有機磷類化合物(OPs)的急性毒性半致死量(LD50)預(yù)測的工作流程圖注：GA-MLR表示基于遺傳算法的多元線性回歸模型；SW-MLR表示基于逐步算法的多元線性回歸模型；R2表示決定系數(shù)；MAE表示平均絕對誤差。Fig. 1 Workflow diagram for quantitative structure-activity relationship (QSAR) models construction and acute toxicity semi-lethal dose (LD50) prediction of organophosphorus compounds (OPs)Note: GA-MLR represents the multiple linear regression model based on genetic algorithm; SW-MLR represents the multiple linear regression model based on stepwise algorithm; R2 represents the coefficient of determination; MAE represents mean absolute error.

1.4 模型評估

1.5 模型應(yīng)用域

使用Williams圖(又稱杠桿方法)[30]對建立的QSAR模型進行應(yīng)用域分析?；诿弊泳仃?一類投影矩陣)得出的值(h)與標準化殘差(δ)的比較。h大于警告杠桿值(h*)的化合物被認為是結(jié)構(gòu)異常值，標準化殘差>2.5的化合物則被認為是計算異常值。Williams圖的公式定義如下:

(1)

表1 本研究中用于模型評價的公式及參數(shù)Table 1 Equations and parameters used for model evaluation in this study

(2)

(3)

1.6 分子對接

我們先前的實驗研究證實了稀有鮈鯽的丁酰膽堿酯酶(BChE)能夠被一些OPs顯著抑制，從而引起神經(jīng)毒性[31]。這項研究中，借助分子對接工具探究預(yù)測集中9種OPs對哺乳動物類的BChE的結(jié)合與抑制能力。從蛋白質(zhì)數(shù)據(jù)庫中檢索到BChE的人類蛋白質(zhì)結(jié)構(gòu)(PDB ID: 5K5E)作為對接蛋白。通過AutoDockTools4程序進行預(yù)處理，去除水分、添加氫并分配電荷。對接網(wǎng)格的大小設(shè)置為52 ?×52 ?×52 ?，并使用AutoDock Vina程序[32]進行分子對接。對接前，從BChE蛋白中提取了共結(jié)晶化合物并重新對接至原蛋白進行對接程序的驗證。最終使用PLIP工具[33]分析了對接后化合物與蛋白質(zhì)的結(jié)合模式。

2 結(jié)果(Results)

2.1 預(yù)測模型

逐步算法的參數(shù)設(shè)置為進入F值：<0.010，逐出F值：>0.050。建立的SW-MLR模型如下：

logLD50=-0.154(±1.5150)+0.008(±0.0018)ATS3s-2.783(±0.5247)GATS3e+2.054(±0.4083)MATS8i-88.853(±13.1459)JGI7-2.871(±0.6230)MATS3i+1.257(±0.2907)maxssO-(±0.0553)0.180topoDiameter

遺傳算法的參數(shù)設(shè)置為最大迭代次數(shù)：100，方程參數(shù)數(shù)量：7，初始方程數(shù)：100，每代方程數(shù)：30，適應(yīng)函數(shù)：MAE。建立的GA-MLR模型如下；

logLD50=5.608(±0.7354)+0.0067(±0.0049)AATSC7m-0.012(±0.0088)ATSC3i+0.0133(±0.0087)VR3_Dzp-3.1038(±0.9454)AVP_4-85.9765(±14.6601)JGI7+0.263(±0.1021)nsssCH-1.2293(±0.3834)GATS4p

對于每個模型還執(zhí)行了2 000次Y-隨機測試。隨機測試結(jié)果顯示SW-MLR和GA-MLR的最大偽R2值分別為0.488和0.504，最大偽Q2值分別為0.299和0.234。所有的偽R2和偽Q2都很低，這表明原始模型的結(jié)果并不是偶然關(guān)聯(lián)的。SW-MLR和GA-MLR模型的擬合關(guān)系結(jié)果如圖2所示。

圖2 SW-MLR(a)和GA-MLR(b)模型中OPs實驗值與預(yù)測值擬合關(guān)系Fig. 2 Fitting relationship between experimental and predicted values of OPs in SW-MLR (a) and GA-MLR (b) models

2.2 分子描述符解釋

在SW-MLR方程中，分子描述符ATS3s、MATS8i和maxssO與毒性終點logLD50正相關(guān)。ATS3s、MATS8i分別是基于原子的內(nèi)在狀態(tài)和電離電位來表征化合物拓撲結(jié)構(gòu)的分子描述符[34]。maxssO則代表分子結(jié)構(gòu)中—O—基團的最大電性拓撲狀態(tài)，前人的研究中也選用maxssO來描述一些無毒的抗瘧化合物[35]，這與我們研究中maxssO描述符在SW-MLR方程中的正系數(shù)一致。分子描述符GATS3e、JGI7、MATS3i和topoDiameter則與毒性終點logLD50負相關(guān)。GATS3e和JGI7描述符與化合物中高電負性的原子相關(guān)[36-37]。MATS3i和topoDiameter是特定拓撲結(jié)構(gòu)的分子描述符。

在GA-MLR方程中，分子描述符AATSC7m、nsssCH和VR3_Dzp與毒性終點logLD50正相關(guān)。AATSC7m和nsssCH分別對應(yīng)OPs分子結(jié)構(gòu)中R基上支鏈的長度[38]、支鏈或取代基烷基化的程度[39]。VR3_Dzp則是比較抽象的二維描述符，代表極化率加權(quán)的Barysz矩陣的特征向量[40]。分子描述符ATSC3i、AVP_4、JGI7和GATS4p則與毒性終點logLD50負相關(guān)。這4個描述符與分子的疏水性、極性和原子的電負性相關(guān)[36,41]。對于2個模型分子描述符的解釋從一定程度上揭示了OPs的化學(xué)結(jié)構(gòu)與毒性相關(guān)的內(nèi)在因素，對于識別OPs的毒性起到關(guān)鍵作用。

2.3 模型應(yīng)用域

建立的SW-MLR模型和GA-MLR模型應(yīng)用域表征Willimas圖如圖3所示。SW-MLR模型中，所有OPs類化合物都很好地落入了應(yīng)用域范圍內(nèi)。GA-MLR模型中，僅有一種化合物(dimethoate, CAS: 60-51-5,h=0.570>h*=0.558)為離群點。結(jié)論表明本研究的2種模型能夠較好地預(yù)測OPs類化合物對大鼠的急性毒性數(shù)據(jù)。表2中給出了使用2種QSAR模型預(yù)測的沒有實驗數(shù)據(jù)的9種OPs的急性毒性結(jié)果。

圖3 SW-MLR(a)和GA-MLR(b)的Williams應(yīng)用域圖Fig. 3 Williams plotting for SW-MLR (a) and GA-MLR (b)

2.4 模型比較

2.5 分子對接

AutoDock Vina重新對接BChE(5K5E)的原配體的結(jié)合能結(jié)果為-11.9 kcal·mol-1，均方根偏差(RMSD)為0.89 ?，說明對接程序適用于進一步的對接與研究。使用AutoDock Vina軟件對接預(yù)測集(9種OPs)與BChE(5K5E)的最低對接能量結(jié)果如表2所示。PLIP工具對OPs與BChE的結(jié)合模式可視化如圖4所示。

AutoDock Vina對接的結(jié)果顯示其中8個OPs的對接結(jié)合能低于閾值[43]-6 kcal·mol-1，表明了潛在的結(jié)合能力以及進一步可能造成的神經(jīng)毒性。此外，取代基或支鏈烷基化程度較高的化合物(如磷酸三(3-氯丙基)酯)對接結(jié)合能也較高，這表明更低的結(jié)合能力與更低的神經(jīng)毒性，這一結(jié)果與GA-MLR模型中分子描述符nsssCH的含義一致。因此我們認為OPs的R基及取代基的烷基化能夠降低毒性。由圖4可知，BChE與OPs結(jié)合的關(guān)鍵氨基酸殘基有Gly116和Ser198(氫鍵)；His428、Trp332和Phe329(鹽橋)；Trp231、Phe329、Trp82、Tyr332和His438(π-堆疊)。這些關(guān)鍵氨基酸殘基與前人的研究較為一致[44-45]。

3 討論(Discussion)

研究采用了SW和GA算法，基于PaDEL分子描述符，對OPs化合物的大鼠口服急性毒性建立了多元線性回歸QSAR模型，SW-MLR模型具有更好的擬合能力，GA-MLR模型有著更好的預(yù)測能力。對模型分子描述符的分析顯示了OPs的R基的高烷基化和支鏈的延長能夠降低毒性。此外，通過模型預(yù)測了數(shù)據(jù)集中沒有實驗數(shù)據(jù)的9種OPs的毒性，并以分子對接方法檢測其與BChE的潛在結(jié)合能力。所構(gòu)建模型能夠方便有效地預(yù)測OPs類化合物的毒性，并為新化學(xué)品的合成與環(huán)境保護政策的制定提供幫助。

表2 SW-MLR和GA-MLR對9種OPs的LD50預(yù)測值、沒有實驗值的9種OPs與丁酰膽堿酯酶(BChE)對接的能量與關(guān)鍵殘基Table 2 Predicted LD50 values by SW-MLR and GA-MLR for nine OPs, docking energy and key residue results by docking to butyrylcholinesterase (BChE)

圖4 9種OPs與BChE(5K5E)分子對接結(jié)果的可視化注：A～I對應(yīng)表2中化合物序號；疏水性相互作用為灰色虛線，平行π-堆疊為亮綠色虛線，垂直π-堆疊為深綠色虛線，氫鍵為藍色實線，鹽橋為黃色虛線。Fig. 4 Ligand interaction diagram of the main interactions of nine OPs in the active site of BChE (5K5E)Note: A～I represents the compounds in Table 2; hydrophobic interactions are represented by gray dashes; π-stackings are represented by bright green (parallel) and dark green (perpendicular) dashes; hydrogen bonds are represented by blue solid lines; salt bridges are represented by yellow dashes.

通訊作者簡介：查金苗(1975—)，男，博士，研究員，主要研究方向為水生模型生物體系的構(gòu)建與發(fā)展、水環(huán)境生物毒性測試方法、環(huán)境內(nèi)分泌干擾物的篩選技術(shù)研究、環(huán)境污染物對水生生物分子毒理機制和水生態(tài)系統(tǒng)完整性評估方法等。

共同通訊作者簡介：閆賽紅(1988—)，女，博士，助理研究員，主要研究方向為水生態(tài)毒理學(xué)。