String方法在丙氨酸多肽链中的应用

贾红霞，靳聪明

(浙江理工大学理学院 杭州 310018)



String方法在丙氨酸多肽链中的应用

贾红霞，靳聪明

(浙江理工大学理学院 杭州 310018)

蛋白质的构象变化在生命过程中起着重要的作用，蛋白质在发挥其功能时常伴有构象变化，有些疾病也是由蛋白质的构象发生变化引起的，典型的例子如老年痴呆、疯牛病等就是某些蛋白质的构象由水溶性的蛋白变成了不溶性的纤维性蛋白引起的。知道构象变化的过程对理解蛋白质相关性质起着重要作用，但从原子尺度研究构象变化是极具挑战的问题。文章利用String方法从原子尺度研究了多肽链的310-螺旋和α-螺旋之间的构象变化。文中用零温度String方法研究多肽链310-螺旋和α-螺旋之间的构象变化，从不同的初始路径得到了不同的最优路径和过渡态，其中有显著势垒的过渡态是与310-螺旋相似的构型；利用有限温度String方法研究多肽的构象变化，采用了Langevin动力学和分子动力学两种采样法，得到了它的转移路径系综，通过与零温度String方法结果比较说明了模拟结果的有效性。

String方法；多肽链；转移路径；构象变化

0　引　言

蛋白质构象变化决定了蛋白质的性质，蛋白质在发挥其功能时常伴有构象变化。研究蛋白质构象变化对研究蛋白质的性质和功能非常重要。蛋白质的配体与受体之间对接时常伴有分子的构象变化，知道了蛋白质的构象变化过程即可设计出更好的配体来与受体作用。再如有些疾病是由蛋白质的构象发生变化引起的，如老年痴呆、疯牛病等就是Prion蛋白质[1-2]的构象由螺旋结构占优的水溶性的蛋白变成了β-折叠占优的不溶性纤维性蛋白引起的。弄清楚这个变化过程，可以指导药物设计来阻止病变甚至消除引起病变的蛋白质。再如血清白蛋白是人和动物内血浆中含量最丰富的蛋白质，在生理上有着十分重要的载运功能，可作为金属离子、脂肪酸、氨基酸、代谢物、胆红素、酶、药物和激素的载运体，研究血清白蛋白的构象变化是探讨其生物功能奥秘的重要途径。

蛋白质构象变化就是蛋白质从一个状态变化到另一个状态，蛋白质在行使其生理功能时往往会发生构象变化。蛋白质构象变化发生的时间尺度大多在10-6～10-3s，运动幅度一般达到5～10Å(Å =10-10m)。而原子的振动时间尺度为10-15s。在蛋白质构象变化过程中，原子尺度的振动和大幅度的结构变化耦合在一起，使构象变化的过程非常复杂，是一个多尺度问题。计算机模拟在蛋白质构象变化研究中起着重要的作用，可以处理各种条件下的构象变化，不受实验条件的限制。而且构象变化过程中的关键状态—过渡态—很难在实验中观察到，因为出现的概率很低，且寿命很短，而计算机模拟可以抓住这些状态。但目前传统的计算模拟方法主要是分子动力学方法，能模拟的时间尺度为10-12～10-9s，因此很难处理蛋白质构象变化。针对这个问题人们发展了一些新技术。例如靶向分子动力学[3-4]、路径采样法(transition path sample)[5]、伞形抽样方法[6]、NEB(nudged elastic band)方法[7]等，这些方法往往计算量很大。

由于蛋白质构象变化相对原子的振动时间长得多，因此在原子尺度考虑蛋白质构象变化时，蛋白质构象变化可以看作小概率事件。由鄂维南教授等提出的String方法[8-10]是处理小概率事件的有效方法，被广泛用于位错运动[11-12]、相变中的成核[13]、固体融化机制[14]、生物分子构象变化[15]等的研究。α-螺旋与310-螺旋之间的变化是蛋白质中重要的构象变化，与很多生物过程有密切的联系，比如一些配体与受体的结合中有α-螺旋与310-螺旋之间的转换，它是一些酶发生作用时的重要过程[16-17]。丙氨酸多肽链是生物分子模拟中常用的代表性生物分子模型[18-20]。本文利用String方法研究了丙氨酸多肽链的310-螺旋和α-螺旋之间的构象变化。在计算中采用了AMBER(assisted model building and energy refinement)力场参数，由于计算量很大，在运算中我们采用了并行计算。

String方法包括零温度String(zero temperature string,ZTS)方法[9]和有限温度String(finite temperature string，FTS)方法[8]。本研究用ZTS方法研究了多肽链310-螺旋和α-螺旋之间的构象变化，从不同的初始路径得到了不同的最优转移路径以及多个过渡态。接着用FTS方法研究多肽的构象变化，并且得到了它的转移路径系综，并与ZTS方法的结果进行了比较。

1　String方法

1.1零温度String方法

假设A和B是系统的两个稳态。令φ是连接A和B的一条光滑曲线，称之为String。φ有参数表示形式φ(α)，α是它的内部参数。φ是从A到B的最小能量路径(minimalenergypath，MEP)，如果φ满足：

(▽V)⊥(φ)=0

(1)

(2)

其中：γ≡γ(α,t)是一个Lagrange乘子，由参数α的选取来决定，最简单的可以选α为φ的归一化弧长，在A有α=0，在B有α=1。这时φ满足下面的约束条件：

(|φα|)α=0

(3)

γ由此约束条件给出。其它的参数选取可以通过改变式(3)得到。假设已知V(X)的两个稳态A和B，下面是寻找MEP的零温度String方法的算法。

b)解微分方程。由第n次的迭代结果，利用下面的差分方程求第n+1次的路径：

(4)

c)重新参数化String。利用插值使路径上的点均匀分布。

d)转向b)。

1.2有限温度String方法

零温度String方法只能处理光滑势能面问题，如何计算粗糙势能面系统的最优转移路径一直是人们关心的问题，例如蛋白质的构象变化和蛋白质折叠等。目前处理这类问题的方法主要有转移路径采样法[21]、有限温度String方法[8]、NEB方法[7]等。

在粗糙势能面中，转移路径不是简单的一条曲线，而是多条转移概率差不多的路径组成的集合，称为路径系综，有限温度String方法可以找到路径系综的平均值φ(α)，也就是有限温度String方法中的String。下面是有限温度String方法算法。

a)给出初始String；

e)重新初始化动力学，转到b)。

2　String方法在丙氨酸多肽链中的应用

α-螺旋是蛋白质中常见的二级结构，其中第i个氨基酸的C=O中的O原子与第i+4个氨基酸的N-H中的H原子形成氢键。二面角{φ,φ}在{-60°,-50°}附近。310-螺旋中第i个氨基酸的C=O中的O原子与第i+3个氨基酸的N-H中的H原子形成氢键，在氢键间有10个原子，因此得名，二面角{φ,φ}在{-60°,-70°}附近。

本文把String方法应用到丙氨酸多肽链ACE-(ALA)20-NME的α-螺旋和310-螺旋间的变化，该多肽链由20个丙氨酸组成，两头的分子基团是ACE(CH3CO-)和NME(-NHCH3)。ACE一端被称为C端，NME一端被称为N端。

310-螺旋比α-螺旋长，在从α-螺旋到310-螺旋的变化过程中，它可能从一头开始伸长，从C端或者从N端，也可能整条链均匀伸长。从这3种机制可以给出3条初始路径，对α-螺旋从C端每次改变两个氨基酸的二面角φ和φ使其变化成310-螺旋的二面角的值，直到整条肽链都变成310-螺旋的二面角的值，这样得到一条初始路径。同理，从N端进行这个过程，得到第2条初始路径。同时改变所有的二面角使得α-螺旋逐渐变为310-螺旋，得到第3条初始路径。

2.1零温度String方法对3条路径的计算

首先分析零温度String方法对这3种初始路径的收敛情况，定义：

a)均匀拉伸的初始路径

当String上有26个点时，Fmax在10-3附近震荡，不能进一步收敛到0；当String上有102个点时，Fmax<10-5可以达到。

b)从C端拉伸的路径

路径上有38个点和102个点时，Fmax<10-5可以达到。

c)从N端拉伸的路径

当路径上有38个点和102个点时，Fmax都只能在10-3附近震荡。

得到的MEP的性质在图1中给出，可以看出能量相对光滑，每条路径上都存在一些中间稳态，即局部极小值。

另外如果计算得到的转移路径是MEP，它应该满足方程：

(6)

其中，M为质量矩阵。即当用如下方程

(7)

优化MEP上的系统时，系统应该沿着MEP达到邻近的极小势能状态。通过分析，当优化零温度String方法的数值解MEP上的系统时，系统确实沿着得到的MEP到达邻近的极小能量状态。这进一步说明得到的最优转移路径确实是MEP。

记C38为从C端拉伸的38个点的MEP，C101为从C端拉伸的101个点的MEP，U102为均匀拉伸的有102个点的MEP。中间稳态的信息如表1所示。

表1　MEP上的稳态能量

注：符号∃表示MEP上有这个中间稳态，但优化MEP上的系统时没有系统收敛到该状态；符号√表示MEP上有这个状态，同时优化MEP上的系统时有系统收敛到该状态；符号—表示无此状态。

从不同的路径我们求得了不同的过渡态和中间稳态，而且路径收敛情况比较复杂。

表2中列出了各条MEP上的过渡态以及能量值。

表2　MEP上过渡态的能量

注：括号中的数字表示MEP上过渡态的编号，符号—表示无更多状态。

图1　丙氨酸链能量氢键示意图注：○表示对应的氨基酸的氧原子生成的是α-类型的氢键，□表示对应的氨基酸的氧原子没有生成氢键，◇表示对应的氨基酸的氧原子生成的是310-类型的氢键。

2.2有限温度String方法对转移路径的计算

对从C端拉伸、上面有38个点的初始路径，在有限温度String方法中用了两种采样方法，Langevin动力学和分子动力学，分别记为FTS-LD和FTS-MD。模拟的系统是正则系综，温度为272 ℃。目前为了保证结果的稳定性，超平面内的采样都是限制在一个管子里面，即在约束条件下进行的。在FTS-LD方法中，约束管内侧半径为1.2 Å，外侧半径为3.0 Å。在FTS-MD方法中，约束管内侧半径为0.8 Å，外侧半径为2.0 Å。在计算中进行了100步迭代，每步迭代中在垂直于路径的超平面中采样进行了1000步动力学，积分步长取为10-15s。

在重新参数化String时，用了String上相邻两个多肽的主链上的原子间距离作为参数。所谓主链上的原子就是形成二面角φ和φ的那些碳原子。

(8)

i=0,1,2,…,N

(9)

另外分析了由FTS-MD方法和FTS-LD方法得到的最优转移路径氢键的类型，可以看到α-类型的氢键基本上是从C端开始断裂并接着形成310-类型的氢键，与ZTS方法的结果所不同的是，氢键在两种类型氢键间可能存在反复变化，尤其在多肽链的两头，这正是动力学的特点。

图2　丙氨酸链的均方根偏差

3　讨　论

本文用零温度String方法和有限温度String方法研究多肽链310-螺旋和α-螺旋之间的构象变化，从不同的初始路径得到了不同的最优路径和最优路径上的多个稳态和过渡态，给出了最优转移路径系综。结果显示不同转移路径的势垒相近，因此3种转移机制存在竞争关系。从过渡态出发进行能量优化，可以发现系统会沿着所得到的最优转移路径到达邻近的稳态或另一个结构非常相似的稳态，从而说明了得到的最优转移路径的合理性。通过比较零温度String方法和有限温度String方法得到的转移路径的均方根偏差和氢键模式变化情况，可以说明两种方法得到的最优转移路径相似，同时也说明了有限温度String方法结果的合理性。今后通过把String方法应用到实际的大分子中，如蛋白质或DNA的构象变化，并考虑各种系综下的构象变化，希望对实际问题(例如药物设计等)的研究有所帮助。

[1]HYAREH,YOUSRYT.Humanpriondiseases[J].HumanBrainMapping, 2015, 3: 683-691.

[2]GOOLDR,MCKINNONC,TABRIZISJ.Priondegradationpathways:potentialfortherapeuticintervention[J].MolecularandCellularNeuroscience, 2015, 66(A): 12-20.[3]FRIEDRICHSMS,EASTMANP,VAIDYANATHANV.Acceleratingmoleculardynamicsimulationongraphicsprocessingunits[J].JournalofComputationalChemistry, 2009, 30(6): 864-872.

[4]SCHLITTERJ,ENGELSM,KRÜGERP.Targetedmoleculardynamics:anewapproachforsearchingpathwaysofconformationaltransitions[J].JournalofMolecularGraphics, 1994, 12(2): 84-89.

[5]BOLHUISPG,DELLAGOC,CHANDLERD.Reactioncoordinatesofbiomolecularisomerization[J].ProceedingsoftheNationalAcadenyofSciences, 2000, 97(11): 5877-5882.

[6]MAN,HUAC,VANDERVAARTA.Freeenergysimulationofhelicaltransitions[J].JournalofComputationalChemistry, 2013, 34(8): 640-645.

[7]MATHEWSDH,CASEDA.NudgedelasticbandcalculationofminimalenergypathsfortheconformationalchangeofaGGnon-canonicalpair[J].JournalofMolecularBiology, 2006, 357(5): 1683-1693.

[8]EWN,RENWQ,VANDEN-EIJNDENE.Finite-temperaturestringmethodforthestudyofrareevents[J].TheJournalofPhysicalChemistryB, 2005, 109(14): 6688-6693.

[9]EWN,RENWQ,VANDEN-EIJNDENE.Stringmethodforthestudyofrareevents[J].PhysicalReviewBCondensedMatter, 2002, 66(5): 052301.

[10]RENWQ,VANDEN-EIJNDENE,MARAGAKISP,etal.Transitionpathwaysincomplexsystems:applicationofthefinite-temperaturestringmethodtothealaninedipeptide[J].JournalofComputationalChemistry, 2005,123(13): 134109.

[11]JINCM,RENWQ,XIANGY.Computingtransitionratesofthermallyactivatedeventsindislocationdynamics[J].ScriptaMaterialia, 2010, 62(4): 206-209.

[12]JINCM,XIANGY,LUG.Cross-slipofthescrewdislocationinaluminium[J].PhilosophicalMagazine, 2011, 91(32): 4109-4125.

[13]LITJ,ZHANGPW,ZHANGW.Nucleationratecalculationforthephasetransitionofdiblockcopolymersunderstochasticcahn-hilliarddynamics[J].MultiscaleModeling&Simulation, 2013, 11(1): 385-409.

[14]SAMANTAA,TUCKERMANME,YUT,etal.Microscopicmechanismsofequilibriummeltingofasolid[J].Science, 2014, 346(6210): 729-732.

[15]OVCHINNIKOVV,KARPLUSM.Investigationsofalpha-helix:beta-sheettransitionpathwaysinaminiproteinusingthefinite-temperaturestringmethod[J].JournalofChemicalPhysics, 2014, 140(17): 175103.

[16]HUSTONSE,MARSHALLlGR.Alpha/3(10)-helixtransitionsinalpha-methylalaninehomopeptides:conformationaltransitionpathwayandpotentialofmeanforce[J].Biopolymers, 1994, 34(1): 75-90.

[17]NELLASRB,JOHNSONQR,ShenT.Solvent-inducedα-to3(10)-helixtransitionofanamphiphilicpeptide[J].JournalofInorganicBiochemistry, 2013, 52 (40): 7137-7144.

[18]BIELERNS,HUNENBERGERPH.Ontheambiguityofconformationalstates:AB&S-LEUSsimulationstudyofthehelicalconformationsofdecaalanineinwater[J].JournalofChemistryPhysical, 2015, 142(16): 165102.

[19]HAZELA,CHIPOTC,GUMBARTJC.Thermodynamicsofdeca-alaninefoldinginwater[J].JournalofChemicalTheoryandComputation, 2014, 10(7): 2836-2844.

[20]KOKUBOH,HUCH,PETTITTBM.Peptideconformationalpreferencesinosmolytesolutions:transferfreeenergiesofdecaalanine[J].JournaloftheAmericanChemicalSociety, 2011,133(6): 1849-1858.

[21]BOLHUISPG,CHANDLERD,DellagoC.Transitionpathsampling:throwingropesoverroughmountainpasses,inthedark[J].AnnualReviewofPhysicalChemistry, 2002, 53(1): 291-318.

(责任编辑： 许惠儿)

The Application of String Method in Alanine Polyalanine Chain

JIAHongxia，JINCongming

(Department of Mathematics, Zhejiang Sci-Tech University, Hangzhou 310018,China)

Conformational changes of proteins play an important role in life processes. The protein plays its role often accompanied by a conformational change. Some diseases are caused by conformational change of protein, such as Alzheimer’s disease andBSE (bovine spongiform encephalopathy ), where a water-soluble protein becomes an insoluble protein fiber form. Knowing the process of conformational change of protein is very important to understand the properties of protein. However, finding the transition path of the protein conformational change from atomic scale is a big challenge. In this paper, we applied the zero-temperature string methods in protein conformational changes of polypeptide chain between 310-helix and α-helix. From different initial paths, we got different optimal paths and different transition states; the transition state with significant barriers is the configuration similar to the 310-helix. And we also applied the finite temperature string method to study the conformational change and got the transition path ensemble by using two sampling methods, namely Langevin dynamics and molecular dynamics. The result indicates the effectiveness of the finite tem perature string method when compared with the zero temperature string method.

String method; polypeptide chain; transition path;conformational changes

10.3969/j.issn.1673-3851.2016.03.023

2015-07-03

浙江省自然科学基金项目(LY13A010014)；国家自然科学基金项目(11571314)； 浙江理工大学521人才培养计划

贾红霞(1991-)，女，山西晋中人，硕士研究生，主要从事蛋白质构象变化方面的研究。

靳聪明，E-mail：jincm@lsec.cc.ac.cn

O242.1

A

1673- 3851 (2016) 02- 0297- 07 引用页码： 030705