改进的模拟退火算法在蛋白质结构预测中的应用

2013-12-24朱均燕温永仙

福建农林大学学报（自然科学版） 2013年2期

朱均燕，温永仙

(福建农林大学计算机与信息学院，福建福州350002)

蛋白质是一种生物大分子，是由20种氨基酸以肽键连接成的肽链［1］.蛋白质结构预测对于探索蛋白质空间结构与功能的关系，以及进行蛋白质结构设计、突变体设计等具有重要意义.蛋白质结构预测是指直接从氨基酸序列推断某一蛋白质的功能位点或对其三维结构进行预测(包括二级和三级结构)，是目前分子生物学研究中迫切需要解决的问题［2］.蛋白质二级结构预测不仅是联系蛋白质一级结构和三级结构的纽带，还是从一级结构预测其三级结构的关键步骤［3］.

迄今为止，对蛋白质结构预测问题已提出了一些简化模型.目前有2种典型的蛋白质折叠的简化模型，一种是由Dill et al提出的HP格点模型(HP lattice model)［4，5］，但该模型仅考虑了蛋白质残基之间的疏水性，忽略了残基之间的亲水性，且相邻残基之间的夹角只能是直角或平角;另一种是由Stillinger et al提出的Toy模型(AB off-lattice model)［6］，该模型同时考虑了蛋白质残基之间的疏水性和亲水性，而且相邻残基之间的夹角可是任意的.研究［7］表明，Toy模型比HP模型更接近真实蛋白质，由此，本文选用Toy模型作为蛋白质结构预测的研究对象.

1995年Stillinger对Toy模型所体现出来的蛋白质折叠的共有特点进行了研究，并提出了用于Toy模型的2种氨基酸序列，即Center-doped序列和Fibonacci序列［8］.由于Center-doped序列考虑的情况比较特殊，主链上的氨基酸序列必须是对称的，这与多数天然蛋白质序列结构不符，所以本文选用Fibonacci序列.

1 改进的模拟退火算法

模拟退火算法最初的思想是Metropolis在1953年提出的，Kirkpatrick于1983年将其成功地应用在组合最优化问题中［9］.其主要特点是在搜索的过程中，不仅接受优化解，也以一定的概率接受恶化解，这样就容易陷入局部最优解.文献［10］针对蛋白质结构预测模型中目标函数的多极值多变量的特点，给出了一种新解的产生方法，本文对这一算法进行了改进.

1.1 新解边界值处理

在蛋白质Toy二维模型中，通过找到最佳的折叠角度获得最小能量值，针对这一问题，本文在模拟退火算法中对于产生新解边界值的处理给出一种新方法.对于产生的一个新解X中的任意一个分量xr，按如下规则进行边界处理:

经式(1)处理后，新解的任意一个分量的值都在规定的范围内，且分量值符号保持不变.

1.2 初始温度的选择

在模拟退火算法中，温度是一个比较关键的参数.如果初始温度不够高或者退火时间不够长会使得搜索过快，从而导致算法陷入局部最优解;但是如果初始温度过高或者退火时间过长，则会造成因算法的运行时间过长.温度参数的设定将影响到模拟退火算法是否能收敛到全局最优解.

虽然模拟退火算法的求解不依赖于初始值，但初始温度的选择也不是随机的.初始温度的选择与临界温度有关［11］.临界温度是指液体变成固体时的温度，从模拟退火算法来说是指目标函数值开始变化时的温度［12］.采用文献［12］中区间分半搜寻临界温度的方法来确定初始温度.

1.3 算法的具体步骤

(1)用区间分半法确定初始温度tmax，设置同一温度下的迭代次数k、邻域规模因子λ、温度下降因子dt、折叠角度的最大值和最小值.

(2)利用xr=rand*(max－min)+min产生初始解X0(x01，x02，…，x0n)，计算目标函数值f0=f(X0).

(3)判断是否满足程序终止条件，如果满足就结束程序;否则，令t=t*dt，k=1，转下一步.

(4)从{1，2，…，n}中随机选一数 r，xr=x0r+λ*rand*(max－min)，rand产生－1到1之间的随机数.如果xr超出上下边界最大和最小值，则按式(1)进行边界处理后，计算产生的新解的目标函数f=f(X).

(5)计算△f=f(X)－f(X0)，根据Metropolis法则来判断是否接受新解，若接受，则用X代替X0，用f(X)代替f(X0).

(6)判断是否达到最大迭代值k，若达到则转步骤(3);否则迭代次数加1，转步骤(4).

在算法中有以下几个可调参数.

(1)初始温度Tmax:主要根据目标函数取值范围来确定，取值应保证初始接受率足够高.

(2)终止条件:在最优值未知的情况下，终止条件比较难确定.一般采取2种准则:一是给定终止温度;二是连续多次降温，直到能量函数的值不再下降为止.

(3)温度下降因子 dt:一般取0.95－0.98.

(4)同一温度迭代次数(Markov链长度):Lmax的选取与问题规模和解空间大小有关.

(5)邻域规模因子scale:邻域规模因子的取值和解空间有直接关系，取0.2－0.5.

(6)最大、最小值与角度的取值坐标有关，选逆时针方向为正方向，最大值取π，最小值取－π.

2 结果与分析

本文数据是选用Fibonacci序列形式.首先对于N＜8序列进行了测试，结果表明本文的算法能迅速获得最小能量值，并与Stillinger所获得的最小能量值［6，7］完全一致.随后又对N值分别为13、21、34、55的序列进行了测试，并与利用其它算法所得的最低能量值进行比较，结果如表1所示.

表1 序列及最小能量值1)Table 1 Sequences and the minimum values of energy

由表1可知:采用本文的模拟退火算法获得的最小能量值都优于Stillinger获得的最小能量值，也优于文献［10］中改进的模拟退火算法获得的最小能量值.对于N值分别为13、21、55的序列，采用本文算法获得的最小能量值优于用PERM方法得到的最小能量值;对于N值为34的序列，采用本文算法获得的最小能量值与采用PERM方法得到的最小能量值有差异，但两者的结果还比较接近.

图1是采用本文改进的模拟退火算法获得的表1中各序列的最小能量构象.

图1 各序列最小能量构象Fig.1 The minimum values of energy conformation form

3 小结与讨论

本文对模拟退火算法中产生的新解边界值的处理给出一种新方法.在初始温度选择方面，采用李丽等提出的区间分半搜寻临界温度的方法来确定算法的初始温度［12］.将改进后的算法应用到蛋白质二维的Toy模型中，通过对序列长度分别是13、21、34、55的Fibonacci序列测试结果的分析，证明该算法可行有效.但是由于天然蛋白质结构远比Toy模型复杂得多，所以模拟退火算法在蛋白质结构预测中的应用有待更深入研究.

［1］陈清西，李松刚.KCIO3诱导龙眼成花及其叶片碳水化合物与蛋白质的变化［J］.福建农林大学学报:自然科学版，2004，33(2):182－185.

［2］马栋苹，阮晓钢.基于改进BP神经网络预测蛋白质二级结构［J］.北京联合大学学报:自然科学版，2005，19(2):70－73.

［3］WANG Z X.The current situation and prospect of protein structure prediction［J］.Chemistry of Life，1998，18(6):19－22.

［4］SHORTLE D，CHAN H S，DILL K A.Modeling the effects of mutations on the denatured states of proteins［J］.Protein Science，1992(1):201－205.

［5］DILL K A，BROMBERG S，YUE K，et al.Principles of protein folding-A perspective from simple exact models［J］.Protein Science，1995(4):561－602.

［6］FRANK H S，TERESA H G，CATHERINE L H.Toy model for protein folding［J］.Physical Review E，1993，48(2):1469－1477.

［7］张晓龙，李婷婷，芦进.基于Toy模型蛋白质折叠预测的多种群微粒群优化算法研究［J］.计算机科学，2008，35(10):230－235.

［8］FRANK H S.Collective aspects of protein folding illustrated by a toy model［J］.Physical Review E，1995，52(3):2872－2877.

［9］王翼飞，史定华.生物信息学——智能化算法及其应用［M］.北京:化学工业出版社，2006:195－197.

［10］张红娟.基于非格点模型的蛋白质结构预测研究［C］.大连:大连理工大学，2005:30－31.

［11］姚姚.地球物理非线性反演模拟退火算法的改进［J］.地球物理学报，1995，38(5):643－650.

［12］李丽，朱国同，陈秀娟，等.模拟退火算法的改进及在静校正中的应用［J］.大庆石油地质与开发，2008，10，27(5):120－123.