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

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

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

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

1 理論和計(jì)算方法

1.1 計(jì)算方法

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

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

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

(1)

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

(2)

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

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

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

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

表2 POPE分子的鍵長(zhǎng)伸縮振動(dòng)參數(shù)

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

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

2 結(jié)果與討論

2.1 小分子的幾何構(gòu)型和電荷分布

應(yīng)用公式(1)計(jì)算模型小分子MAS、DMP和ETMA的結(jié)構(gòu)和電荷分布,獲得的穩(wěn)定結(jié)構(gòu)與QM結(jié)果進(jìn)行比較.表5給出了3種分子在兩種方法下的鍵長(zhǎng)、鍵角以及二面角的平均偏差(AAD)、均方根偏差(RMSD)和相對(duì)均方根偏差(RRMSD).其中,偏差最大的DMP的鍵長(zhǎng)、鍵角、二面角的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力場(chǎng)方法計(jì)算結(jié)果與QM方法結(jié)果有較好的一致性.

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

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

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

2.2 POPE分子的結(jié)構(gòu)和電荷分布

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

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

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

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

3 結(jié)論與展望

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