結合分子模擬和網絡方法分析β2腎上腺素受體激活過程的構象變化

肖秀嬋+華勇攀+高楠+　江源遠+蒲雪梅+李夢龍

摘 要 順應蛋白質結構功能研究的發展對構象信息的需求，結合靶向分子動力學模擬和蛋白質網絡方法建立了對β2腎上腺素受體激活信號傳導相關的構象變化的解析方法，通過分子動力學模擬獲取激活過程的構象集合，并進一步采用網絡分析方法去識別構象變化過程中的關鍵殘基及關鍵路徑。結果表明，β2腎上腺素受體構象變化是各個區域之間信號傳遞的綜合，跨膜螺旋的連接區域是信號傳遞的樞紐，螺旋兩端包括胞內外Loop區域是信號傳導的核心區域。關鍵路徑分析發現，信號傳遞通路不止一條，都是從配體結合口袋殘基Ser204附近開始，分別傳遞到NPxxY區域和ICL2區域，其中螺旋Ⅲ，Ⅴ，Ⅵ是信號傳遞的必經區域，這些研究成果為β2AR激活信號傳導過程的闡明提供了重要的分子信息。

關鍵詞 β2腎上腺素受體； 激活過程； 構象變化； 靶向分子動力學模擬； 蛋白質網絡； 關鍵殘基

1 引 言

蛋白質具有柔性，其結構在功能發揮過程中會發生相應變化，而蛋白質的晶體結構只是一種靜態的平均結構，并不能完全代表其發揮功效的構象，因此對蛋白質構象變化的分析已成為蛋白質結構功能研究的一個重要發展趨勢[1]。Lanucara等[2]利用離子遷移質譜解析蛋白質的構象變化動力學，為人們理解其功能提供了蛋白質柔性和折疊機制的構象信息；Burmann等[3]將高分辨率NMR用于大腸桿菌探針與配體作用的構象變化動力學研究。此外，Kahsai等[4]提出了結合同位素標記和質譜法研究蛋白質構象變化的分析方法，并用該方法成功揭示了一種重要靶標蛋白受藥物配體誘導所引起的激活功能的構象變化。上述工作清楚地展示了儀器分析方法用于蛋白質構象分析的發展新趨勢，以及構象解析對一些新分析方法的需求。

G蛋白偶聯受體（GPCRs）是最大的膜蛋白家族和藥物靶標之一，由7個跨膜螺旋構成，在信號傳導過程中具有重要作用，是藥物研發中主要靶標蛋白，目前市場上30%的藥物是以GPCRs為靶標[5]。GPCRs主要是通過從非活性態到活性態的激活來發揮其信號功能，在其激活過程中，伴隨著明顯的構象變化，因此，對激活過程中構象變化的了解是闡明GPCRs激活機制這一信號功能的關鍵。然而，對于膜蛋白的解析難度導致GPCRs高分辨率晶體結構實驗缺乏，更加限制了其激活過程中構象變化的分析。因此，有必要引入一些先進的其它學科技術，協助實驗解析其構象變化的結構信息。

分子動力學方法可以模擬出物質在原子尺度上的結構和動力學行為的時間變化過程，因此已成為輔助實驗研究蛋白質動態結構（構象）變化的一種有效方法[6，7]。實驗研究已揭示GPCRs激活過程時間在ms～s級[8]，現有的計算機能力還難于采用平衡的分子動力學方法獲取GPCR的激活過程。靶向分子動力學模擬（TMD）是一種非平衡的分子動力學模擬，可以借助外力加速蛋白質從初始結構到目標結構的構象變化，在有限的時間內能獲取蛋白質功能過程中構象變化空間，已成功應用于一些涉及大的構象變化的蛋白質結構性能的研究中[9，10]。此外，分子動力學模擬軌跡中包含了眾多構象，如何對這些大量構象進行有效的分析，以獲取其重要的關鍵區域及通路是解析與功能相關的構象變化的關鍵。

復雜網絡方法可以利用網絡的原理通過對節點網絡特征值和邊的分析給出體系中信號變化的重要信息，蛋白質網絡就是采用復雜網絡方法對蛋白質拓撲結構進行分析，其中以氨基酸為頂點，以氨基酸之間的相互作用為邊構建的拓撲結構，可以將復雜的蛋白質結構以一種全新的網絡形式描述，而且也可將分子動力學模擬軌跡所得的每個構象視為節點，而構象之間的轉換關系就是網絡的邊，由此構建的蛋白質網絡可以有效并快速地獲取蛋白質動態變化中的氨基酸特征以及變化通路，這類蛋白質網絡分析方法已成功應用于蛋白質的穩定性[11]，蛋白質折疊機制[12，13]以及蛋白質相互作用[14]等方面的研究。

基于上面的研究背景，選擇β2腎上腺素能受體（β2AR）作為研究對象，該受體是2007年解析出的人類第一個GPCR受體，其對肌肉平滑肌、支氣管哮喘和喘息性支氣管炎都有重要的調節作用[15，16]。首先采用靶向分子動力學模擬方法來獲取β2AR從非活性狀態到活性狀態的激活過程的構象變化過程，再結合蛋白質網絡的方法來識別和揭示β2AR的激活過程中的關鍵殘基以及關鍵通路，為闡明GPCR基于激活的信號傳導機制提供有價值的分子信息。

2 實驗部分

2.1 初始結構的構建

非活性和活性狀態的初始坐標分別來自X-ray解析所得的β2AR晶體結構（PDB ID分別為2RH1和3SN6[17，18]。為了保證模擬環境真實性將去掉了所有非受體分子得到的2個apo受體[19]放入已構建好的磷脂雙層膜POPC中，添加TIP3P水分子，分別構建了非活性態和活性態兩個體系。

2.2 靶向分子動力學模擬

β2AR是由AMBER03力場描述，磷脂分子由GAFF（General amber force filed）力場描述[20]。TMD模擬采用的是NPT系綜，在Amber12程序[21]的Sander模塊中實現的，該方法用了一個標準的分子動力學勢能和一個與時間相關的限制勢能UTMD：

UTMD=12NK（RMSD-p（t）2）（1）

其中，N是RMSD值計算中所包含的原子數量，K是力常數，RMSD等于模擬結構初始和最終狀態之間的RMSD值，ρ（t）是模擬體系在時間t時相對于目標構象的RMSD值。在本研究中，TMD模擬的初始坐標分別是上面所得非活性態和活性態的坐標。此外，NPT系綜的模擬在310 K下，壓力采用了1個大氣壓的周期邊界條件，SHAKE算法被用于限制所有與H相連的鍵長，其中限制的閾值為1.0×10

5[22]。非鍵相互作用的截斷值采用10，Particle-Mesh-Ewald （PME） 方法[23]被用于處理長程的靜電相互作用。采用1 kcal mol1 2限制力和1 ns的模擬時間。MD模擬中步長為2 fs，且每隔1 ps保存了軌跡用于分析。endprint

2.3 復雜網絡方法

2.3.1 氨基酸網絡的構建 基于靶向分子模擬所得的5000個構象，本研究采用現有研究中相同的方法（如下公式2和2）構建氨基酸網絡[24]。在氨基酸網絡MCα中，以氨基酸Cα原子代表整個氨基酸并作為網絡中的節點。邊的構建方法：當任意一對氨基酸的Cα原子之間的歐幾里得距離dCαij≤7時[24]，兩個節點之間添加一條邊。

dCαij=3（2）

MCαij={1， dCαij≤7， i≠j

0 其它（3）

2.3.2 網絡特征的挑選

主要分析了氨基酸在激活過程中度、接近中心性、介數和K-核的期望值。在各種不同的網絡拓撲特征中， 這4個特征常用來衡量網絡中節點的重要性和權威性。

度表示網絡中節點連接的邊數，用k表示。度越大的節點說明該節點的鄰居比較多，在傳出信號時能夠快速地將信號轉發給多個接收者，并且能夠從多個方向接收信號。

接近中心性表示特定節點i到網絡中每個節點的最短平均距離的倒數，簡稱接近數，用符號CCi表示，公式如下：

Li=1N（4）

CCi=LLi=N∑Nj=1（Lij）（5）

其中，Lij表示節點i和節點j之間的最短路徑長度。Li表示從節點i到網絡中所有節點最短距離的平均值。Li值的大小能反映節點i在網絡中的相對重要性。Li越小，說明節點i距離網絡中所有節點的平均距離更近，信號由節點i向網絡中其它節點傳播更為容易，速度更快，信號衰減更小。

介數g（i）是網絡中經過節點i的最短路徑個數所占比例，介數值越大表明該節點是多條最短路徑的必經之處，即該節點的樞紐通道能力很強。

K-核衡量了節點位于網絡中的層次，能夠從一定程度上衡量節點在氨基酸網絡中信號傳導的重要性。

3 結果與討論

3.1 靶向分子動力學模擬所揭示與激活信號傳導過程相關的構象集合的分子信息

β2AR的實驗研究表明從非活性態轉換到活性態是伴隨復雜的構象變化的過程[15，25]。利用TMD方法模擬了β2AR從非活性態到活性態的變化過程（具體過程參照2.1和2.2節）。為了觀察受體整體結構的變化情況分析了整個模擬過程中以非活性態晶體結構為參照的結構的RMSD值波動，結果如圖1所示，受體在激活過程中其構象發生了很大的變化，在400 ps內達到了活性態的結構；在起始的100 ps內，是接近非活性態的結構，而在隨后的時間內逐漸向活性態變化，最后達到完全活性態。

圖1 β2AR靶向分子動力學模擬過程中RMSD值隨時間的變化（）

Fig.1 Changes in RMSD values （in ） of β2 adrenergic receptor （β2AR） as a function of time in the 1 ns target molecular dynamics （TMD） simulation （）

3.2 基于復雜網絡分析的激活過程中關鍵區域和關鍵殘基的識別

采用氨基酸網絡方法將靶向分子動力學模擬β2AR得到的5000個構象構建成相應的5000個氨基酸網絡，分析了每個氨基酸在氨基酸網絡變化過程中網絡拓撲特征的變化，發現氨基酸的變化情況，進而挖掘β2AR激活過程的信號傳遞中構象變化的關鍵區域及關鍵殘基。

計算TMD構象所構建的5000個氨基酸網絡中每一個氨基酸的度期望值、接近中心性期望值，介數期望值和K-核期望值，將結果映射到β2AR氨基酸網絡中，圖2是網絡特征值的氨基酸網絡熱點映射圖。

圖2 β2AR氨基酸網絡。度期望值（a），接近中心性期望值（b），介數期望值（c），K-核期望值（d）的分布，紅色的點代表值較小，綠色的點代表值較大

Fig.2 Amino acid network of β2 adrenergic receptor. the degree expected value （a）， Closeness expected value （b）， betweenness expected value （c）， K-core expected value（d）. The distribution of red dots represents small values， green dots represent big values

從圖2a可見，在整個網絡分析中度期望值較大的節點遍及所有的區域，包括跨膜螺旋區域以及胞內外Loop區，而度值越大的氨基酸的網絡鄰居越多，更容易將信號傳遞給更多的節點，圖中各個區域都有度值較大的殘基，說明受體各個區域之間信號傳遞有很高的相互關聯性。

圖2b是接近中心性的期望值，此值衡量了網絡中每一個節點與其余所有節點距離平均值的倒數，該網絡參數越大傳遞信號的能力越強，體現出殘基在構象變化中的重要性。圖2b可見，接近中心性較大的殘基大多位于螺旋的中間部分，即連接區域，揭示了該區域可作為激活信號從配體結合口袋傳遞到G蛋白結合區域的樞紐區域。

介數描述的是一個網絡中節點可能承載的信息流量，如果節點的介數越大，流經它的信息就越多。從圖2c可見，螺旋Ⅱ，Ⅲ，V，Ⅵ及Ⅶ的介數值較大，配體結合口袋處殘基（Ala200，Ser203，Ser204，Ser207）以及NPxxY區域較為明顯，說明這些區域在信號傳遞中是必經的氨基酸節點。

圖2d給出是K-核期望值，此值衡量了氨基酸在網絡中核心層次的位置，越接近核心位置的節點，相對信號越容易擴散到網絡中的每一個節點，螺旋兩端的K-核最大，尤其是胞內外loop區域（ECL2，ECL3，ICL2，靠近ICL3的殘基），說明在β2AR激活過程中胞內外區域是核心，為受體激活中信號從胞外loop區傳到胞內loop區的功能過程提供了進一步的分子證據。

綜上所述， 在受體激活過程中度值，接近中心性值以及介數值最大的區域在配體結合口袋Ser204附近。此外，在胞內區域網絡特征值最大的氨基酸節點并不明顯，而在受體激活過程中胞內區域變化較大的螺旋Ⅲ和Ⅵ的K-核值特征明顯，體現了激活過程中螺旋Ⅲ和螺旋Ⅵ構象變化的核心作用。這與一些實驗研究[26]所指出的螺旋Ⅲ和螺旋Ⅵ是G蛋白偶聯受體激活的一個重要特征是一致的，為實驗現象提供了分子水平上的構象證據，同時也表明了分子模擬的可靠性。

3.3 構象變化的關鍵路徑分析

上面的分析已揭示了β2AR激活過程中的一些重要區域，本部分進一步采用WebPSN[27]方法對整個激活過程的關鍵路徑進行了分析。關鍵路徑分析算法的核心思想：在網絡中，存在一個起始狀態類、一個終點狀態類和若干中間狀態類（按照靶向動力學變化方向排列）。而每一條關鍵路徑必須從起點狀態到達終點狀態類，途中經過每一個中間狀態的最短路徑。圖3顯示了信號從胞外傳遞到胞內的激活通路的3種情況：

圖3 β2AR激活的關鍵通路及關鍵殘基。黃色小球代表氨基酸，綠色線代表路徑

Fig.3 The key residues and the critical pathways of β2AR during activation process. Yellow balls represent amino acids， green lines represent pathway

圖3a中可以看出信號從配體結合口袋處殘基Ser204開始，經過螺旋Ⅴ→Ⅵ→Ⅲ→Ⅱ→Ⅳ，最后到達DRY區域以及胞內ICL2區域；圖3b信號從配體結合口袋處殘基Ser203開始，經過螺旋Ⅴ→Ⅲ→Ⅱ→Ⅳ，最后到達DRY區域以及胞內ICL2區域；而圖3c信號從配體結合口袋處殘基Ser204開始，經過螺旋Ⅴ→Ⅵ→Ⅲ→Ⅶ，最后到達NPxxY區域。

由此可知，信號從胞外傳遞到胞內的過程中配體結合口袋為必經的區域，且Ser204附近的殘基為關鍵氨基酸，與3.2節結果相同。從圖3還可見，受體激活通路并不止一條，圖3a和圖3b所展示的通路都是由配體結合口袋經過連接區域不同的氨基酸到達高保守的DRY區域[28]，最終為ICL2區域上的Tyr141，而由活性態的晶體結構可知ICL2區域正是G蛋白結合的關鍵區域[28]，驗證了結果的可靠性；圖3c所展示的是由配體結合口袋經過連接區域到達NPxxY區域的，而NPxxY區域被實驗研究指出是衡量G蛋白結合區域活性的關鍵區域[29]，與實驗結果相符。這些研究結果表明， 受體激活過程中會激活胞內不同區域，這應該是導致受體可以與不同的G蛋白以及β-抑制蛋白結合的原因。此外，圖3還顯示，通路中間殘基大多位于螺旋Ⅲ，Ⅴ和Ⅵ上，與3.2節介數分布所給出的重要殘基相吻合。

4 結 論

實驗采用靶向分子動力學模擬的計算手段獲取了重要的藥物靶標蛋白β2AR在其信號傳導激活過程中的構象變化的分子信息，并進一步采用復雜網絡方法對信號傳導過程中構象變化的關鍵區域和關鍵路徑進行了識別分析，研究結果發現β2AR激活是各個區域之間信號傳遞的綜合，跨膜螺旋的連接區域是信號傳遞的樞紐，螺旋兩端包括胞內外Loop區域是信號傳導的核心區域。關鍵路徑分析發現信號傳遞通路不止一條，都是從配體結合口袋殘基Ser204附近開始，分別傳遞到NPxxY區域和ICL2區域，其中螺旋Ⅲ，Ⅴ，Ⅵ是信號傳遞的必經區域。這些研究結果為闡明G蛋白偶聯受體的功能機制提供了重要的分子信息，并進一步展示了計算機模擬和網絡分析技術在蛋白質結構功能分析中的重要作用和發展潛力。

References

1 Vardy E， Roth B L. Cell， 2013， 152（3）： 385-386

2 Lanucara F， Holman S W， Gray C J， Eyers C E. Nat. Chem.， 2014， 6（4）： 281-294

3 Burmann B M， Wang C， Hiller S. Nat. Struct. Mol. Biol.， 2013， 20（11）： 1265-1272

4 Kahsai A W， Rajagopal S， Sun J， Xiao K. Nat. Protoc.， 2014， 9（6）： 1301-1319

5 Gether U， Asmar F， Meinild A K， Rasmussen S G. Pharm. Toxicol.， 2002， 91（6）： 304-312

6 Cui Y L， Zhang J L， Zheng Q C， Niu R J， Xu Y， Zhang H X， Sun C C. Chem. Eur. J.， 2013， 19（2）： 549-557

7 Sun H Y， Ji F Q. Biochem. Biophys.Res. Co.， 2012， 423（2）： 319-324

8 Swift R V， McCammon J A. Biochemistry， 2008， 47（13）： 4102-4111

9 Zhang J， Li C， Chen K， Zhu W， Shen X， Jiang H. Proc. Natl. Acad. Sci. USA， 2006， 103（36）： 13368-13373

10 Weng J， Fan K， Wang W. PloS one， 2012， 7（1）： e30465

11 Vendruscolo M， Paci E， Dobson CM， Karplus M. Nature， 2001， 409（6820）： 641-645endprint

12 Dokholyan N V， Li L， Ding F， Shakhnovich E I. Proc. Natl. Acad. Sci. USA， 2002， 99（13）： 8637-8641

13 Vendruscolo M， Dokholyan N， Paci E， Karplus M. Phys. Rev. E.， 2002， 65（6）： 1-4

14 Thibert B， Bredesen DE， del Rio G. BMC Bioinformatics， 2005， 6（213）： 1-15

15 Feng Z， Hou T， Li Y. J. Chem. Inf Model.， 2012， 52（4）： 1005-1014

16 Litonjua A A. Curr. Opin. Pulm. Med.， 2006， 12（1）： 12-17

17 Rasmussen S G， Choi H J， Rosenbaum D M， Kobilka T S， Thian F S， Edwards P C， Burghammer M， Ratnala V R， Sanishvili R， Fischetti R F. Nature， 2007， 450（7168）： 383-387

18 Rasmussen S G， DeVree B T， Zou Y， Kruse A C， Chung K Y， Kobilka T S， Thian F S， Chae P S， Pardon E， Calinski D. Nature， 2011， 477 （7366）： 549-555

19 Rubenstein R C， Wong S， Ross E. J. Biol. Chem.， 1987， 262（34）： 16655-16662

20 Pearlman D A， Case D A， Caldwell J W， Ross W S， Cheatham III T E， DeBolt S， Ferguson D， Seibel G， Kollman P. Comput.Phys. Commun.， 1995， 91（1）： 1-41

21 Dickson C J， Rosso L， Betz R M， Walker R C， Gould I R. Soft Matter.， 2012， 8（37）： 9617-9627

22 Ryckaert J P， Ciccotti G， Berendsen H J. J. Comput. Phys.， 1977， 23（3）： 327-341

23 Essmann U， Perera L， Berkowitz M L， Darden T， Lee H， Pedersen L G. J. Chem. Phys.， 1995， 103（19）： 8577-8593

24 Yan W， Zhou J， Sun M， Chen J， Hu G， Shen B. Amino Acids， 2014， 46（6）： 1419-1439

25 Chung K Y， Rasmussen S G， Liu T， Li S， DeVree B T， Chae P S， Calinski D， Kobilka B K， Woods V L， Jr. Sunahara R K. Nature， 2011， 477（7366）： 611-615

26 Dror R O， Arlow D H， Borhani D W， Jensen M O， Piana S， Shaw D E. Proc. Natl. Acad. Sci. USA， 2009， 106（12）： 4689-4694

27 Seeber M， Felline A， Raimondi F， Mariani S， Fanelli F. Bioinformatics， 2015， 31（5）： 779-781

28 Szalai B， Barkai L， Turu G， Szidonya L， Varnai P， Hunyady L. Biochem. Pharmacol.， 2012， 84（4）： 477-485

29 Fritze O， Filipek S， Kuksa V， Palczewski K， Hofmann K P， Ernst O P. Proc. Natl. Acad. Sci. USA.， 2003， 100（5）： 2290-2295endprint