基于圖卷積的離子液體CO2溶解度可解釋性預測

張茜茜,陳 平

(1.中北大學 信息與通信工程學院,山西 太原 030051;2.信息探測與處理山西省重點實驗室(中北大學),山西 太原 030051)

0 引 言

離子液體是一種理想的綠色CO2吸收劑[1-2],具有化學穩定性好、腐蝕性低、從中分離CO2需要更少的能量、可供選擇的種類比較多等優點[3-5]。CO2在離子液體中的溶解度值是重要信息,可以為篩選、設計滿足工業需求的離子液體提供重要指導[6]。目前,獲取離子液體中CO2溶解度的主要方法包括實驗測量和建模。而通過實驗測量大量離子液體的CO2溶解度并不現實,實驗過程不僅耗時費力且成本高,故借助模型實現準確的預測具有重要意義?；跓崃W或動力學方法構建的模型難以處理大量的性質數據,同時預測準確度不夠高、通用性差。因此,根據現有的實驗數據開發機器學習模型實現準確、快速、高效的離子液體的CO2溶解度預測是十分必要的。

構建預測模型過程中離子液體的表征方式選取尤為重要,當前的模型構建選擇的分子表征方式主要包括經驗描述符、簡單熱物理參數、基團信息以及包含多種物理化學性質值、拓撲性質的分子描述符。

經驗描述符需要化學背景以及相關先驗知識。簡單熱物理參數多借助現有的數據集,數據集中包含離子液體的偏心因子、臨界性質(臨界溫度與臨界壓力)以及分子量等。Moosanezhad-Kermani H等采用GMDH分組數據處理方法,提出了一種數學關聯形式的白盒模型,模型輸入參數為傳統的熱力學參數[7]?；鶊F信息是指將離子液體手動分解為若干基團的組合,Song Z等使用基團信息表征離子液體,使用人工神經網絡(Artificial Neural Network,ANN)和支持向量機(Support Vector Machines,SVM)構建模型預測離子液體中CO2溶解度[8]。這兩種表征方式都受限于數據集,對新的離子液體應用困難。Aghaie M等引入了幾何、拓撲、信息指數、基于3D矩陣的描述符和量子化學特征等作為輸入參數,用決策樹(Decision Tree,DT)與隨機森林(Random Forest,RF),并與熱力學參數模型進行對比,結果證明相比于熱力學參數,結構性質與離子液體的CO2溶解度有更強的相關性[9]。然而大量的分子描述符需要借助專業的軟件獲取,計算成本高。分子指紋作為容易獲取的表征方式,常見于藥物分子的虛擬篩選模型中,分子指紋是由結構信息編碼而成的高維、稀疏向量,這一過程會造成分子的結構信息表征能力受限,且高維的向量在模型訓練中易過擬合。

而圖卷積神經網絡(Graph Convolutional Neural networks,GCN)進行離子液體的結構特征學習可以解決上述分子表征方式存在的問題。GCN可以從分子圖中捕獲結構特征信息,已經被應用于生物化學領域的活性分析、性能預測等任務中[10],在這些領域的成功應用證明GCN能夠有效學習分子圖的結構數據,挖掘化合物分子圖中的結構信息實現與分子性質的關聯。這些任務與離子液體中CO2溶解度預測任務的共性是要實現結構性質與性質的關聯,說明通過GCN實現離子液體的預測任務同樣是可行的。

該文以離子液體中的CO2溶解度預測為研究目的,提出了基于注意力圖卷積神經網絡的全流程建模方法,主要貢獻概括如下:

(1)在圖嵌入模塊引入注意力機制,使模型關注原子的非局部特征信息。

(2)將圖卷積神經網絡與XGBoost有效組合為預測模型,解決了預測模型容易過擬合的問題。

(3)借助注意力池化層與SHAP對預測模型進行解釋,能夠篩選出對于CO2溶解度預測重要的離子液體結構片段。

1 CO2溶解度預測模型

1.1 總體框架流程

整體流程見圖1。首先將離子液體的分子式轉化為GCN所需要的圖數據,原子的特征向量經過圖卷積后生成新的原子特征向量,然后經過池化轉換為整個分子圖的特征向量,也就是圖嵌入網絡生成的離子液體的圖嵌入向量。

圖1 二氧化碳溶解度預測模型構建及應用流程

1.2 圖數據獲取

要對分子進行表征,通過軟件ChemDraw[11](大多數化學結構軟件都支持)獲取陽離子、陰離子的簡化分子線性輸入規范(Simplified Molecular Input Line Entry System,SMILES)。SMILES字符串是一種類似語言的描述,可以通過利用化學信息學工具方便地表達分子二維和三維結構特征,從而用于表征分子特征。然后使用Python中的化學信息學工具包RDKit將SMILES字符串轉換為所需要的圖數據與Morgan分子指紋。RDKit獲取的離子液體的SMILES表達式所包含的圖數據包含節點特征(如原子的元素類型、價態等)與結構特征(如鍵的種類、鍵是否在環中等)。

1.3 圖卷積

訓練過程旨在學習這些特征來獲得分子的整體特征信息,最終表示為一個固定長度的向量[12]。假設離子液體中包含n個原子,包含所有原子的特征向量組成的n×d維特征矩陣X,以及原子之間的n×n維的鄰接矩陣A,X與A也就是GCN模型的輸入,H(l)表示所有原子在第l層的特征向量矩陣,H(l+1)表示所有原子經過一次卷積操作之后的特征向量矩陣,GCN層與層之間的傳播方式如式1所示:

(1)

在每一層中,將原子的特征與鄰域原子的特征信息聚合,生成新的原子表示,對原子的特征向量進行更新,作為下一層的輸入。隨著網絡層數的增加,逐漸獲取所有的原子特征信息與結構信息生成新的原子特征向量。

1.4 注意力池化層

為了實現性質預測任務,希望從離子液體的分子圖生成圖級表示,大多數先前的研究用簡單的最大值池化、平均池化等[14-16]來聚合分子圖中所有的原子特征來生成整個分子圖的表示,但可能會導致不必要的合并、不完整的解釋等問題,因此使用注意力池化層作為傳統池化方式的替代。目標是將分子圖編碼成固定大小的嵌入矩陣,同時最大化原子特征向量所隱含的信息,學習原子的全局信息。

注意力機制允許方法關注神經網絡的任務相關部分。對具有序列結構數據的任務應用注意力機制已成為一種慣例,以使模型能夠集中注意輸入的最相關部分并實現更好的預測[17]。

通過注意力池化層生成注意力系數矩陣,使用注意力機制,將從圖卷積模塊學習的原子特征向量作為輸入來輸出權重向量。注意力系數也就是輸出的權重,代表不同特征的重要性,權重越大,特征越重要,根據注意力系數大小篩選出對于溶解度預測相對重要的特征,這里的特征是分子圖中的化學結構片段,經過圖卷積模塊的卷積、合并,每個原子的特征向量不僅包含原子本身的信息,還包含原子的鄰域信息。

經過圖卷積模塊,離子液體的結構特征被編碼為每個原子的特征向量組合,此時的圖表示為G,αvn為原子的特征向量,G的每一行表示一個原子的特征向量,如下式:

G=(αv1,αv2,…,αvn)

(2)

通過注意力池化:

A=softmax(β2tanh(β1GT))

(3)

其中,β1,β2表示權重矩陣,則最終輸出的圖的表達形式變為:

G'=AG

(4)

加入注意力池化層的圖卷積神經網絡(Attention Pool Artificial Neural Network,APGCN)所構成的圖嵌入生成模塊流程如圖2所示,通過圖卷積神經網絡獲取分子圖中每個原子的特征向量。

圖2 APGCN模塊示意圖

1.5 預測模塊算法選擇

目前,離子液體的CO2溶解度預測模型構建中ANN,SVM,RT模型等被廣泛應用。在輸入數據為分子描述符、高維向量時機器學習算法綜合表現更好,訓練效率高、不易過擬合[18]。通過6折交叉驗證選擇最合適的機器學習算法進行建模,在數據集上對機器學習算法XGBoost(e-Xtreme Gradient Boosting)[19]、梯度提升決策樹(Gradient Boosting Decision Tree,GBDT)[20]、LightGBM(Light Gradient Boosting Machine)[21]進行6折交叉驗證實驗對比,評估指標選擇決定系數(R2)、平均絕對誤差(MAE)、均方誤差(MSE)。表1中結果表明,經過驗證,在選定的數據集上使用XGBoost構建的模型表現更好,更能避免過擬合,且相比于其他機器學習算法有較高的準確性,同時計算成本低、效率高,故選擇XGBoost算法來構建CO2的預測模型。

表1 圖嵌入和Morgan指紋交叉驗證對比

此預測任務的輸入不僅有圖數據,還有對于離子液體的CO2溶解度不可忽略的重要熱力學參數T,P,因此模型要實現基于圖數據與數值型混合數據的預測。綜上所述,充分考慮到離子液體的CO2預測模型需要輸入數據的特點,經過圖嵌入生成網絡將圖數據轉化為一個低維向量,將其與對應的T,P條件結合作為XGBoost預測模型的輸入,輸出最終的CO2溶解度預測值。

XGBoost是基于決策樹的集成模型,每次訓練時對數據集進行采樣,隨機抽取一定數量的樣本作為每棵樹的根節點。每個樣本所包含的分子指紋所表示的離子液體結構特征、溫度T、壓力P即為待分裂的特征。每一棵樹生成過程如圖3所示,溫度T、壓力P、結構特征都是影響離子液體CO2溶解度的重要因素,從根節點開始分裂,圖3中的樹根節點先對壓力P進行判斷,中間節點有溫度T、分子的結構特征向量,對于APGCN-XGBoost模型來說,分子的結構特征向量為APGCN生成的圖嵌入向量,對于APGCN-XGBoost模型來說,分子的結構特征向量為Morgan分子指紋,中間節點不斷進行分裂生長,直到葉子節點結束。

圖3 決策樹生成

XGBoost將增益作為特征選擇和分裂點選擇的準則,增益值越大,說明分裂后能使目標函數減小越多,就越好,用增益值來尋找每次分裂的最優特征和最優分裂點。通過不斷地進行特征分裂來生成一棵樹,每次添加一棵樹,其實是學習一個新函數,去擬合上次預測的殘差。在目標函數中加入了正則項控制模型復雜度,以防止過擬合,并對目標函數進行二階泰勒展開,保留了更多有關目標函數的信息。最終得到XGBoost的最優目標函數如式5所示,

(5)

其中,Gj,Hj分別為損失函數對前t-1輪迭代的一階導數和二階導數和,γ和λ為加權因子,ωj為葉子節點權重,T為葉子的個數。

當訓練完成后得到了多棵樹,對一個樣本進行預測時,每棵樹的節點對相應的離子液體的分子指紋和溫度T、壓力P進行判斷,最終每棵樹都會落到對應的一個葉子節點,模型最終預測的CO2溶解度為每棵樹對應的葉子節點之和。

2 實驗及結果分析

2.1 實驗環境及參數設置

設計APGCN為4層、最終生成的圖嵌入維度為50的向量構建APGCN-XGBoost預測模型進行實驗。選擇半徑為8、陽離子和陰離子指紋長度均為512的Morgan分子指紋[22-23]來構建Morgan-XGBoost預測模型進行對比實驗,RDKit將SMILES轉化Morgan分子指紋,也就是指定長度的二進制比特串。分子結構、T、P對于離子液體的CO2溶解度均為重要的因素,故最終不管是哪種分子表示方式,都需要結合相應的T、P作為預測模型的輸入。

超參數優化使用了優化工具Optuna,對基于樹的超參數搜索進行了優化,它使用被稱為TPE(“Tree-structured Parzen Estimator”)Sampler的方法,這種方法依靠貝葉斯概率來確定哪些超參數選擇是最有希望的并迭代調整搜索[24]。反復調用目標函數,找到使目標函數最小化的值。最終優化后模型的超參數見表2。

表2 模型的主要超參數

2.2 實驗環境及參數設置

使用的數據來源于Lei等早期的工作[25],其中包含在不同溫度和壓力下多種離子液體的CO2溶解度數據。共有118種離子液體9 897組CO2溶解度數(0.000 064 8-0.951 6 mol·mol-1),溫度243.2-453.15 K,壓力0.798 bar-49 990 kPa。陽離子有[BMIM],[EMIM],[MMIM],[HMIM][OMIM],[BMPYR],[N1,4,4,4],[P1,4,4,4]等,陰離子有[BF4][Cl],[DCA],[NO3],[PF6],[SCN],[C(CN)3],[Tf2N],[MeSO4]等。將數據集劃分出外部測試集,包含8種離子液體的1 040組溶解度數據,這8種離子液體未被包含在訓練集中。將訓練集中的實驗數據隨機劃分出20%作為測試集,測試集與外部測試集中的數據均不參與模型的訓練。

2.3 評價指標

(6)

(7)

2.4 對比分析

不同模型的相關參數與訓練集和測試集的預測性能如表3所示。其中Daryayehsalameh B僅僅對[BMIM][BF4]中CO2溶解度預測建模[26],而APGCN-XGBoost模型可以對多種離子液體實現更準確的預測。

表3 模型的預測結果

由預測結果的對比可以看出,APGCN獲取的圖嵌入構建的模型測試集的預測結果與實驗值的AARD最小,R2最大,能夠實現離子液體中CO2溶解度更為準確的預測。APGCN-XGBoost與Morgan-XGBoost,GCN-XGBoost預測模塊均采用XGBoost模型,三個模型外部測試集決定系數R2分別為0.943 3,0.729 9,0.924 1,外部測試集的結果說明APGCN-XGBoost模型對于未知的離子液體也能實現更為準確的預測。

APGCN-XGBoost模型不僅在預測準確度方面表現優秀,同時相比于量子化學、分子軌道、拓撲信息、基團信息等描述符的計算,APGCN獲取的分子圖嵌入無需復雜的計算、專業的軟件或者物理化學先驗知識,特征工程成本低,易于推廣到新的離子液體。綜合表明APGCN實現了更好的分子表征,注意力池化層的加入提高了離子液體結構信息的提取效率,提高了預測模型的準確性、泛化性。

訓練集與測試集的實驗值與APGCN-XGBoost模型預測值進行比較,圖4(a)實驗值與預測值的擬合效果表明模型能夠實現較為準確的預測;圖4(b)可以看出測試集大多數誤差落在[-0.1,0.1]之間,表明模型的預測值與實驗值具有較高的對應性。

圖4 實驗值和APGCN-XGBoost模型預測的CO2溶解度值比較

實驗結果驗證了模型的有效性,證明模型可以很好地提取離子液體中的特征信息用于CO2溶解度預測,不僅在準確度上具有明顯的優勢,且通過圖卷積獲取的圖嵌入維度更低,這能夠降低下游的計算量,為預測CO2在離子液體中的溶解度提供了一種便捷、準確的方法。

2.5 模型學習效果

原子的特征向量在學習過程中不斷變化,可通過計算不同層中不同原子的特征向量之間的Pearson相關系數(取值為-1到1之間,接近-1或1被稱為具有強負相關或正相關),并繪制原子相似性矩陣的熱圖,以觀察學習過程中的變化。以[BBIM][MeSO3]為例,原子的標號為0-17,第0-4層的原子相關性熱圖見圖5,具體標號對應位置如圖5(f)所示。在訓練之前,如圖5(a),不同原子之間相似性混亂,并沒有明顯的區分,而在訓練過程中,具有相似化學特性的原子的特征向量越來越相似,最終整個離子液體按照原子的化學環境分為幾個區域,每個區域內的原子的化學特性相似,如圖5(e)所示,原子0-3在同一支鏈上,具有相似的節點特征且隨著APGCN層數的增加,幾個原子特征向量之間的相關系數增大,最終為1,也就是說原子特征向量之間相似性也隨之提高。

圖5 [BBIM][MeSO3]的熱圖

具有相似的節點特征并且在拓撲上彼此接近的原子在最終訓練后的嵌入空間中具有更高的相似性。說明APGCN通過學習成功提取了原子化學特性相關的信息,并且能夠捕獲到原子的非局部效應。

2.6 模型的可解釋性

模型不僅能夠進行準確的預測,還能對模型的預測機理進行解釋,實現結構-性質之間的關聯,借助注意力池化層計算的注意力系數與SHAP(SHapley Additive exPlanation)[27]實現模型的解釋。首先通過分析溫度、壓力對模型溶解度預測的影響,驗證模型中特征重要性的可信度。SHAP的部分依賴圖(partial dependence plot,pdp)可以顯示特征與預測值之間的關系。如圖6所示,pdp圖顯示了溫度、壓力的不同取值與模型預測的溶解度數值之間的關系?？梢钥闯?離子液體中的CO2溶解度隨著溫度的升高而降低,隨著壓力的增加而增加,結果與溶解度隨著溫度和壓力變化的理論是一致的,說明模型的預測和解釋是可信的。

圖6 溫度和壓力的pdp

用注意力池化層得到的注意力系數來衡量離子液體中結構片段對于溶解度預測的重要性,最終導出對溶解度預測影響較大的結構信息,并將其與已知的結論進行對比。相應的結構信息見圖7,分子圖中的淺色部分為標注的重要結構片段。從結果可以看出:標注結構均為陰離子,S=O基團與含氟陰離子對溶解度影響較大。已知:離子液體CO2溶解度中陰離子起關鍵作用[28];S=O基團能夠提高離子液體CO2溶解度[29];陰陽離子氟化都能使溶解度升高,陰離子氟化作用更為明顯[30]。通過對離子液體CO2溶解度重要結構的可視化,尋找溶解度預測中重要的子結構。將化學常識與模型自動篩選得到的重要特征對比,有利于進一步研究結構與性質的關系,對碳捕集過程中的理想離子液體的篩選與設計具有重要意義。

圖7 對模型預測貢獻較大的結構片段

從模型的結果可以看出,APGCN能夠學習到離子液體分子圖中更全局的特征,基于APGCN的離子液體表征方式足夠“學習”化學結構與性質之間的關系,并且APGCN-XGBoost模型的開發和使用不依賴于化學反應、化學性質和描述符計算的任何高級先驗知識,準確預測的同時還可以指出哪些化學結構與化學特性最相關,可以將這些結構篩選出來,此項工作旨在通過深度學習來解釋結構與性質的關聯,作為相關離子液體的設計和篩選的輔助工具。這種方法還可以很容易地擴展到許多其他領域,用于多種化學分子的性質預測以及篩選、結構優化工作。

3 結束語

為實現離子液體中CO2溶解度的高效預測,通過圖卷積神經網絡從離子液體的分子圖中捕獲有效的特征信息,將APGCN與XGBoost組成高效的預測模型,圖數據經過APGCN轉化為包含結構特征的低維向量,結合其他參數數據作為輸入特征,APGCN-XGBoost模型輸出溶解度的預測值。最終結果表明,對于新的離子液體,APGCN能夠實現比其他分子表征方式更為準確的預測,且對比于傳統的熱力學、量子化學、拓撲信息等表征方式特征工程成本低,更容易推廣到新的離子液體,說明APGCN實現了更好的分子表征。一方面,APGCN-XGBoost模型能夠實現在較寬的溫度、壓力范圍內離子液體CO2溶解度準確的預測,解決了實驗測量耗時、成本高、傳統模型結構復雜、通用性差、計算成本高等問題,為預測CO2在離子液體中的溶解度提供了一種高效的方法。另一方面,APGCN成功學習到了各個原子和結構的特征信息,這些特征信息不僅可以用于預測離子液體的性質,還可以用于探索不同的原子和結構是如何影響最終的性質預測的,有利于進一步研究結構與性質之間的關系。