基于各向异性网络模型研究δ阿片受体的动力学与关键残基*

陈 磊 巩卫康 李春华

（北京工业大学环境与生命学部，北京 100124）

G蛋白偶联受体（GPCRs）是最大的膜蛋白家族，调节许多生理过程中的信号通路，如行为、认知和免疫反应［1］。它们在人体生理方面的中心调节作用使其成为关键的药理学靶点［2］。阿片受体是一种A 类GPCR，包括μ（MOP）、δ（DOP）、κ（KOP）和孤啡肽FQ（OFQ）成员［3］，与疼痛控制相关。其中，DOP 还与情绪控制有关，因此作用于DOP 的药物显示出额外的抗焦虑和抗抑郁作用［4］。DOP 的结构和变构动力学引起了广泛的关注。

在实验研究方面，2012年，Granier 等［5］利用X射线去探索阿片类配体识别的保守片段，揭示了对配体亚型选择性重要的结构特征。2014年，Fenalti等［6］以1.8 Å的高分辨率解析了人源DOP与拮抗剂naltrindole7 的复合物结构。2019年，Claff等［4］研究了与DOP 结合的激动剂，并确定了激动剂结合以及受体激活的关键决定因素。在理论方面，Shang 等［7］使 用 分 子 动 力 学（molecular dynamics，MD）模拟研究了正构激动剂SNC‑80存在的情况下变构剂BMS‑986187 与DOP 的结合过程。利用多尺度模拟，Wang 等［8］分别研究了MOP 和DOP 二聚反应及其激活过程中的协同机制。此外，对于重要的变构调节剂钠离子，Shang等［3］利用全原子MD 模拟观察到钠离子与MOP、KOP 和DOP 的结合途径类似，还发现钠离子的结合降低了特异性激动剂的结合水平。值得注意的是，结合MD 模拟和实验，Sun 等［9］探索了DOP中钠离子的变构调节机制，揭示了钠离子利用残基Trp2746.48的独特构象将信号传递到TM5和TM6。

MD模拟是一种耗时的方法，特别是对于生物大分子。为了解决这个问题，科学家们提出粗粒化模型。其中，弹性网络模型（elastic network model，ENM）是研究蛋白质内在动力学和功能相关运动的一种特别有效的模型［10］。高斯网络模型（Gaussian network model，GNM）［11］和各向异性网络模型（anisotropic network model，ANM）［12］是两种广泛使用的ENM，它们可以计算蛋白质结构所具有的动力学信息。相比而言，ANM 包含了残基运动的方向性，所以能够揭示的动力学信息更多。通常，ENM 获得的低频运动模式代表与蛋白质功能相关的大规模集体运动，而高频模式反映蛋白质结构的几何不规则性，其下的高涨落残基被认为是对蛋白质稳定性重要的残基［13］。由于ENM的高效性，基于它的微扰响应方法被发展用于研究蛋白质的变构特性。2009年，Atilgan等［14］提出了微扰响应扫描（perturbation‑response scanning，PRS）方法探索蛋白质残基对外界微扰的响应程度，目前该方法已被广泛应用于识别变构调控中的关键残基［15］。

本研究基于ANM探索了DOP的动力学与功能特性。除此之外，结合PRS方法识别出对蛋白质的变构通信具有重要作用的关键残基。

Fig.1 Crystal structure of δ opioid receptor（DOP）（PDB ID：4N6H）The seven transmembrane helixes (TM1-7) are connected by 3 extracellular loops(ECL1-3)and 3 intracellular loops(ICL1-3).

1 研究体系与方法

1.1 研究体系

人源δ阿片受体（DOP）结构是由X射线衍射实验解析出来的（PDB ID：4N6H）（图1），共有303 个残基，包括N 端（残基36～38），7 个跨膜（transmembrane helixes，TM）螺旋（残基39～77、82～112、117～152、161～187、205～243、249～287和293～321）、3 个 胞 内 环（intracellular loops，ICLs）（残基78～81、153～160 和244～248）、3 个胞外环（extracellular loops，ECLs）（残基113～116、188～204和288～292）、1个螺旋H8（残基322～335）以及C端（残基336～338）。

这项工作采用了Ballesteros和Weinstein定义的GPCRs 残基编号方案［16］。简单地说，每个残基的上标都被标记为X.YY，其中数字X 对应于残基存在的TM 螺旋（从1 到7），YY 为相对于螺旋中最保守的残基（YY 值为50）的位置。例如，如果TM2 中最保守的残基是一个天冬氨酸（残基号95），则其标识符为2.50，即Asp952.50。对于一个残基Asn3107.45，表示在TM7 中最保守残基Pro3157.50之前的第5个残基处。

1.2 各向异性网络模型（ANM）

在ANM 中，蛋白质被抽象为一个弹性网络，其中氨基酸用一个节点来表示（通常选用Cα原子作为节点），当两个残基节点间的距离小于某一截断半径时（本文选取17 Å），它们之间用一个弹簧相连，所有弹簧的弹性系数相同。在ANM中，整个体系的势能函数为：

其中，Rij和Rij0分别是节点i 和j 之间的瞬时和平衡距离，γ是弹簧的弹性系数。蛋白质的运动模式是由Hessian 矩阵H 所决定，矩阵的各元素为体系势函数对位置的二阶偏导，具体形式如下：

其中，每一个元素hij是3×3 的子矩阵。当i ≠j 时，hij可以表示为：

当i=j时，hii可以表示为：

对Hessian矩阵进行分解可获得3N个本征值和本征向量。本征值小的本征向量对应蛋白质的慢运动模式，与功能有关；本征值大的对应蛋白质的快运动模式，与结构稳定有关。

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

其中，λ-1k是H 矩阵的第k 个本征值的倒数，[uk]3i是H矩阵的第k个本征向量中的第3i个元素。

归一化的交叉相关性为：

该值的取值范围为-1～1。正负值分别表示残基沿相同和相反的方向运动。绝对值越高，两个残基的相关性越强。零值意味着残基的运动是完全不相关的。

1.3 微扰响应扫描模型（PRS）

线性响应理论（LRT）是一种线性扰动蛋白质残基后，预测构象变化的方法［17］。基于LRT 的微扰响应扫描（PRS）方法被提出［15］，用于研究蛋白质的变构调控并识别关键残基。基于胡克定律F= H ΔR，在受到外力扰动时，残基位置的改变量为ΔR，其中H 是ANM 中3N × 3N 的Hessian 矩阵。当对网络中的一个残基i施加外力时：

观测到其他残基的响应：

其中，ΔRi是在外力对残基i 扰动后，所有残基位置的改变量。

对蛋白质的残基i 施加一个力Fi，可以得到1个其他残基产生的3N维的位移响应向量ΔRi，依次扰动所有残基可以得到位移响应矩阵ΔR。具体实施中，在每个残基上施加7个单位向量的外力，即方向为（1，0，0）、（0，1，0）、（0，0，1）、（1，1，0）、（1，0，1）、（0，1，1）和（1，1，1），之后对7 次微扰得到的7 个响应矩阵ΔR 求平均，然后对这个平均矩阵每一行中的元素相对于对角线元素进行归一化，得到一个N×N 阶的响应矩阵P，对其数字进行颜色编码即为热图，热图中第ij个元素表示微扰节点i，节点j产生的平均响应。归一化后的P矩阵中每一列的平均值代表了该列对应残基的敏感性，即残基响应其他微扰残基的倾向性。所有残基的敏感性就组成了敏感性特征曲线。该曲线峰值对应的残基被认为是蛋白质变构中的功能性运动位点。相似的，每一行的平均值代表了该行对应残基的效应性，即衡量微扰残基影响其他残基动力学的能力。相应曲线的峰值被认为对信号传输起重要作用的热点残基［18］。

基于上述方法的描述，本文使用MATLAB R2019a 对ANM 和PRS 的相关程序进行编写并在DOP这个体系上进行研究。

2 结果与讨论

2.1 基于ANM的实验与理论B-factor的比较

温度因子（B‑factor）可以反映生物大分子的局部结构的柔性。本文构建了DOP 的ANM 模型，选用残基的Cα原子作为节点。传统的ANM只有一个截断半径参数，对截断半径的选取方法分为使用经验值［19］和在一定范围内通过最大化理论值与实验值的相关性进行寻优［12］。后一种方法能够提高ANM 构建蛋白质模型的准确性，所以在这个工作中，采用寻优的方法对截断半径进行优化，其寻优范围为：6～24 Å，步长为1 Å。在截断半径为17 Å时，ANM 计算的理论B‑factor 与实验B‑factor 的皮尔森相关系数（Pearson’s correlation coefficient，PCC）可以达到最大值（0.64），因此，本研究选择17 Å。图2 显示了在最优的参数下ANM 得到的理论和实验B‑factor 比较图。从图2 可以看出，ANM 可以较好地再现DOP 的柔性信息，TM 区域柔性最小，ECL 和ICL 柔性最大，这与实验结果一致［20］。

Fig.2 Comparison between the theoretical（dotted line）and experimental（solid line）B-factor of DOP

2.2 基于ANM的DOP慢运动模式分析

全局慢运动模式代表大范围的集体运动，通常与功能有关［13］。为研究DOP 的功能动力学特性，计算了DOP在前两个慢运动模式下的涨落（图3a，c）。7个局部最小值分别对应于TM1～7，峰值对应于ECL和ICL区域，中间区域的局部最小值对应于ECL2 中的β 片层，表明ANM 对DOP 的结构有较好的识别能力。

有趣的是，发现其中几个涨落较小的区域与DOP的功能相关（图3a中圆圈）。第1个区域（残基82～97）位于TM2，具有一定的保守性，其中Asp952.50是一个关键的钠结合位点，实验发现其突变可以破坏“钠效应”。第2 个区域（残基134～144）位于TM3，其中残基Ser1353.39、Asn1313.35和Asp952.5形成了第1 个钠离子结合位点。位于TM7中的第3 个区域（残基301～320）也在钠离子转运中起着重要作用，其中残基Asn3107.45和Asn3117.46形成了钠离子的第2个配位壳，位于NP7.50xxY基序的残基Asn3147.49和Tyr3187.53可以调节变构钠离子结合，与GPCR的激活有关［21］。

2.3 基于ANM的DOP的快运动模式分析

快运动模式对应的是蛋白质局部结构的不规则性［22］。在快运动模式下，运动幅度较大的残基通常是热点残基，对稳定蛋白质结构至关重要［23］。图3b 显示了前两个快运动模式的涨落分布图。由图3b 可以看出，11 个峰值区域对应的中心残基分别 是 Asn671.50、 Ile882.43、 Leu912.46、 Asp952.50、Asp1283.32、 Asn1313.35、 Ser1353.39、 Val2666.40、Yyr3087.43、Ser3117.46、Asn3147.49。这些残基分布在DOP跨膜螺旋的中心区域。

接下来，将根据现有的实验和理论数据，对识别出的关键残基功能进行讨论。高度保守的残基Asn671.50是DOP I‑I 二聚体的一个结合热点，对于稳定DOP的二级结构具有重要作用［24］。在钠离子转运的过程中，它除了与两个钠离子配位壳残基结合之外，还与高度保守的残基如Leu912.46形成氢键［25］。残 基Asp952.50、Asn1313.35和Ser1353.39与 两个保守的水分子组成钠离子的第1个结合口袋［21］。残基Asp1283.32是DOP 与抑制剂naltrindole7 的结合位点［5］。在激动剂DPI‑287 与DOP 对接的过程中，Yyr3087.43与DPI‑287形成氢键，使复合物的结构更稳定［26］。Ser3117.46和Asn3147.49是第2 个钠离子配位口袋处残基，在DOP 激活的过程中，Asn3147.49的位移会导致活性GPCR 的变构钠结合口袋崩塌［26］。Ile882.43和Val2666.40分别位于TM2 和TM6，目前还没有关于它们对DOP 变构作用的研究，值得被进一步探讨。

综上所述，识别的关键残基在快运动模式下具有高度的活性，在稳定DOP 结构方面发挥着重要作用或者在DOP与配体相互作用中扮演重要角色。

2.4 运动相关性分析

为了研究蛋白质运动过程中残基之间的关联性，通过式（7）计算了所有残基之间的交叉相关系数，并用二维彩图形式显示出来（图4a）。因为低频运动模式常常与蛋白质功能有关，选择对残基涨落贡献刚超过50%的低频运动模式进行运动相关性的分析是一种通常的做法，可以提高信噪比［27］，所以本文选取对涨落贡献刚大于50%的运动模式计算运动相关性。从图4a 中可以看出，对角线区域被分成7部分，垂直于对角线的部分显示相邻TM 螺旋之间具有强的正相关性，表明ANM可以对DOP 的结构具有较好的识别能力。除此之外，胞外/胞内跨膜螺旋与环状区域之间也具有一定的正相互作用（图中橙红色区域），可能是因为DOP 的跨膜螺旋部分区域在空间上是相邻的。特别是位于TM2 的残基Asp952.50、位于TM3 的Asn1313.35和Ser1353.39与 位 于TM6 和TM7 的Trp2746.48、Ser3117.46和Asn3147.49呈较强的正相关，研究表明这些残基位于两个钠离子结合口袋处［21］。在DOP 非活性构象中，Leu246ICL3和Val243ICL3与Val1503.54形 成 疏 水 簇。ICL3 还 通 过Leu246ICL3和Val1503.54和Arg2395.66之间的水介导氢键网络与TM3 相互作用。这些结果表明，ICL3、TM3 和TM5 存在相互作用，进一步稳定蛋白质结构［21］。Yyr1293.33、Met1323.36、Ile3047.39和Yyr3087.43形成一个结合口袋，它与DPI‑287 形成了疏水相互作用。Lys2145.39、His2786.52和Yyr1293.33与KGCHM07形成一个极性相互作用网络［26］。总之，DOP 的跨膜螺旋之间的相互运动不仅可以稳定构象状态，还能促进与配体之间的结合。

Fig.3 Residue mean square fluctuation （MSF）profiles for the DOP(a) Fluctuation distribution of the first two slow motion modes. (b) The fluctuation distribution of the first two fast motion modes. (c) First two slowest motion modes (depicted with a cone model) mapped on the DOP structure. The cone’s length is proportional to the motion magnitude, and the cone’s orientation indicates the motion direction.

2.5 微扰响应扫描分析

使用PRS 方法研究DOP 中变构信号的可能途径。该方法是基于LRT，已成功地用于揭示变构信号转导过程中的关键区域［15］。本文对DOP中的残基逐个微扰可以得到一个响应矩阵P（研究体系与方法1.3）。使用矩阵P 直接可以画出DOP 的PRS响应热图（图4b）。在PRS 图谱中，第ij 个元素表示是扰动残基i 对j 产生的影响，它衡量的是蛋白质结构的敏感性和效应性。

经过分析，发现DOP 中具有大的敏感性残基（变构信号的潜在接收器）主要分布在胞外环和胞内环区域（图4c），然而具有大的效应性残基（传播变构信号）主要分布在跨膜螺旋区域（图4d）。图4b的顶部是敏感性曲线图。敏感性高的残基有6簇（图4c），中心残基是Thr78ICL1、Glu1122.67、Asp158ICL2、 Asp193ECL2、 Arg244ICL3、 Asp290ECL3。这些残基在前面的慢运动模式涨落分析中位于大的柔性区域，因此对于蛋白质的功能性运动具有一定作用。其中残基Lys792.49在形成DOP 的I‑I 二聚体时发挥作用［24］。在正变构调节剂BMS‑986187 与DOP结合过程中，残基Glu1122.67与BMS‑986187形成疏水相互作用［28］。口袋残基Trp114ECL1、Ile289ECL3和Arg291ECL3分别与配体KGCH07 形成ππ、疏水和阳离子‑π 相互作用，可以进一步稳定DOP的激活状态［26］。由于ECL2是β链折叠，是所有阿片受体亚型的典型［29］，负责识别多种配体，因此位于此处的残基如Asp193ECL2可能参与配体的结合。残基Arg244ICL3通过与TM6和TM7的其他残基形成广泛的氢键网络，在稳定ICL3 方面起着关键作用，使DOP 稳定在非活性状态［21］。在研究钠离子与DOP 结合的MD 模拟实验中发现，Asp293ECL3周围的钠密度较高，因此这个位点可能与钠离子结合有关［25］。

图4b 的右边是效应性曲线图。效应性高的残基有12 簇（图4d），中心残基分别是Ser421.25、Gly631.46、 Asp952.50、 Ala1172.53、 Asn1313.35、Cys1513.36、 Asp2105.35、 Met2365.61、 Leu2405.67、Leu246ICL3、Arg2576.31和Ser3117.46。残基Ser421.25是DOP 形成A‑I 二聚体的界面残基，残基Gly631.46是DOP 形成I‑I 二聚体时的界面残基，这些残基可以促进二聚体的形成［24］。在这些效应残基中，Asp952.50和Asn1313.35是组成第1 个钠离子结合口袋的残基，Asn3107.45是组成第2个钠离子结合口袋的残基［21］。Ala1172.53是激活状态下的MOP 单体形成A‑I 二聚体时的界面残基［24］。MD 模拟研究发现，Cys1513.36对阿片受体与芬太尼的识别和相互作用至关重要［30］。Asp2105.35可以与激动剂KGCHM07形成水介导的盐桥相互作用，稳定DOP 的激活构象［26］。Leu246ICL3位于双亮氨酸基序（Leu245ICL3～Leu246ICL3）中，实验发现基序的缺失或Leu245 的突变会减缓DOP 的溶酶体靶向性［31］。残基Arg2576.31与Leu2405.67、Arg244ICL3和Val243ICL3形成氢键网络，与Asp2536.27形成盐桥，它连接TM5 和TM6的胞内结构，将DOP稳定在非活性状态［21］。

Fig.4 Dynamics identified by ANM and PRS for DOP(a)The cross‑correlation map for DOP.The regions marked by black oval box in the figure are the parts of key analysis.(b)PRS heat map for DOP.The strongest disturbance (dark red) is shown in the figure, which describes the dynamic influence of disturbance residue i on residue j. Sensitivity and effectiveness are shown along the top and right,and their values are mapped to(c)and(d),respectively.The highest value on the sensitivity curve corresponds to the strongest sensor,while the highest value of effectiveness corresponds to the strongest effector.

进一步分析发现效应性高的残基中也有部分残基处在柔性较高的区域，这些区域与敏感性高的区域共同发挥功能性作用，如效应性强的残基Ser421.25、Gly631.46和敏感性强的残基Lys792.49共同耦合促进DOP二聚体的形成［24］。效应性强的残基Asp952.50、Asn1313.35、Ser3117.46和敏感性强的残基Asp293ECL3共同作用组成DOP 中钠离子结合口袋［21］。

基于以上分析发现，通过PRS 确定的关键残基在DOP 的二聚化、功能状态稳定以及与变构配体（如钠离子和激动剂）的相互作用中起着重要作用，这些都参与了DOP的变构调节。

3 结 论

本工作探索了人DOP 的动力学特性，以及对蛋白质变构具有重要作用的关键残基。首先，根据DOP的结构构建出最优参数下的ANM，实验与理论B‑factor的PCC值达到0.64，并且识别出了DOP中的柔性区域，表明ANM可以较好的研究DOP的动力学特性。然后，从慢运动模式的涨落图可以很好地识别DOP 的结构域，包括7 个TM 螺旋、3 个ICL 和3 个ECL。有趣的是，发现钠离子结合口袋处具有最小的运动性，这有助于分析钠离子的变构调节。此外，从快运动模式的涨落曲线来看，模式中的活性位点主要分布在TM1～3、TM6 和TM7中，这些区域在稳定DOP 结构或与变构配体（如拮抗剂和钠离子）的复合物结构方面起着重要作用。接着对蛋白质的运动相关性的分析可以看出胞外/胞内跨膜螺旋与环状区域之间存在较强的相互作用，这有助于蛋白质的构象稳定以及与配体的结合。最后采用PRS 方法对DOP 进行微扰，识别出的关键残基已经在实验或理论上发现直接参与DOP 二聚化和与配体结合，认为这些关键残基对于DOP的变构通信至关重要。这项工作有助于对δ阿片受体变构调节的物理机制的理解，并为药物设计提供有价值的信息。