仓库源文，站点原文

layout: post title: AMBER高级教程A12：使用Amber11处理聚脯氨酸 categories:

摘要

本文介绍了几种利用sander程序(Amber 11中的MPI程序包)评估一条短聚脯氨酸链构象是否达到平衡的方法. 特别地, 本文还阐述了所谓的"拉伸分子动力学"和著名的"副本交换"[1]协议. 这篇文章的意图是公开技术细节而不是描述基础知识.

引言

当准备这个教程的时候, 我面临着一个不可避免的选择: 使用一个在1个小时之内就很容易运行的简单而又相对乏味的例子, 还是考虑一个更大的, 需要更长的时间运行的, 且稍微有趣一点的例子. 这次后者更受我的青睐, 因此这篇文章描述了在Amber 11中模拟Ace-(Pro)5-NHe肽的两种方法. 重复这篇教程需要的所有文件和脚本以及这篇文档都可以在以下地址下载: polyproline-tutorial-files.tar.bz2(提示: 文件包大小大约为60M)

初步准备

我们使用teLeap程序获得prmtop和inpcrd文件, 输入文件如下

<div class="highlight"><pre style="line-height:125%"><span style="color:#A2F">source</span> leaprc.ff99SB <span style="color:#A2F">m</span> = sequence {ACE PRO PRO PRO PRO PRO NHE} <span style="color:#A2F">set</span> default PBradii mbondi2 <span style="color:#A2F">saveamberparm</span> m prmtop inpcrd <span style="color:#A2F">savepdb</span> m inpcrd.pdb <span style="color:#A2F">quit</span> </pre></div>

图1 生成Ace-(Pro)5-NHe肽的prmtop和inpcrd文件所需的LEaP程序输入文件

我们准备使用文献[2, 3, 4]中的GB/SA溶剂模型(igb=2, gbsa=1), 因此第五行令teLeap程序使用合适的半径. 其余内容应该是不言自明的.

用常规分子动力学方法观察室温下的遍历性是有指导意义的(比如在300K条件下). 鉴于此, 我们可以使用mdin文件如下所示.

<div class="highlight"><pre style="line-height:125%"><span style="color:#A2F">****</span> this line is required **** <span style="color:#A2F">&cntrl</span> <span style="color:#A2F">ntwx</span> = 0, ntpr = 5000, ntwr = 50000, <span style="color:#A2F">ntt</span> = 3, temp0 = 300.0, gamma_ln = 1.0, <span style="color:#A2F">igb</span> = 2, gbsa = 1, dielc = 1.0, cut = 18.0, <span style="color:#A2F">ntb</span> = 0, ntc = 2, ntf = 2, tol = 0.000001, <span style="color:#A2F">nstlim</span> = 500000, dt = 0.002, ntp = 0, ibelly = 0, <span style="color:#A2F">ntr</span> = 0, imin = 0, irest = 0, ntx = 1, ig = 27606 <span style="color:#A2F">/</span> <span style="color:#A2F">ncsu_abmd</span> <span style="color:#A2F">mode</span> = ANALYSIS <span style="color:#A2F">monitor_file</span> = 'monitor.txt' <span style="color:#A2F">monitor_freq</span> = 500 ! 1 ps <span style="color:#A2F">variable</span> <span style="color:#A2F">type</span> = COS_OF_DIHEDRAL <span style="color:#A2F">i</span> = ( 2, 5, 7, 17, ! :1@CH3 == :1@C == :2@N == :2@CA <span style="color:#A2F">17,</span> 19, 21, 31, ! :2@CA == :2@C == :3@N == :3@CA <span style="color:#A2F">31,</span> 33, 35, 45, ! :3@CA == :3@C == :4@N == :4@CA <span style="color:#A2F">45,</span> 47, 49, 59, ! :4@CA == :4@C == :5@N == :5@CA <span style="color:#A2F">59,</span> 61, 63, 73) ! :5@CA == :5@C == :6@N == :6@CA <span style="color:#A2F">end</span> variable <span style="color:#A2F">end</span> ncsu_abmd </pre></div>

图2 室温下做常规分子动力学的sander程序输入文件

除了常规的&cntrl参数外, 还包括ncsu_abmd部分(定义了集合变量的类型COS_OF_DIHEDRAL, 17-25行), 命令sander程序每500步 (1 ps) 保存一次变量的值到monitor.txt文件中 (12,14和15行). 集合变量由连续连接脯氨酸残基的肽键的二面角的余弦值之和来确定(见图2中20-24行的注释). 对于顺式构象中的所有链, 变量值可以设成5, 而对于反式构象中, 值应该等于-5. 一旦模拟完成, 检查monitor.txt文件会发现在1 ns内没有反式/顺式构象转变发生. 当然也可以尝试长一点的时间和更高的温度, 这些都是有意义的. 总结而言, 在此情况下, 常规分子动力学并不是具有遍历性的: 轨迹没有达到所有的具有玻尔兹曼概率的平衡构象状态; 它卡在初始构象附近, 因为在分开的势能面极小值之间存在着高高的自由能势垒. 尽管有这些势垒存在, 但这个教程的目的是演示两个可能的方法去学习研究平衡态系综. 在下一个环节, 我们将考虑所谓的拉伸分子动力学, 它可以用来计算不同构象之间的自由能差. 然后, 我们将讨论一个更苛刻(但更强大)的"副本-交换"技术.

拉伸分子动力学

所谓的拉伸分子动力学的意思是, 在一个或一组集合变量上增加一个非稳态的谐振约束. 这种方法在文献中被广泛讨论, 因此我们跳过所有细节, 将读者直接引向文献[5, 6, 7]和这三篇文献中引用的文献. 为想偷懒的人简单总结一下: 这种方法可以将系统在不同态之间"拉来拉去", 衡量非平衡功的执行, 然后推导出相应的平衡自由能差.

Sander中的几个子程序可以实现"拉伸"功能, 本教程只描述一个子程序(ncsu_smd), 对这个子程序或者其它相关子程序有疑问的读者可以参看用户手册.

我们对COS_OF_DIHEDRAL变量进行拉伸, 改变它的值在-5到5之间(也就是说, 让分子在全反式构象和顺式构象之间变换). 为了获得更加准确的Δf值, 我们进行了100次向前拉伸(从-5到+5)以及100次向后拉伸(从+5到-5). 每个用于ncsu_smd操作的初始坐标都在相应的状态下达到平衡. (全反式的向前拉, 全顺式的向后拉). Mdin文件的相关部分如图3所示, 内容也应该是不言自明的. 每10ps会程序会将拉伸的过程写进work.txt文件. 一旦模拟完成, 总功将会在这个文件的末尾显示. 所有两百次拉伸的文件可以在polyproline-tutorial-files/3.smd文件夹下找到, 我们鼓励读者检查或者重复运行它们. 文件夹下还有一个perl脚本可以用来计算自由能差.

<table class="highlighttable"><th colspan="2" style="text-align:left">bash</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%"><span></span>**** this line is required **** &cntrl ... / ncsu_smd <span style="color: #B8860B">output_file</span> <span style="color: #666666">=</span> <span style="color: #BB4444">'</span>output/00<span style="color: #BB4444">'</span> <span style="color: #B8860B">output_freq</span> <span style="color: #666666">=</span> <span style="color: #666666">5000</span> ! <span style="color: #666666">10</span> ps variable <span style="color: #AA22FF">type</span> <span style="color: #666666">=</span> COS_OF_DIHEDRAL <span style="color: #B8860B">i</span> <span style="color: #666666">=</span> <span style="color: #666666">(</span> 2, 5, 7, 17, ! :1@CH3 <span style="color: #666666">==</span> :1@C <span style="color: #666666">==</span> :2@N <span style="color: #666666">==</span> :2@CA 17, 19, 21, 31, ! :2@CA <span style="color: #666666">==</span> :2@C <span style="color: #666666">==</span> :3@N <span style="color: #666666">==</span> :3@CA 31, 33, 35, 45, ! :3@CA <span style="color: #666666">==</span> :3@C <span style="color: #666666">==</span> :4@N <span style="color: #666666">==</span> :4@CA 45, 47, 49, 59, ! :4@CA <span style="color: #666666">==</span> :4@C <span style="color: #666666">==</span> :5@N <span style="color: #666666">==</span> :5@CA 59, 61, 63, 73<span style="color: #666666">)</span> ! :5@CA <span style="color: #666666">==</span> :5@C <span style="color: #666666">==</span> :6@N <span style="color: #666666">==</span> :6@CA <span style="color: #B8860B">path</span> <span style="color: #666666">=</span> <span style="color: #666666">(</span>-5.0, 5.0<span style="color: #666666">)</span> <span style="color: #B8860B">harm</span> <span style="color: #666666">=</span> <span style="color: #666666">(</span>100.0<span style="color: #666666">)</span> end variable end ncsu_smd </pre></div> </td></tr></table>

图3 Ace-(Pro)5-NHe肽的从全反式到全顺式的构象的拉伸动力学输入文件(&cntrl参数没有显示)

向前和向后拉伸的非平衡功数值可以根据Crooks关系式[5]用来估算Δf值, 估算的方程与早期Bennet的接受率公式一致, 公式如下(1)式所示:

$$\Sum_{i=1}^{n_F} {1 \over 1+n_F/n_R \exp(Wi^F-\D F)} -\Sum{i=1}^{n_R} { 1 \over 1+n_R/n_F \exp(W_i^R+\D F) }=0 \tag 1$$

$n_F, n_R$ 分别指向前和向后模拟的数字序号, 我们计算得到的Δf值为4.58kcal/mol. 将这个估计量与Jarzynski恒等式计算的只有向前过程的10.65kcal/mol和只有向后过程的-1.49kcal/mol相比较是很有意义的.

(哈密顿)副本交换

副本交换方法是一个强有力的构想, 它由Geyer在1991年的会议论文集中提出[1], 随后在不同的科学领域被重新发现[10]. 该技术的基本思想是在几种模拟之间交换构象, 以保持精细的平衡, 从而使每个模拟都能从其目标分布中取样. 一个副本(根据势垒被称为"冷"副本)中高的势垒可能在另一个副本("热"副本)中会明显变低, 从而跨越势垒的能力就大大加强了.

该方法的成功与否, 关键取决于"热"副本的选择. 因此我们将在"热"副本中使用平均力势(PMF)来偏置动力学, 这种方法优于传统的平行回火方法在于, 温度不会在整个复制过程中发生变化, 从而使这种设置更适合于显式的溶剂模拟. 我们(厚着脸皮地)向有兴趣的读者推荐我们的预印版论文来获得更加全面的介绍.

这个想法很简单: 假设一个"慢模式"可以由一个集合变量σ(r)来描述, 也就是说与这个集合变量有关的平均力势(也称为朗道自由能)是已知的. 我们通过对PMF取负值来偏置最初的动力学, 可以建立一个"热"副本, 因而能够致使不同σ(r)值对应的状态是等概率出现的. 实际上, PMF通常是未知的, 但幸运的是, 偏置的势能在几个kBT的实际PMF值之内, 所以是可行的. PMF值可以由各种各样的方法得到. 本文中我们将使用所谓的ABMD11, 它是早期的局部高程方法12的另一个变体, 我们努力使得所谓的meta动力学方法13,14更加符合实际, 从而发展出这种方法.

对于Ace-(Pro)5-NHe肽, 它很容易辨别出肽键的反式/顺式构象, 并转变为"慢模式". 一共有5个肽键, 所以我们目标是建立一个被PMFs以及相应的二面角值偏置过的有五个"热"副本的副本-交换组合, 按比例缩小的PMFs还可以偏置出五个副本(作为在副本之间进行更有效的随机游走的代理), 还有一个未被偏置的副本, 可以从平衡态样品中"读取". 我们使用ABMD计算(近似的)PMFs值. 首先, 粗略的近似(为要求更高的副本-交换操作提供初始条件), 我们假设为Ace-(Pro)5 -NHe的:1@CH3 == :1@C == :2@N == :2@CA二面角计算出的PMF与为Ace-(Pro)5-NHe计算出的PMF, 在Pro残基之间的二面角计算的PMF值与Ace-(Pro)2 -NHe二聚体计算的相近.

我们从单体开始, 进行两个阶段的"漂移", 从更加有侵略性的开始(不平衡), 之后会更平坦(接近于平衡). 请在帮助文档中的4.hremd/1.preliminaries/0.monomer子目录寻找相应的文件. "漂移"这里指的是"na_ve", 线性增长的偏置势旨在补偿真正的PMF. 详情请参阅我们的预印本. 最后, 我们通过做短平衡偏置运行来检验PMF的质量(文件在0.monomer/3.umbrella文件夹下). 事实证明, 角度确实在短期内历经了所有可能的值. 这就意味着, 偏置势(如预期的那样)是在真正的PMF的几个kBT范围内的. 如果事实并非如此, 我们将在需要的时候用尽可能更慢的漂移.

下一步我们以相同的步骤计算Ace-(Pro)2-NHe二聚体:2@CA == :2@C == :3@N == :3@CA二面角的PMF值(文件在4.hremd/1.preliminaries/0.dimer文件夹下).

为了平衡初始副本-交换设置和检查这些用较小分子近似计算的PMF值是否接近真实值, 误差是否小于几个kBT, 我们运行了一个短副本-交换模拟, 包含十个副本(五个充分偏置过的和五个代理副本). 这些文件位于4.hremd/2.hremd.eq文件夹下, 我们鼓励读者去检查他们. 我们利用熟悉的multisander原理, 还使用mdin文件中ncsu_bbmd部分激活的另一个(非主流)代码路径. 一切都应该是不言自明的(如果不是的话, 可以通过邮件联系我). 从这段短的运行结果来看, 最初的PMF与真实值不是很相近, 由于二面角的值分布不是很均匀. 要纠正这个问题, 我们使用平稳的"代理"偏执势进行了一次"漂移" (4. hremd / 3. hremd.flooding). 然后我们通过短的平衡模拟(大概)检查了更准确PMF值(4. hremd/4. hremd.umbrella), 之后进入生产阶段(以防我们不满意PMF的质量,我们将做更多的缓慢的"漂移").

生产阶段的运行文件在4.hremd/5.hremd.production文件夹下, 它包含11个副本(五个"热"副本, 五个"代理"副本和一个"冷"副本), 每个运行50ns. 我们尝试每91步交换5个随机的副本. 交换的统计数据在exchange_log_file文件中, 可以看到结果非常好.

在完成"生产"运行后, 我们得到了" (1)所有11个副本的监控文件: 10个持有偏置的二面角样本, 带有11个COS_OF_DIHEDRAL集合变量个平衡样本; (2)没有偏置的构象, 保存为来自11个副本的轨迹, 我们将分析轨迹, 在传递中指出可以用来提高PMF以及相关二面角的精度的偏置样本. 我们使用ptraj从轨迹中提取变化慢的角的值. 可以看出, 这些角度确实非常集中于顺式和反式值, 因此它们可以被明确地分为顺式或反式. 因此, 使用由T (对于反式)和C (对于顺式)组成的五字母串来描述每个轨迹框架是很自然的. 一个简单的perl代码(补充文档的脚本子目录下的omegas.pl文件)用来比较不同构象出现的频率. 我们发现TTTTT(全反式)确实是最有可能出现的构象(全局最小值), 它的概率有33%. 我们也可以计算TTTTT和CCCCC状态之间的自由能差(可以根据设定COS_OF_DIHEDRAL值分别为-5和+5). 我们得到这个自由能差为4.42kcal/mol, 这与上述拉伸动力学描述的结果相吻合.

参考文献

  1. C. J. Geyer. Markov chain monte carlo maximum likelihood. In Computing Science and Statistics: The 23rd symposium on the interface, pages 156–163, Fairfax, 1991. Interface Foundation.
  2. A. Onufriev, D. Bashford, and D. A. Case. Modification of the generalized Born model suitable for macromolecules. J. Phys. Chem. B, 104:3712–3720, 2000.
  3. A. Onufriev, D. Bashford, and D. A. Case. Exploring protein native states and large-scale conformational changes with a modified generalized Born model. Proteins, 55:383–394, 2004.
  4. J. Weiser, P. S. Shenkin, and W. C. Still. Approximate Atomic Surfaces from Linear Combinations of Pairwise Overlaps (LCPO). J. Comp. Chem., 20:217–230, 1999.
  5. Gavin E. Crooks. Path-ensemble averages in systems driven far from equilibrium. Phys. Rev. E, 61(3):2361–2366, Mar 2000.
  6. Michael R. Shirts, Eric Bair, Giles Hooker, and Vijay S. Pande. Equilibrium free energies from nonequilibrium measurements using maximum-likelihood methods. Phys. Rev. Lett., 91(14):140601, 2003.
  7. Ioan Kosztin, Bogdan Barz, and Lorant Janosi. Calculating potentials of mean force and diffusion coefficients from nonequilibrium processes without jarzynski’s equality. J. Chem. Phys., 124(6):064106, 2006.
  8. C. H. Bennett. Efficient estimation of free energy differences from monte carlo data. J. Comp. Phys., 22:245–268, 1976.
  9. C. Jarzynski. Nonequilibrium equality for free energy differences. Phys. Rev. Lett., 78(14):2690–2693, 1997.
  10. K Hukushima and K Nemoto. Exchange monte carlo method and application to spin glass simulations. Journal of the Physical Society of Japan, 65:1604–1608, 1996.
  11. Volodymyr Babin, Christopher Roland, and Celeste Sagui. Adaptively biased molecular dynamics for free energy calculations. J. Chem. Phys., 128(13):134101, 2008.
  12. Thomas Huber, Andrew E. Torda, and Wilfred F. van Gunsteren. Local elevation: a method for improving the searching properties of molecular dynamics simulation. J. Comput. Aided. Mol. Des., 8:695–708, 1994.
  13. A. Laio and M. Parrinello. Escaping free-energy minima. Proc. Natl. Acad. Sci., 99:12562–12566, 2002.
  14. M. Iannuzzi, A. Laio, and M. Parrinello. Efficient exploration of reactive potential energy surfaces using car-parrinello molecular dynamics. Phys. Rev. Lett., 90:238302–1, 2003.