1、分子模拟与药物设计,药物设计,研究蛋白质与小分子药物之间的相互作用,药物设计,基于配体的药物设计,QSAR with CoMFA可以建立分子性质(包括生物活性)与它们结构之间相关的统计学与图形化模型。这些模型可用于预测新化合物的性质或活性。,定量构效关系(QSARs)将分子的化学性质或生物活性与它的结构相关以便设计出活性更好的新化合物,凝血酶抑制剂的CoMSIA分子场等势图(左:立体场,绿色有利黄色不利;右:疏水场,蓝色有利红色不利),药物设计,基于配体的药物设计,叠合一系列分子从而共享出体现生物活性的公共模式,生长出相应的药效团假设,采用先进的遗传算法和多元化的打分函数,考虑分子能量、立体相
2、似性和药效位点重合,同时还考虑构象的柔性,不确定的立体化学性质,可变的环结构,药物设计,基于受体的药物设计,高通量虚拟筛选 大量化合物数据库的筛选是一项昂贵费时的任务。计算高通量筛选可以增加筛选数据集中合适化合物的占有率并且能降低先导物发现的成本,打分函数是一个蛋白-配体原子表面距离的非线性函数的线性组合。蛋白-配体相互作用包括立体作用、极性作用、熵和溶剂化作用。,药物设计,基于受体的药物设计,同源模建是一个根据模板蛋白将一级序列转成3D结构的技术总称,包括了threading和homology两种,人类Xa因子的晶体结构(PDBID:2BOK,橙色)与同源模建结构(白色)的比较,分子动力学G
3、romacs介绍,主要内容:,Gromacs简介 Gromacs原理 Gromacs运算步骤 Gromacs结果分析,GROMACS,Fast Flexible Free,Lysozyme in Water,RMSD,什么是分子模拟,分子模拟是在分子模型的基础上用计算机做实验,“计算机实验”通过模拟微观粒子的运动来计算宏观性质,温度压力黏度传递性质表面张力,分子间的作用模型,牛顿力学量子力学统计力学等,分子模拟的双重性质,分子模拟具有理论和实验的双重性质,分子模拟不能完全取代实验,理论,实验,模拟,理论的正确性,模拟参数的正确性,模拟方法的选择理论的更新,模拟分子在特定环境下一定时间内的构象能
4、量变化趋势;分子动力学模拟是根据分子力学性质建立的适用于生化体系,聚合物,金属或非金属材料的力场和牛顿的运动力学原理发展出来的计算方法;分子力场(force field)中的各项参数,包括键长,键角,电 荷分布等都可以通过量子化学计算得到;与量化计算相比,分子力学在速度方面具有明显的优势;最早出现在上世纪70年代;,分子动力学模拟的基本原理,势能模型,分子动力学对势能函数的依赖性:所有从分子动力学计算出来得到的宏观性质最终都取决于势能模型,分子动力学的核心:牛顿运动方程,势能(位能)模型:,i=1,2,3,N,简单分子的势能模型,r,U,r,例:甲烷,某些惰性气体,质点处理,U,r,方阱模型,
5、U,r,阶梯模型,复杂分子的势能模型,键的振动,键角,扭矩,分子内部各原子(基团)之间的范德华力、静电力一般要计算1-4(相隔超过两个键的原子或基团对),1,5,4,3,2,复杂分子的势能模型,q,q,q,分子之间的范德华力,分子之间的静电力,例子:丙烷,C,C,C,H,H,H,H,H,H,H,H,10根键长作用18个键角作用8个扭矩作用27个范德华力作用27个静电作用,分子动力学程序的一般步骤,初始化,能量优化,平衡,数据产出,避免局部分子重叠,并不是动力学模拟,根据所有分子的当前坐标计算个分子的受力(位能函数)根据受力更新分子的坐标在此过程中收集用来计算宏观性质的有关信息,读入模型参数,模
6、拟控制参数,初始能量优化方法,去除某些可能存在的原子重叠去除某些严重扭曲的键长、键角、扭矩等方法最速下降法牛顿拉夫森方法其他一般优化几千到几万步,分子模拟的体系分类和方法,简单小型体系大型(复杂)体系和并行计算,简单小型体系,气体的模拟小分子体系,不需要复杂的势能模型几百到几千个分子,分子分布稀疏,大部分是短程作用一般用一台微机就可以处理,计算时间几分钟几小时 简单的液体,不涉及太多的界面性质小分子体系,势能模型不是很复杂几百个分子,可能涉及到静电作用,可能需要长程校正用微机也可以处理,计算时间一般几小时几天,大型(复杂)体系和并行算法,必要性体系越来越大模拟时间越来越长 解决办法制造更快的处
7、理器并行计算机,例子:50000原子的生物体系,1ns模拟单个处理器:12天16个并行处理器:1天,或者,MPI,Message Passing Interface90年代初制定和完善的一套并行语法支持Fortran, C, C+简单易学,几种常见的分子动力学软件,NAMDAMBERCHARMMTINKERLAMMPSDL-POLYGROMACS,GROMACS,主要针对生物体系,也适当照顾一般化学体系优点算法好,计算效率高界面友好维护服务好免费软件 缺点兼容性不好,http:/www.gromacs.org/,Gromacs力场,Select the Force Field: From /u
8、sr/local/gromacs/share/gromacs/top: 1: AMBER03 force field (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003) 2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995) 3: AMBER96 force field (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996) 4: AMBER99 force field (Wang et al., J. Comp. Che
9、m. 21, 1049-1074, 2000) 5: AMBER99SB force field (Hornak et al., Proteins 65, 712-725, 2006) 6: AMBER99SB-ILDN force field (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010) 7: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002) 8: CHARMM27 all-atom force field (with CMAP) - versi
10、on 2.0,Gromacs立场,9: GROMOS96 43a1 force field10: GROMOS96 43a2 force field (improved alkane dihedrals) 11: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205) 12: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)13: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656) 14: OPLS-AA/L all-atom forc
11、e field (2001 aminoacid dihedrals)15: DEPRECATED Encad all-atom force field, using full solvent charges16: DEPRECATED Encad all-atom force field, using scaled-down vacuum charges17: DEPRECATED Gromacs force field (see manual) 18: DEPRECATED Gromacs force field with hydrogens for NMR,GROMACS运算流程,pdb2
12、gmx -f 1AKI.pdb -o 1AKI_processed.gro -water spce-f 蛋白质名称 -o坐标 water 水分子类型,Step One: Prepare the Topology,GROMACS运算流程,#include “oplsaa.ff/forcefield.itp“ 立场文件; Name nrexcl 包含的分子名称Protein_A 3 atoms 每个原子的坐标、电荷、质量、类型信息; nr type resnr residue atom cgnr charge mass typeB chargeB massB; residue 1 LYS rtp
13、LYSH q +2.0 1 opls_287 1 LYS N 1 -0.3 14.0067 ; qtot -0.3 2 opls_290 1 LYS H1 1 0.33 1.008 ; qtot 0.03 3 opls_290 1 LYS H2 1 0.33 1.008 ; qtot 0.36 4 opls_290 1 LYS H3 1 0.33 1.008 ; qtot 0.69 5 opls_293B 1 LYS CA 1 0.25 12.011 ; qtot 0.94,Step Two: Examine the Topology,GROMACS运算流程,Step Three: Defin
14、ing the Unit Cell & Adding Solvent,1. Define the box dimensions using editconf.2. Fill the box with water using genbox.editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic-c 居中 d 盒子边缘离质心的距离 bt 盒子类型genbox -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top-cs 水分子类型,GROMA
15、CS运算流程,Step Four: Adding Ions,grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr-f 输入mdp文件 o生成tpr文件mdp 即molecule dynamics parameter文件,起控制动力学模拟过程的作用genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -nn 8-pname 阳离子类型 nname 阴离子类型 nn 离子数目,GROMACS运算流程,Step Five: Energy M
16、inimization,grompp -f minim.mdp -c 1AKI_solv_ions.gro -p topol.top -o em.tprmdrun -v -deffnm em-v 输出每步的运算信息 deffnm 定义输出文件的名称em.log: EM过程的记录文件em.edr: 二进制能量文件em.trr: 二进制轨迹文件em.gro: 能量最小化后的结构,GROMACS运算流程,Step Five: Energy Minimization,能量分析:g_energy -f em.edr -o potential.xvg,GROMACS运算流程,Step Six: Equil
17、ibration,目的:使蛋白质分子周围的溶剂和离子达到平衡状态,防止整个体系坍塌;先等温使得溶剂分子达到既定的模拟温度;再等压使得整个体系的密度达到均衡一致;grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr(等温恒容过程)mdrun -deffnm nvtg_energy -f nvt.edr,GROMACS运算流程,Step Six: Equilibration,grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -o npt.tpr(等温恒压过程)mdrun -deffnm npt
18、g_energy -f npt.edr -o pressure.xvg(分析压力)g_energy -f npt.edr -o density.xvg(分析密度),GROMACS运算流程,Step Eight: Production MD,grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tprmdrun -deffnm md_0_1mdrun -deffnm md_0_1mpirun -np X mdrun_mpi -deffnm md_0_1(并行),GROMACS运算流程,Step Nine: Analysis,t
19、rjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -ur compact-pbc周期性边界条件 -xtc文件压缩后的坐标文件g_rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns-tu 时间单位nsg_rms -s em.tpr -f md_0_1_noPBC.xtc -o rmsd_xtal.xvg -tu ns与结晶结构相比整个体系的RMSD值g_gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xv
20、g回转半径计算,代表整个蛋白质体系的紧致程度,评价体系是否稳定的指标。,GROMACS运算流程,Step Nine: Analysis,GROMACS运算流程,Step Nine: Analysis,ngmx f md.trr (or md.xtc) s md.tpr(GROMACS自带图形分析查看模块),GROMACS运算流程,Step Nine: Analysis,make_ndx f xxx.pdb o xxx.ndx(体系分组)g_confrms f1 1OMB.pdb f2 md.gro o fit.pdb(叠合结构并输出RMSD值)g_rmsf -f xxx.xtc s xxx.t
21、pr b e o xxx.xvg ox xxx.pdb(均方根波动)-b开始时间 e结束时间 ox 输出开始时间到结束时间的平均PDB结构g_hbond f xxx.xtc s xxx.tpr num xxx.xvg输出氢键的数目 r0.35nm 角度小于30trjconv f xxx.xtc s xxx.tpr o xxx.pdb dump xxx在特定的时间点输出该时间点的pdb结构文件,GROMACS运算流程,Step Nine: Analysis,g_dist f md.xtc s md.tpr o dist.xvg(距离计算命令)g_angle f md.tpr n angle.nd
22、x o angledist.xvg(角度计算命令) g_helix1. Helix radius (file radius.xvg). -螺旋半径变化情况,理想值0.23 nm2. Rise per residue (file rise.xvg). 每个螺旋上升的距离为0.15 nm3. Total helix length (file len-ahx.xvg). 整个蛋白质中-螺旋的长度4. Number of helical residues (file n-ahx.xvg). 组成-螺旋的残基数目5. RMS deviation from ideal helix, calculated for the Calpha atoms only 与理想状态下蛋白质结构的RMSD值。,GROMACS运算流程,Step Nine: Analysis,如何继续crashed动力学模拟tpbconv f pre.xtc s pre.tpr e pre.edr o restart.tprmdrun s restart.tpr deffnm 如何进行并行计算lambootgrompp np # f md.mdp c b4md.gro p md.top o md.tprmpirun np # mdrun deffnm md-np CPU个数,预测化合物的ADME/T性质,