高精度终点自由能方法的开发与测试。4. 扩展溶剂适用范围及logBB预测能力

《The Journal of Physical Chemistry B》：Development and Test of Highly Accurate End Point Free Energy Methods. 4. Expanding Solvents Capability and logBB Prediction

　　高分辨率图像下载 MS PowerPoint 幻灯片

我们报告了一个扩展的泊松-玻尔兹曼表面积（PBSA）溶剂化模型，用于高通量计算多种有机溶剂中的溶剂化自由能（SFEs）和分配系数（logP）。此外，还开发了一个预测血脑屏障（BBB）渗透性（logBB）的模型。利用明尼苏达溶剂化（MNSol）数据库中的1246个实验SFE测量值，我们为27种有机溶剂参数化了溶剂可及表面积（SASA）的非极性项，从而能够在广泛的溶剂环境中进行PBSA计算。在MNSol数据集中，PBSA的整体均方根误差（RMSE）为1.10 kcal/mol，而基于B3LYP/6-31G*水平的量子力学SMD计算的RMSE为0.98 kcal/mol。为了评估参数的转移性，我们计算了248种化合物在五种溶剂系统中的水-有机溶剂logP值；PBSA的整体RMSE为2.02 log单位，略高于SMD的1.5 log单位。最后，我们基于PBSA得出的水合自由能和选定的有机溶剂SFE开发了一个预测BBB渗透性（logBB）的定量模型。该模型在训练集（N = 865）和测试集（N = 97）上均表现出稳定的预测性能，RMSE约为0.65 log单位，MAE约为0.5 log单位，并且在测试集上显示出改进的排序能力（Kendall’s tau = 0.45）。总体而言，这项工作将PBSA的适用范围扩展到了广泛的溶剂，并建立了物理上有根据且计算效率高的方法来预测一些关键的ADMET描述符。

特刊

作为《物理化学B杂志》“分子模拟与计算化学：Peter A. Kollman的遗产”特刊的一部分发表。

引言

在电子结构计算、(1)生物分子静电学、(2,3)和计算机辅助药物设计中，模拟溶剂效应是必不可少的。显式溶剂模型可以精确模拟热力学系综并分析多体相互作用。然而，它们的广泛应用受到高计算成本的限制。相比之下，隐式溶剂模型通过将溶剂视为可极化的介电连续体来近似溶质-溶剂静电相互作用。当与适当的非极性模型结合使用时，隐式溶剂模型为计算溶剂化自由能（SFE）提供了一个高效的框架。(4) SFE是药物发现中的一个关键物理化学描述符。药物的水溶性直接影响制剂设计和生物利用度，而分配系数（logP）是评估药物的吸收、分布、代谢、排泄和毒性（ADMET）特性的核心指标。(5) 值得注意的是，从溶剂j到i的溶质logP可以通过SFE计算得出：logP??/??=?Δ????→??????ln10logPi/j=?ΔGj→iRTln?10(1)Δ????→??=Δ?????Δ????ΔGj→i=ΔGi?ΔGj(2) 已经开发了多种理论模型用于计算SFE。基于密度的溶剂化模型（SMD）是一种代表性的隐式溶剂模型，它将可极化连续体模型（PCM）与精细调整的非极性项结合起来。(6) 由于其设计，一旦有了七个特定于溶剂的参数，SMD就可以扩展到任意溶剂。与之密切相关的模型，包括COSMO-RS、(7,8)嵌入式簇参考相互作用位点模型（EC-RISM），以及Miertus–Scrocco–Tomasi（MST）模型(10)，也在预测SFE方面表现出高准确性，这从它们在SAMPL系列盲测挑战中的出色表现中得到了证明。(11?19) 另一种广泛采用的用于模拟反应场静电学的框架基于泊松-玻尔兹曼（PB）方程。PB模型在准确性和计算效率之间提供了良好的平衡，因此特别适用于通常需要密集计算和处理大分子的药物发现。(20) 广义玻恩（GB）模型进一步提高了计算效率，使得在分子动力学（MD）模拟中可以频繁进行计算，尽管其公式通常不如PB严格。(21) 因此，人们投入了大量努力来发展和改进基于PB的溶剂化模型，包括在数值PB求解器(22?25)和计算协议方面的进步。例如，Swanson及其同事提出了两组针对突变和光滑介电边界的优化原子半径(26)，Venken等人提出了一个用于结合亲和力计算的异质介电常数方案(27)，并总结了用于蛋白质-配体结合亲和力计算的PB模型的系统研究(20)。最近，我们通过拟合实验水合自由能（HFE)(28)开发了一组基于GAFF2原子类型的半径，并将其应用于logP和logD计算，证明了出色的准确性(29,30)。深度神经网络的发展使得模拟溶剂效应成为可能。图神经网络（GNN）已被用于从显式或隐式溶剂计算的数据中学习溶剂势的平均力（PMF），Chen等人以及Airas等人基于SchNet架构开发了用于MD模拟的隐式溶剂模型。(31,32) 更近期，Wu等人引入了一个GNN框架PBGNN来预测PB能量。(33) 在实际的药物设计中，预测溶解度和分配系数需要跨多种溶剂进行高通量计算。从头算方法如SMD和COSMO-RS通常计算要求过高，而基于神经网络的溶剂模型往往受到有限的训练数据的限制，并且在实现广泛的溶剂转移性方面面临挑战。在这项工作中，我们为27种有机溶剂提供了泊松-玻尔兹曼表面积（PBSA）模型的参数集，从而能够高效计算不同有机溶剂中的SFE和分配系数。此外，我们还利用PBSA计算的SFE通过多元线性回归开发了一个预测血脑屏障（BBB）渗透性（logBB）的模型，假设logBB在多种溶剂中对SFE具有线性依赖性。LogBB是中枢神经系统（CNS）药物发现的关键指标，因为它反映了化合物到达人体CNS目标位点的能力。(34) 最近，logBB预测模型主要使用QSAR和机器学习方法进行分类或量化开发。(35?42) 然而，数据可用性的限制意味着定量logBB预测通常依赖于更简单的QSAR方法或经典机器学习算法，这些方法在小数据集上通常更加稳健。(42,43)

方法

数据准备和处理

PBSA非极性参数是通过拟合明尼苏达溶剂化（MNSol）数据库（v2012）中的实验SFE数据来开发的。(44) 共使用了1246个SFE数据点，涵盖27种溶剂进行参数化。PBSA计算的溶质几何结构取自MNSol XYZ文件，并使用Schr?dinger Maestro（v11.2）转换为MOL2格式。(45) 使用AmberTools中的Antechamber（49）模块指定了GAFF2（46）原子类型和ABCG2（47,48）原子电荷。(50) 从Leo等人(51)收集了一个独立的logP数据集来评估结果的PBSA溶剂化模型。Leo及其同事数据集中带有脚注的logP条目被排除，因为这些值是在特定或非标准实验条件下测量的。此外，报告了多个值且差异大于0.1 log单位的条目也被排除。对于有多个值的条目，如果它们仍在这一范围内，则采用中位数值作为代表logP。对于这部分数据，从PubChem（https://pubchem.ncbi.nlm.nih.gov/）检索了2D结构信息，并使用Open Babel 3.1.0（52）通过“--gen3d”选项生成了3D几何结构，并保存为MOL2格式。每个分子只使用了一种构象。原子类型和部分电荷的分配遵循了用于MNSol化合物的相同协议。

对于logBB预测模型的构建，最初从Shaker等人(43)处获得了1000个具有定量测量logBB值的化合物。值得注意的是，这个数据集中的实验方法多种多样，从体内技术如放射性标记到体外方法如固定化人工膜测定不等。虽然体内实验通常提供最准确的logBB值，但它们的技术复杂性和成本限制了数据的可用性。为了克服这一限制并提高计算模型的泛化能力，将体外logBB数据纳入了建模过程。移除了分子量>600、盐类和多组分（多含量）的记录，得到了963个经过整理的logBB数据点。随机选取了97个分子作为测试集，其余866个分子用于训练。训练集和测试集分子的分子量和logBB值的分布分别绘制在图1中。logBB数据集的分子结构准备、原子类型和部分电荷分配遵循了与logP数据集相同的程序。

图1

图1. 训练集和测试集中的分子量（Da）和logBB值分布。

高分辨率图像下载 MS PowerPoint 幻灯片

PBSA计算和参数化

PB能量使用Delphi V4版本1.1进行计算。(53) 使用了Sun等人的精细调整后的原子半径(28)，介电边界定义为具有1.4 ?探测半径的库仑边界。溶剂可及表面积（SASA）是使用内部的MS程序计算的。PB计算的额外设置与我们之前的出版物相同。(30)

SFE计算如下：??solv=Δ??PB+Δ??SASA(3)Δ??SASA=??SASA+??ΔGSASA=γSASA+b(4)其中γ和b是特定于溶剂的参数。在这个PBSA模型中，每种溶剂由三个量参数化：介电常数（ε）、γ和b。对于水，γ = 0.005，b = 1.12，来自Sun等人(28)。使用公式1计算logP。

进行系统搜索以确定每种溶剂的最佳γ和b值，同时保持介电常数不变。

从头算计算溶剂化自由能

使用SMD（6）隐式溶剂模型进行了SFE和测试集logP的基准计算。所有计算都使用Gaussian 16（54）在B3LYP/6-31G*理论水平上进行。对于每种溶质，几何结构在同一水平上使用SMD进行优化。然后通过SMD能量与使用相同几何结构的气相单点能量之差获得SFE。

用化学方法计算溶剂化自由能

Another基准计算使用热力学积分（TI）对测试集logP值进行了处理。在AMBER 2022中实现的化学增强采样（ACES）方法用于溶剂化自由能计算，而GPU-TI用于水合自由能计算。(55) 使用AmberTools 2022中的TLEAP模块构建了溶剂-溶剂立方箱，最终尺寸约为30 ?。模拟协议与我们之前的工作相同。(30) 进行了CPU-TI以平衡系统，并选择了九个λ窗口（0.01592, 0.08198, 0.19331, 0.33787, 0.5, 0.66213, 0.80669, 0.91802, 和0.98408）用于GPU-TI计算，初始配置从CPU-TI的λ = 0.01592采样。GPU-TI用于ACES模拟的配置采样。TI模拟（CPU-TI, GPU-TI和ACES）在周期性边界条件和等温-等压NPT系综下进行。系统温度通过Langevin动力学调节，碰撞频率为2.0 ps–1。使用蒙特卡洛压力计进行压力控制，并为溶质打开了SHAKE约束。GPU-TI在λ = 0.5时的最后一个快照作为ACES的初始配置，溶剂箱尺寸扩大到约46 ?。对于溶液相和气相，分别使用了九个λ窗口（0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 和0.9）。每种溶质总共进行了2 ns的模拟，其中前0.5 ns的模拟在自由能计算期间被舍弃。

结果和讨论

PB溶剂化模型的开发

我们首先报告了MNSol数据库中可用溶剂的拟合PBSA非极性参数以及PBSA和基准SMD计算的误差评估（表1）。在整个数据集中，PBSA对SFE的整体RMSE为1.10 kcal/mol，而SMD基准的RMSE为0.98 kcal/mol。每种溶剂的PBSA计算SFE值提供在支持信息中。总体而言，PBSA在SFE方面的准确性略低于SMD，但所需的时间显著减少。值得注意的是，这里使用的实验SFE值不仅用于PBSA参数化，还用于SMD模型的参数化。在我们之前的工作中，我们研究了PBSA准确性的构象依赖性，并表明使用高质量的构象集合可以提高与实验的一致性。然而，为了保持参数开发和高通量任务应用之间的一致性，我们使用每种分子的单一构象进行参数拟合和所有后续的PBSA计算。这种策略还避免了耗时的构象采样步骤，并提高了计算吞吐量。值得注意的是，在开发ABCG2电荷模型期间，也对MNSol数据集进行了TI计算，取得了更好的准确性，SFE计算的RMSE为0.89 kcal/mol。表1显示了使用PBSA和SMD计算得到的溶解自由能与实验测量值的误差指标。表1. 使用PBSA和SMD计算得到的溶解自由能与实验测量值的误差指标

| 溶剂 | PBSA | B3LYP/6-31G* | MSE | MUE | MSE |
| -------------- | ------------------ | ---------------- | -------------- | ---------------- |
| 正戊烷 | 26 | 1.83 | 1–0.01 | 90 | 3.40 | 0.03 |
| | | | 0.44 | 0.64 | –0.09 |
| 正己烷 | 59 | 1.88 | 19–0.01 | 90 | 2.20 | 0.30 |
| | | | 0.89 | 1.11 | 0.08 |
| 正庚烷 | 69 | 1.91 | 13–0.01 | 90 | 2.00 | 0.53 |
| | | | 1.04 | 1.40 | 0.33 |
| 2,4-三甲基戊烷 | 32 | 1.93 | 58–0.02 | 94 | 4.20 | 0.11 |
| | | | 0.53 | 0.63 | –1.68 |
| | | | 1.68 | 1.73 |
| 正辛烷 | 38 | 1.94 | 0.60 | 92 | 4.50 | 0.10 |
| | | | 0.43 | 0.60 | 0.20 |
| | | | 0.40 | 0.50 |
| 正壬烷 | 26 | 1.96 | 0.50 | 90 | 0.35 |
| | | | 0.02 | 95 | 85 | 0.01 |
| | | | 0.34 | 0.51 |
| | | | 0.21 | 0.35 |
| 正癸烷 | 39 | 1.98 | 0.66 | 94 | 4.50 | 0.15 |
| | | | 0.89 | 1.04 | 0.68 |
| | | | 0.02 | 97 | 4.50 |
| | | | 0.92 | 1.40 | 0.33 |
| | | | 0.54 | 0.89 |
| 环己烷 | 83 | 2.01 | 65–0.02 | 97 | 6.00 | 0.47 |
| | | | 0.88 | 1.11 | 0.64 |
| | | | 0.63 | 0.81 |
| 正十六烷 | 19 | 2.04 | 0.02 | 92 | 13.70 |
| | | | 0.25 | 0.95 | 1.23 |
| | | | 0.41 | 0.68 | 0.97 |
| 辛烷混合物 | 27 | 2.19 | 0.60 | 91 | 1.80 |
| | | | 0.42 | 0.85 | 1.06 |
| | | | 0.65 | 0.73 |
| | | | 0.89 | |
| 四氯化碳 | 78 | 2.22 | 0.01 | 85 | 2.80 | 0.29 |
| | | | 1.05 | 0.13 |
| | | | 0.49 | 0.68 |
| 苯 | 74 | 2.27 | 0.06 | 74 | 2.30 | 0.16 |
| | | | 0.80 | 1.07 |
| | | | 0.88 | 1.19 |
| 苯二甲醚 | 47 | 2.38 | 0.02 | 26 | 94 | 0.12 |
| | | | 0.51 | 0.74 |
| | | | 0.62 | 0.69 |
| | | | 0.81 |
| 乙苯 | 28 | 2.43 | 0.02 | 57 | 35 |
| | | | 0.11 | 0.50 |
| | | | 0.77 |
| | | | 0.41 | 0.49 |
| | | | 0.63 |
| 二异丙基醚 | 21 | 3.38 | 0.00 | 60 | 0.07 |
| | | | 0.81 | 1.00 |
| | | | 0.52 | 0.82 |
| | | | 1.14 |
| 二乙醚 | 71 | 4.24 | 0.00 | 0.03 | 0.60 |
| | | | 1.01 | 1.29 |
| | | | 0.39 |
| | | | 0.86 |
| | | | 1.23 |
| 碘苯 | 20 | 4.54 | 0.02 | 79 | 66 |
| | | | 0.01 | 0.62 |
| | | | 0.76 | –0.04 |
| | | | 0.32 |
| | | | 0.44 |
| 氯仿 | 10 | 74.71 | 13–0.02 | 258 | 4.70 |
| | | | 0.30 | 1.02 |
| | | | 1.40 |
| | | | 0.45 |
| | | | 0.89 |
| 丁酸乙酯 | 21 | 4.99 | 0.02 | 31 | 5.00 |
| | | | 0.27 |
| | | | 1.39 | –1.08 |
| | | | 1.31 |
| | | | 1.60 |
| | | | 0.44 |
| 溴苯 | 27 | 5.39 | 0.02 | 67 | 6.60 |
| | | | 0.09 |
| | | | 0.68 |
| | | | 0.98 | –0.28 |
| | | | 0.41 |
| | | | 0.50 |
| 氯苯 | 38 | 5.69 | 0.02 | 87 | 4.70 |
| | | | 0.09 |
| | | | 0.57 |
| | | | 0.80 | –0.40 |
| | | | 0.50 |
| | | | 0.61 |
| 乙酸乙酯 | 23 | 5.98 | 0.02 | 57 | 5.6 |
| | | | 0.03 |
| | | | 0.85 |
| | | | 1.14 |
| –1-二氯乙烷 | 39 | 10.12 | 0.02 | 37 |
| | | | 0.02 |
| | | | 0.58 |
| | | | 0.85 |
| | | | –0.05 |
| | | | 0.58 |
| | | | 0.87 |
| 1-戊醇 | 22 | 15.13 | 0.00 | 24 |
| | | | 0.03 |
| | | | 0.74 |
| | | | 0.27 |
| | | | 0.94 |
| 1-丁醇 | 21 | 17.33 | 0.02 | 29 |
| | | | 0.83 |
| | | | 0.15 |
| | | | 0.83 |
| | | | 0.05 |
| | | | 0.98 |
| 乙腈 | 73 | 5.68 | 0.02 | 80 |
| | | | 0.02 |
| | | | 0.87 |
| | | | 0.02 |
| | | | 0.49 |
| 二甲基亚砜 | 74 | 6.82 | 0.02 | 48 | 6.60 |
| | | | 0.04 |
| | | | 0.93 |
| | | | 1.03 |
| | | | 0.81 |
| | | | 1.07 |
| 总计 | 124 | 0.13 | 1.10 | 0.09 | 0.71 |

a N代表与溶剂相关的溶质数量。ε是介电常数。λ和b是特定于溶剂的参数。MSE：带符号的平均误差。MUE：无符号的平均误差。RMSE：均方根误差。MSE、MUE和RMSE的单位是kcal/mol。

表2总结了使用PBSA、SMD和TI计算得到的logP与实验测量值的RMSE、MAE和MUE。所有溶质的计算logP值都提供在支持信息中。在涵盖五种水-有机溶剂系统的249种溶质中，PBSA的总体RMSE为2.02 log单位，而SMD和TI的RMSE分别为1.55和1.59 log单位。总体而言，TI和SMD在传递自由能（以及logP）的准确性上表现相当，而PBSA的RMSE大约大0.5 log单位。

在溶剂系统层面，PBSA的表现优于SMD，特别是在水-苯（PBSA的RMSE = 1.59，SMD的RMSE = 1.90）和水-四氯化碳（PBSA的RMSE = 1.69，SMD的RMSE = 1.76）系统中。PBSA在水-丁醇系统中的异常值较大，其RMSE为3.10。排除这四种溶质后，PBSA的RMSE降至2.46，SMD为1.08，TI为1.57。

注意，logP预测的RMSE值分别为1-丁醇的3.10 log单位和1-戊醇的1.81 log单位。在我们之前的工作中，我们报告了水/1-辛醇系统中一大组类似药物分子的logP预测的RMSE为0.91 log单位。这些结果表明，我们的PBPK模型对更疏水的溶剂表现更好，这与疏水性排名1-辛醇 > 1-戊醇 > 1-丁醇一致。有趣的是，使用PB极性能量和SMD非极性能量的组合计算的logP值显示了显著的预测准确性提高（1-丁醇的RMSE = 1.52 log单位，1-戊醇的RMSE = 0.92 log单位），这表明用于估计溶解自由能非极性组分的SASA缩放方法需要针对更极性的溶剂进行改进。

接下来，我们构建了一些基于溶解自由能（logBB）的预测模型，通过将线性回归拟合到物理化学描述符（如水/1-辛醇logP、水/环己烷logP、溶质表面积和水合自由能）上来实现。为了保持模型的简洁性并降低过拟合的风险，我们将有机溶剂的数量限制在最多五种。该模型通过不同溶剂的组合来模拟分子在不同环境中的内在转移能力。然后，我们对本研究中参数化的27种溶剂以及我们之前工作中使用的1-辛烷、甲苯和环己烷，共174,436种溶剂组合进行了全面的搜索。（29,30）HFE始终用于多变量线性回归。对于每种组合，我们进行了10次自助法迭代以确定最佳组合。在每次迭代中，训练集被随机分为9:1的比例，较小的子集用于验证。训练集中的一个化合物甲氨蝶呤因异常HFE而被排除在模型拟合之外。训练集上的最佳模型使用水（HFE）、1-辛烷、四氯化碳、正戊烷、正辛烷和二甲苯混合物的SFE作为描述符：

```python
logBB = -0.1116 × HFE + 4.5656 × SFEocl?1 + 245.4075 × SFEcar?1 + 170.7768 × SFEpen + 1506.9692 × SFEoce + 891.4166 × SFEyl?? + 30.4952
```

其中HFE表示溶质在水中的水合自由能，SFE表示溶质在有机溶剂中的溶解自由能。下标表示溶剂：ocl为1-辛烷；car为四氯化碳；pen为正戊烷；oce为正辛烷；xyl为二甲苯混合物。最佳模型是根据10次自助法迭代的平均均方根误差选择的。随后，使用完整的训练集拟合了每个SFE的系数。

logBB预测模型的性能在表3中进行了总结，并在图2中进行了说明。线性模型在训练集和测试集上都显示出一致的准确性。在训练集上，它的RMSE为0.65，MAE为0.49。在测试集上，RMSE仍为0.65，MAE为0.50。此外，Kendall’s Tau在测试集上相对于训练集有所提高，表明排名顺序性能有所改善。水、四氯化碳和戊烷的负系数反映了与过度亲水性或非生产性疏水捕获相关的惩罚，而辛烷、辛烷和二甲苯的正系数表明在类似膜和芳香环境中的有利溶解促进了BBB穿透。总体而言，这些结果表明模型在化学多样化的化合物中具有良好的转移能力，并具有稳定的预测性能。对于分类任务，应用了-1.0的logBB阈值来确定分子是否能够穿过BBB。使用这一标准，模型的准确率为0.84，敏感性为0.98，特异性为0.07。考虑到跨膜转运机制的复杂性，评估我们的模型在涉及主动转运的化合物上的性能非常重要。我们通过查询ChEMBL数据集识别出38种P-糖蛋白（P-gp）底物，其中7种包含在