多尺度粗粒化建模方法及分子动力学模拟方法

摘要

一种多尺度粗粒化建模方法及分子动力学模拟方法，该建模方法包括将由多个单体键合而成的生物大分子按照结构特征划分为多个结构区域；对所述多个结构区域进行不同精度的粗粒化建模，分别得到与多个所述结构区域对应的多个粗粒化模型，所述多个粗粒化模型包括以单个或多个球体表征相互作用位点的相互作用模型；将与多个所述结构区域对应的多个所述粗粒化模型以弹簧键链接，共同组成所述生物大分子的多尺度粗粒化模型。本发明基于该多尺度粗粒化模型进行大规模分子动力学模拟，可比全原子分子动力学或单体层次粗粒化模型的分子动力学提高数倍至数百倍的模拟效率，能够实现对大规模统计行为更大时间和空间尺度的模拟。

说明书

技术领域

本发明涉及分子动力学技术领域，具体涉及一种多尺度粗粒化建模方法及分子动力学模拟方法。

技术背景

细胞由大量生物大分子构成，如蛋白质、核酸、糖类等，它们组成细胞结构、实现细胞内的各种功能。然而大量的生物大分子如何共同实现结构和功能、如何组织起复杂细胞内的生化反应依然是重要的科学问题，受到了人们的广泛关注。例如为创造适合特定生化反应的间隔区域，生物分子会区域化形成无膜细胞器，例如核仁、P颗粒、卡哈尔体、应激颗粒等等，他们的形成依赖于生物大分子的液液相分离过程；此外中心粒生长、微观生长、染色质区域化等过程电在细胞中发挥重要作用。这些细胞内生命过程不断得到新的认识，但是如何从本质上理解这些过程的形成机理、如何理解突变、环境因素对各种过程的影响乃至对功能的作用，仍然在实验上具有巨大的挑战。

分子动力学模拟是理论研究的重要手段，已成为研究细胞内分子机理的关键步骤。在过去的十年中，随着计算机水平的迅速发展，理论化学家已经可以进行全原子动力学模拟以研究分子量达上百万的蛋白的构象转变、受体识别等等过程。因此，分子动力学模拟为理解由围观生物分子到宏观细胞内现象的关系提供了一种有前景的方法。

然而对于上述相分离、聚集等统计行为，由于需涉及到成千上万个生物大分子，若使用全原子动力学模拟在，当前技术下计算成本极高且难以实现。因此，采用更粗粒化的建模模拟十分必要。发明人发现，现有技术对生物大分子的粗粒化方法一般只涉及单一尺度，主要分为两类：一类是将整个蛋白或RNA等粗粒化为一个球体，这种处理无法反映分子内更多结构信息，和反映相互作用位点，难以在研究突变等问题时实现准确刻画；另一种是将大分子内每个单体粗粒化为一个球体，从而构成一个软的球珠链，这种处理方法能够反映序列信息，但对于大量的局部有序的大分子的高级结构刻画不足，且面临计算量依然巨大的问题。因此在本发明中，通过对生物大分子多尺度粗粒化建模，可以在抓住大分子主要结构特征和作用位点分布的同时，实现大规模高效动力学模拟，在研究生物大分子统计行为的机理中发挥重要作用。

发明内容

有鉴于此，本发明的主要目的在于提供了一种多尺度粗粒化建模方法及分子动力学模拟方法，以期至少部分地解决上述提及的技术问题中的至少之一。

为实现上述目的，本发明的技术方案如下：

根据本发明的一个方面，提供了一种生物大分子的多尺度粗粒化建模方法，包括：

将由多个单体键合而成的生物大分子按照结构特征划分为多个结构区域；对所述多个结构区域进行不同精度的粗粒化建模，分别得到与多个所述结构区域对应的多个粗粒化模型，所述多个粗粒化模型包括以单个或多个球体表征相互作用位点的相互作用模型；将与多个所述结构区域对应的多个所述粗粒化模型以弹簧键链接，共同组成所述生物大分子的多尺度粗粒化模型。

根据本发明的另一个方面，提供了一种基于多尺度粗粒化建模的分子动力学模拟方法，包括：采用如上所述的多尺度粗粒化建模方法建立多尺度粗粒化模型；基于所述多尺度粗粒化模型和模拟体系大小生成动力学模拟所需的体系初始构型；基于所述体系初始构型采用分子动力学运动方程进行体系演化，得到关于所述多尺度粗粒化模型的模拟结果。

基于上述技术方案，本发明的多尺度粗粒化建模方法及分子动力学模拟方法至少具有以下有益效果其中之一或其中一部分：

(1)选用该多尺度粗粒化模型进行大规模分子动力学模拟，将生物大分子划分成不同结构区域，针对不同结构区域的特点进行不同精度的粗粒化建模，可比全原子分子动力学、或单体层次粗粒化模型的分子动力学提高数倍至数百倍的模拟效率，能够实现对大规模统计行为更大时间和空间尺度的模拟。

(2)本发明基于多尺度粗粒化模型生成模拟体系初始构型，通过选用模型适配的动力学方程对体系演化进行模拟，从而可得到相变、聚集等统计形成的形成机理及实验难以得到的分子间统计信息，并通过参数依赖为相应生物体系实验提供高效指导和方案。

附图说明

图1是本发明实施例1中的EB1蛋白序列结构信息及各结构区域的粗粒化示意图；

图2是本发明实施例1中的EB1蛋白的多尺度粗粒化模型结构图；

图3是本发明实施例1中的基于EB1蛋白的多尺度粗粒化模型在分子动力学模拟的平衡态的体系状态图；

图4A是本发明实施例1中的基于EB1蛋白的多尺度粗粒化模型在分子动力学模拟过程中的液滴生长动力学统计图；

图4B是本发明实施例1中的基于EB1蛋白的多尺度粗粒化模型在分子动力学模拟过程中的分子内相互作用占比统计图；

图5是本发明实施例2中的G3BP1蛋白序列结构信息及各结构区域的粗粒化示意图；

图6是本发明实施例2中的G3BP1蛋白的多尺度粗粒化模型结构图；

图7是本发明实施例2中的基于G3BP1蛋白的多尺度粗粒化模型在分子动力学模拟的平衡态的体系状态图；

图8是本发明实施例2中的基于G3BP1蛋白的多尺度粗粒化模型在分子动力学模拟的平衡态时与RNA的平均价态统计。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

本发明针对当前生物分子粗粒化模拟中模拟效率低、对生物大分子难以实现准确刻画的问题，本发明提供了一种多尺度粗粒化建模方法以及基于多尺度粗粒化建模的分子动力学模拟方法。具体而言，以生物大分子的粗粒化模拟为核心，根据生物分子的结构特点对分子或其多聚体的不同结构区域进行不同精度的粗粒化建模，再在与结构区域对应的粗粒化模型中建立相互作用位点的模型，为动力学模拟提供最高效的模型。并进一步地，基于多尺度粗粒化模型生成确定模拟体系初始构型，通过选用模型适配的动力学方程对体系演化进行模拟，从而得到相变、聚集等统计行为的形成机理及实验难以得到的分子间统计信息，并通过参数依赖为相应生物体系实验提供高效指导和方案。

具体而言，根据本发明的实施例，提供了一种生物大分子的多尺度粗粒化建模方法，包括步骤A～C。

在步骤A中，将由多个单体键合而成的生物大分子按照结构特征划分为多个结构区域。

在一些实施例中，多个结构区域独立地为无序软链区域或含二级或高级结构的有序结构区域。

以蛋白质为例，其为由氨基酸单体键合而成的生物大分子，蛋白质分子的部分区域自身无法折叠为明确唯一的空间结构，即为无序软链区域，部分区域为包含二级或高级结构的有序结构区域，可折叠为有序结构，其中二级结构例如α螺旋、β折叠、β转角，三级结构例如在二级结构基础上多肽链进一步盘绕、折叠形成的特定空间结构。

在步骤B中，对多个结构区域进行不同精度的粗粒化建模，分别得到与多个结构区域对应的多个粗粒化模型，多个粗粒化模型包括以单个或多个球体表征相互作用位点的相互作用模型。

在一些实施例中，在结构区域为无序软链区域时，以单个球体表征单体，例如蛋白质的氨基酸单体，形成多个球体相连而成的球珠软链模型；在结构区域为有序结构区域时，以单个球体或多个球体形成与有序结构区域形态近似的刚体模型。

需要说明的是，球珠软链模型或者刚体模型是对于分子动力学计算而言，可基于相互作用势以及动力学方程的计算得以区分，具体而言，球珠软链模型内的两个球体之间的相互作用势设置为弹簧，表明两个球体之间可相对移动，从而形成软链，而刚体模型内两个球体之间不设相互作用或合力为0，并且在分子动力学方式求解时，刚体模型内每个球体的位置是根据其相对于刚体模型质心的相对位置来进行更新。

更具体地，在结构区域为有序结构区域时，若有序结构区域在近似为单个球体时能够反应主要物理性质，则建模为单个球形刚体模型，否则根据有序结构区域的形状特点用多个球体勾勒有序结构区域的轮廓，并组成刚体。

在一些实施例中，可基于分子动力学模拟过程中所关注的作为重要因素的相互作用来建立相互作用模型。其中，用于表征相互作用位点的单个或多个球体位于球珠软链模型内或位于刚体模型表面，其中每个球体用于表征单个相互作用位点。

需要说明的是，用于表征相互作用位点的单个或多个球体位于球珠软链模型内，表示球珠软链模型内中用于表征单体的单个或多个球体同时用于表征相互作用位点。这些用于表征相互作用位点的小球根据所关注的相互作用可以附上电荷或者相互作用势等，以便进行分子动力学计算。

在一些实施例中，在存在反映生物大分子的单体之间相互作用强度的核磁数据时，则基于生物大分子的晶体衍射结构及核磁数据确定单个或多个球体的位置及大小；举例而言，以蛋白质分子为例，基于核磁数据可以确定蛋白质分子中能够产生关键相互作用的氨基酸单体，从而结合蛋白质分子的晶体结构将这些氨基酸单体以单个或多个球体建立相互作用模型。

或者，在不存在核磁数据时，基于生物大分子的晶体衍射结构以及模拟实验参数确定单个或多个球体的位置及大小。举例而言，结合晶体结构和模拟实验参数例如pH和离子浓度可用于确定蛋白质中氨基酸单体的带电量，基于这些氨基酸单体的带电量来确定表征相互作用位点的球体位置和大小。

在此基础上，在以单个球体作为生物分子建模的案例中，模型对生物分子反映不足，能够给出的统计信息较少，能与实验对照的因素较定性，预言能力较差，直接表现为模拟结果与实验误差较大。基于核磁、晶体衍射数据确定的粗粒化相互作用位点，与实际生物实验结果符合较好，说明该多尺度模拟过程准确度高。有助于为今后研究该基因突变对疾病发展机制的影响，以及以该基因为作用靶点治疗疾病等奠定基础。

在一些实施例中，若分子动力学模拟涉及细胞内过程时，若在细胞内过程发生前生物大分子已形成寡聚，则可以寡聚体为单位进行建模，具体而言，将两个以上所述生物大分子的多尺度粗粒化模型中与发生聚合的所述结构区域对应的粗粒化模型共同组成单个刚体模型。

在步骤C中，将与多个结构区域对应的多个粗粒化模型以弹簧键链接，共同组成生物大分子的多尺度粗粒化模型。

根据本公开的实施例，还提供了一种基于多尺度粗粒化建模的分子动力学模拟方法，包括步骤a～c。

在步骤a中，采用如上所述的多尺度粗粒化建模方法建立多尺度粗粒化模型。

在步骤b中，基于多尺度粗粒化模型和模拟体系大小生成动力学模拟所需的体系初始构型。

在一些实施例中，步骤b包括子步骤b1～b2：

在子步骤b1中，根据模拟体系大小建立具有预设周期性边界的模拟盒子，模拟盒子内包含预设数量的生物大分子的多尺度粗粒化模型；

在子步骤b2中，采用与生物大分子的预设数量相同的循环次数，循环确定每个生物大分子的多尺度粗粒化模型内每个球体的初始坐标。

具体而言，首先确定多尺度粗粒化模型中首个球的初始坐标，然后基于各结构区域所对应的粗粒化模型的连接顺序，通过随机行走(randomwalk)得到多尺度粗粒化模型中其余球体的初始坐标。

在一些实施例中，位于刚体模型上的每个球体的初始坐标是基于刚体模型质心的初始坐标及其与刚体模型质心的相对位置而确定的。

在步骤c中，基于体系初始构型采用分子动力学运动方程进行体系演化，得到关于多尺度粗粒化模型的模拟结果。

在一些实施例中，步骤c包括子步骤c1～c6：

在子步骤c1中，根据多尺度粗粒化模型内每个球体的初始坐标、初始速度和相互作用势截断范围，建立模拟体系内所有球体的全邻居列表粒子关系矩阵。

在子步骤c2中，根据全邻居列表粒子关系矩阵，依次计算多尺度粗粒化模型中存在相互作用关系的球体的受力，从而累计计算每个球体所受的合力；

在子步骤c3中，将与有序结构区域相对应的刚体模型中所有球体的合力进行四元数运算，得到刚体模型所受的合力；

在子步骤c4中，根据刚体模型所受合力，及球珠软链模型中每个球体所受的合力，计算所述刚体模型以及球珠软链模型中每个球体所受到的朗之万热浴的随机力以及随机转矩；

在子步骤c5中，根据刚体模型及球珠软链模型中所有球体所受的合力、随机力以及随机转矩代入朗之万方程中，求解得到多尺度粗粒化模型内所有球体的位移和旋转速度；

在子步骤c6中，利用Verlet速度积分算法计算和更新多尺度粗粒化模型内所有球体的速度和位置信息，刚体模型内每个球体的位置根据其与所述刚体模型质心的相对位置进行更新。

在上述子步骤c1～c6中，全邻居列表离子关系矩阵的建立、朗之万热浴下随机力及随机转矩的计算、朗之万方程的求解过程已为本领域技术人员所熟知，故在此不作赘述。本发明与常规动力学模拟中体系演化的区别主要在于，基于多尺度粗粒化模型的体系演化过程中同时包含了刚体模型和球珠软链模型关于球体所受合力以及位置坐标的计算。

以下列举多个具体实施例来对本发明的技术方案作详细说明。需要说明的是，下文中的具体实施例仅用于示例，并不用于限制本发明。在具体实施例中模拟所使用的软件为自主编程大规模粗粒化并行平台。所有晶体结构来源于Protein Data Bank (PDB)。

实施例1：

本实施例对EB1蛋白及其突变体相分离行为进行多尺度建模并进行分子动力学模拟，主要包括以下步骤：

(1)获取EB1的晶体结构，根据结构特点进行粗粒化建模。EB1蛋白由C端至N端主要包含4个结构区域，如图1A所示：CH结构区域(CH domain)、无序的Linker区域、卷曲螺旋(C-coil)区域和无序的c-tail区域。CH结构区域(CH domain)可被近似为一个球状的刚体模型，如图1C所示；C-coil结构区域根据其棒状的结构特点，用约10个球体勾勒其轮廓，如图1D所示，这些球体共同构成分子动力学计算中的一个刚体模型；Linker无序区域(共70个氨基酸)为体现其内禀无序性特征，采用以单个氨基酸为精度的球珠软链模型，如图1B所示；c-tail无序区域采用与Linker无序区域相类似的建模方式，采用单个氨基酸为精度的球珠软链模型。由于在关注的相分离过程发生时，EB1通常首先发生二聚，故进一步将EB1二聚体作为粗粒化模拟的建模单位，将二聚体的两个单体中c-coil部分共同组成单个刚体，以反映其二聚特点。

(2)确定所关注科学问题起重要因素的相互作用，建立相互作用位点的模型。以EB1体内发生相分离过程为例，其分子间相互作用位点分类：C端结构域及中间的无序区，主导相互作用为静电相互作用，则在CH结构区域的球状刚体模型上，参考表面带电的氨基酸相对位置(共四个正电荷位点，3个负电荷位点)，根据晶体结构中这些位点的相对位置确定粗粒化相互作用位点的球体在球状刚体模型表面的相对位置。对于无序结构区域内的带电氨基酸，则在模型中将对应于带电氨基酸的粗粒化球体带上相应的电荷。C-coil结构区域主要参与疏水相互作用，则在C-coil结构区的棒状刚体模型上每个球体作为相互作用位点并设定两EB1间的疏水相互作用势。由此，将对应于各结构区域的各粗粒化模型之间用弹簧键链接，共同组成多尺度粗粒化模型，如图2所示。

(3)确定模拟体系大小，生成动力学模拟所需的初始构型。以模拟600条EB1蛋白，浓度为1.0Mm为例，则需要周期性边界为100nm的模拟盒子，模拟盒子内含300个EB1蛋白二聚体的多尺度粗粒化模型。基于软件编程“for”循环，每次循环确定一个二聚体分子的全部坐标：循环时首先确定首个球的坐标，再根据步骤(1)、(2)中建模方法所建立的多尺度粗粒化模型，按照各结构区域的连接顺序通过随机行走(random walk)得到其余粗粒化结构中所有球体的坐标。其中位于刚体模型上的球体的初始坐标基于刚体质心及相对位置确定。

(4)对体系进行分子动力学模拟，对粗粒化的生物大分子采用布朗动力学，运用朗之万方程进行体系演化，直到体系演化到平衡态或所需时间。具体而言，首先根据多尺度粗粒化模型内每个球体的初始坐标、初始速度和相互作用势截断范围，建立模拟体系内所有球体的全邻居列表粒子关系矩阵；然后，根据全邻居列表粒子关系矩阵，依次计算多尺度粗粒化模型中存在相互作用关系的球体的受力，从而累计计算每个球体所受的合力；然后，将与有序结构区域相对应的刚体模型中所有球体的合力进行四元数运算，得到所述刚体模型所受的合力；然后，根据刚体模型所受合力及球珠软链模型中每个球体所受的合力，计算刚体模型以及球珠软链模型中每个球体所受到的朗之万热浴的随机力以及随机转矩；然后，根据刚体模型及球珠软链模型中所有球体所受的合力、随机力以及随机转矩代入朗之万方程中，求解得到多尺度粗粒化模型内所有球体的位移和旋转速度；最后利用Verlet速度积分算法计算和更新多尺度粗粒化模型内所有球体的速度和位置信息，刚体模型内每个球体的位置根据其与刚体模型质心的相对位置进行更新。

(5)基于步骤(4)可得到模拟过程的粗粒化模型中各球体随时间运动轨迹，其到达平衡态的体系截图如图3所示。由此，可对模拟结果进行分析，统计EB1相分离的液滴生长动力学过程(如图4A所示)、EB1野生型和各突变体之间分子内相互作用比例差别(如图4B所示)等等以阐释相关机理。

实施例2：

本发明对G3BP1蛋白进行多尺度建模并进行分子动力学模拟，步骤为：

(1)获取G3BP1的晶体结构，根据结构特点进行粗粒化建模。G3BP1蛋白主要包含一个结构域：RNA识别结构域(RNArecognition motif，RRM)可被近似为一个球形刚体模型，如图5A所示；在RRM后面链接一段无序结构区域(RGG)，该段无序结构区域采用以单个氨基酸为精度的球珠软链模型，如图5B所示。

(2)确定所关注科学问题起重要因素的相互作用，建立相互作用位点的模型。以G3BP1结合RNA形成应激颗粒过程为例，其与RNA的主导相互作用为静电相互作用。根据实验的核磁数据，能够确定G3BP1和RNA结合RRM上有6个作用位点、在无序区有两个相互作用位点，则在RRM的球状刚体模型上，参考表面带电氨基酸相对位置(共6个正电荷位点)，根据晶体结构中这些位点的相对位置确定粗粒化相互作用位点的球体在球状刚体模型表面的相对位置。对于无序结构区域内的带电氨基酸，对核磁数据确定的参与相互作用氨基酸的粗粒化球体上带上相应的电荷。由此，将对应于各结构区域的各粗粒化模型之间用弹簧键链接，共同组成多尺度粗粒化模型，如图6所示。

(3)确定模拟体系大小，生成动力学模拟所需的提初始构型。以模拟600条G3BP1蛋白，浓度为1.0Mm为例，则需要周期性边界为100nm的模拟盒子。基于软件编程“for”循环，每次循环确定一个分子的全部坐标：循环首先确定首个球的坐标，再根据步骤(1)、(2)中建模方法所建立的多尺度粗粒化模型，按照各结构区域连接顺序通过随机行走(random walk)得到其余粗粒化结构中所有球体的坐标。其中位于刚体模型上的球体的初始坐标基于刚体质心及相对位置确定。

(4)对体系进行分子动力学模拟，对粗粒化的生物大分子采用布朗动力学，运用朗之万方程进行体系演化，直到体系演化到平衡态或所需时间。具体而言，首先根据多尺度粗粒化模型内每个球体的初始坐标、初始速度和相互作用势截断范围，建立模拟体系内所有球体的全邻居列表粒子关系矩阵；然后，根据全邻居列表粒子关系矩阵，依次计算多尺度粗粒化模型中存在相互作用关系的球体的受力，从而累计计算每个球体所受的合力；然后，将与有序结构区域相对应的刚体模型中所有球体的合力进行四元数运算，得到所述刚体模型所受的合力；然后，根据刚体模型所受合力及球珠软链模型中每个球体所受的合力，计算刚体模型以及球珠软链模型中每个球体所受到的朗之万热浴的随机力以及随机转矩；然后，根据刚体模型及球珠软链模型中所有球体所受的合力、随机力以及随机转矩代入朗之万方程中，求解得到多尺度粗粒化模型内所有球体的位移和旋转速度；最后利用Verlet速度积分算法计算和更新多尺度粗粒化模型内所有球体的速度和位置信息，刚体模型内每个球体的位置根据其与刚体模型质心的相对位置进行更新。

(5)基于步骤(4)可得到模拟过程的粗粒化模型中各球体随时间运动轨迹，对模拟结果进行分析，其到达平衡态的体系截图如图7所示。由此，可根据模拟结果进行机理分析，如统计其与RNA的平均价态，如图8所示。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。