2016年10月5日 星期三

模擬Tryptophan cage的分子動態

Tryptophan cage是個人造的蛋白質,長度只有20個氨基酸,序列如下:
NLYIQ WLKDG GPSSG RPPPS

Tryptophan cage在一般的環境下會摺疊成固定的立體結構、稱之為Trp-cage motif。立體結構可以從蛋白質結構資料庫的 1L2Y 下載。由於此結構是用NMR解出來的,會有許多相似的結構。下面的Tryptophan cage的分子動態模擬(molecular dynamics simulation)用NMR的第一個結構當作是初始結構(檔案名稱叫做1L2Ym1.pdb),參數設定參考Justin A. Lemkul寫的《GROMACS Tutorial:Lysozyme in Water》,在Mac上面跑模擬。指令如下:
  1. sudo port install gromacs
  2. port info gromacs
  3. gmx pdb2gmx -f 1L2Ym1.pdb -o processed.gro -ff oplsaa -water tip4p
  4. gmx editconf -f processed.gro -o newbox.gro -c -d 1.0 -bt cubic
  5. gmx solvate -cp newbox.gro -cs tip4p.gro -o solv.gro -p topol.top
  6. gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
  7. gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -nname CL -nn 1    
    選13 SOL
  8. gmx grompp -f minim.mdp -c solv_ions.gro -p topol.top -o em.tpr
  9. gmx mdrun -nt 4 -v -deffnm em
  10. gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr
  11. gmx mdrun -nt 4 -v -deffnm nvt
  12. gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
  13. gmx mdrun -nt 4 -v -deffnm npt
  14. gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
  15. gmx mdrun -nt 4 -v -deffnm md
  16. gmx trjconv -s md.tpr -f md.xtc -o md_0-1ns_noPBC.xtc -pbc mol -ur compact
    選1 Protein
  17. gmx rms -s em.tpr -f md_0-1ns_noPBC.xtc -tu ns -o rmsd.xvg
    選4 Backbone、再選4 Backbone
  18. gmx gyrate -s md.tpr -f md_0-1ns_noPBC.xtc -o gyrate.xvg
    選1 Protein
  19. gmx trjconv -s md.tpr -f md_0-1ns_noPBC.xtc -o md_0-1ns.pdb
    選1 Proteins

上述指令的說明如下:
  1. 在Mac上用Macports安裝GROMACS
  2. 顯示安裝的GROMACS版本,目前安裝的版本是5.1.4
  3. 分子模擬用到的force field與water model,參考《分子模擬選擇的水分子模型》
    • 1L2Ym1.pdb 是初始結構。注意!此初始結構事先把所有的氫原子都去除了,讓force field根據此結構來加入氫原子
    • force field是 OPLS-AA/L
    • 水分子模型是 TIP4P,GROMACS推薦oplsaa用tip4p
    • 跑完這個指令後會會出現這個系統總共有20 residues with 304 atoms、total mass 2170.437 a.m.u.、total charge 1.000e。系統帶的電荷後面會說如何加入離子來平衡
  4. 把上述的分子放到一個盒子裡面,之後跑的分子模擬會用到
    • -c 代表把分子放在盒子的中央
    • -d 1.0 代表分子到盒子邊界的距離是1.0 nm
    • -bt cubic 代表這個盒子是一個cubic
  5. 在盒子裡面注入水分子
  6. 準備加入離子來平衡電性,這邊的ions.mdp從Justin Lemkul的網站用取得
    • wget http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/Files/ions.mdp
  7. 參考步驟3的結果,加入一個氯離子來平衡系統的電性
    • -nname CL 代表陰離子選用氯離子
    • -nn 1 代表要一個陰離子,如果加入的是一個陽離子則用 -np 1
    • 輸入指令後選13 SOL,意思是說把隨機拿掉一個水、代換成上面的離子
  8. 設定 energy minimization 需要的參數,minim.mdp從Justin Lemkul的網站用取得
    • wget http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/Files/minim.mdp
    • 這邊用的演算法是 steepest descent minimization
    • 跑到maximum force < 1000 kJ/mol/nm,或是跑50000 步
    • -nt 4 代表用4個核心來跑
    • -defnm em 代表跑出來的檔案都用 em.* 來儲存
  9. 跑 energy minimization
  10. 設定 equilibration under an NVT ensemble (constant Number of particles, Volume, and Temperature are all constant)。nvt.mdp設定檔從Justin Lemkul的網站用取得,總共跑 100 ps
    • wget http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/Files/nvt.mdp
    • 設定檔裡面 ref_t = 300 代表說溫度設定是300 K ≈ 27°C
    • 設定檔裡面 gen_temp = 300 代表根據300 K的Maxwell distribution用給予每個原子一個初始速度
  11. 跑 NVT equilibration
  12. 設定 equilibration under an NPT ensemble (Number of particles, Pressure, and Temperature are all constant)。npt.mdp設定檔從Justin Lemkul的網站用取得,總共跑 100 ps
    • wget http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/Files/npt.mdp
    • 設定檔裡面 continuation = yes 代表原子初始速度延續 NVT equilibration 最後的結果
    • 設定檔裡面 gen_vel = no 代表說不設定初始速度
    • 設定檔裡面的 ref_p = 1.0 代表壓力設定成 1.0 bar ≈ 1大氣壓
  13. 跑 NPT equilibration
  14. 設定分子模擬的參數。md.mdp設定檔從Justin Lemkul的網站用取得,設定成 300 K, 1.0 bar,約莫就是27 °C、1大氣壓的環境
  15. 跑Tryptophan cage的分子動態模擬,總共跑 1000 ps = 1 ns
    • wget http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/Files/md.mdp
    • 設定檔裡面 dt = 0.002 代表每一步的大小是0.002 ps = 2 fs
    • 設定檔裡面 constraint_algorithm = lincs 
    • 設定檔裡面 coulombtype = PME 代表long-range electrostatics用Particle Mesh Ewald這個演算法
    • 設定檔裡面 tcoupl = V-rescale 代表用modified Berendsen thermostat
    • 設定檔裡面 pcoupl = Parrinello-Rahman 代表用此演算法設定壓力
    • 跑完以後會出現Performance,這個例子跑了 1h00:56,23.627 ns/day、1.016 hour/ns
  16. 把模擬的軌跡檔(trajectory file)存出來,參考trjconv的說明
    • -pbc mol 意思是把分子的coler of mass放在盒子的正中央
    • -ur compact 代表將所有的原子放在離正中央最近的位置
  17. 用RMSD計算模擬過程中每個時間點對初始結構的差異
    • rmsd.xvg 這個檔案可以用xmgrace去看
  18. 計算模擬過程中radius of gyration的改變
    • gyrate.xvg 這個檔案可以用xmgrace去看
  19. 用PDB格式輸出軌跡檔,可以直接用PyMOL觀看
如果覺得1 ns模擬太短,想要繼續跑下去,請參考《GROMACS續跑》的說明

_EOF_

沒有留言:

張貼留言