结合分子模拟和网络方法分析β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保存了轨迹用于分析。
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-645
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-2295
摘 要 顺应蛋白质结构功能研究的发展对构象信息的需求，结合靶向分子动力学模拟和蛋白质网络方法建立了对β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保存了轨迹用于分析。
2.3 复杂网络方法
2.3.1 氨基酸网络的构建 基于靶向分子模拟所得的5000个构象，本研究采用现有研究中相同的方法（如下公式2和2）构建氨基酸网络[24]。在氨基酸网络MCα中，以氨基酸Cα原子代表整个氨基酸并作为网络中的节点。边的构建方法：当任意一对氨基酸
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-645
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-2295
