Effect of Charge Regulation and Conformational Equilibria in the Stretching Properties of Weak Polyelectrolytes

Weak polyelectrolytes can modulate their charge in response to external perturbations, such as changes in the pH, ionic strength (I), or electrostatic interactions with other charged species, a phenomenon known as charge regulation (CR). On the other hand, it is well established that CR is highly coupled with the conformational degrees of freedom. In this paper, the influence of CR in the stretching properties of weak polyelectrolytes is analyzed, and the possibility of CR induced by mechanical stretching is explored. With this aim, we make use of a minimal model, which captures the fundamental aspects present in the stretching of a flexible weak linear polyelectrolyte: internal angle rotation, bond stretching, bond bending, and proton binding, which is the paradigmatic mechanism of CR. The angle rotation is described by using the rotational isomeric state approximation, while for protonation, the site binding model is assumed. Mechanical stretching is studied by performing semi-grand canonical Monte Carlo simulations at different pH and ionic strength conditions. The simulations simultaneously provide both conformational (bond state probabilities, persistence length lP, and chain elongation) and protonation properties (degree of protonation θ and the effective protonation constant Kc). The obtained force−extension curves suggest that the pH value and the ionic strength I have a significant effect on polyelectrolyte stretching. Three different force regimes can be observed. For large forces (F > 100 pN for typical force constants), the force−extension curve is almost independent of the pH and I. For low forces, the persistence length lP is force-independent, although it strongly increases with the pH value. Under this regime, linear and Pincus scaling behaviors are observed. Finally, in the intermediate-force regime, both rotational and protonation degrees of freedom are mechanically activated, and the picture becomes more complicated. It is found that lP increases with F and, under certain conditions, a significant increase of θ with F is observed, indicating that CR could in principle be induced by means of mechanical stretching. This fact can be explained by analyzing the coupling between θ and the probability of a bond to be in the gauche state P(g). P(g) decreases with F as the bonds adopt the trans conformation so that the electrostatic repulsion is reduced and θ increases. Finally, the intricate interplay between short-range and long-range interactions is analyzed, leading to apparently contradictory behaviors (P(g) and lp simultaneously decrease with I), which can only be explained by CR and the presence of complex spatial correlations. ■ INTRODUCTION In the last two decades, the development of single-molecule force spectroscopy has led to an extraordinary expansion of the field of mechanochemistry. By applying a controlled external force to a molecule that is chemically attached to a surface, a wide range of mechanically induced physicochemical events can take place. Just to mention a few examples, AFM has been used to mechanically induce cis-to-trans isomerization of carbon−carbon double bonds, prolyl cis−trans isomerization, or conformational chair−boat transitions or hydrogen bond breaking in polysaccharides. Some ring-opening reactions, normally forbidden by orbital symmetry, become possible if a tensile force is applied to the polymer chain. Single-molecule AFM experiments have been recently used in Received: June 12, 2019 Revised: September 16, 2019 Published: October 16, 2019 Article pubs.acs.org/Macromolecules Cite This: Macromolecules 2019, 52, 8017−8031 © 2019 American Chemical Society 8017 DOI: 10.1021/acs.macromol.9b01160 Macromolecules 2019, 52, 8017−8031 D ow nl oa de d vi a U N IV O F C IN C IN N A T I on F eb ru ar y 12 , 2 02 1 at 0 3: 51 :2 7 (U T C ). Se e ht tp s: //p ub s. ac s. or g/ sh ar in gg ui de lin es f or o pt io ns o n ho w to le gi tim at el y sh ar e pu bl is he d ar tic le s. monitoring force-dependent enzyme catalysis and surface desorption of polypeptides; the characterization of new supramolecular polymers based on host-enhanced π−π interaction; or the design of mechanophores embedded in macrocycles, which allows pinpointing of the mechanochemical bond rupture. Optical tweezers have also been used in the study of the elastic properties of biomacromolecules such as single-stranded DNA (ss-DNA). In parallel to the experimental work, several theoretical approaches have been developed, differing in the detail of description of the macromolecular structure. In the freely jointed chain (FJC) model, the polymer chain is represented at the coarse-grain level by a set of rigid links joined with fully random orientations. Although able to account for the stretching properties of a wide variety of synthetic polymers with different structures and solvents, this model was shown to present clear deviations from the elastic response of many other macromolecules of interest, such as double-stranded DNA (ds-DNA). Aiming at overcoming these limitations, Marko and Siggia modeled the polymer as a worm-like chain (WLC), which, assuming exponential decaying correlations between chain segments, accounted for the capability of the chain to deform on short-length scales. The resulting highforce regime matched very well with a variety of polymers for which electrostatic interactions can be neglected, such as some polypeptides or ds-DNA. Models including freely rotating bonds, bond elasticity, or ligand−receptor equilibria have also been proposed, leading to theoretical predictions of new force−extension regimes. As an alternative to these coarse-graining approaches, theoretical methods based on first principles, which account for the detailed atomistic structure of the macromolecular backbone, have also been proposed. In most of these studies, ab initio calculations are first performed in order to detect the more stable conformational states of the interacting monomers at different elongations of the bonds. Once the structural microscopic information is available, the necessary thermal averages are performed by using Monte Carlo (MC) or transfer matrix techniques. The resulting scheme has been successful in reproducing the experimental force−extension curves of several polymers. In particular, the stretching behavior of poly(ethylene glycol) (PEG) has been analyzed in detail. In essence, this methodology can be regarded as a generalization of the rotational isomeric state (RIS) model developed mainly by Flory to study the conformational properties of linear chains in which only the rotational states of minimum energy (commonly trans, gauche+, and gauche−) are taken into account in the computation of the thermal averages. The methods mentioned above only account for short-range interactions and cannot thus be applied to charged macromolecules for which the long-range Coulombic forces cannot be neglected. The presence of self-avoiding electrostatic forces produces, however, new elastic regimes, which strongly deviate from the ideal, non-interacting FJC and WLC models. The resulting stretching behavior, mainly studied in single-stranded nucleic acids (ss-DNA and ss-RNA), is extremely dependent on the valence and concentration of the counterions and seems to be well explained by the recently proposed “snake chain model”. This model was motivated by recent MC simulations with explicit ions, which suggested two elastic regimes. At low forces, the polyelectrolyte behaves as a set of swollen electrostatic blobs on a long-length scale, while at high forces, a short-length, ion-stabilized, crumpling structure is detected. Since ss-DNA and ss-RNA are strong polyelectrolytes, in those studies, the macromolecular charge is considered constant and independent of the degree of stretching. However, this could not be the case of weak polyelectrolytes for which the charge is in general a dynamical and fluctuating variable. This fact leads to the phenomenon of charge regulation (CR), defined as the capability of weak polyelectrolytes to modulate their ionization state as a response to some physicochemical perturbation. The main aim of the present work is to study the influence of CR in the stretching properties of weak polyelectrolytes. The possibility of mechanically induced CR will also be explored. CR is ubiquitous in a wide range of processes of biological, environmental, and technological interest. A few examples are the stability of colloidal systems and nanoparticle coatings, receptor−ligand interactions in biochemical systems, and protein−protein and protein−surface interactions, among many others, which can be found in ref 52 and references quoted therein. The paradigmatic mechanism for CR is the binding of protons and other small ions present in the backward medium. Although CR can take place in rigid structures such as nanoparticles or surfaces, most polyelectrolytes are flexible so that alterations in the ionization state induce changes in the rotational states of the bonds. Sometimes, this can even produce dramatic conformational transitions in the global macromolecular structure, such as the helix−coil transitions of polypeptides or the abrupt swelling of poly(methacrylic acid) in a very narrow range of pH values. The mechanochemistry of weak polyelectrolytes is still a fairly unexplored area from the experimental, theoretical, and computational point of view. Although some AFM experiments have been performed on weak polyelectrolytes (such as hyaluronic acid), they have been either focused on the temperature effect or carried out at pH conditions where CR is negligible. As a result, the effect of other environmental variables such as pH or salt concentration is still unknown. In this work, we introduce a minimal model that captures the fundamental aspects present in the stretching of a flexible weak linear polyelectrolyte: internal angle rotation, bond stretching, bond bending, and proton binding. The model presented is based on the site binding rotational isomeric state (SBRIS) model. The model is implemented in a Monte Carlo simulation scheme in the semi-grand canonical ensemble (SGCMC) widely used in computational modeling of CR phenomena. An outline of the used methodology is reported in Model and Simulations, while the Results and Discussion section is devoted to analyzing the behavior of both the conformational (bond state probabilities, persistence length, and chain extension) and protonation properties (degree of protonation θ). ■ MODEL AND SIMULATIONS Minimal Site Binding Rotational Isomeric State (SBRIS) Model of Stretched Weak Polyelectrolytes. In this work, we will make use of a model, which, containing a minimum number of parameters, still captures the fundamental aspects present in the stretching of a flexible weak linear polyelectrolyte: internal angle rotation, bond stretching, bond bending, and proton binding. The polyelectrolyte, outlined in Figure 1a, can be considered a simplification of a previously Macromolecules Article DOI: 10.1021/acs.macromol.9b01160 Macromolecules 2019, 52, 8017−8031 8018 proposed model for linear poly(ethylenimine) (LPEI). A similar model has been recently proposed to study the role of long-range interactions in the conformational/protonation coupling. Let us assume that the chain is symmetric (i.e., the chain has a plane of symmetry when it is completely elongated), thus avoiding the question of tacticity, and contains a protonating site situated every three chain positions. In Figure 1, inert and protonating sites are depicted in gray and blue, respectively. A macromolecule with N protonating sites thus contains M = 3N − 3 bonds. The protonation equilibria are treated using the site binding (SB) model for which the protonating sites can adopt two possible states: protonated (dark blue) and deprotonated (cyan). Within the SB approach, the ionization state of the macromolecule can be characterized by a set of variables s = {si}, i = 1, ..., N, with values 0 (deprotonated) or 1 (protonated). On the other hand, the conformational degrees of freedom are treated assuming the rotational isomeric state (RIS) approximation; that is, only the rotational states corresponding to local energy minima are taken into account, typically trans (t), gauche+ (g), and gauche− (g−). A conformational state of the macromolecule can now be defined by a set of variables c = {cj}, j = 1, ..., M where cj is the row vector with as many components as the number of states can adopt the bond j. In our case, each cj can only take three different values: c (1, 0, 0) j = if bond j is in trans, c (0, 1, 0) j = if it is in gauche+, and c (0, 0, 1) j = if it is in gauche−. Since the chain is symmetric, gauche+ and gauche− states have the same energy. For simplicity, let us also assume that only the bonds with adjacent protonating sites (c bonds in Figure 1a) are allowed to rotate, while the rest of the bonds (a and b) are forced to be in the trans state. This approximation has been previously used in the modeling of stretching properties of poly(ethylene glycol) (PEG), the neutral counterpart of the model proposed here. As a result, in our model, only N − 1 bonds from the total M bonds are allowed to rotate. Finally, we introduce the possibility of elastic bond stretching and bending. Combining the SB and RIS approaches, we obtain the SBRIS model, which deals with ionization and conformational equilibria simultaneously. The resulting free energy can be expressed as the sum of five contributions W length angle SR LR = + + + + (1)


■ INTRODUCTION
In the last two decades, the development of single-molecule force spectroscopy has led to an extraordinary expansion of the field of mechanochemistry. 1,2 By applying a controlled external force to a molecule that is chemically attached to a surface, a wide range of mechanically induced physicochemical events can take place. Just to mention a few examples, AFM has been used to mechanically induce cis-to-trans isomerization of carbon−carbon double bonds, 3 prolyl cis−trans isomeriza-tion, 4,5 or conformational chair−boat transitions or hydrogen bond breaking in polysaccharides. 6,7 Some ring-opening reactions, normally forbidden by orbital symmetry, become possible if a tensile force is applied to the polymer chain. 8,9 Single-molecule AFM experiments have been recently used in monitoring force-dependent enzyme catalysis 10 and surface desorption of polypeptides; 11,12 the characterization of new supramolecular polymers based on host-enhanced π−π interaction; 13 or the design of mechanophores embedded in macrocycles, which allows pinpointing of the mechanochemical bond rupture. 14 Optical tweezers have also been used in the study of the elastic properties of biomacromolecules such as single-stranded DNA (ss-DNA). 2,15 In parallel to the experimental work, several theoretical approaches have been developed, differing in the detail of description of the macromolecular structure. 1,2,16 In the freely jointed chain (FJC) model, the polymer chain is represented at the coarse-grain level by a set of rigid links joined with fully random orientations. Although able to account for the stretching properties of a wide variety of synthetic polymers with different structures and solvents, 1 this model was shown to present clear deviations from the elastic response of many other macromolecules of interest, such as double-stranded DNA (ds-DNA). Aiming at overcoming these limitations, Marko and Siggia modeled the polymer as a worm-like chain (WLC), which, assuming exponential decaying correlations between chain segments, accounted for the capability of the chain to deform on short-length scales. 17 The resulting highforce regime matched very well with a variety of polymers for which electrostatic interactions can be neglected, such as some polypeptides 18 or ds-DNA. 19 Models including freely rotating bonds, 20 bond elasticity, 21 or ligand−receptor equilibria 22 have also been proposed, leading to theoretical predictions of new force−extension regimes.
As an alternative to these coarse-graining approaches, theoretical methods based on first principles, which account for the detailed atomistic structure of the macromolecular backbone, have also been proposed. 23 In most of these studies, ab initio calculations are first performed in order to detect the more stable conformational states of the interacting monomers at different elongations of the bonds. 24 Once the structural microscopic information is available, the necessary thermal averages are performed by using Monte Carlo (MC) or transfer matrix techniques. 25, 26 The resulting scheme has been successful in reproducing the experimental force−extension curves of several polymers. 20,27−29 In particular, the stretching behavior of poly(ethylene glycol) (PEG) has been analyzed in detail. 27,30,31 In essence, this methodology can be regarded as a generalization of the rotational isomeric state (RIS) model developed mainly by Flory to study the conformational properties of linear chains 32,33 in which only the rotational states of minimum energy (commonly trans, gauche+, and gauche−) are taken into account in the computation of the thermal averages.
The methods mentioned above only account for short-range interactions and cannot thus be applied to charged macromolecules for which the long-range Coulombic forces cannot be neglected. The presence of self-avoiding electrostatic forces produces, however, new elastic regimes, which strongly deviate from the ideal, non-interacting FJC and WLC models. 34 The resulting stretching behavior, mainly studied in single-stranded nucleic acids (ss-DNA and ss-RNA), is extremely dependent on the valence and concentration of the counterions 2,15,35−39 and seems to be well explained by the recently proposed "snake chain model". 40,41 This model was motivated by recent MC simulations with explicit ions, 42,43 which suggested two elastic regimes. At low forces, the polyelectrolyte behaves as a set of swollen electrostatic blobs on a long-length scale, while at high forces, a short-length, ion-stabilized, crumpling structure is detected. 41−44 Since ss-DNA and ss-RNA are strong polyelectrolytes, in those studies, the macromolecular charge is considered constant and independent of the degree of stretching.
However, this could not be the case of weak polyelectrolytes for which the charge is in general a dynamical and fluctuating variable. This fact leads to the phenomenon of charge regulation (CR), defined as the capability of weak polyelectrolytes to modulate their ionization state as a response to some physicochemical perturbation. 45,46 The main aim of the present work is to study the influence of CR in the stretching properties of weak polyelectrolytes. The possibility of mechanically induced CR will also be explored. CR is ubiquitous in a wide range of processes of biological, environmental, and technological interest. A few examples are the stability of colloidal systems and nanoparticle coatings, 47,48 receptor−ligand interactions in biochemical systems, 49 and protein−protein 50 and protein−surface interactions, 51 among many others, which can be found in ref 52 and references quoted therein. The paradigmatic mechanism for CR is the binding of protons and other small ions present in the backward medium. Although CR can take place in rigid structures such as nanoparticles or surfaces, most polyelectrolytes are flexible so that alterations in the ionization state induce changes in the rotational states of the bonds. Sometimes, this can even produce dramatic conformational transitions in the global macromolecular structure, such as the helix−coil transitions of polypeptides 53 or the abrupt swelling of poly(methacrylic acid) in a very narrow range of pH values. 54 The mechanochemistry of weak polyelectrolytes is still a fairly unexplored area from the experimental, theoretical, and computational point of view. Although some AFM experiments 7,55 have been performed on weak polyelectrolytes (such as hyaluronic acid), they have been either focused on the temperature effect 7 or carried out at pH conditions where CR is negligible. 55 As a result, the effect of other environmental variables such as pH or salt concentration is still unknown. In this work, we introduce a minimal model that captures the fundamental aspects present in the stretching of a flexible weak linear polyelectrolyte: internal angle rotation, bond stretching, bond bending, and proton binding. The model presented is based on the site binding rotational isomeric state (SBRIS) model. 52,56−58 The model is implemented in a Monte Carlo simulation scheme in the semi-grand canonical ensemble (SGCMC) widely used in computational modeling of CR phenomena. 51,59−66 An outline of the used methodology is reported in Model and Simulations, while the Results and Discussion section is devoted to analyzing the behavior of both the conformational (bond state probabilities, persistence length, and chain extension) and protonation properties (degree of protonation θ).

■ MODEL AND SIMULATIONS
Minimal Site Binding Rotational Isomeric State (SBRIS) Model of Stretched Weak Polyelectrolytes. In this work, we will make use of a model, which, containing a minimum number of parameters, still captures the fundamental aspects present in the stretching of a flexible weak linear polyelectrolyte: internal angle rotation, bond stretching, bond bending, and proton binding. The polyelectrolyte, outlined in Figure 1a, can be considered a simplification of a previously  57 A similar model has been recently proposed to study the role of long-range interactions in the conformational/protonation coupling. 52 Let us assume that the chain is symmetric (i.e., the chain has a plane of symmetry when it is completely elongated), thus avoiding the question of tacticity, and contains a protonating site situated every three chain positions. In Figure 1, inert and protonating sites are depicted in gray and blue, respectively. A macromolecule with N protonating sites thus contains M = 3N − 3 bonds.
The protonation equilibria are treated using the site binding (SB) model for which the protonating sites can adopt two possible states: protonated (dark blue) and deprotonated (cyan). Within the SB approach, the ionization state of the macromolecule can be characterized by a set of variables s = {s i }, i = 1, ..., N, with values 0 (deprotonated) or 1 (protonated). On the other hand, the conformational degrees of freedom are treated assuming the rotational isomeric state (RIS) approximation; 32,33 that is, only the rotational states corresponding to local energy minima are taken into account, typically trans (t), gauche+ (g + ), and gauche− (g − ). A conformational state of the macromolecule can now be defined by a set of variables c = {c j }, j = 1, ..., M where c j is the row vector with as many components as the number of states can adopt the bond j. In our case, each c j can only take three different values: Since the chain is symmetric, gauche+ and gauche− states have the same energy. For simplicity, let us also assume that only the bonds with adjacent protonating sites (c bonds in Figure 1a) are allowed to rotate, while the rest of the bonds (a and b) are forced to be in the trans state. This approximation has been previously used in the modeling of stretching properties of poly(ethylene glycol) (PEG), 27 the neutral counterpart of the model proposed here. As a result, in our model, only N − 1 bonds from the total M bonds are allowed to rotate. Finally, we introduce the possibility of elastic bond stretching and bending.
Combining the SB and RIS approaches, we obtain the SBRIS model, which deals with ionization and conformational equilibria simultaneously. 52,56,57 The resulting free energy can be expressed as the sum of five contributions The term represents the mechanical work exerted by the applied force F, which is considered to act in the z axis direction, as shown in Figure 1, and r is the polyelectrolyte chain end-to-end vector. length and angle quantify the elastic deformation of the length and the angles of the M bonds, respectively, which can be important at large forces. In this work, they are represented by the harmonic potentials where l j , α j , l j,0 , and α j,0 represent, respectively, the length, the bending angle, the equilibrium length, and the equilibrium bending angle of bond j. Finally, k length, j and k angle, j denote the bond stretching and bending force constants. Note that the geometrical parameters and the force constants in the potentials (eqs 3 and 4) do not depend on the ionization state of the sites. At this level of description, this is a reasonable approximation, according to quantum-mechanical computations, which show only small variations in the bond lengths (see, for instance, the results for LPEI at different degrees of ionization 67 ). On the other hand, as will be shown in the next section, the bond bending and bond stretching will be essentially induced by the mechanical work at high enough forces rather than by electrostatic repulsions. The electrostatic/conformational interaction free energy has been split into short-range (SR) SR and long-range (LR) LR contributions, as depicted in Figure 1. This distinction becomes necessary due to the fundamental differences in the physical chemistry of SR and LR interactions. It is well established that LR interactions are chemically unspecific, Only the rotational states corresponding to minimum energy (rotational isomeric state approximation), that is, trans (t), gauche+ (g + ), and gauche− (g − ), are taken into account. In order to minimize the number of parameters, only the bonds holding protonable sites (c bonds) are allowed to rotate, while the rest of the bonds (a and b) are forced to be in the trans state. The macromolecular chain is considered symmetric so that g + and g − have the same energy ε σ . Two kinds of sites are considered: inert sites (gray) and protonating sites. The latter can adopt two possible states: protonated (dark blue) or deprotonated (cyan), with protonation constant K (site binding model). Long-range (LR) and short-range (SR) interactions are treated in a different way. LR interactions (mediated by the solvent) are described by the Debye−Huckel potential. Conversely, neighboring protonated sites linked by a c bond in the trans state interact with energy ε u, t since SR interactions are mainly mediated by the macromolecular skeleton. Two neighboring sites linked by a c bond in the gauche state cannot be simultaneously protonated due to the huge electrostatic repulsion (ε u, g → ∞). (b) Snapshot from a semi-grand canonical Monte Carlo (SGCMC) simulation of the stretching of a weak polyelectrolyte with pK = 9, pH = 6, I = 10 −3 M, F = 10 pN, ε σ = − 1, ε u, t = 1, l 0 = 1.5 Å, and α 0 = 120°. Both the extension and the degree of protonation depend on the applied force in the z axis direction. mediated by the solvent, and can be reasonably approximated by a simple pair-interaction continuous force field. Conversely, SR interactions between neighboring sites and bonds are mediated by the macromolecular skeleton so that they depend on the detailed chemical environment of the site. As a result, they cannot be described by simple potentials, and specific interaction parameters must be used. These parameters will depend on the particular rotational state of the bond connecting the two protonating sites (in our case, the c bonds). SR is the result of two contributions , corresponding to the classical RIS model, 32,33 represents the conformational energy of the bonds for a given conformational state c = {c j } when the polyelectrolyte is uncharged. On the other hand, s c ( , ) p includes the binding free energy and the SR electrostatic interaction between charged sites and accounts for the coupling between the conformational and ionization degrees of freedom. The term rot can be expressed as where β = 1/k B T is the inverse of the thermal energy, ϵ j is a row vector whose elements are the free energies associated with the rotational state of bond j, and E j is a square matrix containing the interaction energies between the neighboring bonds j and j + 1. ϵ and E are expressed in thermal units and divided by a factor ln 10 in order to be compared with the pH scale. The sum in eq 6 could be extended to take into account three-bond interactions, four-bond interactions, and so on. In this work, however, we only considered SR interactions involving the first neighbor bonds. Following Flory, 32,33 we choose as the ground state the configuration with all the bonds in the trans state. In this way, the rotational parameters ϵ j and E j can be expressed as where ε σ, j is the free energy of the gauche states, while ε ψ, j and ε ω, j are related to the interaction energies between two consecutive gauche states with the same and different orientation, respectively. The Boltzmann factor corresponding to ε α is precisely α, that is, ε α = − log α (α = σ, ψ, ω). In choosing this notation, the resulting Boltzmann factors are denoted by the same symbols used in previous works, 33,56−58 which make use of the transfer matrix approach. For instance, the transfer matrix corresponding to the free energy (eq 6) is given by On the other hand, the second term in eq 5, p , can be expressed in terms of the ionization state s = {s i } by means of the cluster expansion 52 where μ i = pH − pK i = − log (Ka H ) is the reduced chemical potential of the ionizable site i, which depends on the proton activity, a H , and the pK values of the protonation constant of site i, pK i . ϵ u is a row vector whose components correspond to the electrostatic interaction energy between the neighboring charged sites i and i + 1. The intensity of the interaction is determined by the rotational state of the c bond between the two sites (bond number 3i − 1 in the chain) In eq 10, we have assumed that only the bonds of type c are able to rotate. Note that in the SBRIS model, the SR electrostatic interactions depend on the conformation of the bond linking of the sites, which couples the ionization and the conformational degrees of freedom. As in eq 6, the sum in eq 9 could be extended to take into account triplet interactions, quadruplet interactions, and so on. In this work, however, only neighboring pair SR interactions will be taken into account.
Finally, as in most of the previous literature, 57,68 LR electrostatic interactions will be described by the Debye− Huckel (DH) potential where B ≈ 0.7 nm is the Bjerrum length in water at 298.15 K, d i j is the distance between sites i and j, and is the Debye length for water 298.15 K at ionic strength I.
Since we are interested in a model with the minimum number of parameters in order to analyze the fundamental aspects of the stretching, in this work, we will restrict ourselves to the special situation in which all the bonds have the same length, bond angle (l j,0 = l 0 and α j,0 = α 0 ), and force constants (k length, j = k length and k angle, j = k angle ). Moreover, we consider that all the protonating sites are identical (pK i = pK) and the possible end effects of the chain are neglected so that ϵ i = ϵ and ϵ u, i = ϵ u . It is also assumed that when two neighboring sites are charged, the very strong SR repulsion hinders the gauche conformation of the c bond so that u g = 0 or ε u, g → ∞.
Since one of every two consecutive bonds is always in the trans state, the interaction terms ε ψ and ε ω become irrelevant, and they can be taken as zero without loss of generality. To summarize, the model presented involves the following assumptions: 1. The SBRIS model is used to describe the conformational and protonation equilibria on the same foot. 2. The molecule contains one protonating site every two inert groups, as shown in Figure 1. 3. The ionizable sites are identical, with the same pK value, and the bonds have the same length, bending angle, and constant forces. 4. Only the bonds of type c are allowed to rotate, while bonds of types a and b are constrained to be in the trans state. In practice, this implies that the rotation of the bonds is independent when the macromolecule is fully uncharged. 5. LR interactions are described by the DH potential, which accounts for screening effects so that co-and counterions intervene only in an effective way. Only excluded volume effects induced by electrostatics are taken into account. 6. Specific parameters are used to describe interactions between neighboring sites. Moreover, when two neighboring sites are simultaneously protonated, the c bond linking the sites cannot adopt the gauche state. 7. As a result, the parameters involved in the model are ε σ (free energy difference between gauche and trans states), interaction energy between neighboring sites (ε u, t ) when the c bonds are in the gauche state, equilibrium length and equilibrium angle of the bonds (l 0 and α 0 ), and constant forces for the bending and bond stretching (k length and k angle ). The control variables are the reduced chemical potential of the protonating sites (μ = pH − pK) and the ionic strength (I). Finally, note that in the absence of LR interactions ( 0 LR = ), that is, at high enough ionic strengths, the model can be exactly solved by using the transfer matrix (TF) method. 52,57,58 When applied to calculate stretching properties, the resulting TF combine conformational energies, by means of TF of the type (eq 8), and geometrical restrictions imposed by the macromolecular skeleton. In the absence of ionization processes, this method has been used to study the stretching of chains with freely rotating bonds 20 and the stretching properties of POE. 20,27 However, in this work, we are especially interested in the effect of electrostatic interactions that are long ranged, and the LR term (eq 11) cannot be neglected. As a consequence, the transfer matrix approach would be too restrictive since only the cases of high ionic strength could be analyzed. For this reason, a Monte Carlo computational code has been developed in order to implement the model, which is described in the next subsection.
Monte Carlo Simulations at Constant pH Value. The proposed SBRIS model is analyzed by means of simulations in the semi-grand canonical Monte Carlo (SGCMC); that is, the pH value is the control variable, and it is kept constant along the computation. The SGCMC code is a modification of the one previously developed by our group to compute conformational and ionization properties of linear polyelectrolytes. 52,57 In particular, it has been extended in order to include the effect of mechanical work. As a result, bending and stretching of the bonds have also been implemented. The resulting program is rather general since it allows working with sites and bonds of different pK values, interaction energies, conformational energies, and so on. Excluded volume effects can also be taken into account. Moreover, the code can deal with any arbitrary distribution of the sites along the chain, which is chosen by the user. However, in this paper, we restrict its use to the assumptions detailed previously. A snapshot of one of the SGCMC realizations is shown in Figure 1.
The Metropolis algorithm 68,69 generates new states at constant pH in a chain with N = 50 ionizable sites (i.e., 148 nodes or M = 147 bonds), a number which is large enough to avoid end effects and ensures the reproducibility of the intensive properties of the polymer, such as bond state probabilities or degree of protonation. An outline of the algorithm is depicted in Figure 2. In each new MC configuration, the polyelectrolyte can change either (A) the rotational state of a bond, (B) the length or bending angle of a bond, (C) the ionization state of a binding site, or (D) the spatial orientation of the polyelectrolyte chain in the laboratory reference frame, with trial probabilities of 0.88, 0.1, 0.01, and 0.01, respectively. These values allow us to obtain a good equilibration of the conformational structure for a given ionization state so that the system does not get trapped in local minima. Each change in the rotational state of a bond j implies a ±120°rotation of its dihedral angle and the recalculation of distances among the sites situated before and after the rotating bond. The changes in the stretching and bending states of bond j are Δl j = ± 0.01 Å and Δα j = ± 0.5°, respectively. These variations provide an average acceptance ratio of ∼20%, which is an acceptable value to make proper statistics. The global spatial orientation of the polymeric chain is altered by changing the polar angle θ and the azimuth angle ϕ of the first and second bonds of the polyelectrolyte with respect to the laboratory coordinate frame by amounts Δθ = ± 2°and Δcos (ϕ) = ± 0.1, respectively. The latter change is important to avoid preferred orientations in the space at zero force. Once the free energy difference (ΔF of eq 1) between trial and current conformations is calculated, the new configuration is always accepted if 0 Δ < and accepted with a probability exp( The results presented represent the average over eight different SGCMC simulations. Each simulation has been equilibrated in the first 5 × 10 7 configurations, and the thermal averages have been computed in the following 4.5 × 10 8 realizations. The simulations were performed using a parallel code developed in C++ on a 126 CPU cluster. For each pH, ionic strength, and force, typical jobs were run using eight CPUs for 1 to 2 h. The chosen parameters are pK = 9 and u t = 10, similar to those corresponding to LPEI. Note, however, that the reduced free energy only depends on the reduced chemical potential μ = pH − pK. This means that the results and conclusions taken from the simulations are the same for any pK value by choosing a suitable pH value for which the difference pH − pK is the same. The simulations are performed at room temperature T = 298.15 K. Free protons, co-ions, and counterions are not explicit in the simulations, and the screening effects are taken into account via the Debye length parameter, κ −1 , in eq 11. The chosen values for the parameters in the stretching and bending potentials are l 0 = 1.5 Å, α 0 = 120°, k length = 300 kcal mol −1 Å −2 , and k angle = 0.01 kcal mol −1 deg −2 , which are typical values used in molecular dynamics force fields for C−C bonds. 70 The average degree of protonation (θ) is computed as where ⟨N + ⟩ is the thermal average number of protonated sites. Note that since the simulations are performed at constant pH, N + is a fluctuating quantity, different in each new accepted configuration. Another interesting quantity is the effective protonation constant (K c ), which provides information about the average affinity of the macromolecular sites for the protons. 62,71,72 In general, K c depends on the charge of the macromolecule, different at each pH value. It is defined as The probability of a rotating c bond to be in the gauche state, P(g), is calculated as where ⟨N g ⟩ is the thermal average number of rotating bonds in a gauche state. Other quantities, such as the probability of having two neighboring c bonds in given conformations, (e.g., tt, tg + , g + g + , etc.) can be calculated in the same way. The extension of the polyelectrolyte chain (L z ) in the direction of the mechanical force, that is, the z axis, is obtained as where z i is the z coordinate of site i in the laboratory coordinate frame. Finally, a very useful quantity to understand where r i is the position of site i. As a consequence of eq 17, the Kuhn length can be directly related to the persistence length by l K = 2l P − l 0 .

■ RESULTS AND DISCUSSION
In this section, we will discuss the effect of the pH value and the ionic strength on the force−extension curves by simultaneously analyzing the dependence of conformational (chain elongation, bond state probabilities, and persistence length) and protonation properties (degree of protonation and effective protonation constant). As commented above, the microscopic pK value of the protonating sites will be fixed at pK = 9 throughout this section without loss of generality. Concerning the energy of the gauche state in the absence of charge, two cases physically relevant are considered. In the first case, the polyelectrolyte exhibits free rotation (i.e., gauche and trans states have the same energy, ε σ = 0, σ = 1). In the second case, the gauche states of the c bonds are favored, for instance, because of the existence of hydrogen bonding between two consecutive protonating sites, which means that σ > 1 (we take ε σ = 1, σ = 10). This phenomenon has been observed in LPEI and POE. 57 The case σ < 1 is not much interesting since most of the bonds are in the trans state so the chain is basically extended even in the absence of force. Finally, the interaction energy between two charged neighboring sites through a trans c bond is fixed at ε u, t = 1 (u t = 0.1) for all the studied cases, which is the order of magnitude found in a number of weak polyelectrolytes by using potentiometry. These works indicate that, for the same molecule, neighboring interactions are rather independent of the ionic strength. 56,57,69 Since charge regulation is a key ingredient of the model, let us first analyze the behavior of the degree of protonation θ when no mechanical force is applied, which will be useful in the foregoing discussion. In Figure 3, the titration curves for the cases σ = 1 (a) and σ = 10 (b) are shown at four different ionic strengths (I): 1, 0.1, 0.01, and 0.001 M, from top to bottom. It can be observed that in both cases, lowering the ionic strength results in a decrease of the degree of protonation for all of the pH values, which is explained by the increase in the LR electrostatic repulsion. Note that this effect is larger in the case σ = 10, which can be explained by the fact that gauche states are energetically favored, which hinders the possibility of having two neighboring sites charged (since u g = 0). The effective protonation constant K c is depicted in Figure 3c,d for σ = 1 and σ = 10, respectively. Clearly, K c presents two asymptotic values. At high pH, the macromolecule is not charged, electrostatic interactions are absent, and the K c value corresponds to the microscopic pK value of the ionizable sites (pK = 9). However, as pH decreases, sites get ionized and the work needed to protonate an empty site increases due to electrostatic repulsion. This results in a decrease in K c , which, at low enough pH values, reaches a new asymptotic value. This decrease in the affinity for the proton is especially relevant at low ionic strengths for which the LR interactions are stronger since screening is weaker. Finally, note that, for the same pH and I values, the decrease in pK c in the case σ = 10 is more pronounced than in the case σ = 1. For σ = 10, the gauche states are energetically promoted, the chain is more folded, and the distance between charger sites is shorter, which leads to larger LR interactions.
Effect of the pH Value on the Force−Extension Curves. The force−extension curves are shown in Figure 4a,b for the cases σ = 1 and σ = 10, respectively, for pH values ranging from 4 to 10 (top to bottom). The chain extension is normalized to the contour length  Figure 4g,h. The images in the left side always correspond to the case σ = 1, while those in the right side refer to the case σ = 10. It is important to stress that stretching, conformational, and ionization properties are highly coupled, so we will discuss them all at once. Owing to the lack of space, in the main document of this work, we just present and discuss the force−extension curves, which are more relevant for the purpose of this study. However, for the lecturers interested, the full casuistry, covering the complete range of pH values and ionic strengths, is reported in the Supporting Information. Let us first analyze the different force regimes in terms of the progressive activation of the degrees of freedom of the bonds, that is, rotational, bending, and stretching. First of all, note that, for low enough forces and for all the pH values, both the gauche state probability P(g) and the normalized persistence length l P/ l 0 remain constant. Despite this fact, the chain is considerably extended, around 10% for pH = 10 and a remarkable 50% for pH = 4. This indicates the existence of a low-force regime (corresponding to F < F E = k B T/l P ≈ 1 pN) 20 under which the chain behaves as a structureless set of segments with characteristic length l P .
The dependence of P(g) and l P on the pH value is shown in Figure 4. As observed in Figure 4c,d, P(g) is strongly affected by the pH value at low forces. This fact can be better explained by observing the behavior of the average degree of protonation θ (Figure 4g,h) for F < F E . For σ = 1, θ increases from θ ≈ 0.05 at pH = 10 to θ ≈ 0.7 at pH = 4, while for σ = 10, θ increases from θ ≈ 0.05 at pH = 10 to θ ≈ 0.6 at pH = 4. In both cases, when the adjacent sites are simultaneously protonated, the electrostatic repulsion is so strong that the gauche states are forbidden (u g = 0). As a result, the rotational properties change with the pH value, resulting in two limit behaviors. At high pH, when the polymer is discharged, the bonds present free rotation (when σ = 1) or preference for the gauche state (when σ = 10). On the contrary, at low enough pH values, when the macromolecule is almost completely charged, the restriction u g = 0 forces all the bonds to adopt the trans conformation. The increase of the number of bonds in the trans state due to lowering the pH value can be clearly observed in Figure 4c,d. As a result, the polymer chain gets stiffer, and the persistence length increases, as can be observed in Figure 4e,f. In turn, this fact leads to the increase in the chain elongation observed in Figure 4a,b, which is more marked for σ = 1 than for σ = 10. In the latter case, the gauche state is energetically favored so that a larger charge (i.e., a lower pH) is required to obtain the same number of trans bonds and to increase the persistence length.
The elongation versus force curves in the low force regime are shown in Figure 5 where the simulation values are depicted as markers. As can be observed, the low-force regime can be divided into two subregimes. For very low forces (F < 0.1 pN), the chain behaves as an entropic spring, and the elongation responds linearly to the applied force (dashed lines in the figure). The relation between elongation and force is given by This expression is independent of the structure of the chain since it comes directly from the fluctuation−dissipation theorem. 73 Under this subregime, the mechanical work is much smaller than the thermal energy. However, for larger forces (0.1 < F < 1 pN), the situation becomes more complex. It is found that the extension follows the Pincus scaling law 73,74 (continuous lines) Note that ν is found to range from ν = 1/2 at high pH (uncharged chain and corresponding to the linear regime) to ν = 3/5 at low pH (at which the chain is almost fully charged). The latter value was first predicted by Pincus 74 for strong polyelectrolytes, and it can be explained as the effect of electrostatic excluded volume interactions and the corresponding swelling of the macromolecule. Interestingly, both limiting values are obtained with great accuracy from the simulations, which nicely confirms Pincus theory. For pH values ranging from 6 to 8, intermediate values of ν are obtained, indicating that a partially charged weak polyelectrolyte can be seen as an intermediate situation between the neutral chain and a strong polyelectrolyte. For the case σ = 10, the transition between the two limiting cases is more gradual than for σ = 1 due to the fact that lower pH values are necessary to fully charge the chain (see Figure 3).
So far, we have analyzed the low-force regime for which the persistence length remains constant with the applied force. For larger forces, however, the rotational degrees of freedom are activated. This fact makes P(g) decrease with the force. In this new situation, the bonds, which were initially in the gauche state, gradually adopt the trans state by effect of the pulling force. The stretching mechanism is no longer entropic, but it depends on the free energy of the trans/gauche transition, and the persistence length becomes force-dependent. This fact is in contrast with DNA, with a much more rigid structure and for which charge regulation can be neglected. 35 The new characteristic force F R for which the rotational degrees of freedom are activated can be roughly estimated by equating the mechanical work per monomer to the free energy of the bond state transition ΔF t → g so that The F R resulting values are 30 and 70 pN for σ = 1 and σ = 10, respectively, in agreement with the observed range of forces (1−100 pN) for which the variation of P(g) is more pronounced. Moreover, the conformational degrees of freedom are coupled with proton binding due to the term (eq 9) in the reduced free energy. It is observed that, in increasing the force value, the change in the rotational states from gauche to trans is simultaneous with the increase in the macromolecular charge. Two asymptotic behaviors are found again. At low forces, the protonation state is the same as the one of the non-stretched molecule. As commented above, θ is lower for σ = 10 than for σ = 1. At large enough forces, however, a new plateau arises, and θ is the same value for both σ values. The gap between the θ value at low and high forces is thus larger for σ = 10 and depends on the pH value. For instance, for σ = 10 and pH = 4, θ ranges from less than 0.58 in the linear regime to 0.83 for the larger forces. At large pH values, however, the change in θ is smaller in absolute terms although it can be significant in relative terms. For instance, for σ = 10 and pH = 8, θ ranges from 0.2 to 0.32, which means an increase of more than 50%. Using the definition (eq 13), the effective protonation constant K c can be calculated as a function of the force. The obtained effective pK value, which is not shown here but can be found in the Supporting Information, slightly increases with F for both σ values. However, this effect is very weak since the pK value at large forces exceeds the low force limit by, at most, 0.5 pK units. However, this variation seems to be enough to cause significant changes in the macromolecular charge in applying an external force. This point will be discussed in more detail in the next subsection for the full range of ionic strengths. Finally, when the force is large enough to deform bond angles and lengths, a third situation arises. The average bond length ⟨l⟩ (green squares) and bond angle ⟨α⟩ (black dots), normalized by their respective equilibrium values, are shown in Figure 6 as functions of the applied force. It can be observed that the bond length and angle only start to be significantly elongated for values larger than the characteristic force F S > 300 pN. A rough estimation of F S using the Hooke law also confirms this value. Note that the bond angle is slightly easier to deform than the bond length. Although it is not shown in the figure, bond stretching and bending have been found to be independent of the pH value and ionic strength. This fact explains why the force−extension curves in Figure 4a,b also become almost independent of the pH for forces larger than F S . At this point, most of the bonds are in the trans state. Since the force constants are independent of the ionization state, the response to the applied force is also the same for any degree of protonation. Finally, for forces around and larger than F S , the bending potential becomes anharmonic, finally leading to bond breaking, as reported in several AFM single-molecule stretching experiments. 1 We conclude that, for intermediate forces and suitable pH and ionic strength values, CR can be induced by an applied force when the mechanism of CR is proton binding. For other types of binding mechanisms, such as metal binding, the ionic charge and binding constants are much larger, and the binding mechanism strongly depends on the conformational state (for instance, because of the presence of chelate complexes). As a result, CR could be significantly enhanced. In those cases, which are out of the scope of the present work, the interplay between stretching and CR could be of technological interest.
Effect of the Ionic Strength. Let us investigate the effect of the ionic strength in the force−extension curves, which are represented in Figure 7a,b for ionic strengths (from bottom to top) of 1, 0.1, 0.01, and 0.001 M. Now, the pH value is fixed to pH = 6. For this pH value, the macromolecule is approximately half charged, so it is a suitable value to discuss the influence of  finally the large-force regime, corresponding to deformations in the bond angle and length. However, unlike the effect of the pH value, the effect of ionic strength on the conformational properties is more complicated due to protonation and the complex interplay between SR and LR interactions.
Let us first analyze the dependence of the binding properties on the applied force, depicted in Figure 8. In all the cases, θ increases with F for forces larger than the low-force regime F > F E ≈ 1 pN. As a general trend, the polyelectrolyte chain is on average more elongated as F increases so that the mean distance between sites increases and the LR electrostatic repulsion decreases, allowing more sites to be protonated. Concerning the SR interactions, note that P(g) experiences an important decrease in the interval F = 10−100 pN (see Figure  7c,d), which is especially dramatic in the case σ = 10. This fact implies a drastic change in the chemical environment of the ionizable sites, which become separated by trans bonds, through which the repulsion is much weaker. Charge regulation is clearly induced by the mechanical force. Interestingly, the larger the ionic strength, the more intense charge regulation is. Especially remarkable is the case I = 1 M and σ = 10 for which the charge is almost doubled at high forces. This indicates that the charging process is basically a local phenomenon, which is essentially driven by SR interactions and the conformational state of the c bonds and is rather independent of the ionic strength. Conversely, LR interactions, which increase on lowering the ionic strength, weaken charge regulation because they discharge the molecule in all the force regimes. In the same way, the effective pK value also increases with the stretching process for F > F E ≈ 1 pN, as can be observed in Figure 8c,d, so that a larger affinity for the protons is induced by the applied force. Again, this effect is especially relevant for σ = 10 and at high ionic strengths.
Concerning the dependence of the gauche probabilities on the ionic strength, note that, as observed in Figure 7c,d, P(g) seems to be in contradiction with the behavior of the persistence length. On increasing the ionic strength, P(g) decreases and so does the number of gauche bonds, and apparently, the chain should be stiffer. However, for forces below F = F R ≈ 20 pN, the persistence length also decreases, so actually the chain gets more folded. This effect can be observed for the two σ values although it is especially relevant for σ = 10. This apparent paradox can be explained by taking a look at the degree of protonation. As commented above, θ increases with the ionic strength. As a result, the probability of having two charged neighboring sites increases. Since they cannot be both protonated and linked by a bond in the gauche state at the same time (u g = 0), the number of gauche states decreases. Note that this effect is due to action of the SR interactions. However, unlike the gauche probability, the persistence length is a "global" property, the result of the simultaneous action of many bonds and sites. As a consequence, LR interactions play a more important role in the behavior of l P . The same argument can be used for the chain elongation, which decreases with the ionic strength in all the curves below for F < 20 pN.
For F > 20 pN, the elongation becomes almost independent of I, but on examining in detail the curves, we can observe an unexpected result: the extension slightly decreases with I. This intriguing trend of the elongation is consistent with the behavior of the persistence length: for F < F R , the chain gets stiffer on lowering the ionic strength, but, for larger forces, it gets slightly more flexible. For lower pH values, this effect is more visible (as in the curves for pH = 4 shown in the Supporting Information). This certainly small effect seems to be apparently irrelevant but again points out the non-trivial connection between SR and LR interactions. We suspect that this result is related to spatial correlations of the states of neighboring bonds although probably a deeper insight is necessary.
As an example of the formation of spatial patterns, the probabilities of having two consecutive c bonds in a given conformation P(c i c i + 1 ) as a function of F have been calculated at two different ionic strengths I = 1 M (Figure 9a,b) and I = 0.001 M (Figure 9c,d). Again, the images on the left side and on the right side correspond to the cases σ = 1 and σ = 10, respectively. Due to the polyelectrolyte symmetry, there are only four different combinations of bond states: both bonds in trans (with a probability P(tt)), one bond in trans and its neighbor in gauche (P(tg)), both bonds in gauche with the same orientation (P(g + g + ) + P(g − g − )), and both bonds with different orientations (P(g + g − )). As expected, P(tt) monotonically increases as the mechanical force increases and, for large enough forces, P(tt) tends to 1 independently of the value of σ. For σ = 1 and I = 0.001 M, the preferred combination is trans− gauche at low forces, but there is a crossover, which can be observed at F ≈ 50 pN, curiously at the same force interval at which the change of tendency in l P and in the elongation is observed. For σ = 10, the preferred combination at low forces is g ± g ± . However, the most intriguing point is the observed maximum in P(tg) at F ≈ 50 pN. This maximum implies the existence of an intermediate situation where the mechanical work contribution, which promotes the trans state, competes with the energetic stabilization of the gauche states due to the fact that σ > 1. In this force regime, the polyelectrolyte adopts a highly ordered structure, which alternates bonds in the trans state with bonds in the gauche state. Again, the presence of this maximum coincides with the force interval for which l P and the elongation switch their dependence on the ionic strength. We would like to highlight that, even for the very simple model of polyelectrolyte presented here, one finds a rather rich physical−chemical behavior, which includes charge regulation, complex conformational transitions, and highly correlated spatial structures.

■ CONCLUSIONS
In this work, the influence of charge regulation, highly coupled with the conformational degrees of freedom, in the stretching Figure 9. Probability of having two neighbor bonds in given conformations P(c i c i + 1 ) versus F at two different ionic strengths of (a, b) I = 1 M and (c, d) I = 0.001 M at constant pH = 6. The combinations shown are both bonds in trans P(tt) (blue squares), one bond in trans and the neighbor in gauche P(tg) (orange circles), both bonds in gauche with the same orientation P(g + g + ) + P(g − g − ) (purple upward triangles), and both bonds in gauche with the opposite orientation P(g + g − ) (green downward triangles). The images on the left side (a and c) correspond to σ = 1, and those on the right side (b and d) refer to σ = 10. The rest of the parameters have the same values as those in Figure 1b. properties of weak polyelectrolytes is studied. With this aim, we propose a model, which captures the fundamental aspects present in a flexible weak linear polyelectrolyte (internal angle rotation, bond stretching, bond bending, and proton binding) with a minimum number of parameters. The model was inspired by the structure of linear poly(ethylenimine) (LPEI), a symmetric weak polyelectrolyte with an ionizable site every three chain positions. It is based on the site binding rotational isomeric state (SBRIS) model, which allows studying conformational and ionization properties on the same foot. Short-range (SR) and long-range (LR) electrostatic interactions are treated separately. LR interactions are chemically unspecific and can be reasonably implemented using the meanfield Debye−Huckel potential. Conversely, SR interactions between neighboring sites and bonds are mediated by the macromolecular skeleton so that they depend on the detailed chemical environment of the site. As a result, specific energetic parameters are used to describe them. Bond stretching and bending are included by means of harmonic potentials. The resulting scheme is used to perform semi-grand canonical Monte Carlo (SGCMC) simulations at constant pH and applied force. Concerning the energy of the gauche state, two situations are studied, controlled by the energy of the gauche state (corresponding to the Boltzmann factor σ): when trans and gauche states have the same energy (σ = 1) and when the gauche state is energetically stabilized, for instance, by hydrogen bonding (we take σ = 10). The influence of the pH value and the ionic strength in the force−extension curves is analyzed. In order to understand the different mechanisms of chain stretching, the degree of protonation θ, bond state probabilities of the gauche state P(g), and persistence length l P as functions of the force are also analyzed.
As a general trend, three force regimes are found. In the lowforce regime, the persistence length is force-independent, and two subregimes are identified. Up to 0.1 pN linear behavior is found, as demanded by the fluctuation−dissipation theorem, for all the pH values. From 0.1 to 1 pN, however, the chain presents Pincus scaling behavior depending on the pH value. For high pH values (i.e., neutral chain), the elongation is still linear with Pincus scaling exponent ν = 1/2, while for low pH values (fully charged chain), ν = 3/5, the theoretically predicted value for strong polyelectrolytes. For intermediate pH values, ν has been found to present a gradual transition between the two limiting values. In the large-force regime, most of the bonds are in the trans state so that the stretching becomes approximately independent of the pH and the ionic strength. Finally, there is an intermediate regime, between 1 and 100 pN, for which the rotational and protonation degrees of freedom, which are highly coupled, are activated. This force regime is the most interesting one since conformational transitions, charge regulation, and spatial correlations are observed. It is in this force regime that the pH value and the ionic strength maximally influence the chain elongation. Mechanically induced charge regulation is mainly driven by SR interactions. When the macromolecule is elongated, the trans states are promoted so that the electrostatic interaction between neighboring sites decreases, favoring the affinity for the protons. This effect seems to be larger at large ionic strengths and pH values for which the molecule is partially charged.
The role of the pH value is relatively straightforward to understand. On lowering the pH value, the macromolecule gets charged, promoting trans states and larger distances between sites, thus reducing the electrostatic repulsion. This results in an increase in the persistence length, a reduction in the number of gauche states, and an easier extension of the chain. Therefore, a significant influence of the pH value on the curve−extension curves is found both for σ = 1 and σ = 10.
On the other hand, the effect of the ionic strength for a fixed pH value is more complicated since it depends on the complex interplay between SR phenomena (bond conformations and protonation) and the LR interactions. It is found that the exhibited tendency of P(g) seems to be in contradiction with the behavior of the persistence length. On increasing the ionic strength, P(g) decreases and so does the number of gauche bonds, and apparently, the chain should be stiffer. However, for forces below F = F R ≈ 20 pN, the persistence length also decreases, so the chain gets more folded. This apparent paradox can be solved by observing that the charge decreases on increasing the ionic strength while the intensity of LR interactions is enhanced. Ionization equilibria therefore play a fundamental role in the stretching properties of weak polyelectrolytes. Finally, it is found that in the intermediateforce regime, spatial correlations are formed, which also determine some subtle aspects of the stretching process.
We would like to highlight that, even for the very simple model of polyelectrolyte presented here, one finds a rather rich physical−chemical behavior, which includes charge regulation, complex conformational transitions, and highly correlated spatial structures. To our knowledge, this work is the first attempt to study, at least by means of computational simulation, the mechanical stretching of a weak polyelectrolyte including the coupling of charge regulation and conformational equilibria. Notes HORIZON 2020 (692146). P.M.B. also acknowledges the financial support from a grant (FI-2017) of Generalitat de Catalunya. We thank Professor Michal Borkovec for introducing us to the topic, helpful discussions, and suggestions.