考慮選擇性和反應速率的多目標制藥反應溶劑設計

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

（大連理工大學化工學院，化工系統工程研究所，遼寧大連116024）

引 言

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

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

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

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

1 QM-CAMD設計框架

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

圖1 QM-CAMD設計方法Fig.1 QM-CAMD design method

1.1 反應體系的選擇

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

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

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

1.2 反應機理研究及反應速率常數計算

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

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

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

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

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

（1）使用Gaussian 09[25]軟件對反應物分子進行結構優化和振動分析，并基于合理的過渡態初猜結構搜索真實的過渡態結構，然后進行內稟反應坐標（intrinsic reaction coordinate,IRC）[26]驗證確保過渡態結構正確；

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

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

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

1.3 構建反應動力學模型

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

1.4 反應溶劑設計模型

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

（1）目標函數（其中lgk2與lgk4的表達式由1.2 節及1.3節中描述的方法回歸得到）

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

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

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

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

1.5 分解式算法

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

（1）由選擇的初始基團根據結構約束通過數學算法生成大量結構可行的分子結構（數量為N1）；

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

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

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

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

2 結果分析與討論

2.1 反應機理

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

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

2.2 反應動力學模型

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

其中r2=0.934。

圖4 真空中的反應過渡態結構Fig.4 Transition state structure of reaction in vacuum

表1 真空中反應物、過渡態與產物的相對自由能Table 1 Relative Gibbs free energies of reactants，transition states，and products in a vacuum

其中r2=0.855。

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

2.3 CAMD設計結果

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

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

表2 MINLP模型中的結構和性質約束Table 2 Structure and property constraints in MINLP

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

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

表3 反應溶劑設計結果Table 3 Reaction solvents design results

3 結 論

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

符 號 說 明

f(ni)——基團結構約束

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——普朗克常數

k——反應速率常數，L·mol-1·s-1

kB——Boltzmann常數

ni——基團結構

PLh——溶劑性質h約束下限

PUh——溶劑性質h約束上限

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

R——普適氣體常數

r2——決定系數

T——溫度，K

κ——傳輸系數