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

Abstract

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 for analyzing 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, which we will refer to as the Local Effective Interaction Parameters (LEIP) method, is illustrated with a minimal model of a flexible linear polyelectrolyte containing only one type of rotating bond. 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 explicitly taken 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.

1. 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, the importance of the ionization configuration in the conformational properties of intrinsically disordered proteins, whose function-structure relationship still remains a controversial matter, has been recognized [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 Figure 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 is parametrized, the machinery of statistical mechanics can be used in order to quantify the relevant physical properties such as titration curves, site-specific binding probabilities, macroscopic protonation constants, microscopic protonation enthalpies, site–site binding correlation, etc. [5,15,21,22,23,24,25]. For systems with a small number of sites ($N\le 20$), the necessary thermal averages can be performed by direct enumeration, while, for a large number of sites, Monte Carlo (MC) simulations become necessary [26,27,28,29,30,31,32,33,34,35,36,37].

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. The 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 a 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 Figure 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 a 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, semi-grand canonical Monte Carlo simulations are introduced as a tool to test LEIP accuracy when applied to flexible polyelectrolytes. In Section 5, LEIP 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 bond.

2. 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=\left\{s_i\right\}$, 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]

$$\frac{H\left(s\right)}{ln10}=\sum {\mu }_i{s}_i+\sum _{i>j}{\varphi }_{ij}{s}_i{s}_j+\sum _{i>j>k}{\tau }_{ijk}{s}_i{s}_j{s}_k+\dots ,$$

(1)

where ${\mu }_i=\mathrm{pH}-\mathrm{p}{K}_i=-log\left(K_i{a}_{\mathrm{H}}\right)$ is the reduced chemical potential, which depends on the proton activity, $a_{\mathrm{H}}$, and the protonation pK-value of the ionizable site i, $\mathrm{p}{K}_i$; ${\varphi }_{ij}$ represents the interaction energy of the sites i and j; ${\tau }_{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., $\beta =1/k_{\mathrm{B}}T=1$, and divided by a factor $ln10$ in order to be compared in the pH scale. Note that the conformation degrees of freedom are omitted in Equation (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 the 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\left(s\right)$ by means of the semi-grand canonical partition function

$$\mathsf{\Xi }=\sum _s{e}^{-H\left(s\right)}.$$

(2)

The average degree of protonation of a particular site i is related to $\mathsf{\Xi }$ as

$${\theta }_i=〈s_i〉=-\frac{\partial log\mathsf{\Xi }}{\partial {\mu }_i}=\frac{1}{ln10}\frac{\partial \mathsf{\Omega }}{\partial {\mu }_i},$$

(3)

where $\mathsf{\Omega }=-ln\mathsf{\Xi }$ is the thermodynamic potential associated with the semi-grand canonical ensemble. The average number of bound protons is given by

$$\nu =〈\sum _i{s}_i〉=\frac{\partial \mathsf{\Omega }}{\partial ln{a}_{\mathrm{H}}}.$$

(4)

The correlation of the protonation degrees of two sites i and j, a quantity which will be used later, can be expressed as

$$h_{ij}=〈s_i{s}_j〉=-\frac{1}{{\left(ln10\right)}^2}\frac{\partial \mathsf{\Omega }}{\partial {\mu }_i\partial {\mu }_j}.$$

(5)

As can be seen, the quantification of all the relevant physical quantities relies in the accurate determination of the partition function $\mathsf{\Xi }$. If the number of sites is small ($N\le 20)$, $\mathsf{\Xi }$ 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 Figure 1a and assuming only nearest neighbour interactions, the partition function can be expressed as [47]

$$\mathsf{\Xi }=\mathbf{q}{\mathbf{T}}^{\mathit{N}}{\mathbf{p}}^{\mathbf{T}},$$

(6)

where $\mathbf{T}$ is the transfer matrix

$$\mathbf{T}=\left(\begin{array}{cc}1& z\\ 1& zu\end{array}\right).$$

(7)

$z=K{a}_{\mathrm{H}}$ represents the reduced activity and $ϵ=-logu$ is the interaction free energy between neighbouring sites. $\mathbf{q}=(1,0)$ and $\mathbf{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 for including the LR interactions in a very accurate way. In this approach, the full free energy Equation (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 Reference [56]. Although the following arguments can be readily generalized to the general form of the free energy Equation (1), let us consider the simplest case of a rigid linear chain with identical sites, such as the one sketched in Figure 1a. For this system, ${\mu }_1={\mu }_2=\dots =\mu$ and triplet interactions are omitted, i.e., ${\tau }_{ijk}=0$. The reduced free energy H can be split into two contributions $H=H_0\left(x\right)+\Delta H\left(x\right)$ such as

$$\begin{array}{c}\frac{H_0}{ln10}=\left(\mu -x\right){\displaystyle \sum _i}{s}_i+ϵ{\displaystyle \sum _i}{s}_i{s}_{i+1},\\ \frac{\Delta H}{ln10}={\displaystyle \sum _{j>i+1}}{\varphi }_{ij}{s}_i{s}_j+x{\displaystyle \sum _i}{s}_i,\end{array}$$

(8)

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 $ϵ={\varphi }_{i,i+1}$, which can be exactly solved by using the transfer matrix (7). Now, we can use the Gibbs–Bogoliubov variational principle [39,57]

$$\mathsf{\Omega }\le \tilde{\mathsf{\Omega }}={\mathsf{\Omega }}_0\left(x\right)+{〈\Delta H\left(x\right)〉}_0$$

(9)

to determine the optimal value of x, where ${\mathsf{\Omega }}_0\left(x\right)=-ln{\mathsf{\Xi }}_0$ and ${〈\cdots 〉}_0$ represent the free energy and the thermal average corresponding to $H_0$, respectively. Minimizing $\tilde{\mathsf{\Omega }}$ with respect to x, it is found that x fulfills the equation [56]

$$x=\frac{d{\phi }_0/dx}{d{\nu }_0/dx}=\frac{d{\phi }_0}{d{\nu }_0},$$

(10)

where

$${\phi }_0={〈\sum _{j>i+1}{\varphi }_{ij}{s}_i{s}_j〉}_0=\sum _{j>i+1}{\varphi }_{ij}{h}_{ij}^0$$

(11)

is the LR energy averaged over the unperturbed free energy $H_0$, whose correlation function $h_{ij}^0$, 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 $〈\cdots 〉$ by ${〈\cdots 〉}_0$, which can be exactly determined since only SR interactions are involved. Equation (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 the 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 $\mu$ in order to account for the LR interactions but in a local effective way. By the definition of $\mu =\mathrm{pH}-\mathrm{p}\mathit{K}$, x can also be understood as the correction to the site pK-value, so that $\mathrm{p}{K}^{\mathrm{eff}}=\mathrm{p}K-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 Equation (11), although in an approximate way. This approximation, however, results in being extremely accurate, as can be observed in Figure 2a, where the titration curves corresponding to a rigid linear chain with identical sites are depicted. The chosen parameters are $\mathrm{p}K=9$ and $ϵ=1.5$. In this model, the LR interactions between distant sites are described by the Debye–Hückel potential

$${\varphi }_{ij}=\frac{1}{ln10}\frac{{\ell }_{\mathrm{B}}{e}^{-\kappa d_{ij}}}{d_{ij}};j>i+1,$$

(12)

where ${\ell }_{\mathrm{B}}\simeq 0.7$ nm is the Bjerrum length in water at 298 K, $d_{ij}$ is the distance between the sites i and j, and ${\kappa }^{-1}\left(\mathrm{nm}\right)=0.304/\sqrt{I\left(\mathrm{M}\right)}$ is the Debye length at the ionic strength I. For a rigid linear chain, as the one shown in Figure 1, $d_{ij}=\left|j-i\right|b$, 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 Equations (10) and (11) (continuous lines) for all the range of ionic strengths and b = 0.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 many orders of magnitude lower than that required in MC simulations, allowing the fitting of parameters to experimental titration curves. The correction to the pK, x, shown in Figure 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 ($\mathrm{p}K\to \mathrm{p}K-x$), but also the nearest neighbour interaction energy ($ϵ\to ϵ+x_{ϵ}$). Proceeding in the same way, it can be shown that x and $x_{ϵ}$ fulfill the nonlinear system of equations [56]

$$\mathbf{J}\left(\begin{array}{c}x\\ x_{ϵ}\end{array}\right)=\left(\begin{array}{c}{\left(\frac{\partial {\phi }_0}{\partial x}\right)}_{x_{ϵ}}\\ {\left(\frac{\partial {\phi }_0}{\partial x_{ϵ}}\right)}_x\end{array}\right);\mathbf{J}=\left(\begin{array}{cc}{\left(\frac{\partial {\nu }_0}{\partial x}\right)}_{x_{ϵ}}& {\left(\frac{\partial D_0}{\partial x}\right)}_{x_{ϵ}}\\ {\left(\frac{\partial {\nu }_0}{\partial x_{ϵ}}\right)}_x& {\left(\frac{\partial D_0}{\partial x_{ϵ}}\right)}_x\end{array}\right),$$

(13)

where ${\nu }_0\left(x,x_{ϵ}\right)$ and $D_0={〈s_i{s}_{i+1}〉}_0=h_{12}^0\left(x,x_{ϵ}\right)$ represent the average number of protons and the average number of nearest neighbour interactions, respectively, which can be exactly calculated using $H_0$. Solving Equation (13), the correction to the pK and $ϵ$ are obtained as functions of the pH. The physical meaning of x and $x_{ϵ}$ becomes transparent if Equation (13) are rewritten in terms of ${\nu }_0$ and $D_0$ as independent variables. After some elementary algebra, x and $x_{ϵ}$ adopt the much simpler form

$$x={\left(\frac{\partial {\phi }_0}{\partial {\nu }_0}\right)}_{D_0};x_{ϵ}={\left(\frac{\partial {\phi }_0}{\partial D_0}\right)}_{{\nu }_0}.$$

(14)

Equation (14) states that x represents the increase in ${\phi }_0$ for a constant number of interactions $D_0$, while $x_{ϵ}$ can be interpreted as the change in ${\phi }_0$ in creating a nearest neighbour interaction, keeping constant the number of bound protons ${\nu }_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 Figure 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 Figure 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 Figure 2b is no longer present in Figure 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 type (14) can be generalized in a straightforward manner. The same treatment leads to very good results for heterogeneous polyelectrolytes and polyampholytes, by correcting the pK-values of the different kind of sites ($\mathrm{p}{K}_i\to \mathrm{p}{K}_i-x_i$) [56].

3. 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=\left\{c_{\alpha }\right\}$, j = 1, ..., M. The variables $c_{\alpha }$ can adopt several values corresponding to the rotational angles of the bond $\alpha$. 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 ${\mathsf{\Xi }}_{\mathrm{rot}}$, which can be expressed as

$${\mathsf{\Xi }}_{\mathrm{rot}}=\mathbf{q}{\mathbf{U}}_1{\mathbf{U}}_2\dots {\mathbf{U}}_{\mathit{M}-1}{\mathbf{U}}_{\mathit{M}}{\mathbf{p}}^{\mathbf{T}},$$

(15)

where ${\mathbf{U}}_{\alpha }$ is the transfer matrix corresponding to the bond $\alpha$. For a symmetric chain, for which the states gauche+ and gauche− have the same energy and identical bonds, the transfer matrices are of the form

$$\mathbf{U}=\left(\begin{array}{ccc}1& \sigma & \sigma \\ 1& \sigma \psi & \sigma \omega \\ 1& \sigma \omega & \sigma \psi \end{array}\right),$$

(16)

where $\sigma$, $\psi$ and $\omega$ are the Boltzmann factors associated with the conformational energies of the bonds: $-k_{\mathrm{B}}Tln\sigma$ is the free energy of the gauche states while $\psi$ ($\omega$) are related to the interaction energies between two consecutive gauche states of different (same) orientation. $\psi$ and $\omega$ equate one if the rotation of the bonds is independent. $\mathbf{q}=\left(1,0,0\right)$ and $\mathbf{p}=\left(1,1,1\right)$ 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]

$$g=\frac{\partial ln{\mathsf{\Xi }}_{\mathrm{rot}}}{\partial ln\sigma }.$$

(17)

The 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 $\left(s,c\right)$ defines a roto-microstate with reduced free energy

$$F\left(s,c\right)=F_{\mathrm{rot}}\left(c\right)+F_{\mathrm{p}}\left(s,c\right).$$

(18)

$F_{\mathrm{rot}}\left(c\right)$ is the free energy corresponding to the fully deprotonated state of each conformer, while $F_{\mathrm{p}}\left(s,c\right)$ represents the reduced free energy due to the protonation process, which, for a given conformation, can be expressed as

$$\frac{\beta F_{\mathrm{p}}\left(s,c\right)}{ln10}=\sum _i{\mu }_i\left(c\right)s_i+\sum _{j>i}{\varphi }_{ij}\left(c\right)s_i{s}_j,$$

(19)

where triplet interactions have been neglected. Note that the cluster parameters now depend on the conformational state c. The reduced free energy Equations (18) and (19) combines the RIS and the SB model and defines the SBRIS model, which allows for studying conformational and ionization properties on the same foot. In recent publications, 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

$$p\left(s,c\right)=\frac{e^{-\beta F\left(s,c\right)}}{{\mathsf{\Xi }}_{\mathrm{SBRIS}}},$$

(20)

where the SBRIS partition function ${\mathsf{\Xi }}_{\mathrm{SBRIS}}$ is defined as

$${\mathsf{\Xi }}_{\mathrm{SBRIS}}=\sum _{s,c}{e}^{-\beta F\left(s,c\right)}.$$

(21)

The SBRIS partition function can alternatively be expressed in the fashion

$${\mathsf{\Xi }}_{\mathrm{SBRIS}}=\sum _s{\mathsf{\Xi }}_{\mathrm{rot}}\left(s\right),$$

(22)

where ${\mathsf{\Xi }}_{\mathrm{rot}}\left(s\right)$ denotes the rotational partition function for the macromolecule in a ’frozen’ binding configuration $s=\left\{s_1,s_2,\cdots ,s_N\right\}$. ${\mathsf{\Xi }}_{\mathrm{rot}}\left(s\right)$ can then be calculated as a RIS partition function as in Equation (15), 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 (Equation (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

$$\mathbf{U}\to \mathbf{B}=\left(\begin{array}{cc}\mathbf{U}& \mathbf{U}z\\ \mathbf{U}& \mathbf{U}\mathbf{u}z\end{array}\right),$$

(23)

where u is a diagonal matrix containing the Boltzmann factors corresponding to the short-range interactions

$$\mathbf{u}=\left(\begin{array}{ccc}{u}_{\mathrm{t}}& 0& 0\\ 0& u_{\mathrm{g}}& 0\\ 0& 0& u_{\mathrm{g}}\end{array}\right).$$

(24)

In matrix (24), $-k_{\mathrm{B}}Tln{u}_{\mathrm{t}}$ and $-k_{\mathrm{B}}Tln{u}_{\mathrm{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

$$\mathbf{U}\to \mathbf{B}=\left(\begin{array}{cc}\mathbf{U}& \mathbf{0}\\ \mathbf{0}& \mathbf{U}\end{array}\right).$$

(25)

The resulting SBRIS partition function reads

$${\mathsf{\Xi }}_{\mathrm{SBRIS}}=\mathbf{r}{\mathbf{B}}_1{\mathbf{B}}_2\dots {\mathbf{B}}_{\mathit{M}-1}{\mathbf{B}}_{\mathit{M}}{\mathbf{t}}^{\mathbf{T}},$$

(26)

where $\mathbf{r}=\left(\mathbf{q}\mathbf{q}\right)$ and $\mathbf{t}=\left(\mathbf{p}\mathbf{p}\right)$ 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 Equation (26). It can also be 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 on 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 $\mathrm{p}K\to \mathrm{p}K-x$, but also the conformational parameters must be corrected as outlined in Figure 3. In the simplest case, with only one kind of rotating bonds, the substitution $\mathrm{p}\sigma \to \mathrm{p}\sigma +x_{\sigma }$ where $\mathrm{p}\sigma =-log\sigma$ will be necessary. The treatment is almost identical to the one used for rigid molecules. Now, the “unperturbed” free energy is ${\mathsf{\Omega }}_0=-ln{\mathsf{\Xi }}_{\mathrm{SBRIS}}\left(x,x_{\sigma }\right).$ It can be easily shown that the corrections x for the pK and $x_{\sigma }$ for $\mathrm{p}\sigma$ fulfill equivalent equations to (13)

$$\mathbf{J}\left(\begin{array}{c}x\\ x_{\sigma }\end{array}\right)=\left(\begin{array}{c}{\left(\frac{\partial {\phi }_0}{\partial x}\right)}_{x_{\sigma }}\\ {\left(\frac{\partial {\phi }_0}{\partial x_{\sigma }}\right)}_x\end{array}\right);\mathbf{J}=\left(\begin{array}{cc}{\left(\frac{\partial {\nu }_0}{\partial x}\right)}_{x_{\sigma }}& {\left(\frac{\partial g_0}{\partial x}\right)}_{x_{\sigma }}\\ {\left(\frac{\partial {\nu }_0}{\partial x_{\sigma }}\right)}_x& {\left(\frac{\partial g_0}{\partial x_{\sigma }}\right)}_x\end{array}\right),$$

(27)

where ${\nu }_0\left(x,x_{\sigma }\right)$ represents the average number of bound protons (Equation (4)) and $g_0\left(x,x_{\sigma }\right)$ the average number of bonds in the gauche state (Equation (17)), calculated using the unperturbed free energy. If we use ${\nu }_0$ and $g_0$ as independent variables, instead of x and $x_{\sigma }$, Equation (27) can be rewritten in a similar fashion as Equation (14)

$$x={\left(\frac{\partial {\phi }_0}{\partial {\nu }_0}\right)}_{g_0};x_{\sigma }={\left(\frac{\partial {\phi }_0}{\partial g_0}\right)}_{{\nu }_0},$$

(28)

which essentially tell us that x represents the average change in ${\phi }_0$ when a new proton is bound (keeping constant the number of bonds in gauche) while $x_{\sigma }$ is the the average change in ${\phi }_0$ when a bond is brought to its gauche state (keeping constant the number of bound protons). Note that the LEIP method always leads to expressions for the interaction corrections of the same type of Equations (14) and (28).

As in the previous section, LR interactions are described by the Debye–Hückel potential, although the method could in principle be applied to other kind of interactions such as van der Waals interactions. Moreover, by including convenient “hard core” terms in the interaction potentials, the excluded volume effect could in principle be taken into account. The study of this effect, however, is not trivial and it is out of the scope of this work. Unlike rigid molecules, for flexible molecules, the average LR interaction energy ${\phi }_0$ for the unperturbed free energy can only be approximately calculated. In this work, as a first approximation, we assume that

$${\phi }_0={〈\sum _{ij}\varphi \left(d_{ij}\right)s_i{s}_j〉}_0\simeq \sum _{ij}\varphi \left(\sqrt{{〈d_{ij}^2〉}_0}\right){〈s_i{s}_j〉}_0.$$

(29)

This approximation could in principle be improved by using higher order moments of $d_{ij}$. Matricial expressions for ${〈d_{ij}^2〉}_0$ and higher moments where derived by Flory and Jernigan [47,58] 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.

4. Monte Carlo Simulations

In order to estimate the accuracy of the LEIP method when applied to flexible polyelectrolytes, we compare the theoretical values with those resulting from MC simulations. Two main MC techniques have been previously proposed: the Reaction Ensemble approach, for which the pH is a calculated quantity [59,60], and the constant pH method, corresponding to the semi-grand canonical ensemble [33,61,62,63]. Since the control variable in the LEIP method is the pH-value, as indicated by the reduced free energies in Equations (1) and (18), the constant pH method has been chosen here. In previous studies about polyelectrolyte ionization properties, both Reaction Ensemble and constant pH methods have been coupled to Molecular Dynamics schemes in order to deal with explicit ions. In this work, free protons, co- and counter-ions are not explicit in the simulations and the screening effects are taken into account via the Debye length parameter, ${\kappa }^{‒1}$. The MC code generalises the one previously used in the computation of conformational and ionization properties of linear poly(ethylenimine) [46]. The polyelectrolyte is modeled as a linear chain with rigid bond lengths and angles. Bonds can adopt one of the three states of minimum energy (trans, gauche+ or gauche−). Each change of a bond state implies a $120°$ rotation of its dihedral angle and the recalculation of distances among the sites situated before and after the rotating bond. The linear chain is composed of interacting nodes which can correspond to inert or protonating groups. In Figure 4, two snapshots of Monte Carlo simulations at ionic strength 0.001 M and two pH-values (four in Figure 4a and eight in Figure 4b) are presented. As in Figure 1 and Figure 3, the ionizable sites are depicted in blue (dark blue if they are protonated and cyan otherwise). It can be observed that a decrease in the pH-value promotes the elongation of the chain, and the consequent reduction of the electrostatic repulsion, by increasing the number of bonds in the trans conformations.

In the MC simulations, the free energy of the system is divided into SR and LR terms

$$F(s,c)=F_{\mathrm{SR}}+F_{\mathrm{LR}}+\sum _i(\mathrm{pH}-\mathrm{p}{K}_i)s_i.$$

(30)

The SR term is computed using SBRIS free energy (Equation (18)) which involves the energies present in the transfer matrices ($\sigma$, $\psi$, $\omega$, $u_{\mathrm{t}}$ and $u_{\mathrm{g}}$), while the LR contribution is calculated using the Debye–Hückel potential (Equation (12)). If $F_{LR}$ is set to zero, the obtained results coincide, within the numerical error, with those obtained using the transfer matrix method. This was one of the tests used to check the reliability of the Monte Carlo code. A Metropolis algorithm [15,27] is used to generate roto-microstates at constant pH in a chain with 50 ionizable sites (i.e., 148 nodes or 147 bonds). In each new MC configuration, the polyelectrolyte can change either (i) the conformational state of a rotating bond or (ii) the ionization state of a binding site, with trial probabilities of 0.999 and 0.001, respectively. These trial probabilities allow us to obtain a good equilibration of the conformational structure for each ionization state and the system does not become trapped in local minima. The probability to accept a new configuration is obtained by computing the free energy difference ($\Delta F(s,c)$) between trial and actual conformations. When the state of the bond $\alpha$ is changed, the following free energy differences must be calculated: (i) the conformational energy of bond $\alpha$ and its interaction with bonds $\alpha \pm 1$ (corresponding to the parameters $\sigma$, $\psi$ and $\omega$); (ii) the electrostatic SR interaction between the two sites bound to $\alpha$ when they are charged (corresponding to $u_{\mathrm{t}}$ and $u_{\mathrm{g}}$, which depend on the new conformation of $\alpha$ ); and (iii) the change in the LR Debye–Hückel interaction among sites before and after $\alpha$, which involves the recalculation of the distances between the charged sites. On the other hand, a change in the ionization state of a site $s_i$ implies to recalculate: (i) the reduced chemical potential of the site i by an amount $\Delta {\mu }_i=(\mathrm{pH}-\mathrm{p}{K}_i)\Delta n$, where $\Delta n=\pm 1$ is the variation in the number of protons; (ii) the SR repulsive interaction between $s_i$ and $s_{i\pm 1}$; and (iii) the LR Debye–Hückel interactions between the trial protonating site and the rest of ionized sites. Once $\Delta F(s,c)$ is computed, the new configuration is always accepted if $\Delta F(s,c)<0$ and accepted with a probability $exp(-\beta \Delta F(s,c\left)\right)$ if $\Delta F(s,c)>0$. The values presented are the average over eight different MC simulations. Each MC simulation has been equilibrated in the first $5\times {10}^7$ configurations and the thermal averages have been computed in the following $4.5\times {10}^8$ realizations. The simulations were performed using a parallel code developed in C++ on a 126 CPU cluster. For each pH and ionic strength (one point of the curves), typical jobs were run using 8 CPUs during 1 to 2 h.

5. Results and Discussion

As a model system, we use the linear polyelectrolyte outlined in Figure 1b, with protonating sites situated every three chain positions. Only 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 $\mathrm{p}\sigma =-log\sigma$. The c bonds are assumed to rotate independently when the macromolecule is uncharged ($\omega =\psi =1$ in Equation (16)). The protonating sites are considered to be identical with the same protonation pK-value. The interactions between protonated sites are characterized by the energies ${ϵ}_{\mathrm{t}}=-log{u}_{\mathrm{t}}$ (when the bond c is in trans state) and ${ϵ}_{\mathrm{g}}=-log{u}_{\mathrm{g}}$ (when the bond c is in gauche state). In the computations, the used values of the bond length and the bond angle are 0.2 nm and 120°, respectively. This model can be regarded as a minimal model of a flexible polyelectrolyte, with only four energetic parameters involved ($\sigma$, ${ϵ}_{\mathrm{t}}$, ${ϵ}_{\mathrm{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.

Let us firstly consider the case for which c bonds can freely rotate when the adjacent sites are deprotonated (i.e., $\sigma =1$). When both sites are charged, however, the very strong SR repulsion hinders the gauche conformation, so that we take $u_{\mathrm{g}}=0$ (i.e., ${ϵ}_{\mathrm{g}}\to \infty$). The resulting titration curves are shown in Figure 5a for ionic strengths ranging from 1 M to 0.001 M. The chosen parameters are $\mathrm{p}K=9$ and ${ϵ}_{\mathrm{t}}=1$. The black continuous lines represent the average protonation degree $\theta$ calculated using the LEIP method correcting both the pK-value ($\mathrm{p}K\to \mathrm{p}K-x$) and the conformational energy of c bonds ($\mathrm{p}\sigma \to \mathrm{p}\sigma +x_{\sigma }$), 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.001 M if only the pK-value is corrected, while $\sigma$ remains constant. Although relative good prediction of the MC simulations is obtained, the quality of the titration curve clearly improves if $\sigma$ 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 behaves 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 when the gauche states of the c bonds are favored, for instance, because of the existence of hydrogen bond, which means that $\sigma >1$. When the adjacent sites are both protonated, on the contrary, the electrostatic repulsion is so strong that the gauche states are forbidden ($u_{\mathrm{g}}=0$). In this case, the conformational propensity changes when the ionization state of the sites change. Figure 5b compares the titration curves obtained using LEIP correcting pK and $\mathrm{p}\sigma$ and MC simulations for $\mathrm{p}K=9$, ${ϵ}_{\mathrm{t}}=1$ and $\sigma =10.$ As can be observed, LEIP method and MC simulations yield to almost identical titration curves. In this case, however, the correction of $\mathrm{p}\sigma$ becomes compulsory. If only $\mathrm{p}K$ is corrected, the titration curve obtained by LEIP at 0.001 M (green dashed line) exhibits a phase transition-like behaviour at $\mathrm{pH}\approx 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 Figure 6a,b. Markers correspond to MC simulations while black lines represent the theoretical values at ionic strengths 1 M, 0.01 M and 0.001 M, for $\mathrm{p}K=9$, ${ϵ}_{\mathrm{t}}=1$ and $u_{\mathrm{g}}=0$. Figure 6a corresponds to the case with $\sigma =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 the 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 = 1 M, the LR interactions can be neglected and total correspondence between simulated and calculated values is found. At higher ionic strengths (0.01 M and 0.001 M), 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.

Figure 6b corresponds to the case with $\sigma =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\sigma /(2\sigma +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 ($\mathrm{p}\sigma \to \mathrm{p}\sigma +x_{\sigma }$). 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 Figure 5b (green dashed line).

Figure 7 shows the LEIP method correction $x_{\sigma }$ to the bond conformational energy ($\mathrm{p}\sigma \to \mathrm{p}\sigma +x_{\sigma }$) for $\sigma =1$ (Figure 7a) and $\sigma =10$ (Figure 7b). In both cases, it is observed that $x_{\sigma }$ 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_{\sigma }$ 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 $\sigma =10$, a wavy behaviour is observed for ionic strengths 0.1 M and 0.01 M and $x_{\sigma }$ exhibits a smooth maximum at p$H\simeq 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 included 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_{\mathrm{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_{\mathrm{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 Figure 8. Figure 8a shows the computed titration curves with ${ϵ}_{\mathrm{g}}=2$ at ionic strengths ranging, from top to bottom, from 1 M to 0.001 M. Excellent agreement between the theoretical predictions and simulations is obtained for all the ionic strengths, so that the relaxation of the condition $u_{\mathrm{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 Figure 8b: 1 M (circles and continuous line), 0.01 M (squares and dotted line) and 0.001 M (triangles and dashed line). As expected, even at low pH-values, some bonds can remain in the gauche state due to the finite value of $u_{\mathrm{g}}$. Despite the complexity of the obtained profiles for this case, LEIP is able to accurately reproduce the MC simulations.

6. 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 charged sites interact by means of the Debye–Hückel potential, which accounts for the electrostatic screening in an average way, while excluded volume effects are neglected. 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. The 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 of 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 using LEIP to directly fit parameters to experimental information. The LEIP method could also represent a complementary tool to the study of other aspects of the polyelectrolyte physical chemistry, such as the dependence of the molecular size on the pH, the influence of excluded volume interactions, the presence of attractive hydrophobic interactions or the competitive binding of metal ions. The clarification of these points, which have not been the subject of the present study, would be desirable in order to extend the applicability of the LEIP method. We think that the ideas presented here could be useful in the design of pH-dependent force fields based on experimental ionization and conformational properties.