【2.5.3.2】amberTI结果分析
Small molecule binding to T4-lysozyme L99A
一、设置系统
绘制了两个配体，并用gaff atom types进行了参数化，并using antechamber on gaussian03 output files生成了相应的电荷。 （请参阅教程4，了解如何使用antechamber工具对配体进行参数化。苯和酚分子被保存在两个OFF库中（benz.lib和phen.lib）以备将来使用。截图显示选择了C6作为在酚中带有羟基的位置。
我们将使用来自pdb的T4-L99A的X射线结构（从中去除水分子和不需要的杂原子：pdb文件） 作为建立模拟文件的基础。我们将使用两次leap来生成四组参数并重新启动文件，其中包含处于蛋白质结合状态和溶剂化状态的两个配体。
- 第一个leap运行（输入文件 ) 将生成溶剂化和中和的苯配合物以及水中苯配体（complex.pdb和ligand.pdb ）的pdb文件。通过将BNZ分子重命名为PHN并删除H6 （t4_phn.pdb 和phn.pdb） ，可以从这两个附加的pdb文件中创建文件。
- 然后，将这四个pdb文件用于第二次leap运行 （输入文件）， 以生成* .prm和* .rst文件。这将产生4个parameter和4个rst文件
我们使用这个令人困惑的两步过程来构建参数文件是有原因的。在以后的TI运行中，将根据组合电势V（λ）耦合和传播不同配体的两个参数文件。这两个过程将具有一组独特的原子（苯中的H6和苯酚中的O1和H6）以及一组公共原子（水，芳环的其余部分和蛋白质）绑定的系统）。要求两个输入文件中的所有公共原子都具有相同的起始坐标并以相同的顺序出现，并且如果我们仅在两个不同的配体上使用了溶剂盒（即使它们仅相差几个原子），水分子的添加方式也将大不相同。
或者，一个人可以迅速建立一个系统，然后通过添加/删除原子和键将其修改为另一个系统。唯一重要的是，在最终输入文件中，两个过程（以及框的大小）所共有的原子的所有坐标都相同（或接近相同，由于舍入误差等原因，您可以避免出现较小的坐标偏差，将在sander启动期间通过警告消息进行更正）。请注意，设置不需要任何特殊类型的prmtop文件（如Amber7）或原子数与虚拟原子的平衡（如Amber9）。
二、输入文件准备
需要模拟两种不同的配体转化，即水中的苯转化为苯酚并与溶菌酶结合（无需模拟“反向”转化苯酚转化为苯，TI计算没有方向，而是每个λ点都代表苯的静态中间体）。出于模拟一致性的原因，将这两个转换进一步细分为三个子步骤：
1。 首先，去除氢原子上从苯中消失的原子部分电荷， 2. 其次进行配体vdW转换（从周围环境中消除或脱除苯氢，同时出现酚羟基） 3. 最后苯酚氢基团的原子部分电荷被打开。
这样做的原因是，当vdW与周围环境的相互作用变弱时，原子上具有非零电荷会导致模拟不稳定。同样出于模拟一致性的原因，在范德瓦尔斯变换中，我们将使用软核势，这是一种改良的LJ方程，可防止发生自由能发散的原点奇异类型（see Simonson T, Mol. Phys 1993: 80:441pp for more information）。
用于这些转换的mdin文件包含执行最小化/ MD模拟的常规Amber参数，除了一些用于执行TI计算的独特参数外，我将不作进一步讨论。 输入文件的非TI部分将执行以下步骤：
将2 fs的时间步与SHAKE算法，9 A和截止，各向同性压力换算和Langevin thermostat一起使用。 模板输入文件中有TI命令的占位符，其中包含以下关键字：
icfe = 1 告诉Sander这将是TI运行
clambda = XXX 设置为该运行中要使用的λ的实际值
ifsc = [0 | 1] 决定是否使用软核电势( softcore potentials)，指示输入文件中出现和消失的原子
scmask =='XXX' 将掩码字符串设置为此过程唯一的原子（即从V0消失或出现在V1中，另请参见下表）
crgmask ='XXX' 这是一个掩码字符串，指示在模拟开始之前某些原子应除去其原子部分电荷（这只是为了方便起见，可以构建单独的输入文件，其中电荷实际上设置为零，这将完成 与crgmask相同）
在此示例中，参数clambda设置为0.1、0.2，…，0.9。 根据转换的子步骤，ifsc，crgmask和scmask具有不同的值，下面对此进行了概述。 下表中的参数还将用于水中苯向苯酚的转化。
仅对于步骤2，激活软核电位（通过设置ifsc = 1）。 这使用了另一种形式的LJ方程，专门设计用于在出现或分散原子的情况下更好地收敛TI计算：
设置ifsc = 1 还会通知sander，每个prmtop文件（由scmask指定的原子）中的许多原子应被视为对此潜在函数唯一。在第1步和第3步的转换中，我们不需要使用软核电位，因为V0和V1之间的唯一区别是存在或不存在某些部分电荷，而起始态和终止态都具有相同数量的原子。
实际运行模拟将使用sander.MPI的multisander功能，例如：‘mpirun -np 2 sander.MPI -ng 2 -groupfile group’。由于输入文件众多，建议使用某种脚本自动生成它们。例如，像这样一个简单的shell脚本可以作为起点，该脚本可以设置针对配合物中苯的TI去除电荷的运行。该脚本会生成一个pbs输入文件，您可以根据自己的品牌的超级计算机排队系统对其进行调整。因此，最后，将执行总共两个转换，每个转换具有三个子步骤，每个转换具有九个λ窗口，每个具有三个 模拟，总共对sander模块进行162次调用（每个调用至少使用两个CPU）。真正的TI研究可能涉及更多的λ窗口和显着更长的模拟时间，这使其成为了计算量很大的技术。
三、运行和分析MD模拟
运行上一页概述的所有模拟之后，我们获得了大量的输出数据。 本节说明如何进行分析，以便最终获得我们两个配体的相对结合自由能的估计值。 第一步将是查看amber的mdout文件，以查看我们的物理参数是否有任何严重错误。 下表概述了六个不同模拟的输出文件：
Transformation: Remove charges from BNZ H6 - Medium: Water
Lambda DVDL rms ksteps TEMP rms PRESS rms Density rms Filename
0.100 -5.424 1.193 100 299.55 5.96 2.8 360.7 0.9842 0.0038 bnz_prod_v0_l1.out
0.200 -5.774 1.200 100 299.80 5.94 0.3 361.2 0.9833 0.0033 bnz_prod_v0_l2.out
0.300 -5.968 1.199 100 300.07 6.23 -0.7 358.2 0.9834 0.0042 bnz_prod_v0_l3.out
0.400 -6.205 1.213 100 299.31 5.99 1.4 356.7 0.9839 0.0038 bnz_prod_v0_l4.out
0.500 -6.488 1.250 100 299.50 5.95 1.9 347.2 0.9820 0.0035 bnz_prod_v0_l5.out
0.600 -6.681 1.261 100 299.46 5.90 -0.7 354.5 0.9841 0.0037 bnz_prod_v0_l6.out
0.700 -6.977 1.281 100 299.53 6.23 0.7 357.2 0.9844 0.0032 bnz_prod_v0_l7.out
0.800 -7.279 1.299 100 300.38 6.06 1.0 354.6 0.9818 0.0043 bnz_prod_v0_l8.out
0.900 -7.445 1.308 100 300.14 6.07 4.9 357.4 0.9835 0.0043 bnz_prod_v0_l9.out
Transformation: Transform BNZ into PHN - Medium: Water
Lambda DVDL rms ksteps TEMP rms PRESS rms Density rms Filename
0.100 5.188 3.076 100 300.01 5.94 0.6 354.1 0.9824 0.0040 bnz_prod_v0_l1.out
0.200 4.404 3.033 100 300.29 5.82 2.3 354.0 0.9832 0.0037 bnz_prod_v0_l2.out
0.300 3.660 3.010 100 299.67 5.97 -0.4 359.0 0.9835 0.0036 bnz_prod_v0_l3.out
0.400 2.523 3.112 100 299.80 6.06 -0.8 357.1 0.9817 0.0045 bnz_prod_v0_l4.out
0.500 1.367 3.188 100 300.24 6.09 -0.1 357.6 0.9830 0.0038 bnz_prod_v0_l5.out
0.600 0.602 3.468 100 299.52 6.04 -1.2 362.6 0.9820 0.0036 bnz_prod_v0_l6.out
0.700 -0.598 3.275 100 299.97 6.13 -0.3 359.5 0.9831 0.0037 bnz_prod_v0_l7.out
0.800 -1.757 3.102 100 300.21 6.15 4.1 368.2 0.9841 0.0041 bnz_prod_v0_l8.out
0.900 -2.695 2.989 100 299.84 6.01 2.7 357.4 0.9817 0.0037 bnz_prod_v0_l9.out
Transformation: Add charges to PHN O1,H6 - Medium: Water
Lambda DVDL rms ksteps TEMP rms PRESS rms Density rms Filename
0.100 -29.498 2.153 100 299.78 6.07 2.3 356.8 0.9850 0.0032 phn_prod_v0_l1.out
0.200 -30.213 2.233 100 299.92 6.05 0.0 361.6 0.9834 0.0038 phn_prod_v0_l2.out
0.300 -31.081 2.330 100 300.03 6.08 0.9 361.2 0.9835 0.0040 phn_prod_v0_l3.out
0.400 -32.142 2.486 100 300.34 5.95 1.3 350.8 0.9831 0.0033 phn_prod_v0_l4.out
0.500 -32.923 2.598 100 300.49 5.94 0.1 356.1 0.9824 0.0035 phn_prod_v0_l5.out
0.600 -34.052 2.774 100 300.45 6.02 0.8 352.7 0.9833 0.0041 phn_prod_v0_l6.out
0.700 -35.768 3.075 100 300.12 6.17 1.3 357.9 0.9826 0.0042 phn_prod_v0_l7.out
0.800 -37.675 3.415 100 300.21 5.93 -1.1 357.9 0.9826 0.0036 phn_prod_v0_l8.out
0.900 -39.873 3.655 100 300.49 6.06 3.0 364.9 0.9831 0.0041 phn_prod_v0_l9.out
Transformation: Remove charges from BNZ H6 - Medium: Protein
Lambda DVDL rms ksteps TEMP rms PRESS rms Density rms Filename
0.100 -5.380 0.697 100 299.99 1.67 0.2 135.8 1.0126 0.0011 t4_bnz_prod_v0_l1.out
0.200 -6.381 0.726 100 300.04 1.70 2.1 134.0 1.0125 0.0012 t4_bnz_prod_v0_l2.out
0.300 -5.967 0.823 100 299.87 1.70 0.9 136.0 1.0130 0.0012 t4_bnz_prod_v0_l3.out
0.400 -6.368 0.831 100 299.90 1.67 0.6 136.8 1.0127 0.0011 t4_bnz_prod_v0_l4.out
0.500 -6.502 0.733 100 300.02 1.68 1.6 133.1 1.0131 0.0010 t4_bnz_prod_v0_l5.out
0.600 -6.887 0.712 100 300.19 1.63 0.1 137.5 1.0127 0.0013 t4_bnz_prod_v0_l6.out
0.700 -6.747 0.764 100 300.09 1.71 1.1 135.7 1.0130 0.0011 t4_bnz_prod_v0_l7.out
0.800 -6.885 0.803 100 300.02 1.67 0.8 136.0 1.0128 0.0011 t4_bnz_prod_v0_l8.out
0.900 -7.104 0.766 100 300.15 1.68 1.5 137.6 1.0131 0.0013 t4_bnz_prod_v0_l9.out
Transformation: Transform BNZ into PHN - Medium: Protein
Lambda DVDL rms ksteps TEMP rms PRESS rms Density rms Filename
0.100 2.213 2.198 100 299.99 1.67 1.9 135.1 1.0129 0.0012 t4_bnz_prod_v0_l1.out
0.200 2.511 2.078 100 299.97 1.66 1.1 134.5 1.0128 0.0012 t4_bnz_prod_v0_l2.out
0.300 3.504 2.000 100 299.87 1.70 0.8 134.6 1.0130 0.0011 t4_bnz_prod_v0_l3.out
0.400 3.975 2.314 100 300.05 1.68 1.3 134.7 1.0130 0.0011 t4_bnz_prod_v0_l4.out
0.500 4.160 2.997 100 300.16 1.71 1.5 137.4 1.0126 0.0014 t4_bnz_prod_v0_l5.out
0.600 3.330 3.413 100 300.29 1.72 2.0 133.4 1.0126 0.0011 t4_bnz_prod_v0_l6.out
0.700 2.957 3.247 100 300.14 1.72 1.1 135.5 1.0128 0.0013 t4_bnz_prod_v0_l7.out
0.800 0.988 2.943 100 300.16 1.69 1.8 133.3 1.0128 0.0013 t4_bnz_prod_v0_l8.out
0.900 3.490 3.524 100 300.03 1.67 1.2 135.8 1.0124 0.0011 t4_bnz_prod_v0_l9.out
Transformation: Add charges to PHN O1,H6 - Medium: Protein
Lambda DVDL rms ksteps TEMP rms PRESS rms Density rms Filename
0.100 -32.831 1.170 100 300.00 1.69 1.9 132.4 1.0129 0.0011 t4_phn_prod_v0_l1.out
0.200 -33.164 1.169 100 300.09 1.66 0.4 134.2 1.0126 0.0013 t4_phn_prod_v0_l2.out
0.300 -33.385 1.322 100 300.00 1.72 0.7 133.9 1.0126 0.0012 t4_phn_prod_v0_l3.out
0.400 -33.828 1.271 100 299.98 1.74 1.1 135.0 1.0127 0.0012 t4_phn_prod_v0_l4.out
0.500 -33.840 1.247 100 300.06 1.69 0.6 136.6 1.0128 0.0012 t4_phn_prod_v0_l5.out
0.600 -34.164 1.189 100 300.00 1.68 1.3 134.9 1.0133 0.0011 t4_phn_prod_v0_l6.out
0.700 -34.454 1.293 100 299.98 1.69 1.9 137.1 1.0130 0.0013 t4_phn_prod_v0_l7.out
0.800 -34.659 1.341 100 300.01 1.70 1.1 137.7 1.0127 0.0015 t4_phn_prod_v0_l8.out
0.900 -34.388 1.318 100 299.93 1.67 1.1 133.9 1.0131 0.0011 t4_phn_prod_v0_l9.out
我们发现水中的配体温度波动在5K范围内，而配合物的温度波动在2K范围内，与人们预期的一样。 当使用ntp = 2时，压力会像往常一样剧烈波动，但密度保持相对恒定。 dV /dλrms值在几千卡/摩尔的范围内。 在实际研究中，对结果的快速浏览将需要用对所有轨迹的仔细分析来代替，以查看是否发生了有趣的构象变化或配体运动，但这对于本教程而言已足够。 因此，我们可以继续分析所得的dV /dλ与λ曲线。
下图显示了由9个λ点组成的dV /dλ曲线。 我们可以看到，水中的转化曲线比蛋白质中的配体变化曲线更平滑，这是可以预期的，这是由于与各向同性溶液相比，蛋白质-配体复合物的构象更为复杂。 复合物中的vdW转换显示 dV /dλ 的波动最大，而复合物中步骤2的自由能曲线是最“颠簸”的曲线。 由于TI的自由能曲线几乎总是平滑的曲线，这很可能是由于采样不足以及对复杂的转换进行更长的模拟而导致的自由能曲线很可能会被平滑化。
当以络合物或溶液形式进行时，除去氢原子上的电荷的热力学作用很小，并且非常相似。 将电荷添加到苯基上（步骤3）会比其他两个步骤产生更大的能量变化，但请记住，单个转换的结果并不是物理上有意义的数字，只有两个转换结果之间的差异。 与每个步骤相关的自由能变化可以通过对 λ= 0-1 右边的ΔG（复合水）曲线进行数值积分来计算。 集成可以像这样，可以在具有λ，dV /dλ和rms（dV /dλ）值的列表上使用该脚本完成此操作）。 梯形法则（从最接近的两个值线性推断出λ= 0和1的值）得出：
有趣的是，对结合自由能差的主要贡献来自vdW相互作用，而不是静电。可以预料到，极性羟基上的溶剂化作用的丧失可能会强烈影响相对的结合自由能，但是根据这一计算，主要区别似乎是在静电作用下，使结合位点中较大的羟基适应了空间效应。扮演次要角色。但是，如果进一步完善计算，则此结果可能会更改。
我们的预测是，苯酚对T4-溶菌酶的结合能力比苯弱1.8 kcal / mol。尽管这与实验数据并不完全吻合，但它抓住了苯成为更强粘合剂的趋势。为了改善结果，可以在不同的λ点处添加更多的模拟（这确实是TI的强项之一，您可以根据需要添加尽可能多的其他数据点以精炼结果，而不必重做初始计算）或可以使用更复杂的数值积分方案。最终结果的均方根值不应被误认为是结果的统计不确定性。为了获得良好的误差估计，需要一些附加的分析（参见例如Steinbrecher T.等人，J Chem.Phys.127，214108，2008）。
参考资料
