基于随机模拟的折叠病致病机理的分析方法

摘要

本发明公开了一种基于随机模拟的折叠病致病机理的分析方法。该方法应用ABEEMσπ浮动电荷力场能量模型进行研究，模拟及计算蛋白质系统的态密度的主从模式采用MPI并行程序算法进行计算；所述主从模式是指在N个分进程中，其中分进程1为主进程，其余都为子进程。本发明采用主从模式的MPI并行程序算法，有效加快了模拟和计算蛋白质系统态密度的速度；在计算过程中，通过引入具有适应特点的更新修正因子，进一步提高了计算的精度和速度，可应用于研究蛋白质折叠热力学过程及其致病机理。

说明书

【技术领域】

本发明属于生物信息学、生物化学和随机模拟技术领域，特别涉及一种基 于随机模拟的折叠病致病机理的分析方法。

【背景技术】

除了编码基因组和确定蛋白质中氨基酸的序列，生物信息学和生物化学的 另一个越来越重要的任务是预测蛋白质结构和功能。这个任务关键需要理解蛋 白质折叠成其天然结构的作用机制，也即蛋白质折叠研究。折叠病是一类因蛋 白质构象发生改变所引起疾病。由于蛋白质的功能与其三维结构是密切相关的， 故蛋白质折叠研究对探索各种折叠病的致病机理意义重大。

WangLandau算法就是随机模拟领域最吸引人最有发展情景的新算法，它能 解决生物信息学、生物化学等多个领域的很多复杂的问题。比如在蛋白质折叠 研究，该算法有两个最显著的优点：第一，蛋白质模拟不会局限在局部最小能 量状态，因而能较好的在整个能量区间进行自由行走；第二，通过该算法可模 拟和计算出蛋白质系统态密度，因而就能进一步求解得到宽广温度范围内的很 多热动力学量如比热等，这样就能高效的分析和研究蛋白质折叠的整个热力学 过程。

【发明内容】

本发明的目的在于提供一种基于随机模拟的折叠病致病机理的分析方法。

本发明的目的通过以下技术方案实现：一种基于随机模拟的折叠病致病机 理的分析方法，包括如下步骤：

（1）确定蛋白质能量模型和能量区间：

应用ABEEMσπ浮动电荷力场能量模型进行研究，然后使用全局优化算法获 得蛋白质能量区间；在ABEEMσπ模型中，分子总能量表示为：

E_ABEEMσπ＝E_b+E_θ+E_φ+E_imptors+E_vdW+E_elec

其中，E_b为键伸缩振动势能；E_θ为键角弯曲振动势能；E_φ和Ei_mptors分别为二 面角扭转势能和非共面扭转势能；E_vdW、E_elec均为非键作用势能；

（2）确定蛋白质能量区间的分段方式：对步骤（1）的蛋白质能量区间进 行离散化处理，若取k个能量bin区间值，则对[E_min,E_max]平均划分k个bin区间， 用每个bin区间中间的一个能量值代表能量区间值；

（3）模拟及计算蛋白质系统态密度：

所述的模拟及计算蛋白质系统的态密度的主从模式采用MPI并行程序算法 进行计算；所述主从模式是指在N个分进程中，其中分进程1为主进程，其余 都为子进程。

步骤（1）中：

所述ABEEMσπ模型中，将分子体系划分为原子区域、化学键区域及孤对电 子区域；孤对电子位点处于距离原子共价半径处；

所述化学键区域分为σ键和π键区域；其中σ键位点处于两个成键原子共价 半径之比处；π键位点处于垂直于双键所在的平面上，且置于双键原子上下两侧 共价半径处；

所述分子总能量通过如下公式进行计算：

$E_{\mathrm{ABEEM\sigma \pi }}\left(\begin{array}{c}=\underset{\mathrm{Bonds}}{\Sigma }{k}_r{(r-r_{\mathrm{eq}})}^2+\underset{\mathrm{angles}}{\Sigma }{k}_{\theta }{(\theta -{\theta }_{\mathrm{eq}})}^2+\underset{\mathrm{imptors}}{\Sigma }v[1-\cos \left(2\phi \right)]\\ +\underset{\mathrm{torsions}}{\Sigma }\left\{\frac{V_1}{2}\right[1+\cos \left(\phi \right)]+\frac{V_2}{2}[1-\cos \left(2\phi \right)]+\frac{V_3}{2}[1+\cos \left(3\phi \right)\left]\right\}\\ +\underset{i<j}{\Sigma 4f_{\mathrm{ij}}{ϵ}_{\mathrm{ij}}({\sigma }_{\mathrm{ij}}^{12}/r_{\mathrm{ij}}^{12}-{\sigma }_{\mathrm{ij}}^6/r_{\mathrm{ij}}^6)+\underset{i<j}{\Sigma }}{k}_{\mathrm{ij}}{q}_i{q}_j/r_{\mathrm{ij}}\end{array}\right)$

其中，k_r和k_θ分别为键伸缩和键角弯曲势能的力常数，r和θ分别为实际的 键长和键角值，r_eq和θ_eq分别为平衡键长和平衡键角值，V₁、V₂、V₃和v分别为二 面角扭转势能项及非共面扭转势能项的展开力常数；对于非键相互作用项的范 德华相互作用，采用标准的联合规则，即ε_ij＝(ε_iiε_jj)^1/2,σ_ij＝(σ_iiσ_jj)^1/2；对分子A和B 间或分子A和A内的相互作用，求和遍及所有i＜j的原子对；对于任何1-2和1-3 关系的原子对，f_ij＝0；对于任何1-4关系的原子对，f_ij＝0.5；其余情况下f_ij＝1.0； q_i和q_j分别为位点i和j的电荷数，r_ij是位点i和j之间的距离，当i和j之间的间 隔为两原子以下时，k_ij＝0；当i和j在氢键相互作用区域时，k_ij＝k_H-bond(氢键拟 合函数)；其它所有情况下，k_ij＝0.57。

键长伸缩振动能和键角弯曲振动能项中的参数被认为是“硬自由度”参数， 它们具有较好的可转移性，因此直接采用了OPLS-AA固定电荷力场的相应参数； 其它参数[包括电荷参数(x*,η*)、非共面扭转势能项的力常数、二面角力场数和 范德华参数]都是经ABEEMσπ模型重新调节的。

所述全局优化算法优选为模拟退火算法。

步骤（3）中：

所述主进程的算法采用如下步骤进行：

S1:初始化系统的态密度的对数S(E)＝lng(E)＝0，直方图H(E)＝(， E_min≤E≤E_max；修正因子df＝1(＝lnf＝lne)，总步数计数器s＝1;

S2:t＝1;

S3:在主进程中，对步骤（1）的ABEEMσπ模型进行随机变动，得到新的模 型，计算能量E_new，根据Metropolis准则确定新模型被接受的概率：简称为MCS 步，t＝t+1；

$P(\mathrm{old}\to \mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

若接受新模型，则：

S(E_new)＝S(E_new)+df,H(E_new)＝H(E_new)+1；

否则就：

S(E_old)＝S(E_old)+df,H(E_old)＝H(E_old)+1；

S4:经过tmax次（如100次）MCS后（即S3循环tmax次），所有进程间相互 通信，主进程收集所有从进程的S_tmp(E)和H_tmp(E)并累加计算出全局的S(E)和 H(E)，即S(E)＝S(E)+所有从进程的S_tmp(E)，H(E)＝H(E)+所有从进程的H_tmp(E)， 然后将全局的S(E)和H(E)的广播给所有从进程；

S5:改变修正因子df＝ln(s+tmax)-ln(s)、总步数计数器s＝s+tmax，再返 回到S2继续迭代，直到进程终止条件（也即如可取 0.0001），这样就能求得S(E)进而得到蛋白质系统的相对的态密度g(E)＝eS(E)。

所述从进程的算法采用如下步骤进行：

_S1:初始化系统的态密度的对数S(E)＝lng(E)＝0，S_tmp(E)＝lng_tmp(E)＝0，直方 图H(E)＝0，H_tmp(E)＝0，E_min≤E≤E_max；修正因子df＝1(＝lnf＝lne),总步数计数器 s＝1;

S2:t＝1;

S3:在从进程中，对步骤（1）的ABEEMσπ模型进行随机变动，产生新的模 型，计算能量E_new，根据Metropolis准则确定新模型被接受的概率：简称为MCS 步，t＝t+1；

$P(\mathrm{old}\to \mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

若接受新模型，则：

S(E_new)＝S(E_new)+df,H(E_new)＝H(E_new)+1，

S_tmp(E_new)＝S_tmp(E_new)+df,H_tmp(E_new)＝H_tmp(E_new)+1；

否则就：

S(E_old)＝S(E_old)+df,H(E_old)＝H(E_old)+1，

S_tmp(E_old)＝S_tmp(E_old)+df,H_tmp(E_old)＝H_tmp(E_old)+1；

S4:经过tmax次（如100次）MCS后（即S3循环tmax次），所有进程间相互 通信，从进程将S_tmp(E)和H_tmp(E)发送给主进程，然后接收经主进程计算得到的全 局的S(E)和H(E)更新原来的S(E)和H(E)，再将S_tmp(E)和H_tmp(E)初始化为0；

S5:改变修正因子df＝ln(s+tmax)-ln(s)、总步数计数器s＝s+tmax，再返回到 S2继续迭代，直到进程终止条件（也即可取0.0001）。

本发明相对于现有技术具有如下的优点及有益效果：

本发明采用主从模式的MPI并行程序算法，有效加快了模拟和计算蛋白质 系统态密度的速度；在计算过程中，通过引入具有适应特点的更新修正因子， 进一步提高了计算的精度和速度，可以应用于研究蛋白质折叠热力学过程及其 致病机理。

【附图说明】

图1是实施例1的模拟及计算蛋白质系统的态密度的并行算法的流程详 图。

图2是实施例1的基于随机模拟的折叠病致病机理的分析方法的流程图。

【具体实施方式】

下面结合实施例和附图对本发明作进一步详细的描述，但本发明的保护范 围并不限于此。

实施例1

一种基于随机模拟的折叠病致病机理的分析方法，包括如下步骤：

（1）确定蛋白质能量模型和能量区间：

应用ABEEMσπ浮动电荷力场能量模型进行研究，然后使用全局优化算法获 得蛋白质能量区间；在ABEEMσπ模型中，分子总能量表示为：

E_ABEEMσπ＝E_b+E_θ+E_φ+E_imptors+E_vdW+E_elec

其中，E_b为键伸缩振动势能；E_θ为键角弯曲振动势能；E_φ和E_imptors分别为二 面角扭转势能和非共面扭转势能；E_vdW、E_elec均为非键作用势能；

（2）确定蛋白质能量区间的分段方式：对步骤（1）的蛋白质能量区间进 行离散化处理，若取k个能量bin区间值，则对[E_min,E_max]平均划分k个bin区间， 用每个bin区间中间的一个能量值代表能量区间值；

（3）模拟及计算蛋白质系统态密度：

所述的模拟及计算蛋白质系统的态密度的主从模式采用MPI并行程序算法 进行计算；所述主从模式是指在N个分进程中，其中分进程1为主进程，其余 都为子进程。

步骤（1）中：

所述ABEEMσπ模型中，将分子体系划分为原子区域、化学键区域及孤对电 子区域；孤对电子位点处于距离原子共价半径处；

所述化学键区域分为σ键和π键区域；其中σ键位点处于两个成键原子共价 半径之比处；π键位点处于垂直于双键所在的平面上，且置于双键原子上下两侧 共价半径处；

所述分子总能量通过如下公式进行计算：

$E_{\mathrm{ABEEM\sigma \pi }}\left(\begin{array}{c}=\underset{\mathrm{Bonds}}{\Sigma }{k}_r{(r-r_{\mathrm{eq}})}^2+\underset{\mathrm{angles}}{\Sigma }{k}_{\theta }{(\theta -{\theta }_{\mathrm{eq}})}^2+\underset{\mathrm{imptors}}{\Sigma }v[1-\cos \left(2\phi \right)]\\ +\underset{\mathrm{torsions}}{\Sigma }\left\{\frac{V_1}{2}\right[1+\cos \left(\phi \right)]+\frac{V_2}{2}[1-\cos \left(2\phi \right)]+\frac{V_3}{2}[1+\cos \left(3\phi \right)\left]\right\}\\ +\underset{i<j}{\Sigma 4f_{\mathrm{ij}}{ϵ}_{\mathrm{ij}}({\sigma }_{\mathrm{ij}}^{12}/r_{\mathrm{ij}}^{12}-{\sigma }_{\mathrm{ij}}^6/r_{\mathrm{ij}}^6)+\underset{i<j}{\Sigma }}{k}_{\mathrm{ij}}{q}_i{q}_j/r_{\mathrm{ij}}\end{array}\right)$

其中，k_r和k_θ分别为键伸缩和键角弯曲势能的力常数，r和θ分别为实际的 键长和键角值，r_eq和θ_eq分别为平衡键长和平衡键角值，V₁、V₂、V₃和v分别为二 面角扭转势能项及非共面扭转势能项的展开力常数；对于非键相互作用项的范 德华相互作用，采用标准的联合规则，即ε_ij＝(ε_iiε_jj)^1/2,σ_ij＝(σ_iiσ_jj)^1/2；对分子A和B 间或分子A和A内的相互作用，求和遍及所有i＜j的原子对；对于任何1-2和1-3 关系的原子对，f_ij＝0；对于任何1-4关系的原子对，f_ij＝0.5；其余情况下f_ij＝1.0； q_i和q_j分别为位点i和j的电荷数，r_ij是位点i和j之间的距离，当i和j之间的间 隔为两原子以下时，k_ij＝0；当i和j在氢键相互作用区域时，k_ij＝k_H-bond(氢键拟 合函数)；其它所有情况下，k_ij＝0.57。

键长伸缩振动能和键角弯曲振动能项中的参数被认为是“硬自由度”参数， 它们具有较好的可转移性，因此直接采用了OPLS-AA固定电荷力场的相应参数； 其它参数[包括电荷参数(x*,η*)、非共面扭转势能项的力常数、二面角力场数和 范德华参数]都是经ABEEMσπ模型重新调节的。

所述全局优化算法优选为模拟退火算法。

步骤（3）中：

如图1所示，所述主进程的算法采用如下步骤进行：

S1:初始化系统的态密度的对数S(E)＝lng(E)＝0，直方图H(E)＝（， E_min≤E≤E_max；修正因子df＝1(＝lnf＝lne)，总步数计数器s＝1;

S2:t＝1;

S3:在主进程中，对步骤（1）的ABEEMσπ模型进行随机变动，得到新的模 型，计算能量E_new，根据Metropolis准则确定新模型被接受的概率：简称为MCS 步，t＝t+1；

$P(\mathrm{old}\to \mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

若接受新模型，则：

S(E_new)＝S(E_new)+df,H(E_new)＝H(E_new)+1；

否则就：

S(E_old)＝S(E_old)+df,H(E_old)＝H(E_old)+1；

S4:经过tmax次（如100次）MCS后（即S3循环tmax次），所有进程间相互 通信，主进程收集所有从进程的S_tmp(E)和H_tmp(E)并累加计算出全局的S(E)和 H(E)，即S(E)＝S(E)+所有从进程的S_tmp(E)，H(E)＝H(E)+所有从进程的H_tmp(E)， 然后将全局的S(E)和H(E)的广播给所有从进程；

S5:改变修正因子df＝ln(s+tmax)-ln(s)、总步数计数器s＝s+tmax，再返 回到S2继续迭代，直到进程终止条件（也即如可取 0.0001），这样就能求得S(E)进而得到蛋白质系统的相对的态密度g(E)＝e^S(E)。

如图1所示，所述从进程的算法采用如下步骤进行：

S1:初始化系统的态密度的对数S(E)＝lng(E)＝0，S_tmp(E)＝lng_tmp(E)＝0，直方 图H(E)＝0，H_tmp(E)＝0，E_min≤E≤E_max；修正因子df＝1(＝lnf＝lne),总步数计数器 s＝1;

S2:t＝1;

S3:在从进程中，对步骤（1）的ABEEMσπ模型进行随机变动，产生新的模 型，计算能量E_new，根据Metropolis准则确定新模型被接受的概率：简称为MCS 步，t＝t+1；

$P(\mathrm{old}\to \mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

若接受新模型，则：

S(E_new)＝S(E_new)+df,H(E_new)＝H(E_new)+1，

S_tmp(E_new)＝S_tmp(E_new)+df,H_tmp(E_new)＝H_tmp(E_new)+1；

否则就：

S(E_old)＝S(E_old)+df,H(E_old)＝H(E_old)+1，

S_tmp(E_old)＝S_tmp(E_old)+df,H_tmp(E_old)＝H_tmp(E_old)+1；

S4:经过tmax次（如100次）MCS后（即S3循环tmax次），所有进程间相互 通信，从进程将S_tmp(E)和H_tmp(E)发送给主进程，然后接收经主进程计算得到的全 局的S(E)和H(E)更新原来的S(E)和H(E)，再将S_tmp(E)和H_tmp(E)初始化为0；

S5:改变修正因子df＝ln(s+tmax)-ln(s)、总步数计数器s＝s+tmax，再返回到 S2继续迭代，直到进程终止条件（也即可取0.0001）。

以上所述本发明的具体实施方式，并不构成对本发明保护范围的限定。任 何根据本发明的技术构思所作出的各种其他相应的改变与变形，均应包含在本 发明权利要求的保护范围内。