计算机科学 ›› 2013, Vol. 40 ›› Issue (9): 212-215.
程正华,张贵军,邓勇跃,金媚媚
CHENG Zheng-hua,ZHANG Gui-jun,DENG Yong-yue and JIN Mei-mei
摘要： 针对现阶段药物设计中对于蛋白质结构多模态的需求,提出了一种基于排挤差分进化策略的多模态优化算法。为了降低蛋白质构象空间求解的复杂度,算法采用能量极小化过程,有效缩小了可行域的搜索空间；同时,为了有效地平衡多模态优化问题的局部收敛性和模态多样性,在排挤差分进化算法的框架下,在保证算法收敛速度的前提下,算法采用空间局部性原理,同时随机选取不同交叉策略的集结思想又有效改善了种群的多样性。以脑啡肽为例,算法不仅得到了其全局最稳定结构,还获得了一系列局部最优结构。
| [3] 许忠能．生物信息学[M].北京:清华大学出版社,2008:361-363 [4] Unger R,Moult J．Genetic algorithms for protein folding simula-tions[J]．Journal of Molecular Biology,1993,1(1):75-81 [5] Kalegari D H,Lopes H S．A differential evolution approach for protein structure optimization using a 2D off-lattice model[J]．International Journal of Bio-Inspired Computation,2010,2(3/4):242-250 [6] 。 3)执行交叉操作前,随机从M个种子中选取一个做交叉操作中的种子个体。 4)在交叉操作过程中,若选取的是小组或集合组,将种子中和选取的集合组或小组相对应的局部片段直接复制给测试个体。 4 案例研究 脑啡肽(Try1-Gly2-Gly3-Phe4-Met5)是由5个氨基酸组成的,它作为一种理想的评估新的蛋白质预测算法效率的基准,现在被广泛应用在各种蛋白质结构预测算法中。脑啡肽分子由75个原子组成,可用24个独立的主-侧链二面角描述,公认的脑啡肽稳定结构能量值为－11.7073kcal/mol,其势能面上至少包含1011个以上的局部最优解,优化难度极大,为了说明该优化问题的复杂性,我们使用多起点局部搜索方法对模型(2)进行优化:首先,随机产生100000个初始点；然后以这些起点进行能量极小化处理得到各模态的局部最优解；最后将得到的100000个局部最优解与公认的标准解进行结构对比,并计算Ca-RMSD。从图1中可以清楚地看到,100000个随机结构的局部优化都没有收敛到公认的脑啡肽稳定结构解。 图1 100000局部稳定结构与全局稳定结构对比结果 算法以脑啡肽为例,将脑啡肽的二面角搜索范围定在－180°到180°之间,并将其对应的24个二面角分为8个小组,如图2所示,小组中的φ、ψ、ω代表脑啡肽主链中的二面角,χi代表脑啡肽侧链中的二面角。在算法中,我们进一步将8个小组分为7个集合组,如表1所列。小组和集合组中的成员分别对应24个二面角中的某些片段,这些小组和集合组类似于原始集结过程中两种不同的片段。 图2 将24个二面角分为1－8个小组表1 将8个小组划分为7个集合组 小组集合组 11,2,3 21,2,4 33,4,5 44,5,6 55,6,7 65,6,8 76,7,8 设计一种蛋白质结构多模态预测的优化算法必须要解决的难点是:由于蛋白质势能曲面上拥有大量的局部极小值,虽然算法得到大量的蛋白质结构,但是其中大部分是蛋白质的重叠结构,或者不是蛋白质稳定的模态结构。针对实验中所得的重叠结构,本文设定两个标准(式(4)和式(5))来区分不同的模态:1)两个蛋白质结构解中二面角的dij不大于24°(dij中参数参考文献 [7] )；2)两个蛋白质结构解中任意二面角θ相差不大于5°。如果满足1)、2)中任意一个,则认为这两个蛋白质结构为同一个模态。针对蛋白质的伪模态结构,算法同时采用能量极小化对这些解进行处理,使其尽可能地逼近离它们最近的模态。 dij=∑nk=1min[mod{(θik－θjk),sym(k)},{sym(k)－mod{(θik－θjk),sym(k)}}]≤24°(4) ｜θ1i－θ2j｜≤5°(i≠j=1,2,…,24)(5) 5 算法测试及其结果分析 为了同BECDE-SL算法进行对比,文中给出了其它3种同类的算法。算法1为DE算法+能量极小化,算法2为DE算法+集结过程+能量极小化,算法3为CrowdingDE算法+集结过程+能量极小化,算法4为BECDE-SL算法。算法1－4各自针对模型(2)进行优化运算,并对实验结果进行对比。为了保证实验结果对比的公平性,算法1和2、3、4的实验参数保持一致,设定如下:popSize=100,mutatorFactor=0.9,crossFactor=0.1,迭代次数为400次,每种算法都独立运行50次,实验中选取的种子个数M=10,随机概率m=0.5,n=0.2,p=0.3。针对ECEPP/3势能模型,当脑啡肽的能量值收敛至－11.7070kcal/mol时,我们认定其达到最稳态结构。同时在对数据进行处理时,当迭代至400次仍然没有达到稳态值时,我们认定算法在400次迭代时达到稳态值。实验以文献[13]中所得到的脑啡肽全局稳定结构解为公认稳定结构,实验结果以此解为标准进行对比。 表2 4种算法各方面性能的对比算法可靠性(%)达到稳态的平均迭代次数每次所得模态的平均 个数50次所得模态总数 (＜－10kcal/mol)最低稳定 结构(kcal/mol) 172265.6648.1630－11.7073 274215.451.125－11.7073 388181.0695.2867－11.7073 492188.6896.6484－11.7073 从表2运行结果中可以看出,算法1－4虽然都能搜索到脑啡肽能量最小值－11.7073kcal/mol对应的结构,但是加入了集结过程和排挤因子的算法3、4运行结果明显较好。和算法3比较,空间局部性原理的加入虽然增加了算法4的复杂度,但是算法4在单个模态上的收敛速度(达到稳态的平均迭代次数/50次所得到模态总数)比算法3要好(即表2的数据188.68/84<181.06/67)；同时在保证收敛速度的前提下,算法4在50次运行中有46次能够搜索到脑啡肽的全局最稳态结构,可靠性相对最好。如图3所示的是4种算法50次运行的种群平均能量分布图,算法4其种群的平均能量曲线整体较为平缓,并且误差波动小,表明算法4拥有更好的性能,得到了更多的脑啡肽模态结构。 图3 种群运行50次的平均能量图 由于蛋白质结构的多重对称性和实验的误差,多模态优化算法得到了大量的蛋白质结构,但其中有蛋白质的重叠结构和不稳定的伪模态结构。图4所示的是算法1、2、3、4在50次运行后得到的所有能量值小于－10kcal/mol的种群个体的聚类曲线,从图中可以直观地看出,算法1是DE+能量极小化,导致种群收敛特性不是很好,400次迭代后,种群中也出现了一些脑啡肽的结构。从图4中看出,尽管算法3和算法4都是多模态优化算法,但是算法4种群个体的聚类曲线明显更加曲折,表明其有更好的种群多样性。算法4加入了空间局部原理,虽然增加小部分收敛速度,但是其平均每次运行就能得到96.64个不同的蛋白质结构,50次运行总共得到了679个能量值小于－10kcal/mol的脑啡肽结构,经过dij等限制条件的筛选,得到了84个独立的高质量脑啡肽稳定结构,明显比算法3所得到的运算结果要好。 图4 各算法50次运行所得的模态 图5 算法4所得到的6个高质量稳定结构的三维图 (下转第229页)(上接第215页) 表3是算法4所得到的全局最稳态结构的二面角解。表4给出的是算法4同CSA等算法所得的全局最稳态结构的精度对比。尽管这几种算法都能得到脑啡肽的全局最稳态结构,但是算法4所求的最优解同公认稳定结构之间的标准差为0.606°,最接近公认的脑啡肽稳定结构。图5所示的是算法4所得到的脑啡肽全局稳定结构和其它5个高质量的局部稳定结构的三维图。从图中可以直观地看出,尽管能量值相差很小,但6个三维结构图有很大的差别。图5中还给出了这6种高质量脑啡肽稳定结构与公认稳定结构相比较的Ca-RMSD。如图5中所示,能量值为－10.8748kcal/mol所对应的脑啡肽结构尽管能量值相比较大,但是其结构和公认稳定结构的误差为0.471。由此表明,优化算法所得到的脑啡肽结构虽然其能量值较小,但是其三维结构和脑啡肽的全局最稳定结构可能有较大的误差。 表3 算法4运行所得最稳定结构解 氨基酸符号二面角(°)氨基酸符号二面角(°) 酪氨酸 Try1 甘氨酸 Gly2 甘氨酸 Gly3 1 φ1 ω1χ1χ2χ3 2 φ2 ω2 3 φ3 ω3 －83.499 155.793 －177.127 －173.177 79.384 －166.329 －154.261 85.817 168.503 82.918 －75.029 －169.97 苯丙 氨酸 Phe4 甲硫 氨酸 Met5 4－136.847 φ419.097 ω4－174.09 χ158.857 χ2－85.476 5－163.446 φ5160.899 ω5－179.79 χ152.863 χ2175.293 χ3－179.868 χ4－58.589表4 各算法所得最优解的精度对比算法CSA1234 方差()0.6720.6250.6180.6280.606 标准差()0.8120.7910.7860.7930.779结束语 本文针对现在药物设计对蛋白质结构多模态的需求,提出了一种新的多模态混合算法BECDE-SL,新算法可较好地解决多模态问题,能优化存在的无法有效地平衡收敛性和多样性的问题,特别是针对蛋白质结构预测这种高维非凸函数的优化。算法在CrowdingDE-SL算法的框架下,对种群个体进行能量极小化处理,同时集结过程和空间局部原理大大加快了算法的收敛速度,不同的交叉策略又增加了种群的多样性。仿真结果表明,算法BECDE-SL同其它算法相比,不仅具有较好的收敛特性,而且保持了种群的多样性,能够得到众多高质量的脑啡肽稳定结构。 Bradley P,Misura K M S,Baker D．Toward high-resolution de novo structure prediction for small proteins[J]．Science,2005,9(5742):1868-1871[3]许忠能．生物信息学[M].北京:清华大学出版社,2008:361-363[4]Unger R,Moult J．Genetic algorithms for protein folding simula-tions[J]．Journal of Molecular Biology,1993,1(1):75-81[5]Kalegari D H,Lopes H S．A differential evolution approach for protein structure optimization using a 2D off-lattice model[J]．International Journal of Bio-Inspired Computation,2010,2(3/4):242-250[6]Wong K C,Wu C H,Peng C B,et al．Evolutionary multimodal optimization using the principle of locality[J]．Information Sciences,2012,4:138-170[7]Jooyoung L,Harold A S,Shalom R．New optimization methodfor conformational energy calculations on polypeptides:conformational space annealing[J]．Journal of Computational Chemistry,1997,8(9):1222-1232 [8] Klepeis J L,Pieja M J,Floudas C A．A new class of hybrid global optimization algorithms for peptide structure prediction:integrated hybrids[J]．Computer Physics Communications,2003(151):121-140 [9] Wong K C,Leung K S,Wong M H．Protein structure prediction on a lattice model via multimodal optimization techniques[C]∥GECCO’10Genetic and Evolutionary Coputation Conference．Portland,USA:ACM,2010:155-162 [10] ): Step1算法初始化:设置种群规模popSize,变异因子mutatorFactor,交叉因子crossFactor,初始种群POP=x1,x2,…,xpopSize｜xi=(xi1,xi2,…,xiN)i∈I,计算f2(xi),i∈I,并设f*2=mini∈If2(xi),I∈(1,2,…,popSize)。 Step2对种群进行能量极小化处理,并对其排序,选取种群前M个个体作为种子。 Step3 对于每个目标xi做以下处理: 3a)设置i=1； 3b)计算种群中非亲代和亲代的距离,并由大到小排序,再通过高斯转换函数将距离转化为轮盘形式。 3c)首先选取距亲代最近的个体作为变异的基准矢量xa；然后在轮盘中随机选取个体xb,xc,并保证a≠b≠c；最后对{xa,xb,xc}执行变异操作。 3d)为保证种群多样性,算法以不同的概率(m、n、p)从下列3种策略中选取一种执行交叉操作。 ①以概率m执行基本的DE交叉策略。 ②以概率n随机选取一个小组执行交叉操作。 ③以概率p随机选取一个集合组执行交叉操作。 3e)对所得的测试个体进行能量极小化处理。 3f)对所得的测试个体执行基本DE的选择操作,如果测试个体比亲代好,则替换亲代个体,否则保持种群不变。 Step4 如果i [11] 来鲁华．蛋白质结构预测与分子设计[M]．北京:北京出版社,1993 [12] Anfinsen C B．Principles that govern the folding of protein[J]．Science,1973,1(96):223-230 [13] 。 对模型(1)进行平移及旋转变换后,优化问题维数为3N－6,考虑到肽链的键长、键角均固定在平衡状态,设置Ebond=Eangle=Eother=0,可将维数降至3N－Ebond－Eangle－6。同时,为了消除转换过程中高阶的非线性等式约束条件,设rij＝ζ(τ1,τ2,…,τN),i,j=(1,2,…,N),i≠j带入式(1)可得到(具体转换过程参考文献 [14] ): min f2(τ)≡f2(i,ψi,ωi,χki) =f1(ζ(τ1,τ2,…,τN))｜Ebond=Eangle=Eother=0(2) subject to －π≤i≤π,i=1,…,NRES －π≤ψi≤π,i=1,…,NRES －π≤ωi≤π,i=1,…,NRES －π≤χki≤π,i=1,…,NRES,k=0,1,…,Ki 式中,τ=(τ1,τ2,…,τN)≡{i,ψi,ωi,χki |i=1,…,NRES,k=0,1,…,Ki}为肽链二面角向量；N为肽链中二面角的自由度,NRES表示肽链长度(或残基)个数,Ki为第i个残基侧链二面角的个数,且满足3NRES+∑NRESi=1Ki=N。 3 算法描述 3.1 基于空间局部原理的排挤差分进化算法 差分进化算法(DE)是1995年Storn等人提出的一种全局优化算法 [15] ,该算法采用基于种群的全局搜索策略和竞争生存策略对多个非劣解并行搜索,容易被改造成多模态混合优化算法。基本的DE算法是以变异的方式产生新个体,其数学表达式为: vt+1i,j=xta,j+mutatorFactor·(xtb,j－xtc,j)(3) 用差分进化这种群体优化算法处理多模态优化问题时,由于使用了全局选择因子,算法只能收敛到全局最优解,而忽略了众多局部极值解；其次,模型的复杂性造成这些算法极易陷入某个局优解；同时差分这种随机算法因缺乏全局收敛理论依据及解的不确定性而进一步限制了它们在实际问题中的应用。因此,Thomsen 2004年将排挤因子融合到差分进化算法中,得到排挤差分进化算法(CrowdingDE) |
| No related articles found! |
|
||
首页