基于弹性和复杂网络模型的内向整流型钾通道Kir3.2的动力学及关键位点研究

Study on the Kinetics and Key Sites of Inward Rectifier Potassium Channel Kir3.2 Based on the Elastic Network and Complex Network Models

摘要: 内向整流型钾通道蛋白Kir3.2能够调节钾离子的跨膜流动，具有重要的生理功能。由于其与神经和心肌系统的疾病相关，因此是一个重要的药物靶点。为研究该通道的动力学和功能性关键位点，我们首先利用高斯网络模型(GNM)对该通道进行建模，分析其结构的柔性，及各功能区域间的运动相关性；然后，构建运动相关性加权的复杂网络模型，通过依次敲除节点获得网络特征路径长度的变化来识别关键位点。结果发现，GNM重现了通道各区域的柔性；界面螺旋、跨膜螺旋和CTD区域间存在显著的运动正相关，这促进了通道与配体的结合和结构的稳定；所识别的关键残基对通道的激活、结构的稳定和变构信号交流起重要作用。本工作有助于加强对Kir3.2通道工作机制的理解，所识别的功能位点可为药物设计提供有价值的信息。

Abstract: The inward rectifier potassium channel Kir3.2 can regulate the flow of potassium ions across the cell membrane, playing important physiological functions. Due to its association with diseases of the nervous and myocardial systems, it is an important target for drug design. To study its kinetics and functional key sites, we first construct its GNM model, analyze the structural flexibility and the mo-tion correlations between different functional regions. Then, the motion correlation-weighted com-plex network model is constructed, and the key sites are identified in terms of the changes of char-acteristic path lengths upon the nodes knocked one by one. The results found that the GNM model reproduces the flexibilities of different functional regions; evident positive correlations exist among the interfacial helices, transmembrane helices and CTD regions, which promotes channel-ligand binding and structural stability; the identified key residues play an important role in channel acti-vation, structural stability and signal communication. This work helps to enhance the understand-ing of Kir3.2 channel’s working mechanism, and the identified key sites can provide valuable in-formation for drug design.

1. 引言

内向整流型钾通道Kir3.2能够调节钾离子的跨膜流动，是内向型整流钾离子通道家族的一员，在许多生理过程中发挥重要作用 [1] 。Kir3.2通道的突变可引起许多疾病，包括原发性醛固酮增多症、安德森综合征、巴特综合征和先天性高胰岛素血症 [2] [3] 。实验发现，该通道主要由4个相同的亚基组成，它们构成了两个重要的功能区域：跨膜结构域(transmembrane domain, TMD)和胞质结构域(cytoplasmic domain, CTD) [4] 。目前认为，离子跨膜流动可能是由两个串联的门来调节的：由TMD的内部螺旋形成的内螺旋门 [5] ，由CTD顶端的G环形成的G环门 [6] ，阴离子脂质磷脂酰肌醇4，5-二磷酸(phosphatidylinositol 4, 5-diphosphate, PIP2)是该通道的激活剂 [2] [7] [8] 。目前，尽管有一些工作在探讨该离子通道的门控机制，并有所发现 [9] [10] ，但对其内在动力学如何调控通道的打开与关闭，哪些残基在门控中发挥重要作用还不完全清楚 [11] 。

分子动力学与复杂网络模型方法结合已被广泛用于蛋白质变构信号交流、蛋白质折叠、关键残基识别和蛋白质相互作用的研究 [12] [13] 。Sethi等人将全原子分子动力学(Molecular dynamics, MD)模拟与复杂网络模型结合，构建了残基间运动相关性加权的动态复杂网络模型 [14] 。由于MD模拟非常耗时，随后人们又发展了粗粒化的弹性网络模型 (Elastic network model, ENM)方法，用该方法获得的残基运动相关性进行加权构建的动态复杂网络模型已被用于一些大的分子机器的变构动力学研究。我们组曾利用该策略研究了U1A-snRNA结合偶联变构的动力学过程 [15] 。

在本工作中，我们基于Kir3.2通道的三维结构构建了GNM模型，预测了通道结构的柔性，分析了残基间的运动偶联。之后，构建了一个动态复杂网络模型来识别Kir3.2中参与变构和信号交流的关键残基。

Figure 1. Cartoon diagram of Kir3.2 channel protein structure with the important functional regions labelled

图1. Kir3.2通道蛋白结构的卡通图，标注了Kir3.2的重要功能区

2. 方法和体系

2.1. 研究体系

内向整流型钾通道Kir3.2的高分辨率结构来自蛋白质数据库(Protein Data Bank, PDB) [16] ，其PDB ID为3SYA，该体系是PIP2结合态结构，如图1所示。Kir3.2通道主要由4个相同的亚基组成，共有1292个残基，每个亚基包含323个氨基酸。

2.2. 高斯网络模型(GNM)

高斯网络模型(Gaussian Network Model, GNM)是粗粒化弹性网络模型中的一种，它将蛋白质的三维结构模拟成一个弹性的网络。其中，氨基酸残基的Cα原子被简化为网络中的节点，节点间距离小于某一截断半径时，认为它们之间存在相互作用，并用弹性系数统一的弹簧进行连接。基于上述简化，弹性网络模型总的势能为 [17] ：

$V_{\text{GNM}}=\frac{1}{2}\gamma \left[\Delta R^{\text{T}}\left(\Gamma \otimes E\right)\Delta R\right]$ (1)

其中，γ为弹簧的弹性系数，列向量ΔR代表蛋白质中Cα原子的位置涨落，N为氨基酸个数，上标T表示向量的转置，Γ表示N × N的基尔霍夫对称矩阵，其矩阵元定义为：

${\Gamma }_{ij}=\{\begin{array}{ll}-1\hfill & \text{if}i\ne j\text{and}{R}_{ij}\le r_c\hfill \\ 0\hfill & \text{if}i\ne j\text{and}{R}_{ij}>r_c\hfill \\ -{\displaystyle \underset{i,j\ne i}{\sum }{\Gamma }_{ij}}\hfill & \text{if}i=j\hfill \end{array}$ (2)

其中，R_ij是蛋白质中残基i和j间的距离，r_c为截断半径。

残基i的均方涨落(Mean Square Fluctuation, MSF)和残基i与j间涨落的交叉相关性分别为：

$\langle \Delta R_i\cdot \Delta R_i\rangle =\frac{3k_{B}T}{\gamma }{\left[{\Gamma }^{-1}\right]}_{ii}=\frac{3k_{B}T}{\gamma }{\displaystyle \underset{k=2}{\overset{N}{\sum }}{\lambda }_k^{-1}{\left[u_k\right]}_i}{\left[u_k\right]}_i$ (3)

$\langle \Delta R_i\cdot \Delta R_j\rangle =\frac{3k_{B}T}{\gamma }{\left[{\Gamma }^{-1}\right]}_{ij}=\frac{3k_{B}T}{\gamma }{\displaystyle \underset{k=2}{\overset{N}{\sum }}{\lambda }_k^{-1}{\left[u_k\right]}_i{\left[u_k\right]}_j}$ (4)

其中，k_B是波尔兹曼常数，T为绝对温度。Γ矩阵的逆可分解为：

${\Gamma }^{-1}=U{\Lambda }^{-1}{U}^{\text{T}}={\displaystyle \underset{k=2}{\overset{N}{\sum }}{\lambda }_k^{-1}{u}_k{u}_k^{\text{T}}}$ (5)

其中，U为正交矩阵，每一列u_k (1 <k ≤ N)为Γ的本征向量，表示蛋白质的运动模式；Λ为对角矩阵，其矩阵元为Γ的特征值λ_k。根据Debye-Waller理论，B因子可以通过下式计算：

$B_i=\frac{8{\pi }^2}{3}\langle \Delta R_i\cdot \Delta R_i\rangle =\frac{8{\pi }^2{k}_{B}T}{\gamma }{\displaystyle \underset{k=2}{\overset{N}{\sum }}{\lambda }_k^{-1}}{\left[u_k\right]}_i{\left[u_k\right]}_i$ (6)

残基i和j之间归一化的运动交叉相关性为：

$C_{ij}=\frac{\langle \Delta R_i\cdot \Delta R_j\rangle }{{\left[\langle \Delta R_i^2\rangle \times \langle \Delta R_j^2\rangle \right]}^{\frac{1}{2}}}$ (7)

交叉相关性的取值范围为−1到1，正值表示残基的运动方向相同，负值表示运动方向相反。该值的绝对值越大表示两个残基间的运动相关性越强。

2.3. 氨基酸复杂网络模型

氨基酸复杂网络模型 [18] 已被用于研究蛋白质变构交流 。在复杂网络模型中，常用残基Cα原子代表网络中的节点，用一定的截断半径r_c (这里为7.0 Å)来定义网络中节点间的连接 [15] ，当对边进行加权时，复杂网络的邻接矩阵元素为：

$A_{ij}=\{\begin{array}{l}{w}_{ij}i\ne j\text{and}{R}_{ij}<r_c\\ 0i=j\text{and}{R}_{ij}\ge r_c\end{array}$ (8)

这里，r_ij是残基i和j之间的距离，权重来自残基间的运动相关性：

$w_{ij}=-\log \left(\left|C_{ij}\right|\right)$ (9)

网络特征路径长度(Characteristic Path Length, CPL)是网络中所有节点间路径长度的平均值，CPL反映了网络内部节点间的连通性和信息传递效率。k节点对网络通信的贡献可用该节点敲除后网络CPL的变化来度量，通常用Z-score来衡量：

${\text{Z-score}}_k=\frac{\Delta CP{L}_k-\langle \Delta CP{L}_k\rangle }{\sigma }$ (10)

其中，ΔCPL_k表示敲除节点k后特征路径长度的变化， $\langle \Delta CPL\rangle$ 表示所有节点依次敲除后特征路径长度变化的平均值，σ是相应的标准差。Z-score值较大的节点往往在变构信号传导中起关键作用。

3. 结果与讨论

3.1. GNM获得的理论与实验B-Factor的比较

温度因子(B-factor)可以反映生物大分子局部结构的柔性。本文在搭建Kir3.2的GNM模型时，截断半径r_c通过最大化理论和实验B-factor的皮尔森相关系数(Pearson correlation coefficient, PCC)来确定。最大化过程中，截断半径在6.0~13.0 Å内变化，步长为1.0 Å。结果为r_c= 13.0 Å时，理论和实验B-factor的PCC达到了最大0.77。图2显示了在最优的截断半径下的理论和实验B-factor的比较图。由于Kir3.2是由同源四聚体组成的，所以在此只显示了一个亚基的柔性。由图2，GNM可以很好地再现Kir3.2通道结构的柔性。跨膜螺旋(TM1和TM2)和CTD的β片层区域柔性最小，一致于它们紧密堆积的折叠结构，暗示了其高稳定性；选择性过滤器、turret和CTD的β片层环状区域柔性较大，这有助于离子通道的开合运动，与实验结果很好吻合 [1] 。

Figure 2. Comparison between the theoretical (solid line) and experimental (dotted line) B-factors of Kir3.2 chain A

图2. Kir3.2的A链理论(实线)和实验(虚线) B因子的比较

3.2. 基于GNM的Kir3.2残基运动相关性分析

残基涨落的运动相关性可以反映分子功能结构域间的偶合运动，提供对分子功能性运动的洞察。GNM的慢运动模式代表了与蛋白质功能相关的大规模集体运动 [19] ，因此为探究体系功能性运动，选取GNM的前340个慢运动模式(对残基涨落的贡献超过50%)，按照公式(7)计算了残基运动交叉相关性，如图3(b)所示。因Kir3.2是由同源四聚体组成的，所以在此只分析了A和B两个相邻亚基的运动相关性，如图3(a)所示。由图3(a)可见，A链的界面螺旋和B链CTD的多个相邻β片层的环状区域间存在较强的运动正相关性，如图中标号为1的区域所示，这说明它们之间存在强的相互作用，这可促使通道激活过程中CTD的旋转运动；标号为2的区域显示了A和B链的跨膜螺旋间呈现显著的强运动正相关性，这说明不同亚基间的跨膜螺旋存在强相互作用以维持TMD的稳定性，并促进内螺旋的倾斜运动；标号为3的区域反映了结合口袋区域，即A链的中间linker区和B链CTD的3个β片层环的Loop区域之间存在强的运动正相关性，我们推测这有助于两个区域间的变构信号交流，促使两个亚基的CTD产生旋转运动，进而引起通道门的运动 [1] 。

3.3. 基于氨基酸复杂网络模型识别Kir3.2的关键位点

以GNM慢运动模式下残基运动的交叉相关性为边权，对Kir3.2通道建立了加权的氨基酸复杂网络模型。通过依次敲除节点，计算了网络特征路径长度的变化，来识别对蛋白质变构和信号交流起重要作用的关键残基，结果见图4。由图4(a)，有23个残基(功能残基簇的中心)有较高的Z-score值(Z-score ≥ 1.0)，其在结构上的位置分布如图4(b)所示。接下来，将根据现有的实验和理论数据对这些残基的功能进行讨论。

Figure 3. Cross-correlations of Kir3.2 residues calculated based on GNM. (a) Cross-correlations of residues in chain A and chain B. (b) Cross-correlations of the whole Kir3.2 channel residues

图3. 基于GNM计算的Kir3.2残基的运动交叉相关性，(a) 为A链和B链残基的交叉相关性，(b) 为整个Kir3.2通道残基的运动交叉相关性

位于界面螺旋及其附近的功能残基簇的中心为Val72、Thr74、Tyr78、Thr85、Val87、Leu89和Asn94。Tyr78与相邻亚基β_C-β_D环上的Leu229形成氢键，稳定了TMD-CTD结合界面，对通道的生物学功能实现有重要作用 [20] 。位于界面螺旋附近的Asn94与界面螺旋Leu86形成了氢键，稳定了界面螺旋的位置，在控制通道活性中起重要作用 [1] 。其余的残基位于PIP2结合口袋和界面螺旋的中间位置，我们认为它们对配体结合很重要，在促使界面螺旋运动和通道的激活中起关键作用。

位于内螺旋上的功能残基簇的中心为Gly180、Ala185、Val188、Phe192和Lys194。Val188的侧链在通道激活时发生轻微偏转，通过溶剂化作用参与了钾离子的转运 [11] 。Phe192侧链的苯环参与了内螺旋门的形成，处于通道的狭窄部分，对钾离子的跨膜流动有重要作用 [11] 。Lys194参与了与PIP2的结合，对识别PIP2和发送变构信号到TMD起重要作用 [1] 。另外两个残基Gly180和Ala185位于TMD的内螺旋上，我们推测其可能作为变构信号传导的枢纽，对信号向上传到选择性过滤器和turret区域发挥重要作用。

位于孔螺旋上的中心残基是Glu152，该残基紧邻选择性滤器，可能对稳定选择性滤器结构，并帮助其发挥生物学功能有重要的作用。

位于TMD-CTD铰链区附近的中心残基为Gln197和Thr204。这两个残基位于铰链区两端的loop上且在PIP2结合位点附近，可能作为变构信号传输的重要枢纽对Kir3.2的激活和CTD的扭转运动有重要作用。

位于G环门附近的功能残基簇的中心为Gly318。该残基紧邻负责G环门的Met319，可能作为信号交流的枢纽对信号传到Met319并激活G环门起重要作用。另外，Gly318主链羰基和Thr320侧链氧原子位于通道孔，对形成亲水孔允许水合K离子通过起重要作用 [1] [21] 。

最后，识别到的其他关键残基Ser208、Asp217、Asp261、Asn263、Asp270、His283、Gln322和Ser326位于CTD上。目前尚未有实验数据证实它们在Kir3.2变构中的重要作用。CTD是通道功能发挥不可缺少的结构域，这些残基潜在的功能还需实验证实。

(a) (b)

Figure 4. Identification of key residues in chain A based on Z-score. (a) Z-score values of residues; (b) Distribution of the 23 central residue clusters on Kir3.2 channel

图4. 基于Z-score识别到的A链上的关键残基。(a) 残基的Z-score值；(b) 23个中心残基簇在Kir3.2结构上的分布

4. 结论

Kir3.2是一种内向整流型钾离子通道，在神经和心肌系统中发挥重要作用。本工作首先利用GNM对Kir3.2进行建模，成功再现了通道结构的柔性。结果发现，跨膜螺旋和CTD的β片层区域结构比较刚性，一致于其折叠紧密的结构；位于TMD顶端的turret和CTD上的β片层环柔性较大，这有助于离子通道的开合运动，与实验结果吻合。接着，计算了Kir3.2的残基运动相关性，发现不同亚基的跨膜螺旋之间以及CTD的β片层环和linker区域间存在较强的正相关，这暗示了这些区域间的强相互作用，进而说明它们可能对结构的稳定和通道的激活有重要作用。最后，基于加权的氨基酸网络模型预测到了Kir3.2上与实验数据高度吻合的功能性关键残基，另外也识别到了一些尚未被实验证实的关键位点，它们可为未来实验研究提供重要信息。本工作有助于深入理解Kir3.2通道的工作机制，对药物设计具有指导意义。

基金项目

国家自然科学基金项目(31971180, 32271294)。

NOTES

^*通讯作者。