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

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

引 言

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

1 蛋白質結構

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

2 本文方法

2.1 高斯網絡模型

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

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

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

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

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

蛋白質中兩個殘基均方漲落的互相關性計算可表示為

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

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

在高斯網絡模型中，歸一化的互相關性系數(shù)可寫成［15］

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

2.2 去折疊原理

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

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

蛋白質結構去折疊過程的模擬方案如下：

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

（2）斷開距離均方漲落值最大的殘基對之間的天然接觸，得到對應新Γ矩陣的結構拓撲；

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

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

（5）綜合由以上步驟得到的所有結構拓撲信息，以獲取蛋白質的去折疊過程。

3 結果與討論

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

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

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

3.2 蛋白質的快運動與慢運動模式分析

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

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

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

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

3.3 蛋白質去折疊過程分析

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

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

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

3.4 蛋白質去折疊過程中的互相關分析

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

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

4 結束語

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