基于深度學習的金屬離子-有機配體配位穩定常數的預測

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

（大連理工大學化工學院，化工系統工程研究所，遼寧 大連 116024）

引 言

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

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

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

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

1 基于深度學習的配合物穩定常數預測模型框架

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

圖1 預測金屬離子-有機配體穩定常數的深度學習模型框架Fig.1 A deep learning model framework for predicting the stability constants of metal ions-organic ligands

1.1 數據準備

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

圖2 數據分布圖Fig.2 The distribution of data points

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

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

1.2 模型訓練

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

深度學習模型主要由多頭自注意力特征提取層和全連接層兩部分組成。其中GAT 層作為第一層次讀入分子屬性圖，將其按照圖神經網絡節點信息更新規則對分子屬性圖進行迭代，使得模型認識到這個有機配體；全連接層作為第二層次讀入經過抽提的分子屬性圖、金屬離子以及實驗條件，目的是讓在模型認識到這個有機配體的情況下再進一步獲取其他必要信息，以便進行lgK的預測。之所以將分子屬性圖單獨作為GAT 層的輸入而不結合其他信息（如金屬離子或某些實驗條件），是因為希望所建立的模型先通過屬性圖了解到分子的信息后大致確定一個預測范圍（也就是配體絡合能力的強弱），再通過后續輸入的其他必要信息縮小預測范圍，直至實現精準預測。綜上，通過建立雙層次的圖神經網絡，模型具有了學習預測lgK所需要的物理知識的能力。多頭自注意力特征提取層是通過multi-head GAT 實現的，它背后的關鍵思想是注意力機制，可以將不同的權重分配給鄰域中的不同節點，特征提取過程后得到的信息對穩定常數預測更為有效。計算如式（1）所示。

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

2 訓練結果與案例分析

2.1 模型訓練結果

由于本文的穩定常數預測是一項回歸任務，故選用均方誤差（mean square error, MSE）作為損失函數，如式（2）所示。

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

圖3 模型訓練結果Fig.3 Model training results

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

2.2 案例分析

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

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

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

表2 模型預測結果和DFT計算結果與實驗值的比較Table 2 Comparison of model prediction results and DFT calculation results with experimental values

3 結 論

本文利用基于多頭圖注意力網絡對29 種金屬離子和1371 個有機配體構建配合物穩定常數預測模型，該模型綜合考慮了多種影響因素，可對多種配體、多種金屬離子、多種實驗條件進行穩定常數的高通量預測。本文的主要貢獻如下：（1）創新性地提出了一種配合物的表征新范式，并將multihead GAT 用于配合物的性質預測；（2）本文模型涉及的金屬離子和有機配體的多樣性遠超以往研究，通過建立一個模型完成了對過渡金屬離子（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+）與有機配體的穩定常數預測，模型適用范圍廣泛。

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