基于弹性网络模型的蛋白质变构路径与关键残基识别研究

李 娇 王韦卜 苏计国

（1）燕山大学理学院，河北省微结构材料物理重点实验室，秦皇岛 066004；2）国药中生生物技术研究院有限公司第六研究室，北京 101111）

变构调节是蛋白质功能调节的一种重要机制，变构效应因子在变构位点的结合，可以通过特定的变构信号传输路径传递到远处的活性位点，使活性部位的局部结构或构象运动发生改变，进而对蛋白质的功能进行调节［1-2］。变构相互作用最早由Umbarger等［3］在1956年提出，他们发现了左旋异亮氨酸对左旋苏氨酸脱氨酶的反馈抑制作用。之后，Changeux等［4-6］在含有几个亚单位的蛋白质体系中，观察到配体结合引起蛋白质结构发生一定的构象变化，由此提出了变构效应。尽管蛋白质的变构调节在结构生物学领域的研究进行了大半个世纪，但是蛋白质变构调节的物理机制以及相应变构路径的有效识别，仍是一个有待深入研究的重要科学问题。变构调节在蛋白质功能反馈和调控机制中起着关键性作用，许多蛋白质的活性及其功能的调节都是通过变构效应来实现的。因此，蛋白质变构路径的有效识别，不仅有助于揭示蛋白质功能调控的分子机制，还可以为变构药物的开发或纳米生物传感器件［7-9］的开发提供理论依据。

分子动力学模拟（molecular dynamics，MD）已广泛应用于蛋白质的变构效应研究中，通过对构象运动和运动关联性的分析等来揭示蛋白质体系的变构路径和关键残基［10-15］。Roy 等［16］运用MD 模拟结合运动相关性分析，对人类鼻病毒（HRⅤ）的变构调控过程进行了研究，他们通过对不同残基之间径向运动关联性的分析，揭示了抗病毒化合物WⅠN 52084 对HRⅤ衣壳蛋白“呼吸”运动的长程变构调控机制。Bowerman 和Wereszczynski［17］基于MD模拟轨迹，利用残基运动相关性分析、互信息关联分析以及图论的方法来识别蛋白质体系的变构网络，利用该方法他们详细分析了凝血酶变构位点与催化中心之间的长程变构信息传播路径。MD模拟虽然在蛋白质变构效应研究中得到了成功应用，但是，很多情况下，蛋白质的变构调控并不一定伴随明显的构象运动，而是通过刚性区域来耗散变构所产生的应力，这些刚性区域不会发生大的构象变化，而是产生较大的力信号［18］，常规的平衡态MD模拟很难识别到这些力信号。因此，一些研究者通过对体系施加外界的扰动，利用MD模拟结合力分布分析（force distribution analysis，FDA）方法来揭示蛋白质内力信号的长程传播过程，进而获得体系的变构调控路径［19-23］。Stacklies 等［24］对肌联蛋白Ⅰ27结构域的长程力学变构调控进行了研究，对Ⅰ27的N端、C端施加300 pN的恒力进行拉伸，运用MD 结合FDA 的方法，计算了不同残基间的力分布，得到了蛋白质力信号的传导路径。随后，该小组运用同样的方法，对MetJ 蛋白进行了研究，识别了力信号从变构位点传递到活性位点的变构路径［18］。2013年，Palmai等［25］利用MD结合FDA 的方法，研究了3-磷酸还原酶（hPGK）催化功能执行过程中的变构调控机制，揭示了变构信号从底物结合位点到功能性运动铰链区的传导路径。这些研究表明，基于MD模拟轨迹，通过对蛋白质体系内部力分布的分析来研究蛋白质变构调控路径是一种有效的分析方法。但是，利用MD模拟来计算蛋白质内的力分布存在一定的局限性。首先，由于蛋白质体系的热涨落，力学信号常常被热噪声所淹没，需要设计一定的去噪音手段才能获取有价值的力学信号。其次，为了去除热噪声，获得具有统计意义的力分布结果，需要进行大量轨迹的模拟，计算非常耗时。而弹性网络模型（elastic network model，ENM）不需要考虑热噪声的影响，并且对于体系的配分函数和热力学性质可以进行解析求解，大大的减少计算量，提高计算效率。

ENM 方法是研究蛋白质结构和动力学性质的有效方法［26-27］，在蛋白质结构-功能关系研究中得到了广泛应用，包括蛋白质结构稳定性的分析、关键残基的识别、蛋白质结构域的分解、大幅度功能性构象运动的分析等［28-35］。Wang等［36］基于ENM，通过力分布的计算来分析蛋白质体系对外力的响应过程，研究了GB1 蛋白质力学性能的变构调控机制。对GB1蛋白的N端、C端施加20 pN的外力进行拉伸，分别研究结合配体与未结合配体情况下蛋白质体系的形变情况和力分布情况。研究发现，配体的结合能够通过长程变构效应调控抗力部位的机械性能，显著增强GB1的力学稳定性。

本工作将基于ENM 的力分布计算方法用于蛋白质变构路径和关键残基位点的识别研究，基于该方法，分析了人类磷酸甘油酸激酶（human phosphoglycerate kinase，hPGK）和蛋白质酪氨酸磷酸酶（protein tyrosine phosphatase，PTP）PDZ2结构域的变构路径和关键残基。对两种蛋白质的变构位点施加外力进行拉伸，通过计算外力作用下蛋白质内部的力分布，来识别与力承载区域形变相耦合的残基位点，并研究力信号在蛋白质结构内的传播路径。

1 理论与方法

ENM 将蛋白质的空间结构简化为一个由大量节点和节点之间的弹簧所构成的弹性网络。将蛋白质体系中的每一个残基简化为网络中的一个节点，用Cα 原子来代替，其他的原子忽略。选取一个几何距离作为截断半径rc［26］。如果两个残基间的距离小于该截断半径，认为它们之间存在相互作用，则用弹簧连接；若两残基间距离大于该截断半径，则认为不存在相互作用。所有连接的弹簧的弹性系数均相同。在本工作中，弹簧的弹性系数设为1.0 N/m，选取的截断半径为8.0 Å。由整个体系中弹簧的连接情况，可得到体系的势能为

其中，γ为弹簧的弹性系数；ΔR为蛋白质内各个氨基酸偏离其平衡位置的位移，是3N维列向量，N为残基的个数；右上角的T表示向量的转置；H为3N×3N的Hessian矩阵，该矩阵可表示为

上式中的hij为3×3的子矩阵，i、j分别表示第i、j个残基。

当i≠j时，有

当i=j时，有

其中，矩阵元为势能函数V对坐标x、y求二阶偏导，表达式为

在本工作中，由于对蛋白质关键位点施加外力进行拉伸，体系的总势能Vtotal由原来的弹性势能V变为弹性势能与外力做功之和［33，37］，即：

其中，F表示施加在蛋白质上的外力，为3N维列向量。对i、j残基施加外力进行拉伸，则 外 力F矩 阵 的 构 建 为 ：F=(000…Fix Fiy Fiz…Fjx Fjy Fjz…000)T1×3N。将体系的总势能在平衡位置求梯度，可得到体系达到平衡后，各个残基偏离平衡位置的位移为：

上式中，H-1矩阵的矩阵元为[H-1]ij=表示H-1矩阵的i行j列的矩阵 元，λkk为Hessian 矩 阵 的 第k个 本 征 值，U为Hessian 矩阵的本征向量，[Uk]i表示第k个本征向量中的第i个元素。残基位移所导致的体系内弹簧形变可以表示为［38-39］：

ΔrM×1为残基间弹簧的形变，M表示组成弹性网络的弹簧总数，AT是残基偏移平衡位置位移在弹簧形变矢量上的投影。残基间相应的内力分布可以表示为：

K是弹簧弹性系数所构成的M维的对角矩阵，f是体系中所有弹簧承受的力，为M×1维的向量。

2 结果与讨论

2.1 hPGK结构域中的关键残基及变构路径

hPGK 在生物体内参与糖酵解过程，是催化1，3 二磷酸甘油酸转化为3-磷酸甘油酸的关键酶，可产生1分子的ATP，在糖酵解的过程中发挥着不可或缺的作用［25，40］。hPGK 的三维结构（PDB 代码2XE7）由两个大小相近的结构域组成（图1a）。研究证实，hPGK两个结构域的开合运动在催化功能的发挥中起着重要作用，而Arg38、Glu192、Asn336、Gly394 残基是两个结构域打开与闭合过程中的关键残基［25］。其中，残基Arg38 和Asn336直接参与了底物的结合。残基Glu192 和Gly394 位于体系两个结构域之间的铰链区，直接调控着两个结构域的相对运动。大量的实验结果表明，底物的结合会引起铰链区的构象调整，从而导致hPGK的两个结构域发生大幅度的构象变化，从开放构象转变到闭合构象，进而催化底物的化学反应，完成其生物学功能［41-44］。底物的结合如何调控体系铰链区的构象运动，变构信号如何从底物结合位点有效传递到体系的铰链区是值得研究的重要问题，变构路径的识别研究有助于更好地揭示hPGK发挥功能的分子机制。

在本工作中，首先基于hPGK的三维结构构建弹性网络，然后对残基Arg38 和Asn336 施加1 pN的外力进行拉伸，拉伸方向为两残基位点连线的方向。进而，通过公式（9）和（10）计算蛋白质在外力作用下的形变情况和相应的内力分布情况。计算得到体系内力的分布情况如图1b 所示。图中蓝色到红色表示力的分布从小到大，受力越大表明相应残基对于力信号的传导越重要，相应位点为力信号变构传导的关键位点。

Fig.1 Three-dimensional structure of hPGK(a)and internal force distribution of hPGK in response to the exertion of external forces(b)

将图1b 中不同残基之间弹簧的受力情况映射到蛋白质的三级结构上，并根据不同残基间力的分布情况来识别力信号的传播路径（图2）。将外力所施加的两个残基位点Arg38 和Asn336 在图中用紫色球标出，为了更清晰的显示力信号的传播路径，将受力大于0.05 pN 的弹簧显示在图2 中。图中残基间连线的粗细表示其承受的内力大小，连线越粗，表示所承受的力越强，对力信号传导也越重要；反之，连线越细表示受力越弱。由残基间内力分布情况，将底物结合位点与铰链区之间受力最大的弹簧相连（图2 红色连线），得到力信号的传导路径。根据该方法，从Arg38出发，可以得到两条到达铰链区的变构路径：Arg38→Thr393→Gly394（铰 链 区） 和Arg38→Thr393→Ala41→Leu188→Glu192（铰链区）。这两条变构路径用红色箭头显示在了图2中，其中所涉及到的关键残基用黄色球标出。采用类似方法，从Asn336 施力位点出发，也可以得到两条传播到铰链区的力信号传导路径：Asn336→Gly371→Thr393→Gly394 （铰 链 区） 和Asn336→Gly371→Thr393→Ser392→Glu192 （ 铰链），传导路径在图2 中用蓝色箭头标出。变构路径中残基间的内力值见表1。

Fig.2 Allosteric pathways in hPGK identified by our method

Table 1 Internal forces between residues in the allosteric pathway of hPGK in response to the exertion of external forces

Szabo 等［45］通过残基突变、酶反应动力学实验以及结构分析，研究了底物结合对hPGK功能性开合构象运动的变构调控过程，并识别了变构信号从底物结合位点到构象运动铰链区传播过程中的关键残基，发现Arg38、Lsy219、Asn336、Glu343等残基在结构域开合过程中发挥着重要作用。另外，Palmai 等［25］利用MD 模拟与FDA 相结合的方法，研究了底物结合对铰链区构象转变的变构调控作用，得到了两条长程变构路径，分别是从结合位点残基Arg38 传递到Thr393 再到铰链区残基Glu192、Gly394，以及从结合位点残基Asn336 到Gly371、Ser398、Ser392 再传递到铰链位区Glu192 和Gly394。本文的计算结果与上述实验和MD模拟结果能够很好地吻合。

2.2 PTP PDZ2结构域的关键残基及变构路径

PTP是生物体内的关键酶［46-48］。它与蛋白酪氨酸激酶协同作用，共同维持着酪氨酸磷酸化的平衡，进而参与体内的细胞生长、新陈代谢等活动。PDZ 结构域存在于许多蛋白质内，被认为具有蛋白质定位、信号传输等功能［46，49-50］。PTP 包含有5个PDZ 结构域，其中，第二个结构域由6 条β 链（β1～β6）和2条α螺旋（α1、α2）组成，含有96个氨基酸，属于典型的PDZ 结构［51-52］。图3a 中不同的颜色表示不同的局部二级结构。β2与α2螺旋之间存在一个配体结合口袋，实验研究表明，配体的结合能够长程变构调控远端的α1和β1区域与其他蛋白质的相互作用，变构信号如何从配体结合口袋传递到远端区域是有待深入研究的重要问题，有助于更好揭示PDZ2功能的分子机制。

Fig.3 The three-dimensional structure of PTP PDZ2(a)and internal force distribution of PTP PDZ2 in response to the exertion of external forces(b)

本研究基于PTP PDZ2 的三级结构（PDB 代码3PDZ）构建相应的弹性网络，然后对配体结合口袋区域的变构位点残基Gly19和His71施加1 pN外力进行拉伸，拉伸方向为两残基的连线方向。运用公式（10）计算蛋白质体系在外力作用下的内力分布情况（图3b）。图3b中的红色区域对应于承受内力较大的残基，是力信号传导过程中的关键残基。

为了更清晰地分析力信号的传导路径，将计算得到的力分布显示到蛋白质的三级结构上（图4）。图中仅画出了受力大于0.01 pN 的弹簧，显示为绿色连线。连线的粗细表示相应弹簧承受的内力情况，连线越粗表示弹簧所承受的内力越强。由残基间内力分布情况，将施力位点与远端β1区域间受力最大的弹簧相连，得到两条特异性变构路径（图4红色连线）。一条变构路径从β2的施力位点开始，传递到α1螺旋的N 端，最终到达β1-β2之间的loop区，即Gly19、Ser21（β2）→Ⅰle41、Gln43、Gly45（α1）→Ser17（β1-β2之间的loop区），如图4右侧上图中红色箭头所示，变构路径中的关键残基用黄色球显示在图中。另一条力信号的传导路径是α2螺旋的施力位点开始，穿过β2、β3、β4、β6，最终到达β1，即His71 （α2） →Ⅴal22、Ⅰle20 （β2） →Ⅴal37（β3）→Gly55（β4）→Gln93、Lys91（β6）→Gly4（β1），如图4右侧下图中红色箭头所示，变构路径中的关键残基用橙色球显示在图中。变构路径中残基间的内力值见表2。

Fig.4 Allosteric pathways in PTP PDZ2 identified by our method

Table 2 Internal forces between residues in the allosteric pathway of PTP PDZ2 in response to the exertion of external forces

Kong等［48］基于MD模拟，采用残基间相互作用的关联性分析，研究了PTP PDZ2体系的长程变构调控机制，识别出两条变构路径，一条路径始于β2上的变构位点长程传递到α1螺旋，另一条路径始于α2上的配体结合位点，垂直穿过β2、β3、β4、β6最终到达β1链的N 端。本研究利用基于ENM 的力分布分析方法得到的变构路径与Kong 等利用MD模拟得到的结果能够很好地吻合。相对于MD 模拟，基于ENM 的力分布分析方法计算更为简单、快速。

3 结 论

本文利用基于ENM 的力分布计算方法，通过分析蛋白质对外力的响应，研究了hPGK、PTP PDZ2 两个蛋白质体系力信号的变构调控机制，成功识别出体系中特异的变构路径和关键残基，与实验和其他计算模拟结果能够很好地吻合。该方法将蛋白质的三级结构简化为弹性网络，通过对变构位点施加外力，计算蛋白质体系在外力作用下的形变情况和体系内部的力分布情况，来识别力信号传播过程中起重要作用的关键残基和相应的变构传播路径。相对于MD 模拟方法，基于ENM 的力分布分析方法具有明显的优势，MD模拟需要耗费大量的模拟时间，并且需要利用有效的统计分析方法，去除噪音的干扰，来识别变构信号，而ENM 方法可以通过解析求解来获得体系的力分布情况，不需要进行数值模拟和轨迹采样，也不包含噪声，计算更加简洁快速，能够用于超大蛋白质体系的研究。因此，本文为蛋白质的变构效应研究提供了一种简单、有效的方法。但是，需要指出的是，该方法也存在一定的不足，由于ENM 是一种对蛋白质结构进行高度简化的粗粒化模型，每个氨基酸简化为一个点，很难考虑不同配体结合、残基的不同突变、翻译后修饰等相互作用的精细差别，要更为精确的研究这些不同类型相互作用的影响，还需要结合更精细的全原子模型。