靜電作用力算法對(duì)分子動(dòng)力學(xué)模擬結(jié)果的影響分析

李 浩，白 姝

（天津大學(xué)化工學(xué)院，天津 300072；天津大學(xué)系統(tǒng)生物工程教育部重點(diǎn)實(shí)驗(yàn)室，天津 300072）

分子模擬是基于原子的結(jié)構(gòu)、勢(shì)能等信息，利用統(tǒng)計(jì)學(xué)原理對(duì)系統(tǒng)宏觀性質(zhì)進(jìn)行預(yù)測(cè)的一種非量子化計(jì)算方法，它具有足夠小的時(shí)間和空間尺度，能讓研究者方便地解析動(dòng)力學(xué)過(guò)程的微觀信息。然而，由于在處理作用力時(shí)所用的算法都只是近似處理，不能完全真實(shí)的描述作用力，因此非自然模擬結(jié)果時(shí)有發(fā)生［1-3］。近期，Wang等［4］研究發(fā)現(xiàn)同電荷交換介質(zhì)可以顯著提高帶同種電荷的蛋白質(zhì)復(fù)性收率，在研究中他們沒(méi)有檢測(cè)到蛋白質(zhì)被吸附到介質(zhì)表面，因此推測(cè)蛋白質(zhì)和平板之間的作用力為靜電排斥作用。但是，他們并沒(méi)有從實(shí)驗(yàn)的角度考察此過(guò)程的機(jī)理。而在利用分子模擬研究此課題時(shí)發(fā)現(xiàn)，在實(shí)驗(yàn)所用的體系中，蛋白質(zhì)和色譜基質(zhì)間相互作用力類(lèi)型很大程度上取決于所選取的力場(chǎng)參數(shù)，因此本研究就這個(gè)現(xiàn)象及其產(chǎn)生的原因進(jìn)行了分析。

分子動(dòng)力學(xué)模擬用于考察蛋白質(zhì)之間［5］以及蛋白質(zhì)和界面之間的相互作用已有很多報(bào)道。在以溶菌酶為目標(biāo)蛋白的研究中，也有利用不同的模擬方法研究其在帶電表面的行為［6-8］。但是，對(duì)于蛋白質(zhì)在由帶同種電荷的離子交換介質(zhì)介導(dǎo)的靜電作用過(guò)程鮮有研究者關(guān)注。同時(shí)，對(duì)于靜電作用力的處理大部分采用的Cut-off算法，未考慮遠(yuǎn)程靜電力可能對(duì)模擬結(jié)果帶來(lái)的影響。

因此，本研究首先利用前期構(gòu)建的靜電作用模型來(lái)模擬帶同電荷的離子交換介質(zhì)。通過(guò)分子動(dòng)力學(xué)模擬，考察不同的靜電參數(shù)對(duì)蛋白質(zhì)和帶電表面作用過(guò)程的影響，展示了蛋白質(zhì)在帶電表面的行為，計(jì)算了此過(guò)程的構(gòu)象和能量變化，初步揭示了溶菌酶體系下靜電作用力參數(shù)的選取對(duì)分子動(dòng)力學(xué)模擬的影響，揭示了蛋白質(zhì)在界面行為的微觀機(jī)理。為進(jìn)一步研究蛋白質(zhì)在界面表面的行為奠定了一定的理論基礎(chǔ)。

1 模擬方法

1.1 模擬體系構(gòu)建

根據(jù)蛋白質(zhì)數(shù)據(jù)庫(kù)（http：//www.rcsb.org/pdb/）中的晶體衍射結(jié)構(gòu)（PDB ID：3LYZ）構(gòu)建溶菌酶的全原子模型［9］，如圖1a）所示。在pH值為7的條件下，溶菌酶帶8個(gè)正電荷。由于實(shí)驗(yàn)中所用色譜介質(zhì)孔道較為復(fù)雜，為節(jié)省計(jì)算量，前人將其簡(jiǎn)化為兩互相平行的平板，包括模型化的帶電表面以及惰性平板，平板之間距離為12nm，如圖1b）所示。色譜介質(zhì)由平板及配基構(gòu)成，其中色譜介質(zhì)的拓?fù)浣Y(jié)構(gòu)由 Dundee PRODRG 2.5 server（beta）［10］生成。根據(jù)實(shí)驗(yàn)情況，每個(gè)色譜配基帶電量設(shè)置為1個(gè)正電荷，而對(duì)于模擬中配基的帶電基團(tuán)分配及每個(gè)基團(tuán)的帶電量通過(guò)以下方法獲得，首先通過(guò)Gaussian 03軟件進(jìn)行量化計(jì)算，然后通過(guò)查閱相關(guān)包含類(lèi)似結(jié)構(gòu)的分子的文獻(xiàn)［11］對(duì) Gaussian軟件計(jì)算結(jié)果進(jìn)行調(diào)整。

圖1 全原子模擬中溶菌酶和色譜孔道模型Fig.1 All-atom model of lysozyme and chromatographic pore structure.a）The lysozyme is shown by a New Cartoon model using the VMD program；b）The chromatographic pore is shown by LINES model using the VMD program

最終，1個(gè)溶菌酶分子被放置在帶電表面一側(cè)，距離平板約2nm（蛋白質(zhì)主軸與帶電平板平行）。然后將整個(gè)盒子居中放入12×12×50nm3的體系中，即在平板上下各加入19nm的真空層以防止周期性邊界條件在z方向上對(duì)模擬結(jié)果的影響。旋轉(zhuǎn)蛋白質(zhì)的方位以獲得不同的初始構(gòu)象，共進(jìn)行了3次模擬，同時(shí)進(jìn)行3次將帶電介質(zhì)所帶電荷重置為0的空白對(duì)照（包括以Cut-off和PME算法為靜電作用參數(shù)的兩組空白對(duì)照）模擬。每次模擬時(shí)間均為10 ns。

1.2 分子動(dòng)力學(xué)模擬

分子動(dòng)力學(xué)模擬利用GROMACS 4.0.5軟件［12，13］（http：//www.gromacs.org/）完成，使用的力場(chǎng)為 GROMOS96 43A1力場(chǎng)［14-15］。模擬體系為NVT系綜，通過(guò)Berendsen方法控制溫度為298.15 K，時(shí)間步長(zhǎng)為0.5 ps。勢(shì)能分析中 LJ勢(shì)能計(jì)算采用Cut-off算法，靜電作用力則分別采用Cut-off算法和PME算法，其中Cut-off算法表示在截?cái)喟霃絻?nèi)按庫(kù)侖公式計(jì)算，在截?cái)喟霃酵鈩t進(jìn)行補(bǔ)償，雖然理論而言，越大的截?cái)喟霃侥塬@到越精確的結(jié)果，但是在 Yonetani［2］的研究中發(fā)現(xiàn)，使用 Cut-off算法時(shí)，截?cái)喟霃皆O(shè)為1.8nm就能使得水分子出現(xiàn)錯(cuò)誤的分層，因此本研究的截?cái)喟霃礁鶕?jù)Gromacs手冊(cè)［14］設(shè)為1.2nm，此算法由于機(jī)理簡(jiǎn)單且運(yùn)算速度快，廣泛用于早期分子動(dòng)力學(xué)模擬；PME算法對(duì)截?cái)喟霃絻?nèi)的靜電作用力同樣采用庫(kù)侖計(jì)算，但是在計(jì)算截?cái)嗤忪o電作用力時(shí)則是先將體系分為多個(gè)小格，后利用三維正向傅里葉變換計(jì)算每個(gè)格子內(nèi)的靜電作用力，最后直接加和得到靜電作用長(zhǎng)程部分的能量［16-17］，從理論上而言計(jì)算結(jié)果更加準(zhǔn)確，同時(shí)模擬所需時(shí)間也更長(zhǎng)。

所有模擬均使用曙光TC2600刀片服務(wù)器完成。文中蛋白質(zhì)構(gòu)象圖利用VMD軟件［18］繪制。

1.3 數(shù)據(jù)分析方法

1.3.1 最小距離分析

計(jì)算蛋白質(zhì)和帶電平板質(zhì)心之間的距離（dcom）隨時(shí)間變化，質(zhì)心距離計(jì)算使用的是Gromacs自帶的g_dist程序。

1.3.2 二級(jí)結(jié)構(gòu)分析

為了描述模擬過(guò)程中蛋白質(zhì)的構(gòu)象轉(zhuǎn)化，使用GROMACS軟件包的do_dssp程序完成其二級(jí)結(jié)構(gòu)分析，采用的是Define Secondary Structure of Proteins（DSSP）方法。

1.3.3 分子間勢(shì)能分析

利用GROMACS軟件包的g_energy程序進(jìn)行勢(shì)能分析，包括 Lennard-Jones（LJ）勢(shì)能和庫(kù)侖勢(shì)能。計(jì)算庫(kù)侖勢(shì)能時(shí)，Cut-off算法直接利用庫(kù)侖公式得到所需靜電作用能和補(bǔ)償項(xiàng)相加即所得勢(shì)能；而PME算法對(duì)截?cái)嗤獾拈L(zhǎng)程靜電作用力（Coul.recip）采用傅里葉變換進(jìn)行計(jì)算，但是由于算法本身限制，這部分的靜電作用力只能以總體值來(lái)表示，因此我們采用數(shù)學(xué)方法對(duì)此總體值進(jìn)行分解，計(jì)算出蛋白質(zhì)和帶電平板之間的長(zhǎng)程靜電作用力，并和截?cái)鄡?nèi)作用力加和，進(jìn)而獲得修正庫(kù)侖勢(shì)能（Coul.fix）。

2 結(jié)果與討論

2.1 最小距離分析

利用g_dist程序計(jì)算對(duì)照組蛋白質(zhì)和帶電平板之間的距離變化，如圖2所示。

從圖2中可以看出，在空白對(duì)照組［平板電荷被重置為0，圖2a）和圖2b）］中，無(wú)論是是采用 Cutoff參數(shù)或PME參數(shù)，蛋白質(zhì)整體質(zhì)心變化并不大，蛋白質(zhì)沒(méi)有發(fā)生明顯的位移，僅僅在水的擾動(dòng)下在其初始位置附近震蕩。在平板與蛋白質(zhì)帶同種電荷的情況下，模擬結(jié)果出現(xiàn)了差別：首先，蛋白質(zhì)在兩套體系中均發(fā)生了較大位移，蛋白質(zhì)和帶電平板之間發(fā)生了相互作用，在采用 Cut-off參數(shù)的體系中，蛋白質(zhì)被吸附到平板表面，蛋白質(zhì)和平板之間最小距離隨著時(shí)間迅速減小，并穩(wěn)定在1.5nm左右；而在采用PME參數(shù)的體系中，蛋白質(zhì)被較為快速的排斥遠(yuǎn)離了帶電平板表面，最終穩(wěn)定在約7nm左右。以上結(jié)果直接表明蛋白質(zhì)在不同參數(shù)的體系中會(huì)發(fā)生不同的行為，說(shuō)明參數(shù)的選取對(duì)蛋白質(zhì)的移動(dòng)方式有決定作用。根據(jù) Wang［4］的研究結(jié)果，在實(shí)驗(yàn)過(guò)程中并未檢測(cè)到蛋白質(zhì)被吸附殘留在色譜介質(zhì)上，因此以Cut-off為靜電作用力參數(shù)的模擬體系的結(jié)果是不符合實(shí)際情況的，而以PME為靜電作用力參數(shù)的模擬結(jié)果則和實(shí)驗(yàn)結(jié)果相符。

圖2 溶菌酶和帶電平板之間質(zhì)心距離隨模擬時(shí)間的變化Fig.2 The dcom values for protein and the charged surface as a function of simulation time The average of three independent simulation trajectories is used.（A）with the cut-off parameters in control.（B）with the PME parameters in control.（C）with the cut-off parameters in charged system.（D）with the PME parameters in charged system.

2.2 蛋白質(zhì)二級(jí)結(jié)構(gòu)分析

利用DSSP方法確定蛋白質(zhì)在空白對(duì)照組和帶電體系中的二級(jí)結(jié)構(gòu)隨模擬時(shí)間的變化，如圖3所示。

結(jié)果表明，在空白對(duì)照中，無(wú)論是采用 Cut-off參數(shù)或PME參數(shù)，經(jīng)過(guò)10 ns的分子動(dòng)力學(xué)模擬，蛋白質(zhì)的二級(jí)結(jié)構(gòu)均保持的較為良好，蛋白質(zhì)沒(méi)有發(fā)生明顯的變性，這也是和實(shí)際相吻合的，即蛋白質(zhì)在純水環(huán)境中不會(huì)發(fā)生變性。而在帶電平板的體系中，情況發(fā)生了較為明顯的變化，在采用 Cut-off為靜電作用力參數(shù)的體系中，蛋白質(zhì)很快失去了大部分的二級(jí)結(jié)構(gòu)，且沒(méi)有自動(dòng)恢復(fù)部分二級(jí)結(jié)構(gòu)，說(shuō)明平板和蛋白質(zhì)分子之間的作用力使得蛋白質(zhì)在靠近平板的同時(shí)由活性態(tài)轉(zhuǎn)變?yōu)榱俗冃詰B(tài)，蛋白質(zhì)被吸附在平板表面且二級(jí)結(jié)構(gòu)被破壞；而在采用PME參數(shù)的體系中，蛋白質(zhì)的二級(jí)結(jié)構(gòu)雖然也有一定程度的改變，但并未出現(xiàn)各顏色條帶長(zhǎng)時(shí)間消失或增加的情況，這也意味著蛋白質(zhì)的活性沒(méi)有受到影響，并且僅有的二級(jí)結(jié)構(gòu)變化發(fā)生在模擬后期，結(jié)合蛋白質(zhì)質(zhì)心變化和二級(jí)結(jié)構(gòu)變化圖可以推測(cè)這是由于部分氨基酸殘基被帶電表面排斥并觸碰到了惰性平板造成的。根據(jù) Wang等［4］的研究結(jié)果，同電荷色譜介質(zhì)對(duì)蛋白質(zhì)的復(fù)性是促進(jìn)作用而非抑制作用，因此以Cut-off為靜電作用力參數(shù)的模擬體系所得結(jié)果是不符合實(shí)際情況的。

圖3 分子動(dòng)力學(xué)模擬中溶菌酶的二級(jí)結(jié)構(gòu)變化Fig.3 Secondary structure evolution of lysozyme during the molecular dynamics simulation（A）with cut-off parameters in control.（B）with PME parameters in control.（C）with cut-off parameters in charged system.（D）with PME parameters in charged system.

2.3 分子間勢(shì)能分析

通過(guò)考察蛋白質(zhì)和表面間相互作用勢(shì)能分析（在PME體系中所有庫(kù)侖勢(shì)能均為修正后的數(shù)值），考察蛋白質(zhì)在不同參數(shù)體系中構(gòu)象及位移出現(xiàn)不同的因素，如圖4所示。

圖4 蛋白質(zhì)和帶電平板間作用勢(shì)能隨模擬時(shí)間的變化Fig.4 Time evolution of the potential energy between the protein and the charged surface Average of three independent simulation trajectories is used.（A）with cut-off param eters in control.（B）with PME parameters in control.（C）with cut-off parameters in charged system.（D）with PME parameters in charged system.

圖4顯示，在空白對(duì)照中，由于蛋白質(zhì)和平板之間距離較遠(yuǎn)，故LJ勢(shì)能基本為0，由于平板電荷被重置為0，故 Coul勢(shì)能也為0，也就是說(shuō)，在這個(gè)體系中，蛋白質(zhì)和平板之間沒(méi)有相互作用。而在平板帶上正電荷之后，在采用 Cut-off參數(shù)的體系中，蛋白質(zhì)和帶電平板之間表現(xiàn)為強(qiáng)烈的靜電吸引作用（約-8 MJ/mol），且隨著時(shí)間的延長(zhǎng)，越來(lái)越多的氨基酸殘基和帶電平板之間發(fā)生了作用，能量呈小幅度波動(dòng)狀態(tài)；同時(shí)，由于靜電吸引導(dǎo)致蛋白質(zhì)分子和平板之間距離迅速減小，蛋白質(zhì)和平板之間的LJ勢(shì)能也表現(xiàn)出來(lái)，蛋白質(zhì)和平板之間出現(xiàn)了親水作用，然而其數(shù)量級(jí)僅為 kJ/mol，不起主導(dǎo)作用。而在采用PME參數(shù)的體系中，當(dāng)?shù)鞍滋幱诔跏嘉恢脮r(shí)，蛋白質(zhì)和帶電平板之間的庫(kù)侖作用能很大，但在短時(shí)間內(nèi)庫(kù)侖作用能下降到較低的水平，蛋白質(zhì)被迅速的排斥遠(yuǎn)離了帶電表面，然后呈趨向穩(wěn)定的變化趨勢(shì)，而自始至終由于蛋白質(zhì)距離平板較遠(yuǎn)，蛋白質(zhì)和帶電平板之間未產(chǎn)生明顯的LJ勢(shì)能，這說(shuō)明蛋白質(zhì)的排斥是由靜電作用力所介導(dǎo)，而且主要是遠(yuǎn)程靜電作用力。

對(duì)于采用Cut-off參數(shù)處理靜電作用力的體系，計(jì)算每個(gè)氨基酸殘基和帶電平板之間的作用能表明，較多的帶負(fù)電的天冬氨酸暴露在溶菌酶表面，使得天冬氨酸和配基發(fā)生作用的幾率很大，而帶正電的氨基酸則較多包埋在氨基酸內(nèi)部，并且Cut-off算法只計(jì)算其截?cái)喟霃絻?nèi)的靜電作用力，對(duì)于遠(yuǎn)程靜電作用力則采取簡(jiǎn)單的補(bǔ)償法進(jìn)行修正，這樣帶正電的氨基酸和配基作用的幾率就較小，因而即便溶菌酶分子總體是帶的正電荷，仍然發(fā)生了蛋白質(zhì)分子吸附到帶電平板的現(xiàn)象。反之，由于PME算法對(duì)遠(yuǎn)程作用力計(jì)算的較為精確，溶菌酶分子作為一個(gè)帶正電的整體能被帶電平板排斥遠(yuǎn)離平板表面，這也是兩種算法導(dǎo)致不同模擬結(jié)果的最直接原因。但是，需要指出的是，靜電作用參數(shù)選取對(duì)模擬結(jié)果的影響和蛋白質(zhì)本身以及其所處的體系的性質(zhì)存在很大關(guān)系。但隨著計(jì)算機(jī)硬件和軟件技術(shù)的飛速發(fā)展，使用更加精細(xì)的算法來(lái)進(jìn)行分子動(dòng)力學(xué)模擬必定是大勢(shì)所趨。

3 結(jié)論

以同電荷介質(zhì)促進(jìn)蛋白質(zhì)復(fù)性的實(shí)驗(yàn)結(jié)果為背景，采用了溶菌酶分子在靜電表面的全原子模型，通過(guò)使用不同的靜電作用力算法來(lái)模擬同電荷介質(zhì)和帶電蛋白之間的作用，考察了不同參數(shù)下模擬結(jié)果的差異及原因。結(jié)果表明，在空白對(duì)照組中，采用Cut-off算法和PME算法對(duì)模擬結(jié)果的影響較小，蛋白質(zhì)沒(méi)有失去其活性，蛋白質(zhì)均未發(fā)生明顯位移。但是，在平板帶上同種電荷之后，模擬結(jié)果出現(xiàn)了較大差異，采用Cut-off算法的體系中蛋白質(zhì)被吸附到帶電表面，蛋白的結(jié)構(gòu)及活性被強(qiáng)大的靜電吸引力破壞，這是與實(shí)驗(yàn)結(jié)論是不相符的［4］；而在采用 PME算法的體系中，蛋白質(zhì)被遠(yuǎn)程靜電作用力排斥遠(yuǎn)離了帶電平板表面，雖然蛋白質(zhì)的結(jié)構(gòu)出現(xiàn)過(guò)一部分的破壞，但是其骨架結(jié)構(gòu)依舊保存完好，且隨著時(shí)間的推移，蛋白質(zhì)的結(jié)構(gòu)出現(xiàn)了一定程度的恢復(fù)，這與 Wang［4］的實(shí)驗(yàn)結(jié)論是相吻合的，而溶菌酶分子的帶電特性以及不同算法對(duì)長(zhǎng)程靜電力的處理是造成這個(gè)現(xiàn)象的主要因素。因此，對(duì)于本體系而言，采用PME算法所得模擬結(jié)果明顯優(yōu)于采用 Cut-off算法所得模擬結(jié)果?？傮w而言，本論文通過(guò)分子動(dòng)力學(xué)模擬闡述了采取不同靜電作用力參數(shù)對(duì)模擬結(jié)果的影響，考察了此過(guò)程的微觀細(xì)節(jié)，確定了靜電作用力是導(dǎo)致蛋白質(zhì)被吸附或排斥的主導(dǎo)因素。

［1］Mark P，Nilsson L.Structure and dynamics of liquid water with different long-range interaction truncation and temperature controlmethods in molecular dynamics simulations［J］.Journal of Computational Chemistry，2002，23（13）：1211-1219

［2］Yonetani Y.A severe artifact in simulation of liquid water using a long cut-off length：Appearance of a strange layer structure［J］.Chemical Physics Letters，2005，406（1/3）：49-53

［3］Bergdorf M，Peter C.Influence of cut-off truncation and artificial periodicity of electrostatic interactions in molecular simulations of solvated ions：A continuum electrostatics study［J］.Journal of Chemical Physics，2003，119（17）：9129-9144

［4］Wang G，Dong X，Sun Y.Ion-Exchange resins greatly facilitate refolding of like-charged proteins at high concentrations［J］.Biotechnology and Bioengineering，2011，108：1068-1077

［5］白紅軍，來(lái)魯華.蛋白質(zhì)相互作用：界面分析，結(jié)合自由能計(jì)算與相互作用設(shè)計(jì) ［J］.物理化學(xué)學(xué)報(bào)，2010，26（7）：1988-1997 Bai Hongjun，Lai Luhua.Protein-Protein interactions：Interface analysis，binding free energy calculation and interaction design［J］.Acta Physico-Chimica Sinica，2010，26（7）：1988-1997（in Chinese）

［6］Ravichandran S，Madura JD，Talbot J.A brownian dynamics study of the initial stages of hen egg-white lysozyme adsorption at a solid interface［J］.The Journal of Physical Chemistry B，2001，105（17）：3610-3613

［7］Kubiak-Ossowska K，Mulheran P A.Mechanism of hen egg white lysozyme adsorption on a charged solid surface［J］.Langmuir，2010，26（20）：15954-15965

［8］Zhang L，Zhao G，Sun Y.Effects of ligand density on hydrophobic charge induction chromatography：Molecular dynamics simulation ［J］.Journal of Physical Chemistry B，2010，114（6）：2203-2211

［9］Diamond R.Real-Space refinement of the structure of hen egg-white lysozyme［J］.Journal of Molecular Biology，1974，82（3）：371-391

［10］Schuttelkop f A W，Van Aalten D M F.PRODRG：A tool for high-throughput crystallography of protein-ligand complexes［J］.Acta Crystallographica Section D：Biological Crystallography，2004，60：1355-1363

［11］Chandrasekhar I，Kastenholz M，Lins R D，et al.A consistent potential energy parameter set for lipids：Dipalmitoylphosphatidylcholine as a benchmark of the GROMOS96 45A3 force field ［J］.European Biophysics Journal with Biophysics Letters，2003，32（1）：67-77

［12］Berendsen H J C，Van Der Spoel D，Van Drunen R.GROMACS：A message-passing parallel molecular dynamics implementation ［J］.Computer Physics Communications，1995，91（1/3）：43-56

［13］Lindahl E，Hess B，Van Der Spoel D.GROMACS 3.0：A package formolecular simulation and trajectory analysis［J］.Journal of Molecular Modeling，2001，7（8）：306-317

［14］Van Gunsteren W，Billeter SR，Eising A A，et al.Biomolecular simulation：The GROMOS96 manual and user guide［B］.1996

［15］Schuler L D，Daura X，Van Gunsteren W F.An improved GROMOS96 force field for aliphatic hydrocarbons in the condensed phase［J］.Journal of Computational Chemistry，2001，22（11）：1205-1218

［16］Essmann U，Perera L，Berkowitz M L，et al.A smooth particle mesh Ewald method ［J］.The Journal of Chemical Physics，1995，103（19）：8577-8593

［17］Darden T，York D，Pedersen L.Particle mesh Ewald：An N［center-dot］log（N）method for Ewald sums in large systems［J］.The Journal of Chemical Physics，1993，98（12）：10089-10092

［18］Humphrey W，Dalke A，Schulten K.VMD：Visualmolecular dynamics［J］.Journal of molecular graphics，1996，14（1）：33-38