GIRK2通道外流门控机制的分子动力学模拟

王 瑞，赵伊婷，石定媛，李岱霖，林 嫣

（厦门理工学院环境科学与工程学院，福建 厦门 361024）

G 蛋白门控内向整流钾离子（G Protein-Gated Inwardly Rectifying K+，GIRK）通道通过细胞兴奋性调节在神经系统或组织器官的G蛋白信号通路调控方面发挥重要作用，在神经递质释放、成瘾药物研制表现出极具前景的靶向治疗潜力［1-3］。细胞兴奋性调节通过与多种耦联受体结合介导胞内外K+传递实现［4-5］。该传导过程受限于通道状态，当通道受到激活而使门控呈现打开状态时，允许K+实现选择性穿越迁移。因此，探究门控调节机制对于深入理解相关病理、药效机制，开发特异性靶向治疗调节剂具有重要意义。

作为GIRK 家族成员，GIRK2通道的X 射线晶体结构于2011年被首个解析，得到结构生物学研究广泛关注。GIRK2 结构组成包括跨膜区（Transmembrane domain，TM）和细胞溶质区（Cytosolic domain，CTD）。K+的外向迁移须穿越三个限制性结构域：CTD 顶端的G-loop 门控、螺旋束交叉门控（Helix Bundle Crossing，HBC）、高度保守的选择性过滤器（Selectivity Filter，SF）［6］。门控激活依赖于细胞膜上磷脂酰肌醇（Phosphatidylinositol（4，5）Bisphosphate，PIP2）的稳定作用和激活因子调控［7-8］。目前，GIRK2 通道已发现的激活因子包括Na+［9］、G 蛋白βγ 亚单位（Gβγ）［10］、乙醇［11］等。考虑到内向整流特性，关于GIRK2 通道的门控机制研究均以K+内流为前提。通过晶体结构叠合提出了CTD 结构域相对TM 的逆向旋转机制可激活GIRK2 通道［12］。此外，分子动力学研究结果表明，GIRK2 通道通过CTD 摇摆、逆向旋转、TM2 倾斜运动被激活，其中HBC 门控激活是CTD 的摇摆和TM2的倾斜导致，而G-loop门控由于CTD结构域的逆向旋转和向内摆动激活［6，13-14］。尽管GIRK2通道的内流机制得到深入研究，但与之相对的外流门控机制未见报道。已有研究表明，利用离子通道外向电流阻滞或激活已成功研制出众多药物分子，如dofetilide［15］、cromakalim、nicorandil［16］等，可见外流机制研究对于探讨药效机制是必要的。本文借助分子力学、分子动力学理论，以GIRK2 通道为对象，在双层脂膜环境下探究不同因子调控下的构效关系，揭示GIRK2的外流门控机制，为GIRK 家族的通道激活调节机制研究提供结构生物学信息，以助于深入理解外流调控型病理或药效产生机制。

1 GIRK2体系分子动力学模型搭建

为了探究饱和浓度配比下GIRK2 通道的外流门控调节机制，选取GIRK2 通道（蛋白质数据库代码：4KFM［17］）晶体结构为初始结构模板，如图1所示。

图1 GIRK2通道晶体结构图[17]Fig.1 Crystal structure of GIRK2 channel[17]

对于未饱和浓度配比下，考虑GIRK2 通道空间位阻和双层脂膜的锚定效应等因素，使得PIP2、Gβγ 与GIRK2 的结合位点固定、单一，致使其实验条件众多，计算和分析量巨大，本文仅讨论外向电场作用下GIRK2 的门控机制［17］。通过化学修饰增删内源性调节因子PIP2、Na+、Gβγ 搭建3 个初始模型：GIRK2-PIP2-Na+（GIRK2-Na+）、GIRK2-PIP2-Gβγ（GIRK2-Gβγ）、GIRK2-PIP2-Gβγ -Na+（GIRK2-Na+-Gβγ），GIRK2 的模型搭建与计算具体如图2 所示。使用H++服务器（http：//biophysics.cs.vt.edu/）添加残基支链缺失的氢原子［18］，在生理条件（pH=7.0）下通过计算各残基的pKa 确定残基质子化状态。将GIRK2复合物浸没在1-棕榈酰-2-油酰-Sn-甘油-3-磷酸胆碱、1-棕榈酰-2-油酰-Sn-甘油-3-磷脂酰乙醇胺、1-棕榈酰-2-油酰-Sn-甘油-3-磷脂酰丝氨酸和胆固醇构成的显式双层脂膜（分子比为25∶5∶5∶1）［19］。添加显示水溶剂（GIRK2-Na+、GIRK2-Gβγ、GIRK2-Na+-Gβγ 溶剂盒子尺寸分别为116 Å×116 Å×150 Å、162 Å× 162 Å× 170 Å、162 Å× 162 Å× 170 Å）和150 mmol·L-1KCl溶液的生理环境。据此搭建的3 个模拟体系，调节因子PIP2、Gβγ、Na+的初始浓度约2～ 4 mmol·L-1，因子与GIRK2 的配位比4∶1。针对复合物模型中的蛋白质、脂膜、调节因子，分别采用FF14SB、Lipid17 和GAFF2力场［20］描述原子受力和运动情况。据此搭建的模拟体系包含160 000到400 000个原子。

图2 GIRK2的模型搭建与计算Fig.2 Model building and calculation of GIRK2

本文以外向电场下分子动力学平衡轨迹为对象，进行钾离子运动轨迹、通道距离、结构域几何参数分析。运动轨迹主成分分析和几何参数计算分别借助AmberTools17 和MDAnalysis 0.18 软件［21］实现。通道结构绘制借助VMD1.9.3软件包完成［22］。所有计算在配备英伟达GTX 1080图形卡的服务器中进行。计算环境为Community Enterprise Operating System（CentOS）7 和NVIDIA Compute Unified Device Architecture（CUDA）8.0。

2 微秒尺度全原子分子动力学模拟

对搭建好的3个初始模型分别进行几何优化及微秒尺度全原子分子动力学模拟（将体系每个原子当做一个小球，求解体系牛顿力学方程）计算。首先，对模型体系进行能量最小化计算；其次，对体系依次进行升温、平衡动力学计算；最后，在外向电场下进行平衡动力学等。能量最小化采取最速下降算法（10 000步）和共轭梯度算法（10 000步）。体系升温采用Langevin 恒温算法，时间步长0.5 fs以避免内部扰动。平衡模拟首先在无电场条件下，采取恒温恒压系综进行0.2 μs，而后在外向电场（跨膜电压约200 mV）［23-24］及恒温恒容系综进行0.5 μs。该电压在GIRK2 二级结构保持完好情况下用以促进K+外流，更高的外电场会导致结构不稳定。长程静电作用采取Particle-Mesh Ewald 方法，设定截断值为10 Å。借助溶质的氢质量重分配算法［25］，以4 fs 为时间步长加速分子动力学模拟，并对溶剂采用SHAKE算法。模拟工作借助AMBER16中的PMEMD.CUDA程序完成。

3 实验结果与讨论

3.1 Na+和Gβγ亚基协同作用激活GIRK2通道

为了探究外流传导门控机制，分别对GIRK2-Na+、GIRK2-Gβγ和GIRK2-Na+-Gβγ体系在外向电场下进行分子动力学模拟。均方根偏差结果如图3显示，所有模拟体系在200 ns左右达到平衡态，适于结果与讨论分析。

图3 GIRK2通道3个体系的Cα原子均方根偏差Fig.3 Root mean square deviation of Cα atom in three systems of GIRK2 channel

外流穿越各门控的钾离子数和消耗时间如图4所示。

图4 外流穿越各门控的钾离子数和消耗时间Fig.4 Number of outflowing potassium ions and consumption time through each gate

结果显示调节因子Na+和Gβγ共同作用可激活通道，但是各自调控GIRK2通道时，并未检测到K+流。从图4（c）可见，在GIRK2-Na+-Gβγ 体系，有6 个K+实现外流穿越G-loop、HBC 门控，平均时间分别为77、1 ns，其中1 个K+进一步外向穿越了SF，消耗时间262 ns。从图4（a）可见，在GIRK2-Na+体系，仅有2 个K+通过G-loop 门，平均通过时间为180 ns，未检测到HBC 和SF 穿越事件。从图4（b）可见，在GIRK2-Gβγ体系，穿越G-loop门控的K+为1个，通过时间为134 ns，外流通过HBC门消耗时间为1 ns，但未检测到SF穿越事件。简言之，GIRK2-Na+-Gβγ体系，K+容易穿越三个门控迁移至胞外，即Na+和Gβγ 亚基协同作用可激活GIRK2 通道。与之不同的是，在GIRK2-Na+体系，K+仅穿过了G-loop 门；在GIRK2-Gβγ 体系，穿越G-loop 门的K+减少，通过HBC 门的增多，即两个体系未出现通道完全穿越事件。造成该现象的原因可能是由于通道内向整流特性（K+在既定驱动电压下由细胞外向细胞内的内向迁移强于反向驱动的外向渗透），使K+外向迁移相较内向困难得多，所需时长增加。已有实验结果有力地支持了这一观点。以GIRK2-Gβγ 体系为例，K+内向穿越GIRK2 通道所需时间约为200 ns［6］，而外向穿越需1 435 ns，表现出显著的内向整流特性。此外，本文的外向电流顺序（GIRK2-Na+-Gβγ ＞ GIRK2-Gβγ ＞ GIRK2-Na+）表明，在调节因子的激活能力上，Na+-Gβγ ＞ Gβγ ＞Na+。该结论与电生理实验结果表现出良好一致性，显示出本文计算结果的可靠性［6，14］。

3.2 门控因子Na+主要调节GIRK2的G-loop门控

已有报道显示Na+通过与GIRK2 或GIRK4 的CD 环结合稳定其与G-loop 门控的相互作用［26］。本文分子动力学模拟结果表明Na+通过稳定G-loop 门控激活GIRK2，支持上述结论。3 个体系的G-loop 和HBC门控最小距离分布，如图5所示。

图5 GIRK2通道G-loop和HBC门控最小距离分布直方图Fig.5 Histogram of minimum distance distribution of G-loop and HBC gates

向GIRK2 加入Na+，G-loop 门控最小距离分布相较于Gβγ 体系明显右移（图5（a）），且超过K+穿越所需最小距离5.69 Å［6］，而HBC最小距离分布未受到显著影响（图5（b））。表1是GIRK2通道的G-loop、HBC 门控打开概率，其中GIRK2-Na+体系的G-loop 门控打开概率高达94.6%，即长时间处于激活状态，而HBC门控打开概率仅为61.0%。这解释了GIRK2-Na+体系中外流K+容易穿越G-loop门控而很难通过HBC的原因。综上，Na+主要调节GIRK2的G-loop门控。

表1 GIRK2通道的G-loop和HBC门控打开概率Table 1 Open probability of G-loop and HBC gates of GIRK2 channel 单位：%

3.3 调节因子Gβγ亚基主要调节HBC门控

当Gβγ加入GIRK2体系时，HBC门控最小距离分布发生明显右移（图5（b）），而G-loop门控最小距离分布未发生显著移动（图5（a））。表1 显示该体系HBC 门控打开概率（73.6%）高于G-loop门控（15.8%），表明HBC更易激活。该结果解释了GIRK2-Gβγ体系外流K+一旦穿越G-loop门控，即容易迅速通过下一门控HBC。相较于GIRK2-Na+体系的G-loop门控打开概率94.6%，GIRK2-Gβγ体系则低得多。该结果与外流K+数目相符，即在同样的模拟时间（500 ns），有2个K+外流穿越GIRK2-Na+体系的G-loop门控，而只有1个穿越GIRK2-Gβγ体系的G-loop门控。综上，Gβγ主要调节GIRK2通道的HBC门控。

3.4 G-loop和HBC门控机制

图6显示了不同因子作用下GIRK2通道的G-loop和HBC门控激活机制。

图6 不同因子作用下GIRK2通道的G-loop和HBC门控激活机制Fig.6 G-loop and HBC gated activation mechanisms of GIRK2 channel under different factors

从图6 可见，当通道与Na+作用时，主成分分析结果显示存在CTD 结构域摇摆和跨膜螺旋2（TM2）相对CTD的逆向旋转（胞外到胞内视角）运动，本文利用TM2-CTD二面角和TM2-CTD旋转角描述上述两种运动，具体如图6（a）～6（c）所示。图6（d）表明随着TM2-CTD 二面角减小（CTD 平均内向摆动148.8°）和逆向旋转角增加（平均逆向旋转11.9°），G-loop 门控尺寸相应增加，这与GIRK2-Na+体系G-loop 门控拥有最大尺寸（图5（a））和打开概率（表1）结果一致。由此，CTD 结构域内向摆动和TM2-CTD逆向旋转可导致G-loop门激活。

当通道与Gβγ 作用时，主成分分析结果表现出CTD 摆动和TM2 倾斜两种运动模式［12］。随着TM2-CTD 二面角增加，导致TM2 外向倾斜约38.3°，HBC 的Cα 距离相应增加（GIRK2-Na+-Gβγ），即CTD结构域的外向摆动和TM2的外向倾斜可导致HBC激活（图6（e））。

4 结论

本文借助分子力学、分子动力学理论，采用全原子分子动力学模拟方法，在外向电场和双层脂膜环境下，探究GIRK2-Na+、GIRK2-Gβγ、GIRK2-Na+-Gβγ三个体系的外流门控机制。结果表明，调节因子Na+和Gβγ 亚基协同作用可激活GIRK2 通道，其中Na+使得GIRK2 通道G-loop 门控打开概率高达94.6%，显现G-loop门控调控作用，而Gβγ主要影响HBC门控，使其打开概率高达73.6%；CTD结构域内向摆动约148.8°和TM2-CTD逆向旋转约11.9°可导致G-loop激活，而CTD外向摆动和TM2外向倾斜约38.3°可导致HBC 激活。考虑到GIRK 家族成员关键残基保守性和序列相似性，GIRK2 通道的门控调节机制可为家族的通道激活机制探究提供理论依据和参考，为深入理解外流调控型病理或药效产生机制提供结构生物学信息。