The effect of chain stiffness and salt on the elastic response of a polyelectrolyte

Some of the seminal theoretical papers on polyelectrolytes treat the chain stiffness as a combination of intrinsic and electrostatic contributions.^1,2 Many polyelectrolytes fit into this categorization having significant intrinsic stiffness due to their molecular structure. One prominent example is double-stranded DNA (dsDNA), which has a large intrinsic persistence length. The force-extension curve for semiflexible dsDNA has been well characterized and there are simple formulas for the curve.^3,4 On the other hand, single-stranded DNA (ssDNA) and single-stranded RNA (ssRNA) are flexible, strongly-charged polyelectrolytes, and have yielded fundamentally different force-extension curves in single-molecule measurements.^5–10 In particular, experiments discovered logarithmic scaling of the extension with force for monovalent salt at high forces.^7–9 A fundamental issue in polyelectrolytes is how and to what degree intrinsic stiffness alters the chain conformational statistics and correspondingly the force-extension curves.

We have been able to characterize the fully flexible polyelectrolyte using a coarse-grained bead-spring model that treats the electrostatics, connectivity, and entropy.^11,12 This model reproduces the logarithmic dependence of the chain extension as a function of applied force in the single-molecule force-extension measurements and has given key insight into the chain structure as a function of applied force. The logarithmic elastic regime corresponds to the straightening of the short length scale “wrinkles” that occur in the local structure of the flexible polyelectrolyte due to neither electrostatic nor entropic interactions being dominant. More recently, the effect of the charge spacing on the force-extension curve has been studied in simulations.¹³ The transition at high forces from the logarithmic scaling for fully charged chains to the neutral scaling for sufficiently large charge separation was observed. The variation of the tensile screening length as a function of the applied force was explicitly exhibited in the single chain structure factor. Interestingly, it was also discovered that the excluded volume without electrostatic interactions can also produce the logarithmic dependence. In a similar vein, we are interested in examining the effect of chain stiffness on the force-extension curves and determining how stiffness alters the fully flexible force-extension curves. Recent experiments on hyaluronic acid (HA) which has a stiffness intermediate between flexible and semiflexible provide data for comparison.¹⁴ 

The theory of strong polyelectrolytes like HA and ssDNA has proven to be difficult because the electrostatic and entropic interactions are comparable, which eliminates perturbative methods. For strong polyelectrolytes, the Bjerrum length l_B, the separation at which the Coulomb energy is k_BT, is close to the charge separation distance in the polymer a. The screening of the electrostatic interactions by the salt brings in another length scale, the Debye length $lD=1/(8πlBcs)1/2$⁠, where c_s is the salt concentration. A broad theory literature has treated the interplay of electrostatics and chain stiffness in the flexible and semi-flexible limits, and has found two differing scaling relations for the electrostatic persistence length, either as l_D or $lD2$⁠.^1,2,15–19 Because the calculations are at or beyond their limit of application for strong polyelectrolytes, simulations that do not use additional approximations should provide missing insight into the structure of strong polyelectrolytes.

In this paper, we focus on simulations treating the effect of varying intrinsic stiffness for the difficult strong polyelectrolyte regime. In addition, the simulations will determine whether a simple bead-spring model is sufficient to describe hyaluronic acid force-extension data. In order to treat chains sufficiently long enough to have distinct elastic regimes, we use the screened Coulomb potential primarily at a relatively high salt concentration. This allows us to use Monte Carlo simulations with the pivot algorithm which is very efficient and can equilibrate very long chains in contrast to molecular dynamics simulations.¹² As in our previous work, we characterize the chain structure as a function of applied force.

The general problem we treat is a charged chain of N monomers with bond length b and a chain stiffness k_a being stretched by an applied force f. The force is applied in the z-direction and the chain extension L in the z-direction are calculated. The simulations treat N = 5000, which is sufficiently long that the persistence length and l_D are far shorter than the chain length. Consequently, on long length scales, the chain is a random walk and obeys the standard random walk (i.e., neutral) polymer statistics.

The system consists of a single polyelectrolyte chain in the simulation cell. Monte Carlo simulations were performed only using the pivot algorithm.²⁰ The polyelectrolyte chain is modeled with the potentials used in previous molecular dynamics (MD) simulations that used a standard bead-spring chain.^11,21 The same bond potential is used as in the MD simulations, but the bond length is fixed, since only the pivot move is used. Two models are treated. The first model is a generic polyelectrolyte with added stiffness that is based on previous work examining the variation of charge separation in the chain.¹³ The second model specifically treats hyaluronic acid (HA) for comparison to new experimental data.¹⁴ 

The excluded volume interaction between two monomer beads is given by the Lennard-Jones (LJ) potential,

where d is the bead diameter in units of σ, ϵ is the LJ unit of energy, and r is the distance between the two beads. We set the LJ cutoff as r_c = 2^1/6d for all interactions. The value of d is set to 4 Å for all particles. This creates a simple polyelectrolyte model with all ions of the same size. This value of d is a typical distance for the closest approach between two hydrated ions. In terms of interaction strengths at the closest ion-ion separations, 4 Å is a typical value. The focus of this work is on the scaling behavior and should not depend noticeably on the value of d.

We treat the electrostatic interactions using the screened Coulomb potential

where $κ=(8πlBcs)1/2$ is the inverse Debye length, c_s is the salt concentration, and the Bjerrum length l_B = e²/εk_BT is 7.1 Å. For these simulations, we use a cutoff such that U_sc(r_c) = 0.001ϵ.

As noted above, a fixed bond length of b = 0.96σ was used to match the average bond length in the previous MD simulations. The bond potential is the finitely extensible, non-linear elastic potential (FENE),

where the bond spring constant is k = 30ϵ/σ² and the maximum extent is given by R₀ = 7σ. The contour length of the polyelectrolyte is L₀ = b(N − 1), where N is the number of beads.

We introduce chain stiffness by using the standard three-body angle term in the potential that matches the stiffness term in the free energy of theoretical calculations,¹⁹ 

where the angle spring constant k_a is a variable and the θ is the angle between three consecutive beads. This angle potential has an equilibrium angle equal to 180°, which is special because it affects the tangent of the chain’s space curve. Thus, for sufficiently large k_a, the persistence length for a neutral chain with this angle potential is L_p = bk_a/k_bT.

In the generic polyelectrolyte model, the monomer is a single charged bead and using d = 4 Å, the charge separation corresponds to a = 3.67 Å and l_B = 1.86σ to match the previous setup. For this system, the spring constant for the angle potential is systematically increased from 0. The salt concentration is fixed at 200 mM, where the screened Coulomb approximation matched the all ion simulations. In a second system, which is setup to match hyaluronic acid, there are two beads per monomer, with one bead charged and the other neutral, because the charge spacing is larger in HA and, in fact, HA is a disaccharide with only one charged monomer. In this case, there are 10 000 total beads with every other one charged. The bond separation is b = 5.1 Å and l_B = 1.34σ. We chose k_a = 15 corresponding to an intrinsic persistence length of 8 nm, which is in the measured range.¹⁴ For the HA system, the salt concentration is varied from 1 to 500 mM. We note that at 1 mM, the Debye length is 25σ, which requires a long cutoff length, but it is small relative to the chain length and to the end-to-end distance at f = 0 which is 1406σ.

The force f is applied to both the terminal beads in opposite directions along the z-axis. We report the z component of the end-to-end distance as the extension L of the chain under applied force. In the plots, the forces are scaled by k_BT/l_B. Taking the temperature to be room temperature, the LJ unit of force ϵ/σ maps to 9.1 pN for system 1 and 3.2 pN for the HA system. However, comparisons are best made between dimensionless quantities such as fℓ_B/k_BT.

We first address the generic system with varying chain stiffness. The normalized force-extension plot for this system is given in Fig. 1. The solid points are for a fully flexible chain (k_a = 0), and the line at large f shows the logarithmic dependence found previously in experiments and simulations.^7–9,11,12 As k_a increases from 0, the extension naturally is larger at a given f, since the intrinsic or f = 0 conformation is more extended by the chain stiffness. In the large f regime, small deviations from the logarithmic scaling are noticeable even with k_a = 1. By k_a = 5, the deviations are large and clearly a different behavior occurs. However, for k_a ≤ 10, comparison with the Marko-Siggia formula for the wormlike chain shows noticeable differences especially at smaller forces. In particular, Fig. 1 shows the Marko-Siggia line using L_p = 12.5 for the k_a = 10 case. The Marko-Siggia values of L are accurate for large f (fℓ_B/k_BT ≳ 1) but deviate strongly at low f (fℓ_B/k_BT ≲ 0.1). For more stiff chains, k_a = 50, the Marko-Siggia formula fits well over the whole range of forces. We also include the more recent formula by Carrillo and Dobrynin, which is slightly better at intermediate forces.⁴ Thus, for sufficiently large k_a (≃42k_BT), the force-extension behavior corresponds to that of a semiflexible, wormlike chain.

We characterize the local structure of the chains through calculations of the structure factor S(q). The single chain structure factor is calculated using the formula

where r_j is the position of the j-th bead in the polymer chain and the normalization is S(0) = N. First, we examine S(q) at f = 0 for the varying k_a in Fig. 2 to see the structural variation as a function of k_a before tension is applied. In general, S(q) scales as q^−1/ν, where ν is the Flory exponent, and for these systems, we expect two scaling regimes. At low q corresponding to large separations much greater than L_p, the chains are self-avoiding random walks and have ν = 3/5, which is clearly seen in the figure. At large q, the intrinsic stiffness determines the structure. There is a change in the slope at log qσ = −1.2 for k_a = 50; this corresponds to a separation of 2π/q = 100σ ≈ 2 L_p, which is the Kuhn length. For q larger than this crossover point, the scaling switches to ν = 1 corresponding to a straight, stiff chain as expected due to inclusion of the angle potential. For the smaller values of k_a, the description of S(q) is similar. At low q, there is a ν = 3/5 scaling regime for all k_a. At high q, the value of ν decreases below 1 with decreasing k_a as shown in the inset of Fig. 2. The transition to the high q behavior occurs at higher q as k_a decreases.

At k_a = 1, the stiffness due to the bending potential is small and the electrostatic contribution to L_p is significant. In this case, the high q scaling is close to the fully flexible value of ν = 0.85 as shown in the inset. The logarithmic regime in the force-extension at k_a = 0 occurs because of the structure corresponding to this high q regime. While S(q) at k_a = 1 is close to that at k_a = 0 at high q, the difference is sufficient to yield a noticeable difference in the scaling of L. At stronger intrinsic stiffness (k_a > 1), the short range wrinkles in the fully flexible chain that yield the logarithmic regime are removed and the high q regime locally has ν = 1. Thus, the logarithmic regime does not occur.

The effect of the applied force on the chain structure is shown through the structure factor in Fig. 3 for k_a = 10 by plotting qS(q). The regions with S(q) ∼ 1/q are horizontal lines in the plot. Specifically at high q ≳ 2π/(2L_p) where the angle potential yields the persistence length, all the curves are flat. At small q where S(q) ∼ q^−5/3, the effect of the applied force is to progressively flatten the peak present at f = 0. The structural change corresponding to this peak flattening at small q is a combination of stretching the chain (increasing L) and progressively aligning more of the chain to the force direction at shorter and shorter length scales. Thus, the range of q for which S(q) is horizontal at low q increases with f. The minimum q of the range is decreased because the extension is increasing. The maximum q of the range is increasing because the alignment is progressing to shorter lengths as the Pincus tension length shrinks with f. The points in the figure at q^* are estimates of where the difference between S(q) at f = 0 and the particular value of f occurs. At a given f, the chain structure for q > q^* is unchanged with respect to the unstrained chain. This behavior corresponds to the tension screening as originally noted by Pincus.²² 

We now have further confirmation that the logarithmic scaling seen only with a fully flexible polyelectrolyte is due to the short-ranged, local structure of the chain. Our previous work showed that the logarithmic regime involves the stretching out of the chain at the short separations (or large q). In the chains with added stiffness, this local structure is altered. Thus, the force-extension plot is different at the large forces, as clearly seen in Fig. 1. While adding a small angle potential (k_a = 1) leads to a small difference compared to k_a = 0 in global structural quantities such as R_G, the local structure is sufficiently altered to yield a nonlogarithmic force-extension scaling. As k_a is increased, the local structural difference increases and the force-extension behavior tends toward the semiflexible form. By k_a ≃ 42k_BT, the large q behavior matches that of a semiflexible chain.

For the HA system, we set k_a = 15 to match the intrinsic persistence length of HA. From Sec. III A, for chains with every bead charged, we know this is an intermediate value of k_a that would yield behavior close to but not quite equal to that of the wormlike chain. For the HA system, we have calculated the force-extension curves for a range of salt concentrations (c_s = 1, 3, 10, 50, 100, 200, and 500 mM) to compare with recent experiments.¹⁴ We performed the simulations at the lower c_s to see how well the screened Coulomb model would reproduce the scaling observed in the experiments. While previous comparison to MD simulations treating all the ions explicitly for a ssDNA model determined that below 100 mM the screened Coulomb approximation has noticeable differences,¹² these differences will be smaller in HA. HA is not as strongly charged as ssDNA, which will make the screened Coulomb approximation better for HA. In addition, the charge separation in HA is larger than the Bjerrum length, which eliminates counterion condensation effects that are important in ssDNA.

To validate that the simulation captures the salt-dependent elasticity of HA, we compared the simulated force-extension curves, at three different salts (1, 10, and 500 mM), to curves measured in recent experiments.¹⁴ Comparison was carried out by finding multiplicative factors in both the force and length that would best transform the experimental data so as to overlay the simulation data; the best transformation was judged by minimizing the sum of least-squared differences between data and simulations. The length factor is needed to account for contour length differences between the experimental and simulated chains; however, since the extension is expected to be simply proportional to the contour length, this factor is not very meaningful. The force factor is more important: a factor near unity would indicate that the simulation is properly representing the real, physical forces (electrostatic and tensile) present in the experiment, while a factor far from unity would indicate an error in the microscopic representation.

We find that the best-fit force factors are indeed near unity when comparing experiments and simulations of equivalent ionic strength. Particularly, data acquired on HA at 10.4 mM ionic strength are well-fit to the simulation performed at 10 mM, using a force factor of 1.05 (green curve, Fig. 4); similarly, the 500 mM data match the 500 mM simulation using a factor of 1.14 (red curve). Data acquired at 1.44 mM ionic strength fit to the 1 mM simulation with a force factor of 1.42 (blue curve, Fig. 4); the slightly larger factor here is likely due to the slight mismatch of ionic strength. However, this difference is minor relative to comparisons involving strongly divergent ionic strengths. In such cases, the data could still be partially overlaid onto the simulation, but the force factor was far from unity. For example, the 1.44 mM data overlay the 500 mM simulation only when using a force factor of 0.16; vice versa, the 500 mM data can be made to overlay the 1 mM simulation with a force factor of 8.0. These results confirm the validity of the simulation in representing the electrostatic and conformational behavior of HA and allow us to interpret HA behavior in more detail through analysis of the simulation.

The salt dependence of the force-extension curves in the simulation is shown in Fig. 5(a). The extension decreases as c_s increases due to the increased screening. At the largest f, the curves merge, where the chain is highly stretched and the applied force is the dominant interaction. The plot shows that the Marko-Siggia curve fits well at large f and high salt concentrations. At large c_s and low f, the Marko-Siggia curve is slightly below the simulation data. This is consistent with the above result for k_a = 10 (and c_s = 200 mM) in Fig. 1 that found larger k_a values that were needed to fully reach the wormlike chain regime at low f. However, at low salt, the Marko-Siggia curve can only fit the high force regime and is progressively worse as the salt concentration decreases.

The experimental data on HA have been shown to have a master curve by rescaling each curve by salt dependent parameters f_c and L_c.¹⁴ We apply the same analysis to the simulation data by fitting each data set to the 100 mM data and plot the scaled data in Fig. 5(b). The overlap is very good, reproducing another characteristic of the experiments and demonstrating that the scaling does occur within the screened Coulomb approximation. Furthermore, even though HA is a somewhat complex polymer (local secondary structure; charges slightly off of the backbone), the simple simulated bead-spring model reproduces the experimental results. This validates the use of HA as a model polyelectrolyte with intermediate stiffness, where chemical details are unimportant for the questions we are exploring.

One of the interesting results of the experiments is the dependence of f_c on the ionic strength. The expectation is that f_c is proportional to the crossover force where the Pincus regime ends which should scale as k_BT/L_p. Thus, one can obtain the scaling behavior of L_p from f_c. The experimental data were fit as a function of ionic strength I to

where I_ref is the reference I (100 mM) used in the fitting, I_c defines the plateau ionic strength, and α is the scaling exponent for the electrostatic persistence length, L_e(I) ∼ I^−α. Least squares fitting to the simulation data with α fixed at 1/2 gives I_c = 12.5 mM and relaxing the constraint gives α = 0.71 and I_c = 9.5 mM. In Fig. 6, we plot f_c as a function I with the fit curves. This result is similar to the experimental result obtained with α = 0.65 as a best fit with I_c = 21 mM. The experimental data reach lower I where the difference between the predictions of the Odijk-Skolnick-Fixman (OSF) model^1,2 as well as others^15,18 and the Barrat-Joanny (BJ) model²³ as well as others^16,24,25 is largest. The simulation data are not as strongly distinguishing as the experimental data. As in the experimental data, we find the salt dependence on the simulation values for the L_c to be weak.

The value of α for the simulations is close to the experimental value and both are clearly larger than α = 0.50. Prior studies on flexible polyelectrolytes found α = 0.5.⁷ This result corresponds to a linear scaling with the Debye length and matches many theoretical predictions for flexible polyelectrolytes.^16,19,23,25 The new result for HA, a system with intermediate stiffness, being distinct from the flexible case raises some questions. As noted previously,¹⁴ the value α = 0.7 was predicted theoretically in a model that accounted for salt-dependent counterion condensation.²⁶ However, the present results rule out this mechanism, since the screened-Coulomb simulations do not contain counterion condensation yet reproduce the experiments. Furthermore, the HA charge density is below the Manning limit, 1/l_B, so strong counterion condensation is not expected. We conclude that the agreement of α with the prior theory²⁶ is a coincidence. But, the correct explanation is not obvious. We suggest that the fact that experiments and simulations match so well implies that the simple bead-spring model is a good starting point for theory. Indeed, the underlying simulation Hamiltonian is the one that many calculations start from, but, as noted in the Introduction, does present challenges especially in the strong polyelectrolyte regime; some approximations may not be valid. For the moment, we only conclude that the present result implies that there is something missing in the theories.

We have included intrinsic stiffness through a bending potential with spring constant k_a in simulations of the force-extension of a single polyelectrolyte chain. Even for small values of k_a, the logarithmic regime in the force-extension curve that occurs at large applied forces for fully flexible chains disappears. This result further confirms the picture that the logarithmic regime is a consequence of the short-ranged structure, which is only stretched at large applied forces. In the systems with an intrinsic stiffness, it is this short-ranged structure that is altered resulting in the loss of the logarithmic regime and in a different functional dependence of the force-extension curves. As k_a increases, the functional form of the force-extension curves approaches the wormlike chain behavior as represented by the Marko-Siggia formula. At k_a = 50, the force-extension data are well represented by the formula. In calculations of the single chain structure factor, we explicitly see the tension screening behavior predicted by Pincus.²² That is, at a given applied force, there is a q^* where the structure factor is identical to that at zero force (and thus screened) for q > q^*. Unfortunately, determining the dependence of q^* on the applied force f requires simulations of a higher fidelity than in the present set.

Simulations specifically modeling hyaluronic acid have also been performed to compare to recent experiments studying the salt dependence. Hyaluronic acid has an intermediate stiffness between a fully flexible and a wormlike chain. Using k_a = 15 in the simulations, the overlap for individual salt concentrations of the force-extension data between the simulation and experiment is very good. Thus, the simple model incorporating bonds, an angle term, and the screened Coulomb interactions works well for hyaluronic acid. Also matching the experiments,¹⁴ the simulation data for varying salt concentration can be scaled to a universal curve. One interesting but unexplained result is that the scaling exponent for the salt concentration dependence is 0.71, which is close to the experimental value 0.65 ± 0.02. In neither case is the value consistent with one of the two commonly predicted values, 0.5 or 1.0. This raises the question of what happens for larger k_a.

O.A.S. acknowledges support from NSF Grant No. DMR-1611497. J.P.B. acknowledges support from NSF MRSEC Grant No. DMR-1428302. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under Contract No. DE-NA0003525. This work was performed at the U.S. Department of Energy, Center for Integrated Nanotechnologies, at Los Alamos National Laboratory (Contract No. DE-AC52-06NA25396) and Sandia National Laboratories. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.