Adsorption Isotherm and Mechanism of Ca2+ Binding to Polyelectrolyte

Polyelectrolytes, such as poly(acrylic acid) (PAA), can effectively mitigate CaCO3 scale formation. Despite their success as antiscalants, the underlying mechanism of binding of Ca2+ to polyelectrolyte chains remains unresolved. Through all-atom molecular dynamics simulations, we constructed an adsorption isotherm of Ca2+ binding to sodium polyacrylate (NaPAA) and investigated the associated binding mechanism. We find that the number of calcium ions adsorbed [Ca2+]ads to the polymer saturates at moderately high concentrations of free calcium ions [Ca2+]aq in the solution. This saturation value is intricately connected with the binding modes accessible to Ca2+ ions when they bind to the polyelectrolyte chain. We identify two dominant binding modes: the first involves binding to at most two carboxylate oxygens on a polyacrylate chain, and the second, termed the high binding mode, involves binding to four or more carboxylate oxygens. As the concentration of free calcium ions [Ca2+]aq increases from low to moderate levels, the polyelectrolyte chain undergoes a conformational transition from an extended coil to a hairpin-like structure, enhancing the accessibility to the high binding mode. At moderate concentrations of [Ca2+]aq, the high binding mode accounts for at least one-third of all binding events. The chain’s conformational change and its consequent access to the high binding mode are found to increase the overall Ca2+ ion binding capacity of the polyelectrolyte chain.


Introduction
Divalent metal ions, such as Ca 2+ , exhibit a pronounced affinity for dissolved anions like carbonate (CO 2 - 3 ) in an aqueous solution.The ion pairs, which readily associate, precipitate out of solution and form solid deposits (scale) due to their low solubility limit.2][3] These metal ions also form complexes with household products, such as detergents, and disrupt their efficacy.
5][6][7][8] Although it is not fully understood what makes an effective antiscalant polyelectrolyte, a few mechanisms for scale prevention have been proposed.The polyelectrolytes could prevent nucleation via chelating metal ions from the solution. 5,9The chelation reduces the concentration of metal ion and the likelihood of their association with the dissolved anions.Simultaneously, polyelectrolytes could also adsorb onto the surfaces of scale crystals and prevent further growth or deposition. 6,7,10,113][14][15][16][17][18][19][20][21][22][23] Cloud point measurements revealed instances of phase separation into a polymer-poor (supernatant) liquid and a polymer-rich liquid with addition of divalent ions. 24Boisvert et al. performed osmotic pressure measurements and showed that the solution behavior of polyelectrolytes in divalent salt solutions is primarily influenced by a proposed site-binding mechanism of divalent cations 25 .The site-binding mechanism facilitates bridging of non-adjacent repeat units by the divalent cations.Through a Fourier Transform Infrared (FTIR) dialysis technique, 26 Fantinel et al. identified monodentate, bidentate, and bridging modes when Ca 2+ ions bind to carboxylate groups of a polyacrylate chain.However, the relative importance of each of these binding modes-particularly that of the bridging mode-on the ability of polyelectrolyte to chelate Ca 2+ ion has not been elucidated.Furthermore, the influence of Ca 2+ ion concentration on these binding modes remains elusive.
8][29][30][31][32][33] In addition to establishing conditions for the precipitation behaviours, these theoretical models have predicted a large reduction in polymer size with addition of divalent ions, 34 beyond what is expected from electrostatic screening.The chain collapse was primarily attributed to ion bridging between non-neighboring repeat units.8][39][40][41][42][43] These investigations have each reported that the Ca 2+ ion is strongly coordinated with the polyelectrolyte chain, resulting in highly coiled conformations with a chain rigidity reminiscent of crystal-like structures.Due to the strong Ca 2+ -polyelectrolyte interactions and the number of Ca 2+ ions binding to the polyelectrolyte chain, these models hint at the overcharging of the Ca 2+ -polyelectrolyte complex.However, recent potentiometric titration assays 9 suggest that only 1 3 of the binding sites are occupied before a fully charged polyelectrolyte chain reaches its saturated value of Ca 2+ -binding capacity.
The saturation value in the adsorption isotherm, which describes the maximum amount of Ca 2+ that can be bound to a polyelectrolyte chain, reflects the chelating capacity of the polyelectrolyte chain.However, the molecular principles that govern the corresponding adsorption isotherm have not been addressed.Specifically, the interplay between the site-binding nature of Ca 2+ -polyelectrolyte interactions and the conformational transitions of the polyelectrolyte, aimed at enhancing both the chelating capacity and the solubility of the Ca 2+ -polyelectrolyte complex, has not been explored.
In this study, we address the mechanism of Ca 2+ adsorption onto a polyacrylate chain.We construct an adsorption isotherm to describe the binding behaviors of Ca 2+ ions to polyacrylate.We then determine the different binding modes accessible for Ca 2+ binding to the chain and quantify their impact on the adsorption isotherm.In our follow-up manuscript, 44 we address questions related to Ca 2+ ion-mediated association between polyelectrolyte chains.The rest of the manuscript is organized as follows.
In Section 2, we present a Hamiltonian Replica Exchange Molecular Dynamics (HREMD) protocol to selectively bias Ca 2+ -polyelectrolyte interactions and efficiently sample the configurational space.We then introduce a free energy perturbation approach coupled with molecular dynamics to compute the adsorption isotherm describing Ca 2+ binding to the polyelectrolyte chain.We discuss the results of these calculations in Section 3 and present the conclusions in Section 4.

Models and Methods
We investigated the mechanism of calcium ion binding to polyacrylate chain using all-atom molecular dynamics (MD) simulations.Our simulation system consists of a single polyacrylic acid (PAA) chain with 32 repeat units, solvated in a cubic water box with an edge length of 12 nm.All repeat units on the polymer are charged, consistent with the solution conditions for antiscalant activity(i.e., solution pH ∼ 10).Sodium ions were added for electroneutrality.The average end-to-end distance of such a polymer was ∼ 5 nm.We chose our box dimensions so that the polymer does not interact with its periodic image.During our modeling of calcium ions, we observed that the "full" electrostatic charge force field parameters overestimated calcium ion binding to carboxylate groups on the polyacrylate chain, resulting in a charge inversion of the polymer chain inconsistent with experimental reports (see Supporting Information).Duboué-Dijon et al., who investigated calcium ion binding to insulin, reported that such an overestimation was due to inadequate treatment of electronic dielectric screening when using full charges on ions with non-polarizable forcefields. 46The correct energetics of ion-pair formation could be captured by molecular models with explicit polarization. 47However, a general-purpose force field with explicit polarization, tested to reproduce the properties of polyacrylates, was not readily available, and polarizeable force fields introduce large computational expense.
Alternatively, in our models, we included electronic dielectric screening by uniformly scaling the charges of all solute atoms.This approach, known as the electronic continuum charge correction (ECC), is a mean-field method that attempts to mimic charge-carrying species within an electronic dielectric continuum. 48,49While the ECC scheme is a physically meaningful concept, it is primarily used as an ad hoc solution to incorporate electronic polarization into an otherwise non-polarizable force field.ECC schemes tend to fail in systems with a discontinuity in the high-frequency dielectric constant. 50Moreover, force field parameterizations implicitly account for electronic polarization effects to some extent.Introducing additional ECC correction may overly compensate for electronic polarization effects, resulting in a significant underestimation of cohesive energy density and leading to unphysical solution behavior. 51evertheless, the ECC scheme yields meaningful observations in common electrolyte systems where the high-frequency dielectric constant remains uniform throughout.
Within the scope of our study, the interactions of interest involve calcium ions and carboxylate groups on the polyacrylate chain.3][54][55] We particularly employ the ECC scheme reported by Jungwirth and coworkers, who provided parameters for modeling calcium and other ions in an aqueous solution. 53In their work, the authors uniformly scaled the charges on ions by a factor of 0.75.Lennard-Jones parameters describing dispersion interactions of corresponding ions were optimized to reproduce ab initio molecular dynamics results for ion-pairing and neutron scattering experiments.Such a parameter set accurately described the properties of calcium ions in an aqueous solution and their association with carboxylate groups on amino acids.We used the same parameter set to model calcium and other ions in our simulation system.
When modeling the electronic polarization effects of the polyacrylate chain in an aqueous solution, we scaled the partial charges of all the atoms on the polyacrylate chain by 0.75 and used Lennard-Jones parameters reported by Mintis et al. 45 to describe dispersion interactions.Although our approach is a commonly accepted practice, 55 we emphasize that it is not rigorous.However, different properties of polymers strongly depend on their chain length, and finding the right strategy to optimize their dispersion interaction parameters is not obvious.
Since the "full-charge" parameters by Mintis et al. 45 reproduce the structural properties of the polyacrylate chain in a salt-free solution, we validated the predictions of the "scaled-charge" model against the former.Even though we did not re-optimize the Lennard-Jones interaction parameters of the polyelectrolyte chain, the conformational flexibility of a polyelectrolyte chain in a salt-free solution did not change.Additionally, we observed identical distributions for structural properties when compared to the full-charge" model (see Supporting Information).
Simulation methodology: We used GRO-MACS 2022.3 [56][57][58] patched with PLUMED 2.8.1 [59][60][61] and employed the following protocol to conduct our molecular dynamics simulations of the systems reported in Table 1.First, we minimized the energy of the system using the steepest descent algorithm until the maximum force on any atom in the system was smaller than 100 kJ mol −1 nm −1 .Next, we equilibrated the system at constant NPT conditions with pressure and temperature set to 1 atm and 300 K, respectively.While fluctuations in the energy and box size minimized within a few nanoseconds of the equilibration run, it took tens of nanoseconds to equilibrate the chain structure.We used the Berendsen barostat 62 with the velocity-rescaling stochastic thermostat during the first 10 ns of the equilibration.For the remainder of the equilibration run, we switched to the Parrinello-Rahman barostat 63 with the Nosé-Hoover chain thermostat for the accurate reproduction of the thermodynamic properties of the system. 64,65Following the equilibration run, we conducted production MD simulations of these systems in the NVT ensemble using the Nosé-Hoover chain thermostat.
We employed the leap-frog time integration algorithm with a finite time step of 2 fs to integrate the equations of motion.Additionally, we utilized the LINCS constraint algorithm to convert all bonds with hydrogen atoms into constraints. 66We applied periodic boundary conditions along all three spatial axes and used the particle mesh Ewald (PME) method with a minimum Fourier spacing of 0.12 nm to calculate the long-ranged electrostatic interactions. 67,68We applied a cutoff distance of 1.2 nm for computing van der Waals interactions.We used the same distance as the real-space cut-off value while computing the PME electrostatics.
Need for enhanced sampling molecular simulations: Although the ECC scheme with the non-polarizable force field greatly reduced the PAA-Ca 2+ binding/unbinding relaxation times, regular MD simulations were unable to efficiently sample the polymer conformational space in an aqueous CaCl 2 solution.Even a microsecond long trajectory was not sufficient to sample polymer conformations in any of the systems listed in Table 1 with calcium numbers higher than 4 Ca ions.We direct readers to the Supporting Information document for the relevant data and discussion.
To address the challenges associated with polymer conformational sampling in an aqueous CaCl 2 solution, we employed Hamiltonian Replica Exchange Molecular Dynamics (HREMD). 69,70Our HREMD framework is based on the flexible implementation of the REST2 variant, 69 as previously reported by Bussi. 70We introduced a parameter  to selectively bias the interactions between the polymer and ions, as well as the dihedral potential components of the Hamiltonian.The charges of the ions and the polymer were scaled by a factor of √ , while their Lennard-Jones interaction parameter () was scaled by .Similarly, the polymer dihedral potential was also scaled by .With this scheme for -parameterized Hamiltonian,  = 1 corresponds to the system of interest with full-scale interactions.We determined that the parameterized Hamiltonian with  = 0.67 rapidly sampled the polymer conformational space.Coordinate exchange between neighboring replicas were attempted every 500 steps.We utilized 16 replicas, with  values ranging from 1 to 0.67 (geometrically spaced), to simulate polymer conformations in an aqueous solution containing 8 CaCl 2 or 16 CaCl 2 .For higher Ca 2+ numbers, we increased the number of replicas to 24.This combination of number of replicas and the -range yielded acceptable exchange probabilities (∼ 0.3) between neighboring replicas.All the relevant results reported in the subsequent sections were obtained by averaging over a 250 ns production HREMD run, conducted at constant volume and a temperature of 300 K.
Computing ion adsorption isotherm from molecular dynamics simulations: The primary objective of the current work is to investigate Ca 2+ chelation onto a model polyacrylate chain.In an aqueous CaCl 2 solution containing a polyacrylate chain, a dynamic equilibrium exists between calcium ions that are freely dispersed in the solution (Ca 2+ aq.free ) and the calcium ions adsorbed per monomer of the polyacrylate chain (AA − Ca 2+ ).
Here, AA represents the concentration of re-peat units on the polymer chain that are not bound to any calcium ions.An adsorption isotherm, quantifying Equation (1), describes calcium ion chelating ability of a model polyacrylate chain.We constructed the isotherm by plotting the number of calcium ions adsorbed per monomer (AA − Ca 2+ ) against the concentration of calcium ions that are freely dispersed in the solution (Ca 2+ aq.free ).Computing AA − Ca 2+ from the simulation trajectory is straightforward.We calculated the number of calcium ions that are within 0.7 nm (see SI) of the polymer atoms at each frame.This quantity was then divided by the number of repeat units per chain (i.e.32 in this case) and ensemble averaged over the trajectory.
We .][73][74] In brief,  CaCl 2 represents the change in free energy when adding (or removing) a unit of CaCl 2 , This is expressed in Equation (2) as the sum of the corresponding ideal gas component ( id CaCl 2 ) and the residual component ( R CaCl 2 ).
Here,  0 Cl − and  0 Ca 2+ are the respective chemical potentials of Cl -and Ca 2+ at the reference pressure of 1 bar.We used the tabulated value of  0 Cl − from the NIST-JANAF thermochemical tables 75 .Moučka et al. estimated the value for  0 Ca 2+ 76 , and we employed the same value in our calculations.N CaCl 2 represents the num-ber of CaCl 2 units in the simulation box,   denotes the Boltzmann constant,  0 = 1, bar is the reference pressure, and ⟨⟩ stands for the average volume obtained from NPT simulations of a system with a specific composition.
To calculate  R CaCl 2 , we gradually decoupled a Ca 2+ and two Cl -from the system.The initial configuration for this study came from the most probable polymer conformation identified from the previously described enhanced sampling simulations.We then tagged a Ca 2+ ion that was adsorbed onto the polymer backbone.In a solution without the polyacrylate chain, a randomly selected Ca 2+ ion served as the tagged ion, and a comparison point between the two systems.Regardless of the system, two randomly chosen Cl -ions were tagged.
The electrostatic interactions of the tagged Ca 2+ ion with other particles in the system were gradually turned off in 11 stages ( = 1, 0.9, . . ., 0).Here,  = 1 represents the system with fully activated electrostatic interactions of the tagged calcium ion, while  = 0 corresponds to complete deactivation.Each stage consisted of a series of molecular dynamics simulation steps, beginning with energy minimization, followed by 5 ns of equilibration and a 10 ns production simulation run at a pressure of 1 bar and a temperature of 300 K.The output from the production step of the th stage (  ) served as the initial configuration for the MD simulation steps of the ( +1)th stage ( +1 ).The same methodology was then applied to deactivate the electrostatic interactions of each tagged chloride ion and also the van der Waals interactions of the tagged calcium ion and the two tagged chloride ions.We then used the Bennett Acceptance Ratio (BAR) 77,78 approach, natively implemented in the GROMACS MD package 56 , to estimate the residual chemical potential ( R CaCl 2 ) from these different stages of MD simulations.Notably, given the positiondependent nature of the unbinding free energy of a Ca 2+ ion from the polymer backbone, we calculated the unbinding free energy for each adsorbed Ca 2+ ion individually and used their average to estimate  R CaCl 2 .

Results and Discussion
The ability of a polyacrylate chain to sequester Ca 2+ ions is intricately tied to the bulk solution concentration of Ca 2+ , denoted as Ca 2+ aq.free .In Figure 1  From the plateau in the adsorption isotherm, we note that the maximum binding capacity of a model polyacrylate chain is 0.40 ± 0.05 Ca 2+ ions per monomer, which aligns with recent potentiometric titration experiments 9 .Although the adsorption isotherm bears some resemblance to a Langmuir model, we observe that the Ca 2+ binding does not obey the model assumptions.Notably, the conformational flexibility of the polyelectrolyte chain results in distinct chemical environments around each binding site, rendering the binding sites nonequivalent.
We demonstrate this phenomenon in a polyelectrolyte solution corresponding to Ca 2+ aq.free = 0.026mol/kg.First, we identify the system configuration that corresponds to the polymer chain with the most probable radius of gyration (Figure 2a).Then, we independently unbind each of the bound Ca 2+ ions.We employ the free energy perturbation approach described in Section 2 to compute the unbinding free energy, and report the corresponding values in Figure 2b.
The unbinding free-energies of the 11 Ca 2+ ions broadly fall into two classes.We categorize these binding sites into a "low" binding mode (≤ 818 kJ/mol) and "high" binding mode (> 818 kJ/mol).From the binding sites depicted in Figure 2a along with the accompanying free energies in Figure 2b, we identify that the high binding mode is approximately 5-10 kJ/mol energetically more favorable than the low binding mode and facilitates a Ca 2+ ion-bridge between non-neighboring carboxylate groups on the polyacrylate chain.
The various binding modes and their coupling to polyelectrolyte conformation could impact the number of Ca 2+ ions sequestered.We determine a binding mode by tracking the number of carboxylate oxygen atoms that are around a Ca 2+ ion.From the radial distribution of carboxylate oxygen atoms around a Ca 2+ ion, we identify that their most probable separation is about 0.35 nm (see Supporting Information).We compute the probability of finding a varying number of carboxylate oxygen atoms within 0.35 nm from a Ca 2+ ion and report this as a function of Ca 2+ aq.free in Figure 3.We observe that the low binding mode, in which two carboxylate oxygen atoms bind to a given Ca 2+ ion, remains dominant at all concentrations of Ca 2+ aq.free studied.However, with an increase in Ca 2+ aq.free , the high binding mode corresponding to ion-bridging becomes more favorable.Interestingly, when Ca 2+ aq.free exceeds 2.18 × 10 −2 mol/kg, this trend reverses: further increases in Ca 2+ aq.free promote the low binding mode once more.We hypothesize that this nonmonotonic trend in the population of different binding modes arises from the enhanced electrostatic screening at higher Ca 2+ concentrations.
Nevertheless, at moderate and high concentrations of [Ca 2+ aq.free , nearly 1 3 of the binding events are due to the high binding mode.This high coordination binding environment is seen to have a large impact on the polymer size and conformation.We investigate the coupling between the binding modes and the polyelectrolyte chain conformation by tracking the chain radius of gyration (R g ) as a function of Ca 2+ aq.free .We report these results in Figure 4. We note from Figure 4(a) that at low concentrations of Ca 2+ aq.free , where the population of the high binding mode is insignificant, the polymer chain adopts an extended conformation with an R g of approximately 1.65-1.9nm (Figure 4(b)).At these concentrations, Ca 2+ ions bind to at most one carboxylate group on the polyacrylate chain.As Ca 2+ aq.free increases and the population of high coordination binding sites subsequently increases, the polymer conformation transitions to a collapsed state (Figure 4(c)).Here, the high coordination binding sites, which bridge two strands of the polyelectrolyte chain, are nearly as prominent as the low coordination binding sites located on the solvent-exposed side of each strand.Intriguingly, in solutions with high concentrations of Ca 2+ ions, although some of the bridging events are disrupted, the conformation of a polyelectrolyte chain still resembles that of a collapsed state (Figure 4(d)).The disruption of the bridging events due to enhanced screening from the other ions in the system increases the conformational flexibility and hence the average chain size.
Since the high coordination binding sites are energetically more favorable, we anticipate that access to a larger population of these binding sites would enhance the ability of a polyelectrolyte chain to sequester Ca 2+ ions.To explore this, we investigate the inverse problem; we limit the polyelectrolyte chain to only binding sites with low coordination environments and  construct the adsorption isotherm.We achieve this by first identifying the polymer conformation that corresponds to the most probable R g in a system with no added Ca 2+ ions.Utilizing this chain conformation, with harmonic restraints imposed, we prepare polyelectrolyte solutions with added Ca 2+ ions.We compute the Ca 2+ adsorption isotherm from these systems and compare it to that obtained from unrestrained simulations reported earlier.We show these results in Figure 5.
We find that at low concentrations, the number of Ca 2+ ions adsorbed on the polyacrylate chain in both the restrained and unrestrained simulations is indistinguishable.This is because, at lower concentrations of Ca 2+ ions in the solution, unrestrained polyelectrolyte chains can only access the low coordination binding sites.However, we see a noticeable difference at moderate and high concentrations, where the unrestrained polyacrylate chain favors the high coordination binding sites.The restrained simulations, which do not have access to high binding modes, consistently adsorb fewer Ca 2+ ions per chain when compared to the unrestrained simulations.
These observations indicate that a polyelec-   trolyte chain's ability to efficiently chelate Ca 2+ ions is impacted by its access to a large number of high coordination binding sites.These sites are energetically more favorable for ion binding, and as such, they provide a more stable environment for the ions, leading to a more effective sequestration.This insight could prove useful in applications where selective and efficient ion capture is paramount, such as in water treatment processes or biomedical applications.

Summary
Using all-atom molecular simulations coupled with a free energy perturbation approach, we constructed an adsorption isotherm to describe the binding of Ca 2+ ions to a model polyacrylate chain.Analysis of the adsorption isotherm revealed that the per-monomer Ca 2+ ion binding capacity of a fully charged polyelectrolyte chain saturates at a value of 0.40 ± 0.05.This saturation value correlates with the binding modes accessible to Ca 2+ ions as they bind to the polyelectrolyte chain.Two predominant binding modes were identified: one mode involves Ca 2+ ions binding to at most two carboxylate oxygen atoms on a polyacrylate chain, and the other involves Ca 2+ ions binding to four or more carboxylate oxygen atoms.
The population of low binding mode sites remains high across all concentrations of Ca 2+ in the solution.Nevertheless, at least one-third of the binding events at moderate and high concentrations of Ca 2+ ions in the solution are defined by the high binding mode events.These binding events, while responsible for enhancing the Ca 2+ ion binding capacity of a polyelectrolyte chain, also lead to the collapse of the polyelectrolyte chain's conformation.The solution concentration of Ca 2+ ions corresponding to the conformational transition falls within the range close to the supernatant (polymer-poor) side of the phase behavior reported by Sabbagh et al. 24 The adsorption isotherm constructed in this study suggests that, at these concentrations, the collapsed polyelectrolyte chain had attained the saturation value of the Ca 2+ binding capacity.
The collapsed polyelectrolyte chain, with all available binding sites saturated, may indicate the onset of putative phase separation.This observation carries important implications for the chelating capacity of a polyelectrolyte towards Ca 2+ ions, and consequently, for its antiscalant activity.Yet, polyelectrolyte-divalent ion complex falling out of equilibrium is inherently a multi-chain problem.In our follow up work, 44 we tackle this problem and establish the conditions under which a polyelectrolyte in divalent salt solution falls out of equilibrium.

Figure 1 .
Figure 1.Isotherm describing Ca 2+ ion binding to a polyacrylate chain with 32 repeat units.Error bars represent one standard deviation around the sample mean.Note: The apparent positive slope of the last two data points is a consequence of large uncertainties in determining corresponding solution concentration of free Ca 2+ in the system.The differences in their vertical axis values are very minimal.

Figure 4 .
Figure 4. a) Radius of gyration of 32-mer polyacrylate chain at different concentrations of Ca 2+ aq.free .Dominant binding modes at b) low, c) moderate, and d) high Ca 2+aq.free concentrations.The conformation in (c) corresponds to the R g minimum in (a)

Figure 5 .
Figure 5. Ca 2+ ion binding isotherm to a polyacrylate chain with 32 repeat units restrained to an extended coil conformation and unrestricted.

Table 1 . Composition of the cubic simulation box with an edge length of 12 nm, containing a single sodium polyacrylate chain with 32 repeat units, and at various
numbers of CaCl 2 in water.
determined Ca 2+ aq.free by equating the chemical potential of the Ca 2+ ions in the system ( System , we present the adsorption isotherm and describe the extent of Ca 2+ binding onto a 32-mer polyacrylate chain as a function ofCa 2+aq.free .At low Ca 2+ concentrations, most of the binding sites on the polyacrylate chain are available, and an increase in Ca 2+ aq.free results in a rapid increase in the number of Ca 2+ ions sequestered by the polyacrylate chain.However, at moderate concentrations, the accessible binding sites become occupied, and the polyacrylate chain saturates with Ca 2+ .