Collapse–Swelling Transitions of a Thermoresponsive, Single Poly(N

Dec 5, 2016 - B , 2016, 120 (51), pp 13184–13192 ... Three independent simulations with durations of 1 μs are performed at each temperature ... by ...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Collapse−Swelling Transitions of a Thermoresponsive, Single Poly(N‑isopropylacrylamide) Chain in Water Yunwon Kang, Heesun Joo, and Jun Soo Kim* Department of Chemistry and Nanoscience, Ewha Womans University, Seoul 03760, Republic of Korea S Supporting Information *

ABSTRACT: We present molecular dynamics (MD) simulations of a single poly(N-isopropylacrylamide) (PNIPAM) chain in explicit water at temperatures between 270 and 320 K near the lower critical solution temperature (LCST). The force-fields of OPLS-AA and TIP4P/2005 are used for a PNIPAM chain and water molecules, respectively. Three independent simulations with durations of 1 μs are performed at each temperature for a 30-mer PNIPAM chain starting with three distinct conformations: extended, loosely collapsed, and tightly collapsed states. The simulation trajectories exhibit reversible conformational transitions between swollen- and collapsed-chain conformations, which has rarely been reported in previous simulation studies, with the overall transition occurring at different temperatures depending on the initial conformation. The inconsistency of the transition temperatures depending on the initial conformation implies that, in spite of the simulation duration of 1 μs distinctly longer than that in previous simulation studies, the conformational sampling from the MD simulations is not enough to draw conclusions on equilibrium properties. Instead of evaluating average properties, therefore, the focus is on dynamic changes in the chain conformation during reversible collapse−swelling transitions at each temperature. The simulation trajectories are analyzed in terms of the radius of gyration, intrachain distances, hydrophobic contacts, and chain−water and intrachain hydrogen bonding. In particular, the formation of stable intrachain hydrogen bonds is a signature of the tightly collapsed-chain conformations that persist, once formed, for the entire simulation duration.

■

tacticity32 of a PNIPAM chain as well as the salt effect26 on the collapse transition were also investigated. However, when the MD simulations were started with a precollapsed-chain conformation, the swelling of the precollapsed chain at temperatures below the estimated LCST was not observed within the simulation time of up to 300 ns in previous MD studies.25,31 This may be either due to the slow relaxation of the chain conformation at low temperatures or due to the collapsed conformation being trapped in a local energy minimum, possibly implying the insufficient conformational sampling in either way. In this work, we perform 3 μs MD simulations of a single, 30mer PNIPAM chain at different temperatures ranging from 270 to 300 K in the explicit presence of water molecules. The first 1 μs simulation is run starting from an extended-chain conformation, whereas the other two 1 μs simulations are run from two distinct collapsed-chain conformations, referred to as loosely collapsed and tightly collapsed conformations. Distinction of the two collapsed conformations is based on the stability of the intrachain hydrogen bonding. For the system with an extended conformation, we perform additional 1 μs

INTRODUCTION Poly(N-isopropylacrylamide) (PNIPAM) is a smart polymer that gained huge popularity in the last three decades due to its temperature-responsive coil−globule transition or volume phase transition near physiological temperature with the lower critical solution temperature (LCST) around 32 °C.1−4 At temperatures below the LCST, PNIPAM is soluble in water with swollen polymer conformations, whereas it becomes collapsed or aggregated at temperatures above the LCST. The transition of PNIPAM near the physiological temperature enables the versatile utility of PNIPAM-based materials, in particular, in various biomedical applications.5−9 Temperature-responsive coil−globule transitions or volume phase transition of PNIPAM and PNIPAM-containing copolymers have been examined using various experimental methods including turbidimetry,3,10,11 light scattering,11−14 calorimetry,15,16 and nuclear magnetic resonance,17,18 infrared,19,20 and fluorescence21 spectroscopies. Complementary to the experimental observation of the LCST, there have been recent attempts to characterize the collapse−swelling transition of a single PNIPAM chain in water using molecular dynamics (MD) simulations.22−32 These studies concluded that the collapse transition of a stretched PNIPAM chain at various temperatures provided information on the transition temperature in qualitative agreement with the experimentally estimated LCST. The influences of the polymer length28 and © XXXX American Chemical Society

Received: September 10, 2016 Revised: December 1, 2016 Published: December 5, 2016 A

DOI: 10.1021/acs.jpcb.6b09165 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

simulations). Summary and conclusions are presented at the end.

simulations at 310 and 320 K to confirm the collapse transition of the extended chain. The simulation duration of 1 μs in this work is the longest of the MD simulation trajectories with the atomistic resolution of a PNIPAM chain. If the failure to observe the swelling of a precollapsed chain even at temperatures below the estimated LCST in the previous simulation studies arises from the slow conformational relaxation, then extension of the simulation duration to 1 μs should allow the observation of the reversible swelling−collapse transitions near the LCST. Indeed, the repetitive transitions between swollen and collapsed conformations are observed, although limited to a few, in our simulations starting with extended and loosely collapsed initial conformations. This implies that the conformational relaxation is so slow that the reversible collapse−swelling transitions of a single, 30mer PNIPAM chain are observed only when the simulation duration is extended, at least, to 1 μs. However, the swelling of a precollapsed-chain conformation is not observed within the simulation duration of 1 μs in our simulations with the tightly collapsed initial conformations. Thus, it may be interpreted that the tightly collapsed-chain conformations are trapped in a local energy minimum. In addition, the overall transition between swollen and collapsed conformations occurs at different temperatures, depending on different initial conformations, specifically, at temperatures between 300 and 310 K for the extended-chain simulations, between 280 and 290 K for the loosely collapsed-chain simulations, and below 270 K for the tightly collapsed-chain simulations. As a whole, these results imply that the simulation duration of 1 μs is not long enough to provide meaningful information on the equilibrium behavior of a PNIPAM chain, and the transition temperatures mentioned just above should not be taken as a thermodynamic transition temperature for PNIPAM. Nevertheless, our simulations are long enough to exhibit the reversible collapse−swelling transitions in the simulation trajectories that started with extended and loosely collapsed conformations at temperatures below the experimentally estimated transition temperature of 305 K. During these transitions, we examine the dynamic changes in intrachain distances, hydrophobic contacts, and hydrogen bonding. The changes in intrachain distances provide better information about the conformational states than the mere change in the radius of gyration. In previous MD studies, the hydrogenbonding and hydrophobic interactions were compared at different temperatures above and below the estimated LCST.22,24,27,29−32 In this work, however, we compare the changes in the interactions during the reversible collapse− swelling at the same temperature as well as at different temperatures. The results show that the collapse (or swelling) of a PNIPAM chain is correlated with an increase (or decrease) in intrachain hydrogen bonding and hydrophobic contacts but with a decrease (or increase) in chain−water hydrogen bonding. In particular, it is found that the stability of the intrachain hydrogen bonding can vary for the collapsed chains, by which the collapsed conformations are distinguished into loosely and tightly collapsed conformations. The rest of this article is organized as follows. We describe the simulation model and method in the next section. Simulation results are then presented and discussed, where the time-dependent changes in polymer conformations are investigated at temperatures ranging between 270 and 300 K (additionally at 310 and 320 K for the extended-chain

■

SIMULATION METHODS The collapse−swelling transition of a single 30-mer PNIPAM chain is investigated by performing MD simulations in the temperature range from 270 to 320 K. An isotactic arrangement of side groups is considered for chain tacticity. Three sets of 1 μs simulations are run with different PNIPAM chain conformations: extended, loosely collapsed, and tightly collapsed conformations, as shown in Figure 1 and summarized

Figure 1. Structure and conformations of a PNIPAM chain. (a) Chemical structure of a PNIPAM chain with n indicating the number of monomers. (b1−b3) Three initial conformations of a PNIPAM chain (with n = 30) in extended, loosely collapsed, and tightly collapsed states, respectively. Water molecules are omitted for visual clarity, and the snapshots are created using visual molecular dynamics.33

Table 1. Average Radius of Gyration (nm), Average Number of Intrachain Hydrogen Bonds, and Average Number of Hydrophobic Contacts of Three Initial Conformations, Estimated in the First 1 ns Simulations at 300 Ka

a

initial chain conformation

radius of gyration (nm)

intrachain hydrogen bonds

hydrophobic contacts

extended loosely collapsed tightly collapsed

1.91 (0.05) 1.13 (0.02) 1.04 (0.01)

1.3 (1.0) 4.2 (1.1) 5.2 (0.8)

20.9 (2.6) 31.9 (5.4) 36.9 (3.1)

The values in parentheses are the standard deviation.

in Table 1. An extended-chain conformation is obtained from the Gaussview software,34 and the 1 μs simulations are run at each temperature between 270 and 320 K. A loosely collapsedchain conformation is obtained from the simulation of the extended chain at 300 K, for which the chain collapses significantly at times, and the 1 μs simulations are run with the loosely collapsed conformation at each temperature between 270 and 300 K. Through the hydrogen-bonding analysis of the simulation trajectory of a loosely collapsed-chain conformation at 300 K, as will be shown in the next section, it is found that B

DOI: 10.1021/acs.jpcb.6b09165 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B the collapsed conformations can be distinguished into two different ones based on the stability of the intrachain hydrogen bonding, namely, one with loosely collapsed structures and another with tightly collapsed structures. Taken from the simulation of the loosely collapsed chain at 300 K, a tightly collapsed-chain conformation is used as an input for the 1 μs simulations performed at temperatures between 270 and 300 K. We refer to each set of these simulations as extended-chain, loosely collapsed-chain, and tightly collapsed-chain simulations hereafter. It is noted that the loosely collapsed-chain and tightly collapsed-chain simulations are run at temperatures from 270 and 300 K, whereas the extended-chain simulations are run at 310 and 320 K as well as at temperatures between 270 and 300 K, to confirm the collapse transition of the extended-chain conformation at higher temperatures. The simulation temperatures are chosen below or near the experimentally estimated LCST of 305 K such that the observation of reversible collapse−swelling transitions becomes more probable. We employ the OPLS-AA force-field25,35,36 for PNIPAM chains and the TIP4P/2005 force-field37 for water molecules. The Lennard-Jones (LJ) parameters for a pair of different types of atoms are determined using a geometric average of the LJ parameters of both atoms, including those between PNIPAM atoms and a water oxygen. Several force-fields have been used in previous MD studies, and we use the OPLS-AA force-field for PNIPAM chains because the OPLS force-field was shown to predict the temperature-dependent change in average chain sizes qualitatively well in combinations with several water models.25,28,31 The TIP4P/2005 force-field is employed for water because it describes the water properties relatively well at low temperatures,38 compared to other water models that have been employed in combination with a PNIPAM chain. The TIP4P/2005 water model, for instance, predicts the melting point at 249 K,38,39 whereas the SPC/E and the TIP4P models estimate the melting temperature at 215 and 232 K, respectively, implying that the hydrogen-bonding interaction is better described with the TIP4P/2005 model at low temperatures. All NPT simulations are performed with GROMACS v 5.1.2.40 All simulations are performed using a time step of 2 fs, and the temperature is set using a V-rescale method41 with a time constant of 0.1 ps. The pressure of the system is set at 1 bar using a Parrinello−Rahman algorithm42,43 with a time constant of 2.0 ps. A cutoff distance of 1 nm is used for LJ interactions, and the long-range electrostatic interactions are calculated using the Particle-Mesh Ewald method44 with a real space cutoff distance of 1 nm. A periodic boundary condition is applied with the simulation box size of roughly 7 × 7 × 7 nm3, with a small deviation due to pressure coupling. The system contains a single 30-mer PNIPAM chain with 11 262 water molecules. Due to the use of the periodic boundary condition, one end of a PNIPAM chain in the simulation box may interact with the other end of the chain image assumed on its neighboring periodic box if the box size is not large enough. The possibility of the artificial interaction of a PNIPAM chain with its periodic images is tested by calculating the closest distance of one chain end from the other end of all the periodic images as well as from that of the chain itself, as explained in Figure S1a. It is shown that the closest distance in Figure S1b is rarely less than 1.5 nm unless the chain folds such that the intrachain ends get closer (verified by comparison with the change in the radius of gyration shown in Figure 2 that is presented in the next section). Therefore, it is concluded that

Figure 2. Changes in the radius of gyration with time at temperatures between 270 and 320 K. The three lines at each temperature between 270 and 300 K are obtained from three sets of 1 μs simulations starting from extended (in orange), loosely collapsed (in blue), and tightly collapsed conformations (in gray), respectively. Simulations are performed at 310 and 320 K only for the sets with an extended conformation, and the corresponding data are presented at 310 and 320 K.

the direct contact between the original chain and its periodic images is highly unlikely during the simulation trajectories obtained in this work.

■

RESULTS AND DISCUSSION Collapse−Swelling Transition in Terms of the Radius of Gyration. The collapse−swelling transitions of a single PNIPAM chain are characterized by the change in the radius of gyration of the PNIPAM chain, as shown in Figure 2. In the figure, the change in the radius of gyration is presented for three sets of simulations that started with extended (in orange), loosely collapsed (in blue), and tightly collapsed conformations (in gray), respectively, at temperatures between 270 and 300 K. Simulations are performed at 310 and 320 K only for the system with an initially extended conformation. The distribution of the radius of gyration of the PNIPAM chain is presented in Figures S2, S3, and S7 for the extended, loosely collapsed, and tightly collapsed conformations, respectively. For the system with an extended conformation (orange curves in Figure 2), it is shown that the radius of gyration fluctuates most significantly at 300 K, indicating reversible collapse−swelling transitions. The chain collapses after the initial 150 ns, remains mostly collapsed for the next 350 ns with transient extensions, and finally swells back again at 500 ns with repeated transient collapses afterward. It is noticed that the swelling of the collapsed chain at 500 ns would have not been observed if the simulation duration were less than 500 ns. At temperatures below 300 K, the chain remains extended for most of the simulation duration. In contrast, the chain remains mostly collapsed at 310 and 320 K after the initial collapse at around 100 ns. To determine the overall transition between C

DOI: 10.1021/acs.jpcb.6b09165 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

limited to a few times, as shown in Figure 2. Therefore, it is emphasized that the conformational sampling is not sufficient, and the transition temperature estimated above does not represent the thermodynamic transition temperature for the collapse and swelling of a PNIPAM chain. However, 1 μs-long simulation trajectories in this work allow for the observation of reversible changes with regard to the chain collapse and swelling, although limited to a few times. The transient collapse−swelling transitions are most clearly observed at 300 K for the extended-chain simulation and at 270 and 280 K for the loosely collapsed-chain simulations. Such reversible collapse−swelling transitions were rarely reported in previous MD studies,23 probably due to the limited simulation time or due to the trapping of collapsed conformations in a local energy minimum. In the previous MD studies, the average value and the distribution of the radius of gyration were often presented at various temperatures and used to infer the LCST from the MD simulations of a single PNIPAM chain, eliminating the dynamic details regarding the collapse−swelling transitions. Instead of calculating the average properties due to the insufficiency of conformational sampling, however, we analyze the structural changes occurring during the reversible collapse−swelling transitions over the entire trajectory, presenting more dynamic information hidden in the average values of the radius of gyration. Dynamical Changes of Chain Conformations during the Collapse−Swelling Transitions. The reversible collapse−swelling transition occurs most clearly in the loosely collapsed-chain simulations. Therefore, the structural changes are analyzed for the loosely collapsed-chain simulations in terms of the intrachain distances, hydrophobic contacts, and hydrogen bonds. The same analyses are performed for the extended-chain and tightly collapsed-chain simulations, and the results are presented in the Supporting Information (Figures S4−S6 for the extended-chain simulations and Figures S8−S10 for the tightly collapsed-chain simulations). Intrachain Spatial Distances. The time-averaged intrachain distance and segment−segment contact probability were calculated in previous MD studies to reveal structural differences between the chain conformations at temperatures below and above the LCST.30,32 Here, we present instantaneous intrachain distances to show the change in the chain conformations during the course of reversible collapse−swelling transitions, as shown in Figure 3. Figure 3a shows the intrachain spatial distances of three initial chain conformations as a function of the sequential difference between a pair of chain monomers, |i − j|. The spatial distance at |i − j| = 29 indicates the chain end-to-end distance, whereas the distances at |i − j| < 29 are the average distance between all pairs of monomers that are separated by |i − j| bonds. For instance, there are three different pairs of monomers, (i, j) = (1, 28), (2, 29), and (3, 30), for which |i − j| = 27, and the intrachain distance at |i − j| = 27 is the average distance for the three pairs. As expected, the extendedchain conformation has an intrachain spatial distance that increases with an increase in |i − j|. On the other hand, the intrachain distances for the collapsed conformations increase initially but decrease at higher values of |i − j| due to the chain collapse. Figure 3b shows the instantaneous intrachain distances as a function of time in the simulation trajectories starting with a loosely collapsed-chain conformation. Only the results of a subset of the |i − j| values, that is, |i − j| = 13, 17, 21, 25, and 29,

collapsed conformations at lower temperatures and extended conformations at higher temperatures, we calculate the probability that the radius of gyration is less than a critical distance of 1.1 nm. The critical distance is chosen based on the distribution of the radius of gyration from the tightly collapsedchain simulations in which the radius of gyration is less than 1.1 nm (in Figure S7). The probability is provided in Table 2. Table 2. Probability That the Radius of Gyration is Less than a Cutoff Distance rcut = 1.1 nm for Different Initial Conformations and Temperatures initial chain conformation extended

loosely collapsed

tightly collapsed

temp. (K)

prob. (Rg < rcut; rcut = 1.1 nm)

270 280 290 300 310 320 270 280 290 300 270 280 290 300

0.00 0.01 0.27 0.30 0.71 0.89 0.35 0.26 0.59 0.93 1.00 1.00 1.00 1.00

Large probability implies that the collapsed conformation with the radius of gyration less than 1.1 nm is more abundant, whereas small probability implies that the extended conformation is more abundant. Assuming that the chain is in the collapsed state for probability ≥0.5, the overall transition temperature for the extended-chain simulations is then estimated to be between 300 and 310 K. For the system with a loosely collapsed conformation shown by blue curves in Figure 2, the overall transition temperature is estimated to be between 280 and 290 K based on the probability ≥0.5 in Table 2. At 270 and 280 K, the swelling of the chain from a collapsed conformation is clearly observed with subsequent collapse−swelling transitions. At 290 K, however, the chain is slightly extended in the first 100 ns, collapses again for the next 400 ns, and repeats the swelling and collapse afterward. Although repeated collapse−swelling transitions are observed in the last 500 ns, the chain remains collapsed more probably, as shown in Table 2. At 300 K, the chain remains collapsed except in the first 50 ns. The simulations of the tightly collapsed conformation at temperatures between 270 and 300 K, as shown by gray curves in Figure 2, do not show any deviation of the radius of gyration from the initial value, suggesting that the swelling transition of the tightly collapsed conformation does not occur during the simulation duration of 1 μs in the temperature range considered in this work. Therefore, it is concluded that the overall transition temperature is below 270 K. From the comparison of the three sets of simulations, we find that the overall transition temperature for the collapse−swelling transition determined by the unbiased 1 μs simulations depends on the initial chain conformation. Unfortunately, this implies that the simulations do not sample PNIPAM conformations enough to draw meaningful conclusions on the equilibrium properties. In addition, the number of reversible structural transitions between the collapsed and swollen conformations is D

DOI: 10.1021/acs.jpcb.6b09165 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

extended for a significant amount of time, and then collapsed again at the end. The same interpretation is possible for 300 K, where the chain is extended for a short time at the beginning, collapsed rapidly at 50 ns, extended again at 350 ns, and then collapsed at 550 ns for the rest of the simulation. This representation of the distinction between collapsed and extended conformations may provide a slightly more adequate depiction of the chain collapse−swelling transitions than that based on the mere change in the radius of gyration. For the loosely collapsed-chain simulation at 300 K, the chain is mostly collapsed based on the radius of gyration except in the first 50 ns and for a very short time at around 550 ns, as shown in Figure 2. Based on the intrachain distances in Figure 3b, however, the chain is extended for about 200 ns between 350 ns and 550 ns in the middle of the trajectory as well as in the first 50 ns, as described above. In this period between 350 and 550 ns, the radius of gyration ranges between 0.93 and 1.20 (implying a collapsed conformation) before it suddenly increases to ∼1.5 at 550 ns. In contrast, the curves of the intrachain distances have an order of |i − j| = 29 > 25 > 21 > 17 ≈ 13 between 350 and 550 ns, indicating that the chain is extended in this time range. Hydrophobic Contacts. The LCST behavior of PNIPAM with swelling at low temperature and collapse at high temperature has been interpreted as a consequence of the competitive roles of hydrogen-bonding and hydrophobic interactions.15,19−21,45 The chain conformations are determined both by the hydrophilic and the hydrophobic interactions at all temperature ranges, and it was shown that an earlier notion that the chain is hydrophilic at low temperatures and becomes hydrophobic at high temperatures is incorrect.45 At low temperatures, hydrophilic interaction through hydrogen bonding between amide groups of the PNIPAM chain and water molecules is dominant and hence the chain is hydrated by water molecules, keeping it swollen in water. As the temperature increases, the hydrophobic interaction becomes more dominant, resulting in the dehydration and collapse of the chain, inhibiting the hydrogen bond formation between chain segments and water molecules. In this work, we analyze the changes in hydrophobic interactions and hydrogen-bonding interactions during the collapse−swelling transitions of the loosely collapsed-chain simulations. The analysis on the hydrophobic interactions is presented in this subsection and those on the hydrogen bonding are in the following subsection. The change in the hydrophobic interactions is quantified in terms of the contact numbers between isopropyl CH3 groups in different monomers of a PNIPAM chain. Although the backbone of the chain is also considered hydrophobic, it is not counted here for the hydrophobic contact for simplicity. For this purpose, a quantitative criterion is required for the hydrophobic contacts between isopropyl CH3 groups, which is provided by analyzing the radial distribution functions (RDFs) between isopropyl CH3 groups in different monomers at each temperature. The RDFs between CH3 groups in Figure 4a show a temperature dependence with the overall larger values at 290 and 300 K and with smaller values at 270 and 280 K. This temperature dependence is in agreement with our estimation of the transition temperature based on the radius of gyration. At higher temperatures, the chain is mostly collapsed and thus, contact between isopropyl CH3 groups is more frequent than that at lower temperatures where the chain is largely swollen. Here, we attribute the first peak at around 0.40 nm to the hydrophobic contacts. Therefore, the distance

Figure 3. (a) Intrachain distances in the three initial chain conformations. Distances are presented as a function of the number of in-between bonds separating a pair of monomers i and j, represented as |i − j|. (b) Instantaneous intrachain distances shown for entire trajectories of loosely collapsed-chain simulations. Data are presented only for a subset of the intrachain distances with |i − j| = 13, 17, 21, 25, and 29 in black, green, red, orange, and purple, respectively.

are presented for visual clarity. The maximum distance for the collapsed chains in Figure 3a is found at around |i − j| = 15. Therefore, the curves for |i − j| = 13 and 17 (in black and green in Figure 3b, respectively) either have the largest values if the chain is collapsed or have the smallest values if the chain is extended. The relative values of the curves with |i − j| = 13 and 17 to those with |i − j| > 17 can be used as an indicator of the chain collapse and swelling. The instantaneous intrachain distances as a function of time in Figure 3b clearly reveals the reversible collapse−swelling transitions. The collapsed conformations are characterized by the distances in the order |i − j| = 13 > 17 > 21 > 25 > 29, and the extended conformations are characterized by the distances in the opposite order. At 280 K, for instance, the data for |i − j| = 13 (black curve) have the highest values at the beginning and end of the simulation, whereas they have the smallest values during the rest of the simulation trajectory. This implies that the chain is collapsed at the beginning of the simulation, E

DOI: 10.1021/acs.jpcb.6b09165 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

number of NIPAM monomers, which is 30 in this work. This indicates that the hydrophobic interaction is not negligible even at low temperatures. The non-negligible values of the hydrophobic contacts in the extended conformations are expected to originate from the local interactions, whereas larger contact numbers in the collapsed conformations are from the global interactions that collapse the chain. To confirm the local and global hydrophobic interactions in the extended and collapsed conformations, respectively, we perform cluster analysis of the monomers involved in the hydrophobic contact formation. Each cluster is defined as a group of monomers mutually connected through hydrophobic contacts between isopropyl CH3 groups. For each isopropyl C atom, we calculate the distance with all other isopropyl C atoms on different monomers. If the distance is less than a critical distance of 0.49 nm, the monomers bearing these isopropyl C atoms are considered to belong to the same cluster. The cluster analysis on the hydrophobic contacts is shown in Figure 4c for 280 and 300 K, representing the most extended and collapsed conformations, respectively. In the figure, the monomers involved in hydrophobic contact formation and mutually connected to each other are depicted in the same color, but those not mutually connected are colored differently, as shown in the color bar on the right-hand side of the figure. In the figure for 280 K, different colors appear over the course of the simulation, implying that there are several clusters that form hydrophobic contacts. In particular, at specific times during the period between 200 and 800 ns, where the chain is mostly extended, the yellow-to-red colors are frequently shown together with blue and green, meaning that there can be up to 10 separate clusters forming hydrophobic contacts. In such cases, the same colors are shown for the monomers close to each other, and this proves that the hydrophobic interactions in the extended-chain conformations are mostly local. On the other hand, at 300 K where the chain is largely collapsed, only the cluster in blue is dominant, with much less frequent appearance of green and yellow. This means that the monomers with hydrophobic contacts are mostly connected to each other and hydrophobic interaction is global throughout the chain in the collapsed conformations. Hydrogen Bonding. In addition to the hydrophobic interactions, the hydrogen-bonding interactions between water molecules and chain segments as well as those between chain segments are important for determining the chain conformation. The hydrogen-bonding analysis is performed based on the definition by the donor−acceptor distance less than 0.35 nm and the hydrogen−donor−acceptor angle less than 30°. The number of hydrogen bonds between water molecules and chain segments is presented in Figure 5a and that between intrachain segments is in Figure 5b. The numbers of hydrogen bonds given in the figures are the average values over the intervals of each 1 ns. The changes in the number of each type of hydrogen bonds are also in qualitative agreement with the chain conformational change. When the chain conformation is extended, the number of chain−water hydrogen bonds reaches up to 72. When the chain is collapsed, the number decreases to 52. The number of intrachain hydrogen bonds shows the opposite trend, with higher numbers for collapsed conformation and lower numbers for extended conformation. However, the number of intrachain hydrogen bonds is fairly small. It reaches, at maximum, ∼5 when the chain is mostly collapsed. Therefore, it can be concluded that the formation of intrachain hydrogen bonds is

Figure 4. Analysis of hydrophobic contacts in the loosely collapsedchain simulations. (a) Changes in the RDF g(r) between isopropyl CH3 groups in different monomers of a chain. g(r) decreases with the decrease in temperature in the order 300, 290, and 270 ≈ 280 K. (b) Changes in the number of hydrophobic contacts are defined as a pair of isopropyl CH3 groups within a distance of 0.49 nm. This is the distance corresponding to the first minimum of g(r), as indicated by the dashed line in (a). (c) Cluster analysis of the monomers involved in the hydrophobic contact formation. The cluster analysis data are presented at every 10 ns. The color bar on the right indicates different clusters of monomers mutually connected to each other.

of the first minimum at 0.49 nm is used as a criterion to calculate the time evolution of the number of hydrophobic contacts, as shown in Figure 4b,c. The contact number between isopropyl CH3 in Figure 4b shows the evolution of the hydrophobic contacts during the collapse−swelling transitions of a single PNIPAM chain. The change in the hydrophobic contact numbers is consistent with the change in the chain conformation. For instance, at 280 K, the chain is swollen at the beginning and end of the simulation, as shown in Figures 2 and 3, and collapsed in the rest of the simulation. The contact number shown in Figure 4b is larger at the beginning and end of the trajectory but smaller elsewhere. At 300 K, the chain is swollen before 50 ns and between 350 and 550 ns and is collapsed in the rest of the trajectory based on the intrachain distances of Figure 3. The contact number in Figure 4b also shows the same behavior: smaller values at the beginning and larger values afterward except a slight decrease in the contact numbers between 350 and 550 ns. Here, we emphasize that the hydrophobic contact number is large even in the swollen conformations at low temperatures. In the middle of the trajectory at 280 K, where the chain is mostly extended, the contact number is estimated to be 15−30. Although this contact number is less than 30−45 for the most collapsed conformations, it is still significant considering the F

DOI: 10.1021/acs.jpcb.6b09165 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

trajectory except in the first 50 ns. However, the number of intrachain hydrogen bonds changes significantly earlier and becomes stable after 700 ns. We conclude that the collapsedchain conformations with similar radii of gyration can be distinguished into two collapsed conformations, one with highly stable hydrogen bonding and the other with less stable hydrogen bonding. Throughout the article, we distinguished the collapsed conformations into tightly collapsed and loosely collapsed conformations, based on such differences in the stability of the hydrogen bonds, as was briefly explained in the “Simulation Methods” section. The stability of intrachain hydrogen bonds was discussed above, solely in terms of very small variation in the number of intrachain hydrogen bonds, as shown during the last 300 ns at 300 K in Figure 5b. The stability of the intrachain hydrogen bonding is also demonstrated by the cluster analysis in Figure 5c. In the figure, the monomers connected by intrachain hydrogen bonding are represented with the same color. In the cluster analysis data for 300 K, multiple bands representing the connected monomers are maintained during the last 300 ns of the simulation, implying that intrachain hydrogen bonds connecting the same pairs of monomers persist for a significant period of time. The formation of stable intrachain hydrogen bonds is more probable at higher temperatures, as shown in Figure 5. The same temperature dependence is observed in the extendedchain simulations (Figure S6), in which the formation of stable intrachain hydrogen bonds is observed only at 310 and 320 K. Interestingly, however, the stable intrachain hydrogen bonds, once formed, persist longer at lower temperatures. The intrachain hydrogen bonds in the tightly collapsed-chain simulations (Figure S10) remain stable for several hundreds of nanoseconds, but become less stable as the temperature increases from 270 to 300 K, especially, in the final 300 ns. This is due to weakening of hydrogen-bonding interactions at high temperatures. At the same time, however, the number of hydrophobic contacts increases at higher temperatures, as shown in the last 300 ns of Figure S9, whereas the chain radius of gyration remains about the same in Figure 2. Therefore, it can be conceived that the depressed intrachain hydrogen bonding at high temperatures is compensated by the increase in the hydrophobic interactions to keep the chain collapsed during conformational reorganization, emphasizing the role of hydrophobic interactions in chain collapsing at high temperatures. Finally, the formation of stable intrachain hydrogen bonds may have prevented the tightly collapsed chain from swelling within the simulation duration of 1 μs at temperatures between 270 and 300 K by trapping the collapsed conformations in a local energy minimum. It was reported in previous experimental studies that the coil-to-globule (collapse) transition during heating occurs at higher temperatures than that at which the globule-to-coil (swelling) transition occurs during cooling.13,14,19 This thermal hysteresis was interpreted as a consequence of the persistence of intrachain structures formed in the globule state,13 which is important in the globuleto-coil (swelling) transition during cooling but not in the coilto-globule (collapse) transition during heating. Therefore, it is not absurd to presume that the formation of stable intrachain hydrogen bonds at high temperatures maintains the collapsed state of a PNIPAM chain for a long period of time even when the temperature is lowered during cooling.

Figure 5. Analysis of hydrogen bonding in the loosely collapsed-chain simulations. Changes in the number of hydrogen bonds (a) between chain amide groups and water molecules and (b) between intrachain amide groups (CO···H−N). The numbers given in (a) and (b) are the average values over the intervals of each 1 ns. The dashed lines at (a) 55.3 and (b) 5.2 are guides to the eye for the tightly collapsed conformations, averaged from 700 to 1000 ns during the simulation of a precollapsed chain at 300 K. (c) Cluster analysis of the monomers involved in the intrachain hydrogen bonding. The data are presented at every 10 ns. The color bar on the right indicates different clusters of monomers mutually connected to each other.

not a major driving force for the chain collapse. As in the case of the hydrophobic contacts, we perform cluster analysis to identify the monomers connected through intrachain hydrogen bonds, as shown in Figure 5c. It is determined that the intrachain hydrogen bonding is truly scarce, with only a few monomers mutually connected to each other via hydrogen bonding. It is obvious that there are more clusters of the monomers that are mutually connected through hydrogen bonding in the collapsed conformation at 300 K than those in the overall extended conformation at 280 K. In addition, for collapsed conformations at 300 K, the intrachain hydrogen bonds between distant monomers (presented in the same color) are shown in Figure 5c. Closer observation of the change in the hydrogen bonds reveals that there are two regimes in the simulation trajectory that have different variations in the number of hydrogen bonds, indicating different levels of stability for the hydrogen-bonding interactions. The difference in stability is more pronounced in the case of the intrachain hydrogen bonds at 300 K, as shown in Figure 5b. The number of intrachain hydrogen bonds is almost constant near 5.2 after 700 ns, whereas it fluctuates significantly before 700 ns. Looking back at Figure 2, the radius of gyration at 300 K remains almost constant during the entire simulation G

DOI: 10.1021/acs.jpcb.6b09165 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

■

CONCLUSIONS Simulation results in this work reveal that the overall transition temperatures determined by the unbiased 1 μs simulations of a single PNIPAM chain depend on the initial chain conformation. The inconsistency of the estimated transition temperatures depending on the initial conformation implies a deficiency in conformation sampling, and it is not convincing to discuss the equilibrium properties of the chain. Therefore, this work suggests that caution should be exercised when conformational transition of a polymer chain that involves the formation and breakage of strong intersegment interactions is studied by conventional MD simulations. However, the simulation duration of 1 μs was long enough for each trajectory to show repeated collapse−swelling transitions of a PNIPAM chain, which has not been reported previously. In this work, we analyzed the structural changes that occur during the reversible collapse−swelling transitions of a single PNIPAM chain. Intrachain spatial distances over the entire trajectory show that the chain conformation can still be open and extended even when the radius of gyration is comparable to that of tightly collapsed conformations. The analysis of hydrophobic contacts over collapse−swelling transitions at various temperatures showed that the number of hydrophobic contacts varies significantly depending on the chain conformation, with large values for collapsed conformations and smaller values for extended conformations. Nevertheless, the number of hydrophobic contacts is non-negligible even when the chain is extended. In such cases, hydrophobic contacts are formed locally among monomers in close proximity. For the collapsed conformations, hydrophobic interaction is global over the entire chain. The number of hydrogen-bonding interactions changes according to the chain conformation. The number of hydrogen bonds between the chain and water decreases with chain collapse but increases with chain swelling. The number of hydrogen bonds between intrachain monomers has opposite dependence on the chain conformation. The most notable finding in the analysis of hydrogen bonding is that the collapsed conformations may be distinguished into two: one with very stable intrachain hydrogen bonding and another with less stable hydrogen bonding. Throughout this article, we refer to these two types of conformations as tightly collapsed and loosely collapsed conformations, respectively. These conformations have similar radii of gyration but distinct dynamic behaviors. It is concluded that intrachain hydrogen bonding may play a significant role in keeping the collapsed structure tight once the chain collapses, although it may not be a major driving force for the chain collapse considering the small number of intrachain hydrogen bonds. Finally, although an isotactic PNIPAM chain was employed in this work, the tacticity of the chain may require more careful consideration. In the literature, it was shown that the LCST depends on the tacticity of the chain.32,46,47 An increase in the isotacticity of the polymer decreases the solubility of the polymer in water, and the polymer with very large isotacticity does not dissolve in water even at low temperatures without a phase transition across the LCST.46,47 However, our model system with a single, 30-mer PNIPAM chain is clearly different from the experimental PNIPAM solution, with significantly higher degrees of polymerization (>300); 46 thus, the observation of the LCST behavior for the isotactic chain in this work may still be plausible due to much shorter polymer

length, as has been shown by a significant increase in the transition temperature with decrease in the polymer length, particularly, for polymers with low molecular weights (≤∼100).11 In addition, because the isotactic PNIPAM is known to have the greatest tendency for chain collapse, the use of the isotactic conformation in this work can be rationalized to accelerate the collapse transition of an extended conformation in silico.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b09165. Complete data of the extended-chain and tightly collapsed-chain simulations and additional data for the loosely collapsed-chain simulations (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jun Soo Kim: 0000-0001-8113-0605 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the National Research Foundation of Korea (NRF) under Grant Nos. NRF-20110024621 and NRF-2015R1D1A1A01059005.

■

REFERENCES

(1) Heskins, M.; Guillet, J. E. Solution Properties of Poly(Nisopropylacrylamide). J. Macromol. Sci., Chem. 1968, 2, 1441−1455. (2) Bae, Y. H.; Okano, T.; Kim, S. W. Temperature Dependence of Swelling of Crosslinked Poly(N,N′-alkyl substituted acrylamides) in Water. J. Polym. Sci., Part B: Polym. Phys. 1990, 28, 923−936. (3) Schild, H. G. Poly(N-isopropylacrylamide): Experiment, Theory and Application. Prog. Polym. Sci. 1992, 17, 163−249. (4) Yoshida, R.; Uchida, K.; Kaneko, Y.; Sakai, K.; Kikuchi, A.; Sakurai, Y.; Okano, T. Comb-type Grafted Hydrogels with Rapid Deswelling Response to Temperature Changes. Nature 1995, 374, 240−242. (5) Okano, T.; Yamada, N.; Okuhara, M.; Sakai, H.; Sakurai, Y. Mechanism of Cell Detachment from Temperature-Modulated, Hydrophilic-Hydrophobic Polymer Surfaces. Biomaterials 1995, 16, 297−303. (6) Hoffman, A. S.; Stayton, P. S.; Bulmus, V.; Chen, G.; Chen, J.; Cheung, C.; Chilkoti, A.; Ding, Z.; Dong, L.; Fong, R.; et al. Really Smart Bioconjugates of Smart Polymers and Receptor Proteins. J. Biomed. Mater. Res. 2000, 52, 577−586. (7) Jeong, B.; Kim, S. W.; Bae, Y. H. Thermosensitive Sol-Gel Reversible Hydrogels. Adv. Drug Delivery Rev. 2002, 54, 37−51. (8) Schmaljohann, D. Thermo- and PH-responsive Polymers in Drug Delivery. Adv. Drug Delivery Rev. 2006, 58, 1655−1670. (9) Cobo, I.; Li, M.; Sumerlin, B. S.; Perrier, S. Smart Hybrid Materials by Conjugation of Responsive Polymers to Biomacromolecules. Nat. Mater. 2015, 14, 143−159. (10) Fujishige, S.; Kubota, K.; Ando, I. Phase Transition of Aqueous Solutions of Poly(N-isopropylacrylamide) and Poly(N-isopropylmethacrylamide). J. Phys. Chem. 1989, 93, 3311−3313. (11) Pamies, R.; Zhu, K.; Kjøniksen, A.-L.; Nyström, B. Thermal Response of Low Molecular Weight Poly-(N-isopropylacrylamide)Polymers in Aqueous Solution. Polym. Bull. 2009, 62, 487. H

DOI: 10.1021/acs.jpcb.6b09165 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(32) Chiessi, E.; Paradossi, G. Influence of Tacticity on Hydrophobicity of Poly(N-isopropylacrylamide): A Single Chain Molecular Dynamics Simulation Study. J. Phys. Chem. B 2016, 120, 3765−3776. (33) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (34) Dennington, R.; Keith, T.; Millam, J. GaussView, version 5; Semichem Inc., 2009. (35) Jorgensen, W. L.; Tirado-Rives, J. The OPLS potential Functions for Proteins. Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (36) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (37) Abascal, J. L.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, No. 234505. (38) Vega, C.; Abascal, J. L. F. Physics and Chemistry of Water and Ice. Phys. Chem. Chem. Phys. 2011, 13, 19660−19662. (39) Fernández, R. G.; Abascal, J. L. F.; Vega, C. The Melting Point of Ice l-h for Common Water Models Calculated from Direct Coexistence of the Solid-Liquid Interface. J. Chem. Phys. 2006, 124, No. 144506. (40) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations Through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1−2, 19−25. (41) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, No. 014101. (42) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (43) Nosé, S.; Klein, M. L. Constant Pressure Molecular-Dynamics for Molecular-Systems. Mol. Phys. 1983, 50, 1055−1076. (44) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577. (45) Pelton, R. Poly(N-isopropylacrylamide) (PNIPAM) Is Never Hydrophobic. J. Colloid Interface Sci. 2010, 348, 673−674. (46) Ray, B.; Okamoto, Y.; Kamigaito, M.; Sawamoto, M.; Seno, K.; Kanaoka, S.; Aoshima, S. Effect of Tacticity of Poly(N-isopropylacrylamide) on the Phase Separation Temperature of Its Aqueous Solutions. Polym. J. 2005, 37, 234−237. (47) Biswas, C. S.; Vishwakarma, N. K.; Patel, V. K.; Mishra, A. K.; Saha, S.; Ray, B. Synthesis and Study of the Properties of Stereocontrolled Poly(N-isopropylacrylamide) Gel and its Linear Homopolymers Prepared in the Presence of a Y(OTf)3 Lewis Acid: Effect of the Composition of Methanol-Water Mixtures as Synthesis Media. Langmuir 2012, 28, 7014−7022.

(12) Kubota, K.; Fujishige, S.; Ando, I. Single-Chain Transition of Poly(N-isopropylacrylamide) in Water. J. Phys. Chem. 1990, 94, 5154− 5158. (13) Wu, C.; Wang, X. Globule-to-Coil Transition of a Single Homopolymer Chain in Solution. Phys. Rev. Lett. 1998, 80, 4092− 4094. (14) Lu, Y.; Zhou, K.; Ding, Y.; Zhang, G.; Wu, C. Origin of Hysteresis Observed in Association and Dissociation of Polymer Chains in Water. Phys. Chem. Chem. Phys. 2010, 12, 3188−3194. (15) Otake, K.; Inomata, H.; Konno, M.; Saito, S. Thermal Analysis of the Volume Phase Transition with N-Isopropylacrylamide Gels. Macromolecules 1990, 23, 283−289. (16) Schild, H.; Tirrell, D. Microcalorimetric Detection of Lower Critical Solution Temperatures in Aqueous Polymer Solutions. J. Phys. Chem. 1990, 94, 4352−4356. (17) Ohta, H.; Ando, I.; Fujishige, S.; Kubota, K. Molecular Motion and 1H NMR Relaxation of Aqueous Poly(N-isopropylacrylamide) Solution Under High Pressure. J. Polym. Sci., Polym. Phys. Ed. 1991, 29, 963−968. (18) Tokuhiro, T.; Amiya, T.; Mamada, A.; Tanaka, T. NMR Study of Poly(N-isopropylacrylamide) Gels near Phase Transition. Macromolecules 1991, 24, 2936−2943. (19) Maeda, Y.; Higuchi, T.; Ikeda, I. Change in Hydration State during the Coil-Globule Transition of Aqueous Solutions of Poly(Nisopropylacrylamide) as Evidenced by FTIR Spectroscopy. Langmuir 2000, 16, 7503−7509. (20) Meersman, F.; Wang, J.; Wu, Y.; Heremans, K. Pressure Effect on the Hydration Properties of Poly(N-isopropylacrylamide) in Aqueous Solution Studied by FTIR Spectroscopy. Macromolecules 2005, 38, 8923−8928. (21) Winnik, F. Fluorescence Studies of Aqueous Solutions of Poly(N-isopropylacrylamide) below and above Their LCST. Macromolecules 1990, 23, 233−242. (22) Longhi, G.; Lebon, F.; Abbate, S.; Fornili, S. L. Molecular Dynamics Simulation of a Model Oligomer for Poly(N-isopropylamide) in Water. Chem. Phys. Lett. 2004, 386, 123−127. (23) Gangemi, F.; Longhi, G.; Abbate, S.; Lebon, F.; Cordone, R.; Ghilardi, G. P.; Fornili, S. L. Molecular Dynamics Simulation of Aqueous Solutions of 26-unit Segments of p(NIPAAm) and of p(NIPAAm) “Doped” with Amino Acid Based Comonomers. J. Phys. Chem. B 2008, 112, 11896−11906. (24) Deshmukh, S.; Mooney, D. A.; McDermott, T.; Kulkarni, S.; MacElroy, J. M. D. Molecular Modeling of Thermo-responsive Hydrogels: Observation of Lower Critical Solution Temperature. Soft Matter 2009, 5, 1514−1521. (25) Walter, J.; Ermatchkov, V.; Vrabec, J.; Hasse, H. Molecular Dynamics and Experimental Study of Conformation Change of Poly(N-isopropylacrylamide) Hydrogels in Water. Fluid Phase Equilib. 2010, 296, 164−172. (26) Du, H.; Wickramasinghe, R.; Qian, X. Effects of Salt on the Lower Critical Solution Temperature of Poly(N-isopropylacrylamide). J. Phys. Chem. B 2010, 114, 16594−16604. (27) Pang, J.; Yang, H.; Ma, J.; Cheng, R. Solvation Behaviors of Nisopropylacrylamide in Water/Methanol Mixtures Revealed by Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 8652− 8658. (28) Tucker, A. K.; Stevens, M. J. Study of the Polymer Length Dependence of the Single Chain Transition Temperature in Syndiotactic Poly(N-isopropylacrylamide) Oligomers in Water. Macromolecules 2012, 45, 6697−6703. (29) Alaghemandi, M.; Spohr, E. Molecular Dynamics Investigation of the Thermo-Responsive Polymer Poly(N-isopropylacrylamide). Macromol. Theory Simul. 2012, 21, 106−112. (30) Abbott, L. J.; Tucker, A. K.; Stevens, M. J. Single Chain Structure of a Poly(N-isopropylacrylamide) Surfactant in Water. J. Phys. Chem. B 2015, 119, 3837−3845. (31) Boţan, V.; Ustach, V.; Faller, R.; Leonhard, K. Direct Phase Equilibrium Simulations of NIPAM Oligomers in Water. J. Phys. Chem. B 2016, 120, 3434−3440. I

DOI: 10.1021/acs.jpcb.6b09165 J. Phys. Chem. B XXXX, XXX, XXX−XXX