QuanticationofthedierencesbetweenquenchedandannealedaveragingforRNAs - PDF document

QuanticationofthedierencesbetweenquenchedandannealedaveragingforRNAsecondarystructuresTsunglinLiuandRalfBundschuhDepartmentofPhysics,OhioStateUniversity,191WWoodruAv.,ColumbusOH43210-1117Theanalyticalstudyofdisorderedsystemisusuallydicultduetothenecessitytoperformaquenchedaverageoverthedisorder.Thus,onemayresorttotheeasierannealedensembleasanapproximationtothequenchedsystem.InthestudyofRNAsecondarystructures,weexplicitlyquantifythedeviationofthisapproximationfromthequenchedensemblebylookingatthecorrela-tionsbetweenneighboringbases.Thisquantieddeviationthenallowsustoproposeaconstrainedannealedensemblewhichpredictsphysicalquantitiesmuchclosertotheresultsofthequenchedensemblewithoutbecomingtechnicallyintractable.PACSnumbers:87.14.Gg,87.15.v,05.70.FhI.INTRODUCTIONHeteropolymerfoldingisofcrucialsignicanceinmolecularbiology.Itisthebasisforthemechanismwithwhichcellscanproducethreedimensionalbuild-ingblocksoutoftheone-dimensionalinformationstoredintheirgenome.Cellsachievethisbyforming(stillone-dimensional)polymers(proteinsandRNA)bystringingtogetherdierentmonomerswithcovalentbonds.Allmonomersshareacompatiblebackbonebuttheyhavedierentsidechainsandoccurinapredenedorderalongthesequence.Physicalinteractionsbetweenthesemonomersforcethepolymertostablyfoldintoathreedimensionalstructure.Thisstructureiscrucialforthefunctionofthemolecule;itisdeterminedbythespecicsequenceofthepolymer[1{4].Inadditiontoitsbiologicalrelevance,heteropolymerfoldingisalsoaveryinterestingproblemofstatisticalmechanics[5{17].Thecompetitionbetweenthecongu-rationalentropyofthepolymer,theoveralltendencyofthemonomerstosticktoeachother,thesequencedisor-der,andthepreferenceoffoldingtowardabiologicallyactivenativestate,leadstoaveryrichthermodynamicphasediagram.WhilethesamequalitativebehaviorisexpectedforproteinsandRNA,wewillhereconcentrateonRNAsinceRNAfoldingismoreamenabletoanalyt-icalandnumericalapproachesthanproteinfolding.TherelativesimplicityoftheRNAfoldingproblemcomparedtotheproteinfoldingproblemdoesnotstemfromthefactthatRNAconsistsofonlyfourdierentbasesversusthetwentyaminoacidstheproteinsarecomposedof,butitcomesfromthesimplerinteractionrules:Thedomi-nantinteractionbetweenthefourbasesA,U,G,andCofanRNAmoleculeisWatson-Crick(G{CandA{U)pairformation,i.e.,iftwobaseshaveformedapairtheytorstorderdonottakepartinanyfurtherinteractions.Everyaminoacidofaproteinonthecontraryinteractswithallitsspatialneighbors,i.e.,withontheorderoftenotheraminoacidsatatime.Fromastatisticalphysicspointofview,thepossibilityofaglassphaseatlowtemperaturesdrivenbysequencedisorder,isofspecialinterestintheheteropolymerfold-ingproblem[6{12].Unfortunately,evenforthecaseofRNAfoldingananalyticquantitativedescriptionoftheglassphaseisstilloutstanding.Thus,quantitativestud-ieshavetoeitherrelyonnumericsortheyhavetousewhatisknownastheannealedaverage.Intheannealedaverage,thefreeenergyofthesystemisapproximatedasthelogarithmoftheensembleaveragedpartitionfunc-tion(insteadoftakingtheensembleaverageoverthelogarithmofthepartitionfunctioncalledthequenchedaverage).Physically,thisapproximationcorrespondstotreatingthesequencedegreesoffreedomasdynamicalinsteadoffrozenvariables.Thus,theannealedsystemrepresentsasequenceensemblethatiscoupledtothestructuralensemblebywayoftheinteractionenergies.Thissequenceensemblemaybedierentfromtheorigi-nalsequenceensembleofuncorrelatedrandomsequencesoverwhichthefreeenergyissupposedtobeaveraged.Duetothesedierencesbetweentheannealedandthequenchedsequenceensembletheannealedfreeenergyisonlyanapproximationtothetrue(quenched)freeenergyofadisorderedsystem.Thepurposeofthismanuscriptistorstquantifythedierencesbetweentheannealedandthequenchedse-quenceensembles.Specically,wewilllookatcorrelationbetweenneighboringbases.Weshowthatwhilethiscor-relationisstrictlyzerointhecorrect(quenched)sequenceensemble,theyarenon-zerointheannealedsequenceen-sembleandincreasewithdecreasingtemperature-uptocompletecorrelationincertainmodelsofRNAfolding.ThisclearlyunderlinesandquantiesthefundamentalshortcomingsoftheannealedaverageintheRNAfold-ingproblematlowtemperatures.Basedonthequantiednon-zeronearestneighborcor-relations,wethentrytodiminishthedierencesbe-tweentheannealedandquenchedensemblesbyforcingtheannealedensembletopresentzeroneighboringcor-relation.Thisconstrainedannealedensemblebehavesmuchmoresimilartothequenchedensemblethantheannealedensemble.Althoughtheglassphaseitselfcannotbeidentiedusingtheconstrainedannealedensem-blewhichonlypartiallycorrectstheoverallnon-randomcorrelations,onecanobtainthermodynamicquantities 2whicharemuchclosertothequenchedresultsthantheannealedonesusingthismethod.Thispaperisorganizedasfollows:InSec.II,webrie\ryreviewtheRNAsecondarystructureandintroducethegeneralRNAfoldingproblemwithsequencedisorder.InSec.III,wequantifythedeviationofnearestneighborcorrelationsoftheannealedensemble.Finally,weim-provethepureannealedensemblebyapplyingacon-straintofrandomcorrelationsinSec.IV.II.RNAFOLDINGPROBLEMWITHSEQUENCEDISORDERA.RNAsecondarystructuresRNAisasingle-strandedbiopolymeroffourdierentbasesA,U,C,andG.ThestrandcanfoldbackontoitselfandformhelicesconsistingofstacksofstableWatson-Crickpairs(AwithUorGwithC).ThiscomparativelysimpleinteractionschememakestheRNAfoldingprob-lemveryamenabletotheoreticalapproacheswithoutlos-ingtheoverall\ravorofthegeneralbiopolymerfoldingproblem[5].AnRNAsecondarystructureSischaracterizedbyitssetofWatson-Crickbasepairs(i;j)whereiandjdenotetheithandjthbaseoftheRNApolymerrespectively(conventionallyij).Here,wefollowmanypreviousstudies[5{17]andapplythereasonableapproximationtoexcludeso-calledpseudoknots[18],i.e.,fortwoWatson-Crickpairs(i;j)2Sand(k;l)2S,congurationswithikjlarenotallowed.Thisapproximationisjus-tied,becauseshortpseudoknotsdonotcontributemuchtotheoverallenergyandlongpseudoknotsarekineticallydiculttoform.B.QuenchedaveragingThepropertiesofRNAfolding,especiallythepossibil-ityofaglassphasedrivenbythesequencedisorder,havebeenachallengingproblemfromthestatisticalphysicspointsofview.Tounderstandthestatisticsofthisdis-orderedsystem,onersthastoassignanenergyE(;S)toeverysecondarystructureSforagivensequence.Thiscould,e.g.,simplybethenegativeofthetotalnum-berofWatson-Crickbasepairs.ThisthenallowsustocalculatethepartitionfunctionZ()=XS\(;S)eE(;S)=T(1)foragivensequencewhere\(;S)isonewhenthesecondarystructureSiscompatiblewiththesequenceandzerootherwise.Finally,onehastocalculatethequenchedaverageFq=kBThlnZ()i(2)overallsequences.C.AnnealedaveragingUnfortunately,thequenchedfreeenergyFqisverydif-culttocalculate.Thus,onecantrytoapproximatethequenchedfreeenergybythemucheasiercomputedan-nealedfreeenergy,whichtreatsthedisorderedsequencesasdynamicvariables.Thisannealedfreeenergyisonlyalowerboundofthequenchedfreeenergy,Fa=kBTlnhZ()iFq:(3)Itcanbequitedierentfromthequenchedfreeen-ergysincetheannealedensemblefavorsthosesequenceswheremorebindingpairsareallowed.Moreimportantly,physicalquantitiesderivedfromthisannealedfreeenergycanbeverydierentfromtheirquenchedcounterpartsaswewillshowexplicitlyinthefollowingsections.Tobespecic,wewillmeasurethecorrelationbetweenneigh-boringbaseswhichareknowntovanishinthequenchedcase.D.EnergymodelsInthispaper,westudythesimplestmodelofdisor-deredRNAsequenceswhichcontainonlythetwobasesAandU.Inassigningfreeenergiestosecondarystructures,weneglectanyloopentropiesandfocusonthebasepairsalone.Besides,formostpartsofthismanuscript,wedonotconsidertheminimalhairpinlengthconstraintwhichrequiresthetwobasesofabindingpairtobeseparatedbyatleastthreebasesinarealRNAmolecule.Withintheseapproximationswedoconsidertwodierentenergymodels.Inthebindingenergymodel,wesimplyassignanen-ergy=1toeachAU(orUA)bindingpair.WedenotethecorrespondingBoltzmannfactorbyq=e1=T.Thismodelcapturesthemainfeaturesoftheenergeticsandissimpleenoughforanalyticalandnumericalstudies.Wealsostudyasomewhatmorerealisticenergymodel,namelythestackingenergymodel.Inthismodel,weas-signenergiestothestackingoftwobasepairsratherthantoindividualbasepairs.Thisstackingenergydependsinrealityontheidentitiesofallfourbasesinvolved.WeimplementthiseectbyassociatingaBoltzmannfactors1withstackingsoftypesAAUUandUUAAwhileassociatingadierentBoltzmannfactors2withstackingsoftypesAUAandUAAU.Tobespecic,wewillchoosetheseBoltzmannfactorsass1=e2=Tands2=e1=Tfortheremainderofthiscommunication.Themainreasontostudythestackingenergymodelisthatthesimplebindingenergymodelisknowntobepathologicalwithoutaglassphaseatlowtemperatureinthedisorderedsequenceensemble[7{9].Asimplereasonisthatwhateverthesequence,eachbaseAcanalwaysndanotherbaseUtopairwithprovidedwehavethesameamountofbasesAandU.Thus,sequencesdisorder 3doesnotcausefrustration.Incontrast,theenergydis-tributionofthestackingenergymodelisgreatlyaectedbysequences,andastructureinwhichallbasepairsarestackedcaningeneralnotbefoundforeverysequence.Thus,sequencedisorderisexpectedtocausefrustration,andaglassphaseisexpectedinthisenergymodelforlowenoughtemperature.III.NEARESTNEIGHBORCORRELATIONSOFTHEANNEALEDENSEMBLEInthissection,wecalculatequantitativelyhowthenearestneighborcorrelationsintheannealedensembledeviatefromtheirtruevaluesintherandomsequenceensemble.Tothisend,wehavetocalculatetheannealedpartitionfunctionforsequenceswithlengthN-1,whichisdenedasZa(N)=12N1XS X\(;S)eE(;S)=T!:(4)1NNN-11Nk1Sk=1N-1AAUFIG.1:RecursiverelationexploringallpossiblesecondarystructuresforahomogeneoussequenceoflengthN.Thewavylinesstandsforcontributionfromallpossiblestructuresandsequences.Thestraightlinestandsfornon-pairedbases.Forthebindingenergymodel,thisannealedpartitionfunctioncanbeeasilyobtainedviatherecursiverelationshowninFig.1alongthelinesofpreviousstudies[10,19{22]buttakingthesequencesintoaccountexplicitly.Theideaistoseparatethetwocasesforthelastbase,whichiseitherunboundorboundtoacertainbasek,andthenrelatethepartitionfunctiontotheshorterlengthoneasZa(N+1;q)=Za(N;q)(5)+q2N1X=1Za(k;q)Za(Nk;q):Withthisrelation,onecanobtainananalyticalformulafortheannealedpartitionfunctioninthelargeNlimitbyperformingthez-transform,whichisdenedascZa(z;q)=1XN=1Za(N;q)zN;(6)ontherecursiverelation.AftersolvingtheresultingquadraticequationforcZa(z;q),wecanobtaintheparti-tionfunctionbydoingtheinversez-transform,Za(N;q)=12iIcZa(z;q)zN1dz:(7)Thisapproachcanbeeasilygeneralizedtothestackingenergymodel.Inordertokeeptrackofthecorrelationsbythean-nealedensemble,weassignanadditionalBoltzmannfac-torLtoallAAandUUneighborswithinthesequence.ThemodiedannealedpartitionfunctionisthenZa(N;q;L)=12N1XS X\(;S)qnq(S)LnL()!;(8)wherenq(S)isthenumberofbindingpairsinasecondarystructureS,andnL()isthenumberofconjugateneigh-bors,i.e.,AAandUUneighborsinthesequence.TheadditionalBoltzmannfactorcomplicatesthecal-culationofthepartitionfunctionsincedierentbasesAandUcontributedierently.However,wecanstillformu-laterecursiverelationsbynoticingthatthetwoendbasesofasequencepiecedeterminethecorrelationswithotherpieces.Thus,wecanseparateasequenceintotwocaseswheretheendbasesareeitherofthesametypeornot,andformulatetherecursiverelationforeachcaseinde-pendently.TheannealedpartitionfunctionZa(N;q;L)isthenobtainedviaz-transformasbefore.Sincethefor-mationoftherecursiverelationsisquitetechnical,weonlyaddresstheresulthere,anddeferthedetailstoAp-pendixA.Fromthepartitionfunctionwecanobtainthenear-estneighborcorrelationsbylookingatthedeviationoftheaveragedfractionofAU(orUA)neighborsfromtheexpectedvalue1/2inthedisorderedsequenceensemble.Thisdeviationisobtainedbytakingthederivativeas=121NL@Lln(Za(N;q;L))jL=1:(9)A.BindingenergymodelFig.2showstheneighborcorrelationsforthebindingenergymodel.Wendthatthedeviationmovesfurtherawayfromzeroastemperaturedecreases.Thisisadi-rectresultfromthefactthatatlowtemperature,themaincontributionstotheannealedpartitionfunctioncomefromthosesequenceswhichallowalotofbind-ingpairs,unlikethequenchedcasewheresequencesareequallyweighted.Theexactwaythattheneighborcorrelationsarebi-asedcanbeunderstoodasfollows.Inthisbindingenergymodel,theonlythingthatbiasesthenearestneighborcorrelationsistheformationofminimalhairpinssincetheyenforcetheneighboringbasestobedierent,whichareeitherAUorUA.Thus,thedegreeofbiasisdirectlycoupledtothefractionofsmallesthairpinsinasequence.Thisassertioncanbeveriedbystudyingthefractionofminimalhairpins.Asanexample,westudythezerotemperaturecasewhereallthebasesareexpectedtobepaired.Amongallpossiblepairingstructures,weexplic-itlycalculatethefractionofsmallesthairpins(withthe 400.511.522.53T00.0250.050.0750.10.125 dFIG.2:DeviationofthefractionofAU(orUA)nearestneigh-bors.Thedeviationisplottedasafunctionoftemperatureinunitsofthebindingenergyforthebindingenergymodel.Noticethatthedeviationmovesfurtherawayfromzeroandstopsataxedconstantastemperaturedecreases.Italsoapproachesalimitlargerthanzeroathightemperatureindi-catedbythedashedline.detailsshowninAppendixB).Asaresult,everyfourthbaseispartofaminimalsizehairpin.Thus,wehave1=4AU(orUA)nearestneighborsfromthesehairpinsandanother1=23=4=3=8fromtherestofthebasessincetheydonotshownearestneighborcorrelationbias.ThedeviationofthefractionofAU(orUA)neighborsisthenexpectedtobe5=81=2=1=8,whichmatchesexactlythezerotemperaturelimitinFig.2.Inthiscase,thesequence,asadynamicvariable,adjustsitselftoallthebindingpairs.Eveninthehightemperaturelimit,althoughallal-lowedsequencesareequallyweighted,therestillexistsanitefractionofminimalsizehairpinsonaverage.Asaresult,thedeviationofneighborcorrelationsapproachesaconstantlargerthanzero.Theassertionthatthedeviationiscoupledtotheformationofminimalsizehairpinsisagainveriedasweadditionallyrequireallthehairpinsbeingoflengthlargerthanone.Inthiscase,thecorrelationbetweennearestneighborsbecomesrandomatalltemperatures.How-ever,thesecondnearestneighborcorrelationsbecomenon-trivial.Thissimplebindingenergymodelgivesusatastehowthenearestneighborcorrelationsarecoupledwiththeen-ergythroughthestructure,i.e.,theformationofminimalhairpins.Thiscorrelationisbiasedsincetheannealeden-sembleputsmoreweightonlowerenergysequences.B.StackingenergymodelFollowingthesameapproach,wecheckthesamedevi-ationasafunctionoftemperatureinthemorerealisticstackingenergymodel.Again,onlytheresultisquotedhereinFig.3(interestedreaderscancheckthedetailedcalculationsinAppendixC).0.511.52T-0.5-0.4-0.3-0.2-0.100.1 dFIG.3:DeviationofthefractionofAU(orUA)nearestneigh-borsfortheenergymodelinvolvingstackingenergies.Unlikeinthecaseofthebindingenergymodel,theAU(orUA)neighborcorrelationsarecompletelybiasedatzerotempera-tureinthestackingenergymodel.Athightemperature,thisdeviationapproachesthesamelimitasthebindingenergymodel.Unlikethebindingenergymodel,atzerotemperature,thenearestneighborcorrelationsofthestackingenergymodelarecompletelybiased.AlmostnoAU(orUA)neighborscanbefoundinthisannealedsystem.Thiscanbeunderstoodsinceatzerotemperature,theonlydominatingstructureisalongsteminwhichallstack-ingloopsareoftypes1.Thus,theonlytwoimportantsequencesaretheonesmadeofhalfconsecutiveAbasesandtheotherhalfofUbases.Toverifythisstructure,weadditionallyintroducean-otherBoltzmannfactorhforeachhairpinloopformation.WiththisBoltzmannfactorwecankeeptrackofthefrac-tionofhairpinsfhintheannealedsystembycalculatingfh=1Nh@hln(Za(N;s1;s2;h;L=1))jh=1:(10)FromFig.4,wedoseethatthefractionofhairpinsofthisannealedsystemindeedgoestozeroastemperaturegoestozero,whichisafeatureofthelongstemstructure.Athightemperature,however,theenergymodeldoesnotmattersinceentropydominates.Thus,theAU(orUA)fractionapproachesthesamelimitasinthebindingenergymodel.Fromthisstackingenergymodel,welearnthatthestrongertheenergyiscoupledtothenearestneighborcorrelations,thelargerdeviationinnearestneighborcor-relationsoftheannealedsystemwillbepresentatlowtemperature.IV.CONSTRAINEDANNEALINGSofarwehaveonlyobservedthesequencecorrela-tionsarticiallyintroducedthroughtheannealedensem- 50.10.20.30.40.50.60.70.80.91T00.050.10.150.2fhannealedconstrained annealedquenchedFIG.4:Fractionofhairpinsinthestackingenergymodelforthreedierentensembles.ble.However,ourapproachcaninfactbeusedtogener-atemorerealisticensembleswithintheannealedframe-work.Theideaistoforcethenearestneighborcorrela-tionstoberandomwhenperformingtheannealedaver-age[23,24].Wesimplyenforcethisrandomdisorderconstraint,i.e.,thefractionofAU(orUA)neighborsbeingonehalfbysettingtheBoltzmannfactorL,whichcontrolsthenear-estneighborcorrelations,towhatevervalueitneedstohaveforthecorrelationsoftheannealedensembletovan-ish.Thisconstrainedannealingturnsouttopredictther-modynamicquantitiesmuchclosertothequenchedre-sults.Anditcanbedoneimmediatelyfollowingourquantieddeviationsindisorder.A.BindingenergymodelTheconstraintforthebindingenergymodelisreadas1NL@Lln(Za(N;q;L))jL=Lc=12:(11)Inthisenergymodel,weexpectthesequenceswithmoreAU(orUA)nearestneighborstobesuppressedsincetheannealedsystemfavorsthoseneighbors.Asaresult,Lc,whichfavorsAA(orUU)neighbors,isexpectedtobelargerthanoneinordertomeettheconstraint.Further-more,Lcshouldbelargeratlowertemperaturessincetheneighborcorrelationismorebiasedatlowertemper-atures.Oneimportantnoteisthattheresultingfreeenergyisonlydeneduptoanadditiveconstant,i.e.,addingaconstantbackgroundpotentialdoesnotchangethere-sultatall.Thus,theabsolutevalueofthisconstrainedannealedfreeenergyaswellastheBoltzmannfactorLchasnorealmeaning.Forexample,onecouldassigntheBoltzmannfactorLtoAU(orUA)neighborsinsteadofAA(orUU)neighbors.Theresultingchemicalpotentialwouldthenchangeasignandthefreeenergywoulddif-ferbyaconstantamount.However,thethermodynamicquantities,whicharecalculatedbytakingderivativesoftheconstrainedfreeenergy,willnotseethisconstantandareexpectedtobeclosertothequenchedresult.Toverifythisassertion,wearegoingtocomputetheaveragefractionofbindingpairsforthebindingenergymodelviaq=N@qln(Za(N;q;L))asafunctionoftem-perature.Then,wecomparethecasesoftheannealed(L=1),theconstrainedannealed(L=Lc)andthequenchedensembles.Astothequenchedresult,wenumericallycalculatethepartitionfunctiongivenrandomsequencesoflength1280andcollectthedatafrom1000randomsequences.Inor-dertoavoidthetrivialnitesizeeectsdueto\ructuationofthefractionofAbasesawayfromitsexpectedvalue1/2,weonlychoosesequencesthatcontainexactly640A'sand640U's.TheresultisshowninFig.5.Thesta-tisticalerrorsofthequenchedresultsarealwayssmallerthanthesizeofthecorrespondingsymbol,suchthatwithintheerrorbarsthequenchedresultsneveroverlapothercurves.Thisconditionholdsforallotherquenchedresultsinthismanuscript.0.20.30.40.50.60.70.80.91T0.340.360.380.40.420.440.460.48fqannealedconstrained annealedquenchedFIG.5:Fractionofbindingpairsinthebindingenergymodel.Theconstraintofrandomnearestneighborsbringstheannealedquantityclosertothequenchednumericalesti-mate.Thestatisticalerrorsofthequenchedresultsarealwayssmallerthanthesizeofthesymbol.Weseethattheconstrainedannealedresultisindeedveryclosetothequenchednumericalestimate.However,allthreeresultsareratherclosetoeachotheranyway.Thereasonforthesethreecasesbeingsoclosetoeachotherissimplythatunderthisenergymodelthesystemisnotglassy,andeverybaseisabletondanotherbaseforpairinginthisbindingenergymodel.Thus,atzerotemperature,allthebasesarepairedinallthreesystems.Thefactthatthenearestneighborcorrelationsarenotbi-asedalotcanalsobeveriedaswendthatatT=0.1,Lctobejust1.59.Thus,thechemicalpotentialintroducedfromtheconstraintiscomparativelysmallanddoesnotaecttheresulttoomuch. 6B.stackingenergymodelThesituationforthestackingenergymodelisverydierentfromthatofthebindingenergymodel.Here,wefollowthesameapproachandcomputetheaveragedfractionofstackingloopsoftypeAAUU(orUUAA)andAU(orUAUA)asafunctionoftemperatureundertheconstraint,1NL@Lln(Za(N;s1;s2;h=1;L))jL=Lc=12:(12)Similarly,inordertoavoidthetrivialnitesizeeectsforthequenchednumericalestimate,wexthenumberofAA,AU,UA,UUneighborsintherandomlychosensequencestobe320each[25].0.10.20.30.40.50.60.70.80.91T00.10.20.30.40.50.6fs1annealedconstrained annealed2 constrained annealedquenchedFIG.6:FractionofstackingloopsAAUU(orUUAA)inthestack-ingenergymodel.Theconstraintofrandomnearestneigh-borsxesthisquantitymuchbetterthanaveragednumberofpairsinthebindingenergymodel.Thephenomenolog-icalconstraint,i.e.,axedfractionofhairpins,bringsthisquantityonlyabitclosertothequenchedresult.FromFigs.6and7,weseethattheconstrainedan-nealedresultsaregreatlyimprovedovertheplainan-nealedresults.Thisveriestheideathatlargerdevia-tionsfromtherandomdisorderresultinabettercor-rectionviatheconstraintoftherandomdisorder.Forthisstackingenergymodel,atT=0.1,Lc=0.0067ismuchmoredierentfrom1thaninthebindingenergymodel.Fromtheseresults,wecanseethattheconstrainedannealedensembleofthestackingenergymodelbehavesinthefollowingway.SincetheensembleisforcedtohavethesamenumberofAA(orUU)andAU(orUA)neighbors,atzerotemperature,thedominatingstructureisstillalongstemstructure,butwithhalfthestackingloopsbeingoftypes1andtheotherhalfbeings2.ThisisconsistentwiththefactthatfractionofhairpinsgoingtozeroastemperaturegoestozerofortheconstrainedannealedsystemasshowninFig.4.Onedierencebetweenthequenchedensembleandtheconstrainedannealedensembleisthatnotallthebasesofarandomsequencecanformstackingloops.Thus,wehaveanitefractionofhairpinsinthequenchedensem-ble(Fig.4).Thisdierencecanusedasanadditional0.10.20.30.40.50.60.70.80.91T00.050.10.150.20.25fs2annealedconstrained annealed2 constrained annealedquenchedFIG.7:FractionofstackingloopsAUA(orUAAU)inthestack-ingenergymodel.Again,theconstraintofrandomnearestneighborsgreatlyimprovestheresult.However,unlikethecaseinFig.6,theconstraintofaxedfractionofhairpinsalsocontributesinbringingtheannealedquantityclosertothequenchedresult.phenomenologicalconstrainttoimprovetheconstrainedannealedsystemevenfurther.Weapplythisadditionalphenomenologicalconstraintbyrequiringthefractionofhairpinsfhtotthequenchednumericalestimatesandneighboringbasestobeuncor-relatedatthesametime,i.e.,toenforceLN@Lln(Za(N;s1;s2;h;L))jh=hc;L=Lc=12(13)hN@hln(Za(N;s1;s2;h;L))jh=hc;L=Lc=fh(T);(14)wherefh(T)isthequenchednumericalestimateinthisequation.FromFigs.6and7,weseethatthisadditionalcon-straintslightlyimprovethefractionofstackingloopss1,butsignicantlyimprovesthefractionofstackingloopss2.Thiscanbeunderstoodsincetheexistenceofhair-pinsintroducesAU(orUA)neighbors,ifthefractionofAU(orUA)neighborsisalsorequiredtobeonehalf,itwilldecreasethefractionofstackingloopss2amongthestemstructures.V.CONCLUSIONWeconcludethatthedeviationoftheannealeden-semblefromthequenchedensembleisstronglyrelatedtotheenergymodelandcanbecompletelybiasedwhenthecorrelationisstronglycoupledtotheenergyofthesystem.Quantifyingthisdeviationallowsustodocon-strainedannealingwhichbringsthepredictionsofther-modynamicquantitiesmuchclosertotherealvaluesinthequenchedensemble.Asthedeviationislarger,theconstraintisstrongerandthusbringstheannealeden-sembleevenclosertothequenchedresults.Unfortu-nately,thebiasingtowardthequenchedensembleisnot 7strongenoughtoactuallydrivethesystemintotheglasstransition.Besidesthenearestneighborcorrelations,onecouldalsoconsiderthecorrelationsfornextnearestneighborsoreventwobasesseparatedbyarbitrarydistances.Inprinciple,allthesecorrelationstogetherwouldbringustotheexactquenchedresultsandthustotheglasstransi-tion.However,thecalculationsbecomemuchmorecum-bersomeasoneincreasesthedistancebetweenthetwobases,andareleftforfuturework.VI.APPENDIXA.AnnealedpartitionfunctionforthebindingenergymodelTheannealedpartitionfunctionisobtainedbyrstsummingoverallcompatiblesequencesgivenasecondarystructureSandthensummingoverallpossiblestructuresS,whichcanbedoneviatherecursiverelationinFig.1.WedenetheannealedpartitionfunctionforasequenceoflengthNasZa(N+1).Inaddition,theannealedpartitionfunctionforasequenceoflengthNwithitstwoendbasespairedisdenedasAe(N1).TherecursiverelationinFig.1isthenreadasZa(N+1)=1+L2Za(N)+N1Xk=11+L4Za(k1)Ae(Nk):(15)Thefactor(1+L)=2forthersttermontherighthandsidecomesfromthecontributioninnearestneighborcor-relationsbetweenthefreebaseNandbaseN-1,andthe2takescareofaveragingoverthenumberofsequences.Wehaveasimilarfactorinthesecondtermcomingfromthecorrelationbetweenbasekofthearchandbasek-1.InthelaterpartwewillshowthatthebehavioroftheannealedpartitionfunctionismainlydeterminedbythearchtermAe,sowewillonlylookatthisquantityhere.TherstbaseofAeisalsospeciedtobeAandthelastbasetobetheconjugatebaseU.1NAUEE1NN1AUEE1kN1NAUSk=2N2FIG.8:Recursiverelationfortheannealedpartitionfunc-tionoverheterogeneoussequenceswheretherstandthelastbasesformaconjugatepair.Aletter'E'isusedtodenotethatthetwobasesattheendsofthearchareconjugatebases.Again,theannealedpartitionfunctionforthearchcanbeobtainedthroughasimilarrecursiverelation(Fig.8).Thetwotermsontherighthandsidearefurtherdecom-posedinFigs.9and10.Intheserelations,weneedtokeeptrackoftwofactors:theenergycontributionsandthenearestneighborcorre-E1NN1AU1EUN1N1N1NAOFIG.9:Decompositionofanarchwithitslastbaseinsidebeingafreebase,whichcanbeeitherAorU,intotwocases.Theletter'O'isusedtodenotethatthetwobasesattheendsofthearcharenon-conjugate.1NN1kEEAUk+11kOAAkN1EAUFIG.10:Separationofarches.lations.Fromtheenergeticpointofview,anarchcanbethoughtofsimplycontributingaBoltzmannfactorqandneednotstandforarealbindingpair,eventhoughinitiallyitisusedtorepresentarealbindingpair.Thus,inFig.9,aswetrytorelatetheannealedpartitionfunc-tiontoitsshorterlengthones,weassumeaneectivebindingpairbetweenbases1andN-1simplytoconservetheenergycontribution.Inthiscase,thetwobasesarenotreallypaired.Inordertokeeptrackofthecorrectnearestneigh-borcorrelations,weusealetterEonanarchtodenotethatthetwobasesattheendsofthearchareconju-gatebases.Similarly,aletterOisusedtorepresenttwonon-conjugatebasesattheendsofthearch.Thus,inFig.9,thetwocaseswherebaseN-1iseitherAorUareseparatedandaredenotedbyletterEandO,whichisdeterminedbywhetherthebases1andN-1areconjugateornot.Thesenotationsenableustoconnectthedecom-posedtermsrecursivelybacktotherelationinFig.8.InFig.10,aninnerarchcanbetreatedasafreebaseinconsideringtheenergyandcorrelationsfortherestofthebasesoutsidetheinnerarch.However,thereisadier-enceincountingneighborcorrelationsforthistreatmentbecausethefreebaselooksasabaseAfromtheright,butasabaseUfromtheleft.Thecorrectcorrelationscanbeobtainedifweshiftthisdiscrepancytothelastbaseand\ripitfromUtoA.Thus,thelasttermcarriesaletterOonthearchinsteadofE.TheserecursiverelationsarethenreadasAe(N1)=L2Ae(N2)+12Ao(N2)+(16)14N2X=2Ae(Nk1)[LAo(k1)+Ae(k1)];Ao(N1)=L2Ao(N2)+12Ae(N2)+(17)14N2X=2Ae(Nk1)[LAe(k1)+Ao(k1)]: 8Togetherwiththeinitialconditions,Ae(1)=q,Ae(2)=qL,Ao(1)=qL,Ao(2)=q(1+L)=2,onecansolveforAe(N)byperformingthez-transformcAe(z)=1XN=1Ae(N)zN;(18)cAo(z)=1XN=1Ao(N)zN;(19)ontherecursiverelations.AftersolvingforcAe(z),Ae(N)canbeobtainedthroughtheinversetransformAe(N)=12iIcAe(z)zN1dz:(20)Frompreviousstudies[10],weknowthatinthethermo-dynamiclimit,thepartitionfunctionhasananalyticalformasAe(N)/N3=2zc(q;L)N,wherezc(q;L)isthegreatestrealpartamongthebranchpointsobtainedfromthesolutionofcAe(z).Similarly,ifweperformz-transformonequation15,wecanrelatethez-transformoftheannealedpartitionfunctioncZa(z)tothatofthearchcAe(z).Sincethesetwosharethesamebranchpoints,theasymptoticbehavioroftheannealedpartitionfunctionisdierentfromtheaboveformulaforthearchbyjustadierentprefactor,whichdoesnotplayaroleinthethermodynamiclimit.ThefractionofAA(orUU)neighboringbasesperbaseoftheannealedsystemistheneasilycalculatedasL@Lln(zc(q;L))jL=1.Unfortunately,theanalyticalsolu-tionofthissetofpolynomialequationsistoocumber-sometoconveyanyusefulinformation.Thus,weresorttonumericalevaluationofthisanalyticalsolutioninthispaper.B.FractionofminimalhairpinsatzerotemperatureAsdiscussedinthemaintext,thefractionofminimalsizehairpinscanbeeasilyobtainedoncewegureoutthepartitionfunction.Atzerotemperature,thepartitionfunctionissimplerthanthenitetemperatureonesinceweonlyneedtoconsiderthegroundstateswhereallbasesarepaired.ThispartitionfunctionisobtainedthroughtherecursiverelationinasimilarwayasshowninFig.11.12NSN12k12N112N12NFIG.11:Recursiverelationforthepartitionfunctionwhereallthebasesarepaired.Wedenethepartitionfunctionforasequenceoflength2(N-1)asZm(N;h),wherehistheBoltzmannfactorforaminimalsizehairpin.TherecursiverelationisthenreadasZm(N+1)=N1Xk=1Zm(k)Zm(Nk+1)+hZm(N):(21)Togetherwiththeinitialconditions,Zm(1)=1andZm(2)=h,onecanobtaintheasymptoticbehaviorthroughz-transform.Aftersimplealgebra,wehavethelargestpolezc(h)=h+2ph+1forthez-transformofpartitionfunctioncZm(z;h).ThepartitionfunctionZm(N)isthenproportionaltozc(h)N.Thefractionofminimalsizehairpinspertwobasesistheneasilycalculatedas@hlnzc(h)jh=1=1+1=phh+2ph+1jh=1=1=2:(22)Thus,thefractionofminimalsizehairpinsperbaseis1/4.C.AnnealedpartitionfunctionforthestackingenergymodelThecalculationforthestackingenergymodelfollowsthesameapproache.However,itisabitmorecompli-catedsinceweneedtokeeptrackofstackingloopsin-volvingfourbaseswhichleadsustotherecursiverelationdepictedinFig.12.E1NN1AU1NAUES1AU2N1NAUUASSEE1kN1NAUSESk=2N2FIG.12:Recursiverelationforthestackingenergymodel.Intheserecursiverelations,weuseanadditionalletterSonthearchtodenotethefactthatweconsiderthestackingenergyofthestackingloopformedpartlybythatbindingpair.Independentofthetypeofthearch,allthestackingenergiesinsidethearchesarestillconsideredinallcases.Thus,thersttermontherighthandsideinFig.12doesnotcontainanSbecauseitsbaseN-1isunbound,andnostackingloopcanbeformedwiththebindingpairofthearch.E1NN1AU1EUN1N1N1NAO1N1AUES12N2N1AAUUAUESEESL1N1A12N2N1AAUUAESOSOSOAAFIG.13:Decompositionoftheannealedpartitionfunctionwhichlastbaseinsidethearchisafreebase.Similartotherecursiverelationinprevioussection,wethendiscardthelastbaseasafreebaseasshown 9inFig.13.Again,thearchesontherighthandsidearemeanttopreservetheenergycontributionsonly.Inthesecondlineoftherelation,wefurtherdecomposethetermsinordertorelatethesetermswiththerstrecur-siverelationinFig.12.ESAkN1E1kUESAN1N(1/h1)k+11kO1k+1OAAFIG.14:Separationofarchesconsideringthehairpincontri-bution.InFig.14,wealsoseparatethecontributionsoftheinnerarchfromtherestpartasinFig.10.Onedier-enceisthatweconsiderthecontributionfromthehairpinloopsinthiscase.Thus,thehairpinloopcontainedintheouterarchtermisnotarealhairpinloopbecauseoftheexistenceoftheinnerarch.Thecorrectresultisobtainedbyaddingthelasttermintherelation.Inthisstackingenergymodel,wedenotetheannealedpartitionfunctionforanarchoflengthN-1asAes(N).TherecursiverelationsarethenreadasfollowsAes(N1)=L2Aes(N2)(s11)(1+L2)4Aes(N4)+12Aos(N2)(s21)2L4Aes(N4)+s1L2+s24Aes(N3)(23)+14N2Xk=3Aes(Nk1)LAos(k1)(s11)2L4Aes(k3)+Aes(k1)(s21)1+L24Aes(k3)+2(1h1)Ho(k);Aos(N1)=L2Aos(N2)(s11)2L4Aes(N4)+12Aes(N2)(s21)1+L24Aes(N4)+s1L+s2L4Aes(N3)(24)+14N2Xk=3Aes(Nk1)LAes(k1)(s11)1+L24Aes(k3)+Aos(k1)(s21)2L4Aes(k3)+2(1h1)He(k);wherethetermsHeandHostandforthecontributionfromahairpin.TheyareobtainedseparatelyfromarecursiverelationsimilartotheoneinFig.9,byjustreplacingthewavylinebyastraightline,whichmeansthatbasesarenotbound.OnecantheneasilyformulatetherecursiverelationsforHeandHo.Togetherwiththeinitialconditions:Aes(1)=h,Aes(2)=hL,Aes(3)=h(1+3L2+s2+s1L2)=4,Aos(1)=hL,Aos(2)=h(1+L2)=2,Aos(3)=h(3L+L3+s1L+s2L)=4,wecanperformz-transformtoob-taintheasymptoticbehavioroftheannealedpartitionfunctionforthestackingenergymodel.[1]K.A.Dilletal.,ProteinSci.4,561(1995).[2]J.N.Onuchic,Z.Luthey-Schulten,andP.G.Wolynes,Annu.Rev.Phys.Chem.48,545(1997).[3]T.Garel,H.OrlandandE.Pitard,J.Phys.I(France)7,1201(1997).[4]E.I.Shakhnovich,Curr.Opin.Struct.Biol.7,29(1997).[5]P.G.Higgs,Q.Rev.BioPhys.33,199(2000).[6]P.G.Higgs,Phys.Rev.Lett.76,704(1996).[7]A.Pagnani,G.ParisiandF.Ricci-Tersenghi,Phys.Rev.Lett.84,2026(2000).[8]A.K.Hartmann,Phys.Rev.Lett.86,1382(2001).[9]A.Pagnani,G.ParisiandF.Ricci-Tersenghi,Phys.Rev.Lett.86,1383(2001).[10]R.BundschuhandT.Hwa,Phys.Rev.E65,031903(2002).[11]F.Krzakala,M.MezardandM.Muller,EuroPhys.Lett.57,752(2002)[12]E.Marinari,A.PagnaniandF.Ricci-Tersenghi,Phys.Rev.E.65,041919(2002)[13]M.Muller,F.KrzakalaandM.MezardEuro.Phys.J.E.9,67(2002).[14]H.Orland,andA.Zee,Nucl.Phys.B.620,456(2002)[15]R.Mukhopadhyay,E.Emberly,C.TangandN.S.Wingreen,Phys.Rev.E.68,041904(2003)[16]M.Baiesi,E.OrlandiniandA.L.Stella,Phys.Rev.Lett.91,198102(2003)[17]P.LeoniandC.Vanderzande,Phys.Rev.E.68,051904(2003)[18]I.TinocoJr.andC.Bustamante,J.Mol.Biol.293,271(1999),andreferencestherein.[19]P.G.deGennes,Biopolymers6,175(1968).[20]M.S.Waterman,AdvancesinMathematics,Supplemen-tarystudies,editedbyG.-C.Rota(Academic,NewYork,1978),pp.167-212.[21]J.S.McCaskill,Biopolymers29,1105(1990)[22]M.ZukerandD.Sanko,Bull.Math.Biol.46,591(1984)[23]T.Morita,J.Math.Phys.5,1401(1964)[24]E.Orlandini,A.RechnitzerandS.G.Whittington,J.Phys.A.35,7729(2002)[25]S.F.AltschulandB.W.Erickson,Mol.Biol.Evol.2,526(1985)