Coil-to-globule transitions of homopolymers and multiblock copolymers Wei Wang, Peng Zhao, Xi Yang, and Zhong-Yuan Lu Citation: The Journal of Chemical Physics 141, 244907 (2014); doi: 10.1063/1.4904888 View online: http://dx.doi.org/10.1063/1.4904888 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/24?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Collapse transitions in thermosensitive multi-block copolymers: A Monte Carlo study J. Chem. Phys. 140, 204904 (2014); 10.1063/1.4875694 Coil-to-globule transition by dissipative particle dynamics simulation J. Chem. Phys. 134, 244904 (2011); 10.1063/1.3604812 Globular structures of a helix-coil copolymer: Self-consistent treatment J. Chem. Phys. 126, 034902 (2007); 10.1063/1.2403868 A lattice model Monte Carlo study of coil-to-globule and other conformational transitions of polymer, amphiphile, and solvent J. Chem. Phys. 112, 7711 (2000); 10.1063/1.481363 Observation of the molten globule state in a Monte Carlo simulation of the coil-to-globule transition of a homopolymer chain J. Chem. Phys. 110, 10212 (1999); 10.1063/1.478893
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
THE JOURNAL OF CHEMICAL PHYSICS 141, 244907 (2014)
Coil-to-globule transitions of homopolymers and multiblock copolymers Wei Wang, Peng Zhao, Xi Yang, and Zhong-Yuan Lua) Institute of Theoretical Chemistry, State Key Laboratory of Theoretical and Computational Chemistry, Jilin University, Changchun 130023, China
(Received 1 October 2014; accepted 9 December 2014; published online 30 December 2014) We study the coil-to-globule transitions of both homopolymers and multiblock copolymers using integrated tempering sampling method, which is a newly proposed enhanced sampling method that can efficiently sample the energy space with low computational costs. For homopolymers, the coil-to-globule structure transition temperatures (Ttr) are identified by the radius of gyration of the chain. The transition temperature shows a primary scaling dependence on the chain length (N) with Ttr ∼ N −1/2. For multiblock copolymers, the coil-to-globule transition can be identified as first order, depending on the block size and the difference in attractive interactions of blocks. The influence of mutating a small portion of strongly attractive blocks to weakly attractive blocks on the coil-to-globule transition is found to be related to the position of the mutation. C 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4904888]
I. INTRODUCTION
The change of a single polymer chain from an extended coil in good solvent to a collapsed globule in poor solvent is called coil-to-globule (C2G) transition.1 This collapse typically corresponds to a second-order phase transition for homopolymers with suitable lengths.2,3 The importance of the C2G transition is due to its inherent connection with the problem of tuning function of a protein.4,5 The key factor influencing protein function strongly relates to a single polypeptide chain folding from coil state to globule state with specific three-dimensional structures. In this process, the protein, in which some amino acid residues are hydrophilic or charged while others are hydrophobic,4,5 undergoes a firstorder folding transition to reach the native state.4–6 Thus, heteropolymers or even block copolymer models (resembling protein structures) are widely used to study the protein folding for a long time, and the transition process together with C2G transition of both homopolymer and heteropolymer chain had been explained theoretically.4,7–11 Computer simulations can provide more structural, thermodynamic, and dynamic information for C2G transition of a polymer chain in dilute solution. Both Monte Carlo (MC) and molecular dynamics (MD) methods with either implicit or explicit solvent models have been used to elucidate the influence of chain length, chain stiffness, confinement, etc., on the C2G transition temperature, order of transition, and even transition dynamics.12–26 But for long chains, the energy landscape of the system is always too rugged with many local minima and barriers, thus, the ergodicity may not be guaranteed in normal simulations. To efficiently sample the configuration space, enhanced sampling methods, such as Wang-Landau methods24,25 and multicanonical ensemble methods,26–28 are usually adopted in the studies of C2G a)Electronic address: [email protected]
0021-9606/2014/141(24)/244907/7/$30.00
transitions of polymer chains. With the aid of enhanced sampling methods, C2G transitions of homopolymers and several sequences of heteropolymers were studied for giving hints of the formation of protein native structures.27,29 The multiblock copolymers, which can be tuned to homopolymers if all the blocks have the same property or to polypeptides if the block hydrophobicities are designed to mimic specific protein sequence, can serve as a versatile system bridging homopolymer model and heteropolymer model in the studies of C2G transition. The influence of inherent driving force for protein folding can be well illustrated in multiblock copolymer models by designing block sizes, positions, and interactions in a chain. The most prominent examples are the lattice HP model30–32 and the off-lattice AB model,33,34 which had been exhaustively investigated.24,27,35–37 Moreover, previously overlooked influence of mutation38 in strongly attractive blocks can be systematically studied with a regularly distributed multiblock copolymer model. Therefore, focusing on these issues, we study the C2G transitions of regularly distributed block copolymers and compare the results to those of homopolymers with a novel enhanced sampling method–integrated tempering sampling (ITS)39,40 implemented in GPU MD code GALAMOST41 in this work. This paper is organized as follows: in Sec. II, the models of homopolymer and block copolymers are described. The ITS method will be briefly introduced in Sec. III, in which the main idea and the most important formulae are given. The MD simulation details are also illustrated in this section. The results and discussion are shown in Sec. IV, followed by Sec. V for conclusions.
II. MODELS
Common bead-spring models are used in this study to describe homopolymers and block copolymers. The homopolymer is a linear chain with one type of beads
141, 244907-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-2
Wang et al.
J. Chem. Phys. 141, 244907 (2014)
connected by bonds with finite extensible nonlinear elastic (FENE) potential between connected beads
representing the contribution from each temperature Tk . The biased potential U ′ is defined through
1 r2 − Kr 02 ln *1.0 − 2 + +UWCA r ≤ r 0 UFENE(r) = 2 , (1) r0 , ∞ r > r0 where ( 1 σWCA ) 12 ( σWCA ) 6 4ϵ WCA − r ≤ 2 6 σWCA r r UWCA(r) = , 1 0 6 σWCA r > 2 (2)
e−βU (r ) ≡ W (r),
(5)
1 U ′(r) = − ln nk e−β kU (r ). β k
(6)
and r 0 is the bond extension parameter, K is the force constant, and r is the instantaneous bond length. To be consistent with previous work,40 we set σWCA = 1.05, ϵ WCA = 1.0, r 0 = 1.5, and K = 40. The Lennard-Jones potential is used between non-bonded beads with σ = 1.0 and ϵ = 1.0 (which, together with bead mass m ≡ 1.0, are taken as the units in this study), ( ) 12 ( ) 6 σ σ − r ≤ r cut, 4ϵ r r , (3) ULJ (r) = 0 r > r cut with r cut being set to 3.0 in our simulations. We consider homopolymers with chain lengths (N) between 13 and 800. The solvent effects are treated by implicit solvent model. For block copolymers, we use two types of beads denoted by A and B, respectively. The linear block copolymer chain with blocks of A and B type beads are connected by FENE potential. To represent the incompatibility between A and B beads, we use Lennard-Jones potential but with different energy parameter ϵ to characterize the strength of short-range attractive interaction. To describe the amphiphilicity of the block copolymer chain, we set ϵ AA = 1.0 and ϵ BB = 0.2 in most of the cases. The interaction parameter between A and B beads is set the same as that between B type beads, i.e., ϵ AB = ϵ BB = 0.2, which is consistent with the “AB” model used in the literature.27 To represent the effect of changing block copolymer amphiphilicity, we keep the value of ϵ AB and ϵ BB the same and adjust the value in the range of 0.2–0.8.
III. SIMULATION METHOD
The integrated tempering sampling method is a novel easy-to-use enhanced sampling method based on generalized ensemble to generate a distribution covering a broad range of energies.39,40 Here, we only give a brief introduction to this method, since details of the calculation process have been introduced in previous work.40 The main idea in ITS method is to use a biased potential generated from the combination of a series of canonical energy distributions. This allows efficient sampling in a desired temperature and energy range without requiring a predefined reaction coordinate. The generalized distribution W (r) is defined as W (r) = nk e−β kU (r ), (4) k
where U is potential energy, βk = 1/k BT (k B is Boltzmann constant and T is temperature), and nk is the weighting factor
′
thus,
The biased force Fb is simply the force in conventional MD simulations, F, carrying a pre-factor ∂U ′ ∂U ′ ∂U nk βk e−β kU (r ) F, (7) Fb = − =− = k ∂r ∂U ∂r β k nk e−β kU (r ) thus, ITS method can be readily implemented in any MD code by slightly modifying the force subroutine. The thermodynamic properties of any canonical ensemble whose temperature (T j ) is in the desired range can be calculated from the generalized ensemble by using reweighting A(r )e −β j U (r ) W (r)dr A(r)e−β jU (r )dr W (r ) ⟨A⟩ β j = −β U (r ) = −β U (r ) j e e j dr W (r ) W (r)dr −β j U (r ) =
A(r )e W (r )
e
−β j U (r )
W (r )
W
.
(8)
W
Typical computational procedure of ITS simulation is therefore as follows: 1. Determine the desired temperature range. 2. Choose a set of temperatures in the desired range and run several short time conventional MD simulations. Then, the relation between potential energy and temperature can be obtained by interpolation. Here, we define an overlap factor t, which can actually be determined by the requirement that the ratio between energy probability density functions at two adjacent temperatures is a constant. Giving a proper value of the overlap factor t, sufficient overlap of the energy distribution between adjacent temperatures can be guaranteed. 3. Determine the ITS temperature distribution and the corresponding weighting factors nk as described in Ref. 40. 4. Use the parameters generated in steps 2 and 3 to perform ITS simulation, which is essentially nothing but a conventional MD simulation using biased force calculated by Eq. (7). 5. After ITS simulation, the canonical ensemble properties can be calculated by Eq. (8). Nosè-Hoover thermostat is adopted in our ITS simulations. The potential energies at different temperatures are first obtained by 1.0 × 106 steps MD simulations at 8 temperatures in the range of 1.0–5.0 for homopolymers and in the range of 0.5–5.0 for block copolymers. The dependence of potential energy on temperature is used to estimate nk .40 The overlap factor t is set to e0.5, thus, 8–83 temperatures are generated in the range of 1.0–4.0 for systems with different chain lengths.40 The specific numbers of temperatures are shown in Table I. For homopolymers with different chain lengths, we then perform 1.0 × 109 steps ITS simulations with GALAMOST package41 at each chain length, and the data are recorded every 103 steps. Using a single ITS simulation
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-3
Wang et al.
J. Chem. Phys. 141, 244907 (2014)
TABLE I. The number of temperatures in ITS simulations. Chain length
Number of temperatures
13 25 50 100 200 400 800
8 12 18 27 40 58 83
trajectory, we can obtain thermodynamic properties of the system at any temperature in the temperature range. For multiblock copolymers, we exclusively focus on the influence of block size, strength of attraction, and mutation on the C2G transitions, so only chain with length N = 100 is considered, and the numbers of temperatures to perform ITS simulations are 25.
IV. RESULTS AND DISCUSSION A. Homopolymers
We calculate the radius of gyration (Rg ) of the chain and its derivative and the specific heat (Cv ) as structural and thermodynamic indicators, respectively, for C2G transitions for homopolymer with N = 13 ∼ 800. The radius of gyration (Rg ) is defined as
Rg2
N N 1 2 = r , 2N 2 i=1 j=1 i j
and the specific heat (Cv ) is defined as
2 E − ⟨E⟩2 ∂U . = Cv = ∂T k BT 2
(9)
(10)
The dependence of both Rg and Cv on temperature for homopolymers with different lengths is shown in Fig. 1. In Fig. 1(a), we can see the peaks in the derivative of Rg for all the chains, which clearly indicate strong structure changes corresponding to C2G transition. In Fig. 1(b), we can see for N < 200, there is no clear signal in the Cv curves for C2G transition, which implies strong finite size effect. While for N > 200, Cv curves have shoulders at the temperatures related to the peaks in the curves of the first derivative of Rg in Fig. 1(a), which can be recognized as the energetic signal of the C2G transition, similar to those reported in Ref. 24. No sharp peak in the Cv curves implies that the C2G transition for homopolymers is not thermodynamically first-order, which is consistent with previous conclusions.6 We can also clearly see the transition temperatures (Ttr) identified from Fig. 1(a) move to higher value with increasing chain length. To further elucidate the dependence of Ttr on the chain length, we fit Ttr versus N with the following equation: √ (11) Ttr(N) −TΘ = a1/ N + a2/N, where Ttr(N) is the C2G transition temperature for each chain length and TΘ is the Θ temperature characterizing the turning
FIG. 1. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the specific heat C v per monomer versus temperature.
point of structure transition of a polymer chain.3 Here, the √ N term corresponds to a mean field fitting, and the N term is adopted for finite chain length correction. As shown in Fig. 2, the fitting is quite acceptable with parameters a1 = −12.43, a2 = 9.58, and TΘ = 3.696, which are in agreement with those reported by using Wang-Landau method.24 Combining the results shown above for homopolymers, we can verify that ITS is an effective method for elucidating C2G transitions. Fig. 1(b) also shows rugged regions at low temperatures especially when the chain length is larger than 200. It may be
FIG. 2. The dependence of the transition temperature on homopolymer chain length with N from 13 to 800. The dashed line shows fitting with Eq. (11).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-4
Wang et al.
attributed to the existence of metastable intermediate states for long chains that have been identified in experiments42 and in simulations,43 but it also may be attributed to the inefficiency of sampling with ITS method for long chains. However, since we are focusing C2G transition (which takes place at comparatively higher temperature) only, we will not pay much attention to these rugged regions at low temperatures in the following studies. B. Heteropolymers
In this section, we exclusively study the influence of block size, amphiphilicity between blocks, and mutation at specific sites on the C2G transitions of multiblock copolymers. Since changing chain length only results in shift of C2G transition temperature, in the following, we will keep the chain length constant at N = 100 and consider four types of multiblock copolymers with different block sizes: (A2 B2)25, (A5 B5)10, (A10 B10)5, and (A25 B25)2. We keep ϵ AA = 1.0 to manifest A blocks strongly attractive and ϵ BB = ϵ AB = 0.2 to manifest B blocks weakly attractive and incompatible with A blocks. We calculate the first derivative of radius of gyration and Cv at different temperatures, as shown in Figs. 3(a) and 3(c), respectively. To further clarify the structure transition characteristics of multiblock copolymer, we also calculate the shape factor ρ, defined as ρ = Rg /Rh , where Rg is the radius of gyration, and Rh is the strongly attractive radius N N
−1 1 r−1 . (12) Rh = 2 N i=1 j=1(i, j) i j The shape factor ρ was widely used in light scattering techniques to probe the shape, size, or structure of certain molecule systems.44,47 Thus, it is used here to measure the conformational change when C2G transitions of single polymer chains occur. Generally, the value of ρ is about 1.5 for flexible chains in good solvent, while ρ = 0.77 for a solid sphere.44,47 In Fig. 3(b), we can see that the shape factor
J. Chem. Phys. 141, 244907 (2014)
is around 1.2 for chains at high temperature (corresponding to good solvent condition), which is a little lower than the theoretical value for coil state in good solvent. It is mainly attributed to the strongly attractive parts of the multiblock copolymer chain, which make the conformation of the chain not extended as much as an ideal coil. However, the C2G transitions can be clearly discerned from the turning points in the curves. To further justify the order of the transition other than using Cv , we also calculate the Binder cumulant VE 45 for all the chains at different temperatures. The Binder cumulant has an isolated minimum VE < 2/3 at the transition point,46 thus can be a good energy indicator. VE is defined as
4 E (13) VE = 1 − 2 . 3 E2 We can see from Figs. 3(a) and 3(b) that both the first derivative of mean square of radius of gyration and the shape factor give clear signals of the C2G transition, and the transition temperature increases with increasing block size. As compared to homopolymers, Fig. 3(c) shows a clear difference: obvious peaks from the Cv curves of multiblock copolymers can be observed, while there are only shoulders in the cases of homopolymers (Fig. 1(b)). This implies that C2G transition is first order for certain multiblock copolymers studied here. The Binder cumulant in Fig. 3(d) also gives consistent results, showing valleys at the same positions as compared to the peaks in Cv curves. From Fig. 3, we can also find that, the transition temperatures indicated by the peaks of the structural indicators are a little bit higher than those from energetic indicators. It may be due to the fact that the response of energetic properties is not sensitive as compared to the structural properties. The peaks in the structural curves correspond to the maximum changes of the molecular structures, while the energy changes may show apparent signal when the polymer chains are compact enough (at lower temperatures). It should be noted that the shape factor of (A25 B25)2 has two turning points, as shown in Fig. 3(b). From Fig. 3(a),
FIG. 3. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the shape factor ρ versus temperature; (c) the specific heat C v versus temperature; (d) the Binder cumulant versus temperature for multiblock copolymer chains with different block sizes. In all cases, ϵ AA = 1.0 and ϵ AB = ϵ BB = 0.2.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-5
Wang et al.
FIG. 4. The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩ for (A25 B25)2. The results for homopolymer chains A 50 and B50 are also shown for comparison. In all cases, ϵ AA = 1.0 and ϵ AB = ϵ BB = 0.2.
we can see that the Rg derivative of (A25 B25)2 increases again at low temperature, implying another peak in even low temperature region. To check the existence of this peak and elucidate its physical origin, we simulate this system with ITS method by including lower temperatures. The result of Rg derivative is shown in Fig. 4, and the results for strongly attractive A50 (i.e., ϵ AA = 1.0) and for weakly attractive B50 (i.e., ϵ BB = 0.2) are also shown for comparison. Now, we can clearly see two peaks for (A25 B25)2. These two peaks correspond to the C2G transitions of strongly attractive A50 and weakly attractive B50, respectively, since the positions of the peaks of multiblock copolymers are just the same as those of the respective homopolymer. Note that the transition temperature for the strongly attractive parts of (A25 B25)2 is a little lower than that of A50, which may be attributed to the fact that the extension of the hydrophilic blocks prevents the collapse of strongly attractive blocks, thus, the C2G transition occurs at relatively lower temperature. Therefore, for multiblock copolymers with large enough block sizes (such as (A25 B25)2), the strongly attractive blocks first go through a C2G transition and collapse into a compact core as the temperature decreases, but at this temperature, the weakly attractive blocks remain extended since the temperature is not low enough for them to collapse. When the
J. Chem. Phys. 141, 244907 (2014)
temperature decreases further, the weakly attractive blocks will also go through a C2G transition and form a globule state. Such two-step structural transition is not observed for multiblock copolymers with smaller block sizes. For multiblock copolymers with small block sizes, the strongly attractive blocks go through a C2G transition first with decreasing temperature. But after the strongly attractive core is formed, the weakly attractive blocks will be anchored locally on the surface of the strongly attractive core, which prevents further aggregation of weakly attractive blocks. Now, we focus on the influence of amphiphilicity between blocks on the order of C2G transitions. Two representative cases for varying the amphiphilicity between blocks of multiblock copolymers are considered. In the first case, we keep ϵ AA = 1.0 (manifesting A blocks are always more strongly attractive) and adjust ϵ BB and ϵ AB from 0.2 to 0.9. With increasing ϵ BB and ϵ AB, the whole polymer chain is changing gradually from an amphiphilic multiblock copolymer to a strongly attractive chain. In the second case, we keep ϵ BB = ϵ AB = 0.2 and adjust ϵ AA from 0.9 to 0.3. With decreasing ϵ AA, the whole polymer chain gradually changes from an amphiphilic multiblock copolymer to a weakly attractive chain. In both cases, we calculate Cv and the Binder cumulant to identify the order of the C2G transitions. Therefore, we can draw phase diagrams showing the influence of block size and amphiphilicity between blocks on the order of C2G transitions. The results are shown in Fig. 5, in which letter “F” means there is a first-order C2G transition, while letter “N” means the transition is not first-order. We can see that the C2G transition is strongly related to the block size and interaction strength of the block: only for suitable block size and apparent amphiphilicity of the multiblock copolymers, the C2G transitions are of first order. It implies that to mimic protein folding with multiblock copolymers, we need to design suitable block size and chain amphiphilicity, since the folding transition is first-order for typical proteins. Besides the multiblock copolymers with equal A and B compositions, we also study the C2G transitions of multiblock copolymers with different compositions with ITS method. The sequences of the chains are (A2 B8)10, (A3 B7)10, (A4 B6)10, (A6 B4)10, (A7 B3)10, and (A8 B2)10, i.e., the proportions of strongly attractive beads change from 0.2 to 0.8. The results are shown in Fig. 6. We can see that with the increase of the proportion of strongly attractive beads, the transition temperature increases. This is because with less
FIG. 5. The phase diagrams referring to the order of the C2G transitions of multiblock copolymers. “F” means first order and “N” means not. For the systems shown in the left diagram, ϵ AA = 1.0. For the systems shown in the right diagram, ϵ AB = ϵ BB = 0.2.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-6
Wang et al.
J. Chem. Phys. 141, 244907 (2014)
FIG. 6. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the specific heat C v versus temperature for multiblock copolymers with different proportions of hydrophobic beads.
FIG. 7. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the specific heat C v versus temperature for multiblock copolymers with different mutation sequences as listed in Table II. Sequence 0 corresponds to (A 25 B 25)2 without mutation.
weakly attractive components, the interactions of the beads in the whole chain turns to be more attractive, thus, the transition temperature approaches to that of the corresponding strongly attractive homopolymer A100. Meanwhile, as shown in Fig. 6(b), the C2G transition is no more of first order for multiblock copolymers with large strongly attractive proportion, since there is no obvious peak in Cv curves when the proportion of strongly attractive beads is larger than 0.6. We notice that the multiblock copolymer chains with short strongly attractive blocks always show a clear first-order C2G transition. As we know, the protein may not fold to its correct state due to the mutation of one single amino acid, which may further cause serious diseases.48 As the strongly attractive parts of protein play the dominant role in the folding process, we suppose mutation in strongly attractive blocks in our multiblock copolymers may influence the C2G transitions. Thus, we choose (A5 B5)10 as the standard sequence, and consider changing one (or several) strongly attractive bead(s) into weakly attractive, to study the influence of mutation on
C2G transitions. Specifically, we consider the sequences as shown in Table II after mutation. From sequence 1 to sequence 4, one, two, five, and ten A beads in the middle of each A blocks are changed to B type. The C2G transition results are shown in Fig. 7. We can see that a single mutation in one block can lead to the decrease of transition temperature; with the increase of the number of mutations, the decrease of transition temperature is more apparent. As the mutation occurs, the attractive interactions of the whole chain turn to be weaker, thus, lower temperature is needed to observe C2G transitions. We further study the influence of the mutation position in strongly attractive blocks on C2G transitions of (A5 B5)10. The specific sequences after mutations are shown in Table III. The proportions of strongly attractive beads for all the sequences are the same. The C2G transition results are shown in Fig. 8. It is clear that the transition temperatures for the first and the third sequences are basically the same, whereas the transition temperature for the second sequence is a little lower than the others. Comparing the sequences listed in Table III, we can
TABLE II. The sequences of multiblock copolymers with different mutations. Sequence index 1 2 3 4
TABLE III. The sequences of multiblock copolymers with different mutation positions.
Sequences (AABAABBBBB)1(AAAAABBBBB)9 (AABAABBBBB)2(AAAAABBBBB)8 (AABAABBBBB)5(AAAAABBBBB)5 (AABAABBBBB)10
Sequence index
Repeating blocks
1 2 3
(AAABABBBBB)10 (AABAABBBBB)10 (ABAAABBBBB)10
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-7
Wang et al.
J. Chem. Phys. 141, 244907 (2014)
ACKNOWLEDGMENTS
This work is subsidized by the National Basic Research Program of China (973 Program, 2012CB821500) and supported by National Science Foundation of China (21025416). We also thank the supports from Jilin Province Science and Technology Development Plan (20140519004JH). 1I.
Nishio, S.-T. Sun, G. Swislow, and T. Tanaka, Nature 281, 201 (1979). A. Moore, J. Phys. A: Math. Gen. 10, 305 (1977). 3I. Lifshitz, A. Y. Grosberg, and A. Khokhlov, Rev. Mod. Phys. 50, 683 (1978). 4V. S. Pande, A. Y. Grosberg, and T. Tanaka, Rev. Mod. Phys. 72, 259 (2000). 5E. Sherman and G. Haran, Proc. Natl. Acad. Sci. U. S. A. 103, 11539 (2006). 6A. Gutin, V. Abkevich, and E. Shakhnovich, Phys. Rev. Lett. 80, 208 (1998). 7A. Y. Grosberg and E. I. Shakhnovich, Sov. Phys. JETP 64, 1284 (1986). 8A. R. Khokhlov and P. G. Khalatur, Phys. A 249, 253 (1998). 9T. M. Birshtein and V. A. Pryamitsyn, Macromolecules 24, 1554 (1991). 10D. Yang and Q. Wang, ACS Macro Lett. 2, 952 (2013). 11R. Wang and Z.-G. Wang, Macromolecules 47, 4094 (2014). 12P. Grassberger and R. Hegger, J. Chem. Physica. 102, 6881 (1995). 13A. M. Rubio, J. J. Freire, J. H. R. Clarke, C. W. Yong, and M. Bishop, J. Chem. Phys. 102, 2277 (1995). 14A. Milchev, W. Paul, and K. Binder, J. Chem. Phys. 99, 4786 (1993). 15J. Ma, J. E. Straub, and E. I. Shakhnovich, J. Chem. Phys. 103, 2615 (1995). 16Y. Zhou, M. Karplus, J. M. Wichert, and C. K. Hall, J. Chem. Phys. 107, 10691 (1997). 17G. Tanaka and W. L. Mattice, Macromolecules 28, 1049 (1995). 18J. P. K. Doye, R. P. Sear, and D. Frenkel, J. Chem. Phys. 108, 2134 (1998). 19A. Irbärk and E. Sandelin, J. Chem. Phys. 110, 12256 (1999). 20H. Liang and H. Chen, J. Chem. Phys. 113, 4469 (2000). 21H. Noguchi and K. Yoshikawa, J. Chem. Phys. 109, 5070 (1998). 22J. M. Polson and N. E. Moore, J. Chem. Phys. 122, 024905 (2005). 23J. K. C. Suen, F. A. Escobedo, and J. J. de Pablo, J. Chem. Phys. 106, 1288 (1997). 24D. T. Seaton, T. Wüst, and D. P. Landau, Phys. Rev. E 81, 011802 (2010). 25Z. Wang, L. Wang, and X. He, Soft Matter 9, 3106 (2013). 26M. Bachmann and W. Janke, J. Chem. Phys. 120, 6779 (2004). 27M. Bachmann, H. Arkin, and W. Janke, Phys. Rev. E 71, 031906 (2005). 28S. Schnabel, M. Bachmann, and W. Janke, J. Chem. Phys. 131, 124904 (2009). 29J. E. Magee, J. Warwicker, and L. Lue, J. Chem. Phys. 120, 11285 (2004). 30K. A. Dill, Biochemistry 24, 1501 (1984). 31K. M. Fiebig and K. A. Dill, J. Chem. Phys. 98, 3475 (1993). 32K. Yue and K. A. Dill, Phys. Rev. E 48, 2267 (1993). 33F. Stillinger, T. Head-Gordon, and C. Hirshfeld, Phys. Rev. E 48, 1469 (1993). 34F. Stillinger and T. Head-Gordon, Phys. Rev. E 52, 2872 (1995). 35A. Irbäck, C. Peterson, F. Potthast, and O. Sommelius, J. Chem. Phys. 107, 273 (1997). 36Y. W. Li, T. Wust, and D. P. Landau, Commun. Comput. Phys. 182, 1896 (2011). 37T. Wust and D. P. Landau, J. Chem. Phys. 137, 064903 (2012). 38G. Shi, T. Vogel, T. Wüst, Y.-W. Li, and D. P. Landau, Phys. Rev. E 90, 033307 (2014). 39L. Yang, Q. Shao, and Y. Q. Gao, J. Chem. Phys. 130, 124111 (2009). 40P. Zhao, L. J. Yang, Y. Q. Gao, and Z.-Y. Lu, Chem. Phys. 415, 98 (2013). 41Y.-L. Zhu, H. Liu, Z.-W. Li, H.-J. Qian, G. Milano, and Z.-Y. Lu, J. Comput. Chem. 34, 2197 (2013). 42C. Wu and X. Wang, Phys. Rev. Lett. 80, 4092 (1998). 43W. Paul, T. Strauch, F. Rampf, and K. Binder, Phys. Rev. E 75, 060801 (2007). 44W. Burchard, M. Schmidt, and W. H. Stockmayer, Macromolecules 13, 580 (1980). 45B. Binder, Z. Phys. B: Condens. Matter 43, 119 (1981). 46W. Janke, Phys. Rev. B 47, 14757 (1993). 47W. Burchard, Adv. Polym. Sci. 48, 1 (1983). 48T. K. Chaudhuri and S. Paul, FEBS J. 273, 1331 (2006). 2M.
FIG. 8. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the specific heat C v versus temperature for multiblock copolymers with different mutation positions as shown in Table III.
find that the most influential position of mutation corresponds to the middle of each strongly attractive blocks. This result may be helpful for the modification of peptide sequence and further control its folding structures.
V. CONCLUSIONS
In summary, we find that the ITS method is an appropriate tool to study the C2G transitions of polymer chains. Especially for multiblock copolymers, we study the block size, amphiphilicity between blocks, and composition ratio between the two blocks on the C2G transition temperatures and the transition order. The transition temperature decreases with decreasing block size at a fixed chain length, and the transition temperature also decreases with increasing proportion of weakly attractive blocks. For homopolymers, the C2G transition is not of first order, but by choosing suitable block size and chain amphiphilicity of multiblock copolymer, the C2G transition can be tuned to be first-order. We also find that mutating the beads at middle positions of strongly attractive blocks to weakly attractive can largely change the C2G transition temperature. These results systematically illustrate the characteristics of C2G transitions of multiblock copolymers and may be helpful to design new functions in applications of multiblock copolymers.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
THE JOURNAL OF CHEMICAL PHYSICS 141, 244907 (2014)
Coil-to-globule transitions of homopolymers and multiblock copolymers Wei Wang, Peng Zhao, Xi Yang, and Zhong-Yuan Lua) Institute of Theoretical Chemistry, State Key Laboratory of Theoretical and Computational Chemistry, Jilin University, Changchun 130023, China
(Received 1 October 2014; accepted 9 December 2014; published online 30 December 2014) We study the coil-to-globule transitions of both homopolymers and multiblock copolymers using integrated tempering sampling method, which is a newly proposed enhanced sampling method that can efficiently sample the energy space with low computational costs. For homopolymers, the coil-to-globule structure transition temperatures (Ttr) are identified by the radius of gyration of the chain. The transition temperature shows a primary scaling dependence on the chain length (N) with Ttr ∼ N −1/2. For multiblock copolymers, the coil-to-globule transition can be identified as first order, depending on the block size and the difference in attractive interactions of blocks. The influence of mutating a small portion of strongly attractive blocks to weakly attractive blocks on the coil-to-globule transition is found to be related to the position of the mutation. C 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4904888]
I. INTRODUCTION
The change of a single polymer chain from an extended coil in good solvent to a collapsed globule in poor solvent is called coil-to-globule (C2G) transition.1 This collapse typically corresponds to a second-order phase transition for homopolymers with suitable lengths.2,3 The importance of the C2G transition is due to its inherent connection with the problem of tuning function of a protein.4,5 The key factor influencing protein function strongly relates to a single polypeptide chain folding from coil state to globule state with specific three-dimensional structures. In this process, the protein, in which some amino acid residues are hydrophilic or charged while others are hydrophobic,4,5 undergoes a firstorder folding transition to reach the native state.4–6 Thus, heteropolymers or even block copolymer models (resembling protein structures) are widely used to study the protein folding for a long time, and the transition process together with C2G transition of both homopolymer and heteropolymer chain had been explained theoretically.4,7–11 Computer simulations can provide more structural, thermodynamic, and dynamic information for C2G transition of a polymer chain in dilute solution. Both Monte Carlo (MC) and molecular dynamics (MD) methods with either implicit or explicit solvent models have been used to elucidate the influence of chain length, chain stiffness, confinement, etc., on the C2G transition temperature, order of transition, and even transition dynamics.12–26 But for long chains, the energy landscape of the system is always too rugged with many local minima and barriers, thus, the ergodicity may not be guaranteed in normal simulations. To efficiently sample the configuration space, enhanced sampling methods, such as Wang-Landau methods24,25 and multicanonical ensemble methods,26–28 are usually adopted in the studies of C2G a)Electronic address: [email protected]
0021-9606/2014/141(24)/244907/7/$30.00
transitions of polymer chains. With the aid of enhanced sampling methods, C2G transitions of homopolymers and several sequences of heteropolymers were studied for giving hints of the formation of protein native structures.27,29 The multiblock copolymers, which can be tuned to homopolymers if all the blocks have the same property or to polypeptides if the block hydrophobicities are designed to mimic specific protein sequence, can serve as a versatile system bridging homopolymer model and heteropolymer model in the studies of C2G transition. The influence of inherent driving force for protein folding can be well illustrated in multiblock copolymer models by designing block sizes, positions, and interactions in a chain. The most prominent examples are the lattice HP model30–32 and the off-lattice AB model,33,34 which had been exhaustively investigated.24,27,35–37 Moreover, previously overlooked influence of mutation38 in strongly attractive blocks can be systematically studied with a regularly distributed multiblock copolymer model. Therefore, focusing on these issues, we study the C2G transitions of regularly distributed block copolymers and compare the results to those of homopolymers with a novel enhanced sampling method–integrated tempering sampling (ITS)39,40 implemented in GPU MD code GALAMOST41 in this work. This paper is organized as follows: in Sec. II, the models of homopolymer and block copolymers are described. The ITS method will be briefly introduced in Sec. III, in which the main idea and the most important formulae are given. The MD simulation details are also illustrated in this section. The results and discussion are shown in Sec. IV, followed by Sec. V for conclusions.
II. MODELS
Common bead-spring models are used in this study to describe homopolymers and block copolymers. The homopolymer is a linear chain with one type of beads
141, 244907-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-2
Wang et al.
J. Chem. Phys. 141, 244907 (2014)
connected by bonds with finite extensible nonlinear elastic (FENE) potential between connected beads
representing the contribution from each temperature Tk . The biased potential U ′ is defined through
1 r2 − Kr 02 ln *1.0 − 2 + +UWCA r ≤ r 0 UFENE(r) = 2 , (1) r0 , ∞ r > r0 where ( 1 σWCA ) 12 ( σWCA ) 6 4ϵ WCA − r ≤ 2 6 σWCA r r UWCA(r) = , 1 0 6 σWCA r > 2 (2)
e−βU (r ) ≡ W (r),
(5)
1 U ′(r) = − ln nk e−β kU (r ). β k
(6)
and r 0 is the bond extension parameter, K is the force constant, and r is the instantaneous bond length. To be consistent with previous work,40 we set σWCA = 1.05, ϵ WCA = 1.0, r 0 = 1.5, and K = 40. The Lennard-Jones potential is used between non-bonded beads with σ = 1.0 and ϵ = 1.0 (which, together with bead mass m ≡ 1.0, are taken as the units in this study), ( ) 12 ( ) 6 σ σ − r ≤ r cut, 4ϵ r r , (3) ULJ (r) = 0 r > r cut with r cut being set to 3.0 in our simulations. We consider homopolymers with chain lengths (N) between 13 and 800. The solvent effects are treated by implicit solvent model. For block copolymers, we use two types of beads denoted by A and B, respectively. The linear block copolymer chain with blocks of A and B type beads are connected by FENE potential. To represent the incompatibility between A and B beads, we use Lennard-Jones potential but with different energy parameter ϵ to characterize the strength of short-range attractive interaction. To describe the amphiphilicity of the block copolymer chain, we set ϵ AA = 1.0 and ϵ BB = 0.2 in most of the cases. The interaction parameter between A and B beads is set the same as that between B type beads, i.e., ϵ AB = ϵ BB = 0.2, which is consistent with the “AB” model used in the literature.27 To represent the effect of changing block copolymer amphiphilicity, we keep the value of ϵ AB and ϵ BB the same and adjust the value in the range of 0.2–0.8.
III. SIMULATION METHOD
The integrated tempering sampling method is a novel easy-to-use enhanced sampling method based on generalized ensemble to generate a distribution covering a broad range of energies.39,40 Here, we only give a brief introduction to this method, since details of the calculation process have been introduced in previous work.40 The main idea in ITS method is to use a biased potential generated from the combination of a series of canonical energy distributions. This allows efficient sampling in a desired temperature and energy range without requiring a predefined reaction coordinate. The generalized distribution W (r) is defined as W (r) = nk e−β kU (r ), (4) k
where U is potential energy, βk = 1/k BT (k B is Boltzmann constant and T is temperature), and nk is the weighting factor
′
thus,
The biased force Fb is simply the force in conventional MD simulations, F, carrying a pre-factor ∂U ′ ∂U ′ ∂U nk βk e−β kU (r ) F, (7) Fb = − =− = k ∂r ∂U ∂r β k nk e−β kU (r ) thus, ITS method can be readily implemented in any MD code by slightly modifying the force subroutine. The thermodynamic properties of any canonical ensemble whose temperature (T j ) is in the desired range can be calculated from the generalized ensemble by using reweighting A(r )e −β j U (r ) W (r)dr A(r)e−β jU (r )dr W (r ) ⟨A⟩ β j = −β U (r ) = −β U (r ) j e e j dr W (r ) W (r)dr −β j U (r ) =
A(r )e W (r )
e
−β j U (r )
W (r )
W
.
(8)
W
Typical computational procedure of ITS simulation is therefore as follows: 1. Determine the desired temperature range. 2. Choose a set of temperatures in the desired range and run several short time conventional MD simulations. Then, the relation between potential energy and temperature can be obtained by interpolation. Here, we define an overlap factor t, which can actually be determined by the requirement that the ratio between energy probability density functions at two adjacent temperatures is a constant. Giving a proper value of the overlap factor t, sufficient overlap of the energy distribution between adjacent temperatures can be guaranteed. 3. Determine the ITS temperature distribution and the corresponding weighting factors nk as described in Ref. 40. 4. Use the parameters generated in steps 2 and 3 to perform ITS simulation, which is essentially nothing but a conventional MD simulation using biased force calculated by Eq. (7). 5. After ITS simulation, the canonical ensemble properties can be calculated by Eq. (8). Nosè-Hoover thermostat is adopted in our ITS simulations. The potential energies at different temperatures are first obtained by 1.0 × 106 steps MD simulations at 8 temperatures in the range of 1.0–5.0 for homopolymers and in the range of 0.5–5.0 for block copolymers. The dependence of potential energy on temperature is used to estimate nk .40 The overlap factor t is set to e0.5, thus, 8–83 temperatures are generated in the range of 1.0–4.0 for systems with different chain lengths.40 The specific numbers of temperatures are shown in Table I. For homopolymers with different chain lengths, we then perform 1.0 × 109 steps ITS simulations with GALAMOST package41 at each chain length, and the data are recorded every 103 steps. Using a single ITS simulation
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-3
Wang et al.
J. Chem. Phys. 141, 244907 (2014)
TABLE I. The number of temperatures in ITS simulations. Chain length
Number of temperatures
13 25 50 100 200 400 800
8 12 18 27 40 58 83
trajectory, we can obtain thermodynamic properties of the system at any temperature in the temperature range. For multiblock copolymers, we exclusively focus on the influence of block size, strength of attraction, and mutation on the C2G transitions, so only chain with length N = 100 is considered, and the numbers of temperatures to perform ITS simulations are 25.
IV. RESULTS AND DISCUSSION A. Homopolymers
We calculate the radius of gyration (Rg ) of the chain and its derivative and the specific heat (Cv ) as structural and thermodynamic indicators, respectively, for C2G transitions for homopolymer with N = 13 ∼ 800. The radius of gyration (Rg ) is defined as
Rg2
N N 1 2 = r , 2N 2 i=1 j=1 i j
and the specific heat (Cv ) is defined as
2 E − ⟨E⟩2 ∂U . = Cv = ∂T k BT 2
(9)
(10)
The dependence of both Rg and Cv on temperature for homopolymers with different lengths is shown in Fig. 1. In Fig. 1(a), we can see the peaks in the derivative of Rg for all the chains, which clearly indicate strong structure changes corresponding to C2G transition. In Fig. 1(b), we can see for N < 200, there is no clear signal in the Cv curves for C2G transition, which implies strong finite size effect. While for N > 200, Cv curves have shoulders at the temperatures related to the peaks in the curves of the first derivative of Rg in Fig. 1(a), which can be recognized as the energetic signal of the C2G transition, similar to those reported in Ref. 24. No sharp peak in the Cv curves implies that the C2G transition for homopolymers is not thermodynamically first-order, which is consistent with previous conclusions.6 We can also clearly see the transition temperatures (Ttr) identified from Fig. 1(a) move to higher value with increasing chain length. To further elucidate the dependence of Ttr on the chain length, we fit Ttr versus N with the following equation: √ (11) Ttr(N) −TΘ = a1/ N + a2/N, where Ttr(N) is the C2G transition temperature for each chain length and TΘ is the Θ temperature characterizing the turning
FIG. 1. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the specific heat C v per monomer versus temperature.
point of structure transition of a polymer chain.3 Here, the √ N term corresponds to a mean field fitting, and the N term is adopted for finite chain length correction. As shown in Fig. 2, the fitting is quite acceptable with parameters a1 = −12.43, a2 = 9.58, and TΘ = 3.696, which are in agreement with those reported by using Wang-Landau method.24 Combining the results shown above for homopolymers, we can verify that ITS is an effective method for elucidating C2G transitions. Fig. 1(b) also shows rugged regions at low temperatures especially when the chain length is larger than 200. It may be
FIG. 2. The dependence of the transition temperature on homopolymer chain length with N from 13 to 800. The dashed line shows fitting with Eq. (11).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-4
Wang et al.
attributed to the existence of metastable intermediate states for long chains that have been identified in experiments42 and in simulations,43 but it also may be attributed to the inefficiency of sampling with ITS method for long chains. However, since we are focusing C2G transition (which takes place at comparatively higher temperature) only, we will not pay much attention to these rugged regions at low temperatures in the following studies. B. Heteropolymers
In this section, we exclusively study the influence of block size, amphiphilicity between blocks, and mutation at specific sites on the C2G transitions of multiblock copolymers. Since changing chain length only results in shift of C2G transition temperature, in the following, we will keep the chain length constant at N = 100 and consider four types of multiblock copolymers with different block sizes: (A2 B2)25, (A5 B5)10, (A10 B10)5, and (A25 B25)2. We keep ϵ AA = 1.0 to manifest A blocks strongly attractive and ϵ BB = ϵ AB = 0.2 to manifest B blocks weakly attractive and incompatible with A blocks. We calculate the first derivative of radius of gyration and Cv at different temperatures, as shown in Figs. 3(a) and 3(c), respectively. To further clarify the structure transition characteristics of multiblock copolymer, we also calculate the shape factor ρ, defined as ρ = Rg /Rh , where Rg is the radius of gyration, and Rh is the strongly attractive radius N N
−1 1 r−1 . (12) Rh = 2 N i=1 j=1(i, j) i j The shape factor ρ was widely used in light scattering techniques to probe the shape, size, or structure of certain molecule systems.44,47 Thus, it is used here to measure the conformational change when C2G transitions of single polymer chains occur. Generally, the value of ρ is about 1.5 for flexible chains in good solvent, while ρ = 0.77 for a solid sphere.44,47 In Fig. 3(b), we can see that the shape factor
J. Chem. Phys. 141, 244907 (2014)
is around 1.2 for chains at high temperature (corresponding to good solvent condition), which is a little lower than the theoretical value for coil state in good solvent. It is mainly attributed to the strongly attractive parts of the multiblock copolymer chain, which make the conformation of the chain not extended as much as an ideal coil. However, the C2G transitions can be clearly discerned from the turning points in the curves. To further justify the order of the transition other than using Cv , we also calculate the Binder cumulant VE 45 for all the chains at different temperatures. The Binder cumulant has an isolated minimum VE < 2/3 at the transition point,46 thus can be a good energy indicator. VE is defined as
4 E (13) VE = 1 − 2 . 3 E2 We can see from Figs. 3(a) and 3(b) that both the first derivative of mean square of radius of gyration and the shape factor give clear signals of the C2G transition, and the transition temperature increases with increasing block size. As compared to homopolymers, Fig. 3(c) shows a clear difference: obvious peaks from the Cv curves of multiblock copolymers can be observed, while there are only shoulders in the cases of homopolymers (Fig. 1(b)). This implies that C2G transition is first order for certain multiblock copolymers studied here. The Binder cumulant in Fig. 3(d) also gives consistent results, showing valleys at the same positions as compared to the peaks in Cv curves. From Fig. 3, we can also find that, the transition temperatures indicated by the peaks of the structural indicators are a little bit higher than those from energetic indicators. It may be due to the fact that the response of energetic properties is not sensitive as compared to the structural properties. The peaks in the structural curves correspond to the maximum changes of the molecular structures, while the energy changes may show apparent signal when the polymer chains are compact enough (at lower temperatures). It should be noted that the shape factor of (A25 B25)2 has two turning points, as shown in Fig. 3(b). From Fig. 3(a),
FIG. 3. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the shape factor ρ versus temperature; (c) the specific heat C v versus temperature; (d) the Binder cumulant versus temperature for multiblock copolymer chains with different block sizes. In all cases, ϵ AA = 1.0 and ϵ AB = ϵ BB = 0.2.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-5
Wang et al.
FIG. 4. The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩ for (A25 B25)2. The results for homopolymer chains A 50 and B50 are also shown for comparison. In all cases, ϵ AA = 1.0 and ϵ AB = ϵ BB = 0.2.
we can see that the Rg derivative of (A25 B25)2 increases again at low temperature, implying another peak in even low temperature region. To check the existence of this peak and elucidate its physical origin, we simulate this system with ITS method by including lower temperatures. The result of Rg derivative is shown in Fig. 4, and the results for strongly attractive A50 (i.e., ϵ AA = 1.0) and for weakly attractive B50 (i.e., ϵ BB = 0.2) are also shown for comparison. Now, we can clearly see two peaks for (A25 B25)2. These two peaks correspond to the C2G transitions of strongly attractive A50 and weakly attractive B50, respectively, since the positions of the peaks of multiblock copolymers are just the same as those of the respective homopolymer. Note that the transition temperature for the strongly attractive parts of (A25 B25)2 is a little lower than that of A50, which may be attributed to the fact that the extension of the hydrophilic blocks prevents the collapse of strongly attractive blocks, thus, the C2G transition occurs at relatively lower temperature. Therefore, for multiblock copolymers with large enough block sizes (such as (A25 B25)2), the strongly attractive blocks first go through a C2G transition and collapse into a compact core as the temperature decreases, but at this temperature, the weakly attractive blocks remain extended since the temperature is not low enough for them to collapse. When the
J. Chem. Phys. 141, 244907 (2014)
temperature decreases further, the weakly attractive blocks will also go through a C2G transition and form a globule state. Such two-step structural transition is not observed for multiblock copolymers with smaller block sizes. For multiblock copolymers with small block sizes, the strongly attractive blocks go through a C2G transition first with decreasing temperature. But after the strongly attractive core is formed, the weakly attractive blocks will be anchored locally on the surface of the strongly attractive core, which prevents further aggregation of weakly attractive blocks. Now, we focus on the influence of amphiphilicity between blocks on the order of C2G transitions. Two representative cases for varying the amphiphilicity between blocks of multiblock copolymers are considered. In the first case, we keep ϵ AA = 1.0 (manifesting A blocks are always more strongly attractive) and adjust ϵ BB and ϵ AB from 0.2 to 0.9. With increasing ϵ BB and ϵ AB, the whole polymer chain is changing gradually from an amphiphilic multiblock copolymer to a strongly attractive chain. In the second case, we keep ϵ BB = ϵ AB = 0.2 and adjust ϵ AA from 0.9 to 0.3. With decreasing ϵ AA, the whole polymer chain gradually changes from an amphiphilic multiblock copolymer to a weakly attractive chain. In both cases, we calculate Cv and the Binder cumulant to identify the order of the C2G transitions. Therefore, we can draw phase diagrams showing the influence of block size and amphiphilicity between blocks on the order of C2G transitions. The results are shown in Fig. 5, in which letter “F” means there is a first-order C2G transition, while letter “N” means the transition is not first-order. We can see that the C2G transition is strongly related to the block size and interaction strength of the block: only for suitable block size and apparent amphiphilicity of the multiblock copolymers, the C2G transitions are of first order. It implies that to mimic protein folding with multiblock copolymers, we need to design suitable block size and chain amphiphilicity, since the folding transition is first-order for typical proteins. Besides the multiblock copolymers with equal A and B compositions, we also study the C2G transitions of multiblock copolymers with different compositions with ITS method. The sequences of the chains are (A2 B8)10, (A3 B7)10, (A4 B6)10, (A6 B4)10, (A7 B3)10, and (A8 B2)10, i.e., the proportions of strongly attractive beads change from 0.2 to 0.8. The results are shown in Fig. 6. We can see that with the increase of the proportion of strongly attractive beads, the transition temperature increases. This is because with less
FIG. 5. The phase diagrams referring to the order of the C2G transitions of multiblock copolymers. “F” means first order and “N” means not. For the systems shown in the left diagram, ϵ AA = 1.0. For the systems shown in the right diagram, ϵ AB = ϵ BB = 0.2.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-6
Wang et al.
J. Chem. Phys. 141, 244907 (2014)
FIG. 6. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the specific heat C v versus temperature for multiblock copolymers with different proportions of hydrophobic beads.
FIG. 7. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the specific heat C v versus temperature for multiblock copolymers with different mutation sequences as listed in Table II. Sequence 0 corresponds to (A 25 B 25)2 without mutation.
weakly attractive components, the interactions of the beads in the whole chain turns to be more attractive, thus, the transition temperature approaches to that of the corresponding strongly attractive homopolymer A100. Meanwhile, as shown in Fig. 6(b), the C2G transition is no more of first order for multiblock copolymers with large strongly attractive proportion, since there is no obvious peak in Cv curves when the proportion of strongly attractive beads is larger than 0.6. We notice that the multiblock copolymer chains with short strongly attractive blocks always show a clear first-order C2G transition. As we know, the protein may not fold to its correct state due to the mutation of one single amino acid, which may further cause serious diseases.48 As the strongly attractive parts of protein play the dominant role in the folding process, we suppose mutation in strongly attractive blocks in our multiblock copolymers may influence the C2G transitions. Thus, we choose (A5 B5)10 as the standard sequence, and consider changing one (or several) strongly attractive bead(s) into weakly attractive, to study the influence of mutation on
C2G transitions. Specifically, we consider the sequences as shown in Table II after mutation. From sequence 1 to sequence 4, one, two, five, and ten A beads in the middle of each A blocks are changed to B type. The C2G transition results are shown in Fig. 7. We can see that a single mutation in one block can lead to the decrease of transition temperature; with the increase of the number of mutations, the decrease of transition temperature is more apparent. As the mutation occurs, the attractive interactions of the whole chain turn to be weaker, thus, lower temperature is needed to observe C2G transitions. We further study the influence of the mutation position in strongly attractive blocks on C2G transitions of (A5 B5)10. The specific sequences after mutations are shown in Table III. The proportions of strongly attractive beads for all the sequences are the same. The C2G transition results are shown in Fig. 8. It is clear that the transition temperatures for the first and the third sequences are basically the same, whereas the transition temperature for the second sequence is a little lower than the others. Comparing the sequences listed in Table III, we can
TABLE II. The sequences of multiblock copolymers with different mutations. Sequence index 1 2 3 4
TABLE III. The sequences of multiblock copolymers with different mutation positions.
Sequences (AABAABBBBB)1(AAAAABBBBB)9 (AABAABBBBB)2(AAAAABBBBB)8 (AABAABBBBB)5(AAAAABBBBB)5 (AABAABBBBB)10
Sequence index
Repeating blocks
1 2 3
(AAABABBBBB)10 (AABAABBBBB)10 (ABAAABBBBB)10
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
244907-7
Wang et al.
J. Chem. Phys. 141, 244907 (2014)
ACKNOWLEDGMENTS
This work is subsidized by the National Basic Research Program of China (973 Program, 2012CB821500) and supported by National Science Foundation of China (21025416). We also thank the supports from Jilin Province Science and Technology Development Plan (20140519004JH). 1I.
Nishio, S.-T. Sun, G. Swislow, and T. Tanaka, Nature 281, 201 (1979). A. Moore, J. Phys. A: Math. Gen. 10, 305 (1977). 3I. Lifshitz, A. Y. Grosberg, and A. Khokhlov, Rev. Mod. Phys. 50, 683 (1978). 4V. S. Pande, A. Y. Grosberg, and T. Tanaka, Rev. Mod. Phys. 72, 259 (2000). 5E. Sherman and G. Haran, Proc. Natl. Acad. Sci. U. S. A. 103, 11539 (2006). 6A. Gutin, V. Abkevich, and E. Shakhnovich, Phys. Rev. Lett. 80, 208 (1998). 7A. Y. Grosberg and E. I. Shakhnovich, Sov. Phys. JETP 64, 1284 (1986). 8A. R. Khokhlov and P. G. Khalatur, Phys. A 249, 253 (1998). 9T. M. Birshtein and V. A. Pryamitsyn, Macromolecules 24, 1554 (1991). 10D. Yang and Q. Wang, ACS Macro Lett. 2, 952 (2013). 11R. Wang and Z.-G. Wang, Macromolecules 47, 4094 (2014). 12P. Grassberger and R. Hegger, J. Chem. Physica. 102, 6881 (1995). 13A. M. Rubio, J. J. Freire, J. H. R. Clarke, C. W. Yong, and M. Bishop, J. Chem. Phys. 102, 2277 (1995). 14A. Milchev, W. Paul, and K. Binder, J. Chem. Phys. 99, 4786 (1993). 15J. Ma, J. E. Straub, and E. I. Shakhnovich, J. Chem. Phys. 103, 2615 (1995). 16Y. Zhou, M. Karplus, J. M. Wichert, and C. K. Hall, J. Chem. Phys. 107, 10691 (1997). 17G. Tanaka and W. L. Mattice, Macromolecules 28, 1049 (1995). 18J. P. K. Doye, R. P. Sear, and D. Frenkel, J. Chem. Phys. 108, 2134 (1998). 19A. Irbärk and E. Sandelin, J. Chem. Phys. 110, 12256 (1999). 20H. Liang and H. Chen, J. Chem. Phys. 113, 4469 (2000). 21H. Noguchi and K. Yoshikawa, J. Chem. Phys. 109, 5070 (1998). 22J. M. Polson and N. E. Moore, J. Chem. Phys. 122, 024905 (2005). 23J. K. C. Suen, F. A. Escobedo, and J. J. de Pablo, J. Chem. Phys. 106, 1288 (1997). 24D. T. Seaton, T. Wüst, and D. P. Landau, Phys. Rev. E 81, 011802 (2010). 25Z. Wang, L. Wang, and X. He, Soft Matter 9, 3106 (2013). 26M. Bachmann and W. Janke, J. Chem. Phys. 120, 6779 (2004). 27M. Bachmann, H. Arkin, and W. Janke, Phys. Rev. E 71, 031906 (2005). 28S. Schnabel, M. Bachmann, and W. Janke, J. Chem. Phys. 131, 124904 (2009). 29J. E. Magee, J. Warwicker, and L. Lue, J. Chem. Phys. 120, 11285 (2004). 30K. A. Dill, Biochemistry 24, 1501 (1984). 31K. M. Fiebig and K. A. Dill, J. Chem. Phys. 98, 3475 (1993). 32K. Yue and K. A. Dill, Phys. Rev. E 48, 2267 (1993). 33F. Stillinger, T. Head-Gordon, and C. Hirshfeld, Phys. Rev. E 48, 1469 (1993). 34F. Stillinger and T. Head-Gordon, Phys. Rev. E 52, 2872 (1995). 35A. Irbäck, C. Peterson, F. Potthast, and O. Sommelius, J. Chem. Phys. 107, 273 (1997). 36Y. W. Li, T. Wust, and D. P. Landau, Commun. Comput. Phys. 182, 1896 (2011). 37T. Wust and D. P. Landau, J. Chem. Phys. 137, 064903 (2012). 38G. Shi, T. Vogel, T. Wüst, Y.-W. Li, and D. P. Landau, Phys. Rev. E 90, 033307 (2014). 39L. Yang, Q. Shao, and Y. Q. Gao, J. Chem. Phys. 130, 124111 (2009). 40P. Zhao, L. J. Yang, Y. Q. Gao, and Z.-Y. Lu, Chem. Phys. 415, 98 (2013). 41Y.-L. Zhu, H. Liu, Z.-W. Li, H.-J. Qian, G. Milano, and Z.-Y. Lu, J. Comput. Chem. 34, 2197 (2013). 42C. Wu and X. Wang, Phys. Rev. Lett. 80, 4092 (1998). 43W. Paul, T. Strauch, F. Rampf, and K. Binder, Phys. Rev. E 75, 060801 (2007). 44W. Burchard, M. Schmidt, and W. H. Stockmayer, Macromolecules 13, 580 (1980). 45B. Binder, Z. Phys. B: Condens. Matter 43, 119 (1981). 46W. Janke, Phys. Rev. B 47, 14757 (1993). 47W. Burchard, Adv. Polym. Sci. 48, 1 (1983). 48T. K. Chaudhuri and S. Paul, FEBS J. 273, 1331 (2006). 2M.
FIG. 8. (a) The first derivative of the mean square radius of gyration versus temperature ⟨dR 2g /dT ⟩; (b) the specific heat C v versus temperature for multiblock copolymers with different mutation positions as shown in Table III.
find that the most influential position of mutation corresponds to the middle of each strongly attractive blocks. This result may be helpful for the modification of peptide sequence and further control its folding structures.
V. CONCLUSIONS
In summary, we find that the ITS method is an appropriate tool to study the C2G transitions of polymer chains. Especially for multiblock copolymers, we study the block size, amphiphilicity between blocks, and composition ratio between the two blocks on the C2G transition temperatures and the transition order. The transition temperature decreases with decreasing block size at a fixed chain length, and the transition temperature also decreases with increasing proportion of weakly attractive blocks. For homopolymers, the C2G transition is not of first order, but by choosing suitable block size and chain amphiphilicity of multiblock copolymer, the C2G transition can be tuned to be first-order. We also find that mutating the beads at middle positions of strongly attractive blocks to weakly attractive can largely change the C2G transition temperature. These results systematically illustrate the characteristics of C2G transitions of multiblock copolymers and may be helpful to design new functions in applications of multiblock copolymers.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Sat, 28 Feb 2015 13:31:37
