Selective molecular transport in thermo-responsive polymer membranes: role of nanoscale hydration and fluctuations
For a wide range of modern soft functional materials the selective transport of sub-nanometer-sized molecules (‘penetrants’) through a stimuli-responsive polymeric membrane is key to the desired function. In this study, we investigate the diffusion properties of penetrants ranging from non-polar to polar molecules and ions in a matrix of collapsed Poly(N-isopropylacrylamide) (PNIPAM) polymers in water by means of extensive molecular dynamics simulations. We find that the water distributes heterogeneously in fractal-like cluster structures embedded in the nanometer-sized voids of the polymer matrix. The nano-clustered water acts as an important player in the penetrant diffusion, which proceeds via a hopping mechanism through ‘wet’ transition states: the penetrants hop from one void to another via transient water channels opened by rare but decisive polymer fluctuations. The diffusivities of the studied penetrants extend over almost five orders of magnitude and thus enable a formulation of an analytical scaling relation with a clear non-Stokesian, exponential dependence of the diffusion coefficient on the penetrant’s radius for the uncharged penetrants. Charged penetrants (ions) behave differently as they get captured in large isolated water clusters. Finally, we find large energetic activation barriers for hopping, which significantly depend on the hydration state and thereby challenge available transport theories.
KEYWORDS: polymers, thermo-responsiveness, hydration, penetrant diffusion, energy barrier, molecular dynamics simulation
Thermo-responsive polymer systems, typically realized as polymeric liquids, networks, or hydrogels in solutions, have become integral components in the engineering of soft functional materials. Their attractiveness stems primarily from their responsive behavior, high water content, and rubbery nature, being similar to biological tissue Peppas (1997). These versatile and useful features triggered many applications in material science, spanning from drug delivery Peppas (1997); Stuart et al. (2010); Kabanov and Vinogradov (2009); Oh et al. (2008); Guan and Zhang (2011), catalysis Zhang et al. (2010); Lu et al. (2009); Wu et al. (2012), biosensing Guan and Zhang (2011); Lapeyre et al. (2006), thin-film techniques Guan and Zhang (2011), as well as in environmental science Parasuraman and Serpe (2011), including nanofiltrationShannon et al. (2008), water purification, and desalination applications Nghiem et al. (2004); Geise et al. (2011). One of the most commonly studied responsive polymer is poly(N-isopropylacrylamide) (PNIPAM), which exhibits its volume transition in water close to the human body temperature Pelton (2000); Gil and Hudson (2004); Halperin et al. (2015). As a versatile model component, it has been prototypical for many developments of soft responsive materials Stuart et al. (2010); Lu et al. (2006); Ballauff and Lu (2007).
In most of the above mentioned applications, the selective diffusive transport of small solutes (‘penetrants’) through the polymer matrix is key to the desired function. For instance, in nanocarrier particles, the permeability of the polymer shell plays a crucial role in uptake and release rates of functional penetrants Stuart et al. (2010). In particular, hydrogels can be used to release small ligands and drugs in a controlled way over time, or as a response to a local chemical or physical stimulus Hoffman (1987); Gehrke et al. (1997). In desalination, selective salt flux and permeability is key to be optimized, which strongly couples to water permeability of the polymer membrane Geise et al. (2011). Other examples include responsive ‘nanoreactors’ Lu et al. (2006); Wu et al. (2012); Contreras-Cáceres et al. (2008); Horecha et al. (2014), where a responsive hydrogel shell around the catalysts is used to tune the reactor’s selective permeability of reactants and by that the reaction rate Herves et al. (2012); Roa et al. (2017).
Due to its fundamental importance, the study of the penetrant transport through polymer meshes and gels has a long history. While the diffusion in swollen (dilute or semi-dilute) polymer regimes has been tackled successfully by a large body of different theories Masaro and Zhu (1999); Amsden (1998), a much more challenging problem represents the diffusion in collapsed, highly concentrated polymers, such as melts, glasses, and gels. Indeed, due to material-specific complexity on nanoscales, such as specific solute–polymer interactions in aqueous solution Walter et al. (2010); Algaer and van der Vegt (2011); Tucker and Stevens (2012); Mukherji and Kremer (2013); Chiessi and Paradossi (2016); Rodriguez-Ropero and van der Vegt (2014); Micciulla et al. (2016); Lesch et al. (2016); Adroher-Benítez et al. (2017); Kanduč et al. (2017), the full understanding of the problem demands insights by rigorous atomistic modeling. Pioneering simulations since the early ’90s Takeuchi (1990); Müller-Plathe (1991) have shown that transport in a restricting dense polymer matrix can intrinsically differ from those in liquids and dilute systems. In fact, the diffusion of low-mass particles follows a hopping motion Müller-Plathe (1991); Müller-Plathe et al. (1992); Gusev et al. (1994); Müller-Plathe (1998a); Fritz and Hofmann (1997); Hahn et al. (1999); Hofmann et al. (2000); van der Vegt (2000); Neyertz et al. (2010), where the penetrant is constrained within a microscopic cavity for most of the time and only occasionally hops into a neighboring cavity through a transient void that opens between the chains Sok et al. (1992); Müller-Plathe et al. (1993a); Hofmann et al. (2000). Despite the aforementioned insights, the long-time diffusivity in dense, atom-resolved polymer models is often so low that its characterization poses a serious challenge for atomistic simulations even nowadays Xi et al. (2013). Hence, most useful insight on scalings (e.g., the diffusion versus temperature or solute size) have become just recently available from still quite coarse-grained simulations without solvent Zhang and Kumar (2017); Zhang et al. (2018) and generic statistical mechanics theories on activated hopping Cai et al. (2015); Zhang and Schweizer (2017).
Hence, molecular insights into penetrant diffusion in concentrated polymers including chemical specificity and in particular the influence of hydration has been very scarce. What is of particular concern in responsive polymer solutions or hydrogels is that the hydration state and polymer packing fraction are temperature dependent. As such, temperature interpolates between very wet swollen and much less hydrated collapsed states of the polymer, and it is unclear what role the water plays in the particular regimes. Is there sufficient water such that diffusion is simply Stokes-like as in a homogeneous fluid around obstructions Masaro and Zhu (1999); Amsden (1998)? Or is the diffusion dominated by hopping and the structural rearrangements of polymer segments as in the statistical mechanics pictures of very dense and ‘dry’ polymer melts Cai et al. (2015); Zhang and Schweizer (2017)? If so, does water contribute to the temperature dependence and activation (free) energies of diffusion, and how?
Experiments indicate that the distribution of water in dense polymers may be highly non-trivial and may crucially contribute to transport and thermodynamic properties of penetrant molecules Müller-Plathe (1998b). For instance, it has been indicated that, depending on the polarity of the host polymer matrix, water can either be uniformly distributed Karlsson et al. (2004) or tend to cluster Fukuda (1998). The existence of disconnected water clusters in collapsed PNIPAM gels was also revealed by dielectric dispersion spectroscopy Sasaki et al. (1999). Apart from these and other fragmented experimental facts, a bigger picture of the interplay between the polymers, water, and diffusing particles is still missing and is difficult to obtain from experiments. Namely, most experimental techniques average over the behavior of a large number of molecules, and thus can elucidate only particular facets of the entire problem. On the contrary, detailed and trustworthy atomistic simulations can offer at least qualitative insights of a broader picture, which we pursue in this work.
In this work, we employ for the first time extensive molecular dynamics simulations to examine the diffusion of small and medium-sized penetrant molecules of various kinds in a collapsed thermo-responsive PNIPAM polymer in water above its lower critical solution temperature (LCST). We first analyze sorbed water content in the collapsed polymer and its structure. After that we focus on diffusion of selected penetrant molecules in the polymer at different temperatures and water contents. The resulting diffusivities span almost over five orders of magnitude and hence enable us to extract reliable scaling laws to be compared to previous work. Here, sorbed nano-clustered water and matrix fluctuations play an important role in the diffusion by creating pathways for penetrant molecules. We finally analyze the hopping diffusion mechanisms in terms of energy barriers in relation to solute size and hydration state of the polymer.
Results and discussions
Water content in the collapsed phase
To set up a system of collapsed PNIPAM polymers (20-monomer-long atactic chains) above the LCST with the same chemical potential of water in as in bulk water, we first construct a system with two distinct phases Boţan et al. (2016); Adroher-Benítez et al. (2017): a collapsed amorphous polymer phase in one part of the box, forming a membrane, and a polymer-free water reservoir in the other, as shown in Figure 1a. We use a novel, recently introduced, OPLS-based force field Palivec et al. (2018) for the PNIPAM polymers, which better captures the thermo-responsiveness than the standard OPLS-AA Jorgensen and Tirado-Rives (1988). The simulation details and equilibration procedures are described in the Methods section.
Figure 1b shows equilibrated mass density profiles of water (blue shades) and the polymer (orange shades) at three different temperatures above the LCST. The water density in the polymer phase is drastically reduced compared with the bulk phase. Upon heating, the polymer expels even more water into bulk, which reflects its increasing hydrophobic character with rising temperature. In Figure 1c, the densities inside the polymer phase have been converted into mass water fraction, , and are plotted as a function of temperature (red circles). As seen, the mass fraction roughly linearly decreases with temperature from around 0.22 near the LCST to around 0.14 at 370 K. However, the quantitative consensus regarding the experimental PNIPAM/water phase diagram is limited due to various issues facing experimental measurements, such as differences in synthesis protocols, dependence on polymer chain length, and choices of criteria for the onset of demixing Van Durme et al. (2004); Halperin et al. (2015). In Figure 1c, for instance, we also show measurements by Van Durme et al. Van Durme et al. (2004) (triangles) for two different molecular weights of PNIPAM (18 and 78 kDa), which differ considerably in the water content, thereby demonstrating the subtlety and sensitivity of the system on various parameters. The water content in the collapsed PNIPAM phase can also be compared to cross-linked PNIPAM hydrogels above the transition temperature. The reported values for hydrogels Dong and Hoffman (1990); Sasaki et al. (1999) shown in the figure tend to be a bit higher than the MD values, with 0.3 near the LCST. It should be noted that an exact comparison to the hydrogel experiments turns out to be difficult due to possible additional effects of crosslinkers (around 1 wt% in both experiments), which make the whole network structure less flexible and likely more porous Raccis et al. (2011); Kaneko et al. (1995); Dong and Hoffman (1986). However, within these uncertainties, we conclude that the used novel PNIPAM modelPalivec et al. (2018) captures the overall features of the simulated system reasonably well. On the same plot, we also show MD results obtained with the standard OPLS-AA Jorgensen and Tirado-Rives (1988) force field of the same system by blue squares. With far lower values of than reported by the experiments, we do not find it suitable for our system.
Water structure in the collapsed polymer
From the above introduced two-phase system (Figure 1a) we now switch to a system with a collapsed polymer phase occupying the entire simulation box (Figure 2a). Here, the water amount is set in accordance with the outcomes of the two-phase model. Due to periodic boundary conditions, it mimics an infinite ‘bulk’ of collapsed PNIPAM polymers in water or a collapsed weakly cross-linked gel. The simulation details are described in the Methods section.
A typical simulation snapshot in various representations is shown in Figure 2. Panel (a) displays all atoms in the space-filling representation, with PNIPAM chains in blue and water in white–red. Panel (b) highlights the PNIPAM backbones, thereby revealing a disordered amorphous structure with voids between the chains that can accommodate water molecules, displayed in (c) and (d). As can be noted, water molecules are not uniformly distributed throughout the polymer phase. Instead, they tend to structure into heterogeneous formations, that is, into nanosized water clusters and pockets, as featured in panel (d), where individual water clusters (see further below for the definition of a cluster) are depicted in different colors.
In Figure 3a we plot the static structure factors of water oxygen atoms (blue shades) and PNIPAM heavy atoms (orange shades) that result from the simulations at different temperatures. For comparison, the structure factor of bulk (SPC/E) water at 340 K is plotted by a green dash-dotted line. The structure factor of water in the polymer differs only slightly from that of bulk water for higher wavenumbers, nm. Only the disappearance of the characteristic shoulder for water at nm indicates some distortion of the tetrahedral hydrogen-bond network Sorenson et al. (2000); Sedlmeier et al. (2011). At lower values, particularly for nm (i.e., for length scales nm), the structure factor in the polymer phase reaches much higher values than in bulk, which must be ascribed to large-scale spatial correlations due to the clustering. Here, also the structure factor of the polymer shows large peaks and reflects the mesh and void structure on the scale of 1 nm, which accommodates the water clusters. Increasing temperature mostly affects the signal at low -values and signifies changes of the nanosized cluster distribution.
For purposes of further analysis, we define a cluster as a group of the water molecules that are mutually separated by less than 0.35 nm (corresponding to the size of the first water hydration shell Berendsen et al. (1987)). Figure 3b shows the cluster-size distribution , defined as the fraction of clusters in terms of their number size , in the log–log scale at different temperatures. As can be seen, the clusters do not have a characteristic size but are extremely polydisperse. The size distribution of smaller clusters (with below 20) can be roughly described by the power law . Larger clusters (with ) tend to be progressively less present than assumed by the power law. Note a non-monotonicity of the distribution of larger clusters with rising temperature, which could be due to competing effects of varying the water content and increasing hydrophobicity of the polymer matrix.
A glance at the snapshot in Figure 2d reveals that the clusters are far from being compact structures, but rather of ‘lacy’ forms. It is known that various aggregation processes (e.g., in colloidal systems) can lead to random cluster formations describable as fractals Forrest and Witten Jr (1979); Smirnov (1990); Sorensen and Roberts (1997). A consequence of the fractal morphology is that the cluster size, typically expressed in terms of its radius of gyration , scales with the number of particles as , where is fractal dimension Smirnov (1990). Lower values of are associated with more open structures of clusters, and higher values with more compact clusters. The plot in Figure 3c, showing the evaluated mean square of versus the cluster size , reveals a clear linear relation for not too large clusters,
with nm (determined by the fit to the MD values for ), shown by a dashed line. The fractal dimension of the clusters is hence , that is, essentially the same as for the random walk. Larger clusters, however, deviate from the above scaling and tend to be more compact, namely, smaller than predicted by the scaling relation. In the limiting case of completely compact (i.e., spherical) clusters, we expect to scale as .
Molecular penetrant hydration
We examine different types of penetrant molecules, categorized into three major groups: non-polar molecules, polar molecules, and ions (shown in Figure 4). The non-polar molecules comprise noble-gas atoms, simple alkanes, and simple aromatics. The polar molecules are water, methanol, and 4-nitrophenol (NP). Finally, we treat three monatomic ions and 4-nitrophenolate (NP), which is the deprotonated form of NP. The used force fields are described in the Methods section.
We first explore the hydration environment in which the penetrants diffuse. For that, we monitor the average hydration number , defined as the number of water molecules that reside within a spherical shell of radius 0.54 nm (corresponding to the first hydration shell of CH in pure water, see Supporting Information) around any of the penetrant’s atoms. In Figure 5 we show the probability distributions of the hydration numbers for several penetrants in the collapsed polymer: The non-polar solutes (Me and NB) exhibit a rather weak hydration, with a significant probability of even a completely dry environment (). On the other hand, a polar NP is notably more hydrated than the similarly sized NB. An extreme case of the hydration behavior is exhibited by ions (shown for Cl), which are ‘wrapped’ in a considerable hydration layer that they not give up readily. From experimental conductivity measurements, it was already speculated that small ions in a collapsed PNIPAM gel get captured in isolated water clusters Sasaki et al. (1999). The above results also suggest that since ions tend to preferably reside well inside the water clusters and non-polar molecules in dry regions, the polar molecules may then preferentially reside at the boundaries of the water clusters with their non-polar parts pointing towards the polymer and the polar parts towards water. For comparison, we plot in the inset also the hydration numbers of the same molecules in bulk water. The comparison indeed reveals a much larger dehydration of non-polar molecules when transferred from bulk water into the polymer phase.
Molecular penetrant diffusion
An example of a diffusion trajectory is shown in Figure 6a (projected on the plane) for the case of NP. The trajectory exhibits individual localized states, where the particle dwells for a sufficient amount of time without considerable migration. After some time the particle suddenly performs a larger jump to a different location where it gets ‘localized’ into another dwelling state. This is a characteristics of hopping diffusion Takeuchi (1990); Müller-Plathe (1991), which has been identified also experimentally in collapsed PNIPAM-based gels Ghugare et al. (2010); Philipp et al. (2014) and is observed quite generally in simulations of amorphous melts and glassy polymer matrices Müller-Plathe (1991); Müller-Plathe et al. (1992); Gusev et al. (1994); Müller-Plathe (1998a); Fritz and Hofmann (1997); Hahn et al. (1999); Hofmann et al. (2000); Kucukpinar and Doruker (2003); Mozaffari et al. (2010).
To quantitatively characterize the hopping events, we define the displacement of the particle from the preceding mean position as
The first term in the brackets is the time-averaged position of the particle during the time window , whereas the second term is the position of the particle at time . The quantity thus measures the square displacement of the particle from its mean position during the preceding time . We plot the displacement with 100 ns for the NP trajectory [presented in panel (a)] in Figure 6b. During the dwelling states, fluctuates around small values, reflecting the degree of the particle’s fluctuations within a cavity. As the particle hops from one cavity to another, exhibits an abrupt jump, followed by a decay on the time scale of . Note that a hopping transition is much shorter than a typical residence time in the cavities. In a case the particle performs only a larger displacement fluctuation and promptly jumps back into the previous dwelling location, which is not considered as a hopping transition, exhibits a -like spike. This property enables us to distinguish hopping events from temporal displacement fluctuations. By using simple numerical criteria (see Methods), we isolate the hopping transitions from the trajectory, which are indicated by green vertical lines and annotated by letters A–D in Figure 6b, and correspondingly denoted also in panel (a).
Having identified and characterized the hopping transitions, we are now able to track the changes in the hydration, , and the polymer coordination, , numbers during the transitions. Here, counts, analogously as , the number of the polymer heavy atoms within the cut-off nm. Figure 7 shows the resulting averaged time evolution of the changes in both numbers during the hopping events. As representative examples we show cases of (a) methane (small hydrophobe), (b) nitrobenzene (large hydrophobe), (c) nitrophenol (hydrophilic solute), and chloride (charged solute) in (d). Universally, in all the cases, the penetrants get more hydrated during a hopping transition (observed jump in ), which is for larger penetrants accompanied by a ‘detachment’ from the polymer (drop in ). Both quantities relax back to their respective mean values thereafter on a time scale of around 10 ns. These main features are not qualitatively affected by the particular choice of the cut off (see Supporting Information). The results convey a clear three-phase pattern of the hopping mechanism: (i) by thermal fluctuations, a transient channel opens, creating a water pathway from one cavity to an adjacent one, (ii) the penetrant hops through the channel, (iii) the channel closes. This pattern, which appears to be quite general, has been found in other amorphous polymers as well Takeuchi (1990); Sok et al. (1992); Müller-Plathe et al. (1993b, a); Gusev and Suter (1993); Gusev et al. (1994); Fritz and Hofmann (1997); Müller-Plathe (1998b). The results imply that particles in our system universally perform hops not through dry regions of the polymer phase (even if they are non-polar), but through transient water channels instead.
We thus arrive at a central question of this study: How does the long-time diffusion coefficient depend on the particle size and type? As is often the practice in the literature Amsden (1998), we express the sizes of the penetrants in terms of their hydrodynamic Stokes radii in pure water, (see Methods). We plot the long-time self-diffusion coefficients in the collapsed polymer at 340 K versus the size of the penetrants in Figure 8. The results display a dramatic, five orders of magnitude large decrease in the diffusion coefficients as the penetrant size increases by a factor of 7. Orders of magnitude decrease of diffusion coefficients of polar molecular penetrants in a collapsed PNIPAM has been observed also in experiments Zhang et al. (2002).
The results show that the diffusion coefficients of the non-polar and polar neutral molecules roughly follow the same trend, which can be described by an exponential function
The best fit of this equation to the MD results for the neutral penetrants yields the decay length 0.019 nm. Based on these results, the behavior of the polar molecules does not differ significantly from the non-polar ones, indicating that the role of polarity in the diffusion plays a minor role. On the other hand, the diffusion of ions (shown by red triangles) clearly deviates from the trend of the neutral molecules of similar size. Notably, in the latter case, strong Coulomb interaction and the hydration seem to add an important contribution to the diffusion.
The significance of the fluctuating polymer matrix can be furthermore demonstrated by additional simulations of the same system but with rigid (i.e., completely immobilized) polymer chains Sok et al. (1992); Sentjabrskaja et al. (2016). The resulting diffusivities of selected penetrants are shown by filled purple symbols in Figure 8. In the rigid case, the diffusion plummets by an order of magnitude for the smallest penetrant helium, and even more for larger penetrants! Hence, the Boltzmann probability of crossing a barrier in the rigid matrix is far lower than in the case of mobile chains, which is additionally confirming that particle jumps are exceedingly facilitated by an elastic response of the matrix, which enables the openings of channels.
Activated hopping: Analysis of temperature and hydration effects
We now address the influences of temperature and hydration on the diffusion. Both effects are demonstrated on the case of ethane by two scenarios in Figure 9a in a form of an Arrhenius plot. In the first scenario (blue circles), the water content in the polymer is fixed at 0.14 at all temperatures. Here, the expected diffusivity increment with temperature can be entirely ascribed to thermal effects. The second scenario (orange squares) corresponds to the polymer phase with the water component in chemical equilibrium with bulk water at each temperature, where the water fraction equals the values shown in Figure 1c. Also in this case, the diffusion coefficient rises with temperature, albeit significantly less than in the first scenario. Namely, the thermal effects in this case mix with the effects of temperature-dependent (de)hydration. Hence, at given temperature, a more hydrated polymer provides higher diffusivity than a less hydrated one.
The temperature dependence of activated diffusion can be analyzed by the Arrhenius relation
Assuming that the prefactor is independent of temperature, the energy barrier then follows as
where the partial derivative is taken at constant water fraction . Figure 9b shows evaluated energy barriers for selected uncharged penetrants versus their sizes. Clearly, the results suggest a linear dependence
plotted by a dashed line with the fitting parameter kJ molnm.
In the other scenario, where the water component is in equilibrium with bulk reservoir, the resulting temperature dependence of the diffusivity is more gradual (Figure 9a), with the slope
This quantity represents an apparent energy barrier that lumps together the contributions of the thermal influence and the variable hydration into a single macroscopic parameter. The values of , shown in Figure 9b, are much less sensitive to the particle size than . Consequently, an important finding is that the variable hydration contributes a strongly compensating component to the apparent energy barrier for diffusion.
Activated hopping: Discussion
The rate-determining step in the hopping diffusion is the opening of a channel, which is associated with a free energy barrier and can be via Boltzmann probability related to the diffusion coefficient as . In conjunction with the empirically obtained diffusion relation (eq 3), this implies
This linear dependence of the free energy barrier on the particle size represents a special case of possible scenarios predicted by various theories. Most of the now classical theories that are based on activated diffusion describe regimes with either square Peppas and Reinhart (1983); Amsden (1999, 2002) or cubic Vrentas and Duda (1977); Vrentas et al. (1985); Fujita (1991); Dell and Schweizer (2013) scaling. However, these theories have been developed assuming a water-swollen network and bigger penetrants, such as enzymatic drugs or nanoparticles. A possible linear dependence of the free energy barrier has recently been theoretically envisioned in scaling theories for particle mobility in dense polymer solutions Cai et al. (2011, 2015) and in dense liquids by using a self-consistent cooperative hopping theory Zhang and Schweizer (2017). According to these studies, the actual linear scaling regime seems to depend in a very complex way on various system parameters Cai et al. (2015) that are still under debate Zhang and Kumar (2017); Zhang et al. (2018). Note also that these statistical mechanics theories based on ideal-chain models of polymers and no explicit hydration effects predict purely entropic free energy barriers for hopping, in stark contrast to the strong -dependence of diffusion observed in our simulations.
Previous atomistic simulation studies Kucukpinar and Doruker (2003); Mozaffari et al. (2010) support a diffusion relation that decays exponentially with the square of the effective penetrant size, , however, one should keep in mind the significance of a fit to given numerical data. If a span of particle sizes is not large enough, several functional forms can be fitted to the same data points within a given accuracy. In our results in Figure 8 the ratio between the smallest and the largest penetrant is around 7, which is much more than in older atomistic studies, thus providing higher significance to the relation given by eq 3.
In addition, our study points to the crucial role of the water component in the polymer. We found that, at given temperature, a more hydrated polymer provides higher diffusivity than a less hydrated one. A similar effect is known in some hydrophilic synthetic polymers (e.g., polyvinyl alcohol, used in food industry for packaging materials) that possess high barrier resistance for diffusion of various gases under dry conditions, but become drastically more permeable in humid environments Hernandez (1994); Zhang et al. (2001); Karlsson et al. (2004). The phenomenon is commonly ascribed to the sorbed water, which acts as a plasticizer and reduces the tensile strength of the polymer by softening interactions between the chains Müller-Plathe (1998b); Tamai and Tanaka (1998); Karlsson et al. (2004). This means that less hydrated polymer systems have higher free energy barriers for channel formations, and according to eq 8 this implies that they should exhibit a smaller diffusion decay length . Indeed, as we show in Supporting Information, the diffusion coefficients of selected penetrants in a less hydrated polymer decay faster with their size. This can be viewed within a simplified picture as follows: Firstly, in more hydrated polymers, the average separation between the chains is larger (seen as a decreased PNIPAM partial density in Figure 1b), which have therefore more freedom for fluctuations and thereby creating channel openings more readily. Secondly, in more hydrated cases with larger and more abundant water clusters, the separations between neighboring clusters are smaller on average, consequently the activated transient channels can be shorter.
Finally, our analysis conveyed another important message: The energy barrier is proportional to the particle size. This has an interesting consequence that larger diffusing molecules are more affected by temperature changes at fixed hydration of the polymer than smaller ones. This effect can be used as an additional means in tailoring and tuning the selectivity of a molecular transport through hydrogels by external stimuli.
We have presented a detailed all-atom simulation analysis of penetrant diffusion in a collapsed PNIPAM polymer phase in water. Water has been found to structure in fractal-like clusters sorbed between the voids made by the polymeric chains. The diffusion advances via the hopping mechanism, in which a penetrant resides for longer time in a local cavity and suddenly performs a longer jump into a neighboring cavity through a transient water channel that forms between the chains.
We found that the diffusion heavily depends on the temperature and the hydration level of the polymer. These two effects are typically coupled, since the thermo-responsive nature of PNIPAM directly impacts the affinity to water, thereby regulating its hydration through an outer water reservoir. By systematic and careful simulation approaches, we were able to separate both effects and demonstrate plasticizing effects of water on the diffusion.
Furthermore, the penetrant molecules in our study extend almost over an order of magnitude in their sizes and almost over five orders of magnitude in the resulting diffusivities. This allowed us to formulate a reliable and statistically significant relation between the diffusion coefficients and the sizes of the penetrants. We find that the diffusion of non-ionic penetrants follow a universal, exponential dependence on the size. The diffusion of ions, on the other hand, deviates from this relation, probably due to additional Coulomb interactions with the polymer chains.
The outcome of this work seriously challenges theoretical understanding of diffusion in such systems. As we have seen in our simulations, the dynamics of small and mid-sized penetrants is coupled to the structural relaxation of the polymer, which has far-reaching consequences and it complicates theoretical modeling. That means, it is not sufficient to treat a collapsed polymer with rigid obstruction models, which are very popular due to their simplicity Amsden (1998); Masaro and Zhu (1999). Moreover, the hopping cannot be explained by hydrodynamic forces. In addition, the diffusion in our system is thermally activated, which manifests in large energy barriers that scale linearly with the size of the penetrant. As we have demonstrated, the hydration effects stemming from the clustered water, together with the thermo-responsive nature of the polymer, are crucial elements that have to be included in future theoretical considerations.
We believe that our findings will stimulate new theoretical efforts into this problem. Understanding the transport mechanisms inside not only PNIPAM but also other responsive hydrogels is important for the rational design of novel materials. We are currently performing a follow-up study of solvation properties of penetrants in these systems, which will complement our understanding of penetrant interactions with collapsed hydrogels.
We utilize an atomistic model of PNIPAM polymer in the presence of explicit water. The PNIPAM chains are composed of 20 monomeric units with atactic stereochemisty (i.e., with random distribution of monomeric enantiomers along the chain). To describe water in the simulations we use the SPC/E water model Berendsen et al. (1987).
For PNIPAM polymers we first tested the standard OPLS-AA Jorgensen and Tirado-Rives (1988) force field, which is among the most popular ones for PNIPAM simulations Walter et al. (2010); Algaer and van der Vegt (2011); Tucker and Stevens (2012); Mukherji and Kremer (2013); Chiessi and Paradossi (2016); Rodriguez-Ropero and van der Vegt (2014); Adroher-Benítez et al. (2017). Since this force field did not yield satisfactory hydration results in the collapsed state (see Figure 1c) and also due to revealed issues with its thermo-responsive properties in the recent literature Kang et al. (2016); Boţan et al. (2016), we use its recent modification with recalculated partial charges by Palivec et al. Palivec et al. (2018). As has been demonstrated, the novel forcefield exhibits thermo-responsive properties of a single PNIPAM polymer much closer to experimental observations. For the neutral penetrant molecules, we use the OPLS-AA force field Jorgensen and Tirado-Rives (1988); Price et al. (2001), which keeps our model on the generic level and which sufficiently captures the hydration properties in combination with the SPC/E water model Hess and van der Vegt (2006). For the deprotonated ion NP, which is not provided within OPLS-AA, we used the partial charge parameterization ‘OPLS/QM1’ from Ref. dccclxxx(36). For the monatomic ions we employ the parameters from the Jorgensen force fields Chandrasekhar et al. (1984); McDonald et al. (1998). For comparison, we additionally use the Åqvist Åqvist (1990) force field for Na and Dang et al. Dang (1992); Rajamani et al. (2004) for I.
In the first part of the study, we assembled 48 polymer chains in a slab-like structure that extended across the and box dimensions with a finite thickness (approx. 3–4 nm) in -direction. Initially, the polymers were only loosely arranged in a slab-like assembly and solvated with water. The equilibration steered the assembly into a more compact structure, thereby making two well separated phases of water-only (supernatant) and the polymer-rich (precipitant) slab along -axis as shown in the snapshot in Figure 1a. For equilibration purposes, we performed simulated annealing where the temperature was linearly decreased by the thermostat from 450 K down to a target temperature (i.e., 310 K, 340 K, or 370 K) on a time interval of 100 ns. After that, the equilibration continued at the target temperature until the water amount in the gel phase, which we monitored, equilibrated (see Supporting Information). The necessary simulation times for equilibration depend significantly on temperature and are around 2,000 ns, 1,000 ns, and 500 ns for 310 K, 340 K, and 370 K, respectively. Notably, much longer equilibration times are needed at lower temperatures due to slower kinetics. In order to improve sampling statistics, we averaged the results over four independent sets of simulations for 310 K and 340 K, and over two sets for 370 K.
In the second part we set up single PNIPAM-phase simulations. Here, 48 atactic PNIPAM polymers in a cubic box were mixed with a certain number of water molecules. We equilibrated the systems via the annealing simulations with linearly decreasing temperatures from 900 K down to target temperatures on a time interval of 100 ns. After that, the equilibration continued for another 50–100 ns. The resulting equilibrated structures (with an example shown in Figure 2) were then used as initial configurations for all further analyses in this work.
To study diffusion, we inserted 10–15 penetrant molecules of the same kind at random positions into the equilibrated polymer structures (using 2–4 independent replicas). When simulating ions, we estimate the Bjerrum length to 5–6 nm (assuming the dielectric permittivity of a collapsed PNIPAM in water to be around 10 Sasaki et al. (1999)), which could lead to mutual influences between the ions due to strong Coulomb interactions. To reduce the effects, we restricted the number of ions in a simulation box to three, which on the other hand compromised the statistics. The net charge was compensated by applying a uniform neutralizing background charge. The simulation production runs spanned from 300 ns for the smallest penetrants (He, Ne), to around 1,000 ns for medium-sized (methane, ethane), and up to 8,000 ns for the largest ones (NB, NP) and ions. Such long simulation times for the latter penetrants were required in order to observe at least several hopping events.
The molecular dynamics (MD) simulations were performed with the GROMACS 5.1 simulation suite van der Spoel et al. (2005); Pronk et al. (2013) in the constant-pressure (NPT) ensemble, where the box sizes are independently adjusted in order to maintain the external pressure of 1 bar via Berendsen barostat Berendsen et al. (1984) with the time constant of 1 ps. The system temperature was maintained by the velocity-rescaling thermostat Bussi et al. (2007) with a time constant of 0.1 ps. The Lennard-Jones (LJ) interactions were truncated at 1.0 nm. Electrostatics was treated using Particle-Mesh-Ewald (PME) methods Darden et al. (1993); Essmann et al. (1995) with a 1.0 nm real-space cutoff.
Determination of hopping transitions
Identifying the hopping transitions from is related to the step detection problem in a noisy signal and is thus a matter of subjective criteria. Due to the noisy nature, we scanned through by defining two moving average values, a lagging average , which is the mean value of in the time window , and similarly an advancing average in the window . Obviously, a step in occurs when the ratio exhibits a local peak. We specified two criteria that have to be fulfilled in order to consider an identified peak as a hopping event. First, the local peak in should correspond to the maximal value in the time window . This assumes that at most one hopping event can occur in the given time window of length . Second, the jump should fulfill the condition , that is, the particle has to jump by more than a threshold distance . For all the analyzed penetrants in Figure 7, we chose averaging interval length 5 ns (which is considerably longer than temporal displacement fluctuations and much shorter than a typical residence time of a penetrant in a cavity) and the threshold 0.75 nm (which roughly corresponds to an estimated inter-chain distance, thus to a typical cavity size).
A conventional way to analyze the particle dynamics in computer simulations, is to evaluate mean square displacement (MSD) of the particles,
where is the position of the particle at time . When evaluating an MSD, the center-of-mass motion of the whole system should be removed from the displacement vector . However, the diffusion of particles in crowded environments often behaves anomalously at short time scales, as it is coupled to the segmental dynamics of polymers, and MSD follows a more general power-law pattern , where the scaling exponent can deviate from unity Bouchaud and Georges (1990); Müller-Plathe et al. (1992); Metzler and Klafter (2000); Sokolov (2012); Ernst et al. (2014); Metzler et al. (2014); Ghosh et al. (2015). Only at sufficiently long observation times the MSD behavior crosses over to normal (Brownian) diffusion with . Since we are interested solely in the normal (long-time) diffusion, it is important to assess the crossover time for each individual MSD. As we show in Supporting Information, the crossover time depends on the size of the particle, spanning from around 1 ns for the smallest one (helium) and up to several 100 ns for larger molecules. This analysis indicates that for the largest considered penetrant particles, trajectory lengths of several s are needed in order to properly evaluate the diffusion coefficients. As noted before Gusev et al. (1994), this stringent verification of the crossover time was often lacking in many previous evaluations of MD results. Once the crossover time has been determined, the diffusion coefficient can be calculated from the Einstein relation, viz.
which we achieve by a linear fit of the MSD in the long-time limit.
In our simulations, we also evaluated the MSD of individual polymer chains. On the time scale of the simulations, the diffusion of the polymer chains is negligible compared to the diffusion of the penetrants, which indicates that the penetrant diffusion is not related to the diffusion of the whole network.
Stokes radii in water
The dominant factor that governs the diffusion is the size of the particle Gusev et al. (1994), whereas other factors, such as the shape and polarity, are typically of secondary importance. As is often the practice in the literature Amsden (1998), we express the sizes of the penetrants by their Stokes hydrodynamic radii in pure water, , defined via the Stokes–Einstein equation,
where is the diffusion coefficient in water and the water viscosity. The Stokes radius represents a suitable measure of the effective particle size, as it captures the ‘bare’ particle size together with its hydration shell. The latter one is relevant in our case, since the penetrants, as we have seen, diffuse predominantly through transient water channels. In addition, the Stokes radius is a well-defined quantity in experiments and simulations. Therefore, we used the Stokes–Einstein equation 11 to determine Stokes hydrodynamic radii of the molecules in pure water, which required the evaluation of the and . The latter was computed with the standard procedure of transverse current correlation function Palmer (1994); Hess (2002) from an independent simulation of pure water at 340 K. The obtained value mPas agrees very well with previously reported 0.38 mPas for the SPC/E water at this temperature Markesteijn et al. (2012). The diffusion coefficients were obtained from the MSDs (eq 10) of the molecules in a water box. Since hydrodynamic interactions are long range (), they lead to effective coupling between the molecule, the solvent, and the periodic images. The evaluated apparent diffusion coefficients are therefore box-size dependent. The finite-size effects on the diffusivity can be estimated by the Ewald summation of the Oseen tensor over the periodic images, which enables the evaluation of the actual diffusion coefficient as Yeh and Hummer (2004a, b)
Here, stands for the box length and the constant stems from the summation over all periodic images.
Equation 12 is based on the first-order correction and therefore applicable only for Hasimoto (1959); Yeh and Hummer (2004b). In order to determine what is the minimal box size that enables the application of the correction given by eq 12, we performed a size-scaling analysis of for water, NB, and I (see Supporting Information). Based on the outcomes, we used box sizes of nm for the larger (aromatic) molecules and – nm for smaller molecules.
The results for Stokes radii do not significantly depend on temperature in the range 310–370 K (see Supporting Information). Some differences are observed only for He and Ne. However, we use the values of obtained at 340 K for all the particles and regard them as fixed molecular parameters.
Conflict of Interest
The authors declare no competing financial interest.
Supporting Information Available
Equilibration; hydration shell of methane; finite-size correction of diffusion in water; Stokes radii at different temperatures; short-diffusion diffusion analysis; diffusion at different temperatures and hydrations
The authors thank Richard Chudoba, Vladimir Palivec, Jan Heyda, Sebastian Milster, Yan Lu, and Matthias Ballauff for useful discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement n 646659-NANOREACTOR). The simulations were performed with resources provided by the North-German Supercomputing Alliance (HLRN).
- Peppas (1997) N. A. Peppas, Curr. Opin. Colloid Interface Sci. 2, 531 (1997).
- Stuart et al. (2010) M. A. C. Stuart, W. T. Huck, J. Genzer, M. Müller, C. Ober, M. Stamm, G. B. Sukhorukov, I. Szleifer, V. V. Tsukruk, M. Urban, F. Winnik, S. Zauscher, I. Luzinov, and S. Minko, Nat. Mater. 9, 101 (2010).
- Kabanov and Vinogradov (2009) A. Kabanov and S. Vinogradov, Angew. Chem. Int. 48, 5418 (2009).
- Oh et al. (2008) J. K. Oh, R. Drumright, D. J. Siegwart, and K. Matyjaszewski, Prog. Polym. Sci. 33, 448 (2008).
- Guan and Zhang (2011) Y. Guan and Y. Zhang, Soft Matter 7, 6375 (2011).
- Zhang et al. (2010) J.-T. Zhang, G. Wei, T. F. Keller, H. Gallagher, C. Stötzel, F. A. Müller, M. Gottschaldt, U. S. Schubert, and K. D. Jandt, Macromol. Mater. Eng. 295, 1049 (2010).
- Lu et al. (2009) Y. Lu, S. Proch, M. Schrinner, M. Drechsler, R. Kempe, and M. Ballauff, J. Mater. Chem. 19, 3955 (2009).
- Wu et al. (2012) S. Wu, J. Dzubiella, J. Kaiser, M. Drechsler, X. Guo, M. Ballauff, and Y. Lu, Angew. Chem. Int. 51, 2229 (2012).
- Lapeyre et al. (2006) V. Lapeyre, I. Gosse, S. Chevreux, , and V. Ravaine, Biomacromolecules 7, 3356 (2006).
- Parasuraman and Serpe (2011) D. Parasuraman and M. J. Serpe, ACS Appl. Mater. Interfaces 3, 4714 (2011).
- Shannon et al. (2008) M. A. Shannon, P. W. Bohn, M. Elimelech, J. G. Georgiadis, B. J. Mariñas, and A. M. Mayes, Nature 452, 301 (2008).
- Nghiem et al. (2004) L. D. Nghiem, A. I. Schäfer, and M. Elimelech, Environ. Sci. Technol. 38, 1888 (2004).
- Geise et al. (2011) G. M. Geise, H. B. Park, A. C. Sagle, B. D. Freeman, and J. E. McGrath, J. Membr. Sci. 369, 130 (2011).
- Pelton (2000) R. Pelton, Adv. Colloid Interface Sci. 85, 1 (2000).
- Gil and Hudson (2004) E. S. Gil and S. M. Hudson, Prog. Polym. Sci. 29, 1173 (2004).
- Halperin et al. (2015) A. Halperin, M. Kröger, and F. M. Winnik, Angew. Chem. Int. Ed. 54, 15342 (2015).
- Lu et al. (2006) Y. Lu, Y. Mei, M. Drechsler, and M. Ballauff, Angew. Chem. Int. 45, 813 (2006).
- Ballauff and Lu (2007) M. Ballauff and Y. Lu, Polymer 48, 1815 (2007).
- Hoffman (1987) A. S. Hoffman, J. Control. Release 6, 297 (1987).
- Gehrke et al. (1997) S. H. Gehrke, J. P. Fisher, M. Palasis, and M. E. Lund, Ann. N. Y. Acad. Sci. 831, 179 (1997).
- Contreras-Cáceres et al. (2008) R. Contreras-Cáceres, A. Sánchez-Iglesias, M. Karg, I. Pastoriza-Santos, J. Pérez-Juste, J. Pacifico, T. Hellweg, A. Fernández-Barbero, and L. M. Liz-Marzán, Adv. Mater. 20, 1666 (2008).
- Horecha et al. (2014) M. Horecha, E. Kaul, A. Horechyy, and M. Stamm, J. Mater. Chem. A 2, 7431 (2014).
- Herves et al. (2012) P. Herves, M. Perez-Lorenzo, L. M. Liz-Marzan, J. Dzubiella, Y. Lu, and M. Ballauff, Chem. Soc. Rev. 41, 5577 (2012).
- Roa et al. (2017) R. Roa, W. K. Kim, M. Kanduč, J. Dzubiella, and S. Angioletti-Uberti, ACS Catal. 7, 5604 (2017).
- Masaro and Zhu (1999) L. Masaro and X. Zhu, Prog. Polym. Sci. 24, 731 (1999).
- Amsden (1998) B. Amsden, Macromolecules 31, 8382 (1998).
- Walter et al. (2010) J. Walter, V. Ermatchkov, J. Vrabec, and H. Hasse, Fluid Phase Equilib. 296, 164 (2010).
- Algaer and van der Vegt (2011) E. A. Algaer and N. F. A. van der Vegt, J. Phys. Chem. B 115, 13781 (2011).
- Tucker and Stevens (2012) A. K. Tucker and M. J. Stevens, Macromolecules 45, 6697 (2012).
- Mukherji and Kremer (2013) D. Mukherji and K. Kremer, Macromolecules 46, 9158 (2013).
- Chiessi and Paradossi (2016) E. Chiessi and G. Paradossi, J. Phys. Chem. B 120, 3765 (2016).
- Rodriguez-Ropero and van der Vegt (2014) F. Rodriguez-Ropero and N. F. van der Vegt, J. Phys. Chem. B 118, 7327 (2014).
- Micciulla et al. (2016) S. Micciulla, J. Michalowsky, M. A. Schroer, C. Holm, R. von Klitzing, and J. Smiatek, Phys. Chem. Chem. Phys. 18, 5324 (2016).
- Lesch et al. (2016) V. Lesch, A. Heuer, C. Holm, and J. Smiatek, ChemPhysChem 17, 387 (2016).
- Adroher-Benítez et al. (2017) I. Adroher-Benítez, A. Moncho-Jordá, and G. Odriozola, J. Chem. Phys. 146, 194905 (2017).
- Kanduč et al. (2017) M. Kanduč, R. Chudoba, K. Palczynski, W. K. Kim, R. Roa, and J. Dzubiella, Phys. Chem. Chem. Phys. 19, 5906 (2017).
- Takeuchi (1990) H. Takeuchi, J. Chem. Phys. 93, 2062 (1990).
- Müller-Plathe (1991) F. Müller-Plathe, J. Chem. Phys. 94, 3192 (1991).
- Müller-Plathe et al. (1992) F. Müller-Plathe, S. C. Rogers, and W. F. van Gunsteren, Chem. Phys. Lett. 199, 237 (1992).
- Gusev et al. (1994) A. Gusev, F. Müller-Plathe, W. Van Gunsteren, and U. Suter, in Atomistic Modeling of Physical Properties (Springer, 1994) pp. 207–247.
- Müller-Plathe (1998a) F. Müller-Plathe, Ber. Bunsenges. Phys. Chem. 102, 1679 (1998a).
- Fritz and Hofmann (1997) L. Fritz and D. Hofmann, Polymer 38, 1035 (1997).
- Hahn et al. (1999) O. Hahn, D. A. Mooney, F. Müller-Plathe, and K. Kremer, J. Chem. Phys. 111, 6061 (1999).
- Hofmann et al. (2000) D. Hofmann, L. Fritz, J. Ulbrich, and D. Paul, Comput. Theor. Polym. Sci. 10, 419 (2000).
- van der Vegt (2000) N. van der Vegt, Macromolecules 33, 3153 (2000).
- Neyertz et al. (2010) S. Neyertz, D. Brown, S. Pandiyan, and N. F. van der Vegt, Macromolecules 43, 7813 (2010).
- Sok et al. (1992) R. Sok, H. Berendsen, and W. Van Gunsteren, J. Chem. Phys. 96, 4699 (1992).
- Müller-Plathe et al. (1993a) F. Müller-Plathe, L. Laaksonen, and W. F. van Gunsteren, J. Mol. Graph. 11, 118 (1993a).
- Xi et al. (2013) L. Xi, M. Shah, and B. L. Trout, J. Phys. Chem. B 117, 3634 (2013).
- Zhang and Kumar (2017) K. Zhang and S. K. Kumar, ACS Macro Lett. 6, 864 (2017).
- Zhang et al. (2018) K. Zhang, D. Meng, F. Müller-Plathe, and S. K. Kumar, Soft Matter 14, 440 (2018).
- Cai et al. (2015) L.-H. Cai, S. Panyukov, and M. Rubinstein, Macromolecules 48, 847 (2015).
- Zhang and Schweizer (2017) R. Zhang and K. S. Schweizer, J. Chem. Phys. 146, 194906 (2017).
- Müller-Plathe (1998b) F. Müller-Plathe, J. Membr. Sci. 141, 147 (1998b).
- Karlsson et al. (2004) G. E. Karlsson, U. W. Gedde, and M. S. Hedenqvist, Polymer 45, 3893 (2004).
- Fukuda (1998) M. Fukuda, J. Chem. Phys. 109, 6476 (1998).
- Sasaki et al. (1999) S. Sasaki, S. Koga, and H. Maeda, Macromolecules 32, 4619 (1999).
- Jorgensen and Tirado-Rives (1988) W. L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc. 110, 1657 (1988).
- Van Durme et al. (2004) K. Van Durme, G. Van Assche, and B. Van Mele, Macromolecules 37, 9596 (2004).
- Dong and Hoffman (1990) L.-C. Dong and A. S. Hoffman, J. Control. Release 13, 21 (1990).
- Boţan et al. (2016) V. Boţan, V. Ustach, R. Faller, and K. Leonhard, J. Phys. Chem. B 120, 3434 (2016).
- Palivec et al. (2018) V. Palivec, D. Zadrazil, and J. Heyda, arXiv:1806.05592 (2018).
- Raccis et al. (2011) R. Raccis, R. Roskamp, I. Hopp, B. Menges, K. Koynov, U. Jonas, W. Knoll, H.-J. Butt, and G. Fytas, Soft Matter 7, 7042 (2011).
- Kaneko et al. (1995) Y. Kaneko, R. Yoshida, K. Sakai, Y. Sakurai, and T. Okano, J. Membr. Sci. 101, 13 (1995).
- Dong and Hoffman (1986) L. C. Dong and A. S. Hoffman, J. Control. Release 4, 223 (1986).
- Sorenson et al. (2000) J. M. Sorenson, G. Hura, R. M. Glaeser, and T. Head-Gordon, J. Chem. Phys. 113, 9149 (2000).
- Sedlmeier et al. (2011) F. Sedlmeier, D. Horinek, and R. R. Netz, J. Am. Chem. Soc. 133, 1391 (2011).
- Berendsen et al. (1987) H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
- Forrest and Witten Jr (1979) S. Forrest and T. Witten Jr, J. Phys. A 12, L109 (1979).
- Smirnov (1990) B. M. Smirnov, Phys. Rep. 188, 1 (1990).
- Sorensen and Roberts (1997) C. M. Sorensen and G. C. Roberts, J. Colloid Interface Sci 186, 447 (1997).
- Ghugare et al. (2010) S. V. Ghugare, E. Chiessi, M. T. Telling, A. Deriu, Y. Gerelli, J. Wuttke, and G. Paradossi, J. Phys. Chem. B 114, 10285 (2010).
- Philipp et al. (2014) M. Philipp, K. Kyriakos, L. Silvi, W. Lohstroh, W. Petry, J. K. Krüger, C. M. Papadakis, and P. Müller-Buschbaum, J. Phys. Chem. B 118, 4253 (2014).
- Kucukpinar and Doruker (2003) E. Kucukpinar and P. Doruker, Polymer 44, 3607 (2003).
- Mozaffari et al. (2010) F. Mozaffari, H. Eslami, and J. Moghadasi, Polymer 51, 300 (2010).
- Müller-Plathe et al. (1993b) F. Müller-Plathe, S. C. Rogers, and W. F. van Gunsteren, J. Chem. Phys. 98, 9895 (1993b).
- Gusev and Suter (1993) A. A. Gusev and U. W. Suter, J. Chem. Phys. 99, 2228 (1993).
- Zhang et al. (2002) W. Zhang, I. Gaberman, and M. Ciszkowska, Anal. Chem. 74, 1343 (2002).
- Chandrasekhar et al. (1984) J. Chandrasekhar, D. C. Spellmeyer, and W. L. Jorgensen, J. Am. Chem. Soc. 106, 903 (1984).
- Åqvist (1990) J. Åqvist, J. Phys. Chem. 94, 8021 (1990).
- McDonald et al. (1998) N. A. McDonald, E. M. Duffy, and W. L. Jorgensen, J. Am. Chem. Soc. 120, 5104 (1998).
- Dang (1992) L. X. Dang, J. Phys. Chem. 97, 2659 (1992).
- Rajamani et al. (2004) S. Rajamani, T. Ghosh, and S. Garde, J. Phys. Chem. 120, 4457 (2004).
- Sentjabrskaja et al. (2016) T. Sentjabrskaja, E. Zaccarelli, C. De Michele, F. Sciortino, P. Tartaglia, T. Voigtmann, S. U. Egelhaaf, and M. Laurati, Nat. Commun. 7 (2016).
- Peppas and Reinhart (1983) N. A. Peppas and C. T. Reinhart, J. Membr. Sci. 15, 275 (1983).
- Amsden (1999) B. Amsden, Macromolecules 32, 874 (1999).
- Amsden (2002) B. Amsden, Polymer 43, 1623 (2002).
- Vrentas and Duda (1977) J. Vrentas and J. Duda, J. Polym. Sci. Part B Polym. Phys. 15, 403 (1977).
- Vrentas et al. (1985) J. Vrentas, J. Duda, and H.-C. Ling, J. Polym. Sci. Part B Polym. Phys. 23, 275 (1985).
- Fujita (1991) H. Fujita, Polym. J 23, 1499 (1991).
- Dell and Schweizer (2013) Z. E. Dell and K. S. Schweizer, Macromolecules 47, 405 (2013).
- Cai et al. (2011) L.-H. Cai, S. Panyukov, and M. Rubinstein, Macromolecules 44, 7853 (2011).
- Hernandez (1994) R. J. Hernandez, J. Food Eng. 22, 495 (1994).
- Zhang et al. (2001) Z. Zhang, I. J. Britt, and M. A. Tung, J. Appl. Polym. Sci. 82, 1866 (2001).
- Tamai and Tanaka (1998) Y. Tamai and H. Tanaka, Fluid Ph. Equilibria 144, 441 (1998).
- Kang et al. (2016) Y. Kang, H. Joo, and J. S. Kim, J. Phys. Chem. B 120, 13184 (2016).
- Price et al. (2001) M. L. P. Price, D. Ostrovsky, and W. L. Jorgensen, J. Comp. Chem. 22, 1340 (2001).
- Hess and van der Vegt (2006) B. Hess and N. F. van der Vegt, J. Phys. Chem. B 110, 17616 (2006).
- van der Spoel et al. (2005) D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C. Berendsen, J. Comput. Chem. 26, 1701 (2005).
- Pronk et al. (2013) S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess, and E. Lindahl, Bioinformatics 29, 845 (2013).
- Berendsen et al. (1984) H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984).
- Bussi et al. (2007) G. Bussi, D. Donadio, and M. Parrinello, J. Chem. Phys. 126, 014101 (2007).
- Darden et al. (1993) T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993).
- Essmann et al. (1995) U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995).
- Bouchaud and Georges (1990) J.-P. Bouchaud and A. Georges, Phys. Rep. 195, 127 (1990).
- Metzler and Klafter (2000) R. Metzler and J. Klafter, Phys. Rep. 339, 1 (2000).
- Sokolov (2012) I. M. Sokolov, Soft Matter 8, 9043 (2012).
- Ernst et al. (2014) D. Ernst, J. Köhler, and M. Weiss, Phys. Chem. Chem. Phys. 16, 7686 (2014).
- Metzler et al. (2014) R. Metzler, J.-H. Jeon, A. G. Cherstvy, and E. Barkai, Phys. Chem. Chem. Phys. 16, 24128 (2014).
- Ghosh et al. (2015) S. K. Ghosh, A. G. Cherstvy, and R. Metzler, Phys. Chem. Chem. Phys. 17, 1847 (2015).
- Palmer (1994) B. J. Palmer, Phys. Rev. E 49, 359 (1994).
- Hess (2002) B. Hess, J. Chem. Phys. 116, 209 (2002).
- Markesteijn et al. (2012) A. Markesteijn, R. Hartkamp, S. Luding, and J. Westerweel, J. Chem. Phys. 136, 134104 (2012).
- Yeh and Hummer (2004a) I.-C. Yeh and G. Hummer, Biophys. J. 86, 681 (2004a).
- Yeh and Hummer (2004b) I.-C. Yeh and G. Hummer, J. Phys. Chem. B 108, 15873 (2004b).
- Hasimoto (1959) H. Hasimoto, J. Fluid Mech. 5, 317 (1959).