生命体系中等位基因竞争性反应QM/MM方法

摘要

本发明涉及一种判断生命体系致病基因反应速率算法，具体为一种生命体系中等位基因竞争性反应QM/MM方法。（1）从蛋白质、核酸等数据库下载相应的分子晶体结构文件，从中选择若个比较合理的分子结构，（2）先用Tinker程序计算MM部分的能量、电荷以及能量梯度，（3）用Gaussian98程序对QM部分进行优化，（4）用tinker中的Newton程序对MM部分进行构型的优化，（5）选取合适的反应坐标，（6）得到了反应路径，进行计算。通过应用本发明中提出的计算方法，可以为诸如镰刀形红细胞贫血症等基因疾病找到合适的治疗方案，即寻找使致病基因不能发生反应或者反应速率慢于正常基因的条件。这对医学是非常有意义的。

说明书

技术领域

本发明涉及一种判断生命体系致病基因反应速率算法，具体为一种生命体系中等位基因竞争性反应QM/MM方法。

背景技术

《哈尔滨师范大学学报》.31卷6期，2015.QM/MM方法的研究。在该文中指出：分子动力学方法是较成熟的模拟方法之一，其可以对由千万个原子构成的生物大分子体系进行高效地模拟计算．但由于其是以经典力学为基础的，无法充分地描述电子的运动，在二十世纪初诞生的量子力学能够充分的地描述电子运动，现今应用量子力学方法已经可以对电子运动进行精确的计算，但是由于其需要进行大量的并且及其复杂的积分运算，其结果就是产生了巨大的运算量，并且严重的增加了计算成本．即使在现今最优秀的计算环境下，也无法很好地进行这种大规模的运算，这种大规模的运算在化学、生物学等多门学科的发展中占有重要的地位，为了可以高效地进行这种大规模的运算，既能得到精确的结果，同时可以有效地降低计算成本，则应运而生了被称作QM／MM方法的混合动力学计算方法，这是采用量子力学(Quantum Mechanics，QM)与分子动力学(Molecular Mechanics，MM)相结合的方法，对生物大分子体系中需要仔细观察的重点部位中所包含的少数的原子使用QM方法进行精确的计算，而体系中剩余的部分采用MM方法进行计算．这种方法提高计算结果的精确性，同时可以有效地降低计算成本，因此正被更多的科研人员所采用。目前，QM／MM方法已经用于溶液体系中的化学变化、生物大分子的处理、酶催化原理等方面，得到了很好的发展。

生命体系中存在着许多化学反应，而且这些化学反应无时无刻不在进行之中，正是由于这些化学反应的进行，生命才得以维持。其中许多反应都是具有竞争性的。例如生命体内酶的催化作用， 有些酶在遇到化学结构与底物相似的分子时，这些分子与底物竞争结合酶的活性中心，就会表现出酶的活性降低（抑制）。抑制剂通过与酶结合而使其无法再与反应底物结合，竞争性抑制剂占据了酶的活性位点而使底物无法进入。竞争的结果就是使反应速率快的化学反应占主导地位，而使体系呈现出该化学反应所表达出的现象。等位基因之间的竞争性表达往往起了决定性的作用。

以镰刀型细胞贫血症为例，镰刀形红细胞贫血(Sickle Cell Anemia SCA)是一种经常出现疼痛现象的严重贫血病。自1910年发现该病，直到1949年Pauling等才证实以是一种血红蛋白异常的“分子病” 。其特点是在氧分压低的情况下，红细胞由正常的双凹圆饼状变成镰刀形。产生镰刀形细胞贫血症的根本原因是基因突变，导致转录、翻译错误,即编码谷氨酸的密码子GAG突变为编码缬氨酸GUG（GAG→GUG），实际上指决定β多肽链那条DNA分子上一碱基发生了改变，即A→T。从而使正常血红蛋白分子（HbA）的β多肽链上的第6位的谷氨酸被缬氨酸代替，造成异常血红蛋白（HbS），导致镰刀状细胞贫血症。

HbA和HbS的氨基酸残基顺序如下：

HbA：缬一组一亮一酪一脯一谷一谷一赖⋯ ⋯

HbS：缬一组一亮一酪一脯一缬一谷一赖⋯ ⋯

谷氨酸的侧链集团带负电荷，缬氨酸的侧链集团不带电荷，为非极性氨基酸残基，见下式：

谷氨酸

缬氨酸

在血红蛋白分子的三维结构中，β链第六位的氨基酸残基位于血红蛋白分子的表面，而非极性缬氨酸残基的存在，使血红蛋白分子的表面出现了一个粘结中心，这样在去氧状态和pH值降低的情况下，血红蛋白另一侧本身存在一个与该区互补的部位。因此，一个Hbs的粒结中心与第二个Hbs的互补部位结合，第二个Hbs的粘结中心与第三个Hbs的互补部位结合，循环往复后形成钢性很强、螺旋管样的长链巨型缔合分子。在Hbs的氧和状态下，互补部位消散或内陷，仅在氧分压低或PH值低的情况下才暴露在分子表面，所以仅在去氧状态下，Hbs才发生线型缔合，撑住红细胞使其镰刀型细胞化。红细胞由于镰刀化而不能正常携带氧气，其寿命又特别短，血中的红细胞数量因此很少，而且极易破裂，从而导致了溶血性贫血症状。

我们猜测，在低氧或pH值低状态下，缬氨酸与生命体内物质发生反应的速率比谷氨酸的反应速率快，导致了红细胞的镰刀形细胞化，从而使生命体呈现出贫血的症状。

我们认为，在生命体系中，发生着多种反应。其中红细胞内存在的缬氨酸与生命体内物质发生反应的速率、谷氨酸与生命体内物质发生反应的速率，前者大于后者，则红细胞形成镰刀型细胞；后者大于前者，形成的红细胞形态正常。

那么我们可以通过对两个反应过程的速率的进行计算比较，就可以预判出生命体系内由于基因变化导致的疾病是否发生。然后我们可以通过一些手段来影响致病基因的反应速录，从而防止致病基因对机体产生破坏。

发明内容

本发明提供了一种生命体系中等位基因竞争性反应QM/MM方法。

本发明的技术方案是，一种生命体系中等位基因竞争性反应QM/MM方法，步骤是：

（1）从蛋白质、核酸等数据库下载相应的分子晶体结构文件，从中选择若个比较合理的分子结构，对每一个分子进行下一步骤计算，

（2）先用Tinker程序计算MM部分的能量、电荷以及能量梯度，

（3）用Gaussian98程序对QM部分进行优化，在这个优化过程中MM部分是固定不动的，并且这个过程中Yps原子用F原子来表示，它的赝势场目前只有3-21G*和6-31G*两个基组下的参数。

（4）用tinker中的Newton程序对MM部分进行构型的优化，同样在这个过程中QM部分是固定不动的，循环进行上述三个过程直到系统总的能量以及RMSD (均方根偏差)均收敛为止，其中RMSD的计算表达式如下:

（5）选取合适的反应坐标，反应坐标可以是两个原子间的距离、两组原子的距离差、三个原子间的夹角以及四个原子间的二面角等参数，反应坐标确定后,就可以对反应路径进行研究，在这个过程中对反应坐标上的每一个构型都要进行上述（2）～（4）步骤的优化。

在优化的过程中为了让反应坐标稳定在预先设定的某个值上，必须在体系的能量函数中引入约束能这一能量项，它的表达式如下：

其中，k为力常数，R为反应坐标，s为反应坐标变化的步长。这样做的目的就是迫使分子体系向预定反应路径的方向进行优化。通过如上的步骤就可以得到整个反应路径的能量-反应坐标曲线。

（6）得到了反应路径，就相当于掌握了反应的整个过程，就可以对这个反应路径进行一系列的计算。

比如可以计算整个反应过程中的自由能垒，反应中的自由能变化，这个可以用自由能微扰来得到。自由能微扰就是将整个反应路径分成很多的点，然后计算相邻两点间的自由能，这样反应路径中任何两点间的自由能变化等于它们之间所有相邻两点间的自由能变化的叠加，对反应路径划分得越细，得到的自由能也越准确。

在步骤（5）中，我们刚好得到了反应路径-反应坐标的曲线，通过对这个曲线进行自由能微扰的计算，就可以得到整个反应过程的自由能变化和自由能垒等重要参数。而这些参数是有实验数据可以考证的，所以这样反过来又可以检验得到的这个反应路径是否合理。在进行自由能微扰的计算过程中值得注意的一点就是，对反应路径上的每一个分子，假定QM部分对整个分子体系自由能的贡献是固定的，它们有前面得到的固定的分子几何构型和电荷分布。

等位基因竞争反应至今尚未在理论上进行研究，本专利申请将弥补理论上的空白，提出对反应过程的合理认识。正是由于等位基因在竞争表达中的化学反应速率不同，才造成了不同的蛋白质表达结果。QM/MM方法具有准确性和高效性，适合处理生物大分子体系，产生重要的理论研究结论，对理论和实验都有非常重要的意义，可以促进学科之间的交叉发展。

通过应用本专利中提出的计算方法，可以为诸如镰刀形红细胞贫血症等基因疾病找到合适的治疗方案，即寻找使致病基因不能发生反应或者反应速率慢于正常基因的条件。这对医学是非常有意义的。

附图说明

图1为本发明QM/MM组合理论中体系的划分示意图。

具体实施方式

使用QM/MM方法，第一步是将整个体系划分割成三个不同的区域 ，即包括反应活性中心的Ⅰ区，反应活性中心周围的Ⅱ区，及所有其他原子构成的Ⅲ区（如图1所示）：

Ⅰ区是考察的关键部位。因为涉及到化学键形成或断裂，这一区域需要用量子力学处理，所以又称 QM区域。Ⅱ区是活性中心之外的所有原子构成的外围环境，不直接涉及化学键的断裂形成，且一般涉及数千或数万个原子，故用分子力学方法来处理，所以又称MM区域。对于不属于Ⅰ区和Ⅱ区的原子则通称为边界区（Ⅲ区）。所以整个体系的有效Hamiltonian 算符可表示为：

(1)

其中等号右侧三项分别是量子体系、量子/分子力学体系及整个分子力学体系的Hamiltonian 算符。如果Ⅰ区包括M 个原子和N个电子，则：

(2)

其中是电子的Laplsoian算符，分别是电子—电子、电子—原子核及核—核之间的距离，Zm是原子核电荷。由于Ⅱ区采用分子力学方法，体系的势能函数仅与原子坐标有关，且包括键伸缩Eb、键弯曲Ea、二面角扭曲Ed和非键作用（库仑Eele、范德华力Evdw）等能量项。

>

>

显然，Ⅰ区和Ⅱ区的相互作用构成了QM/MM组合理论的核心。这种作用主要包括库仑作用，范德华势和极化作用，其Hamiltonian算符可表达为：

>

其中qs是MM区s原子上的部分电荷S则是总的MM作用点，r 和R则代表了QM电子和原子核与MM作用点的距离，是MM 原子s上的诱导偶极矩：

（6）

其中和 各代表原子极化率和原子s感受到外部电场，包括MM和QM两方面的贡献。[44-47]与反应场理论相似，溶剂极化效应的计算须通过洽场方法进行，这一过程计算量很大，故在多数情况下这种极化效应被平均化处理，并直接反应在分子力学参数上。这样，在具体计算时极化作用一般不涉及其余各项[44-46]，体系的总能量为体系有效Hamiltonian对QM波函数的期望值：

（7）

其中前两项包括Ⅰ区中的电子坐标从而须用量子化学方法求得。

Pseudobond QM/MM计算的具体步骤如下：

（1）从蛋白质、核酸等数据库下载相应的分子晶体结构文件，然后通过一系列的能量最小化、平衡步骤，从中选择几个比较合理的分子结构，对每一个分子进行如下的步骤。

（2）先用Tinker程序计算MM部分的能量、电荷以及能量梯度等数据，它们将在下面的计算过程中起到非常重要的作用。

（3）用Gaussian98程序对QM部分进行优化，在这个优化过程中MM部分是固定不动的。并且这个过程中Yps原子用F原子来表示，它的赝势场目前只有3-21G*和6-31G*两个基组下的参数。

（4）用tinker中的Newton程序对MM部分进行构型的优化，同样在这个过程中QM部分是固定不动的。循环进行上述三个过程直到系统总的能量以及RMSD (均方根偏差)均收敛为止，其中RMSD的计算表达式如下:

它表明不同构象的两个分子间的构象差异。至此，对其中L为基组中的最大λ值，ai和bi分别为适当的参数。有效核电势的参数化过程就是让Yps- X键的参数和原来的X- Y键的参数拟合得最好，这样在QM的计算过程中使得Yps和MM部分的作用是相当的。当得到这些参数后再对其他的化学体系进行检验。

值得注意的是有效核电势的参数化过程不是在分子力场的计算中得到的，而是在纯QM计算过程中得到的。这种算法的和其他方法相比，其明显优势就是参数的通用性很强，而且对不同体系计算得到的结果与用纯QM计算得到的结果相近。此构象分子的优化过程结束。

（5）要研究一个反应的机理，也就是反应路径的搜索步骤，必须选取合适的反应坐标。反应坐标可以是两个原子间的距离、两组原子的距离差、三个原子间的夹角以及四个原子间的二面角等参数。一旦反应坐标确定后,就可以对反应路径进行一系列的研究。在这个过程中对反应坐标上的每一个构型都要进行上述（2）～（4）步骤的优化。

在优化的过程中为了让反应坐标稳定在我们预先设定的某个值上，必须在体系的能量函数中引入约束能这一能量项，它的表达式如下：

其中，k为力常数，R为反应坐标，s为反应坐标变化的步长。这样做的目的就是迫使分子体系向预定反应路径的方向进行优化。通过如上的步骤就可以得到整个反应路径的能量-反应坐标曲线。

（6）得到了反应路径，就相当于掌握了反应的整个过程，还可以对这个反应路径进行一系列的计算。比如可以计算整个反应过程中的自由能垒，反应中的自由能变化，这个可以用自由能微扰来得到。自由能微扰就是将整个反应路径分成很多的点，然后计算相邻两点间的自由能，这样反应路径中任何两点间的自由能变化等于它们之间所有相邻两点间的自由能变化的叠加，对反应路径划分得越细，得到的自由能也越准确。

在步骤（5）中，我们刚好得到了反应路径-反应坐标的曲线，通过对这个曲线进行自由能微扰的计算，就可以得到整个反应过程的自由能变化和自由能垒等重要参数。而这些参数是有实验数据可以考证的，所以这样反过来又可以检验得到的这个反应路径是否合理。在进行自由能微扰的计算过程中值得注意的一点就是，对反应路径上的每一个分子，假定QM部分对整个分子体系自由能的贡献是固定的，它们有前面得到的固定的分子几何构型和电荷分布。