简化线形柔性大分子动力学的仿真方法.pdf

一种简化线形柔性大分子动力学的仿真方法，从药物分子的药效基团或受体的活性部位入手，建立分子简化线形部分的具有“球－虚弹簧”结构的示意模型，模型中将线形部分一端固定在分子母体上，线形部分由球和虚弹簧相间连接而成，其中的球是将重原子及与之相连的轻原子的质量和转动惯量等效在一个原子上而得到的等效原子，等效原子间连接的虚弹簧同时具备弹簧和棒的特性。根据“球－虚弹簧”模型中原子间非共价键势及原子间作用力，建立柔性线形分子的横向伸缩、轴向扭转振动及线形摆动模型的运动方程，模拟线形生物受体大分子低能构象时存在的形态，可用于计算机辅助药物设计。

1、  一种简化线形柔性大分子动力学的仿真方法，其特征在于包括如下步骤：
1)将药物分子的药效基团或受体的活性部位中近似线形部分简化成线形，建立分子简化线形部分的具有“球-虚弹簧”结构的示意模型，模型中将线形部分一端固定在分子母体上，线形部分由球和虚弹簧相间连接而成，其中的球是将重原子及与之相连的轻原子的质量和转动惯量等效在一个原子上而得到的等效原子，等效原子间用虚弹簧连接，虚弹簧同时具备“弹簧和棒”的特性；
2)分别确定模型中原子间非共价键势能，包括范德华能量项E_vdw，静电作用能量项E_elec，氢键作用能量项E_H；
3)根据非共价键势能计算原子间作用力，根据分子动力学构象仿真分析的方法，用势能函数的梯度_iE来计算力Q_i，随机产生一个初速度，以原子的初始坐标为起点，计算原子在t时刻的新位置和速度。其中，Q_i(t)沿三个方向的分力为：沿轴向的F_i(t)，沿径向的F_j(t)和F_k(t)，F_i(t)，F_j(t)，F_k(t)之间符合右手定律；
4)根据“球-虚弹簧”模型建立伸缩模型的振动方程式，振动方程式用等效原子的质量矩阵、虚弹簧的弹性系数矩阵和轴向力向量来表示，轴向力向量等于质量矩阵与加速度向量的乘积加上弹性系数矩阵与位移向量的积，其中，质量矩阵是对角阵，每个等效原子质量m_i是重原子连同相连的轻原子的质量和，加速度向量是等效原子轴向位移向量的两次导数；
5)根据“球-虚弹簧”模型建立扭转模型的振动方程式，振动方程式用等效原子的转动惯量矩阵、虚弹簧的扭转系数矩阵和扭转力矩向量来表示，扭转力矩向量为转动惯量矩阵与角加速度向量的乘积加上扭转弹性系数矩阵与角位移向量的积，这里扭转力矩向量主要是由作用在与重原子相连的轻原子或基团上的与轴向力方向垂直的力与它们到相应重原子距离乘积，其中：转动惯量矩阵是对角阵，每个对角元素J_i是等效转动惯量，大小等于重原子连同相连轻原子对线形分子轴线的转动惯量和，角加速度向量是等效原子轴向位移向量的两次导数；
6)根据“球-虚弹簧”模型建立摆动模型的摆动方程式，摆动方程式用等效原子的质量矩阵、虚弹簧的弯曲系数矩阵和径向力向量来表示，其中，径向力向量包括作用于每个等效原子上的垂直于线分子方向的分力和固定架作用于原子上的力；
7)分别对伸缩、扭转及摆动三种模型的运动方程式求解，然后用计算机高级语言和三维图形仿真工具对运动表达式仿真，得到线形生物受体大分子低能构象时存在的形态，包括瞬时位置，速度，加速度。

简化线形柔性大分子动力学的仿真方法
技术领域
本发明涉及一种简化线形柔性大分子动力学的仿真方法，通过对线形大分子的活性部位的振动和摆动的仿真，给出了线形大分子在低能状态下存在的形态(位置，速度，加速度)，属于计算机系统仿真技术领域。
背景技术
目前，大分子结构研究已经成为一个热点。人们希望通过对大分子结构的研究为发现新药物找到捷径，这种基于结构的寻找药物的方法就是合理药物设计，其中重要的一点是靶点(即受体)的三维结构仿真问题。目前的技术手段很难捕捉大分子某一瞬间的状态，即使现有的X射线技术和NMR核磁共振也只能检测到分子在液态和结晶时的构象，但受体在与药物作用时构象会变化，所以人们力图采用分子动力学方法来解决这个问题。目前分子动力学的研究方法很多，基本上都是通过改变分子温度来找到分子的低能构象，但从微观结构上研究分子的动力学特性较少。为了得到受体分子的三维构象，研究人员提出了一些方法，大体上可以分为如下几种：
(1)能量最小化方法，它是运用分子力学知识找到分子能量最小时的构象，即药效构象，其缺点是计算量很大而且容易遗漏许多潜在的先导分子；
(2)遗传算法，它是把遗传学引入药物设计的成功例子，它希望通过对初始值进行杂交变异处理从而找到较优化的分子构象，特别适应于传统方法解决不了的复杂和非线性问题，但由于存在初始群体难以确立和经常会出现在一个结果上左右振动无法迅速平衡等缺点，其应用受到限制；
(3)计算显微镜方法，它是基于不同的电子等密度图反映分子结构的不同信息，来勾勒分子的骨架、键和外形特征，但该方法还不能进行体系的能量计算和构型优化。
已有的简化分子模型有多种，如球-棒结构，带状结构等，但它们都是适应某个专门的研究领域。球-棒结构在三维结构分析中应用广泛，但它们被看成刚性分子不适合受力分析特别是分子动力学仿真分析；而带状结构是为了从整体上研究大分子如蛋白质的二级结构而提出的，也不适于分子动力学的仿真分析。
发明内容
本发明的目的在于针对现有线形大分子动力学仿真方法的不足与局限性，提出一种简化线形柔性大分子动力学的仿真方法，能有效模拟线形生物受体大分子低能构象时存在的形态，包括瞬时位置，速度，加速度，而且能避免大量的繁复运算。
为实现这样的目的，本发明着重从药物分子的药效基团或受体的活性部位入手，即从大分子的局部而不是整体进行简化，将分子一端固定在分子母体上，建立大分子活性部位低能构象的伸缩、扭转振动和摆动模型。本发明通过建立“球-虚弹簧”结构模型来研究柔性线形分子的横向伸缩、轴向扭转振动及线形摆动，在“球-虚弹簧”结构中将重原子及与之相连的轻原子的质量、转动惯量等效在球上，这个球即为等效原子，球与球之间相连的虚弹簧具有弹簧和棒的特征。当只考虑等效原子之间的受力时，建立伸缩振动模型；当考虑扭转力对分子构象影响时需要建立扭转振动模型；而分子之间有径向力或发生碰撞时，这时需要建立摆动模型。这样可以模拟线形生物受体大分子低能构象时存在的形态，用于寻找潜在药物。
本发明的方法具体包括如下步骤：
1、将药物分子的药效基团或受体的活性部位中近似线形部分简化成线形，建立分子简化线形部分的具有“球-虚弹簧”结构的示意模型。模型中将线形部分一端固定在分子母体上，线形部分由球和虚弹簧相间连接而成。其中的球是将重原子(如N，C，O等)及与之相连的轻原子(如H)的质量和转动惯量等效在一个原子上而得到的等效原子，等效原子间用虚弹簧连接。虚弹簧同时具备“弹簧和棒”的特性。
本发明将分子局部简化成线形主要考虑到以下几点：分子中存在近似线形局部的大分子，而且链状分子也可近似看成线形的；一般的大分子结构很复杂，分子间动力学难以计算，需要从一个简化的模型出发来研究复杂问题的本质。
线形大分子中原子之间以共价键结合，结合力基本符合经典牛顿力学，故可以简化成弹性体如“弹簧”；考虑到弹簧对模拟原子间伸缩振动很有效但不能象“棒”一样很好的模拟扭转振动，而实际分子中的原子存在于不同的空间位置，为了既能描述它们不同的空间位置，又能符合它们之间的力学原理，本发明构建了“虚弹簧”这一模型，它同时具备“弹簧和棒”的特性。本发明中，不考虑电子和原子核，把原子看成大小不同的柔性球体。
2、分别确定模型中原子间非共价键势能。本发明中原子间非共价键势能E包括以下几个主要部分：范德华能量项E_vdw，静电作用能量项E_elec，氢键作用能量项E_H。
由于原子之间的非共价键力只有在距离很近时才会很明显，本发明忽略了远距离原子间作用力产生的能量。
3、根据步骤2所得的非共价键势能E，计算原子间作用力。根据分子动力学构象仿真分析地方法，用势能函数的梯度_iE来计算力Q_i，随机产生一个初速度，以原子的初始坐标为起点，计算原子在t时刻的新位置和速度。其中，Q_i(t)沿三个方向的分力为：沿轴向的F_i(t)，沿径向的F_j(t)和F_k(t)，F_i(t)，F_j(t)，F_k(t)之间符合右手定律。
在步骤2中忽略了某些能量，故由它计算的力中：药物作用时药物原子与受体分子的原子之间有力的作用，本发明只考虑一个原子对一个原子的作用力而不计它同时对另外原子的作用力。另外当发生碰撞时可以认为在相应原子上作用了力脉冲，与对能量取导得到的力一起对原子的运动起作用。
4、根据步骤1中的“球-虚弹簧”模型建立伸缩模型的振动方程式，振动方程式用等效原子的质量、虚弹簧的弹性系数和轴向力向量来表示。轴向力向量等于质量矩阵与加速度向量的乘积加上弹性系数矩阵与位移向量的积。其中，质量矩阵是对角阵，每个等效原子质量m_i是重原子连同相连的轻原子的质量和，加速度向量是等效原子轴向位移向量的两次导数。
其中虚弹簧的伸缩弹性系数和计算等效原子质量时用到的各个参数可以从相关的化学参考书中查到。由于原子处在不同位置时受力不同，需要分开考虑，然后将分开的振动方程式写成矩阵形式。从弹性系数矩阵可以看出某个等效原子受力与它前后相连的两个弹簧有关。以上是过程力的作用，当考虑碰撞时，只要在相应等效原子所受外力中加入力脉冲即可。
之所以可以用经典动力学来处理分子中原子间的作用是考虑到分子动力学适用于复杂分子及生物大分子的构象分析，而量子力学却适用于与电子运动有关的性质，当温度不是很低时误差不大。故原子之间受共价键力，符合虎克定律。
5、根据“球-虚弹簧”模型建立扭转模型的振动方程式，振动方程式用等效原子的转动惯量矩阵、虚弹簧的扭转系数矩阵和扭转力矩向量来表示。扭转力矩向量为转动惯量矩阵与角加速度向量的乘积加上扭转弹性系数矩阵与角位移向量的积。这里扭转力矩向量主要是由作用在与重原子相连的轻原子或基团上的与轴向力方向垂直的力与它们到相应重原子距离乘积。其中：转动惯量矩阵是对角阵，每个对角元素J_i是等效转动惯量，大小等于重原子连同相连轻原子对线形分子轴线的转动惯量和，角加速度向量是等效原子轴向位移向量的两次导数。
扭转振动的建模与伸缩振动的建模相似，只是一些参数选取不同，其中虚弹簧的扭转弹性系数和计算等效原子转动惯量时用到的各个参数可以从相关的化学参考书中查到。另外当构象发生改变时也会产生振动，当然这相当于母体提供扭转力矩的冲击振动会逐渐衰退的，其实这和伸缩振动中的碰撞类似。
6、根据“球-虚弹簧”模型建立摆动模型的摆动方程式，摆动方程式用等效原子的质量，虚弹簧的弯曲系数和径向力向量来表示。其中，径向力向量包括作用于每个等效原子上的垂直于线分子方向的分力和固定架作用于原子上的力。其中虚弹簧的弯曲系数和计算等效原子质量时用到的各个参数可以从相关的化学参考书中查到。
7、分别对伸缩、扭转及摆动三种模型的运动方程式求解，然后用计算机高级语言和三维图形仿真工具对运动表达式仿真，得到线形生物受体大分子低能构象时存在的形态，包括瞬时位置，速度，加速度。
由于三种运动：轴向伸缩，环轴扭转和径向摆动是互不干涉的，故可以分别仿真或同时仿真而不必考虑互相影响。
本发明考虑到线形大分子只有某局部而不是整体才对药理起作用，针对药物分子的药效基团和或受体的活性部位进行线形简化，节省了大量的计算量，很好的模拟出线形分子局部在低能任何时间和状态时的存在形态。这较别的方法更适合于展现大分子近似线形局部瞬时的三维结构，当然就更适合于基于结构的药物设计的方法，不象以前的方法只能给出静态时的结构，而实际情况是药物分子与受体分子作用时结构都会发生变化，这样有些潜在药物分子在静态时是合适的但相互作用时候就不行，本发明避免这种情况的发生。
附图说明
图1是本发明中横向伸缩振动简化模型示意图。
图2是本发明中轴向扭转振动简化模型示意图。
图3是本发明中摆动简化模型示意图。
图4是将具体参数带入振动方程式得到的第1、4、7等效原子的伸缩振动运动曲线和10^-3倍作用力的变化曲线。
图5是将具体参数带入振动方程式得到的第1、4、7等效原子的扭转振动运动曲线和10倍扭转力矩的变化曲线。
具体实施方式
本发明的线形分子模型是大分子的一部分而不是整体，可以看成是药物分子的药效基团或受体的活性部位，故简化时可以将分子一端固定在分子母体上。目前，大多数的分子动力学计算中，都是考虑整个分子的力学性能和能量，认为蛋白质分子的功能受大分子整体结构的决定。但是研究分子结构主要的目的是模拟疾病的靶点(受体)的活性部位进而找到治疗的药物，较之现有整体研究的方法，本发明重在从大分子的局部入手进行简化，提出了伸缩、扭转振动模型和摆动模型。
以下结合附图及实施例对本发明的技术方案作进一步描述，具体步骤如下：
1、建立分子简化线形部分的具有“球-虚弹簧”结构的示意模型。在这里选取药物分子的一小段近似线形部分对它进行简化然后建立模型，建模时将线形部分一端固定在分子母体上，首先将虚弹簧与母体相连，其余部分由球和虚弹簧相间连接成线形。其中的球是将重原子(如N，C等)及与之相连的轻原子(如H)的质量和转动惯量等效在一个原子上而得到的等效原子，等效原子间用虚弹簧连接，虚弹簧同时具备“弹簧和棒”的特性。
其中用到的参数有：N、C、H的原子量分别为14、12、1，半径分别为0.75埃米、0.91埃米、0.79埃米，原子间距离C-H为109pm、C-C为134pm、C-N为147pm、N-H为101pm。
2、分别确定模型中原子间非共价键势能。本发明中原子间非共价键势能E包括以下几个主要部分：范德华能量项E_vdw，静电作用能量项E_elec，氢键作用能量项E_H。
E vdw = Σ ( A ij r ij 12 - B ij r ij 6 ) ]]>
其中：A_ij是排斥项系数；B_ij是吸引项系数；r_ij是非键原子距离。
E elec = Σ 1 ϵ Q i Q j r ij ]]>
其中：ε是介电常数；Q_i，Q_j是原子电荷；r_ij是原子之间距离。
E H = M r m - N r n ]]>
其中：通常m＝12，n＝10；r是氢键“给体原子”D与“受体原子”A之间的距离；
M，N依赖于成键原子的化学性质。
                        E＝E_vdw+E_elec+E_H
3、根据步骤2所得的非共价键势能E，确定原子间作用力。根据分子动力学构象仿真分析的方法，用势能函数的梯度_iE来计算力，其中，Q_i(t)沿三个方向的分力为：沿轴向的F_i(t)，沿径向的F_j(t)和F_k(t)，F_i(t)、F_j(t)、F_k(t)之间符合右手定律。
其中的轴向分力{F(t)}是{Q_i}在轴向的投影，这里F_i(t)近似等于N_i乘上cos(w₀*t)，它是将算出的轴向力(很复杂所以要转化成易于处理的三角函数表示的力)通过傅立叶变换成不同频率的三角函数力之和时忽略很小的高次阶项而得到的。这里w₀是基频取10¹⁴Hz，由于变化的原子位置导致变化的力，所以驱动力的频率几乎与位置变化频率相同，同时有[N]＝[4 -5 6 -3 6 -6 -8]^T×10^-8(N)。
将径向力与其到相应重原子距离相乘得到扭转力矩，同样道理M_i(t)近似等于M_i乘以sin(w₀*t)，这里w0＝10^-22Hz，写成矩阵形式有[M]＝[4；-5；6；-3；6；-6；-8]*10^-11(N.m)。
4、根据步骤1中的“球-虚弹簧”模型建立伸缩模型的振动方程式。图1是本发明中横向伸缩振动简化模型示意图。为了建立伸缩振动模型，这里只给出等效原子的质量和虚弹簧的伸缩弹性系数，图1所示模型其实是“球-虚弹簧”模型在分析伸缩振动时的具体应用。这里的虚弹簧同时具备弹簧和棒的特征，而在建立轴向伸缩振动模型只用到它的弹簧特性，所以模型中的虚弹簧用弹簧来表示，这样易于看出球间的轴向力符合虎克规律。图1中：mi是第i个等效原子的质量；ri是各个等效原子的相对距离；xi是各个等效原子的相对位移(选取向右为正方向)；Fi、Fj、Fk是各个等效原子受力在坐标系下的分力；ki是第i个虚弹簧的弹性系数；Qi是由_iE计算出来的矢量力。
为了得到矩阵表示的伸缩振动方程式，首先对处于不同位置的等效原子列出运动方程式：
m 1 x · · 1 + ( k 1 + k 2 ) x 1 - k 2 x 2 = F 1 ( t ) ]]>
m i x · · i + ( k i + k i + 1 ) x i - k i x i - 1 - k i + 1 x i + 1 = F i ( t ) ]]>
m n x · · n + k n x n - k n x n - 1 = F n ( t ) ]]>
将以上的微分方程组用矩阵表示为： [ M ] { x · · ( t ) } + [ K ] { x ( t ) } = { F ( t ) } ]]>
在伸缩振动中需要的参数有{F(t)}，[M]，[K]。
[M]是等效原子质量矩阵。线形大分子的链状结构是由重原子(N，C等)连同轻原子(如H)组成，计算得到矩阵形式的等效质量矩阵[M]＝[20 21 22 23 24 25 26]^T×10^-27(kg)。
[K]是等效原子之间的虚弹簧的弹性系数矩阵。由于不同重原子之间的共价键的弹性系数是不同的，即使相同重原子之间的键也有几种连接形式，如单键、双键、三键，它们的弹性系数也不同，所以弹性系数矩阵中各元素值是不同的，其表达式及带入具体值时如下：
[ K ] = 30 - 12 - 12 28 - 16 - 16 30 - 14 - 14 26 - 12 - 12 20 - 8 - 8 14 - 6 - 6 6 × 100 ]]>
3)当有碰撞存在时，只要在轴向力中加入冲击力脉冲即可(本实例没有)。
5、和步骤4雷同，根据“球-虚弹簧”模型建立扭转模型的振动方程式。图2是本发明中轴向扭转振动简化模型示意图，这个模型其实是“球-虚弹簧”模型在分析扭转振动时的具体应用。为了建立扭转振动模型，这里只给出等效原子的转动惯量和虚弹簧的扭转弹性系数。如前面提到，这里的虚弹簧同时具备弹簧和棒的特征，而在建立扭转振动模型只用到它棒的特性，所以模型中的虚弹簧用棒来表示。图2中：M_i为等效扭转力矩；J_i是等效转动惯量，是重原子和相连轻原子对线形分子轴线的转动惯量之和；k_i是扭转弹性系数。
振动方程式用等效原子的转动惯量，虚弹簧的扭转系数和扭转力矩向量来表示。这里扭转力矩向量主要是由作用在与重原子相连的轻原子或基团上的与轴向力方向垂直的力与它们到相应重原子距离乘积。其中：转动惯量矩阵是对角阵，每个对角元素J_i是等效转动惯量，大小等于重原子连同相连轻原子对线形分子轴线的转动惯量和，扭转力矩向量是作用在每个等效原子上力矩也包括冲击力矩，加速度向量和扭转弹性系数矩阵与步骤4中相似，只是用扭转角代替伸缩位移。同理得到扭转振动矩阵形式的微分方程式： [ J ] { x · · ( t ) } + [ K ] { x ( t ) } = { M ( t ) } ]]>这里[J]＝diag([13，4，7，3，9，15，7])*10^-45(Kg.m²)；而[K]值如下：
[ K ] = 6 - 3 - 3 7 - 4 - 4 10 - 6 - 6 9.5 - 3.5 - 3.5 11.5 - 8 - 8 15 - 7 - 7 7 ( N . m ) ]]>
6、同上，根据“球-虚弹簧”模型建立扭转模型的振动方程式。图3是本发明中摆动简化模型示意图，这个模型其实是“球-虚弹簧”模型在分析扭转振动时的具体应用。为了建立扭转振动模型，这里只给出等效原子的质量和虚弹簧的弯曲系数。图3中：m_i是各个等效原子质量，计算同上；F_ji是第i个等效原子上受到的力F_j(这里考虑的是原子在力F_j作用下的摆动，同理也可以仿照它得到原子在F_k作用下的摆动)。
7、分别对三种模型的运动方程式求解，然后用计算机高级语言和三维图形处理工具对运动表达式仿真，得到线形生物受体大分子低能构象时存在的形态。
这里只给出轴向伸缩和环轴向扭转振动的分立的运动曲线，如图4和图5。
当仅考虑7个等效原子时，通过对伸缩振动矩阵方程的求解得到其运动轨迹，其曲线如图4所示，图4中给出1、4、7三个处于不同位置的等效原子的运动曲线，可以看出第一个振动幅值很小大概只有原子间距离的十分之一，但第4、7个就累积很大了，更何况分子链很长其积累效应就不能忽视了，而且变化规律近似三角函数。而对扭转振动曲线，由图5可以看出其扭转振动幅度在10^-10米数量级，同样它的累积也会很大，但与分子力学研究中的扭转角度相比却为小。