一种稳定蛋白质的全部原子结构预测和折叠模拟

AMBER教程8：研究案例——一种稳定蛋白质的全部原子结构预测和折叠模拟这段教程展示的是一个研究实例，像您演示如何重现下述文章中的研究工作：Simmerling, C., Strockbine, B., Roitberg, A.E., J. Am. Chem. Soc., 2002, 124, 11258-11259(/10.1021/ja0273851)我们建议您在开始本教程前首先阅读上述文章，获得该蛋白的氨基酸序列及其他有用信息。

警告1: 本教程中的一些计算耗时很长，我使用了由16个1.3GHz cup的SGI Altix进行了27小时计算才完成整个工作，因此如果您没有足够的计算能力，我强烈建议您在重复本教程的过程中使用我为您提供的out文件，以使得您能够流畅地完成整个教程。

警告2: 如果您重复本教程，我们并不能保证您能够精确地重现我的计算结果，在计算过程中，不同结构的计算机会产生不同的近似误差，从而使得计算过程搜索的是相空间的不同部位，但是模拟的平均结果是大致相同的。

另外，尽管您完全重复了本教程也有可能无法获得论文中给出的结果，而且即便是我们自己也无法保证论文中的结果能够重现，这可能是因为我模拟的时间不够长，获取的仅仅是一个局部最小点，但是尽管如此，本教程的工作还是展示了蛋白折叠中一些有趣的行为。

背景这篇论文应用AMBER FF99力场和经典的全原子动力学对一个肽的折叠过程进行了模拟。

模拟的对象"trpcage"是一个由20个氨基酸构成的小肽，华盛顿大学的Andersen已经对这个蛋白做过了结构优化，它是现在已知最小的能够显示两种不同折叠状态的蛋白，而且这个蛋白在室温下可以稳定存在。

该蛋白的小身量使得它成为模拟蛋白质折叠的绝嘉对象。

当最早的关于这个蛋白的折叠的计算结果出炉时，对这个蛋白结构的实验测定还没有完成，所以整个模拟过程是在没有实验数据作为指导的情况下完成的。

当蛋白的结构经由实验手段测定之后，人们惊喜地发现，计算机模拟的结果与实验测定的数值之间的RMSD值仅为1.4A。

考虑到整个模拟过程是从蛋白的一级结构开始并且完全没有同源蛋白作为参考，这样的一个计算结果是非常精确的。

本教程中，我们试图重复论文中的结果，计算的设定都与论文非常接近，只是由于计算能力的限制，在教程中我们只进行一个50ns级的模拟。

这已经足够重见蛋白质折叠的结果了。

在这里必须提醒的是，由于模拟过程的长度所限，在不同的计算机，或在处理器数量不同的情况下，计算的结果将会是不同的。

这是由分子动力学模拟的方法决定的，实施过程的细微变化或者浮点计算中舍入的变化都意味着由不同的计算机进行采样的动力学轨迹会随着时间的流逝发生不可预知的分化。

这并非误差或者程序的bug，也并不意味着某一个模拟过程比其他的过程更合理。

这仅仅意味着不同的模拟过程搜索的是相空间的不同区域，如果我们平均一下模拟的结果，或者运行更长时间的动力学过程，我们会在不同的机器上得到完全相同的结果，他们之间仅仅在过程上有所不同。

因而我们说在本教程中我们很难精确的再现论文中的结果，但是我们试图重新创造那个重要的结果，即用AMBER程序来预测一个20氨基酸的小蛋白的空间结构是可以完成的。

那么记住这一点，让我们开始吧第一步：构建起始结构在以往的教程中，我们要么有一个可用的晶体结构，要么可以通过程序生成一个已经初步优化的结构。

而在这个教程中我们要用的结构太复杂，没法通过手画的办法完成，同时我们也没有一个可用的PDB结构，因此我们就需要构建一个线形的肽链，非常幸运的是，在LEAP中有一个命令可以完成这个工作，就是sequence。

蛋白的一级结构序列在所列论文中可以查到，如下所示:下面使用sequence 命令来建立蛋白的起始结构(如需了解sequence 命令的详细情况可以在Leap中键入: help sequence). 注意：为了版面设计的需要，下面将命名分为三行显示，实际上您必须将所有内容在一行内输入，其间不能回车。

>TC5b = sequence { NASN LEU TYR ILE GLN TRP LEU LYSASP GLY GLY PRO SER SER GLY ARGPRO PRO PRO CSER }我们需要的起始结构就放在TC5b中我们可以使用edit命令来观察这个结构。

>edit TC5b现在我们获得了一个线形的蛋白质序列作为起始结构，但是在这个起始结构中很多原子是相互抵触的，所以在进行分子动力学模拟之前我们要对这个结构首先进行短时间的优化。

我们暂时将Unit中的这个结构存成一个.lib文件，这样在之后的操作中，我们只要调用这个lib就可以简单地取出起始结构，同时我们还要将这个结构存成一个PDB文件，以便直观地进行观察。

>saveoff TC5b TC5b_linear.lib>savepdb TC5b TC5b_linear.pdb(TC5b_linear.lib, TC5b_linear.pdb)第二步：创建prmtop和inpcrd文件我们已经有了起始结构，下一步的工作是创建prmtop以及inpcrd文件。

在进行这一步之前我们需要首先确认我们使用的参数和文献中报道的是一样的，在论文的第三段讲到：We initiated our simulations using only the trpcage TC5b2 amino acid sequence(N20LYIQWLKDGGPSSGRPPPS39), with an extended initial conformation built by the LEaP module of AMBER version 6.0.4 All molecular dynamics (MD) simulations were fully unrestrained and carried out in the canonical ensemble using the SANDER module, which we modified to improve performance on the Linux/Intel PC cluster that was used for all calculations. The ff99 force field5 was employed, with the exception of [phi/psi] dihedral parameters which were refit6 (see Supporting Information) to improve agreement with ab initio relative energies7 of alanine tetrapeptide conformations. Parameters were not fit to data for the trpcage. Solvation effects were incorporated using the Generalized Born model,8 as implemented9 in AMBER.Input Files: TC5b.prmtop, TC5b.inpcrd, min1.inOutput Files: min1.out, min1.rst在16个1.3GHzCPU的SGI Altix上这个过程需要3.5秒完成为了直观的比较优化前后的结构，我们生成一个pdb文件:$AMBERHOME/exe/ambpdb -p TC5b.prmtop < min1.rst > min1.pdb将优化前后的两个文件打开(min1.pdb and TC5b_linear.pdb)你可以选择任何可用的显示软件，比如VMD起始结构用蓝色显示，优化后的结构用黄色显示。

如你所见，优化过程并未造成主链结构太大的变化但是色氨酸和酪氨酸残基发生了比较明显的移动，这些能量热点集中的区域有可能在我们开始分子动力学模拟之后带来麻烦，如果你不相信，可以用未经优化的结构跑一个动力学过程看看，肯定飞！第四步：体系加热.接下来我们要在这个体系中正式开始分子动力学模拟，首先我们要分7步花费50ps时间对体系进行升温模拟。

将升温过程分为7步完成可以在每一步升温之后维持一段时间，以免一次升温造成体系能量聚集最终跑飞，另外一种可行的方法是对升温过程加一个权重限制。

您可以参阅AMBER用户手册以获取更多信息。

一般而言我们升温的最终目标是室温即300K但是为了重复文献的运算我们选择325K：MD simulations of 100 ns were performed at 300 K, but all were kinetically trapped on this time scale, showing strong dependence on initial conditions and failing to converge to similar conformational ensembles. We therefore increased the temperature to 325 K.文献认为必须将体系加温到325K进行模拟，否则有可能使模拟的结果最终落入局部最小点，所以我们也做同样的设定。

但是你必须牢记更高的模拟温度会导致体系中各化学键发生更加显著的振动，这意味着如果你打算做一个600K，以2fs为步长的动力学模拟，你就要考虑一个应用了shaken的300k效果会与之相同，但600K的模拟却要临步长过大的问题，过大的步长会导致体系不稳定。

还好325K不算太高，还setenv AMBERHOME /usr/people/rcw/amber9cd ~rcw/initial_heatingmpirun -np 16 $AMBERHOME/exe/sander -O -i heat1.in -p TC5b.prmtop -c min1.rst -r heat1.rst -o heat1.out -x heat1.mdcrd gzip -9 heat1.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i heat2.in -p TC5b.prmtop -c heat1.rst -r heat2.rst -o heat2.out -x heat2.mdcrdgzip -9 heat2.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i heat3.in -p TC5b.prmtop -c heat2.rst -r heat3.rst -o heat3.out -x heat3.mdcrdgzip -9 heat3.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i heat4.in -p TC5b.prmtop -c heat3.rst -r heat4.rst -o heat4.out -x heat4.mdcrdgzip -9 heat4.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i heat5.in -p TC5b.prmtop -c heat4.rst -r heat5.rst -o heat5.out -x heat5.mdcrdgzip -9 heat5.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i heat6.in -p TC5b.prmtop -c heat5.rst -r heat6.rst -o heat6.out -x heat6.mdcrdgzip -9 heat6.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i heat7.in -p TC5b.prmtop -c heat6.rst -r heat7.rst -o heat7.out -x heat7.mdcrdgzip -9 heat7.mdcrdecho "DONE"译者提供的bash脚本如下：#!/bin/bash#heatingsander -O -i heat1.in -p TC5b.prmtop -c min1.rst -r heat1.rst -o heat1.out -x heat1.mdcrdgzip -9 heat1.mdcrdsander -O -i heat2.in -p TC5b.prmtop -c heat1.rst -r heat2.rst -o heat2.out -x heat2.mdcrdgzip -9 heat2.mdcrdsander -O -i heat3.in -p TC5b.prmtop -c heat2.rst -r heat3.rst -o heat3.out -x heat3.mdcrdgzip -9 heat3.mdcrdsander -O -i heat4.in -p TC5b.prmtop -c heat3.rst -r heat4.rst -o heat4.out -x heat4.mdcrdgzip -9 heat4.mdcrdsander -O -i heat5.in -p TC5b.prmtop -c heat4.rst -r heat5.rst -o heat5.out -x heat5.mdcrdgzip -9 heat5.mdcrdsander -O -i heat6.in -p TC5b.prmtop -c heat5.rst -r heat6.rst -o heat6.out -x heat6.mdcrdgzip -9 heat6.mdcrdsander -O -i heat7.in -p TC5b.prmtop -c heat6.rst -r heat7.rst -o heat7.out -x heat7.mdcrdgzip -9 heat7.mdcrdmkdir initial_heatingcp heat1.out initial_heatingStage 2 equilibration 1 0-5ns&cntrlimin=0, irest=1, ntx=5,nstlim=2500000, dt=0.002,ntc=2, ntf=2,ntt=1, tautp=0.5,tempi=325.0, temp0=325.0,ntpr=500, ntwx=500,ntb=0, igb=1,cut=999.,rgbmax=999./每一个阶段的模拟都会进行250,000步(由nstlim的取值决定)，步长为2fs(由dt的取值决定) 即总共进行 5 ns的模拟。