磷脂分子POPE的浮動電荷分子力場

宮利東, 趙力儐, 田 博, 楊忠志

(遼寧師范大學 化學化工學院，遼寧 大連 116029)

生物膜是細胞的重要組成成分,它作為細胞內(nèi)和細胞外環(huán)境之間的屏障,提供專門的物理環(huán)境實現(xiàn)物質轉運和信號傳導,這些功能需要動態(tài)的流體環(huán)境,允許膜蛋白的旋轉和平移擴散[1]. 細胞膜主要由磷脂構成,在原子和分子水平上理解磷脂雙層的結構和動態(tài)特性是至關重要的[2-3].對于脂質雙層的研究已有近百年的歷史[4],主要實驗技術有量熱法(用于研究膜相行為)、熒光共振能量轉移(FRET)、原子力顯微鏡(AFM)、核磁共振(NMR)、小角中子散射(SANS)、小角X射線散射(SAXS)、準彈性中子散射(QENS)和質譜等.隨著計算機技術的迅速發(fā)展,分子動力學模擬已被廣泛應用于探討各種復雜體系的結構和性質,與實驗技術相比,分子動力學模擬能夠得到體系的原子級別的詳細結構信息以及實驗中難以捕獲的微觀細節(jié).但其計算的成敗很大程度上取決于力場的適用性、計算速度的快慢和計算方法的正確性.大多數(shù)的脂質力場采用固定電荷模型處理靜電相互作用,比如AMBER[5]、CHARMM[6]以及OPLS-AA[7]和GROMOS[8]等力場.雖然固定電荷力場的計算成本相對較低,可以對較大體系進行長時間尺度模擬,但忽略了極化效應[9].由于磷脂分子的特殊結構,其嵌入式通過膜的外部和內(nèi)部的環(huán)境是不同的,需要明確考慮極化作用[10].與傳統(tǒng)力場相比,可極化力場模型可以表現(xiàn)體系的電荷分布或偶極矩等隨環(huán)境和結構的改變而發(fā)生的相應變化,因此，建立磷脂分子的可極化力場近年來受到關注.目前，主要有3種可極化模型:誘導偶極模型[9-11]、Drude諧振子模型[1,10,12-13]和浮動電荷模型[14-17].李國輝等人報道了基于誘導偶極模型的AMOEBA力場、磷脂雙層的1,2-二油酰基-磷酸膽堿(DOPC)、1-棕櫚酰-2-油?；?磷脂酰乙醇胺(POPE)、1-棕櫚酰-2-油?；?sn-磷脂酰膽堿(POPC)和1,2-二肉豆蔻?；?sn-磷脂酰膽堿(DMPC)的可極化力場[9,18];Roux課題組[1,13]開發(fā)了1,2-二棕櫚酰基-sn-甘油-3-磷酸膽堿(DPPC)、二月桂酰磷脂酰膽堿(DLPC)、二棕櫚酰磷脂酰乙醇胺(DPPE)和1,2-二油?；?sn-甘油基-3-磷酸乙醇胺(DOPE)、DMPC、POPE、POPC、DOPC等多種飽和及不飽和的兩性磷脂分子體系的Drude可極化力場.

在密度泛函理論下的電負性均衡方法的基礎上,楊忠志等創(chuàng)建和發(fā)展了原子-鍵電負性均衡浮動電荷分子力場方法(ABEEMσπ/MM),在對純水、有機分子與聚合物、電解質溶液、核酸以及多肽蛋白質等體系的研究中,均獲得了優(yōu)于傳統(tǒng)力場的結果[17].在本課題組之前建立的DPPC分子的ABEEMσπ/MM力場的基礎上[19],繼續(xù)開展對POPE的研究,為建立廣泛適用于磷脂體系ABEEMσπ/MM力場打下基礎.

1 理論和計算方法

1.1 計算方法

量子化學計算采用Gaussian 09[20]程序.將POPE分子的極性頭部拆分為乙基三甲基銨(ETMA)、乙酸甲酯(MAS)、磷酸二甲酯(DMP)等小分子片段,選用MP2/6-31G*水平對中性和帶正電的片段進行結構優(yōu)化,而對于帶負電的片段使用MP2/6-31+G(d)方法.對于POPE分子,為節(jié)省計算量,選用M062X泛函.其中,極性頭部和尾部烴鏈分布選用6-31++G(d,p)和6-31G(d)基組進行結構優(yōu)化和頻率計算.對于電荷分布的計算,基于本課題組之前的研究[16],應用較大基組會高估兩原子間極化作用,采用HF/STO-3G的方法.

1.2 POPE的ABEEMσπ/MM浮動電荷勢能函數(shù)

根據(jù)ABEEMσπ/MM模型,構建了POPE分子的浮動電荷勢能函數(shù):

(1)

前2項分別為磷脂分子的鍵長伸縮和鍵角的彎曲振動勢能,采用諧振子勢能函數(shù)形式;第3項為二面角的扭轉振動勢能,采用Fourier級數(shù)形式;第4項為非共面二面角勢能;第5項為范德華相互作用勢,采用Lennard-Jones12-6形式;最后1項為靜電相互作用.其中,i、j表示不同分子,a、b表示不同原子,kr和fθ為力常數(shù),r、θ和φ為體系的真實鍵長、鍵角和二面角值,req、θeq為平衡時的鍵長、鍵角值.V1、V2、V3表示二面角的Fourier展開系數(shù).σia,jb、εia,jb分別表示第i個分子中原子a和第j個分子中原子b之間的碰撞直徑和勢阱深度,fij為常數(shù).qm和qn分別表示第m和第n個位點所帶電荷.kmn為靜電相互作用校正系數(shù),①當位點m和n處于1-2或1-3關系時,kmn=0;②當兩位點處于氫鍵作用區(qū)域(HBIR)時,kmn=kHB,kHB為氫鍵擬合函數(shù),引入它可以防止2個原子之間的電荷過度極化,其具體公式如下所示:

(2)

其中，A、B、C、D為相關參數(shù),Rlp,H為孤對電子和氫原子間的距離.③當兩位點在其他區(qū)域時,kmn=0.57.ABEEMσπ/MM浮動電荷力場將分子劃分為原子、鍵以及孤對電子位點,并且各位點電荷會隨環(huán)境的改變而發(fā)生變化.

1.3 ABEEMσπ/MM參數(shù)

ABEEMσπ/MM參數(shù)包括力場參數(shù)和電荷參數(shù).力場參數(shù)中的鍵長伸縮振動、鍵角彎曲振動等參數(shù)來自CHARMM36力場.電荷參數(shù)包含價態(tài)硬度和價態(tài)電負性.根據(jù)QM計算得到的模型分子ETMA、MAS、DMP的穩(wěn)定結構和電荷分布,調節(jié)優(yōu)選相關參數(shù),使由ABEEMσπ/MM勢能函數(shù)公式計算得到的模型分子結構、電荷分布等與QM結果吻合,優(yōu)選后的參數(shù)列在表1中.進一步應用ABEEM/MM勢能函數(shù)計算磷脂分子POPE的結構和電荷分布,檢驗參數(shù)的合理性.優(yōu)選得到的POPE力場參數(shù)分別列在表2～表4中.

表1 小分子MAS、DMP、TMA的靜電參數(shù)

表2 POPE分子的鍵長伸縮振動參數(shù)

表3 POPE分子的鍵角彎曲振動參數(shù)

表4 POPE分子的二面角扭轉勢能參數(shù)

2 結果與討論

2.1 小分子的幾何構型和電荷分布

應用公式(1)計算模型小分子MAS、DMP和ETMA的結構和電荷分布,獲得的穩(wěn)定結構與QM結果進行比較.表5給出了3種分子在兩種方法下的鍵長、鍵角以及二面角的平均偏差(AAD)、均方根偏差(RMSD)和相對均方根偏差(RRMSD).其中,偏差最大的DMP的鍵長、鍵角、二面角的AAD、RMSD、RRMSD分別為0.021 9、0.003 5 nm和2.59%,2.39°、3.48°和3.16%,2.61°、3.26°和2.87%.總體來說,ABEEMσπ/MM力場方法計算結果與QM方法結果有較好的一致性.

表5 ABEEMσπ/MM與QM方法優(yōu)化的模型小分子的結構對比

圖1為ABEEMσπ/MM方法和QM方法計算得到的3種分子電荷分布的線性相關圖,可以看出電荷分布的線性相關系數(shù)為0.972 8,斜率接近于1,截距接近于0,兩者獲得的電荷分布具有良好的一致性.

圖1 ABEEMσπ/MM與QM計算的DMP、MAS和 ETMA的電荷分布的線性相關圖Fig.1 The linear correlation of charge distributions of DMP,MAS and ETMA between ABEEMσπ/MM and QM

2.2 POPE分子的結構和電荷分布

表6為ABEEMσπ/MM與QM方法得到的POPE分子的鍵長R(nm)、鍵角θ(°)、二面角τ(°)的比較.結果表明,ABEEMσπ/MM與QM結果符合很好.其中，二面角的AAD最大為0.344 0°,RMSD為0.507 2°,RRMSD為0.452 0%.

表6 ABEEMσπ/MM與QM方法得到的POPE分子穩(wěn)定結構的對比

圖2為ABEEMσπ/MM與QM方法計算得到的POPE分子的電荷分布的線性相關圖.從圖中可以看出二者線性相關很好,斜率為1.024 0,線性相關系數(shù)為0.981 3,截距趨近于0.

圖2 ABEEMσπ/MM與QM計算POPE電荷分布的線性相關圖Fig.2 The linear correlation of charge distributions of POPE between ABEEMσπ/MM and QM

3 結論與展望

建立了POPE分子的ABEEMσπ/MM浮動電荷勢能函數(shù),基于QM計算結果,優(yōu)選確定相關參數(shù).應用該勢能函數(shù)計算獲得的POPE分子的結構和電荷分布與QM結果符合很好,表明本文初步建立的POPE分子的ABEEMσπ/MM力場的合理性和可靠性.此項研究為進一步將ABEEMσπ/MM推廣到其他磷脂分子體系,建立可用于復雜生物分子體系的磷脂分子的浮動電荷分子力場打下了基礎.