本發(fā)明涉及一種判斷生命體系致病基因反應(yīng)速率算法，具體為一種生命體系中等位基因競爭性反應(yīng)QM/MM方法。

背景技術(shù)：

《哈爾濱師范大學(xué)學(xué)報》.31卷6期，2015.QM/MM方法的研究。在該文中指出：分子動力學(xué)方法是較成熟的模擬方法之一，其可以對由千萬個原子構(gòu)成的生物大分子體系進(jìn)行高效地模擬計算．但由于其是以經(jīng)典力學(xué)為基礎(chǔ)的，無法充分地描述電子的運動，在二十世紀(jì)初誕生的量子力學(xué)能夠充分的地描述電子運動，現(xiàn)今應(yīng)用量子力學(xué)方法已經(jīng)可以對電子運動進(jìn)行精確的計算，但是由于其需要進(jìn)行大量的并且及其復(fù)雜的積分運算，其結(jié)果就是產(chǎn)生了巨大的運算量，并且嚴(yán)重的增加了計算成本．即使在現(xiàn)今最優(yōu)秀的計算環(huán)境下，也無法很好地進(jìn)行這種大規(guī)模的運算，這種大規(guī)模的運算在化學(xué)、生物學(xué)等多門學(xué)科的發(fā)展中占有重要的地位，為了可以高效地進(jìn)行這種大規(guī)模的運算，既能得到精確的結(jié)果，同時可以有效地降低計算成本，則應(yīng)運而生了被稱作QM／MM方法的混合動力學(xué)計算方法，這是采用量子力學(xué)(Quantum Mechanics，QM)與分子動力學(xué)(Molecular Mechanics，MM)相結(jié)合的方法，對生物大分子體系中需要仔細(xì)觀察的重點部位中所包含的少數(shù)的原子使用QM方法進(jìn)行精確的計算，而體系中剩余的部分采用MM方法進(jìn)行計算．這種方法提高計算結(jié)果的精確性，同時可以有效地降低計算成本，因此正被更多的科研人員所采用。目前，QM／MM方法已經(jīng)用于溶液體系中的化學(xué)變化、生物大分子的處理、酶催化原理等方面，得到了很好的發(fā)展。

生命體系中存在著許多化學(xué)反應(yīng)，而且這些化學(xué)反應(yīng)無時無刻不在進(jìn)行之中，正是由于這些化學(xué)反應(yīng)的進(jìn)行，生命才得以維持。其中許多反應(yīng)都是具有競爭性的。例如生命體內(nèi)酶的催化作用， 有些酶在遇到化學(xué)結(jié)構(gòu)與底物相似的分子時，這些分子與底物競爭結(jié)合酶的活性中心，就會表現(xiàn)出酶的活性降低（抑制）。抑制劑通過與酶結(jié)合而使其無法再與反應(yīng)底物結(jié)合，競爭性抑制劑占據(jù)了酶的活性位點而使底物無法進(jìn)入。競爭的結(jié)果就是使反應(yīng)速率快的化學(xué)反應(yīng)占主導(dǎo)地位，而使體系呈現(xiàn)出該化學(xué)反應(yīng)所表達(dá)出的現(xiàn)象。等位基因之間的競爭性表達(dá)往往起了決定性的作用。

以鐮刀型細(xì)胞貧血癥為例，鐮刀形紅細(xì)胞貧血(Sickle Cell Anemia SCA)是一種經(jīng)常出現(xiàn)疼痛現(xiàn)象的嚴(yán)重貧血病。自1910年發(fā)現(xiàn)該病，直到1949年P(guān)auling等才證實以是一種血紅蛋白異常的“分子病” 。其特點是在氧分壓低的情況下，紅細(xì)胞由正常的雙凹圓餅狀變成鐮刀形。產(chǎn)生鐮刀形細(xì)胞貧血癥的根本原因是基因突變，導(dǎo)致轉(zhuǎn)錄、翻譯錯誤,即編碼谷氨酸的密碼子GAG突變?yōu)榫幋a纈氨酸GUG（GAG→GUG），實際上指決定β多肽鏈那條DNA分子上一堿基發(fā)生了改變，即A→T。從而使正常血紅蛋白分子（HbA）的β多肽鏈上的第6位的谷氨酸被纈氨酸代替，造成異常血紅蛋白（HbS），導(dǎo)致鐮刀狀細(xì)胞貧血癥。

HbA和HbS的氨基酸殘基順序如下：

HbA：纈一組一亮一酪一脯一谷一谷一賴? ?

HbS：纈一組一亮一酪一脯一纈一谷一賴? ?

谷氨酸的側(cè)鏈集團(tuán)帶負(fù)電荷，纈氨酸的側(cè)鏈集團(tuán)不帶電荷，為非極性氨基酸殘基，見下式：

谷氨酸

纈氨酸

在血紅蛋白分子的三維結(jié)構(gòu)中，β鏈第六位的氨基酸殘基位于血紅蛋白分子的表面，而非極性纈氨酸殘基的存在，使血紅蛋白分子的表面出現(xiàn)了一個粘結(jié)中心，這樣在去氧狀態(tài)和pH值降低的情況下，血紅蛋白另一側(cè)本身存在一個與該區(qū)互補的部位。因此，一個Hbs的粒結(jié)中心與第二個Hbs的互補部位結(jié)合，第二個Hbs的粘結(jié)中心與第三個Hbs的互補部位結(jié)合，循環(huán)往復(fù)后形成鋼性很強、螺旋管樣的長鏈巨型締合分子。在Hbs的氧和狀態(tài)下，互補部位消散或內(nèi)陷，僅在氧分壓低或PH值低的情況下才暴露在分子表面，所以僅在去氧狀態(tài)下，Hbs才發(fā)生線型締合，撐住紅細(xì)胞使其鐮刀型細(xì)胞化。紅細(xì)胞由于鐮刀化而不能正常攜帶氧氣，其壽命又特別短，血中的紅細(xì)胞數(shù)量因此很少，而且極易破裂，從而導(dǎo)致了溶血性貧血癥狀。

我們猜測，在低氧或pH值低狀態(tài)下，纈氨酸與生命體內(nèi)物質(zhì)發(fā)生反應(yīng)的速率比谷氨酸的反應(yīng)速率快，導(dǎo)致了紅細(xì)胞的鐮刀形細(xì)胞化，從而使生命體呈現(xiàn)出貧血的癥狀。

我們認(rèn)為，在生命體系中，發(fā)生著多種反應(yīng)。其中紅細(xì)胞內(nèi)存在的纈氨酸與生命體內(nèi)物質(zhì)發(fā)生反應(yīng)的速率、谷氨酸與生命體內(nèi)物質(zhì)發(fā)生反應(yīng)的速率，前者大于后者，則紅細(xì)胞形成鐮刀型細(xì)胞；后者大于前者，形成的紅細(xì)胞形態(tài)正常。

那么我們可以通過對兩個反應(yīng)過程的速率的進(jìn)行計算比較，就可以預(yù)判出生命體系內(nèi)由于基因變化導(dǎo)致的疾病是否發(fā)生。然后我們可以通過一些手段來影響致病基因的反應(yīng)速錄，從而防止致病基因?qū)C體產(chǎn)生破壞。

技術(shù)實現(xiàn)要素：

本發(fā)明提供了一種生命體系中等位基因競爭性反應(yīng)QM/MM方法。

本發(fā)明的技術(shù)方案是，一種生命體系中等位基因競爭性反應(yīng)QM/MM方法，步驟是：

（1）從蛋白質(zhì)、核酸等數(shù)據(jù)庫下載相應(yīng)的分子晶體結(jié)構(gòu)文件，從中選擇若個比較合理的分子結(jié)構(gòu)，對每一個分子進(jìn)行下一步驟計算，

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

（3）用Gaussian98程序?qū)M部分進(jìn)行優(yōu)化，在這個優(yōu)化過程中MM部分是固定不動的，并且這個過程中Yps原子用F原子來表示，它的贗勢場目前只有3-21G*和6-31G*兩個基組下的參數(shù)。

（4）用tinker中的Newton程序?qū)M部分進(jìn)行構(gòu)型的優(yōu)化，同樣在這個過程中QM部分是固定不動的，循環(huán)進(jìn)行上述三個過程直到系統(tǒng)總的能量以及RMSD (均方根偏差)均收斂為止，其中RMSD的計算表達(dá)式如下:

（5）選取合適的反應(yīng)坐標(biāo)，反應(yīng)坐標(biāo)可以是兩個原子間的距離、兩組原子的距離差、三個原子間的夾角以及四個原子間的二面角等參數(shù)，反應(yīng)坐標(biāo)確定后,就可以對反應(yīng)路徑進(jìn)行研究，在這個過程中對反應(yīng)坐標(biāo)上的每一個構(gòu)型都要進(jìn)行上述（2）～（4）步驟的優(yōu)化。

在優(yōu)化的過程中為了讓反應(yīng)坐標(biāo)穩(wěn)定在預(yù)先設(shè)定的某個值上，必須在體系的能量函數(shù)中引入約束能這一能量項，它的表達(dá)式如下：

其中，k為力常數(shù)，R為反應(yīng)坐標(biāo)，s為反應(yīng)坐標(biāo)變化的步長。這樣做的目的就是迫使分子體系向預(yù)定反應(yīng)路徑的方向進(jìn)行優(yōu)化。通過如上的步驟就可以得到整個反應(yīng)路徑的能量-反應(yīng)坐標(biāo)曲線。

（6）得到了反應(yīng)路徑，就相當(dāng)于掌握了反應(yīng)的整個過程，就可以對這個反應(yīng)路徑進(jìn)行一系列的計算。

比如可以計算整個反應(yīng)過程中的自由能壘，反應(yīng)中的自由能變化，這個可以用自由能微擾來得到。自由能微擾就是將整個反應(yīng)路徑分成很多的點，然后計算相鄰兩點間的自由能，這樣反應(yīng)路徑中任何兩點間的自由能變化等于它們之間所有相鄰兩點間的自由能變化的疊加，對反應(yīng)路徑劃分得越細(xì)，得到的自由能也越準(zhǔn)確。

在步驟（5）中，我們剛好得到了反應(yīng)路徑-反應(yīng)坐標(biāo)的曲線，通過對這個曲線進(jìn)行自由能微擾的計算，就可以得到整個反應(yīng)過程的自由能變化和自由能壘等重要參數(shù)。而這些參數(shù)是有實驗數(shù)據(jù)可以考證的，所以這樣反過來又可以檢驗得到的這個反應(yīng)路徑是否合理。在進(jìn)行自由能微擾的計算過程中值得注意的一點就是，對反應(yīng)路徑上的每一個分子，假定QM部分對整個分子體系自由能的貢獻(xiàn)是固定的，它們有前面得到的固定的分子幾何構(gòu)型和電荷分布。

等位基因競爭反應(yīng)至今尚未在理論上進(jìn)行研究，本專利申請將彌補理論上的空白，提出對反應(yīng)過程的合理認(rèn)識。正是由于等位基因在競爭表達(dá)中的化學(xué)反應(yīng)速率不同，才造成了不同的蛋白質(zhì)表達(dá)結(jié)果。QM/MM方法具有準(zhǔn)確性和高效性，適合處理生物大分子體系，產(chǎn)生重要的理論研究結(jié)論，對理論和實驗都有非常重要的意義，可以促進(jìn)學(xué)科之間的交叉發(fā)展。

通過應(yīng)用本專利中提出的計算方法，可以為諸如鐮刀形紅細(xì)胞貧血癥等基因疾病找到合適的治療方案，即尋找使致病基因不能發(fā)生反應(yīng)或者反應(yīng)速率慢于正常基因的條件。這對醫(yī)學(xué)是非常有意義的。

附圖說明

圖1為本發(fā)明QM/MM組合理論中體系的劃分示意圖。

具體實施方式

使用QM/MM方法，第一步是將整個體系劃分割成三個不同的區(qū)域 ，即包括反應(yīng)活性中心的Ⅰ區(qū)，反應(yīng)活性中心周圍的Ⅱ區(qū)，及所有其他原子構(gòu)成的Ⅲ區(qū)（如圖1所示）：

Ⅰ區(qū)是考察的關(guān)鍵部位。因為涉及到化學(xué)鍵形成或斷裂，這一區(qū)域需要用量子力學(xué)處理，所以又稱 QM區(qū)域。Ⅱ區(qū)是活性中心之外的所有原子構(gòu)成的外圍環(huán)境，不直接涉及化學(xué)鍵的斷裂形成，且一般涉及數(shù)千或數(shù)萬個原子，故用分子力學(xué)方法來處理，所以又稱MM區(qū)域。對于不屬于Ⅰ區(qū)和Ⅱ區(qū)的原子則通稱為邊界區(qū)（Ⅲ區(qū)）。所以整個體系的有效Hamiltonian 算符可表示為：

(1)

其中等號右側(cè)三項分別是量子體系、量子/分子力學(xué)體系及整個分子力學(xué)體系的Hamiltonian 算符。如果Ⅰ區(qū)包括M 個原子和N個電子，則：

(2)

其中是電子的Laplsoian算符，分別是電子—電子、電子—原子核及核—核之間的距離，Zm是原子核電荷。由于Ⅱ區(qū)采用分子力學(xué)方法，體系的勢能函數(shù)僅與原子坐標(biāo)有關(guān)，且包括鍵伸縮Eb、鍵彎曲Ea、二面角扭曲Ed和非鍵作用（庫侖Eele、范德華力Evdw）等能量項。

（3）

（4）

顯然，Ⅰ區(qū)和Ⅱ區(qū)的相互作用構(gòu)成了QM/MM組合理論的核心。這種作用主要包括庫侖作用，范德華勢和極化作用，其Hamiltonian算符可表達(dá)為：

(5)

其中qs是MM區(qū)s原子上的部分電荷S則是總的MM作用點，r 和R則代表了QM電子和原子核與MM作用點的距離，是MM 原子s上的誘導(dǎo)偶極矩：

（6）

其中和 各代表原子極化率和原子s感受到外部電場，包括MM和QM兩方面的貢獻(xiàn)。[44-47]與反應(yīng)場理論相似，溶劑極化效應(yīng)的計算須通過洽場方法進(jìn)行，這一過程計算量很大，故在多數(shù)情況下這種極化效應(yīng)被平均化處理，并直接反應(yīng)在分子力學(xué)參數(shù)上。這樣，在具體計算時極化作用一般不涉及其余各項[44-46]，體系的總能量為體系有效Hamiltonian對QM波函數(shù)的期望值：

（7）

其中前兩項包括Ⅰ區(qū)中的電子坐標(biāo)從而須用量子化學(xué)方法求得。

Pseudobond QM/MM計算的具體步驟如下：

（1）從蛋白質(zhì)、核酸等數(shù)據(jù)庫下載相應(yīng)的分子晶體結(jié)構(gòu)文件，然后通過一系列的能量最小化、平衡步驟，從中選擇幾個比較合理的分子結(jié)構(gòu)，對每一個分子進(jìn)行如下的步驟。

（2）先用Tinker程序計算MM部分的能量、電荷以及能量梯度等數(shù)據(jù)，它們將在下面的計算過程中起到非常重要的作用。

（3）用Gaussian98程序?qū)M部分進(jìn)行優(yōu)化，在這個優(yōu)化過程中MM部分是固定不動的。并且這個過程中Yps原子用F原子來表示，它的贗勢場目前只有3-21G*和6-31G*兩個基組下的參數(shù)。

（4）用tinker中的Newton程序?qū)M部分進(jìn)行構(gòu)型的優(yōu)化，同樣在這個過程中QM部分是固定不動的。循環(huán)進(jìn)行上述三個過程直到系統(tǒng)總的能量以及RMSD (均方根偏差)均收斂為止，其中RMSD的計算表達(dá)式如下:

它表明不同構(gòu)象的兩個分子間的構(gòu)象差異。至此，對其中L為基組中的最大λ值，ai和bi分別為適當(dāng)?shù)膮?shù)。有效核電勢的參數(shù)化過程就是讓Yps- X鍵的參數(shù)和原來的X- Y鍵的參數(shù)擬合得最好，這樣在QM的計算過程中使得Yps和MM部分的作用是相當(dāng)?shù)?。?dāng)?shù)玫竭@些參數(shù)后再對其他的化學(xué)體系進(jìn)行檢驗。

值得注意的是有效核電勢的參數(shù)化過程不是在分子力場的計算中得到的，而是在純QM計算過程中得到的。這種算法的和其他方法相比，其明顯優(yōu)勢就是參數(shù)的通用性很強，而且對不同體系計算得到的結(jié)果與用純QM計算得到的結(jié)果相近。此構(gòu)象分子的優(yōu)化過程結(jié)束。

（5）要研究一個反應(yīng)的機理，也就是反應(yīng)路徑的搜索步驟，必須選取合適的反應(yīng)坐標(biāo)。反應(yīng)坐標(biāo)可以是兩個原子間的距離、兩組原子的距離差、三個原子間的夾角以及四個原子間的二面角等參數(shù)。一旦反應(yīng)坐標(biāo)確定后,就可以對反應(yīng)路徑進(jìn)行一系列的研究。在這個過程中對反應(yīng)坐標(biāo)上的每一個構(gòu)型都要進(jìn)行上述（2）～（4）步驟的優(yōu)化。

在優(yōu)化的過程中為了讓反應(yīng)坐標(biāo)穩(wěn)定在我們預(yù)先設(shè)定的某個值上，必須在體系的能量函數(shù)中引入約束能這一能量項，它的表達(dá)式如下：

其中，k為力常數(shù)，R為反應(yīng)坐標(biāo)，s為反應(yīng)坐標(biāo)變化的步長。這樣做的目的就是迫使分子體系向預(yù)定反應(yīng)路徑的方向進(jìn)行優(yōu)化。通過如上的步驟就可以得到整個反應(yīng)路徑的能量-反應(yīng)坐標(biāo)曲線。

（6）得到了反應(yīng)路徑，就相當(dāng)于掌握了反應(yīng)的整個過程，還可以對這個反應(yīng)路徑進(jìn)行一系列的計算。比如可以計算整個反應(yīng)過程中的自由能壘，反應(yīng)中的自由能變化，這個可以用自由能微擾來得到。自由能微擾就是將整個反應(yīng)路徑分成很多的點，然后計算相鄰兩點間的自由能，這樣反應(yīng)路徑中任何兩點間的自由能變化等于它們之間所有相鄰兩點間的自由能變化的疊加，對反應(yīng)路徑劃分得越細(xì)，得到的自由能也越準(zhǔn)確。

在步驟（5）中，我們剛好得到了反應(yīng)路徑-反應(yīng)坐標(biāo)的曲線，通過對這個曲線進(jìn)行自由能微擾的計算，就可以得到整個反應(yīng)過程的自由能變化和自由能壘等重要參數(shù)。而這些參數(shù)是有實驗數(shù)據(jù)可以考證的，所以這樣反過來又可以檢驗得到的這個反應(yīng)路徑是否合理。在進(jìn)行自由能微擾的計算過程中值得注意的一點就是，對反應(yīng)路徑上的每一個分子，假定QM部分對整個分子體系自由能的貢獻(xiàn)是固定的，它們有前面得到的固定的分子幾何構(gòu)型和電荷分布。