基于各向异性网络模型分析A_2A腺苷受体变构过程中的关键残基

Analysis of key residues in A_2A adenosine receptor metamorphosis based on the anisotropic network model

摘要

A_2A腺苷受体（A_2A Adenosine receptor，A_2AAR）是A类G蛋白偶联受体的重要成员，在内源性激动剂腺苷介导的跨膜信号传导中发挥关键作用，参与许多重要的信号调节过程。基于各向异性网络模型，分析了A_2AAR非激活态、中间态以及激活态的快、慢运动模式残基的均方涨落，识别了A_2AAR变构过程中与结构稳定性及功能相关的关键残基。利用残基对间运动的相关性与相互作用，分析了A_2AAR变构过程中二级结构间的相关性和维持结构稳定的重要残基，研究结果为相关疾病与药物开发提供理论帮助。

Abstract

A_2A adenosine receptor (A_2AAR), an important member of class A G-protein-coupled receptors, plays a critical role in transmembrane signaling mediated by endogenous agonist adenosine and is involved in many important signal-regulatory processes. Using the anisotropic network model, we analyzed the mean-square fluctuations of residues in the inactive, intermediate, and active states of A_2AAR in both fast and slow-motion modes, and found the key residues that contribute to the structural stability and function of A_2AAR during its conformational changes. By examining the correlation and interaction of inter-residue pair motions, we investigated the correlation between the secondary structures and the important residues to maintain the structural stability in the A_2AAR allosteric process. The findings of this study provide theoretical guidance for the study of related diseases and drug development.

Graphical abstract

关键词

A_2A腺苷受体 / 各向异性网络模型 / 均方涨落 / 关键残基

Key words

A_2A adenosine receptor / anisotropic network model / mean-square fluctuations / key residues

G蛋白偶联受体(G-protein coupled receptors，GPCRs)是负责细胞膜信号传导的最大类膜蛋白家族。GPCRs介导了多数细胞对激素、神经递质、光与气味等外界环境的生理反应，具有作为重要药理学治疗相关疾病靶点的巨大潜力^[1]。腺苷受体(Adenosine receptors，ARs)是A类G蛋白偶联受体的重要亚家族成员，嘌呤腺苷是其内源性配体，通常ARs分为四种亚型受体，即A₁、A_2A、A_2B和A₃受体^[2]。ARs在身体各种组织中的广泛分布和表达，以及其在控制身体不同功能中的作用，使其成为治疗各种病理状况的潜在药物靶点，如心脏病、癌症、帕金森等其他疾病^[3]。A_2A腺苷受体(A_2A Adenosine receptor，A_2AAR)已被广泛地研究，并发现了多种与其结合的配体，包括激动剂、拮抗剂、反向激动剂和变构调节剂^[4]。

在实验方面，Lebon等^[5]通过比对A_2AAR与不同激动剂结合的晶体结构，发现配体的腺苷部分以相同的方式与受体结合，并指出与不同配体结合过程中影响构象变化的关键残基及相互作用方式。Jaakola等^[6]解析出A_2A与拮抗剂ZM241385结合的复合物晶体结构，基于结构特性和残基间的相互作用方式确定了配体结合过程中的关键残基。随着A_2AAR晶体结构被解析出，已成功地基于结构进行分子建模为A_2AAR设计新的药物结构^[7]。在理论方面，分子动力学(Molecular dynamics，MD)模拟在解决GPCRs里的配体结合和激活过程关键问题中取得了巨大贡献^[8-9]。MARTÍNEZ等^[10]应用MD方法对活性和非活性的A_2AAR进行模拟，研究激活或失活过程中受体的结构特征，并揭示了受体激活过程中的关键残基。Lovera等^[11]应用自适应采样的MD模拟与目标导向的评分函数建立马尔科夫状态模型(Markov state modelling，MSM)重构了A_2AAR的激活途径同时分析了激活过程中关键残基的特征变化。MD模拟虽然能实现蛋白质发生功能性运动时的构像变化但计算量巨大，为解决这个问题，Brooks等^[12]发展了粗粒化的正则模分析(Normal model analysis，NMA)，并成功用于蛋白质运动模式的研究。Tirion^[13]基于NMA提出了弹性网络模型(Elastic network model，ENM)，并表明其计算精度不输于传统NMA方法且速度有所提升。Haliloglu等^[14]对ENM做了进一步简化，提出了每一个氨基酸只考虑Cα原子的高斯网络模型(Gaussian network model，GNM)，该模型将氨基酸简化为一个点，大大降低了计算量，能够计算体系更大的蛋白质。

由于GNM认为蛋白质中所有残基的运动是各向同性的，没有考虑残基运动方向问题，根据这一问题，Atilgan等^[15]提出了各向异性网络模型(Anisotropic network model，ANM)，ANM不仅能够得到残基的运动幅度，还能得到各残基的运动方向，残基的运动是各向异性的，这一模型相比于GNM能够更好地模拟蛋白质功能性运动。Yang等^[16]使用ANM来预测蛋白质的温度因子与实验值比较。陈磊等^[17]使用ANM研究了\mathrm{\delta }阿片受体的动力学特征以及对变构具有重要作用的关键残基。Su等^[18]使用ANM研究了维生素B₁₂转运蛋白BtuCD-BtuF的构象运动，并预测了与功能运动相关的关键残基。

本研究基于ANM探索了A_2AAR激活过程中的动力学与功能性残基，根据模型快慢模式运动分析蛋白质功能运动相关的残基，同时根据原子的均方涨落振幅大小及运动相关性划分蛋白质柔性结构域和关键残基。

1 材料与方法

1.1 数据

从GPCRdb数据库下载了同源建模的未激活态(PDB代码：4EIY)、中间态(PDB代码：4UG2)和激活态(PDB代码：5G53)A_2AAR晶体结构数据。A_2AAR由412个残基组成，其序列见图1。7个跨膜螺旋对应的残基片段分别是TM1(Met1^1.27-Asn34^1.60)、TM2(Asn39^2.37-Gly69^2.66)、TM3(Ala73^3.21-Ile108^3.56)、TM4(Thr117^4.38-Gly142^4.63)、TM5(Pro173^5.36-Ser213^5.74)、TM6(Glu219^6.21-Cys259^6.61)、TM7(Pro266^7.30-Arg291^7.56)。1个不穿过细胞膜的短螺旋H8(Ile292^8.47-Phe314^8.69)位于膜-细胞质界面。残基编号采用的是Ballesteros-Weinstein通用编号方案^[19]，该方案提供了每个残基的相对位置，如Asp52^2.50的上标2代表天冬氨酸位于TM2，50表示最保守的残基，其他残基编号根据最保守残基的位置确定。在4EIY的结构数据中Lys209^5.70-Gly218^ICL3结构缺失，4UG2中Lys209^5.70-Arg220^6.22结构缺失，5G53中Gly147^ECL2-Gly158^ECL2和Glu212^5.73-Ser223^6.25结构缺失，因此在模型计算时不考虑这些结构的残基。

1.2 各向异性网络模型

ANM将蛋白质的三维空间结构简化为一个弹性网络。蛋白质中的每个残基简化为一个点，用其Cα原子代替，残基之间的相互作用简化为谐振势，用弹簧来代替。当残基之间的距离小于截断半径Rc时存在相互作用，体系的势能V等于所有弹簧的势能之和：

V=\frac{\mathrm{1}}{\mathrm{2}}\gamma {\sum }_{ij}(∆R_{ij}{)}^{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}\gamma \left\{∆{\vec{R}}^{\mathrm{T}}\right\}\mathrm{H}\left\{∆\vec{R}\right\}(1)

其中：\gamma代表体系中弹簧的劲度系数，\left\{∆\vec{R}\right\}是残基波动的3N维列向量，T表示向量转置，H是海森矩阵，其形式为

\mathrm{H}=\left[\begin{array}{cc}\begin{array}{cc}{\mathrm{H}}_{\mathrm{11}}& {\mathrm{H}}_{\mathrm{12}}\\ {\mathrm{H}}_{\mathrm{21}}& {\mathrm{H}}_{\mathrm{22}}\end{array}& \begin{array}{cc}\cdots & {\mathrm{H}}_{\mathrm{1}N}\\ \cdots & {\mathrm{H}}_{\mathrm{2}N}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ {\mathrm{H}}_{N\mathrm{1}}& {\mathrm{H}}_{N\mathrm{2}}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \cdots & {\mathrm{H}}_{NN}\end{array}\end{array}\right](2)

其中，H _ij 是三维子矩阵，其元素给定为

{\mathrm{H}}_{ij}=\left[\begin{array}{ccc}\frac{{\partial }^{\mathrm{2}}V}{\partial {\mathrm{x}}_i\partial {\mathrm{x}}_j}& \frac{{\partial }^{\mathrm{2}}V}{\partial {\mathrm{x}}_i\partial {\mathrm{y}}_j}& \frac{{\partial }^{\mathrm{2}}V}{\partial {\mathrm{x}}_i\partial {\mathrm{z}}_j}\\ \frac{{\partial }^{\mathrm{2}}V}{\partial {\mathrm{y}}_i\partial {\mathrm{x}}_j}& \frac{{\partial }^{\mathrm{2}}V}{\partial {\mathrm{y}}_i\partial {\mathrm{y}}_j}& \frac{{\partial }^{\mathrm{2}}V}{\partial {\mathrm{y}}_i\partial {\mathrm{z}}_j}\\ \frac{{\partial }^{\mathrm{2}}V}{\partial {\mathrm{z}}_i\partial {\mathrm{x}}_j}& \frac{{\partial }^{\mathrm{2}}V}{\partial {\mathrm{z}}_i\partial {\mathrm{y}}_j}& \frac{{\partial }^{\mathrm{2}}V}{\partial {\mathrm{z}}_i\partial {\mathrm{z}}_j}\end{array}\right](3)

海森矩阵H可以根据相似对角化进行分解

{\mathrm{H}}^{-\mathrm{1}}=\mathrm{U}{\mathrm{\Lambda }}^{-\mathrm{1}}{\mathrm{U}}^{\mathrm{T}}={\sum }_{i=\mathrm{2}}^N{\mathrm{\lambda }}_i^{-\mathrm{1}}{\mathrm{u}}_i{\mathrm{u}}_i^{\mathrm{T}}(4)

其中，U是正交矩阵，其列u _i 代表H的第i个本征向量，Λ是对角矩阵，其对角元素λ _i 是H的本征值，N为蛋白质体系内残基的个数。

残基i在第k个运动模式下的均方涨落可以写为

\begin{array}{l}{〈{\left(\Delta {\vec{R}}_i\right)}^{\mathrm{2}}〉}_k=\frac{k_{B}T}{\gamma }({\mathrm{\lambda }}_k^{-\mathrm{1}}[{\mathrm{u}}_k{]}_{\mathrm{3}i-\mathrm{2}}[{\mathrm{u}}_k{]}_{\mathrm{3}i-\mathrm{2}}\\ +{\mathrm{\lambda }}_k^{-\mathrm{1}}[{\mathrm{u}}_k{]}_{\mathrm{3}i-\mathrm{1}}[{\mathrm{u}}_k{]}_{\mathrm{3}i-\mathrm{1}}+{\mathrm{\lambda }}_k^{-\mathrm{1}}[{\mathrm{u}}_k{]}_{\mathrm{3}i}[{\mathrm{u}}_k{]}_{\mathrm{3}i})\end{array}(5)

上式[{\mathrm{u}}_k{]}_{\mathrm{3}i}是H矩阵的第k个本征向量中的第3i个元素，T为温度。

在ANM中，残基的均方涨落以及残基间涨落的交叉关联性为

〈{\left(\Delta {\vec{R}}_i\right)}^{\mathrm{2}}〉=\frac{k_{B}T}{\gamma }\left(H_{\mathrm{3}i-\mathrm{2},\mathrm{3}i-\mathrm{2}}^{-\mathrm{1}}+H_{\mathrm{3}i-\mathrm{1},\mathrm{3}i-\mathrm{1}}^{-\mathrm{1}}+H_{\mathrm{3}i,\mathrm{3}i}^{-\mathrm{1}}\right)(6)

〈\Delta {\vec{R}}_i\cdot \Delta {\vec{R}}_j〉=\frac{k_{B}T}{\gamma }\left(H_{\mathrm{3}i-\mathrm{2},\mathrm{3}j-\mathrm{2}}^{-\mathrm{1}}+H_{\mathrm{3}i-\mathrm{1},\mathrm{3}j-\mathrm{1}}^{-\mathrm{1}}+H_{\mathrm{3}i,\mathrm{3}j}^{-\mathrm{1}}\right)(7)

归一化的互相关系数为

C_{ij}=\frac{〈\mathrm{\Delta }{\vec{R}}_i\cdot \mathrm{\Delta }{\vec{R}}_j〉}{{\left[〈\mathrm{\Delta }{\vec{R}}_i^{\mathrm{2}}〉\times 〈\mathrm{\Delta }{\vec{R}}_j^{\mathrm{2}}〉\right]}^{\mathrm{1}/\mathrm{2}}}(8)

第i个残基的温度因子(B-factor)可以根据下式计算

B_i=\frac{\mathrm{8}{\mathrm{\pi }}^{\mathrm{2}}}{\mathrm{3}}〈{\left(\mathrm{\Delta }{\vec{R}}_i\right)}^{\mathrm{2}}〉(9)

残基i的变形能为

E_i=\frac{\mathrm{1}}{\mathrm{2}}{\sum }_{j=\mathrm{1}}^{N}k\left(R_{ij}^{\left(\mathrm{0}\right)}\right)\frac{{\left|\left(d_i-d_j\right)R_{ij}^{\left(\mathrm{0}\right)}\right|}^{\mathrm{2}}}{{\left|R_{ij}^{\left(\mathrm{0}\right)}\right|}^{\mathrm{2}}}(10)

k\left(r\right)=c\cdot \mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\left|r\right|}^{\mathrm{2}}}{r_{\mathrm{0}}^{\mathrm{2}}}\right)(11)

其中，d_i 是待分析模式下原子i的位移量，k(r)代表成对力常数，r_{\mathrm{0}}给定的距离截止值，R_{ij}^{\left(\mathrm{0}\right)}表示初始距离向量，残基i到j的距离。

ANM模型对蛋白质数据的分辨率要求较高，但可以计算体系很大的蛋白质，其计算精度不输MD，模拟且计算速度很快。

1.3 残基的功能性运动柔性量度

采样残基的均方涨落和变形能的组合来衡量残基的功能性运动，但由于残基的均方涨落与变形能的量纲不一致，故使用归一化的无量纲数进行线性加权来定义残基的功能性运动柔性量度Δ：

x^s=\frac{x-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{x_{\mathrm{m}\mathrm{a}\mathrm{x}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}(12)

因此第i个残基的柔性量Δ定义为

{∆}_i=\mathrm{0.5}{〈{\left(\mathrm{\Delta }{\vec{R}}_i\right)}^{\mathrm{2}}〉}^s+\mathrm{0.5}{E}_i^s(13)

2 结果与讨论

2.1 基于ANM预测A_2AAR各残基B-factor值

在对蛋白质结构的研究中，B-factor可被看作是结构中不同部位的相对振动，当B-factor较小时可以体现出该原子处于晶格结构的稳定有序部分，而B-factor较大时则表现出该原子处于柔性易变的晶格结构中^[20]。本文通过对A_2A腺苷受体建立粗粒化弹性网络，选取残基的Cα原子为节点建立ANM对A_2AAR的B-factor进行预测。根据Debye-Waller理论，用残基的均方涨落预测B-factor^[21]。模型的可调参数截断半径(Rc)越大，每个残基相连的弹簧数越多，受到的作用也随之增加，因此残基在力场中产生的波动与参数截断半径相关^[16]。不同的截断半径直接影响预测结果，首先对实验B-factor值进行Shapiro-Wilk正态性检验，表1中的检验结果p值远小于0.05，故B-factor不服从正态分布，使用Spearman相关系数作为衡量标准。截断半径的最优取值是由B-factor的ANM预测值与实验值之间的斯皮尔曼相关系数的最大化来确定，截断半径的寻优范围为7~50 Å，步长为1 Å。

图2显示了B-factor的ANM预测值与实验值之间的Spearman相关系数随截断半径的变化曲线。由图2可以看出，对于A_2AAR非激活态(4EIY)、中间态(4UG2)和激活态(5G53)其最优截断半径分别为Rc=30、26、28 Å，相应的Spearman相关系数分别为ρ=0.811、0.732、0.674。不难发现同一截断半径下ANM对B-factor预测值的准确性随着受体的激活状态依次降低，激活态的受体柔性波动最大，其结构最为复杂。最优截断半径下A_2AAR的三种状态温度因子的预测值与实验值的比较见图3。由图3可知ANM能较好地再现腺苷受体的柔性信息，TM1-TM7刚性较强，残基的波动较小，其预测值与实验值吻合度较高，而未激活态的ICL3、ECL2预测值波动最大，其次是ICL2，代表其结构的柔性最大易发生形变。由于缺失结构具有较大的柔性，B-factor值也相对较大致使结构不易测得，而在ANM的预测中邻近的残基也具有较高的B-factor值，根据趋势能推断出缺失结构同样具有较高的B-factor，因此缺失部位对模型的影响较小。在A_2A激活过程中其结构发生形变，激活态的ECL2波动与实验值较好的重合，而ICL3的B-factor值在整个结构中峰值最大，而在激活时ICL3的结构产生相应的偏移运动，这为使用ANM研究A_2AAR提供了有效性。

2.2 腺苷受体激活过程的动力学分析

为了探究A_2AAR激活过程中结合配体时受体蛋白的动力学影响，基于式(6)和式(10)计算了10~35 Å截断半径下A_2AAR激活态、中间态和未激活态下残基的均方涨落和变形能，通过归一化处理消除量纲，利用式(12)计算每个残基不同截断半径下的柔性量并取均值代表此残基在变构时的柔性量度，各个残基运动的柔性量度见图4(a)。结构越粗形变量越大，与受体激活过程中的功能性运动可能相关。由图4(a)可知，ECL的柔性较跨膜螺旋束而言更强，ECL在三种形态模型中的柔性都较大，根据三种形态的柔性量可明显发现ECL中的残基柔性量随着受体的激活不断改变并发生偏移，这是因为ECL具有受体变构调节剂可以结合的位点，当受体激活时胞外环会产生开口，利于配体结合，且ECL和TM顶部之间的二硫键有效地将ECL提供构象约束，二硫键网络形成一个刚性的开放结构，使配体结合腔暴露在溶剂中，允许小分子配体进入，受体与配体结合后的相互作用可使结构更加稳定^[6,22]。事实上A_2A腺苷受体的细胞外表面特征在很大程度上取决于其第二细胞外环(ECL2)，这与β₁AR、β₂AR和视紫红质有很大不同^[6,23]。观察图4的ICL3结构不难发现其柔性量较大，并且A_2AAR在RCSB PDB数据库中的ICL3结构未有完整的结构数据，在变构过程中ICL3结构逐渐弯曲变窄仍具有较大的柔性量，这是因为ICL3作为细胞内配体结合的控制开关，通过特定的相互作用促进G蛋白激活^[24]。在A_2AAR与CGS21680结合时(pdb代码：4UG2)，CGS21680中的苯乙胺基（2-羧乙基）的延伸在由跨膜TM2和TM7中，以及在ECL2和ECL3形成的延伸前庭中结合，与残基Glu169^ECL2、His264^ECL3、Leu267^7.31和Ile274^7.38的侧链形成范德华力作用，胺基与Ser67^2.64的侧链形成氢键，观察柔性量图4(a)可知，ECL3的柔性量随激活程度有所下降，这可能是配体结合后的氢键对结构起稳定作用，使其柔性有所损失，实验结果表明，ECL2和ECL3中His264^ECL3和Glu169^ECL2之间存在氢键且距离相对较近(见图4(b))，致使结构趋于稳定状态^[5]。螺旋束在未激活态与中间态均显示较强的刚性，而在激活态中螺旋束发生明显的位移形变，但其柔性没有显著增强，表明在激活过程中螺旋束发生整体性移动改变受体的结构实现结合配体。由实验结果可知，受体结合内源性激动剂会导致A_2A在变构过程中TM3、TM5、TM7向内运动发生形变，使配体结合口袋收缩。当核糖基团深入配体结合口袋后会与TM7中的保守残基Ser277^7.41和His278^7.42发生极性相互作用，以及与TM3中的残基产生范德华力，并导致TM5发生隆起使Cys185^5.46、Val186^5.47向内移动，α螺旋束的柔性也相应地产生损失^[25-26]。未激活态TM6的下端与TM7的上端柔性量较大，与其相连的ECL3与ICL3柔性量也较大，而在激活态处的柔性量则相对减弱，可知TM6、TM7参与激活过程，影响激活途径。根据实验结果可知，激动剂结合引起的构象变化，特别是TM7的旋转和向内移动以及TM6的向外移动，导致Na⁺结合口袋塌陷，增强Na⁺的位移，从而使激动剂克服Na⁺的负变构效应并激活受体^[27]，而短螺旋H8的柔性较大且发生移动，对受体激活时的作用还有待进一步研究。

2.3 快运动模式分析

在研究A_2AAR中对结构起稳定性作用的动力学关键残基时，快运动模式涨落幅度的极大值点残基可能是对构象变化起延缓作用，并维持结构稳定状态的重要残基^[28]。图5绘制了A_2AAR前两个快运动模式的均方涨落，其中黑色圈内的残基是已被实验确定的结合位点，而红色圈内的残基是模型预测的新的关键残基，这些残基可能是变构过程中对结构起稳定作用的关键残基。

由图5发现，TM1~TM7呈峰状，而ICL、ECL呈谷状，表明快运动模式能有效识别A_2AAR的结构特征。通过分析峰值区域的残基，整理可得到表2中的残基。由快运动模式识别A_2AAR的关键残基共有27个，其中有19个是已被实验证实的配体结合位点（在表2中白体显示的是残基），8个是模型预测的新关键残基（在表2中粗体显示的残基），可能是潜在的配体结合位点。这些关键残基大部分都是保守残基且邻近跨膜螺旋的中心处。缺失部位附近的残基涨落值较低，且相似结构同样具有较低的均方涨落，在快运动分析中需要识别较高涨落值的残基，故ICL3结构缺失对快运动模式分析影响较小，本小节不考虑ICL3结构处的关键残基。

在识别的这些关键残基中已有学者发现Na⁺与受体结合时会通过相互作用或间接作用达到结合的稳定状态，如Ser91^3.39的侧链氧原子会与Na⁺发生相互作用且Na⁺与Tpr246^6.48间存在水分子介导的氢键，在Na⁺结合过程中His278^7.42和Asn284^7.49之间的区域会经历向外扩张的位移运动，驱使TM7远离TM3，扩大了囊腔，利于钠离子结合^[29]。在变体结构实验中，Asn92^3.40的侧链可与Cys185^5.46的羰基形成氢键调节受体平衡^[30]。与激活相关的构象变化涉及“核糖袋”的扩大，是由残基Thr88^3.36、Ser277^7.41和His278^7.42决定，所有与A_2AAR共结晶的激动剂通过Thr88^3.36与TM3结合以及TM7的Ser277^7.41和His278^7.42结合，而这3个残基在快运动模式中均具有较大的涨落值^[24]。图5中Thr88^3.36位于波峰，对应的均方涨落值仅次于Ser91^3.39和Asn92^3.40，其对A_2AAR的结构至关重要。在突变实验的研究中Thr88A^3.36突变可将受体的活性状态转变为非活性状态，受体构象的稳定性也随之降低，而Ser281^7.46与Thr88^3.36都与钠离子配位的结构化水分子具有相互作用，对Thr88^3.36的研究值得进一步探索^[29,31]。研究表明，配体的核糖部分与Ser277^7.41、His278^7.42间的氢键相互作用可影响TM7的稳定性，这两个残基对激动剂的识别也至关重要^[26]，但在快运动模式中Ser277^7.41并不处于峰值，而His278^7.42处于峰值处，通过快运动模式的峰值识别关键残基具有一定的局限性，仍需结合其他方法识别关键残基。识别的非结合位点残基Tyr197^5.58在激活过程中其侧链采用新的旋转异构体来填补Val239^6.41的移动，稳定激活过程中的结构，且TM6移动时末端与螺旋TM5、TM7会产生明显的堆积变化^[24]。其余识别的非结合位点残基当前未发现其在变构过程中的作用，有待进一步深入研究。

2.4 慢运动模式分析

为了更深入了解激活A_2AAR过程中的功能性残基，计算了A_2AAR结构的前两个慢运动模式的均方涨落，结果如图6所示，黑色圈的残基是已被实验确定的结合位点，红色圈的残基是快运动模式识别的处于峰值处的非结合位点残基。图6中TM1和H8的涨落波动较大是因为其位于结构两端的位置，在模拟时其他结构对其约束相对较小就容易发生波动，在利用ANM分析其它具有N端与C端的GPCR时，TM1和H8的涨落一般较小。通过分析激活前、后三种蛋白质在慢运动模式的涨落值，发现受体的结合位点均具有较小的均方涨落，且配体结合位点在激活态中的涨落值相较于未激活态明显降低，而中间态与未激活态的涨落值比较接近且基本相似。缺失结构根据趋势具有较高的涨落值，故在慢运动模式分析中不用考虑缺失结构处的残基。

关键残基在慢运动模式中通常位于谷底处，观察图6的结合位点发现其大多数都处于谷底，个别位点在谷底残基的邻近处。如Asn24^1.50、Val27^1.53都处于谷底，在激活态的均方涨落高于另外两种状态，且在快运动模式中处于峰处。Ser281^7.46、Asn284^7.49、Pro285^7.50、Tyr288^7.53在中间态的慢运动模式涨落值最低，但Tyr288^7.53在快运动模式中并不处于峰处，事实上Tyr288^7.53位于NPxxY基序中，在结合激动剂的过程中与TM2、TM6和TM7中的残基产生相互作用，是重要的结合位点^[32]。Trp246^6.48在图6的中间态与激活态位于谷底，在快运动模式中也处于峰处，而Trp246^6.48位于保守的CWxP位点，即是Na⁺结合位点，又是配体结合位点，且是受体激活的旋转扭结开关并可以与进入受体的水分子产生相互作用^[2,33]。已有研究发现在激活过程中Ile92^3.40和Trp246^6.48、Cys185^5.46、Asn280^7.45之间的距离缩小，且它们之间存在氢键及弱相互作用，图6中激活态的这4个残基处于谷底处且涨落低于其他态^[30]。虽然胞外环在慢运动模式中总体的均方涨落较大，但重要位点残基的涨落仍在谷底，如Phe168^ECL2，激活过程中Phe168^ECL2的芳香族可与配体发生π-堆积相互作用或疏水接触^[6,34]。Asn253^6.55的涨落较其他配体结合位点稍大，但在局部仍处于谷底处，然而Phe168^ECL2和Asn253^6.55直接影响受体与激动剂或拮抗剂的结合，已作为重要的正构结合靶点设计药物开发^[31]。

综上所述，重要的残基位点在慢运动模式中具有较低的均方涨落，一般处于谷底且在激活态的涨落值低于未激活态或中间态。

2.5 互相关分析

为了直观地分析腺苷受体在发生功能性运动时残基间的相互联系，计算了各个残基间的交叉关联性系数，其反应运动过程中联动性较强的残基，结合GPCRdb数据库记录的残基相互作用散点图可知，A_2AAR各个结构之间产生相互作用的残基。分析残基间相互作用能判断受体结构中的重点残基，其对变构过程中的运动至关重要。从残基相关性热图(图7(a))可以看出相邻螺旋的相关性较高且每个螺旋中残基之间关联性很强，ECL2与ECL1的关联性也较强。根据实验结果可知，在受体激活过程中螺旋II、III和V位置会发生相应变化重新定义配体结合口袋的位置，A_2A腺苷受体的结合口袋移至靠近螺旋VI和VII结构^[6]。由图7(a)发现TM7与TM2、TM3具有较强的正相关性，TM5、TM6与TM3也具有较强的正相关性，所以螺旋区域的相互作用影响受体激活时各个部位的局部运动，如Arg102^3.50在稳定相邻Asp101^3.49的去质子化状态方面发挥作用，使该残基与TM2和ICL2形成更强的氢键相互作用，且在热图中Arg102^3.50与TM2和ICL2也呈现较强的正相关性^[6]。此外，在GPCRdb数据库记录的残基对间的相互作用(图7(b))可知各个结构与相邻结构存在相互作用的情况，结合相关性热图可发现相邻结构间的残基更容易产生相互作用，空间距离较近的残基也容易产生相互作用。TM1与TM7、H8相互作用较多且强，TM1与TM2之间、TM2与TM3之间、TM3与TM4之间存在较强的相互作用，图中红色区域较为明显。形成保守序列基序Asp/Glu-Arg-Tyr(D/ERY)的TM3(Asp101^3.49-Arg102^3.50-Tyr103^3.51)和TM6(Glu 228^6.30)的细胞质末端之间的相互作用构成一个“离子锁”，且它们呈现较强的负相关性，可能起到抑制受体激活的作用。

另外Glu169^ECL2和His264^ECL3这两个残基之间的相关系数为0.7，大于邻近残基的相关系数，在局部最大，而这两个残基通过形成的盐桥稳定结合口袋中的配体^[27]，能发现互相关性较强的残基对更容易产生相互作用，在GPCR变构过程中对配体的结合及结构的稳定起着重要作用，而负相关较强的残基则可能抑制受体的激活过程。

3 总结

研究结果揭示了A_2AAR激活过程中所涉及的影响蛋白质构象变化的重点残基及与功能性运动相关的关键残基。首先对A_2AAR构建了参数最优时的ANM模型预测B-factor，然后根据A_2AAR残基的均方涨落和变形能进行线性加权得到残基在变构时的柔性量。其次根据快慢运动模式识别腺苷受体变构时的关键残基分析了其在蛋白质结构中的分布及作用。最后计算了残基之间运动的交叉关联系数，结合残基对间存在的相互作用，识别对受体结构起稳健功能的重点残基。通过识别关键残基可以更好地了解蛋白质如何发生功能。本文研究有助于理解A_2AAR激活时的功能运动和变构动力学，根据关键残基设计特定药物干预或促进蛋白质实现生物功能，有助于蛋白质的功能改造和靶向药物设计。

参考文献

[1] WANG W J, QIAO Y H, LI Z J. New insights into modes of GPCR activation[J]. Trends in Pharmacological Sciences, 2018, 39(4): 367-386.

[2] CALIMAN A D, MIAO Y L, MCCAMMON J A. Mapping the allosteric sites of the A_2A adenosine receptor[J]. Chemical Biology & Drug Design, 2018, 91(1): 5-16.

[3] DE LERA RUIZ M, LIM Y H, ZHENG J Y. Adenosine A_2A receptor as a drug discovery target[J]. Journal of Medicinal Chemistry, 2014, 57(9): 3623-3650.

[4] MÜLLER C E, JACOBSON K A. Recent developments in adenosine receptor ligands and their potential as novel drugs[J]. Biochimica et Biophysica Acta (BBA)-Biomembranes, 2011, 1808(5): 1290-1308.

[5] LEBON G, EDWARDS P C, LESLIE A G W, et al. Molecular determinants of CGS21680 binding to the human adenosine A_2A receptor[J]. Molecular Pharmacology, 2015, 87(6): 907-915.

[6] JAAKOLA V P, GRIFFITH M T, HANSON M A, et al. The 2.6 angstrom crystal structure of a human A_2A adenosine receptor bound to an antagonist[J]. Science, 2008, 322(5905): 1211-1217.

[7] RODRÍGUEZ D, GAO Z G, MOSS S M, et al. Molecular docking screening using agonist-bound GPCR structures:probing the A_2A adenosine receptor[J]. Journal of Chemical Information and Modeling, 2015, 55(3): 550-563.

[8] MCROBB F M, NEGRI A, BEUMING T, et al. Molecular dynamics techniques for modeling G protein-coupled receptors[J]. Current Opinion in Pharmacology, 2016, 30: 69-75.

[9] CIANCETTA A, SABBADIN D, FEDERICO S, et al. Advances in computational techniques to study GPCR-Ligand recognition[J]. Trends in Pharmacological Sciences, 2015, 36(12): 878-890.

[10] MARTÍNEZ A M, CORREA B J. Molecular dynamics simulations reveal initial structural and dynamic features for the A₂AR as a result of ligand binding[J]. Molecular Simulation, 2014, 40(13): 996-1014.

[11] LOVERA S, CUZZOLIN A, KELM S, et al. Reconstruction of apo A_2A receptor activation pathways reveal ligand-competent intermediates and state-dependent cholesterol hotspots[J]. Scientific Reports, 2019, 9(1): 14199.

[12] BROOKS B, KARPLUS M. Harmonic dynamics of proteins: normal modes and fluctuations in bovine pancreatic trypsin inhibitor[J]. Proceedings of the National Academy of Sciences of the United States of America, 1983, 80(21): 6571-6575.

[13] TIRION M M. Large amplitude elastic motions in proteins from a Single-Parameter, atomic analysis[J]. Physical Review Letters, 1996, 77(9): 1905-1908.

[14] HALILOGLU T, BAHAR I, ERMAN B. Gaussian dynamics of folded proteins[J]. Physical Review Letters, 1997, 79(16): 3090-3093.

[15] ATILGAN A R, DURELL S R, JERNIGAN R L, et al. Anisotropy of fluctuation dynamics of proteins with an elastic network model[J]. Biophysical Journal, 2001, 80(1): 505-515.

[16] YANG L, SONG G, JERNIGAN R L. Comparisons of experimental and computed protein anisotropic temperature factors[J]. Proteins, 2009, 76(1): 164-175.

[17] 陈磊, 巩卫康, 李春华. 基于各向异性网络模型研究δ阿片受体的动力学与关键残基[J]. 生物化学与生物物理进展, 2022, 49(6): 1146-1154.

[18] SU J G, ZHANG X, ZHAO S X, et al. Conformational motions and functionally key residues for vitamin B₁₂ transporter BtuCD-BtuF revealed by elastic network model with a function-related internal coordinate[J]. International Journal of Molecular Sciences, 2015, 16(8): 17933-17951.

[19] BALLESTEROS J A, WEINSTEIN H. Integrated methods for the construction of three-dimensional models and computational probing of structure-function relations in G protein-coupled receptors[J]. Methods in Neurosciences, 1995, 25: 366-428.

[20] SUN Z T, LIU Q, QU G, et al. Utility of B-Factors in protein science:interpreting rigidity, flexibility, and internal motion and engineering thermostability[J]. Chemical Reviews, 2019, 119(3): 1626-1665.

[21] BAHAR I, JERNIGAN R L. Inter-residue potentials in globular proteins and the dominance of highly specific hydrophilic interactions at close separation[J]. Journal of Molecular Biology, 1997, 266(1): 195-214.

[22] WHEATLEY M, WOOTTEN D, CONNER M T, et al. Lifting the lid on GPCRs: the role of extracellular loops[J]. British Journal of Pharmacology, 2012, 165(6): 1688-1703.

[23] SCHOLL D J, WELLS J N. Serine and alanine mutagenesis of the nine native cysteine residues of the human A₁ adenosine receptor[J]. Biochemical Pharmacology, 2000, 60(11): 1647-1654.

[24] CARPENTER B, NEHMÉ R, WARNE T, et al. Structure of the adenosine A_2A receptor bound to an engineered G protein[J]. Nature, 2016, 536(7614): 104-107.

[25] LEBON G, WARNE T, EDWARDS P C, et al. Agonist-bound adenosine A_2A receptor structures reveal common features of GPCR activation[J]. Nature, 2011, 474(7352): 521-525.

[26] XU F, WU H X, KATRITCH V, et al. Structure of an agonist-bound human A_2A adenosine receptor[J]. Science, 2011, 332(6027): 322-327.

[27] CARPENTER B, LEBON G. Human adenosine A_2A receptor: molecular mechanism of ligand binding and activation[J]. Frontiers in Pharmacology, 2017, 8: 898.

[28] BAHAR I, DEMIREL M C, ERMAN B, et al. Vibrational dynamics of folded proteins: significance of slow and fast motions in relation to function and stability[J]. Physical Review Letters, 1998, 80: 2733-2736.

[29] GUTIÉRRÉZ D T H, MASSINK A, RODRÍGUEZ D, et al. The role of a sodium ion binding site in the allosteric modulation of the A_2A adenosine G protein-coupled receptor[J]. Structure, 2013, 21(12): 2175-2185.

[30] CUI M, ZHOU Q T, XU Y M, et al. Crystal structure of a constitutive active mutant of adenosine A_2A receptor[J]. IUCrJ, 2022, 9, Part 3: 333-341.

[31] JESPERS W, HEITMAN L H, IJZERMAN A P, et al. Deciphering conformational selectivity in the A_2A adenosine G protein-coupled receptor by free energy simulations[J]. PLoS Computational Biology, 2021, 17(11): e1009152.

[32] CALIMAN A D, SWIFT S E, WANG Y, et al. Investigation of the conformational dynamics of the apo A_2A adenosine receptor[J]. Protein Science, 2015, 24(6): 1004-1012.

[33] OLIVELLA M, CALTABIANO G, CORDOMÍ A. The role of Cysteine 6.47 in class A GPCRs[J]. BMC Structural Biology, 2013, 13: 3.

[34] WINK L H, BAKER D L, COLE J A, et al. A benchmark study of loop modeling methods applied to G protein-coupled receptors[J]. Journal of Computer-aided Molecular Design, 2019, 33(6): 573-595.

基金资助

内蒙古自治区自然科学基金项目(2022LHMS03014)

内蒙古自治区直属高校基本科研业务费项目(JY20220069)