Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (6): 784-788

The article information

Xin-ke Zhang, Jia-ye Su
张鑫柯, 苏加叶
Monte Carlo Simulation of Coil-to-Globule Transition of Compact Polymer Chains: Role of Monomer Interacting
紧密高分子链从无规线团到球形相变的蒙特卡洛模拟:链单体的相互作用
Chinese Journal of Chemical Physics, 2018, 31(6): 784-788
化学物理学报, 2018, 31(6): 784-788
http://dx.doi.org/10.1063/1674-0068/31/cjcp1801002

Article history

Received on: January 16, 2018
Accepted on: April 22, 2018
Monte Carlo Simulation of Coil-to-Globule Transition of Compact Polymer Chains: Role of Monomer Interacting
Xin-ke Zhang, Jia-ye Su     
Dated: Received on January 16, 2018; Accepted on April 22, 2018
Department of Applied Physics, Nanjing University of Science and Technology, Nanjing 210094, China
*Author to whom correspondence should be addressed. Jia-ye Su, E-mail: jysu@iccas.ac.cn, jysu@njust.edu.cn
Abstract: Coil-to-globule transitions are fundamental problems existing in polymer science for several decades; however, some features are still unclear, such as the effect 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 effects 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    
Ⅰ. INTRODUCTION

Half century ago Stockmayer predicted the coil-to-globule (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 affecting 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 effect of chain length, chain stiffness 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 effects, 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 different residue interactions.

Ⅱ. 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 $\varepsilon$. We first consider the chain length of $N$=30, 60, 90, at $\varepsilon$=-1. Then, we study the chain length of $N$=90 at $\varepsilon$=-1, -2, -3. In the MC simulation, the system is gradually to reach the equilibrium states by the metropolis sampling algorithm: if $E_{\textrm{new}}$$ <$$E_{\textrm{old}}$, the new conformation of the chain will be accepted, and if $E_{\textrm{new}}$$\geq$$E_{\textrm{old}}$, the new conformation will be accepted with the probability of exp[-($E_{\textrm{new}}$-$E_{\textrm{old}}$)/$k_\textrm{B}T$]. In this way, the polymer chain will reach low energy states for a given temperature.

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 $\varepsilon$=-1.0.

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 $t_0$=500000 MC steps for the chain to reach equilibrium and also another $t_0$ for data collection. As the chain moves slowly at low temperatures, for the other three zones, we use 2$t_0$, 3$t_0$, 4$t_0$, 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],

$ \begin{eqnarray} C_\textrm{v}=\frac{\langle E^2\rangle-\langle E\rangle^2}{k_\textrm{B} T^2} \end{eqnarray} $ (1)

where $E$ is the contact energy. In general, the conformational changes of polymer chains will be reflected by peaks in the $C_\textrm{v}$ curve, 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 configuration can be obtained by evaluating the principal components $L_1$$^2$$\geq$$ L_2$$^2$$\geq$$L_3$$^2$ of the squared radius of gyration $S^2$=$L_1$$^2$+$L_2$$^2$+$L_3$$^2$ of individual configurations taken along the principal axes of inertia [39]. Another normalized parameter, the shape factor $\langle \delta^*\rangle$ can well describe the shapes of polymer chains with regardless of chain lengths,

$ \begin{eqnarray} \langle \delta^*\rangle=1-3\left\langle\frac{L_1{^2} L_2{^2}+L_2{^2} L_3{^2}+L_{1^2} L_3{^2}}{(L_1{^2}+L_2{^2}+L_3{^2})^2} \right\rangle \end{eqnarray} $ (2)

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

Ⅲ. RESULTS AND DISCUSSION

We first consider the CG transition of different chain lengths with $\varepsilon$=-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 liquid-like 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 liquid-solid transition is suggested to be a first-order transition [36, 37], since the energy change is discontinuous, confirmed 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. 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 $\varepsilon$=-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 $\langle$$S^2$$\rangle$ mostly happens in the first $C_\textrm{v}$ peak 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 stiffness of polyethylene chains, our previous work shows a value of 0.69 for the linear chain [35]. In fact, only for a sufficient stiff chain, 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 different 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 $\varepsilon$=-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 effect 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 $\varepsilon$=-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 $\varepsilon$=-0.5, -1.0, the two critical temperatures are at $T_1$=0.75, $T_2$=0.25, and $T_1$=1.5, $T_2$=0.5, respectively. Specifically, for $\varepsilon$=-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 $\varepsilon$=-1.0 is about twice larger than that of $\varepsilon$=-0.5. Thus, the contact number should increase almost linearly as the monomer interaction increases.

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 $\varepsilon$=-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 $\varepsilon$=-0.1. Obviously, the reduction region shifts to lower temperatures with the interaction decreasing and the ultimate values for $\varepsilon$=-0.1 are larger than those for $\varepsilon$=-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 $\varepsilon$=-0.1 is about 0.2, which is about ten times larger than that for $\varepsilon$=-0.5, -1.0. Thus, for $\varepsilon$=-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.

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 $\varepsilon$=-0.1, -0.5, -1.0, respectively.
Ⅳ. CONCLUSION

We have systematically investigated the coil-to-globule transition of compact polymer chains by using the Monte Carlo simulations. We first consider the chain length effect and find that the short chain exhibits a one-step collapse while long chains demonstrate a two-step collapse, indicated by the specific heat. We then focus on the effect 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 different residue interactions.

Ⅴ. 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.

Reference
[1] W. H. Stockmayer, Makromol. Chem. 35 , 54 (1960). DOI:10.1002/macp.1960.020350103
[2] I. M. Lifshitz, A. Yu. Grosberg, and A. R. Khokhlov, Rev. Mod. Phys. 50 , 683 (1978). DOI:10.1103/RevModPhys.50.683
[3] I. C. Sanchez, Macromolecules 12 , 980 (1979). DOI:10.1021/ma60071a040
[4] D. S. Simmons, and I. C. Sanchez, Macromolecules 41 , 5885 (2008). DOI:10.1021/ma800151p
[5] D. S. Simmons, and I. C. Sanchez, Macromolecules 46 , 4691 (2013). DOI:10.1021/ma400338d
[6] Yu. A. Budkov, and M. G. Kiselev, J. Phys.: Condens. Matter 30 , 043001 (2018). DOI:10.1088/1361-648X/aa9f56
[7] D. Mukherji, C. M. Marques, and K. Kremer, Nat. Commun. 5 , 4882 (2014). DOI:10.1038/ncomms5882
[8] Yu. A. Budkov, A. L. Kolesnikov, N. N. Kalikin, and M. G. Kiselev, EPL 114 , 46004 (2016). DOI:10.1209/0295-5075/114/46004
[9] Yu. A. Budkov, and A. L. Kolesnikov, Soft Matter 13 , 8362 (2017). DOI:10.1039/C7SM01637A
[10] A. M. Tom, S. Vemparala, R. Rajesh, and V. Nikolai Brilliantov, Phys. Rev. Lett. 117 , 147801 (2016). DOI:10.1103/PhysRevLett.117.147801
[11] Yu. A. Budkov, A. L. Kolesnikov, and M. G. Kiselev, J. Chem. Phys. 143 , 201102 (2015). DOI:10.1063/1.4936661
[12] A. L. Kolesnikov, Yu. A. Budkov, E. A. Basharovac, and M. G. Kiselev, Soft Matter 13 , 4363 (2017). DOI:10.1039/C7SM00417F
[13] Y. A. Budkov, and A. L. Kolesnikov, Eur. Phys. J. E 39 , 110 (2016). DOI:10.1140/epje/i2016-16110-x
[14] A. Gutin, V. Abkevich, and E. Shakhnovich, Phys. Rev. Lett. 80 , 208 (1998). DOI:10.1103/PhysRevLett.80.208
[15] V. S. Pande, A. Y. Grosberg, and T. Tanaka, Rev. Mod. Phys. 72 , 259 (2000). DOI:10.1103/RevModPhys.72.259
[16] E. Sherman, and G. Haran, Proc. Natl. Acad. Sci. USA 103 , 11539 (2006). DOI:10.1073/pnas.0601395103
[17] W. Wang, P. Zhao, X. Yang, and Z. Y. Lu, J. Chem. Phys. 141 , 244907 (2014). DOI:10.1063/1.4904888
[18] D. Yang, and Q. Wang, ACS Macro Lett. 2 , 952 (2013). DOI:10.1021/mz400457h
[19] R. Wang, and Z. G. Wang, Macromolecules 47 , 4094 (2014). DOI:10.1021/ma5003968
[20] H. S. Chan, and K. A. Dill, Macromolecules 22 , 4559 (1989). DOI:10.1021/ma00202a031
[21] A. M. Rubio, J. J. Freire, J. H. R. Clarke, C. W. Yong, and M. Bishop, J. Chem. Phys. 102 , 2277 (1995). DOI:10.1063/1.468748
[22] G. Tanaka, and W. L. Mattice, Macromolecules 28 , 1049 (1995). DOI:10.1021/ma00108a036
[23] H. Liang, and H. Chen, J. Chem. Phys. 113 , 4469 (2000). DOI:10.1063/1.1288186
[24] H. Noguchi, and K. Yoshikawa, J. Chem. Phys. 109 , 5070 (1998). DOI:10.1063/1.477121
[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). DOI:10.1039/c3sm27454c
[27] M. Bachmann, and W. Janke, J. Chem. Phys. 120 , 6779 (2004). DOI:10.1063/1.1651055
[28] Z. T. Jiang, W. H. Dou, Y. Shen, T. T. Sun, and P. Xu, Chin. Phys. B 24 , 116802 (2015). DOI:10.1088/1674-1056/24/11/116802
[29] Z. T. Jiang, W. H. Dou, T. T. Sun, Y. Shen, and D. Cao, J. Polym. Res. 22 , 236 (2015). DOI:10.1007/s10965-015-0875-3
[30] Z. Y. Yang, A. H. Chai, P. Li, and Y. F. Yang, Chin. J. Chem. Phys. 29 , 381 (2016). DOI:10.1063/1674-0068/29/cjcp1511231
[31] P. Zhao, L. J. Yang, Y. Q. Gao, and Z. Y. Lu, Chem. Phys. 415 , 98 (2013). DOI:10.1016/j.chemphys.2012.12.043
[32] S. Schnabel, M. Bachmann, and W. Janke, J. Chem. Phys. 131 , 124904 (2009). DOI:10.1063/1.3223720
[33] E. Shakhnovich, G. Farztdinov, A. M. Gutin, and M. Karplusm, Phys. Rev. Lett. 67 , 1665 (1991). DOI:10.1103/PhysRevLett.67.1665
[34] L. Zhang, and J. Su, Polymer 47 , 735 (2006). DOI:10.1016/j.polymer.2005.11.081
[35] J. Su, L. Zhang, and H. Liang, J. Chem. Phys. 129 , 044905 (2008). DOI:10.1063/1.2957486
[36] Y. Zhou, C. K. Hall, and M. Karplus, Phys. Rev. Lett. 77 , 2822 (1996). DOI:10.1103/PhysRevLett.77.2822
[37] Y. Zhou, M. Karplus, J. M. Wichert, and C. K. Hall, J. Chem. Phys. 107 , 10691 (1997). DOI:10.1063/1.474186
[38] G. Zifferer, J. Chem. Phys. 102 , 3720 (1995). DOI:10.1063/1.468554
[39] J. Xu, Z. Y. Zhu, S. Z. Luo, C. Wu, and S. Y. Liu, Phys. Rev. Lett. 96 , 027802 (2006). DOI:10.1103/PhysRevLett.96.027802
紧密高分子链从无规线团到球形相变的蒙特卡洛模拟:链单体的相互作用
张鑫柯, 苏加叶     
南京理工大学应用物理系,南京 210094
摘要: 本文采用蒙特卡洛模拟研究紧密高分子单链的无规线团到球形转变过程.首先,研究某一固定链单体的相互作用下的有效尺寸效应,并发现短链表现出一次的坍塌相变,而长链则出现两次,这从热熔曲线上可以明显看出.随着链单体相互作用的减小,热熔曲线上峰值对应的转变温度朝着低温转移.从体系的能量、均方回转半径和形状因子可以进一步看出,高分子相变发生在更低的温度区间.这些结果有助于加深理解高分子坍塌相变过程中链单体相互作用的影响.
关键词: 无规线团到球形的转变    高分子    蒙特卡洛