Atomistic Molecular Dynamics Simulations of the Lower Critical

Subscriber access provided by Technical University of Munich University Library

B: Fluid Interfaces, Colloids, Polymers, Soft Matter, Surfactants, and Glassy Materials

Atomistic Molecular Dynamics Simulations of The Lower Critical Solution Temperature Transition of Poly(N-Vinylcaprolactam) in Aqueous Solutions Xiaoquan Sun, and Xianghong Qian J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b01711 • Publication Date (Web): 24 May 2019 Downloaded from http://pubs.acs.org on May 27, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Atomistic Molecular Dynamics Simulations of The Lower Critical Solution Temperature Transition of Poly(N-vinylcaprolactam) in Aqueous Solutions Xiaoquan Sun and Xianghong Qian* Department of Biomedical Engineering, University of Arkansas, Fayetteville, AR 72701

Abstract Poly (N-vinylcaprolactam) (PVCL) is a thermo-responsive polymer which exhibits a lower critical solution temperature (LCST) in aqueous solution. The LCST of this hydrophilic-to-hydrophobic transition is found to be strongly dependent on the salt type and salt concentration as well as on the molecular weight and concentration of polymer. Here atomistic molecular dynamics (MD) simulations have been successfully conducted for the first time to investigate the LCST transition of a 100 degree of polymerization (DP) PVCL chain in water, 1 M NaCl, 3.5 M NaCl and 0.5 M CaCl 2 solutions. Our results show that steric hinderance resulting from the bulky 7-member ring on the PVCL chain plays a critical role in conformational transition. Moreover, the degrees of hydration and dehydration below or above the transition temperature are highly dependent on the specific solution condition and temperature. Water molecules are found to be trapped inside the collapsed polymer chains leading to the varying degrees of hydration and dehydration of the polymer chain in different solutions. Calculated water diffusion coefficients for both trapped and free water molecules agree very well with experimental measurements.

*Corresponding author, Email: [email protected]; Tel: 479-575-8401

Keywords: Lower Critical Solution Temperature, Molecular Dynamics, Hydration-Dehydration, Thermo-responsive Polymer

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. INTRODUCTION Thermo-responsive polymers such as poly(N-isopropylacrylamide) (PNIPAM) and poly(Nvinylcaprolactam) (PVCL) exhibit lower critical solution temperature (LCST), above which the polymers become hydrophobic and adopt collapsed conformations, below which the polymers are hydrophilic and adopt extended coil-like conformations1-2. As a result, the LCST transition is also often called hydrophilic-to-hydrophobic transition. The LCST of a thermos-responsive polymer depends strongly on the balance of hydrophobic and hydrophilic groups on the polymer chain. Incorporating more hydrophobic groups will lead to a decrease in the LCST of the polymer. Similarly, increasing the hydrophilicity of the functional groups on the polymer chain tends to increase the LCST. The LCST transition is considered to be entropic in nature due to the competition between the hydration and dehydration free energies of the polymer chain. This hydrophilic-to-hydrophobic transition is not only affected by the temperature, but also other environmental factors such as the presence of salt. Salt ions tend to reduce the transition temperature generally. The reduction of LCST is found to be strongly dependent on the salt type and salt concentration. The higher the ionic strength is, the larger the reduction will be3-5. It is also ion specific with cations and anions each following a different order depending on the ionic charge and size. Our earlier studies1-2,

6

investigating the effects of salt ions on LCST transition of

PNIPAM show that the cations tend to interact directly with the amide O. Anionic interaction with the polymer is mediated by the cations even though larger size anions form hydrophobic interaction with the isopropyl group on PNIPAM. The LCST of PNIPAM in water is around 305 K and only depends slightly on the molecular weight of the polymer chain and polymer concentration. In contrast, the LCST of PVCL depends strongly on the polymer molecular weight and concentration7-8. For very dilute PVCL

2


Page 2 of 30

Page 3 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


concentrations, the LCST of PVCL varies from about 305 K for long polymer chains to over 323 K for short chains. Moreover, the effects of salt ions on the LCST transition of PNIPAM and PVCL are different for these two poly(amide) chains. Even though deep insights into the LCST transition of PNIPAM have been obtained from both theoretical and experimental investigations, the hydrophilic-to-hydrophobic transition of PVCL is largely unexplored, particularly at the atomistic level. This is partly due to the lack of accurate force field parameters for PVCL polymer chains. These temperature- and salt ion-responsive polymers have a wide range of applications5,

9-10

including as smart carriers for drug delivery, as responsive hydrophobic interaction chromatography ligands for protein purifications and as self-cleaning anti-fouling surfaces. Elucidating the nature of the PVCL transition has tremendous implications for these applications. Hofmeister series refers to the ability of different ions to denature protein11. Cations and anions have their respective orders. The effects of salt ions on the LCST transition of polymers as well as protein denaturation follow the Hofmeister series12-13 or inverse Hofmeister series

14-15

depending on the charge and polarity of the surfaces16. So far, the exact nature of ion specificity on the physical and biological phenomena remains elusive and under considerable debate. It is not sure whether it is caused by the changes in the bulk water structure in different salt solutions or by the direct ion-protein interactions. Previously it was thought that ions in solution can affect the bulk water structure thus affecting the hydrophilic-to-hydrophobic transition of responsive polymers in aqueous solution1718

. Ions are considered either chaotropic (structure breaking) or kosmotropic (structure making).

Previous spectroscopic and calorimetry studies show that ions only affect the closest hydration shells and the bulk water structure is not altered by the presence of the salt ions19-21. Recent theoretical and experimental work has focused on specific ion−protein interactions. Our previous

3



simulations1, 6, 22 show that this hydrophilic-to-hydrophobic transition of PNIPAM are affected by both cations and anions. Moreover, transition dynamics, polymer swellability and degree of hydration are found to be salt ion specific6. Here, more accurate force field parameters for PVCL are developed. Classical atomistic molecular dynamics (MD) simulations are subsequently performed to investigate the LCST transition of PVCL in water and several salt solutions. The salt solutions investigated include 1 M and 3.5 M NaCl and 0.5 M CaCl 2 . Other saline solutions were not investigated as these simulations are rather expensive. For the systems investigated, insights on effects of salt ions on the degree of hydration and dehydration of the PVCL polymer, its transition dynamics, water dynamics associated with the polymer and water diffusion are obtained. The radius of gyration, end-to-end distance of the polymer chain, the water molecule numbers in the first hydration shell, the residence time of the water molecules associated with the polymer as well as water diffusion coefficient are calculated. 1. COMPUTATIONAL METHODS In order to conduct the MD simulations of PVCL in water and salt solutions at different temperatures, accurate force field parameters are essential. Since, no existing force field parameters available for the VCL monomer especially the interaction parameters involving the 7member ring, quantum mechanical (QM) calculations were used to map the force field parameters prior to conducting the MD simulations. In order to determine the force field parameters, VCL monomer was terminated by two methyl groups to mimic the polymeric environment as shown in Figure 1. The modified VCL monomer was optimized at the MP2/6-31G* level in Gaussian09 23. Electrostatic potential (ESP24) was calculated at the HF/6-31G* level, and the atomic charges were

4


Page 4 of 30

Page 5 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


then derived based on restrained electrostatic potential protocol (RESP25-28) by AMBERTOOLS package29.

Figure 1

The structure of modified VCL terminated by two methyl groups mimicking the polymeric environment (C, N, and O atoms are labeled).

The bonding force field parameters were determined using SCAN in Gaussian09 at the MP2/6-31G* level. However, the seven-member ring would break during the SCAN of angles and dihedrals. PARAMFIT30 was conducted to develop the force field parameters of angles and dihedrals on the seven-member ring. PARAMFIT can adjust multiple force field parameters of angles and dihedrals simultaneously, and then fit the potentials calculated from AMBER31 using the adjusted force field parameters to the result calculated from QM calculations in Gaussian09 by linear regression. The coefficient fitting parameters with R2 larger than 0.98 from the least square regression were adopted as the force field parameters for the MD simulations of PVCL polymer chain. A PVCL chain with degree of polymerization (DP) of 100 was constructed by AMBERTOOLS package. The initial conformation of the 100 DP polymer chain was obtained by energy minimization in vacuum. The polymer chain is isotactic. The LCST of the PVCL chain with DP 100 in water was estimated based on the previous experimental studies32. The results indicate that LCST decreases linearly with the increase of molecular weight (MW) when the 5



Page 6 of 30

polymer MW is less than 22100. Since DP 100 PVCL has a MW of ~13950, its LCST in water is estimated to be around 317 K. The LCSTs of the polymer chain in salt solutions were estimated based on previous experimental results 7, 33-34. Table 1 shows experimentally estimated LCST for the 100 DP polymer in water and different salt solutions as well as the simulation temperatures adopted. The simulation temperatures chosen were above and below the LCSTs of the respective systems. For PVCL in 3.5 M NaCl, one temperature at 283 K close to the transition temperature was also investigated. Table 1

PVCL systems investigated using atomistic MD simulations

System

H2O

PVCL w/w Concentration (%) Estimated LCST (K) from Experiments Simulation Temperature (K) Number of Cations

2.19

0.5 M CaCl 2

1 M NaCl

3.5 M NaCl

2.12

2.12

1.87

304-314

300-310

295-305

273-283

283 & 353

273 & 333

270 & 350

263, 283 & 353

0

300

600

2244

Number of Anions

0

600

600

2244

Number of Water Molecules Unit Cell Size (Å3)

34647

33850

33873

33403

100×112×91

99×111×91

99×111×91

100×95×111

The simulations were conducted with AMBER package35. TIP3P was used as the water model36. The unit cell dimensions, the number of water molecules and the corresponding cation and anion numbers in each system, the weight-based polymer concentration and the simulation temperatures were also listed in Table 1. The polymer chain was placed diagonally to minimize its the potential interactions between the polymer and its neighboring images. The initial distance between the PVCL chain and the edge of the unit cells was kept at 10 Å. Before the production

6


Page 7 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


run, a total of 1000 steps of energy minimization were performed by allowing the water molecules to hydrate the restrained polymer. Additional 2500 steps of energy minimization were conducted by allowing all the atoms to relax and equilibrate. This was followed by increasing the temperature to a targeted value in 20 picosecond (ps) of time while the polymer was being restraining. Finally all simulations were run at the constant temperature and pressure (1 bar) (NPT) under the periodic boundary condition with the Langevin dynamics thermostat

37

. Ewald summation was used to

determine the electrostatic interactions. A 10 Å cutoff was used for calculating both the electrostatic and van der Waals (VdW) interactions. The time step was chosen to be 1 femtosecond (fs). The force field parameters of Na+, Cl- were from ion0838 in AMBER. Improved force field parameters were used for Ca2+, in which an additional term was added to Lennard-Jones (LJ) potential to construct LJ12-6-4 VdW energies considering the charge-induced dipole interaction, which has been validated by experimental results of hydration free energy, Ca2+-O (H 2 O) distance and coordination number39. The number of water molecules in the first hydration shell is counted if the water molecules are within 3.3 Å from the surface atoms on the PVCL chain. This cut-off distance is based on the first minimum on the amide O-O (H 2 O) radial distribution function (RDF). RDFs were calculated between the water molecules with various heavy atoms on the PVCL chain as well as between the salt ions and the polymer chain. Since these simulations are computational expensive, the production run time for each system is somewhat different due to its different LCST transition kinetics. At the temperatures above the respective LCSTs, at least 50 ns of simulations was performed after the polymer chain collapsed. As a result, about 100-200 ns of simulations was conducted for each system.

7



2. RESULTS AND DISCUSSION

Figure 2

The radius of gyration (top left), end-to-end distance (top right), the water number

in the first hydration shell (bottom left) and the conformations of PVCL at the end of the simulations at temperatures above (353 K) and below (283 K) its LCST in water.

A total of about 200 ns simulations were conducted for the PVCL chain in water at 283 K and 353 K below and above its LCST respectively. Figure 2 shows clockwise the radius of gyration (R g ), end-to-end distance, the final conformations at the end of the simulation, and the water number in the first hydration shell respectively for the two simulations. The initial conformations for both simulations are stretched and coil-like. At high and low temperatures, both the R g values start to decrease from the initial value of 33 Å to about 25 Å as the simulations continue until about

8


Page 8 of 30

Page 9 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


100 ns following a similar trend. Thereafter, the R g value for the PVCL chain at 283 K remains more or less the same or even increases slightly during the remaining 100 ns of the simulation. On the other hand, for the simulations performed at 353 K, the R g value of the PVCL chain continues to decrease reaching ~18 Å at 150 ns and stabilizing at this value for the rest of the simulation period. The end-to-end distances for the two systems follow a similar trend. The end-to-end distance of the PVCL chain at 353 K increases initially due to the high thermal energies. However, it decreases steadily from the maximum value of about 120 Å to about 50 Å at the end of the 200 ns simulation. At 283 K, the end-to-end distance also decreases steadily from the initial value of about 80 Å to about 30 Å at the end of the simulation. The end-to-end distance is not necessarily accurate to describe the stretched or collapsed state of the polymer as can be seen from the conformations shown at the end of the simulations. The polymer chain at 283 K is clearly only partially folded whereas the polymer is completely collapsed at 353 K above its LCST. The best indication of the LCST phase transition is the water numbers in the first hydration shell of the polymers above and below its LCST as also shown in Figure 2. The water molecules hydrating the PVCL chain is significantly higher at 283 K (~950 at the end of 200 ns simulation) than those at 253 K (~700) indicating the different degrees of hydration of the two polymer chains. Even though the polymer conformations at the beginning of the simulations at two temperatures are the same, the water molecule numbers in the first hydration shell are different even at the beginning of the production run. This is due to the energy minimization steps and 20 ps of pre-equilibration at different targeted temperatures before the production run. Our MD simulations with PVCL force field parameters developed in-house can clearly demonstrate the hydrophilic-to-hydrophobic transition of PVCL in water when the temperature increases above its LCST in water.

9



Figure 3

The ring distribution in the stretched region (left) and folded region (right) for a

100 DP PVCL chain. The positions of the 7-member VCL ring around the backbone are also shown for the two regions.

As discussed, the hydrated PVCL has a coil-like conformation below its LCST. When the temperature increases above its LCST, it becomes dehydrated leading to the collapse of the polymer chain. As this hydrophilic-to-hydrophobic transition involves dramatic conformational changes, the VCL 7-member rings along the backbone also need to undergo dramatic rearrangement due to the steric hinderance effects. The VCL ring distribution along the backbone during the LCST phase transition is shown in Figure 3. It can be seen that the VCL ring distributions along the backbones of the stretched and folded fragments have significant differences. In the stretched region, 5 to 20 VCL rings are distributed along the backbone evenly to form a helix structure due to the bulkiness of the 7-member VCL ring. The backbone maintains the rod-like linear conformation as shown in the left panel of Figure 3. The folded fragments on the chain are much shorter. These VCL rings all face the water molecules while the backbone collapses onto each other from the opposite sides to form a folded local structure as shown in the

10


Page 10 of 30

Page 11 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


right panel of Figure 3. The VCL rings in the folded region need to be arranged in a specific way as described above to drive the hydrophilic-to-hydrophobic conformational change. A completely folded globular PVCL chain is comprised of several folded regions that are formed by alternating longer linear and shorter curved fragments.

Figure 4


in the first hydration shell (bottom left) and the conformations of PVCL at the end of the 200 ns simulations at temperatures above (330 K) and below (270 K) its LCST in 1 M NaCl solution.

The MD simulations of the PVCL chain in 1 M NaCl solution were conducted for about 195 ns at 270 K and 330 K below and above its LCST respectively. The four panels in Figure 4 show in clockwise R g , the end-to-end distance, the polymer conformation at the end of the 11



simulation and the water molecule number in the first hydration shell respectively. In 1 M NaCl solution, the R g values of the two systems follow more or less the same trend. The R g values exhibit a rapid increase from 30 to 35 Å during the first few ns of simulation followed by a steady decrease to around 22 Å at about 80 ns. The R g values keep more or less the same for the remaining simulation period. The R g value at 330 K was observed to be slightly higher than the corresponding value at 270 K. Both R g value and end-to-end distance of a polymer chain are two different commonly used measures of the polymer conformation. However, they are not necessarily the most reliable measure for the degree of hydration or dehydration. Here the end-to-end distances of the PVCL chains have substantial differences after 60-70 ns of simulations. The two ends of the PVCL chain at 270 K are far away from each other at around 40 Å, whereas the end-to-end distance at 330 K is around 10 Å indicating the formation of the dehydrated state above LCST. At the end of the simulations, both chains show a more or less a folded conformation. However, the degree of hydration for the two temperatures are rather different. At 270 K, the water number in the first hydration shell only decreases slightly from 1100 to about 1050 during the first 80 ns of simulations whereas the R g values and the end-to-end distances demonstrate significant decrease. During the subsequent simulations at 270 K, the PVCL chain becomes partially folded with the water number in the first hydration shell decreases from 1050 to 900. On the other hand, at 330 K, the water numbers in the first hydration shell decreases from the initial 900 to about 750 after 70 ns of simulations and keeps the same number till the end of the simulation. This result indicates that the PVCL chain is much more hydrated below its LCST even though it is also partially folded with a similar R g value and end-to-end distance to the polymer chain above its LCST. Our results also indicate that the water number in the first hydration shell is a reliable measure of the degree of hydration for PVCL.

12


Page 12 of 30

Page 13 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


From the simulations of PVCL in water and 1 M NaCl solution above their respective LCSTs, it can be seen that the dynamics of the LCST transition from a coil-like conformation to a collapsed globule conformation is different under these two solutions. The LCST transition takes much longer in water than in 1 M NaCl salt solution. In addition, the degrees of hydration or dehydration of the PVCL chain is also different in water and 1 M NaCl solution. Earlier studies have also shown that the LCST transition of thermo-responsive PNIPAM is strongly ion specific1, 6, 40

. In order to further investigate the ion-specificity of the LCST transition, MD simulations were

conducted for PVCL chain in 0.5 M CaCl 2 solution for a total of about 120 ns simulations at 273 K and 333 K above and below its estimated experimental LCST respectively. Figure 5 shows in clockwise the radius of gyration, the end-to-end distance, the final conformation of the polymer and the water number in the first hydration shell during or after about 120 ns simulations. It can be seen that the LCST transition at 333 K occurs rapidly in less than 40 ns with the R g value decreased from the initial 30 Å to below 20 Å whereas the PVCL chain remains stretched when the simulation was conducted at 273 K. The R g value even increases to above 40 Å at 273 K. Similarly, the end-to-end distances exhibit significant differences above and below its LCST after 50 ns simulations. The water molecule number in the first hydration shell is over 1000 for PVCL chain in 0.5 M CaCl 2 during the entire simulation period at 333 K whereas it decreases to just about 700 at 333 K. The conformation of the PVCL demonstrates a completely folded structure at 333 K whereas it is linear and stretched out at 273 K at the end of the 120 ns simulations. Our simulations results show that the transition dynamics and degree of hydration/dehydration indeed depend on the specific salt ions.

13



Figure 5


in the first hydration shell (bottom left) and the conformations of PVCL at the end of the 150 ns simulations at temperatures above (333 K) and below (273 K) its LCST in 0.5 M CaCl 2 solution. In addition to the ion type, ion concentration also affects LCST transition significantly5, 10. For most of the cations and anions, the higher the concentration is, the lower the transition temperature becomes. In order to investigate the effects of ionic concentration on the LCST transition of PVCL, MD simulations of a 100 DP PVCL in 3.5 M NaCl solution at temperatures above (353 K), below (263 K) and close to 283 K its LCST were conducted. Figure 6 shows similarly clockwise the radius of gyration, the end-to-end distance, the final conformation and the water number in the first hydration shell for the simulation period conducted. It is found that 3.5 M NaCl salt solution speeds up the transition process significantly. At 353 K above the LCST, the hydrophilic-to-hydrophobic transition occurred during the first 40 ns. The R g decreased to 15 Å 14


Page 14 of 30

Page 15 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from the initial value of 30 Å. The R g value remains more or less the same during the subsequent 50 ns simulations. At 283 K which is close to the LCST transition temperature, the hydrophilicto-hydrophobic collapsing occurred after about 90 ns of simulation reaching a R g value of about 16 Å, similar to the corresponding value at 353 K. On the other hand, the PVCL chain at 263 K maintained a coil-like conformation with a R g value of around 35 Å throughout the simulation. Both PVCL chains at 283 K and 353 K went through the LCST transition and exhibited a collapsed conformation at the end of the simulation with their end-to-end distances decreasing from the initial ~60 Å to ~30 Å.

Figure 6


in the first hydration shell (bottom left) and the conformations of PVCL at the end of simulation times at temperatures above (353 K), below (263 K) and close to 283 K its LCST in 3.5 M NaCl solution.

15



On the other hand, the three systems demonstrate substantial differences in their water numbers in the first hydration shell of the PVCL chain. At 353 K, the water number decreased from the initial 900 to about 630 during the first 40 ns of simulations when the polymer goes through the hydrophilic-to-hydrophobic transition. At 283 K, the water number in the first hydration shell only shows a slight but steady decrease over the simulation period from the initial 1000 to about 800 at the end of the simulation even though the R g values exhibit a clear conformational transition during the first 90 ns of simulation. At 263 K, the water number in the first hydration shell of PVCL decreases slightly but is always higher than the corresponding water number at 283 K. It appears that even though PVCL goes through a hydrophilic-to-hydrophobic transition forming a folded hydrophobic conformation at both 353 K and 283 K, the degree of hydration reflected by the water molecule numbers associated with the polymer remains to be significantly different. The LCST transition is not strictly a two-state transition, i.e. either hydrophobic or hydrophilic. At collapsed hydrophobic conformation, there could still be substantial differences in their hydrophobicity or water molecule numbers associated with the polymer. Similarly, at hydrophilic state, water number hydrating the polymer can also vary. More details about different degrees of hydration and dehydration will be discussed below.

16


Page 16 of 30

Page 17 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7

The combined radii of gyration (left) and water numbers in the first hydration water

shells (right) for the nine (9) systems investigated. The R g values and water numbers in the first hydration shell of the PVCL chain in all 9 systems are plotted together and compared as shown in the left and right panels of Figure 7 respectively. The R g values for the coil-like PVCL chain below its LCST are around 23, 43, 20, and 37 Å respectively in water (283 K), 0.5 M CaCl 2 (273 K), 1 M NaCl (270 K) and 3.5 M NaCl (263 K). The corresponding R g values for the collapsed PVCL above its LCST are around 16, 17, 23, and 15 Å respectively in water (353 K), 0.5 M CaCl 2 (333 K), 1 M NaCl (330 K) and 3.5 M NaCl (353 K). Clearly, the size of the PVCL chain is different in different salt solutions. The numbers of water molecules in the first hydration shell below LCST are shown in the right panel of Figure 7. It can be seen that water molecule number in the first hydration shell varies when the PVCL chain is either in the hydrated state or in the dehydrate state. There is no single number to represent the number of water molecules associated with the polymer chain. As discussed previously, the dynamics of the hydrophilic-to-hydrophobic transition is strongly salt type and salt concentration dependent. Our results further demonstrate that the degree of hydration and dehydration is also salt type and salt concentration dependent for both hydrophilic and hydrophobic states. 17



Figure 8

Radial distribution functions between the water molecules and various atoms on the

7-member ring of PVCL chain for the simulations conducted in water at 283 K.

Radial distribution functions (RDFs) between the polymer and H 2 O were calculated over a period of 10 ns at 283 K after the system reached quasi-equilibrium. The RDFs between O on the H 2 O molecule and the C, N, O atoms on the 7-member ring of the PVCL chain are shown in Figure 8. A strong correlation peak at 2.8 Å is observed for O-O1 between the oxygen atom on the amide bond and the oxygen atom on the water molecule indicating the formation of strong cooperative hydrogen bonds41 between the polymer and surrounding water molecules. The correlation O-C5 observed with a peak at 3.6 Å between O (H 2 O) and amide C5 atom is likely due to the proximity of the C5 atom to O1. Other atoms on the PVCL ring have rather weak correlations with the solvent water molecules. Other weak correlations observed are also likely due to the strong hydrogen bonding formation between the water molecule and amide O1. Even though amide N is negatively charged, its interaction with water molecules is really weak with a correlation peak at 4.8 Å due to the proximity of the N1 to the amide O1. This agrees well with our earlier observation of PNIPAM in water1.

18


Page 18 of 30

Page 19 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9

Radial distribution functions between the cations (Na+, Ca2+) and amide O1 on the

7-member ring of the PVCL chain in 1 M NaCl, 3.5 M NaCl and 0.5 M CaCl 2 salt solutions.

The RDFs between the cations (anions) and the amide O were also calculated following the same procedure as that for the RDFs between water molecules and atoms on the polymer chain. Figure 9 shows the cation-amide O correlation (top) and anion-amide O (bottom) respectively. A correlation peak at about 2.3 Å was observed for the Na+–amide O interaction in both 1 M and 3.5 M NaCl solutions at 263, 270, 283, 330 and 353 K. A correlation peak at about 2.4 Å was observed for the Ca2+–amide O interaction at 353 K, but not at 273 K. Our previous theoretical studies1-2, 6 19



on the effects of salt ions on PNIPAM demonstrate that the hydration of ions tend to decrease the LCST whereas the direct cation–amide O binding tends to increase the LCST. The strength of the cation–amide O interaction is dictated by the competition between the electrostatic and the hydration forces. For the singly charged alkali cations (Li+, Na+, K+ and Rb+), electrostatic interactions dominate. The larger the cation is, the weaker the binding interaction with the amide O. For the doubly charged cations (Mg2+ and Ca2+) , the strong hydration of these divalent ions overcomes the electrostatic attraction leading to the very weak binding between the cations and the amide O. For Ca2+–amide O interaction, the correlation tends to be stronger at higher temperature due to the reduced hydration free energy at elevated temperatures. On the contrary, electrostatic interaction is less affected by the temperature. For the anion-amide O interaction as seen from the bottom panel of Figure 9, no strong correlation was observed for all the systems. This observation was consistent with our previous results1, 6.

Figure 10

The residence time of water molecules interacting with PVCL chain during 10 ns

interval at the start and end of the simulations for PVCL chain in water at 283 K.

20


Page 20 of 30

Page 21 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Estimating water residence time is critical to understand the water dynamics. The residence time of water molecules in the first hydration shell surrounding each residue was calculated taking an average of 10 ns simulation of the stretched or folded PVCL chain in water at 283 and 353 K at the beginning and end of the simulations as shown in Figures 10 and 11 respectively. The water molecules are considered in the first hydration shell when the distance between the O on the H 2 O molecule and amide O is less than the first minimum at 3.3 Å based on the correlation in Figure 8. A detailed procedure for calculating the water residence time was described in our previous work42. At 283 K below its LCST, the PVCL chain is coil-like during the initial 10 ns period. The water residence time around each residue fluctuates in a narrow range between 55 to 80 pico-second (ps) as shown in Figure 10. At the end of the simulation, fragments surrounding residues 30, 50, and 90 have been found to be partially folded even though the chain remains to be in a hydrated state. The water residence times around these residues increase substantially reaching 100-200 ps indicating that some water molecules have been trapped in these folded fragments. On the other hand, the water residence times surrounding the other unfolded residues remain in the range 5080 ps.

Figure 11

The residence time of water molecules interacting with PVCL chain during 10 ns

interval at the start and end of the simulations for PVCL chain in water at 353 K.

21



Figure 11 shows that water dynamics at temperature above its LCST at 353 K is rather different from the temperature below its LCST. Water movement is much faster at 353 K with a water residence time in the range of 20-35 ps at the beginning of the simulation when the polymer is still coil-like as seen from the left panel of Figure 11. However, towards the end of the simulation during 148-158 ns and when the polymer has gone through the hydrophilic-to-hydrophobic transition to form a collapsed folded conformation, water residence times surrounding the residues show a very different characteristic. For water molecules surrounding some residues, particularly for residue number 60-80, water movement was found to be drastically slower than other residues. Water residence time increases to 100-240 ps similar to the lower temperature case when water molecules were being trapped in the folded chain. In a collapsed conformation, these residues were buried inside the folded polymer chain. Water dynamics are found to be much slower in other regions as well. However, water residence times surrounding some other residues which are exposed to the bulk water remain fast in the range of 20-30 ps. Not only for PVCL in water, PVCL in salt solutions exhibits similar behavior. The residence time of water molecules around coiled PVCL always fluctuates at a faster speed, whereas water dynamics around folded PVCL chain becomes much slower at some residues indicating that water molecules have been trapped in certain folded pockets. Here if the residence time of a water molecule exceeds 1 ns, it is considered as trapped water. Water molecules with residence time less than 1 ns are considered free water. Self-diffusion coefficients of trapped and free water molecules were calculated based on the mean square displacement of the molecules during the simulations.

22


Page 22 of 30

Page 23 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12

The conformation of a partially folded PVCL chain and the trapped water

molecules inside the polymer for PVCL in water.

The trapped water molecules are found in hydrophobic collapsed structures or partially folded fragments as shown in Figure 12. During the transition from the hydrated coil-like conformation to the dehydrated folded conformation, water molecules surrounding the PVCL chain are being removed. Depending on the salt type, salt concentration and temperature, the degree of dehydration or the number of water molecules removed appears to be somewhat different. Trapped water molecules form fewer hydrogen bonds with the neighboring water molecules, instead they form hydrogen bonds with the PVCL chain. These water molecules can still enter or leave the trapped regions, but their motion is substantially restricted. To better understand these structured water molecules, the nearest neighbor coordination numbers with other water molecules for both free and trapped water molecules were calculated as shown in Table 2. The nearest neighbor coordination numbers of free and trapped water molecules are 3.62 and 2.22 respectively, which indicates that trapped water molecules probably do not form a 3-D tetrahedral water network, instead they likely form a 2-D water structure.

23



Table 2 The nearest neighbor coordination numbers for free and trapped water molecules calculated from 1 ns simulation of PVCL chain in 3.5 M NaCl at 263 K. Water Position

Coordination Number

Trapped Water

2.22

Free Water

3.62

Previous studies43 used 1H-NMR to measure the diffusion coefficients of trapped water molecules in hydrogels made of two compositions of PVCL and PNIPAM. Water diffusion coefficients in hydrogels of 54/46 and 50/50 of PVCL/PNIPAM were measured. Both experimental diffusion coefficients and calculated values are shown in Table 3. Since the experimental measurements were performed at 320 K, water diffusion coefficient was scaled to the simulation temperature of 283 K based on Stokes–Einstein equation. The viscosity used was from IAPWS44. The measured diffusion coefficient for the free water is 1.12~1.18×10-5 cm2/s at 283 K. The calculated value from MD simulations is 2.27×10-5 cm2/s in very good agreement with the experimental value. The measured diffusion coefficient for the trapped water molecules is 0.03~0.35×10-5 cm2/s at 283 K. The calculated diffusion coefficient from current MD simulations is 0.08×10-5 cm2/s also in very good agreement with the experimental value. It can be seen that the diffusion coefficient of trapped water is about 1~2 orders of magnitude smaller than that of free water indicating the movement of trapped water molecules is significantly restricted. This result agrees with the previous water residence time and coordination number calculations.

Table 3 The comparison of the diffusion coefficient of free and trapped water between experiment and simulation. The system of simulation result is in water at 283 K. 24


Page 24 of 30

Page 25 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


System

Temperature (K)

D (Free Water) (cm2/s)

D (Trapped Water) (cm2/s)

Experiment62

320

2.85~3.02×10-5

0.08~0.89×10-5

Experiment (scaled)

283

1.12~1.18×10-5

0.03~0.35×10-5

MD Simulation

283

2.27×10-5

0.08×10-5

3. CONCLUSIONS The MD simulations of hydrophilic-to-hydrophobic transition of a 100 DP responsive PVCL chain were successfully conducted for the first time in water, 1 M NaCl, 3.5 M NaCl as well as in 0.5 M CaCl 2 salt solutions based on in-house developed force field parameters from ab initio calculations. The 7-member rings on the PVCL chain are found to arrange differently in the coillike and collapsed regions. The transition dynamics in different salt solutions are ion specific. The radius of gyration, the water molecule number in the first hydration shell, the degree of hydration and dehydration for both the hydrophilic and hydrophobic states are found to be salt ion type and salt concentration dependent. This difference in hydration or dehydration degrees is likely due to the presence and number of trapped water molecules inside the collapsed polymer chains or partially folded regions.

4. ACKNOWLEDGEMENT Funding from University of Arkansas is gratefully acknowledged. Computations were performed on Razor2 from the Arkansas High Performance Computing Center (AHPCC) from the University of Arkansas.

25



5. REFERENCE

1.

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 (49), 16594-16604. 2.

Du, H.; Qian, X., The Interactions between salt ions and thermo-responsive poly (N-

Isopropylacrylamide) from molecular dynamics simulations. In Responsive Membranes and Materials, John Wiley & Sons, Ltd.: 2012; pp 229-242. 3.

Zhang, Y.; Cremer, P. S., Interactions between macromolecules and ions: the Hofmeister

series. Curr. Opin. Chem. Biol. 2006, 10 (6), 658-663. 4.

Hiruta, Y.; Nagumo, Y.; Suzuki, Y.; Funatsu, T.; Ishikawa, Y.; Kanazawa, H., The effects

of anionic electrolytes and human serum albumin on the LCST of poly(N-isopropylacrylamide)based temperature-responsive copolymers. Colloids and Surfaces B: Biointerfaces 2015, 132 (0), 299-304. 5.

Liu, Z.; Wickramasinghe, S. R.; Qian, X., Ion-specificity in protein binding and recovery

for the responsive hydrophobic poly (vinylcaprolactam) ligand. RSC Advances 2017, 7 (58), 36351-36360. 6.

Du, H.; Wickramasinghe, S. R.; Qian, X., Specificity in cationic interaction with poly(N-

isopropylacrylamide). J. Phys. Chem. B 2013, 117 (17), 5090-5101. 7.

Meeussen, F.; Nies, E.; Berghmans, H.; Verbrugghe, S.; Goethals, E.; Du Prez, F., Phase

behaviour of poly(N-vinyl caprolactam) in water. Polymer 2000, 41 (24), 8597-8602. 8.

Spevacek, J.; Dybal, J.; Starovoytova, L.; Zhigunov, A.; Sedlakova, Z., Temperature-

induced phase separation and hydration in poly(N-vinylcaprolactam) aqueous solutions: a study by NMR and IR spectroscopy, SAXS, and quantum-chemical calculations. Soft Matter 2012, 8 (22), 6110-6119. 9.

Vu, A.; Qian, X.; Wickramasinghe, S. R., Membrane-based hydrophobic interaction

chromatography. Separation Science and Technology 2017, 52 (2), 287-298. 10.

Liu, Z.; Wickramasinghe, S. R.; Qian, X., The architecture of responsive polymeric ligands

on protein binding and recovery. RSC Advances 2017, 7 (44), 27823-27832. 11.

Hofmeister, F., Zur Lehre von der Wirkung der Salze. Archiv f. experiment. Pathol. u.

Pharmakol 1888, 24 (4-5), 247-260.

26


Page 26 of 30

Page 27 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


12.

Melander, W.; Horváth, C., Salt effects on hydrophobic interactions in precipitation and

chromatography of proteins: An interpretation of the lyotropic series. Arch. Biochem. Biophys. 1977, 183 (1), 200-215. 13.

Nfor, B. K.; Hylkema, N. N.; Wiedhaup, K. R.; Verhaert, P. D. E. M.; van der Wielen, L.

A. M.; Ottens, M., High-throughput protein precipitation and hydrophobic interaction chromatography: Salt effects and thermodynamic interrelation. Journal of Chromatography A 2011, 1218 (49), 8958-8973. 14.

Boström, M.; Parsons, D. F.; Salis, A.; Ninham, B. W.; Monduzzi, M., Possible origin of

the inverse and direct Hofmeister series for lysozyme at low and high salt concentrations. Langmuir 2011, 27 (15), 9504-9511. 15.

Paterová, J.; Rembert, K. B.; Heyda, J.; Kurra, Y.; Okur, H. I.; Liu, W. R.; Hilty, C.;

Cremer, P. S.; Jungwirth, P., Reversal of the hofmeister series: specific ion effects on peptides. J. Phys. Chem. B 2013, 117 (27), 8150-8158. 16.

Schwierz, N.; Horinek, D.; Netz, R. R., Reversed anionic Hofmeister series: the interplay

of surface charge and surface polarity. Langmuir 2010, 26 (10), 7370-7379. 17.

Vanzi, F.; Madan, B.; Sharp, K., Effect of the protein denaturants urea and guanidinium on

water structure: a structural and thermodynamic study. Journal of the American Chemical Society 1998, 120 (41), 10748-10753. 18.

Zou, Q.; Bennion, B. J.; Daggett, V.; Murphy, K. P., The molecular mechanism of

stabilization of proteins by TMAO and its ability to counteract the effects of urea. Journal of the American Chemical Society 2002, 124 (7), 1192-1202. 19.

Omta, A. W.; Kropman, M. F.; Woutersen, S.; Bakker, H. J., Negligible effect of ions on

the hydrogen-bond structure in liquid water. Science 2003, 301 (5631), 347-349. 20.

Batchelor, J. D.; Olteanu, A.; Tripathy, A.; Pielak, G. J., Impact of protein denaturants and

stabilizers on water structure. Journal of the American Chemical Society 2004, 126 (7), 1958-1961. 21.

Omta, A. W.; Kropman, M. F.; Woutersen, S.; Bakker, H. J., Influence of ions on the

hydrogen-bond structure in liquid water. Journal of Chemical Physics 2003, 119 (23), 1245712461. 22.

Du, H.; Qian, X., Molecular dynamics simulations of PNIPAM-co-PEGMA copolymer

hydrophilic to hydrophobic transition in NaCl solution. Journal of Polymer Science Part B: Polymer Physics 2011, 49 (15), 1112-1122.

27



23.

Frisch, M.; Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Scalmani, G.;

Barone, V.; Mennucci, B.; Petersson, G., Gaussian 09, revision D. 01. Gaussian, Inc., Wallingford CT: 2009. 24.

Singh, U. C.; Kollman, P. A., An approach to computing electrostatic charges for

molecules. Journal of Computational Chemistry 1984, 5 (2), 129-145. 25.

Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A., A well-behaved electrostatic

potential based method using charge restraints for deriving atomic charges: the RESP model. The Journal of Physical Chemistry 1993, 97 (40), 10269-10280. 26.

Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollmann, P. A., Application of RESP charges to

calculate conformational energies, hydrogen bond energies, and free energies of solvation. Journal of the American Chemical Society 1993, 115 (21), 9620-9631. 27.

Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A., Application of the multimolecule

and multiconformational RESP methodology to biopolymers: Charge derivation for DNA, RNA, and proteins. Journal of Computational Chemistry 1995, 16 (11), 1357-1377. 28.

Fox, T.; Kollman, P. A., Application of the RESP methodology in the parametrization of

organic solvents. The Journal of Physical Chemistry B 1998, 102 (41), 8070-8079. 29.

Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A., Automatic atom type and bond type

perception in molecular mechanical calculations. Journal of molecular graphics and modelling 2006, 25 (2), 247-260. 30.

Betz, R. M.; Walker, R. C., Paramfit: automated optimization of force field parameters for

molecular dynamics simulations. Journal of computational chemistry 2015, 36 (2), 79-87. 31.

Salomon‐Ferrer, R.; Case, D. A.; Walker, R. C., An overview of the Amber biomolecular

simulation package. Wiley Interdisciplinary Reviews: Computational Molecular Science 2013, 3 (2), 198-210. 32.

Lyon, L. A.; Serpe, M. J., Hydrogel micro and nanoparticles. John Wiley & Sons: 2012.

33.

Maeda, Y.; Nakamura, T.; Ikeda, I., Hydration and phase behavior of poly (N-

vinylcaprolactam) and poly (N-vinylpyrrolidone) in water. Macromolecules 2002, 35 (1), 217222. 34.

Yuan, G.; Zhao, Z.; Li, W.; Yin, L.; Wang, L., Correlation between temperature-dependent

dissolution of poly (N-vinylcaprolactam)/poly (Nl-(1-hydroxymethyl) propylmethacrylamide) layer-by-layer films and cloud point of mixed solutions of the two polymers in the presence of

28


Page 28 of 30

Page 29 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


chloride salts. Colloids and Surfaces A: Physicochemical and Engineering Aspects 2013, 436, 1186-1191. 35.

Case, D.; Darden, T.; Cheatham, T.; Simmerling, C.; Wang, J., PA AMBER 11. San

Francisco: University of California 2010. 36.

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L.,

Comparison of simple potential functions for simulating liquid water. The Journal of chemical physics 1983, 79 (2), 926-935. 37.

Chandrasekhar, S., Stochastic problems in physics and astronomy. Reviews of modern

physics 1943, 15 (1), 1. 38.

Joung, I. S.; Cheatham III, T. E., Determination of alkali and halide monovalent ion

parameters for use in explicitly solvated biomolecular simulations. The journal of physical chemistry B 2008, 112 (30), 9020-9041. 39.

Li, P.; Merz Jr, K. M., Taking into account the ion-induced dipole interaction in the

nonbonded model of ions. Journal of chemical theory and computation 2013, 10 (1), 289-297. 40.

Patra, L.; Vidyasagar, A.; Toomey, R., The effect of the Hofmeister series on the

deswelling isotherms of poly(N-isopropylacrylamide) and poly(N,N-diethylacrylamide). Soft Matter 2011, 7 (13), 6061-6067. 41.

Qian, X., The effect of cooperativity on hydrogen bonding interactions in native cellulose

Iβ from ab initio molecular dynamics simulations. Molecular Simulation 2008, 34 (2), 183-191. 42.

Du, H.; Qian, X., The hydration properties of carboxybetaine zwitterion brushes. Journal

of computational chemistry 2016, 37 (10), 877-885. 43.

Hanyková, L.; Spěváček, J.; Radecki, M.; Zhigunov, A.; Šťastná, J.; Valentová, H.;

Sedláková, Z., Structures and interactions in collapsed hydrogels of thermoresponsive interpenetrating polymer networks. Colloid and Polymer Science 2015, 293 (3), 709-720. 44.

Wagner, W.; Cooper, J.; Dittmann, A.; Kijima, J.; Kretzschmar, H.-J.; Kruse, A.; Mares,

R.; Oguchi, K.; Sato, H.; Stocker, I., The IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam. Journal of engineering for gas turbines and power 2000, 122 (1), 150-184.

29



TOC Graphic

30


Page 30 of 30

Atomistic Molecular Dynamics Simulations of the Lower Critical

Recommend Documents