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

1、AMBER教程8：研究案例一種穩(wěn)定蛋白質(zhì)的全部原子結(jié)構(gòu)預測和折疊模擬這段教程展示的是一個研究實例，像您演示如何重現(xiàn)下述文章中的研究工作： 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小時計算才完成整個工作，因

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

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

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

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

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

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

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

9、SN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO CSER 我們需要的起始結(jié)構(gòu)就放在 TC5b中我們可以使用 edit命令來觀察這個結(jié)構(gòu)。 edit TC5b現(xiàn)在我們獲得了一個線形的蛋白質(zhì)序列作為起始結(jié)構(gòu)，但是在這個起始結(jié)構(gòu)中很多原子是相互抵觸的，所以在進行分子動力學模擬之前我們要對這個結(jié)構(gòu)首先進行短時間的優(yōu)化。我們暫時將Unit中的這個結(jié)構(gòu)存成一個.lib文件，這樣在之后的操作中，我們只要調(diào)用這個lib就可以簡單地取出起始結(jié)構(gòu)，同時我們還要將這個結(jié)構(gòu)存成一個PDB文件，以便直觀地進行觀察。 s

10、aveoff 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文件。在進行這一步之前我們需要首先確認我們使用的參數(shù)和文獻中報道的是一樣的，在論文的第三段講到：We initiated our simulations using only the trpcage TC5b2 amino acid sequence (N20LYIQWLKDGGPSSGRPPPS39), w

11、ith 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

12、 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. S

13、olvation effects were incorporated using the Generalized Born model,8 as implemented9 in AMBER.文獻顯示，他們首先建立了一個線形的起始結(jié)構(gòu)，這一步我們已經(jīng)完成了，之后他們運行了沒有限制的恒溫動力學模擬過程（即正則系綜中的模擬）在動力學過程中他們使用了廣義波恩近似來模擬溶劑效應的影響。AMBER程序可以支持很多不同的廣義波恩模型，在這些模型中最先進的是由A. Onufriev, D. Bashford 和 D.A. Case等人開發(fā)的改良GB模型，這個GB模型使用了模型II的半徑 (IGB=5) 具體可

14、以參考AMBER用戶手冊的GB模型一章。在論文中，他們沒有使用特殊的GB模型，這是因為在那時AMBER程序中只有IGB=1這個設(shè)定可用。為了使我們的教程盡可能接近文獻報道，我們也使用IGB=1的設(shè)定。由于Leap默認的設(shè)定就是IGB=1，所以我們無需專門對此作出設(shè)定。論文中還聲明他們使用了FF99力場，這與我們之前設(shè)定的是一樣的，但是他們的立場有改進的 phi/psi二面角參數(shù)，這是對FF99立場中phi/psi二面角參數(shù)的一種校正，可以更好的模擬蛋白質(zhì)中 alpha螺旋的結(jié)構(gòu)。為了更好地重復文獻中的工作，我們需要建立一個包含上述修正的參數(shù)文件。但是比較麻煩的是，文獻中并沒有明確給出那些參數(shù)做

15、了如何的改變，僅僅給出了一個修正后的parm99.dat文件。出現(xiàn)這種情況的原因我認為可能是 AMBER6本身不帶FF99力場，在當時存在很多不同的版本，文獻的作者為了讓人們了解他們使用的是官方版本的FF99力場所以在論文中展示了 parm99.dat文件。但很不幸，ACS以PDF文件格式給出了這個文件，這使得我們很難直接使用這個文件。幸運的是，在AMBER8版本中，給出了這個修正的力場，位于下述路徑：$AMBERHOME/dat/leap/parm/ as frcmod.mod_phipsi.1. 下面我也列出了文件的內(nèi)容，以備不時：frcmod.mod_phipsi.1from Simme

16、rling, 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如你所見，只有三個二面角參數(shù)發(fā)生了變化，所以我們只需要打開xLEaP，讀取這個文件，其中的參數(shù)就會自動覆蓋原有的參數(shù)。如果現(xiàn)在你已經(jīng)關(guān)閉了xLEaP，可以重新打開并調(diào)入蛋白質(zhì)結(jié)構(gòu)： $AMBERHOME/exe/xleap -s -f $AMBE

17、RHOME/dat/leap/cmd/leaprc.ff99 loadoff TC5b_linear.lib然后調(diào)入修正的二面角參數(shù)： loadamberparams frcmod.mod_phipsi.1下面就可以存儲我們的 prmtop和 inpcrd文件: saveamberparm TC5b TC5b.prmtop TC5b.inpcrd下面是生成的輸入文件 TC5b.prmtop, TC5b.inpcrd第三步：預優(yōu)化蛋白質(zhì).在運行分子動力學模擬之前我們必須對起始結(jié)構(gòu)進行短時間的優(yōu)化，這樣的話我們的體系就不會因為局部能量的聚集而在動力學過程中出現(xiàn)問題。結(jié)構(gòu)優(yōu)化的過程會整理整個分子結(jié)構(gòu)

18、，重新修整氫的位置，這樣的過程之后我們的動力學模擬會比較穩(wěn)定 。下面是結(jié)構(gòu)優(yōu)化的輸入文件:min1.inStage 1 - minimisation of TC5b &cntrl imin=1, maxcyc=1000, ncyc=500, cut=999., rgbmax=999.,igb=1, ntb=0, ntpr=100 /我們總共運行1000步優(yōu)化過程，其中500步為最陡下降法(ncyc=500)，然后緊跟500步共軛梯度法(maxcyc-ncyc)。這樣的設(shè)置已經(jīng)足夠充分地釋放聚集在起始結(jié)構(gòu)中的能量。需要提醒的是我在輸入文件中設(shè)置了非常大的截斷值(cut=999. angstrom

19、s)，這樣設(shè)置是因為我們使用了非周期模擬(ntb=0) ，故而我們沒有使用PME方法，也就不會出現(xiàn)長程的靜電相互作用。如果使用了 PME，推薦的截斷值是8埃，在這樣一個范圍內(nèi)實際上模擬的主要是范德華相互作用。 但是如果不使用PME而設(shè)置截斷值的話，范德華相互作用和靜電相互作用都在截斷值的范圍內(nèi)被截斷了，所以在沒有使用PME方法的狀況下，最好要將截斷值盡可能設(shè)大。需要提醒的是，模擬的耗時是與截斷值的平方成正比的，（參見 教程一 的 5.1.2節(jié)）所幸我們模擬的體系非常小，足夠承受沒有截斷值(cut=999)的計算?；谕瑯拥脑砦覀儗gbmax也設(shè)置為 999埃，這個參數(shù)控制了在計算非鍵相互作

20、用過程中列用于計算有效波恩半徑的粒子對的最遠間距。這個值設(shè)定的越大，計算的結(jié)果就越好，當然也就需要花費越多的計算時間?？紤]到我們面對的體系只有20個氨基酸殘基我們可以把所有的粒子都納入到有效波恩半徑中來，所以我們設(shè)定的rgbmax遠遠大于計算的尺度。.下面開始運行優(yōu)化過程: $AMBERHOME/exe/sander -O -i min1.in -o min1.out -p TC5b.prmtop -c TC5b.inpcrd -r min1.rstInput Files: TC5b.prmtop, TC5b.inpcrd, min1.inOutput Files: min1.out, min

21、1.rst在16個1.3GHzCPU的 SGI Altix上這個過程需要3.5秒完成為了直觀的比較優(yōu)化前后的結(jié)構(gòu)，我們生成一個pdb文件: $AMBERHOME/exe/ambpdb -p TC5b.prmtop min1.pdb將優(yōu)化前后的兩個文件打開(min1.pdb and TC5b_linear.pdb)你可以選擇任何可用的顯示軟件，比如VMD起始結(jié)構(gòu)用藍色顯示，優(yōu)化后的結(jié)構(gòu)用黃色顯示。如你所見，優(yōu)化過程并未造成主鏈結(jié)構(gòu)太大的變化但是色氨酸和酪氨酸殘基發(fā)生了比較明顯的移動，這些能量熱點集中的區(qū)域有可能在我們開始分子動力學模擬之后帶來麻煩，如果你不相信，可以用未經(jīng)優(yōu)化的結(jié)構(gòu)跑一個動力學過

22、程看看，肯定飛！第四步：體系加熱.接下來我們要在這個體系中正式開始分子動力學模擬，首先我們要分7步花費50ps時間對體系進行升溫模擬。將升溫過程分為7步完成可以在每一步升溫之后維持一段時間，以免一次升溫造成體系能量聚集最終跑飛，另外一種可行的方法是對升溫過程加一個權(quán)重限制。您可以參閱AMBER用戶手冊以獲取更多信息。一般而言我們升溫的最終目標是室溫即300K但是為了重復文獻的運算我們選擇325K：MD simulations of 100 ns were performed at 300 K, but all were kinetically trapped on this time scal

23、e, showing strong dependence on initial conditions and failing to converge to similar conformational ensembles. We therefore increased the temperature to 325 K.文獻認為必須將體系加溫到325K進行模擬，否則有可能使模擬的結(jié)果最終落入局部最小點，所以我們也做同樣的設(shè)定。但是你必須牢記更高的模擬溫度會導致體系中各化學鍵發(fā)生更加顯著的振動，這意味著如果你打算做一個600K，以2fs為步長的動力學模擬，你就要考慮一個應用了shaken的300k

24、效果會與之相同，但600K的模擬卻要臨步長過大的問題，過大的步長會導致體系不穩(wěn)定。還好325K不算太高，還比較接近常用的300K，2fs的步長可以處理含氫的鍵的振動?？墒羌偃缥覀円?00K的條件下運行動力學模擬的話，那模擬的步長就要縮減到1.5fs。我們的起始結(jié)構(gòu)是手工搭建的，不是通常常見的來自實驗的晶體結(jié)構(gòu)，所以我們的體系在模擬的開始階段要面臨不如晶體結(jié)構(gòu)穩(wěn)定的問題。為了讓我們的體系能夠在可控制的狀況下來釋放能量，在模擬起始的升溫階段我們選擇步長為0.5fs，進入相對穩(wěn)定的生成相之后，我們再選擇常規(guī)的2fs步長0.5fs的步長確實有些矯枉過正，但是保證體系的安全畢竟還是最重要的。我們進行升

25、溫模擬的方案如下： 第一階段 - 10,000 步, 步長 0.5fs (共5ps), 起始溫度 0.0K, 結(jié)束溫度 50.0K, 溫度耦合系數(shù) 1.0ps 第二階段 - 10,000 步, 步長 0.5fs (共5ps), 結(jié)束溫度 100.0K, 溫度耦合系數(shù) 1.0ps 第三階段 - 10,000 步, 步長 0.5fs (共5ps), 結(jié)束溫度 150.0K, 溫度耦合系數(shù) 1.0ps 第四階段 - 10,000 步, 步長 0.5fs (共5ps), 結(jié)束溫度 200.0K, 溫度耦合系數(shù) 1.0ps 第五階段 - 10,000 步, 步長 0.5fs (共5ps), 結(jié)束溫度 2

26、50.0K, 溫度耦合系數(shù) 1.0ps 第六階段 - 10,000 步, 步長 0.5fs (共5ps), 結(jié)束溫度 300.0K, 溫度耦合系數(shù) 1.0ps 第七階段 - 40,000 步, 步長 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

27、.0, ntpr=50, ntwx=50, ntb=0, igb=1, cut=999.,rgbmax=999. /其他六個階段的輸入文件與之非常接近，只需要改變一下相應的溫度就可以了，可以從此處下載現(xiàn)成的輸入文件： (heat2.in, heat3.in, heat4.in, heat5.in, heat6.in, heat7.in).下面是一個運行升溫模擬的PBS腳本，你也可以根據(jù)你的系統(tǒng)自己寫一個腳本。#PBS -l ncpus=16#PBS -l walltime=500:00:00#PBS -l cput=2000:00:00#PBS -j oesetenv AMBERHOME /u

28、sr/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.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.mdcrdgz

29、ip -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

30、-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/sa

31、nder -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 hea

32、t1.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 hea

33、t5.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.mdcrd gzip -9

34、 heat7.mdcrdmkdir initial_heatingcp heat1.out initial_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 he

35、at3.mdcrd.gz initial_heatingcp heat4.mdcrd.gz initial_heatingcp heat5.mdcrd.gz initial_heatingcp heat6.mdcrd.gz initial_heatingcp heat7.mdcrd.gz initial_heatingecho DONE在16個1.3GHz CPU的SGI Altix上，7個步驟全部完成共耗時7分鐘。下面提供了全部過程的輸出文件，你可以分別下載，也可以下載最終的一個壓縮文件。Heating StageOutput FileRestrt FileMdcrd FileStage 1

36、heat1.outheat1.rstheat1.mdcrd.gzStage 2heat2.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

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

38、明之前的計算獲得的就是一個穩(wěn)定的結(jié)果。 盡管文獻的作者發(fā)現(xiàn)在模擬的最初520ns中蛋白就已經(jīng)充分折疊，我們在本教程中仍然進行50ns的動力學計算，以重復文獻的報道。Two independent simulations converged to essentially identical families of structures after 5 and 20 ns.我們將這個總時間長度為50ns的模擬分為10個階段，每段5ns，這樣做是為了一旦系統(tǒng)崩潰，我們不會損失已經(jīng)進行的所有工作。而且這樣分開處理還可以保證每個輸出文件和軌跡文件的大小都適合處理。這10個階段的模擬我們會使用相同的輸入文

39、件，文件如下所示:equil.inStage 2 equilibration 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. /每一個階段的模擬都會進行250,000步(由nstlim的取值決定)，步長為2fs(由dt的取值決定) 即總共進行 5 ns的模擬。請注意我們在模擬全過程中使用了S

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

41、您系統(tǒng)的狀況修改腳本。#PBS -l ncpus=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 $A

42、MBERHOME/exe/sander -O -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/s

43、ander -O -i equil.in -p 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 eq

44、uil.in -p TC5b.prmtop -c 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

45、.prmtop -c equil7.rst -r 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 equ

46、il9.rst -r equil10.rst -o equil10.out -x equil10.mdcrdgzip -9 equil10.mdcrdecho DONE在16個1.3GHz CPU的SGI Altix上全部的十段模擬共耗時27小時。下面是每一步模擬的結(jié)果輸出，您可以下載這些結(jié)果，但是您需要注意每一個軌跡文件在壓縮之后都有13Mb 左右你可以將軌跡文件載入到VMD等看圖軟件中，提醒一下，你要首先解壓縮這個軌跡文件，然后載入它，整個過程有可能需要花費一天時間，然后你就可以欣賞整個動力學過程了，如果你用飄帶模型來顯示的話會非常好看。下面是一個抓圖：第六步：分析模擬結(jié)果下面我們要開始分析模擬的結(jié)果，我們要繪制模擬過程的溫度、總能量、動能、勢能變化曲線，這些信息可以從輸出文件中獲取，檢查這些數(shù)據(jù)也可以讓我們知道在動力學模擬過程中有沒有發(fā)生什么異常狀況。我們可以看到，溫度的曲線非常平滑，并且維持在325K。動能和勢能的曲線也一樣，會在平衡