高性能分子动力学模拟在生物大分子结构和功能研究中的应用

APPLICATIONe-Science应用60e-Science2010年11月高性能分子动力学模拟在生物大分子结构和功能研究中的应用基于经典物理学理论的分子动力学模拟是研究生物大分子结构和功能的重要工具。然而，分子动力学模拟对计算资源的巨大需求，限制了其在生物领域的应用。近几年来，随着计算机科学的迅猛发展以及计算化学理论和方法的完善，进行生物大分子在溶液状态下微秒时间尺度的模拟具备了可行性。本篇综述主要介绍近期高性能分子动力学模拟的发展，及在蛋白质三维结构、动态构象变化和生物学功能上的具体应用实例。分子力学；力场；分子动力学模拟；蛋白结构功能关系摘要：关键词：StudyingStructureandFunctionRelationshipofBiomoleculesviaHighPerformanceMolecularDynamicsSimulationsStateoftheartmoleculardynamicssimulationtechniquesapplythelawsofphysicstostudystructurefunctionrelationshipsofbiologicalmacromoleculesofinterest,whichenableustoacquireinformationnotobtainablebyexperimentaltechniquesandtoinvestigatestructuralphenomenaattheatomiclevel.Recentadvancesincomputerscienceandalgorithmdevelopmentmakethemicrosecondtimescalesimulationsontensofthousandsofatomspractical.Here,wehighlightthemostrecentprogressesinhighperformancemoleculardynamicssimulationsandbrieflyintroduceseveralapplicationprojectsinstudyingproteinstructure,dynamicsandbiologicalfunction.Molecularmechanics;Forcefield;Moleculardynamicssimulation;ProteinstructurefunctionrelationshipAbstract:Keywords:吴超，马路遥，黄牛*北京生命科学研究所，北京102206WuChao,MaLuyao,HuangNiu*NationalInstituteofBiologicalSciences,Beijing102206,ChinaAPPLICATIONe-Science应用e-Science2010年11月611.引言近年来，得益于计算资源的迅猛增长及生物大分子体系模拟理论和算法的日趋成熟，计算化学方法逐渐被广泛应用于生物学领域。其主要通过经典分子力学计算实时地模拟生物体系中各个原子的微观运动，结合统计力学原理研究体系的宏观物理化学性质，从而允许我们获得目前实验条件下难以观测到的微观动态情况[1,2]；预测可能进行的实验过程以及在原子水平上分析和解释实验结果[3,4,5,6,7,8,9]；设计和开发新型结构的活性化合物[10,11]。“生命在于运动”在分子水平上具有了全新的内涵。生物大分子功能的发挥与其三维结构在生物环境中的构象动态变化紧密相关[12,13,14,15]，RichardFeynman曾在其著名的费曼物理学讲义中表述过这一观点，“毫无疑问，没有哪个学科或领域能像生物学一样，可以同时在如此多的前沿有如此多的进展。在人们前赴后继地试图理解生命现象的过程中，如果要选出一个最基本的假设，那就是所有的物体都是由原子组成的，所有的生命现象都可以用这些原子的运动来解释。”[16]分子动力学模拟（MolecularDynamicsSimulation）方法以经典物理学为理论基础，在原子水平上研究大分子在三维空间中随时间的动态变化过程。1977年McCammon等人第一次利用全原子模型结合分子动力学模拟来研究蛋白质在真空状态下的动态构象变化[17]，预示着这一全新研究领域的诞生。三十多年来，随着模拟理论算法的改进和计算科学的发展，分子动力学模拟在生命现象的研究中得到了越来越多的应用[18]，其重要性开始被实验生物学家所认可。本篇综述通过具体例子的阐述，主要介绍高性能分子动力学模拟在蛋白质分子结构和功能领域中的应用。全文的组织结构如下：首先介绍分子动力学模拟的基本原理，以及其高性能计算中相关的软硬件发展；然后以具体例子的形式说明高性能分子动力学模拟在研究蛋白质结构和动态变化及其生物学功能中所起的重要作用；最后讨论高性能分子动力学模拟的发展方向。2.分子动力学模拟2.1分子动力学模拟的基本原理分子动力学模拟的核心是分子力场（Forcefield），包括势能量函数的形式及其相应的力场参数。力场势能量函数通常包括键长、键角、二面角等反映分子内成键原子间作用的共价键项，以及反映非成键原子间相互作用的范德华、静电等非共价键项[19]。一个分子系统的总势能V(R)等于所有原子的共价作用能量V(R)bonded和非共价作用能量V(R)non-bonded之和。V(R)=V(R)bonded+V(R)non−bonded(1)V(R)bonded=Kb(b−b0)2+(2)Kθ(θ−θ0)2+Kχ[1+cos(nχ−σ)]V(R)non−bonded=(3)(εij[(Rmin,ij/rij)12−(Rmin,ij/rij)6]+qiqj/εDrij)2其中b为键长，θ为键角，χ为二面角，rij为原子间的距离，如图1所示。在分子动力学模拟中，每个原子位置的变化遵守经典力学规律。Fi=miai(4)Fi=−DiV(5)公式（4）中Fi为作用在原子i上的力，mi为它的质量，ai为产生的加速度，其中Fi可以用公式（5）中的势函数梯度来计算。每隔一个模拟步幅，分子中的各个原子移动到新的位置，由此时的结构可以根据力场求出势能梯度，进而计算每个原子所受到的来自体系brijθ12345χ图1用于描述公式（1）-（3）中各个能量项的虚拟分子，包括共价连接的四个原子（1-4）及非共价作用的原子5。共价键能量项包括键长项b（1-2，2-3，3-4），键角项θ（1-2-3，2-3-4），二面角项χ（1-2-3-4）；非共价键能量项包括非共价原子对（1-5，2-5，3-5，4-5）之间的范德华项和静电项。APPLICATIONe-Science应用62e-Science2010年11月中其他原子的力，然后求出加速度，最后求出新的位置。如此循环下去，就可以得到描述体系各个原子位置、速度和加速度随时间变化的运动轨迹。根据各向异性假设，也就是系统在长时间内的取样平均等价于系综平均，通过统计力学方法计算，可以将分子动力学轨迹中包含的微观信息同宏观物理化学性质联系起来。这样，分子动力学模拟不但可以提供实验中很难测定的高精度动态微观结构信息，还可以通过理论计算得到宏观热力学量，从而可以用于预测实验测量数值。2.2分子动力学模拟的高性能计算在生物大分子的分子动力学模拟中，生物大分子连同其水溶液环境（或者膜环境）通常包括几万到上百万的原子数，在每步的模拟中，最耗计算机时间的势函数求解理论上与原子数的平方成正比，因此每步要做几亿到几万亿次计算。同时，分子中包含的原子的最快振动周期要求分子动力学模拟的步幅在1~2飞秒（10-15秒）范围内，而生物大分子行使功能所需的时间则多在微秒-秒范围内，甚至更长时间尺度，因此为得到有意义的相空间取样，需要至少109到1015动力学步骤。每步计算的规模之大加上需要计算的步骤之多，要求必须具有极高性能的超级计算机，才能真正意义上实现生物学有意义的时间尺度的模拟工作。以我们在中科院超算中心深腾7000超级计算机上进行的流感血凝素蛋白（HA）模拟为例，为研究该蛋白在介导膜融合过程中受低pH环境调控的构象变化，从而指导流感膜融合抑制剂的设计，我们构建了H3和H5亚型HA在不同pH环境中的数十个模拟体系（原子数在200，000～400，000），运用高度并行化的Desmond分子模拟软件[20]，进行了最高1024核的动力学模拟（128，256，512，1024核并行效率分别为95%，91%，75%，65%），模拟时间在100ns级别，总模拟计算量达到70CPU年，通过模拟揭示了在低pH调控的构象变化早期阶段，不同亚型HA蛋白构象变化的机制可能有明显区别，并鉴定了构象变化早期可能存在的能量亚稳态，而单独通过比较蛋白序列或晶体结构是无法获得这种基于蛋白动态构象变化得到的信息（未发表结果）。近年来，高性能分子动力学的模拟取得了长足的发展[21]。虽然在过去三年中，单机的计算速度只增长了50%，但模拟速度却几乎增长了一个数量级，其中很大一部分归功于分子动力学软件的并行化程度的提高，如采用Dynamicloadbalancing技术的NAMD模拟软件[22]，在300个核的计算集群上模拟一个30，000原子的体系速度达到了100纳秒/天[23]，而采用Neutralterritory并行技术的Desmond软件[20]，在1024核的通用CPU集群上模拟一个23，000原子的体系速度达到了471纳秒/天[24]。在计算硬件上也有一定的发展，以D.E.Shaw研究小组为例，他们设计了专门用于分子动力学模拟的机器Anton[25]。另外，用于图形加速的GPU的硬件的兴起，一些分子动力学软件如NAMD、AMBER在GPU上的重编程[26,27,28]，使得单芯片模拟速度有了几个数量级的提高，虽然节点间通讯的问题尚未解决，但这对于像Folding@Home这样的分配式计算模拟项目仍有较大的促进[26]。3.高性能分子动力学模拟在生物大分子研究中的应用3.1蛋白质折叠对蛋白质如何从一级序列折叠成具有生物学功能的特定三维结构问题（蛋白质折叠问题，见图2）的理解，是研究蛋白质序列和结构关系的关键步骤。对此机理的研究，有助于预测蛋白质结构，研究蛋白质复合化等基础生物学问题，促进药物分子的设计[29]。虽然不同的蛋白质序列不同，但折叠背后的物理化学原理却是一样的，因此基于物理学原理的理论方法对研究蛋白质折叠有着广泛适用性[30]。随着高性能计算的发展，基于物理学原理的蛋白质折叠将不再是一个遥不可及的梦想。1998年，UCSF的YongDuan和PeterKollman在Cray3E超级计算机上开展了第一个微秒级的蛋白质折叠分子动力学模拟[29]，有力APPLICATIONe-Science应用e-Science2010年11月63地推动了分子动力学的并行计算研究，堪称分子动力学历史上的一次壮举。在这项工作中，作者模拟了一个36个氨基酸的亚结构域VillinHeadpiece在水溶液中的折叠过程。计算结果显示，在该蛋白折叠的最初阶段，a天然螺旋含量逐渐上升直到60%，原子间的相互接触也从起始的3%逐渐升到45%。通过各种能量项的相关分析，推断这一阶段主要是疏水相互作用驱动的，可以视为初始疏水坍塌阶段。从开始到达到50%α螺旋含量大约需要50纳秒，同另外两个小蛋白的动力学实验数据较为符合。紧接着是一个较慢的调整阶段，天然α螺旋含量迅速下降到20%，然后和原子间接触一样又开始缓慢增长，这一时期为三级结构形成期。作者对所得到的构象进行了聚类分析，发现最大的一类构象同天然构象较类似，关键区域的RMSD为4.0Å，但侧链的构象并不接近天然构象。作者认为这并不奇怪，因为该蛋白的折叠时间下限在10微秒，而模拟时间仅为1微秒。另一方面，该类亚稳构象持续了150纳秒，时间相较天然构象折叠只小了2个数量级，又和天然构象较为接近，因此作者推断这类构象可能是折叠过程的其中一个中间状态，在折叠的过程还会