基于各项异性网络模型研究μ阿片受体的动力学及关键位点

Study on the Dynamics and Key Sites of the μ-Opioid Receptor Based on the Anisotropic Network Model

摘要: 阿片类药物主要通过选择性作用于μ阿片受体(μ-opioid receptor, MOR)，激活抑制性G蛋白来实现镇痛的效果。由于其使用常伴随成瘾、呼吸抑制等严重的不良反应，针对阿片受体的药物研发成为一大挑战。为研究MOR的结构动力学和功能性关键位点，首先基于各向异性网络模型(Anisotropic network model, ANM)分析低频、高频运动模式下的残基涨落以及残基间的交叉相关涨落；然后，构造氨基酸复杂网络识别与变构信号传导相关的关键残基。结果表明，慢运动模式可有效识别结构中的功能区域，快运动模式能够识别对稳定MOR非活性结构重要的热点残基，运动相关性分析则揭示了跨膜螺旋区域变构信号的转导机制。另外，复杂网络识别出的关键残基在配体结合和信号传递中起重要作用。本工作有助于理解μ阿片受体的动力学特性和关键位点，并为药物设计提供有益的参考信息。

Abstract: Opioids mainly achieve analgesic effects by selectively acting on μ-opioid receptors (MOR) and activating inhibitory G proteins. Since its use is often accompanied by serious adverse reactions such as addiction and respiratory depression, the development of drugs targeting opioid receptors has become a major challenge. To investigate the structural dynamics and functional key sites of MOR, we first analyzed the residue fluctuations under low-frequency and high-frequency motion modes and the cross-correlation fluctuations between residues based on the Anisotropic network model (ANM). Then, a complex network of amino acids is constructed to identify key residues associated with allosteric signaling. The results showed that the lowest mode effectively identified functional regions within the structure, the fastest mode identified hotspot residues important for stabilizing the MOR inactive structure, and motion correlation analysis revealed the transduction mechanism of allosteric signals in the transmembrane helical region. Additionally, the complex network identified key residues that play important roles in ligand binding and signaling. This work helps to understand the dynamics characteristics and key sites of μ-opioid receptors, and provides useful reference information for drug design.

1. 引言

阿片受体(opioid receptors, ORs)负责调控镇痛、奖赏[1]以及情绪控制[2]等重要生理活动的信号通路。其中，μ阿片受体(μ-opioid receptor, MOR)是吗啡等阿片类镇痛药物的主要作用靶点[3]。阿片类药物和MOR结合可以产生止咳、镇痛的作用，但也会导致严重的不良反应，包括成瘾、呼吸抑制和便秘[4] [5]，从而限制了其临床应用。探究MOR的动力学特性，并明确参与变构信号传导的关键残基，对于了解其变构效应和开发镇痛药物具有重要的意义。

2022年，Robertson等人通过冷冻电镜技术解析了与拮抗剂Alvimopan结合的非活性MOR (残基65-352)结构[6]。与其他GPCRs类似，MOR的结构由7个跨膜(transmembrane, TM1-7) α螺旋、3个细胞内环(ECL1-3)、3个细胞外环(ICL1-3)以及1个C末端的第八螺旋(Helix 8, H8)组成，如图1所示。TM3通过C140^3.25 (上标表示Ballesteros-Weinstein编号)和C217之间的保守二硫键连接到ECL2。TM螺旋区域作为蛋白质的结构核心参与信号传导：当配体结合后会引起配体结合域(ligand-binding domain, LBD)微妙的构象重排，这种重排通过连接区域(connectorregion, CR)以变构和动态的方式触发长程构象转变，直至影响细胞内耦合域（intracellular coupling domains, ICDs)。随后，激活的ICD与信号伴侣蛋白(如G蛋白或β-抑制蛋白)结合，诱导ICD进一步大规模开放，从而达到完全激活状态[7]。

目前，许多计算模拟方法被应用于MOR变构机制的研究。Yuan等使用长时间尺度的分子动力学(molecular dynamics, MD)模拟探究了水分子和钠离子在MOR跨膜信号传导过程中的作用[8]。Kapoor等利用自适应采样模拟分析了MOR的激活态和失活态结构之间的配体特异性构象转变，同时，运用马尔科夫状态模型揭示了吗啡(平衡性激动剂)和TRV-130 (偏向性激动剂)对受体构象空间的不同调节方式[9]。Cong等采用核磁共振与溶质回火副本交换(replica exchange with solute tempering, REST2) MD模拟相结合的方法，通过比较不同激动剂对MOR构象动力学的影响，提出了MOR功能选择性可能的分子机制[10]。考虑到MD模拟的耗时性，人们发展了各向异性网络模型(Anisotropic network model, ANM) [11]，能够揭示与蛋白功能密切相关的动力学特性。另外，将复杂网络模型与ANM相结合可以用于研究蛋白质的变构通讯[12]，Del等发现网络中对特征路径长度(characteristic path length, CPL)影响大的残基通常对变构信号的传导至关重要[13]。

本工作基于ANM研究MOR的结构动力学与功能间的关系，并结合复杂网络模型识别蛋白质中与变构通信相关的关键位点。

Figure 1. Crystal structure of μ opioid receptor (MOR), PDB code: 7UL4

图1. μ阿片受体(MOR)的晶体结构，PDB代码为7UL4

2. 研究方法

2.1. 各向异性网络模型（ANM）

从蛋白质结构数据库 (Protein data bank, PDB) [14]中下载MOR的晶体结构(PDB ID: 7UL4)，其分辨率为2.8 Å。将受体的三维结构抽象为一个弹性网络，每个氨基酸简化为网络中的一个节点(通常用Cα原子表示)，若节点间的距离小于某一特定的截断半径，则用力常数统一的弹簧进行连接，整个体系的总势能为：

$V=\frac{1}{2}\gamma \left[\Delta R^{\text{T}}H\Delta R\right]$ (1)

其中，γ是弹簧的力常数，$\Delta R={\left[\Delta X_{\text{1}},\Delta Y_{\text{1}},\Delta Z_{\text{1}},\Delta X_{\text{2}},\Delta Y_{\text{2}},\Delta Z_{\text{2}},\cdots ,\Delta X_N,\Delta Y_N,\Delta Z_N\right]}^{\text{T}}$ ，N为残基数，上标T表示向量的转置，H表示3N × 3N的海森矩阵，其元素为：

$H=\left[\begin{array}{cccc}{H}_{11}& H_{12}& \cdots & H_{1N}\\ H_{21}& H_{22}& \cdots & H_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ H_{N1}& H_{N2}& \cdots & H_{NN}\end{array}\right]$ (2)

当$i\ne j$ 时，子矩阵H_ij为：

$H_{ij}=\left\{\begin{array}{l}\frac{{\partial }^{2}V}{\partial X_i\partial X_j}\frac{{\partial }^{2}V}{\partial X_i\partial Y_j}\frac{{\partial }^{2}V}{\partial X_i\partial Z_j}\\ \frac{{\partial }^{2}V}{\partial Y_i\partial X_j}\frac{{\partial }^{2}V}{\partial Y_i\partial Y_j}\frac{{\partial }^{2}V}{\partial Y_i\partial Z_j}\\ \frac{{\partial }^{2}V}{\partial Z_i\partial X_j}\frac{{\partial }^{2}V}{\partial Z_i\partial Y_j}\frac{{\partial }^{2}V}{\partial Z_i\partial Z_j}\end{array}\right\}$ (3)

当$i=j$ 时，子矩阵H_ii为：

$H_{ii}=-{\displaystyle \sum _{i\ne j}{H}_{ij}}$ (4)

海森矩阵分解得到的特征向量表示蛋白质的运动模式，残基的均方涨落(mean square fluctuations, MSF)和残基间的交叉相关涨落(correlated fluctuation, CoF)为：

$\langle \Delta R_i\cdot \Delta R_i\rangle =\frac{k_{B}T}{\gamma }\left(H_{3i-2,3i-2}^{-1}+H_{3i-1,3i-1}^{-1}+H_{3i,3i}^{-1}\right)$ (5)

$\langle \Delta R_i\cdot \Delta R_j\rangle =\frac{k_{B}T}{\gamma }\left(H_{3i-2,3j-2}^{-1}+H_{3i-1,3j-1}^{-1}+H_{3i,3j}^{-1}\right)$ (6)

其中，k_B是Boltzmann常数，T为绝对温度。归一化的运动交叉相关性为：

$C_{ij}=\frac{\langle \Delta R_i\cdot \Delta R_j\rangle }{\sqrt{\langle {\left(\Delta R_i\right)}^2\rangle \cdot \langle {\left(\Delta R_j\right)}^2\rangle }}$ (7)

C_ij的取值范围在[−1, 1]之间，其中正值表示残基间的运动方向相同，而负值则表示运动方向相反。C_ij的绝对值越大，意味着残基间的运动相关性越强。

2.2. 氨基酸复杂网络模型

蛋白质的结构和功能依赖于氨基酸残基之间复杂的相互作用网络。在加权的复杂网络模型[12]中，用残基中的Cα原子代表网络中的节点，截断半径内的节点用边连接，边权重基于ANM得到的残基间的运动相关性计算得出：

$w_{ij}=-\log \left(\left|C_{ij}\right|\right)$ (8)

特征路径长度(Characteristic Path Length, CPL)是网络中所有节点对之间最短路径的平均长度：

$CPL=\frac{1}{N_p}{\displaystyle \sum _{j>i}^N{d}_{ij}}$ (9)

其中N和N_p分别是节点数和节点对数，d_{i j}表示从节点i到j累计边权重最小的路径长度。节点k对网络内信息通信的贡献是通过从网络中移除节点k后CPL的变化(ΔCPL_k)来衡量[15]，用Z-score表示：

$Z{\text{-score}}_k=\frac{\Delta CP{L}_k-\langle \Delta CP{L}_k\rangle }{\sigma }$ (10)

其中，ΔCPL_k是移除节点k后特征路径长度的变化，$\langle \Delta CPL\rangle$ 是依次移除各节点后特征路径长度变化的平均值，σ是标准差。

3. 结果与讨论

3.1. 温度因子的理论值与实验值比较

温度因子(B-factor)可以体现分子的局部柔性。对MOR构建模型时，通过最大化理论B-factor和实验之间的皮尔森相关系数(Pearson correlation coefficient, PCC) [16]来优化ANM中的截断半径，其寻优范围是：6.0~30.0 Å，步长为1.0 Å。当r_c = 23.0 Å时，理论B-factor和实验的PCC值达到最大，为0.79 (图2)。ANM可以较好地模拟MOR中各结构域的波动。TM区域的涨落幅度较小，而ECLs和ICLs区域的涨落则更为明显，符合跨膜结构的特点。另外，实验发现ECLs负责多种配体的识别，ICLs可以与G蛋白等效应蛋白相互作用，这些区域的柔性变化在信号转导中至关重要[17] [18]。

Figure 2. Comparison between the theoretical B-factors (dotted line) and experimental (solid line) of MOR

图2. MOR的理论B-factors(虚线) 和实验 (实线)的比较

3.2. 慢运动模式分析

慢运动模式通常代表与蛋白功能相关的大幅度集合运动，图3显示了在第一慢运动模式下MOR中各残基的涨落。峰值对应重要的功能区域——ECL2和ICL3，ECL2的β发夹结构在配体结合时会出现轻微的移动，该区域残基的突变对配体结合和受体激活有不利影响[19]；ICL3通过阻断或暴露G蛋白结合位点的动态构象平衡来调节受体活性[20]。MOR中涨落较小的区域同样具有特定的功能，第1个区域中的I77^1.41和V80^1.44以及第2个区域中的P122^2.58参与构成MOR的二聚体界面；位于TM3的第3个区域包含保守的DRY基序[Asp164^3.49-Arg165^3.50-Tyr166^3.51]，基序中的局部螺旋内接触以及和T279^6.34之间的极性相互作用在阿片受体中普遍存在，研究表明消除残基间的极性相互作用会使受体具有组成性活性[21]。

Figure 3. Mean square fluctuation distribution of the residues in the first slow motion mode

图3. 第一慢运动模式下残基的均方涨落分布

3.3. 快运动模式分析

在快运动模式下，表现出较大运动幅度的残基往往被视为动力学上的热点残基，这些残基在稳定蛋白质空间结构中起着至关重要的作用[22]。图4中显示了5个涨落峰值，对应的残基分别是：D114^2.50、S154^3.39、L158^3.43、V282^6.37和S329^7.46。其中，D114^2.50、S154^3.39、S329^7.46和N150^3.35、N328^7.45组成钠离子结合位点，并且形成广泛的极性相互作用网络，钠离子结合到该位点，通过变构效应可以稳定受体的非活性构象。将D114^2.50位点上的带电残基替换为不带电残基后，一些GPCRs的激动剂依赖性信号传导能力会显著下降，而配体结合能力和基础信号传导水平则得以保持，揭示了该位点在激动剂引起的GPCR信号传递过程中的特定作用[23]。保守的疏水残基L158^3.43和V286^6.41组成疏水锁，其闭合会压缩TM6和TM3的堆积，有利于TM6细胞内端的向内移动，这是GPCR失活所必需的[24]。V282^6.37可以与G蛋白的Gα亚基上的保守残基相互作用，在MOR激活过程中，胞质侧偶联区域残基的构象变化会提高受体对激动剂的结合亲和力[25]。

Figure 4. Mean square fluctuation distribution of the residues in the first fast motion mode

图4. 第一快运动模式下残基的均方涨落分布

3.4. 残基运动相关性分析

为研究MOR中残基之间的功能性运动耦合，基于公式(7)计算了残基涨落之间的动态互相关矩阵 (Dynamic cross-correlation matrix, DCCM)，图5是前208个低频运动模式(对残基涨落贡献大于50%)获得的结果。涨落之间的交叉相关性能够提供与空间运动和相互作用相关的信息，图中橙红色区域表示相邻TM螺旋以及空间上紧邻的TM螺旋之间具有一定的相互作用。TM2、TM3、TM6和TM7螺旋内以及各自和其余螺旋之间存在强烈的正相关，归因于配体结合在这四个螺旋的胞外侧，进而触发了螺旋之间的相互作用和构象变化。具体的说，TM6和TM2、TM3、以及TM7呈正相关，对应受体的“传输开关”W293^6.48和钠离子结合位点(D114^2.50、N150^3.35、S154^3.39、N328^7.45和S329^7.46)的相互作用，实验表明保守残基W293^6.48的侧链朝向变化会导致钠离子结合口袋的重排[24]。TM2、TM3和TM7上部分残基通过水分子形成极性网络，稳定了MOR的非活性构象[25]。其次TM5和TM3、TM6之间有较强的正相关，是由于W293^6.48和关键的PIF基序[Pro244^5.50-Ile155^3.40-Phe289^6.44]共同构成所谓的CR，转导了LBD和ICD之间的变构通讯[26]。

3.5. 关键位点识别

接收和传递变构信号的重要残基通常是蛋白质相互作用网络的中心，位于大多数残基之间的最短路径上。通过依次移除氨基酸网络中的节点及其相连通路，计算特征路径长度变化的Z-score分数，来量化

Figure 5. Dynamic cross-correlation map between MOR residues

图5. MOR残基之间的动态互相关图

残基对网络内通信的影响，结果如图6(a)所示，共识别出22个关键残基簇。对氨基酸网络影响比较大的残基主要分布在TM2、TM3、TM6和TM7上，表明这些区域在变构通信中起重要作用，与运动相关性分析的结果(图5)一致。将Z-score分数映射到MOR结构上，并标注值高的簇对应的中心残基(Z-score ≥ 1.5，图6(b))。按照信号从胞外侧向胞质侧传递的顺序逐个分析：W318^7.35参与形成疏水口袋，很大程度上有助于与配体的相互作用。Y166^3.51通过与其它残基形成广泛的氢键网络来稳定TM3和TM6之间的连接，从而稳定MOR的非活性结构[25]。A337^7.54与K344^8.51的螺旋内接触则有利于维持TM7的外移[27]。R277^6.32与R280^6.35在G蛋白偶联区域，该区域与配体结合口袋存在弱偶联，被认为是潜在的变构调节剂结合口袋[28]。对于残基M90^1.54，尚未有相关研究发现它在变构通信中的作用，其潜在的功能还需实验证实。

Figure 6. (a) The Z-score value for the difference in CPL. (b) MOR structures are color-coded by Z-score values

图6. (a) 特征路径长度变化的Z-score分数；(b) MOR结构用Z-score分数进行颜色编码

4. 结论

本工作研究了μ阿片受体 (MOR) 的动力学特性，以及对蛋白变构信号转导起重要作用的残基。在最佳参数下构建的受体的ANM模型，计算得到的B-factor与实验值的PCC达到0.79，较好地再现了各结构域的柔性，证明ANM可用于研究MOR的结构动力学与功能的关系。在慢运动模式下，ECL2和ICL3有大规模的集合运动，猜测ECL2的移动有助于配体结合口袋的打开，ICL3的构象变化负责调节受体活性，涨落比较小的残基与MOR的二聚化、激活密切相关。快运动模式下识别出了可以稳定受体内外两侧非活性结构的关键残基。运动相关性分析表明变构信号主要是通过TM2、TM3、TM5、TM6和TM7的耦合构象变化传导的。最后，加权的氨基酸网络模型预测了MOR的功能性关键残基，大部分结果与实验数据吻合，个别残基还未被证实，可为未来的实验研究提供参考。本工作有助于理解MOR的变构通讯机制，对相关药物设计具有指导意义。