F–水合过程中共价相互作用的量子化学拓扑分析以及NPL4锌指结构的应用：用于界定量子力学/分子力学（QM/MM）方法的边界

《ACS Omega》：Quantum Chemical Topology Analysis of Covalent Interactions in the Hydration of F– along with the Zinc Finger of NPL4 and Its Application to the Delimitation of QM/MM Boundaries

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

金属蛋白的精确建模仍然是计算化学中的一个核心挑战，尤其是在经典力场无法描述具有显著共价特性的金属-配体相互作用时。在此，我们提出了一种诊断方法，用于在混合量子力学/分子力学（QM/MM）模拟中划分量子力学区域，以F-离子的水合作用和核蛋白定位4（NPL4）的锌指（ZF）结构域作为案例研究。我们的方法基于量子化学拓扑学中的波函数分析技术。具体来说，我们依赖于分子中原子的量子理论以及相互作用量子原子方法，将F-阴离子与溶剂水分子之间的相互作用能量，以及NPL4中Zn2+与其周围氨基酸之间的相互作用能量，划分为经典（离子）和交换-相关（共价）两部分。所提出的方法还提供了一个合理的标准，以确定是否将水分子包含在量子力学区域内，因为电子结构几何优化表明Zn2+···OH2相互作用可能是相关的。我们的诊断方法提供了一种系统的方法来定义合适的QM/MM边界，并且在共价贡献不确定或仅部分存在时，有助于决定是否在混合模拟的量子力学区域内包含溶剂或其他物种。分析显示，F-离子的第一层溶剂化层和Zn2+的第一层配位壳层在与F-和Zn2+离子的相互作用中表现出显著的共价性，这证明了NPL4中ZF的量子力学区域应至少包括三个半胱氨酸、一个组氨酸和金属中心。对于F-阴离子的水合作用，QM/MM混合模拟与经典分子力学在F-···O峰的配位数和径向分布函数的最大值方面存在重要差异。使用基于共价性的NPL4 ZF区域的混合QM/MM模拟能够恢复准确的配位几何结构和键长，纠正了经典力场引入的错误，包括虚假的水合作用和Zn2+中心的超配位。扩大或缩小量子力学区域都会导致不正确的配位几何结构和不太准确的Zn2+–配体距离分布。本文提出的诊断方法表明，基于化学键合情景的指标可以为离子和金属蛋白的溶剂化构建可转移的量子力学区域提供信息。总体而言，我们的结果建立了一种将波函数分析与混合建模相结合的通用协议，并突显了此类分析提高QM/MM模拟准确性的潜力。

1. 引言
分子动力学（MD）模拟已被广泛用于研究金属蛋白，(1?7) 然而经典力场往往无法充分描述金属-配体配位，特别是对于过渡金属。(1?3,8,9) 这些限制主要是由于固定电荷模型和球对称相互作用势不能考虑所研究金属中心与其配体之间的配位键的方向性或共价特性。(10,11) 混合量子力学/分子力学（QM/MM）计算通常比经典MD模拟提供了更准确的金属配位键研究方法。例如，QM/MM计算改进了对于Cd2+、Zn2+、Cu2+和Hg2+等过渡金属的配位环境的描述。(12?15) 另一个关于QM/MM计算相对于经典分子模拟更准确的情况是，尽管带有三体校正的经典MD模拟可以提供水合过渡金属的结构数据的适当评估，但它们不足以描述这些系统的某些特性。(12?15) 例如，正确描述水合能和离子-配体振动需要量子力学校正。(13,14) 考虑到Jahn–Teller畸变等效应在水合Cu2+中的作用也需要量子力学处理。(12) 此外，与MM模拟不同，QM/MM计算在建模和热力学研究方面非常有效。例如，QM/MM模拟对Zn2+阳离子的水合作用描述与实验确定的结构和热力学数据非常吻合。(14) 此外，与经典模型相比，QM/MM计算更准确地预测了相对溶剂化自由能。(14) 最后一个例子是，QM/MM再现了实验结果，即Zn2+的溶剂化自由能明显低于Mg2+，尽管这两种阳离子具有相同的净电荷和相似的溶剂结构。(14) 使用QM/MM模拟时，通常仅包含第一层配位壳层就足以准确描述上述效应，尽管包含第二层配位壳层可以进一步改进这种描述。(12?14) 实际上，混合QM/MM模拟的可靠性和准确性在很大程度上取决于量子力学区域的定义方式。(16?19) 为了在混合模拟中划分量子力学（QM）和分子力学（MM）区域，一些方法利用某些性质的收敛行为来定义最佳的QM/MM边界，旨在平衡准确性和计算效率。虽然在某些情况下增加量子力学区域的大小可以获得更准确的结果，但随之而来的计算成本也会显著增加。因此，需要确定一个合理大小的量子力学区域，同时仍能捕捉到所研究系统的基本物理和化学特性。众所周知，只有当使用非常大的量子力学区域（约500–1000个原子）时，能量和结构特性才会趋近于渐近极限。(20?24) 这种缓慢收敛的主要因素之一是量子力学区域与周围蛋白质环境之间的电荷转移。(25) 因此，需要一种有效的策略来确定最佳的量子力学区域，同时深入了解只有大尺度量子力学模型才能捕捉到的关键相互作用。Karelina和Kulik提出了一种方法，(26) 他们结合了两种互补技术，即电荷转移分析（CSA）和Fukui位移分析（FSA）。CSA方法考虑了在两种不同情况下（即考虑完整酶（holo）和缺少活性位点的蛋白质（apo）时残基（RES）上的原子部分电荷变化，从而计算电荷转移：ΔqRES = qholoRES ? qapoRES。(1) 类似地，FSA基于在片段上对Fukui函数求和：fRES+ = ∑i∈RES qi(N+1) ? qi(N) 和 fRES? = ∑i∈RES qi(N) ? qi(N?1)。(2) 在两种不同情况下，即添加额外残基与否的情况下。尽管CSA和FSA具有有前景的特点，但它们也存在一些不便之处。首先，我们指出CSA和FSA方法的基础物理原理，即通过原子电荷变化和Fukui函数变化来确定量子力学区域，最终体现在电子密度的极化上。这种极化可以在一定程度上通过经典极化力场来解释。然而，这些方法没有考虑共价作用的强度。我们强调这一点的重要性，因为共价性是一种量子力学效应，而经典力场无法完全考虑这一点。即使在没有共价键形成或破坏的系统中，正确描述共价作用也很重要，因为适当描述由共价作用控制的区域内的电荷密度可以影响并纠正MM区域的描述。例如，正确描述Cu2+离子水合作用中的第一层溶剂化层（量子力学区域）可以改善对第二层溶剂化层（MM区域）的描述，而这在量子力学处理中并未明确包括。(12) CSA和FSA方法的另一个不便之处是它们依赖于电荷分布分析方法，如Mulliken提出的方法，这些方法可能强烈依赖于计算的技术方面（例如基组的选择）。(27) 此外，为了效率和数值稳健性，CSA和FSA的实施基于Voronoi变形密度（VDD）部分电荷。(28) 虽然VDD电荷对基组的选择相对不敏感，(28) 但Voronoi单元仅是几何构造，没有任何理论或方法的支持。(29) 关于FSA，Fukui函数还有其他限制，(30) 这对于确定混合QM/MM模拟中的量子力学区域尤为重要。例如，FSA方法需要使用具有渐近正确交换的密度泛函近似。此外，离域误差可能导致Kohn–Sham分子轨道的不准确性，从而影响Fukui函数的计算，导致误判。(26) 对于这项研究来说，Fukui函数的另一个重要缺点是它不适合处理过渡金属复合物。(31) 混合QM/MM模拟中量子力学区域和分子力学区域之间的主要区别在于化学相互作用的处理方式。因此，利用最先进的技术来检查化学键合情景以确定这些区域似乎是合理的。更具体地说，由于共价性本质上是一种量子力学效应，我们考虑在混合QM/MM模拟的量子力学区域内对感兴趣的相互作用的共价贡献作为划分量子力学和分子力学区域的基础。与其它确定QM/MM混合模拟中QM子系统的方法（如CSA和FSA）相比，这项研究的主要特点是考虑了共价性。在这项工作中，我们开发并应用了一种基于量子化学拓扑学（QCT）领域的波函数分析方法来定义（i）F-阴离子的水合作用和（ii）核蛋白定位4（NPL4）的锌指（ZF）的量子力学区域，我们将这些系统作为案例研究。我们选择这些系统作为示例，以说明这里提出的方法。实际上，在这些案例研究中确定量子力学区域涉及选择（i）根据该系统的外部特性，包括F-离子水合和锌指的溶剂化水分子，以及（ii）围绕Zn2+中心的氨基酸。我们利用分子中原子的量子理论（QTAIM）来划分系统内的原子系统边界，因为从理论角度来看，QTAIM分区无疑具有最多的优点，(29) 例如，在计算原子动能的期望值方面。(32) 我们还使用了相互作用量子原子（IQA）对电子能量的分区来量化F-···H–OH和Zn2+···配体相互作用的共价特性。有了这些信息，我们划定了量子力学子系统，并使用GFN2-xTB方法(33)对量子力学区域进行了QM/MM模拟，对系统的其余部分使用标准力场进行了F-离子水合作用的分析和NPL4中锌指的研究。具体来说，我们实现了一个自适应的QM/MM方案，以动态考虑水溶液中F-阴离子的量子力学和分子力学子系统之间的水分子交换。对于NPL4中的ZF，我们基于QTAIM和IQA分析分析了两种代表性的量子力学区域选择方案。首先，(a) 通过仅包括与Zn2+相互作用具有强共价成分的残基来定义量子力学区域。该模型进一步用于进行QM/MM计算，通过扩展和缩小量子力学区域来评估相应模拟的结果。其次，(b) 使用(a)中定义的区域以及两个额外的水分子扩展量子力学系统，在Zn2+中心周围形成八面体配位。使用QCT确定的量子力学区域进行的混合模拟得到的配位几何结构和键长分布与电子结构理论和实验数据更为吻合。这些发现表明，通过电子密度的拓扑分析来检查化学键合情景可以作为定义混合QM/MM模拟中量子力学区域的严格、系统和明确的诊断工具，以研究离子和金属蛋白的水合作用。理论框架
我们通过密度泛函理论（DFT）计算研究了F–离子的水合过程以及NPL4的ZF（锌结合域）的电子密度ρ(r)，并使用了QTAIM方法。QTAIM利用这个标量场将电子系统的三维（3D）空间划分为原子盆地。这种方法能够将分子、分子簇和晶体固体分割成原子和官能团，从而恢复经典的化学概念，如化学键合和键级。(35) 另一方面，IQA是理论化学中一个有价值且多用途的工具，可以分析电子系统中的分子内和分子间相互作用，而无需参考系统。(36?41) IQA和QTAIM都被证明在以相同严格的基础上表征不同性质的化学相互作用方面非常重要。(40,42?47) 特别是，这些方法允许详细评估电子系统中原子间相互作用的离子和共价贡献，如下所述。
关于这些方法的一般特征，QTAIM已成为分子量子化学中的标准概念工具。(36) QTAIM定义了一个理论框架，该框架将分子片段分解为具有明确定义的量子力学可观测量期望值的量子子系统。例如，可以计算原子盆地ΩA内的净原子电荷，(48)
Q(ΩA)=Z(ΩA)??N(ΩA)???(ΩA)=??(ΩA)????(ΩA)?(4)
其中Z(ΩA)和?N(ΩA)?分别是盆地ΩA中的原子序数和预期的电子数。?N(ΩA)?的值是通过在相应盆地的体积上积分ρ(r)来计算的。QTAIM电荷的计算提供了关于电荷转移过程的有价值信息。(29,49) 此外，QTAIM还提供了关于化学相互作用性质的其他有价值见解。例如，QTAIM定义了两个原子之间的非局域电子数，称为非局域化指数（DI），(48)
DI(ΩA,ΩB)=?2cov(N(ΩA),N(ΩB))
DI(ΩA,ΩB)=?2cov(??(ΩA),??(ΩB))(5)
其中cov(x, y)表示随机变量x和y之间的协方差。这个指数代表了共价性对ΩA和ΩB之间相互作用的贡献；也就是说，它是共价性对相互作用影响的一个度量。QTAIM指标对于所研究系统的电子结构近似或基组的选择相对稳定。这些指标也是轨道不变的，因为它们是从电子密度的分析中得出的，而电子密度对应于狄拉克可观测量的期望值，不同于CSA和FSA。(48)
IQA分析通过引入将电子系统的总能量分解为分子内和分子间贡献的概念框架，扩展了QTAIM的空间划分方法，(34,50)
E=∑AEAnet+∑A∑A此外，iqa能量划分还允许我们通过将相互作用能量分为经典（vclab）和交换-相关（vxcab）项来描述分子内和分子间贡献，(8)
eabint=vabcl+vabxc(8)
vclab和vxcab项分别对应于原子ωa和ωb之间相互作用的离子和共价贡献。对于组g和h也有类似的公式。鉴于vxcab和di(ωa, ωb)都与共价键合有关，可以预期这两个量之间存在相关性，事实确实如此。更具体地说，vxcab的多极展开中的主要项直接涉及指数di(ωa, ωb)，即
vabxc≈?di(ωa,ωb)2r,
(对于较大的r)(9)
其中r是原子盆地ωa和ωb内核之间的距离。总之，iqa相互作用能量eintab可以分解为离子和共价成分，从而能够定量评估电子系统中每对相互作用中的离子性和共价性程度。
iqa分析最显著的计算需求来自于在qtaim原子边界定义的通常不规则体积上评估单电子和双电子性质所需的三维和六维积分。(37) 然而，方程9由于不涉及六维积分中的权重因子r12–1，因此显著节省了计算时间，并且di(ωa, ωb)的计算可以简化为在盆地ωa和ωb内的重叠积分计算。(37)

3. 计算细节
3.1. 氟离子的水合结构
由于x···(h2o)n簇的势能面上存在大量可访问的构型，确定局部最小值具有挑战性。为了高效地采样f–···(h2o)n复合物的这种景观，我们使用了alexandrova及其同事开发的梯度嵌入遗传算法（gega），并在cluster中实现了该算法。(54) 在gega中，遗传算法操作符生成新的候选结构，但每个候选结构随后都使用从头算梯度沿着势能面进行优化，以确保进化池中保留的每个结构都对应于一个真正的局部最小值。这一条件通过几何优化和随后的谐波频率分析来强制执行，从而过滤掉鞍点并确保在真正的局部最小值处比较异构体。这种梯度嵌入的工作流程减少了传统方案通常所需的遗传算法迭代次数，同时保持了结果几何结构和相对异构体能量的从头算质量。gega搜索和局部优化是使用orca 4.2和m06 2x/def2-tzvp近似方法进行的。(56,57) gega程序使用了10个候选结构的初始种群，并经过8代进化。然后使用b3lyp/def2-tzvp方法对低能量结构进行了重新优化，(57?59) 并进行了谐波频率计算，以确认最终几何结构代表相应势能面上的局部最小值。
3.2. npl4系统的模型构建和系统准备
在rcsb-pdb数据库中，我们选择了pdb id为6jwh的结构(60)（对应于酵母npl4的晶体结构)作为npl4 zf的模板，因为其分辨率很高（1.72 ?）。我们关注了酵母npl4蛋白的锌结合域（残基129–151），我们将其称为npl4(zn2+)129–151。(1) 使用pdb ids 6jwi (60)和6jwj (60)的额外模型作为模板来重建缺失的环，从而获得了用于经典和混合qm/mm模拟的完整且稳健的npl4 zf模型。结构完成使用了modeler v.10.6程序。(61) 使用charmm-gui (62?64)和h++ (65)服务器分配了质子化状态。在n端应用了?nh2基团，在c端应用了中性羧基。his139被建模为nδ1–h互变异构形式。天冬氨酸、谷氨酸、精氨酸和赖氨酸残基被视为带电的，而与锌结合的半胱氨酸被视为去质子的。所有其余的半胱氨酸残基都被建模为带电的。
3.3. npl4 zf的量子力学计算和口袋定义
基于npl4(zn2+)129–151的x射线结构（pdb id: 6jwh (60)）进行了几何优化。协调残基在cα位置被截断，并添加了甲基自由基作为封端基团。原子坐标使用b3lyp (58,59)交换-相关泛函和def2-tzvp (66)基组进行了优化。频率分析确认优化后的结构对应于相关势能面上的局部最小值。根据它们与zn2+中心的距离，将金属离子的协调环境分解为空间口袋。这些口袋（图s1）不应与第4.4节中定义的区域i–iii混淆。口袋1包含了直接与zn2+结合的残基，即cys137、his139、cys145和cys148。为了考虑溶剂水分子可能包含在qm区域中，我们从使用charmm36m力场的npl4的经典md模拟中提取了一个代表性快照。在这种结构中，金属中心表现出由三个半胱氨酸、一个组氨酸和两个水分子组成的八面体协调几何结构。使用半经验方法gfn2-xtb提取并优化了zf原子的坐标。(33) 关于这种方法的准确性，参考文献(67)报告了一项关于gfn1-xtb和gfn2-xtb在几何优化中的性能与混合meta-gga泛函tpssh-d3(bj)-atm以及def2-tzvpp基组相比的基准研究。使用这两种紧密结合方法计算的过渡金属复合物和有机金属超分子结构（直到hg (z = 80)）与使用dft计算的结果非常吻合。(67) 基于gfn2-xtb的几何优化，我们将口袋2定义为口袋1的扩展，包括相同的协调残基，并增加了两个与锌中心结合的水分子。口袋3包含了口袋1中的残基，以及另外两个相邻的残基，即met144和tyr147。最后，口袋4包括了口袋3中的氨基酸以及ser134和lys138。为了进一步评估上述gfn2-xtb几何优化的质量，我们将使用这种方法优化的结构与b3lyp/def2-tzvp优化的结构进行了对比。比较显示，b3lyp和gfn2-xtb对于zn–s距离的差异范围为0.03至0.07 ?，而对于zn–n距离，差异约为0.12 ?。两种结构的相似性还通过相应的均方根偏差（rmsd）0.20 ?以及系统的众多自由度得到了体现。使用b3lyp/tzvp (58,59,68)近似对每个口袋进行了单点dft计算，以计算电子密度。这些计算是使用gaussian v.16程序(69)进行的，随后使用aimall (70)程序套件对得到的电子密度进行了qtaim和iqa波函数分析。我们想指出，我们基于效率原因选择了gfn2-xtb方法，因为它在过渡金属复合物的几何结构方面表现良好。(67) 不幸的是，基于这种电子结构方法无法进行iqa能量划分，即使是近似方式也不可行。因此，我们使用了dft（更具体地说是b3lyp）密度来进行本研究中呈现的iqa分析。
3.4. 分子动力学模拟
使用li和merz开发的12-6-4 lennard-jones参数(71,72)以及tip3p水模型(73)进行了f–阴离子的水合经典分子动力学模拟。(73) 为了进一步研究f–阴离子的水合过程，采用了两种自适应qm/mm方案进行相应的模拟。第一种方案在qm区域使用了一个移动的水分子和四个固定的水分子，因此qm区域最多可以有五个h2o分子。第二种方案采用了类似的设置，但它扩展到了包括一个移动的六个固定的h2o分子，从而有效地探索了多达七个h2o单体的协调环境。对于npl4中的zf，对完整的npl4蛋白和截短的npl4(zn2+)129–151域进行了经典分子动力学模拟，而仅对后者进行了混合qm/mm分子动力学计算。这些模拟使用了三种力场：charmm36m、(74) ff99sb-ildn、(75) 和ff19sb。(76) 对于混合qm/mm模拟，qm区域的定义基于对(i) f–和h2o分子在水合过程中的相互作用以及(ii) zn2+中心及其周围配体（即第3.3节中定义的口袋1–4中的氨基酸和溶剂水分子）的化学键合情况的分析。特别是，我们关注了第4节中描述的这些相互作用的共价贡献。量子力学子系统使用gfn2-xtb (33)紧密结合方法进行了处理。 此外，iqa能量划分还允许我们通过将相互作用能量分为经典（vclab）和交换-相关（vxcab）项来描述分子内和分子间贡献，(8) eabint=VABcl+VABxc(8) vclab和vxcab项分别对应于原子ωa和ωb之间相互作用的离子和共价贡献。对于组g和h也有类似的公式。鉴于vxcab和di(ωa, ωb)都与共价键合有关，可以预期这两个量之间存在相关性，事实确实如此。更具体地说，vxcab的多极展开中的主要项直接涉及指数di(ωa, ωb)，即 vabxc≈?di(ωa,ωb)2r, (对于较大的r)(9) 其中r是原子盆地ωa和ωb内核之间的距离。总之，iqa相互作用能量eintab可以分解为离子和共价成分，从而能够定量评估电子系统中每对相互作用中的离子性和共价性程度。 iqa分析最显著的计算需求来自于在qtaim原子边界定义的通常不规则体积上评估单电子和双电子性质所需的三维和六维积分。(37) 然而，方程9由于不涉及六维积分中的权重因子r12–1，因此显著节省了计算时间，并且di(ωa, ωb)的计算可以简化为在盆地ωa和ωb内的重叠积分计算。(37) 3. 计算细节 3.1. 氟离子的水合结构 由于x···(h2o)n簇的势能面上存在大量可访问的构型，确定局部最小值具有挑战性。为了高效地采样f–···(h2o)n复合物的这种景观，我们使用了alexandrova及其同事开发的梯度嵌入遗传算法（gega），并在cluster中实现了该算法。(54) 在gega中，遗传算法操作符生成新的候选结构，但每个候选结构随后都使用从头算梯度沿着势能面进行优化，以确保进化池中保留的每个结构都对应于一个真正的局部最小值。这一条件通过几何优化和随后的谐波频率分析来强制执行，从而过滤掉鞍点并确保在真正的局部最小值处比较异构体。这种梯度嵌入的工作流程减少了传统方案通常所需的遗传算法迭代次数，同时保持了结果几何结构和相对异构体能量的从头算质量。gega搜索和局部优化是使用orca 4.2和m06 2x def2-tzvp近似方法进行的。(56,57) gega程序使用了10个候选结构的初始种群，并经过8代进化。然后使用b3lyp def2-tzvp方法对低能量结构进行了重新优化，(57?59) 并进行了谐波频率计算，以确认最终几何结构代表相应势能面上的局部最小值。 3.2. npl4系统的模型构建和系统准备 在rcsb-pdb数据库中，我们选择了pdb id为6jwh的结构(60)（对应于酵母npl4的晶体结构)作为npl4 zf的模板，因为其分辨率很高（1.72 ?）。我们关注了酵母npl4蛋白的锌结合域（残基129–151），我们将其称为npl4(zn2+)129–151。(1) 使用pdb ids 6jwi (60)和6jwj (60)的额外模型作为模板来重建缺失的环，从而获得了用于经典和混合qm mm模拟的完整且稳健的npl4 zf模型。结构完成使用了modeler v.10.6程序。(61) 使用charmm-gui (62?64)和h++ (65)服务器分配了质子化状态。在n端应用了?nh2基团，在c端应用了中性羧基。his139被建模为nδ1–h互变异构形式。天冬氨酸、谷氨酸、精氨酸和赖氨酸残基被视为带电的，而与锌结合的半胱氨酸被视为去质子的。所有其余的半胱氨酸残基都被建模为带电的。 3.3. npl4 zf的量子力学计算和口袋定义 基于npl4(zn2+)129–151的x射线结构（pdb id: 6jwh (60)）进行了几何优化。协调残基在cα位置被截断，并添加了甲基自由基作为封端基团。原子坐标使用b3lyp (58,59)交换-相关泛函和def2-tzvp (66)基组进行了优化。频率分析确认优化后的结构对应于相关势能面上的局部最小值。根据它们与zn2+中心的距离，将金属离子的协调环境分解为空间口袋。这些口袋（图s1）不应与第4.4节中定义的区域i–iii混淆。口袋1包含了直接与zn2+结合的残基，即cys137、his139、cys145和cys148。为了考虑溶剂水分子可能包含在qm区域中，我们从使用charmm36m力场的npl4的经典md模拟中提取了一个代表性快照。在这种结构中，金属中心表现出由三个半胱氨酸、一个组氨酸和两个水分子组成的八面体协调几何结构。使用半经验方法gfn2-xtb提取并优化了zf原子的坐标。(33) 关于这种方法的准确性，参考文献(67)报告了一项关于gfn1-xtb和gfn2-xtb在几何优化中的性能与混合meta-gga泛函tpssh-d3(bj)-atm以及def2-tzvpp基组相比的基准研究。使用这两种紧密结合方法计算的过渡金属复合物和有机金属超分子结构（直到hg (z=80)）与使用DFT计算的结果非常吻合。(67) 基于gfn2-xtb的几何优化，我们将口袋2定义为口袋1的扩展，包括相同的协调残基，并增加了两个与锌中心结合的水分子。口袋3包含了口袋1中的残基，以及另外两个相邻的残基，即met144和tyr147。最后，口袋4包括了口袋3中的氨基酸以及ser134和lys138。为了进一步评估上述gfn2-xtb几何优化的质量，我们将使用这种方法优化的结构与b3lyp def2-tzvp优化的结构进行了对比。比较显示，b3lyp和gfn2-xtb对于zn–s距离的差异范围为0.03至0.07 ?，而对于zn–n距离，差异约为0.12 ?。两种结构的相似性还通过相应的均方根偏差（rmsd）0.20 ?以及系统的众多自由度得到了体现。使用b3lyp tzvp (58,59,68)近似对每个口袋进行了单点dft计算，以计算电子密度。这些计算是使用gaussian v.16程序(69)进行的，随后使用aimall (70)程序套件对得到的电子密度进行了qtaim和iqa波函数分析。我们想指出，我们基于效率原因选择了gfn2-xtb方法，因为它在过渡金属复合物的几何结构方面表现良好。(67) 不幸的是，基于这种电子结构方法无法进行iqa能量划分，即使是近似方式也不可行。因此，我们使用了dft（更具体地说是b3lyp）密度来进行本研究中呈现的iqa分析。 3.4. 分子动力学模拟 使用li和merz开发的12-6-4 lennard-jones参数(71,72)以及tip3p水模型(73)进行了f–阴离子的水合经典分子动力学模拟。(73) 为了进一步研究f–阴离子的水合过程，采用了两种自适应qm mm方案进行相应的模拟。第一种方案在qm区域使用了一个移动的水分子和四个固定的水分子，因此qm区域最多可以有五个h2o分子。第二种方案采用了类似的设置，但它扩展到了包括一个移动的六个固定的h2o分子，从而有效地探索了多达七个h2o单体的协调环境。对于npl4中的zf，对完整的npl4蛋白和截短的npl4(zn2+)129–151域进行了经典分子动力学模拟，而仅对后者进行了混合qm mm分子动力学计算。这些模拟使用了三种力场：charmm36m、(74) ff99sb-ildn、(75) 和ff19sb。(76) 对于混合qm mm模拟，qm区域的定义基于对(i) f–和h2o分子在水合过程中的相互作用以及(ii) zn2+中心及其周围配体（即第3.3节中定义的口袋1–4中的氨基酸和溶剂水分子）的化学键合情况的分析。特别是，我们关注了第4节中描述的这些相互作用的共价贡献。量子力学子系统使用gfn2-xtb>此外，iqa能量划分还允许我们通过将相互作用能量分为经典（vclab）和交换-相关（vxcab）项来描述分子内和分子间贡献，(8)
eabint=vabcl+vabxc(8)
vclab和vxcab项分别对应于原子ωa和ωb之间相互作用的离子和共价贡献。对于组g和h也有类似的公式。鉴于vxcab和di(ωa, ωb)都与共价键合有关，可以预期这两个量之间存在相关性，事实确实如此。更具体地说，vxcab的多极展开中的主要项直接涉及指数di(ωa, ωb)，即
vabxc≈?di(ωa,ωb)2r,
(对于较大的r)(9)
其中r是原子盆地ωa和ωb内核之间的距离。总之，iqa相互作用能量eintab可以分解为离子和共价成分，从而能够定量评估电子系统中每对相互作用中的离子性和共价性程度。
iqa分析最显著的计算需求来自于在qtaim原子边界定义的通常不规则体积上评估单电子和双电子性质所需的三维和六维积分。(37) 然而，方程9由于不涉及六维积分中的权重因子r12–1，因此显著节省了计算时间，并且di(ωa, ωb)的计算可以简化为在盆地ωa和ωb内的重叠积分计算。(37)

3. 计算细节
3.1. 氟离子的水合结构
由于x···(h2o)n簇的势能面上存在大量可访问的构型，确定局部最小值具有挑战性。为了高效地采样f–···(h2o)n复合物的这种景观，我们使用了alexandrova及其同事开发的梯度嵌入遗传算法（gega），并在cluster中实现了该算法。(54) 在gega中，遗传算法操作符生成新的候选结构，但每个候选结构随后都使用从头算梯度沿着势能面进行优化，以确保进化池中保留的每个结构都对应于一个真正的局部最小值。这一条件通过几何优化和随后的谐波频率分析来强制执行，从而过滤掉鞍点并确保在真正的局部最小值处比较异构体。这种梯度嵌入的工作流程减少了传统方案通常所需的遗传算法迭代次数，同时保持了结果几何结构和相对异构体能量的从头算质量。gega搜索和局部优化是使用orca 4.2和m06 2x/def2-tzvp近似方法进行的。(56,57) gega程序使用了10个候选结构的初始种群，并经过8代进化。然后使用b3lyp/def2-tzvp方法对低能量结构进行了重新优化，(57?59) 并进行了谐波频率计算，以确认最终几何结构代表相应势能面上的局部最小值。
3.2. npl4系统的模型构建和系统准备
在rcsb-pdb数据库中，我们选择了pdb id为6jwh的结构(60)（对应于酵母npl4的晶体结构)作为npl4 zf的模板，因为其分辨率很高（1.72 ?）。我们关注了酵母npl4蛋白的锌结合域（残基129–151），我们将其称为npl4(zn2+)129–151。(1) 使用pdb ids 6jwi (60)和6jwj (60)的额外模型作为模板来重建缺失的环，从而获得了用于经典和混合qm/mm模拟的完整且稳健的npl4 zf模型。结构完成使用了modeler v.10.6程序。(61) 使用charmm-gui (62?64)和h++ (65)服务器分配了质子化状态。在n端应用了?nh2基团，在c端应用了中性羧基。his139被建模为nδ1–h互变异构形式。天冬氨酸、谷氨酸、精氨酸和赖氨酸残基被视为带电的，而与锌结合的半胱氨酸被视为去质子的。所有其余的半胱氨酸残基都被建模为带电的。
3.3. npl4 zf的量子力学计算和口袋定义
基于npl4(zn2+)129–151的x射线结构（pdb id: 6jwh (60)）进行了几何优化。协调残基在cα位置被截断，并添加了甲基自由基作为封端基团。原子坐标使用b3lyp (58,59)交换-相关泛函和def2-tzvp (66)基组进行了优化。频率分析确认优化后的结构对应于相关势能面上的局部最小值。根据它们与zn2+中心的距离，将金属离子的协调环境分解为空间口袋。这些口袋（图s1）不应与第4.4节中定义的区域i–iii混淆。口袋1包含了直接与zn2+结合的残基，即cys137、his139、cys145和cys148。为了考虑溶剂水分子可能包含在qm区域中，我们从使用charmm36m力场的npl4的经典md模拟中提取了一个代表性快照。在这种结构中，金属中心表现出由三个半胱氨酸、一个组氨酸和两个水分子组成的八面体协调几何结构。使用半经验方法gfn2-xtb提取并优化了zf原子的坐标。(33) 关于这种方法的准确性，参考文献(67)报告了一项关于gfn1-xtb和gfn2-xtb在几何优化中的性能与混合meta-gga泛函tpssh-d3(bj)-atm以及def2-tzvpp基组相比的基准研究。使用这两种紧密结合方法计算的过渡金属复合物和有机金属超分子结构（直到hg (z = 80)）与使用dft计算的结果非常吻合。(67) 基于gfn2-xtb的几何优化，我们将口袋2定义为口袋1的扩展，包括相同的协调残基，并增加了两个与锌中心结合的水分子。口袋3包含了口袋1中的残基，以及另外两个相邻的残基，即met144和tyr147。最后，口袋4包括了口袋3中的氨基酸以及ser134和lys138。为了进一步评估上述gfn2-xtb几何优化的质量，我们将使用这种方法优化的结构与b3lyp/def2-tzvp优化的结构进行了对比。比较显示，b3lyp和gfn2-xtb对于zn–s距离的差异范围为0.03至0.07 ?，而对于zn–n距离，差异约为0.12 ?。两种结构的相似性还通过相应的均方根偏差（rmsd）0.20 ?以及系统的众多自由度得到了体现。使用b3lyp/tzvp (58,59,68)近似对每个口袋进行了单点dft计算，以计算电子密度。这些计算是使用gaussian v.16程序(69)进行的，随后使用aimall (70)程序套件对得到的电子密度进行了qtaim和iqa波函数分析。我们想指出，我们基于效率原因选择了gfn2-xtb方法，因为它在过渡金属复合物的几何结构方面表现良好。(67) 不幸的是，基于这种电子结构方法无法进行iqa能量划分，即使是近似方式也不可行。因此，我们使用了dft（更具体地说是b3lyp）密度来进行本研究中呈现的iqa分析。
3.4. 分子动力学模拟
使用li和merz开发的12-6-4 lennard-jones参数(71,72)以及tip3p水模型(73)进行了f–阴离子的水合经典分子动力学模拟。(73) 为了进一步研究f–阴离子的水合过程，采用了两种自适应qm/mm方案进行相应的模拟。第一种方案在qm区域使用了一个移动的水分子和四个固定的水分子，因此qm区域最多可以有五个h2o分子。第二种方案采用了类似的设置，但它扩展到了包括一个移动的六个固定的h2o分子，从而有效地探索了多达七个h2o单体的协调环境。对于npl4中的zf，对完整的npl4蛋白和截短的npl4(zn2+)129–151域进行了经典分子动力学模拟，而仅对后者进行了混合qm/mm分子动力学计算。这些模拟使用了三种力场：charmm36m、(74) ff99sb-ildn、(75) 和ff19sb。(76) 对于混合qm/mm模拟，qm区域的定义基于对(i) f–和h2o分子在水合过程中的相互作用以及(ii) zn2+中心及其周围配体（即第3.3节中定义的口袋1–4中的氨基酸和溶剂水分子）的化学键合情况的分析。特别是，我们关注了第4节中描述的这些相互作用的共价贡献。量子力学子系统使用gfn2-xtb (33)紧密结合方法进行了处理。>相应地，MM区域也使用了各种力场进行了研究。所有的分子动力学（MD）模拟之前都进行了能量最小化步骤，并且使用AMBER 24软件进行。模拟中使用了边长为24.6 ?的立方盒子以及499个显式的TIP3P水分子，在周期性边界条件下模拟F–阴离子的水合过程。优化过程进行了5000步，其中前2500步采用最速下降法（steepest descent）。所有相互作用都设置了10.0 ?的非键合截断距离。另一方面，NPL4中的ZF在由显式TIP3P水分子组成的立方盒子中进行了溶剂化处理，确保任何蛋白质原子与盒子边缘之间的最小距离为10.0 ?。所有方向都应用了周期性边界条件。为了中和NPL4中ZF的净电荷并再现生理离子强度，系统中加入了Na+和Cl–离子。NPL4中ZF的能量最小化过程最多进行了10,000步：前1000步使用最速下降法，其余步骤采用共轭梯度最小化法。电静力和范德华相互作用都设置了12 ?的非键合截断距离，在能量最小化过程中对蛋白质主链施加了1.0 kcal mol–1 ?–2的位置约束。然后，在NVT系综下，使用Langevin动力学，以2 fs的时间步长，将NPL4中的ZF和水合后的F–离子逐渐加热100 ps，碰撞频率γ = 2.0 ps–1。在整个加热过程中，主链原子受到10.0 kcal mol–1 ?–2的力约束。温度从10 K线性升高到298 K。加热后，在298 K和1 atm的NPT条件下进行了20 ps的平衡过程，使用了Berendsen压强计（τp = 1.0 ps）和Langevin恒温器（γ = 2.0 ps–1）。在此阶段，主链原子上的力约束降低到1.0 kcal mol–1 ?–2。最后，对于NPL4中的ZF，在相同的恒温器和压强计设置下进行了无约束的生产MD模拟。关于F–阴离子的水合过程，模拟是在NVT系综下进行的。QM区域的电子结构使用GFN2-xTB近似方法计算，而系统的其余部分则使用Li和Merz参数（71,72）描述离子，以及TIP3P模型描述水分子（73）。为了确保QM区域的数值稳定性，使用了0.5 fs的减小时间步长，并且没有施加键合约束。通过Langevin恒温器维持碰撞频率γ = 5.0 ps–1。经典非键合相互作用使用10.0 ?的截断距离处理，而QM/MM界面上的电静力相互作用则使用粒子-网格Ewald方法处理，特定的QM截断距离为5.0 ?。经典和自适应QM/MM生产运行都进行了500 ps。

对于NPL4中ZF的所有混合QM/MM MD模拟，都采用了SHAKE算法来约束穿过QM/MM边界的键。非键合相互作用设置了10 ?的截断距离，QM/MM界面上的电静力接触使用粒子-网格Ewald方法处理。生产运行对于经典MD进行了10 ns，对于混合QM/MM模拟进行了5 ns，两种情况下的积分时间步长都是2 fs，温度保持在298 K恒定。使用SHAPE v2.1（80,81）和AMBERTOOLS v24.77软件进行了配位几何和金属-配体距离分布的分析。

4. 结果与讨论
4.1 F–阴离子的水合
首先，我们以氟化物阴离子的水合为例，验证本文提出的方法论。之后，我们使用NPL4中的ZF作为测试系统来评估我们的诊断方法。为了量化氟化物水合过程中涉及的共价贡献，我们利用IQA分区分析来评估F–···H2O相互作用中的共价和离子成分。如表1所示，[OH2···F···H2O]?簇中的交换-相关贡献最大，其值为?39.8 kcal mol–1，而经典贡献为?28.1 kcal mol–1。换句话说，F–···H2O相互作用中的共价成分比其离子成分强得多，这一点之前已有Sauza-de la Vega及其同事的报道（82）。我们强调，这一结果与化学直觉相反，从而展示了本文提出的诊断方法的一个重要应用领域。由于F–的电荷和H2O分子的大偶极矩，化学直觉会认为这种相互作用应该主要由强经典力（即库仑力）主导，例如电荷-偶极相互作用。更具体地说，F–与水分子之间的相互作用能量中，共价贡献分别为[OH2···F···H2O]?簇中的39.80 kcal mol–1和36.40 kcal mol–1。值得注意的是，这种交换-相关成分比(H2O)6中水分子之间的相互作用更强，如表1所示。

鉴于F–···H2O相互作用的明显共价特性，我们进行了QM/MM计算，在QM区域内包含了第一层溶剂化壳层的水分子，并将相应的结果与经典MD模拟进行了比较。我们通过F–···O和F–···H的径向分布函数（RDFs）及其相应的配位数来表征水合F–的结构特性。在这方面，经典MD模拟显示第一个F–···O峰位于2.60 ?（图1）。第一层水合壳层与第二层分离明显，平均配位数为6。在QM/MM模拟中，观察到一个更宽且不那么明显的第一个峰，其F–···O距离为2.69 ?，对应的第一层水合壳层的平均配位数为7（图2）。混合QM/MM计算中F–···O峰比经典MD计算中的峰更远，这与之前参考文献（82）中描述的F–···H2O相互作用的反合作用一致。与经典MD模拟相比，混合QM/MM模拟中F–的第一层水合壳层更大，这与Tongraar等人的先前结果（83）一致。这些研究者报告的QM/MM（Hartree–Fock (HF)/BJH-CF2）的第一峰距离约为2.68 ?，而经典MD的第一峰距离为2.60 ?。然而，他们报告的配位数（NC）为4.6，与我们在图2中报告的NC = 7形成鲜明对比。这种差异主要是由于HF理论中缺乏电子相关性。在没有色散力的情况下，短程Pauli排斥力没有得到吸引性范德华相互作用的平衡（84），导致了对离子的人为硬化描述，实际上将溶剂分子推开。相比之下，GFN2-xTB势提供了更可极化和稍微柔和的F–阴离子描述。当结合溶剂灵活性（O–H键拉伸）时，系统减少了空间障碍，允许第七个水分子在动态的第一层内稳定存在。这一特性通常被HF近似或经典模型的刚性所抑制。QM/MM结果进一步表明，第一层水合壳层中的水分子相对移动，促进了第一层与外部区域之间的更容易交换。在经典和混合QM/MM模拟中，都可以弱识别出第二层水分子（图3）。经典MD模拟显示第一个F–···H峰位于1.70 ?，随后是一个明显的第二个峰位于2.99 ?，对应于第一层水分子中的两个氢原子。第一个F–···H峰和F–···O峰之间的大约0.9 ?差异表明了线性氢键形成的偏好。另一方面，前两个F–···H峰移动到了更短的距离（1.61和2.84 ?），并且它们比经典MD模拟中的相应峰更宽且不那么明显。在混合QM/MM模拟中，第一个F–···H和F–···O峰之间的距离增加到1.08 ?，表明与经典MD计算中观察到的氢键线性不同。这种变化，加上观察到的O–H键拉伸和H–O–H角松弛，反映了更灵活和可变形的溶剂化几何结构。QM/MM（GFN2-xTB）描述捕捉到了更大的角度位移和更无序、更不稳定的水合复合物。GFN2-xTB势的非刚性特性有效地解释了离子与其灵活水环境之间的复杂相互作用，与HF预测的更刚性和无结构的壳层形成对比。

图1
(a) F–···O和(b) F–···H的径向分布函数及其相应的积分数。
高分辨率图像
下载MS PowerPoint幻灯片

图2
F–···H2O相互作用中F–···O RDF的第一最小值处的配位数分布。
高分辨率图像
下载MS PowerPoint幻灯片

图3
从自适应(a)混合QM/MM和(b)经典MD模拟获得的F–阴离子与其前七个配位水分子之间的距离分布。实线代表第三到第五近的水分子，而虚线黑色线条表示第六和第七近的分子。垂直虚线表示F–···O径向分布函数的第一最小值。虽然经典MD显示第七个分子明显倾向于位于第二层溶剂化壳层，导致CN = 6，但QM/MM结果显示分布更加分散，CN = 7。
高分辨率图像
下载MS PowerPoint幻灯片

我们采用了自适应QM/MM模拟，最初在QM区域内保持四个固定的水分子，以及一个在过渡区域内的移动水分子。在这个活跃区域内，基于距离的标准决定了一个分子是在QM级别还是MM级别处理。为了验证第一配位球内的额外分子是否是模拟设置的伪影，我们将QM区域扩展以包括一个移动的加上六个固定的H2O分子。这些扩展计算在F–···O和F–···H RDFs中产生了一致的结构结果，收敛到NC = 7。正如Hofer之前指出的（85），第一层配体分布的明显尾部可能会使得仅基于径向截断确定配位数变得困难。在这项工作中，我们分析了距离F–中心的七个最近水分子的距离分布。我们观察到连续的水分子邻居分布逐渐变宽和峰值强度减小。虽然最内层的分子显示出尖锐且狭窄的轮廓，但从第五个邻居开始，分布明显向更分散的曲线过渡。值得注意的是，第六和第七个邻居的大部分概率密度仍然位于F–···O RDF的第一最小值内。这些观察表明第一层水合壳层具有高度动态性，其中紧密结合的水分子核心被额外的H2O分子包围，这些分子虽然仍在配位范围内，但表现出更不稳定的特性，并且容易与其他水合壳层中的水分子交换（见图3）。图4显示了计算到F–···O RDF第一最小值的O···F–···O角分布。在经典势模拟中，可以清楚地识别出八面体排列，通过70–110°和150–180°之间的两个明显峰。QM/MM计算的···F–···O角显示第一个峰在50–100°之间，而在120–140°附近有一个更宽、不那么明显的峰。这些角分布的最大值比经典MD中的参考八面体配置要小，反映了配位数的增加，如前文所述。此外，混合计算导致第一层水分子配置更加无序，如图4所示。

总之，经典和含水QM/MM模拟对水合F–阴离子的不同结构描述强调了对该系统进行量子力学处理的必要性。虽然经典的分子动力学（依赖于固定点电荷和刚性的原子尺寸）强加了严格的八面体几何结构，但量子力学/分子力学（QM/MM）方法揭示了一个更加灵活的水合壳层，能够容纳第七个水分子。这种结构塑性与我们通过IQA分析发现的F–···H2O相互作用中的显著共价成分密切相关，这是经典势能本质上无法再现的特性。因此，包含显式的电子结构效应对于准确捕捉氟离子中共价性和水合动力学之间的微妙相互作用是不可或缺的。4.2. NPL4中的ZF关于NPL4中ZF的结果讨论分为四个部分，每个部分分别讨论了本文研究的金属蛋白的混合QM/MM模拟中量子力学处理区域的构建、验证和评估的不同阶段。首先，我们报告了通过QTAIM和IQA分析对NPL4蛋白中Zn2+配位所涉及的共价和离子贡献的表征。如前所述，这些分析构成了本文提出的用于界定QM/MM计算中QM区域的方法论基础。其次，我们讨论了使用不同QM区域的混合模拟性能。通过这些方法，我们能够表征基于共价性界定的QM区域，该区域与电子结构理论和结构实验数据相比提供了最佳结果。第三，我们将经典力场模拟的结果与混合QM/MM计算的结果进行了比较，并指出了经典MD计算在配位几何和金属-配体距离分布方面的系统性偏差。第四，我们评估了即使使用由经典MD得出的错误八面体几何结构进行初始化，混合模拟也能恢复正确的配位数，最后总结了我们的方法论工作流程，并讨论了该方法对其他金属蛋白的潜在通用性。4.3. QCT分析和QM区域的界定为了确定用于ZF的混合QM/MM模拟的化学上有意义的量子区域，我们应用了QCT理论框架来研究NPL4(Zn2+)129–151结构域中锌配位环境内的化学键合情况。我们使用QTAIM和IQA波函数分析方法将Zn2+与其周围氨基酸和溶剂化水分子之间的相互作用能量分解为经典（即离子）和交换-相关（即共价）贡献。为此，我们考虑了第3.3节中定义的四个口袋。如图5所示，两个能量分量在口袋1内都非常吸引人，该口袋包括金属中心以及Cys137、His139、Cys145和Cys148氨基酸残基。与共价特性相关的IQA交换-相关项对于Zn–S和Zn–N相互作用分别达到了-60到-80 kcal mol–1的值。虽然经典静电相互作用在这个内部区域也具有吸引力，但其范围更广，最小值约为-120 kcal mol–1。图5（左）IQA相互作用能量中的交换-相关部分和（右）经典部分，涉及NPL4蛋白ZF中口袋1–4（在第3.3节末尾定义）内的Zn2+与氨基酸残基以及溶剂化水分子之间的相互作用。高分辨率图像下载MS PowerPoint幻灯片在口袋2中，Zn2+···OH2相互作用的共价成分明显较弱（-15.82和-19.35 kcal mol–1），不到与His139观察到的最弱Zn–配体相互作用（-39.74 kcal mol–1）的一半。尽管Zn2+···OH2接触保留了部分共价特性，但其单独的贡献不足以稳定八面体配位环境。表1报告称，对于由六个水分子完全水合的Zn2+离子，每个Zn2+···OH2相互作用的平均共价贡献显著更强（高达-58 kcal mol–1）。此外，在环状水六聚体中，单个H2O···OH2相互作用可以贡献高达-27.53 kcal mol–1。换句话说，锌指环境中Zn2+···OH2相互作用的共价贡献在幅度上远小于表1中显示的所有其他相互作用，表明水与Zn2+的相互作用中的共价贡献远小于在水溶液中的贡献，甚至小于水簇中H2O分子的共价贡献。我们想强调的是，尽管水六聚体中的O···O距离（2.68 ?）比ZF中的Zn2+···O距离（2.47和2.50 ?）长，但ZF中的Zn2+···OH2相互作用仍然显示出明显的共价特性较弱。这些发现表明，在QM区域中包含水分子应该由定量相互作用能量而不是仅由几何标准来指导，因为这些指导原则可能会导致不同的结论。表1中的数据分析表明，在MD QM/MM模拟中应排除水分子。尽管希望可以预见排除显式溶剂化水分子的效果可以转移到其他锌指蛋白上，但我们认为这种情况不幸并非如此。如表1的第一行所示，Zn2+离子与水分子之间的相互作用可能表现出显著的共价成分。因此，我们认为必须针对每个单独的情况来确定其他锌指蛋白中是否包含水分子。当我们考虑那些不在口袋1和2中的氨基酸残基，即Lys138、Tyr147、Ser134和Met144时，我们注意到Zn2+中心与这些氨基酸之间的IQA交换-相关贡献迅速降至零，这表明这些物种之间的共价性可以忽略不计。通过考虑方程9中的VxcAB近似，我们在图S2中也观察到了同样的趋势。同样，Zn2+与口袋1内残基之间的IQA相互作用能量的经典贡献是吸引人的，而对于口袋3和4中的Lys138、Tyr147和Met144，同一项则变得排斥。这些结果表明，只有没有溶剂化水分子的Zn2+的第一配位壳层在其与金属中心的键合中表现出显著的共价特性，而与更远氨基酸残基的相互作用主要是静电的。基于这一分析，我们选择了口袋1中的原子作为NPL4中ZF的混合QM/MM模拟的合适量子力学子系统。该区域包括所有具有不可忽视的IQA交换-相关相互作用的残基，从而确保了键合的准确量子力学描述，同时通过排除可以经典处理的第二层氨基酸和溶剂化水分子来保持计算效率。因此，这个口袋被采纳为整个工作中分析NPL4中ZF的参考QM区域。在检查NPL4中ZF的混合QM/MM和经典MD模拟之间的比较之前，我们指出表1和图5中的数据表明，将一个部分包含在QM区域中的合理标准是（i）其与QM系统最相关部分的共价贡献与其静电贡献相似或更显著，即|VGHxc|?|VGHcl|，以及（ii）这种共价成分在能量上有重要贡献。表1和图5中的数据表明|VGHxc|?35kcalmol?1。4.4. 使用不同QM区域的混合模拟性能为了评估这种选择作为我们QM/MM计算中QM子系统的适用性，我们定义了三个候选量子区域，如图6所示。这些区域不应与第3.3节中定义的口袋1–4混淆。区域I是与口袋1相比的简化模型，它排除了His139残基，即区域I包括金属中心以及Cys137、Cys145和Cys148氨基酸。另一方面，区域II对应于基于上述Zn2+配位球共价性分析提出的最佳QM区域。换句话说，区域II包括Zn2+离子以及Cys137、His139、Cys145和Cys148。最后，区域III通过包括第一配位壳层之外的额外邻近残基来扩展区域II。也就是说，区域III包含了区域II中的物种以及Leu136、Met144、Glu146和Ser149。我们使用这三个区域进行了混合QM/MM模拟，以确定包含或排除特定氨基酸如何影响Zn2+周围的配位环境结构。第4.4节讨论了包含显式溶剂化水分子的QM区域的混合QM/MM模拟，以评估所选QM区域的适用性。图6为本次研究中考虑的NPL4锌指的混合QM/MM模拟定义的量子力学区域。（a）区域I包括Zn2+金属中心和三个配位半胱氨酸残基Cys137、Cys145和Cys148，但不包括His139。（b）通过QTAIM和IQA波函数分析方法检查化学键合情况确定的区域II包括区域I中的物种以及His139。（c）区域III通过包括第一配位壳层之外的额外邻近残基来扩展区域II，即Leu136、Met144、Glu146和Ser149。原子颜色表示代码：硫显示为黄色，氮显示为蓝色，氢显示为白色，碳显示为灰色，氧显示为红色。高分辨率图像下载MS PowerPoint幻灯片图7显示了图6中三个QM区域中Zn2+与其配位配体之间的距离分布。我们强调，只有使用区域II的模拟，该区域包括Zn2+金属中心以及通过QTAIM和IQA分析确定为与金属离子共价结合的Cys137、His139、Cys145和Cys148，才能得到与使用B3LYP/Def2-TZVP近似优化的几何结构相匹配的Zn–S和Zn–N相互作用距离分布。SHAPE度量使用一个无量纲参数（即连续形状度量（CShM）来量化一组原子与理想参考几何结构（例如IUPAC表示为T-4的四面体）的偏差。该参数评估了原子位置与理想形状位置、大小和方向最佳对齐后的总距离。S = 0.00的值表示与参考形状完全对应。S值较小（S < 1.0）对应于轻微的偏差，而S值超过1则表示与参考几何结构的显著偏离。（86）在这项工作中，我们获得了四面体几何结构的CShM得分为S = 0.81，而空的三棱双锥体、跷跷板和正方形几何结构的得分分别为S = 2.43、6.70和27.76。这些结果表明，QM/MM模拟中的Zn2+配位结构不仅再现了QM计算预测的原子间距离，还展示了与晶体学数据一致的四面体配位几何结构。（60）此外，三个硫原子之间的Zn–S距离在2.3 ?左右有明显的峰值（图7a–c），而Zn–N分布在一个接近2.0 ?的明确最大值附近（图7d），这与电子结构理论和实验确定的值一致。（60）图7图7. 从QM/MM模拟获得的NPL4蛋白ZF中Zn2+离子与其配位残基之间的距离分布，对于图6中定义的三个QM区域，使用CHARMM36m力场和GFN2-xTB半经验方法作为量子化学近似。电子结构计算（B3LYP/Def2-TZVP）的参考值用虚线绿色显示。面板（a）–（c）分别针对Cys137、Cys145和Cys148的Zn2+–S相互作用。面板（d）和His139的Zn2+–N相互作用也是如此。区域I中的Zn–N距离约为7.5 ?，这太长而无法维持Zn–N配位，因此根据将区域I视为该QM子系统，His139无法与Zn2+形成化学键。高分辨率图像下载MS PowerPoint幻灯片相反，排除His139的区域I未能正确描述Zn–N配位。该区域的Zn–N距离分布较宽，并向更长的距离偏移，如图7d所示，第一个峰值接近7.5 ?。经典力场无法准确描述配位相互作用，因为它们缺乏对电子转移、诱导极化以及金属-配体键的显著共价性质的明确描述。最重要的是，这一结果反映了在没有量子描述的情况下，分子力学（MM）方法无法稳定Zn–N相互作用。另一方面，Zn–S距离的描述相对较好，表明一个最小的量子力学（QM）区域就足以捕捉His139对Zn2+的硫配位。然而，同一个区域不足以维持氮对同一金属中心的正确配位。因此，区域I无法在NPL4中正确再现Zn2+离子周围的配位球。同样，考虑量子子系统区域III的QM/MM模拟（该区域包括额外的第二壳层残基）与使用电子结构理论（B3LYP/Def2-TZVP）计算的结构参数相比，并没有改善一致性。相反，区域III的Zn–S（Cys137）距离分布变得更宽，并且与区域I和II计算的结果相比有所偏离。这种结果质量的下降表明，将QM区域扩展到第一配位壳层之外可能会引入系统中的过度灵活性或由于QM/MM界面处电子相互作用不平衡而导致的不稳定边界效应。使用区域III作为量子力学子系统的QM/MM计算结果表明，如果没有系统和合理的选择标准来定义QM区域，整个计算容易出错。这一观察并不意味着扩大QM区域本身是有害的，而是表明不明确的边界定义可能会引入人为误差，无论QM区域的大小如何。实际上，如果没有适当的采样或精心定义的选择标准来扩展QM区域，可能会夸大MM区域中远处带电残基对QM子系统的影响。最后，仅基于距离标准选择QM原子通常是不足的，因为通常不使用化学直觉选择的残基可能对于能量和结构特性的正确收敛是必要的。为了进一步探讨这个问题，我们回顾了AMBER 24中默认实现的H-link原子方案，其中在QM/MM界面处放置一个额外的氢原子来封闭最终断裂的共价键。引入这样的辅助原子可能会导致几个众所周知的问题：(a) 每个链接原子引入了实际系统中不存在的结构自由度；(b) 会导致电荷不平衡；(c) 链接原子以及QM电子密度被放置在非常接近前沿MM原子（M）的位置，因此位于M上的点电荷可能会过度极化QM密度。为了缓解问题(a)，链接原子的坐标默认以父原子的坐标表示，这样就不会有额外的内部自由度显式地贡献到系统的能量或力中。关于问题(b)，即电荷不平衡，我们进行了采用电荷守恒方案的混合QM/MM计算。在这种方法中，当定义了QM区域时，QM原子和任何参与链接键的MM原子的原子电荷被设置为零，并且在整个模拟过程中明确约束和保持QM区域的净电荷。因此，QM区域在QM/MM计算过程中不能获得或失去净电荷，从而防止QM和MM子系统之间出现任何人为的电荷不平衡。然而，关于问题(c)，QM密度有被固定的MM点电荷过度极化的趋势。当点电荷与可极化的电子分布相互作用时，这种人为效应总是存在到某种程度，但当点电荷接近QM密度时，尤其是当QM密度可极化时，这种效应变得更加严重。当存在链接原子时，这个问题在边界处尤为关键，因为它们位于靠近前沿MM原子的位置。由于MM原子带有固定的部分电荷，靠近边界的带电氨基酸的出现可能会引入强烈的静电相互作用，这些相互作用会传播到QM区域。在我们的系统中，区域III包括一个位于QM/MM边界附近的Glu残基。因此，上述Zn–Cys137距离分布的缺陷行为可以归因于与该Glu–残基在界面处相关的强静电相互作用对QM区域几何的额外扰动。总的来说，这些结果表明区域II在计算成本和物理准确性之间提供了最佳平衡，它保持了正确的配位数和几何形状以及距离分布，同时避免了在QM区域中包含或排除残基所导致的人为误差。此外，这一分析表明基于IQA和QTAIM波函数分析方法得出的共价性评估来选择QM区域的有效性，并确立了区域II作为后续NPL4中ZF的混合QM/MM模拟最可靠的参考QM区域。为了评估经典分子力学在描述NPL4中Zn2+配位方面的性能，我们比较了使用三种标准力场（CHARMM36m、ff99SB-ILDN和ff19SB-GAFF）进行的经典MD模拟结果与使用区域II作为QM子系统的混合QM/MM模拟结果。我们在相同的NPL4(Zn2+)129–151域上进行了所有模拟，并关注了Zn–S和Zn–N距离分布以及Zn2+离子周围的几何配位。在测试的力场中，CHARMM36m与电子结构计算（B3LYP/Def2-TZVP）和实验参考数据的偏差最大。特别是如图8a所示，CHARMM36m模拟表现出将持续将两个水分子纳入Zn2+的配位球中的趋势，从而促进了Zn2+周围的假想八面体几何结构，而不是预期的四面体配置。这种效应在模拟的最初几纳秒后一致出现，并伴随着Zn–S键距离的系统性延长。从四面体到八面体的配位几何变化与Scrima等人的先前研究结果一致，他们利用CHARMM36m、ff99SB-ILDN和ff99SB*-ILDN力场研究了NPL4中的ZF。同样，参考文献(1)报告了使用CHARMM36m力场时，锌域和整个蛋白质中的配位数人为增加以及Zn–S和Zn–N键长度的高估。这些效应特别是在CHARMM36m（一个非键合模型）中，与Zn–O相互作用的过度稳定有关，这有利于添加氧化的额外配体（例如水），并改变了Zn2+金属中心的配位。CHARMM36m力场对Zn–O相互作用的过度稳定反映了其参数的制定，这在很大程度上依赖于对[Zn(H2O)6]2+复合物的考虑。表1报告称，[Zn(H2O)6]2+复合物中的Zn2+···OH2相互作用比上述NPL4中的ZF中的相互作用强得多。相比之下，使用ff99SB-ILDN力场得到的Zn2+周围有一个稳定的配位球。这一观察与之前使用ff99SB*-ILDN进行的NPL4中ZF的经典MD模拟结果一致，这些模拟显示配位键长度与参考QM距离的偏差很小（小于0.2 ?）。在这项研究中，ff99SB-ILDN和ff19SB与QM距离的偏差较大（分别约为0.5 ?和0.2 ?）。如图S3底部所示，AMBER力场倾向于产生狭窄的距离分布，这表明了一个深而刚性的势阱，限制了热波动。这样的结果可能表明模型过度拟合，迫使系统趋向于固定的、更短的平衡距离，可能会以牺牲结构真实性为代价。计算出的Zn–S距离始终比从晶体学数据获得的距离短（2.3 ?）。这一观察表明，力场更倾向于整体结构稳定性而不是局部化学准确性，可能导致Zn2+配位球的几何形状出现系统误差。总之，使用ff99SB-ILDN力场进行的模拟，在图8的(c)和(d)面板中显示的几何形状更接近使用区域II进行的混合QM/MM计算，避免了虚假的水合作用，与使用CHARMM36m进行的模拟相反。然而，ff99SB-ILDN力场倾向于低估Zn–配体键长度。虽然不如CHARMM36m的缺陷严重，但使用ff99SB-ILDN和ff19SB-GAFF计算的偏差表明，力场难以同时保持正确的配位几何和化学键指标。图8显示了使用CHARMM36m（a, b）和ff99SB-ILDN（c, d）力场进行的经典MD模拟得到的NPL4中ZF的Zn2+–配体相互作用距离分布。使用经典分子动力学计算得到的具有虚假水合八面体结构的参考混合GFN2-xTB/CHARMM36m模拟显示在(e)和(f)面板中。左侧面板显示了His139的Zn2+–N相互作用，而右侧面板分别讨论了Cys137、Cys145和Cys148的Zn2+–S接触。绿色虚线垂直线表示来自电子结构计算（B3LYP/Def2-TZVP）的参考距离。红色虚线垂直线和晶体学数据也是如此。CHARMM36m模拟显示出将水分子纳入Zn2+配位球的趋势，如图(a)所示。现在我们考虑Zn2+复合物的化学键合情景以及本研究中使用的力场。Zn2+离子是一个d10阳离子，具有完全填充的d壳层，表明电子分布是球对称的，并且倾向于近似对称的配位。尽管如此，配体原子的电负性和供体性质（即组氨酸的氮和半胱氨酸的硫）的差异会在锌中心周围的电子密度中引起强烈的极化效应。这些效应在IQA分析中得到了体现，其中Zn–His相互作用的共价部分的大小被发现为-60 kcal mol–1，而Zn–Cys相互作用的大小约为-80 kcal mol–1。这种共价贡献的不对称性加剧了Zn–配体距离的差异。尽管标准经典力场通常使用成对的、球对称的相互作用势来促进均匀的配位几何，但CHARMM36m能够再现Zn–配体距离的不对称性。然而，如前所述，这个力场还会在Zn2+的配位球中引入额外的水分子。相比之下，ff99-ILDN倾向于将Zn–N和Zn–S相互作用视为相似，从而在Zn2+离子周围产生更对称的配位环境。此外，QTAIM分析表明Zn2+离子带有大约1 au的部分电荷；因此，大量的电子密度从周围的配体捐赠给锌中心。总的来说，这些观察突出了使用标准经典力场难以捕捉的两个现象：由具有不同供体性质的配体引起的电子密度的各向异性重新分布以及由此产生的配位几何的不对称性。如前所述，AMBER力场再现了锌域的正确配位数，但显示了狭窄的距离分布。这种效应导致相应模拟中的动态受限，这可能是因为这些力场被参数化为再现固定的参考值。相比之下，QM/MM计算显示了更宽的距离分布，其峰值接近参考值。这些更宽的位移在化学反应中尤为重要，因为在化学反应中原子运动对于化学键的形成和断裂至关重要。虽然CHARMM36m力场可以通过其参数化策略近似Zn2+离子周围的不对称几何，但只有QM/MM混合计算才能准确描述Zn2+中心周围第一配位球的电子结构和几何形状。混合量子力学/分子力学（QM/MM）模拟，以经典分子动力学（MD）结构为初始条件

为了评估本文提出的NPL4中ZF（锌指结构）的QM区域的稳健性，我们检验了当模拟从使用经典MD模拟计算得到的畸变结构开始时，该区域是否能够恢复Zn2+中心周围的正确配位数（从6变为4）。具体来说，我们使用了经典MD CHARMM36m计算的最终快照作为混合QM/MM计算的起点，在这些快照中，有两个水分子进入了金属配位球，导致Zn2+周围形成了一个人为的八面体排列，并且另外两个水分子作为QM子系统被纳入计算。我们在图8的(e)和(f)面板中报告了相应的结果。在所有情况下，混合模拟在最初的几纳秒内都恢复到了四配位数。Zn–O距离分布显著扩大，配位水分子远离了金属中心，而Zn–S和Zn–N距离则恢复到了与DFT电子结构计算和实验数据一致的值。这种行为证实了混合方案能够纠正由力场引起的一些人为误差，选择性地稳定具有化学意义的相互作用，同时抑制虚假的配位事件。然而，配位几何结构从八面体变为一个空的三棱双锥体，其CShM得分为0.95，而四面体几何结构的得分为4.48。

为了进一步研究这个问题，我们使用GFN2-xTB近似方法对Region III和两个水分子进行了几何优化。我们观察到，当水分子靠近金属中心时，最小能量结构确实可能呈现八面体几何形状。因此，水分子的溶剂化过程确实可以扭曲金属离子周围的几何结构。然而，正如QM/MM混合模拟所显示的，水分子在系统演化过程中被逐出。简而言之，当溶剂化H2O分子出现在Zn2+中心附近时，上述从四面体几何结构的畸变在物理上是合理的，但这些分子最终会被迫离开金属中心的配位球。

此外，从人为构建的八面体配置开始的混合QM/MM模拟由于水分子接近QM/MM边界而表现出不稳定性。在整个模拟过程中，H2O单体逐渐移动到QM/MM边界区域，导致QM和水分子之间的交换。这种原子在穿过边界时的物理描述变化会在势能面上引入显著的不连续性，严重影响混合QM/MM模拟中量子力学计算的收敛性。为了减轻这种影响，我们专注于水分子与QM/MM边界保持足够距离的短时间段进行分析。在这些精心选择的时间内，混合模拟保持稳定，并一致地纠正了Zn2+中心的配位数。

这些发现表明，混合QM/MM模型不仅能够保持准确的金属-配体相互作用，还能够动态纠正结构畸变，例如由经典MD计算引入的配位数。然而，最终的配位几何结构仍然对初始配置敏感。虽然经典力场在金属-配体共价性强的系统中具有有限的适用性，常常导致错误的配位几何结构或虚假的配体结合，但使用Region II作为QM子系统的混合QM/MM模拟至少在某种程度上能够纠正Zn2+的结构和键合模式。这些观察结果强调了划定基于共价性信息的QM区域的重要性，并突出了在具有复杂配位化学的金属蛋白建模中结合电子结构理论和波函数分析数据的重要性。

5. 结论

本研究提出了一种系统化的协议，用于构建和验证涉及离子和水溶性金属蛋白的QM/MM模拟中的QM区域。我们的方法依赖于量子化学拓扑分析，特别是QTAIM和IQA理论框架，来量化不同化学相互作用的共价性。通过将这一协议应用于F-阴离子的水合过程，我们证明了F-···H2O相互作用具有重要的共价贡献。因此，在混合QM/MM模拟的QM区域中包含F-的第一层溶剂化分子，会导致配位数和径向分布函数的重要差异。同样，将该协议应用于NPL4蛋白的ZF结构域，使我们能够证明基于共价性贡献选定的QM区域在计算效率和化学准确性之间达到了最佳平衡。通过检查金属离子Zn2+的化学键合情况来确定QM区域的混合QM/MM模拟，能够保持距离分布，准确再现实验和电子结构计算得到的几何形状和结构参数。该区域排除了具有可忽略共价性的氨基酸残基和溶剂化水分子，证实了在仔细分析系统化学键合情况后，使用最小限度的QM描述是足够的。此外，具有我们QCT分析确定的QM区域的混合QM/MM模拟纠正了一些由经典力场引入的结构畸变，如配位数，这突显了传统MD在金属-配体共价性强系统中的局限性。本文提出的诊断方法遵循一个工作流程（共价性分析、QM区域划定以及通过混合QM/MM模拟进行验证），为模拟不同于NPL4的金属蛋白提供了一种可转移的策略。这种方法特别适用于涉及边界共价配位或难以表征的配体环境的情况。基于化学键合情况分析来划定QM区域，而不是通过任意距离截止或残基身份来确定相同区域，有助于提高混合QM/MM模拟的可靠性和可解释性。未来的应用可能支持在无机化学以及生物无机和药物化学中设计更准确的QM/MM模拟。