基于定量构效关系（QSAR）建模、分子对接和分子动力学模拟的计算方法，用于筛选新型二肽基肽酶-4（DPP-4）抑制剂

《Letters in Drug Design & Discovery》：Computational identification of novel DPP-4 inhibitors based on QSAR modeling, molecular docking, and molecular dynamics simulations

　　谢世建波；郑景燕；曾彦荣；张雅坤
中国西安西北大学附属医院核医学系，710000

**摘要**
背景：二肽基肽酶-4（DPP-4）抑制剂是治疗2型糖尿病的重要药物，其疗效取决于对DPP-4酶活性的有效抑制。
**目的**：本研究旨在构建可靠的QSAR模型，并设计新型尿嘧啶–苯甲酸衍生的DPP-4抑制剂，以提高预测活性和优良的药代动力学特性。
**方法**：利用25种已报道的DPP-4抑制剂，开发了HQSAR和Topomer CoMFA模型。通过原子贡献和等高线分析确定关键结构片段，并通过虚拟筛选重组生成20个新候选分子。进行分子对接和分子动力学模拟以评估结合亲和力和稳定性，并预测ADMET特性以评估药代动力学潜力。
**结果**：HQSAR和Topomer CoMFA模型表现出良好的预测性能（验证系数分别为0.778/0.915和0.719/0.824）。这20种新设计的化合物的预测活性高于模板分子。对接结果显示与关键活性位点残基有良好的氢键和疏水相互作用，分子动力学模拟证实了蛋白-配体结合的稳定性。ADMET预测表明具有可接受的药代动力学特性。
**结论**：这些综合计算分析为2型糖尿病新型DPP-4抑制剂的合理设计提供了理论支持。

1. **引言**
糖尿病是一组慢性代谢紊乱疾病，由胰岛素分泌不足、胰岛素抵抗不足或两者共同作用引起，最终导致脂肪、蛋白质和碳水化合物代谢紊乱。未经控制的糖尿病可引发转移性癌症、肾病、神经病变、心血管疾病、视网膜病变和肥胖。随着生活水平的提高，糖尿病发病率逐年上升。截至2021年，20–79岁成人糖尿病患者总数达到5.37亿。国际糖尿病联合会（IDF）的最新数据显示，到2030年全球糖尿病患者将增至6.43亿，到2045年将增至7.83亿。这已成为增长最快的公共卫生威胁。在临床实践中，1型糖尿病（胰岛素依赖型糖尿病，T1DM）和2型糖尿病（非胰岛素依赖型糖尿病，T2DM）是两种主要类型，其中T2DM占全球糖尿病病例的90%。目前，T2DM的治疗是现代社会面临的最昂贵的医疗问题之一。尽管已有治疗方法，但血糖控制仍对T2DM患者构成重大挑战。新的T2DM治疗手段包括胰高血糖素样肽-1（GLP-1）激动剂、α-葡萄糖苷酶抑制剂、钠依赖性葡萄糖转运蛋白-2（SGLT-2）抑制剂和二肽基肽酶-4（DPP-4）抑制剂。目前，大多数研究认为DPP-4抑制剂是T2DM药物研究的新热点，可能成为糖尿病治疗的革命性进展。DPP-4是一种能够降解葡萄糖依赖性胰岛素促分泌多肽（GIP）和胰高血糖素释放激素（GLP-1）的酶。GIP和GLP-1通过促进胰岛素分泌和抑制胰高血糖素释放来降低血糖水平。DPP-4抑制剂的优势在于耐受性良好，但也存在一些缺点，如可能引起头痛和胃肠道问题。因此，迫切需要开发具有新靶点的新药物或补充现有治疗方案。
本研究采用计算机辅助药物设计（CADD）方法，设计出活性高、疗效好且副作用少的新型DPP-4抑制剂。CADD方法的主要步骤包括建立定量结构-活性关系（QSAR）模型、设计新分子、进行分子对接和分子动力学（MD）模拟。该方法因开发周期短、风险低而在药物研发领域得到广泛应用。QSAR利用化学计量技术建立结构与功能之间的计算或数学模型，主要揭示药物分子与生物分子之间的关系。分子对接是分析配体-蛋白相互作用的有效工具，有助于深入理解小配体分子与蛋白质大分子之间的相互作用。ADMET预测确保新设计分子在人体内的良好药代动力学特性，因此在体外和体内测试前评估候选分子的安全性和毒性风险至关重要。分子动力学模拟可以直观展示药物分子与蛋白质的结合动态。

2. **材料与方法**
2.1 **数据集**：从文献中筛选出25种具有相似骨架但生物活性不同的尿嘧啶苯甲酸及其酯类衍生物作为DPP-4抑制剂。将它们的半数抑制浓度（IC50）值转换为负对数常数pIC50值。表S1展示了采用这些化合物的化学结构和生物活性值的QSAR模型分析结果。排除了一些活性较低的分子，选择具有多样化学结构和广泛生物活性的化合物作为训练集和测试集。
2.2 **HQSAR建模**：HQSAR建模是一种2D-QSAR方法，将分子的生物活性与其结构片段相结合，并根据不同原子团对生物活性的贡献来优化分子结构。首先将分子分解为不同大小的片段，片段大小和长度是关键参数。片段区分基于原子（A）、键（B）、连接（C）、手性（CH）、氢原子（H）以及供体和受体（DA）等基本参数。这些参数被认为是影响QSAR结果的最小单元。片段编码为分子全息图，并应用偏最小二乘（PLS）回归算法基于全息图描述符和生物活性数据构建QSAR模型。
2.3 **Topomer CoMFA建模**：Topomer CoMFA是一种新型3D-QSAR方法，可预测化合物的生物活性。首先将化合物切割成两个或多个小片段，同时保留尽可能多的共同骨架。然后通过自动系统运行生成所有片段的立体和静电场势能值，随后进行偏最小二乘回归分析，建立定量结构-活性关系。
2.4 **偏最小二乘（PLS）分析**：利用偏最小二乘算法构建QSAR模型。该算法量化模型描述符与生物活性之间的关系。
2.5 **QSAR模型验证**：理想的QSAR模型不仅具有良好的内部验证结果（高q2和r2值），外部验证也同样重要，用于评估模型的预测能力。通过计算未用于模型构建的化合物的生物活性来验证模型的预测能力。外部测试参数rpred2用于评估模型预测能力。
2.6 **Y-随机化测试和适用范围验证**：为评估所构建QSAR模型的稳定性并排除偶然性，进行了Y-随机化测试。Y-随机化验证基于随机扰动实验值、保持变量不变后再构建模型的原理，比较结果与初始建模结果。适用范围评估是判断模型能否准确预测新药性质或活性的关键指标。目前，已经开发了多种方法来评估模型的适用范围。在本研究中，采用了广泛使用的杠杆值方法来确定QSAR模型的适用范围，其定义如公式（13）所示：(13)hi=xiT(XTX)?1xi，其中h是杠杆值，X是描述符矩阵，XT是X的转置，xi是化合物的描述符行向量。通过与阈值（h*）进行比较，如果hi < h*，则表示该模型在该区域的预测性是稳定的；相反，则表明该模型的预测能力较差。阈值定义如下：(14)h*=3p′/n，其中p′是描述符的数量；n是训练集化合物的数量。

2.7 分子筛选和活性评估
虚拟筛选是一种用于药物发现中识别先导化合物的方法。拓扑搜索是一种基于片段的筛选工具，能够快速从化学数据库中搜索和组装小分子。在本研究中，选择了包含大量生物活性小分子的ZINC药物数据库作为基于R基团的拓扑搜索的源数据库。在分子设计过程中，首先根据预定义的碎片规则将参考活性化合物系统地分解为核心骨架和相应的R基团片段。随后，从数据库中检索出符合目标要求的拓扑特征的R基团片段，并使用建立的Topomer CoMFA模型对其取代基的贡献进行定量评估和排序。在此基础上，通过将高贡献的R基团片段与选定的通用核心骨架重组来构建新分子。这种基于片段的重组策略有效地固定了中心骨架，同时在特定位置系统地变化取代基，从而能够高效地探索化学空间并进行结构优化，生成具有潜在增强生物活性的候选化合物库。最后，使用经过验证且稳健的Topomer CoMFA模型对新设计的分子的生物活性进行预测。
值得注意的是，我们选择满足以下两个条件的分子片段来设计新分子：1）分子片段的拓扑距离应接近185；2）R基团的贡献值大于Topomer CoMFA模型中相应R基团的贡献值。从相应的高贡献基团组合中设计出新分子后，使用建立的稳定可靠的Topomer CoMFA模型对其生物活性进行了预测。

2.8 分子对接
分子对接是一种理论模拟方法。该方法能够显示配体小分子与蛋白质之间的相互作用，并预测它们之间的结合模式和亲和力。存在多种分子对接方法，不同方法得到的对接结果也不同。因此，只有选择最合适的方法才能显示小分子配体与蛋白质之间最真实的相互作用。有三种常见的对接方法，本研究选择了最合适的方法进行对接分析：
Surflex-Dock：这是SYBYL-X 2.0软件中对接模块使用的方法。它利用独特的经验评分函数和基于分子相似性的搜索引擎将配体分子与蛋白质结合位点进行对接。该方法基于使用原型分子表示蛋白质的结合口袋，并使用探针来探测蛋白质口袋的疏水性、氢键和静电性质，从而生成蛋白质活性口袋的图像。
AutoDock：AutoDock基于经验自由能场和快速拉马克遗传算法。其中有两个主要部分起关键作用。首先，AutoGrid通过预处理目标蛋白质并定义结合位点来生成网格图；随后，AutoDock在准备好的网格内进行配体对接以评估结合相互作用。
AutoDock Vina：Vina是一个基于AutoDock 4的即用型分子对接程序。其主要特点是使用快速梯度优化的概念搜索算法和简单的评分函数。研究表明AutoDock Vina提高了组合模型的平均预测准确性。
综合所有因素，本研究选择了AutoDock Vina方法来对接分子。该方法用于将半柔性配体与刚性蛋白质结构进行对接。第一步是对蛋白质和配体进行预处理，即最小化配体的能量并从蛋白质中去除水分子。这些操作的参数设置如下：使用Kollman联合力场，非键合截止值设为8.0，介电常数设为1.0。使用AutoGrid模块生成尺寸为50 × 50 × 50 ?的结合能网格，网格间距为0.375 ?，中心坐标为X = 51.099, Y = 49.464, Z = 36.186。对接后检查计算出的结合能以获得最佳结构构象。

2.9 分子动力学模拟
为了进一步研究配体分子与目标蛋白质之间的相互作用以及复合物在生理动态条件下的稳定性，对选定的复合物进行了分子动力学（MD）模拟。MD模拟是一种成熟且高效的计算方法，广泛用于分析生物相关分子的动态行为并阐明结构-功能关系。所有分子动力学模拟均使用GROMACS 5.15进行。蛋白质的拓扑文件和力场参数使用pdb2gmx模块和CHARMM（2021年7月）力场生成。配体参数基于CHARMM通用力场（CGenFF）进行赋值。积分时间步长设置为2 fs。所有涉及氢原子的共价键均使用LINCS算法进行约束。系统置于周期性边界条件（PBC）十二面体盒子中，盒子内填充TIP3P水模型，保持蛋白质与盒子边缘之间的最小距离为1.0 nm。添加钠离子和氯离子以中和系统，达到0.15 M的生理离子强度。长程静电相互作用使用Particle Mesh Ewald（PME）方法处理，实空间截止距离为1.0 nm。范德华相互作用也使用1.0 nm的截止距离进行计算，并应用了力切换修改以平滑处理截止附近的相互作用。能量最小化使用最速下降算法进行了50,000步，直到最大力小于1000 kJ·mol?1·nm?1，确保起始构象的稳定性。随后系统在NVT和NPT系综下分别平衡100 ps，温度保持在300 K，使用V-rescale恒温器，压力使用Parrinello–Rahman压力计保持为1 atm。最后，在恒定温度和压力条件下对每个选定的复合物进行了100 ns的生产模拟。模拟过程中生成的拓扑和轨迹文件使用GROMACS的内置程序进行分析，包括均方根偏差（RMSD）、均方根波动（RMSF）、回转半径（Rg）、可溶剂表面积（SASA）、氢键数量以及吉布斯自由能景观（FEL）的分析。基于MD模拟获得的结果文件，并使用分子力学/泊松-玻尔兹曼表面积（MM/PBSA）方法，我们可以计算新设计分子与蛋白质之间的结合自由能。每个系统的结合自由能通过MM/PBSA方法使用gmx_MMPBSA计算MD模拟期间最后5 ns的稳态轨迹。结合自由能ΔGbind根据公式（15）计算：(15)ΔGbind=ΔEvdW+ΔEelec+ΔEpolar+ΔEnonpolar?TΔS，其中ΔGbind表示结合自由能，ΔEvdW表示范德华相互作用，ΔEelec表示静电相互作用，ΔEpolar表示极性溶剂化-溶剂化相互作用，ΔEnonpolar表示非极性溶剂化-溶剂化相互作用，TΔS表示温度T下的构象熵贡献。由于TΔS对复合物的结合自由能影响较小，我们在计算中将其省略。

2.10 ADMET预测
ADMET药代动力学性质的预测是当前药物设计和开发过程中的关键步骤。在临床前阶段，ADMET预测可以克服物种差异的问题，减少药物毒性和副作用，从而更好地指导药物在临床中的合理使用。这一操作的结果可以提高药物开发的成功率，降低开发成本，并避免不必要的时间和资源浪费。因此，我们对设计的新分子进行了ADMET性质预测。预测包括人体肠道吸收（HIA）、口服生物利用度（HOA）、血脑屏障渗透（BBB）和总清除率。新设计分子的ADMET性质使用在线平台ADMET Lab和admetSAR进行预测。

3. 结果与讨论
3.1 HQSAR模型结果与分析
影响HQSAR模型构建的关键因素是全息长度（HL）、片段区分（FD）和片段长度（FL）。因此，我们需要选择和设置这些参数。在本研究中，最初应用了默认的片段长度（FL，4–7），并通过不同的HL和FD组合获得了具有更高预测性能的模型。通过组合各种片段参数，共构建了53个HQSAR模型，详细结果总结在表S2中。研究发现，通过考虑参数A、B、C、H和CH获得了最佳的HQSAR模型。其参数为q2 = 0.778, r2 = 0.915, SEE = 0.180, HL = 199, N = 4。随后，保持最佳的HQSAR片段区分类型不变，评估不同的片段长度以获得具有改进预测性能的模型。表S3列出了最佳模型50–4的计算参数，结果显示q2 = 0.778, r2 = 0.915, SEE = 0.180, HL = 199, N = 4。表S4列出了HQSAR模型50–4为每种化合物计算的预测值，结果显示大多数化合物的实际值与预测值一致。图1a显示了基于HQSAR模型的实验值和预测值之间的线性回归图，显示出良好的线性关系。图1b显示了实验值和HQSAR模型预测值之间的残差图。从图中可以看出，所有残差值都小于±0.6，这表明模型的结果没有系统误差。

3.2 Topomer CoMFA模型结果与分析
Topomer CoMFA模型基于将分子切割成不同大小的片段，计算每个片段的贡献值，该贡献值与分子的活性密切相关，最终计算活性与QSAR之间的关系。训练集中最活跃的分子（化合物14）被选为碎片化的模板。表1中显示的两种碎片化方案被应用，剩余的分子使用相同的切片策略自动进行碎片化处理，从而得到了一个更加稳健且具有预测性的Topomer CoMFA模型。表1的结果表明，通过模型2的切割方法构建的QSAR模型具有良好的适应性和预测能力。该模型的统计参数为：q2 = 0.719，r2 = 0.824，F = 37.449，N = 2，SEE = 0.240，rpred2 = 0.783。这些模型参数符合稳定模型的标准，说明我们选择的切割方法构建的Topomer CoMFA模型更为有效。利用碎片化模型2，获得了Topomer CoMFA模型的统计参数并随后进行了评估。基于这些结果，生成了Topomer CoMFA模型的回归线和残差图。图3a显示了Topomer CoMFA模型预测的分子活性与其实验活性之间的关系。从图中可以看出，大多数分子在回归线的两侧对称分布。图3b显示了分子实验值与预测值之间的残差。图中零线两侧的残差范围较小，表明残差值中没有系统性错误。因此，可以得出结论，Topomer CoMFA模型表现出很强的预测性能。

表1. Topomer CoMFA模型的计算分子切割和参数。
| No. | Cutting Model | NSE | EFE | q2 | r2 | rpred2 |
|-------------|-----------------|--------------------|------------------|------------------|------------------|
| 1 | 2 | 0.230 | 46.355 | 0.70 | 0.827 |
| 2 | 2 | 0.240 | 37.449 | 0.719 | 0.824 |
| Ra fragments | Red | Rc fragments | Blue | Common backbone | Black |
| Rb fragments | Green | | | |
| | | | | |

下载：下载高分辨率图像（108KB）
下载：下载全尺寸图像

图3. 基于Topomer CoMFA模型的实验值与预测值的线性回归图（a）以及残差与实验活性的关系图（b）。

根据Topomer CoMFA模型，为每个片段组生成了等高线图。等高线图分为立体场和电场。在Topomer CoMFA的等高线图中，结构被划分为有利区域和不利区域，这为抑制剂设计的进一步优化提供了参考。在立体场等高线图中，绿色区域表示较大取代基的活性较高；在电场等高线图中，红色区域表示带负电荷的取代基对活性有利；蓝色区域表示带正电荷的取代基更活跃。

活性最强的化合物14（pIC50 = 9.097）被选为分析参考。图4展示了化合物14的Ra片段空间场。氮原子周围的黄色区域表明添加小而 bulky 的取代基可以增强分子的活性，这解释了化合物2（pIC50 = 8.036）的活性高于化合物4（pIC50 = 7.695）。图4d显示了化合物14的Ra片段的电场，其周围的红色大块表示添加带负电荷的片段可以增加分子的活性。图4b显示了化合物14的Rb片段的空间场，溴原子周围的绿色区域表明添加该片段可以增强分子的活性。图4e显示了化合物14的Rc片段的电场，溴原子周围的红色区域表明添加带正电荷的片段会降低分子的活性。图4c显示了化合物14的Rc片段的空间场，羧基周围的绿色区域表明增加该片段的体积会降低分子的活性。图4f显示了化合物14的Rc片段的电场，羧基周围的红色区域表明添加电负性的片段会增强分子的活性。

下载：下载高分辨率图像（62KB）
下载：下载全尺寸图像

图4. 化合物14 Topomer CoMFA模型的三维等高线图：(a-c) 立体场；(d-f) 电场。

3.3. QSAR模型的验证
已建立的QSAR模型的q2和r2参数具有很强的内部预测能力。使用测试集计算了HQSAR和Topomer CoMFA模型的外部验证参数，相应结果展示在表2中。两个模型的参数结果均符合标准，表明这两个模型是稳定的，并且具有一定的预测能力。

表2. QSAR模型验证方法的统计参数
| Parameter | Topomer CoMFA | HQSAR |
|------------------|------------------|----------------------|
| q2 | > 0.5 | 0.719 |
| r2 | > 0.6 | 0.824 |
| SEE | 0.224 | 0.177 |
| N | 24 |
| External validation | r2pred | > 0.6 | 0.783 |
| Golbraikh-Tropsha method | R2 | > 0.6 | 0.778 |
| R02 | 0.776 | 0.808 |
| R'02 | 0.675 | 0.782 |
| (R2 - R02) / R2 | < 0.1 | 0.003 |
| (R2 - R′02) / R2 | < 0.1 | 0.097 |
| k | 0.85 < k < 1.15 | 1.000 |
| k' | 0.85 < k' < 1.15 | 0.999 |
| rm2 | rm2 > 0.5 | 0.738 |
| ? rm2 | ? rm2 < 0.2 | 0.192 |
| RMSE | Close to 0 | 0.255 |
| MAE | 0.244 | 0.179 |
| CCCC | > 0.85 | 0.867 |

3.4. 随机化测试和适用域验证
对HQSAR和Topomer CoMFA模型进行了Y-随机化测试。以生物活性作为因变量，并随机重新分配其值以构建新的QSAR模型进行验证。20次随机化测试的结果显示，两个模型的q2和r2值均较低，表明所得模型并非基于偶然性。

在本研究中，使用杠杆方法检测了计算模型的应用领域，结果显示在图5中。对于训练集和测试集中的25种有机物，所有分子的杠杆值均在指定阈值范围内，这表明该方法构建的模型具有广泛的应用范围。

下载：下载高分辨率图像（124KB）
下载：下载全尺寸图像

图5. HQSAR模型适用域地图（a）和Topomer CoMFA模型适用域地图（b）。

3.5. 分子设计与活性预测
利用Topomer CoMModel的结果进行分子设计。通过Topomer搜索模块从ZINC数据库中识别并选择贡献值较高的分子片段，然后将某些最优分子片段与共同骨架结合形成新的分子。每个基序的最终筛选片段及其对应的贡献值显示在图S1中，包括贡献值高且结构匹配的Ra片段、Rb片段和Rc片段。通过将这些筛选出的单个基序的分子片段与共同骨架结合，我们获得了20种新的抑制剂分子。表S6列出了它们的结构以及基于Topomer CoMModel的预测活性值。所有新设计分子的预测pIC50值均高于模板分子14，这表明所构建的模型具有说服力，并为DPP-4抑制剂的研发提供了有用信息。然而，需要明确指出的是，目前的发现仅基于计算预测，尚未经过实验验证。因此，这些新设计的分子仍需要进一步的化学合成以及全面的体外和体内生物活性评估以确认其实际的抑制效果。

3.6. 分子对接
分子对接提供了关于设计化合物在目标蛋白（PDB:2RGU45）活性位点结合模式的见解，包括关键的氢键、疏水作用和静电相互作用。这有助于与已知的DPP-4抑制剂进行比较，并支持机制解释和结构-活性关系分析。原配体与重新对接配体之间的RMSD值是确保对接可靠性的关键指标。当RMSD ≤ 2.0 ?时，我们认为对接是可靠的。原配体是从蛋白质结构中提取的，然后使用AutoDock Vina进行重新对接。重新对接的结果显示在图6中，其中原配体显示为绿色，重新对接的配体显示为蓝色。重新对接的配体与原配体几乎完全相同，它们的旋转趋势也相同。原配体与重新对接配体之间的RMSD为0.444 ?，表明所选择的对接方法是可行的，配体在对接前后的构象没有显著变化，且配体与蛋白质的结合模式是合理的。

图7. 模板化合物14与蛋白质的对接结果。配体与Tyr547、Glu206和His740氨基酸形成氢键相互作用，并与Tyr631形成疏水相互作用。对接过程的结合能为为-8.6 kcal/mol。

图6. 模板化合物14与蛋白质的对接结果显示，配体与Tyr547、Glu206和His740氨基酸形成氢键，与Tyr631形成疏水相互作用。通过AutoDock Vina方法对接所有新设计分子的结合能得分显示在表S6中，所有化合物的总得分范围为-8.5至-10.3 kcal/mol，表明新设计的分子能够很好地结合到蛋白质的结合位点，从而增加了抑制DPP-4的可能性。基于良好的对接结果，选择了三种新分子N08、N16和N20。它们与蛋白质的对接结果如图8所示。N08与蛋白质的对接结果显示，Asp545残基与Rb取代基上的H原子形成氢键。此外，Trp629和His740氨基酸残基分别与Rc取代基上的O和H原子形成氢键。N16与蛋白质的对接结果显示，N16主要通过氢键与蛋白质稳定结合。N20与蛋白质的对接结果显示，Tyr666与Rc取代基上的H原子形成氢键。Val546与Rb取代基上的H原子形成氢键。N20与蛋白质的对接结果显示，Tyr666与Rb取代基上的H原子形成氢键。

图8. 新分子N08与蛋白质对接后的结合构象。
图9. 新分子N16与蛋白质对接后的结合构象。
图10. 新分子N20与蛋白质对接后的结合构象。

配体-蛋白质结合的稳定性与某些残基密切相关。从各个分子的结合模式图中可以看出，配体分子与Tyr666周围的残基形成疏水相互作用，而蛋白质也与Tyr547和His740等重要残基形成氢键。与模板分子14相比，新分子与蛋白质的结合更为稳定。这些结果表明我们新设计的分子与蛋白质的结合更为牢固。此外，还表明DPP-4抑制剂的活性结合位点的主要氨基酸残基是Tyr547、Tyr666和His740。这些理论发现为新型DPP-4抑制剂的开发提供了有价值的见解。

总之，化合物N08、N16和N20在DPP-4活性位点显示出有利的结合构象，与Tyr547、Tyr666、His740等关键残基形成了多个氢键，并具有稳定的疏水相互作用。它们的结合能表明它们对目标蛋白质有很强的亲和力。与模板化合物14相比，这些分子展示了更强的相互作用网络和更稳定的结合模式，显示出作为有前景的DPP-4抑制剂候选物的潜力。

3.7. 分子动力学模拟分析
进行了分子动力学模拟，以研究N08、N16和N20与目标蛋白质的结合稳定性。在模型构建过程中还进行了模板分子复合物的比较模拟。在进行了100纳秒的分子动力学模拟后，使用GROMACS内置程序分析了轨道和坐标文件。为了探索受体-配体复合物的性质，如柔韧性和稳定性，采用了RMSD、RMSF、Rg、SASA和氢键数量等客观指标。如图11所示，通过RMSD分析来评估100纳秒分子动力学模拟过程中配体-蛋白质复合物的稳定性。所有系统在0-30纳秒内都表现出初始波动，这对应于构象调整，随后逐渐达到平衡。30纳秒之后，RMSD曲线稳定下来并保持相对恒定。所有复合物的主链RMSD值（图11d）都保持在较低水平，表明整体结构稳定。具体来说，化合物N08（图11a）在最初30纳秒内仅有轻微波动，并在剩余的模拟时间内保持稳定。化合物N16（图11b）显示出相对较高的RMSD值，表明其具有更大的柔韧性，但在大约36.2纳秒时达到稳定，平均RMSD为0.263纳米。化合物N20（图11c）在0-10纳秒内RMSD略有增加，之后在30纳秒时轨迹稳定。

相比之下，模板化合物14在60至80纳秒之间持续波动，显示出较弱的动态稳定性。总体而言，新设计的配体（N08、N16和N20）比模板化合物14表现出更稳定的结合构象，这支持了它们的良好结合特性。

图12展示了复合物的RMSF曲线，RMSF表征了配体结合对蛋白质链每个氨基酸残基的影响。较大的RMSF值表明蛋白质结构的柔韧性和松散的键合；相反，较小的值表明更高的稳定性和二级结构的存在。在所有复合物中，氨基酸残基在较低范围内以相似的模式波动，这表明配体分子在结合位点的作用方式相似。从图中可以看出，复合物结构保持高度稳定，主要是由于配体分子与两个关键氨基酸残基Tyr547和Tyr666之间的氢键和疏水相互作用。此外，Tyr629和Tyr666是产生氢键相互作用的重要氨基酸，配体与它们的结合使复合物的RMSF保持较低。在蛋白质结合位点内，662至666之间的氨基酸形成了重要的疏水腔，与配体分子形成强疏水键。这些关键相互作用使复合物的RMSF保持在较低水平。Asp243至Gln247之间的氨基酸显示出最高的波动。这些波动是因为它们位于蛋白质结构的边缘，远离结合位点，并且不与小分子形成明确的相互作用。

Rg也是在动力学模拟过程中生成的文件中的一个重要参数。它可以确定模拟过程中复合物整体构象的稳定性和松弛程度。一般来说，较大的Rg值表示蛋白质结构较为松散；而较小的Rg值表示蛋白质结构更为紧凑。

SASA指的是直接与溶剂分子接触并产生效应的复合物表面积。SASA的增加表明复合物的结构已经扩展，增加了与溶剂分子接触的面积。SASA的减少表明由于结构紧缩，与溶剂接触的表面积减小。如图13b所示，所有复合物的溶剂可接触表面积（SASA）在模拟过程中保持相对恒定。它们的SASA值在385至400纳米2的范围内。这表明复合物结构在模拟过程中相对稳定和平衡。

配体与蛋白质的结合受到配体分子与蛋白质之间氢键存在的影响。因此，可以通过分析模拟过程中配体与蛋白质之间形成的氢键数量来评估配体结合的稳定性。图14d显示了所有复合物的氢键数量。图14a、b和c分别显示了化合物N08、N16和N20在结合位点产生的氢键数量。在整个模拟过程中，N08有3至4个氢键以维持稳定的极性相互作用，并在模拟期间最多生成7个氢键。化合物N08、N16和N20在100纳秒的模拟时间内分别与蛋白质生成1.906、1.911和1.936个平均氢键。相比之下，模板分子14在模拟期间仅生成1至2个氢键，平均氢键数为1.186，低于新化合物的平均氢键数。蛋白质-配体相互作用的氢键分析表明，新设计的分子与蛋白质的相互作用比模板分子更强。

为了理解复合物构象的结构变化并识别最低能量的构象，可以使用FEL（自由能景观分析）。深蓝色区域表示复合物处于最低能量状态。颜色块的分散程度表示不同能量状态下的构象灵活性。如果配体-蛋白质相互作用弱或不稳定，则会显示多个最小能量簇；而更强和更稳定的相互作用则在能量分布中显示一个单一的构象簇。如图15、图16、图17所示，图a和b分别展示了2D和3D自由能图。可以看出，复合物在图中显示出一个显著的蓝色区域，表明其具有强且稳定的结合构象。

图15和图16、图17分别展示了新设计化合物N08的自由能景观，这些自由能景观是在RMSD和Rg坐标之间绘制的。图a和b展示了2D结果图，图c和d展示了最低能量构象相互作用的结果。

图c和d展示了配体分子在最低能量结合位点和蛋白质静电势面上的结合构象，稳定的结合姿势与我们之前的研究结果一致。化合物N08的自由能范围为0至12.628千焦/摩尔，N16和N20的自由能范围为0至12.967千焦/摩尔。

3.9 MM/PBSA结合自由能计算
结合自由能是配体-蛋白质结合稳定性的另一个关键参数。它通常包括ΔEelec、ΔEpolar和ΔEnonpolar。新分子（N08、N16和N20）以及模板分子14与蛋白质结合的结合自由能分别显示在图18中。结合能量的负值表明对稳定性有益，而正值则具有不利影响。化合物N08、N16、N20和模板分子4的结合能量分别为-84.129千焦/摩尔、-79.215千焦/摩尔、-70.270千焦/摩尔和-48.119千焦/摩尔（ΔGbind）。新分子的结合能量总是高于模板分子的结合能量。范德华相互作用对结合能量的贡献最大，其次是静电能和SASA能量，其中负值对结合能量有益。

3.10 ADMET预测分析
对20种新设计的分子进行了全面的体外ADMET预测（表S7）。所有化合物都表现出高的人体肠道吸收（HIA，88.83%–99.19%）和中等程度的口服生物利用度（HOB），表明具有良好的口服药物潜力。理想的Caco-2渗透性进一步支持了它们的高效肠道传输能力。预测的血液-脑屏障（BBB）渗透性较低，表明中枢神经系统暴露有限，这可能降低外周作用DPP-4抑制剂的CNS相关副作用的风险。几种化合物显示出相对较高的血浆蛋白结合（PPB），这可能有助于延长系统保留时间，前提是结合是可逆的。所有分子都被预测为潜在的CYP3A4底物，表明肝脏代谢是主要的消除途径。大多数化合物显示出中等的预测清除率，反映了代谢稳定性和消除率之间的合理平衡。然而，应进一步考虑CYP3A4介导的药物-药物相互作用的可能性。毒性预测表明，大多数化合物符合评估的安全标准，表明理论上的副作用风险相对较低。尽管如此，这些发现仅基于计算模型，需要实验验证。总体而言，ADMET谱支持这些分子的药物特性，并为它们的进一步合成、生物评估和作为潜在DPP-4抑制剂的优化提供了合理的依据。

4. 讨论
在这项研究中，建立了具有强大内部和外部验证的预测模型，并应用于新型DPP-4抑制剂的合理设计。与之前的QSAR研究相比，HQSAR和Topomer CoMFA的整合提供了更全面的原子和片段级别的结构-活性关系描述。从贡献图和等高线分析中识别出的关键结构决定因素与已知的DPP-4抑制特征一致，特别是疏水相互作用和取代基效应。在经过验证的Topomer CoMFA模型的指导下，片段重组能够高效生成结构新颖的候选物，并提高预测的活性。分子对接揭示了与DPP-4抑制剂的已知相互作用模式一致的关键活性位点残基和结合模式。分子动力学模拟进一步确认了配体-蛋白质复合物在动态条件下的结构稳定性。ADMET预测表明，大多数设计化合物表现出符合一般药物特性标准的可接受药代动力学性质。

然而，这项研究仅依赖于计算机模拟方法，缺乏实验验证，这给预测的生物活性带来了不确定性。此外，QSAR建模、分子对接和分子动力学模拟的固有局限性，包括数据集依赖性、评分不准确性以及力场近似可能会影响预测的可靠性。因此，ADMET结果应被视为初步估计。

未来的工作将侧重于实验验证，包括酶学、细胞和体内研究。在实验反馈的指导下，结合更严格的计算方法（如自由能计算），将提高所提出化合物的稳健性和转化潜力。

5. 结论
总之，这项研究建立了可靠的预测模型，并实施了用于新型DPP-4抑制剂合理设计的综合计算策略。确定了控制抑制活性的关键结构特征，并生成了一系列具有改进预测活性和有利ADMET特性的候选化合物。尽管目前缺乏实验验证，但所提出的框架为DPP-4抑制剂的设计提供了一种系统和高效的方法，并为后续的糖尿病药物开发实验研究提供了理论基础。所有作者均阅读并批准了最终稿件。

作者贡献声明：
曾雁蓉（Yan-Rong Zeng）：撰写、审稿与编辑、数据可视化、结果验证。
闫静（Jing Yan）：研究方法论设计、形式化分析。
张亚坤（Ya-Kun Zhang）：撰写初稿。
史乐（Le Shi）：审稿与编辑。
童建波（Jian-Bo Tong）：审稿与编辑。

伦理审查与参与同意：
不适用。

资金支持：
本研究得到了国家自然科学基金（项目编号：22373062）和陕西科技大学研究生创新基金的支持。