氧气和一氧化碳在人血红蛋白迁移过程研究

彭雪，芦琛璘，卢滇楠

（化学工程联合国家重点实验室，清华大学化学工程系，北京100084）

引 言

气体分子与生物大分子的结合与解离对于生命活动具有重要的意义。血红素分子和氧气（O2）的结合与解离是血红蛋白发挥生理功能的物理基础[1]。一般认为，一氧化碳（CO）和氧气具有共同的结合位点，从而形成竞争机制。然而，这两种气体如何从血液进入血红蛋白内部的血红素结合位点，又是如何脱离该位点而重新回到血液之中一直存在着比较大的争论。因此，从原子水平上去理解血红蛋白的载氧功能及其机制是一个巨大的挑战[2-5]。

人血红蛋白是由两条α 链和两条β 链组成的四聚体蛋白。氧合血红蛋白为R 态，脱氧血红蛋白为T 态[6]。每个亚基上均有一个包埋在蛋白质内部的血红素分子，可以与O2或者CO 结合。血红素分子被包埋在蛋白内，因此解析气体分子从溶剂扩散到结合位点的迁移路径成为理解血红蛋白功能机制的关键[7]，同时也是人工血液替代品理性设计的重点。1966 年，Perute 等[8]依据人血红蛋白晶体结构，发现组氨酸可能是影响气体分子迁移的关键氨基酸（histine gate），提出了一条气体迁移的组氨酸通道。在高压Xe 气氛下，用X 射线晶体衍射确定了可能与Xe 结合的位点。这些位点周围氨基酸形成的空穴可能组成了气体分子到达活性位点所经历的迁移途径[9]。近年来，随着高分辨瞬时检测技术的发展，研究者对于气体在血红蛋白内部迁移的动力学进行了解析。例如，Lepeshkevich 等采用激光动力学光谱对人体R 态氧化血红蛋白α 单链进行了分析[10]，进而又采用了纳秒激光闪光光解技术研究了O2迁移过程[11]，得出了O2扩散进入以及离开蛋白的完整动力学路径，并且观察到了除Xe位点之外的其他O2停留位点。然而目前的实验研究仅能够给出人血红蛋白内部空腔中气体分子的大致分布，无法对扩散通道进行准确辨认。更为重要的是，目前尚无法实验观察到气体迁移过程途径中因蛋白质柔性产生的瞬时通道，更无法对其功能予以描述[10]。人血红蛋白的瞬时通道主要存在于蛋白质内部，是由疏水性氨基酸组成。迄今为止，实验唯一观察到的是远端组氨酸侧链因运动而打开的组氨酸大门[12-13]，以及由Leu-29、Phe-33、Phe-43、Phe-46、His-64和Val-68组成的通道，研究表明该扩散距离最短的通路可能是气体分子迁移的主要通道[9-13]。

作为实验研究的有益补充，分子动力学模拟日益得到研究者的关注。由于人血红蛋白结构复杂，目前研究主要集中在变构效应[14]。相比而言，肌红蛋白因其结构简单，是目前O2和CO 迁移路径研究的热点[15-24]。例如，Yu 等[17]在2015 年提出肌红蛋白内完整的CO 迁移动力学示意图，这与实验上普遍认为的经过血红素上方组氨酸大门的单一路径不符。随着计算能力增强，近年来出现了对人血红蛋白的分子模拟研究。例如，采用全原子分子动力学模拟对O2离开R 态人血红蛋白的α 单链进行分析，提出四条没有受到其他亚基阻挡的通道，但该研究并没有发现近端空腔[11]。Lepeshkevich 等[10]采用PELE（protein energy landscape exploration）方法，研究了CO 离开T 态及R 态人血红蛋白的过程，发现T态和R 态的主要区别是T 态时E 螺旋的倒塌部分堵住了组氨酸大门，而在肌红蛋白中发现的XE1 位点并没有出现在人血红蛋白中。接着研究者对T态人血红蛋白进行全原子分子动力学模拟，找到了O2离开以及进入α 链和β 链的多条通路，而不再仅仅是实验观察到的以组氨酸为出口的唯一通道[25-26]。除此之外，其他类血红蛋白也多涉及此过程[27]，如Bustamante 等[28]曾对O2在截断血红蛋白的扩散行为进行详细分析。

马尔科夫状态模型是一种后处理统计学模型，其通过构建平衡概率矩阵和转移概率矩阵，可以从有限的分子动力学模拟数据中演化出与蛋白结构变化相关的热力学和动力学信息，基于大规模分子动力学模拟构建马尔科夫状态模型已逐渐成为研究生物大分子构象变化的重要工具[17,24,27]。总体来看，人体血红蛋白是最复杂、含量最多、最有研究意义的血红蛋白[29]，而目前还没有对O2和CO扩散进入R态迁移过程的报道。因此本论文的目的是运用分子动力学模拟结合马尔科夫状态模型的方法对该过程进行分析，以加深对人体血红蛋白生理功能的认识，为后续人工血液替代品的理性设计提供参考。

1 模型和方法

1.1 模型

血红蛋白的晶体结构源自PDB 数据库（PDB：1hho），为R 态人血红蛋白[30]，去掉结构中的结晶水分子得到独立的α 链和β 链。分子动力学模拟在Gromacs5.12 平台上进行，采用GROMOS96 43a1 力场，体系压力为0.1 MPa，温度为293 K。向边长为6.5 nm 的正方形水盒子中加入溶剂水分子8629 个，水分子采用SPC/E 模型[10]，O2和CO 的力场参数分别源自文献[25]和文献[31]，忽略气体分子的四极矩运动[21]。为了提高模拟体系取样精度，模拟盒子中O2和CO 的分子数分别为10 个和50 个，高于O2和CO在人体内的溶解度[32-33]，类似的处理方式在其他分子模拟研究中也经常使用[31,34-35]。将人血红蛋白置于正方形盒子中心，溶剂化后将气体分子随机分布在蛋白周围的溶液中，添加离子使得模拟体系保持电中性。x、y和z方向均采用周期性边界条件。

1.2 模拟方法

首先使用最陡下降法对系统进行能量优化，然后在正则系综（NVT）下固定蛋白运行模拟1 ns，最后在平衡温度和压力下（NPT 系综）进行多次不同初始构象的分子动力学模拟。分子动力学积分采用跳蛙（leap-frog）积分的方法，积分步长为2 fs，每10 步更新一次邻里列表（neighbor list）。库仑和范德华截断半径均为1.4 nm。采用PME 方法处理静电相互作用，采用SHAKE 算法固定键长。温度和压力分别采用V-rescale 和Parrinello-Rahman 方法来控制。

1.3 马尔科夫状态模型和转移路径理论

为了构建马尔科夫状态模型，对气体在α 链和β 链中的扩散进行多次模拟，累计O2和CO 的模拟时长分别为450 ns 和700 ns。以处理后的分子动力学轨迹作为输入，蛋白碳骨架扭转运动为特征值，利用PyEMMA 2.5.5 建立马尔科夫状态模型。为了捕捉到过程中最慢的时间尺度，采用时滞独立成分分析（TICA）方法将输入数据的高维空间线性降维为描述最慢过程的低维空间。在低维空间中，采用K 均值聚类算法获得不同的微观状态。利用最大似然估计算法，在微观状态下生成过渡矩阵P。用不同滞后时间下的隐含时间尺度和查普曼-科尔莫戈罗（Chapman-Kolmogorow）检验对模型进行验证。在微观状态的弛豫时间尺度中找到跳跃间隙，从而确定不同平衡分布的亚稳状态，再由主成分聚类分析(perron analysis)方法将微观状态划分为相对应的亚稳状态。因为模型本身的复杂性，很难了解不同状态之间的平衡和转换，因此对于不同的亚稳状态采用过渡路径理论（TPT）来计算初始状态到终止状态的转移通量，其满足转移概率矩阵的细致平衡条件（micro-reversibility）和各态遍历性（ergodicity）。

1.4 人血红蛋白结构分析

迴转半径（radius of gyration,Rg）反映了蛋白质分子结构紧密程度，其值越小则蛋白质结构越致密[36]，其计算公式为

均方根偏差（root mean square deviation，RMSD)反映了蛋白质结构与特定结构（这里选择晶体结构）的偏离程度，其值越大则蛋白质结构变化越显著[37]，其计算公式为

均 方 根 涨 落（root mean square fluctuation，RMSF）反映了模拟过程中氨基酸残基的涨落情况，其值越大则蛋白质中该氨基酸涨落越明显，局部结构柔性越强，其计算公式为

2 结果及讨论

2.1 迁移路径分析

图1～图4 分别给出了O2和CO 分子从溶液扩散进入人血红蛋白α 和β 链血红素中心的路径和停留/结合位点。停留/结合位点的定义是，在气体扩散过程中气体分子经常出现（长时间停留）的氨基酸位点。目前，通过实验手段已知的主要有六个Xe位点和一个近端位点（DP）[11]。图1～图4 的结果表明，分子动力学模拟研究可以捕获到全部实验观测的气体结合位点。更为重要的是，通过分子动力学模拟，还捕获到多个瞬时空穴（Ph 位点），这是实验无法观测到的。瞬时空穴是由于氨基酸柔性侧链发生瞬时运动而产生，可以充当O2和CO 分子的临时停留位点[38-39]。这些瞬时空穴主要由亮氨酸、苯丙氨酸及缬氨酸等疏水氨基酸组成。图1～图4中同时将肌红蛋白的Xe 结合位点也标记出来。结果表明，尽管人血红蛋白和肌红蛋白功能类似，但是气体分子停留/结合位点存在显著差异。

图1 O2在血红蛋白α链中的结合位点和六条主要迁移路径(蓝色球：O2位点;箭头连线：O2迁移轨迹;青色带：血红蛋白α链骨架;深棕色分子：血红素;Xe1～Xe6：实验确定Xe结合位点;DP：实验确定近端位点;Ph1～Ph5：瞬时空穴;MbXe1～MbXe4：实验确定肌红蛋白Xe结合位点)Fig.1 O2binding sites and 6 major migration pathways in α chain of human hemoglobin(blue ball：O2 sites;line with arrow：migration trajectories of O2;cyan belt：backbone of α chain of human hemoglobin;brown stick：heme;Xe1—Xe6：Xe binding sites determined by experiments;DP：distal path site determined by experiemnts;Ph1—Ph5：instant cavities;MbXe1—MbXe4：Xe binding sites of myoglobin determined by experiments)

下面将详细阐述O2和CO 在人血红蛋白α 和β链中的迁移特点。

图1 给出了O2在血红蛋白α 链中的结合位点和主要迁移路径。除了实验观测的位点之外，模拟还捕获到5 个瞬时空穴，标记为Ph1～Ph5，这些空穴捕获迁移路径中的O2分子。O2扩散的六条路径及其在分子模拟中出现次数和对应的概率如表1 所示。表1 结果表明，路径1（水溶液→Xe6→Ph3→Xe5→DP）为O2分子在人血红蛋白α 链中迁移扩散的主要通道，而Xe1 和Xe5 结合位点是O2迁移路径中的关键位点。值得一提的是，尽管路径5（水溶液→DP）迁移距离最短，但是其在可能迁移路径中的比例仅为6.0%。这说明O2迁移并不是想象中的最短迁移途径，而是经历复杂的动态迁移网络。

图2 给出了CO 在血红蛋白α 链中的结合位点和主要迁移路径。与O2结果类似，所有实验观测得到的Xe 结合位点和近端位点均可以结合CO。此外，还存在6个瞬时空穴，标记为Ph1～Ph6。CO扩散的六条路径及其在分子模拟中出现的次数和对应的概率如表2 所示。结果表明，路径4（水溶液→DP）和路径5（水溶液→Ph1→Xe3→DP）是CO 在α链中较为重要的两条通道。这说明尽管CO 和O2具有相似的迁移路径，但是这些路径的优先级是不一样的。除此之外，与O2不同，Ph1 和Xe3 对于CO 的停留和扩散影响更为显著，即迁移路径中的关键位点也是不一样的。CO 和O2迁移路径选择概率上的区别也许是CO失活人血红蛋白的重要原因之一。

表1 O2在α链中迁移路径的次数及概率Table 1 Frequency and probability of O2 migration pathways in α chain

图2 CO在血红蛋白α链中的结合位点和六条主要迁移路径(红色球：CO位点;箭头连线：CO迁移轨迹;青色带：血红蛋白α链骨架;深棕色分子：血红素;Xe1～Xe6：实验确定Xe结合位点;DP：实验确定近端位点;Ph1～Ph6：瞬时空穴;MbXe1～MbXe4：实验确定肌红蛋白Xe结合位点)Fig.2 CO binding sites and 6 major migration pathways in α chain of human hemoglobin(red ball：CO sites;line with arrow：migration trajectories of CO;cyan belt：backbone of α chain of human hemoglobin;brown stick：heme;Xe1—Xe6：Xe binding sites determined by experiments;DP：distal path site determined by experiemnts;Ph1—Ph6：instant cavities;MbXe1—MbXe4：Xe binding sites of myoglobin determined by experiment)

表2 CO在α链中迁移路径的次数及概率Table 2 Frequency and probability of CO migration pathways in α chain

类似地，可以分析O2和CO 在血红蛋白β 链中的结合位点和主要迁移路径，以及各个路径所占比例，结果如图3 和图4、表3 和表4 所示。结果表明，β 链中存在8 个O2停留位点以及6 条主要通道，由Xe1 作为入口的通道3 是O2的主要通道。CO 的停留位点有10 个，相比O2来说增加了血红素附近的Ph2 以及Ph3 位点，其迁移通道也存在6 条，主要是直接从水溶液经His 残基和血红素结合。CO 迁移扩散途径较短，尽管其结合常数（热力学参数）要低于O2[40]，但从动力学上会造成蛋白β 链丧失功能。因此CO 造成人血红蛋白丧失功能是一个动力学控制过程。在α 链以及β 链中，O2主要都是由内部空腔组成的通道进行迁移，CO则主要经由路径较短的通道进入。而β 链相对于α 链来说，没有发现经由Xe3 直接进行迁移的通道存在，而这在α 链中却是较为重要的气体迁移通道。需要强调的是，气体分子在α 链和β 链中迁移途径存在差异，表明尽管两个亚基功能相同，但是发挥作用的手段却不同，这可能是生物在进化过程中为了应对环境复杂性而逐渐形成的。

表3 O2在β链中迁移路径的次数及概率Table 3 Frequency and probability of O2 migration pathways in β chain

图3 O2在血红蛋白β链中的结合位点和六条主要迁移路径(蓝色球：O2位点;箭头连线：O2迁移轨迹;粉色带：血红蛋白β链骨架;深棕色分子：血红素;Xe1～Xe6：实验确定Xe结合位点;DP：实验确定近端位点;Ph1：瞬时空穴;MbXe1～MbXe4：实验确定肌红蛋白Xe结合位点)Fig.3 O2binding sites and 6 major migration pathways in β chain of human hemoglobin(blue ball：O2 sites;line with arrow：migration trajectories of O2;pink belt：backbone of β chain of human hemoglobin;brown stick：heme;Xe1—Xe6：Xe binding sites determined by experiments;DP：distal path site determined by experiemnts;Ph1：instant cavities;MbXe1—MbXe4：Xe binding sites of myoglobin determined by experiments)

表4 CO在β链中迁移路径的次数及概率Table 4 Frequency and probability of CO migration pathways in β chain

2.2 血红蛋白结构变化与气体迁移之间的关系

前述分析可知，气体分子迁移进入血红素结合位点时，人血红蛋白由于柔性会产生瞬时空穴，即人血红蛋白结构变化与气体分子迁移存在一定的相关性。因此，在体系达到平衡之后，进一步分析了气体迁移过程中人血红蛋白结构的变化，结果如图5和图6所示（最后1 ns的统计数据）。

图5(a)和(b)分别给出了气体分子迁移过程中Rg和RMSD 变化。结果表明，相对于β 链，气体在α 链中的扩散行为对蛋白结构的影响较大。图5(c)和(d)分别给出了O2和CO 扩散过程中，人血红蛋白α 链和β 链RMSF 变化。这里以0.3 nm 作为截断判据，如果RMSF 值大于0.3 nm，则该氨基酸残基及其周围结构柔性较强，在迁移过程中会发生显著变化。研究结果表明，气体分子扩散过程中人血红蛋白α链的波动明显高于β 链，这与图5(b)所示结果一致。除去蛋白前后端之外，可以看到位于近血红素的氨基酸45，81，92(95)附近均有较大波动，其分别对应于α 链中的Ph1,Xe3,Ph5,Ph6 和β 链中的Xe3,Ph2,Ph3 位点，这些关键位点恰恰位于前述气体迁移中的主要通道，即气体分子迁移过程与人血红蛋白关键氨基酸运动相关联。

图4 CO在血红蛋白β链中的结合位点和六条主要迁移路径(红色球：CO位点;箭头连线：CO迁移轨迹;粉色带：血红蛋白β链骨架;深棕色分子：血红素;Xe1～Xe6：实验确定Xe结合位点;DP：实验确定近端位点;Ph1～Ph3：瞬时空穴;MbXe1～MbXe4：实验确定肌红蛋白Xe结合位点)Fig.4 CO binding sites and 6 major migration pathways in β chain of human hemoglobin(red ball：CO sites;line with arrow：migration trajectories of CO;pink belt：backbone of β chain of human hemoglobin;brown stick：heme;Xe1—Xe6：Xe binding sites determined by experiments;DP：distal path site determined by experiemnts;Ph1—Ph3：instant cavities;MbXe1—MbXe4：Xe binding sites of myoglobin determined by experiments)

图5 四个模型蛋白迴转半径以及RMSD/RMSFFig.5 Gyration and RMSD/RMSF in four protein models

除此之外，α 链的氨基酸17和72附近的蛋白片段显现出较为明显的波动，其分别对应于α 链中的Ph3和MbXe3位点。这表明尽管在晶体结构中无法观测到Ph3 瞬时空穴，但其所在的α 链氨基酸17 是气体分子迁移主通路中的关键位点，这是分子动力学模拟的独特发现。而MbXe3 位点是肌红蛋白的O2主要结合位点，其柔性也较强，但是在人血红蛋白中气体分子却很少在此占据，主要原因在于其并不处于气体分子的扩散路径上。值得一提的是，人血红蛋白β 链的氨基酸17 波动非常小，但其处于Ph1 瞬时空穴，这源于β 链该氨基酸位点本身具备一定容纳O2分子能力[26]，这也充分说明了气体在人血红蛋白中迁移的复杂性。

对比CO 和O2迁移过程中，人血红蛋白α 链和β链的RMSF 可知，CO 对人血红蛋白结构影响要高于O2，即CO 迁移过程中会造成人血红蛋白结构的变化，这可能会进一步影响人血红蛋白的载氧能力。

进而，分析了气体迁移过程中人血红蛋白二级结构的变化，结果如图6所示。

图6 蛋白二级结构变化Fig.6 Protein secondary structure changes

天然人血红蛋白晶体结构由α螺旋和无规卷曲组成，而气体迁移过程中则会造成α 螺旋含量的降低，无规卷曲含量的增加，甚至会出现少量β折叠结构，即迁移过程与人血红蛋白结构变化相关联，这与前面的结论一致。总体而言，相对于β链，迁移过程中气体分子对α 链影响要更大一些。因此，用马尔科夫状态模型（Markov state model, MSM）进一步研究气体迁移过程中人血白蛋白α 链结构变化，结果如图7和图8所示。

首先，利用PyEMMA 2.5.5建立马尔科夫状态模型，通过PCCA 聚类分析得到了不同平衡分布的亚稳状态。在O2扩散的过程中，血红蛋白的结构显现出4 种亚稳态构象，结果如图7 所示。需要说明的是，这里依据不同状态中气体分子和血红素距离的远近为依据，气体和血红素结合最少且距离最远的状态作为初始状态，反之则作为终止状态，对分出的亚稳状态进行路径迁移理论分析，以95%为截断，比例太小的路径忽略不计。

结果表明，由于模拟过程中全部O2分子均进入α 链内部，因此O2在血红素附近的结合状态占据的比例最大，其余状态平衡时占据的比例较小。比较O2_MS1 与晶体结构可知，O2_MS1 状态结构变化较大的区域是对应路径1 的入口氨基酸（黑色虚线圆圈所示），此处α 链由α 螺旋变成了无规卷曲。该结构状态直接转换到O2_MS4 状态的转移概率最大，而该变化正好对应了路径1（水溶液→Xe6→Ph3→Xe5→DP），成 为O2迁 移 的 主 要 路 径。O2_MS1 结构也可转换成O2_MS2 和O2_MS3，分别对应路径6（水溶液→Ph1→Xe3→DP）中的Ph1 瞬时空穴结构和路径4（水溶液→Ph5→Xe5→Xe4→DP）中的Ph5 瞬时空穴结构。相对于晶体结构，结合位点之间的缝隙均变窄，和血红素之间的距离变小（黑色虚线圆圈所示），即O2分子与α 链特定氨基酸结合造成蛋白质结构的变化，从而形成瞬时空穴以结合气体分子。此时，在O2_MS2 状态下路径1 的入口氨基酸结构又由无规卷曲变回α 螺旋结构，说明气体分子在蛋白质内部迁移过程中，蛋白质结构始终处于动态变化状态以实现气体分子的定向输运。相比于CO，O2经由蛋白尾部附近Ph5 位点的路径比例较大，故在O2_MS4 结构中尾部氨基酸波动较剧烈。

类似地，对CO 扩散过程进行分析，结果如图8所示。由结果可知，血红蛋白的结构显现出5 种亚稳态构象。由于模拟过程中CO 扩散到血红素结合位点的数量相对O2较少（即血红素与CO 结合要低于与O2的结合常数，这与文献报道一致[40]），因此中间状态CO_MS2 所占据的比例最大。比较CO_MS1与晶体结构可知，CO_MS1 状态结构变化较大的区域是对应比例最大的路径5 以及路径1，3，6 的入口氨基酸（黑色虚线圆圈所示）。和O2不同，由于CO分子较多，在未结合状态时多条通道被CO 分子占据，故CO_MS1 结构变化较剧烈。在CO_MS1 结构中，血红素分子上方空间变窄，该变化正好对应了路径5（水溶液→Ph1→Xe3→DP），成为CO 迁移的主要路径。同时Ph5 位点也由α 螺旋变成了无规卷曲，对应于路径3（水溶液→Ph5→Ph6→DP）。Ph4和Xe6 处的α 螺旋和邻近螺旋距离变小，其对应于路径6（水溶液→Ph4→Xe2→Xe4→DP）和路径1（水溶液→Xe6→Ph3→Xe4→DP）。CO_MS1 结构也可转换成CO_MS2，在CO_MS2 状态处Ph4 位点由α 螺旋变成了无规卷曲（黑色虚线圆圈所示），其对应于路径6结构变化的中间状态。同时位于背后的两条α 螺旋形成的空间变小，这在一定程度上可以阻止溶液中的CO 分子进入蛋白，起到关闭通路的作用。而CO_MS3 在CO_MS1 转化到CO_MS4 的过程中起到连接作用，其结构变化基于两者结构之间。和O2分子不同的是，在这个过程出现了独立的CO_MS0状态。CO_MS0 状态的结构变化是对应于路径最短的通道4，由图8中CO_MS0结构变化可知，在通道4周围的蛋白结构变为了无规卷曲，这在之前的各种状态中没有出现过，且蛋白空腔变小，柔性增大，均有利于调节气体分子进入蛋白。CO 可以通过该通道直接进入血红素上方。这个通道对应的蛋白质构象占据很大的比例，因此形成一个独立的状态。由前面的结论可知，CO 相对O2来说更倾向于直接由路程最短的通道和血红素直接结合，因此形成了其特有的CO_MS0 状态结构，进一步说明了不同种类的气体分子直接影响人血红蛋白的结构，进而影响其功能。

图7 O2迁移过程中人血红蛋白α链结构转换路径(灰色带状模型：人血红蛋白α链晶体结构;O2_MS1(红色)：初始状态结构,即O2未结合态;O2_MS2(黄色)和O2_MS3(粉色)：中间状态结构,即迁移路径关键结构;O2_MS4(橙色)：最终结合状态结构;含箭头连线：不同状态迁移路径,其粗细表明迁移概率;不同颜色圆：不同状态,其大小表明状态出现概率;不同颜色带状模型：对应状态人血红蛋白α链结构;绿色球状模型：血红素分子)Fig.7 Tructural transition pathways of α chain of human hemoglobin during O2 migration processes(gray ribbon model：crystal structure of α chain of human hemoglobin;O2_MS1 in red：initial structure,i.e.,O2 unbinding state;O2_MS2 in yellow and O2_MS3 in pink：middle states,i.e.,key structures during O2 migration;O2_MS4 in orange：final O2 binding state;lines with arrow：transition pathways among different states,their thickness reflects transition probability;circles with different color：different states,their sizes reflect appearance probability;ribbon model with different color：corresponding structures of α chain of human hemoglobin;green ball model：heme molecule)

图8 CO迁移过程中人血红蛋白α链结构转换路径(灰色带状模型：人血红蛋白α链晶体结构;CO_MS1(红色)：初始状态结构,即CO未结合态;CO_MS2(黄色)和CO_MS3(粉色)：中间状态结构,即迁移路径关键结构;CO_MS0(紫色)：独立状态结构,即不参加转换的结构;CO_MS4(橙色)：最终结合状态结构;含箭头连线：不同状态迁移路径,其粗细表明迁移概率;不同颜色圆：不同状态,其大小表明状态出现概率;不同颜色带状模型：对应状态人血红蛋白α链结构;绿色球状模型：血红素分子)Fig.8 Structural transition pathways of α chain of human hemoglobin during CO migration processes(gray ribbon model：crystal structure of α chain of human hemoglobin;CO_MS1 in red：initial structure,i.e.,CO unbinding state;CO_MS2 in yellow and CO_MS3 in pink：middle states,i.e.,key structures during CO migration;CO_MS0 in purple：independent state,i.e.,a structure not in transition network;CO_MS4 in orange：final CO binding state;lines with arrow：transition pathways among different states,their thickness reflects transition probability;circles with different color：different states,their sizes reflect appearance probability;ribbon model with different color：corresponding structures of α chain of human hemoglobin;green ball model：heme molecule)

综上可知，CO 在人血红蛋白迁移过程中，不仅仅选取与O2相同的迁移途径，而且也有独特的“便捷通道”，同时造成人血红蛋白结构的更大变化，这就是CO 造成人血红蛋白载氧能力丧失的重要原因。

3 结 论

本文采用分子动力学模拟结合马尔科夫状态模型的方法对O2和CO 从溶剂外部进入人血红蛋白内部结合位点的迁移过程进行了研究，分别得出了两种气体在人血红蛋白的α 链以及β 链中的主要结合/停留位点和主要迁移通道，以及迁移过程中蛋白结构的变化。人血红蛋白α 链中存在12 个O2结合位点以及6 条主要迁移通道，其中由Xe6 作为大门的通道1 是主要的进入通道。而CO 的结合位点有13 个，相比O2来说多了一个Ph6 位点，其主要的迁移通道也有6个，但是路径较短的通道4和5是其主要的进入通道。β 链中存在8 个O2停留位点以及6条主要通道，由Xe5 作为入口的通道3 是O2的主要通道。CO 的停留位点有10 个，相比O2来说增加了血红素附近的Ph2 以及Ph3 位点，其迁移通道也存在6 条，主要由路径较短的通道6 进入蛋白内部。在α 链以及β 链中，O2主要都是由内部空腔组成的通道进行迁移，CO 则主要经由路径较短的通道进入。而β 链相对于α 链来说，没有发现由Xe3 直接进行迁移的通道存在。除了蛋白前后端，四种蛋白模型在血红素附近的氨基酸都波动较大，主要是受气体迁移的影响，特别的对于α 链还存在两个波动大的蛋白区域。由马尔科夫状态模型分析得出蛋白结构变化的详细转移图，可知为了实现气体分子的定向输运，蛋白质结构始终处于动态的变化状态。本文尝试对R 态下O2和CO 进入人血红蛋白迁移过程的模拟研究，对血液替代品的分子设计有重要的基础研究价值，在原子水平上加深对气体分子在人血红蛋白中迁移动力学过程的理解。在未来工作中，可以更进一步指认各个关键蛋白区域调控气体分子迁移的关键氨基酸及其结构，通过饱和突变和理性设计等方式筛选和设计载氧能力更强的人工血液替代品。