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

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

引 言

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

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

本研究選擇分析的蛋白質(zhì)GB1（PDB代碼：6CHE）如圖1所示［12］。GB1是一種小球狀蛋白，由β折疊和α螺旋組成，共有56個殘基。8個殘基與W43形成天然接觸：其中4個殘基（F52、T53、V54和T55）位于相鄰的β折疊中，并與W43的骨架形成天然接觸，而其他4個殘基（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ò)模型中，每個蛋白質(zhì)的三維結(jié)構(gòu)可以簡化為一個彈性網(wǎng)絡(luò)，其中每個氨基酸（殘基）被看作為該網(wǎng)絡(luò)中的頂點，如果兩個頂點間距離小于截止距離，則用一根彈簧將其連接，所有彈簧的彈性系數(shù)都相同［14］?；谠摼W(wǎng)絡(luò)模型，網(wǎng)絡(luò)的總能量可以寫成

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

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

N階對稱矩陣Γ的逆矩陣可表示為

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

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

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

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

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

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

2.2 去折疊原理

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

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

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

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

（2）斷開距離均方漲落值最大的殘基對之間的天然接觸，得到對應(yīng)新Γ矩陣的結(jié)構(gòu)拓撲；

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

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

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

3 結(jié)果與討論

3.1 實驗與模擬所得的B因子分析

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

圖2 實驗與模擬所得的B因子對比Fig.2 Comparison of B factor between experiment and simulation

3.2 蛋白質(zhì)的快運動與慢運動模式分析

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

圖3 GB1快運動模式結(jié)果Fig.3 The fastest mode shapes of GB1

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

圖4 GB1慢運動模式結(jié)果圖Fig.4 The slowest mode shapes of GB1

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

為了詳細說明展開模擬過程中自然接觸的損失，構(gòu)建了不同快照中構(gòu)象的接觸圖，結(jié)果如圖5所示。圖5（a）顯示了GB1蛋白天然結(jié)構(gòu)的接觸圖，即當(dāng)兩個殘基之間的距離小于7.4?時，兩個殘基被定義為相互接觸。如果兩個殘基直接有接觸，則用*表示，圖5（b～f）分別展示了GB1蛋白的非共價接觸損失數(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蛋白的展開有一個優(yōu)先的過程，它顯示了一系列事件。

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

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

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

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

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

4 結(jié)束語

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