Aqueous Nanoclusters Govern Ion Partitioning in Dense Polymer Membranes.

The uptake and sorption of charged molecules by responsive polymer membranes and hydrogels in aqueous solutions is of key importance for the development of soft functional materials. Here, we investigate the partitioning of simple monatomic (Na+, K+, Cs+, Cl-, I-) and one molecular ion (4-nitrophenolate; NP-) within a dense, electroneutral poly(N-isopropylacrylamide) membrane using explicit-water molecular dynamics simulations. Inside the predominantly hydrophobic environment, water distributes in a network of polydisperse water nanoclusters. The average cluster size determines the mean electrostatic self-energy of the simple ions, which preferably reside deeply inside them; we therefore find substantially larger partition ratios K ≃10-1 than expected from a simple Born picture using a uniform dielectric constant. Despite their irregular shapes, we observe that the water clusters possess a universal negative electrostatic potential with respect to their surroundings, as is known for aqueous liquid-vapor interfaces. This potential, which we find concealed in cases of symmetric monatomic salts, can dramatically impact the transfer free energies of larger charged molecules because of their weak hydration and increased affinity to interfaces. Consequently, and in stark contrast to the simple ions, the molecular ion NP- can have a partition ratio much larger than unity, K ≃10-30 (depending on the cation type) or even 103 in excess of monovalent salt, which explains recent observations of enhanced reaction kinetics of NP- reduction catalyzed within dense polymer networks. These results also suggest that ionizing a molecule can even enhance the partitioning in a collapsed, rather hydrophobic gel, which strongly challenges the traditional simplistic reasoning.

The principal factor in the design of polymer membranes and hydrogels used in applications for water purification, pervaporation, drug delivery, nanocarriers, etc. is the capability to control the uptake and permeability of different molecular species.[1−9] For instance, the ability to control the permeability of ions makes such materials interesting for water purification technologies, such as reverse and forward osmosis.[10,11] A separation is achieved between different penetrants because of the distinction in the amount of material that dissolves in the membrane and the rate at which the material diffuses through the membrane.[2,12] A recent modification of forward osmosis, where the hydrogel acts as the draw solute (osmotic agent) and the separation membrane at the same time,[13−15] exploits the difference of the salt concentration inside the hydrogel in its swollen and collapsed (or compressed) state and can be used to yield low-salinity water upon repeating the cycle.[13]

As an alternative to mechanical compression, thermoresponsive shrinking of a hydrogel can be utilized for desalination purposes.[12,16,17] Thermoresponsive polymers, such as poly(N-isopropylacrylamide) (PNIPAM), undergo a transition from good to poor solvent conditions upon an increase in temperature or triggered by some other environmental stimulus, which results in substantial shrinkage of the gel.[18−21] An advantage of responsive hydrogels in purification lies in simple separation by filtration methods and their reusability. Thus, polymer membranes are expected to continue to dominate the water purification technologies owing to their energy efficiency.[12,22]

Hydrogels are also particularly appealing materials for drug delivery systems as they can be designed in numerous ways in order to selectively encapsulate and release particular types of pharmacologically relevant molecules in a controllable manner.[8] For charged drugs, the electrostatic interactions, facilitated by accompanying ions, are an important, if not dominating, driving force for the transport through the gel.[23] Related phenomena are exploited in adaptive nanocatalysis where catalytic nanoparticles are confined in a permeable hydrogel that can be potentially employed as a programmable nanoreactor to shelter and control the catalytic process.[24−28] Since most of the penetrating molecules in these applications, such as salt, ligands, and reactants, are charged, the understanding and control of the behavior of simple and molecular ions inside hydrogels are critical for the development of high-performance materials for the applications mentioned above.

Two basic quantities that need first attention here are the solubility and the mobility (or diffusivity) of the molecules inside the polymer. The solubility is usually simply expressed by the partition ratio K, which describes the ratio between mean solute concentrations inside the membrane and in bulk.[29] The product of partitioning and diffusivity defines the membrane permeability within the linear response solution–diffusion model.[2] In fact, because of a lower water content, a collapsed hydrogel exhibits an apparently more hydrophobic environment than in the swollen state; therefore, the partitioning of ions is typically below unity (that is, below bulk concentration).[30] In a simple Born solvation picture, one could argue that the mean dielectric constant in a dense hydrogel is an order of magnitude lower than in bulk water, and thus, the transfer from bulk is strongly penalized by a lack of polarization. This penalty (see eq later) would lead to partition ratios orders of magnitude smaller than observed for simple salts such as NaCl,[30] a fact that is essentially unexplained. More strikingly, more complex, molecular ions exhibit a much larger span of partitioning values, in some cases (e.g., small ionized drugs) even exceeding unity by orders of magnitude, K ≃ 103, in collapsed hydrogels.[23,31] This clearly eludes any simple electrostatic and mean-field mechanisms. Hence, despite a rich history of research on synthetic gels and ion exchange polymers, the present knowledge still does not suffice for quantitative, predictive connections between the polymer structure and the thermodynamics of ion solvation in dense membranes.[1,10,22]

The reason for the lack of understanding must be sought in the complex interplay between water, ions, and the polymer matrix on the molecular level, particularly in dense polymers and collapsed hydrogels, where a very crowded and heterogeneous molecular environment can be expected.[32] In fact, it has been shown by us recently using molecular simulations that the water–polymer spatial heterogeneity in a collapsed PNIPAM polymer affects the distribution of neutral solutes;[33] therefore, one may anticipate something similar for ions. We studied previously also the diffusion of ions in the same dense polymer matrix where magnitude and scaling of diffusion with ionic size were very different from nonpolar solutes,[34] which could be clearly traced back to local hydration effects. When comparing the interaction and adsorption of simple monatomic ions and molecular ions[35] to a single PNIPAM chain, we observed that the former were repelled from the chain and stayed nicely hydrated in bulk, while the latter could easily adsorb to the chain. Hence, ionic solvation in a dense, weakly hydrated polymer matrix is not easily accessible and is challenging to categorize without a deeper molecular level understanding.

In this study, we aim for such a detailed understanding of the problem of ion partitioning in the dense, weakly hydrated polymers by using classical molecular dynamics simulations. As a case study we choose the popular PNIPAM polymer in the collapsed state, above the transition temperature, where the water amount is low (20 wt %). We perform a detailed analysis of ionic solvation free energies in bulk water and in the membrane as well as of hydration structure in the polymer matrix. We discuss in detail the electrostatic contributions to the salt and ionic partition ratios and devise a simple solvation theory that explains the simulation outcomes. At the end we devote special attention to the comparison of simple monatomic ions to the molecular ion 4-nitrophenolate (NP–), which behaves markedly different because of weaker hydration and its affinity to the heterogeneous nanosized interfaces within the membrane. This study thus provides molecular insight into ion partitioning in a dense polymer and will be helpful for the interpretation of various experimental data, directly or indirectly related to ionic concentrations inside the polymeric functional material.

Results and Discussion

Background

The uptake of salt ions and larger charged (druglike) molecules by a hydrogel depends crucially on the molecular details (the “chemistry”) of the polymer matrix. It is characterized by the partition ratio K(29) and is for low enough concentrations independent of the concentration.

In Figure we show the correlation between the NaCl salt (K(salt)) and water (Kw) partition ratios (the latter defined as the ratio of water densities inside and outside the gel) for several uncharged polymers reported in the literature. Here, the primary polymer [methacrylate,[36] cellulose acetate,[37] silicone,[30] PNIPAM,[38] or polyethylene glycol (PEG)[39−41]] was cross-linked and either exposed to different temperatures or copolymerized with other monomers in order to tune the equilibrium water content. Quite generally, the uptake of salt by a polymer depends to some extent on water amount; polymers that sorb more water often tend to sorb more salt than those polymers that sorb less water.[42] Also, most of the data fall below the dashed line, depicting an apparent limiting case of K(salt) = Kw where the salt concentration in the sorbed water is equal to that in bulk water.[36] However, since mostly K(salt) < Kw, this indicates that both polymer–ion and polymer–water interactions influence salt partitioning.

Figure 1

Correlation between NaCl partition ratio, K(salt), and the water partition ratio, Kw, for several polymers reported in the literature: methacrylate,[36] cellulose acetate,[37] silicone,[30] PNIPAM,[38] and PEG[39−41] at different temperatures or with different degrees of copolymerization.

On the other hand, larger charged molecules (e.g., ionized drugs) exhibit a much larger span of partitioning values, in some cases even exceeding K = 103 in collapsed hydrogels.[23,31] This clearly eludes the empirical observations for monatomic salts and calls for additional mechanisms, which we will address in this study.

Transfer Free Energies

In this computer simulation study, we use a model of a collapsed PNIPAM polymer phase above the transition temperature with the water partition ratio Kw = 0.2 (see Figure A for a single extended polymer chain). This is also a good model of a collapsed hydrogel with very low (few percent) cross-linker concentrations. A typical simulation snapshot is shown in Figure B with 48 polymer chains in blue and 1325 water molecules in red–white. Removing the polymer component from the plot (Figure C) reveals that water molecules are very nonuniformly distributed, forming irregular, lacy-like clusters, which engender a heterogeneous polar–nonpolar environment. Such water clustering and aggregation are otherwise well-known to occur in various amorphous polymer structures.[43−47] The radius of gyration Rg of the clusters in our system scales with the number of water molecules Nw as Rg ∼ Nw1/2 and hence retains some characteristic of the random walk. More details about the structure can be found elsewhere.[34]

Figure 2

(A) PNIPAM polymer chain (in a stretched conformation) and ions in this study. Snapshots of condensed polymers and water containing three chloride ions, showing (B) all components and (C) water and ions only. The blue bubble zooms into the proximity of one of the Cl– ions, featuring its hydration shell.

For simplicity, we restrict ourselves to the solvation of monovalent ions: three alkali metal cations (sodium, potassium, and cesium), two halides (chloride and iodide), and a molecular ion 4-nitrophenolate (NP–), which is 4-nitrophenol (NP0) with a deprotonated hydroxyl group, see Figure A. Nitrophenol became popular in model reactions in nanocatalytic benchmark experiments and was also used in a connection with PNIPAM hydrogels.[27]

In order to analyze the ionic partitioning we resort to free energy calculations (see the Methods section). Note that a direct evaluation of K in our system is not tractable due to a very slow diffusion of ions (see the Methods section). Figure shows the evaluated solvation free energies of single ions in water (Gw, blue triangle symbols) and in the gel (Gg, orange square symbols). They are negative and comparable in size (of the order of several −100 kJ/mol). Their difference corresponds to the central quantity of interest in this study, namely, to the transfer free energy from water into the gel, ΔG = Gg – Gw. Note that ΔG represents the real single-ion solvation free energy because it is the one encompassing the totality of the reversible work associated with the physical process of transferring the ion from water into the gel. In general, the experimental determination of single-ion solvation free energies is complicated because of the electroneutrality of macroscopic matter.[54] Note that the gel–water interface potential is not included in ΔG because it is screened by ions sufficiently away from the interface compared with the Debye screening length (see the Supporting Information).

Figure 3

Free energies of ions: solvation free energy in water (Gw), solvation free energy in PNIPAM gel (Gg), and the transfer free energy from water into the gel (ΔG = Gg – Gw). For all monatomic ions we use the Jorgensen force field[48−50] and additionally *Åqvist[51] for Na+ and †Dang[52,53] for I–.

The computed values of ΔG are depicted by red circles in Figure . Evidently, we are confronting a very surprising and important observation: there is a stark contrast between the cations and anions in terms of their transfer free energies, ΔG! Systematically, all the monatomic anions span in the narrow window 50–70 kJ/mol and are thus intrinsically repelled from the gel, whereas cations cover the negative range between −30 and −40 kJ/mol and should be attracted by the gel. This difference stretches far beyond ion-specific effects (interpreted as the variation among different ions and force fields of the same valency). Note that NP–, as a molecular ion, includes additional contributions and will be discussed separately.

This is surprising since in the traditional view one expects that charged entities should be universally repelled from a less polarizable medium. Furthermore, also atomistic simulation studies of isolated PNIPAM polymers and swollen networks utilizing similar or even the same force fields[35,55,56] show a universal repulsion of small cations and anions from the chains. Somehow, the aggregation of the polymers and expulsion of water (from excess of water down to 20 wt % in our case) inverts the scenario for cations but not for anions. Because the effect is triggered by reducing the water amount, this also suggests that the cation–anion asymmetry does most probably not stem from polymer–ion interactions but rather has something to do with the water structure and its amount.

Monatomic ions

We will first focus on the mechanisms operative in the solvation of monatomic ions. For that, we begin by taking a closer look into the structural properties of the solvation. It has already been shown that the water–polymer spatial heterogeneity in this collapsed polymer affects the distribution of neutral solutes;[33] therefore, one can expect a related behavior for ions. Indeed, already glancing at the representative snapshot in Figure C reveals that ions (Cl– in this case) enclose themselves with clusters of water molecules. To put this into a quantitative perspective, we analyze the local water densities around the ions in bulk water and in the gel phase, see Figure A. The main peak, corresponding to the first hydration shell, does not weaken that much upon entering from water into the gel. Instead, the ions preserve their first hydration shell even though the mean water density in PNIPAM is only around one-fifth that in the bulk. We conclude that the ions mostly distribute within the aqueous nanoclusters.

Figure 4

(A) Water density profiles around ions in bulk water (dotted lines) and PNIPAM (solid lines). (B) Ion-hydrating cluster size distribution for different ions (symbols) and size distribution of all water clusters in the gel (dashed line). The green-shaded region indicates the mean ion-hydrating cluster size. (C) Continuum picture of an ion encapsulated by a water cluster of an effective radius rw embedded into a hydrogel environment.

Further understanding of the extent of the surrounding water can be obtained by analyzing the hydrating clusters around the ions. We define a water cluster as the group of the water molecules that are mutually separated by less than 0.35 nm. An ion resides in a cluster if it is separated from any of the water molecules of the cluster by less than 0.4 nm. Figure B shows the size distributions of all water clusters (dashed line) in the gel as well as ion-hydrating clusters (symbols), which are the clusters formed around ions. The distribution of water clusters first roughly follows a power-law[34]P(Nw) ∼ Nw–1.74 and crosses over into a roughly exponential decay for larger clusters (Nw > 100). Similarly, the distribution of ion-hydrating clusters also approximately follows an exponential decay for Nw > 10. The mean ion-hydrating cluster size is ⟨Nw⟩ = 120–160, indicated by a green-shaded stripe in Figure B. However, the probability of very small ion-hydrating clusters rapidly decays to zero for Nw < 10, indicating that an ion does not very likely become dehydrated. Another important observation is that the results do not reveal any significant difference between the cations and anions within the accuracy of the data—all the ions behave very similarly.

For the purposes of a simple analysis, we convert the mean ion-hydrating cluster size ⟨Nw⟩ into an equivalent size of a spherical water droplet of the bulk water density (ρw = 32.2 nm–3)which results in the effective droplet size of rw = 0.96–1.06 nm.

A traditional mode of thinking about ionic solvation is based on the continuum picture of the Born free energy,[57] which is particularly powerful for estimating partitioning in homogeneous liquids.[58] In this framework, the water and the gel phases are regarded as homogeneous background media with relative permittivities εw and εg, respectively. The transfer free energy is then associated with the change of the Born free energy of an ion with charge q as[57,59−61]

Here, aB is the effective Born radius of the ion,[62,63] which characterizes its effective size, usually regarded as the location of the first hydration layer. The first hydration layer can be deduced from water density distributions around the ions (aB = 0.23–0.36 nm, cf.Figure A). Using evaluated values for the relative permittivities in our model, εw = 67 for bulk water and εg = 8.5 for the gel (see the Methods section), we obtain ΔGB = 20–30 kJ/mol. However, since in our system the ions in the gel are encapsulated by water clusters (Figure C shows a simplified picture), the effective Born radius should rather reflect the size of the ion along with its hydration shell of effective size rw. In other words, upon the transfer of an ion from water into the gel, the dielectric environment around the ion changes only beyond the radial distance of r = aB. Therefore, we can conveniently assume aB = 1 nm in eq , which gives us ΔGB = +7 kJ/mol for both cations and anions. Note that eq is symmetric on the sign of the ionic charge q and strictly positive as long as εw > εg.

The analysis thus far demonstrated no differences in the structure and distribution of monatomic ions inside the gel. Since ions are enclosed by rather bulky water clusters, this raises the question of the influence of these clusters on the containing charge. An effect that is not accounted for in the implicit Born solvation model is the polarization of water clusters, which comes about due to preferential orientation of water molecules at the cluster interface. The polarization at the cluster interface arises from the fact that the water molecules in this region are subject to an anisotropic environment.

To answer this question, we evaluate the electrostatic potential inside water clusters of various sizes (Figure A). The radial dependence of the potential from the center of a cluster is plotted in Figure B for three different cluster sizes. Interestingly, in spite of an irregular, not well-defined shape, the interior of a cluster seems to universally acquire a well-defined potential of ψcl ≈ −0.5 V with respect to the surrounding “dry” part of the gel. In the same plot we also indicate the potential drop across the planar macroscopic PNIPAM–water boundary (ψs = −0.53 V, as obtained from separate simulations, see the Supporting Information). Already a few water molecules (Nw ∼ 10) are enough to generate the potential drop in the center that is almost equivalent to the drop across macroscopically large interface. This means that an ion enclosed in a cluster is directly subjected to the potential drop at the interface. Cations gain a favorable negative contribution of eψcl = −49 kJ/mol, whereas anions are penalized by this same amount, +49 kJ/mol.

Figure 5

(A) Snapshots of water clusters of various sizes Nw in PNIPAM. (B) Radial dependence of the electrostatic potential from the center of a water cluster. The horizontal dashed line shows the potential drop across a planar PNIPAM–water interface, ψs, where the thickness of the shaded region represents the numerical uncertainty (see the Supporting Information). (C) Snapshots of water droplets of various sizes Nw in vapor. (D) Radial dependence of the electrostatic potential from the center of a water droplet. The horizontal dashed line shows the potential drop across a planar vapor–water interface.

The phenomenon of the interface potential is well-known from the context of the water–vapor and water–oil interface.[54,64−71] For the sake of comparison we also show a similar analysis for water nanodroplets in vapor (Figure C,D). In all the droplets, even the smallest one (composed of 10 molecules), the potential reaches around −0.6 V with respect to vapor, which is the same as the potential drop at a planar water–vapor interface (horizontal dashed line). The analysis of potential inside the clusters and droplets shows a similar behavior and reveals essentially the same physics. Note that the potential profiles tend to smear out for larger clusters because the center of a larger cluster may reside outside the cluster owing to its random-walk shape.

In Figure A we again plot the transfer free energies from the simulations (red circles) along with the contribution from the cluster potential, ±eψcl (blue dashed lines). Clearly, the cluster potential alone already roughly captures the difference of simulation data between the monatomic cations and anions. The remaining deviations can be reconciled by adding the Born solvation component,

Figure 6

(A) Transfer free energies of ions from simulations (red circles); the same as in Figure . The blue dashed lines depict the values ±eψcl, the green solid lines are the predictions of eq for monatomic ions, and the orange lines are the predictions of eq with α = 0.7 and 1 for NP–. (B) Partition ratios K from simulations computed viaeq assuming 1:1 electrolytes for different cation–anion combinations (symbols). The green line is the prediction of eq , and the orange lines are the predictions of eq with α = 0.73 and 1 for combinations involving NP–. The results are based on the Jorgensen force field for ions,[48−50] and where marked additionally *Åqvist[51] for Na+ and †Dang[52,53] for I–.

Accepting the Born part of 7 kJ/mol (assuming aB = 1 nm in eq ), the predictions of eq are plotted by solid green lines in Figure A and capture the data very well.

With the transfer free energy at hand, we are finally in the position to tackle the ion partitioning. The individual partition ratio of ion i is defined as the Boltzmann factor of its transfer free energywhich does not include collective electrostatic effects from other ions and should not be confused with the partition ratio K( of ion i. Because of the electroneutrality requirements in the gel, the partitionings K( of different ion species are coupled and depend on the electrolyte composition (see the Supporting Information). In the simplest case of a 1:1 electrolyte, the concentrations of cations (+) and anions (−) are equal, and their collective partition ratio is given as the geometric mean of both individual partition ratios

In Figure B we plot the salt partitioning obtained from the simulation transfer free energy data (symbols) for all combinations of cation and anion species of a 1:1 salt. For the monatomic combinations, the partition ratios fall inside the window K(salt) = 0.04–0.2, which is also in the range of reported experimental values for numerous polymers (including PNIPAM) with a water content of Kw = 0.2 (Figure ).

An important observation is also that partitioning does not significantly depend on the type of salt (i.e., ion-specific effects are below the numerical resolution), which is in line with simulation studies of core–shell PNIPAM membranes[72] as well as with experimental observations.[38] An easy explanation for this weak ion specificity based on our picture (Figure C) is that their specific character is shielded by the surrounding strongly bound hydration shell; that is, the microenvironment does not specifically change upon insertion into the polymer. However, we can further speculate that the ion specificity may come forth once the hydration cluster starts to vanish, which may happen at even lower hydrations or for larger and more chaotropic ions.

For further discussions, it is useful to define the cluster-potential partition ratiowhich represents a hypothetical partition ratio of a cation solely due to the cluster potential, and conversely Kcl–1 ≃ 3.0 × 10–8 for an anion. However, these enormous values are not disclosed in the total salt partitioning, since they exactly cancel out (cf.eq ). This is because of equal amounts of cations and anions that are enclosed by the clusters, and therefore, the cluster potential remains “concealed”. Namely, the theoretical approximation given by eq together with eq for a monovalent salt leads towhere only the contribution from the Born solvation survives and yields K(salt) ≃ 0.08. This value, indicated by a green solid line in Figure B, matches the simulation values for the monatomic ions considerably well. For a comparison, the standard homogeneous Born solvation approach, where one uses the position of the first hydration shell as the effective Born radius, which then yields the Born part of 20–30 kJ/mol (see above), results in substantially too low values of K(salt) ≃ 10–5–10–3.

The modified Born solvation model (eq ) furthermore predicts that the ionic solvation should decrease if the ion hydration layer (and with that the Born radius aB) decreases, which happens, for instance, when the polymer is further dehydrated upon heating. Indeed, simulations at higher temperatures (thereby at lower equilibrium water amount, Kw) show that NaCl partitioning goes down, which is quantitatively well-captured by the Born solvation model (see the Supporting Information). In the opposite limit of a very swollen gel at low temperatures, the Born model correctly predicts that when aB tends to infinity, and εg tends to εw, then ΔGB vanishes, and K(salt) approaches unity.

Finally, our conclusions for salt partitioning are based on single-ion transfer free energies, valid at low enough concentrations. Certainly, one can expect that, at higher salt concentrations, the ion–ion interactions could add an essential contribution to their solvation. We tested the free energy calculations also at finite salt concentrations and concluded that the single-ion results should be valid up to at least several 100 mM of bulk salt concentrations (see the Supporting Information).

Molecular Ions

Quite generally, large molecular ions (e.g., ionized medicinal molecules, some pharmacological molecules) that possess a considerable electroneutral part can behave very differently from small monatomic ions.[23,31,73] Also, the physics involved in their hydration is expected to be comparatively more complex. Their geometry and charge density are not spherically symmetric. Our representative for a molecular ion is NP–, where the deprotonation of the OH group for pH > 7.15 creates an ionized oxygen center.[74] As we have already seen, the transfer free energy ΔG clearly deviates from the trends of the monatomic ions (see Figure A).

We first take a look at the hydration of the oxygen atoms in NP–: the deprotonated one O– and the two in the nitro NO2 group, see Figure A. For a comparison, the hydration of monatomic ions (similar as in Figure ) is shown by gray symbols. It can be seen that O– is not hydrated to that extent as compared to the monovalent ions. For instance, the probability that the hydrating cluster consists of more than 5 water molecules (by summing up the probabilities P(Nw) for Nw ≥ 6) is only around 0.8, whereas it approaches 1 in cases of the monatomic ions. On the other hand, the NO2 group is only weakly hydrated (with 0.66 probability of the hydration shell smaller than 5 water molecules). The incomplete hydration of NP– can be attributed to a lowered charge density at the O– site, which exerts a more moderate electric field on the surrounding water molecules, resulting in a more labile and shorter-ranged hydration structure. This is because the lone pair from the oxygen delocalizes via conjugation to the benzene ring and the nitro group, which in turn makes the NO2 group slightly charged (see the Supporting Information for the distribution of charges).

Figure 7

(A) Ion-hydrating cluster size distribution P(Nw) around the oxygen atoms in NP– in comparison with monatomic ions (shown by gray symbols in the log–linear plot in the inset). The dashed line shows the size distribution of all water clusters in the gel. (B) Schematic continuum representation of a large molecular ion that partially sticks out of the hydrating cluster (partial dehydration).

In a simple mechanistic picture that transpires from the discussion above, the charged O– and NO2 centers of NP– manage to occasionally escape the water cluster, as schematically depicted in Figure B. With this in mind, we construct a phenomenological description for the free energy of a molecular ion (M) aswhich is composed of the contribution of the neutral, nonionized form of the molecule (first term), the modified cluster contribution (second term), and the Born contribution (last term). The latter is for simplicity assumed to be the same as for monatomic ions. In our case, the first term is the transfer energy of nitrophenol, NP0 (nonionized NP–), which is ΔGneut(M) = −22(1) kJ/mol as obtained previously.[33] As has been shown, the transfer free energy of a neutral molecule scales to a very good extent with the molecular surface area, ΔGneut(M) ∝ Am(M).[33] The second term is the cluster-potential contribution, now multiplied by a hydration parameter α, a phenomenological parameter that accounts for incomplete hydration. If the entire charged part lies inside the cluster, then α = 1 as for monatomic ions, whereas α < 1 represents situations where the charge is partially dehydrated and partially evading the influence of the cluster potential, which decreases the free energy of a negative charge. The fit of eq to the simulation data point of NP– in Figure A gives α = 0.73. For a comparison we show also the prediction for α = 1, which yields a too high value.

Moving on to the partitioning, we show in Figure B the partition ratios of the 1:1 salts with NP– as the anionic component. The nitrophenol salts partition in the gel much more than monatomic salts do. This can now be easily understood by using the approximate expressions for a monatomic cation ΔG(+) and for the molecular ion ΔG(M) (eqs and 8, respectively) in eq , which results in the expressionwhere Kneut(M) = exp(−ΔGneut(M)/kBT) is the partition ratio of the nonionized form of the molecule, and K(salt) = 0.08 is the theoretical prediction for the partitioning of monatomic salts. Notably, because of an incomplete hydration (α < 1) of one of the salt components (NP– in our case), the cluster potential (expressed as Kcl, eq ) influences the partitioning explicitly, in contrast to the case of monatomic ions. The predictions of eq for α = 0.73 and α = 1 are depicted in Figure B by orange lines and demonstrate that the cluster potential has an important contribution to the partitioning.

Moreover, in many practical scenarios, larger molecular ions are not a component of the electrolyte but are instead typically introduced in trace (submillimolar) amounts into a system containing a simple, monatomic monovalent salt (e.g., NaCl). Such are typical cases of charged drugs/pharmaceuticals in a physiological solution or reactants in catalytic experiments. In the case of a 1:1 salt of a much larger concentration c0 than the concentration cM of the molecular ions, the expression for the partitioning of the molecular ions reads (see the Supporting Information)

For our system, we obtain the partitioning of NP– ions in excess of NaCl in the hydrogel K(M) = 4 × 103, which is an enormous value compared to the partitioning of NaCl salt! A full calculation for the cases with two anions of arbitrary concentrations shows that the drastic consequence of the large partitioning of species M is that only relatively small amounts of it (∼mM) are needed to almost completely exchange the smaller anion; i.e., the partitioning of the latter decreases by orders of magnitude (see Figure S5 in the Supporting Information).

Furthermore, the nonionized form of the molecule, NP0, partitions with a very similar ratio, Kneut(M) = 3 × 103.[33] This may seem counterintuitive and against common knowledge, which would anticipate that charged molecules should be significantly less attracted to weakly hydrated and hydrophobic gels than their neutral counterparts. In order to understand this nonapparent and surprising outcome, we use our phenomenological theoretical framework, which helps us to explain the main qualitative trends. Inserting the theoretical expressions (eqs and 8) into eq gives us

In other words, the partitioning of a charged molecule in excess of salt is proportional to the partitioning of its nonionized form, Kneut(M), which is then modified by the factor Kcl1−αK(salt) due to ionization. In case the charged molecular centers are fully hydrated (α = 1), the ionization decreases the partitioning only by the factor of K(salt) ≃ 0.08. On the other hand, if the charged centers become partially dehydrated (α < 1), this blows up the partitioning in an exponential manner—in our case with α = 0.73, this adds a factor of Kcl1−α ≃ 102. One needs to be aware that eq is very rough and cannot yield reliable quantitative results, but it nevertheless gains important qualitative insights into the partitioning governed by water clusters due to incomplete hydration.

Our simulations as well as the theoretical analysis show that large ionized molecules can exhibit enormous partitioning (K ∼ 103) in weakly hydrated gels, which was reported in numerous experiments of charged pharmaceuticals in PNIPAM hydrogels.[23,31,73] Notably, the charge does not necessarily impede the partitioning but can also enhance it.

Conclusion

The understanding of the solvation and structure of charged molecules within dense polymer systems remains rudimentary owing to a high complexity of the molecular mechanisms involved. Resorting to a classical atomistic simulation framework, we investigate partitioning of small monatomic and larger molecular ions in dense, weakly hydrated neutral polymers. The simulations reveal that ions are enclosed by hydrating, nanosized water clusters that accommodate them within a principally hydrophobic surrounding.

The estimates for partitioning of monatomic salts between 0.04 and 0.2 from simulations agree well with experiments for a wide range of polymers of a similar water content. This partitioning is much larger than inferred from a simple Born solvation model with an effective homogeneous dielectric constant. The hydration of monatomic cations and anions is structurally almost indistinguishable. However, the polarization stemming from the water–cluster interfaces induces a potential drop of around −0.5 V and renders a substantial difference in the transfer free energies of ions from water into the gel: negative for cations and positive for anions. Nevertheless, the explicit impact of the cluster potential on the thermodynamic partitioning disappears for monatomic salts because of overall electroneutral sets of ions that are subjected to this potential. This fact makes the modified Born solvation model that accounts for a 1 nm thick hydration shell sufficient.

The story becomes more intricate for larger molecular ions, in our study represented by 4-nitrophenolate. Because of its labile hydration shell, the charged parts of the molecule can occasionally “step out” of the cluster, thereby escaping the influence of its electrostatic potential. This has far-reaching consequences on the thermodynamics. The cluster potential is no longer completely compensated, and it thus reveals its direct influence on the total partitioning of the ions, which is by itself an interesting phenomenon. Namely, the cluster potential is a special type of an interface potential (the water–air interface potential being the most prominent example), and as such it may be expected to be thermodynamically inaccessible and nullified in the overall salt partitioning due to the net charge in action. However, the fragmented water clusters together with the asymmetry in hydration strength of ion species involved elude this paradigm. A direct consequence of the incomplete cancellation of the cluster potential is a higher ionic partitioning. In fact, ionizing the molecule may even enhance its solubility in the gel, which opposes the standard picture of ion solvation. Indeed, high partitioning has been observed in experiments with charged organic molecules.[23,31] Also, the high partitioning of 4-nitrophenolate explains its fast reduction kinetics in collapsed PNIPAM-based nanoreactors.[27]

Our study provides important and very fundamental insights into the microscopic mechanisms behind the solvation of ions and charged molecules not only in dense hydrogels but also in other amorphous matter with water microphases as well. The notion of the surface potential is extremely vital for various research directions, ranging from desalination membranes to biomedical applications and pharmacokinetics.

Methods

Model and Force Fields

The atomistic computer model of the collapsed PNIPAM polymers with sorbed water is adopted from our previous works.[33,34] In this model, 48 atactic PNIPAM chains built out of 20 monomers are condensed together with 1325 water molecules (corresponding to the water amount in equilibrium with bulk water) in a cubic simulation box (of size ∼6 nm) at 340 K and isotropic pressure of 1 bar.

For PNIPAM polymers, we use the recently improved OPLS-based model by Palivec et al.(75) with an ad hoc parametrization of partial charges that reproduces the thermoresponsive properties better than the standard OPLS-AA force field. For water we use the SPC/E model.[76] We use the force field by Jorgensen and co-workers[48−50] for all monatomic ions, and additionally by Åqvist[51] for Na+ and Dang[52,53] for I–. For the nitrophenolate ion NP–, we use the OPLS-based force field with the partial charge parametrization “OPLS/QM1” from ref (35).

Simulations

All atomistic MD simulations are performed using the GROMACS simulation package.[77,78] Electrostatics is treated using the particle-mesh-Ewald method[79,80] with a 1 nm real-space cutoff. The Lennard-Jones (LJ) interactions are truncated at 1 nm. The simulations are performed with an integration time step of 2 fs in the constant-pressure (NPT) ensemble with periodic boundary conditions. Temperature was maintained at 340 K, using the velocity-rescale thermostat[81] with the time constant of 0.1 ps. The isotropic pressure was maintained at 1 bar using the Berendsen barostat[82] with the time constant of 1 ps.

Free Energy Calculations

The solvation free energies of ions are calculated using thermodynamic integration (TI).[83] To reduce mutual ion–ion effects (note that the Bjerrum length is 6 nm), we restricted the number of ions in the simulation box to three. The net charge is compensated by applying a uniform neutralizing background charge.

We first insert three ions of the same type at random positions into an equilibrated polymer system and equilibrate it further with fully interacting ions for at least 100 ns. The necessary equilibration time is estimated based on the crossover time of the ion to reach the normal diffusion.[34]

During the TI simulations, the Coulombic and LJ interactions of the equilibrated ions are gradually switched off. Introducing a coupling parameter λ ∈ [0, 1] that continuously switches the interactions in the Hamiltonian U(λ) between the original interactions (for λ = 1) and a noninteracting ion (for λ = 0), the solvation free energy is computed as[83]

The integration is performed in two stages: We first linearly scale down the charge, while keeping the LJ interactions intact. In the second stage, we scale down the LJ interactions using the “soft-core” LJ potentials as implemented in GROMACS in order to avoid singularity problems when the interactions are about to vanish (λ → 0).[84] The entire TI procedure is composed of 24 individual simulation windows with equidistant λ values for the Coulomb part and likewise 24 windows for the LJ part. The simulation time of each individual λ window is 4 ns where the first 3 ns is discarded from the sampling to allow an equilibration of the ion. All the TI calculations are performed with 5–6 independently equilibrated systems for the Coulomb part and 1–2 systems for the LJ part (note that the LJ part converges considerably faster than the Coulomb part, allowing also a much more accurate evaluation). The final results are averaged over all three ions in the simulation box and over all the systems. Even though each ion samples only a local phase space during a short TI time window, the used procedure assures adequately sampled values of the free energy.

Because of the Ewald summation with periodic boundary conditions, the free energy of solvation depends on the size of the simulation box, and we have to apply the finite-size correction[85]where q is the ion charge, ε the relative permittivity of medium k [water (w) or gel (g)], and L the size of the simulation box. The correction in our case accounts for around 0.5 kJ/mol in water and 3.8 kJ/mol in the gel, which is comparable to the uncertainties of the evaluated free energies. Therefore, higher-order corrections to eq (e.g., orientational polarization of the solvent due to periodic images of the ion[68]) are not necessary. Finally, the solvation free energies for transferring an ion from vacuum into the medium are obtained as a difference of the values obtained from TI (corrected for eq )Note that GvacTI may be nonzero for molecules (NP–) due to intramolecular Coulombic and LJ contributions.

The evaluated free energies represent the intrinsic single-ion solvation free energies, arising exclusively from the interaction between the ion and its local solution environment, and take neither the surface potential ψs[68] nor the compression free energy[86] into account. We also verified that performing TI on three ions of the same kind with a uniform background charge yields the same results as performing TI on ion pairs (see the Supporting Information).

Need for Free Energy Calculations

In principle, an alternative simulation setup containing a polymer membrane in water provides a possibility for a direct evaluation of K from the equilibrated ion concentrations inside and outside the membrane. The minimal width of both the water and polymer slabs would need to be at least d = 3 nm in order to overcome the interface regions of 1 nm (Debye length at 100 mM). An estimated time for an ion to diffuse from one end of the membrane slab to the other is d2/Dg, where Dg ∼ 10–4 nm2/ns[34] is the diffusion coefficient in the gel membrane. Furthermore, the statistics is hindered by the low ion partitioning K, such that the necessary simulation time for a meaningful sampling would be tsim ≫ d2/(KDg) ∼ 106 ns. This is 3 orders of magnitude longer than for the free-energy calculations and hence far beyond the current rational capabilities. The direct method lends itself, for example, for higher polymer hydrations[72] and implicit-solvent models.[87]

Relative Permittivity

The static relative permittivities of bulk water and PNIPAM gel are calculated from the fluctuations of the total dipole moment M of the system as[88]where V is the volume of the simulation box. The fluctuations are evaluated from the systems without ions.

Electrostatic Potential

In the case of water droplets and clusters, we compute the electrostatic potential along the radial distance from the center of mass (COM) of each target droplet or cluster aswhere the summation runs over all the partial charges q of the system, located at positions r. For each radial distance r = |r| from the COM, we average the results over the entire solid angle around the droplet or the cluster. In the case of clusters, we also average the results over all the clusters of a target size Nw in each stored simulation frame.

The potentials across planar water–vapor and PNIPAM–water interfaces are calculated from independent simulations that consist of two separated phases by double integration of the Poisson equation. For a detailed analysis, see the Supporting Information.