Article pubs.acs.org/Macromolecules
Simulations of Stretching a Strong, Flexible Polyelectrolyte Mark J. Stevens* Sandia National Laboratories, P.O. Box 5800, MS 1315, Albuquerque, New Mexico 87185-1315, United States
Dustin B. McIntosh Physics Department, University of California, Santa Barbara, Santa Barbara, California 93106, United States
Omar A. Saleh Materials Department and BMSE Program, University of California, Santa Barbara, Santa Barbara, California 93106, United States ABSTRACT: We present simulations of the stretching of a single strongly charged, flexible polyelectrolyte chain. Over a range of salt concentrations and valences, our results quantitatively match recent single-molecule measurements of the elastic response of single-stranded DNA. Of particular interest is a heretofore unexplained high-force regime of logarithmic elasticity seen in experiments with monovalent salts. We investigate this by calculating the forcedependent structure factor of the charged chain. We find distinct structures at long and short length scales and show that the compliance of the short-scale crumpled structure underlies the logarithmic elasticity regime. This is corroborated by results in divalent salts, in which stronger electrostatics leads to more short-scale crumpling and a corresponding increase in the high-force compliance. Thus, we show that the high-force elasticity of a charged chain is a signature of a structural regime that exists in flexible polyelectrolytes, but not neutral polymers.
trolytes amenable to single-molecule testing.11−15 Dessinges et al.11 showed that the elastic behavior of chemically denatured ssDNA in monovalent salt is quite distinct from that predicted by simple polymer models and were the first to note that ssDNA displays an extensive logarithmic elastic regime. Seol et al.12 measured the force−extension curve for poly(U) ssRNA, which has negligible base pairing and stacking interaction, finding that the extension also seems to increase approximately logarithmically with force. Saleh et al. measured the force− extension behavior of chemically denatured ssDNA across a wide range of salt conditions, including various concentrations of both monovalent and divalent cations.13−15 They found a salt-dependent transition from a low-force regime that matches scaling theory for long chains to a high-force, valence-specific behavior. Specifically, over 2 decades of monovalent salt, the extension again scales as a logarithm at high force, while the divalent data display deviations from the logarithmic behavior
I. INTRODUCTION A full understanding of the structure of flexible polyelectrolytes is a major outstanding problem in the theory of macromolecules, with importance both to technological and medical uses of charged polymers and to understanding of charged biopolymer (e.g., RNA, protein) behavior. Accordingly, it has been an active area of research over the past two decades.1−8 The theory of flexible polyelectrolytes is particularly difficult in the regime where the interaction strengths of the entropic and electrostatic interactions are comparable, and mean-field methods break down. There has been significant progress, but a comprehensive theory for the physics of flexible polyelectrolytes is yet to be established.5,9,10 In this context, the generation of direct and quantitative data on flexible polyelectrolytes plays a critical role in advancing the field. An important development of the past decade is the development of single molecule experiments that measure polyelectrolyte force/extension relations, as these can be directly connected to theory. In particular, single-stranded nucleic acids (ssDNA and ssRNA) that lack secondary structure have been established as flexible, strongly charged polyelec© 2012 American Chemical Society
Received: May 3, 2012 Revised: June 27, 2012 Published: July 12, 2012 5757
dx.doi.org/10.1021/ma300899x | Macromolecules 2012, 45, 5757−5765
Macromolecules
Article
dominates locally, causing flexible polymers to coil up. In this case the applied force orients the (neutral) chain into a straight set of “tension blobs”.22 Inside the tension blobs there is a random walk structure. The standard scaling theory then predicts
that the authors speculated may be due to local ion-wrapping by the polymer. The logarithmic behavior, observed by multiple experiments on ssDNA and ssRNA, over a large monovalent salt concentration range, has yet to be theoretically derived. To better understand these issues, we have performed molecular dynamics (MD) simulations of the stretching of ssDNA using a bead−spring model for flexible polyelectrolytes. With this model, we treat the competition between the entropy, electrostatics, and the stretching force. We study both monovalent and divalent salt systems, as in the experiments. We show that this model can reproduce key aspects of the measurements, including the valence-specific universal response and the logarithmic elastic regime. Additionally, the simulations permit direct examination of the structure of the stretched polyelectrolyte, which we quantify by examining changes in the structure factor as force is applied. Further, the ions are explicitly modeled, enabling direct examination of ioncoordinated structure. In particular, we find that the polymer forms short-length-scale, counterion-associated bent structures that we term “wrinkles” and that the force-induced straightening of these structures dominates the high-force elastic response, including the logarithmic regime seen in monovalent salt. Our data thus permit a microscopic-level understanding of the physics underlying polyelectrolyte elastic response.
R ∼ fγ
(1)
where γ = 1/ν − 1 with ν the Flory exponent. In a good solvent, ν ≈ 3/5 and γ ≈ 2/3. Analysis of polyelectrolyte behavior is complicated by the introduction of other length scales including the Debye length SD and the electrostatic blob size ξe, the length scale over which unscreened intrapolymer electrostatic repulsion becomes of order the thermal energy.23 At typical salt concentrations, experimental polymers can be much larger than SD (e.g., in 100 mM monovalent salt, SD = 9.6 Å) and, on long length scales, random walk scaling is recovered. The standard theory then applies to the experiments in the low force regime. However, the simple scaling prediction can be expected to break down when multiples of the relevant length scales are comparable. The simulation data presented here are in the regime kBT/f ∼ SD ∼ ξe ∼ l, consistent with experimental data on single-stranded nucleic acids. In this complex regime, theoretical predictions are lacking; a logarithmic force law is not unexpected13 but has not been derived. With the ability of the simulations to examine the chain structure and the ionic interactions simultaneously, insight into these issues may be gained.
II. REVIEW OF POLYELECTROLYTE ELASTICITY THEORY Most of the theoretical analysis of single polymer stretching has focused on semiflexible polymers, particularly dsDNA. dsDNA’s elastic behavior is well-described by the Marko− Siggia worm-like chain model.16 In this model, electrostatics can largely be ignored due to dsDNA’s large intrinsic persistence length. An exception is at very low salt concentrations, where short-range electrostatic repulsion stiffens the chain. Marko and Siggia accounted for this effect at the mean-field level by including the electrostatic kernel of Barrat and Joanny.16,17 Similarly, most of the simulations have focused on semiflexible polymers; of particular relevance is recent work that examined the effect of electrostatic interactions on semiflexible polymer elasticity.18 For flexible polymers, the results of the semiflexible theories are limited, although their methods could be applicable with suitable adjustment. The ssRNA poly(U) force−extension data12 was partially reconciled with the Marko−Siggia/Barrat− Joanny approach,16,17 but the logarithmic dependence in the data persists to lower forces than the theory predicts. Recently, Dobrynin et al. have proposed a unified chain deformation model for flexible, neutral chains,19 in analogy with previous studies of semiflexible polymers.20,21 This model is analytic for all forces and successfully fits the ssDNA force−extension behavior at high salt, where electrostatics are largely screened away. However, the model fails, expectedly, at lower salt concentrations where electrostatics become more relevant. To our knowledge, a comparable theoretical analysis that includes electrostatic effects on flexible polyelectrolytes has not been performed. Monte Carlo simulations were presented in the early experimental paper on ssDNA.11 They used a screened Coulomb model, which qualitatively produces the same shape as the experimental force−extension curves, but differs quantitatively. Another factor for flexible polymers is the effect of longrange excluded volume interactions which occur for forces f ≪ kBT /l, where l is the Kuhn length; at such low forces, entropy
III. METHOD The system consists of a single polyelectrolyte chain in the simulation cell with counterions and salt ions. The chain has 200 beads, each of which has a single charge of negative sign. The counterions are monovalent. The valence of the positively charged salt ion is varied as in the experiments. We treat 1:1 and 2:1 added salt, which we refer to as monovalent and divalent (salt) systems. The polyelectrolytes chains were modeled using a standard bead−spring model.2,24 We used a Langevin thermostat with temperature 1.2ε with damping constant of 1.0τ−1, where τ is the LJ time unit. The time step is 0.010τ. The force f is applied to both terminal beads in opposite directions along the z-axis. We report the z component of the end-to-end distance R as the extension Rz of the chain under applied force. The simulation cell is longer in the stretch direction than the lateral two dimensions in order to encompass a fully stretched chain. We chose Lz = 300σ, which is much longer than the contour length of the chain. Most of the simulations were performed with the lateral dimensions being Lx = Ly = 84.85σ. The lateral dimension must at least be greater than the equilibrium chain size and sufficiently large that the ion distribution about the chain can decay to the bulk value. In some of the simulations at the lowest divalent salt concentrations (and low forces), this system width was too small, and the lateral dimension was increased by a factor of 1.5−3. The system width was determined to be sufficiently large when the value of Rz became constant with increasing system width. The excluded volume interactions between monomer beads is given by the Lennard-Jones (LJ) potential: 5758
dx.doi.org/10.1021/ma300899x | Macromolecules 2012, 45, 5757−5765
Macromolecules
Article
⎧ ⎡⎛ ⎞12 ⎛ ⎞6 d d 1⎤ ⎪ ⎪ 4ε⎢⎜ ⎟ − ⎜ ⎟ + ⎥ , r < rc ⎝r⎠ 4 ⎥⎦ ULJ(r ) = ⎨ ⎢⎣⎝ r ⎠ ⎪ ⎪ 0, r > rc ⎩
Usc = kBT ∑ zizj i,j
where κ = (4πlBcs(1 + z+)) is the inverse Debye length, cs is the salt concentration, and z+ is the valence of the positive salt ion. We only do the screened Coulomb simulations for the 1:1 salt as we do not expect the screened Coulomb to be reasonable for divalent ions. We present data here using for monovalent salt concentrations of 100 and 200 mM. For these simulations we use a cutoff rc, where Usc(rc) = 0.001ε.
(2)
IV. RESULTS A. Force−Extension Curves. We calculated the extension Rz for a range of forces f for 20−200 mM monovalent and 0.2− 10 mM divalent salt concentrations. Notably, surprisingly small concentrations of divalent ions are required to produce similar effects on the elastic behavior of the polymer in comparison to the large monovalent concentrations, an effect also noted in the experimental data.15 The range of f starts at small values where the energy of the pulling force is less than the thermal energy, i.e., f lB/kBT ≪ 1. At the high end of the force range, the chain is becoming taut, but only at the largest force do the bonds begin to stretch. We are not interested in treating (with a coarsegrained model) bond stretching and so did not enter that regime. The plot of the force−extension data for the monovalent and divalent systems in Figure 1 shows the
(3)
where the spring constant is k = 30ε/σ2 and the maximum extent is given by R0 = 7σ. The average bond length, which is also the average distance between charges, is a = 0.96σ. Taking the ssDNA monomer separation to be 6.4 Å gives σ = 6.67 Å. The contour length of the polyelectrolyte is L0 = a(N − 1) = 192.0σ. The LJ unit of force ε/σ maps to 5.2 pN, taking the temperature to be room temperature. In the plots the forces are scaled by kBT/lB = 1.13ε/σ, which maps to 5.8 pN at room temperature. All the charges interact through a medium of dielectric constant ε that fills all space. Equivalently, for water the Bjerrum length, lB = e2/εkBT is 7.1 Å. In this medium, any two charged particles i and j interact via the Coulomb potential, and the total Coulombic energy is given by UCoulomb = kBT ∑ zizj i,j
lB rij
(5)
1/2
where d is the bead diameter in units of σ, ε is the LJ unit of energy, and r is the distance between two beads. We set the LJ cutoff as rc = 21/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 closest approach between two 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. In fact, the experimental results did not depend strongly on ion types which have different sizes.15 The bonds of the chain are modeled using a finitely extensible, nonlinear elastic potential (FENE) Ubond(r ) = −1/2kR 0 2 ln(1 − r 2/R 0 2)
lB exp( −κrij) rij
(4)
where z is the valence, kB is Boltzmann’s constant, T is temperature, and rij is the separation between particles i and j. The sum of all the Coulombic interactions is calculated using the PPPM method, described elsewhere.25 Periodic boundary conditions in all directions are used. Modeling ions explicitly introduces several computational issues. First, particularly in monovalent salt, the number of particles in the simulation cell is large, requiring heavy computational resources. We have thus been limited to four monovalent salt concentrations. The divalent salt concentrations are much smaller, and as a result, we obtain data at six divalent salt concentrations. However, in 0.2 and 0.5 mM divalent salt, the Debye length is large (25.3σ and 16.6σ, respectively). We found the distance for the ion distribution about the polymer to decay was large in these cases requiring an increase in the simulation cell to obtain accurate data at low forces as discussed above. In all cases, fluctuations in the extension are large at low forces and computational times required to reduce the error can be excessive. As a result, the error in the extension at some of the low force data presented here is substantial, and we are unable to do simulations at lower forces. We also performed some simulations using the screened Coulomb potential
Figure 1. Plot of the extension Rz normalized by the contour length L0 versus the applied force f normalized by kBT/lB for different salt concentrations cs. (a) Monovalent salt at cs = 20 mM (magenta triangle), 50 mM (blue square) 100 mM (red circle), and 200 mM (green star). (b) Divalent salt at cs = 0.2 mM (purple up triangle), 0.5 mM (brown star), 1 mM (magenta square), 2 mM (red circle), 5 mM (green diamond), and 10 mM (blue down triangle).
extension normalized by the contour length L0 and the pulling force normalized by kBT/lB. In the weak force regime the uncertainty in the simulation data is larger, since the fluctuations in Rz are large relative to its magnitude. The error bars are plotted for the low force cases when the uncertainty is greater than the point size. At a given force, the extension decreases with increasing salt concentration as predicted due to electrostatic screening of the 5759
dx.doi.org/10.1021/ma300899x | Macromolecules 2012, 45, 5757−5765
Macromolecules
Article
Table 1. Salt Concentration cs, Scaling Parameters fc and Rz(fc), and Structure Factor Exponents at f = 0 for Monovalent (top) and Divalent (bottom) Casesa
monomer−monomer repulsion of the charged backbone. At the lowest forces the extension depends on the equilibrium (f = 0) extension, which decreases with increasing salt concentration. As the force increases, the different data sets merge and then overlap at the highest forces. In this regime the pulling force is dominating over entropy and electrostatics. The largest forces yield an extension of almost 180σ which is close, but less than the contour length L0 = 192σ. We did not go to larger f because the bonds are beginning to be stretched at the highest force (f = 12.8ε/σ). One of the major results of the experiments is that suitable scaling results in overlap of the data.13,15 We have employed similar procedures to produce scaled plots of the force− extension curves in Figure 2. Specifically, we extract a salt-
cs 20 50 100 200 0.5 1 2 5 10 a
fc
Rz( fc)
νl
νh
0.077 0.122 0.203 0.360
44.3 47.4 52.2 62.1
0.82 0.73 0.67 0.64
0.79 0.78 0.76 0.75
0.050 0.055 0.080 0.150 0.250
28.4 27.9 30.2 36.6 41.3
0.83 0.78 0.76 0.71 0.67
0.70 0.69 0.67 0.66 0.65
The salt concentration is in mM, and the fc and Rz( fc) are in LJ units.
mM, which has the same values as 0.5 mM in the scaling plot (Figure 2). However, since the force−extension data for 0.2 mM do not appear to reach forces below fc, the uncertainty is too large, and we have excluded the point in the figure and the table. In the experiments, the divalent salt is in competition with an 8 mM monovalent buffer which the authors originally associated with the plateau of fc at low concentrations. The simulations indicate that there may be another cause of this effect, but presently we have not been able to ascertain the source of this behavior. In the high salt regime, the power-law exponent of the fc line in Figure 3 is 0.66, and the exponent for
Figure 2. Scaled plots of the extension Rz versus the scaled applied force f for different (top) monovalent and (bottom) divalent salt concentrations cs. Simulation data points are the same as in Figure 1. In addition, there are black solid points from the experimental data for NaCl at 20 mM (triangles) and for CaCl2 at 0.2 (triangles) and 5 mM (squares).
dependent crossover point ( fc, Rz( fc)) for each force−extension curve. For the monovalent salt curves this is done by a fit to the phenomenological function from ref 13, fixing γ = 0.6.26 For the divalent data, the crossover point was extracted by hand, overlapping the data by eye. Figure 2 shows not only that the simulation data scales but also that simulation data can be overlapped with the experimental data. The overlap between the simulation and experimental data for both monovalent and divalent salt is excellent.27 The values of fc and Rz( fc) are given in Table 1 and plotted in Figure 2. For the monovalent data, we find that fc scales as a power law in salt concentration cs, fc ∼ cαs , with α = 0.54, which agrees with the experimental result of α = 0.54 for the NaCl system and α = 0.47 for the KCl system. We find that Rz( fc) = Lc also scales as a power law in cs with exponent 0.09, which also agrees with the experimental values of 0.11 for NaCl and 0.088 for KCl. We note a log−log plot of Lc versus the Debye length indicates that the monovalent and divalent data lie on a single line. For the divalent systems, both the simulation and experimental data indicate that the low salt regime is distinct, close to being constant as a function of cs.13,15 This leveling off at low cs would be more clear if we included the point at 0.2
Figure 3. Plot of the scaling values fc and Lc as a function of salt concentration for the monovalent (squares) and divalent (circles) systems. The solid points are for the all ion simulations, and the open points are for the screened Coulomb simulations. The numbers indicate the slopes of the fit lines.
the Lc line is 0.18. These values are larger than the monovalent ones, which is in agreement with the experiment, but the magnitudes appear to be smaller than the experimental ones. One of the intriguing results of the experiments is that the large f regime has a logarithmic dependence in monovalent salt. A semilog plot (Figure 4) of the scaled monovalent data from the simulations clearly shows a logarithmic dependence. The straight line is a least-squares fit to the data for log f/f > 0.4. Varying the range of the fit data set slightly shifts the line but 5760
dx.doi.org/10.1021/ma300899x | Macromolecules 2012, 45, 5757−5765
Macromolecules
Article
strong applied force, the chain conformation is almost a linear, taut chain. The extension is 85% and 93% of contour length for the monovalent and divalent cases, respectively. There remain some amount of short length scale bends (“wrinkles”) at this force. The chain is substantially stretched, but as stated earlier, the bonds are not being strained at this force. At the intermediate force, the extension is a little less than half the extension at the very high force. There are large fluctuations in the conformation lateral to the applied force direction. In the weak force case, the conformations are only stretched slightly compared to the zero force conformations. To examine the structure at different length scales, we calculate the single chain structure factor N
Figure 4. Semilog plot of the scaled extension Rz versus the scaled applied force f, demonstrating a logarithmic elasticity regime for the monovalent data (same point types as in Figure 1). A single divalent case at 2 mM (black squares) shows the difference between the valences. The line is a least-squares fit to the data.
S( k ) =
1 |∑ exp(ik·rj)|2 N j=1
(6)
where rj is the position of the jth bead in the polymer chain and the normalization is S(0) = N. Our calculations are for the time average and the spherical average of k. For equilibrium polymers S(k) ∼ k1/ν. A good solvent neutral polymer has ν = 0.588, and we can obtain the three digit accuracy for the neutral systems.28 For sufficiently long polyelectrolytes such as the ones in the experiments the electrostatic interactions are fully screened at long length scales and the chain is a random walk. In this case, there will be at least two distinct regimes with different values of ν. The random walk structure at large length scales would yield neutral scaling at small k. The stretched polyelectrolyte structure at shorter length scales which is not (fully) screened would yield a larger ν at large k. For a polyelectrolyte under tension, the scaling exponent will change as a function of the applied force. As the force increases ν will tend toward 1, which is the value for a straight linear conformation. For a fully taut chain, the ν = 1 scaling would hold on the range 2π/L0 ≲ k ≲ 2π/a. We show an example of the effect of the tensile force on S(k) for the 1 mM divalent system in Figure 6. The range of k displayed goes from the k → 0 limit where S(k) = N to near 2π/a, where the rise of the bond peak is visible. For reference log(2πσ/L0) = −1.5. The dashed lines are fits to different sets of data indicating the range of ν as a function of f. As a baseline, we show the line for ν = 0.59 that comes from a fit to the corresponding neutral system. The line with νh = 0.68 is a fit to the high k regime (−0.3 < log kσ < 0.3) for the f = 0.0ε/σ case. The figure shows that f = 0.4ε/σ has the same S(k) in this high k regime as f = 0. In the force range between these values, the changes occur only in the low k regime. The line with νl = 0.98 is a fit to the low k regime (−0.9 < log kσ < 0.4) for the f = 6.4ε/σ case. We see that the chain does become rodlike (i.e., ν → 1) at this high force. In addition, we see that at large enough forces ( f ≳ 0.4ε/σ) the short-range structure determined by high k is altered.
has no effect on the conclusion that the data have a logarithmic regime. The simulation data are predominately in the logarithmic regime. Only the lowest f values show a deviation above the least-squares fit line. Thus, we have shown that this simple bead−spring polyelectrolyte model has recovered a number of results from force−extension experiments on ssDNA: (1) A vastly smaller concentration of divalent ions, compared with monovalent ions, is required to modulate the elastic behavior of the polymer. (2) The force−extension data can be rescaled by a salt-dependent crossover (fc, Rz(fc)) which collapses the data to a single curve overlapping with the experiments. Notably, this curve indicates that Rz scales logarithmically with force for a broad range of forces and monovalent salt concentrations. (3) The crossover forces fc and extensions Rz( fc) scale with salt concentration as power laws that compare favorably with those measured in experiments. These results indicate that the bead−spring model incorporated sufficient features to capture the experimental data. The molecular details of the monomer are not important for this scaling result. The main interactions are entropy, electrostatics, and the extension force. Any other interactions (e.g., chemical) are negligible. B. Force-Dependent Structure. In addition to the end-toend extension, the simulations give us access to the full structure of the polyelectrolyte as well as the explicitly modeled ions as force is applied. Images of the polymer conformations along with (condensed) nearby ions are shown in Figure 5. Because ssDNA is a strong polyelectrolyte, there are many condensed ions. In the divalent systems, almost all the condensed ions are divalent, outcompeting the monovalent counterions, as expected. The conformations are shown at very weak, intermediate, and very strong applied forces. At the very
Figure 5. Images of the stretched polyelectrolyte for the monovalent salt case at 50 mM for (a) f = 6.4ε/σ (scaled 50%), (b) f = 0.40ε/σ, and (c) f = 0.05ε/σ and the divalent salt case at 2 mM for (d) f = 12.8ε/σ (scaled 50%), (e) f = 0.40ε/σ, and (f) f = 0.05ε/σ. The monomers are blue. All ions within 2σ are shown. The positive salt ions are red; the equivalent counterions (if any) are green; negative salt ions are cyan. 5761
dx.doi.org/10.1021/ma300899x | Macromolecules 2012, 45, 5757−5765
Macromolecules
Article
Figure 6. Single chain structure factor for the divalent system at 1 mM salt concentration is shown for f = 0.0, 0.4 1.6, and 6.4ε/σ (right to left). The dashed lines are fits that give the value of the exponent ν; each is labeled with the best-fit value of ν. The inset shows the f = 0 case with the low and high k fits (dotted and dashed lines, respectively).
Figure 7. Scaling exponents from S(k) for the monovalent systems (top) and divalent systems (bottom). The point types indicate the salt concentrations and are the same as in Figure 1. Each S(k) has been fit as described in the text over two regimes (high and low k), yielding two ν. The points connected by the solid line designate νh for the high k regime, and the dotted line designates νl for the low k regime.
In general, we find that a single power law does not fit S(k). This contrasts with the neutral polymer case, which can be fit with one exponent over the above ranges (data not shown). Two lines representing two regimes with distinct values of ν are required even at f = 0 as is shown in the inset of Figure 6, which shows the low and high k fits for the f = 0 case. This is an intrinsic feature of a polyelectrolyte chain. The equilibrium system at 1 mM has νl = 0.68 and νh = 0.78. In general, the transition point depends little on salt concentration and valence. The transition occurs at about k = 0.5/σ, which corresponds to a length r = 2π/k = 12.6σ. The high k fit is good for all cases at least for k > 0.5/σ, implying that this regime occurs for lengths r up to at least 12σ, which is an appreciable length. The need for (at least) two length scales to describe polyelectrolyte structure is not new. In calculations of the bond−bond correlation function for the screened Coulomb simulations of polyelectrolytes, Dobrynin et al. found it necessary to use two length scales to obtain good fits.29,30 The bond−bond correlation function and the structure factor are exhibiting the same length-scale-dependent structural variation. We can use the two regime fits to examine the full salt concentration, valence, and force dependence of S(k). The data are shown in Figure 7. We use the same ranges as given above for fitting the low and high k region in Figure 6 for all the data. The uncertainty in ν is ≈±0.01. Considering first the νh values, the dependence on f has the same form for all cs and salt valence. At low f, νh is constant, equal to the equilibrium f = 0 value and then increases toward 1. (The f = 0 values are listed in Table 1.) The differences in the values of νh at low f is the expected effect of screening by the salt yielding smaller values for larger concentrations. The transition between the constant and increasing regions occurs at about f = kBT/lB for the monovalent salt and a slightly smaller force for the divalent salt. This force magnitude is the point where the applied force is comparable to the electrostatic and entropic interactions on the short length scales, which is roughly where such a transition should occur. In contrast to the high k behavior, the νl increases monotonically for the range of f calculated. Once f is about 3ε/ σ, then νl has reached 1 (within uncertainty) and remains 1 for
larger f. There is a much larger variation at low f in νl than in νh. The main source of this variation is that the equilibrium value of νl depends much more strongly on the salt concentration (see Table 1). For most cases, νl is always greater than νh. However, for the monovalent system at 100 mM νl is smaller than νh, which makes sense since at this salt concentration the screening is rather strong. In the divalent systems, νl > νh for all salt concentrations, though the trend is such that this would likely change at higher concentrations. Overall, the picture of structural change as a function of f that appears from Figure 7 is the following. The low k value of ν increases immediately upon application of a nonzero force. This implies that on long length scales the chain is being oriented and stretched out. In examining the images in Figure 5, there is a bending and orienting of long chain segments, particularly when comparing parts b and c. The short-ranged, local structure with its bends and ion aggregation is not strongly affected in the low force regime. Once the segments have been oriented and stretched close to linear (i.e., νl reaches ≈0.9) the short-range structure begins to straighten. As f increases in this regime, most of compliance is due to the short-range structure. This large f regime corresponds to the logarithmic scaling regime for Rz (cf. Figure 4). Consequently, the logarithmic regime corresponds to a straightening of the short-ranged structure of the polyelectrolyte chain. As Figure 4 shows, most of the monovalent system data are in the logarithmic regime. In order to get better data particularly at low forces, we have have performed simulations using the screened Coulomb potential (eq 5) for 1:1 salt systems. Since no free ions are treated in this approximation, the simulations are 2−3 orders of magnitude faster. We expect the screened Coulomb interaction to be an acceptable approximation for monovalent salt at high salt. The screened Coulomb data will differ from the explicit ion data at sufficiently low screening primarily because counterion condensation requires a renormalization of the Debye screening length used in the screened Coulomb data. We find that at 200 5762
dx.doi.org/10.1021/ma300899x | Macromolecules 2012, 45, 5757−5765
Macromolecules
Article
“unfolding”; i.e., long segments of the chain are bending apart and straightening. Figure 9 shows images of the (explicit salt)
mM the difference between the screened Coulomb data and the explicit ion data is negligible. At 100 mM there is a noticeable difference between the two methods, but the scaling behavior is matched and the force−extension data have a similar logarithmic regime. With the screened Coulomb systems, we are able to obtain data at lower forces and see the transition to the logarithmic regime more clearly. In Figure 8, we show the data at cs = 100 and 200 mM. The extension data have been fit using least-squares at low f to Rz ∼
Figure 9. Images of the chain for the monovalent salt case at 50 mM and f = 0.05ε/σ showing the folds and wrinkles. The numbers label the intermediate scale folds in the conformation in the right image, which are correlated with ion clustering in the left image. All ions within 2σ are shown in the left image. The monomers are blue and the positive salt ions are red; the equivalent counterions (if any) are green; negative salt ions are cyan.
monovalent 50 mM chain. The bottom part shows only the chain conformation which has several local “folds” labeled by numbers. In these folds monomers separated by at least several other monomers along the backbone come almost into contact. The top image showing the nearby ions indicates that these folds are related to clustering of ions along the chain which enable large bends. At sufficiently large forces νh starts to increase; this corresponds to an “unwrinkling” of the chain in which the short length scale, small bends begin to be pulled taut. The image of the chain conformation in Figure 9 clearly shows these wrinkles in the conformation. The transition between algebraic scaling and logarithmic dependence of Rz on f is bracketed by the transition in the ν values. The increase in νl occurs a bit below the transition, and the increase in νh occurs a bit above the transition. The primary dynamics occurring in the logarithmic regime is the unwrinkling. While there is some increase in νl, by the point where νh starts to increase, νl has reached 0.9. Thus, most of the unfolding has already occurred and the chain is rather straight on long length scales.
Figure 8. Force−extension data (top) and the scaling exponents from S(k) (bottom) for the screened Coulomb simulations at cs = 100 mM (squares) and cs = 200 mM (circles). In the top panel, the solid points are Rz and the open points are the full end-to-end distance R. The red lines are least-squares fits to the low f elastic region Rz ∼ f 0.80, and the green lines are least-squares fits of a logarithm to the high f region. In the bottom panel, the points connected by the solid line designate νh for the low k regime and the dotted line designates νl for the high k regime.
fγ and at large f as log f. The figure shows line fits for both regimes, which clearly display that two functional forms are required to fit the full range of data. In the top panel, the full end-to-end distance data are displayed. The important feature is that R is a constant equal to the equilibrium f = 0 value for a range of small f. In addition, at small f both ν are also constant. Both R and νl start to increase at about f = 0.05ε/σ. νh remains constant to much larger f, about 0.40ε/σ. These results imply the chain is not being stretched by the applied force in the low f limit. Since Rz is changing at small force, all that is occurring is an orientation of the set of chain conformations by the applied force. An unperturbed polymer has an anisotropic shape on average. For a neutral chain the three eigenvalues of the radius of gyration tensor ei have ratios e1/e2 ≈ 2.5 and e1/e3 ≈ 11.8.31 This implies a flattened cigar shape for the chain. For polyelectrolytes, the anisotropy is even larger.2 The low force data imply that the statistical ensemble of configuration is being oriented by the force, but not being stretched. Direct examination of the ei at f > 0 confirms this. Beyond this low-force, orientation regime, the chain does begin to stretch as indicated by the rise in R and νl. The values of νh remain constant to larger forces, indicating that the polymer is not being stretched on short length scales. This difference in behavior of νl and νh implies that the chains is
V. DISCUSSION The simulations have produced a picture of the effect of tensile force on a single flexible, strongly charged polyelectrolyte. There are three stages as force increases: (i) At the lowest forces, we find that the zero-force chain conformations are simply aligning with the force. (ii) At intermediate force, the long-range structure of the chain is affected via unfolding. (iii) As the force grows larger and the chain approaches its contour length, the force serves to straighten local wrinkles in the polymer. These stages correspond to structures on different length scales, with the larger length scales affected by the weaker forces. Here we discuss the nature of each of these three regimes and their connection to polymer theory and previously observed experimental results. At the smallest forces, the zero-force conformation of the polymer is simply being reoriented in the direction of the force. This regime is only fully observable in the simulations with electrostatics represented by the Yukawa potential as large fluctuations of the polymer at these low forces impose computationally expensive simulations in the explicit ion case. In this regime we observe that the force−extension curves appear to follow a power law steeper than both that observed in 5763
dx.doi.org/10.1021/ma300899x | Macromolecules 2012, 45, 5757−5765
Macromolecules
Article
experiment,13 which is predicted with the scaling theory for long polymers.22 The global reorientation of the chain may be chain length dependent.32 The (zero force) simulation systems lack the random walk regime because of short chain length. For sufficiently long polyelectrolytes the electrostatic interactions should be screened on the long length scales, and we should find ν = 0.588 in the single chain structure factor. While we find this scaling for neutral chains of the same length, for charged chains even at high salt concentrations with short screening lengths, we find that ν > 0.6. We have performed preliminary simulations using the Yukawa potential for larger chain lengths and find that neutral scaling is present for N ≳ 1500. It remains to be determined for such long chains with a neutral scaling regime if the low force limit will be the reorientation or the scaling found in experiment. At intermediate forces, the chain on long length scales (tens of monomers) is unfolded. This long length scale structure depends strongly on the salt concentration, since the salt strongly screens the interactions on this length scale. At low salt, the folds occur from the bending of the relatively straight segments, which are related to the standard picture of the persistence length. As the salt increases, these segments become progressively more compact, but for the salt concentration range studied, the segment structure is extended in comparison to neutral structure. Thus, the degree of folding depends on the salt concentration. The key aspect with respect to the pulling force is that these segments have a net electrostatic repulsion due to the screened Coulomb potential. Thus, the pulling force is actually aided by electrostatics in this regime and is predominantly overcoming entropy in elongating the folds. In the experiments the longer chain length has the additional neutral scaling regime in the structure, and these segments should be even weaker to unfold and occur at lower forces. At the largest forces, wrinkles that occur on the length scale of several monomers are being straightened. Indicators of wrinkles have been observed in structure factors of early polyelectrolyte simulations,2,33 which had the high k signature of the wrinkles (ν ≠ 1), but have not been a focus of consideration. More recent calculations of bond−bond correlation functions are also indicative of the wrinkles.29,30,34 These calculations of bond−bond correlation for flexible chains show a large drop in the first several bond lengths, which corresponds to the wrinkled structure on short length scales. We note that the effect on the bond−bond correlation function causes a problem with the standard picture of persistence length. The definition of Lp is the decay length of the bond− bond correlation function, but that is very short due to the large drop. On a longer, intermediate length scale there is a correlation that better corresponds to the standard picture of the persistence length.30 The actual picture of chain segments of length Lp is not one of a simple, straight segment, but one of neighboring bonds zigzagging in a wrinkled structure. The source of the wrinkling is within the combination of electrostatic and entropic interactions. In strong, flexible polyelectrolytes there are a large number of condensed counterions. The electrostatic energy is locally reduced by putting a counterion as close to multiple monomers, which is achieved by wrapping the chain about a counterion (cf. Figure 9). In addition, there is a complex interplay with entropy. On one hand localizing counterions lowers entropy; this effect limits the amount of screening that counterions can do. On the other hand, entropy is increased by having as many chain
conformations as possible, which is promoted by the wrinkling relative to a straight chain. The full combination of repulsive (monomer−monomer) and attractive (monomer−counterion) Coulomb interactions with ion and chain entropy yield the short-range structure. The screened Coulomb simulations are quite accurate for the monovalent system at 200 mM, but there are no counterions in these simulations for the chain to wrap around. Without a wrapping effect, how can the screened Coulomb interaction yield the correct wrinkled structure, i.e., the correct value of νh and the logarithmic elastic regime? The counterion wrapping effectively reduces the direct monomer−monomer Coulomb repulsion. In the screened Coulomb simulations, the monomer−monomer repulsion is reduced by the exponential screening factor. For the monovalent system at 200 mM, this screening produces an effective monomer−monomer interaction that matches the explicit ion monomer−monomer interactions. As the salt concentration is lowered, a noticeable and increasing difference occurs between the two models with the explicit ion simulations having a smaller νh or less wrinkled structure. Thus, the screened Coulomb approximation captures the wrinkled structure over a limited range of salt concentrations. Given the presence of the folds and wrinkles, why does unfolding occur before unwrinkling in the force−extension dynamics? The fold structure occurs on length scales where screening of the electrostatic interactions is strong. There is not a segment−segment attraction that needs to be overcome to unfold. The wrinkle structure occurs where screening is weak. Fundamentally, the wrinkles are due to the aggregation of oppositely charged ions, which have strong Coulomb interactions. To unwrinkle requires altering the attractions among the ions and therefore requires a larger force. Evidence for this can be found in the difference between the monovalent and divalent systems at large forces. The reason for this difference is that the local ion binding depends strongly on the salt ion valence. This is partly seen in the values of νh which are different for the two valences. The values of νh for the divalent systems are smaller than for the monovalent systems, indicating that the short-range scale structure is more compact and wrinkled. Calculations of the average charge density about a monomer also show a difference between the monovalent and divalent systems. As expected, the divalent systems have a larger opposite charge neighboring the monomer. These differences are expected given that the divalent ions will pack more charged monomers around them to lower the electrostatic energy and thereby produce a more wrinkled structure. The force to unwrinkle the divalent system is thus commensurately larger.
VI. CONCLUSION We have used coarse-grained molecular dynamics to study the stretching of a single strongly charge, flexible polyelectrolyte under an applied force. Both monovalent and divalent added salt have been treated. The simulation data overlaps experimental data on ssDNA,13,15 reproducing in particular the logarithmic elastic regime in monovalent salt. In addition, the simulations recover the difference in mono- and divalent concentration required to modulate ssDNA elasticity: vastly smaller concentrations of divalent ions produce extensions similar to monovalent ions at similar forces. Thus, we have demonstrated that the simulation has all of the relevant physics to capture the effect of ions on ssDNA elasticity and that the experimental force−extension data should be interpreted as a 5764
dx.doi.org/10.1021/ma300899x | Macromolecules 2012, 45, 5757−5765
Macromolecules
Article
result of entropic and electrostatic effects of strong, flexible polyelectrolytes rather than a chemistry-specific property of the ion−polymer interaction. The conformational data available from the simulations reveal the underlying structural changes as a function of applied force. We attribute the logarithmic elastic regime to the straightening of the short length scale wrinkles that are associated with the bends in the local structure of the flexible polyelectrolyte about the condensed ions. This local, wrinkled structure has been found in previous simulations,30,33,34 yet it has not been shown to be significant in experiments until now. Finally, we note the recent appearance of an article in which the authors used the correct (nonexponential) bond−bond correlation function to analytically calculate a logarithmic force−extension relation.35
■
(19) Dobrynin, A.; Carrillo, J. Y.; Rubinstein, M. Macromolecules 2010, 9181−9190. (20) Livadaru, L.; Netz, R.; Kreuzer, H. Macromolecules 2003, 36, 3732−3722. (21) Toan, N.; Thirumalai, D. Macromolecules 2010, 43, 4394−4400. (22) Pincus, P. Macromolecules 1976, 9, 386−388. (23) Everaers, R.; Milchev, A.; Yamakov, V. Eur. Phys. J. E 2002, 8, 3−14. (24) Hehmeyer, O.; Stevens, M. J. Chem. Phys. 2005, 122, 134909. (25) Pollock, E. L.; Glosli, J. Comput. Phys. Commun. 1996, 95, 93. (26) As the simulation data were primarily focused on the high-force regime, we were forced to fix the low-force power law γ, choosing to fix it to the value measured for ssDNA.13 It should be noted that we do not measure a power law Rz ∼ fγ over a significant range of forces in the simulation data. (27) The experimental data extend to lower applied forces, which is important for obtaining the low force scaling regime. (28) Grest, G. S.; Kremer, K. Phys. Rev. 1986, A33, 3628. (29) Gubarev, A.; Carrillo, J.-M.; Dobrynin, A. Macromolecules 2009, 42, 5815−5860. (30) Carrillo, J.-M.; Dobrynin, A. Macromolecules 2011, 44, 5798− 5816. (31) Bishop, M.; Michele, J. M. J. J. Chem. Phys. 1986, 84, 444. (32) Morrison, G.; Hyeon, C.; Toan, N.; Ha, B.; Thirumalai, D. Macromolecules 2007, 40, 7343−7353. (33) Stevens, M.; Kremer, K. Macromolecules 1993, 26, 4717. (34) Nguyen, T.; Shklovskii, B. I. Phys. Rev. E 2002, 66, 021801. (35) Toan, N.; Thirumalai, D. J. Chem. Phys. 2012, 136, 235103.
AUTHOR INFORMATION
Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DE-AC0494AL85000. This work was performed at the US Department of Energy, Center for Integrated Nanotechnologies, at Los Alamos National Laboratory (Contract DE-AC52-06NA25396) and Sandia National Laboratories. O.S. and D.M. were partially supported by the NSF through Grant DMR-1006737. This work was started while M.S. and O.S. were attending a program at the Aspen Center for Physics supported by the NSF Grant 1066293, and the collaboration was furthered by M.S. visiting the Kavli Institute for Theoretical Physics, which was supported in part by the National Science Foundation under Grant NSF PHY11-25915.
■
REFERENCES
(1) Khokhlov, A.; Khachaturian, K. Polymer 1982, 23, 1742. (2) Stevens, M.; Kremer, K. J. Chem. Phys. 1995, 103, 1669. (3) Dobrynin, A. V.; Colby, R. H.; Rubinstein, M. Macromolecules 1995, 28, 1859−1871. (4) Barrat, J.-L.; Joanny, J.-F. Adv. Chem. Phys. 1996, 94, 1. (5) Grosberg, A.; Nguyen, T.; Shklovskii, B. Rev. Mod. Phys. 2002, 74, 329−345. (6) Netz, R.; Andelman, D. Phys. Rep. 2003, 380, 1−95. (7) Holm, C.; Joanny, J.; Kremer, K.; Netz, R.; Reineker, P.; Seidel, C.; Vilgis, T.; Winkler, R. Adv. Polym. Sci. 2004, 166, 67−111. (8) Dobrynin, A. V.; Rubinstein, M. Prog. Polym. Sci. 2005, 30, 1049− 1118. (9) Stevens, M. Biophys. J. 2001, 80, 130. (10) Moreira, A. G.; Netz, R. Europhys. Lett. 2000, 52, 705−711. (11) Dessinges, M.; Maier, B.; Zhang, Y.; Peliti, M.; Bensimon, D.; Croquette, V. Phys. Rev. Lett. 2002, 24, 248102. (12) Seol, Y.; Skinner, G.; Visscher, K. Phys. Rev. Lett. 2004, 93, 118102. (13) Saleh, O. A.; McIntosh, D. B.; Pincus, P.; Ribeck, N. Phys. Rev. Lett. 2009, 102, 068301. (14) McIntosh, D.; Ribeck, N.; Saleh, O. A. Phys. Rev. E 2009, 80, 041803. (15) McIntosh, D.; Saleh, O. A. Macromolecules 2011, 44, 2328− 2333. (16) Marko, J.; Siggia, E. Macromolecules 1995, 28, 8759−8770. (17) Barrat, J.-L.; Joanny, J.-F. Europhys. Lett. 1993, 3, 343. (18) Carrillo, J. Y.; Dobrynin, A. Macromolecules 2010, 43, 2589− 2604. 5765
dx.doi.org/10.1021/ma300899x | Macromolecules 2012, 45, 5757−5765
