基于增強采樣構建的隱式馬爾可夫狀態模型分析GLP-1R激動劑對GLP-1R激活機制

劉一卜，湯磊，范菊娣

1 貴州醫科大學藥學院藥物化學教研室，貴陽 550004；2 貴州省化學合成藥物研發利用工程技術研究中心

蛋白質的結構與功能之間的關系密不可分，觀察蛋白質構象變化有助了解其生理功能和作用機制。胰高血糖素樣肽1 受體（GLP-1R）是B1 類G 蛋白偶聯受體，具有細胞外結構域及由α-螺旋束構成的7 次跨膜核心域（跨膜結構域），GLP-1R 是Ⅱ型糖尿病的主要治療靶點［1］。PF-06882961 是激活GLP-1R 正構位點的激動劑，激活效力類似內源性多肽GLP-1，對包含cAMP、pERK1/2 等多個下游信號均有激活作用［1-2］。然而目前對于小分子激活GLP-1R的研究僅僅停留在濕實驗層面。在原子水平，PF06882961 激活GLP-1R 的具體機制仍未見報道。2020 年TIMPER 等［2］首次揭示了PF06882961 與GLP-1R 共晶結構，同時首次提出位于GLP-1R 跨膜結構域內由關鍵氨基酸之間的非鍵相互用形成的極性網絡對于GLP-1R 的激活具有重要意義，但小分子PF-06882961 與GLP-1R 結合后引起GLP-1R 發生怎樣的結構變化未見探討。本研究在增強采樣的基礎上構建了PF06882961與GLP-1R 復合物體系的隱式馬爾可夫狀態模型（HMM），基于此模型分析了PF06882961 激活GLP-1R 的機制，為靶向GLP-1R 的藥物設計提供結構上的參考。

1 材料與方法

1.1 研究體系 下載RCSB（www. pdbus. org）數據庫中全長GLP-1R 晶體結構6X1A（PDBID）［1］，Gal?axyFill（版本2.0）［3］進行氨基酸殘基的修補。使用CHARMM-GUI 構建PF06882961 與GLP-1R 復合物拓撲，補充Cys46-Cys71、Cys62-Cys104、Cys85-Cys126、Cys226-Cys296 之間的二硫鍵［3］，同時刪除該體系的水分子、納米體（鏈代碼D）等冗余分子，供后續構建常規動力學模擬（cMD）模型。

1.2 PF06882961激活GLP-1R的機制分析

1.2.1 PF-06882961 與GLP-1R 結合的cMD 體系構建 使用CHARMM-GUI（https：//www. CHARMMGUI. org/）中的 Charmm36 力場輔助構建PF06882961-GLP-1R 復合物拓撲，構建過程在軟件Gromacs（版本：2020.7）中運行［4］。復合物體系中插入100*100 ?2，厚度為74.84?的棕櫚酰油酰磷脂酰膽堿（POPC）雙層膜包裹跨膜結構域，溶劑化蛋白—配體—生物膜體系后，在蛋白質頂部和底端（雙側）插入22.5?的TIP3P水分子，加入Na+和Cl–維持電中性。采用線性約束求解器（LINCS）算法約束含氫原子的化學鍵鍵長，快速平滑粒子網絡Ewald（SPME）法計算靜電相互作用的長程部分，應用周期性邊界條件（PBC）消除溶劑盒子的邊緣效應。設定vdw 截斷值為1.2 nm，處理長、短程靜電相互作用的分界為1.2 nm。能量體系采用最陡下降法優化5 000步，當復合物體系最大的力小于1 000 kJ/（mol·nm）時，停止收斂。使用Nose-Hoover 擴展系綜進行溫度耦合（303.15 K），耦合兩次，每次按照蛋白質骨架、側鏈、水分子、脂質、離子的順序施加位置限制約束（kJ/（mol·?2））：10.0，5.0，2.5，2.5，20.0；5.0，2.5，2.5，2.5，0 ，步長為1 fs，單次平衡時長為25 ps。使用擴展系綜Parrinello-Rahman 構建壓力耦合，使用Semiisotropic 協議耦合膜體系，耦合4 次，每次按照蛋白質骨架、側鏈、水分子、脂質、離子的順序施加位置限制約束（kJ/（mol·?2））：2.5，1.0，1.0，1.0，0；1.0，0.5，0.5，0.5，0；0.5，0.1，0.1，0.1，0；0.1，0，0，0，0 ，運行10 ns 模擬，步長為1 fs，首次平衡時長為25 ps，其余3次為100 ps。后去除限制平衡100 ns，采用蛙跳式算法積分牛頓運動方程，時間步長為2 fs，最終構建了全長GLP-1R的動力學模型，包括細胞外結構域、跨膜結構域、細胞內結構域三個部分，修復四個二硫鍵以維持受體激活狀態。通過100 ns的常規動力學演化，復合物已達到相對穩定的結構，以此為基準構建高斯加速動力學（GaMD）體系［5］。

1.2.2 PF-06882961 與GLP-1R 結合的GaMD 體系構建 使用CHARMM-GUI 輔助構建PF06882961 與GLP-1R 結合的復合物拓撲，在軟件NAMD 中運行［5］。使用CHARMM36 力場構建復合物體系，在GLP-1R 的細胞外結構域和細胞內結構域插入厚度為22.5? 的TIP3P水分子層，加入Na+和Cl-以平衡電荷。采用朗之萬動力學控制溫度耦合，將復合物體系升溫至303.15 K。使用周期性邊界條件減小邊緣效應，在設定周期性邊界條件下，使用PME 計算庫侖相互作用（Coulombic potential），在恒溫條件下加壓，NAMD 將根據體系原子之間的力和動能來計算系統壓力。使用Langevin Piston 控制系統壓力，維持壓力在1 個大氣壓，按照默認振蕩時間常數和衰減時間常數，設置此時系統噪音溫度為環境溫度，此時不考慮氫鍵的增壓。在正式執行模擬之前，復合物體系將進行預模擬以確定升壓勢。最終進行50 ns GaMD模擬，并平行運行8次。

1.2.3 PF-06882961 結合 GLP-1R GaMD 軌跡的HMM 構建及驗證 將8 段平行運行的50 ns GaMD動力學軌跡作為分析對象，利用PyEMMA 工具包構建［6］HMM 模型，從一級結構［關鍵氨基酸殘基間的αC 間距（Glu247-His180；Glu364-Arg190）］、二級結構 ［關鍵α 螺旋間扭轉角（以關鍵氨基酸之間的夾角測算：Val365-Pro358-Ala350；Arg380-Phe390-Met397）］ 兩個層面對HMM 模型中的PF-06882961與GLP-1R 若干復合物構象進行聚類分析，獲得5 類復合物宏觀態構象（S1、2、3、4、5）。HMM 模型構建的大致流程見圖1。采用C-K 檢驗（Chapman-Kol?mogorov 檢驗）［7-8］驗證復合物宏觀態之間是否具有馬爾可夫性，結果顯示從一級結構及二級結構兩個層面將HMM 模型中的PF-06882961 與GLP-1R 若干復合物構象歸類為5 個宏觀態時具有可靠性（符合馬爾可夫性），可以用于觀察PF-06882961 與GLP-1R結合時的復合物構象。

圖1 PyEMMA工作流程示意圖

1.2.4 PF-06882961 與 GLP-1R 復合物宏觀態構象的差異分析 使用Pymol（open-source）對上述聚類分析所得的5 類復合物宏觀態構象（S1、2、3、4、5）進行可視化分析。使用RING（版本2.0，https：//ring.biocomputingup. it/submit）［9］將 GLP-1R 中跨膜結構域關鍵氨基酸之間的非鍵相互作用網絡進行可視化分析。

2 結果

S4、S5 為優勢態穩定存在，S1、S2、S3 為過渡態，不穩定，分布概率較小。由于S1、S2、S3 同屬于過渡態，且S1、S2、S3 在關鍵氨基酸和螺旋束夾角之間結構差異較小，故下文僅展示S4、S5 與S1 之間的構象差異。

2.1 從二級結構層面聚類的S1、S4 中GLP-1R 細胞外結構域及下游G 蛋白的構象變化情況 S1 與S4 在疊合后的構象差異見圖2A。從二級結構層面聚類的復合物宏觀態構象中，GLP-1R 細胞外結構域部分的α螺旋頂部發生傾斜，底部發生平移，使得小分子PF-06882961可以與Trp33發生更加直接的氫鍵相互作用和疏水相互作用。S4中PF-06882961 與GLP-1R 的結合模式見圖2B。相較于S1，S4 在跨膜螺旋束1 和跨膜螺旋束2（TM1、2）頂部發生內旋，使得PF-06882961 與Glu138、Lys197、Trp203 發生疏水相互作用、氫鍵以及p-π 共軛。同時TM7 頂部發生大幅度內旋，利于Phe381 與PF-06882961 之間發生π-π 堆積。G 蛋白的構象變化見圖2C。Gs 末端的阿爾法螺旋發生大幅度扭轉，Gi/Go 部分發生平移。

圖2 從二級結構層面聚類的復合物宏觀態構象變化圖

2.2 3 個復合物宏觀態中GLP-1R 跨膜結構域中關鍵氨基酸構成的極性網絡形成情況

2.2.1 從一級結構層面聚類的3 個復合物宏觀態中GLP-1R 跨膜結構域中關鍵氨基酸構成的極性網絡形成情況 S1（圖3A）中， GLP-1R 跨膜結構域中關鍵氨基酸Glu364-Arg190之間存在離子鎖，這使得二者間距較小，Tyr152-Tyr241-Glu364-Arg190 重排并形成非鍵相互作用組成的極性網絡，而Tyr402-Thr353-His180-Glu247 并未形成此類極性網絡。S4（圖3B）中，Glu364-Arg190 與His180-Glu247 之間均存在離子鎖，氨基酸距離的減小使得Tyr402-His180-Thr353-Glu247、Arg190-Glu364-Tyr241 重排構成了極性網絡，氨基酸殘基之間通過大量的氫鍵網絡相互作用，體系較為穩定。而同時，Glu364-Tyr241-His180-Glu247 之間通過離子鎖以及氫鍵相互作用形成了串聯，構建了穩定存在的新的極性網絡。S5（圖3C）中，關鍵氨基酸Trp33［1］與其他氨基酸通過多種非鍵相互作用穩定，Phe385-Tyr203-Tyr148之間穩定存在一個π-π堆疊網絡，Glu247-His180之間存在的離子鎖并未使此處的極性網絡重排，Glu364-Arg190 之間既沒有構成氫鍵相互作用也無離子鎖，Glu364-Tyr241 之間構成了氫鍵相互作用。Glu364-Tyr241-His180-Glu247之間存在氫鍵相互作用。

圖3 從一級結構層面聚類的3個復合物宏觀態中GLP-1R跨膜結構域中關鍵氨基酸構成的極性網絡圖

2.2.2 從二級結構層面聚類的3 個復合物宏觀態中GLP-1R 跨膜結構域中關鍵氨基酸構成的極性網絡形成情況 S1 中GLP-1R 跨膜結構域中氨基酸殘基間相互作用網絡見圖4A。Trp33 與π-π 堆疊網絡（Phe385-Trp203-Tyr148）之間通過氫鍵相互作用結合。Glu364-Tyr241-Arg190 之間的極性網絡通過氫鍵相互作用與His180-Glu247-Tyr402 之間發生串聯，而Glu247-His180 之間存在離子鎖。S4中GLP-1R 跨膜結構域中氨基酸殘基間相互作用網絡見圖4B。S4 中，Phe385-Trp203-Tyr148 之間存在π-π 堆疊，而His180-Glu247 之間的離子鎖打開，以氫鍵相互作用相連。Glu364-Tyr241 之間存在氫鍵相互作用，并與His180-Glu247 形成了一個新的重排網絡。S5 中GLP-1R 跨膜結構域中氨基酸殘基間相互作用網絡見圖4C。S5 中Trp203-Tyr148 之間存在π-π 堆疊，兩個離子鎖均呈現出打開的狀態；Tyr402-Glu247-Arg180 重排并形成極性網絡，通過氫鍵相互作用與Glu364-Tyr241 之間發生串聯，形成新重排網絡Glu364-Tyr241-His180-Glu247。

圖4 從二級結構層面聚類的3個復合物宏觀態中GLP-1R跨膜結構域中關鍵氨基酸構成的極性網絡圖

3 討論

HMM 模型和全原子分子動力學模擬提供了觀察生物學過程的原子結構的視野［10-12］。全原子分子動力學可以還原生物膜、體液環境等復雜的環境，而增強采樣動力學則在分子動力學的基礎上，人為縮短了構象從局部最小值跨越到下一個局部最小值之間所需的時間（更快地翻越能壘），從而在有限時間內更高效地探索研究體系的自由能景觀［13］。HMM 模型可以更有效地解讀分子動力學過程獲得的結果，并根據研究目的、不同的聚類特征處理研究體系的動力學軌跡。對于PF-06882961 激活GLP-1R 的研究具有重要意義。PF-06882961作為目前惟一的GLP-1R 完全激動劑，對于GLP-1R 下游多條通路都具有不同程度的激活效力，與其他的部分激動劑之間的激活機制之間存在的差異性有重要的研究價值［1，13］。

本研究構建了包含磷脂雙層膜和下游信號蛋白（G 蛋白）的GLP-1R 的全原子模型，并基于此實現了100 ns 經典動力學模擬，同時構建了GaMD 體系，實現不受集體變量約束的50 ns 增強采樣全原子分子動力學模擬。基于8 條平行GaMD 軌跡，構建了PF-06882961 與GLP-1R 復合物體系的HMM 模型，并采用C-K檢驗驗證了宏觀態之間是否具有馬爾可夫性［10］。在驗證HMM 模型的科學性后，繪制了在兩種聚類特征下［二級結構層面 （α 螺旋間扭轉角）和一級結構層面］復合物構象之間的時間概率轉移矩陣，有助于幫助捕捉PF-06882961 在50 ns 下誘導GLP-1R發生重要轉變的關鍵PF-06882961與GLP-1R復合物構象及相應的分布概率，分析PF06882961結合GLP-1R 后，GLP-1R 的細胞外結構域部分以及下游蛋白（G蛋白）的構象轉變情況。

文獻［14-15］報道，小分子PF-06882961 作為GLP-1R 激動劑，在激活GLP-1R 的過程中形成了兩對關鍵的極性網絡（Tyr152-Tyr241-Glu364-Arg190；Tyr402-Thr353-His180-Glu247），并猜測由于這兩對極性網絡的出現使得GLP-1R 穩定在激活狀態，從而引導下游G 蛋白發生信號傳導。但通過對PF-06882961 與GLP-1R 結合的前期研究中，我們觀察到小分子PF-06882961 在100 ns 常規動力學期間，其中一對關鍵氨基酸間的離子鎖（Arg190-Glu364）并不穩定，同時獲得了存在“打開”和“閉合”兩種構象，由于模擬時間尚短，無法對這樣的現象進行定性。故而本研究選擇采用增強采樣動力學進行更深入的模擬，同時構建HMM 模型以更科學地采樣動力學模擬軌跡。

GLP-1R跨膜區極性網絡的形成對GLP-1R發生構象轉變具有重要影響［1］。本研究從兩個層面分別觀察了PF-06882961 結合GLP-1R 后，GLP-1R 的構象轉變：①以跨膜結構域中形成極性網絡的關鍵氨基酸之間的間距（Glu247-His180；Glu364-Arg190）作為觀察指標，觀察在PF-06882961 結合GLP-1R 的動力學過程中極性網絡的變化情況；以GLP-1R 的細胞外結構域中的關鍵α 螺旋的扭轉角度（Val365-Pro358-Ala350；Arg380-Phe390-Met397）作為觀察指標分析HMM 模型，以觀察在動力學過程中GLP-1R細胞外結構域中部分構象的變化情況。

從兩個層面進行聚類分析時，PF-06882961 與GLP-1R 復合物體系動力學軌跡的采樣具有差異。從二級結構層面進行聚類時，在50 ns GaMD 過程中對PF-06882961 與GLP-1R 復合物體系的采樣更加連貫；從一級結構層面進行聚類時，PF-06882961 與GLP-1R 復合物體系的采樣并不連貫。本研究結果顯示，相較于從一級結構層面獲得的復合物宏觀態構象，從二級結構層面獲得的復合物宏觀態優勢構象更加穩定（由過渡態轉變至優勢構象時轉變概率更高）。這表明，以更加宏觀的（二級結構或三級結構）角度聚類所得的復合物構象更加能夠反映研究體系的真實生物過程。

本研究中， S4、S5具有較低的自由能，二者結構差異較小（RMSD=0.802），在這兩個復合物宏觀態中，PF-06882961 中的環氧丁基與GLP-1R 的Gln22形成了氫鍵相互作用。S1、S2、S3為過渡態，這3個復合物構象中PF-06882961的差異較小（RMSD分別為0.940、0.715），說明此時PF-06882961 在與GLP-1R結合時具有相似的結合模式，這種結合模式下，由于PF-06882961 中環氧丁基構象的改變形成空間位阻，GLP-1R 的Gln22 與PF-06882961 的環氧丁基之間的氫鍵消失，表現為結合力較弱的范德華力相互作用，導致整個體系的能量變高。

從二級結構層面聚類所得宏觀態中，優勢構象態S4 中可以觀察到GLP-1R 經過PF-06882961 激活后，下游的G 蛋白（Gs部分）發生了偏轉，Gi/Go部分發生平移。這或許是PF-06882961 能夠成功激活GLP-1R 的有力證據，為后續開發PF-06882961 類似物提供下游構象轉變方面的參考。同時在過渡態S1向S4轉變的過程中的細胞外結構域部分閉合，可以視為PF-06882961 激活GLP-1R 過程中，發揮激活效力的一個重要二級結構層面的信號。

從兩個層面進行聚類分析時，通過觀察GLP-1R跨膜結構域的關鍵氨基酸極性網絡形成情況，我們得到了一些有關PF-06882961 激活GLP-1R 的重要結構性啟示。一個關鍵的π-π 堆疊網絡（Phe385-Tyr203-Tyr148）穩定存在于由PF-06882961 激活GLP-1R的激活狀態中，同時通過氫鍵網絡與關鍵氨基酸Trp33 結合，對GLP-1R 的細胞外結構域具有穩定作用。觀察從兩個層面聚類所得優勢宏觀態，其中兩個極性網絡：Tyr152-Tyr241-Glu364-Arg190；Tyr402-Thr353-His180-Glu247 并不穩定。從一級結構層面進行聚類時，兩個優勢宏觀態（S4、S5）中GLP-1R的跨膜區氨基酸只能重排出部分極性網絡，且關鍵氨基酸中存在離子鎖。從二級結構層面聚類所得的優勢宏觀態（S4、S5）中，我們發現，構成GLP-1R 細胞外結構域和下游G 蛋白發生構象改變的優勢宏觀態，形成了一個新的重排網絡（S5）：Glu364-Tyr241-His180-Glu247。這個新的極性網絡相比文獻報道的兩個網絡更穩定地存在于聚類所得的優勢宏觀態中，在PF-06882961 激活GLP-1R 過程中的穩定存在，穩定了GLP-1R的跨膜結構域。

總之，PF-06882961 激活GLP-1R 具有2 個關鍵的結構改變：一個由Phe385-Tyr203-Tyr148 組成的π-π 堆疊網絡穩定GLP-1R 的細胞外結構域；一個由Glu364-Tyr241-His180-Glu247 重排而成的新的極性網絡。在其他先導化合物激活GLP-1R 過程中，出現這兩個關鍵網絡將大概率預示著該先導化合物具有較高的激活潛力，從而加速激活GLP-1R 的小分子藥物的研發。GLP-1R 隸屬于G 蛋白偶聯受體B家族，后者在人體內具有廣泛的分布與相當類似的三級結構［16-17］，有關于GLP-1R 激活機制的研究將一定程度上給予其他G蛋白偶聯受體機制研究一定的啟發和參考。