靜電作用力算法對分子動力學模擬結果的影響分析

李 浩，白 姝

（天津大學化工學院，天津 300072；天津大學系統生物工程教育部重點實驗室，天津 300072）

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

分子動力學模擬用于考察蛋白質之間［5］以及蛋白質和界面之間的相互作用已有很多報道。在以溶菌酶為目標蛋白的研究中，也有利用不同的模擬方法研究其在帶電表面的行為［6-8］。但是，對于蛋白質在由帶同種電荷的離子交換介質介導的靜電作用過程鮮有研究者關注。同時，對于靜電作用力的處理大部分采用的Cut-off算法，未考慮遠程靜電力可能對模擬結果帶來的影響。

因此，本研究首先利用前期構建的靜電作用模型來模擬帶同電荷的離子交換介質。通過分子動力學模擬，考察不同的靜電參數對蛋白質和帶電表面作用過程的影響，展示了蛋白質在帶電表面的行為，計算了此過程的構象和能量變化，初步揭示了溶菌酶體系下靜電作用力參數的選取對分子動力學模擬的影響，揭示了蛋白質在界面行為的微觀機理。為進一步研究蛋白質在界面表面的行為奠定了一定的理論基礎。

1 模擬方法

1.1 模擬體系構建

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

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

1.2 分子動力學模擬

分子動力學模擬利用GROMACS 4.0.5軟件［12，13］（http：//www.gromacs.org/）完成，使用的力場為 GROMOS96 43A1力場［14-15］。模擬體系為NVT系綜，通過Berendsen方法控制溫度為298.15 K，時間步長為0.5 ps。勢能分析中 LJ勢能計算采用Cut-off算法，靜電作用力則分別采用Cut-off算法和PME算法，其中Cut-off算法表示在截斷半徑內按庫侖公式計算，在截斷半徑外則進行補償，雖然理論而言，越大的截斷半徑能獲到越精確的結果，但是在 Yonetani［2］的研究中發現，使用 Cut-off算法時，截斷半徑設為1.8nm就能使得水分子出現錯誤的分層，因此本研究的截斷半徑根據Gromacs手冊［14］設為1.2nm，此算法由于機理簡單且運算速度快，廣泛用于早期分子動力學模擬；PME算法對截斷半徑內的靜電作用力同樣采用庫侖計算，但是在計算截斷外靜電作用力時則是先將體系分為多個小格，后利用三維正向傅里葉變換計算每個格子內的靜電作用力，最后直接加和得到靜電作用長程部分的能量［16-17］，從理論上而言計算結果更加準確，同時模擬所需時間也更長。

所有模擬均使用曙光TC2600刀片服務器完成。文中蛋白質構象圖利用VMD軟件［18］繪制。

1.3 數據分析方法

1.3.1 最小距離分析

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

1.3.2 二級結構分析

為了描述模擬過程中蛋白質的構象轉化，使用GROMACS軟件包的do_dssp程序完成其二級結構分析，采用的是Define Secondary Structure of Proteins（DSSP）方法。

1.3.3 分子間勢能分析

利用GROMACS軟件包的g_energy程序進行勢能分析，包括 Lennard-Jones（LJ）勢能和庫侖勢能。計算庫侖勢能時，Cut-off算法直接利用庫侖公式得到所需靜電作用能和補償項相加即所得勢能；而PME算法對截斷外的長程靜電作用力（Coul.recip）采用傅里葉變換進行計算，但是由于算法本身限制，這部分的靜電作用力只能以總體值來表示，因此我們采用數學方法對此總體值進行分解，計算出蛋白質和帶電平板之間的長程靜電作用力，并和截斷內作用力加和，進而獲得修正庫侖勢能（Coul.fix）。

2 結果與討論

2.1 最小距離分析

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

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

圖2 溶菌酶和帶電平板之間質心距離隨模擬時間的變化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 蛋白質二級結構分析

利用DSSP方法確定蛋白質在空白對照組和帶電體系中的二級結構隨模擬時間的變化，如圖3所示。

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

圖3 分子動力學模擬中溶菌酶的二級結構變化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 分子間勢能分析

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

圖4 蛋白質和帶電平板間作用勢能隨模擬時間的變化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顯示，在空白對照中，由于蛋白質和平板之間距離較遠，故LJ勢能基本為0，由于平板電荷被重置為0，故 Coul勢能也為0，也就是說，在這個體系中，蛋白質和平板之間沒有相互作用。而在平板帶上正電荷之后，在采用 Cut-off參數的體系中，蛋白質和帶電平板之間表現為強烈的靜電吸引作用（約-8 MJ/mol），且隨著時間的延長，越來越多的氨基酸殘基和帶電平板之間發生了作用，能量呈小幅度波動狀態；同時，由于靜電吸引導致蛋白質分子和平板之間距離迅速減小，蛋白質和平板之間的LJ勢能也表現出來，蛋白質和平板之間出現了親水作用，然而其數量級僅為 kJ/mol，不起主導作用。而在采用PME參數的體系中，當蛋白處于初始位置時，蛋白質和帶電平板之間的庫侖作用能很大，但在短時間內庫侖作用能下降到較低的水平，蛋白質被迅速的排斥遠離了帶電表面，然后呈趨向穩定的變化趨勢，而自始至終由于蛋白質距離平板較遠，蛋白質和帶電平板之間未產生明顯的LJ勢能，這說明蛋白質的排斥是由靜電作用力所介導，而且主要是遠程靜電作用力。

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

3 結論

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

［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］白紅軍，來魯華.蛋白質相互作用：界面分析，結合自由能計算與相互作用設計 ［J］.物理化學學報，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