Simulation of weak polyelectrolytes: A comparison between the constant pH and the reaction ensemble method

TL;DR

This study compares the reaction ensemble and constant pH methods for simulating protonation in weak polyelectrolytes, highlighting their similarities and differences, especially in systems with strong acids or confined geometries.

Contribution

It demonstrates the conditions under which both methods agree and identifies the advantages of the reaction ensemble approach for complex systems with diverse pKa values.

Findings

Good agreement between methods when reaction constant is swept

Pronounced differences for strong acids/bases due to explicit protons

Reaction ensemble better for confined geometries and diverse pKa values

Abstract

The reaction ensemble and the constant pH method are well-known chemical equilibrium approaches to simulate protonation and deprotonation reactions in classical molecular dynamics and Monte Carlo simulations. In this article, we show similarity between both methods {under certain conditions}. We perform molecular dynamics simulations of a weak polyelectrolyte in order to compare the titration curves obtained by both approaches. Our findings reveal a good agreement between the methods when the reaction ensemble is used to sweep the reaction constant. Pronounced differences between the reaction ensemble and the constant pH method can be observed for stronger acids and bases in terms of adaptive pH values. These deviations are due to the presence of explicit protons in the reaction ensemble method which induce a screening of electrostatic interactions between the charged titrable groups of the polyelectrolyte. The outcomes of our simulation hint to a better applicability of the reaction ensemble method for systems in confined geometries and titrable groups in polyelectrolytes with different pK$_\text{a}$ values.

Equations

$$\text{nˉ}=\frac{N_{\text{HA}}}{N_0}$$

$$K_a=\frac{c({\text{A}}^{-})c(\text{H+})}{c(\text{HA})}$$

$$p(r)t(l|r)=p(l)t(r|l),$$

$$t(r|l)=g(r|l)\text{acc}(r|l)$$

$$\frac{\text{acc}(l|r)}{\text{acc}(r|l)}=\frac{p(l)g(r|l)}{p(r)g(l|r)}.$$

$$\text{acc}(l|r)=\min \left(1,\frac{p(l)g(r|l)}{p(r)g(l|r)}\right),$$

$$\text{acc}(l|r)=\min \left(1,\frac{p(l)}{p(r)}\right)$$

$$\sum _{i=1}^z{\nu }_i{s}_i=0$$

$$\text{acc}^{\text{RE}}(l|r)=\min\bigg{\{}1,(\beta P^{0}V)^{\overline{\nu}\xi}K^{\xi}\prod_{i=1}^{z}\left[\frac{N_{i}^{0}!}{(N_{i}^{0}+\xi\nu_{i})!}\right]\exp(-\beta\Delta E_{\text{pot}})\bigg{\}},$$

$$\text{acc}(l|r)=\min\left(1,\exp\left[-\beta\left(\Delta E_{\text{pot}}+\bigg{(}\pm(\ln(10)/\beta)\left(\text{pH}_{\text{in}}-\text{pK}_{a}\right)\bigg{)}\right)\right]\right)$$

$$Z_{\text{pH}}=\sum _{\text{nˉ}}\left(\genfrac{}{}{0ex}{}{N_0}{(1-\text{nˉ})N_0}\right)x^{N_0(1-\text{nˉ})}\sum _{i(\text{nˉ})}\exp (-\beta E_{\text{pot,i}}),$$

$$p(\text{nˉ},E_{\text{pot, i}})=\left(\genfrac{}{}{0ex}{}{N_0}{(1-\text{nˉ})N_0}\right)x^{N_0(1-\text{nˉ})}\exp (-\beta E_{\text{pot,i}})$$

$$\widetilde{\text{acc}}(\text{diss}|\text{ass})=\min \left(1,\frac{p(\text{nˉ}-\Delta \text{nˉ},E_{\text{pot,diss}})}{p(\text{nˉ},E_{\text{pot,ass}})}\right)$$

$$\widetilde{\text{acc}}(\text{diss}|\text{ass})=\min \left(1,\frac{\left(\genfrac{}{}{0ex}{}{N_0}{(1-\text{nˉ}+\Delta \text{nˉ})N_0}\right)}{\left(\genfrac{}{}{0ex}{}{N_0}{(1-\text{nˉ})N_0}\right)}{x}^{N_0\Delta \text{nˉ}}\exp (-\beta \Delta E_{\text{pot}})\right)$$

$$\widetilde{\text{acc}}(\text{diss}|\text{ass})=\min \left(1,\frac{\left(\genfrac{}{}{0ex}{}{N_0}{(1-\text{nˉ}+1/N_0)N_0}\right)}{\left(\genfrac{}{}{0ex}{}{N_0}{(1-\text{nˉ})N_0}\right)}x\exp (-\beta \Delta E_{\text{pot}})\right)$$

$$\frac{\left(\genfrac{}{}{0ex}{}{N_0}{(1-\text{nˉ}+1/N_0)N_0}\right)}{\left(\genfrac{}{}{0ex}{}{N_0}{(1-\text{nˉ})N_0}\right)}=\frac{N_0\text{nˉ}}{N_0(1-\text{nˉ})+1}\stackrel{N_0\to \infty }{\sim }\frac{N_0\text{nˉ}}{N_0(1-\text{nˉ})}=\frac{N_{\text{HA}}}{N_{{\text{A}}^{-}}}$$

$$\widetilde{\text{acc}}(\text{diss}|\text{ass})=\min \left(1,\frac{N_{\text{HA}}}{N_{{\text{A}}^{-}}}10^{{\text{pH}}_{\text{in}}-{\text{pK}}_{\text{a}}}\exp (-\beta \Delta E_{\text{pot}})\right)$$

$${\text{acc}}^{\text{RE}}(\text{diss}|\text{ass})=\min \left(1,K_a\frac{N_{\text{HA}}}{{\text{N}}_{{\text{A}}^{-}}{c}^{\ast }({\text{H}}^{+})}\exp (-\beta \Delta E_{\text{pot}})\right)$$

$$K_a\frac{N_{\text{HA}}}{N_{{\text{A}}^{-}}{c}^{\ast }({\text{H}}^{+})}=K_a\frac{N_{\text{HA}}}{N_{{\text{A}}^{-}}c({\text{H}}_{\text{in}}^{+})}$$

$${\text{pK}}_{\text{a}}-{\text{pH}}^{\ast }={\text{pK}}_{\text{a}}-{\text{pH}}_{\text{in}}$$

$$U_{\text{FENE}}(r)=-\frac{1}{2}k{r}_{\text{max}}^2\log \left(1-{\left(\frac{r-r_0}{r_{\text{max}}}\right)}^2\right)$$

$$m_i{\ddot{\vec{r}}}_i=-\zeta {\dot{\vec{r}}}_i+{\vec{R}}_i+{\vec{F}}_i$$

Full text

1 Introduction

Weak polyelectrolytes, like polyacrylic acid or most proteins [1] have titrable groups which can be either in a protonated or a deprotonated state depending on the pH value of the solution [2]. The influence of the pH value gives rise to phenomena like protonation-configuration [3, 4], or charge regulation effects [5, 6] as they are known for weak polyelectrolytes like proteins. A protonation/deprotonation reaction of a titrable group in a weak polyelectrolyte can be written as

\ce

HA ¡=¿ A^- + H^+

where $\rm{HA}$ denotes the protonated form of the titrable group, $\rm{A^{-}}$ the deprotonated form and $\rm{H^{+}}$ the dissociated proton. The presence of water as a proton acceptor or proton donor is omitted for the sake of simplicity. Depending on the pH value defined by $\rm{pH=-\log_{10}(c(H^{+})/(\text{mol/L)})}$, the degree of association can be calculated via

[TABLE]

with the number of associated titrable groups $N_{\textrm{HA}}$ divided by the total number of titrable groups N_{0}=N_{\textrm{HA}}+N_{\textrm{A{}^{-}}}, where N_{\textrm{A{}^{-}}} denotes the number of deprotonated units [7, 8]. The equilibrium concentration of each species is steered by an apparent reaction constant in accordance to the law of mass action

[TABLE]

where $c(\cdot)$ denotes the concentration of the individual species [9]. More conveniently, the logarithmic reaction constant is defined by $\rm{pK_{a}=-\log_{10}(K_{a}/(\text{mol/L)})}$, where typical values for $\rm{pK_{a}}$ vary between $\rm{pK_{a}=2-10}$ for common weak polyelectrolytes [10].

The simulation of weak polyelectrolytes in classical molecular dynamics (MD) or Monte Carlo (MC) simulations is a challenging task. Over the last decades, the constant pH and the reaction ensemble method were the most frequently used algorithms to perform these simulations. The constant pH method [11] implements a constant and global pH value as an input parameter which balances the probability of protonation and deprotonation reactions. The method was originally developed to simulate linear polyelectrolytes in presence and absence of excess salt and counterions [11, 12, 13, 14] and reveals a good agreement with analytical results [15]. More effort was additionally spent on the study of different solvent conditions and their influence on polyelectrolyte conformations [16, 17] and the properties of polyampholytes [18, 5]. Moreover, the constant pH method can also be employed in slightly modified versions [19] in order to analyze shifts in the apparent reaction constant for proteins.

A different approach was introduced by the reaction ensemble (RE) method which also provides the possibility to model arbitrary chemical reactions in classical simulations [20, 21, 22]. The reaction ensemble method was used to study chemical and phase equilibria in different systems and to investigate the influence of high pressure and high temperature on chemical reactions [20]. Furthermore, it can be applied for acidic molecules in confined geometries, for the study of interface effects [20] and for the simulation of acid-base reactions in weak polyelectrolyte systems [23, 24, 25, 26]. In contrast to the constant pH method, the reaction ensemble method considers an adaptive pH value which significantly differs from the global interpretation of the predefined implicit pH value in the constant pH method.

In this article, we elucidate the main properties of the constant pH and the reaction ensemble method. Moreover, we also demonstrate similarity between the reaction ensemble method and the constant pH method under certain conditions. Furthermore, a new interpretation of the reaction ensemble method is proposed in order to reproduce a real titration experiment. Our MD simulations verify that the explicit presence of protons in the reaction ensemble method induces electrostatic screening effects between the charged titrable groups of a polyelectrolyte in contrast to the constant pH simulations. Hence, differences in the results obtained by both methods can be assigned to the different implementations of predefined and adaptive pH values.

The article is organized as follows. In the next section, we discuss the main properties of the reaction ensemble and the constant pH method. Furthermore, we present a novel interpretation of the reaction ensemble method which is useful for the study of titration curves. In section 3, we show similarity between the constant pH and the reaction ensemble method under certain conditions. After the presentation of the numerical details in section 4, we compare the titration curves obtained by both methods for weak polyelectrolytes in coarse-grained MD simulations. We briefly conclude and summarize in the last section.

2 Properties of the constant pH and the reaction ensemble method

The reaction ensemble and the constant pH method rely on Monte Carlo techniques and can be implemented in terms of the Metropolis-Hastings algorithm [27, 28]. The Metropolis-Hastings algorithm is an importance sampling technique which fosters the population of states with a certain probability, for example the canonical probability $p_{i}=\exp(-\beta E_{i})/Z$ with the energy of a state $i$, the inverse thermal energy $\beta=1/k_{B}T$ with the Boltzmann constant $k_{B}$, the temperature $T$ and the canonical partition sum $Z$. A detailed balance condition between two states $r$ and $l$ is defined by

[TABLE]

with the probability $p$ and the associated transition probability $t$ between the two states. The transition probability itself can be calculated by the definition of a proposal probability $g$ and an acceptance probability acc according to

[TABLE]

which can be inserted into Eqn. (3) to yield

[TABLE]

A standard choice [28] for the acceptance probability $\text{acc}(l|r)$ reads

[TABLE]

which fulfills the requirements of the detailed balance condition (Eqn. (3)). In the original Metropolis algorithm [29], the proposal probability is symmetric ($g(l|r)=g(r|l)$) and therefore the acceptance probability simplifies to

[TABLE]

for arbitrary choices of $g$. Specific expressions of $g$ are mostly intended to achieve an efficient and effective sampling of the phase space.

2.1 The reaction ensemble method

In presence of chemical equilibrium, the reaction ensemble method proposes changes in the particle numbers of the reacting species by forward (deprotonation) and backward (protonation) reactions [21, 20, 22]. As outlined by Turner et al. [20], the definition of the reaction ensemble with fluctuating particle numbers can be derived from the grand canonical ensemble via the separation of the kinetic and the configurational canonical partition sum.

In general, a chemical reaction can be written as

[TABLE]

for $z$ chemical species of type $s_{i}$ with stoichiometric coefficients $\nu_{i}$ [9]. The acceptance probability in the reaction ensemble for an arbitrary forward reaction from state $r$ to $l$ is defined as

[TABLE]

where $N_{i}$ is the number of particles after a reaction, $N_{i}^{0}$ the number of particles prior to a reaction and $\xi$ the “extent” of the reaction which is selected randomly with $\xi\pm 1$ [20]. A deprotonation (forward) reaction is defined by $\xi=+1$ and a protonation (backward) reaction by $\xi=-1$. Additional parameters are the dimensionless reaction constant $K$ for each titrable group, which is proportional to the apparent reaction constant in Eqn. (2), the standard pressure $P^{0}$, the potential energy difference with $\Delta E_{\text{pot}}=E_{\text{pot,r}}-E_{\text{pot,l}}$, the volume of the system $V$ and the total change in the number of particles $\overline{\nu}=\sum_{i}\nu_{i}$. The corresponding protonation and deprotonation reactions are usually performed after a predefined number of MD simulation steps with constant particle numbers [20].

2.2 Operation modes of the reaction ensemble method

Two different options on how to apply the reaction ensemble method can be defined. We coin the first operation mode “sweeping the reaction constant” and our proposed second operation mode “real titration”.

The “sweeping” operation mode, which was for example applied in Refs. [23, 24, 25, 26] considers multiple independent simulation runs in which the individual titrable units are characterized by various dimensionless input reaction constants $K=K_{a}/(\beta P_{0})^{\overline{\nu}}$ [26]. This definition has a direct chemical interpretation: in each independent simulation, titrable units with a fixed and arbitrary reaction constant $K$ are inserted into initially neutral water solution. After equilibration of the reaction ensemble, a certain pH value in the simulation box is adjusted, which typically deviates from the neutral pH value. The corresponding pH value can be regarded as the equilibrium pH value for the given choice of the reaction constant $K$ and the given surrounding. In the following, we denote this pH value as the “eigen pH value”. By choosing various values of $K$, the resulting degrees of association as well as the corresponding (eigen) pH values can be obtained and therefore also the titration curves as a function of the pKa-pH value.

As a second operation mode, we propose to directly imitate a “titration” experiment [9]. In a real titration experiment, a certain substance is titrated by injecting a strong acid or a strong base into the system. It is important to note, that the pH value changes in contrast to the reaction constant, which is a fixed substance property. Therefore, in the real titration mode that we propose, we fix the intrinsic reaction constant and add certain amounts of a strong acid (*e. g. * HCl with species \ceH+ and \ceCl-) or a strong base (*e. g. * NaOH with species \ceNa+ and \ceOH-) to the system. Due to their chemical properties [9], strong acids and bases reveal a very high dissociation constant and therefore remain dissociated even at extremely high or low pH values, respectively. Moreover, also the autoprotolysis reaction of water \ce2 H2O ¡=¿ H_3O+ +OH- with an apparent reaction constant of $K_{w}=10^{-14}\text{mol}^{2}/\text{L}^{2}$ is explicitly taken into account. Therefore, the protons and hydroxide ions can react to neutral water molecules which are neglected in our implicit solvent approach. It has to be noticed that hydroxide ions and protons differ from ordinary counterions (\ceNa+ and \ceCl-) only due to their properties in our algorithm, which enables them to participate in chemical reactions (autoprotolysis and deprotonation/protonation reactions of the polyelectrolyte), whereas all other properties are identical.

In order to approach higher or lower pH values than the eigen pH value of the weak acid, a strong base or a strong acid, respectively, are injected into the system. Based on the autoprotolysis reaction of water, the concentration of the deprotonated or protonated titrable units of the weak polyacid directly adapts to the concentration of the excess H+ or OH- species in the solution. In fact, this operational mode resembles an experimental titration procedure due to the explicit presence of protons or hydroxide ions and cannot be imitated by the constant pH method without further effort.

2.3 The constant pH method

In the constant pH method, the protonation or deprotonation probability for a titrable group is determined after two steps. First, a random titrable group is chosen. If it is of type A- or HA, the group or particle, respectively and the corresponding properties are exchanged. Thus, dissociated protons are randomly placed in the simulation box and charge neutrality is fulfilled. The trial move is accepted with a probability

[TABLE]

where $\Delta E_{\text{pot}}$ is the potential energy change due to the exchange of chemical species, pH${}_{\text{in}}$ is an input parameter which determines the implicit pH value of the solution and $\rm{pK_{a}}$ is the negative common logarithm of the apparent reaction constant, which is also a simulation parameter with a predefined value. The expression $\pm(1/\beta)\left(\text{pH}_{\text{in}}-\text{pK}_{a}\right)$ can be interpreted as a change of the chemical potential [11]. A negative prefactor defines a deprotonation reaction (diss) and a positive prefactor, vice versa assigns a protonation reaction (ass). Moreover, it has to be noted that the proposal probability for a protonation or deprotonation reaction in the constant pH method is asymmetric [11, 30, 31, 15, 16, 17, 14, 18]. This can be shown by comparing the proposal probability for a deprotonation reaction which reads $g(\text{diss}|\text{ass})=N_{\text{HA}}/N_{0}$, with the proposal probability for a protonation reaction defined as $g(\text{ass}|\text{diss})=N_{\text{A}^{-}}/N_{0}$, which implies $g(\text{ass}|\text{diss})\neq g(\text{diss}|\text{ass})$.

3 The constant pH method with a symmetric proposal probability

As we have discussed above, the constant pH method relies on asymmetric proposal probabilities. In order to show similarity between the constant pH and the reaction ensemble method under certain conditions, we develop an expression for the constant pH method with a symmetric proposal probability which is then compared to the reaction ensemble method. Noteworthy, the acceptance probability can be evaluated by the original Metropolis algorithm [29] with symmetric proposal probabilities per definition. In order to follow this approach, we use an expression for the partition sum of the constant pH ensemble which was proposed in Ref. [11] and reads

[TABLE]

as a sum over all degrees of association and over all corresponding configurational microstates $i$ of the system. The individual probability for a microstate with a certain degree of association reads

[TABLE]

with $x=10^{\text{pH}_{\text{in}}-\text{pK}_{a}}$ and predefined and fixed values for $\text{pH}_{\text{in}}$ and $\text{pK}_{a}$. A deprotonation step for a single titrable group can be expressed by a change of the degree of association \Delta\text{\bar{n}} in order to describe the transition from (\text{\bar{n}},E_{\text{pot,ass}}) to (\text{\bar{n}}-\Delta\text{\bar{n}},E_{\text{pot,diss}}). Thus, the Metropolis acceptance probability [29] for this Monte Carlo move reads

[TABLE]

which yields

[TABLE]

after inserting Eqn. (12) into Eqn. (13). The equation above can be reformulated for a single deprotonation step in order to read

[TABLE]

with \Delta\text{\bar{n}}=1/N_{0}. By using the relation

[TABLE]

in the thermodynamic limit for an infinite number of titrable groups $N_{0}$, we finally obtain a simple expression for the acceptance probability in the constant pH method with a symmetric proposal probability according to

[TABLE]

which can be also derived for a protonation reaction. One has to notice that the so derived acceptance probability $\tilde{\text{acc}}$ in Eqn. (17) differs from the standard acceptance probability in Eqn. (10) with regard to the prefactor $N_{\text{HA}}/N_{\text{A}^{-}}$ which accounts for the usage of a symmetric proposal probability.

3.1 Similarity between the reaction ensemble and the constant pH method under certain conditions

As it was discussed in the introduction, one can either use the reaction ensemble or the constant pH method for the simulation of weak polyelectrolytes. In this section, we demonstrate similarity between the constant pH method and the reaction ensemble method in the sweeping mode, as it was introduced in section 2.2 under certain conditions. In terms of dissociation reactions, the reaction ensemble yields the acceptance probability

[TABLE]

which is a simplified version of Eqn. (9) with the apparent reaction constant $K_{a}=K\beta P^{0}$ and the currently present proton concentration $c^{*}(\text{H}^{+})=N_{\text{H}^{+}}/V$ in the simulation box. A comparison between the acceptance probability $\tilde{\text{acc}}$ in Eqn. (17) for the constant pH method with a symmetric proposal probability and Eqn. (18) yields that both acceptance probabilities are equal if the following relation hold

[TABLE]

with $10^{\text{pH}_{\text{in}}-\text{pK}_{\text{a}}}=K_{\text{a}}/c(\text{H}_{\text{in}}^{+})$ that can be also expressed by

[TABLE]

with the implicit and predefined $\text{pH}_{\text{in}}$ value as used in the constant pH method denoted by $\text{pH}_{\text{in}}=-\log_{10}(c(\text{H}_{\text{in}}^{+})/(\text{mol/L}))$ including the virtual proton concentration $c(\text{H}_{\text{in}}^{+})$ and $\text{pH}^{*}=-\log_{10}(c^{*}(H^{+})/(\text{mol/L}))$ the current pH in the simulation box. Eqn. (20) is valid for the sweeping operational mode in the reaction ensemble method and for $\text{pH}^{*}=\text{pH}_{\text{in}}$. If these requirements are fulfilled, the reaction ensemble method and the reformulated constant pH method reveal equal acceptance probabilities. Since the current ${\text{pH}^{*}}$ in the reaction ensemble simulation fluctuates around the $\text{pH}:=-\log_{10}(\langle c(\text{H}^{+})\rangle)$ in the reaction ensemble the average particle number in the constant pH method and the reaction ensemble in the sweeping mode are the same if $\text{pH}=\text{pH}_{\text{in}}$. However the variance of the particle number is typically different in both methods.

At this point, it is important to note that the pH value in the reaction ensemble method (pH) is measured via the actual proton concentration whereas the pH value in the constant pH method (pH${}_{\text{in}}$) represents a constant input parameter. Hence, a change of the box volume $V$ in the reaction ensemble method induces a variation of the measured pH value in contrast to the constant pH method where the pH is a fixed number. It thus follows, that the reaction ensemble method allows the study of concentration dependent effects in terms of the law of dilution, which enforces a more pronounced deprotonation behavior for lower concentrations of titrable groups [9]. In fact, the constant pH method can be interpreted as a coupling scheme to an implicit proton bath of infinite dimensions which fixes the pH value of the solution. In contrast to the reaction ensemble method in the real titration mode, the absence of all free protons in the constant pH method reduces a screening of electrostatic interactions between the charged titrable groups. Thus, the resulting deprotonation behavior differs between the methods which can be recognized by differences in the titration curves as it will be discussed in the next sections.

4 Simulation details

We study the properties of weak polyelectrolytes in terms of a coarse-grained bead-spring model with $N_{0}=50$ beads. All titrable groups (beads) repel each other by a truncated and shifted Lennard-Jones potential [32] with amplitude $\epsilon=1k_{B}T$ and range $1\sigma$ yielding a cutoff radius $r_{c}=2^{1/6}\sigma$. Electrostatic interactions were calculated by the P3M method [33] with a Bjerrum length $\lambda_{B}=e^{2}/4\pi\varepsilon_{0}\varepsilon_{\rm r}k_{\rm B}T=2\sigma$ including the dielectric constant $\epsilon_{\rm r}$ and the elementary charge $e$. In comparison to an aqueous solution at room temperature yielding $\lambda_{B}=0.71$ nm, we thus identify $\sigma=0.355$ nm. The cubic simulation box with periodic boundary conditions in all three dimensions has a box length of $b=56.3124\sigma$ with a monomer or titrable group concentration of $c_{0}=0.00028\sigma^{-3}$ and a polymer concentration of $c_{p}=5.6\cdot 10^{-6}\sigma^{-3}$. With the Avogadro constant $N_{A}=6.022\cdot 10^{23}$ mol*-1* and $\sigma=0.355$ nm, these values correspond to concentrations of $c_{0}\approx 0.01$ mol/L and $c_{p}\approx 2\cdot 10^{-3}$ mol/L. Bonds between adjacent beads of the polyelectrolyte are modeled by a FENE potential [34] according to

[TABLE]

with the spring constant $k=10\epsilon/\sigma^{2}$, the maximum elongation $r_{\text{max}}=1.5\sigma$ and an equilibrium length $r_{0}=2^{1/6}\sigma\approx 1.12\sigma$. We perform Langevin Dynamics simulations according to

[TABLE]

with the mass $m_{i}=1m$ for each particle, the conservative force $\vec{F}_{i}$, the friction force $-\gamma\dot{\vec{r}}_{i}$ and the random force ${\vec{R}}_{i}$. The random force acts on each particle independently and obeys the fluctuation-dissipation theorem $\langle R_{ik}\rangle=0$ and $\langle R_{ik}(t)R_{jl}(t^{\prime})\rangle=2\gamma k_{B}T\delta_{ij}\delta_{kl}\delta(t-t^{\prime})$ which ensures the presence of Gaussian white noise for particles $i$ and $j$ in the spatial directions $k$ and $l$. The friction coefficient has a value of $\gamma=1\sigma^{-1}(m\epsilon)^{1/2}$. The temperature is $T=1\epsilon/k_{B}$ and the Langevin equation is integrated by a Velocity Verlet algorithm with a time step of $\delta t=0.01\sigma(m/\epsilon)^{1/2}$.

The apparent reaction constant $K_{a}=K\beta P^{0}$ in the reaction ensemble depends on the standard pressure $P^{0}=1\text{ bar}=0.00108\epsilon/\sigma^{3}$ and on $\beta=1\epsilon^{-1}$. For the constant pH method, we varied the values for $\rm{pK_{a}-pH}$ between -4 and 2. In the sweeping mode of the reaction ensemble, we choose dimensionless reaction constants $K$ between the values -8 and 2 in logarithmically equidistant intervals. In contrast, we consider values of $\text{pK}_{a}=0.49$ and $\text{pK}_{a}=3$ in the real titration mode and vary the number of negatively and positively charged excess protons or hydroxide ions at specific pH values. Moreover, we also take the autoprotolysis of water into account by adding the apparent autoprotolysis reaction constant $K_{\text{w}}=10^{-14}\rm{mol^{2}/L^{2}}=10^{-14}\cdot(0.02694/\sigma^{3})^{2}$ [9]. All simulations are performed with the software package ESPResSo [35, 36].

5 Numerical Results

In Fig. 1, we compare the titration curves between the constant pH and the reaction ensemble method in the sweeping mode for a flexible weak polyelectrolyte with Bjerrum length $\lambda_{B}=2\sigma$. As it was discussed in Sec. 3.1, a nearly perfect agreement between the results for the reaction ensemble and the constant pH method can be observed. Moreover, strong deviations to an ideal titration curve with $\Delta E_{\text{pot}}=0$ can be seen.

The results of the reaction ensemble method in the real titration mode and the constant pH method for different pKa values are presented in Fig. 2. Depending on the intrinsic $\text{pK}_{a}$ value of the titrable units and the pH value of the solution, the titration curves of the reaction ensemble deviate significantly from those obtained by the constant pH method. As an example, for a moderately strong acid with $\rm{pK_{a}=0.49}$, one can observe pronounced differences between both curves at $\text{pK}_{a}-\text{pH}\geq-1.67$ corresponding to low pH values.

The differences to the constant pH method can be mainly attributed to the additionally occurring electrostatic screening effects in the real titration mode. In order to achieve higher association degrees $\bar{n}$ and lower pH values, a strong acid is injected into the system whose chemical species \ceH+ and \ceCl- induce a screening of electrostatic interactions between the deprotonated titrable groups [37]. The corresponding points of eigen pH values for the given choice of the reaction constant, as introduced in section 2.2, are denoted by orange stars in Fig. 2. Hence, for very low pH values and high concentrations of the excess strong acid, the charges of the titrable units are screened such that the polyelectrolyte becomes more and more ideal for decreasing pH values. Furthermore, it can be observed that for polyelectrolytes with pK${}_{a}\leq 2$, the degrees of association for $\text{pK}_{a}-\text{pH}\geq 0.25$ are significantly smaller compared with the ideal titration curve. We can attribute this finding to the explicit presence of charged species in the solution. It was experimentally found and in depth discussed in Refs. [38, 23] that the presence of ions, like in a salty solution favors a stronger dissociation of polyelectrolytes. Thus, the charged species stemming from the injected strong acid resemble a salty solution such that lower degrees of association at specific pH values in comparison to ideal titration curves can be observed.

In order to verify the presence of electrostatic screening effects, we studied a polyelectrolyte with pK${}_{\text{a}}=2$ in presence of explicit salt ions and by using the constant pH method (gray triangles in Fig. 2). More specifically, we calculated the Debye Hückel length $\lambda_{D}=\sqrt{1/(4\pi\lambda_{B}\sum_{i}c_{i}z_{i}^{2})}$ [37] where $c_{i}$ is the concentration of charged species with valency $z_{i}$ in the reaction ensemble in the real titration mode at distinct pH values (yellow squares) and added the corresponding concentrations of chemically inert salt anions and cations to the constant pH simulations. The coincidence between the curves verifies our assumption that electrostatic screening effects are mainly responsible for the differences between the constant pH and the reaction ensemble method.

Moreover, due to the increasing number of protons from the injected strong acid at pH values lower than the point of the eigen pH value, a significant decrease of the electrostatic Debye-Hückel screening length can be observed at the right side of Fig. 3. For low pH values, the results for the constant pH method indicate that the Debye-Hückel length diverges in comparison to the reaction ensemble method in the real titration mode, which can be attributed to the above discussed absence of explicit excess free protons in the constant pH method. Moreover, after a comparison between Fig. 2 and Fig. 3, one can observe that most pronounced differences between the methods are evident for $\lambda_{D}\leq 10\sigma$ which is in the order of the polyelectrolyte size. Vice versa, in order to simulate higher pH values, a strong base is added to the system whose hydroxide ions annihilate with the free protons of the weak polyelectrolyte in terms of autoprotolysis reactions (left side of Fig. 3). Thus, the Debye-Hückel lengths and the titration curves are identical in the constant pH and the reaction ensemble method for high pH values until the point of the eigen pH value is reached. The Debye-Hückel length also decreases for higher pH values, due to a significantly more pronounced dissociation of the weak acid aresulting in a high amount of free charged species. Due to these reasons, we conclude that electrostatic screening effect impose a significant influence on weak polyelectrolytes which complicates the applicability of the constant pH method at low pH values. However for weak polyelectrolytes which have a high enough $\text{pK}_{a}$ value (e.g. $\text{pK}_{a}=3$), the difference between the constant pH titration curve and the reaction ensemble titration curve in the real titration mode practically vanish (see Fig. 2).

In summary, pH dependent screening effects are not adequately reproduced by the standard constant pH method or the reaction ensemble in the sweeping mode due to the fact that both approaches consider the pH value implicitly. Nevertheless, it can be assumed that the explicit treatment of pH dependent screening effects might be relevant for proteins, based on the findings that the individual amino acids strongly differ in their deprotonation/protonation properties [2].

6 Summary and conclusion

In this article, we demonstrated the similarity between the reaction ensemble method in the sweeping mode and the constant pH method under certain conditions. Both methods can be used to study the dissociation properties of titrable groups in weak polyelectrolytes. Noteworthy, the implicit interpretation of the pH value in the constant pH method inhibits electrostatic screening effects due to the absence of all explicit free protons according to the pH value of the solution around the charged groups of the polyelectrolyte. This finding points at certain complications of the constant pH method. It is evident that these effects are mostly important at extremely high or low pH values. In fact, for moderate pH values, the constant pH and the reaction ensemble method reveal comparable results. Moreover, we proposed a new operational mode for the reaction ensemble method, which can be used to study the behavior of polyelectrolytes according to real titration procedures. Based on our findings, we conclude that the usage of the reaction ensemble method in the real titration mode is specifically preferred for polyelectrolytes with different functional groups, under confinement and for the simulation of acids and bases with moderate pKa values.

7 Acknowledgments

We thank Tobias Richter, Peter Kosovan and Kai Szuttor for fruitful discussions and helpful ideas. Financial funding by the Deutsche Forschungsgemeinschaft through the grants AR 593/7-1 and HO 1108/26-1 is gratefully acknowledged.

We dedicate this work to Wolfhard Janke on the occasion of his 60th birthday and wish him an equally scientific fruitful future.