賴自成，張玉萍，馬 燕

（上海師范大學信息與機電工程學院，上海 201418）

0 引言

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

1 有機化學反應預測模型研究現(xiàn)狀

計算機輔助化學反應預測有著豐富的研究歷史。從1980年起，Salatin等［9］便提出了計算機輔助評價有機化學反應機理算 法 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ū)θ我庖粋€普通的化學反應進行反應機理的預測分析，但是最后都沒有能夠得到廣泛的應用。目前有機化學反應預測模型主要有基于規(guī)則的模型、神經(jīng)序列到序列的模型以及圖卷積神經(jīng)網(wǎng)絡模型。

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

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

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

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

1.3 圖卷積網(wǎng)絡

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

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

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

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

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

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

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

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

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

算法1 GGCN算法。

輸入 原子鄰接表，化學鍵鄰接表，原子輸入特征，化學鍵輸入特征；

輸出 原子隱藏特征。

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)絡獲取得到的上下文原子隱藏特征，cv=。原子u和原子v之間的反應活性值suv可由下列公式得到，其中U1、U2、Pa、Pb、Ma、Mb、Mc為矩陣變量，buv是一個額外的特征向量用來編碼關(guān)于原子對的輔助備用信息，例如該原子對中的原子是否屬于同一分子或者它們之間的化學鍵類型。

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

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

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

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

3 實驗和分析

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

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

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

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

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

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

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

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

表1 反應產(chǎn)物預測準確率Tab.1 Prediction accuracy of reaction product

4 模型可解釋性

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

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

圖5 中，圖（a）～（f）中的橫縱坐標均表示原子序號，為圖6中原子所對應的序號，圖（a）、圖（b）、圖（d）、圖（e）右側(cè)條形坐標為原子相似度值，圖（c）、圖（f）右側(cè)條形坐標為原子注意力值。圖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ù)反應物序號大體分為兩塊，其中14，15，16，17 號原子屬于環(huán)氧丙烷，ASGGCN 的劃分結(jié)果更為明顯。圖5（c）為WLDN 的全局注意力矩陣，圖5（f）為ASGGCN的全局注意力矩陣。

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

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

5 結(jié)語

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