CN103413067A - 一种基于抽象凸下界估计的蛋白质结构预测方法 - Google Patents

Abstract

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

Description

一种基于抽象凸下界估计的蛋白质结构预测方法

技术领域

本发明涉及一种生物信息学、计算机应用领域，尤其涉及的是，一种基于抽象凸下界估计的蛋白质结构预测方法。

背景技术

生物信息学是生命科学和计算机科学交叉领域的一个研究热点。生物信息学研究成果已经被广泛应用于基因发现和预测、基因数据的存储管理、数据检索与挖掘、基因表达数据分析、蛋白质结构预测、基因和蛋白质同源关系预测、序列分析与比对等。而蛋白质三维结构预测又是生物信息学领域的一个重要分支。著名的Anfinsen实验表明蛋白质的一级结构决定其三维结构，即肽链中的氨基酸残基序列决定其空间结构。从热力学的角度来讲，肽链分子天然构型对应于其物理能量的最小状态。目前，根据Anfinsen假设，直接从氨基酸序列出发，基于势能模型，采用全局优化方法，搜索分子系统的最小能量状态，从而能够高通量、廉价地预测肽链的天然构像，已经成为生物信息学最重要的研究课题之一。经过40多年的发展，对于序列相似度较高的情况(>50%)，预测精度能够达到

左右；但是序列相似度<30%的情况下并不理想；尤其是对于序列相似度低、或多肽（<10个残基的小蛋白）来说，从头预测方法则是唯一的选择。

从头预测方法是选择一种力场势能模型作为目标函数，利用全局优化算法在势能曲面上搜索全局最小能量解，在数学上属于一类极其复杂的非凸多极值优化问题。针对力场模型优化问题，1987年Li开创性地提出了一种蒙特卡洛能量最小化方法(Monte Carlo Minimization,MCM)，采用Metropolis Monte Carlo采样策略，结合局部最小化方法，优化得到了脑啡肽(包含5个残基)的最低能量结构；1997年Lee提出的构象空间退火(Conformational SpaceAnnealing,CSA)算法，基于“种子”概念生成构象群体，结合遗传算法和构象集结过程，除了得到脑啡肽的最小能量构象之外，还得到了一系列次稳态构象；1999年Wales提出的Basin-Hopping(BH)方法本质上就是结合局部优化方法的蒙特卡洛退火优化，2005年Bradley在结合蒙特卡洛优化方法开发了分子建模Rosetta程序，成功预测CASP6中T0281目标蛋白(包含70个残基)的C_α-RMSD为

对某些小目标蛋白(<85个残基)，甚至得到误差精度小于

的高分辨率三维解析结构；2006年，Zhan采用Basin-Paving(BP)方法得到了脑啡肽在ECEPP/2、ECEPP/3两种力场模型下最小能量构象；2009年Roosi提出几种BP改进策略，并给出了比较结果。另外，其它随机全局优化算法，如遗传算法(Genetic Algorithms,GA)、差分进化(Differential Evolution,DE)算法也应用于力场模型的优化，但是由于模型的复杂性，这些算法极容易陷入某个局优解；同时由于随机算法缺乏全局收敛理论依据及解的不确定性，进一步限制了它们在实际问题中的应用。

因此，现有的技术在确定蛋白质稳定构象方面存在着缺陷，需要改进。

发明内容

为了克服已有方法复杂度高，易陷入局部极值点而不易得到全局最优稳定构象的不足，本发明在降低力场模型复杂性的同时，结合抽象凸理论，提出了一种基于抽象凸下界估计的蛋白质结构预测方法。

本发明解决其技术问题所采用的技术方案是：一种基于抽象凸下界估计的蛋白质结构预测方法，包括以下步骤：

1）选取合适的力场模型。

力场模型是依赖于原子三维坐标的经验势能函数，由于其忽略了电子的相互作用，使得分子力场模型结构相对简单，计算速度较快，辅以成熟的力场参数，可达到很高的精度，能够用于生物大分子的结构预测。力场模型总势能通常被经验性的划分成若干项，本发明采用ECEPP/3力场模型能量函数的表示形式如下：

$f_1({\stackrel{\&OverBar;}{x}}^1,{\stackrel{\&OverBar;}{x}}^2,\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{x}}^{\stackrel{\&OverBar;}{N}})=E_{\mathrm{bond}}+E_{\mathrm{angle}}+E_{\mathrm{torsion}}+E_{\mathrm{electrostatic}}+E_{\mathrm{vdw}}+E_{\mathrm{hydrogen}}+E_{\mathrm{other}}$

$=\underset{b\&Element;\mathrm{BOND}}{\mathrm{\&Sigma;}}\frac{k_b}{2}{(b-b_0)}^2+\underset{a\&Element;\mathrm{ANGLE}}{\mathrm{\&Sigma;}}\frac{k_a}{2}{(a-a_0)}^2+\underset{\mathrm{\&tau;}\&Element;\mathrm{TOR}}{\mathrm{\&Sigma;}}\underset{m\&Element;\mathrm{MUL}}{\mathrm{\&Sigma;}}{V}_{m,\mathrm{\&tau;}}[1+\cos (\mathrm{m\&tau;}-{\mathrm{\&gamma;}}_{m,\mathrm{\&tau;}})]---\left(I\right)$

$+\underset{i,j\&Element;\mathrm{ES}}{\mathrm{\&Sigma;}}\frac{q_i{q}_j}{\mathrm{\&epsiv;}\&CenterDot;r_{\mathrm{ij}}}+\underset{i,j\&Element;\mathrm{VDW}}{\mathrm{\&Sigma;}}\frac{A_{\mathrm{ij}}}{{r_{\mathrm{ij}}}^{12}}-\frac{B_{\mathrm{ij}}}{r_{\mathrm{ij}}^6}-\underset{i,j\&Element;\mathrm{HB}}{\mathrm{\&Sigma;}}\frac{C_{\mathrm{ij}}}{{r_{\mathrm{ij}}}^{12}}-\frac{D_{\mathrm{ij}}}{r_{\mathrm{ij}}^{10}}+E_{\mathrm{other}}$

式中

表示肽链中原子个数，为第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）将力场模型转换为单位单纯形约束下的目标函数，并通过局部优化的方法获得其简化的势能模型。

模型（I）中b，a，τ，r_ij等变量均为

个原子坐标

的函数，经过平移及旋转变换后问题维数为

给定肽链中所有原子坐标，可计算得到能量值。通常在生物学条件下，键长偏差小于键角偏差小于2°。因此，考虑肽链分子的键长、键角均固定在平衡状态，设置E_bond=E_angle=E_other=0，可将其维数降至

其中

为肽链分子中共价键总数，

为肽链分子中独立键角总数。与此同时，该简化过程也引入高价（最高12价）的非线性等式约束条件。为了消除高价非线性等式约束条件，设 $r_{\mathrm{ij}}=\mathrm{\&zeta;}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,...,{\mathrm{\&tau;}}_N),i,j=\mathrm{1,2},...,\stackrel{\&OverBar;}{N},i\&NotEqual;j,$ 代入模型（I），可得到：

$-\mathrm{\&pi;}\&le;{\mathrm{\&psi;}}_i\&le;\mathrm{\&pi;},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{RES}}$

$-\mathrm{\&pi;}\&le;{\mathrm{\&omega;}}_i\&le;\mathrm{\&pi;},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{RES}}---\left(\mathrm{II}\right)$

$-\mathrm{\&pi;}\&le;{\mathrm{\&chi;}}_i^{\stackrel{\&OverBar;}{k}}\&le;\mathrm{\&pi;},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{RES}},\stackrel{\&OverBar;}{k}=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{K}}_i$

式中

为肽链二面角向量；N为肽链中二面角的自由度（即优化问题维数），N_RES表示肽链长度（或残基）个数，

为第i个残基侧链二面角的个数；且满足

为第i个残基主链C–N–C_α–C四个原子之间的二面角，ψ_i为第i个残基主链N–C_α–C–N四个原子之间的二面角，ω_i为第i个残基主链C_α–C–N–C_α四个原子之间的二面角，

为第i个残基侧链第

个二面角变量。

2.1）将力场模型转换为单位单纯形。

设

常数

其中τ_t(t=1,2,...,N)为模型（II）中二面角变量，N为二面角自由度（即优化问题维数）；则由模型（II）约束条件可知： ${\mathrm{\&Sigma;}}_{t=1}^N({\mathrm{\&tau;}}_i+\mathrm{\&pi;})={\mathrm{\&Sigma;}}_{t=1}^N{\mathrm{\&tau;}}_i+\mathrm{\&pi;N}\&le;\mathrm{\&Gamma;};$ ${\mathrm{\&tau;}}_t^{\&prime;}\&GreaterEqual;0,t=\mathrm{1,2},...,N.$ 采用下列投影变换：

$\left\{\begin{array}{cc}{x}_t\&equiv;\frac{{\mathrm{\&tau;}}_t^{\&prime;}}{\mathrm{\&Gamma;}}=\frac{{\mathrm{\&tau;}}_t+\mathrm{\&pi;}}{2\mathrm{\&pi;N}},& t=\mathrm{1,2},...,N\\ x_{N+1}\&equiv;1-{\mathrm{\&Sigma;}}_{t=1}^N{x}_t=\frac{\mathrm{\&pi;N}-{\mathrm{\&Sigma;}}_{t=1}^N{\mathrm{\&tau;}}_t}{2\mathrm{\&pi;N}}.& \end{array}\right.---\left(2\right)$

将模型（II）边界约束可行域一对一地映射成单位单纯形将式（2）代入模型（II）可得：

min f₃(x)=f₂(2πNx₁-π,2πNx₂-π,...,2πNx_N-π),x∈S.   （III）

近一步，应用局部优化方法，得到以下松弛模型：

$\min f_4\left(x\right)=f_3\left(\tilde{x}\right);x,\tilde{x}\&Element;S.---\left(\mathrm{IV}\right)$

其中

为以x为初始点在模型（III）势能曲面上应用局部最小化算法得到的局优解。

显然，模型（IV）为模型（III）松弛势能曲面，即min f₃(x)=min f₄(x)。图1给出了脑啡肽以ψ₁(即τ₂,酪氨酸残基中N–C_α–C–N四个原子形成的二面角)为自变量，其余23个二面角固定在平衡态时模型(III)和模型(IV)的势能曲线。

3）参数初始化:设置群体规模popSize，变异因子F为0.5，交叉因子CR为0.1，低估概率underFactor，常数M，支撑向量规模K=N+1，支撑向量阈值K_T，初始群体 $P=\{x^1,x^2,...,x^{\mathrm{popSize}}|x^i=(x_1^i,x_2^i,...,x_{N+1}^i)\&Element;S,i\&Element;I\},$ 其中 $S=\{x\&Element;R^{N+1}:x_i\&GreaterEqual;0,{\mathrm{\&Sigma;}}_{i=1}^{N+1}{x}_i=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) $\mathrm{trial}\left[j\right]=\left\{\begin{array}{cccc}{\hat{x}}_j^i& \mathrm{if}(\mathrm{randb}\left(j\right)\&le;\mathrm{CR}& \mathrm{or}& j=\mathrm{rnbr}\left(i\right)\\ x_j^i& \mathrm{if}(\mathrm{randb}\left(j\right)>\mathrm{CR}R& \mathrm{or}& j\&NotEqual;\mathrm{rnbr}\left(i\right)\end{array}\right.$ 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\left({\stackrel{\&OverBar;}{x}}^i\right):={\tilde{H}}^K\left({\stackrel{\&OverBar;}{x}}^i\right)=\max \{{\stackrel{\&OverBar;}{l_1}}^{k_1}{\stackrel{\&OverBar;}{x}}_1^i,{{\stackrel{\&OverBar;}{l}}_2}^{k_2}{\stackrel{\&OverBar;}{x}}_2^i,...,{{\stackrel{\&OverBar;}{l}}_{N+1}}^{k_{N+1}}{\stackrel{\&OverBar;}{x}}_{N+1}^i\};$

5.7）如果转至5.9）;

5.8）计算

如果

则置

K=K+1，更新树结构T^K;

5.9）i=i+1；如果i≤popSize，转至5.2）；

6)置 $f_5^{*}=\underset{i\&Element;I}{\min }{f}_5\left(x^i\right);$

7）判断是否满足终止条件，如不满足转至步骤3）；

8）输出结果，退出。

附图说明

图1、脑啡肽模型III与模型IV势能曲线比较（其余23个二面角固定在平衡状态）；

图2、脑啡肽稳定构象（能量值-11.7073Kcal/mol）；

图3、脑啡肽最优构象解对应的24个二面角的值；

具体实施方式

下面结合附图对本发明作进一步描述，参照图1～图3，以脑啡肽为实施例，一种基于抽象凸下界估计的蛋白质结构预测方法，其中包含以下步骤：

1）选取合适的力场模型：力场模型是依赖于原子三维坐标的经验势能函数，由于其忽略了电子的相互作用，使得分子力场模型结构相对简单，计算速度较快，辅以成熟的力场参数，可达到很高的精度，能够用于生物大分子的结构预测。力场模型总势能通常被经验性的划分成若干项，本发明采用ECEPP/3力场模型能量函数的表示形式如下：

$f_1({\stackrel{\&OverBar;}{x}}^1,{\stackrel{\&OverBar;}{x}}^2,\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{x}}^{\stackrel{\&OverBar;}{N}})=E_{\mathrm{bond}}+E_{\mathrm{angle}}+E_{\mathrm{torsion}}+E_{\mathrm{electrostatic}}+E_{\mathrm{vdw}}+E_{\mathrm{hydrogen}}+E_{\mathrm{other}}$

$=\underset{b\&Element;\mathrm{BOND}}{\mathrm{\&Sigma;}}\frac{k_b}{2}{(b-b_0)}^2+\underset{a\&Element;\mathrm{ANGLE}}{\mathrm{\&Sigma;}}\frac{k_a}{2}{(a-a_0)}^2+\underset{\mathrm{\&tau;}\&Element;\mathrm{TOR}}{\mathrm{\&Sigma;}}\underset{m\&Element;\mathrm{MUL}}{\mathrm{\&Sigma;}}{V}_{m,\mathrm{\&tau;}}[1+\cos (\mathrm{m\&tau;}-{\mathrm{\&gamma;}}_{m,\mathrm{\&tau;}})]---\left(I\right)$

$+\underset{i,j\&Element;\mathrm{ES}}{\mathrm{\&Sigma;}}\frac{q_i{q}_j}{\mathrm{\&epsiv;}\&CenterDot;r_{\mathrm{ij}}}+\underset{i,j\&Element;\mathrm{VDW}}{\mathrm{\&Sigma;}}\frac{A_{\mathrm{ij}}}{{r_{\mathrm{ij}}}^{12}}-\frac{B_{\mathrm{ij}}}{r_{\mathrm{ij}}^6}-\underset{i,j\&Element;\mathrm{HB}}{\mathrm{\&Sigma;}}\frac{C_{\mathrm{ij}}}{{r_{\mathrm{ij}}}^{12}}-\frac{D_{\mathrm{ij}}}{r_{\mathrm{ij}}^{10}}+E_{\mathrm{other}}$

式中

表示肽链中原子个数，

为第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）将力场模型转换为单位单纯形约束下的目标函数，并通过局部优化的方法获得其简化的势能模型。

模型（I）中b，a，τ，r_ij等变量均为

个原子坐标

的函数，经过平移及旋转变换后问题维数为

给定肽链中所有原子坐标，可计算得到能量值。通常在生物学条件下，键长偏差小于

键角偏差小于2°。因此，考虑肽链分子的键长、键角均固定在平衡状态，设置E_bond=E_angle=E_other=0，可将其维数降至

其中

为肽链分子中共价键总数，

为肽链分子中独立键角总数。与此同时，该简化过程也引入高价（最高12价）的非线性等式约束条件。为了消除高价非线性等式约束条件，设 $r_{\mathrm{ij}}=\mathrm{\&zeta;}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2,...,{\mathrm{\&tau;}}_N),i,j=\mathrm{1,2},...,\stackrel{\&OverBar;}{N},i\&NotEqual;j,$ 代入模型（I），可得到：

$-\mathrm{\&pi;}\&le;{\mathrm{\&psi;}}_i\&le;\mathrm{\&pi;},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{RES}}$

$-\mathrm{\&pi;}\&le;{\mathrm{\&omega;}}_i\&le;\mathrm{\&pi;},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{RES}}---\left(\mathrm{II}\right)$

$-\mathrm{\&pi;}\&le;{\mathrm{\&chi;}}_i^{\stackrel{\&OverBar;}{k}}\&le;\mathrm{\&pi;},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{RES}},\stackrel{\&OverBar;}{k}=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{K}}_i$

式中为肽链二面角向量；N为肽链中二面角的自由度（即优化问题维数），N_RES表示肽链长度（或残基）个数，

为第i个残基侧链二面角的个数；且满足

为第i个残基主链C–N–C_α–C四个原子之间的二面角，ψ_i为第i个残基主链N–C_α–C–N四个原子之间的二面角，ω_i为第i个残基主链C_α–C–N–C_α四个原子之间的二面角，

为第i个残基侧链第

个二面角变量。

2.1）将力场模型转换为单位单纯形。

设

其中τ_t(t=1,2,...,N)为模型（II）中二面角变量，N为二面角自由度（即优化问题维数）；则由模型（II）约束条件可知：

采用下列投影变换：

$\left\{\begin{array}{cc}{x}_t\&equiv;\frac{{\mathrm{\&tau;}}_t^{\&prime;}}{\mathrm{\&Gamma;}}=\frac{{\mathrm{\&tau;}}_t+\mathrm{\&pi;}}{2\mathrm{\&pi;N}},& t=\mathrm{1,2},...,N\\ x_{N+1}\&equiv;1-{\mathrm{\&Sigma;}}_{t=1}^N{x}_t=\frac{\mathrm{\&pi;N}-{\mathrm{\&Sigma;}}_{t=1}^N{\mathrm{\&tau;}}_t}{2\mathrm{\&pi;N}}.& \end{array}\right.---\left(2\right)$

将模型（II）边界约束可行域一对一地映射成单位单纯形将式（2）代入模型（II）可得：

min f₃(x)=f₂(2πNx₁-π,2πNx₂-π,...,2πNx_N-π),x∈S.   （III）

近一步，应用局部优化方法，得到以下松弛模型：

$\min f_4\left(x\right)=f_3\left(\tilde{x}\right);x,\tilde{x}\&Element;S.---\left(\mathrm{IV}\right)$

其中

为以x为初始点在模型（III）势能曲面上应用局部最小化算法得到的局优解。

显然，模型（IV）为模型（III）松弛势能曲面，即min f₃(x)=min f₄(x)。

为了获取势能函数低估模型，在抽象凸理论的基础上，引入以下定理。

定理1：设

使得模型（III）目标函数f₃:S→R满足：

$\stackrel{\&OverBar;}{L}=\underset{x^1\&NotEqual;x^2}{\inf }\frac{|f_3\left(x^1\right)-f_3\left(x^2\right)|}{{\left|\right|x^1-x^2\left|\right|}_1},{\&ForAll;x}^1,x^2\&Element;S\&equiv;\{x\&Element;R_{+}^{N+1}:{\mathrm{\&Sigma;}}_{i=1}^{N+1}{x}_i=1\}.---\left(3\right)$

其中 ${\left|\right|x^1-x^2\left|\right|}_1\&equiv;{\max }_{i=\mathrm{1,2},...,N+1}|x_i^1-x_i^2|;$ 如果 ${\min }_{x\&Element;S}{f}_3\left(x\right)>2\stackrel{\&OverBar;}{L},$ 则f₃:S→R₊在y∈S处的支撑函数h^y(x)为:

$h^y\left(x\right)=\underset{i\&Element;\mathrm{\&Theta;}\left(l\right)}{\min }{l}_i{x}_i,\&ForAll;x\&Element;S.---\left(4\right)$

其中l=f₃(y)/y，

索引集合Θ(l)={i:l_i>0}。

证明：首先证明f₃:S→R₊为严格递增射线凸函数（SICAR）。

考虑下列函数

$g\left(z\right)=\left\{\begin{array}{cc}{f}_3\left(\frac{z}{{\mathrm{\&Sigma;}}_{i=1}^{N+1}{z}_i}\right)\left({\mathrm{\&Sigma;}}_{i=1}^{N+1}{z}_i\right)& z\&NotEqual;0\\ 0& z=0\end{array}\right.---\left(5\right)$

（i）根据引理3， $\underset{x\&Element;S}{\min }{f}_3\left(x\right)>2\stackrel{\&OverBar;}{L}\&DoubleLeftRightArrow;2\stackrel{\&OverBar;}{L}/\underset{x\&Element;S}{\min }{f}_3\left(x\right)<1;$ 故 ${\&ForAll;x}^1,x^2\&Element;R_{+}^{N+1},$ 如果x¹>x²，则g(x¹)>g(x²)，即

为严格递增函数；设

由式（5）可知g(x)=f₃(x)，即f₃:S→R₊为严格递增函数；

（ii）定义f_3y(λ)≡{f₃(λy)|y∈S，λ∈(0，+∞)}，设

有g(y)=f₃(y)，由式（5）可得出：

$f_{3y}\left(\mathrm{\&lambda;}\right)=f_3\left(\mathrm{\&lambda;y}\right)=f_3\left(\frac{\mathrm{\&lambda;y}}{\mathrm{\&lambda;}{\mathrm{\&Sigma;}}_{t=1}^{N+1}{y}_i}\right)\mathrm{\&lambda;}\left({\mathrm{\&Sigma;}}_{t=1}^{N+1}{y}_i\right)=\mathrm{\&lambda;}{f}_3\left(y\right)---\left(6\right)$

即(f₃)_y(λ)是在y∈S点关于λ∈(0,+∞)的凸函数。

根据定义1，由（i）、（ii）可知：f₃:S→R₊为严格递增射线凸函数（SICAR）；

根据引理1，可知f₃:S→R₊在y∈S点的H-次微分为：

$h^y\left(x\right)=l\left(x\right)-l\left(y\right)+f\left(y\right)=\underset{i\&Element;\mathrm{\&Theta;}\left(l\right)}{\min }{l}_i{x}_i-\underset{i\&Element;\mathrm{\&Theta;}\left(l\right)}{\min }{l}_i{y}_i+f\left(y\right),\&ForAll;x\&Element;S.---\left(7\right)$

$l\&Element;U=\{\frac{v}{y}:v\&Element;{\&PartialD;f}_{3y}\left(1\right)\}---\left(8\right)$

基于凸分析理论可知 ${\&PartialD;f}_{3y}\left(\mathrm{\&lambda;}\right),\&ForAll;\mathrm{\&lambda;}\&Element;(0,\&infin;)$ 非空，且 ${\&PartialD;f}_{3y}\left(\mathrm{\&lambda;}\right)\&Element;[f_{3y}^{-}\left(\mathrm{\&lambda;}\right),f_{3y}^{+}\left(\mathrm{\&lambda;}\right)],$ 其中：

$f_{3y}^{+}\left(\mathrm{\&lambda;}\right)=\underset{\mathrm{\&beta;}\&DownArrow;0}{\lim }\frac{f_{3y}(\mathrm{\&lambda;}+\mathrm{\&beta;})-f_{3y}\left(\mathrm{\&lambda;}\right)}{\mathrm{\&beta;}}=\underset{\mathrm{\&beta;}\&DownArrow;0}{\lim }\frac{(\mathrm{\&lambda;}+\mathrm{\&beta;})f_3\left(y\right)-{\mathrm{\&lambda;f}}_3\left(y\right)}{\mathrm{\&beta;}}=f_3\left(y\right);---\left(9\right)$

$f_{3y}^{-}\left(\mathrm{\&lambda;}\right)=\underset{\mathrm{\&beta;}\&DownArrow;0}{\lim }\frac{f_{3y}\left(\mathrm{\&lambda;}\right)-f_{3y}(\mathrm{\&lambda;}-\mathrm{\&beta;})}{\mathrm{\&beta;}}=\underset{\mathrm{\&beta;}\&DownArrow;0}{\lim }\frac{{\mathrm{\&lambda;f}}_3\left(y\right)-(\mathrm{\&lambda;}-\mathrm{\&beta;})f_3\left(y\right)}{\mathrm{\&beta;}}=f_3\left(y\right);---\left(10\right)$

故：

${\&PartialD;f}_{3y}\left(\mathrm{\&lambda;}\right){|}_{\mathrm{\&lambda;}-1}={\&PartialD;f}_{3y}\left(1\right)=f_3\left(y\right)---\left(11\right)$

将上式代入（8），得到在y点的支撑向量：

$l=\frac{f_3\left(y\right)}{y}=(\frac{f_3\left(y\right)}{y_1},\frac{f_3\left(y\right)}{y_2},...,\frac{f_3\left(y\right)}{y_{N+1}})---\left(12\right)$

根据式（7）、（12），h^y(x)为关于y点的支撑函数为：

$h^y\left(x\right)=\underset{i\&Element;\mathrm{\&Theta;}\left(l\right)}{\min }{l}_i{x}_i,\&ForAll;x\&Element;S---\left(13\right)$

证毕。

推论1：设h^y(x)为模型（III）f₃(x)关于y点的支撑函数，则可得到：

$f_3\left(x\right)\&GreaterEqual;h^y\left(x\right),\&ForAll;x\&Element;S.---\left(14\right)$

f₃(x)=h^y(x),x=y.   (15)

证明：由式（13）可知，

$h^y\left(x\right)=\underset{i\&Element;\mathrm{\&Theta;}\left(l\right)}{\min }{l}_i{x}_i.$

不失一般性，假定

根据式（12）可得出：

$h^y\left(x\right)=f_3\left(y\right)\min \{\frac{x_1}{y_1},\frac{x_2}{y_2},...,\frac{x_{N+1}}{y_{N+1}}\}---\left(16\right)$

假设：

$\stackrel{\&OverBar;}{\mathrm{\&lambda;}}=\min \{\frac{x_1}{y_1},\frac{x_2}{y_2},...,\frac{x_{N+1}}{y_{N+1}}\}=\frac{x_{\mathrm{\&xi;}}}{y_{\mathrm{\&xi;}}},\mathrm{\&xi;}\&Element;\{\mathrm{1,2},...,N+1\}.---\left(17\right)$

则对于 $\&ForAll;t\&Element;\{\mathrm{1,2},...,N+1\},$ 满足：

$\frac{x_t}{y_t}\&GreaterEqual;\stackrel{\&OverBar;}{\mathrm{\&lambda;}},t\&Element;\{\mathrm{1,2},...,N+1\}.\&DoubleLeftRightArrow;x_t\&GreaterEqual;\stackrel{\&OverBar;}{\mathrm{\&lambda;}}{y}_t,t\&Element;\{\mathrm{1,2},...,N+1\}.\&DoubleLeftRightArrow;x\&GreaterEqual;\stackrel{\&OverBar;}{\mathrm{\&lambda;}}y.$

由于f₃:S→R₊为严格递增射线凸函数（ICAR），由函数的递增性可知：

$f_3\left(x\right)\&GreaterEqual;f_3\left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}y\right)=\stackrel{\&OverBar;}{\mathrm{\&lambda;}}{f}_3\left(y\right)=h^y\left(x\right),\&ForAll;x\&Element;S.$

特别地：当x=y时，由式（16）可知：f₃(x)=h^y(x).

证毕。

推论2：设y¹,y²,...,y^K∈S，则

$H^K\left(x\right)=\underset{k=1,...,K}{\max }{h}^k\left(x\right)=\underset{k=1,...,K}{\max }\underset{i\&Element;\mathrm{\&Theta;}\left(l\right)}{\min }{l}_i^k{x}_i---\left(18\right)$

为f₃:S→R₊的支撑函数族；则：

$H^K\left(x\right)\&le;f_3\left(x\right),\&ForAll;x\&Element;S.---\left(19\right)$

$H^K\left(x\right)=f_3\left(x\right),\&ForAll;x\&Element;\{y^1,y^2,...,y^K\}.---\left(20\right)$

证明：由推论1可知：

$f_3\left(x\right)\&GreaterEqual;h^k\left(x\right),\&ForAll;k\&Element;\{\mathrm{1,2},...,k\},\&ForAll;x\&Element;S.$

故可得到：

$f_3\left(x\right)\&GreaterEqual;\underset{k=1,...,K}{\max }{h}^k\left(x\right)=H^K\left(x\right),\&ForAll;x\&Element;S.---\left(21\right)$

设x=y^ξ,ξ∈{1,2,...,K}，由式（15）可知：

f₃(y^ξ)=h^ξ(y^ξ).   (22)

另外，由式（21）可知：

$f_3\left(y^{\mathrm{\&xi;}}\right)\&GreaterEqual;\underset{k=1,...,K}{\max }{h}^k\left(y^{\mathrm{\&xi;}}\right)=H^K\left(y^{\mathrm{\&xi;}}\right),\&ForAll;x\&Element;S.---\left(23\right)$

由式（22）、（23）可知：

f₃(y^ξ)=H^K(y^ξ),ξ∈{1,2,...,K}.   (24)

证毕。

定理2：设模型（III）目标函数f₃:S→R满足式（3）。假设

则给定y∈S，

${\stackrel{\&OverBar;}{h}}^y\left(x\right)=f_4\left(y\right)\underset{i\&Element;\mathrm{\&Theta;}\left(y\right)}{\min }\frac{x_i}{y_i},\&ForAll;x\&Element;S.---\left(25\right)$

为模型（III）目标函数f₃:S→R在

处的支撑弱函数。其中f₄:S→R为模型（IV）目标函数，Θ(y)={i:y_i>0}。

证明：由于f₄(x)≤f₃(x),

故：

${\stackrel{\&OverBar;}{h}}^y\left(x\right)=f_4\left(y\right)\underset{i\&Element;\mathrm{\&Theta;}\left(y\right)}{\min }\frac{x_i}{y_i}\&le;f_3\left(y\right)\underset{i\&Element;\mathrm{\&Theta;}\left(y\right)}{\min }\frac{x_i}{y_i}\&le;f_3\left(x\right),\&ForAll;x\&Element;S.---\left(26\right)$

证毕。

注：定理2表明虽然不是松弛曲面f₄:S→R的支撑弱函数，但仍然为模型（III）目标函数f₃:S→R的下界估计曲面。

定理3：给定下列势能模型：

$\underset{x\&Element;S}{\min }{f}_5\left(x\right)=f_4\left(x\right)+M.---\left(V\right)$

其中f₄:S→R为模型（IV）目标函数，M>0为常数。如果选取

式中

满足式（3），f₃:S→R为模型（III）目标函数，则给定y∈S：

${\tilde{h}}^y\left(x\right)=(f_4\left(y\right)+M)\underset{i\&Element;\mathrm{\&Theta;}\left(y\right)}{\min }\frac{x_i}{y_i},\&ForAll;x\&Element;S.---\left(27\right)$

为模型（V）的支撑弱函数，式中Θ(y)={i:y_i>0}。

注：由于f₅:S→R₊仅仅是在模型（IV）f₄:S→R上增加了常数项，故二者的最优解是等价的，又因为

故模型（III）、（IV）、（V）的最优解是等价的。

推论3：设x¹,x²,...,x^K∈S，则

${\tilde{H}}^K\left(x\right)=\underset{k=1,...,K}{\max }{\tilde{h}}^k\left(x\right)=\underset{k=1,...,K}{\max }\underset{i\&Element;\mathrm{\&Theta;}\left(l\right)}{\min }{l}_i^k{x}_i---\left(28\right)$

为模型（V）f₅:S→R₊的支撑函数族。式中

Θ(l)={i:l_i>0}。

3）参数初始化:设置群体规模popSize，变异因子F为0.5，交叉概率CR为0.1，低估概率underFactor，常数M，支撑向量规模K=N+1，支撑向量阈值K_T，初始群体 $P=\{x^1,x^2,...,x^{\mathrm{popSize}}|x^i=(x_1^i,x_2^i,...,x_{N+1}^i)\&Element;S,i\&Element;I\},$ 计算f₅(xⁱ),i∈I。并设 $f_5^{*}=\underset{i\&Element;I}{\min }{f}_5\left(x^i\right),$ 其中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） $\mathrm{trial}\left[j\right]=\left\{\begin{array}{cccc}{\hat{x}}_j^i& \mathrm{if}(\mathrm{randb}\left(j\right)\&le;\mathrm{CR}& \mathrm{or}& j=\mathrm{rnbr}\left(i\right)\\ x_j^i& \mathrm{if}(\mathrm{randb}\left(j\right)>\mathrm{CR}R& \mathrm{or}& j\&NotEqual;\mathrm{rnbr}\left(i\right)\end{array}\right.$ 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，转至步骤2.2）；

5）对目标个体xⁱ∈S和测试个体

逐个更新操作：

5.1）设置i=1；

5.2）查询包含

的子区域

计算

其中

为低估支撑面

在

子区域中唯一的最优解，

为唯一对应于子区域的支撑向量矩阵对角项。

为了查询给定区域的最优解，引入以下定理。

考虑N+1维支持向量矩阵L：

假定L满足以下两个条件：

(I) $\&ForAll;i,j\&Element;\{k_1,k_2,\&CenterDot;\&CenterDot;\&CenterDot;k_{N+1}\},i\&NotEqual;j:l_i^{k_i}<l_i^{k_j}$

(II) $\&ForAll;v\&Element;\{l^1,l^2,...,l^K\}\backslash \{l^{k_1},l^{k_2},...,l^{k_{N+1}}\},$ 则 $\&Exists;i\&Element;\{k_1,k_2,\&CenterDot;\&CenterDot;\&CenterDot;k_{N+1}\}$ 使得 $l_i^{k_i}\&GreaterEqual;{\mathrm{\&nu;}}_i$

所有的局优解对应于满足上述两个条件的所有矩阵的集合，对于上述L的局优值d(L)和局优解x_min(L)可分别表示为：

$d\left(L\right)=1/{\mathrm{\&Sigma;}}_{i=1}^{N+1}(1/l_i^{k_i})---\left(30\right)$

x_min(L)=d(L)Trace(L)   (31)

枚举所有满足条件(I)、(II)的矩阵L，可以得到

所有的局优解，其计算复杂度为

而且满足上述两个条件的矩阵L位于一个定向的无环图上，通过检查图的叶节点V^K-1（对应于的所有局优解）可发现满足条件(I)、(II)的所有组合，这样计算复杂度简化至O(N|V^K-1|))，|V^K-1|表示

局优解的数量；同时新的支撑向量l^K且满足(I)、(II)的矩阵L可从V^K-1中刚删掉的叶节点得到（v=l^K时可能导致V^K-1某些叶节点不满足(II)条件）。通过从旧节点中生成新的组合可创建一个树图T^K，假设T^K为平衡树，在算法中保存整棵树结构可使计算复杂度降至log_N+1(|V^K|)。而且，在可行域S上每一个局优解

唯一确定一个子区域S^u(u=1,2,...,U)，且子区域之间互不交叉，即：

$S=\underset{u=1}{\overset{U}{\&cup;}}{S}^u,S^{u_i}\&cap;S^{u_j},u_i\&NotEqual;u_j,u_i,u_j\&Element;\{\mathrm{1,2},...,U\}.---\left(32\right)$

对每一个子区域S^u(u=1,2,...,U)，有

${\tilde{H}}^K\left(x\right)=\max \{l_1^{k_1}{x}_1,l_2^{k_2}{x}_2,...,l_{N+1}^{k_{N+1}}{x}_{N+1}\},\&ForAll;x\&Element;S^u.---\left(33\right)$

${\tilde{H}}^K\left(x^{u^{*}}\right)=1/{\mathrm{\&Sigma;}}_{i=1}^{N+1}(1/l_i^{k_i}),\&ForAll;x\&Element;S^u.---\left(34\right)$

其中

为满足条件(I)、(II)的支撑向量矩阵。

显然， ${\tilde{H}}^K\left(x^{u^{*}}\right)\&le;{\tilde{H}}^K\left(x\right),$ $\&ForAll;x\&Element;S^u.$ 设 $\&ForAll;y\&Element;S\backslash S^u,$ 如果 $f_5\left(y\right)<{\tilde{H}}^K\left(x^{u^{*}}\right),$ 则子区域S^u肯定不会包含全局最优解，从而不必计算f₅(x),x∈S^u而将该区域可靠排除；另一方面式（32）也给出了f₅(x),x∈S^u相对廉价的下界估计值，当K→∞时，

5.3）如果

转至5.9）；

5.4）如果K<K_T，转至5.8）；

5.5）如果random（0,1）<underFactor，转至5.8）；

5.6）计算 $d\left({\stackrel{\&OverBar;}{x}}^i\right):={\tilde{H}}^K\left({\stackrel{\&OverBar;}{x}}^i\right)=\max \{{{\stackrel{\&OverBar;}{l}}_1}^{k_1}{\stackrel{\&OverBar;}{x}}_1^i,{{\stackrel{\&OverBar;}{l}}_2}^{k_2}{\stackrel{\&OverBar;}{x}}_2^i,...,{{\stackrel{\&OverBar;}{l}}_{N+1}}^{k_{N+1}}{\stackrel{\&OverBar;}{x}}_{N+1}^i\},$ 其中

为给定区域的有效下界估计值；

5.7）如果

转至5.9）；

5.8）计算

如果

则置

K=K+1，更新树结构T^K；

5.9）i=i+1；如果i≤popSize，转至5.2）；

6）置 $f_2^{*}=\underset{i\&Element;I}{\min }{f}_5\left(x^i\right);$

7）判断是否满足终止条件，如不满足转至步骤3）；

8）输出结果，退出。

以脑啡肽为实施例，运用以上方法得到了其最优构象值为-11.7073Kcal/mol。脑啡肽稳态结构如图2所示，稳态结构下对应的24个二面角如图3所示。

以上阐述的是本发明给出的一个实施例表现出来的优良优化效果，显然本发明不仅适合上述实施例，在不偏离本发明基本精神及不超出本发明实质内容所涉及内容的前提下可对其做种种变化加以实施。

Claims (3)

1.一种基于抽象凸下界估计的蛋白质结构预测方法，包括以下步骤：

1）选取合适的力场模型，本发明采用ECEPP/3力场模型能量函数的表示形式如下：

$f_1({\stackrel{\&OverBar;}{x}}^1,{\stackrel{\&OverBar;}{x}}^2,\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{x}}^{\stackrel{\&OverBar;}{N}})=E_{\mathrm{bond}}+E_{\mathrm{angle}}+E_{\mathrm{torsion}}+E_{\mathrm{electrostatic}}+E_{\mathrm{vdw}}+E_{\mathrm{hydrogen}}+E_{\mathrm{other}}$

$=\underset{b\&Element;\mathrm{BOND}}{\mathrm{\&Sigma;}}\frac{k_b}{2}{(b-b_0)}^2+\underset{a\&Element;\mathrm{ANGLE}}{\mathrm{\&Sigma;}}\frac{k_a}{2}{(a-a_0)}^2+\underset{\mathrm{\&tau;}\&Element;\mathrm{TOR}}{\mathrm{\&Sigma;}}\underset{m\&Element;\mathrm{MUL}}{\mathrm{\&Sigma;}}{V}_{m,\mathrm{\&tau;}}[1+\cos (\mathrm{m\&tau;}-{\mathrm{\&gamma;}}_{m,\mathrm{\&tau;}})]---\left(I\right)$

$+\underset{i,j\&Element;\mathrm{ES}}{\mathrm{\&Sigma;}}\frac{q_i{q}_j}{\mathrm{\&epsiv;}\&CenterDot;r_{\mathrm{ij}}}+\underset{i,j\&Element;\mathrm{VDW}}{\mathrm{\&Sigma;}}\frac{A_{\mathrm{ij}}}{{r_{\mathrm{ij}}}^{12}}-\frac{B_{\mathrm{ij}}}{r_{\mathrm{ij}}^6}-\underset{i,j\&Element;\mathrm{HB}}{\mathrm{\&Sigma;}}\frac{C_{\mathrm{ij}}}{{r_{\mathrm{ij}}}^{12}}-\frac{D_{\mathrm{ij}}}{r_{\mathrm{ij}}^{10}}+E_{\mathrm{other}}$

式中

表示肽链中原子个数，为第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=\{x^1,x^2,...,x^{\mathrm{popSize}}|x^i=(x_1^i,x_2^i,...,x_{N+1}^i)\&Element;S,i\&Element;I\},$ 其中 $S=\{x\&Element;R^{N+1}:x_i\&GreaterEqual;0,{\mathrm{\&Sigma;}}_{i=1}^{N+1}{x}_i=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） $\mathrm{trial}\left[j\right]=\left\{\begin{array}{cccc}{\hat{x}}_j^i& \mathrm{if}(\mathrm{randb}\left(j\right)\&le;\mathrm{CR}& \mathrm{or}& j=\mathrm{rnbr}\left(i\right)\\ x_j^i& \mathrm{if}(\mathrm{randb}\left(j\right)>\mathrm{CR}R& \mathrm{or}& j\&NotEqual;\mathrm{rnbr}\left(i\right)\end{array}\right.$ 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\left({\stackrel{\&OverBar;}{x}}^i\right):={\tilde{H}}^K\left({\stackrel{\&OverBar;}{x}}^i\right)=\max \{{\stackrel{\&OverBar;}{l_1}}^{k_1}{\stackrel{\&OverBar;}{x}}_1^i,{{\stackrel{\&OverBar;}{l}}_2}^{k_2}{\stackrel{\&OverBar;}{x}}_2^i,...,{{\stackrel{\&OverBar;}{l}}_{N+1}}^{k_{N+1}}{\stackrel{\&OverBar;}{x}}_{N+1}^i\},$ 其中

为给定区域中的唯一最优解；

5.7）如果

转至5.9）；

5.8）计算

如果则置K=K+1，更新树结构T^K;

5.9）i=i+1；如果i≤popSize，转至5.2）；

6）置 $f_5^{*}=\underset{i\&Element;I}{\min }{f}_5\left(x^i\right);$

7）判断是否满足终止条件，如不满足转至步骤3）；

8）输出结果，退出。

2.如权利要求1所述，步骤2）中模型变换方法为，模型（I）中b，a，τ，r_ij等变量均为

个原子坐标

的函数，经过平移及旋转变换后问题维数为

给定肽链中所有原子坐标，可计算得到能量值；通常在生物学条件下，键长偏差小于

键角偏差小于2°；因此，考虑肽链分子的键长、键角均固定在平衡状态，设置E_bond=E_angle=E_other=0，可将其维数降至

其中

为肽链分子中共价键总数，

为肽链分子中独立键角总数；与此同时，该简化过程也引入高价（最高12价）的非线性等式约束条件；为了消除高价非线性等式约束条件，设

代入模型（I），可得到：

$-\mathrm{\&pi;}\&le;{\mathrm{\&psi;}}_i\&le;\mathrm{\&pi;},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{RES}}$

$-\mathrm{\&pi;}\&le;{\mathrm{\&omega;}}_i\&le;\mathrm{\&pi;},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{RES}}---\left(\mathrm{II}\right)$

$-\mathrm{\&pi;}\&le;{\mathrm{\&chi;}}_i^{\stackrel{\&OverBar;}{k}}\&le;\mathrm{\&pi;},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{RES}},\stackrel{\&OverBar;}{k}=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,{\stackrel{\&OverBar;}{K}}_i$

式中

为肽链二面角向量；N为肽链中二面角的自由度（即优化问题维数），N_RES表示肽链长度（或残基）个数，

为第i个残基侧链二面角的个数；且满足

为第i个残基主链C–N–C_α–C四个原子之间的二面角，ψ_i为第i个残基主链N–C_α–C–N四个原子之间的二面角，ω_i为第i个残基主链C_α–C–N–C_α四个原子之间的二面角，

为第i个残基侧链第个二面角变量。

3.如权利要求1和权利要求2所述，将给定的力场模型转换为单位单纯形，并获得简化的力场模型：

设常数

其中τ_t(t=1,2,...,N)为模型（II）中二面角变量，N为二面角自由度（即优化问题维数）；则由模型（II）约束条件可知： ${\mathrm{\&Sigma;}}_{t=1}^N({\mathrm{\&tau;}}_i+\mathrm{\&pi;})={\mathrm{\&Sigma;}}_{t=1}^N{\mathrm{\&tau;}}_i+\mathrm{\&pi;N}\&le;\mathrm{\&Gamma;};{\mathrm{\&tau;}}_t^{\&prime;}\&GreaterEqual;0,t=\mathrm{1,2},...,N;$ 采用下列投影变换：

$\left\{\begin{array}{cc}{x}_t\&equiv;\frac{{\mathrm{\&tau;}}_t^{\&prime;}}{\mathrm{\&Gamma;}}=\frac{{\mathrm{\&tau;}}_t+\mathrm{\&pi;}}{2\mathrm{\&pi;N}},& t=\mathrm{1,2},...,N\\ x_{N+1}\&equiv;1-{\mathrm{\&Sigma;}}_{t=1}^N{x}_t=\frac{\mathrm{\&pi;N}-{\mathrm{\&Sigma;}}_{t=1}^N{\mathrm{\&tau;}}_t}{2\mathrm{\&pi;N}}.& \end{array}\right.---\left(2\right)$

将模型（II）边界约束可行域一对一地映射成单位单纯形

将式（2）代入模型（II）可得：

min f₃(x)=f₂(2πNx₁-π,2πNx₂-π,...,2πNx_N-π),x∈S.   （III）

近一步，应用局部优化方法，得到以下松弛模型：

$\min f_4\left(x\right)=f_3\left(\tilde{x}\right);x,\tilde{x}\&Element;S.---\left(\mathrm{IV}\right)$

其中

为以x为初始点在模型（III）势能曲面上应用局部最小化算法得到的局优解；

显然，模型（IV）为模型（III）松弛势能曲面，即minf₃(x)=minf₄(x)。

Images