江建敏　李娟　方慧生

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

關(guān)鍵詞：分子動(dòng)力學(xué)模擬；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ì)的空間結(jié)構(gòu)在行使其功能時(shí)起著重要的作用，其功能與動(dòng)力學(xué)特征密切相關(guān)，不同結(jié)構(gòu)決定了蛋白質(zhì)的不同功能。蛋白質(zhì)分子內(nèi)部運(yùn)動(dòng)和構(gòu)象變化對其正常功能有重要作用，這些過程包括蛋白質(zhì)的折疊和拉伸，配體和底物的結(jié)合等。要了解蛋白質(zhì)的特性并找到其發(fā)揮生物功能的分子基礎(chǔ)模型，只有氨基酸順序和三維結(jié)構(gòu)是不夠的，還必須研究蛋白質(zhì)的結(jié)構(gòu)和功能的關(guān)系等。

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

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

1 GROMACS簡介

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

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

2 模擬流程

分子動(dòng)力學(xué)模擬的步驟主要包括四步。

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

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

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

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

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

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

首先可以使用GROMACS的g_rms程序進(jìn)行計(jì)算蛋白質(zhì)各個(gè)部分的RMSD值隨模擬時(shí)間的變化，查看結(jié)構(gòu)的穩(wěn)定性。RMSD值可以反應(yīng)出模擬過程中體系的穩(wěn)定情況，當(dāng)體系的RMSD趨于穩(wěn)定時(shí)，說明動(dòng)力學(xué)模擬達(dá)到了平衡。同樣可以利用g_gyrate程序計(jì)算體系回旋半徑的變化，從而度量結(jié)構(gòu)的緊密度?；匦霃绞敲枋龅鞍踪|(zhì)結(jié)構(gòu)緊密性的一個(gè)物理量，半徑越小說明致密性越好，即蛋白質(zhì)結(jié)構(gòu)就越穩(wěn)定。利用g_energy程序獲得體系的勢能隨時(shí)間的變化情況。隨著模擬時(shí)間的不斷增加，體系也會(huì)趨于平衡，而勢能也是趨于平衡的。

3 模擬軟件的比較

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

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

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

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

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

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

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

NAMD是并行分子動(dòng)力學(xué)軟件，它被設(shè)計(jì)用于生物大分子的快速仿真，適用性極其廣泛?？梢允褂肁mber，CHARMM，GROMACS等多種力場，適合模擬蛋白質(zhì)、核酸等生物大分子。CHARMM力場適用范圍包括化學(xué)小分子到生物大分子，同時(shí)也是NAMD的默認(rèn)力場。

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

4 制約因素

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

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

5 前景展望

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

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

參考文獻(xiàn)：

[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.

編輯/成森