GROMACS分子動力學教程：水中溶菌酶的分子動力學模擬

溶菌酶 （Lysozyme）

又稱

胞壁質酶 （Muramidase）

，存在於卵清、唾液等生物分泌液中，能夠催化細菌細胞壁肽聚糖N-乙醯氨基葡糖和 N-乙醯胞壁酸之間的β-1，4 糖苷鍵的水解。當與帶負電荷的病毒蛋白相互結合時，溶菌酶能夠與 DNA、RNA、脫輔基蛋白形成複鹽，這種複鹽可使病毒失去活力。蛋清溶菌酶是目前研究最為透徹的一種溶菌酶。本例中就以

雞蛋清溶菌酶

為例，為大家演示使用

Gromacs

完成其在水中的動力學模擬。考慮到許多同學不熟悉Linux系統的操作，本例所有操作均可在Windows系統下完成，關於Windows系統下Groamcs的安裝，大家可參考下面這篇博文：

http://sobereva.com/458

01結構處理

首先在PDB資料庫中下載溶菌酶蛋白質的晶體結構，PDB ID為1AKI；

下載好結構後，使用PyMOL檢查結構，刪除雜原子和水分子。對於水的處理，刪除所有的水分子並非普適規則，特別是在某些情況下，緊密結合的水可能具有某種功能意義，這時必須進行適當的建模。如在研究蛋白質和配體之間的親和作用時，結合在活性位點附近的比如當水參與了小分子結合作用時就不能刪除。

下圖為處理完成的結構，將其另存為“

1AKI_clean.pdb

”檔案。

02產生拓撲檔案

下面我們需要使用

pdb2gmx

模組產生拓撲檔案。在工作目錄中開啟命令視窗，輸入以下命令：

gmx pdb2gmx -f 1AKI_clean。pdb -o 1AKI_processed。gro -water spce -ignh

其中“

-f

”指定座標檔案；“

-o

”輸出檔案；“

-water

”指定使用的水模型，本例中使用SPC/E模型的水；“

-ignh

”設定忽略PDB檔案中的氫原子。

此時會出現選擇力場選擇的提示，如下：

1： AMBER03 protein， nucleic AMBER94 （Duan et al。， J。 Comp。 Chem。 24， 1999-2012， 2003）2： AMBER94 force field （Cornell et al。， JACS 117， 5179-5197， 1995）3： AMBER96 protein， nucleic AMBER94 （Kollman et al。， Acc。 Chem。 Res。 29， 461-469， 1996）4： AMBER99 protein， nucleic AMBER94 （Wang et al。， J。 Comp。 Chem。 21， 1049-1074， 2000）5： AMBER99SB protein， nucleic AMBER94 （Hornak et al。， Proteins 65， 712-725， 2006）6： AMBER99SB-ILDN protein， nucleic AMBER94 （Lindorff-Larsen et al。， Proteins 78， 1950-58， 2010）7： AMBERGS force field （Garcia & Sanbonmatsu， PNAS 99， 2782-2787， 2002）8： CHARMM27 all-atom force field （CHARM22 plus CMAP for proteins）9： GROMOS96 43a1 force field10： GROMOS96 43a2 force field （improved alkane dihedrals）11： GROMOS96 45a3 force field （Schuler JCC 2001 22 1205）12： GROMOS96 53a5 force field （JCC 2004 vol 25 pag 1656）13： GROMOS96 53a6 force field （JCC 2004 vol 25 pag 1656）14： GROMOS96 54a7 force field （Eur。 Biophys。 J。 （2011）， 40，， 843-856， DOI： 10。1007/s00249-011-0700-9）15： OPLS-AA/L all-atom force field （2001 aminoacid dihedrals）

力場的選擇至關重要，如果對於力場的瞭解不夠，可選擇與相關文獻一致的力場。本例中選擇全原子OPLS-AA力場，輸入15，回車即可。

注意，此時工作目錄產生三個新的檔案：

1AKI_processed.gro

、

topol.top

以及

posre.itp

；

1AKI_processed.gro

是

Gromacs

的結構檔案格式，包含力場中定義的所有原子，

topol.top

是體系的拓撲檔案，

posre.itp

中包含了用於限制重原子位置的資訊。

03定義盒子並溶劑化

我們模擬一個簡單的蛋白水溶液體系，可在適當的情況下模擬其他不同溶劑中的蛋白質或其他分子，只要有合適的力場引數即可。

首先使用

editconf

模組定義盒子的大小，輸入以下命令：

gmx editconf -f 1AKI_processed。gro -o 1AKI_newbox。gro -c -d 1。0 -bt cubic

上述命令會將蛋白置於立方體盒子的中心

（-bt cubic）

，此外，“

-d

”指定分子到盒子邊緣的距離，“1。0”表示1nm。

下面使用solvate模組給盒子中填充溶劑（水），輸入以下命令：

gmx solvate -cp 1AKI_newbox。gro -cs spc216。gro -o 1AKI_solv。gro -p topol。top

其中“

spc216.gro

”為溶劑構型檔案，為Gromacs自帶的。如果我們使用SPC、TIP3P等三點水模型，均可使用“

spc216.gro

”作為溶劑構型。

下面我們使用視覺化軟體VMD檢視一下溶劑化的體系。

04新增抗衡離子

溶劑化的體系包含了一個帶電荷的蛋白，如果大家注意了

pdb2gmx

模組的輸出結果，可以看到蛋白質帶有+8e的淨電荷。分子動力學模擬通常是在電中性的條件下進行的，因此我們需要在體系中新增抗衡離子，來使體系處於中性電荷的狀態。

新增抗衡離子需要兩步，第一步需要一個“

ions.mdp

”的引數檔案，本例中使用的“

ions.mdp

”檔案內容如下：

； ions。mdp - used as input into grompp to generate ions。tpr； Parameters describing what to do， when to stop and what to saveintegrator = steep ； Algorithm （steep = steepest descent minimization）emtol = 1000。0 ； Stop minimization when the maximum force < 1000。0 kJ/mol/nmemstep = 0。01 ； Minimization step sizensteps = 50000 ； Maximum number of （minimization） steps to perform； Parameters describing how to find the neighbors of each atom and how to calculate the interactionsnstlist = 1 ； Frequency to update the neighbor list and long range forcescutoff-scheme = Verlet ； Buffered neighbor searchingns_type = grid ； Method to determine neighbor list （simple， grid）coulombtype = cutoff ； Treatment of long range electrostatic interactionsrcoulomb = 1。0 ； Short-range electrostatic cut-offrvdw = 1。0 ； Short-range Van der Waals cut-offpbc = xyz ； Periodic Boundary Conditions in all 3 dimensions

下面輸入以下命令來使用

grompp

模組生成

.tpr

檔案：

gmx grompp -f ions。mdp -c 1AKI_solv。gro -p topol。top -o ions。tpr

這時產生了

.tpr

結尾的檔案，該檔案提供了體系的原子級別描述，還包含了模擬引數說明。

第二步就是使用

genion

模組新增抗衡離子，將一些溶劑水分子替換為離子，輸入以下命令：

gmx genion -s ions。tpr -o 1AKI_solv_ions。gro -p topol。top -pname NA -nname CL -neutral

出現提示時，選擇溶劑組“

13（SOL）

”。

05能量最小化

在開始動力學模擬之前，我們需要保證體系的結構正常，幾何構型合理，原子沒有衝撞，此時就需要先將溶劑化體系弛豫至能量較低的狀態，即能量最小化。與上一步類似，我們也需要一個“

minim.mdp

”檔案，本例中使用的“

minim.mdp

”檔案內容如下：

； minim。mdp - used as input into grompp to generate em。tpr； Parameters describing what to do， when to stop and what to saveintegrator = steep ； Algorithm （steep = steepest descent minimization）emtol = 1000。0 ； Stop minimization when the maximum force < 1000。0 kJ/mol/nmemstep = 0。01 ； Minimization step sizensteps = 50000 ； Maximum number of （minimization） steps to perform； Parameters describing how to find the neighbors of each atom and how to calculate the interactionsnstlist = 1 ； Frequency to update the neighbor list and long range forcescutoff-scheme = Verlet ； Buffered neighbor searchingns_type = grid ； Method to determine neighbor list （simple， grid）coulombtype = PME ； Treatment of long range electrostatic interactionsrcoulomb = 1。0 ； Short-range electrostatic cut-offrvdw = 1。0 ； Short-range Van der Waals cut-offpbc = xyz ； Periodic Boundary Conditions in all 3 dimensions

下面我們仍然是使用

grompp

模組將結構、拓撲和模擬引數整合到

.tpr

檔案中，輸入以下命令：

gmx grompp -f minim。mdp -c 1AKI_solv_ions。gro -p topol。top -o em。tpr

下面使用

mdrun

模組來執行能量最小化過程，輸入以下命令：

gmx mdrun -v -deffnm em

06預平衡

在開始真正的動力學模擬之前，我們需要對蛋白質周圍的溶劑和離子進行預平衡。預平衡通常分兩個階段進行，第一個階段在NVT系綜（也稱正則系綜，等溫等容系綜）下進行，其中體系的粒子數（N），體積（V） 和平均溫度（T） 保持不變。本例中進行100 ps的NVT預平衡，需要一個“

nvt.mdp

”檔案，本例中“

nvt.mdp

”檔案的內容如下：

title = OPLS Lysozyme NVT equilibrationdefine = -DPOSRES ； position restrain the protein； Run parametersintegrator = md ； leap-frog integratornsteps = 50000 ； 2 * 50000 = 100 psdt = 0。002 ； 2 fs； Output controlnstxout = 500 ； save coordinates every 1。0 psnstvout = 500 ； save velocities every 1。0 psnstenergy = 500 ； save energies every 1。0 psnstlog = 500 ； update log file every 1。0 ps； Bond parameterscontinuation = no ； first dynamics runconstraint_algorithm = lincs ； holonomic constraintsconstraints = h-bonds ； bonds involving H are constrainedlincs_iter = 1 ； accuracy of LINCSlincs_order = 4 ； also related to accuracy； Nonbonded settingscutoff-scheme = Verlet ； Buffered neighbor searchingns_type = grid ； search neighboring grid cellsnstlist = 10 ； 20 fs， largely irrelevant with Verletrcoulomb = 1。0 ； short-range electrostatic cutoff （in nm）rvdw = 1。0 ； short-range van der Waals cutoff （in nm）DispCorr = EnerPres ； account for cut-off vdW scheme； Electrostaticscoulombtype = PME ； Particle Mesh Ewald for long-range electrostaticspme_order = 4 ； cubic interpolationfourierspacing = 0。16 ； grid spacing for FFT； Temperature coupling is ontcoupl = V-rescale ； modified Berendsen thermostattc-grps = Protein Non-Protein ； two coupling groups - more accuratetau_t = 0。1 0。1 ； time constant， in psref_t = 300 300 ； reference temperature， one for each group， in K； Pressure coupling is offpcoupl = no ； no pressure coupling in NVT； Periodic boundary conditionspbc = xyz ； 3-D PBC； Velocity generationgen_vel = yes ； assign velocities from Maxwell distributiongen_temp = 300 ； temperature for Maxwell distributiongen_seed = -1 ； generate a random seed

下面輸入以下命令執行NVT平衡：

gmx grompp -f nvt。mdp -c em。gro -r em。gro -p topol。top -o nvt。tprgmx mdrun -deffnm nvt

上面階段的NVT預平衡穩定了體系的溫度，下面我們對體系進行第二階段的預平衡，在NPT系綜（也稱等溫等壓系綜）下進行，其中粒子數（N），平均壓力（P） 和平均溫度（T） 保持不變，穩定體系的壓力，進而穩定密度。本例中進行100 ps的NPT預平衡，需要一個“

npt.mdp

”檔案，本例中“

npt.mdp

”檔案的內容如下：

title = OPLS Lysozyme NPT equilibrationdefine = -DPOSRES ； position restrain the protein； Run parametersintegrator = md ； leap-frog integratornsteps = 50000 ； 2 * 50000 = 100 psdt = 0。002 ； 2 fs； Output controlnstxout = 500 ； save coordinates every 1。0 psnstvout = 500 ； save velocities every 1。0 psnstenergy = 500 ； save energies every 1。0 psnstlog = 500 ； update log file every 1。0 ps； Bond parameterscontinuation = yes ； Restarting after NVTconstraint_algorithm = lincs ； holonomic constraintsconstraints = h-bonds ； bonds involving H are constrainedlincs_iter = 1 ； accuracy of LINCSlincs_order = 4 ； also related to accuracy； Nonbonded settingscutoff-scheme = Verlet ； Buffered neighbor searchingns_type = grid ； search neighboring grid cellsnstlist = 10 ； 20 fs， largely irrelevant with Verlet schemercoulomb = 1。0 ； short-range electrostatic cutoff （in nm）rvdw = 1。0 ； short-range van der Waals cutoff （in nm）DispCorr = EnerPres ； account for cut-off vdW scheme； Electrostaticscoulombtype = PME ； Particle Mesh Ewald for long-range electrostaticspme_order = 4 ； cubic interpolationfourierspacing = 0。16 ； grid spacing for FFT； Temperature coupling is ontcoupl = V-rescale ； modified Berendsen thermostattc-grps = Protein Non-Protein ； two coupling groups - more accuratetau_t = 0。1 0。1 ； time constant， in psref_t = 300 300 ； reference temperature， one for each group， in K； Pressure coupling is onpcoupl = Parrinello-Rahman ； Pressure coupling on in NPTpcoupltype = isotropic ； uniform scaling of box vectorstau_p = 2。0 ； time constant， in psref_p = 1。0 ； reference pressure， in barcompressibility = 4。5e-5 ； isothermal compressibility of water， bar^-1refcoord_scaling = com； Periodic boundary conditionspbc = xyz ； 3-D PBC； Velocity generationgen_vel = no ； Velocity generation is off

下面輸入以下命令執行NPT平衡：

gmx grompp -f npt。mdp -c nvt。gro -r nvt。gro -t nvt。cpt -p topol。top -o npt。tprgmx mdrun -deffnm npt

07分子動力學模擬

預平衡階段完成後，體系處於適當的溫度和壓力下，現在可以放開蛋白質的位置限制，進行無限制分子動力學模擬了。這一過程預前面的步驟類似，需要的引數檔案“

md.mdp

”內容如下：

title = OPLS Lysozyme NPT equilibration； Run parametersintegrator = md ； leap-frog integratornsteps = 500000 ； 2 * 500000 = 1000 ps （1 ns）dt = 0。002 ； 2 fs； Output controlnstxout = 0 ； suppress bulky 。trr file by specifyingnstvout = 0 ； 0 for output frequency of nstxout，nstfout = 0 ； nstvout， and nstfoutnstenergy = 5000 ； save energies every 10。0 psnstlog = 5000 ； update log file every 10。0 psnstxout-compressed = 5000 ； save compressed coordinates every 10。0 pscompressed-x-grps = System ； save the whole system； Bond parameterscontinuation = yes ； Restarting after NPTconstraint_algorithm = lincs ； holonomic constraintsconstraints = h-bonds ； bonds involving H are constrainedlincs_iter = 1 ； accuracy of LINCSlincs_order = 4 ； also related to accuracy； Neighborsearchingcutoff-scheme = Verlet ； Buffered neighbor searchingns_type = grid ； search neighboring grid cellsnstlist = 10 ； 20 fs， largely irrelevant with Verlet schemercoulomb = 1。0 ； short-range electrostatic cutoff （in nm）rvdw = 1。0 ； short-range van der Waals cutoff （in nm）； Electrostaticscoulombtype = PME ； Particle Mesh Ewald for long-range electrostaticspme_order = 4 ； cubic interpolationfourierspacing = 0。16 ； grid spacing for FFT； Temperature coupling is ontcoupl = V-rescale ； modified Berendsen thermostattc-grps = Protein Non-Protein ； two coupling groups - more accuratetau_t = 0。1 0。1 ； time constant， in psref_t = 300 300 ； reference temperature， one for each group， in K； Pressure coupling is onpcoupl = Parrinello-Rahman ； Pressure coupling on in NPTpcoupltype = isotropic ； uniform scaling of box vectorstau_p = 2。0 ； time constant， in psref_p = 1。0 ； reference pressure， in barcompressibility = 4。5e-5 ； isothermal compressibility of water， bar^-1； Periodic boundary conditionspbc = xyz ； 3-D PBC； Dispersion correctionDispCorr = EnerPres ； account for cut-off vdW scheme； Velocity generationgen_vel = no ； Velocity generation is off

本例中執行1ns（500000 步×0。002 ps/步= 1000 ps = 1 ns）的分子動力學模擬，輸入以下命令執行計算：

gmx grompp -f md。mdp -c npt。gro -t npt。cpt -p topol。top -o md_0_1。tprgmx mdrun -deffnm md_0_1 -nb gpu

下面我們使用視覺化軟體看一下動力學軌跡動畫。

參考資料：1。http：//www。mdtutorials。com/gmx/lysozyme/index。html2。https：//jerkwin。github。io/3。 http：//sobereva。com/458

版權資訊

本文系AIDD Pro接受的外部投稿，文中所述觀點僅代表作者本人觀點，不代表AIDD Pro平臺，如您發現釋出內容有任何版權侵擾或者其他資訊錯誤解讀，請及時聯絡AIDD Pro （請新增微訊號plgrace）進行刪改處理。