基于QSPR方法的活性化合物熱穩(wěn)定性預(yù)測研究

高 月，厲 鵬，萬可風(fēng)，張子炎，苑 媛，張宏哲

(1.中石化安全工程研究院有限公司化學(xué)品安全控制國家重點(diǎn)實(shí)驗(yàn)室，山東 青島 266000；2.應(yīng)急管理部化學(xué)品登記中心，山東 青島 266000)

0 引言

活性化合物易發(fā)生劇烈反應(yīng)使化學(xué)能轉(zhuǎn)化為熱能或者機(jī)械能，是具有潛在危害的不穩(wěn)定物質(zhì)。在工業(yè)生產(chǎn)及儲運(yùn)過程中，以有機(jī)過氧化物和硝基化合物為代表的活性化合物引發(fā)許多嚴(yán)重的熱失控事故，在全球范圍內(nèi)造成大量人員傷亡[1]。活性化合物的熱穩(wěn)定性是固有屬性，評判熱穩(wěn)定性的常見參數(shù)包括起始放熱溫度(To)和分解熱(-ΔHd)，兩者分別反映熱危害發(fā)生的可能性和嚴(yán)重性[1]。傳統(tǒng)獲取活性化合物熱穩(wěn)定性參數(shù)的方法是實(shí)驗(yàn)方法，如采用差示掃描量熱法(Differential Scanning Calorimetry,DSC)等，然而由于實(shí)驗(yàn)方法具有滯后性、工作量大、費(fèi)用昂貴且危險(xiǎn)的特點(diǎn)，有必要探究1種簡便、快速的方法來評估活性化合物的熱穩(wěn)定性。

定量結(jié)構(gòu)-性質(zhì)相關(guān)性(Quantitative Structure-Property Relationship,QSPR)是通過將分子結(jié)構(gòu)參數(shù)與目標(biāo)性質(zhì)之間的內(nèi)在定量關(guān)系關(guān)聯(lián)，建立構(gòu)效關(guān)系模型，基于分子結(jié)構(gòu)信息預(yù)測化合物性質(zhì)的1種方法[2-3]。QSPR方法可以彌補(bǔ)實(shí)驗(yàn)方法的不足，是化學(xué)品熱穩(wěn)定性評估的重要理論預(yù)測方法，也是當(dāng)前研究熱點(diǎn)[4]。近年來，已有國內(nèi)外學(xué)者針對有機(jī)過氧化物和硝基化合物的起始放熱溫度和分解熱建立QSPR預(yù)測模型，Prana等[5]基于38種有機(jī)過氧化物的熱穩(wěn)定性實(shí)驗(yàn)數(shù)據(jù)，建立四元和三元線性回歸(MLR)預(yù)測模型，分別預(yù)測To和-ΔHd；Zhang等[1]基于63種硝基芳香族化合物和16種有機(jī)過氧化物，綜合遺傳算法和多元線性回歸(GA-MLR)分別建立硝基芳香族化合物和有機(jī)過氧化物的To預(yù)測模型?？傮w來說，上述研究中模型均取得了較好的預(yù)測效果，但是針對活性化合物的熱穩(wěn)定性的QSPR預(yù)測研究仍較少，且起始放熱溫度和分解熱的預(yù)測仍以線性模型為主，在表征分子結(jié)構(gòu)與熱穩(wěn)定性之間的復(fù)雜關(guān)系方面存在不足。針對上述問題，本文采用遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)，建立QSPR預(yù)測模型，對有機(jī)過氧化物和硝基化合物分子結(jié)構(gòu)及其起始放熱溫度、分解熱的內(nèi)在非線性關(guān)系進(jìn)行探討。

1 數(shù)據(jù)樣本選取

為獲得性能良好的QSPR模型，最標(biāo)準(zhǔn)的數(shù)據(jù)樣本應(yīng)選自同一數(shù)據(jù)源，且在實(shí)驗(yàn)室條件、測試人員和測試方法不變的情況下獲得[6]。通過嚴(yán)格篩選，本文的數(shù)據(jù)樣本均來自同一數(shù)據(jù)源。有機(jī)過氧化物起始放熱溫度和分解熱的實(shí)驗(yàn)數(shù)據(jù)來自Prana等[5]研究中測得的38種有機(jī)過氧化物的實(shí)驗(yàn)數(shù)據(jù)；硝基化合物起始放熱溫度和分解熱的實(shí)驗(yàn)數(shù)據(jù)來自Ando等[7]開展的DSC實(shí)驗(yàn)研究中104種硝基化合物的實(shí)驗(yàn)數(shù)據(jù)。以上數(shù)據(jù)樣本分別用于建立有機(jī)過氧化物和硝基化合物熱穩(wěn)定性預(yù)測模型，其中硝基化合物熱穩(wěn)定性的數(shù)據(jù)樣本比Zhang等研究[1]的范圍更廣。采用隨機(jī)原則，劃分?jǐn)?shù)據(jù)樣本的80%(30種有機(jī)過氧化物和83種硝基化合物的起始放熱溫度和分解熱)作為訓(xùn)練集，用來篩選分子描述符和建立模型；劃分?jǐn)?shù)據(jù)樣本的20%(8種有機(jī)過氧化物和21種硝基化合物的起始放熱溫度和分解熱)作為測試集，用來驗(yàn)證模型。數(shù)據(jù)樣本信息如表1所示。

表1 數(shù)據(jù)樣本信息Table 1 Data sample statistics

2 分子描述符計(jì)算與篩選

本文在計(jì)算分子描述符之前，首先采用HyperChem 8.0軟件[8]繪制有機(jī)過氧化物和硝基化合物的3D分子結(jié)構(gòu)，并采用VB編程，利用DDE接口操控HyperChem的分子力學(xué)(MM+)和量子化學(xué)半經(jīng)驗(yàn)方法(AM1)對繪制的分子結(jié)構(gòu)進(jìn)行批量優(yōu)化，獲得能量最低的穩(wěn)定構(gòu)型，將優(yōu)化后分子結(jié)構(gòu)導(dǎo)入E-Dragon在線計(jì)算工具[9]計(jì)算數(shù)據(jù)樣本中各化合物的分子描述符。

經(jīng)過計(jì)算，數(shù)據(jù)樣本中每個化合物均獲得1 664個分子描述符。為提取預(yù)測模型建模所需的特征變量，需要對大量的分子描述符進(jìn)行降維、篩選。首先利用MATLAB進(jìn)行分子描述符初篩，剔除大量冗余的描述符。初篩的原則包括：1)剔除對所有化合物來說數(shù)值為常數(shù)或近似常數(shù)的描述符；2)對于兩兩相關(guān)系數(shù)大于0.95的描述符，剔除與熱穩(wěn)定性參數(shù)相關(guān)性較差的一方；3)剔除與熱穩(wěn)定性參數(shù)相關(guān)系數(shù)小于0.1的描述符。初篩后，有機(jī)過氧化物與硝基化合物的起始放熱溫度數(shù)據(jù)樣本對應(yīng)的分子結(jié)構(gòu)描述符分別由1 664個減少至477個和396個；分解熱數(shù)據(jù)樣本對應(yīng)的分子結(jié)構(gòu)描述符分別由1 664個減少至480個和479個。

基于初篩的描述符，采用Materials Studio軟件的遺傳函數(shù)算法(GFA)，按照以下步驟篩選描述符：1)設(shè)置描述符篩選目標(biāo)數(shù)量；2)由GFA算法篩選出與熱穩(wěn)定性參數(shù)相關(guān)的對應(yīng)目標(biāo)數(shù)量的分子描述符組合10組，并建立篩選的分子描述符與熱穩(wěn)定性參數(shù)的10個回歸方程；3)比較10個回歸方程的決定系數(shù)R2，R2最接近1的方程的變量即為目標(biāo)數(shù)量下最佳的分子描述符組合；4)重復(fù)步驟1)～3)，獲得目標(biāo)數(shù)量為1～10的最佳分子描述符組合；5)分析不同目標(biāo)數(shù)量下最佳分子描述符組合的R2值隨描述符目標(biāo)數(shù)量的變化規(guī)律，如圖1所示，并依據(jù)Katritzky等[10]提出的“斷點(diǎn)原則”，由“斷點(diǎn)”位置確定最佳的分子描述符數(shù)量，該數(shù)量下的描述符組合即為建模的最優(yōu)描述符，如表2所示。

表2 最優(yōu)分子描述符組合Table 2 Combination of optimal molecular descriptors

圖1 描述符組合的R2值隨數(shù)量的變化規(guī)律Fig.1 Change law of R2 value of descriptors combination with number

3 QSPR預(yù)測模型的建立

BP神經(jīng)網(wǎng)絡(luò)是1種單向傳播的具有3層或3層以上的前向神經(jīng)網(wǎng)絡(luò)，包括輸入層、隱含層和輸出層，上下層之間實(shí)現(xiàn)全連接。若網(wǎng)絡(luò)輸出與期望輸出值不一致，則將其誤差信號反向傳播，并在傳播過程中對加權(quán)系數(shù)不斷修正，使在輸出層節(jié)點(diǎn)上得到的輸出結(jié)果盡可能接近期望輸出值[6]。為優(yōu)化BP神經(jīng)網(wǎng)絡(luò)的預(yù)測效果，采用遺傳算法(GA)優(yōu)化網(wǎng)絡(luò)的權(quán)值和閾值，優(yōu)化后的BP神經(jīng)網(wǎng)絡(luò)能得到更好的預(yù)測輸出結(jié)果。本文基于GFA篩選出的描述符，在MATLAB中實(shí)現(xiàn)BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練與測試，同時(shí)采用GA算法優(yōu)化BP網(wǎng)絡(luò)的權(quán)值和閾值，提升BP模型的性能，最終獲得有機(jī)過氧化物、硝基化合物的起始放熱溫度、分解熱的GA-BP預(yù)測模型。同時(shí)，為對比模型預(yù)測效果，以相同的描述符作為輸入?yún)?shù)，采用偏最小二乘法(PLS)建立線性預(yù)測模型。GA-BP模型和PLS模型預(yù)測結(jié)果分別如圖2、圖3所示，對比2圖中GA-BP模型和PLS模型的離散性可知，GA-BP模型預(yù)測結(jié)果明顯優(yōu)于PLS模型。以有機(jī)過氧化物的分解熱預(yù)測模型為例，列舉部分樣本的預(yù)測結(jié)果，見表3。

圖2 有機(jī)過氧化物熱穩(wěn)定性GA-BP預(yù)測模型和PLS預(yù)測模型預(yù)測結(jié)果對比Fig.2 Comparison on prediction results of GA-BP prediction model and PLS prediction model for thermal stability of organic peroxides

圖3 硝基化合物熱穩(wěn)定性GA-BP預(yù)測模型和PLS預(yù)測模型結(jié)果對比Fig.3 Comparison on prediction results of GA-BP prediction model and PLS prediction model for thermal stability of nitro compounds

表3 部分有機(jī)過氧化物的分子描述符與分解熱數(shù)據(jù)Table 3 Molecular descriptors and heat of decomposition of some organic peroxides

4 模型驗(yàn)證與討論

4.1 模型驗(yàn)證與對比

表4 熱穩(wěn)定性預(yù)測模型驗(yàn)證結(jié)果Table 4 Validation results of thermal stability prediction models

同時(shí)，為驗(yàn)證本文建立的熱穩(wěn)定性參數(shù)預(yù)測的GA-BP模型的應(yīng)用效果，額外選取數(shù)據(jù)樣本集以外的1種硝基化合物，驗(yàn)證預(yù)測方法的合理性。4,4-二硝基二苯二硫醚(CAS：100-32-3)起始放熱溫度預(yù)測值為258.12 ℃，分解熱預(yù)測值為-2019.23 J/g；起始放熱溫度實(shí)驗(yàn)值為248.00 ℃，分解熱實(shí)驗(yàn)值為-1955.24 J/g。通過應(yīng)用效果驗(yàn)證，本文預(yù)測方法獲得熱穩(wěn)定性參數(shù)預(yù)測值與實(shí)驗(yàn)值較為接近。

4.2 GA-BP模型的應(yīng)用域評價(jià)

模型應(yīng)用域的評價(jià)對于預(yù)測模型的應(yīng)用是非常重要的，只有當(dāng)所預(yù)測的化合物處于模型的應(yīng)用域范圍內(nèi)，預(yù)測模型對性質(zhì)的預(yù)測才是有效的[12]。前文中，模型驗(yàn)證與對比證明了GA-BP模型的優(yōu)越性，本文采用Williams圖進(jìn)一步驗(yàn)證模型的應(yīng)用域。

若Williams圖中所有數(shù)據(jù)均在±3倍標(biāo)準(zhǔn)偏差和臨界臂比值(h*)構(gòu)成的范圍之內(nèi)(圖4中3條虛線與縱坐標(biāo)軸圍成區(qū)域)，則說明所有的化合物均在模型的應(yīng)用域范圍之內(nèi)。從圖4中可以看出，4個模型的數(shù)據(jù)均落在±3倍標(biāo)準(zhǔn)偏差之內(nèi)；其中圖4(b)～圖4(d)中少數(shù)數(shù)據(jù)點(diǎn)的臂比值(h)超過了h*，說明對應(yīng)化合物的某些分子結(jié)構(gòu)對于樣本集整體而言比較特殊，預(yù)測結(jié)果是模型外推得來的。雖然部分臂比值超過了臨界值，但是標(biāo)準(zhǔn)化殘差值落在±3倍標(biāo)準(zhǔn)偏差范圍之內(nèi)，有助于使模型更加穩(wěn)定并具有一定的外推能力，因此可以看作是“良性異常點(diǎn)”[13]。綜上，通過應(yīng)用域驗(yàn)證了模型的有效性以及可靠性。

圖4 GA-BP模型的Williams圖Fig.4 Williams diagrams of GA-BP model

4.3 GA-BP神經(jīng)網(wǎng)絡(luò)模型機(jī)理淺析

BP神經(jīng)網(wǎng)絡(luò)是1種“黑箱”模型，因此難以對GA-BP模型進(jìn)行機(jī)理解釋。本文采用描述符重要度衡量方法(MMDI)[14]分析模型描述符對模型的影響程度，淺析影響熱穩(wěn)定性的主要結(jié)構(gòu)因素。MMDI方法的基本原理是，以GA-BP模型的AARD值作為基數(shù)AARD值，依次打亂模型中各分子描述符的數(shù)值順序，利用GA-BP模型對新的數(shù)據(jù)樣本進(jìn)行預(yù)測，新預(yù)測結(jié)果的AARD值與基數(shù)AARD值之差作為各描述符對應(yīng)的AARD差值，比較其大小即可獲得各描述符對模型的影響程度。本文中有機(jī)過氧化物起始放熱溫度GA-BP模型的描述符重要度排序?yàn)椋篗s>H7m>R5v+>L2v>Mor08m，其中Ms和H7m對模型的影響最大，分別是組成描述符和GETAWAY描述符，反映了分子的平均電拓?fù)錉顟B(tài)和中心原子坐標(biāo)；有機(jī)過氧化物分解熱GA-BP模型的描述符重要度排序?yàn)椋篟7e+>Hy>O-063>Mor08m>AAC>Mor32m，其中R7e+和Hy對模型的影響最大，分別是GETAWAY描述符和分子特性描述符，反映了分子的中心原子坐標(biāo)和親水性能；硝基化合物起始放熱溫度GA-BP模型的描述符重要度排序?yàn)椋篨1Av>RDF010m>HATS7p>H4m>R7e>EEig12r，其中X1Av和RDF010m對模型的影響最大，分別是連接性指數(shù)和RDF描述符，反映了分子圖中的原子頂點(diǎn)度和原子間距離；硝基化合物分解熱GA-BP模型的描述符重要度排序?yàn)椋篗or30u>Mor03e>Mor24m>GGI2>ATS3m>Mor21v，其中Mor30u和Mor03e對模型的影響最大，兩者均是3D-MoRSE描述符，反映了原子的3D排列。綜上，利用MMDI方法明確各模型中影響化合物熱穩(wěn)定性參數(shù)的主要結(jié)構(gòu)因素。

5 結(jié)論

1)采用遺傳算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò)參數(shù)，建立預(yù)測有機(jī)過氧化物和硝基化合物的起始放熱溫度與分解熱的GA-BP模型，其擬合能力、穩(wěn)定性、預(yù)測能力以及應(yīng)用域得到驗(yàn)證。其中，硝基化合物熱穩(wěn)定性預(yù)測模型建立在比前人的研究更廣泛的數(shù)據(jù)樣本基礎(chǔ)之上，具有更強(qiáng)的代表性和適用性。

2)本文構(gòu)建的GA-BP模型在熱穩(wěn)定性預(yù)測方面更具優(yōu)勢，說明非線性模型更能表征分子結(jié)構(gòu)與熱穩(wěn)定性之間的復(fù)雜關(guān)系，從一定程度上證明熱穩(wěn)定性參數(shù)與分子結(jié)構(gòu)之間存在非線性關(guān)系。

3)采用MMDI方法對GA-BP模型的機(jī)理進(jìn)行淺析，得出各模型中影響化合物熱穩(wěn)定性參數(shù)的主要結(jié)構(gòu)因素。