β-環(huán)糊精抑制α-半乳糖苷酶的動力學模擬

鄒 偉，王 敏，陶佳妮，姚 娣，李 赟，王雪松，沈汪洋,2*

(1.武漢輕工大學 食品科學與工程學院，湖北 武漢 430023； 2.大宗糧油精深加工省部共建教育部重點實驗室，湖北 武漢 430023)

β-環(huán)糊精抑制α-半乳糖苷酶的動力學模擬

鄒 偉1，王 敏1，陶佳妮1，姚 娣1，李 赟1，王雪松1，沈汪洋1,2*

(1.武漢輕工大學 食品科學與工程學院，湖北 武漢 430023； 2.大宗糧油精深加工省部共建教育部重點實驗室，湖北 武漢 430023)

為進一步了解β-環(huán)糊精與α-半乳糖苷酶的作用機制，在前、中期采用同源建模、分子對接的基礎上，采用對復合物體系進行分子動力學模擬，并計算出均方根偏差、溫度、密度以及總能量四個衡量指標，對體系動力平衡的穩(wěn)定程度進行評價，得到了動力學模擬軌跡是穩(wěn)定的結果，再采用MMPBSA能量計算及分解的方法，對β-環(huán)糊精與α-半乳糖苷酶之間的相互作用進行能量計算，確定了兩者的結合是一個自發(fā)的過程，且能夠達到一個穩(wěn)定的狀態(tài)，并得到了β-環(huán)糊精對α-半乳糖苷酶的抑制作用主要是由范德華相互作用能和靜電作用能產生的結論。

β-環(huán)糊精；α-半乳糖苷酶；分子動力學模擬；MMPBSA能量計算

1 引言

β-環(huán)糊精(β-cyclodextrin)是由7個α-D-吡喃葡萄糖單元通過α-1,4糖苷鍵連接而成的環(huán)狀低聚糖，具有一個環(huán)外親水、環(huán)內疏水且具有一定尺寸的立體手性空腔[1]。α-半乳糖苷酶(α-galactosidase, EC 3.2.1.22)是一種外切糖苷酶，主要作用對象是具有α-半乳糖苷結構的碳水化合物[2-3]。

咖啡豆α-半乳糖苷酶是α-半乳糖苷酶的一種,從咖啡豆中提取得到[4-7]，咖啡豆α-半乳糖苷酶能把α-半乳糖基直接轉移到β-環(huán)糊精的環(huán)上，形成新的分支環(huán)糊精：α-半乳糖基-β-環(huán)糊精[8-9](α-galactosyl-β-cyclodextrin)。1992年首次報道酶法合成了這種分支環(huán)糊精，目前咖啡豆α-半乳糖苷酶酶法合成是制備α-半乳糖基-β-環(huán)糊精主要的手段。α-半乳糖基-β-環(huán)糊精具有較高的水溶性，其分支上的半乳糖基在動物組織中具有特殊的識別受體，具有廣闊的應用前景。

目前有關α-半乳糖基-β-環(huán)糊精的合成及應用報道并不多，主要原因是酶法合成α-半乳糖基-β-環(huán)糊精的產率不高(22 %左右)[10-12]，對分子機制了解的缺乏嚴重制約了該領域的深入研究。有研究發(fā)現(xiàn)，環(huán)糊精可以改變酶促反應的進程，可以與酶分子的特殊基團形成復合物，激活或抑制酶的活性，從而影響反應速率或反應平衡[13-14]。而酶法合成α-半乳糖基-β-環(huán)糊精的方法中，β-環(huán)糊精為合成反應的底物，α-半乳糖苷酶(咖啡豆來源)為反應酶。在反應過程中，環(huán)糊精和咖啡豆α-半乳糖苷酶共存于同一個反應體系中，可能是環(huán)糊精自身抑制了咖啡豆α-半乳糖苷酶的活性，從而導致酶法合成產率降低。

本實驗從這一現(xiàn)象入手，前期采用同源建模的方法[15-17]，合理構建了咖啡豆α-半乳糖苷酶的結構模型，中期采用分子對接的方法，將α-環(huán)糊精、β-環(huán)糊精及γ-環(huán)糊精與咖啡豆α-半乳糖苷酶分別進行對接，結果表明環(huán)糊精抑制了酶的活性，β-環(huán)糊精的抑制作用較為明顯。本文使用MMPBSA方法計算β-環(huán)糊精與α-半乳糖苷酶的結合自由能，研究二者的結合模式，探索β-環(huán)糊精抑制α-半乳糖苷酶的作用機制，可以為降低抑制效應，提高合成效率提供理論依據(jù)。

2 材料與方法

2.1 材料與儀器

分子動力學軟件包Amber 12；分子顯示軟件PyMOL 1.7；計算機Intel Xeon E5-2650 v2

2.2 實驗方法

2.2.1 分子動力學模擬

動力學模擬的蛋白質初始構象來自實驗前期同源建模的咖啡豆α-半乳糖苷酶的結構模型，在AmberTools 12的Tleap模塊中補全缺失的氫原子，由ff 10力場產生立場參數(shù)。小分子的初始構象來自蛋白質數(shù)據(jù)庫的晶體結構(3CGT)中的配體β-環(huán)糊精, 在實驗的前期工作中已對接完畢，用sqm模塊中的AM1-BCC計算分配原子電荷。將復合物置于截角八面體水盒子中心，水分子采用TIP3P模型，水盒子壁到復合物最近原子的距離為10.0 ?。加入5個Na+保持模擬體系呈電中性。

采用Amber 12的PMEMD模塊對模擬體系進行四個步驟的操作: (1)能量最小化:每一步優(yōu)化先執(zhí)行10000步的最陡下降優(yōu)化，再執(zhí)行20000步的共軛梯度優(yōu)化；(2)升溫:在50 ps內將體系溫度從0K升到300 K；(3)平衡:執(zhí)行50 ps的常溫300 K的溫度平衡；(4)動力學模擬:對體系進行6 ns的無約束分子動力學模擬，時間步長為2 fs，坐標文件每1 ps記錄一次，用于后續(xù)分析。

2.2.2 計算結合自由能

自由能計算由Amber 12的MMPBSA模塊完成，從分子動力學模擬的后2 ns每10 ps取一次構象，分子動力學軌跡共200個構象用來做自由能計算。

受體和配體的結合自由能由如下方程計算：

ΔGBind=ΔGComplex-ΔGReceptor-ΔGLigand

其中ΔGBind是結合自由能，ΔGComplex,ΔGReceptor和ΔGLigand分別是復合物，受體和配體的自由能。用MMPBSA方法把它們分解為三個部分：

ΔG=ΔEgas+ΔEsol

其中為ΔEgas是氣相能量焓，ΔEsol是溶解自由能。

氣相能量焓ΔEgas可分解為如下兩項：

ΔEgas=ΔEvdw+ΔEele

其中為ΔEvdw是范德華相互作用能，ΔEele是靜電相互作用能。

溶解自由能ΔEsol可分解為如下兩項：

ΔEsol=ΔEcal+ΔEsurf

其中為ΔEcal是極性溶解自由能，ΔEsurf是非極性溶解自由能。

ΔEvdw和ΔEele即MM部分；ΔEcal可通過泊松-玻爾茲曼(Poisson-Bohzmann)方程求解，即PB部分；ΔEsurf由如下經驗公式計算：

ΔEsurf=γSASA +β

其中SASA是溶劑可接觸表面積(solvent-accessible surface)，即SA部分。

3 結果與討論

3.1 動力學模擬體系的平衡

對復合物體系進行了6 ns的分子動力學模擬，為了評價體系動力學平衡的穩(wěn)定程度，計算了均方根偏差(Root-mean-square deviation,RMSD)、溫度、密度以及總能量這四個衡量指標,分別見圖1、圖2、圖3及圖4。

圖1 MD模擬中體系的均方根偏差隨時間變化圖

圖2 MD模擬中體系的溫度隨時間變化圖

圖3 MD模擬中體系的密度隨時間變化圖

圖4 MD模擬中體系的總能量隨時間變化圖

從3.5 ns后，體系趨于穩(wěn)定，在3.5～6 ns內均方根偏差RMSD未見明顯波動，收斂于2.74 ?，漲落范圍低于0.31 ?，體系達到了動態(tài)平衡。在整個模擬過程中體系的溫度、密度及總能量隨時間的變化曲線較為穩(wěn)定，無異常變化。上述結果表明用于下一步計算分析的動力學模擬軌跡的穩(wěn)定性是可靠的。

3.2 MMPBSA能量計算及分解

采用MMPBSA方法，使用動力學模擬軌跡后2 ns的200幀構象，對β-環(huán)糊精與α-半乳糖苷酶之間的相互作用進行能量計算，結合自由能各成分的貢獻均在表1中列出。

表1 MMPBSA計算結合自由能

組成成分復合物受體配體差值均值均值均值均值范德華作用能-1681．2489-1683．4017-9．8159-58．0313靜電作用能-15302．5695-15063．6189-183．8766-55．0741氣象能量焓-12312．0024-12144．1514-53．2335-144．6175非極性溶劑化能60．779861．95104．2188-5．3899極性溶劑化能-3581．3175-3622．9516-44．415286．0492溶解自由能-3520．5377-3561．0006-40．196480．6593結合自由能(PB法)-15832．5401-15705．1520-93．4299-33．9582結合自由能-33．9582

VDW為范德華作用能；ELE為靜電作用能；GAS為氣象能量焓，GAS=VDW+ELE；PBSUR為非極性溶劑化能；PBCAL為極性溶劑化能；PBSOL為溶解自由能，PBSOL=PBSUR+PBCAL；PBTOT為用PB方法計算的結合自由能，PBTOT= GAS+ PBSOL；ΔGBind為受體與配體的結合自由能，ΔGBind= PBTOT-TSTOT。因為該體系是一個能夠達到動態(tài)平衡的狀態(tài)，反應的熵變?yōu)榱悖宰詈蟮玫降慕Y合自由能ΔGBind即為PBTOT的值-33.9582。結合自由能是一個負值，表明兩者結合的過程為一自發(fā)過程，且能夠達到一個穩(wěn)定的狀態(tài)。從表1中可以看到體系中的范德華相互作用能有利于體系的穩(wěn)定，并且相對于其它的能量項數(shù)值最小，因此為包合的主要驅動力。體系的靜電作用能也有利于二者結合的形成，但相對于范德華相互作用要弱一些，可能是極性的溶解自由能抵消了靜電能，有利于包合的這種趨勢。非極性溶解自由能為負值，說明疏水效應也有利于兩者結合的穩(wěn)定性，但其絕對值較小，因此也不是形成結合的主要驅動力。

4 結論

在本研究中對β-環(huán)糊精與α-半乳糖苷酶復合物體系進行了分子動力學模擬，得到了動力學模擬軌跡是穩(wěn)定的結果，再采用MMPBSA能量計算及分解的方法，對β-環(huán)糊精與α-半乳糖苷酶之間的相互作用進行能量計算，得到PBTOT的值-33.9582，確定了兩者的結合是一個自發(fā)的過程，且能夠達到一個穩(wěn)定的狀態(tài)，并得到了β-環(huán)糊精對α-半乳糖苷酶的抑制作用主要是由范德華相互作用能和靜電作用能產生的。這為β-環(huán)糊精與α-半乳糖苷酶的相互作用機制的研究和開發(fā)該酶的一些特殊作用提供了相應的依據(jù)。

[1] 宮斌, 黃據(jù)芹, 王懷友. β-環(huán)糊精對中性紅及熒光素包合性能的研究 [J].化學分析計算, 2010,02:19-23.

[2] 郝桂娟, 張凱, 王學智. α-半乳糖苷酶的研究進展 [J].中國畜牧獸醫(yī), 2013, 03:149-154.

[3] 陶佳妮, 李赟, 姚娣. 咖啡豆α-半乳糖苷酶的同源建模及分析 [J].食品工業(yè)科技, 2016,01:326-330.

[4] Marraccini P,Rogers W J. Biochemical and molecular characterization of from coffee beans[J].Plant Physiology and Biochemistry, 2005,43: 909-920.

[5] 沈汪洋, 金征宇. 發(fā)芽咖啡豆α-半乳糖苷酶的性質研究 [J]. 安徽農業(yè)科學, 2011,19:11958-11960.

[6] 侯廷軍, 徐筱杰. 基于分子動力學模擬和連續(xù)介質模型的自由能計算方法 [J]. 化學進展, 2004,02:153-158.

[7] 沈汪洋, 金征宇. 植物多酚對α-半乳糖苷酶提取的影響 [J]. 食品工業(yè)科技, 2008,04:90-92.

[8] 沈汪洋, 金征宇, 高虹. α-半乳糖苷酶法合成α-半乳糖基-β-環(huán)糊精的研究 [J]. 湖北農業(yè)科學, 2011,19:4041-4043.

[9] 沈汪洋, 金征宇. 半乳糖基-β-環(huán)糊精的合成與結構研究 [J]. 食品科學, 2011,12:57-59.

[10] 崔波, 金征宇. 麥芽糖基(α-1-6)β-環(huán)糊精的酶法合成和結構鑒定 [J]. 高等學?；瘜W學報, 2007,02:283-285.

[11] Harak , Fujitak , Kuwaharad, et al. Galaetosylation of cyclodextrins and branched cyclodextrins by α-galactosidase [J]. Bioscienee, Bioteelmology and Biochem istry, 1994, 4:652—659.

[12] Koizumi K, Tanimoto T, Fujita K. et al. Preparation, iIsolation, and characterization of novel heterogeneous branched cyclomalto-oligosaccharides having β-galactosyl residue(s) on the side chain [J].Carbohydrate Research, 1993,238:75一91.

[13] Seosj, Kimsh, Sasagawat, et al. Delivery of all trans-retinoic acid(RA) to hepatocyte cell line from RA/galactosyl α-cyclodextrin inclusion complex [J]. European Journal of Pharm aceutics and Biopharm aceutics, 2004,58: 681-687.

[14] 于博, 張煥新, 金征宇等. 環(huán)糊精及其衍生物對普魯蘭酶抑制機理的研究 [J]. 食品工業(yè)科技, 2011,03:109-111+115.

[15] 洪瑋. 糖苷水解酶Fncel5A和GB1多肽的分子動力學模擬 [D].合肥：中國科學技術大學生物信息學，2012.

[16] 王加磊, 孟現(xiàn)美, 張慶剛. HIV-1蛋白酶異位抑制劑體系的長時問分子動力學模擬 [J]. 山東科學, 2013,02:1-6.

[17] 柯國濤. 用分子模擬方法研究HIV-1整合酶與病毒DNA的相互作用 [D].北京：北京工業(yè)大學生物物理學.2008.

Dynamics simulations of inhibition of β-cyclodextrin on α-galactosidase

ZOUWei1,WANGMin1,TAOJia-ni1,YAODi1,LIYun1,WANGXue-song1,SHENWang-yang1,2*

(1.School of Food Science and Engineering, Wuhan Polytechnic University, Wuhan 430023,China; 2. Key Laboratory of the Deep Processing of Bulk Grain and Oil Authorized by Ministry of Education, Wuhan 430023,China)

In order to further understand the mechanism of action of β-cyclodextrin and α-galactosidase, this test used homology modeling and molecular docking based on the molecular dynamics simulation of complex system in the early and middle stage, and calculated four indicators include rmsd deviation, temperature, density and total energy, and evaluated the degree of stability of the system dynamic equilibrium, aquired a stable result of dynamics simulation trajectory , then by using the MMPBSA energy calculation and decomposition method for calculating the energy of interaction between β-cyclodextrin and α-galactosidase, determined the combination of β-cyclodextrin and α-galactosidase was a spontaneous process, and could achieve a steady state, and aquired the result that the inhibitory effect of β-cyclodextrin on α-galactosidase was mainly composed of Fan Dehua interaction energy and the electrostatic interaction.

β-cyclodextrin；α-galactosidase；molecular dynamics simulation；MMPBSA energy calculation

2016-12-23

鄒偉(1992-)，男，碩士研究生。E-mail：530294379@qq.com.

*通信作者：沈汪洋(1978-)，男，副教授，博士。E-mail: whwangyangshen@126.com

國家自然科學基金委員會；國家自然科學基金青年科學基金項目(31301415)

2095-7386(2017)01-0022-04

10.3969/j.issn.2095-7386.2017.01.004

TS 201.2

A