基于門控圖卷積神經(jīng)網(wǎng)絡(luò)的有機(jī)化學(xué)反應(yīng)預(yù)測

賴自成，張玉萍，馬 燕

（上海師范大學(xué)信息與機(jī)電工程學(xué)院，上海 201418）

0 引言

健康是人類永恒的目標(biāo)，而疾病和病毒往往對人的生命構(gòu)成挑戰(zhàn)，采用人工智能技術(shù)輔助設(shè)計藥物來加速藥物的研發(fā)進(jìn)度成為研究熱點［1］。機(jī)器學(xué)習(xí)在化學(xué)領(lǐng)域已展現(xiàn)出巨大的應(yīng)用潛力，在量子化學(xué)、密度泛函計算、藥物發(fā)現(xiàn)、分子從頭設(shè)計、反應(yīng)預(yù)測與逆合成分析等領(lǐng)域均已有重要應(yīng)用［2］。四川大學(xué)的李秉軻等［3］使用分子描述符和機(jī)器學(xué)習(xí)方法預(yù)測和篩選基質(zhì)金屬蛋白酶的選擇性抑制劑。美國密西根大學(xué)研究團(tuán)隊利用Synthia?人工智能軟件對12種尚在研究階段的抗新冠藥物展開了逆合成研究，嘗試使用廉價易得的原料進(jìn)行藥物的合成［4］。傳統(tǒng)的新藥研發(fā)面臨著投入高、耗時長、成功率低的困境［5］。開發(fā)新藥的平均成本約為28 億美元，新藥研發(fā)的平均時長約為14 a［6］。Stokes 等［7］基于圖卷積神經(jīng)網(wǎng)絡(luò)從1億多個分子中發(fā)現(xiàn)了強(qiáng)大的新型廣譜抗生素halicin。在實驗中，其他抗生素化合物的耐藥性通常在一兩天內(nèi)出現(xiàn)，但經(jīng)過30 天的檢測，也沒有發(fā)現(xiàn)細(xì)菌對halicin 有任何耐藥性。雖然人工智能技術(shù)在藥物發(fā)現(xiàn)方面取得初步進(jìn)展，但由于化合物搜索空間太過龐大、化學(xué)反應(yīng)規(guī)則太過龐雜，所以通過人工智能技術(shù)設(shè)計藥物合成路線仍然是一個巨大的挑戰(zhàn)。目前主流的單步逆合成反應(yīng)模型的預(yù)測準(zhǔn)確率都比較低，例如GLN（Graph Logic Network）模型［8］在USPTO-full 測試集上的Top-1準(zhǔn)確率為39.3%。在難以提高單步逆合成反應(yīng)模型性能的情況下，使用一套可靠的正向有機(jī)化學(xué)反應(yīng)預(yù)測算法用以篩除單步逆合成反應(yīng)模型中推薦的不合理的反應(yīng)變得尤為關(guān)鍵。

1 有機(jī)化學(xué)反應(yīng)預(yù)測模型研究現(xiàn)狀

計算機(jī)輔助化學(xué)反應(yīng)預(yù)測有著豐富的研究歷史。從1980年起，Salatin等［9］便提出了計算機(jī)輔助評價有機(jī)化學(xué)反應(yīng)機(jī)理算 法 CAMEO（Computer Assisted Mechanistic Evaluation of Organic reactions），更早的類似的專家系統(tǒng)算法還有EROS（Elaboration of Reactions for Organic Synthesis）［10］，IGOR（Interactive Generation of Organic Reactions）［11］，SOPHIA（System for Organic Reaction Prediction by Heuristic Approach）［12］，和ROBIA（Reaction Outcomes By Informatics Analysis）［13］。這些算法在設(shè)計之初都希望能夠?qū)θ我庖粋€普通的化學(xué)反應(yīng)進(jìn)行反應(yīng)機(jī)理的預(yù)測分析，但是最后都沒有能夠得到廣泛的應(yīng)用。目前有機(jī)化學(xué)反應(yīng)預(yù)測模型主要有基于規(guī)則的模型、神經(jīng)序列到序列的模型以及圖卷積神經(jīng)網(wǎng)絡(luò)模型。

1.1 基于規(guī)則的模型

基于規(guī)則的模型主要分為模板法和通式規(guī)則兩大類，均是先通過子圖匹配找到相應(yīng)規(guī)則，然后根據(jù)制定的反應(yīng)規(guī)則來預(yù)測反應(yīng)。兩者在使用上大致相同，主要區(qū)別在于通式規(guī)則是由化學(xué)專家總結(jié)得來，而模板規(guī)則是由算法從反應(yīng)中自動提取得到。模板規(guī)則比較呆板，因為只能關(guān)注反應(yīng)中心周圍的環(huán)境情況所以泛化性能一般。通式規(guī)則比較靈活，但是化學(xué)反應(yīng)規(guī)則龐雜，人工無法列舉所有的反應(yīng)規(guī)則。Coley等［14］從一個擁有100 萬化學(xué)反應(yīng)的數(shù)據(jù)庫中提取了1.4×105個反應(yīng)模板，由于模板法涉及圖匹配，這使得測試大量模板的時間變得難以承受。因此模板法被限制在較小的數(shù)據(jù)集以及有限的反應(yīng)類型上。此外新藥研發(fā)極可能涉及新的化學(xué)反應(yīng)，基于規(guī)則的模型難以預(yù)測規(guī)則之外的化學(xué)反應(yīng)，所以基于規(guī)則的模型在新藥研發(fā)的使用上會有所限制。

1.2 神經(jīng)序列到序列

由于三維空間結(jié)構(gòu)的化合物可以依照簡化分子線性輸入規(guī)范（Simplified Molecular Input Line Entry Specification）轉(zhuǎn)換為一維的SMILES表達(dá)式，因此也可用神經(jīng)序列到序列的方法對化合物進(jìn)行分析。Goh 等［15］利用長短時記憶（Long Short-Term Memory，LSTM）網(wǎng)絡(luò)將SMILES 表達(dá)式編碼為分子向量，并借以預(yù)測化合物的細(xì)胞毒性和溶解度。Schwaller 等［16］利用LSTM 和Luong 注意力機(jī)制進(jìn)行化學(xué)反應(yīng)預(yù)測，利用LSTM 對反應(yīng)物的SMILES 表達(dá)式進(jìn)行編碼，再通過LSTM 解碼得到生成物的SMILES 表達(dá)式。序列到序列的方法比較靈活，能夠預(yù)測數(shù)據(jù)庫之外的化學(xué)反應(yīng)，但缺點是不易遵守原子守恒，有時預(yù)測的生成物會缺失或多出幾個原子。

1.3 圖卷積網(wǎng)絡(luò)

傳統(tǒng)圖像領(lǐng)域的卷積是在歐氏空間進(jìn)行的操作，圖像和卷積核都是規(guī)則的矩陣形狀。卷積核在圖像上滑動，對應(yīng)位置相乘求和即可。一個3×3大小的卷積核進(jìn)行的運算，可以看成為一個像素點與周圍9 個相鄰像素點之間的信息整合。而在非歐氏空間中每個節(jié)點的相鄰節(jié)點數(shù)并不是固定的，所以需要與之相應(yīng)的圖卷積網(wǎng)絡(luò)來進(jìn)行信息整合。圖卷積網(wǎng)絡(luò)是一種在非歐氏空間中利用局部節(jié)點信息獲取全局節(jié)點信息的方法［17］?；衔镏忻恳粋€原子可以視為一個節(jié)點，每個節(jié)點將周圍節(jié)點信息進(jìn)行聚集再整合，如此循環(huán)往復(fù)使得局部節(jié)點信息能在整張圖上進(jìn)行流動。Jin 等［18］基于WLDN（Weisfeiler-Lehman Difference Network）模型預(yù)測有機(jī)化學(xué)反應(yīng)，該模型能夠更好地處理有機(jī)化學(xué)分子式這類節(jié)點和邊均有不同形式的異構(gòu)圖，預(yù)測準(zhǔn)確率達(dá)到了85.6%，該模型采用的圖卷積網(wǎng)絡(luò)如圖1所示，圖中NN 表示帶sigmoid 激活函數(shù)的全連接神經(jīng)網(wǎng)絡(luò)（Neural Network with sigmoid activation function）。圖卷積網(wǎng)絡(luò)目前主要通過不斷迭代遞推公式來獲取更加全局的信息，這將導(dǎo)致所有節(jié)點的信息趨于平均。根據(jù)熱力學(xué)第二定律，一個孤立系統(tǒng)的熵值不會減小，最終系統(tǒng)將趨于無序的狀態(tài)。圖卷積中節(jié)點信息聚集加和再取平均的過程可以看成熱力學(xué)模型中熱能從高溫物體流向低溫物體，最終溫度達(dá)到一致的過程，這同樣印證了傳統(tǒng)圖卷積網(wǎng)絡(luò)中節(jié)點信息易趨于平均。

圖1 圖卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 Structure of graph convolutional neural network

2 門控圖卷積神經(jīng)網(wǎng)絡(luò)模型

為了解決現(xiàn)有模型中出現(xiàn)的模板規(guī)則無法窮舉、預(yù)測結(jié)果原子不守恒、無法有效獲取全局信息的問題，提出了門控圖卷積神經(jīng)網(wǎng)絡(luò)模型，取得了一定效果。

解決方案有兩大基本步驟：其一輸入化學(xué)反應(yīng)物通過門控圖卷積神經(jīng)網(wǎng)絡(luò)（Gated Graph Convolutional Neuralnetwork，GGCN）以及注意力機(jī)制預(yù)測候選反應(yīng)中心所在位置，即哪些原子對間可能發(fā)生化學(xué)鍵的斷裂和組合；其二，根據(jù)化學(xué)約束條件和候選反應(yīng)中心枚舉生成候選產(chǎn)物，通過門控圖卷積差分網(wǎng)絡(luò)對候選產(chǎn)物進(jìn)行評估排序，得分最高的候選產(chǎn)物即為最終反應(yīng)產(chǎn)物。

2.1 候選產(chǎn)物的生成

將給定化學(xué)反應(yīng)視為一對分子圖(Gr，Gp)，亦可定義為G=(V，E)。其中反應(yīng)物定義為Gr，反應(yīng)產(chǎn)物定義為Gp，原子集合定義為V={a1，a2，…，an}?；瘜W(xué)鍵集合定義為E={b1，b2，…，bm}，化學(xué)鍵包括單鍵、雙鍵、芳香鍵等類型。反應(yīng)中心即化學(xué)鍵發(fā)生改變的位置，是從反應(yīng)物到產(chǎn)物之間進(jìn)行圖修改所需要的最小集合。每一個Gr中的原子對(au，av)與一個二值反應(yīng)標(biāo)簽yuv相關(guān)聯(lián)，若原子對間的關(guān)系發(fā)生改變則yuv為真。原子au和原子av的節(jié)點特征定義為fu以及fv，節(jié)點之間邊buv的特征定義為fuv。

在步驟一中首先通過門控圖卷積神經(jīng)網(wǎng)絡(luò)（GGCN）獲取節(jié)點的隱藏特征，然后通過注意力機(jī)制預(yù)測原子對間的反應(yīng)活性值，反應(yīng)中心即在反應(yīng)活性值較高的原子對處。門控圖卷積神經(jīng)網(wǎng)絡(luò)算法框圖見圖2，門控圖卷積網(wǎng)絡(luò)偽代碼如下：

圖2 GGCN算法框圖Fig.2 Block diagram of GGCN algorithm

算法1 GGCN算法。

輸入 原子鄰接表，化學(xué)鍵鄰接表，原子輸入特征，化學(xué)鍵輸入特征；

輸出 原子隱藏特征。

GGCN 通過以下公式來更新原子的特征。式中Wr、Wz、Wh是跨層共享的矩陣變量，1≤l≤3，，σ(?) 為sigmoid激活函數(shù)，N(v)為原子v的所有鄰居原子組成的集合。

定義αuv為原子v分配到原子u上的注意力值，αuv越高，原子v與原子u越可能發(fā)生關(guān)聯(lián)。cv、cu為通過門控圖卷積神經(jīng)網(wǎng)絡(luò)獲取得到的上下文原子隱藏特征，cv=。原子u和原子v之間的反應(yīng)活性值suv可由下列公式得到，其中U1、U2、Pa、Pb、Ma、Mb、Mc為矩陣變量，buv是一個額外的特征向量用來編碼關(guān)于原子對的輔助備用信息，例如該原子對中的原子是否屬于同一分子或者它們之間的化學(xué)鍵類型。

反應(yīng)活性值suv可以預(yù)測原子和原子之間化學(xué)鍵斷裂與連接的可能性。選取反應(yīng)活性值較高的原子對，枚舉出這些原子對間可能的化學(xué)鍵組合方式，并根據(jù)化學(xué)約束條件用以生成候選產(chǎn)物。在步驟2）中通過門控圖卷積差分網(wǎng)絡(luò)從候選產(chǎn)物中挑選出真正的反應(yīng)產(chǎn)物。

2.2 候選產(chǎn)物的篩選

步驟二中首先將候選產(chǎn)物pi輸入到GGCN 中，得到候選產(chǎn)物pi中原子v的隱藏特征。由于原子守恒，所以反應(yīng)物和反應(yīng)產(chǎn)物中的原子是一一映射的，故通過式（9）定義差分向量來關(guān)注原子隱藏特征的變化。

差分向量僅當(dāng)靠近反應(yīng)中心時才會偏離零值，因此它會關(guān)注于處理反應(yīng)中心和它鄰近的信息。將差分向量輸入到GGCN 中，經(jīng)過L次循環(huán)，得到差分向量的隱藏特征，對其進(jìn)行加和池化得到反應(yīng)分?jǐn)?shù)，取反應(yīng)分?jǐn)?shù)最高的候選產(chǎn)物為最后結(jié)果。反應(yīng)分?jǐn)?shù)由式（10）計算得到，其中U3及M為矩陣變量。

3 實驗和分析

本實驗采用經(jīng)麻省理工大學(xué)整理的美國專利局?jǐn)?shù)據(jù)集USPTO［19］，有機(jī)化學(xué)反應(yīng)方程式在訓(xùn)練集多達(dá)409 035 條、測試集有為4 萬條。模型使用主動采樣方法進(jìn)行訓(xùn)練，即學(xué)習(xí)完一輪完整的訓(xùn)練集數(shù)據(jù)后，采集loss 值較高的前40%數(shù)據(jù)進(jìn)行訓(xùn)練，再在此基礎(chǔ)上選擇loss 值高的數(shù)據(jù)再進(jìn)行訓(xùn)練，如此迭代3 次后再學(xué)習(xí)一輪完整的訓(xùn)練集數(shù)據(jù)進(jìn)行下一輪主動采樣訓(xùn)練。假設(shè)訓(xùn)練集有40萬條，則先學(xué)習(xí)完整的40萬條訓(xùn)練樣本，然后取loss值較高的前16萬條訓(xùn)練樣本再進(jìn)行訓(xùn)練，記錄訓(xùn)練過程中每條訓(xùn)練樣本對應(yīng)的loss 值，再取16 萬條中l(wèi)oss 值較高的6.4 萬條訓(xùn)練樣本再訓(xùn)練。完成上述訓(xùn)練后繼續(xù)學(xué)習(xí)完整的40萬條訓(xùn)練樣本，開始下一輪主動采樣訓(xùn)練。

模型基于化學(xué)工具包rdkit 進(jìn)行輸入特征的提取。原子輸入特征有97 個維度，其中63 個維度用于表示原子序數(shù)，其他維度用于表示元素周期、原子連接數(shù)、原子顯示化合價、原子隱式化合價、原子芳香性等。化學(xué)鍵輸入特征共有6 個維度，分別是單鍵、雙鍵、三鍵、芳香性、共軛、成環(huán)。門控圖卷積神經(jīng)網(wǎng)絡(luò)中所有隱藏層均為300維，采用Adam 優(yōu)化器進(jìn)行訓(xùn)練，初始學(xué)習(xí)率為0.001，每一萬次迭代學(xué)習(xí)率衰減為原來的90%。模型的損失函數(shù)為交叉熵?fù)p失函數(shù)加上L2 范數(shù)正則化損失。訓(xùn)練過程中采用了梯度截斷技巧，梯度裁剪率設(shè)置為5可以有效避免發(fā)生梯度爆炸。

WLDN模型、門控圖卷積神經(jīng)網(wǎng)絡(luò)（GGCN）模型和主動采樣訓(xùn)練下的門控圖卷積網(wǎng)絡(luò)（Active Sampling-training Gated Graph Convolutional Neural-network，ASGGCN）模型在USPTO測試集中的覆蓋率實驗結(jié)果如圖3。圖3 中覆蓋率指所有反應(yīng)中心均被預(yù)測到的準(zhǔn)確率，定義原子對計數(shù)符號為K，定義反應(yīng)活性值最高的前K個原子對為集合TOPK。若所有反應(yīng)中心原子對（反應(yīng)中心即化學(xué)鍵發(fā)生改變的位置）組成的集合包含于集合TOPK，且預(yù)測的化學(xué)鍵變化符合實際產(chǎn)物，則稱集合TOPK覆蓋了該有機(jī)化學(xué)反應(yīng)的反應(yīng)中心。由圖可見每個模型隨著K的增大覆蓋率均有提高，收斂后ASGGCN 模型在TOP12、TOP16和TOP20集合中均優(yōu)于其他兩個模型。

圖3 模型的反應(yīng)中心覆蓋率Fig.3 Reaction center coverage rates of models

在ASGGCN 模型的收斂過程存在一定的波動，如第3×104次迭代時覆蓋率出現(xiàn)較大跌幅，隨著迭代次數(shù)的增加波動逐漸趨緩，在1.3×105次迭代后覆蓋率穩(wěn)居最高。原因是采用了主動采樣的方式進(jìn)行訓(xùn)練，訓(xùn)練初期對較差樣本進(jìn)行了重點學(xué)習(xí)，由于訓(xùn)練參數(shù)尚未穩(wěn)定，所以對普通樣本的分析能力有部分的遺忘，而較差樣本所占比例較低造成了初期整體覆蓋率的下跌。訓(xùn)練后期，訓(xùn)練參數(shù)逐漸穩(wěn)定，在對較差樣本進(jìn)行學(xué)習(xí)時普通樣本的分析能力仍能較好地保留。在1.6×105次迭代時覆蓋率未出現(xiàn)明顯下跌，表明ASGGCN 模型能夠兼顧較差樣本和普通樣本的分析能力。

ASGGCN 模型只考慮前16 個候選反應(yīng)中心進(jìn)行預(yù)測，此時ASGGCN 模型與WLDN 模型在USPTO 測試集上的實驗比較結(jié)果如圖4。圖4中Top-1準(zhǔn)確率為預(yù)測一次便成功預(yù)測的概率，將模型預(yù)測的產(chǎn)物與實際產(chǎn)物的SMILES字符串表達(dá)式進(jìn)行一一比較，任意一位匹配不成功，則該次預(yù)測失敗。橫坐標(biāo)為反應(yīng)方程式的模板匹配數(shù)，若反應(yīng)越罕見則模板匹配數(shù)越低，所有模型的準(zhǔn)確率隨模板匹配數(shù)的降低而下降。ASGGCN 模型在任意模板匹配數(shù)時均優(yōu)于WLDN 模型。在模板匹配數(shù)較低，為5～49 及小于5 時，ASGGCN 的優(yōu)勢更加明顯，二者Top-1的領(lǐng)先幅度分別為2.8%及3.0%。

圖4 預(yù)測準(zhǔn)確率與反應(yīng)模板數(shù)的關(guān)聯(lián)性Fig.4 Correlation between prediction accuracy and the number of reaction templates

各模型進(jìn)行有1、2、3 和5 次機(jī)會預(yù)測時，在USPTO 測試集中反應(yīng)產(chǎn)物預(yù)測成功的概率如表1 所示。實驗結(jié)果表明ASGGCN 模型參數(shù)規(guī)模最小，僅有2.4×106，但在Top-1、Top-2、Top-5指標(biāo)上明顯優(yōu)于其他模型。

表1 反應(yīng)產(chǎn)物預(yù)測準(zhǔn)確率Tab.1 Prediction accuracy of reaction product

4 模型可解釋性

以環(huán)氧丙烷與吡咯并喹啉的反應(yīng)為例，輸入反應(yīng)物的SMILES表達(dá)式，基于rdkit的預(yù)處理程序?qū)⒆詣犹崛≡虞斎胩卣骱突瘜W(xué)鍵輸入特征，模型將輸出反應(yīng)產(chǎn)物的SMILES表達(dá)式，并基于式（6）計算得到全局注意力分?jǐn)?shù)矩陣。圖5 為WLDN模型與ASGGCN模型的對比分析。

圖5 模型可解釋性分析Fig.5 Model interpretability analysis

圖5 中，圖（a）～（f）中的橫縱坐標(biāo)均表示原子序號，為圖6中原子所對應(yīng)的序號，圖（a）、圖（b）、圖（d）、圖（e）右側(cè)條形坐標(biāo)為原子相似度值，圖（c）、圖（f）右側(cè)條形坐標(biāo)為原子注意力值。圖5 中圖（a）和圖（d）為輸入層原子的皮爾遜相關(guān)系數(shù)矩陣，可以看見子圖（a）和（d）結(jié)果一致，且由于1，10，17 號原子為非C 原子，與其他原子的相關(guān)系數(shù)在0.4 以下。圖5 中圖（b）和（e）為圖2 圖卷積操作之后的隱藏層中原子的皮爾遜相關(guān)系數(shù)矩陣，可以看出矩陣根據(jù)反應(yīng)物序號大體分為兩塊，其中14，15，16，17 號原子屬于環(huán)氧丙烷，ASGGCN 的劃分結(jié)果更為明顯。圖5（c）為WLDN 的全局注意力矩陣，圖5（f）為ASGGCN的全局注意力矩陣。

圖6為WLDN 模型下該反應(yīng)的反應(yīng)中心14號C 原子對其他原子的注意力情況，原子顏色越深注意力值越大，其中反應(yīng)中心原子被標(biāo)注為綠色。注意力值可以反映化學(xué)鍵的連接情況，反應(yīng)產(chǎn)物中14號C 原子與1號N 原子相連，故圖5（f）中14號C 原子對1 號N 原子分配了較多的注意力。而WLDN 模型對反應(yīng)中心的判斷能力稍弱，故在圖5（c）中14 號C 原子對1號N原子和10號N原子都分配了較多的注意力。

圖6 環(huán)氧丙烷與吡咯并喹啉反應(yīng)Fig.6 Propylene oxide reaction with pyrroloquinoline

5 結(jié)語

本文提出了一種基于主動采樣訓(xùn)練及門控圖卷積網(wǎng)絡(luò)的有機(jī)化學(xué)反應(yīng)預(yù)測模型。給定化學(xué)反應(yīng)物，模型可預(yù)測相應(yīng)產(chǎn)物。ASGGCN 擁有三個權(quán)重參數(shù)矩陣并通過門控對信息加以融合，所以ASGGCN 能夠獲取更加豐富的原子隱藏特征信息；同時ASGGCN通過主動采樣的方式進(jìn)行訓(xùn)練，能夠兼顧較差樣本和普通樣本的分析能力。最終ASGGCN 的反應(yīng)產(chǎn)物Top-1 預(yù)測準(zhǔn)確率達(dá)到了87.2%，相較于WLDN 模型提高了1.6個百分點。