Tuning the selective permeability of polydisperse polymer networks

Abstract

We study the permeability and selectivity (‘permselectivity’) of model membranes made of polydisperse polymer networks for molecular penetrant transport, using coarse-grained, implicit-solvent computer simulations. In our work, permeability [scr P, script letter P] is determined on the linear-response level using the solution–diffusion model, [scr P, script letter P] = [scr K, script letter K]D_in, i.e., by calculating the equilibrium penetrant partition ratio [scr K, script letter K] and penetrant diffusivity D_in inside the membrane. We vary two key parameters, namely the network–network interaction, which controls the degree of swelling and collapse of the network, and the network–penetrant interaction, which tunes the selective penetrant uptake and microscopic energy landscape for diffusive transport. We find that the partitioning [scr K, script letter K] covers four orders of magnitude and is a non-monotonic function of the parameters, well interpreted by a second-order virial expansion of the free energy of transferring one penetrant from a reservoir into the membrane. Moreover, we find that the penetrant diffusivity D_in in the polydisperse networks, in contrast to highly ordered membrane structures, exhibits relatively simple exponential decays. We propose a semi-empirical scaling law for the penetrant diffusion that describes the simulation data for a wide range of densities and interaction parameters. The resulting permeability [scr P, script letter P] turns out to follow the qualitative behavior (including maximization and minimization) of partitioning. However, partitioning and diffusion are typically anti-correlated, yielding large quantitative cancellations, controlled and fine-tuned by the network density and interactions, as rationalized by our scaling laws. We finally demonstrate that even small changes of network–penetrant interactions, e.g., by half a k_BT, modify the permselectivity by almost one order of magnitude.

---

1 Introduction

Being a key transport property in materials science, the permeability of membranes has been excessively studied for more than a century.^1–5 The permeability determines the fundamental ability of functional solutes, such as ions, ligands, proteins, and reactants to penetrate and be transported through dense but permeable membranes of various kinds. Membranes, constituting typically quite crowded environments, are mostly polymer-based, and are ubiquitous in soft matter applications, materials science, and naturally in biological systems. In the latter, bio-hydrogels such as cytoskeletons, mucus gels, and the extracellular matrix (ECM) are complex molecular assemblies composed of hydrated polymer networks.^4,6–11 In general, they function as selectively permeable barriers for solutes to penetrate.⁹ For instance, the ECM constructs a selective barrier around the cells, thereby regulating the transport of signaling molecules.^9,12–18 Hence, the permeability of bio-hydrogels plays a decisive role in maintaining life.

Other important examples of polymer-network-based membranes can be found in functional soft matter composed of synthetic hydrogels, such as cross-linked poly(N-isopropylacrylamide) (PNIPAM).¹⁹ Due to their thermoresponsiveness and relatively sharp volume transition, they are widely used as representative and promising components in emerging material technologies for stimuli-responsive carrier particles, actuators, sensors, or responsive nanoreactors.^20–33 In the latter, for instance, the hydrogel embeds nano-sized enzymes or metal nanoparticles catalyzing chemical reactions, which are ultimately controlled by responsive membrane permeability.^34–36 In general, responsive polymeric matrices can be expected to control the permeation of solutes (penetrants) in a selective manner, modulated by external stimuli such as temperature, pH, and salinity. The tunable selectivity of permeability (‘permselectivity’)⁴ thus bears enormous potential for the development of ‘intelligent’, programmable and adaptive membranes for diverse applications ranging from gas separation,^37–43 water purification, and filtration^44–50 to dialysis and drug delivery.^51,52

Typically, the permeability of dense membranes is quantified by the so-called solution–diffusion model on the linear-response level, via^{3,5,37,38,43,53–60}

[scr P, script letter P] = [scr K, script letter K]D_in, (1)

which is the product of two key quantities measured in equilibrium, namely the partitioning (partition ratio) [scr K, script letter K] = c_in/c₀, simply given by the ratio between the penetrant concentrations inside the membrane c_in and in the bulk reservoir c₀, and the diffusivity (diffusion coefficient) D_in of the penetrants inside the membrane. In this definition, the permeability of the bulk reservoir reference is thus equal to the free penetrant diffusivity in the bulk reservoir, i.e., [scr P, script letter P]₀ = D₀. Definition (1) is based on linear Fick's type of laws,^5,56,61 where penetrant flux j is generated by a concentration difference Δc, i.e., j = −[scr P, script letter P]Δc/L across the membrane of thickness L, and the proportionality constant is given by eqn (1). Note that in the gas separation field, the permeability is often very similarly defined as the product of a solubility coefficient and a diffusion coefficient,³⁷ which may lead to different physical units of [scr P, script letter P], depending on the exact definition of solubility.⁶² Other definitions of permeability exist depending on length scales and transport mechanisms, e.g., Knudsen diffusion or convection, which are not applicable to our simulations.⁶³ The elegance of the model represented by eqn (1) is that it is simply based on two intuitive and fundamental properties of a medium in equilibrium, which should be easily accessible in experiments and theoretically tractable.

However, there is still a growing number of theoretical studies pursuing a better understanding of partitioning^{39–41,64–71} and diffusivity^{37,53,55,56,72–94,131} in polymer-based membranes and hydrogels. It is the complexity arising from diverse molecular interactions (e.g., excluded volume and attraction) and conformational structures (cross-linked, ordered, polydisperse) inside the membrane that renders the problem very challenging. Since the pioneering attempts^95–99 to tackle this problem, computer simulations have become powerful, requisite tools used to understand the permeability in polymers.^100–110 In this context, for instance, we recently presented a simple coarse-grained (CG) simulation model of penetrant transport across a rigid immobile lattice-based membrane, pursuing a better comprehension of permeability, particularly in dense and attractive systems.⁸⁷ Despite the simplicity of that model, we demonstrated a very intricate behavior of permeability: it varied over many orders of magnitude, and could even be minimized or maximized by tailoring the potential energy landscape for the diffusing penetrants through small variations of membrane attraction, structure, and density. Supported by limited scaling theories, we showed that the possible occurrence of extreme values is far from trivial, being evoked by a strong anti-correlation and substantial (orders of magnitude) cancellation between penetrant partitioning and diffusivity, especially in the case of dense and highly attractive membranes.

In this work, we extend the previous study of a fixed, ordered membrane topology to a more complex and more realistic case of a membrane composed of fluctuating and cross-linked polydisperse polymers to study the transport of diffusive penetrants. For this, we consider a polydisperse tetra-functional network, i.e., each cross-linker connects four polymer strands, which have a polydisperse length distribution. As considered previously,^69,87 the system includes the network region and the bulk reservoir region, enabling a direct calculation of partitioning, diffusivity, and thus permeability. We focus on two important control parameters: the polymer network density ϕ_n (volume fraction), tuned by internal interactions, and the interaction between the network monomers and the penetrants. We calculate the linear-response permeability [scr P, script letter P] according to eqn (1) and systematically analyze and rationalize our findings by presenting semi-empirical scaling laws. Finally, we demonstrate how minute changes of the interactions can modify the permselectivity of the membrane substantially.

2 Methods

2.1 Simulation model

2.1.1 Network structure and setup. We performed implicit-solvent Langevin dynamics computer simulations of a model membrane made of a polydisperse, tetra-functional polymer network^111–114 including diffusive penetrants (see Fig. 1a). As described in our previous studies on penetrant partitioning in regular polymer networks,⁶⁹ and permeability of highly ordered membranes,⁸⁷ we used LAMMPS software¹¹⁵ with the stochastic Langevin integrator. The iteration time step δτ = 5 × 10⁻³τ was used with the time units , where m is the unit mass, and k_BT = 1/β is the thermal energy. The friction coefficient γ was chosen such that the momentum relaxation time τ_γ = m/γ = τ, thus the free penetrants’ motion becomes diffusive after 200 time steps.

Fig. 1 (a) Simulation snapshot of the polydisperse tetra-functional polymer network in the swollen state with diffusive penetrants (blue). The polymer segments (red beads) are connected by tetra-functional cross-linkers (yellow beads) and have a random length distribution (see the text in the Methods section). (b) Various network conformations depending on ε_nn and ε_np. The network collapses as the the network–network interaction parameter ε_nn increases (i.e., lowering the solvent quality to poor solvent conditions). (c) Polymer network volume fractions ϕ_nvs. the solvent quality parameter ε_nn at different values of the network–penetrant parameters ε_np (see the text for details). The solid line interpolates between the mean over the volume fraction averaged over all simulated ε_np values.

For the initial configuration of the network, we considered 4 × 4 × 4 unit cells of a diamond cubic lattice, where N_xlink = 64 × 8 = 512 cross-linkers were located on the lattice points. The unit cell length is 25σ (with σ defining the penetrant size and our length scale) with the initial normalized cross-linker positions per cell, {(0,0,0),(0,0.5,0.5),(0.5,0,0.5),(0.5,0.5,0),(0.25,0.25,0.25), (0.25,0.75,0.75),(0.75,0.25,0.75),(0.75,0.75,0.25)}. In the tetra-functional network, each cross-linker connects four polymer chains. The number of polymer monomers between the (closest neighboring) cross-linkers was randomly drawn from a uniform distribution between 2 and 18, thereby resulting in the polydisperse structure with an average chain length of 10 monomers, and a standard deviation of about 5. With the above construction we ended up with N_m = 10 364 monomers in the network, yielding a cross-linker fraction of f_xlink = 4.7%. This cross-linker fraction is in the range of typical experimental values for tetra-functional polymer networks, such as cross-linked PNIPAM hydrogels.^88–92

For initial equilibration, the membrane was placed in the middle of a simulation box of lateral lengths L_x = L_y = 100σ and longitudinal length L_z = 300σ, with periodic boundary conditions in all three Cartesian directions. For a small molecule of typical size σ = 0.4 nm, the longitudinal box length then corresponds to L_z = 120 nm. Note that the periodic boundary conditions lead to an infinitely large membrane in the x and y directions, while the membrane-reservoir regions periodically repeat in the z direction with a spatial period L_z, reminiscent of a stack of lamellar hydrogel layers.¹¹⁶ The membrane was first equilibrated in the NVT ensemble using the force-field described below. We then added N_p = 1000 penetrant particles to the bulk reservoir region, and equilibrated the whole system.

In the next step, the longitudinal box length L_z was kept fixed, while L_x and L_y could adjust equivalently according to the NpT ensemble with a given particle number N = N_m + N_xlink + N_p = 11 876, pressure p_x = p_y = p, and temperature T. To maintain the target pressure p on an average, we used the Berendsen barostat with semi-isotropic pressure coupling.¹¹⁷ Since we consider the implicit-solvent model, the anisotropic pressure couples only to the virial coefficients in the x and y directions of the polymer monomers and the penetrants. The value of the lateral pressure defines the network reference states without penetrants and was chosen to be p = 1 bar ≈10⁻²k_BT nm⁻³. The latter is in fact negligibly small (essentially p ≈ 0) compared to the interactions in our system,^67,118 which are of the order of k_BT nm⁻³ ≈ 41 bar ≫ p. The pressure relaxation time τ_p and bulk modulus K_b for the barostat were chosen in the range of 1 ≤ τ_p/τ ≤ 2 and 1 ≤ K_b/p ≤ 10, respectively, to maintain reasonable dynamic pressure relaxation comparable to the Langevin momentum relaxation, which depends on the interaction parameters.

After an equilibration time of 1.5 × 10⁵τ, we performed the production simulations typically up to 10⁷δτ = 5 × 10⁴τ. As the finite network membrane is connected to the large bulk reservoir region of solute penetrants, the penetrants can always equilibrate their partitioning between the large bulk reservoir and the (responsive) membrane. The resulting equilibrium density profiles of the monomers as well as the penetrants within the membrane and in the bulk reservoir are considerably homogeneous along the z direction, cf. Fig. S1 in the ESI† demonstrates that the systems are well equilibrated. For the chosen typical molecule size σ = 0.4 nm, the corresponding membrane widths, determined from the monomer concentration profiles, equilibrate then typically from ≃4 to 16 nm. The lateral equilibrated lengths of the simulation box L_x (= L_y) are shown in Table S2 in the ESI.†

Selected two-dimensional radial density distribution functions between the cross-linkers, g^2D_xlink(r), shown in Fig. 2, also demonstrate that the equilibration procedure leads to reasonable and homogeneous network structures in the lateral directions (note that g^2D_xlink(r) is averaged over thin two-dimensional membrane slabs in the x and y directions, see Fig. S2 in the ESI† for details). Especially in the dense state, apart from some short-ranged packing effects, for r ≳ 3σ, g^2D_xlink(r) becomes homogeneous in the lateral directions. For the swollen network, g^2D_xlink(r) reveals some more structures with a local peak at 3σ ≲ r ≲ 4σ, reflecting short-range correlations between the crowded cross-linker regions, and a second peak close to the average chain length (i.e., average mesh size in the swollen case) of r ≃ 10σ. Changing the network–penetrant interaction affects these distributions only slightly in the dense systems, while in the swollen case few structures are observed for large attractions between the network and the penetrants (see Fig. S2 in the ESI†), indicating some amount of penetrant-induced homogenization.

Fig. 2 Two-dimensional radial distribution function between the cross-linkers g^2D_xlink(r) for three different polymer volume fractions, from swollen (ϕ_n = 0.06) to collapsed (ϕ_n = 0.36) states. The network–penetrant interaction parameter is βε_np = 0.1. The distribution function g^2D_xlink(r) is averaged over thin 2D membrane slabs in the x and y directions, see the ESI† for details and more data.

2.1.2 Force-field. For the non-bonded interactions, all particles (i.e., monomers, crosslinkers, and penetrants) interact via the generic Lennard-Jones (LJ) potential U^ij_LJ for i, j = n or p, where n denotes the network particles (polymer monomers and cross-linkers), and p denotes the penetrant. The LJ size is σ for all particles and sets the length unit. (Note that our model is quite generic and, while we have molecular sizes σ ≲ 1 nm in mind, it can be applied also to, e.g., macromolecular penetrants of larger size.) The LJ energy ε_pp = 0.1k_BT of the LJ potential U^pp_LJ is fixed such that the mutual penetrant interaction is essentially repulsive^69,87 (see also the positive second virial coefficient of the LJ interaction in Table S1 in the ESI,† and the following sections including Fig. 4 for details of the virial coefficients). In this work, we vary two interaction parameters, the network–network interaction ε_nn, and the network–penetrant interaction ε_np between 0.1 and 1.5k_BT. Here we should mention that, in fact, we considered even higher attractions, βε_np > 1.5, but the simulations revealed metastable non-equilibrium phases where penetrants were kinetically trapped within the membrane, evading equilibration. For details, see Section D and Fig. S3 in the ESI.†

The intra-network interaction ε_nn is interpreted as a measure of solvent quality,^69,87,119 thereby controlling the network volume fraction ϕ_n = (N_m + N_xlink)v₀/V_n, where v₀ = πσ³/6 is the monomer volume with diameter σ = σ_nn = σ_np, and V_n is the network volume. As discussed in previous works,^69,87,119 small/high ε_nn corresponds to a good/poor solvent leading to a small/high volume fraction, respectively. The network–penetrant interaction ε_np governs the strength of attraction between the polymers and the penetrants.

For the bonded interactions of the (semi-flexible) polymers, we employed harmonic stretching (bond) and bending (angle) potentials.⁶⁹ The bonded polymer parameters were determined via coarse-graining from explicit-water, all-atom simulation results of cross-linked PNIPAM chains, utilizing a force-field from our group's work.¹²⁰ Since the cross-linker connects monomers of four polymer chains, the network is tetra-functional, and in addition to the m-m-m bending, there are six bending potentials for the m-xlink-m arrangement. Therefore, we have nine different bonded (7 bending (angles), 2 stretching (bonds)) potentials in total and we determined eighteen bonded parameters K^ij_r, r^ij₀, K^ijk_θ, and θ^ijk₀ by fitting harmonic potentials to the free energies obtained from the all-atom simulations. The details of all the bonded interactions, that is, their calculation from the all-atom (explicit-water) simulations of PNIPAM and their final definition, can be found in the last section in the ESI.†

2.2 Analysis

In order to determine the solution–diffusion permeability given in eqn (1), we need to individually evaluate the partition ratio and the penetrant diffusion inside the network. Note again that we analyze the simulation results in equilibrium, not by measuring any fluxes of the penetrants.

The partition ratio, [scr K, script letter K] = c_in/c₀, was computed by counting and averaging the equilibrium number density of penetrants inside the network and bulk reservoir, as similarly done in our previous works.^36,69,87 We carefully divided the simulation box into three regions (inner membrane, membrane surface, and bulk reservoir) to sample the concentrations without any effects of the surface and the finite membrane width. See Fig. S1 in the ESI† for details.

For the calculation of the penetrant diffusivity in the network, D_in, to avoid the surface effect, we generated 20 auxiliary simulation boxes of diamond unit cells of the polydisperse tetra-functional networks periodically repeated in all three dimensions. We included the penetrants for each parameter set of ε_nn and ε_np and we performed additional simulations of these periodic cells (see Fig. S4 in the ESI†). To determine the cell size and the number of the penetrants in the cell, we used the equilibrium values of the penetrant density and the polymer density obtained from the main simulation data. We computed the mean-squared-displacement (MSD) of the penetrants in the networks, averaged over time and particles,¹²¹ as shown in the ESI,† Fig. S5 (upper panels), within the dimensionless simulation time range from t = 100 to t = 1000 to obtain diffusivity via MSD = 6D_int, ensuring the normal diffusion,¹²¹ which fulfills (see Fig. S5 (lower panels) in the ESI†). This auxiliary setup is fully isotropic and can possibly differ from the network structure in the main simulations, where the membrane deformations can be globally anisotropic. We therefore checked the isotropy of the membrane simulations using several shape descriptors.^122,123 Our examinations with results shown in Fig. S6 in the ESI,† confirm that the polymeric structure is highly isotropic in the central membrane region throughout the whole range of interaction parameters.

Note that in our implicit-solvent simulations in the NpT ensemble, the penetrant's concentration in the bulk reservoir is not a conserved quantity among the systems (as it would be in a grand canonical ensemble), as the penetrants partition differently in the membrane. Therefore, we intentionally employed small penetrant concentrations in the reservoir and only weak interactions (ε_pp = 0.1k_BT) to reduce non-ideality corrections. To check the effects of a varying reservoir concentration, we compare the partitioning results for two different numbers of penetrants (N_p) for all polymer densities at the lowest and highest network–penetrant interaction in Fig. S7 in the ESI.† This comparison indeed reveals a small quantitative difference in the partitioning, but only for large polymer packing fractions and without changing the qualitative picture. Hence, the conclusions based on our generic simulation model should not depend on the ensemble. Actually, our periodic isothermal–isobaric representation (with fixed N) could be adequate to describe permeability experiments using lamellar hydrogel layers.¹¹⁶

We note here that there are other simulation approaches to study the solvent permeability in polymeric networks: Masoud et al. developed DPD (dissipative particle dynamics) CG simulations to investigate the solvent permeability in random networks under external flow.¹⁰³ Wang et al. developed explicit-solvent simulations to study pressure-driven water transport in a polyamide membrane in the isotropic NpT ensemble.¹⁰⁷ A similar all-atom approach was presented by the same authors to study pressure-driven water transport through a CNT (carbon nanotube) membrane pore in the presence of moving walls in the longitudinal direction.¹⁰⁴ This method as an alternative to the grand canonical ensemble¹²⁴ was originally proposed by Ghoufi et al.¹²⁵ and used to study the water transport through a polyamide membrane.¹⁰⁸ These approaches successfully elucidated the permeability of solvents, especially under non-equilibrium situations. However, our approach is different since we focus on the permeability of solutes in an implicit-solvent model in equilibrium, in the absence of pressure gradients. In particular, the chosen solution–diffusion approach has the advantage that effects of partitioning and diffusion on permeability can be studied separately.

3 Results and discussion

3.1 Network density response to solvent quality and penetrants

Six representative simulation snapshots of the system are shown in Fig. 1b for different values of the solvent quality parameter ε_nn and the network–penetrant interaction parameter ε_np. The most swollen state is shown by the lower left snapshot, whereas the most compact state is depicted by the upper right snapshot. The polymer network collapses due to strong network–network attractions ε_nn (poor solvent), otherwise it swells (good solvent). In addition, upon changing the network–penetrant interaction, particularly while using intermediate solvent quality (ε_nn = 0.5k_BT), we note that the larger the attraction ε_np, the more packed is the network (lower volume). This is due to the bridging effects of highly attractive penetrants, contracting the network to maximize favorable interaction contacts.⁶⁹ See Fig. S8 in the ESI† for details of the network volume change depending on the interactions.

In this respect, note that for the same solvent quality βε_nn = 1, the membrane width is increasing with growing network–penetrant interaction βε_np, i.e., there is a penetrant-induced deformation. Recall that we apply isotropic pressure with independent scaling of the box in the x and y directions with essentially vanishing target pressure, hence the membrane can simply freely deform based on its internal virial. For large network–penetrant attractions, the polymer tends to maximize contacts with penetrants (a volume effect) and therefore the membrane deforms to increase the volume-to-interface ratio.

The global effects of the two interaction parameters are summarized in Fig. 1c, which depicts the network volume fraction ϕ_n, i.e., the ratio of the volume occupied by the polymers to the entire network volume, as a function of the solvent quality parameter ε_nn for different values of the network–penetrant interaction ε_np. The network undergoes a typical collapse transition as ε_nn increases at small and intermediate values of ε_np, while the transition becomes more gradual when the penetrants are very attractive (ε_np = 1.5k_BT). This is probably due to local monomer clustering and less structured networks, cf. Fig. S1 (ESI†), indicating a smooth transition. The solid line interpolates between the mean values over the volume fraction averaged over all simulated ε_np values. Also note that the grey data for ϕ_n(βε_np = 0.1) almost coincide with the purple data for ϕ_n(βε_np = 0.7) as shown in Fig. 1c.

3.2 Penetrant partitioning, diffusivity, and permeability

3.2.1 Partitioning. In Fig. 3 we show the partitioning [scr K, script letter K], the penetrant diffusion inside the network D_in, and the permeability [scr P, script letter P], as a function of the network volume fraction (a–c) and the network–penetrant interaction parameter (d–f).

Fig. 3 Simulation results (symbols) and theoretical fits (solid lines). (a) Partitioning [scr K, script letter K] as a function of the network volume fraction for different values of network–penetrant interaction ε_np. The solid lines are the fits from virial expansion eqn (2) (see the text and Table S1 in the ESI† and Fig. 4 for details). (b) Penetrant diffusion in the polymer network D_in as a function of ϕ_n. The solid lines are the fits (also depicted in (e)) from the free-volume approach eqn (4) (see Table S3 in the ESI† for the details). The dashed line is e^−ϕ_n as a reference. (c) Permeability [scr P, script letter P] as a function of ϕ_n. The solid lines are the predictions from eqn (5) using the fitting values from (a) and (b). (d) Partitioning [scr K, script letter K] as a function of ε_np. The solid lines are the prediction according to the virial expression eqn (2) using the B^np₂ and B^nnp₃ values from the fit in (a) and the mean density . The data scale roughly as [scr K, script letter K] ∼ e^7βε_np for dense and attractive polymers (the red dashed line). (e) Penetrant diffusion D_in as a function of ε_np. The solid lines are the fits according to the exponential (Kramers’) scaling D_in/D₀ ∼ e^−βε_np from eqn (4). (f) Permeability [scr P, script letter P] as a function of ε_np. The solid lines are the predictions from eqn (5) using the fitting values from (a) (see also panels (d) and (e)). The cross symbols in each panel are corresponding reference values in the bulk reservoir.

Fig. 4 Virial coefficients (a) B^np₂ and (b) B^nnp₃ in eqn (2) obtained as fitting parameters in Fig. 3a, shown by symbols for different values of ε_np. The exact values of B^LJ₂ and B^LJ₃ for LJ potential are shown by the solid lines for comparison. (c) Fitting parameter values b and c in the free-volume scaling theory eqn (3) (see also Fig. 3b). See the ESI† for details.

The partitioning as a function of the network volume fraction, [scr K, script letter K](ϕ_n), exhibits diverse behavior, ranging over four orders of magnitude depending on the interactions, as shown in Fig. 3a. For low network–penetrant interaction parameters ε_np, [scr K, script letter K] decreases monotonically with increasing network density, since the essentially repelled penetrants are excluded by highly packed polymers (see the second virial coefficient of the LJ system as shown in Fig. 4). For higher values of the LJ potential depth ε_np, the penetrants are increasingly more attracted to the network. The partitioning [scr K, script letter K], however, becomes non-monotonic and reaches a maximum around ϕ_n ≃ 0.3. This partitioning maximization is due to the volume exclusion of the penetrants, which wins over the attraction at high densities.^70,87

The cross-over from penetrant exclusion to enrichment for increasing ε_np at fixed polymer density ϕ_n becomes obvious in Fig. 3d, where we plot [scr K, script letter K](ε_np). At around βε_np ≃ 0.5–0.7 (depending on details of polymer density) the attraction outvalues the steric obstruction and penetrants are on average preferentially adsorbed rather than being in the bulk reservoir, i.e., [scr K, script letter K] > 1. We also observe that partitioning [scr K, script letter K](ε_np) roughly exhibits an exponential increase with a larger slope as ϕ_n increases. The exponential increase of partitioning is also found in ordered membranes,⁸⁷ reflecting that the overall scaling behavior of partitioning (upon changing the interactions) is rather insensitive to the regularity of the network. For dense and attractive polymer networks, we empirically find that [scr K, script letter K] ∼ e^7βε_np, as depicted in Fig. 3d. The prefactor 7 reflects the total mean attraction in the dense systems, where the potential wells of many attractive monomers densely overlap.

In order to gain more theoretical insight and develop an analytical framework for describing the data, we perform a virial expansion of the transfer free energy , and apply it to the partition coefficient [scr K, script letter K] = exp(−βΔG),⁸⁷ as

where B^np₂ is the second virial coefficient, B^nnp₃ is the third virial coefficient, and v₀ = πσ³/6 is the network monomer volume with the diameter σ = σ_nn = σ_np. The expansion eqn (2) is compared with the simulation data by fitting the parameters B^np₂ and B^nnp₃. The final best fits are depicted by the solid curves in Fig. 3a and are in very good agreement. The comparison implies the pronounced contribution of many-body (B^nnp₃) correlations, which are responsible for the non-monotonicity in the attractive and dense regimes.

The fitted B^np₂ and B^nnp₃ parameters can be found in Fig. 4a and b. We find that the second virial coefficients B^np₂ agree well with the values from the explicit relation for LJ particles, cf. the solid line in Fig. 4a. However, as shown in Fig. 4b, the third virial coefficient B^nnp₃ from the fitting deviates from the explicitly computed values of the LJ fluid. This implies that as the polymer density increases, many-body interactions, in particular at the crowded cross-linker regions,¹²⁰ play a major role, which is beyond the correlations of a simple LJ liquid. In fact, the fitted B^nnp₃ values are always positive, i.e., the many-body effect can be identified as, on average, a repulsive contribution.

The data in Fig. 3d are also well described by the virial form, eqn (2), where the solid lines agree with the simulation data. For this, we use eqn (2) with the same virial coefficients obtained from the result shown in Fig. 3a, and assume , which is in fact a good approximation particularly for low and high polymer densities. The dependence of the partitioning on the network volume fraction can thus again be explained by a balance between the network–penetrant attraction and exclusion, which is particularly important for high volume fractions.

3.2.2 Diffusivity. In Fig. 3b, penetrant diffusivity D_in in the network is shown versus the polymer packing fraction ϕ_n. Note that the diffusivity is rescaled by the diffusivity in the bulk reservoir D₀. The diffusivity decreases monotonically and tends to decay rapidly as the network volume fraction increases.^{40,75,76,126,127} The dashed line depicts e^−ϕ_n as a simple reference. We furthermore compare the simulation results with the ‘free-volume’ theory,^{53,72,73,76,128–130}The solid lines show the fitting with the prefactor b and the exponent c, which perform in an excellent fashion. The fitting values of b and c are shown in Fig. 4c. We note that b decays exponentially with ε_np, while c is rather independent. This is physically reasonable if we regard diffusion for large attractions as an activated process, in which the penetrants have to escape from locally bound states (‘traps’). Therefore, here we present a semi-empirical scaling expression for the penetrant diffusivity,

In Fig. 3e we confirm that D_in(ε_np) indeed tends to exponentially decrease. Hence, the Kramers’ type scaling D_in ∝ e^−βε_np for the diffusion limited escape from a single attractive well⁷⁵ fits well, such that our prediction from eqn (4) holds. It is interesting that the energy barrier in the dense systems (i.e., the micro-roughness of the energy landscape) is simply described by ε_np and not by multiples of it, as we observed in the more ordered systems.⁸⁷ Apparently, the random structure (e.g., polydispersity of the network) in thermal fluctuations smoothens out the roughness significantly. Note again that the overall mean attraction (i.e., the mean of the landscape in contrast to its roughness) is much higher than ε_np, since we needed 7ε_np to fit the partition ratio above. We remark here that our semi-empirical scaling law eqn (4) has limitations since it does not behave well when ϕ_n → 0, where D_in/D₀ should go up to unity. However, this highly dilute membrane limit with little influence on transport is not interesting anyway for applications and selectivity-control.

We recall that in the literature there are in fact various proposals for scaling theories of diffusivity, depending on microscopic diffusion mechanisms and length scales.^75,76 In Fig. S9 in the ESI,† we present several appropriate scaling theories for the diffusivity, which are compared to our simulation results. We concluded that for our systems, eqn (4) performs the best throughout the range of ϕ_n, including the dense regimes.

It is interesting that the diffusivity is a simple monotonic function of ϕ_n. In fact, this result is very different from our previous finding for regular topologies, that is, membranes made of a fixed (static) fcc (face-centered-cubic) or simple-cubic lattice of LJ spheres.⁸⁷ There, we found that the diffusivity is rather a complex function of the density of the membranes. We rationalized the effect by the roughness of a potential landscape, which for ordered potential wells on a regular lattice can be a very rapidly changing function of membrane density in certain density regions.⁸⁷ But in the case here, the fluctuations and the polydispersity of the polymer network smoothen out the sharp density effects on the energy landscape and all diffusivities scale similarly exponentially, qualitatively almost independent of parameter ε_np.

3.2.3 Permeability. In Fig. 3c we present permeability [scr P, script letter P] = [scr K, script letter K]D_inversus the packing fraction. The permeability varies by about 4 orders of magnitude in our parameter range. Due to the generic monotonic behavior of the diffusion, the nature of the permeability depending on the volume fraction reflects essentially the one of the partition ratio [scr K, script letter K], while the diffusivity only quantitatively scales the results. Hence, we find that for small interactions ε_np, the permeability monotonically decreases with density, whereas for stronger interactions, it becomes a non-monotonic function of density. Therefore, as an important finding, the permeability can be maximized in our network model system. For the largest network–penetrant attraction, the permeability is maximized at around ϕ_n ≃ 0.28 by a factor of around 20 when compared to the bulk reservoir reference permeability [scr P, script letter P] = D₀ (the cross symbol).

Having well-performing scaling laws for [scr K, script letter K] and D_in from eqn (2) and (4), now we attempt to empirically construct a scaling law for the permeability, via their product, eqn (1),

which comprises the attractive contribution as a function of the network–penetrant interaction ε_np, and the exclusion contribution as a function of the packing fraction. The maximization of [scr P, script letter P] can therefore be understood viaeqn (2) and (4). The solid lines in Fig. 3c are the predictions from eqn (5) using the fitting parameters determined already in panels a and b, showing very good agreement with the simulation results.

The permeability as a function of the network–penetrant interaction, [scr P, script letter P](ε_np), as shown in Fig. 3f, is an increasing function from the global minimum at around βε_np = 0.1, which substantially depends on membrane density. Here, the selective tuning of [scr P, script letter P] is mainly controlled by the penetrant's excluded volume. The prediction from the empirical scaling eqn (5) indeed agrees well with the simulation data, in particular, capturing the competition and cancellation between the exponentially growing partitioning and the exponentially decreasing diffusion.

3.3 Correlations between [scr K, script letter K] and D_in and tuning of permselectivity

The diagram in Fig. 5 plots the partitioning [scr K, script letter K] versus the diffusivity D_in and thus presents a landscape visualizing how they are correlated, i.e., a partitioning–diffusivity correlation diagram. The plot shows a wide landscape of the permeability spanning over several orders of magnitude. The black dashed line depicts the iso-permeability line of the bulk reservoir permeability [scr P, script letter P]/D₀ = 1, where the two contributions exactly cancel out. The data in Fig. 5 at low and intermediate polymer densities lead to final permeabilities close to the iso-permeability line, hence exhibiting clear anti-correlations and cancellations. Such cancellations were also observed, even more massively, in membranes constructed of static regular obstacles.⁸⁷

Fig. 5 Partitioning–diffusion [scr K, script letter K]–D_in/D₀ correlation diagram. As depicted in the legend, symbols of the same color have the same network-penetrant interaction, but different polymer densities, i.e., the network-network interaction and hence the polymer volume fraction increases (for each color individually) from right (high diffusion) to left (low diffusion). The black dashed line depicts the iso-permeability line [scr P, script letter P]/D₀ = [scr K, script letter K]D_in/D₀ = 1 (reference bulk reservoir permeability), where the actions of [scr K, script letter K] and D_in on the membrane permeability exactly cancel each other. The arrow connects two states II (βε_np = 1.2) and I (βε_np = 0.7) at packing fraction ϕ_n ≃ 0.17, featuring the selectivity ratio α_II,I ≡ [scr P, script letter P]_II/[scr P, script letter P]_I ≃ 6 (see the text).

The anti-correlations can be understood by glancing back at our empirical scaling law for permeability, eqn (5). The attraction between monomers and penetrants increases the uptake of penetrants in the membrane roughly exponentially. However, at the same time the attraction enhances the microscopic roughness and deepens local traps, thereby impeding the thermally activated (Kramers) escape, which in turn also leads to an exponential decrease of diffusion. In many regimes, these two effects cancel out, but the exact behavior depends on the details of the variation of the energy landscape.⁸⁷

In fact, the data do not simply follow one universal trend with some statistical scatter. On the contrary, we observe a set of distinct non-monotonic trends, each of which has its own physical explanation. For example, for the most repulsive penetrants (gray upside triangles) the data diverge substantially downwards, as both partitioning and diffusion decrease with polymer density due to excluded-volume (steric) constraints. In contrast, for the most attractive penetrants (red circles) the permeability stays closer to the bulk permeability [scr P, script letter P]/D₀ = 1 but shows a distinct maximum. The reason for this is that partitioning first increases stronger (beyond unity) than diffusion decreases for increasing polymer density, but then steric exclusion limits the uptake. Diffusion apparently does not decrease as much as partitioning increases because in our fluctuating, polydisperse network the energy landscape for diffusion stays relatively smooth and independent of density, in contrast to highly ordered membranes.⁸⁷ This knowledge can be harvested to tune and optimize the selectivity of a polymer membrane. In contrast to the ordered membranes,⁸⁷ however, this work indicates that the diffusivity in fluctuating, disordered networks only rescales the permeability, while the shape of the dependence is dictated by the partitioning behavior. Hence, membrane design for polydisperse networks should focus more on the behavior of solubility in the network than on the mobility.

Hence, the diagram in Fig. 5 presents non-trivial pathways of the permeability [scr P, script letter P] along the two variable parameters, density and network–penetrant attraction. This clearly shows how the permeability can be tuned substantially over several orders of magnitude already by a relatively small material parameter space. With this, a significant selective permeability (permselectivity) can be demonstrated depending on the interaction parameter, ε_np (which in reality is different for various chemically specific penetrants). For instance, defined as α_II,I ≡ [scr P, script letter P]_II/[scr P, script letter P]_I,⁴² the selectivity for the states II (βε_np = 1.2) and I (βε_np = 0.7) depicted by the arrow in Fig. 5 at a packing fraction ϕ_n ≃ 0.17 amounts to α_II,I ≈ 6, which is large. Hence, a small difference in interactions of half a k_BT results already in a permeability ratio of almost one order of magnitude.

We note here that there are various mechanisms that govern the selectivity. We limit our study to varying the interaction strengths (ε_nn and ε_np) at fixed penetrant size (σ). Varying the size, in fact, leads to even more complex behavior, such as molecular sieving, once the penetrant size becomes larger than the mesh size of the network.^109,110 In this work, we focus on the interaction-dependent permselectivity, where the sieving can be effectively tuned by changing the interactions.

4 Conclusion

We presented extensive (implicit-solvent) coarse-grained simulations and scaling theories for penetrant transport through semi-flexible, cross-linked, and polydisperse polymer networks with a focus on the linear-response permeability, calculated by the equilibrium partitioning and diffusion of the penetrants inside the network. The permeability has been found to be largely tunable by varying the polymer network density and the microscopic interactions between the network and the diffusive penetrants. In particular, significant maximization and minimization of the permeability were found, fine-tuned by the solvent quality and the network–penetrant interactions. The results were rationalized by phenomenological scaling theories that we developed in this work, which include a virial expansion with two-body attractions and many-body exclusion effects for the partitioning, and a combination of the free-volume and Kramers’ escape scaling laws for the diffusivity. Our presented laws, despite their simplicity, capture salient features of the system, showing good agreement with the simulation results.

The penetrant diffusivity turned out to be a rather smooth function of the network density, implying substantial effects of the fluctuation and randomness of the polymer network. The highly fluctuating and polydisperse nature of the network averages out the roughness of the energy landscape, which was more pronounced and sensitive to parameter changes in highly ordered, lattice-based and static membrane systems.⁸⁷ However, it is not clear whether smoothing-out is more a consequence of polydispersity, or just a consequence of the irregular distribution of polymer segments in space, which might also occur in monodisperse fluctuating networks that may be studied in the future.

Nevertheless, permeability revealed a rather intricate, non-monotonic behavior over several orders of magnitude, originating from the complex nature of partitioning, while quantitatively and substantially modified by the anti-correlated and canceling contributions of the diffusion. As a consequence, only small changes of interactions, e.g., by half a k_BT can already modify the selectivity of the membrane by a factor of 6. Our study provides a further step in the fundamental understanding and development of a minimal theory to characterize better the permeability in flexible and fluctuating polymer-based membrane systems.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank Matthias Ballauff, Benjamin Rotenberg, Arturo Moncho-Jordá and Changbong Hyeon for fruitful 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 No. 646659). W. K. K. acknowledges the support by a KIAS Individual Grant (CG076001) at Korea Institute for Advanced Study. M. K. acknowledges the financial support from the Slovenian Research Agency (research core funding No. P1-0055). The simulations were performed with resources provided by the North-German Supercomputing Alliance (HLRN). We acknowledge the Center for Advanced Computation at Korea Institute for Advanced Study for providing computing resources for this work.

References

1. T. Graham, Philos. Mag., 1866, 32, 401–420.
2. A. Finkelstein, Current Topics in Membranes and Transport, Elsevier, 1984, vol. 21, pp. 295–308.
3. Q. Al-Awqati, Nat. Cell Biol., 1999, 1, E201.
4. K. P. Lee, T. C. Arnot and D. Mattia, J. Membr. Sci., 2011, 370, 1–22.
5. R. M. Venable, A. Krämer and R. W. Pastor, Chem. Rev., 2019, 119, 5954–5997.
6. D. Shasby, S. Shasby, J. Sullivan and M. Peach, Circ. Res., 1982, 51, 657–661.
7. J. Wingender, T. R. Neu and H.-C. Flemming, Microbial extracellular polymeric substances, Springer, 1999, pp. 1–19.
8. E. D. Hay, Cell biology of extracellular matrix, Springer Science & Business Media, 2013.
9. J. Witten, K. Ribbeck, V. V. Khutoryanskiy, L. Wu, M. Liu, W. Shan, X. Zhu, L. Li, Z. Zhang and Y. Huang, Nanoscale, 2017, 9, 8080.
10. C. P. Goodrich, M. P. Brenner and K. Ribbeck, Nat. Commun., 2018, 9, 4348.
11. G. Fuhrmann, Nat. Nanotechnol., 2020, 15, 168–169.
12. J. Taipale and J. Keski-Oja, FASEB J., 1997, 11, 51–59.
13. C. J. Dowd, C. L. Cooney and M. A. Nugent, J. Biol. Chem., 1999, 274, 5236–5244.
14. E. W. Raines, Int. J. Clin. Exp. Pathol., 2000, 81, 173–182.
15. A. M. Garcia, N. Szasz, S. B. Trippel, T. I. Morales, A. J. Grodzinsky and E. H. Frank, Arch. Biochem. Biophys., 2003, 415, 69–79.
16. R. G. Thorne, A. Lakkaraju, E. Rodriguez-Boulan and C. Nicholson, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 8416–8421.
17. L. Zhang, B. S. Gardiner, D. W. Smith, P. Pivonka and A. J. Grodzinsky, J. Theor. Biol., 2010, 263, 20–29.
18. A. D. Theocharis, S. S. Skandalis, C. Gialeli and N. K. Karamanos, Adv. Drug Delivery Rev., 2016, 97, 4–27.
19. A. Halperin, M. Kröger and F. M. Winnik, Angew. Chem., Int. Ed., 2015, 54, 15342–15367.
20. D. M. Vriezema, M. Comellas Aragonès, J. A. A. W. Elemans, J. J. L. M. Cornelissen, A. E. Rowan and R. J. M. Nolte, Chem. Rev., 2005, 105, 1445–1490.
21. S. Carregal-Romero, N. J. Buurma, J. Pérez-Juste, L. M. Liz-Marzán and P. Hervés, Chem. Mater., 2010, 22, 3051–3059.
22. M. A. C. Stuart, W. T. S. 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., 2010, 9, 101–113.
23. Y. Lu and M. Ballauff, Prog. Polym. Sci., 2011, 36, 767–792.
24. K. Renggli, P. Baumann, K. Langowska, O. Onaca, N. Bruns and W. Meier, Adv. Funct. Mater., 2011, 21, 1241–1259.
25. P. Tanner, P. Baumann, R. Enea, O. Onaca, C. Palivan and W. Meier, Acc. Chem. Res., 2011, 44, 1039–1049.
26. Y. Guan and Y. Zhang, Soft Matter, 2011, 7, 6375–6384.
27. P. Hervés, M. Pérez-Lorenzo, L. M. Liz-Marzán, J. Dzubiella, Y. Lu and M. Ballauff, Chem. Soc. Rev., 2012, 41, 5577–5587.
28. S. Wu, J. Dzubiella, J. Kaiser, M. Drechsler, X. Guo, M. Ballauff and Y. Lu, Angew. Chem., Int. Ed., 2012, 51, 2229–2233.
29. J. Gaitzsch, X. Huang and B. Voit, Chem. Rev., 2015, 116, 1053–1093.
30. S. Campisi, M. Schiavoni, C. Chan-Thaw and A. Villa, Catalysts, 2016, 6, 185.
31. G. Prieto, H. Tüysüz, N. Duyckaerts, J. Knossalla, G.-H. Wang and F. Schüth, Chem. Rev., 2016, 116, 14056–14119.
32. S. H. Petrosko, R. Johnson, H. White and C. A. Mirkin, J. Am. Chem. Soc., 2016, 138, 7443–7445.
33. H. Jia, R. Roa, S. Angioletti-Uberti, K. Henzler, A. Ott, X. Lin, J. Möser, Z. Kochovski, A. Schnegg, J. Dzubiella, M. Ballauff and Y. Lu, J. Mater. Chem. A, 2016, 4, 9677–9684.
34. S. Angioletti-Uberti, Y. Lu, M. Ballauff and J. Dzubiella, J. Phys. Chem. C, 2015, 119, 15723–15730.
35. R. Roa, W. K. Kim, M. Kanduč, J. Dzubiella and S. Angioletti-Uberti, ACS Catal., 2017, 7, 5604–5611.
36. M. Kanduč, W. K. Kim, R. Roa and J. Dzubiella, Mol. Syst. Des. Eng., 2020, 5, 602–619.
37. L. M. Robeson, J. Membr. Sci., 1991, 62, 165–185.
38. P. Pandey and R. Chauhan, Prog. Polym. Sci., 2001, 26, 853–893.
39. E. Atci, I. Erucar and S. Keskin, J. Phys. Chem. C, 2011, 115, 6833–6840.
40. K. Falk, B. Coasne, R. Pellenq, F.-J. Ulm and L. Bocquet, Nat. Commun., 2015, 6, 6949.
41. A. Obliger, R. Pellenq, F.-J. Ulm and B. Coasne, J. Phys. Chem. Lett., 2016, 7, 3712–3717.
42. B. D. Freeman, Macromolecules, 1999, 32, 375–380.
43. H. B. Park, J. Kamcev, L. M. Robeson, M. Elimelech and B. D. Freeman, Science, 2017, 356, 1137.
44. M. A. Shannon, P. W. Bohn, M. Elimelech, J. G. Georgiadis, B. J. Mariñas and A. M. Mayes, Nature, 2008, 452, 301–310.
45. G. M. Geise, H.-S. Lee, D. J. Miller, B. D. Freeman, J. E. McGrath and D. R. Paul, J. Polym. Sci., Part B: Polym. Phys., 2010, 48, 1685–1718.
46. G. M. Geise, H. B. Park, A. C. Sagle, B. D. Freeman and J. E. McGrath, J. Membr. Sci., 2011, 369, 130–138.
47. D. Menne, F. Pitsch, J. E. Wong, A. Pich and M. Wessling, Angew. Chem., Int. Ed., 2014, 53, 5706–5710.
48. B. Tansel, J. Sager, T. Rector, J. Garland, R. F. Strayer, L. Levine, M. Robert, M. Hummerick and J. Bauer, Sep. Purif. Technol., 2006, 51, 40–47.
49. Z. Tan, S. Chen, X. Peng, L. Zhang and C. Gao, Science, 2018, 360, 518–521.
50. W. Hyk and K. Kitka, J. Environ. Chem. Eng., 2018, 6, 6108–6117.
51. C. S. Brazel and N. A. Peppas, Polymer, 1999, 40, 3383–3398.
52. D. F. Stamatialis, B. J. Papenburg, M. Girones, S. Saiful, S. N. Bettahalli, S. Schmitmeier and M. Wessling, J. Membr. Sci., 2008, 308, 1–34.
53. H. Yasuda, A. Peterlin, C. Colton, K. Smith and E. Merrill, Die Makromol. Chem., 1969, 126, 177–186.
54. D. R. Paul, Sep. Purif. Methods, 1976, 5, 33–50.
55. J. Williams and R. W. Baker, J. Membr. Sci., 1995, 107, 1–21.
56. S. Gehrke, J. Fisher, M. Palasis and M. E. Lund, Ann. N. Y. Acad. Sci., 1997, 831, 179–207.
57. S. C. George and S. Thomas, Prog. Polym. Sci., 2001, 26, 985–1017.
58. M. Ulbricht, Polymer, 2006, 47, 2217–2262.
59. A. Missner and P. Pohl, Chem. Phys. Chem., 2009, 10, 1405–1414.
60. R. W. Baker and B. T. Low, Macromolecules, 2014, 47, 6999–7013.
61. P. Nelson, Biological physics, WH Freeman, New York, 2004.
62. R. Sander, Atmos. Chem. Phys., 2015, 15, 4399–4981.
63. R. W. Baker, Membrane Technology and Applications, John Wiley & Sons, Ltd, 2004.
64. A. Obliger, M. Jardat, D. Coelho, S. Bekri and B. Rotenberg, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 043013.
65. A. Moncho-Jordá and I. Adroher-Bentez, Soft Matter, 2014, 10, 5810–5823.
66. I. Adroher-Bentez, S. Ahualli, A. Martn-Molina, M. Quesada-Pérez and A. Moncho-Jordá, Macromolecules, 2015, 48, 4645–4656.
67. A. Erbas and M. Olvera de la Cruz, Macromolecules, 2016, 49, 9026–9034.
68. B. Rotenberg, J.-F. Dufreche, B. Bagchi, E. Giffaut, J.-P. Hansen and P. Turq, J. Chem. Phys., 2006, 124, 154701.
69. W. K. Kim, A. Moncho-Jordá, R. Roa, M. Kanduč and J. Dzubiella, Macromolecules, 2017, 50, 6227–6237.
70. L. Pérez-Mas, A. Martín-Molina, M. Quesada-Pérez and A. Moncho-Jordá, Phys. Chem. Chem. Phys., 2018, 20, 2814–2825.
71. M. Kanduč, W. K. Kim, R. Roa and J. Dzubiella, ACS Nano, 2019, 13, 11224–11234.
72. H. Yasuda, C. Lamaze and L. D. Ikenberry, Die Makromol. Chem., 1968, 118, 19–35.
73. H. Yasuda, L. Ikenberry and C. Lamaze, Die Makromol. Chem., 1969, 125, 108–118.
74. I. C. Kim and S. Torquato, J. Chem. Phys., 1992, 96, 1498–1503.
75. L. Masaro and X. Zhu, Prog. Polym. Sci., 1999, 24, 731–775.
76. B. Amsden, Macromolecules, 1998, 31, 8382–8395.
77. A. P. Chatterjee, J. Phys.: Condens. Matter, 2011, 23, 375103.
78. Y. Jiao and S. Torquato, Phys. Biol., 2012, 9, 036009.
79. M. Spanner, S. K. Schnyder, F. Höfling, T. Voigtmann and T. Franosch, Soft Matter, 2013, 9, 1604–1611.
80. A. Godec, M. Bauer and R. Metzler, New J. Phys., 2014, 16, 092002.
81. H. Liasneuski, D. Hlushkou, S. Khirevich, A. Höltzel, U. Tallarek and S. Torquato, J. Appl. Phys., 2014, 116, 034904.
82. X. Zhang, J. Hansing, R. R. Netz and J. E. DeRouchey, Biophys. J., 2015, 108, 530–539.
83. J. Hansing, C. Ciemer, W. K. Kim, X. Zhang, J. E. DeRouchey and R. R. Netz, Eur. Phys. J. E: Soft Matter Biol. Phys., 2016, 39, 53.
84. J. Hansing, J. R. Duke III, E. B. Fryman, J. E. DeRouchey and R. R. Netz, Nano Lett., 2018, 18, 5248–5256.
85. J. Hansing and R. R. Netz, Macromolecules, 2018, 51, 7608–7620.
86. J. Hansing and R. R. Netz, Biophys. J., 2018, 114, 2653–2664.
87. W. K. Kim, M. Kanduč, R. Roa and J. Dzubiella, Phys. Rev. Lett., 2019, 122, 108001.
88. P. Jha, J. Zwanikken, F. Detcheverry, J. de Pablo and M. Olvera de la Cruz, Soft Matter, 2011, 7, 5965–5975.
89. M. Quesada-Pérez, J. Ramos, J. Forcada and A. Martín-Molina, J. Chem. Phys., 2012, 136, 244903.
90. P. Košovan, T. Richter and C. Holm, Macromolecules, 2015, 48, 7698–7708.
91. H. Kobayashi and R. G. Winkler, Sci. Rep., 2016, 6, 19836.
92. A. Schmid, J. Dubbert, A. A. Rudov, J. Pedersen, P. Lindner, M. Karg, I. Potemkin and W. Richtering, Sci. Rep., 2016, 22736.
93. I. Zadok and S. Srebnik, J. Phys. Chem. B, 2018, 122, 7091–7101.
94. M. Kanduč, W. K. Kim, R. Roa and J. Dzubiella, Macromolecules, 2018, 51, 4853–4864.
95. H. Takeuchi and K. Okazaki, J. Chem. Phys., 1990, 92, 5643–5652.
96. J. Sonnenburg, J. Gao and J. Weiner, Macromolecules, 1990, 23, 4653–4657.
97. F. Müller-Plathe, J. Chem. Phys., 1991, 94, 3192–3199.
98. R. Sok, H. Berendsen and W. Van Gunsteren, J. Chem. Phys., 1992, 96, 4699–4704.
99. F. Müller-Plathe, Acta Polym., 1994, 45, 259–293.
100. P. A. Netz and T. Dorfmüller, J. Chem. Phys., 1997, 107, 9221–9233.
101. D. Hofmann, L. Fritz, J. Ulbrich, C. Schepers and M. Böhning, Macromol. Theory Simul., 2000, 9, 293–327.
102. H. Zhou and S. B. Chen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 021801.
103. H. Masoud and A. Alexeev, Macromolecules, 2010, 43, 10117–10122.
104. L. Wang, R. S. Dumont and J. M. Dickson, J. Chem. Phys., 2013, 138, 124701.
105. H. Li, A. Erbas, J. Zwanikken and M. Olvera de la Cruz, Macromolecules, 2016, 49, 9239–9246.
106. D. Sandrin, D. Wagner, C. Sitta, R. Thoma, S. Felekyan, H. Hermes, C. Janiak, N. de Sousa Amadeu, R. Kühnemuth and H. Löwen, et al. , Phys. Chem. Chem. Phys., 2016, 18, 12860–12876.
107. L. Wang, R. S. Dumont and J. M. Dickson, RSC Adv., 2016, 6, 63586–63596.
108. M. Ding, A. Szymczyk and A. Ghoufi, J. Membr. Sci., 2016, 501, 248–253.
109. K. Zhang and S. K. Kumar, ACS Macro Lett., 2017, 6, 864–868.
110. K. Zhang, D. Meng, F. Müller-Plathe and S. K. Kumar, Soft Matter, 2018, 14, 440–447.
111. P. Higgs and R. Ball, J. Phys. France, 1988, 49, 1785–1811.
112. E. Geissler, F. Horkay and A.-M. Hecht, Phys. Rev. Lett., 1993, 71, 645.
113. G. Glatting, R. Winkler and P. Reineker, Macromolecules, 1995, 28, 5906–5909.
114. J. S. Soares and P. Zunino, Biomaterials, 2010, 31, 3032–3042.
115. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
116. Y. F. Yue, M. A. Haque, T. Kurokawa, T. Nakajima and J. P. Gong, Adv. Mater., 2013, 25, 3106–3110.
117. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
118. A. Erbas and M. Olvera de la Cruz, ACS Macro Lett., 2015, 4, 857–861.
119. J. Heyda, A. Muzdalo and J. Dzubiella, Macromolecules, 2013, 46, 1231–1238.
120. S. Milster, R. Chudoba, M. Kanduč and J. Dzubiella, Phys. Chem. Chem. Phys., 2019, 21, 6588–6599.
121. J. Shin, A. G. Cherstvy, W. K. Kim and V. Zaburdaev, Phys. Chem. Chem. Phys., 2017, 19, 18338–18347.
122. D. N. Theodorou and U. W. Suter, Macromolecules, 1985, 18, 1206–1214.
123. P. Fu and Y. F. Dafalias, Int. J. Solids Struct., 2015, 63, 68–81.
124. H. Eslami, F. Mozaffari, J. Moghadasi and F. Müller-Plathe, J. Chem. Phys., 2008, 129, 194702.
125. A. Ghoufi, D. Morineau, R. Lefort, I. Hureau, L. Hennous, H. Zhu, A. Szymczyk, P. Malfreyt and G. Maurin, J. Chem. Phys., 2011, 134, 074104.
126. J. Haus and K. Kehr, Phys. Rep., 1987, 150, 263–406.
127. S. K. Ghosh, A. G. Cherstvy and R. Metzler, Phys. Chem. Chem. Phys., 2014, 17, 1847–1858.
128. N. A. Peppas and C. T. Reinhart, J. Membr. Sci., 1983, 15, 275–287.
129. C. T. Reinhart and N. A. Peppas, J. Membr. Sci., 1984, 18, 227–239.
130. S. R. Lustig and N. A. Peppas, J. Appl. Polym. Sci., 1988, 36, 735–747.
131. P. Kumar, L. Theeyancheri, S. Chaki and R. Chakrabarti, Soft Matter, 2019, 15, 8992–9002.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sm01083a

---