Interactions between Polyelectrolyte Gel Surfaces - Macromolecules

Nov 17, 2016 - Cosolute Partitioning in Polymer Networks: Effects of Flexibility and Volume Transitions. Won Kyu Kim , Arturo Moncho-Jordá , Rafael R...
1 downloads 9 Views 4MB Size
Download PDF
Article pubs.acs.org/Macromolecules

Interactions between Polyelectrolyte Gel Surfaces Aykut Erbaş† and Monica Olvera de la Cruz*,‡ †

Department of Materials Science and Engineering and ‡Department of Chemistry, Department of Chemical and Biological Engineering, and Department of Physics, Northwestern University, Evanston, Illinois 60208, United States ABSTRACT: Surfaces formed by charged polymeric species are highly abundant in both synthetic and biological systems, for which maintaining an optimum contact distance and a pressure balance is paramount. Here we investigate interactions between surfaces of two same-charged and highly swollen polyelectrolyte gels, using extensive molecular dynamic simulations and minimal analytical methods. The external-pressure responses of the gels and the polymer-free ionic solvent layer separating two surfaces are considered. Simulations confirmed that the surfaces are held apart by osmotic pressure resulting from excess charges diffusing out of the network. At a threshold pressure, the counterion-induced osmotic pressure in the gel interiors becomes nearly equal to that in solvent gap, and the gels begin to deform. At the threshold pressure, the distance between surfaces is equivalent to the electrostatic screening length imposed by excess charges. Both the solvent layer and pressure dependence are well described by an analytical model based on the Poisson−Boltzmann solution for low and moderate electrostatic strengths. Scaling descriptions for the threshold pressure and critical distance compare well with the simulations. The threshold pressure decreases with increasing polymerization degree of network chains, similar to sparsely grafted polyelectrolyte brushes. As the electrostatic coupling strength is increased, the deviations from the analytical model also increase due to condensed counterions. Our results are of great importance for systems where charged gels or gel-like structures interact in various solvents, including systems encapsulated by gels and microgels in confinement.

■

INTRODUCTION Equilibrium swelling of polyelectrolyte hydrogels is described by the balance of elastic energy of the network chains constituting the gels and the osmotic pressure of ionized counterions. This is a Donnan equilibrium as long as the electrostatic interactions are not stronger than the thermal energy.1,2 The Donnan equilibrium limits the translational entropy of counterions and confines the ions in the gel, preserving local electroneutrality. For a finite-size gel, the concentration difference of ions inside and outside of the gel creates a potential difference and an outward counterion layer at the gel−solvent interface. This is mainly due to dominance of the counterions’ entropy over electrostatic attractions between backbone charges and ionized counterions near the interface. If two finite-size gel or gel coated-surfaces are brought into close proximity, the surfaces are physically kept apart by the osmotic pressure of the excess ions. Thus, a polymer-free ionic solvent layer occupying the space between two surfaces3−5 can support relatively high external pressures until an osmotic-pressure balance between interior and exterior of the gels is reached. If counterions are confined between two charged impenetrable plates, any decrease in interplate distance (i.e., compression) increases the counterion concentrations in the space bounded by the plates. The increase in ionic concentration gives rise to a pressure. If solid plates are replaced by two gel surfaces, counterions can diffuse in and out of the gels. As the distance between surfaces is decreased, a © XXXX American Chemical Society

different picture can emerge as a result of the intricate balance between entropic and electrostatic components; compression can deform the gels (if the solvent content in the gels is allowed to squeezed out), or the counterion concentration throughout the entire system can adapt a new state. Indeed, interactions between two hydrogel surfaces under normal loads have been considered in various theoretical frameworks4−7 and also in the context of polyelectrolyte brushes.8−10 The osmotic pressure and the thickness of the solvent layer have been calculated as a function of external load at the Debye−Huckel level. Indirect experimental comparisons have also been achieved through sliding experiments.4,11 Indeed, a direct validation of the analytical descriptions is possible via computational simulations. Surprisingly, however, although equilibrium properties of hydrogels have been studied extensively,12−17 and direct validation of the analytical descriptions is possible via computational simulations, only one computational study has been reported on the interaction between gel surfaces, to our knowledge.18 However, that study lacks of any comparison. Moreover, a detailed picture of electrostatic strength on the interactions between gel surfaces and the resulting pressure regimes are missing in the literature even though relevant systems are quite abundant in synthetic and biological systems. Received: July 1, 2016 Revised: November 1, 2016

A

DOI: 10.1021/acs.macromol.6b01416 Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

■

RESULTS AND DISCUSSION Theoretical Considerations. Consider two gel surfaces separated by a polymer-free gap with thickness δ as shown in Figure 1. The center of the gap is located at z = 0. The gap and

Apart from the relevance of surfaces coated with charged gels to characterize the electrostatic-originated interactions in soft systems, such surfaces have great potential relevance in various material applications: for instance, in biomimetic systems that include soft components operating in solvents19−21 or in microgel or gel coated capsules confined at high concentrations or interacting with identical surfaces.22−24 Particularly at small separations, these systems can be modeled as two planar surfaces, for which curvature effects can be negligible. Furthermore, cracks observed in multicomponent brittle hydrogels are also due to formation of similar surfaces in bulk structures under various stimuli.25,26 Surfaces coated with gel-like structures are also abundant in biological machinery, particularly in connective and intercellular tissues, where optimum contact between surfaces is paramount. Cartilage tissue and its components,27,28 mucin bilayers bridging the cornea stroma and epithelial surfaces,29,30 and a mucin-brush layer coating the interior surfaces of the mammalian respiratory system31 are among some examples. It is thus desirable to understand the interaction between these soft, charged surfaces and possibly relate them to the chemical properties of the polymer network (e.g., swelling ratio, charge fraction) and electrostatic strength of solvent. To better understand the interactions between these soft, charged surfaces, in this work we revisit a previous theory based on the solution of Poisson−Boltzmann (PB) equation.4 This theory is extended with scaling arguments to relate the electrostatic strength of the solvent and polymeric properties of gels. Polymerization degree, the Kuhn segment size (i.e., persistence length), and electrostatic strength of solution are related to the pressure regimes and a threshold pressure separating these two regimes. At this threshold pressure, the distance between the two gel surfaces reaches an effective Gouy−Chapman length. The osmotic pressures of the polymer-free layer separating the two gels and that of the gel interiors become equal. At higher pressures, the elastic energy of constituting chains begin to deform the gels. These theoretical descriptions are compared to the strain-controlled deformation simulations of hydrogel bilayers composed of two identical hydrogel slabs. The simulations were performed for various network-chain sizes and electrostatic strengths by employing a coarse-grained molecular dynamics (MD) model. We found that interactions between the same-charge gels are mainly dominated by counterion induced repulsive forces at noncontact distances for weak electrostatic couplings. The separation distance can be described well by the PB-based solution in a pressure range spanning 3 orders of magnitude. We also found that the distance between two gel surfaces at a given pressure is strongly affected by the electrostatic properties of the medium and exhibits two distinct regimes in accord with the theory. The scaling arguments can predict the polymerization-degree dependence of the threshold pressure (intersurface distance) separating low- and high-pressure regimes well within the parameter ranges considered here. In the following, we first present a theoretical description of pressure dependence of the solvent layer confined by surfaces of two highly swollen gels. After describing the simulation methodology, we then compare the analytical model to the coarse-grained simulations of strain-controlled gel deformations. Finally, we discuss the applicability and the limits of the PB solution. Hierarchical deformation of gels and the separation layer and the effect of dangling ends are also considered in the Discussion section.

Figure 1. Separation gap between two highly swollen hydrogel surfaces. Effective surface charge density σeff due to charged network chains (green beads) and condensed counterions can be defined at the interface. The ρgap denotes the counterion (purple beads) density at the center of the two gels at z = 0. The contributions to the osmotic pressure are shown within the corresponding volume.

gels are in mechanical equilibrium. The chemical details of flexible network chains are described by a Kuhn segment size (persistence length) b and a charge fraction factor f. Thus, the distance between two charged monomers along the backbone is d = b/f. The number of Kuhn segments between two crosslinks is N. The concentration of counterions within the gel is cci ≈ Nf/L03, where the swollen size of where the swollen size of an inter-cross-link chain in chain in good solvent is given by L0 ≈ bNf1−ν 12,13 with exponent ν ≈ 3/5. The surface normal of gel slabs are taken parallel to the z-direction. If the polymer networks are compressed in the z-direction, the diagonal component of the stress tensor in the deformation direction, referred to here as Pz, will be measured throughout the system. In mechanical equilibrium, the pressure balance is

Pz ≈ Πgap ≈ Πgel

(1)

where Πgap and Πgel are free energy densities in the gap and inside the gels, respectively. Ideally, Πgap and Πgel should include entropic contributions, steric interactions, and electrostatic energy resulting from charge−charge correlations. However, for highly swollen gels (i.e., in a dilute regime) and for low electrostatic strengths, we may write Pz ≈ Πgap ≈ Πci − Πel

(2)

where Πci and Πel are the counterion osmotic pressure and elastic energy density of network chains, respectively (see Figure 1). In eq 2, while Πel favors the collapse of network chains, Πci favors the expansion of gels. The pressure inside the gels is balanced by the outer osmotic pressure, Πgap. In dilute regime, the main contribution to the free energy in the polymer-free gap comes from the translational free energy of counterions Pz ≈ Πgap ≈ kBTρgap B

(3) DOI: 10.1021/acs.macromol.6b01416 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

center of the gap becomes ρgap = ρ*gap. At smaller distances b ≪ δ < δ*, the entire gel system (gels + layer) enters a semidilute regime, in which osmotic pressure is mainly dominated by counterion concentration.32 By using eqs 9 and 10, a more compact relation for eq 8 can be obtained as

where kB is the Boltzmann constant, T is temperature and ρgap is the counterion density at the geometrical center of the separation gap (Figure 1). Indeed, eq 3 is in agreement with the simulations that we will discuss in the next section as well as with the previous work.13 An analytical expression for ρgap can be obtained by solving the Poisson−Boltzmann (PB) equation for the separation gap confined by two gel surfaces (see Figure 1 for the schematic description). The PB equation for the counterion density between two charged planar surfaces can be solved if the charge of surface-exposed chains is replaced by an effective surface-charge density σeff.4 In this case, the PB equation for the interior of the gap reads ∇Ψ̃(z) = −4π SBρgap exp[−Ψ̃(z)]

δ = δ*

where Ψ̃(z) ≡ eΨ/kBT is the rescaled electrostatic potential defined within −δ/2 < z < δ/2, and the elementary charge is e. In eq 4, the strength of electrostatic interactions is characterized by the Bjerrum length SB = e2/4πε0εkBT, where ε0 is the vacuum permittivity and ε is the relative permittivity of the medium. One of the solutions for eq 4, with the boundary condition at the center of the gap ∂zΨ̃(z)|z=0 = 0, is (5)

σeff ≈

The total number of counterions trapped between two gel surfaces is determined by the effective surface charge density induced by the chains adjacent to the interface. Hence, within −Ψ̃(z) the gap, this charge, σeff = ρgap∫ δ/2 dz, is −δ/2e σeff =

2ρgap π SB

⎛ Sρ π B gap tan⎜⎜δ 2 ⎝

⎞ ⎟ ⎟ ⎠

(6)

δ=

⎞ ⎟⎟ ⎠

δ* ≈

Pz* ≈

(8)

* = kBTσeff 2 SB Pz* ≈ kBTρgap (9)

and similarly, a crossover value for the separation gap is δ* =

1 σeff SB

Nb3/2 SB1/2

(13)

kBT N 2b3

(14)

Note that eq 14 does not have any explicit SB dependence for weak coupling (i.e., f SB/b < 1). However, for strong electrostatic couplings (i.e., f SB/b > 1), the concentration of ionized counterions can decrease. As a result, the counterion-induced osmotic pressure Πci decreases, and the elastic energy of network chains Πel shrinks the gel. Hence, for systems, in which charge condensation is strong (f SB/b ≫ 1), due to electrostatic correlations and resulting steric interactions, eq 7 will not describe the corresponding systems. Alternatively, in obtaining eq 14, the PB equation can be solved for the electrostatic potential inside the hydrogels to obtain a relation for σeff.5 However, if the overall height of the hydrogel slabs is greater than the thickness of the gap, i.e., h ≫ δ, then σeff can be approximated successfully by eq 12 as we will show in the next section. In the following section, we will compare the results derived in this section to the coarse-grained simulations with explicit electrostatic interactions.

By equating two limits, the crossover values for the pressure is

* = σeff 2 SB ρgap

(12)

Equation 13 predicts that δ* decreases with increasing SB and increases with N, which hint that uncondensed counterions is the determining factor for threshold thickness. A scaling relation for P*z , by using eqs 9 and 12, can be obtained as

(7)

Limits of eq 7 for weak and high normal pressures for a fixed σeff and SB are after neglecting prefactors. 1/2 ⎧ for Pz ≪ Pz* ⎪ (Pz*/ Pz) δ ≃⎨ δ* ⎪ for Pz ≫ Pz* ⎩ (Pz*/Pz)

1 Nb3/2 SB1/2

In eq 12, the term for the charge fraction is dropped for simplicity since 0 < f ≲ 1. However, the lower limit of f should be sufficiently large enough to provide stretched chains (i.e., L0 ∼ N). In eq 12, σeff decreases with increasing SB as a result of electrostatic attractions between the counterions and backbone charges. Using eqs 10 and 12, one can also write

Extracting ρgap from eq 6 and plugging the result into eq 3 leads to a relation for the thickness of the separation gap as a function of the normal pressure, Pz.4 ⎛ πkBT SB 2kBT arctan⎜⎜σeff πPz SB 2Pz ⎝

(11)

What remains is to relate Pz* and δ* in eqs 9 and 10 to the polymeric properties of hydrogels (i.e., N, f, b). As the hydrogel surfaces are brought closer, the density of counterions inside the separation gap increases. At a threshold value of the gap thickness, the osmotic pressures of counterions inside and outside of hydrogels become nearly equal (i.e., Πci ≈ Πgap). Note that this ansatz is confirmed in the simulations and discussed in the next section. According to eq 2, above this threshold, the elastic energy Πel dominates, and the gels begin to deform under the external pressure. Under the assumption that at the onset of deformation, the counterion osmotic pressure inside gap is given by Πci ≈ kBTcci ≈ kBTNf/L03, the pressure balance in eq 2 with the help of eq 9 leads to

(4)

Ψ̃(z) = 2 log[cos(z 2π SBρgap )]

⎛ P* ⎞ Pz* arctan⎜⎜ z ⎟⎟ Pz ⎝ Pz ⎠

(10)

Here, δ* is analogous to the Gouy−Chapman length that characterizes the width of the counterion layer adjacent to a charged surface. The electrostatic forces due to the surface are screened beyond the Gouy−Chapman length. At intersurface separations δ > δ*, two counterion layers do not interact. At separations, δ = δ* (i.e., at Pz = Pz*), the Gouy−Chapman layers of two surfaces overlap, and the counterion density at the C

DOI: 10.1021/acs.macromol.6b01416 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules Comparison with Simulations. Simulation Methods. In the coarse-grained simulations with explicit counterions in implicit good solvent, one hydrogel slab is placed on a second identical slab, as shown in Figure 2. Each hydrogel is

UC(r ) = −ϵ

for r < re

The electrostatic cutoff distance is re = 6σ, above which longerrange electrostatic interactions are calculated via the particle− particle−particle mesh (PPPM) Ewald solver with an error tolerance 10−5 and 10−6 for N = 32, 64 and N = 128 hydrogel bilayers, respectively, for computational efficiency. The Bjerrum length quantifies the length scale, above which Coulomb energy between two charges is less than thermal energy kBT. The strength of the electrostatic interactions in simulations is adjusted by tuning the dielectric constant of the medium to obtain SB ≃ 1b, 4b, and 8b. For a flexible polymer, the case with SB ≃ 1b can correspond to water (i.e., SB ≈ 7 Å). The cases SB ≃ 4b, 8b may represent ethyl alcohol at room temperature (SB ≈ 30 Å) and butyl alcohol at 320 K (SB ≈ 60 Å), respectively.35 All MD simulations were run using the Lammps MD package.36 During the relaxation of initial configurations, a twostep procedure was followed: at the initial stage, all electrostatic interactions were switched off, and each system was simulated for 105−107 MD steps at a constant pressure P = 0, particle number Nm, and temperature T. The system temperature was set to T = 1.0ϵ/kB via a Nosé−Hoover thermostat. For the second stage, electrostatic interactions were switched on, and the systems were allowed to relax another 105−107 MD steps before the data production runs are performed. Simulations were run with a time step of Δt = 0.002τ. Here the unit time scale τ = σ m/ϵ is the LJ time step, where monomeric LJ mass m is unity. Monomers of the network chains and counterions were of the same size. To construct the hydrogel slabs, a single polymeric unit lattice was repeated three times in the ẑ direction and four times in x̂ and ŷ directions as shown in Figure 1. While the slabs are finite in the ẑ-direction, in the lateral directions each terminal monomer is bonded to the terminal monomer at the other side of the same slab to ensure lateral periodicity of the hydrogel networks. Periodic boundary conditions for the simulation boxes were introduced in all directions. Hence, our gel slabs can be considered to extend to infinity in the lateral directions, whereas in the vertical direction, there is a stack of hydrogel slabs separated by polymer-free gaps. Alternatively, our model system can be considered as a hydrogel system enclosed by a semipermeable membrane that can allow solvent to flow in and out of the gels, but not the counterions. Equilibrium Swelling of Hydrogels. This is defined at zero pressure so that the hydrogels neither swell nor collapse. The hydrogel bilayers are composed of two unconnected hydrogel slabs (Figure 1), and there is no force holding two slabs at close proximity. Ideally hydrogel slabs separate until the value of vertical pressure component reaches Pz = 0. In simulations, we allowed two slabs to separate in the vertical direction until value of each pressure component reaches the order of error bars that we can achieve here, i.e., P ≈ Px ≈ Py ≈ Pz ≈ 0.0 ± δP, where the error tolerance δP ≈ 10−5−10−6kBT/σ3 for N = 32, 64 and N = 128 systems, respectively. This criteria resulted in vertical box heights of approximately D0 ≈ 12L0 ≈ 4h0, where h0 is the initial height of each slab. At this box height, each hydrogel was isotropically swollen within the error bars. Hence, we accepted D0 ≈ 12L0 as a reference point for our deformation simulations. The largest system size considered here is over 1000σ for N = 128 gels.

Figure 2. Illustration of the simulated polyelectrolyte−gel systems composed of two identical slabs with f = 0.5 and N = 64 at the initial state. The chain and counterion monomers are rendered as green and purple beads, respectively. For clarity, on the left, only network chains chains are shown; on the right, both network and counterions are shown. At the initial configuration, two gel slabs are separated by a gap with a thickness δ0 ≈ h0 (in simulations, it is larger than illustrated here and δ0 ≈ h0). The dashed lines show the actual size of the simulation box. The snapshot is obtained via VMD.

constructed as a defect-free cubic polymer lattice. The lattice junction points (i.e., cross-links) are permanent. The chains connecting two cross-links of the hydrogel networks are modeled, using the coarse-grained “Kremer−Grest (KG)” bead−spring model.33,34 Nonbonded interactions are accounted for a shifted Lennard-Jones (LJ) potential ⎧ ⎡⎛ ⎞12 ⎛ ⎞6 1⎤ σ σ ⎪ ⎪ 4ϵ⎢⎜ ⎟ − ⎜ ⎟ + ⎥ for r < rc ⎝r⎠ 4⎦ ULJ(r ) = ⎨ ⎣⎝ r ⎠ ⎪ ⎪0 for r > rc ⎩

where the energy unit ϵ = kBT is the strength of LJ potential, and the unit length σ is the monomer size. In eq 15, r = |r| is the distance between two monomers, and rc = 21/6σ is the cutoff distance, which provides a repulsive LJ force between monomers. Two adjacent chain monomers are bonded via a finite extensible nonlinear elastic (FENE) potential ⎛ 1 r2 ⎞ UFENE(r ) = − kr0 2 log⎜1 − 2 ⎟ 2 r0 ⎠ ⎝

SB r

for r < r0

where k = 30ϵ/σ2 is the spring constant; the maximum bond stretching length is r0 = 1.5σ. The numbers of monomers per linear network chain connecting two cross-links are N = 32, 64, and 128 to obtain various gel bilayers with equilibrium volume fractions ϕ ∼ N/ L03. In the flexible-chain model, the Kuhn size is b ≈ σ, and the charge fraction is set to f = 0.5. The overall system is electroneutral. That is, for each positive charge on the backbone, one monovalent counterion is added into simulation box at a random position. In the simulations, short-range Coulomb interactions between all charges are calculated via D

DOI: 10.1021/acs.macromol.6b01416 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules Simulations of Strain-Controlled Gel Deformations. The vertical height of the simulation box, D0, was reduced to a prescribed height D to obtain a deformation ratio D/D0. The subscript “0” is dropped to donate all quantities at D < D0. During the compression process, no pressure coupling in the corresponding direction was applied while for lateral, i.e., x̂ and ŷ, directions, the components of pressure barostat were set to Px = Py = 0. Once the prescribed strain was reached, the hydrogel bilayers were equilibrated for another 105−106 MD steps before collecting data for analyses. The pressure-coupling constant of Nosé−Hoover thermostat was set to τp = 1.0τ−1. In our deformation simulations, the normal pressure component was always Pz ≫ Px ≈ Py ≈ 0.0. Unless noted otherwise, all results presented in this paper were averaged over time. Error bars are not shown if they are smaller than the size of the corresponding data point. Various numerical analyses on the MD trajectories were done by employing in-house analysis programs. Strain-Free Swelling. In the simulations, each identical slab is equally swollen. For the Bjerrum length SB = b, the swelling ratios of the gels obey a power law λ0 ≡ L0/Rdry ∼ N1/2, where the size of a network chain in the unswollen state is taken as Rdry = bN1/2 (see Table 1 for swelling ratios). Since longer Table 1. Table of Equilibrium Swelling Ratios, λ0, and All Threshold Values Used To Rescale Figures 4 and 6a N

SB /b

λ0b

ρgap * [σ−3]

δ* [σ]

Pz* [ϵ/σ3]

σeff [σ−2]

32 32 32 64 64 128 128

1 4 8 1 4 1 4

3.5(3) 2.9(3) 1.8(7) 5.0(0) 3.8(7) 7.0(6) 5.6(0)

0.0064 0.0110 0.0003 0.0013 0.0027 0.000 26 0.000 57

10.1 3.4 21 27 13 59 29

0.0099(6) 0.0213(8) 0.0002(9) 0.0014(1) 0.0014(3) 0.00028(8) 0.00028(2)

0.099(8) 0.073(1) 0.006(0) 0.0376(9) 0.0188(7) 0.0169(7) 0.0084(4)

Figure 3. Normalized density profiles of counterions of a uniaxially deformed (compressed) N = 32 (left column) and N = 128 (right column) hydrogel bilayers as a function of rescaled vertical distance for f = 0.5 and SB /b = 1, 4, and 8 cases. Note that each system has a different D value.

respect to the cases with lower SB values since most ions are attracted by backbone charges. The increasing counterion accumulation near the network chains suggests a Manning-condensation type behavior, which is valid for infinitely long charged rods at infinite dilution.37,38 Accordingly, at high electrostatic strengths (f b/SB > 1), a fraction of ionized counterions “condenses” on the chains and reduces the 1D charge density of the chains from Γ = 1/d = f/b to ΓOM = 1/SB. As the SB is increased, more counterions are accumulated near the network chains, and there are fewer free ions contributing to the counterion osmotic pressure. Hence, the swelling ratio of gels is lower for SB > b. Therefore, specifically for the strongest electrostatic strength SB = 8b, the gels are weakly swollen, and the network chains are much less stretched.14,39 Strain-Controlled Deformation. Using the configurations obtained from P ≈ 0.0 (i.e., D = D0) simulations, a set of uniaxial deformation simulations were performed as described above. As the compression ratio D/D0 is decreased, the counterion concentration between two gels increases as shown in Figure 3. To quantify the density changes, in Figure 4a we plot the counterion density at the center of the gap (filled symbols) and inside the hydrogels (open symbols) as a function of D/D0 for N = 32 gels. For gel interiors, the counterion concentrations are obtained by averaging the minima at the centers of the gel slabs in the density profiles given in Figure 3. Each density is rescaled by its corresponding value of ρgap * to obtain a crossover at ρgap/ρgap * = 1 (see Table 1). Figure 4a reveals that at around compression ratios D/D0 ≈ 0.5 for N = 32 gels the counterion concentration in the gap and inside the gels becomes roughly equal. If the gel system is further compressed (i.e., D/D0 < 0.5), the concentration of

a The values for δ* and Pz* are obtained via eqs 9 and 10 by using fitted σ*eff values, whereas ρgap is obtained via concentration profiles. bλ0 ≡ L0/Rdry, where Rdry = bN1/2 is dry size of a chain in melt conditions, and the Kuhn size is b = 1σ.

chains are more flexible (i.e., the elastic energy density, on the leading order, scales as Πel ∼ N−1), the counterion osmotic pressure can stretch the longer network chains further. For moderate strength of electrostatics (SB = 4b), the swelling ratios also obey λ0 ∼ N1/2; however, they are lower than the SB = b case. For the shortest chain size we have here, N = 32, deviation from the scaling is observed due to the prevalence of elastic energy over counterion osmotic pressure. At the initial configurations (i.e., D = D0), the majority of the counterions are located in the gels. This situation is more noticeable for N = 128 hydrogels as can be seen in the top frames of Figure 3. Here we show the normalized density profiles of counterions throughout the entire simulation box, ρ̃ci, as a function of rescaled vertical distance for N = 32 and N = 128 gels. All densities are normalized so that D−1∫ ρ̃ci(z) dz = 1. In Figure 3, the center of the gap corresponds to z/D ≈ 0.5. The undulations in the density profiles in Figure 3 are due to the counterions near the network chains and cross-links. A visual inspection reveals that the difference between maxima and minima in the profiles is higher for SB = 4b (blue curves in Figure 3) compared to SB = b cases (red curves). For the strongest electrostatic strength we consider here (SB = 8b), the undulations in the density profile are significantly irregular with E

DOI: 10.1021/acs.macromol.6b01416 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

before two gels are at physical contact, which was also predicted for polyelectrolyte brushes of flexible chains.9,40,41 The assumption in eq 3 that counterions in the gap can be treated as ideal gas can be assessed by the data in Figure 4a,b. In Figure 5, we show the ratio Pz/ρgap for various gel systems. This

Figure 5. Ratio of the ẑ-component of the stress tensor, Pz, and the counterion density at the center of the solvent free gap, ρgap, for various gels and electrostatic strengths. The arrow indicates the direction of increasing electrostatic strengths. The shaded area highlights the data for SB /b = 1 case.

ratio corresponds to the osmotic coefficient in polymer solutions that quantifies the deviations from ideal gas behavior. For weak electrostatic strengths (i.e., SB/b = 1), the ratios are nearly independent of deformation ratio for various gels (see the data in the rendered area in Figure 5). Indeed, this behavior at δ < δ* (or for Pz > P*z ) is analogous to that observed in semidilute solutions of polyelectrolytes, in which the osmotic coefficient exhibits a plateau with polymer concentration.32,42 As the gels and the solvent-free layer are strongly compressed, the average monomer concentration in the overall system increases toward a highly concentrated regime. Consistently, at strong deformations (D/D0 → 0), regardless of the system size or electrostatic strength, the steric interactions begin to contribute to the pressure (Figure 5). As for SB = 4b and 8b cases, however, Pz/ρgap deviates from unity: this indicates that electrostatic correlations start contributing to the pressure negatively as the gels are drawn closer. In Figure 6a, the rescaled thickness of the gap, δ/δ*, is shown as a function of the rescaled normal pressure, Pz/P*z , for various gels and electrostatic strengths. The thickness of the separation gap, δ, is defined as the root-mean-square distance between the center-of-masses of the cross-links adjacent to the separation gap (see also Figure 1 for a schematic definition). Almost all data points overlap, particularly for low-electrostatic cases (i.e., SB/b = 1, 4) at around (Pz/Pz* ≈ 1, δ/δ* ≈ 1). In the simulations, the general trend is that the higher the electrostatic strength is the narrower the separation layer is. The smaller values of δ with increasing SB values are due to the counterions that evacuate the separation gap to condense on the chains: as more counterions leave the gap, the counterion-induced osmotic pressure in the gap decreases, and the gap shrinks. Note that the leading order terms in eq 8 include either SB or σeff. In the inset of Figure 6a, the same data are shown in log−log scale to reveal two pressure regimes discussed in eq 8 for low-

Figure 4. (a) Average counterion densities inside the gap (filled symbols) and inside the gels (open symbols) as a function of compression ratio D/D0 for N = 32 hydrogel bilayer for various SB values. The densities are rescaled to obtain ρgap/ρ*gap = 1 at the * ). (b) Measured normal crossover (see Table 1 for values of ρgap pressures of a uniaxially deformed N = 32 and N = 128 (inset) hydrogels systems. (c) Height of the hydrogel slabs normalized by their values at D = D0 (see Figure 2 for schematic definition). Vertical lines highlight the deformation ratio at which gels start to deform under external pressure for N = 32 gels.

counterions in and outside of the polymer-free gap is indistinguishable. This behavior is persistent for various electrostatic strengths (Figure 4a). We observe similar regimes for N = 64, 128 gels, however, at smaller compression ratios due to lower volume fractions of counterions. As the counterion density in the gap approaches the crossover density, ρ*gap, according to eqs 8 and 9 a crossover for pressure sets in. Indeed, as shown in Figure 4b, the crossovers for densities observed in simulations also correspond to the crossovers between two distinct pressure regimes (see vertical dashed lines in Figure 4a,b). Moreover, at this compression ratio we observe crossovers for the density and pressure, the gels themselves begin to deform significantly as shown in Figure 4c, where the normalized heights of various gels are plotted as a function of compression ratio D/D0. This confirms that the onset of gel deformation is set by the ionic concentration difference inside and outside of gels (i.e., Πci ≈ Πgap). Importantly, deformation of the gel networks is initiated F

DOI: 10.1021/acs.macromol.6b01416 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

be stronger and render the constant surface charge approximation invalid. In such situations, the PB equation fails to describe the regimes, where the pressure regime is determined by charge correlations. In Figure 6b, the threshold pressure and effective surface charge density obtained from the fits are plotted as a function of N for various electrostatic strengths. For weak electrostatic coupling (red data in Figure 6b), the scaling predictions in eqs 12 and 14 capture the decreases of σeff and Pz* ≈ σeff2SB with increasing N (solid lines in Figure 6b). The slopes P*z ∼ N−2.2±0.1 and σeff ∼ N−1.1±0.05 seem to describe the data slightly better. To validate the fitting further, the values of the crossover density ρgap * obtained directly from the concentration profiles are compared to the density values imposed by fitting eq 7 to the data via P*z ≈ ϵρ*gap with ϵ = 1. While for the case of SB/b = 1, agreement is very good (see Table 1), for the moderate electrostatic strength (i.e., SB/b = 4) both the P*z and σeff are slightly underestimated by eq 7. As mentioned earlier, for SB/b = 4, the ideal gas approximation already underestimates the pressure (see the inset of Figure 4a) due to electrostatic interactions between counterions and backbone charges.14 Also, note that in more collapsed gels, the surface charge density can recover its mean-field value and results in a steeper slope (i.e., σ2d ∼ cci2/3 ∼ N−4/3).

■

DISCUSSION In our simulations, we observed an nonaffine deformation throughout the simulation boxes below the threshold pressures. That is, for weak pressures (Pz < P*z ), the deformation of the gels is insignificant (Figure 4c), but the polymer-free gap deforms considerably. At around Pz ≈ P*z , both the gels and the gap begin to deform. Indeed, the osmotic modulus is the change in osmotic pressure with respect to monomeric density. At weak compression, the solvent gap with a lower monomer density deform substantially with compression (i.e., δ ∼ (D/ D0)2.3±0.1), but the gel slabs are not deformed (Figure 4c). As the gel surfaces are further squeezed, the counterion density inside the gap increases. At the threshold pressure, osmotic pressures of the gap and gels approach each other. Consequently, at higher pressures, both gel networks and polymer-free gap deform simultaneously. In the gelation process of hydrogels, gel networks are not defect free. Inhomogeneous cross-linking along the chains results in “dangling” chains with at least one end unbound to the network. Any dangling chain can have a monomer number 1 ≤ Nd ≤ N. At interfaces, these free ends can modify the interactions between two gels as well. An estimation for the contribution of the dangling chains to the normal pressure can be obtained by treating these chains similar to sparsely grafted polyelectrolyte brushes. If the distance between grafting points is to be d ∼ N for swollen gels (i.e., the dangling ends only spring out from cross-links), the disjoining pressure at separations, below which dangling chains of two opposing brushes overlap, is Pz ∼ f Nd/d2δ.8 Indeed, this scaling coincides with the high pressure limit given in eq 8 as Nd → N and d → N (i.e., Pz ≃ δ*P*z /δ ∼ 1/Nδ). This upper limit indicates that the contributions due to the dangling ends do not change the results presented here in a qualitative manner for swollen gels. However, they are rather expected to modify prefactors in the scaling arguments. If the distance between two hydrogels is larger than the threshold thickness δ > δ*, two surfaces interact via the tails of

Figure 6. (a) Rescaled thickness of the gap between two gels as a function of rescaled normal pressure for different compression ratios obtained from simulations for various gel systems. The values for rescaling pressure and thickness are given in eqs 9 and 10, respectively (see also Table 1). The dashed line is the PB solution given in eq 11. (b) Threshold pressure as a function of number of Kuhn monomers between two cross-linkers in log−log scale. Inset shows the values of σeff obtained for various gels via fitting the simulation data to eq 7.

pressure (weak compression) and high-pressure (strong compression) limits: as the pressure is increased, the δ/δ* decreases with a slope of −1 as predicted (see eq 8). Below the threshold pressure, the dependence of the thickness on the pressure is weaker (δ/δ* ∼ Pz−0.5) as shown in the inset of Figure 6a. The simulation results are compared to the analytical model (eq 7) in Figure 6a. To fit the data, effective surface charge density, σeff, in eq 7 is released as a fit parameter. The fit values for N = 32, 64, and 128 gels are shown in Table 1. The values for Pz* and δ* are obtained using fitted σeff values according to eqs 9 and 10 for each case separately. In Figure 6a, the gap thickness, δ, can be described well for a wide range of pressures by eq 7 (dashed curve in Figure 6a). Equation 7 describes the pressure dependence of thickness of the separation layer well, again specifically at low and moderate pressures. However, eq 7 fails to describe the data at high pressures for SB = 8b (green triangles in Figure 6a). At high electrostatic strengths, almost all counterions condense on the chains (i.e., ρgap → 0), and the two gels interact directly chain by chain. Indeed, in Figure 6a, a plateau appears for the SB = 8b case at around Pz/P*z ≈ 1. At this pressure, the network chains of the two gels are in contact (see the density profiles in Figure 3), and the polymer-free solvent gap is ill-defined. Also, note that if the gels are nearly at contact, particularly at SB ≫ b/f, attractive electrostatic interactions can arise between the chains of opposite slabs.43 Condensation can G

DOI: 10.1021/acs.macromol.6b01416 Macromolecules XXXX, XXX, XXX−XXX

Macromolecules the counterion diffuse layers. Thus, the degree of localization of the counterions inside the hydrogels (i.e., extend of the counterion layers outside of the gels) can alter the pressure response significantly. There are two limiting cases for the pressure response; if all counterions are localized inside the gels, two gels interact only at physical contact. Consequently, the corresponding pressure drops to vanishing values at noncontact distances. Contrarily, if all counterions evacuate the gels to fill the gap, the surfaces can interact via their extended counterion layers at large intersurface separation distances δ ≫ δ*. Indeed, in Figure 4b we observe an osmotic pressure behavior in between these two limits at D/D0 > 0.5 due to the partially localized counterions inside the gels. Such behavior was also predicted for polyelectrolyte brushes by considering a parabolic form for the electrostatic potential.40 In this regard, on the basis of Figure 4a,b, one can conclude that the hydrogel surfaces considered here behave similar to polyelectrolyte brushes in the “weak osmotic regime”. This is the situation between the osmotic-brush regime, where the polyelectrolyte brushes behave like neutral brushes due to the localized counterions within the brushes, and nonosmotic brush, where condensation is weak (see Figure 4 of ref 40). Interestingly, the pressure drop in Figure 4b corresponds to an abrupt increase in the thickness of the gap at D/D0 > 0.5 as discussed above (in the Discussion section). Surprisingly, in Figure 6, where we show the rescaled pressure versus the rescaled thickness of the gap, the cancellation of these two effects leads to a smooth transition and perfect overlap of the data (i.e., δ2Pz ∼ constant at δ > δ*). Note that this steep decrease in the pressure (or in the counterion concentration in the gap) can be even more prominent for strong electrostatic couplings (i.e., SB > 1), for which ionic condensation onto chains is even stronger. At intermediate deformations (i.e., δ* > δ ≫ b), the osmotic coefficient is independent of counterion concentration (i.e., Pz/ kBTρgap ∼ constant) in high-dielectric media, similar to semidilute polyelectrolyte solutions. This indicates that ionic condensation with compression is balanced by the volume decrease upon compression. At strong compressions, however, condensation onto network chains can decrease the effective surface charge and lead to lower disjoining pressures. These effects can be incorporated by using cylindrical cell models suggested for polyelectrolyte solutions; however, applicability of those theories to polymer networks requires a detailed study. The relations presented in eqs 9 and 10 for the threshold pressure and the intersurface distance can be estimated for some relevant systems. If we consider a flexible backbone for the chains constituting the gels, a Kuhn size of b = 1 nm with SB ≈ 1 nm and 1kBT = 4.1 pN × nm give threshold pressure values of P*z ≈ 104 Pa and P*z ≈ 102 Pa for N ≈ 10 and N ≈ 100, respectively. Indeed, these pressures are normal loads measured in various soft tissues19 or in collapse of microgel particles.44 Above these values of pressure tensor, the gels enter a regime of deformation and squeeze out their solvent contents. For much stiffer polymeric systems (e.g., b ≈ 10 nm) and smaller N, due to decreasing entropic elasticity of chains, the threshold pressure also decreases. However, for stiff systems, the intrinsic mechanical response of the network itself should also be taken into account.

Article

■

CONCLUSION

■

AUTHOR INFORMATION

Coarse-grained molecular dynamics simulations and Poission− Boltzmann-based scaling analysis demonstrate that interactions between two polylelectrolyte gel surfaces are mainly dominated by repulsive forces. These forces are due to osmotic pressure of excess counterions occupying the volume between two surfaces. As two surfaces are brought closer, the osmotic pressure in the separation gap increases. At a critical pressure, the counterion osmotic pressure inside the gels and within the separation layer becomes equal, and the gels begin to deform under external pressure (Figure 4). The threshold pressure that separates weak and strong deformation regimes is lower for swollen gels (i.e., gels composed of high molecular weight chains) due to the decreasing effective surface charge density. The scaling prediction for the threshold pressure is consistent with simulations (Figure 6b). In the scaling limit, the pressure response is on the same order of magnitude as that of sparsely grafted polyelectrolyte brushes for surface separations smaller than Gouy−Chapman length. The distance between two gel surfaces is well described by the solution of the PB equation for a wide range of polymerization degrees at weak electrostatic interaction strengths, which are relevant to aqueous solvents. In the simulations, we also observe that the gel deformation starts before two gels are at physical contact, possibly due to the longrange nature of electrostatic interactions This indicates that for more swollen gels there is a wider range of pressure before two surfaces come to contact. This may be an important parameter in structures abundant in biological systems to reduce friction and wear. Our findings are also relevant to the studies of dense suspensions of microgels and their mechanical applications in anisotropic confinement, as well as crack formation in brittle hydrogels.

Corresponding Author

*E-mail [email protected] (M.O.d.l.C.) Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A.E. acknowledges Dr. Ozan S. Sariyer for his critical reading of the manuscript. This work was supported by the Center for Bio-Inspired Energy Science (CBES), which is an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award DE-SC0000989.

■

REFERENCES

(1) Guggenheim, E. A.; Donnan, F. G. Die genaue thermodynamik der membrangleichgewichte. Z. Phys. Chem. 1932, 162, 346−360. (2) Katchalsky, A.; Lifson, S.; Exsenberg, H. Equation of swelling for polyelectrolyte gels. J. Polym. Sci. 1951, 7 (5), 571−574. (3) Klein, J.; Kumacheva, E.; Mahalu, D.; Perahia, D.; Fetters, L. J Reduction of frictional forces between solid surfaces bearing polymer brushes. Nature 1994, 370 (6491), 634−636. (4) Gong, J. P.; Kagata, Go; Osada, Y. Friction of Gels. 4. Friction on Charged Gels. J. Phys. Chem. B 1999, 103, 6007−6014. (5) Sokoloff, J. B. Theory of the observed ultralow friction between sliding polyelectrolyte brushes. J. Chem. Phys. 2008, 129 (1), 014901. (6) Sokoloff, J. B. Theory of hydrostatic lubrication for two likecharge polymer hydrogel coated surfaces. Soft Matter 2010, 6 (16), 3856.

H

DOI: 10.1021/acs.macromol.6b01416 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules (7) Wynveen, A.; Likos, C. N. Interactions between planar polyelectrolyte brushes: effects of stiffness and salt. Soft Matter 2009, 6 (1), 163. (8) Pincus, P. Colloid stabilization with grafted polyelectrolytes. Macromolecules 1991, 24 (10), 2912−2919. (9) Zhulina, E. B.; Borisov, O. V. Structure and interaction of weakly charged polyelectrolyte brushes: Self-consistent field theory. J. Chem. Phys. 1997, 107, 5952. (10) Zhulina, E. B.; Borisov, O. V.; Pryamitsyn, V. A.; Birshtein, T. M. Coil-globule type transitions in polymers. 1. Collapse of layers of grafted polymer chains. Macromolecules 1991, 24 (1), 140−149. (11) Gong, J. P.; Kagata, G.; Iwasaki, Y.; Osada, Y. Surface friction of polymer gels: 1. Effect of interfacial interaction. Wear 2001, 251, 1183. (12) Mann, B. A; Holm, C.; Kremer, K. Swelling of polyelectrolyte networks. J. Chem. Phys. 2005, 122 (15), 154903. (13) Mann, B. A.; Everaers, R.; Holm, C.; Kremer, K. Scaling in polyelectrolyte networks. Europhys. Lett. 2004, 67, 786−792. (14) Erbaş, A.; Olvera de la Cruz, M. Energy Conversion in Polyelectrolyte Hydrogels. ACS Macro Lett. 2015, 4, 857−861. (15) Yin, D.-W.; Olvera de la Cruz, M.; de Pablo, J. J. Swelling and collapse of polyelectrolyte gels in equilibrium with monovalent and divalent electrolyte solutions. J. Chem. Phys. 2009, 131 (19), 194907. (16) Escobedo, F. A.; de Pablo, J. J Molecular simulation of polymeric networks and gels: phase behavior and swelling. Phys. Rep. 1999, 318 (3), 85−112. (17) Yin, D.-W.; Yan, Q.; de Pablo, J. J Molecular dynamics simulation of discontinuous volume phase transitions in highlycharged crosslinked polyelectrolyte networks with explicit counterions in good solvent. J. Chem. Phys. 2005, 123 (17), 174909. (18) Ou, Y.; Sokoloff, J. B.; Stevens, M. J. Comparison of the kinetic friction of planar neutral and polyelectrolyte polymer brushes using molecular dynamics simulations. Phys. Rev. E 2012, 85, 011801 DOI: 10.1103/PhysRevE.85.011801. (19) Dunn, A. C.; Cobb, J. A.; Kantzios, A. N.; Lee, S. J.; Sarntinoranont, M.; Tran-Son-Tay, R.; Sawyer, W. G. Friction Coefficient Measurement of Hydrogel Materials on Living Epithelial Cells. Tribol. Lett. 2008, 30, 13−19. (20) Gong, J. P. Friction and lubrication of hydrogels?its richness and complexity. Soft Matter 2006, 2 (7), 544. (21) Rennie, A. C.; Dickrell, P. L.; Sawyer, W. G. Friction coefficient of soft contact lenses: measurements and modeling. Tribol. Lett. 2005, 18, 499−509. (22) Heyes, D. M.; Brańka, A. C. Interactions between microgel particles. Soft Matter 2009, 5 (14), 2681. (23) Lyon, L. A.; Debord, J. D.; Debord, S. B.; Jones, C. D.; McGrath, J. G.; Serpe, M. J. Microgel Colloidal Crystals. J. Phys. Chem. B 2004, 108, 19099−19108. (24) Denton, A. R. Counterion penetration and effective electrostatic interactions in solutions of polyelectrolyte stars and microgels. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 67, 011804. (25) Gong, J. P. Why are double network hydrogels so tough? Soft Matter 2010, 6 (12), 2583. (26) Brown, H. R. A Model of the Fracture of Double Network Gels. Macromolecules 2007, 40, 3815−3818. (27) Chen, M.; Briscoe, W. H.; Armes, S. P.; Klein, J. Lubrication at physiological pressures by polyzwitterionic brushes. Science 2009, 323, 1698−1701. (28) Klein, J. Repair or Replacement A Joint Perspective. Science 2009, 323 (5910), 47−48. (29) Roba, M.; Duncan, E. G.; Hill, G. A.; Spencer, N. D.; Tosatti, S. G. P. Friction Measurements on Contact Lenses in Their Operating Environment. Tribol. Lett. 2011, 44, 387−397. (30) Meyer, A. E.; Baier, R. E.; Chen, H.; Chowhan, M. Tissue-onTissue Testing of Dry Eye Formulations for Reduction of Bioadhesion. J. Adhes. 2006, 82, 607−627. (31) Button, B.; Cai, L.-H.; Ehre, C.; Kesimer, M.; Hill, D. B.; Sheehan, J. K.; Boucher, R. C.; Rubinstein, M. A periciliary brush

promotes the lung health by separating the mucus layer from airway epithelia. Science 2012, 337, 937−941. (32) Liao, Q.; Dobrynin, A. V.; Rubinstein, M. Molecular Dynamics Simulations of Polyelectrolyte Solutions: Osmotic Coefficient and Counterion Condensation. Macromolecules 2003, 36, 3399−3410. (33) Kremer, K.; Grest, G. S. Dynamics of entangled linear polymer melts: A moleculardynamics simulation. J. Chem. Phys. 1990, 92 (8), 5057. (34) Grest, G. S.; Murat, M. Structure of grafted polymeric brushes in solvents of varying quality: a molecular dynamics study. Macromolecules 1993, 26 (12), 3108−3117. (35) Akerlof, G. Dielectric constants of some organic solvent-water mixtures at various temperatures. J. Am. Chem. Soc. 1932, 54 (11), 4125−4139. (36) Plimpton, S. Journal of Computational Physics. J. Comput. Phys. 1995, 117 (1), 1−19. (37) Manning, G. S. Limiting laws and counterion condensation in polyelectrolyte solutions I. Colligative properties. J. Chem. Phys. 1969, 51 (3), 924−933. (38) Schiessel, H.; Pincus, P. Counterion-Condensation-Induced Collapse of Highly Charged Polyelectrolytes. Macromolecules 1998, 31 (22), 7953−7959. (39) Sing, C. E.; Zwanikken, J. W.; Olvera de la Cruz, M. Effect of IonIon Correlations on Polyelectrolyte Gel Collapse and Reentrant Swelling. Macromolecules 2013, 46, 5053−5065. (40) Zhulina, E. B.; Boulakh, A. B.; Borisov, O. V. Repulsive Forces between Spherical Polyelectrolyte Brushes in Salt-Free Solution. Z. Phys. Chem. 2012, 226, 625−643. (41) Zhulina, E. B.; Rubinstein, M. Lubrication by Polyelectrolyte Brushes. Macromolecules 2014, 47, 5825−5838. (42) Fuoss, R. M.; Katchalsky, A. The potential of an infinite rod-like molecule and the distribution of the counter ions. PNAS 1951, 37, 579. (43) Naji, A.; Jungblut, S.; Moreira, A. G.; Netz, R. R. Electrostatic interactions in strongly coupled soft matter. Phys. A 2005, 352, 131− 170. (44) Lietor-Santos, J.-J.; Sierra-Martin, B.; Vavrin, R.; Hu, Z.; Gasser, U.; Fernandez-Nieves, A. Deswelling Microgel Particles Using Hydrostatic Pressure. Macromolecules 2009, 42, 6225−6230.

I

DOI: 10.1021/acs.macromol.6b01416 Macromolecules XXXX, XXX, XXX−XXX