DRUGAI
靶标-药物的解离速率常数(koff)对药效与给药频率影响重大。本文提出了一种耦合深度学习与数学规划法的小分子药物从头设计框架，可高效优化设计koff较小的慢解离“长效药”。首先，创建了一个包含2,773个koff实验数据的高质量数据库，并基于此，开发了一种基于混合专家架构的新型解离动力学通用模型，可实现高效、精准预测koff的目标。其次，耦合解离动力学模型与数学规划法，优化设计具有低koff的潜在药物候选物。最后，利用随机加速分子动力学模拟方法严格验证所设计的潜在药物候选物。以HSP90候选抑制剂设计为例，本文成功设计出许多全新的潜在HSP90候选抑制剂，与已知性能优异的HSP90抑制剂相比，潜在HSP90候选抑制剂的驻留时间(τ = 1/koff)最多提高了45.7%。上述结果验证了本文所提框架在设计慢解离药物候选物方面的可行性和有效性。
引言
近年来，药物发现领域发生了范式转变，强调了动力学性质的重要性，如靶标-药物的解离速率常数(koff)和驻留时间(τ, τ = 1/koff)。解离动力学性质在确定药效方面发挥着关键作用。例如，一些抑制剂与其靶标表现出微摩尔水平的结合活性，而在细胞培养和动物模型中均表现出纳摩尔水平的抗增殖效力和抗肿瘤效力，这一现象可归结于药物对靶标的持久结合，这也会影响药物的给药频率。因此，在设计和开发小分子药物候选物时，考虑解离动力学性质是非常重要的。
因此，本文开发了一种基于深度学习的新型解离动力学通用模型，该模型采用混合专家(MoE)架构，并基于新建立的解离动力学数据库进行训练，从而实现高效精准预测koff的目标。随后，将所开发的深度学习模型集成于一个基于优化的从头药物设计框架，以反向设计在特定靶标上具有长驻留时间的潜在药物候选物。
方法
本文所提出的慢解离药物从头设计框架如图1所示。首先，开发了基于MoE的解离动力学通用模型，利用易于获取的配体SMILES字符串和蛋白质FASTA字符串，来高效精准预测pkoff(-logkoff)(图1(a))。其次，集成解离动力学模型与数学规划法，将慢解离药物设计问题表述为由目标函数(pkoff)、结构约束、性质约束组合而成的混合整数非线性规划模型(图1(b))。再次，采用分解式算法求解数学规划模型，从而高通量优化设计得到针对某一靶标的按pkoff降序排序的定制虚拟候选药物库(图1(c))。最后，潜在药物候选物通过严格的随机加速分子动力学模拟(τ-RAMD)方法进行验证(图1(d))。
图1. 慢解离药物从头设计框架
基于MoE的解离动力学通用模型的架构如图2所示。
图 2. 基于MoE的解离动力学通用模型的架构
结果与讨论
基于MoE的解离动力学通用模型性能
经过100次随机训练，基于MoE的解离动力学通用模型在训练集、验证集和测试集上的平均MSE和r值分别为0.200±0.012、0.300±0.045、0.288±0.028和0.890±0.007、0.821±0.025、0.835±0.017，测试集的预测结果表明，所建立的解离动力学模型在评估靶标-药物解离动力学方面具有可接受的预测能力，并能够对新样品进行满意的泛化预测。若不使用MoE层，解离动力学通用模型在训练集、验证集和测试集上的平均MSE和r值分别为0.245±0.018、0.320±0.047、0.307±0.034和0.864±0.010、0.808±0.028、0.825±0.022，说明引入MoE层的重要性。图3展示了100个随机训练的基于MoE的解离动力学通用模型在整个数据集上的预测结果。通过比较pkoff预测值(以100次预测的均值±标准偏差表示)和pkoff实验值来评估预测误差，在整个数据集上的R2、MAE和MAPE分别为0.793、0.324和38.4%。pkoff预测值(以100次预测的均值±标准偏差表示)和pkoff实验值之间的绝对预测误差范围为0到2.473，大约96.5%的样本显示绝对预测误差小于或等于1。因此，在使用这100个随机训练的MoE模型时，预测误差较小。
图 3. 100个随机训练的基于MoE的解离动力学通用模型的预测结果
为了评估特征提取层(MoE层之前的层)在提取靶标-药物特征方面的能力，以及MoE层在提高解离动力学模型预测准确性中的作用机制，本文利用提取的靶标-药物特征和t-SNE聚类方法创建了化学空间，如图4所示。图4(a)中的不同颜色代表按蛋白样本大小降序排列的前十种蛋白质。图4(a)显示，几乎所有样本都根据其蛋白质类型得到了良好的聚类，这表明特征提取层通过完全数据驱动的方式成功区分了靶标-药物特征。在图4(b)中，相同的化学空间用专家类别重新标注，标注的颜色通过MoE层中门控网络对每个样本的最大权重确定的。一方面，某些蛋白质(例如蛋白质1和蛋白质7)的预测结果主要受到特定专家(例如专家1)的影响。此外，具有相似靶标-药物特征的蛋白质(例如蛋白质4-6)与同一专家(例如专家1)关联。这些观察结果证实了一些专家特别擅长处理特定蛋白质，尤其是那些具有相似靶标-药物特征的样本(例如蛋白质4-6)。另一方面，一些蛋白质(例如蛋白质8)的预测结果涉及多个专家(例如专家3和专家4)，突显了MoE系统协作机制在处理展现独特靶标-药物特征的特定蛋白质复杂样本时的优越性。以上发现表明，门控网络可根据专家各自擅长处理的靶标-药物特征来分配专家权重，从而更好地预测koff大小。
图 4. 由靶标-药物特征和t-SNE聚类方法创建的化学空间
此外，本文使用100个随机训练的MoE模型对来自Amangeldiuly等人和Liu等人的样本进行pkoff预测（以100次预测的均值±标准偏差表示）。最终，对Amangeldiuly等人和Liu等人样本的预测结果的r值分别为0.870和0.728，表示本文的MoE模型在预测koff方面的高稳健性和高准确性。这一结果可能归结于本文构建了一个经过良好清理且样本多样的高质量数据库。此外，相较于依赖靶标-药物三维结构的深度学习模型，本文模型所利用的文本描述符容易获取，因此计算效率更高。
4.2 数学规划法为HSP90蛋白设计定制虚拟抑制剂库
针对HSP90靶标，本文集成深度学习模型与数学规划法设计慢解离HSP90候选抑制剂，设计流程如图5所示。
图 5. 慢解离HSP90候选抑制剂设计流程
首先，收集具有慢解离动力学的HSP90抑制剂及其koff实验值，并将其视为参考化合物。其次，采用Bemis–Murcko算法和基于骨架的相似性算法提取参考化合物的骨架，并从骨架数据库中识别出120个相似骨架，删除重复骨架后剩下37个骨架。基于38个骨架(包括参考化合物的骨架)和25个常用基团，HSP90候选抑制剂设计任务可表述为由目标函数(pkoff)、药物结构约束和药物性质约束(里宾斯基五规则、合成可及性和深度学习模型预测的结合亲和力性质)组合而成的混合整数非线性规划模型。最后，采用分解式算法求解高度非线性的混合整数非线性规划模型，以生成可行解池(即按pkoff降序排序的HSP90候选抑制剂库)，这些可行解是通过在模型约束下组合骨架和基团生成的。本文成功设计出261,736个SMILES字符串，这些SMILES字符串表示满足里宾斯基五规则、合成可及性和结合亲和力性质且结构合理的化合物。在这261,736个SMILES字符串中，70,360个SMILES字符串的pkoff值高于参考化合物的pkoff(pkoff = 2.598±0.281)，表明本文所提框架在从头设计具有慢解离动力学和满意结合亲和力的药物候选物方面具有强大的能力。整个设计过程在台式计算机(Intel(R) Core (TM) i7-10700F CPU @ 2.90GHz 24.0 GB RAM)上大约耗时5.4小时。通过爬虫，发现261,736个候选分子中只有96个化合物(0.037%)在PubChem数据库中存在，这表明所提框架不仅能够找到现有的药物候选物，而且还有能力设计从未被研究过的新化合物。
此外，采用ECFP描述符和t-SNE方法创建化学空间(图6(a))，用以表征261,736个设计的HSP90候选抑制剂的结构多样性，其中不同颜色代表不同的骨架。如图6(a)所示，大量设计的HSP90候选抑制剂分布在靠近参考化合物(由图6(a)中的五角星表示)的位置，表明所提框架在设计与参考化合物相似的药物候选物方面具有强大的能力。此外，远离参考化合物的众多化合物的存在突显了所提框架在设计显著不同于参考化合物的HSP90候选抑制剂方面的能力。这种能力可以归因于基于骨架的相似性算法，该算法利用24个不同的评估标准来扩大药物候选物的设计范围。
图 6. HSP90候选抑制剂的设计结果。(a)使用ECFP描述符和t-SNE方法设计的HSP90候选抑制剂的化学空间。(b)参考化合物的二维分子结构。(c)潜在HSP90抑制剂的二维分子结构。(d)具有全新骨架的化合物的二维分子结构。
此外，对排名结果的分析显示，PubChem数据库中可用的化合物的pkoff值普遍低于参考化合物的值。这一发现表明，目前已知的化合物不太可能在解离动力学方面表现优越，突显了设计新药候选物的迫切需要。关于其他尚未上市且已知合成路径有限的HSP90候选抑制剂，本文采用严格的τ-RAMD方法对参考化合物(预测pkoff=2.598±0.281，如图6(b)所示)、潜在HSP90候选抑制剂(预测pkoff=2.608±0.284，如图6(c)所示)和具有全新骨架的化合物(预测pkoff=2.480±0.363，如图6(d)所示)进行严格验证。
使用τ-RAMD方法验证HSP90候选抑制剂
τ-RAMD方法可通过增强采样分子动力学模拟来计算化合物在蛋白质中的相对驻留时间(τ)。计算得到的潜在HSP90候选抑制剂的τ值为1.02±0.20 ns(基于六次重复实验的平均值)。相比之下，参考化合物的τ值为0.70±0.14 ns(基于六次重复实验的平均值)。这个结果表明，潜在HSP90候选抑制剂的τ值比参考化合物大，即解离速率较低。这证明了所提框架在设计具有慢解离动力学的药物候选物方面的强大能力。此外，几个具有全新骨架的HSP90候选抑制剂也被发现展现出缓慢的解离动力学。例如，具有全新骨架的化合物(预测pkoff=2.480±0.363，如图6(d)所示)的τ值为0.50±0.16 ns(基于六次重复实验的平均值)。尽管其τ值未超过参考化合物，但其骨架与潜在HSP90候选抑制剂和参考化合物完全不同，表明通过所提框架有设计出具有慢解离动力学的新药物候选物的潜力。
HSP90是一种在预测解离动力学方面具有挑战性的蛋白质，因为其灵活的腺苷三磷酸(ATP)结合位点由不稳定的α-helix3组成。这个α-helix3片段(残基101-122)可以根据所结合的配体采取“helix”(图7(a))或“loop”(图7(b))构象。“helix”构象有一个额外的瞬态亚口袋，位于ATP结合位点旁边。这个亚口袋提供了在“helix结合位点”上的取代基空间。参考化合物、潜在HSP90候选抑制剂和具有全新骨架的化合物均属于“helix”构象。图8(a-b)中红色和黑色圆圈突出显示的配体片段代表了潜在HSP90候选抑制剂与参考化合物之间的差异。配体占据了靠近α-helix3的瞬态亚口袋，如图8(c-d)所示。最近的研究表明，瞬态亚口袋与HSP90抑制剂的解离动力学之间存在相关性。研究发现，增加R1取代基的大小能够有效稳定α-helix3，从而降低结合和解离速率，尤其是具有更大R1取代基的配体通常表现出更慢的结合和解离动力学。显然，潜在HSP90候选抑制剂的R1取代基大于参考化合物，因而提高了τ值。
图 7. HSP90结合位点的两种主要构象
图 8. 潜在HSP90候选抑制剂和参考化合物结合状态的分子动力学模拟结果
结论
本文建立了一个慢解离药物从头设计框架，用于开发“长效药”。创建了一个高质量的解离动力学数据库，包含了不同靶标-药物的2,773个koff实验值。基于该数据库，随机训练了100次基于MoE的解离动力学通用模型，以高效精准预测koff。解离动力学模型在训练、验证和测试集上的MSE和r值分别为0.200±0.012、0.300±0.045、0.288±0.028和0.890±0.007、0.821±0.025、0.835±0.017。随后，将所开发的解离动力学模型集成到基于数学规划法的从头药物设计框架中，成功设计了一个定制的虚拟HSP90抑制剂库，包含261,736个具有满意结合亲和力的药物候选分子，其中70,360个相较于参考化合物具有更慢的解离动力学。通过τ-RAMD方法严格验证发现，潜在HSP90候选抑制剂的模拟τ值(1.02±0.20 ns)比参考化合物(0.70±0.14 ns)高出45.7%，证明了所提框架在设计具有长驻留时间的全新候选药物方面的可行性和有效性。
第一作者为大连理工大学生物工程学院，赵雨靓博士后。
通讯作者为大连理工大学化工学院，刘奇磊副教授。
个人主页：
https://faculty.dlut.edu.cn/liuqilei/zh_CN/index.htm。
参考资料
Zhao, Y., Zhang, L., Du, J., Meng, Q., Zhang, L., Wang, H., Sun, L. and Liu, Q., 2024. Mixture-of-Experts Based Dissociation Kinetic Model for De Novo Design of HSP90 Inhibitors with Prolonged Residence Time. Journal of Chemical Information and Modeling.
