考虑选择性和反应速率的多目标制药反应溶剂设计

赵红庆，刘奇磊，张磊，董亚超，都健

（大连理工大学化工学院，化工系统工程研究所，辽宁大连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

κ——传输系数