Monte Carlo Simulation of Coil-to-Globule Transition of Compact Polymer Chains:Role of Monomer Interacting

Xin-ke Zhang,Jia-ye Su

Department of Applied Physics,Nanjing University of Science and Technology,Nanjing 210094,China

Coil-to-globule transitions are fundamental problems existing in polymer science for several decades;however,some features are still unclear,such as the eflect of chain monomer interaction.Herein,we use Monte Carlo simulation to study the coil-to-globule transition of simple compact polymer chains.We first consider the finite-size eflects for a given monomer interaction,where the short chain exhibits a one-step collapse while long chains demonstrate a two-step collapse,indicated by the specific heat.More interestingly,with the decrease of chain monomer interaction,the critical temperatures marked by the peaks of heat capacity shift to low values.A closer examination from the energy,mean-squared radius of gyration and shape factor also suggests the lower temperature of coil-to-globule transition.

Key words:Coil-to-globule transition,Polymer,Monte Carlo

I.INTRODUCTION

Half century ago Stockmayer predicted the coil-toglobule(CG)transition of polymer chains which has attracted extensive subsequent studies[1].In general,polymer chains will exhibit extended coils in good solvent or at high temperatures,and collapsed globules in poor solvent or at low temperatures.Thus,the change of solvent quality or temperature should be a direct way to drive the CG transition.In some early work,while the computer simulation was still luxury,scientists proposed pioneering theory for the CG transition,such as the work by Lifshitz[2]and Sanchez[3].Later,the theory for CG transition has been further explored extensively[4,5].Recently,Budkov and Kiselev presented a comprehensive review on the Flory-type theories of polymer collapse under various external stimuli[6],including the polymer collapse in miscible good solvents[7],the statistical theory of a polymer chain in a mixture of solvents[8,9],the collapse of strongly charged polyelectrolytes[10],as well as the polarizable polymer under external electric field[11?13].Remarkably,the CG transition is crucial to the physical properties of polymer materials and thus important for material design.Actually,the CG transition is also a key factor aflecting the function of proteins,where polypeptide chains can fold from coil to globule so as to demonstrate specific functions[14?16].Normally,these polypeptide chains contain both hydrophilic and hydrophobic residues,and thus the CG transition of heteropolymer chains has been considered to be the key role of protein folding[17?19].Simply,the principle conformations of the native states of globular proteins can be well described by compact polymer chain models[20].

Through Monte Carlo(MC)and molecular dynamics(MD)simulations,previous studies have extensively investigated the eflect of chain length,chain stiflness and confinement on the dynamic and thermodynamic information of CG transition[21?30].Meanwhile,due to local energy minima for long polymer chains,it is difficult to obtain their equilibrium conformations especially at low temperatures.Thus,some novel sampling methods,such as the Wang-Landau[25,26],integrated tempering sampling[31]and multicanonical ensemble methods[32]are developed to capture the CG transition.

It is believed that the folding and structures of proteins should highly depend on the interaction between polypeptide residues.Thus,it will be interesting and necessary for us to consider the CG transition of general model chains with varying monomer-monomer interactions.Herein,we use Monte Carlo simulations to investigate the CG transition of compact polymer chains,focusing on the role of chain monomer interaction.For a given monomer interaction,we study the finite-size eflects,and with the increase of chain length,the CG transition varies from one-step to two-step,indicated by the specific heat.We then consider the monomer interaction for a given chain length,and with the increase of interaction,the heat capacity shows a broad plateau instead of peaks.The mean-squared radius of gyration and shape factor suggest a smoother and less intense collapse.As the compact polymer chains are simple models of globular proteins,the results are helpful for us to understand the folding of proteins with diflerent residue interactions.

II.MODEL AND SIMULATION METHOD

In this work,we use lattice Monte Carlo simulations to study the CG transition of compact polymer chains.Actually,in an early work,similar model and method were used to study the folding of proteins[33].FIG.1 shows three typical chain structures for N=60 monomers from high to low temperatures,where the chain exhibits an extended coil,a globule,and a compact globule conformations,respectively.For the compact polymer chains,only nearest-neighbor interaction,e.g.,contact energy is considered,which is based on a self-avoiding walk.For simplicity,we do not consider the HP sequence[34],i.e.,all the chain monomers are hydrophobic and if two monomers are unbonded but with one lattice distance,they will have a contact energy of ε.We first consider the chain length of N=30,60,90,at ε=?1.Then,we study the chain length of N=90 at ε=?1,?2,?3.In the MC simulation,the system is gradually to reach the equilibrium states by the metropolis sampling algorithm:if Enew

In addition,similar to our previous work of a MD simulation[35],we use the thermal annealing technique for the chain to be relaxed sufficiently.In total,we consider 29 temperatures from T=10 to T=0.15,and we roughly divide them into four zones.For the first zone,we use t0=500000 MC steps for the chain to reach equilibrium and also another t0for data collection.As the chain moves slowly at low temperatures,for the other three zones,we use 2t0,3t0,4t0,respectively.We use our local computer code to run the simulation in a single CPU,and as there are so many temperatures,for the chain of N=90,the actual running time will be larger than two weeks.Thus,we do not consider larger chains further.In addition,we run ten independent simulations for each system and estimate the standard error bars.

In order to capture the thermodynamic features of the polymer chains,we calculate the contact energy and heat capacity.The heat capacity as a function of temperature is defined by the standard energy fluctuation[35],

where E is the contact energy.In general,the conformational changes of polymer chains will be reflected by peaks in the Cvcurve,as exemplified by many previous work[17,23,35?37].From thermodynamics point of view,isochoric thermal capacity can adequately describe the phase behavior of many-particle systems only in thermodynamic limit,i.e.in case of very big number of particles.However,for single polymer chains from several tens to hundreds of monomers,the heat capacity can well determine the coil-to-globule transition[17,23,35?37].To capture the polymer size during the CG transition,we calculate the radius of gyration,which is defined as the average distance of chain monomers from the center of mass.Clearly,it demonstrates how tightly packed a chain is.The“equivalent ellipsoid” of a confi guration can be obtained by evaluating the principal components L12≥L22≥L32of the squared radius of gyration S2=L12+L22+L32of individual configurations taken along the principal axes of inertia[39].Another normalized parameter,the shape factor ?δ??can well describe the shapes of polymer chains with regardless of chain lengths,

Where the shape factor changes between 0(sphere)and 1(rod).

III.RESULTS AND DISCUSSION

We first consider the CG transition of diflerent chain lengths with ε=?1.0,and the heat capacity and contact energy are shown as a function of temperature in FIG.2.From FIG.2(a)we can see clearly that there is a single peak for N=30,while there are two peaks with the same locations for N=60 and 90.Thus,the CG transition for short chains is a one-step collapse and for long chains is a two-step collapse.Actually,previous studies have well established the CG transition for square-well chains,and they suggest a three-step collapse,i.e.,gas-liquid,liquid-solid,and solid-solid ordering transitions[36,37].However,it is still difficult to observe the solid-solid transition,perhaps due to the ultralow temperatures.Thus,the extended coil structure in FIG.1(a)should be a gas-like state,while the globules in FIG.1(b)and(c)are corresponding to liquidlike and solid-like states,respectively.Generally,during the gas-liquid transition the average distance between monomers will have an intensive drop,and the energy should be changed similarly,as shown in FIG.2(b)from T=10.0 to T=1.0.While from T=1.0 to T=0.15,the energy change is relatively small.However,the liquidsolid transition is suggested to be a first-order transition[36,37],since the energy change is discontinuous,confi rmed by our contact energy in FIG.2(b).We also note that at low temperatures,for N=30 the contact number(energy)is close to the chain length;while it is clearly larger than the chain length for N=60 and 90.This indicates that the energy and conformation change of N=30 should be less intensive,which may be relevant to the single peak behavior of heat capacity.

FIG.1 Snapshot of the lattice polymer chain(N=60)at temperatures of(a)T=10.0,(b)T=1.0,(c)T=0.15 and contact energy of ε=?1.0.

FIG.2(a)The heat capacity and(b)average contact energy as a function of the temperature for compact polymer chains with length of N=30,60,90 and monomer interaction of ε=?1.0,respectively.Error bars are shown for ten data points.

During the CG transition,not only the energy but also the chain size has a remarkable change with the decrease of temperature,as shown in FIG.3.Notably,the decreasing behavior of mean-squared radius of gyration in FIG.3(a)should be another proof for the CG transition.We also note that the reduction of ?S2?mostly happens in the first Cvpeak region during the gas-liquid transition,in consistence with the change of contact energy;while the chain size only changes slightly for the liquid-solid transition.Thus,the chain should shrink to a globule first and then further tune some local position to a more compact structure.FIG.3(b)further shows the shape factor and at high temperatures,the values are slightly above 0.4,which is close to 0.39 of linear random walks[38].Due to the stiflness of polyethylene chains,our previous work shows a value of 0.69 for the linear chain[35].In fact,only for a sufficient sti ffchain,the value of shape factor can be close to 1.At low temperatures,the shape factor is below 0.1,implying a globule conformation.As a whole,the CG transition of present compact polymer chain is in good agreement with previous simulations in diflerent models[17,23?31,35?37]and experiments[39].

FIG.3(a)The mean-squared radius of gyration and(b)average shape factor as a function of the temperature for compact polymer chains with length of N=30,60,90 and monomer interaction of ε=?1.0,respectively.

FIG.4(a)The heat capacity and(b)average contact energy as a function of the temperature for compact polymer chains with length of N=60 and monomer interaction of ε=?0.1,?0.5,?1.0,respectively.

As we have well established the CG transition behaviors of compact polymer chains with a given monomer interaction,we then consider the eflect of monomer interaction.Similar to the CG transition of heteropolymer chains[17],the role of monomer interaction can help us to understand the folding of proteins with different residue interactions.FIG.4 shows the heat capacity and contact energy for N=60 with ε=?0.1,?0.5,?1.0.Remarkably,as seen in FIG.4(a),with the decrease of monomer interaction,the peak location shifts to lower temperatures.For ε=?0.5,?1.0,the two critical temperatures are at T1=0.75,T2=0.25,and T1=1.5,T2=0.5,respectively.Specifically,for ε=?0.1,there is no peak and the peak should exist in ultralow temperatures beyond our simulation range.When the chain monomer interaction is very low,it should be like a pure self-avoiding chain.Larger compact energy will lead to more compact conformation of chains and then the polymer collapse should happen at high temperatures,confirmed by the contact energy in FIG.4(b).For larger monomer interaction,the contact energy becomes to decrease at higher temperatures,corresponding to the heat capacity behavior.We also note that from the contact energy we can estimate the average contact number,and clearly at high or low temperatures the contact number of ε=?1.0 is about twice larger than that of ε=?0.5.Thus,the contact number should increase almost linearly as the monomer interaction increases.

FIG.5 (a)The mean-squared radius of gyration and(b)average shape factor as a function of the temperature for compact polymer chains with length of N=60 and monomer interaction of ε=?0.1,?0.5,?1.0,respectively.

As shown in FIG.5(a),for larger monomer interaction,the radius of gyration has smaller values at high temperatures,and also it decreases smoothly with the decrease of temperature rather than a sudden decrease for ε=?0.1.Obviously,the reduction region shifts to lower temperatures with the interaction decreasing and the ultimate values for ε=?0.1 are larger than those for ε=?0.5,?1.0.The change of shape factor in FIG.5(b)exhibits similar behaviors to the radius of gyration.Specifically,the final value at the lowest temperature for ε=?0.1 is about 0.2,which is about ten times larger than that for ε=?0.5, ?1.0.Thus,for ε=?0.1 the chain does not still collapses into a final globule structure.As a whole,the CG transition for chains with small monomer interaction will happen at lower temperatures,indicated by the shifting heat capacity.Considering the competition between temperature and monomer interaction,these results can be well understood.For the same temperature,large monomer interaction should lead to more compact chain conformation.Remarkably,for a given temperature if we change the monomer interaction step by step,the CG transition should also happen.

IV.CONCLUSION

We have systematically investigated the coil-toglobule transition of compact polymer chains by using the Monte Carlo simulations.We first consider the chain length eflect and find that the short chain exhibits a one-step collapse while long chains demonstrate a twostep collapse,indicated by the specific heat.We then focus on the eflect of chain monomer interaction.With the decrease of chain monomer interaction,the peaks of heat capacity shift to low temperatures and even disappear,confirmed by the similar change of contact energy,mean-squared radius of gyration and shape factor.These results can help us to understand the folding of proteins with diflerent residue interactions.

V.ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China(No.21574066 and No.21204093).The allocated computer time at the National Supercomputing Center in Shenzhen is gratefully acknowledged.

[1]W.H.Stockmayer,Makromol.Chem.35,54(1960).

[2]I.M.Lifshitz,A.Yu.Grosberg,and A.R.Khokhlov,Rev.Mod.Phys.50,683(1978).

[3]I.C.Sanchez,Macromolecules 12,980(1979).

[4]D.S.Simmons and I.C.Sanchez,Macromolecules 41,5885(2008).

[5]D.S.Simmons and I.C.Sanchez,Macromolecules 46,4691(2013).

[6]Yu.A.Budkov and M.G.Kiselev,J.Phys.:Condens.Matter 30,043001(2018).

[7]D.Mukherji,C.M.Marques,and K.Kremer,Nat.Commun.5,4882(2014).

[8]Yu.A.Budkov,A.L.Kolesnikov,N.N.Kalikin,and M.G.Kiselev,EPL 114,46004(2016).

[9]Yu.A.Budkov and A.L.Kolesnikov,Soft Matter 13,8362(2017).

[10]A.M.Tom,S.Vemparala,R.Rajesh,and Nikolai V.Brilliantov,Phys.Rev.Lett.117,147801(2016).

[11]Yu.A.Budkov,A.L.Kolesnikov,and M.G.Kiselev,J.Chem.Phys.143,201102(2015).

[12]A.L.Kolesnikov,Yu.A.Budkov,E.A.Basharovac,and M.G.Kiselev,Soft Matter 13,4363(2017).

[13]Y.A.Budkov and A.L.Kolesnikov,Eur.Phys.J.E 39,110(2016).

[14]A.Gutin,V.Abkevich,and E.Shakhnovich,Phys.Rev.Lett.80,208(1998).

[15]V.S.Pande,A.Y.Grosberg,and T.Tanaka,Rev.Mod.Phys.72,259(2000).

[16]E.Sherman and G.Haran,Proc.Natl.Acad.Sci.USA 103,11539(2006).

[17]W.Wang,P.Zhao,X.Yang,and Z.Y.Lu,J.Chem.Phys.141,244907(2014).

[18]D.Yang and Q.Wang,ACS Macro Lett.2,952(2013).

[19]R.Wang and Z.G.Wang,Macromolecules 47,4094(2014).

[20]H.S.Chan and K.A.Dill,Macromolecules 22,4559(1989).

[21]A.M.Rubio,J.J.Freire,J.H.R.Clarke,C.W.Yong,and M.Bishop,J.Chem.Phys.102,2277(1995).

[22]G.Tanaka and W.L.Mattice,Macromolecules 28,1049(1995).

[23]H.Liang and H.Chen,J.Chem.Phys.113,4469(2000).

[24]H.Noguchi and K.Yoshikawa,J.Chem.Phys.109,5070(1998).

[25]D.T.Seaton,T.Wust,and D.P.Landau,Phys.Rev.E 81,011802(2010).

[26]Z.Wang,L.Wang,and X.He,Soft Matter 9,3106(2013).

[27]M.Bachmann and W.Janke,J.Chem.Phys.120,6779(2004).

[28]Z.T.Jiang,W.H.Dou,Y.Shen,T.T.Sun,and P.Xu,Chin.Phys.B 24,116802(2015).

[29]Z.T.Jiang,W.H.Dou,T.T.Sun,Y.Shen,and D.Cao,J.Polym.Res.22,236(2015).

[30]Z.Y.Yang,A.H.Chai,P.Li,and Y.F.Yang,Chin.J.Chem.Phys.29,381(2016).

[31]P.Zhao,L.J.Yang,Y.Q.Gao,and Z.Y.Lu,Chem.Phys.415,98(2013).

[32]S.Schnabel,M.Bachmann,and W.Janke,J.Chem.Phys.131,124904(2009).

[33]E.Shakhnovich,G.Farztdinov,A.M.Gutin,and M.Karplusm,Phys.Rev.Lett.67,1665(1991).

[34]L.Zhang and J.Su,Polymer 47,735(2006).

[35]J.Su,L.Zhang,and H.Liang,J.Chem.Phys.129,044905(2008).

[36]Y.Zhou,C.K.Hall,and M.Karplus,Phys.Rev.Lett.77,2822(1996).

[37]Y.Zhou,M.Karplus,J.M.Wichert,and C.K.Hall,J.Chem.Phys.107,10691(1997).

[38]G.Ziflerer,J.Chem.Phys.102,3720(1995).

[39]J.Xu,Z.Y.Zhu,S.Z.Luo,C.Wu,and S.Y.Liu,Phys.Rev.Lett.96,027802(2006).