Coupling of Charge Regulation and Conformational Equilibria in Linear Weak Polyelectrolytes : Treatment of Long Range Interactions via Effective Short-Ranged and pH-Dependent Interaction Parameters

The classical Rotational Isomeric State (RIS) model, originally proposed by Flory, has been used to rationalize a wide range of physicochemical properties of neutral polymers. However, many weak polyelectrolytes of interest are able to regulate their charge depending on the conformational state of the bonds. Recently, it has been shown that the RIS model can be coupled with the Site Binding (SB) model, for which the ionizable sites can adopt two states: protonated or deprotonated. The resulting combined scheme, the SBRIS model, allows to analyse ionization and conformational equilibria on the same foot. In the present work this approach is extended to include pH-dependent electrostatic Long Range (LR) interactions, ubiquitous in weak polyelectrolytes at moderate and low ionic strengths. With this aim the original LR interactions are taken into account by defining effective Short Range (SR) and pH-dependent parameters, such as effective microscopic protonation constants and rotational bond energies. The new parameters are systematically calculated using variational methods. The machinery of statistical mechanics for SR interactions, including the powerful and fast transfer matrix methods, can then be applied. The resulting technique, to which we will refer as Local Effective Interaction Parameters (LEIP) method, is illustrated with a minimal model of a flexible linear polyelectrolyte containing only one type of rotating bonds. LEIP reproduces very well the pH dependence of the degree of protonation and bond probabilities obtained by semi-grand canonical Monte Carlo simulations, where LR interactions are taken explicitly into account. The reduction in the computational time in several orders of magnitude suggests that the LEIP technique could be useful in a range of areas involving linear weak polyelectrolytes, allowing direct fitting of the relevant physical parameters to the experimental quantities.


Introduction
The ionization state of charged macromolecules in solution is regulated by the binding of small ions (protons, metal ions, etc.) present in the backward medium.In particular, acid-basic equilibria in weak polyelectrolytes represent the paradigmatic mechanism of charge regulation due to the ubiquitous presence of proton ions in aqueous solution.These processes are of paramount importance to understand the physicochemical behavior of charged macromolecules in a wide range of situations.Just to mention a few examples, charge regulation plays a fundamental role in receptor-ligand equilibria in biochemical systems [1][2][3][4], supramolecular chemistry [5][6][7], the role of natural organic matter in geochemical cicle of metal ions [8], wastewater treatment [9], stability of colloidal systems [10], advanced coating in material science [11][12][13], or drug delivery [14].Charge regulation can take place on rigid structures, such as surfaces or nano-particles [15], but in general polyelectrolytes are flexible and conformational and ionization degrees of freedom are strongly coupled.This fact can result in dramatic structural changes in the macromolecule.Classical examples are the helix-coil transitions of poly(peptides) [16], the swelling of poly(methacrylic) acid in a very narrow range of pH [17] or the strong influence of ionization in the folding of proteins [18].More recently, it has been recognized the importance of the ionization configuration in the conformational properties of intrinsically disordered proteins, whose function-structure relationship still remains a controversial matter [19,20].
The understanding of the ionization processes has been mainly based on the so-called Site Binding (SB) model.In this approach, the ionization configuration of the macromolecule is defined as a set of sites which can be in two possible states, i.e. protonated or deprotonated, as outlined in Fig. 1a.The free energy is then parametrized by site-specific microscopic protonation constants and interaction energies between sites.Triplet or higher-order interactions among sites can also be considered.Once the system Figure 1.a) Outline of the Site Binding (SB) model for a linear polyelectrolyte represented as a linear chain of ionizable sites.The ionization state of the macromolecule is characterized by a set of variables s = {s i }, i=1,... , N, which can adopt two possible values: s i = 1 if the site i is protonated (dark blue circles) and s i = 0 if it is deprotonated (cyan circles).b) Sketch of a linear chain joining ionizable sites by means of rotating bonds as an example of the Site Binding-Rotational Isomeric State (SBRIS) model.Both ionization and conformation degrees of freedom are now taken into account.In the depicted chain only the bonds type c are able to rotate, which can take three possible states of minimum energy, i.e., trans, gauche+ and gauche-.
In the case of linear polyelectrolytes, the transfer matrix method can be used to compute the relevant thermal averages [15].This powerful and elegant technique was originally designed to solve the classical Ising model of ferromagnets.It is based on the fact that the partition function of a system with N+1 sites can be related in a recursive way to the one of a system with N sites.The resulting recursive relationship can be expressed in terms of the transfer matrix, whose elements represent the contributions of the new site to the partition function, for a given state of the preceding site [38,39].The technique is very versatile and can be generalized to systems composed with repetitive units (spins, bonds or binding sites) which can take in principle more than two states.When applied to the binding of ions to polyelectrolytes, the method can be adapted to include a wide range of phenomena such as triplet interactions between sites [21], chelate complexation of metal ions [23], proton binding to polyampholytes [40,41], protein-DNA binding [42], super-capacitator charging [43] or coupling between ionization and conformational degrees of freedom [44][45][46].
Probably the most productive application of transfer matrices was proposed by Flory in the context of the Rotational Isomeric State (RIS) model [47,48], aiming to compute conformational properties of neutral linear molecules.RIS model relies on the observation that, although a particular bond can adopt in principle any rotation angle, only those of minimum energy (typically trans, gauche+ and gauche-) are significantly populated.As a consequence, each bond can be regarded as a 'unit' of the system adopting three possible states.The corresponding partition function and the necessary thermal averages (bond probabilities, end-to-end distance, radius of gyration, etc.) can be calculated using proper product of transfer matrices.In recent works [45,46], it has been shown that SB and RIS models can be combined in a unique scheme so that conformational and ionization equilibria can be analyzed on the same foot.It has been shown that all the matricial expressions of RIS can be systematically extended to account for the ionization degrees of freedom.The resulting SBRIS model, outlined in Fig. 1b, has been recently applied to the detailed characterization of the conformational and ionization properties of linear poly(ethylene)imine [46].
The main limitation of the transfer matrices used in SB, RIS and SBRIS models is that they can only deal with Short Range (SR) interactions [49,50].SR interactions are chemically specific and can produce important correlations between neighbouring sites and bonds.They cannot be modeled by simple continuous force fields (such as van der Waals or Debye-Hückel potentials) [51] but, in exchange, they can be easily implemented in a transfer matrix scheme.For polyelectrolytes, however, this is an important restriction, due to the Long Range (LR) nature of coulombic interactions, which severely restricts the range of application of the transfer matrix approach.In practice, the possibility of neglecting LR coulombic interactions must be restricted to high ionic strengths, an important limitation specially for polyelectrolytes which become insoluble under such conditions [52][53][54][55].
In a recent paper [56], the SB model has been extended to include LR interactions by introducing a modified free energy involving Local Effective Interaction Parameters (LEIP) which account for the LR interactions in an effective way.The LR force field is thus replaced by a short ranged effective one.The new local effective parameters, i.e. effective protonation free energies, effective pair interactions and so on, can be systematically calculated by using the Gibbs-Bogoliubov variational principle [39].The resulting modified free energy converges very fast to the exact free energy.It was found that the correction to the site protonation pK (first order correction) is enough to obtain an excellent, exact from the practical point of view, agreement between theory and MC simulations.This previous study, however, was restricted to rigid molecules, and conformational degrees of freedom were not taken into account.The main goal of the present work is to extend the LEIP method to account for the coupling between charge regulation and conformational equilibria involving LR interactions.In addition to allow much faster computations of ionization/conformational properties (computational times are reduced in orders of magnitude), the methodology here presented adds new physical insight in the interplay of conformational and ionization degrees of freedom in polymeric structures.For instance, the energy of the gauche state of a bond will now depend on the pH and the ionic strength, even if such a bond does not hold any ionizable group.The use of the LEIP technique in the SB model is reviewed in section 2. In section 3, the technique is generalized in order to include conformational equilibria coupled to LR coulombic interactions represented by the Debye-Hückel potential.In section 4 the theoretical results are compared to MC simulations.The new ideas here introduced are illustrated with a minimal model of a flexible linear weak polyelectrolyte containing only one type of rotating bonds.

Simultaneous treatment of short and long range interactions in rigid molecules
The ionization state of a macromolecule with N ionizable sites can be characterized by a set of variables s = {s i }, i=1,... , N, which can adopt two possible values: s i = 1 if the site i is protonated, and s i = 0 if it is deprotonated.The corresponding reduced free energy can be expressed in terms of the variables s i by means of the so-called cluster expansion [24] H (s) ln 10 where ) is the reduced chemical potential, which depends on the proton activity, a H , and the protonation pK-value of the ionizable site i, pK i ; φ ij represents the interaction energy of the sites i and j; τ ijk accounts for possible triplet interactions among sites i, j and k, and so on.The term "reduced" refers to the fact that the chemical potential incorporates both the pH and the protonation pK, which simplifies the subsequent expressions.The interaction (or cluster) parameters are expressed in thermal units, i.e. β = 1/k B T = 1, and divided by a factor ln 10 in order to be compared in the pH scale.Note that the conformation degrees of freedom are omitted in 1, so that the interaction parameters should be understood as proper averages over the conformational states.The mathematical form of these averages is not trivial and expressions for them are given in [45].Throughout this work we will assume that a site is charged when it is protonated, i.e. we are dealing with poly-cations.However, the subsequent arguments are also applicable to poly-anions with a suitable change in the protonation variables [15].The expansion of the free energy (1) usually converges very fast to the exact free energy, and for most of cases, the inclusion of triplet interactions is enough to accurately reproduce the measurable quantities, such as the degree of ionization of the individual sites [22].These can be obtained from H (s) by means of the semi-grand canonical partition function The average degree of protonation of a particular site i is related to Ξ as where Ω = − ln Ξ is the thermodynamic potential associated to the semi-grand canonical ensemble.The average number of bound protons is given by The correlation of the protonation degrees of two sites i and j, a quantity which will be used later, can be expressed as As can be seen, the quantification of all the relevant physical quantities relies in the accurate determination of the partition function Ξ.If the number of sites is small (N ≤ 20), Ξ can be evaluated by direct enumeration of all the possible ionization states.Otherwise, Monte Carlo (MC) simulations must be performed.In some cases, however, methods borrowed from Statistical Mechanics can be used.Among them, probably the most elegant one is the transfer matrix method, consisting of relating the partition function for a system with N+1 sites with that with N sites in a recursive way.This method was firstly used in the exact solution of the Ising model of ferromagnets [38,39].The link between both partition functions is the transfer matrix whose elements are the Boltzmann factors corresponding to the increase in the reduced free energy.For instance, for the linear polyelectrolyte sketched in Fig. 1, and assuming only nearest neighbour interactions, the partition function can be expressed as [47] where T is the transfer matrix z = Ka H represents the reduced activity and = − log u is the interaction free energy between neighbouring sites.q = (1, 0) and p = (1, 1) are the initiating and terminating vectors.This would be the simplest use of the transfer matrix.
The main limitation of the transfer matrix methods is that they can be only used when Long Range (LR) interactions are neglected, since the size of the transfer matrices grows exponentially with the range of the interactions [50].This is an important limitation of the method when dealing with polyelectrolytes, since it can be only used at high enough ionic strengths, for which the screening is enough to avoid the LR interactions.In a recent paper [56] we introduced a method which allows to include the LR interactions in a very accurate way.In this approach, the full free energy (1) is replaced by a new one involving only short range (SR) interaction parameters, accounting for the LR interactions in an effective way.The resulting formalism deals with both SR and LR interactions simultaneously.The method can be used for any kind of molecular or surface geometry, but it is restricted to rigid structures, so that the conformational degrees of freedom are not explicitly taken into account.Since the main goal of this work is to extend this formalism to flexible molecules and polyelectrolytes, we briefly outline the basic ideas of the method.The details of the derivations are given in Ref. [56].Although the following arguments can be readily generalized to the general form of the free energy (1), let us consider the simplest case of a rigid linear chain with identical sites, such as the one sketched in Fig. 1a.For this system, µ 1 = µ 2 = ... = µ and triplet interactions are omitted, i.e. τ ijk = 0.The reduced free energy H can be split into two contributions H = H 0 (x) + ∆H (x) such as where x is a parameter to be determined.Note that H 0 corresponds to a reduced free energy containing only nearest neighbour interactions of energy = φ i,i+1 , which can be exactly solved by using the transfer matrix (7).Now we can use the Gibbs-Bogoliubov variational principle [39] Ω ≤ Ω = Ω 0 (x) + ∆H (x) 0 ( 9) to determine the optimal value of x, where Ω 0 (x) = − ln Ξ 0 and • • • 0 represent the free energy and the thermal average corresponding to H 0 , respectively.Minimizing Ω with respect to x, it is found that x fulfills the equation [56] x = dϕ 0 /dx dν 0 /dx = dϕ 0 dν 0 (10) where is the LR energy averaged over the unperturbed free energy H 0 , whose correlation function h 0 ij , can be exactly evaluated using (5).If the optimal value for x is used in the computations, the variational principle (9) implies that all the thermal averages (degree of protonation, correlation functions, etc.) can be obtained replacing the average • • • by • • • 0 , which can be exactly determined since only SR interactions are involved.Eqn.(10) provides a transparent physical interpretation of x: it is the average change in the LR interaction energy when a new proton is bound to the molecule at a given pH-value.As expected, x vanishes in absence of LR interactions and the nearest neighbour interaction model becomes exact.Therefore, x can be interpreted as the necessary correction to the reduced chemical potential µ in order to account for the LR interactions but in a local effective way.By the definition of µ = pH − pK, x can also be understood as the correction to the site pK-value, so that pK eff = pK − x is the effective pK-value, and it represents the extra energetic cost of the site protonation due to the presence of LR interactions.We will refer to this procedure as the Local Effective Interaction Parameters (LEIP) method.It is important to highlight that LEIP, unlike other approaches involving some mean-field approximation (such as the Bragg-Williams approximation in Ising models), includes the correlations via Eqn.11, although in an approximate way.This approximation, however, results to be extremely accurate, as can be observed in Fig. 2a, where the titration curves corresponding to a rigid linear chain with identical sites are depicted.The chosen parameters are pK = 9 and = 1.5.In this model, the LR interactions between distant sites are described by the Debye-Hückel potential where B 0.7 nm is the Bjerrum length in water at 298 K , d ij is the distance between the sites i and j, and κ −1 (nm) = 0.304/ I (M) is the Debye length at the ionic strength I.For a rigid linear chain, as the one shown in Fig 1, , where b is the separation between consecutive protonating sites.
We have plotted the titration curves obtained by Monte Carlo (MC) simulations in the semi-grand canonical ensemble, i.e. at constant pH (blue circles), together with the ones calculated using Eqns.10 and 11 (continuous lines) for all the range of ionic strengths and b=2 nm.Surprisingly, simulated and calculated curves overlap, so that, for this model, the LEIP solution can be regarded as exact from the practical point of view.The computational cost of LEIP methods is orders of magnitude much lower than that required in MC simulations, allowing the fitting of parameters to experimental titration curves.The correction to the pK, x, shown in Fig. 2b, increases in lowering the pH (i.e.increasing the charge), and in decreasing the ionic strength (lower electrostatic screening), since the energetic cost to protonate a site increases with the macromolecular charge and with the intensity of the LR interactions.Another advantage of the LEIP method is that it can be systematically improved by selectively correcting other cluster parameters.For instance one could decide to correct, not only the pK-value (pK → pK − x), but also the nearest neighbour interaction energy ( → − x ).Proceeding in the same way, it can be shown that x and x fulfill the non-linear system of equations [56]  where ν 0 (x, x ) and D 0 = s i s i+1 0 = h 0 12 (x, x ) represent the average number of protons and the average number of nearest neighbour interactions respectively, which can be exactly calculated using H 0 .Solving Eqns.(13) the correction to the pK and are obtained as functions of the pH.The physical meaning of x and x becomes transparent if the Eqns.(13) are rewritten in terms of ν 0 and D 0 as independent variables.After some elementary algebra, x and x adopt the much simpler form Eqns. ( 14) stay that x represents the increase in ϕ 0 for a constant number of interactions D 0 , while x can be interpreted as the change in ϕ 0 in creating a nearest neighbour interaction, keeping constant the number of bound protons ν 0 .Intuitively, one can guess that x is much smaller than x, so that the correction to the pK is enough to reproduce almost exactly the exact free energy, generating physical properties almost indistinguishable from the MC simulations.In Fig. 2a the titration curves have been recalculated using the correction to .As expected, no significant improvement is obtained.x and x as functions of the pH are shown in Fig. 2c, where it is clearly observed that x is much lower than x for all the ionic strengths.Note that the wavy behaviour of x in Fig. 2b is no longer present in Fig. 2c, and seems to be replaced by the contribution x .Using the same procedure, corrections to higher order interactions, such as triplet or next-nearest neighbour interactions, can be calculated until the desired accuracy is obtained, and expressions of the type ( 14) can straightforwardly generalized.The same treatment leads to very good results for heterogeneous polyelectrolytes and polyampholytes, by correcting the pK-values of the different kind of sites (pK i → pK i − x i ) [56].

Coupling of ionization and conformational equilibria
For a linear macromolecule composed by M bonds, a particular conformational state is denoted by a set of variables c = {c α }, j = 1,... , M. The variables c α can adopt several values corresponding to the rotational angles of the bond α.The possible states of the bonds are usually chosen as those of minimum energy, three in the simplest situation: trans, gauche+ and gauche-The selection of a finite number of rotational states instead of working with the full continuous rotational potential greatly simplifies the statistical mechanics treatment, and constitutes the basis of the Rotational Isomeric State (RIS) model, mainly developed by Flory [47].In the case of linear polymers, the transfer method can be used to determine the conformational partition function Ξ rot , which can be expressed as where U ff is the transfer matrix corresponding to the bond α.For a symmetric chain, for which the states gauche+ and gauchehave the same energy, and identical bonds, the transfer matrices are of the form where σ, ψ and ω are the Boltzmann factors associated to the conformational energies of the bonds: −k B T ln σ is the free energy of the gauche states while ψ (ω) are related to the interaction energies between two consecutive gauche states of different (same) orientation.ψ and ω equate one if the rotation of the bonds is independent.q = (1, 0, 0) and p = (1, 1, 1) are the initiating and terminating vectors.As in the SB model, the necessary thermal averages can be obtained by performing proper derivatives of the partition function.For instance, the average number of bonds in the gauche state is given by [47] RIS model can be generalized in order to take into account the protonation degrees of freedom.If the macromolecule is in a protonation state s, the pair (s, c) defines a roto-microstate with reduced free energy F rot (c) is the free energy corresponding to the fully deprotonated state of each conformer, while F p (s, c) represents the reduced free energy due to the protonation process which, for a given conformation, can be expressed as βF p (s, c) ln 10 where triplet interactions have been neglected.Note that the cluster parameters now depend on the conformational state c.The reduced free energy (18)(19) combines the RIS and the SB model and defines the SBRIS model, which allows to study conformational and ionization properties on the same foot.
In recent publicacions, the SBRIS model has been used to explain conformational transitions in weak linear polyelectrolytes [45] and in the characterization of ionization/conformational properties of linear poly(ethylenimine) [46].The probability of a specified roto-microstate is given by where the SBRIS partition function Ξ SBRIS is defined as The SBRIS partition function can alternatively be expressed in the fashion where Ξ rot (s) denotes the rotational partition function for the macromolecule in a 'frozen' binding configuration s = {s 1 , s 2 , • • • , s N }.Ξ rot (s) can then be calculated as a RIS partition function in terms of ( 16), but now decorating the transfer matrices with the suitable binding parameters.The sum over the protonation states can be performed by using proper matricial methods described elsewhere [46].They are outlined as supplementary information and here we just comment the final results.The SBRIS partition function is obtained by replacing the conformational RIS transfer matrices U (Eqn. 16) for suitable super-matrices.The rule is that, if a bond is holding at its ends two ionization groups, U must be replaced by where u is a diagonal matrix containing the Boltzmann factors corresponding to the short range interactions In (24) −k B T ln u t and −k B T ln u g represent the short range interaction energy between two protonated sites separated by a bond in trans and gauche conformation, respectively.For the bonds which do not hold ionization sites, the substitution is The resulting SBRIS partition function reads where r = (q q) and t = (p p) are the initiating and terminating vectors, respectively.The average number of bound protons and the bond state probabilities can be again obtained by proper derivatives of (26).It can be also shown that matricial expressions for other physical quantities derived in the context of the RIS model can also be generalized to ionizable molecules by performing suitable substitutions by super-matrices [45].For instance, there are available matricial expressions for the average square distance between two sites of the chain.These expressions are used in this work to estimate the average distance between charged sites and the corresponding LR interaction energy.As commented in the preceding section, transfer matrix methods can only be applied if only SR interactions are taken into account.In this work we propose to use the LEIP technique to include the LR interactions via local parameters, as done in the case of rigid molecules.Now, however, not only the ionization parameters, such as pK → pK − x, but also the conformational parameters must be corrected as outlined in Fig. 3.In the simplest case, with only one kind of rotating bonds, the substitution pσ → pσ − x σ where pσ = − log σ will be necessary.The treatment is almost identical to  the one used for rigid molecules.Now the "unperturbed" free energy is Ω 0 = − ln Ξ SBRIS (x, x σ ) .It can be easily shown that the corrections x for the pK and x σ for pσ fulfill equivalent equations to (13) where ν 0 (x, x σ ) represents the average number of bound protons (Eqn. 4)and g 0 (x, x σ ) the average number of bonds in the gauche state (Eqn.17), calculated using the unperturbed free energy.In the next section the case in which the LR interactions are described by the Debye-Hückel potential is discussed, although the method could in principle be applied to other kind of interactions such as van der Waals interactions.For flexible molecules, however, the average LR interaction energy ϕ 0 for the unperturbed free energy can only be calculated approximately.In this work, as a first approximation, we assume that This approximation could in principle be improved by using higher order moments of d ij .Matricial expressions for d 2 ij 0 and higher moments where derived by Flory and Jernigan [47,57] for neutral chains.Here these expressions are modified in order to account for the protonation degrees of freedom.An outline of the derivations is provided as supplementary information.

Results and Discussion
In order to estimate the accuracy of the LEIP method when applied to flexible polyelectrolytes, we compare the theoretical results with those obtained from Monte Carlo simulations in the semi-grand canonical ensemble.As a model system, we use the linear polyelectrolyte outlined in Fig. 1b, with protonating sites situated every three chain positions.Only the bonds between non-ionisable sites (c bonds) are allowed to rotate and they can adopt the three states of minimum energy, i.e. trans, gauche+ or gauche-.The conformation of c bonds determines the intensity of the SR interactions between neighbouring protonated sites.The rest of bonds (a and b in the figure) are forced to be in the trans state.The energy of the gauche state of the c bonds is denoted as pσ = − log σ.Two consecutive c bonds are assumed to rotate independently when the macromolecule is uncharged (ω = ψ = 1 in Eqn.16) .The protonating sites are considered to be identical with the same protonation pK-value.In the computations the used values of the bond length and the bond angle are 0.2nm and of 120 o , respectively.In the MC simulations, the free energy (F total ) is divided into a SR and a LR energetic terms The SR term involves the reduced chemical potential µ, the conformational energy of the bonds pσ, and the interactions between two consecutive protonated sites.The latter are characterized by the energies t = − log u t (when the bond c is in trans state) and g = − log u g (when the bond c is in gauche state).The LR contribution is explicitly calculated using the Debye-Hückel potential (Eq.12).This model can be regarded as a minimal model of a flexible polyelectrolyte, with only four energetic parameters involved (σ, t , g and pK), and it is here used to illustrate the application of the LEIP method to the analysis of the interplay of conformational and protonation degrees of freedom.
A Metropolis algorithm [15,27] is used to generate roto-microstates at constant pH with probability given by Eq. 20 in a chain with N=50.In each MC new configuration, the polyelectrolyte can change either (i) the conformational state of a rotating bond c or (ii) the ionization state of a binding site, with trial probabilities of 0.999 and 0.001 respectively.In doing so, good equilibration of the conformational structure for each new binding configuration is achieved and the system does not become trapped in local minima.Conformational changes of the bonds c imply a ±120 o rotation of its dihedral angle and a recalculation of the distance between all ionizable sites is performed.The results presented represent the average over 10 different MC simulations.Each MC 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.
Let us firstly consider the case for which c bonds can freely rotate when the adjacent sites are deprotonated (i.e.σ = 1) .When both sites are charged, however, the very strong SR repulsion hinders  2 with pK = 9 , t = 1, u g = 0 (i.e.g → ∞) and σ = 1 (a) ; σ = 10 (b).The chosen ionic strengths are, from top to bottom, 1M, 0.1M, 0.01M and 0.001M.Black continuous lines represent calculations using the LEIP method in which two effective short range parameters (pK and pσ) has been corrected.Red circles depict the MC values.The green dashed line corresponds to the LEIP values for I = 0.001M when only the pK-value is corrected, while pσ is kept constant.the gauche conformation, so that we take u g = 0 (i.e.g → ∞).The resulting titration curves are shown in Fig. 4a for ionic strengths ranging from 1M to 0.001M.The chosen parameters are pK = 9 and t = 1.The black continuous lines represent the average protonation degree θ calculated using the LEIP method correcting both the pK-value (pK → pK − x) and the conformational energy of c bonds (pσ → pσ − x σ ), while the red circles represent the results of the MC simulations.It is observed that the LEIP method reproduces very accurately the MC simulations for all the range of pH-values and ionic strengths.The dashed line depicts the values provided by LEIP for I=0.001M if only the pK-value is corrected, while σ remains constant.Although relative good prediction of the MC simulations is obtained, the quality of the titration curve clearly improves if σ is corrected.This means that the rotational energy of the c bonds is affected by LR interactions even if their pendant sites are not charged, as a result of the tendency of the chain to separate the rest of charged groups.Actually, the system behave as if c bonds "feel" the LR interactions in an effective way.This effect is more remarkable in the subsequent case.
Let us check the accuracy of LEIP method where the gauche states of the c bonds are favored, for instance, because of the existence of hydrogen bond, which means that σ > 1.When the adjacent sites are both protonated, however, the electrostatic repulsion is so strong that the gauche state states forbidden (u g = 0).In this case the rotational properties are opposite when the ionization states of the sites change.Fig. 4b compares the titration curves obtained using LEIP correcting pK and pσ and MC simulations for pK = 9, t = 1 and σ = 10.As can be observed, LEIP method and MC simulations yield to almost identical titration curves.In this case, however, the correction of pσ becomes compulsory.If only pK is corrected, the titration curve obtained by LEIP at 0.001M (green dashed line) exhibits a phase transition-like behaviour at pH ≈ 5.This is an artifact resulting from the impossibility to explain the complex interplay of charge regulation and conformational transition without taking into account the influence of LR interactions in the effective energy of the gauche state.
For the two cases commented above, the gauche state probabilities versus the pH are shown in Fig. 5a and Fig. 5b.Markers correspond to MC simulations while black lines represent the theoretical values at ionic strengths 1M, 0.01M and 0.001M, for pK = 9 , t = 1 and u g = 0. Fig. 5a corresponds to the case with σ = 1.A good correspondence between simulated and theoretical profiles is obtained for all the ionic strengths.Since at low pH-values the polyelectrolyte is almost fully protonated, the gauche state probability tends to zero because of the high electrostatic repulsion between nearest charged sites in the gauche position (u g = 0).On the other hand, at high pH-values the macromolecule is completely uncharged and the c bonds are freely to rotate.As a result, the probability of the two gauche conformers tends to 2/3.For I=1M, the LR interactions can be neglected and total correspondence between simulated and calculated values is found.At higher ionic strengths (0.01M and 0.001M), for  which the Debye-Hückel potential is not screened enough, some small differences arise.However, still now, good agreement between the LEIP method and MC simulations is observed.Fig. 5b corresponds to the case with σ = 10.Simulated and theoretical probabilities are also in good agreement.Again, the gauche state probability tends to zero for low pH-values, while at high pH-values the population of the gauche conformer is 2σ/(2σ + 1) = 0.95.A continuous transition from gauche to trans conformations as the pH decreases is observed.This transition is sharper than in the previous case.Note that, from the LEIP point of view, this transition occurs because of a double effect.On the one hand, there is the charging process, so that two adjacent sites tend to minimize the repulsion when the bonds adopt the trans conformation.This effect is present even when LR interactions are not present.On the other hand, the effective gauche state energy is increasing due to the average effect of the LR interactions (pσ → pσ − x σ ).Both effects are important to correctly reproduce the MC simulations.Otherwise, the lack of flexibility in the conformational energy leads to the spurious phase transition observed in Fig. 4b (green dashed line).
Fig. 6 shows the LEIP method correction x σ to the bond conformational energy (pσ → pσ − x σ ) for σ = 1 (Fig. 6a) and σ = 10 (Fig. 6b).In both cases, it is observed that x σ tends to zero at high pH-values, since the molecule is uncharged and no LR interactions are present, so no correction is necessary.As a general tendency, x σ tends to increase as the pH decreases due to the charging process and the corresponding increase in the LR interactions.This effect is larger at low values of the ionic strength since the Debye-Hückel potential is less screened.For the case σ = 10 a wavy behaviour is observed for ionic strengths 0.1M and 0.01M and x σ exhibits a smooth maximum at pH 4, which coincides with the pH regime where the trans to gauche transition is sharper.This fact could be due to correlations between the rotation of neighbouring bonds or because part of the correction is effectively include in the pK-correction x.Further analysis would probably be necessary in order to clarify this point.
In the two cases discussed above, we have taken u g = 0, which means that bonds between two protonated sites cannot be in the gauche state.Let us now relax this condition and take a finite value for u g , so that the electrostatic interaction between two charged sites in gauche is not forbidden but only penalized.LEIP predictions (black lines) and MC simulations (red markers) are plotted in Fig. 7. Fig. 7a shows the computed titration curves with g = 2 at ionic strengths ranging, from top to bottom, from 1M to 0.001M.Excellent agreement between the theoretical predictions and simulations is obtained for all the ionic strengths, so that the relaxation of the condition u g = 0 does not seem to affect the accuracy of the LEIP approach.gauche state probabilities versus pH at three ionic strengths are depicted in Fig. 7b: 1M (circles and continuous line), 0.01M (squares and dotted line) and 0.001M (triangles and dashed line).As expected, even at low pH-values some bonds can remain the the gauche state due to the finite value of u g .Despite the complexity of the obtained profiles for this case, LEIP is able to accurately reproduce the MC simulations.

Conclusions
The ionization and conformational properties of polyelectrolytes are determined by a combination of Short Range (SR) and Long Range (LR) interactions between bonds and ionizable sites.In particular, electrostatic LR interactions can only be neglected at high enough ionic strengths, which is an important limitation for many macromolecular systems of interest.The present work explores the possibility of defining local, short ranged, interaction parameters which are corrected to account for the LR interactions in an effective way.The new parameters are systematically calculated by using variational methods, and equations for them are provided.The resulting approach, the Local Effective Interaction Parameters (LEIP) method, was firstly developed to study the binding properties of rigid polyelectrolytes.In this paper these ideas are extended to flexible polyelectrolytes, for which conformational and ionization equilibria (charge regulation) are strongly coupled.With this aim, LEIP is combined with the Site Binding-Rotational Isomeric State (SBRIS) model in order to deal simultaneously with conformational and protonation degrees of freedom for the full range of ionic strengths.
The LEIP method is illustrated by using a model of a linear symmetric polyelectrolyte containing protonating sites situated regularly along the polymer backbone.The bonds linking the ionizable sites can be in three possible states, i.e. trans, gauche+ and gauche-.This model, with only four relevant parameters (protonation pK-value, gauche state energy and SR electrostatic interactions between neighbouring sites through bonds in trans and gauche states) can be regarded as a minimal model of a flexible polyelectrolyte where conformational and binding equilibria are strongly coupled.LEIP method is applied to correct both the protonation pK-values and the gauche state energy.As a result, local pH dependent rotational potentials are obtained.The correction to the gauche energy represents the contribution of the LR interactions in rotating a bond to its gauche state.The degree of protonation and the gauche state probabilities obtained using the LEIP method are compared with those computed using semi-grand canonical Monte Carlo (MC) simulations.In all the studied cases, the agreement between LEIP and MC simulations is excellent.The computational cost, however, is orders of magnitude lower in the LEIP method.This fact allows to use LEIP to directly fit parameters to experimental information.We think that the ideas here presented could be useful in the design of pH-dependent force fields based in experimental ionization and conformational properties.

Figure 2 .
Figure 2.(a) Titration curves corresponding to a rigid linear chain with interacting ionizable sites separated by a distance b=0.2 nm obtained using MC simulations (blue circles), LEIP method correcting only the pK-value (continuous line) and LEIP method correcting the pK-value and the the nearest neighbour interaction energy (red triangles).The chosen parameters are pK=9 and = 1.5.The LR interactions are calculated using the Debye-Hückel potential.The dashed line represents the titration curve in absence of LR interactions.Note that the correction to the pK-value is enough to reproduce almost exactly the MC simulations, and no significant improvement is obtained in correcting .(b) Correction x to the pK-value using the LEIP method.(c) Corrections x (black lines) and x (blue lines) to the pK-value and the nearest neighbour interaction energy , respectively.In all the figures, from top to bottom, the ionic strengths are 1M, 0.5M, 0.1M, 0.05M, 0.01M and 0.001M.

Figure 3 .
Figure 3. Outline of the Local Effective Interaction Parameters (LEIP) method.LEIP accounts for the SR interactions exactly, but LR interactions are replaced by effective SR free energies.In the resulting scheme only SR interaction are present, which considerably simplify the theoretical treatment.The bonds "feel" the presence of the LR interactions in an effective way, which leads to an apparent pH-dependent rotational energy.Other parameters, such as the protonation pK-values, also become pH-dependent due to the presence of LR interactions.

Figure 4 .
Figure 4. Titrations curve for the model polyelectrolyte depicted in Fig.2with pK = 9 , t = 1, u g = 0 (i.e.g → ∞) and σ = 1 (a) ; σ = 10 (b).The chosen ionic strengths are, from top to bottom, 1M, 0.1M, 0.01M and 0.001M.Black continuous lines represent calculations using the LEIP method in which two effective short range parameters (pK and pσ) has been corrected.Red circles depict the MC values.The green dashed line corresponds to the LEIP values for I = 0.001M when only the pK-value is corrected, while pσ is kept constant.

Figure 5 .
Figure 5. gauche state probabilities of bond c versus the pH-value computed by means of MC simulations (red markers) and using the LEIP method (black lines) for ionic strengths: 1M (circles and continuous line), 0.01M (squares and dotted line) and 0.001M (triangles and dashed lines).The chosen parameters are pK = 9 , t = 1, u g = 0 and σ = 1 (a) ; σ = 10 (b).

Figure 6 .
Figure 6.Correction to the gauche state energy x σ versus the pH for pK = 9 , t = 1, u g = 0 and σ = 1 (a) ; σ = 10 (b).From bottom to top the ionic strengths are 1M, 0.1M, 0.01M and 0.001M.x σ represents the average effective energy "felt" by c bonds as a result of the LR interactions.