发展可极化键偶极模型快速预测环肽分子构象稳定性

Extending the Polarizable Bond-dipole Model to Enable the Rapid Prediction of the Conformational Stability of Cyclic Peptides

摘要

环肽分子具有独特的结构稳定性、 多样的生物活性和良好的靶向性， 是药物研发的重要先导化合物. 快速准确预测环肽构象稳定性， 不仅有助于揭示蛋白质错误折叠的分子机制， 还可为靶点识别与干预提供理论依据， 对理性设计结构稳定、 活性优良的环肽类药物具有重要意义. 本文将环肽分子中主链上C=O, N—H, C ^α —H化学键和侧链上C—O， O—H化学键视为键偶极， 使用固有偶极-固有偶极作用描述体系中的静电作用， 使用固有偶极-诱导偶极作用以及诱导偶极-诱导偶极作用描述极化作用， 并引入键长、 键角和二面角作用项描述成键作用， 从而将可极化偶极-偶极作用模型发展为可预测环肽分子不同构象相对稳定性的分子力场势函数. 将所得势函数应用于数据库中9个环肽分子共计33个不同构象， 快速计算其构象能， 并与高精度DLPNO-MP2/aug-cc-pVTZ方法和可极化分子力场AMOEBA方法计算所得构象能进行比较. 结果表明， 相对于DLPNO-MP2/aug-cc-pVTZ方法的环肽分子构象能， 本文方法的线性相关系数R²=0.9784， 均方根偏差RMSE=13.43 kJ/mol， 略优于AMOEBA方法（R²=0.9682， RMSE=16.28 kJ/mol）. 对环肽分子的结构优化和频率计算结果进一步表明了本文方法的合理性. 相较于AMOEBA可极化分子力场方法， 本文方法显著降低了静电作用项数.

Abstract

Cyclic peptides possess unique conformational stability， diverse biological activities， and favorable target specificity， making them important lead compounds in drug development. Rapid and accurate prediction of their conformational stability not only aids in uncovering the molecular mechanisms of protein misfolding， but also provides a theoretical basis for target identification and intervention. This is of great significance for the rational design of structurally stable and highly active cyclic peptide-based drugs. In this paper， the polar chemical bonds C=O， N—H， C ^α —H and C—O， O—H in cyclic peptides are regarded as bond dipoles. The permanent dipole- permanent dipole interaction is used to describe the electrostatic interaction in the system， and the permanent dipole-induce dipole interaction and induce dipole-induced dipole interaction are used to describe the polarization. The bonded terms， including the bond-stretching， angle-bending， and dihedral torsion， are also introduced. The polarizable dipole-dipole interaction model is thus developed into a potential function that can be used to rapidly calculate the relative energies of different conformations of cyclic peptides. The potential function is applied to 9 cyclic peptides， total 33 different conformations to rapidly predict the conformational energies of these conformations. The conformational energies of these conformations are also calculated using the AMOEBA and DLPNO-MP2/aug-cc-pVTZ methods. The calculation results show that， compared with the DLPNO-MP2/aug- cc-pVTZ conformational energies， the linear correlation coefficient R² of our model is 0.9784， and the root mean square deviation is 13.43 kJ/mol， slightly better than the linear correlation coefficient 0.9682 and root mean square deviation 16.28 kJ/mol of the AMOEBA method. The results of structural optimization and frequency calculation further suggest the rationality of our model. Furthermore， compared with the AMOEBA polarizable force field， our polarizable model significantly reduces the number of electrostatic terms. The model proposed in this paper may provide a new tool for the research and development of novel cyclic peptides as drug candidate molecules.

Graphical abstract

关键词

环肽 / 构象能 / 化学键偶极 / 诱导键偶极 / 可极化势函数

Key words

Cyclic peptide / Conformational energy / Chemical bond dipole / Induced bond dipole / Polarizable potential energy function

环肽具有独特的结构稳定性、 多样的生物活性和良好的靶向性， 是药物研发的重要先导化合物^［1~8］. 理性从头设计在药物研发过程中起着越来越重要的作用， 它能够针对特定靶点进行精准设计， 提高药物发现的成功率， 缩短研发周期， 降低研发成本^［8~14］. 发展高精度高效率理论计算方法， 准确预测环肽分子不同构象的构象能， 有助于理解蛋白质错误折叠机制， 寻找干预靶点， 对理性从头设计药物候选物有重要价值. 准确的分子构象能通常需要使用高精度量子化学方法计算获得^［15~24］. 通常将CCSD（T）/CBS方法计算得到的分子构象能作为检测其它计算方法计算精度的标准数据^［15~19］. 但是 CCSD（T）/CBS方法十分昂贵， 仅适用于较小分子体系. 对于较大分子体系， 通常使用二阶微扰理论MP2/aug-cc-pVTZ方法以及高精度密度泛函理论B2PLYP-D3BJ/aug-cc-pVTZ或LC-ωPBE-XDM/aug-cc-pVTZ等方法计算得到较为可靠的分子构象能^［20~24］， 但是， 针对诸如蛋白质核酸等含有成千上万个原子的大分子体系， 这几种方法的计算成本仍旧很高.

为了模拟研究蛋白质等大分子体系的结构和性能， 建立和发展了固定电荷分子力场方法^［25~29］， 该力场将原子分数电荷作为固定参数. 因为分子中各原子上电荷数值的大小不能随环境的变化而改变， 所以固定电荷力场无法描述体系中存在的极化作用. 随着研究的深入， 人们发现极化作用对于正确描述生物分子水体系、 电解质溶液和离子液体等体系至关重要. 由此， 可极化分子力场方法应运而生^［30~34］， 通过在固定电荷分子力场基础上引入诱导点偶极^［30~32］、 Drude Oscillator^［33］或浮动电荷等手段^［34］对体系中的极化效应进行处理. 相对于固定电荷分子力场， 可极化分子力场方法极大地提高了预测的准确性， 同时也大幅度降低了计算效率. AMOEBA可极化分子力场在原子点电荷的基础上， 又为每个原子引入3个偶极矩和6个四极矩^［30］. 这显著增加了AMOEBA可极化分子力场需要计算的静电作用项数目， 从而导致其计算效率大幅降低， 也极大地限制了可极化分子力场的应用.

近年来， 机器学习力场迅速发展， 成为传统分子力场之外的新兴模拟手段. 通过深度神经网络、 图神经网络、 高斯过程回归等算法， 使其能够对高精度量子化学数据进行拟合， 在分子结构优化、 势能面重构、 构象采样等方面表现出极高的准确性^{［35，36］}. 在理想情况下， 这些模型可以在计算成本远低于高阶量子化学方法的同时， 实现密度泛函理论甚至耦合簇方法级别的预测精度. 然而， 机器学习力场在广泛应用中仍面临若干挑战. 首先， 其性能高度依赖于训练数据的质量与覆盖范围^［37］. 若训练集不能涵盖目标体系的构象空间或化学空间， 模型在预测未见样本时往往出现泛化能力下降的问题. 尤其在涉及分子间弱相互作用、 多体效应或远离平衡态的反应路径时， 模型可能产生非物理的结果. 此外， 构建高质量训练数据集通常依赖于高阶量子化学计算， 这对大分子或多构象体系而言计算代价极高， 从而限制了机器学习力场的数据获取与训练效率. 其次， 机器学习力场在物理解释性方面仍然存在不足^［38］. 与传统基于明确物理模型（如电荷、 电荷转移、 极化）的力场不同， 机器学习力场的黑箱性质使得其在能量分解、 结构驱动机制分析等任务中缺乏透明性， 这对材料设计和药物筛选等需要 “结构-性质关系”明确反馈的应用场景带来限制.

针对于蛋白质多肽和核酸碱基等生物分子， 本课题组^［39~41］基于化学键偶极建立了可极化偶极-偶极作用模型. 大量测试表明， 可极化键偶极模型能够较准确地计算多肽分子间、 碱基分子间、 多肽碱基分子间相互作用强度， 计算结果与CCSD（T）， MP2等高精度量子化学方法的计算结果相等^［39~41］. 可极化键偶极模型还能够较准确地计算纯水团簇以及小分子-水团簇体系中的多体极化作用强度^［42~45］. 为了进一步拓宽可极化偶极作用模型的应用范围， 本文通过引入成键作用项， 将可极化偶极作用模型发展为可极化分子势函数， 并用于快速计算模拟环肽分子构象稳定性. 从文献中选取若干环肽分子构象， 将环肽分子中的N—H和C=O等极性化学键视为键偶极， 使用化学键固有偶极-固有偶极作用描述体系中的静电作用， 使用化学键固有偶极-诱导偶极作用以及诱导偶极-诱导偶极作用描述极化作用. 通过拟合训练集分子的高精度量子化学方法计算得到的二体和三体作用能随分子间距离变化的势能曲线确定所需参数. 将本文方法和参数应用于环肽分子构象能的快速计算， 并与高精度量子力学DLPNO-MP2/aug-cc-pVTZ方法的计算结果进行比较， 以检验本文势函数的合理性和参数的可转移性.

1 理论模型和参数化

1.1　理论模型

从数据库中提取9个环肽分子共计33个不同构象作为研究对象（表S1， 见本文支持信息）. 这些环肽分子由甘氨酸、 丙氨酸、 脯氨酸、 苯丙氨酸、 亮氨酸、 苏氨酸、 色氨酸、 络氨酸、 缬氨酸等氨基酸残基构成， 含有C=O， N—H， C—O， O—H等极性化学键和C—C， C=C， C—H等非极性或弱极性化学键. 仅将这些环肽分子中的C=O， N—H， C—O， O—H和C ^α —H 5种极性化学键视为键偶极矩， 不考虑其它非极性或弱极性化学键， 使用偶极-偶极作用描述体系中的静电和极化作用. 可极化偶极作用模型下的环肽分子势函数（U）如下式：

\begin{array}{l}U=\sum _{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}}{k}_b{\left(b-b_{\mathrm{0}}\right)}^{\mathrm{2}}+\sum _{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}{k}_{\theta }{\left(\theta -{\theta }_{\mathrm{0}}\right)}^{\mathrm{2}}+\sum _{\mathrm{d}\mathrm{i}\mathrm{h}\mathrm{e}\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{l}}\frac{V_n}{\mathrm{2}}\left[\mathrm{1}+\mathrm{c}\mathrm{o}\mathrm{s}\left(n\varphi -{\gamma }_n\right)\right]+\\ \sum _{l\ne m}\frac{{\mu }_l{\mu }_m}{r_{lm}^{\mathrm{3}}}\left(\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }^{\text{'}}+\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}{\alpha }^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{s}\beta \right)+\sum _{i\ne j}\left(\frac{A_{ij}}{R_{ij}^{\mathrm{12}}}-\frac{B_{ij}}{R_{ij}^{\mathrm{6}}}\right)\end{array}（1）

式中： 右侧前3项分别是键长、 键角和二面角作用项； b（nm）， θ（°）和φ（°）分别表示键长、 键角和 二面角； b₀（nm）， θ₀（°）和γ_n （°）分别表示平衡键长、 平衡键角和相角； k_b ［kJ·mol^-1·nm^-2）］， k_θ （kJ·mol^-1·degree^-2）和V_n （4.18 kJ/mol）分别为键长伸缩振动力常数、 键角弯曲振动力常数和傅里叶展开系数. 键长b、 键角θ和二面角φ均易由体系的笛卡尔坐标计算得到. 与传统分子力场方法相同， 本文将平衡键长b₀、 平衡键角θ₀、 相角γ_n 以及键长伸缩振动力常数k_b 、 键角弯曲振动力常数k_θ 、 傅里叶展开系数V_n 作为参数处理.

式（1）中右侧最后一项是范德华作用项， 其中， R_ij （nm）为第i个原子和第j个原子间的距离， A_{ij}=\sqrt[]{{\epsilon }_i{\epsilon }_j}{\left(R_i^{\mathrm{*}}+R_j^{\mathrm{*}}\right)}^{\mathrm{12}}， B_{ij}=\mathrm{2}\sqrt[]{{\epsilon }_i{\epsilon }_j}{\left(R_i^{\mathrm{*}}+R_j^{\mathrm{*}}\right)}^{\mathrm{6}}. R_i^{\mathrm{*}}（nm）和R_j^{\mathrm{*}}（nm）分别为第i个和第j个原子的范德华半径参数； {\epsilon }_i（kJ/mol）和{\epsilon }_j（kJ/mol）分别为第i个和第j个原子的范德华势阱参数.

式（1）中右侧第4项是可极化偶极作用模型框架下的静电作用项（包含静电极化作用）， 其中， μ_l （C·m）和μ_m （C·m）表示体系中的第l和第m个化学键偶极， r_lm （nm）， α（°）， α′（°）， β（°）是描述μ_l 和μ_m 两个键偶极矩空间相对关系的几何参量. 图1给出了这些几何参量的具体定义， A原子和B原子形成化学键AB， C原子和D原子形成化学键CD. E点为原子A和原子B连线的中点， F点为原子C和原子D连线的中点. 化学键A—B为极性化学键， 将之视为键偶极， 化学键偶极位于原子A和原子B连线的中点E， 其偶极矩大小使用符号μ_l 表示， 其方向由电负性大的原子A指向电负性小的原子B. 化学键C—D也为极性化学键， 化学键偶极位于原子C和原子D连线的中点F， 其偶极矩大小使用符号μ_m 表示， 其方向由电负性大的原子C指向电负性小的原子D. 如图1所示， r_lm 表示键偶极矩μ_l 和键偶极矩μ_m 间的距离（即点E和点F间的距离）， α和α′是键偶极矩正方向与线段EF间的夹角， β是平面BEF与平面EFD间的二面角. 通过A， B， C和D 4个原子的笛卡尔坐标， 容易计算得到几何参量r_lm， α，α′， β的值.

体系中任意一个化学键偶极μ（μ_l 或μ_m ）可进一步表示如下：

μ=μ₀+δμ（2）

式中： μ₀（C·m）为固有键偶极矩； δμ（C·m）为诱导键偶极矩.

在可极化偶极作用模型框架下， μ₀的大小不随环境的改变或分子构象的改变而改变， 可作为常数处理.δμ的大小可以随环境的改变或分子构象的改变而改变， 因而本文模型式（1）可以描述体系中的极化作用.

体系中任意一个化学键的诱导键偶极矩δμ的数值由下式计算：

δμ=c·(q-q₀)·d（3）

式中： d（nm）为化学键键长； q（e）为原子电荷； q₀（e）为参考电荷； c为校正系数. 如在本文理论模型框架下， 将环肽分子中的N—H， C=O等极性化学键视为键偶极. 对于极性化学键N—H， d为N—H化学键键长， q为H原子电荷， q₀为H原子的参考电荷； 对于极性化学键C=O， d为C=O化学键键长， q为O原子电荷， q₀为O原子的参考电荷. 化学键键长容易通过分子的笛卡尔坐标计算得到. 原子电荷q通过半经验AM1方法计算获得. 参考电荷q₀和校正系数c作为参数处理.

1.2　原子类型和参数化

使用本文模型式（1）计算任意一个环肽分子结构的能量， 需要首先确定平衡键长b₀、 平衡键角θ₀、 相角γ_n、 键长伸缩振动力常数k_b 、 键角弯曲振动力常数k_θ 、 傅里叶展开系数V_n 、 范德华半径、 范德华势阱、 固有键偶极矩μ₀、 参考电荷q₀和校正系数c等11种参数. 本文所研究的环肽分子含有脯氨酸（PRO）、 甘氨酸（GLY）、 丙氨酸（ALA）、 缬氨酸（VAL）、 亮氨酸（LEU）、 苯丙氨酸（PHE）、 苏氨酸（THR）、 色氨酸（TRP）和络氨酸（TYR）等9种氨基酸残基. 图2给出了这9种氨基酸残基的原子类型.

由图2可见， 在可极化偶极作用模型框架下， 除PRO， GLY， ALA， VAL， LEU， PHE， THR， TRP和TYR等每个氨基酸残基均含有N， H， C， O， CT和H1等6个原子类型（图2中第一行第一个结构式， 其中， R表示氨基酸侧链）和N—H， C=O， CT—H1等3类极性化学键偶极（此处碳原子CT为α碳原子， 将与α碳原子相键连的氢原子用H1表示）. 研究表明， 多肽分子中的α碳原子和氢原子形成的碳氢化学键可与多肽分子中的C=O键形成一定强度的氢键作用^［46］， 所以本文模型将α碳原子和氢原子H1形成的CT—H1键视为化学键偶极）. 除了N， H， C， O， CT和H1这6个原子类型之外， 依据侧链不同， 每个氨基酸残基还会存在不同的原子类型和化学键类型. 如ALA的侧链只有CT和HC两种原子类型和CT—HC一种化学键； TYR的侧链有CT， HC， CA， HA， OH， HO等6种原子类型和CT—HC， CA—HA，CA—OH， OH—HO， CT—CT， CA—CA等6种化学键； TRP的侧链有CT， HC， CA， CB， CN， CW， C*， NA， HA， H4， HNA等11种原子类型和CT—HC， CA—HA， CW—H4， NA—HNA等多种化学键. 相较于其它氨基酸残基， PRO没有原子类型H和化学键N—H（图2中第一行第二个结构式）. 本文仅将主链上的N—H， C=O和C ^α —H化学键以及侧链上的CA—OH， OH—HO， CT—OH， NA—HNA化学键作为键偶极处理.

环肽分子主链上的N—H， C=O和C ^α —H化学键以及侧链上的CA—OH， OH—HO， CT—OH， NA—HNA化学键所对应的固有键偶极矩μ₀、 参考电荷q₀和校正系数c等静电作用项参数的确定过程如下： 如图S1（见本文支持信息）所示， 设计二肽分子和两个水分子形成的氢键三聚体， 固定其中两个分子的位置不变， 改变第3个分子的位置， 使用MP2/aug-cc-pVTZ方法计算得到三体作用能随分子间距离变化的势能曲线. 对N—H， C=O和C ^α —H化学键的固有键偶极矩使用实验值为其初始值^［47］， 水分子O—H键的固有偶极矩和诱导偶极矩相关参数取自本课题组的前期工作^［43］， 使用本文模型计算二肽-水三聚体的三体作用能势能曲线. 调整优化N—H， C=O和C ^α —H化学键的固有键偶极矩μ₀和校正系数c， 直至使本文模型能够较好地重复MP2/aug-cc-pVTZ方法的三体作用能势能曲线， 从而确定了静电作用项参数μ₀和c. 对于TRP侧链所涉及的NA—HNA极性化学键， 其固有键偶极矩μ₀和校正系数c的取值与N—H的相同. 对于THR和TYR侧链所涉及的CT—OH， CA—OH， OH—HO极性化学键， 其固有键偶极矩μ₀的值直接取自实验值^［47］. 使用式（3）计算诱导偶极矩所需的参考电荷q₀的确定方法如下： 使用AM1方法计算多个具有铺展结构的多肽分子的原子电荷分布， 对于C=O键， 取其O原子的AM1电荷的平均值作为参考电荷； 对于N—H键， 取其H原子的AM1电荷的平均值作为参考电荷； 对于 NA—HNA键， 取其HNA原子的AM1电荷的平均值作为参考电荷； 对于CT—H1键， 取其H1原子的AM1电荷的平均值作为参考电荷； 对于OH—HO键， 取其HO原子的AM1电荷的平均值作为参考电荷. 表1列出了本文确定的静电作用项参数. 可以看出， 极性化学键C=O所对应的参考电荷q₀为负值， 而其它化学键所对应的参考电荷q₀为正值， 这与本文模型使用O原子的电荷作为C=O键的参考电荷和H原子的电荷作为N—H， O—H， C ^α —H键的参考电荷相一致. 表1中CT—OH和CA—OH化学键所对应的参考电荷q₀和校正系数c为零， 表示本文模型对这两个化学键的极化不作考虑.

本文所需范德华参数的确定方法如下： 如图S2（见本文支持信息）所示， 设计由两个多肽分子形成的氢键二聚体和堆积二聚体， 计算其二体作用能随分子间距离变化的势能曲线. 所有计算均在DLPNO-MP2/aug-cc-pVTZ方法下进行， 并采用TightPNO设置， 以确保计算结果的准确性. 使用AMBER99sb力场中的范德华参数作为本文所需范德华参数的初始值， 使用本文模型计算得到二聚体的二体作用能随分子间距离变化的势能曲线. 调整优化范德华参数， 使本文模型能够较好地重现DLPNO-MP2/aug-cc-pVTZ方法的二体作用能势能曲线， 从而确定了本文模型所需要的范德华参数. 表S2（见本文支持信息）列出了本文确定的范德华作用项参数. 对于范德华参数， 表S2中黑体数字表示该范德华参数与AMBER99sb所用参数不同， 正体数字表示该范德华参数与AMBER99sb所用参数相同.

本文模型所需的所有键长、 键角、 二面角作用项参数（即平衡键长b₀、 平衡键角θ₀、 相角γ_n 、 键长伸缩振动力常数k_b 、 键角弯曲振动力常数k_θ 、 傅里叶展开系数V_n 等）均直接取自AMBER99sb^［25］. 键长、 键角、 二面角作用项的所有参数列于表S3~表S5（见本文支持信息）.

2 应 用

图3为本文研究的9个环肽分子结构示意图， 其中， 大写字母表示L型氨基酸残基、 小写字母表示D型氨基酸残基， 圆括号里第一个数字表示氨基酸残基个数、 第二个数字表示原子个数、 第三个数字表示构象数. 如环肽分子FlGLFG中的大写字母F， L和G分别表示L型苯丙氨酸、 L型亮氨酸和甘氨酸， 小写字母l则表示D型亮氨酸. 因此， 环肽分子FlGLFG由2个L型苯丙氨酸、 2个甘氨酸、 1个D型亮氨酸和1个L型亮氨酸共计6个氨基酸残基构成， 含有92个原子和2种不同构象. 同理， 环肽分子PYPV由2个L型脯氨酸残基、 1个L型络氨酸残基和1个L型缬氨酸残基共计4个氨基酸残基构成， 含有65个原子和10种不同构象， 其余类推.

2.1　构象能

图3所示的9个环肽分子共计含有33个不同构象结构， 表S6（见本文支持信息）给出了这33个构象结构的笛卡尔坐标. 使用高精度DLPNO-MP2/aug-cc-pVTZ方法计算得到了33个构象结构的总能量. 对同一个环肽分子的不同构象结构， 将DLPNO-MP2/aug-cc-pVTZ方法计算得到的最低能量选为参考能量（即将该能量设为零. 本文将参考能量所对应的构象称为参考构象）， 从而得到DLPNO-MP2/aug-cc-pVTZ水平下环肽分子不同构象的构象能. 使用本文方法和AMOEBAbio18极化力场方法计算这33个不同构象结构的总能量. 将DLPNO-MP2/aug-cc-pVTZ方法的参考构象所对应的本文方法能量及

AMOEBAbio18方法能量取为零， 得到本文方法和AMOEBAbio18方法的环肽分子不同构象的构象能. 表2列出DLPNO-MP2/aug-cc-pVTZ方法、 AMOEBAbio18方法和本文方法计算得到的这9个环肽分子33个不同构象的构象能（RE）及AMOEBAbio18方法和本文方法的构象能相对于DLPNO-MP2/aug-cc-pVTZ方法的构象能的偏差（Δ）. 表2中第一列FlGLFG-1表示环肽分子FlGLFG的第一个构象，FlGLFG-2表示环肽分子FlGLFG的第二个构象； aGPFaGPF-1表示环肽分子aGPFaGPF的第一个构象， aGPFaGPF-2表示环肽分子aGPFaGPF的第二个构象， aGPFaGPF-3表示环肽分子aGPFaGPF的第三个构象， aGPFaGPF-4表示环肽分子aGPFaGPF的第四个构象， 其余类推.

如表2所示， 所有环肽分子的第一个构象被选为参考构象， 其相对能量设为零， 其余构象的能量相对于该参考构象进行比较. 以aGPFaGPF分子为例， 该分子由8个氨基酸残基组成， 共包含4个不同构象， 其中参考构象（aGPFaGPF-1）的构象能设为零， 其余3个构象的相对构象能分别被计算和比较. DPLNO-MP2/aug-cc-pVTZ方法的计算结果表明， 相对于环肽分子aGPFaGPF的参考构象（aGPFaGPF-1）， 其它3个构象的构象能分别为10.96， 53.81和391.04 kJ/mol. AMOEBAbio18预测的构象能分别为47.15， 105.60和436.56 kJ/mol， 相对于DLPNO-MP2方法构象能的误差分别为36.19， 51.80和45.52 kJ/mol. 本文方法预测的构象能分别为29.66， 64.94和361.04 kJ/mol， 相对于DLPNO-MP2方法构象能的误差分别为18.70， 11.13和-30.00 kJ/mol. 虽然两种方法均能正确预测构象稳定性的相对顺序， 但本文方法在构象能的绝对数值上更接近基准值， 显示出更高的精度. 对于另一个分子PYPV（由4个氨基酸残基组成）， 共存在10种不同构象. 相对于DLPNO-MP2/aug-cc-pVTZ方法的结果， 本文方法的预测误差在-3.68~0.08 kJ/mol之间， 而AMOEBAbio18方法的误差范围为-3.89~4.77 kJ/mol， 表明两者均具备良好的预测能力. 对本文研究的9个环肽分子共计24个构象能， AMOEBAbio18的计算结果 相对于DLPNO-MP2/aug-cc-pVTZ结果的均方根偏差（RMSE）为19.12 kJ/mol， 最大绝对偏差（MAE） 为51.80 kJ/mol， 本文方法的计算结果相对于DLPNO-MP2/aug-cc-pVTZ结果的RMSE为12.97 kJ/mol， MAE为33.85 kJ/mol， 说明本文方法对环肽分子构象能的预测精度整体略优于AMOEBAbio18力场方法.

尽管本文方法在大多数环肽分子的构象能预测中表现出较高的精度， 但在个别体系中仍存在一定的偏差， 可能源于环肽分子主链构象的高度刚性. 由于主链呈闭环结构， 构象自由度受到显著限制， 微小的键长、 键角或二面角的变化即可引发较大的能量波动， 从而使构象能计算对结构扰动更加敏感. 以pFTFWF分子为例， 在其构象pFTFWF-2中， 极化和范德华作用的能量贡献较小， 构象能主要由固有静电与成键作用主导. 值得注意的是， 本文模型对成键项的能量贡献存在一定高估， 其数值高达40.38 kJ/mol（表S6）， 这可能是导致该构象能预测偏差较大的主要原因. 这一结果提示， 对于构象高度受限、 非键作用贡献较小而成键能主导的体系， 本文势函数中成键项的建模仍有优化空间. 未来工作可进一步改进成键项的物理描述， 以提升模型在此类高刚性分子中的适用性与泛化能力.

众所周知， 对如蛋白质核酸等含有成千上万原子的生物大分子体系， 获得其高精度量子化学计算方法的构象能十分困难. 生物大分子诸多不同构象的高精度量子化学计算方法的构象能已不再是已知条件， 所以如表2那样选择量子化学方法计算得到的能量最低的构象为参考构象就变得不再可行. 为了更加合理地评估本文模型的计算精度， 将同一个分子的每一个构象的能量均取为参考能量. 以环肽分子aGPFaGPF为例， 该分子包含4个构象（aGPFaGPF-1~aGPFaGPF-4）. 通过依次将每个构象设为能量零点， 即参考构象， 可计算其余3个构象的相对构象能. 如， 以aGPFaGPF-1为参考构象时， aGPFaGPF-2， aGPFaGPF-3和aGPFaGPF-4 的相对构象能列于表2； 若以aGPFaGPF-2为参考构象， 则 aGPFaGPF-1， aGPFaGPF-3和aGPFaGPF-4的构象能分别为-10.96， 42.84和380.07 kJ/mol. 通过这种方式， 可全面获得所有构象之间的相对能量关系. 对于aGPFaGPF分子， 总计得到12组构象能数据（参考构象能量不计入）. 同理对环肽分子PYPV， 可得到10×9=90个构象能. 推而广之， 对于任意一个分子， 可以得到构象能的个数为M×（M-1）， 其中， M为该分子的构象数. 对本文研究的9个环肽分子的33个构象， 依次将同一个分子的每一个构象的能量取为参考能量， 共计得到144个构象能. 图4（A）和（B）分别为使用AMOEBAbio18方法和本文方法计算得到的环肽分子的144个构象能与使用DLPNO-MP2/aug-cc-pVTZ方法计算得到的144个构象能的线性相关图.

由图4可见， 本文方法构象能与DLPNO-MP2/aug-cc-pVTZ方法构象能的线性相关系数（R²）为0.9784， RMSE为13.43 kJ/mol， MAE为48.70 kJ/mol. 而AMOEBAbio18可极化分子方法的构象能与DLPNO-MP2/aug-cc-pVTZ方法的构象能的R²为0.9682， RMSE为16.28 kJ/mol， MAE为61.59 kJ/mol. 说明本文模型对环肽分子构象能的预测结果略优于AMOEBAbio18方法.

2.2　结构优化、 频率计算和计算量比较

通过对环肽分子的结构优化和频率计算进一步对本文方法进行评估. 分别使用密度泛函理论B3LYP/6-31G（d，p）方法、 本文方法和AMOEBAbio18方法对环肽分子的结构进行优化， 并将优化结构与实验结构进行比较. 图5给出了部分环肽分子结构的优化结果， 其中， 黑色是实验结构， 红色是 使用B3LYP/6-31G（d，p）方法优化得到的结构， 绿色是使用本文方法优化得到的结构， 黄色是使用 AMOEBAbio18方法优化得到的结构. 括弧中3个数字依次为使用AMOEBAbio18方法、 本文方法和B3LYP/6-31G（d，p）方法优化得到的结构相对于实验结构的RMSE.

由图5可知， 对于环肽分子FPaFPa的能量最低构象FPaFPa-1， B3LYP/6-31G（d，p）方法、 本文方法和AMOEBAbio18方法优化得到的结构相对于实验结构的RMSE分别为0.12， 0.10和0.12 nm. 对于环肽分子lLGlLG的能量最低构象lLGlLG-1， B3LYP/6-31G（d，p）方法、 本文方法和AMOEBAbio18方法优化得到的结构相对于实验结构的RMSE分别为0.05， 0.03和0.04 nm. 以上结果说明本文方法、 AMOEBAbio18方法与密度泛函理论B3LYP/6-31G（d，p）方法优化得到的环肽分子结构均与实验结构符合得很好.

使用AMOEBAbio18方法和本文方法对环肽分子的振动频率进行了计算， 并将计算结果与密度泛函理论B3LYP/6-31G（d，p）方法的结果进行比较［研究表明， 使用B3LYP/6-31G（d，p）方法计算得到的振动频率与实验值相符^［48］］. 图6（A）和（B）分别给出了环肽分子FPaFPa-1和lLGlLG-1的振动频率计算结果. 环肽分子FPaFPa-1含有88个原子， 共计有258个振动模式. lLGlLG-1含有90个原子， 共计有264个振动模式. 由图6可以看出， 对于环肽分子FPaFPa-1， 相较于密度泛函理论B3LYP/6-31G（d，p）方法的振动频率， AMOEBAbio18方法和本文方法的RMSE分别为85.39和82.88 cm^-1. 对于环肽分子 lLGlLG-1， 相较于B3LYP/6-31G（d，p）方法的振动频率， AMOEBAbio18方法和本文方法的RMSE分别为73.71和75.51 cm^-1. 这些结果说明， 本文方法以及AMOEBAbio18方法均能够较好地重现密度泛函理论B3LYP方法的振动频率.

多体极化作用对于正确描述蛋白分子水体系、 电解质溶液和离子液体等体系的性质至关重要^［49］. 非极化力场虽然简单快捷， 但是由于其使用固定电荷模型， 所以无法描述体系中存在的多体极化作用. 可极化力场能够较好地描述多体极化作用， 但是其所需完成的计算量相对于非极化力场又增加得太显著. 这极大限制了可极化力场的应用. 寻求建立计算量较小的可极化力场方法十分必要.

本文构建的可极化分子势函数使用极性化学键间的偶极-偶极作用描述静电作用， 而可极化分子力场AMOEBA方法使用原子点电荷、 原子偶极矩、 原子四极矩描述静电作用. 因而相比于AMOEBA可极化分子力场方法， 本文模型需要计算的静电作用项更少. 以环肽分子LFGFLG为例， 本文模型仅将该环肽分子中的C=O， N—H和C ^α —H化学键视为键偶极， 所以LFGFLG分子共计20个键偶极， 存在20×19÷2=190个偶极-偶极作用. 扣除间隔两个及两个以下化学键的偶极-偶极作用， 本文模型仅需计算144对偶极-偶极相互作用. 而AMOEBAbio18方法在每个原子位点放1个电荷、 3个偶极矩和6个四极矩. LFGFLG共计有92个原子， 所以在AMOEBAbio18框架下该分子共计有92个原子点电荷、 92×3=276个原子偶极矩、 以及92×6=552个原子四极矩. 92个原子点电荷间存在92×91÷2=4186个电荷-电荷作用， 扣除间隔一个化学键和间隔两个化学键的原子间的电荷-电荷相互作用， AMOEBAbio18需要计算3924个电荷-电荷相互作用. 92个原子点电荷和276个原子偶极矩间有92×276=25392个 电荷-偶极相互作用， 扣除间隔一个化学键和间隔两个化学键的原子间的电荷-偶极相互作用， AMOEBAbio18需要计算25392－92×3－262×6=23544个电荷-偶极相互作用. 276个原子偶极矩间 有276×275÷2=37950个偶极-偶极相互作用， 扣除间隔一个化学键和间隔两个化学键的原子间的 偶极-偶极相互作用， AMOEBAbio18需要计算37950－92×3－262×9=35316个偶极-偶极相互作用. AMOEBAbio18还需要计算若干项四极-四极相互作用、 电荷-四极相互作用和偶极-四极相互作用.

对于环肽分子pFTFWF， 本文模型仅将分子中的C=O， N—H， C ^α —H， C—O， O—H化学键视为键偶极， 所以该分子共计有23个键偶极， 存在253个偶极-偶极作用. 扣除间隔两个及两个以下化学键的偶极-偶极作用， 本文模型仅需计算168对偶极-偶极相互作用. pFTFWF有112个原子， 所以在 AMOEBAbio18框架下， 该分子共计有112个点电荷、 112×3=336个偶极矩、 以及112×6=672个四极矩. 112个原子点电荷间存在6216个电荷-电荷作用， 扣除间隔一个化学键和间隔两个化学键的原子间的电荷-电荷相互作用， AMOEBAbio18需要计算5890个电荷-电荷相互作用. 112个点电荷和336个偶极矩间有112×336=37632个电荷-偶极相互作用. 扣除间隔一个化学键和间隔两个化学键的原子间的电荷-偶极相互作用， AMOEBAbio18需要计算37632－112×3－326×6=35340个电荷-偶极相互作用. 336个偶极矩间有336×335÷2=56280个偶极-偶极相互作用. 扣除间隔一个化学键和间隔两个化学键的原子间的偶极-偶极相互作用， AMOEBAbio18需要计算56280－112×3－326×9=53010个偶极-偶极相互作用.

由92个原子的环肽分子LFGFLG到112个原子的环肽分子pFTFWF， AMOEBAbio18方法需要计算的电荷-电荷相互作用由3924项增加到5890项， 电荷-偶极作用由23544项增加到35340项， 偶极-偶极作用由35316项增加到53010项. 而本文模型需要计算的静电作用（即偶极-偶极作用）由144项增加到168项. 这说明， 随着体系原子数目的增多， AMOEBAbio18方法需要计算的静电作用项数急剧增加（AMOEBAbio18的静电作用计算量增加50%）， 而本文模型需要计算的静电作用项数仅稍有增加（计算量增加16%）.

由以上讨论可以看出， 相较于著名的可极化分子力场AMOEBA方法， 本文模型显著降低了需要计算的静电作用项数. 因此， 当应用于含有成千上万个原子的复杂体系时， 本文模型将会极大地降低计算量， 提高可极化力场的计算效率， 拓展可极化力场的应用范围.

为了进一步评估本文方法的计算效率， 分别对9个环肽参考构象（即各分子中能量最低的构象）进行单点能计算， 比较了本文方法与AMOEBAbio18力场在相同硬件条件下的CPU时间. 测试环境为 Intel（R） Xeon（R）CPU E5-2650 v4@2.20GHz， 采用单线程运行（表S7， 见本文支持信息）. 统计结果显示， 在未经程序优化的前提下， 本文方法与使用Tinker 8.10.1程序计算AMOEBAbio18力场在大多数体系中的计算时间相当， 部分体系中甚至略优， 说明所提出的模型在保持较高精度的同时， 也具有良好的计算效率潜力. 这为未来面向大规模生物体系的高效模拟打下了基础.

3 结 论

基于化学键偶极构建了适用于环肽分子的可极化分子力场势函数， 并通过构象能预测、 结构优化及振动频率计算等测试验证了其准确性与适用性. 相比于AMOEBA等典型可极化力场， 本文方法在保持较高精度的同时， 显著减少了静电作用项的计算数量， 具有较好的计算效率. 目前， 该方法已通过自研程序实现， 支持单点能、 多体相互作用及结构优化计算， 尚未与主流分子动力学软件直接兼容， 未来将致力于接口开发与程序优化， 以实现其在复杂体系中的应用拓展. 本文的方法为构建具有物理基础的高效可极化力场提供了新思路， 并有望应用于新型环肽分子的理性设计与筛选.

支持信息见http： //www.cjcu.jlu.edu.cn/CN/10.7503/cjcu20250173.

参考文献

[1] Zorzi A.， Deyle K.， Heinis C.， Curr. Opin. Chem. Biol.， 2017， 38， 24—29

[2] Jing X.， Jin K.， Med. Res. Rev.， 2019， 40（2）， 1—58

[3] Ojeda P. G.， Cardoso M. H.， Franco O. L.， Drug Discov. Today， 2019， 24（11）， 2152—2161

[4] Saunders G. J.， Yudin A. K.， Chem. Sci.， 2022， 13（44）， 12942—12944

[5] Zhang H.， Chen S.， RSC Chem. Biol.， 2022， 3（1）， 18—31

[6] Jin K.， Future Med. Chem.， 2020， 12（19）， 1687—1690

[7] Ji X.， Nielsen A. L.， Heinis C.， Angew. Chem. Int. Ed.， 2023， 63（3）， e202308251

[8] Li T.， Jiang S.， Li T.， Xu H.， Zhang X.， Yan R.， Wu X.， Jin Y.， Wang Z.， J. Med. Chem.， 2024， 67（14）， 11789—11813

[9] Frecer V.， Ho B.， Ding J. L.， Antimicrob. Agents Chemother.， 2004， 48（9）， 3349—3357

[10] Santini B. L.， Zacharias M.， Front. Chem.， 2020， 8， Article573259

[11] Xue X.， Wang X.， Ye C.， Gao M.， Li P.， Yu K.， Chen G.， Chem. Res. Chinese Universities， 2024， 40（6）， 1298—1310

[12] Zhang J.， Lv L.， Zhu H.， Zhang Y.， Xu X.， Long L.， Fu W.， Chem. Res. Chinese Universities， 2024， 40（6）， 1201—1211

[13] Zhang B.， Li Y.， Shi W.， Wang T.， Zhang F.， Liu L.， Chem. Res. Chinese Universities， 2020， 36（5）， 733—747

[14] Yin J.， Wang J.， Niu R.， Ren S.， Wang D.， Chao J.， Chem. Res. Chinese Universities， 2020， 36（2）， 219—226

[15] Ehlert S.， Grimme S.， Hansen A.， J. Phys. Chem. A， 2022， 126（22）， 3521—3535

[16] Plett C.， Grimme S.， Hansen A.， J. Comput. Chem.， 2024， 45（7）， 419

[17] Stylianakis I.， Zervos N.， Lii J.， Pantazis D.A.， Kolocouris A.， J. Comput. Aided Mol. Des.， 2023， 37（12）， 607—656

[18] Prasad V. K.， AOtero⁃de⁃la⁃Roza A.， DiLabio G. A.， Sci. Data， 2019， 6， 180310

[19] Řezáč J.， Bím D.， Gutten O.， Rulíšek L.， J. Chem. Theory Comput.， 2018， 14（3）， 1254—1266

[20] Mladek A.， Banas P.， Jurecka P.， Otyepka M.， Zgarbova M.， Sponer J.， J. Chem. Theory Comput.， 2014， 10， 463—480

[21] Wang P.， Shu C.， Ye H.， Biczysko M.， J. Phys. Chem. A， 2021， 125（45）， 9826—9837

[22] Sameera W. M. C.， Pantazis D. A.， J. Chem. Theory Comput.， 2012， 8（8）， 2630—2645

[23] Robertson M. J.， Tirado⁃Rives J.， Jorgensen W. L.， J. Chem. Theory Comput.， 2015， 11（7）， 3499—3509

[24] Yan X. C.， Robertson M. J.， Tirado⁃Rives J.， Jorgensen W. L.， J. Phys. Chem. B， 2017， 121（27）， 6626—6636

[25] Hornak V.， Abel R.， Okur A.， Strockbine B.， Roitberg A.， Simmerling C.， Proteins， 2006， 65（3）， 712—725

[26] Mackerell A. D. Jr.， Feig M.， Brooks C. L. III， J. Comput. Chem.， 2004， 25（11）， 1400—1415

[27] Dodda L. S.， Cabeza D. V. I.， Tirado⁃Rives J.， Jorgensen W. L.， Nuc. Acids Res.， 2017， 45（1）， W331—W336

[28] Geng H.， Jiang F.， Wu Y. D.， J. Phys. Chem. Lett.， 2016， 7（10）， 1805—1810

[29] Chen J. N.， Jiang F.， Wu Y. D.， J. Chem. Theory Comput.， 2022， 18（10）， 6386—6395

[30] Ren P.， Ponder J. M.， J. Phys. Chem. B， 2003， 107（24）， 5933—5947

[31] Shi Y.， Xia Z.， Zhang J.， Best R.， Wu C.， Ponder J. M.， Ren P.， J. Chem. Theory Comput.， 2013， 9（9）， 4046—4063

[32] Zhang C.， Lu C.， Jing Z.， Wu C.， Piquemal J.， Ponder J. M.， Ren P.， J. Chem. Theory Comput.， 2018， 14（4）， 2084—2108

[33] Lopes P. E. M.， Huang J.， Shim J.， Luo Y.， Li H.， Roux B.， MacKerell A. D. Jr.， J. Chem. Theory Comput.， 2013， 9， 5430—5449

[34] Liu C.， Li Y.， Han B. Y.， Gong L. D.， Lu L. N.， Yang Z. Z.， Zhao D. X.， J. Chem. Theory Comput.， 2017， 13（5）， 2098—2111

[35] Zheng T.， Wang A.， Han X.， Xia Y.， Xu X.， Zhan J.， Liu Y.， Chen Y.， Wang Z.， Wu X.， Gong S.， Yan W.， Chem. Sci.， 2025， 16， 2730—2740

[36] Tang H.， Xiao B.， He W.， Subasic P.， Harutyunyan A.， Wang Y.， Liu F.， Xu H.， Li J.， Nat. Comput. Sci.， 2025， 5， 144—154

[37] Wang Z.， Wu H.， Sun L.， He X.， Liu Z.， Shao B.， Wang T.， Liu T.， J. Chem. Phys.， 2023， 159， 035101

[38] Hu J.， Xu L.， Jiang J.， Sci. Sin. Chim.， 2025， 55（6）， 1734—1750

[39] Sun C. L.， Jiang X. N.， Wang C. S.， J. Comput. Chem.， 2009， 30， 2567—2575

[40] Li S. S.， Huang C. Y.， Hao J. J.， Wang C. S.， J. Comput. Chem.， 2014， 35（6）， 415—426

[41] Gao X. C.， Hao Q.， Wang C. S.， J. Chem. Theory Comput.， 2017， 13（6）， 2730—2741

[42] Li X. L.， Sun Y. J.， Tang Y.， Wang C. S.， Chem. J. Chinese Universities， 2021， 42（12）， 3664—3671

[43] 李晓蕾， 孙云娇， 唐颖， 王长生. 高等学校化学学报， 2021， 42（12）， 3664—3671

[44] Li X. L.， Li C. M.， Zhu J. Y.， Zhou Z.， Hao Q.， Wang C. S.， J. Comput. Chem.， 2023， 44（5）， 677—686

[45] Bai M. Y.， Gao S. S.， Wang L.， Jiang X. Y.， Wang Z. M.， Hao Q.， Gong L. D.， Wang C. S.， Chem. Phys. Lett.， 2025， 876， 142230

[46] Zhu J. Y.， Jiang X. N.， Zheng X. H.， Hao Q.， Wang C. S.， Chem. J. Chinese Universities， 2025， 46（4）， 20240518

[47] 祝佳怡， 姜笑楠， 郑笑函， 郝强， 王长生. 高等学校化学学报， 2025， 46（4）， 20240518

[48] Vargas R.， Garza J.， Dixon D. A.， Hay B. P.， J. Am. Chem. Soc.， 2000， 122（19）， 4750—4755

[49] Dean J. A.， Lange’s Handbook of Chemistry， 15th Edition， McGraw⁃Hill， Inc.， New York， 1998

[50] Szczepaniak K.， Szczesniak M. M.， Person W. B.， J. Phys. Chem. A， 2000， 104（16）， 3852—3863

[51] Bedrov D.， Piquemal J. P.， Borodin O.， MacKerell A. D. Jr.， Roux B.， Schroder C.， Chem. Rev.， 2019， 119（13）， 7940—7995

基金资助

国家自然科学基金(21773102)

辽宁省教育厅科学研究项目(LJKMZ20221411)