Table Of Content

Freezing andfoldingbehavior insimpleoff-latticeheteropolymers J. E. Magee,1 J. Warwicker,2 and L. Lue3 1Department of Chemical Engineering 2Department of Biomolecular Sciences 3Department of Chemical Engineering, UMIST, PO Box 88, Manchester, M60 1QD, UK (Dated:January11,2010) We have performed parallel tempering Monte Carlo simulations using a simple continuum heteropolymer modelforproteins. Alltenheteropolymersequenceswhichwehavestudiedhaveshownfirst-ordertransitions at low temperature to ordered states dominated by single chain conformations. These results are in contrast 4 withthetheoreticalpredictionsoftherandomenergymodelforheteropolymers,fromwhichwewouldexpect 0 0 continuoustransitionstoglassybehavioratlowtemperatures. 2 n I. INTRODUCTION a J 6 Manybiologicalheteropolymers(inparticularproteins)canadoptasingle,sequencedependentconformation(thefold,orna- 1 tivestate)underphysiologicalconditions.Thisspecificityallowsgroupsofatomswithinsuchamoleculetobeplacedaccurately inspace,allowingthevastrangeofprecisebiochemicalfunctionswithinlivingorganisms.Uponalteringtheenvironmentaround ] theprotein(e.g.pH,solventquality,temperature,etc.),themoleculecanbemadetolosethissingleconformation(i.e.unfoldor t f denature),adoptingafluctuating,disorderedstate. Thisfoldingtransitionhasbeenshowntobebothreversibleandfirst-order o formanyproteins[1],andisstillpoorlyunderstood. s . The modern theoretical viewpoint on the thermodynamics of protein folding (recently reviewed by Pande et al. [2]) uses t a conceptsfromthetheoryofspinglasses. ThisapproachwasintroducedbyBryngelsonandWolynes[3],whousedtherandom- m energymodel(REM)[4]todescribethelowtemperaturebehaviorofheteropolymers. Thissimplemodelassumestheenergies - ofthepossibleconformationsofaheteropolymertobeindependentanduncorrelated,andpredictsasecondordertransition(no d latentheat)toaglassystatedominatedbyalowenergyconformation. ThemorecomplicatedtheoryofShakhnovitchetal.[5] n usesanothertrickfromthespinglassrepertoire,thereplicaapproach,toarriveatthesameconclusions. o c The low energy conformation occupied in the glassy state is not necessarily the native state. Depending upon details of [ conditions, sequenceandsample history, the heteropolymermay well become“misfolded”; thatis, lockedin to a non-native, lowenergyconformation. Indeed,mostrandomlyselectedsequencesshouldhaveseveralclosely-spacedconformationsatthe 1 low end of their energy spectrum, and will not reliably or reproducibly fold into only one of these. Some annealing of the v 1 heteropolymersequence(i.e.designwork)shouldbeneededto“drivedown”theenergyofachosen(target)conformationwell 9 belowtheclosely-spacedregionoftheenergyspectrum,preventingsamplingofcompetingconformationsandallowingrobust 2 folding[2]. Simulations[2, 6]andnumericalstudies[7]usinglatticeheteropolymermodelssupportthesepredictionsandthe 1 approximationsbehindthemodels. 0 However,therandomenergyapproximationisessentiallymeanfieldinnature,andcannotaccuratelytreatsystemswithper- 4 sistent correlations [8]. Further, lattice models artificially constrain the configurations available to a polymer. While lattice 0 homopolymerresultshavebeeninstructiveathightemperatures[9,10],theydonotpredictthelowtemperaturefreezingtran- / at sition seenin thesimulationsofZhouetal. [11, 12], whereanisolated homopolymerfreezesinto a single, apparentlyunique m crystallinestateviaafirst-ordertransition. We have performedMonte Carlo simulationsusing ten randomsequence realizations of a simple off-lattice heteropolymer - d modeltotestthesepredictionsincontinuumsystems. Wereportthat,ratherthanexhibitingglassybehavioratlowtemperature, n alltensequencesstudiedpassthroughafirst-ordertransitiontoanorderedlowenergystate,dominatedbyasingleconformation, o withoutneedforsequencedesignwork. c : v Xi II. MODEL r a Our heteropolymersare modeledas freely-jointedchainsofN hard-spheremonomersof diameters , each assigneda type; we denote the type of monomer i by s. The model is a simple extension of that used by Zhou et al. [11, 12] and is shown i schematicallyinFig.1. Bondedmonomersinteractviathepotential ¥ r < (1−d )s u (r)= 0 (1−d )s < r < (1+d )s (1) i,i+1 ¥ (1+d )s < r  The parameter d is a bond-lengthflexibility parameter; this allows for conversionto molecular dynamics simulation and has beenshowntomakenophenomenologicaldifferencetohomopolymerfoldingresults[11,12]. Weused =0.1,allowingaten percentvariationinbondlength. Ifmonomersiand jarenotdirectlybonded(thatis,|i−j|>1),thentheyinteractviaasquare-wellpotential: ¥ r < s uij(r)= es0isj lss << rr < ls (2)  wherer denotesthe distancebetweenmonomersiand j, ls isthesquarewellwidth,ande is aninteractionmatrix giving ss′ square-wellinteractionstrengths(square-welldepths)betweenthemonomertypes.Inthiswork,weuseamodifiedhydrophobic- polar,ormHPmatrix[36]. Thisisgivenby: −1 −1 e mHP=e ′ 2 (3) −1 0 (cid:18) 2 (cid:19) wheree ′isthecharacteristicenergyscale.Withthismodel,monomerswiths =0arehydrophobic-like(H-type),preferringtype- i symmetric interactions, while type s =1 monomersare polar-like (P-type), preferringtype-asymmetricinteractions. Similar i interactionmatricesarewidelyusedinon-latticeproteinfoldingsimulationsandnumericalwork[13,14,15]andarebelieved to be good approximations to real intraprotein interactions [16]. In the dense, globular heteropolymer state, we expect this interactionmatrix to lead to configurationswith a core of hydrophobicmonomerssurroundedby polarmonomers, as seen in globularproteins. Weusereducedunitsscaledbythehard-sphereradiuss andtheenergyscalee ′,suchthat r∗ = r/s E∗ = E/e ′ (4) T∗ = k T/e ′ B Polypeptidesoflengthgreaterthanapproximately30aminoacidscanformstablefoldeddomains(e.g.the35aminoacidWW domain[17]). WehavechosentostudyN=64lengthheteropolymers;thislengthisshortenoughforsimulationinreasonable time,whilelongenoughthatequivalentrealproteinscanfold. WehavestudiedtenrandomsequencesinteractingviathemHP interaction matrix, each with a 1:1 H:P composition ratio. If we consider P-type monomersto map onto charged and polar amino acid residues, and H-type monomers to map to all other amino acid residues, this composition ratio is similar to that foundinrealproteins[18]. III. SIMULATIONMETHOD Monte Carlo simulations were performed for a single polymer in isolation. As such, we have not used periodic boundary conditions.Wehaveimplementedcrankshaft,pivotandcontinuumconfigurationalbias(CCB)polymermoves[19].Bondlength fluctuation is allowed both in the CCB moveimplementation, as well as throughstandard Monte Carlo particle displacement moves performedon the monomers. CCB moves regrow end segments of up to four monomers. Moves are performed with relativeprobabilitiesintheratioN:(1/2):(1/4):(1/4)forparticledisplacement,crankshaft,CCBandpivotmovesrespectively. Initialequilibrationrunslast3×109MonteCarlomoveattempts;datacollectionrunslast1×109moveattempts. TheMonteCarlo simulationsalso useparalleltempering(also knownasreplicaexchange)[20]. SeveralMonteCarlo sim- ulationsofthesamesystematdifferingtemperaturesareruninparallel; thesesimulationsattempttoexchangeconfigurations at regular intervals. This helps to maintain the ergodicity of the simulations; an individual system which has become stuck at low temperature in a metastable “trap” may be promoted to higher temperatures, where it can explore a wider range of phasespace,aidingitsequilibration. Ourparalleltemperingsetupusestenindividualsimulationsacrossthetemperaturerange T∗=0.15...1.5,distributedsuchthattheratiobetweenconsecutivetemperaturesisconstant[21]. Configurationswapmoves areattemptedevery1000MonteCarlomoveattempts. The simulations sample the configurational energy E, the radius of gyration r , the two-particle density r (2)(r), and the g instantaneouscontactmapcforthesystemevery105MonteCarlomoveattempts. Thetwo-particledensityisdefinedas r (2)(r)= 1 (cid:229) (cid:229) d (r-r ) (5) ij N * i j6=i + wherer isthedistancebetweenmonomersiand j[22]. ij The instantaneouscontactmap c reflectsa heteropolymerconformation;it is an N×N matrix, with each elementc equal ij tooneifmonomersiand j arewithineachotherssquarewell,andzerootherwise. Theinstantaneouscontactmapistherefore givenby: 0 r <ls c = ij (6) ij 1 r >ls ij (cid:26) Theaveragecontactmap c givesusinformationonthestructureadoptedbythesystematgiventemperature. Wecanuse ij ittocalculateanormalizedthermalaverageoverlap: (cid:10) (cid:11) (cid:229) 2 c ij hQni= i6=(cid:229) j(cid:10) (cid:11) (7) c ij i6=j (cid:10) (cid:11) This normalized overlap parameter can take values between zero and one. Values close to one indicate that the system is dominated by only a very small number of very similar conformations; small values indicate that the system is exploring a varietyofdifferentconformations. Energyandradiusofgyrationdatafromtheindividualsimulationswithinaparalleltemperingrunaretiedtogetherusingself consistentmultiplehistogramextrapolation[23]. Thisallowsthe histogramsofenergyE andotherobservables{A}recorded fromtheindividualsimulationsto be combinedto estimate W (E,{A}), the densityofstates ofthe system. Fromthis, we can interpolatethe behaviorof the system across, and to a certain extentbeyond,the simulated temperatureand parameterrange. Weusetheextrapolateddistributionatgiventemperaturetocalculatetheconstantvolumespecificheatc∗,givenby v E∗2 −hE∗i2 c∗= (8) v T∗2 (cid:10) (cid:11) Histogramextrapolationmethodsrelyongoodsamplingacrosstherangeofinterestbytheindividualsimulations. Assuch, simulationsaroundfirst-ordertransitionswithasignificantfreeenergybarrierbetweenphasescanbeproblematic.Ifsimulations on either side of the transition temperaturesample onlythe phase which is moststable at that temperature,but no simulation actually observes the system passing back and forth between phases, the extrapolation process can “blur out” the transition, leavingasignatureonlyinanincreaseintheerrorestimateandtheresponsefunctions(e.g.specificheat)acrosstheregionof thetransition.Worsestill,ifthesystemexhibitsstrongmetastability,thiscanintroducealargesystematicerrorintoanyestimate forthepositionofthetransition. Toavoidthisproblem,we haveusedmulticanonicalsampling[24]to directlyconfirmordenyphasecoexistencewherewe suspecttheexistenceofafirst-ordertransition. Wehaveusedaniterativemulticanonicalscheme[25]. RepeatedbiasedMonte Carlosimulationsrefineourestimateforthedensityofstatesofthesystem, withtheaimofsamplingevenlyinenergyacross thesuspectedrangeofcoexistence.Thebiascanthenberemovedfromtheresultsofthesesimulations,andtherefinedestimate forthedensityofstatescanbesplicedbackintothewidersimulationresults. Theindividualsimulationsusedinthisprocedure areperformedusingthedatacollectionrunparametersstatedabove. IV. SIMULATIONRESULTS Plots of E∗, c∗ and r∗ for our sequences, interpolatedacross the studied temperature range, are shown in Fig. 2. The phe- v g nomenologydoes notvary considerablywith sequence. Above a temperatureof T∗ ≈0.9, the heteropolymerexists in a dis- orderedcoilstate -a highenergystructurethatrapidlyincreasesitsradiusofgyrationanddecreasesits heatcapacity(energy fluctuations)withtemperature. Between T∗ ≈0.3andT∗ ≈0.9, wesee theexpectedlowtemperaturecollapsedglobulestate with a small radius of gyration (order of 0.05s ) and with relatively high heat capacity (order of 2.0k ). Around T∗ ≈0.26, B however, we see sharp spikes in the specific heat, associated with sudden drops in both conformationalenergyand radius of gyrationwithdecreasingtemperature,suggestiveoffirst-ordertransitions.Wehaveusedmulticanonicalsamplingtostudythese apparenttransitions;forallsequences,wehavefoundtemperaturesatwhichdoublepeakedprobabilitydistributionsoccurfor conformationalenergy,indicatingtwo-phasecoexistence. TheseprobabilitydistributionsareshowninFig.3(a). Thesequence dependenceoftheprobabilitydistributionsissmall. In Fig. 3 we plot coexistence temperatures for the sequences studied against the energy difference D E∗/T∗ between the coexisting phases. The quantitiesare clustered around T∗ ≈0.27 and D E/k T ≈56, but otherwise appear uncorrelated. For B comparison,thetypicallatentheatD H ofthefoldingtransitionforaglobularproteinisoforderof100k T [26]. B A. Structuralcharacterofthelowenergyphase We presenttwo-particledensitiesr (2) data, calculatedusingEq.(5), for sequence4 in Fig. 4; these plotsare typicalof the resultsforalltensequencesstudied. Thesedatawerecollectedbeforethemulticanonicalalgorithmwasused,andmayexhibit somedegreeofmetastability. Wefirstconsiderr (2) (r∗),thehydrophobic-hydrophobictwo-particledensities,inFig.4(a). We HH see thatbelow the transitiontemperature,thefirst threeneighborshells areclear anddistinct, separatedby regionswith near- zero hydrophobicparticledensity. The nextthreeneighborshells are much moresparsely populated(limited bythe diameter of the collapsed state), but remain distinct, separated by sharp minima. These data show a low energy phase structure with localtranslationalorderbetweenhydrophobicmonomers. Forr (2)(r∗)abovethetransition,thefirstthreeneighborshellsare HH still clear up to T∗ =0.7, but not distinctly separated. The next three neighbor shells have blurred together into one single peak. Above the transition temperature, the hydrophobicmonomers have no translational order, as expected in a disordered globularphase.Weshowtypicallow-energyphaseandglobularphaseconfigurationsinFig.5;theorientationofthelow-energy configuration has been chosen to clearly show the hydrophobicmonomersarranged in hexagonalcrystal planes, whereas no crystalstructureexistsintheglobularconfiguration. Figures4(b)and(c)showhydrophobic-polarandpolar-polartwo-particledensitiesr (2)(r∗)andr (2)(r∗),respectively. Nei- HP PP thershowanysignoftranslationalorder,eitheraboveorbelowthetransition. Thepolarmonomersremaininadisorderedstate acrossthetemperaturerangestudied. WeshowthethermalaveragecontactmapsinFig.6. Notethatthesedataarecollectedfromequilibratedparalleltempering runsbeforethemulticanonicalalgorithmisused,andassuchexhibitmetastability. Atlowtemperatures(belowT∗≈0.32),we seethatthecontactmapsaredominatedbyverypersistentcontacts,suggestingthatthelowtemperaturestructureisdominated byveryfewconfigurations.Contactmapscollectedathighertemperaturesshowlargelowpersistenceregions,asexpectedfora disorderedglobularstate. Toquantifythis,plotsofthenormalizedoverlapparameter,calculatedfromthesecontactmaps,areshowninFig.7. Aswith thecontactmaps,theseresultsexhibitmetastabilityaroundthetransition. Forallsequences,weseealargeincreaseoncooling in hQ i from values of hQ i≈0.3 in the disordered globule phase to values of hQ i≈0.85 in the low energy phase. This n n n confirmsthatthelowenergyphaseisdominatedbysmallfluctuationsaroundasingleconformation. V. DISCUSSIONANDCONCLUSIONS Ourresultsindicatethat,forallofthetenrandomheteropolymersequenceswhichwehavestudied,thereexistsalowtemper- aturetransitionbetweenthefamiliarcollapsedliquid-likeglobulephaseandalowenergyphase. Thesetransitionsdonothave theglassycharacterexpectedfromreplicatheory,showinginsteaddouble-peakedenergydistributionfunctionswithsubstantial associated latent heat. In each case, the low energy phase is more compact than the disordered globule, having a lower r , g exhibitslocaltranslationalorderandcrystalplanes,asshownbythetwo-particledensityfunctionr (2)(similartothelowenergy phasesseeninoff-latticehomopolymersimulations[11,12]),andisnotglassy,beingdominatedbyasingleconformationofthe polymerbackbone,asshownbythenormalizedthermalaverageoverlapparameterhQ i. Inanalogytothemorefamiliarphase n behaviorofbulksystems,wedescribethetransitionasafreezingtransition,andthelow-energystateasafrozenstate. Wenote theobviousparallelswithproteinfoldingtransitions, wherepolypeptideheteropolymersundergoa first-orderphasetransition toacollapsedsingleconformation,exhibitingrepeatingstructuralunits. Thelargecomputationaleffortassociatedwithsimulatinghigh-density,lowtemperaturecontinuumpolymerstateshaslimited thenumberofsequenceswe havebeenabletostudy. Nevertheless,wecanusetheNewcombe-Wilsonscoremethod(method 3 in Ref. [27]) to estimate with 95% confidence that at least 72.25% of 64 length 1:1 HP composition ratio sequences will demonstrateanequivalentfirst-orderfreezingtransitiontoanon-glassysingleconformationstate. Wesuggestthatthereasonthatweseedifferentbehaviorfromequivalentlatticemodels,whichundergolowtemperatureglass transitions,isamatterofsymmetries.Asevidencedbyther (2)data(seeFig.4),thedisorderedglobularstateforoursequences has no translational order, with hydrophobic monomers able to occupy any position, within the bounds of bond lengths and hard-sphereoverlap.Passingbelowthefreezingtemperature,thefrozenstatesofoursequencesgaincrystalline-likeorder[37]; hydrophobicmonomers must occupy lattice points, or points very near them. The spatial symmetry of the system has been brokeninamannerwhich,inthreedimensions,requiresafirst-orderphasetransition.Inthelatticemodel,theapproximationis thatspatialsymmetryisalreadybrokenforthesystematalltemperatures,andthereisnopossibilityofsuchanorder-disorder transition. Further,thesymmetryofthelatticeisfixedinadvance;ifthefrozenstate oftheoff-latticesystemisonadifferent lattice,orifsomedistortionofthelatticeisnecessaryforstability,thelatticemodelcannotaccuratelyportraytheorderedfrozen state. Withcooling,thelatticesystemwillremaininadisorderedstate,eventuallyformingaglass. Thecrystallinenatureofthefrozenstatesseenheredonotresembletheorderedstructuresseeninbiologicalheteropolymers (i.e., proteins); we see no signature of large scale helices or b -sheets in the contact maps. We believe this to be due to the completeflexibilityofthechainsinthemodel. Ithasbeensuggested[28]thathelicesandb -sheetlike“saddle”structuresare optimalshapes for a stiff stringlike object, in the same manner thatthe ubiquitousface centeredcubic crystal structure is the optimalarrangementforsetsofspheres.Wearecurrentlyinvestigatingthismatterforbothhomo-andheteropolymericsystems. Itcanbeseenfromthecontactmaps(seeFig.6)thatthefrozenconfigurationsarestronglysequencedependent.However,the thermodynamicbehavioroftheheteropolymers(Fig.2)hasverylittle(thoughstillstatisticallysignificant)sequencedependence. We believe that the thermodynamicbehavior, in the HP model at least, has stronger composition dependence than sequence dependence. Our1:1H:Pcompositionratio(withitsfreezingtemperaturesbunchedaroundT∗≈0.27)liesbetweenthelimits ofthepurehydrophobichomopolymer(withafreezingtemperatureofT∗=0.33[11]),andthepurepolarhomopolymerwith itseffectivelyhard-spheremonomers,whichshouldnotfreezeinisolation. Itisgenerallysuggestedthatthefractionofheteropolymericsequencescapableoffreezingtouniqueconformationsshould be small, around10−9 [29]. Dill and coworkersestimate the chanceof finding a proteinsequence which will fold to a given approximatestructureasbeingoforder10−10 [15,30,31]. Recentexperimentalwork,usingalargerandom-sequenceprotein librarytofindsequenceswhichbindtoATP,estimatesthefractionofproteinsequenceswhichare“functional”asaround10−11 [32]. For our system, we have estimated that over 70% of 64-length 1:1 composition ratio HP sequences freeze to unique conformations. Similar over-estimation is also seen for lattice models; the impressive direct enumeration performed by Li et al. [14] on an HP lattice model suggests that as many as 4.75% of all 27-length HP sequences have unique ground state configurations. Webelievesuchover-estimationtobeduetoalackofsourcesoffrustrationinourmodel,ascomparedtorealproteinsystems. Ourmodelhasonlytwomonomertypes,nosizepolydispersityorsidechains,completelyflexiblebonds,andnoorientational dependencein its interaction potential. Any of these factors could significantly lower the numberof model sequenceswhich exhibitfolding-likebehavior. In summary, we have demonstratedthat a simple continuumheteropolymermodelcan exhibit a first-order transition to an orderedlowtemperaturestatedominatedbyasingleconformation,withoutglassybehavior.Thisbehaviorissimilartothatseen inproteinfolding.Itwouldbeofinteresttoseewhetherextensionstothismodelyieldfurtherprotein-likecharacter,forexample secondarystructure,orinsteadintroducefrustrationswhichpreventtheobservedbehaviors. Acknowledgments ThisworkissupportedbytheBBSRC(grantreferenceB17005). [1] C.Anfinsen,Science181,223(1973). [2] V.Pande,A.Grosburg,andT.Tanaka,Rev.Mod.Phys.72,259(2000). [3] J.BryngelsonandP.Wolynes,ProceedingsoftheNationalAcademyofSciencesoftheUSA84,7524(1987). [4] B.Derrida,PhysicalReviewLetters45,79(1980). [5] E.ShakhnovichandA.Gutin,BiophysicalChemistry34,187(1989). [6] U.Bastolla,H.Frauenkron,andP.Grassberger,JournalofMolecularLiquids84,111(2000). [7] E.ShakhnovichandA.Gutin,JournalofChemicalPhysics93,5967(1990). [8] V.Pande,A.Grosberg,C.Joerg,andT.Tanaka,PhysicalReviewLetters76,3987(1996). [9] A.GrosbergandA.Khokhlov,StatisticalPhysicsofMacromolecules(Woodbury,1994). [10] P.Flory,Statisticalmechanicsofchainmolecules(Interscience,1969). [11] Y.Zhou,C.Hall,andM.Karplus,Phys.Rev.Lett.77,2822(1996). [12] Y.Zhou,M.Karplus,J.Wichert,andC.Hall,J.Chem.Phys.107,10691(1997). [13] V.Pande,A.Grosberg,andT.Tanaka,Biophys.J.73,3192(1997). [14] H.Li,R.Helling,C.Tang,andN.Wingreen,Science273,666(1996). [15] K.LauandK.Dill,ProceedingsoftheNationalAcademyofSciencesUSA87,638(1990). [16] H.Li,C.Tang,andN.Wingreen,Phys.Rev.Lett.79,765(1997). [17] J.Ilsley,M.Sudol,andA.Winder,CellularSignalling14,183(2002). [18] S.FukuchiandK.Nishikawa,JournalofMoleularBiology309,835(2001). [19] R.Sadus,MolecularSimulationofFluids(Elsevier,2002). [20] B.SmitandD.Frenkel,UnderstandingMolecularSimulation(AcademicPress,2002),2nded. [21] D.Kofke,JournalofChemicalPhysics117,6911(2002). [22] J.HansenandI.McDonald,TheoryofSimpleLiquids(AcademicPress,1986),2nded. [23] A.FerrenbergandR.Swendsen,Phys.Rev.Lett.63,1195(1989). [24] B.BergandT.Neuhaus,Phys.Rev.Lett.68,9(1992). [25] G.SmithandA.Bruce,J.Phys.A28,6623(1995). [26] A.Cooper,C.Johnson,J.Lakey,andM.Nöllmann,BiophysicalChemistry93,215(2001). [27] R.Newcombe,StatisticsinMedicine17,857(1998). [28] A.Maritan,C.Micheletti,A.Trovato,andJ.Banavar,Nature406,287(2000). [29] B.Alberts,A.Johnson,J.Lewis,M.Raff,K.Roberts,andP.Walter,Molecularbiologyofthecell(GarlandScience,2002),4thed. [30] H.ChanandK.Dill,JournalofChemicalPhysics95,3775(1991). [31] K.Dill,ProteinScience8,1166(1999). [32] A.KeefeandJ.Szostak,Nature410,715(2001). [33] A.Chertovich,V.Ivanov,B.Zavin,andA.Khoklov,Macromol.TheorySimul.11,751(2002). [34] L.Laaksonen,JournalofMolecularGraphics10,33(1992). [35] B.Bergman,L.Laaksonen,andA.Laaksonen,JournalofMolecularGraphicsandModelling15,301(1997). [36] Anunmodifiedhydrophobic-polarmatrixhaspolar-typemonomersinteractingashardspheres.ThisformfortheinteractionleadstoP- typerichsectionsoftheheteropolymersequenceforminglongloopsoffofacentralsphericalnucleusrichinhydrophobic-typemonomers [33];thisbehaviorisnotseeninrealproteins. [37] Wesaycrystalline-likeratherthancrystallinesinceitispossiblethat,asincaseoftwodimensionalcrystals,thistranslationalorderisnot long-ranged,andinthelargeN limitwouldbedestroyedbyfluctuations.However,thiswillbeunimportantforpolymerswhichfreeze intoobjectswithshortlengthscales,whichwenoteincludesrealproteinsaswellasthesimplemodelheteropolymersstudiedhere. Figures Figure1: Schematicdrawingofmolecularmodel;asystemofsphericalmonomers(filledcircles),diameters .Monomersarejoinedtogether inanon-branchinglinearsequencewithbondsoflengths (1−d )<l<s (1+d ).Non-bondedmonomersinteractiftheyarewithindistance ls ofeachother(shownwithdottedcircles);inthefigure,monomersi+1andi+2areinteracting. (a) -0.5 -1 -1.5 N *E/ -2 -2.5 -3 -3.5 20 (b) 15 *v10 c 5 0 (c) 0.064 0.062 0.25 0.2 *rg 0.06 0.15 0.1 0.058 0.05 0.056 0.15 0.5 1 1.5 0.15 0.5 1 1.5 * T Figure2: (a)Energyperunitmonomer,(b)heatcapacityand(c)radiusofgyrationasafunctionoftemperatureforthetensequencesstudied. Thescaleofplot(c)hasbeenchosentoshowthelowtemperaturebehavior;theinsetshowstherg behavioroverthefulltemperaturerange studied.CurvesarecalculatedatT∗intervalsof0.01usingW E,rg calculatedfromparalleltemperingrunsasdescribedinthetext.Estimated errorsaresmallerthansymbolsizeforradiusofgyrationandenergydata. (cid:0) (cid:1) (a) (b) 70 1.4 65 1.2 1 60 x e o N) 0.8 *Tc *E/ >/ 55 p( *E 0.6 D< 50 0.4 45 0.2 0 40 -3.2 -3 -2.8 -2.6 0.25 0.255 0.26 0.265 0.27 0.275 0.28 * * E /N T coex Figure3: (a)Coexistenceenergydistributions p(E∗/N)forthesequencesstudied. (b)CoexistencetemperatureT∗ plottedagainstenergy coex differenceD E/kBT betweenthecoexistingphases.Colordenotessequenceinbothgraphs.