基于各向異性網(wǎng)絡模型研究δ阿片受體的動力學與關鍵殘基*

陳 磊 鞏衛(wèi)康 李春華

（北京工業(yè)大學環(huán)境與生命學部，北京 100124）

G蛋白偶聯(lián)受體（GPCRs）是最大的膜蛋白家族，調節(jié)許多生理過程中的信號通路，如行為、認知和免疫反應［1］。它們在人體生理方面的中心調節(jié)作用使其成為關鍵的藥理學靶點［2］。阿片受體是一種A 類GPCR，包括μ（MOP）、δ（DOP）、κ（KOP）和孤啡肽FQ（OFQ）成員［3］，與疼痛控制相關。其中，DOP 還與情緒控制有關，因此作用于DOP 的藥物顯示出額外的抗焦慮和抗抑郁作用［4］。DOP 的結構和變構動力學引起了廣泛的關注。

在實驗研究方面，2012年，Granier 等［5］利用X射線去探索阿片類配體識別的保守片段，揭示了對配體亞型選擇性重要的結構特征。2014年，F(xiàn)enalti等［6］以1.8 ?的高分辨率解析了人源DOP與拮抗劑naltrindole7 的復合物結構。2019年，Claff等［4］研究了與DOP 結合的激動劑，并確定了激動劑結合以及受體激活的關鍵決定因素。在理論方面，Shang 等［7］使 用 分 子 動 力 學（molecular dynamics，MD）模擬研究了正構激動劑SNC?80存在的情況下變構劑BMS?986187 與DOP 的結合過程。利用多尺度模擬，Wang 等［8］分別研究了MOP 和DOP 二聚反應及其激活過程中的協(xié)同機制。此外，對于重要的變構調節(jié)劑鈉離子，Shang等［3］利用全原子MD 模擬觀察到鈉離子與MOP、KOP 和DOP 的結合途徑類似，還發(fā)現(xiàn)鈉離子的結合降低了特異性激動劑的結合水平。值得注意的是，結合MD 模擬和實驗，Sun 等［9］探索了DOP中鈉離子的變構調節(jié)機制，揭示了鈉離子利用殘基Trp2746.48的獨特構象將信號傳遞到TM5和TM6。

MD模擬是一種耗時的方法，特別是對于生物大分子。為了解決這個問題，科學家們提出粗?；Ｐ?。其中，彈性網(wǎng)絡模型（elastic network model，ENM）是研究蛋白質內(nèi)在動力學和功能相關運動的一種特別有效的模型［10］。高斯網(wǎng)絡模型（Gaussian network model，GNM）［11］和各向異性網(wǎng)絡模型（anisotropic network model，ANM）［12］是兩種廣泛使用的ENM，它們可以計算蛋白質結構所具有的動力學信息。相比而言，ANM 包含了殘基運動的方向性，所以能夠揭示的動力學信息更多。通常，ENM 獲得的低頻運動模式代表與蛋白質功能相關的大規(guī)模集體運動，而高頻模式反映蛋白質結構的幾何不規(guī)則性，其下的高漲落殘基被認為是對蛋白質穩(wěn)定性重要的殘基［13］。由于ENM的高效性，基于它的微擾響應方法被發(fā)展用于研究蛋白質的變構特性。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 個 胞 內(nèi) 環(huán)（intracellular loops，ICLs）（殘基78～81、153～160 和244～248）、3 個胞外環(huán)（extracellular loops，ECLs）（殘基113～116、188～204和288～292）、1個螺旋H8（殘基322～335）以及C端（殘基336～338）。

這項工作采用了Ballesteros和Weinstein定義的GPCRs 殘基編號方案［16］。簡單地說，每個殘基的上標都被標記為X.YY，其中數(shù)字X 對應于殘基存在的TM 螺旋（從1 到7），YY 為相對于螺旋中最保守的殘基（YY 值為50）的位置。例如，如果TM2 中最保守的殘基是一個天冬氨酸（殘基號95），則其標識符為2.50，即Asp952.50。對于一個殘基Asn3107.45，表示在TM7 中最保守殘基Pro3157.50之前的第5個殘基處。

1.2 各向異性網(wǎng)絡模型（ANM）

在ANM 中，蛋白質被抽象為一個彈性網(wǎng)絡，其中氨基酸用一個節(jié)點來表示（通常選用Cα原子作為節(jié)點），當兩個殘基節(jié)點間的距離小于某一截斷半徑時（本文選取17 ?），它們之間用一個彈簧相連，所有彈簧的彈性系數(shù)相同。在ANM中，整個體系的勢能函數(shù)為：

其中，Rij和Rij0分別是節(jié)點i 和j 之間的瞬時和平衡距離，γ是彈簧的彈性系數(shù)。蛋白質的運動模式是由Hessian 矩陣H 所決定，矩陣的各元素為體系勢函數(shù)對位置的二階偏導，具體形式如下：

其中，每一個元素hij是3×3 的子矩陣。當i ≠j 時，hij可以表示為：

當i=j時，hii可以表示為：

對Hessian矩陣進行分解可獲得3N個本征值和本征向量。本征值小的本征向量對應蛋白質的慢運動模式，與功能有關；本征值大的對應蛋白質的快運動模式，與結構穩(wěn)定有關。

在ANM中，殘基的均方漲落以及殘基間漲落的交叉相關性為：

其中，λ-1k是H 矩陣的第k 個本征值的倒數(shù)，[uk]3i是H矩陣的第k個本征向量中的第3i個元素。

歸一化的交叉相關性為：

該值的取值范圍為-1～1。正負值分別表示殘基沿相同和相反的方向運動。絕對值越高，兩個殘基的相關性越強。零值意味著殘基的運動是完全不相關的。

1.3 微擾響應掃描模型（PRS）

線性響應理論（LRT）是一種線性擾動蛋白質殘基后，預測構象變化的方法［17］?；贚RT 的微擾響應掃描（PRS）方法被提出［15］，用于研究蛋白質的變構調控并識別關鍵殘基?；诤硕蒄= H ΔR，在受到外力擾動時，殘基位置的改變量為ΔR，其中H 是ANM 中3N × 3N 的Hessian 矩陣。當對網(wǎng)絡中的一個殘基i施加外力時：

觀測到其他殘基的響應：

其中，ΔRi是在外力對殘基i 擾動后，所有殘基位置的改變量。

對蛋白質的殘基i 施加一個力Fi，可以得到1個其他殘基產(chǎn)生的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，對其數(shù)字進行顏色編碼即為熱圖，熱圖中第ij個元素表示微擾節(jié)點i，節(jié)點j產(chǎn)生的平均響應。歸一化后的P矩陣中每一列的平均值代表了該列對應殘基的敏感性，即殘基響應其他微擾殘基的傾向性。所有殘基的敏感性就組成了敏感性特征曲線。該曲線峰值對應的殘基被認為是蛋白質變構中的功能性運動位點。相似的，每一行的平均值代表了該行對應殘基的效應性，即衡量微擾殘基影響其他殘基動力學的能力。相應曲線的峰值被認為對信號傳輸起重要作用的熱點殘基［18］。

基于上述方法的描述，本文使用MATLAB R2019a 對ANM 和PRS 的相關程序進行編寫并在DOP這個體系上進行研究。

2 結果與討論

2.1 基于ANM的實驗與理論B-factor的比較

溫度因子（B?factor）可以反映生物大分子的局部結構的柔性。本文構建了DOP 的ANM 模型，選用殘基的Cα原子作為節(jié)點。傳統(tǒng)的ANM只有一個截斷半徑參數(shù)，對截斷半徑的選取方法分為使用經(jīng)驗值［19］和在一定范圍內(nèi)通過最大化理論值與實驗值的相關性進行尋優(yōu)［12］。后一種方法能夠提高ANM 構建蛋白質模型的準確性，所以在這個工作中，采用尋優(yōu)的方法對截斷半徑進行優(yōu)化，其尋優(yōu)范圍為：6～24 ?，步長為1 ?。在截斷半徑為17 ?時，ANM 計算的理論B?factor 與實驗B?factor 的皮爾森相關系數(shù)（Pearson’s correlation coefficient，PCC）可以達到最大值（0.64），因此，本研究選擇17 ?。圖2 顯示了在最優(yōu)的參數(shù)下ANM 得到的理論和實驗B?factor 比較圖。從圖2 可以看出，ANM 可以較好地再現(xiàn)DOP 的柔性信息，TM 區(qū)域柔性最小，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區(qū)域，中間區(qū)域的局部最小值對應于ECL2 中的β 片層，表明ANM 對DOP 的結構有較好的識別能力。

有趣的是，發(fā)現(xiàn)其中幾個漲落較小的區(qū)域與DOP的功能相關（圖3a中圓圈）。第1個區(qū)域（殘基82～97）位于TM2，具有一定的保守性，其中Asp952.50是一個關鍵的鈉結合位點，實驗發(fā)現(xiàn)其突變可以破壞“鈉效應”。第2 個區(qū)域（殘基134～144）位于TM3，其中殘基Ser1353.39、Asn1313.35和Asp952.5形成了第1 個鈉離子結合位點。位于TM7中的第3 個區(qū)域（殘基301～320）也在鈉離子轉運中起著重要作用，其中殘基Asn3107.45和Asn3117.46形成了鈉離子的第2個配位殼，位于NP7.50xxY基序的殘基Asn3147.49和Tyr3187.53可以調節(jié)變構鈉離子結合，與GPCR的激活有關［21］。

2.3 基于ANM的DOP的快運動模式分析

快運動模式對應的是蛋白質局部結構的不規(guī)則性［22］。在快運動模式下，運動幅度較大的殘基通常是熱點殘基，對穩(wěn)定蛋白質結構至關重要［23］。圖3b 顯示了前兩個快運動模式的漲落分布圖。由圖3b 可以看出，11 個峰值區(qū)域對應的中心殘基分別 是 Asn671.50、 Ile882.43、 Leu912.46、 Asp952.50、Asp1283.32、 Asn1313.35、 Ser1353.39、 Val2666.40、Yyr3087.43、Ser3117.46、Asn3147.49。這些殘基分布在DOP跨膜螺旋的中心區(qū)域。

接下來，將根據(jù)現(xiàn)有的實驗和理論數(shù)據(jù)，對識別出的關鍵殘基功能進行討論。高度保守的殘基Asn671.50是DOP I?I 二聚體的一個結合熱點，對于穩(wěn)定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形成氫鍵，使復合物的結構更穩(wěn)定［26］。Ser3117.46和Asn3147.49是第2 個鈉離子配位口袋處殘基，在DOP 激活的過程中，Asn3147.49的位移會導致活性GPCR 的變構鈉結合口袋崩塌［26］。Ile882.43和Val2666.40分別位于TM2 和TM6，目前還沒有關于它們對DOP 變構作用的研究，值得被進一步探討。

綜上所述，識別的關鍵殘基在快運動模式下具有高度的活性，在穩(wěn)定DOP 結構方面發(fā)揮著重要作用或者在DOP與配體相互作用中扮演重要角色。

2.4 運動相關性分析

為了研究蛋白質運動過程中殘基之間的關聯(lián)性，通過式（7）計算了所有殘基之間的交叉相關系數(shù)，并用二維彩圖形式顯示出來（圖4a）。因為低頻運動模式常常與蛋白質功能有關，選擇對殘基漲落貢獻剛超過50%的低頻運動模式進行運動相關性的分析是一種通常的做法，可以提高信噪比［27］，所以本文選取對漲落貢獻剛大于50%的運動模式計算運動相關性。從圖4a 中可以看出，對角線區(qū)域被分成7部分，垂直于對角線的部分顯示相鄰TM 螺旋之間具有強的正相關性，表明ANM可以對DOP 的結構具有較好的識別能力。除此之外，胞外/胞內(nèi)跨膜螺旋與環(huán)狀區(qū)域之間也具有一定的正相互作用（圖中橙紅色區(qū)域），可能是因為DOP 的跨膜螺旋部分區(qū)域在空間上是相鄰的。特別是位于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之間的水介導氫鍵網(wǎng)絡與TM3 相互作用。這些結果表明，ICL3、TM3 和TM5 存在相互作用，進一步穩(wěn)定蛋白質結構［21］。Yyr1293.33、Met1323.36、Ile3047.39和Yyr3087.43形成一個結合口袋，它與DPI?287 形成了疏水相互作用。Lys2145.39、His2786.52和Yyr1293.33與KGCHM07形成一個極性相互作用網(wǎng)絡［26］?？傊珼OP 的跨膜螺旋之間的相互運動不僅可以穩(wěn)定構象狀態(tài)，還能促進與配體之間的結合。

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，已成功地用于揭示變構信號轉導過程中的關鍵區(qū)域［15］。本文對DOP中的殘基逐個微擾可以得到一個響應矩陣P（研究體系與方法1.3）。使用矩陣P 直接可以畫出DOP 的PRS響應熱圖（圖4b）。在PRS 圖譜中，第ij 個元素表示是擾動殘基i 對j 產(chǎn)生的影響，它衡量的是蛋白質結構的敏感性和效應性。

經(jīng)過分析，發(fā)現(xiàn)DOP 中具有大的敏感性殘基（變構信號的潛在接收器）主要分布在胞外環(huán)和胞內(nèi)環(huán)區(qū)域（圖4c），然而具有大的效應性殘基（傳播變構信號）主要分布在跨膜螺旋區(qū)域（圖4d）。圖4b的頂部是敏感性曲線圖。敏感性高的殘基有6簇（圖4c），中心殘基是Thr78ICL1、Glu1122.67、Asp158ICL2、 Asp193ECL2、 Arg244ICL3、 Asp290ECL3。這些殘基在前面的慢運動模式漲落分析中位于大的柔性區(qū)域，因此對于蛋白質的功能性運動具有一定作用。其中殘基Lys792.49在形成DOP 的I?I 二聚體時發(fā)揮作用［24］。在正變構調節(jié)劑BMS?986187 與DOP結合過程中，殘基Glu1122.67與BMS?986187形成疏水相互作用［28］。口袋殘基Trp114ECL1、Ile289ECL3和Arg291ECL3分別與配體KGCH07 形成ππ、疏水和陽離子?π 相互作用，可以進一步穩(wěn)定DOP的激活狀態(tài)［26］。由于ECL2是β鏈折疊，是所有阿片受體亞型的典型［29］，負責識別多種配體，因此位于此處的殘基如Asp193ECL2可能參與配體的結合。殘基Arg244ICL3通過與TM6和TM7的其他殘基形成廣泛的氫鍵網(wǎng)絡，在穩(wěn)定ICL3 方面起著關鍵作用，使DOP 穩(wěn)定在非活性狀態(tài)［21］。在研究鈉離子與DOP 結合的MD 模擬實驗中發(fā)現(xiàn)，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是激活狀態(tài)下的MOP 單體形成A?I 二聚體時的界面殘基［24］。MD 模擬研究發(fā)現(xiàn)，Cys1513.36對阿片受體與芬太尼的識別和相互作用至關重要［30］。Asp2105.35可以與激動劑KGCHM07形成水介導的鹽橋相互作用，穩(wěn)定DOP 的激活構象［26］。Leu246ICL3位于雙亮氨酸基序（Leu245ICL3～Leu246ICL3）中，實驗發(fā)現(xiàn)基序的缺失或Leu245 的突變會減緩DOP 的溶酶體靶向性［31］。殘基Arg2576.31與Leu2405.67、Arg244ICL3和Val243ICL3形成氫鍵網(wǎng)絡，與Asp2536.27形成鹽橋，它連接TM5 和TM6的胞內(nèi)結構，將DOP穩(wěn)定在非活性狀態(tài)［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.

進一步分析發(fā)現(xiàn)效應性高的殘基中也有部分殘基處在柔性較高的區(qū)域，這些區(qū)域與敏感性高的區(qū)域共同發(fā)揮功能性作用，如效應性強的殘基Ser421.25、Gly631.46和敏感性強的殘基Lys792.49共同耦合促進DOP二聚體的形成［24］。效應性強的殘基Asp952.50、Asn1313.35、Ser3117.46和敏感性強的殘基Asp293ECL3共同作用組成DOP 中鈉離子結合口袋［21］。

基于以上分析發(fā)現(xiàn)，通過PRS 確定的關鍵殘基在DOP 的二聚化、功能狀態(tài)穩(wěn)定以及與變構配體（如鈉離子和激動劑）的相互作用中起著重要作用，這些都參與了DOP的變構調節(jié)。

3 結 論

本工作探索了人DOP 的動力學特性，以及對蛋白質變構具有重要作用的關鍵殘基。首先，根據(jù)DOP的結構構建出最優(yōu)參數(shù)下的ANM，實驗與理論B?factor的PCC值達到0.64，并且識別出了DOP中的柔性區(qū)域，表明ANM可以較好的研究DOP的動力學特性。然后，從慢運動模式的漲落圖可以很好地識別DOP 的結構域，包括7 個TM 螺旋、3 個ICL 和3 個ECL。有趣的是，發(fā)現(xiàn)鈉離子結合口袋處具有最小的運動性，這有助于分析鈉離子的變構調節(jié)。此外，從快運動模式的漲落曲線來看，模式中的活性位點主要分布在TM1～3、TM6 和TM7中，這些區(qū)域在穩(wěn)定DOP 結構或與變構配體（如拮抗劑和鈉離子）的復合物結構方面起著重要作用。接著對蛋白質的運動相關性的分析可以看出胞外/胞內(nèi)跨膜螺旋與環(huán)狀區(qū)域之間存在較強的相互作用，這有助于蛋白質的構象穩(wěn)定以及與配體的結合。最后采用PRS 方法對DOP 進行微擾，識別出的關鍵殘基已經(jīng)在實驗或理論上發(fā)現(xiàn)直接參與DOP 二聚化和與配體結合，認為這些關鍵殘基對于DOP的變構通信至關重要。這項工作有助于對δ阿片受體變構調節(jié)的物理機制的理解，并為藥物設計提供有價值的信息。