齊書平，王文龍，張磊，都健

（大連理工大學(xué)化工學(xué)院，化工系統(tǒng)工程研究所，遼寧 大連 116024）

引 言

金屬離子（Mn+）和有機(jī)配體（L）的絡(luò)合過程是配位化學(xué)中應(yīng)用最廣泛的過程之一[1]，涉及催化[2-4]、濕法冶金[5]、醫(yī)療[6]、生物醫(yī)學(xué)與環(huán)境監(jiān)測(cè)[7]等多種領(lǐng)域。該過程在這些領(lǐng)域的應(yīng)用都離不開配合物的穩(wěn)定性問題，而配合物的穩(wěn)定常數(shù)是衡量溶液中有機(jī)配體與金屬離子間結(jié)合穩(wěn)定性的重要指標(biāo)。目前，獲取配合物穩(wěn)定常數(shù)的傳統(tǒng)方法是通過實(shí)驗(yàn)方法測(cè)定或通過分子模擬計(jì)算。實(shí)驗(yàn)方法包括電位法、分光光度法、色譜法等，這類方法通常煩瑣復(fù)雜且耗時(shí)，特別是對(duì)大環(huán)化合物需要幾天或幾周的時(shí)間才能達(dá)到平衡。傳統(tǒng)的分子模擬計(jì)算方法主要是利用密度泛函理論（density functional theory,DFT）或分子動(dòng)力學(xué)（molecular dynamics, MD）進(jìn)行幾何優(yōu)化和能量計(jì)算，通常費(fèi)時(shí)費(fèi)力且很難達(dá)到模型精度要求。

由于上述原因，利用實(shí)驗(yàn)數(shù)據(jù)建立回歸模型預(yù)測(cè)金屬離子和配體的穩(wěn)定常數(shù)成為一種很有潛力的方法。此前，文獻(xiàn)中報(bào)道過許多試圖預(yù)測(cè)配體和金屬離子間穩(wěn)定常數(shù)的機(jī)器學(xué)習(xí)模型，涉及的機(jī)器學(xué)習(xí)方法包括多元線性回歸（multiple linear regression, MLR）、支 持 向 量 機(jī)（support vector machine, SVM）、K-近鄰（K-nearest neighbor, KNN）等[8-18]。這些模型中的絕大多數(shù)都是在某種確定的實(shí)驗(yàn)條件下，局限于構(gòu)建單一金屬離子與不同配體的定量結(jié)構(gòu)性質(zhì)關(guān)系（quantitative structure-property relationship, QSPR）模型[8-16]或某一類有機(jī)配體及其衍生物（如席夫堿及其衍生物）與不同金屬的QSPR模型[17]，配體使用的大都是子結(jié)構(gòu)片段（substructural molecular fragments，SMF）描述符，每個(gè)模型所用的數(shù)據(jù)集很?。◣资缴习伲?，配體的類型十分有限，也沒有考慮金屬離子和配體間的相互作用。針對(duì)這些模型缺陷，Chaube 等[18]提出了考慮不同實(shí)驗(yàn)條件的預(yù)測(cè)鑭系金屬離子-配體間穩(wěn)定常數(shù)的機(jī)器學(xué)習(xí)模型，使用的描述符為物理化學(xué)描述符和拓?fù)涿枋龇?，該工作的不足之處在于模型只適用于性質(zhì)相近的鑭系金屬離子，且沒有考慮多配位的情況。綜上所述，目前提出的模型僅針對(duì)有限類別的配體或金屬離子而構(gòu)建，且由于這些限制，這些模型的錯(cuò)誤率較高，泛化能力較弱，限制了它們?cè)诓町惡艽蟮呐潴w化學(xué)空間中的應(yīng)用。

此外，傳統(tǒng)的QSPR 模型通?；诔砂偕锨€(gè)已知描述符的組合來提取分子特征[18]，缺少對(duì)配合物的針對(duì)性。至于配合物穩(wěn)定常數(shù)lgK的大小，主要與有機(jī)配體和金屬離子形成的配合物的結(jié)構(gòu)有關(guān)。而不同的金屬離子和配體由于配位數(shù)和配位原子數(shù)的不同，難以表示整個(gè)配合物的分子圖。因此為了便于建模，把每種配合物的表示分成了兩步，第一步通過有機(jī)配體的分子圖及原子和鍵的性質(zhì)，準(zhǔn)確地描述有機(jī)配體（包括順反異構(gòu)、手性等），模型認(rèn)識(shí)了這個(gè)有機(jī)分子的特征，就有了預(yù)測(cè)的基礎(chǔ)；第二步將多頭圖注意力網(wǎng)絡(luò)（multi-head GAT）[19]提取的分子特征拼接金屬離子的獨(dú)熱編碼，考慮有機(jī)配體和金屬離子的相互作用。

由于穩(wěn)定常數(shù)還受溫度、介質(zhì)、離子強(qiáng)度等實(shí)驗(yàn)條件的影響[20]，因此本研究除了考慮有機(jī)配體和金屬離子的相互作用，還編碼了穩(wěn)定常數(shù)測(cè)定時(shí)所處的實(shí)驗(yàn)條件。本文提出的深度學(xué)習(xí)模型以配體的簡化線性輸入規(guī)范（simplified molecular input line entry system, SMILES）轉(zhuǎn)換得到的分子屬性圖作為輸入，經(jīng)圖注意力特征提取層，識(shí)別對(duì)配合物穩(wěn)定常數(shù)有重要貢獻(xiàn)的節(jié)點(diǎn)特征，得到分子圖的特征表示并將提取的分子特征與使用獨(dú)熱編碼的金屬離子和實(shí)驗(yàn)條件進(jìn)行拼接后送入全連接層，預(yù)測(cè)金屬離子與有機(jī)配體的穩(wěn)定常數(shù)。

1 基于深度學(xué)習(xí)的配合物穩(wěn)定常數(shù)預(yù)測(cè)模型框架

本文提出的模型框架由數(shù)據(jù)準(zhǔn)備和模型訓(xùn)練兩部分組成，如圖1 所示。首先，從IUPAC（the international union of pure and applied chemistry,IUPAC）穩(wěn)定常數(shù)數(shù)據(jù)庫[21]的子庫（mini stability constant database）中獲取建立模型所需的實(shí)驗(yàn)數(shù)據(jù)，利用 DGL[22]（deep graph library，DGL）對(duì)有機(jī)配體生成分子屬性圖，其中各原子的屬性值由RDKit[23]計(jì)算得到。其次，使用獨(dú)熱編碼對(duì)實(shí)驗(yàn)數(shù)據(jù)中的金屬離子和實(shí)驗(yàn)條件部分進(jìn)行編碼，獲得條件特征向量。最后，將分子屬性圖送入建立的多頭圖注意力網(wǎng)絡(luò)進(jìn)行分子特征的抽提，將所得特征與條件特征向量進(jìn)行拼接，送入全連接層進(jìn)行穩(wěn)定常數(shù)的預(yù)測(cè)。

圖1 預(yù)測(cè)金屬離子-有機(jī)配體穩(wěn)定常數(shù)的深度學(xué)習(xí)模型框架Fig.1 A deep learning model framework for predicting the stability constants of metal ions-organic ligands

1.1 數(shù)據(jù)準(zhǔn)備

本文使用IUPAC SCDatabase 的子庫mini-SCDatabase作為訓(xùn)練所需的實(shí)驗(yàn)數(shù)據(jù)來源。數(shù)據(jù)的提取過程為：首先，根據(jù)所研究的金屬離子的種類選定金屬離子下載一個(gè)文本文件，文件中每條數(shù)據(jù)包括配體名稱（俗名和配體的大致分類，部分?jǐn)?shù)據(jù)包含CAS 號(hào)，分子式）、金屬離子、介質(zhì)的種類和濃度、實(shí)驗(yàn)溫度、穩(wěn)定常數(shù)；然后，編寫python腳本自動(dòng)整理數(shù)據(jù)，刪去不含碳和碳數(shù)大于20 的配體，從文本文件中提取建模所需的信息并導(dǎo)出到excel 表格中，包括配體名稱、CAS 號(hào)、實(shí)驗(yàn)條件、單配位 （Mn+∶L = 1∶1）穩(wěn)定常數(shù)lgK1和雙配位（Mn+∶L=1∶2）穩(wěn)定常數(shù)lgβ2；最后，根據(jù)配體的CAS 號(hào)在PubChem 數(shù)據(jù)庫中爬取對(duì)應(yīng)的SMILES，對(duì)CAS 號(hào)錯(cuò)誤或缺失的配體，根據(jù)分子式和俗名手動(dòng)寫出SMILES。對(duì)于本文所涉及的29種金屬離子，經(jīng)篩選和整理后得到7127條實(shí)驗(yàn)數(shù)據(jù)，每條數(shù)據(jù)由有機(jī)配體、金屬離子、實(shí)驗(yàn)條件（實(shí)驗(yàn)介質(zhì)、支持電解質(zhì)的離子強(qiáng)度或混合溶劑的體積比、配位狀態(tài)以及實(shí)驗(yàn)溫度）和穩(wěn)定常數(shù)（配位比為Mn+∶L = 1∶1 的lgK1和配位比為Mn+∶L=1∶2的lgβ2）組成。數(shù)據(jù)集中絡(luò)合物的金屬離子和穩(wěn)定常數(shù)的數(shù)據(jù)點(diǎn)分布如圖2所示。從圖2（a）中可以看出除少數(shù)金屬離子外，大多數(shù)金屬離子的數(shù)據(jù)量都大于100，從圖2（b）中可以看出穩(wěn)定常數(shù)大都處在（1，15）的取值范圍中。

圖2 數(shù)據(jù)分布圖Fig.2 The distribution of data points

本文使用DGL 對(duì)所有實(shí)驗(yàn)數(shù)據(jù)中涉及的1371個(gè)有機(jī)分子生成分子屬性圖，屬性圖中原子及鍵的特征與編碼方式如表1所示。分子屬性圖不僅考慮了原子層面的特征（如原子種類、形式電荷、原子雜化方式等）和分子層面的特征（如Gasteiger 電荷、Crippen 摩爾折射率、拓?fù)錁O性表面積TPSA 等），還考慮了各種化學(xué)鍵的特征（如鍵類型、鍵方向、鍵數(shù)量等）。原子特征的選取參考了Coley 等[24]在分子屬性圖方面的工作，并在原子特征矩陣中引入鍵特征參數(shù)，既考慮了分子的空間信息，也能更有效地表示原子所處的化學(xué)環(huán)境。

表1 分子屬性圖中原子與鍵的特征及編碼方式Table 1 Features and coding methods of atoms and bonds in molecular attribute diagrams

1.2 模型訓(xùn)練

本文在模型訓(xùn)練過程中采用的是留出驗(yàn)證，即把數(shù)據(jù)劃分為訓(xùn)練集、驗(yàn)證集和測(cè)試集，劃分比例為90∶5∶5，其中驗(yàn)證集用來選擇模型的超參數(shù)（如學(xué)習(xí)率、層數(shù)等），測(cè)試集用來評(píng)估模型性能。此外，本文在構(gòu)建訓(xùn)練集時(shí)充分考慮了有機(jī)配體的多樣性，具體做法為：首先，在隨機(jī)打亂的數(shù)據(jù)池中，順序依次取出所有第一次出現(xiàn)的有機(jī)配體SMILES所對(duì)應(yīng)的數(shù)據(jù)，共1371 條數(shù)據(jù)，然后在剩余的數(shù)據(jù)中隨機(jī)抽取5044條數(shù)據(jù)共同組成訓(xùn)練集。最后，隨機(jī)抽取驗(yàn)證集數(shù)據(jù)356條，剩下的則放入測(cè)試集。

深度學(xué)習(xí)模型主要由多頭自注意力特征提取層和全連接層兩部分組成。其中GAT 層作為第一層次讀入分子屬性圖，將其按照?qǐng)D神經(jīng)網(wǎng)絡(luò)節(jié)點(diǎn)信息更新規(guī)則對(duì)分子屬性圖進(jìn)行迭代，使得模型認(rèn)識(shí)到這個(gè)有機(jī)配體；全連接層作為第二層次讀入經(jīng)過抽提的分子屬性圖、金屬離子以及實(shí)驗(yàn)條件，目的是讓在模型認(rèn)識(shí)到這個(gè)有機(jī)配體的情況下再進(jìn)一步獲取其他必要信息，以便進(jìn)行l(wèi)gK的預(yù)測(cè)。之所以將分子屬性圖單獨(dú)作為GAT 層的輸入而不結(jié)合其他信息（如金屬離子或某些實(shí)驗(yàn)條件），是因?yàn)橄Ｍ⒌哪Ｐ拖韧ㄟ^屬性圖了解到分子的信息后大致確定一個(gè)預(yù)測(cè)范圍（也就是配體絡(luò)合能力的強(qiáng)弱），再通過后續(xù)輸入的其他必要信息縮小預(yù)測(cè)范圍，直至實(shí)現(xiàn)精準(zhǔn)預(yù)測(cè)。綜上，通過建立雙層次的圖神經(jīng)網(wǎng)絡(luò)，模型具有了學(xué)習(xí)預(yù)測(cè)lgK所需要的物理知識(shí)的能力。多頭自注意力特征提取層是通過multi-head GAT 實(shí)現(xiàn)的，它背后的關(guān)鍵思想是注意力機(jī)制，可以將不同的權(quán)重分配給鄰域中的不同節(jié)點(diǎn)，特征提取過程后得到的信息對(duì)穩(wěn)定常數(shù)預(yù)測(cè)更為有效。計(jì)算如式（1）所示。

本文中多頭圖注意力機(jī)制的實(shí)現(xiàn)使用了DGL框架（后端為PyTorch[25]）。DGL 是一款易于使用的開源圖神經(jīng)網(wǎng)絡(luò)計(jì)算框架，它通過與主流的深度學(xué)習(xí)框架（目前支持Pytorch、TensorFlow[26]等）集成，實(shí)現(xiàn)了從傳統(tǒng)的張量運(yùn)算到圖運(yùn)算的自由轉(zhuǎn)換。本文使用DGL 庫內(nèi)置的GATConv 模塊搭建8 層的多頭圖注意力網(wǎng)絡(luò)對(duì)分子屬性圖進(jìn)行特征提取，將提取的分子特征與使用獨(dú)熱編碼的金屬離子與實(shí)驗(yàn)條件進(jìn)行拼接，送入3 層的全連接層進(jìn)行配合物穩(wěn)定常數(shù)的預(yù)測(cè)。全連接層的大小分別為128、32 和1。訓(xùn)練過程中的批次大小設(shè)置為128，使用初始學(xué)習(xí)率為0.0003 的Adam 優(yōu)化器，在訓(xùn)練1500 代后停止。

2 訓(xùn)練結(jié)果與案例分析

2.1 模型訓(xùn)練結(jié)果

由于本文的穩(wěn)定常數(shù)預(yù)測(cè)是一項(xiàng)回歸任務(wù)，故選用均方誤差（mean square error, MSE）作為損失函數(shù)，如式（2）所示。

式中，Pi和Yi分別是第i個(gè)數(shù)據(jù)的預(yù)測(cè)值和實(shí)驗(yàn)值；n為樣本總數(shù)。模型的訓(xùn)練結(jié)果如圖3 所示，圖3（a）為損失函數(shù)下降曲線，其中紅色曲線表示訓(xùn)練集的損失函數(shù)曲線，藍(lán)色曲線表示驗(yàn)證集的損失函數(shù)曲線。從圖3（a）中可以看出在第1500 代時(shí)訓(xùn)練集與驗(yàn)證集損失函數(shù)趨于穩(wěn)定，故停止訓(xùn)練。模型的預(yù)測(cè)能力可用回歸決定系數(shù)R2和均方根誤差（root mean square error, RMSE）進(jìn)行評(píng)價(jià)。R2計(jì)算如式（3）所示。RMSE計(jì)算如式（4）所示。

圖3 模型訓(xùn)練結(jié)果Fig.3 Model training results

圖3（b）、（c）、（d）分別為訓(xùn)練集、驗(yàn)證集和測(cè)試集的實(shí)驗(yàn)值和預(yù)測(cè)值之間的結(jié)果差異圖，其中位于兩條虛線之內(nèi)的點(diǎn)表示其與實(shí)驗(yàn)值的絕對(duì)誤差不超過2。訓(xùn)練集、驗(yàn)證集和測(cè)試集的R2與RMSE 分別為0.995/0.446、0.956/1.316 和0.956/1.251。數(shù)據(jù)集中共有88 個(gè)點(diǎn)在誤差線外，其中訓(xùn)練集中有29個(gè)在誤差線外的數(shù)據(jù)點(diǎn)，占訓(xùn)練集總樣本數(shù)的0.45%；驗(yàn)證集中有30 個(gè)點(diǎn)在誤差線外，涉及13 種金屬離子，占驗(yàn)證集樣本總數(shù)的8.42%；測(cè)試集中有29 個(gè)點(diǎn)在誤差線外，涉及14 種金屬離子，占測(cè)試集樣本總數(shù)的8.14%。在這88 個(gè)點(diǎn)中，穩(wěn)定常數(shù)的實(shí)驗(yàn)值大于15 的點(diǎn)有41 個(gè)，占這個(gè)區(qū)間所有數(shù)據(jù)點(diǎn)的3.90%，實(shí)驗(yàn)值小于15 的點(diǎn)占同區(qū)間所有數(shù)據(jù)點(diǎn)的0.77%，這是由穩(wěn)定常數(shù)的數(shù)據(jù)分布導(dǎo)致的。

2.2 案例分析

近年來，研究者們廣泛使用DFT 計(jì)算[27-28]、MD模擬[29-30]來獲得配合物在溶液中的絡(luò)合自由能ΔG，絡(luò)合自由能和穩(wěn)定常數(shù)間的轉(zhuǎn)化如式（5）所示[27]。

式中，R為熱力學(xué)常數(shù)，R= 8.314 J·mol-1·K-1；T為熱力學(xué)溫度。這些方法雖然對(duì)比黑箱的深度學(xué)習(xí)模型解釋性更強(qiáng)，但想要得到可靠的結(jié)果，需要選擇合適的基組或力場(chǎng)，并設(shè)置超參數(shù)等，對(duì)不同的體系計(jì)算方法也不同，因此找到合適的計(jì)算方法是一個(gè)難題。以DFT 計(jì)算為例，定量計(jì)算配合物在溶液中的絡(luò)合自由能需要使用計(jì)算步驟煩瑣的熱力學(xué)循環(huán)。此外，即使使用精度較高的基組（如PBE/TPZ），計(jì)算結(jié)果往往也不能盡如人意，大多情況下只能定性比較[27]。

因此，本節(jié)使用前文建立的深度學(xué)習(xí)模型來預(yù)測(cè)文獻(xiàn)中螯合物的穩(wěn)定常數(shù)，表2 為模型預(yù)測(cè)結(jié)果和文獻(xiàn)中DFT 計(jì)算結(jié)果的對(duì)比，結(jié)果展示了模型強(qiáng)大的預(yù)測(cè)能力，所選用的配體和金屬離子對(duì)在數(shù)據(jù)集中沒有出現(xiàn)過。預(yù)測(cè)結(jié)果表明，本文建立的深度學(xué)習(xí)模型計(jì)算成本更低，泛化程度更高，更容易建立與應(yīng)用，且在大多數(shù)情況下模型預(yù)測(cè)結(jié)果更為準(zhǔn)確。表2所示的預(yù)測(cè)情況中，第三種配體與Ni2+的單配位穩(wěn)定常數(shù)預(yù)測(cè)誤差最大，絕對(duì)誤差為4.71，這是因?yàn)镾CDatabase 數(shù)據(jù)庫不再提供維護(hù)，所以本文所用為其子數(shù)據(jù)庫mini-SCDatabase，數(shù)據(jù)量不夠豐富，穩(wěn)定常數(shù)較大的數(shù)據(jù)點(diǎn)較少，存在數(shù)據(jù)不平衡的問題。以單配位態(tài)的Ni2+為例，數(shù)據(jù)集中僅有4%的數(shù)據(jù)lgK1＞17，模型學(xué)習(xí)到螯合能力強(qiáng)的配體的特征較少，以致于模型對(duì)螯合能力太強(qiáng)的配體預(yù)測(cè)偏差較大。

表2 模型預(yù)測(cè)結(jié)果和DFT計(jì)算結(jié)果與實(shí)驗(yàn)值的比較Table 2 Comparison of model prediction results and DFT calculation results with experimental values

3 結(jié) 論

本文利用基于多頭圖注意力網(wǎng)絡(luò)對(duì)29 種金屬離子和1371 個(gè)有機(jī)配體構(gòu)建配合物穩(wěn)定常數(shù)預(yù)測(cè)模型，該模型綜合考慮了多種影響因素，可對(duì)多種配體、多種金屬離子、多種實(shí)驗(yàn)條件進(jìn)行穩(wěn)定常數(shù)的高通量預(yù)測(cè)。本文的主要貢獻(xiàn)如下：（1）創(chuàng)新性地提出了一種配合物的表征新范式，并將multihead GAT 用于配合物的性質(zhì)預(yù)測(cè)；（2）本文模型涉及的金屬離子和有機(jī)配體的多樣性遠(yuǎn)超以往研究，通過建立一個(gè)模型完成了對(duì)過渡金屬離子（Sc3+,V2+, V3+, Mn2+, Fe2+, Fe3+, Co2+, Ni2+）、后過渡金屬離子（Al3+,Ga3+, In3+, Tl3+, Pb2+, Tb3+, Bi3+）、鑭系和錒系金屬離子（La3+, Ce3+, Pr3+, Nd3+, Sm3+, Eu3+, Gd3+, Tb3+, Dy3+,Ho3+, Er3+, Tm3+, Yb3+, Th4+）與有機(jī)配體的穩(wěn)定常數(shù)預(yù)測(cè)，模型適用范圍廣泛。

模型首先基于實(shí)驗(yàn)數(shù)據(jù)涉及的有機(jī)配體生成分子屬性圖，然后使用獨(dú)熱編碼對(duì)金屬離子和實(shí)驗(yàn)條件進(jìn)行編碼，最后將分子屬性圖送入特征提取層進(jìn)行特征提取，將提取出的特征拼接上金屬離子和實(shí)驗(yàn)條件的特征編碼，送入全連接層預(yù)測(cè)穩(wěn)定常數(shù)。模型在測(cè)試集上的RMSE 為1.251。與DFT 計(jì)算結(jié)果相比，本文提出的模型更為高效。而且值得注意的是，與僅針對(duì)螯合劑建模的模型不同，本文中參與建模的有機(jī)配體不僅有螯合劑（多齒配體），還有非螯合劑（單齒配體），這在很大程度上拓展了模型的可用范圍，可用于對(duì)特異金屬螯合劑的高通量篩選。此外，在本模型中使用多頭圖自注意力網(wǎng)絡(luò)直接處理有機(jī)配體的分子圖，可自動(dòng)提取與所預(yù)測(cè)性質(zhì)相關(guān)的分子特征，因此相信該方法還可以應(yīng)用于穩(wěn)定常數(shù)以外的其他配合物性質(zhì)預(yù)測(cè)（如磁性、顏色等）。