How the Shape and Chemistry of Molecular Penetrants Control Responsive Hydrogel Permeability

The permeability of hydrogels for the selective transport of molecular penetrants (drugs, toxins, reactants, etc.) is a central property in the design of soft functional materials, for instance in biomedical, pharmaceutical, and nanocatalysis applications. However, the permeation of dense and hydrated polymer membranes is a complex multifaceted molecular-level phenomenon, and our understanding of the underlying physicochemical principles is still very limited. Here, we uncover the molecular principles of permeability and selectivity in hydrogel permeation. We combine the solution–diffusion model for permeability with comprehensive atomistic simulations of molecules of various shapes and polarities in a responsive hydrogel in different hydration states. We find in particular that dense collapsed states are extremely selective, owing to a delicate balance between the partitioning and diffusivity of the penetrants. These properties are sensitively tuned by the penetrant size, shape, and chemistry, leading to vast cancellation effects, which nontrivially contribute to the permeability. The gained insights enable us to formulate semiempirical rules to quantify and extrapolate the permeability categorized by classes of molecules. They can be used as approximate guiding (“rule-of-thumb”) principles to optimize penetrant or membrane physicochemical properties for a desired permeability and membrane functionality.


SIZE SCALING OF DIFFUSIVITY
We fit the diffusivity data in the collapsed state (D c ) with three different exponential forms (for n = 1, 2, 3): where D 0 and λ are fitting parameters. The fits for the three different penetrant categories (spherical, planar + Me, linear + Me) are plotted in Figure S1A (for n = 1; the same as Figure 2A in the main text), Figure S1B (for n = 2), and Figure S1C (for n = 3).
We compute the sum of squared residuals, which serves us to compare the quality of the fits:   Figure S1: Diffusion coefficients of penetrants in the collapsed PNIPAM polymer (same data as in Figure 2A in the main text) plotted versus different powers of the penetrants' Stokes radii in water: (a) aw, (b) aw 2 , and (c) aw 3 . The dashed lines are fits of eq S1 to the three different categories of penetrants.
Visual inspection of the fits in Figure S1 as well as noticeably smaller values of S (1) than S (2) and S (3) in Table S1 clearly show that the linear form (i.e., n = 1, which is eq 8 in the main text) provides by far the best description of the data among the three different popular possibilities suggested in the literature.

GLASS TRANSITION TEMPERATURE
We use two approaches to estimate the glass transition temperature, T g , of the collapsed PNIPAM model at 340 K (containing 20 wt% of water): (i) by cooling simulations and (ii) by monitoring water diffusivity at different temperatures.

Cooling simulations:
In the first approach, we perform cooling simulations from an equilibrated structure at 400 K down to 100 K and monitor the density. The glass transition temperature is considered as the temperature at which the density vs. temperature slope undergoes a gradient change. The crossover is sensitive to the analysis-it depends on the cooling technique, statistical variation, cooling rate (Γ), and the fitting range. 1 To probe the latter two dependencies, we perform cooling simulations on three independent samples and at three different cooling rates, Γ = 5, 20, and 100 K ns −1 . For fitting each branch, we use two different temperature intervals, ∆T = 75 K and 100 K (i.e., 100-175 K and 100-200 K for the glassy branch and 325-400 K and 300-400 K for the rubbery branch).
An example of such a density-temperature plot is shown in Figure S2A, where Γ = 5 K ns −1 . The resulting T g averaged over the three independent samples are shown Figure S2B.
Water diffusivity: In the other approach, we perform longer simulations (1000 ns each) at different temperatures (ranging from 160 to 400 K) of the system (20 wt% water) and evaluate the diffusion coefficient of water. For even lower temperatures, the diffusion does not reach the normal diffusion regime within the simulation time, therefore the evaluation is not possible. Figure S3 shows an Arrhenius plot of the resulting diffusion coefficients, where a change in the slope is nicely visible. The crossover temperature obtained from linear fits is T * = 275 K. This agrees well with the T g values obtained from the cooling method in Figure S2B.
From both analyses, we can thus confidently conclude that the glass transition temperature of our system is around T g = 275 ± 10 K, and the temperature to the glass transition temperature ratio is

CHAIN RELAXATION TIME
For analyzing chain relaxation in the collapsed state at 340 K, we use a self-intermediate scattering which is based on the atomic displacements during the time interval t. The angle brackets ... i,t denote the average over all chiral carbon atoms in the PNIPAM backbone and over time t . Our simulated systems are isotropic, therefore, the scattering function only depends on the absolute value of the scattering wave vector q = |q|, determining the length scale at which dynamics is monitored. Evaluated F s (q, t) is plotted in Figure S4A as a function of time for different values of q. The time dependence can be described by a two-step relaxation 3 The first term reflects the fast dynamics modes of the polymer, and the second term describes the main relaxation mode of the polymer matrix. Here, β stands for 'stretching exponent', empirically found β < 1 for most cases. For simplicity, we will assume β = 1 in our fits, which yields smaller values for the relaxation time τ than if assuming β < 1. Since the fast relaxation dies out on the order much smaller than 100 ns, we fit the second term of eq S5 to the data in Figure S4A in the range 100-500 ns. The fitted relaxation times τ as a function of wave vector are shown in Figure S4B. A relevant length scale for diffusion processes is the distance between neighboring polymer chains. From the effective PNIPAM radius R 0 = 0.5 nm and the fact that PNIPAM is a densely packed major component of the system (80 wt%), the distance between chains is around 2R 0 = 1 nm, and the corresponding wave vector is q = 2π/2R 0 ≈ 6 nm −1 . The chain relaxation time at this length scale is around 2 × 10 3 ns (see Figure S4B).

COMPARISON WITH THE SELF-CONSISTENT COOPERATIVE HOPPING THEORY
We compare the scaling of diffusion resulting from our simulations with predictions of the self-consistent cooperative hopping theory (SCCHT) by Zhang and Schweizer. 4 The SCCHT is namely one of a few theories that predict a near exponential decrease of penetrant diffusivity with its size. The theory is based on the activated diffusion of spherical penetrants in dense melts and glasses of matrix spheres. The predictions for diffusion can be approximately described with the equation where σ is the diameter of the matrix sphere, d the diameter of the penetrant, and the dimensionless parameter ξ depends on the matrix volume fraction φ m . By fitting eq S6 to the predictions for various φ m in Figure 3 in ref 4 on the interval of d/σ from 0.3 to 0.6, we obtain coefficients ξ shown in Table S2. As it turns out, the penetrant-matrix interaction does not significantly change the slope (i.e., ξ) but instead modifies only the prefactor D 0 in eq S6 (see Figures 8 and 9 in ref 4).
To quantitatively compare the predictions of the SCCHT with our simulations, we have to map our molecular model to an effective mixture of spheres. To do that, we interpret the matrix sphere radius σ in the SCCHT as the Kuhn length of the polymer chains in our model, σ = L K , as proposed in ref 4. With that, we can rewrite eq S6 directly to eq 8 in the main text by matching the exponents, ξ(d/σ) = (d/2)/λ, which provides the parameter λ in our notation as This allows us to evaluate the effective λ from the SCCHT and compare it to our simulations. Before doing so, we have to extract the Kuhn length of PNIPAM chains from our simulations, which is defined as where r 2 ee is the mean square distance between the ends of the chain (in our case defined by the chiral carbon atoms in the PNIPAM backbone) and L c = 5.0 nm is the contour length of the backbone of the chain. Averaging r 2 ee over time and over all chains in the system (denoted by angle brackets), gives us L K = 1.14 nm.
The calculated values for λ from the SCCHT, obtained from eq S7, is given in Table S2. In the end, we want to find a reasonable estimate for the matrix volume fraction φ m of the hard-sphere system that corresponds best to our molecular model. The volume fraction of the close-packed spheres equals φ max ≈ 0.74, which can be considered as a completely dehydrated, densely packed polymer system. Our system at 340 K contains 20 wt% of water and 80 wt% of the polymer, thus a reasonable mapping to the hard-sphere model is φ m = φ max × 0.80 ≈ 0.59. As seen from Table S2, at this particular volume fraction the value for λ from the SCCHT is 0.020 nm, whereas our simulations give λ = 0.0186 nm for spherical penetrants ( Table 1 in the main text), which is very good agreement.

COMPRESSIBILITY OF THE COLLAPSED PNIPAM STATE
In order to evaluate the compressibility of the collapsed PNIPAM state, we perform additional simulations (three independent samples) of the PNIPAM at 340 K at the pressure of 1000 bar for 70 ns. We discard the first 10 ns for relaxation reasons and evaluate the volume, shown in Table S3. The isothermal compressibility is then computed as which gives χ p = 3.3(1) × 10 −5 bar −1 . This is smaller than the compressibility of bulk water at 340 K, χ w = 4.65(1) × 10 −5 bar −1 .

SIZE SCALING OF TRANSFER FREE ENERGY
We fit the transfer free energy data from water into the collapsed PNIPAM state (∆G) with two different forms which are shown in Figure S5.  Figure S5: Transfer free energies of penetrants from water into the collapsed PNIPAM (same data as in Figure 2B in the main text) plotted versus penetrant (A) surface area and (B) volume. The dashed lines are fits of eqs S10 and S11 to the different categories of penetrants in panels A and B, respectively.
We asses the quality of the fits by computing the sum of squared residuals, (S13) Table S4: Sum of squared residuals (eqs S12 and S13, in the units of kJ 2 mol −2 ) of the fits in Figure S5. Comparing the fits in Figure S5 by visual inspection, it is difficult to discern the difference in the fitting quality. However, the squared residuals in Table S4 demonstrate better fits with surface area scaling (eq S10) for all three groups.