基于增强采样构建的隐式马尔可夫状态模型分析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蛋白偶联受体机制研究一定的启发和参考。