一种利用计算机模拟计算蛋白质与DNA之间作用力的方法

摘要

本发明公开了一种利用计算机模拟计算蛋白质与DNA之间相互作用力的方法。该方法利用NAMD_2.6计算软件和VMD－1.8.6分析软件进行分子动力学模拟，首先将PDB数据库的蛋白质结构制成适当的研究体系，然后对该体系进行能量最小化处理，将体系加热到设定温度后，对体系进行多步平衡处理，接着对体系实施分子动力学模拟，然后标记体系中目标研究对象并对体系进行SMD模拟，最后处理和分析数据得到蛋白质与DNA之间的相互作用力变化情况，并求得它们之间的最大相互作用力。该方法测定对象范围广泛、实验材料易于制备、设备要求低并且利用率高，便于在与分子识别和基因表达行为相关的生命科学和物理化学以及医药领域广泛应用。

说明书

技术领域

本发明涉及一种利用计算机模拟计算蛋白质与DNA之间相互作用力大小的方法，尤其涉及一种利用NAMD_2.6和VMD-1.8.6软件定量研究蛋白质与DNA之间相互作用力的方法。

背景技术

蛋白质与DNA相互作用是生命活动过程的关键步骤，如转录因子对DNA上特定位点的识别、启动或关闭某个特定基因的转录，调控DNA的复制等。常规的实验方法采用测定反应平衡常数来检测蛋白质与DNA之间的结合强度，也可通过原子力显微镜，对单分子进行操作，将蛋白质与DNA分开为解离过程提供信息。但是由于实验材料难以制备，无论是测定反应平衡常数还是用原子力显微镜进行操作都不能很好的揭示蛋白质与DNA之间相互作用力的大小，而分子动力学模拟的方法可以弥补上述不足。不仅如此，通过分子动力学模拟还能得到蛋白质与DNA相互作用的具体过程(每个原子的运动轨迹)以及该过程中的许多重要的动力学参数。

自1977年McCammon J.A.等首次用分子动力学的方法研究蛋白质至今，分子动力学模拟的计算方法迅速发展，模拟的体系从最初的数百个原子几皮秒的尺度到今天260多万个原子几百个纳秒的尺度，模拟的对象包括蛋白质、DNA、RNA等，应用范围涉及制药和生命科学的众多领域。但是，利用NAMD_2.6计算软件和VMD-1.8.6分析软件通过对蛋白质与DNA进行分子动力学模拟从而计算彼此之间相互作用力大小的方法，目前国内外未见报道。

发明内容

针对现有技术的不足，本发明要解决的问题是提供一种利用NAMD_2.6计算软件和VMD-1.8.6分析软件(http://www.ks.uiuc.edu/Research/namd/)通过模拟蛋白质与DNA的分子动力学行为从而计算得到它们之间相互作用力大小的方法。

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

(1)从PDB数据库获得含有由蛋白质和DNA所形成复合物的结构数据文件；

(2)利用VMD-1.8.6软件提取欲研究的目标蛋白质和DNA片断，将其溶解到盛有0.9％NaCl溶液的长方体水箱(Water Box)中，得到蛋白质-DNA-盐溶液体系；所述体系中含有非氢原子空间位置数据的结构，命名为ionized.pdb，含有氢原子以及氢键参数的结构，命名为ionized.psf；

其中：上述水箱内尺寸以蛋白质和DNA片断距离水箱各壁的最小距离为12，且蛋白质沿DNA-蛋白质质心连线方向能移动不少于50为准；

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

(4)用20kcal/mol·²的和谐约束作用(Harmonic Constraint)约束蛋白质的α碳原子(CA atoms)，并用10kcal/mol·²的和谐约束作用约束DNA骨架原子(nucleic backboneatoms)，采用与步骤(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，并且使用Langevin动力学(Langevin Dynamics)控温方法和Nose-Hoover Langevin piston控压方法将温度控制在280K-320K，压强控制在1atm，对步骤(3)能量最小化处理后的体系加热100皮秒(ps)，使体系逐渐升温并稳定到280K-320K；

(5)对步骤(4)加热后的体系，用10kcal/mol·²的和谐约束作用约束蛋白质的α碳原子，其它条件与步骤(4)相同，进行100皮秒(ps)的分子动力学模拟；

(6)去除对蛋白质的α碳原子的约束，其它条件同步骤(4)，对步骤(5)模拟后的体系进行100皮秒(ps)的分子动力学模拟；

(7)用10kcal/mol·²的和谐约束作用约束DNA片段两个末端各一个碱基对的所有原子，固定体系中所有的氢氧键长，将PME(Particle Mesh Ewald)格子边长设置为1±0.1，采用与步骤(4)相同的力场、周期性边界条件、范德华力和电场力计算方法、控温控压方法、PME长程电场力计算方法，对步骤(6)模拟后的体系进行2-6纳秒(ns)的分子动力学模拟，使得体系的RMSD(Root Mean Square Distance)值达到水平；

(8)将步骤(7)模拟后体系中的DNA骨架原子固定，以蛋白质的α碳原子为SMD(steered molecular dynamics)原子，将SMD系数(SMDk)设定为0.1-1kcal/mol·²的任意数值，SMD移动速度(SMDVel)设定为0.0005/fs，SMD记录频率(SMDOutputFreq)设定为100，其它条件与步骤(7)相同，对步骤(7)模拟后的体系进行100皮秒(ps)的SMD模拟，得到该过程中各记录点的热力学参数，命名为smd.log；

(9)用VMD-1.8.6软件处理文件smd.log，得到各记录点的力，命名为ft.dat；

(10)以各记录点的力为纵坐标、各记录点的时刻为横坐标作图，得到将蛋白质从DNA表面匀速拉开过程中力随时间的变化曲线，曲线中纵坐标的峰值是蛋白质与DNA之间的最大相互作用力。

上述的利用计算机模拟蛋白质与DNA之间相互作用力的方法中：

步骤(4)中所述的温度控制优选为310K；

步骤(7)中所述的模拟时间优选为4纳秒；

步骤(8)中所述的SMD系数优选为0.7kcal/mol·²。

利用本发明所述的方法进行蛋白质与DNA之间相互作用力的计算与传统实验方法相比具有明显的优越性：

(1)测定对象范围广泛，在数据库中依靠非共价键相互作用结合形成的开放式蛋白质-DNA复合物之间的力都能测定；

(2)实验材料易于制备，利用适当软件对蛋白质-DNA复合物的晶体结构稍加处理即可；

(3)设备要求低，在普通PC机或由普通PC机组成的机群(Cluster)上计算即可，不需要购买原子力显微镜等大型设备；

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

(5)便于在与分子识别和基因表达行为相关的生命科学和物理化学以及医药领域广泛应用。

附图说明

图1POU_HD结构域和与之相互作用的DNA片断溶解在0.9％NaCl溶液中。

图中POU_HD结构域和DNA片断位于水箱的右侧，左侧有足够的空间供SMD模拟中POU_HD结构域运动。其中：A是POU_HD结构域，B是与POU_HD结构域相互作用的DNA片断。

图2SMD模拟后，POU_HD结构域远离DNA片断。

图中POU_HD结构域沿DNA-POU_HD质心连线方向运动到水箱左侧，DNA位置不变。其中：A是POU_HD结构域，B是与POU_HD结构域相互作用的DNA片断。

图3POU_HD结构域运动过程中受力的变化情况。

在100ps时POU_HD结构域受到的力最大，约为1542pN。

具体实施方式

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

实施例1：

(1)从PDB数据库(http://www.rcsb.org/pdb/home/home.do)获得代码为1o4x的复合物数据文件，命名为1o4x.pdb；

(2)利用VMD-1.8.6软件提取欲研究的目标蛋白质和DNA片断，将其溶解到盛有0.9％NaCl溶液的长方体水箱(84*70*61)中，得到蛋白质-DNA-盐溶液体系(见图1)；所述体系中含有非氢原子空间位置数据的结构，命名为1o4x_ionized.pdb，含有氢原子以及氢键参数的结构，命名为1o4x_ionized.psf；

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

(4)用20kcal/mol·²的和谐约束作用(Harmonic Constraint)约束蛋白质的α碳原子(CA atoms)，并用10kcal/mol·²的和谐约束作用约束DNA骨架原子(nucleic backboneatoms)，采用与步骤(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，并且使用Langevin动力学(Langevin Dynamics)控温方法和Nose-Hoover Langevin piston控压方法将温度控制在310K，压强控制在1atm，对步骤(3)能量最小化处理后的体系加热100皮秒(ps)，使体系逐渐升温并稳定到310K；

(5)对步骤(4)加热后的体系，用10kcal/mol·²的和谐约束作用约束蛋白质的α碳原子，其它条件与步骤(4)相同，进行100皮秒(ps)的分子动力学模拟；

(6)去除对蛋白质的α碳原子的约束，其它条件同步骤(4)，对步骤(5)模拟后的体系进行100皮秒(ps)的分子动力学模拟；

(7)用10kcal/mol·²的和谐约束作用约束DNA片段两个末端各一个碱基对的所有原子，固定体系中所有的氢氧键长，将PME(Particle Mesh Ewald)格子边长设置为1±0.1，采用与步骤(4)相同的力场、周期性边界条件、范德华力和电场力计算方法、控温控压方法、PME长程电场力计算方法，对步骤(6)模拟后的体系进行4纳秒(ns)的分子动力学模拟，使得体系的RMSD(Root Mean Square Distance)值达到水平；

(8)将步骤(7)模拟后体系中的DNA骨架原子固定，以蛋白质的α碳原子为SMD(steered molecular dynamics)原子，将SMD系数(SMDk)设定为0.7kcal/mol·²的任意数值，SMD移动速度(SMDVel)设定为0.0005/fs，SMD记录频率(SMDOutputFreq)设定为100，其它条件与步骤(7)相同，对步骤(7)模拟后的体系进行100皮秒(ps)的SMD模拟(结果见图2)，得到该过程中各记录点的热力学参数，命名为1o4x_smd.log；

(9)用VMD-1.8.6软件处理文件1o4x_smd.log，得到各记录点的力，命名为1o4x_ft.dat；

(10)以各记录点的力为纵坐标、各记录点的时刻为横坐标作图，得到将蛋白质从DNA表面匀速拉开过程中力随时间的变化曲线(见图3)，曲线中纵坐标的峰值是蛋白质与DNA之间的最大相互作用力，约为1542pN。

实施例2：

(1)从PDB数据库(http://www.rcsb.org/pdb/home/home.do)获得代码为1wtq的复合物数据文件，命名为1wtq.pdb；

(2)利用VMD-1.8.6软件提取欲研究的目标蛋白质和DNA片断，将其溶解到盛有0.9％NaCl溶液的长方体水箱中，得到蛋白质-DNA-盐溶液体系；所述体系中含有非氢原子空间位置数据的结构，命名为1wtq_ionized.pdb，含有氢原子以及氢键参数的结构，命名为1wtq_ionized.psf；

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

(4)用20kcal/mol·²的和谐约束作用(Harmonic Constraint)约束蛋白质的α碳原子(CA atoms)，并用10kcal/mol·²的和谐约束作用约束DNA骨架原子(nucleic backboneatoms)，采用与步骤(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，并且使用Langevin动力学(Langevin Dynamics)控温方法和Nose-Hoover Langevin piston控压方法将温度控制在280K，压强控制在1atm，对步骤(3)能量最小化处理后的体系加热100皮秒(ps)，使体系逐渐升温并稳定到280K；

(5)对步骤(4)加热后的体系，用10kcal/mol·²的和谐约束作用约束蛋白质的α碳原子，其它条件与步骤(4)相同，进行100皮秒(ps)的分子动力学模拟；

(6)去除对蛋白质的α碳原子的约束，其它条件同步骤(4)，对步骤(5)模拟后的体系进行100皮秒(ps)的分子动力学模拟；

(7)用10kcal/mol·²的和谐约束作用约束DNA片段两个末端各一个碱基对的所有原子，固定体系中所有的氢氧键长，将PME(Particle Mesh Ewald)格子边长设置为1±0.1，采用与步骤(4)相同的力场、周期性边界条件、范德华力和电场力计算方法、控温控压方法、PME长程电场力计算方法，对步骤(6)模拟后的体系进行2纳秒(ns)的分子动力学模拟，使得体系的RMSD(Root Mean Square Distance)值达到水平；

(8)将步骤(7)模拟后体系中的DNA骨架原子固定，以蛋白质的α碳原子为SMD(steered molecular dynamics)原子，将SMD系数(SMDk)设定为0.1kcal/mol·²的任意数值，SMD移动速度(SMDVel)设定为0.0005/fs，SMD记录频率(SMDOutputFreq)设定为100，其它条件与步骤(7)相同，对步骤(7)模拟后的体系进行100皮秒(ps)的SMD模拟，得到该过程中各记录点的热力学参数，命名为1wtq_smd.log；

(9)用VMD-1.8.6软件处理文件1wtq_smd.log，得到各记录点的力，命名为1wtq_ft.dat；

(10)以各记录点的力为纵坐标、各记录点的时刻为横坐标作图，得到将蛋白质从DNA表面匀速拉开过程中力随时间的变化曲线，曲线中纵坐标的峰值是蛋白质与DNA之间的最大相互作用力。

实施例3：

(1)从PDB数据库(http://www.rcsb.org/pdb/home/home.do)获得代码为1hry的复合物数据文件，命名为1hry.pdb；

(2)利用VMD-1.8.6软件提取欲研究的目标蛋白质和DNA片断，将其溶解到盛有0.9％NaCl溶液的长方体水箱中，得到蛋白质-DNA-盐溶液体系；所述体系中含有非氢原子空间位置数据的结构，命名为1wtq_ionized.pdb，含有氢原子以及氢键参数的结构，命名为1wtq_ionized.psf；

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

(4)用20kcal/mol·²的和谐约束作用(Harmonic Constraint)约束蛋白质的α碳原子(CA atoms)，并用10kcal/mol·²的和谐约束作用约束DNA骨架原子(nucleic backboneatoms)，采用与步骤(3)中相同的力场、周期性边界条件、范德华力和电场力计算方法，并且使用Langevin动力学(Langevin Dynamics)控温方法和Nose-Hoover Langevin piston控压方法将温度控制在320K，压强控制在1atm，对步骤(3)能量最小化处理后的体系加热100皮秒(ps)，使体系逐渐升温并稳定到320K；

(5)对步骤(4)加热后的体系，用10kcal/mol·²的和谐约束作用约束蛋白质的α碳原子，其它条件与步骤(4)相同，进行100皮秒(ps)的分子动力学模拟；

(6)去除对蛋白质的α碳原子的约束，其它条件同步骤(4)，对步骤(5)模拟后的体系进行100皮秒(ps)的分子动力学模拟；

(7)用10kcal/mol·²的和谐约束作用约束DNA片段两个末端各一个碱基对的所有原子，固定体系中所有的氢氧键长，将PME(Particle Mesh Ewald)格子边长设置为1±0.1，采用与步骤(4)相同的力场、周期性边界条件、范德华力和电场力计算方法、控温控压方法、PME长程电场力计算方法，对步骤(6)模拟后的体系进行6纳秒(ns)的分子动力学模拟，使得体系的RMSD(Root Mean Square Distance)值达到水平；

(8)将步骤(7)模拟后体系中的DNA骨架原子固定，以蛋白质的α碳原子为SMD(steered molecular dynamics)原子，将SMD系数(SMDk)设定为1kcal/mol·²的任意数值，SMD移动速度(SMDVel)设定为0.0005/fs，SMD记录频率(SMDOutputFreq)设定为100，其它条件与步骤(7)相同，对步骤(7)模拟后的体系进行100皮秒(ps)的SMD模拟，得到该过程中各记录点的热力学参数，命名为1hry_smd.log；

(9)用VMD-1.8.6软件处理文件1hry_smd.log，得到各记录点的力，命名为1hry_ft.dat；

(10)以各记录点的力为纵坐标、各记录点的时刻为横坐标作图，得到将蛋白质从DNA表面匀速拉开过程中力随时间的变化曲线，曲线中纵坐标的峰值是蛋白质与DNA之间的最大相互作用力。