大豆中一种富含半胱氨酸蛋白性质的分子动力学模拟

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

(1. 吉林工商学院 工学院， 长春 130507； 2. 吉林大学 化学学院理论化学研究所， 长春 130021)

大豆是人类获取优质蛋白的主要来源， 但中国大豆产量较低. 因此， 发展具有优越性能的大豆分子设计育种技术非常重要[1]. 在大豆育种和生产过程中， 如何获得营养丰富、 功能特异的蛋白是重要目的之一. 以大豆中富含半胱氨酸的蛋白为例， 其主要用于医药、 化妆品和生化研究等方面， 可用于美容水、 烫发液、 防晒霜等产品中， 并可促进面包发酵、 防止维生素C氧化、 解丙烯腈中毒、 预防放射线损伤人体、 治疗支气管炎等[2-10]. 实验结果表明， 富含半胱氨酸蛋白还具有配合金属离子[11]等功能. 具有疏水性的富含半胱氨酸蛋白的晶体结构已被解析[12]， 但通过实验研究富含半胱氨酸蛋白与溶剂相互作用的文献较少[13-17].

利用计算机可模拟研究富含半胱氨酸蛋白与溶剂相互作用. 计算机模拟具有成本低、 速度快等优点[18-25]， He等[26]构建了植物类金属硫蛋白半胱氨酸富含区结构的模型， 但未进行详细的动力学描述. 本文利用分子动力学方法， 研究大豆中一种富含半胱氨酸蛋白的疏水性质以及蛋白与溶剂间的弱相互作用， 并给出其自由能变化规律.

1 研究方法及模型结构

分子动力学方法是研究化学小分子和生物大分子的有效方法， 在研究蛋白质的热力学和动力学性质等方面具有一定优势[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两个状态， 其中λ=0表示考虑了范德华相互作用，λ=1表示未考虑范德华相互作用，λ为耦合参数没有实际意义， 但通过∂H/∂λ曲线可计算其在水溶液中的自由能变化ΔG. 富含半胱氨酸蛋白在模拟中采用OPLS-AA力场， 水分子采用SPC(simple point charge)模型. 为更清晰地观察蛋白质， 未画水分子. 模拟在具有周期性边界条件、 尺寸为5.7 nm3的立方箱中进行， 全部计算使用GROMACS软件包完成[38-39]. 对每个λ状态均进行了预处理和足够长时间的平衡， 流程如下：

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组螺旋结构和扭曲的β链. 由于溶剂的环境不同， 蛋白与溶剂间的弱相互作用也不同， 因此蛋白的构象有细微差别.

图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时， 螺旋结构基本消失， 蛋白质无法保持折叠结构， 这是由于蛋白的疏水性和溶剂环境变差共同作用的结果.当λ=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状态到λ=1状态是一个自由能升高、 不易发生的过程. 其中λ=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力场的分子动力学方法研究了大豆中富含半胱氨酸蛋白的疏水性质以及蛋白与溶剂间的弱相互作用. 结果表明， 半胱氨酸对蛋白的性质具有决定性影响， 溶剂主要影响蛋白的构象. 并从热力学和动力学上研究了该蛋白在不同λ下的自由能变化规律, 从分子层面给出了蛋白质构象的变化细节， 为大豆分子设计育种提供了新思路.