一種穩(wěn)定蛋白質(zhì)的全部原子結(jié)構(gòu)預(yù)測(cè)和折疊模擬

1、AMBER教程8：研究案例一種穩(wěn)定蛋白質(zhì)的全部原子結(jié)構(gòu)預(yù)測(cè)和折疊模擬這段教程展示的是一個(gè)研究實(shí)例，像您演示如何重現(xiàn)下述文章中的研究工作： Simmerling, C., Strockbine, B., Roitberg, A.E., J. Am. Chem. Soc., 2002, 124, 11258-11259 (/10.1021/ja0273851)我們建議您在開(kāi)始本教程前首先閱讀上述文章，獲得該蛋白的氨基酸序列及其他有用信息。警告1: 本教程中的一些計(jì)算耗時(shí)很長(zhǎng)，我使用了由16個(gè)1.3GHz cup的SGI Altix進(jìn)行了27小時(shí)計(jì)算才完成整個(gè)工作，因

2、此如果您沒(méi)有足夠的計(jì)算能力，我強(qiáng)烈建議您在重復(fù)本教程的過(guò)程中使用我為您提供的out文件，以使得您能夠流暢地完成整個(gè)教程。警告2: 如果您重復(fù)本教程，我們并不能保證您能夠精確地重現(xiàn)我的計(jì)算結(jié)果，在計(jì)算過(guò)程中，不同結(jié)構(gòu)的計(jì)算機(jī)會(huì)產(chǎn)生不同的近似誤差，從而使得計(jì)算過(guò)程搜索的是相空間的不同部位，但是模擬的平均結(jié)果是大致相同的。另外，盡管您完全重復(fù)了本教程也有可能無(wú)法獲得論文中給出的結(jié)果，而且即便是我們自己也無(wú)法保證論文中的結(jié)果能夠重現(xiàn)，這可能是因?yàn)槲夷M的時(shí)間不夠長(zhǎng)，獲取的僅僅是一個(gè)局部最小點(diǎn)，但是盡管如此，本教程的工作還是展示了蛋白折疊中一些有趣的行為。背景這篇論文應(yīng)用AMBER FF99力場(chǎng)和經(jīng)典的

3、全原子動(dòng)力學(xué)對(duì)一個(gè)肽的折疊過(guò)程進(jìn)行了模擬。模擬的對(duì)象"trpcage"是一個(gè)由20個(gè)氨基酸構(gòu)成的小肽，華盛頓大學(xué)的 Andersen已經(jīng)對(duì)這個(gè)蛋白做過(guò)了結(jié)構(gòu)優(yōu)化，它是現(xiàn)在已知最小的能夠顯示兩種不同折疊狀態(tài)的蛋白，而且這個(gè)蛋白在室溫下可以穩(wěn)定存在。該蛋白的小身量使得它成為模擬蛋白質(zhì)折疊的絕嘉對(duì)象。當(dāng)最早的關(guān)于這個(gè)蛋白的折疊的計(jì)算結(jié)果出爐時(shí)，對(duì)這個(gè)蛋白結(jié)構(gòu)的實(shí)驗(yàn)測(cè)定還沒(méi)有完成，所以整個(gè)模擬過(guò)程是在沒(méi)有實(shí)驗(yàn)數(shù)據(jù)作為指導(dǎo)的情況下完成的。當(dāng)?shù)鞍椎慕Y(jié)構(gòu)經(jīng)由實(shí)驗(yàn)手段測(cè)定之后，人們驚喜地發(fā)現(xiàn)，計(jì)算機(jī)模擬的結(jié)果與實(shí)驗(yàn)測(cè)定的數(shù)值之間的RMSD值僅為1.4A?？紤]到整個(gè)模擬過(guò)程是從蛋白的一級(jí)結(jié)構(gòu)

4、開(kāi)始并且完全沒(méi)有同源蛋白作為參考，這樣的一個(gè)計(jì)算結(jié)果是非常精確的。本教程中，我們?cè)噲D重復(fù)論文中的結(jié)果，計(jì)算的設(shè)定都與論文非常接近，只是由于計(jì)算能力的限制，在教程中我們只進(jìn)行一個(gè)50ns級(jí)的模擬。這已經(jīng)足夠重見(jiàn)蛋白質(zhì)折疊的結(jié)果了。在這里必須提醒的是，由于模擬過(guò)程的長(zhǎng)度所限，在不同的計(jì)算機(jī)，或在處理器數(shù)量不同的情況下，計(jì)算的結(jié)果將會(huì)是不同的。這是由分子動(dòng)力學(xué)模擬的方法決定的，實(shí)施過(guò)程的細(xì)微變化或者浮點(diǎn)計(jì)算中舍入的變化都意味著由不同的計(jì)算機(jī)進(jìn)行采樣的動(dòng)力學(xué)軌跡會(huì)隨著時(shí)間的流逝發(fā)生不可預(yù)知的分化。這并非誤差或者程序的bug，也并不意味著某一個(gè)模擬過(guò)程比其他的過(guò)程更合理。這僅僅意味著不同的模擬過(guò)程搜索的

5、是相空間的不同區(qū)域，如果我們平均一下模擬的結(jié)果，或者運(yùn)行更長(zhǎng)時(shí)間的動(dòng)力學(xué)過(guò)程，我們會(huì)在不同的機(jī)器上得到完全相同的結(jié)果，他們之間僅僅在過(guò)程上有所不同。因而我們說(shuō)在本教程中我們很難精確的再現(xiàn)論文中的結(jié)果，但是我們?cè)噲D重新創(chuàng)造那個(gè)重要的結(jié)果，即用AMBER程序來(lái)預(yù)測(cè)一個(gè)20氨基酸的小蛋白的空間結(jié)構(gòu)是可以完成的。那么記住這一點(diǎn)，讓我們開(kāi)始吧第一步：構(gòu)建起始結(jié)構(gòu)在以往的教程中，我們要么有一個(gè)可用的晶體結(jié)構(gòu)，要么可以通過(guò)程序生成一個(gè)已經(jīng)初步優(yōu)化的結(jié)構(gòu)。而在這個(gè)教程中我們要用的結(jié)構(gòu)太復(fù)雜，沒(méi)法通過(guò)手畫(huà)的辦法完成，同時(shí)我們也沒(méi)有一個(gè)可用的PDB結(jié)構(gòu)，因此我們就需要構(gòu)建一個(gè)線形的肽鏈，非常幸運(yùn)的是，在LEAP中

6、有一個(gè)命令可以完成這個(gè)工作，就是 sequence。蛋白的一級(jí)結(jié)構(gòu)序列在所列論文中可以查到，如下所示:NLYIQWLKDGGPSSGRPPPS這是用單字母符號(hào)顯示的蛋白質(zhì)一級(jí)結(jié)構(gòu)序列，在Leap中使用之前我們需要將其轉(zhuǎn)換成標(biāo)準(zhǔn)的三字母表示下面的表格給出了單字母表示和三字母表示之間的轉(zhuǎn)換關(guān)系：?jiǎn)巫帜概c三字母的轉(zhuǎn)換 conversionGPAVLIMCFYWHKRQNEDST甘氨酸 (Gly)脯氨酸 (Pro)丙氨酸(Ala)纈氨酸(Val)亮氨酸 (Leu)異亮氨酸 (Ile)蛋氨酸 (Met)半胱氨酸 (Cys)苯丙氨酸 (Phe)酪氨酸 (Tyr)色氨酸 (Trp)組氨酸 (His)賴氨酸

7、(Lys)精氨酸 (Arg)谷氨酸鹽 (Glu)天冬酰氨 (Asn)谷氨酸 (Glu)天冬氨酸 (Asp)絲氨酸 (Ser)蘇氨酸 (Thr)那么上述序列可以轉(zhuǎn)寫(xiě)為：ASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO SER但是這還沒(méi)有結(jié)束，LEaP 不能自動(dòng)識(shí)別序列的兩端，所以我們必須手工為這個(gè)序列標(biāo)定N末端和C末端，標(biāo)定的方法就是在N末端氨基酸前方加上N，C末端氨基酸前方加上字母C。最終在 LEaP中使用的序列如下：NASN LEU TYR ILE GLN TRP LEU LYS ASP G

8、LY GLY PRO SER SER GLY ARG PRO PRO PRO CSER下面啟動(dòng)xleap并調(diào)用ff99力場(chǎng):    $AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99(使用xLeap的時(shí)候一定要記住要關(guān)閉Num Lock鍵！否則工具欄會(huì)無(wú)法使用)下面使用 sequence 命令來(lái)建立蛋白的起始結(jié)構(gòu) (如需了解 sequence 命令的詳細(xì)情況可以在Leap中鍵入: help sequence). 注意： 為了版面設(shè)計(jì)的需要，下面將命名分為三行顯示，實(shí)際上您必須將所有內(nèi)容在一行

9、內(nèi)輸入，其間不能回車。    >TC5b = sequence NASN LEU TYR ILE GLN TRP LEU LYS                        ASP GLY GLY PRO SER SER GLY ARG        

10、                PRO PRO PRO CSER 我們需要的起始結(jié)構(gòu)就放在 TC5b中我們可以使用 edit命令來(lái)觀察這個(gè)結(jié)構(gòu)。    >edit TC5b現(xiàn)在我們獲得了一個(gè)線形的蛋白質(zhì)序列作為起始結(jié)構(gòu)，但是在這個(gè)起始結(jié)構(gòu)中很多原子是相互抵觸的，所以在進(jìn)行分子動(dòng)力學(xué)模擬之前我們要對(duì)這個(gè)結(jié)構(gòu)首先進(jìn)行短時(shí)間的優(yōu)化。我們暫時(shí)將Unit中的這個(gè)結(jié)構(gòu)存成一個(gè).lib文件，這樣在之后的操作中，我們只要調(diào)用這

11、個(gè)lib就可以簡(jiǎn)單地取出起始結(jié)構(gòu)，同時(shí)我們還要將這個(gè)結(jié)構(gòu)存成一個(gè)PDB文件，以便直觀地進(jìn)行觀察。    >saveoff TC5b TC5b_linear.lib    >savepdb TC5b TC5b_linear.pdb(TC5b_linear.lib, TC5b_linear.pdb)第二步：創(chuàng)建prmtop和inpcrd文件我們已經(jīng)有了起始結(jié)構(gòu)，下一步的工作是創(chuàng)建prmtop以及inpcrd文件。在進(jìn)行這一步之前我們需要首先確認(rèn)我們使用的參數(shù)和文獻(xiàn)中報(bào)道的是一樣的，在論文的第三段講到：We initiated

12、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 SANDE

13、R 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 energ

14、ies7 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.文獻(xiàn)顯示，他們首先建立了一個(gè)線形的起始結(jié)構(gòu)，這一步我們已經(jīng)完成了，之后他們運(yùn)行了沒(méi)有限制的恒溫動(dòng)力學(xué)模擬過(guò)程（即正則系綜中的模擬）在動(dòng)力學(xué)過(guò)程中他們使用了廣義波恩近似來(lái)模擬溶劑效應(yīng)的影響。AMBER程序可以支持很多不同的廣

15、義波恩模型，在這些模型中最先進(jìn)的是由A. Onufriev, D. Bashford 和 D.A. Case等人開(kāi)發(fā)的改良GB模型，這個(gè)GB模型使用了模型II的半徑 (IGB=5) 具體可以參考AMBER用戶手冊(cè)的GB模型一章。在論文中，他們沒(méi)有使用特殊的GB模型，這是因?yàn)樵谀菚r(shí)AMBER程序中只有IGB=1這個(gè)設(shè)定可用。為了使我們的教程盡可能接近文獻(xiàn)報(bào)道，我們也使用IGB=1的設(shè)定。由于Leap默認(rèn)的設(shè)定就是IGB=1，所以我們無(wú)需專門(mén)對(duì)此作出設(shè)定。論文中還聲明他們使用了FF99力場(chǎng)，這與我們之前設(shè)定的是一樣的，但是他們的立場(chǎng)有改進(jìn)的 phi/psi二面角參數(shù)，這是對(duì)FF99立場(chǎng)中phi/p

16、si二面角參數(shù)的一種校正，可以更好的模擬蛋白質(zhì)中 alpha螺旋的結(jié)構(gòu)。為了更好地重復(fù)文獻(xiàn)中的工作，我們需要建立一個(gè)包含上述修正的參數(shù)文件。但是比較麻煩的是，文獻(xiàn)中并沒(méi)有明確給出那些參數(shù)做了如何的改變，僅僅給出了一個(gè)修正后的parm99.dat文件。出現(xiàn)這種情況的原因我認(rèn)為可能是 AMBER6本身不帶FF99力場(chǎng)，在當(dāng)時(shí)存在很多不同的版本，文獻(xiàn)的作者為了讓人們了解他們使用的是官方版本的FF99力場(chǎng)所以在論文中展示了 parm99.dat文件。但很不幸，ACS以PDF文件格式給出了這個(gè)文件，這使得我們很難直接使用這個(gè)文件。幸運(yùn)的是，在AMBER8版本中，給出了這個(gè)修正的力場(chǎng)，位于下述路徑：$AM

17、BERHOME/dat/leap/parm/ as frcmod.mod_phipsi.1. 下面我也列出了文件的內(nèi)容，以備不時(shí)：frcmod.mod_phipsi.1from Simmerling, Strockbine, Roitberg, JACS 124:11258, 2002. Modifies parm99.MASSBONDANGLDIHEDRALN -CT-C -N 1 0.700 180.000 -1.N -CT-C -N 1 1.100 180.000 2.C -N -CT-C 1 1.000 0.000 1.NONB如你所見(jiàn)，只有三個(gè)二面角參數(shù)發(fā)生了變化，所以我們只需要打開(kāi)

18、xLEaP，讀取這個(gè)文件，其中的參數(shù)就會(huì)自動(dòng)覆蓋原有的參數(shù)。如果現(xiàn)在你已經(jīng)關(guān)閉了xLEaP，可以重新打開(kāi)并調(diào)入蛋白質(zhì)結(jié)構(gòu)：     $AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99    >loadoff TC5b_linear.lib然后調(diào)入修正的二面角參數(shù)：    >loadamberparams frcmod.mod_phipsi.1下面就可以存儲(chǔ)我們的 prmtop和 inpcrd文件:  

19、  >saveamberparm TC5b TC5b.prmtop TC5b.inpcrd下面是生成的輸入文件 TC5b.prmtop, TC5b.inpcrd第三步：預(yù)優(yōu)化蛋白質(zhì).在運(yùn)行分子動(dòng)力學(xué)模擬之前我們必須對(duì)起始結(jié)構(gòu)進(jìn)行短時(shí)間的優(yōu)化，這樣的話我們的體系就不會(huì)因?yàn)榫植磕芰康木奂趧?dòng)力學(xué)過(guò)程中出現(xiàn)問(wèn)題。結(jié)構(gòu)優(yōu)化的過(guò)程會(huì)整理整個(gè)分子結(jié)構(gòu)，重新修整氫的位置，這樣的過(guò)程之后我們的動(dòng)力學(xué)模擬會(huì)比較穩(wěn)定 。下面是結(jié)構(gòu)優(yōu)化的輸入文件:min1.inStage 1 - minimisation of TC5b &cntrl imin=1, maxcyc=1000, ncyc=5

20、00, cut=999., rgbmax=999.,igb=1, ntb=0, ntpr=100 /我們總共運(yùn)行1000步優(yōu)化過(guò)程，其中500步為最陡下降法(ncyc=500)，然后緊跟500步共軛梯度法(maxcyc-ncyc)。這樣的設(shè)置已經(jīng)足夠充分地釋放聚集在起始結(jié)構(gòu)中的能量。需要提醒的是我在輸入文件中設(shè)置了非常大的截?cái)嘀?cut=999. angstroms)，這樣設(shè)置是因?yàn)槲覀兪褂昧朔侵芷谀M(ntb=0) ，故而我們沒(méi)有使用PME方法，也就不會(huì)出現(xiàn)長(zhǎng)程的靜電相互作用。如果使用了 PME，推薦的截?cái)嘀凳?埃，在這樣一個(gè)范圍內(nèi)實(shí)際上模擬的主要是范德華相互作用。 但是如果不使用PME而設(shè)

21、置截?cái)嘀档脑?，范德華相互作用和靜電相互作用都在截?cái)嘀档姆秶鷥?nèi)被截?cái)嗔?，所以在沒(méi)有使用PME方法的狀況下，最好要將截?cái)嘀当M可能設(shè)大。需要提醒的是，模擬的耗時(shí)是與截?cái)嘀档钠椒匠烧鹊?，（參?jiàn) 教程一 的 5.1.2節(jié)）所幸我們模擬的體系非常小，足夠承受沒(méi)有截?cái)嘀?cut=999)的計(jì)算?；谕瑯拥脑砦覀儗gbmax也設(shè)置為 999埃，這個(gè)參數(shù)控制了在計(jì)算非鍵相互作用過(guò)程中列用于計(jì)算有效波恩半徑的粒子對(duì)的最遠(yuǎn)間距。這個(gè)值設(shè)定的越大，計(jì)算的結(jié)果就越好，當(dāng)然也就需要花費(fèi)越多的計(jì)算時(shí)間?？紤]到我們面對(duì)的體系只有20個(gè)氨基酸殘基我們可以把所有的粒子都納入到有效波恩半徑中來(lái)，所以我們?cè)O(shè)定的rgbmax遠(yuǎn)遠(yuǎn)

22、大于計(jì)算的尺度。.下面開(kāi)始運(yùn)行優(yōu)化過(guò)程:    $AMBERHOME/exe/sander -O -i min1.in -o min1.out -p TC5b.prmtop                             -c TC5b.inpcrd -r min1.rstInpu

23、t Files: TC5b.prmtop, TC5b.inpcrd, min1.inOutput Files: min1.out, min1.rst在16個(gè)1.3GHzCPU的 SGI Altix上這個(gè)過(guò)程需要3.5秒完成為了直觀的比較優(yōu)化前后的結(jié)構(gòu)，我們生成一個(gè)pdb文件:    $AMBERHOME/exe/ambpdb -p TC5b.prmtop < min1.rst > min1.pdb將優(yōu)化前后的兩個(gè)文件打開(kāi)(min1.pdb and TC5b_linear.pdb)你可以選擇任何可用的顯示軟件，比如VMD起始結(jié)構(gòu)用藍(lán)色顯示，優(yōu)化后的結(jié)構(gòu)用

24、黃色顯示。如你所見(jiàn)，優(yōu)化過(guò)程并未造成主鏈結(jié)構(gòu)太大的變化但是色氨酸和酪氨酸殘基發(fā)生了比較明顯的移動(dòng)，這些能量熱點(diǎn)集中的區(qū)域有可能在我們開(kāi)始分子動(dòng)力學(xué)模擬之后帶來(lái)麻煩，如果你不相信，可以用未經(jīng)優(yōu)化的結(jié)構(gòu)跑一個(gè)動(dòng)力學(xué)過(guò)程看看，肯定飛！第四步：體系加熱.接下來(lái)我們要在這個(gè)體系中正式開(kāi)始分子動(dòng)力學(xué)模擬，首先我們要分7步花費(fèi)50ps時(shí)間對(duì)體系進(jìn)行升溫模擬。將升溫過(guò)程分為7步完成可以在每一步升溫之后維持一段時(shí)間，以免一次升溫造成體系能量聚集最終跑飛，另外一種可行的方法是對(duì)升溫過(guò)程加一個(gè)權(quán)重限制。您可以參閱AMBER用戶手冊(cè)以獲取更多信息。一般而言我們升溫的最終目標(biāo)是室溫即300K但是為了重復(fù)文獻(xiàn)的運(yùn)算我們選

25、擇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.文獻(xiàn)認(rèn)為必須將體系加溫到325K進(jìn)行模擬，否則有可能使模擬的結(jié)果最終

26、落入局部最小點(diǎn)，所以我們也做同樣的設(shè)定。但是你必須牢記更高的模擬溫度會(huì)導(dǎo)致體系中各化學(xué)鍵發(fā)生更加顯著的振動(dòng)，這意味著如果你打算做一個(gè)600K，以2fs為步長(zhǎng)的動(dòng)力學(xué)模擬，你就要考慮一個(gè)應(yīng)用了shaken的300k效果會(huì)與之相同，但600K的模擬卻要臨步長(zhǎng)過(guò)大的問(wèn)題，過(guò)大的步長(zhǎng)會(huì)導(dǎo)致體系不穩(wěn)定。還好325K不算太高，還比較接近常用的300K，2fs的步長(zhǎng)可以處理含氫的鍵的振動(dòng)?？墒羌偃缥覀円?00K的條件下運(yùn)行動(dòng)力學(xué)模擬的話，那模擬的步長(zhǎng)就要縮減到1.5fs。我們的起始結(jié)構(gòu)是手工搭建的，不是通常常見(jiàn)的來(lái)自實(shí)驗(yàn)的晶體結(jié)構(gòu)，所以我們的體系在模擬的開(kāi)始階段要面臨不如晶體結(jié)構(gòu)穩(wěn)定的問(wèn)題。為了讓我們的體

27、系能夠在可控制的狀況下來(lái)釋放能量，在模擬起始的升溫階段我們選擇步長(zhǎng)為0.5fs，進(jìn)入相對(duì)穩(wěn)定的生成相之后，我們?cè)龠x擇常規(guī)的2fs步長(zhǎng)0.5fs的步長(zhǎng)確實(shí)有些矯枉過(guò)正，但是保證體系的安全畢竟還是最重要的。我們進(jìn)行升溫模擬的方案如下：    第一階段 - 10,000 步, 步長(zhǎng) 0.5fs (共5ps), 起始溫度 0.0K, 結(jié)束溫度 50.0K, 溫度耦合系數(shù) 1.0ps    第二階段 - 10,000 步, 步長(zhǎng) 0.5fs (共5ps), 結(jié)束溫度 100.0K, 溫度耦合系數(shù) 1.0ps   

28、第三階段 - 10,000 步, 步長(zhǎng) 0.5fs (共5ps), 結(jié)束溫度 150.0K, 溫度耦合系數(shù) 1.0ps    第四階段 - 10,000 步, 步長(zhǎng) 0.5fs (共5ps), 結(jié)束溫度 200.0K, 溫度耦合系數(shù) 1.0ps    第五階段 - 10,000 步, 步長(zhǎng) 0.5fs (共5ps), 結(jié)束溫度 250.0K, 溫度耦合系數(shù) 1.0ps    第六階段 - 10,000 步, 步長(zhǎng) 0.5fs (共5ps), 結(jié)束溫度 300.0K, 溫度耦合系數(shù) 1.0ps &

29、#160;  第七階段 - 40,000 步, 步長(zhǎng) 0.5fs (共20ps),結(jié)束溫度 325.0K, 溫度耦合系數(shù) 1.0ps下面是第一階段的輸入文件：heat1.inStage 1 heating of TC5b 0 to 50K &cntrl imin=0, irest=0, ntx=1, nstlim=10000, dt=0.0005, ntc=2, ntf=2, ntt=1, tautp=1.0, tempi=0.0, temp0=50.0, ntpr=50, ntwx=50, ntb=0, igb=1, cut=999.,rgbmax=999. /其他六個(gè)階段

30、的輸入文件與之非常接近，只需要改變一下相應(yīng)的溫度就可以了，可以從此處下載現(xiàn)成的輸入文件： (heat2.in, heat3.in, heat4.in, heat5.in, heat6.in, heat7.in).下面是一個(gè)運(yùn)行升溫模擬的PBS腳本，你也可以根據(jù)你的系統(tǒng)自己寫(xiě)一個(gè)腳本。#PBS -l ncpus=16#PBS -l walltime=500:00:00#PBS -l cput=2000:00:00#PBS -j oesetenv AMBERHOME /usr/people/rcw/amber9cd rcw/initial_heatingmpirun -np 16 $AMBERHO

31、ME/exe/sander -O -i heat1.in -p TC5b.prmtop -c min1.rst -r heat1.rst -o heat1.out -x heat1.mdcrdgzip -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

32、.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 he

33、at4.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

34、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

35、-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 he

36、at5.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.mdcrd gzip -9 heat7.mdcrdmkdir initial_heatingcp heat1.out initia

37、l_heatingcp heat2.out initial_heatingcp heat3.out initial_heatingcp heat4.out initial_heatingcp heat5.out initial_heatingcp heat6.out initial_heatingcp heat7.out initial_heatingcp heat1.mdcrd.gz initial_heatingcp heat2.mdcrd.gz initial_heatingcp heat3.mdcrd.gz initial_heatingcp heat4.mdcrd.gz initia

38、l_heatingcp heat5.mdcrd.gz initial_heatingcp heat6.mdcrd.gz initial_heatingcp heat7.mdcrd.gz initial_heatingecho "DONE"在16個(gè)1.3GHz CPU的SGI Altix上，7個(gè)步驟全部完成共耗時(shí)7分鐘。下面提供了全部過(guò)程的輸出文件，你可以分別下載，也可以下載最終的一個(gè)壓縮文件。Heating StageOutput FileRestrt FileMdcrd FileStage 1heat1.outheat1.rstheat1.mdcrd.gzStage 2h

39、eat2.outheat2.rstheat2.mdcrd.gzStage 3heat3.outheat3.rstheat3.mdcrd.gzStage 4heat4.outheat4.rstheat4.mdcrd.gzStage 5heat5.outheat5.rstheat5.mdcrd.gzStage 6heat6.outheat6.rstheat6.mdcrd.gzStage 7heat7.outheat7.rstheat7.mdcrd.gzComplete file setheating.tar.gz (5.2 Mb)將軌跡文件用VMD打開(kāi)就可以看到在升溫過(guò)程中究竟發(fā)生了什么。你可以看

40、到體系隨著溫度升高開(kāi)始折疊，我們對(duì)這一階段的軌跡并不關(guān)心，觀看升溫過(guò)程主要的目的在于確認(rèn)整個(gè)升溫過(guò)程一切OK，沒(méi)有發(fā)生什么意外。下圖顯示了升溫過(guò)程結(jié)束后肽鏈的結(jié)構(gòu):從這個(gè)結(jié)構(gòu)我們看出，一些alpha螺旋已經(jīng)形成了，但是這個(gè)結(jié)構(gòu)距離最終的穩(wěn)定折疊構(gòu)像還有很長(zhǎng)的路要走。第五步：生產(chǎn)相動(dòng)力學(xué)模擬本教程分子模擬部分的最后一步是在325K條件下進(jìn)行一個(gè)時(shí)間非常長(zhǎng)的動(dòng)力學(xué)模擬。在文獻(xiàn)中他們做了50ns的動(dòng)力學(xué)模擬，但是實(shí)際上我們看到的結(jié)果在模擬進(jìn)行了20ns之后就已經(jīng)呈現(xiàn)在人們面前了，之后繼續(xù)進(jìn)行的30ns模擬的意義僅僅在于說(shuō)明之前的計(jì)算獲得的就是一個(gè)穩(wěn)定的結(jié)果。 盡管文獻(xiàn)的作者發(fā)現(xiàn)在模擬的最初520ns

41、中蛋白就已經(jīng)充分折疊，我們?cè)诒窘坛讨腥匀贿M(jìn)行50ns的動(dòng)力學(xué)計(jì)算，以重復(fù)文獻(xiàn)的報(bào)道。"Two independent simulations converged to essentially identical families of structures after 5 and 20 ns."我們將這個(gè)總時(shí)間長(zhǎng)度為50ns的模擬分為10個(gè)階段，每段5ns，這樣做是為了一旦系統(tǒng)崩潰，我們不會(huì)損失已經(jīng)進(jìn)行的所有工作。而且這樣分開(kāi)處理還可以保證每個(gè)輸出文件和軌跡文件的大小都適合處理。這10個(gè)階段的模擬我們會(huì)使用相同的輸入文件，文件如下所示:equil.inStage 2 equ

42、ilibration 1 0-5ns &cntrl imin=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. /每一個(gè)階段的模擬都會(huì)進(jìn)行250,000步(由nstlim的取值決定)，步長(zhǎng)為2fs(由dt的取值決定) 即總共進(jìn)行 5 ns的模擬。請(qǐng)注意我們?cè)谀M全過(guò)程中使用了SHAKE(ntc=2, ntf=2)，我們使用了

43、Berendsen 恒溫器來(lái)保證模擬過(guò)程中體系溫度恒定(ntt=1)，另外由于我們的體系已經(jīng)完成升溫，并且保持穩(wěn)定，所以我們選擇了耦合地更加緊密的恒溫器(tautp=0.5)，這個(gè)恒溫器的作用在于保持我們模擬的對(duì)象始終保持325K的溫度。如同文獻(xiàn)中所列，我們將模擬的溫度設(shè)定在325K，每隔500步書(shū)寫(xiě)一次輸出文件和軌跡文件. 過(guò)于頻繁地寫(xiě)入這些文件會(huì)產(chǎn)生非常巨大的文件，這是不容易處理的。按照我們進(jìn)行的設(shè)定，每5ns的模擬會(huì)產(chǎn)生一個(gè) 35 Mb的軌跡文件，對(duì)于我們的研究來(lái)說(shuō)，500步記錄一次已經(jīng)足夠了。下面是進(jìn)行生產(chǎn)相模擬的PBS 腳本，您可以根據(jù)您系統(tǒng)的狀況修改腳本。#PBS -l ncpus

44、=16#PBS -l walltime=500:00:00#PBS -l cput=8000:00:00#PBS -j oesetenv AMBERHOME /usr/people/rcw/amber9cd rcw/productionmpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c heat7.rst -r equil1.rst -o equil1.out -x equil1.mdcrdgzip -9 equil1.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -

45、i equil.in -p TC5b.prmtop -c equil1.rst -r equil2.rst -o equil2.out -x equil2.mdcrdgzip -9 equil2.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil2.rst -r equil3.rst -o equil3.out -x equil3.mdcrdgzip -9 equil3.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p

46、TC5b.prmtop -c equil3.rst -r equil4.rst -o equil4.out -x equil4.mdcrdgzip -9 equil4.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil4.rst -r equil5.rst -o equil5.out -x equil5.mdcrdgzip -9 equil5.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c

47、 equil5.rst -r equil6.rst -o equil6.out -x equil6.mdcrdgzip -9 equil6.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil6.rst -r equil7.rst -o equil7.out -x equil7.mdcrdgzip -9 equil7.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil7.rst -r

48、 equil8.rst -o equil8.out -x equil8.mdcrdgzip -9 equil8.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil8.rst -r equil9.rst -o equil9.out -x equil9.mdcrdgzip -9 equil9.mdcrdmpirun -np 16 $AMBERHOME/exe/sander -O -i equil.in -p TC5b.prmtop -c equil9.rst -r equil10.rst -

49、o equil10.out -x equil10.mdcrdgzip -9 equil10.mdcrdecho "DONE"在16個(gè)1.3GHz CPU的SGI Altix上全部的十段模擬共耗時(shí)27小時(shí)。下面是每一步模擬的結(jié)果輸出，您可以下載這些結(jié)果，但是您需要注意每一個(gè)軌跡文件在壓縮之后都有13Mb 左右你可以將軌跡文件載入到VMD等看圖軟件中，提醒一下，你要首先解壓縮這個(gè)軌跡文件，然后載入它，整個(gè)過(guò)程有可能需要花費(fèi)一天時(shí)間，然后你就可以欣賞整個(gè)動(dòng)力學(xué)過(guò)程了，如果你用飄帶模型來(lái)顯示的話會(huì)非常好看。下面是一個(gè)抓圖：第六步：分析模擬結(jié)果下面我們要開(kāi)始分析模擬的結(jié)果，我們要繪制

50、模擬過(guò)程的溫度、總能量、動(dòng)能、勢(shì)能變化曲線，這些信息可以從輸出文件中獲取，檢查這些數(shù)據(jù)也可以讓我們知道在動(dòng)力學(xué)模擬過(guò)程中有沒(méi)有發(fā)生什么異常狀況。我們可以看到，溫度的曲線非常平滑，并且維持在325K。動(dòng)能和勢(shì)能的曲線也一樣，會(huì)在平衡位置比較平滑。這些曲線任何的突刺都暗示了我們的模擬過(guò)程發(fā)生了一些意外，我們就需要看看是不是用了不好的起始結(jié)構(gòu)、太長(zhǎng)的動(dòng)力學(xué)步長(zhǎng)或者不合適的參數(shù)等等。為了從輸出文件中提取我們需要的數(shù)據(jù)，我們可以使用下述 perl腳本： process_mdout.perl。用下面的命令我們可以把10個(gè)輸出文件中的信息提取到一個(gè)文件中：mkdir analysiscd analysis.

51、/process_mdout.perl ./heat1.out ./heat2.out ./heat3.out ./heat4.out ./heat5.out    ./heat6.out ./heat7.out ./equil1.out ./equil2.out ./equil3.out ./equil4.out ./equil5.out ./equil6.out ./equil7.out ./equil8.out ./equil9.out ./equil10.out下面這個(gè)文件是運(yùn)行process_mdout.perl提取出來(lái)的信息： analysis.tar.gz (2.3 Mb)下面我們可以用一些圖形化的數(shù)據(jù)處理軟件如xmgr來(lái)分析動(dòng)力學(xué)數(shù)據(jù)。xmgr我們