文章_可极化力场（PPC）

摘要

对溶液中的蛋白质进行了从头量子力学计算，以产生极化的特定蛋白质电荷（PPC），用于蛋白质的分子动力学（MD）模拟。通过开发基于片段的量子化学方法和隐式的连续溶剂模型，使蛋白质的量子计算成为可能。蛋白质的计算电子密度被用来推导PPCs，这些PPCs代表了蛋白质在原始结构附近的极化静电状态。这些PPCs与标准力场中的PPCs一样是以原子为中心的，因此在计算上对蛋白质的分子动力学模拟具有吸引力。为了研究电子极化对硫氧还蛋白的结构和动力学的影响，已经进行了广泛的MD模拟。我们的研究表明，通过明确比较使用PPC和AMBER指控的MD结果，硫氧还蛋白的动力学被电子极化所稳定。特别是，使用PPCs的MD自由能计算准确地再现了埋在硫氧还蛋白内的可电离残基Asp26的pKa移位的实验值，而以前使用标准力场的计算则高估了pKa移位的两倍之多。通过严格的MD自由能模拟准确预测埋藏在蛋白质内部的可电离残基的pKa移位，几十年来一直是计算生物学的一个重大挑战。这项研究提出了强有力的证据，证明蛋白质的电子极化在蛋白质动力学中起着重要作用。对溶液中的蛋白质进行了无源量子力学计算，以产生极化的蛋白质特定电荷（PPC），用于蛋白质的分子动力学（MD）刺激。通过开发基于片段的量子化学方法和隐含的连续溶剂模型，使蛋白质的量子计算成为可能。蛋白质的计算电子密度被用来推导PPCs，这些PPCs代表了蛋白质在原始结构附近的极化静电状态。这些PPCs与标准力场中的PPCs一样是以原子为中心的。为了研究电子极化对硫氧还蛋白的结构和动力学的影响，已经进行了广泛的MD模拟。我们的研究表明，通过明确比较使用PPC和AMBER电荷的MD结果，硫氧还蛋白的动力学被电子极化所稳定。特别是，使用PPCs的MD自由能计算准确地再现了埋在硫氧还蛋白内的可电离残基Asp26的pKa移位的实验值，而以前使用标准力场的计算则高估了pKa移位的两倍之多。通过严格的MD自由能模拟准确预测埋藏在蛋白质内部的可电离残基的pKa移位，几十年来一直是计算生物学的一个重大挑战。这项研究提出了强有力的证据，证明蛋白质的电子极化在蛋白质动力学中起着重要作用。

引言

蛋白质中相互作用的一个重要组成部分是静电相互作用。蛋白质的折叠、蛋白质与配体的结合、蛋白质与蛋白质的相互作用、电子转移、质子的结合和释放、酶的反应等过程主要是由静电相互作用驱动的。特别是，质子结合是一个依赖于pH值的过程，并受到与蛋白质和局部环境的静电作用的强烈影响。蛋白质侧链之间的静电作用不仅取决于它们的距离，而且还取决于它们在蛋白质中的位置和它们的局部溶剂环境。在质子结合中，蛋白质和溶剂经历了涉及电子极化和原子团位移的介电松弛。蛋白质内部的局部静电环境是不均匀的和疏水的，这通常对电离不太有利。因此，与蛋白质表面或溶液中的孤立氨基酸相比，蛋白质内部的可电离（或带电）残基可能有很大的pKa偏移。

目前的标准力场，如CHARMM和AMBER，在应用上取得了巨大的成功，但存在着根本的局限性。具体来说，这些FFs是氨基酸特异性的或类似于均值场的，因此不能准确表示高度不均匀和蛋白质特异性的特定蛋白质环境的静电作用。例如，在相同或不同的蛋白质中的两个相同类型的氨基酸，由于其不同的静电环境，应该有相当不同的电荷状态。而目前基于氨基酸的FFs无法描述特定蛋白质结构的极化状态，例如，原生蛋白质结构。为了克服标准FFs的这一基本缺陷，人们努力开发肽和蛋白质的可极化力场模型。

在这篇文章中，我们提出了一个包含极化的特定蛋白电荷(The polarized protein-specific charge) 的替代力场，用于蛋白质的动力学研究。极化特定蛋白质电荷（PPC）在原子电荷中建立了蛋白质的极化。PPC本身并不是 "可极化的"，但它是由第一原理的量子溶解计算得出的原生（或特定）结构中的蛋白质。因此，PPCs正确地代表了蛋白质的电子极化状态，因此提供了接近原生结构的准确的静电作用。在我们的方法中，蛋白质溶解的量子化学计算是通过将最近开发的基于片段的方案，即molecular fractionation with conjugate caps （MFCC），与连续溶剂模型相结合而实现的。在这篇文章中，我们采用线性化泊松-波尔兹曼方法来解决自洽反应场方程，同时使用MFCC方案对溶质进行量子化学计算 。使用RESP方法对蛋白质片段（氨基酸）的收敛电子密度进行拟合，生成蛋白质中每个氨基酸的部分电荷。拟合出的原子偏电荷是protein-specific的，它们正确地代表了蛋白质在原生（或其他给定）结构中的极化电子状态。由于PPC是以原子为中心的，并保持了与AMBER中标准电荷相同的简单性，它可以很容易地应用于MD模拟，而不需要任何额外的复杂化。因此，我们期望PPCs在蛋白质接近其原始结构的MD模拟中，在结构和动力学方面都能提供更好的静电相互作用。为了在MD模拟中应用PPC，人们只需用PPC取代AMBER力场中的标准电荷，同时保持其余的力参数不变。

为了研究PPCs在蛋白质的MD模拟中可能产生的影响，我们进行了一系列的MD研究，研究了硫氧还蛋白的各种结构和动力学特性。使用PPCs的新的MD结果与相应的使用了AMBER电荷MD研究中得到的结果进行了明确的比较。特别是，进行了MD自由能模拟以预测埋藏Asp26的硫氧还蛋白的pKa移动。以前使用AMBER和CHARMM进行的MD自由能计算给出的pKa偏移是实验值的两倍。为了与这些结果进行直接比较，我们按照Simonson等人的计算方法和程序来计算pKa移位，只是用我们计算的硫氧还蛋白PPC取代标准AMBER电荷。这将消除结果比较中任何可能的动态不确定性。

理论方法

在连续体-溶剂模型中，溶质(蛋白质)由电荷分布ρ(r)表示，它嵌入到一个空腔中，周围是介电常数为ϵ的可极化介质。溶质电荷分布ρ(r)使介电介质极化，并产生一个反应场，该反应场反作用于溶质，直到达到平衡。根据经典的静电理论，作用于溶质的反应场可以有效地由腔体表面的感应电荷来表示。通过离散化空腔表面的诱导电荷并迭代求解表面电荷产生的外部反应场中的溶质的量子化学方程，可以得到流行的PCM方法，最近被推广到蛋白质的溶解。然而，对于大的蛋白质，PCM模型需要许多离散的表面电荷，从而使线性方程的解决在计算上变得困难。在目前的方法中，我们通过数值求解泊松-波尔兹曼（PB）方程来获得反应场。这避免了大型线性方程的求解。

在我们的方法中，拟合蛋白质原子电荷的基本程序可以描述如下。首先，用MFCC方法进行蛋白质的气相计算，以获得给定结构的蛋白质的初始电子密度。使用RESP程序将计算出的电子密度用于拟合原子电荷。这里使用的电荷拟合理念与AMBER力场中使用的相同，这就保证了PPC与AMBER力场的其他参数一致。然后进行PB方程的求解，得到反应场，从中产生腔体表面的离散表面电荷。

蛋白质的MFCC/RESP电荷拟合方法与Poison-Boltzmann连续溶剂的耦合应提供蛋白质在水中的溶解自由能的准确估计。由MFCC/RESP程序产生的蛋白质每个原子上的Partial charges被传递给PB求解器DELPHI以确定自洽反应场。在电介质边界上得出一组诱导表面电荷q^ind，代表溶剂分子的反应场效应。介质溶质/溶剂边界由AMBER范德瓦尔斯半径定义，溶质分子原子的探测半径为1.4 Å。内部介电常数，表示为ϵsolute，被设置为unity，因为分子的极化性明确地包含在量子力学计算中。溶剂介电常数，表示为ϵsolvent，被设置为80。网格密度被设置为4.0网格/Å。然后，在接下来的量子力学计算中，每个被盖住的片段（CF）的表面电荷被添加为背景电荷。如上所述，蛋白质的部分电荷和筛选的表面电荷相互极化，直到达到收敛。只有当蛋白质的偶极子和表面电荷都收敛到一定的数值精度内时，这个循环才会停止。

计算公式见原文。

从图1的流程图中可以很容易理解生成PPCs的过程。对于本文报道的结果，单个蛋白质片段的电子密度的量子化学计算是在B3LYP/6-31G*水平进行的。

结果与讨论

蛋白质的溶化和极化

上一节所述的MFCC-PB计算协议需要进行广泛的数值测试以验证其功效。特别是，我们需要研究由PPC代表的蛋白质极化对蛋白质溶解的影响。我们进行了基准研究，考察了一些蛋白质系统的电子极化对溶解能的贡献。一些蛋白质系统的静电溶解能的计算贡献列于表1，为了比较，之前用MFCC-CPCM方法计算的结果也列于表中。与经典的PB或广义Born（GB）方法不同，溶质极化能在计算的静电溶存自由能中是相当重要的。我们的计算表明，如表1所示，这个项对总的静电溶解能的贡献是7∼19%。极化能的这一重要部分表明，在目前的力场中需要包括极化效应。

为了研究蛋白质极化对各种蛋白质特性的影响，我们还将目前PPC计算出的溶液中蛋白质的偶极矩与AMBER03力场中的偶极矩进行了比较。由于大多数蛋白质具有非零净电荷，而计算出的偶极矩取决于带电实体在坐标系中的位置，因此，蛋白质首先以其电荷中心为中心。电荷中心定义为rc=（∑ri∥qi∥）/∑∥qi∥，rc=（∑ri∥qi∥）/∑∥qi∥，其中ri是蛋白质中第i个原子的原始坐标，qi是其部分电荷。然后，整个蛋白质的偶极子由D=∑（ri-rc）qi.D=∑（ri-rc）qi.为了研究更详细的局部静电环境，我们计算了蛋白质单个残基的偶极矩。图2显示了一些蛋白质系统中由PPC计算的单个残基的偶极矩与AMBER03计算的偶极矩的对比图。比较结果表明，从PPCs计算出的大多数残基的偶极矩一般都比从AMBER电荷的大。这可能表明，在MFCC-PB电荷中，残基的极化程度更高。

Protein-specific的原子电荷

由于MFCC-PB电荷正确地描述了蛋白质的静电极化，所以在蛋白质不同位置的同一类型的氨基酸上partial charge通常是不同的。一个特定残基上的partial charge是由它的特定构象和由于蛋白质的其他残基而产生的化学环境决定的。为了说明这一特点，我们在表2中列出了溶菌酶的一些局部电荷。如同在溶菌酶中，我们发现同一类型的氨基酸上的电荷相互之间有很大的不同。例如，ASP中骨干原子N的原子电荷值从-0.083到-0.642，Ser中骨干原子O从-0.727到-0.556，Ile中原子CG1从0.025到0.304。这是PPC最重要的特点，它与标准的基于氨基酸的对一个给定的残基有相同的固定原子电荷的力场不同。

硫氧还蛋白中Asp26/Asp20的pKa移位计算（过）

尽管标准的基于氨基酸的部分电荷模型在MD模拟中对蛋白质的许多宏观特性建模时效果很好，但在模拟对局部静电环境更敏感的特性时预计会有困难。这是因为标准力场电荷是类似于均值场的，它们不包含蛋白质的极化和上述其他特定的蛋白质静电信息。例如，准确预测蛋白质中埋藏残基的pKa位移是目前MD模拟中的一个难题。准确计算pKa位移对于我们理解酶的酸碱催化机制至关重要。这种酶的活性需要催化残基以适当的质子化状态存在。pKa计算中常用的方法是基于求解PB方程，其中溶剂被视为连续介质。PB方法是一种平均场理论，对于一些简单的情况可以给出有用的见解。然而，PB方法在更复杂的情况下失败了，因为它不是一个微观模型，因此不能说明详细的分子过程。为了准确地预测pKa移动，人们需要正确地说明影响质子结合过程的分子因素。因此，需要有明确包括水分子的微观方法来从第一原理上正确预测pKa位移。然而，以前使用微观方法的一些尝试未能准确预测pKa位移。

Simonson等人最近进行了一项分子动力学自由能（MDFE）研究，计算硫氧还蛋白中Asp26的pKa位移，证明了目前标准AMBER和CHARMM力场的问题。用两种力场进行严格的MDFE模拟，不同的运行结果都高估了Asp26的pKa偏移量，误差为4-5千卡/摩尔。Simonson等人对pKa的MDFE计算可能是迄今为止最严格和最广泛的研究，其结果可以作为所研究系统的一个基准。Simonson等人的结论是，这个问题一定部分来自于所采用的力场的系统误差。在这篇文章中，遵循Simonson等人的严格的MDFE程序来计算埋藏的Asp26和表面暴露的Asp20的pKa偏移，但使用了本文之前描述的新的PPC。

我们的MDFE模拟的数值细节与Simonson等人的模拟密切相关。硫氧还蛋白的核磁共振结构（Protein Data Bank，即PDB，代码1XOA；并见图3）被作为初始结构。模拟是用AMBER程序进行的，使用了具有周期性边界条件的显式水模型（TIP3P）。首先，系统在NPT集合下平衡，在298K下运行500-ps，然后在NVE集合下继续模拟。SHAKE[44.](44)被用来约束与氢相连的共价键。长程静电相互作用是由粒子网格Ewald方法处理的。时间步长为一飞秒。我们使用的模型化合物是带有n-乙酰基和n-甲基酰胺阻断基的天冬氨酸。质子化和电离状态下的电荷分布见图4。

 表3给出了埋藏的Asp26以及表面暴露的Asp20和相应模型系统的计算结果。我们计算的Asp26去质子化的自由能转移是5.0千卡/摩尔，这非常接近实验值4.8千卡/摩尔。仿真一直进行到18 ns，为蛋白质的所有中间状态提供足够的采样。图5显示了电位相对于系统参数λ=0.5的能量导数与模拟时间的关系。这相当于Asp26的质子化和d质子化状态之间的一个虚构的中间状态。为了验证我们的计算，我们还计算了Asp20的pKa位移，一个表面暴露的残基。我们的模拟结果给出了1.0千卡/摩尔的自由能差，考虑到模拟的统计不确定性，这与实验值的零是相当合理的一致。我们的结果表明，使用正确描述蛋白质中带电残基周围的极化静电环境的PPC，可以从第一原理的严格的MDFE模拟中得到pKa移动的准确预测。相比之下，同样的数值模拟但使用AMBER和CHARMM电荷产生的pKa位移是Simonson等人中实验值的两倍。

为了估计统计误差，我们进行了块平均分析，其中每个持续时间为6ns的轨迹被分为四个块。导数的统计不确定性被估计为块平均的标准偏差的两倍。从表3可以看出，自由能的统计误差为1千卡/摩尔，计算的理论值在Asp26和Asp20的实验数据的统计误差范围内。由于我们的MDFE计算与Simonson等人[30.](30)的计算完全相同，只是我们用PPC取代了AMBER原子电荷，同时保持力场中的所有其他参数，所以不同的结果一定来自电荷的影响。换句话说，PPC中包含的蛋白质的极化效应是准确预测Asp26的pKa转变的原因。这一结果真正证明了电子极化在蛋白质微妙的静电相互作用中发挥的重要作用。

上述结果意味着电子极化的缺乏可能是早期对埋藏在蛋白质内部的可电离残基的pKa位移的一些MD自由能计算的主要罪魁祸首。当然，需要对不同的蛋白质系统进行更多这样的计算来证实这一结论。当然，还有其他重要的动力学效应被提出来解释pKa移动的MD模拟的失败。例如，Kato和Warshel提出，决定埋藏基团pKa值的一个主要因素是内部基团被电离时发生的结构变化。然而，由于有限的MD采样所带来的问题，大型结构重组的可能性并不容易在计算上进行研究。例如，到目前为止，葡萄球菌核酸酶中Glu66的pKa移动的计算是一个巨大的挑战。Kato和Warshel开发了一种特殊的过度充电方法，允许更多的结构重组，并获得了葡萄球菌核酸酶中Glu66的令人鼓舞的pKa位移。在现实中，我们预计电子极化和结构重组都应该在决定埋藏残基的pKa移位方面发挥一些作用，而且它们很可能相互交织在一起。我们相信，特定效应的实际贡献水平将取决于具体情况。显然，需要进一步的计算研究来澄清这种情况。

由于PPC代表了蛋白质在给定（原生）结构中的极化静电状态，我们预计在使用PPC的MD模拟中，与AMBER中的平均场状电荷相比，蛋白质的结构细节会有一些差异。我们使用PPC和AMBER电荷对硫氧还蛋白进行了进一步的MD计算，以比较结构细节。首先，我们检查了分别使用PPC和AMBER03进行的两组MD模拟的RMSD的特征。如图6所示，它们明显表现出不同的特征。由于两个计算中使用的力场中的所有其他参数都是相同的，除了电荷，图6的结果清楚地显示了蛋白质极化对详细的蛋白质动力学的重要影响，这在本质上是静电的。由于蛋白质内部氢键的形成使供体和受体极化，这种极化效应用PPCs比用AMBER电荷更好地表示。因此，在使用PPC的MD计算过程中，蛋白质中的氢键通常应该更加稳定，这可能解释了图6所示的硫氧还蛋白在MD模拟中的RMSD的不同特征

 总结

 我们提出了一种新的计算方案，通过结合蛋白质电子结构的线性缩放片段方案和溶剂的连续介质模型进行自洽处理，对溶液中的蛋白质进行量子力学计算。计算出的蛋白质原生结构中的电子密度被用来产生极化蛋白质力场（电荷）。派生的蛋白质特定的PPC正确地代表了蛋白质在其原始结构附近的极化静电状态，可以很容易地用于蛋白质运动的MD模拟。我们对硫氧还蛋白的MD模拟的计算研究揭示了PPC的以下特性。

1.
PPC是特定于蛋白质的，因此正确地描述了蛋白质在接近原始结构附近的极化静电状态。同一类型的残基上的电荷通常是不同的，这取决于它们特定的局部静电环境。这与基于氨基酸的电荷相反，后者是均值场的，无论在蛋白质中的具体位置如何，都保持不变。
2.
通过采用PPCs，使用热力学整合方法的MD自由能模拟准确地再现了硫氧还蛋白中埋藏的Asp26的实验性pKa偏移，而采用CHARMM和AMBER的标准平均场力场的相同计算则高估了pKa偏移，大约是两倍。

指出PPC的局限性也很重要。由于PPC是基于一个蛋白质的给定结构（大多数情况下是原生结构），用它来描述远离给定结构的结构可能不合适，也没有好处。幸运的是，大多数蛋白质的MD研究都是为了模拟接近原始结构的蛋白质运动，PPC在这种蛋白质的计算研究中应该是非常有吸引力的。