融合多模態數據的藥物合成反應的虛擬篩選

孫曉飛，朱靜遠，陳 斌，游恒志*

（1.中國科學院 成都計算機應用研究所，成都 610041；2.中國科學院大學 計算機科學與技術學院，北京 100049；3.哈爾濱工業大學（深圳）理學院，廣東 深圳 518055；4.哈爾濱工業大學（深圳）人工智能研究院，廣東 深圳 518055）

0 引言

藥物合成中的有機化學反應的發現依賴于實踐經驗和化學機理支配的“化學直覺”，實驗人員試圖定性地識別有機化學反應中的模式，以確定反應產物和反應效率。然而，這種方法受到了很多因素的限制，包括反應的復雜性、活性懸崖、對機理理解的缺乏以及人工處理大數據的艱難。基于計算機的虛擬篩選［1-5］已經成為吸引化學家的重要方案，主要因為它不需要機理的理解，化合物結構可以用分子性質的數值表示來表征，從而量化數以千計的候選分子的化學性質。在實驗和文獻數據的基礎上，虛擬篩選可以通過計算機模型來對藥物合成反應的結果和催化劑的選擇性進行量化。

機器學習在化學領域已成功應用于藥物虛擬篩選［3-5］、分子生成［6］、有機反應預測［7-8］、催化劑篩選［9-10］、材料發現［11］、計算機輔助合成設計［4，12］和反應條件優化。線性回歸是傳統的反應預測和分析工具［13-14］，它假設反應物的物理特征和反應性之間存在線性關系，可以根據反應的機制假設人工對輸入變量進行選擇，所以非常符合數據科學家的思維和統計方式。Hammett［15］在線性自由能關系的推斷中使用線性回歸擬合化合物描述符和輸出是一個代表性工作。長期以來，由于分子特征的多維性和反應空間的復雜性，很難生成足夠完整和一致的數據，從而限制了機器學習的發展［11］。如今，高通量實驗（High Throughput Experimentation，HTE）已經成為逐步掃除這一障礙的有效手段［8，16-20］。Ahneman 等［21］使用了支持向量機（Support Vector Machine，SVM）和隨機森林（Random Forest，RF）等方法在4 000 多個高通量實驗數據中預測了Buchwald-Hartwig 偶聯反應的產率。此外，Zahrt等［22］通過RF 對1 000 多個反應中的手性磷酸（Chiral Phosphoric Acid，CPA）催化劑的對映選擇性進行了預測。

使用計算機對化學反應進行虛擬篩選的流程如圖1 所示，首先從已有的化學反應數據庫和文獻中提取分子的簡化分子線性輸入（Simplified Molecular-Input Line-Entry System，SMILES）或分子指紋，或者使用Gaussian 等密度泛函（Density Functional Theory，DFT）工具對這些分子進行結構優化并計算出與反應有關的性質；然后，用這些物理和化學性質構建出分子描述符，再選用合適的機器學習方法進行建模；最后對數據集中待分析的反應進行篩選。這種方法對于數據科學家來說直觀有效，不需要關注反應機理的理解，已經成為化學合成預測的標準流程。該流程的成功與否，取決于兩個關鍵因素：1）所選的DFT 特征或者分子指紋，以及用它們構建的描述符是否準確；2）機器學習方法是否有效。藥物合成相關的有機反應預測經過幾十年發展，在這兩個方面仍受制約。下面對這兩個問題進行詳細描述：

圖1 對藥物合成反應進行虛擬篩選的流程Fig.1 Flow of virtual screening of drug synthesis reactions

1）對于基于量子力學的DFT 特征來說，針對不同反應進行特征的選擇一直是藥物合成相關的預測需要面對的難題，特別是對反應產率和選擇性預測的特征選擇往往有很大差別。若能降低特征選擇的難度，將為藥物合成相關的反應預測帶來促進作用。對于使用SMILES 和分子指紋的序列特征來說，對于三維（3D）結構信息表達不足是一直存在的困難。這是由SMILES 作為一種簡化的分子結構線性表示以及分子指紋的算法本質決定的。

2）在藥物合成相關的化學反應預測中，傳統的機器學習方法如SVM 和RF 等，甚至是線性回歸方法一直占據主流。由于“維度災難”的存在，隨著特征維度的上升，所需的化學反應數據急劇上升，大幅超出了人工實驗的工作量。而高通量實驗的出現，使這一問題逐步得到緩解。如今，如何將高通量化學反應數據應用于深度學習，已經擺在數據科學家面前。然而，由于對化學知識的缺乏和反應數據相對稀少（相對于傳統深度學習應用領域，如圖像、視頻、音頻和文本），深度學習方法在這一領域的研究仍然罕見。雖然已經有一些工作對文獻中的反應數據提取SMILES 后使用深度學習方法進行預測［23］，但如何使用深度學習方法對日漸累積的DFT反應數據進行虛擬篩選仍亟待研究［4，24］。

針對上述問題，本文主要工作如下：1）提出基于加權平均占位和分子DFT 性質的描述符，并用于4-甲基苯胺與芳基鹵化物的C-N 偶聯反應的預測。這種描述符以構象能量計算權重，得出多個構象的平均值來近似反應發生時的構象概率分布。2）提出一種針對DFT 計算性質的圖卷積網絡，用量子力學（Quantum Mechanics，QM）計算的原子性質作為圖節點的輸入特征，并融合分子指紋特征進行預測。使用這種網絡在CPA 催化的N，S-縮醛反應上進行驗證。據了解，現在尚無這種使用DFT 計算的性質構建的圖卷積網絡，這是一種將QM 特征引入深度學習方向的有益嘗試。3）針對現有工作對多模態數據的表示和利用不足，采用基于卷積神經網絡（Convolutional Neural Network，CNN）或圖卷積神經網絡（Graph Convolutional Network，GCN）［25］等深度學習手段，將所提出的3D描述符與其他來源的分子描述符融合起來，應用于反應產率和對映選擇性的預測中。

在C-N 偶聯反應產率預測中，采用格點精度為1 埃（0.1 nm）的加權平均占位描述符，相較于文獻［21］和文獻［18］提高了3和2個百分點，對僅使用DFT計算的描述符有明顯優勢。在N，S-縮醛反應的對映選擇性預測中，在隨機劃分的數據集上誤差從文獻［22］的0.152降到0.147，驗證了方法的有效性。

1 相關工作

1.1 三維分子描述符

比較分子場分析（Comparative Molecular Field Analysis，CoMFA）作為一種使用了三維空間信息的定量構效關系（Quantitative Structure-Activity Relationship，QSAR）方法，在不對稱催化相關問題上已有近20年的歷史。該方法用分子力學力場近似范德華相互作用，用庫侖勢確定靜電相互作用。這種分子描述符也考慮了分子的三維結構，被認為具有探索化學空間的潛力。該方法被首次用于分析環戊二烯和3-乙烯基惡唑烷-2-酮的Diels-Alder反應中含有膦惡唑啉或雙惡唑啉配體的催化劑［26］的對映選擇性，以及在PhCHO 中添加Et2Zn［27］時氨基醇催化劑的對映選擇性。后來的工作采用了類似的基于CoMFA 的方法，結合了半經驗［28］和量子力學相互作用能，并用于不同的不對稱催化反應，如氧烯丙基陽離子［29］環加成反應和手性sparteine 替代物的不對稱鋰化反應［30］。在最近的杰出工作中，Zahrt等［22］介紹了一種與CoMFA完全不同的新方法，該方法關注催化劑分子的眾多構象，并將構象的平均空間占有率（Average Steric Occupancy，ASO）作為描述符。

1.2 藥物合成反應篩選中的多模態數據融合

隨著高通量實驗和機器學習在藥物合成反應預測中的應用和發展，不同來源、類型和分布的多模態反應數據被產生出來，每種數據都包含特定的信息，如SMILES、分子指紋、量子力學性質和各種人工設計的2D 和3D 描述符。藥物合成反應篩選中的多模態數據具有維度差異大、信息容量差異大和信息類別多樣的特點。由于不同模態數據的維度和物理意義不同，對它們的利用非常困難?，F有的工作都使用單一模態的反應數據，如Ahneman等［21］使用量子力學計算的分子、原子和振動描述符，采用CoMFA 的一系列研究使用分子力學力場近似的相互作用力為基礎構建的描述符。

多模態融合負責聯合多個模態的數據進行分類或者回歸，在深度學習的很多領域應用廣泛。它還有其他常見的別名，如多傳感器融合和多源信息融合。單模態的機器學習將信息表示為計算機可以處理的數值向量并進一步抽象為更高層的特征向量，而多模態機器學習可以通過利用多模態數據之間的互補性，學習到更好的特征表示，如圖2 所示。常見的機器學習模型都可以用于多模態融合［31］。

圖2 多模態數據融合Fig.2 Multimodal data fusion

1.3 圖卷積神經網絡

可以使用圖卷積神經網絡來建模分子相關的問題，將一個化合物構建為一個圖形，它的節點代表原子，連接它們的邊代表鍵。在圖卷積（圖3（a））中，對于一個節點，將特征和鄰居饋送到兩個密集層中，然后添加密集層的輸出作為節點的新特征。對具有相同度的節點計算新特征時共享權重。在一個化合物中，如果一個原子a總共有n個鄰居，那么它經過圖卷積后的新特征可以表示為：

其中：Wa為節點a的權重；Wr是相鄰節點的權重；b是偏差，σ是激活函數。從節點a出發的箭頭表示原子a及其相鄰原子的密集層，權重分別為Wa和Wr。

與CNN 中的池化層類似，圖池化層（圖3（b））也被用于給化合物分子編碼。圖池化是返回原子及其鄰居中最大或平均特征的操作，可以在不增加其他參數的情況下增加感受野：

經過圖卷積和圖池化，每個原子都有一個特征向量；但為了進行最終預測，需要為整個分子圖提供一個固定大小的特征向量。圖聚集層（圖3（c））將化合物分子中所有原子的特征向量相加，以獲得分子的特征向量：

將這三種操作按圖3 中黑色箭頭所示的順序進行組合，就可以構成不對稱反應催化劑篩選網絡中的圖卷積網絡模塊。

圖3 圖卷積神經網絡Fig.3 Graph convolution neural network

2 本文方法

本文提出了一種基于加權平均占位的3D描述符，它在三維空間中整合分子級別性質和原子級別性質，為機器學習模型提供良好的輸入特征，并將這種方法用于鈀催化的Buchwald-Hartwig偶聯反應的篩選；還提出一種基于圖卷積的多模態模型，通過將量子力學性質和分子指紋融合起來來篩選N，S-縮醛反應中的手性磷酸催化劑。下面介紹分子描述符合網絡的結構的一些實現細節。

2.1 基于Boltzmann分布的加權平均占位描述符

在特征選擇和模型建立過程中，需要考慮描述符對有機反應的適應性，即要尋找一組不受限于特定機理假設的描述符，以區分和描述反應之間的差異?；谶@種考慮，在Buchwald-Hartwig 偶聯反應的產率預測中引入一種3D 描述符（圖4（a）），它首先根據結構優化后的分子3D坐標構建出一個三維網格，將分子置于網格中心并采用自開發的分子對齊算法進行對齊，然后統計每個網格點上原子的量子力學性質?？紤]到單一的構象不能很好地表示反應發生時構象的復雜情況，使用Boltzmann 分布計算不同構象所占比例，生成加權平均的3D構象描述符。

其中：p是構象所占比例，k、i是構象的序號，E是構象的能量，R是理想氣體常數，T是開爾文溫度。di是單個構象的描述符，d是經過加權計算后的分子描述符。

3D 坐標不僅表示原子的空間信息，而且對描述符中其他物理性質起到明確的位置編碼作用。其中的分子和原子性質是由Gaussian 計算得到，如最高占據分子軌道（Highest Occupied Molecular Orbital，HOMO）、最低占據分子軌道（Lowest Occupied Molecular Orbital，LUMO）、偶極矩、容積、密立根電荷等（圖4（a））。這些性質的都是易于計算的且經過初步篩選以確保它們的有效性和通用性。

圖4 基于加權平均占位的3D描述符及CNN特征融合模型Fig.4 3D descriptor based on weighted average occupancy and CNN feature fusion model

2.2 反應產率篩選網絡架構

融合多種來源的分子特征數據，對藥物合成相關的有機反應進行篩選，是本次工作的重要目標。基于加權平均占位的3D 描述符可以將原子級別的量子力學特征與立體空間特征融合起來。分子整體的量子力學性質如HOMO、LOMO 等也是與反應機理相關的重要特征，在3D描述符中并不能很好表達這種特征，在此將多個分子級別特征表示為一維向量描述符。3D 描述符數據的維度巨大，采用卷積操作構建網絡以減少過擬合，而對于分子級別特征直接使用密集層構建網絡，然后將兩種特征進行融合（圖4（b））。損失函數的形式如下：

其中：x是樣本，y是化學實驗得到的觀察值。

2.3 不對稱反應中催化劑篩選的網絡架構

使用DFT 計算的量子力學性質具有化學意義且精度高，但由于DFT 計算及其數據的復雜性，至今未見到有工作將量子力學性質用圖卷積的方式來表示和建模。這里從Gaussian的優化結果中提取每個原子的3D 坐標，使用RDKit 提取分子的鄰接矩陣以構建圖結構，并將DFT 計算的量子力學性質作為每個圖節點的屬性向量。這樣構建起來的圖卷積網絡，接受的原子性質和提取的特征具有化學意義，符合化學家的化學直覺。接受SMILES 或分子指紋作為輸入的圖卷積網絡已經被應用于藥物虛擬篩選中［25，32］。構建一個網絡將這些不同種類的特征融合起來，可以發揮優勢互補的作用（圖5）。將這個模型應用于手性磷酸催化的N，S-縮醛反應。

圖5 基于圖卷積的量子力學特征與分子指紋融合模型Fig.5 Quantum mechanical feature and molecular fingerprint fusion model based on graph convolution

3 實驗與結果分析

3.1 有機反應數據和評價指標

本文在兩個藥物相關的有機合成反應上進行了實驗，所采用的數據劃分方法與文獻［21］和文獻［22］等保持一致，以便于比較。第一個數據集是Buchwald-Hartwig 偶聯反應，該反應的底物含有雜原子-雜原子鍵的五元雜環（例如異惡唑），這些雜環化合物具有藥物樣特征，但被成功篩選為候選藥物的數量仍然不足［33］。因此，使用人工智能（Artificial Intelligence，AI）預測異惡唑存在下Buchwald-Hartwig 反應的性能是很有必要的。該反應數據集包含3 960 個高通量實驗的反應物、催化劑、添加劑、溶劑、產物和產率。將反應按添加劑分成4 個測試集，從而可以測試添加劑對反應帶來的影響（表1）。

表1 Buchwald-Hartwig反應的數據集劃分Tab.1 Splitting of Buchwald-Hartwig reaction dataset

第二個虛擬篩選數據集是手性磷酸催化的N，S-縮醛反應。該數據集中反應數量稀少，且該數據集的任務之一是探究底物和催化劑在不對稱催化反應預測中的重要性。因此在數據劃分時，首先將催化劑和底物按訓練用途和測試用途進行劃分，然后將它們組合以形成1 個訓練集和3 個測試集（表2）。實驗結果證明這種劃分可以反映出不同模型對催化劑和底物在反應預測中的表達能力的高低。

表2 N，S-縮醛反應的數據集劃分Tab.2 Splitting of N，S-acetal formation dataset

評價指標采用與已報道工作相同的R2（R-squared）和RMSE（Root Mean Square Error），其中R2指標計算方式為：

其中：分子部分表示觀察值y與預測值的平方差之和；分母部分表示觀察值y與均值的平方差之和。

3.2 實驗細節

在對藥物合成相關反應的數據進行處理時，出于數據一致性的考慮，對于分子統一使用了Gaussian 計算出所需要的性質。為了減少人工數據分析，開發了自動化軟件將分子提交給Gaussian 進行計算，并提取和解析結果。在計算催化劑的性質時，使用ChemDraw、Gaussian 創建初始的分子結構，或者從Cambridge Structural Database［34］導入。經過構象搜索從催化劑的多個構象中獲得能量最低的構象以及不同能量的構象用于生成加權平均構象描述符。計算中未考慮溶劑影響，停止優化的條件為Gaussian 的默認收斂標準。使用自研的程序對催化劑進行對齊，以提供更靈活和自動化的對齊方式。

在AI 模型設計和訓練部分，所有的程序都是基于Pytorch 和Scikit-learn 實現，并將模型預測結果與已報道文獻結果進行對比。鑒于數據規模較小，模型訓練中使用了1 個NVIDIA GeForce GTX 1070，顯存為8 GB。在數據預處理中，由于分子中原子數目不相同，故按列（特征）對數據進行統一長度和縮放處理，這有利于減少無效的數據、縮短描述符長度從而減少網絡參數和提升效果。

3.3 主要結果

首先采用了基于加權平均占位的3D 描述符和分子級別的量子力學特征融合的模型對異惡唑存在下的Buchwald Hartwig 偶聯反應進行篩選。在已有工作中，都使用傳統機器學習算法來評判一個反應描述符的好壞，因此本文也做了類似的比較，如圖6 所示。

圖6 使用基于加權平均占位的3D描述符的6種方法在C-N偶聯反應上的性能對比Fig.6 Performance comparison of six methods using 3D descriptors based on weighted average occupancy in C-N coupling reaction

在圖6 所示的分析中，可以看到預測產率和觀察產率的散點，其擬合出的直線應該滿足如下規律：1）斜率為1。兩種產率結果良好的擬合必然是擬合直線斜率為1。2）截距為0。截距能在一定程度上反映出預測結果相對于觀察值之間的偏移。3）良好的散點疏密程度。真正有效的擬合必然是散點大部分聚集于擬合直線周圍。

在隨機劃分（70/30，即訓練和測試數據分別占70%和30%）數據集上，所提出的基于Boltzmann 分布的加權平均占位描述符可以使隨機森林等機器學習方法得到較高的準確率，證明它對催化劑選擇性的表達良好。其中，線性回歸（圖7（b））和SVM（圖7（c））的截距大，斜率也不接近1，因而效果較差。K 近鄰的斜率和截距都很好，但其中的散點大部分分散于距離擬合直線很遠的區域，因而效果很差，自適應增強算法（Adaptive Boosting，AdaBoost）也存在同樣的問題。表現好的是隨機森林（圖7（e））和本文的CNN 特征融合模型（圖7（f）），而從散點的疏密程度來看，CNN 特征融合模型是最好的。圖7 中使用了對映體過量（Enantiomeric Excess，ee）值指標來評估模型性能。

圖7 使用基于量子力學特征與分子指紋的描述符的6種方法在N，S-縮醛反應上的性能對比Fig.7 Performance comparison of six methods using descriptors based on quantum mechanical features and molecular fingerprints in N，S-acetal formation

CNN 特征融合模型在隨機劃分數據上的R2比Ahneman［21］提高了3 個百分點，在additive test 2 數據集上比Schwaller［23］提高了1 個百分點，在additive test 4 數據集上比Ahneman［21］提高了2 個百分點（表3）。在4 個添加劑數據集上的平均性能比MFF［18］提高了8.5 個百分點，這說明了模型在添加劑水平不同的數據集中具有適用性。

表3 C-N偶聯反應上的測試結果（R2）Tab.3 Prediction results（R2）on C-N coupling reaction

在手性磷酸催化的N，S-縮醛反應篩選中，使用的反應選擇性評價指標是結合自由能差ΔΔG，實踐中發現它比ee值的預測更有挑戰性。ΔΔG的表達式為：

其中：R是理想氣體常數，T是常溫（取23℃）。

與前面的討論類似，使用結合量子力學性質與分子指紋的描述符后，隨機森林等傳統機器學習方法表現出令人鼓舞的效果，這證明了描述符的有效性。在與已發表工作對比中，GCN 特征融合模型在多個數據集上的平均絕對誤差（Mean Absolute Error，MAE）比MFF 降低1.2 個百分點，接近Zahrt 等［22］。在大多數數據集上取得了比基線方法更小的MAE（表4）。

表4 N，S-縮醛反應上的測試結果（MAE）Tab.4 Prediction results（MAE）on N，S-acetal formation

從已發表工作［18，21-22，35-36］可知，選擇和構建分子描述符在有機反應預測中非常重要。實驗結果表明，通過合適的特征融合模型將多種來源的描述符如量子力學性質、3D 空間性質、SMILES 和分子指紋等綜合起來，在產率預測和選擇性預測等篩選工作中都可以獲得出色的預測性能，也是將深度學習中一些模型引入該領域的可選方案。

3.4 催化劑篩選結果

有效的催化劑的篩選可以大幅加快催化劑優化的進程。在手性磷酸催化的N，S-縮醛反應篩選中，通過對訓練集之外的反應進行篩選，發現雖然訓練集中包含大量的中、低反應選擇性的催化劑和反應，高于95%的反應和催化劑很少，但GCN 特征融合模型卻能把高選擇性的催化劑大部分都預測出來，表現出令人鼓舞的外推能力。如圖8 所示，GCN 特征融合模型可以篩選出有價值和高選擇性的催化劑，且預測的ee 值與觀察值非常接近。這種能力對于不對稱催化的藥 物合成反應的優化很有意義。

圖8 使用GCN特征融合模型篩選出的一部分催化劑Fig.8 Some catalysts screened by GCN feature fusion model

4 結語

本文通過分析與實驗證明了使用深度學習方法融合多模態數據對藥物合成反應進行虛擬篩選是可行的，并提出一種基于Boltzmann 分布的3D 描述符，它按構象能量計算權重從而將多個構象加權平均，來近似反應發生時構象的概率分布。針對單一模態分子特征表示的不足，提出一種融合3D描述符和DFT 分子性質的CNN 模型，用于4-甲基苯胺與芳基鹵化物的C-N 偶聯反應的產率預測；同時還提出一種基于DFT 性質的圖卷積網絡，該網絡將原子的QM 性質作為節點的輸入特征，使用RDKit 產生的鄰接矩陣來構建圖結構，并利用分支結構將分子指紋特征融合進來。最后，在不對稱N，S-縮醛反應上驗證了這種基于圖卷積的特征融合方法的有效性。因此，得出結論：使用量子力學性質構建的分子描述符，特別是3D 描述符是預測藥物合成反應的有效特征表示，而且基于多模態融合的深度學習可以將它與SMILES 以及分子指紋等結合起來，協同發揮作用。