考慮選擇性和反應(yīng)速率的多目標制藥反應(yīng)溶劑設(shè)計

趙紅慶，劉奇磊，張磊，董亞超，都健

（大連理工大學化工學院，化工系統(tǒng)工程研究所，遼寧大連116024）

引 言

藥物研發(fā)是一個風險大、周期長、成本高的過程，在這個以間歇合成為主要方法的過程中溶劑的使用量占總物耗的80%～90%[1]。過去，藥物化學家們致力于開發(fā)出可靠的合成過程，希望可以以高純度、高收率生產(chǎn)活性藥物成分（active pharmaceutical ingredient, API），同時保證效率、可操作性、安全性以及對環(huán)境的低影響性。然而，溶劑的作用一直沒有受到重視[2]。溶劑在不同階段可以起到不同的作用：①作為反應(yīng)介質(zhì)或同時作為試劑參與反應(yīng)；②影響反應(yīng)的選擇性，從而減少副產(chǎn)物的形成，實現(xiàn)高質(zhì)量合成；③改善產(chǎn)物的分離效果[2]。因此，為不同的階段尋找到合適的溶劑對于實現(xiàn)藥物的清潔高效生產(chǎn)具有重要的意義。

傳統(tǒng)的溶劑篩選過程一般采用實驗方法實現(xiàn)。然而，此類方法需依托實驗平臺，搜索范圍有限且費時費力。為優(yōu)化實驗，F(xiàn)isher[3]提出通過對實驗進行合理安排，以較小的實驗規(guī)模、較短的實驗周期和較低的成本獲得理想的實驗結(jié)果，并得出科學結(jié)論。但是，這種方法并不能從根本上解決問題。計算機輔助分子設(shè)計（computer-aided molecular design,CAMD）方法的提出，為以上問題提供了一個更有前景和系統(tǒng)化的解決辦法[4]。CAMD 方法經(jīng)過多年的發(fā)展，目前已應(yīng)用于工業(yè)中分離過程的溶劑設(shè)計，如液液萃取[4-6]、萃取蒸餾[7]等。

近年來，CAMD 的應(yīng)用領(lǐng)域進一步拓展到了反應(yīng)溶劑設(shè)計[8-12]。在液-液均相有機反應(yīng)中，溶劑作為介質(zhì)與反應(yīng)中各物質(zhì)產(chǎn)生相互作用，進而改變化學反應(yīng)的反應(yīng)速率、反應(yīng)平衡、溶解度、穩(wěn)定性甚至反應(yīng)機理。換言之，化學反應(yīng)中的熱力學和動力學可通過選擇適當溶劑來進行控制[13-15]。近年來，多位學者根據(jù)溶劑對反應(yīng)速率的影響，成功地利用CAMD 方法進行了反應(yīng)溶劑的設(shè)計。Foli? 等[8]通過多參數(shù)溶劑變色方程將少量溶劑性質(zhì)（溶劑變色參數(shù)、內(nèi)聚能密度參數(shù)等）和反應(yīng)速率常數(shù)實驗值的對數(shù)值進行關(guān)聯(lián)，并以速率常數(shù)最大化為目標函數(shù)進行了反應(yīng)溶劑設(shè)計。該方法只需使用少量溶劑（6～8 種）進行實驗，利用基團貢獻法（group contribution, GC）即可實現(xiàn)反應(yīng)溶劑的設(shè)計。Zhou等[9]使用由量化計算得到的溶劑描述符和少量實驗數(shù)據(jù)，回歸得到與反應(yīng)速率常數(shù)關(guān)聯(lián)的代理模型，并成功應(yīng)用在Diels-Alder（DA）反應(yīng)中。在構(gòu)建反應(yīng)動力學模型時，反應(yīng)速率常數(shù)除了通過實驗獲取，還可以通過量化計算獲得。2013 年Struebing等[10]提出了一種結(jié)合QM-CAMD 的反應(yīng)溶劑設(shè)計方法，其中反應(yīng)速率常數(shù)使用QM（quantum mechanics）與過渡態(tài)理論計算得到。該方法被成功應(yīng)用于雙分子親核取代（SN2）反應(yīng)中，避免了煩瑣耗時的實驗操作，節(jié)省了大量人力物力。然而，以上方法在描述符的選擇上尚存在一些問題：①部分描述符仍為實驗測量值，如溶劑化變色參數(shù)[8]等；②盡管量化計算得到的描述符避免了煩瑣的實驗，但是使用此描述符的模型應(yīng)用范圍有限，仍需通過耗時的計算獲取。因此，Liu 等[11]基于傳統(tǒng)過渡態(tài)理論（transition state theory,TST），結(jié)合類導體屏蔽片段活度系數(shù)模型（COSMO-SAC），通過篩選算法得到新的描述符，包括無限稀釋活度系數(shù)γ∞、氫鍵供體XDON、氫鍵受體XACC、表面張力σ 等，這些描述符可以通過量化計算或基團貢獻法快速獲得。同時，由于該方法從過渡態(tài)理論出發(fā)，經(jīng)過嚴格的熱力學推導，因此具有更為廣泛的適用范圍和更高的預測精度。但是，以上研究均是將反應(yīng)速率最大化作為優(yōu)化目標，并未考慮反應(yīng)的選擇性。近年來，有學者開始考慮使用CAMD 方法對反應(yīng)的選擇性進行研究，即以最大化主產(chǎn)物與副產(chǎn)物產(chǎn)率的差值為優(yōu)化目標。Adjiman等[12]以最大化反應(yīng)速率與選擇性為目標，提出了一種基于溶劑變色方程的CAMD 方法，并通過芳香親核反應(yīng)進行了驗證。但是此方法中所使用的描述符和反應(yīng)速率常數(shù)仍需要實驗數(shù)據(jù)進行擬合。

因此，本文提出了QM-CAMD 結(jié)合的反應(yīng)溶劑設(shè)計方法，基于Liu等[11]提出的描述符構(gòu)建反應(yīng)動力學模型，并通過QM 計算代替實驗來獲取不同溶劑作用下的反應(yīng)速率常數(shù)，從而對模型中的參數(shù)進行擬合，進而以最大化反應(yīng)速率常數(shù)和選擇性為目標建立包括目標函數(shù)和溶劑結(jié)構(gòu)及性質(zhì)約束的混合整數(shù)非線性規(guī)劃（MINLP）模型進行反應(yīng)溶劑設(shè)計，實現(xiàn)反應(yīng)溶劑設(shè)計的多目標優(yōu)化。本文以嘧啶類化合物的合成過程為例，建立了溶劑設(shè)計模型，采用分解式算法對該模型進行求解，并對結(jié)果進行了分析。

1 QM-CAMD設(shè)計框架

本文的研究框架如圖1所示。第一步選擇反應(yīng)體系，同時根據(jù)介電常數(shù)大小選擇初始溶劑用于后續(xù)反應(yīng)速率常數(shù)及描述符計算；第二步對于選定的反應(yīng)體系經(jīng)過反應(yīng)物結(jié)構(gòu)優(yōu)化、過渡態(tài)搜索、內(nèi)稟反應(yīng)坐標（IRC）驗證過程進行反應(yīng)機理研究，然后根據(jù)過渡態(tài)理論計算反應(yīng)速率常數(shù)同時計算反應(yīng)物和過渡態(tài)分子在溶劑環(huán)境下的無限稀釋活度系數(shù)以及溶劑對應(yīng)的描述符；第三步反應(yīng)速率常數(shù)與描述符經(jīng)過多元線性回歸方法回歸得到反應(yīng)動力學模型；第四步結(jié)合目標函數(shù)與約束條件共同組成MINLP 模型，并采用分解式算法對模型進行優(yōu)化求解得到可行分子列表；最后對評價好的溶劑進行實驗驗證，最終得到性能更優(yōu)的反應(yīng)溶劑。

圖1 QM-CAMD設(shè)計方法Fig.1 QM-CAMD design method

1.1 反應(yīng)體系的選擇

嘧啶環(huán)是生物分子脫氧核糖核酸（DNA）和三磷酸腺苷（ATP）的重要結(jié)構(gòu)單元，含有嘧啶結(jié)構(gòu)單元的化合物通常具有廣泛的藥理活性，因此對含嘧啶環(huán)分子的研究一直是藥物化學工作者的關(guān)注重點[16]。其中，2,4-二氯-5-硝基嘧啶是一種非常有用的原料，可以用來合成（或經(jīng)進一步合成轉(zhuǎn)化）大量具有生物活性的含嘧啶結(jié)構(gòu)的藥物，例如腺苷A2A受體拮抗劑[17]、人類激酶抑制劑[18]、二芳基嘧啶（DAPYs）[19]等。本文選擇DAPYs合成的第一步反應(yīng)為研究對象。DAPYs 是HIV-1 非核苷逆轉(zhuǎn)錄酶（RT）抑制劑中的一種，可用于治療艾滋病，該合成過程的反應(yīng)方程式如圖2所示。

圖2 2,4-二氯-5-硝基嘧啶與對氨基苯腈的化學反應(yīng)方程式Fig.2 Reaction of 2,4-dichloro-5-nitropyrimidine and 4-aminobenzonitrile

由于2,4-二氯-5-硝基嘧啶存在兩個親核反應(yīng)位點，當胺類與2,4-二氯-5-硝基嘧啶發(fā)生芳香親核反應(yīng)（SNAr）時，C-4 和C-2 產(chǎn)物比例約為9∶1～19∶1[20]。其中C-4 產(chǎn)物為所需產(chǎn)物，然而產(chǎn)物中依然存在著少量C-2副產(chǎn)物。因此對于具有多個親核反應(yīng)位點的反應(yīng)，提高選擇性有利于增加產(chǎn)率并減少副產(chǎn)物，實現(xiàn)更加清潔高效的合成過程。

1.2 反應(yīng)機理研究及反應(yīng)速率常數(shù)計算

1935 年，Eyring[21]提出了過渡態(tài)理論，其基本假設(shè)是由反應(yīng)物產(chǎn)生的過渡態(tài)（即活化絡(luò)合物）處于一種準平衡狀態(tài)，并且任何反應(yīng)體系沿著反應(yīng)坐標一次通過過渡態(tài)，永不復返。過渡態(tài)理論經(jīng)過多年的發(fā)展，已經(jīng)取得了巨大成就。通過過渡態(tài)理論研究化學反應(yīng)需要確定反應(yīng)路徑，即反應(yīng)機理。針對本文所研究的反應(yīng)體系，其存在兩種可能的反應(yīng)機理：第一種是自1870 年以來，長期以來被廣泛接受的機理，即反應(yīng)過程中生成Meisenheimer 絡(luò)合物中間體的兩步反應(yīng)機理[22-23]；第二種是一步反應(yīng)的連接機理，即反應(yīng)過程不產(chǎn)生中間體，親核試劑的進攻和離去基團的離去這兩者同時發(fā)生[24]。

當反應(yīng)機理確定后，即可根據(jù)過渡態(tài)理論進行反應(yīng)速率常數(shù)的計算，對于氣相雙分子反應(yīng)方程，如式(1)所示，其反應(yīng)速率常數(shù)可根據(jù)艾林方程（Eyring equation）計算，如式(2)所示，其中ΔG?為氣相活化Gibbs自由能。

對于液相雙分子反應(yīng)，ΔG?為液相中的活化Gibbs 自由能，其值等于真空下的Gibbs 自由能與溶劑化自由能之和，如式(3)所示。

將式(3)代入式(2)可得式(4)：

在式(4)中，在標準狀態(tài)下，除真空下反應(yīng)能壘ΔGIG和反應(yīng)的溶劑化自由能Δ?ΔGsolv外均為常數(shù)，因此只要獲得反應(yīng)物和過渡態(tài)分子在真空和溶劑下的能量值即可獲得反應(yīng)速率常數(shù)k。反應(yīng)機理及反應(yīng)速率常數(shù)計算步驟可總結(jié)如下。

（1）使用Gaussian 09[25]軟件對反應(yīng)物分子進行結(jié)構(gòu)優(yōu)化和振動分析，并基于合理的過渡態(tài)初猜結(jié)構(gòu)搜索真實的過渡態(tài)結(jié)構(gòu)，然后進行內(nèi)稟反應(yīng)坐標（intrinsic reaction coordinate,IRC）[26]驗證確保過渡態(tài)結(jié)構(gòu)正確；

（2）基于第一步優(yōu)化得到的反應(yīng)物和過渡態(tài)的分子結(jié)構(gòu)，對反應(yīng)物和過渡態(tài)分子使用ORCA[27]量化計算軟件，采用更高級別的基組與泛函PWPB95-D3(BJ)/def2-TZVPP 進行單點能計算，獲得它們在真空下更高精度的電子能量，隨后將過渡態(tài)與反應(yīng)物分子能量加上對應(yīng)的Gibbs 熱校正量后作差得到真空下反應(yīng)能壘ΔGIG；

（3）基于第一步的結(jié)構(gòu)在Gaussian 09 中使用M052X/6-31G(d)和基于密度的溶劑模型（solvation model based on density,SMD）[28]方法分別計算分子結(jié)構(gòu)在真空和溶劑下的電子能量，隨后將二者作差獲得反應(yīng)的溶劑化自由能Δ?ΔGsolv。

（4）更換溶劑，重復步驟(3)，獲得不同溶劑下的溶劑化自由能，并分別代入式(4)獲得不同溶劑下的反應(yīng)速率常數(shù)k。

1.3 構(gòu)建反應(yīng)動力學模型

通過以上計算可以得到不同溶劑作用下的反應(yīng)速率常數(shù)k。在此基礎(chǔ)上，通過選取適當?shù)娜軇┟枋龇?，利用?shù)據(jù)回歸方法建立溶劑描述符與反應(yīng)速率常數(shù)間的代理模型。本文建立的反應(yīng)動力學代理模型選取的溶劑描述符包括溶劑環(huán)境下反應(yīng)物與過渡態(tài)分子的無限稀釋活度系數(shù)γ∞、溶劑分子的氫鍵供體XDON、氫鍵受體XACC以及表面張力σ。其中，描述符分別采用以下兩種策略獲?。孩僭诮⒎磻?yīng)動力學模型時為了提高模型的準確性，反應(yīng)物和過渡態(tài)分子的γ∞通過SCM-ADF 軟件[29]基于2016-ADF Chen 方 法 的COSMO-SAC 模 型[30]計 算 得到，其余描述符均來自數(shù)據(jù)庫；②在CAMD 設(shè)計過程中，無限稀釋活度系數(shù)γ∞通過GC-COSMO[11]方法獲得，XDON、XACC、σ 等描述符通過基團貢獻法[4]進行快速預測。隨后對已獲得的反應(yīng)速率常數(shù)的對數(shù)值和描述符以R2最大為目標進行多元線性回歸得到反應(yīng)動力學模型，如式(5)所示。

1.4 反應(yīng)溶劑設(shè)計模型

本文利用MINLP 優(yōu)化模型建立了反應(yīng)溶劑設(shè)計模型。模型包括兩個目標函數(shù)，分別為4-位取代的反應(yīng)速率常數(shù)f1= lg k4以及選擇性f2= lg k4-lg k2，具體模型如式(6)～式（11）所示。

（1）目標函數(shù)（其中l(wèi)gk2與lgk4的表達式由1.2 節(jié)及1.3節(jié)中描述的方法回歸得到）

（2）分子結(jié)構(gòu)約束（八隅體規(guī)則、價鍵規(guī)則等[31],詳細方程見文獻[31]中方程(5)～方程(22)）

（3）溶劑性質(zhì)約束（線性性質(zhì)）：如XDON、XACC、σ、熔沸點、毒性、溶解度系數(shù)等

（4）溶劑性質(zhì)約束（非線性性質(zhì)）：γ∞[32]（詳細見文獻[32]中方程(1)～方程(12)）

其中，熔點和沸點約束是為了保證溶劑在反應(yīng)溫度下為液態(tài)，毒性約束了設(shè)計得到的溶劑的綠色安全性。由于該優(yōu)化模型中存在大量非線性方程，給模型的求解帶來了很大的困難。因此有必要尋找一種有效的求解方法。

1.5 分解式算法

本文采用分解式算法[33]（decomposition-based algorithm,DA），將難以求解的MINLP 問題分解為多個子問題，在GAMS[34]軟件中使用BARON 求解器[35]編程求解。分解算法的具體求解原理如圖3所示。

（1）由選擇的初始基團根據(jù)結(jié)構(gòu)約束通過數(shù)學算法生成大量結(jié)構(gòu)可行的分子結(jié)構(gòu)（數(shù)量為N1）；

（2）對上一步生成的分子結(jié)構(gòu)進行線性性質(zhì)的約束，通過基團貢獻法計算熔點、沸點、毒性等性質(zhì)，篩選得到滿足線性性質(zhì)的分子結(jié)構(gòu)（數(shù)量為N2，N2＜N1）。

圖3 分解式算法原理Fig.3 Principle of decomposition-based algorithm

（3）使用GC-COSMO 方法分別計算N2個溶劑分子與溶質(zhì)（反應(yīng)物與過渡態(tài)）的γ∞，通過反應(yīng)動力學模型[式(5)]得到目標函數(shù)f1和f2的值。

（4）在經(jīng)過分解式算法計算后，按照目標函數(shù)f1以及f2的取值獲得可行溶劑分子列表，并進行評價。

2 結(jié)果分析與討論

2.1 反應(yīng)機理

由于本文研究體系存在兩種可能的機理，因此本文在真空條件以及不同溶劑作用下通過密度泛函理論（density functional theory, DFT）計算來探索在2-位和4-位可能發(fā)生的兩種反應(yīng)機理。經(jīng)過計算發(fā)現(xiàn)，按照兩步反應(yīng)機理反應(yīng)時，基團離去過程的過渡態(tài)無法定位。而使用一步反應(yīng)機理可以準確定位過渡態(tài)并成功通過IRC 驗證，其過渡態(tài)結(jié)構(gòu)如圖4所示。

根據(jù)以上結(jié)果，本文對于速率常數(shù)的計算均基于一步反應(yīng)機理。反應(yīng)過程中的相對Gibbs 自由能如表1所示?？梢钥闯?-位取代的活化自由能低于2-位，從而證明同樣條件下4-位發(fā)生親核取代的可能性更高，這與實驗中的主產(chǎn)物為4-位取代物相符合。為得到不同溶劑作用下的速率常數(shù)值，本文從Gaussian 09 內(nèi)置溶劑中選取了80種溶劑，覆蓋了多種極性和非極性溶劑，分別進行了速率常數(shù)的計算。

2.2 反應(yīng)動力學模型

在DFT 計算的基礎(chǔ)上，對得到的反應(yīng)速率常數(shù)的對數(shù)值和溶劑描述符進行多元線性回歸，分別得到4-位取代（k4）和2-位取代（k2）的反應(yīng)動力學模型，如式(12)、式(13)所示。

其中r2=0.934。

圖4 真空中的反應(yīng)過渡態(tài)結(jié)構(gòu)Fig.4 Transition state structure of reaction in vacuum

表1 真空中反應(yīng)物、過渡態(tài)與產(chǎn)物的相對自由能Table 1 Relative Gibbs free energies of reactants，transition states，and products in a vacuum

其中r2=0.855。

以上兩個模型回歸效果如圖5所示。兩個模型預測值的絕對誤差絕大部分在±1 以內(nèi)，表明該反應(yīng)動力學模型預測效果良好，預測結(jié)果可靠性強，可以用于后續(xù)的CAMD設(shè)計。

2.3 CAMD設(shè)計結(jié)果

根據(jù)圖1 的計算流程，本文預選了52 個基團，通過表2所示的分子結(jié)構(gòu)及溶劑性質(zhì)約束建立溶劑設(shè)計模型并進行求解，得到了可行分子列表，最終的結(jié)果如圖6 所示。從圖6 可知使主反應(yīng)的反應(yīng)速率常數(shù)增幅大的溶劑，會使反應(yīng)選擇性下降；選擇性大幅度提升的溶劑，其對應(yīng)的反應(yīng)速率常數(shù)下降。因此，兩個目標函數(shù)存在矛盾。

對設(shè)計結(jié)果中的四種溶劑作詳細分析，如表3所示。由表3 的結(jié)果可知，0 號溶劑（即圖6 中放大的“×”所代表的溶劑）四氫呋喃（tetrahydrofuran,THF）[19]為文獻中已經(jīng)報道的溶劑，1～4 號溶劑為CAMD 設(shè)計得到的溶劑。1 號溶劑對反應(yīng)的選擇性雖然提高了兩個數(shù)量級，但是在速率常數(shù)上起到了抑制作用；3 號溶劑對速率常數(shù)提升明顯，但是同時對2-位取代也起到了促進作用，導致對反應(yīng)選擇性的提升遜于2 號溶劑；而4 號溶劑雖然對主副反應(yīng)的反應(yīng)速率常數(shù)提升明顯，但對副反應(yīng)提升更大，導致選擇性下降。因此，對反應(yīng)速率和選擇性進行綜合考慮后，可以認為2號溶劑效果最好，其反應(yīng)速率常數(shù)相對THF 提升了62.8%，同時選擇性也提高了一個數(shù)量級。在后續(xù)的工作中，將對以上設(shè)計結(jié)果進行進一步實驗驗證。

表2 MINLP模型中的結(jié)構(gòu)和性質(zhì)約束Table 2 Structure and property constraints in MINLP

圖5 lgk量化計算結(jié)果與回歸模型預測結(jié)果的比較Fig.5 Comparison of regression model and DFT results for lgk

圖6 反應(yīng)溶劑設(shè)計結(jié)果（×代表80種初始選取溶劑，代表設(shè)計結(jié)果）Fig.6 Reaction solvents design results(×are 80 initial solvents,are design results)

表3 反應(yīng)溶劑設(shè)計結(jié)果Table 3 Reaction solvents design results

3 結(jié) 論

本文將QM-CAMD 方法進行集成，同步考慮了反應(yīng)速率常數(shù)和選擇性兩個目標函數(shù)的反應(yīng)溶劑設(shè)計，并以藥物合成過程中2,4-二氯-5-硝基嘧啶與對氨基苯腈發(fā)生的SNAr 反應(yīng)為例進行了反應(yīng)溶劑設(shè)計。首先，通過QM 計算的反應(yīng)速率常數(shù)數(shù)據(jù)與COSMO-SAC 等模型得到溶劑描述符進行回歸，得到了該反應(yīng)的反應(yīng)動力學模型，隨后以反應(yīng)速率常數(shù)與選擇性最大化為目標函數(shù)通過CAMD 方法進行反應(yīng)溶劑設(shè)計，最終對溶劑進行評價排名，得到了速率與選擇性更優(yōu)的溶劑。相較于文獻報道中采用的THF 溶劑，其主反應(yīng)的反應(yīng)速率常數(shù)提升了62.8%，副反應(yīng)的反應(yīng)速率常數(shù)降低了85.1%，反應(yīng)選擇性提高了一個數(shù)量級。在未來的工作中，將對本文的理論設(shè)計結(jié)果進行進一步實驗驗證。

符 號 說 明

f(ni)——基團結(jié)構(gòu)約束

G——52種預選基團集合

ΔGp→c,IG——Gibbs自由能變，kJ·mol-1

ΔGsolv——分子的溶劑化自由能，kJ·mol-1

ΔG?——氣相活化Gibbs自由能，kJ·mol-1

ΔG?,L——液相活化Gibbs自由能，kJ·mol-1

h——普朗克常數(shù)

k——反應(yīng)速率常數(shù)，L·mol-1·s-1

kB——Boltzmann常數(shù)

ni——基團結(jié)構(gòu)

PLh——溶劑性質(zhì)h約束下限

PUh——溶劑性質(zhì)h約束上限

Ph(ni)——基團ni對溶劑性質(zhì)h的貢獻值

R——普適氣體常數(shù)

r2——決定系數(shù)

T——溫度，K

κ——傳輸系數(shù)