基于GROMACS的分子動力學模擬的研究

江建敏　李娟　方慧生

摘要：分子動力學模擬可以通過多款軟件進行，根據(jù)需要不同可以選擇不同的軟件，但GROMACS是應用比較廣泛的。生物大分子的結構動力學決定其生物學功能，而分子動力學模擬是研究其結構動力學的有效手段，因此采用GROMACS進行分子動力學計算分析從而指導進一步的實驗工作。本文重點介紹了GROMACS的特點、基本流程以及制約因素等，并展望了未來的應用前景。

關鍵詞：分子動力學模擬；GROMACS；力場

Abstract：Molecular dynamics simulation can be carried out through a variety of software， depending on the need to choose a different software， but GROMACS is widely used. biological functions of biological macromolecules is determined by its structural dynamics， and molecular dynamics simulation is an effective means to study its structural dynamics. Therefore，GROMACS is used to carry out molecular dynamics analysis to guide the further experimental work. This review mainly introduces the characteristics， basic process and the control of GROMACS， and the prospect of future prospects.

Key words：Molecular dynamics simulation；GROMACS；Field

蛋白質(zhì)的空間結構在行使其功能時起著重要的作用，其功能與動力學特征密切相關，不同結構決定了蛋白質(zhì)的不同功能。蛋白質(zhì)分子內(nèi)部運動和構象變化對其正常功能有重要作用，這些過程包括蛋白質(zhì)的折疊和拉伸，配體和底物的結合等。要了解蛋白質(zhì)的特性并找到其發(fā)揮生物功能的分子基礎模型，只有氨基酸順序和三維結構是不夠的，還必須研究蛋白質(zhì)的結構和功能的關系等。

目前，通過測定蛋白質(zhì)空間結構的方法可以觀察到當?shù)鞍滋幱诓煌臉嬒髸r，其對應的功能狀態(tài)也不同，但在這些不同的構象變化的細節(jié)仍不清楚。因此利用計算機研究蛋白質(zhì)的動態(tài)結構和動力學信息，在一定程度上彌補了實際實驗上的不足。在過去的幾十年中，隨著計算機軟硬件技術的快速提升，在生物大分子等方面的模擬逐漸發(fā)展為在原子水平上研究蛋白質(zhì)、脂類等生物大分子的結構和動力學性質(zhì)的一個很重要工具，從而輔助實驗進行一定的解析。計算機對生物大分子的模擬研究使人們進一步了解生命現(xiàn)象的基本過程，對生物大分子結構和功能的研究更加深入和系統(tǒng)。利用分子動力學模擬計算生物大分子的空間結構，從而得到系統(tǒng)內(nèi)部的構象隨著時間變化的各種信息，有助于我們進一步理解蛋白質(zhì)、核酸等生物大分子的運動特征。

1957年，Alder和Wainwright[1]最早利用分析動力學及硬球模型實現(xiàn)了第一次模擬。1977年，McCammon[2]等人發(fā)表了第一篇蛋白分子模擬研究論文。雖然模擬時間很短，但是使人們看到分子動力學技術對研究蛋白質(zhì)分子功能的重要作用。1980年，Andersen[3]等人創(chuàng)造了一個恒定的壓力分子動力學方法。分子動力學模擬在誕生至今的幾十年中不斷發(fā)展完善，已經(jīng)成為應用最廣泛的分子模擬方法之一。隨著分子動力學模擬的快速發(fā)展，也逐漸形成了一些比較完善的應用于生物領域的模擬軟件，例如NAMD[4]、GROMACS[5]等。GROMACS是一個非常有效并且免費的的軟件，遵循GNU許可，可以從網(wǎng)站http：//www.gromacs.org/上進行下載得到。本文主要介紹GROMACS的基本過程。

1 GROMACS簡介

GROMACS[6]，一款開源的分子動力學仿真軟件包，最初是由荷蘭格羅寧根大學（Groningen University）開發(fā)。目前主要由瑞典烏普薩拉大學、斯德哥爾摩大學和德國的馬普學會高分子研究所維護和管理。它開創(chuàng)了分子動力學模擬的先河，成為迄今為止應用很廣泛的動力學模擬軟件之一。

GROMACS設計的初衷主要是應用于生物大分子的分析計算。其模擬程序包除了包含Gromacs力場，還包含了其他動力學軟件的力場如CHARMM力場[7]、AMBER力場[8，9]，這些力場的加入也讓其應用性更加廣泛。GROMACS軟件主要是用來研究分子體系的動力學，還可用于分子動力學，隨機動力學等方法來模擬自然的分子動力學，分子能量優(yōu)化，構象分析，等等。相比于其他的軟件GROMACS提供了極高的性能，而且它允許使用標準的MPI通訊進行并行運算。同時GROMACS支持很多常用的算法，與其他動力學模擬軟件相比還具有獨特的優(yōu)勢，使其脫穎而出。

2 模擬流程

分子動力學模擬的步驟主要包括四步。

2.1初始結構處理 在分子動力學模擬中，分子結構文件需要轉(zhuǎn)化成分子的拓撲文件以供后續(xù)模擬使用。運行pdb2gmx程序后會要求我們選擇合適的分子力場。運行pdb2gmx程序這一步會將原始的pdb文件轉(zhuǎn)化成結構文件以及相應的拓撲文件。同時考慮到體系在真實環(huán)境中的狀態(tài)，首先需要采用editconf程序給體系添加一個盒子，然后采用genbox程序選擇合適的溶劑化模型。溶劑化模型一般分為顯性溶劑化模型和隱形溶劑化模型。同時為了使體系在分子動力學模擬過程保持中性，我們還需要采用genion程序添加抗衡離子。

然而在準備模擬的時候，體系中有很多原子距離太近，或者新加入的粒子位置不一定合理，導致局部能量過高。如果能量過高可能導致運動速度過快，從而使體系不穩(wěn)定，最終導致模擬失敗。因此在開始模擬之前，首先需要進行能量最優(yōu)化，使系統(tǒng)達到穩(wěn)定的狀態(tài)。能量優(yōu)化的過程是改變系統(tǒng)中局部高能量原子的位置，降低這些部位的能量。

2.2體系優(yōu)化 進行能量最優(yōu)化的程序是mdrun，它是GROMACS的核心程序。能量最優(yōu)化會消除高能量原子之間不合理的結構關系，使能量降低到局部極小值。所有的模擬都是采用 mdrun程序進行的，類似于amber軟件中的sander程序，就好像GROMACS的心臟一樣。

進行能量最小化后，就可以在控制壓強和溫度的前提下對體系進行平衡模擬，使體系能量達到全局最小值時，整個體系達到平衡態(tài)。首先在NVT系統(tǒng)中，該系統(tǒng)稱為等溫等容系。然后經(jīng)過NVT平衡，溫定了系統(tǒng)的溫度，還需要穩(wěn)定系統(tǒng)的壓強。壓強平衡在NPT系統(tǒng)中進行，該系統(tǒng)被稱為等溫等容系。經(jīng)過兩個階段的完成，系統(tǒng)已經(jīng)在需要的溫度和壓強下平衡良好，可以進行正式模擬了。

2.3正式模擬 經(jīng)過上一步的平衡相過程后，分子系統(tǒng)已經(jīng)可以進行生產(chǎn)相過程即完全無限制的分子動力學模擬。然后再次調(diào)用grompp程序?qū)⑤斎胛募D(zhuǎn)化成tpr預處理文件，進行真正的分子動力學模擬，從而得到按時間排列的所有原子的運動情況。為了在一定的計算資源和時間下，盡可能增長實際模擬時間。

2.4分析結果 模擬完成后就可以對體系進行分析，從模擬的軌跡中提取有用的信息，并加以分析解釋才是最重要的。而作為開源的分子動力學軟件的GROMACS，另一個主要優(yōu)勢是提供了一系列分析軌跡文件的小程序。

首先可以使用GROMACS的g_rms程序進行計算蛋白質(zhì)各個部分的RMSD值隨模擬時間的變化，查看結構的穩(wěn)定性。RMSD值可以反應出模擬過程中體系的穩(wěn)定情況，當體系的RMSD趨于穩(wěn)定時，說明動力學模擬達到了平衡。同樣可以利用g_gyrate程序計算體系回旋半徑的變化，從而度量結構的緊密度。回旋半徑是描述蛋白質(zhì)結構緊密性的一個物理量，半徑越小說明致密性越好，即蛋白質(zhì)結構就越穩(wěn)定。利用g_energy程序獲得體系的勢能隨時間的變化情況。隨著模擬時間的不斷增加，體系也會趨于平衡，而勢能也是趨于平衡的。

3 模擬軟件的比較

目前，用于分子動力學模擬的軟件越來越成熟。廣泛應用于生物大分子研究的除了Gromacs之外，常用的還有NAMD、CHARMM、AMBER、TINKER、LAMMPS等。

3.1應用范圍 分子動力學模擬[11]作為一種理想的計算機模擬方法，可以用來模擬兩相之間的相互作用，也是對計算和實驗的補充。目前，隨著計算機的快速發(fā)展，分子動力學模擬在化學工程、生命科學和醫(yī)學等領域的應用也比較成熟。

GROMACS是通用的分子動力學模擬軟件，尤其適用于模擬生物大分子體系。同時由于其速度快，在非生物體系統(tǒng)中也得到了廣泛的應用。其中GROMACS軟件的一個關鍵優(yōu)勢是用戶界面友好。而且GROMACS模擬中的參數(shù)條件和基本功能也已經(jīng)趨于成熟。此外，我們可以在本地系統(tǒng)上運行GROMACS而不用占用并行系統(tǒng)。

NAMD[15]是一套分子動力學模擬軟件，由美國伊利諾伊大學理論與計算生物物理小組開發(fā)。它包括Amber、CHARMM、OPLS、GROMACS等多種力場，適用于模擬計算蛋白質(zhì)、核酸等生物大分子體系，而且并行計算效率非常高。為了配合NAMD進行模擬結果分析，還需要VMD[6]作為分子可視化和輔助分析軟件。

AMBER[13]是由美國加利福尼亞大學舊金山分校的Peter A.Kollman小組開發(fā)的軟件包?，F(xiàn)在由羅格斯大學、猶他大學、美國國立環(huán)境衛(wèi)生科學研究所等機構進行維護管理。但是它只支持AMBER力場，其它力場都不適用。另外它的力場是免費的，而其分子模擬軟件需要付費使用，能夠模擬蛋白質(zhì)、核酸等生物大分子。

3.2力場 分子力場，有時被稱為勢函數(shù)，對于一個特定的目的，力場可以分為不同的形式，但也有不同的適用范圍和局限性。在分子力學和分子動力學的計算中，適當?shù)牧鲞x取是非常重要的，這是決定計算結構的重要條件。

GROMACS[12]是一款功能強大的分子動力學模擬軟件，它包括一系列的GROMOS力場、AMBER力場和CHARMM力場等，能普遍應用于處理生物大分子。GROMOS力場[14]誕生后不斷被改進，現(xiàn)階段廣泛使用的版本有GROMOS96 43a1、GROMOS96 43a2、GROMOS96 45a3等。GROMACS進行分子動力學模擬時首先采用pdb2gmx程序，然后會出現(xiàn)界面，要求我們選擇合適的分子力場，關于力場的選擇可以根據(jù)不同需要進行不同的選擇。

NAMD是并行分子動力學軟件，它被設計用于生物大分子的快速仿真，適用性極其廣泛。可以使用Amber，CHARMM，GROMACS等多種力場，適合模擬蛋白質(zhì)、核酸等生物大分子。CHARMM力場適用范圍包括化學小分子到生物大分子，同時也是NAMD的默認力場。

AMBER包括一系列用于分子模擬的分子力學力場和分子模擬的軟件包[13]。其分子力場在計算生物學領域得到了廣泛的應用。它是一種分子動力學模擬和分析程序，可以模擬蛋白質(zhì)、糖類、核酸等生物大分子。它的勢能函數(shù)很簡單，所需的參數(shù)不多，并且參數(shù)是從計算結果和實驗值的比對中得到。

4 制約因素

4.1兼容性 在進行分子動力學模擬時首先需要處理輸入的pdb文件，而GROMACS只能識別pdb初始文件，其他格式的文件均不能被其轉(zhuǎn)化從而進行模擬。而且如果初始文件含有缺失的原子或者雜原子，還需要采用其他軟件進行填補或者刪除。而且NAMD具有一款很好的分析輔助軟件VMD，并能夠兼容多種文件格式。

4.2結合自由能計算 在分子動力學模擬領域，準確計算結合自由能雖取得一定進展，但仍然是一個挑戰(zhàn)?，F(xiàn)在開發(fā)了多種計算結合自由能[17]的方法，可分為三類：第一類為自由能微擾，第二類為熱力學積分和使用最廣泛的第三類MM-PBSA/GBSA。MM-PB/GBSA是 Srinivasan 等人提出的一種自由能計算方法，用來計算復合物體系的絕對和相對結合自由能。此外，Amber軟件包中專門提供了MM-PBSA/GBSA程序，可以非常方便的計算出受體與配體間的結合自由能。而GROMACS中并沒有像AMBER中那樣單獨的模塊，所以計算比較麻煩。

5 前景展望

隨著計算機技術的迅速發(fā)展，人們對于生物大分子結構和功能的逐步認識，為分子動力學的進一步發(fā)展和完善提供了良好的氛圍。本文簡要介紹了分子動力學模擬的一款軟件GROMACS的基本流程和應用。傳統(tǒng)實驗需要大量的人力物力而且耗時，所以采用計算機模擬的方法能給予一定的指導，省時省力，提高效率。

目前，分子動力學模擬能在很大程度上進行預測指導實驗，因此模擬與實驗的對照是將來研究的主要方向之一。對模擬得到的結構進行分析的方法雖然已經(jīng)提出了多種，然而不能認為它們已經(jīng)構成了結構分析方法的完整集合。分子動力學模擬不僅廣泛應用于生物大分子中，而且在藥物設計和開發(fā)中具有廣闊的應用前景。分子動力學模擬研究藥物中配體受體間相互作用的結構信息，為篩選藥物和藥物設計提供了重要線索。藥物的開發(fā)需要大量的資金，需要投入大量的人力和物力，而分子動力學模擬可以減少實驗的盲目性，從而節(jié)省人力物力，提高效率。

參考文獻：

[1] Alder B J，Wainwright T E. Phase transition for a hard sphere system[J].The Journal of chemical physics， 1957， 27（5）： 1208.

[2] McCammon J A， Karplus M. Internal motions of antibody molecules [J]. Nature， 1977， 268（5622）： 765-766.

[3]ANDENSONHC.Molecular dynamics simulations at constant press and/or temperature[J].Journal of Chemical Physics，1980， 72：2384-2393.

[4] Phillips J C， Braun R， Wang W， et al. Scalable molecular dynamics with NAMD[J]. Journal of computational chemistry，2005， 26（16）： 1781-1802.

[5]Hess B，Kutzner C， Van Der Spoel D， et al.GROMACS 4： algorithms for highly efficient， load-balanced， and scalable molecular simulation[J]. Journal of chemical theory and computation， 2008，4（3）：435-447.

[6] H.J.C. Berendsen， l. van der Spoe， R. van Drunen， Gromacs： A Message-Passing Parallel Molecular-Dynamics Implementation[J].Computer Physics Communications，1995， 91（1-3）： 43-56.

[7] Momany F.A.， Rone R. Validation of the general purpose QUANTA 3.2/CHARMm force field[J]. Journal of computational chemistry， 1992， 13（7）： 888-900.

[8] Weiner S.J.， Kollman P.A.， Case D.A.， et al. A new force field for molecular mechanical simulation of nucleic acids and proteins[J]. Journal of the American Chemical Society， 1984， 106（3）： 765-784.

[9] Weiner S.J.， Kollman P.A.， Nguyen D.T.， et al. An all atom force field for simulations of proteins and nucleic acids[J]. Journal of computational chemistry， 1986， 7（2）： 230-252.

[10] MACKAY D H J， CROSS A J， HAGLER A T. Prediction of protein structure and the principles of protein conformation [M]. New York： Plenum Press， 1989.

[11] Alder B.J.， Wainwright T.Studies in molecular dynamics. I. General method[J]. The Journal of Chemical Physics， 2004， 31（2）： 459-466.

[12] E. Lindahl， B. Hess， D. van der Spoel， Gromacs 3.0： A Package for Molecular Simulation and Trajectory Analysis[J].Journal of Molecular Modeling，2001， 7（8）： 306-317.

[13] D. A. Case， T. E. C. III， T. Darden， H. Gohlke， et al. Woods， The Amber Biomolecular Simulation Programs.[J].Journal of Computational Chemistry，2005， 26（16）： 1668-1688.

[14] Oostenbrink C， Villa A， Mark A E， et al. A biomolecular force field based on the free enthalpy of hydration and solvation： The GROMOS force‐field parameter sets 53A5 and 53A6[J]. Journal of computational chemistry， 2004， 25（13）： 1656-1676.

[15] L. Kale， R. Skeel， M. Bhandarkar， R. Brunner， A. Gursoy， N. Krawetz， J. Phillips， A. Shinozaki， K. Varadarajan， K. Schulten， Namd2： Greater Scalability for Parallel Molecular Dynamics[J].J Comput Phys， 1999， 151： 283-312.

[16] W.Humphrey， A. Dalke， K. Schulten， Vmd： Visual Molecular Dynamics[J]. J Mol Graph，1996， 14（1）： 33-38.

[17] Chipot C.， Pohorille A. Free energy calculations： theory and applications in chemistry and biology[M]. Springer， 2007.

編輯/成森