有機(jī)磷類化合物大鼠急性毒性QSAR模型構(gòu)建與毒性機(jī)制研究

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

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

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

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

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

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

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

1 材料與方法(Materials and methods)

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

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

1.2 分子描述符的計(jì)算

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

1.3 模型的構(gòu)建

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

圖1 定量結(jié)構(gòu)活性關(guān)系(QSAR)模型構(gòu)建及有機(jī)磷類化合物(OPs)的急性毒性半致死量(LD50)預(yù)測(cè)的工作流程圖注：GA-MLR表示基于遺傳算法的多元線性回歸模型；SW-MLR表示基于逐步算法的多元線性回歸模型；R2表示決定系數(shù)；MAE表示平均絕對(duì)誤差。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 模型評(píng)估

1.5 模型應(yīng)用域

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

(1)

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

(2)

(3)

1.6 分子對(duì)接

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

2 結(jié)果(Results)

2.1 預(yù)測(cè)模型

逐步算法的參數(shù)設(shè)置為進(jìn)入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

對(duì)于每個(gè)模型還執(zhí)行了2 000次Y-隨機(jī)測(cè)試。隨機(jī)測(cè)試結(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實(shí)驗(yàn)值與預(yù)測(cè)值擬合關(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與毒性終點(diǎn)logLD50正相關(guān)。ATS3s、MATS8i分別是基于原子的內(nèi)在狀態(tài)和電離電位來表征化合物拓?fù)浣Y(jié)構(gòu)的分子描述符[34]。maxssO則代表分子結(jié)構(gòu)中—O—基團(tuán)的最大電性拓?fù)錉顟B(tài)，前人的研究中也選用maxssO來描述一些無毒的抗瘧化合物[35]，這與我們研究中maxssO描述符在SW-MLR方程中的正系數(shù)一致。分子描述符GATS3e、JGI7、MATS3i和topoDiameter則與毒性終點(diǎn)logLD50負(fù)相關(guān)。GATS3e和JGI7描述符與化合物中高電負(fù)性的原子相關(guān)[36-37]。MATS3i和topoDiameter是特定拓?fù)浣Y(jié)構(gòu)的分子描述符。

在GA-MLR方程中，分子描述符AATSC7m、nsssCH和VR3_Dzp與毒性終點(diǎn)logLD50正相關(guān)。AATSC7m和nsssCH分別對(duì)應(yīng)OPs分子結(jié)構(gòu)中R基上支鏈的長(zhǎng)度[38]、支鏈或取代基烷基化的程度[39]。VR3_Dzp則是比較抽象的二維描述符，代表極化率加權(quán)的Barysz矩陣的特征向量[40]。分子描述符ATSC3i、AVP_4、JGI7和GATS4p則與毒性終點(diǎn)logLD50負(fù)相關(guān)。這4個(gè)描述符與分子的疏水性、極性和原子的電負(fù)性相關(guān)[36,41]。對(duì)于2個(gè)模型分子描述符的解釋從一定程度上揭示了OPs的化學(xué)結(jié)構(gòu)與毒性相關(guān)的內(nèi)在因素，對(duì)于識(shí)別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)為離群點(diǎn)。結(jié)論表明本研究的2種模型能夠較好地預(yù)測(cè)OPs類化合物對(duì)大鼠的急性毒性數(shù)據(jù)。表2中給出了使用2種QSAR模型預(yù)測(cè)的沒有實(shí)驗(yàn)數(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 分子對(duì)接

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

AutoDock Vina對(duì)接的結(jié)果顯示其中8個(gè)OPs的對(duì)接結(jié)合能低于閾值[43]-6 kcal·mol-1，表明了潛在的結(jié)合能力以及進(jìn)一步可能造成的神經(jīng)毒性。此外，取代基或支鏈烷基化程度較高的化合物(如磷酸三(3-氯丙基)酯)對(duì)接結(jié)合能也較高，這表明更低的結(jié)合能力與更低的神經(jīng)毒性，這一結(jié)果與GA-MLR模型中分子描述符nsssCH的含義一致。因此我們認(rèn)為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分子描述符，對(duì)OPs化合物的大鼠口服急性毒性建立了多元線性回歸QSAR模型，SW-MLR模型具有更好的擬合能力，GA-MLR模型有著更好的預(yù)測(cè)能力。對(duì)模型分子描述符的分析顯示了OPs的R基的高烷基化和支鏈的延長(zhǎng)能夠降低毒性。此外，通過模型預(yù)測(cè)了數(shù)據(jù)集中沒有實(shí)驗(yàn)數(shù)據(jù)的9種OPs的毒性，并以分子對(duì)接方法檢測(cè)其與BChE的潛在結(jié)合能力。所構(gòu)建模型能夠方便有效地預(yù)測(cè)OPs類化合物的毒性，并為新化學(xué)品的合成與環(huán)境保護(hù)政策的制定提供幫助。

表2 SW-MLR和GA-MLR對(duì)9種OPs的LD50預(yù)測(cè)值、沒有實(shí)驗(yàn)值的9種OPs與丁酰膽堿酯酶(BChE)對(duì)接的能量與關(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)分子對(duì)接結(jié)果的可視化注：A～I(xiàn)對(duì)應(yīng)表2中化合物序號(hào)；疏水性相互作用為灰色虛線，平行π-堆疊為亮綠色虛線，垂直π-堆疊為深綠色虛線，氫鍵為藍(lán)色實(shí)線，鹽橋?yàn)辄S色虛線。Fig. 4 Ligand interaction diagram of the main interactions of nine OPs in the active site of BChE (5K5E)Note: A～I(xiàn) 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.

通訊作者簡(jiǎn)介：查金苗(1975—)，男，博士，研究員，主要研究方向?yàn)樗Ｐ蜕矬w系的構(gòu)建與發(fā)展、水環(huán)境生物毒性測(cè)試方法、環(huán)境內(nèi)分泌干擾物的篩選技術(shù)研究、環(huán)境污染物對(duì)水生生物分子毒理機(jī)制和水生態(tài)系統(tǒng)完整性評(píng)估方法等。

共同通訊作者簡(jiǎn)介：閆賽紅(1988—)，女，博士，助理研究員，主要研究方向?yàn)樗鷳B(tài)毒理學(xué)。