李 洋， 王科飛， 胡海燕， 王志超， 李延春

(1. 吉林工商學院 工學院， 長春 130507； 2. 吉林大學 化學學院理論化學研究所， 長春 130021)

大豆是人類獲取優(yōu)質蛋白的主要來源， 但中國大豆產量較低. 因此， 發(fā)展具有優(yōu)越性能的大豆分子設計育種技術非常重要[1]. 在大豆育種和生產過程中， 如何獲得營養(yǎng)豐富、 功能特異的蛋白是重要目的之一. 以大豆中富含半胱氨酸的蛋白為例， 其主要用于醫(yī)藥、 化妝品和生化研究等方面， 可用于美容水、 燙發(fā)液、 防曬霜等產品中， 并可促進面包發(fā)酵、 防止維生素C氧化、 解丙烯腈中毒、 預防放射線損傷人體、 治療支氣管炎等[2-10]. 實驗結果表明， 富含半胱氨酸蛋白還具有配合金屬離子[11]等功能. 具有疏水性的富含半胱氨酸蛋白的晶體結構已被解析[12]， 但通過實驗研究富含半胱氨酸蛋白與溶劑相互作用的文獻較少[13-17].

利用計算機可模擬研究富含半胱氨酸蛋白與溶劑相互作用. 計算機模擬具有成本低、 速度快等優(yōu)點[18-25]， He等[26]構建了植物類金屬硫蛋白半胱氨酸富含區(qū)結構的模型， 但未進行詳細的動力學描述. 本文利用分子動力學方法， 研究大豆中一種富含半胱氨酸蛋白的疏水性質以及蛋白與溶劑間的弱相互作用， 并給出其自由能變化規(guī)律.

1 研究方法及模型結構

分子動力學方法是研究化學小分子和生物大分子的有效方法， 在研究蛋白質的熱力學和動力學性質等方面具有一定優(yōu)勢[27-29]. 在分子動力學中， 體系中每個分子的運動均遵守牛頓方程

(1)

OPLS-AA(optimized potentials for liquid simulations all atom)力場是能準描述生物分子(如蛋白質)的力場， 其形式為

V(rN)=Vbonds+Vangles+Vdihedrals+Vnonbonded，

(2)

其中:

OPLS-AA力場包括分子間和分子內相互作用兩部分， 其中分子內相互作用包括化學鍵、 鍵角和二面角三部分； 分子間相互作用包括范德華和靜電相互作用.

圖1 自由能計算模型Fig.1 Free energy calculation models

文獻[12]研究了富含半胱氨酸蛋白的晶體結構(PDB： 1HYP)， 其分子由四組螺旋結構和扭曲的β鏈組成， 共含有75個氨基酸， 如圖1所示. 其中纈氨酸(Val)、 亮氨酸(Leu)、 異亮氨酸(Ile)是典型的疏水性氨基酸， 天冬氨酸(Asp)、 天冬酰胺(Asn)、 谷氨酸(Glu)、 谷氨酰胺(Gln)、 精氨酸(Arg)是典型的親水性氨基酸， 脯氨酸(Pro)、 絲氨酸(Ser)、 半胱氨酸(Cys)、 甘氨酸(Gly)、 酪氨酸(Tyr)、 丙氨酸(Ala)既不親水也不疏水. 分子中4個精氨酸帶正電荷， 4個天冬氨酸和1個谷氨酸帶負電荷， 因此整個蛋白分子帶1個負電荷. 通過加入1個帶正電的鈉離子使整個模擬體系保持電中性. 富含半胱氨酸蛋白中含有24個疏水性氨基酸， 17個親水性氨基酸和34個既不親水也不疏水的氨基酸， 整體呈疏水性. 氨基酸序列為： Pro Ser Cys Pro Asp Leu Ser Ile Cys Leu Asn Ile Leu Gly Gly Ser Leu Gly Thr Val Asp Asp Cys Cys Ala Leu Ile Gly Gly Leu Gly Asp Ile Glu Ala Ile Val Cys Leu Cys Ile Gln Leu Arg Ala Leu Gly Ile Leu Asn Leu Asn Arg Asn Leu Gln Leu Ile Leu Asn Ser Cys Gly Arg Ser Tyr Pro Ser Asn Ala Thr Cys Pro Arg Thr.

采用Bennett接受率(Bennett acceptance ratio, BAR)方法[30-37]計算蛋白的自由能差值. 由于分子間相互作用由兩部分組成， 一部分是靜電相互作用， 另一部分是范德華相互作用， 如圖1所示， 因此構建了該蛋白的λ=0和λ=1兩個狀態(tài)， 其中λ=0表示考慮了范德華相互作用，λ=1表示未考慮范德華相互作用，λ為耦合參數沒有實際意義， 但通過?H/?λ曲線可計算其在水溶液中的自由能變化ΔG. 富含半胱氨酸蛋白在模擬中采用OPLS-AA力場， 水分子采用SPC(simple point charge)模型. 為更清晰地觀察蛋白質， 未畫水分子. 模擬在具有周期性邊界條件、 尺寸為5.7 nm3的立方箱中進行， 全部計算使用GROMACS軟件包完成[38-39]. 對每個λ狀態(tài)均進行了預處理和足夠長時間的平衡， 流程如下：

1) 利用最陡下降法， 進行能量最小化；

2) 利用L-BFGS(limited-memoy Broyden-Fletcher-Goldfarb-Shanno)法進行能量最小化；

3) 采用正則(NVT)系綜平衡0.5 ns；

4) 采用等溫-等壓(NPT)系綜平衡0.5 ns；

5) 進行動力學模擬5 ns.

2 結果與討論

2.1 動力學過程

選取λ=0,0.05,0.1,…,1等21個數值， 按上述流程準備21個平衡的初始構象， 共進行21個5 ns的分子動力學模擬.圖2給出了λ=0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.75時的9個平衡構象.由圖2可見， 當λ=0～0.75時， 由于存在疏水性氨基酸， 蛋白質整體呈疏水性.因此， 整個蛋白質基本保持了4組螺旋結構和扭曲的β鏈. 由于溶劑的環(huán)境不同， 蛋白與溶劑間的弱相互作用也不同， 因此蛋白的構象有細微差別.

圖2 當λ=0～0.75時的平衡結構Fig.2 Equilibrium structures when λ=0—0.75

圖3給出了λ=0.8,0.85,0.9,1時的4個平衡構象.由圖3可見， 當λ=0.8～1時， 螺旋結構基本消失， 蛋白質無法保持折疊結構， 這是由于蛋白的疏水性和溶劑環(huán)境變差共同作用的結果.當λ=1時， 蛋白與溶劑間的非鍵相互作用為靜電相互作用， 范德華相互作用完全消失， 這種情況下蛋白質會失去活性， 也不利于實際應用.

圖3 當λ=0.8～1時的平衡結構Fig.3 Equilibrium structures whenλ=0.8—1

2.2 自由能變化

圖4為不同λ下的自由能變化曲線， 總的自由能變化ΔG=108.41 kT=270.41 kJ/mol， 表明從λ=0狀態(tài)到λ=1狀態(tài)是一個自由能升高、 不易發(fā)生的過程. 其中λ=0.75→λ=0.8時， ΔG=-3.96 kT=-9.87 kJ/mol；λ=0.8→λ=0.85時， ΔG=-32.20 kT=-80.32 kJ/mol；λ=0.85→λ=0.9時， ΔG=-64.19 kT=-160.10 kJ/mol；λ=0.9→λ=0.95時， ΔG=-28.03 kT=-69.92 kJ/mol. 只有這4個轉變過程的自由能為負值； 其他16個轉變過程的自由能均為正值. 這4個轉變過程對應蛋白質從螺旋結構到解螺旋結構的轉變.

圖5為不同λ下的累積自由能變化曲線. 由圖5可見： 當λ=0時， ΔG(λ=0)=0 kT=0 kJ/mol； 當λ=1時， ΔG(λ=1)=G(λ=0)+G(λ=0.05)+…+G(λ=0.95)=108.41 kT=270.41 kJ/mol. 因此， 從累積的自由能變化也可得到總的自由能變化ΔG=270.41 kJ/mol. 累積的自由能變化曲線從λ=0.8開始有一個下降過程， 直到λ=1時， 略有上升.

圖4 不同λ下的自由能變化曲線Fig.4 Free energy curves with different λ

圖5 不同λ下的累積自由能變化曲線Fig.5 Cumulative free energy curves with different λ

圖6為不同λ樣本的Hamilton ΔH變化曲線. 由圖6可見， ΔH主要分布在0～1 000間， 樣本基本分布在1 000～1 500間， 表明計算結果不受采樣影響. 圖7為兩個不同λ樣本的Hamilton ΔH變化曲線.N(ΔH(λ=0.05)|λ=0)表示從λnative=0的軌跡中算出λforeign=0.05的ΔH. 計算結果表明， ΔH主要分布在-100～100間， 樣本基本分布在1 000～2 000間， 同樣表明計算結果不受采樣影響.

圖6 ?H/?λ在每個λ數值下的分布Fig.6 Distributions of ?H/?λfor each λ value

圖7 ?H/?λ在兩個 λ數值下的分布Fig.7 Distributions of ?H/?λfortwoλvalues

綜上所述， 本文利用基于OPLS-AA力場的分子動力學方法研究了大豆中富含半胱氨酸蛋白的疏水性質以及蛋白與溶劑間的弱相互作用. 結果表明， 半胱氨酸對蛋白的性質具有決定性影響， 溶劑主要影響蛋白的構象. 并從熱力學和動力學上研究了該蛋白在不同λ下的自由能變化規(guī)律, 從分子層面給出了蛋白質構象的變化細節(jié)， 為大豆分子設計育種提供了新思路.