Thermosensitive Hydration of Four Acrylamide-Based Polymers in Coil and Globule Conformations.

To characterize the thermosensitive coil-globule transition in atomistic detail, the conformational dynamics of linear polymer chains of acrylamide-based polymers have been investigated at multiple temperatures. Therefore, molecular dynamic simulations of 30mers of polyacrylamide (AAm), poly-N-methylacrylamide (NMAAm), poly-N-ethylacrylamide (NEAAm), and poly-N-isopropylacrylamide (NIPAAm) have been performed at temperatures ranging from 250 to 360 K for 2 μs. While two of the polymers are known to exhibit thermosensitivity (NEAAm, NIPAAm), no thermosensitivity is observed for AAm and NMAAm in aqueous solution. Our computer simulations consistently reproduce these properties. To understand the thermosensitivity of the respective polymers, the conformational ensembles at different temperatures have been separated according to the coil-globule transition. The coil and globule conformational ensembles were exhaustively analyzed in terms of hydrogen bonding with the solvent, the change of the solvent accessible surface, and enthalpic contributions. Surprisingly, independent of different thermosensitive properties of the four polymers, the surface affinity to water of coil conformations is higher than for globule conformations. Therefore, polymer-solvent interactions stabilize coil conformations at all temperatures. Nevertheless, the enthalpic contributions alone cannot explain the differences in thermosensitivity. This clearly implies that entropy is the distinctive factor for thermosensitivity. With increasing side chain length, the lifetime of the hydrogen bonds between the polymer surface and water is extended. Thus, we surmise that a longer side chain induces a larger entropic penalty due to immobilization of water molecules.

Introduction

Since their first discovery, thermosensitive polymers (TSPs) have been the subject of interest in many different fields of research.[1−4] Next to medical applications, such as tissue engineering, bioseparation, and drug delivery,[5−10] they have been proven to be highly valuable as programmable materials[11] and for gel actuators amongst other things. Their development has advanced to the point that it is now possible to use TSPs as a remote-controlled drug delivery platform in cancer therapy.[12] Crucial for the majority of the applications is the extraordinary phase behavior of these polymers. Counterintuitively, depending on solvent conditions, TSPs undergo a liquid–gel-phase transition with a lower critical solution temperature (LCST). Thus, above the associated phase transition temperature, TSPs form a viscoelastic gel, whereas below, they are found to be in a liquid mixture with the solvent.[13]

The above-mentioned macroscopic liquid-gel phase transition of TSPs has been linked to a microscopic conformational rearrangement, i.e., the coil–globule transition (CGT).[14−16] While the liquid–gel-phase transition can be measured by clear changes in the properties of the bulk fluid, such as the change of opacity, the conformational CGT can be reproduced in experiments in dilute solution, without necessarily exhibiting a phase transition of the liquid mixture.[15,17−20] Since detailed information about polymer conformations is extremely challenging to be obtained experimentally, this transition is solely defined by the change in the size distribution of the polymer chains in solution. At “low temperatures”, i.e., below the CGT temperature T < T*, TSPs exist in extended conformations, indicated by a large radius of gyration (Rg). However, at “high temperatures”, T > T*, the polymer chains collapse, exhibiting a clearly lower distribution of Rg. Conventionally, the extended conformational state is called coil (C), whereas the collapsed state is called globule (G).[13,21,22]

Experimental studies have been complemented by computational investigations to study the conformational changes of polymers on an atomistic level. Molecular dynamic (MD) simulations have been established as the state-of-the-art method.[23,24] Nevertheless, there are multiple reasons, why characterizing the CGT is computationally challenging: One important aspect is that none of the two conformational states, i.e., neither C nor G, are defined by distinct structural descriptors. Consequently, the conformational space is large in both conformational states, the C and the G subensembles.[25] Hence, the CGT is fundamentally different from protein-folding processes in that respect.[26] The absence of a well-defined fold leads to an insuperable uncertainty of whether the captured conformational subspace comprises the most favorable conformational states, i.e., the lowest free energy minima. Thus, both states can be very diverse conformationally. In fact, there is experimental evidence for conformational substates within the G, which can be separated kinetically.[27]

Due to the slow timescale of the transition between C and G, it is a challenge to achieve sufficient conformational sampling. Hence, it is not trivial to assess how large the variance of thermodynamic and structural properties within C and G may be. Furthermore, depending on the simulation temperature, it is very challenging to sample a large number of transitions. Therefore, reliable ensemble averages can only be estimated from long simulations, in particular, at low simulation temperatures.[28−30] The slow timescale of the conformational transition translates to a large barrier in free energy, which can often only be overcome by sophisticated advanced sampling methods.[31] Albeit often compared with protein folding/denaturation processes, the conformational collapse of TSPs has often been studied with comparably short simulation time.[32,33] We schematically depict the CGT in Figure .

Figure 1

Schematic illustration of the coil–globule transition (CGT). Exemplary structures of coil and globule conformations of NIPAAm are shown.

To explain the thermosensitivity of the CGT, the free energy of the process needs to be evaluated at different temperatures.[31] Generally, the free energy can be separated into enthalpic and entropic contributions. Whereas the entropic contributions to the free energy can be a challenge to obtain, the enthalpic contributions can be estimated by calculating the internal energy difference between the two states, i.e., C and G. To our knowledge, these contributions to the free energy of the CGT have not been assessed in detail before. In this study, we aim to systematically compare C and G conformations at different temperatures with respect to their difference in enthalpy. Furthermore, we performed a detailed analysis of the hydration of the C and G subensembles since we expect the interactions with water to be of crucial importance for the energy balance of the transition. Moreover, to facilitate the understanding of the thermosensitive character of the transition, we investigated a set of polymers, which are partly TSP and partly non-TSP (see below).

To compare the two conformational states, i.e., C and G, at different temperatures, we separated the conformational ensembles accordingly on the basis of established structural descriptors. Since the C and G ensemble are defined by the size distribution of the polymer chains in experiments, we calculate the radius of gyration (Rg) and the hydrodynamic radius (Rh) of the conformational ensembles, which we obtained from our simulations at different temperatures. Furthermore, we calculated the solvent accessible surface area (σ), which has been shown to be a useful additional indicator to distinguish C and G.[31,34] In addition, we calculated the persistence length of the polymers, as a measure for the general stiffness of the polymer chains.[35,36]

Furthermore, we rigorously quantified the contributions to the enthalpy of the CGT process, which originate from polymer–solvent, solvent–solvent, and internal interactions of the polymer, respectively. Moreover, we calculated the lifetimes of hydrogen bonds between the polymer and solvent. A comparison of these lifetimes is useful to evaluate the hydration entropy qualitatively. As already mentioned before, to obtain reasonable estimates for the mean internal energies, sufficient sampling is crucial. Therefore, we invested a significantly longer simulation time than previous studies to minimize the bias introduced by arbitrarily chosen starting structures.

To better understand the onset of the thermosensitivity of the CGT, we compared a set of chemically closely related polymers. Despite their close relationship, they differ in their thermosensitive behavior, one half being TSPs and the other half being non-TSPs. Comparing more than one representative of both TSP and non-TSP opens up the possibility to make more general statements about the differences between these two classes of polymers. To interpret the results of simulations of TSPs consistently, it is crucial to validate the force field with non-TSPs also. Furthermore, examination of additional TSPs next to NIPAAm, a prominent model system for the thermosensitive CGT,[22,37] is highly beneficial to determine the decisive mechanisms for the thermosensitivity. To our knowledge, TSPs and non-TSPs have very rarely been compared systematically before.[38]

Model Systems

We selected four acrylamide-based polymers as model systems: besides poly-N-isopropylacrylamide (NIPAAm), we simulated three additional closely related polymers. We systematically decreased the length of the substituent at the nitrogen, obtaining poly-N-ethylacrylamide (NEAAm), poly-N-methylacrylamid (NMAAm), and polyacrylamide (AAm). These polymers are displayed in Figure . NEAAm and NIPAAm are known to undergo a thermosensitive transition with a LCST of T* ≈ 347 K,[4,22] and T* ≈ 305 K,[1,18,21] respectively. AAm on the other hand is an example for a non-TSP in watery solution.[39] NMAAm has been conjectured to exhibit an LCST above 373 K, in aqueous solution.[4,22] To our knowledge, this assertion has hitherto not been confirmed experimentally. Therefore, we expect it to be nonthermosensitive in the investigated range of temperatures.

Figure 2

Set of closely related acrylamide-based polymers. We introduce colors for easier identification in figures below. From left to right: in blue, polyacrylamide (AAm); in orange, poly-N-methylacrylamide (NMAAm); in green, poly-N-ethylacrylamide (NEAAm); and in red, poly-N-isopropylacrylamide (NIPAAm).

Generally, the existence of the CGT and the actual transition temperature of linear polymer chains also depend on the chain length.[40−42] We chose to simulate polymer chains with 30 monomer units since the 30mer is a well-established model system for the CGT of NIPAAm.[41,43,44]

Computational Methods

Simulation Setup

As starting structures for the MD simulations, we prepared extended conformations of syntactic 30mers of AAm, NMAAm, NEAAm, and NIPAAm, making use of the Maestro software package.[45] We solvated these structures in cubic boxes with a side length of 8 nm with SPC/E water.[46,47] Prior to the MD simulations, we minimized the initial configurations using the steepest descent method. Before the production runs, we equilibrated the system in short NVT simulations. Except for the preparation of the initial polymer configuration, we used the GROMACS MD-simulation software package.[48] For all simulations, we used the OPLS2005 force field,[49,50] which has been established for simulations of NIPAAm in previous publications.[31,43,44,51−53] We show force-field parameters in the Supporting Information. In our production runs, we applied the Parinello–Rahman barostat,[54,55] with a reference pressure of 1 bar and the velocity-rescaling thermostat[56] at respective simulation temperatures ranging from 250 to 360 K. We used the LINCS algorithm to constrain the bonds involving hydrogen atoms and used a timestep of 2 fs for our MD integration.[57,58] All production runs were performed with a simulation length of 2 μs. Throughout all simulations, we applied periodic boundary conditions and used the particle mesh Ewald method for treating long-range electrostatic interactions.[59]

Conformational Analysis

We calculated the radius of gyration (Rg) and the solvent accessible surface area (σ) of conformations sampled every 20 ps in our trajectory. For these analyses, we used the GROMACS tools.[60,61] Following our previously published approach, we defined the conformational subensembles of C and G in the joint Rg-σ plane.[31] Thus, we separated populations in the two-dimensional histograms in this space. We verified this conformational distinction by visual inspection of the structures. Furthermore, we validated the ratio of the mean Rg within both subensembles. This method led to a temperature-independent conformational criterion for every polymer to identify C or G, respectively. Therewith, we were able to separate C and G conformations at all applicable temperatures, which facilitated the systematic comparison of these two conformational states.

As an additional metric for the dimensions of the polymers, we calculated the hydrodynamic radius (Rh) of snapshots along our trajectory. Therefore, we used HYDROpro software to process conformations, which we sampled every 10 ns.[62] By calculating the ratio of Rg over Rh,

we were able to evaluate the compactness of the conformations of the polymer chains. Typically, we find ρ ≈ 0.7 for globular structures, whereas ρ > 1 for extended conformations, i.e., the C.[63−66]

Furthermore, to quantify the polymer stiffness, we calculated the persistence length. An introduction into the methodology, as well as all related results and discussion, can be found in the Supporting Information.

Hydration Analysis

To quantify polymer–polymer and polymer–water interactions, we performed a detailed analysis of hydrogen bonds. Therefore, we considered oxygen atoms of the polymer to be hydrogen bond acceptors and nitrogen atoms of the polymers and oxygen atoms of water molecules to be able to accept and donate hydrogen bonds. We applied a combined criterion for hydrogen bonds: first, donor and acceptor atoms need to be 3.5 Å or closer, and second, the angle between the atoms of hydrogen-donor–acceptor needs to be 30° or smaller. We separately counted the hydrogens between the following pairs of sets of atoms: polymer–polymer, polymer–water, polymer oxygen–water, and polymer nitrogen–water. To compare the polymers with respect to their interaction with water, we calculated the mean number of hydrogen bonds between the polymers and water (ν). For this analysis, we separated the structural subensembles of C and G (see above) and calculated the mean within these subensembles for all polymers at every applicable temperature. As a measure for the affinity of the surface of the polymers to water, we calculated the number of hydrogen bonds per solvent accessible surface area νσ = ν/σ. Likewise, we calculated these quantities for C and G separately.

Furthermore, we estimated the lifetime of the hydrogen bonds of the polymers with water from the autocorrelation of the existence function of the hydrogen bonds.[67] Therefore, we separately analyzed the hydrogen bonds with the nitrogen of the polymer and the oxygen of the polymer and water at different temperatures. To judge the lifetimes of the hydrogen bonds at different temperatures, we performed analogous analyses of pure water simulations and used the lifetime of the hydrogen bonds between water molecules in bulk water at different temperatures as reference.

Thermodynamic Analysis

To quantify the enthalpy of the CGT at different temperatures, we separately calculated the following contributionswhere ΔCGHSol is the difference in the enthalpy of solvent–solvent interactions, ΔCGHPol is the difference in enthalpy of the polymer, and ΔCGHPol – Sol is the difference in enthalpy of polymer–solvent interactions. All these differences ΔCG are calculated as mean changes with the conformational transition, at the respective temperature. Under the assumption that no work is done (ΔCGV = 0), we estimate the enthalpy of the CGT by the internal energy. Therefore, we calculated the enthalpy of conformations as the sum of the nonbonded potential energy terms of the respective groups, i.e., Lennard–Jones and electrostatic potentials. These potential energies have been obtained with the energy-group feature of the rerun functionality of GROMACS. In the framework of this analysis, we also took the change in the potential energy of the torsional degrees of freedom of the polymer into account. Since we ascertained that the contribution of this term to the enthalpy difference is two magnitudes smaller than the other contributions, we neglected it in the following.

Since the solvent accessible surface area decreases with the CGT, the number of water molecules bound to G is expected to be smaller than for C. Therefore, the number of polymer–solvent interactions is expected to decrease, with the collapse of the polymer chain. At the same time, the number of solvent–solvent interactions is expected to increase. Furthermore, the number of polymer–polymer interactions is expected to increase with the CGT. Assuming the number of interaction sites of solvent and polymers to be conserved, we reason that for every two water molecules, which unbind from the polymer, a water–water and a polymer–polymer bond are formed. Furthermore, assuming all bonds to be energetically equivalent (solvent–solvent, polymer–solvent, and polymer–polymer), in a crude approximationfollows. Thus, after a small rearrangement we obtain

These assumptions may be very inaccurate for other compounds, depending on the affinity of the accessible chemical groups to water. We will test the validity of these assumptions based on our simulation data.

Results

Conformational Analysis and Thermosensitivity

To facilitate the comparison of C and G, we separated our conformational ensembles at different temperatures accordingly. Therefore, we projected all sampled polymer conformations onto the Rg-σ plane for each polymer species (see Figure S5). Since we aimed to achieve a temperature-independent classification, we used the conformations from simulations at all temperatures to make our state definition. In these two-dimensional histograms, we separated the conformational subensembles C and G for every polymer in the series. We note that this separation becomes less distinct with decreasing side chain length. In particular, for the polymers, which do not exhibit a thermosensitive CGT, i.e., AAm and NMAAm, the C ensemble is less pronounced. For comparison, the criteria for the state definition are shown in Figure S6. Furthermore, we show the defined state borders for every polymer and the mean properties within the C and G subensembles in Figures S7 and S8, respectively. There, we can observe that the state borders follow a consistent trend and that Rg(coil) over Rg(globule) is between 1.2 and 1.4 for all polymers. Additionally, we note that the latter relation appears to be systematically lower for non-TSP than for TSP.

With this state definition, we identified the predominant conformational state in the simulations at different temperatures. In Figure , we show the mean conformational state at different temperatures for every polymer. Here, the eventual temperature sensitivity is visible. It is apparent that only NEAAm and NIPAAm show a thermosensitive transition. From the turning point of the logistic fit of the mean state, we obtain an estimate for the CGT temperature. According to these fits, we report T* = 281 K for NEAAm and T* = 272 K for NIPAAm, respectively. We show the underlying timeseries of Rg, σ, and the conformational state we derived from these two quantities, for all temperatures and all polymers in Figures S1–S4.

Figure 3

Mean conformational state at different simulation temperatures (T), calculated with the state definition in two-dimensional Rg-σ histograms (Figure S5). Here, we show the arithmetic mean of the assigned states over all frames in our 2 μs simulations. The error bars represent the standard error of the mean. From top to bottom, in respective colors and symbols: in blue triangles, we show data points from simulations of AAm; in orange squares, NMAAm; in green diamonds, NEAAm; and in red crosses, NIPAAm. Where applicable, we show a logistic fit of the data and depict the turning point of the curve with a dashed line, which gives an estimate for the CGT temperature (T*).

To validate the definition of conformational states, we evaluated the compactness of the polymer conformations. Therefore, we characterized the polymers’ ratio of radius of gyration over hydrodynamic radius, ρ = Rg / Rh, Figure S9. There, we can generally not only judge the conformational change with the temperature but also note that ρ spreads around 0.7 ± 0.05 for the collapsed conformations of all polymers in our set.

In Figure S12, we show the mean number of internal hydrogen bonds at different temperatures for the series of polymers. We observed consistent trends as for the state definition with other structural descriptors. Generally, hardly any internal hydrogen bonds exist in the C. These internal hydrogen bonds are formed when the polymer chains collapse, which is in line with the picture of the CGT outlined in earlier sections. The number of internal hydrogen bonds in the G state is higher for AAm than for the other polymers due to the primary amide function. This major difference between AAm and the other polymers influences basically all other analyses.

Hydration and Enthalpic Differences of Coils and Globules

Above, we visualize the hydrogen bonds between the polymers and water at different temperatures. We show the number of hydrogen bonds (ν) in Figure and the number of hydrogen bonds per solvent accessible surface area (νσ) in Figure , respectively. At every temperature, we analyze the C and G subensembles separately. We note that the C exhibits consistently higher ν but a lower νσ than the G. This holds true for every polymer at all temperatures. We note that νσ is smaller for the TSPs in our set, namely, NIPAAm and NEAAm, than for the non-TSPs. As a trend, νσ increases with decreasing sidechain length. Furthermore, AAm generally exhibits higher ν than the other polymers in the set, hence also much higher νσ. This is due to the fact that the primary amine can form two hydrogen bonds (as already mentioned before). Except for AAm, ν is very similar within the set. In Figure S14, we show the naïve mean over the whole trajectories, without separation of the subensembles, which is generally less straightforward to interpret. This is due to the fact that the naïve mean represents a weighted average of C and G for which the weights change with temperature.

Figure 4

Number of hydrogen bonds (ν) between polymers and water. We show data related to the respective polymers as follows: (a) AAm, (b) NMAAm, (c) NEAAm, and (d) NIPAAm. We calculated the mean and standard deviation of these quantities within the C (purple spheres) and the G (dark red diamonds), respectively, and performed a linear fit for these points. The opacity of the symbols represents the state predominantly existing at a given temperature. In addition, we plot small dots to make it easier to identify the position of the highly transparent points.

Figure 5

Number of hydrogen bonds between polymers and water per solvent accessible surface area (νσ). We show data related to these respective polymers as follows: (a) AAm, (b) NMAAm, (c) NEAAm, and (d) NIPAAm. We calculate the mean and standard deviation of these quantities within the C (purple spheres) and the G (dark red diamonds), respectively, and perform a linear fit for these points. The opacity of the symbols represents the state predominantly existing at a given temperature. In addition, we plot small dots to make it easier to identify the position of the highly transparent points. Please note the different y-axis scale for AAm.

In Figure S15, we show the mean number of hydrogen bonds between water and nitrogen or oxygen of the polymers, respectively. There, we observe no significant difference between the number of hydrogen bonds at the nitrogen and oxygen of the polymer with respect to the comparison of the C and G. This observation holds true for almost all polymers at all temperatures, with AAm being the exception (see above).

To compare the enthalpic contributions with the CGT in accordance with eq , we plotted the enthalpy difference of the polymer ΔCGHPol against the enthalpy difference of the water–water interactions ΔCGHSol, Figure . There, we generally notice that the data scatters around the diagonal,

Figure 6

Comparison of enthalpic contributions to the free energy of the CGT. ΔCGHPol is the difference in internal energy of the polymer; ΔCGHSol is the difference in internal energy of the water–water interactions. We show data related to the respective polymers as follows: (a) AAm, (b) NMAAm, (c) NEAAm, and (d) NIPAAm. In all panels, we plot the mean difference of these two quantities at different temperatures. The corresponding simulation temperature is color coded in accordance with the color bar. To facilitate the comparison of the magnitude of these two properties, we plot x = y as a dashed grey line.

x = y. From that, we conclude that the approximation ΔCGHSol ≈ ΔCGHPol is valid within a certain variance (eq ). Thus, both these quantities are of very similar magnitude at all temperatures for all polymers. Taking a closer look, we note that the shorter the sidechain, the more the trendline is shifted in the positive direction (linear fit in green in the respective panels). Accordingly, as a trend, ΔCGHSol > ΔCGHPol for these polymers. Generally, since both ΔCGHPol and ΔCGHSol are negative, the G is favored at all temperatures for all polymers in respect to these quantities. Once more, both these quantities are lower for AAm in comparison to the other polymers because it can form an additional hydrogen bond per monomer unit (see above).

In Figure , we compare ΔCGHPol – Sol with the sum ΔCGHPol + ΔCGHSol. Here, we can see that

Figure 7

Comparison of the enthalpic contributions to the free energy of the CGT. ΔCGHPol is the difference in internal energy of the polymer, ΔCGHSol is the difference in internal energy of the water–water interactions, and ΔCGHPol – Sol is the difference in internal energy of the polymer–water interactions. We show data related to the respective polymers as follows: (a) AAm, (b) NMAAm, (c) NEAAm, and (d) NIPAAm. In all panels, we plot the sum, ΔCGHPol + ΔCGHSol vs ΔCGHPol – Sol at different temperatures. The corresponding simulation temperature is color coded in accordance with the color bar. To facilitate the comparison of the magnitude of the two axes, we plot y = −x as a dashed grey line within the respective panels. Furthermore, we mark the area where x + y > 0 in grey. In this region, the total enthalpy of the transition (ΔCGHTot = ΔCGHSol + ΔCGHPol + ΔCGHPol – Sol) would be unfavorable. In addition, we plot a linear fit of the data with the function: y = −x + b in green.

ΔCGHPol – Sol > 0 for all polymers at all temperatures. Therefore, the enthalpy of the polymer–solvent interactions generally disfavors the CGT, which agrees with the expectation. In accordance with eq , we plot the linear function y = – x in a dashed grey line, as an orientation. In addition, we fit

the data of all polymers with the linear function y = – x + b. The resulting fits confirm that

ΔCGHPol + ΔCGHSol ≈ – ΔCGHPol – Sol is generally a valid approximation (eq ). We note that the shorter the sidechain, the further is the trendline shifted in negative direction. We show the differences in ΔCGHTot in Figure S18. There, we can see that the mean ΔCGHTot shows a clear trend to be higher for longer sidechains. While on average ΔCGHTot > 0 for NIPAAm, for AAm, almost all values are below 0.

Furthermore, we compared ΔCGHPol and ΔCGHSol with −ΔCGHPol – Sol/2 in Figures S16 and S17, respectively. According to eq , these three quantities are expected to be approximately equal. In Figure S16, we can ascertain that as a trend ΔCGHSol and −ΔCGHPol – Sol/2 are indeed approximately equal for all polymers, whereas the shorter the sidechain, the systematically lower is ΔCGHPol in comparsion to −ΔCGHPol – Sol/2.

In Figure , we depict the mean lifetimes of hydrogen bonds between the polymers and water. We display the ratio of the lifetimes of hydrogen bonds between polymer and water (⟨τPol – Sol ⟩) over the lifetimes of hydrogen bonds in bulk water at the same temperature (⟨τSol – Sol ⟩). The lifetimes of the hydrogen bonds with oxygen and nitrogen are displayed in separate panels. Generally, ⟨τPol – Sol ⟩ is longer with increasing sidechain length. This trend is more prominent for hydrogen bonds involving the polymer’s oxygen atoms of the amide groups than for the hydrogen bonds at the nitrogen atoms. Nevertheless, the trend is also visible there. We note that at low temperatures, ⟨τPol – Sol ⟩ is significantly longer than ⟨τSol – Sol ⟩, for both nitrogen and oxygen. Moreover, the higher the temperature, the smaller is the ratio of ⟨τPol – Sol ⟩ and ⟨τSol – Sol ⟩. In Figure S19, we show the unscaled values of ⟨τPol – Sol ⟩, and as a reference, we show ⟨τSol – Sol ⟩ in bulk solvent in Figure S20.

Figure 8

Mean lifetimes of hydrogen bonds between the solvent and amide group of the different polymers (⟨τPol – Sol ⟩) at different temperatures. Scaled by the mean lifetime of hydrogen bonds in bulk water at the same temperature (⟨τSol – Sol ⟩). In the upper panel, we show the lifetimes between the solvent and nitrogen atoms of the polymers, and in the lower panel, we show the lifetimes between the solvent and oxygen atoms of the polymers. In blue triangles, we show data points from simulations of AAm; in orange squares, NMAAm; in green diamonds, NEAAm; and in red crosses, NIPAAm. We show a logistic fit of the data in the respective color.

Discussion

Conformational States and Thermosensitivity

Our estimates of the CGT temperatures are in excellent qualitative agreement with experimental references, i.e., they reproduce the existence and ordering of CGT temperatures within the set of polymers. In experiments, the CGT temperatures of NIPAAm and NEAAm have been determined as approximately 305[1,18,21] and 347 K,[4,22] respectively. We estimate these CGT temperatures to be 272 and 281 K, respectively. Therefore, despite being qualitatively correct, the obtained transition temperatures are systematically lower than the experimentally obtained values. We believe that this shift originates from small inaccuracies of the force field and water model (see below). However, our estimate for the CGT temperature of NIPAAm is in perfect agreement with prior computational publications. Generally, this temperature has been estimated to lie within 270–280 K with this combination of the force field and water model.[31,52] Hence, considering our superior sampling efforts combined with our improved temperature resolution, we are confident to estimate the CGT temperature with higher accuracy than before.

For the first time, a set of closely related polymers has been simulated for such a long simulation time, i.e., 2 μs, at different temperatures. Due to the high barrier between C and G, the expected timescales of the transition are generally long.[31] In Figure S4, we show the timeseries of Rg, σ, and the resulting conformational state of NIPAAm at different temperatures. There, we can assess that, at 270 K, the first transition to a stable globule only occurs after almost 1 μs. The accuracy of MD-based predictions inherently relies on the extent of the captured conformational ensemble. Insufficient sampling may provide unphysical state populations, which inevitably result in erroneous predictions and unreliable observations. Therefore, simulations of only a few 100 ns may not be enough to determine the preferred conformational state at low temperatures. Nevertheless, we are aware that also 2 μs of simulation time will most likely not suffice to obtain convergence in all degrees of freedom. However, we do observe back and forth transitions between the conformational subensembles. This indicates, that we, in fact, sample close-to-equilibrium distributions and not only kinetically trapped conformations.

We further report that the previously published method for separating C and G in the Rg-σ plane is applicable to other polymers besides NIPAAm. Despite this method being increasingly challenging with less pronounced populated C states, we were able to define distinct state borders, leading to thoroughly consistent results in all aspects. As an additional verification, we evaluated the ratio of the mean Rg(coil) over Rg(globule), which agrees with the expected values and published size distributions for NIPAAm 30mers.[68] Furthermore, having analyzed the ratio of Rg over Rh, i.e., ρ, we can consistently confirm the globularity of the structures at the respective temperatures. In summary, we are confident in the validity of the separation of the conformational subensembles. This method enabled us to systematically compare C and G at different temperatures.

Based on the aforementioned state separation, we characterized the hydration of the C and G subensembles of the polymers at different temperatures. Therefore, we profiled the interaction of each polymer with the surrounding water molecules. We report that our results for ν are not only well in line, with similar computational studies on NIPAAm,[69] but also qualitatively agree with experimental measurements of the hydration number of NIPAAm.[70,71] Comparing ν and νσ, we notice that the comparison between C and G gives qualitatively consistent results for all polymers in our set at all temperatures, no matter whether thermosensitive or not. Generally, ν is higher in C than in G. This is due to the fact that the polymer forms internal hydrogen bonds with the structural collapse. In contrast, νσ is generally lower in C than in G. Thus, the general affinity to water of the surface changes. In addition, the volume of the solvation shell, therefore, the number of nonbulk waters, decreases with the conformational transition to the G. Furthermore, we find that νσ systematically decreases with increasing sidechain length. Generally, these results go in line, with the quantification of the enthalpic contributions of the CGT.

We were able to confirm our initial assumptions: On the one hand ΔCGHSol ≈ ΔCGHPol (eq ) and on the other hand ΔCGHSol + ΔCGHPol ≈ – ΔCGHPol – Sol (eq ). Since these equations are approximately true, we reason that the general picture of the release of water molecules from the solvation shell to bulk water, with polymer bonds forming is valid. Nevertheless, we noticed systematic deviations from these approximations. Comparing the values for ΔCGHtot of the different polymers, we find the following trend: With increasing sidechain length, on average the globule becomes less favored. Furthermore, according to our results, the G is not generally enthalpically favored. Additionally, the enthalpy of the CGT is dominated by ΔCGHPol – Sol, which is generally large in magnitude and positive thus unfavorable. As a trend, with decreasing size of the substituent ΔCGHSol is more positive than ΔCGHPol. Under the assumption that the number of binding sites, of both solvent and polymer molecules, is conserved, we infer that the potential energy per bond varies. Furthermore, we found that ΔCGHPol is systematically shifted in comparison to −ΔCGHPol – Sol/2. In contrast to that ΔCGHSol ≈ – ΔCGHPol – Sol/2, for all polymers. Therefore, we assume that solvent–solvent bonds and polymer–solvent bonds are approximately equivalent in potential energy for all polymers. However, the shorter the sidechain, the more favorable are polymer–polymer interactions in comparison to polymer–solvent interactions.

In summary, in terms of the potential energy, both polymer–polymer and solvent–solvent interactions favor the CGT. In contrast, the polymer–solvent interactions, which dominate the enthalpy, disfavor the CGT. Furthermore, the longer the sidechain, the weaker are the polymer–polymer interactions in comparison to the polymer–solvent interactions. Therefore, with increasing sidechain, the globule is in sum enthalpically less favorable. We assume electrostatic interactions to be of major importance for the potential energy. Therefore, deficiencies in the determination of the partial charges of the polymers may affect the energetic balance of the CGT. This may lead to shifted thermosensitive behavior, as already reported by the referenced studies.[72,73] We furthermore believe, that the strength of the polymer–polymer and polymer–solvent bonds may be influenced by sterical hindrance, i.e., binding partners with high affinity may not get into close contact for geometrical reasons. This would not be contradicted by the known effect of the tacticity on the thermosensitivity.

Lastly, we evaluated the lifetime of hydrogen bonds between the polymers and solvent. Generally, ⟨τPol – Sol ⟩ is longer than ⟨τSol – Sol ⟩ for all polymers in our set at every temperature we simulated. With increasing temperature, the ratio of ⟨τPol – Sol ⟩ over ⟨τSol – Sol ⟩ decreases. Furthermore, we note that this ratio is systematically higher with longer sidechains. From significantly longer lifetimes of hydrogen bonds, we reason an increased immobility of water molecules at the surface of the polymers, which is consistent with recently published experimental results.[74] We hypothesize that the trend of increasing ⟨τPol – Sol ⟩ with increasing sidechain is partly due to steric hindrance of water molecules, which intercalate between neighboring oxygen and nitrogen atoms. Therefore, we expect the sterical barrier for water molecules in the solvation shell to diffuse to the bulk to grow with the increasing size of the group at the nitrogen. Moreover, we hypothesize that large rearrangements in the solvation shell are necessary for the CGT to occur, which agrees with recent studies.[53,75] Therefore, the deceleration of the water dynamics at low temperatures significantly increases transition times.

Furthermore, we believe that the immobilization of water molecules in the solvation shell translates into an entropic penalty in comparison to bulk water. Therefore, as it has been stated before in literature, we expect the thermosensitivity to originate from an entropic effect.[76,77] With increasing temperature, the effect of this penalty on the free energy increases. Furthermore, the larger the solvation shell volume, the larger the impact of this penalty. This hypothesis is well in line with recently published experimental data[74,76] and computational studies.[78,79] Furthermore, this theory is not contradicted by the fact that the CGT temperature shifts in water–alcohol mixtures or in deuterated water.[76,80,81] Accordingly, we expect the water model in simulations to be of great importance for the balance of the CGT.[76] Hence, we surmise that the temperature shift for the CGT temperature may find its origin in the interplay between deficiencies of the water model and inaccuracies in the force field. According to our results, a purely enthalpic contemplation of the energetic balance of this transition cannot fully explain the thermosensitivity of the process.

Conclusions

The separation of the conformational ensembles at different temperatures enabled us to make some general statements about the hydration of C and G. Comparing the number of hydrogen bonds per solvent accessible surface area (νσ) for TSPs and non-TSPs, we can identify clear distinctions. We conclude that not only the affinity of the polymers’ surface to water but also the decrease in the solvation shell volume with the collapse of the chain play an important role in the stabilization of C. We want to emphasize that the differences in hydration between C and G can be identified at temperatures above and below the CGT temperature. Therefore, merely describing the hydration of the conformations, without making considerations about the effect on the free energy of the whole system (including solvent), may not lead to an exhaustive explanation of the thermosensitivity.

We want to emphasize that the non-TSPs in our set collapse at every simulated temperature, whereas the TSPs do not collapse below the transition temperature. Therefore, to explain the thermosensitivity, it is important to understand what stabilizes the C at lower temperatures, rather than what drives the collapse at higher temperatures. Furthermore, we report that it is not sufficient to simulate only a few hundred nanoseconds (or even less), especially at low temperatures, since the timescales of the transition are of a similar magnitude. Depending on the starting structure, the simulation results may be largely biased if the simulation time is too short.

According to our results, the G is not generally enthalpically favored. We found that with increasing size of the substituent at the nitrogen, the polymer–polymer interactions with forming the G become weaker in comparison to the polymer–water interactions. This leads to the G being enthalpically less favored with longer sidechains. Furthermore, we ascertained the lifetimes of the hydrogen bonds between polymer and water to be generally longer than in bulk water. In addition, we note that these lifetimes are systematically higher in TSPs than in non-TSPs. Therefore, we hypothesize that the higher affinity of water to the surface of C leads to an entropic penalty for the systems due to tightly bound waters. As a result, at higher temperatures, the C will not be stabilized by the surrounding solvent anymore.

In conclusion, we believe that the thermosensitive CGT can only be fully explained with a detailed consideration of the entropic contributions to the free energy of this transition. We believe that the entropy of the solvation is a crucial quantity for the free-energy balance of the CGT. We aim at quantifying this property in future publications.