一种利用计算机模拟蛋白质相互作用的方法

摘要

本发明公开了一种利用计算机模拟蛋白质相互作用的方法。该方法利用NAMD_2.6_Linux－i686计算软件和VMD－1.8.5分析软件在PC机上进行分子动力学模拟，首先将PDB数据库的蛋白质结构制成适当的研究体系，然后对该体系进行能量最小化处理，将体系加热到设定温度后，再对体系进行多步平衡处理，接着对体系实施分子动力学模拟，最后研究蛋白质在设定条件下所有原子的运动轨迹、构象变化和变化过程中起关键作用的氨基酸残基。该方法不需借助于大型计算机和高性能计算中心，设备要求低，CPU利用率高，所用机时短并且结果准确可靠，便于在与分子识别和蛋白质动力学行为相关的生命科学和物理化学以及医药领域广泛应用。

说明书

技术领域

本发明涉及一种利用计算机模拟蛋白质相互作用的方法，尤其涉及一种利用NAMD_2.6和VMD-1.8.5软件研究蛋白质的转录因子之间相互作用的方法。

背景技术

蛋白质相互作用是生命活动过程的关键步骤，如多转录因子的协同调控、蛋白质的折叠及酶的合成、分子识别及信号分子的转导、免疫反应等。常规的实验方法采用测定反应平衡常数来检测蛋白质之间的结合强度，也可通过原子力显微镜，对单分子进行操作，将蛋白质彼此分开为解离过程提供信息。但是无论是测定反应平衡常数还是用原子力显微镜进行操作都不能很好的揭示蛋白质间的非键相互作用，而分子动力学模拟的方法可以弥补上述不足。不仅如此，通过分子动力学模拟还能得到蛋白质相互作用的具体过程(每个原子的运动轨迹)以及该过程中的许多重要的动力学参数。

自1977年McCammon J.A.等首次用分子动力学的方法研究蛋白质至今，分子动力学模拟的计算方法迅速发展，模拟的体系从最初的数百个原子几皮秒的尺度到今天260多万个原子几百个纳秒的尺度，模拟的对象包括蛋白质、DNA、RNA等。但是，这些分子动力学模拟的过程大都是在大型计算机或者高性能计算中心完成的，而利用NAMD_2.6_Linux-i686计算软件和VMD-1.8.5分析软件在PC机上进行分子动力学模拟，目前国内外未见报道。

发明内容

针对现有技术的不足，本发明要解决的问题是提供一种利用NAMD_2.6_Linux-i686计算软件和VMD-1.8.5分析软件(http：//www.ks.uiuc.edu/Research/namd/)在普通PC机上模拟蛋白质相互作用的方法。

本发明所述的方法，按如下步骤进行：

(1)从PDB数据库获得含有一条或多条肽链的蛋白质结构；

(2)利用VMD-1.8.5软件提取欲研究的目标肽链，将其溶解到盛有0.9％NaCl溶液的长方体水箱(Water Box)中，并且使其距离水箱各壁的最小距离为10，得到肽-盐溶液体系，将该体系制成含有非氢原子空间位置数据的结构，命名为ionized.pdb，同时将该体系制成含有氢原子以及氢键参数的结构，命名为ionized.psf；

(3)利用NAMD_2.6_Linux-i686软件，采用Charmm27力场，用Steepest Descent方法，在周期性边界条件(Periodic Boundary Conditions PBC)下对上述肽-盐溶液体系进行总时间长度为20皮秒(ps)的能量最小化(Minimize)处理，在此过程中每2飞秒(fs)计算一次以各原子为球心、半径为10的球形空间内的其它原子对该原子的范德华力，每4fs计算一次以各原子为球心、半径为12的球形空间内的其它原子对该原子的电场力，计算完成后得到各原子的运动轨迹，命名为minimize.dcd；

(4)然后，固定肽链中的所有原子，采用与(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，使用Langevin动力学(Langevin Dynamics)控温方法加热20ps的时间，结果体系逐渐升温并稳定到280-320K，得到该过程中各原子的运动轨迹，命名为heat.dcd；

(5)依次用100kcal/mol·²、50kcal/mol·²、20kcal/mol·²和10kcal/mol·²的和谐约束作用(Harmonic Constraint)约束上述肽链，采用与(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，使用Langevin动力学(Langevin Dynamics)控温方法和Nose-Hoover Langevin piston控压方法将温度控制在300K，压强控制在1atm，分别平衡体系20ps，依次分别得到各步中所有原子的运动轨迹，依次分别命名为equ100.dcd、equ50.dcd、equ20.dcd和equ10.dcd；

(6)去除(5)中对肽链的和谐约束，采用与(5)相同的力场、周期性边界条件、范德华力和电场力计算方法、控温控压方法，将PME grid(Particle Mesh Ewald grid)边长设置为1±0.1，利用PME长程电场力计算方法，固定体系中所有的氢氧键长，对整个体系进行2纳秒(ns)的分子动力学模拟，得到该过程中各原子的运动轨迹，命名为simulate.dcd；

(7)将上述(3)-(6)中的计算结果minimize.dcd、heat.dcd、equ100.dcd、equ50.dcd、equ20.dcd、equ10.dcd和simulate.dcd载入VMD-1.8.5软件观察所有原子的运动轨迹；

(8)通过上述运动轨迹分析蛋白质的构象变化过程；

(9)找出蛋白质构象变化过程中起关键作用的氨基酸残基。

上述的利用计算机模拟蛋白质相互作用的方法中：步骤(4)中所述的体系稳定温度优选为310K。

利用本发明所述的方法进行蛋白质相互作用的模拟与传统方法相比具有明显的优越性：

(1)对范德华力和电场力的计算精度高，计算范围分别为10和12，而传统方法大都使用9；

(2)不需借助于大型计算机和高性能计算中心，设备要求低，在普通PC上即可进行；

(3)CPU利用率高，在没有干扰的情况下利用率可达到95％以上；

(4)所用机时短，对于3万到4万个原子的体系在双核Pentium4 3.0G Hz处理器、1G内存、Redfleg Desktop 5.0操作系统的PC机上约一周时间即可；

(5)结果准确可靠；

(6)便于在与分子识别和蛋白质动力学行为相关的生命科学和物理化学以及医药领域广泛应用。

附图说明

图1 HMG结构域和POUs结构域的初始态和终态构象

图中A为模拟之前Sox的HMG结构域和Oct4的POUs结构域的相对位置和构象，B为模拟之后两者的相对位置和构象(HMG结构域在上面，POUs结构域在下面)。

具体实施方式

下面结合实施例对本发明作进一步说明：

实施例1：

(1)从PDB数据库(http：//www.rcsb.org/pdb/home/home.do)获得代码为1o4x的蛋白质结构命名为1o4x.pdb；

(2)利用VMD-1.8.5软件提取欲研究的目标肽链，将其溶解到盛有0.9％NaCl溶液的长方体水箱(Water Box)中，并且使其距离水箱各壁的最小距离为10，得到肽-盐溶液体系，将该体系制成含有非氢原子空间位置数据的结构，命名为1o4x_ionized.pdb，同时将该体系制成含有氢原子以及氢键参数的结构，命名为1o4x_ionized.psf；

(3)利用NAMD_2.6_Linux-i686软件，采用Charmm27力场，用Steepest Descent方法，在周期性边界条件(Periodic Boundary Conditions PBC)下对上述肽-盐溶液体系进行总时间长度为20皮秒(ps)的能量最小化(Minimize)处理，在此过程中每2飞秒(fs)计算一次以各原子为球心、半径为10的球形空间内的其它原子对该原子的范德华力，每4fs计算一次以各原子为球心、半径为12的球形空间内的其它原子对该原子的电场力，计算完成后得到各原子的运动轨迹，命名为1o4x_minimize.dcd；

(4)然后，固定肽链中的所有原子，采用与(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，使用Langevin动力学(Langevin Dynamics)控温方法加热20ps的时间，结果体系逐渐升温并稳定到310K，得到该过程中各原子的运动轨迹，命名为1o4x_heat.dcd；

(5)依次用100kcal/mol·²、50kcal/mol·²、20kcal/mol·²和10kcal/mol·²的和谐约束作用(Harmonic Constraint)约束上述肽链，采用与(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，使用Langevin动力学(Langevin Dynamics)控温方法和Nose-Hoover Langevin piston控压方法将温度控制在300K，压强控制在1atm，分别平衡体系20ps，依次分别得到各步中所有原子的运动轨迹，依次分别命名为1o4x_equ100.dcd、1o4x_equ50.dcd、1o4x_equ20.dcd和1o4x_equ10.dcd；

(6)去除(5)中对肽链的和谐约束，采用与(5)相同的力场、周期性边界条件、范德华力和电场力计算方法、控温控压方法，将PME grid(Particle Mesh Ewald grid)边长设置为1±0.1，利用PME长程电场力计算方法，固定体系中所有的氢氧键长，对整个体系进行2纳秒(ns)的分子动力学模拟，得到该过程中各原子的运动轨迹，命名为1o4x_simulate.dcd；

(7)将上述(3)-(6)中的计算结果1o4x_minimize.dcd、1o4x_heat.dcd、1o4x_equ100.dcd、1o4x_equ50.dcd、1o4x_equ20.dcd、1o4x_equ10.dcd和1o4x_simulate.dcd载入VMD-1.8.5软件观察所有原子的运动轨迹；

(8)通过上述运动轨迹分析蛋白质的构象变化过程，初始态和终态构象见图1；

(9)找到蛋白质构象变化过程中起关键作用的氨基酸残基HMG的ASP274和POUs的LYS18。

实施例2：

(1)从PDB数据库(http：//www.rcsb.org/pdb/home/home.do)获得代码为1JLU的蛋白质结构命名为1JLU.pdb；

(2)利用VMD-1.8.5软件提取欲研究的目标肽链，将其溶解到盛有0.9％NaCl溶液的长方体水箱(Water Box)中，并且使其距离水箱各壁的最小距离为10，得到肽-盐溶液体系，将该体系制成含有非氢原子空间位置数据的结构，命名为1JLU_ionized.pdb，同时将该体系制成含有氢原子以及氢键参数的结构，命名为1JLU_ionized.psf；

(3)利用NAMD_2.6_Linux-i686软件，采用Charmm27力场，用Steepest Descent方法，在周期性边界条件(Periodic Boundary Conditions PBC)下对上述肽-盐溶液体系进行总时间长度为20皮秒(ps)的能量最小化(Minimize)处理，在此过程中每2飞秒(fs)计算一次以各原子为球心、半径为10的球形空间内的其它原子对该原子的范德华力，每4fs计算一次以各原子为球心、半径为12的球形空间内的其它原子对该原子的电场力，计算完成后得到各原子的运动轨迹，命名为1JLU_minimize.dcd；

(4)然后，固定肽链中的所有原子，采用与(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，使用Langevin动力学(Langevin Dynamics)控温方法加热20ps的时间，结果体系逐渐升温并稳定到320K，得到该过程中各原子的运动轨迹，命名为1JLU_heat.dcd；

(5)依次用100kcal/mol·²、50kcal/mol·²、20kcal/mol·²和10kcal/mol·²的和谐约束作用(Harmonic Constraint)约束上述肽链，采用与(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，使用Langevin动力学(Langevin Dynamics)控温方法和Nose-Hoover Langevin piston控压方法将温度控制在300K，压强控制在1atm，分别平衡体系20ps，依次分别得到各步中所有原子的运动轨迹，依次分别命名为1JLU_equ100.dcd、1JLU_equ50.dcd、1JLU_equ20.dcd和1JLU_equ10.dcd；

(6)去除(5)中对肽链的和谐约束，采用与(5)相同的力场、周期性边界条件、范德华力和电场力计算方法、控温控压方法，将PME grid(Particle Mesh Ewald grid)边长设置为1±0.1，利用PME长程电场力计算方法，固定体系中所有的氢氧键长，对整个体系进行2纳秒(ns)的分子动力学模拟，得到该过程中各原子的运动轨迹，命名为1JLU_simulate.dcd；

(7)将上述(3)-(6)中的计算结果1JLU_minimize.dcd、1JLU_heat.dcd、1JLU_equ100.dcd、1JLU_equ50.dcd、1JLU_equ20.dcd、1JLU_equ10.dcd和1JLU_simulate.dcd载入VMD-1.8.5软件观察所有原子的运动轨迹；

(8)通过上述运动轨迹分析蛋白质的构象变化过程；

(9)找到蛋白质构象变化过程中起关键作用的氨基酸残基。

实施例3：

(1)从PDB数据库(http：//www.rcsb.org/pdb/home/home.do)获得代码为1HN2的蛋白质结构命名为1HN2.pdb；

(2)利用VMD-1.8.5软件提取欲研究的目标肽链，将其溶解到盛有0.9％NaCl溶液的长方体水箱(Water Box)中，并且使其距离水箱各壁的最小距离为10，得到肽-盐溶液体系，将该体系制成含有非氢原子空间位置数据的结构，命名为1HN2_ionized.pdb，同时将该体系制成含有氢原子以及氢键参数的结构，命名为1HN2_ionized.psf；

(3)利用NAMD_2.6_Linux-i686软件，采用Charmm27力场，用Steepest Descent方法，在周期性边界条件(Periodic Boundary Conditions PBC)下对上述肽-盐溶液体系进行总时间长度为20皮秒(ps)的能量最小化(Minimize)处理，在此过程中每2飞秒(fs)计算一次以各原子为球心、半径为10的球形空间内的其它原子对该原子的范德华力，每4fs计算一次以各原子为球心、半径为12的球形空间内的其它原子对该原子的电场力，计算完成后得到各原子的运动轨迹，命名为1HN2_minimize.dcd；

(4)然后，固定肽链中的所有原子，采用与(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，使用Langevin动力学(Langevin Dynamics)控温方法加热20ps的时间，结果体系逐渐升温并稳定到280K，得到该过程中各原子的运动轨迹，命名为1HN2_heat.dcd；

(5)依次用100kcal/mol·²、50kcal/mol·²、20kcal/mol·²和10kcal/mol·²的和谐约束作用(Harmonic Constraint)约束肽链，采用与(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，使用Langevin动力学(Langevin Dynamics)控温方法和Nose-Hoover Langevin piston控压方法将温度控制在300K，压强控制在1atm，分别平衡体系20ps，依次分别得到各步中所有原子的运动轨迹，依次分别命名为1HN2_equ100.dcd、1HN2_equ50.dcd、1HN2_equ20.dcd和1HN2_equ10.dcd；

(6)去除(5)中对肽链的和谐约束，采用与(5)相同的力场、周期性边界条件、范德华力和电场力计算方法、控温控压方法，将PME grid(Particle Mesh Ewald grid)边长设置为1±0.1，利用PME长程电场力计算方法，固定体系中所有的氢氧键长，对整个体系进行2纳秒(ns)的分子动力学模拟，得到该过程中各原子的运动轨迹，命名为1HN2_simulate.dcd；

(7)将上述(3)-(6)中的计算结果1HN2_minimize.dcd、1HN2_heat.dcd、1HN2_equ100.dcd、1HN2_equ50.dcd、1HN2_equ20.dcd、1HN2_equ10.dcd和1HN2_simulate.dcd载入VMD-1.8.5软件观察所有原子的运动轨迹；

(8)通过上述运动轨迹分析蛋白质的构象变化过程；

(9)找出蛋白质构象变化过程中起关键作用的氨基酸残基。