張莉莉，蔣益鋒，謝良旭，孔 韌，常 珊

（江蘇理工學(xué)院電氣信息工程學(xué)院生物信息與醫(yī)藥工程研究所，常州 213001）

引 言

在生物體內(nèi)，由蛋白質(zhì)組成的各種類型的分子機(jī)器驅(qū)動(dòng)著各種各樣的生命活動(dòng)所必需的化學(xué)反應(yīng)［1］。深入理解和分析蛋白質(zhì)的結(jié)構(gòu)特征，可以闡明蛋白質(zhì)功能，解釋蛋白質(zhì)錯(cuò)誤折疊引起的相關(guān)疾病起源，以及對(duì)于藥物設(shè)計(jì)工作具有重要意義［2］。蛋白質(zhì)GB1參與許多生理信號(hào)的檢測(cè)，包括激素、神經(jīng)遞質(zhì)和各種感覺(jué)刺激（光、氣味等物質(zhì)）［3］。此外，它還與疾?。ɡ珉貌《竞桶柎暮Ｄ。┫嚓P(guān)的錯(cuò)誤折疊狀態(tài)的存在以及β聚集（淀粉樣疾?。┑难芯肯嚓P(guān)［4］，因此關(guān)于GB1蛋白的研究具有重要意義。理解蛋白質(zhì)GB1結(jié)構(gòu)的折疊機(jī)制和穩(wěn)定性是治療人類疾病的重要基礎(chǔ)，也有助于蛋白質(zhì)的開(kāi)發(fā)設(shè)計(jì)。近年來(lái)，關(guān)于蛋白質(zhì)結(jié)構(gòu)折疊的計(jì)算機(jī)模擬研究方法主要有兩種，即分子動(dòng)力學(xué)模擬和彈性網(wǎng)絡(luò)模型。分子動(dòng)力學(xué)模擬方法是一種細(xì)粒度方法，可以觀察到蛋白質(zhì)的折疊路徑、過(guò)渡態(tài)等，但是該方法計(jì)算復(fù)雜、耗時(shí)較長(zhǎng)，就目前計(jì)算機(jī)的模擬水平，僅僅只對(duì)一些小蛋白質(zhì)分子的折疊結(jié)構(gòu)模擬效果較好［5］。而彈性網(wǎng)絡(luò)模型關(guān)鍵是給出適合簡(jiǎn)化模型的勢(shì)函數(shù)，計(jì)算簡(jiǎn)單、耗時(shí)短，可模擬時(shí)間跨度大的去折疊過(guò)程，相對(duì)分子動(dòng)力學(xué)模擬方法而言效率較高［6-7］，能夠很好地再現(xiàn)蛋白質(zhì)的低頻運(yùn)動(dòng)（長(zhǎng)時(shí)間動(dòng)力學(xué)），提供關(guān)于它們的平衡動(dòng)力學(xué)、天然結(jié)構(gòu)拓?fù)鋵?duì)它們穩(wěn)定性的影響、蛋白質(zhì)波動(dòng)的定位特性或蛋白質(zhì)結(jié)構(gòu)域的定義的信息［8］。彈性網(wǎng)絡(luò)模型（Elastic network model，ENM）在蛋白質(zhì)結(jié)構(gòu)-功能關(guān)系研究中得到了廣泛應(yīng)用。先前有研究通過(guò)應(yīng)用ENM來(lái)評(píng)估生物分子整體編碼、蛋白質(zhì)功能性運(yùn)動(dòng)分析和關(guān)鍵位點(diǎn)識(shí)別等［9］，此外，還有研究結(jié)果表明，應(yīng)用ENM有助于更好地理解轉(zhuǎn)運(yùn)體系發(fā)揮生物學(xué)功能的分子機(jī)制。經(jīng)典的彈性網(wǎng)絡(luò)模型能夠提供蛋白質(zhì)在平衡態(tài)（通常為原生態(tài)）附近的動(dòng)態(tài)特性，因此它們被廣泛應(yīng)用于許多蛋白質(zhì)的系統(tǒng)比較。然而，蛋白質(zhì)折疊通常遠(yuǎn)離平衡態(tài)，所以一般的ENM不適合蛋白質(zhì)折疊的研究。蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)一直是生命科學(xué)里的一個(gè)重要問(wèn)題，研究蛋白質(zhì)序列和結(jié)構(gòu)間關(guān)系的蛋白質(zhì)折疊問(wèn)題是生物物理領(lǐng)域最重要的基礎(chǔ)問(wèn)題之一。在2020年舉辦的第14屆蛋白質(zhì)結(jié)構(gòu)預(yù)測(cè)競(jìng)賽CASP14（Critical assessment of protein structure prediction）中，Google DeepMind團(tuán)隊(duì)使用AlphaFold2預(yù)測(cè)了多個(gè)物種中共30余萬(wàn)個(gè)無(wú)實(shí)驗(yàn)結(jié)構(gòu)的蛋白質(zhì)的結(jié)構(gòu)模型，并聯(lián)手EBI建立了結(jié)構(gòu)預(yù)測(cè)數(shù)據(jù)庫(kù)AFDB［10］。這一系列成果的出現(xiàn)吸引了科學(xué)界的大量關(guān)注。AlphaFold2等結(jié)構(gòu)預(yù)測(cè)方法目前僅能預(yù)測(cè)特定氨基酸序列的靜態(tài)構(gòu)象。蛋白質(zhì)在行使生物學(xué)功能時(shí)往往需要發(fā)生構(gòu)象變化。比如酶從失活狀態(tài)轉(zhuǎn)變?yōu)榛钚誀顟B(tài)、膜轉(zhuǎn)運(yùn)蛋白需要通過(guò)構(gòu)象變化交替接觸膜兩側(cè)的溶液、蛋白和配體結(jié)合時(shí)發(fā)生構(gòu)象變化等等。高斯網(wǎng)絡(luò)模型（Gaussian network model，GNM）是經(jīng)典ENM方法的發(fā)展，是一種基于拓?fù)涞摹⒉灰蕾囆蛄刑禺愋缘拇至６饶Ｐ?。高斯網(wǎng)絡(luò)模型可以從晶體結(jié)構(gòu)提供蛋白質(zhì)構(gòu)象轉(zhuǎn)變的信息，不需要分子動(dòng)力學(xué)模擬的高計(jì)算成本，是一種基于正態(tài)模式計(jì)算的迭代方法，被提出來(lái)用于研究蛋白質(zhì)折疊/去折疊過(guò)程。多年來(lái)，GB1在蛋白質(zhì)折疊的計(jì)算和實(shí)驗(yàn)研究中被廣泛用作模型系統(tǒng)［11］。本文主要就是通過(guò)利用彈性網(wǎng)絡(luò)模型，模擬GB1蛋白結(jié)構(gòu)的展開(kāi)過(guò)程，再現(xiàn)GB1的快運(yùn)動(dòng)與慢運(yùn)動(dòng)模式，同時(shí)研究它的拓?fù)浣Y(jié)構(gòu)對(duì)自身穩(wěn)定性的影響。

1 蛋白質(zhì)結(jié)構(gòu)

本研究選擇分析的蛋白質(zhì)GB1（PDB代碼：6CHE）如圖1所示［12］。GB1是一種小球狀蛋白，由β折疊和α螺旋組成，共有56個(gè)殘基。8個(gè)殘基與W43形成天然接觸：其中4個(gè)殘基（F52、T53、V54和T55）位于相鄰的β折疊中，并與W43的骨架形成天然接觸，而其他4個(gè)殘基（L5、F30、K31和M34）與W43的側(cè)鏈相互 作 用。在GB1中，殘 基2-19形 成N端β折 疊，殘 基23-36形成α螺旋，殘基42-55形成C端β折疊［13］。

2 本文方法

2.1 高斯網(wǎng)絡(luò)模型

在高斯網(wǎng)絡(luò)模型中，每個(gè)蛋白質(zhì)的三維結(jié)構(gòu)可以簡(jiǎn)化為一個(gè)彈性網(wǎng)絡(luò)，其中每個(gè)氨基酸（殘基）被看作為該網(wǎng)絡(luò)中的頂點(diǎn)，如果兩個(gè)頂點(diǎn)間距離小于截止距離，則用一根彈簧將其連接，所有彈簧的彈性系數(shù)都相同［14］?；谠摼W(wǎng)絡(luò)模型，網(wǎng)絡(luò)的總能量可以寫(xiě)成

式中：V為所有接觸殘基的總能量；γ為彈性系數(shù)；{ΔR}為殘基漲落的N維列向量；Γ為N階對(duì)稱矩陣，在對(duì)稱矩陣中的元素可寫(xiě)為

式中：Rij為蛋白質(zhì)中第i個(gè)和第j個(gè)殘基之間的距離；Γc是截止距離（本研究中采用的截止距離為7.4?）。

N階對(duì)稱矩陣Γ的逆矩陣可表示為

式中：U為正交矩陣，其列向量Ui（1≤i≤N）是Γ的特征向量；Λ為對(duì)角矩陣，其對(duì)角線上的元素是Γ的特征值。

蛋白質(zhì)中兩個(gè)殘基均方漲落的互相關(guān)性計(jì)算可表示為

式中：i和j分別表示蛋白質(zhì)中第i個(gè)和第j個(gè)殘基；kB為玻爾茲曼常數(shù)；T為絕對(duì)溫度。當(dāng)i=j時(shí)，第i個(gè)殘基的均方漲落計(jì)算式可表示為

根據(jù)Debye-Waller理論，第i個(gè)殘基的B因子計(jì)算式可表示為

在高斯網(wǎng)絡(luò)模型中，歸一化的互相關(guān)性系數(shù)可寫(xiě)成［15］

高斯網(wǎng)絡(luò)模型是建立在多聚體網(wǎng)絡(luò)的波動(dòng)動(dòng)力學(xué)基礎(chǔ)之上的，彈性網(wǎng)絡(luò)模型可以是原子層次上的粗?；Ｐ停部梢允菤埢鶎哟紊系拇至；Ｐ?。高斯網(wǎng)絡(luò)模型模擬方法可以把蛋白質(zhì)的功能性運(yùn)動(dòng)分解成為各個(gè)不同種運(yùn)動(dòng)模式的疊加，在不同種運(yùn)動(dòng)模式中，慢運(yùn)動(dòng)模式為對(duì)應(yīng)著與蛋白質(zhì)功能相關(guān)的大幅度集合運(yùn)動(dòng)［16］。通過(guò)與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比可以發(fā)現(xiàn)這種方法所得到的數(shù)據(jù)結(jié)果是可靠且有效的。

2.2 去折疊原理

為了研究蛋白質(zhì)的去折疊過(guò)程，本文提出了一種基于高斯網(wǎng)絡(luò)模型的迭代方法。所有殘基對(duì)之間距離的均方漲落都是基于高斯網(wǎng)絡(luò)模型計(jì)算的，第i個(gè)殘基和第j個(gè)殘基之間距離的均方漲落可表示為［17］

式中：Rij和分別為殘基i和殘基j之間的瞬時(shí)和平衡分離向量。

蛋白質(zhì)結(jié)構(gòu)去折疊過(guò)程的模擬方案如下：

（1）基于式（8）和蛋白質(zhì)的天然拓?fù)浣Y(jié)構(gòu)計(jì)算出結(jié)構(gòu)中所有殘基對(duì)之間距離的均方漲落值；

（2）斷開(kāi)距離均方漲落值最大的殘基對(duì)之間的天然接觸，得到對(duì)應(yīng)新Γ矩陣的結(jié)構(gòu)拓?fù)洌?/p>

（3）基于新的Γ矩陣，利用式（8）重新計(jì)算所有殘基對(duì)之間距離的均方漲落值；

（4）重復(fù)上述兩個(gè)步驟，直到蛋白質(zhì)中所有的非共價(jià)接觸被斷開(kāi)；

（5）綜合由以上步驟得到的所有結(jié)構(gòu)拓?fù)湫畔?，以獲取蛋白質(zhì)的去折疊過(guò)程。

3 結(jié)果與討論

3.1 實(shí)驗(yàn)與模擬所得的B因子分析

為了評(píng)價(jià)高斯網(wǎng)絡(luò)模型方法在本研究中應(yīng)用的可行性，計(jì)算了B因子，并與X射線（X-RAY）實(shí)驗(yàn)數(shù)據(jù)對(duì)比。根據(jù)GNM模擬所得的數(shù)據(jù)與X-RAY實(shí)驗(yàn)結(jié)果對(duì)比結(jié)果如圖2所示，其中，紅色曲線對(duì)應(yīng)基于GNM模擬所得的數(shù)據(jù)，綠色曲線對(duì)應(yīng)X-RAY實(shí)驗(yàn)數(shù)據(jù)?？梢钥闯觯瑑蓷l曲線的峰值和谷值出現(xiàn)的位置幾乎相同，模擬所得的數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)之間的相關(guān)系數(shù)為0.70。綜合以往的文獻(xiàn)研究得到，一般模擬數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)之間的相關(guān)系數(shù)取值為0.53～0.89［6，14-15，18］，本 次 實(shí) 驗(yàn)結(jié) 果 所得 值 為0.70，在 此 范圍內(nèi)，可見(jiàn)該方法是適用的，表明該模型適用于研究GB1蛋白的固有動(dòng)力學(xué)。

圖2 實(shí)驗(yàn)與模擬所得的B因子對(duì)比Fig.2 Comparison of B factor between experiment and simulation

3.2 蛋白質(zhì)的快運(yùn)動(dòng)與慢運(yùn)動(dòng)模式分析

運(yùn)動(dòng)的快模式對(duì)應(yīng)于局部結(jié)構(gòu)中的幾何不規(guī)則性。以前的研究發(fā)現(xiàn)，高頻波動(dòng)殘基被認(rèn)為是動(dòng)力學(xué)關(guān)鍵殘基，對(duì)三級(jí)折疊的穩(wěn)定性至關(guān)重要［19］。圖3顯示了GB1蛋白的快運(yùn)動(dòng)模式。圖3中，橫坐標(biāo)表示殘基序號(hào)，縱坐標(biāo)表示殘基自身距離的均方漲落值，基于式（5）求得，單位為平方埃（?2）。從圖3可以看出，殘基Lys4、Ala26、Thr51和Val54（圖中已標(biāo)注）是曲線中的峰值。本文結(jié)果與以前的研究［20］一致，表明這些殘基在蛋白質(zhì)的穩(wěn)定性中起著關(guān)鍵作用。

圖3 GB1快運(yùn)動(dòng)模式結(jié)果Fig.3 The fastest mode shapes of GB1

在蛋白質(zhì)的研究過(guò)程中，慢運(yùn)動(dòng)模式代表著蛋白質(zhì)結(jié)構(gòu)中編碼的長(zhǎng)程集體運(yùn)動(dòng)，同時(shí)相關(guān)研究認(rèn)為那些慢運(yùn)動(dòng)模式就相當(dāng)于大幅度的集體運(yùn)動(dòng)，而大幅度集體運(yùn)動(dòng)往往與蛋白質(zhì)運(yùn)動(dòng)相關(guān)［21］。圖4顯示了基于高斯網(wǎng)絡(luò)模型計(jì)算的GB1蛋白的最慢模式。從圖4可以看出，大多數(shù)殘基波動(dòng)值較高，這意味著這些結(jié)構(gòu)相對(duì)而言不是很穩(wěn)定。同時(shí)，還可以從圖中看出，殘基Gln2、Tyr3和Thr18的波動(dòng)值保持較低。

圖4 GB1慢運(yùn)動(dòng)模式結(jié)果圖Fig.4 The slowest mode shapes of GB1

3.3 蛋白質(zhì)去折疊過(guò)程分析

為了詳細(xì)說(shuō)明展開(kāi)模擬過(guò)程中自然接觸的損失，構(gòu)建了不同快照中構(gòu)象的接觸圖，結(jié)果如圖5所示。圖5（a）顯示了GB1蛋白天然結(jié)構(gòu)的接觸圖，即當(dāng)兩個(gè)殘基之間的距離小于7.4?時(shí)，兩個(gè)殘基被定義為相互接觸。如果兩個(gè)殘基直接有接觸，則用*表示，圖5（b～f）分別展示了GB1蛋白的非共價(jià)接觸損失數(shù)（Loss number of noncovalent contact，LNNC）分別為20、50、100、130和170的結(jié)果。圖5（a）天然狀態(tài)下的接觸呈現(xiàn)結(jié)果與之前的相關(guān)研究一致［22］。結(jié)果表明，GB1蛋白的展開(kāi)有一個(gè)優(yōu)先的過(guò)程，它顯示了一系列事件。

圖5 GB1天然結(jié)構(gòu)以及非共價(jià)接觸損失數(shù)分別為20、50、100、130和170的接觸圖Fig.5 Contact maps of native conformation and conformations with LNNC of 20,50,100,130,170 for GB1

由圖5（a）可以看出，在GB1的天然結(jié)構(gòu)中，碳末端折疊比氮末端折疊有更多更強(qiáng)的接觸（圖1），這可能導(dǎo)致碳末端區(qū)域更快的折疊。從圖5（b，c）的實(shí)驗(yàn)結(jié)果可以看出，隨著殘基對(duì)之間非共價(jià)接觸損失個(gè)數(shù)的增加，GB1蛋白一開(kāi)始主要是從β2折疊部分的殘基對(duì)之間的接觸先斷開(kāi)，此外，從圖5（c）也可以看出β4在非共價(jià)接觸損失數(shù)為50左右的時(shí)候開(kāi)始斷開(kāi)了。繼而如圖5（d，e），α螺旋部分的殘基對(duì)之間的接觸再逐漸斷開(kāi)，直至如圖5（f），最終幾乎所有接觸斷開(kāi)，即GB1蛋白完全展開(kāi)。該過(guò)程顯示了GB1蛋白的展開(kāi)是從大量的α螺旋和β2折疊結(jié)構(gòu)元素的接觸損失開(kāi)始，同時(shí)先保持了大部分其他β結(jié)構(gòu)的完整。本模擬結(jié)果與之前的實(shí)驗(yàn)研究結(jié)果一致［13］。

此外，折疊協(xié)同性被認(rèn)為是蛋白質(zhì)折疊動(dòng)力學(xué)的一個(gè)重要行為［23］。在本研究模型中，展開(kāi)路徑是連續(xù)的，很難直接觀察展開(kāi)過(guò)程中的協(xié)同性。事實(shí)上，這些高度合作的行為發(fā)生在這個(gè)迭代展開(kāi)模型的近鄰步驟中。結(jié)果表明，解折疊路徑主要由其自身的拓?fù)浣Y(jié)構(gòu)決定，迭代解折疊方法可以合理地描述GB1的去折疊過(guò)程。

3.4 蛋白質(zhì)去折疊過(guò)程中的互相關(guān)分析

此外，本文還研究了在GB1蛋白去折疊過(guò)程中殘基波動(dòng)之間的相關(guān)性的變化。殘基波動(dòng)之間的互相關(guān)用式（7）計(jì)算?；ハ嚓P(guān)值的取值范圍為-1到1。其中，正值表示殘基間運(yùn)動(dòng)方向相同，負(fù)值則表示它們之間運(yùn)動(dòng)方向相反。絕對(duì)互相關(guān)值越高，兩個(gè)殘基越相關(guān)（或反相關(guān)）。另外，互相關(guān)值0意味著殘基的運(yùn)動(dòng)完全不相關(guān)［14］。圖6顯示了GB1蛋白的互相關(guān)圖。

圖6 GB1天然結(jié)構(gòu)以及展開(kāi)過(guò)程中非共價(jià)接觸損失數(shù)分別為20、50、100、130和170時(shí)的殘基互相關(guān)圖Fig.6 Cross-correlation maps calculated using all modes for native conformation and conformations with LNNC of 20,50,100,130,170 during the unfolding process of GB1

如圖6（a）所示，沿著圖的對(duì)角線，有一些正相關(guān)的光塊，對(duì)應(yīng)α螺旋和β折疊的二級(jí)結(jié)構(gòu)。隨著殘基對(duì)之間非共價(jià)損失個(gè)數(shù)的增加，即隨著GB1蛋白的逐漸展開(kāi)，如圖6（b，c），當(dāng)α螺旋和β折疊中的天然觸點(diǎn)開(kāi)始丟失時(shí)，α螺旋和β折疊之間負(fù)相關(guān)，β折疊之間的正相關(guān)性提高；隨著天然觸點(diǎn)丟失個(gè)數(shù)的增加，如圖6（d，e），α螺旋和β折疊僅部分保留，最后，如圖6（f），蛋白質(zhì)的結(jié)構(gòu)似乎被分成兩個(gè)方向相反的方向波動(dòng)。該圖反映的是去折疊的最后狀態(tài)，即蛋白質(zhì)結(jié)構(gòu)展開(kāi)回到了最初未折疊的多肽鏈結(jié)構(gòu)。根據(jù)先前的研究發(fā)現(xiàn)［14］，當(dāng)去折疊模擬到最后的階段時(shí)，蛋白質(zhì)的結(jié)構(gòu)也似乎被分為兩部分，上下波動(dòng)方向相反，與本次實(shí)驗(yàn)結(jié)果一致。

4 結(jié)束語(yǔ)

本研究基于GB1的拓?fù)浣Y(jié)構(gòu)，采用高斯網(wǎng)絡(luò)模型模擬了GB1的快運(yùn)動(dòng)與慢運(yùn)動(dòng)模式，并對(duì)其做了相應(yīng)的結(jié)果分析；同時(shí)，對(duì)其拓?fù)浣Y(jié)構(gòu)做了展開(kāi)過(guò)程的路徑研究；此外，還研究了GB1蛋白在去折疊過(guò)程中殘基波動(dòng)之間相關(guān)性的變化。與相關(guān)實(shí)驗(yàn)和分子動(dòng)力學(xué)模擬數(shù)據(jù)吻合良好，表明彈性網(wǎng)絡(luò)模型的計(jì)算效率高，能夠準(zhǔn)確模擬蛋白質(zhì)的動(dòng)態(tài)和結(jié)構(gòu)特性，能夠很好地再現(xiàn)蛋白質(zhì)的運(yùn)動(dòng)特性，提供關(guān)于它們的平衡動(dòng)力學(xué)、天然結(jié)構(gòu)拓?fù)鋵?duì)其穩(wěn)定性的影響、蛋白質(zhì)波動(dòng)的定位特性或蛋白質(zhì)結(jié)構(gòu)域的信息，適用于對(duì)蛋白質(zhì)的研究。