生物分子功能性动力学多转变路径的计算方法及系统

摘要

本发明提出了一种生物分子功能性动力学多条转变路径的计算方法及系统，方法包括：确定起始态结构和目标态结构，并在隐含水条件下由起始态结构出发得到M个初始结构；进行一轮或多轮采样；基于每轮采样中的模拟轨迹构建M条转变路径，对M条转变路径进行聚类并选取每一类中的代表转变路径，得到N个不同的代表转变路径；以隐含水条件下的N个代表转变路径为参照，利用预设算法重建显含水条件下的N条转变路径。本发明的方案，创新性的采用了隐含水的分子动力学模拟，可大幅降低计算量并提高路径搜索效率。基于隐含水的多条转变路径，实现显含水转变路径的重建以及优化，并最终得到目标生物分子功能性动力学的转变机制。

说明书

技术领域

本发明涉及生物分子系统计算模拟研究领域，特别涉及生物分子功能性动力学多转变路径的计算方法及系统。

背景技术

生物分子实现其功能时，常伴随着结构的复杂转变（亦被称为生物分子功能性动力学），这种结构转变往往存在多种不同的转变方式，而每种转变方式的具体过程可被称为转变路径。对于每条转变路径，还需进一步优化以搜寻到相应的最小自由能路径（MinimumFree Energy Path），而最小自由能路径所反映的是与实际情况最相近的微观转变机制。因而，对于复杂的生物分子系统的功能性动力学研究，要实现对其微观机制的全面深入理解，则必需设法获取其主要的多条转变路径，并优化至最小自由能路径，而基于分子动力学（Molecular Dynamics，MD）模拟进行增强采样算法开发成为解决上述挑战的可行方案之一。

常规MD模拟是一种基于牛顿动力学来研究生物分子微观性质的计算模拟研究方法，其目标为模拟生物分子在热扰动下随时间变化的结构涨落行为，因而常规MD模拟中结构转变的难易程度与能垒差别对应。对于生物分子功能性动力学，其所涉及的结构转变较为复杂，且转变过程中需翻越较高能垒，这使得由常规MD模拟和使用通用计算设备在有限的时间内观测到上述过程的可能性极低；同时，生物分子功能性动力学转变过程还存在着多种转变路径；因此，已有的方案都试图通过将常规MD模拟结合增强采样算法，来解决上述挑战，以获取生物分子功能性动力学的多条最优转变路径（最小自由能路径）。

现有技术中，第一种方案，通过得到一条结构转变路径，沿着此参考转变路径进行分割，得到多个采样空间，随后从每个采样空间中随机挑选初始结构，重复完成多个短时常规MD采样，直至采样结构所分属于任意采样空间的概率都趋于稳定；从所有的采样数据中，通过提取和整合采样空间相互之间的大量转变路径，进而筛选出生物分子功能性动力学相应的多种不同转变路径。

现有技术中，第二种方案，从已有的大量路径中选取任意转变路径，从路径所有的结构出发，完成短时MD采样，并随机的在不同转变路径之间交换路径的某部分结构，或者随机生成某部分结构的维度信息，从而使转变路径的作用量发生改变，通过筛选去除冗余，最终获取生物分子功能性动力学的多条转变路径。

首先，现有技术的方案需在显含水分子条件下完成大量常规MD模拟采样，而MD模拟采样中计算资源主要便是消耗于处理水分子的动力学演化计算，而这对于较大的生物分子系统，研究其功能性动力学复杂转变过程，便要消耗巨量的计算成本和时间成本。

其次，对于研究生物分子功能性动力学复杂的结构转变：若采用第一种方案，由于转变路径中生物分子的结构自由度过多，则将导致沿转变路径划分的采样空间数量过多，这也极大增加了其最终实现足够采样的难度。若采用第二种方案，考虑到生物分子结构中可调整的变量过多，使得该方法通过部分结构随机互换或者随机产生部分结构信息可获取的新转变路径数量大大增加。

发明内容

有鉴于此，本发明提出了一种生物分子功能性动力学多转变路径的计算方法及系统，具体方案如下：

一种生物分子功能性动力学多转变路径的计算方法，包括如下：

确定生物分子功能性动力学的起始态结构和目标态结构，并在隐含水条件下由所述起始态结构出发得到M个初始结构；其中，M≥1；

进行一轮或多轮采样，每一轮采样包括：从M个初始结构出发，产生预设长度的M条模拟轨迹；随机生成当前概率，通过判断当前概率是否大于预设截断概率以从M条模拟轨迹中挑选M个构象；对挑选的M个构象进行预处理，并将预处理后的结果作为下一轮采样中的M个初始结构，完成一轮采样；

基于每轮采样中的模拟轨迹构建M条转变路径，对M条转变路径进行聚类并选取每一类中的代表转变路径，得到N个不同的代表转变路径；

以隐含水条件下的N个代表转变路径为参照，利用预设算法重建显含水条件下的N条转变路径。

在一个具体实施例中，在隐含水条件下得到M个初始结构，具体包括：

获取转变路径的期望数量M；

由所述起始态结构出发，于隐含水条件下产生初始模拟轨迹；

从所述初始模拟轨迹的所有构象中，挑选M个与所述目标态结构距离最近的结构，得到M个初始结构。

在一个具体实施例中，“随机生成当前概率，通过判断当前概率是否大于预设截断概率以从M条模拟轨迹中挑选M个构象”具体包括：

随机生成当前概率，判断当前概率是否大于预设截断概率：

若是，则从M条模拟轨迹的所有构象中筛选出

若否，则从M条模拟轨迹的所有构象中筛选出势能大于预设截断势能的高势能构象，并通过预设第一算法从所有高势能构象中挑选M个构象。

在一个具体实施例中，所述预处理具体包括：对挑选的M个构象进行排序，以使得累加的第i个构象与前一轮采样中第i条模拟轨迹的最后一帧构象之间RMSD距离和为最小；

其中，i=1,2,3……，M。

在一个具体实施例中，按照从第一轮到最后一轮的采集顺序，将每轮采样中第i条模拟轨迹整合为一条完整的转变路径，最终得到M条转变路径；

计算M条转变路径中任意两条转变路径的相似度，当两条转变路径之间的相似度大于相似度阈值时，此两条路径被归为一类；

完成所有转变路径之间的相似度计算，完成对M条转变路径的聚类。

在一个具体实施例中，所述代表转变路径的获取包括：

对于被聚为一类的所有转变路径，通过计算该类中所有转变路径的OM作用量，并筛选出OM作用量最小的转变路径作为该类的代表转变路径。

在一个具体实施例中，重建得到显含水条件下的N条转变路径后，还包括：

利用基于旅行商问题的自动路径搜索方法优化显含水条件下的N条转变路径，得到最小自由能转变路径；

利用伞形采样方式，计算最小自由能转变路径的自由能面。

在一个具体实施例中，显含水条件下的N条转变路径的获取过程包括：

以隐含水条件下的N个代表转变路径为参照，从显含水的初始结构系统出发，分别沿着生物分子功能性动力学的不同转变路径，利用定向分子动力学模拟算法重建得到显含水条件下的N条转变路径。

在一个具体实施例中，在进行一轮采样过程中，当挑选的M个构象中，与目标态结构的最大RMSD小于预设目标截断距离时，则当前采样为最后一轮采样。

一种生物分子功能性动力学多转变路径的计算系统，包括如下：

初始单元，用于确定生物分子功能性动力学的起始态结构和目标态结构，并在隐含水条件下，由所述起始态结构出发得到M个初始结构；

采样单元，用于进行一轮或多轮采样，每一轮采样包括：从M个初始结构出发，产生预设长度的M条模拟轨迹；随机生成当前概率，通过判断当前概率是否大于预设截断概率以从M条模拟轨迹中挑选M个构象；对挑选的M个构象进行预处理，并将预处理后的结果作为下一轮采样中的M个初始结构，完成一轮采样；

聚类单元，用于基于每轮采样中的模拟轨迹构建M条转变路径，对M条转变路径进行聚类并选取每一类中的代表转变路径，得到N个不同的代表转变路径；

重建单元，用于以隐含水条件下的N个代表转变路径为参照，利用预设算法重建显含水条件下的N条转变路径。

在一个具体实施例中，所述采样单元具体包括：

随机概率模块，用于随机生成当前概率，判断当前概率是否大于预设截断概率；

第一挑选模块，用于在当前概率大于截断概率时，从M条模拟轨迹的所有构象中筛选出

第二挑选模块，用于在当前概率小于等于截断概率时，从M条模拟轨迹的所有构象中筛选出势能大于预设截断势能的高势能构象，并通过预设第一算法从所有高势能构象中挑选M个构象；

预处理模块，用于对挑选的M个构象进行排序，以使得累加的第i个构象与前一轮采样中第i条模拟轨迹的最后一帧构象之间RMSD距离和为最小；其中，i=1,2,3……，M。

在一个具体实施例中，所述聚类单元具体包括：

路径转换模块，用于按照从第一轮到最后一轮的采集顺序，将每轮采样中第i条模拟轨迹整合为一条完整的转变路径，最终得到M条转变路径；

路径聚类模块，用于计算M条转变路径中任意两条转变路径的相似度，当两条转变路径之间的相似度大于相似度阈值时，此两条路径被归为一类；完成所有转变路径之间的相似度计算，完成对M条转变路径的聚类；

筛选代表模块，用于对被聚为一类的所有转变路径，通过计算该类中所有转变路径的OM作用量，并筛选出OM作用量最小的转变路径作为该类的代表转变路径。

在一个具体实施例中，还包括：

优化单元，用于利用基于旅行商问题的自动路径搜索方法优化显含水条件下的N条转变路径，得到最小自由能转变路径；利用伞形采样方式，计算最小自由能转变路径的自由能面。

有益效果：

本发明提供了一种生物分子功能性动力学多转变路径的计算方法及系统，创新性的采用了隐含水的分子动力学模拟，可大幅降低计算量并提高路径搜索效率，采样难度低。基于隐含水的多转变路径，实现显含水转变路径的重建以及优化，并最终得到目标生物分子功能性动力学的转变机制，能够稳定应用于复杂的生物分子功能性动力学机制研究。

附图说明

图1为本发明实施例的计算方法流程图；

图2为本发明实施例的具体计算方法流程示意图；

图3为本发明实施例的基于路径相似度的路径聚类示意图；

图4为本发明实施例的OM作用量计算示意图；

图5为本发明实施例的计算系统模块示意图。

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

附图标记：1-初始单元；2-采样单元；3-聚类单元；4-重建单元；5-优化单元；21-随机概率模块；22-第一挑选模块；23-第二挑选模块；24-预处理模块；31-路径转换模块；32-路径聚类模块；33-筛选代表模块。

具体实施方式

在下文中，将更全面地描述本发明公开的各种实施例。本发明公开可具有各种实施例，并且可在其中做出调整和改变。然而，应理解：不存在将本发明公开的各种实施例限于在此公开的特定实施例的意图，而是应将本发明公开理解为涵盖落入本发明公开的各种实施例的精神和范围内的所有调整、等同物和/或可选方案。

构象本质上就是一种结构，本申请将构象和结构在不同阶段进行了适应性的区分。实质上，构象可理解为结构。

在本发明公开的各种实施例中使用的术语仅用于描述特定实施例的目的并且并非意在限制本发明公开的各种实施例。如在此所使用，单数形式意在也包括复数形式，除非上下文清楚地另有指示。除非另有限定，否则在这里使用的所有术语(包括技术术语和科学术语)具有与本发明公开的各种实施例所属领域普通技术人员通常理解的含义相同的含义。所述术语(诸如在一般使用的词典中限定的术语)将被解释为具有与在相关技术领域中的语境含义相同的含义并且将不被解释为具有理想化的含义或过于正式的含义，除非在本发明公开的各种实施例中被清楚地限定。

实施例1

本发明实施例1公开了一种生物分子功能性动力学多转变路径的计算方法，采用了隐含水的分子动力学模拟，基于隐含水的多转变路径，实现显含水转变路径的重建以及优化，并最终得到目标生物分子功能性动力学的转变机制。计算方法流程图如说明书附图1所示，具体方案如下：

一种生物分子功能性动力学多转变路径的计算方法，包括如下步骤：

101、确定生物分子功能性动力学的起始态结构和目标态结构，并在隐含水条件下，由起始态结构出发得到M个初始结构；

102、进行一轮或多轮采样，每一轮采样包括：

从M个初始结构出发，产生预设长度的M条模拟轨迹；随机生成当前概率，通过判断当前概率是否大于预设截断概率以从M条模拟轨迹中挑选M个构象；对挑选的M个构象进行预处理，并将预处理后的结果作为下一轮采样中的M个初始结构，完成一轮采样；

103、基于每轮采样中的模拟轨迹构建M条转变路径，对M条转变路径进行聚类并选取每一类中的代表转变路径，得到N个不同的代表转变路径；

104、以隐含水条件下的N个代表转变路径为参照，利用预设算法重建显含水条件下的N条转变路径；

105、利用基于旅行商问题的自动路径搜索方法优化显含水条件下的N条转变路径，得到最小自由能转变路径；利用伞形采样方式，计算最小自由能转变路径的自由能面。

其中，M一般会小于N，由于路径转变过程较为复杂，需要经历多个（数目为N）连续的结构转变，而搜索多条转变路径时，是从N个结构中用聚类方法找到有代表性的M个初始结构（代表性是指尽可能的挑选有着足够差异的M个结构）。说明书附图2提供了具体的搜索流程图。

计算方法的初始准备，即产生短的隐水MD模拟轨迹并挑选M个距目标态结构RMSD最小的构象，将这些构象作为初始结构。获取转变路径的期望数量M；由起始态结构出发，于隐含水条件下产生初始模拟轨迹；从初始模拟轨迹的所有构象中，挑选M个与目标态结构距离最近的结构，得到M个初始结构。

确定生物分子功能性动力学的两个关键结构，分别作为起始态结构和目标态结构。由起始结构出发于隐含水条件下产生初始模拟轨迹。隐含水条件，即水分子的贡献由广义波恩模型近似描述。利用更加节省计算资源的隐含水MD模拟，且在隐含水条件下，生物分子更易发生结构转变，使得其采样效率更高。相较于已有的计算方法，本实施例创新性的采用了隐含水的分子动力学模拟，此举可提高隐水条件下构象转变的概率，大幅降低计算量并提高路径搜索效率。

从初始轨迹的所有构象中直接挑选M个与目标态结构差异性最小的初始结构，评估生物分子不同结构间的差异性则是通过计算构象与目标态结构之间的均方根偏差（RootMean Square Deviation，RMSD）来得到。

由M个结构产生M条隐水MD模拟轨迹，从M个初始结构出发，分别产生一定长度的模拟轨迹，模拟轨迹一般为20-200ps。

参考提前设置的截断概率P

若当前概率大于截断概率，则从M条模拟轨迹的所有构象中筛选出

若当前概率小于等于截断概率，则从M条模拟轨迹的所有构象中筛选出势能大于预设截断势能的高势能构象，并通过预设第一算法从所有高势能构象中挑选M个构象。优选地，第一算法为K-center算法，从M条模拟轨迹的所有高势能构象中由K-center算法挑选M个构象。

多转变路径搜索中，从每轮采样结果中筛选下一轮的初始结构时，即以P

当M个构象中与目标态结构的最大RMSD距离小于目标截断距离r

预处理具体包括：对挑选的M个构象进行排序，以使得累加的第i个构象与前一轮采样中第i条模拟轨迹的最后一帧构象之间RMSD距离和为最小，以公式表示为

其中，i=1,2,3……，M。

每一轮采样都可设置编号，按照从第一轮到最后一轮的采集顺序，将每轮采样中第i条模拟轨迹整合为一条完整的转变路径，最终得到M条转变路径。每轮采集都会得到M个构象，对每轮的M个构象进行排序，将所有编号相同的构象整合为一条转变路径。

得到M条转变路径之后，需要对其进行聚类。聚类的方式包括多种，如基于预定的反应坐标，评估结构间相似性。本实施例提供了一种聚类的方案，通过计算自定义的路径相似度实现快速聚类。具体包括：通过计算M条转变路径中任意两条转变路径的相似度，当两条转变路径之间的相似度大于相似度阈值时，此两条路径被归为一类。优选地，相似度阈值为0.6，即两条转变路径之间的相似度超过60%，则认定此两条转变路径为一类。完成所有转变路径之间的相似度计算，完成对M条转变路径的聚类。基于路径相似度的路径聚类如说明书附图3所示。

完成路径聚类后，需要进行最优路径筛选。在本实施例中，通过计算重新定义的OM作用量准确筛选出每类路径中的代表路径。具体包括：对于被聚为一类的所有转变路径，通过计算该类中所有转变路径的OM作用量，并筛选出OM作用量最小的转变路径作为该类的代表转变路径。代表转变路径即为该类中的最优路径。OM作用量计算示意图如说明书附图4所示。

进一步，基于隐含水条件下的路径计算获取显含水条件下的路径。以隐含水条件下的N个代表转变路径为参照，从显含水的初始结构系统出发，分别沿着生物分子功能性动力学的不同转变路径，利用标靶动力学模拟重建得到显含水条件下的N条转变路径。优选地，利用定向分子动力学模拟算法（Targeted Molecular Dynamics）重建显水条件下的N条转变路径。

重建得到显含水条件下的N条转变路径后，利用TAPS方法完成最终的转变路径优化。TAPS路径优化算法，即为基于旅行商问题的自动路径搜索方法，能够实现并行化和GPU支持，这使得mTAPS算法可快速基于隐含水的多转变路径，实现显含水转变路径的重建以及优化，并最终得到目标生物分子功能性动力学的转变机制。

此外，针对多条代表性路径的路径优化，可使用其他路径优化方法进行优化。需要说明的是，针对HePaCS的多转变路径采样阶段的MD采样条件可进行调整，如可直接于显含水条件下进行采样、或隐含水显含离子条件下进行采样、或粗颗粒化的系统进行采样。

本实施例方案中所构建的生物分子路径采样，相较于现有的计算方法，创新性的采用了隐含水的分子动力学模拟，由于隐水条件下构象转变的概率提高，可大幅降低计算量并提高路径搜索效率。正因如此，本实施例方案有着更好的实用价值，可用来稳定研究复杂的复合体系统（包括蛋白质、核酸和脂质的生物分子系统）的多条转变路径。

经过实验论证，利用本实施例的方案已完成针对5残基的多肽分子（Met-Enkephalin, Met-En）、164残基的T4溶菌酶（T4 Lysozyme L99A，T4L L99A）结构的功能性动力学测试研究。其中，对于Met-En，得到了其由3

进一步，利用本实施例的方案实现了对蛋白质核酸复合体的功能动力学研究，即685残基长度的嗜热古菌Argonaute蛋白与引导DNA单链（21-碱基）的复合体结构，识别载入将被编辑的靶标DNA单链（21-碱基）过程采样，得到了六条代表性的转变路径。

本实施例提供了一种生物分子功能性动力学多转变路径的计算方法，创新性的采用了隐含水的分子动力学模拟，可大幅降低计算量并提高路径搜索效率。基于隐含水的多转变路径，实现显含水转变路径的重建以及优化，并最终得到目标生物分子功能性动力学的转变机制，能够稳定应用于复杂的生物分子功能性动力学机制研究。

实施例2

本发明实施例2公开了一种生物分子功能性动力学多转变路径的计算系统。在实施例1的基础上，将实施例1的方法系统化，具体结构如说明书附图5所示，具体方案如下：

一种生物分子功能性动力学多转变路径的计算系统，包括如下：

初始单元1，用于确定生物分子功能性动力学的起始态结构和目标态结构，并在隐含水条件下，由起始态结构出发得到M个初始结构；

采样单元2，用于进行一轮或多轮采样，每一轮采样包括：从M个初始结构出发，产生预设长度的M条模拟轨迹；随机生成当前概率，通过判断当前概率是否大于预设截断概率以从M条模拟轨迹中挑选M个构象；对挑选的M个构象进行预处理，并将预处理后的结果作为下一轮采样中的M个初始结构，完成一轮采样；

聚类单元3，用于基于每轮采样中的模拟轨迹构建M条转变路径，对M条转变路径进行聚类并选取每一类中的代表转变路径，得到N个不同的代表转变路径；

重建单元4，用于以隐含水条件下的N个代表转变路径为参照，利用预设算法重建显含水条件下的N条转变路径。

优化单元5，用于利用基于旅行商问题的自动路径搜索方法优化显含水条件下的N条转变路径，得到最小自由能转变路径；利用伞形采样方式，计算最小自由能转变路径的自由能面。

其中，采样单元2具体包括：

随机概率模块21，用于随机生成当前概率，判断当前概率是否大于预设截断概率；

第一挑选模块22，用于在当前概率大于截断概率时，从M条模拟轨迹的所有构象中筛选出

第二挑选模块23，用于在当前概率小于等于截断概率时，从M条模拟轨迹的所有构象中筛选出势能大于预设截断势能的高势能构象，并通过预设第一算法从所有高势能构象中挑选M个构象；

预处理模块24，用于对挑选的M个构象进行排序，以使得累加的第i个构象与前一轮采样中第i条模拟轨迹的最后一帧构象之间RMSD距离和为最小；其中，i=1,2,3……，M。

其中，聚类单元3具体包括：

路径转换模块31，用于按照从第一轮到最后一轮的采集顺序，将每轮采样中第i条模拟轨迹整合为一条完整的转变路径，最终得到M条转变路径；

路径聚类模块32，用于计算M条转变路径中任意两条转变路径的相似度，当两条转变路径之间的相似度大于相似度阈值时，此两条路径被归为一类；完成所有转变路径之间的相似度计算，完成对M条转变路径的聚类；

筛选代表模块33，用于对被聚为一类的所有转变路径，通过计算该类中所有转变路径的OM作用量，并筛选出OM作用量最小的转变路径作为该类的代表转变路径。

本实施例提供了一种生物分子功能性动力学多转变路径的计算系统，将实施例1的方法系统化，使其更具实用性。

本发明提供了一种生物分子功能性动力学多转变路径的计算方法及系统，创新性的采用了隐含水的分子动力学模拟，可大幅降低计算量并提高路径搜索效率。基于隐含水的多转变路径，实现显含水转变路径的重建以及优化，并最终得到目标生物分子功能性动力学的转变机制，能够稳定应用于复杂的生物分子功能性动力学机制研究。

本领域技术人员可以理解附图只是一个优选实施场景的示意图，附图中的模块或流程并不一定是实施本发明所必须的。本领域技术人员可以理解实施场景中的装置中的模块可以按照实施场景描述进行分布于实施场景的装置中，也可以进行相应变化位于不同于本实施场景的一个或多个装置中。上述实施场景的模块可以合并为一个模块，也可以进一步拆分成多个子模块。上述本发明序号仅仅为了描述，不代表实施场景的优劣。以上公开的仅为本发明的几个具体实施场景，但是，本发明并非局限于此，任何本领域的技术人员能思之的变化都应落入本发明的保护范围。