基于分子动力学模拟和自由能计算的非小细胞肺癌对克唑替尼(crizotinib)耐药性研究

桑鹏 李智 杨力权 宋玲

引用本文:
Citation:

基于分子动力学模拟和自由能计算的非小细胞肺癌对克唑替尼(crizotinib)耐药性研究

    作者简介: 桑 鹏(1983– ),男,山东人,博士,讲师,主要从事结构生物信息学方面的研究. E-mail:speng431@163.com;
    通讯作者: 宋玲, songl0870@yeah.net

Insight into crizotinib resistance of non-small cell lung cancer derived from molecular dynamics simulations and free energy analysis

    Corresponding author: SONG Ling, songl0870@yeah.net
  • 摘要: 为了研究间变性淋巴瘤激酶(anaplastic lymphoma kinase,ALK)F1174V突变引起非小细胞肺癌患者对crizotinib的耐药性机制,使用分子动力学模拟(molecular dynamics simulations, MD). 本质动力学(essential dynamics, ED)分析及分子力学泊松-玻尔兹曼比表面积(molecular mechanics Poisson-Boltzmann surface area, MM-PBSA)结合自由能计算等方法研究了野生型及F1174V突变型ALK与克唑替尼(crizotinib)结合模式的差异. MD模拟和ED分析结果表明,F1174V突变导致ALK的crizotinib结合口袋部位构象柔性降低,由此可能造成药物进出ALK速度减慢,影响药物发挥正常作用. MM-PBSA结合自由能计算结果表明,F1174V突变造成ALK与crizotinib的结合能力降低. 进一步构建了野生型和突变型ALK的自由能图谱(free energy landscape,FEL),并选取了二者最主要自由能井的代表性结构进行比较. 结果表明,F1174V突变可以显著改变ALK的构象,特别是可以通过收缩P-环和缩小结合口袋来减弱其与crizotinib的相互作用. 研究有助于揭示基于ALK的靶向药物耐药机制,并为开发抗非小细胞肺癌药物提供帮助.
  • 图 1  模拟过程中ALK野生型和突变型骨架原子相对初始结构的RMSD值.

    Figure 1.  The RMSD Values of backbone atoms from the initial structures of WT and mutant ALK during simulation

    图 2  野生型和突变体残基在模拟过程中的RMSF值

    Figure 2.  The RMSF Value of eacd residue of WT and mutant during the simulation

    图 3  野生型和突变体前20个本征向量对应的本征值

    Figure 3.  Eigenvalues as a function of the first 20 eigenvector of WT and mutant

    图 4  野生型和突变体ALK的FEL及自由能最小化区域的代表性结构.自由能图谱通过主成分分析方法进行构建,其中野生型和突变体的分子结构用 Pymoe软件[28]作图,ALK和crizontinib分子之间的相互作用由Ligplot软件[29]进行分析.

    Figure 4.  FEL and the representative structure in the lowest energy basin. The FEL were constructed based on the principal component analysis. The surface and cartoon representation of crizotinib were generated by Pymol[28]. Residue interations in the ALK-crizotinib interface was analyzed by Ligplot[29].

    表 1  组成crizotinib结合部位的残基RMSF值

    Table 1.  RMSF values of the crizotinib binding sites

    残基均方根波动/nm
    野生型突变型
    11220.23460.2108
    11300.24420.0967
    11480.26150.1123
    11970.20960.0911
    11980.25030.2442
    11990.18100.1784
    12560.22190.2087
    下载: 导出CSV

    表 2  野生型和突变体ALK在模拟过程中的几何属性

    Table 2.  Structural properties of the WT and mutant of ALK

    蛋白质NNCNHBRg/nmSASA/nm2
    野生型334530 (1094)227 (7.6)1.96 (0.02) 140.82 (2.82)
    突变型340622 (1125)232 (7.5)1.94 (0.01)135.26(2.75)
    下载: 导出CSV

    表 3  野生型和突变体ALK与crizotinib之间的结合自由能

    Table 3.  Binding free energy and its components of the WT and mutant KJ/mol

    结合自由能组分野生型突变型
    ΔEvdw−41.52−39.56
    ΔEele−242.37−211.05
    ΔGpolar173.42165.31
    ΔGnonpolar−8.76−4.35
    ΔGMM−283.89−250.61
    ΔGsol164.66160.96
    ΔG−119.23−89.65
    TΔS24.8519.36
    ΔGbind−94.38−70.29
    下载: 导出CSV
  • [1] Lynch T J, Bell D W, Sordella R, et al. Activating mutations in the epidermal growth factor receptor underlying responsiveness of non–small-cell lung cancer to gefitinib[J]. New England Journal of Medicine, 2004, 350(21): 2 129-2 139. DOI:  10.1056/NEJMoa040938.
    [2] Wang D D, Zhou W, Yan H, et al. Personalized prediction of EGFR mutation-induced drug resistance in lung cancer[J]. Scientific Reports, 2013, 3: 2 855. DOI:  10.1038/srep02855.
    [3] Kumar A, Ramanathan K. Analyzing resistance pattern of non-small cell lung cancer to crizotinib using molecular dynamic approaches[J]. Indian J Biochem Biophys, 2015, 52(1): 23-28.
    [4] Schlessinger J. Cell signaling by receptor tyrosine kinases[J]. Cell, 2000, 103(2): 211-225. DOI:  10.1016/S0092-8674(00)00114-8.
    [5] Ni Z, Zhang T C. Computationally unraveling how ceritinib overcomes drug-resistance mutations in ALK-rearranged lung cancer[J]. Journal of Molecular Modeling, 2015, 21(7): 1-10.
    [6] Lee C C, Jia Y, Li N, et al. Crystal structure of the ALK (anaplastic lymphoma kinase) catalytic domain[J]. Biochemical Journal, 2010, 430(3): 425-437. DOI:  10.1042/BJ20100609.
    [7] Katayama R, Shaw A T, Khan T M, et al. Mechanisms of acquired crizotinib resistance in ALK-rearranged lung cancers[J]. Science Translational Medicine, 2012, 4(120): 120ra17.
    [8] Rodig S J, Shapiro G I. Crizotinib, a small-molecule dual inhibitor of the c-Met and ALK receptor tyrosine kinases[J]. Current Opinion in Investigational Drugs, 2010, 11(12): 1 477-1 490.
    [9] Thomas R K. Overcoming drug resistance in ALK-rearranged lung cancer[J]. New England Journal of Medicine, 2014, 370(13): 1 250-1 251. DOI:  10.1056/NEJMe1316173.
    [10] Ou S HI, Azada M, Hsiang D J, et al. Next-generation sequencing reveals a Novel NSCLC ALK F1174V mutation and confirms ALK G1202R mutation confers high-level resistance to alectinib (CH5424802/RO5424802) in ALK-rearranged NSCLC patients who progressed on crizotinib[J]. Journal of Thoracic Oncology, 2014, 9(4): 549-553. DOI:  10.1097/JTO.0000000000000094.
    [11] Friboulet L, Li N, Katayama R, et al. The ALK inhibitor ceritinib overcomes crizotinib resistance in non–small cell lung cancer[J]. Cancer Discovery, 2014, 4(6): 662-673. DOI:  10.1158/2159-8290.CD-13-0846.
    [12] Šali A, Blundell T L. Comparative protein modelling by satisfaction of spatial restraints[J]. Journal of Molecular Biology, 1993, 234(3): 779-815. DOI:  10.1006/jmbi.1993.1626.
    [13] Pronk S, Páll S, Schulz R, et al. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit[J]. Bioinformatics, 2013, 29(7): 845-854.
    [14] Duan Y, Wu C, Chowdhury S, et al. A point‐charge force field for molecular mechanics simulations of proteins based on condensed‐phase quantum mechanical calculations[J]. Journal of Computational Chemistry, 2003, 24(16): 1 999-2 012. DOI:  10.1002/jcc.10349.
    [15] Hornak V, Abel R, Okur A, et al. Comparison of multiple Amber force fields and development of improved protein backbone parameters[J]. Proteins: Structure, Function, and Bioinformatics, 2006, 65(3): 712-725. DOI:  10.1002/prot.21123.
    [16] SchuÈttelkopf A W, van Aalten D M. PRODRG: a tool for high-throughput crystallography of protein–ligand complexes[J]. Acta Crystallographica Section D: Biological Crystallography, 2004, 60(8): 1 355-1 363. DOI:  10.1107/S0907444904011679.
    [17] Jorgensen W L, Chandrasekhar J, Madura J D, et al. Comparison of simple potential functions for simulating liquid water[J]. The Journal of Chemical Physics, 1983, 79(2): 926-935. DOI:  10.1063/1.445869.
    [18] Essmann U, Perera L, Berkowitz M L, et al. A smooth particle mesh Ewald method[J]. The Journal of Chemical Physics, 1995, 103(19): 8 577-8 593. DOI:  10.1063/1.470117.
    [19] Hess B, Bekker H, Berendsen H J, et al. LINCS: a linear constraint solver for molecular simulations[J]. Journal of Computational Chemistry, 1997, 18(12): 1 463-1 472. DOI: 10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H.
    [20] Amadei A, Linssen A, Berendsen H J. Essential dynamics of proteins[J]. Proteins: Structure, Function, and Bioinformatics, 1993, 17(4): 412-425. DOI:  10.1002/prot.340170408.
    [21] Jin H, Zhou Z, Wang D, et al. Molecular dynamics simulations of acylpeptide hydrolase bound to chlorpyrifosmethyl oxon and dichlorvos[J]. International Journal of Molecular Sciences, 2015, 16(3): 6 217-6 234.
    [22] Gohlke H, Kiel C, Case D A. Insights into protein–protein binding by binding free energy calculation and free energy decomposition for the Ras–Raf and Ras–RalGDS complexes[J]. Journal of Molecular Biology, 2003, 330(4): 891-913. DOI:  10.1016/S0022-2836(03)00610-7.
    [23] Kollman P A, Massova I, Reyes C, et al. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models[J]. Accounts of Chemical Research, 2000, 33(12): 889-897. DOI:  10.1021/ar000033j.
    [24] Xu Y, Wang R. A computational analysis of the binding affinities of FKBP12 inhibitors using the MM‐PBSA method[J]. Proteins: Structure, Function, and Bioinformatics, 2006, 64(4): 1058-1068. DOI:  10.1002/prot.21044.
    [25] Kumari R, Kumar R, Lynn A. g_mmpbsa—A GROMACS Tool for High-Throughput MM-PBSA Calculations[J]. Journal of Chemical Information and Modeling, 2014, 54(7): 1951-1962. DOI:  10.1021/ci500020m.
    [26] Gilson M K, Honig B. Calculation of the total electrostatic energy of a macromolecular system: solvation energies, binding energies, and conformational analysis[J]. Proteins: Structure, Function, and Bioinformatics, 1988, 4(1): 7-18. DOI:  10.1002/prot.340040104.
    [27] Kokkinidis M, Glykos N, Fadouloglou V. Protein flexibility and enzymatic catalysis[J]. Adv Protein Chem Struct Biol, 2012, 87: 181-218. DOI:  10.1016/B978-0-12-398312-1.00007-X.
    [28] DeLano W L. The PyMOL molecular graphics system. [CP]http://www.pymol.org,2002.
    [29] Laskowski R A, Swindells M B. LigPlot+: multiple ligand–protein interaction diagrams for drug discovery[J]. Journal of Chemical Information and Modeling, 2011, 51(10): 2 778-2 786. DOI:  10.1021/ci200227u.
    [30] Stefl S, Nishi H, Petukh M, et al. Molecular mechanisms of disease-causing missense mutations[J]. Journal of Molecular Biology, 2013, 425(21): 3 919-3 936. DOI:  10.1016/j.jmb.2013.07.014.
    [31] Merlino A, Mazzarella L, Carannante A, et al. The importance of dynamic effects on the enzyme activity X-ray structure and molecular dynamics of Onconase mutants[J]. Journal of Biological Chemistry, 2005, 280(18): 17 953-17 960. DOI:  10.1074/jbc.M501339200.
    [32] Chen J, Chen H, Shi Y, et al. Probing the effect of the non-active-site mutation Y229W in New Delhi metallo-β-lactamase-1 by site-directed mutagenesis, kinetic studies, and molecular dynamics simulations[J]. PLoS One, 2013, 8(12): e82080. DOI:  10.1371/journal.pone.0082080.
    [33] Liu S Q, Tao Y, Meng Z H, et al. The effect of calciums on molecular motions of proteinase K[J]. Journal of Molecular Modeling, 2011, 17(2): 289-300. DOI:  10.1007/s00894-010-0724-6.
  • [1] 桑鹏杨丕仁许丹朱月勋沈建新杨力权 . 碱性和中性食线虫真菌丝氨酸蛋白酶动力学行为差异研究*. 云南大学学报(自然科学版), doi: 10.7540/j.ynu.20180135
    [2] 佘连兵张文林李扬荣 . 非自治Kuramoto-Sivashinsky方程的半一致动力学. 云南大学学报(自然科学版), doi: 10.7540/j.ynu.20170626
    [3] 蒋丽红王亚明陕绍云贾庆明 . 非贵金属新型催化剂催化松香歧化动力学研究. 云南大学学报(自然科学版),
    [4] 杨汉春孙越梁洁 . 气体动力学零压流的狄拉克激波和真空相互作用的数值分析. 云南大学学报(自然科学版), doi: 10.7540/j.ynu.20130075
    [5] 张为杨斌蔡伊易俊年张杰 . 水溶液中NaF缔合行为的分子动力学研究. 云南大学学报(自然科学版), doi: 10.7540/j.ynu.20180726
    [6] 张耀宇张芳薛喜昌贾利群 . 相对运动动力学Chetaev型非完整系统Appell方程Mei对称性共形不变性和守恒量. 云南大学学报(自然科学版), doi: 10.7540/j.ynu.20170188
    [7] 彭守礼 . 符号动力学,星花积及其他. 云南大学学报(自然科学版),
    [8] 李彦敏 . 变质量非完整力学系统的共形不变性. 云南大学学报(自然科学版),
    [9] 肖雪王刚陈进 . 鸡嗉子榕果实和妃延腹榕小蜂性状的地理马赛克分化动力初析. 云南大学学报(自然科学版), doi: 10.7540/j.ynu.20170196
    [10] 赵晓华 . Lotka-Volterra方程:约化、分类及动力学性质. 云南大学学报(自然科学版),
    [11] 张捍卫许厚泽郑勇张超 . 地球内核的动力学理论及章动效应. 云南大学学报(自然科学版),
    [12] 张达敏蔡绍洪周海平郭长睿 . 规则网格中带人工免疫SIRS模型的动力学行为. 云南大学学报(自然科学版),
    [13] 张一方李艳梅 . 双星形成的非线性动力学机制和定性分析理论. 云南大学学报(自然科学版),
    [14] 戴正德 . (1+2)维斑图方程的动力学性态. 云南大学学报(自然科学版),
    [15] 张捍卫王庆林郭增长 . 三轴液核地球自转的动力学方程. 云南大学学报(自然科学版),
    [16] 陈晓萍包桂蓉王华李法社 . 蓝藻在超临界乙醇中液化反应动力学的研究. 云南大学学报(自然科学版),
    [17] 陈兴利 . 催化速差动力学分析法同时测定钼、钨. 云南大学学报(自然科学版), doi: 10.7540/j.ynu.20160194
    [18] 黄潇芮伟国 . 分数阶广义Bagley-Torvik方程的各种精确解及其动力学性质. 云南大学学报(自然科学版), doi: 10.7540/j.ynu.20170350
    [19] 金蕊钱正强刘飞虎杨明挚 . 工业大麻的营养吸收及其动力学特征. 云南大学学报(自然科学版), doi: 10.7540/j.ynu.20140426
    [20] 张通谷应丽杨汉春 . 气体动力学等熵流2个疏散波的相互作用. 云南大学学报(自然科学版),
  • 加载中
图(4)表(3)
计量
  • 文章访问数:  311
  • HTML全文浏览量:  243
  • PDF下载量:  14
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-01-04
  • 录用日期:  2020-05-15
  • 网络出版日期:  2020-06-23

基于分子动力学模拟和自由能计算的非小细胞肺癌对克唑替尼(crizotinib)耐药性研究

    作者简介:桑 鹏(1983– ),男,山东人,博士,讲师,主要从事结构生物信息学方面的研究. E-mail:speng431@163.com
    通讯作者: 宋玲, songl0870@yeah.net
  • 1. 大理大学 农学与生物科学学院,云南 大理 671003
  • 2. 昭通卫生职业学院,云南 昭通 657000

摘要: 为了研究间变性淋巴瘤激酶(anaplastic lymphoma kinase,ALK)F1174V突变引起非小细胞肺癌患者对crizotinib的耐药性机制,使用分子动力学模拟(molecular dynamics simulations, MD). 本质动力学(essential dynamics, ED)分析及分子力学泊松-玻尔兹曼比表面积(molecular mechanics Poisson-Boltzmann surface area, MM-PBSA)结合自由能计算等方法研究了野生型及F1174V突变型ALK与克唑替尼(crizotinib)结合模式的差异. MD模拟和ED分析结果表明,F1174V突变导致ALK的crizotinib结合口袋部位构象柔性降低,由此可能造成药物进出ALK速度减慢,影响药物发挥正常作用. MM-PBSA结合自由能计算结果表明,F1174V突变造成ALK与crizotinib的结合能力降低. 进一步构建了野生型和突变型ALK的自由能图谱(free energy landscape,FEL),并选取了二者最主要自由能井的代表性结构进行比较. 结果表明,F1174V突变可以显著改变ALK的构象,特别是可以通过收缩P-环和缩小结合口袋来减弱其与crizotinib的相互作用. 研究有助于揭示基于ALK的靶向药物耐药机制,并为开发抗非小细胞肺癌药物提供帮助.

English Abstract

  • 肺癌是威胁人类健康的主要疾病,可分为小细胞肺癌(small cell lung cancer,SCLC)和非小细胞肺癌(non-small cell lung cancer,NSCLC)两种类型,而后者可占全世界肺癌死亡人数的85%左右[1-3]. 目前研究表明,间变性淋巴瘤激酶(anaplastic lymphoma kinase,ALK)的重排是除表皮生长因子受体(epidermal growth factor receptor,EGFR)激活外另一引起非小细胞肺癌的重要原因. 因此,ALK成为了目前研发抗非小细胞肺癌药物的重要靶点.

    ALK是一种细胞质酪氨酸蛋白激酶,主要在中枢和外周神经系统中表达[4-5]. ALK属于胰岛素受体超家族成员,其与包括糖尿病和癌症在内的多种疾病相关. 与大多数蛋白激酶类似,ALK的结构采取了一种双叶折叠的模式,即包括了一个较小的氮末端叶(N-lobe)结构域和一个较大的碳末端叶(C-lobe)结构域. 以上两个结构域通过一个较小的环区相连接,该环区同时也是ALK与ATP分子及其他竞争性抑制剂进行结合的口袋. N-lobe结构域由5个反平行β-折叠及相邻的1个α-螺旋构成,而C-lobe则主要由α-螺旋所构成[5-6]. 研究表明,ALK基因可发生染色体重排并引起ALK融合基因的表达,而后者可以使ALK蛋白发生非配体依赖性的聚合导致其一直处于激活状态,从而引发包括非小细胞肺癌和大细胞淋巴瘤在内的多种癌症[5, 7]. 因此,有效抑制ALK的活性成为一种新的治疗非小细胞肺癌方式. 克唑替尼(crizotinib)是被发现的一种可以高效抑制ALK活性的化合物,并且很快被美国FDA批准为治疗非小细胞肺癌的药物而被广泛使用[8]. 虽然crizotinib在治疗ALK相关的非小细胞肺癌时疗效显著,但较好的疗效会随着肿瘤耐药性的产生而迅速减弱[9]. 研究表明,大约三分之一的ALK相关非小细胞肺癌患者在使用crizotinib 12个月内都会出现耐药性,从而会导致癌症的复发[5]. 因此,耐药性的产生已经成为制约开发ALK抑制剂和靶向药物的主要障碍[7]. 目前,人们已经在发生crizotinib耐药性的病人ALK基因上发现了多个突变,在这些突变中最为常见的是L1196M 和 G1269A突变,而且这两种突变的空间位阻耐药机制已经研究得相对比较清楚[7].

    Ou等在非小细胞肺癌患者中发现了另一种crizotinib耐药性ALK突变F1174V[10]. 与常见的L1196M 和 G1269A突变不同,F1174V突变导致crizotinib耐药性的机制仍不清楚. 为了揭示F1174V突变对ALK蛋白结构和动力学性质以及其与crizotinib相互作用的影响,我们分别使用了MD模拟,ED分析,MM-PBSA结合自由能计算以及FEL构建等方法对野生型和F1174V突变型ALK- crizotinib复合物进行了比较研究. 研究在蛋白质原子结构水平上揭示了突变是如何影响ALK与crizotinib相互作用的,并为突变导致的crizotinib耐药性机制提供合理解释.

    • 野生型ALK-crizotinib复合物三维空间结构下载自PDB蛋白质结构数据库(PDB检索号:2XP2[11]). 使用MODELLE软件包通过同源模建技术构建F1174V突变型ALK-crizotinib复合物三维空间模型[12]. 使用分子概率密度函数对所构建的结构进性评估,选择函数值最小的结构作为最终的突变体结构模型.

    • MD模拟及分子力学优化用GROMACS软件完成[13],选择的力场为Amber03[14-15]. Crizotinib的力场由PRODRG 2.5产生[16]. MD模拟之前,将能量最小化的ALK-crizotinib复合物野生型和突变体结构模型置于包含周期边界条件的立方体盒子中央,并保证蛋白质的每个原子与盒子壁的距离不小于1.1 nm. 之后,向盒子中加入TIP3P水分子[17]. 随后分别使用最陡下降法和共轭梯度法对系统进行能量最小化处理. 在MD模拟之前,对野生型和突变体模拟系统分别进行了100 ps的位置限制性模拟以促使水分子充分浸润蛋白质分子. 最后对两个模拟体系分别进行100 ns的生产MD模拟.

      MD模拟的详细参数如下:积分时间步长为2 fs;质量中心平移和旋转的消除频率为1步/次;非键相互作用的更新频率为10步/次;静电相互作用采用PME(Partial-Mesh Ewald)算法[18],其内插阶(interpolation order)设为4;傅立叶空间(Fourierspacing)和库仑(Coulomb)半径分别设为0.135 nm和1.0 nm;范德华(van der Waals,VDW)相互作用的截断半径为1.4 nm;蛋白质和非蛋白质组分(如溶剂和离子)的热浴温度均为300 K,耦合时间常数τ_t设为0.1 ps;系统压强为1 atm,耦合时间常数τ_p设为0.5 ps;原子的初始速度根据300 K条件下的麦克斯韦分布随机产生;用阶次(order)为4的LINCS算法[19]将键长约束在它们的平衡位置附近;记录结构框架的频率为10 ps/次. 模拟过程中的均方根偏差(root mean square deviation, RMSD)和均方根波动(root mean square fluctuation, RMSF)分别由GROMACS软件包中的g_rms和g_rmsf进行计算.

    • 野生型和突变体ALK-crizotinib复合物在模拟过程中的各种几何性质,如氢键数量(NHB)、溶剂可及性表面面积(SASA)、范德华接触数量(NNC)、分子回旋半径(Rg)分别由GROMACS软件包中g_hbond,g_sas, g_mindist, g_gyrate 等程序进行计算.

    • ED分析的基本原理及详细过程可以参考文献[20]. 使用GROMACS软件包中的g_covar程序构建野生型和突变型模拟体系的Cα原子协方差矩阵并进行对角线化,进一步使用g_anaeig程序将动力学模拟轨迹投射到对应的本征向量上进行ED分析.

      FEL可以用来表征一个模拟体系蛋白质的动态变化过程[21]. 利用ED分析的结果分别对野生型和突变型ALK模拟体系进行FEL构建. 自由能的计算公式为:

      $ \Delta G\left( X \right) = - {K_{\rm{B}}}T\; \ln P\left( X \right), $

      其中反应坐标X代表模拟轨迹在某个本征向量上的投射,KB代表玻尔兹曼常数,T代表系统的温度,P(X)代表某个本征向量上的构象分布概率. 在本研究中,最终的FEL由SigmaPlot软件作图显示. 为了得到每个自由能井的代表性结构,我们对自由能井里的所有构象根据RMSD值进行聚类分析,其中聚类的阈值为0.2 nm,使用每一个聚类的中间结构作为代表性结构.

    • 目前最为广泛使用的结合自由能算法为MM-PBSA算法[22]. 该算法通过对组成结合自由能的各组分进行分解和计算,可以为解析复合物分子之间复杂的相互作用过程提供帮助[23-24]. 在本研究中,野生型和突变体ALK 与crizotinib的结合自由能使用GROMACS软件包的插件g_mmpbsa程序来完成[25]. 在MM-PBSA算法中,结合自由能定义为[23, 25-26]:

      $ \Delta {G_{{\rm{binding}}}} = \Delta {G_{{\rm{complex}}}} - \left( {\Delta {G_{{\rm{protein}}}} + \Delta {G_{{\rm{ligand}}}}} \right), $

      其中ΔGcomplex, ΔGprotein 和 ΔGligand分别代表蛋白质–配体复合物,蛋白质及配体三者各自在溶剂环境中的总自由能. 对于上式中的每个组分,即自由能G,可以表示为:

      $ G = {E_{{\rm{MM}}}} + {G_{{\rm{sol}}}} - TS, $

      其中EMM代表在真空中的分子力学势能,Gsol定义了溶解自由能. TS代表自由能中的熵贡献部分,其中T代表温度,S代表熵. EMM包括了键作用和非键相互作用两类,即键角,二面角,静电相互作用(Eelec)及范德华相互作用(Evdw). 由于在单轨迹方法中键作用的贡献默认为0,因此g_mmpbsa在计算结合自由能时只计算EMM部分,即静电相互作用和范德华相互作用部分,即

      $ {E_{{\rm{MM}}}} = {E_{{\rm{ele}}}} + {E_{{\rm{vdw}}}}. $

      溶解自由能包括极性部分(Gpolar)和非极性部分(Gnonpolar),即:

      $ {G_{{\rm{solvation}}}} = {G_{{\rm{polar}}}} + {G_{{\rm{nonpolar}}}}. $

      虽然g_mmpbsa计算的结合自由能只是相对自由能,但这个软件还是非常适合比较配体在与不同蛋白受体结合时的自由能差异的. 本研究使用20~100 ns的模拟轨迹用于结合自由能的计算,其中每100 ps取出1个结构快照,即共使用800个结构快照用于结合自由能计算.

    • ALK-crizotinib复合物在模拟过程中的结构稳定性由随时间变化的碳骨架RMSD值来表征(图1). 如图1,对于野生型和突变体两个模拟系统,它们的RMSD值在模拟的前20 ns均有一个快速上升的过程,然后则变得较为平稳,因此我们选取最后80 ns模拟轨迹用于后续的结构属性分析.

      图  1  模拟过程中ALK野生型和突变型骨架原子相对初始结构的RMSD值.

      Figure 1.  The RMSD Values of backbone atoms from the initial structures of WT and mutant ALK during simulation

      为了比较野生型和突变体ALK-crizotinib复合物在模拟过程中构象柔性上的差异,我们计算了ALK各个残基在模拟过程中的RMSF值. 图2显示了以残基序号为函数的野生型和突变体ALK-crizotinib复合物分子的RMSF值. 由图2可以看出,除了数量有限的区域之外,突变体ALK在整体上比野生型蛋白具有更低的柔性(或更高的刚性). 野生型和突变体ALK所有残基RMSF的平均值分别为0.195,0.158 nm. 进一步观察图2显示突变体结构柔性显著降低的区域主要集中于暴露于蛋白分子表面的一些外部环区,特别是crizotinib的结合口袋区域. 例如,F1174V突变导致ALK构象柔性降低的区域包括残基1121~1129, 1135~1145, 1187~1192 及1197~1203. 在这些区域中,残基1121~1129区域是富含甘氨酸的环区(P-loop),且组成了crizotinib结合口袋的一部分,而残基1197~1203则组成了结合口袋的另外一部分.

      图  2  野生型和突变体残基在模拟过程中的RMSF值

      Figure 2.  The RMSF Value of eacd residue of WT and mutant during the simulation

      为了研究F1174V突变对crizotinib结合部位构象柔性的影响,我们比较了野生型和突变体ALK在相应残基位置上的RMSF差异. 如表1所示,与野生型相比,所有组成突变体crizotinib结合口袋的残基RMSF值均降低.

      残基均方根波动/nm
      野生型突变型
      11220.23460.2108
      11300.24420.0967
      11480.26150.1123
      11970.20960.0911
      11980.25030.2442
      11990.18100.1784
      12560.22190.2087

      表 1  组成crizotinib结合部位的残基RMSF值

      Table 1.  RMSF values of the crizotinib binding sites

      研究表明,蛋白质柔性在其发挥生物学功能时起着至关重要的作用[27]. 例如,较高的柔性会导致蛋白质底物结合口袋增大,因此有利于底物的进出. 另外,蛋白质构象柔性的增加也有利于提高其与底物的结合能力. 因此,F1174V突变导致的蛋白整体构象柔性特别是crizotinib结合部位构象柔性的降低,可能会影响ALK与crizotinib的相互作用.

      蛋白质构象柔性的降低往往是由于分子内部一些原子间相互作用及非键相互作用的增强所致的. 为了探究突变造成ALK蛋白构象柔性降低的原因,使用两个模拟系统的模拟平衡轨迹部分计算了ALK的各种动态几何属性. 表2描述了在模拟过程中野生型和突变体ALK各种动态几何属性的差异. 这些比较的几何属性包括氢键数量(Number of the H-bonds, NHB)、溶剂可及性表面面积(Solvent accessible surface area, SASA)、范德华接触的数量(Number of native contacts, NNC)以及旋转半径(Radius of gyration, Rg). 由表2可知,F1174V突变体比野生型ALK具有更多的NHB和NNC,但是具有较小的Rg和SASA,由此表明突变体比野生型ALK具有更多分子内相互作用以及更加紧密的结构包装.

      蛋白质NNCNHBRg/nmSASA/nm2
      野生型334530 (1094)227 (7.6)1.96 (0.02) 140.82 (2.82)
      突变型340622 (1125)232 (7.5)1.94 (0.01)135.26(2.75)

      表 2  野生型和突变体ALK在模拟过程中的几何属性

      Table 2.  Structural properties of the WT and mutant of ALK

    • 野生型和突变体ALK模拟轨迹的本质动力学分析表明只有极少数本征向量具有显著的本征值(图3). 对角线化处理野生型,F1174V突变体的Cα原子协方差矩阵后,得到的总均方波动值(total mean square fluctuation, TMSF)分别为5.71和 4.81 nm2,表明野生型ALK在模拟过程中经历了更大的构象波动. 此外,野生型ALK的前6个本征向量的本征值明显高于突变体对应的本征值(图3).

      图  3  野生型和突变体前20个本征向量对应的本征值

      Figure 3.  Eigenvalues as a function of the first 20 eigenvector of WT and mutant

      为了进一步比较野生型和突变体ALK在模拟过程中可以采到的构象空间大小,我们使用模拟轨迹在第1个本征向量(PC1)和第2个本征向量(PC2)的投射构建了二者的FEL(图4). 另外,为了研究F1174V突变对ALK结构的影响,我们使用聚类的方法分别抽提了野生型和突变体FEL中最主要自由能阱的代表性结构,并进行了详细的比较. 如图4所示,野生型ALK的FEL相比F1174V突变体更加宽阔,而且具有更少的自由能最小化区域,由此可以说明野生型ALK在模拟过程中相比突变体采集到了更加丰富构象空间.

      图  4  野生型和突变体ALK的FEL及自由能最小化区域的代表性结构.自由能图谱通过主成分分析方法进行构建,其中野生型和突变体的分子结构用 Pymoe软件[28]作图,ALK和crizontinib分子之间的相互作用由Ligplot软件[29]进行分析.

      Figure 4.  FEL and the representative structure in the lowest energy basin. The FEL were constructed based on the principal component analysis. The surface and cartoon representation of crizotinib were generated by Pymol[28]. Residue interations in the ALK-crizotinib interface was analyzed by Ligplot[29].

      为了探究构象柔性降低对ALK和crizotinib之间相互作用的影响,提取了野生型和突变体模拟系统的代表性结构,并进行比较. 野生型和F1174V突变体ALK与crizotinib之间的结合模式如图4所示. 对于野生型来讲,crizotinib的结合口袋大而宽阔,且存在5个相互作用残基和1对氢键. 对于图4进一步观察表明野生型ALK的P-loop采取了更加开放和伸展的构象. 与此相反,突变体ALK的crizotinib结合口袋小而狭窄,并且只存在3个相互作用残基,且并不存在氢键. 另外,突变体ALK采取了更加封闭和紧缩的构象. 突变体与crizotinib之间数量较少的相互作用残基以及氢键的缺失可能会降低二者之间的作用力,且较狭窄的结合口袋及更加封闭和紧缩的P-loop构象也不利于crizotinib进出ALK分子. 由此可以推测,ALK在发生F1174V突变后会显著减弱其与crizotinib的相互作用.

    • 为了定量描述F1174V突变对ALK与crizotinib之间相互作用的影响,我们使用MM-PBSA算法计算了二者之间的结合自由能,并详细解析了各自由能组成成分的差异. 如表3所示,野生型和突变体ALK与crizotinib之间的结合自由能分别为−94.38 kJ/mol和−70.29 kJ/mol,表明crizotinib与野生型结合的更加紧密,因此F1174V突变发生不利于ALK与crizotinib之间的结合. 进一步对结合自由能进行详细的分解表明,驱动ALK与crizotinib结合的主要作用力是真空中的势能(ΔGMM)和非极性的溶解自由能(ΔGnonpolar). 对于野生型ALK-crizotinib复合物,ΔGMM和ΔGnonpolar分别为−283.89和−8.76 kJ/mol,而对于突变型ΔGMM和ΔGnonpolar则分别为−250.61和−4.35 kJ/mol. 因此,ALK突变后会导致驱动其与crizotinib相互作用的作用力ΔGMM和ΔGnonpolar均降低. 与ΔGMM和ΔGnonpolar对ALK与crizotinib之间的结合自由能正贡献相反,极性的溶解自由能(ΔGpolar)和熵作用(−TΔS)均不利于ALK与crizotinib之间的结合. 对于野生型ALK-crizotinib复合物,ΔGpolar和−TΔS分别为173.42 kJ/mol和24.85 kJ/mol,而对于突变型ΔGpolar和−TΔS则分别为165.31 kJ/mol和19.36 kJ/mol. 值得注意的是,虽然突变型ALK较低的ΔGpolar和−TΔS值对于其与crizotinib之间的结合自由能负贡献更小,然而由于野生型ALK与crizotinib之间结合时具有更低的ΔGMM和ΔGnonpolar值,因此使其最终的结合自由能低于突变型,即ALK发生F1174V突变后会使其与crizotinib的结合能力降低,从而可能导致耐药性的产生.

      结合自由能组分野生型突变型
      ΔEvdw−41.52−39.56
      ΔEele−242.37−211.05
      ΔGpolar173.42165.31
      ΔGnonpolar−8.76−4.35
      ΔGMM−283.89−250.61
      ΔGsol164.66160.96
      ΔG−119.23−89.65
      TΔS24.8519.36
      ΔGbind−94.38−70.29

      表 3  野生型和突变体ALK与crizotinib之间的结合自由能

      Table 3.  Binding free energy and its components of the WT and mutant KJ/mol

    • 蛋白质氨基酸序列的变异往往会对其结构、动力学性质及功能产生重要的影响[30],特别是对于药物分子的靶蛋白来讲,错义突变的出现往往预示着其耐药性的产生. 为了研究F1174V突变影响ALK与crizotinib之间相互作用及产生耐药性的机制,分别对野生型和F1174V突变型ALK-crizotinib复合物进行了MD模拟,ED分析及MM-PBSA结合自由能计算比较研究.

      通过对MD模拟轨迹进行结构稳定性和构象柔性分析并结合ED分析结果表明,F1174V突变型ALK结构相比野生型更加紧凑,且刚性更强. 进一步对组成crizotinib结合口袋的氨基酸残基分析表明,突变型ALK中所有对应的残基柔性均有所降低. 为了探究ALK在发生F1174V突变后蛋白质构象柔性降低的原因,分别计算了野生型和突变型ALK的各种分子内相互作用,结果表明ALK在发生突变后氢键数量及范德华接触数量等表征分子间相互作用强弱的指标均有所增加,进而导致ALK分子的构象柔性降低.

      目前研究表明,蛋白质功能的发挥依赖于其动力学性质及构象柔性[31]. 较高的构象柔性可以使蛋白质能够以不同的构象状态存在,例如可以很容易实现底物结合口袋的舒张,使底物可以更方便地进出结合口袋,从而提高了与底物结合的效率[32]. 另外,现有研究表明蛋白质构象柔性的增加也可以增强其与底物之间的相互作用,进而使二者之间的结合能力增强[33]. 使用计算生物学的方法对野生型和F1174V突变型ALK分别与crizotinib之间的相互作用进行了比较研究,研究结果表明ALK在发生突变后构象柔性降低,并且与野生型相比具有更狭窄的crizotinib结合口袋,因此降低了ALK与crizotinib之间的结合效率. 另外,通过计算ALK与crizotinib之间的结合自由能表明,ALK在发生F1174V突变后导致其与crizotinib之间的结合自由能升高,因此结合能力降低. 综上所述,研究结果从结构生物学信息学的角度证实了 F1174V突变发生后不利于ALK与crizotinib之间的相互作用,从而解释了当非小细胞肺癌患者ALK蛋白发生F1174V突变后产生crizotinib耐药性的机制.

参考文献 (33)

目录

    /

    返回文章
    返回