Modeling Gel Swelling Equilibrium in Mean-Field: From explicit Models to Poisson-Boltzmann

TL;DR

This paper introduces a double mean-field model for charged macrogels that accurately predicts swelling behavior by combining explicit polymer sampling with a Poisson-Boltzmann approximation, extending to pH-sensitive gels.

Contribution

It develops a computationally efficient double mean-field approach that bridges explicit polymer models and Poisson-Boltzmann theory for swelling predictions.

Findings

Accurate swelling predictions with explicit coarse-grained models.

Excellent agreement with molecular dynamics simulations.

Effective extension to pH-sensitive gel systems.

Abstract

We develop a double mean-field theory for charged macrogels immersed in electrolyte solutions in the spirit of the cell model approach. We first demonstrate that the equilibrium sampling of a single explicit coarse-grained charged polymer in a cell yields accurate predictions of the swelling equilibrium if the geometry is suitably chosen and all pressure contributions have been incorporated accurately. We then replace the explicit flexible chain by a suitably modeled penetrable charged rod that allows to compute all pressure terms within the Poisson-Boltzmann approximation. This model, albeit computationally cheap, yields excellent predictions of swelling equilibria under varying chain length, polymer charge fraction, and external reservoir salt concentrations when compared to coarse-grained molecular dynamics simulations of charged macrogels. We present an extension of the model to the experimentally relevant cases of pH-sensitive gels.

Figures

Equations

$$P_{\mathrm{in}}=\frac{1}{3}{P}_{\mathrm{cap}}+\frac{2}{3}{P}_{\mathrm{side}},$$

$${\Pi }_{(z,z)}=\frac{{\sum }_i{m}_i{v}_i^{(z)}{v}_i^{(z)}}{V}+\frac{{\sum }_{j>i}{\vec{F}}_{ij}^{(z)}\cdot {\vec{r}}_{ij}^{(z)}}{V}+{\Pi }_{(z,z)}^{\text{Coulomb}}.$$

$${\nabla }^2\psi =-\frac{1}{{ϵ}_0{ϵ}_{\mathrm{r}}}\left(\sum _i{q}_i{c}_i(\vec{r})+{\rho }_f(\vec{r})\right).$$

$$c_i(r)=c_i^{\mathrm{b}}\exp \left(-\frac{q_i\psi (r)}{k_{\mathrm{B}}T}\right).$$

$$P_{\mathrm{side}}$$

$$P_{\mathrm{cap}}^{\text{ions}}$$

$$\begin{array}{rl}& P_{\mathrm{cap}}^{\mathrm{str}}(R_{\mathrm{e}})=P^{\mathrm{conf}}-{\sigma }^{\mathrm{chain}}\\ & =\frac{1}{\pi R_{\mathrm{out}}^2}\left(\frac{k_{B}T}{b}\frac{R_0^3}{R_{\mathrm{e}}^3}{\mathcal{L}}^{-1}\left(\frac{R_0}{R_{\max }}\right)-\frac{k_{B}T}{b}{\mathcal{L}}^{-1}\left(\frac{R_{\mathrm{e}}}{R_{\max }}\right)\right),\end{array}$$

$$c(\text{ce}A-,\vec{r})=\frac{c_0(\vec{r})K_{\mathrm{a}}}{c^{\mathrm{b}}(\text{ce}H+)\exp (-e_0\psi (\vec{r})/(k_{\mathrm{B}}T))+K_{\mathrm{a}}},$$

$$c^{\mathrm{b}}(\text{ce}H+)+c^{\mathrm{b}}(\text{ce}Na+)=c^{\mathrm{b}}(\text{ce}Cl-)$$

Full text

Polyelectrolyte gels consist of crosslinked charged polymers (polyelectrolytes) that can be synthesized with various topologies and are produced in sizes ranging from nanometers (nanogels) up to centimeters (macrogels) [1, 2]. They show a large, reversible uptake of water that is exploited in numerous daily-life products, such as in superabsorbers, cosmetics, pharmaceuticals [3, 4, 5], agriculture [6, 7], or quite recently water desalination [8, 9]. Tailoring polyelectrolyte gels to their applications requires a sufficiently accurate prediction of their swelling capabilities and elastic responses, a task that still goes beyond analytical approaches [10, 11, 12, 13, 14, 15, 16, 17, 18]. So far only all-atom simulations of short single chains in the bulk (not of whole hydrogels) with explicit water have been performed [19, 20, 21, 22]. On the other hand, coarse-grained polyelectrolyte network models have demonstrated their ability to amend analytical approaches, showing that structural microscopic details can have noticeable effects on the macroscopic properties such as the swelling [23, 24, 25, 26, 27, 28, 29, 30, 31, 32]. Macroscopic gels with monodisperse chain length can be simulated with microscopic detail using molecular dynamics (MD) simulations with periodic boundary conditions (PBCs) (cf. periodic gel model) where a unit gel section is connected periodically to yield an infinite gel without boundaries. However, even MD simulations of periodic gels remain computationally very expensive due to the many particles and the slow relaxation times of the involved polymers. Thus, the development of computationally efficient mean-field models capable of predicting swelling equilibria have been of scientific interest in the last years [33, 15, 31, 32]. First ideas of using a Poisson-Boltzmann (PB) cell model under tension were put forward by Mann for salt-free gels, with moderate success [33].

About sixty years ago Katchalsky and Michaeli [11] suggested a free energy model that has recently been shown to predict swelling equilibria reasonably well [31] when compared to MD simulations of charged bead-spring gels. This model has been applied to explore a wide parameter space in search of optimal desalination conditions [9]. However, the Katchalsky model fails [31] for Manning parameters $\xi=\lambda_{B}/\langle d\rangle>1$ [34], where $\lambda_{B}$ denotes the Bjerrum length, and $\langle d\rangle$ the average distance between polymer backbone charges. This is presumably due to the usage of the Debye-Hückel approximation.

In this letter we describe two successive mean-field approaches to render the determination of swelling equilibria of polyelectrolytes accurately and efficiently. Figure 1 displays our construction scheme of the two different models. First, we describe a single-chain MD cell model, that reproduces results similar to those obtained from expensive MD simulations of multiple crosslinked chains. This reduces the many-body problem of the macroscopic gel to one of computing the pressure exerted within a cell containing a single polyelectrolyte chain under varying environmental conditions. The single-chain cell model can thus be viewed as a mean-field attempt to factorize the many-body partition function of the macrogel into a product state of suitable identical subunits [35]. We then show that the single-chain cell model can further be simplified in a second mean-field step using a PB description of the chain with appropriate boundary conditions. The PB cell description has been successful in describing a variety of polyelectrolyte phenomena [36, 37, 38, 39, 40, 41] and is here applied to macroscopic polyelectrolyte gels for the first time. The quality of our two mean-field models is gauged by comparing them to 60 data points for the swelling equilibrium of periodic monodisperse gel MD simulations obtained within a wide range of system parameters. We want to emphasize that none of our models assumes a specific stretching state of the chains: they are constructed to incorporate the main physical effects which happen during stretching (at high chain extensions) and compression of a polyelectrolyte gel (at low chain extensions).

Finally, we generalize the PB cell model to account for the effect of weak groups along the chain backbone. From this we can efficiently predict swelling equilibria for the experimentally relevant cases of weak polyelectrolyte gels.

We specifically compare our results to MD data for the periodic gel model obtained by Košovan et. al [31] where a continuum solvent is used with standard charged bead-spring polymers connected in a diamond lattice as in the work of Ref. [27] together with explicit salt and counterions. Here a perfect tetrafunctional gel is described by the chain length $N$ and the monomers charge fraction $f$. This explicit particle based model uses monodisperse chain lengths which is in contrast to the heterogeneity observed in a wide variety of synthesized gels [42]. It would be computationally very costly to introduce chain length heterogeneity into this model since it would require to simulate a much larger representative volume element. However, these periodic gel simulations of a monodisperse gel are sufficient to test the validity of our two consecutive mean-field approaches.

All MD simulations (namely the periodic gel and the single-chain model described later) are performed with PBCs using the simulation package ESPResSo [43, 44]. All particles interact via WCA interactions [45, 46]. Monomers are connected via FENE bonds (including the ends of the single periodic chain) with Kremer-Grest parameters [47]. We employ the Langevin thermostat [48], and all electrostatic interactions between particles are calculated with the P3M method [49] tuned to an accuracy in the root mean squared error of the electrostatic force of at least $10^{-3}$ in electrostatic simulation units [50].

In addition, salt ion pair exchanges between the simulation volume and an external reservoir are performed using grand canonical Monte Carlo moves [48]. The equilibrium pressure inside the gel and the electro-chemical potentials of all $i$ species balance out with that of the reservoir: $P_{\mathrm{in}}(V_{\mathrm{eq}})=P_{\mathrm{res}}$ and $\mu^{\mathrm{gel}}_{i}=\mu^{\mathrm{res}}_{i},$ respectively. The simulations are performed at different imposed volumes and the internal pressure is measured after chemical equilibrium is reached. The reservoir pressure is approximated by the ideal gas expression $P_{\mathrm{res}}=\sum_{i}k_{\mathrm{B}}Tc_{i}^{\mathrm{b}}$, with the Boltzmann constant $k_{\mathrm{B}}$, temperature $T$ and bulk ion concentrations $c_{i}^{\mathrm{b}}$, which are chosen such that the bulk is electroneutral.

The above approach requires multiple simulations at different imposed volumes until the condition $P_{\mathrm{in}}(V_{\mathrm{eq}})=P_{\mathrm{res}}$ can be narrowed down to a satisfactory small interval. The equilbrium volume is found at the intersection of $P_{\mathrm{in}}(V_{\mathrm{eq}})$ and $P_{\mathrm{res}}$ by using linear interpolation. The errorbar is given by the width of the interval. Under equilibrium conditions, the end-to-end distance $R_{\mathrm{e}}$ is equal to the equilibrium chain extension $R_{\mathrm{eq}}$.

The first model for complexity reduction is the mean-field single-chain model. Like in the cylindrical cell model used to describe solutions of polyelectrolytes [38, 40, 41], we propose constructing the many-body partition function of a periodic gel as a suitable product of individual cylindrical cells containing a single polyelectrolyte chain with added salt. Since the main physical principle is the balance between the polyelectrolyte chain tension and the remaining pressure contributions (mainly the ionic ones), these cylindrical cells have an axial length chosen such as to represent the polymer chain extension between gel crosslinks. Like in the periodic gel model, we perform MD simulations for a single-chain in cylindrical confinement allowing explicit salt ion pairs to enter the cell volume and reach chemical equilibrium with an external reservoir [48]. For a perfect affine compression of a (fully stretched) tetrafunctional gel (built in a diamond cubic lattice) the volume per chain is given by: $V_{\text{chain}}=R_{\mathrm{e}}^{3}/A$, with the geometrical prefactor $A=\sqrt{27}/4$ [31]. Using the volume of the cylindrical cell as the volume per chain, we arrive at a constant aspect ratio $R_{\mathrm{out}}/R_{\mathrm{e}}=1/\sqrt{\pi A}\approx 0.49$, where $R_{\mathrm{out}}$ denotes the radius of the cylindrical cell, see Fig. 1. This is in contrast to the deformation of a pure isolated chain, which does not occupy a cylindrical volume of constant aspect ratio upon a stretching deformation. The simulation volume is completely defined by the length of the cylinder, or equivalently $R_{\mathrm{e}}$ 111In the single chain MD simulations the cylinder height is $R_{\mathrm{e}}+b$ (with the average bond length $b\approx 0.966\sigma$) in order to be able to satisfy PBCs.. Note that the single chain under confinement sees its images in the axial direction (due to PBCs) which in a simplified way mimics the electrostatic environment in a gel, where the end of a single chain sees the next chain. For cylindrical geometries under affine compression, the pressure inside the volume is given by [41, 51]:

[TABLE]

where the total internal pressure $P_{\mathrm{in}}$ is split into the two contributions from, $P_{\mathrm{cap}}$, the cylinder end caps and, $P_{\mathrm{side}}$, the side wall. The latter is mainly dominated by collisions between mobile ions and the boundary whereas the cap pressure is given as the $(z,z)$ component of the pressure tensor $P_{\mathrm{cap}}=\Pi_{(z,z)}$:

[TABLE]

Here $V=\pi R_{\mathrm{out}}^{2}R_{\mathrm{e}}$ is the effective available volume, $m_{i}$ ($\vec{v}_{i}$) is the mass (velocity) of particle $i$, and $\vec{F}_{ij}$ ($\vec{r}_{ij}$) the force (connection vector) between particles $i$ and $j$. The last term represents the Coulomb contribution to the pressure tensor and is calculated according to Ref. [52]. The side contribution $P_{\mathrm{side}}$ is obtained directly by measuring the average normal force on the constraint and dividing by its area. Having expressions for $P_{\mathrm{cap}}$, $P_{\mathrm{side}}$ and $P_{\mathrm{res}}$, we determine the equilibrium volume using $P_{\mathrm{in}}(V_{\mathrm{eq}})=P_{\mathrm{res}}$. To check the accuracy of the single-chain cell model we compare our equilibrium chain extensions to the ones obtained via the periodic-cell model, cf. Fig. 2 of Ref. [31]. We will discuss the results after describing the second mean-field approximation.

Since the single-chain cell model uses a cylindrical cell a further reduction of the model complexity is to construct an adapted PB description of the polyelectrolyte in the salt solution. As before, the PB model uses a semi-infinite cylinder having an external radius $R_{\mathrm{out}}$ and conceptual length $R_{\mathrm{e}}$, as shown in Fig. 1. The electrostatics for an infinite rod is solved, again mimicking in a simplified way the electrostatic situation in a gel. We also depict how the explicit single chain is now modeled as a penetrable concentric charged rod of radius $a$, characterized by a prescribed charge distribution.

The radius $a$ of the charged rod is chosen such that the polymer density implies the same average distance from the end-to-end vector as obtained from single-chain MD simulations. This amounts to finding the length scale equivalent to the tensional blob size. The numerical value is obtained via fitting a second degree polynomial $a_{\text{MD}}(R_{\mathrm{e}})=N\sigma\left(C_{1}(R_{\mathrm{e}}/(N\sigma))^{2}+C_{2}(R_{\mathrm{e}}/(N\sigma))+C_{3}\right)$, with $C_{1}=-0.17$, $C_{2}=0.14$, $C_{3}=0.03$ based on our single-chain MD data. Enforcing $\langle r\rangle=\int_{V}p(\vec{r})r\,\mathrm{d}V\overset{!}{=}a_{\text{MD}}$ we obtain for the rectangular distribution function $a=3/2a_{\text{MD}}$. The monomer charges thus have a probability density $p(\vec{r})=\mathcal{N}H(-(r-a))$, where $H(x)$ is the Heaviside function and $\mathcal{N}$ a normalization such that $p(\vec{r})$ is a probability density. For a strong polyelectrolyte, the charge is homogeneously distributed in the rod with $\rho_{f}(\vec{r})=-Nfe_{0}p(\vec{r})$.

The Poisson equation describes the electrostatic interaction and ionic distribution in the system:

[TABLE]

The number densities of the ions $c_{i}$ are related to the charge densities via $\rho_{i}(\vec{r})=q_{i}c_{i}(\vec{r})$ given by standard PB theory (with the bulk potential $\psi^{\mathrm{b}}=0$):

[TABLE]

Water is modeled implicitly via a relative dielectric permittivity of $\epsilon_{\mathrm{r}}\approx 80$.

The PB pressure inside the cell $P_{\mathrm{in}}$ has two contributions:

1. The combined ideal and Maxwell pressure [53] which yields

[TABLE]

where $E_{r}=-\partial_{r}\Psi(r)$ is the electric field in radial direction and \langle\mathcal{A}\rangle_{z}=\int_{0}^{2\pi}\int_{0}^{\text{R_{\mathrm{out}}}}r\mathcal{A}(r)dr/(\pi R_{\mathrm{out}}^{2}) denotes the average over all radii. And 2) the stretching pressure $P_{\mathrm{cap}}^{\mathrm{str}}$ (acting only on the cap) which we define to be the pressure due to confinement $P^{\mathrm{conf}}$ (in the spirit of [54, p. 115]) minus the tensile stress $\sigma^{\mathrm{chain}}$ of the chain in order to have a finite extension for a neutral polymer gel:

[TABLE]

where $\mathcal{L}^{-1}$ is the inverse Langevin function. The stretching pressure is constructed such that $P_{\mathrm{cap}}^{\mathrm{str}}(R_{0})=0$, where $R_{0}=1.2bN^{0.588}$ is the average end-to-end distance of an unconfined neutral chain  [31]. For a neutral gel, only the stretching pressure will determine the swelling equilibrium $P_{\mathrm{in}}(R_{\mathrm{eq}})=P_{\mathrm{res}}$ which is found at $R_{\mathrm{eq}}=R_{0}$. The added confinement pressure dominates at low extensions.

Therefore the pressure inside the gel can be obtained via Eq. (1) where the cap pressure is given by $P_{\mathrm{cap}}=P_{\mathrm{cap}}^{\text{ions}}+P_{\mathrm{cap}}^{\mathrm{str}}$. Applying the equilibrium condition $P_{\mathrm{in}}(R_{\mathrm{eq}})=P_{\mathrm{res}}$, we obtain the equilibrium end-to-end distance $R_{\mathrm{eq}}$. All equations were solved with a finite element solver [55].

In the following, we compare the obtained swelling equilibria for both models to the periodic gel model for different charge fractions, chain lengths, and reservoir salt concentrations. For selected parameters, Fig. 2 demonstrates that for all models the gel swells: i) more with increased charge fraction $f$; ii) less with higher salt concentration in the reservoir $c_{s}^{b}$, in good agreement with the data of the periodic gel model. Both models also work for Manning parameters larger than unity, contrary to the Katchalsky model [31]. Further, at high charge fractions ($f>0.5$) they show better agreement with the periodic gel model than the self-consistent field theory presented in Ref. [32]. The PB model has basically the same accuracy as the single-chain MD model for the selected parameter regions. Enlarging the comparison across a wide parameter range to all available data, yielding 60 data points, Fig. 3 shows an excellent agreement of both models against the periodic gel model used here as the reference standard. Our PB model has the known limitations [56, 39]: multivalent ions, high charge densities (e.g. at high compressions of the gel) or high ionic concentrations lead to deviations due to neglecting ionic and excluded volume correlations which also exist in polyelectrolyte gels [57] or charged rod systems [39, 58]. However, these limitations do not apply to our single chain MD model.

We now generalize the PB model to account for weak groups similar to using the charge regulating boundary condition [59]. However, we use charge regulation as a way to determine the space charge density of the penetrable rod. For a weak polyelectrolyte where monomers may be neutral or charged (\ceHA ¡=¿ A- + H+) the titratable monomers (\ceA- or \ceHA) are distributed with $p(\vec{r})$ resulting in a concentration $c_{0}(\vec{r})=Np(\vec{r})$. The dissociation constant is given by $K_{\mathrm{a}}=c(\ce{A-})c(\ce{H+})/c(\ce{HA})=10^{-4}$\mathrm{mol}\text{,}{\mathrm{L}}^{-1}$$. Chemical equilibrium results in:

[TABLE]

and therefore $\rho_{f}(\vec{r})=-e_{0}c(\ce{A-},\vec{r})$. In the case of charge regulation we also explicitly model pH (while neglecting any small \ceOH- concentration). Bulk charge neutrality implies that the sum of the product of all species multiplied with their valency needs to be zero (or equivalently):

[TABLE]

Note that the bulk salt concentration $c_{\mathrm{salt}}^{\mathrm{b}}$ is related to $c^{\mathrm{b}}_{\ce{Na+}}=c_{\mathrm{salt}}^{\mathrm{b}}$ and $c^{\mathrm{b}}_{\ce{Cl-}}=c_{\mathrm{salt}}^{\mathrm{b}}+c^{\mathrm{b}}_{\ce{H+}}$ to ensure charge neutrality.

In Fig. 4, we show the equilibrium extension $R_{\mathrm{eq}}$ as a function of the $\mathrm{p}K_{\mathrm{a}}-\text{pH}$. The gel swells less with lower $\text{pH}=-\log_{10}(c^{\mathrm{b}}(\ce{H+})/$\mathrm{(}\mathrm{mol}\mathrm{/}\mathrm{L}\mathrm{)}$)$ ($\mathrm{p}K_{\mathrm{a}}-\text{pH}$ becoming larger) since the acid becomes less dissociated (less charged). The gel also swells less with higher salt concentration due to increased screening and a higher pressure exerted by the salt reservoir. These findings are in good qualitative agreement with experiments [60, 61].

In summary, we presented two successively coarsening mean-field models that can predict swelling equilibria for charged macrogels. The first one, the single-chain cell model, is a charged bead-spring model with explicit salt ions confined within a cylindrical cell which can undergo affine volume changes. The computational cost for solving the single-chain cell model is at least an order of magnitude lower than for the periodic gel model. In the next model we replaced the charged single-chain and all ions by suitable charge distributions and use the PB framework to derive the equilibrium cylindrical cell length. This model can be solved numerically, i.e., with standard finite element solvers, and is yet at least another order of magnitude faster than the single-chain cell model. Since both models can predict the swelling equilibria in similar good agreement for a wide parameter range with the more elaborate periodic gel model, we can use the extremely efficient PB model for predictions about gel swelling equilibria. The PB model was further generalized to account for charged gels containing weak groups. We find that the gel behavior under different $\mathrm{p}K_{\mathrm{a}}-\text{pH}$ conditions as well as salt reservoir concentrations qualitatively agrees with experimental findings and theoretical expectations.

1 Acknowledgments

The authors acknowledge inspiring discussions with T. Richter and want to thank G. Rempfer and F. Weik for fruitful discussions regarding the pressure. Funding from the DFG through the SFB 716 and Grants HO 1108/26-1 and AR 593/7-1 is gratefully acknowledged.