基于萘的微管蛋白抑制剂的计算机模拟评估：作为潜在的抗乳腺癌药物

《Scientific African》：In Silico Evaluation of Naphthalene-Based Tubulin Inhibitors as Potential Anticancer Agents Targeting Breast Cancer Cells

　　萨吉尔·哈姆扎·阿卜杜拉伊|雷纳托·阿劳霍·达科斯塔|亚当乌·乌扎伊鲁|穆罕默德·图库尔·易卜拉欣|马加吉·拉丹|阿尤巴·阿卜杜拉伊·穆罕默德|巴巴恩吉达·阿卜杜拉伊·萨杰
艾哈迈杜·贝洛大学物理科学学院化学系，扎里亚

摘要
乳腺癌仍然是全球主要的健康问题，女性中的发病率和死亡率都很高。现有化疗药物的局限性，包括耐药性和不良副作用，突显了需要更安全、更有效的替代方案。在这项研究中，使用多元线性回归（MLR）方法开发了一个统计上稳健的二维定量结构-活性关系（2D-QSAR）模型，并通过内部和外部统计指标进行了验证，如：R2 = 0.951、Q2 = 0.921、R2ext = 0.722、cRp2 = 0.713以及适用性域分析。该模型促进了虚拟筛选，并被用于设计七种新型萘基 tubulin 抑制剂。与 tubulin 目标蛋白（PDB ID：1SA0）的分子对接显示，化合物 1d 的预测结合亲和力为 –142.358，化合物 1e 的结合亲和力为 –147.636，高于秋水仙素（–131.453）。通过 300-ns 分子动力学模拟和 MM/GBSA 计算进一步验证了稳定性和结合强度，化合物 1e 显示出有利的 ΔG 为 –27.96 kcal/mol。药代动力学和药效学分析确认其符合 Lipinski 规则，具有良好的口服生物利用度和低毒性以及有希望的药理特性。这些计算结果表明，所设计的萘衍生物可以作为进一步实验研究的潜在骨架，作为乳腺癌治疗的 tubulin 靶向药物。

引言
乳腺癌是全球影响女性最普遍和致命的恶性肿瘤之一，每年报告的新病例数达数百万例 [1,2]。在亚型中，雌激素受体阳性（ER+）乳腺癌（如 MCF-7 细胞系）由于其对抗化疗的响应性，仍然是治疗干预的重点 [3], [4], [5]。尽管治疗策略有所进步，但对传统化疗药物的耐药性和不良副作用仍突显了开发更有效、更具选择性的抗癌药物的需求 [6], [7], [8]。Tubulin 是抗癌药物设计中最成熟和最有吸引力的靶点之一，它是一种对微管组装和有丝分裂纺锤体形成至关重要的结构蛋白。破坏 tubulin 的聚合会导致细胞周期停滞和凋亡，使 tubulin 成为抗癌治疗的关键分子靶点 [9]。特别是 tubulin 上的秋水仙素结合位点因其独特的作用机制及其克服癌细胞多重耐药性的潜力而受到广泛关注 [10]。最近的研究强调了基于萘的化合物（尤其是苯甲酮、查尔酮和吡唑衍生物）作为有效 tubulin 聚合抑制剂的潜力。这些分子通过与其 tubulin 结合位点的相互作用来抑制微管动态，显示出良好的抗癌活性 [11]。

在这项研究中，采用基于配体的药物设计（LBDD）策略来探索和优化萘衍生的 tubulin 抑制剂，目标是提高它们的生物效力和结合亲和力。利用计算技术为新型治疗药物的研发提供了节省时间和经济可行的途径 [12], [13], [14]。结合定量结构-活性关系（QSAR）建模、分子对接和分子动力学（MD）模拟等计算机方法来评估研究分子的效力、选择性和类药物特性。这些技术有助于识别适合后续实验室验证的有希望的先导候选物 [15]。虽然计算药物发现方法提供了许多优势，如快速筛选大量化合物库和能够探索多种化学骨架，但其使用也存在某些局限性，包括计算机预测与实验结果之间可能存在的不一致性、对目标蛋白准确结构数据的依赖性以及计算资源的需求 [16], [17]。然而，配体-靶标相互作用的模拟提供了对指导先导优化和推进早期药物发现至关重要的机制洞察 [18]。

材料与方法
本研究使用的数据集来自 Wang 等人报告的实验结果 [19], [20], [21]，其中包括评估了 tubulin 抑制活性的化合物。这些活性最初以微摩尔（μM）浓度的 IC?? 值表示。为了提高统计可靠性并确保建模的线性关系，IC?? 值被转换为其负对数形式（pIC??）[22]。分子的二维（2D）化学结构最初使用 Perkin-Elmer ChemDraw 绘制（补充文件表 S1），然后使用 Spartan v14.0 软件套件转换为三维（3D）分子模型。随后，每个 3D 结构在 Spartan 环境中通过密度泛函理论（DFT）在 B3LYP/6-31G* 理论水平下进行几何优化，以获得适合进一步计算分析的能量稳定构象 [23], [24], [25]。

生成分子描述符
在优化数据集的分子几何结构后，使用 PaDEL-Descriptor 软件计算与生物活性相关的一组分子描述符 [26]。为了提高模型的质量和预测能力，通过手动预处理仔细过滤这些描述符，以去除冗余和高度相关的变量。使用 Kennard-Stone 算法通过专用数据预处理工具进一步细化数据集，确保样本的代表性及均匀分布 [27]。最终的数据集被划分为包含 32 个样本的训练集和包含 13 个样本的测试集，以便进行模型构建和外部验证。

提取最佳描述符
在这项研究中，使用 Material Studio v8.0 软件实施遗传算法（GA）从初始描述符池中选择最具信息量的描述符，用于 QSAR 模型的开发 [28]。用于优化的目标函数是留一法交叉验证相关系数（cvR2）[29]。GA 系统地搜索最大化模型预测性能的描述符组合（通过较高的 cvR2 值表示），同时尽量减少使用的描述符数量，从而提高模型的简单性和泛化能力 [30]。

模型开发
在 Material Studio v8.0 套件中应用多元线性回归（MLR）框架生成 QSAR 模型，建立因变量（pIC??）与选定的分子描述符之间的直接线性相关性 [31]。所得回归方程的一般形式为：
(1)
pIC50 = β0 + β1X1 + β2X2 + … + βnXn
其中，pIC50 是因变量，β 是相应自变量 X 的回归系数，β0 是截距。

QSAR 模型质量评估与验证
QSAR 模型的验证在确定其预测可靠性方面起着关键作用，这对于准确计算和后续筛选研究至关重要。留一法交叉验证（LOOCV）是一种稳健的内部验证方法，其中模型在所有数据点上迭代训练，除了一个数据点，通过排除该点来评估预测准确性 [32]。这个循环重复进行，直到每个数据点都被测试，从而全面衡量模型的泛化能力。LOOCV 的迭代性质确保每个数据点都被系统地排除并用作测试案例，从而能够严格和全面地评估模型的预测准确性 [33]。该方法通过将每个观测值都视为潜在的验证实例来最小化偏差。采用了一整套验证指标来评估模型质量，包括 Friedman Lack-of-Fit（LOF）、决定系数（R2）、调整后的 R2（adjR2）、交叉验证 R2（cvR2）以及回归显著性 [34]。还分析了其他统计指标，如回归显著性 F 值、95% 置信水平下的临界 SOR F 值、重复点、实验误差、缺乏拟合的数据点以及非显著 LOF（95%）所需的最小实验误差。除了内部验证外，还通过独立测试集进行模型随机化以进一步评估模型的稳健性。这一过程涉及评估随机数据的平均相关系数（rR）、随机数据的平均决定系数（rR2）、随机数据的平均留一法交叉验证决定系数（rQ2）以及 Y-随机化系数（cR2p）[35]。此外，还使用独立测试集进行了外部验证。外部评估的关键指标包括外部验证系数（extR2）、观察值与预测值的斜率（k）以及预测值与观察值的斜率（k’）等 [36]。

适用性域评估
通过使用 William 图来评估 QSAR 模型的适用性域是一种系统的方法，用于确定对新分子的预测可靠性。该过程从计算训练集和测试集的分子描述符开始，然后使用这些描述符计算每个样本的杠杆值。这些杠杆值构成了 William 图的基础 [37]。一个警告杠杆阈值（h*）有助于识别位于模型适用性域内的化合物。杠杆值低于此阈值的化合物被认为能产生可靠的预测 [38]。相反，杠杆值超过定义阈值的化合物被认为是有影响力的，表明它们的预测响应可能不太可靠，需要进一步审查 [39]。这种结构化的评估通过确认模型对于具有与训练集中化合物相似结构和物理化学特征的化合物的预测可靠性，增强了 QSAR 模型的稳健性，从而验证其在未来预测任务中的适用性。

基于配体的新型 tubulin 抑制剂设计
采用基于配体的药物设计策略生成具有潜在 tubulin 抑制活性的新型候选分子。lead化合物的选择完全基于开发的 QSAR 模型的预测。选择过程侧重于识别具有最高预测抑制效力（pIC??）、最小残差误差、位于模型定义的适用性域（AD）内以及符合药理学特性标准（如 Lipinski 规则五项）的化合物 [40]。因此，从测试集中识别出化合物 1 作为 lead 候选物。对该 lead 的结构优化涉及在模板分子上的特定位置进行有针对性的官能团添加和替换。这些修改由模型进化中使用的选定描述符的平均效应值系统指导。

设计的 tubulin 抑制剂的对接模拟
从蛋白质数据库检索 tubulin 蛋白（pdb id：1SA0）的晶体结构，并使用 Discovery Studio 进行分子对接 [41]。为了验证对接协议的可靠性，将共结晶配体秋水仙素重新对接到其天然结合位点 [42]。随后，使用 Molegro Virtual Docker（MVD）软件将设计的 tubulin 靶向化合物对接到最合适的结合口袋中，该口袋位于 XYZ 坐标 125.40、95.22 和 13.13，半径为 15 ? 的定义球形区域内。每个对接模拟进行了最多 1500 次迭代，种群大小固定为 50 [42]。选择显示最佳结合亲和力的配体构象进行结构解释和进一步分析。

最佳复合物的分子动力学模拟
使用 Amber18 软件包进行分子动力学（MD）模拟，以研究最稳定的蛋白质-配体复合物的动态行为。这些模拟用于阐明从对接研究中识别出的两个顶级抑制剂复合物之间的结合相互作用 [43]。配体参数使用 General Amber 力场（GAFF）推导，而蛋白质则使用 ff14SB 力场进行建模 [44,45]。使用在高斯 09 套件中实现的受限电静势（RESP）方法在 HF/6-31G* 理论水平下获得配体的部分原子电荷 [46]。H++ 服务器用于预测生理 pH（7.00）下可电离残基的质子化状态，估计它们的 pKa 值。系统准备使用 tLeap 模块进行，其中每个复合物都溶解在 TIP3P 水分子的立方周期性盒子中，并用钠离子中和。之后，在多个阶段进行能量最小化以消除空间冲突。水分子和反离子的初始能量最小化包括 2000 步的最陡下降步骤，随后是 3000 步的共轭梯度最小化。然后将类似的协议应用于放松氢原子和水分子。整个系统在 200 ps 的分子动力学（MD）过程中逐渐从 0 加热到 300 K，然后在恒容条件下进行 300 ps 的密度平衡，同时对蛋白质原子施加位置约束。随后进行的平衡包括 10 ns 的无约束 MD 模拟，允许蛋白质-配体复合物完全放松。为了将系统温度保持在300 K，使用了具有2 ps?1碰撞频率的Langevin恒温器进行温度调节。应用了10 ?的非键合相互作用截止距离，并使用Particle Mesh Ewald（PME）方法处理长程静电作用。涉及氢原子的键长通过SHAKE算法进行了约束。最终，在300 K下进行了300 ns的生产性分子动力学（MD）模拟。分析了所得的MD轨迹，以评估蛋白质-配体复合物的结构稳定性和动态行为。计算了包括均方根偏差（RMSD）、均方根波动（RMSF）、回转半径（Rg）和溶剂可及表面积（SASA）在内的关键指标，以评估整个模拟过程中自由蛋白质和结合复合物的构象偏差和灵活性[46]。相互作用自由能是基于分子力学中的广义Born表面积（MM/GBSA）技术计算的。总结合自由能、极性和非极性溶剂化能以及势能组分（静电和范德华相互作用）是从MD轨迹中估算出的[47,48]。

药理学研究
利用了两个可访问的在线工具：http://www.swissadme.ch/ 和 https://biosig.lab.uq.edu.au/pkcsm/ 来预测设计的Tubulin抑制剂的药理学和药效学特性，同时使用PreADMET服务器来预测分子的毒性特性[49,50]。

结果与讨论
通过结合理论和体外方法，开发了一个稳健的2D-QSAR模型。数据集使用Kennard-Stone算法被划分为训练集（n = 32）和测试集（n = 13）。使用Material Studio v8.0进行模型开发，并通过训练集上的多元线性回归（MLR）技术进行内部验证。随后使用独立的测试集进行外部验证，以评估模型的预测可靠性[51]。为了进一步研究分子结构和理化特性对抗Tubulin活性的贡献，应用了遗传函数算法（GFA），提供了驱动生物活性的关键描述符的机制洞察。最终模型结合了自相关和拓扑描述符，包括AATSC1c、ATSC1v、AATS6v、ATS2i、SpMax6_Bhe和SpMin2_Bhp，每个描述符都对观察到的活性有独特的贡献。分析回归系数以确定每个描述符的影响方向和程度。值得注意的是，AATSC1c、ATSC1v、AATS6v和SpMax6_Bhe显示出正系数，表明它们的值增加与预测的抗Tubulin活性增强相关[52]。相比之下，ATS2i和SpMin2_Bhp显示出负系数，表明这些描述符的较低值有利于提高生物活性。此外，计算了方差膨胀因子（VIF）以评估所选描述符之间的潜在多重共线性。所有VIF值（如补充文件表2所示）均低于5的阈值，表明没有显著的相互相关性，从而确认了模型的统计稳健性。VIF超过10的阈值通常表示存在问题性的多重共线性，这可能会影响模型的有效性并需要修订[53]。从开发的QSAR模型中获得的统计指标总结在表1中。内部验证显示Friedman LOF得分为0.007，R2值为0.951，调整后的R2值为0.939，Q2值为0.921，表明模型性能良好。Y-随机化系数（cR2p = 0.713）进一步证实了模型的稳健性，并表明其预测能力并非偶然。外部验证显示预测R2（extR2）为0.722，超过了推荐的阈值[53]。Kaushal等人基于喹诺酮衍生物开发了一个2D-QSAR模型作为潜在的抗癌剂，报告了训练集的相关系数r2=0.78和q2=0.71。外部验证的预测相关系数（pred-r2）为0.68。这些统计参数比本研究中获得的要低，从而突显了我们模型的优越稳健性和预测可靠性[54]。

表1. 开发模型的统计参数。
QSAR模型：pIC50 = 0.000109817 * ATS2i + 0.007430667 * AATS6v + 0.005431168 * ATSC1v - 124.915500296 * AATSC1c + 1.337502926 * SpMax6_Bhe - 2.944015958 * SpMin2_Bhp +3.922015155
/ 无
验证指标
可接受值
模型值
备注
1. Friedman LOF
流动值
0.007
通过
2. R-squared
≥ 0.600
0.951
通过
3. 调整后的R-squared
≥ 0.600
0.939
通过
4. 交叉验证的R-squared
≥ 0.600
0.921
通过
5. 回归显著性F值
80.254
通过
6. 临界SOR F值（95%）
≥ 2.092.503
通过
7. 复制点
0.00
通过
8. 计算的实验误差
0.00
通过
9. 无显著LOF的最低实验误差（95%）
0.031
通过
10. 外部验证
12. 外部验证系数（ext R2）
≥ 0.600
0.722
通过
13. 观察到的活性与预测的活性值在零截距处的图表斜率（k）
0.850.991
通过
14. 预测的活性与观察到的活性在零截距处的图表斜率（k′）
0.851.009
通过
15. r02-r0’2/< 0.300
0.045
通过
16. r2?r0’2
r2<0.100
0.011
通过
Y-随机化参数
17. 平均r<0.500
0.364
通过
18. 平均r2<0.500
0.139
通过
19. 平均Q2<0.500
-0.378
通过
20. cRp2≥ 0.600
0.713
通过
关键：r2是实际与预测活性图表的相关系数，r20是实际与零截距处预测活性值图表的相关系数，r’20是预测与实际活动图表的相关系数，且截距为零。

通过训练集和测试集样本的观察到的pIC??值与预测值的散点图进一步评估了所提出模型的预测性能，如图1所示。数据点的分布显示出预测活性和观察活性之间的强相关性，大多数分子紧密地沿着回归线排列。这表明模型具有高水平的准确性和内部一致性[54]。值得注意的是，训练样本与回归线的偏差最小，表明模型成功捕获了训练数据集中的潜在结构-活性关系。围绕线的紧密聚集反映了强大的内部预测能力，这与模型验证期间获得的高决定系数（R2 = 0.951）和调整后的R2（adjR2 = 0.939）值一致。测试样本也与回归线相当吻合，从而确认了模型的外部预测能力。尽管一些测试分子显示出轻微的偏差，但它们的分布并没有表现出任何系统性偏差或显著的过拟合。这一观察结果与稳健的留一法交叉验证系数（Q2 = 0.921）和令人满意的外部预测平方相关系数（extR2 = 0.722）一致，后者超过了QSAR模型通常接受的阈值[55]。

下载：下载高分辨率图像（1MB）
下载：下载全尺寸图像
图1. 实验pIC50与预测pIC50的图表。
图2显示了残差与实验pIC??值的图表。这种分析作为评估预测误差均匀性和模型中任何潜在偏见的诊断工具[56]。如图所示，训练集和测试集的残差随机分布在零线周围，没有明显的系统趋势或模式。这种随机性表明模型没有受到显著的结构偏差，且预测误差的方差在整个pIC??范围内保持恒定。大多数残差值落在±0.4的狭窄范围内，表明预测活性和实验活性之间有很好的一致性。只有一个测试化合物的残差略高，这种偏差可能归因于其内在的结构多样性或生物数据的实验变异性。尽管如此，这种残差的幅度相对较低，其分布并不影响模型的整体预测性能[56]。

下载：下载高分辨率图像（1MB）
下载：下载全尺寸图像
图2. 实验值与残差的图表。

适用性领域
William’s图用于界定QSAR模型的适用性领域（AD），定义模型预测被认为是可靠和有效的化学空间[58]。该图评估两个关键参数：杠杆值（h）及其对应的阈值（h*），后者表示化合物对模型预测产生显著影响的临界值。图3展示了所开发模型的William’s图，提供了对有影响力的化合物和潜在异常值的视觉评估。所有数据点都落在可接受的标准化残差范围±3内，表明没有响应异常值。然而，根据杠杆标准，一个训练化合物和五个测试化合物被识别为有影响力的化合物，它们的杠杆值超过了阈值（h* = 0.656）。这些化合物虽然不一定错误，但可能会对模型性能产生显著影响，在未来的应用中应谨慎解释。尽管如此，大多数化合物都位于AD范围内，确认所开发的模型具有稳健且明确的化学空间，因此可以期望在未来设计工作中为结构相关的类似物提供可靠的预测[58]。

下载：下载高分辨率图像（2MB）
下载：下载全尺寸图像
图3. 所提模型的William’s图。

描述符平均效应评估
平均效应（ME）指标定量反映了每个描述符对建模化合物的生物活性（pIC??）的平均影响，从而表明它们的相对重要性[59]。图4和补充文件表2展示了模型描述符及其相应的ME值。SpMin2_Bhp显示出最高的负绝对ME值（-4.473），表明它与分子的抗Tubulin活性有最显著的负相关关系。这意味着其较高的值可能会降低分子的效力，可能是由于该描述符所代表的电子或空间因素。另一方面，SpMax6_Bhe具有最大的正ME值（+3.124），表明其对抗Tubulin活性有很强的直接贡献；增加其值可能会提高生物活性。描述符如AATS6v（+1.026）和ATS2i（+0.894）也对模型有积极的贡献，尽管程度小于SpMax6_Bhe。它们的正ME值支持了它们在提高预测性能方面的相关性。描述符AATSC1c（+0.256）和ATSC1v（+0.173）的ME值相对较低，表明它们对模型的贡献适度但有意义。这种平衡的ME值分布反映了模型整合了多样化的结构和理化特征[59]。

下载：下载高分辨率图像（212KB）
下载：下载全尺寸图像
图4. 模型描述符的平均效应评估。

表2. 设计的Tubulin抑制剂及其与Tubulin蛋白目标的结合亲和力（MolDock单位）。

基于使用报告模型的虚拟筛选，化合物1作为先导分子出现，其pIC50值高，模型预测准确（残差值低），并且位于定义的模型范围内。由于SpMin2_Bhp和SpMax6_Bhe具有最大的负和正ME值，因此被选为设计目的，如图4所示。SpMin2_Bhp是按疏水性（hp）加权的Burden矩阵的最小特征值。该描述符的高度负ME值表明较低的值有利于活性。为了提高活性，设计了在特定位置减少疏水性的新分子，可能是通过引入极性官能团或减少庞大的非极性基团来实现的。同时，SpMax6_Bhe指的是按电负性（he）加权的Burden矩阵的最大特征值，在路径长度为6的情况下。该描述符的大正ME值表明较高的值可以增强活性。因此，在2-甲氧基苯酚基团周围添加活性基团（如-OH、-NHCH3、-OCH3和-NH2）可以设计出七种（7）种比先导分子和标准药物（Colchicine）活性更好的新型抗Tubulin候选物，如表2和表3所示。

下载：下载高分辨率图像（335KB）
下载：下载全尺寸图像
图5. 设计模板。
表3. 描述符的值及设计分子的预测活性。

基于所报告模型的虚拟筛选，化合物1作为先导分子出现，其pIC50值高，模型预测准确（残差值低），并且在定义的模型范围内。SpMin2_Bhp和SpMax6_Bhe因具有最大的负和正ME值而被选为设计目的，如图4所示。SpMin2_Bhp是按疏水性（hp）加权的Burden矩阵的最小特征值。该描述符的高度负ME值表明较低的值有利于活性。为了提高活性，设计了在特定位置减少疏水性的新分子，可能是通过引入极性官能团或减少庞大的非极性基团来实现的。同时，SpMax6_Bhe指的是按电负性（he）加权的Burden矩阵的最大特征值，在路径长度为6的情况下。该描述符的大正ME值表明较高的值可以增强活性。因此，在2-甲氧基苯酚基团周围添加如-OH、-NHCH3、-OCH3和-NH2等活性基团可以设计出七种（7）种活性比先导分子和标准药物（Colchicine）更好的新型抗Tubulin候选物，如表2和表3所示。

下载：下载高分辨率图像（335KB）
下载：下载全尺寸图像
图5. 设计模板。
表3. 描述符的值及设计分子的预测活性。

空单元
ID ATS2i AATS6v ATSC1v AATSC1c SpMax6_Bhe SpMin2_Bhp
Pred pIC50 110 64 6.213 17 6.645 61.069 -0.005 3.069 1.947 5.793 1a 122 48.000 154.960 67.928 -0.005 83.135 91.970 5.910 1b 134 96.585 153.070 38.617 -0.007 3.147 1.986 5.950 1c 120 37.319 156.002 68.723 -0.008 3.093 1.969 6.089 1d 113 12.012 160.129 189.839 -0.009 3.071 1.972 6.829 1e 124 01.565 154.160 114.113 -0.008 3.104 1.978 6.353 1f 134 96.585 149.693 38.617 -0.007 3.116 1.987 5.879 1g 120 37.319 156.002 68.723 -0.008 3.093 1.969 6.089 Cholcine 100 21.108 145.112 -22.317 -0.009 3.185 1.864 5.861

分子对接模拟
使用MVD软件通过分子对接研究了Tubulin蛋白（PDB ID：1SA0）与化合物（1）、设计的类似物以及参考药物（Colchicine）之间的相互作用，并利用MolDock评分函数来估计结合能。算法验证是分子对接的关键方面，用于确定该方法准确的姿态预测和结合亲和力估计的能力。如图6所示，准备好的Tubulin目标与共结晶配体一起显示，突出了关键活性位点残留物和共结晶配体的叠加姿态。为了确认配子在特定活性位点的结合精度，将Colchicine重新对接到其原始结合口袋中。重新对接的结果是RMSD为1.357?，远低于≤ 2.0 ?的接受阈值，从而验证了对接协议的精度[60]。

下载：下载高分辨率图像（600KB）
下载：下载全尺寸图像
图6. Colchicine在Tubulin目标活性位点的相互作用及其叠加姿态。

设计的候选物与Tubulin蛋白目标之间的对接模拟显示了显著的相互作用潜力，这从它们的结合亲和力（表2，图7）可以得到证明。值得注意的是，主要候选物的结合亲和力得分为-131.874，而设计出的类似物的亲和力范围在-132.669到-147.636之间。设计化合物1d和1e显示出最佳的结合亲和力，分别为-142.358和-147.636，它们与目标活性位点残基的2D和3D残基相互作用分别展示在图7d和7e中。下载：下载高分辨率图像（536KB）下载：下载全尺寸图像下载：下载高分辨率图像（481KB）下载：下载全尺寸图像下载：下载高分辨率图像（528KB）下载：下载全尺寸图像下载：下载高分辨率图像（487KB）下载：下载全尺寸图像图7. 设计化合物的2D残基相互作用：(A) 1a；(B) 1b；(C) 1c；(D) 1d；(E) 1e；(F) 1f；(G) 1g；以及(H) 秋水仙碱与Tubulin目标蛋白的活性位点残基（PDB ID: 1SA0）。化合物1d与Thr145形成了常规的氢键，这被认为是固定配体的关键。与Thr145的相互作用特别显著，因为该残基靠近秋水仙碱结合区域，并且经常参与与抑制剂的氢键结合。此外，还观察到了配体与Asp98、Gly144和Ile171残基之间的碳-氢键，进一步稳定了配体在结合口袋内的位置。Tyr224（一种已知能支持π电子云相互作用的芳香残基）表现出π-π堆叠，增强了活性位点内的范德华力。Ala12和Ile171之间的π-烷基相互作用也很明显，表明配体1g能够很好地适应Tubulin结合位点的疏水环境。另一方面，1e与几个关键氨基酸残基（尤其是Gln15和Gln11）形成了常规的氢键，表明在结合裂隙内存在强大的极性相互作用网络。特别是Gln15与配体的羟基和胺基形成了双重氢键，表明其在配体固定中起着关键作用。还观察到了与Asn101、Glu183、Gly143和Asn206之间的碳-氢键，这进一步支持了整体的结合亲和力。Gly144形成了π-供体氢键，通过静电作用进一步加强了配体的稳定性。与Tyr224之间的显著π-π堆叠相互作用也表明配体1g能够很好地适应Tubulin结合位点。秋水仙碱作为参考药物，其对Tubulin目标的结合亲和力为-131.453。配体-蛋白质相互作用分析显示，它与Tyr224和Ser140形成了两个常规的氢键，与Ser178和Ile171形成了碳-氢键，并且与Ile171形成了疏水烷基相互作用，如图7g所示。与设计化合物相比，秋水仙碱在结合口袋内建立的稳定相互作用较少，这与它较低的结合亲和力一致。这一观察结果突显了设计配体的增强相互作用特性，表明它们具有更高的结合稳定性和作为有前景的Tubulin抑制剂的潜力，从而支持了它们作为潜在药物候选物的新颖性。

分子动力学（MD）模拟旨在评估Tubulin在其apo形式下以及与化合物1d和1e结合时的稳定性和动态行为，模拟时间为300纳秒。分析了关键动态参数，如RMSD、RMSF、回转半径（Rg）和溶剂可及表面积（SASA），以量化结构波动和构象稳定性。

RMSD分析
RMSD值提供了整个模拟过程中的整体构象稳定性估计（图8a）。apo Tubulin的平均RMSD为1.82 ± 0.12 ?，表明在没有配体的情况下结构是稳定的。1d-Tubulin复合物的平均RMSD为2.45 ± 0.18 ?，大约在100纳秒后达到平衡，这反映了配体结合后的适度构象调整。相比之下，1e-Tubulin复合物的平均RMSD为2.00 ± 0.10 ?，与apo形式非常相似，表明形成了稳定的复合物。这些结果表明两种复合物都达到了平衡，其中化合物1e引起的结构偏差小于1d，意味着蛋白质-配体相互作用更加稳定。

RMSF分析
计算了每个残基的RMSF值以评估局部灵活性（图8b）。apo蛋白质显示了最高的平均残基波动（2.10 ± 0.42 ?），尤其是在靠近残基40和310的-loop区域。配体结合导致整体灵活性降低：1d-Tubulin的平均RMSF为1.78 ± 0.35 ?，而1e-Tubulin为1.85 ± 0.38 ?。1e在残基270–290和350–370附近观察到轻微的灵活性增加，表明有局部的构象适应。总体而言，配体结合稳定了蛋白质 backbone，化合物1d显示出比1e更大的刚性增强。

回转半径（Rg）分析
描述结构紧凑性的Rg值在各个系统中保持稳定（图8c）。apo蛋白质的平均Rg为21.25 ?，与良好折叠的状态一致。1d-Tubulin复合物在早期模拟中的平均Rg稍高（21.32 ± 0.11 ?），表明在稳定之前有轻微的膨胀。1e-Tubulin复合物的平均Rg较低（21.18 ± 0.08 ?），表明结构更加紧凑。这些数据表明配体结合并没有破坏整体折叠，而且化合物1e可能有助于增强结构的紧凑性。

溶剂可及表面积（SASA）分析
SASA分析（图8d）显示，配体结合相对减少了与apo形式的溶解剂接触。apo Tubulin的平均SASA为16,420 ± 95 ?2，而1d-和1e-结合的复合物的SASA值分别降低到15,980 ± 88 ?2和15,760 ± 81 ?2。配体结合系统的较低SASA值表明了更好的紧凑性和减少的灵活性，这与Rg和RMSD的发现一致。1e-Tubulin复合物表现出最低的SASA，似乎形成了更紧密和更稳定的蛋白质-配体组装。

4.7. MM/GBSA结合自由能分析
为了进一步评估配体1d和1e对Tubulin蛋白活性位点残基的结合亲和力和能量贡献，对300纳秒MD模拟轨迹进行了分子力学/广义Born表面面积（MM/GBSA）计算。表4展示了能量组成部分，如范德华（VDWAALS）、静电（EEL）、溶剂化自由能（ΔGsolv）、气相能量（ΔGgas）和总相互作用自由能（ΔGTOTAL）。两种配体都表现出有利的结合亲和力，这通过它们的负ΔGTOTAL值体现出来，表明与Tubulin结合位点的自发和热力学上是有利的相互作用。值得注意的是，配体1e显示出更有利的总结合自由能（–27.96 kcal/mol），而1d为–22.13 kcal/mol），表明化合物1e与蛋白质形成的复合物更为稳定且能量上更有利。1e-Tubulin复合物中的范德华（VDWAALS）贡献更为显著（–31.58 kcal/mol），相对于1d复合物（–16.59 kcal/mol），这突显了配体1e在结合腔内的强疏水和非极性相互作用。相反，静电成分（EEL）在1d复合物中更为有利（–38.62 kcal/mol），而在1e中为–18.74 kcal/mol），表明配体1d可能与Tubulin残基形成更极性或带电的相互作用，如氢键或盐桥。溶剂化能量惩罚（ΔGsolv）对于1e（22.37 kcal/mol）低于1d（33.08 kcal/mol），表明配体1e在结合时经历了较小的脱溶成本，这可能归因于其增加的疏水性质或与活性位点的更好兼容性。总体气相能量（ΔGgas）结合VDWAALS和EEL项，对于两种复合物来说都是可比的（1d为–55.22 kcal/mol，1e为–50.32 kcal/mol），进一步证实了两种配体在无溶剂效应的情况下都表现出强的内在结合力。总体而言，尽管配体1d更多地受益于静电相互作用，但配体1e实现了更优的整体结合亲和力，主要得益于增强的范德华相互作用和减少的溶剂化惩罚。这些发现与SASA分析结果一致，并提供了额外的能量证据，表明配体1e可能是一种更有效且结构上更兼容的Tubulin抑制剂。

药物相似性和ADMET分析
设计Tubulin抑制剂的药物相似性评估表明它们完全符合Lipinski的五条规则（见补充文件表3 [60]），该规则定义了口服活性药物的物理化学标准，如分子量（MW）≤ 500 g/mol、分配系数（Log P）≤ 5、氢键受体（nHA）≤ 10和氢键供体（nHD）≤ 5 [61]。没有一种设计分子违反了一个以上的标准，表明它们具有有利的物理化学特性，适合口服吸收和渗透。所有类似物的生物利用度得分为0.55，这表示化合物在大鼠中至少具有10%的口服生物利用度或在人体中具有可测量的Caco-2渗透性 [62]。实际上，0.55的得分被认为是中等偏高的，意味着这些化合物可能具有口服生物利用性，尽管需要进行实验确认。此外，合成可行性得分在2.8到4.7之间，表明这些设计类似物在合成上是可行的，因为低于5的值表明合成过程相对简单。设计化合物的ADME分析显示了整体上有利的药代动力学特性（见补充文件表4）。预测的人体肠道吸收（HIA）值在88.17%到96.92%之间，超过了良好吸收的30%最低阈值 [64]。稳态分布体积（VDss）值相对较低，表明化合物在组织中的分布有限，倾向于主要留在血液中。关于代谢稳定性，所有设计的抑制剂都被预测为CYP3A4的底物和抑制剂，CYP3A4是负责外源物质代谢的主要酶之一。这表明这些化合物可能会经历显著的肝脏代谢，可能会影响给药频率或潜在的药物-药物相互作用。总清除（TC）值处于可接受范围内，表明消除过程平衡，可能降低累积或快速清除的风险。补充文件表5中的PreADMET毒性预测显示，衍生物1a–1g的安全性得到了显著改善，与母体化合物相比。结构优化成功消除了所有衍生物的Ames突变性和大鼠致癌性预测。化合物1b–1g表现出最有利的安全性特征，没有突变性和致癌性警告。虽然所有化合物都显示出中等的hERG抑制风险，但没有观察到高风险的细胞毒性。然而，需要注意的是，这些ADMET和药物相似性预测仅基于计算模型，尽管对于初步筛选很有价值，但这些模型存在固有的局限性。SwissADME和pkCSM中的算法依赖于从已知数据集得出的QSAR统计相关性，可能无法完全捕捉生物系统的复杂性。因此，对代谢稳定性、酶抑制和毒性的实验（体外和体内）验证对于确认这些预测并准确建立设计Tubulin抑制剂的药代动力学和安全性特征是必要的。

结论
在这项研究中，采用了一种基于配体的设计方法来识别和设计新型的萘衍生物Tubulin抑制剂。开发并使用多重线性回归验证了一个统计上稳健的2D-QSAR模型，展示了强大的内部和外部预测能力。描述符平均效应分析指导了七个新型类似物的设计，这些类似物显示出增强的预测活性。分子对接结果显示出强烈的结合亲和力，特别是化合物1d和1e，它们的表现优于参考化合物秋水仙碱。300纳秒MD模拟证实了这些复合物的构象稳定性，其中化合物1e形成了最紧凑和稳定的复合物。MM/GBSA结合自由能计算进一步证实了这些观察结果，突显了化合物1e作为最具能量优势的配体。此外，药物相似性和ADMET分析显示出了有利的药代动力学和毒性特征，所有设计化合物都符合Lipinski的规则，并表现出良好的口服吸收和低毒性。总体而言，这些结果表明，尤其是化合物1e，代表了作为进一步临床前评估的有希望的Tubulin靶向抗癌剂的候选物。然而，重要的是要注意，目前的工作完全是计算性的。因此，为了验证这些发现，未来的研究需要进行体外实验，如细胞毒性和Tubulin聚合实验，以确认设计化合物的预测抑制潜力和生物学效果。穆罕默德·图库尔·易卜拉欣（Muhammad Tukur Ibrahim）：结果解释与审阅。
马加吉·拉丹（Magaji Ladan）和阿尤巴·阿卜杜拉希·穆罕默德（Ayuba Abdullahi Muhammad）：审阅与编辑工作。
巴班吉达·阿卜杜拉希·萨杰（Babangida Abdullahi Saje）：审阅与编辑工作。

**资金情况**
作者未为此项工作获得任何资助。

**数据与材料的可用性**
手稿已包含所有关键信息；不过，如果通过期刊及所属机构提出请求，并遵循标准流程，作者愿意提供进一步的数据。

**伦理审批与参与同意**
不适用。

**发表同意**
作者已正式同意文章的发表。

**人类与动物权利**
不适用。

**数据可用性声明**
本文描述的研究未使用任何数据。

**未被引用的参考文献**
[57, 63]