NAMD入门教程(三)

4 受控分子动力学模拟(Steered Molecular Dynamics)所谓受控分子动力学模拟（Steered Molecular Dynamics，SMD），就是指在进行分子动力学模拟时，人为地给分子中的某个或某几个原子施加一个假想的外力，或者人为地固定某个或某几个原子的位置。从而控制整个分子的行为。我们前面进行的动力学模拟都是通过各种参数设置，尽量逼近蛋白质分子在溶液体系中的真实状态，以研究分子的各种性质和行为。但受控分子动力学模拟却要用一个假想的外力干扰控制生物大分子的行为，这是为什么呢？（需要查一下更多的SMD的应用实例，详细了解SMD究竟可以用于研究什么问题，才能做出回答）在下面的例子中，我们将首先固定泛素分子中一个原子的位置，然后用一个假想的外力牵拉另一个原子。球形的泛素分子会因此而被逐渐拉开，最终成为伸展状态的肽链。进行SMD时，我们需要用已经完成能量最小化和能量平衡，达到稳定状态的蛋白结构。如果用含有扭曲、拉伸、变形构象的原始蛋白结构，我们将无法分清蛋白各部分的运动是由于蛋白内部的形变张力引起的，还是由于SMD实验附加的外力造成的。因此，进行SMD之前必须预先进行一次平衡态分子动力学模拟，获得稳定结构。这里，我们使用2.4节球状水体分子动力学模拟输出的恢复文件（.restart）进行SMD。恢复文件输出时，泛素已在水体中完成了能量最小化和能量平衡，达到了稳定状态。4.1 除去水分子为了节省计算时间，我们在进行本次SMD之前将除去体系中所有的水分子。但读者必须要注意：在真正进行SMD实验的时候决不可以将水分子除去！1、打开VMD，选择FileNew Molecule菜单项，载入common目录下的文件ubq_ws.psf。不要关闭窗口，此时窗口“Load file for:”一项应该显示“0:ubq_ws.psf”。再次单击按钮Browser，找到1-2-sphere目录下的文件ubq_ws_eq.restart.coor，载入该文件。关闭Molecule File Browser窗口。现在我们载入了2.4节球状水体动力学模拟时输出的恢复文件(.restart.coor)。2、选择Extension tk菜单项，打开tk，首先用cd命令改变当前目录至common目录下，然后输入：set selprotein atomselect top protein$selprotein writepdb ubq_ww_eq.pdb这样我们便在common文件夹下新建了一个pdb文件ubq_ww_eq.pdb，储存已经达到平衡态的蛋白质分子，但没有水分子。3、删除当前分子，但不要关闭VMD。4.2 恒速牵拉图 恒速牵拉的图示。图中SMD原子是蓝色的，假原子（dummy atom）是红色的，二者之间有一根假想的弹簧相连。假原子以恒定的速度运动，弹簧发生形变后作用于SMD原子，于是SMD原子开始运动。在这一节我们将使用恒速牵拉进行SMD动力学模拟。进行恒速牵拉时，我们需要首先设定一个假想的原子作为施力原子，这个实际并不存在的原子就叫做假原子（dummy atom）。假原子通过一根假想的弹簧与系统中真实存在的某个原子（称为SMD原子）相连，牵拉时，假原子以恒定的速度运动，因此称这种牵拉方式为恒速牵拉。实际上被牵拉原子（SMD原子）的速度不是恒定的。SMD原子所受拉力由弹簧的形变和弹性系数决定，满足胡克定律：。又有：其中是SMD原子所受的力，是假想弹簧的形变量。它等于假原子的位移减去SMD原子的位移，其中是SMD原子的初始位置，是SMD原子的当前位置。因此，SMD原子的受力可以计算出来，然后根据牛顿定律计算它的运动情况。4.2.1 设定SMD原子和固定原子（Fixed atom）NAMD使用pdb文件中的B因子一栏区别哪些原子是被固定的，哪些原子未被固定：如果某个原子在该栏中的对应值是1（非0值即可），那么该原子被固定，在整个动力学模拟过程中，它的空间坐标不会变化；反之如果为0，那么该原子不受影响。知识链接：PDB文件的格式使用写字板打开common目录下我们刚刚制作的ubq_ww_eq.pdb，可以看到文件的内容分成几栏（column）：在上图中，给我们提供有用信息的共有10栏数据，已在图中标注出：（1）表示所指示为原子（2）该原子序列号（3）IUPAC标准格式的原子名称（4）残基名称（5）残基序列号（6）原子的X坐标（7）原子的Y坐标（8）原子的Z坐标（9）位置（Occupancy）（10）温度因子或B因子（beta factor）其中后两栏在动力学模拟的时候被NAMD用于标记两类特殊的原子：位置（Occupancy）栏用于标记SMD原子，B因子一栏用于标记固定原子（fixed atom）。4.2.1 设定固定原子（fixed atom）和SMD原子下面我们将使用VMD创建我们所需的pdb文件。我们在本例中将使用pdb文件中B因子一栏（即温度因子）指定哪个原子将被固定，使用位置（occupancy）一栏指定哪个原子作为SMD原子被牵拉。1、在VMD中选择FileNew Molecule 菜单项，单击Browse按钮找到common目录下的文件ubq.psf，载入它。不要关闭Molecule File Browser窗口，注意Load file for 一项应当显示1:ubq.psf。再次单击Browse按钮，载入common目录下的ubq_ww_eq.pdb。关闭Molecule File Browse窗口，在VMD图形窗口中应当可以看到没有水分子的泛素分子（图）。图 载入除去水分子后的泛素分子2、在VMD tk中输入以下命令，固定第一个氨基酸残基的碳：set allatoms atomselect top all$allatoms set beta 0set fixedatom atomselect top “resid 1 and name CA”$fixedatom set beta 1下面是对命令的解释： set allatoms atomselect top all 新建了一个变量allatoms，它代表体系中的所有原子。 $allatoms set beta 0 所有原子在pdb文件中的B栏（即beta）全部设为0。 set fixedatom atomselect top “resid 1 and name CA” 新建一个变量fixedatom，代表第一个氨基酸残基的碳。 $fixedatom set beta 1 fixedatom所代表的原子第一个氨基酸残基的碳的B栏被设置为1。这样NAMD会在分子动力学模拟时保持该原子固定。同样地，我们还需要设定哪一个原子将被牵拉。前面提到过，被牵拉的原子就叫做SMD原子。对于这一种原子，NAMD使用pdb文件中的“位置”（occupancy）一栏进行标记。同样地，0代表不受影响的原子，标记为1的原子将被牵拉。设置方法仍然是在tk中输入命令进行。$allatoms set occupancy 0set smdatom atomselect top “resid 76 and name CA”$smdatom set occupancy 1以下是对上述命令的解释： $allatoms set occupancy 0 当前pdb文件中所有原子的“occupancy”一栏设置为0。 set smdatom atomselect top “resid 76 and name CA” 新建了一个变量smdatom，代表第76个氨基酸（即最后一个氨基酸）的碳。 $smdatom set occupancy 1 变量smdatom所代表的原子的occupancy 一栏设定为1，这样NAMD会在分子动力学模拟时将假想的力作用于该原子上。完成上述设定之后，我们需要让VMD输出我们设置好的文件。如果当前目录不是common，用cd命令改变当前目录到common中，然后输入：$allatoms writepdb ubq_ww_eq.ref 回车后会在当前目录生成文件ubq_ww_eq.ref，这就是我们进行SMD所需的文件，它记录了泛素的结构，以及哪一个原子作为固定原子（fixed atom），哪一个作为牵拉原子。不要关闭VMD，将所有窗口最小化。使用写字板打开common目录下的文件ubq_ww_eq.ref，注意打开时“文件类型”依然要选择“所有文件（*.*）”否则将看不到该文件。在文件开头第六行即可找到固定原子第一个氨基酸甲硫氨酸的碳。这个原子的B因子一栏为1.00，而其它原子都为0.00（图）。拖动滚动条向下，一直到倒数第四行，可以看到我们定义的SMD原子最后一个氨基酸甘氨酸的碳。这个原子的Occupancy一栏是1.00，而其它原子都是0.00（图）。图 固定原子（Fixed atom）的设置图 SMD原子的设置需要说明一点：在上面两张图中出现的原子坐标数值（第7、8、9栏数值）可能和读者看到的不一致。这是因为动力学模拟中，能量平衡（equilibration）一步中各个原子的初始速度是随机选择的。我们使用的pdb文件是各个原子能量平衡结束后的末速度，因为初速度是随机的，末速度也不一定相同。4.2.2 设定拉力的方向现在，关闭写字板。我们已经设定了固定原子（fixed atom）和SMD原子。下面我们需要做的是设定拉力的方向。我们定义拉力的方向为沿固定原子（fixed atom）和SMD原子的连线方向（即由固定原子到SMD原子的矢量的方向）。在VMD tk中输入：set smdpos lindex $smdatom get x y z 0set fixedpos lindex $fixedatom get x y z 0vecnorm vecsub $smdpos $fixedpos输入以上命令后，tk中给出三个数值，这三个值就是固定原子（fixed atom）到SMD原子的矢量坐标，分别为x，y，z坐标。记下这三个值（图），我们后面还会用到。然后关闭VMD即可。图 假想的拉力矢量的坐标4.2.3 配置文件现在我们已经得到了文件：ubq_ww_eq.ref，这个文件储存了原子的位置，以及哪个原子是固定原子，哪个原子是SMD原子。下一步我们将制作NAMD配置文件。注意在下面的步骤中一定要避免输入错误，否则很可能导致动力学模拟无法正常进行。我们将对一个样本配置文件进行修改，制作我们所需要的配置文件。在实际工作中，读者也可以修改NAMD教程中提供的样本配置文件，得到自己所需要的文件。1、使用Windows资源管理器获得我们所需的样本配置文件sample.conf。文件在common目录下。将这一文件粘贴到3-1-pullcv目录下,然后改名为ubq_ww_pcv.conf以和其他文件相区分。2、使用写字板打开ubq_ww_pcv.conf。现在我们看到的是一个配置文件的样本（图）。浏览一下可以发现很多参数是用户需要根据实际情况进行设定的。事实上，这个样本配置文件就是NAMD动力学模拟配置文件的基本蓝本，用户可以在此基础上量身定做符合自己需要的配置文件。图 NAMD配置文件样本3、在Job Description 一栏加入:# N- C- Termini Constant Velocity Pulling当然这只是一个注释， 不加也不会影响我们的动力学模拟。4、在Adjustable Parameters 一部分进行如下改变，以符合我们本次动力学模拟的实际条件：structure mypsf.psf structure ./common/ubq.psfcoordinates mypdb.pdb coordinates ./common/ubq_ww_eq.pdboutputName myoutput outputName ubq_ww_pcv以上各项的含义参见2.3.1 配置文件5、在Input一部分进行如下改变，以符合我们本次动力学模拟的实际条件：parameters par_all27_prot_lipid.inp parameters ./common/par_all27_prot_lipid.inp我们没有使用周期性边界条件，因此不需要设置Periodic Boundary Conditions一部分。此外Force-Field Parameters,Integrator Parameters和PME三项也不需改变。6、Constant Temperature Control 一项应该关闭，因为温度控制可能干扰原子的运动：langevin on langevin off7、Constant Pressure Control不需要进行改变，默认是不启用。8、Fixed Atoms Constraint 需要启用。请对前两行作出如下改变：if 0 if 1 fixedAtomsFile myfixedatoms.pdb fixedAtomsFile ./common/ubq ww eq.ref（注：和C语言中判断表达式真假类似，if1表示条件成立，if0表示条件不成立）第三行不需要改变，因为我们是使用B值这一栏设定固定原子的。如果在实际工作中使用了其他栏，就需要在这里作出相应改变以通知NAMD那一栏用于标记固定原子。9、IMD Settings 一项不需要改变，默认是不启用10、在Extra Parameters 一部分输入以下内容：SMD onSMDFile ./common/ubq_ww_eq.refSMDk 7SMDVel 0.005上面的内容的作用是：通知NAMD我们将使用SMD原子。SMD原子和假原子之间的弹簧弹性常数为 7 kcal/mol/2 (1 kcal/mol = 69.479 pN /)。SMDVel一项设定的就是假原子的速度。这里设定为 0.005/timestep。因为我们设定 1timestep = 2fs，因此假原子的速度为 2.5/ps。在实际动力学模拟过程中，假原子将以该恒定速度运动，并通过一个假想的弹簧同SMD原子相连，牵拉SMD原子运动。11、在Extra Parameters部分我们还需要制定拉力的方向。请输入SMDDir nx ny nz注意nx，ny，nz 三个参数就是我们在4.2.2一节求出的矢量的坐标。12、接下来继续输入：SMDOutputFreq 10这个值指定每隔多少步输出一次SMD过程产生的数据13、最后，在Execution Script部分将动力学模拟时间改为20000步，相当于40ps：run 50000 run 20000然后保存文件，关闭写字板。我们已经完成了配置文件的设定。提示请再次注意，我们的设置仅仅是为了使得动力学模拟能够迅速完成，因此采取了许多不恰当的体系设定。比如，具有科学意义的动力学模拟应当至少计算几个ns（几百万步）,牵拉速度应该尽可能慢一些，在0.1/ps左右。并且正如我们提到的，绝不可以忽略水分子。4.2.4 SMD模拟下面我们可以进行动力学模拟了。为了保证模拟正常进行，请确认3-1-pullcv目录下是否有文件ubq_ww_pcv.conf，并确认common目录下有以下文件： ubq.psf ubq_ww_eq.pdb ubq_ww_eq.ref par_all27_prot_lipid.inp下面，打开terminal，使用cd命令改变当前目录到namd-tutorial/namd目录下（即NAMD所在目录），然后输入：namd2 ./3-1-pullcv/ubq_ww_pcv.conf ./3-1-pullcv/ubq_ww_pcv.log实测：CPU：Pentium M740 1.73GHz，内存：512M 可以在10min之内完成。结果的获得与分析将在4.4节讲解。4.3 恒力牵拉下面我们将进行另一种类型的SMD模拟：恒力牵拉。我们仍需要固定一个原子并指定另一个原子被牵拉，但是拉力是恒定的，直接作用于SMD原子上。因此不需要指定假原子以及假想弹簧的弹性系数。4.3.1 设定固定原子（fixed atom）和SMD原子我们在本例中仍将使用pdb文件中B因子一栏（即温度因子）指定哪个原子将被固定，使用位置（occupancy）一栏指定哪个原子作为SMD原子被牵拉，并设定拉力的大小。1、在VMD中选择FileNew Molecule 菜单项，单击Browse按钮找到common目录下的文件ubq.psf，载入它。不要关闭Molecule File Browser窗口，注意Load file for 一项应当显示0:ubq.psf。再次单击Browse按钮，载入common目录下的ubq_ww_eq.pdb。关闭Molecule File Browse窗口，在VMD图形窗口中应当可以看到没有水分子的泛素分子2、在VMD tk中输入以下命令，固定第一个氨基酸残基的碳：set allatoms atomselect top all$allatoms set beta 0set fixedatom atomselect top “resid 1 and name CA”$fixedatom set beta 13、输入以下命令，设定SMD原子和拉力的大小：$allatoms set occupancy 0set smdatom atomselect top “resid 76 and namd CA”$smdatom set occupancy 11.54这样第76个残基的碳的Occupancy一栏被设定为11.54，其他原子都为0。这样不仅设定了SMD原子，而且设定了拉力大小为11.54 kcal/mol/。相当于800pN（8.0 x 10-10N）4、下面设定拉力的方向。$smdatom set x nx$smdatom set y ny$smdatom set z nz注意nx，ny，nz 三个参数就是我们在4.2.2一节求出的矢量的坐标（0.3671,0.3783,0.8498）。设定好之后图形窗口显示一个极度变形的分子。不过没有关系，这是因为VMD将我们设定的拉力坐标按照原子的坐标显示出来了。图 设定好拉力坐标后图形窗口显示的分子5、下面在tk中首先使用cd命令改变当前目录到common目录下，然后输入$allatoms writepdb ubq_ww_eq2.ref 保存我们的设定。然后打开Windows资源管理器，找到common目录下，检查一下是否已经有ubq_ww_eq2.ref文件？如果有了，就可以关闭VMD。4.3.2 配置文件和上一节我们进行的过程一样，我们还需要修改样本配置文件sample.conf，制作我们需要的配置文件。首先打开资源管理器，将common目录下的sample.conf拷贝到3-2-pullcf目录下，并改名为ubq_ww_pcf.conf，然后用写字板打开该文件开始编辑：3、在Job Description 一栏加入:# N- C- Termini Constant Force Pulling当然这只是一个注释，不加也不会影响我们的动力学模拟。4、在Adjustable Parameters 一部分进行如下改变，以符合我们本次动力学模拟的实际条件：structure mypsf.psf structure ./common/ubq.psfcoordinates mypdb.pdb coordinates ./common/ubq_ww_eq.pdboutputName myoutput outputName ubq_ww_pcf以上各项的含义参见2.3.1 配置文件5、在Input一部分进行如下改变，以符合我们本次动力学模拟的实际条件：parameters par_all27_prot_lipid.inp parameters ./common/par_all27_prot_lipid.inp同上一节，我们不需要设置Periodic Boundary Conditions一部分。此外Force-Field Parameters,Integrator Parameters和PME三项也不需改变。6、同上一节，Constant Temperature Control 一项应该关闭，因为温度控制可能干扰原子的运动：langevin on langevin off7、Constant Pressure Control不需要进行改变，默认是不启用。8、Fixed Atoms Constraint 需要启用。请对前两行作出如下改变：if 0 if 1 fixedAtomsFile myfixedatoms.pdb fixedAtomsFile ./common/ubq ww eq2.ref9、IMD Settings 一项不需要改变，默认是不启用10、在Extra Parameters 一部分输入以下内容：constantforce yesconsforcefile ./common/ubq_ww_eq2.ref上面的内容的作用是：通知NAMD我们将对SMD原子施加一个恒定的拉力。拉力的大小就是Occupancy一栏我们设定的值11.54 kcal/mol/。拉力的方向不需要在此指定，因为我们刚才已经记录到了ubq_ww_eq2.ref文件中。11、最后，在Execution Script部分将动力学模拟时间改为20000步，相当于40ps：run 50000 run 20000然后保存文件，关闭写字板。我们已经完成了配置文件的设定。4.3.3 动力学模拟下面我们可以进行动力学模拟了。为了保证模拟正常进行，请确认3-2-pullcf目录下是否有文件ubq_ww_pcf.conf，并确认common目录下有以下文件： ubq.psf ubq_ww_eq.pdb ubq_ww_eq2.ref par_all27_prot_lipid.inp下面，打开terminal，使用cd命令改变当前目录到namd-tutorial/namd目录下（即NAMD所在目录），然后输入：namd2 ./3-2-pullcf/ubq_ww_pcf.conf ./3-2-pullcf/ubq_ww_pcf.log实测：CPU：Pentium M740 1.73GHz，内存：512M 可以在10min之内完成。4.4 结果分析4.4.1 恒速拉伸SMD结果分析我们将在本节中分析恒速拉伸SMD的结果。如果读者没有成功完成动力学模拟，可以使用3-1-pullcv/example-output中的结果文件。4.4.1.1 日志文件首先我们看一下输出的日志文件。打开写字板，找到日志文件3-1-pullcv/ubq_ww_pcv.log打开它。向下拖动滚动条，是否发现比我们第一次动力学模拟的日志文件多了一些内容？读者可以找到1-2-sphere/ubq_ws_eq.log对照一下增加的内容（图）。图 输出文件中新增的内容可以看到，在每次输出ENERGY之后，紧接着增加了以SMD开头的几行文字（图），这是NAMD输出的关于SMD原子（即我们定义的第76氨基酸的碳）的数据。图中（1）记录的是步数，（2）记录的是输出时SMD原子的坐标，（3）记录的是输出时SMD原子所受拉力的大小（x,y,z三个方向上的分力）。4.4.1.2 轨迹动画下面我们将看一看动力学模拟生成的原子运动轨迹，直观地感受一下泛素肽链是如何被拉力拉开的。1、打开VMD，选择File New Molecule 菜单项，单击按钮Browse找到common目录下的ubq.psf，载入该文件；不要关闭Molecule File Browser窗口，注意Load file for 一项应当显示0:ubq.psf。再次单击Browse按钮，载入3-1-pullcv目录下的轨迹文件ubq_ww_pcv.dcd。关闭Molecule File Browse窗口，在VMD图形窗口中应当可以看到泛素分子。2、为了更形象的观察蛋白质分子的变化，我们需要改变一下显示模式。选择GraphicsRepresentation菜单项，打开显示模式设置，将Drawing Method一项由Lines改为Cartoon（图）以便直观地观察到二级结构的变化。然后调节主窗口右下角的速度滑块至中央，单击最右边的播放按钮，就可以看到整个过程的动画了：注意Cartoon模式下二级结构的表示：螺旋为圆柱，片层为片层箭头，loop为线绳。但是因为VMD不会每一帧都将二级结构重新渲染，因此我们看到，在动画的整个过程中肽链都是线绳状（图）。这是因为我们从最后一帧开始渲染。最后一帧时，整个肽链都处于伸展状态，相当于全是loop；因此VMD在渲染其它帧的时候也认为整条肽链都是loop，从而以线绳表示3、单击播放控制滑块，将它向右推到第一帧，这时的蛋白还是以loop形式表示的（图）。我们需要让VMD重新渲染一下以正确地看到二级结构。方法是打开tk，输入命令：mol ssrecalc top图：伸展状态下VMD认为整个肽链都是loop图 将播放控制滑块推到第一帧 图 输入命令mol ssrecalc top之后重新渲染，可以看到二级结构的表示可以看到VMD重新对二级结构进行了渲染（图）。按照此方法可以逐帧地观察肽链被拉开的过程中二级结构的变化, 研究一下在拉伸的过程中哪一种二级结构先消失，哪一种后消失。下面显示的是调整其他参数后渲染出的三帧。图中左下灰色圆球为固定原子，右上角的灰色圆球是SMD原子。知识链接：蛋白质的解折叠途径（unfolding pathway）做到这一步，可能有的读者会问，我们用一个假想的拉力将蛋白质拉开解折叠有什么实际意义呢？ 事实上，蛋白质的解折叠途径是生命体中很重要的一个过程。我们知道，除少数在线粒体和叶绿体中所合成的蛋白质之外，游离核糖体和内置网膜上的核糖体所合成的蛋白质，一般都含有分拣信号（sorting signal），决定它们的最终去向和定位。分拣型号可以引导蛋白质跨膜转运，将蛋白质从细胞质溶质装移到细胞内的不同部位。在跨膜的过程中，蛋白质必须解折叠，以便于跨膜转运。解折叠的过程是蛋白质信号序列或前导序列一端固定在膜上，整个肽链signal）序列。这一信号序列决定了蛋白质的趋向和最终定位（470页）4.4.1.3 牵引力-时间作图下面我们将分析牵引拉力随时间的变化情况，并使用Excel进行作图。1、在tk中使用cd命令改变当前目录到3-1-pullcv，然后输入：mkdir analysis 这样我们在3-1-pullcv目录下新建了一个文件夹analysis，用于储存我们分析时生成的数据文件。2、下面我们将使用一个脚本文件从日志文件中提取NAMD输出文件中关于SMD受力大小的信息。脚本文件就在3-1-pullcv目录下，读者可以输入dir命令查看当前目录下的文件，应该有一个文件为ft.tcl。下面在tk中输入：source ft.tcl这样就调用了脚本文件ft.tcl。调用之后，脚本会要求用户输入三个值：n_x，n_y，n_z。分别输入我们在4.2.2一节计算得到的拉力矢量坐标（0.3671,0.3783,0.8498）（图）。该脚本最后会在3-1-pullcv/analysis目录下新建文件ft.dat储存提取出的数据。图 ft.tcl脚本的调用3、关闭VMD，下面我们将使用