孫曉飛，朱靜遠(yuǎn)，陳 斌，游恒志*

（1.中國(guó)科學(xué)院 成都計(jì)算機(jī)應(yīng)用研究所，成都 610041；2.中國(guó)科學(xué)院大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，北京 100049；3.哈爾濱工業(yè)大學(xué)（深圳）理學(xué)院，廣東 深圳 518055；4.哈爾濱工業(yè)大學(xué)（深圳）人工智能研究院，廣東 深圳 518055）

0 引言

藥物合成中的有機(jī)化學(xué)反應(yīng)的發(fā)現(xiàn)依賴于實(shí)踐經(jīng)驗(yàn)和化學(xué)機(jī)理支配的“化學(xué)直覺”，實(shí)驗(yàn)人員試圖定性地識(shí)別有機(jī)化學(xué)反應(yīng)中的模式，以確定反應(yīng)產(chǎn)物和反應(yīng)效率。然而，這種方法受到了很多因素的限制，包括反應(yīng)的復(fù)雜性、活性懸崖、對(duì)機(jī)理理解的缺乏以及人工處理大數(shù)據(jù)的艱難?；谟?jì)算機(jī)的虛擬篩選［1-5］已經(jīng)成為吸引化學(xué)家的重要方案，主要因?yàn)樗恍枰獧C(jī)理的理解，化合物結(jié)構(gòu)可以用分子性質(zhì)的數(shù)值表示來表征，從而量化數(shù)以千計(jì)的候選分子的化學(xué)性質(zhì)。在實(shí)驗(yàn)和文獻(xiàn)數(shù)據(jù)的基礎(chǔ)上，虛擬篩選可以通過計(jì)算機(jī)模型來對(duì)藥物合成反應(yīng)的結(jié)果和催化劑的選擇性進(jìn)行量化。

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

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

圖1 對(duì)藥物合成反應(yīng)進(jìn)行虛擬篩選的流程Fig.1 Flow of virtual screening of drug synthesis reactions

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

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

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

在C-N 偶聯(lián)反應(yīng)產(chǎn)率預(yù)測(cè)中，采用格點(diǎn)精度為1 埃（0.1 nm）的加權(quán)平均占位描述符，相較于文獻(xiàn)［21］和文獻(xiàn)［18］提高了3和2個(gè)百分點(diǎn)，對(duì)僅使用DFT計(jì)算的描述符有明顯優(yōu)勢(shì)。在N，S-縮醛反應(yīng)的對(duì)映選擇性預(yù)測(cè)中，在隨機(jī)劃分的數(shù)據(jù)集上誤差從文獻(xiàn)［22］的0.152降到0.147，驗(yàn)證了方法的有效性。

1 相關(guān)工作

1.1 三維分子描述符

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

1.2 藥物合成反應(yīng)篩選中的多模態(tài)數(shù)據(jù)融合

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

多模態(tài)融合負(fù)責(zé)聯(lián)合多個(gè)模態(tài)的數(shù)據(jù)進(jìn)行分類或者回歸，在深度學(xué)習(xí)的很多領(lǐng)域應(yīng)用廣泛。它還有其他常見的別名，如多傳感器融合和多源信息融合。單模態(tài)的機(jī)器學(xué)習(xí)將信息表示為計(jì)算機(jī)可以處理的數(shù)值向量并進(jìn)一步抽象為更高層的特征向量，而多模態(tài)機(jī)器學(xué)習(xí)可以通過利用多模態(tài)數(shù)據(jù)之間的互補(bǔ)性，學(xué)習(xí)到更好的特征表示，如圖2 所示。常見的機(jī)器學(xué)習(xí)模型都可以用于多模態(tài)融合［31］。

圖2 多模態(tài)數(shù)據(jù)融合Fig.2 Multimodal data fusion

1.3 圖卷積神經(jīng)網(wǎng)絡(luò)

可以使用圖卷積神經(jīng)網(wǎng)絡(luò)來建模分子相關(guān)的問題，將一個(gè)化合物構(gòu)建為一個(gè)圖形，它的節(jié)點(diǎn)代表原子，連接它們的邊代表鍵。在圖卷積（圖3（a））中，對(duì)于一個(gè)節(jié)點(diǎn)，將特征和鄰居饋送到兩個(gè)密集層中，然后添加密集層的輸出作為節(jié)點(diǎn)的新特征。對(duì)具有相同度的節(jié)點(diǎn)計(jì)算新特征時(shí)共享權(quán)重。在一個(gè)化合物中，如果一個(gè)原子a總共有n個(gè)鄰居，那么它經(jīng)過圖卷積后的新特征可以表示為：

其中：Wa為節(jié)點(diǎn)a的權(quán)重；Wr是相鄰節(jié)點(diǎn)的權(quán)重；b是偏差，σ是激活函數(shù)。從節(jié)點(diǎn)a出發(fā)的箭頭表示原子a及其相鄰原子的密集層，權(quán)重分別為Wa和Wr。

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

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

將這三種操作按圖3 中黑色箭頭所示的順序進(jìn)行組合，就可以構(gòu)成不對(duì)稱反應(yīng)催化劑篩選網(wǎng)絡(luò)中的圖卷積網(wǎng)絡(luò)模塊。

圖3 圖卷積神經(jīng)網(wǎng)絡(luò)Fig.3 Graph convolution neural network

2 本文方法

本文提出了一種基于加權(quán)平均占位的3D描述符，它在三維空間中整合分子級(jí)別性質(zhì)和原子級(jí)別性質(zhì)，為機(jī)器學(xué)習(xí)模型提供良好的輸入特征，并將這種方法用于鈀催化的Buchwald-Hartwig偶聯(lián)反應(yīng)的篩選；還提出一種基于圖卷積的多模態(tài)模型，通過將量子力學(xué)性質(zhì)和分子指紋融合起來來篩選N，S-縮醛反應(yīng)中的手性磷酸催化劑。下面介紹分子描述符合網(wǎng)絡(luò)的結(jié)構(gòu)的一些實(shí)現(xiàn)細(xì)節(jié)。

2.1 基于Boltzmann分布的加權(quán)平均占位描述符

在特征選擇和模型建立過程中，需要考慮描述符對(duì)有機(jī)反應(yīng)的適應(yīng)性，即要尋找一組不受限于特定機(jī)理假設(shè)的描述符，以區(qū)分和描述反應(yīng)之間的差異?；谶@種考慮，在Buchwald-Hartwig 偶聯(lián)反應(yīng)的產(chǎn)率預(yù)測(cè)中引入一種3D 描述符（圖4（a）），它首先根據(jù)結(jié)構(gòu)優(yōu)化后的分子3D坐標(biāo)構(gòu)建出一個(gè)三維網(wǎng)格，將分子置于網(wǎng)格中心并采用自開發(fā)的分子對(duì)齊算法進(jìn)行對(duì)齊，然后統(tǒng)計(jì)每個(gè)網(wǎng)格點(diǎn)上原子的量子力學(xué)性質(zhì)。考慮到單一的構(gòu)象不能很好地表示反應(yīng)發(fā)生時(shí)構(gòu)象的復(fù)雜情況，使用Boltzmann 分布計(jì)算不同構(gòu)象所占比例，生成加權(quán)平均的3D構(gòu)象描述符。

其中：p是構(gòu)象所占比例，k、i是構(gòu)象的序號(hào)，E是構(gòu)象的能量，R是理想氣體常數(shù)，T是開爾文溫度。di是單個(gè)構(gòu)象的描述符，d是經(jīng)過加權(quán)計(jì)算后的分子描述符。

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

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

2.2 反應(yīng)產(chǎn)率篩選網(wǎng)絡(luò)架構(gòu)

融合多種來源的分子特征數(shù)據(jù)，對(duì)藥物合成相關(guān)的有機(jī)反應(yīng)進(jìn)行篩選，是本次工作的重要目標(biāo)?；诩訖?quán)平均占位的3D 描述符可以將原子級(jí)別的量子力學(xué)特征與立體空間特征融合起來。分子整體的量子力學(xué)性質(zhì)如HOMO、LOMO 等也是與反應(yīng)機(jī)理相關(guān)的重要特征，在3D描述符中并不能很好表達(dá)這種特征，在此將多個(gè)分子級(jí)別特征表示為一維向量描述符。3D 描述符數(shù)據(jù)的維度巨大，采用卷積操作構(gòu)建網(wǎng)絡(luò)以減少過擬合，而對(duì)于分子級(jí)別特征直接使用密集層構(gòu)建網(wǎng)絡(luò)，然后將兩種特征進(jìn)行融合（圖4（b））。損失函數(shù)的形式如下：

其中：x是樣本，y是化學(xué)實(shí)驗(yàn)得到的觀察值。

2.3 不對(duì)稱反應(yīng)中催化劑篩選的網(wǎng)絡(luò)架構(gòu)

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

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

3 實(shí)驗(yàn)與結(jié)果分析

3.1 有機(jī)反應(yīng)數(shù)據(jù)和評(píng)價(jià)指標(biāo)

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

表1 Buchwald-Hartwig反應(yīng)的數(shù)據(jù)集劃分Tab.1 Splitting of Buchwald-Hartwig reaction dataset

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

表2 N，S-縮醛反應(yīng)的數(shù)據(jù)集劃分Tab.2 Splitting of N，S-acetal formation dataset

評(píng)價(jià)指標(biāo)采用與已報(bào)道工作相同的R2（R-squared）和RMSE（Root Mean Square Error），其中R2指標(biāo)計(jì)算方式為：

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

3.2 實(shí)驗(yàn)細(xì)節(jié)

在對(duì)藥物合成相關(guān)反應(yīng)的數(shù)據(jù)進(jìn)行處理時(shí)，出于數(shù)據(jù)一致性的考慮，對(duì)于分子統(tǒng)一使用了Gaussian 計(jì)算出所需要的性質(zhì)。為了減少人工數(shù)據(jù)分析，開發(fā)了自動(dòng)化軟件將分子提交給Gaussian 進(jìn)行計(jì)算，并提取和解析結(jié)果。在計(jì)算催化劑的性質(zhì)時(shí)，使用ChemDraw、Gaussian 創(chuàng)建初始的分子結(jié)構(gòu)，或者從Cambridge Structural Database［34］導(dǎo)入。經(jīng)過構(gòu)象搜索從催化劑的多個(gè)構(gòu)象中獲得能量最低的構(gòu)象以及不同能量的構(gòu)象用于生成加權(quán)平均構(gòu)象描述符。計(jì)算中未考慮溶劑影響，停止優(yōu)化的條件為Gaussian 的默認(rèn)收斂標(biāo)準(zhǔn)。使用自研的程序?qū)Υ呋瘎┻M(jìn)行對(duì)齊，以提供更靈活和自動(dòng)化的對(duì)齊方式。

在AI 模型設(shè)計(jì)和訓(xùn)練部分，所有的程序都是基于Pytorch 和Scikit-learn 實(shí)現(xiàn)，并將模型預(yù)測(cè)結(jié)果與已報(bào)道文獻(xiàn)結(jié)果進(jìn)行對(duì)比。鑒于數(shù)據(jù)規(guī)模較小，模型訓(xùn)練中使用了1 個(gè)NVIDIA GeForce GTX 1070，顯存為8 GB。在數(shù)據(jù)預(yù)處理中，由于分子中原子數(shù)目不相同，故按列（特征）對(duì)數(shù)據(jù)進(jìn)行統(tǒng)一長(zhǎng)度和縮放處理，這有利于減少無效的數(shù)據(jù)、縮短描述符長(zhǎng)度從而減少網(wǎng)絡(luò)參數(shù)和提升效果。

3.3 主要結(jié)果

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

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

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

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

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

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

表3 C-N偶聯(lián)反應(yīng)上的測(cè)試結(jié)果（R2）Tab.3 Prediction results（R2）on C-N coupling reaction

在手性磷酸催化的N，S-縮醛反應(yīng)篩選中，使用的反應(yīng)選擇性評(píng)價(jià)指標(biāo)是結(jié)合自由能差ΔΔG，實(shí)踐中發(fā)現(xiàn)它比ee值的預(yù)測(cè)更有挑戰(zhàn)性。ΔΔG的表達(dá)式為：

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

與前面的討論類似，使用結(jié)合量子力學(xué)性質(zhì)與分子指紋的描述符后，隨機(jī)森林等傳統(tǒng)機(jī)器學(xué)習(xí)方法表現(xiàn)出令人鼓舞的效果，這證明了描述符的有效性。在與已發(fā)表工作對(duì)比中，GCN 特征融合模型在多個(gè)數(shù)據(jù)集上的平均絕對(duì)誤差（Mean Absolute Error，MAE）比MFF 降低1.2 個(gè)百分點(diǎn)，接近Zahrt 等［22］。在大多數(shù)數(shù)據(jù)集上取得了比基線方法更小的MAE（表4）。

表4 N，S-縮醛反應(yīng)上的測(cè)試結(jié)果（MAE）Tab.4 Prediction results（MAE）on N，S-acetal formation

從已發(fā)表工作［18，21-22，35-36］可知，選擇和構(gòu)建分子描述符在有機(jī)反應(yīng)預(yù)測(cè)中非常重要。實(shí)驗(yàn)結(jié)果表明，通過合適的特征融合模型將多種來源的描述符如量子力學(xué)性質(zhì)、3D 空間性質(zhì)、SMILES 和分子指紋等綜合起來，在產(chǎn)率預(yù)測(cè)和選擇性預(yù)測(cè)等篩選工作中都可以獲得出色的預(yù)測(cè)性能，也是將深度學(xué)習(xí)中一些模型引入該領(lǐng)域的可選方案。

3.4 催化劑篩選結(jié)果

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

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

4 結(jié)語

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