基于弹性网络模型和微扰响应扫描的PKC-θ的动力学和关键残基研究

Study on the Dynamics and Key Residues of PKC-θ Based on Elastic Network Model and Perturbation Response Scanning

摘要: PKC-θ (protein kinase C-theta)激酶能够对底物蛋白中的丝氨酸/苏氨酸进行磷酸化修饰，并参与T细胞活化、信号转导及免疫应答等重要生理过程，与许多免疫疾病相关，因此研究PKC-θ的结构动力学及其关键位点具有重要意义。在本工作中，我们对PKC-θ构建了弹性网络模型，并在慢运动模式下识别出了结构中重要的功能区域。通过残基的运动相关性分析发现：N端小叶和C端大叶的运动是负相关的，这一致于两叶间的“开合运动”，有助于结合配体的稳定；此外，催化相关残基与DFG基序的运动呈正相关，这有利于催化功能的发挥。接着，利用微扰响应扫描方法来识别变构关键残基，结果发现，高敏感性残基对配体结合及酶的催化活性很重要。这项工作有助于加强对PKC-θ结构动力学的理解，并为设计激酶抑制剂提供重要信息。

Abstract: PKC-θ (Protein Kinase C-theta) can phosphorylate the serine/threonine of substrate protein, participating in important physiological processes such as T cell activation, signal transduction and immune response, and related to many immune diseases. Thus, it is of great significance to study its structural dynamics and key residues. In this work, we construct PKC-θ’s elastic network model, and the important functional areas are identified based on GNM’s slow motion modes. According to interresiduemotion correlations, it is found that there exists a negative motion correlation between N-lobe and C-lobe, consistent with their open-closed motion, which helps for the bound ligand’s stability; in addition, there is a positive motion correlation between DFG motif and catalysis related residues, which is conducive to the catalytic function exertion. Then, the perturbation response scanning method is used to identify the key residues for allostery, and the results suggest that the key residues of high sensitivity are important for ligand binding and enzyme’s catalytic activity. This work helps to strengthen the understanding of structural dynamics of PKC-θ and provides important information for kinase inhibitor design.

1. 引言

新型蛋白激酶C (novel protein kinase C, nPKC)是一类能够磷酸化修饰底物丝氨酸/苏氨酸的激酶[1] [2]。nPKC包括δ、ε、η、υ和θ等成员，参与多种细胞信号转导途径和生理过程[3]。其中，PKC-θ主要分布在免疫细胞中并参与T细胞的活化和免疫应答，其已被证明与多种疾病的发生和发展相关[4] [5]，是各种T细胞介导的疾病的有吸引力的治疗靶点。

PKC-θ由N端调节域和C端激酶结构域(催化域)组成。激酶结构域自身的磷酸化对其催化活性至关重要[6]，一系列磷酸化事件使其成熟[7]，进而具有使底物蛋白磷酸化的功能。如图1所示，激酶结构域为双叶折叠，N端小叶(N-lobe)和C端大叶(C-lobe)间结合一个ATP分子，该ATP在磷酸化过程中充当磷酸供体[8]。激活环(Active-loop)中含有保守的催化组件DFG基序(Asp-Phe-Gly)，其作用是将ATP定位并进行磷酸基转移。激活环通过αC螺旋与N-lobe相连，环上Thr536的磷酸化会导致激酶构象转变。αC螺旋是N-lobe疏水口袋的一部分，在所有激酶中发挥重要的催化和调节功能[8] [9]。富含甘氨酸的环(G-rich loop)，即三个高度保守的甘氨酸GxGxxG基序，是激酶核心中最动态的环。该环的主要功能是定位腺嘌呤环，并在磷酸基的转移中发挥重要作用，其突变与多种疾病有关[10] [11]。由于PKC-θ是很多疾病的治疗靶点，所以其结构动力学及关键位点研究受到了广泛关注。

实验方面，Graham等人用单克隆抗体和磷酸化特异性抗体血清方法确定了PKC-θ的磷酸化位点Thr538、Ser676和Ser695，它们分别位于激活环、转角基序和疏水基序上，其中激活环的磷酸化至关重要[12]。理论方面，Seco等人利用分子动力学(molecular dynamics, MD)模拟发现激活环的磷酸化诱导了ATP结合位点的打开，揭示了激活环的构象与激酶的催化间的关系[6]。Iqbal和Hatmal等人利用MD模拟和分子对接筛选出针对不同疾病的有效激酶抑制剂[13] [14]。众所周知，对大分子进行MD模拟是非常耗时的，因此人们提出了粗粒化的弹性网络模型(Elastic Network Model, ENM)用于研究生物大分子结构动力学与功能的关系[15]。两种传统的ENM模型是高斯网络模型(Gaussian Network Model, GNM) [16]和各向异性网络模型(Anisotropic Network Model, ANM) [17]。Atilgan等人提出了微扰响应扫描(perturbation response scanning, PRS) [18]方法，该方法与ENM相结合来评估蛋白质残基因外界扰动所产生的响应，目前该方法已被证明在蛋白质变构调控位点识别中非常有效[19]。

本研究使用弹性网络模型和微扰响应扫描方法分析了PKC-θ的结构动力学与功能间的关系，并识别出了在变构通信中发挥重要作用的关键残基，加深了对PKC-θ的激活及催化机制背后动力学的理解。

Figure 1. Structure of the kinase domain of PKC-θ. The active-loop (orange), αC and αG helix (violet), G-rich loop (yellow) turn motif (bright orange, the phosphorylation site shown in stick model), hydrophobic motif (blue, the phosphorylation site shown in stick model) and DFG (orange stick at active-loop) are all labeled

图1. PKC-θ激酶结构域的结构。激活环(橙色)、αC和αG螺旋(紫罗兰色)、G-rich loop (黄色) turn motif (亮橙色，磷酸化位点显示为棍棒)、hydrophobic motif (蓝色，磷酸化位点显示为棍棒)和DFG (橙色棍棒)均已标注

2. 研究方法

2.1. 弹性网络模型

从蛋白质数据库(Protein data bank, PDB)中获取PKC-θ激酶结构域的晶体结构(PDB ID: 4Q9Z)，其分辨率为2.6 Å，选取A链作为研究对象，用于搭建ENM模型。在ENM中，该激酶的结构被粗粒化为一个由弹簧连接的节点网络，其中残基被粗粒化为由Cα原子表示的节点，如果两个节点间的距离小于某一截断半径，则用弹簧连接，这表示残基之间有相互作用。在GNM中，网络的总势能为：

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

其中，γ是弹簧的弹性系数，∆R代表N个Cα原子偏离其平衡位置的位移，N是该蛋白中残基的个数，E是单位矩阵，$\otimes$ 为矩阵的直积，Γ为N × N阶的Kirchhoff矩阵，其元素为：

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

其中，R_ij是节点i和j间的距离，r_c为截断半径。残基的均方涨落和残基间涨落的交叉相关性分别为：

$\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 \sum _{k=2}^N{\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 \sum _{k=2}^N{\lambda }_k^{-1}{\left[u_k\right]}_i}{\left[u_k\right]}_j$ (4)

其中，k_B为Boltzmann常数，T为绝对温度。

与GNM模型相比，ANM模型不但可以提供残基运动的幅度信息，而且还能够提供残基运动的方向信息。在ANM中，网络的总势能为：

$V_{\text{ANM}}=\frac{1}{2}\gamma {\displaystyle \sum _{i,j}^N{\left(R_{ij}-R_{ij}^0\right)}^2}$ (5)

其中，R_ij表示残基i和j之间的瞬时距离，$R_{ij}^0$ 表示在平衡位置时残基i和j之间的距离。对公式(5)求二阶偏导数计算出Hessian矩阵H，该矩阵分解得到的特征向量表示蛋白质的运动模式。H的元素h_ij是3 × 3的子矩阵。当i≠ j时，h_ij为：

$h_{ij}=\left[\begin{array}{ccc}\frac{{\partial }^{\text{2}}V}{\partial X_i\partial X_j}& \frac{{\partial }^{\text{2}}V}{\partial X_i\partial Y_j}& \frac{{\partial }^{\text{2}}V}{\partial X_i\partial Z_j}\\ \frac{{\partial }^{\text{2}}V}{\partial Y_i\partial X_j}& \frac{{\partial }^{\text{2}}V}{\partial Y_i\partial Y_j}& \frac{{\partial }^{\text{2}}V}{\partial Y_i\partial Z_j}\\ \frac{{\partial }^{\text{2}}V}{\partial Z_i\partial X_j}& \frac{{\partial }^{\text{2}}V}{\partial Z_i\partial Y_j}& \frac{{\partial }^{\text{2}}V}{\partial Z_i\partial Z_j}\end{array}\right]$ (6)

当i = j时，h_ij为：

$h_{ii}=-{\displaystyle \sum _{i\ne j}{h}_{ij}}$ (7)

在ANM中，残基的均方涨落和残基间涨落的交叉相关性为：

$\langle \Delta R_i\cdot \Delta R_i\rangle =\frac{k_{B}T}{\gamma }\left(H_{3i-2,3i-2}^{-1}+H_{3i-1,3i-1}^{-1}+H_{3i,3i}^{-1}\right)$ (8)

$\langle \Delta R_i\cdot \Delta R_j\rangle =\frac{k_{B}T}{\gamma }\left(H_{3i-2,3j-2}^{-1}+H_{3i-1,3j-1}^{-1}+H_{3i,3j}^{-1}\right)$ (9)

归一化的交叉相关性为：

$C_{ij}=\frac{\langle \Delta R_i\cdot \Delta R_j\rangle }{\sqrt{\langle {\left(\Delta R_i\right)}^2\rangle \cdot \langle {\left(\Delta R_j\right)}^2\rangle }}$ (10)

C_ij的取值范围为[−1, 1]，1表示残基的运动方向完全相同，−1表示运动方向完全相反，0则表示两残基的运动不相关。

2.2. 微扰响应扫描方法

为了探究影响蛋白质构象变化的关键残基，Atilgan等人根据线性响应理论(linear response theory, LRT) [20]提出了PRS方法[11]。该方法可以根据胡克定律计算出所有残基在受到外力时的位置变化量。在PRS中，依次对残基施加外力F，通过观察所有残基的位置变化来判断该残基的重要性。其中施加在残基i上的力F表示为

$F^i={\left(000\cdots \Delta F_x^i\Delta F_y^i\Delta F_z^i\cdots 000\right)}^{\text{T}}$ (11)

其他残基的响应：

$\Delta R^i={{\rm H}}^{-1}{F}^i$ (12)

其中，∆Rⁱ表示所有残基在残基i被扰动时的位置变化量。

对每个残基分别施加x方向、y方向、z方向、x和y方向、x和z方向、y和z方向、x和y和z方向的7个外力。PRS方法通过对所有残基进行扰动，得到一个N  ×  N阶的位移响应矩阵P，矩阵的第k列表示扰动残基k时其他残基的响应。利用对角线元素对矩阵P中的每一列进行归一化，归一化后的列平均值被定义为相应残基的敏感性。由此可以得到所有残基的敏感性曲线，曲线中的极大值所对应的残基被认为在蛋白质变构中起关键作用[21]。

3. 结果与讨论

3.1. 基于GNM的实验和理论B因子的比较

在X-ray结构中，Debye-Waller因子(又称温度因子或B因子)常被用于衡量原子的柔性。本文用Cα原子的B因子来表示其所在残基的柔性信息。通过最大化理论和实验B因子之间的皮尔森相关系数(Pearson correlation coefficient, PCC)来优化GNM中的截断半径，以构建最佳模型。在本工作中，截断半径r_c的寻优范围为：6.0~13.0 Å，步长为1.0 Å。当r_c = 8.0 Å时，GNM计算的理论B因子与实验B因子的PCC达到最大值0.71，故本工作中r_c = 8.0 Å。图2显示了GNM在截断半径为8.0 Å时得到的理论和实验B因子的对比图，可以看出GNM能较好的再现PKC-θ的柔性。

Figure 2. Comparison between the experimental (solid line) and computed (dotted line) B factor of PKC-θ

图2. PKC-θ的实验 (实线)和理论 (虚线) B因子的比较

3.2. 基于GNM的PKC-θ慢运动模式分析

研究表明，与蛋白质功能有关的慢运动模式通常是大规模的集体运动[22]。为研究PKC-θ的功能动力学特性，本文计算了PKC-θ在第一慢运动模式下的涨落，如图3所示。图中峰值主要对应于G-rich loop和C末端的转角基序及疏水基序，这些区域多为loop结构，故柔性较大。有趣的是，几个涨落较小的区域与PKC-θ的功能密切相关。区域1 (残基421-441)对应于αC螺旋，该螺旋具有催化和调节功能[8]，其中保守的谷氨酸(Glu428)与赖氨酸(Lys429)形成的氢键有利于催化功能的发挥；区域2 (残基461-473)是连接N-lobe和C-lobe的铰链区，其中含有一个ATP结合位点Leu461，该铰链区对PKC-θ的催化活性和结构稳定起重要作用[8]；区域3 (残基506-526)的部分残基位于激活环中，该环在激酶结构域中较为保守，其中包含了重要的催化组件DFG基序，另外该区域中的Ala521和Asp522是ATP结合位点[23]。

Figure 3. Mean square fluctuation distribution of the residues in the first slow motion mode

图3. 第一慢运动模式下残基的均方涨落分布

3.3. 基于ANM的PKC-θ运动相关性分析

为了研究PKC-θ在运动过程中不同残基的相关性，本文通过公式(10)计算出所有残基对的交叉相关性，如图4所示。因为慢运动模式对应着蛋白质的功能性运动，所以选取ANM (r_c = 15.0 Å，理论与实验B因子之间的PCC为0.69)的前19个低频运动模式(对残基涨落贡献大于50%)进行运动相关性分析。从图4看出，N-lobe (左下角虚线方框)和C-lobe (右上角虚线方框)内的残基运动多呈正相关，两叶间的残基运动多呈负相关，这与其“开合运动”相一致[24]。图中区域1显示αC螺旋与Ala483和Phe484等残基具有强烈的运动正相关性，它们均与ATP的结合有关[25]。此外，对催化起重要作用的Asp466和Asn471与催化组件DFG基序也存在强的运动相关性，它们一起保证催化功能的稳定发挥。区域2表明C末端的疏水基序与αC螺旋呈正相关，两者的相互作用对激酶活性具有调节作用[26]。区域3反应了转角基序与β1间的运动相关性，由于其空间上的近邻，所以它们具有强正相关性。

3.4. 微扰响应扫描分析

PRS方法结合了ENM和LRT来评估在蛋白质上施加外力扰动对其构象的影响。本工作对PKC-θ中的残基逐个微扰得到响应热图和残基敏感性曲线，如图5(a)所示；残基敏感性值映射在PKC-θ结构上，如图5(b)所示。

经过分析，发现PKC-θ中敏感性较高的残基大多与ATP或底物结合有关。在图5(b)中，区域1是G-rich loop，该环的功能是定位腺嘌呤环和ATP的磷酸基，敏感性较高。区域2是位于αC螺旋和β4间

Figure 4. The cross-correlation map for PKC-θ

图4. PKC-θ的残基运动相关性热图

(a) (b)

Figure 5. (a) PRS heat map of PKC-θ, with sensitivity curve at the top. (b) Mapping of sensitivity values on the structure of PKC-θ

图5. (a) PKC-θ的PRS热图，顶部是敏感性曲线；(b) 敏感性值在PKC-θ结构上的映射

的小螺旋和一段loop (αC-β4 loop)，它们是激酶发挥功能的多种调节的枢纽，也是调节激酶活性的复杂相互作用中的热点。Yeung等人的实验表明，该位置的突变可改变家族保守的和特异性的调节机制[26]。αC-β4 loop被认为是真核激酶中调节相互作用的保守但研究不足的热点，本文对其敏感性的探究为理解其结构动力学–功能间的关系提供了一个新见解，并为研究新型蛋白激酶抑制剂提供重要信息。区域3是αG螺旋，虽然它远离底物结合位点，但仍与底物相互作用有关[27]。敏感性较高的区域4中包含了转角基序的磷酸化位点和疏水基序的磷酸化位点，对激酶催化活性的调节非常重要。总的来说，识别出的敏感性高的残基与激酶的变构和催化活性相关，这也说明了PRS方法的有效性。

4. 结论

PKC-θ是一种能够使底物蛋白发生磷酸化的激酶，在蛋白质的翻译后修饰中发挥重要作用。首先，构建了PKC-θ在最优截断半径下的高斯网络模型，实验和理论B因子的PCC达到0.71，并且在第一慢运动模式下较好的识别出了PKC-θ中的功能性铰链区和柔性区域。对PKC-θ的运动相关性分析发现N-lobe和C-lobe之间存在较强的相互作用，一致于激酶的“开合运动”，这有助于结合配体的稳定。接着，通过PRS方法对PKC-θ进行微扰，可以识别出关键残基区域，这些区域被证明与该蛋白质的配体结合及催化活性相关。这些关键残基的识别可以为进一步的实验或理论研究提供重要的方向和依据。本工作有利于加深人们对PKC-θ催化活性成熟及其变构机制的理解，并为设计新的蛋白激酶抑制剂提供有用信息。

基金项目

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

NOTES

^*通讯作者。