基于分子模擬技術研究來源于嗜熱地衣芽孢桿菌葡聚糖酶的活性位點

禹亞男，李 婷

(重慶醫科大學附屬第二醫院 醫學重癥監護科，重慶 400000)

嗜熱酶是從嗜熱微生物中分離得到的一類熱穩定性酶，具有化學催化劑無法比擬的優點，尤其是在高溫條件下保持的極好穩定性，克服了常溫酶在應用中出現的化學性質不穩定的現象，使很多高溫化學反應得以實現，從而將極大的促進生物技術產業的發展。葡聚糖酶(glucanase)能專一性催化葡聚糖或地衣多糖中臨近β-1,3鍵、β-1,4鍵水解，產生3- 5個葡萄糖單位的低聚糖及葡萄糖，從而降低葡聚糖的親水性和粘性[1]，因此現被很多國家用來作為飼料添加劑[2]及啤酒生產的酶制劑[3]。對于來自不同物種的糖苷水解，其葡聚糖酶的氨基酸組成和的活性部位存在差異。例如，姜彤等研究來源于柄孢霉的β-葡聚糖酶晶體屬于GH131糖苷酶家族，由外圍的9個β折疊和內部的6個β折疊構成具有典型的β-jelly-roll結構，表面呈凹槽結構，易于和底物結合。其中E99和E139作為酸堿供體，是保守的催化結合位點[4]。來源于嗜熱擬青霉的Ptlic16A晶體屬于糖苷水解酶GH16家族，其催化結構為β-jelly-roll，其中E113作為親核劑，E118作為質子供體，W97、W101、W253和底物的識別有關[5]。來源于黃孢原毛平革菌的Lam55A葡聚糖酶屬于GH55家族，其具有兩個平行的右手β折疊結構域構成的一個類似于胸腔的結構，其催化位點位于兩個結構域之間[6]。植物葡聚糖酶主要屬于GH17家族，其結構具有典型的TIM桶結構(TIM-barrel fold)，其酸堿催化、親核體分別位于第4和第7 β折疊鏈上[7]。來源于GH64家族的LPHase葡聚糖酶晶體其催化部位由一個β折疊和一個α/β結構域構成，在兩個結構域之間形成了寬廣的溝槽[8]。目前尚沒有關于來源嗜熱地衣芽孢桿菌葡聚糖酶晶體結構的報道。有研究表明，葡聚糖酶的最適溫度為55℃，在溫度100℃下，30 min仍有70%的生物活性并且可以催化多種糖苷鍵水解[9]，是一種新型的葡聚糖酶。本研究在此研究的基礎之上，通過利用分子對接與分子動力學模擬技術研究葡聚糖酶與二聚糖UNK之間分子的識別機制，以期為嗜熱地衣芽孢桿菌葡聚糖酶分子結構后續的理性改造和高效利用提供有力依據。

1 材料與方法

1.1 葡聚糖酶結構的制備

利用在線服務器I-TASSER(http://zhanglab.ccmb.med.umich.edu/I-TASSER/)對嗜熱地衣芽孢桿菌葡聚糖酶(登錄號：Y_006712752.1)進行同源建模獲的試驗模板，通過GROMACS4.6.5軟件[10]對其進行結構優化，獲得最佳的構象。

1.2 配體分子的制備

利用ChemOffice 3D軟件[11]構建二聚糖分子結構，用ChemOffice 3D的MM2模板對底物進行能量最小化，保存結構為pdbqt格式，以便進行進一步的對接工作。

1.3 分子對接

采用AutoDock4.2程序包[12]對二聚糖和葡聚糖酶進行預處理。對接時，在格點盒子中，X維、y維、z維格點數目分別設為52、54、54.該格點盒子的中心坐標為-24.855，14.414，93.059，格點間距為0.0392 nm。采用拉馬克遺傳算法(LGA)來尋找配體和受體最佳的結合位置，將局部能量搜索與遺傳算法相結合，以半經驗勢函數為打分函數，對小分子構象和位置進行全局搜索。為了使結果更加準確，采用遺傳算法中較大的運算參數，種群個數設為150，能量評估的最大數目設為1×107，運算循環數設為1000，其他的使用該程序默認的參數。同時設立2組平行組。最后依據獲得的1000個配體小分子構象進行聚類分析(參數為0.2 nm)，根據聚類的情況、結合能以及構象的合理性選取最佳構象進行分析并作為分子動力學模擬的初始構象，同時使用acpype程序[13]生成小分子配體的拓撲文件，用于構建分子動力學模擬中的復合物拓撲文件。

1.4 分子動力學模擬

本實驗采用gromacs4.6.5軟件，使用gromacs43 al力場[14]，水分子采用SPCE模型[15]，遠程靜電作用采用PME方法[16]，系統的溫度和壓力通過V-rescale和Parrinel-Rahman實現控制[17]。將復合物置于一個立方體水盒子中，使用周期性邊界條件以消除邊界效應，加入離子中和電荷使整個體系處于中性環境，使用5000步最陡下降法進行能量最小化，然后進行400 ps的限制性動力學模擬，使得溫度分別達到300 K，壓力為1個大氣壓，溫度和壓力耦合系數為0.2 ps[18]，待體系平衡后，進行50 ns的自由動力學模擬，模擬中采用的積分步長為2 ps，每2 ps保存一次軌跡結構用于后期數據的處理。通過分析復合物root mean squared deviation (rmsd)值，系統能量來判斷復合物構象的穩定性。

1.4.1 底物和酶之間接觸數和接觸距離的計算

底物和葡聚糖酶之間的接觸數和接觸距離均采用GROMACS自帶程序g_mindist進行計算，其中計算接觸數時所用的截斷距離為0.6 nm[19]。

1.4.2 復合物自由能計算

采用單軌跡方案的MM-PBSA方法計算葡聚糖酶與底物小分子的結合自由能，此方法是從復合物的MD模擬平衡的后5 ns軌跡中每隔一定的時間間隔取出體系結構，通過下式計算出平均結合自由能。

△Gbind=△Emm+△Gsol-T△S

其中△Gbind是結合自由能，△Emm、△Gsol和T△S分別是氣相中的分子力學能、溶劑自由能和熵變對結合自由能的貢獻。

△Emm=△Eele+△Evdw

式中△Eele和△Evdw分別表示氣相中靜電相互作用和范德華能。溶劑自由能也由如下的兩部分組成：

△Gsol=△Gel+△Gnonel

式中△Gel和△Gnonel分別代表極性溶劑自由能和非極性溶劑自由能。前者的成分可以通過求解泊松—波爾茲曼方程獲得[20]，溶質和溶劑的介電常數分別設為1.0和80.0。后一項由下列的經驗方程求解：△Gnonel=γSASA+β

式中△Gnonel代表非極性溶劑自由能，SASA代表溶劑可及接觸面積，γ和β的值分別取為2.27 kJ/(mol·nm2)和3.85 kJ/mol。至于熵效應，可用簡正模分析來計算，此項的計算非常耗時且對結果影響較小，因此一般情況下，可忽略熵變效應[21]。

1.4.3 氫鍵分析

本文中小分子配體和葡聚糖酶之間的氫鍵利用GROMACS自帶程序g_hbond進行分析。當供體原子D和受體原子A之間的距離小于0.35 nm，且D-H-A之間的夾角大于120°時，即認為供體原子D和受體原子A之間形成氫鍵。

2 結果與分析

2.1 分子對接

圖1 復合物三維結構和平面圖

通過利用分子對接方法，對二聚糖小分子配體與葡聚糖酶進行對接，選取自由能結合最小的構象，其試驗結果如圖1所示，圖1a中，我們可以看出，Thr126、Thr124、Asp71、Ala72、His73、Gly46和Lys47殘基參與二聚糖小分子的結合，并且其結合部位與Cheng等人研究的來自嗜熱擬青霉的Ptlic16A晶體結合部位一致[5]。通過使用LigPlot軟件[22]，獲得復合物結合部位平面圖1b，從圖中，我們可以看出，Ala72、Gly46和Gly125殘基分別與小分子配體形成疏水相互作用，Thr124、Lys47、His73、Asp71和Lys126殘基分別參與復合物氫鍵的形成。

2.2 分子動力學模擬

2.2.1 均方根偏差分析

從圖2中，我們可以看出，葡聚糖酶和二聚糖分子結構從5 ns后趨于平衡，其中分別圍繞著0.3 nm和0.15 nm值上下波動，在45～50 ns的時間內，二聚糖小分子結構的rmsd波動較小，并且系統總能量的平均波動量低于0.8%，因此可以認為系統在MD模擬后比較穩定，用此軌跡分析得出的結果是可靠的。

圖2 酶和底物的RMSD值的變化

2.2.2 二聚糖和葡聚糖酶之間的接觸數和分子間的距離

從圖3A中我們可以看出，小分子配體與蛋白質之間的接觸數在前10 ns內，隨著時間的增加，其接觸數目增多，在30～50 ns內，其值趨于平衡，總體波動范圍為700～950，并且在794附件集中；圖3B中為蛋白質與小分子之間的接觸距離的頻率分布，從圖中我們可以看出，其值分布在0.1～0.3 nm之間，并且集中于0.22 nm，因此我們可以認為葡聚糖酶與二聚糖分子的結合較為緊密。

圖3 A酶與底物的接觸數，B酶與底物之間的距離

2.2.3 復合物結合自由能

通過上文對系統中RMSD的分析，我們選取了45～50 ns時間內的軌跡作為研究對象，每隔100 ps提取構象，共獲得51個構象，通過利用MM-PBSA方法計算結合自由能，其結果如表1所示，從表1我們可以看出，小分子UNK與葡聚糖酶的結合自由能為-54.718 kJ/mol，這個結果表明小分子UNK與葡聚糖酶產生了比較強的結合能力。范德華作用能(△Evdw)為-88.473 kJ/mol，該結果表明范德華作用可以有效的促進小分子與蛋白質的結合。靜電相互作用能(△Eele)和非極性溶劑自由能(△Gnonel)分別為-7.919 kJ/mol和-8.164 kJ/mol，這表明該作用為小分子的結合提供了較小的貢獻。極性溶劑自由能(△Gel)為49.838 kJ/mol，這說明，該相互作用極大的衰減了靜電相互作用，不利于小分子與蛋白質的結合。

表1 MM-PBSA計算所得的能量(kJ/mol)

2.2.4 活性位點氨基酸殘基的貢獻

為了進一步研究小分子UNK與葡聚糖酶分子間的相互作用模式，通過利用MM-PBSA程序計算每個氨基酸殘基對結合小分子的貢獻，試驗結果如圖4所示，Gly46、Leu63、Gln67、Tyr68、His73、Thr124、Gly125和Gln129對結合小分子的貢獻較大，其中His73和Lys126殘基可能作為親核劑，為催化部位的活性位點，Gly46、Leu63、Gln67、Tyr68、Thr124、Gly125和Gln129可能與小分子結合識別有關。從表2 我們也可以看出，His73為堿性氨基酸，帶有正電荷，與小分子之前形成的靜電相互作用和范德華作用力較強，結合能為-4.4088 kJ/mol，Tyr68為疏水性氨基酸，極性溶劑自由能對靜電作用影響較小，其對小分子的結合自由能為-4.6058 kJ/mol，其他氨基酸殘基對小分子的結合自由能貢獻均低于-1 kJ/mol。與上文對接試驗結果較為一致。

圖4 底物與酶結合的相互作用譜

表2 能量分解的主要貢獻殘基各部分能量值

2.2.5 氫鍵分析

葡聚糖酶與小分子UNK之間形成的氫鍵數由圖5所示，從圖5可知，氫鍵數目隨時間的變化而變化，在0～10 ns內，氫鍵數目呈遞增的趨勢，最大值為21；在20～40 ns時間段內，其氫鍵數目處于動態平衡，在3～15范圍內波動；其整體的平均值為7。通過利用VMD軟件分析，我們可以得出表3， 從表3我們可以看出His73和Thr124殘基與小分子之間的氫鍵占有率均大于50%，這說明這些氨基酸殘基與小分子之間形成的氫鍵較為穩定。

圖5 底物與酶之間氫鍵數隨時間的變化

表3 主要殘基的氫鍵分析

3 結論

本文利用分子對接與動力學模擬的方法，研究葡聚糖酶結合底物的活性位點。通過分子對接技術，選取最佳的分子構象用于分子動力學模擬，用MM-PBSA方法計算了小分子UNK和葡聚糖酶之間的結合能，計算結果顯示范德華相互作用、靜電相互作用和非極性溶劑化自由能是形成復合物的主要驅動力，通過對各個氨基酸殘基的貢獻研究發現，Gly46、Leu63、Gln67、Tyr68、His73、Thr124、Gly125和Gln129是對小分子UNK的結合起關鍵作用的氨基酸殘基。通過進一步研究小分子與葡聚糖酶之間的氫鍵數目及其對關鍵氨基酸的氫鍵分析，我們可以得出，小分子UNK與葡聚糖酶的結合較為緊密，同時也可以得出，His73和Thr124殘基與小分子所形成的氫鍵較為穩定，對結合能貢獻較大。