基于拉伸分子动力学模拟的黄嘌呤氧化酶抑制剂解离途径的理论研究

齐人睿，李明昊，常 浩，付学奇，高 波，韩葳葳，韩 璐，李婉南

（1.吉林大学分子酶学工程教育部重点实验室,生命科学学院,长春130012；2.吉林省特医食品生物科技有限公司,长春130052）

痛风的临床症状有高尿酸血症、肾尿酸结石、慢性肾炎以及关节畸形和慢性关节炎，给患者造成极大的痛苦［1，2］.痛风是由嘌呤代谢紊乱导致尿酸沉积在组织中而引发的慢性疾病［3］.嘌呤代谢过程中产生的黄嘌呤和次黄嘌呤在黄嘌呤氧化酶（XO）的氧化作用下生成尿酸，是导致尿酸在体内堆积的主要原因之一.目前，黄嘌呤氧化酶的抑制剂是用于治疗痛风及高尿酸血症的主要药物，因此黄嘌呤氧化酶与抑制剂的相互作用是一个研究热点.

XO是一种复合型的核黄素蛋白，由两个完全重复对称的结构单元组成，每个结构单元由3条异源的肽链组成［图1（A）］，其催化中心包括一个黄素腺嘌呤二核苷酸（Flavin adenine dinucleotide，FAD）、两个铁-硫中心（Iron-sulfur centers，Fe/S I，II）以及一个钼蝶呤中心（Molybdopterin，Mo-pt）.图1（B）给出了黄素腺嘌呤二核苷酸（C27H33N9O15P2，FAD）、MTE（Phosphonic acidmono ester，C10H14N5O6PS2）、鸟嘌呤（Guanine，C5H5N5O）、MOS（Dioxothiomolybdenum ion，HMoO2S）和铁硫簇（Fe2S2cluster，FES）5种小分子配体的三维（3D）结构，其中鸟嘌呤，MOS及MTE结合在由氨基酸残基Phe767，Phe798，Gly800，Glu802，Gly876，Arg880，Ala910，Phe911，Phe914，Thr1077，Ala1078，Ala1079，Ser1080，Val1081及Ser1082组成的钼蝶呤中心口袋［图1（C）］，FES结合在由氨基酸残基Lys40，Leu41，Gly42，Cys43，Gly44，Gly46，Ala50，Cys51，Asn71，Cys73，Ser111，Gln112，Cys113，Phe115，Cys116，Cys143，Cys148，Cys150及Thr151组成的铁-硫中心口袋［图1（D）］，FAD结合在由残基Asn251，Lys256，Leu257，Val258，Val259，Gly260，Asn261，Thr262，Val342，Ser347，Asp360，Leu398，Glu402，Ile403，Leu404及Arg426组成的黄素腺嘌呤二核苷酸中心口袋［图1（E）］.而钼蝶呤中心是黄嘌呤氧化酶催化底物产生尿酸的关键位点［4，5］，在黄嘌呤的氧化过程中，Mo（Ⅵ）被还原为Mo（Ⅳ），产生的电子经由FAD和Fe/S传递给氧分子，与氢离子结合形成过氧化氢，Mo（Ⅳ）则再氧化为Mo（Ⅵ）［6］.

Fig.1 Overview of xanthine oxidase(XO,PDB ID:3NVW)

别嘌呤醇（Allopurinol）的结构与黄嘌呤相似，可以竞争性抑制黄嘌呤氧化酶［7］，使黄嘌呤氧化酶不能催化次黄嘌呤和黄嘌呤生成尿酸，从而降低了尿酸在体内的沉积，是治疗高尿酸血症的主要药物，但其具有恶心、头痛、白细胞减少或增多等一系列副作用［8～10］.有研究表明，属于异黄酮类化合物的葛根素（Puerarin）对黄嘌呤氧化酶表现出较好的抑制效果，属于温和性的竞争型黄嘌呤氧化酶抑制剂，能够明显降低高尿酸血症患者血清中的尿酸水平［11～13］，有望被开发为新的治疗痛风的药物，但是关于抑制剂与黄嘌呤氧化酶解离途径之间的相互作用却鲜有报道.本文通过对XO/allopurinol和XO/puerarin 2种复合物进行拉伸分子动力学模拟，从计算生物学的角度解释黄嘌呤氧化酶与抑制剂的结合方式，以及抑制剂从黄嘌呤氧化酶上解离的过程.

1 计算模拟方法

1.1 模拟体系准备

蛋白质的初始结构为XO（PDB ID：3NVW）［14］.从Chemspider数据库中下载Allopurinol和Puerarin的几何结构，底物初始结构采用Gaussian 09软件包，在B3LYP 6-31+G*基组水平上进行优化［15］.优化后结构的分子轨道使用Multiwfn程序［16］显示.采用AutoDock Vina［17，18］对接构建XO/allopurinol和XO/puerarin 2种复合物的结构.

1.2 XO蛋白通道预测

采用CAVER 3.0软件，分别对2种复合物的10 ns常规分子动力学模拟轨迹产生的1000个快照进行配体解离通道的预测［19］.设置每个构象中识别出的通道瓶颈半径（Bottleneck radiu）≥0.12 nm，聚类阈值（Clustering_threshold）为7.5.由于在同一构象中可能会识别出多条几乎同轴线的通道，为了去除多余通道，则保留每个构象中能量最低通道（One_tunnel_in_snapshot cheapest）.

1.3 拉伸分子动力学模拟

依据CAVER 3.0得到的结果，确定配体解离通道及其拉伸方向，蛋白采用GROMACS软件中AMBER99SB-ILDN力场，小分子力场参数由GAFF程序生成，进而实施拉伸分子动力学模拟［20，21］.在进行拉伸分子动力学模拟之前，将分子对接得到的2个复合物分别放入到两个周期性边界为1 nm的立方体盒子中，盒子中除蛋白外的位置采用TIP3P［22］水模型填充，其中在含有别嘌呤醇的体系中添加57349个水分子，在含葛根素的体系中添加57340个水分子.为了保证电荷的平衡，向2个模拟体系中均添加5个Cl-.在进行拉伸分子动力学模拟时，首先使用最陡下降法对模拟体系进行50000步的能量最小化，然后再进行100 ps的恒温恒容（NVT）模拟和100 ps的恒温恒压（NPT）模拟.在进行NVT模拟时，2个体系均采用LINCS算法［23］约束所有共价键的振动.采用Berendsen算法［24］控温，将体系的温度保持在300 K，温度耦合常数为0.1 ps.采用SETTLE算法限制水分子的几何构象，范德华相互作用的双程截断半径设为1.0和1.2 nm，长程静电相互作用使用Particle mesh Ewald（PME）［25］方法处理.在进行NPT模拟时，采用Berendsen算法［24］控压，压力耦合的时间常数设为2.0 ps，压力耦合的参考压力设为1×105Pa，原子的各种参数限制及温度设置均与NVT模拟相同，NPT模拟之后体系的密度达到最优状态.限制之后对2个体系进行步长为1 fs、牵引力常数为1000 kJ·mol-1·nm-2、时长为4000 ps的恒速拉伸分子动力学模拟（0.001 nm/ps）.

1.4 平均力势计算

为了研究别嘌呤醇、葛根素与黄嘌呤氧化酶之间的结合能，需要计算一系列伞形采样模拟的平均力势（PMF）［26，27］.使用伞形偏离势将黄嘌呤氧化酶的质心限制在窗口中，以质心与抑制剂之间的距离作为反应坐标，沿反应轴，每隔约0.1～0.2 nm设置一个伞形采样窗口，其中XO/allopurinol体系取9个窗口，XO/puerarin体系取7个窗口，每个窗口独立模拟，首先进行100 ps的NPT模拟，再进行10 ns的采样，最终通过使用GROMACS软件的gmx wham工具获得PMF曲线.

2 结果与讨论

2.1 抑制剂结构的优化及与黄嘌呤氧化酶的相互作用

为了得到模拟复合物中合理的抑制剂结构，对别嘌呤醇和葛根素2种抑制剂进行结构优化，研究了它们的化学性质.最高占据分子轨道（HOMO）与最低未占据分子轨道（LUMO）之间的能量差称为HOMO-LUMO gap，可用于预测电子转移的强度和稳定性.其中别嘌呤醇的化学结构及其LUMO轨道分别见图2（A）和（B）；葛根素的化学结构及其LUMO轨道分别见图2（C）和（D）.2种抑制剂LUMO轨道能量贡献主要位于2个芳香环上，结构与黄嘌呤相似，这可能是2种化合物能竞争性抑制XO的原因.

Fig.2 Chemical structures(A,C)and population analysis of LUMO orbits(B,D)for Allopurinol(A,B)and Puerarin(C,D)

采用AutoDock Vina软件，将别嘌呤醇和葛根素分别对接到黄嘌呤氧化酶中.对接结果显示，XO的Glu802，Leu873，Arg880，Phe914，Phe1009，Thr1010，Val1011，Ala1078，Ala1079及Glu1261氨基酸残基将别嘌呤醇包裹在中心［图3（A）］，XO的Leu648，Phe649，Lys771，Glu802，Leu873，His875，Ser876，Arg880，Ser1008，Phe1009，Thr1010，Val1011，Phe1013，Leu1014，Pro1076及Ala1078等氨基酸残基将葛根素包裹在中心［图3（B）］.其结果与图1（C）中鸟嘌呤的结合位置极其相似，由此可见，两者均成功地对接到了XO的MO-pt中心.此外，别嘌呤醇和葛根素与活性中心氨基酸形成的氢键可以帮助其更好地结合在XO的钼蝶呤中心.

Fig.3 Allopurinol(A)and Puerarin(B)docking to MO-pt pocket of XO

丙氨酸扫描是指将活性中心的氨基酸突变成丙氨酸，用以分析单个氨基酸残基对蛋白质与配体结合的影响.XO的丙氨酸扫描的结果显示，来自于MO-pt中心的Val789，Arg880，Phe911，Phe914和Val1081突变成丙氨酸之后，XO与别嘌呤醇之间的结合自由能分别下降了-11.13，-9.67，-16.50，-13.45和-7.60 kJ/mol［图4（A）］，与葛根素的结合自由能分别下降了-11.97，-10.50，-14.41，-15.12和-8.44 kJ/mol［图4（D）］，这种结合自由能下降较为明显的结果说明以上氨基酸在XO与抑制剂的结合中起到非常重要的作用.此外，来自于Fe/S结合中心［图4（B）和（E）］和FAD结合中心［图4（C）和（F）］氨基酸的突变也能在一定程度上影响XO与抑制剂的结合.

Fig.4 Computational scanning of the binding site residues in the XO/allopurinol(A―C)and XO/puearin(D―F)complexes

2.2 XO蛋白通道的分析

在进行拉伸分子动力学模拟之前，利用CAVER 3.0软件预测抑制剂从XO蛋白上离去的通道，图5显示了通道的形状及其在XO中的相对位置.抑制剂从XO的离去通道位于MO-pt中心，主要由Leu648，Phe649，Lys771，Ser876，Val1011和Phe1013组成，这些氨基酸在空间结构上形成一个较为规则的近似圆形的通道，且XO通道的形状也较为规则.蛋白通道的瓶颈半径（Bottleneck radius）、长度（Length）和弯曲率（Curvature）是用于定性分析蛋白质通道的指标，瓶颈半径越大，表示蛋白通道的空间维度越大；长度越短，表示配体从蛋白上解离需要经过的路程越短；弯曲率越小，则配体从蛋白解离的通道形状越平坦.

Fig.5 CAVER 3.0 analysis of the shape and size of tunnel for XO

依据XO/allopurinol复合物的分子动力学模拟轨迹，得到XO蛋白离去通道的平均瓶颈半径（Average bottleneck radius）、平均长度（Average length）和平均弯曲率（Average curvature）分别为0.26，1.24和0.12 nm；依据XO/puerarin复合物的分子动力学模拟轨迹，得到XO蛋白离去通道的平均瓶颈半径、平均长度和平均弯曲率分别为0.25，0.58和0.12 nm.由此可见，抑制剂从XO蛋白上解离时会经过一个相对较为顺畅的通道，接下来的拉伸分子动力学模拟即基于这条通道进行.从计算结果可见，两个蛋白通路的瓶颈半径和弯曲率相差无几，而通路的长度相差较大，这可能导致相比于葛根素，别嘌呤醇更不易从XO蛋白中解离.

2.3 别嘌呤醇和葛根素从XO蛋白上解离的过程

为了研究别嘌呤醇和葛根素从XO蛋白上解离的过程，分别对XO/allopurinol和XO/puerarin进行了6次时长为4000 ps的拉伸分子动力学模拟，表1显示了抑制剂从XO解离所需要的拉力和时间，发现别嘌呤醇从XO中解离所需要的拉力值明显高于葛根素.同样的，别嘌呤醇从XO中解离所需要的时间也明显高于葛根素.这说明相比葛根素，别嘌呤醇在解离过程中受到的阻碍更大，更不容易从XO中解离，即XO与别嘌呤醇之间的结合更加紧密，也可能是别嘌呤醇对XO具有较强抑制效果的原因.

Table 1 Time and largest applied force of ligand dissociating from XO during the SMD simulations

Fig.6 Steered molecular dynamics(SMD)simulation of typical force profiles of allopurinol and purearin pulled out of the MO-pt pocket of XO along the unbinding pathway

图6 给出了别嘌呤醇在解离过程中所需外力随模拟时间的变化.为验证该曲线具有代表性，共进行了6次重复性实验，结果见图S1（A）（见本文支持信息）.在0～700 ps时间内，别嘌呤醇所受的外力沿直线上升至700.54 pN.外力达到最大值之后逐渐下降，直至920 ps时，外力出现小幅度的上升.在2526 ps时到达平衡，所受外力在0 pN处波动，说明此时别嘌呤醇已经成功从XO上解离.图6还给出了葛根素在解离过程中所需外力随模拟时间的变化.为验证该曲线具有代表性，同样进行了6次重复性实验，结果见图S1（B）.在0～425 ps时间内，葛根素所受的外力沿直线上升到428 pN.外力达到最大值之后逐渐下降，直至706 ps时，所受外力产生小幅度的上升.在1089 ps时达到平衡，所受外力在0 pN处波动，说明此时葛根素已经成功从XO上解离.

为了分析在拉伸过程中别嘌呤醇与XO成键即相互作用的变化，选取了700和920 ps时的构象进行分析.在700 ps时，外力达到峰值，别嘌呤醇与Arg880和Thr1010之间存在氢键作用；与Ala1079之间存在碳氢键作用和疏水作用；与S原子产生的π-Sulfur作用对蛋白质稳定起重要作用；与Ser1008，Phe1009，Val1011，Leu1014和Glu1261之间存在范德华作用［图7（A）］.在920 ps时，别嘌呤醇与Phe649和Ser876之间存在范德华作用；与His875之间新产生的π-π堆积类型的疏水作用，可能是导致外力上升的原因［图7（B）］.外力到达峰值后，随着模拟时间的延长，XO与别嘌呤醇形成的非键相互作用减少，所受外力减小，说明非键相互作用是阻碍别嘌呤醇从XO上解离的主要因素.因此这些残基在拉伸过程中对维持XO与别嘌呤醇之间的相对稳定起到重要作用.

Fig.7 Non-bond interactions between ligands and residues

为了分析在拉伸过程中葛根素与XO成键即相互作用的变化，分别选取了425和706 ps时的构象进行分析.在425 ps时，外力达到峰值，葛根素与Ser876和Ala1079之间存在氢键作用；与Leu648，Leu873，Phe914，Phe1009，Val1011，Leu1014和Ala1078之间存在疏水作用，其中与Leu648之间的π-σ相互作用是氢与sp2杂化原子组成平面环之间的弱相互作用；与Phe649，Lys771，His875，Arg880，Ala910，Thr1010，Phe1013和Pro1076之间存在范德华作用［图7（C）］.在706 ps时，葛根素与Leu648，Phe649，Phe1013和Leu1014之间存在范德华作用；与Val1011之间依旧存在疏水作用，与His875之间新产生的碳氢键即一种弱的氢键作用，可能是导致外力小幅度上升的原因［图7（D）］.外力到达峰值后，随着模拟时间的延长，XO与葛根素之间的非键相互作用减少，所受外力减小，显然，这些残基也在拉伸过程中对维持XO与葛根素之间的相对稳定也起到重要作用.

通过分析别嘌呤醇、葛根素从XO解离路径离去过程中所受外力随时间变化的曲线可见，相比于葛根素，别嘌呤醇需更大的外力与更长的时间才能从XO中解离，这可能是由别嘌呤醇与XO的活性中心残基存在更强的相互作用所致；从图7（A）与（C）可见，拉力最大时，Phe1009，Arg880，Ala1079，Val1011和Thr1010均对2种复合物的结构稳定起到重要作用.随着模拟时间增加，由图7（B）与（D）可见，His875在两者从XO的解离过程中均起到阻碍作用，由于别嘌呤醇与His875之间新产生的π-π堆积类型的疏水作用要强于葛根素与His875之间产生的碳氢键（即一种弱的氢键）作用，所以，与葛根素相比，别嘌呤醇更不易从XO中离去.

拉伸分子动力学模拟和CAVER的分析结果表明，Phe649和Phe1013在抑制剂离去以及组成通道的过程中起到非常重要的作用.由图8（A）可见，Phe649与Phe1013位于XO通道的两边，犹如两扇门控制着通道的大小.因此，分析了Phe649与Phe1013苯环之间的距离在6次拉伸模拟过程中随时间的变化.由图8（B）和（C）可见，在别嘌呤醇存在时，Phe649与Phe1013苯环之间的距离收敛为1.25 nm，小于葛根素存在时的1.30 nm［图8（C）］.虽然图8（C）中蓝色和粉色曲线分别在2000和3000 ps开始从1.30 nm处下降，但此时葛根素已从XO中离去（见表1），所以可忽略不计.这种距离上差异也说明别嘌呤醇在解离时遇到的阻碍比较大，不容易从XO蛋白中解离.

2.4 抑制剂解离过程中PMF计算

为了进一步了解别嘌呤醇、葛根素与黄嘌呤氧化酶解离途径之间的相互作用，通过计算PMF，得到了抑制剂与黄嘌呤氧化酶之间的结合能.图9（A）和（B）分别为别嘌呤醇与葛根素从黄嘌呤氧化酶中解离过程的PMF曲线.起初能量值迅速增加，这是由抑制剂与活性位点残基之间相互作用导致的，说明2个抑制剂初始位置与XO的结合状态非常稳定，且相比于葛根素，别嘌呤醇与XO结合更稳定.从计算结果分析，别嘌呤醇从XO中解离的能垒为41.46 kJ/mol，葛根素从XO中解离的能垒为21.70 kJ/mol，说明别嘌呤醇与XO结合更加紧密，这也是别嘌呤醇对XO具有更强抑制效果的原因.

Fig.9 PMF profile along the unbinding pathway of XO/allopurinol(A)and XO/puerarin(B)

综上所述，通过对XO/allopurinol与XO/puerarin 2个复合物的拉伸分子动力学模拟和PMF的计算，从蛋白质构象和能量差异两个角度研究了抑制剂别嘌呤醇和葛根素分别从XO蛋白上解离的过程.通过分子对接结果显示，两者均处于XO的活性中心位置；丙氨酸扫描的结果显示，Val789，Arg880，Phe911，Phe914和Val1081在XO与别嘌呤醇和葛根素的结合中均起到非常重要的作用；拉伸分子动力学模拟结果显示，随着模拟时间的增加，所受外力值的变化趋势基本一致，且发现在拉伸过程中Arg880，Phe1009，Thr1010，Val1011和Ala1079对维持两种复合物的结构稳定均起到重要作用，Phe649与Phe1013在其解离的过程中起到门控的作用，His875对2个抑制剂均起到阻碍解离的作用，但相比于葛根素，别嘌呤醇则需更大的外力和更长的时间才能从XO的MO-pt中心解离.研究结果阐释了别嘌呤醇、葛根素与黄嘌呤氧化酶解离途径之间有着非常相似的相互作用，在解离过程中，两者均受到许多相同残基的阻碍作用，但发现这些残基对别嘌呤醇的阻碍作用更大，这也体现了葛根素会较易从XO中解离，该结果也为研究抑制剂与XO的相互作用提供了一定的理论指导.

支持信息见http：//www.cjcu.jlu.edu.cn/CN/10.7503/cjcu20200752.