基于随机优化的老年痴呆致病机理的分析方法

摘要

本发明涉及一种基于随机优化的老年痴呆致病机理的分析方法，包括以下步骤：S1、确定能量模型：采用ECEPP能量力场模型以及角度坐标系；S2、搜索最小自由能蛋白质构型；S3、计算蛋白质的态密度。本发明的基于随机优化的老年痴呆致病机理的高效分析方法，与经典的WangLandau算法相比，使用具有全局更新特点的分段函数形式的修正因子能加快搜索和模拟速度，利用退火机制的灵活的更新修正因子方式可提高模拟精度和速度，且该方法的并行方式可大大加快搜索和模拟速度。

说明书

【技术领域】

本发明涉及生物信息学，特别涉及一种基于随机优化的老年痴呆致病机理 的分析方法。

【背景技术】

蛋白质折叠研究蛋白质在短时间内从一级结构（亦即一维多肽链）折叠成 天然三维结构，形成具有生命功能的大分子。生物体的遗传信息（DNA）通过 RNA转录和翻译过程传递给蛋白质，因此蛋白质折叠也被称为第二遗传密码。 通过蛋白质折叠研究可以揭示生命遗传信息的表达和功能传递的奥秘。蛋白质 在从一级结构折叠到天然三维结构的过程中，可能发生误折叠或聚集，其结构 和功能也受到破坏，从而引起‘折叠病’，例如老年痴呆症。

老年痴呆症又叫阿尔茨海默病，是一种中枢神经系统变性病。老年痴呆症 起病隐袭，病程呈慢性进行性，是老年期痴呆常见的一种类型，主要表现为渐 进性记忆障碍、认知功能障碍、人格改变及语言障碍等神经精神症状，严重影 响社交、职业与生活功能。目前，美国500多万老年痴呆症患者每年的花费是 1830亿美元；中国的老年痴呆症患者也有500-600万。蛋白质折叠研究有助于 蛋白质分子药物的设计，因此，蛋白质折叠研究对探索老年痴呆症致病机理以 及对老年痴呆症的预防和治疗将起重大的帮助作用。

对于一个有100个氨基酸的蛋白质，假设每个氨基酸有10个构象，其构象 数目即有10¹⁰⁰；进一步假设搜索每个构象需要10-15秒，按照目前最快的千万 亿次计算机处理速度，要搜索所有的蛋白质结构空间至少需3×10⁷⁷年。因此， 搜索蛋白质全部结构空间不现实，这也对蛋白质结构高效搜索算法提出了非常 高的要求。

蛋白质折叠算法依赖于一个基本假设，即蛋白质的天然结构是蛋白质自由 能最低的结构。传统的分子动力学模拟和蒙特卡洛模拟方法在蛋白质折叠研究 时，低温下的模拟会‘陷’在蛋白质能量的局部最优区间，很难‘跳’出来从 而找到全局最优解。

经典的WangLandau算法属于随机优化领域，WangLandau算法能解决生物 医学、统计物理等多个领域复杂的问题。例如，蛋白质折叠研究，该算法有两 个最显著的优点：第一，该算法不会局限在局部最小能量状态，因而能搜索到 全局最小能量状态；第二，该算法可模拟和计算出蛋白质系统的态密度，因而 可进一步求解得到宽广温度范围内的热动力学量，如比热，实现对蛋白质系统 进行定量分析。

【发明内容】

本发明要解决的技术问题在于提高随机优化的老年痴呆致病机理的高效方 法精度和速度。

为此，本发明提供一种基于随机优化的老年痴呆致病机理的分析方法，包 括以下步骤：

S1、确定能量模型：采用ECEPP能量力场模型以及角度坐标系；

S2、搜索最小自由能蛋白质构型：包括

主进程，

S211、局部最小化方法获得能量上界E_min和能量下界E_max，初始化蛋白质系 统的态密度对数S(E)、修正因子df，df|_E＝(κΘ(E₀-E)+1)lnf，

其中，Θ(E₀-E)为Heaviside分段函数，κ、E₀、f是模型的参数；

S212、进行Metropolis式的随机变动，构建新构型，计算能量E_new、并修改 E_min或E_max为E_new；

S213、重复步骤S212，完成进程间通信，计算全局S(E)；广播全局S(E)至 所有从进程；

S214、执行步骤S212至步骤S213的迭代，E_min三次保持，迭代终止；

从进程，

S221、局部最小化方法获得能量上界E_min和下界E_max，初始化蛋白质系统的 态密度对数S(E)、修正因子df，df|_E＝(κΘ(E₀-E)+1)lnf，

其中，Θ(E₀-E)为Heaviside分段函数，κ、E₀、f是模型的参数；

S222、进行Metropolis式的随机变动，构建新构型，计算能量E_new、并修改 E_min或E_max为E_new；

S223、重复步骤S222，完成进程间通信，接收步骤S213计算的全局S(E)；

S224、执行步骤S222、步骤S223，接收步骤S214迭代终止信息，停机；

S3、计算蛋白质的态密度：

主进程，

S311、初始化蛋白质系统的态密度对数S(E)，直方图H(E)，修正因子df， df|_E＝(κΘ(E₀-E)+1)lnf，

其中，Θ(E₀-E)为Heaviside分段函数，κ、E₀、f是模型的参数；

S312、进行Metropolis式的随机变动，构建新构型，计算能量E_new；

S313、重复步骤S312，完成进程间通信，计算全局S(E)、H(E)；判断直方 图平缓阈值，当满足是执行步骤S314，否则重复执行步骤S312、步骤S313迭 代；

S314、改变df，并执行步骤S312至步骤S313的迭代，直至df小于第二阈 值，获得蛋白质系统的相对的态密度

从进程，

S321、初始化蛋白质系统的态密度对数S(E)，直方图H(E)，修正因子df， df|_E＝(κΘ(E₀-E)+1)lnf，

其中，Θ(E₀-E)为Heaviside分段函数，κ、E₀、f是模型的参数；

S322、进行Metropolis式的随机变动，构建新构型，计算能量E_new；

S323、重复步骤S322，完成进程间通信，接收步骤S313计算的全局S(E)、 H(E)，即更新原S(E)、H(E)；判断直方图平缓阈值，当满足是执行步骤S324， 否则重复执行步骤S322、步骤S323迭代；

S324、改变df，并执行步骤S322至步骤S323的迭代，直至df小于第二阈 值。

本发明的基于随机优化的老年痴呆致病机理的高效分析方法，包括确定能 量模型、搜索自由能最小的蛋白质构型、模拟及计算蛋白质系统的态密度等步 骤，其与经典的WangLandau算法相比，使用具有全局更新特点的分段函数形式 的修正因子能加快搜索和模拟速度，利用退火机制的灵活的更新修正因子方式 可提高模拟精度和速度，且该方法的并行方式可大大加快搜索和模拟速度。

【附图说明】

图1示出本发明分析方法的流程图。

图2示出本发明分析方法的搜索最小自由能蛋白质构型的流程图。

图3示出本发明分析方法的计算蛋白质的态密度的流程图。

【具体实施方式】

下面结合附图和具体实施例对本发明作进一步详细说明。

如图1所示，本发明提供一种基于随机优化的老年痴呆致病机理的分析方 法，包括以下步骤：

步骤S1、确定能量模型：采用ECEPP能量力场模型以及角度坐标系；

所述ECEPP能量力场模型为：

E_ECEPP＝E_C+E_LJ+E_HB+E_Tor

其中，E_C为两原子的电荷之间的库伦作用力；E_LJ为两原子之间的兰纳 -琼斯作用力；E_HB为氢键作用力；E_Tor为两面角旋转作用力(模型参数，与 蛋白质的多肽链结构有关。)。

所述$E_C={\Sigma }_{(i,j)}\frac{{332q}_i{q}_j}{{\mathrm{ϵr}}_{\mathrm{ij}}};$

所述$E_{\mathrm{LJ}}={\Sigma }_{(i,j)}(\frac{A_{\mathrm{ij}}}{r_{\mathrm{ij}}^{12}}-\frac{B_{\mathrm{ij}}}{r_{\mathrm{ij}}^6});$

所述$E_{\mathrm{HB}}={\Sigma }_{(i,j)}(\frac{C_{\mathrm{ij}}}{r_{\mathrm{ij}}^{12}}-\frac{D_{\mathrm{ij}}}{r_{\mathrm{ij}}^{10}});$

所述E_Tor＝∑_lU_l(1±cos(n_lξ_l))；

其中，r_ij为原子i和j之间的距离；ξ_l为第l个两面角，σ_ij＝0。基于角度坐 标系的ECEPP能量力场，计算效率高于基于笛卡尔坐标系的能量力场。同时， 为便于模拟，本发明进一步对能量值E进行适当划分的离散化处理，如n个能 量区间值，对[Emin，Emax]平均划分n个区间，用每个区间中间的一个能量值 代表该能量区间值。

如图2所示，步骤S2、搜索最小自由能蛋白质构型：自由能最小所对应的 构型采用主从模式的MPI并行算法，具体步骤包括:

主进程，

步骤S211、局部最小化方法获得能量上界E_min和能量下界E_max，初始化蛋白 质系统的态密度对数S(E)、修正因子df，df|_E＝(κΘ(E₀-E)+1)lnf，

其中，Θ(E₀-E)为Heaviside分段函数，κ、E₀、f是模型的参数；本发明 实施例中κ＝5，E₀＝-2，f＝e。初始化t=0。其中，κ是与经典的wanglandau算 法中的修正因子lnf调整倍数比；E₀用于能量分段，使得不同的分段区间具有不 同的搜索和模拟速率。

初始化蛋白质系统的态密度对数S(E)为S(E)＝ln g(E)＝0，（E_min≤E≤E_max）。

步骤S212、进行Metropolis式的随机变动，构建新构型，计算能量E_new、并 修改E_min或E_max为E_new；具体为：

若，E_new＜E_min，则，E_new～E_min；初始化态密度对数S(E)为0，并修改E_min为E_new；

若，E_max＜E_new，则，E_max～E_new；初始化态密度对数S(E)为0，并修改E_max为E_new；

并根据Metropolis准则确定新构型被接受的概率（Metropolis式随机游动）：

$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；否则S(E_old)＝S(E_old)+df。对t累加， 即t=t+1。判定：t mod 1000=0?

步骤S213、重复步骤S212，完成进程间通信，计算全局S(E)；主进程和所 有的从进程通信，传递S_tmp(E)，H_tmp(E)信息以便计算、更新得到全局的S(E)，H(E)。

本发明实施中重复次数为100，也可为80或200，具体根据测试数据调整 重复次数；上述计算全局S(E)为：主进程收集所有从进程的S_tmp(E)并累加计算出 全局S(E)，亦即S(E)＝S(E)+所有从进程S_tmp(E)，并对选取E_min和E_max分别为所有 进程中（主进程和从进程）的最小值或最大值。

广播全局S(E)至所有从进程。

步骤S214、执行步骤S212至步骤S213的迭代，E_min三次保持，迭代终止。 具体为重复步骤S212到步骤S213的迭代直到连续三次迭代E_min都保持原值；此 时E_min及其构型E_min自适应并有效获得。

从进程，

步骤S221、局部最小化方法获得能量上界E_min和下界E_max，初始化蛋白质系 统的态密度对数S(E)、修正因子df，df|_E＝(κΘ(E₀-E)+1)lnf，

其中，Θ(E₀-E)为Heaviside分段函数，κ、E₀、f是模型的参数；本发明 实施例中κ＝5，E₀＝-2，f＝e。

初始化蛋白质系统的态密度对数S(E)为S(E)＝ln g(E)＝0，（E_min≤E≤E_max）。

步骤S222、进行Metropolis式的随机变动，构建新构型，计算能量E_new、并 修改E_min或E_max为E_new；具体为：

若，E_new＜E_min，则，E_new～E_min；初始化态密度对数S(E)为0，并修改E_min为E_new；

若，E_max＜E_new，则，E_max～E_new；初始化态密度对数S(E)为0，并修改E_max为E_new；

并根据Metropolis准则确定新构型被接受的概率：

$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，S_tmp(E_new)＝S_tmp(E_new)+df；

否则S(E_old)＝S(E_old)+df，S_tmp(E_old)＝S_tmp(E_old)+df。

判定：t mod 1000=0?

步骤S223、重复步骤S222，完成进程间通信，接收步骤S213计算的全局 S(E)，即更新原S(E)，即用计算的全局S(E)代原来从进程的S(E)。；其还包括在 完成进程间通信后，将从进程的S_tmp(E)发送至主进程。对选取E_min和E_max分别为 所有进程中的最小值或最大值，并对S_tmp(E)初始化为0。

步骤S224、执行步骤S222、步骤S223，接收步骤S214迭代终止信息，停 机；亦即继续步骤S222到步骤S223的迭代直到主进程满足停机条件。

以上为本发明一主进程以及从进程1的流程，其还包括从进程的并行进程， 进程2、。。。进程N。

如图3所示，步骤S3、计算蛋白质的态密度，亦即计算蛋白质系统的态密 度的主从模式的MPI并行程序算法步骤为：

主进程，

步骤S311、初始化蛋白质系统的态密度对数S(E)，直方图H(E)，修正因子 df，df|_E＝(κΘ(E₀-E)+1)lnf，

其中，Θ(E₀-E)为Heaviside分段函数，κ、E₀、f是模型的参数；本发明 实施例中κ＝5，E₀＝-2，f＝e。初始化t=0。

初始化蛋白质系统的态密度对数S(E)为S(E)＝lng(E)＝0，H(E)=0； （E_min≤E≤E_max）。

步骤S312、进行Metropolis式的随机变动，构建新构型，计算能量E_new；

并根据Metropolis准则确定新构型被接受的概率（Metropolis式随机游动）：

$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。对t累加，即t=t+1。判定：t mod 1000=0?

步骤S313、重复步骤S312，完成进程间通信，计算全局S(E)、H(E)；判断 直方图平缓阈值，当满足是执行步骤S314，否则重复执行步骤S312、步骤S313 迭代；

具体为：主进程收集所有从进程的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)至所有从进程，并求得

判断直方图平缓阈值<φ，即

$\frac{\max \left(H_{\mathrm{real}}\left(E\right)\right)-\min \left(H_{\mathrm{real}}\left(E\right)\right)}{\max \left(H_{\mathrm{real}}\left(E\right)\right)+\min \left(H_{\mathrm{real}}\left(E\right)\right)}<\phi (0<\phi <1),$此时执行步骤S314，本发明实施例 中φ可取0.2；

当判断直方图平缓阈值>φ，重复执行步骤S312至步骤S313的迭代。

步骤S314、改变df，并执行步骤S312至步骤S313的迭代，直至df小于第 二阈值，获得蛋白质系统的相对的态密度

其中，第二阈值为即本发明实施例中可取 0.0001；

且S_real(E)＝S(E)+lnk×Θ(E₀-E)，求得g(E)。

具体的改变df方式是：先连续进行N次迭代f＝f^α(0＜α＜1)，再进行1次迭 代并反复重复前述迭代方式。

从进程，

步骤S321、初始化蛋白质系统的态密度对数S(E)，直方图H(E)，修正因子 df，df|_E＝(κΘ(E₀-E)+1)lnf，

其中，Θ(E₀-E)为Heaviside分段函数，κ、E₀、f是模型的参数；本发明 实施例中κ＝5，E₀＝-2，f＝e。初始化t=0。

初始化蛋白质系统的态密度对数S(E)为S(E)＝lng(E)＝0，H(E)=0； （E_min≤E≤E_max）。

步骤S322、进行Metropolis式的随机变动，构建新构型，计算能量E_new；

并根据Metropolis准则确定新构型被接受的概率：

$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。

对t累加，即t=t+1。判定：t mod 1000=0?

步骤S323、重复步骤S322，完成进程间通信，接收步骤S313计算的全局 S(E)、H(E)，即更新原S(E)、H(E)；判断直方图平缓阈值，当满足是执行步骤 S324，否则重复执行步骤S322、步骤S323迭代；

具体为：所有从进程发送S_tmp(E)和H_tmp(E)至主进程，并接收主进程的全局 S(E)和H(E)，并更新原S(E)、H(E)，并将S_tmp(E)和H_tmp(E)初始化为0

并求得

判断直方图平缓阈值<φ，即

$\frac{\max \left(H_{\mathrm{real}}\left(E\right)\right)-\min \left(H_{\mathrm{real}}\left(E\right)\right)}{\max \left(H_{\mathrm{real}}\left(E\right)\right)+\min \left(H_{\mathrm{real}}\left(E\right)\right)}<\phi (0<\phi <1),$此时执行步骤S324，本发明实施例 中φ可取0.2；

当判断直方图平缓阈值>φ，重复执行步骤S322至步骤S323的迭代。

步骤S324、改变df，并执行步骤S322至步骤S323的迭代，直至df小于第 二阈值；

其中，第二阈值为即本发明实施例中可取 0.0001；

具体的改变df方式是：先连续进行N次迭代f＝f^α(0＜α＜1)，再进行1次迭 代并反复重复前述迭代方式。

以上为本发明一主进程以及从进程1的流程，其还包括从进程的并行进程， 进程2、。。。进程N。

本发明的基于随机优化的老年痴呆致病机理的高效分析方法，包括确定能 量模型、搜索自由能最小的蛋白质构型、模拟及计算蛋白质系统的态密度等步 骤，其与经典的WangLandau算法相比，使用具有全局更新特点的分段函数形式 的修正因子能加快搜索和模拟速度，利用退火机制的灵活的更新修正因子方式 可提高模拟精度和速度，且该方法的并行方式可大大加快搜索和模拟速度。

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