A self-consistent mean-field model for polyelectrolyte gels

Abstract

We present a novel approach to modeling polyelectrolyte gels, exploiting the analogy between star-branched polymers and polymer networks as a computationally inexpensive yet reliable alternative to full-scale simulations. In the numerical mean-field model of a star-like polymer we modify the boundary conditions to represent an infinite network. We validate the predictions of our new model against a coarse-grained simulation model. We also validate it against a phenomenological analytical model which has been previously shown to agree with simulations in a limited range of parameters. The mean-field model explicitly considers local density gradients and agrees with the simulation results in a broad range of parameters, beyond that of the analytical model. Finally, we use the mean-field model for predictions of the swelling behaviour of weak polyelectrolyte gels under different pH conditions. We demonstrate that the local density gradients are important and that the ionization of the weak polyelectrolyte gel is significantly suppressed. Under the studied conditions the effective pK_A is about one unit higher than that of the free monomer. This shift in the effective pK_A stems from the different pH values inside and outside the gel.

---

1 Introduction

The swelling capacity of polyelectrolyte gels is much greater than that of neutral polymer gels because of the additional osmotic pressure contributions of small ions confined in the gel and mutual repulsion of polymer-bound charges. While modeling of the swelling behaviour of neutral gels can be considered well understood in terms of the Flory–Rehner theory and the related later developments,^1,2 this is not true for polyelectrolyte gels. On the one hand, phenomenological analytical theories exist, which provide relatively simple explanations of the swelling behaviour in terms of different contributions to the free energy.^3–6 These theories typically assume a homogeneous gel phase in equilibrium with the bulk solution that can be characterized by a few macroscopic parameters. On the other hand, particle-based coarse-grained models have been used in several simulation studies.^7–15 The simulations sample particle configurations and therefore explicitly account for correlations between different species. Their disadvantage is that they are much more computationally expensive when compared to analytical theories. When it comes to modeling pH-sensitive (weak) polyelectrolyte gels, where the dissociation of gel segments is a result of chemical equilibrium determined by solution pH, theoretical and simulation literature is quite scarce. On the side of simulations, an additional limiting factor is the necessity to explicitly account for the dissociation equilibrium. On the theoretical side, it is the lack of suitable analytical models that would account for the different ionizations of polyelectrolytes, as compared to free monomers. Therefore the existing literature employs a number of simplifications, e.g., the use of ideal ionization equilibrium to describe the ionization inside the gel⁶ in the analytical models or fixed ionization which is assumed to increase linearly with solution pH.^16,17 An interesting approach to swelling of nanogels has been proposed by Jha et al.,^18,19 which combines explicit particle simulations with a density functional approach used earlier for grafted polyelectrolyte layers.^20,21 Using a similar approach, Longo and coworkers^14,22 combined coarse-grained simulations with local ionization equilibria to study the swelling of pH-sensitive hydrogel layers. Most of the above simulation works have employed a diamond-like topology for their model networks.^7,9–13,15 This is an idealized representation of the most common experimental situation with four chains emanating from each node of the network, e.g. in gels obtained by photocrosslinking. However, model topologies derived from a simple cubic lattice with six chains attached to a node have been used as well.^8,14,22 Several other theoretical works on polyelectrolyte microgels explicitly considered density gradients across the gel–solution interface.^23,24 However, the above approaches, which require explicit particle simulations, are computationally rather expensive. Apparently, an inexpensive method, which would provide a first approximation to the local density variations inside polyelectrolyte gels and to their ionization response, seems to be lacking.

In this work we present a novel approach to model the swelling behaviour of polyelectrolyte hydrogels, based on the numerical self-consistent field model due to Scheutjens and Fleer (SF-SCF),²⁵ which has been used for studies of star-like polyelectrolytes.^26–30 We propose that a covalently cross-linked polyelectrolyte hydrogel can be viewed as a collection of spherically symmetric star-like polyelectrolytes, mutually connected together by end segments of the arms, as schematically shown in Fig. 1.

Fig. 1 (A) Network of polyelectrolyte chains is considered as the set of stars. The arm ends, representing midpoints of the gel strands, are fixed at a certain distance, R_star, from the centre. (B) Representation of the star as random walks on the one-dimensional array of concentric spherical layers with the end segments pinned in a particular layer given by R_star. The box size R_box is chosen such that the volume of the coloured regions in (A) and (B) is equal. (C) Simulation snapshot showing the simulation model of the gel. Periodic images of the gel outside the unit cell are shown as well to illustrate the quasi-infinite connectivity of the network.

We demonstrate the capabilities of the numerical SF-SCF approach by comparing results to recent simulation data and to an analytical model that has been shown to successfully predict the swelling of polyelectrolyte gels within a limited range of parameters.¹³

The newly proposed approach has the advantage that it explicitly treats the concentration gradients of all the involved species in the radial direction with respect to the node of the hydrogel network (the centre of the star in the mean-field representation). Once validated against the simulation data on strong polyelectrolyte gels, the SF-SCF approach can be easily extended to predict the swelling of weak polyelectrolyte gels with pH-dependent ionization and to poor solvent conditions, where radial variations of the properties significantly affect the behaviour that result in non-trivial features. The spherical symmetry approximation used in the model provides a correction to the most relevant density gradients in the radial direction, but neglects the angular correlations. At the same time, it allows us to reduce the three-dimensional problem to a one-dimensional one, which saves a tremendous amount of computational cost. In addition, our approach provides a guideline to extrapolate the existing theoretical knowledge about star polyelectrolytes to polyelectrolyte gels.

2 Models

The target system represented by our models is a hydrogel network made up by polyelectrolyte chains. All chains consist of the same number of segments N = 79, and are connected to a diamond-like network with f = 4 chains connected to each node. We consider two situations: (i) a strong polyelectrolyte in which some segments are permanently ionized while others are not; (ii) a weak polyelectrolyte in which all segments are subject to the ionization reaction

HA ⇋ A⁻ + H⁺ (1)

where HA and A⁻ denote the protonated and deprotonated form of the acidic monomer unit A. We denote by α the fraction of ionized groups on the chains, such that on average there are αN charged segments per chain. The gel is in equilibrium with a bulk reservoir of monovalent salt solution at concentration c^b_s. The salt concentration inside the gel, c^g_s, is the result of our calculation, determined by the equivalence of chemical potentials of mobile ions in the bulk and in the gel phase. For simplicity of notation, we refer to the salt as ‘NaCl’.

In the strong polyelectrolyte case we fix the degree of ionization α = 1, 1/2, 1/4, and 1/8 with regularly spaced charges, such that one ionized segment on the polymer is followed by (1/α − 1) non-ionized ones. In the weak polyelectrolyte case, the reaction in eqn (1) is characterized by an equilibrium constant, K_A. Then pK_A and pH are the input parameters of the model, and α is the result of the calculation. We choose pK_A = −log₁₀ K_A = 5 as a typical value for monomers containing weak carboxylic acids, such as acrylic acid.

In the weak polyelectrolyte case, we use α_s to denote the probability that the segment with ranking number s (counted from the network node) is ionized. Then we define the overall degree of ionization of the gel, α as the average over individual segments:

Similarly we define the local degree of ionization, α(r), as a fraction for ionized segments out of all segments at given distance r from the network node.

For the bulk salt concentration we use the values 0.01 M ≤ c^b_s ≤ 0.2 M in order to cover the range from potable water to seawater, motivated by the potential application of charged hydrogels to water desalination.^12,31 The analytical calculations and the numerical mean-field calculations were performed also for a broader range of values of α and c^b_s than those mentioned above.

We use three different models to study the properties of the gels, each of them employing a different level of approximation.

• The analytical model based on our modification of the theory of Katchalsky and Michaeli,^3,13 hereafter referred to as the analytical model.

• The coarse-grained simulation model of the Kremer–Grest type with explicit particles and charges, and implicit solvent, identical to our preceding publication,¹³ hereafter referred to as the simulation model.

• The mean-field model based on the Scheutjens–Fleer self-consistent field (SF-SCF) model for star-like polyelectrolytes, and on the analogy between gels and stars, hereafter referred to as the mean-field model.

Although the above-mentioned models all represent the same hydrogel, due to the different approximations they differ in some specific details of the representation. In the following paragraphs we briefly outline each of these representations, focusing on the differences. For the analytical model and for the simulation model we only provide a brief description, and refer to the preceding publication for additional details.¹³

2.1 The analytical model

The analytical model of Katchalsky and Michaeli,^3,13,32 that we extended recently to account for finite extensibility of the chains, is based on the analogy between the swelling of a polyelectrolyte gel and stretching of an individual strand between two crosslinks. The end-to-end distance of the strand, R_e, and the volume per chain of the gel, V, are coupled assuming affine deformation:

R_e³ = AV, (3)

where the coefficient A is defined by the topology of the network. For the diamond network used here (see Fig. 1), . The affine deformation holds well close to the full extension but deteriorates for less extended chains, which results in systematic errors of up to 30% in the predicted end-to-end distance.¹³ The mechanical equilibrium is found by balancing different contributions to the pressure (each of them being the derivative of the corresponding free energy contribution):The free swelling equilibrium is attained at the external pressure P_ext = 0, while P_ext > 0 corresponds to a gel under compression, and P_ext < 0 corresponds to a gel under tension. The short-range interaction contribution is assumed to be negligible under good solvent conditions, P_int ≈ 0. The electrostatic contribution, P_el, the osmotic contribution, P_osm, and the chain stretching contribution, P_str, attain the form:

P_osm = k_BT(2c^g_s + αc_p − 2c^b_s), (7)

where a is the segment size, c_p is the polymer concentration in the gel, λ_B is the Bjerrum length and [script L]⁻¹(R_e/aN) is the inverse Langevin function to account for finite extensibility. The term P_el in eqn (4) accounts for the intra-chain electrostatic repulsions. Parameter ξ in eqn (5) has been introduced to simplify the notation:where κ is the inverse Debye screening length inside the gel, V₀ = R₀³/A is the volume of the corresponding neutral gel, and R₀ = 1.18N^0.588 is the end-to-end distance of the corresponding neutral chain in athermal solvent. The prefactor 1.18 has been determined numerically for our simulation model.¹³

With V₀ as the reference state, we define the swelling ratio as

Q = V/V₀. (9)

The salt concentration in the gel, c^g_s, which is needed as an input for the calculation of P_osm in eqn (7), and of the Debye screening length in eqn (8), couples the swelling equilibrium and salt partitioning:

With C = 1 eqn (10) reduces to the Donnan equation for partitioning of charged solute across a semi-permeable membrane. The Donnan equilibrium explicitly considers electroneutrality but otherwise assumes ideal behaviour of the solute. Deviations from the ideal behaviour due to interactions between the polymer and small ions are comprised in the term C. In the numerical implementation, we calculate the external pressure viaeqn (4) for a series of values of R_e to obtain the pressure–extension relation. Interpolating between these values, we obtain the gel swelling ratio at the desired pressure, e.g. P_ext = 0 for free swelling equilibrium.

2.2 The coarse-grained (CG) simulation model

Here, we restrict ourselves to a brief summary, since this model is identical to the one used in our previous publications,^12,13 where a full-detailed description can be found. Similar models have been used by other authors to study polyelectrolyte gels.^{9–11,33–35} In this work, we extended our earlier simulation results beyond the Manning condensation threshold³⁶ to α = 1.0. In brief, the gel network consists of polymer chains interconnected in a diamond topology via the finite-extensible nonlinear elastic potential, with excluded volume modeled by a purely repulsive Weeks–Chandler–Andersen potential. Charged segments, their counterions, as well as salt ions are modeled explicitly as charged particles with excluded volume, while solvent is treated as a dielectric continuum with the relative permittivity ε_r = 80. Equilibrium with the salt reservoir is ensured by insertion/deletion of salt ions in a semi-grand canonical Monte Carlo manner. The effective particle size and bond length are set to a ≈ 0.35 nm, which ensures that the linear charge density at full ionization is λ_B/a ≈ 2, where λ_B ≈ 0.71 nm is the Bjerrum length in water at 300 K. The same ratio of λ_B/a yields the concentration dependence of the chemical potential of the coarse-grained salt solution which is consistent with experimental data on NaCl up to c^b_s ≈ 1 M.¹³

We carry out a series of molecular dynamics simulations with the Langevin thermostat for fixed sizes of the simulation box to obtain the pressure–extension relation and interpolate to find the gel volume at a given pressure.

2.3 The numerical mean-field model

Here we exploit the idea originally suggested to us by A. Polotsky, that a gel can be viewed as a collection of interconnected star-like polymers under tension, as illustrated in Fig. 1. As discussed below, this provides some advantages over the traditional depiction as a collection of linear chains.

Our mean-field model of a crosslinked gel is derived from the model of a star-like polyelectrolyte,²⁷ based on the Scheutjens–Fleer self-consistent field (SF-SCF) approach described in more detail in Section 2.4. The star is placed in the centre of the spherical coordinate system, and spherical symmetry is assumed, explicitly accounting for the radial gradients in the densities and neglecting all angular correlations.

To map the star to the gel, we use f = 4 arms, each of which has half the length of the strand connecting two nodes of the network, N/2. The key aspect of representing the gel as a star consists of fixing the end segment of each arm at a particular distance,

R_star = R_e/2, (12)

where R_e is the extension of the corresponding gel strand. We then fix the radius of the simulation box to

R_box = (Aπ/12)^−1/3R_star, (13)

in order to preserve the volume per chain given by eqn (3), as illustrated in Fig. 1.

Besides the star we introduce salt ions, Na⁺, Cl⁻, and water. The bulk concentrations of the salt ions and the H⁺ ions are fixed in the outermost layer of the box, which defines the boundary condition for the semi-grand canonical calculation. The concentration of OH⁻ ions is determined by explicitly considering the autoprotolysis reaction of water. The salt concentration in the gel is then obtained as the total number of salt ions in the box, divided by the box volume:

where r is the distance from the centre of the box. In the same manner we define the concentration of H⁺ ions in the gel, .

In order to get the pressure from the mean field model, we use the partial open Helmholtz free energy, defined as . This quantity is analogous to the Grand potential for open systems. The main difference is that our system is open to all particles but the macromolecule, which is fixed at the centre, hence the term partial open. We calculate pressure as a numerical derivative using the seven point stencil³⁷

Finally, we determine the swelling ratio from the mean-field model by linear interpolation between the values of R_star which are closest to P_ext = 0 usingwhere R^ref_star corresponds to the free swelling equilibrium of the same star but with α = 0.

2.4 The Scheutjens–Fleer self-consistent field implementation

In the Scheutjens–Fleer self-consistent field method, SF-SCF,^25,38 the free energy is expressed as a function of the density profiles of all the system components, i.e. polymer segments, free ions, and the solvent. This method has been described many times in the literature. Therefore we only point out the essential parts of it and refer the reader to the existing literature for further details.^25,27,29 The one-dimensional spherical approximation employed here reduces the three-dimensional problem to a one-dimensional spherically symmetric problem, which greatly reduces the number of iteration variables and speeds up the search for the optimum solution. Assuming spherical symmetry, the central segment of the star is placed at the centre of a system of M concentric shells (layers) of constant thickness a = 0.35 nm, which also defines the effective size of the monomer and of all other species in the system. The mean-field potential u_x(z) experienced by segments of type x in layer z can be expressed aswhere e is the elementary charge, ν_x the valency of segment x. The term u_x′(z) is the Lagrange field due to the incompressibility constraint,²⁵ and the second term accounts for the electrostatic interactions. The last term accounts for the short-ranged interactions by means of the Flory–Huggins parameter χ. The notation 〈φ(z)〉 stands for the density averaged over layers adjacent to z. The suitable value of χ is determined by quantitative comparison with the simulations as explained in Section 3.1.1. The local electrostatic potential, ψ(z), is obtained by solving the discretised Poisson equation in the spherically symmetric system, with the density profiles of all charged species as input parameters. As is implied by eqn (17), angular density gradients within a given layer are neglected, while density gradients in the radial direction are explicitly considered in the potentials. All components in the system follow the Boltzmann distribution with respect to the potentials, and therefore the statistical weight of finding component x in layer z is given by

G_x(z) = exp(−u_x(z)) (18)

The SF-SCF formulation introduces a propagator scheme, which allows us to compute the density profiles from the given potentials, taking into account the connectivity of chain-like molecules and all their possible conformations which satisfy the given constraints – in our case it is fixation of the first segment at z* = 1. Then the propagator scheme is as follows:

G_x(z,s|z*,1) = G_x(z)〈G_x(z,s − 1|z*,1)〉 (19)

G_x(z,s|N) = G_x(z)〈G_x(z,s + 1|N)〉 (20)

where G_x(z,s|z*,1) denotes the statistical weight of a segment with ranking number s in layer z for a conformation with segment 1 in layer z*, and G_x(z,s|N) is the statistical weight of segment s in layer z for a conformation which starts with segment N in any layer. With the initial conditions G_x(z,N|N) = G_x(z) and G_x(z*,1|z*,1) = G_x(z*) eqn (19) and (20) can be solved recursively. The volume fraction profiles, φ_x(z,s), then follow aswhere C_x is the normalization constant. The volume fractions of each species are related to molar concentrations as c = φ/(1000N_Aa³) mol dm⁻³.

For weak polyelectrolytes it is necessary to consider different states of the ionizable species (polymer segments and water molecules). Then the above equations are modified such that the fraction of species of type j in state k is given by³⁹

α_j,k(z) = α^b_j,k(z)G_j,k(z)/G_j(z) (22)

where α^b_j,k(z) is the corresponding fraction of species j in state k in the bulk, i.e., in the absence of any field. Because the ionized and non-ionized states have different charges, they experience different potentials and attain different statistical weights in each layer. At the same time, the local concentrations of ionized and non-ionized polymer segments are coupled by the relationwhere α(z) is the local degree of ionization in layer z. Similarly, the autoprotolysis of water, determined by its ionic product, K_w = 10⁻¹⁴ mol² dm⁻⁶, is considered locally within each layer:

c_H⁺(z)c_OH⁻(z) = K_w (24)

In this way, the local variations of pH and polymer ionization are explicitly included in the mean-field model.^39,40

The solution of the set of equations starts from an initial guess of the density profiles, which yield the respective potentials and the statistical weights. Application of the propagator scheme provides a new set of density profiles. The self-consistent solution is reached when the old and the new volume fractions match within a pre-defined accuracy threshold. The final set of density profiles and potentials yield an accurate prediction of the free energy within the mean-field approximation, which consists of replacing explicit inter-particle interactions with a mean-field potential as given in eqn (17). Due to the variational principle, the mean-field free energy is an upper bound to the true free energy with all correlations included. This difference is inherent to the involved approximation and is not related to the above mentioned accuracy threshold.

The SF-SCF formulation is a discrete analogue of the continuous Edwards diffusion equation.^41,42 The discrete version implemented on a set of concentric shells (“the lattice”) with lattice spacing a allows for the full enumeration of all the chain configurations starting in the first layer and ending in the layer number given by R_star, see figure Fig. 1B. This automatically accounts also for finite chain extensibility. The discretised lattice has an important consequence for our depiction of the gel as a collection of stars, sketched in Fig. 1. The SF-SCF approach requires that both R_box and R_star are integer multiples of lattice spacing a. If we choose R_star to be an integer multiple of a, then R_box obtained from eqn (13) usually is not an integer multiple of a. Therefore we always perform two calculations, with R_box rounded to the nearest lower and the nearest higher integer multiple of a. These provide us with an upper and lower estimate of the desired result. Therefore, although the SF-SCF calculation can be performed to an arbitrary accuracy, the discretization limits the accuracy of the representation of the gel. All results that we report are arithmetic means of the two calculations, while the error bars represent the upper and lower bound to indicate discretization artifacts.

3 Results and discussion

3.1 Model validation: strong polyelectrolyte gels

3.1.1 The pressure–extension relation. As a starting point, we calibrate the excluded volume interactions in the mean-field model such that its predictions of the pressure as a function of chain extension (gel swelling) for the neutral gel quantitatively agree with the simulation model. As deduced from Fig. S7 in the ESI,† the quantitative agreement is obtained with χ = −1.5. With this calibration our model is ready to be used for the prediction of the swelling of charged gels. Fig. 2 shows the gel pressure as a function of chain extension, comparing the three different models for a selected set of system parameters (see the ESI† Fig. S4 for analogous results for other parameter sets). At α = 1/4, both the analytical model (dashed line) and the mean-field model (empty points connected with lines) reasonably well reproduce the simulation results (full points). This description holds also for lower α. The agreement between the simulations and both the analytical and the mean-field model improves at lower c^b_s. An increase to α = 1/2, which is the Manning condensation threshold (λ_B = a/α), results in poorer agreement. While the analytical model fails to reproduce even the qualitative trend from the simulations, the mean-field model succeeds qualitatively, but at the quantitative level overestimates the pressure and predicts higher chain extension at free swelling equilibrium (P_ext = 0). The analytical model exhibits an undulation close to P_ext = 0, which is not indicated by the simulation data. Beyond the Manning condensation threshold, at α = 1, the undulation predicted by the analytical model grows out of the plotting bounds. This is mainly due to the electrostatic pressure term, which loses its validity in this range of parameters. If this term is omitted, the predicted pressure–extension curves retain the qualitatively correct shape, but lose quantitative agreement in the whole range of parameters. In contrast, the mean-field prediction is qualitatively correct, but further deteriorates in quantitative terms. Interestingly, at higher α, we observe pronounced differences in performance of the models with varying c^b_s (see the ESI† Fig. S4 for the corresponding plots). At α = 1/2 the analytical model works well for short chains, and quantitatively outperforms the mean-field model at N = 39. It exhibits initial signs of failure at higher N and high c^b_s, although at high c^b_s it still works comparably well with the mean-field model. At α = 1 the analytical model fails qualitatively in the whole range of studied parameters, while the mean-field model exhibits the same behaviour as for low α: overall qualitative agreement and improving quantitative agreement with increasing N or c^b_s. The improvement of the mean-field approach with increasing N can be understood in terms of discretization effects, which are more pronounced at small N. The addition of salt screens the direct interactions between neighbouring groups on the chains, which we have identified in previous work as the key component of the discrepancy between the mean-field and explicit particle treatment.²⁹ Therefore it is not surprising that the mean-field model improves with increasing c^b_s. In our previous work we have shown that for star polyelectrolytes the agreement between mean-field predictions and simulations improves with an increasing number of arms.²⁹ In this context f = 4 arms in a star is the lower limit for the mean-field prediction to be actually applicable. A higher value of f would presumably lead to better agreement between the mean-field and the simulation model; however, such connectivities are not commonly encountered in polymer gels.

Fig. 2 Pressure–extension curves for a gel with strand length N = 79, salt concentration c^b_s = 0.1 M, and three different degrees of ionization: α = 1/4 (green), 1/2 (blue) and 1 (red). Empty points connected with lines represent the mean-field results, full points of matching shape and colour are the simulation results, and the dashed lines of matching colour are the corresponding predictions of the analytical model.

Fig. 3 shows the salt partitioning, expressed as the ratio of salt concentration in the gel to that in the bulk, as a function of the total concentration of charged groups bound to the gel, αc_p. Following the simple Donnan prediction, the ratio c^g_s/c^b_s should be a universal function of αc_p, independent of other parameters. The Donnan prediction does not explicitly account for any inter-particle interactions (except for electroneutrality), and in this respect plays the role of the ideal gas approximation. Deviations of the simulation data (empty points in Fig. 3) from the Donnan theory provide a measure of non-ideality, which is well captured by the analytical model at low c^b_s and low α (see also ESI† Fig. S5 for additional data). This agreement deteriorates with increasing c^b_s and α. In particular at α = 1 and high c^b_s the analytical model fails even at the qualitative level, predicting that initially c^g_s/c^b_s increases (see also ESI† Fig. 5), which is not observed in the simulation data. This failure is presumably related to the convex shape of the predicted pressure–extension curves in Fig. 2 in the same range of parameters. On the other hand, the mean-field model again retains its qualitative agreement throughout the whole range of parameters. The rather big estimated error in the mean-field model, which increases with α, stems from discretization, as discussed in Section 2.4. Within the accuracy limits set by the discretization, the mean-field model agrees with the simulation results and with the analytical model up to α = 1/2 and higher salt concentrations. At the lowest salt concentrations and α = 1/2, the mean-field model systematically underestimates c^g_s, and is outperformed by the analytical model. At α = 1, the mean-field model provides poor quantitative predictions, but qualitatively the predicted trend does not contradict the simulation data as it occurred in the case of the analytical model. Last but not least, we should mention that at the highest c^b_s the uncertainty due to discretization becomes bigger than the difference between simulation results and the Donnan prediction, which makes the mean-field results useless.

Fig. 3 The gel to bulk ratio of salt concentration, c^g_s/c^b_s, as a function of the concentration of bound charges in the gel for a system with strand length N = 79, salt concentration c^b_s = 0.1 M, and three different degrees of ionization: α = 1/4 (green), 1/2 (blue) and 1 (red). The black solid curve displays the Donnan prediction as a reference. The empty symbols show the results of the mean-field calculations, full symbols of matching shape and colour are the corresponding simulation results, and the dashed lines of matching colour represent the corresponding analytical predictions. Error bars of the mean-field results have been omitted for clarity. See the ESI† Fig. S5 for the same plots with error bars.

At the qualitative level, the mean-field model captures all the trends that emerge from the simulation data, in particular the variation of salt content and the pressure–extension relation with the crosslink length, N, the degree of ionization of the polymer, α, and the salt content of the bulk solution, c^b_s.

At the quantitative level, the mean-field model underestimates the salt content in the gel, and consistent with that it underestimates the pressure inside the gel, and hence overestimates the swelling ratio. This discrepancy produces a systematic deviation which is consistent throughout the whole parameter range. The pressure–extension relation predicted by the SF-SCF approach remains consistent with the simulations (apart from the above-mentioned systematic deviation) even at high c^b_s or high α, where the analytical model fails. At high α the analytical model predicts a non-monotonic pressure–extension relation, while both the mean-field model and the simulations predict a monotonic dependence.

3.1.2 Swelling of the strong gels. In Fig. 4 we show the swelling ratio of the gel as a function of the degree of ionization and for various salinities of the bulk solution. At low α the analytical model very well matches the simulation results, while the SF-SCF model tends to underestimate the swelling. With increasing α both models follow very well the trend of the simulation results. Around and beyond the Manning condensation threshold, α = 1/2, the analytical model fails to qualitatively reproduce the pressure–extension curves from the simulations, and the predictions of Q(α) from the analytical model cannot be considered reliable. Therefore we plot the corresponding curves in Fig. 4 only in the range where the qualitative agreement is still retained. This is expected since the Katchalsky model does not account at all for the counterion condensation. In contrast, the mean-field model overshoots the predicted value of Q from the simulations but does not lose qualitative agreement. The spherically averaged mean-field model cannot fully account for the counterion condensation either, but it explicitly accounts for the radial correlations between the local densities of charged polymer segments and the counterions. If we consider predictions of Q by different models as a function of salt concentration at fixed α (see ESI,† Fig. S2), it turns out that both the analytical model and the numerical mean-field model reproduce comparably well the gel collapse with increasing c^b_s. The agreement is not quantitative but, as discussed in our previous publication,¹³ it can be partly attributed to the failure of the assumption of affine deformation.

Fig. 4 The swelling ratio of the strong gel as a function of α at various salt concentrations. In this figure, the colour code distinguishes different salt concentrations in the bulk. Full points display simulation data, lines show the numerical mean-field results, and dashed lines represent the analytical predictions.

3.1.3 Density profiles. The advantage of the mean-field model is that it explicitly considers inhomogeneous distributions of all species in the radial direction with respect to the nodes of the gel network (centre of the star). Two examples of such profiles are shown in Fig. 5 for the case of moderate salinity, and two different degrees of ionization. Additional profiles for other studied systems can be found in the ESI,† Fig. S3. The density profiles of monomers and counterions at a high degree of ionization, α = 1, both follow the power law decay ρ(r) ∼ r⁻², theoretically predicted for strongly charged star-like polyelectrolytes.^26,43 On the other hand, the profiles of weakly ionized gels, α = 1/4, follow the same power law at low salinity, while at high salinity they are close to the ρ(r) ∼ r^−4/3 decay of neutral star-like polymers.^43–45 Comparing the profiles of counterions and coions, we observe, that up to c^b_s ≈ 0.1 M, the charge density in the whole gel is dominated by the counterions of the polymer, while the salt concentration (given by the concentration of coions at a given point) is much lower. Only at c^b_s ≳ 0.1 M the situation reverses, and the charge density is dominated by the salt, apart from the region close to the node. Comparing the profiles from the simulations (data points in Fig. 5) with the mean-field predictions (corresponding lines in Fig. 5) we notice nearly quantitative agreement in the case of polymer segments and counterions. This agreement deteriorates at very short distances from the node, where the excluded volume contribution is overrated by the mean-field. A more significant disagreement shows up at r > R_star, where the mean-field density of monomers sharply drops. Consequently, the mean-field counterion density near R_star drops more dramatically than in the simulation data. The range and magnitude of this drop might seem to be insignificant on the logarithmic scale. However, when we realize that it represents a major part of the total gel volume, we may conclude that the systematic underestimation of the salinity inside the gel by the mean-field approach stems from the same roots. Apparently, the drop in the counterion concentration in the distant part of the mean-field profile is due to the missing charged polymer segments beyond R_star, while in a real gel the segment density beyond R_star is nearly constant, passing through a shallow minimum and then increasing again at separations where the nearest-neighbour node can be found.

Fig. 5 Density profiles from simulations (points) and from the SF-SCF method (lines) at intermediate chain length and salt concentration, at the two extreme ionization degrees: α = 1/4 (top panel), and α = 1 (lower panel). Different particle types (monomers, counterions, coions) are distinguished by colour. The two gray lines with slopes −4/3 and −2 show the scaling prediction for the neutral star and a highly charged star, respectively. See the ESI† Fig. S3 for plots with other combinations of parameters.

3.2 pH-Sensitive gels

3.2.1 Pressure–extension curves and salinity. Having established the range of reliability of our mean-field model in the previous section, we now demonstrate its capabilities for the predictions of the swelling of weak polyelectrolyte gels with variable ionizations which depends on the local properties of the system. For this purpose we choose a gel with all segments being weak acids with pK_A = 5. With pH ≳ pK_A, the ionization of such a weak polyelectrolyte gel should vary with the monomer position in a similar manner as it does in a star-like polyelectrolyte in solution.^25,27–29 In addition, its degree of ionization should vary with compression, as is shown in Fig. 6. At pH = 7 = pK_A + 2 we observe α ≈ 0.9 at full chain extension, which slowly decreases down to 0.6 upon compression. Note that at this pH, a low-molecular acid would be fully ionized at almost any concentration. At lower pH values we observe the same trend: α at full extension decreases with pH, and for a given pH the ionization decreases further with compression until eventually, at pH = 4, we get α ≈ 0 throughout the whole range.

Fig. 6 Degree of ionization of the weak polyelectrolyte gels as a function of chain extension at a fixed salt concentration of c^b_s = 0.01 M and for different pH values.

This variation of α with chain extension has a pronounced effect on the pressure–extension curves, shown in Fig. 7. With decreasing pH the curves become less steep (the gels are softer) at a given value of R_e, and consequently the free swelling equilibrium is attained at a lower chain extension and a lower degree of ionization. Interestingly, if we compare the pressure–extension curves at a given pressure (free swelling equilibrium, P = 0) rather than at a given compression, it turns out that the more ionized gels are softer than the less ionized ones. The trend in the pressure–extension curves with decreasing pH is qualitatively similar to the trend with decreasing ionization in the strong polyelectrolyte gels (Fig. 2). The difference becomes apparent by comparing the data for weak gel at pH 5.5 and a strong gel at α = 0.35. This is the average value of α for the weak polyelectrolyte gel at free swelling equilibrium (zero applied pressure). For the strong polyelectrolyte this was realized by making all segments ionized with a valency of 0.35. In the weak case α decreases with compression, which makes the pressure–extension curve less steep compared to that of the corresponding strong gel.

Fig. 7 The pressure–extension curves for the weak polyelectrolyte gels at a fixed salt concentration and different pH values. The data for the strong gel correspond to the same α as the weak gel at pH 5.5 and P_ext = 0.

Predictions of the SF-SCF model for the salinity inside the weak gels exhibit only weak deviations from the Donnan prediction, very similar to the strong case (see the ESI,† Fig. S1). Since the Donnan prediction is a universal function of αc_p, the fact that in the weak gel α varies with compression does not significantly change the situation. However, as deduced from the discussion of the strong case, the actual deviations from the Donnan model should be stronger than those predicted by the SF-SCF model.

3.2.2 Density profiles and local properties of weak gels. Unlike in a strong polyelectrolyte gel, in a weak one the degree of ionization is a locally varying property. Fig. 8 reveals that the local degree of ionization, α(r), is very low near the node (low r), and monotonically increases with r. This leads to a different distribution of charges in the weak case, and therefore density profiles of the ionized polymer segments at a given value of R_star are different from those of the strong gel at the same α and R_star (see also ESI,† Fig. S6 for density profiles of all species in the weak gels). In particular, the density profile of the H⁺ ions, shown in Fig. 8, follows the same (Boltzmann) distribution as the salt cations, Na⁺: c(r) ∼ exp(−ψ(r)/k_BT), where ψ(r) is the local electrostatic potential (see the plot in the ESI† Fig. S6). However, they differ in the pre-factor, which is proportional to the bulk concentration of H⁺ and Na⁺, respectively. From Fig. 8 one can see that up to r = R_star the profile of H⁺ ions approximately follows the profile of the ionized monomers. Beyond R_star the value of c_H⁺(r) in our model slowly decays to its bulk value at R_box. This is the same situation as discussed in the strong case (Fig. 5): beyond R_star our model predicts a decrease of the counterion concentration down to its bulk value, while the simulation results show that it should be roughly constant in this range. Irrespective of the above-mentioned artifact, Fig. 8 reveals that both c_H⁺(r) and α(r) inside the gel significantly differ from their bulk values. Because the magnitude of the difference is probably underestimated, all pertinent deviations predicted by our model should be considered as a lower bound to the actual ones.

Fig. 8 Density profiles of polymer segments (top), H⁺ ions (middle) and the degree of ionization of the polymer (bottom) for weak polyelectrolyte gels at free swelling equilibrium, fixed salinity and various pH values.

3.2.3 Titration curves and local pH. Fig. 9 shows that the increase in the swelling ratio of the weak gel with increasing pH follows the trend of increasing ionization of the gel. In the first approximation, the degree of ionization of a weak acid follows the ideal titration curve

1/α = 1 + 10^pK_A–pH. (25)

Polyelectrolyte chains generally exhibit deviations from this ideal behaviour, and their ionization is suppressed.²⁹ A comparison with the ideal titration curve of the free monomer in Fig. 9 reveals that the ionization of the gel is shifted about one unit towards higher pH values. At the mean-field level, this shift can be explained by the different c_H⁺ concentrations in the gel and in the bulk. For this purpose, let us define the “Local pH” inside the gel aswhere is the average concentration of H⁺ inside the gel, defined in analogy with eqn (14). As follows from eqn (23), the ideal titration curve is satisfied locally within each individual layer, if the local pH is used instead of the bulk pH. However, because α and the local pH are coupled by a non-linear relation, if we plot the average α over the whole box as a function of the average pH inside the gel in Fig. 9, this curve does not fully coincide with the ideal titration curve.

Fig. 9 Left axis: Degree of ionization of the gel, α, as a function of pH in the bulk (squares), and the same quantity as a function of the local pH in the gel (circles). The red solid curve displays the ideal titration curve for comparison. Right axis: Equilibrium swelling ratio, Q_eq, (stars), as a function of pH. Inset: Local pH in the gel as a function of pH in the bulk.

As discussed in detail in our previous publication,²⁹ the shift and deformation of titration curves of star-branched weak polyelectrolytes have two important contributions: (1) inter-arm correlations, which lead to higher c_H⁺(r) and lower α(r) inside the gel as compared to the bulk; (2) intra-arm correlations due to neighbouring groups along the chain. The one-dimensional model of a polyelectrolyte star well accounts for the first but underestimates the second contribution. The latter cannot be explained in terms of the (average) local pH in the gel and becomes important especially at α > 1/2 when the distance between the nearest ionized groups along the chain becomes smaller than the Bjerrum length. Therefore, our mean-field model should reasonably well describe weak polyelectrolyte gels at low α but it underestimates the shift of the titration curve at high α.

4 Conclusions

In this work, we introduced a numerical mean-field model for the swelling of polyelectrolyte gels, based on the analogy between gels and star-like polymers. This model employs the Scheutjens–Fleer self-consistent field approach and assumes spherically symmetric density profiles around the node of the gel network, in analogy with the central segment of star-like polymers. We validated this new model, comparing its predictions with our earlier simulation results and with the analytical model for strong polyelectrolyte gels. Below the Manning condensation threshold, the mean-field model provides quantitatively more accurate predictions of gel swelling than the analytical model. Beyond the condensation threshold the analytical model fails even at the qualitative level, while the mean-field model retains a qualitative agreement with the simulations. For the prediction of salt concentration in the gel, the simulation data exhibit positive deviations from the ideal-gas Donnan approximation. While both models capture the general trend of the deviations, the analytical model tends to overshoot them and the mean-field model tends to underestimate them.

Next, we demonstrate a straightforward extension of the new mean-field model to weak polyelectrolyte gels and provide predictions of their swelling as a function of pH. In the weak case, the local coupling between the conformation and ionization is included in the mean-field model and proves to be essential to capture important effects. The degree of ionization of the weak gel changes with compression, and therefore the pressure–extension curves at a given pH are less steep than those for the strong gels at a fixed degree of ionization. The Donnan-like partitioning of the mobile ions between the gel and the bulk also applies to H⁺ and OH⁻ ions. This results in a lower pH inside the gel than in the bulk, and in the suppressed ionization of the gel as compared to the free monomer. This can be interpreted in terms of effective pK_A which is about one unit higher than that of the single monomer. Within the mean-field approximation this suppression of ionization is dominated by the difference between the local and bulk pH, but there is also a clear contribution from the inhomogeneous distribution of charges. The latter is presumably underestimated by the mean-field and should be even more significant in a real system. Some of the above effects were neglected in earlier analytical theories which assumed a homogeneous density inside the gels.

In terms of computational effort, the mean-field calculations are more expensive than solving the analytical model. However, in both cases one calculation of the pressure–extension curve completes within few seconds, so that it can be performed interactively. In contrast, the explicit simulations, which provided our reference data, require several days per each data point on the pressure–extension curve. From this point of view, the new mean-field model provides a promising and cheap alternative for the reliable prediction of the swelling of polyelectrolyte gels in cases where local density gradients are important, simple analytical approaches fail, and explicit simulations are too demanding even at the coarse-grained level.

Acknowledgements

We gratefully acknowledge Dr Lucie Nová for reading the manuscript and providing valuable comments. P. K. and O. R. acknowledge support from the Czech Science Foundation grant (P208/14-23288J) and the Ministry of Education, Youth and Sports (CUCAM CZ.02.1.01/0.0/0.0/15\_003/0000417). P. K., O. R. and O. B. acknowledge support from the Marie Curie International Research Staff Exchange Scheme Fellowship (PIRSES-GA-2013-612562) within the 7th European Community Framework Programme. C. H. and T. R. were supported by the DFG grant HO 1108/26-1 and AR 593/7-1. Computational resources were provided by the CESNET LM2015042 and the CERIT Scientific Cloud LM2015085, provided under the programme “Projects of Large Research, Development, and Innovations Infrastructures”. This work has been partially supported by the Government of Russian Federation, Grant 074-U01, and by the Russian Foundation of Basic Research, Grant 14-03-00372.

References

1. P. J. Flory and J. Rehner, J. Chem. Phys., 1943, 11, 512.
2. M. Doi, Introduction to Polymer Physics, Clarendon Press, Oxford, 1996.
3. A. Katchalsky and I. Michaeli, J. Polym. Sci., 1955, 15, 69.
4. Responsive Gels: Volume Transitions I, ed. K. Dušek, Springer-Verlag, New York, 1993, vol. 109.
5. Y. Huang, I. Szleifer and N. A. Peppas, Macromolecules, 2002, 35, 1373–1380.
6. A. A. Polotsky, F. A. Plamper and O. V. Borisov, Macromolecules, 2013, 46, 8702–8709.
7. S. Schneider and P. Linse, Eur. Phys. J. E: Soft Matter Biol. Phys., 2002, 8, 457–460.
8. Z.-Y. Lu and R. Hentschke, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 061807.
9. S. Edgecombe, S. Schneider and P. Linse, Macromolecules, 2004, 37, 10089–10100.
10. B. A. Mann, R. Everaers, C. Holm and K. Kremer, Europhys. Lett., 2004, 67, 786–792.
11. B. A. Mann, C. Holm and K. Kremer, J. Chem. Phys., 2005, 122, 154903.
12. J. Höpfner, T. Richter, P. Košovan, C. Holm and M. Wilhelm, Intelligent Hydrogels, Springer International Publishing, 2013, vol. 140, pp. 247–263.
13. P. Košovan, T. Richter and C. Holm, Macromolecules, 2015, 48, 7698–7708.
14. G. S. Longo, M. O. de la Cruz and I. Szleifer, Macromolecules, 2011, 44, 147–158.
15. I. Adroher-Bentez, S. Ahualli, A. Martn-Molina, M. Quesada-Pérez and A. Moncho-Jordá, Macromolecules, 2015, 48, 4645–4656.
16. R. Schroeder, A. a. Rudov, L. A. Lyon, W. Richtering, A. Pich and I. I. Potemkin, Macromolecules, 2015, 48, 5914–5927.
17. Z. Posel, Z. Limpouchová, K. Šindelka, M. Lísal and K. Procházka, Macromolecules, 2014, 47, 2503–2514.
18. P. K. Jha, J. W. Zwanikken, J. J. de Pablo and M. O. de la Cruz, Curr. Opin. Solid State Mater. Sci., 2011, 15, 271–276.
19. P. K. Jha, J. W. Zwanikken, F. A. Detcheverry, J. J. de Pablo and M. Olvera de la Cruz, Soft Matter, 2011, 7, 5965–5975.
20. P. Gong, J. Genzer and I. Szleifer, Phys. Rev. Lett., 2007, 98, 018302.
21. R. Nap, P. Gong and I. Szleifer, J. Polym. Sci., Part B: Polym. Phys., 2006, 44, 2638–2662.
22. G. S. Longo, M. O. de la Cruz and I. Szleifer, J. Chem. Phys., 2014, 141, 124909.
23. T. Colla, C. N. Likos and Y. Levin, J. Chem. Phys., 2014, 141, year.
24. A. Moncho-Jorda and J. Dzubiella, Phys. Chem. Chem. Phys., 2016, 18, 5372–5385.
25. G. J. Fleer, M. A. Cohen Stuart, T. Cosgrove and B. Vincent, Polymers at interfaces, Chapman and Hall, London, 1993.
26. J. Klein Wolterink, F. A. M. Leermakers, G. J. Fleer, L. K. Koopal, E. B. Zhulina and O. V. Borisov, Macromolecules, 1999, 32, 2365–2377.
27. J. Klein Wolterink, J. van Male, M. A. Cohen Stuart, L. K. Koopal, E. B. Zhulina and O. V. Borisov, Macromolecules, 2002, 35, 9176–9190.
28. O. V. Rud, A. A. Mercurieva, F. A. M. Leermakers and T. M. Birshtein, Macromolecules, 2012, 45, 2145–2160.
29. F. Uhlík, P. Košovan, Z. Limpouchová, K. Procházka, O. V. Borisov and F. A. M. Leermakers, Macromolecules, 2014, 47, 4004–4016.
30. F. Uhlík, P. Košovan, E. B. Zhulina and O. V. Borisov, Soft Matter, 2016, 4846–4852.
31. J. Höpfner, C. Klein and M. Wilhelm, Macromol. Rapid Commun., 2010, 31, 1337.
32. A. Katchalsky, S. Lifson and J. Mazur, J. Polym. Sci., 1953, 11, 409–423.
33. B. A. Mann, C. Holm and K. Kremer, Macromol. Symp., 2006, 237, 90–107.
34. B. A. Mann, O. Lenz, K. Kremer and C. Holm, Macromol. Theory Simul., 2011, 20, 721–734.
35. P. Košovan, T. Richter and C. Holm, Intelligent Hydrogels, Springer International Publishing, 2013, vol. 140, pp. 205–221.
36. G. S. Manning, J. Chem. Phys., 1969, 51, 924–933.
37. A. Griewank and A. Walther, Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, SIAM, 2nd edn, 2008.
38. J. Scheutjens and G. Fleer, J. Phys. Chem., 1979, 83, 1619–1635.
39. R. Israëls, F. A. M. Leermakers and G. J. Fleer, Macromolecules, 1994, 27, 3087–3093.
40. R. Israëls, F. A. M. Leermakers, G. J. Fleer and E. B. Zhulina, Macromolecules, 1994, 27, 3249–3261.
41. M. Doi and S. F. Edwards, The theory of polymer dynamics, Oxford Science Publications, 1986.
42. P. G. de Gennes, Scaling Concepts in Polymer Physics, Cornell University Press, Ithaca, NY, 1979.
43. O. V. Borisov, E. B. Zhulina, F. A. Leermakers, M. Ballauff and A. H. E. Müller, Self Organized Nanostructures of Amphiphilic Block Copolymers I, Springer Berlin Heidelberg, 2011, vol. 241, pp. 1–55.
44. T. M. Birshtein, A. A. Mercurieva, F. A. M. Leermakers and O. V. Rud, Polym. Sci., Ser. A, 2008, 50, 992–1007.
45. M. Daoud and J. Cotton, J. Phys., 1982, 43, 531–538.

---

Footnote

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

---