We present simulations of the force-extension curves of strong polyelectrolytes with varying intrinsic stiffness as well as specifically treating hyaluronic acid, a polyelectrolyte of intermediate stiffness. Whereas fully flexible polyelectrolytes show a high-force regime where extension increases nearly logarithmically with force, we find that the addition of even a small amount of stiffness alters the short-range structure and removes this logarithmic elastic regime. This further confirms that the logarithmic regime is a consequence of the short-ranged “wrinkles” in the flexible chain. As the stiffness increases, the force-extension curves tend toward and reach the wormlike chain behavior. Using the screened Coulomb potential and a simple bead-spring model, the simulations are able to reproduce the hyaluronic acid experimental force-extension curves for salt concentrations ranging from 1 to 500 mM. Furthermore, the simulation data can be scaled to a universal curve like the experimental data. The scaling analysis is consistent with the interpretation that, in the low-salt limit, the hyaluronic acid chain stiffness scales with salt with an exponent of −0.7, rather than either of the two main theoretical predictions of −0.5 and −1. Furthermore, given the conditions of the simulation, we conclude that this exponent value is not due to counterion condensation effects, as had previously been suggested.

## I. INTRODUCTION

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.^{13} 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.^{14}

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*_{B}*T*, 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\pi 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.^{12} As in our previous work, we characterize the chain structure as a function of applied force.

## II. METHOD

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.^{20} 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.^{13} The second model specifically treats hyaluronic acid (HA) for comparison to new experimental data.^{14}

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/6}*d* 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 $\kappa =(8\pi lBcs)1/2$ is the inverse Debye length, *c*_{s} is the salt concentration, and the Bjerrum length *l*_{B} = *e*^{2}/*εk*_{B}*T* 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*ϵ*/*σ*^{2} and the maximum extent is given by *R*_{0} = 7*σ*. The contour length of the polyelectrolyte is *L*_{0} = *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,^{19}

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*_{b}*T*.

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.^{14} 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*_{B}*T*/*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*_{B}*T*.

## III. RESULTS

### A. General effect of stiffness

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*_{B}*T* ≳ 1) but deviate strongly at low *f* (*fℓ*_{B}/*k*_{B}*T* ≲ 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.^{4} Thus, for sufficiently large *k*_{a} (≃42*k*_{B}*T*), 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*π*/(2*L*_{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.^{22}

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} ≃ 42*k*_{B}*T*, the large *q* behavior matches that of a semiflexible chain.

### B. Hyaluronic acid model and salt dependence

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.^{14} 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,^{12} 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.^{14} 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}.^{14} 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*_{B}*T*/*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^{23} 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.^{7} 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,^{14} the value *α* = 0.7 was predicted theoretically in a model that accounted for salt-dependent counterion condensation.^{26} 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^{26} 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.

## IV. CONCLUSION

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.^{22} 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,^{14} 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}.

## ACKNOWLEDGMENTS

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.