基于彈性網(wǎng)絡(luò)模型的蛋白質(zhì)變構(gòu)路徑與關(guān)鍵殘基識別研究

李 嬌 王韋卜 蘇計國

（1）燕山大學(xué)理學(xué)院，河北省微結(jié)構(gòu)材料物理重點實驗室，秦皇島 066004；2）國藥中生生物技術(shù)研究院有限公司第六研究室，北京 101111）

變構(gòu)調(diào)節(jié)是蛋白質(zhì)功能調(diào)節(jié)的一種重要機制，變構(gòu)效應(yīng)因子在變構(gòu)位點的結(jié)合，可以通過特定的變構(gòu)信號傳輸路徑傳遞到遠(yuǎn)處的活性位點，使活性部位的局部結(jié)構(gòu)或構(gòu)象運動發(fā)生改變，進而對蛋白質(zhì)的功能進行調(diào)節(jié)［1-2］。變構(gòu)相互作用最早由Umbarger等［3］在1956年提出，他們發(fā)現(xiàn)了左旋異亮氨酸對左旋蘇氨酸脫氨酶的反饋抑制作用。之后，Changeux等［4-6］在含有幾個亞單位的蛋白質(zhì)體系中，觀察到配體結(jié)合引起蛋白質(zhì)結(jié)構(gòu)發(fā)生一定的構(gòu)象變化，由此提出了變構(gòu)效應(yīng)。盡管蛋白質(zhì)的變構(gòu)調(diào)節(jié)在結(jié)構(gòu)生物學(xué)領(lǐng)域的研究進行了大半個世紀(jì)，但是蛋白質(zhì)變構(gòu)調(diào)節(jié)的物理機制以及相應(yīng)變構(gòu)路徑的有效識別，仍是一個有待深入研究的重要科學(xué)問題。變構(gòu)調(diào)節(jié)在蛋白質(zhì)功能反饋和調(diào)控機制中起著關(guān)鍵性作用，許多蛋白質(zhì)的活性及其功能的調(diào)節(jié)都是通過變構(gòu)效應(yīng)來實現(xiàn)的。因此，蛋白質(zhì)變構(gòu)路徑的有效識別，不僅有助于揭示蛋白質(zhì)功能調(diào)控的分子機制，還可以為變構(gòu)藥物的開發(fā)或納米生物傳感器件［7-9］的開發(fā)提供理論依據(jù)。

分子動力學(xué)模擬（molecular dynamics，MD）已廣泛應(yīng)用于蛋白質(zhì)的變構(gòu)效應(yīng)研究中，通過對構(gòu)象運動和運動關(guān)聯(lián)性的分析等來揭示蛋白質(zhì)體系的變構(gòu)路徑和關(guān)鍵殘基［10-15］。Roy 等［16］運用MD 模擬結(jié)合運動相關(guān)性分析，對人類鼻病毒（HRⅤ）的變構(gòu)調(diào)控過程進行了研究，他們通過對不同殘基之間徑向運動關(guān)聯(lián)性的分析，揭示了抗病毒化合物WⅠN 52084 對HRⅤ衣殼蛋白“呼吸”運動的長程變構(gòu)調(diào)控機制。Bowerman 和Wereszczynski［17］基于MD模擬軌跡，利用殘基運動相關(guān)性分析、互信息關(guān)聯(lián)分析以及圖論的方法來識別蛋白質(zhì)體系的變構(gòu)網(wǎng)絡(luò)，利用該方法他們詳細(xì)分析了凝血酶變構(gòu)位點與催化中心之間的長程變構(gòu)信息傳播路徑。MD模擬雖然在蛋白質(zhì)變構(gòu)效應(yīng)研究中得到了成功應(yīng)用，但是，很多情況下，蛋白質(zhì)的變構(gòu)調(diào)控并不一定伴隨明顯的構(gòu)象運動，而是通過剛性區(qū)域來耗散變構(gòu)所產(chǎn)生的應(yīng)力，這些剛性區(qū)域不會發(fā)生大的構(gòu)象變化，而是產(chǎn)生較大的力信號［18］，常規(guī)的平衡態(tài)MD模擬很難識別到這些力信號。因此，一些研究者通過對體系施加外界的擾動，利用MD模擬結(jié)合力分布分析（force distribution analysis，F(xiàn)DA）方法來揭示蛋白質(zhì)內(nèi)力信號的長程傳播過程，進而獲得體系的變構(gòu)調(diào)控路徑［19-23］。Stacklies 等［24］對肌聯(lián)蛋白Ⅰ27結(jié)構(gòu)域的長程力學(xué)變構(gòu)調(diào)控進行了研究，對Ⅰ27的N端、C端施加300 pN的恒力進行拉伸，運用MD 結(jié)合FDA 的方法，計算了不同殘基間的力分布，得到了蛋白質(zhì)力信號的傳導(dǎo)路徑。隨后，該小組運用同樣的方法，對MetJ 蛋白進行了研究，識別了力信號從變構(gòu)位點傳遞到活性位點的變構(gòu)路徑［18］。2013年，Palmai等［25］利用MD結(jié)合FDA 的方法，研究了3-磷酸還原酶（hPGK）催化功能執(zhí)行過程中的變構(gòu)調(diào)控機制，揭示了變構(gòu)信號從底物結(jié)合位點到功能性運動鉸鏈區(qū)的傳導(dǎo)路徑。這些研究表明，基于MD模擬軌跡，通過對蛋白質(zhì)體系內(nèi)部力分布的分析來研究蛋白質(zhì)變構(gòu)調(diào)控路徑是一種有效的分析方法。但是，利用MD模擬來計算蛋白質(zhì)內(nèi)的力分布存在一定的局限性。首先，由于蛋白質(zhì)體系的熱漲落，力學(xué)信號常常被熱噪聲所淹沒，需要設(shè)計一定的去噪音手段才能獲取有價值的力學(xué)信號。其次，為了去除熱噪聲，獲得具有統(tǒng)計意義的力分布結(jié)果，需要進行大量軌跡的模擬，計算非常耗時。而彈性網(wǎng)絡(luò)模型（elastic network model，ENM）不需要考慮熱噪聲的影響，并且對于體系的配分函數(shù)和熱力學(xué)性質(zhì)可以進行解析求解，大大的減少計算量，提高計算效率。

ENM 方法是研究蛋白質(zhì)結(jié)構(gòu)和動力學(xué)性質(zhì)的有效方法［26-27］，在蛋白質(zhì)結(jié)構(gòu)-功能關(guān)系研究中得到了廣泛應(yīng)用，包括蛋白質(zhì)結(jié)構(gòu)穩(wěn)定性的分析、關(guān)鍵殘基的識別、蛋白質(zhì)結(jié)構(gòu)域的分解、大幅度功能性構(gòu)象運動的分析等［28-35］。Wang等［36］基于ENM，通過力分布的計算來分析蛋白質(zhì)體系對外力的響應(yīng)過程，研究了GB1 蛋白質(zhì)力學(xué)性能的變構(gòu)調(diào)控機制。對GB1蛋白的N端、C端施加20 pN的外力進行拉伸，分別研究結(jié)合配體與未結(jié)合配體情況下蛋白質(zhì)體系的形變情況和力分布情況。研究發(fā)現(xiàn)，配體的結(jié)合能夠通過長程變構(gòu)效應(yīng)調(diào)控抗力部位的機械性能，顯著增強GB1的力學(xué)穩(wěn)定性。

本工作將基于ENM 的力分布計算方法用于蛋白質(zhì)變構(gòu)路徑和關(guān)鍵殘基位點的識別研究，基于該方法，分析了人類磷酸甘油酸激酶（human phosphoglycerate kinase，hPGK）和蛋白質(zhì)酪氨酸磷酸酶（protein tyrosine phosphatase，PTP）PDZ2結(jié)構(gòu)域的變構(gòu)路徑和關(guān)鍵殘基。對兩種蛋白質(zhì)的變構(gòu)位點施加外力進行拉伸，通過計算外力作用下蛋白質(zhì)內(nèi)部的力分布，來識別與力承載區(qū)域形變相耦合的殘基位點，并研究力信號在蛋白質(zhì)結(jié)構(gòu)內(nèi)的傳播路徑。

1 理論與方法

ENM 將蛋白質(zhì)的空間結(jié)構(gòu)簡化為一個由大量節(jié)點和節(jié)點之間的彈簧所構(gòu)成的彈性網(wǎng)絡(luò)。將蛋白質(zhì)體系中的每一個殘基簡化為網(wǎng)絡(luò)中的一個節(jié)點，用Cα 原子來代替，其他的原子忽略。選取一個幾何距離作為截斷半徑rc［26］。如果兩個殘基間的距離小于該截斷半徑，認(rèn)為它們之間存在相互作用，則用彈簧連接；若兩殘基間距離大于該截斷半徑，則認(rèn)為不存在相互作用。所有連接的彈簧的彈性系數(shù)均相同。在本工作中，彈簧的彈性系數(shù)設(shè)為1.0 N/m，選取的截斷半徑為8.0 ?。由整個體系中彈簧的連接情況，可得到體系的勢能為

其中，γ為彈簧的彈性系數(shù)；ΔR為蛋白質(zhì)內(nèi)各個氨基酸偏離其平衡位置的位移，是3N維列向量，N為殘基的個數(shù)；右上角的T表示向量的轉(zhuǎn)置；H為3N×3N的Hessian矩陣，該矩陣可表示為

上式中的hij為3×3的子矩陣，i、j分別表示第i、j個殘基。

當(dāng)i≠j時，有

當(dāng)i=j時，有

其中，矩陣元為勢能函數(shù)V對坐標(biāo)x、y求二階偏導(dǎo)，表達(dá)式為

在本工作中，由于對蛋白質(zhì)關(guān)鍵位點施加外力進行拉伸，體系的總勢能Vtotal由原來的彈性勢能V變?yōu)閺椥詣菽芘c外力做功之和［33，37］，即：

其中，F(xiàn)表示施加在蛋白質(zhì)上的外力，為3N維列向量。對i、j殘基施加外力進行拉伸，則 外 力F矩 陣 的 構(gòu) 建 為 ：F=(000…Fix Fiy Fiz…Fjx Fjy Fjz…000)T1×3N。將體系的總勢能在平衡位置求梯度，可得到體系達(dá)到平衡后，各個殘基偏離平衡位置的位移為：

上式中，H-1矩陣的矩陣元為[H-1]ij=表示H-1矩陣的i行j列的矩陣 元，λkk為Hessian 矩 陣 的 第k個 本 征 值，U為Hessian 矩陣的本征向量，[Uk]i表示第k個本征向量中的第i個元素。殘基位移所導(dǎo)致的體系內(nèi)彈簧形變可以表示為［38-39］：

ΔrM×1為殘基間彈簧的形變，M表示組成彈性網(wǎng)絡(luò)的彈簧總數(shù)，AT是殘基偏移平衡位置位移在彈簧形變矢量上的投影。殘基間相應(yīng)的內(nèi)力分布可以表示為：

K是彈簧彈性系數(shù)所構(gòu)成的M維的對角矩陣，f是體系中所有彈簧承受的力，為M×1維的向量。

2 結(jié)果與討論

2.1 hPGK結(jié)構(gòu)域中的關(guān)鍵殘基及變構(gòu)路徑

hPGK 在生物體內(nèi)參與糖酵解過程，是催化1，3 二磷酸甘油酸轉(zhuǎn)化為3-磷酸甘油酸的關(guān)鍵酶，可產(chǎn)生1分子的ATP，在糖酵解的過程中發(fā)揮著不可或缺的作用［25，40］。hPGK 的三維結(jié)構(gòu)（PDB 代碼2XE7）由兩個大小相近的結(jié)構(gòu)域組成（圖1a）。研究證實，hPGK兩個結(jié)構(gòu)域的開合運動在催化功能的發(fā)揮中起著重要作用，而Arg38、Glu192、Asn336、Gly394 殘基是兩個結(jié)構(gòu)域打開與閉合過程中的關(guān)鍵殘基［25］。其中，殘基Arg38 和Asn336直接參與了底物的結(jié)合。殘基Glu192 和Gly394 位于體系兩個結(jié)構(gòu)域之間的鉸鏈區(qū)，直接調(diào)控著兩個結(jié)構(gòu)域的相對運動。大量的實驗結(jié)果表明，底物的結(jié)合會引起鉸鏈區(qū)的構(gòu)象調(diào)整，從而導(dǎo)致hPGK的兩個結(jié)構(gòu)域發(fā)生大幅度的構(gòu)象變化，從開放構(gòu)象轉(zhuǎn)變到閉合構(gòu)象，進而催化底物的化學(xué)反應(yīng)，完成其生物學(xué)功能［41-44］。底物的結(jié)合如何調(diào)控體系鉸鏈區(qū)的構(gòu)象運動，變構(gòu)信號如何從底物結(jié)合位點有效傳遞到體系的鉸鏈區(qū)是值得研究的重要問題，變構(gòu)路徑的識別研究有助于更好地揭示hPGK發(fā)揮功能的分子機制。

在本工作中，首先基于hPGK的三維結(jié)構(gòu)構(gòu)建彈性網(wǎng)絡(luò)，然后對殘基Arg38 和Asn336 施加1 pN的外力進行拉伸，拉伸方向為兩殘基位點連線的方向。進而，通過公式（9）和（10）計算蛋白質(zhì)在外力作用下的形變情況和相應(yīng)的內(nèi)力分布情況。計算得到體系內(nèi)力的分布情況如圖1b 所示。圖中藍(lán)色到紅色表示力的分布從小到大，受力越大表明相應(yīng)殘基對于力信號的傳導(dǎo)越重要，相應(yīng)位點為力信號變構(gòu)傳導(dǎo)的關(guān)鍵位點。

Fig.1 Three-dimensional structure of hPGK(a)and internal force distribution of hPGK in response to the exertion of external forces(b)

將圖1b 中不同殘基之間彈簧的受力情況映射到蛋白質(zhì)的三級結(jié)構(gòu)上，并根據(jù)不同殘基間力的分布情況來識別力信號的傳播路徑（圖2）。將外力所施加的兩個殘基位點Arg38 和Asn336 在圖中用紫色球標(biāo)出，為了更清晰的顯示力信號的傳播路徑，將受力大于0.05 pN 的彈簧顯示在圖2 中。圖中殘基間連線的粗細(xì)表示其承受的內(nèi)力大小，連線越粗，表示所承受的力越強，對力信號傳導(dǎo)也越重要；反之，連線越細(xì)表示受力越弱。由殘基間內(nèi)力分布情況，將底物結(jié)合位點與鉸鏈區(qū)之間受力最大的彈簧相連（圖2 紅色連線），得到力信號的傳導(dǎo)路徑。根據(jù)該方法，從Arg38出發(fā)，可以得到兩條到達(dá)鉸鏈區(qū)的變構(gòu)路徑：Arg38→Thr393→Gly394（鉸 鏈 區(qū)） 和Arg38→Thr393→Ala41→Leu188→Glu192（鉸鏈區(qū)）。這兩條變構(gòu)路徑用紅色箭頭顯示在了圖2中，其中所涉及到的關(guān)鍵殘基用黃色球標(biāo)出。采用類似方法，從Asn336 施力位點出發(fā)，也可以得到兩條傳播到鉸鏈區(qū)的力信號傳導(dǎo)路徑：Asn336→Gly371→Thr393→Gly394 （鉸 鏈 區(qū)） 和Asn336→Gly371→Thr393→Ser392→Glu192 （ 鉸鏈），傳導(dǎo)路徑在圖2 中用藍(lán)色箭頭標(biāo)出。變構(gòu)路徑中殘基間的內(nèi)力值見表1。

Fig.2 Allosteric pathways in hPGK identified by our method

Table 1 Internal forces between residues in the allosteric pathway of hPGK in response to the exertion of external forces

Szabo 等［45］通過殘基突變、酶反應(yīng)動力學(xué)實驗以及結(jié)構(gòu)分析，研究了底物結(jié)合對hPGK功能性開合構(gòu)象運動的變構(gòu)調(diào)控過程，并識別了變構(gòu)信號從底物結(jié)合位點到構(gòu)象運動鉸鏈區(qū)傳播過程中的關(guān)鍵殘基，發(fā)現(xiàn)Arg38、Lsy219、Asn336、Glu343等殘基在結(jié)構(gòu)域開合過程中發(fā)揮著重要作用。另外，Palmai 等［25］利用MD 模擬與FDA 相結(jié)合的方法，研究了底物結(jié)合對鉸鏈區(qū)構(gòu)象轉(zhuǎn)變的變構(gòu)調(diào)控作用，得到了兩條長程變構(gòu)路徑，分別是從結(jié)合位點殘基Arg38 傳遞到Thr393 再到鉸鏈區(qū)殘基Glu192、Gly394，以及從結(jié)合位點殘基Asn336 到Gly371、Ser398、Ser392 再傳遞到鉸鏈位區(qū)Glu192 和Gly394。本文的計算結(jié)果與上述實驗和MD模擬結(jié)果能夠很好地吻合。

2.2 PTP PDZ2結(jié)構(gòu)域的關(guān)鍵殘基及變構(gòu)路徑

PTP是生物體內(nèi)的關(guān)鍵酶［46-48］。它與蛋白酪氨酸激酶協(xié)同作用，共同維持著酪氨酸磷酸化的平衡，進而參與體內(nèi)的細(xì)胞生長、新陳代謝等活動。PDZ 結(jié)構(gòu)域存在于許多蛋白質(zhì)內(nèi)，被認(rèn)為具有蛋白質(zhì)定位、信號傳輸?shù)裙δ埽?6，49-50］。PTP 包含有5個PDZ 結(jié)構(gòu)域，其中，第二個結(jié)構(gòu)域由6 條β 鏈（β1～β6）和2條α螺旋（α1、α2）組成，含有96個氨基酸，屬于典型的PDZ 結(jié)構(gòu)［51-52］。圖3a 中不同的顏色表示不同的局部二級結(jié)構(gòu)。β2與α2螺旋之間存在一個配體結(jié)合口袋，實驗研究表明，配體的結(jié)合能夠長程變構(gòu)調(diào)控遠(yuǎn)端的α1和β1區(qū)域與其他蛋白質(zhì)的相互作用，變構(gòu)信號如何從配體結(jié)合口袋傳遞到遠(yuǎn)端區(qū)域是有待深入研究的重要問題，有助于更好揭示PDZ2功能的分子機制。

Fig.3 The three-dimensional structure of PTP PDZ2(a)and internal force distribution of PTP PDZ2 in response to the exertion of external forces(b)

本研究基于PTP PDZ2 的三級結(jié)構(gòu)（PDB 代碼3PDZ）構(gòu)建相應(yīng)的彈性網(wǎng)絡(luò)，然后對配體結(jié)合口袋區(qū)域的變構(gòu)位點殘基Gly19和His71施加1 pN外力進行拉伸，拉伸方向為兩殘基的連線方向。運用公式（10）計算蛋白質(zhì)體系在外力作用下的內(nèi)力分布情況（圖3b）。圖3b中的紅色區(qū)域?qū)?yīng)于承受內(nèi)力較大的殘基，是力信號傳導(dǎo)過程中的關(guān)鍵殘基。

為了更清晰地分析力信號的傳導(dǎo)路徑，將計算得到的力分布顯示到蛋白質(zhì)的三級結(jié)構(gòu)上（圖4）。圖中僅畫出了受力大于0.01 pN 的彈簧，顯示為綠色連線。連線的粗細(xì)表示相應(yīng)彈簧承受的內(nèi)力情況，連線越粗表示彈簧所承受的內(nèi)力越強。由殘基間內(nèi)力分布情況，將施力位點與遠(yuǎn)端β1區(qū)域間受力最大的彈簧相連，得到兩條特異性變構(gòu)路徑（圖4紅色連線）。一條變構(gòu)路徑從β2的施力位點開始，傳遞到α1螺旋的N 端，最終到達(dá)β1-β2之間的loop區(qū)，即Gly19、Ser21（β2）→Ⅰle41、Gln43、Gly45（α1）→Ser17（β1-β2之間的loop區(qū)），如圖4右側(cè)上圖中紅色箭頭所示，變構(gòu)路徑中的關(guān)鍵殘基用黃色球顯示在圖中。另一條力信號的傳導(dǎo)路徑是α2螺旋的施力位點開始，穿過β2、β3、β4、β6，最終到達(dá)β1，即His71 （α2） →Ⅴal22、Ⅰle20 （β2） →Ⅴal37（β3）→Gly55（β4）→Gln93、Lys91（β6）→Gly4（β1），如圖4右側(cè)下圖中紅色箭頭所示，變構(gòu)路徑中的關(guān)鍵殘基用橙色球顯示在圖中。變構(gòu)路徑中殘基間的內(nèi)力值見表2。

Fig.4 Allosteric pathways in PTP PDZ2 identified by our method

Table 2 Internal forces between residues in the allosteric pathway of PTP PDZ2 in response to the exertion of external forces

Kong等［48］基于MD模擬，采用殘基間相互作用的關(guān)聯(lián)性分析，研究了PTP PDZ2體系的長程變構(gòu)調(diào)控機制，識別出兩條變構(gòu)路徑，一條路徑始于β2上的變構(gòu)位點長程傳遞到α1螺旋，另一條路徑始于α2上的配體結(jié)合位點，垂直穿過β2、β3、β4、β6最終到達(dá)β1鏈的N 端。本研究利用基于ENM 的力分布分析方法得到的變構(gòu)路徑與Kong 等利用MD模擬得到的結(jié)果能夠很好地吻合。相對于MD 模擬，基于ENM 的力分布分析方法計算更為簡單、快速。

3 結(jié) 論

本文利用基于ENM 的力分布計算方法，通過分析蛋白質(zhì)對外力的響應(yīng)，研究了hPGK、PTP PDZ2 兩個蛋白質(zhì)體系力信號的變構(gòu)調(diào)控機制，成功識別出體系中特異的變構(gòu)路徑和關(guān)鍵殘基，與實驗和其他計算模擬結(jié)果能夠很好地吻合。該方法將蛋白質(zhì)的三級結(jié)構(gòu)簡化為彈性網(wǎng)絡(luò)，通過對變構(gòu)位點施加外力，計算蛋白質(zhì)體系在外力作用下的形變情況和體系內(nèi)部的力分布情況，來識別力信號傳播過程中起重要作用的關(guān)鍵殘基和相應(yīng)的變構(gòu)傳播路徑。相對于MD 模擬方法，基于ENM 的力分布分析方法具有明顯的優(yōu)勢，MD模擬需要耗費大量的模擬時間，并且需要利用有效的統(tǒng)計分析方法，去除噪音的干擾，來識別變構(gòu)信號，而ENM 方法可以通過解析求解來獲得體系的力分布情況，不需要進行數(shù)值模擬和軌跡采樣，也不包含噪聲，計算更加簡潔快速，能夠用于超大蛋白質(zhì)體系的研究。因此，本文為蛋白質(zhì)的變構(gòu)效應(yīng)研究提供了一種簡單、有效的方法。但是，需要指出的是，該方法也存在一定的不足，由于ENM 是一種對蛋白質(zhì)結(jié)構(gòu)進行高度簡化的粗?；Ｐ停總€氨基酸簡化為一個點，很難考慮不同配體結(jié)合、殘基的不同突變、翻譯后修飾等相互作用的精細(xì)差別，要更為精確的研究這些不同類型相互作用的影響，還需要結(jié)合更精細(xì)的全原子模型。