【2.5.2】amber使用热力学积分进行pKa计算
本教程将重现蛋白质硫氧还(thioredoxin)蛋白中ASP残基的pKa值的计算，如下文所述：
Simonson, T., Carlsson, J., Case, D.A., “Proton Binding to Proteins: pKa Calculations with Explicit and Implicit Solvent Models”, JACS 2004, 126, pp4167-4180.
在尝试本教程之前，您应确保已阅读该文献。如果您在2006年Amber高级研修班中，则这是必备阅读材料的一部分。
本文的作者使用隐式和显式溶剂以及AMBER和CHARMM力场进行了许多模拟。就本教程而言，我们将仅针对蛋白质硫氧还蛋白中的单个天冬氨酸残基（ASP 26）设置隐式溶剂模拟。欢迎您自行尝试和/或其他天冬氨酸残基或其他蛋白质的溶剂实例。
pKa值的计算基于以下公式：
pKa,model的值是模型系统的实验pKa。 在这种情况下，ACE/NME阻止了天冬氨酸。 我们实际上对计算一部分蛋白质基质中的天冬氨酸残基的pKa位移感兴趣。 在这种情况下，我们不需要知道模型系统的实验pKa。 我们感兴趣的是获得上述方程中涉及delta-delta G的后半部分的值。这表示pKa漂移，这就是我们要计算的结果。 我们将使用本页顶部图像中显示的热力学循环来执行此操作。
为此，我们需要使用热力学积分来计算将模型AspH转换为Asp-的deltaG和将硫氧还蛋白-AspH，转换为硫氧还蛋白-Asp-的delG。 我们将从对模型化合物的deltaG计算开始。
本教程共包括4个步骤，如下所示：
- 为AspH和Asp-模型设置prmtop和inpcrd文件。
- 对模型化合物运行热力学积分，以找到delta-G（model）。
- 设置AspH和Asp-硫氧还蛋白蛋白的prmtop和inpcrd文件。
- 对硫氧还蛋白进行热力学整合。 计算deltadeltaG和pKa偏移。
一、设置AspH和Asp-模型的prmtop和inpcrd文件。
我们将使用Amber的IGB = 5模型和FF94力场进行计算。较新的力场可用，也许更合适，但出于本教程的目的，我们将使用与本文相同的方法。 Amber 9中进行热力学积分的方法已从早期版本更改。您不再需要创建包含扰动数据的prmtop文件，而是创建代表热力学循环两端的两个单独的prmtops。要求是这两个prmtop文件具有完全相同的原子数和原子顺序。这意味着我们不能简单地创建ACE-ASH-NME和ACE-ASP-NME，而是需要遵循它们在本文中使用的方法来创建两个ACE-ASH-NME系统，但是在第二种情况下，我们将零使质子上的电荷带电，并使其他原子上的电荷重新分配，使其与ASH残基的定义一致。我们用这种方法会遇到的一个问题是如何使涉及质子的VDW术语消失。幸运的是，在FF94力场中，所有与氧原子相连的质子的范德华半径为零。这意味着与羟基质子不存在VDW相互作用，因此将其电荷清零将消除与该质子的所有非键相互作用。我们不必担心删除键，角度和二面角项，因为它们都是残基的内部，因此这些项的平均属性在模型和蛋白质系统中应该相同。这些术语对lambda势的导数也没有贡献。
因此，我们将建立两个ACE-ASH-NME系统，并收取以下电荷，如图2所示。注意到本文的工作时，他们假设仅质子化时羟基质子上的局部电荷才会发生变化，因此对于ASP残基，他们将两个氧原子的电荷设置为与FF94力场值相匹配，但随后将所有其他原子变成C（γ）原子。因此，这两个系统要使用的电费为：
最高电荷与ASH的力场中定义的电荷相同。 您可以根据需要在Xleap中自行检查。
我创建了一个简单的跳转脚本来为我们进行设置。 这里是：
create_ash.leaprc
source leaprc.ff94
set default PBradii mbondi2
mol=sequence{ACE ASH NME}
saveamberparm mol asph_model.prmtop asph_model.inpcrd
desc mol.2
set mol.2.OD2 charge -0.8014
set mol.2.OD1 charge -0.8014
set mol.2.CG charge 0.5307
set mol.2.HD2 charge 0.0000
saveamberparm mol asp_model.prmtop asp_model.inpcrd
exit
终端运行
$AMBERHOME/exe/tleap -s -f create_ash.leaprc
您应该检查输出以确保其正常工作。 即 ASPH模型电荷为-1.0吗？ 等等
这将创建：asph_model.prmtop，asph_model.incprd，asp_model.prmtop和asp_model.inpcrd
下一步是对该模型化合物进行热力学积分计算，以找到ΔG（模型）。
二、在模型化合物上进行热力学积分，以找到delta-G（模型）
现在，我们将对模型化合物进行热力学积分计算。在Amber 9中，热力学整合是使用multi-sander以并行方式完成的。这是对并行版本的sander的包装，该版本在单个并行运行中运行多个sander副本。由于这个原因，您必须在系统上安装并测试过并行版本的sander（sander.MPI）。您还必须知道如何在系统上启动任何必需的MPI demons，并成功执行在两个cpus上运行的并行sander。如果您在工作室中进行本教程，则应该已经在路径中安装了必需的软件，并且正在运行lamd mpi demon。 （使用’ps -aux’进行检查-如果看不到lamd运行，请使用lamboot启动它）。
我们将在模型系统上总共运行5个lambda值，每个3ns。这与他们在论文中所做的相同。我们将使用igb = 5，摇动所有氢（ntc = 2，ntf = 2），使用1fs的时间步长（dt = 0.001）。在本文中，他们使用Berendsen thermostat进行温度控制。但是，不再建议这样做。现在，温度控制的首选方法是使用Langevin动力学（ntt = 3）。由于这只是一个小型模型系统，因此在执行MD之前不必最小化我们的系统就可以摆脱困境。
我们将在此处运行5个模拟。第一个不会重新启动（irest = 0，ntx = 1）。其他4个将重新启动。每种情况下的差异将是clambda的值。它在amber手册中对此进行了更详细的说明。第一个模拟（irest = 0）将使用clambda = 0.0。其他4个将分别使用clambda = 0.11270、0.5000、0.88729和1.0000。以下是第二步的示例mdin文件。您应该设置5个mdin文件（记住step1的ntx = 1和irest = 0）。 ntave = 1000使sander每1000步写一次运行平均值。这使得更容易检查平衡。
model_step2.mdin
model ASH to ASP for pKa calculation
&cntrl
nstlim =3000000, nscm=2000,
ntx=5, irest=1, ntpr=1000,
tempi=300.0, temp0=300.0, ntt=3, gamma_ln=5.0,
ntb=0, igb=5, cut=999.0,
dt=0.001,
ntc=2, ntf=2,
ntwr = 10000, ntwx=1000, ntave=1000
icfe=1, clambda=0.11270,
/
对于这5次运行，我们还将需要一个组文件。 该文件为将要运行的每个sander线程指定命令行选项。 这是步骤2的示例：
model_step2.group
-O -i model_step2.mdin -p asph_model.prmtop -c asph_model_step1.rst -o asph_model_step2.out -inf asph.mdinfo -x asph_model_step2.mdcrd -r asph_model_step2.rst
-O -i model_step2.mdin -p asp_model.prmtop -c asp_model_step1.rst -o asp_model_step2.out -inf asp.mdinfo -x asp_model_step2.mdcrd -r asp_model_step2.rst
创建所有5个文件后，您可以在一个运行脚本中运行所有5个作业。 例如。：
run_model.x
mpirun -np 2 $AMBERHOME/exe/sander.MPI -ng 2 -groupfile model_step1.group
mpirun -np 2 $AMBERHOME/exe/sander.MPI -ng 2 -groupfile model_step2.group
mpirun -np 2 $AMBERHOME/exe/sander.MPI -ng 2 -groupfile model_step3.group
mpirun -np 2 $AMBERHOME/exe/sander.MPI -ng 2 -groupfile model_step4.group
mpirun -np 2 $AMBERHOME/exe/sander.MPI -ng 2 -groupfile model_step5.group
-np指定2个cpus。 原则上，这可以是任何偶数，并且会将处理器分成两组。 例如。 如果您拥有4个处理器的计算机或群集，则可以执行-np 4以使用所有4个处理器。
此计算大约需要4400 s（1.22小时）才能在奔腾D 3.2GHz机器上运行所有5个lambda。
这是一个包含输出文件的tar文件：model.tar.gz（12.25 MB）
2.1 检查收敛 Check for convergence
接下来，我们需要检查我们最初的3ns运行是否足以平衡我们的系统。 实际上，我们应该对本教程中的内容进行更深入的分析，但是在此所做的操作对于我们的目的应该足够了。 我们应该做的第一件事是绘制5个lambda值的1000步瞬时平均值和累积平均值。 我们可以编写一个简单的awk脚本来从输出文件中提取相关数据：
dvdl.awk
BEGIN{ isd=0; epsum=0; npts=0 }
/Molecular dynamics:/{
while( $1 != "nstlim") getline;
split($3,a,",");
nstlim = a[1];
}
/DV\/DL,/{ isd=1}
/NSTEP/{ if(isd == 1 ){nstep = $3 ; tim = $6; }}
/EPtot/ {
if(isd == 1 ) {
if( nstep >= nstart ) {
npts += 1;
eptot = $9;
epsum = epsum + $9;
print tim,eptot,epsum/npts ;
}
isd = 0;
}
}
您应该在5个lambda值输出文件中的每一个上运行此awk脚本（选择ASP或ASPH都没关系，因为它们应该相同）以获得1000步瞬时 DV/DL值和累积平均值的列表 DV/DL值。 例如。：
awk -f dvdl.awk asp_model_step1.out > step1.dat
首先，您应该在每个stepx.dat文件的末尾尾部检查最后一行上的两个DV / DL值是否匹配。 这是由sander计算的总体平均值和由awk脚本计算的全部累计平均值。 如果值不匹配，则说明发生了严重错误，您应该仔细检查输出文件中是否有错误迹象。 例如。
tail -n 5 step1.dat
2997.000 -58.0952 -57.8712
2998.000 -58.0155 -57.8712
2999.000 -58.2206 -57.8713
3000.000 -58.0181 -57.8714
3000.000 -57.8714 -57.8714
注意，由于不同体系结构之间的舍入差异，您得出的结果可能略有不同。 但是，如果我们已经足够长时间地进行了平衡，那么我们获得的最终结果应该与采样误差内的结果相同。
现在，您可以打开xmgrace（或您自己喜欢的制图程序），并检查5个stepx.dat文件中的每个文件是否收敛。 从命令行运行xmgrace，然后选择 Data->import->ascii ，选择文件名， for the Load as box select NXY。 （单击图像可查看大图）
从这些图像可以看出，所有5次运行的累计平均值（以红色显示）已经收敛到可以接受的公差范围。 注意瞬时平均值中的大量跳跃。 这就是为什么必须进行长时间的大规模模拟并收集足够的数据的原因。 由于运行似乎已经收敛，因此我们现在可以使用它们来提取计算delta-G（model）所需的数据。
2.2 计算delta-G（模型）
现在我们有了输出文件，可以通过使用论文中的方程5来计算delta-G（model）以找到 dG/dL：
然后，我们将其在从零到一的lambda范围内进行数值积分，以得到delta-G（model）。
我们使用 dU/dL（DV/DL）的平均值的事实就是为什么我们需要确认数据的收敛性。 现在，我们需要选择要提取的数据部分，以供上式使用。 从上面的图表中，我们将舍弃每3 ns运行的第一个纳秒作为平衡阶段。 因此，我们需要再次处理输出文件，并且仅从每个文件中提取最后2,000,000步（1 fs时间步）的数据。 上面的awk脚本可以为我们做到这一点。 我们只需要告诉它开始的步骤：
awk -f dvdl.awk nstart=1000000 asp_model_step1.out > prod1.dat
文件:
您应该再次在xmgrace中可视化这些文件，以验证累积平均值看起来不错。
用高斯求积法计算deltaG（model）
现在，我们将根据上述计算中的数据估算以下积分，从而计算deltaG（model）的值：
我们使用下面显示的公式来计算DV / DL下的面积。 我们使用高斯正交进行此操作，这比使用梯形方法要精确得多：
下表总结了我们需要的数据，这是prod2-4.dat的最后一行，它是每2,000个步骤块的累积平均值。 有关Lambda和Wi值的来源的更多信息，请参考Amber 9手册的第155页：
由于对于0.0和1.0的lambda值，Wi的值为零，因此这两个步骤对积分的值没有贡献。 这样，我们实际上可以通过不计算终点来节省自己的模拟时间。 但是，这是一种很好的做法，如果我们在模拟方面遇到问题，它将使调试工作变得更加容易。 因此，上述积分的值为：
deltaA = (-58.3322 * 0.27777) + (-59.6629 * 0.44444) + (-61.7819 * 0.27777) = -59.8807 kcal / Mol
在这种情况下，deltaA = deltaG，因此我们的deltaG（model）值为-59.8807 kcal / Mol。 不幸的是，由于这不是在恒定温度和压力下的deltaG值，因此我们无法使用它直接计算模型的pKa。 但是，我们可以将其与蛋白质模拟结合使用以获得deltadeltaG，从中我们可以计算出蛋白质基质中残基的delta pKa（pKa位移）。
现在，我们准备转到第3部分，在这里我们将为AspH和Asp-thioredoxin蛋白构建输入文件。
三、设置AspH和Asp-thioredoxin蛋白的prmtop和inpcrd文件。
现在，我们已经计算出模型化合物的delta-G值，我们必须对硫氧还蛋白进行重复此过程。 这是蛋白质的pdb文件。
该蛋白质的原始晶体结构（2TRX）包含1条以上的链，并且某些残基还具有许多其他构象。 在上面给出的pdb文件中，我已经提取了链A，并且还为所有重复的残基选择了构象A。 但是，我们仍然需要做一些初步的准备。
首先，我们需要检查任何二硫键。 我们可以在VMD中直观地执行此操作。 打开pdb文件，并显示球和棒中的所有半胱氨酸残基：
如您所见，此蛋白质中只有两个半胱氨酸残基，并且它们之间显然存在二硫键，因此我们需要编辑pdb文件并将它们都更改为CYX。 创建prmtop和inpcrd文件时，我们还必须在Leap中手动添加此绑定。
接下来，我们需要指定任何组氨酸残基的质子化状态。 正如您在下面看到的那样，该蛋白质中只有一个组氨酸残基，它暴露于溶剂中，距离我们将计算的pKa p26残基很远。
因此，该组氨酸的质子化状态不太可能对我们从计算中获得的结果产生很大影响。 出于本教程的目的，我们将遵循它们在本教程所基于的论文中指出的内容：
“The unique histidine (His6) was doubly protonated. This residue is solvent-exposed, 10 � away from Asp26, so that its protonation state is not expected to affect Asp26 strongly.”
因此，通过将其名称从HIS更改为HIP，我们将使HIS6倍增质子化。
这是修改后的pdb文件：thioredoxin_chaina_modified.pdb
接下来，我们需要为Asp26的两个质子化状态创建prmtop和inpcrd文件。 我们将使用IGB = 5进行此计算，因此必须记住将默认的PBradii设置为mbondi2集。 这是用于创建prmtop的leap文件：
thio_leap.in
source leaprc.ff94
set default PBradii mbondi2
thio = loadpdb thioredoxin_chaina_modified.pdb
#Add the disulphide bond
bond thio.32.SG thio.35.SG
saveamberparm thio thio_ash.prmtop thio_ash.inpcrd
#Adjust the charges to match the Asp- charges from the paper
set thio.26.OD2 charge -0.8014
set thio.26.OD1 charge -0.8014
set thio.26.CG charge 0.5307
set thio.26.HD2 charge 0.0000
saveamberparm thio thio_asp.prmtop thio_asp.inpcrd
quit
终端运行
$AMBERHOME/exe/tleap -s -f thio_leap.in
您应该检查输出是否有任何错误。 尤其要检查以确保带-1电子的ASP的电荷比ASH的“更多”。
以下是文件：
现在，我们准备进行下一步，在此步骤中将结构最小化和平衡，然后进行热力学积分。
四、对硫氧还蛋白进行热力学整合。
在我们真正开始进行热力学积分运行之前，我们应该迅速最小化我们的系统，以消除Leap氢化中产生的任何不良接触。但是要注意的一件事是，因为我们正在进行热力学积分，所以这两个状态的端点的坐标文件必须相同。因此，目前我们有一个thio_ash.inpcrd和thio_asp.incrd，它们的名称不同。现在，如果我们分别将这两者最小化，我们将得到两个不同的结构，这些结构将无法用于进行热力学积分。但是，由于我们的目的只是在启动MD之前清理结构，因此我们可以只运行一种状态，然后将结果结构用于ASP和ASH起始结构。在本文中，他们没有提及所使用的nonbond cut off的值，因此我们必须假设他们没有使用 cut off。通常，在GB模拟中，如果可能的话，希望避免使用截止值。如果出于计算费用的原因必须使用截止值，则理想情况下应为18埃或更大。我们还将半径计算的截止值设置为也大于系统。
thio_min.in
minimise thioredoxin
&cntrl
imin=1, maxcyc=500, ncyc=250,
igb=5, ntb=0, ntf=2, ntc=2,
cut=999.0, rgbmax=999.0,
ntpr=50
/
终端运行
$AMBERHOME/exe/sander -O -i thio_ash_min.in -o thio_ash_min.out -p thio_ash.prmtop -c thio_ash.inpcrd -r thio_ash_min.rst
输出文件：
下一步是运行热力学积分。 我们将以与模型系统相同的方式执行此操作。 在本文中，它们运行了一个加热阶段和一个平衡阶段，并对晶体结构进行了谐波(harmonic)约束，并在平衡阶段逐渐降低了晶体结构。 总体而言，它们运行了大约6 ns的模拟。 在本教程中，我们无力承担此费用，因此出于我们的目的，我们将以300K热启动，并运行5个lambda值模拟，每个模拟总共1 ns。 然后，我们将查看每个模拟，并可能舍弃前200ps左右的平衡。 在真实的计算中，人们希望更加仔细地平衡，并且还要进行更长的生产。
thio_step2.in
thioredoxin ASH to ASP for pKa calculation
&cntrl
nstlim = 1000000, nscm=2000,
ntx=5, irest=1, ntpr=1000,
tempi=300.0, temp0=300.0, ntt=3, gamma_ln=5.0,
ntb=0, igb=5, cut=999.0,
dt=0.001,
ntc=2, ntf=2,
ntwr = 10000, ntwx=1000, ntave=1000
icfe=1, clambda=0.11270,
/
thio_groupfile.2
-O -i thio_step2.in -p thio_ash.prmtop -c thio_step1_ash.rst \
-o thio_step2_ash.out -inf ash.mdinfo -x thio_step2_ash.mdcrd
-r thio_step2_ash.rst
-O -i thio_step2.in -p thio_asp.prmtop -c thio_step1_asp.rst \
-o thio_step2_asp.out -inf asp.mdinfo -x thio_step2_asp.mdcrd \
-r thio_step2_asp.rst
您应该对所有5个lambda值运行以上计算。 您可以编写脚本。 例如。：
mpirun -v -machinefile $PBS_NODEFILE -np 128 /usr/local/apps/amber9/exe/sander.MPI \
-ng 2 -groupfile thio_groupfile.1
mpirun -v -machinefile $PBS_NODEFILE -np 128 /usr/local/apps/amber9/exe/sander.MPI \
-ng 2 -groupfile thio_groupfile.2
mpirun -v -machinefile $PBS_NODEFILE -np 128 /usr/local/apps/amber9/exe/sander.MPI
-ng 2 -groupfile thio_groupfile.3
mpirun -v -machinefile $PBS_NODEFILE -np 128 /usr/local/apps/amber9/exe/sander.MPI
-ng 2 -groupfile thio_groupfile.4
mpirun -v -machinefile $PBS_NODEFILE -np 128 /usr/local/apps/amber9/exe/sander.MPI
-ng 2 -groupfile thio_groupfile.5
请注意，与Myrinet连接的1.7GHz Itanium 2群集的128 cpus大约需要25个小时。 现在也许您知道为什么我们只运行1 ns。 您可以在此处下载输出的tar文件：thio_calc_results.tar.gz。
我们再次需要检查 DV/DL累积值的收敛性。 我们会将每次运行的前200ps数据视为平衡。 与之前一样，我们实际上只需要点2,3和4，您可以使用前面的awk脚本来获取它们：
从上图可以看出，特别是对于λ= 0.5和0.872的情况，我们尚未达到累积DV / DL平均的收敛。 实际上，如本文所述，存在多个状态是显而易见的。 因此，为了获得准确的结果，我们将需要进行更多的平衡。 但是，出于本教程的目的，我们将继续进行操作。
用高斯求积法计算deltaG（model）
对于我们的模型系统，我们现在将通过估算以下内容来计算deltaG（model）的值：
根据我们从上面得出的结果。 那就是我们对以下项进行高斯求积：
下表总结了我们需要的数据，该数据是 2-4_800_dvdl.dat 的最后一行，这是每800个步骤块的累计平均值。 有关Lambda和Wi值的来源的更多信息，请参考Amber 9手册的第155页：
因此，上述积分的值为：
deltaA = (-49.8394 * 0.27777) + (-53.6001 * 0.44444) + (-62.1712 * 0.27777) = -55.0886 kcal / Mol
在这种情况下，deltaA = deltaG，因此我们对deltaG（蛋白质）的值为-55.0886 kcal / Mol。
deltadeltaG和pKa偏移的计算
现在我们有了deltaG（模型）和deltaG（蛋白质），我们可以通过参考本教程开始时显示的热力学循环来计算deltadeltaG：
从而：
deltadeltaG = -55.0886--59.8807 = 4.7921 kcal / mol
如果将其与文献进行比较，您会发现他们在隐式溶剂中计算出的ASP26的deltadeltaG为7.1。 这显然比我们的大，这可能是因为我们没有足够长的采样时间。 但是，有趣的是，实验值是4.8，因此我们很满意。 尽管我认为这更多是靠运气而不是判断。
从这个deltadeltaG中，我们可以使用以下方程式计算pKa位移（显然只是忽略了pKa（model）项）：
= (1/2.303kT)*4.7921
= 1/(2.303*0.001987*300)*4.7921
pKa shift (ASP26) = 3.49
因此，我们计算出由硫氧还蛋白蛋白基质引起的ASP26的pKa位移为3.49。 这在实验中很明显。 虽然我们在这里可能只是幸运的。
参考资料
