刘永刚 杨云帆 唐 妍 徐建洁 刘 冰
(1. 西南科技大学分析测试中心 四川绵阳 621010; 2. 西南科技大学材料科学与工程学院 四川绵阳 621010; 3. 国民核生化灾害防护国家重点实验室 北京 102205)
随着我国经济的发展,城市人口越来越密集,以投放生物化学毒剂为主的恐怖袭击将成为社会稳定和经济发展的巨大隐患[1]。因此,从理论上研究毒剂分子结构和发展生化毒剂检测技术成为重大事故后果分析的重要内容。定性地描述一个可能发生或正在发生的重大化学事故对人员和环境造成危害的程度,预测和评估危害后果,对突发危害性公共卫生事件的应急处置具有重要意义。
由分子光谱特性可知,不同的物质有不同的能级轨道,有其特定的精细结构,其表现出来的光谱性质是不同的,因此可利用分子光谱的特性来对物质进行甄别和鉴别,从而进行定性定量分析,故分子光谱亦称为“指纹性光谱”。采用光谱学检测方法检测危险化学品具有采样灵活、样品无需前处理、分析速度快、甚至可以原位检测等优点。这些特点使光谱法在检测技术中占有重要地位,成为一种应用最广泛的方法[2]。红外光谱具有“灵敏度高”“指纹性强”“原位测定”等优点,在分析复杂化学物质组成及结构信息等方面有着广泛的应用,特别是利用“红外遥感”来实时监测大气环境污染、无人区油气田泄露等方面表现出红外光谱技术的优越性,可对有毒有害物质进行“指认”[3-5]。随着红外光谱技术的发展,该技术在光谱学界越来越受到重视。
生物体系的毒剂由于其庞大且复杂的结构信息,仅通过实验测定其光谱容易忽略一些重要的结构信息,计算模拟是对传统实验方法的一种补充,能够了解一些通过传统实验手段不能得到的信息[6-8]。量子化学计算方法是在量子力学基础上发展起来的与量子化学相关的计算方法[9-11]。随着计算理论和计算机技术的不断发展,计算模拟在光谱模拟分析领域的精确度越来越高,应用范围也越来越广,已经成为一种非常重要的研究光谱产生机理的方法[12]。
本文以葡萄球菌肠毒素蛋白为研究对象,通过分子动力学计算方法模拟该蛋白在不同温度条件下的动态过程,并通过傅里叶变换方法建立大体系毒素分子的红外光谱计算模型,模拟得到葡萄球菌肠毒素的红外光谱,并分析其振动光谱与温度之间的关系。
分子动力学模拟计算充当了微观长度、时间尺度和宏观世界之间的桥梁(见图1),根据分子间相互作用提出假设,预测分子的“真实”综合特性。在分子动力学模拟过程中可以根据“实际意想”的环境调整系统性质的准确性,此预测具有一定的“真实”性。模拟仿真方法也充当理论与实验之间的桥梁。通过使用相同的模型进行模拟,可以测试理论的正确性,还可以模拟在实验室中难以执行或无法执行的实验(例如,需要在极端温度或压力下完成的实验)[13-15]。对于大体系物质葡萄球菌毒素可采用分子动力学方法计算能量与时间的关系,并利用傅里叶变换得到红外光谱[16-17]。
葡萄球菌肠毒素 B 蛋白的初始构象取自蛋白质结构数据库(PDB号为3SEB)[18],该结构共有 238 个氨基酸,使用 Discovery Studio 软件补充氢原子后共有3 921个原子。通过添加合适的水分子和其他离子来保证整个系统稳定且为电中性。在分子动力学模拟的整个过程中,仅记录目标蛋白的运动和受力等信息,系统内添加的水分子和其他离子用于维持系统的稳定和目标蛋白的稳定构型,不会对计算目标蛋白的红外光谱产生影响。对葡萄球菌肠毒素 B 蛋白系统采用Charmm 36 m力场的方法,共有66 080个原子,系统结构如图 1 所示。
图1 葡萄球菌肠毒素模型示意图Fig.1 Schematic diagram of the staphylococcal enterotoxin model
对于分子动力学运动的计算,在给出所有的原子坐标ri和势能函数U(rN)的情况下,接下来要做的是计算原子力fi,其定义式为:
(1)
首先通过能量最小化进行结构弛豫以保证葡萄球菌肠毒素B系统里面没有立体几何冲突或不当,在能量最小化的过程中采用最陡下降法(steep)进行5 000步能量优化,对于氢键位置使用 Lincs 方法进行约束。通过Nosw-Hoover的热浴方法进行NVT(NVT为正则系综简写,表示具有确定的粒子数(N)、体积(V)、温度(T))系综 125 ps 模拟,将体系温度分别设定在254,264,274,299和314 K,模拟稳定系统的温度,在每个稳定的温度下均进行 125 ps 的分子动力学模拟。此时通过Parrinello-Rahman耦合的方法将系统的压强稳定在1.05×105Pa,还需要在NPT(NPT为等温等压系综简写,表示具有确定的粒子数(N)、压强(P)、温度(T))条件下进行模拟稳定系统的压强,模拟时长为125 ps, 将系统的压强稳定在1.05×105Pa。随着两个平衡阶段的完成,系统现在已经处于目标温度和压强下,再进行时长为 1 ns 的分子动力学模拟(MD),模拟步长为 2 fs,每隔 0.2 ps 记录一次数据(坐标、能量、力和速度)。
记录根据时间变化的总偶极矩信号,做自相关函数,然后做傅里叶变换,得到红外光谱。公式为:
I(ω)n(ω)=量子校正系数×
(8)
式中:I(ω)为吸收强度;ω为频率;M为体系的总偶极矩;〈M(t)·M(0)〉是体系总偶极矩的自相关函数;n(ω)是折射率。
对目标蛋白能量最小化后,将系统的温度设定到各目标温度并维持稳定。在最终 1 ns 的分子动力学模拟中,监测各温度下系统的多种参数,结果如图 2、图 3 所示。
图2 葡萄球菌肠毒素B系统温度随时间的变化Fig.2 Temperature changes of staphylococcal enterotoxin B system over time
图3 不同温度下葡萄球菌肠毒素B系统能量随时间的变化Fig.3 Energy changes of staphylococcal enterotoxin B system over time at different temperatures
从图 2 可知,葡萄球菌肠毒素 B 系统处于比较完美的目标温度,其平均值和波动均在合理范围内。图 3 是不同温度下葡萄球菌肠毒素 B 系统能量随时间的变化曲线图,系统在各温度条件下能量都较为稳定,波动处于合理范围内,并且从图 3 可知平均能量与温度成正比,即温度越高,系统的能量越高。为研究葡萄球菌肠毒素 B 蛋白在不同温度下的稳定性,分析了在模拟计算过程中蛋白质骨架相对于平衡结构的均方根偏差( Root Mean Square Deviation,RMSD)以及蛋白质回旋半径(Rg),结果如图 4 和图 5 所示[19]。
图4 不同温度下葡萄球菌肠毒素B系统RMSD随时间的变化Fig.4 Changes of RMSD of staphylococcal enterotoxin B system over time at different temperatures
图5 不同温度下葡萄球菌肠毒素B蛋白质Rg随时间的变化Fig.5 Changes of Rg of staphylococcal enterotoxin B protein over time at different temperatures
RMSD 反映了两个结构间的偏离程度,该值越大,表明结构间差异越大;Rg表征了蛋白质的致密程度,在一定程度上体现了蛋白质在空间分布的延展性。图4显示,在不同温度下,蛋白质骨架的 RMSD 在模拟初始阶段(0~0.5 ns)均显著增加,表明结构发生较大变化。随着模拟时间进行,在 254,264,274,314 K温度下,系统的 RMSD 均保持在较小恒定值(0.15 nm)附近。但是在 299 K 时,RMSD 最终保持在 0.275 nm 附近,说明温度对葡萄球菌肠毒素B 蛋白结构的稳定性有很大影响。这也与图 3 中温度为 299 K 时系统的能量较高相呼应。在图 5 中,温度较低时(254 ,264 ,274 K),Rg保持在一个恒定值,说明此时的蛋白质折叠稳定。在 314 K 时,其Rg也保持在一个相对恒定值,但是其波动涨落比低温下更大。在 299 K 时,其Rg发生明显变化,说明此时蛋白质处于一个展开不稳定状态,这与前面分析一致。综上所述,温度对蛋白质结构稳定性有明显影响,尤其是温度为 299 K 时变化最明显。
对上述系统再次进行不采用任何约束力的模拟,时长 200 ps,步长为 2 fs,每隔1步记录数据(位置、力、能量和速度),结果共有 200 000 个构象。然后使用 Gromacs 程序和傅里叶变换,记录系统内目标蛋白的运动信息和偶极矩随时间变化的曲线,然后做自相关函数的傅里叶变换,即将时间-偶极信号转变为频率-吸收强度信号,计算不同温度下的红外光谱,结果如图 6 和表 1 所示。
对葡萄球菌肠毒素蛋白进行红外光谱测定,对比实验和理论结果。在干燥条件下称取葡萄球菌肠毒素蛋白样品 2 mg,按照质量比 1∶100 加入烘干的 KBr,用玛瑙研钵研磨至细粉状,压片(将混合研磨后的样品置于压片夹内,压力约为 10 MPa,压片时间约为 60 s)后, 用 Frontier(美国PE公司)红外光谱仪分析,分辨率 4 cm-1,以 KBr 为空白,并去除背景值,收集 4 000 ~ 500 cm-1范围内光谱。
图6 不同温度下葡萄球菌肠毒素B蛋白质的红外光谱与实验红外光谱Fig.6 Computational and experimental infrared spectroscopy of staphylococcal enterotoxin B protein at different temperatures
表1 不同温度下葡萄球菌肠毒素B蛋白质的红外光谱的峰值Table 1 Infrared spectrum peaks of staphylococcal enterotoxin B protein at different temperatures
表2 实验与计算得到的葡萄球菌肠毒素红外特征峰及其振动模式Table 2 Experimental and calculated infrared characteristic peaks and vibration modes of staphylococcal enterotoxin
采用分子动力学方法对葡萄球菌肠毒素蛋白进行非简谐振动分析计算,建立了生物大分子红外光谱的理论计算模型,计算其特征峰位并研究了温度对其稳定性及振动情况的影响,与实验数据对比得到如下结论:通过分子动力学模拟了葡萄球菌肠毒素蛋白的动态过程,得到了不同温度条件下系统的能量、均方根偏差(RMSD)和蛋白质回旋半径(Rg)随时间变化图。分析得知葡萄球菌毒素平均能量与温度成正比,即温度越高,系统的能量越高,温度对蛋白质结构稳定性有明显影响。计算表明温度对其高频振动光谱影响更大。
利用此数据和理论实践方法,后续拟将该模型应用于其他化学战剂和危化品的近红外光谱计算中,将大量的毒素低频振动光谱收集成库,即可以实现利用其光谱特征信号来实现对该类物质的分子识别和精确预警。