[发明专利]一种基于抽象凸下界估计的蛋白质结构预测方法

摘要

一种基于抽象凸下界估计的蛋白质结构预测方法，包括以下步骤：首先针对蛋白质高维构象空间采样难题，采用一系列变换方法将ECEPP/3力场模型转换为单位单纯性约束条件下的递增射线凸函数；基于抽象凸理论，证明并分析给出了递增射线凸函数的支撑超平面集；然后，在差分进化群体算法框架下，基于群体极小化构象次微分知识构建下界低估支撑面；进而，通过低估支撑面极值点快速枚举方法，逐步减小构象采样空间以提高采样效率;同时，利用下界低估支撑面快速廉价地估计原势能模型能量值，有效减少势能模型目标函数的评价次数；最后，甲硫氨酸—脑啡肽(TYR1-GLY2-GLY3-PHE4-MET5)构象空间优化实例验证了本发明的有效性。本发明提供一种可靠性高、复杂性较低、计算效率高的基于抽象凸下界估计的蛋白质结构预测方法。

主权项

1.一种基于抽象凸下界估计的蛋白质结构预测方法，包括以下步骤：1）选取合适的力场模型，本发明采用ECEPP/3力场模型能量函数的表示形式如下：f1(x‾1,x‾2,···,x‾N‾)=Ebond+Eangle+Etorsion+Eelectrostatic+Evdw+Ehydrogen+Eother]]>=Σb∈BONDkb2(b-b0)2+Σa∈ANGLEka2(a-a0)2+Στ∈TORΣm∈MULVm,τ[1+cos(mτ-γm,τ)]---(I)]]>+Σi,j∈ESqiqjϵ·rij+Σi,j∈VDWAijrij12-Bijrij6-Σi,j∈HBCijrij12-Dijrij10+Eother]]>式中表示肽链中原子个数，为第i个原子的坐标E_bond为键长贡献项（1-2相互作用），BOND为键长集合，b为1-2原子之间的距离，b₀为1-2原子之间平衡状态下键长，k_b为键长强度；E_angle为键角贡献项（1-3相互作用），ANGLE为键角集合，a为两键矢量之间的夹角，a₀为平衡状态下键角，k_a为键角的强度；E_torsion为正常二面角贡献项（1-4相互作用），采用余弦函数的形式描述，TOR为正常二面角集合，MUL为二面角多样性集合，τ为正常二面角，m为多样性，V_m,τ为二面角τ对应于多样性m的势垒高度，γ_m,τ为二面角τ对应于多样性m的初始相位；E_{electrostatic}为静电力（库仑力）贡献项，ES为静电作用力原子集合，q_i为原子i的部分电荷，q_j为原子j的部分电荷，ε为介电常数；E_vdw为范德华力相互作用贡献项（6-12作用力），采用Lennard-Jones势描述，VDW为范德华作用力集合，r_ij为原子i与原子j之间的欧氏距离，参数A_ij和B_ij依赖于特定原子类型和相互作用的特征；E_hydrogen为氢键相互作用贡献项（10-12作用力），HB为氢键作用力集合，C_ij和D_ij依赖于相互作用特征；E_other为其它额外的能量贡献项；2）将力场模型转换为单位单纯形约束下的目标函数，并通过局部优化的方法获得其简化的势能模型；3）参数初始化：设置群体规模popSize，变异因子F为0.5，交叉因子CR为0.1，低估概率underFactor，常数M，支撑向量规模K=N+1，支撑向量阈值K_T，初始群体P={x1,x2,...,xpopSize|xi=(x1i,x2i,...,xN+1i)∈S,i∈I},]]>其中S={x∈RN+1:xi≥0,Σi=1N+1xi=1}]]>称为单位单纯形，计算f₅(xⁱ),i∈I，其中f₅为简化势能模型对应的目标函数并设为区域I上的最小值，其中I={1,2,...,popSize}；4）对每一个目标个体xⁱ∈S(i=1,2,…,popSize)作如下处理：4.1）任意选取三个个体{x^a,x^b,x^c|a,b,c∈{1,2,...,popSize},a≠b≠c≠i}；4.2）对{x^a,x^b,x^c}执行变异操作生成变异个体4.3）对目标个体xⁱ和变异个体执行交叉操作，生成测试个体4.3.1）设置j=1;4.3.2）trial[j]=x^jiif(randb(j)≤CRorj=rnbr(i)xjiif(randb(j)>CRRorj≠rnbr(i)]]>j=1,2,...,N+1；其中randb(j)产生0到1之间的随机数；rnbr(i)产生1到N+1之间的任一个整数；4.3.3）j=j+1;4.3.4）如果j<N+1；转至步骤4.3.2）；4.4）i=i+1；如果i<popSize，转至步骤4.1）；5）对目标个体xⁱ∈S和测试个体逐个更新操作：5.1）设置i=1；5.2）查询包含的子区域计算其中为低估支撑面在子区域中唯一的最优解，为唯一对应于子区域的支撑向量矩阵对角项；5.3）如果转至5.9）；5.4）如果K<K_T，转至5.8）；5.5）如果random（0,1）<underFactor，转至5.8）；5.6）计算d(x‾i):=H~K(x‾i)=max{l1‾k1x‾1i,l‾2k2x‾2i,...,l‾N+1kN+1x‾N+1i},]]>其中为给定区域中的唯一最优解；5.7）如果转至5.9）；5.8）计算如果则置K=K+1，更新树结构T^K;5.9）i=i+1；如果i≤popSize，转至5.2）；6）置f5*=mini∈If5(xi);]]>7）判断是否满足终止条件，如不满足转至步骤3）；8）输出结果，退出。