CN105787292B - 蛋白质折叠的并行预测方法 - Google Patents

Abstract

本发明涉及蛋白质折叠的并行预测方法，包括以下步骤：不同的蛋白质结构通过分子动力学模拟得到；设定不同温度，通过副本互换得到更多不同构象的蛋白质结构；采用共轭帽分子分割法对不同构象的蛋白质结构进行分割，将每个氨基酸从蛋白质中分解出来，并且给每个氨基酸分子碎片带上共轭帽；通过量子力学软件计算分子碎片的静电势，从而拟合出每个氨基酸分子碎片中各原子的电荷，并更新分子力场中的原子电荷。本发明为了研究蛋白质折叠问题，将上述考虑蛋白质极化效应的方案和基于副本互换的方法有机结合起来，从而克服动力学模拟中的原子相互作用势和取样技术这两大难题，以期能在蛋白质折叠模拟中获得比较好的效果。

Description

蛋白质折叠的并行预测方法

技术领域

本发明涉及的是一种蛋白质折叠预测的并行方法，属于生物信息技术、计算方法与计算机虚拟现实技术。

背景技术

蛋白质的结构决定功能，仅仅知道基因组序列并不能使我们充分了蛋白质的功能，更无法知道它是如何工作的。X射线衍射法(X-ray diffraction)和核磁共振技术(nuclear magnetic resonance)是目前获得蛋白质三维空间结构最有效的实验方法，虽然目前蛋白质晶体数据库中已经出现了非常多的晶体结构，但是由于实验的难度，很多重要的蛋白质目前并没有晶体结构的。而且更为关键的是蛋白质的功能发挥，不仅和其静态结构相关，而且和其结构的动态变化密切相关，以往的单纯的结构决定功能的研究已经转变成蛋白质结构决定蛋白质的动态行为，从而决定蛋白质的功能发挥的研究思路。

而蛋白质如何从合成出来的无规则结构，转变成具有完整的结构并发挥生物功能的蛋白质折叠过程就成为研究蛋白质科学的重要问题。蛋白质折叠问题被列为“21世纪的生物物理学”的重要课题。从生物信息学的统计方法，到基于物理的方法，人们提出了非常多的方案，而目前，比较行之有效的方法就是利用基于分子力场的动力学方法，此种方法已经在比较小的蛋白质研究中获得了比较大的成功。

目前蒙特卡洛模拟(MC)分子动力学(MD)是考察蛋白质动态行为及折叠的模拟方法中最重要的手段。蒙特卡洛模拟是建立在统计力学的基础上，使用真实的分子模型，用真实分子的键长、键角，根据实验的各种外部、内部条件，以及化学反应、物质变化的各种物理化学定律，来考察计算模型体系的各种统计性质的变化及其对所研究的问题给出统计参数。如今，蒙特卡洛模拟已经被广泛用于生命体系的研究，并取得了一系列进展。分子动力学模拟则是在牛顿力学基础上描述模拟体系运动随时间演化的方法。近年来，计算机技术和算法的发展使得分子动力学模拟所能研究体系的空间和时间尺度都大大增加。如今，分子动力学模拟已经能在较大的体系规模和较长的时间尺度上模拟蛋白质的动力学行为，模拟体系的规模也从最初BPTI的500多个原子发展到今天的百万个原子以上，时间尺度也从最初的皮秒发展到现在的微秒级。决定分子动力学模拟准确性和高效性主要有原子分子之间相互作用势、溶剂模型和取样技术这三个方面。

目前广泛应用于生物大分子模拟的分子力场(如CHARMM,AMBER,OPLS和GROMOS等力场)，然而由于它们都采用固定的原子点电荷计算体系的静电势能，而这种描述方式有一个很大的弊端，即固定点电荷的使用意味着这种模型是不能对于外界环境做出反应，而调整电荷分布。这样的力场就忽略了外势所引起的电荷转移和静电极化所带来的影响，而这一影响对于模拟体系的结构和能量性质起着非常重要的作用。因此，极化力场的研究和发展成为当今分子力场研究的热点之一。

目前引入极化作用主要通过三种模型，诱导偶极模型，浮动电荷模型和经典德鲁特振子(Drude oscillator)模型。在诱导偶极模型中每个原子具有一定的可极化率，外场影响是通过诱导偶极来实现的。在浮动点电荷模型中，这一电荷可以在外界环境的影响下随时调整分布，通过原子的电负性来作为收敛的判据。经典鲁特振子模型，即引入一个虚质量的电荷点，此电荷点通过弹簧与原子中心相连。而对外界环境的响应可以通过调整这一电荷点到原子中心的距离，故可以考虑外界环境对此原子的影响。国内Yang等发展的原子-键电负性均衡方法即属于此类方法，可以快速准确计算体系随外势变化的电荷分布，并将其应用到了有机分子和多肽的研究中去。

不管是分子动力学模拟还是蒙特卡洛模拟，得出有意义的物理量都是建立在对体系构象取样的充分性基础上。目前主要有两种解决途径，一是增加模拟时间，随着计算能力的提高和算法的改进，目前的动力学模拟已经能够达到微秒数量级。但这一数量级对于许多重要的蛋白质发挥功能的行为来说，还是太短，而且当所研究的体系处于一个能垒较大的局域极小值点时，通过常规的延长分子动力学模拟时间或者蒙特卡洛步数的作用是不大的。另一种方法是改进取样方法，可以帮助在较短的模拟时间内观察到需要长时间模拟才能观察到的现象。为了解决普通的模拟方法不能有效跨过能垒的问题，近年来很多基于广义系统统计力学的取样算法被建立和发展，如副本互换(replica exchange method)，其中包括温度副本互换(temperature replica exchange method TREM)和哈密顿副本互换(Hamiltonian replica exchange,HREM),其基本思想是利用处于高温或者改变体系势能的方式得到高能状态的构象，然后通过构象交换机制，将限于局部极小值的构象交换出来，从而达到取样的充分性。目前这一方法已经被广泛用于蛋白质折叠的研究，NMR所得蛋白质结构的修正，蛋白结合机制等的研究中，并取得了很大的进展。

除了广义系综外还有一类另外一种基于改变作用势的方法。此种加强局部取样的方法一般是需要预先定义一个反应坐标，然后在这个反应坐标上给定一个势，从而加强这个方向上的取样。如soft-core势的应用，通过改变势能函数的形式，从而使能垒降低，从而达到在模拟中容易跨过能垒的，伞形取样就是这样一类方法。

在蛋白质折叠过程中经常使用的是副本互换分子动力学模拟(REMD)的取样技术，目前在常规的小肽的折叠过程中已经取得了非常好的效果，可是目前由于计算资源的限制，蛋白质本身力场的限制，目前能够获得的折叠是非常有限的。另外，目前这种温度副本互换的方法，也存在如下问题：在实际模拟过程中，蛋白体系是处在水溶液体系中，而要取得较好的取样效果，由于所需的副本数目和体系总原子数目的平方根成正比，这样由于水分子数目非常巨大，模拟所需要的副本数目就基本取决于水分子的数目，所以所需的计算量也大大增大。实际上我们关注的是蛋白分子的构象空间，基于此，我们提出了一种新的理论计算模型，通过此计算方法，理论上，所需副本数目就只和蛋白的原子数目相关，而与水分子数目无关，从而使所需的副本数目大大减少，而且我们关心的蛋白质的构象取样效果和普通的温度副本互换没什么区别。

发明内容

针对上述技术不足，本发明的目的提供一种考虑蛋白质极化效应和副本互换的方法有机结合起来的蛋白质折叠预测的并行方法。

本发明解决其技术问题所采用的技术方案是：蛋白质折叠的并行预测方法，包括以下步骤：

1)通过分子动力学模拟得到蛋白质结构；

2)设定不同温度，通过副本互换得到在不同温度下的蛋白质不同构象；

3)采用共轭帽分子分割法对蛋白质不同构象进行分割，将每个氨基酸从蛋白质中分解出来，并且每个氨基酸的分子碎片带有共轭帽；

4)通过量子力学软件计算分子碎片的静电势，通过拟合得到每个氨基酸分子碎片的原子电荷；通过更新分子力场中相应的原子电荷后返回至步骤1)，直到达到设定时间为止。

所述设定不同温度，通过副本互换得到在不同温度下的蛋白质不同构象包括以下步骤：

对每个分子动力学模拟设定N个不同的温度，也就是N个副本；

在分子动力学模拟过程中，在N个副本中任意相邻的两个副本之间进行通讯，判断是否需要交换温度；如需要交换温度，则两个副本交换温度；否则两个副本不交换温度；得到交换温度后不同温度下的蛋白质不同构象。

所述判断是否需要交换温度包括以下步骤：

在两个具有不同温度的分子动力学模拟中，分别得到在两个副本下分子动力学模拟后具有不同构象的蛋白质的能量，使用波尔兹曼权重因子将二者作差；然后随机产生一个阈值，如果该差值大于阈值，则需要交换温度，否则不需要交换温度。

所述采用共轭帽分子分割法对不同构象的蛋白质进行分割包括以下步骤：

将不同温度下的蛋白质均在肽键处切开，在每个切开的地方引进一对共轭的帽子，该共轭帽子为与之相邻的两个或四个氨基酸残基，或者N端使用CH₃CO-、C端使用NCH₃的共轭帽子结构

本发明具有以下有益效果及优点：

1.本发明为了更好的研究蛋白质折叠动态过程，蛋白质的极化作用不可忽略，采用基于分块思想的共轭帽分子分割法(MFCC)结合量子力学方法来计算可极化电荷,考虑蛋白质周围环境的电性影响，从而考察到蛋白的极化效应，使目前的常规动力学模拟充分考虑到蛋白质的动态变化所引起的极化效应，从而可以得到更加准确的结果。

2.本发明为了研究蛋白质折叠问题，将上述考虑蛋白质极化效应的方案和基于副本互换的方法有机结合起来，从而克服动力学模拟中的原子相互作用势和取样技术这两大难题，以期能在蛋白质折叠模拟中获得比较好的效果。

3、拟采用量子力学的方法计算蛋白质的电荷性质，并且为了实时反应出蛋白的周围环境对结构的影响，在模拟过程中实时考虑蛋白质周围的环境，不断更新蛋白质的电荷性质，从而考虑蛋白质的极化效应。使得能量计算更加逼近真实情况。

4、针对蛋白质构象取样问题，为了降低常规副本互换所需副本数目，提高效率，采用基于广义系综的高效蛋白质构象取样方法。并将这种极化电荷和新型采样技术有机地结合到动力学模拟中去，最终构成一个高准确度、高效率的分子模拟算法，为研究蛋白质折叠提供理论计算手段。

附图说明

图1是本发明的方法流程图；

图2是分子碎片分割法示意图；

图3是实现副本互换与分子分割法方法示意图；

图4是AMNW软件在超算中心部署结构图；

图5是AMNW蛋白质折叠软件平台在中科院北京深腾超算中心并行效率测试结果图；

图6a是AMNW蛋白质折叠软件平台在天河II号上并行效率的测试结果图；

图6b是AMNW蛋白质折叠软件平台在天河II号上加速比的测试结果图。

具体实施方式

下面结合实施例对本发明做进一步的详细说明。

本发明涉及研究生命体系中的蛋白质折叠的预测方法，包括以下内容：针对极化效应这一问题，拟采用量子力学的方法计算蛋白质的电荷性质，并且为了实时反应出蛋白的周围环境对结构的影响，在模拟过程中实时考虑蛋白质周围的环境，不断更新蛋白质的电荷性质，从而考虑蛋白质的极化效应。针对蛋白质构象取样问题，为了降低常规副本互换所需副本数目，提高效率，拟采用最近发展的基于广义系综的高效蛋白质构象取样方法。并将这种极化电荷和新型采样技术有机地结合到动力学模拟中去，最终构成一个高准确度、高效率的分子模拟算法，为研究蛋白质折叠提供理论计算手段。

蛋白质折叠预测的并行方法，包括以下步骤：

针对极化效应这一问题，采用量子力学的方法计算蛋白质的电荷性质；在目前的方案中使用MFCC的方式来计算可极化的蛋白质的原子电荷。首先采用MFCC的方式将蛋白质划分成N个片段，然后利用量化方法计算分别计算蛋白质所有残基的量化性质，利用RESP的方法拟合出此时的电荷。

为了实时反应出蛋白的周围环境对结构的影响，在模拟过程中实时考虑蛋白质周围的环境，不断更新蛋白质的电荷性质，从而考虑蛋白质的极化效应包括以下步骤：

在计算完体系每个原子的电荷后，利用软件Delphi来计算这些蛋白质所产生的溶剂化电荷，然后再将这些溶剂化电荷代入下一轮的量化计算中，经过迭代自洽后我们就可以得到这种量化计算的考虑了蛋白质具体环境的可极化电荷了。

针对蛋白质构象取样问题，为了降低常规副本互换所需副本数目，提高效率，采用最近发展的基于广义系综的高效蛋白质构象取样方法。并将这种极化电荷和新型采样技术有机地结合到动力学模拟中去，最终构成一个高准确度、

高效率的分子模拟算法，为研究蛋白质折叠提供理论计算手段。包括以下步骤：对于一个新的蛋白体系，首先为了确定所选的温度范围和副本个数，对体系进行一系列不同温度下的全原子动力学模拟，以确定蛋白质所需跨过的能垒，从而确定每个副本的温度。

针对目前不能直接使用量化方法(如密度泛函方法或解波函数方法)对整个蛋白质进行计算，我们采用MFCC的分块方法来进行量化计算更新电荷。即针对给定结构的具有N个氨基酸的蛋白质可以用如下来表示P＝NA1-A2-A3-...-A_NC,(其中Ai(i＝1,2,3...N)是每个氨基酸单元，N和C分别是蛋白质的氮端和碳端)。为了要使用量子力学的方法计算蛋白质的电子结构等性质，采用分块的方法，即在肽键处切开，在每个切开的地方，我们引进一对共轭的帽子，分别为Ccap1和Ccap2表示。共轭的帽子的引进基于两个方法，一是为了保持原来切开前的共价键的性质，和尽量模拟原来被切掉的片段的电子结构的性质，二是直接使用传统的N端和C端小共轭帽(如在N端使用CH₃CO-，在C端使用NCH₃)。这样蛋白质的整体性质我们可以通过将其分解成N个片段的方法来进行，同时为了保证所计算性质的准确，我们根据自洽理论，将进行多次迭代循环，直至电荷的变化小与0.001,认为收敛完成。本实施例的共轭帽子可为与之相邻的两个或四个氨基酸残基，也可使用简单的共轭帽子结构，如CH₃CO-(N端)和NCH₃(C-端)。共轭帽子大小的选择决定其精度和速度，使用大的共轭帽子，精度会提高，但是计算量就会增加；而使用小的共轭帽子，计算量会降低，但精度就会降低。

在目前的方案中我们使用MFCC的方式来计算可极化的蛋白质的原子电荷，计算路线如图2所示:首先我们采用MFCC的方式将蛋白质划分成N个片段，然后利用量化方法分别计算蛋白质所有残基的量化性质(如静电势)，利用限制性静电势拟合(RESP)的方法拟合出分子碎片中的原子电荷。为了考察溶剂化效应，在计算完体系每个原子的电荷后，利用软件Delphi来计算这些蛋白质所产生的溶剂化电荷，然后再将这些溶剂化电荷代入下一轮的量化计算中，经过迭代自洽后就可以得到这种量化计算的考虑了蛋白质具体环境的可极化电荷了。

为了能获得更加具有统计效应的模拟，使取样更加充分，采用最近发展起来的副本互换的方法来进行动力学模拟。对于一个新的蛋白体系，首先为了确定所选的温度范围和副本个数，对体系进行一系列不同温度下的全原子动力学模拟，以确定蛋白质所需跨过的能垒，从而确定每个副本的温度。

具体实现方案如图1所示。首先将整个蛋白质体系处于N个不同的温度(N个副本)下进行分子动力学模拟，每隔一段动力学模拟时间比如1ps，相邻的两副本之间进行通讯，在两个副本(不同温度)的分子动力学模拟中，分别得到在高温和低温下分子动力学模拟后的不同蛋白质构象的能量，使用波尔兹曼权重将二者作差；如果该差值大于阈值(随机产生的一个数值)，则需要交换温度，否则不需要交换温度。

判断是否需要交换温度具体为：在两个具有不同温度的分子动力学模拟中，分别得到在两个副本下分子动力学模拟后具有不同构象的蛋白质的能量，将两个波尔兹曼权重因子分别乘以两个不同构象的蛋白质的能量后作差；然后根据模特卡洛的大都会(Metropolis)准则随机产生一个阈值，如果该差值大于阈值，则需要交换温度，否则不需要交换温度。判断结束后，无论是否交换温度，都更新各副本的分子力场文件中的蛋白质各原子电荷信息。

将不同温度下所获取的具有不同构象的蛋白质结构均在肽键处切开，在每个切开的地方引进一对共轭的帽子，然后将每个氨基酸分子碎片进行量化计算拟合量化计算的静电势，从而得到分子碎片中的每个原子的电荷，通过更新旧分子力场中相应的原子电荷后，新的分子力场可以作为下一次的分子动力学模拟的输入，继续动力学模拟，从而不断循环下去，直至整个分子动力学模拟结束。而在我们的方案中为了考虑蛋白质极化效应，我们拟在判断完温度是否交换后，马上更新电荷，从而让接下来的动力学模拟在这个更加真实的电性环境下进行模拟，可以获得更加准确的模拟数据。

此种方案(图3)将所需要的计算资源分成两部分，一部分为跑动力学部分，另一部分是用来做量化计算的部分。假设每个副本进行动力学计算需要K个核，那么纯粹的在副本互换的动力学模拟中所需要的核数就有K*N个，由于我们需要计算每个蛋白质的极化电荷，那么通过MFCC方案，假设有L个残基的蛋白质，每个残基的量化计算我们通过M个核来计算，那么每个蛋白质构象所需要的量化计算量就是M*L,而我们一共有N个副本，则花在量化计算中的核数就有N*L*M，综上，对于我们具有L个残基的蛋白质，使用N个副本，我们所需要的计算的核数就是N*(K+L*M)。

如图1所示，为了将共轭帽分子分割法(MFCC)和基于副本互换的方法(REMD)，将分子动力学软件(AMBER),共轭帽分子分割法软件(MFCC)和量子力学软件(NWChem)有机结合起来，发展出AMNW软件平台。根据超级计算中心的具体架构以及需要的计算核数、具体的模拟时间长度、具体的副本数目以及温度范围(一般为100K–600K)，AMNW软件平台设置相应的提交分子动力学模拟和量子力学计算的壳(shell)脚本，并且产生分子动力学模拟软件(AMBER)所需要的输入文件，然后处理动力学模拟的输出的蛋白质的结构信息文件，通过共轭帽分子分割法软件(MFCC)将蛋白质的结构信息文件分解成氨基酸分子碎片文件作为量子力学计算软件(NWchem)所需要的输入文件，量子力学计算完成后，AMNW软件平台从量子力学计算的结果中提取出各原子的电荷数据，最后将这些新原子电荷数据更新到原有的分子力场文件中。为了使得共轭帽分子分割法软件(MFCC)与量子力学软件(NWChem)能够结合使用，我们对共轭帽分子分割法软件(MFCC)中调用量化计算的模块进行了相应的代码修改。我们将AMNW软件平台部署在超算中心进行了并行效率的测试，如图4所示，以克服分子动力学模拟中的原子相互作用势和取样技术这两大难题，能在蛋白质折叠模拟中获得比较好的效果。

首先，从AMNW软件平台可以简单地调用分子动力学软件(AMBER)对蛋白质执行副本互换的分子动力模拟，在不同温度(副本)下，副本互换可以帮助在更有效地搜索蛋白质的构象空间，从而达到尽可能多获取不同的蛋白质构象，增加蛋白质构象取样的完整性。其次，AMNW软件平台调用共轭帽分子分割法软件(MFCC)可对副本互换动力学模拟产生的蛋白质结构进行分子碎片分割，为下一步量子力学的计算准备输入文件。然后，AMNW软件平台调用量子力学软件(NWchem)对各蛋白质结构的分子碎片进行量化计算，最终得到每个原子的具体电荷，这样可以充分体现每个原子的电荷随着环境变化而变化，从而比固定点电荷模型更能准确描述原子间的静电相互作用。最后，AMNW软件平台可以通过调用我们自主完成的电荷更新程序将分子动力模拟软件(AMBER)的旧分子力场文件中各原子的电荷进行更新，然后AMNW软件平台再一次调用分子动力学软件(AMBER)可对蛋白质执行副本互换的分子动力模拟，如此循环，直到预先设定的整个分子动力学模拟(如几个纳秒)。

AMNW平台的并行效率测试：

现已将AMNW蛋白质折叠软件平台部署在各超算中心，包括中科院北京深腾超算中心，上海魔方超算中心以及国防科技大设计天河II号超算中心。首先，在中科院北京深腾超算中心上对AMNW蛋白质折叠软件平台进行并行效率的测试。在测试中我们选用20个氨基酸的Trpcage蛋白作为测试体系，并使用8个副本(不同温度)，每500步进行副本互换，然后AMNW软件平台可以调用量子力学软件(NWchem)更新每个构象的电荷，如此循环5次，整个分子动力学模拟进行了5ps。在并行效率测试中，我们分别使用了128个核，256个核，512个核和1024个核来计算整个分子动力学模拟所需要的时间，我们发现千核并行效率在80％以上(见图5)。

然后在国防科技大设计天河II号超算中心对AMNW蛋白质折叠软件平台进行并行效率的测试。在测试中我们选用的测试体系是具有106个氨基酸的Fe(II)cytochrome b562蛋白质体系，并使用32个温度(副本)，每500步进行副本互换，然后AMNW软件平台可以通过相关模块调用量子力学软件(NWchem)更新每个构象的电荷，整个分子动力学模拟进行了5ps。在并行效率测试中，我们分别使用了3392个核，6784个核，13568个核,27136个核,54274个核,108544个核,217088个核和284928个核来计算整个分子动力学模拟所需要的时间，发现三十万核量级的并行效率可以达到30％，符合我们预期的结果。如图6a和6b所示，分别表示AMNW蛋白质折叠软件平台在天河II号上并行效率和加速比的测试结果。

Claims (3)

1.蛋白质折叠的并行预测方法，其特征在于包括以下步骤：

1)通过分子动力学模拟得到蛋白质结构；

2)设定不同温度，通过副本互换得到在不同温度下的蛋白质不同构象；

3)采用共轭帽分子分割法对蛋白质不同构象进行分割，将每个氨基酸从蛋白质中分解出来，并且每个氨基酸的分子碎片带有共轭帽；

4)通过量子力学软件计算分子碎片的静电势，通过拟合得到每个氨基酸分子碎片的原子电荷；通过更新分子力场中相应的原子电荷后返回至步骤1)，直到达到设定时间为止；

在计算完体系每个原子的电荷后，利用软件Delphi来计算这些蛋白质所产生的溶剂化电荷，然后再将这些溶剂化电荷代入下一轮的量化计算中，经过迭代自洽后就得到这种量化计算的考虑了蛋白质具体环境的可极化电荷；

所述设定不同温度，通过副本互换得到在不同温度下的蛋白质不同构象包括以下步骤：

对每个分子动力学模拟设定N个不同的温度，也就是N个副本；

在分子动力学模拟过程中，在N个副本中任意相邻的两个副本之间进行通讯，判断是否需要交换温度；如需要交换温度，则两个副本交换温度；否则两个副本不交换温度；得到交换温度后不同温度下的蛋白质不同构象。

2.根据权利要求1所述的蛋白质折叠的并行预测方法，其特征在于所述判断是否需要交换温度包括以下步骤：

在两个具有不同温度的分子动力学模拟中，分别得到在两个副本下分子动力学模拟后具有不同构象的蛋白质的能量，使用波尔兹曼权重因子将二者作差；然后随机产生一个阈值，如果差值大于阈值，则需要交换温度，否则不需要交换温度。

3.根据权利要求1所述的蛋白质折叠的并行预测方法，其特征在于所述采用共轭帽分子分割法对蛋白质不同构象进行分割包括以下步骤：

将不同温度下的蛋白质均在肽键处切开，在每个切开的地方引进一对共轭帽子，该共轭帽子为与之相邻的两个或四个氨基酸残基，或者N端使用CH₃CO-、C端使用N(CH₃)₂的共轭帽子结构。