Multiphase Simulated Annealing Based on Boltzmann and Bose-Einstein Distribution Applied to ...

16 Pages·2016·1.28 MB·English

Hindawi Publishing Corporation
Advances in Bioinformatics
Volume 2016, Article ID 7357123, 16 pages
http://dx.doi.org/10.1155/2016/7357123
Research Article
Multiphase Simulated Annealing Based on Boltzmann and
Bose-Einstein Distribution Applied to Protein Folding Problem
JuanFrausto-Solis,1ErnestoLiñán-García,2JuanPauloSánchez-Hernández,3
J.JavierGonzález-Barbosa,1CarlosGonzález-Flores,2andGuadalupeCastilla-Valdez1
1InstitutoTecnolo´gicodeCiudadMadero,Tecnolo´gicoNacionaldeMe´xico,AvenidaSorJuanaIne´sdelaCruzs/n,
ColonialosMangos,89440CiudadMadero,TAMPS,Mexico
2UniversidadAuto´nomadeCoahuila,CiudadUniversitaria,25280Arteaga,COAH,Mexico
3UPEMOR,BoulevardCuauhna´huac566,Jiutepec,62550MorMe´xico,CP,Mexico
CorrespondenceshouldbeaddressedtoJuanFrausto-Solis;[email protected]
Received24November2015;Revised5April2016;Accepted19April2016
AcademicEditor:F.Fdez-Riverola
Copyright©2016JuanFrausto-Solisetal. This is an open access article distributed under the Creative Commons Attribution
License,whichpermitsunrestricteduse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperly
cited.
A new hybrid Multiphase Simulated Annealing Algorithm using Boltzmann and Bose-Einstein distributions (MPSABBE) is
proposed.MPSABBEwasdesignedforsolvingtheProteinFoldingProblem(PFP)instances.Thisnewapproachhasfourphases:
(i)MultiquenchingPhase(MQP),(ii)BoltzmannAnnealingPhase(BAP),(iii)Bose-EinsteinAnnealingPhase(BEAP),and(iv)
DynamicalEquilibriumPhase(DEP).BAPandBEAParesimulatedannealingsearchingproceduresbasedonBoltzmannandBose-
Einsteindistributions,respectively.DEPisalsoasimulatedannealingsearchprocedure,whichisappliedatthefinaltemperature
ofthefourthphase,whichcanbeseenasasecondBose-Einsteinphase.MQPisasearchprocessthatrangesfromextremely
hightohightemperatures,applyingaveryfastcoolingprocess,andisnotveryrestrictivetoacceptnewsolutions.However,BAP
andBEAPrangefromhightolowandfromlowtoverylowtemperatures,respectively.Theyaremorerestrictiveforaccepting
newsolutions.DEPusesaparticularheuristictodetectthestochasticequilibriumbyapplyingaleastsquaresmethodduringits
execution.MPSABBEparametersaretunedwithananalyticalmethod,whichconsidersthemaximalandminimaldeterioration
ofprobleminstances.MPSABBEwastestedwithseveralinstancesofPFP,showingthattheuseofbothdistributionsisbetterthan
usingonlytheBoltzmanndistributionontheclassicalSA.
1.Introduction geneticalgorithms,andsomedeterministicmethods[3,4].In
fact,therearesomeexamples,suchasinsulinalphalytic[5,6]
IngeneticsDNA,RNA,andproteinsarethebasicelements with natural conformations whose energy is not minimal.
formanyresearches.DNAisamoleculethatcontainsgenetic This structure is usually named Native Structure (NS). In
instructions,whichareinvolvedinproteinsynthesisprocess addition,thefreeenergyofanNSconformationdependson
[1]. This molecule represents a complete set of heredi- theinteractionamongtheatomsandtheirrelativepositions.
tary information of any organism. DNA has four differ- ProteinFoldingProblem(PFP)isanenormouschallenge
ent nucleotides, which are adenine, cytosine, guanine, and and important problem in bioinformatics, medicine, and
thymine. This molecule is divided into genes, and a gene is otherareas[7].Thefunctionofaproteinisdirectlyrelatedto
asequenceofnucleotidesthatexpressaprotein.Afunctional its three-dimensionalstructure, and misfolded proteinscan
proteinisconformedinanapproximatedgeometricalmodel causeavarietyofdiseases.Theaimofthisproblemistofind
oftheglobalminimumenergy[2,3].Thisisadinamicprocess thenaturaltertiarystructureofaproteinusingonlyatarget
wherethelowestfreeenergyoftheproteinplusthesolvent sequence.Aproteincantakeahighnumberofdifferentcon-
canbereasonablyapproximatedbytheminimumfreeenergy formationalstructures fromits primarystructureto itsNS.
found by Monte Carlo, conformational space annealing, ThecomputationalprobleminvolvedtofindtheNSisknown
2 AdvancesinBioinformatics
as Protein Folding Problem. Because PFP is an NP-hard thatallconformationsofaprotein(exceptthenativestate)are
problem[8],heuristicmethodsavoidingthegenerationofall equallylikely.Thus,itismoreefficienttofindthenativestate
possible states of the protein are commonly used. In order byarandomsearch.PFPisaninterdisciplinaryproblemthat
tofindanNS,computationalmethodssearchstructureson involves molecular biology, biophysics, computational biol-
ahugespaceofpossiblesolutions.Thesemethodscanobtain ogy,andcomputerscience.Intheabinitiocase,NSprediction
severalstructuresveryclosetotheNS.Aparticularclassof requiresdifferentmechanismsthatleadthesearchingprocess
thesemethodsisknowntobeabinitiowhichlooksfortheNS toabiologicalthree-dimensionalstructure.Aswaspreviously
usingonlytheprotein’saminoacidsequence. mentioned, this process requires only the amino acids’
As a consequence, to solve PFP, new metaheuristics are sequence.PFPisanenormouschallengeandisveryhardto
applied, where simulated annealing (SA) [9, 10] is one of findtheNSofaproteinbecausethespaceofpossibleconfor-
the most successful [11–13]. Currently, classical SA applies mations of the protein is in general extremely large. For all
a Boltzmann distribution in order to accept bad solutions practicalpurposes,PFPcanbedefinedasfollows.
and escape from local minima. However, to generate high- Given
quality solutions for PFP, new and more efficient SA have
beendesigned;oneofthem,namedChaoticMultiquenching (i)asequenceof𝑛aminoacids𝑎1,𝑎2,...,𝑎𝑛 thatrepre-
Annealing Algorithm (CMQA), has obtained very good sentstheprimarystructureofaprotein,
results for proteins such as Met5-enkephalin, proinsulin, (ii)anenergyfunction𝑓∗(𝜎1,𝜎2,...,𝜎𝑛),wherethevari-
T0549,T0335,andT0281or1PLXW,1T0C,2K5E,SR384,and ables𝜎1,𝜎2,...,𝜎𝑛represent𝑛dihedralangles,
1A19, in PDB format, respectively. There are three central
phases of this algorithm [14]: (i) Multiquenching Phase findthefollowing:
(MQP), (ii) Annealing Phase (AP), and (iii) Dynamical
EquilibriumPhase(DEP).Allofthesephasesareexplained (i)theNativeStructuresuchthat𝑓∗(𝜎1,𝜎2,...,𝜎𝑛)rep-
inthepaper;forthisintroductionallweneedtoknowisthat resentsthelowestenergyvalue,where
eachphaseisdesignedwithanannealingapproachlooking (ii)thesolution𝜎∗ =𝜎1,𝜎2,...,𝜎𝑛definesthebestthree-
forfindingthebestconfigurationofthepreviousone.Atthe dimensionalconfiguration.
beginningoftheprocess,MQPimprovesarandomconfigu-
rationthroughanannealingprocedureexecutedatextremely Force fields are used to represent the energy of a protein;
very high temperatures; AP searches for a better solution someofthemostcommonareAMBER[17],CHARMM[18],
than that of MQP with an annealing search applied at high ECEPP/2[19–21],ECEPP/3[22],andGROMACS[23].These
temperatures,and,finally,DEPisappliedatlowtemperatures force fields compute energy components, for instance, the
lookingforabettersolutionthanthatobtainedbyAP.Asthe electrostatic energy, the torsion energy, the hydrogen bond
classicalSA,allofthesephasesapplyBoltzmanndistribution energy,andtheLennard-Jonesenergy.InthispaperECEPP/2
foracceptbadsolutions.However,Bose-Einsteindistribution forcefieldisused.
can also be used for escape from local minima [15]. Never- The atoms of a protein are represented in three-dimen-
theless,algorithmsusingthesetwodistributionsindifferent sionalcartesiancoordinates.Therearefourtypesoftorsion
rangesoftemperatureshavenotbeenpublishedforPFP. anglesordihedralanglesasfollows:
In this paper, a new SA algorithm named MPSABBE
(Multiphase Simulated Annealing based on Boltzmann and (i)The angle between the amino group and the alpha
Bose-Einstein distributions) is introduced. MPSABBE ap- carbonisreferredtoasPhi(𝜙).Thisanglerepresents
pliestheBoltzmannandBose-Einsteindistributionsathigh theanglebetweentheaminogroup(orNH2)ofthe
and low temperatures, respectively. The paper shows that aminoacid𝑖andthealphaCarbonC𝑖inthesequence;
usingbothdistributionsthequalitysolutionisimproved.This specifically, it represents the bond angle between N𝑖
paperisorganizedasfollows.InSection2,PFPisdescribed. atomofaminogroupandthecentralcarbon(𝛼C𝑖).
InSection3,theclassicalSAandMPSABBEalgorithmsare (ii)Thedihedralanglebetweenthealphacarbonandthe
explained. In Section 4, the SA applied for solving PFP is carboxylgroupisreferredtoasPsi(𝜓).Psirepresents
detailed. In Section 4, all the four MPSABBE’s phases are theanglebetweenthecarboxyl(COOH𝑖)groupofthe
presented. In Section 5, analytical tuning methods SA and aminoacid𝑖andthecentralcarbon𝑖(C𝑖)ofthesame
MPSABBEaredescribed.InSection6,experimentalresults amino acid. In particular, Psi measures the angle of
are shown. Finally, in Section 7, the conclusions of this thecovalentbondbetweenC𝑖 ofthecarboxylgroup
researcharediscussed. andthecentralcarbon(𝛼C𝑖).
(iii)Foreveryaminoacidssequence,anomegaangle(𝜔)
2.ProteinFoldingProblem is defined for each two consecutive amino acids 𝑖 −
1,𝑖; specifically, it is the angle of the covalent bond
PFPisrelatedtothequestionsofhowandwhyaproteinis
foldedintoitsNS.Theproteinsadoptanextremenumberof betweentheatomN𝑖ofaminoacid𝑖andcarbonC𝑖−1
ofthecarboxylgroupaminoacid𝑖−1.
possibleconformations[16],whichdependsonthenumberof
aminoacidsandthenumberofconformationsbyeachamino (iv)And, finally, each Chi angle (𝜒) is defined between
acid.TheessentialconceptintroducedbyLevinthalisthatthe thetwoplanesconformedbytwoconsecutivecarbon
PFP is a random search problem. This general idea means atomsintheradicalgroup.
AdvancesinBioinformatics 3
(1) Settinginitialtemperature(𝑇 )
initial
(2) Settingfinaltemperature (𝑇 )
final
(3) Generate𝑆 fromInitialSolution (𝑆 )
current initial
(4) 𝑇=𝑇
initial
(5) While(𝑇>𝑇 )do//TemperatureCycle
final
(6) While(stopcondition)//MetropolisCycle
(7) Generate𝑆 byapplyingaperturbationto𝑆
new current
(8) Obtaindifferencebetween𝑆 and𝑆
new current
(9) If(difference≤0)then
(10) Accept𝑆
new
(11) else
(12) BoltzmannProbability=exp(−difference/𝑇)
(13) If(BoltzmannProbability>random(0,1))then
(14) Accept𝑆
new
(15) endif
(16) endif
(17) endwhile
(18) Decrease𝑇byacoolingfunction
(19)endwhile
(20)Shownbettersolution(𝑆 )
better
Algorithm1:Pseudocodeofclassicalsimulatedannealing.
Thevariablesoftheproblemareallofthesefourangleswhich usuallyisveryclosetozero[9,10].Thereareseveralcooling
areintherange[0,360].Inthesimulationsconductedinthis functionsusedinSA[31–36],forexample,
researchwork,theseanglesaresetwithdiscretevalues.Some
variables have well-defined ranges as is the case of Psi and 𝑇𝑘+1 =𝛼𝑇𝑘 (1)
Pplhoita[n2g4l]e.sThwheoPsheiraannggleesiasrededfienfiendedinbtyhtehreaRngaemsa[c1h8a0n,d3r0a0n] 𝑇𝑘+1 =𝑒−𝛼𝑇𝑘 (2)
and[45,60].ThePsiangleisdefinedinthreeranges[20,180], 𝑇 = 𝑇𝑘 .
[300,330],and[180,205].Finally,theomegaangleisfixedat 𝑘+1 1+𝜂𝑇 (3)
𝑘
180degrees.
Themostcommonfunctionis(1).Thisfunctionreduces
thetemperatureparameterby𝛼factor,whichiscommonlyin
3.SimulatedAnnealingAlgorithm therangeof0.7≤𝛼<1.0.Aslowcoolingisappliedwhen𝛼is
verycloseto1,whileafastcoolingisappliedwhen𝛼isaround
3.1. Simulated Annealing Based on Boltzmann Distribution.
0.70.
Simulated Annealing (SA) Algorithm is a probabilistic
TheclassicalSAhastwocyclesasisshowninAlgorithm1;
methodproposedbyKirkpatricketal.[9]andCˇerny´[10]and
the first one is named temperature cycle and is used to
is an adaptation of the Metropolis algorithm, which is a
decreasethevalueofthetemperaturewithaspecificcooling
MonteCarlomethod[25].SAisbasedonthegradualmetal
function.Thesecondcycleisnamedmetropoliscycleandit
coolingforcrystallization.Thisalgorithmworksbyemulating
generates, accepts, or rejects solutions of the problem to be
the physical process where a metal is heating at very high optimized.Theinitialandfinaltemperaturevaluesareset(see
temperatureandthencooledveryslowlyuntilitsfrozenstate. lines(1)-(2)).Thesevaluesareobtainedbyananalytical(see
When this process happens, the metal is crystallized with Section5)orexperimentalway: 𝑇 shouldbeashighas
initial
thelowestenergyconfiguration.SAisanalgorithmthathas possible,while𝑇 shouldbeclosetozero.Aninitialsolu-
final
beenusedforfindingtheoptimalsolutionorclosetoitfor tion(𝑆 )isrequiredinSA;thissolutionisgenerated(see
initial
different NP-hard problems including biological problems line(3))andissetto𝑆 .Atthebeginningoftheprocess,
current
suchassequencealignment[26–28],phylogenetictrees[29], theparameter𝑇issetattheinitialtemperature(seeline(4)).
andPFP[30].Fromatheoreticalpointofview,SAconverges The temperature cycle is executed from 𝑇 to 𝑇 (see
initial final
totheoptimalsolutionorclosetothelowestfreeenergy[31]. lines(5)–(19)).Thenthemetropoliscycleisrepeated(seelines
However, classical SA is not able to find the lowest energy (6)–(17)) a certain number of times until a stop condition,
becauseenergybarriersaretoohighforSAandcannotescape which is explained later in this paper. A new solution
fromlocalminima.Asaconsequence,variantsofthismethod (𝑆 )isgeneratedwithinthemetropoliscyclebyapplyinga
new
areproposed[14,30]. small perturbation to the current solution 𝑆 (see line
current
Simulated annealing usually starts at a very high initial (7)). The difference between these two solutions (𝑆 and
new
temperature (𝑇 ). Through a cooling function, the tem- 𝑆 ) is calculated (see line (8)). In practice, SA can be
initial current
peraturevalueisgraduallyreducedfrom𝑇 to𝑇 ,which stoppedwhentheprobabilityofacceptinganewsolutionis
initial final
4 AdvancesinBioinformatics
negligible. For a minimization problem, if this difference is 3.3. Simulated Annealing Applied to Solve Protein Folding
less than or equal to zero (see line (9)), the new solution is Problem. The classical Simulated Annealing Algorithm can
accepted(seeline(10)).Whenthisdifferenceisgreaterthan beimplementedtosolvetheProteinFoldingProblem[40]as
zero, the Boltzmann distribution is applied. Then, a Boltz- isshowninthepseudocodeofAlgorithm2.Theinitialand
mann probability is calculated using (4) in line (12). If this finaltemperature(seelines(1)-(2))canbecalculatedaccord-
probabilityishigherthanarandomvaluebetween0and1(see ingtotheinstanceoftheproblembyapplyingtheanalytical
line (13)), then the new solution 𝑆 is accepted (see line methodparametersofSection5;thatmeansthattheprotein
new
(14)): shouldbepreprocessed.
𝑃(𝑆new)=𝑒(−Δ𝑆/𝑇). (4) 𝛼isArepqpulyirinedg.thTheecoteomlinpgerfuatnucrteiovnal(u1)e,tisherecdouocliendgvfearcytosrlovwalluye;
Afterthemetropoliscycleiscompleted,thetemperature thus,𝛼mustbeverycloseto1(seeline(3)).Inordertoreduce
value is reduced by a cooling function (see line (18)). For veryfastthetemperature,thecoolingfactor𝛼issetveryclose
a maximization problem, if the difference of 𝑆 − 𝑆 to0.70.AninitialsolutionofPFPiscreated,whichissetto
new current
is greater than zero, the new solution 𝑆 is accepted; else thecurrentsolution𝑆 (seeline(4)).Theinternalangles
new current
𝑆 canberejectedoraccepteddependingontheBoltzmann oftheinitialsolutionaremodifiedatrandom.Atthispoint,
new
probabilityvalue. thebestsolution𝑆 is𝑆 (seeline(5)).Theenergyof
better current
𝑆 iscalculatedbyapplyingaforcefieldfunction(seeline
current
3.2.SimulatedAnnealingBasedonBose-EinsteinDistribution. (6)).Beforestartingthetemperaturecycle,theinitialisloaded
Statistical Mechanics (SM) study the overall behavior of a into𝑇variableinline(7).Thenthetemperaturecyclestarts
systemconsistingofalargenumberofparticleswhosebehav- (see lines (8)–(26)) with a logic condition (𝑇 greater than
iorisunpredictable.SMusesstatisticsandprobabilitytheory 𝑇 inline(8)).Insideoftemperaturecycle,themetropolis
final
andthermodynamicprinciples.AccordingtoSM,theoccur- cycle is executed (see lines (9)–(24)). After this cycle is
rence of each future result is determined by a probabilistic completed,thevalueofthetemperatureisdecreased(seeline
functionsuchasBoltzmannandBose-Einsteindistributions. (25)).
Inaddition,onlythemostprobablebehaviorofthesystemin Inside the metropolis cycle, a new solution of Protein
thermalequilibriumatagiventemperatureisobserved[37]. FoldingProblem𝑆 isgeneratedbymodifyingtheprevious
new
Bose-Einstein distribution is obtained by finding the most solution𝑆 .Thisisdonebymodifyingtheinternalangles
current
probable distribution, that is, solving the problem defined of the protein (see line (10)). The energy of the protein is
bymaximizingthemostprobabledistribution,subjecttothe calculated(seeline(11)),andthedifferenceofenergies(i.e.,
followingconstraints:(ℎ1)thenumberofparticles(definedby between𝑆new and𝑆current)isdetermined(seeline(12)).This
thesummationofparticlesineachmicrostate)isconstantand difference is denoted by Δ𝑆 = 𝑆 − 𝑆 . The new
current new
(ℎ2)thetotalenergy(definedbythesummationofindividual solutionisacceptedwhenthenewsolutionisbetterthanthe
energiesofeachmicrostate)isconstant.Theproblemissolved previous one; thus, the current solution 𝑆 is replaced
current
using Lagrange multipliers. The parameters 𝜆 and 𝛽 are by 𝑆 (see line (14)). When a new solution is worse than
new
definedaslagragemultiplierofℎ1 andℎ2,respectively[38]. thecurrentsolution,itcanbeacceptedusingtheBoltzmann
ThentheBose-Einsteindistributionappliedforlowandvery distribution(seeline(21)).Theprobabilityofthisdistribution
lowtemperaturesisdefinedby (or acceptance probability)is directly related to the current
1 valueofthetemperatureandthedifferenceofenergybetween
ℎ(Δ𝐸)= . (5) 𝑆 and 𝑆 . This probability is calculated by (4). As
(𝑒𝜆+𝛽𝑒𝑖−1) new current
thetemperaturevalueisreduced,theacceptanceprobability
ThenparticlesbehaviorcanbemodeledbyBose-Einstein 𝑃(𝑆 )decreases.
new
distributiondefinedin(6).Thisequationdefinestheaccep-
tanceprobabilitydistributionofanewconfigurationofparti- 4.MPSABBE Algorithm
cles:
1 4.1. General Description. MPSABBE is a hybrid algorithm,
ℎ(Δ𝐸)= ,
(𝑒𝜆𝑒(Δ𝐸/𝐾𝑇)−1) (6) which has four phases (see Figure 1). These phases are (i)
MultiquenchingPhase (MQP) applied fromextremelyhigh
where 𝑇 is the temperature parameter, 𝜆 is related to the tohightemperatures,(ii)BoltzmannAnnealingPhase(BAP),
constraint of the total of particles in the system, and 𝐾 is whichisexecutedfromhightolowtemperatures,(iii)Bose-
theBoltzmannconstant.However,atveryhightemperatures Einstein Annealing Phase (BEAP) from low to very low
Bose-Einstein distribution practically becomes the Boltz- temperatures,andfinally(iv)DynamicalEquilibriumPhase
manndistribution.Nevertheless,atlowandverylowtemper- (DEP)whichappliesanannealingprocessatextremelylow
atures,theparticlesbehavedifferentlyandtheytendtocon- temperaturesusingBose-Einsteindistribution.
gregateatthesamelowestenergystate;theresultisknownas Inordertoacceptworsesolutions,BAPandBEAPapply
aBose-Einsteincondensate[39];asaconsequence,thesystem Boltzmann and Bose-Einstein distributions, respectively.
can be modeled by Bose-Einstein distribution. Section 4 Thisisdonewiththeaimofescapingfromlocalminima.DEP
presents a new SA applying both Boltzmann and Bose- is an extension of BEAP, where the stochastic equilibrium
Einstein distributions for accepting bad solutions for high is dynamically detected. This is done by using a regression
andlowtemperatures. method into the metropolis cycle; the iterations’ number is
AdvancesinBioinformatics 5
(1) Tuneinitialtemperature(𝑇 )byanalyticalmethod
initial
(2) Tunefinaltemperature(𝑇 )byanalyticalmethod
final
(3) Settingcoolingfactor(𝛼)
(4) 𝑆 iscreatedbymodifyingtheinternalanglesofprotein
current
(5) 𝑆 =𝑆
better current
(6) CalculateenergyofproteinapplyingaForceField
(7) 𝑇=𝑇
initial
(8) While(𝑇>𝑇 )do//TemperatureCycle
final
(9) While(stopcondition)//MetropolisCycle
(10) Createnewsolution(𝑆 )bymodifyinginternalanglesoftheprotein
new
(11) CalculateEnergyofproteinsusingaforcefieldfuntion
(12) Obtaindifferenceofenergiesbetweenthesetwoproteins
(13) If(difference≤0)then
(14) 𝑆 =𝑆
current new
(15) Ifenergy(𝑆 )>energy(𝑆 )then
new better
(16) 𝑆 =𝑆
better current
(17) endif
(18) else
(19) BoltzmannProbability=exp(−difference/𝑇)
(20) If (BoltzmannProbability>random(0,1))then
(21) 𝑆 =𝑆
current new
(22) endif
(23) endif
(24) endwhile
(25) Decrease𝑇bycoolingfunction(𝑇𝑘+1=𝛼∗𝑇𝑘)
(26)endwhile
(27)ShowthebestsolutionofPFP(𝑆 )
better
Algorithm2:PseudocodeSAappliedtoproteinfoldingproblem.
T
initial
High
temperatures
e
ur
at
er
p
m
Te
T
fMQP
Low
temperatures
T
Very low fBAP
temperatures T
fBEAP
T
f
Time
MQP BAP BEAP DEP
phase phase phase phase
Figure1:MPSABBEphases.
consideredastheindependentvariableandtheenergyvalue initial and final temperatures 𝑇initial and 𝑇𝑓 are determined
ofeachiterationasthedependentvariable.Theequilibrium using the analytical tuning method of Section 6. The other
detectioncriterionistheslopeoftheenergyfunctionintothe temperatures are determined using a variability criterion,
metropoliscycle.ThefourphasesMQP,BAP,BEAPandDEP suchasthevariabilitybeinglargerwherethetemperatureis
areexecutedinthetemperaturesrangeshowninTable1.The higher.
6 AdvancesinBioinformatics
Table1:TemperaturesrangesofMPSABBE.
Phase Initialtemperature Finaltemperature
MQP(fromveryhightohightemperatures) 𝑇initial 𝑇𝑓MQP
BAP(fromhightolowtemperatures) 𝑇𝑓MQP 𝑇𝑓BAP
BEAP(fromlowtoverylowtemperatures) 𝑇𝑓BAP 𝑇𝑓BEAP
DEP(fromverylowtoextremelylowtemperatures) 𝑇𝑓BEAP 𝑇𝑓
Tinitial In Algorithm 3, the MQP pseudocode of MPSABBE
is shown. At setting section (see lines (4)–(6)), the initial
temperatureiscalculatedbyananalyticalmethod.Thefinal
temperature of this phase (𝑇𝑓MQP) is set to an initial value,
determinedinanexperimentalway.Inlinefive,thevariable
ure 𝜏≈0 𝑇 is set to the initial temperature. The factors 𝛼Quenching
erat and 𝜏 are set to their initial values. The initial solution
mp 𝑆 isgenerated(seeline(8)).Theenergyofthissolution
Te Ecnuerrregnty(𝑆 )iscalculated,and𝐸(𝑆 )and𝑆 areset
current current min
toEnergy(𝑆 )and𝑆 ,respectively.
current current
T
fMQP Theexternalcycleisstartedatline(10),andthisisfinished
atline(31).ThisinternalcyclegeneratessolutionsofPFPand
acceptsorrejectssolutionsusingtheBoltzmanndistribution.
MQP The temperature parameter is decreased into this cycle by
phase applyingacoolingfunction(seeline(30)).Inthiscycle,𝜏is
Figure2:MQPphaseofMPSABBEalgorithm. set by (9) (see line (26)). When 𝜏 is very close to zero, this
variableissettoitsinitialvalue(seeline(28)).TheTempera-
turevalueiscalculatedby(7).
4.2.MQPPhaseofMPSABBE. MQPhasseveralsubphases. Aftertheexternalcycleisstarted,themetropoliscycleis
It starts at an extremely high initial temperature (𝑇 ), startedtoo.ThiscyclegeneratesnewsolutionsofPFP.Anew
which is obtained by an analytical method [41]. This ipnhitiaalse solution 𝑆new is obtained by applying a small perturbation
MisfiQnPisuhseedswthheecnoaoltihnrgesfhuonlcdtitoenmgpievreantubrye(𝑇𝑓MQP)isreached. tboettwheeecnutrhreenetnseorgluietisoonf𝑆𝑆ncuewrreanntd(s𝑆eceurlrienneti(s1c2a))lc.uThlateeddi(ffseereelnincee
(13)). If this difference is less than zero (see line (14)), then
𝑇𝑘+1 =𝛼Quenching𝛾𝑘𝑇𝑘, (7) (thseeenlienwes(1o5lu))t.io𝐸n(𝑆𝑆new is)aiscrceeppltaecde.d𝑆bcuyrr𝐸en(t𝑆isre)p(lsaeceeldinbey(1𝑆6n)e)w.
current new
where 𝛼 is a decrement factor of the temperature Ifthedifferenceofenergiesbetweenthesesolutionsislarger
Quenching
parameter,intherange[0.7,1.0],anddefineshowfasteach than zero, then the Boltzmann probability is applied (see
MQPsubphaseisdecreased.Averylow𝛼 valuewill line(17)).Ifthisprobabilityislargerthanarandomnumber
decreasethetemperatureveryfast.BesidesQ,u𝛾e𝑘ncihsindgefinedas between0and1(seeline(17)),thenthenewsolution𝑆new is
accepted (see line (18)). The 𝑆 is replaced by 𝑆 (see
current new
𝛾𝑘 =1−𝜏𝑘. (8) line(19)).If𝐸(𝑆current)islessthan𝐸(𝑆new)(seeline(21))then
𝑆 issetto𝑆 (seeline(22)).The𝐸(𝑆 )isreplacedby
min current min
The𝜏parameterisdefinedby(9),where0<𝜏<1,andit 𝐸(𝑆 )(seeline(23)).
current
definesaquadraticdecrementofthetemperature.Noticethat
𝜏convergestozeroand(7)isequivalentto(10):
4.3. BAP Phase of MPSABBE. In Algorithm 4, pseudocode
𝜏 =𝜏2 (9) ofBAPisshown.BAPisbasedonsimulatedannealing.The
𝑘 𝑘−1 temperatureparameterisdecreasedby(𝑇𝑘+1 = 𝛼Annealing𝑇𝑘)
𝑇𝑘+1 =𝛼Quenching𝑇𝑘. (10) or (𝑇𝑘+1 = 𝑒−𝛼Annealing𝑇𝑘). On the other hand, the length of
metropoliscycleisdeterminedby(21)or(27),respectively.In
The transition between two subphases is based on 𝜏 theinternalcycleoftheBAP,newsolutionsfortheinstance
parameter;itoccurswhen𝜏convergestozero(𝜏≈0).When aregenerated.Inthiscycle,abettersolutionthanaprevious
𝜏 is very close to zero, a new MQP subphase is started and oneisalwaysaccepted.However,worsesolutionsareaccepted
𝜏 is set to its initial value. This process continues until the orrejectedbyapplyingtheBoltzmanndistribution(4).The
temperature 𝑇𝑓MQP is reached. In Figure 2, the MQP phase length of the Markov chain (i.e., the internal cycle length)
isshown.Inthisphase,severalsubphasesareshown.Whena is determined by (21), where the increment 𝛽 is calculated
subphaseisstarted,theparameter𝜏issettoitsinitialvalue. with(22).Theinitialtemperaturewassettoathresholdvalue,
AdvancesinBioinformatics 7
(1) MQPProcedure()
(2) Begin
(3) //Settingsection
(4) 𝑇 =InitialTemperaturecalculatedbyanalyticalmethod
initial
(5) 𝑇𝑓MQP=Initialvalue,𝑇=𝑇initial
(6) 𝛼 =initialvalue,𝜏=Initialvalue
Quenching
(7) //Creationofinitialsolution
(8) 𝑆 =Createtheinitialsolution,𝐸(𝑆 )=Energy(𝑆 )
current current current
(9) 𝑆 =𝑆 ,𝐸(𝑆 )=𝐸(𝑆 )
min current min current
(10) Repeat//ExternalCycle
(11) Repeat//InternalCycle(MetropolisCycle)
(12) 𝑆 =Perturbation(𝑆 )//Uniformperturbation
new current
(13) Difference=𝐸(𝑆 )−𝐸(𝑆 )
new current
(14) IfDifference≤0Then
(15) 𝑆 =𝑆
current new
(16) 𝐸(𝑆 )=𝐸(𝑆 )
current new
(17) elseifexp(−Difference/𝑇)>random[0,1]Then//BoltzmannProbability
(18) 𝑆 =𝑆
current new
(19) 𝐸(𝑆 )=𝐸(𝑆 )
current new
(20) endif
(21) If𝐸(𝑆 )<𝐸(𝑆 )then//save𝑆
current min min
(22) 𝑆 =𝑆
min current
(23) 𝐸(𝑆 )=𝐸(𝑆 )
min current
(24) endif
(25) UntilMetropolisCycleisFinished
(26) 𝜏=𝜏∗𝜏
(27) If𝜏verycloseto0Then
(28) 𝜏=initialvalue
(29) endif
(30) 𝑇=𝛼 ∗(1−𝜏)∗𝑇
Quenching
(31) Until𝑇>𝑇𝑓MQP//ExternalCycle
(32)Endprocedure
Algorithm3:MQPpseudocodeofMPSABBEalgorithm.
which was the final temperature of MQP phase. The final whichbecomes
temperatureofBAPphaseisveryclosetozero.
𝑛 𝑛
𝑚=𝑘 ∑𝑖𝐸 −𝑘 ∑𝐸,
1 𝑖 2 𝑖 (12)
𝑖=1 𝑖=1
4.4.BEAPPhaseofMPSABBE. InAlgorithm5,pseudocode
of BEAP is shown. Again the external cycle decreases its
where
temperature value according to the cooling functions(1) or
(2). This time, the metropolis cycle length is constant, and 12
𝑘 = ,
it is equal to the maximum length of the last metropolis 1 𝑛3−𝑛
cycle in BAP phase. In this second cycle, the Bose-Einstein (13)
6
distributionisappliedforacceptingworsesolutions. 𝑘 = .
2 𝑛2+𝑛
4.5.DEPPhaseofMPSABBE. InAlgorithm6,theDEPgoalis Noticethatthecomplexityofthecomputationof(12)is
todetectthestochasticequilibriumbydeterminingtheitera- 𝑂(𝑛).Thisequationcontainsonlysummations;thus,itisless
tionwheretheslopeoftheenergyfunctionremainsveryclose complex than (11). These summations are computed using
tozero.Inordertodothat,letusdefinethenextvariables: simple data structures. 𝑘1 and 𝑘2 are only constants for a
(a) 𝑥𝑖 the number of the iterations in the metropolis cycle particular𝑛value.
(1,2,...,𝑛)and(b)𝐸𝑖theenergyfoundforthealgorithmin
iteration𝑥𝑖.Usingastandardleastsquaresmethod,theslope 5.AnalyticalTuningMethod
for𝑛iterationsisdefinedby
5.1. Parameters Setting Based on Boltzmann Distribution.
𝑛∑𝑛 𝑥𝐸 −(∑𝑛 𝑥)(∑𝑛 𝐸) ParametersofMPSABBEaretunedbytheanalyticalmethod
𝑚= 𝑖=1 𝑖 𝑖 𝑖=1 𝑖 𝑖=1 𝑖 , [42]. The initial temperature is defined by the maximum
𝑛∑𝑛 𝑥2−(∑𝑛 𝑥)2 (11) difference named maximum decrement Δ𝑍 , which is
𝑖=1 𝑖 𝑖=1 𝑖 max
8 AdvancesinBioinformatics
(1) BAPPhase()
(2) Begin
(3) 𝑇=𝑇𝑓MQP(FinaltemperatureofMQPPhase)
(4) 𝑇𝑓BAP=FinalTemperatureofthisphase
(5) 𝛼=initialvalue(veryclosetoone)
(6) 𝛽=Valuecalculatedbyanalyticalmethod
(7) CM=Initialvalue
(8) While(𝑇>𝑇𝑓BAP)do
(9) 𝑘=1
(10) while(𝑘≤CM)do
(11) 𝑆 =perturbationsystem(𝑆 )
new current
(12) Difference=𝐸 −𝐸
new current
(13) If(Difference≤0)then
(14) 𝑆 =𝑆
current new
(15) 𝐸(𝑆 )=𝐸(𝑆 )
current new
(16) ElseIf(exp(−Difference/𝑇))>random[0,1])then
(17) 𝑆 =𝑆
current new
(18) 𝐸(𝑆 )=𝐸(𝑆 )
current new
(19) Endif
(20) If𝐸(𝑆 )<𝐸(𝑆 )then//save𝑆
current min min
(21) 𝑆 =𝑆
min current
(22) 𝐸(𝑆 )=𝐸(𝑆 )
min current
(23) endif
(24) 𝑘=𝑘+1
(25) endwhile
(26) 𝑇=𝛼∗𝑇or𝑇=(exp(−𝛼))∗𝑇
(27) CM=𝛽∗CM
(28) Endwhile
(29)End
Algorithm4:PseudocodeofBAPphaseofMPSABBE.
(1) BEAPPhase()
(2) Begin
(3) 𝑇=Thresholdvalue;Determineexp(𝜆)
(4) 𝑇 =valueveryclosetozero
final
(5) 𝛼=initialvalue
(6) 𝛽=Valuecalculatedbyanalyticalmethod
(7) CM=Initialvalue
(8) While(𝑇>𝑇 )do
final
(9) 𝑘=1
(10) while(𝑘≤CM)do
(11) 𝑆𝑗=perturbationsystem(𝑆𝑖)
(12) Δ𝐸=𝐸𝑗−𝐸𝑖
(13) If(Δ𝐸≤0)then
(14) 𝑆𝑖=𝑆𝑗
(15) ElseIf((1/(exp(𝜆)∗exp(Δ𝐸/𝑇)−1)))>random[0,1])then
(16) 𝑆𝑖=𝑆𝑗
(17) Endif
(18) 𝑘=𝑘+1
(19) endwhile
(20) 𝑇=𝛼∗𝑇or𝑇=(exp(−𝛼))∗𝑇
(21) Endwhile
(22)End
Algorithm5:PseudocodeofBEAPphase.
AdvancesinBioinformatics 9
(1) DEPPhase()
(2) Begin
(3) While(𝑚>0)do
(4) 𝑛=1
(5) 𝐸 =0,Summary=0
summary
(6) while(𝑘≤CM)do
(7) 𝑆𝑗=perturbation system(𝑆𝑖)
(8) If𝐸(𝑆𝑗)=total clausulethenstop()
(9) Δ𝐸=𝐸𝑗−𝐸𝑖
(10) If(Δ𝐸≥0)then
(11) 𝑆𝑖=𝑆𝑗
(12) ElseIf((1/(exp(𝜆)∗exp(Δ𝐸/𝑇)−1))>random[0,1])then
(13) 𝑆𝑖=𝑆𝑗
(14) Endif
(15) 𝑛=𝑛+1
(16) Summary=Summary+𝑛∗𝐸𝑖
(17) 𝐸summary =𝐸summary+𝐸𝑖
(18) endwhile
(19) 𝑇=𝛼∗𝑇or𝑇=(exp(−𝛼))∗𝑇
(20) 𝑘1 =12/(𝑛∗𝑛∗𝑛−𝑛)
(21) 𝑘2=6/(𝑛∗𝑛+𝑛)
(22) 𝑚=𝑘1∗Summary+𝑘2∗𝐸summary
(23) Endwhile
(24)End
Algorithm6:PseudocodeofDEPphase.
calculatedusingasampleofrandomproteinstructuresatthe theminimumdeterioration.Theequationforcalculatingthe
highesttemperaturerange.Inthissample,theenergyoftwo finaltemperatureisdefinedby
consecutiveproteinstructuresdefinesasimpledecrementof
energy Δ𝑍𝑖,𝑗, and Δ𝑍max is the maximum difference in the 𝑇 = −Δ𝑍min . (16)
sample.Ontheotherhand,thefinaltemperatureiscalculated final ln(𝑃(Δ𝑍 ))
min
by applying the minimum deterioration (i.e., minimum
decrement)Δ𝑍 ofasampleofproteinstructurestakenat ThereareotherparametersofMPSABBEthatarecalculated
min
lowtemperatures.AnalyticaltuningbasedonBoltzmanndis- by applying a particular cooling function; for example, the
tributioncanbehelpfulforsettinguptheinitialtemperature. metropoliscyclelengthiscalculatedbyapplying
Theprobabilityofacceptinganynewsolution𝑆 isnearto
one(𝑃(𝑆 ) ≈ 1)athightemperatures,sothedneewcrementof 𝑇𝑘+1 =𝛼𝑇𝑘. (17)
new
thecostfunctionismaximal.Theinitialtemperature(𝑇 )
initial The analytical method determines the metropolis cycle
isassociatedwiththemaximumdeteriorationadmittedand
thedefinedacceptanceprobability𝑃(𝑆 ). length 𝐿𝑘 with a simple Markov model [42]; at high tem-
Let 𝑆 be the current solutionneawnd 𝑆 a new pro- peratures,onlyafewiterationsarerequiredbecause,inthis
posedonceu,rarenndt 𝑍(𝑆 )and𝑍(𝑆 )arethecnoewstsassociated condition, the stochastic equilibrium is reached very fast.
to 𝑆 and 𝑆 c,urrreenstpectively.nThewe maximum and mini- Nevertheless,atlowtemperatures,amoreexhaustiveexplo-
mumcurdreentterioratinoenwsareΔ𝑍 andΔ𝑍 ,respectively;then rationisneededand𝐿𝑘shouldbeaslargestaspossible.Let𝐿1
𝑃(Δ𝑍 )probabilityofaccmeapxtinganemwinsolution𝑆 with be𝐿𝑘atthetemperature𝑇initialandlet𝐿maxbethemaximum
themmaxaximumdeteriorationisdefinedby new metropoliscyclelength.Letthetemperature𝑇𝑘bedecreased
bythecoolingfunction(17)andlet𝐿𝑘+1becalculatedby
−Δ𝑍
𝑃(Δ𝑍)=exp( 𝑇 ). (14) 𝐿𝑘+1 =𝛽𝐿𝑘, (18)
ThisequationisbasicallytheBoltzmanndistribution,which where𝛽istheincrementcoefficientofmetropoliscycle(𝛽>
is applied for calculating 𝑇initial. This temperature value is 1),so𝐿𝑘+1 >𝐿𝑘and𝐿1istheinitialvalue.Themarkovchain
definedby lengthofthelastmetropoliscycleisequalto𝐿 .Functions
max
−Δ𝑍 (17)and(18)areconsecutivelyappliedinsimulatedannealing
𝑇initial = ln(𝑃(Δ𝑍max )). (15) from𝑇initial to𝑇final;consequently𝑇𝑛 and𝐿max areobtained
max by(19)and(20),respectively,
Similarly,thefinaltemperature(𝑇 )isestablishedaccord-
ingtotheprobabilityofacceptingfinaalnewsolution𝑆new with 𝑇𝑛 =𝛼𝑛𝑇initial (19)
10 AdvancesinBioinformatics
5
Table2:AverageresultsofMet -enkephalinwithMPSABBEalgorithm.
𝛼 Averageenergy(kcal/mol) Processingtime(minutes) AverageRMSD(A˚)
Annealing
0.75 −3.0836 0.1252 0.4517
0.80 −4.3025 0.1701 0.4327
0.85 −4.4093 0.2023 0.3510
0.90 −4.6493 0.3384 0.5097
0.95 −5.0634 0.8427 0.3610
𝐿max =𝛽𝑛𝐿1, (20) areavailable[42].Asiswellknowntheformercansimplybe
setclosetoone,whiletheseconddependsontheexploration
wheNreo𝑛ticisetthheatsttehpesinnucrmembeernftrocomeffi𝑇icniiteianltto𝛽𝑇cfiannalb.ecalculated lteiovnel,ethsteabvalilsuheedofin𝐿thedaelgpoernitdhsmo.nThtheeresifzoereo,ffonreiagnhyb𝑆o𝑖rhsooloud-
aifrethaevainiliatbialle.leAnsgtishw𝐿e1llaknndotwhnemthaexfiomrmumerlceanngtshimvaplluyeb𝐿emseaxt |𝑉𝑠𝑖|.Thus,𝐿max =m𝐶a|x𝑉𝑠𝑖|and𝐶 = −ln(𝑃𝑟(𝑆𝑖)),where𝑃𝑟(𝑆𝑖)
close to one, while the second depends on the exploration is the rejection probability for a solution 𝑆𝑖. The parameter
𝐶rangesfrom1to4.6;thelargervalueof𝐶assuresagood
levelestablishedinthealgorithmasfollows.
Thus,thenumberoftimesthatthemetropoliscycleisexe- explorationlevelintheneighborhoodof𝑆𝑖 atthefinaltem-
cutedcanbesimplyobtainedbyusing(21).Once𝑛isdeter- perature.Hence,differentexplorationlevelscanbeapplied.
mined the increment of the metropolis cycle length can be When we explore with 𝑃𝑟(𝑆𝑖) values of 63%, 86%, 95%, or
99%,theexplorationlevelsare𝐶=1,2,3,or4.6,respectively.
calculatedby(22):
Because𝐿 canbeverylargeforPFPinstances,itisimpor-
max
ln(𝑇 )−ln(𝑇 ) tanttoapplyaparticularprocessfordetectingthestochastic
𝑛= final initial (21)
ln(𝛼) equilibrium; this is done in DEP phase of MPSABBE that
detectsefficientlythestochasticequilibrium.Thenextsection
𝛽=exp(ln(𝐿max)−ln(𝐿1)). (22) explainsallMPSABBEphasesandtheperformanceofusing
𝑛 BoltzmannandBose-Einsteindistribution.
5.2. Parameters Setting Based on Bose-Einstein Distribution. 6.ExperimentalResults
The initial and final temperatures can be calculated by
applyingtheBose-Einsteindistribution.Then,theprobability MPSABBE is tested with five instances of PFP, which are
5
ofacceptinganewsolutionwiththemaximumdeterioration Met -enkephalin,proinsulin,T0549,T0335,andT0281.These
𝑃(Δ𝑍 ) is defined by (23). Consequently, the initial and instances have different sequence’s length and a different
max
finaltemperaturesarecalculatedwith(24)and(25),respec- numberofvariables(dihedralangles).Thesmallestsequence
5
tively, is Met -enkephalin, which has five amino acids and 19
variables. The largest sequence is a hypothetical protein
1
𝑃(Δ𝑍)= (23) (CASPT0281),whichhas90aminoacidsand458variables.
𝑒Δ𝑍/𝑇−1 Theproinsulininstancehas31aminoacidsand132variables;
Δ𝑍 the2K5E(CASPT0549)has73aminoacidsand343variables.
𝑇 = max (24) The instance Bacillus subtilis (CASP T0335) has 85 amino
initial ln((𝑃(Δ𝑍 )+1)/𝑃(Δ𝑍 ))
max max acidsand450variables.Thedihedralanglesusedinthesim-
Δ𝑍 ulationswerephi(Φ),psi(Ψ),omega(𝜔),andChi(𝜒).The
𝑇 = min .
(25) initialandfinaltemperaturearetunedanalytically.InMQP,
final ln((𝑃(Δ𝑍min)+1)/𝑃(Δ𝑍min)) parameters 𝛼 and 𝜏 are set with 0.85 and 0.999,
Quenching
Let𝑇𝑘bedecreasedbythecoolingfunction(2).Thus,𝑇𝑛 respectively.IneachsubphaseofMQPthefinalvalueof𝜏is
iscalculatedby setto0.001.
5
InTable2,theresultsofMet -enkephalinobtainedwith
𝑇𝑛 =𝑒−𝑛𝛼𝑇initial. (26) MPSABBE algorithm are shown. In this table, we show
the traditional average energy, processing time in minutes,
Asaconsequence,𝑛and𝛽arecalcultedby
and the average of the traditional RMSD (Root-Mean-
ln(𝑇 )−ln(𝑇 ) SquareDeviation)[43].TheRMSDwascalculatedusingTM-
𝑛= final initial (27) Align [44]. The best average solution for Met5-enkephalin
−𝛼
is−5.0634kcal/molwith0.8427minutesofprocessingtime,
𝛽=exp(ln(𝐿max)−ln(𝐿1)). (28) and the average RMSD obtained was 0.361A˚ (Angstroms).
𝑛 The RMSD is a measure which represents a structural
alignment between two proteins (target and solution). The
Noticethattheincrementcoefficient𝛽canbecalculated targetusedinthispaperwastakenfromProteinDataBank
if the initial and maximum metropolis length 𝐿1 and 𝐿max (PDB).AnRMSDneartozeroistakenasaperfectstructural